Skip to content

Commit

Permalink
Merge pull request erf-model#1275 from AMLattanzi/MOST_Terrain
Browse files Browse the repository at this point in the history 
Fix MOST with grid stretching in z.
  • Loading branch information
@ewquon
ewquon authored Oct 26, 2023
2 parents f7c653c + bfcaf08 commit 5d224cf
Show file tree
Hide file tree
Showing 6 changed files with 51 additions and 34 deletions.
5 changes: 4 additions & 1 deletion Source/BoundaryConditions/ABLMost.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,13 +132,16 @@ public:
void
impose_most_bcs (const int lev,
const amrex::Vector<amrex::MultiFab*>& mfs,
amrex::MultiFab* eddyDiffs);
amrex::MultiFab* eddyDiffs,
amrex::MultiFab* z_phys);

template<typename FluxCalc>
void
compute_most_bcs (const int lev,
const amrex::Vector<amrex::MultiFab*>& mfs,
amrex::MultiFab* eddyDiffs,
amrex::MultiFab* z_phys,
const amrex::Real& dz_no_terrain,
const FluxCalc& flux_comp);

void
Expand Down
53 changes: 33 additions & 20 deletions Source/BoundaryConditions/ABLMost.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -108,15 +108,16 @@ ABLMost::compute_fluxes (const int& lev,
void
ABLMost::impose_most_bcs (const int lev,
const Vector<MultiFab*>& mfs,
MultiFab* eddyDiffs)
MultiFab* eddyDiffs,
MultiFab* z_phys)
{
const int zlo = 0;
if (flux_type == FluxCalcType::DONELAN) {
donelan_flux flux_comp(zlo,m_geom[lev].CellSize(2));
compute_most_bcs(lev,mfs,eddyDiffs,flux_comp);
donelan_flux flux_comp(zlo);
compute_most_bcs(lev, mfs, eddyDiffs, z_phys, m_geom[lev].CellSize(2), flux_comp);
} else {
moeng_flux flux_comp(zlo,m_geom[lev].CellSize(2));
compute_most_bcs(lev,mfs,eddyDiffs,flux_comp);
moeng_flux flux_comp(zlo);
compute_most_bcs(lev, mfs, eddyDiffs, z_phys, m_geom[lev].CellSize(2), flux_comp);
}
}

Expand All @@ -134,17 +135,20 @@ void
ABLMost::compute_most_bcs (const int lev,
const Vector<MultiFab*>& mfs,
MultiFab* eddyDiffs,
MultiFab* z_phys,
const Real& dz_no_terrain,
const FluxCalc& flux_comp)
{
const int zlo = 0;
const int icomp = 0;
for (MFIter mfi(*mfs[0]); mfi.isValid(); ++mfi)
{
// Get field arrays
const auto cons_arr = mfs[Vars::cons]->array(mfi);
const auto velx_arr = mfs[Vars::xvel]->array(mfi);
const auto vely_arr = mfs[Vars::yvel]->array(mfi);
const auto eta_arr = eddyDiffs->array(mfi);
const auto cons_arr = mfs[Vars::cons]->array(mfi);
const auto velx_arr = mfs[Vars::xvel]->array(mfi);
const auto vely_arr = mfs[Vars::yvel]->array(mfi);
const auto eta_arr = eddyDiffs->array(mfi);
const auto zphys_arr = (z_phys) ? z_phys->const_array(mfi) : Array4<const Real>{};

// Get average arrays
const auto *const u_mean = m_ma.get_average(lev,0);
Expand All @@ -168,36 +172,45 @@ ABLMost::compute_most_bcs (const int lev,
auto dest_arr = (*mfs[var_idx])[mfi].array();

if (var_idx == Vars::cons) {
amrex::Box b2d = bx; // Copy constructor
Box b2d = bx;
b2d.setBig(2,zlo-1);
int n = Cons::RhoTheta;

ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
flux_comp.compute_t_flux(i,j,k,n,icomp,cons_arr,velx_arr,vely_arr,eta_arr,
umm_arr,tm_arr,u_star_arr,t_star_arr,t_surf_arr,dest_arr);
Real dz = (zphys_arr) ? ( zphys_arr(i,j,zlo) - zphys_arr(i,j,zlo-1) ) : dz_no_terrain;
flux_comp.compute_t_flux(i, j, k, n, icomp, dz,
cons_arr, velx_arr, vely_arr, eta_arr,
umm_arr, tm_arr, u_star_arr, t_star_arr, t_surf_arr,
dest_arr);
});

} else if (var_idx == Vars::xvel || var_idx == Vars::xmom) { //for velx
} else if (var_idx == Vars::xvel || var_idx == Vars::xmom) {

amrex::Box xb2d = surroundingNodes(bx,0); // Copy constructor
Box xb2d = surroundingNodes(bx,0);
xb2d.setBig(2,zlo-1);

ParallelFor(xb2d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
flux_comp.compute_u_flux(i,j,k,icomp,var_idx,cons_arr,velx_arr,vely_arr,eta_arr,
umm_arr,um_arr,u_star_arr,dest_arr);
Real dz = (zphys_arr) ? ( zphys_arr(i,j,zlo) - zphys_arr(i,j,zlo-1) ) : dz_no_terrain;
flux_comp.compute_u_flux(i, j, k, icomp, var_idx, dz,
cons_arr, velx_arr, vely_arr, eta_arr,
umm_arr, um_arr, u_star_arr,
dest_arr);
});

} else if (var_idx == Vars::yvel || var_idx == Vars::ymom) { //for vely
} else if (var_idx == Vars::yvel || var_idx == Vars::ymom) {

amrex::Box yb2d = surroundingNodes(bx,1); // Copy constructor
Box yb2d = surroundingNodes(bx,1);
yb2d.setBig(2,zlo-1);

ParallelFor(yb2d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
flux_comp.compute_v_flux(i,j,k,icomp,var_idx,cons_arr,velx_arr,vely_arr,eta_arr,
umm_arr,vm_arr,u_star_arr,dest_arr);
Real dz = (zphys_arr) ? ( zphys_arr(i,j,zlo) - zphys_arr(i,j,zlo-1) ) : dz_no_terrain;
flux_comp.compute_v_flux(i, j, k, icomp, var_idx, dz,
cons_arr, velx_arr, vely_arr, eta_arr,
umm_arr, vm_arr, u_star_arr,
dest_arr);
});
}
} // var_idx
Expand Down
3 changes: 1 addition & 2 deletions Source/BoundaryConditions/ERF_FillPatch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,7 +158,6 @@ ERF::FillIntermediatePatch (int lev, Real time,
const Vector<MultiFab*>& mfs,
int ng_cons, int ng_vel, bool cons_only,
int icomp_cons, int ncomp_cons,
MultiFab* eddyDiffs,
bool allow_most_bcs)
{
BL_PROFILE_VAR("FillIntermediatePatch()",FillIntermediatePatch);
Expand Down Expand Up @@ -260,7 +259,7 @@ ERF::FillIntermediatePatch (int lev, Real time,

// MOST boundary conditions
if (!(cons_only && ncomp_cons == 1) && m_most && allow_most_bcs)
m_most->impose_most_bcs(lev,mfs,eddyDiffs);
m_most->impose_most_bcs(lev,mfs,eddyDiffs_lev[lev].get(),z_phys_nd[lev].get());
}

//
Expand Down
18 changes: 10 additions & 8 deletions Source/BoundaryConditions/MOSTStress.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -587,9 +587,8 @@ private:
*/
struct moeng_flux
{
moeng_flux (int l_zlo,
amrex::Real l_dz)
: zlo(l_zlo), dz(l_dz) {}
moeng_flux (int l_zlo)
: zlo(l_zlo) {}


AMREX_GPU_DEVICE
Expand All @@ -600,6 +599,7 @@ struct moeng_flux
const int& k,
const int& n,
const int& icomp,
const amrex::Real& dz,
const amrex::Array4<const amrex::Real>& cons_arr,
const amrex::Array4<const amrex::Real>& velx_arr,
const amrex::Array4<const amrex::Real>& vely_arr,
Expand Down Expand Up @@ -663,6 +663,7 @@ struct moeng_flux
const int& k,
const int& icomp,
const int& var_idx,
const amrex::Real& dz,
const amrex::Array4<const amrex::Real>& cons_arr,
const amrex::Array4<const amrex::Real>& velx_arr,
const amrex::Array4<const amrex::Real>& vely_arr,
Expand Down Expand Up @@ -725,6 +726,7 @@ struct moeng_flux
const int& k,
const int& icomp,
const int& var_idx,
const amrex::Real& dz,
const amrex::Array4<const amrex::Real>& cons_arr,
const amrex::Array4<const amrex::Real>& velx_arr,
const amrex::Array4<const amrex::Real>& vely_arr,
Expand Down Expand Up @@ -781,7 +783,6 @@ struct moeng_flux

private:
int zlo;
amrex::Real dz;
};


Expand All @@ -790,9 +791,8 @@ private:
*/
struct donelan_flux
{
donelan_flux (int l_zlo,
amrex::Real l_dz)
: zlo(l_zlo), dz(l_dz) {}
donelan_flux (int l_zlo)
: zlo(l_zlo) {}


AMREX_GPU_DEVICE
Expand All @@ -803,6 +803,7 @@ struct donelan_flux
const int& k,
const int& n,
const int& icomp,
const amrex::Real& dz,
const amrex::Array4<const amrex::Real>& cons_arr,
const amrex::Array4<const amrex::Real>& /*velx_arr*/,
const amrex::Array4<const amrex::Real>& /*vely_arr*/,
Expand Down Expand Up @@ -849,6 +850,7 @@ struct donelan_flux
const int& k,
const int& icomp,
const int& var_idx,
const amrex::Real& dz,
const amrex::Array4<const amrex::Real>& cons_arr,
const amrex::Array4<const amrex::Real>& velx_arr,
const amrex::Array4<const amrex::Real>& vely_arr,
Expand Down Expand Up @@ -915,6 +917,7 @@ struct donelan_flux
const int& k,
const int& icomp,
const int& var_idx,
const amrex::Real& dz,
const amrex::Array4<const amrex::Real>& cons_arr,
const amrex::Array4<const amrex::Real>& velx_arr,
const amrex::Array4<const amrex::Real>& vely_arr,
Expand Down Expand Up @@ -975,6 +978,5 @@ struct donelan_flux

private:
int zlo;
amrex::Real dz;
};
#endif
2 changes: 1 addition & 1 deletion Source/ERF.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -369,7 +369,7 @@ private:
void FillIntermediatePatch (int lev, amrex::Real time,
const amrex::Vector<amrex::MultiFab*>& mfs,
int ng_cons, int ng_vel, bool cons_only, int icomp_cons, int ncomp_cons,
amrex::MultiFab* eddyDiffs, bool allow_most_bcs = true);
bool allow_most_bcs = true);

// Fill all multifabs (and all components) in a vector of multifabs corresponding to the
// grid variables defined in vars_old and vars_new just as FillCoarsePatch.
Expand Down
4 changes: 2 additions & 2 deletions Source/TimeIntegration/TI_utils.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,7 @@

FillIntermediatePatch(level, time_for_fp,
{&S_data[IntVar::cons], &xvel_new, &yvel_new, &zvel_new},
ng_cons_to_use, 0, cons_only, scomp_cons, ncomp_cons, eddyDiffs);
ng_cons_to_use, 0, cons_only, scomp_cons, ncomp_cons);

// Here we don't use include any of the ghost region because we have only updated
// momentum on valid faces
Expand Down Expand Up @@ -132,7 +132,7 @@
FillIntermediatePatch(level, time_for_fp,
{&S_data[IntVar::cons], &xvel_new, &yvel_new, &zvel_new},
ng_cons_to_use, ng_vel, cons_only, scomp_cons, ncomp_cons,
eddyDiffs, allow_most_bcs);
allow_most_bcs);

// Now we can convert back to momentum on valid+ghost since we have
// filled the ghost regions for both velocity and density
Expand Down

0 comments on commit 5d224cf

Please sign in to comment.