ERF_FireUtils.H File Reference

Utility functions for fire model atmospheric state extraction. More...

#include <AMReX_MultiFab.H>
#include <AMReX_Geometry.H>
#include <AMReX_MFIter.H>
#include <ERF_MicrophysicsUtils.H>
#include <ERF_EOS.H>
#include <ERF_IndexDefines.H>
#include <ERF_Constants.H>

Include dependency graph for ERF_FireUtils.H:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions
void	compute_rh_from_conservative (amrex::MultiFab &RH_out, const amrex::MultiFab &S_atm, const amrex::Geometry &)
	Compute relative humidity from conserved atmospheric variables at k=0. More...
void	compute_t_from_conservative (amrex::MultiFab &T_out, const amrex::MultiFab &S_atm)
	Air temperature [K] at k = 0 from the conserved state. More...
void	fire_copy_k0_plane (amrex::MultiFab &dst_2d, const amrex::MultiFab &src_3d)
void	fire_host_copy (amrex::MultiFab &dst, const amrex::MultiFab &src, int ncomp=1, int nghost=0)
	Copy a MultiFab between device and host-accessible memory, for an engine whose fire physics runs on the host. More...

Detailed Description

Utility functions for fire model atmospheric state extraction.

Function Documentation

compute_rh_from_conservative()

void compute_rh_from_conservative ( amrex::MultiFab &  RH_out, const amrex::MultiFab &  S_atm, const amrex::Geometry &    )

inline

Compute relative humidity from conserved atmospheric variables at k=0.

Computes RH from pressure, temperature, and water vapor mixing ratio at the first atmospheric level (k=0). Uses the ERF microphysics saturation vapor pressure functions.

Parameters

[out]	RH_out	Output relative humidity MultiFab [fraction 0-1] Must be allocated with ncomp=1 on same grid as inputs
[in]	S_atm	Conservative state at k=0 Contains RhoTheta (comp 1) and RhoQ1 (comp 4)
[in]	geom	Geometry at level 0

Note

• Output is clamped to [0, 1] range
• Assumes S_atm contains standard conserved variables (Rho, RhoTheta, RhoQ1)
• Temperature extracted from RhoTheta using standard EOS
• Pressure computed from RhoTheta and qv using getPgivenRTh

{
    using namespace amrex;
 
    // Precompute rdOcp = R_d / Cp_d for use in Exner function
    // Cp_d = 1004.5 J/(kg-K), R_d = 287.0 J/(kg-K)
    constexpr Real rdOcp = R_d / Cp_d;
 
    // Iterate over all boxes
    for (MFIter mfi(RH_out); mfi.isValid(); ++mfi) {
        const Box& bx = mfi.tilebox();
        Array4<Real>       rh_arr = RH_out.array(mfi);
        Array4<const Real> S_arr  = S_atm.const_array(mfi);
 
        // Kernel to compute RH at k=0
        ParallelFor(bx, [=] AMREX_GPU_DEVICE(const IntVect& iv) {
            int i = iv[0];
            int j = iv[1];
            int k = 0;  // First atmospheric level
 
            // Extract conserved variables at k=0
            Real rho     = S_arr(i, j, k, Rho_comp);           // Density
            Real rho_th  = S_arr(i, j, k, RhoTheta_comp);     // RhoTheta
 
            // Compute mixing ratios (guard against moisture not being present)
            Real qv = 0.0_rt;  // Water vapor mixing ratio [kg/kg]
            if (RhoQ1_comp < S_arr.nComp()) {
                Real rho_qv  = S_arr(i, j, k, RhoQ1_comp);    // RhoQ1 (water vapor)
                qv = rho_qv / rho;
            }
 
            // Compute potential temperature
            Real theta = rho_th / rho;  // Potential temperature [K]
 
            // Compute pressure from RhoTheta and qv using standard EOS
            Real p = getPgivenRTh(rho_th, qv);  // Pressure [Pa]
 
            // Convert pressure to mbar for saturation vapor pressure functions
            Real p_mbar = p * 0.01_rt;  // Pa to mbar (1 mbar = 100 Pa)
 
            // Compute temperature from potential temperature
            // T = theta * (p / p_0)^(R_d / c_p)
            // Use Exner function: pi = (p / p_0)^(R_d / c_p)
            Real exner = getExnergivenP(p, rdOcp);
            Real T_K = theta * exner;
 
            // Compute saturation vapor pressure [mbar]
            Real e_sat = erf_esatw(T_K, false);  // Use Flatau polynomial
 
            // Compute saturation mixing ratio [kg/kg]
            Real q_sat = erf_qsat_from_esat(e_sat, p_mbar);
 
            // Compute relative humidity as RH = qv / qsat
            Real RH = 0.0_rt;
            if (q_sat > 1.0e-10_rt) {
                RH = qv / q_sat;
            }
 
            // Clamp RH to valid range [0, 1]
            RH = amrex::min(1.0_rt, amrex::max(0.0_rt, RH));
 
            rh_arr(i, j, k) = RH;
        });
    }
}

Referenced by FireCoupler::advance().

Here is the call graph for this function:

Here is the caller graph for this function:

compute_t_from_conservative()

void compute_t_from_conservative ( amrex::MultiFab &  T_out, const amrex::MultiFab &  S_atm  )

inline

Air temperature [K] at k = 0 from the conserved state.

T = theta * (p / p_0)^(R_d / c_p), with p from RhoTheta and qv as in compute_rh_from_conservative. Fuel-moisture curves and ignition criteria take an air temperature, not the potential temperature of the first cell, which at 840 hPa is 14-16 K above it.

Parameters

[out]	T_out	2D MultiFab (one cell in z) on the level's boxes
[in]	S_atm	Conserved state (Rho, RhoTheta, RhoQ1 when moist)

{
    using namespace amrex;
    constexpr Real rdOcp = R_d / Cp_d;
    for (MFIter mfi(T_out); mfi.isValid(); ++mfi) {
        const Box& bx = mfi.tilebox();
        Array4<Real>       t_arr = T_out.array(mfi);
        Array4<const Real> S_arr = S_atm.const_array(mfi);
        ParallelFor(bx, [=] AMREX_GPU_DEVICE(const IntVect& iv) {
            const int i = iv[0];
            const int j = iv[1];
            const int k = 0;
            const Real rho    = S_arr(i, j, k, Rho_comp);
            const Real rho_th = S_arr(i, j, k, RhoTheta_comp);
            Real qv = 0.0_rt;
            if (RhoQ1_comp < S_arr.nComp()) {
                qv = S_arr(i, j, k, RhoQ1_comp) / rho;
            }
            const Real p     = getPgivenRTh(rho_th, qv);
            const Real exner = getExnergivenP(p, rdOcp);
            t_arr(i, j, k)   = (rho_th / rho) * exner;
        });
    }
}

Referenced by FireCoupler::advance().

Here is the call graph for this function:

Here is the caller graph for this function:

fire_copy_k0_plane()

void fire_copy_k0_plane ( amrex::MultiFab &  dst_2d, const amrex::MultiFab &  src_3d  )

inline

Copy the k = 0 plane of a cell-centred atmospheric field into a MultiFab on the same boxes coarsened to one cell in z. A namespace-scope function because nvcc rejects a device kernel inside a private member such as ERF::Advance.

{
    for (amrex::MFIter mfi(dst_2d); mfi.isValid(); ++mfi) {
        amrex::Array4<amrex::Real> d = dst_2d.array(mfi);
        amrex::Array4<const amrex::Real> s = src_3d.const_array(mfi);
        amrex::ParallelFor(mfi.tilebox(), [=] AMREX_GPU_DEVICE (const amrex::IntVect& iv) {
            d(iv[0], iv[1], 0) = s(iv[0], iv[1], 0);
        });
    }
}

Here is the call graph for this function:

fire_host_copy()

void fire_host_copy ( amrex::MultiFab &  dst, const amrex::MultiFab &  src, int  ncomp = 1, int  nghost = 0  )

inline

Copy a MultiFab between device and host-accessible memory, for an engine whose fire physics runs on the host.

The memory contract of ERF_FireEngine.H: the coupling's fields are in the default arena (device memory on a GPU build) and are only valid during the call that hands them over; an engine's returned fields must be readable there too. dst and src have the same BoxArray and DistributionMapping; ncomp components and nghost ghost cells are copied from component 0. dst belongs to the caller (allocate a host copy once with amrex::MFInfo().SetArena(amrex::The_Pinned_Arena()), or The_Managed_Arena(), and reuse it). The copy completes before the call returns: the GPU stream is synchronised, so the host may read dst, or the device read it, right away. On a CPU build it is a plain copy.

{
    AMREX_ALWAYS_ASSERT(dst.boxArray() == src.boxArray() && dst.DistributionMap() == src.DistributionMap());
    amrex::MultiFab::Copy(dst, src, 0, 0, ncomp, nghost);
    amrex::Gpu::streamSynchronize();
}

Here is the call graph for this function: