From 72628167ff53cc144a2faf1374979469d1436034 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Thu, 31 Oct 2024 11:35:30 -0700 Subject: [PATCH] fix issue at coarse/fine boundaries (#1920) * fix issue at coarse/fine boundaries * remove unused --- Source/ERF.H | 4 +- Source/ERF.cpp | 4 +- Source/ERF_make_new_arrays.cpp | 7 +-- Source/TimeIntegration/ERF_Advance.cpp | 51 ++++++++++---------- Source/TimeIntegration/ERF_TI_fast_headers.H | 1 - Source/TimeIntegration/ERF_TI_fast_rhs_fun.H | 2 - Source/TimeIntegration/ERF_fast_rhs_N.cpp | 39 ++------------- Source/TimeIntegration/ERF_slow_rhs_pre.cpp | 26 +++++----- 8 files changed, 48 insertions(+), 86 deletions(-) diff --git a/Source/ERF.H b/Source/ERF.H index 9027c6dc7..a914b926f 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -731,8 +731,8 @@ private: amrex::Vector rW_old; amrex::Vector rW_new; - amrex::Vector xmom_crse_rhs; - amrex::Vector ymom_crse_rhs; + // amrex::Vector xmom_crse_rhs; + // amrex::Vector ymom_crse_rhs; amrex::Vector zmom_crse_rhs; std::unique_ptr micro; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index c812c9b11..65e216654 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -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) { diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index 90ab59dff..d1a856c19 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/Source/ERF_make_new_arrays.cpp @@ -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 diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index f634580f0..197edaddd 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -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; + + // 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 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); - } } // *********************************************************************************************** diff --git a/Source/TimeIntegration/ERF_TI_fast_headers.H b/Source/TimeIntegration/ERF_TI_fast_headers.H index 8a244dbcb..5c470d9da 100644 --- a/Source/TimeIntegration/ERF_TI_fast_headers.H +++ b/Source/TimeIntegration/ERF_TI_fast_headers.H @@ -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& S_data, amrex::Vector& S_scratch, const amrex::Geometry geom, diff --git a/Source/TimeIntegration/ERF_TI_fast_rhs_fun.H b/Source/TimeIntegration/ERF_TI_fast_rhs_fun.H index 152c60227..1213b8cf6 100644 --- a/Source/TimeIntegration/ERF_TI_fast_rhs_fun.H +++ b/Source/TimeIntegration/ERF_TI_fast_rhs_fun.H @@ -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], @@ -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], diff --git a/Source/TimeIntegration/ERF_fast_rhs_N.cpp b/Source/TimeIntegration/ERF_fast_rhs_N.cpp index b11e62254..a77fda150 100644 --- a/Source/TimeIntegration/ERF_fast_rhs_N.cpp +++ b/Source/TimeIntegration/ERF_fast_rhs_N.cpp @@ -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& S_data, // S_sum = most recent full solution Vector& S_scratch, // S_sum_old at most recent fast timestep for (rho theta) const Geometry geom, @@ -281,8 +280,6 @@ void erf_fast_rhs_N (int step, int nrk, const Array4& slow_rhs_cons = S_slow_rhs[IntVars::cons].const_array(mfi); const Array4& slow_rhs_rho_w = S_slow_rhs[IntVars::zmom].const_array(mfi); - const Array4& rho_w_crse_arr = (level > 0) ? zmom_crse_rhs->const_array(mfi) : Array4{}; - const Array4& prev_zmom = S_prev[IntVars::zmom].const_array(mfi); const Array4< Real>& cur_zmom = S_data[IntVars::zmom].array(mfi); @@ -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 diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index b35f3329e..1a42951f0 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -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, @@ -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 dxInv = geom.InvCellSizeArray(); @@ -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); @@ -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");