What's the best way to avoid CFL-related instabilities in the early simulation steps of LES

[tomchor]

I'm running a LES of the upper ocean with a Smagorinsky closure. In the beginning of the simulations the Courant number varies a lot (probably because of the small-scale noise) so I have to set the CFL to something really small to get passed that. (Something like 0.02!) The problem is that eventually the Courant number stabilizes and I get stuck with a simulation that takes much smaller steps than it could, and therefore doesn't use computational resources efficiently.

In other codes that I know this is generally counter-acted by starting the first couple hundred steps of an LES with a constant eddy diffusivity, and only then turning on the LES closure. Is there an easy way to do that with Oceananigans? Sorry if I missed something, but I searched the docs and couldn't find anything.

Ideally I'd wanna avoid starting a very brief simulation with a constant viscosity and then pickup that simulation with another script to run it "for real". Alternatively having a time-dependent CFL limit might work, but that seems unnecessarily convoluted...

Thoughts?

[ali-ramadhan]

I've used a time-dependent CFL in the past with some good results. It's easy to implement in the simulation progress function but a bit awkward. This is how I did it for the stratified Couette flow validation experiment (the script needs some real TLC):

Oceananigans.jl/validation/stratified_couette_flow/stratified_couette_flow.jl

Lines 237 to 244 in cadfdae

	# We will ramp up the CFL used by the adaptive time step wizard during spin up.
	cfl(t) = min(0.01t, 0.1)
	function print_progress(simulation)
	model = simulation.model
	clock = model.clock
	wizard.cfl = cfl(model.clock.time)

Switching model properties is not allowed since that would modify the model's Julia type. You also might not want to create two models as that would allocate twice as much memory.

Might not be too hard to write a function like

new_model = switch_out_model(old_model, closure=new_closure, ...)

that constructs and returns a new model using all the arrays from the old model to avoid allocating any new memory.

It would perform a somewhat similar job to the checkpointer set! function: https://github.com/CliMA/Oceananigans.jl/blob/master/src/OutputWriters/checkpointer.jl#L236

[glwagner]

Order shouldn't matter if we are using kwargs though.

@ali-ramadhan fair point that similar (the function I was thinking of) returns an uninitialized object. Perhaps copy_model? swap_model_properties could work too.

@tomchor I think the problem is that the constructor assumes background_fields kwarg / user input is "flat", but then gathers the background fields into velocities and tracers categories. (This is a good lesson that its best for user input to correspond as closely as possible to the object that's created.) So you'd need to write:

new_model = copy_model(old_model; closure=SmagorinskyLilly(),
                       background_fields=merge(old_model.background_fields.velocities, old_model.background_fields.tracers))

Smoothing out the user experience could be part of the PR for this feature. There will still be some bumps, especially due to the fact that model.diffusivities depends on model.closure. This would be fixed by storing the eddy diffusivities in model.closure rather than model.diffusivities (which we know how to do now since we now know how to use adapt_structure properly), which hopefully we'll have before too long.

[tomchor]

@glwagner Just FYI, that didn't work. Here's an example with a more simple code:

julia> new_model  = swap_model(model; closure=closure_LES,
                          background_fields=(u=old_model.background_fields.velocities.u, b=old_model.background_fields.tracers.b))
ERROR: MethodError: no method matching regularize_background_field(::Type{Face}, ::Type{Center}, ::Type{Center}, ::Oceananigans.Fields.FunctionField{Face,Center,Center,Clock{Float64},Nothing,typeof(u_g),RegularCartesianGrid{Float64,Periodic,Bounded,Bounded,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}}, ::RegularCartesianGrid{Float64,Periodic,Bounded,Bounded,OffsetArrays.OffsetArray{Float64,1,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}, ::Clock{Float64})
Closest candidates are:
  regularize_background_field(::Any, ::Any, ::Any, ::BackgroundField{var"#s64",P} where P where var"#s64"<:Function, ::Any, ::Any) at /glade/u/home/tomasc/.julia/packages/Oceananigans/iSJKi/src/Fields/background_fields.jl:59
  regularize_background_field(::Any, ::Any, ::Any, ::Function, ::Any, ::Any) at /glade/u/home/tomasc/.julia/packages/Oceananigans/iSJKi/src/Fields/background_fields.jl:62
...

So we apparently need to pass the functions:

julia> new_model  = swap_model(model; closure=closure_LES,
                          background_fields=(u=old_model.background_fields.velocities.u.func, b=old_model.background_fields.tracers.b.func))

Although this little quirk can be dealt with inside the function I guess, so that the user doesn't need to worry about it.

This would be fixed by storing the eddy diffusivities in model.closure rather than model.diffusivities (which we know how to do now since we now know how to use adapt_structure properly), which hopefully we'll have before too long.

I upvote that!

[glwagner]

We can hack around it, but I think a more robust and long-term improvement is to change the way that background_fields are parsed. Fixing that up might have additional benefits down the line.

[tomchor]

That sounds okay with me. Just for the record, this worked, and now my simulation runs much smoother. Should I do something with the code for now? I feel that until we make it work seamlessly (meaning that no special treatment of the background fields is required of the user) I shouldn't add it to Oceananigans. But I don't know if everyone agrees.

[glwagner]

I think adding

BackgroundFields(background_fields::NamedTuple{(:velocities, :tracers)}, args...) = background_fields

to

https://github.com/CliMA/Oceananigans.jl/blob/df9b830e574261efb1b7b0ee43b783be2381ceaf/src/Fields/background_fields.jl

will make copy_model work without needing special translation for background_fields. This makes is so that this line:

Oceananigans.jl/src/Models/IncompressibleModels/incompressible_model.jl

Line 140 in df9b830

background_fields = BackgroundFields(background_fields, tracernames(tracers), grid, clock)

in the constructor for IncompressibleModel returns the right thing when background_fields is already "regularized" and organized into a named tuple with fields velocities and tracers.

There are other bumps that might arise since we not covering all cases in this discussion.

More broadly this suggests a flaw in the design of BackgroundFields, but fixing that deeper issue is for another time / PR.