Skip to content

Commit

Permalink
Tweaks (#1264)
Browse files Browse the repository at this point in the history 
* update particle functionality

* some minor cleanup
  • Loading branch information
@asalmgren
asalmgren authored Oct 11, 2023
1 parent 0db2007 commit 522382e
Show file tree
Hide file tree
Showing 8 changed files with 69 additions and 84 deletions.
42 changes: 25 additions & 17 deletions Source/ERF.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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") {
Expand Down
39 changes: 14 additions & 25 deletions Source/IO/ERF_WriteScalarProfiles.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
Expand Down
11 changes: 7 additions & 4 deletions Source/IO/NCColumnFile.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
});
}
}
Expand Down
19 changes: 13 additions & 6 deletions Source/Initialization/ERF_init_custom.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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());
Expand Down
10 changes: 0 additions & 10 deletions Source/Microphysics/Init.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand All @@ -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();
Expand Down Expand Up @@ -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);
Expand Down
3 changes: 0 additions & 3 deletions Source/Microphysics/Microphysics.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -159,8 +159,5 @@ class Microphysics {
// data (output)
amrex::TableData<amrex::Real, 1> qifall;
amrex::TableData<amrex::Real, 1> tlatqi;

amrex::TableData<amrex::Real, 1> qpsrc; // source of precipitation microphysical processes
amrex::TableData<amrex::Real, 1> qpevp; // sink of precipitating water due to evaporation
};
#endif
10 changes: 0 additions & 10 deletions Source/Microphysics/Precip.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Expand All @@ -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);
Expand Down Expand Up @@ -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) {

Expand Down Expand Up @@ -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;
}
}
Expand Down
19 changes: 10 additions & 9 deletions Source/TimeIntegration/ERF_slow_rhs_inc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand All @@ -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();
Expand Down

0 comments on commit 522382e

Please sign in to comment.