diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 2d056fbbd..573d8c7bf 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -494,7 +494,14 @@ ERF::post_timestep (int nstep, Real time, Real dt_lev0) // We need to do this before anything else because refluxing changes the // values of coarse cells underneath fine grids with the assumption they'll // be over-written by averaging down - AverageDownTo(lev,0,ncomp); + int src_comp; + if (solverChoice.anelastic[lev]) { + src_comp = 1; + } else { + src_comp = 0; + } + int num_comp = ncomp - src_comp; + AverageDownTo(lev,src_comp,num_comp); } } diff --git a/Source/IO/ERF_Checkpoint.cpp b/Source/IO/ERF_Checkpoint.cpp index 4a27fe16d..72c3a9da1 100644 --- a/Source/IO/ERF_Checkpoint.cpp +++ b/Source/IO/ERF_Checkpoint.cpp @@ -425,13 +425,16 @@ ERF::ReadCheckpointFile () } MultiFab base(grids[lev],dmap[lev],ncomp_base,ng_base); VisMF::Read(base, MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "BaseState")); + MultiFab::Copy(base_state[lev],base,0,0,ncomp_base,ng_base); + if (read_old_base_state) { for (MFIter mfi(base_state[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) { - const Box& bx = mfi.growntilebox(1); + // We only compute theta_0 on valid cells since we will impose domain BC's after restart + const Box& bx = mfi.tilebox(); Array4 const& fab = base_state[lev].array(mfi); - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { fab(i,j,k,BaseState::th0_comp) = getRhoThetagivenP(fab(i,j,k,BaseState::p0_comp)) / fab(i,j,k,BaseState::r0_comp); diff --git a/Source/SourceTerms/ERF_buoyancy_utils.H b/Source/SourceTerms/ERF_buoyancy_utils.H index e8db1015e..fd18827f3 100644 --- a/Source/SourceTerms/ERF_buoyancy_utils.H +++ b/Source/SourceTerms/ERF_buoyancy_utils.H @@ -7,61 +7,32 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real -buoyancy_dry_anelastic (int& i, - int& j, - int& k, +buoyancy_dry_anelastic (int& i, int& j, int& k, amrex::Real const& grav_gpu, const amrex::Array4& r0_arr, - const amrex::Array4& p0_arr, + const amrex::Array4& th0_arr, const amrex::Array4& cell_data) { // Note: this is the same term as the moist anelastic buoyancy when qv = qc = qt = 0 - amrex::Real theta_d_lo = cell_data(i,j,k-1,RhoTheta_comp)/r0_arr(i,j,k-1); - amrex::Real theta_d_hi = cell_data(i,j,k ,RhoTheta_comp)/r0_arr(i,j,k); + amrex::Real theta_d_lo = cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp); + amrex::Real theta_d_hi = cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp); amrex::Real theta_d_wface = amrex::Real(0.5) * (theta_d_lo + theta_d_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_d_0_lo = th0_arr(i,j,k-1); + amrex::Real theta_d_0_hi = th0_arr(i,j,k ); amrex::Real theta_d_0_wface = amrex::Real(0.5) * (theta_d_0_lo + theta_d_0_hi); - return (-grav_gpu * (theta_d_wface - theta_d_0_wface) / theta_d_0_wface * 0.5 * (r0_arr(i,j,k) + r0_arr(i,j,k-1))); -} - -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 rho0_wface = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1)); - 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); + return (-rho0_wface * grav_gpu * (theta_d_wface - theta_d_0_wface) / theta_d_0_wface); } - AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real -buoyancy_moist_anelastic (int& i, - int& j, - int& k, +buoyancy_moist_anelastic (int& i, int& j, int& k, amrex::Real const& grav_gpu, amrex::Real const& rv_over_rd, const amrex::Array4& r0_arr, @@ -87,15 +58,15 @@ buoyancy_moist_anelastic (int& i, amrex::Real theta_v_0_wface = amrex::Real(0.5) * (theta_v_0_lo + theta_v_0_hi); - return (-grav_gpu * (theta_v_wface - theta_v_0_wface) / theta_v_0_wface * 0.5 * (r0_arr(i,j,k) + r0_arr(i,j,k-1))); + amrex::Real rho0_wface = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1)); + + return (-rho0_wface * grav_gpu * (theta_v_wface - theta_v_0_wface) / theta_v_0_wface); } AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real -buoyancy_dry_default (int& i, - int& j, - int& k, +buoyancy_dry_default (int& i, int& j, int& k, amrex::Real const& grav_gpu, amrex::Real const& rd_over_cp, const amrex::Array4& r0_arr, @@ -118,13 +89,11 @@ buoyancy_dry_default (int& i, AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real -buoyancy_type1 (int& i, - int& j, - int& k, - const int& n_qstate, - amrex::Real const& grav_gpu, - const amrex::Array4& r0_arr, - const amrex::Array4& cell_data) +buoyancy_rhopert (int& i, int& j, int& k, + const int& n_qstate, + amrex::Real const& grav_gpu, + const amrex::Array4& r0_arr, + const amrex::Array4& cell_data) { amrex::Real rhop_hi = cell_data(i,j,k ,Rho_comp) - r0_arr(i,j,k ); amrex::Real rhop_lo = cell_data(i,j,k-1,Rho_comp) - r0_arr(i,j,k-1); @@ -139,9 +108,7 @@ buoyancy_type1 (int& i, AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real -buoyancy_type2 (int& i, - int& j, - int& k, +buoyancy_type2 (int& i, int& j, int& k, const int& n_qstate, amrex::Real const& grav_gpu, amrex::Real* rho_d_ptr, @@ -191,9 +158,7 @@ buoyancy_type2 (int& i, AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real -buoyancy_type3 (int& i, - int& j, - int& k, +buoyancy_type3 (int& i, int& j, int& k, const int& n_qstate, amrex::Real const& grav_gpu, amrex::Real* rho_d_ptr, @@ -235,9 +200,7 @@ buoyancy_type3 (int& i, AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real -buoyancy_type4 (int& i, - int& j, - int& k, +buoyancy_type4 (int& i, int& j, int& k, const int& n_qstate, amrex::Real const& grav_gpu, amrex::Real* rho_d_ptr, @@ -272,4 +235,56 @@ buoyancy_type4 (int& i, return (-qavg * r0avg * grav_gpu); } + +// ************************************************************************************** +// Routines below this line are not currently used +// ************************************************************************************** + +AMREX_GPU_DEVICE +AMREX_FORCE_INLINE +amrex::Real +buoyancy_dry_Tpert (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& th0_arr, + const amrex::Array4& cell_data) +{ + 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 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; + + 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 +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); +} + #endif diff --git a/Source/SourceTerms/ERF_make_buoyancy.cpp b/Source/SourceTerms/ERF_make_buoyancy.cpp index e66f11ab0..416c79a71 100644 --- a/Source/SourceTerms/ERF_make_buoyancy.cpp +++ b/Source/SourceTerms/ERF_make_buoyancy.cpp @@ -70,7 +70,6 @@ void make_buoyancy (Vector& S_data, // Base state density and pressure 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) { @@ -81,7 +80,7 @@ void make_buoyancy (Vector& S_data, // buoyancy_fab(i, j, k) = buoyancy_dry_anelastic(i,j,k, grav_gpu[2], - r0_arr,p0_arr,cell_data); + r0_arr,th0_arr,cell_data); }); } else { // NOTE: For decomposition in the vertical direction, klo may not @@ -127,7 +126,7 @@ void make_buoyancy (Vector& S_data, ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - buoyancy_fab(i, j, k) = buoyancy_type1(i,j,k,n_q_dry,grav_gpu[2],r0_arr,cell_data); + buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,n_q_dry,grav_gpu[2],r0_arr,cell_data); }); } // mfi @@ -202,7 +201,7 @@ void make_buoyancy (Vector& S_data, ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - buoyancy_fab(i, j, k) = buoyancy_type1(i,j,k,n_qstate, + buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,n_qstate, grav_gpu[2],r0_arr,cell_data); }); } // mfi diff --git a/Source/Utils/ERF_AverageDown.cpp b/Source/Utils/ERF_AverageDown.cpp index 6280169a2..c080d1e1b 100644 --- a/Source/Utils/ERF_AverageDown.cpp +++ b/Source/Utils/ERF_AverageDown.cpp @@ -16,10 +16,17 @@ void ERF::AverageDown () { AMREX_ALWAYS_ASSERT(solverChoice.coupling_type == CouplingType::TwoWay); - int src_comp = 0; - int num_comp = vars_new[0][Vars::cons].nComp(); + + int src_comp, num_comp; for (int lev = finest_level-1; lev >= 0; --lev) { + // If anelastic we don't average down rho because rho == rho0. + if (solverChoice.anelastic[lev]) { + src_comp = 1; + } else { + src_comp = 0; + } + num_comp = vars_new[0][Vars::cons].nComp() - src_comp; AverageDownTo(lev,src_comp,num_comp); } } @@ -28,8 +35,13 @@ ERF::AverageDown () void ERF::AverageDownTo (int crse_lev, int scomp, int ncomp) // NOLINT { - AMREX_ALWAYS_ASSERT(scomp == 0); - AMREX_ALWAYS_ASSERT(ncomp == vars_new[crse_lev][Vars::cons].nComp()); + if (solverChoice.anelastic[crse_lev]) { + AMREX_ALWAYS_ASSERT(scomp == 1); + } else { + AMREX_ALWAYS_ASSERT(scomp == 0); + } + + AMREX_ALWAYS_ASSERT(ncomp == vars_new[crse_lev][Vars::cons].nComp() - scomp); AMREX_ALWAYS_ASSERT(solverChoice.coupling_type == CouplingType::TwoWay); // ******************************************************************************************