diff --git a/Exec/CMakeLists.txt b/Exec/CMakeLists.txt index a5b384ebd..046865f20 100644 --- a/Exec/CMakeLists.txt +++ b/Exec/CMakeLists.txt @@ -19,7 +19,6 @@ elseif (ERF_ENABLE_REGRESSION_TESTS_ONLY) add_subdirectory(DevTests/MovingTerrain) else () add_subdirectory(ABL) - add_subdirectory(SuperCell) add_subdirectory(SquallLine_2D) add_subdirectory(RegTests/Bubble) add_subdirectory(RegTests/Couette_Poiseuille) diff --git a/Exec/RegTests/DensityCurrent/GNUmakefile b/Exec/RegTests/DensityCurrent/GNUmakefile index 75690762b..cda34e9af 100644 --- a/Exec/RegTests/DensityCurrent/GNUmakefile +++ b/Exec/RegTests/DensityCurrent/GNUmakefile @@ -24,6 +24,8 @@ DEBUG = FALSE TEST = TRUE USE_ASSERTION = TRUE +USE_POISSON_SOLVE = TRUE + # GNU Make Bpack := ./Make.package Blocs := . diff --git a/Exec/RegTests/DensityCurrent/inputs_amr b/Exec/RegTests/DensityCurrent/inputs_amr index 11a5373d4..173529de8 100644 --- a/Exec/RegTests/DensityCurrent/inputs_amr +++ b/Exec/RegTests/DensityCurrent/inputs_amr @@ -36,8 +36,8 @@ amr.max_level = 1 # maximum level number allowed amr.n_cell = 128 1 32 # dx=dy=dz=100 m, Straka et al 1993 / Xue et al 2000 erf.fixed_dt = 2.0 # fixed time step [s] -- Straka et al 1993 erf.fixed_fast_dt = 0.5 # fixed time step [s] -- Straka et al 1993 -erf.plot_int_1 = 100 # number of timesteps between plotfiles -erf.check_int =-1000 # number of timesteps between checkpoints +erf.plot_int_1 = 150 # number of timesteps between plotfiles +erf.check_int =-100 # number of timesteps between checkpoints # DIAGNOSTICS & VERBOSITY erf.sum_interval =-1 # timesteps between computing mass @@ -73,7 +73,6 @@ prob.U_0 = 0.0 ################################ MULTILEVEL ################################ amr.ref_ratio_vect = 2 1 2 -erf.coupling_type = "OneWay" erf.regrid_int = 10 erf.refinement_indicators = lo_theta diff --git a/Exec/SuperCell/CMakeLists.txt b/Exec/SuperCell/CMakeLists.txt deleted file mode 100644 index 5fc424233..000000000 --- a/Exec/SuperCell/CMakeLists.txt +++ /dev/null @@ -1,12 +0,0 @@ -set(erf_exe_name erf_super_cell) - -add_executable(${erf_exe_name} "") -target_sources(${erf_exe_name} - PRIVATE - ERF_prob.cpp -) - -target_include_directories(${erf_exe_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) - -include(${CMAKE_SOURCE_DIR}/CMake/BuildERFExe.cmake) -build_erf_exe(${erf_exe_name}) diff --git a/Exec/SuperCell/ERF_prob.H b/Exec/SuperCell/ERF_prob.H deleted file mode 100644 index 94f58ed36..000000000 --- a/Exec/SuperCell/ERF_prob.H +++ /dev/null @@ -1,63 +0,0 @@ -#ifndef ERF_PROB_H_ -#define ERF_PROB_H_ - -#include - -#include "AMReX_REAL.H" - -#include "ERF_prob_common.H" - -struct ProbParm : ProbParmDefaults { - amrex::Real T_0 = 300.0; // surface temperature == mean potential temperature - amrex::Real U_0 = 10.0; - amrex::Real V_0 = 0.0; - amrex::Real x_c = 0.0; // center of thermal perturbation - amrex::Real z_c = 3200.0; - amrex::Real x_r = 1000.0; - amrex::Real z_r = 1000.0; - amrex::Real T_pert = 5.0; // perturbation temperature - // overridden physical constants - amrex::Real C_p = 1004.0; -}; // namespace ProbParm - -class Problem : public ProblemBase -{ -public: - Problem(); - -#include "Prob/ERF_init_density_hse_dry.H" - - void init_custom_pert ( - const amrex::Box& bx, - const amrex::Box& xbx, - const amrex::Box& ybx, - const amrex::Box& zbx, - amrex::Array4 const& state, - amrex::Array4 const& state_pert, - amrex::Array4 const& x_vel_pert, - amrex::Array4 const& y_vel_pert, - amrex::Array4 const& z_vel_pert, - amrex::Array4 const& r_hse, - amrex::Array4 const& p_hse, - amrex::Array4 const& z_nd, - amrex::Array4 const& z_cc, - amrex::GeometryData const& geomdata, - amrex::Array4 const& mf_m, - amrex::Array4 const& mf_u, - amrex::Array4 const& mf_v, - const SolverChoice& sc) override; - - void erf_init_rayleigh ( - amrex::Vector >& rayleigh_ptrs, - amrex::Geometry const& geom, - std::unique_ptr& z_phys_nd, - amrex::Real zdamp) override; - -protected: - std::string name() override { return "Supercell"; } - -private: - ProbParm parms; -}; - -#endif diff --git a/Exec/SuperCell/ERF_prob.cpp b/Exec/SuperCell/ERF_prob.cpp deleted file mode 100644 index 41ed24774..000000000 --- a/Exec/SuperCell/ERF_prob.cpp +++ /dev/null @@ -1,233 +0,0 @@ -#include "ERF_prob.H" -#include "ERF_EOS.H" - -using namespace amrex; - -std::unique_ptr -amrex_probinit( - const amrex_real* /*problo*/, - const amrex_real* /*probhi*/) -{ - return std::make_unique(); -} - -Problem::Problem() -{ - // Parse params - amrex::ParmParse pp("prob"); - pp.query("T_0", parms.T_0); - pp.query("U_0", parms.U_0); - pp.query("x_c", parms.x_c); - pp.query("z_c", parms.z_c); - pp.query("x_r", parms.x_r); - pp.query("z_r", parms.z_r); - pp.query("T_pert", parms.T_pert); - - init_base_parms(parms.rho_0, parms.T_0); -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -amrex::Real -init_supercell_temperature(amrex::Real z, amrex::Real z_0, amrex::Real z_trop, amrex::Real z_top, - amrex::Real T_0, amrex::Real T_trop, amrex::Real T_top) -{ - if (z <= z_trop) { - amrex::Real lapse = - (T_trop - T_0) / (z_trop - z_0); - return T_0 - lapse * (z - z_0); - } else { - amrex::Real lapse = - (T_top - T_trop) / (z_top - z_trop); - return T_trop - lapse * (z - z_trop); - } -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -amrex::Real -init_supercell_pressure(amrex::Real z, amrex::Real z_0, - amrex::Real z_trop, amrex::Real z_top, - amrex::Real T_0, amrex::Real T_trop, amrex::Real T_top) -{ - if (z <= z_trop) { - amrex::Real lapse = - (T_trop - T_0) / (z_trop - z_0); - amrex::Real T = init_supercell_temperature(z, z_0, z_trop, z_top, T_0, T_trop, T_top); - return p_0 * std::pow( T / T_0 , CONST_GRAV/(R_d*lapse) ); - } else { - // Get pressure at the tropopause - amrex::Real lapse = - (T_trop - T_0) / (z_trop - z_0); - amrex::Real p_trop = p_0 * std::pow( T_trop / T_0 , CONST_GRAV/(R_d*lapse) ); - // Get pressure at requested height - lapse = - (T_top - T_trop) / (z_top - z_trop); - if (lapse != 0) { - amrex::Real T = init_supercell_temperature(z, z_0, z_trop, z_top, T_0, T_trop, T_top); - return p_trop * std::pow( T / T_trop , CONST_GRAV/(R_d*lapse) ); - } else { - return p_trop * std::exp(-CONST_GRAV*(z-z_trop)/(R_d*T_trop)); - } - } -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -amrex::Real -init_supercell_sat_mix(amrex::Real press, amrex::Real T ) { - return 380./(press) * std::exp(17.27*(T-273.)/(T-36.)); -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -amrex::Real -init_supercell_relhum(amrex::Real z, amrex::Real z_trop) -{ - if (z <= z_trop) { - return 1. - 0.75 * pow(z / z_trop , 1.25); - } else { - return 0.25; - } -} - -void -Problem::init_custom_pert( - const Box& bx, - const Box& xbx, - const Box& ybx, - const Box& zbx, - Array4 const& /*state*/, - Array4 const& state_pert, - Array4 const& x_vel_pert, - Array4 const& y_vel_pert, - Array4 const& z_vel_pert, - Array4 const& r_hse, - Array4 const& p_hse, - Array4 const& /*z_nd*/, - Array4 const& /*z_cc*/, - GeometryData const& geomdata, - Array4 const& /*mf_m*/, - Array4 const& /*mf_u*/, - Array4 const& /*mf_v*/, - const SolverChoice& sc) -{ - const bool use_moisture = (sc.moisture_type != MoistureType::None); - - const int khi = geomdata.Domain().bigEnd()[2]; - - AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); - - // This is what we do at k = 0 -- note we assume p = p_0 and T = T_0 at z=0 - const amrex::Real& dz = geomdata.CellSize()[2]; - const amrex::Real& prob_lo_z = geomdata.ProbLo()[2]; - const amrex::Real& prob_hi_z = geomdata.ProbHi()[2]; - - const amrex::Real rdOcp = sc.rdOcp; - - // const amrex::Real thetatr = 343.0; - // const amrex::Real theta0 = 300.0; - const amrex::Real ztr = 12000.0; - const amrex::Real Ttr = 213.0; - const amrex::Real Ttop = 213.0; - const amrex::Real deltaz = 1000.*0.0; - const amrex::Real zs = 5000.; - const amrex::Real us = 30.; - const amrex::Real uc = 15.; - - amrex::ParallelForRNG(bx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept - { - // Geometry (note we must include these here to get the data on device) - const auto prob_lo = geomdata.ProbLo(); - const auto dx = geomdata.CellSize(); - const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2]; - - amrex::Real relhum = init_supercell_relhum(z, ztr); - amrex::Real temp = init_supercell_temperature(z, prob_lo_z, ztr, prob_hi_z, parms_d.T_0, Ttr, Ttop); - amrex::Real press = init_supercell_pressure(z, prob_lo_z, ztr, prob_hi_z, parms_d.T_0, Ttr, Ttop); - amrex::Real qvs = init_supercell_sat_mix(press, temp); - -#if 0 - amrex::Real thetaeq; - amrex::Real faceq; - if (z <= ztr) { - thetaeq = theta0 + (thetatr - theta0)*std::pow(z/ztr, 5./4.); - } - else { - thetaeq = thetatr*std::exp(CONST_GRAV*(z-ztr)/c_p*Ttr); - } - - faceq = 1.; - for (int km = 0; km < k; ++km) { - amrex::Real zloc = prob_lo[2]+(km+0.5)*dx[2]; - amrex::Real theq; - if (zloc <= ztr) { - theq = theta0 + (thetatr - theta0)*std::pow(zloc/ztr, 5./4.); - } else { - theq = thetatr*std::exp(CONST_GRAV*(zloc-ztr)/c_p*Ttr); - } - faceq -= CONST_GRAV/(c_p*theq)*dz; - } - - temp = faceq*thetaeq; -#endif - - if (relhum*qvs > 0.014) relhum = 0.014/qvs; - amrex::Real qvapor = std::min(0.014, qvs*relhum); - amrex::Real rho = press/(R_d+qvapor*R_v)/temp; - - // perturb theta - amrex::Real rand_double = amrex::Random(engine) - 1.0; // Random number in [-1,1] - amrex::Real scaling = (khi-static_cast(k))/khi; // Less effect at higher levels - amrex::Real deltaT = parms_d.T_pert*scaling*rand_double; - - amrex::Real theta = getThgivenRandT(rho, temp+deltaT, rdOcp); - - // This version perturbs rho but not p - state_pert(i, j, k, RhoTheta_comp) = rho*theta - getRhoThetagivenP(p_hse(i,j,k)); - state_pert(i, j, k, Rho_comp) = rho - r_hse(i,j,k); - - // Set scalar = 0 everywhere - state_pert(i, j, k, RhoScalar_comp) = 0.0; - - // mean states - if (use_moisture) { - state_pert(i, j, k, RhoQ1_comp) = rho*qvapor; - state_pert(i, j, k, RhoQ2_comp) = 0.0; - } - }); - - // Set the x-velocity - amrex::ParallelFor(xbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - const amrex::Real z = prob_lo_z + (k+0.5) * dz; - if (z < zs-deltaz) { - x_vel_pert(i, j, k) = us*(z/zs) - uc; - } else if (std::abs(z-zs) < deltaz) { - x_vel_pert(i, j, k) = (-0.8+3.*(z/zs)-1.25*(z/zs)*(z/zs))*us-uc; - } else { - x_vel_pert(i, j, k) = us-uc; - } - }); - - // Set the y-velocity - amrex::ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - y_vel_pert(i, j, k) = 0.0; - }); - - // Set the z-velocity - amrex::ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - z_vel_pert(i, j, k) = 0.0; - }); - - amrex::Gpu::streamSynchronize(); -} - -void -Problem::erf_init_rayleigh( - amrex::Vector >& rayleigh_ptrs, - amrex::Geometry const& geom, - std::unique_ptr& /*z_phys_nd*/, - amrex::Real /*zdamp*/) -{ - const int khi = geom.Domain().bigEnd()[2]; - - // We just use these values to test the Rayleigh damping - for (int k = 0; k <= khi; k++) - { - rayleigh_ptrs[Rayleigh::ubar][k] = 2.0; - rayleigh_ptrs[Rayleigh::vbar][k] = 1.0; - rayleigh_ptrs[Rayleigh::wbar][k] = 0.0; - rayleigh_ptrs[Rayleigh::thetabar][k] = parms.T_0; - } -} diff --git a/Exec/SuperCell/GNUmakefile b/Exec/SuperCell/GNUmakefile deleted file mode 100644 index 968d704e3..000000000 --- a/Exec/SuperCell/GNUmakefile +++ /dev/null @@ -1,32 +0,0 @@ -# AMReX -COMP = gnu -PRECISION = DOUBLE - -# Profiling -PROFILE = FALSE -TINY_PROFILE = FALSE -COMM_PROFILE = FALSE -TRACE_PROFILE = FALSE -MEM_PROFILE = FALSE -USE_GPROF = FALSE - -# Performance -USE_MPI = TRUE -USE_OMP = FALSE - -USE_CUDA = FALSE -USE_HIP = FALSE -USE_SYCL = FALSE - -# Debugging -DEBUG = FALSE - -TEST = TRUE -USE_ASSERTION = TRUE - -# GNU Make -Bpack := ./Make.package -Blocs := . -ERF_HOME := ../.. -ERF_PROBLEM_DIR = $(ERF_HOME)/Exec/SuperCell -include $(ERF_HOME)/Exec/Make.ERF diff --git a/Exec/SuperCell/Make.package b/Exec/SuperCell/Make.package deleted file mode 100644 index 5fc21f61c..000000000 --- a/Exec/SuperCell/Make.package +++ /dev/null @@ -1,2 +0,0 @@ -CEXE_headers += ERF_prob.H -CEXE_sources += ERF_prob.cpp diff --git a/Exec/SuperCell/README b/Exec/SuperCell/README deleted file mode 100644 index 9242958db..000000000 --- a/Exec/SuperCell/README +++ /dev/null @@ -1,2 +0,0 @@ -This problem setup is the evolution of a supercell, which primarily tests the ability -of ERF to model moisture physics. diff --git a/Exec/SuperCell/inputs_moisture b/Exec/SuperCell/inputs_moisture deleted file mode 100644 index 92e812fe7..000000000 --- a/Exec/SuperCell/inputs_moisture +++ /dev/null @@ -1,63 +0,0 @@ -# ------------------ INPUTS TO MAIN PROGRAM ------------------- -max_step = 1000 -stop_time = 90000.0 - -amrex.fpe_trap_invalid = 1 - -fabarray.mfiter_tile_size = 2048 1024 2048 - -# PROBLEM SIZE & GEOMETRY -geometry.prob_lo = -25600. 0. 0. -geometry.prob_hi = 25600. 400. 12800. -amr.n_cell = 128 4 32 # dx=dy=dz=100 m - -# periodic in x to match WRF setup -# - as an alternative, could use symmetry at x=0 and outflow at x=25600 -geometry.is_periodic = 1 1 0 -zlo.type = "SlipWall" -zhi.type = "SlipWall" - -# TIME STEP CONTROL -erf.fixed_dt = 1.0 # fixed time step [s] -- Straka et al 1993 -erf.fixed_fast_dt = 0.25 # fixed time step [s] -- Straka et al 1993 - -# DIAGNOSTICS & VERBOSITY -erf.sum_interval = 1 # timesteps between computing mass -erf.v = 1 # verbosity in ERF.cpp -amr.v = 1 # verbosity in Amr.cpp - -# REFINEMENT / REGRIDDING -amr.max_level = 0 # maximum level number allowed - -# CHECKPOINT FILES -amr.check_file = chk # root name of checkpoint file -amr.check_int = 10000 # number of timesteps between checkpoints -#amr.restart = chk01000 - -# PLOTFILES -erf.plot_file_1 = plt # root name of plotfile -erf.plot_int_1 = 1 # number of timesteps between plotfiles -erf.plot_vars_1 = density rhotheta rhoQ1 rhoQ2 rhoQ3 x_velocity y_velocity z_velocity pressure theta temp qt qp qv qc qi - -# SOLVER CHOICE -erf.use_gravity = true -erf.use_coriolis = false - -erf.moisture_model = "SAM" - -erf.les_type = "Deardorff" -erf.KE_0 = 0.1 # for Deardorff -#erf.les_type = "None" -# -# diffusion coefficient from Straka, K = 75 m^2/s -# -#erf.molec_diff_type = "ConstantAlpha" -erf.molec_diff_type = "None" -erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities -erf.dynamicViscosity = 75.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s -erf.alpha_T = 75.0 # [m^2/s] - -# PROBLEM PARAMETERS (optional) -prob.T_0 = 300.0 -prob.U_0 = 0 -prob.T_pert = 3 diff --git a/Source/BoundaryConditions/ERF_FillIntermediatePatch.cpp b/Source/BoundaryConditions/ERF_FillIntermediatePatch.cpp index 6b8c013e2..a494132f1 100644 --- a/Source/BoundaryConditions/ERF_FillIntermediatePatch.cpp +++ b/Source/BoundaryConditions/ERF_FillIntermediatePatch.cpp @@ -121,10 +121,6 @@ ERF::FillIntermediatePatch (int lev, Real time, Vector ctime = {t_old[lev-1], t_new[lev-1]}; Vector ftime = {time,time}; - // Impose physical bc's on coarse data (note time and 0 are not used) - (*physbcs_cons[lev-1])(vars_old[lev-1][Vars::cons],icomp_cons,ncomp_cons,IntVect{ng_cons},time,BCVars::cons_bc,true); - (*physbcs_cons[lev-1])(vars_new[lev-1][Vars::cons],icomp_cons,ncomp_cons,IntVect{ng_cons},time,BCVars::cons_bc,true); - // Impose physical bc's on fine data (note time and 0 are not used) (*physbcs_cons[lev])(*mfs_vel[Vars::cons],icomp_cons,ncomp_cons,IntVect{ng_cons},time,BCVars::cons_bc,true); @@ -257,10 +253,6 @@ ERF::FillIntermediatePatch (int lev, Real time, fmf = {&mfu,&mfu}; cmf = {&vars_old[lev-1][Vars::xvel], &vars_new[lev-1][Vars::xvel]}; - // Impose physical bc's on coarse data (note time and 0 are not used) - (*physbcs_u[lev-1])(vars_old[lev-1][Vars::xvel],0,1,IntVect{ng_vel},time,BCVars::xvel_bc,true); - (*physbcs_u[lev-1])(vars_new[lev-1][Vars::xvel],0,1,IntVect{ng_vel},time,BCVars::xvel_bc,true); - // Call FillPatchTwoLevels which ASSUMES that all ghost cells have already been filled FillPatchTwoLevels(mfu, IntVect{ng_vel}, IntVect(0,0,0), time, cmf, ctime, fmf, ftime, @@ -275,10 +267,6 @@ ERF::FillIntermediatePatch (int lev, Real time, fmf = {&mfv,&mfv}; cmf = {&vars_old[lev-1][Vars::yvel], &vars_new[lev-1][Vars::yvel]}; - // Impose physical bc's on coarse data (note time and 0 are not used) - (*physbcs_v[lev-1])(vars_old[lev-1][Vars::yvel],0,1,IntVect{ng_vel},time,BCVars::yvel_bc,true); - (*physbcs_v[lev-1])(vars_new[lev-1][Vars::yvel],0,1,IntVect{ng_vel},time,BCVars::yvel_bc,true); - // Call FillPatchTwoLevels which ASSUMES that all ghost cells have already been filled FillPatchTwoLevels(mfv, IntVect{ng_vel}, IntVect(0,0,0), time, cmf, ctime, fmf, ftime, @@ -293,16 +281,6 @@ ERF::FillIntermediatePatch (int lev, Real time, fmf = {&mfw,&mfw}; cmf = {&vars_old[lev-1][Vars::zvel], &vars_new[lev-1][Vars::zvel]}; - // Impose physical bc's on coarse data (note time and 0 are not used) - (*physbcs_w[lev-1])(vars_old[lev-1][Vars::zvel], - vars_old[lev-1][Vars::xvel], - vars_old[lev-1][Vars::yvel], - IntVect{ng_vel},time,BCVars::zvel_bc,true); - (*physbcs_w[lev-1])(vars_new[lev-1][Vars::zvel], - vars_new[lev-1][Vars::xvel], - vars_new[lev-1][Vars::yvel], - IntVect{ng_vel},time,BCVars::zvel_bc,true); - // Call FillPatchTwoLevels which ASSUMES that all ghost cells have already been filled FillPatchTwoLevels(mfw, IntVect{ng_vel}, IntVect(0,0,0), time, cmf, ctime, fmf, ftime, diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index 645b5c7a5..0eff71724 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -339,7 +339,7 @@ struct SolverChoice { // Which type of multilevel coupling - coupling_type = CouplingType::OneWay; // Default + coupling_type = CouplingType::TwoWay; // Default pp.query_enum_case_insensitive("coupling_type",coupling_type); // Which type of windfarm model diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index 3f0fd47cf..377eb09be 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -332,39 +332,40 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp // ***************************************************************************************************** make_physbcs(lev); - // ************************************************************************************************* - // 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 - // ************************************************************************************************* - 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); - // ******************************************************************************************** // 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; + 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}; + 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); + // 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()); + // 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]); - } + // ************************************************************************************************* + // 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); + + // Now swap the pointers since we needed both old and new in the FillPatch + std::swap(temp_base_state, base_state[lev]); // ******************************************************************************************** // Copy from new into old just in case @@ -445,6 +446,10 @@ ERF::ClearLevel (int lev) rW_new[lev].clear(); rW_old[lev].clear(); + if (lev > 0) { + zmom_crse_rhs[lev].clear(); + } + if (solverChoice.anelastic[lev] == 1) { pp_inc[lev].clear(); } diff --git a/Source/TimeIntegration/ERF_TI_fast_headers.H b/Source/TimeIntegration/ERF_TI_fast_headers.H index 5c470d9da..8a244dbcb 100644 --- a/Source/TimeIntegration/ERF_TI_fast_headers.H +++ b/Source/TimeIntegration/ERF_TI_fast_headers.H @@ -26,6 +26,7 @@ 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 1213b8cf6..152c60227 100644 --- a/Source/TimeIntegration/ERF_TI_fast_rhs_fun.H +++ b/Source/TimeIntegration/ERF_TI_fast_rhs_fun.H @@ -157,6 +157,7 @@ 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], @@ -166,6 +167,7 @@ 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_TI_slow_headers.H b/Source/TimeIntegration/ERF_TI_slow_headers.H index 86dfb38b5..14ab70cde 100644 --- a/Source/TimeIntegration/ERF_TI_slow_headers.H +++ b/Source/TimeIntegration/ERF_TI_slow_headers.H @@ -67,8 +67,6 @@ 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, diff --git a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H index 092f87589..b96c2e096 100644 --- a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H @@ -125,8 +125,6 @@ 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(), @@ -238,8 +236,6 @@ 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(), @@ -449,8 +445,6 @@ 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(), diff --git a/Source/TimeIntegration/ERF_fast_rhs_N.cpp b/Source/TimeIntegration/ERF_fast_rhs_N.cpp index 5bafd34bf..b11e62254 100644 --- a/Source/TimeIntegration/ERF_fast_rhs_N.cpp +++ b/Source/TimeIntegration/ERF_fast_rhs_N.cpp @@ -41,6 +41,7 @@ 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, @@ -191,17 +192,17 @@ void erf_fast_rhs_N (int step, int nrk, // ********************************************************************* // Define updates in the RHS of {x, y, z}-momentum equations // ********************************************************************* - if (nrk == 0 and step == 0) { + if (nrk == 0 and step == 0) { // prev == stage ParallelFor(tbx, tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real new_drho_u = prev_xmom(i,j,k) - stage_xmom(i,j,k) + dtau * slow_rhs_rho_u(i,j,k); + Real new_drho_u = dtau * slow_rhs_rho_u(i,j,k); avg_xmom(i,j,k) += facinv*new_drho_u; temp_cur_xmom_arr(i,j,k) = stage_xmom(i,j,k) + new_drho_u; }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real new_drho_v = prev_ymom(i,j,k) - stage_ymom(i,j,k) + dtau * slow_rhs_rho_v(i,j,k); + Real new_drho_v = dtau * slow_rhs_rho_v(i,j,k); avg_ymom(i,j,k) += facinv*new_drho_v; temp_cur_ymom_arr(i,j,k) = stage_ymom(i,j,k) + new_drho_v; }); @@ -220,11 +221,10 @@ void erf_fast_rhs_N (int step, int nrk, } Real pi_c = 0.5 * (pi_stage_ca(i-1,j,k,0) + pi_stage_ca(i,j,k,0)); - Real fast_rhs_rho_u = -Gamma * R_d * pi_c * gpx; Real new_drho_u = prev_xmom(i,j,k) - stage_xmom(i,j,k) - + dtau * fast_rhs_rho_u + dtau * slow_rhs_rho_u(i,j,k); + + dtau * fast_rhs_rho_u + dtau * slow_rhs_rho_u(i,j,k); avg_xmom(i,j,k) += facinv*new_drho_u; @@ -243,11 +243,10 @@ void erf_fast_rhs_N (int step, int nrk, } Real pi_c = 0.5 * (pi_stage_ca(i,j-1,k,0) + pi_stage_ca(i,j,k,0)); - Real fast_rhs_rho_v = -Gamma * R_d * pi_c * gpy; Real new_drho_v = prev_ymom(i,j,k) - stage_ymom(i,j,k) - + dtau * fast_rhs_rho_v + dtau * slow_rhs_rho_v(i,j,k); + + dtau * fast_rhs_rho_v + dtau * slow_rhs_rho_v(i,j,k); avg_ymom(i,j,k) += facinv*new_drho_v; @@ -282,6 +281,8 @@ 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); @@ -422,13 +423,41 @@ 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) - + dtau * slow_rhs_rho_w(i,j,lo.z); + // 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); - // 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) - + dtau * slow_rhs_rho_w(i,j,hi.z+1); +#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 + + // 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 }); // b2d #ifdef AMREX_USE_GPU diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index c92ed95cb..b35f3329e 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -32,6 +32,7 @@ using namespace amrex; * @param[in] xmom_src source terms for x-momentum * @param[in] ymom_src source terms for y-momentum * @param[in] zmom_src source terms for z-momentum + * @param[in] zmom_crse_rhs update term from coarser level for z-momentum; non-zero on c/f boundary only * @param[in] Tau11 tau_11 component of stress tensor * @param[in] Tau22 tau_22 component of stress tensor * @param[in] Tau33 tau_33 component of stress tensor @@ -81,9 +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* xmom_crse_rhs, - const MultiFab* ymom_crse_rhs, - 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, @@ -746,6 +745,34 @@ void erf_slow_rhs_pre (int level, int finest_level, } }); +#if 0 + auto const lo = lbound(bx); + auto const hi = ubound(bx); + + // Note: the logic below assumes no tiling in z! + if (level > 0) { + + const Array4& rho_w_rhs_crse = zmom_crse_rhs->const_array(mfi); + + 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.) { + 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.) { + rho_w_rhs(i,j,hi.z+1) = rho_w_rhs_crse(i,j,hi.z+1); + } + }); + } +#endif + { BL_PROFILE("slow_rhs_pre_fluxreg"); // We only add to the flux registers in the final RK step