Create problem to study He flame movement in a double det (#2870)

Begin creating a problem to study the lateral movement of a He flame in the accreted layer of a white dwarf.

We want to determine how various shock detection thresholds and burning treatments affect the initial flame movement in a double detonation model. In particular, we seek to see if the results of disabling/enabling burning in detected shock zones converges as resolution becomes finer.

Exec/science/subch_planar/GNUmakefile

@@ -0,0 +1,33 @@
+PRECISION        = DOUBLE
+PROFILE          = FALSE
+DEBUG            = FALSE
+DIM              = 2
+
+COMP	         = gnu
+
+USE_MPI          = TRUE
+USE_OMP          = FALSE
+
+USE_GRAV         = TRUE
+USE_REACT        = TRUE
+
+USE_SHOCK_VAR = TRUE
+
+USE_MODEL_PARSER = 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/subch_planar/Make.package

@@ -0,0 +1 @@
+

Exec/science/subch_planar/README.md

@@ -0,0 +1 @@
+We want to study the impact of shock detection and burning treatment in the context of a double detonation SNe Ia by simulating a lateral flame in the accreted He layer on the white dwarf.

Exec/science/subch_planar/_prob_params

@@ -0,0 +1,17 @@
+H_min                 real           1.e-4_rt      y
+
+model_name            string         ""            y
+
+# the distance from the base of the He layer star to put the perturbation
+R_pert                real           1.e7_rt       y
+
+# location of base of He layer
+R_He_base             real           0.0_rt
+
+# the amplitude of the temperature perturbation (above background)
+pert_temp_factor      real           10.0_rt       y
+
+# the physical scale of the perturbation
+pert_rad_factor       real            2.0_rt       y
+
+ihe4                  integer         -1

Exec/science/subch_planar/inputs_2d

@@ -0,0 +1,86 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 1000000
+stop_time =  1.0
+
+# PROBLEM SIZE & GEOMETRY
+geometry.is_periodic = 0       0
+geometry.coord_sys   = 1                  # 0 => cart, 1 => RZ  2=>spherical
+geometry.prob_lo     = 0.0     0.0
+geometry.prob_hi     = 4.e8    4.e8
+amr.n_cell           = 200     200
+
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+# 0 = Interior           3 = Symmetry
+# 1 = Inflow             4 = SlipWall
+# 2 = Outflow            5 = NoSlipWall
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+castro.lo_bc       =  3   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 = 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.466e9
+
+# TIME STEP CONTROL
+castro.cfl            = 0.7     # 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 = 5.0e4
+castro.sponge_lower_density = 8.0e1
+castro.sponge_timescale     = 1.0e-3
+
+# 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   = 64 
+amr.n_error_buf     = 2 2 2 2 # number of buffer cells in error est
+
+# CHECKPOINT FILES
+amr.check_file      = planar_chk    # root name of checkpoint file
+amr.check_int       = 200       # number of timesteps between checkpoints
+
+# PLOTFILES
+amr.plot_file       = planar_plt    # root name of plotfile
+amr.plot_int        = -1       # number of timesteps between plotfiles
+amr.plot_per        = 0.002   #simulation time (s) between plotfiles
+amr.derive_plot_vars = ALL
+
+# PROBLEM PARAMETERS
+problem.model_name =  "toy_subch.hse.tanh.delta_50.00km.dx_20.00km"
+
+problem.pert_temp_factor = 20.0
+problem.pert_rad_factor = 0.5
+problem.R_pert = 1.e7
+
+# MICROPHYSICS
+integrator.jacobian = 1
+
+integrator.use_burn_retry = 1
+integrator.retry_swap_jacobian = 1
+
+integrator.atol_spec = 1.e-6
+integrator.rtol_spec = 1.e-6

Exec/science/subch_planar/problem_initialize.H

@@ -0,0 +1,73 @@
+#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();
+
+    problem::ihe4 = network_spec_index("helium-4");
+
+    if (problem::ihe4 < 0) {
+        amrex::Error("Error: no He4 in the network");
+    }
+
+    // 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 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;
+
+// find the distance (starting from the center) where the He layer begins
+
+    problem::R_He_base = 0.0;
+
+    for (int n = 0; n < model::npts; n++) {
+        if (model::profile(0).state(n, model::ispec+problem::ihe4) > 0.5_rt) {
+            problem::R_He_base = model::profile(0).r(n);
+            break;
+        }
+    }
+
+    amrex::Print() << "base of He layer found at " << problem::R_He_base << std::endl;
+}
+#endif

Exec/science/subch_planar/problem_initialize_state_data.H

@@ -0,0 +1,126 @@
+#ifndef problem_initialize_state_data_H
+#define problem_initialize_state_data_H
+
+#include <prob_parameters.H>
+#include <eos.H>
+#include <model_parser.H>
+
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void problem_initialize_state_data (int i, int j, int k,
+                                    Array4<Real> const& state,
+                                    const GeometryData& geomdata)
+{
+
+    const Real* dx = geomdata.CellSize();
+    const Real* problo = geomdata.ProbLo();
+
+    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
+
+    state(i,j,k,URHO) = interpolate(height, model::idens);
+    state(i,j,k,UTEMP) = interpolate(height, model::itemp);
+    for (int n = 0; n < NumSpec; n++) {
+        state(i,j,k,UFS+n) = interpolate(height, model::ispec+n);
+    }
+
+    eos_t eos_state;
+    eos_state.rho = state(i,j,k,URHO);
+    eos_state.T = state(i,j,k,UTEMP);
+    for (int n = 0; n < NumSpec; n++) {
+        eos_state.xn[n] = state(i,j,k,UFS+n);
+    }
+
+    eos(eos_input_rt, eos_state);
+
+    state(i,j,k,UEINT) = state(i,j,k,URHO) * eos_state.e;
+    state(i,j,k,UEDEN) = state(i,j,k,UEINT);
+
+    Real sumX{0.0_rt};
+
+    for (int n = 0; n < NumSpec; n++) {
+        state(i,j,k,UFS+n) = state(i,j,k,URHO) * state(i,j,k,UFS+n);
+        sumX += state(i,j,k,UFS+n);
+    }
+
+    // normalize
+    for (int n = 0; n < NumSpec; ++n) {
+        state(i,j,k,UFS+n) /= sumX;
+    }
+
+    burn_t burn_state;
+
+    burn_state.rho = state(i,j,k,URHO);
+    burn_state.T = state(i,j,k,UTEMP);
+    for (int n = 0; n < NumSpec; n++) {
+        burn_state.xn[n] = state(i,j,k,UFS+n);
+    }
+
+    eos(eos_input_rt, burn_state);
+
+    state(i,j,k,UEINT) = state(i,j,k,URHO) * burn_state.e;
+    state(i,j,k,UEDEN) = state(i,j,k,URHO) * burn_state.e;
+
+    for (int n = 0; n < NumSpec; n++) {
+        state(i,j,k,UFS+n) = state(i,j,k,URHO) * state(i,j,k,UFS+n);
+    }
+
+    // Initial velocities = 0
+    state(i,j,k,UMX) = 0.0_rt;
+    state(i,j,k,UMY) = 0.0_rt;
+    state(i,j,k,UMZ) = 0.0_rt;
+
+// add a perturbation at the north pole
+
+    Real T0 = state(i,j,k,UTEMP);
+
+    // perturbation is on the vertical-axis
+
+    Real pert_center = problem::R_pert + problem::R_He_base;
+
+    #if AMREX_SPACEDIM == 1
+    Real r1 = std::sqrt((x - pert_center) * (x - pert_center)) /
+        (2.5e6_rt * problem::pert_rad_factor);
+    #elif AMREX_SPACEDIM == 2
+    Real r1 = std::sqrt(x * x + (y - pert_center) * (y - pert_center)) /
+        (2.5e6_rt * problem::pert_rad_factor);
+    #else
+    Real r1 = std::sqrt(x * x + y * y + (z - pert_center) * (z - pert_center)) /
+        (2.5e6_rt * problem::pert_rad_factor);
+    #endif
+
+    // convolve the temperature perturbation with the amount of He
+    Real X_he = burn_state.xn[problem::ihe4];
+
+    Real Tpert = T0 * (1.0_rt + X_he * problem::pert_temp_factor *
+                       (0.150e0_rt * (1.0_rt + std::tanh(2.0_rt - r1))));
+
+    Real dT = Tpert - T0;
+
+    burn_state.rho = state(i,j,k,URHO);
+    burn_state.T = T0 + dT;
+
+    // we don't need to refill xn, since it still holds unchanged from above
+
+    eos(eos_input_rt, burn_state);
+
+    // the internal energy changed
+
+    state(i,j,k,UEINT) = burn_state.e * state(i,j,k,URHO);
+    state(i,j,k,UEDEN) = burn_state.e * state(i,j,k,URHO);
+}
+#endif