CloudChamberBudget Class Reference

#include <ERF_CloudChamberBudget.H>

Collaboration diagram for CloudChamberBudget:

Public Types
enum	Scalar : int { RhoTheta = 0 , RhoQv = 1 , RhoQc = 2 , NumScalars = 3 }

Public Member Functions
	CloudChamberBudget (int interval)
bool	enabled () const noexcept
bool	due (int step) const noexcept
void	set_initial_state (const amrex::MultiFab &cons, const amrex::Geometry &geom, int step=0, double time=0.0)
void	capture_stage (int scalar, int nrk, amrex::Real dt, const amrex::MultiFab &xflux, const amrex::MultiFab &yflux, const amrex::MultiFab &zflux, const amrex::Geometry &geom, int flux_comp=0, bool trapezoidal=true)
std::array< amrex::Real, NumScalars >	state_integrals (const amrex::MultiFab &cons, const amrex::Geometry &geom) const
void	record_internal_source (const std::array< amrex::Real, NumScalars > &before, const std::array< amrex::Real, NumScalars > &after)
void	report (int step, double time, const amrex::MultiFab &cons, const amrex::Geometry &geom, bool cloudy)

Static Public Member Functions
static const char *	row_status (Scalar scalar, bool cloudy, bool closes) noexcept

Static Public Attributes
static constexpr int	NumFaces = 2 * AMREX_SPACEDIM

Private Attributes
int	m_interval = 0
int	m_interval_start_step = 0
double	m_interval_start_time = 0.0
bool	m_have_initial = false
std::array< bool, NumScalars >	m_have_stage0 {}
std::array< amrex::Real, NumScalars >	m_initial_state {}
std::array< amrex::Real, NumScalars >	m_internal_source {}
std::array< std::array< amrex::Real, NumFaces >, NumScalars >	m_stage0 {}
std::array< std::array< amrex::Real, NumFaces >, NumScalars >	m_cumulative {}

Detailed Description

Low-overhead Stage 1 conserved-scalar budget accumulator.

Fluxes are positive in the positive coordinate direction. The reported boundary balance is therefore low-minus-high for each coordinate pair. Local face sums are held until a configured report interval, then one packed MPI reduction produces the six-face values. For the anelastic RK2 path, nrk=0 is retained and nrk=1 applies the trapezoidal stage weight; this avoids counting both RHS evaluations as two physical timesteps. For the anelastic MidPoint path there is nothing to average: the first stage advances only half a timestep and the second advances the whole timestep with the source evaluated there, so only nrk=1 contributes.

Member Enumeration Documentation

Scalar

enum CloudChamberBudget::Scalar : int

Enumerator
RhoTheta
RhoQv
RhoQc
NumScalars

: int { RhoTheta = 0, RhoQv = 1, RhoQc = 2, NumScalars = 3 };

Constructor & Destructor Documentation

CloudChamberBudget()

CloudChamberBudget::CloudChamberBudget ( int  interval )

inlineexplicit

: m_interval(interval) {}

Member Function Documentation

capture_stage()

void CloudChamberBudget::capture_stage ( int  scalar, int  nrk, amrex::Real  dt, const amrex::MultiFab &  xflux, const amrex::MultiFab &  yflux, const amrex::MultiFab &  zflux, const amrex::Geometry &  geom, int  flux_comp = 0, bool  trapezoidal = true  )

inline

Integrate boundary fluxes for a specific Runge-Kutta stage.

Parameters

scalar	Index of the scalar being tracked
nrk	Runge-Kutta stage index
dt	Timestep size
xflux	Fluxes in the x-direction
yflux	Fluxes in the y-direction
zflux	Fluxes in the z-direction
geom	Geometry used for face area calculation
flux_comp	Component index of the flux
trapezoidal	True for the RK2 stages, whose two sources are averaged; false for the MidPoint stages, where only nrk=1 counts

    {
        if (!enabled() || scalar < 0 || scalar >= NumScalars) { return; }
        std::array<amrex::Real, NumFaces> current{};
        for (int d = 0; d < AMREX_SPACEDIM; ++d) {
            amrex::Box low = geom.Domain();
            amrex::Box high = low;
            low.surroundingNodes(d);
            high.surroundingNodes(d);
            low.setSmall(d, geom.Domain().smallEnd(d));
            low.setBig(d, geom.Domain().smallEnd(d));
            high.setSmall(d, geom.Domain().bigEnd(d)+1);
            high.setBig(d, geom.Domain().bigEnd(d)+1);
            const amrex::MultiFab* flux = (d == 0) ? &xflux : ((d == 1) ? &yflux : &zflux);
            const amrex::Real area = (d == 0) ? geom.CellSize(1)*geom.CellSize(2) :
                                     ((d == 1) ? geom.CellSize(0)*geom.CellSize(2) :
                                                 geom.CellSize(0)*geom.CellSize(1));
            current[2*d]   = flux->sum(low, flux_comp, true) * area * dt;
            current[2*d+1] = flux->sum(high, flux_comp, true) * area * dt;
        }
 
        if (!trapezoidal) {
            // MidPoint: the first stage spans only dt/2 and exists solely to make the
            // midpoint state, so the second stage alone carries the step's flux.
            if (nrk == 1) {
                for (int face = 0; face < NumFaces; ++face) {
                    m_cumulative[scalar][face] += current[face];
                }
            }
            m_have_stage0[scalar] = false;
        } else if (nrk == 0) {
            m_stage0[scalar] = current;
            m_have_stage0[scalar] = true;
        } else if (nrk == 1 && m_have_stage0[scalar]) {
            for (int face = 0; face < NumFaces; ++face) {
                m_cumulative[scalar][face] +=
                    amrex::Real(0.5) * (m_stage0[scalar][face] + current[face]);
            }
            m_have_stage0[scalar] = false;
        }
    }

Referenced by erf_slow_rhs_post(), and erf_slow_rhs_pre().

Here is the call graph for this function:

Here is the caller graph for this function:

due()

bool CloudChamberBudget::due ( int  step ) const

inlinenoexcept

Determine if a budget report is due for the given timestep.

Parameters

step	Current timestep

Returns

True if a report is due, false otherwise

                                       {
        return enabled() && step > 0 && (step % m_interval == 0);
    }

Referenced by report().

Here is the call graph for this function:

Here is the caller graph for this function:

enabled()

bool CloudChamberBudget::enabled ( ) const

inlinenoexcept

{ return m_interval > 0; }

Referenced by capture_stage(), and due().

Here is the caller graph for this function:

record_internal_source()

void CloudChamberBudget::record_internal_source ( const std::array< amrex::Real, NumScalars > &  before, const std::array< amrex::Real, NumScalars > &  after  )

inline

Accumulate internal source terms for the tracked scalars.

Parameters

before	State of scalars before source application
after	State of scalars after source application

    {
        for (int s = 0; s < NumScalars; ++s) {
            m_internal_source[s] += after[s] - before[s];
        }
    }

report()

void CloudChamberBudget::report ( int  step, double  time, const amrex::MultiFab &  cons, const amrex::Geometry &  geom, bool  cloudy  )

inline

Perform global reduction and write the budget report to file.

Parameters

step	Current timestep
time	Current simulation time
cons	Current conserved variable MultiFab
geom	Geometry used for state integration
cloudy	Flag indicating if the state is cloudy

    {
        if (!due(step) || !m_have_initial) { return; }
        const auto current = state_integrals(cons, geom);
        std::array<amrex::Real, NumScalars*NumFaces + NumScalars*3> packed{};
        int index = 0;
        for (int s = 0; s < NumScalars; ++s) {
            for (int f = 0; f < NumFaces; ++f) { packed[index++] = m_cumulative[s][f]; }
        }
        for (int s = 0; s < NumScalars; ++s) { packed[index++] = m_initial_state[s]; }
        for (int s = 0; s < NumScalars; ++s) { packed[index++] = current[s]; }
        for (int s = 0; s < NumScalars; ++s) { packed[index++] = m_internal_source[s]; }
        amrex::ParallelDescriptor::ReduceRealSum(packed.data(), static_cast<int>(packed.size()));
 
        std::array<std::array<amrex::Real, NumFaces>, NumScalars> global_faces{};
        index = 0;
        for (int s = 0; s < NumScalars; ++s) {
            for (int f = 0; f < NumFaces; ++f) { global_faces[s][f] = packed[index++]; }
        }
        std::array<amrex::Real, NumScalars> initial{};
        std::array<amrex::Real, NumScalars> global_current{};
        index = NumScalars*NumFaces;
        for (int s = 0; s < NumScalars; ++s) { initial[s] = packed[index++]; }
        for (int s = 0; s < NumScalars; ++s) { global_current[s] = packed[index++]; }
        std::array<amrex::Real, NumScalars> global_internal{};
        for (int s = 0; s < NumScalars; ++s) { global_internal[s] = packed[index++]; }
 
        if (amrex::ParallelDescriptor::IOProcessor()) {
            std::ofstream out("cloud_chamber_budget.dat", std::ios::app);
            if (out.tellp() == std::streampos(0)) {
                out << "# start_step start_time end_step end_time scalar units "
                       "xlo xhi ylo yhi zlo zhi net_boundary volume_change "
                       "internal_source residual tolerance status\n";
            }
            out << std::setprecision(17);
            const char* names[NumScalars] = {"rhoTheta", "water_vapor", "cloud_water"};
            const char* units[NumScalars] = {"conserved_scalar", "mixing_ratio_density", "mixing_ratio_density"};
            for (int s = 0; s < NumScalars; ++s) {
                std::array<amrex::Real, NumFaces> faces = global_faces[s];
                amrex::Real net = faces[0] - faces[1] + faces[2] - faces[3] + faces[4] - faces[5];
                const amrex::Real change = global_current[s] - initial[s];
                const amrex::Real internal = global_internal[s];
            const amrex::Real residual = change - net - internal;
            const amrex::Real scale = std::max({amrex::Real(1.0),
                    std::abs(initial[s]), std::abs(global_current[s]),
                    std::abs(net), std::abs(internal)});
                const amrex::Real effective_operations =
                    static_cast<amrex::Real>(geom.Domain().numPts()) +
                    static_cast<amrex::Real>(NumFaces);
                amrex::Real boundary_scale = amrex::Real(0.0);
                for (const auto value : faces) { boundary_scale += std::abs(value); }
                const amrex::Real tolerance = amrex::Real(128.0) *
                    std::numeric_limits<amrex::Real>::epsilon() *
                    (effective_operations * scale + boundary_scale + std::abs(internal));
                const bool status = std::abs(residual) <= tolerance;
                const char* status_text = row_status(
                    static_cast<Scalar>(s), cloudy, status);
                out << m_interval_start_step << ' ' << m_interval_start_time << ' '
                    << step << ' ' << time << ' ' << names[s] << ' ' << units[s];
                for (const auto value : faces) { out << ' ' << value; }
                out << ' ' << net << ' ' << change << ' ' << internal << ' '
                    << residual << ' ' << tolerance << ' ' << status_text << '\n';
            }
            std::array<amrex::Real, NumFaces> total_faces{};
            for (int f = 0; f < NumFaces; ++f) {
                total_faces[f] = global_faces[RhoQv][f] + global_faces[RhoQc][f];
            }
            const amrex::Real total_net = total_faces[0] - total_faces[1] +
                total_faces[2] - total_faces[3] + total_faces[4] - total_faces[5];
            const amrex::Real total_change = global_current[RhoQv] + global_current[RhoQc] -
                initial[RhoQv] - initial[RhoQc];
            const amrex::Real total_internal = global_internal[RhoQv] + global_internal[RhoQc];
            const amrex::Real total_residual = total_change - total_net - total_internal;
            const amrex::Real total_scale = std::max({amrex::Real(1.0), std::abs(total_change),
                std::abs(total_net), std::abs(total_internal)});
            amrex::Real total_boundary_scale = amrex::Real(0.0);
            for (const auto value : total_faces) { total_boundary_scale += std::abs(value); }
            const amrex::Real total_tolerance = amrex::Real(128.0) *
                std::numeric_limits<amrex::Real>::epsilon() *
                (static_cast<amrex::Real>(geom.Domain().numPts() + NumFaces) * total_scale +
                 total_boundary_scale + std::abs(total_internal));
            const bool total_status = std::abs(total_residual) <= total_tolerance;
            out << m_interval_start_step << ' ' << m_interval_start_time << ' '
                << step << ' ' << time << " total_nonprecipitating_water mixing_ratio_density";
            for (const auto value : total_faces) { out << ' ' << value; }
            out << ' ' << total_net << ' ' << total_change << ' ' << total_internal << ' '
                << total_residual << ' ' << total_tolerance << ' '
                << (total_status ? "PASS" : "FAIL") << '\n';
        }
        // The next report is an interval-local balance beginning at the
        // state just reported.  Keep this reference rank-local: assigning the
        // MPI-global value here would make every rank contribute the same
        // global reference at the next reduction.
        m_initial_state = current;
        m_cumulative = {};
        m_internal_source = {};
        m_stage0 = {};
        m_have_stage0 = {};
        m_interval_start_step = step;
        m_interval_start_time = time;
    }

Here is the call graph for this function:

row_status()

static const char* CloudChamberBudget::row_status ( Scalar  scalar, bool  cloudy, bool  closes  )

inlinestaticnoexcept

    {
        if (cloudy && scalar == RhoTheta) { return "UNSUPPORTED_SOURCE"; }
        return closes ? "PASS" : "FAIL";
    }

Referenced by report().

Here is the caller graph for this function:

set_initial_state()

void CloudChamberBudget::set_initial_state ( const amrex::MultiFab &  cons, const amrex::Geometry &  geom, int  step = 0, double  time = 0.0  )

inline

Initialize the base state for the budget accumulation period.

Parameters

cons	Conserved variable MultiFab
geom	Geometry used for volume calculation
step	Timestep at the start of the interval
time	Simulation time at the start of the interval

    {
        if (m_have_initial) { return; }
        const amrex::Real volume = geom.CellSize(0) * geom.CellSize(1) * geom.CellSize(2);
        for (int s = 0; s < NumScalars; ++s) {
            const int comp = (s == RhoTheta) ? RhoTheta_comp : RhoQ1_comp + (s-1);
            m_initial_state[s] = (comp < cons.nComp()) ? cons.sum(comp, true) * volume : amrex::Real(0.0);
        }
        m_interval_start_step = step;
        m_interval_start_time = time;
        m_have_initial = true;
    }

state_integrals()

std::array<amrex::Real, NumScalars> CloudChamberBudget::state_integrals ( const amrex::MultiFab &  cons, const amrex::Geometry &  geom  ) const

inline

    {
        const amrex::Real volume = geom.CellSize(0) * geom.CellSize(1) * geom.CellSize(2);
        std::array<amrex::Real, NumScalars> result{};
        result[RhoTheta] = cons.sum(RhoTheta_comp, true) * volume;
        result[RhoQv] = (RhoQ1_comp < cons.nComp()) ? cons.sum(RhoQ1_comp, true) * volume : amrex::Real(0.0);
        result[RhoQc] = (RhoQ2_comp < cons.nComp()) ? cons.sum(RhoQ2_comp, true) * volume : amrex::Real(0.0);
        return result;
    }

Referenced by report().

Here is the caller graph for this function:

Member Data Documentation

m_cumulative

std::array<std::array<amrex::Real, NumFaces>, NumScalars> CloudChamberBudget::m_cumulative {}

private

Referenced by capture_stage(), and report().

m_have_initial

bool CloudChamberBudget::m_have_initial = false

private

Referenced by report(), and set_initial_state().

m_have_stage0

std::array<bool, NumScalars> CloudChamberBudget::m_have_stage0 {}

private

Referenced by capture_stage(), and report().

m_initial_state

std::array<amrex::Real, NumScalars> CloudChamberBudget::m_initial_state {}

private

Referenced by report(), and set_initial_state().

m_internal_source

std::array<amrex::Real, NumScalars> CloudChamberBudget::m_internal_source {}

private

Referenced by record_internal_source(), and report().

m_interval

int CloudChamberBudget::m_interval = 0

private

Referenced by due(), and enabled().

m_interval_start_step

int CloudChamberBudget::m_interval_start_step = 0

private

Referenced by report(), and set_initial_state().

m_interval_start_time

double CloudChamberBudget::m_interval_start_time = 0.0

private

Referenced by report(), and set_initial_state().

m_stage0

std::array<std::array<amrex::Real, NumFaces>, NumScalars> CloudChamberBudget::m_stage0 {}

private

Referenced by capture_stage(), and report().

NumFaces

constexpr int CloudChamberBudget::NumFaces = 2 * AMREX_SPACEDIM

staticconstexpr

Referenced by capture_stage(), and report().

---

The documentation for this class was generated from the following file:

• Source/Prob/ERF_CloudChamberBudget.H