Details on the Level Set Method

reinitialize_ϕ_HCR!(ϕ, dom::Domain; maxsteps = 50, tol=1e-4, err_reg=:B)

Reinitialize ϕ throughout the domain.

Implementation of Eq. 22 in Hartmann 2010, scheme HCR-2.

Checks L∞ error (of |∇ϕ|=1) against tol either in band around interface (err_reg=:B) or throughout domain (err_reg=:all), and ends iteration early if tolerance is met.

reinitialize_ϕ_HCR(ϕ, dom::Domain)

Thin wrapper on reinitialize_ϕ_HCR! to avoid mutating.

extrap_v_fastmarch(v_front, u, dom::Domain)

Compute an extrapolated velocity field in place. v_front should be nonzero only on Γ⁺, positive side of interface.

Using fast marching, instead of the PDE-based approach, to get second order accuracy more easily.

function identify_Γ(ϕ, dom::Domain)

Identify cells on the sublimation front (interface), returning as a Matrix::Bool.

Γ_cells(ϕ, dom::Domain)

Compute findall(identify_Γ(ϕ, dom)). (That's the whole implementation.)

𝒢_weno(ϕ, ir::Int, iz::Int, dom::Domain)
𝒢_weno(ϕ, ind::CartesianIndex{2}, dom::Domain)
𝒢_weno(ϕ, ir::Int, iz::Int, ϕ0sign, dom::Domain)
𝒢_weno(ϕ, ind::CartesianIndex{2}, ϕ0sign, dom::Domain)

Compute the norm of the gradient by Godunov's scheme with WENO differences (wenodiffs_local).

If supplied, ϕ0sign is used in Godunov's scheme, rather than the current sign of ϕ. Described in (Hartmann et al., 2009), eq. 6 to eq. 9. Boundary treatment of ghost cells handled by get_or_extrapolate_ϕ.

𝒢_weno_all(ϕ, [dϕdx_all,] dom::Domain; signs=nothing)

Compute the norm of the gradient of ϕ throughout domain by Godunov's scheme with WENO derivatives. If signs of ϕ should refer to some prior state, can be provided. If dϕdx_all is provided, then is not recomputed internally.

Internally, calls 𝒢_weno on all computational cells.

dϕdx_all_WENO(ϕ, dom)

Compute (∂ϕ/∂r west, ∂ϕ/∂r east, ∂ϕ/∂z south, ∂ϕ/∂z north) using WENO derivatives.

Implemented by computing WENO derivatives for each cell separately, which is a little wasteful. Beyond the boundaries of domain, ϕ is extrapolated according to get_or_extrapolate_ϕ.

wenodiffs_local(u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, dx)

Compute one-sided finite differences, using Jiang and Peng's WENO approximation (Jiang and Peng, 2000).

A relatively easy-to-read reference is §3.2 of (Hartmann et al., 2009) .

identify_B(Γc::Vector{CartesianIndex{2}}, dom::Domain; extra=0)
identify_B(Γ_field::Matrix{Bool}, dom::Domain)
identify_B(ϕ::Matrix{Float64}, dom::Domain)
identify_B(ϕ::Any, dom::Domain)

Return a field of bools identifying the band around the interface.

The width in the band around Γ is specified by the fields bwr and bwz, which represent number of cells in the 𝑟 and 𝑧 directions respectively. extra will tack on extra cells, if in some (but not all) places you need a larger band than the default.

sdf_err_L1(ϕ, dom; region=:B)

Compute the L1 norm of the error in the Eikonal equation |∇ϕ| = 1, on the given region (either :B or :all).

sdf_err_L∞(ϕ, dom; region=:B)

Compute the L∞ norm of the error in the Eikonal equation |∇ϕ| = 1, on the given region (either :B or :all).

compute_icevol_H(ϕ, dom)

Compute the volume of the ice, using a discrete Heaviside function. Calls compute_discrete_H, which is implemented according to Min and Gibou 2008.

compute_discrete_H(ϕ, dom)

Compute the discrete Heaviside H across the domain, for use in volume integrals.

Implements the discrete Heaviside of Min & Gibou, 2008.

compute_local_H(ij::CartesianIndex, ϕ, dom)

Return the discrete Heaviside for level set ϕ at location ij.

Implementation of the final expression for a discrete Heaviside in Min & Gibou, 2008.

compute_icesurf_δ(ϕ, dom)

Compute the surface area of the ice, using a discrete Dirac delta function. Calls compute_discrete_δ, which is implemented according to Min and Gibou 2008.

compute_discrete_δ(ϕ, dom)

Compute the discrete Dirac δ throughout the domain, for use in surface integrals. Implements the discrete delta of Min and Gibou 2008.

compute_local_δ(ij::CartesianIndex, ϕ, dom)

Return the discrete delta for level set ϕ at location ij.

Implementation of the final expression for a discrete delta in Min and Gibou 2008.

make_ϕ0(ϕtype::Symbol, dom::Domain; ϵ=1e-4)

Return a ϕ0 with appropriate size for the passed Domain.

Parameter ϵ * sqrt(dom.rmax*dom.zmax) is added to ensure that interface is within domain. Currently allowed setups:

• :top, :flat – interface at zmax*(1 - ϵ)
• :rad, :cyl – interface at rmax*(1 - ϵ)
• :box – interface at both zmax(1-ϵ) and rmax(1-ϵ)
• :cone – interface a line decreasing in r
• :midflat – interface at zmax*0.5
• :ell_bub – ellipse in center of vial, separated from boundaries
• :circ – circle at r=0, z=0
• :tinycirc – circle at r=0, z=0, very small radius