add the simple_convection problem

Exec/hydro_tests/simple_convection/GNUmakefile

@@ -0,0 +1,31 @@
+PRECISION        = DOUBLE
+PROFILE          = FALSE
+DEBUG            = FALSE
+DIM              = 2
+
+COMP	         = gnu
+
+USE_MPI          = TRUE
+USE_OMP          = FALSE
+
+USE_GRAV         = TRUE
+USE_REACT        = FALSE
+
+USE_MODEL_PARSER = TRUE
+USE_RNG_STATE_INIT = TRUE
+
+CASTRO_HOME ?= ../../..
+
+# This sets the EOS directory in $(MICROPHYSICS_HOME)/EOS
+EOS_DIR     := gamma_law
+
+# This sets the network directory in $(MICROPHYSICS_HOME)/networks
+NETWORK_DIR := general_null
+NETWORK_INPUTS = gammalaw.net
+
+PROBLEM_DIR ?= ./
+
+Bpack   := $(PROBLEM_DIR)/Make.package
+Blocs   := $(PROBLEM_DIR)
+
+include $(CASTRO_HOME)/Exec/Make.Castro

Exec/hydro_tests/simple_convection/Make.package

@@ -0,0 +1,3 @@
+CEXE_headers += initial_model.H
+
+

Exec/hydro_tests/simple_convection/README.md

@@ -0,0 +1,5 @@
+# simple_convection
+
+This is an implementation of the convection problem from pyro.
+
+

Exec/hydro_tests/simple_convection/_prob_params

@@ -0,0 +1,12 @@
+dens_base      real     10.0    y
+
+scale_height   real      4.0    y
+
+y_height       real      2.0    y
+
+thickness      real      0.25   y
+
+e_rate         real      0.1    y
+
+low_density_cutoff   real    0.01    y
+

Exec/hydro_tests/simple_convection/initial_model.H

@@ -0,0 +1,94 @@
+#ifndef INITIAL_MODEL_H
+#define INITIAL_MODEL_H
+
+#include <prob_parameters.H>
+#include <network.H>
+#include <eos.H>
+
+struct model_t {
+    amrex::Real dens_base = -1;
+    amrex::Real scale_height = -1;
+    Real xn[NumSpec] = {0.0};
+};
+
+
+
+///
+/// construct an initial model in HSE.  Note: this does not return
+/// anything, but rather updates the model_parser globals with the
+/// model information.
+///
+AMREX_INLINE
+void
+generate_initial_model(const int npts_model, const Real xmin, const Real xmax,
+                       const model_t model_params) {
+
+    model::npts = npts_model;
+    model::initialized = true;
+
+    if (npts_model > NPTS_MODEL) {
+        amrex::Error("Error: model has more than NPTS_MODEL points,  Increase MAX_NPTS_MODEL");
+    }
+
+    // compute the pressure scale height (for an isothermal, ideal-gas
+    // atmosphere)
+
+    amrex::Real pres_base = model_params.scale_height * model_params.dens_base * std::abs(gravity::const_grav);
+
+    // create the grid -- cell centers
+
+    Real dx = (xmax - xmin) / npts_model;
+
+    for (int i = 0; i < npts_model; i++) {
+        model::profile(0).r(i) = xmin + (static_cast<Real>(i) + 0.5_rt) * dx;
+    }
+
+    for (int i = 0; i < npts_model; i++) {
+
+        amrex::Real density;
+        amrex::Real pressure;
+
+        amrex::Real profile = 1.0 - (eos_rp::eos_gamma - 1.0) / eos_rp::eos_gamma * model::profile(0).r(i) / model_params.scale_height;
+
+        if (profile > 0.0) {
+            density = std::max(model_params.dens_base * std::pow(profile, 1.0/(eos_rp::eos_gamma - 1.0)),
+                               problem::low_density_cutoff);
+        } else {
+            density = problem::low_density_cutoff;
+        }
+
+        if (i == 0) {
+            pressure = pres_base;
+        } else if (density <= problem::low_density_cutoff + 1.e-30_rt) {
+            pressure = model::profile(0).state(i-1, model::ipres);
+        } else {
+            pressure = pres_base *
+                std::pow(density / model_params.dens_base, eos_rp::eos_gamma);
+        }
+
+        // initial guess
+
+        amrex::Real temp = T_guess;
+
+        eos_t eos_state;
+        eos_state.p = pressure;
+        eos_state.T = temp;
+        eos_state.rho = density;
+        for (int n = 0; n < NumSpec; n++) {
+            eos_state.xn[n] = model_params.xn[n];
+        }
+
+        eos(eos_input_rp, eos_state);
+
+        model::profile(0).state(i, model::idens) = density;
+        model::profile(0).state(i, model::ipres) = pressure;
+        model::profile(0).state(i, model::itemp) = eos_state.T;
+        for (int n = 0; n < NumSpec; n++) {
+            model::profile(0).state(i, model::ispec+n) = model_params.xn[n];
+        }
+
+    }
+
+}
+
+#endif

Exec/hydro_tests/simple_convection/inputs_2d

@@ -0,0 +1,78 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 1000000
+stop_time =  10.0
+
+# PROBLEM SIZE & GEOMETRY
+geometry.is_periodic = 0       0
+geometry.coord_sys   = 0                  # 0 => cart, 1 => RZ  2=>spherical
+geometry.prob_lo     = 0       0
+geometry.prob_hi     = 4.0     12.0
+amr.n_cell           = 128     384
+
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+# 0 = Interior           3 = Symmetry
+# 1 = Inflow             4 = SlipWall
+# 2 = Outflow            5 = NoSlipWall
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+castro.lo_bc       =  2   3
+castro.hi_bc       =  2   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 = 0
+castro.add_ext_src = 1
+castro.do_grav = 1
+castro.do_sponge = 1
+
+castro.ppm_type = 1
+castro.grav_source_type = 2
+castro.use_pslope = 1
+
+gravity.gravity_type = ConstantGrav
+gravity.const_grav   = -2.0
+
+# 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 = 0.02
+castro.sponge_lower_density = 0.002
+castro.sponge_timescale     = 1.e-2
+
+# 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       = 0       # maximum level number allowed
+amr.ref_ratio       = 2 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   = 128
+amr.n_error_buf     = 2 2 2 2 # number of buffer cells in error est
+
+# CHECKPOINT FILES
+amr.check_file      = chk        # root name of checkpoint file
+amr.check_int       = 1000       # number of timesteps between checkpoints
+
+# PLOTFILES
+amr.plot_file        = plt        # root name of plotfile
+amr.plot_int         = 1000       # number of timesteps between plotfiles
+amr.derive_plot_vars = ALL
+
+# PROBLEM PARAMETERS
+problem.scale_height = 2.0
+problem.dens_base = 1000.0
+problem.low_density_cutoff = 1.e-3
+
+problem.e_rate = 0.5
+
+

Exec/hydro_tests/simple_convection/inputs_2d.big

@@ -0,0 +1,78 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 1000000
+stop_time =  10.0
+
+# PROBLEM SIZE & GEOMETRY
+geometry.is_periodic = 0       0
+geometry.coord_sys   = 0                  # 0 => cart, 1 => RZ  2=>spherical
+geometry.prob_lo     = 0       0
+geometry.prob_hi     = 8.0     12.0
+amr.n_cell           = 512     768
+
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+# 0 = Interior           3 = Symmetry
+# 1 = Inflow             4 = SlipWall
+# 2 = Outflow            5 = NoSlipWall
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+castro.lo_bc       =  2   3
+castro.hi_bc       =  2   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 = 0
+castro.add_ext_src = 1
+castro.do_grav = 1
+castro.do_sponge = 1
+
+castro.ppm_type = 1
+castro.grav_source_type = 2
+castro.use_pslope = 1
+
+gravity.gravity_type = ConstantGrav
+gravity.const_grav   = -2.0
+
+# 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 = 0.02
+castro.sponge_lower_density = 0.002
+castro.sponge_timescale     = 1.e-2
+
+# 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       = 0       # maximum level number allowed
+amr.ref_ratio       = 2 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   = 128
+amr.n_error_buf     = 2 2 2 2 # number of buffer cells in error est
+
+# CHECKPOINT FILES
+amr.check_file      = chk        # root name of checkpoint file
+amr.check_int       = 1000       # number of timesteps between checkpoints
+
+# PLOTFILES
+amr.plot_file        = plt        # root name of plotfile
+amr.plot_int         = 1000       # number of timesteps between plotfiles
+amr.derive_plot_vars = ALL
+
+# PROBLEM PARAMETERS
+problem.scale_height = 2.0
+problem.dens_base = 1000.0
+problem.low_density_cutoff = 1.e-3
+
+problem.e_rate = 0.5
+
+

Exec/hydro_tests/simple_convection/problem_initialize.H

@@ -0,0 +1,81 @@
+#ifndef problem_initialize_H
+#define problem_initialize_H
+
+#include <prob_parameters.H>
+#include <eos.H>
+#include <model_parser.H>
+#include <initial_model.H>
+#include <global.H>
+#include <ambient.H>
+
+AMREX_INLINE
+void problem_initialize ()
+{
+
+    const Geometry& dgeom = DefaultGeometry();
+
+    const Real* problo = dgeom.ProbLo();
+    const Real* probhi = dgeom.ProbHi();
+
+    problem::center[0] = 0.5_rt * (problo[0] + probhi[0]);
+#if AMREX_SPACEDIM >= 2
+    problem::center[1] = 0.5_rt * (problo[1] + probhi[1]);
+#endif
+#if AMREX_SPACEDIM == 3
+    problem::center[2] = 0.5_rt * (problo[2] + probhi[2]);
+#endif
+
+    // first make a 1D initial model for the entire domain
+    // we'll create it ourselves, but hook it into the model_parser
+
+    // we use the fine grid dx for the model resolution
+    auto fine_geom = global::the_amr_ptr->Geom(global::the_amr_ptr->maxLevel());
+
+    auto dx = fine_geom.CellSizeArray();
+    auto dx_model = dx[AMREX_SPACEDIM-1];
+
+    int nx = (2.0_rt * problem::center[AMREX_SPACEDIM-1] + 1.e-8_rt) /
+        dx_model;
+
+    model_t model_params;
+    model_params.dens_base = problem::dens_base;
+    model_params.scale_height = problem::scale_height;
+    for (int n = 0; n < NumSpec; ++n) {
+        model_params.xn[n] = 0.0_rt;
+    }
+    // only initialize the first species
+    model_params.xn[0] = 1.0_rt;
+
+    generate_initial_model(nx, problo[AMREX_SPACEDIM-1], probhi[AMREX_SPACEDIM-1],
+                           model_params);
+
+
+    // set the ambient state for the upper BC
+    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