GitHub

A travelling top hat: incremental remapping against flux-form advection

This example advects a top hat of sea ice concentration through a periodic channel at a uniform speed and watches what happens to the thickness recovered from the transported quantities. It demonstrates:

  • why the sea ice model transports concentration ℵ and content 𝓋 = ℵh rather than h itself,
  • how the recovered thickness h = 𝓋/ℵ behaves at an ice edge,
  • the difference between flux-form advection and IncrementalRemapping.

The exact answer is a rigid translation: the hat keeps its shape, and the thickness stays at 1.5 m everywhere inside it. Nothing in the problem creates or destroys ice.

Install dependencies

using Pkg
pkg"add Oceananigans, ClimaSeaIce, CairoMakie"
using ClimaSeaIce
using Oceananigans
using Oceananigans.BoundaryConditions
using Printf
using CairoMakie

The domain

The problem is one-dimensional. We give the grid a few cells in y only because the model is two-dimensional; nothing varies across the channel. Incremental remapping reads one cell beyond its stencil in both directions, so the halo has to be at least two.

Nx = 200

grid = RectilinearGrid(size = (Nx, 8, 1),
                       x = (0, 1),
                       y = (0, 1),
                       z = (-1, 0),
                       halo = (5, 5, 5),
                       topology = (Periodic, Periodic, Bounded))
200×8×1 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 5×5×5 halo
├── Periodic x ∈ [0.0, 1.0)  regularly spaced with Δx=0.005
├── Periodic y ∈ [0.0, 1.0)  regularly spaced with Δy=0.125
└── Bounded  z ∈ [-1.0, 0.0] regularly spaced with Δz=1.0

The flow

A uniform current carries the ice to the right. One revolution returns the hat to where it started.

U = 1

u = Field{Face, Center, Center}(grid)
v = Field{Center, Face, Center}(grid)

set!(u, U)
set!(v, 0)

fill_halo_regions!(u)
fill_halo_regions!(v)

The initial condition

A top hat of fully packed ice, 1.5 m thick, surrounded by open water. The concentration is discontinuous; the thickness is uniform wherever there is ice.

iced(x) = 0.3 < x < 0.6

ℵᵢ(x, y) = iced(x) ? 1.0 : 0.0
hᵢ(x, y) = iced(x) ? 1.5 : 0.0
hᵢ (generic function with 1 method)

Three configurations

We compare flux-form WENO with both of the model's timesteppers against incremental remapping. Incremental remapping transports the material swept through each face over the whole step, so it carries its own time integration and is a Forward-Euler scheme by construction.

configurations = [(name = "WENO5, RK3",            advection = WENO(order=5),          timestepper = :SplitRungeKutta3),
                  (name = "WENO5, Forward Euler",  advection = WENO(order=5),          timestepper = :ForwardEuler),
                  (name = "Incremental remapping", advection = IncrementalRemapping(), timestepper = :ForwardEuler),
                  (name = "UpwindBiased 1",        advection = UpwindBiased(order=1),  timestepper = :ForwardEuler)]
4-element Vector{NamedTuple{(:name, :advection, :timestepper)}}:
 (name = "WENO5, RK3", advection = WENO{3, Float64, Nothing}(order=5)
├── buffer_scheme: WENO{2, Float64, Nothing}(order=3)
│   └── buffer_scheme: Centered(order=2)
└── advecting_velocity_scheme: Centered(order=4), timestepper = :SplitRungeKutta3)
 (name = "WENO5, Forward Euler", advection = WENO{3, Float64, Nothing}(order=5)
├── buffer_scheme: WENO{2, Float64, Nothing}(order=3)
│   └── buffer_scheme: Centered(order=2)
└── advecting_velocity_scheme: Centered(order=4), timestepper = :ForwardEuler)
 (name = "Incremental remapping", advection = IncrementalRemapping(quadrature_nodes=3), timestepper = :ForwardEuler)
 (name = "UpwindBiased 1", advection = UpwindBiased(order=1)
├── buffer_scheme: Nothing
└── advecting_velocity_scheme: Centered(order=2), timestepper = :ForwardEuler)

Running

The Courant number is 0.25 and we integrate for one revolution, saving a snapshot every eight steps.

courant = 0.25
Δt = courant / (U * Nx)
steps = round(Int, 1 / (U * Δt))
save_interval = 8

function run(configuration)
    model = SeaIceModel(grid; velocities = (; u, v),
                        advection = configuration.advection,
                        timestepper = configuration.timestepper,
                        ice_thermodynamics = nothing,
                        dynamics = nothing)

    set!(model, ℵ = ℵᵢ, h = hᵢ)

    ℵ = [Array(interior(model.ice_concentration))[:, 1, 1]]
    h = [Array(interior(model.ice_thickness))[:, 1, 1]]

    for n in 1:steps
        time_step!(model, Δt)

        if n % save_interval == 0
            push!(ℵ, Array(interior(model.ice_concentration))[:, 1, 1])
            push!(h, Array(interior(model.ice_thickness))[:, 1, 1])
        end
    end

    @printf("%-22s peak h = %12.4f m\n", configuration.name, maximum(maximum, h))

    return (; ℵ, h)
end

solutions = map(run, configurations)
4-element Vector{@NamedTuple{ℵ::Vector{Vector{Float64}}, h::Vector{Vector{Float64}}}}:
 (ℵ = [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]  …  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]], h = [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]  …  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]])
 (ℵ = [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]  …  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]], h = [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]  …  [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]])
 (ℵ = [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  1.0727429585864132e-83, 8.043058733543874e-87, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]  …  [1.7676784845398517e-25, 7.869275634756825e-25, 3.488587610299336e-24, 1.5405069528040195e-23, 6.777804387411167e-23, 2.971897098652847e-22, 1.298976309343864e-21, 5.660954469103872e-21, 2.4603214162489464e-20, 1.0665835020906227e-19  …  4.3199274702622644e-32, 2.036365256800712e-31, 9.532525207211897e-31, 4.432509312941503e-30, 2.0478563077870868e-29, 9.403253116250031e-29, 4.292475313226813e-28, 1.948550751106111e-27, 8.798598246847189e-27, 3.953088171908666e-26], [9.35800169007482e-27, 4.1928000704723263e-26, 1.8700352887296025e-25, 8.304967342365889e-25, 3.673532418974382e-24, 1.6188247039453983e-23, 7.108765416689509e-23, 3.111509040751486e-22, 1.357781269539001e-21, 5.908328035263796e-21  …  2.1032185251450326e-33, 1.0012472719790196e-32, 4.731975884002064e-32, 2.220743031617563e-31, 1.0351942168078976e-30, 4.7943387544247626e-30, 2.206666367917121e-29, 1.0096367481024292e-28, 4.593405630985811e-28, 2.0785810929873477e-27], [4.9062284906255e-28, 2.2134343463717466e-27, 9.936944454596525e-27, 4.4404274092155936e-26, 1.9756010677638556e-25, 8.753687374645708e-25, 3.8637772800731926e-24, 1.699303703908604e-23, 7.448562318787276e-23, 3.2547428167948133e-22  …  1.010096364396926e-34, 4.857641197829222e-34, 2.3185395091739767e-33, 1.098586316682483e-32, 5.168822891702724e-32, 2.415443408408987e-31, 1.1214046556673754e-30, 5.173719551349487e-30, 2.3726577209248552e-29, 1.0818749143387882e-28], [2.545883561644052e-29, 1.157058014705773e-28, 5.231009914232335e-28, 2.3531364318248456e-27, 1.0535536840488564e-26, 4.6960200820250745e-26, 2.0843991679977642e-25, 9.215546164305772e-25, 4.059374222528806e-24, 1.7819681750392048e-23  …  4.784187516103846e-36, 2.3248070122605913e-35, 1.1209471897959102e-34, 5.364149367850114e-34, 2.548209120912466e-33, 1.2019648937807729e-32, 5.63090071275961e-32, 2.620612452740112e-31, 1.2119324029963496e-30, 5.570811883585264e-30], [1.3068222814937226e-30, 5.98576488300525e-30, 2.726394385515042e-29, 1.2352038815914672e-28, 5.567817192505789e-28, 2.4977147958348863e-27, 1.1153898814178878e-26, 4.9596330202248244e-26, 2.1964553059650397e-25, 9.690664578533365e-25  …  2.2343032921812057e-37, 1.0973134810509705e-36, 5.3461656945133284e-36, 2.584455334122143e-35, 1.2399537008112449e-34, 5.905405152227359e-34, 2.7925610072616138e-33, 1.311499664558057e-32, 6.118609140642089e-32, 2.836384806437985e-31], [6.632318000696434e-32, 3.062884886387257e-31, 1.4061159064959025e-30, 6.418718872353408e-30, 2.914238872531703e-29, 1.3163303062479206e-28, 5.9167200575304995e-28, 2.6471988021125423e-27, 1.1792161219260338e-26, 5.231326528292835e-26  …  1.0287482514262521e-38, 5.10723524567602e-38, 2.5147756761921008e-37, 1.2283869145486578e-36, 5.9536402146622203e-36, 2.863739034866953e-35, 1.3673595749075773e-34, 6.482289650997123e-34, 3.0519073255121643e-33, 1.427299046193494e-32], [3.326539564085678e-33, 1.549463725669305e-32, 7.172369812749578e-32, 3.3002281246936464e-31, 1.5098522219324327e-30, 6.869807545502494e-30, 3.1094647274807834e-29, 1.400455345998803e-28, 6.277791745050902e-28, 2.8016200757464437e-27  …  4.669542723323841e-40, 2.3437068656381162e-39, 1.1665186844835644e-38, 5.758639405107455e-38, 2.8201425078617876e-37, 1.370349566778791e-36, 6.608303130326312e-36, 3.1632894411321006e-35, 1.5033958019782514e-34, 7.09563145582191e-34], [1.6482807007903047e-34, 7.746208013629305e-34, 3.6167295688758456e-33, 1.6780872353817137e-32, 7.739082606307619e-32, 3.5485222305084366e-31, 1.6180680297685028e-30, 7.339160067696146e-30, 3.312119470668537e-29, 1.4875976066095498e-28  …  2.0893996509174973e-41, 1.060360851338486e-40, 5.335520162561939e-40, 2.6623541459341753e-39, 1.317645726542395e-38, 6.469278830411349e-38, 3.1515263249937964e-37, 1.5236325928428463e-36, 7.311780136423824e-36, 3.4837080441289173e-35], [8.065627718219828e-36, 3.825566137065958e-35, 1.8022200922984585e-34, 8.434747541045195e-34, 3.92273074597557e-33, 1.8132565779399061e-32, 8.332752839030908e-32, 3.807868455581039e-31, 1.73079850696587e-30, 7.826903078023804e-30  …  9.216171485405112e-43, 4.729609438360341e-42, 2.406191945429452e-41, 1.2137680895997998e-40, 6.071756636017576e-40, 3.012600646006346e-39, 1.4828439199649907e-38, 7.241972585014038e-38, 3.5100201120970143e-37, 1.6886522430382952e-36], [3.896684012178807e-37, 1.8658089479486962e-36, 8.871331309238707e-36, 4.1894048468030666e-35, 1.965407616917856e-34, 9.161931424360413e-34, 4.244779411451108e-33, 1.9550528880864993e-32, 8.953658372650116e-32, 4.0783628474852235e-31  …  4.0075191516881777e-44, 2.0798009268802023e-43, 1.0699037870086317e-42, 5.456450208186343e-42, 2.7592193958524944e-41, 1.3837044428232268e-40, 6.882631945897021e-40, 3.396217932857324e-39, 1.662818120023425e-38, 8.079478055235599e-38]], h = [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  1.4999999999999991, 1.5000000000000018, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]  …  [1.5000000000000007, 1.5000000000000002, 1.4999999999999996, 1.5000000000000002, 1.5000000000000007, 1.5000000000000016, 1.5000000000000024, 1.5000000000000024, 1.5000000000000029, 1.5000000000000024  …  1.4999999999999998, 1.4999999999999996, 1.5000000000000004, 1.5000000000000007, 1.4999999999999998, 1.5, 1.5, 1.4999999999999996, 1.4999999999999998, 1.5000000000000002], [1.4999999999999993, 1.5000000000000002, 1.5000000000000004, 1.4999999999999998, 1.5000000000000002, 1.5000000000000016, 1.5000000000000016, 1.5000000000000016, 1.5000000000000022, 1.5000000000000027  …  1.5000000000000002, 1.4999999999999998, 1.5000000000000007, 1.5000000000000007, 1.5000000000000007, 1.5000000000000007, 1.5000000000000007, 1.5, 1.4999999999999996, 1.4999999999999996], [1.4999999999999998, 1.5000000000000002, 1.5, 1.5, 1.5, 1.5000000000000007, 1.5000000000000018, 1.5000000000000016, 1.5000000000000016, 1.5000000000000024  …  1.5000000000000002, 1.5, 1.4999999999999998, 1.5000000000000007, 1.5000000000000007, 1.5000000000000004, 1.5000000000000009, 1.5000000000000007, 1.5000000000000004, 1.5000000000000004], [1.5, 1.5000000000000004, 1.5000000000000004, 1.5000000000000007, 1.4999999999999998, 1.4999999999999998, 1.5, 1.5000000000000009, 1.5000000000000013, 1.5000000000000022  …  1.5, 1.5000000000000004, 1.5000000000000009, 1.5000000000000007, 1.5000000000000009, 1.5000000000000007, 1.5000000000000004, 1.5000000000000004, 1.5000000000000004, 1.5000000000000004], [1.5, 1.5000000000000002, 1.5000000000000007, 1.5000000000000004, 1.5000000000000004, 1.5000000000000007, 1.5000000000000004, 1.5000000000000004, 1.5000000000000004, 1.500000000000001  …  1.4999999999999993, 1.5000000000000002, 1.5000000000000002, 1.5000000000000004, 1.5000000000000007, 1.5000000000000009, 1.5000000000000004, 1.5000000000000004, 1.5000000000000002, 1.5000000000000002], [1.5000000000000002, 1.5000000000000002, 1.5, 1.5, 1.4999999999999998, 1.5, 1.5000000000000002, 1.500000000000001, 1.5000000000000013, 1.5000000000000004  …  1.4999999999999993, 1.5, 1.4999999999999998, 1.5000000000000002, 1.5000000000000004, 1.5000000000000013, 1.5000000000000007, 1.5000000000000004, 1.5000000000000013, 1.5000000000000002], [1.4999999999999998, 1.5000000000000002, 1.5, 1.5000000000000007, 1.5000000000000004, 1.5000000000000002, 1.4999999999999996, 1.5000000000000004, 1.5000000000000007, 1.5000000000000004  …  1.499999999999999, 1.4999999999999991, 1.4999999999999996, 1.5000000000000004, 1.4999999999999993, 1.5000000000000004, 1.5000000000000002, 1.5000000000000004, 1.5000000000000009, 1.5000000000000009], [1.5000000000000002, 1.4999999999999996, 1.5, 1.5000000000000007, 1.5, 1.5000000000000002, 1.4999999999999998, 1.5000000000000007, 1.5000000000000007, 1.5000000000000004  …  1.4999999999999991, 1.4999999999999993, 1.4999999999999996, 1.4999999999999996, 1.4999999999999998, 1.4999999999999998, 1.5000000000000009, 1.5000000000000004, 1.5000000000000004, 1.5000000000000009], [1.5000000000000009, 1.4999999999999993, 1.5000000000000007, 1.5000000000000004, 1.5, 1.5, 1.5, 1.4999999999999991, 1.4999999999999996, 1.5000000000000009  …  1.4999999999999987, 1.499999999999999, 1.4999999999999991, 1.499999999999999, 1.4999999999999998, 1.5, 1.5, 1.5, 1.5000000000000002, 1.5000000000000002], [1.5000000000000009, 1.5, 1.5, 1.5, 1.5000000000000002, 1.4999999999999993, 1.4999999999999993, 1.4999999999999993, 1.4999999999999993, 1.5  …  1.4999999999999987, 1.4999999999999987, 1.499999999999999, 1.4999999999999987, 1.4999999999999982, 1.4999999999999993, 1.5, 1.5, 1.5000000000000002, 1.5000000000000002]])
 (ℵ = [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]  …  [0.00013285133689763936, 0.00018837334860680806, 0.0002648790312067532, 0.0003693832967551969, 0.0005108994894978833, 0.000700889044630984, 0.0009537746851922591, 0.001287514918144405, 0.0017242335561452402, 0.0022908930862785017  …  2.5204185967666594e-6, 3.898294310559031e-6, 5.975546057111719e-6, 9.078408038043536e-6, 1.367095224117653e-5, 2.0406683476604686e-5, 3.01966063182632e-5, 4.429790612907177e-5, 6.442802229164784e-5, 9.290948138531075e-5], [7.050828866377653e-5, 0.00010125420309379661, 0.00014419744017069743, 0.00020365730668653513, 0.00028527698436726906, 0.00039635479839855183, 0.0005462337664611077, 0.0007467521135767589, 0.0010127550558612162, 0.0013626649663499235  …  1.1778678247271828e-6, 1.8451114763709734e-6, 2.8645113426952075e-6, 4.4076627908564e-6, 6.7223897046544025e-6, 1.0163045245920208e-5, 1.5231265521651485e-5, 2.2630126980554345e-5, 3.333524270717165e-5, 4.868694319756331e-5], [3.671801563631085e-5, 5.3397141084805945e-5, 7.70063254708524e-5, 0.00011013622968146392, 0.00015622667491149172, 0.0002197997026115454, 0.0003067408978316628, 0.00042463394694310314, 0.0005831523033357914, 0.0007945101116826865  …  5.40824050910432e-7, 8.579214258304755e-7, 1.3487823179407241e-6, 2.101681870068094e-6, 3.246014194492151e-6, 4.969576274513579e-6, 7.542245909337972e-6, 1.1348044045890013e-5, 1.6928042507244054e-5, 2.503710476767576e-5], [1.8769283920236967e-5, 2.7637494720237913e-5, 4.035692358563738e-5, 5.8442824537070206e-5, 8.393892079278125e-5, 0.0001195748436189245, 0.0001689606999298265, 0.0002368242573806481, 0.0003292959965061543, 0.00045424662368339166  …  2.4407140384410835e-7, 3.920229643042608e-7, 6.240416307342177e-7, 9.845737828543749e-7, 1.5397224228083156e-6, 2.386838441494838e-6, 3.6678938137500243e-6, 5.587918279958581e-6, 8.440128334268384e-6, 1.2639785109761327e-5], [9.42122444311968e-6, 1.4044835321869834e-5, 2.0763180474927386e-5, 3.0441355773443477e-5, 4.426411979662349e-5, 6.383842448387763e-5, 9.132289317299543e-5, 0.00012958925072631376, 0.00018242095588132551, 0.0002547542969689847  …  1.0832317803254931e-7, 1.76115901432981e-7, 2.8381230086611615e-7, 4.5332819793503943e-7, 7.177272270581747e-7, 1.1264058093754888e-6, 1.7524465629204176e-6, 2.7029323331883177e-6, 4.133252854551828e-6, 6.266728489249752e-6], [4.645329654268873e-6, 7.01020100019593e-6, 1.0490891813288884e-5, 1.556994162830022e-5, 2.2918103761629896e-5, 3.3458915423033414e-5, 4.845189616588457e-5, 6.959845835048195e-5, 9.917506875605625e-5, 0.00014019855501425972  …  4.7382613717679086e-8, 7.786469037692502e-8, 1.269532783910335e-7, 2.0523260811282905e-7, 3.2890196609223535e-7, 5.225089470644243e-7, 8.228906525078614e-7, 1.2847948255088658e-6, 1.9888099211703318e-6, 3.05243066403107e-6], [2.250765536728034e-6, 3.437916946053215e-6, 5.207497602170109e-6, 7.822682214147513e-6, 1.16546540646764e-5, 1.7222024516063433e-5, 2.5242596822485933e-5, 3.670055704662502e-5, 5.2932667144534773e-5, 7.573751064317242e-5  …  2.073143565754871e-8, 3.407067281373876e-8, 5.597975651037492e-8, 9.145261048164917e-8, 1.4825830845673055e-7, 2.3834569173663128e-7, 3.7990351166839453e-7, 6.003478131312386e-7, 9.406065692610427e-7, 1.461198620732435e-6], [1.0720160159432735e-6, 1.6571567738490504e-6, 2.5403651886538673e-6, 3.862091234930112e-6, 5.8232641938588125e-6, 8.70866837500816e-6, 1.2918196428514548e-5, 1.9008171156215352e-5, 2.7745364462291157e-5, 4.0176807089477034e-5  …  1.0002969427821593e-8, 1.533727161759615e-8, 2.4697001050814237e-8, 4.034114824317163e-8, 6.589067930234742e-8, 1.070277963447405e-7, 1.7254392577092675e-7, 2.7588743618726494e-7, 4.37421416855086e-7, 6.876802598882871e-7], [5.021029742893842e-7, 7.853961000770366e-7, 1.2183278148831826e-6, 1.8742918078249042e-6, 2.8597600652218125e-6, 4.327770017540596e-6, 6.496269780761731e-6, 9.672802776014216e-6, 1.428736098121199e-5, 2.093562824472219e-5  …  7.903532915295078e-9, 8.815423235504647e-9, 1.19968723835996e-8, 1.8306381093832385e-8, 2.9306234580565583e-8, 4.758690860435066e-8, 7.727763946266267e-8, 1.2482309443363473e-7, 2.0014109916008042e-7, 3.18323299115033e-7], [2.3140727208608616e-7, 3.661532135462114e-7, 5.7463968252801e-7, 8.944408964111791e-7, 1.3808198551368044e-6, 2.114310176664285e-6, 3.211201219088133e-6, 4.837885841088648e-6, 7.230286024970139e-6, 1.0719872651654277e-5  …  1.3960690268471914e-8, 1.0531870437388545e-8, 9.416053329606381e-9, 1.054557843746341e-8, 1.4330227106684214e-8, 2.1764622951534455e-8, 3.464477508523095e-8, 5.593129972670833e-8, 9.03178821505272e-8, 1.4509139563986718e-7]], h = [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]  …  [1.4999999999999996, 1.5000000000000002, 1.5000000000000004, 1.5000000000000004, 1.5000000000000002, 1.5000000000000007, 1.5000000000000007, 1.5000000000000009, 1.500000000000001, 1.500000000000001  …  1.4999999999999984, 1.499999999999999, 1.4999999999999996, 1.4999999999999993, 1.4999999999999993, 1.4999999999999991, 1.4999999999999993, 1.4999999999999998, 1.4999999999999996, 1.4999999999999991], [1.4999999999999996, 1.4999999999999996, 1.5, 1.5000000000000002, 1.5, 1.5000000000000002, 1.5000000000000002, 1.5, 1.5000000000000007, 1.500000000000001  …  1.4999999999999993, 1.499999999999999, 1.499999999999999, 1.4999999999999991, 1.4999999999999993, 1.4999999999999991, 1.4999999999999996, 1.4999999999999998, 1.4999999999999998, 1.4999999999999998], [1.4999999999999998, 1.4999999999999993, 1.4999999999999996, 1.4999999999999998, 1.5000000000000004, 1.5000000000000004, 1.5, 1.5, 1.5000000000000004, 1.5000000000000004  …  1.4999999999999996, 1.4999999999999998, 1.4999999999999996, 1.499999999999999, 1.499999999999999, 1.499999999999999, 1.4999999999999991, 1.4999999999999996, 1.5, 1.5], [1.5, 1.4999999999999998, 1.4999999999999998, 1.4999999999999996, 1.5, 1.4999999999999998, 1.4999999999999998, 1.5, 1.5000000000000002, 1.5000000000000002  …  1.4999999999999993, 1.4999999999999993, 1.4999999999999998, 1.4999999999999996, 1.4999999999999993, 1.4999999999999984, 1.4999999999999991, 1.4999999999999993, 1.4999999999999991, 1.4999999999999998], [1.4999999999999993, 1.5000000000000002, 1.5000000000000002, 1.5, 1.4999999999999998, 1.4999999999999998, 1.4999999999999998, 1.4999999999999998, 1.4999999999999996, 1.4999999999999998  …  1.4999999999999996, 1.4999999999999996, 1.4999999999999996, 1.4999999999999996, 1.4999999999999993, 1.4999999999999991, 1.4999999999999987, 1.4999999999999991, 1.4999999999999993, 1.4999999999999993], [1.4999999999999993, 1.4999999999999998, 1.4999999999999998, 1.4999999999999998, 1.4999999999999998, 1.4999999999999998, 1.4999999999999998, 1.5, 1.5, 1.5000000000000002  …  1.4999999999999987, 1.4999999999999991, 1.4999999999999993, 1.4999999999999998, 1.5, 1.4999999999999991, 1.4999999999999991, 1.4999999999999984, 1.4999999999999993, 1.4999999999999993], [1.4999999999999993, 1.4999999999999996, 1.4999999999999991, 1.4999999999999996, 1.4999999999999996, 1.4999999999999996, 1.4999999999999996, 1.5, 1.5000000000000002, 1.5000000000000002  …  1.499999999999999, 1.4999999999999991, 1.4999999999999993, 1.4999999999999993, 1.4999999999999998, 1.4999999999999998, 1.4999999999999993, 1.499999999999999, 1.4999999999999991, 1.4999999999999993], [1.4999999999999991, 1.4999999999999993, 1.4999999999999993, 1.4999999999999996, 1.4999999999999996, 1.4999999999999998, 1.4999999999999991, 1.4999999999999993, 1.4999999999999998, 1.5  …  1.4999999999999991, 1.499999999999999, 1.499999999999999, 1.4999999999999991, 1.4999999999999991, 1.4999999999999993, 1.4999999999999998, 1.4999999999999993, 1.499999999999999, 1.4999999999999987], [1.4999999999999991, 1.4999999999999991, 1.499999999999999, 1.4999999999999991, 1.4999999999999993, 1.4999999999999991, 1.4999999999999996, 1.4999999999999991, 1.4999999999999993, 1.4999999999999996  …  1.5, 1.4999999999999991, 1.4999999999999991, 1.499999999999999, 1.4999999999999987, 1.4999999999999991, 1.4999999999999993, 1.4999999999999993, 1.4999999999999998, 1.4999999999999991], [1.4999999999999993, 1.499999999999999, 1.499999999999999, 1.4999999999999991, 1.4999999999999991, 1.4999999999999991, 1.4999999999999996, 1.4999999999999993, 1.4999999999999996, 1.4999999999999998  …  1.5, 1.5, 1.4999999999999998, 1.4999999999999993, 1.4999999999999993, 1.4999999999999993, 1.4999999999999993, 1.4999999999999991, 1.4999999999999991, 1.4999999999999993]])

Visualizing

The top panel shows the concentration, the bottom panel the recovered thickness where the concentration is above 10⁻³. The thickness axis is clipped at 3 m: the exact answer never leaves 1.5 m, so anything that runs off the top of the panel is the h = 𝓋/ℵ recovery failing at the ice edge.

x = Array(xnodes(grid, Center()))

frames = length(solutions[1].ℵ)
iter = Observable(1)

colors = [:seagreen, :crimson, :royalblue, :darkorange]

fig = Figure(size = (1000, 620))

axℵ = Axis(fig[1, 1], ylabel = "concentration ℵ", limits = ((0, 1), (-0.05, 1.15)))
axh = Axis(fig[2, 1], ylabel = "thickness h (m)", xlabel = "x", limits = ((0, 1), (0, 3)))

hlines!(axh, [1.5], color = :black, linestyle = :dash, linewidth = 1)

# The exact solution is the initial hat, translated and wrapped around the periodic domain.
time(i) = (i - 1) * save_interval * Δt
ℵexact = @lift([iced(mod(ξ - U * time($iter), 1)) ? 1.0 : 0.0 for ξ in x])
hexact = @lift([iced(mod(ξ - U * time($iter), 1)) ? 1.5 : 0.0 for ξ in x])

lines!(axℵ, x, ℵexact, color = (:black, 0.35), linewidth = 6, label = "exact")
lines!(axh, x, hexact, color = (:black, 0.35), linewidth = 6)

# `h` is masked where the ice vanishes, so the schemes' different trailing tails do not show up as shifts
ice_threshold = 1e-3

for (n, configuration) in enumerate(configurations)
    ℵn = @lift(solutions[n].ℵ[$iter])
    hn = @lift([solutions[n].ℵ[$iter][k] > ice_threshold ? solutions[n].h[$iter][k] : NaN
                for k in eachindex(x)])

    lines!(axℵ, x, ℵn, color = colors[n], linewidth = 2, label = configuration.name)
    lines!(axh, x, hn, color = colors[n], linewidth = 2)
end

axislegend(axℵ, position = :rt, framevisible = false)

title = @lift(@sprintf("Travelling top hat — t = %.2f revolutions, peak h = %.1f m",
                       U * time($iter),
                       maximum(maximum(solutions[n].h[$iter]) for n in 1:length(configurations))))

Label(fig[0, 1], title, fontsize = 18, tellwidth = false)

CairoMakie.record(fig, "one_dimensional_ice_advection.mp4", 1:frames, framerate = 12) do i
    iter[] = i
end

The concentration behaves as expected in every case: the hat translates, its edges smear — a lot for first-order upwinding — and it stays between zero and one.

The thickness does not. Flux-form advection reconstructs ℵ and 𝓋 = ℵh independently, so their ratio is a weighted average whose numerator and denominator come from different reconstructions. Where the concentration collapses at the ice edge, the denominator vanishes but the numerator does not vanish with it, and the recovered thickness leaves the physical range — dramatically so under Forward Euler, and the same mechanism is what produces spurious thick ice in a full simulation.

Incremental remapping integrates the transported area and the transported content over the same swept region, so

\[h⁺ = \frac{Σ_k w_k ℵ̃_k h̃_k}{Σ_k w_k ℵ̃_k}\]

is a weighted average of the reconstructed thickness with non-negative weights. It is bounded by construction, and here it holds 1.5 m to round-off.

This page was generated using Literate.jl.