Example using nonlinear equation of state

[glwagner]

We should add an example that uses TEOS10EquationOfState from SeawaterPolynomials.

cc @BrodiePearson

[BrodiePearson]

Sure, although I ran into an issue after our discussion on this the other day. EDIT: I just realized that I had an erroneous factor of density in there - so the below error may just be because I was imposing a very high momentum flux, but I'm not sure why that manifests as a negative salinity.

When I increase the momentum flux with the non-linear EOS switched on, the model crashes with a DomainError as the term in the following square root is negative (presumably the absolute salinity is the problem). Reducing the momentum flux from 1e-3 to 1e-4 allows the model to run.

https://github.com/CliMA/SeawaterPolynomials.jl/blob/9e3342ec405bff01dee23e6309dd9de6fe342aa8/src/TEOS10.jl#L73

Here is a script that reproduces the error: https://github.com/BrodiePearson/Oceananigans.jl/blob/Nonlinear_EOS/examples/nonlinear_EOS_crash.jl

[glwagner]

I guess this is not the issue (since Qᵇ = 0), but just be aware that this line:

https://github.com/BrodiePearson/Oceananigans.jl/blob/838cc38552ba8c4acff583aef840157cf7f145da/examples/nonlinear_EOS_crash.jl#L104

is not correct if Qᵇ is buoyancy flux (as the numbers of the commented out line suggest it is)--- it is only correct is Qᵇ = 0 is heat flux. Which may be what you want anyways for this problem...

Note also that TEOS10 requires a reference density. Thus this line:

https://github.com/BrodiePearson/Oceananigans.jl/blob/838cc38552ba8c4acff583aef840157cf7f145da/examples/nonlinear_EOS_crash.jl#L93

introduces a second reference density that isn't equivalent to the one used by TEOS10 (the default for TEOS10 in SeawaterPolynomials.jl is 1020 kg / m^3). The two can be synced by instantiating the TEOS10 equation of state object first:

eos = TEOS10EquationOfState(FT)

and then using eos.refererence_density. You can also adjust the reference density for the equation of state by writing

eos = TEOS10EquationOfState(FT, reference_density=1027)

As for negative salinity... the only thing I can think of is that its a numerical error due to a violation of a CFL condition (or something similar). Does reducing the CFL used in adaptive time stepping to 0.1 or 0.05 help? You might also reduce the initial time-step. With a strong momentum flux, you'll get strong surface velocities initially, and it may be best to start time-stepping with a very small time-step.

Yo could try to replace this line:

https://github.com/BrodiePearson/Oceananigans.jl/blob/838cc38552ba8c4acff583aef840157cf7f145da/examples/nonlinear_EOS_crash.jl#L223

with

wizard = TimeStepWizard(cfl=0.1, Δt=0.1, max_change=1.1, max_Δt=10.0)

[glwagner]

Side note: hopefully we will have RK3 time stepping at some point, and then we will be able to dramatically improve our adaptive time-stepping as there will be no penalty in updating the step size every time step.

[ali-ramadhan]

Not really a full example but the model setup docs now show how to use the equations of state from SeawaterPolynomials.jl: https://clima.github.io/OceananigansDocumentation/stable/model_setup/buoyancy_and_equation_of_state/#Idealized-nonlinear-equations-of-state

[glwagner]

Ah nice, that's great.

I do think a full example would be nice since that's the entry point for a lot of people. @BrodiePearson might have a good script to start an example from that showcases a few features at the same time, including spatially-dependent boundary conditions and a nonlinear equation of state... ?