search

CliMA / Oceananigans.jl

🌊 Julia software for fast, friendly, flexible, ocean-flavored fluid dynamics on CPUs and GPUs
https://clima.github.io/OceananigansDocumentation/stable
MIT License
1k stars 196 forks source link

Laplacian on `ImmersedBoundaryGrid` is incorrect if using `∂x`, `∂y`, and `∂z` composition #2915

Open xkykai opened 1 year ago

xkykai commented 1 year ago

Computing the Laplacian on ImmersedBoundaryGrid using different methods give different results. Here's a MWE:

julia> using Oceananigans

julia> using Oceananigans.Operators: ∇²ᶜᶜᶜ, ∂²xᶜᶜᶜ, ∂²yᶜᶜᶜ, ∂²zᶜᶜᶜ, ∂xᶠᶜᶜ, ∂yᶜᶠᶜ, ∂zᶜᶜᶠ, ∂xᶜᶜᶜ, ∂yᶜᶜᶜ, ∂zᶜᶜᶜ, δxᶜᵃᵃ, δyᵃᶜᵃ, δzᵃᵃᶜ, Δxᶜᶜᶜ, Δyᶜᶜᶜ, Δzᶜᶜᶜ

julia> grid = RectilinearGrid(size = (5, 1, 5),
                              halo = (4, 4, 4),
                              x = (0, 1),
                              y = (0, 1),
                              z = (-1, 0),
                              topology = (Bounded, Periodic, Bounded))

julia> slope(x, y) = -0.8

julia> grid = ImmersedBoundaryGrid(grid, GridFittedBottom(slope))

julia> field = CenterField(grid)

julia> rand_field(x, y, z) = rand()

julia> set!(field, rand_field)

julia> ∇² = [∇²ᶜᶜᶜ(i, 1, k, grid, field) for i in 0:6, k in 0:6]
7×7 Matrix{Float64}:
 0.0  0.0    0.0        0.0        0.0        0.0      0.0
 0.0  0.0   33.3801   -51.1901    44.2001     7.36811  0.0
 0.0  0.0   -3.44097  -29.2144    -1.68447  -29.1814   0.0
 0.0  0.0  -42.4603    61.5155   -12.8827    -7.86909  0.0
 0.0  0.0   16.4908     6.82589  -41.1333    52.9584   0.0
 0.0  0.0   -5.28556   23.4819   -21.1286   -20.672    0.0
 0.0  0.0    0.0        0.0        0.0        0.0      0.0

julia> ∂² = [∂²xᶜᶜᶜ(i, 1, k, grid, field) +
             ∂²yᶜᶜᶜ(i, 1, k, grid, field) +
             ∂²zᶜᶜᶜ(i, 1, k, grid, field) for i in 0:6, k in 0:6]
7×7 Matrix{Float64}:
 0.0  0.0    0.0        0.0        0.0        0.0      0.0
 0.0  0.0   28.7172   -76.0349    43.2687    -7.21832  0.0
 0.0  0.0    5.28475  -29.2144    -1.68447  -50.788    0.0
 0.0  0.0  -59.1201    61.5155   -12.8827   -23.0124   0.0
 0.0  0.0   26.3115     6.82589  -41.1333    51.5877   0.0
 0.0  0.0   -2.94252   19.9348   -42.7884   -62.6937   0.0
 0.0  0.0    0.0        0.0        0.0        0.0      0.0

julia> ∂∂ = [∂xᶜᶜᶜ(i, 1, k, grid, ∂xᶠᶜᶜ, field) +
             ∂yᶜᶜᶜ(i, 1, k, grid, ∂yᶜᶠᶜ, field) +
             ∂zᶜᶜᶜ(i, 1, k, grid, ∂zᶜᶜᶠ, field) for i in 0:6, k in 0:6]
7×7 Matrix{Float64}:
 0.0  0.0    0.0        0.0        0.0        0.0      0.0
 0.0  0.0   28.7172   -76.0349    43.2687    -7.21832  0.0
 0.0  0.0    5.28475  -29.2144    -1.68447  -50.788    0.0
 0.0  0.0  -59.1201    61.5155   -12.8827   -23.0124   0.0
 0.0  0.0   26.3115     6.82589  -41.1333    51.5877   0.0
 0.0  0.0   -2.94252   19.9348   -42.7884   -62.6937   0.0
 0.0  0.0    0.0        0.0        0.0        0.0      0.0

julia> δ∂ = [δxᶜᵃᵃ(i, 1, k, grid, ∂xᶠᶜᶜ, field) / Δxᶜᶜᶜ(i, 1, k, grid) +
             δyᵃᶜᵃ(i, 1, k, grid, ∂yᶜᶠᶜ, field) / Δyᶜᶜᶜ(i, 1, k, grid) +
             δzᵃᵃᶜ(i, 1, k, grid, ∂zᶜᶜᶠ, field) / Δzᶜᶜᶜ(i, 1, k, grid) for i in 0:6, k in 0:6]
7×7 Matrix{Float64}:
 0.0  0.0    0.0        0.0        0.0        0.0      0.0
 0.0  0.0   33.3801   -51.1901    44.2001     7.36811  0.0
 0.0  0.0   -3.44097  -29.2144    -1.68447  -29.1814   0.0
 0.0  0.0  -42.4603    61.5155   -12.8827    -7.86909  0.0
 0.0  0.0   16.4908     6.82589  -41.1333    52.9584   0.0
 0.0  0.0   -5.28556   23.4819   -21.1286   -20.672    0.0
 0.0  0.0    0.0        0.0        0.0        0.0      0.0

julia> ∇² .≈ ∂²
7×7 BitMatrix:
 1  1  1  1  1  1  1
 1  1  0  0  0  0  1
 1  1  0  1  1  0  1
 1  1  0  1  1  0  1
 1  1  0  1  1  0  1
 1  1  0  0  0  0  1
 1  1  1  1  1  1  1

julia> ∇² .≈ δ∂
7×7 BitMatrix:
 1  1  1  1  1  1  1
 1  1  1  1  1  1  1
 1  1  1  1  1  1  1
 1  1  1  1  1  1  1
 1  1  1  1  1  1  1
 1  1  1  1  1  1  1
 1  1  1  1  1  1  1

julia> ∂² .≈ ∂∂
7×7 BitMatrix:
 1  1  1  1  1  1  1
 1  1  1  1  1  1  1
 1  1  1  1  1  1  1
 1  1  1  1  1  1  1
 1  1  1  1  1  1  1
 1  1  1  1  1  1  1
 1  1  1  1  1  1  1

From above, ∇² is consistent with δ∂, but they do not agree with ∂∂. The deviations appear at the regular and immersed boundaries. All 3 methods of calculating Laplacian yield the same result in a non-immersed RectilinearGrid.

It appears that ∇² and δ∂ are the correct solutions, so the problem is likely to be in the composing of 2 derivative operations on an ImmersedBoundaryGrid. I tried to fix it myself, but I couldn't as I am unfamiliar with the metaprogramming aspects of the code.

navidcy commented 1 year ago

cc @simone-silvestri

simone-silvestri commented 1 year ago

Hmmm, ∇² and δ∂ are executing the same code (under the hood the ∇²ᶜᶜᶜ is calculated as δ∂ for a constant spacing rectilinear grid). The difference between the two calculations yielding different results is how many times you account for an immersed boundary: ∂ operators are boundary-aware (i.e. yield a zero if calculated across an immersed boundary), while δ operators are not. Despite this, they should yield the same result since the second pass of the derivative should not cross immersed boundaries. I ll take a look

navidcy commented 1 year ago

This

using Oceananigans
using Oceananigans.Operators

grid = RectilinearGrid(size = (5, 1, 5),
                       halo = (4, 4, 4),
                       x = (0, 1),
                       y = (0, 1),
                       z = (-1, 0),
                       topology = (Bounded, Periodic, Bounded))

slope(x, y) = -0.8

grid = ImmersedBoundaryGrid(grid, GridFittedBottom(slope))

field = CenterField(grid)

rand_field(x, y, z) = rand()

set!(field, rand_field)

∂x = [∂xᶜᶜᶜ(i, 1, k, grid, field) for i in 0:6, k in 0:6]

∂δx = [δxᶜᵃᵃ(i, 1, k, grid, field) / Δxᶜᶜᶜ(i, 1, k, grid) for i in 0:6, k in 0:6]

∂x .≈ ∂δx

would give 

7×7 BitMatrix:
 1  0  0  0  0  0  1
 1  0  1  1  1  1  1
 1  0  1  1  1  1  1
 1  0  1  1  1  1  1
 1  0  1  1  1  1  1
 1  0  1  1  1  1  1
 1  1  1  1  1  1  1
navidcy commented 1 year ago

Though I'm not sure if we'd expect this to work since δxᶜᵃᵃ are not immersed-boundary aware, right?

simone-silvestri commented 1 year ago

yeah, the difference is expected. Now we have to understand if the immersed version is correct

xkykai commented 1 year ago

∇²ᶜᶜᶜ gives the right solution for the pressure solver (correct as in no flow divergence at the boundary) that we are developing for the ImmersedPoissonSolver here: https://github.com/xkykai/Oceananigans.jl/blob/df579ed61935378273f2f0e3df0706554be22c30/validation/immersed_boundaries/immersed_pressure_solver.jl#L114

simone-silvestri commented 1 year ago

Remember to fill the halos before computing the laplacian. Anyways, it looks like ∂² gives the "correct" answer when I look at analytical solutions. Still, these operators should not give different results

julia> set!(field, (x, y, z) -> z^2)
5×1×5 Field{Center, Center, Center} on ImmersedBoundaryGrid on CPU
├── grid: 5×1×5 ImmersedBoundaryGrid{Float64, Bounded, Periodic, Bounded} on CPU with 4×4×4 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: ZeroFlux, east: ZeroFlux, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 13×9×13 OffsetArray(::Array{Float64, 3}, -3:9, -3:5, -3:9) with eltype Float64 with indices -3:9×-3:5×-3:9
    └── max=0.49, min=0.01, mean=0.21

julia> fill_halo_regions!(field)

julia> ∇² = [∇²ᶜᶜᶜ(i, 1, k, grid, field) for i in 1:5, k in 1:5]
5×5 Matrix{Float64}:
 0.0  -6.0  2.0  2.0  2.0
 0.0  -6.0  2.0  2.0  2.0
 0.0  -6.0  2.0  2.0  2.0
 0.0  -6.0  2.0  2.0  2.0
 0.0  -6.0  2.0  2.0  2.0

julia> ∂² = [∂²xᶜᶜᶜ(i, 1, k, grid, field) +
             ∂²yᶜᶜᶜ(i, 1, k, grid, field) +
             ∂²zᶜᶜᶜ(i, 1, k, grid, field) for i in 1:5, k in 1:5]
5×5 Matrix{Float64}:
 0.0  2.0  2.0  2.0  2.0
 0.0  2.0  2.0  2.0  2.0
 0.0  2.0  2.0  2.0  2.0
 0.0  2.0  2.0  2.0  2.0
 0.0  2.0  2.0  2.0  2.0

Edit: actually ∇² seems to be correct because of the immersed boundary

simone-silvestri commented 1 year ago

Ok, I see why this is happening. conditional_∂z_c "throws away" the immersed boundary when it passes through. Since a second derivative is just a first derivative of a first derivative, the immersed condition is tested only on the "outer" derivative, which doesn't see the immersed boundary. The inner derivative is then just called on the underlying grid and does not satisfy immersed boundary conditions. This does not happen in the ∇² operator since the "inner" derivative correctly calls the conditional_∂z_c, thus being aware of the immersed boundaries. Good catch, I ll solve this issue

glwagner commented 1 year ago

there it is

glwagner commented 1 year ago

For documentation:

https://github.com/CliMA/Oceananigans.jl/blob/f70d0949f7688df89a7f8a9e7d44177cde27091a/src/Operators/laplacian_operators.jl#L36-L40