CFTransport

godunov! = GodunovScheme(kind, dim, rank) <: OneDimFV{dim, rank, kind}
godunov!(backend, newtransported, transported, flux, mass)

Returns the callable godunov! that applies a one-dimensional Godunov scheme to dimension dim among rank. Computations are offloaded to backend.

If kind==:density then transported is a density. If kind==:scalar then transported is a scalar whose density is mass*transported.

If newtransported is the same as transported, it is updated in-place.

The Godunov scheme has 2, possibly 3 steps : 0 - (for densities) : compute scalar from density 1 - compute fluxes (upwind/downwind) 2 - update scalar or density

Boundary conditions are left to the user. If BCs must be enforced between steps 1-2 (e.g. zero boundary fluxes), the user should rather call individual steps concentrations!, fluxes! and FV_update! .

abstract type OneDimFV{kind,dim,rank} end

One-dimensional finite volume transport operator acting on dimension dim of arrays of rank rank. If kind==:density the operator transports a density field (e.g. in kg/m3) If kind==:scalar the operator transports a scalar field (e.g. in kg/kg)

abstract type OneDimOp{dim,rank} end

One-dimensional operator acting on dimension dim of arrays of rank rank.

abstract type Stencil <: Function end

A one-dimensional stencil is a function/closure/callable of the form

function stencil((i,j), params, arrays)
    a, b, ... = arrays
    a[i,j] = expression(params, b[i,j], b[i+1,j], b[i-1,j], ...)
end

vanleer! = VanLeerScheme(kind, limiter, dim, rank) <: OneDimFV{dim, rank, kind}
vanleer!(backend, newtransported, transported, flux, mass)

Returns the callable vanleer that applies a one-dimensional Van Leer scheme with limiter to dimension dim among rank. Computations are offloaded to backend.

If kind==:density then transported is a density. If kind==:scalar then transported is a scalar whose density is mass*transported.

If newtransported is the same as transported, it is updated in-place.

The VanLeer scheme has 3, possibly 4 steps : 0 - (for densities) : compute scalar from density 1 - compute slopes (with slope limiter) 2 - compute fluxes (upwind/downwind) 3 - update scalar or density

Boundary conditions are left to the user. If BCs must be enforced between steps 1-2 (e.g. zero boundary slopes) and 2-3 (e.g. zero boundary fluxes), the user should rather call individual steps concentrations!, slopes!, fluxes! and FV_update! .

stencil = expand_stencil(dim, rank, stencil1)

Return a stencil of rank rank that applies the one-dimensional stencil stencil1 to dimension dim. For example:

function stencil((i,j), params, arrays)
    a, b = arrays
    a[i,j] = expression(params, b[i,j], b[i+1,j], b[i-1,j])
end

stencil2 = expand_stencil(1, 2, stencil1)
for i in axes(a,1), j in axes(a,2)
    stencil2((i,j), params, (a,b))
end

is equivalent to:

for i in axes(a,1), j in axes(a,2)
    a[i,j] = expression(params, b[i,j], b[i,j+1], b[i,j-1])
end

See also Stencil

slope = minmod(slope1, slope2)

Minmod limiter

slope = minmod_simd(slope1, slope2)

Minmod limiter. This implementation avoids branching and may be more suitable for SIMD vectorization.

remap_fluxes!(mgr, mcoord::CFDomains.MassCoordinate, layout, flux, newmass, mass)

Computes the target (pseudo-)mass distribution newmass prescribed by mass coordinate mcoord and current mass distribution mass, as well as the vertical (pseudo-)mass flux flux needed to remap from current mass distribution mass to target newmass. layout specifies the data layout, see CFDoamins.VHLayout and CFDomains.HVLayout.

lim = ML()

Return a limiter that is is less diffusive than SL() but requires a method-of-lines time integration.

lim = SL()

Return a slope limiter which guarantees positivity with a semi-Lagrangian time integration.

lim = SL_simple()

Return a simplified version of SL(), slightly more diffusive and more efficient.

qflux, gradq, q = SLFV_flux!(qflux, gradq, q, mgr, lim, vsphere, disp, mass, mflux, qmass)

For each edge of vsphere, compute the time-integrated scalar flux qflux following the SLFV scheme, using the displacement disp computed by backwards_trajectories, the carrier mass mass and time-integrated mass flux mflux and the scalar mass qmass. Additionnally, for each primal cell, compute the slope-limited gradient gradq, and mixing ratio q

qflux, gradq and q may be ::Void, in which case they will be appropriately allocated. vsphere must be a VoronoiSphere or another struct or named tuple with the necessary fields. mgr is nothing or a LoopManager. In the latter case, mgr manages the computational loops.

mass, qmass and q are scalar fields known at primal cells. mflux and qflux are vector fields known by their contravariant components (integrals across primal cell edges). disp is a 3d-vector-valued field known by its values at edges. gradq is a 3d-vector-valued field known at primal cells.

If mflux is an AbstractVector representing a 2D field, disp and gradq are vectors of 3-uples. If mflux is an AbstractMatrix representing a 3D field, it must have layout VHLayout{1}, i.e. nz=size(mflux,1) is the number of layers. disp and gradq are then arrays of size (nz, 3, size(mflux,2)) and (nz, 3, size(qmass, 2)). This layout favors SIMD vectorization on CPUs and merged memory accesses on GPUs.

disp, dx = backwards_trajectories!(disp, dx, mgr, vsphere, mass, mflux)

For each edge of vsphere, compute the 3D displacement disp from the center of the upwind cell to the midpoint of the backwards trajectory starting from the edge midpoint. The backward trajectory is deduced from the mass and the time-integrated mass flux mflux. dx is the normal (to the edge) time-integrated velocity.

disp and dx may be ::Void, in which case they will be appropriately allocated. vsphere must be a VoronoiSphere or another struct or named tuple with the necessary fields. mgr is nothing or a LoopManager. In the latter case, mgr manages the computational loops.

mass is a scalar field known at primal cells. mflux is a vector field known by its contravariant components (integrals across primal cell edges). dx is a vector field known by its normal (to primal cell edges) components. disp is a 3d-vector-valued field known by its values at edges.

mass and mflux must be arrays of the same type. Udt is similar to mflux. If mflux is an AbstractVector representing a 2D field, disp is a similar vector, but of 3-uples. If mflux is an AbstractMatrix representing a 3D field, it must have layout VHLayout{1}, i.e. nz=size(mflux,1) is the number of layers. disp is then an array of size (nz, 3, size(mflux,2)). This layout favors SIMD vectorization on CPUs and merged memory accesses on GPUs.