Compare ImmersedBoundary with tilted geometry

[francispoulin]

To better understand Wenegrat and Thomas, JPO 2020 (WT2020) I modifed the following examples:

1. Internal tide by a seamount and

2. Tilted bottom boundary layer example ,

to study the dynamics of constant flow over a linearly slowing bottom. Note that I tried to pick the parameters as close to one of the cases in WT2020, however I did make some modifications.

1. I changed the direction of the constant horizontal flow (along bathymetry) to avoid inertial instability. This produces simpler dynamics.
2. The code uses free-slip boundary conditons by default, but I believe I have no slip working as well.
3. Currently diffusion and viscosity is larger by a factor of 100 as it's easier to resolve the boundary layer that develops.

The examples were very easy to modify, which shows that people have done a great job in making them friendly, as promised. Thanks everyone for that! You can find my examples on branch fjp/Test_ImmersedBoundaryMethod.

The first code, what I call take1, uses the immersed boundary method.

Thet second code, what I call take2 used the tilted geometry approach.

Below I will share some results that people might find interest, maybe @tomchor in particular?

The first plot shows at heatmap of the velocity profile and a contour plot of the isopycnals for the case of an immersed boundary (take1). We see the isopycnals incrop (reach the bathymetry) at a 90 degree angle, which is consistent with what is said in the paper. It seems odd that slip boundary conditions yields a flow that gets faster over the topography, but that is not a problem.

The second plot shows at heatmap of the velocity profile and a contour plot of the isopycnals for the case of a tilted geometry (take2). This case generates a boundary layer at the top as well, but it's less well resolved because I am using a stretched grid at the bottom. Not a problem, just an observation. This also shows the speed up over the topography, so it must be physical.

There is good qualitative agreement between these two and I'm happy with the results.

Next you will find animations for the two cases where we show v, w b, in that order. First, for the immersed boundary. What's note worthy is that there are relatively strong vertical velocities that occur over the immersed boundary.

linear_bathymetry_slip_take1.mp4

This is the same simulation but for the case of a tilted geometry, and we don't see the same spikes, equally spaced based on the immersed boundary.

linear_bathymetry_slip_take2.mp4

Is the vertical velocity generated over the piecewise constant topogrpahy something that is of a concern? When I increased the resolution I believe I saw the vertical velocity generated increased as well.

Any other thoughts about this comparison?

If people want to see the results with no slip boundary conditions, just let me know and I'm happy to share those too.

[glwagner]

Is the vertical velocity generated over the piecewise constant topogrpahy something that is of a concern?

I think so! Is this with GridFittedBottom or PartialCellBottom? I feel like this motivates cut / shaved cells.

Should we convert this to a discussion?

[tomchor]

Is the vertical velocity generated over the piecewise constant topogrpahy something that is of a concern? When I increased the resolution I believe I saw the vertical velocity generated increased as well.

I agree with @glwagner that it's concerning at least for your current configuration. My guess is that with bottom drag and a more turbulent simulation that might not be a huge problem. But as it stands it's definitely affecting the dynamics.

(Btw, I say that because I've run a fair amount of turbulent simulations with immersed boundaries at this point and I've never seen this being an issue in my simulations. But they are always turbulent.)

I think so! Is this with GridFittedBottom or PartialCellBottom? I feel like this motivated cut / shaved cells.

@glwagner reading the code here it seems that it's a GridFittedBottom.

@francispoulin Like Greg mentioned, shaved cells would be ideal here, but since we don't have that yet, I'd suggest trying out the PartialCellBottom. I expect some improvement there.

[hdrake]

I did almost the exact comparison as @francispoulin has done here early on in my PhD, but using the MITgcm: comparing flat-bottom tilted-geometry vs. tilted-bottom.

I also want to point out that the circulation actually is qualitatively different between the flat-bottom , because of the nature of the bounded vs. periodic lateral boundary conditions (see MITgcm simulations in Ruan et al. 2020 for the untilted case and Drake et al. 2022 for the tilted case) and the inluence of barotropic pressure constraints (see Peterson and Callies 2022). So I would not read too much into the differences between the two configurations–they are mostly due to other factors besides the cell steppiness.

By the way, there's nothing stopping you from putting topography on top of a tilted bottom! @liuchihl and I have successfully done this, both in MITgcm and in Oceananigans (building off of @tomchor's tilted BBL example).

In any case, my MITgcm experiments demonstrated that the shaved cells can help a lot to reduce the grid-scale noise in the vertical velocity.

[tomchor]

I also want to point out that the circulation actually is qualitatively different between the flat-bottom , because of the nature of the bounded vs. periodic lateral boundary conditions

I agree with that, but just for clarification, I think the tilted set-up shown here also uses bounded lateral BCs. @francispoulin can you confirm?

[hdrake]

Oh, looks like you're right @tomchor. I guess this also addresses the barotropic pressure constraint issue as well, but I'd need to think more about it. I still would not expect the solutions to exactly converge (e.g. because the upper boundary condition is still slightly different), but I guess they might be pretty close.

[francispoulin]

Yes, I did use bounded walls on the sides. It seems natural and easy.

[hdrake]

I think that's a sensible approach and wish I'd thought of that back in the day!

[glwagner]

We've seen issues noise in the vertical velocity along immersed boundaries in other configurations too. I've seen at least one result in which the noise is substantially mitigated by using the PCG solver, which avoids the approximations of the "naive" FFT solver. Unfortunately, we don't yet have a performant PCG-based solver cc @simone-silvestri @xkykai .

I think it's interesting and also convenient that the noise is mitigated in turbulent cases or by the inclusion of bottom drag.

[navidcy]

This is an interesting discussion!

(Thus, I think we should convert this into a "Discussion".)

[francispoulin]

Good suggestion @glwagner (and second by @navidcy) and I now know how to convert issues to discussions! I'll try and plan things better in the future.

As has been said, I am using GridFittedBottom but I will try PartialCellBottom and let you know what I find.

I agree with @tomchor that when you add bottom drag, this will likely kill the spurious vertical velocities that are generated, but I hope we can simulate this without the drag.

I will try and give an update soon after I relearn how partial cells work.

Cheers!

[francispoulin]

I am trying to run resting_stratified_bumpy_ocean.jl to relearn how to use PartialCellBottom, unfortunately, it doesn't run. There were some minor fixes I was able to make but there is a bigger problem that I don't understand. I think this has to do with there being a Flat dimension as I can define the model in heterogeneous_windy_convection.jl, but that's what fails here.

I updated it on this branch but am including the error message below. Any suggestions on how we can fix this?

julia> include("resting_stratified_bumpy_ocean.jl")
┌ Warning: Assignment to `grid` in soft scope is ambiguous because a global variable by the same name exists: `grid` will be treated as a new local. Disambiguate by using `local grid` to suppress this warning or `global grid` to assign to the existing global variable.
└ @ ~/Software/TMP/Oceananigans.jl/validation/immersed_boundaries/resting_stratified_bumpy_ocean.jl:41
┌ Warning: Assignment to `model` in soft scope is ambiguous because a global variable by the same name exists: `model` will be treated as a new local. Disambiguate by using `local model` to suppress this warning or `global model` to assign to the existing global variable.
└ @ ~/Software/TMP/Oceananigans.jl/validation/immersed_boundaries/resting_stratified_bumpy_ocean.jl:45
┌ Warning: Assignment to `N²` in soft scope is ambiguous because a global variable by the same name exists: `N²` will be treated as a new local. Disambiguate by using `local N²` to suppress this warning or `global N²` to assign to the existing global variable.
└ @ ~/Software/TMP/Oceananigans.jl/validation/immersed_boundaries/resting_stratified_bumpy_ocean.jl:52
grid = 1×128×64 ImmersedBoundaryGrid{Float64, Flat, Periodic, Bounded} on CPU with 0×3×3 halo:
├── immersed_boundary: PartialCellBottom(mean(z)=-0.894182, min(z)=-1.0, max(z)=-0.500488, ϵ=0.2)
├── underlying_grid: 1×128×64 RectilinearGrid{Float64, Flat, Periodic, Bounded} on CPU with 0×3×3 halo
├── Flat x
├── Periodic y ∈ [-1.0, 1.0)      regularly spaced with Δy=0.015625
└── Bounded  z ∈ [-1.0, 0.0]      regularly spaced with Δz=0.015625
ERROR: LoadError: BoundsError: attempt to access Tuple{Float64, Float64} at index [3]
Stacktrace:
  [1] indexed_iterate
    @ ./tuple.jl:92 [inlined]
  [2] Δzᶜᶜᶜ
    @ ~/Software/TMP/Oceananigans.jl/src/ImmersedBoundaries/partial_cell_bottom.jl:111 [inlined]
  [3] Δzᶠᶜᶜ
    @ ~/Software/TMP/Oceananigans.jl/src/ImmersedBoundaries/partial_cell_bottom.jl:139 [inlined]
  [4] getindex
    @ ~/Software/TMP/Oceananigans.jl/src/AbstractOperations/grid_metrics.jl:138 [inlined]
  [5] getindex
    @ ~/Software/TMP/Oceananigans.jl/src/AbstractOperations/conditional_operations.jl:96 [inlined]
  [6] _getindex
    @ ./abstractarray.jl:1338 [inlined]
  [7] getindex
    @ ./abstractarray.jl:1288 [inlined]
  [8] macro expansion
    @ ./reducedim.jl:309 [inlined]
  [9] macro expansion
    @ ./simdloop.jl:77 [inlined]
 [10] _mapreducedim!(f::typeof(identity), op::typeof(Base.add_sum), R::SubArray{…}, A::Oceananigans.AbstractOperations.ConditionalOperation{…})
    @ Base ./reducedim.jl:308
 [11] mapreducedim!(f::Function, op::Function, R::SubArray{…}, A::Oceananigans.AbstractOperations.ConditionalOperation{…})
    @ Base ./reducedim.jl:324
 [12] sum!(f::Function, r::SubArray{…}, A::Oceananigans.AbstractOperations.ConditionalOperation{…})
    @ Base ./reducedim.jl:1034
 [13] sum!(r::Field{…}, a::Oceananigans.AbstractOperations.GridMetricOperation{…}; condition::Nothing, mask::Int64, kwargs::@Kwargs{})
    @ Oceananigans.Fields ~/Software/TMP/Oceananigans.jl/src/Fields/field.jl:683
 [14] sum!(r::Field{…}, a::Oceananigans.AbstractOperations.GridMetricOperation{…})
    @ Oceananigans.Fields ~/Software/TMP/Oceananigans.jl/src/Fields/field.jl:677
 [15] Oceananigans.Models.HydrostaticFreeSurfaceModels.SplitExplicitAuxiliaryFields(grid::ImmersedBoundaryGrid{Float64, Flat, Periodic, Bounded, RectilinearGrid{…}, PartialCellBottom{…}, Nothing, Nothing, CPU})
    @ Oceananigans.Models.HydrostaticFreeSurfaceModels ~/Software/TMP/Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface.jl:234
 [16] materialize_free_surface(free_surface::SplitExplicitFreeSurface{…}, velocities::@NamedTuple{…}, grid::ImmersedBoundaryGrid{…})
    @ Oceananigans.Models.HydrostaticFreeSurfaceModels ~/Software/TMP/Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface.jl:121
 [17] HydrostaticFreeSurfaceModel(; grid::ImmersedBoundaryGrid{…}, clock::Clock{…}, momentum_advection::Centered{…}, tracer_advection::Centered{…}, buoyancy::BuoyancyTracer, coriolis::FPlane{…}, free_surface::SplitExplicitFreeSurface{…}, tracers::Symbol, forcing::@NamedTuple{}, closure::Nothing, boundary_conditions::@NamedTuple{}, particles::Nothing, biogeochemistry::Nothing, velocities::Nothing, pressure::Nothing, diffusivity_fields::Nothing, auxiliary_fields::@NamedTuple{})
    @ Oceananigans.Models.HydrostaticFreeSurfaceModels ~/Software/TMP/Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl:177
 [18] top-level scope
    @ ~/Software/TMP/Oceananigans.jl/validation/immersed_boundaries/resting_stratified_bumpy_ocean.jl:45
 [19] include(fname::String)
    @ Base.MainInclude ./client.jl:489
 [20] top-level scope
    @ REPL[58]:1
 [21] top-level scope
    @ ~/.julia/packages/CUDA/htRwP/src/initialization.jl:206
in expression starting at /home/fpoulin/Software/TMP/Oceananigans.jl/validation/immersed_boundaries/resting_stratified_bumpy_ocean.jl:40
Some type information was truncated. Use `show(err)` to see complete types.

[navidcy]

In src/ImmersedBoundaries/partial_cell_bottom.jl, change

x, y, z = node(i, j, k+1, underlying_grid, c, c, f)

to

z = node(i, j, k+1, underlying_grid, c, c, f)

From a quick look I had, x, y are not used anyways....

[navidcy]

wait, which branch? I can have a look..

[navidcy]

seems like z and h are not both Floats?

@show z

@show h

will reveal the truth :)

[francispoulin]

fjp/Test_ImmersedBoundaryMethod

[navidcy]

z = znode(i, j, k+1, underlying_grid, c, c, f)

not

z = node(i, j, k+1, underlying_grid, c, c, f)

!

(node returns a Tuple... znode returns the Float value of z)

[navidcy]

I pushed some edits (9b03066 and 706c80c)
have a look now

[francispoulin]

Thanks @navidcy for finding the fix and fixing it!

(edited)

Sorry that I said things were running before. I had mistakenly copied your version of my script, which was using GridFittedBottom.

Since the problem seems to be subtracting a tuple and a real I replaced - with .- and it runs!

[francispoulin]

Thanks to @navidcy I was able to get the parital cell method working. Below are the new results.

First, this is a plot of the final velocity and buoyancy fields. Unfortunately, it's much messier at the bottom, when compared with GridFittedBottom.

The simulation shows that the vertical velocities on the bottom seem to pulsate, unlike the previous case where they were localized. Different but not sure if this is any better.

linear_bathymetry_PartialCell_take1.mp4

I remember that last year, when talking to @glwagner , the partial cell method was created but not validated. Maybe there is still something that needs to be fixed? I'd be happy to work on this with someone if there was interest in improving it.

[glwagner]

I don't think PartialCellBottom has been tested much.

To me this points towards there being something wrong with PartialCellBottom since its whole point is to make the jumps smoother.

I think this is a good hypothesis, but it's not certainly true. The problem could be with the pressure solver / near-boundary divergence and thus could be relatively unaffected by the differences between the two immersed boundary methods.

PartialCellBottom is really simple and there are no mysteries. Here's the implementation:

https://github.com/CliMA/Oceananigans.jl/blob/main/src/ImmersedBoundaries/partial_cell_bottom.jl

The main subtlety is how we treat the integral for hydrostatic pressure. In order to satisfy the "resting stratified ocean remains at rest" criteria, we do not use the partial metrics in the hydrostatic pressure integral:

Oceananigans.jl/src/Models/NonhydrostaticModels/update_hydrostatic_pressure.jl

Lines 27 to 32 in fc871ef

	# Partial cell "algorithm"
	const PCB = PartialCellBottom
	const PCBIBG = ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:Any, <:PCB}
	update_hydrostatic_pressure!(pHY′, arch, ibg::PCBIBG, buoyancy, tracers; parameters = p_kernel_parameters(ibg.underlying_grid)) =
	update_hydrostatic_pressure!(pHY′, arch, ibg.underlying_grid, buoyancy, tracers; parameters)

So, if there are problems with horizontal gradients of hydrostatic pressure, those problems will be similar between grid fitted and partial cell bottoms. I guess this is yet another reason to omit the hydrostatic pressure integral for NonhydrostaticModel...

[francispoulin]

I can return to our validation script and rerun it to see whether it generates any spurious modes over the topography. If yes, then we will have a better clue as to what's going on.

I will try and look at this later today and let you know what I find.

[glwagner]

The simulation shows that the vertical velocities on the bottom seem to pulsate, unlike the previous case where they were localized.

To me it looks like there is a large scale internal wave that is driving those oscillations. It's not necessarily a bad thing / sign of an issue.

Huh, that's weird. I count 18 w "bottom sources" on the video with GridFittedBottom, and 17 on the video with PartialCellBottom. So they're distributed approximately at the same locations, and with the same distance from each other.

The criteria that determines which cells are immersed for PartialCellBottom is different than for GridCellBottom, so the shape / total volume of the domain can be slightly different between the two.

All in all I'm not sure there's a reason to suspect an issue with PartialCellBottom. I think it's more likely that the artifacts in the vertical velocity are either 1) present whenever a discontinuous representation is used (eg perhaps cut cells will solve this) or 2) associated with the near-boundary divergence. Note that 2) isn't easily dismissed; the boundary normal velocity is of course tightly associated with the near boundary divergence.

[glwagner]

@francispoulin @tomchor I'm not sure the code still works... but you could try running this case with the pressure solver that's introduced on #3188. If you find issues, you could report them on #3188 since it will help us move forward with that PR...

[hdrake]

I agree that the domain-scale pulsating is due to a transient initial adjustment. I saw a lot of this in bounded domains with a tilted bottom, and was one of the reasons it was much cleaner to just go into the tilted geometry coordinate system (w/ periodic cross-slope and along-slope BCs).

The spikes are definitely due to the cell discretization near the bottom though.

[navidcy]

PR #3530 is related

[francispoulin]

@glwagner The fact that we don't get strong spikes of vertical velocities in the case of a tilted gravity suggests to me that these spikes are not physical. But it would be nice to have a more convincing argument.

I modified resting_stratified_bumpy_ocean.jl so that it runs and pushed the changes. Below is a plot of the buoyancy fields for each of the two methods and the difference. For some reason the plot of velocity is not working, but I will try and get that working tomorrow.

It is easy to add a constant velocity in the orthogonal direction here but I think we want a simpler test. Any suggestions?

[glwagner]

Hmm, it looks like this case is behaving correctly. The tilted bottom case you set up in the original post is very revealing though, I think we should keep investigating that.

[simone-silvestri]

@glwagner The fact that we don't get strong spikes of vertical velocities in the case of a tilted gravity suggests to me that these spikes are not physical. But it would be nice to have a more convincing argument.)

In my opinion, the vertical velocity spikes are physical, most likely the effect of the steps and the boundary conditions. Or at least they would be physical if a flow would move parallel to a steppy bathymetry like that subject to those boundary conditions.
The physical domain is quite different in the immersed case vs the tilted case. The immersed case approaches the tilted case for infinite resolution.
Also, the boundary conditions are numerically different in the immersed boundary case where the "no-slip" boundary condition is modeled with a ValueBoundaryCondition for the y-velocity on top of the steps combined with a no penetration on the vertical side of the step, compared to the tilted domain where no-slip is modeled with an actual no-slip condition for the y-velocity that is directed along the sloping boundary.
It seems quite plausible to me that the combination of the two boundary conditions in the former case leads to a localized divergence of the y-velocity that generates a vertical flow. The same should not necessarily happen in a tilted domain where the y-velocity changes are smooth.

[glwagner]

Probably by "physical", @francispoulin may have meant "reflecting the solution in the case of infinite resolution"...

@francispoulin you can perhaps gain some clarity by running with different resolutions (refining both the vertical and and the horizontal resolution).

I think @simone-silvestri is thinking about the problem correctly though --- we have to interpret the discrete solution and the steps are part of that. I would only add the nuance that the magnitude of the vertical velocity features could depend on the pressure solver (something we have seen before).

[francispoulin]

Thanks @simone-silvestri and @glwagner for your responses.

I agree that flow over bathymetry that is piecewise constant is very different than a linear slope, and as the resolution increases they should converge. I'll do some tests now and see how that goes.

I will add that with the partial cell method, my understanding is we are making the bathymetry piecewise linear. In that case, it should be able to do a linear bathymetry very well, however we still get strong vertical velocities. I'll investigate this as well.

More results soon, I hope!

[hdrake]

I think seeing how the immersed boundary configuration converges with higher horizontal / vertical resolution would be useful.

As mentioned above, be careful interpreting differences beyond the vertical velocity spikes: the immersed boundary configuration will not approach the exactly approach the tilted geometry configuration even with infinite resolution!

[francispoulin]

Below are plots of the vertical velocity in the simulations of the GridFittedBottom cases after 2 hours for different resolutions. Unfortunately, in these simulations the vertical velocities are not going away and in the last simulation they increase by an order of magnitude. Indeed, the last case has a spike that forms, which I have seen previously.

I will now look at the partial cell method but thoughts would be greatly appreciated?

1. 100x100 with a maximum w of 3.4e-4.

2. 200x200 with a maximum w of 3.6e-4.

3. 400x400 with a maximum w of 3.5e-4.

4. 800x800 with a maximum w of 2.9e-3!!!.

linear_bathymetry_take1.mp4

[tomchor]

For the record: I've seen this issue of a random velocity "spike" that happens in some simulations, but only at specific resolutions. In my simulations that happens because, due to the immersed criterion and depending on the resolution, two "steps" can be in very close proximity to each other, causing flow constriction that accentuates any small divergences that happen near the immersed boundary.

@glwagner @xkykai idk if you remember but I discussed that previously with you guys on slack, although those messages have been lost since we have a free account. In my case this issue went away when I removed the hydrostatic pressure separation.

[glwagner]

Another reason to remove the separation

[francispoulin]

Below are some results for the PartialCellBottom and they look worst, maybe not surprising as this hasn't been tested very carefully.

The irregularity of the buoyancy contours near the bottom are a cause for concern.

1. 100x100 with a maximum w of 9.4e-4.

2. 200x200 with a maximum w of `1.1e-3'.

[hdrake]

Following up, I was looking back at some of my MITgcm analysis comparing the rotated and non-rotated configurations and I have never seen spikes in vertical velocity this large. The caveat is that my simulations were all much lower aspect ratios and therefore run in hydrostatic mode with dx=500 m and dz=5 m.

Here is a notebook comparing simulations with sinusoidal hills (in x) on top of either a flat bottom or a tilted-geometry bottom, both using MITgcm's partial cell implementation.

The figure at the bottom is the main result (left untilted vs. right tilted; both periodic in x and y):

A minor caveat is that these simulations were 3D, but at this point in the spinup have not yet become turbulent and are essentially uniform in y.

And here is a notebook comparing the flat-bottom case (analagous to your experiments, except the tilted one is periodic in x while the untilted one is bounded in x):

The noise in the vertical velocity you see here in the untilted case is due to how the noise I added to the initial temperature field impacts hydrostatic pressure, not the numerical adjustment near bottom cells.

[glwagner]

@hdrake that is super useful! It might be nice to have both cases in hand, one hydrostatic with thin cells, more traditional, and another nonhydrostatic with O(1) aspect ratio cells, which maybe the ocean community has less experience with but we think will likely be important in the future.

Perhaps @hdrake's results suggest that improvements to PartialCellBottom are the thing to focus on in the near term, not sure (I think we would to see the exact same case with Oceananigans to be sure)

[hdrake]

It would definitely be useful to see these vertical velocity spikes are affected by grid aspect ratio and hydrostatic vs. non-hydrostatic!

Overall, I think PartialCellBottom will only make a difference on the margins (when topographic slopes are comparable to the grid aspect ratio, not much greater or much less). In global hydrostatic z-coordinate ocean models (e.g. MITgcm or MOM5) with grid cell aspect ratios dz/dx << 1, you inevitably will have steppy vertical velocities along basically any topographic features (even large-scale ones like mid-ocean riges and continental shelves). I've never seen it quite as bad as in @francispoulin's animations, however, so there could still be a bug here or else this is just a feature of O(1) aspect ratio simulations.

[francispoulin]

@hdrake , your notebooks are great for plotting but don't have much information in terms of what parameters were used in the simulations. Is there somewhere else where I could learn more about your configuration?

[hdrake]

Sorry; these runs were abandoned a long time ago so the setup files don't live on github but I've been able to track them down and will try to upload them when I can.

[francispoulin]

That is very kind of you and please take your time.

[francispoulin]

I did some simulations with GridFittedBottom for the same resolutions as before but now with no slip boundary conditions. Below are the plots for the vertical and horizontal velocities in the finest case with 800x800.

The good news is that he horizontal velocity looks smooth and do go to zero by the bottom.

The not so good news is that the spurious vertical velocities persist.

Things that I tried that didn't do much are the following:

1. change aspect ratio
2. increase the vertical diffusivities
3. change advection schemes
4. use nonhydrostatic model.

Next, I might try adding some damping like @tomchor has done in the tilted gravity example to see what effect that has.

[simone-silvestri]

I have done some thinking and I maintain that those "spikes" are completely consistent with the discretized equations and a consequence of (1) the vertical grid (2) the bathymetric representation.

In this setup, the flow is pushed uphill (against a buoyancy gradient but that is not the point). It's only logical that the flow has to rise somewhere. The absence of a positive vertical velocity somewhere would be unphysical, not the presence of it.
Now, the question is: where will the flow rise? Aka, where will the positive updraft be?

In my opinion, we should expect it to happen when the flow hits a vertical surface, and this is exactly what you see in your simulation. Think about climbing stairs - you don't step up on flat ground, but when there's a step in front of you.
If you want to be a little more precise, the presence of a vertical step forces the zonal velocity to zero, which leads to a horizontal convergence and consequently a local vertical flow. If nothing, if we find that the strength of these localized updrafts is too large it is a free surface issue. However, I am pretty sure that implementing the same setup in MITgcm would yield the same result.
I would expect this phenomenon to reduce in a non-hydrostatic flow with a correct pressure gradient where a flow convergence is mitigated by a nonhydrostatic pressure.

I suspect we do not see those vertical velocities in @hdrake's setup due to a couple of reasons (or the combination of those) that can mask this feature of a steppy bathymetry:
(1) a larger horizontal viscosity
(2) a ratio of vertical to horizontal surfaces which is closer to one
All this paired with a much lower forcing.
@francispoulin, note that increasing the vertical resolution alongside the horizontal one like you did does not remove these "spikes" because it does not change the ratio of horizontal to vertical steps.

It is possible to "reduce" or even "get rid" of these spikes if you play with one (or both) of the above.
As an example, here is the results of a simulation with the same setting but a steeper slope with a larger horizontal-to-vertical surface ratio. The vertical velocity "spikes" are completely removed because the updraft happens in all cells.

higher_vertical_to_horizontal.mp4

This does not mean that increasing horizontal viscosity or increasing the slope are a "solution" to this problem. I think we should acknowledge that a steppy bathymetry leads to a solution that necessarily has a "discontinuity" in the vertical velocity. But this example should be a further motivation to implement cut cells.

[francispoulin]

Thank you @simone-silvestri for your very thoughtful reply. I agree! At first I thought the spikes should go away under higher resolution, as did some others, but they don't and you have nicely explained why they have to be there. Essentially, there must be vertical motion and unlike a tilted bottom where it can happy gradually, here it can only happen in between cells. Thanks everyone for a very stimulating discussion!

I'll have one more comment below and then mabye close this discussion, or I can keep it open depending on how others feel.

[glwagner]

I have done some thinking and I maintain that those "spikes" are completely consistent with the discretized equations and a consequence of (1) the vertical grid (2) the bathymetric representation.

I think we all agree about this, don't we?

[glwagner]

I think we should acknowledge that a steppy bathymetry leads to a solution that necessarily has a "discontinuity" in the vertical velocity. But this example should be a further motivation to implement cut cells.

We also would like to understand thoroughly exactly the physics / geometry that produces this problem, as well as other aspects of a setup that affect this behavior (horizontal viscosity, presence of ambient turbulence, etc). I'd also like to understand how the pressure solver interacts with these near-boundary flows.

[francispoulin]

Again, thanks everyone for joining in this discussion!

Below i have plots of v, w, u, for the case of GridFittedBottom with a resolution of 400x400. The new part is the u profile, and following @simone-silvestri 's argument, this must be zero between the cells. I think it's tehre but unfortuantely, it's hard to see with the color scheme.

I am convinced the spikes are physical in the GridFittedBottom approach but given that PartlalCellBottom is noisy, it seems to me that there is probably a bug there. It should be similar in nature to the previous case. Do people agree with me?

Also, a silly question, if the MITgcm does shaved cells, can't we just translate that process into julia? A naive question but I am happy to give tihs a try with someone if people think this is doable and someone else wants to join in the fun.

[francispoulin]

I will close this discussion and suggest the following:

1. I create an issue saying we need to fix PartialCellBottom.
2. I create another issue saying we build on Implementing a CutCellBottom immersed boundary #3123 and work on a shaved/cut cell method.

That okay with everyone?

[glwagner]

No need to close the discussion! (I don't think we can close discussions anyways)

[glwagner]

Might make sense to continue investigating PartialCellBottom here. There's a lot of info already on this discussion that is helpful.

I create another issue saying we build on #3123 and work on a shaved/cut cell method.

Since #3123 already exists I don't think we need another issue to discuss it.

[navidcy]

I agree with @glwagner
We can continue discussing in this Discussion and, if we need to discuss about the #3123 then what makes most sense is to continue in #3123 rather than starting a new issue.

[francispoulin]

Good ideas and thanks for the insights!

In the next week I'll revisit the PartialCell method and see if I can find anythind odd in the code. I think this example is still a good test problem and I will look at it on a very coarse grid.

[navidcy]

Also, a silly question, if the MITgcm does shaved cells, can't we just translate that process into julia? A naive question but I am happy to give tihs a try with someone if people think this is doable and someone else wants to join in the fun.

#3123 is related; I don't think anyone is actively working on it atm

[glwagner]

I don't think there was ever progress on that PR, was there? We may want to close it so as not to mislead.

[navidcy]

oh I see, PR #3146
yes might be good idea to close

or also we could have a label "stale" and add to PRs/issues...?

[francispoulin]

All sound good to me.

If there is desire and time to push this forward soon I can create an issue or draft PR, but it's probably good to see this will happen soonish before we make that step.

[glwagner]

Should probably just close the stale ones, they are still in history and can be found but don't clutter our "to-do list". @francispoulin not sure what you are saying. There is already an issue. Why do you think we need a new one? As for a PR, you are definitely free to take up the implementation of CutCellBottom. Whether or not you use the existing PR, or a new PR, is entirely your choice in this case. I think there are some challenges though so you have to be ready for a long road!

[francispoulin]

I'll just look at the partial cell method for now and go from there.

[jm-c]

An ultra simple illustration (with MITgcm and partial cell and bottom friction), top is with usual viscosity and below is when I increase vertical viscosity just above bottom cell when the bottom cell thickness decreases.

[francispoulin]

Thank you/merci @jm-c!

Any chance I could find the parameters for your set up to try and reproduce the same thing in Oceananigans?

[glwagner]

Interesting that vertical viscosity solves the problem... is that because the vertical viscosity acts to decrease the along slope flow?

[francispoulin]

I am looking at PartialCellBottom and have a question.

I see that the grid stores the function to evaluate the bathymetry but what do I use to find the actual height used in the simulations? I know there is something called is_immersed but not sure if that's what we should use as I didn't find any examples.

[simone-silvestri]

This is the function that discriminates between immersed cells and fluid cells:

Oceananigans.jl/src/ImmersedBoundaries/partial_cell_bottom.jl

Line 97 in f7042b7

@inline function _immersed_cell(i, j, k, underlying_grid, ib::PartialCellBottom)

The dz is either full or "partial" based on these functions:

Oceananigans.jl/src/ImmersedBoundaries/partial_cell_bottom.jl

Line 107 in f7042b7

@inline function Δzᶜᶜᶜ(i, j, k, ibg::PCBIBG)

Oceananigans.jl/src/ImmersedBoundaries/partial_cell_bottom.jl

Line 127 in f7042b7

@inline function Δzᶜᶜᶠ(i, j, k, ibg::PCBIBG)

[francispoulin]

Thanks @simone-silvestri . I have seen these functions and good to know these are the ones I should focus on.

How do I import them directly into my scripts?

[navidcy]

using Oceananigans.ImmersedBoundaries: Δzᶜᶜᶠ ?

[francispoulin]

Thanks @navidcy !

After I also added immersed_cell that worked nicely. Now to enjoy these functions.

[francispoulin]

One thing that I have learned is that when I set

minimum_fractional_cell_height=0.0

in the partial cell method I don't recover the same thing as in the grid fitted bottom. Next step is to figure out why.

[hdrake]

Oh I see, we that makes sense.

[francispoulin]

Thanks @simone-silvestri . I was able to get the InterfaceImmersedCondition() working and it's similar to the standard case but slightly different.

I and do find that the PartialCellMethod is the same. This is good!

Now to figure out what goes wrong when the the parameter changes from 1. More soon I hope.

[francispoulin]

This is a result from GridFittedBottom or PartialCellBottom with minimum_fractional_cell_height=1 or even 0.99.

This is a result from PartialCellBottom with minimum_fractional_cell_height=0.95. Only a 5 percent difference yields a rather large difference. I find this surprising but I'm not sure if others do too?

[glwagner]

You do seem to be on to something. Is it possible that somewhere in the code, the metrics are used inconsistently?

Can you remind me about the physics of this case? Is there bottom drag?

[hdrake]

Can you plot them both with symmetric color maps? White does not correspond to zero so it's hard to read into color changes.

[francispoulin]

Sorry for the plots before.

Here are the same plots but with a balanced colorbar, the same for each.

Physically, this is initialized with a constant geostrophic velocity going along isobars into the page.

I agree something is odd here.

[francispoulin]

Thanks @simone-silvestri for the suggestion!

I tried that and the results don't change.

You can find my script for the GridFittedBottom here.

And my script for the PartialCellBottom here.

Any other suggestions?

[glwagner]

Can you let us know the vertical and horizontal CFL just so we are sure about what we are seeing?

[francispoulin]

I am using TimeStepWizard and setting cfl=0.1. I don't how how to find the horizontal and vertical CFL numbers but happy to do so if you let me know how.

[glwagner]

You can write cfl = AdvectiveCFL(dt) for a time-step dt and then call cfl(model).

I believe the wizard should work correctly because cell advection timescale uses the right metrics:

Oceananigans.jl/src/Advection/cell_advection_timescale.jl

Lines 19 to 27 in 1775f2b

	@inline function cell_advection_timescaleᶜᶜᶜ(i, j, k, grid, u, v, w)
	Δx = Δxᶠᶜᶜ(i, j, k, grid)
	Δy = Δyᶜᶠᶜ(i, j, k, grid)
	Δz = Δzᶜᶜᶠ(i, j, k, grid)
	inverse_timescale = @inbounds abs(u[i, j, k]) / Δx + abs(v[i, j, k]) / Δy + abs(w[i, j, k]) / Δz
	return 1 / inverse_timescale
	end

Another useful check is to make sure that minimum_zspacing(grid) indeed changes when you use PartialCellBottom.

[glwagner]

I guess another cool analysis would be to directly compute the difference between the two solutions. That'll help point out exactly what's changing...

[francispoulin]

This is a plot that compares the linear topography function I am using, along with what is used in both the GridFittedBottom and PartialCellBottom with ϵ = 0.1. The good news is that everything appears to be correct.

[The script can be found here.]

Maybe the next step is to look to see how the fluxes are computed in PartialCellBottom.

[glwagner]

I made a sketch of your situation:

On the underlying grid, the bottom cell spans -200 to -160. The bottom height, indicated by the black bar, is -190. Therefore, we should get a bottom tracer cell that spans from -190 to -160 --- 30 units high. The bottom Δzᶜᶜᶜ should be 30. Is this what you get?

The center of that cell is -160 - 15 = -175 --- 15 units below the next cell interface above at -160. Is this what you get?

Turning to Δzᶜᶜᶠ, the height of the first cell (ie k=1) is not important physically (this represents the height of a cell that spans the boundary). But this should span from -220 (the center of the halo cell below the boundary) to -175 --- Δzᶜᶜᶠ(k=1) = 45. The second cell, k=2, spans from -175 to -140 --- so Δzᶜᶜᶠ(k=2) = 35. It sounds like this is what you get?

Not sure what was expected, or if you think we are conceiving of this grid incorrectly...

[francispoulin]

Thanks @glwagner for your detaled response, with a nice picture!

I am happy to say that we both agree on the first part, what we get using the Δzᶜᶜᶜ functionand I understand this much better.

I didn't appreciate that we need to define the dual grid based on the reduced cell sizes, but I agree that makes sense.

I do find this to be true, Δzᶜᶜᶠ(k=2) = 35, however for the first one I get Δzᶜᶜᶠ(k=1) = 40. Do you think this value is used anywhere?

[glwagner]

I don't think it's used since that cell is not part of the domain. But I am not totally sure. We could try modifying the metric to be correct.

[francispoulin]

I also don't think this is going to make any difference but I admit that I am running out of ideas of things to test.

I suppose I could look at the difference in the immersed boundary solutions after one or two time steps might help point to where the differences occur.

Other ideas are very welcome.

[glwagner]

I would try reducing the resolution