add a compressible version of the xrb_layered problem

.github/workflows/good_defines.txt

@@ -34,6 +34,7 @@ NSE_TABLE
 RADIATION
 RAD_INTERP
 REACTIONS
+RNG_STATE_INIT
 ROTATION
 SCREENING
 SHOCK_VAR

Exec/Make.Castro

@@ -172,6 +172,10 @@ ifeq ($(USE_MODEL_PARSER), TRUE)
   DEFINES += -DMODEL_PARSER -DNPTS_MODEL=$(MAX_NPTS_MODEL) -DNUM_MODELS=$(NUM_MODELS)
 endif
 
+ifeq ($(USE_RNG_STATE_INIT), TRUE)
+  DEFINES += -DRNG_STATE_INIT
+endif
+
 # add / define any special physics we need
 
 ifeq ($(USE_MHD), TRUE)

Exec/science/xrb_layered/GNUmakefile

@@ -0,0 +1,28 @@
+PRECISION        = DOUBLE
+PROFILE          = FALSE
+DEBUG            = FALSE
+DIM              = 2
+
+COMP	         = gnu
+
+USE_MPI          = TRUE
+USE_OMP          = FALSE
+
+USE_GRAV         = TRUE
+USE_REACT        = TRUE
+
+USE_CXX_MODEL_PARSER = TRUE
+USE_RNG_STATE_INIT = TRUE
+
+CASTRO_HOME = ../../..
+
+# This sets the EOS directory in $(MICROPHYSICS_HOME)/eos
+EOS_DIR     := helmholtz
+
+# This sets the EOS directory in $(MICROPHYSICS_HOME)/networks
+NETWORK_DIR := CNO_extras
+
+Bpack   := ./Make.package
+Blocs   := .
+
+include $(CASTRO_HOME)/Exec/Make.Castro

Exec/science/xrb_layered/Make.package

@@ -0,0 +1 @@
+

Exec/science/xrb_layered/Problem.H

@@ -0,0 +1,2 @@
+amrex::Real maxVal (const std::string& name, amrex::Real time);
+

Exec/science/xrb_layered/Problem_Derive.H

@@ -0,0 +1,5 @@
+void ca_dertpert
+  (const amrex::Box& bx, amrex::FArrayBox& derfab, int dcomp, int /*ncomp*/,
+   const amrex::FArrayBox& datfab, const amrex::Geometry& geomdata,
+   amrex::Real /*time*/, const int* /*bcrec*/, int /*level*/);
+

Exec/science/xrb_layered/Problem_Derive.cpp

@@ -0,0 +1,54 @@
+#include <AMReX_REAL.H>
+
+#include <Derive.H>
+#include <Problem_Derive_F.H>
+#include <Castro.H>
+#include <Castro_F.H>
+#include <model_parser.H>
+
+using namespace amrex;
+
+
+void ca_dertpert(const Box& bx, FArrayBox& derfab, int dcomp, int /*ncomp*/,
+                 const FArrayBox& datfab, const Geometry& geomdata,
+                 Real /*time*/, const int* /*bcrec*/, int /*level*/)
+{
+
+    // derive the temperature perturbation
+
+    const auto dx = geomdata.CellSizeArray();
+    const auto problo = geomdata.ProbLoArray();
+
+    auto const dat = datfab.array();
+    auto const der = derfab.array();
+
+    amrex::ParallelFor(bx,
+    [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
+    {
+
+        Real x = problo[0] + dx[0] * (static_cast<Real>(i) + 0.5_rt);
+
+        Real y = 0.0;
+#if AMREX_SPACEDIM >= 2
+        y = problo[1] + dx[1] * (static_cast<Real>(j) + 0.5_rt);
+#endif
+
+        Real z = 0.0;
+#if AMREX_SPACEDIM == 3
+        z = problo[2] + dx[2] * (static_cast<Real>(k) + 0.5_rt);
+#endif
+
+#if AMREX_SPACEDIM == 2
+        Real height = y;
+#else
+        Real height = z;
+#endif
+
+        Real temp = interpolate(height, model::itemp);
+
+        der(i,j,k,0) = dat(i,j,k,UTEMP) - temp;
+
+    });
+
+}
+

Exec/science/xrb_layered/Problem_Derives.H

@@ -0,0 +1,5 @@
+    //
+    // tpert
+    //
+    derive_lst.add("tpert",IndexType::TheCellType(),1,ca_dertpert,the_same_box);
+    derive_lst.addComponent("tpert",desc_lst,State_Type,URHO,NUM_STATE);

Exec/science/xrb_layered/README.md

@@ -0,0 +1,3 @@
+# xrb_layered
+
+The initial model comes from Simon's MAESTROeX setup

Exec/science/xrb_layered/_prob_params

@@ -0,0 +1,11 @@
+H_min                 real           1.e-4_rt      y
+
+cutoff_density        real           50.0_rt       y
+
+model_name            character      ""            y
+
+apply_perturbation    integer        1             y
+
+pert_density          real           1.0           y
+
+model_shift           real           0.0           y

Exec/science/xrb_layered/inputs_2d

@@ -0,0 +1,88 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 1000000
+stop_time =  1000.0
+
+# PROBLEM SIZE & GEOMETRY
+geometry.is_periodic = 1       0
+geometry.coord_sys   = 0                  # 0 => cart, 1 => RZ  2=>spherical
+geometry.prob_lo     = 0.0     0.0
+geometry.prob_hi     = 3072    3072
+amr.n_cell           = 128     128
+
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+# 0 = Interior           3 = Symmetry
+# 1 = Inflow             4 = SlipWall
+# 2 = Outflow            5 = NoSlipWall
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+castro.lo_bc       =  0   3
+castro.hi_bc       =  0   2
+
+castro.fill_ambient_bc = 1
+castro.ambient_fill_dir = 1
+castro.ambient_outflow_vel = 1
+
+
+# WHICH PHYSICS
+castro.do_hydro = 1
+castro.do_react = 1
+castro.add_ext_src   = 0
+castro.do_grav = 1
+castro.do_sponge = 1
+
+castro.ppm_type = 1
+castro.grav_source_type = 2
+castro.use_pslope = 1
+castro.pslope_cutoff_density = 1.e4
+
+gravity.gravity_type = ConstantGrav
+gravity.const_grav   = -1.29e14
+
+# burning
+castro.react_rho_min = 1.e2
+castro.react_T_min = 5.e6
+
+
+# TIME STEP CONTROL
+castro.cfl            = 0.8     # cfl number for hyperbolic system
+castro.init_shrink    = 0.1     # scale back initial timestep
+castro.change_max     = 1.1     # max time step growth
+
+# SPONGE
+castro.sponge_upper_density = 1.e3
+castro.sponge_lower_density = 1.e2
+castro.sponge_timescale     = 1.0e-7
+
+# DIAGNOSTICS & VERBOSITY
+castro.sum_interval   = 1       # timesteps between computing mass
+castro.v              = 1       # verbosity in Castro.cpp
+amr.v                 = 1       # verbosity in Amr.cpp
+
+# REFINEMENT / REGRIDDING 
+amr.max_level       = 1       # maximum level number allowed
+amr.ref_ratio       = 4 2 2 2 # refinement ratio
+amr.regrid_int      = 2 2 2 2 # how often to regrid
+amr.blocking_factor = 8       # block factor in grid generation
+amr.max_grid_size   = 256
+amr.n_error_buf     = 2 2 2 2 # number of buffer cells in error est
+
+# CHECKPOINT FILES
+amr.check_file      = xrb_chk    # root name of checkpoint file
+amr.check_int       = 1000       # number of timesteps between checkpoints
+
+# PLOTFILES
+amr.plot_file       = xrb_plt    # root name of plotfile
+amr.plot_int        = 1000       # number of timesteps between plotfiles
+amr.derive_plot_vars = ALL
+
+# PROBLEM PARAMETERS
+problem.model_name =  "toy_atm_hot_castro_3cm.hse"
+
+problem.apply_perturbation = 1
+
+problem.model_shift = 300.0
+
+# MICROPHYSICS
+integrator.jacobian = 1
+
+integrator.atol_spec = 1.e-6
+integrator.rtol_spec = 1.e-6

Exec/science/xrb_layered/problem_diagnostics.H

@@ -0,0 +1,73 @@
+#ifndef problem_diagnostics_H
+#define problem_diagnostics_H
+
+AMREX_INLINE
+void
+Castro::problem_diagnostics ()
+{
+    int finest_level = parent->finestLevel();
+    Real time        = state[State_Type].curTime();
+    Real mass        = 0.0;
+    Real rho_E       = 0.0;
+    Real Tmax        = -std::numeric_limits<Real>::max();
+    Real MachMax     = -std::numeric_limits<Real>::max();
+    Real Tmax_level;
+    Real MachMax_level;
+
+    for (int lev = 0; lev <= finest_level; lev++)
+    {
+        Castro& ca_lev = getLevel(lev);
+
+        mass     += ca_lev.volWgtSum("density", time);
+
+        rho_E    += ca_lev.volWgtSum("rho_E", time);
+
+        auto temp_mf = ca_lev.derive("Temp", time, 0);
+        Tmax_level = temp_mf->max(0, 0);
+        if (Tmax_level > Tmax)
+        {
+            Tmax = Tmax_level;
+        }
+
+        auto mach_mf = ca_lev.derive("MachNumber", time, 0);
+        MachMax_level = mach_mf->max(0, 0);
+        if (MachMax_level > MachMax)
+        {
+            MachMax = MachMax_level;
+        }
+
+    }
+
+
+    if (ParallelDescriptor::IOProcessor())
+    {
+        if (verbose > 0) {
+            std::cout << '\n';
+            std::cout << "TIME= " << time << " MASS  = "   << mass  << '\n';
+            std::cout << "TIME= " << time << " RHO*E = "   << rho_E     << '\n';
+        }
+
+        std::ostream& datalog = *Castro::problem_data_logs[0];
+
+        if (time == 0.0) {
+            datalog << std::setw(14) <<  "#     time    ";
+            datalog << std::setw(14) <<  "       max(T) ";
+            datalog << std::setw(14) <<  "       max(M) ";
+            datalog << std::setw(14) <<  "        rho_E ";
+            datalog << std::endl;
+        }
+
+        // Write the quantities at this time
+        datalog << std::setw(14) <<  time;
+        datalog << std::setw(14) <<  std::setprecision(6) << Tmax;
+        datalog << std::setw(14) <<  std::setprecision(6) << MachMax;
+        datalog << std::setw(14) <<  std::setprecision(6) << rho_E;
+        datalog << std::endl;
+
+        if (verbose > 0) {
+            std::cout<<'\n';
+        }
+    }
+}
+
+#endif

Exec/science/xrb_layered/problem_initialize.H

@@ -0,0 +1,70 @@
+#ifndef problem_initialize_H
+#define problem_initialize_H
+
+#include <prob_parameters.H>
+#include <eos.H>
+#include <model_parser.H>
+#include <ambient.H>
+
+AMREX_INLINE
+void problem_initialize ()
+{
+
+    const Geometry& dgeom = DefaultGeometry();
+
+    const Real* problo = dgeom.ProbLo();
+    const Real* probhi = dgeom.ProbHi();
+
+    // Read initial model
+    read_model_file(problem::model_name);
+
+    if (ParallelDescriptor::IOProcessor()) {
+        for (int i = 0; i < model::npts; i++) {
+            std::cout << i << " " << model::profile(0).r(i) << " " << model::profile(0).state(i,model::idens) << std::endl;
+        }
+    }
+
+    // Set up Castro data logs for this problem
+
+    if (castro::sum_interval > 0 && amrex::ParallelDescriptor::IOProcessor()) {
+
+        Castro::problem_data_logs.resize(1);
+
+        Castro::problem_data_logs[0].reset(new std::fstream);
+        Castro::problem_data_logs[0]->open("xrb_diag.out", std::ios::out | std::ios::app);
+        if (!Castro::problem_data_logs[0]->good()) {
+            amrex::FileOpenFailed("xrb_diag.out");
+        }
+
+    }
+
+    // set the ambient state for the upper boundary condition
+
+    ambient::ambient_state[URHO] = model::profile(0).state(model::npts-1, model::idens);
+    ambient::ambient_state[UTEMP] = model::profile(0).state(model::npts-1, model::itemp);
+    for (int n = 0; n < NumSpec; n++) {
+        ambient::ambient_state[UFS+n] =
+            ambient::ambient_state[URHO] * model::profile(0).state(model::npts-1, model::ispec+n);
+    }
+
+    ambient::ambient_state[UMX] = 0.0_rt;
+    ambient::ambient_state[UMY] = 0.0_rt;
+    ambient::ambient_state[UMZ] = 0.0_rt;
+
+    // make the ambient state thermodynamically consistent
+
+    eos_t eos_state;
+    eos_state.rho = ambient::ambient_state[URHO];
+    eos_state.T = ambient::ambient_state[UTEMP];
+    for (int n = 0; n < NumSpec; n++) {
+        eos_state.xn[n] = ambient::ambient_state[UFS+n] / eos_state.rho;
+    }
+
+    eos(eos_input_rt, eos_state);
+
+    ambient::ambient_state[UEINT] = eos_state.rho * eos_state.e;
+    ambient::ambient_state[UEDEN] = eos_state.rho * eos_state.e;
+
+}
+#endif
+