-
Notifications
You must be signed in to change notification settings - Fork 195
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Include background tracer and velocities in closure tendency terms #3646
base: main
Are you sure you want to change the base?
Conversation
Probably before merging this we also want to use the total velocities in the stress computation |
But I wonder also if we should add a feature to |
…liMA/Oceananigans.jl into glw/background-flux-divergence
Ok @hdrake @liuchihl I have added the option to include or not the background field when computing closure fluxes. They are noit included by default. If you want to include them you need to build the background_fields = BackgroundFields(; background_closure_fluxes=true, b=B) where Let me know if this seems like a good interface and also if it works. |
Looks like a good interface to me. But is it on purpose that there is only support for background fields in the @liuchihl will test it in our configurations. |
Well yes, it's substantial effort to support background fields. So we implemented it in the nonhydrostatic model first. Nobody has requested having background fields for the hydrostatic model. It's not impossible but might require some thinking if it's going to work with the more complicated turbulence closures (like CATKE or k-epsilon) that sometimes get used for hydrostatic applications. Since the nonhydrostatic model is fast (at least on one GPU) the hydrostatic model is mostly important for simulations on the sphere (although this statement needs to be evaluated more carefully for complex domains when we have a proper nonhydrostatic solver). |
@glwagner, can you help us understand how the diffusive flux calculations here interact with gradient boundary conditions? We've run some simulations and seem to be getting the expected behavior in the interior, where buoyancy tendencies are due to the convergence of both the perturbation and background diffusive fluxes. However, we are not seeing the expected behavior from the flux convergence due to boundary conditions (on the domain, not even considering immersed boundary conditions yet). My expectation was that the
which in our case should be the sum of the perturbation and background tracer fields. Instead, it seems that our solutions are behaving as though the |
Oceananigans.jl/src/ImmersedBoundaries/immersed_boundary_condition.jl Lines 100 to 101 in 738d172
Contrast this with the routine for Oceananigans.jl/src/ImmersedBoundaries/immersed_boundary_condition.jl Lines 103 to 107 in 738d172
These gradients are then used to compute the flux, for example Oceananigans.jl/src/ImmersedBoundaries/immersed_boundary_condition.jl Lines 117 to 122 in 738d172
I don't know exactly what it means for a gradient to be applied to the field. Can you please clarify? |
Additionally can you help me understand why you are using |
@glwagner Thanks for implementing the total tracer diffusive flux at a high level. After running several tests, I found it to work exceptionally well! I conducted a series of tests: 1) comparing 1D vs 3D, 2) with and without the Coriolis force, and 3) with and without the immersed boundary. Everything looks great! Here are some simple examples on a rotated coordinate:
nonconstantdiffusivity250days-theta.0.002_Nx4_Ny4_smallf_zlargerf.mp4
nonconstantdiffusivity8days-theta.0.2_Nx4_Ny4_immersed_3Dfields_withcrossflux.mp4The only caveat mentioned by @hdrake is that |
Flux boundary conditions are both conceptually and practically easier to implement when the boundary fluxes are zero or constant. They can be trickier when they depend on interior flow variables. In our case, for example, the boundary condition on the perturbation variable Imposing a flux boundary condition requires knowing the diffusivity |
I just meant where in the code the gradient boundary conditions get imposed, which you've shown us is in the calculation of the gradients that feed into the downgradient diffusive fluxes that are used in the diffusive flux divergence contribution to the tracer tendencies. Thanks! |
Okay great. I would only add, I think it's clearer to think of the gradient as being used to diagnose the cross-boundary flux (rather than imposed). I guess the point here is that there is actually an apparent flux of tracer into the perturbation field because of the presence of the background. So we are imagining that the background is being maintained by some large scale circulation which is ultimately the source of this apparent flux. |
Thanks @liuchihl ! |
Hi @hdrake @liuchihl please review this PR and let me know if things look good. I've tried to consistently add the background fields to the fields that go into computing both tracer fluxes and subgrid stresses. It makes things a little more complicated, but the code is more general now. I think with this change, things like a background shear will be used in computing a subgrid stress divergence. But the shear is still not used when computing a nonlinear diffusivity (this could be added though...) |
Thanks @glwagner! We haven't been using background velocities in our setups, so won't have as much to say on that yet. But we'll continue testing the background tracer fluxes. @samlewin, are you using background shear in addition to background tracer fields in your configurations, or is the shear just in your initial conditions? This PR might be relevant. |
Should we add the test @hdrake suggests in this PR or leave it for later? Happy either way... I will elaborate a bit my train of thoughts and where my question came before. I opened the PR in an attempt to review. I looked at the file changes, I like the new feature... if I had seen that some test was added in the CI to test the new feature and the test suite was passing I'd feel very comfortable with it and would be more keen to approve. Now that I don't see that I need to look into the details of what happened here before I feel comfortable to approve (which I'm happy to do, but perhaps on Fri or next week) |
+1 for CI tests in this PR! Also, can someone please update the first comment with a brief example/description of the usage? I'm having a bit of trouble visualizing it just from the source code. |
agreed for the test. The rest looks good |
@glwagner
|
What code are you running? |
@liuchihl Have you tried running on the CPU to see if you get a more useful error message? I think
could just be the result of a typo in an expression like |
I was testing my internal tide example, where the code is written as a function, and the run script is used to adjust the parameters. |
Thanks, I'll give that a try. I could be wrong, but my guess is that the issue might be related to the use of GPU. |
@ali-ramadhan background_fields = Oceananigans.BackgroundFields(; background_closure_fluxes=true, b=B̄_field) But it runs fine like the usual way background_fields = (; b=B̄_field)
|
@liuchihl I think the best way forward is to write a simple test that illustrates the error. Then I can help fix the error to make the test pass. Once that is done, we may be ready to merge this PR. What do you think? PS it is always best to work with minimal examples, and to paste code directly into a discussion stream (rather than providing links). This will help us keep up an efficient workflow. |
Sure, I agree with that! I will work on that and let you know how it goes.
For sure, sorry about that, I haven't been able to create an MWE for this specific issue because I don't understand the problem yet. |
This does it: using Oceananigans
grid = RectilinearGrid(size=1, x=(0, 1), topology=(Periodic, Flat, Flat))
B(args...) = 0
background_fields = Oceananigans.BackgroundFields(; background_closure_fluxes=true, b=B)
model = NonhydrostaticModel(; grid, background_fields) You can use this in a test, something like B(args...) = 0
function time_step_background_fields_with_closure_fluxes(arch)
grid = RectilinearGrid(arch, size=1, x=(0, 1), topology=(Periodic, Flat, Flat))
background_fields = Oceananigans.BackgroundFields(; background_closure_fluxes=true, b=B)
model = NonhydrostaticModel(; grid, background_fields)
time_step!(model, 1)
return true
end
@test time_step_background_fields_with_closure_fluxes(arch) If it throws an error, the test fails. This is a basic unit test which is always a good idea to write first. The next step is to write a test that confirms the functionality works correctly. This is an example of a hierarchy of tests proceeding from simple to complex. |
Since the code is just a few lines, it's easy to figure out what's going on in the REPL. The top of the error message says:
So let's go look at this line. It is: Oceananigans.jl/src/Models/NonhydrostaticModels/nonhydrostatic_tendency_kernel_functions.jl Line 35 in f89445d
Ok. Does julia> background_fields.u
ERROR: type BackgroundFields has no field u
Stacktrace:
[1] getproperty(x::Oceananigans.Models.NonhydrostaticModels.BackgroundFields{true, @NamedTuple{…}, @NamedTuple{…}}, f::Symbol)
@ Base ./Base.jl:37
[2] top-level scope
@ REPL[7]:1
Some type information was truncated. Use `show(err)` to see complete types. Ok, we found the error. But where is julia> background_fields.
tracers
velocities Ok, so maybe we want Oceananigans.jl/src/Models/NonhydrostaticModels/background_fields.jl Lines 32 to 38 in f89445d
confirms what we found, that julia> background_fields.velocities.u
ZeroField{Int64} so yes, using This is an easy fix. Just to go through the motions and get a hang for how this works, let's add the test first, and then make the fix. |
I fixed the conflicts so this PR should be ready for work |
Thanks for explaining these in detail! |
CliMA#3646 (comment) @glwagner implementing the fix suggested by Greg
#3646 (comment) @glwagner implementing the fix suggested by Greg Co-authored-by: Chih-Lun Liu <[email protected]>
To summarize where we're at: we need at least two unit tests:
B(args...) = 0
function time_step_background_fields_with_closure_fluxes(arch)
grid = RectilinearGrid(arch, size=1, x=(0, 1), topology=(Periodic, Flat, Flat))
background_fields = Oceananigans.BackgroundFields(; background_closure_fluxes=true, b=B)
model = NonhydrostaticModel(; grid, background_fields)
time_step!(model, 1)
return true
end @test time_step_background_fields_with_closure_fluxes(arch)
@liuchihl, can you create a PR into Greg's branch that adds both of these tests to the test suite? |
Sure, I will work on the unit tests. |
And note that we only really need a unit test to merge this, same for many things. Correctness is a high bar and it's ok if we can't come up wtih someone right away. Sometimes we don't really have a way to do correctness and our best option is something like a regression test (eg we verified it worked at one point, so we just make sure that it keeps returning that same result). One way to go partway towards functional test but not all the way to "correctness" is to simply test that a simple set up returns a different result when closure fluxes are included vs not (for example). |
…ergence! Passed the WindowedTimeAverage test and seemed to be okay for the coarse-resolution internal tide test, but still need to run a high-resolution test case for a few tidal cycles, with picking up from a checkpoint, to be sure! See PRs: - CliMA/Oceananigans.jl#3721 - CliMA/Oceananigans.jl#3646
Is there still interest in this PR? |
@glwagner Thank you for following up. Yes, we’re still interested in this great PR. Apologies for the delay—I’m planning to work on the unit test later this month. |
No worries at all. Since you're interested I will try to keep it up to date. |
cc @hdrake @liuchihl