[WIP] create a problem to study impact of species composition gradient on shock burning in 1D

Exec/science/gradient_detonation/GNUmakefile

@@ -0,0 +1,31 @@
+PRECISION  = DOUBLE
+PROFILE    = FALSE
+
+DEBUG      = FALSE
+
+DIM        = 1
+
+COMP	   = gnu
+
+USE_MPI    = TRUE
+
+USE_REACT  = TRUE
+
+USE_SHOCK_VAR = TRUE
+
+USE_SIMPLIFIED_SDC = TRUE
+
+CASTRO_HOME ?= ../../..
+
+# This sets the EOS directory in $(MICROPHYSICS_HOME)/eos
+EOS_DIR     := helmholtz
+
+# This sets the network directory in $(MICROPHYSICS_HOME)/networks
+NETWORK_DIR := subch_simple
+
+PROBLEM_DIR ?= ./
+
+Bpack   := $(PROBLEM_DIR)/Make.package
+Blocs   := $(PROBLEM_DIR)
+
+include $(CASTRO_HOME)/Exec/Make.Castro

Exec/science/gradient_detonation/Make.package

@@ -0,0 +1 @@
+

Exec/science/gradient_detonation/README.md

@@ -0,0 +1,23 @@
+# `Gradient Detonation`
+
+This is a similar setup to `Detonation`, except that we want to study
+the impact of composition gradient on shock detection by applying
+a linear interpolation function to the composition.
+
+A simple (carbon or helium) detonation.  The reaction network should
+be set through the GNUmakefile.
+
+This sets up a domain with a uniform density (dens).  A large
+temperature (``T_ign``) with length ``hotspot_size`` is placed in the accretion region with temperature ``T_l`` on the left side of the domain, and an ambient temperature (``T_r``) representing the WD core is in the right.  The transition from envelope composition to core composition occurs at the center of the domain, with the width of said region specified by ``delta``. Finally, the parameter cfrac_core, nfrac_core, and ofrac_core, are the initial C12, N14, O16 fraction in the core, and the remaining material is He4. This is the same for cfrac_env...etc, which represents the material in the envelope.
+
+The left side and right side can optionally approach each other with a
+non-zero infall velocity, vel.
+
+When run, the large temperature in the left side instantly flashes the
+fuel, generating a lot of energy and a large overpressure.  The
+signature of a propagating detonation is that a narrow energy
+generation zone will keep pace with the rightward propagating
+shockwave (as seen in the pressure field).
+
+Note: if the domain is too small, then the burning will decouple from
+the shock wave, and you will not get a detonation.

Exec/science/gradient_detonation/_prob_params

@@ -0,0 +1,59 @@
+T_l                  real          1.75e8_rt      y
+
+T_r                  real          1.e7_rt        y
+
+T_ign                real          1.75e9_rt      y
+
+hotspot_size         real          40.e5_rt       y
+
+dens                 real          1.e8_rt        y
+
+delta                real          50.e5_rt       y
+
+cfrac_env            real           0.00_rt       y
+
+cfrac_core           real          0.45_rt        y
+
+nfrac_env            real          0.0_rt         y
+
+nfrac_core           real         0.0_rt          y
+
+ofrac_env            real          0.0_rt         y
+
+ofrac_core           real          0.45_rt        y
+
+w_T                  real          5.e-4_rt       y
+
+center_T             real          0.3_rt         y
+
+smallx               real          1.e-12_rt      y
+
+vel                  real          0.0_rt         y
+
+grav_acceleration    real          0.0_rt         y
+
+idir                 integer       1              y
+
+ihe4                 integer       -1
+
+ic12                 integer       -1
+
+in14                 integer       -1
+
+io16                 integer       -1
+
+xn_env               real          0.0_rt         n     nspec
+
+xn_core              real          0.0_rt         n     nspec
+
+ambient_dens         real         0.0_rt
+
+ambient_comp         real         0.0_rt          n     nspec
+
+ambient_e_l          real         0.0_rt
+
+ambient_e_r          real         0.0_rt
+
+mu_p                 real         -6.0_rt         y
+
+mu_n                 real         -11.0_rt        y

Exec/science/gradient_detonation/inputs_grad_zingale_wd

@@ -0,0 +1,95 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 15000
+stop_time =  0.2
+
+# PROBLEM SIZE & GEOMETRY
+geometry.is_periodic = 0 0 0
+geometry.coord_sys   = 0  # 0 => cart, 1 => RZ  2=>spherical
+geometry.prob_lo     = 0     0     0
+geometry.prob_hi     = 4.e7  2500  2500
+amr.n_cell           = 400   16    16
+
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+# 0 = Interior           3 = Symmetry
+# 1 = Inflow             4 = SlipWall
+# 2 = Outflow            5 = NoSlipWall
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+castro.lo_bc       =  2   4   4
+castro.hi_bc       =  2   4   4
+
+# WHICH PHYSICS
+castro.do_hydro = 1
+castro.do_react = 1
+
+castro.ppm_type = 1
+castro.ppm_temp_fix = 0
+
+castro.transverse_reset_density = 1
+
+castro.use_flattening = 1
+
+castro.small_temp     = 1.e7
+
+castro.riemann_solver = 0
+
+castro.disable_shock_burning = 0
+
+# TIME STEP CONTROL
+castro.cfl            = 0.2     # cfl number for hyperbolic system
+castro.init_shrink    = 0.01     # scale back initial timestep
+castro.change_max     = 1.05    # scale back initial timestep
+
+castro.dtnuc_e = 0.1
+
+# DIAGNOSTICS & VERBOSITY
+castro.sum_interval   = 1       # timesteps between computing mass
+castro.v              = 1       # verbosity in Castro.cpp
+amr.v                 = 1       # verbosity in Amr.cpp
+castro.time_integration_method = 3 # use SDC by default
+#amr.grid_log        = grdlog  # name of grid logging file
+
+# REFINEMENT / REGRIDDING
+amr.max_level       = 2       # 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 = 16       # block factor in grid generation
+amr.max_grid_size   = 64
+amr.n_error_buf     = 2 2 2 2 # number of buffer cells in error est
+
+# CHECKPOINT FILES
+amr.check_file      = det_x_chk  # root name of checkpoint file
+amr.check_int       = 1000         # number of timesteps between checkpoints
+
+# PLOTFILES
+amr.plot_file       = det_x_plt  # root name of plotfile
+amr.plot_per = 1.e-3
+amr.derive_plot_vars = ALL
+
+# problem initialization
+
+problem.T_l = 2.e9
+problem.T_r = 7.5e7
+
+problem.dens = 1.12e6
+
+problem.smallx = 1.e-10
+
+problem.idir = 1
+
+problem.w_T = 1.e-3
+problem.center_T = 0.05
+
+problem.nfrac_env = 0.01
+
+# refinement
+
+amr.refinement_indicators = temperr tempgrad
+
+amr.refine.temperr.max_level = 5
+amr.refine.temperr.value_greater = 4.e9
+amr.refine.temperr.field_name = Temp
+
+amr.refine.tempgrad.max_level = 5
+amr.refine.tempgrad.gradient = 1.e8
+amr.refine.tempgrad.field_name = Temp
+

Exec/science/gradient_detonation/problem_bc_fill.H

@@ -0,0 +1,115 @@
+#ifndef problem_bc_fill_H
+#define problem_bc_fill_H
+
+// Given a zone state, fill it with ambient material.
+
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void problem_fill_ambient(int i, int j, int k,
+                          Array4<Real> const& state,
+                          Real x, Real time,
+                          const GeometryData& geomdata)
+{
+    const Real* problo = geomdata.ProbLo();
+    const Real* probhi = geomdata.ProbHi();
+
+    Real c_T = problo[0] + problem::center_T * (probhi[0] - problo[0]);
+
+    state(i,j,k,URHO) = problem::ambient_dens;
+
+    for (int n = 0; n < NumSpec; ++n) {
+        state(i,j,k,UFS+n) = problem::ambient_dens * problem::ambient_comp[n];
+    }
+
+    if (x < c_T) {
+        state(i,j,k,UTEMP) = problem::T_l;
+        state(i,j,k,UEINT) = state(i,j,k,URHO) * problem::ambient_e_l;
+        state(i,j,k,UMX  ) = state(i,j,k,URHO) * (problem::vel + problem::grav_acceleration * time);
+    }
+    else {
+        state(i,j,k,UTEMP) = problem::T_r;
+        state(i,j,k,UEINT) = state(i,j,k,URHO) * problem::ambient_e_r;
+        state(i,j,k,UMX  ) = -state(i,j,k,URHO) * (problem::vel + problem::grav_acceleration * time);
+    }
+
+    state(i,j,k,UMY) = 0.0;
+    state(i,j,k,UMZ) = 0.0;
+    state(i,j,k,UEDEN) = state(i,j,k,UEINT) + 0.5_rt * (state(i,j,k,UMX) * state(i,j,k,UMX) +
+                                                        state(i,j,k,UMY) * state(i,j,k,UMY) +
+                                                        state(i,j,k,UMZ) * state(i,j,k,UMZ)) /
+                                                       state(i,j,k,URHO);
+}
+
+
+
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void problem_bc_fill(int i, int j, int k,
+                     Array4<Real> const& state,
+                     Real time,
+                     const Array1D<BCRec, 0, NUM_STATE-1>& bcs,
+                     const GeometryData& geomdata)
+{
+    // Override the generic routine at the physical boundaries by
+    // setting the material to the ambient state. Note that we
+    // don't want to do this for interior/periodic boundaries
+    // or for reflecting boundaries.
+
+    if (castro::fill_ambient_bc == 1) {
+
+        const int* domlo = geomdata.Domain().loVect();
+        const int* domhi = geomdata.Domain().hiVect();
+
+        bool do_ambient_fill = false;
+
+        // Note: in these checks we're only looking at the boundary
+        // conditions on the first state component (density).
+
+        if (i < domlo[0] && (bcs(0).lo(0) != amrex::BCType::reflect_odd &&
+                             bcs(0).lo(0) != amrex::BCType::int_dir &&
+                             bcs(0).lo(0) != amrex::BCType::reflect_even)) {
+            do_ambient_fill = true;
+        }
+        else if (i > domhi[0] && (bcs(0).hi(0) != amrex::BCType::reflect_odd &&
+                                  bcs(0).hi(0) != amrex::BCType::int_dir &&
+                                  bcs(0).hi(0) != amrex::BCType::reflect_even)) {
+            do_ambient_fill = true;
+        }
+
+        if (AMREX_SPACEDIM >= 2) {
+            if (j < domlo[1] && (bcs(0).lo(1) != amrex::BCType::reflect_odd &&
+                                 bcs(0).lo(1) != amrex::BCType::int_dir &&
+                                 bcs(0).lo(1) != amrex::BCType::reflect_even)) {
+                do_ambient_fill = true;
+            }
+            else if (j > domhi[1] && (bcs(0).hi(1) != amrex::BCType::reflect_odd &&
+                                      bcs(0).hi(1) != amrex::BCType::int_dir &&
+                                      bcs(0).hi(1) != amrex::BCType::reflect_even)) {
+                do_ambient_fill = true;
+            }
+        }
+
+        if (AMREX_SPACEDIM == 3) {
+            if (k < domlo[2] && (bcs(0).lo(2) != amrex::BCType::reflect_odd &&
+                                 bcs(0).lo(2) != amrex::BCType::int_dir &&
+                                 bcs(0).lo(2) != amrex::BCType::reflect_even)) {
+                do_ambient_fill = true;
+            }
+            else if (k > domhi[2] && (bcs(0).hi(2) != amrex::BCType::reflect_odd &&
+                                      bcs(0).hi(2) != amrex::BCType::int_dir &&
+                                      bcs(0).hi(2) != amrex::BCType::reflect_even)) {
+                do_ambient_fill = true;
+            }
+        }
+
+        if (do_ambient_fill) {
+            const Real* dx = geomdata.CellSize();
+            const Real* problo = geomdata.ProbLo();
+
+            Real x = problo[0] + (static_cast<Real>(i) + 0.5_rt) * dx[0];
+
+            problem_fill_ambient(i, j, k, state, x, time, geomdata);
+        }
+
+    }
+}
+
+#endif