2d simulation with variable grid spacings

[parfenyev]

Hello everyone,

I am trying to perform 2d simulation with variable grid spacings in x- and y-directions. The following is a minimal example that gives the error:

using Oceananigans, Statistics, Printf

Nx, Ny = 512, 512;
Lx, Ly = 2π, 2π;

chebychev_spaced_x_faces(i) = Lx/2 - Lx/2 * cos(π * (i - 1) / Nx);
chebychev_spaced_y_faces(j) = Ly/2 - Ly/2 * cos(π * (j - 1) / Ny);

grid = RectilinearGrid(size = (Nx, Ny), 
                              topology=(Periodic, Periodic, Flat),
                              x = chebychev_spaced_x_faces,
                              y = chebychev_spaced_y_faces)

model = NonhydrostaticModel(timestepper = :RungeKutta3,
                              advection = UpwindBiasedFifthOrder(),
                                   grid = grid,
                               buoyancy = nothing,
                                tracers = nothing,
                                closure = ScalarDiffusivity(ν=1e-5))

MethodError: no method matching PressureSolver(::CPU, ::Oceananigans.Grids.ZRegRectilinearGrid{Float64, Periodic, Periodic, Flat, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, Float64, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, CPU})
Closest candidates are:
  PressureSolver(::Any, ::RectilinearGrid{<:Any, <:Any, <:Any, <:Any, <:Number, <:Number, <:Number}) at C:\Users\parfe\.julia\packages\Oceananigans\jmNfq\src\Models\NonhydrostaticModels\NonhydrostaticModels.jl:28
  PressureSolver(::Any, ::RectilinearGrid{<:Any, <:Any, <:Any, <:Any, <:Number, <:Number}) at C:\Users\parfe\.julia\packages\Oceananigans\jmNfq\src\Models\NonhydrostaticModels\NonhydrostaticModels.jl:29
  PressureSolver(::Any, ::Oceananigans.ImmersedBoundaries.ImmersedBoundaryGrid) at C:\Users\parfe\.julia\packages\Oceananigans\jmNfq\src\Models\NonhydrostaticModels\NonhydrostaticModels.jl:32

Stacktrace:
 [1] NonhydrostaticModel(; grid::Oceananigans.Grids.ZRegRectilinearGrid{Float64, Periodic, Periodic, Flat, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, Float64, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, CPU}, clock::Clock{Float64}, advection::UpwindBiasedFifthOrder, buoyancy::Nothing, coriolis::Nothing, stokes_drift::Nothing, forcing::NamedTuple{(), Tuple{}}, closure::ScalarDiffusivity{Oceananigans.TurbulenceClosures.ExplicitTimeDiscretization, Oceananigans.TurbulenceClosures.ThreeDimensionalFormulation, Float64, Float64}, boundary_conditions::NamedTuple{(), Tuple{}}, tracers::Nothing, timestepper::Symbol, background_fields::NamedTuple{(), Tuple{}}, particles::Nothing, velocities::Nothing, pressures::Nothing, diffusivity_fields::Nothing, pressure_solver::Nothing, immersed_boundary::Nothing, auxiliary_fields::NamedTuple{(), Tuple{}})
   @ Oceananigans.Models.NonhydrostaticModels C:\Users\parfe\.julia\packages\Oceananigans\jmNfq\src\Models\NonhydrostaticModels\nonhydrostatic_model.jl:185
 [2] top-level scope
   @ In[45]:14
 [3] eval
   @ .\boot.jl:373 [inlined]
 [4] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base .\loading.jl:1196

How can this be fixed?

[glwagner]

Variable horizontal grid spacing isn't supported by NonhydrostaticModel, because the only pressure solvers that are implemented are fully FFT based (for regular in all directions) and a combined FFT + tridiagonal solve (for regular in horizontal directions, stretched in the vertical).

[glwagner]

True, but there's no closure here! I've only compared cost per time-step, not time-to-solution (effectively). I think if we have explicit diffusion we would indeed require a very small time-step, thus exacerbating the problem of a slow time-step further. One recourse in the finite volume context is an iterative solver for implicit diffusion, which would again increase the computational expense, likely by another factor of 3 in this two-dimensional problem (but much more in three-dimensions). Numerically this is not all that challenging, but implementing it would be a substantial extension of existing functionality.

[glwagner]

Exasim might be able to handle this problem --- it treates the entire tendency implicitly using sophisticated general iterative solvers (that work for nonlinear and linear terms!), so is probably well suited to problems with extreme grid stretching like this one: https://github.com/exapde/Exasim

[parfenyev]

Thank you very much! I have found that people use pseudospectral methods based on Chebyshev polynomials to solve such problems, see e.g. this paper. Can you comment on whether such a method will give a significant increase in the speed of calculations in comparison with the hydrostatic model from Oceananigans.jl?

[glwagner]

I don't know. I think double Chebyshev with implicit diffusion could be a "good choice" for some applications --- like two-dimensional, doubly-bounded Navier-Stokes simulations that resolve viscous boundary layers but have high Reynolds number in the interior. But I'm not an expert.

[francispoulin]

I also can't say much but it depends very much on what you want to do. If you want to resolve boundary layers then perhaps a cheb grid is the way to go, but it will always be more expensive.

[glwagner]

using Oceananigans, Statistics, Printf

Nx, Ny = 512, 512;
Lx, Ly = 2π, 2π;

chebychev_spaced_x_faces(i) = Lx/2 - Lx/2 * cos(π * (i - 1) / Nx);
chebychev_spaced_y_faces(j) = Ly/2 - Ly/2 * cos(π * (j - 1) / Ny);

grid = RectilinearGrid(size = (Nx, Ny, 1), 
                       topology=(Bounded, Bounded, Bounded),
                       # topology=(Periodic, Periodic, Bounded),
                       x = chebychev_spaced_x_faces,
                       y = chebychev_spaced_y_faces,
                       z = (0, 1))

model = HydrostaticFreeSurfaceModel(; grid,
                                    momentum_advection = UpwindBiasedFifthOrder(),
                                    buoyancy = nothing,
                                    tracers = nothing,
                                    free_surface = ImplicitFreeSurface(gravitational_acceleration=10), # c = sqrt(10)
                                    closure = ScalarDiffusivity(ν=1e-3))

ϵ(x, y, z) = 2rand() - 1
set!(model, u=ϵ, v=ϵ)

simulation = Simulation(model, Δt=0.01, stop_time=10)
run!(simulation)

might work

[francispoulin]

I tried changing the model as follows and it didn't complain, but I also didn't go any further.

 model = HydrostaticFreeSurfaceModel(      
                                                     grid = grid,
                                              buoyancy = nothing,
                                               tracers = nothing,
                                               closure = ScalarDiffusivity(ν=1e-5))
HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 512×512×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 1×1×0 halo
├── timestepper: QuasiAdamsBashforth2TimeStepper
├── tracers: ()
├── closure: ScalarDiffusivity{ExplicitTimeDiscretization}(ν=1.0e-5, κ=NamedTuple())
├── buoyancy: Nothing
├── free surface: ExplicitFreeSurface with gravitational acceleration 9.80665 m s⁻²
└── coriolis: Nothing

[glwagner]

Do we really want horizontally-periodic with variable spacing? Or do we want topology = (Bounded, Bounded, Bounded)?

[francispoulin]

I agree that if it's doubly periodic then using a cheb grid might not be very efficient

I would think that topology = (Bounded, Bounded, Flat) would make sense and then you could get boundary layers forming, which is maybe one of the motivations for looking at this?

[glwagner]

We can't do Flat in the vertical with a free surface (just because of the way the algorithm is written). It might not be hard to do that refactor though.

[francispoulin]

Sorry, my mistakes. Thanks for clarifying.

[francispoulin]

Using topology = (Bounded, Bounded, Flat), I was able to create a ShallowWaterModel (see below). This would be maybe a simpler way to study turbulence as there is no inversion to worry about so it should not be too slow. There is only the CFL conditon that could slow things down.

However we have no closer yet. Maybe we could talk about creating closure for this model @glwagner ?

model = ShallowWaterModel(       grid = grid,
                                        gravitational_acceleration = 1.0)
ShallowWaterModel{typename(CPU), Float64}(time = 0 seconds, iteration = 0) 
├── grid: 512×512×1 RectilinearGrid{Float64, Bounded, Bounded, Flat} on CPU with 3×3×0 halo
├── tracers: ()
└── coriolis: Nothing

[glwagner]

We can talk! I may not be able to devote much time to it right now. I think the shallow water model needs a few more things to make this viable: fixing advection operators, and implicit treatment of the h-equation (so we don't have to resolve wave speed).

That said I actually think the shallow water model is a great test bed for idealized stuff, because its simple (and we should keep it that way). We could develop rigid lid, linear-h and implicit-h capabilities and then it might be even slightly more optimized than NonhydrostaticModel for a range of 2D problems.

[francispoulin]

All great stuff and happy to talk when you have time, but know you are very busy on many other fronts right now