Unfiltered List of All Documented Functions

An instance of PhysicalProperties with default values.

Stores information on the grid and domain size for simulation.

Constructors:

Domain(config::Dict) 
Domain(config::NamedTuple) 
Domain(nr, nz, rmax, zmax) 
Domain(nr, nz, rmax, zmax, bwfrac)
Domain(nr, nz, rmin, rmax, zmin, zmax)
Domain(nr, nz, rmin, rmax, zmin, zmax, bwfrac)

If not given, rmin and zmin default to 0.0, while bwfrac defaults to 0.2. Note that nr, nz, bwr, and bwz are assumed to be integers, while rmin, rmax, zmin, and zmax are assumed to be floats

All Fields: (given for 𝑟, likewise named for 𝑧)

• nr - number of grid points in 𝑟
• rmin, rmax - range of domain in 𝑟
• bwfrac - fraction of domain included in band around interface
• rgrid - vector of locations of grid points in 𝑟
• dr - grid spacing in 𝑟
• dr1 - 1/dr
• dr2 - 1/dr^2
• bwr - =ceil(Int, bwfrac*nr) integer width of band around interface in which level set is treated
• ntot - total number of grid points = nr*nz

LD{T}

A "little difference", to make Godunov's scheme in 𝒢_weno easier to read. For a = LD(x::T),

• a.p = max(x, 0.0)
• a.m = min(x, 0.0)

LevelSet(phi, dom::Domain)

A small struct, used in dispatch for plotting.

The plot recipe can be variously called as:

plot(LevelSet(phi, dom), reflect=true)
plot(LevelSet(phi, dom), reflect=false)
plot!(LevelSet(phi, dom), reflect=true)
plot!(LevelSet(phi, dom), reflect=false)

An immutable type for representing a slate of physical properties and constants needed.

A default constructor is provided which uses all defaults (borosilicate glass, water, ice), but thanks to Parameters.jl you can call with just values that need to be changed, e.g.

PhysicalProperties() # all defaults
PhysicalProperties(ρ_vw = 1u"kg/m^3", cp_vw = 1u"J/kg/K", εpp_vw = 1e-1) # adjust vial wall properties

All properties stored here are:

• ρ_vw: Density of vial wall (defaults to borosilicate glass). Default: LevelSetSublimation.ρ_gl

• cp_vw: Heat capacity of vial wall (defaults to borosilicate glass). Default: LevelSetSublimation.cp_gl

• εpp_vw: Dielectric loss coefficient of vial wall (defaults to borosilicate glass). Default: LevelSetSublimation.εpp_gl

• ρf: Density of frozen material (taken as water ice). Default: ρ_ice

• Cpf: Heat capacity of frozen material (taken as water ice). Default: cp_ice

• kf: Thermal conductivity of frozen material. Default: LevelSetSublimation.kf

• Mw: Molecular weight of sublimating species (defaults to water). Default: 0.018 * u"kg/mol"

• μ: Gas phase viscosity (in rarefied regime) of sublimating species (defaults to water). Default: LevelSetSublimation.μ

• ΔH: Heat of sublimation of sublimating species (defaults to water); give as positive number. Default: LevelSetSublimation.ΔH

• εppf: Dielectric loss coefficient of frozen layer (defaults to water ice). Default: LyoPronto.εppf

• εpp_d: Dielectric loss coefficient of dry layer (defaults to 0). Default: 0.0

• R: Universal gas constant. Kept as a constant so that it gets nondimensionalized with the rest. Default: u"R"

• ϵ0: Vacuum permittivity. Kept as a constant so that it gets nondimensionalized with the rest. Default: u"ϵ0" |> u"F/m"

TimeConstantProperties(ϵ, l, κ, Rp0, kd, Kvwf, m_v, A_v, B_d, B_f, B_vw)

A struct for holding physical properties which are likely to change from case to case.

No default constructor is provided by intention–all of these parameters should be at least considered before running a simulation.

• ϵ

• l

• κ

• Rp0

• kd

• Kvwf

• m_v

• A_v

• B_d

• B_f

• B_vw

TimeVaryingProperties(f_RF, P_per_vial, Tsh, pch, Kshf)

A struct for holding controlled parameters that may change over time.

To get the value of all those parameters at a given time, call the struct with a time, and it will evaluate each field at that time and provide a TimeVaryingPropertiesSnapshot.

No default constructor is provided by intention–all of these parameters should be at least considered before running a simulation.

Each should be passed as a callable, which returns the value of the parameter as a function of time, with the exception of Kshf which is a function of pressure. See the RampedVariable and RpFormFit types from LyoPronto, which are intended to make this more convenient.

• f_RF

• P_per_vial

• Tsh

• pch

• Kshf

TimeVaryingPropertiesSnapshot(f_RF, P_per_vial, Tsh, pch, Kshf)

A struct for holding a snapshot of the controlled parameters that may change over time.

This is meant to be constructed by calling an instance of the TimeVaryingProperties type.

• f_RF

• P_per_vial

• Tsh

• pch

• Kshf

arrows(Vf, dom::Domain)

Add "quiver" to current plot, with velocity field Vf of shape (nx, ny, 2).

calc_H0(ϕ0, ϕ1, ϕ2, P0, P1, P2)

Implementation of Min & Gibou 2008, Table 1.

calc_Tf_res(t, sim)

Get Tf at a given time point, regardless of simulation's time integration.

calc_fillvol(dom)

calc_params_at_t(t, config)

calc_psub_highorder(T)

Compute pressure (in Pascals) of sublimation at temperature T in Kelvin

From Feistel and Wagner, 2006

calc_rij_Sij(ϕ, Γ, C, dom::Domain)

Compute rij, neighbors Sij for each cell in Γ.

Implementation of eq. 19b from Hartmann 2010, scheme HCR-2

calc_uTfTp_res(t, sim)

calc_δ0(ϕ0, ϕ1, ϕ2, P0, P1, P2)

Implementation of Min & Gibou 2008, Table 1.

compare_lyopronto_res(ts, sim)

compare_lyopronto_res(sim)

function compute_Qice(ϕ, T, p, dom::Domain, params)

Compute the total heat input into frozen & dried domains. Also passes glass-ice heat as separate return.

See p. 112-113 from paperlike notes. At pseudosteady conditions, all heat getting added to the system goes to the frozen domain, so we don't actually need to treat different areas of the surface. In addition, the sublimation flux is simply evaluated at the top cake surface.

compute_Qice_nodry(u, T, dom::Domain, params)

Compute the total heat input into frozen domain from volumetric, shelf, and glass. Contrast with compute_Qice_noflow and compute_Qice, which include heat to dried domain.

compute_Qice_noflow(u, T, dom::Domain, params)

Compute the total heat input into frozen & dried domains, assuming mass flow is zero. Also passes glass-ice heat as separate return.

compute_Tderiv(u, T, ir::Int, iz::Int, dom::Domain, params)

Compute (∂T/∂r, ∂T/∂z) at point (ir, iz), given system state u, T.

Quadratic ghost cell extrapolation (into frozen domain), second order finite differences, for T.

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_discrete_δ(ϕ, dom)

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

compute_frontvel_fixedspeed(v0, ϕ, dom::Domain)

Compute speed v0 times normal vector on Γ⁺, positive side of interface.

compute_frontvel_heat(u, T, p, dom::Domain, params; debug=false)

Generate an empty velocity field and compute velocity on Γ⁺ (i.e. cells on Γ with ϕ>0).

compute_frontvel_mass(ϕ, Tf, T, p, dom::Domain, params; debug=false)

Generate an empty velocity field and compute velocity on Γ⁺ (i.e. cells on Γ with ϕ>0).

compute_icegl_area(ϕ, dom::Domain)

Compute area of where ice meets radial outer surface of vial.

function compute_iceht_bottopcont(ϕ, dom)

Compute the height of ice at each point radially, as well as identify whether it touches system boundary or not.

Returns (heights, bottom_contact, top_contact), where

• heights is a vector of floats
• bottom_contact, top_contact are vectors of bools

compute_icesh_area(ϕ, dom::Domain)

Compute area of where ice meets bottom surface of vial.

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_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_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_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.

compute_pderiv(u, T, p, ir::Int, iz::Int, dom::Domain, params)

Compute ∂p/∂r, ∂p/∂z at point (ir, iz), given system state u, T, p.

Quadratic ghost cell extrapolation (into frozen domain), second order finite differences, for p.

The distinction in discretization between this and compute_Tderiv is essentially just the boundary treatments.

compute_topmassflux(ϕ, T, p, dom::Domain, params)

Compute total mass flow through top of the cake (that is, mass flux integrated across top surface).

cond_end(u, t, integ)

Changes from positive to negative when the minimum value of ϕ is within min(dom.dz, dom.dr)/4, which is when drying is essentially complete.

dudt_heatmass!(du, u, p, t)

Evaluate local time rate of change for u and put results in du.

Splitting u and du into ϕ, Tf, and Tvw is handled by ϕ_T_from_u and ϕ_T_from_u_view.

Parameters p assumed to be (dom::Domain, params) This is a right-hand-side for ∂ₜϕ = -v⋅∇ϕ, where v = (vr, vz) is evaluated by computing and extrapolating front velocity using compute_frontvel_mass.

dudt_heatmass(u, dom::Domain, params, t)
dudt_heatmass(u, config::Dict, t=0.0)

Compute the time derivative of u with given (nondimensional) parameters.

Wraps a call on dudt_heatmass!, for convenience in debugging and elsewhere that efficiency is less important

doc

doc

dudt_heatmass_params(u, config)

Compute the time derivative of u with given parameters, and also return dom and params associated with the given config.

Wraps a call on dudt_heatmass, for convenience in debugging and elsewhere that efficiency is less important

dudt_heatonly!(du, u, p, t)

Evaluate local time rate of change for u and put results in du.

This function leaves Tf and Tvw untouched, since there isn't a way to govern their dynamics without mass transfer.

u and du are both structured according to

This is a right-hand-side for ∂ₜϕ = -v⋅∇ϕ, where v = (vr, vz) is evaluated by computing and extrapolating front velocity using compute_frontvel_mass.

dudt_heatonly(u, dom::Domain, params)
dudt_heatonly(u, config)

Compute the time derivative of u with given parameters.

Wraps a call on dudt_heatonly!, for convenience in debugging and elsewhere that efficiency is less important

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_ϕ.

dϕdx_all_WENO_loc(ϕ, 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_ϕ.

eval_b(T, p, params)

Compute transport coefficient b as a function of space, given T, p, and params.

params should have the following fields: l, κ, R, Mw, μ. l, and κ may be passed as scalars (and assumed as spatially uniform) or arrays (describing value throughout space, should match Domain dimensions).

When the simulation is started, all these values are converted to SI units and passed accordingly, so in practice there are no units to track.

If κ=0, no spatial variation due to pressure occurs.

eval_b_loc(T, p, ir, iz, params)

Locally evaluate transport coefficient (indexes into spatially varying l and κ if necessary).

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.

fastmarch_v!(vf, acc, locs, ϕ, dom)

Mutate vf and acc to extrapolate vf from interface by fast marching.

Cells are calculated in order of increasing |ϕ|. Uses matrix of bools acc (denoting accepted cells) to determine cells to use in extrapolation, so if you determine one side first, will get used in the derivatives for the other side.

freshplot(dom::Domain)

Generate an empty plot, scaled to show only dom.

gen_anim(config)
gen_anim(config, var)
gen_anim(config, var, casename)

gen_sumplot(config)
gen_sumplot(config, var)
gen_sumplot(config, var, casename)

get_SA(res)

Compute surface area over time for the given simulaiton results.

Returns (ts, SA_t)

get_eff_Rp(sim)

get_or_extrapolate_ϕ(ϕ, ind, stencil)

Either retrieve ϕ at ind on the given stencil, or extrapolate and return domain + extrapolated values.

This extrapolation is on a uniform grid, using a quadratic extrapolant. Define Lagrange interpolant with last three points in domain, then extrapolate outside the domain to make ghost points. Ghost points are then a linear combination of points from domain. A matrix form just keeps this compact and efficient; extrap_ϕ_mat in the source is just the matrix representing this extrapolation.

get_subf_r(ϕ, dom)

Compute the average 𝓇 position of the sublimation front.

get_subf_z(ϕ, dom)

Compute the average 𝑧 position of the sublimation front.

get_t_Tf(sim)

get_t_Tf_subflux(sim)

heat(field, dom::Domain)

Return a heatmap of field, scaled to fit on dom.

This function generates a new plot.

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.

function identify_Γ(ϕ, dom::Domain)

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

interp_saved_Tf(saved_Tf)

Interpolate a set of SavedValues of Tf. Useful so that I use the same interpolation everywhere it shows up.

make_M1_properties()

Returns physical parameters for the mannitol experimental case which was used to originally develop and validate this model.

Useful for testing.

make_u0_ndim(init_prof, Tf0, Tvw0, dom)

Set up a vector of dynamic variables as initial state for simulation.

Structure of vector u: ComponenentArray with fields ϕ, Tf, and Tvw. ϕ is a 2D array of the same size as the domain, Tf is a 1D array of size dom.nr, and Tvw is a scalar.

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

markcells(cells, dom::Domain; lab="", c=:white)

Add a star marker for each CartesianIndex in cells to the current plot.

markfront(phi, dom::Domain; lab="", c=:white)

Add a star marker to Γ (front) cells based on ϕ for current plot, in mutating fashion.

Internally calls Γ_cells.

params_nondim_setup(params_dim)
params_nondim_setup(base_d, tcp_d, tvp_d)

plot_cylheat(T, dom::Domain)

In a mutating fashion, add a "cylindrical" heatmap of T to the current plot. ("Cylindrical" meaning reflected across x=0 axis.)

function plot_frontvel(ϕ, T, dom::Domain)

Calculate, then plot the front velocity given ϕ and T.

Meant for debugging, mostly. Scales all velocity arrows to have length 0.5. Generates a freshplot().

plotframe(t::Float64, sim; maxT=nothing, heatvar=:T)

Unpack simulation results and plot the state at time t.

heatvar = :T or =:ϕ or =:p decides whether temperature, level set function, or pressure is plotted as colors. If given, maxT sets an upper limit for the associated colorbar.

reinit_wrap(integ)

Thin wrapper to reinitialize the state of the level set function.

Uses a tolerance of 0.02 for error in norm of gradient (which should be 1) in the region B (band around the interface). If verbose=true, logs the error in signed distance function before and after reinit, as well as current dried fraction

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.

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

Reinitialize ϕ throughout the domain.

Implementation of Eq. 22 in Hartmann 2010, scheme HCR-2. Augmented by (Della Rocca and Blanquart, 2014)'s treatment of contact line blind spots.

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.

resultsanim(sim, casename; seconds_length=5, heatvar=:T)

Generate a .gif of the given sim, with filename casename_heatvar_evol.gif.

Pass either :p or :T as heatvar. Passing ϕ will probably cause filename problems

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).

sim_from_dict(config; tf=1e5, verbose=false)

Given a simulation configuration config, run a simulation.

Maximum simulation time (not CPU time, but simulation time) is specified by tf, which is in seconds; 1e5 is 270 hours. (No matter how much ice is left at that point, simulation will end.) verbose=true will put out some info messages about simulation progress, i.e. at each reinitialization.

All passed parameters should have Unitful unit annotations, so they can be in whatever units are convenient; inside this function, they will be converted to SI marks then have units stripped off (for easier numerical implementation). The results will be in SI units, so times in seconds, temperature in Kelvin, pressure in Pa, etc. Employ Unitful to do unit conversions after simulation if necessary.

config must have the following fields:

• fillvol, fill volume with Unitful units (e.g. 3u"mL")
• vialsize, a string (e.g. "2R") giving vial size
• paramsd, a Tuple{PhysicalProperties, TimeConstantProperties, TimeVaryingProperties} (i.e. a tuple of the three objects). See PhysicalProperties, TimeConstantProperties, TimeVaryingProperties.

The following fields have default values and are therefore optional:

• simgridsize, a tuple giving number of grid points to use for simulation. Defaults to (51, 51).
• Tf0, an initial ice temperature. Defaults to Tsh(0) if not provided
• Tvw0, an initial glass temperature. Defaults to Tf0.
• time_integ: defaults to Val(:exp_newton); other options are Val(:dae), Val(:dae_then_exp) and Val(:ode_implicit). Among these three, slightly different problem formulations are used (explicit ODE with internal Newton solve, DAE, and implicit ODE).
• init_prof, a Symbol indicating a starting profile (from make_ϕ0)

The three problem formulations have different advantages; in my testing, the DAE formulation and implicit formulation tend to run faster and more closely reflect the problem structure, but they run into instabilities when the sublimation front peels away from the wall, whereas the explicit ODE formulation can jump over that point.

sim_from_u0(u0, t0, config; tf=1e6, verbose=false)

Wrapped by sim_from_dict; useful on its own if you want to start from partway through a simulation.

sim_heatonly(config; tf=1e5, verbose=false)

Given a simulation configuration config, run a simulation.

This simulation is stripped-down: no mass transfer, no variation in ice & glass temperature

Maximum simulation time is specified by tf. verbose=true will put out some info messages about simulation progress, i.e. at each reinitialization.

CURRENTLY OUT OF DATE with parameter input structures.

solve_T(u, Tf, dom::Domain, params)

Compute 2D axisymmetric T profile, returning ghost cell values, for given state u and Tf.

Neumann boundary conditions on all rectangular boundaries; Dirichlet on zero-level set.

This implementation uses second-order finite differences, with linear extrapolation into Ω⁻. Coefficients are all hard-coded here, unfortunately. (For more on extrapolation, see Gibou et al., 2002, "Second-Order-Accurate ... Poisson ... ") Neumann boundaries use a ghost point & BC to define ghost cell, then use same stencil as normal. Coefficients computed in gfm_extrap.ipynb, using Sympy. (For higher order, see Gibou and Fedkiw, 2005, "A fourth order accurate discretization ... Laplace ... ")

solve_p(u, T, dom::Domain, params[, p0]; maxit=20, reltol=1e-6) where G<:AbstractArray

Iteratively compute the pressure profile for system state u, Tf, and T.

There is a weak nonlinearity in the system, since the mass transfer coefficient b depends partially on pressure. To treat this, use a guessed pressure p0 (which, if not provided, is set everywhere to chamber pressure) to compute b, then perform a linear solve for p using solve_p_given_b. At this point, recompute b, then recompute p.

In preliminary testing, this usually converges within 5 or 10 iterations. Usually if it doesn't converge, it is because temperatures are outside the expected range, yielding crazy sublimation pressures. (Occasionally it means I incorrectly wrote a finite difference somewhere.)

solve_p_given_b(ϕ, b, Tf, dom::Domain, params)

Compute 2D axisymmetric pseudosteady pressure profile for given values of level set function ϕ, temperature T, and transport coefficient b.

b is a dusty-gas transport coefficient for the pressure, which can vary spatially. Homogeneous Neumann boundary conditions at r=0, r=R, z=0; Dirichlet on zero-level set (p=p_sub), Robin at top (dp/dz = (pch-p)/Rp0). params should have fields:

• Rp0 : zero-thickness resistance offset, often written R0 in lyo literature
• pch : chamber (or vial) pressure at top surface

This implementation uses second-order finite differences, with linear extrapolation into Ω⁻. (For details, see (Gibou and Fedkiw, 2005).) Coefficients are all hard-coded here, unfortunately. Neumann boundaries use a ghost point & BC to define ghost cell, then use same stencil as normal. Coefficients computed in gfm_extrap.ipynb, using Sympy.

summaryT(sim; layout, clims, tstart, tend)

Plot temperature fields at several time instants throughout the simulation.

summaryplot(sim; layout=(3,2), heatvar=:T)

Return a 2x3 plot of simulation results from start to finish.

sim should have a field "sol" , which is passed to get_ϕ(sol, t, dom::Domain) . heatvar determines what is plotted as a heatmap in the results (:T or :ϕ, currently.)

virtual_thermocouple(sim)

weno_Φ(c, d, e, f)

Return a weighted sum of finite differences for a WENO approximation. Defined in §3.2 of (Hartmann et al., 2009) .

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) .

Γ_cells(ϕ, dom::Domain)

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

𝒢_1st(ϕ, i, j, dom::Domain)

Compute the norm of the gradient of ϕ at point i, j by Godunov's scheme to first-order accuracy.

p. 6830 of Hartmann 2008, 10th page of PDF

𝒢_1st_all(ϕ, dom::Domain)

Compute the norm of the gradient of ϕ throughout domain by Godunov's scheme to first-order accuracy.

Internally, calls 𝒢_1st on all computational cells.

𝒢_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.

𝒢_weno_all_old(ϕ, dom::Domain)

Compute the norm of the gradient of ϕ throughout domain by Godunov's scheme with WENO derivatives.

Internally, calls 𝒢_weno on all computational cells.