diff --git a/Source/ERF_Constants.H b/Source/ERF_Constants.H index a13e5cae2..eda03354f 100644 --- a/Source/ERF_Constants.H +++ b/Source/ERF_Constants.H @@ -15,6 +15,7 @@ constexpr amrex::Real eightthirds = amrex::Real(8.0/3.0); constexpr amrex::Real PI = amrex::Real(3.14159265358979323846264338327950288); constexpr amrex::Real PIoTwo = PI/two; +constexpr amrex::Real four_thirds_pi = four/three*PI; // Physical Constants constexpr amrex::Real R_d = amrex::Real(287.0); // dry air constant for dry air [J/(kg-K)] diff --git a/Source/MaterialProperties/ERF_TerminalVelocity.H b/Source/MaterialProperties/ERF_TerminalVelocity.H index 879611fbf..c17e64e43 100644 --- a/Source/MaterialProperties/ERF_TerminalVelocity.H +++ b/Source/MaterialProperties/ERF_TerminalVelocity.H @@ -8,13 +8,28 @@ template struct TerminalVelocity { - RT m_rho; + RT m_rho_w; /*!< density of water */ + + /*! \brief Compute viscosity coefficient - Pruppacher & Klett, 1997*/ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + auto viscCoeff( const RT a_T /*!< temperature [K] */ ) const noexcept + { + RT T_degC = a_T - RT(273.15); // [K] -> [deg Celsius] + RT viscosity = RT(0.0); + if (T_degC >= RT(0)) { + viscosity = (RT(1.7180) + RT(4.9e-3)*T_degC) * RT(1.0e-4); + } else { + viscosity = ( RT(1.7180) + RT(4.9e-3)*T_degC + - RT(1.2e-5)*T_degC*T_degC ) * RT(1.0e-4); + } + return viscosity; + } /*! \brief Terminal velocity - Rogers & Yau, 1989 */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE RT RogersYau( const RT a_r) const noexcept { - static const RT k1 = amrex::Real(1.233e8); // m^{-1} s^{-1} + static const RT k1 = RT(1.233e8); // m^{-1} s^{-1} return (k1 * a_r * a_r); } @@ -22,9 +37,9 @@ struct TerminalVelocity AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE RT AtlasUlbrich( const RT a_r) const noexcept { - static const RT a = amrex::Real(3.778); - static const RT b = amrex::Real(0.67); - RT D = 2*a_r*1000; // meter -> mm + static const RT a = RT(3.778); + static const RT b = RT(0.67); + RT D = RT(2)*a_r*RT(1000); // meter -> mm return a*(std::exp(b*std::log(D))); } @@ -41,63 +56,56 @@ struct TerminalVelocity const RT a_P, /*!< pressure */ const RT a_T /*!< temperature */) const noexcept { - RT lb4l = amrex::Real(6.62E-6); // base length[cm] for mean free path of air molecules - RT visb4l = amrex::Real(1.818E-4); // base dynamic viscosity[g/(cm*s)] for mean free path of air molecules - RT pb4l = amrex::Real(1013.25); // base pressure[hPa] for mean free path of air molecules - RT tb4l = amrex::Real(293.15); // base temperature[20C] for mean free path of air molecules + RT lb4l = RT(6.62E-6); // base length[cm] for mean free path of air molecules + RT visb4l = RT(1.818E-4); // base dynamic viscosity[g/(cm*s)] for mean free path of air molecules + RT pb4l = RT(1013.25); // base pressure[hPa] for mean free path of air molecules + RT tb4l = RT(293.15); // base temperature[20C] for mean free path of air molecules - RT P_hPa = a_P / amrex::Real(100.0); // [Pa] -> [hPa] - RT T_degC = a_T - amrex::Real(273.15); // [K] -> [deg Celsius] - RT diameter_cm = two*a_r*amrex::Real(100.0); // [m] -> [cm] - RT rho_mat_cgs = m_rho/amrex::Real(1000.0); // [kg/m^3] -> [g/cm^3] - RT rho_air_cgs = a_rho/amrex::Real(1000.0); // [kg/m^3] -> [g/cm^3] - RT grav_cgs = CONST_GRAV * 100; // [m/s^2] -> [cm/s^2] + RT P_hPa = a_P / RT(100.0); // [Pa] -> [hPa] + RT diameter_cm = RT(2.0)*a_r*RT(100.0); // [m] -> [cm] + RT rho_mat_cgs = m_rho_w/RT(1000.0); // [kg/m^3] -> [g/cm^3] + RT rho_air_cgs = a_rho/RT(1000.0); // [kg/m^3] -> [g/cm^3] + RT grav_cgs = CONST_GRAV * RT(100); // [m/s^2] -> [cm/s^2] RT gxdrow = grav_cgs * ( rho_mat_cgs - rho_air_cgs ); - RT viscosity = zero; - if (T_degC >= 0) { - viscosity = (amrex::Real(1.7180) + amrex::Real(4.9e-3)*T_degC) * amrex::Real(1.0e-4); - } else { - viscosity = ( amrex::Real(1.7180) + amrex::Real(4.9e-3)*T_degC - - amrex::Real(1.2e-5)*T_degC*T_degC ) * amrex::Real(1.0e-4); - } + RT viscosity = viscCoeff(a_T); - RT v_terminal = zero; - if (diameter_cm < amrex::Real(1.9e-3)) { + RT v_terminal = RT(0.0); + if (diameter_cm < RT(1.9e-3)) { // small cloud droplets RT sd_l = lb4l * (viscosity/visb4l) * (pb4l/P_hPa) * std::sqrt(a_T/tb4l); - RT c1 = gxdrow / (amrex::Real(18.0)*viscosity); - RT csc = one + amrex::Real(2.510) * (sd_l/diameter_cm); + RT c1 = gxdrow / (RT(18.0)*viscosity); + RT csc = RT(1.0) + RT(2.510) * (sd_l/diameter_cm); v_terminal = c1 * csc * diameter_cm * diameter_cm; - } else if (diameter_cm < amrex::Real(1.07e-1)) { + } else if (diameter_cm < RT(1.07e-1)) { // large cloud droplets and small raindrops RT sd_l = lb4l * (viscosity/visb4l) * (pb4l/P_hPa) * std::sqrt(a_T/tb4l); - RT c2 = rho_air_cgs * (amrex::Real(4.0)*gxdrow)/(three*viscosity*viscosity); + RT c2 = rho_air_cgs * (RT(4.0)*gxdrow)/(RT(three)*viscosity*viscosity); RT nda = c2 * diameter_cm * diameter_cm * diameter_cm; - RT csc = one + amrex::Real(2.510) * ( sd_l / diameter_cm ); + RT csc = RT(1.0) + RT(2.510) * ( sd_l / diameter_cm ); RT nre = csc * std::exp(vz_b_x(std::log(nda))); v_terminal = (viscosity*nre)/(rho_air_cgs*diameter_cm); } else { // large raindrops - RT tau = one - a_T/amrex::Real(647.0960); - RT sigma = amrex::Real(0.2358) * std::exp( amrex::Real(1.2560)*std::log(tau) ) - * ( one - amrex::Real(0.6250)*tau ); - if( a_T g/s^2 - RT c3 = (amrex::Real(4.0)*gxdrow)/(three*sigma); + sigma = sigma * RT(1.E+3); // N/m => g/s^2 + RT c3 = (RT(4.0)*gxdrow)/(RT(three)*sigma); RT bond = c3 * diameter_cm * diameter_cm; RT npp = (rho_air_cgs*rho_air_cgs) * (sigma*sigma*sigma) / (gxdrow*viscosity*viscosity*viscosity*viscosity); - npp = std::exp( std::log(npp)/amrex::Real(6.0) ); + npp = std::exp( std::log(npp)/RT(6.0) ); RT nre = npp * std::exp(vz_c_x(std::log(bond*npp))); v_terminal = (viscosity*nre)/(rho_air_cgs*diameter_cm); } - v_terminal /= amrex::Real(100.0); // [cm/s] -> [m/s] + v_terminal /= RT(100.0); // [cm/s] -> [m/s] return v_terminal; } @@ -110,13 +118,13 @@ struct TerminalVelocity RT x1 = a_x; RT x2 = x1 * x1; RT x3 = x1 * x2; - RT y = - amrex::Real(0.318657E+1) - + amrex::Real(0.9926960) * x1 - - amrex::Real(0.153193E-2) * x2 - - amrex::Real(0.987059E-3) * x3 - - amrex::Real(0.578878e-3) * x2*x2 - + amrex::Real(0.855176E-4) * x3*x2 - - amrex::Real(0.327815E-5) * x3*x3; + RT y = - RT(0.318657E+1) + + RT(0.9926960) * x1 + - RT(0.153193E-2) * x2 + - RT(0.987059E-3) * x3 + - RT(0.578878e-3) * x2*x2 + + RT(0.855176E-4) * x3*x2 + - RT(0.327815E-5) * x3*x3; return y; } @@ -127,12 +135,12 @@ struct TerminalVelocity RT x1 = a_x; RT x2 = x1 * x1; RT x3 = x1 * x2; - RT y = - amrex::Real(0.500015E+1) - + amrex::Real(0.523778E+1) * x1 - - amrex::Real(0.204914E+1) * x2 - + amrex::Real(0.47529400) * x3 - - amrex::Real(0.542819E-1) * x2*x2 - + amrex::Real(0.238449E-2) * x3*x2; + RT y = - RT(0.500015E+1) + + RT(0.523778E+1) * x1 + - RT(0.204914E+1) * x2 + + RT(0.47529400) * x3 + - RT(0.542819E-1) * x2*x2 + + RT(0.238449E-2) * x3*x2; return y; } diff --git a/Source/Microphysics/SuperDropletsMoist/ERF_SuperDropletsMoist.H b/Source/Microphysics/SuperDropletsMoist/ERF_SuperDropletsMoist.H index 265bc0cc5..70230294b 100644 --- a/Source/Microphysics/SuperDropletsMoist/ERF_SuperDropletsMoist.H +++ b/Source/Microphysics/SuperDropletsMoist/ERF_SuperDropletsMoist.H @@ -313,66 +313,36 @@ class SuperDropletsMoist : public NullMoistLagrangian { amrex::MultiFab& a_mf ) const override { AMREX_ASSERT(a_mf.nComp() >= 1); - for (int v = 1; v < m_num_species; v++) { - { - std::string varname = "qv_" + amrex::getEnumNameString(m_species[v]); - if (a_name == varname) { - amrex::MultiFab::Copy( a_mf, - *m_mic_fab_vars[s_qv_idx(v)], - 0, 0, 1, amrex::IntVect::TheZeroVector() ); - return; - } - } - { - std::string varname = "qc_" + amrex::getEnumNameString(m_species[v]); - if (a_name == varname) { - amrex::MultiFab::Copy( a_mf, - *m_mic_fab_vars[s_qc_idx(v)], - 0, 0, 1, amrex::IntVect::TheZeroVector() ); - return; - } - } - { - std::string varname = "qt_" + amrex::getEnumNameString(m_species[v]); - if (a_name == varname) { - amrex::MultiFab::Copy( a_mf, - *m_mic_fab_vars[s_qt_idx(v)], - 0, 0, 1, amrex::IntVect::TheZeroVector() ); - return; - } - } - { - std::string varname = "sat_ratio_" + amrex::getEnumNameString(m_species[v]); - if (a_name == varname) { - amrex::MultiFab::Copy( a_mf, - *m_mic_fab_vars[s_sr_idx(v)], - 0, 0, 1, amrex::IntVect::TheZeroVector() ); - return; - } + + // Helper to copy MultiFab if name matches + auto try_copy = [&](const std::string& prefix, int idx) -> bool { + if (a_name == prefix) { + amrex::MultiFab::Copy(a_mf, *m_mic_fab_vars[idx], + 0, 0, 1, amrex::IntVect::TheZeroVector()); + return true; } - { - std::string varname = "accum_" + amrex::getEnumNameString(m_species[v]); - if (a_name == varname) { - amrex::MultiFab::Copy( a_mf, - *m_mic_fab_vars[s_accum_idx(v)], - 0, 0, 1, amrex::IntVect::TheZeroVector() ); - return; - } + return false; + }; + + // Species variables + for (int v = 1; v < m_num_species; v++) { + std::string sp_name = amrex::getEnumNameString(m_species[v]); + if (try_copy("qv_" + sp_name, s_qv_idx(v)) || + try_copy("qc_" + sp_name, s_qc_idx(v)) || + try_copy("qt_" + sp_name, s_qt_idx(v)) || + try_copy("sat_ratio_" + sp_name, s_sr_idx(v)) || + try_copy("accum_" + sp_name, s_accum_idx(v))) { + return; } } + + // Aerosol variables for (int v = 0; v < m_num_aerosols; v++) { - { - std::string varname = "accum_" + amrex::getEnumNameString(m_aerosols[v]); - if (a_name == varname) { - amrex::MultiFab::Copy( a_mf, - *m_mic_fab_vars[a_accum_idx(m_num_species,v)], - 0, 0, 1, amrex::IntVect::TheZeroVector() ); - return; - } - } + std::string ae_name = amrex::getEnumNameString(m_aerosols[v]); + if (try_copy("accum_" + ae_name, a_accum_idx(m_num_species, v))) { return; } } + amrex::Abort("SuperDropletsMoist::GetPlotVar() called with invalid name"); - return; } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE diff --git a/Source/Microphysics/SuperDropletsMoist/ERF_SuperDropletsMoistUtils.cpp b/Source/Microphysics/SuperDropletsMoist/ERF_SuperDropletsMoistUtils.cpp index d80ab5bf3..c2e6163bc 100644 --- a/Source/Microphysics/SuperDropletsMoist/ERF_SuperDropletsMoistUtils.cpp +++ b/Source/Microphysics/SuperDropletsMoist/ERF_SuperDropletsMoistUtils.cpp @@ -326,14 +326,12 @@ void SuperDropletsMoist::ratioToDensity (MultiFab& a_var, void SuperDropletsMoist::computeQcQrWater () { BL_PROFILE("SuperDropletsMoist::computeQcQrWater()"); - m_super_droplets->speciesMassDensity( *(m_mic_fab_vars[MicVar_SD::q_c]), - m_idx_w, - 0, - m_r_rain ); - m_super_droplets->speciesMassDensity( *(m_mic_fab_vars[MicVar_SD::q_r]), - m_idx_w, - m_r_rain, - one ); + m_super_droplets->cloudRainDensity( *(m_mic_fab_vars[MicVar_SD::q_c]), + 0, + m_r_rain ); + m_super_droplets->cloudRainDensity( *(m_mic_fab_vars[MicVar_SD::q_r]), + m_r_rain, + one ); if (m_dimensionality == SDMSimulationDim::one_d_z) { for ( MFIter mfi(*m_mic_fab_vars[MicVar_SD::q_c]); mfi.isValid(); ++mfi) { diff --git a/Source/Particles/ERFPC.H b/Source/Particles/ERFPC.H index 318d530f1..02917beeb 100644 --- a/Source/Particles/ERFPC.H +++ b/Source/Particles/ERFPC.H @@ -333,6 +333,17 @@ class ERFPC : public amrex::ParticleContainer< ERFParticlesRealIdxAoS::ncomps, // the following functions should ideally be private or protected, but need to be // public due to CUDA extended lambda capture rules + /*! \brief Helper for ParticleToMesh operations with non-uniform z support + * + * Generic CIC deposition from particles to a MultiFab. The caller provides + * a device lambda `value_func(ptd, i)` that returns the per-particle scalar + * to deposit (e.g. mass, number, flux). Handles both uniform and non-uniform + * vertical grids. Implementation in ERFPCParticleToMesh.H. + */ + template + void ERFPCParticleToMesh(amrex::MultiFab& a_mf, int a_lev, int a_comp, + ValueFunc&& value_func) const; + /*! Default particle initialization */ void initializeParticlesUniformDistributionInBox (const std::unique_ptr& a_ptr, const amrex::RealBox& particle_box); diff --git a/Source/Particles/ERFPCParticleToMesh.H b/Source/Particles/ERFPCParticleToMesh.H new file mode 100644 index 000000000..4b43a0524 --- /dev/null +++ b/Source/Particles/ERFPCParticleToMesh.H @@ -0,0 +1,117 @@ +/*! \file ERFPCParticleToMesh.H + * \brief Template implementation of ERFPC::ERFPCParticleToMesh + * + * Separated from ERFPC.H to keep the header lightweight. + * Include this in .cpp files that call ERFPCParticleToMesh(). + */ + +#ifndef ERFPC_PARTICLE_TO_MESH_H_ +#define ERFPC_PARTICLE_TO_MESH_H_ + +#include +#include + +template +void ERFPC::ERFPCParticleToMesh(amrex::MultiFab& a_mf, int a_lev, int a_comp, + ValueFunc&& value_func) const +{ + AMREX_ASSERT(OK()); + AMREX_ASSERT(numParticlesOutOfRange(*this, 0) == 0); + + const auto& geom = Geom(a_lev); + const auto plo = geom.ProbLoArray(); + const auto dxi = geom.InvCellSizeArray(); + + const amrex::Real* zlevels = m_zlevels_d.empty() ? nullptr : m_zlevels_d.data(); + const int nz_levels = static_cast(m_zlevels_d.size()); + + a_mf.setVal(0.0); + + amrex::ParticleToMesh(*this, a_mf, a_lev, + [=] AMREX_GPU_DEVICE (const typename ERFPC::ParticleTileType::ConstParticleTileDataType& ptd, + int i, amrex::Array4 const& rho) + { + auto p = ptd.m_aos[i]; + auto particle_value = value_func(ptd, i); + + if (zlevels) { + amrex::Real lx = (p.pos(0) - plo[0]) * dxi[0] + amrex::Real(0.5); + int ix = static_cast(amrex::Math::floor(lx)) - 1; + amrex::Real wx1 = lx - static_cast(ix + 1); + amrex::Real wx0 = amrex::Real(1.0) - wx1; + + amrex::Real ly = (p.pos(1) - plo[1]) * dxi[1] + amrex::Real(0.5); + int iy = static_cast(amrex::Math::floor(ly)) - 1; + amrex::Real wy1 = ly - static_cast(iy + 1); + amrex::Real wy0 = amrex::Real(1.0) - wy1; + + int kz = p.idata(ERFParticlesIntIdxAoS::k); + int nz = nz_levels - 1; + + amrex::Real dz_k = zlevels[kz+1] - zlevels[kz]; + amrex::Real z_center_k = amrex::Real(0.5) * (zlevels[kz] + zlevels[kz+1]); + + int kz_lo, kz_hi; + amrex::Real wz_lo, wz_hi; + amrex::Real dz_lo, dz_hi; + + if (p.pos(2) >= z_center_k) { + kz_lo = kz; + kz_hi = kz + 1; + dz_lo = dz_k; + if (kz_hi < nz) { + dz_hi = zlevels[kz_hi+1] - zlevels[kz_hi]; + amrex::Real z_center_hi = amrex::Real(0.5) * (zlevels[kz_hi] + zlevels[kz_hi+1]); + amrex::Real span = z_center_hi - z_center_k; + wz_hi = (p.pos(2) - z_center_k) / span; + wz_lo = amrex::Real(1.0) - wz_hi; + } else { + kz_hi = kz; + dz_hi = dz_k; + wz_lo = amrex::Real(1.0); + wz_hi = amrex::Real(0.0); + } + } else { + kz_lo = kz - 1; + kz_hi = kz; + dz_hi = dz_k; + if (kz_lo >= 0) { + dz_lo = zlevels[kz_lo+1] - zlevels[kz_lo]; + amrex::Real z_center_lo = amrex::Real(0.5) * (zlevels[kz_lo] + zlevels[kz_lo+1]); + amrex::Real span = z_center_k - z_center_lo; + wz_lo = (z_center_k - p.pos(2)) / span; + wz_hi = amrex::Real(1.0) - wz_lo; + } else { + kz_lo = kz; + dz_lo = dz_k; + wz_lo = amrex::Real(0.0); + wz_hi = amrex::Real(1.0); + } + } + + amrex::Real inv_vol_lo = dxi[0] * dxi[1] / dz_lo; + amrex::Real inv_vol_hi = dxi[0] * dxi[1] / dz_hi; + + amrex::Real pval_lo = particle_value * wz_lo * inv_vol_lo; + amrex::Real pval_hi = particle_value * wz_hi * inv_vol_hi; + + amrex::Gpu::Atomic::AddNoRet(&rho(ix , iy , kz_lo, a_comp), wx0*wy0*pval_lo); + amrex::Gpu::Atomic::AddNoRet(&rho(ix+1, iy , kz_lo, a_comp), wx1*wy0*pval_lo); + amrex::Gpu::Atomic::AddNoRet(&rho(ix , iy+1, kz_lo, a_comp), wx0*wy1*pval_lo); + amrex::Gpu::Atomic::AddNoRet(&rho(ix+1, iy+1, kz_lo, a_comp), wx1*wy1*pval_lo); + amrex::Gpu::Atomic::AddNoRet(&rho(ix , iy , kz_hi, a_comp), wx0*wy0*pval_hi); + amrex::Gpu::Atomic::AddNoRet(&rho(ix+1, iy , kz_hi, a_comp), wx1*wy0*pval_hi); + amrex::Gpu::Atomic::AddNoRet(&rho(ix , iy+1, kz_hi, a_comp), wx0*wy1*pval_hi); + amrex::Gpu::Atomic::AddNoRet(&rho(ix+1, iy+1, kz_hi, a_comp), wx1*wy1*pval_hi); + } else { + amrex::ParticleInterpolator::Linear interp(p, plo, dxi); + amrex::ParticleReal inv_cell_volume = dxi[0] * dxi[1] * dxi[2]; + interp.ParticleToMesh(p, rho, 0, a_comp, 1, + [=] AMREX_GPU_DEVICE (const ERFPC::ParticleType&, int) { + return particle_value * inv_cell_volume; + }); + } + }); +} + +#endif diff --git a/Source/Particles/ERFPCUtils.cpp b/Source/Particles/ERFPCUtils.cpp index 76ef351d9..3aa974a31 100644 --- a/Source/Particles/ERFPCUtils.cpp +++ b/Source/Particles/ERFPCUtils.cpp @@ -1,6 +1,6 @@ #include +#include #include -#include #include #ifdef ERF_USE_PARTICLES @@ -21,32 +21,10 @@ void ERFPC::massDensity ( MultiFab& a_mf, const int& a_comp ) const { BL_PROFILE("ERFPC::massDensity()"); - - AMREX_ASSERT(OK()); - AMREX_ASSERT(numParticlesOutOfRange(*this, 0) == 0); - - const auto& geom = Geom(a_lev); - const auto plo = geom.ProbLoArray(); - const auto dxi = geom.InvCellSizeArray(); - - const Real inv_cell_volume = dxi[0]*dxi[1]*dxi[2]; - a_mf.setVal(0.0); - - ParticleToMesh( *this, a_mf, a_lev, - [=] AMREX_GPU_DEVICE ( const ERFPC::ParticleTileType::ConstParticleTileDataType& ptr, - int i, Array4 const& rho) - { - auto p = ptr.m_aos[i]; - ParticleInterpolator::Linear interp(p, plo, dxi); - interp.ParticleToMesh ( p, rho, 0, a_comp, 1, - [=] AMREX_GPU_DEVICE ( const ERFPC::ParticleType&, int) - { - auto mass = ptr.m_rdata[ERFParticlesRealIdxSoA::mass][i]; - return mass*inv_cell_volume; - }); + ERFPCParticleToMesh(a_mf, a_lev, a_comp, + [=] AMREX_GPU_DEVICE (const ERFPC::ParticleTileType::ConstParticleTileDataType& ptd, int i) { + return ptd.m_rdata[ERFParticlesRealIdxSoA::mass][i]; }); - - return; } /*! \brief Fix k-indices for all particles after AMR regrid */ @@ -507,4 +485,6 @@ void ERFPC::ExtractAndRouteOORParticles ( int } } + #endif + diff --git a/Source/Particles/ERF_SDInitialization.H b/Source/Particles/ERF_SDInitialization.H index 390e9f436..627632259 100644 --- a/Source/Particles/ERF_SDInitialization.H +++ b/Source/Particles/ERF_SDInitialization.H @@ -10,6 +10,8 @@ #include #include #include +#include +#include #include "ERF_Constants.H" #include "ERF_MaterialProperties.H" @@ -20,6 +22,14 @@ AMREX_ENUM(SDInitShape, null ); +/*! \brief Distribution type enum for GPU compatibility */ +AMREX_ENUM(SDDistributionType, + mass_constant, + mass_exponential, + radius_log_normal, + radius_lognormal_autorange +); + /*! \brief List of super-droplet initializations */ namespace SupDropInit { @@ -28,17 +38,113 @@ namespace SupDropInit /*! Maximum number of aerosols */ const int num_aerosols_max = 8; - - const std::string attrib_init_const = "mass_constant"; - const std::string attrib_init_exp = "mass_exponential"; - const std::string attrib_init_lnr = "radius_log_normal"; - const std::string attrib_init_lnr_auto = "radius_lognormal_autorange"; } AMREX_ENUM(SDMultiplicityType, constant, sampled ); +/*! \brief GPU-compatible structure holding distribution parameters */ +struct SDDistributionParams { + SDDistributionType dist_type; /*!< Type of distribution */ + amrex::Real mass_min; /*!< Minimum mass */ + amrex::Real mass_max; /*!< Maximum mass */ + amrex::Real mass_mean; /*!< Mean mass */ + amrex::Real radius_min; /*!< Minimum radius */ + amrex::Real radius_max; /*!< Maximum radius */ + amrex::Real radius_mean; /*!< Mean radius (mu for log-normal) */ + amrex::Real radius_gstd; /*!< Geometric std dev for radius */ + amrex::Real density; /*!< Material density */ + amrex::Real numdens; /*!< Number density for multiplicity */ + amrex::Real cell_volume; /*!< Cell volume for multiplicity scaling */ + // Pre-computed values for log-normal CDF inversion (truncated distribution) + amrex::Real cdf_min; /*!< CDF at radius_min */ + amrex::Real cdf_max; /*!< CDF at radius_max */ + amrex::Real sigma; /*!< log(radius_gstd) */ + amrex::Real lnrng; /*!< log(radius_max) - log(radius_min) */ + amrex::Real lnmin; /*!< log(radius_min) or log(mass_min) */ + amrex::Real delta; /*!< mass_mean - mass_min for exponential */ + int sampled_mult; /*!< 1 if sampled multiplicity, 0 if constant */ +}; + +/*! \brief Inverse error function approximation for GPU + * \param[in] x Input value in (-1, 1) + * \return Approximate inverse error function value + */ +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +amrex::Real SD_erfinv_gpu(const amrex::Real x) { + amrex::Real a = 0.147; + amrex::Real eps = std::numeric_limits::epsilon(); + amrex::Real term = std::log(1 - x * x + eps); + amrex::Real p1 = 2 / (PI * a) + term / 2.0; + amrex::Real p2 = term / a; + amrex::Real sign = (x >= 0) ? 1.0 : -1.0; + return sign * std::sqrt(std::sqrt(p1 * p1 - p2) - p1); +} + +/*! \brief Generate mass from distribution parameters on GPU + * \param[in] params Distribution parameters + * \param[in] engine Random engine + * \param[out] mult_contribution Contribution to multiplicity (added, not set) + * \return Generated mass value + */ +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +amrex::Real SD_sample_mass_gpu(const SDDistributionParams& params, + const amrex::RandomEngine& engine, + amrex::Real& mult_contribution) +{ + using amrex::Random; + amrex::Real mass = 0.0; + amrex::Real u = Random(engine); + + switch (params.dist_type) { + case SDDistributionType::mass_constant: + mass = params.mass_mean; + mult_contribution = u; // Random [0,1] for later rescaling + break; + + case SDDistributionType::mass_exponential: { + if (params.sampled_mult) { + // Sampled multiplicity: sample uniformly in log-space, weight by exponential PDF + amrex::Real lnval = params.lnmin + u * params.lnrng; + mass = std::exp(lnval); + mult_contribution = (params.numdens * params.cell_volume) * std::exp(-mass / params.delta); + } else { + // Constant multiplicity: true exponential distribution via inverse transform + // If U ~ Uniform(0,1), then -delta*ln(U) ~ Exponential(1/delta) + mass = -params.delta * std::log(u) + params.mass_min; + mult_contribution = 0.0; // Not used for constant multiplicity + } + break; + } + + case SDDistributionType::radius_log_normal: + case SDDistributionType::radius_lognormal_autorange: { + amrex::Real dry_r; + if (params.sampled_mult) { + // Sampled multiplicity: sample uniformly in log-space, weight by log-normal PDF + amrex::Real lnval = params.lnmin + u * params.lnrng; + dry_r = std::exp(lnval); + // Log-normal PDF contribution for multiplicity weighting + amrex::Real term = std::exp(-std::log(dry_r/params.radius_mean)*std::log(dry_r/params.radius_mean) + /(2.0*params.sigma*params.sigma)); + mult_contribution = (params.numdens * params.cell_volume) / (params.sigma * std::sqrt(2.0*PI)) * term; + } else { + // Constant multiplicity: true log-normal via inverse CDF + // Map u from [0,1] to [cdf_min, cdf_max] for truncated distribution + amrex::Real u_trunc = params.cdf_min + u * (params.cdf_max - params.cdf_min); + // Inverse CDF: r = mu * exp(sigma * sqrt(2) * erfinv(2*u - 1)) + amrex::Real z = SD_erfinv_gpu(2.0 * u_trunc - 1.0) * std::sqrt(2.0); + dry_r = params.radius_mean * std::exp(params.sigma * z); + mult_contribution = 0.0; // Not used for constant multiplicity + } + mass = four_thirds_pi * dry_r * dry_r * dry_r * params.density; + break; + } + } + return mass; +} + /*! \brief Super-droplets initial properties */ class SDInitProperties { @@ -72,7 +178,7 @@ class SDInitProperties /*! Standard deviation of species dry radius */ std::vector m_radius_species_geom_std; /*! Initial distribution type for species */ - std::vector m_species_init_type; + std::vector m_species_init_type; /*! Minimum aerosol mass */ std::vector m_mass_aerosol_min; @@ -89,7 +195,7 @@ class SDInitProperties /*! Standard deviation of aerosol dry radius */ std::vector m_radius_aerosol_geom_std; /*! Initial distribution type for aerosol */ - std::vector m_aerosol_init_type; + std::vector m_aerosol_init_type; /*! Number of species */ int m_num_species = 0; @@ -151,16 +257,16 @@ class SDInitProperties * \param[in,out] a_rng Random number generator */ void getDistribution( amrex::Vector& a_mass, - const int a_np, - const amrex::Real a_density, - const std::string& a_init_type, - const amrex::Real a_mass_min, - const amrex::Real a_mass_max, - const amrex::Real a_mass_mean, - const amrex::Real a_radius_min, - const amrex::Real a_radius_max, - const amrex::Real a_radius_mean, - const amrex::Real a_radius_gstd, + int a_np, + amrex::Real a_density, + SDDistributionType a_init_type, + amrex::Real a_mass_min, + amrex::Real a_mass_max, + amrex::Real a_mass_mean, + amrex::Real a_radius_min, + amrex::Real a_radius_max, + amrex::Real a_radius_mean, + amrex::Real a_radius_gstd, std::mt19937& a_rng ) const; /*! \brief Get a distribution with sampled multiplicity @@ -181,17 +287,17 @@ class SDInitProperties */ void getDistribution( amrex::Vector& a_mass, amrex::Vector& a_mult, - const amrex::Real a_dV, - const int a_np, - const amrex::Real a_density, - const std::string& a_init_type, - const amrex::Real a_mass_min, - const amrex::Real a_mass_max, - const amrex::Real a_mass_mean, - const amrex::Real a_radius_min, - const amrex::Real a_radius_max, - const amrex::Real a_radius_mean, - const amrex::Real a_radius_gstd, + amrex::Real a_dV, + int a_np, + amrex::Real a_density, + SDDistributionType a_init_type, + amrex::Real a_mass_min, + amrex::Real a_mass_max, + amrex::Real a_mass_mean, + amrex::Real a_radius_min, + amrex::Real a_radius_max, + amrex::Real a_radius_mean, + amrex::Real a_radius_gstd, std::mt19937& a_rng ) const; /*! \brief Compute the aerosol mass distribution @@ -314,6 +420,87 @@ class SDInitProperties a_rng); } + /*! \brief Get GPU-compatible distribution parameters for a species + * \param[in] a_idx Species index + * \param[in] a_density Density of the species material + * \param[in] a_cell_volume Cell volume for multiplicity scaling + * \param[in] a_sampled_mult Whether using sampled multiplicity mode + * \return SDDistributionParams structure for GPU use + */ + [[nodiscard]] SDDistributionParams getSpeciesDistParams( + int a_idx, + amrex::Real a_density, + amrex::Real a_cell_volume, + bool a_sampled_mult) const + { + return makeDistributionParams( + m_species_init_type[a_idx], + m_mass_species_min[a_idx], + m_mass_species_max[a_idx], + m_mass_species_mean[a_idx], + m_radius_species_min[a_idx], + m_radius_species_max[a_idx], + m_radius_species_mean[a_idx], + m_radius_species_geom_std[a_idx], + a_density, + a_cell_volume, + a_sampled_mult); + } + + /*! \brief Get GPU-compatible distribution parameters for an aerosol + * \param[in] a_idx Aerosol species index + * \param[in] a_density Density of the aerosol material + * \param[in] a_cell_volume Cell volume for multiplicity scaling + * \param[in] a_sampled_mult Whether using sampled multiplicity mode + * \return SDDistributionParams structure for GPU use + */ + [[nodiscard]] SDDistributionParams getAerosolDistParams( + int a_idx, + amrex::Real a_density, + amrex::Real a_cell_volume, + bool a_sampled_mult) const + { + return makeDistributionParams( + m_aerosol_init_type[a_idx], + m_mass_aerosol_min[a_idx], + m_mass_aerosol_max[a_idx], + m_mass_aerosol_mean[a_idx], + m_radius_aerosol_min[a_idx], + m_radius_aerosol_max[a_idx], + m_radius_aerosol_mean[a_idx], + m_radius_aerosol_geom_std[a_idx], + a_density, + a_cell_volume, + a_sampled_mult); + } + + /*! \brief Create GPU-compatible distribution parameters structure + * \param[in] a_init_type Distribution type + * \param[in] a_mass_min Minimum mass + * \param[in] a_mass_max Maximum mass + * \param[in] a_mass_mean Mean mass + * \param[in] a_radius_min Minimum radius + * \param[in] a_radius_max Maximum radius + * \param[in] a_radius_mean Mean radius + * \param[in] a_radius_gstd Geometric standard deviation + * \param[in] a_density Material density + * \param[in] a_cell_volume Cell volume + * \param[in] a_sampled_mult Whether using sampled multiplicity mode + * \return SDDistributionParams structure + */ + [[nodiscard]] SDDistributionParams makeDistributionParams( + SDDistributionType a_init_type, + amrex::Real a_mass_min, + amrex::Real a_mass_max, + amrex::Real a_mass_mean, + amrex::Real a_radius_min, + amrex::Real a_radius_max, + amrex::Real a_radius_mean, + amrex::Real a_radius_gstd, + amrex::Real a_density, + amrex::Real a_cell_volume, + bool a_sampled_mult) const; + /*! \brief Determine whether multiplicity is sampled or constant * \return True if multiplicity is sampled, false if constant */ @@ -336,7 +523,7 @@ class SDInitProperties vol = m_particle_domain.volume(); } else if (m_type == SDInitShape::bubble) { const auto& radius = m_particle_domain.hi(); - vol = (amrex::Real(4.0)/three)*PI*radius[0]*radius[1]*radius[2]; + vol = four_thirds_pi*radius[0]*radius[1]*radius[2]; } return vol; } @@ -372,9 +559,9 @@ class SDInjection : public SDInitProperties /*! \brief Update time-dependent quantities for particle injection * * This function updates: - * one Number density based on injection rate and timestep - * two Super-droplet number density based on SD injection rate - * three Position of the injection domain based on domain velocity + * 1. Number density based on injection rate and timestep + * 2. Super-droplet number density based on SD injection rate + * 3. Position of the injection domain based on domain velocity * * \param[in] a_dt Timestep size in seconds */ @@ -434,8 +621,8 @@ class SDInjection : public SDInitProperties /*! \brief Compute number of super-droplets to inject per grid cell * * Determines the number of super-droplets to inject based on either: - * one The super-droplet number density and cell volume, or - * two A constant number per cell (m_ppc) if number density is not set + * 1. The super-droplet number density and cell volume, or + * 2. A constant number per cell (m_ppc) if number density is not set * * \param[in] a_dv Volume of the grid cell * \return Number of super-droplets to inject in this cell diff --git a/Source/Particles/ERF_SDInitialization.cpp b/Source/Particles/ERF_SDInitialization.cpp index 2216a187e..ac82e9a00 100644 --- a/Source/Particles/ERF_SDInitialization.cpp +++ b/Source/Particles/ERF_SDInitialization.cpp @@ -29,13 +29,22 @@ void SDInitProperties::setDefaults ( const amrex::Geometry& a_geom, for (int i = 0; i < m_num_species; i++) { // default values - m_species_init_type[i] = SupDropInit::attrib_init_const; - m_mass_species_min[i] = amrex::Real(4.1887902e-42); - m_mass_species_max[i] = amrex::Real(4.1887902e-42); - m_mass_species_mean[i] = amrex::Real(4.1887902e-42); - m_radius_species_min[i] = amrex::Real(1.0e-15); - m_radius_species_max[i] = amrex::Real(1.0e-15); - m_radius_species_mean[i] = amrex::Real(1.0e-15); + m_species_init_type[i] = SDDistributionType::mass_constant; + if (a_species_mat[i]->m_is_water) { + m_mass_species_min[i] = amrex::Real(4.1887902e-42); + m_mass_species_max[i] = amrex::Real(4.1887902e-42); + m_mass_species_mean[i] = amrex::Real(4.1887902e-42); + m_radius_species_min[i] = amrex::Real(1.0e-15); + m_radius_species_max[i] = amrex::Real(1.0e-15); + m_radius_species_mean[i] = amrex::Real(1.0e-15); + } else { + m_mass_species_min[i] = zero; + m_mass_species_max[i] = zero; + m_mass_species_mean[i] = zero; + m_radius_species_min[i] = zero; + m_radius_species_max[i] = zero; + m_radius_species_mean[i] = zero; + } m_radius_species_geom_std[i] = two; } @@ -51,7 +60,7 @@ void SDInitProperties::setDefaults ( const amrex::Geometry& a_geom, for (int i = 0; i < m_num_aerosols; i++) { // default values - m_aerosol_init_type[i] = SupDropInit::attrib_init_const; + m_aerosol_init_type[i] = SDDistributionType::mass_constant; m_mass_aerosol_min[i] = zero; m_mass_aerosol_max[i] = zero; m_mass_aerosol_mean[i] = zero; @@ -103,14 +112,6 @@ void SDInitProperties::readInputs ( const std::string& a_prefix, m_particle_domain.setHi(m_init_particle_p2); } - // Backward compatibility - for (int i = 0; i < m_num_species; i++) { - pp.query(std::string(a_key+"condensate_distribution_type").c_str(), m_species_init_type[i]); - pp.query(std::string(a_key+"condensate_mass_min").c_str(), m_mass_species_min[i]); - pp.query(std::string(a_key+"condensate_mass_mean").c_str(), m_mass_species_mean[i]); - pp.query(std::string(a_key+"condensate_min_radius").c_str(), m_radius_species_min[i]); - pp.query(std::string(a_key+"condensate_max_radius").c_str(), m_radius_species_max[i]); - } for (int i = 0; i < m_num_species; i++) { { std::string key = a_key+"species_distribution_type_"+getEnumNameString(a_species_mat[i]->m_name); @@ -260,18 +261,18 @@ void SDInitProperties::printParameters ( const MatVec& a_species_mat, for (unsigned long i=0; i < a_species_mat.size(); i++) { Print() << " " << getEnumNameString(a_species_mat[i]->m_name) - << " (distribution: " << m_species_init_type[i]; - if (m_species_init_type[i] == SupDropInit::attrib_init_const) { + << " (distribution: " << getEnumNameString(m_species_init_type[i]); + if (m_species_init_type[i] == SDDistributionType::mass_constant) { Print() << ", value=" << m_mass_species_mean[i]; - } else if (m_species_init_type[i] == SupDropInit::attrib_init_exp) { + } else if (m_species_init_type[i] == SDDistributionType::mass_exponential) { Print() << ", min=" << m_mass_species_min[i] << ", mean=" << m_mass_species_mean[i] << ", max=" << m_mass_species_max[i]; - AMREX_ALWAYS_ASSERT(m_mass_species_min[i] > zero); + AMREX_ALWAYS_ASSERT(m_mass_species_min[i] >= zero); AMREX_ALWAYS_ASSERT(m_mass_species_max[i] >= m_mass_species_min[i]); AMREX_ALWAYS_ASSERT( (m_mass_species_mean[i] >= m_mass_species_min[i]) && (m_mass_species_mean[i] <= m_mass_species_max[i]) ); - } else if (m_species_init_type[i] == SupDropInit::attrib_init_lnr) { + } else if (m_species_init_type[i] == SDDistributionType::radius_log_normal) { Print() << ", min=" << m_radius_species_min[i] << ", max=" << m_radius_species_max[i] << ", mean=" << m_radius_species_mean[i] @@ -280,7 +281,7 @@ void SDInitProperties::printParameters ( const MatVec& a_species_mat, AMREX_ALWAYS_ASSERT(m_radius_species_max[i] >= m_radius_species_min[i]); AMREX_ALWAYS_ASSERT( (m_radius_species_mean[i] >= m_radius_species_min[i]) && (m_radius_species_mean[i] <= m_radius_species_max[i]) ); - } else if (m_species_init_type[i] == SupDropInit::attrib_init_lnr_auto) { + } else if (m_species_init_type[i] == SDDistributionType::radius_lognormal_autorange) { Print() << ", mean=" << m_radius_species_mean[i] << ", std=" << m_radius_species_geom_std[i]; AMREX_ALWAYS_ASSERT(m_radius_species_mean[i] > zero); @@ -293,11 +294,11 @@ void SDInitProperties::printParameters ( const MatVec& a_species_mat, for (unsigned long i=0; i < a_aerosol_mat.size(); i++) { Print() << " " << getEnumNameString(a_aerosol_mat[i]->m_name) - << " (distribution: " << m_aerosol_init_type[i]; - if (m_aerosol_init_type[i] == SupDropInit::attrib_init_const) { + << " (distribution: " << getEnumNameString(m_aerosol_init_type[i]); + if (m_aerosol_init_type[i] == SDDistributionType::mass_constant) { Print() << ", value=" << m_mass_aerosol_mean[i]; AMREX_ALWAYS_ASSERT(m_mass_aerosol_mean[i] > zero); - } else if (m_aerosol_init_type[i] == SupDropInit::attrib_init_exp) { + } else if (m_aerosol_init_type[i] == SDDistributionType::mass_exponential) { Print() << ", min=" << m_mass_aerosol_min[i] << ", mean=" << m_mass_aerosol_mean[i] << ", max=" << m_mass_aerosol_max[i]; @@ -305,7 +306,7 @@ void SDInitProperties::printParameters ( const MatVec& a_species_mat, AMREX_ALWAYS_ASSERT(m_mass_aerosol_max[i] >= m_mass_aerosol_min[i]); AMREX_ALWAYS_ASSERT( (m_mass_aerosol_mean[i] >= m_mass_aerosol_min[i]) && (m_mass_aerosol_mean[i] <= m_mass_aerosol_max[i]) ); - } else if (m_aerosol_init_type[i] == SupDropInit::attrib_init_lnr) { + } else if (m_aerosol_init_type[i] == SDDistributionType::radius_log_normal) { Print() << ", min=" << m_radius_aerosol_min[i] << ", max=" << m_radius_aerosol_max[i] << ", mean=" << m_radius_aerosol_mean[i] @@ -314,7 +315,7 @@ void SDInitProperties::printParameters ( const MatVec& a_species_mat, AMREX_ALWAYS_ASSERT(m_radius_aerosol_max[i] >= m_radius_aerosol_min[i]); AMREX_ALWAYS_ASSERT( (m_radius_aerosol_mean[i] >= m_radius_aerosol_min[i]) && (m_radius_aerosol_mean[i] <= m_radius_aerosol_max[i]) ); - } else if (m_aerosol_init_type[i] == SupDropInit::attrib_init_lnr_auto) { + } else if (m_aerosol_init_type[i] == SDDistributionType::radius_lognormal_autorange) { Print() << ", mean=" << m_radius_aerosol_mean[i] << ", std=" << m_radius_aerosol_geom_std[i]; AMREX_ALWAYS_ASSERT(m_radius_aerosol_mean[i] > zero); @@ -348,30 +349,30 @@ void SDInjection::printParameters ( const MatVec& a_species_mat, } void SDInitProperties::getDistribution ( amrex::Vector& a_mass, - const int a_np, - const amrex::Real a_density, - const std::string& a_init_type, - const amrex::Real a_mass_min, - const amrex::Real /*a_mass_max*/, - const amrex::Real a_mass_mean, - const amrex::Real a_radius_min, - const amrex::Real a_radius_max, - const amrex::Real a_radius_mean, - const amrex::Real a_radius_gstd, + int a_np, + amrex::Real a_density, + SDDistributionType a_init_type, + amrex::Real a_mass_min, + amrex::Real /*a_mass_max*/, + amrex::Real a_mass_mean, + amrex::Real a_radius_min, + amrex::Real a_radius_max, + amrex::Real a_radius_mean, + amrex::Real a_radius_gstd, std::mt19937& a_rng ) const { a_mass.resize(a_np); - if (a_init_type == SupDropInit::attrib_init_const) { + if (a_init_type == SDDistributionType::mass_constant) { for (int n = 0; n < a_np; n++) { a_mass[n] = a_mass_mean; } - } else if (a_init_type == SupDropInit::attrib_init_exp) { + } else if (a_init_type == SDDistributionType::mass_exponential) { auto delta = a_mass_mean - a_mass_min; std::exponential_distribution ed(one/delta); for (int n = 0; n < a_np; n++) { a_mass[n] = ed(a_rng) + a_mass_min; } - } else if (a_init_type == SupDropInit::attrib_init_lnr) { + } else if (a_init_type == SDDistributionType::radius_log_normal) { std::normal_distribution nrd(std::log(a_radius_mean), std::log(a_radius_gstd)); for (int n = 0; n < a_np; n++) { @@ -382,7 +383,7 @@ void SDInitProperties::getDistribution ( amrex::Vector& a_mass, count++; if (count > 100) { break; } } - a_mass[n] = (amrex::Real(4.0)/three) * PI + a_mass[n] = four_thirds_pi * dry_r * dry_r * dry_r * a_density; } @@ -403,28 +404,28 @@ static amrex::Real SD_erfinv(const amrex::Real x) { void SDInitProperties::getDistribution ( amrex::Vector& a_mass, amrex::Vector& a_mult, - const amrex::Real a_dV, - const int a_np, - const amrex::Real a_density, - const std::string& a_init_type, - const amrex::Real a_mass_min, - const amrex::Real a_mass_max, - const amrex::Real a_mass_mean, - const amrex::Real a_radius_min, - const amrex::Real a_radius_max, - const amrex::Real a_radius_mean, - const amrex::Real a_radius_gstd, + amrex::Real a_dV, + int a_np, + amrex::Real a_density, + SDDistributionType a_init_type, + amrex::Real a_mass_min, + amrex::Real a_mass_max, + amrex::Real a_mass_mean, + amrex::Real a_radius_min, + amrex::Real a_radius_max, + amrex::Real a_radius_mean, + amrex::Real a_radius_gstd, std::mt19937& a_rng ) const { a_mass.resize(a_np); - AMREX_ALWAYS_ASSERT(a_mult.size() == a_np); - if (a_init_type == SupDropInit::attrib_init_const) { + AMREX_ALWAYS_ASSERT(static_cast(a_mult.size()) == a_np); + if (a_init_type == SDDistributionType::mass_constant) { std::uniform_real_distribution<> urd(zero, one); for (int n = 0; n < a_np; n++) { a_mass[n] = a_mass_mean; a_mult[n] += urd(a_rng); // initially this will be a non-integer; later we will rescale to an integer. } - } else if (a_init_type == SupDropInit::attrib_init_exp) { + } else if (a_init_type == SDDistributionType::mass_exponential) { std::uniform_real_distribution<> urd(zero, one); auto delta = a_mass_mean - a_mass_min; auto lnrng = std::log(a_mass_max) - std::log(a_mass_min); @@ -434,7 +435,7 @@ void SDInitProperties::getDistribution ( amrex::Vector& a_mass, a_mass[n] = std::exp(tmp); a_mult[n] += (m_numdens * a_dV) * std::exp(-a_mass[n] / delta); } - } else if (a_init_type == SupDropInit::attrib_init_lnr) { + } else if (a_init_type == SDDistributionType::radius_log_normal) { std::uniform_real_distribution<> urd(zero, one); auto sigma = std::log(a_radius_gstd); auto mu = a_radius_mean; @@ -443,11 +444,11 @@ void SDInitProperties::getDistribution ( amrex::Vector& a_mass, for (int n = 0; n < a_np; n++) { auto tmp = lnmin + urd(a_rng) * lnrng; auto dry_r = std::exp(tmp); - a_mass[n] = (amrex::Real(4.0)/three) * PI * dry_r * dry_r * dry_r * a_density; + a_mass[n] = four_thirds_pi * dry_r * dry_r * dry_r * a_density; auto term = std::exp(-std::log(dry_r/mu)*std::log(dry_r/mu)/(two*sigma*sigma)); a_mult[n] += ( m_numdens * a_dV ) / (sigma*std::sqrt(2*PI)) * term; } - } else if (a_init_type == SupDropInit::attrib_init_lnr_auto) { + } else if (a_init_type == SDDistributionType::radius_lognormal_autorange) { std::uniform_real_distribution<> urd(zero, one); auto sigma = std::log(a_radius_gstd); auto mu = a_radius_mean; @@ -481,7 +482,7 @@ void SDInitProperties::getDistribution ( amrex::Vector& a_mass, for (int n = 0; n < a_np; n++) { auto tmp = lnrmin + urd(a_rng)*dlnr; auto dry_r = std::exp(tmp); - tmp_mass[n] = (amrex::Real(4.0)/three) * PI * dry_r * dry_r * dry_r * a_density; + tmp_mass[n] = four_thirds_pi * dry_r * dry_r * dry_r * a_density; auto term = std::exp(-std::log(dry_r/mu)*std::log(dry_r/mu)/(two*sigma*sigma)); tmp_mult[n] = (m_numdens * a_dV)/ (sigma*std::sqrt(2*PI)) * term; } @@ -494,7 +495,7 @@ void SDInitProperties::getDistribution ( amrex::Vector& a_mass, auto tmp = P_max + (one - P_max) * urd(a_rng); auto tmp2 = SD_erfinv(2 * tmp - 1); auto dry_r = mu * std::exp(sigma * std::sqrt(2) * tmp2); - tmp_mass[sd_id] = (amrex::Real(4.0)/three) * PI * dry_r * dry_r * dry_r * a_density; + tmp_mass[sd_id] = four_thirds_pi * dry_r * dry_r * dry_r * a_density; // set the multiplicity to the same as for the 99th percentile aerosol tmp_mult[sd_id] = (m_numdens * a_dV) * tail_mult; } @@ -508,3 +509,85 @@ void SDInitProperties::getDistribution ( amrex::Vector& a_mass, amrex::Abort("Unknown m_init_type!"); } } + +SDDistributionParams SDInitProperties::makeDistributionParams( + SDDistributionType a_init_type, + amrex::Real a_mass_min, + amrex::Real a_mass_max, + amrex::Real a_mass_mean, + amrex::Real a_radius_min, + amrex::Real a_radius_max, + amrex::Real a_radius_mean, + amrex::Real a_radius_gstd, + amrex::Real a_density, + amrex::Real a_cell_volume, + bool a_sampled_mult) const +{ + SDDistributionParams params; + params.dist_type = a_init_type; + params.mass_min = a_mass_min; + params.mass_max = a_mass_max; + params.mass_mean = a_mass_mean; + params.radius_min = a_radius_min; + params.radius_max = a_radius_max; + params.radius_mean = a_radius_mean; + params.radius_gstd = a_radius_gstd; + params.density = a_density; + params.numdens = m_numdens; + params.cell_volume = a_cell_volume; + params.delta = a_mass_mean - a_mass_min; + params.sampled_mult = a_sampled_mult ? 1 : 0; + + // Pre-compute values for log-normal distribution + params.sigma = std::log(a_radius_gstd); + + if (params.dist_type == SDDistributionType::mass_exponential) { + // For exponential with sampled multiplicity, we sample in log-space + // For constant multiplicity, we use inverse transform sampling + params.lnmin = std::log(a_mass_min); + params.lnrng = std::log(a_mass_max) - params.lnmin; + params.cdf_min = zero; + params.cdf_max = one; + } else if (params.dist_type == SDDistributionType::radius_log_normal || + params.dist_type == SDDistributionType::radius_lognormal_autorange) { + // For log-normal, compute CDF bounds for truncated sampling + // CDF of log-normal: Phi((ln(x) - ln(mu)) / sigma) + // where Phi is the standard normal CDF: Phi(z) = 0.5 * (1 + erf(z/sqrt(2))) + amrex::Real rmin = a_radius_min; + amrex::Real rmax = a_radius_max; + amrex::Real mu = a_radius_mean; + amrex::Real sigma = params.sigma; + + // For LogNormalAuto, find appropriate bounds iteratively + if (params.dist_type == SDDistributionType::radius_lognormal_autorange) { + rmin = amrex::Real(1e-9); + rmax = one; + amrex::Real tol = (m_numdens > 0 && a_cell_volume > 0) ? + one / (m_numdens * a_cell_volume) : amrex::Real(1e-6); + amrex::Real P_min = zero; + amrex::Real P_max = one; + while ((P_max >= one - tol) || (P_min <= tol)) { + if (P_max >= one - tol) rmax *= amrex::Real(0.99); + if (P_min <= tol) rmin *= amrex::Real(1.01); + P_min = myhalf * (one + std::erf(std::log(rmin / mu) / (sigma * std::sqrt(two)))); + P_max = myhalf * (one + std::erf(std::log(rmax / mu) / (sigma * std::sqrt(two)))); + } + // Update the params with computed bounds + params.radius_min = rmin; + params.radius_max = rmax; + } + + params.cdf_min = myhalf * (one + std::erf(std::log(rmin / mu) / (sigma * std::sqrt(two)))); + params.cdf_max = myhalf * (one + std::erf(std::log(rmax / mu) / (sigma * std::sqrt(two)))); + params.lnmin = std::log(rmin); + params.lnrng = std::log(rmax) - params.lnmin; + } else { + // Constant distribution - no special pre-computation needed + params.cdf_min = zero; + params.cdf_max = one; + params.lnmin = zero; + params.lnrng = zero; + } + + return params; +} diff --git a/Source/Particles/ERF_SuperDropletPC.H b/Source/Particles/ERF_SuperDropletPC.H index 71c6b5603..1354dfad0 100644 --- a/Source/Particles/ERF_SuperDropletPC.H +++ b/Source/Particles/ERF_SuperDropletPC.H @@ -7,6 +7,8 @@ #include "ERF_Constants.H" #include "ERF_MaterialProperties.H" #include "ERF_SuperDropletPCDefinitions.H" +#include "ERF_SuperDropletPCProcess.H" +#include class SuperDropletPC : public ERFPC { @@ -55,10 +57,9 @@ class SuperDropletPC : public ERFPC AMREX_FORCE_INLINE virtual void setSpeciesMaterial ( const Species::Name& a_name ) { - if (a_name == Species::Name::H2O) { - m_idx_w = m_species_mat.size(); - } + if (a_name == Species::Name::H2O) { m_idx_w = m_species_mat.size(); } m_species_mat.push_back(std::make_unique(a_name)); + m_device_props_initialized = false; } /*! \brief Set species material */ @@ -84,6 +85,7 @@ class SuperDropletPC : public ERFPC virtual void setAerosolMaterial ( const Species::Name& a_name ) { m_aerosol_mat.push_back(std::make_unique(a_name)); + m_device_props_initialized = false; } /*! \brief Set aerosol material */ @@ -93,6 +95,52 @@ class SuperDropletPC : public ERFPC for (auto& name : a_names) { setAerosolMaterial(name); } } + /*! \brief Update device properties if material properties change */ + void updateDeviceProperties(); + + /*! \brief Setup species and aerosol mass pointers from SOA */ + template + void setupMassPointers(SOAType& soa, SDPCDefn::SDSpeciesMassArr& sp_mass_ptrs, + SDPCDefn::SDAerosolMassArr& ae_mass_ptrs) const + { + for (int i = 0; i < m_num_species; i++) { + sp_mass_ptrs[i] = soa.GetRealData(idx_s(i, m_num_aerosols, m_num_species)).data(); + } + for (int i = 0; i < m_num_aerosols; i++) { + ae_mass_ptrs[i] = soa.GetRealData(idx_a(i, m_num_aerosols, m_num_species)).data(); + } + } + + /*! \brief Update particle radius and total mass from species/aerosol masses */ + AMREX_GPU_DEVICE AMREX_FORCE_INLINE + static void updateParticleAttributes( + int particle_idx, + amrex::ParticleReal* radius_ptr, + amrex::ParticleReal* mass_ptr, + int idx_w, + amrex::ParticleReal rho_w, + int num_sp, int num_ae, + const int* sp_sol_arr, + const int* ae_sol_arr, + const SDPCDefn::SDSpeciesMassArr sp_mass_ptrs, + const SDPCDefn::SDAerosolMassArr ae_mass_ptrs, + const amrex::ParticleReal* sp_rho_arr, + const amrex::ParticleReal* ae_rho_arr) + { + // Update effective radius + radius_ptr[particle_idx] = SD_effective_radius( + particle_idx, idx_w, rho_w, + num_sp, num_ae, + sp_sol_arr, ae_sol_arr, + sp_mass_ptrs, ae_mass_ptrs, + sp_rho_arr, ae_rho_arr); + + // Update total mass + mass_ptr[particle_idx] = SD_total_mass( + particle_idx, num_sp, num_ae, + sp_mass_ptrs, ae_mass_ptrs); + } + /*! \brief Get real-type particle attribute names */ [[nodiscard]] virtual amrex::Vector varNames () const override; @@ -146,18 +194,7 @@ class SuperDropletPC : public ERFPC amrex::Abort("SuperDropletPC::EvolveParticles() is intentionally disabled."); } - /*! \brief Advect particles for one time step - * \param[in] a_lev AMR level - * \param[in] a_time Current simulation time - * \param[in] a_dt Timestep for advection - * \param[in] a_flow_vel Array of face-based velocities - * \param[in] a_density Density field - * \param[in] a_pressure Pressure field - * \param[in] a_temperature Temperature field - * \param[in] a_z_phys_nd Array of terrain heights - * \param[in] a_bctypes Array of boundary condition types - * \param[in] a_recycle Flag to enable particle recycling - */ + /*! \brief Advect particles for one time step */ virtual void AdvectParticles ( int a_lev, amrex::Real a_time, amrex::Real a_dt, @@ -169,17 +206,7 @@ class SuperDropletPC : public ERFPC const BCTypeArr& a_bctypes, const bool a_recycle); - /*! \brief Grow/shrink particles for one time step - * \param[in] a_lev AMR level - * \param[in] a_dt Timestep for mass change - * \param[in] a_vap_name Name of the vapour species undergoing phase change - * \param[in] a_temperature Temperature field - * \param[in] a_pressure Pressure field - * \param[in] a_sat_pressure Saturation pressure field - * \param[in] a_sat_ratio Saturation ratio field - * \param[in] a_z_phys_nd Array of terrain heights - * \param[in] a_is_water Flag indicating if this is water species - */ + /*! \brief Grow/shrink particles for one time step */ virtual void MassChange ( int a_lev, amrex::Real a_dt, const Species::Name& a_vap_name, @@ -190,24 +217,13 @@ class SuperDropletPC : public ERFPC const amrex::Vector& a_z_phys_nd, const bool a_is_water); - /*! \brief Coalescence of super-droplets - * \param[in] a_lev AMR level - * \param[in] a_dt Timestep for coalescence - * \param[in] a_pressure Pressure field used in coalescence calculation - * \param[in] a_temperature Temperature field used in coalescence calculation - */ + /*! \brief Coalescence of super-droplets */ virtual void Coalescence ( int a_lev, amrex::Real a_dt, const amrex::MultiFab& a_pressure, const amrex::MultiFab& a_temperature); - /*! \brief Recycle super-droplets - * \param[in] a_lev AMR level - * \param[in] a_z_phys_nd Array of terrain heights - * \param[in] a_iter Current iteration number - * \param[in] a_dt Timestep size - * \param[in] a_recycle Flag to enable recycling - */ + /*! \brief Recycle super-droplets */ virtual void Recycle ( const int a_lev, const amrex::Vector& a_z_phys_nd, const int a_iter, @@ -249,18 +265,18 @@ class SuperDropletPC : public ERFPC /*! \brief Compute the mass flux component on a mesh */ virtual void massFlux ( amrex::MultiFab&, const int, const int a_comp=0 ) const; - /*! \brief Compute the condensate mass density on a mesh */ - virtual void speciesMassDensity ( amrex::MultiFab& a_mf, const int a_i ) const - { - speciesMassDensity(a_mf, a_i, zero, DBL_MAX); - } - /*! \brief Compute the condensate mass density on a mesh */ + /*! \brief Compute the species mass density on a mesh */ virtual void speciesMassDensity ( amrex::MultiFab&, - const int, - const amrex::Real, - const amrex::Real, + int, const int a_comp = 0) const; - /*! \brief Compute the condensate mass flux component on a mesh */ + + /*! \brief Compute the cloud/rain mass density on a mesh */ + virtual void cloudRainDensity ( amrex::MultiFab&, + amrex::Real, + amrex::Real, + const int a_comp = 0) const; + + /*! \brief Compute the species mass flux component on a mesh */ virtual void speciesMassFlux ( amrex::MultiFab&, const int, const int, const int a_comp=0 ) const; /*! \brief Compute the aerosol mass density on a mesh */ @@ -281,7 +297,7 @@ class SuperDropletPC : public ERFPC SDCoalescenceKernelType m_coalescence_kernel; /*!< Coalescence kernel */ bool m_include_brownian_coalescence; /*!< Include Brownian coalescence? */ - SDTerminalVelocityType m_term_vel_type; /*!< Terminal velocity model */ + SDTerminalVelocityType m_term_vel_type_w; /*!< Terminal velocity model for water droplets */ int m_num_sd_per_cell; /*!< Number of super-droplets per cell */ bool m_density_scaling; /*!< Scale initial number density with air density */ @@ -347,7 +363,7 @@ class SuperDropletPC : public ERFPC std::mt19937 m_rndeng; /*! species index of water */ - int m_idx_w; + int m_idx_w = -1; /*! inactive particles threshold */ amrex::Real m_deac_threshold; @@ -355,6 +371,20 @@ class SuperDropletPC : public ERFPC /*! save inactive particles to file */ bool m_save_inactive; + /*! Persistent device vectors for material properties */ + amrex::Gpu::DeviceVector m_sp_density; + amrex::Gpu::DeviceVector m_sp_solubility; + amrex::Gpu::DeviceVector m_sp_ionization; + amrex::Gpu::DeviceVector m_sp_mol_weight; + + amrex::Gpu::DeviceVector m_ae_density; + amrex::Gpu::DeviceVector m_ae_solubility; + amrex::Gpu::DeviceVector m_ae_ionization; + amrex::Gpu::DeviceVector m_ae_mol_weight; + + /*! Flag to track if device properties are initialized */ + bool m_device_props_initialized = false; + /* recycled particle position bounds */ amrex::Real m_recyc_xmin; amrex::Real m_recyc_xmax; @@ -363,6 +393,133 @@ class SuperDropletPC : public ERFPC amrex::Real m_recyc_zmin; amrex::Real m_recyc_zmax; + /*! Method to initialize device properties */ + void initializeDeviceProperties(); + + /*! \brief Build ProcessContext with geometry and species info */ + SDProcess::ProcessContext buildProcessContext(int a_lev) const + { + SDProcess::ProcessContext ctx; + const amrex::Geometry& geom = m_gdb->Geom(a_lev); + ctx.plo = geom.ProbLoArray(); + ctx.phi = geom.ProbHiArray(); + ctx.dxi = geom.InvCellSizeArray(); + ctx.dx = geom.CellSizeArray(); + ctx.domain = geom.Domain(); + for (int d = 0; d < AMREX_SPACEDIM; d++) { + ctx.is_periodic[d] = geom.isPeriodic(d) ? 1 : 0; + } + const auto cell_size = geom.CellSize(); + ctx.cell_volume = AMREX_D_TERM(cell_size[0], *cell_size[1], *cell_size[2]); + ctx.num_species = m_num_species; + ctx.num_aerosols = m_num_aerosols; + ctx.idx_water = m_idx_w; + ctx.rho_water = m_species_mat[m_idx_w]->m_density; + return ctx; + } + + /*! \brief Setup particle attribute pointers for a tile */ + template + void setupParticlePointers( + SOAType& soa, + AOSType& aos, + SDProcess::ParticlePointers& ptrs) const + { + using namespace SDPCDefn; + + constexpr int rtoff_i = SuperDropletsIntIdxSoA::ncomps; + constexpr int rtoff_r = SuperDropletsRealIdxSoA::ncomps; + + ptrs.num_particles = aos.numParticles(); + ptrs.mass_ptr = soa.GetRealData(SuperDropletsRealIdxSoA::mass).data(); + ptrs.radius_ptr = soa.GetRealData(rtoff_r + SuperDropletsRealIdxSoA_RT::radius).data(); + ptrs.active_ptr = soa.GetIntData(rtoff_i + SuperDropletsIntIdxSoA_RT::active).data(); + ptrs.v_ptr[0] = soa.GetRealData(SuperDropletsRealIdxSoA::vx).data(); + ptrs.v_ptr[1] = soa.GetRealData(SuperDropletsRealIdxSoA::vy).data(); + ptrs.v_ptr[2] = soa.GetRealData(SuperDropletsRealIdxSoA::vz).data(); + ptrs.vterm_ptr = soa.GetRealData(rtoff_r + SuperDropletsRealIdxSoA_RT::term_vel).data(); + ptrs.mult_ptr = soa.GetRealData(rtoff_r + SuperDropletsRealIdxSoA_RT::multiplicity).data(); + setupMassPointers(soa, ptrs.sp_mass_ptrs, ptrs.ae_mass_ptrs); + if (!m_device_props_initialized) { + const_cast(this)->initializeDeviceProperties(); + } + ptrs.sp_rho_arr = m_sp_density.data(); + ptrs.sp_sol_arr = m_sp_solubility.data(); + ptrs.sp_ion_arr = m_sp_ionization.data(); + ptrs.sp_mw_arr = m_sp_mol_weight.data(); + ptrs.ae_rho_arr = m_ae_density.data(); + ptrs.ae_sol_arr = m_ae_solubility.data(); + ptrs.ae_ion_arr = m_ae_ionization.data(); + ptrs.ae_mw_arr = m_ae_mol_weight.data(); + } + + /*! \brief Common body for forEachParticleTile iterations */ + template + void forEachParticleTileBody( + ParIterType& pti, int a_lev, + const SDProcess::ProcessContext& ctx, + TileFunc&& func) + { + int grid = pti.index(); + auto& ptile = ParticlesAt(a_lev, pti); + auto& aos = ptile.GetArrayOfStructs(); + auto& soa = ptile.GetStructOfArrays(); + if (aos.numParticles() == 0) { return; } + auto* p_pbox = aos().data(); + SDProcess::ParticlePointers ptrs; + setupParticlePointers(soa, aos, ptrs); + func(pti, grid, p_pbox, ptrs, ctx); + } + + /*! \brief Iterate over particle tiles with common setup (parallel) */ + template + AMREX_FORCE_INLINE + void forEachParticleTile( + int a_lev, + const SDProcess::ProcessContext& ctx, + TileFunc&& func) + { +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (ParIterType pti(*this, a_lev); pti.isValid(); ++pti) { + forEachParticleTileBody(pti, a_lev, ctx, std::forward(func)); + } + } + + /*! \brief Iterate over particle tiles WITHOUT OpenMP (serial) + * Use when per-tile operations are not thread-safe (e.g., DenseBins). + */ + template + AMREX_FORCE_INLINE + void forEachParticleTileSerial(int a_lev, const SDProcess::ProcessContext& ctx, TileFunc&& func) + { + for (ParIterType pti(*this, a_lev); pti.isValid(); ++pti) { + forEachParticleTileBody(pti, a_lev, ctx, std::forward(func)); + } + } + + /*! \brief Lightweight iteration (no ProcessContext) */ + template + void forEachParticleTile(int a_lev, TileFunc&& func) + { +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (ParIterType pti(*this, a_lev); pti.isValid(); ++pti) { + int grid = pti.index(); + auto& ptile = ParticlesAt(a_lev, pti); + auto& aos = ptile.GetArrayOfStructs(); + auto& soa = ptile.GetStructOfArrays(); + const int num_particles = aos.numParticles(); + if (num_particles == 0) { continue; } + auto* p_pbox = aos().data(); + SDProcess::ParticlePointers ptrs; + setupParticlePointers(soa, aos, ptrs); + func(pti, grid, p_pbox, ptrs, num_particles); + } + } + /*! \brief read inputs from file */ virtual void readInputs () override { diff --git a/Source/Particles/ERF_SuperDropletPCAddParticles.cpp b/Source/Particles/ERF_SuperDropletPCAddParticles.cpp index ff23c260b..fbec48fea 100644 --- a/Source/Particles/ERF_SuperDropletPCAddParticles.cpp +++ b/Source/Particles/ERF_SuperDropletPCAddParticles.cpp @@ -305,112 +305,62 @@ void SuperDropletPC::addParticles ( const MFPtr& a_height_ptr, /*!< terrain */ ae_mass_ptrs[i] = soa.GetRealData(idx_a(i,num_ae,num_sp)).data() + size_old; } - Gpu::DeviceVector species_mass_d(num_sp*np); - { - Vector multiplicity_h(np, zero); - for (int i = 0; i < num_sp; i++) { - Vector species_mass_h; - if (sampled_multiplicity) { - a_init.getSpeciesDistribution( species_mass_h, - multiplicity_h, - cell_volume, - i, - np, - m_species_mat[i]->m_density, - m_rndeng ); - } else { - a_init.getSpeciesDistribution( species_mass_h, - i, - np, - m_species_mat[i]->m_density, - m_rndeng ); - } - Gpu::copy( Gpu::hostToDevice, - species_mass_h.begin(), - species_mass_h.end(), - species_mass_d.begin() + (i*np) ); - } + // GPU-direct distribution generation: no large host memory allocation needed + // Create distribution parameters on host (small, O(num_species + num_aerosols)) + Vector sp_params_h(num_sp); + for (int i = 0; i < num_sp; i++) { + sp_params_h[i] = a_init.getSpeciesDistParams(i, m_species_mat[i]->m_density, cell_volume, sampled_multiplicity); + } + Vector ae_params_h(num_ae); + for (int i = 0; i < num_ae; i++) { + ae_params_h[i] = a_init.getAerosolDistParams(i, m_aerosol_mat[i]->m_density, cell_volume, sampled_multiplicity); } - Gpu::DeviceVector aerosol_mass_d(num_ae*np); - Gpu::DeviceVector multiplicity_d(np); + + // Copy parameters to device + Gpu::DeviceVector sp_params_d(num_sp); + Gpu::DeviceVector ae_params_d(num_ae); + Gpu::copy(Gpu::hostToDevice, sp_params_h.begin(), sp_params_h.end(), sp_params_d.begin()); + Gpu::copy(Gpu::hostToDevice, ae_params_h.begin(), ae_params_h.end(), ae_params_d.begin()); + + // Only need temporary storage for raw multiplicity values (for sampled case) + Gpu::DeviceVector multiplicity_d(np, zero); + + auto* sp_params_ptr = sp_params_d.data(); + auto* ae_params_ptr = ae_params_d.data(); + auto* mult_ptr_init = multiplicity_d.data(); + + // Generate species and aerosol distributions directly into particle SoA + ParallelForRNG(np, [=] AMREX_GPU_DEVICE (int n, const RandomEngine& engine) noexcept { - Vector multiplicity_h(np, zero); - for (int i = 0; i < num_ae; i++) { - Vector aerosol_mass_h; + // Sample species masses directly into particle storage + for (int s = 0; s < num_sp; s++) { + Real sp_mult = zero; + sp_mass_ptrs[s][n] = SD_sample_mass_gpu(sp_params_ptr[s], engine, sp_mult); + } + + // Sample aerosol masses directly into particle storage + Real total_mult = zero; + for (int a = 0; a < num_ae; a++) { + Real ae_mult = zero; + ae_mass_ptrs[a][n] = SD_sample_mass_gpu(ae_params_ptr[a], engine, ae_mult); if (sampled_multiplicity) { - a_init.getAerosolDistribution( aerosol_mass_h, - multiplicity_h, - cell_volume, - i, - np, - m_aerosol_mat[i]->m_density, - m_rndeng ); - } else { - a_init.getAerosolDistribution( aerosol_mass_h, - i, - np, - m_aerosol_mat[i]->m_density, - m_rndeng ); + total_mult += ae_mult; } - Gpu::copy( Gpu::hostToDevice, - aerosol_mass_h.begin(), - aerosol_mass_h.end(), - aerosol_mass_d.begin() + (i*np) ); } if (sampled_multiplicity) { - Gpu::copy( Gpu::hostToDevice, - multiplicity_h.begin(), - multiplicity_h.end(), - multiplicity_d.begin() ); + mult_ptr_init[n] = total_mult; } - } - Gpu::synchronize(); + }); - Gpu::DeviceVector sp_density(num_sp); - Gpu::DeviceVector sp_solubility(num_sp); - { - Vector sp_density_h(num_sp); - Vector sp_solubility_h(num_sp); - for (int i = 0; i < num_sp; i++) { - sp_density_h[i] = m_species_mat[i]->m_density; - sp_solubility_h[i] = static_cast(m_species_mat[i]->m_is_soluble); - } - Gpu::copy( Gpu::hostToDevice, - sp_density_h.begin(), - sp_density_h.end(), - sp_density.begin() ); - Gpu::copy( Gpu::hostToDevice, - sp_solubility_h.begin(), - sp_solubility_h.end(), - sp_solubility.begin() ); - } + Gpu::synchronize(); - Gpu::DeviceVector ae_density(num_ae); - Gpu::DeviceVector ae_solubility(num_ae); - { - Vector ae_density_h(num_ae); - Vector ae_solubility_h(num_ae); - for (int i = 0; i < num_ae; i++) { - ae_density_h[i] = m_aerosol_mat[i]->m_density; - ae_solubility_h[i] = static_cast(m_aerosol_mat[i]->m_is_soluble); - } - Gpu::copy( Gpu::hostToDevice, - ae_density_h.begin(), - ae_density_h.end(), - ae_density.begin() ); - Gpu::copy( Gpu::hostToDevice, - ae_solubility_h.begin(), - ae_solubility_h.end(), - ae_solubility.begin() ); - } + if (!m_device_props_initialized) { initializeDeviceProperties(); } + const ParticleReal* sp_rho_arr = m_sp_density.data(); + const int* sp_sol_arr = m_sp_solubility.data(); + const ParticleReal* ae_rho_arr = m_ae_density.data(); + const int* ae_sol_arr = m_ae_solubility.data(); - auto species_mass = species_mass_d.data(); - auto aerosol_mass = aerosol_mass_d.data(); auto mult_arr = multiplicity_d.data(); - auto sp_rho_arr = sp_density.data(); - auto sp_sol_arr = sp_solubility.data(); - auto ae_rho_arr = ae_density.data(); - auto ae_sol_arr = ae_solubility.data(); auto num_superdroplets_arr = num_superdroplets[mfi].array(); auto random_place = m_place_randomly_in_cells; @@ -501,12 +451,7 @@ void SuperDropletPC::addParticles ( const MFPtr& a_height_ptr, /*!< terrain */ num_to_add -= mult_ptr[n]; if (mult_ptr[n] == 0) { mult_ptr[n] = 1; } - for (int ctr = 0; ctr < num_sp; ctr++) { - sp_mass_ptrs[ctr][n] = species_mass[ctr*np+n]; - } - for (int ctr = 0; ctr < num_ae; ctr++) { - ae_mass_ptrs[ctr][n] = aerosol_mass[ctr*np+n]; - } + // Species and aerosol masses already sampled directly into particle SoA radius_ptr[n] = SD_effective_radius( n, idx_w, rho_w, @@ -523,7 +468,6 @@ void SuperDropletPC::addParticles ( const MFPtr& a_height_ptr, /*!< terrain */ uid_ptr[n] = ParticleReal((pid+n-1)*nprocs + my_proc + 1); } }); - Gpu::synchronize(); const auto height_arr = (*a_height_ptr)[mfi].array(); ParallelFor(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept diff --git a/Source/Particles/ERF_SuperDropletPCAdvection.cpp b/Source/Particles/ERF_SuperDropletPCAdvection.cpp index 8f2fecda1..28930c5b1 100644 --- a/Source/Particles/ERF_SuperDropletPCAdvection.cpp +++ b/Source/Particles/ERF_SuperDropletPCAdvection.cpp @@ -5,10 +5,14 @@ #include #include "ERF_IndexDefines.H" #include "ERF_TerminalVelocity.H" +#include "ERF_InterpolationUtils.H" using namespace amrex; using namespace SDPCDefn; +/*! \brief Field indices for advection interpolation */ +AMREX_ENUM(InterpFieldsAdv, density, pressure, temperature, NUM_FIELDS); + /*! \brief Advect superdroplet particles for one time step * \param[in] a_lev AMR level * \param[in] a_time Current simulation time @@ -45,38 +49,25 @@ void SuperDropletPC::AdvectParticles ( int a_lev, AMREX_ASSERT(!a_flow_vel[1].contains_nan());, AMREX_ASSERT(!a_flow_vel[2].contains_nan());); + const auto ctx = buildProcessContext(a_lev); const Geometry& geom = m_gdb->Geom(a_lev); - const auto plo = geom.ProbLoArray(); - const auto dxi = geom.InvCellSizeArray(); - - const auto is_periodic = geom.isPeriodic(); - auto is_periodic_z = is_periodic[2]; + const auto is_periodic_z = geom.isPeriodic(2); const bool advect_w_flow = m_advect_w_flow; const bool advect_w_gravity = m_advect_w_gravity; const bool prescribed_advection = m_prescribed_advection; - const int num_sp = m_num_species; - const int num_ae = m_num_aerosols; + const auto vterm_type_w = m_term_vel_type_w; - const Real rho_w = m_species_mat[m_idx_w]->m_density; - int idx_w = m_idx_w; + // Terminal velocity calculator (shared across tiles) + TerminalVelocity term_vel { ctx.rho_water }; -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (ParIterType pti(*this, a_lev); pti.isValid(); ++pti) { - int grid = pti.index(); - auto& ptile = ParticlesAt(a_lev, pti); - auto& aos = ptile.GetArrayOfStructs(); - auto& soa = ptile.GetStructOfArrays(); - const int n = aos.numParticles(); - auto *p_pbox = aos().data(); - - Array v_ptr; - v_ptr[0] = soa.GetRealData(SuperDropletsRealIdxSoA::vx).data(); - v_ptr[1] = soa.GetRealData(SuperDropletsRealIdxSoA::vy).data(); - v_ptr[2] = soa.GetRealData(SuperDropletsRealIdxSoA::vz).data(); + forEachParticleTile(a_lev, ctx, + [&](ParIterType& /*pti*/, int grid, ParticleType* p_pbox, + const SDProcess::ParticlePointers& ptrs, + const SDProcess::ProcessContext& ctx) + { + auto zheight = (*z_height)[grid].array(); const FArrayBox* fab[AMREX_SPACEDIM] = { AMREX_D_DECL(&(a_flow_vel[0][grid]), &(a_flow_vel[1][grid]), @@ -90,85 +81,36 @@ void SuperDropletPC::AdvectParticles ( int a_lev, const auto& temperature_arr = a_temperature[grid].array(); const auto& pressure_arr = a_pressure[grid].array(); - auto zheight = (*z_height)[grid].array(); - - int rtoff_i = SuperDropletsIntIdxSoA::ncomps; - auto* active_ptr = soa.GetIntData(rtoff_i+SuperDropletsIntIdxSoA_RT::active).data(); - int rtoff_r = SuperDropletsRealIdxSoA::ncomps; - auto* vterm_ptr = soa.GetRealData(rtoff_r+SuperDropletsRealIdxSoA_RT::term_vel).data(); - - SDSpeciesMassArr sp_mass_ptrs; - Gpu::DeviceVector sp_density(num_sp); - Gpu::DeviceVector sp_solubility(num_sp); - { - Vector sp_density_h(num_sp); - Vector sp_solubility_h(num_sp); - for (int i = 0; i < num_sp; i++) { - sp_mass_ptrs[i] = soa.GetRealData(idx_s(i,num_ae,num_sp)).data(); - sp_density_h[i] = m_species_mat[i]->m_density; - sp_solubility_h[i] = static_cast(m_species_mat[i]->m_is_soluble); - } - Gpu::copy( Gpu::hostToDevice, - sp_density_h.begin(), - sp_density_h.end(), - sp_density.begin() ); - Gpu::copy( Gpu::hostToDevice, - sp_solubility_h.begin(), - sp_solubility_h.end(), - sp_solubility.begin() ); - } - - SDAerosolMassArr ae_mass_ptrs; - Gpu::DeviceVector ae_density(num_ae); - Gpu::DeviceVector ae_solubility(num_ae); - { - Vector ae_density_h(num_ae); - Vector ae_solubility_h(num_ae); - for (int i = 0; i < num_ae; i++) { - ae_mass_ptrs[i] = soa.GetRealData(idx_a(i,num_ae,num_sp)).data(); - ae_density_h[i] = m_aerosol_mat[i]->m_density; - ae_solubility_h[i] = static_cast(m_aerosol_mat[i]->m_is_soluble); - } - Gpu::copy( Gpu::hostToDevice, - ae_density_h.begin(), - ae_density_h.end(), - ae_density.begin() ); - Gpu::copy( Gpu::hostToDevice, - ae_solubility_h.begin(), - ae_solubility_h.end(), - ae_solubility.begin() ); - } - - TerminalVelocity term_vel { m_species_mat[m_idx_w]->m_density }; - auto term_vel_type = m_term_vel_type; - - auto sp_rho_arr = sp_density.data(); - auto sp_sol_arr = sp_solubility.data(); - auto ae_rho_arr = ae_density.data(); - auto ae_sol_arr = ae_solubility.data(); - - ParallelFor(n, [=] AMREX_GPU_DEVICE (int i) + ParallelFor(ptrs.num_particles, [=] AMREX_GPU_DEVICE (int i) { ParticleType& p = p_pbox[i]; if (p.id() <= 0) { return; } - if (active_ptr[i] == 0) { return; } + if (ptrs.active_ptr[i] == 0) { return; } ParticleReal v[AMREX_SPACEDIM]; v[0] = v[1] = v[2] = zero; if (is_periodic_z) { - mac_interpolate(p, plo, dxi, umacarr, v); + mac_interpolate(p, ctx.plo, ctx.dxi, umacarr, v); } else { - mac_interpolate_mapped_z(p, plo, dxi, umacarr, zheight, v); + mac_interpolate_mapped_z(p, ctx.plo, ctx.dxi, umacarr, zheight, v); } + // Interpolate density, pressure, temperature at particle position + constexpr int nf = static_cast(InterpFieldsAdv::NUM_FIELDS); + ParticleReal fv[nf]; + const Array4 fa[nf] = { + density_arr, pressure_arr, temperature_arr + }; + ERF::Interpolation::interpolateFields( + p, ctx.plo, ctx.dxi, fa, fv, nf, + is_periodic_z ? 1 : 0, is_periodic_z ? nullptr : &zheight + ); + const auto density = fv[static_cast(InterpFieldsAdv::density)]; + const auto pressure = fv[static_cast(InterpFieldsAdv::pressure)]; + const auto temperature = fv[static_cast(InterpFieldsAdv::temperature)]; + if (prescribed_advection) { - ParticleReal density; - if (is_periodic_z) { - cic_interpolate( p, plo, dxi, density_arr, &density, 1 ); - } else { - cic_interpolate_mapped_z( p, plo, dxi, density_arr, zheight, &density, 1 ); - } if (a_time < 600) { v[2] = two*sin(PI*a_time/600)/density; } else { @@ -177,40 +119,31 @@ void SuperDropletPC::AdvectParticles ( int a_lev, } // compute effective radius - auto r_eff = SD_effective_radius( i, idx_w, - rho_w, - num_sp, num_ae, - sp_sol_arr, ae_sol_arr, - sp_mass_ptrs, ae_mass_ptrs, - sp_rho_arr, ae_rho_arr ); - + auto r_eff = SD_effective_radius( i, ctx.idx_water, + ctx.rho_water, + ctx.num_species, ctx.num_aerosols, + ptrs.sp_sol_arr, ptrs.ae_sol_arr, + ptrs.sp_mass_ptrs, ptrs.ae_mass_ptrs, + ptrs.sp_rho_arr, ptrs.ae_rho_arr ); ParticleReal terminal_vel = zero; - if (term_vel_type == SDTerminalVelocityType::AtlasUlbrich) { + if (vterm_type_w == SDTerminalVelocityType::AtlasUlbrich) { terminal_vel = term_vel.AtlasUlbrich( r_eff ); - } else if (term_vel_type == SDTerminalVelocityType::RogersYau) { + } else if (vterm_type_w == SDTerminalVelocityType::RogersYau) { terminal_vel = term_vel.RogersYau( r_eff ); - } else if (term_vel_type == SDTerminalVelocityType::CloudRainShima) { - ParticleReal density, pressure, temperature; - if (is_periodic_z) { - cic_interpolate( p, plo, dxi, density_arr, &density, 1 ); - cic_interpolate( p, plo, dxi, pressure_arr, &pressure, 1 ); - cic_interpolate( p, plo, dxi, temperature_arr, &temperature, 1 ); - } else { - cic_interpolate_mapped_z( p, plo, dxi, density_arr, zheight, &density, 1 ); - cic_interpolate_mapped_z( p, plo, dxi, pressure_arr, zheight, &pressure, 1 ); - cic_interpolate_mapped_z( p, plo, dxi, temperature_arr, zheight, &temperature, 1 ); - } + } else if (vterm_type_w == SDTerminalVelocityType::CloudRainShima) { terminal_vel = term_vel.CloudRainShima( r_eff, density, pressure, temperature ); + } else { + amrex::Abort("Invalid option for terminal velocity model"); } for (int dim=0; dim < AMREX_SPACEDIM; dim++) { - v_ptr[dim][i] = v[dim]; + ptrs.v_ptr[dim][i] = v[dim]; } - vterm_ptr[i] = terminal_vel; + ptrs.vterm_ptr[i] = terminal_vel; if (advect_w_flow) { for (int dim=0; dim < AMREX_SPACEDIM; dim++) { @@ -222,14 +155,35 @@ void SuperDropletPC::AdvectParticles ( int a_lev, } // Update z-coordinate carried by the particle - update_location_idata(p,plo,dxi,zheight); + update_location_idata(p,ctx.plo,ctx.dxi,zheight); }); Gpu::synchronize(); - } + }); // end forEachParticleTile applyBoundaryTreatment(a_lev, a_z_phys_nd, a_bctypes, a_recycle); Redistribute(); + + // After redistribution, update k-index for particles that moved to new tiles. + // This is needed because update_location_idata skips particles outside the + // local tile bounds (to avoid out-of-bounds array access). After Redistribute, + // those particles are now on the correct tile and can update their k-index. + forEachParticleTile(a_lev, ctx, + [&](ParIterType& /*pti*/, int grid, ParticleType* p_pbox, + const SDProcess::ParticlePointers& ptrs, + const SDProcess::ProcessContext& ctx) + { + auto zheight = (*z_height)[grid].array(); + + ParallelFor(ptrs.num_particles, [=] AMREX_GPU_DEVICE (int i) + { + ParticleType& p = p_pbox[i]; + if (p.id() > 0) { + update_location_idata(p, ctx.plo, ctx.dxi, zheight); + } + }); + Gpu::synchronize(); + }); } #endif diff --git a/Source/Particles/ERF_SuperDropletPCBoundaries.cpp b/Source/Particles/ERF_SuperDropletPCBoundaries.cpp index a2fb1efb6..8d1619ff2 100644 --- a/Source/Particles/ERF_SuperDropletPCBoundaries.cpp +++ b/Source/Particles/ERF_SuperDropletPCBoundaries.cpp @@ -3,10 +3,12 @@ #ifdef ERF_USE_PARTICLES #include +#include "ERF_InterpolationUtils.H" using namespace amrex; using namespace SDPCDefn; + /*! Handle the boundaries for the particles */ void SuperDropletPC::applyBoundaryTreatment ( int a_lev, const Vector& a_z_phys_nd, @@ -16,24 +18,7 @@ void SuperDropletPC::applyBoundaryTreatment ( int a_lev, BL_PROFILE("SuperDropletPC::applyBoundaryTreatment()"); const MFPtr& z_height = a_z_phys_nd[a_lev]; - const Geometry& geom = m_gdb->Geom(a_lev); - const auto plo = geom.ProbLoArray(); - const auto phi = geom.ProbHiArray(); - const auto dx_h = Geom(m_lev).CellSize(); - const Real cell_volume = dx_h[0]*dx_h[1]*dx_h[2]; - const auto dx = geom.CellSizeArray(); - const auto dxi = geom.InvCellSizeArray(); - const auto domain = geom.Domain(); - const auto is_periodic = geom.isPeriodicArray(); - - const int k_lo = domain.smallEnd(2); - const int k_hi = domain.bigEnd(2); - - const int n_species = m_num_species; - const int n_aerosols = m_num_aerosols; - const Real rho_w = m_species_mat[m_idx_w]->m_density; - const int idx_w = m_idx_w; - + const auto ctx = buildProcessContext(a_lev); const auto save_inac = m_save_inactive; // number of super-droplets per cell @@ -41,88 +26,59 @@ void SuperDropletPC::applyBoundaryTreatment ( int a_lev, // number of physical particles per cell Real num_par_per_cell = zero; for (int i = 0; i < m_num_initializations; i++) { - num_par_per_cell += m_initializations[i]->numParticlesPerCell(cell_volume); + num_par_per_cell += m_initializations[i]->numParticlesPerCell(ctx.cell_volume); } auto multiplicity = (num_sd_per_cell > 0 ? num_par_per_cell / num_sd_per_cell : zero); Long num_deactivated_particles = 0; -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (ParIterType pti(*this, a_lev); pti.isValid(); ++pti) { - int grid = pti.index(); - auto& ptile = ParticlesAt(a_lev, pti); - auto& aos = ptile.GetArrayOfStructs(); - auto& soa = ptile.GetStructOfArrays(); - const int n = aos.numParticles(); - auto *p_pbox = aos().data(); - - Array v_ptr; - v_ptr[0] = soa.GetRealData(SuperDropletsRealIdxSoA::vx).data(); - v_ptr[1] = soa.GetRealData(SuperDropletsRealIdxSoA::vy).data(); - v_ptr[2] = soa.GetRealData(SuperDropletsRealIdxSoA::vz).data(); - auto* mass_ptr = soa.GetRealData(SuperDropletsRealIdxSoA::mass).data(); - + forEachParticleTile(a_lev, ctx, + [&](ParIterType& /*pti*/, int grid, ParticleType* p_pbox, + const SDProcess::ParticlePointers& ptrs, + const SDProcess::ProcessContext& ctx) + { auto zheight = (*z_height)[grid].array(); - int rtoff_i = SuperDropletsIntIdxSoA::ncomps; - auto* active_ptr = soa.GetIntData(rtoff_i+SuperDropletsIntIdxSoA_RT::active).data(); - int rtoff_r = SuperDropletsRealIdxSoA::ncomps; - auto* radius_ptr = soa.GetRealData(rtoff_r+SuperDropletsRealIdxSoA_RT::radius).data(); - auto* vterm_ptr = soa.GetRealData(rtoff_r+SuperDropletsRealIdxSoA_RT::term_vel).data(); - auto* mult_ptr = soa.GetRealData(rtoff_r+SuperDropletsRealIdxSoA_RT::multiplicity).data(); - - SDSpeciesMassArr species_mass_ptrs; - for (int i = 0; i < n_species; i++) { - species_mass_ptrs[i] = soa.GetRealData(idx_s(i,n_aerosols,n_species)).data(); - } - - SDAerosolMassArr aerosol_mass_ptrs; - for (int i = 0; i < n_aerosols; i++) { - aerosol_mass_ptrs[i] = soa.GetRealData(idx_a(i,n_aerosols,n_species)).data(); - } - Gpu::Buffer deactivated_particles({0}); auto* deactivated_particles_ptr = deactivated_particles.data(); - ParallelFor(n, [=] AMREX_GPU_DEVICE (int i) + ParallelFor(ptrs.num_particles, [=] AMREX_GPU_DEVICE (int i) { ParticleType& p = p_pbox[i]; if (p.id() <= 0) { return; } - if (active_ptr[i] == 0) { return; } + if (ptrs.active_ptr[i] == 0) { return; } // check for ground impact { - auto z_ground = plo[2]; + auto z_ground = ctx.plo[2]; { - auto iv = getParticleCell(p, plo, dxi, domain); - z_ground = zheight(iv[0],iv[1],k_lo); + auto iv = getParticleCell(p, ctx.plo, ctx.dxi, ctx.domain); + z_ground = zheight(iv[0],iv[1],ctx.domain.smallEnd(2)); } if (p.pos(2) < z_ground) { - p.pos(2) = z_ground + Real(0.01)*dx[2]; - v_ptr[0][i] = v_ptr[1][i] = v_ptr[2][i] = vterm_ptr[i] = zero; - active_ptr[i] = 0; + p.pos(2) = z_ground + Real(0.01)*ctx.dx[2]; + ptrs.v_ptr[0][i] = ptrs.v_ptr[1][i] = ptrs.v_ptr[2][i] = ptrs.vterm_ptr[i] = zero; + ptrs.active_ptr[i] = 0; if ((!a_recycle) && (!save_inac)) { p.id() = -1; } Gpu::Atomic::Add(deactivated_particles_ptr, Long(1)); - update_location_idata(p,plo,dxi,zheight); + update_location_idata(p,ctx.plo,ctx.dxi,zheight); } } // check for top boundary exits { - auto z_roof = phi[2]; + auto z_roof = ctx.phi[2]; { - auto iv = getParticleCell(p, plo, dxi, domain); - z_roof = zheight(iv[0],iv[1],k_hi+1); + auto iv = getParticleCell(p, ctx.plo, ctx.dxi, ctx.domain); + z_roof = zheight(iv[0],iv[1],ctx.domain.bigEnd(2)+1); } if (p.pos(2) > z_roof) { - p.pos(2) = z_roof - dx[2]; - v_ptr[0][i] = v_ptr[1][i] = v_ptr[2][i] = vterm_ptr[i] = zero; - active_ptr[i] = 0; + p.pos(2) = z_roof - ctx.dx[2]; + ptrs.v_ptr[0][i] = ptrs.v_ptr[1][i] = ptrs.v_ptr[2][i] = ptrs.vterm_ptr[i] = zero; + ptrs.active_ptr[i] = 0; if ((!a_recycle) && (!save_inac)) { p.id() = -1; } Gpu::Atomic::Add(deactivated_particles_ptr, Long(1)); - update_location_idata(p,plo,dxi,zheight); + update_location_idata(p,ctx.plo,ctx.dxi,zheight); } } @@ -130,8 +86,8 @@ void SuperDropletPC::applyBoundaryTreatment ( int a_lev, for (int d = 0; d < 2; d++) { // domain bounds - auto x_min = plo[d]; - auto x_max = phi[d]; + auto x_min = ctx.plo[d]; + auto x_max = ctx.phi[d]; auto bc_lo = a_bctypes[Orientation(d,Orientation::low)]; auto bc_hi = a_bctypes[Orientation(d,Orientation::high)]; @@ -140,9 +96,9 @@ void SuperDropletPC::applyBoundaryTreatment ( int a_lev, if ((bc_lo == ERF_BC::slip_wall) || (bc_lo == ERF_BC::no_slip_wall)) { - p.pos(d) = x_min + Real(0.01)*dx[d]; - v_ptr[0][i] = v_ptr[1][i] = v_ptr[2][i] = vterm_ptr[i] = zero; - active_ptr[i] = 0; + p.pos(d) = x_min + Real(0.01)*ctx.dx[d]; + ptrs.v_ptr[0][i] = ptrs.v_ptr[1][i] = ptrs.v_ptr[2][i] = ptrs.vterm_ptr[i] = zero; + ptrs.active_ptr[i] = 0; if ((!a_recycle) && (!save_inac)) { p.id() = -1; } Gpu::Atomic::Add(deactivated_particles_ptr, Long(1)); @@ -150,25 +106,19 @@ void SuperDropletPC::applyBoundaryTreatment ( int a_lev, auto delta = x_min - p.pos(d); p.pos(d) = x_max - delta; - if (!is_periodic[d]) { - v_ptr[0][i] = v_ptr[1][i] = v_ptr[2][i] = vterm_ptr[i] = zero; + if (!ctx.is_periodic[d]) { + ptrs.v_ptr[0][i] = ptrs.v_ptr[1][i] = ptrs.v_ptr[2][i] = ptrs.vterm_ptr[i] = zero; - for (int ctr = 0; ctr < n_species; ctr++) { - species_mass_ptrs[ctr][i] = zero; + for (int ctr = 0; ctr < ctx.num_species; ctr++) { + ptrs.sp_mass_ptrs[ctr][i] = zero; } - ParticleReal aerosol_mass_total = zero; - for (int ctr = 0; ctr < n_aerosols; ctr++) { - aerosol_mass_total += aerosol_mass_ptrs[ctr][i]; - } + ptrs.mult_ptr[i] = multiplicity; - auto par_radius = Real(1.0e-15); - auto cond_mass = (Real(4.0)/three)*PI - * par_radius*par_radius*par_radius*rho_w; - radius_ptr[i] = par_radius; - species_mass_ptrs[idx_w][i] = cond_mass; - mass_ptr[i] = cond_mass + aerosol_mass_total; - mult_ptr[i] = multiplicity; + SuperDropletPC::updateParticleAttributes( + i, ptrs.radius_ptr, ptrs.mass_ptr, ctx.idx_water, ctx.rho_water, + ctx.num_species, ctx.num_aerosols, ptrs.sp_sol_arr, ptrs.ae_sol_arr, + ptrs.sp_mass_ptrs, ptrs.ae_mass_ptrs, ptrs.sp_rho_arr, ptrs.ae_rho_arr); } } @@ -177,9 +127,9 @@ void SuperDropletPC::applyBoundaryTreatment ( int a_lev, if ((bc_hi == ERF_BC::slip_wall) || (bc_hi == ERF_BC::no_slip_wall)) { - p.pos(d) = x_max - Real(0.01)*dx[d]; - v_ptr[0][i] = v_ptr[1][i] = v_ptr[2][i] = vterm_ptr[i] = zero; - active_ptr[i] = 0; + p.pos(d) = x_max - Real(0.01)*ctx.dx[d]; + ptrs.v_ptr[0][i] = ptrs.v_ptr[1][i] = ptrs.v_ptr[2][i] = ptrs.vterm_ptr[i] = zero; + ptrs.active_ptr[i] = 0; if ((!a_recycle) && (!save_inac)) { p.id() = -1; } Gpu::Atomic::Add(deactivated_particles_ptr, Long(1)); @@ -187,25 +137,19 @@ void SuperDropletPC::applyBoundaryTreatment ( int a_lev, auto delta = p.pos(d) - x_max; p.pos(d) = x_min + delta; - if (!is_periodic[d]) { - v_ptr[0][i] = v_ptr[1][i] = v_ptr[2][i] = vterm_ptr[i] = zero; + if (!ctx.is_periodic[d]) { + ptrs.v_ptr[0][i] = ptrs.v_ptr[1][i] = ptrs.v_ptr[2][i] = ptrs.vterm_ptr[i] = zero; - for (int ctr = 0; ctr < n_species; ctr++) { - species_mass_ptrs[ctr][i] = zero; + for (int ctr = 0; ctr < ctx.num_species; ctr++) { + ptrs.sp_mass_ptrs[ctr][i] = zero; } - ParticleReal aerosol_mass_total = zero; - for (int ctr = 0; ctr < n_aerosols; ctr++) { - aerosol_mass_total += aerosol_mass_ptrs[ctr][i]; - } + ptrs.mult_ptr[i] = multiplicity; - auto par_radius = Real(1.0e-15); - auto cond_mass = (Real(4.0)/three)*PI - * par_radius*par_radius*par_radius - * rho_w; - radius_ptr[i] = par_radius; - mass_ptr[i] = cond_mass + aerosol_mass_total; - mult_ptr[i] = multiplicity; + SuperDropletPC::updateParticleAttributes( + i, ptrs.radius_ptr, ptrs.mass_ptr, ctx.idx_water, ctx.rho_water, + ctx.num_species, ctx.num_aerosols, ptrs.sp_sol_arr, ptrs.ae_sol_arr, + ptrs.sp_mass_ptrs, ptrs.ae_mass_ptrs, ptrs.sp_rho_arr, ptrs.ae_rho_arr); } } @@ -213,12 +157,10 @@ void SuperDropletPC::applyBoundaryTreatment ( int a_lev, } } - - // Update z-coordinate carried by the particle }); Gpu::synchronize(); num_deactivated_particles += *(deactivated_particles.copyToHost()); - } + }); // end forEachParticleTile ParallelDescriptor::ReduceLongSum( &num_deactivated_particles, 1 ); Print() << "SuperDropletPC(" << m_name << "): " diff --git a/Source/Particles/ERF_SuperDropletPCCoalescence.H b/Source/Particles/ERF_SuperDropletPCCoalescence.H index 1a171e4a7..9a44db2c3 100644 --- a/Source/Particles/ERF_SuperDropletPCCoalescence.H +++ b/Source/Particles/ERF_SuperDropletPCCoalescence.H @@ -22,9 +22,9 @@ struct CollisionKernel RT golovin ( const RT a_r1, /*!< radius of droplet 1 */ const RT a_r2 /*!< radius of droplet 2 */ ) const noexcept { - RT b = amrex::Real(1.5e03); - auto X_1 = (amrex::Real(4.0)*PI/three) * a_r1*a_r1*a_r1; - auto X_2 = (amrex::Real(4.0)*PI/three) * a_r2*a_r2*a_r2; + RT b = RT(1.5e03); + auto X_1 = RT(four_thirds_pi) * a_r1*a_r1*a_r1; + auto X_2 = RT(four_thirds_pi) * a_r2*a_r2*a_r2; return b * (X_1 + X_2); } @@ -49,7 +49,7 @@ struct CollisionKernel const RT* const a_v2 ) const noexcept { auto p = std::min(a_r1,a_r2) / std::max(a_r1,a_r2); - auto E = myhalf*p*p / ((one+p)*(one+p)); + auto E = RT(0.5)*p*p / ((RT(1.0)+p)*(RT(1.0)+p)); RT dv = 0; for (int d = 0; d < SPACEDIM; d++) { dv += (a_v1[d]-a_v2[d])*(a_v1[d]-a_v2[d]); @@ -67,7 +67,7 @@ struct CollisionKernel * effective for smaller cloud droplets where hydrodynamic effects are important. * * Reference: Long, A. B., 1974: Solutions to the droplet collection equation - * for polynomial kernels. J. Atmos. Sci., 31, 1040-amrex::Real(1052.) + * for polynomial kernels. J. Atmos. Sci., 31, 1040-RT(1052.) * * \param[in] a_r1 Radius of particle 1 (m) * \param[in] a_r2 Radius of particle 2 (m) @@ -85,9 +85,9 @@ struct CollisionKernel auto r2 = std::min( a_r1, a_r2 ); // small auto sumr = r1 + r2; - RT c_rate = one; - if( r1 <= amrex::Real(5.E-5) ) { - c_rate = amrex::Real(4.5E+8) * ( r1*r1 ) * ( one - amrex::Real(3.E-6)/(std::max(amrex::Real(3.01E-6),r1)) ); + RT c_rate = RT(1.0); + if( r1 <= RT(5.E-5) ) { + c_rate = RT(4.5E+8) * ( r1*r1 ) * ( RT(1.0) - RT(3.E-6)/(std::max(RT(3.01E-6),r1)) ); } c_rate *= (PI*sumr*sumr); @@ -106,21 +106,21 @@ struct CollisionKernel AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE RT Halls_data_r0col (const int a_i /*!< index */) const noexcept { - static const RT vec[15] = { amrex::Real(6.0), - amrex::Real(8.0), - amrex::Real(10.0), - amrex::Real(15.0), - amrex::Real(20.0), - amrex::Real(25.0), - amrex::Real(30.0), - amrex::Real(40.0), - amrex::Real(50.0), - amrex::Real(60.0), - amrex::Real(70.0), - amrex::Real(100.0), - amrex::Real(150.0), - amrex::Real(200.0), - amrex::Real(300.0)}; + static const RT vec[15] = { RT(6.0), + RT(8.0), + RT(10.0), + RT(15.0), + RT(20.0), + RT(25.0), + RT(30.0), + RT(40.0), + RT(50.0), + RT(60.0), + RT(70.0), + RT(100.0), + RT(150.0), + RT(200.0), + RT(300.0)}; return vec[a_i]; } @@ -130,7 +130,7 @@ struct CollisionKernel AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int Halls_data_r0col ( const RT a_r /*!< radius */) const noexcept { - auto radius_microns = a_r*amrex::Real(1.0e6); // [m] -> [mu-m] + auto radius_microns = a_r*RT(1.0e6); // [m] -> [mu-m] int i; for (i = 0; i < 15; i++) { if (radius_microns <= Halls_data_r0col(i)) { return i; } @@ -144,27 +144,27 @@ struct CollisionKernel AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE RT Halls_data_ratcol ( const int a_i /*!< index */) const noexcept { - static const RT vec[21] = {amrex::Real(0.00), - amrex::Real(0.05), - amrex::Real(0.10), - amrex::Real(0.15), - amrex::Real(0.20), + static const RT vec[21] = {RT(0.00), + RT(0.05), + RT(0.10), + RT(0.15), + RT(0.20), fourth, - amrex::Real(0.30), - amrex::Real(0.35), - amrex::Real(0.40), - amrex::Real(0.45), - amrex::Real(0.50), - amrex::Real(0.55), - amrex::Real(0.60), - amrex::Real(0.65), - amrex::Real(0.70), - amrex::Real(0.75), - amrex::Real(0.80), - amrex::Real(0.85), - amrex::Real(0.90), - amrex::Real(0.95), - amrex::Real(1.00)}; + RT(0.30), + RT(0.35), + RT(0.40), + RT(0.45), + RT(0.50), + RT(0.55), + RT(0.60), + RT(0.65), + RT(0.70), + RT(0.75), + RT(0.80), + RT(0.85), + RT(0.90), + RT(0.95), + RT(1.00)}; return vec[a_i]; } @@ -189,27 +189,27 @@ struct CollisionKernel const int a_j /*!< column index */) const noexcept { const RT ecoll[21][15] = - { {amrex::Real(0.0010),amrex::Real(0.0010),amrex::Real(0.0010),amrex::Real(0.0010),amrex::Real(0.0010),amrex::Real(0.0010),amrex::Real(0.0010),amrex::Real(0.0010),amrex::Real(0.0010),amrex::Real(0.0010),amrex::Real(0.0010),amrex::Real(0.0010),amrex::Real(0.0010),amrex::Real(0.0010),amrex::Real(0.0010)}, - {amrex::Real(0.0030),amrex::Real(0.0030),amrex::Real(0.0030),amrex::Real(0.0040),amrex::Real(0.0050),amrex::Real(0.0050),amrex::Real(0.0050),amrex::Real(0.0100),amrex::Real(0.1000),amrex::Real(0.0500),amrex::Real(0.2000),amrex::Real(0.5000),amrex::Real(0.7700),amrex::Real(0.8700),amrex::Real(0.9700)}, - {amrex::Real(0.0070),amrex::Real(0.0070),amrex::Real(0.0070),amrex::Real(0.0080),amrex::Real(0.0090),amrex::Real(0.0100),amrex::Real(0.0100),amrex::Real(0.0700),amrex::Real(0.4000),amrex::Real(0.4300),amrex::Real(0.5800),amrex::Real(0.7900),amrex::Real(0.9300),amrex::Real(0.9600),amrex::Real(1.0000)}, - {amrex::Real(0.0090),amrex::Real(0.0090),amrex::Real(0.0090),amrex::Real(0.0120),amrex::Real(0.0150),amrex::Real(0.0100),amrex::Real(0.0200),amrex::Real(0.2800),amrex::Real(0.6000),amrex::Real(0.6400),amrex::Real(0.7500),amrex::Real(0.9100),amrex::Real(0.9700),amrex::Real(0.9800),amrex::Real(1.0000)}, - {amrex::Real(0.0140),amrex::Real(0.0140),amrex::Real(0.0140),amrex::Real(0.0150),amrex::Real(0.0160),amrex::Real(0.0300),amrex::Real(0.0600),amrex::Real(0.5000),amrex::Real(0.7000),amrex::Real(0.7700),amrex::Real(0.8400),amrex::Real(0.9500),amrex::Real(0.9700),amrex::Real(1.0000),amrex::Real(1.0000)}, - {amrex::Real(0.0170),amrex::Real(0.0170),amrex::Real(0.0170),amrex::Real(0.0200),amrex::Real(0.0220),amrex::Real(0.0600),amrex::Real(0.1000),amrex::Real(0.6200),amrex::Real(0.7800),amrex::Real(0.8400),amrex::Real(0.8800),amrex::Real(0.9500),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000)}, - {amrex::Real(0.0300),amrex::Real(0.0300),amrex::Real(0.0240),amrex::Real(0.0220),amrex::Real(0.0320),amrex::Real(0.0620),amrex::Real(0.2000),amrex::Real(0.6800),amrex::Real(0.8300),amrex::Real(0.8700),amrex::Real(0.9000),amrex::Real(0.9500),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000)}, - {amrex::Real(0.0250),amrex::Real(0.0250),amrex::Real(0.0250),amrex::Real(0.0360),amrex::Real(0.0430),amrex::Real(0.1300),amrex::Real(0.2700),amrex::Real(0.7400),amrex::Real(0.8600),amrex::Real(0.8900),amrex::Real(0.9200),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000)}, - {amrex::Real(0.0270),amrex::Real(0.0270),amrex::Real(0.0270),amrex::Real(0.0400),amrex::Real(0.0520),amrex::Real(0.2000),amrex::Real(0.4000),amrex::Real(0.7800),amrex::Real(0.8800),amrex::Real(0.9000),amrex::Real(0.9400),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000)}, - {amrex::Real(0.0300),amrex::Real(0.0300),amrex::Real(0.0300),amrex::Real(0.0470),amrex::Real(0.0640),amrex::Real(0.2500),amrex::Real(0.5000),amrex::Real(0.8000),amrex::Real(0.9000),amrex::Real(0.9100),amrex::Real(0.9500),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000)}, - {amrex::Real(0.0400),amrex::Real(0.0400),amrex::Real(0.0330),amrex::Real(0.0370),amrex::Real(0.0680),amrex::Real(0.2400),amrex::Real(0.5500),amrex::Real(0.8000),amrex::Real(0.9000),amrex::Real(0.9100),amrex::Real(0.9500),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000)}, - {amrex::Real(0.0350),amrex::Real(0.0350),amrex::Real(0.0350),amrex::Real(0.0550),amrex::Real(0.0790),amrex::Real(0.2900),amrex::Real(0.5800),amrex::Real(0.8000),amrex::Real(0.9000),amrex::Real(0.9100),amrex::Real(0.9500),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000)}, - {amrex::Real(0.0370),amrex::Real(0.0370),amrex::Real(0.0370),amrex::Real(0.0620),amrex::Real(0.0820),amrex::Real(0.2900),amrex::Real(0.5900),amrex::Real(0.7800),amrex::Real(0.9000),amrex::Real(0.9100),amrex::Real(0.9500),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000)}, - {amrex::Real(0.0370),amrex::Real(0.0370),amrex::Real(0.0370),amrex::Real(0.0600),amrex::Real(0.0800),amrex::Real(0.2900),amrex::Real(0.5800),amrex::Real(0.7700),amrex::Real(0.8900),amrex::Real(0.9100),amrex::Real(0.9500),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000)}, - {amrex::Real(0.0370),amrex::Real(0.0370),amrex::Real(0.0370),amrex::Real(0.0410),amrex::Real(0.0750),amrex::Real(0.2500),amrex::Real(0.5400),amrex::Real(0.7600),amrex::Real(0.8800),amrex::Real(0.9200),amrex::Real(0.9500),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000)}, - {amrex::Real(0.0370),amrex::Real(0.0370),amrex::Real(0.0370),amrex::Real(0.0520),amrex::Real(0.0670),amrex::Real(0.2500),amrex::Real(0.5100),amrex::Real(0.7700),amrex::Real(0.8800),amrex::Real(0.9300),amrex::Real(0.9700),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000)}, - {amrex::Real(0.0370),amrex::Real(0.0370),amrex::Real(0.0370),amrex::Real(0.0470),amrex::Real(0.0570),amrex::Real(0.2500),amrex::Real(0.4900),amrex::Real(0.7700),amrex::Real(0.8900),amrex::Real(0.9500),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000),amrex::Real(1.0000)}, - {amrex::Real(0.0360),amrex::Real(0.0360),amrex::Real(0.0360),amrex::Real(0.0420),amrex::Real(0.0480),amrex::Real(0.2300),amrex::Real(0.4700),amrex::Real(0.7800),amrex::Real(0.9200),amrex::Real(1.0000),amrex::Real(1.0200),amrex::Real(1.0200),amrex::Real(1.0200),amrex::Real(1.0200),amrex::Real(1.0200)}, - {amrex::Real(0.0400),amrex::Real(0.0400),amrex::Real(0.0350),amrex::Real(0.0330),amrex::Real(0.0400),amrex::Real(0.1120),amrex::Real(0.4500),amrex::Real(0.7900),amrex::Real(1.0100),amrex::Real(1.0300),amrex::Real(1.0400),amrex::Real(1.0400),amrex::Real(1.0400),amrex::Real(1.0400),amrex::Real(1.0400)}, - {amrex::Real(0.0330),amrex::Real(0.0330),amrex::Real(0.0330),amrex::Real(0.0330),amrex::Real(0.0330),amrex::Real(0.1190),amrex::Real(0.4700),amrex::Real(0.9500),amrex::Real(1.3000),amrex::Real(1.7000),amrex::Real(2.3000),amrex::Real(2.3000),amrex::Real(2.3000),amrex::Real(2.3000),amrex::Real(2.3000)}, - {amrex::Real(0.0270),amrex::Real(0.0270),amrex::Real(0.0270),amrex::Real(0.0270),amrex::Real(0.0270),amrex::Real(0.1250),amrex::Real(0.5200),amrex::Real(1.4000),amrex::Real(2.3000),amrex::Real(3.0000),amrex::Real(4.0000),amrex::Real(4.0000),amrex::Real(4.0000),amrex::Real(4.0000),amrex::Real(4.0000)} }; + { {RT(0.0010),RT(0.0010),RT(0.0010),RT(0.0010),RT(0.0010),RT(0.0010),RT(0.0010),RT(0.0010),RT(0.0010),RT(0.0010),RT(0.0010),RT(0.0010),RT(0.0010),RT(0.0010),RT(0.0010)}, + {RT(0.0030),RT(0.0030),RT(0.0030),RT(0.0040),RT(0.0050),RT(0.0050),RT(0.0050),RT(0.0100),RT(0.1000),RT(0.0500),RT(0.2000),RT(0.5000),RT(0.7700),RT(0.8700),RT(0.9700)}, + {RT(0.0070),RT(0.0070),RT(0.0070),RT(0.0080),RT(0.0090),RT(0.0100),RT(0.0100),RT(0.0700),RT(0.4000),RT(0.4300),RT(0.5800),RT(0.7900),RT(0.9300),RT(0.9600),RT(1.0000)}, + {RT(0.0090),RT(0.0090),RT(0.0090),RT(0.0120),RT(0.0150),RT(0.0100),RT(0.0200),RT(0.2800),RT(0.6000),RT(0.6400),RT(0.7500),RT(0.9100),RT(0.9700),RT(0.9800),RT(1.0000)}, + {RT(0.0140),RT(0.0140),RT(0.0140),RT(0.0150),RT(0.0160),RT(0.0300),RT(0.0600),RT(0.5000),RT(0.7000),RT(0.7700),RT(0.8400),RT(0.9500),RT(0.9700),RT(1.0000),RT(1.0000)}, + {RT(0.0170),RT(0.0170),RT(0.0170),RT(0.0200),RT(0.0220),RT(0.0600),RT(0.1000),RT(0.6200),RT(0.7800),RT(0.8400),RT(0.8800),RT(0.9500),RT(1.0000),RT(1.0000),RT(1.0000)}, + {RT(0.0300),RT(0.0300),RT(0.0240),RT(0.0220),RT(0.0320),RT(0.0620),RT(0.2000),RT(0.6800),RT(0.8300),RT(0.8700),RT(0.9000),RT(0.9500),RT(1.0000),RT(1.0000),RT(1.0000)}, + {RT(0.0250),RT(0.0250),RT(0.0250),RT(0.0360),RT(0.0430),RT(0.1300),RT(0.2700),RT(0.7400),RT(0.8600),RT(0.8900),RT(0.9200),RT(1.0000),RT(1.0000),RT(1.0000),RT(1.0000)}, + {RT(0.0270),RT(0.0270),RT(0.0270),RT(0.0400),RT(0.0520),RT(0.2000),RT(0.4000),RT(0.7800),RT(0.8800),RT(0.9000),RT(0.9400),RT(1.0000),RT(1.0000),RT(1.0000),RT(1.0000)}, + {RT(0.0300),RT(0.0300),RT(0.0300),RT(0.0470),RT(0.0640),RT(0.2500),RT(0.5000),RT(0.8000),RT(0.9000),RT(0.9100),RT(0.9500),RT(1.0000),RT(1.0000),RT(1.0000),RT(1.0000)}, + {RT(0.0400),RT(0.0400),RT(0.0330),RT(0.0370),RT(0.0680),RT(0.2400),RT(0.5500),RT(0.8000),RT(0.9000),RT(0.9100),RT(0.9500),RT(1.0000),RT(1.0000),RT(1.0000),RT(1.0000)}, + {RT(0.0350),RT(0.0350),RT(0.0350),RT(0.0550),RT(0.0790),RT(0.2900),RT(0.5800),RT(0.8000),RT(0.9000),RT(0.9100),RT(0.9500),RT(1.0000),RT(1.0000),RT(1.0000),RT(1.0000)}, + {RT(0.0370),RT(0.0370),RT(0.0370),RT(0.0620),RT(0.0820),RT(0.2900),RT(0.5900),RT(0.7800),RT(0.9000),RT(0.9100),RT(0.9500),RT(1.0000),RT(1.0000),RT(1.0000),RT(1.0000)}, + {RT(0.0370),RT(0.0370),RT(0.0370),RT(0.0600),RT(0.0800),RT(0.2900),RT(0.5800),RT(0.7700),RT(0.8900),RT(0.9100),RT(0.9500),RT(1.0000),RT(1.0000),RT(1.0000),RT(1.0000)}, + {RT(0.0370),RT(0.0370),RT(0.0370),RT(0.0410),RT(0.0750),RT(0.2500),RT(0.5400),RT(0.7600),RT(0.8800),RT(0.9200),RT(0.9500),RT(1.0000),RT(1.0000),RT(1.0000),RT(1.0000)}, + {RT(0.0370),RT(0.0370),RT(0.0370),RT(0.0520),RT(0.0670),RT(0.2500),RT(0.5100),RT(0.7700),RT(0.8800),RT(0.9300),RT(0.9700),RT(1.0000),RT(1.0000),RT(1.0000),RT(1.0000)}, + {RT(0.0370),RT(0.0370),RT(0.0370),RT(0.0470),RT(0.0570),RT(0.2500),RT(0.4900),RT(0.7700),RT(0.8900),RT(0.9500),RT(1.0000),RT(1.0000),RT(1.0000),RT(1.0000),RT(1.0000)}, + {RT(0.0360),RT(0.0360),RT(0.0360),RT(0.0420),RT(0.0480),RT(0.2300),RT(0.4700),RT(0.7800),RT(0.9200),RT(1.0000),RT(1.0200),RT(1.0200),RT(1.0200),RT(1.0200),RT(1.0200)}, + {RT(0.0400),RT(0.0400),RT(0.0350),RT(0.0330),RT(0.0400),RT(0.1120),RT(0.4500),RT(0.7900),RT(1.0100),RT(1.0300),RT(1.0400),RT(1.0400),RT(1.0400),RT(1.0400),RT(1.0400)}, + {RT(0.0330),RT(0.0330),RT(0.0330),RT(0.0330),RT(0.0330),RT(0.1190),RT(0.4700),RT(0.9500),RT(1.3000),RT(1.7000),RT(2.3000),RT(2.3000),RT(2.3000),RT(2.3000),RT(2.3000)}, + {RT(0.0270),RT(0.0270),RT(0.0270),RT(0.0270),RT(0.0270),RT(0.1250),RT(0.5200),RT(1.4000),RT(2.3000),RT(3.0000),RT(4.0000),RT(4.0000),RT(4.0000),RT(4.0000),RT(4.0000)} }; return ecoll[a_i][a_j]; } @@ -241,11 +241,11 @@ struct CollisionKernel RT ek = (one-q) * Halls_data_ecoll(iqq-1,14) + q * Halls_data_ecoll(iqq ,14); - c_rate = std::min( ek, amrex::Real(1.0) ); + c_rate = std::min( ek, RT(1.0) ); } else if ((irr >= 1) && (irr < 15)) { - RT p = ( r1*amrex::Real(1.0e6) - Halls_data_r0col(irr-1) ) + RT p = ( r1*RT(1.0e6) - Halls_data_r0col(irr-1) ) / ( Halls_data_r0col(irr) - Halls_data_r0col(irr-1) ); RT q = ( ratio - Halls_data_ratcol(iqq-1) ) @@ -295,24 +295,24 @@ struct CollisionKernel RT t_crs = a_temperature; // [K] // dynamic viscosity [Pa*s] (Pruppacher & Klett,1997) - RT T_degC = t_crs - amrex::Real(273.15); // [K] => [degC] + RT T_degC = t_crs - RT(273.15); // [K] => [degC] RT vis_crs = zero; if( T_degC >= zero ) { - vis_crs = ( amrex::Real(1.7180) + amrex::Real(4.9E-3)*T_degC ) * amrex::Real(1.E-5); + vis_crs = ( RT(1.7180) + RT(4.9E-3)*T_degC ) * RT(1.E-5); } else { - vis_crs = ( amrex::Real(1.7180) + amrex::Real(4.9E-3)*T_degC -amrex::Real(1.2E-5)*T_degC*T_degC ) * amrex::Real(1.E-5); + vis_crs = ( RT(1.7180) + RT(4.9E-3)*T_degC -RT(1.2E-5)*T_degC*T_degC ) * RT(1.E-5); } //air mean free path [m] - RT lambda_crs = (two*vis_crs) / (p_crs*std::sqrt(amrex::Real(8.0)*mwdair*amrex::Real(1.E-3)/(PI*R_d*t_crs))); + RT lambda_crs = (two*vis_crs) / (p_crs*std::sqrt(RT(8.0)*mwdair*RT(1.E-3)/(PI*R_d*t_crs))); // temporary var RT dtmp = zero; // slip correction of droplets [-] - dtmp = amrex::Real(1.2570) + amrex::Real(0.40) * exp(-amrex::Real(0.550)*diameter_1/lambda_crs); + dtmp = RT(1.2570) + RT(0.40) * exp(-RT(0.550)*diameter_1/lambda_crs); RT slip_corr_1 = one + (two*lambda_crs*dtmp)/(diameter_1); - dtmp = amrex::Real(1.2570) + amrex::Real(0.40) * exp(-amrex::Real(0.550)*diameter_2/lambda_crs); + dtmp = RT(1.2570) + RT(0.40) * exp(-RT(0.550)*diameter_2/lambda_crs); RT slip_corr_2 = one + (two*lambda_crs*dtmp)/(diameter_2); // diffusion term [m*m/s] @@ -321,21 +321,21 @@ struct CollisionKernel RT diff_term_2 = dtmp * (slip_corr_2/diameter_2); // velocity term [m/s] - dtmp = (amrex::Real(8.0)*boltz*t_crs)/PI; + dtmp = (RT(8.0)*boltz*t_crs)/PI; RT vel_term_1 = std::sqrt(dtmp/a_mass_1); RT vel_term_2 = std::sqrt(dtmp/a_mass_2); // mean free path of droplets [m] - dtmp = amrex::Real(8.0)/PI; + dtmp = RT(8.0)/PI; RT lambda_1 = dtmp * (diff_term_1/vel_term_1); RT lambda_2 = dtmp * (diff_term_2/vel_term_2); //length term [m] dtmp = (diameter_1+lambda_1)*(diameter_1+lambda_1)*(diameter_1+lambda_1) - - std::exp(amrex::Real(1.5)*std::log(diameter_1*diameter_1+lambda_1*lambda_1)); + - std::exp(RT(1.5)*std::log(diameter_1*diameter_1+lambda_1*lambda_1)); RT length_term_1 = dtmp/(three*diameter_1*lambda_1) - diameter_1; dtmp = (diameter_2+lambda_2)*(diameter_2+lambda_2)*(diameter_2+lambda_2) - - std::exp(amrex::Real(1.5)*std::log(diameter_2*diameter_2+lambda_2*lambda_2)); + - std::exp(RT(1.5)*std::log(diameter_2*diameter_2+lambda_2*lambda_2)); RT length_term_2 = dtmp/(three*diameter_2*lambda_2) - diameter_2; // Brownian Coagulation Coefficient K12 [m3/s] @@ -344,7 +344,7 @@ struct CollisionKernel RT sumc = std::sqrt( vel_term_1*vel_term_1 + vel_term_2*vel_term_2 ); RT sumg = std::sqrt( two*length_term_1*length_term_1 + two*length_term_2*length_term_2 ); - dtmp = sumdia/(sumdia+two*sumg) + (amrex::Real(8.0)*sumd)/(sumdia*sumc); + dtmp = sumdia/(sumdia+two*sumg) + (RT(8.0)*sumd)/(sumdia*sumc); RT k12 = two*PI * sumdia*sumd/dtmp; return k12; diff --git a/Source/Particles/ERF_SuperDropletPCDefinitions.H b/Source/Particles/ERF_SuperDropletPCDefinitions.H index c66a6b4e8..2687621e4 100644 --- a/Source/Particles/ERF_SuperDropletPCDefinitions.H +++ b/Source/Particles/ERF_SuperDropletPCDefinitions.H @@ -3,10 +3,12 @@ #ifdef ERF_USE_PARTICLES +#include #include #include #include #include "ERF_IndexDefines.H" +#include "ERF_Constants.H" #include "ERFPC.H" #include "ERF_SDInitialization.H" @@ -23,6 +25,9 @@ namespace SDPCDefn { /*! \brief Real-type struct-of-array attributes of a super-droplet */ using SuperDropletsRealIdxSoA = ERFParticlesRealIdxSoA; + /*! \brief Int-type struct-of-array attributes of a super-droplet */ + using SuperDropletsIntIdxSoA = ERFParticlesIntIdxSoA; + /*! \brief Int-type struct-of-array attributes added during runtime */ struct SuperDropletsIntIdxSoA_RT { @@ -47,55 +52,44 @@ namespace SDPCDefn { ncomps }; }; - - /*! \brief Int-type struct-of-array attributes of a super-droplet */ - using SuperDropletsIntIdxSoA = ERFParticlesIntIdxSoA; } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -static int idx_a(const int a_i, /*!< Species index */ - const int a_num_a, /*!< Number of aerosols */ - const int a_num_s /*!< Number of species */) +static int ridx_a(const int a_i, /*!< Species index */ + const int a_num_a, /*!< Number of aerosols */ + const int a_num_s /*!< Number of species */) { amrex::ignore_unused(a_num_a); amrex::ignore_unused(a_num_s); - return SDPCDefn::SuperDropletsRealIdxSoA::ncomps - + SDPCDefn::SuperDropletsRealIdxSoA_RT::ncomps + return SDPCDefn::SuperDropletsRealIdxSoA_RT::ncomps + a_i; } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -static int idx_s(const int a_i, /*!< Species index */ - const int a_num_a, /*!< Number of aerosols */ - const int a_num_s /*!< Number of species */) +static int ridx_s(const int a_i, /*!< Species index */ + const int a_num_a, /*!< Number of aerosols */ + const int a_num_s /*!< Number of species */) { amrex::ignore_unused(a_num_s); - return SDPCDefn::SuperDropletsRealIdxSoA::ncomps - + SDPCDefn::SuperDropletsRealIdxSoA_RT::ncomps + return SDPCDefn::SuperDropletsRealIdxSoA_RT::ncomps + a_num_a + a_i; } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -static int ridx_a(const int a_i, /*!< Species index */ +static int idx_a(const int a_i, /*!< Species index */ const int a_num_a, /*!< Number of aerosols */ const int a_num_s /*!< Number of species */) { - amrex::ignore_unused(a_num_a); - amrex::ignore_unused(a_num_s); - return SDPCDefn::SuperDropletsRealIdxSoA_RT::ncomps - + a_i; + return SDPCDefn::SuperDropletsRealIdxSoA::ncomps + ridx_a(a_i,a_num_a,a_num_s); } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -static int ridx_s(const int a_i, /*!< Species index */ +static int idx_s(const int a_i, /*!< Species index */ const int a_num_a, /*!< Number of aerosols */ const int a_num_s /*!< Number of species */) { - amrex::ignore_unused(a_num_s); - return SDPCDefn::SuperDropletsRealIdxSoA_RT::ncomps - + a_num_a - + a_i; + return SDPCDefn::SuperDropletsRealIdxSoA::ncomps + ridx_s(a_i,a_num_a,a_num_s); } /*! \brief Compute total mass of a droplet from its components */ @@ -155,7 +149,7 @@ static amrex::ParticleReal SD_effective_radius( const int a_idx, /*!< Particle i if (m_p > zero) { rho_p /= m_p; } else { rho_p = one; } auto m_t = m_w + m_s + (a_rho_w/rho_p)*m_p; - auto r_eff = std::cbrt(m_t / (amrex::Real(4.0)/three*PI*a_rho_w)); + auto r_eff = std::cbrt(m_t / (four_thirds_pi*a_rho_w)); return r_eff; } @@ -187,16 +181,16 @@ static amrex::ParticleReal SD_dry_radius( const int a_idx, /*!< Particle index * } if (m_p > zero) { rho_p /= m_p; } else { rho_p = one; } - auto r_dry = std::cbrt(m_p / (amrex::Real(4.0)/three*PI*rho_p)); + auto r_dry = std::cbrt(m_p / (four_thirds_pi*rho_p)); return r_dry; } /*! \brief Types of terminal velocity models */ -enum struct SDTerminalVelocityType { +AMREX_ENUM(SDTerminalVelocityType, RogersYau, AtlasUlbrich, CloudRainShima -}; +); /*! \brief Types of coalescence kernels */ enum struct SDCoalescenceKernelType { @@ -213,4 +207,3 @@ enum struct SDMassChangeTIMethod { #endif #endif - diff --git a/Source/Particles/ERF_SuperDropletPCDiagnostics.cpp b/Source/Particles/ERF_SuperDropletPCDiagnostics.cpp index f07ffe36a..1bc334666 100644 --- a/Source/Particles/ERF_SuperDropletPCDiagnostics.cpp +++ b/Source/Particles/ERF_SuperDropletPCDiagnostics.cpp @@ -2,6 +2,7 @@ #include #include #include +#include #include "ERF_Constants.H" #include "ERF_SuperDropletPC.H" @@ -75,163 +76,154 @@ void SuperDropletPC::Diagnostics( const int a_iter, auto num_total_particles = TotalNumberOfParticles(); auto num_superdroplets = static_cast(NumSuperDroplets()); - auto min_par_radius = ReduceMin( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { return ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::radius][i]; } ); - - auto max_par_radius = ReduceMax( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { return ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::radius][i]; } ); - - auto avg_par_radius = ReduceSum( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { - auto n = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; - auto r = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::radius][i]; - return n*r; - } ); - - auto min_multiplic = ReduceMin( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { return ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; } ); - - auto max_multiplic = ReduceMax( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { return ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; } ); + // Number of base attributes to reduce (radius, multiplicity, mass, vx, vy, vz, term_vel) + constexpr int NUM_BASE_ATTRS = 7; +#ifdef ERF_USE_ML_UPHYS_DIAGNOSTICS + constexpr int NUM_ATTRS = NUM_BASE_ATTRS + 1; // + cond_tendency +#else + constexpr int NUM_ATTRS = NUM_BASE_ATTRS; +#endif - auto avg_multiplic = ReduceSum( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { return ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; } ); + // Fused reduction: compute min, max, and weighted sum for all attributes in a single pass + using ReduceTuple = GpuTuple; - auto min_par_mass = ReduceMin( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { return ptd.m_rdata[SuperDropletsRealIdxSoA::mass][i]; } ); + ReduceOps reduce_ops; - auto max_par_mass = ReduceMax( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { return ptd.m_rdata[SuperDropletsRealIdxSoA::mass][i]; } ); + auto r = ReduceData(reduce_ops); + + for (int lev = 0; lev <= finestLevel(); ++lev) { + for (ParConstIterType pti(*this, lev); pti.isValid(); ++pti) { + const auto& ptd = pti.GetParticleTile().getConstParticleTileData(); + const int np = pti.numParticles(); + + reduce_ops.eval(np, r, + [=] AMREX_GPU_DEVICE (int i) -> ReduceTuple + { + const Real n = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; + const Real radius = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::radius][i]; + const Real mass = ptd.m_rdata[SuperDropletsRealIdxSoA::mass][i]; + const Real vx = ptd.m_rdata[SuperDropletsRealIdxSoA::vx][i]; + const Real vy = ptd.m_rdata[SuperDropletsRealIdxSoA::vy][i]; + const Real vz = ptd.m_rdata[SuperDropletsRealIdxSoA::vz][i]; + const Real term_vel = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::term_vel][i]; +#ifdef ERF_USE_ML_UPHYS_DIAGNOSTICS + const Real cond_t = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::cond_tendency][i]; +#endif - auto avg_par_mass = ReduceSum( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { - auto n = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; - auto m = ptd.m_rdata[SuperDropletsRealIdxSoA::mass][i]; - return n*m; - } ); + return { + // Min values + radius, n, mass, vx, vy, vz, term_vel, + // Max values + radius, n, mass, vx, vy, vz, term_vel, + // Weighted sums (n * value) + n * radius, n, n * mass, n * vx, n * vy, n * vz, n * term_vel, + // Multiplicity sum + n +#ifdef ERF_USE_ML_UPHYS_DIAGNOSTICS + , cond_t, cond_t, n * cond_t // min, max, sum for cond_tendency +#endif + }; + }); + } + } - auto min_par_vx = ReduceMin( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { return ptd.m_rdata[SuperDropletsRealIdxSoA::vx][i]; } ); - - auto max_par_vx = ReduceMax( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { return ptd.m_rdata[SuperDropletsRealIdxSoA::vx][i]; } ); - - auto avg_par_vx = ReduceSum( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { - auto n = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; - auto m = ptd.m_rdata[SuperDropletsRealIdxSoA::vx][i]; - return n*m; - } ); - - auto min_par_vy = ReduceMin( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { return ptd.m_rdata[SuperDropletsRealIdxSoA::vy][i]; } ); - - auto max_par_vy = ReduceMax( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { return ptd.m_rdata[SuperDropletsRealIdxSoA::vy][i]; } ); - - auto avg_par_vy = ReduceSum( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { - auto n = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; - auto m = ptd.m_rdata[SuperDropletsRealIdxSoA::vy][i]; - return n*m; - } ); - - auto min_par_vz = ReduceMin( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { return ptd.m_rdata[SuperDropletsRealIdxSoA::vz][i]; } ); - - auto max_par_vz = ReduceMax( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { return ptd.m_rdata[SuperDropletsRealIdxSoA::vz][i]; } ); - - auto avg_par_vz = ReduceSum( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { - auto n = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; - auto m = ptd.m_rdata[SuperDropletsRealIdxSoA::vz][i]; - return n*m; - } ); - - auto min_term_v = ReduceMin( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { return ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::term_vel][i]; } ); - - auto max_term_v = ReduceMax( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { return ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::term_vel][i]; } ); - - auto avg_term_v = ReduceSum( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { - auto n = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; - auto m = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::term_vel][i]; - return n*m; - } ); + auto hv = r.value(reduce_ops); + + // Extract results from tuple + Real min_par_radius = amrex::get<0>(hv); + Real min_multiplic = amrex::get<1>(hv); + Real min_par_mass = amrex::get<2>(hv); + Real min_par_vx = amrex::get<3>(hv); + Real min_par_vy = amrex::get<4>(hv); + Real min_par_vz = amrex::get<5>(hv); + Real min_term_v = amrex::get<6>(hv); + + Real max_par_radius = amrex::get<7>(hv); + Real max_multiplic = amrex::get<8>(hv); + Real max_par_mass = amrex::get<9>(hv); + Real max_par_vx = amrex::get<10>(hv); + Real max_par_vy = amrex::get<11>(hv); + Real max_par_vz = amrex::get<12>(hv); + Real max_term_v = amrex::get<13>(hv); + + Real avg_par_radius = amrex::get<14>(hv); + Real avg_multiplic = amrex::get<15>(hv); + Real avg_par_mass = amrex::get<16>(hv); + Real avg_par_vx = amrex::get<17>(hv); + Real avg_par_vy = amrex::get<18>(hv); + Real avg_par_vz = amrex::get<19>(hv); + Real avg_term_v = amrex::get<20>(hv); + + // Note: amrex::get<21>(hv) is the multiplicity sum, used for avg_multiplic normalization #ifdef ERF_USE_ML_UPHYS_DIAGNOSTICS - auto min_cond_t = ReduceMin( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { return ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::cond_tendency][i]; } ); - - auto max_cond_t = ReduceMax( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { return ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::cond_tendency][i]; } ); - - auto avg_cond_t = ReduceSum( *this, - [=] AMREX_GPU_HOST_DEVICE (const SDTDType& ptd, const int i) -> Real - { - auto n = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; - auto m = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::cond_tendency][i]; - return n*m; - } ); + Real min_cond_t = amrex::get<22>(hv); + Real max_cond_t = amrex::get<23>(hv); + Real avg_cond_t = amrex::get<24>(hv); #endif - ParallelDescriptor::ReduceRealMin(&min_par_mass,1); - ParallelDescriptor::ReduceRealMin(&min_par_radius,1); - ParallelDescriptor::ReduceRealMin(&min_multiplic,1); - ParallelDescriptor::ReduceRealMin(&min_par_vx,1); - ParallelDescriptor::ReduceRealMin(&min_par_vy,1); - ParallelDescriptor::ReduceRealMin(&min_par_vz,1); - ParallelDescriptor::ReduceRealMin(&min_term_v,1); + // MPI reductions - batch into arrays for efficiency + Real min_vals[NUM_ATTRS] = {min_par_radius, min_multiplic, min_par_mass, + min_par_vx, min_par_vy, min_par_vz, min_term_v #ifdef ERF_USE_ML_UPHYS_DIAGNOSTICS - ParallelDescriptor::ReduceRealMin(&min_cond_t,1); + , min_cond_t #endif - - ParallelDescriptor::ReduceRealMax(&max_par_mass,1); - ParallelDescriptor::ReduceRealMax(&max_par_radius,1); - ParallelDescriptor::ReduceRealMax(&max_multiplic,1); - ParallelDescriptor::ReduceRealMax(&max_par_vx,1); - ParallelDescriptor::ReduceRealMax(&max_par_vy,1); - ParallelDescriptor::ReduceRealMax(&max_par_vz,1); - ParallelDescriptor::ReduceRealMax(&max_term_v,1); + }; + Real max_vals[NUM_ATTRS] = {max_par_radius, max_multiplic, max_par_mass, + max_par_vx, max_par_vy, max_par_vz, max_term_v +#ifdef ERF_USE_ML_UPHYS_DIAGNOSTICS + , max_cond_t +#endif + }; + Real sum_vals[NUM_ATTRS] = {avg_par_radius, avg_multiplic, avg_par_mass, + avg_par_vx, avg_par_vy, avg_par_vz, avg_term_v #ifdef ERF_USE_ML_UPHYS_DIAGNOSTICS - ParallelDescriptor::ReduceRealMax(&max_cond_t,1); + , avg_cond_t #endif + }; + + ParallelDescriptor::ReduceRealMin(min_vals, NUM_ATTRS); + ParallelDescriptor::ReduceRealMax(max_vals, NUM_ATTRS); + ParallelDescriptor::ReduceRealSum(sum_vals, NUM_ATTRS); + + // Unpack back to named variables + min_par_radius = min_vals[0]; min_multiplic = min_vals[1]; min_par_mass = min_vals[2]; + min_par_vx = min_vals[3]; min_par_vy = min_vals[4]; min_par_vz = min_vals[5]; + min_term_v = min_vals[6]; + + max_par_radius = max_vals[0]; max_multiplic = max_vals[1]; max_par_mass = max_vals[2]; + max_par_vx = max_vals[3]; max_par_vy = max_vals[4]; max_par_vz = max_vals[5]; + max_term_v = max_vals[6]; + + avg_par_radius = sum_vals[0]; avg_multiplic = sum_vals[1]; avg_par_mass = sum_vals[2]; + avg_par_vx = sum_vals[3]; avg_par_vy = sum_vals[4]; avg_par_vz = sum_vals[5]; + avg_term_v = sum_vals[6]; - ParallelDescriptor::ReduceRealSum(&avg_par_mass,1); - ParallelDescriptor::ReduceRealSum(&avg_par_radius,1); - ParallelDescriptor::ReduceRealSum(&avg_multiplic,1); - ParallelDescriptor::ReduceRealSum(&avg_par_vx,1); - ParallelDescriptor::ReduceRealSum(&avg_par_vy,1); - ParallelDescriptor::ReduceRealSum(&avg_par_vz,1); - ParallelDescriptor::ReduceRealSum(&avg_term_v,1); #ifdef ERF_USE_ML_UPHYS_DIAGNOSTICS - ParallelDescriptor::ReduceRealSum(&avg_cond_t,1); + min_cond_t = min_vals[7]; + max_cond_t = max_vals[7]; + avg_cond_t = sum_vals[7]; #endif if (num_total_particles > 0) { @@ -399,12 +391,12 @@ void SuperDropletPC::Diagnostics( const int a_iter, auto r_eff_max = max_par_radius; for (int ia = 0; ia < m_num_species; ia++) { const auto rho = m_species_mat[ia]->m_density; - auto r_eff_min_species = std::cbrt( min_mass_species[ia] / ((Real(4.0)/three)*PI*rho) ); + auto r_eff_min_species = std::cbrt( min_mass_species[ia] / (four_thirds_pi*rho) ); if ((r_eff_min_species < r_eff_min) && (r_eff_min_species > Real(1.0e-10))) { r_eff_min = r_eff_min_species; } } for (int ia = 0; ia < m_num_aerosols; ia++) { const auto rho = m_aerosol_mat[ia]->m_density; - auto r_eff_min_aero = std::cbrt( min_mass_aerosols[ia] / ((Real(4.0)/three)*PI*rho) ); + auto r_eff_min_aero = std::cbrt( min_mass_aerosols[ia] / (four_thirds_pi*rho) ); if ((r_eff_min_aero < r_eff_min) && (r_eff_min_aero > Real(1.0e-10))) { r_eff_min = r_eff_min_aero; } } ComputeDistributions( a_iter, r_eff_min, r_eff_max ); @@ -453,7 +445,7 @@ void SuperDropletPC::ComputeDistributions( const int a_iter, } const auto np = NumSuperDroplets(); - const ParticleReal sigma = m_sigma0 * std::exp(-Real(0.2)*std::log(static_cast(np))); + const ParticleReal sigma = m_sigma0 * std::pow(static_cast(np), Real(-0.2)); const ParticleReal lambda = one / (two*sigma*sigma); const ParticleReal gamma = one/(std::sqrt(two*PI)*sigma) * inv_bin_volume; @@ -471,7 +463,7 @@ void SuperDropletPC::ComputeDistributions( const int a_iter, auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; auto ni = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; auto mi = ptd.m_runtime_rdata[ridx_s(idx_w,na,ns)][i]; - auto ri = std::cbrt( mi / ((Real(4.0)/three)*PI*rho_w) ); + auto ri = std::cbrt( mi / (four_thirds_pi*rho_w) ); auto lnRi = std::log(ri); return gamma*ai*ni*mi*std::exp(-lambda*(lnR-lnRi)*(lnR-lnRi)); } ); @@ -483,7 +475,7 @@ void SuperDropletPC::ComputeDistributions( const int a_iter, auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; auto ni = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; auto mi = ptd.m_runtime_rdata[ridx_a(ia,na,ns)][i]; - auto ri = std::cbrt( mi / ((Real(4.0)/three)*PI*rho) ); + auto ri = std::cbrt( mi / (four_thirds_pi*rho) ); auto lnRi = std::log(ri); return gamma*ai*ni*mi*std::exp(-lambda*(lnR-lnRi)*(lnR-lnRi)); } ); @@ -497,7 +489,7 @@ void SuperDropletPC::ComputeDistributions( const int a_iter, auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; auto ni = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; auto mi = ptd.m_runtime_rdata[ridx_s(is,na,ns)][i]; - auto ri = std::cbrt( mi / ((Real(4.0)/three)*PI*rho) ); + auto ri = std::cbrt( mi / (four_thirds_pi*rho) ); auto lnRi = std::log(ri); return gamma*ai*ni*mi*std::exp(-lambda*(lnR-lnRi)*(lnR-lnRi)); } ); @@ -578,7 +570,7 @@ void SuperDropletPC::ComputeBinnedDistributions( const int a_iter) auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; auto ni = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; auto mi = ptd.m_runtime_rdata[ridx_s(idx_w,na,ns)][i]; - auto ri = std::cbrt( mi / ((Real(4.0)/three)*PI*density) ); + auto ri = std::cbrt( mi / (four_thirds_pi*density) ); auto inbin = (r_l <= ri && ri < r_r) ? one : zero; return ai*ni*mi*inbin * inv_cell_volume / dln_R; } ); @@ -588,7 +580,7 @@ void SuperDropletPC::ComputeBinnedDistributions( const int a_iter) auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; auto ni = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; auto mi = ptd.m_runtime_rdata[ridx_s(idx_w,na,ns)][i]; - auto ri = std::cbrt( mi / ((Real(4.0)/three)*PI*density) ); + auto ri = std::cbrt( mi / (four_thirds_pi*density) ); auto inbin = (r_l <= ri && ri < r_r) ? one : zero; return ai*ni*inbin * inv_cell_volume / dln_R; } ); @@ -676,7 +668,7 @@ void SuperDropletPC::ComputeBinnedDistributionsCell( const int a_iter, auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; auto ni = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; auto mi = ptd.m_runtime_rdata[ridx_s(idx_w,na,ns)][i]; - auto ri = std::cbrt( mi / ((Real(4.0)/three)*PI*density) ); + auto ri = std::cbrt( mi / (four_thirds_pi*density) ); auto inbin = (r_l <= ri && ri < r_r) ? one : zero; Gpu::Atomic::AddNoRet(&mf_arr(iv, n), (ai*ni*mi*inbin * inv_cell_volume / dln_R)); @@ -691,7 +683,7 @@ void SuperDropletPC::ComputeBinnedDistributionsCell( const int a_iter, auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; auto ni = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; auto mi = ptd.m_runtime_rdata[ridx_s(idx_w,na,ns)][i]; - auto ri = std::cbrt( mi / ((Real(4.0)/three)*PI*density) ); + auto ri = std::cbrt( mi / (four_thirds_pi*density) ); auto inbin = (r_l <= ri && ri < r_r) ? one : zero; Gpu::Atomic::AddNoRet(&mf_arr(iv, n), (ai*ni*inbin * inv_cell_volume / dln_R) ); diff --git a/Source/Particles/ERF_SuperDropletPCInitializations.cpp b/Source/Particles/ERF_SuperDropletPCInitializations.cpp index 89164304b..db308216c 100644 --- a/Source/Particles/ERF_SuperDropletPCInitializations.cpp +++ b/Source/Particles/ERF_SuperDropletPCInitializations.cpp @@ -81,12 +81,11 @@ void SuperDropletPC::readInputs (const amrex::Real a_dt) m_recyc_zmax = phi[2]; std::string coal_kernel_name = ""; - std::string term_vel_name = ""; m_coalescence_kernel = SDCoalescenceKernelType::sedimentation; coal_kernel_name = "sedimentation"; m_include_brownian_coalescence = false; - term_vel_name = "CloudRainShima"; + m_term_vel_type_w = SDTerminalVelocityType::CloudRainShima; /* read these parameters if specified */ pp.query("density_scaling", m_density_scaling); @@ -153,15 +152,11 @@ void SuperDropletPC::readInputs (const amrex::Real a_dt) amrex::Abort("Error in SuperDropletPC::readInputs() - invalid kernel choice!"); } - pp.query("terminal_velocity_model", term_vel_name); - if (term_vel_name == "AtlasUlbrich") { - m_term_vel_type = SDTerminalVelocityType::AtlasUlbrich; - } else if (term_vel_name == "RogersYau") { - m_term_vel_type = SDTerminalVelocityType::RogersYau; - } else if (term_vel_name == "CloudRainShima") { - m_term_vel_type = SDTerminalVelocityType::CloudRainShima; - } else { - amrex::Abort("Error in SuperDropletPC::readInputs() - invalid terminal velocity choice!"); + { + std::string inp_string = "terminal_velocity_model"; + if (pp.contains(inp_string.c_str())) { + pp.get(inp_string.c_str(), m_term_vel_type_w); + } } { @@ -224,13 +219,30 @@ void SuperDropletPC::define ( const std::vector& a_species_mat, setAerosolMaterial( a_aerosol_mat ); m_num_aerosols = m_aerosol_mat.size(); - AMREX_ASSERT(m_num_species > 0); - AMREX_ASSERT(m_num_species <= SupDropInit::num_species_max); - AMREX_ASSERT(m_num_aerosols <= SupDropInit::num_aerosols_max); + AMREX_ALWAYS_ASSERT(m_num_species > 0); + AMREX_ALWAYS_ASSERT(m_num_species <= SupDropInit::num_species_max); + AMREX_ALWAYS_ASSERT(m_num_aerosols <= SupDropInit::num_aerosols_max); + AMREX_ALWAYS_ASSERT(m_idx_w >= 0); add_superdroplet_attributes(); readInputs(a_dt); + // Initialize device properties for efficient GPU access + initializeDeviceProperties(); + + // Initialize staggered z-levels for non-uniform vertical grids (from terrain_z_levels) + { + ParmParse pp_erf("erf"); + int n_zlevels = pp_erf.countval("terrain_z_levels"); + if (n_zlevels > 0) { + Vector zlevels_h(n_zlevels); + pp_erf.getarr("terrain_z_levels", zlevels_h, 0, n_zlevels); + m_zlevels_d.resize(n_zlevels); + Gpu::copy(Gpu::hostToDevice, zlevels_h.begin(), zlevels_h.end(), + m_zlevels_d.begin()); + } + } + #ifdef AMREX_USE_GPU AMREX_ASSERT(!m_mass_change_logging); #endif @@ -309,14 +321,8 @@ void SuperDropletPC::InitializeParticles (const Real a_t, const MFPtr& a_ptr) Print() << "dirk2"; } Print() << " (cfl = " << m_mass_change_cfl << ")\n"; - Print() << " Terminal velocity model: "; - if (m_term_vel_type == SDTerminalVelocityType::AtlasUlbrich) { - Print() << "AtlasUlbrich" << "\n"; - } else if (m_term_vel_type == SDTerminalVelocityType::RogersYau) { - Print() << "RogersYau" << "\n"; - } else if (m_term_vel_type == SDTerminalVelocityType::CloudRainShima) { - Print() << "CloudRainShima" << "\n"; - } + Print() << " Terminal velocity model: " + << getEnumNameString(m_term_vel_type_w) << "\n"; for (int i = 0; i < m_num_initializations; i++) { Print() << " SuperDropletPC(" << m_name << ") Initialization"; @@ -378,37 +384,13 @@ void SuperDropletPC::SetAttributes (MultiFab& a_rhoc /*!< mass density of conden const Real rho_w = m_species_mat[m_idx_w]->m_density; const int idx_w = m_idx_w; -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (ParIterType pti(*this, m_lev); pti.isValid(); ++pti) { - auto& particle_tile = ParticlesAt(m_lev, pti); - auto& soa = particle_tile.GetStructOfArrays(); - auto& aos = particle_tile.GetArrayOfStructs(); - auto *p_pbox = aos().data(); - const int n = aos.numParticles(); - - /* SoA attributes */ - auto* mass_ptr = soa.GetRealData(SuperDropletsRealIdxSoA::mass).data(); - - /* Runtime-added SoA attributes */ - int rt_off_r = SuperDropletsRealIdxSoA::ncomps; - auto* radius_ptr = soa.GetRealData(rt_off_r+SuperDropletsRealIdxSoA_RT::radius).data(); - auto* mult_ptr = soa.GetRealData(rt_off_r+SuperDropletsRealIdxSoA_RT::multiplicity).data(); - - SDSpeciesMassArr sp_mass_ptrs; - for (int i = 0; i < num_sp; i++) { - sp_mass_ptrs[i] = soa.GetRealData(idx_s(i,num_ae,num_sp)).data(); - } - - SDAerosolMassArr ae_mass_ptrs; - for (int i = 0; i < num_ae; i++) { - ae_mass_ptrs[i] = soa.GetRealData(idx_a(i,num_ae,num_sp)).data(); - } - - auto condensate_mass_density = a_rhoc[pti.index()].array(); + forEachParticleTile(m_lev, + [&](ParIterType& /*pti*/, int grid, ParticleType* p_pbox, + const SDProcess::ParticlePointers& ptrs, int np) + { + auto condensate_mass_density = a_rhoc[grid].array(); - ParallelForRNG(n, [=] AMREX_GPU_DEVICE (int i, const RandomEngine& rnd_engine) + ParallelForRNG(np, [=] AMREX_GPU_DEVICE (int i, const RandomEngine& rnd_engine) { ParticleType& p = p_pbox[i]; if (p.id() <= 0) { return; } @@ -417,32 +399,29 @@ void SuperDropletPC::SetAttributes (MultiFab& a_rhoc /*!< mass density of conden const Real mass_condensate_cell = condensate_mass_density(iv[0],iv[1],iv[2],0) * cell_volume; const Real mass_condensate_sd = mass_condensate_cell / num_sd_per_cell; - Real mult_rnd = -mult_ptr[i]/3 + 2*mult_ptr[i]/3*Random(rnd_engine); - mult_ptr[i] += mult_rnd; + Real mult_rnd = -ptrs.mult_ptr[i]/3 + 2*ptrs.mult_ptr[i]/3*Random(rnd_engine); + ptrs.mult_ptr[i] += mult_rnd; ParticleReal species_mass_total = zero; for (int ctr = 0; ctr < num_sp; ctr++) { if (ctr != idx_w) { - species_mass_total += sp_mass_ptrs[ctr][n]; + species_mass_total += ptrs.sp_mass_ptrs[ctr][np]; } } ParticleReal aerosol_mass_total = zero; for (int ctr = 0; ctr < num_ae; ctr++) { - aerosol_mass_total += ae_mass_ptrs[ctr][n]; + aerosol_mass_total += ptrs.ae_mass_ptrs[ctr][np]; } - const Real mass_particle = mass_condensate_sd / mult_ptr[i] + aerosol_mass_total + species_mass_total; - mass_ptr[i] = mass_particle; + const Real mass_particle = mass_condensate_sd / ptrs.mult_ptr[i] + aerosol_mass_total + species_mass_total; + ptrs.mass_ptr[i] = mass_particle; - Real radius_cubed = mass_particle / ((Real(4.0)/three)*PI*rho_w); + Real radius_cubed = mass_particle / (four_thirds_pi*rho_w); Real radius = (radius_cubed == zero ? zero : std::cbrt(radius_cubed)); - radius_ptr[i] = radius; + ptrs.radius_ptr[i] = radius; }); - - } - - return; + }); // end forEachParticleTile } /*! Scale the multiplicities with density of air */ @@ -455,35 +434,22 @@ void SuperDropletPC::DensityScaling (const MultiFab& a_rho /*!< density of air * const auto dxi = Geom(m_lev).InvCellSizeArray(); const auto domain = Geom(m_lev).Domain(); -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (ParIterType pti(*this, m_lev); pti.isValid(); ++pti) { - auto& particle_tile = ParticlesAt(m_lev, pti); - auto& soa = particle_tile.GetStructOfArrays(); - auto& aos = particle_tile.GetArrayOfStructs(); - auto *p_pbox = aos().data(); - const int n = aos.numParticles(); - - /* Runtime-added SoA attributes */ - int rt_off_r = SuperDropletsRealIdxSoA::ncomps; - auto* mult_ptr = soa.GetRealData(rt_off_r+SuperDropletsRealIdxSoA_RT::multiplicity).data(); - - auto density = a_rho[pti.index()].const_array(); + forEachParticleTile(m_lev, + [&](ParIterType& /*pti*/, int grid, ParticleType* p_pbox, + const SDProcess::ParticlePointers& ptrs, int np) + { + auto density = a_rho[grid].const_array(); - ParallelFor(n, [=] AMREX_GPU_DEVICE (int i) + ParallelFor(np, [=] AMREX_GPU_DEVICE (int i) { ParticleType& p = p_pbox[i]; if (p.id() <= 0) { return; } auto iv = getParticleCell(p, plo, dxi, domain); auto rho_air = density(iv[0],iv[1],iv[2],0); - mult_ptr[i] *= rho_air; + ptrs.mult_ptr[i] *= rho_air; }); - - } - - return; + }); // end forEachParticleTile } #endif diff --git a/Source/Particles/ERF_SuperDropletPCMassChange.H b/Source/Particles/ERF_SuperDropletPCMassChange.H index 8ac7a53a2..640c5aa93 100644 --- a/Source/Particles/ERF_SuperDropletPCMassChange.H +++ b/Source/Particles/ERF_SuperDropletPCMassChange.H @@ -1,11 +1,14 @@ -#ifndef SUPERDROPLET_PC_COALESCENCE_H_ -#define SUPERDROPLET_PC_COALESCENCE_H_ +#ifndef SUPERDROPLET_PC_MASSCHANGE_H_ +#define SUPERDROPLET_PC_MASSCHANGE_H_ #include +#include "ERF_Constants.H" #ifdef ERF_USE_PARTICLES -namespace SDMassChangeUtils +/*! \brief Namespace with classes and functions for condensation/evaporation + * (liquid <--> vapour) */ +namespace SDMassChangeUtils_LV { #ifdef ERF_USE_ML_UPHYS_DIAGNOSTICS /*! \brief Phase change equation (in terms of R) */ @@ -32,13 +35,13 @@ namespace SDMassChangeUtils RT Kn = lambda_v/a_R; RT dcf = (1+Kn) / (1+2*Kn*(1+Kn)); - RT F_k = ( L/(Rv*a_T) - one) * ((L*rho_l) / (K*a_T)); + RT F_k = ( L/(Rv*a_T) - RT(1.0)) * ((L*rho_l) / (K*a_T)); RT F_d = (rho_l*Rv*a_T) / (dcf*a_D*a_e_s); - RT R_inv = one/a_R; + RT R_inv = RT(1.0)/a_R; RT R_inv_cubed = R_inv*R_inv*R_inv; - RT alpha = (a_S-one) / (F_k + F_d); + RT alpha = (a_S-RT(1.0)) / (F_k + F_d); RT retval = alpha*R_inv; RT beta = -(a_a/a_T) / (F_k + F_d); @@ -77,19 +80,19 @@ namespace SDMassChangeUtils RT Kn = lambda_v/std::sqrt(a_R_sq); RT dcf = (1+Kn) / (1+2*Kn*(1+Kn)); - RT F_k = ( L/(Rv*a_T) - one) * ((L*rho_l) / (K*a_T)); + RT F_k = ( L/(Rv*a_T) - RT(1.0)) * ((L*rho_l) / (K*a_T)); RT F_d = (rho_l*Rv*a_T) / (dcf*a_D*a_e_s); - RT R_inv = std::exp(-myhalf*std::log(a_R_sq)); + RT R_inv = RT(1.0)/std::sqrt(a_R_sq); RT R_inv_cubed = R_inv*R_inv*R_inv; - RT alpha = two * (a_S-one) / (F_k + F_d); + RT alpha = RT(2.0) * (a_S-RT(1.0)) / (F_k + F_d); RT retval = alpha; - RT beta = -two * (a_a/a_T) / (F_k + F_d); + RT beta = -RT(2.0) * (a_a/a_T) / (F_k + F_d); retval += beta*R_inv; - RT gamma = two * a_b * a_N_s / (F_k + F_d); + RT gamma = RT(2.0) * a_b * a_N_s / (F_k + F_d); retval += gamma*R_inv_cubed; return retval; @@ -109,19 +112,19 @@ namespace SDMassChangeUtils RT Kn = lambda_v/std::sqrt(a_R_sq); RT dcf = (1+Kn) / (1+2*Kn*(1+Kn)); - RT F_k = ( L/(Rv*a_T) - one) * ((L*rho_l) / (K*a_T)); + RT F_k = ( L/(Rv*a_T) - RT(1.0)) * ((L*rho_l) / (K*a_T)); RT F_d = (rho_l*Rv*a_T) / (dcf*a_D*a_e_s); - RT R_inv = one/std::sqrt(a_R_sq); + RT R_inv = RT(1.0)/std::sqrt(a_R_sq); RT R_inv_3 = R_inv*R_inv*R_inv; RT R_inv_5 = R_inv_3*R_inv*R_inv; - RT retval = zero; + RT retval = RT(0.0); - RT beta = -two * (a_a/a_T) / (F_k + F_d); + RT beta = -RT(2.0) * (a_a/a_T) / (F_k + F_d); retval -= myhalf * beta*R_inv_3; - RT gamma = two * a_b * a_N_s / (F_k + F_d); + RT gamma = RT(2.0) * a_b * a_N_s / (F_k + F_d); retval -= myhalf * three*gamma*R_inv_5; return retval; @@ -167,7 +170,7 @@ namespace SDMassChangeUtils ) const { a_converged = false; - RT res_norm0 = zero; + RT res_norm0 = RT(0.0); for (int k = 0; k < m_maxits; k++) { RT residual = a_mu * a_u @@ -179,7 +182,7 @@ namespace SDMassChangeUtils if (a_res_norm_a > 0) { res_norm0 = a_res_norm_a; } else { - res_norm0 = one; + res_norm0 = RT(1.0); } } a_res_norm_r = a_res_norm_a / res_norm0; @@ -198,7 +201,7 @@ namespace SDMassChangeUtils } RT slope = a_mu - m_ne.rhs_jac( a_u, a_T, a_e_s, a_D, a_a, a_b, a_N_s ); - RT du = zero; + RT du = RT(0.0); du = - residual / slope; RT du_norm = std::sqrt(du*du); @@ -209,7 +212,7 @@ namespace SDMassChangeUtils } a_u += du; - if (a_u <= zero) { + if (a_u <= RT(0.0)) { a_converged = false; break; } @@ -244,13 +247,76 @@ namespace SDMassChangeUtils bool m_adapt_dt; /*!< use error-based adaptive dt? */ bool m_verbose; /*!< verbosity */ + /*! \brief Compute timestep from stiffness estimate */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + RT computeTimestep (const RT& a_u) const + { + RT tau = m_ode.rhs_jac(a_u, m_T, m_e_s, m_D, m_a, m_b, m_N_s); + return m_cfl / std::sqrt(tau*tau); + } + + /*! \brief Compute stiffness estimate (tau) */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + RT computeTau (const RT& a_u) const + { + return m_ode.rhs_jac(a_u, m_T, m_e_s, m_D, m_a, m_b, m_N_s); + } + + /*! \brief Limit timestep to not exceed final time */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + RT limitTimestep (RT dt, RT cur_time) const + { + if ((cur_time + dt) > m_t_final) { + dt = m_t_final - cur_time; + } + return dt; + } + + /*! \brief Check if timestep is too small to continue */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + bool isTimestepTooSmall (RT dt, RT tau) const + { + return (dt < (RT(1.0e-12) * m_cfl / std::sqrt(tau*tau))) + && (dt < (RT(1.0e-12) * m_t_final)); + } + + /*! \brief Evaluate ODE right-hand side */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + RT evalRHS (const RT& a_u) const + { + return m_ode.rhs_func(a_u, m_S, m_T, m_e_s, m_D, m_a, m_b, m_N_s); + } + + /*! \brief Print verbose step information */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void printStepInfo (RT cur_time, RT dt, RT tau, RT radius, RT snorm) const + { + if (m_verbose) { + AMREX_DEVICE_PRINTF("Time %1.2e, dt = %1.2e, cfl = %1.1e, radius = %1.4e, snorm = %1.1e\n", + cur_time, dt, dt * std::sqrt(tau*tau), radius, snorm); + } + } + + /*! \brief Print verbose step info with Newton solver details */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void printStepInfoNewton (RT cur_time, RT dt, RT tau, RT radius, RT snorm, + RT res_norm_a, RT res_norm_r, bool converged) const + { + if (m_verbose) { + AMREX_DEVICE_PRINTF("Time %1.2e, dt = %1.2e, cfl = %1.1e, radius = %1.4e, snorm = %1.1e\n", + cur_time, dt, dt * std::sqrt(tau*tau), radius, snorm); + AMREX_DEVICE_PRINTF(" norms = %1.3e (abs), %1.3e (rel), converged = %s\n", + res_norm_a, res_norm_r, (converged ? "yes" : "no")); + } + } + /*! \brief 3rd-order, 4-stage Bogacki-Shampine explicit RK method * with 2nd order embedded method */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void rk3bs ( RT& a_u, /*!< solution */ bool& a_success /*!< success/failure flag */ ) const { - RT cur_time = zero; + RT cur_time = RT(0.0); a_success = true; RT tau = m_ode.rhs_jac(a_u, m_T, m_e_s, m_D, m_a, m_b, m_N_s); diff --git a/Source/Particles/ERF_SuperDropletPCMassChange.cpp b/Source/Particles/ERF_SuperDropletPCMassChange.cpp index 30f773495..fda30ebba 100644 --- a/Source/Particles/ERF_SuperDropletPCMassChange.cpp +++ b/Source/Particles/ERF_SuperDropletPCMassChange.cpp @@ -4,12 +4,15 @@ #include #include "ERF_Constants.H" #include "ERF_SuperDropletPCMassChange.H" +#include "ERF_InterpolationUtils.H" #ifdef ERF_USE_PARTICLES using namespace amrex; using namespace SDPCDefn; -using namespace SDMassChangeUtils; + +/*! \brief Field indices for liquid-vapour interpolation */ +AMREX_ENUM(InterpFieldsLV, e_sat, sat_ratio, temperature, pressure, NUM_FIELDS); /*! Compute mass change of particles due to evaporation and condensation */ void SuperDropletPC::MassChange ( int a_lev, @@ -22,24 +25,19 @@ void SuperDropletPC::MassChange ( int a_ const Vector>& a_z_phys_nd, const bool a_is_water ) { + using namespace SDMassChangeUtils_LV; + BL_PROFILE("SuperDropletPC::MassChange()"); AMREX_ASSERT( a_lev == m_lev ); - const Geometry& geom = m_gdb->Geom(a_lev); - const auto plo = geom.ProbLoArray(); - const auto dxi = geom.InvCellSizeArray(); + const auto ctx = buildProcessContext(a_lev); - const auto is_periodic = geom.isPeriodic(); - auto is_periodic_z = is_periodic[2]; + const Geometry& geom = m_gdb->Geom(a_lev); + const auto is_periodic_z = geom.isPeriodic(2); const std::unique_ptr& z_height = a_z_phys_nd[a_lev]; - const int num_sp = m_num_species; - const int num_ae = m_num_aerosols; - - const Real rho_w = m_species_mat[m_idx_w]->m_density; - int idx_w = m_idx_w; - + // Find vapour species index int idx_vap = -1; Real mat_density = -1; for (int i = 0; i < m_num_species; i++) { @@ -55,128 +53,50 @@ void SuperDropletPC::MassChange ( int a_ const bool log_unconverged = m_mass_change_logging; [[maybe_unused]] FILE* file_handle = m_mass_change_log; + auto cfl = m_mass_change_cfl; + auto ti_choice = m_mass_change_ti; + + // Solver setup (shared across tiles) + dRsqdt drsqdt{ vapour_mat.m_lat_vap, + therco, + vapour_mat.m_Rv, + mat_density }; + + NewtonSolver< dRsqdt, ParticleReal > newton_solver{ drsqdt, + m_newton_rtol, + m_newton_atol, + m_newton_stol, + m_newton_maxits }; -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) +#ifdef ERF_USE_ML_UPHYS_DIAGNOSTICS + dRdt drdt{ vapour_mat.m_lat_vap, + therco, + vapour_mat.m_Rv, + mat_density }; + constexpr int rtoff_r = SuperDropletsRealIdxSoA::ncomps; #endif - for (ParIterType pti(*this, a_lev); pti.isValid(); ++pti) { - int grid = pti.index(); - auto& ptile = ParticlesAt(a_lev, pti); - auto& aos = ptile.GetArrayOfStructs(); - auto& soa = ptile.GetStructOfArrays(); - const int num_particles = aos.numParticles(); - auto* p_pbox = aos().data(); + forEachParticleTile(a_lev, ctx, + [&](ParIterType& pti, int grid, ParticleType* p_pbox, + const SDProcess::ParticlePointers& ptrs, + const SDProcess::ProcessContext& ctx) + { + amrex::ignore_unused(pti); auto zheight = (*z_height)[grid].array(); - auto* mass_ptr = soa.GetRealData(SuperDropletsRealIdxSoA::mass).data(); - - int rtoff_i = SuperDropletsIntIdxSoA::ncomps; - auto* active_ptr = soa.GetIntData(rtoff_i+SuperDropletsIntIdxSoA_RT::active).data(); - int rtoff_r = SuperDropletsRealIdxSoA::ncomps; - auto* radius_ptr = soa.GetRealData(rtoff_r+SuperDropletsRealIdxSoA_RT::radius).data(); #ifdef ERF_USE_ML_UPHYS_DIAGNOSTICS + auto& soa = ParticlesAt(a_lev, pti).GetStructOfArrays(); auto* condt_ptr = soa.GetRealData(rtoff_r+SuperDropletsRealIdxSoA_RT::cond_tendency).data(); #endif - SDSpeciesMassArr sp_mass_ptrs; - Gpu::DeviceVector sp_ionization(num_sp); - Gpu::DeviceVector sp_mol_weight(num_sp); - Gpu::DeviceVector sp_density(num_sp); - Gpu::DeviceVector sp_solubility(num_sp); - { - Vector sp_ionization_h(num_sp); - Vector sp_mol_weight_h(num_sp); - Vector sp_density_h(num_sp); - Vector sp_solubility_h(num_sp); - for (int i = 0; i < num_sp; i++) { - sp_mass_ptrs[i] = soa.GetRealData(idx_s(i,num_ae,num_sp)).data(); - sp_ionization_h[i] = m_species_mat[i]->m_ionization; - sp_mol_weight_h[i] = m_species_mat[i]->m_mol_weight; - sp_density_h[i] = m_species_mat[i]->m_density; - sp_solubility_h[i] = static_cast(m_species_mat[i]->m_is_soluble); - } - Gpu::copy( Gpu::hostToDevice, - sp_ionization_h.begin(), - sp_ionization_h.end(), - sp_ionization.begin() ); - Gpu::copy( Gpu::hostToDevice, - sp_mol_weight_h.begin(), - sp_mol_weight_h.end(), - sp_mol_weight.begin() ); - Gpu::copy( Gpu::hostToDevice, - sp_density_h.begin(), - sp_density_h.end(), - sp_density.begin() ); - Gpu::copy( Gpu::hostToDevice, - sp_solubility_h.begin(), - sp_solubility_h.end(), - sp_solubility.begin() ); - } - - SDAerosolMassArr ae_mass_ptrs; - Gpu::DeviceVector ae_ionization(num_ae); - Gpu::DeviceVector ae_mol_weight(num_ae); - Gpu::DeviceVector ae_density(num_ae); - Gpu::DeviceVector ae_solubility(num_ae); - { - Vector ae_ionization_h(num_ae); - Vector ae_mol_weight_h(num_ae); - Vector ae_density_h(num_ae); - Vector ae_solubility_h(num_ae); - for (int i = 0; i < num_ae; i++) { - ae_mass_ptrs[i] = soa.GetRealData(idx_a(i,num_ae,num_sp)).data(); - ae_ionization_h[i] = m_aerosol_mat[i]->m_ionization; - ae_mol_weight_h[i] = m_aerosol_mat[i]->m_mol_weight; - ae_density_h[i] = m_aerosol_mat[i]->m_density; - ae_solubility_h[i] = static_cast(m_aerosol_mat[i]->m_is_soluble); - } - Gpu::copy( Gpu::hostToDevice, - ae_ionization_h.begin(), - ae_ionization_h.end(), - ae_ionization.begin() ); - Gpu::copy( Gpu::hostToDevice, - ae_mol_weight_h.begin(), - ae_mol_weight_h.end(), - ae_mol_weight.begin() ); - Gpu::copy( Gpu::hostToDevice, - ae_density_h.begin(), - ae_density_h.end(), - ae_density.begin() ); - Gpu::copy( Gpu::hostToDevice, - ae_solubility_h.begin(), - ae_solubility_h.end(), - ae_solubility.begin() ); - } - const auto& sat_pressure_arr = a_sat_pressure[grid].array(); const auto& sat_ratio_arr = a_sat_ratio[grid].array(); const auto& temperature_arr = a_temperature[grid].array(); const auto& pressure_arr = a_pressure[grid].array(); -#ifdef ERF_USE_ML_UPHYS_DIAGNOSTICS - dRdt drdt{ vapour_mat.m_lat_vap, - therco, /* ERF_Constants.H */ - vapour_mat.m_Rv, - mat_density }; -#endif - - dRsqdt drsqdt{ vapour_mat.m_lat_vap, - therco, /* ERF_Constants.H */ - vapour_mat.m_Rv, - mat_density }; - - NewtonSolver< dRsqdt, ParticleReal > newton_solver{ drsqdt, - m_newton_rtol, - m_newton_atol, - m_newton_stol, - m_newton_maxits }; - Gpu::Buffer unconverged_particles({0}); auto* unconverged_particles_ptr = unconverged_particles.data(); - auto cfl = m_mass_change_cfl; - auto ti_choice = m_mass_change_ti; AMREX_ASSERT_WITH_MESSAGE( ti_choice == SDMassChangeTIMethod::RK4 || ti_choice == SDMassChangeTIMethod::RK3BS || ti_choice == SDMassChangeTIMethod::BE || @@ -184,43 +104,36 @@ void SuperDropletPC::MassChange ( int a_ ti_choice == SDMassChangeTIMethod::DIRK2, "ERROR: invalid time integrator choice!" ); - auto sp_i_arr = sp_ionization.data(); - auto sp_mw_arr = sp_mol_weight.data(); - auto sp_rho_arr = sp_density.data(); - auto sp_sol_arr = sp_solubility.data(); - auto ae_i_arr = ae_ionization.data(); - auto ae_mw_arr = ae_mol_weight.data(); - auto ae_rho_arr = ae_density.data(); - auto ae_sol_arr = ae_solubility.data(); - - ParallelFor(num_particles, [=] AMREX_GPU_DEVICE (int i) + ParallelFor(ptrs.num_particles, [=] AMREX_GPU_DEVICE (int i) { ParticleType& p = p_pbox[i]; if (p.id() <= 0) { return; } - if (active_ptr[i] == 0) { return; } - - ParticleReal sat_ratio, e_sat, temperature, pressure; - if (is_periodic_z) { - cic_interpolate( p, plo, dxi, sat_pressure_arr, &e_sat, 1 ); - cic_interpolate( p, plo, dxi, sat_ratio_arr, &sat_ratio, 1 ); - cic_interpolate( p, plo, dxi, temperature_arr, &temperature, 1 ); - cic_interpolate( p, plo, dxi, pressure_arr, &pressure, 1 ); - } else { - cic_interpolate_mapped_z( p, plo, dxi, sat_pressure_arr, zheight, &e_sat, 1 ); - cic_interpolate_mapped_z( p, plo, dxi, sat_ratio_arr, zheight, &sat_ratio, 1 ); - cic_interpolate_mapped_z( p, plo, dxi, temperature_arr, zheight, &temperature, 1 ); - cic_interpolate_mapped_z( p, plo, dxi, pressure_arr, zheight, &pressure, 1 ); - } + if (ptrs.active_ptr[i] == 0) { return; } + + // Interpolate saturation pressure, saturation ratio, temperature, pressure + constexpr int nf = static_cast(InterpFieldsLV::NUM_FIELDS); + ParticleReal fv[nf]; + const Array4 fa[nf] = { + sat_pressure_arr, sat_ratio_arr, temperature_arr, pressure_arr + }; + ERF::Interpolation::interpolateFields( + p, ctx.plo, ctx.dxi, fa, fv, nf, + is_periodic_z ? 1 : 0, is_periodic_z ? nullptr : &zheight + ); + const auto e_sat = fv[static_cast(InterpFieldsLV::e_sat)]; + const auto sat_ratio = fv[static_cast(InterpFieldsLV::sat_ratio)]; + const auto temperature = fv[static_cast(InterpFieldsLV::temperature)]; + const auto pressure = fv[static_cast(InterpFieldsLV::pressure)]; ParticleReal solute_moles = zero; if (a_is_water) { - for (int j = 0; j < num_sp; j++) { + for (int j = 0; j < ctx.num_species; j++) { if (j != idx_vap) { - solute_moles += (sp_mass_ptrs[j][i]*sp_i_arr[j]/sp_mw_arr[j]); + solute_moles += (ptrs.sp_mass_ptrs[j][i]*ptrs.sp_ion_arr[j]/ptrs.sp_mw_arr[j]); } } - for (int j = 0; j < num_ae; j++) { - solute_moles += (ae_mass_ptrs[j][i]*ae_i_arr[j]/ae_mw_arr[j]); + for (int j = 0; j < ctx.num_aerosols; j++) { + solute_moles += (ptrs.ae_mass_ptrs[j][i]*ptrs.ae_ion_arr[j]/ptrs.ae_mw_arr[j]); } } @@ -230,7 +143,7 @@ void SuperDropletPC::MassChange ( int a_ #ifdef ERF_USE_ML_UPHYS_DIAGNOSTICS if (a_is_water) { - condt_ptr[i] = drdt( radius_ptr[i], + condt_ptr[i] = drdt( ptrs.radius_ptr[i], sat_ratio, temperature, e_sat, @@ -249,12 +162,12 @@ void SuperDropletPC::MassChange ( int a_ solute_moles, cfl, 1e-40, 1e-3, 1e-6, false, false }; - auto r_init = SD_effective_radius( i, idx_w, - rho_w, - num_sp, num_ae, - sp_sol_arr, ae_sol_arr, - sp_mass_ptrs, ae_mass_ptrs, - sp_rho_arr, ae_rho_arr ); + auto r_init = SD_effective_radius( i, ctx.idx_water, + ctx.rho_water, + ctx.num_species, ctx.num_aerosols, + ptrs.sp_sol_arr, ptrs.ae_sol_arr, + ptrs.sp_mass_ptrs, ptrs.ae_mass_ptrs, + ptrs.sp_rho_arr, ptrs.ae_rho_arr ); auto r_sq = r_init*r_init; bool success = false; @@ -275,31 +188,29 @@ void SuperDropletPC::MassChange ( int a_ #ifndef AMREX_USE_GPU fprintf(file_handle, "r=%1.16e, S=%1.16e, T=%1.16e, e=%1.16e, sol_mass=%1.16e\n", - radius_ptr[i], sat_ratio, temperature, e_sat, solute_moles ); + ptrs.radius_ptr[i], sat_ratio, temperature, e_sat, solute_moles ); #endif } Gpu::Atomic::Add(unconverged_particles_ptr, Long(1)); } else { // update particle attributes auto r_new = std::sqrt(r_sq); - auto d_mass = (Real(4.0)/three)*PI*mat_density * (r_new*r_new*r_new - r_init*r_init*r_init); - sp_mass_ptrs[idx_vap][i] += d_mass; + auto d_mass = four_thirds_pi*mat_density * (r_new*r_new*r_new - r_init*r_init*r_init); + ptrs.sp_mass_ptrs[idx_vap][i] += d_mass; // don't let it go negative - sp_mass_ptrs[idx_vap][i] = std::max(sp_mass_ptrs[idx_vap][i],amrex::Real(0)); + ptrs.sp_mass_ptrs[idx_vap][i] = std::max(ptrs.sp_mass_ptrs[idx_vap][i],amrex::Real(0)); - radius_ptr[i] = SD_effective_radius( i, idx_w, - rho_w, - num_sp, num_ae, - sp_sol_arr, ae_sol_arr, - sp_mass_ptrs, ae_mass_ptrs, - sp_rho_arr, ae_rho_arr ); - mass_ptr[i] = SD_total_mass( i, num_sp, num_ae, sp_mass_ptrs, ae_mass_ptrs); + // Update particle attributes (radius and mass) + SuperDropletPC::updateParticleAttributes( + i, ptrs.radius_ptr, ptrs.mass_ptr, ctx.idx_water, ctx.rho_water, + ctx.num_species, ctx.num_aerosols, ptrs.sp_sol_arr, ptrs.ae_sol_arr, + ptrs.sp_mass_ptrs, ptrs.ae_mass_ptrs, ptrs.sp_rho_arr, ptrs.ae_rho_arr); } }); Gpu::synchronize(); m_num_unconverged_particles += *(unconverged_particles.copyToHost()); - } + }); // end forEachParticleTile } diff --git a/Source/Particles/ERF_SuperDropletPCProcess.H b/Source/Particles/ERF_SuperDropletPCProcess.H new file mode 100644 index 000000000..1143b4dba --- /dev/null +++ b/Source/Particles/ERF_SuperDropletPCProcess.H @@ -0,0 +1,82 @@ +/*! \file ERF_SuperDropletPCProcess.H + * \brief Common data structures for SuperDroplet physical processes + * + * This file provides shared infrastructure for physical process implementations: + * 1. ParticlePointers - Bundle of particle attribute pointers for GPU kernels + * 2. ProcessContext - Non-pointer data (geometry, species info) shared across processes + */ + +#ifndef ERF_SUPERDROPLET_PC_PROCESS_H_ +#define ERF_SUPERDROPLET_PC_PROCESS_H_ + +#ifdef ERF_USE_PARTICLES + +#include "ERF_SuperDropletPCDefinitions.H" +#include +#include + +namespace SDProcess { + +using namespace amrex; +using namespace SDPCDefn; + +/*! \brief Bundle of all particle attribute pointers + * + * This struct holds pointers to particle data for use in GPU kernels. + * All pointers are populated by setupParticlePointers(); processes use + * only the pointers they need. + */ +struct ParticlePointers { + // Number of particles in this tile + int num_particles{0}; + + // Core SOA pointers + ParticleReal* mass_ptr{nullptr}; + ParticleReal* radius_ptr{nullptr}; + int* active_ptr{nullptr}; + + // Velocity components + GpuArray v_ptr{nullptr, nullptr, nullptr}; + ParticleReal* vterm_ptr{nullptr}; + ParticleReal* mult_ptr{nullptr}; + + // Species and aerosol mass arrays + SDSpeciesMassArr sp_mass_ptrs{}; + SDAerosolMassArr ae_mass_ptrs{}; + + // Material properties (read-only) + const ParticleReal* sp_rho_arr{nullptr}; + const int* sp_sol_arr{nullptr}; + const ParticleReal* sp_ion_arr{nullptr}; + const ParticleReal* sp_mw_arr{nullptr}; + + const ParticleReal* ae_rho_arr{nullptr}; + const int* ae_sol_arr{nullptr}; + const ParticleReal* ae_ion_arr{nullptr}; + const ParticleReal* ae_mw_arr{nullptr}; +}; + +/*! \brief Context data that doesn't change per-tile */ +struct ProcessContext { + // Geometry + GpuArray plo; + GpuArray phi; + GpuArray dxi; + GpuArray dx; + GpuArray is_periodic; + Box domain; + Real cell_volume{0.0}; + + // Species/aerosol counts and indices + int num_species{0}; + int num_aerosols{0}; + int idx_water{-1}; + + // Material densities (commonly used) + ParticleReal rho_water{0.0}; +}; + +} // namespace SDProcess + +#endif // ERF_USE_PARTICLES +#endif // ERF_SUPERDROPLET_PC_PROCESS_H_ diff --git a/Source/Particles/ERF_SuperDropletPCRecycle.cpp b/Source/Particles/ERF_SuperDropletPCRecycle.cpp index 5f508e70c..830bfdf02 100644 --- a/Source/Particles/ERF_SuperDropletPCRecycle.cpp +++ b/Source/Particles/ERF_SuperDropletPCRecycle.cpp @@ -1,12 +1,14 @@ #include #include #include "ERF_SuperDropletPC.H" +#include "ERF_InterpolationUtils.H" #ifdef ERF_USE_PARTICLES using namespace amrex; using namespace SDPCDefn; + /*! Recycle deactivated particles: particles that have zero multiplicity are * recycled by resetting them to dry aerosol particles and placing them randomly * in the domain. */ @@ -20,7 +22,7 @@ void SuperDropletPC::Recycle ( const int a_lev, const auto num_sd_deactivated = NumSDDeactivated(); const auto num_sd = NumSuperDroplets(); - const auto deac_frac = static_cast(num_sd_deactivated) / static_cast(num_sd); + const auto deac_frac = (num_sd > 0 ? static_cast(num_sd_deactivated)/static_cast(num_sd) : 0.0); int flag = 0; if (deac_frac > m_deac_threshold) { flag = 1; } @@ -94,14 +96,10 @@ void SuperDropletPC::Recycle ( const int a_lev, const auto init_r = *(m_initializations[init_idx]); const auto sampled_multiplicity = init_r.sampledMultiplicity(); + const auto ctx = buildProcessContext(a_lev); const auto dx_h = Geom(m_lev).CellSize(); const Real cell_volume = dx_h[0]*dx_h[1]*dx_h[2]; - const int idx_w = m_idx_w; - const Real rho_w = m_species_mat[m_idx_w]->m_density; - const int num_sp = m_num_species; - const int num_ae = m_num_aerosols; - // number of super-droplets per cell int num_sd_per_cell = m_num_sd_per_cell; // number of physical particles per cell @@ -113,89 +111,28 @@ void SuperDropletPC::Recycle ( const int a_lev, auto mult_avg = num_par_per_cell / num_sd_per_cell; Long np_recycle_total = 0; -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (ParIterType pti(*this, a_lev); pti.isValid(); ++pti) { - auto& ptile = ParticlesAt(a_lev, pti); - auto& aos = ptile.GetArrayOfStructs(); - auto& soa = ptile.GetStructOfArrays(); - const int np = aos.numParticles(); - auto *p_pbox = aos().data(); - - Array v_ptr; - v_ptr[0] = soa.GetRealData(SuperDropletsRealIdxSoA::vx).data(); - v_ptr[1] = soa.GetRealData(SuperDropletsRealIdxSoA::vy).data(); - v_ptr[2] = soa.GetRealData(SuperDropletsRealIdxSoA::vz).data(); - auto* mass_ptr = soa.GetRealData(SuperDropletsRealIdxSoA::mass).data(); - - int rtoff_i = SuperDropletsIntIdxSoA::ncomps; - auto* active_ptr = soa.GetIntData(rtoff_i+SuperDropletsIntIdxSoA_RT::active).data(); - int rtoff_r = SuperDropletsRealIdxSoA::ncomps; - auto* radius_ptr = soa.GetRealData(rtoff_r+SuperDropletsRealIdxSoA_RT::radius).data(); - auto* vterm_ptr = soa.GetRealData(rtoff_r+SuperDropletsRealIdxSoA_RT::term_vel).data(); - auto* mult_ptr = soa.GetRealData(rtoff_r+SuperDropletsRealIdxSoA_RT::multiplicity).data(); - - SDSpeciesMassArr sp_mass_ptrs; - for (int i = 0; i < num_sp; i++) { - sp_mass_ptrs[i] = soa.GetRealData(idx_s(i,num_ae,num_sp)).data(); - } - SDAerosolMassArr ae_mass_ptrs; - for (int i = 0; i < num_ae; i++) { - ae_mass_ptrs[i] = soa.GetRealData(idx_a(i,num_ae,num_sp)).data(); - } - - Gpu::DeviceVector sp_density(num_sp); - Gpu::DeviceVector sp_solubility(num_sp); - { - Vector sp_density_h(num_sp); - Vector sp_solubility_h(num_sp); - for (int i = 0; i < num_sp; i++) { - sp_density_h[i] = m_species_mat[i]->m_density; - sp_solubility_h[i] = static_cast(m_species_mat[i]->m_is_soluble); - } - Gpu::copy( Gpu::hostToDevice, - sp_density_h.begin(), - sp_density_h.end(), - sp_density.begin() ); - Gpu::copy( Gpu::hostToDevice, - sp_solubility_h.begin(), - sp_solubility_h.end(), - sp_solubility.begin() ); - } - - Gpu::DeviceVector ae_density(num_ae); - Gpu::DeviceVector ae_solubility(num_ae); - { - Vector ae_density_h(num_ae); - Vector ae_solubility_h(num_ae); - for (int i = 0; i < num_ae; i++) { - ae_density_h[i] = m_aerosol_mat[i]->m_density; - ae_solubility_h[i] = static_cast(m_aerosol_mat[i]->m_is_soluble); - } - Gpu::copy( Gpu::hostToDevice, - ae_density_h.begin(), - ae_density_h.end(), - ae_density.begin() ); - Gpu::copy( Gpu::hostToDevice, - ae_solubility_h.begin(), - ae_solubility_h.end(), - ae_solubility.begin() ); - } - - auto sp_rho_arr = sp_density.data(); - auto sp_sol_arr = sp_solubility.data(); - auto ae_rho_arr = ae_density.data(); - auto ae_sol_arr = ae_solubility.data(); + const auto x_min = m_recyc_xmin; + const auto x_max = m_recyc_xmax; + const auto y_min = m_recyc_ymin; + const auto y_max = m_recyc_ymax; + const auto z_min = m_recyc_zmin; + const auto z_max = m_recyc_zmax; + + forEachParticleTile(a_lev, ctx, + [&](ParIterType& /*pti*/, int /*grid*/, ParticleType* p_pbox, + const SDProcess::ParticlePointers& ptrs, + const SDProcess::ProcessContext& ctx) + { + const int np = ptrs.num_particles; - // get sampled aerosol mass values based on initialization - Gpu::DeviceVector aerosol_mass_d(num_ae*np); + // Get sampled aerosol mass values based on initialization + Gpu::DeviceVector aerosol_mass_d(ctx.num_aerosols*np); Gpu::DeviceVector multiplicity_d(np); ParticleReal mult_scale = one; { Vector multiplicity_h(np, zero); - for (int i = 0; i < num_ae; i++) { + for (int i = 0; i < ctx.num_aerosols; i++) { Vector aerosol_mass_h; if (sampled_multiplicity) { init_r.getAerosolDistribution( aerosol_mass_h, @@ -239,18 +176,11 @@ void SuperDropletPC::Recycle ( const int a_lev, Gpu::Buffer np_recycle_buf({0}); auto* np_recycle_ptr = np_recycle_buf.data(); - const auto x_min = m_recyc_xmin; - const auto x_max = m_recyc_xmax; - const auto y_min = m_recyc_ymin; - const auto y_max = m_recyc_ymax; - const auto z_min = m_recyc_zmin; - const auto z_max = m_recyc_zmax; - ParallelForRNG(np, [=] AMREX_GPU_DEVICE (int i, const RandomEngine& rnd_engine) noexcept { ParticleType& p = p_pbox[i]; if (p.id() <= 0) { return; } - if (active_ptr[i] > 0) { return; } + if (ptrs.active_ptr[i] > 0) { return; } // Place particle randomly in domain within specified bounds p.pos(0) = x_min + Random(rnd_engine)*(x_max - x_min); @@ -258,101 +188,88 @@ void SuperDropletPC::Recycle ( const int a_lev, p.pos(2) = z_min + Random(rnd_engine)*(z_max - z_min); // Set velocities to zero - v_ptr[0][i] = v_ptr[1][i] = v_ptr[2][i] = vterm_ptr[i] = zero; + ptrs.v_ptr[0][i] = ptrs.v_ptr[1][i] = ptrs.v_ptr[2][i] = ptrs.vterm_ptr[i] = zero; // reset all species masses to zero - for (int ctr = 0; ctr < num_sp; ctr++) { - sp_mass_ptrs[ctr][i] = zero; + for (int ctr = 0; ctr < ctx.num_species; ctr++) { + ptrs.sp_mass_ptrs[ctr][i] = zero; } // Reset water mass auto water_radius = Real(1.0e-15); - auto water_mass = (Real(4.0)/three)*PI - * water_radius*water_radius*water_radius*rho_w; - sp_mass_ptrs[idx_w][i] = water_mass; + auto water_mass = four_thirds_pi + * water_radius*water_radius*water_radius*ctx.rho_water; + ptrs.sp_mass_ptrs[ctx.idx_water][i] = water_mass; // choose a random index auto j = Random_int(np, rnd_engine); // Set aerosol mass - for (int ctr = 0; ctr < num_ae; ctr++) { - ae_mass_ptrs[ctr][i] = aerosol_mass[ctr*np+j]; + for (int ctr = 0; ctr < ctx.num_aerosols; ctr++) { + ptrs.ae_mass_ptrs[ctr][i] = aerosol_mass[ctr*np+j]; } // Set multiplicity to sampled or averaged multiplicity - if (sampled_multiplicity) { mult_ptr[i] = std::ceil(mult_arr[j] * mult_scale); } - else { mult_ptr[i] = mult_avg; } + if (sampled_multiplicity) { ptrs.mult_ptr[i] = std::ceil(mult_arr[j] * mult_scale); } + else { ptrs.mult_ptr[i] = mult_avg; } // compute effective radius and total mass - radius_ptr[i] = SD_effective_radius( i, idx_w, - rho_w, - num_sp, num_ae, - sp_sol_arr, ae_sol_arr, - sp_mass_ptrs, ae_mass_ptrs, - sp_rho_arr, ae_rho_arr ); - mass_ptr[i] = SD_total_mass( i, num_sp, num_ae, sp_mass_ptrs, ae_mass_ptrs); + SuperDropletPC::updateParticleAttributes( + i, ptrs.radius_ptr, ptrs.mass_ptr, ctx.idx_water, ctx.rho_water, + ctx.num_species, ctx.num_aerosols, ptrs.sp_sol_arr, ptrs.ae_sol_arr, + ptrs.sp_mass_ptrs, ptrs.ae_mass_ptrs, ptrs.sp_rho_arr, ptrs.ae_rho_arr); // set as active - active_ptr[i] = 1; + ptrs.active_ptr[i] = 1; // add to count Gpu::Atomic::Add(np_recycle_ptr, Long(1)); }); Gpu::synchronize(); np_recycle_total += *(np_recycle_buf.copyToHost()); - } + }); // end forEachParticleTile ParallelDescriptor::ReduceLongSum( &np_recycle_total, 1 ); Print() << " recycled " << np_recycle_total << " super-droplets.\n"; Redistribute(); + // Update location data for recycled particles const MFPtr& z_height = a_z_phys_nd[a_lev]; - const auto plo = Geom(m_lev).ProbLoArray(); - const auto dxi = Geom(m_lev).InvCellSizeArray(); - for (ParIterType pti(*this, a_lev); pti.isValid(); ++pti) { - int grid = pti.index(); - auto& ptile = ParticlesAt(a_lev, pti); - auto& aos = ptile.GetArrayOfStructs(); - const int n = aos.numParticles(); - auto *p_pbox = aos().data(); + forEachParticleTile(a_lev, ctx, + [&](ParIterType& /*pti*/, int grid, ParticleType* p_pbox, + const SDProcess::ParticlePointers& ptrs, + const SDProcess::ProcessContext& ctx) + { auto zheight = (*z_height)[grid].array(); - ParallelFor(n, [=] AMREX_GPU_DEVICE (int i) + ParallelFor(ptrs.num_particles, [=] AMREX_GPU_DEVICE (int i) { ParticleType& p = p_pbox[i]; if (p.id() <= 0) { return; } - update_location_idata(p,plo,dxi,zheight); - + update_location_idata(p,ctx.plo,ctx.dxi,zheight); }); Gpu::synchronize(); - } + }); // end forEachParticleTile } else { amrex::Print() << "Removing inactive particles.\n"; -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (ParIterType pti(*this, a_lev); pti.isValid(); ++pti) { - auto& ptile = ParticlesAt(a_lev, pti); - auto& aos = ptile.GetArrayOfStructs(); - auto& soa = ptile.GetStructOfArrays(); - const int np = aos.numParticles(); - auto *p_pbox = aos().data(); - - int rtoff_i = SuperDropletsIntIdxSoA::ncomps; - auto* active_ptr = soa.GetIntData(rtoff_i+SuperDropletsIntIdxSoA_RT::active).data(); + const auto ctx = buildProcessContext(a_lev); - ParallelForRNG(np, [=] AMREX_GPU_DEVICE (int i, const RandomEngine& /*rnd_engine*/) noexcept + forEachParticleTile(a_lev, ctx, + [&](ParIterType& /*pti*/, int /*grid*/, ParticleType* p_pbox, + const SDProcess::ParticlePointers& ptrs, + const SDProcess::ProcessContext& /*ctx*/) + { + ParallelForRNG(ptrs.num_particles, [=] AMREX_GPU_DEVICE (int i, const RandomEngine& /*rnd_engine*/) noexcept { ParticleType& p = p_pbox[i]; if (p.id() <= 0) { return; } - if (active_ptr[i] == 0) { p.id() = -1; } + if (ptrs.active_ptr[i] == 0) { p.id() = -1; } }); - } + }); // end forEachParticleTile Redistribute(); } } #endif - diff --git a/Source/Particles/ERF_SuperDropletPCUtils.cpp b/Source/Particles/ERF_SuperDropletPCUtils.cpp index 9503d90f1..b62e834b4 100644 --- a/Source/Particles/ERF_SuperDropletPCUtils.cpp +++ b/Source/Particles/ERF_SuperDropletPCUtils.cpp @@ -1,6 +1,5 @@ -#include -#include "ERF_Constants.H" #include "ERF_SuperDropletPC.H" +#include #ifdef ERF_USE_PARTICLES @@ -8,6 +7,60 @@ using namespace amrex; using namespace SDPCDefn; using SDTDType = typename SuperDropletPC::ParticleTileType::ConstParticleTileDataType; +/*! Initialize device property arrays for species and aerosol materials */ +void SuperDropletPC::initializeDeviceProperties() +{ + if (m_device_props_initialized) return; + + const int num_sp = m_species_mat.size(); + const int num_ae = m_aerosol_mat.size(); + + m_sp_density.resize(num_sp); + m_sp_solubility.resize(num_sp); + m_sp_ionization.resize(num_sp); + m_sp_mol_weight.resize(num_sp); + m_ae_density.resize(num_ae); + m_ae_solubility.resize(num_ae); + m_ae_ionization.resize(num_ae); + m_ae_mol_weight.resize(num_ae); + + amrex::Vector sp_density_h(num_sp), sp_ionization_h(num_sp), sp_mol_weight_h(num_sp); + amrex::Vector sp_solubility_h(num_sp); + amrex::Vector ae_density_h(num_ae), ae_ionization_h(num_ae), ae_mol_weight_h(num_ae); + amrex::Vector ae_solubility_h(num_ae); + + for (int i = 0; i < num_sp; i++) { + sp_density_h[i] = m_species_mat[i]->m_density; + sp_solubility_h[i] = static_cast(m_species_mat[i]->m_is_soluble); + sp_ionization_h[i] = m_species_mat[i]->m_ionization; + sp_mol_weight_h[i] = m_species_mat[i]->m_mol_weight; + } + for (int i = 0; i < num_ae; i++) { + ae_density_h[i] = m_aerosol_mat[i]->m_density; + ae_solubility_h[i] = static_cast(m_aerosol_mat[i]->m_is_soluble); + ae_ionization_h[i] = m_aerosol_mat[i]->m_ionization; + ae_mol_weight_h[i] = m_aerosol_mat[i]->m_mol_weight; + } + + Gpu::copy(Gpu::hostToDevice, sp_density_h.begin(), sp_density_h.end(), m_sp_density.begin()); + Gpu::copy(Gpu::hostToDevice, sp_solubility_h.begin(), sp_solubility_h.end(), m_sp_solubility.begin()); + Gpu::copy(Gpu::hostToDevice, sp_ionization_h.begin(), sp_ionization_h.end(), m_sp_ionization.begin()); + Gpu::copy(Gpu::hostToDevice, sp_mol_weight_h.begin(), sp_mol_weight_h.end(), m_sp_mol_weight.begin()); + Gpu::copy(Gpu::hostToDevice, ae_density_h.begin(), ae_density_h.end(), m_ae_density.begin()); + Gpu::copy(Gpu::hostToDevice, ae_solubility_h.begin(), ae_solubility_h.end(), m_ae_solubility.begin()); + Gpu::copy(Gpu::hostToDevice, ae_ionization_h.begin(), ae_ionization_h.end(), m_ae_ionization.begin()); + Gpu::copy(Gpu::hostToDevice, ae_mol_weight_h.begin(), ae_mol_weight_h.end(), m_ae_mol_weight.begin()); + + m_device_props_initialized = true; +} + +/*! Update device properties if material properties change */ +void SuperDropletPC::updateDeviceProperties() +{ + m_device_props_initialized = false; + initializeDeviceProperties(); +} + /*! Compute mesh variable from particles */ void SuperDropletPC::computeMeshVar( const std::string& a_var_name, MultiFab& a_mf, @@ -15,86 +68,74 @@ void SuperDropletPC::computeMeshVar( const std::string& a_var_name, { BL_PROFILE("SuperDropletPC::computeMeshVar()"); a_mf.setVal(0.0); - if (a_lev == 0) { - if (a_var_name == "number_density") { - numberDensity( a_mf ); - } else if (a_var_name == "sd_number_density") { - SDNumberDensity( a_mf ); - } else if (a_var_name == "mass_density") { - massDensity( a_mf ); - } else if (a_var_name == ("mass_density_"+getEnumNameString(m_species_mat[m_idx_w]->m_name))) { - speciesMassDensity( a_mf, m_idx_w ); - } else if (a_var_name == ("mass_flux_x_" +getEnumNameString(m_species_mat[m_idx_w]->m_name))) { - speciesMassFlux( a_mf, m_idx_w, 0 ); - } else if (a_var_name == ("mass_flux_y_" +getEnumNameString(m_species_mat[m_idx_w]->m_name))) { - speciesMassFlux( a_mf, m_idx_w, 1 ); - } else if (a_var_name == ("mass_flux_z_" +getEnumNameString(m_species_mat[m_idx_w]->m_name))) { - speciesMassFlux( a_mf, m_idx_w, 2 ); - } else if (a_var_name == "mass_flux_x") { - massFlux( a_mf, 0 ); - } else if (a_var_name == "mass_flux_y") { - massFlux( a_mf, 1 ); - } else if (a_var_name == "mass_flux_z") { - massFlux( a_mf, 2 ); - } else if (a_var_name == "radius") { - effectiveRadius( a_mf ); - } else { - for (int i = 0; i < m_num_species; i++) { - std::string var_name = "species_mass_density_"+getEnumNameString(m_species_mat[i]->m_name); - if (a_var_name == var_name) { - speciesMassDensity( a_mf, i ); - break; - } - } - for (int i = 0; i < m_num_species; i++) { - std::string var_name = "species_mass_flux_x_" +getEnumNameString(m_species_mat[i]->m_name); - if (a_var_name == var_name) { - speciesMassFlux( a_mf, i, 0 ); - break; - } - } - for (int i = 0; i < m_num_species; i++) { - std::string var_name = "species_mass_flux_y_" +getEnumNameString(m_species_mat[i]->m_name); - if (a_var_name == var_name) { - speciesMassFlux( a_mf, i, 1 ); - break; - } - } - for (int i = 0; i < m_num_species; i++) { - std::string var_name = "species_mass_flux_z_" +getEnumNameString(m_species_mat[i]->m_name); - if (a_var_name == var_name) { - speciesMassFlux( a_mf, i, 2 ); - break; - } - } - for (int i = 0; i < m_num_aerosols; i++) { - std::string var_name = "aerosol_mass_density_"+getEnumNameString(m_aerosol_mat[i]->m_name); - if (a_var_name == var_name) { - aerosolMassDensity( a_mf, i ); - break; - } - } - for (int i = 0; i < m_num_aerosols; i++) { - std::string var_name = "aerosol_mass_flux_x_"+getEnumNameString(m_aerosol_mat[i]->m_name); - if (a_var_name == var_name) { - aerosolMassFlux( a_mf, i, 0 ); - break; - } - } - for (int i = 0; i < m_num_aerosols; i++) { - std::string var_name = "aerosol_mass_flux_y_"+getEnumNameString(m_aerosol_mat[i]->m_name); - if (a_var_name == var_name) { - aerosolMassFlux( a_mf, i, 1 ); - break; - } - } - for (int i = 0; i < m_num_aerosols; i++) { - std::string var_name = "aerosol_mass_flux_z_"+getEnumNameString(m_aerosol_mat[i]->m_name); - if (a_var_name == var_name) { - aerosolMassFlux( a_mf, i, 2 ); - break; - } - } + + if (a_lev != 0) { return; } + + // Check basic variables + if (a_var_name == "number_density") { + numberDensity(a_mf); return; + } + if (a_var_name == "sd_number_density") { + SDNumberDensity(a_mf); return; + } + if (a_var_name == "mass_density") { + massDensity(a_mf); return; + } + if (a_var_name == "radius") { + effectiveRadius(a_mf); return; + } + + // Check water species shortcuts + const std::string water_name = getEnumNameString(m_species_mat[m_idx_w]->m_name); + if (a_var_name == "mass_density_" + water_name) { + speciesMassDensity(a_mf, m_idx_w); return; + } + if (a_var_name == "mass_flux_x_" + water_name) { + speciesMassFlux(a_mf, m_idx_w, 0); return; + } + if (a_var_name == "mass_flux_y_" + water_name) { + speciesMassFlux(a_mf, m_idx_w, 1); return; + } + if (a_var_name == "mass_flux_z_" + water_name) { + speciesMassFlux(a_mf, m_idx_w, 2); return; + } + + // Check total mass flux + if (a_var_name == "mass_flux_x") { massFlux(a_mf, 0); return; } + if (a_var_name == "mass_flux_y") { massFlux(a_mf, 1); return; } + if (a_var_name == "mass_flux_z") { massFlux(a_mf, 2); return; } + + // Check species variables + for (int i = 0; i < m_num_species; i++) { + const std::string name = getEnumNameString(m_species_mat[i]->m_name); + if (a_var_name == "species_mass_density_" + name) { + speciesMassDensity(a_mf, i); return; + } + if (a_var_name == "species_mass_flux_x_" + name) { + speciesMassFlux(a_mf, i, 0); return; + } + if (a_var_name == "species_mass_flux_y_" + name) { + speciesMassFlux(a_mf, i, 1); return; + } + if (a_var_name == "species_mass_flux_z_" + name) { + speciesMassFlux(a_mf, i, 2); return; + } + } + + // Check aerosol variables + for (int i = 0; i < m_num_aerosols; i++) { + const std::string name = getEnumNameString(m_aerosol_mat[i]->m_name); + if (a_var_name == "aerosol_mass_density_" + name) { + aerosolMassDensity(a_mf, i); return; + } + if (a_var_name == "aerosol_mass_flux_x_" + name) { + aerosolMassFlux(a_mf, i, 0); return; + } + if (a_var_name == "aerosol_mass_flux_y_" + name) { + aerosolMassFlux(a_mf, i, 1); return; + } + if (a_var_name == "aerosol_mass_flux_z_" + name) { + aerosolMassFlux(a_mf, i, 2); return; } } } @@ -130,365 +171,178 @@ Long SuperDropletPC::NumSDDeactivated () } /*! Computes the number density of the SDs over a mesh */ -void SuperDropletPC::SDNumberDensity ( MultiFab& a_mf, /*!< Number density multifab */ - const int a_comp /*!< MultiFab component to fill with number density */) const +void SuperDropletPC::SDNumberDensity ( MultiFab& a_mf, + const int a_comp ) const { BL_PROFILE("SuperDropletPC::SDNumberDensity()"); - - AMREX_ASSERT(OK()); - AMREX_ASSERT(numParticlesOutOfRange(*this, 0) == 0); - - const auto& geom = Geom(m_lev); - const auto plo = geom.ProbLoArray(); - const auto dxi = geom.InvCellSizeArray(); - const ParticleReal inv_cell_volume = dxi[0]*dxi[1]*dxi[2]; - - a_mf.setVal(0.0); - - ParticleToMesh( *this, a_mf, m_lev, - [=] AMREX_GPU_DEVICE ( const SuperDropletPC::ParticleTileType::ConstParticleTileDataType& ptd, - int i, Array4 const& rho) - { - auto p = ptd.m_aos[i]; - ParticleInterpolator::Linear interp(p, plo, dxi); - interp.ParticleToMesh ( p, rho, 0, a_comp, 1, - [=] AMREX_GPU_DEVICE ( const SuperDropletPC::ParticleType&, int) - { - auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; - return ai*inv_cell_volume; - }); + ERFPCParticleToMesh(a_mf, m_lev, a_comp, + [=] AMREX_GPU_DEVICE (const SDTDType& ptd, int i) { + return ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; }); - - return; } /*! Computes the number density of the particles over a mesh */ -void SuperDropletPC::numberDensity ( MultiFab& a_mf, /*!< Number density multifab */ - const int a_comp /*!< MultiFab component to fill with number density */) const +void SuperDropletPC::numberDensity ( MultiFab& a_mf, + const int a_comp ) const { BL_PROFILE("SuperDropletPC::numberDensity()"); - - AMREX_ASSERT(OK()); - AMREX_ASSERT(numParticlesOutOfRange(*this, 0) == 0); - - const auto& geom = Geom(m_lev); - const auto plo = geom.ProbLoArray(); - const auto dxi = geom.InvCellSizeArray(); - - const ParticleReal inv_cell_volume = dxi[0]*dxi[1]*dxi[2]; - - a_mf.setVal(0.0); - - ParticleToMesh( *this, a_mf, m_lev, - [=] AMREX_GPU_DEVICE ( const SuperDropletPC::ParticleTileType::ConstParticleTileDataType& ptd, - int i, Array4 const& rho) - { - auto p = ptd.m_aos[i]; - ParticleInterpolator::Linear interp(p, plo, dxi); - interp.ParticleToMesh ( p, rho, 0, a_comp, 1, - [=] AMREX_GPU_DEVICE ( const SuperDropletPC::ParticleType&, int) - { - auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; - auto num_par = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; - return ai*num_par*inv_cell_volume; - }); + ERFPCParticleToMesh(a_mf, m_lev, a_comp, + [=] AMREX_GPU_DEVICE (const SDTDType& ptd, int i) { + auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; + auto num_par = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; + return ai * num_par; }); - - return; } -/*! Computes the mass density of the particles over a mesh: this does - include the aerosol mass*/ -void SuperDropletPC::massDensity ( MultiFab& a_mf, /*!< Mass density multifab */ - const int a_comp /*!< MultiFab component to fill with mass density */) const +/*! Computes the mass density of the particles over a mesh */ +void SuperDropletPC::massDensity ( MultiFab& a_mf, + const int a_comp ) const { BL_PROFILE("SuperDropletPC::massDensity()"); - - AMREX_ASSERT(OK()); - AMREX_ASSERT(numParticlesOutOfRange(*this, 0) == 0); - - const auto& geom = Geom(m_lev); - const auto plo = geom.ProbLoArray(); - const auto dxi = geom.InvCellSizeArray(); - - const ParticleReal inv_cell_volume = dxi[0]*dxi[1]*dxi[2]; - a_mf.setVal(0.0); - - ParticleToMesh( *this, a_mf, m_lev, - [=] AMREX_GPU_DEVICE ( const SuperDropletPC::ParticleTileType::ConstParticleTileDataType& ptd, - int i, Array4 const& rho) - { - auto p = ptd.m_aos[i]; - ParticleInterpolator::Linear interp(p, plo, dxi); - interp.ParticleToMesh ( p, rho, 0, a_comp, 1, - [=] AMREX_GPU_DEVICE ( const SuperDropletPC::ParticleType&, int) - { - auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; - auto num_par = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; - auto par_mass = ptd.m_rdata[SuperDropletsRealIdxSoA::mass][i]; - return ai*num_par*par_mass*inv_cell_volume; - }); + ERFPCParticleToMesh(a_mf, m_lev, a_comp, + [=] AMREX_GPU_DEVICE (const SDTDType& ptd, int i) { + auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; + auto num_par = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; + auto par_mass = ptd.m_rdata[SuperDropletsRealIdxSoA::mass][i]; + return ai * num_par * par_mass; }); - - return; } /*! Computes the particle velocity components over a mesh */ -void SuperDropletPC::massFlux ( MultiFab& a_mf, /*!< Mass flux multifab */ - const int a_dim, /*!< Flux component */ - const int a_comp /*!< MultiFab component to fill with mass density */) const +void SuperDropletPC::massFlux ( MultiFab& a_mf, + const int a_dim, + const int a_comp ) const { BL_PROFILE("SuperDropletPC::massFlux()"); - - AMREX_ASSERT(OK()); - AMREX_ASSERT(numParticlesOutOfRange(*this, 0) == 0); - - const auto& geom = Geom(m_lev); - const auto plo = geom.ProbLoArray(); - const auto dxi = geom.InvCellSizeArray(); - - const ParticleReal inv_cell_volume = dxi[0]*dxi[1]*dxi[2]; - a_mf.setVal(0.0); - - ParticleToMesh( *this, a_mf, m_lev, - [=] AMREX_GPU_DEVICE ( const SuperDropletPC::ParticleTileType::ConstParticleTileDataType& ptd, - int i, Array4 const& rho) - { - auto p = ptd.m_aos[i]; - ParticleInterpolator::Linear interp(p, plo, dxi); - interp.ParticleToMesh ( p, rho, 0, a_comp, 1, - [=] AMREX_GPU_DEVICE ( const SuperDropletPC::ParticleType&, int) - { - auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; - auto num_par = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; - auto par_mass = ptd.m_rdata[SuperDropletsRealIdxSoA::mass][i]; - auto par_velocity = ptd.m_rdata[SuperDropletsRealIdxSoA::vx+a_dim][i]; - if (a_dim == 2) { - auto term_vel = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::term_vel][i]; - par_velocity -= term_vel; - } - return ai * num_par * par_mass * par_velocity * inv_cell_volume; - }); + ERFPCParticleToMesh(a_mf, m_lev, a_comp, + [=] AMREX_GPU_DEVICE (const SDTDType& ptd, int i) { + auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; + auto num_par = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; + auto par_mass = ptd.m_rdata[SuperDropletsRealIdxSoA::mass][i]; + auto par_velocity = ptd.m_rdata[SuperDropletsRealIdxSoA::vx+a_dim][i]; + if (a_dim == 2) { + par_velocity -= ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::term_vel][i]; + } + return ai * num_par * par_mass * par_velocity; }); - - return; } /*! Computes the aerosol mass density of the particles over a mesh */ -void SuperDropletPC::aerosolMassDensity ( MultiFab& a_mf, /*!< Aerosol mass density multifab */ - const int a_idx, /*!< Aerosol index */ - const int a_comp /*!< MultiFab component to fill */) const +void SuperDropletPC::aerosolMassDensity ( MultiFab& a_mf, + const int a_idx, + const int a_comp ) const { BL_PROFILE("SuperDropletPC::aerosolMassDensity()"); - - AMREX_ASSERT(OK()); - AMREX_ASSERT(numParticlesOutOfRange(*this, 0) == 0); - - const auto& geom = Geom(m_lev); - const auto plo = geom.ProbLoArray(); - const auto dxi = geom.InvCellSizeArray(); - - const ParticleReal inv_cell_volume = dxi[0]*dxi[1]*dxi[2]; - - auto na = m_num_aerosols; - auto ns = m_num_species; - - a_mf.setVal(0.0); - - ParticleToMesh( *this, a_mf, m_lev, - [=] AMREX_GPU_DEVICE ( const SuperDropletPC::ParticleTileType::ConstParticleTileDataType& ptd, - int i, Array4 const& rho) - { - auto p = ptd.m_aos[i]; - ParticleInterpolator::Linear interp(p, plo, dxi); - interp.ParticleToMesh ( p, rho, 0, a_comp, 1, - [=] AMREX_GPU_DEVICE ( const SuperDropletPC::ParticleType&, int) - { - auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; - auto num_par = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; - auto aero_mass = ptd.m_runtime_rdata[ridx_a(a_idx,na,ns)][i]; - return ai*num_par*aero_mass*inv_cell_volume; - }); + const auto na = m_num_aerosols; + const auto ns = m_num_species; + ERFPCParticleToMesh(a_mf, m_lev, a_comp, + [=] AMREX_GPU_DEVICE (const SDTDType& ptd, int i) { + auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; + auto num_par = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; + auto aero_mass = ptd.m_runtime_rdata[ridx_a(a_idx,na,ns)][i]; + return ai * num_par * aero_mass; }); - - return; } -/*! Computes the aerosol mass flux of the partices over a mesh */ -void SuperDropletPC::aerosolMassFlux ( MultiFab& a_mf, /*!< Aerosol mass flux multifab */ - const int a_idx, /*!< Aerosol index */ - const int a_dim, /*!< Flux component */ - const int a_comp /*!< MultiFab component to fill */) const +/*! Computes the aerosol mass flux of the particles over a mesh */ +void SuperDropletPC::aerosolMassFlux ( MultiFab& a_mf, + const int a_idx, + const int a_dim, + const int a_comp ) const { BL_PROFILE("SuperDropletPC::aerosolMassFlux()"); - - AMREX_ASSERT(OK()); - AMREX_ASSERT(numParticlesOutOfRange(*this, 0) == 0); - - const auto& geom = Geom(m_lev); - const auto plo = geom.ProbLoArray(); - const auto dxi = geom.InvCellSizeArray(); - - const ParticleReal inv_cell_volume = dxi[0]*dxi[1]*dxi[2]; - - auto na = m_num_aerosols; - auto ns = m_num_species; - - a_mf.setVal(0.0); - - ParticleToMesh( *this, a_mf, m_lev, - [=] AMREX_GPU_DEVICE ( const SuperDropletPC::ParticleTileType::ConstParticleTileDataType& ptd, - int i, Array4 const& rho) - { - auto p = ptd.m_aos[i]; - ParticleInterpolator::Linear interp(p, plo, dxi); - interp.ParticleToMesh ( p, rho, 0, a_comp, 1, - [=] AMREX_GPU_DEVICE ( const SuperDropletPC::ParticleType&, int) - { - auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; - auto num_par = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; - auto aero_mass = ptd.m_runtime_rdata[ridx_a(a_idx,na,ns)][i]; - auto par_velocity = ptd.m_rdata[SuperDropletsRealIdxSoA::vx+a_dim][i]; - if (a_dim == 2) { - auto term_vel = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::term_vel][i]; - par_velocity -= term_vel; - } - return ai*num_par*aero_mass*par_velocity*inv_cell_volume; - }); + const auto na = m_num_aerosols; + const auto ns = m_num_species; + ERFPCParticleToMesh(a_mf, m_lev, a_comp, + [=] AMREX_GPU_DEVICE (const SDTDType& ptd, int i) { + auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; + auto num_par = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; + auto aero_mass = ptd.m_runtime_rdata[ridx_a(a_idx,na,ns)][i]; + auto par_velocity = ptd.m_rdata[SuperDropletsRealIdxSoA::vx+a_dim][i]; + if (a_dim == 2) { + par_velocity -= ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::term_vel][i]; + } + return ai * num_par * aero_mass * par_velocity; }); - - return; } /*! Computes the species mass density of the particles over a mesh */ -void SuperDropletPC::speciesMassDensity ( MultiFab& a_mf, /*!< Species mass density multifab */ - const int a_idx, /*!< Species index */ - const Real a_rmin, /*!< minimum radius */ - const Real a_rmax, /*!< maximum radius */ - const int a_comp /*!< MultiFab component to fill */) const +void SuperDropletPC::speciesMassDensity ( MultiFab& a_mf, + int a_idx, + const int a_comp ) const { BL_PROFILE("SuperDropletPC::speciesMassDensity()"); - - AMREX_ASSERT(OK()); - AMREX_ASSERT(numParticlesOutOfRange(*this, 0) == 0); - - const auto& geom = Geom(m_lev); - const auto plo = geom.ProbLoArray(); - const auto dxi = geom.InvCellSizeArray(); - - const ParticleReal inv_cell_volume = dxi[0]*dxi[1]*dxi[2]; - - auto na = m_num_aerosols; - auto ns = m_num_species; - - a_mf.setVal(0.0); - - ParticleToMesh( *this, a_mf, m_lev, - [=] AMREX_GPU_DEVICE ( const SuperDropletPC::ParticleTileType::ConstParticleTileDataType& ptd, - int i, Array4 const& rho) - { - auto p = ptd.m_aos[i]; - ParticleInterpolator::Linear interp(p, plo, dxi); - interp.ParticleToMesh ( p, rho, 0, a_comp, 1, - [=] AMREX_GPU_DEVICE ( const SuperDropletPC::ParticleType&, int) - { - auto radius = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::radius][i]; - if ((radius < a_rmin) || (radius >= a_rmax)) { - return zero; - } else { - auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; - auto num_par = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; - auto species_mass = ptd.m_runtime_rdata[ridx_s(a_idx,na,ns)][i]; - return ai*num_par*species_mass*inv_cell_volume; - } - }); + const auto na = m_num_aerosols; + const auto ns = m_num_species; + ERFPCParticleToMesh(a_mf, m_lev, a_comp, + [=] AMREX_GPU_DEVICE (const SDTDType& ptd, int i) { + auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; + auto num_par = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; + auto species_mass = ptd.m_runtime_rdata[ridx_s(a_idx,na,ns)][i]; + return ai * num_par * species_mass; }); +} - return; +/*! Computes the cloud/rain mass density of the particles over a mesh */ +void SuperDropletPC::cloudRainDensity(MultiFab& a_mf, const Real a_rmin, const Real a_rmax, const int a_comp) const +{ + BL_PROFILE("SuperDropletPC::cloudRainDensity()"); + const auto na = m_num_aerosols; + const auto ns = m_num_species; + const auto idx = m_idx_w; + ERFPCParticleToMesh(a_mf, m_lev, a_comp, + [=] AMREX_GPU_DEVICE (const SDTDType& ptd, int i) { + auto radius = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::radius][i]; + if ((radius < a_rmin) || (radius >= a_rmax)) { + return 0.0; + } + auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; + auto num_par = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; + auto species_mass = ptd.m_runtime_rdata[ridx_s(idx,na,ns)][i]; + return ai * num_par * species_mass; + }); } /*! Computes the species mass flux of the particles over a mesh */ -void SuperDropletPC::speciesMassFlux ( MultiFab& a_mf, /*!< Species mass flux multifab */ - const int a_idx, /*!< Species index */ - const int a_dim, /*!< Flux component */ - const int a_comp /*!< MultiFab component to fill */) const +void SuperDropletPC::speciesMassFlux ( MultiFab& a_mf, + const int a_idx, + const int a_dim, + const int a_comp ) const { BL_PROFILE("SuperDropletPC::speciesMassFlux()"); - - AMREX_ASSERT(OK()); - AMREX_ASSERT(numParticlesOutOfRange(*this, 0) == 0); - - const auto& geom = Geom(m_lev); - const auto plo = geom.ProbLoArray(); - const auto dxi = geom.InvCellSizeArray(); - - const ParticleReal inv_cell_volume = dxi[0]*dxi[1]*dxi[2]; - - auto na = m_num_aerosols; - auto ns = m_num_species; - - a_mf.setVal(0.0); - - ParticleToMesh( *this, a_mf, m_lev, - [=] AMREX_GPU_DEVICE ( const SuperDropletPC::ParticleTileType::ConstParticleTileDataType& ptd, - int i, Array4 const& rho) - { - auto p = ptd.m_aos[i]; - ParticleInterpolator::Linear interp(p, plo, dxi); - interp.ParticleToMesh ( p, rho, 0, a_comp, 1, - [=] AMREX_GPU_DEVICE ( const SuperDropletPC::ParticleType&, int) - { - auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; - auto num_par = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; - auto species_mass = ptd.m_runtime_rdata[ridx_s(a_idx,na,ns)][i]; - auto par_velocity = ptd.m_rdata[SuperDropletsRealIdxSoA::vx+a_dim][i]; - if (a_dim == 2) { - auto term_vel = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::term_vel][i]; - par_velocity -= term_vel; - } - return ai*num_par*species_mass*par_velocity*inv_cell_volume; - }); + const auto na = m_num_aerosols; + const auto ns = m_num_species; + ERFPCParticleToMesh(a_mf, m_lev, a_comp, + [=] AMREX_GPU_DEVICE (const SDTDType& ptd, int i) { + auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; + auto num_par = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; + auto species_mass = ptd.m_runtime_rdata[ridx_s(a_idx,na,ns)][i]; + auto par_velocity = ptd.m_rdata[SuperDropletsRealIdxSoA::vx+a_dim][i]; + if (a_dim == 2) { + par_velocity -= ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::term_vel][i]; + } + return ai * num_par * species_mass * par_velocity; }); - - return; } /*! Computes the effective radius of the particles over a mesh */ -void SuperDropletPC::effectiveRadius ( MultiFab& a_mf, /*!< Effective radius multifab */ - const int a_comp /*!< MultiFab component to fill with number density */) const +void SuperDropletPC::effectiveRadius ( MultiFab& a_mf, + const int a_comp ) const { BL_PROFILE("SuperDropletPC::effectiveRadius()"); - AMREX_ASSERT(OK()); - AMREX_ASSERT(numParticlesOutOfRange(*this, 0) == 0); - - const auto& geom = Geom(m_lev); - const auto plo = geom.ProbLoArray(); - const auto dxi = geom.InvCellSizeArray(); - const ParticleReal inv_cell_volume = dxi[0]*dxi[1]*dxi[2]; - - a_mf.setVal(0.0); - - MultiFab number_density( a_mf.boxArray(), a_mf.DistributionMap(), 1, a_mf.nGrowVect() ); + MultiFab number_density(a_mf.boxArray(), a_mf.DistributionMap(), 1, a_mf.nGrowVect()); numberDensity(number_density); - ParticleToMesh( *this, a_mf, m_lev, - [=] AMREX_GPU_DEVICE ( const SuperDropletPC::ParticleTileType::ConstParticleTileDataType& ptd, - int i, Array4 const& rho) - { - auto p = ptd.m_aos[i]; - ParticleInterpolator::Linear interp(p, plo, dxi); - interp.ParticleToMesh ( p, rho, 0, a_comp, 1, - [=] AMREX_GPU_DEVICE ( const SuperDropletPC::ParticleType&, int) - { - auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; - auto num_par = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; - auto radius = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::radius][i]; - return ai*num_par*radius*inv_cell_volume; - }); + ERFPCParticleToMesh(a_mf, m_lev, a_comp, + [=] AMREX_GPU_DEVICE (const SDTDType& ptd, int i) { + auto ai = ptd.m_runtime_idata[SuperDropletsIntIdxSoA_RT::active][i]; + auto num_par = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::multiplicity][i]; + auto radius = ptd.m_runtime_rdata[SuperDropletsRealIdxSoA_RT::radius][i]; + return ai * num_par * radius; }); - for ( MFIter mfi(a_mf); mfi.isValid(); ++mfi) { + for (MFIter mfi(a_mf); mfi.isValid(); ++mfi) { const auto& box = mfi.tilebox(); auto mf_arr = a_mf.array(mfi); const auto nd_arr = number_density.const_array(mfi); @@ -501,11 +355,6 @@ void SuperDropletPC::effectiveRadius ( MultiFab& a_mf, /*!< Effective radius m } } ); } - - - return; } #endif - - diff --git a/Source/Particles/Make.package b/Source/Particles/Make.package index f8c659029..679fce7c2 100644 --- a/Source/Particles/Make.package +++ b/Source/Particles/Make.package @@ -15,8 +15,10 @@ CEXE_sources += ERF_SuperDropletPCUtils.cpp CEXE_headers += ERF_ParticleData.H CEXE_headers += ERFPC.H +CEXE_headers += ERFPCParticleToMesh.H CEXE_headers += ERF_SDInitialization.H CEXE_headers += ERF_SuperDropletPC.H CEXE_headers += ERF_SuperDropletPCDefinitions.H CEXE_headers += ERF_SuperDropletPCCoalescence.H CEXE_headers += ERF_SuperDropletPCMassChange.H +CEXE_headers += ERF_SuperDropletPCProcess.H diff --git a/Source/Utils/ERF_InterpolationUtils.H b/Source/Utils/ERF_InterpolationUtils.H new file mode 100644 index 000000000..3745c534d --- /dev/null +++ b/Source/Utils/ERF_InterpolationUtils.H @@ -0,0 +1,60 @@ +#ifndef ERF_INTERPOLATION_UTILS_H +#define ERF_INTERPOLATION_UTILS_H + +#include +#include +#include +#include +#include "ERF_Constants.H" + +/** + * @file ERF_InterpolationUtils.H + * @brief Utility functions for field interpolation at particle positions + */ + +namespace ERF { +namespace Interpolation { + +/** + * @brief Helper function to interpolate multiple field values at a particle position + * + * This function interpolates values from multiple field arrays at a particle's position + * using the cloud-in-cell (CIC) interpolation method. It supports both periodic and + * non-periodic z-boundaries. + * + * @tparam ParticleType The type of particle + * @param p The particle + * @param plo The domain lower bounds + * @param dxi The inverse of cell spacing + * @param field_arrays Array of field arrays to interpolate from + * @param result_values Output array to store interpolated values + * @param num_fields Number of fields to interpolate + * @param is_periodic_z Whether the z-boundary is periodic + * @param zheight Height field array (only used for non-periodic z-boundary) + */ +template +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void interpolateFields( + const ParticleType& p, + const amrex::GpuArray& plo, + const amrex::GpuArray& dxi, + const amrex::Array4* field_arrays, + amrex::ParticleReal* result_values, + int num_fields, + int is_periodic_z = 1, + const amrex::Array4* zheight = nullptr +) { + // Iterate over each field and interpolate its value at the particle position + for (int i = 0; i < num_fields; ++i) { + if (is_periodic_z != 0) { + cic_interpolate(p, plo, dxi, field_arrays[i], &result_values[i], 1); + } else { + cic_interpolate_mapped_z(p, plo, dxi, field_arrays[i], *zheight, &result_values[i], 1); + } + } +} + +} // namespace Interpolation +} // namespace ERF + +#endif // ERF_INTERPOLATION_UTILS_H