From e76e528905cfc00e0953d5fe241f65ae73e84393 Mon Sep 17 00:00:00 2001
From: Ryan Brady <ryan.m.brady@stonybrook.edu>
Date: Wed, 31 Jul 2024 12:54:38 -0400
Subject: [PATCH 01/14] initial commit

---
 Exec/science/gradient_detonation/GNUmakefile  |  32 +++++
 Exec/science/gradient_detonation/Make.package |   1 +
 Exec/science/gradient_detonation/README.md    |  29 +++++
 Exec/science/gradient_detonation/_prob_params |  51 ++++++++
 .../gradient_detonation/problem_bc_fill.H     | 115 ++++++++++++++++
 .../gradient_detonation/problem_initialize.H  | 123 ++++++++++++++++++
 .../problem_initialize_state_data.H           |  51 ++++++++
 .../gradient_detonation/problem_source.H      |  39 ++++++
 8 files changed, 441 insertions(+)
 create mode 100644 Exec/science/gradient_detonation/GNUmakefile
 create mode 100644 Exec/science/gradient_detonation/Make.package
 create mode 100644 Exec/science/gradient_detonation/README.md
 create mode 100644 Exec/science/gradient_detonation/_prob_params
 create mode 100644 Exec/science/gradient_detonation/problem_bc_fill.H
 create mode 100644 Exec/science/gradient_detonation/problem_initialize.H
 create mode 100644 Exec/science/gradient_detonation/problem_initialize_state_data.H
 create mode 100644 Exec/science/gradient_detonation/problem_source.H

diff --git a/Exec/science/gradient_detonation/GNUmakefile b/Exec/science/gradient_detonation/GNUmakefile
new file mode 100644
index 0000000000..9cfb52c4d2
--- /dev/null
+++ b/Exec/science/gradient_detonation/GNUmakefile
@@ -0,0 +1,32 @@
+PRECISION  = DOUBLE
+PROFILE    = FALSE
+
+DEBUG      = FALSE
+
+DIM        = 1
+
+COMP	   = gnu
+
+USE_MPI    = TRUE
+
+USE_REACT  = TRUE
+
+CASTRO_HOME ?= ../../..
+
+# This sets the EOS directory in $(MICROPHYSICS_HOME)/eos
+EOS_DIR     := helmholtz
+
+# This sets the network directory in $(MICROPHYSICS_HOME)/networks
+ifeq ($(USE_NSE_NET), TRUE)
+	NETWORK_DIR := subch_base
+	SCREEN_METHOD := chabrier1998
+else
+	NETWORK_DIR := aprox19
+endif
+
+PROBLEM_DIR ?= ./
+
+Bpack   := $(PROBLEM_DIR)/Make.package
+Blocs   := $(PROBLEM_DIR)
+
+include $(CASTRO_HOME)/Exec/Make.Castro
diff --git a/Exec/science/gradient_detonation/Make.package b/Exec/science/gradient_detonation/Make.package
new file mode 100644
index 0000000000..8b13789179
--- /dev/null
+++ b/Exec/science/gradient_detonation/Make.package
@@ -0,0 +1 @@
+
diff --git a/Exec/science/gradient_detonation/README.md b/Exec/science/gradient_detonation/README.md
new file mode 100644
index 0000000000..f956ab9290
--- /dev/null
+++ b/Exec/science/gradient_detonation/README.md
@@ -0,0 +1,29 @@
+# `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 sinusoidal 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_l) is placed in the left side of the domain, and an
+ambient temperature (T_r) is in the right.  The parameter center_T
+specifies where, as a fraction of the domain length, to put the
+interface, while width_T determines how wide the transition region is,
+as a fraction of the domain length. Finally, the parameter cfrac, nfrac,
+and ofrac, are the initial C12, N14, O16 fraction, and the remaining
+material is He4.
+
+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.
diff --git a/Exec/science/gradient_detonation/_prob_params b/Exec/science/gradient_detonation/_prob_params
new file mode 100644
index 0000000000..0f94589a81
--- /dev/null
+++ b/Exec/science/gradient_detonation/_prob_params
@@ -0,0 +1,51 @@
+T_l                  real          1.e9_rt        y
+
+T_r                  real          5.e7_rt        y
+
+dens                 real          1.e8_rt        y
+
+cfrac_env           real           0.00_rt         y
+
+cfrac_core           real          0.45_rt         y
+
+nfrac_env            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
\ No newline at end of file
diff --git a/Exec/science/gradient_detonation/problem_bc_fill.H b/Exec/science/gradient_detonation/problem_bc_fill.H
new file mode 100644
index 0000000000..d16f07bf0d
--- /dev/null
+++ b/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
diff --git a/Exec/science/gradient_detonation/problem_initialize.H b/Exec/science/gradient_detonation/problem_initialize.H
new file mode 100644
index 0000000000..1f4165f27e
--- /dev/null
+++ b/Exec/science/gradient_detonation/problem_initialize.H
@@ -0,0 +1,123 @@
+#ifndef problem_initialize_H
+#define problem_initialize_H
+
+#include <prob_parameters.H>
+#include <eos.H>
+#include <network.H>
+
+AMREX_INLINE
+void problem_initialize ()
+{
+
+    // get the species indices
+
+    problem::ihe4 = network_spec_index("helium-4");
+    problem::ic12 = network_spec_index("carbon-12");
+    problem::in14 = network_spec_index("nitrogen-14");
+    problem::io16 = network_spec_index("oxygen-16");
+
+    if (problem::ihe4 < 0 || problem::ic12 < 0 || problem::io16 < 0) {
+        amrex::Error("ERROR: species indices not found");
+    }
+
+    // make sure that the carbon fraction falls between 0 and 1
+    // in the core and the envelope
+
+    if (problem::cfrac_core > 1.0_rt || problem::cfrac_core < 0.0_rt
+        || problem::cfrac_env > 1.0_rt || problem::cfrac_env < 0.0_rt) {
+        amrex::Error("ERROR: cfrac must fall between 0 and 1");
+    }
+
+    // make sure that the nitrogen fraction falls between 0 and 1
+    // in the core and the envelope
+
+    if (problem::nfrac_core > 1.0_rt || problem::nfrac_core < 0.0_rt
+        || problem::nfrac_env > 1.0_rt || problem::nfrac_env < 0.0_rt) {
+        amrex::Error("ERROR: nfrac must fall between 0 and 1");
+    }
+
+    // make sure that the oxygen fraction falls between 0 and 1
+    // in the core and the envelope
+
+    if (problem::ofrac_core > 1.0_rt || problem::ofrac_core < 0.0_rt
+        || problem::ofrac_env > 1.0_rt || problem::ofrac_env < 0.0_rt) {
+        amrex::Error("ERROR: ofrac must fall between 0 and 1");
+    }
+
+    // make sure that the C/O fraction sums to no more than 1
+    // in the core and the envelope
+
+    if (problem::cfrac_core + problem::nfrac_core + problem::ofrac_core > 1.0_rt
+        || problem::cfrac_env + problem::nfrac_env + problem::ofrac_env > 1.0_rt) {
+        amrex::Error("ERROR: cfrac + nfrac + ofrac cannot exceed 1.");
+    }
+
+    // set the default mass fractions
+
+    for (int n = 0; n < NumSpec; n++) {
+        problem::xn_core[n] = problem::smallx;
+        problem::xn_env[n] = problem::smallx;
+    }
+
+    problem::xn_core[problem::ic12] = amrex::max(problem::cfrac_core, problem::smallx);
+    problem::xn_core[problem::io16] = amrex::max(problem::ofrac_core, problem::smallx);
+
+    problem::xn_env[problem::ic12] = amrex::max(problem::cfrac_env, problem::smallx);
+    problem::xn_env[problem::io16] = amrex::max(problem::ofrac_env, problem::smallx);
+
+    if (problem::in14 >= 0) {
+       problem::xn_core[problem::in14] = amrex::max(problem::nfrac_core, problem::smallx);
+       problem::xn_core[problem::ihe4] = 1.0_rt - problem::xn_core[problem::ic12]
+                                           - problem::xn_core[problem::in14]
+                                           - problem::xn_core[problem::io16]
+                                           - (NumSpec - 4) * problem::smallx;
+
+       problem::xn_env[problem::in14] = amrex::max(problem::nfrac_env, problem::smallx);
+       problem::xn_env[problem::ihe4] = 1.0_rt - problem::xn_env[problem::ic12]
+                                           - problem::xn_env[problem::in14]
+                                           - problem::xn_env[problem::io16]
+                                           - (NumSpec - 4) * problem::smallx;
+    }
+    else {
+       problem::xn_core[problem::ihe4] = 1.0_rt - problem::xn_core[problem::ic12]
+                                           - problem::xn_core[problem::io16]
+                                           - (NumSpec - 3) * problem::smallx;
+
+       problem::xn_env[problem::ihe4] = 1.0_rt - problem::xn_env[problem::ic12]
+                                           - problem::xn_env[problem::io16]
+                                           - (NumSpec - 3) * problem::smallx;
+    }
+
+    // Set the ambient material
+
+    problem::ambient_dens = problem::dens;
+    for (int n = 0; n < NumSpec; n++) {
+        problem::ambient_comp[n] = problem::xn_core[n] + problem::xn_env[n];
+    }
+
+    eos_t eos_state;
+    eos_state.rho = problem::ambient_dens;
+    for (int n = 0; n < NumSpec; n++) {
+        eos_state.xn[n]  = problem::ambient_comp[n];
+    }
+
+#ifdef AUX_THERMO
+    // set the aux quantities -- we need to do this if we are using the NSE network
+    set_aux_comp_from_X(eos_state);
+#endif
+
+    eos_state.T = problem::T_l;
+
+    eos(eos_input_rt, eos_state);
+
+    problem::ambient_e_l = eos_state.e;
+
+    eos_state.T = problem::T_r;
+
+    eos(eos_input_rt, eos_state);
+
+    problem::ambient_e_r = eos_state.e;
+
+}
+
+#endif
diff --git a/Exec/science/gradient_detonation/problem_initialize_state_data.H b/Exec/science/gradient_detonation/problem_initialize_state_data.H
new file mode 100644
index 0000000000..76d34c3daa
--- /dev/null
+++ b/Exec/science/gradient_detonation/problem_initialize_state_data.H
@@ -0,0 +1,51 @@
+#ifndef problem_initialize_state_data_H
+#define problem_initialize_state_data_H
+
+#include <prob_parameters.H>
+#include <eos.H>
+#include <network.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();
+    const Real* probhi = geomdata.ProbHi();
+
+    Real width = problem::w_T * (probhi[0] - problo[0]);
+    Real c_T = problo[0] + problem::center_T * (probhi[0] - problo[0]);
+
+    Real xcen = problo[0] + dx[0] * (static_cast<Real>(i) + 0.5_rt);
+
+    state(i,j,k,URHO) = problem::dens;
+
+    Real sigma = 1.0_rt / (1.0_rt + std::exp(-(c_T - xcen)/ width));
+    Real freq = 1.0_rt * 2.0 * 3.1415926535897932384
+
+    state(i,j,k,UTEMP) = problem::T_l + (problem::T_r - problem::T_l) * (1.0_rt - sigma);
+
+    for (int n = 0; n < NumSpec; n++) {
+        state(i,j,k,UFS+n) = state(i,j,k,URHO) * (problem::xn_envelope[n] * std::sin(freq * i) + problem::xn_core[n] * (1.0 - std::sin(freq * i)));
+    }
+
+    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] = problem::xn_core[n] + xn_env[n];
+    }
+
+    eos(eos_input_rt, eos_state);
+
+    state(i,j,k,UMX) = state(i,j,k,URHO) * (problem::vel - 2 * problem::vel * (1.0_rt - sigma));
+    state(i,j,k,UMY) = 0.0_rt;
+    state(i,j,k,UMZ) = 0.0_rt;
+    state(i,j,k,UEINT) = state(i,j,k,URHO) * eos_state.e;
+    state(i,j,k,UEDEN) = state(i,j,k,UEINT) + 0.5_rt * state(i,j,k,URHO) * problem::vel * problem::vel;
+
+}
+
+#endif
diff --git a/Exec/science/gradient_detonation/problem_source.H b/Exec/science/gradient_detonation/problem_source.H
new file mode 100644
index 0000000000..ef3b4572c5
--- /dev/null
+++ b/Exec/science/gradient_detonation/problem_source.H
@@ -0,0 +1,39 @@
+#ifndef problem_source_H
+#define problem_source_H
+
+#include <prob_parameters.H>
+
+using namespace amrex;
+
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void problem_source (int i, int j, int k,
+                     GeometryData const& geomdata,
+                     Array4<const Real> const& state,
+                     Array4<Real> const& src,
+                     const Real dt, const Real time)
+{
+
+    amrex::ignore_unused(dt);
+    amrex::ignore_unused(time);
+
+    // Add a mock gravitational acceleration which points to the center
+    // with uniform magnitude on either side of the center.
+
+    const Real* problo = geomdata.ProbLo();
+    const Real* probhi = geomdata.ProbHi();
+    const Real* dx     = geomdata.CellSize();
+
+    auto c_T = problo[0] + problem::center_T * (probhi[0] - problo[0]);
+
+    auto x = problo[0] + (Real(i) + 0.5_rt) * dx[0] - problem::center[0];
+
+    auto g = (x > c_T) ? -problem::grav_acceleration : problem::grav_acceleration;
+
+    src(i,j,k,UMX) = state(i,j,k,URHO) * g;
+
+    // Energy source is v . momentum source
+
+    src(i,j,k,UEDEN) = state(i,j,k,UMX) * g;
+}
+
+#endif

From c076ced70713961bd9f58e113218d465071519b6 Mon Sep 17 00:00:00 2001
From: Ryan Brady <ryan.m.brady@stonybrook.edu>
Date: Wed, 31 Jul 2024 12:57:49 -0400
Subject: [PATCH 02/14] update README

---
 Exec/science/gradient_detonation/README.md | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/Exec/science/gradient_detonation/README.md b/Exec/science/gradient_detonation/README.md
index f956ab9290..4e78d45ff4 100644
--- a/Exec/science/gradient_detonation/README.md
+++ b/Exec/science/gradient_detonation/README.md
@@ -1,6 +1,6 @@
-# `Detonation`
+# `Gradient Detonation`
 
-This is a similar setup to Detonation, except that we want to study
+This is a similar setup to `Detonation`, except that we want to study
 the impact of composition gradient on shock detection by applying
 a sinusoidal function to the composition.
 
@@ -12,9 +12,9 @@ temperature (T_l) is placed in the left side of the domain, and an
 ambient temperature (T_r) is in the right.  The parameter center_T
 specifies where, as a fraction of the domain length, to put the
 interface, while width_T determines how wide the transition region is,
-as a fraction of the domain length. Finally, the parameter cfrac, nfrac,
-and ofrac, are the initial C12, N14, O16 fraction, and the remaining
-material is He4.
+as a fraction of the domain length. 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.

From b71c961f8cb74efa9b37aecf11500d61da5bffda Mon Sep 17 00:00:00 2001
From: Ryan Brady <ryan.m.brady@stonybrook.edu>
Date: Wed, 31 Jul 2024 14:02:45 -0400
Subject: [PATCH 03/14] use subch_simple and shock var

---
 Exec/science/gradient_detonation/GNUmakefile | 9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/Exec/science/gradient_detonation/GNUmakefile b/Exec/science/gradient_detonation/GNUmakefile
index 9cfb52c4d2..c70577646b 100644
--- a/Exec/science/gradient_detonation/GNUmakefile
+++ b/Exec/science/gradient_detonation/GNUmakefile
@@ -11,18 +11,15 @@ USE_MPI    = TRUE
 
 USE_REACT  = TRUE
 
+USE_SHOCK_VAR = TRUE
+
 CASTRO_HOME ?= ../../..
 
 # This sets the EOS directory in $(MICROPHYSICS_HOME)/eos
 EOS_DIR     := helmholtz
 
 # This sets the network directory in $(MICROPHYSICS_HOME)/networks
-ifeq ($(USE_NSE_NET), TRUE)
-	NETWORK_DIR := subch_base
-	SCREEN_METHOD := chabrier1998
-else
-	NETWORK_DIR := aprox19
-endif
+NETWORK_DIR := subch_simple
 
 PROBLEM_DIR ?= ./
 

From 52d22fb21d96710de646216490c627f42efe67df Mon Sep 17 00:00:00 2001
From: Ryan Brady <ryan.m.brady@stonybrook.edu>
Date: Wed, 31 Jul 2024 14:06:58 -0400
Subject: [PATCH 04/14] minor fixes

---
 Exec/science/gradient_detonation/_prob_params               | 2 ++
 .../gradient_detonation/problem_initialize_state_data.H     | 6 +++---
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/Exec/science/gradient_detonation/_prob_params b/Exec/science/gradient_detonation/_prob_params
index 0f94589a81..f08946e8d1 100644
--- a/Exec/science/gradient_detonation/_prob_params
+++ b/Exec/science/gradient_detonation/_prob_params
@@ -10,6 +10,8 @@ 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
diff --git a/Exec/science/gradient_detonation/problem_initialize_state_data.H b/Exec/science/gradient_detonation/problem_initialize_state_data.H
index 76d34c3daa..61ed21fe0d 100644
--- a/Exec/science/gradient_detonation/problem_initialize_state_data.H
+++ b/Exec/science/gradient_detonation/problem_initialize_state_data.H
@@ -23,19 +23,19 @@ void problem_initialize_state_data (int i, int j, int k,
     state(i,j,k,URHO) = problem::dens;
 
     Real sigma = 1.0_rt / (1.0_rt + std::exp(-(c_T - xcen)/ width));
-    Real freq = 1.0_rt * 2.0 * 3.1415926535897932384
+    Real freq = 1.0_rt * 2.0 * 3.1415926535897932384;
 
     state(i,j,k,UTEMP) = problem::T_l + (problem::T_r - problem::T_l) * (1.0_rt - sigma);
 
     for (int n = 0; n < NumSpec; n++) {
-        state(i,j,k,UFS+n) = state(i,j,k,URHO) * (problem::xn_envelope[n] * std::sin(freq * i) + problem::xn_core[n] * (1.0 - std::sin(freq * i)));
+        state(i,j,k,UFS+n) = state(i,j,k,URHO) * (problem::xn_env[n] * std::sin(freq * i) + problem::xn_core[n] * (1.0 - std::sin(freq * i)));
     }
 
     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] = problem::xn_core[n] + xn_env[n];
+        eos_state.xn[n] = problem::xn_core[n] + problem::xn_env[n];
     }
 
     eos(eos_input_rt, eos_state);

From 925edf3c2df2fc23fb92042cccf180d035ed9432 Mon Sep 17 00:00:00 2001
From: Ryan Brady <ryan.m.brady@stonybrook.edu>
Date: Wed, 31 Jul 2024 14:17:46 -0400
Subject: [PATCH 05/14] add inputs file and use SDC

---
 .../inputs_grad_zingale_wd                    | 96 +++++++++++++++++++
 1 file changed, 96 insertions(+)
 create mode 100644 Exec/science/gradient_detonation/inputs_grad_zingale_wd

diff --git a/Exec/science/gradient_detonation/inputs_grad_zingale_wd b/Exec/science/gradient_detonation/inputs_grad_zingale_wd
new file mode 100644
index 0000000000..a4e5b1cd9b
--- /dev/null
+++ b/Exec/science/gradient_detonation/inputs_grad_zingale_wd
@@ -0,0 +1,96 @@
+# ------------------  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 = 3.e-1
+
+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
+

From 2aabeb2eebd3083cb39989d0bdae8e5346b49d59 Mon Sep 17 00:00:00 2001
From: Ryan Brady <ryan.m.brady@stonybrook.edu>
Date: Wed, 31 Jul 2024 14:18:12 -0400
Subject: [PATCH 06/14] add simplified SDC

---
 Exec/science/gradient_detonation/GNUmakefile | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/Exec/science/gradient_detonation/GNUmakefile b/Exec/science/gradient_detonation/GNUmakefile
index c70577646b..6e5853c5b8 100644
--- a/Exec/science/gradient_detonation/GNUmakefile
+++ b/Exec/science/gradient_detonation/GNUmakefile
@@ -13,6 +13,8 @@ USE_REACT  = TRUE
 
 USE_SHOCK_VAR = TRUE
 
+USE_SIMPLIFIED_SDC = TRUE
+
 CASTRO_HOME ?= ../../..
 
 # This sets the EOS directory in $(MICROPHYSICS_HOME)/eos

From 746ab48de0baaacc54367c9fe0cdf09f5232fabf Mon Sep 17 00:00:00 2001
From: Ryan Brady <ryan.m.brady@stonybrook.edu>
Date: Wed, 31 Jul 2024 15:22:08 -0400
Subject: [PATCH 07/14] add parameter to adjust sine freq.

---
 Exec/science/gradient_detonation/_prob_params          | 10 ++++++----
 .../problem_initialize_state_data.H                    |  2 +-
 2 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/Exec/science/gradient_detonation/_prob_params b/Exec/science/gradient_detonation/_prob_params
index f08946e8d1..136d3a28ee 100644
--- a/Exec/science/gradient_detonation/_prob_params
+++ b/Exec/science/gradient_detonation/_prob_params
@@ -4,17 +4,19 @@ T_r                  real          5.e7_rt        y
 
 dens                 real          1.e8_rt        y
 
-cfrac_env           real           0.00_rt         y
+frequency            real          0.0_rt         y
 
-cfrac_core           real          0.45_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
+nfrac_core           real         0.0_rt          y
 
 ofrac_env            real          0.0_rt         y
 
-ofrac_core           real          0.45_rt         y
+ofrac_core           real          0.45_rt        y
 
 w_T                  real          5.e-4_rt       y
 
diff --git a/Exec/science/gradient_detonation/problem_initialize_state_data.H b/Exec/science/gradient_detonation/problem_initialize_state_data.H
index 61ed21fe0d..c4795512c6 100644
--- a/Exec/science/gradient_detonation/problem_initialize_state_data.H
+++ b/Exec/science/gradient_detonation/problem_initialize_state_data.H
@@ -23,7 +23,7 @@ void problem_initialize_state_data (int i, int j, int k,
     state(i,j,k,URHO) = problem::dens;
 
     Real sigma = 1.0_rt / (1.0_rt + std::exp(-(c_T - xcen)/ width));
-    Real freq = 1.0_rt * 2.0 * 3.1415926535897932384;
+    Real freq = problem::frequency * 2.0 * 3.1415926535897932384;
 
     state(i,j,k,UTEMP) = problem::T_l + (problem::T_r - problem::T_l) * (1.0_rt - sigma);
 

From 3f04cd1e3e17c3c27b13b6516c0c033e7b660cb9 Mon Sep 17 00:00:00 2001
From: Ryan Brady <ryan.m.brady@stonybrook.edu>
Date: Thu, 1 Aug 2024 14:07:53 -0400
Subject: [PATCH 08/14] adjust ambient comp. and ensure eos_state.xn sums to 1

---
 Exec/science/gradient_detonation/problem_initialize.H  |  7 +------
 .../problem_initialize_state_data.H                    | 10 ++++++++++
 2 files changed, 11 insertions(+), 6 deletions(-)

diff --git a/Exec/science/gradient_detonation/problem_initialize.H b/Exec/science/gradient_detonation/problem_initialize.H
index 1f4165f27e..30472d3ed6 100644
--- a/Exec/science/gradient_detonation/problem_initialize.H
+++ b/Exec/science/gradient_detonation/problem_initialize.H
@@ -92,7 +92,7 @@ void problem_initialize ()
 
     problem::ambient_dens = problem::dens;
     for (int n = 0; n < NumSpec; n++) {
-        problem::ambient_comp[n] = problem::xn_core[n] + problem::xn_env[n];
+        problem::ambient_comp[n] = problem::xn_env[n];
     }
 
     eos_t eos_state;
@@ -101,11 +101,6 @@ void problem_initialize ()
         eos_state.xn[n]  = problem::ambient_comp[n];
     }
 
-#ifdef AUX_THERMO
-    // set the aux quantities -- we need to do this if we are using the NSE network
-    set_aux_comp_from_X(eos_state);
-#endif
-
     eos_state.T = problem::T_l;
 
     eos(eos_input_rt, eos_state);
diff --git a/Exec/science/gradient_detonation/problem_initialize_state_data.H b/Exec/science/gradient_detonation/problem_initialize_state_data.H
index c4795512c6..6b8001c2f9 100644
--- a/Exec/science/gradient_detonation/problem_initialize_state_data.H
+++ b/Exec/science/gradient_detonation/problem_initialize_state_data.H
@@ -38,6 +38,16 @@ void problem_initialize_state_data (int i, int j, int k,
         eos_state.xn[n] = problem::xn_core[n] + problem::xn_env[n];
     }
 
+    // normalize eos_state.xn
+    Real sum_xn = 0.0_rt;
+    for (int n = 0; n < NumSpec; n++) {
+        sum_xn += eos_state.xn[n]
+    }
+
+    for (int n = 0; n < NumSpec; n++) {
+        eos_state.xn[n] /= sum_xn;
+    }
+
     eos(eos_input_rt, eos_state);
 
     state(i,j,k,UMX) = state(i,j,k,URHO) * (problem::vel - 2 * problem::vel * (1.0_rt - sigma));

From c6322431e85882d5d6401b87689b53963ac04c1a Mon Sep 17 00:00:00 2001
From: Ryan Brady <ryan.m.brady@stonybrook.edu>
Date: Thu, 1 Aug 2024 14:10:23 -0400
Subject: [PATCH 09/14] forgot ;

---
 .../science/gradient_detonation/problem_initialize_state_data.H | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Exec/science/gradient_detonation/problem_initialize_state_data.H b/Exec/science/gradient_detonation/problem_initialize_state_data.H
index 6b8001c2f9..d4f00bbec2 100644
--- a/Exec/science/gradient_detonation/problem_initialize_state_data.H
+++ b/Exec/science/gradient_detonation/problem_initialize_state_data.H
@@ -41,7 +41,7 @@ void problem_initialize_state_data (int i, int j, int k,
     // normalize eos_state.xn
     Real sum_xn = 0.0_rt;
     for (int n = 0; n < NumSpec; n++) {
-        sum_xn += eos_state.xn[n]
+        sum_xn += eos_state.xn[n];
     }
 
     for (int n = 0; n < NumSpec; n++) {

From 4795e9423229fa5eb5d7eed68b5cf1b274f819e2 Mon Sep 17 00:00:00 2001
From: Ryan Brady <ryan.m.brady@stonybrook.edu>
Date: Thu, 1 Aug 2024 15:56:56 -0400
Subject: [PATCH 10/14] adjust default freq.

---
 Exec/science/gradient_detonation/_prob_params | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Exec/science/gradient_detonation/_prob_params b/Exec/science/gradient_detonation/_prob_params
index 136d3a28ee..0c191ef31f 100644
--- a/Exec/science/gradient_detonation/_prob_params
+++ b/Exec/science/gradient_detonation/_prob_params
@@ -4,7 +4,7 @@ T_r                  real          5.e7_rt        y
 
 dens                 real          1.e8_rt        y
 
-frequency            real          0.0_rt         y
+frequency            real          1.0_rt         y
 
 cfrac_env            real           0.00_rt       y
 

From 033719ca0277ca7110891551ddb130a954026ead Mon Sep 17 00:00:00 2001
From: Ryan Brady <ryan.m.brady@stonybrook.edu>
Date: Mon, 19 Aug 2024 11:47:35 -0400
Subject: [PATCH 11/14] use linear interpolation to transition from envelope to
 core

---
 Exec/science/gradient_detonation/README.md       |  2 +-
 Exec/science/gradient_detonation/_prob_params    |  2 +-
 .../problem_initialize_state_data.H              | 16 ++++++++++++----
 3 files changed, 14 insertions(+), 6 deletions(-)

diff --git a/Exec/science/gradient_detonation/README.md b/Exec/science/gradient_detonation/README.md
index 4e78d45ff4..12a6dd3704 100644
--- a/Exec/science/gradient_detonation/README.md
+++ b/Exec/science/gradient_detonation/README.md
@@ -2,7 +2,7 @@
 
 This is a similar setup to `Detonation`, except that we want to study
 the impact of composition gradient on shock detection by applying
-a sinusoidal function to the composition.
+a linear interpolation function to the composition.
 
 A simple (carbon or helium) detonation.  The reaction network should
 be set through the GNUmakefile.
diff --git a/Exec/science/gradient_detonation/_prob_params b/Exec/science/gradient_detonation/_prob_params
index 0c191ef31f..48d3ef8f40 100644
--- a/Exec/science/gradient_detonation/_prob_params
+++ b/Exec/science/gradient_detonation/_prob_params
@@ -4,7 +4,7 @@ T_r                  real          5.e7_rt        y
 
 dens                 real          1.e8_rt        y
 
-frequency            real          1.0_rt         y
+trans_frac           real          0.5_rt         y
 
 cfrac_env            real           0.00_rt       y
 
diff --git a/Exec/science/gradient_detonation/problem_initialize_state_data.H b/Exec/science/gradient_detonation/problem_initialize_state_data.H
index d4f00bbec2..c0ba6aac3d 100644
--- a/Exec/science/gradient_detonation/problem_initialize_state_data.H
+++ b/Exec/science/gradient_detonation/problem_initialize_state_data.H
@@ -22,20 +22,28 @@ void problem_initialize_state_data (int i, int j, int k,
 
     state(i,j,k,URHO) = problem::dens;
 
-    Real sigma = 1.0_rt / (1.0_rt + std::exp(-(c_T - xcen)/ width));
-    Real freq = problem::frequency * 2.0 * 3.1415926535897932384;
+    Real sigma = 1.0_rt / (1.0_rt + std::exp(-(c_T - xcen) / width));
 
     state(i,j,k,UTEMP) = problem::T_l + (problem::T_r - problem::T_l) * (1.0_rt - sigma);
 
+    // linear interpolation factor
+    Real interp_factor = ((xcen - problo[0]) / (probhi[0] - problo[0])) / problem::trans_frac;
+
+    // clamp interp_factor between 0 and 1
+    interp_factor = std::min(std::max(interp_factor, 0.0_rt), 1.0_rt);
+
+    // create a smooth transition from xn_env (left side of domain)
+    // to xn_core (right side of domain)
     for (int n = 0; n < NumSpec; n++) {
-        state(i,j,k,UFS+n) = state(i,j,k,URHO) * (problem::xn_env[n] * std::sin(freq * i) + problem::xn_core[n] * (1.0 - std::sin(freq * i)));
+        state(i,j,k,UFS+n) = state(i,j,k,URHO) * ((1.0_rt - interp_factor) * problem::xn_env[n] +
+                                                  interp_factor * problem::xn_core[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] = problem::xn_core[n] + problem::xn_env[n];
+        eos_state.xn[n] = state(i,j,k,UFS+n) / state(i,j,k,URHO);
     }
 
     // normalize eos_state.xn

From 7a6a8a57dabd6bc310f67c512cab0e1bf342f50c Mon Sep 17 00:00:00 2001
From: Ryan Brady <ryan.m.brady@stonybrook.edu>
Date: Mon, 26 Aug 2024 21:05:32 -0400
Subject: [PATCH 12/14] modify initialization such that the transition region
 is in the center of the domain, with a width given by problem.delta

---
 Exec/science/gradient_detonation/_prob_params          |  2 +-
 .../science/gradient_detonation/inputs_grad_zingale_wd |  3 +--
 .../problem_initialize_state_data.H                    | 10 ++++++++--
 3 files changed, 10 insertions(+), 5 deletions(-)

diff --git a/Exec/science/gradient_detonation/_prob_params b/Exec/science/gradient_detonation/_prob_params
index 48d3ef8f40..db69f50d0d 100644
--- a/Exec/science/gradient_detonation/_prob_params
+++ b/Exec/science/gradient_detonation/_prob_params
@@ -4,7 +4,7 @@ T_r                  real          5.e7_rt        y
 
 dens                 real          1.e8_rt        y
 
-trans_frac           real          0.5_rt         y
+delta                real          50.e5_rt       y
 
 cfrac_env            real           0.00_rt       y
 
diff --git a/Exec/science/gradient_detonation/inputs_grad_zingale_wd b/Exec/science/gradient_detonation/inputs_grad_zingale_wd
index a4e5b1cd9b..a0fb0c782d 100644
--- a/Exec/science/gradient_detonation/inputs_grad_zingale_wd
+++ b/Exec/science/gradient_detonation/inputs_grad_zingale_wd
@@ -9,7 +9,6 @@ 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
@@ -78,7 +77,7 @@ problem.smallx = 1.e-10
 problem.idir = 1
 
 problem.w_T = 1.e-3
-problem.center_T = 3.e-1
+problem.center_T = 0.05
 
 problem.nfrac_env = 0.01
 
diff --git a/Exec/science/gradient_detonation/problem_initialize_state_data.H b/Exec/science/gradient_detonation/problem_initialize_state_data.H
index c0ba6aac3d..28cbc54612 100644
--- a/Exec/science/gradient_detonation/problem_initialize_state_data.H
+++ b/Exec/science/gradient_detonation/problem_initialize_state_data.H
@@ -26,8 +26,14 @@ void problem_initialize_state_data (int i, int j, int k,
 
     state(i,j,k,UTEMP) = problem::T_l + (problem::T_r - problem::T_l) * (1.0_rt - sigma);
 
+    // define the WD-accretion interface around the center of the domain
+    Real delta = problem::delta;
+    Real center = problo[0] + (probhi[0] - problo[0]) / 2; // center of domain
+    Real x_start = center - delta / 2.0_rt;
+    Real x_end = center + delta / 2.0_rt;
+
     // linear interpolation factor
-    Real interp_factor = ((xcen - problo[0]) / (probhi[0] - problo[0])) / problem::trans_frac;
+    Real interp_factor = (xcen - x_start) / (x_end - x_start);
 
     // clamp interp_factor between 0 and 1
     interp_factor = std::min(std::max(interp_factor, 0.0_rt), 1.0_rt);
@@ -36,7 +42,7 @@ void problem_initialize_state_data (int i, int j, int k,
     // to xn_core (right side of domain)
     for (int n = 0; n < NumSpec; n++) {
         state(i,j,k,UFS+n) = state(i,j,k,URHO) * ((1.0_rt - interp_factor) * problem::xn_env[n] +
-                                                  interp_factor * problem::xn_core[n]);
+                                              interp_factor * problem::xn_core[n]);
     }
 
     eos_t eos_state;

From c8029db58fcd0482ffa0061977dc5148fd3fa904 Mon Sep 17 00:00:00 2001
From: Ryan Brady <ryan.m.brady@stonybrook.edu>
Date: Tue, 27 Aug 2024 14:00:57 -0400
Subject: [PATCH 13/14] separate initial hotspot temperature from envelope
 temperature. make temperature change with species composition.

---
 Exec/science/gradient_detonation/README.md             |  8 +-------
 Exec/science/gradient_detonation/_prob_params          |  6 ++++--
 .../problem_initialize_state_data.H                    | 10 ++++++++--
 3 files changed, 13 insertions(+), 11 deletions(-)

diff --git a/Exec/science/gradient_detonation/README.md b/Exec/science/gradient_detonation/README.md
index 12a6dd3704..a7bb2c57b6 100644
--- a/Exec/science/gradient_detonation/README.md
+++ b/Exec/science/gradient_detonation/README.md
@@ -8,13 +8,7 @@ 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_l) is placed in the left side of the domain, and an
-ambient temperature (T_r) is in the right.  The parameter center_T
-specifies where, as a fraction of the domain length, to put the
-interface, while width_T determines how wide the transition region is,
-as a fraction of the domain length. 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.
+temperature (T_ign) 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.
diff --git a/Exec/science/gradient_detonation/_prob_params b/Exec/science/gradient_detonation/_prob_params
index db69f50d0d..4194964e16 100644
--- a/Exec/science/gradient_detonation/_prob_params
+++ b/Exec/science/gradient_detonation/_prob_params
@@ -1,6 +1,8 @@
-T_l                  real          1.e9_rt        y
+T_l                  real          1.75e8_rt        y
 
-T_r                  real          5.e7_rt        y
+T_r                  real          1.e7_rt        y
+
+T_ign                real          1.75e9_rt      y
 
 dens                 real          1.e8_rt        y
 
diff --git a/Exec/science/gradient_detonation/problem_initialize_state_data.H b/Exec/science/gradient_detonation/problem_initialize_state_data.H
index 28cbc54612..d719e1aa86 100644
--- a/Exec/science/gradient_detonation/problem_initialize_state_data.H
+++ b/Exec/science/gradient_detonation/problem_initialize_state_data.H
@@ -24,8 +24,6 @@ void problem_initialize_state_data (int i, int j, int k,
 
     Real sigma = 1.0_rt / (1.0_rt + std::exp(-(c_T - xcen) / width));
 
-    state(i,j,k,UTEMP) = problem::T_l + (problem::T_r - problem::T_l) * (1.0_rt - sigma);
-
     // define the WD-accretion interface around the center of the domain
     Real delta = problem::delta;
     Real center = problo[0] + (probhi[0] - problo[0]) / 2; // center of domain
@@ -38,8 +36,16 @@ void problem_initialize_state_data (int i, int j, int k,
     // clamp interp_factor between 0 and 1
     interp_factor = std::min(std::max(interp_factor, 0.0_rt), 1.0_rt);
 
+    Real temperature = problem::T_l + (problem::T_r - problem::T_l) * interp_factor;
+
+    if (i < 2) {
+        temperature += problem::T_ign;
+    }
+    
     // create a smooth transition from xn_env (left side of domain)
     // to xn_core (right side of domain)
+    state(i,j,k,UTEMP) = temperature;
+
     for (int n = 0; n < NumSpec; n++) {
         state(i,j,k,UFS+n) = state(i,j,k,URHO) * ((1.0_rt - interp_factor) * problem::xn_env[n] +
                                               interp_factor * problem::xn_core[n]);

From f8a4d87d65516fae1b940378c080dea69cbf4d8e Mon Sep 17 00:00:00 2001
From: Ryan Brady <ryan.m.brady@stonybrook.edu>
Date: Tue, 27 Aug 2024 14:16:38 -0400
Subject: [PATCH 14/14] add new parameter to set size of hotspot independent of
 resolution

---
 Exec/science/gradient_detonation/README.md                    | 2 +-
 Exec/science/gradient_detonation/_prob_params                 | 4 +++-
 .../gradient_detonation/problem_initialize_state_data.H       | 4 ++--
 3 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/Exec/science/gradient_detonation/README.md b/Exec/science/gradient_detonation/README.md
index a7bb2c57b6..837c5efbaa 100644
--- a/Exec/science/gradient_detonation/README.md
+++ b/Exec/science/gradient_detonation/README.md
@@ -8,7 +8,7 @@ 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) 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.
+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.
diff --git a/Exec/science/gradient_detonation/_prob_params b/Exec/science/gradient_detonation/_prob_params
index 4194964e16..19fb582b5c 100644
--- a/Exec/science/gradient_detonation/_prob_params
+++ b/Exec/science/gradient_detonation/_prob_params
@@ -1,9 +1,11 @@
-T_l                  real          1.75e8_rt        y
+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
diff --git a/Exec/science/gradient_detonation/problem_initialize_state_data.H b/Exec/science/gradient_detonation/problem_initialize_state_data.H
index d719e1aa86..646ef6308f 100644
--- a/Exec/science/gradient_detonation/problem_initialize_state_data.H
+++ b/Exec/science/gradient_detonation/problem_initialize_state_data.H
@@ -38,10 +38,10 @@ void problem_initialize_state_data (int i, int j, int k,
 
     Real temperature = problem::T_l + (problem::T_r - problem::T_l) * interp_factor;
 
-    if (i < 2) {
+    if (xcen < problem::hotspot_size) {
         temperature += problem::T_ign;
     }
-    
+
     // create a smooth transition from xn_env (left side of domain)
     // to xn_core (right side of domain)
     state(i,j,k,UTEMP) = temperature;