ERF_MaterialProperties.H

Go to the documentation of this file.

#ifndef MATERIALPROPERTIES_H
#define MATERIALPROPERTIES_H
 
#include <limits>
 
#include <AMReX_Enum.H>
#include <AMReX_REAL.H>
#include <AMReX_Gpu.H>
#include <AMReX_MultiFab.H>
#include "ERF_NumericalConstants.H"
#include "ERF_MicrophysicsConstants.H"
 
/*! Materials */
namespace Species
{
    AMREX_ENUM(Name,
        H2O,
        ice,
        NaCl,
        NH42SO4,
        NH4HSO4,
        soil,
        water, // same as H2O, to test multispecies SDM
        agua,  // same as H2O, to test multispecies SDM
        none
    );
}
 
/*! Functions for computing saturation pressure and vapour fraction */
namespace saturation_funcs
{
    /*! \brief Null function for saturation pressure */
    AMREX_GPU_HOST
    void compute_saturation_pressure_null (  amrex::MultiFab&, const amrex::MultiFab& );
    /*! \brief Compute saturation pressure for moist air */
    AMREX_GPU_HOST
    void compute_saturation_pressure_H2O  (  amrex::MultiFab&, const amrex::MultiFab& );
    /*! \brief Compute saturation pressure for icy air */
    AMREX_GPU_HOST
    void compute_saturation_pressure_ice  (  amrex::MultiFab&, const amrex::MultiFab& );
 
    /*! \brief Null function for saturation vapour fraction */
    AMREX_GPU_HOST
    void compute_saturation_vapfrac_null ( amrex::MultiFab&,
                                           const amrex::MultiFab&,
                                           const amrex::MultiFab& );
 
    /*! \brief Compute saturation vapour fraction for moist air */
    AMREX_GPU_HOST
    void compute_saturation_vapfrac_H2O  ( amrex::MultiFab&,
                                           const amrex::MultiFab&,
                                           const amrex::MultiFab& );
    /*! \brief Compute saturation vapour fraction for icy air */
    AMREX_GPU_HOST
    void compute_saturation_vapfrac_ice  ( amrex::MultiFab&,
                                           const amrex::MultiFab&,
                                           const amrex::MultiFab& );
}
 
/**
 * @brief Base structure containing material property constants and methods for GPU usage.
 */
struct MaterialPropertiesCore { // For GPU
 
    static constexpr amrex::Real s_N_av = amrex::Real(6.02214076e23); /*!< Avogadro's number */
    static constexpr amrex::Real s_kb   = amrex::Real(1.380649e-23); /*!< Boltzmann constant [J K^{-1}] */
 
    /*! Molecular weight of air [kg mol^{-1}] */
    static constexpr amrex::Real s_mol_weight_air = amrex::Real(0.0289647);
    /*! Molecular length scale of air [m] */
    static constexpr amrex::Real s_sigma_air = amrex::Real(3.62e-10);
    /*! Temperature parameter for air [K] */
    static constexpr amrex::Real s_eps_air_k = 97;
 
#ifdef AMREX_USE_FLOAT
    amrex::Real m_density = std::numeric_limits<float>::max(); /*!< density */
    amrex::Real m_ionization = std::numeric_limits<float>::max(); /*!< ionization */
    amrex::Real m_mol_weight = std::numeric_limits<float>::max(); /*!< molecular weight (condensate) */
    amrex::Real m_lat_vap = std::numeric_limits<float>::max(); /*!< latent heat of vaporization */
    amrex::Real m_lat_fus = std::numeric_limits<float>::max(); /*!< latent heat of fusion */
    amrex::Real m_Rv = std::numeric_limits<float>::max(); /*!< gas constant for vapour of this material */
    amrex::Real m_Tc = std::numeric_limits<float>::max(); /*!< critical temperature for surface tension */
    amrex::Real m_Tb = std::numeric_limits<float>::max(); /*!< boiling temperature */
    amrex::Real m_Nav_by_molweight = std::numeric_limits<float>::max(); /*!< Avogadro number by vapour mol. weight */
#else
    amrex::Real m_density = std::numeric_limits<double>::max(); /*!< density */
    amrex::Real m_ionization = std::numeric_limits<double>::max(); /*!< ionization */
    amrex::Real m_mol_weight = std::numeric_limits<double>::max(); /*!< molecular weight (condensate) */
    amrex::Real m_lat_vap = std::numeric_limits<double>::max(); /*!< latent heat of vaporization */
    amrex::Real m_lat_fus = std::numeric_limits<double>::max(); /*!< latent heat of fusion */
    amrex::Real m_Rv = std::numeric_limits<double>::max(); /*!< gas constant for vapour of this material */
    amrex::Real m_Tc = std::numeric_limits<double>::max(); /*!< critical temperature for surface tension */
    amrex::Real m_Tb = std::numeric_limits<double>::max(); /*!< boiling temperature */
    amrex::Real m_Nav_by_molweight = std::numeric_limits<float>::max(); /*!< Avogadro number by vapour mol. weight */
#endif
    /*! vector of molar heat capacity polynomial coeffs */
    amrex::Real m_mol_Cp_coeffs[7] = {  std::numeric_limits<amrex::Real>::max(),
                                        std::numeric_limits<amrex::Real>::max(),
                                        std::numeric_limits<amrex::Real>::max(),
                                        std::numeric_limits<amrex::Real>::max(),
                                        std::numeric_limits<amrex::Real>::max(),
                                        std::numeric_limits<amrex::Real>::max(),
                                        std::numeric_limits<amrex::Real>::max() };
    bool m_is_soluble = false; /*!< is soluble in water? */
    bool m_is_water = false; /*!< is this water? */
    bool m_is_INP = false; /*!< can act as ice nucleating particle? */
 
    /*! \brief Return the molar heat capacity */
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    amrex::Real molarHeatCapacity( const amrex::Real a_T ) const
    {
        auto retval = zero;
        auto term = one;
        for (auto i = 0; i < 7; i++) {
            retval += (m_mol_Cp_coeffs[i] * term);
            term *= a_T;
        }
        return retval;
    }
 
    /*! \brief Compute the surface tension coeff given temperature */
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    amrex::Real surfaceTension( const amrex::Real a_T ) const
    {
        if (m_is_water) {
            return amrex::Real(0.076148325); // N m^{-1}
        } else {
            auto term1 = amrex::Real(6.0) * std::cbrt(s_N_av);
            auto term2 = std::cbrt((m_density/m_mol_weight)*(m_density/m_mol_weight));
            auto retval = molarHeatCapacity(a_T)*(term2/term1)*(m_Tc - a_T);
            return retval;
        }
    }
 
    /*! \brief Return the coeff of curvature given the temperature */
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    amrex::Real coeffCurv( const amrex::Real a_T ) const
    {
        auto surf_tens = surfaceTension(a_T);
        auto retval = 2 * surf_tens / (m_Rv * m_density);
        return retval;
    }
 
    /*! \brief Return the vapour pressure coefficient */
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    amrex::Real coeffVPSolute() const
    {
        auto retval = (amrex::Real(0.75)/PI)*(m_mol_weight/m_density);
        if (m_is_water) {
            retval = amrex::Real(4.3e-06);
        }
        return retval;
    }
 
    /*! \brief Return the molecular diffusion coefficient given temperature and pressure */
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    amrex::Real coeffMolecularDiffusion( const amrex::Real a_T, const amrex::Real a_P ) const
    {
        if (m_is_water) {
 
            return diffelq; // from ERF_Constants.H (amrex::Real(2.21e-05) [m^2 s^{-1}])
 
        } else {
 
            auto mol_weight_v = s_N_av/m_Nav_by_molweight;
            auto sigma = amrex::Real(1.18e-8) * std::cbrt(mol_weight_v/m_density);
            auto sigma_v_air = myhalf * (sigma + s_sigma_air);
 
            auto eps_v = amrex::Real(1.15) * m_Tb * s_kb;
            auto eps_v_air = std::sqrt( eps_v * s_eps_air_k * s_kb);
            auto T_star = a_T / (eps_v_air/s_kb);
 
            auto omegaD =    amrex::Real(1.06036) / std::pow(T_star, amrex::Real(0.1561))
                          +  amrex::Real(0.193)   / std::exp(amrex::Real(0.47635)*T_star)
                          +  amrex::Real(1.03587) / std::exp(amrex::Real(1.52996)*T_star)
                          +  amrex::Real(1.76474) / std::exp(amrex::Real(3.89411)*T_star);
 
            auto term1 = 3 * std::sqrt( s_kb*a_T * s_kb*a_T * s_kb*a_T);
            auto term2 = amrex::Real(8.0) * std::sqrt(2*PI) * a_P * sigma_v_air*sigma_v_air * omegaD;
            auto term3 = std::sqrt( m_Nav_by_molweight + s_N_av/s_mol_weight_air );
 
            auto retval = term1*term3/term2;
            return retval;
        }
    }
 
};
 
/*! Class for material properties */
struct MaterialProperties : public MaterialPropertiesCore {
 
    Species::Name m_name = Species::Name::none; /*!< name */
 
    /*! pointer to function to compute saturation pressure */
    decltype(saturation_funcs::compute_saturation_pressure_null)* m_saturation_pressure_func = nullptr;
 
    /*! pointer to function to compute vapour fraction */
    decltype(saturation_funcs::compute_saturation_vapfrac_null)* m_saturation_vapfrac_func = nullptr;
 
    /*! \brief Constructor */
    AMREX_GPU_HOST_DEVICE
    MaterialProperties ( const Species::Name& a_name );
 
    /*! \brief Copy constructor */
    AMREX_GPU_HOST_DEVICE
    MaterialProperties ( const MaterialProperties& a_matprop );
 
    /*! \brief Default destructor */
    AMREX_GPU_HOST_DEVICE
    ~MaterialProperties () = default;
 
    /*! \brief Print parameters */
    AMREX_GPU_HOST
    void print() const {
        amrex::Print() << "Material properties of " << amrex::getEnumNameString(m_name) << ":\n";
        amrex::Print() << "    density: " << m_density << "\n";
        amrex::Print() << "    ionization: " << m_ionization << "\n";
        amrex::Print() << "    mol. weight: " << m_mol_weight << "\n";
        amrex::Print() << "    latent heat (vap.): " << m_lat_vap << "\n";
        amrex::Print() << "    latent heat (fus.): " << m_lat_fus << "\n";
        amrex::Print() << "    Rv: " << m_Rv << "\n";
        amrex::Print() << "    Tc: " << m_Tc << "\n";
        amrex::Print() << "    Tb: " << m_Tb << "\n";
        amrex::Print() << "    N_av/mol.weight: " << m_Nav_by_molweight << "\n";
        amrex::Print() << "    mol. Cp coeffs: ";
        for (int i = 0; i < 7; i++) { amrex::Print() << m_mol_Cp_coeffs[i] << ", "; }
        amrex::Print() << "\n";
    }
 
    /*! \brief Set this material to H20 */
    AMREX_GPU_HOST_DEVICE
    void setProperties_H2O();
    /*! \brief Set this material to ice */
    AMREX_GPU_HOST_DEVICE
    void setProperties_ice();
    /*! \brief Set this material to water */
    AMREX_GPU_HOST_DEVICE
    void setProperties_water();
    /*! \brief Set this material to agua */
    AMREX_GPU_HOST_DEVICE
    void setProperties_agua();
    /*! \brief Set this material to NaCl */
    AMREX_GPU_HOST_DEVICE
    void setProperties_NaCl();
    /*! \brief Set this material to ammonium sulfate */
    AMREX_GPU_HOST_DEVICE
    void setProperties_NH42SO4();
    /*! \brief Set this material to ammonium bisulfate */
    AMREX_GPU_HOST_DEVICE
    void setProperties_NH4HSO4();
    /*! \brief Set this material to soil */
    AMREX_GPU_HOST_DEVICE
    void setProperties_soil();
 
    /*! \brief Compute saturation pressure */
    AMREX_GPU_HOST AMREX_FORCE_INLINE
    void computeSaturationPressure ( amrex::MultiFab& a_e, const amrex::MultiFab& a_T) const
    {
        m_saturation_pressure_func(a_e, a_T);
    }
 
    /*! \brief Compute saturation vapour fraction */
    AMREX_GPU_HOST AMREX_FORCE_INLINE
    void computeSaturationVapFrac  ( amrex::MultiFab&         a_S,
                                     const amrex::MultiFab&   a_T,
                                     const amrex::MultiFab&   a_p ) const
    {
        m_saturation_vapfrac_func(a_S, a_T, a_p);
    }
};
 
#endif