erf_auxiliary Namespace Reference

Classes
class	AuxiliaryInertTracer
struct	MappedFaceFluxRateTag
struct	IntegratedMappedFaceFluxTag
class	MappedFaceField
struct	ProjectionRule
struct	ProjectionValidation
class	AuxiliaryProjection
struct	ConstTimedFieldView
struct	MutableTimedFieldView
struct	HostCarrierFluxView
struct	AuxiliaryStageRecipe
struct	AuxiliaryStageContext
class	CompletedStepFluxLedger
struct	ComponentDescriptor
class	AuxiliaryStateLayout

Typedefs
using	MappedFaceFluxRate = MappedFaceField< MappedFaceFluxRateTag >
using	IntegratedMappedFaceFlux = MappedFaceField< IntegratedMappedFaceFluxTag >

Enumerations
enum class	ProjectionTargetKind { BulkCoupling , Vapor , Number , AttachedMaterial }
enum class	ProjectionSourceKind { AuxiliaryMass , AuxiliaryNumber , AttachedMaterial }
enum class	HostIntegrator { CompressibleRK3 , AnelasticHeun , AnelasticMidPoint }
enum class	LimiterTrialBase { Anchor , Input }
enum class	AuxiliaryFieldValidationPolicy { Global , AssumeValid }

Functions
bool	SameCellLayout (const amrex::MultiFab &lhs, const amrex::MultiFab &rhs)
bool	BuildMappedCellMeasure (amrex::MultiFab &omega, const amrex::MultiFab &detJ, const amrex::MultiFab &mx, const amrex::MultiFab &my, std::string &diagnostic)
bool	BuildMappedDryAirCarrierFluxRate (MappedFaceFluxRate &rate, const amrex::MultiFab &rho_u, const amrex::MultiFab &rho_v, const amrex::MultiFab &omega, const amrex::MultiFab &ax, const amrex::MultiFab &ay, const amrex::MultiFab &az, const amrex::MultiFab &mf_uy, const amrex::MultiFab &mf_vx, const amrex::MultiFab &mx, const amrex::MultiFab &my, std::string &diagnostic)
bool	CopyNativeMappedDryAirCarrierFluxRate (MappedFaceFluxRate &rate, const amrex::MultiFab &rho_u, const amrex::MultiFab &rho_v, const amrex::MultiFab &rho_w, std::string &diagnostic)
bool	ComputeMaxMappedOutgoingRate (const MappedFaceFluxRate &rate, const amrex::MultiFab &measure, const amrex::MultiFab &conserved, const int density_component, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dx_inv, amrex::Real &max_rate, std::string &diagnostic)
bool	FixedDtExceedsMappedDonorLimit (const double fixed_dt, const amrex::Real max_outgoing_rate, double &hard_limit) noexcept
bool	ValidatePositiveFiniteComponent (const amrex::MultiFab &field, const int component, std::string &diagnostic)
bool	ValidateFiniteComponent (const amrex::MultiFab &field, const int component, std::string &diagnostic)
amrex::Real	MaxFaceFieldDifference (const MappedFaceFluxRate &lhs, const int lhs_comp, const MappedFaceFluxRate &rhs, const int rhs_comp)
void	AccumulateIntegratedFaceFlux (IntegratedMappedFaceFlux &ledger, const MappedFaceFluxRate &rate, const amrex::Real weight)
template<typename Tag >
bool	MappedFaceLayoutMatchesCellLayout (const MappedFaceField< Tag > &faces, const amrex::BoxArray &cell_ba, const amrex::DistributionMapping &dm)
template<typename Tag >
bool	MappedFaceLayoutMatchesCellLayout (const MappedFaceField< Tag > &faces, const amrex::MultiFab &cell)
template<typename LhsTag , typename RhsTag >
bool	SameMappedFaceLayout (const MappedFaceField< LhsTag > &lhs, const MappedFaceField< RhsTag > &rhs)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	ComputationalMappedDivergence (const amrex::Real fx_hi, const amrex::Real fx_lo, const amrex::Real fy_hi, const amrex::Real fy_lo, const amrex::Real fz_hi, const amrex::Real fz_lo, const amrex::Real dx_inv, const amrex::Real dy_inv, const amrex::Real dzeta_inv) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	MappedOutgoingDemandRate (const amrex::Real fx_hi, const amrex::Real fx_lo, const amrex::Real fy_hi, const amrex::Real fy_lo, const amrex::Real fz_hi, const amrex::Real fz_lo, const amrex::Real measure, const amrex::Real density, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dx_inv) noexcept
template<typename Tag >
void	ApplyMappedFluxTendency (const MappedFaceField< Tag > &flux, const int flux_comp, const amrex::MultiFab &omega, const int omega_comp, amrex::MultiFab &tendency, const int tendency_comp, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dx_inv)
const char *	HostIntegratorName (const HostIntegrator method) noexcept
bool	MakeAuxiliaryStageRecipe (const HostIntegrator method, const int stage, const double host_stage_interval, AuxiliaryStageRecipe &recipe, std::string &diagnostic)
bool	BuildAuxiliaryIntensiveState (const ConstTimedFieldView &state_input, const ConstTimedFieldView &rho_input, const double input_time, amrex::MultiFab &intensive, const AuxiliaryFieldValidationPolicy validation, std::string &diagnostic)
bool	AuxiliaryStageTargetIsDisjoint (const AuxiliaryStageContext &context)
void	ApplyAuxiliaryMappedStage (const AuxiliaryStageContext &context, const MappedFaceFluxRate &rate, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &dx_inv, const int rate_component)

Typedef Documentation

IntegratedMappedFaceFlux

using erf_auxiliary::IntegratedMappedFaceFlux = typedef MappedFaceField<IntegratedMappedFaceFluxTag>

MappedFaceFluxRate

using erf_auxiliary::MappedFaceFluxRate = typedef MappedFaceField<MappedFaceFluxRateTag>

Enumeration Type Documentation

AuxiliaryFieldValidationPolicy

enum erf_auxiliary::AuxiliaryFieldValidationPolicy

strong

Enumerator
Global	Perform global positive-density and finite-result validation.
AssumeValid	Skip global checks. The caller guarantees every input state value is finite, density is finite and strictly positive, and each state/density quotient is finite and representable in amrex::Real. No global reduction is performed by this policy.

                                          {
    /** Perform global positive-density and finite-result validation. */
    Global,
    /**
     * Skip global checks. The caller guarantees every input state value is
     * finite, density is finite and strictly positive, and each state/density
     * quotient is finite and representable in amrex::Real. No global
     * reduction is performed by this policy.
     */
    AssumeValid
};

HostIntegrator

enum erf_auxiliary::HostIntegrator

strong

Enumerator
CompressibleRK3
AnelasticHeun
AnelasticMidPoint

                          {
    CompressibleRK3,
    AnelasticHeun,
    AnelasticMidPoint
};

LimiterTrialBase

enum erf_auxiliary::LimiterTrialBase

strong

State view that bounds the low-order limiter trial for this host stage.

Enumerator
Anchor
Input

                            {
    Anchor,
    Input
};

ProjectionSourceKind

enum erf_auxiliary::ProjectionSourceKind

strong

Enumerator
AuxiliaryMass
AuxiliaryNumber
AttachedMaterial

{ AuxiliaryMass, AuxiliaryNumber, AttachedMaterial };

ProjectionTargetKind

enum erf_auxiliary::ProjectionTargetKind

strong

Enumerator
BulkCoupling
Vapor
Number
AttachedMaterial

{ BulkCoupling, Vapor, Number, AttachedMaterial };

Function Documentation

AccumulateIntegratedFaceFlux()

void erf_auxiliary::AccumulateIntegratedFaceFlux	(	IntegratedMappedFaceFlux &	ledger,
		const MappedFaceFluxRate &	rate,
		amrex::Real	weight
	)

Scale every component and direction of a flux rate into an integrated ledger.

{
    AMREX_ALWAYS_ASSERT(ledger.is_defined() && rate.is_defined());
    AMREX_ALWAYS_ASSERT(ledger.nComp() == rate.nComp());
    AMREX_ALWAYS_ASSERT(SameMappedFaceLayout(ledger, rate));
    AMREX_ALWAYS_ASSERT(std::isfinite(weight) && weight >= amrex::Real(0.0));
    for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
        amrex::MultiFab::Saxpy(ledger.dir(dir), weight, rate.dir(dir),
                               0, 0, ledger.nComp(), 0);
    }
}

Here is the call graph for this function:

ApplyAuxiliaryMappedStage()

void erf_auxiliary::ApplyAuxiliaryMappedStage	(	const AuxiliaryStageContext &	context,
		const MappedFaceFluxRate &	rate,
		const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &	dx_inv,
		int	rate_component = 0
	)

Apply H_target = a H_anchor + b H_input - c D_xi(F_rate).

{
    AMREX_ALWAYS_ASSERT(context.level >= 0 && context.stage >= 0);
    AMREX_ALWAYS_ASSERT(std::isfinite(context.step_old_time));
    AMREX_ALWAYS_ASSERT(std::isfinite(context.input_time));
    AMREX_ALWAYS_ASSERT(std::isfinite(context.target_time));
    AMREX_ALWAYS_ASSERT(valid_const_view(context.state_anchor, context.step_old_time));
    AMREX_ALWAYS_ASSERT(valid_const_view(context.state_input, context.input_time));
    AMREX_ALWAYS_ASSERT(context.state_target.field != nullptr);
    AMREX_ALWAYS_ASSERT(context.state_target.component >= 0 &&
                        context.state_target.component < context.state_target.field->nComp());
    AMREX_ALWAYS_ASSERT(std::isfinite(context.state_target.time) &&
                        context.state_target.time == context.target_time);
    AMREX_ALWAYS_ASSERT(valid_const_view(context.rho_anchor, context.step_old_time));
    AMREX_ALWAYS_ASSERT(valid_const_view(context.rho_input, context.input_time));
    AMREX_ALWAYS_ASSERT(valid_const_view(context.rho_target, context.target_time));
    AMREX_ALWAYS_ASSERT(valid_const_view(context.measure_anchor, context.step_old_time));
    AMREX_ALWAYS_ASSERT(valid_const_view(context.measure_input, context.input_time));
    AMREX_ALWAYS_ASSERT(valid_const_view(context.measure_target, context.target_time));
    AMREX_ALWAYS_ASSERT(context.carrier.x != nullptr && context.carrier.y != nullptr &&
                        context.carrier.z != nullptr);
    AMREX_ALWAYS_ASSERT(rate.is_defined() && rate_component >= 0 &&
                        rate_component < rate.nComp());
    AMREX_ALWAYS_ASSERT(AuxiliaryStageTargetIsDisjoint(context));
 
    const auto& target = *context.state_target.field;
    const auto& anchor = *context.state_anchor.field;
    const auto& input = *context.state_input.field;
    const auto& omega_anchor = *context.measure_anchor.field;
    const auto& omega_input = *context.measure_input.field;
    const auto& omega_target = *context.measure_target.field;
    AMREX_ALWAYS_ASSERT(SameCellLayout(target, anchor));
    AMREX_ALWAYS_ASSERT(SameCellLayout(target, input));
    AMREX_ALWAYS_ASSERT(SameCellLayout(target, *context.rho_anchor.field));
    AMREX_ALWAYS_ASSERT(SameCellLayout(target, *context.rho_input.field));
    AMREX_ALWAYS_ASSERT(SameCellLayout(target, *context.rho_target.field));
    AMREX_ALWAYS_ASSERT(SameCellLayout(target, omega_anchor));
    AMREX_ALWAYS_ASSERT(SameCellLayout(target, omega_input));
    AMREX_ALWAYS_ASSERT(SameCellLayout(target, omega_target));
    AMREX_ALWAYS_ASSERT(MappedFaceLayoutMatchesCellLayout(rate, target));
 
    const amrex::Real anchor_weight = static_cast<amrex::Real>(context.recurrence.anchor_weight);
    const amrex::Real input_weight = static_cast<amrex::Real>(context.recurrence.input_weight);
    const amrex::Real face_weight =
        static_cast<amrex::Real>(context.recurrence.face_rate_time_coefficient);
    const amrex::Real dx = dx_inv[0];
    const amrex::Real dy = dx_inv[1];
    const amrex::Real dz = dx_inv[2];
 
    for (amrex::MFIter mfi(target, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
        const amrex::Box bx = mfi.tilebox();
        const auto ua = anchor.const_array(mfi);
        const auto ui = input.const_array(mfi);
        const auto ut = context.state_target.field->array(mfi);
        const auto oa = omega_anchor.const_array(mfi);
        const auto oi = omega_input.const_array(mfi);
        const auto omega_target_array = omega_target.const_array(mfi);
        const auto fx = rate.dir(0).const_array(mfi);
        const auto fy = rate.dir(1).const_array(mfi);
        const auto fz = rate.dir(2).const_array(mfi);
        const int ua_comp = context.state_anchor.component;
        const int ui_comp = context.state_input.component;
        const int ut_comp = context.state_target.component;
        const int oa_comp = context.measure_anchor.component;
        const int oi_comp = context.measure_input.component;
        const int omega_target_component = context.measure_target.component;
        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
            const amrex::Real div = ComputationalMappedDivergence(
                fx(i + 1, j, k, rate_component), fx(i, j, k, rate_component),
                fy(i, j + 1, k, rate_component), fy(i, j, k, rate_component),
                fz(i, j, k + 1, rate_component), fz(i, j, k, rate_component),
                dx, dy, dz);
            const amrex::Real h_target =
                anchor_weight * oa(i, j, k, oa_comp) * ua(i, j, k, ua_comp) +
                input_weight * oi(i, j, k, oi_comp) * ui(i, j, k, ui_comp) -
                face_weight * div;
            ut(i, j, k, ut_comp) = h_target /
                omega_target_array(i, j, k, omega_target_component);
        });
    }
}

Referenced by erf_sbm::SBMTransport::advance_stage(), and erf_auxiliary::AuxiliaryInertTracer::advance_stage().

Here is the call graph for this function:

Here is the caller graph for this function:

ApplyMappedFluxTendency()

template<typename Tag >

void erf_auxiliary::ApplyMappedFluxTendency	(	const MappedFaceField< Tag > &	flux,
		const int	flux_comp,
		const amrex::MultiFab &	omega,
		const int	omega_comp,
		amrex::MultiFab &	tendency,
		const int	tendency_comp,
		const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &	dx_inv
	)

Apply -D_xi(flux)/omega into one selected output component. The supplied face values are already in ERF's mapped coordinate convention.

{
    AMREX_ALWAYS_ASSERT(flux.is_defined());
    AMREX_ALWAYS_ASSERT(flux_comp >= 0 && flux_comp < flux.nComp());
    AMREX_ALWAYS_ASSERT(omega_comp >= 0 && omega_comp < omega.nComp());
    AMREX_ALWAYS_ASSERT(tendency_comp >= 0 && tendency_comp < tendency.nComp());
    AMREX_ALWAYS_ASSERT(SameCellLayout(omega, tendency));
    AMREX_ALWAYS_ASSERT(MappedFaceLayoutMatchesCellLayout(flux, tendency));
 
    for (amrex::MFIter mfi(tendency, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
        const amrex::Box bx = mfi.tilebox();
        const auto fx = flux.dir(0).const_array(mfi);
        const auto fy = flux.dir(1).const_array(mfi);
        const auto fz = flux.dir(2).const_array(mfi);
        const auto measure = omega.const_array(mfi);
        const auto out = tendency.array(mfi);
        const auto dx = dx_inv[0];
        const auto dy = dx_inv[1];
        const auto dz = dx_inv[2];
 
        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
            const amrex::Real div = ComputationalMappedDivergence(
                fx(i + 1, j, k, flux_comp), fx(i, j, k, flux_comp),
                fy(i, j + 1, k, flux_comp), fy(i, j, k, flux_comp),
                fz(i, j, k + 1, flux_comp), fz(i, j, k, flux_comp),
                dx, dy, dz);
            out(i, j, k, tendency_comp) = -div / measure(i, j, k, omega_comp);
        });
    }
}

Here is the call graph for this function:

AuxiliaryStageTargetIsDisjoint()

bool erf_auxiliary::AuxiliaryStageTargetIsDisjoint

(

const AuxiliaryStageContext &

context

)

Whether the mutable stage target is a distinct MultiFab from both inputs.

{
    return context.state_target.field != nullptr &&
           (context.state_anchor.field == nullptr ||
            !storage_overlaps(*context.state_target.field, *context.state_anchor.field)) &&
           (context.state_input.field == nullptr ||
            !storage_overlaps(*context.state_target.field, *context.state_input.field));
}

Referenced by ApplyAuxiliaryMappedStage().

Here is the caller graph for this function:

BuildAuxiliaryIntensiveState()

bool erf_auxiliary::BuildAuxiliaryIntensiveState	(	const ConstTimedFieldView &	state_input,
		const ConstTimedFieldView &	rho_input,
		double	input_time,
		amrex::MultiFab &	intensive,
		AuxiliaryFieldValidationPolicy	validation,
		std::string &	diagnostic
	)

Form z=U_input/rho_input from the explicit input-time views.

{
    diagnostic.clear();
    if (!std::isfinite(input_time) ||
        !valid_const_view(state_input, input_time) ||
        !valid_const_view(rho_input, input_time)) {
        diagnostic = "auxiliary intensive state requires state and density views at input_time";
        return false;
    }
    if (!SameCellLayout(intensive, *state_input.field) ||
        !SameCellLayout(intensive, *rho_input.field) ||
        intensive.nComp() < 1) {
        diagnostic = "auxiliary intensive state fields must share a cell layout";
        return false;
    }
    if (validation == AuxiliaryFieldValidationPolicy::Global) {
        if (!ValidatePositiveFiniteComponent(*rho_input.field, rho_input.component, diagnostic)) {
            return false;
        }
    } else if (validation != AuxiliaryFieldValidationPolicy::AssumeValid) {
        diagnostic = "auxiliary intensive state received an unknown validation policy";
        return false;
    }
    for (amrex::MFIter mfi(intensive, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
        const amrex::Box bx = mfi.tilebox();
        const auto state = state_input.field->const_array(mfi);
        const auto rho = rho_input.field->const_array(mfi);
        const auto out = intensive.array(mfi);
        const int state_comp = state_input.component;
        const int rho_comp = rho_input.component;
        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
            out(i, j, k, 0) = state(i, j, k, state_comp) / rho(i, j, k, rho_comp);
        });
    }
    if (validation == AuxiliaryFieldValidationPolicy::Global) {
        if (!ValidateFiniteComponent(intensive, 0, diagnostic)) {
            return false;
        }
    }
    return true;
}

Referenced by erf_auxiliary::AuxiliaryInertTracer::advance_stage().

Here is the call graph for this function:

Here is the caller graph for this function:

BuildMappedCellMeasure()

bool erf_auxiliary::BuildMappedCellMeasure	(	amrex::MultiFab &	omega,
		const amrex::MultiFab &	detJ,
		const amrex::MultiFab &	mx,
		const amrex::MultiFab &	my,
		std::string &	diagnostic
	)

Build omega=detJ/(mx*my), returning false and a diagnostic on invalid cells.

{
    diagnostic.clear();
    if (omega.nComp() != 1 || detJ.nComp() < 1 || mx.nComp() < 1 || my.nComp() < 1 ||
        !SameCellLayout(omega, detJ) || !same_horizontal_layout(omega, mx) ||
        !same_horizontal_layout(omega, my) ||
        mx.boxArray() != my.boxArray() || mx.DistributionMap() != my.DistributionMap()) {
        diagnostic = "mapped measure requires a cell-centered Jacobian, matching horizontal map-factor layouts, and one output component";
        return false;
    }
 
    amrex::MultiFab invalid(omega.boxArray(), omega.DistributionMap(), 1, 0);
    invalid.setVal(amrex::Real(0.0));
    for (amrex::MFIter mfi(omega, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
        const amrex::Box bx = mfi.tilebox();
        const auto det = detJ.const_array(mfi);
        const auto map_x = mx.const_array(mfi);
        const auto map_y = my.const_array(mfi);
        const auto out = omega.array(mfi);
        const auto bad = invalid.array(mfi);
        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
            const amrex::Real d = det(i, j, k, 0);
            // ERF's horizontal map factors are vertically invariant and all
            // native mapped divergence paths read them from k=0.
            const amrex::Real x = map_x(i, j, 0, 0);
            const amrex::Real y = map_y(i, j, 0, 0);
            const amrex::Real measure = d / (x * y);
            const bool valid = amrex::Math::isfinite(d) && d > amrex::Real(0.0) &&
                               amrex::Math::isfinite(x) && x > amrex::Real(0.0) &&
                               amrex::Math::isfinite(y) && y > amrex::Real(0.0) &&
                               amrex::Math::isfinite(measure) && measure > amrex::Real(0.0);
            out(i, j, k, 0) = valid ? measure : amrex::Real(0.0);
            bad(i, j, k, 0) = valid ? amrex::Real(0.0) : amrex::Real(1.0);
        });
    }
 
    const amrex::Real invalid_count = invalid.sum(0);
    if (invalid_count != amrex::Real(0.0)) {
        diagnostic = "detJ, horizontal map factors, and mapped measure must be finite and positive";
        return false;
    }
    return true;
}

Referenced by erf_auxiliary::AuxiliaryInertTracer::rebuild_static_measure(), and erf_sbm::SBMTransport::rebuild_static_measure().

Here is the call graph for this function:

Here is the caller graph for this function:

BuildMappedDryAirCarrierFluxRate()

bool erf_auxiliary::BuildMappedDryAirCarrierFluxRate	(	MappedFaceFluxRate &	rate,
		const amrex::MultiFab &	rho_u,
		const amrex::MultiFab &	rho_v,
		const amrex::MultiFab &	omega,
		const amrex::MultiFab &	ax,
		const amrex::MultiFab &	ay,
		const amrex::MultiFab &	az,
		const amrex::MultiFab &	mf_uy,
		const amrex::MultiFab &	mf_vx,
		const amrex::MultiFab &	mx,
		const amrex::MultiFab &	my,
		std::string &	diagnostic
	)

Build the production mapped dry-air face rates from face momentum and metrics.

{
    diagnostic.clear();
    if (!rate.is_defined() || rate.nComp() != 1 || rho_u.nComp() < 1 ||
        rho_v.nComp() < 1 || omega.nComp() < 1 || ax.nComp() < 1 ||
        ay.nComp() < 1 || az.nComp() < 1 || mf_uy.nComp() < 1 ||
        mf_vx.nComp() < 1 || mx.nComp() < 1 || my.nComp() < 1 ||
        rho_u.boxArray() != rate.dir(0).boxArray() ||
        rho_u.DistributionMap() != rate.dir(0).DistributionMap() ||
        ax.boxArray() != rate.dir(0).boxArray() ||
        ax.DistributionMap() != rate.dir(0).DistributionMap() ||
        rho_v.boxArray() != rate.dir(1).boxArray() ||
        rho_v.DistributionMap() != rate.dir(1).DistributionMap() ||
        ay.boxArray() != rate.dir(1).boxArray() ||
        ay.DistributionMap() != rate.dir(1).DistributionMap() ||
        omega.boxArray() != rate.dir(2).boxArray() ||
        omega.DistributionMap() != rate.dir(2).DistributionMap() ||
        az.boxArray() != rate.dir(2).boxArray() ||
        az.DistributionMap() != rate.dir(2).DistributionMap() ||
        !same_horizontal_layout(rate.dir(0), mf_uy) ||
        !same_horizontal_layout(rate.dir(1), mf_vx) ||
        mx.boxArray() != my.boxArray() ||
        mx.DistributionMap() != my.DistributionMap() ||
        !same_horizontal_layout(rate.dir(2), mx) ||
        !same_horizontal_layout(rate.dir(2), my)) {
        diagnostic = "mapped dry-air carrier rates require matching staggered momentum, area, map-factor, and output layouts";
        return false;
    }
 
    for (amrex::MFIter mfi(rate.dir(0), amrex::TilingIfNotGPU());
         mfi.isValid(); ++mfi) {
        const amrex::Box bx = mfi.tilebox();
        const auto momentum = rho_u.const_array(mfi);
        const auto area = ax.const_array(mfi);
        const auto map = mf_uy.const_array(mfi);
        const auto out = rate.dir(0).array(mfi);
        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
            out(i, j, k, 0) = area(i, j, k, 0) * momentum(i, j, k, 0) /
                              map(i, j, 0, 0);
        });
    }
    for (amrex::MFIter mfi(rate.dir(1), amrex::TilingIfNotGPU());
         mfi.isValid(); ++mfi) {
        const amrex::Box bx = mfi.tilebox();
        const auto momentum = rho_v.const_array(mfi);
        const auto area = ay.const_array(mfi);
        const auto map = mf_vx.const_array(mfi);
        const auto out = rate.dir(1).array(mfi);
        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
            out(i, j, k, 0) = area(i, j, k, 0) * momentum(i, j, k, 0) /
                              map(i, j, 0, 0);
        });
    }
    for (amrex::MFIter mfi(rate.dir(2), amrex::TilingIfNotGPU());
         mfi.isValid(); ++mfi) {
        const amrex::Box bx = mfi.tilebox();
        const auto contravariant_momentum = omega.const_array(mfi);
        const auto area = az.const_array(mfi);
        const auto map_x = mx.const_array(mfi);
        const auto map_y = my.const_array(mfi);
        const auto out = rate.dir(2).array(mfi);
        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
            out(i, j, k, 0) = area(i, j, k, 0) *
                              contravariant_momentum(i, j, k, 0) /
                              (map_x(i, j, 0, 0) * map_y(i, j, 0, 0));
        });
    }
 
    if (!rate.dir(0).is_finite(0, 1, 0) ||
        !rate.dir(1).is_finite(0, 1, 0) ||
        !rate.dir(2).is_finite(0, 1, 0)) {
        diagnostic = "mapped dry-air carrier face rates are nonfinite";
        return false;
    }
    return true;
}

Referenced by ERF::estTimeStep().

Here is the call graph for this function:

Here is the caller graph for this function:

ComputationalMappedDivergence()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real erf_auxiliary::ComputationalMappedDivergence ( const amrex::Real  fx_hi, const amrex::Real  fx_lo, const amrex::Real  fy_hi, const amrex::Real  fy_lo, const amrex::Real  fz_hi, const amrex::Real  fz_lo, const amrex::Real  dx_inv, const amrex::Real  dy_inv, const amrex::Real  dzeta_inv  )

noexcept

Computational-coordinate divergence only: no metric or physical face area.

{
    return (fx_hi - fx_lo) * dx_inv +
           (fy_hi - fy_lo) * dy_inv +
           (fz_hi - fz_lo) * dzeta_inv;
}

Referenced by erf_sbm::SBMTransport::advance_stage(), erf_auxiliary::AuxiliaryInertTracer::advance_stage(), ApplyAuxiliaryMappedStage(), and ApplyMappedFluxTendency().

Here is the caller graph for this function:

ComputeMaxMappedOutgoingRate()

bool erf_auxiliary::ComputeMaxMappedOutgoingRate	(	const MappedFaceFluxRate &	rate,
		const amrex::MultiFab &	measure,
		const amrex::MultiFab &	conserved,
		int	density_component,
		const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &	dx_inv,
		amrex::Real &	max_rate,
		std::string &	diagnostic
	)

Reduce the largest mapped donor outgoing rate over a cell-centered domain.

{
    diagnostic.clear();
    max_rate = amrex::Real(0.0);
    if (!rate.is_defined() || rate.nComp() != 1 ||
        !MappedFaceLayoutMatchesCellLayout(rate, measure) ||
        !SameCellLayout(measure, conserved) || measure.nComp() < 1) {
        diagnostic = "mapped donor-rate reduction received incompatible cell or face layouts";
        return false;
    }
    if (!ValidatePositiveFiniteComponent(measure, 0, diagnostic)) return false;
    if (!ValidatePositiveFiniteComponent(conserved, density_component,
                                         diagnostic)) return false;
 
    amrex::MultiFab cell_rate(measure.boxArray(), measure.DistributionMap(), 1,
                              0);
    for (amrex::MFIter mfi(cell_rate, amrex::TilingIfNotGPU()); mfi.isValid();
         ++mfi) {
        const amrex::Box bx = mfi.tilebox();
        const auto fx = rate.dir(0).const_array(mfi);
        const auto fy = rate.dir(1).const_array(mfi);
        const auto fz = rate.dir(2).const_array(mfi);
        const auto omega = measure.const_array(mfi);
        const auto rho = conserved.const_array(mfi);
        const auto out = cell_rate.array(mfi);
        const int rho_comp = density_component;
        const auto inv = dx_inv;
        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
            out(i, j, k, 0) = MappedOutgoingDemandRate(
                fx(i + 1, j, k, 0), fx(i, j, k, 0),
                fy(i, j + 1, k, 0), fy(i, j, k, 0),
                fz(i, j, k + 1, 0), fz(i, j, k, 0),
                omega(i, j, k, 0), rho(i, j, k, rho_comp), inv);
        });
    }
    max_rate = cell_rate.max(0);
    if (!std::isfinite(max_rate) || max_rate < amrex::Real(0.0)) {
        diagnostic = "mapped donor outgoing-rate maximum is nonfinite or negative";
        return false;
    }
    return true;
}

Referenced by ERF::estTimeStep().

Here is the call graph for this function:

Here is the caller graph for this function:

CopyNativeMappedDryAirCarrierFluxRate()

bool erf_auxiliary::CopyNativeMappedDryAirCarrierFluxRate	(	MappedFaceFluxRate &	rate,
		const amrex::MultiFab &	rho_u,
		const amrex::MultiFab &	rho_v,
		const amrex::MultiFab &	rho_w,
		std::string &	diagnostic
	)

Copy face fields that are already in ERF's native mapped carrier convention.

{
    diagnostic.clear();
    if (AMREX_SPACEDIM != 3 || !rate.is_defined() || rate.nComp() != 1 ||
        rho_u.nComp() < 1 || rho_v.nComp() < 1 || rho_w.nComp() < 1) {
        diagnostic = "native mapped dry-air carrier requires three dimensions, "
                     "one output component, and one component per input face";
        return false;
    }
 
    const amrex::MultiFab* input_faces[3] = {&rho_u, &rho_v, &rho_w};
    for (int dir = 0; dir < 3; ++dir) {
        const auto& output = rate.dir(dir);
        const auto& input = *input_faces[dir];
        if (output.boxArray() != input.boxArray() ||
            output.DistributionMap() != input.DistributionMap()) {
            diagnostic = "native mapped dry-air carrier input must match the "
                         "output's x/y/z face staggering and distribution map";
            return false;
        }
    }
 
    for (int dir = 0; dir < 3; ++dir) {
        amrex::MultiFab::Copy(rate.dir(dir), *input_faces[dir], 0, 0, 1, 0);
    }
    if (!rate.dir(0).is_finite(0, 1, 0) ||
        !rate.dir(1).is_finite(0, 1, 0) ||
        !rate.dir(2).is_finite(0, 1, 0)) {
        diagnostic = "native mapped dry-air carrier face values are nonfinite";
        return false;
    }
    return true;
}

Referenced by ERF::estTimeStep().

Here is the call graph for this function:

Here is the caller graph for this function:

FixedDtExceedsMappedDonorLimit()

bool erf_auxiliary::FixedDtExceedsMappedDonorLimit ( double  fixed_dt, amrex::Real  max_outgoing_rate, double &  hard_limit  )

noexcept

Test a fixed interval against the hard mapped donor-positivity bound.

{
    hard_limit = std::numeric_limits<double>::infinity();
    if (std::isfinite(max_outgoing_rate) &&
        max_outgoing_rate > amrex::Real(0.0)) {
        hard_limit = 1.0 / static_cast<double>(max_outgoing_rate);
    }
    if (max_outgoing_rate == amrex::Real(0.0) || fixed_dt <= 0.0)
        return false;
    if (!std::isfinite(fixed_dt) ||
        !std::isfinite(max_outgoing_rate) ||
        max_outgoing_rate < amrex::Real(0.0)) {
        hard_limit = 0.0;
        return true;
    }
 
    if (!std::isfinite(hard_limit)) return false;
    const double scale = std::max(std::abs(fixed_dt), std::abs(hard_limit));
    const double tolerance = 64.0 *
        static_cast<double>(std::numeric_limits<amrex::Real>::epsilon()) *
        scale;
    return fixed_dt - hard_limit > tolerance;
}

Referenced by ERF::estTimeStep().

Here is the caller graph for this function:

HostIntegratorName()

const char * erf_auxiliary::HostIntegratorName ( const HostIntegrator  method )

noexcept

{
    switch (method) {
    case HostIntegrator::CompressibleRK3: return "CompressibleRK3";
    case HostIntegrator::AnelasticHeun: return "AnelasticHeun";
    case HostIntegrator::AnelasticMidPoint: return "AnelasticMidPoint";
    }
    return "Unknown";
}

Referenced by erf_auxiliary::AuxiliaryInertTracer::advance_stage().

Here is the caller graph for this function:

MakeAuxiliaryStageRecipe()

bool erf_auxiliary::MakeAuxiliaryStageRecipe	(	HostIntegrator	method,
		int	stage,
		double	host_stage_interval,
		AuxiliaryStageRecipe &	recipe,
		std::string &	diagnostic
	)

Return the explicit ERF recurrence recipe for one supported host stage.

{
    diagnostic.clear();
    recipe = {};
    if (!std::isfinite(host_stage_interval) || !(host_stage_interval > 0.0)) {
        diagnostic = "host stage interval must be finite and positive";
        return false;
    }
    if (method == HostIntegrator::AnelasticMidPoint) {
        diagnostic = "AnelasticMidPoint is represented but is not qualified by the M2 fixture";
        return false;
    }
    if (stage < 0 || stage >= stage_count(method)) {
        diagnostic = "stage index is outside the selected host integrator sequence";
        return false;
    }
 
    if (method == HostIntegrator::CompressibleRK3) {
        recipe.limiter_trial_base = LimiterTrialBase::Anchor;
        recipe.anchor_weight = 1.0;
        recipe.limiter_trial_interval = host_stage_interval;
        recipe.face_rate_time_coefficient = host_stage_interval;
        if (stage == 2) { recipe.completed_ledger_time = host_stage_interval; }
    } else if (method == HostIntegrator::AnelasticHeun) {
        if (stage == 0) {
            recipe.limiter_trial_base = LimiterTrialBase::Anchor;
            recipe.anchor_weight = 1.0;
            recipe.limiter_trial_interval = host_stage_interval;
            recipe.face_rate_time_coefficient = host_stage_interval;
            recipe.completed_ledger_time = 0.5 * host_stage_interval;
        } else {
            recipe.limiter_trial_base = LimiterTrialBase::Input;
            recipe.anchor_weight = 0.5;
            recipe.input_weight = 0.5;
            recipe.limiter_trial_interval = host_stage_interval;
            recipe.face_rate_time_coefficient = 0.5 * host_stage_interval;
            recipe.completed_ledger_time = 0.5 * host_stage_interval;
        }
    }
    return true;
}

Referenced by erf_sbm::SBMTransport::advance_stage(), and erf_auxiliary::AuxiliaryInertTracer::advance_stage().

Here is the caller graph for this function:

MappedFaceLayoutMatchesCellLayout() [1/2]

template<typename Tag >

bool erf_auxiliary::MappedFaceLayoutMatchesCellLayout	(	const MappedFaceField< Tag > &	faces,
		const amrex::BoxArray &	cell_ba,
		const amrex::DistributionMapping &	dm
	)

Check directional face staggering and ownership against a cell layout.

{
    if (!faces.is_defined()) { return false; }
    for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
        const auto& face = faces.dir(dir);
        if (face.boxArray() !=
                amrex::convert(cell_ba, amrex::IntVect::TheDimensionVector(dir)) ||
            face.DistributionMap() != dm) {
            return false;
        }
    }
    return true;
}

Referenced by ApplyAuxiliaryMappedStage(), ApplyMappedFluxTendency(), ComputeMaxMappedOutgoingRate(), and MappedFaceLayoutMatchesCellLayout().

Here is the call graph for this function:

Here is the caller graph for this function:

MappedFaceLayoutMatchesCellLayout() [2/2]

template<typename Tag >

bool erf_auxiliary::MappedFaceLayoutMatchesCellLayout	(	const MappedFaceField< Tag > &	faces,
		const amrex::MultiFab &	cell
	)

{
    return MappedFaceLayoutMatchesCellLayout(faces, cell.boxArray(), cell.DistributionMap());
}

Here is the call graph for this function:

MappedOutgoingDemandRate()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real erf_auxiliary::MappedOutgoingDemandRate ( const amrex::Real  fx_hi, const amrex::Real  fx_lo, const amrex::Real  fy_hi, const amrex::Real  fy_lo, const amrex::Real  fz_hi, const amrex::Real  fz_lo, const amrex::Real  measure, const amrex::Real  density, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > &  dx_inv  )

noexcept

Outgoing donor rate for one cell in mapped coordinates.

{
    const amrex::Real outward =
        (amrex::max(fx_hi, amrex::Real(0.0)) +
         amrex::max(-fx_lo, amrex::Real(0.0))) * dx_inv[0] +
        (amrex::max(fy_hi, amrex::Real(0.0)) +
         amrex::max(-fy_lo, amrex::Real(0.0))) * dx_inv[1] +
        (amrex::max(fz_hi, amrex::Real(0.0)) +
         amrex::max(-fz_lo, amrex::Real(0.0))) * dx_inv[2];
    return outward / (measure * density);
}

Referenced by erf_sbm::SBMTransport::advance_stage(), and ComputeMaxMappedOutgoingRate().

Here is the caller graph for this function:

MaxFaceFieldDifference()

amrex::Real erf_auxiliary::MaxFaceFieldDifference	(	const MappedFaceFluxRate &	lhs,
		int	lhs_comp,
		const MappedFaceFluxRate &	rhs,
		int	rhs_comp
	)

Return the max absolute difference between two component-selected face fields.

{
    AMREX_ALWAYS_ASSERT(lhs.is_defined() && rhs.is_defined());
    AMREX_ALWAYS_ASSERT(lhs_comp >= 0 && lhs_comp < lhs.nComp());
    AMREX_ALWAYS_ASSERT(rhs_comp >= 0 && rhs_comp < rhs.nComp());
    AMREX_ALWAYS_ASSERT(SameMappedFaceLayout(lhs, rhs));
    amrex::Real maximum = amrex::Real(0.0);
    for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
        const auto& a = lhs.dir(dir);
        const auto& b = rhs.dir(dir);
        amrex::MultiFab difference(a.boxArray(), a.DistributionMap(), 1, 0);
        for (amrex::MFIter mfi(difference, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
            const amrex::Box bx = mfi.tilebox();
            const auto av = a.const_array(mfi);
            const auto bv = b.const_array(mfi);
            const auto out = difference.array(mfi);
            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
                out(i, j, k, 0) = amrex::Math::abs(av(i, j, k, lhs_comp) -
                                                    bv(i, j, k, rhs_comp));
            });
        }
        maximum = amrex::max(maximum, difference.norm0(0));
    }
    return maximum;
}

Referenced by erf_auxiliary::AuxiliaryInertTracer::advance_stage().

Here is the call graph for this function:

Here is the caller graph for this function:

SameCellLayout()

bool erf_auxiliary::SameCellLayout	(	const amrex::MultiFab &	lhs,
		const amrex::MultiFab &	rhs
	)

Exact cell-centered BoxArray and DistributionMapping identity.

{
    return lhs.boxArray() == rhs.boxArray() &&
           lhs.DistributionMap() == rhs.DistributionMap();
}

Referenced by erf_sbm::SBMTransport::advance_stage(), ApplyAuxiliaryMappedStage(), ApplyMappedFluxTendency(), BuildAuxiliaryIntensiveState(), BuildMappedCellMeasure(), ComputeMaxMappedOutgoingRate(), erf_sbm::ProlongCarrierRelativeSpectrum(), and erf_sbm::RestrictMappedSpectrum().

Here is the caller graph for this function:

SameMappedFaceLayout()

template<typename LhsTag , typename RhsTag >

bool erf_auxiliary::SameMappedFaceLayout	(	const MappedFaceField< LhsTag > &	lhs,
		const MappedFaceField< RhsTag > &	rhs
	)

Check that two semantically tagged face fields have identical face layouts.

{
    if (!lhs.is_defined() || !rhs.is_defined()) { return false; }
    for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
        if (lhs.dir(dir).boxArray() != rhs.dir(dir).boxArray() ||
            lhs.dir(dir).DistributionMap() != rhs.dir(dir).DistributionMap()) {
            return false;
        }
    }
    return true;
}

Referenced by erf_auxiliary::CompletedStepFluxLedger::accept_stage(), erf_auxiliary::CompletedStepFluxLedger::accumulate_stage_component(), AccumulateIntegratedFaceFlux(), and MaxFaceFieldDifference().

Here is the call graph for this function:

Here is the caller graph for this function:

ValidateFiniteComponent()

bool erf_auxiliary::ValidateFiniteComponent	(	const amrex::MultiFab &	field,
		int	component,
		std::string &	diagnostic
	)

Check that a selected cell-centered component is finite.

{
    diagnostic.clear();
    if (component < 0 || component >= field.nComp()) {
        diagnostic = "finite field validation received an invalid component";
        return false;
    }
    amrex::MultiFab invalid(field.boxArray(), field.DistributionMap(), 1, 0);
    invalid.setVal(amrex::Real(0.0));
    for (amrex::MFIter mfi(field, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
        const amrex::Box bx = mfi.tilebox();
        const auto values = field.const_array(mfi);
        const auto bad = invalid.array(mfi);
        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
            bad(i, j, k, 0) = amrex::Math::isfinite(values(i, j, k, component)) ?
                amrex::Real(0.0) : amrex::Real(1.0);
        });
    }
    if (invalid.sum(0) != amrex::Real(0.0)) {
        diagnostic = "selected field component must be finite in every valid cell";
        return false;
    }
    return true;
}

Referenced by erf_auxiliary::AuxiliaryInertTracer::advance_stage(), BuildAuxiliaryIntensiveState(), and erf_auxiliary::AuxiliaryInertTracer::initialize().

Here is the call graph for this function:

Here is the caller graph for this function:

ValidatePositiveFiniteComponent()

bool erf_auxiliary::ValidatePositiveFiniteComponent	(	const amrex::MultiFab &	field,
		int	component,
		std::string &	diagnostic
	)

Check that a selected cell-centered component is finite and strictly positive.

{
    diagnostic.clear();
    if (component < 0 || component >= field.nComp()) {
        diagnostic = "positive finite field validation received an invalid component";
        return false;
    }
    amrex::MultiFab invalid(field.boxArray(), field.DistributionMap(), 1, 0);
    invalid.setVal(amrex::Real(0.0));
    for (amrex::MFIter mfi(field, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
        const amrex::Box bx = mfi.tilebox();
        const auto values = field.const_array(mfi);
        const auto bad = invalid.array(mfi);
        amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
            const amrex::Real value = values(i, j, k, component);
            const bool valid = amrex::Math::isfinite(value) && value > amrex::Real(0.0);
            bad(i, j, k, 0) = valid ? amrex::Real(0.0) : amrex::Real(1.0);
        });
    }
    if (invalid.sum(0) != amrex::Real(0.0)) {
        diagnostic = "selected field component must be finite and positive in every valid cell";
        return false;
    }
    return true;
}

Referenced by erf_sbm::SBMTransport::advance_stage(), erf_auxiliary::AuxiliaryInertTracer::advance_stage(), BuildAuxiliaryIntensiveState(), ComputeMaxMappedOutgoingRate(), and erf_auxiliary::AuxiliaryInertTracer::initialize().

Here is the call graph for this function:

Here is the caller graph for this function: