Using tripolar grid for quarter-degree simulation

[suhanidalal]

Hi!

I am trying to increase the resolution of my global simulation to quarter-degree. I loosely followed this example for the one-degree version. After increasing the resolution, I run into the following error:

ArgumentError: halo=807 must be ≤ size=720 for coordinate y

I've tried increasing the number of grid points in y, but the error suggests a halo size of 807 in y, which I’m not explicitly setting. My grid setup is:

arch = GPU()
Nx = 1440 
Ny = 720 
Nz = 40
depth = 6000meters 
z = ExponentialCoordinate(Nz, -depth, 0; scale = 0.85*depth)
z = Oceananigans.Grids.MutableVerticalDiscretization(z)
underlying_grid = TripolarGrid(arch; size = (Nx, Ny, Nz), halo = (7,7,4), z) 

Is there something incorrect in my setup / is it possible to run a quarter-degree simulation on the tripolar grid?

Thank you!

[simone-silvestri]

H @suhanidalal. I think you are trying to run the ocean model without an immersed boundary. This will not work because the tripolar grid has singularities which will crash the simulation.

[glwagner]

@simone-silvestri this is an indirect response that does not address the issue, right? @suhanidalal I believe this error arises because you are using a split explicit free surface with a huge number of substeps. What @simone-silvestri is saying is that, there is some default behavior which will estimate the number of substeps needed, and it seems that you probably have a grid with tiny cells. This leads the scheme to estimate that you need O(100)s of free surface substeps.

If you add an immersed boundary, it will mask those tiny cells, and the same estimation procedure will produce a reasonable number of substeps --- probably.

Can you show us more of your code? The error is coming from the model constructor but you are not showing us the model.

We do actually need to address this issue @simone-silvestri --- what we have now is fragile and prone to breakage. We should figure out a more robust way to set up default configurations.

[simone-silvestri]

we could return a warning in the default free surface is someone is trying to run a tripolar grid that does not mask immersed cells or that has very tiny grid cells.

[suhanidalal]

Thanks for replying! I think I’ve figured it out, but I’d love to clarify a few things. I do have an immersed boundary (sorry I didn’t include it in my previous code excerpt):

arch = GPU()
Nx = 1440 
Ny = 720 
Nz = 40

depth = 6000meters 
z = ExponentialCoordinate(Nz, -depth, 0; scale = 0.85*depth)
z = Oceananigans.Grids.MutableVerticalDiscretization(z)

underlying_grid = TripolarGrid(arch; size = (Nx, Ny, Nz), halo = (7,7,4), z) 

bottom_height = regrid_bathymetry(underlying_grid;
                                  minimum_depth = 10,
                                  interpolation_passes = 5, 
                                  major_basins = 3) 

grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height);
                            active_cells_map=true)
tracer_advection = WENO(order=5)
bc = FieldBoundaryConditions()
println("done")
ocean = ocean_simulation(grid; tracers = (:T, :S, :e, :PG, :AC, :IMS, :IBS, :GSA, :GSB, :KC, :SS, :HC, :SFA),  boundary_conditions = (PG = bc, AC = bc, IMS = bc, IBS = bc, GSA = bc, GSB = bc, KC = bc, SS = bc, HC = bc, SFA = bc, ), tracer_advection)

But I had removed:

free_surface = SplitExplicitFreeSurface(grid; substeps=70)

for this simulation, since I didn't see the free surface included in the quarter-degree example. When I added it back in just now, it fixed my error.

For tripolar grids, is the free surface always necessary, and how should I determine the number of substeps? I'm also wondering if I should include other closure/advection schemes:

eddy_closure = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=2e3, κ_symmetric=2e3)
vertical_mixing = ClimaOcean.OceanSimulations.default_ocean_closure()
horizontal_viscosity = HorizontalScalarDiffusivity(ν=4000)
momentum_advection = WENOVectorInvariant(order=5)

Since I increased the resolution, I assumed I could omit these, as smaller scales are better resolved, but I wanted to confirm.

Thank you again!

[glwagner]

The number of substeps must be determined by the gravity wave CFL, eg so that you resolve the gravity wave speed sqrt(g * H) over a substep.

There is an option to use the cfl to determine the number of substeps. I can never remember the syntax, but the docstring may have documentation, cc @simone-silvestri @navidcy

[glwagner]

@suhanidalal there is a quarter degree example here: https://numericalearth.github.io/NumericalEarthDocumentation/dev/literated/near_global_ocean_simulation/

Note, this discussion is on ClimaOcean but development has moved to NumericalEarth (please see the README).

[glwagner]

For tripolar grids, is the free surface always necessary, and how should I determine the number of substeps?

Unfortunately it seems so. There have been efforts in the past to produce a default free surface but apparently we have not yet succeeded and manually configuration is still required. @simone-silvestri may have more information.

[simone-silvestri]

It looks like we must have a bug in our defaults. The free surface is indeed always necessary to have for the ocean_simulation of NumericalEarth but it should not be necessary to specify it explicitly since the defaults that we have set up in place should kick in when the free surface is not explicitly declared.

When it comes to closures and advection schemes, for the quarter degree I would use the default schemes (so do not specify them explicitly). This is the same we do in the quarter degree simulation, just

ocean = ocean_simulation(grid)

The defaults are WENO order 7 for tracers advection, WENOVectorInvariant for momentum advection, and CATKEVerticalDiffusivity as a vertical mixing scheme (you can find the reference for these schemes and parameterizations in the Oceananigans codebase)

[suhanidalal]

Got it, thank you both! This was really helpful.