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

…rf-model#1852)

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
  */