diff --git a/Exec/science/gradient_detonation/GNUmakefile b/Exec/science/gradient_detonation/GNUmakefile new file mode 100644 index 0000000000..6e5853c5b8 --- /dev/null +++ b/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 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..837c5efbaa --- /dev/null +++ b/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. diff --git a/Exec/science/gradient_detonation/_prob_params b/Exec/science/gradient_detonation/_prob_params new file mode 100644 index 0000000000..19fb582b5c --- /dev/null +++ b/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 \ No newline at end of file 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..a0fb0c782d --- /dev/null +++ b/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 + 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 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 const& state, + Real time, + const Array1D& 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(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..30472d3ed6 --- /dev/null +++ b/Exec/science/gradient_detonation/problem_initialize.H @@ -0,0 +1,118 @@ +#ifndef problem_initialize_H +#define problem_initialize_H + +#include +#include +#include + +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_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]; + } + + 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..646ef6308f --- /dev/null +++ b/Exec/science/gradient_detonation/problem_initialize_state_data.H @@ -0,0 +1,81 @@ +#ifndef problem_initialize_state_data_H +#define problem_initialize_state_data_H + +#include +#include +#include + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void problem_initialize_state_data (int i, int j, int k, + Array4 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(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)); + + // 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 - 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); + + Real temperature = problem::T_l + (problem::T_r - problem::T_l) * interp_factor; + + 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; + + 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]); + } + + 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) / state(i,j,k,URHO); + } + + // 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)); + 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 + +using namespace amrex; + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void problem_source (int i, int j, int k, + GeometryData const& geomdata, + Array4 const& state, + Array4 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