Bug in diffusion operators?

thabbott commented 4 years ago

I think the diffusive fluxes in compressible_operators.jl:

@inline diffusive_flux_x(i, j, k, grid, κᶠᶜᶜ, ρ, C) = κᶠᶜᶜ * Axᵃᵃᶠ(i, j, k, grid) * δxᶠᵃᵃ(i, j, k, grid, C_over_ρ, C, ρ)
@inline diffusive_flux_y(i, j, k, grid, κᶜᶠᶜ, ρ, C) = κᶜᶠᶜ * Ayᵃᵃᶠ(i, j, k, grid) * δyᵃᶠᵃ(i, j, k, grid, C_over_ρ, C, ρ)
@inline diffusive_flux_z(i, j, k, grid, κᶜᶜᶠ, ρ, C) = κᶜᶜᶠ * Azᵃᵃᵃ(i, j, k, grid) * δzᵃᵃᶠ(i, j, k, grid, C_over_ρ, C, ρ)

@inline function ∂ⱼκᵢⱼ∂ᵢc(i, j, k, grid, κˣ, κʸ, κᶻ, ρᵈ, C)
    return 1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, diffusive_flux_x, κˣ, ρᵈ, C) +
                                    δyᵃᶜᵃ(i, j, k, grid, diffusive_flux_y, κʸ, ρᵈ, C) +
                                    δzᵃᵃᶜ(i, j, k, grid, diffusive_flux_z, κᶻ, ρᵈ, C))
end

are missing a division by the grid spacing. (Without it the units don't work out.) This might explain when we've been finding that the verification experiments require a smaller diffusivity than is typically reported in publications.

ali-ramadhan commented 4 years ago

The diffusive fluxes have an area factor in them, e.g. Axᵃᵃᶠ, and then a division by a volume factor Vᵃᵃᶜ multiplying the flux divergence which should give the correct 1/distance factor. Does that make sense?

Might be a horrible way to write things since we're just using regular Cartesian grids right now, but I guess I wrote it that way to generalize to more general Cartesian grids.

Also see: https://climate-machine.github.io/Oceananigans.jl/stable/numerical_implementation/spatial_operators/#Discretization-of-isotropic-diffusion-operators-1

thabbott commented 4 years ago

Yes, but even accounting for those geometric factors I think they're missing a factor of 1/length. (In part just from checking units, and also because the A/V factor appears because a given flux on a face produces a smaller volume-average tendency if A/V is smaller---but the diffusive flux itself should be smaller for a given difference between adjacent cells if the centers of those cells are farther apart.) FWIW a 75 m^2/s viscosity gives reasonable results with this (possible) bug fixed.

On Thu, Feb 6, 2020, 8:47 PM Ali Ramadhan notifications@github.com wrote:

The diffusive fluxes have an area factor in them, e.g. Axᵃᵃᶠ, and then a division by a volume factor Vᵃᵃᶜ multiplying the flux divergence which should give the correct 1/distance factor. Does that make sense?

Might be a horrible way to write things since we're just using regular Cartesian grids right now, but I guess I wrote it that way to generalize to more general Cartesian grids.

Also see: https://climate-machine.github.io/Oceananigans.jl/stable/numerical_implementation/spatial_operators/#Discretization-of-isotropic-diffusion-operators-1

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/thabbott/JULES.jl/issues/56?email_source=notifications&email_token=ACZDSTQMS7CRV6Z2FNHYJ7LRBS4UTA5CNFSM4KRC4QJKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOELBN5QI#issuecomment-583196353, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACZDSTQDEUBG2TUOSQP67W3RBS4UTANCNFSM4KRC4QJA .

ali-ramadhan commented 4 years ago

Ah sorry you're right. In the case of Laplacian diffusion you should end up with a 1/length² factor so yeah δxᶠᵃᵃ (just a difference) should be ∂xᶠᵃᵃ (includes a division by a grid spacing like a derivative).

Oceananigans does it right so no idea how I messed up the diffusion operators: https://github.com/climate-machine/Oceananigans.jl/blob/3e803ac0d2e9252654384137aafaef635621285b/src/TurbulenceClosures/diffusion_operators.jl#L5-L7

ali-ramadhan / Atmosfoolery.jl

Bug in diffusion operators? #56