Oceananigans.jl
Oceananigans.jl copied to clipboard
CliMA
Moves masking fields from model update to output writers when using ImmersedBoundary
This PR removes the masking of tracers when using IBM.
At the moment all tracers are masked inside immersed solids. It's my understanding that these masked values aren't being used in any calculations except for the pressure solver.
~The issue is that having active scalars masked can mess up the pressure solver and lead to some numerical artifacts in some cases. I've noticed some small effects like that in my simulation, but @whitleyv had some more serious spurious dynamics pop up as a result in her simulations that only went away by not masking buoyancy. I'll let @whitleyv herself explain/show what was happening.~
~For now @whitleyv and I just commented out the lines doing the masking in the nonhydrostatic model, but if we do decide to move forward with this PR I'll clean the code more properly. I was hoping we could discuss a bit before doing that.~
Updated description:
When I originally opened this PR, the masking was creating an issue with the dynamics. However, this turned out to be a bug and was solved by #2603. Now the current situation is that masking is still done for all tracers and velocities in the main branch, but this has zero effect on the results.
This PR removes masking when updating NonhydrostaticModel and instead masks each output field only when writing to disk, saving computational resources (although just a bit).
CC @simone-silvestri @glwagner @wenegrat
If tracer masking changes the pressure solver, it is likely a bug. The pressure solver finds the solution to the Poisson pressure equation, which is a Poisson equation with a RHS proportional to the divergence of the predictor velocity.
Tracer masking should have no effect on dynamics and is purely a convenience for output (so we should probably move the masking step to output writers to avoid the unnecessary computational expense). However, there is a bug with immersed boundaries and high-order advection schemes (solved by #2603) that means the value of tracer fields in the immersed boundaries does (spuriously and incorrectly) influence the outcome of a simulation.
In addition to that, interpolation / reconstruction across immersed boundaries can in principle touch tracer values within immersed boundaries. However, as far as I know we do not ever reconstruct across immersed boundaries in any core routines except the experimental IsopycnalSkewSymmetricDiffusivity. #2477 solves this problem by implementing "homogeneous" operators that do not touch tracer values within the immersed boundaries (ie differences return 0, reconstruction returns values from "active" cells only).
This PR could move tracer masking to output writing to save the computation. However, if the tracer masking affects the outcome of a simulation, then there is a deeper issue at hand (perhaps the issue that's solved by #2603?)
Perhaps a more specific question is: why does removing the tracer masking step improve simulation outcomes?
In the past, problems seem to arise when masking buoyancy as a tracer. Here is an example of a slope at rest with a uniform stratification in a 2D domain run on the main branch. This is using the nonhydrostatic model, but I've plotted the hydrostatic pressure and its horizontal derivative because it shows how masking buoyancy with a strong stratification results in really strong horizontal pressure gradients. These gradients seem to leak out along the boundary creating the numerical effects you can see in the velocity.
https://user-images.githubusercontent.com/67593861/175611052-1aa81aa0-0aa1-4e9a-bfbd-85369b07d99a.mp4
When I did this analysis a while ago, this problem seemed to be fixed by unmasking the tracers, but on the most up to date version of Oceananigans, it doesn't seem to fix the problem. As @glwagner points out, that makes it sound more like a bug. Running it with https://github.com/CliMA/Oceananigans.jl/pull/2603, seems to fix the issue, at least the issue I noticed before. I believe @tomchor has tried both https://github.com/CliMA/Oceananigans.jl/pull/2603 and unmasking and seemed to see a slight improvement with unmasking?
Okay, thank you very much for the clarification about hydrostatic vs. nonhydrostatic pressure! The hydrostatic pressure integral:
https://github.com/CliMA/Oceananigans.jl/blob/e76deefc1fdf115e05856a18977c593d0b5d3b0d/src/Models/NonhydrostaticModels/update_hydrostatic_pressure.jl#L15
does indeed use tracers; in addition it involves a reconstruction / interpolation that might be affecting results here. Bug perhaps? I say "might" and "perhaps" because the precise mechanism isn't clear to me. Since we integrate downwards, there's only an error in the hydrostatic pressure below the immersed boundary. Hmm...
Note, this particular issue is not fixed by "unmasking" tracers (it's just that the bug / error is different, and perhaps less innocuous). This bug is also not fixed by #2603.
I m wondering if this issue is solved in the new version
It is! Although I'll move the masking from the model update to the output construction since it doesn't change model runs and it saves time.
Tests should be passing, although apparently the GPU runs aren't running currently
@navidcy @simone-silvestri can someone review this PR? I believe it may be ready to merge.
I don't know the details here to provide useful feedback. Sorry... :(
That's a good point. This PR changed from when it was first introduced and it's much simpler in scope, so I updated the PR description. If that's enough for you to give me feedback, feel free to do so :)
I'm only modifying the NonhydrostaticModel for now, as I'm not familiar with the hydrostatic one. (Not sure if eventually the same thing should be that for that...)
We should make the changes needed for both nonhydrostatic and hydrostatic models.
We should make the changes needed for both nonhydrostatic and hydrostatic models.
For sure.
We should make the changes needed for both nonhydrostatic and hydrostatic models.
@glwagner the reason why I hadn't removed masking from the hydrostatic model is that I'm really not familiar with it at all. So it's harder for me to figure out where it's okay to remove it, and testing it.
I just did a best guess and removed it from shallow water and nonhydrostatic models in the places I thought appropriate. It would be good to have some feedback from someone more familiar with the code though on that, though.
@tomchor, if you want I can remove it from the hydrostatic model; it should be the same as for the nonhydrostatic.
Actually, I saw that you already removed it. I am wondering though if we don't need to keep the masking for the velocities, as those need an "impenetrability" condition
We have to keep masking for velocities in the hydrostatic model.
In the nonhydrostatic model, the equivalent masking procedure acts on the predictor velocity field prior to solving for pressure:
https://github.com/CliMA/Oceananigans.jl/blob/ea6826fd2ffaed3f0df330cc952667ec878deb6a/src/Models/NonhydrostaticModels/pressure_correction.jl#L12-L14
We have to keep masking for velocities in the hydrostatic model.
In the nonhydrostatic model, the equivalent masking procedure acts on the predictor velocity field prior to solving for pressure:
Yes, I kept that in the nonhydrostatic model. Although I think I removed it on the hydrostatic one so I'll put it back in.
@simone-silvestri @glwagner I think tests should be passing now. Ready for reviews
Also, I don't understand why we're calling masking events at update for every NonhydrostaticModel regardless of it having an ImmersedGrid or not:
https://github.com/CliMA/Oceananigans.jl/blob/d4a45adf21918f95ca4d23ec7167e720c44c501b/src/Models/NonhydrostaticModels/update_nonhydrostatic_model_state.jl#L14-L19
The same goes for the hydrostatic model
Why not something like
if model.grid isa ImmersedBoundaryGrid
# call masking
end
?
That is what the code does, because
https://github.com/CliMA/Oceananigans.jl/blob/d4a45adf21918f95ca4d23ec7167e720c44c501b/src/ImmersedBoundaries/mask_immersed_field.jl#L8
Note that fetch_and_convert_output is called here:
https://github.com/CliMA/Oceananigans.jl/blob/cb8ac3e0bab2c25ab6c09afd6e73356ce86406a3/src/OutputWriters/jld2_output_writer.jl#L259-L260
for JLD2OutputWriter, and here:
https://github.com/CliMA/Oceananigans.jl/blob/cb8ac3e0bab2c25ab6c09afd6e73356ce86406a3/src/OutputWriters/netcdf_output_writer.jl#L420
for the NetCDF output writer.
@glwagner I tested this locally, tests are passing, and I believe I addressed all your comments, so this should be ready for re-review. That said, github tells me there's a requested revision by you but I can't figure out what to do about it. Could you please take a look?
@simone-silvestri @navidcy also feel free to review
I think we will need to do some more testing to make sure this won't break our near global hydrostatic setups, which are unfortunately in a tenuous position because they rely on some untested features (and we don't have regression tests for some important cases). @simone-silvestri what do you think? We may want to wait for a few more PRs (perhaps containing some of those tests) to go in first.
I think we will need to do some more testing to make sure this won't break our near global hydrostatic setups, which are unfortunately in a tenuous position because they rely on some untested features (and we don't have regression tests for some important cases). @simone-silvestri what do you think? We may want to wait for a few more PRs (perhaps containing some of those tests) to go in first.
@glwagner what is the status of the global hydrostatic setup tests? Do you think we've reached a point where we can merge this PR?
I don't think we have regression tests yet but not sure, @simone-silvestri would know
@tomchor I think we should revisit this PR. Rather than "moving" the masking, I think we should start by adding a feature to both output writers that defines a "mask value" for immersed regions. What do you think about that? I would like to open a new PR for that. And before that, I'm going to open an issue so we can discuss the user interface.
Once we have a nice feature for masking output, we can consider whether or not to also mask (or not) within the time-stepping. These are separate questions, I think.
If you agree with that then we can close this PR and I will open an new issue to discuss the user API.