diff --git a/Source/ERF.H b/Source/ERF.H index 2af6c99f9..9027c6dc7 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -731,6 +731,10 @@ private: amrex::Vector rW_old; amrex::Vector rW_new; + amrex::Vector xmom_crse_rhs; + amrex::Vector ymom_crse_rhs; + amrex::Vector zmom_crse_rhs; + std::unique_ptr micro; amrex::Vector> qmoist; // (lev,ncomp) This has up to 8 components: qt, qv, qc, qi, qp, qr, qs, qg diff --git a/Source/ERF.cpp b/Source/ERF.cpp index a39415562..c812c9b11 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -197,6 +197,10 @@ 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); + zmom_crse_rhs.resize(nlevs_max); + for (int lev = 0; lev < nlevs_max; ++lev) { vars_new[lev].resize(Vars::NumTypes); vars_old[lev].resize(Vars::NumTypes); @@ -460,8 +464,9 @@ ERF::post_timestep (int nstep, Real time, Real dt_lev0) } } // mfi - // This call refluxes from the lev/lev+1 interface onto lev - getAdvFluxReg(lev+1)->Reflux(vars_new[lev][Vars::cons], 0, 0, ncomp); + // This call refluxes all "slow" cell-centered variables + // (i.e. not density or (rho theta) or velocities) from the lev/lev+1 interface onto lev + getAdvFluxReg(lev+1)->Reflux(vars_new[lev][Vars::cons], 2, 2, ncomp-2); // Here we multiply (rho S) by m^2 after refluxing for (MFIter mfi(vars_new[lev][Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi) { diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index a93d83780..90ab59dff 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/Source/ERF_make_new_arrays.cpp @@ -149,6 +149,15 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, rW_old[lev].define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow_vels); 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); + 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 rU_old[lev].setVal(1.2e21); rV_old[lev].setVal(3.4e22); diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index 282114b7d..3f0fd47cf 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -332,40 +332,39 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp // ***************************************************************************************************** make_physbcs(lev); - // ******************************************************************************************** - // Update the base state at this level by interpolation from coarser level AND copy - // from previous (pre-regrid) base_state array - // ******************************************************************************************** - Interpolater* mapper = &cell_cons_interp; - - Vector fmf = {&base_state[lev ], &base_state[lev ]}; - Vector cmf = {&base_state[lev-1], &base_state[lev-1]}; - Vector ftime = {time, time}; - Vector ctime = {time, time}; - - // Call FillPatch which ASSUMES that all ghost cells at lev-1 have already been filled - FillPatchTwoLevels(temp_base_state, temp_base_state.nGrowVect(), IntVect(0,0,0), - time, cmf, ctime, fmf, ftime, - 0, 0, temp_base_state.nComp(), geom[lev-1], geom[lev], - refRatio(lev-1), mapper, domain_bcs_type, - BCVars::base_bc); - - // Impose bc's outside the domain - (*physbcs_base[lev])(temp_base_state,0,temp_base_state.nComp(),base_state[lev].nGrowVect()); - // ************************************************************************************************* // This will fill the temporary MultiFabs with data from vars_new // NOTE: the momenta here are only used as scratch space, the momenta themselves are not fillpatched - // NOTE: we must create the new base state before calling FillPatch because we will - // interpolate perturbational quantities // ************************************************************************************************* FillPatch(lev, time, {&temp_lev_new[Vars::cons],&temp_lev_new[Vars::xvel], &temp_lev_new[Vars::yvel],&temp_lev_new[Vars::zvel]}, {&temp_lev_new[Vars::cons],&rU_new[lev],&rV_new[lev],&rW_new[lev]}, - base_state[lev], temp_base_state, false); + base_state[lev],temp_base_state,false); - // Now swap the pointers since we needed both old and new in the FillPatch - std::swap(temp_base_state, base_state[lev]); + // ******************************************************************************************** + // Update the base state at this level by interpolation from coarser level AND copy + // from previous (pre-regrid) base_state array + // ******************************************************************************************** + if (lev > 0) { + Interpolater* mapper = &cell_cons_interp; + + Vector fmf = {&base_state[lev ], &base_state[lev ]}; + Vector cmf = {&base_state[lev-1], &base_state[lev-1]}; + Vector ftime = {time, time}; + Vector ctime = {time, time}; + + // Call FillPatch which ASSUMES that all ghost cells at lev-1 have already been filled + FillPatchTwoLevels(temp_base_state, temp_base_state.nGrowVect(), IntVect(0,0,0), + time, cmf, ctime, fmf, ftime, + 0, 0, temp_base_state.nComp(), geom[lev-1], geom[lev], + refRatio(lev-1), mapper, domain_bcs_type, + BCVars::base_bc); + + // Impose bc's outside the domain + (*physbcs_base[lev])(temp_base_state,0,temp_base_state.nComp(),base_state[lev].nGrowVect()); + + std::swap(temp_base_state, base_state[lev]); + } // ******************************************************************************************** // Copy from new into old just in case diff --git a/Source/SourceTerms/ERF_Src_headers.H b/Source/SourceTerms/ERF_Src_headers.H index d8ec184ca..7d19eeb0f 100644 --- a/Source/SourceTerms/ERF_Src_headers.H +++ b/Source/SourceTerms/ERF_Src_headers.H @@ -20,8 +20,7 @@ void make_buoyancy (amrex::Vector< amrex::MultiFab>& S_data, amrex::MultiFab& buoyancy, const amrex::Geometry geom, const SolverChoice& solverChoice, - const amrex::MultiFab* r0, - const amrex::MultiFab* p0, + const amrex::MultiFab& base_state, const int n_qstate, const int anelastic); @@ -56,8 +55,7 @@ void make_mom_sources (int level, int nrk, amrex::MultiFab& xmom_source, amrex::MultiFab& ymom_source, amrex::MultiFab& zmom_source, - amrex::MultiFab* r0, - amrex::MultiFab* p0, + const amrex::MultiFab& base_state, const amrex::Geometry geom, const SolverChoice& solverChoice, std::unique_ptr& mapfac_m, diff --git a/Source/SourceTerms/ERF_buoyancy_utils.H b/Source/SourceTerms/ERF_buoyancy_utils.H index c5d5b75be..061240ae5 100644 --- a/Source/SourceTerms/ERF_buoyancy_utils.H +++ b/Source/SourceTerms/ERF_buoyancy_utils.H @@ -11,8 +11,8 @@ buoyancy_dry_anelastic (int& i, int& j, int& k, amrex::Real const& grav_gpu, - const amrex::Array4& p0_arr, const amrex::Array4& r0_arr, + const amrex::Array4& p0_arr, const amrex::Array4& cell_data) { // Note: this is the same term as the moist anelastic buoyancy when qv = qc = qt = 0 @@ -29,6 +29,32 @@ buoyancy_dry_anelastic (int& i, return (-grav_gpu * (theta_d_wface - theta_d_0_wface) / theta_d_0_wface); } +AMREX_GPU_DEVICE +AMREX_FORCE_INLINE +amrex::Real +buoyancy_dry_anelastic_T (int& i, + int& j, + int& k, + amrex::Real const& grav_gpu, + amrex::Real const& rd_over_cp, + const amrex::Array4& r0_arr, + const amrex::Array4& p0_arr, + const amrex::Array4& cell_data) +{ + amrex::Real rt0_hi = getRhoThetagivenP(p0_arr(i,j,k)); + amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), rt0_hi/r0_arr(i,j,k), rd_over_cp); + amrex::Real t_hi = getTgivenPandTh(p0_arr(i,j,k), cell_data(i,j,k,RhoTheta_comp)/r0_arr(i,j,k), rd_over_cp); + amrex::Real qplus = (t_hi-t0_hi)/t0_hi; + + amrex::Real rt0_lo = getRhoThetagivenP(p0_arr(i,j,k-1)); + amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), rt0_lo/r0_arr(i,j,k-1), rd_over_cp); + amrex::Real t_lo = getTgivenPandTh(p0_arr(i,j,k-1), cell_data(i,j,k-1,RhoTheta_comp)/r0_arr(i,j,k-1), rd_over_cp); + amrex::Real qminus = (t_lo-t0_lo)/t0_lo; + + amrex::Real r0_q_avg = amrex::Real(0.5) * (r0_arr(i,j,k) * qplus + r0_arr(i,j,k-1) * qminus); + return (-r0_q_avg * grav_gpu); +} + AMREX_GPU_DEVICE AMREX_FORCE_INLINE @@ -38,28 +64,26 @@ buoyancy_moist_anelastic (int& i, int& k, amrex::Real const& grav_gpu, amrex::Real const& rv_over_rd, - const amrex::Array4& p0_arr, - const amrex::Array4& r0_arr, + const amrex::Array4& r0_arr, + const amrex::Array4& th0_arr, const amrex::Array4& cell_data) { amrex::Real theta_d_lo = cell_data(i,j,k-1,RhoTheta_comp)/r0_arr(i,j,k-1); - amrex::Real qv_lo = cell_data(i,j,k-1,RhoQ1_comp)/r0_arr(i,j,k-1); - amrex::Real qc_lo = cell_data(i,j,k-1,RhoQ2_comp)/r0_arr(i,j,k-1); + amrex::Real qv_lo = cell_data(i,j,k-1,RhoQ1_comp) /r0_arr(i,j,k-1); + amrex::Real qc_lo = cell_data(i,j,k-1,RhoQ2_comp) /r0_arr(i,j,k-1); amrex::Real qt_lo = qv_lo + qc_lo; amrex::Real theta_v_lo = theta_d_lo * (1.0 - (1.0 - rv_over_rd)*qt_lo - rv_over_rd*qc_lo); amrex::Real theta_d_hi = cell_data(i,j,k,RhoTheta_comp)/r0_arr(i,j,k); - amrex::Real qv_hi = cell_data(i,j,k,RhoQ1_comp)/r0_arr(i,j,k); - amrex::Real qc_hi = cell_data(i,j,k,RhoQ2_comp)/r0_arr(i,j,k); + amrex::Real qv_hi = cell_data(i,j,k,RhoQ1_comp) /r0_arr(i,j,k); + amrex::Real qc_hi = cell_data(i,j,k,RhoQ2_comp) /r0_arr(i,j,k); amrex::Real qt_hi = qv_hi + qc_hi; amrex::Real theta_v_hi = theta_d_hi * (1.0 - (1.0 - rv_over_rd)*qt_hi - rv_over_rd*qc_hi); amrex::Real theta_v_wface = amrex::Real(0.5) * (theta_v_lo + theta_v_hi); - amrex::Real theta_d_0_lo = getRhoThetagivenP(p0_arr(i,j,k-1)) / r0_arr(i,j,k-1); - amrex::Real theta_d_0_hi = getRhoThetagivenP(p0_arr(i,j,k )) / r0_arr(i,j,k ); - amrex::Real theta_v_0_lo = theta_d_0_lo * (1.0 - (1.0 - rv_over_rd)*qt_lo - rv_over_rd*qc_lo); - amrex::Real theta_v_0_hi = theta_d_0_hi * (1.0 - (1.0 - rv_over_rd)*qt_hi - rv_over_rd*qc_hi); + amrex::Real theta_v_0_lo = th0_arr(i,j,k-1) * (1.0 - (1.0 - rv_over_rd)*qt_lo - rv_over_rd*qc_lo); + amrex::Real theta_v_0_hi = th0_arr(i,j,k ) * (1.0 - (1.0 - rv_over_rd)*qt_hi - rv_over_rd*qc_hi); amrex::Real theta_v_0_wface = amrex::Real(0.5) * (theta_v_0_lo + theta_v_0_hi); @@ -74,17 +98,16 @@ buoyancy_dry_default (int& i, int& k, amrex::Real const& grav_gpu, amrex::Real const& rd_over_cp, - const amrex::Array4& p0_arr, const amrex::Array4& r0_arr, + const amrex::Array4& p0_arr, + const amrex::Array4& th0_arr, const amrex::Array4& cell_data) { - amrex::Real rt0_hi = getRhoThetagivenP(p0_arr(i,j,k)); - amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), rt0_hi/r0_arr(i,j,k), rd_over_cp); + amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), th0_arr(i,j,k), rd_over_cp); amrex::Real t_hi = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), cell_data(i,j,k ,RhoTheta_comp)); amrex::Real qplus = (t_hi-t0_hi)/t0_hi; - amrex::Real rt0_lo = getRhoThetagivenP(p0_arr(i,j,k-1)); - amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), rt0_lo/r0_arr(i,j,k-1), rd_over_cp); + amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), th0_arr(i,j,k-1), rd_over_cp); amrex::Real t_lo = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp)); amrex::Real qminus = (t_lo-t0_lo)/t0_lo; diff --git a/Source/SourceTerms/ERF_make_buoyancy.cpp b/Source/SourceTerms/ERF_make_buoyancy.cpp index 528c82edf..ca101605e 100644 --- a/Source/SourceTerms/ERF_make_buoyancy.cpp +++ b/Source/SourceTerms/ERF_make_buoyancy.cpp @@ -34,8 +34,7 @@ void make_buoyancy (Vector& S_data, MultiFab& buoyancy, const amrex::Geometry geom, const SolverChoice& solverChoice, - const MultiFab* r0, - const MultiFab* p0, + const MultiFab& base_state, const int n_qstate, const int anelastic) { @@ -50,6 +49,10 @@ void make_buoyancy (Vector& S_data, Real rd_over_cp = solverChoice.rdOcp; Real rv_over_rd = R_v/R_d; + MultiFab r0 (base_state, make_alias, BaseState::r0_comp , 1); + MultiFab p0 (base_state, make_alias, BaseState::p0_comp , 1); + MultiFab th0(base_state, make_alias, BaseState::th0_comp, 1); + if (anelastic == 1) { #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) @@ -66,15 +69,16 @@ void make_buoyancy (Vector& S_data, const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); // Base state density and pressure - const Array4& r0_arr = r0->const_array(mfi); - const Array4& p0_arr = p0->const_array(mfi); + const Array4& r0_arr = r0.const_array(mfi); + const Array4& p0_arr = p0.const_array(mfi); + const Array4& th0_arr = th0.const_array(mfi); if (solverChoice.moisture_type == MoistureType::None) { ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { buoyancy_fab(i, j, k) = buoyancy_dry_anelastic(i,j,k, grav_gpu[2], - p0_arr,r0_arr,cell_data); + r0_arr,p0_arr,cell_data); }); } else { // NOTE: For decomposition in the vertical direction, klo may not @@ -84,7 +88,7 @@ void make_buoyancy (Vector& S_data, { buoyancy_fab(i, j, k) = buoyancy_moist_anelastic(i,j,k, grav_gpu[2],rv_over_rd, - p0_arr,r0_arr,cell_data); + r0_arr,th0_arr,cell_data); }); } } // mfi @@ -113,7 +117,7 @@ void make_buoyancy (Vector& S_data, const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); // Base state density - const Array4& r0_arr = r0->const_array(mfi); + const Array4& r0_arr = r0.const_array(mfi); ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { @@ -142,8 +146,9 @@ void make_buoyancy (Vector& S_data, if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1); // Base state density and pressure - const Array4& r0_arr = r0->const_array(mfi); - const Array4& p0_arr = p0->const_array(mfi); + const Array4& r0_arr = r0.const_array(mfi); + const Array4& p0_arr = p0.const_array(mfi); + const Array4& th0_arr = th0.const_array(mfi); const Array4 & cell_data = S_data[IntVars::cons].array(mfi); const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); @@ -152,7 +157,7 @@ void make_buoyancy (Vector& S_data, { buoyancy_fab(i, j, k) = buoyancy_dry_default(i,j,k, grav_gpu[2],rd_over_cp, - p0_arr,r0_arr,cell_data); + r0_arr,p0_arr,th0_arr,cell_data); }); } // mfi } // buoyancy_type @@ -187,7 +192,7 @@ void make_buoyancy (Vector& S_data, const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); // Base state density - const Array4& r0_arr = r0->const_array(mfi); + const Array4& r0_arr = r0.const_array(mfi); ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { diff --git a/Source/SourceTerms/ERF_make_mom_sources.cpp b/Source/SourceTerms/ERF_make_mom_sources.cpp index b60190bfe..2e9a4d31f 100644 --- a/Source/SourceTerms/ERF_make_mom_sources.cpp +++ b/Source/SourceTerms/ERF_make_mom_sources.cpp @@ -48,7 +48,7 @@ void make_mom_sources (int level, MultiFab & xmom_src, MultiFab & ymom_src, MultiFab & zmom_src, - MultiFab* r0, MultiFab* p0, + const MultiFab& base_state, const Geometry geom, const SolverChoice& solverChoice, std::unique_ptr& mapfac_m, @@ -190,7 +190,7 @@ void make_mom_sources (int level, // ***************************************************************************** // 1. Create the BUOYANCY forcing term in the z-direction // ***************************************************************************** - make_buoyancy(S_data, S_prim, zmom_src, geom, solverChoice, r0, p0, + make_buoyancy(S_data, S_prim, zmom_src, geom, solverChoice, base_state, n_qstate, solverChoice.anelastic[level]); // ***************************************************************************** diff --git a/Source/TimeIntegration/ERF_Advance.cpp b/Source/TimeIntegration/ERF_Advance.cpp index 2a8a73497..f634580f0 100644 --- a/Source/TimeIntegration/ERF_Advance.cpp +++ b/Source/TimeIntegration/ERF_Advance.cpp @@ -241,6 +241,36 @@ 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}); + zmom_crse_rhs[lev+1].mult(1.0/dt_lev,0,1,0); } } diff --git a/Source/TimeIntegration/ERF_TI_slow_headers.H b/Source/TimeIntegration/ERF_TI_slow_headers.H index b4959ac65..86dfb38b5 100644 --- a/Source/TimeIntegration/ERF_TI_slow_headers.H +++ b/Source/TimeIntegration/ERF_TI_slow_headers.H @@ -67,6 +67,9 @@ void erf_slow_rhs_pre (int level, int finest_level, int nrk, const amrex::MultiFab& xmom_src, const amrex::MultiFab& ymom_src, const amrex::MultiFab& zmom_src, + const amrex::MultiFab* xmom_crse_rhs, + const amrex::MultiFab* ymom_crse_rhs, + const amrex::MultiFab* zmom_crse_rhs, amrex::MultiFab* Tau11, amrex::MultiFab* Tau22, amrex::MultiFab* Tau33, diff --git a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H index 73976894a..092f87589 100644 --- a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H @@ -102,18 +102,21 @@ } // mfi - MultiFab r_hse_new (base_state_new[level], make_alias, BaseState::r0_comp, 1); - MultiFab p_hse_new (base_state_new[level], make_alias, BaseState::p0_comp, 1); - MultiFab pi_hse_new(base_state_new[level], make_alias, BaseState::pi0_comp, 1); + MultiFab r_hse_new (base_state_new[level], make_alias, BaseState::r0_comp, 1); + MultiFab p_hse_new (base_state_new[level], make_alias, BaseState::p0_comp, 1); + MultiFab pi_hse_new (base_state_new[level], make_alias, BaseState::pi0_comp, 1); + MultiFab th_hse_new (base_state_new[level], make_alias, BaseState::th0_comp, 1); - MultiFab* r0_new = &r_hse_new; - MultiFab* p0_new = &p_hse_new; + MultiFab* r0_new = &r_hse_new; + MultiFab* p0_new = &p_hse_new; + MultiFab* pi0_new = &pi_hse_new; + MultiFab* th0_new = &th_hse_new; make_mom_sources(level, nrk, slow_dt, old_stage_time, S_data, S_prim, z_phys_nd[level], z_phys_cc[level], xvel_new, yvel_new, xmom_src, ymom_src, zmom_src, - r0_new, p0_new, fine_geom, solverChoice, + base_state_new[level], fine_geom, solverChoice, mapfac_m[level], mapfac_u[level], mapfac_v[level], dptr_u_geos, dptr_v_geos, dptr_wbar_sub, d_rayleigh_ptrs_at_lev, d_sponge_ptrs_at_lev, @@ -122,6 +125,9 @@ erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_old, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, z_t_rk[level], Omega, cc_src, xmom_src, ymom_src, zmom_src, + (level > 0) ? &xmom_crse_rhs[level] : nullptr, + (level > 0) ? &ymom_crse_rhs[level] : nullptr, + (level > 0) ? &zmom_crse_rhs[level] : nullptr, Tau11_lev[level].get(), Tau22_lev[level].get(), Tau33_lev[level].get(), Tau12_lev[level].get(), Tau13_lev[level].get(), Tau21_lev[level].get(), Tau23_lev[level].get(), Tau31_lev[level].get(), Tau32_lev[level].get(), SmnSmn, eddyDiffs, Hfx1, Hfx2, Hfx3, Q1fx1, Q1fx2, Q1fx3, Q2fx3, Diss, @@ -156,6 +162,8 @@ const GpuArray dxInv = fine_geom.InvCellSizeArray(); + const Real l_rdOcp = solverChoice.rdOcp; + #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif @@ -168,8 +176,10 @@ const Array4 r0_new_arr = r0_new->array(mfi); const Array4 r0_tmp_arr = r0_temp.array(mfi); - const Array4 p0_arr = p0->const_array(mfi); - const Array4 p0_new_arr = p0_new->array(mfi); + const Array4 p0_arr = p0->const_array(mfi); + const Array4 p0_new_arr = p0_new->array(mfi); + const Array4 pi0_new_arr = pi0_new->array(mfi); + const Array4 th0_new_arr = th0_new->array(mfi); const Array4& z_t_arr = z_t_rk[level]->array(mfi); @@ -205,7 +215,9 @@ r0_new_arr(i,j,k) = rho0_new / dJ_new_arr(i,j,k); rt0_tmp_new /= dJ_new_arr(i,j,k); - p0_new_arr(i,j,k) = getPgivenRTh(rt0_tmp_new); + p0_new_arr(i,j,k) = getPgivenRTh(rt0_tmp_new); + pi0_new_arr(i,j,k) = getExnergivenRTh(rt0_tmp_new, l_rdOcp); + th0_new_arr(i,j,k) = rt0_tmp_new / r0_new_arr(i,j,k); }); } // MFIter r0_new->FillBoundary(fine_geom.periodicity()); @@ -217,7 +229,7 @@ z_phys_nd[level], z_phys_cc[level], xvel_new, yvel_new, xmom_src, ymom_src, zmom_src, - r0, p0, fine_geom, solverChoice, + base_state[level], fine_geom, solverChoice, mapfac_m[level], mapfac_u[level], mapfac_v[level], dptr_u_geos, dptr_v_geos, dptr_wbar_sub, d_rayleigh_ptrs_at_lev, d_sponge_ptrs_at_lev, @@ -226,6 +238,9 @@ erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_old, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, z_t_rk[level], Omega, cc_src, xmom_src, ymom_src, zmom_src, + (level > 0) ? &xmom_crse_rhs[level] : nullptr, + (level > 0) ? &ymom_crse_rhs[level] : nullptr, + (level > 0) ? &zmom_crse_rhs[level] : nullptr, Tau11_lev[level].get(), Tau22_lev[level].get(), Tau33_lev[level].get(), Tau12_lev[level].get(), Tau13_lev[level].get(), Tau21_lev[level].get(), Tau23_lev[level].get(), Tau31_lev[level].get(), Tau32_lev[level].get(), SmnSmn, eddyDiffs, Hfx1, Hfx2, Hfx3, Q1fx1, Q1fx2, Q1fx3,Q2fx3, Diss, @@ -424,7 +439,7 @@ z_phys_nd[level], z_phys_cc[level], xvel_new, yvel_new, xmom_src, ymom_src, zmom_src, - r0, p0, fine_geom, solverChoice, + base_state[level], fine_geom, solverChoice, mapfac_m[level], mapfac_u[level], mapfac_v[level], dptr_u_geos, dptr_v_geos, dptr_wbar_sub, d_rayleigh_ptrs_at_lev, d_sponge_ptrs_at_lev, @@ -434,6 +449,9 @@ S_rhs, S_old, S_data, S_prim, S_scratch, xvel_new, yvel_new, zvel_new, z_t_rk[level], Omega, cc_src, xmom_src, ymom_src, zmom_src, + (level > 0) ? &xmom_crse_rhs[level] : nullptr, + (level > 0) ? &ymom_crse_rhs[level] : nullptr, + (level > 0) ? &zmom_crse_rhs[level] : nullptr, Tau11_lev[level].get(), Tau22_lev[level].get(), Tau33_lev[level].get(), Tau12_lev[level].get(), Tau13_lev[level].get(), Tau21_lev[level].get(), Tau23_lev[level].get(), Tau31_lev[level].get(), Tau32_lev[level].get(), SmnSmn, eddyDiffs, Hfx1, Hfx2, Hfx3, Q1fx1, Q1fx2, Q1fx3, Q2fx3, Diss, diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index 51890a8cd..c92ed95cb 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -81,6 +81,9 @@ void erf_slow_rhs_pre (int level, int finest_level, const MultiFab& xmom_src, const MultiFab& ymom_src, const MultiFab& zmom_src, + const MultiFab* xmom_crse_rhs, + const MultiFab* ymom_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,