From 522382eb2f46aff1ca9253825d5064d99049a926 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Wed, 11 Oct 2023 16:43:37 -0700 Subject: [PATCH] Tweaks (#1264) * update particle functionality * some minor cleanup --- Source/ERF.cpp | 42 ++++++++++++--------- Source/IO/ERF_WriteScalarProfiles.cpp | 39 +++++++------------ Source/IO/NCColumnFile.cpp | 11 ++++-- Source/Initialization/ERF_init_custom.cpp | 19 +++++++--- Source/Microphysics/Init.cpp | 10 ----- Source/Microphysics/Microphysics.H | 3 -- Source/Microphysics/Precip.cpp | 10 ----- Source/TimeIntegration/ERF_slow_rhs_inc.cpp | 19 +++++----- 8 files changed, 69 insertions(+), 84 deletions(-) diff --git a/Source/ERF.cpp b/Source/ERF.cpp index a5ebc4af1..2e54eded7 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -421,34 +421,42 @@ ERF::InitData () amrex::Abort("We do not currently support init_type = ideal or input_sounding with terrain"); } - // We initialize rho_KE to be nonzero (and positive) so that we end up + // If using the Deardoff LES model, + // we initialize rho_KE to be nonzero (and positive) so that we end up // with reasonable values for the eddy diffusivity and the MOST fluxes // (~ 1/diffusivity) do not blow up Real RhoKE_0 = 0.1; ParmParse pp(pp_prefix); if (pp.query("RhoKE_0", RhoKE_0)) { - // uniform initial rho*e field - int lb = std::max(finest_level-1,0); - for (int lev(lb); lev >= 0; --lev) - vars_new[lev][Vars::cons].setVal(RhoKE_0,RhoKE_comp,1,0); + // Uniform initial rho*e field + for (int lev = 0; lev <= finest_level; lev++) { + if (solverChoice.turbChoice[lev].les_type == LESType::Deardorff) { + vars_new[lev][Vars::cons].setVal(RhoKE_0,RhoKE_comp,1,0); + } else { + vars_new[lev][Vars::cons].setVal(0.0,RhoKE_comp,1,0); + } + } } else { // default: uniform initial e field Real KE_0 = 0.1; pp.query("KE_0", KE_0); - for (int lev = 0; lev <= finest_level; lev++) - { + for (int lev = 0; lev <= finest_level; lev++) { auto& lev_new = vars_new[lev]; - for (MFIter mfi(lev_new[Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Box &bx = mfi.tilebox(); - const auto &cons_arr = lev_new[Vars::cons].array(mfi); - // We want to set the lateral BC values, too - Box gbx = bx; // Copy constructor - gbx.grow(0,1); gbx.grow(1,1); // Grow by one in the lateral directions - amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - cons_arr(i,j,k,RhoKE_comp) = cons_arr(i,j,k,Rho_comp) * KE_0; - }); + if (solverChoice.turbChoice[lev].les_type == LESType::Deardorff) { + for (MFIter mfi(lev_new[Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const Box &bx = mfi.tilebox(); + const auto &cons_arr = lev_new[Vars::cons].array(mfi); + // We want to set the lateral BC values, too + Box gbx = bx; // Copy constructor + gbx.grow(0,1); gbx.grow(1,1); // Grow by one in the lateral directions + amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + cons_arr(i,j,k,RhoKE_comp) = cons_arr(i,j,k,Rho_comp) * KE_0; + }); + } // mfi + } else { + lev_new[Vars::cons].setVal(0.0,RhoKE_comp,1,0); } - } + } // lev } if (coupling_type == "TwoWay") { diff --git a/Source/IO/ERF_WriteScalarProfiles.cpp b/Source/IO/ERF_WriteScalarProfiles.cpp index a88c214e2..78521d620 100644 --- a/Source/IO/ERF_WriteScalarProfiles.cpp +++ b/Source/IO/ERF_WriteScalarProfiles.cpp @@ -293,36 +293,25 @@ ERF::volWgtSumMF(int lev, amrex::MultiFab& ERF::build_fine_mask(int level) { - // Mask for zeroing covered cells - AMREX_ASSERT(level > 0); + // Mask for zeroing covered cells + AMREX_ASSERT(level > 0); - const amrex::BoxArray& cba = grids[level-1]; - const amrex::DistributionMapping& cdm = dmap[level-1]; + const amrex::BoxArray& cba = grids[level-1]; + const amrex::DistributionMapping& cdm = dmap[level-1]; - // TODO -- we should make a vector of these a member of ERF class - fine_mask.define(cba, cdm, 1, 0, amrex::MFInfo()); - fine_mask.setVal(1.0); + // TODO -- we should make a vector of these a member of ERF class + fine_mask.define(cba, cdm, 1, 0, amrex::MFInfo()); + fine_mask.setVal(1.0); - amrex::BoxArray fba = grids[level]; - amrex::iMultiFab ifine_mask = makeFineMask(cba, cdm, fba, ref_ratio[level-1], 1, 0); + amrex::BoxArray fba = grids[level]; + amrex::iMultiFab ifine_mask = makeFineMask(cba, cdm, fba, ref_ratio[level-1], 1, 0); -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for (amrex::MFIter mfi(fine_mask, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) + const auto fma = fine_mask.arrays(); + const auto ifma = ifine_mask.arrays(); + amrex::ParallelFor(fine_mask, [=] AMREX_GPU_DEVICE(int bno, int i, int j, int k) noexcept { - auto& fab = fine_mask[mfi]; - auto& ifab = ifine_mask[mfi]; - const auto arr = fab.array(); - const auto iarr = ifab.array(); - amrex::ParallelFor( - fab.box(), [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { -#ifdef _OPENMP -#pragma omp atomic write -#endif - arr(i, j, k) = iarr(i, j, k); - }); - } + fma[bno](i,j,k) = ifma[bno](i,j,k); + }); return fine_mask; } diff --git a/Source/IO/NCColumnFile.cpp b/Source/IO/NCColumnFile.cpp index 70c6df513..76249da30 100644 --- a/Source/IO/NCColumnFile.cpp +++ b/Source/IO/NCColumnFile.cpp @@ -153,10 +153,13 @@ ERF::writeToNCColumnFile(const int lev, const int idx_vec = k - kstart; const int ialpha = i - iloc; const int jalpha = j - jloc; - ucol[idx_vec] += velx(i+iloc_shift,j,k) * alpha_u(ialpha, jalpha); - vcol[idx_vec] += vely(i,j+jloc_shift,k) * alpha_v(ialpha, jalpha); - thetacol[idx_vec] += state(i,j,k,RhoTheta_comp) / state(i,j,k,Rho_comp) - * alpha_theta(ialpha, jalpha); + auto tmp = velx(i+iloc_shift,j,k) * alpha_u(ialpha, jalpha); + Gpu::Atomic::Add(&(ucol[idx_vec]), tmp); + tmp = vely(i,j+jloc_shift,k) * alpha_v(ialpha, jalpha); + Gpu::Atomic::Add(&(vcol[idx_vec]), tmp); + tmp = state(i,j,k,Temp_comp) / state(i,j,k,Rho_comp) + * alpha_theta(ialpha, jalpha); + Gpu::Atomic::Add(&(thetacol[idx_vec]), tmp); }); } } diff --git a/Source/Initialization/ERF_init_custom.cpp b/Source/Initialization/ERF_init_custom.cpp index 4a6435935..ce9a3a036 100644 --- a/Source/Initialization/ERF_init_custom.cpp +++ b/Source/Initialization/ERF_init_custom.cpp @@ -99,21 +99,28 @@ ERF::init_custom (int lev) const auto &qc_pert_arr = qc_pert.array(mfi); #endif prob->init_custom_pert(bx, xbx, ybx, zbx, cons_pert_arr, xvel_pert_arr, yvel_pert_arr, zvel_pert_arr, - r_hse_arr, p_hse_arr, z_nd_arr, z_cc_arr, + r_hse_arr, p_hse_arr, z_nd_arr, z_cc_arr, #if defined(ERF_USE_MOISTURE) - qv_pert_arr, qc_pert_arr, qi_pert_arr, + qv_pert_arr, qc_pert_arr, qi_pert_arr, #elif defined(ERF_USE_WARM_NO_PRECIP) - qv_pert_arr, qc_pert_arr, + qv_pert_arr, qc_pert_arr, #endif - geom[lev].data(), mf_m, mf_u, mf_v, - solverChoice); + geom[lev].data(), mf_m, mf_u, mf_v, + solverChoice); } //mfi // Add problem-specific perturbation to background flow MultiFab::Add(lev_new[Vars::cons], cons_pert, Rho_comp, Rho_comp, 1, cons_pert.nGrow()); MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoTheta_comp, RhoTheta_comp, 1, cons_pert.nGrow()); MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoScalar_comp,RhoScalar_comp,1, cons_pert.nGrow()); - MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQKE_comp, RhoQKE_comp, 1, cons_pert.nGrow()); + + // RhoQKE is only relevant if using MYNN2.5 or YSU + if (solverChoice.turbChoice[lev].pbl_type == PBLType::None) { + lev_new[Vars::cons].setVal(0.0,RhoQKE_comp,1); + } else { + MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQKE_comp, RhoQKE_comp, 1, cons_pert.nGrow()); + } + #if defined(ERF_USE_MOISTURE) MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQt_comp, RhoQt_comp, 1, cons_pert.nGrow()); MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQp_comp, RhoQp_comp, 1, cons_pert.nGrow()); diff --git a/Source/Microphysics/Init.cpp b/Source/Microphysics/Init.cpp index 9c2cac5de..a349c59cf 100644 --- a/Source/Microphysics/Init.cpp +++ b/Source/Microphysics/Init.cpp @@ -77,9 +77,6 @@ void Microphysics::Init (const MultiFab& cons_in, MultiFab& qmoist, // output qifall.resize({zlo}, {zhi}); tlatqi.resize({zlo}, {zhi}); - - qpsrc.resize({zlo}, {zhi}); - qpevp.resize({zlo}, {zhi}); } auto accrrc_t = accrrc.table(); @@ -98,8 +95,6 @@ void Microphysics::Init (const MultiFab& cons_in, MultiFab& qmoist, auto rho1d_t = rho1d.table(); auto pres1d_t = pres1d.table(); auto tabs1d_t = tabs1d.table(); - auto qpsrc_t = qpsrc.table(); - auto qpevp_t = qpevp.table(); auto gamaz_t = gamaz.table(); auto zmid_t = zmid.table(); @@ -195,11 +190,6 @@ void Microphysics::Init (const MultiFab& cons_in, MultiFab& qmoist, }); #endif - amrex::ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept { - qpsrc_t(k) = 0.0; - qpevp_t(k) = 0.0; - }); - if(round(gam3) != 2) { std::cout << "cannot compute gamma-function in Microphysics::Init" << std::endl; std::exit(-1); diff --git a/Source/Microphysics/Microphysics.H b/Source/Microphysics/Microphysics.H index d6a5cf660..09e538320 100644 --- a/Source/Microphysics/Microphysics.H +++ b/Source/Microphysics/Microphysics.H @@ -159,8 +159,5 @@ class Microphysics { // data (output) amrex::TableData qifall; amrex::TableData tlatqi; - - amrex::TableData qpsrc; // source of precipitation microphysical processes - amrex::TableData qpevp; // sink of precipitating water due to evaporation }; #endif diff --git a/Source/Microphysics/Precip.cpp b/Source/Microphysics/Precip.cpp index 860a85dd3..61f8daaff 100644 --- a/Source/Microphysics/Precip.cpp +++ b/Source/Microphysics/Precip.cpp @@ -29,8 +29,6 @@ void Microphysics::Precip () { auto evaps2_t = evaps2.table(); auto evapg1_t = evapg1.table(); auto evapg2_t = evapg2.table(); - auto qpsrc_t = qpsrc.table(); - auto qpevp_t = qpevp.table(); auto pres1d_t = pres1d.table(); auto qt = mic_fab_vars[MicVar::qt]; @@ -40,11 +38,6 @@ void Microphysics::Precip () { Real dtn = dt; - ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept { - qpsrc_t(k)=0.0; - qpevp_t(k)=0.0; - }); - // get the temperature, dentisy, theta, qt and qp from input for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { auto tabs_array = mic_fab_vars[MicVar::tabs]->array(mfi); @@ -110,7 +103,6 @@ void Microphysics::Precip () { qt_array(i,j,k) = qt_array(i,j,k) - dq; qp_array(i,j,k) = qp_array(i,j,k) + dq; qn_array(i,j,k) = qn_array(i,j,k) - dq; - amrex::Gpu::Atomic::Add(&qpsrc_t(k), dq); } else if(qp_array(i,j,k) > qp_threshold && qn_array(i,j,k) == 0.0) { @@ -140,11 +132,9 @@ void Microphysics::Precip () { dq = std::max(-0.5*qp_array(i,j,k),dq); qt_array(i,j,k) = qt_array(i,j,k) - dq; qp_array(i,j,k) = qp_array(i,j,k) + dq; - amrex::Gpu::Atomic::Add(&qpevp_t(k), dq); } else { qt_array(i,j,k) = qt_array(i,j,k) + qp_array(i,j,k); - amrex::Gpu::Atomic::Add(&qpevp_t(k), -qp_array(i,j,k)); qp_array(i,j,k) = 0.0; } } diff --git a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp index 86f658973..9fbb9d07d 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_inc.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_inc.cpp @@ -107,6 +107,9 @@ void erf_slow_rhs_inc (int /*level*/, int nrk, { BL_PROFILE_REGION("erf_slow_rhs_pre()"); + DiffChoice dc = solverChoice.diffChoice; + TurbChoice tc = solverChoice.turbChoice[level]; + const MultiFab* t_mean_mf = nullptr; if (most) t_mean_mf = most->get_mac_avg(0,2); @@ -121,15 +124,13 @@ void erf_slow_rhs_inc (int /*level*/, int nrk, AMREX_ALWAYS_ASSERT (!l_use_terrain); const bool l_use_ndiff = solverChoice.use_NumDiff; - const bool l_use_diff = ( (solverChoice.molec_diff_type != MolecDiffType::None) || - (solverChoice.les_type != LESType::None) || - (solverChoice.pbl_type != PBLType::None) ); - const bool cons_visc = ( (solverChoice.molec_diff_type == MolecDiffType::Constant) || - (solverChoice.molec_diff_type == MolecDiffType::ConstantAlpha) ); - const bool l_use_turb = ( solverChoice.les_type == LESType::Smagorinsky || - solverChoice.les_type == LESType::Deardorff || - solverChoice.pbl_type == PBLType::MYNN25 || - solverChoice.pbl_type == PBLType::YSU ); + const bool l_use_diff = ( (dc.molec_diff_type != MolecDiffType::None) || + (tc.les_type != LESType::None) || + (tc.pbl_type != PBLType::None) ); + const bool l_use_turb = ( tc.les_type == LESType::Smagorinsky || + tc.les_type == LESType::Deardorff || + tc.pbl_type == PBLType::MYNN25 || + tc.pbl_type == PBLType::YSU ); const amrex::BCRec* bc_ptr = domain_bcs_type_d.data(); const amrex::BCRec* bc_ptr_h = domain_bcs_type.data();