Fix the gridding of level 0 if needed, and add DevTests/FlowInABox

Exec/DevTests/FlowInABox/CMakeLists.txt

@@ -0,0 +1,15 @@
+set(erf_exe_name erf_flow_in_a_box)
+
+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})
+
+#find_package(AMReX)
+#target_link_libraries( ${_target} AMReX::amrex)

Exec/DevTests/FlowInABox/ERF_prob.H

@@ -0,0 +1,79 @@
+#ifndef ERF_PROB_H_
+#define ERF_PROB_H_
+
+#include <string>
+
+#include "AMReX_REAL.H"
+
+#include "ERF_prob_common.H"
+
+struct ProbParm : ProbParmDefaults {
+ amrex::Real rho_0 = 0.0;
+ amrex::Real T_0 = 0.0;
+ amrex::Real A_0 = 1.0;
+ amrex::Real KE_0 = 0.1;
+
+ amrex::Real KE_decay_height = -1;
+ amrex::Real KE_decay_order = 1;
+
+ amrex::Real U_0 = 0.0;
+ amrex::Real V_0 = 0.0;
+ amrex::Real W_0 = 0.0;
+
+ // random initial perturbations (legacy code)
+ amrex::Real U_0_Pert_Mag = 0.0;
+ amrex::Real V_0_Pert_Mag = 0.0;
+ amrex::Real W_0_Pert_Mag = 0.0;
+ amrex::Real T_0_Pert_Mag = 0.0; // perturbation to rho*Theta
+ bool pert_rhotheta = true;
+
+ // divergence-free initial perturbations
+ amrex::Real pert_deltaU = 0.0;
+ amrex::Real pert_deltaV = 0.0;
+ amrex::Real pert_periods_U = 5.0;
+ amrex::Real pert_periods_V = 5.0;
+ amrex::Real pert_ref_height = 100.0;
+
+ // helper vars
+ amrex::Real aval;
+ amrex::Real bval;
+ amrex::Real ufac;
+ amrex::Real vfac;
+}; // namespace ProbParm
+
+class Problem : public ProblemBase
+{
+public:
+ Problem(const amrex::Real* problo, const amrex::Real* probhi);
+
+#include "Prob/ERF_init_constant_density_hse.H"
+#include "Prob/ERF_init_rayleigh_damping.H"
+
+ void init_custom_pert (
+ const amrex::Box& bx,
+ const amrex::Box& xbx,
+ const amrex::Box& ybx,
+ const amrex::Box& zbx,
+ amrex::Array4<amrex::Real const> const& state,
+ amrex::Array4<amrex::Real > const& state_pert,
+ amrex::Array4<amrex::Real > const& x_vel_pert,
+ amrex::Array4<amrex::Real > const& y_vel_pert,
+ amrex::Array4<amrex::Real > const& z_vel_pert,
+ amrex::Array4<amrex::Real > const& r_hse,
+ amrex::Array4<amrex::Real > const& p_hse,
+ amrex::Array4<amrex::Real const> const& z_nd,
+ amrex::Array4<amrex::Real const> const& z_cc,
+ amrex::GeometryData const& geomdata,
+ amrex::Array4<amrex::Real const> const& mf_m,
+ amrex::Array4<amrex::Real const> const& mf_u,
+ amrex::Array4<amrex::Real const> const& mf_v,
+ const SolverChoice& sc) override;
+
+protected:
+ std::string name() override { return "ABL"; }
+
+private:
+ ProbParm parms;
+};
+
+#endif

Exec/DevTests/FlowInABox/ERF_prob.cpp

@@ -0,0 +1,100 @@
+#include "ERF_prob.H"
+#include "AMReX_Random.H"
+
+using namespace amrex;
+
+std::unique_ptr<ProblemBase>
+amrex_probinit(const amrex_real* problo, const amrex_real* probhi)
+{
+ return std::make_unique<Problem>(problo, probhi);
+}
+
+Problem::Problem(const amrex::Real* problo, const amrex::Real* probhi)
+{
+ // Parse params
+ ParmParse pp("prob");
+ pp.query("rho_0", parms.rho_0);
+ pp.query("T_0", parms.T_0);
+ pp.query("A_0", parms.A_0);
+ pp.query("KE_0", parms.KE_0);
+ pp.query("KE_decay_height", parms.KE_decay_height);
+ pp.query("KE_decay_order", parms.KE_decay_order);
+
+ pp.query("U_0", parms.U_0);
+ pp.query("V_0", parms.V_0);
+ pp.query("W_0", parms.W_0);
+ pp.query("U_0_Pert_Mag", parms.U_0_Pert_Mag);
+ pp.query("V_0_Pert_Mag", parms.V_0_Pert_Mag);
+ pp.query("W_0_Pert_Mag", parms.W_0_Pert_Mag);
+ pp.query("T_0_Pert_Mag", parms.T_0_Pert_Mag);
+ pp.query("pert_rhotheta", parms.pert_rhotheta);
+
+ pp.query("pert_deltaU", parms.pert_deltaU);
+ pp.query("pert_deltaV", parms.pert_deltaV);
+ pp.query("pert_periods_U", parms.pert_periods_U);
+ pp.query("pert_periods_V", parms.pert_periods_V);
+ pp.query("pert_ref_height", parms.pert_ref_height);
+ parms.aval = parms.pert_periods_U * 2.0 * PI / (probhi[1] - problo[1]);
+ parms.bval = parms.pert_periods_V * 2.0 * PI / (probhi[0] - problo[0]);
+ parms.ufac = parms.pert_deltaU * std::exp(0.5) / parms.pert_ref_height;
+ parms.vfac = parms.pert_deltaV * std::exp(0.5) / parms.pert_ref_height;
+
+ Real p_0 = 1.e5;
+ Real R_d = 287.0;
+ Real c_p = 1004.5;
+ Real rdOcp = R_d / c_p;
+
+ // Read in T_0 as temperature not theta then convert to theta
+ Real th_0 = getThgivenPandT(parms.T_0, p_0, rdOcp);
+ Real r_0 = p_0 / (R_d * parms.T_0);
+
+ amrex::Print() << "READING IN T0 as " << parms.T_0 << std::endl;
+ amrex::Print() << "COMPUTING THETA to be " << th_0 << std::endl;
+
+ init_base_parms(r_0, th_0);
+}
+
+void
+Problem::init_custom_pert(
+ const amrex::Box& bx,
+ const amrex::Box& xbx,
+ const amrex::Box& ybx,
+ const amrex::Box& zbx,
+ amrex::Array4<amrex::Real const> const& /*state*/,
+ amrex::Array4<amrex::Real > const& state_pert,
+ amrex::Array4<amrex::Real > const& x_vel_pert,
+ amrex::Array4<amrex::Real > const& y_vel_pert,
+ amrex::Array4<amrex::Real > const& z_vel_pert,
+ amrex::Array4<amrex::Real > const& r_hse,
+ amrex::Array4<amrex::Real > const& /*p_hse*/,
+ amrex::Array4<amrex::Real const> const& z_nd,
+ amrex::Array4<amrex::Real const> const& z_cc,
+ amrex::GeometryData const& geomdata,
+ amrex::Array4<amrex::Real const> const& /*mf_m*/,
+ amrex::Array4<amrex::Real const> const& /*mf_u*/,
+ amrex::Array4<amrex::Real const> const& /*mf_v*/,
+ const SolverChoice& sc)
+{
+
+ // Add temperature perturbations
+ ParallelForRNG(bx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept
+ {
+ Real rand_double = amrex::Random(engine); // Between 0.0 and 1.0
+ state_pert(i, j, k, RhoTheta_comp) = (rand_double*2.0 - 1.0)*parms_d.T_0_Pert_Mag;
+ state_pert(i, j, k, RhoTheta_comp) *= r_hse(i,j,k);
+
+ });
+
+ // Set the x-velocity
+ ParallelFor(xbx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+ x_vel_pert(i, j, k) = 0.0;
+ });
+ // Set the y-velocity
+ ParallelFor(ybx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+ y_vel_pert(i, j, k) = 0.0;
+ });
+ // Set the z-velocity
+ ParallelFor(zbx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+ z_vel_pert(i, j, k) = 0.0;
+ });
+}

Exec/DevTests/FlowInABox/GNUmakefile

@@ -0,0 +1,35 @@
+# 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
+
+USE_POISSON_SOLVE = TRUE
+
+# GNU Make
+Bpack := ./Make.package
+Blocs := .
+
+ERF_HOME := ../../..
+ERF_PROBLEM_DIR = $(ERF_HOME)/Exec/DevTests/FlowInABox
+include $(ERF_HOME)/Exec/Make.ERF

Exec/DevTests/FlowInABox/Make.package

@@ -0,0 +1,2 @@
+CEXE_headers += ERF_prob.H
+CEXE_sources += ERF_prob.cpp

Exec/DevTests/FlowInABox/inputs

@@ -0,0 +1,80 @@
+# ------------------ INPUTS TO MAIN PROGRAM -------------------
+stop_time = 370.
+
+amrex.fpe_trap_invalid = 1
+
+erf.anelastic = 1
+
+erf.init_type = "uniform"
+
+fabarray.mfiter_tile_size = 1024 1024 1024
+
+# PROBLEM SIZE & GEOMETRY
+geometry.prob_lo = -1. -1. 0.
+geometry.prob_hi = 1. 1. 1.
+
+#coarse
+amr.n_cell = 64 64 32
+
+#fine
+amr.n_cell = 128 128 64
+
+geometry.is_periodic = 0 0 0
+
+xlo.type = "SlipWall"
+xhi.type = "SlipWall"
+ylo.type = "SlipWall"
+yhi.type = "SlipWall"
+zlo.type = "SlipWall"
+zhi.type = "SlipWall"
+
+xlo.theta = 285.
+xhi.theta = 285.
+ylo.theta = 285.
+yhi.theta = 285.
+zlo.theta = 299.
+zhi.theta = 280.
+
+xlo.density = 1.2225686
+xhi.density = 1.2225686
+ylo.density = 1.2225686
+yhi.density = 1.2225686
+zlo.density = 1.2225686
+zhi.density = 1.2225686
+
+# TIME STEP CONTROL
+#erf.cfl = 0.9
+erf.fixed_dt = 0.025
+
+# 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
+erf.check_file = chk # root name of checkpoint file
+erf.check_int = 1000 # number of timesteps between checkpoints
+
+# PLOTFILES
+erf.plot_file_1 = plt # prefix of plotfile name
+erf.plot_int_1 = 100 # number of timesteps between plotfiles
+erf.plot_vars_1 = density x_velocity y_velocity z_velocity pressure temp theta
+
+erf.use_gravity = true
+
+# SOLVER CHOICE
+erf.molec_diff_type = "None"
+erf.alpha_T = 0.0
+erf.alpha_C = 0.0
+
+erf.les_type = "Smagorinsky"
+erf.Cs = 0.1
+erf.Pr_t = 0.33333333333333
+
+# PROBLEM PARAMETERS
+prob.rho_0 = 1.0
+prob.T_0 = 285.
+prob.T_0_Pert_Mag = 0.1

Source/ERF_make_new_level.cpp

@@ -20,11 +20,27 @@ using namespace amrex;
 // (overrides the pure virtual function in AmrCore)
 // main.cpp --> ERF::InitData --> InitFromScratch --> MakeNewGrids --> MakeNewLevelFromScratch
 // restart --> MakeNewGrids --> MakeNewLevelFromScratch
-void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba,
- const DistributionMapping& dm)
+void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba_in,
+ const DistributionMapping& dm_in)
 {
- // Set BoxArray grids and DistributionMapping dmap in AMReX_AmrMesh.H class
+ BoxArray ba;
+ DistributionMapping dm;
+ if (lev == 0 && ba_in.size() != ParallelDescriptor::NProcs())
+ {
+ // We only decompose in z if max_grid_size_z indicates we should
+ bool decompose_in_z = (max_grid_size[0][2] < Geom(0).Domain().length(2));
+
+ ba = ERFPostProcessBaseGrids(Geom(0).Domain(),decompose_in_z);
+ dm = DistributionMapping(ba);
+ } else {
+ ba = ba_in;
+ dm = dm_in;
+ }
+
+ // Define grids[lev] to be ba
  SetBoxArray(lev, ba);
+
+ // Define dmap[lev] to be dm
  SetDistributionMap(lev, dm);
 
  // amrex::Print() <<" BA FROM SCRATCH AT LEVEL " << lev << " " << ba << std::endl;

Source/ERF_prob_common.H

@@ -453,7 +453,6 @@ public:
  amrex::Real T_0 = base_parms.T_0;
  amrex::Print() << "Initializing uniform fields"
  << " rho=" << rho_0 << " theta=" << T_0
- << " -- this probably only makes sense with gravity turned off"
  << std::endl;
 
  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept

Source/Utils/ERF_ChopGrids.cpp

@@ -2,6 +2,21 @@
 
 using namespace amrex;
 
+BoxArray
+ERFPostProcessBaseGrids (const Box& domain, bool decompose_in_z)
+{
+ //
+ // This is used to avoid the case where the native amrex decomposition makes
+ // too many grids for the number of processors.
+ //
+ // The idea is to not override the user preference if expressed by max_grid_size
+ // but instead to ensure that the default behavior is what we want.
+ //
+ BoxArray ba0 = amrex::decompose(domain, ParallelDescriptor::NProcs(),
+ {true,true,decompose_in_z});
+ return ba0;
+}
+
 void
 ChopGrids2D (BoxArray& ba, const Box& domain, int target_size)
 {

Source/Utils/ERF_Utils.H

@@ -14,6 +14,13 @@
  */
 void ChopGrids2D (amrex::BoxArray& ba, const amrex::Box& domain, int target_size);
 
+/*
+ * Create a new BoxArray with exactly the number of boxes as MPI ranks.
+ * This is designed to work only on the base level which covers the entire domain.
+ */
+amrex::BoxArray
+ERFPostProcessBaseGrids (const amrex::Box& domain, bool decompose_in_z);
+
 /*
  * Create the Jacobian for the metric transformation when use_terrain is true
  */