(0.97.5) Enforce mass-conservation inside domains with Radiation OpenBoundaryConditions

[tomchor]

Recently (I think in #4578) we found that mass wasn't being exactly conserved whenever PerturbationAdvection was used as an open boundary condition. This PR (which is a draft for now) explores adding a correction to normal velocities at the boundaries right before solve_for_pressure!() is called which enforces that no net mass enters or leaves the domain.

At the moment I just focused on something that works to get started and the algorithm is not optimized at all (suggestions on how to do that are welcome) but it works by:

1. Calculating the net mass flux into the domain by integrating the velocities at the boundaries (only for boundaries with an OpenBoundaryCondition with a matching scheme)
2. Dividing the correction equally between all the aforementioned boundaries
3. Adding/subtracting the same correction constant to left/right open boundaries

A MWE measuring the net mass flux into the domain produces this result without the correction:

[ Info: time: 0.000, max|u|: 1.011, CFL: 0.20, Net flux: 1.9306e-03
[ Info: time: 0.100, max|u|: 1.010, CFL: 0.20, Net flux: 4.1407e-03
[ Info: time: 0.200, max|u|: 1.009, CFL: 0.20, Net flux: 4.6124e-03
[ Info: time: 0.300, max|u|: 1.009, CFL: 0.20, Net flux: 4.4737e-03
[ Info: time: 0.400, max|u|: 1.008, CFL: 0.20, Net flux: 4.1404e-03
[ Info: time: 0.500, max|u|: 1.008, CFL: 0.20, Net flux: 3.7634e-03
[ Info: time: 0.600, max|u|: 1.008, CFL: 0.20, Net flux: 3.3962e-03
[ Info: time: 0.700, max|u|: 1.007, CFL: 0.20, Net flux: 3.0557e-03
[ Info: time: 0.800, max|u|: 1.006, CFL: 0.20, Net flux: 2.7461e-03
[ Info: time: 0.900, max|u|: 1.007, CFL: 0.20, Net flux: 2.4666e-03
[ Info: time: 1.000, max|u|: 1.006, CFL: 0.20, Net flux: 2.2152e-03

and this result with the correction:

[ Info: time: 0.000, max|u|: 1.011, CFL: 0.20, Net flux: 3.8326e-16
[ Info: time: 0.100, max|u|: 1.008, CFL: 0.20, Net flux: -2.7974e-18
[ Info: time: 0.200, max|u|: 1.007, CFL: 0.20, Net flux: 3.4705e-16
[ Info: time: 0.300, max|u|: 1.007, CFL: 0.20, Net flux: 5.3560e-16
[ Info: time: 0.400, max|u|: 1.006, CFL: 0.20, Net flux: -1.0820e-16
[ Info: time: 0.500, max|u|: 1.006, CFL: 0.20, Net flux: -1.6355e-16
[ Info: time: 0.600, max|u|: 1.006, CFL: 0.20, Net flux: -4.9130e-16
[ Info: time: 0.700, max|u|: 1.006, CFL: 0.20, Net flux: 2.5534e-16
[ Info: time: 0.800, max|u|: 1.005, CFL: 0.20, Net flux: -3.0518e-16
[ Info: time: 0.900, max|u|: 1.005, CFL: 0.20, Net flux: 4.2174e-16
[ Info: time: 1.000, max|u|: 1.005, CFL: 0.20, Net flux: 1.2774e-16

Notes:

1. This came about in the context of immersed boundaries, but no immersed boundary is necessary to produce these results. Flows with open (radiation) boundaries seem to be divergent even in regular grids, although the FFTSolver seems to handle this well enough that we haven't had problems.
2. @jagoosw I'm curious about your thoughts here. You mentioned before that this might be due to an asymmetry between left and right boundaries when time-stepping the code, but I don't see how the PerturbationAdvection matching scheme (where boundary velocities are computed locally) can be made to always provide velocities at the boundary that are consistent with zero mass flux without having some sort of nonlocal correction.

cc @glwagner @xkykai @ali-ramadhan

[glwagner]

The momentum flux is u*u*A but the mass flux is u*A --- we're concerned with mass flux, right?

[tomchor]

The momentum flux is uuA but the mass flux is u*A --- we're concerned with mass flux, right?

Yes, we're concerned with mass. I meant that I'm measuring velocity at the boundary. I removed the confusing statement from the comment.

[tomchor]

Also for the record, this does seem to allow for the ConjugateGradientSovler to converge when using ImmersedBoundaryGrids + PerturbationAdvectinoOpenBoundaryCondition at least on simple setups. As an example, this code (with a flat immersed boundary bottom) blows up immediately on main, but runs fine on this branch.

On main

[ Info: time: 0.000, max|u|: 1.010, CFL: 0.10, Net flux: 1.1067e-03
[ Info: time: NaN, max|u|: NaN, CFL: NaN, Net flux: NaN
[ Info: time: NaN, max|u|: NaN, CFL: NaN, Net flux: NaN
[ Info: time: NaN, max|u|: NaN, CFL: NaN, Net flux: NaN
[ Info: time: NaN, max|u|: NaN, CFL: NaN, Net flux: NaN
[ Info: time = NaN, iteration = 100: NaN found in field u. Stopping simulation.
[ Info: time: NaN, max|u|: NaN, CFL: NaN, Net flux: NaN

This branch:

...
[ Info: time: 48.497, max|u|: 1.000, CFL: 0.10, Net flux: 5.5511e-17
[ Info: time: 48.997, max|u|: 1.000, CFL: 0.10, Net flux: 4.1633e-17
[ Info: time: 49.497, max|u|: 1.000, CFL: 0.10, Net flux: 6.9389e-17
[ Info: time: 49.997, max|u|: 1.000, CFL: 0.10, Net flux: 9.7145e-17

[glwagner]

Should we make this PR not a draft anymore, or is there a reason to keep it as a draft? It seems like something we want to merge ASAP?

[tomchor]

Should we make this PR not a draft anymore, or is there a reason to keep it as a draft? It seems like something we want to merge ASAP?

I made it a draft because I spent no time optimizing it, so it's very slow (as you can imagine due to the allocations being created every time-step). So that was my attempt as signaling that it's not ready for a full review. I'm okay with making it a full PR though.

[glwagner]

Should we make this PR not a draft anymore, or is there a reason to keep it as a draft? It seems like something we want to merge ASAP?

I made it a draft because I spent no time optimizing it, so it's very slow (as you can imagine due to the allocations being created every time-step). So that was my attempt as signaling that it's not ready for a full review. I'm okay with making it a full PR though.

Are we ready to work on turning this into a proper implementation, or do we want more time to tinker? Happy to help whenever!

[tomchor]

Should we make this PR not a draft anymore, or is there a reason to keep it as a draft? It seems like something we want to merge ASAP?

I made it a draft because I spent no time optimizing it, so it's very slow (as you can imagine due to the allocations being created every time-step). So that was my attempt as signaling that it's not ready for a full review. I'm okay with making it a full PR though.

Are we ready to work on turning this into a proper implementation, or do we want more time to tinker? Happy to help whenever!

Gave it some thought and I think we're ready. The actual algorithm is pretty straightforward. All we need to figure out is the software design.

[jagoosw]

I read a bit more about this, and it turns out open boundary conditions often don't conserve mass, and a common way to correct is to scale the boundary values to correct so this seems reasonable.

[tomchor]

I read a bit more about this, and it turns out open boundary conditions often don't conserve mass, and a common way to correct is to scale the boundary values to correct so this seems reasonable.

Okay, that's good in the sense that it matches my understanding. Could you point to some of the literature you read though? I'm having some trouble finding relevant literature for these kinds of "details".

[tomchor]

These modifications seem to slow down the simulation by a lot. The culprit seems to be the average boundary mass flux calculations. Running this simulation, which has only one radiation open boundary (and therefore needs 1*3 calculations of the mass flux through the boundaries per RK3 time-step) leads to almost a factor ~1.8 slowdown.

Having 4 radiation open boundaries leads to a slowdown of ~6 times...

Not sure what can be done to improve this

(I'm only doing this calculation wherever there is a radiation open boundary (i.e. an OpenBoundaryCondition with a matching scheme))

[glwagner]

I think this needs a good proofread because I found a few details to fix. Please read through. Reduce verbosity if possible, very long lines / variable names can impact the legibility of code

[tomchor]

Addressed all your comments and gave it a thorough proofread afterwards, as per your suggestion.

[tomchor]

#4668 seems to have really sped up things. Here are the benchmarks for a simulation without and then with two open boundaries before #4668 was merged:

julia> b_without_open
BenchmarkTools.Trial: 5785 samples with 1 evaluation per sample.
 Range (min … max):  662.920 μs …   8.982 ms  ┊ GC (min … max): 0.00% … 85.53%
 Time  (median):     871.306 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   861.066 μs ± 559.721 μs  ┊ GC (mean ± σ):  4.45% ±  6.25%

       ▂▄▅▅▃                         ▁▄█▆▄▂▁                     
  ▁▂▃▄▇█████▇▆▅▄▃▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▂▃▄████████▆▅▄▃▃▃▃▂▃▂▂▂▂▂▂▁▁▁▁ ▃
  663 μs           Histogram: frequency by time         1.02 ms <

 Memory estimate: 731.53 KiB, allocs estimate: 2104.

julia> b_with_open
BenchmarkTools.Trial: 5104 samples with 1 evaluation per sample.
 Range (min … max):  754.581 μs …  12.483 ms  ┊ GC (min … max): 0.00% … 89.36%
 Time  (median):     972.037 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   975.329 μs ± 620.934 μs  ┊ GC (mean ± σ):  4.54% ±  6.48%

        ▄███▆▄▁                 ▂▅▇▇▇▇▇▃                         
  ▂▂▃▄▆████████▇▅▅▄▃▃▃▃▂▂▂▂▂▃▄▅▇█████████▆▅▄▄▅▃▃▃▃▃▃▃▃▂▂▃▂▂▂▁▂▂ ▄
  755 μs           Histogram: frequency by time          1.2 ms <

 Memory estimate: 746.67 KiB, allocs estimate: 2382.

and here are the results after it was merged:

julia> b_without_open
BenchmarkTools.Trial: 6946 samples with 1 evaluation per sample.
 Range (min … max):  532.364 μs …   5.758 ms  ┊ GC (min … max): 0.00% … 78.44%
 Time  (median):     626.273 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   716.930 μs ± 406.932 μs  ┊ GC (mean ± σ):  4.45% ±  6.98%

          ▁▅██▆▂                           ▁▃▂                   
  ▁▁▂▂▂▂▃▅███████▆▄▃▃▂▂▂▂▁▁▁▁▁▁▁▂▂▂▂▂▂▃▃▅▇████▇▅▃▂▂▂▂▁▂▁▂▁▁▁▁▁▁ ▃
  532 μs           Histogram: frequency by time          910 μs <

 Memory estimate: 715.31 KiB, allocs estimate: 1776.

julia> b_with_open
BenchmarkTools.Trial: 6801 samples with 1 evaluation per sample.
 Range (min … max):  539.551 μs …   8.671 ms  ┊ GC (min … max): 0.00% … 88.88%
 Time  (median):     661.378 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   731.322 μs ± 411.323 μs  ┊ GC (mean ± σ):  4.48% ±  7.13%

       ▃██▄▂▁                   ▁                                
  ▂▂▃▄▇██████▇▅▄▅▄▄▄▅▄▄▃▃▃▃▃▃▄▆███▆▄▄▃▃▃▂▂▃▃▂▃▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁ ▃
  540 μs           Histogram: frequency by time         1.04 ms <

 Memory estimate: 723.58 KiB, allocs estimate: 1844.

The difference is now much smaller and, perhaps more importantly, both are faster 🎉