Skip to content

Commit

Permalink
fix issue at coarse/fine boundaries
Browse files Browse the repository at this point in the history
  • Loading branch information
@asalmgren
asalmgren committed Oct 31, 2024
1 parent f6b1c87 commit 7947de3
Show file tree
Hide file tree
Showing 8 changed files with 47 additions and 85 deletions.
4 changes: 2 additions & 2 deletions Source/ERF.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -731,8 +731,8 @@ private:
amrex::Vector<amrex::MultiFab> rW_old;
amrex::Vector<amrex::MultiFab> rW_new;

amrex::Vector<amrex::MultiFab> xmom_crse_rhs;
amrex::Vector<amrex::MultiFab> ymom_crse_rhs;
// amrex::Vector<amrex::MultiFab> xmom_crse_rhs;
// amrex::Vector<amrex::MultiFab> ymom_crse_rhs;
amrex::Vector<amrex::MultiFab> zmom_crse_rhs;

std::unique_ptr<Microphysics> micro;
Expand Down
4 changes: 2 additions & 2 deletions Source/ERF.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -197,8 +197,8 @@ ERF::ERF_shared ()
rV_old.resize(nlevs_max);
rW_old.resize(nlevs_max);

xmom_crse_rhs.resize(nlevs_max);
ymom_crse_rhs.resize(nlevs_max);
// xmom_crse_rhs.resize(nlevs_max);
// ymom_crse_rhs.resize(nlevs_max);
zmom_crse_rhs.resize(nlevs_max);

for (int lev = 0; lev < nlevs_max; ++lev) {
Expand Down
7 changes: 2 additions & 5 deletions Source/ERF_make_new_arrays.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,12 +150,9 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
rW_new[lev].define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow_vels);

if (lev > 0) {
xmom_crse_rhs[lev].define(convert(ba, IntVect(1,0,0)), dm, 1, IntVect{0});
xmom_crse_rhs[lev].setVal(3.456e22);
ymom_crse_rhs[lev].define(convert(ba, IntVect(0,1,0)), dm, 1, IntVect{0});
ymom_crse_rhs[lev].setVal(3.456e22);
//xmom_crse_rhs[lev].define(convert(ba, IntVect(1,0,0)), dm, 1, IntVect{0});
//ymom_crse_rhs[lev].define(convert(ba, IntVect(0,1,0)), dm, 1, IntVect{0});
zmom_crse_rhs[lev].define(convert(ba, IntVect(0,0,1)), dm, 1, IntVect{0});
zmom_crse_rhs[lev].setVal(3.456e22);
}

// We do this here just so they won't be undefined in the initial FillPatch
Expand Down
49 changes: 25 additions & 24 deletions Source/TimeIntegration/ERF_Advance.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -241,37 +241,38 @@ ERF::Advance (int lev, Real time, Real dt_lev, int iteration, int /*ncycle*/)
state_new[IntVars::zmom].FillBoundary(geom[lev].periodicity());
FPr_w[lev].RegisterCoarseData({&state_old[IntVars::zmom], &state_new[IntVars::zmom]},
{time, time + dt_lev});
}

//
// Now create a MultiFab that holds (S_new - S_old) / dt from the coarse level interpolated
// on to the coarse/fine boundary at the fine resolution
//
PhysBCFunctNoOp null_bc;

MultiFab tempx(vars_new[lev+1][Vars::xvel].boxArray(),vars_new[lev+1][Vars::xvel].DistributionMap(),1,0);
tempx.setVal(0.0);
xmom_crse_rhs[lev+1].setVal(0.0);
FPr_u[lev].FillSet(tempx , time , null_bc, domain_bcs_type);
FPr_u[lev].FillSet(xmom_crse_rhs[lev+1], time+dt_lev, null_bc, domain_bcs_type);
MultiFab::Subtract(xmom_crse_rhs[lev+1],tempx,0,0,1,IntVect{0});
xmom_crse_rhs[lev+1].mult(1.0/dt_lev,0,1,0);

MultiFab tempy(vars_new[lev+1][Vars::yvel].boxArray(),vars_new[lev+1][Vars::yvel].DistributionMap(),1,0);
tempy.setVal(0.0);
ymom_crse_rhs[lev+1].setVal(0.0);
FPr_v[lev].FillSet(tempy , time , null_bc, domain_bcs_type);
FPr_v[lev].FillSet(ymom_crse_rhs[lev+1], time+dt_lev, null_bc, domain_bcs_type);
MultiFab::Subtract(ymom_crse_rhs[lev+1],tempy,0,0,1,IntVect{0});
ymom_crse_rhs[lev+1].mult(1.0/dt_lev,0,1,0);

MultiFab tempz(vars_new[lev+1][Vars::zvel].boxArray(),vars_new[lev+1][Vars::zvel].DistributionMap(),1,0);
tempz.setVal(0.0);
zmom_crse_rhs[lev+1].setVal(0.0);
FPr_w[lev].FillSet(tempz , time , null_bc, domain_bcs_type);
FPr_w[lev].FillSet(zmom_crse_rhs[lev+1], time+dt_lev, null_bc, domain_bcs_type);
MultiFab::Subtract(zmom_crse_rhs[lev+1],tempz,0,0,1,IntVect{0});
Interpolater* mapper_f = &face_cons_linear_interp;

// MultiFab tempx(vars_new[lev+1][Vars::xvel].boxArray(),vars_new[lev+1][Vars::xvel].DistributionMap(),1,0);
// tempx.setVal(0.0);
// xmom_crse_rhs[lev+1].setVal(0.0);
// FPr_u[lev].FillSet(tempx , time , null_bc, domain_bcs_type);
// FPr_u[lev].FillSet(xmom_crse_rhs[lev+1], time+dt_lev, null_bc, domain_bcs_type);
// MultiFab::Subtract(xmom_crse_rhs[lev+1],tempx,0,0,1,IntVect{0});
// xmom_crse_rhs[lev+1].mult(1.0/dt_lev,0,1,0);

// MultiFab tempy(vars_new[lev+1][Vars::yvel].boxArray(),vars_new[lev+1][Vars::yvel].DistributionMap(),1,0);
// tempy.setVal(0.0);
// ymom_crse_rhs[lev+1].setVal(0.0);
// FPr_v[lev].FillSet(tempy , time , null_bc, domain_bcs_type);
// FPr_v[lev].FillSet(ymom_crse_rhs[lev+1], time+dt_lev, null_bc, domain_bcs_type);
// MultiFab::Subtract(ymom_crse_rhs[lev+1],tempy,0,0,1,IntVect{0});
// ymom_crse_rhs[lev+1].mult(1.0/dt_lev,0,1,0);

MultiFab temp_state(zmom_crse_rhs[lev+1].boxArray(),zmom_crse_rhs[lev+1].DistributionMap(),1,0);
InterpFromCoarseLevel(temp_state, IntVect{0}, IntVect{0}, state_old[IntVars::zmom], 0, 0, 1,
geom[lev], geom[lev+1], refRatio(lev), mapper_f, domain_bcs_type, BCVars::zvel_bc);
InterpFromCoarseLevel(zmom_crse_rhs[lev+1], IntVect{0}, IntVect{0}, state_new[IntVars::zmom], 0, 0, 1,
geom[lev], geom[lev+1], refRatio(lev), mapper_f, domain_bcs_type, BCVars::zvel_bc);
MultiFab::Subtract(zmom_crse_rhs[lev+1],temp_state,0,0,1,IntVect{0});
zmom_crse_rhs[lev+1].mult(1.0/dt_lev,0,1,0);
}
}

// ***********************************************************************************************
Expand Down
1 change: 0 additions & 1 deletion Source/TimeIntegration/ERF_TI_fast_headers.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,6 @@ void erf_fast_rhs_N (int step, int nrk, int level, int finest_level,
const amrex::MultiFab& S_stage_prim,
const amrex::MultiFab& pi_stage,
const amrex::MultiFab& fast_coeffs,
const amrex::MultiFab* zmom_crse_rhs,
amrex::Vector<amrex::MultiFab >& S_data,
amrex::Vector<amrex::MultiFab >& S_scratch,
const amrex::Geometry geom,
Expand Down
2 changes: 0 additions & 2 deletions Source/TimeIntegration/ERF_TI_fast_rhs_fun.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -157,7 +157,6 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk,
// and S_data is the new-time solution to be defined here
erf_fast_rhs_N(fast_step, nrk, level, finest_level,
S_slow_rhs, S_old, S_stage, S_prim, pi_stage, fast_coeffs,
(level > 0) ? &zmom_crse_rhs[level] : nullptr,
S_data, S_scratch, fine_geom, solverChoice.gravity,
dtau, beta_s, inv_fac,
mapfac_m[level], mapfac_u[level], mapfac_v[level],
Expand All @@ -167,7 +166,6 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk,
// and as the new-time solution to be defined here
erf_fast_rhs_N(fast_step, nrk, level, finest_level,
S_slow_rhs, S_data, S_stage, S_prim, pi_stage, fast_coeffs,
(level > 0) ? &zmom_crse_rhs[level] : nullptr,
S_data, S_scratch, fine_geom, solverChoice.gravity,
dtau, beta_s, inv_fac,
mapfac_m[level], mapfac_u[level], mapfac_v[level],
Expand Down
39 changes: 4 additions & 35 deletions Source/TimeIntegration/ERF_fast_rhs_N.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,6 @@ void erf_fast_rhs_N (int step, int nrk,
const MultiFab & S_stage_prim, // Primitive version of S_stage_data[IntVars::cons]
const MultiFab & pi_stage, // Exner function evaluated at last stage
const MultiFab &fast_coeffs, // Coeffs for tridiagonal solve
const MultiFab* zmom_crse_rhs, // Source term from coarser level; non-zero on c/f boundary only
Vector<MultiFab>& S_data, // S_sum = most recent full solution
Vector<MultiFab>& S_scratch, // S_sum_old at most recent fast timestep for (rho theta)
const Geometry geom,
Expand Down Expand Up @@ -281,8 +280,6 @@ void erf_fast_rhs_N (int step, int nrk,
const Array4<const Real>& slow_rhs_cons = S_slow_rhs[IntVars::cons].const_array(mfi);
const Array4<const Real>& slow_rhs_rho_w = S_slow_rhs[IntVars::zmom].const_array(mfi);

const Array4<const Real>& rho_w_crse_arr = (level > 0) ? zmom_crse_rhs->const_array(mfi) : Array4<const Real>{};

const Array4<const Real>& prev_zmom = S_prev[IntVars::zmom].const_array(mfi);
const Array4< Real>& cur_zmom = S_data[IntVars::zmom].array(mfi);

Expand Down Expand Up @@ -424,40 +421,12 @@ void erf_fast_rhs_N (int step, int nrk,
ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int)
{
// w at bottom boundary of grid is 0 if at domain boundary, otherwise w = w_old + dtau * slow_rhs
RHS_a (i,j,lo.z ) = prev_zmom(i,j,lo.z ) - stage_zmom(i,j,lo.z);

#if 1
int k = lo.z;
if (level > 0) {
if (k != 0) {
RHS_a (i,j,lo.z ) += dtau * rho_w_crse_arr(i,j,lo.z);
} else {
RHS_a (i,j,lo.z ) += dtau * slow_rhs_rho_w(i,j,lo.z);
}
} else {
RHS_a (i,j,lo.z ) += dtau * slow_rhs_rho_w(i,j,lo.z);
}
#else
RHS_a (i,j,lo.z ) += dtau * slow_rhs_rho_w(i,j,lo.z);
#endif
RHS_a (i,j,lo.z) = prev_zmom(i,j,lo.z) - stage_zmom(i,j,lo.z)
+ dtau * slow_rhs_rho_w(i,j,lo.z);

// w at top boundary of grid is 0 if at domain boundary, otherwise w = w_old + dtau * slow_rhs
RHS_a (i,j,hi.z+1) = prev_zmom(i,j,hi.z+1) - stage_zmom(i,j,hi.z+1);

#if 1
k = hi.z+1;
if (level > 0) {
if (std::abs(rho_w_crse_arr(i,j,k)) > 0.) {
RHS_a (i,j,hi.z+1) += dtau * rho_w_crse_arr(i,j,hi.z+1);
} else {
RHS_a (i,j,hi.z+1) += dtau * slow_rhs_rho_w(i,j,hi.z+1);
}
} else {
RHS_a (i,j,hi.z+1) += dtau * slow_rhs_rho_w(i,j,hi.z+1);
}
#else
RHS_a (i,j,hi.z+1) += dtau * slow_rhs_rho_w(i,j,hi.z+1);
#endif
RHS_a (i,j,hi.z+1) = prev_zmom(i,j,hi.z+1) - stage_zmom(i,j,hi.z+1)
+ dtau * slow_rhs_rho_w(i,j,hi.z+1);
}); // b2d

#ifdef AMREX_USE_GPU
Expand Down
26 changes: 12 additions & 14 deletions Source/TimeIntegration/ERF_slow_rhs_pre.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ void erf_slow_rhs_pre (int level, int finest_level,
const MultiFab& xmom_src,
const MultiFab& ymom_src,
const MultiFab& zmom_src,
const MultiFab* /*zmom_crse_rhs*/,
const MultiFab* zmom_crse_rhs,
MultiFab* Tau11, MultiFab* Tau22, MultiFab* Tau33,
MultiFab* Tau12, MultiFab* Tau13, MultiFab* Tau21,
MultiFab* Tau23, MultiFab* Tau31, MultiFab* Tau32,
Expand Down Expand Up @@ -169,6 +169,7 @@ void erf_slow_rhs_pre (int level, int finest_level,
AMREX_ALWAYS_ASSERT(!l_use_terrain || !l_anelastic);

const Box& domain = geom.Domain();
const int domlo_z = domain.smallEnd(2);
const int domhi_z = domain.bigEnd(2);

const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();
Expand Down Expand Up @@ -745,7 +746,6 @@ void erf_slow_rhs_pre (int level, int finest_level,
}
});

#if 0
auto const lo = lbound(bx);
auto const hi = ubound(bx);

Expand All @@ -756,22 +756,20 @@ void erf_slow_rhs_pre (int level, int finest_level,

Box b2d = bx; b2d.setRange(2,0);

// Note: we don't need to test on being at the domain boundary
// because tbz is shrunk to not hit top and bottom domain boundaries

ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int ) { // bottom of box but not of domain
if (std::abs(rho_w_rhs_crse(i,j,lo.z)) > 0.) {
if (lo.z > domlo_z) {
ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int ) // bottom of box but not of domain
{
rho_w_rhs(i,j,lo.z) = rho_w_rhs_crse(i,j,lo.z);
}
});
});
}

ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int ) { // top of box but not of domain
if (std::abs(rho_w_rhs_crse(i,j,hi.z+1)) > 0.) {
if (hi.z < domhi_z+1) {
ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int ) // top of box but not of domain
{
rho_w_rhs(i,j,hi.z+1) = rho_w_rhs_crse(i,j,hi.z+1);
}
});
});
}
}
#endif

{
BL_PROFILE("slow_rhs_pre_fluxreg");
Expand Down

0 comments on commit 7947de3

Please sign in to comment.