circular detonation (#2858)

This is useful for testing the multidimensional shock detection for burning

Exec/science/Detonation/nse_runs/detonation.py

@@ -27,11 +27,13 @@ def __init__(self, plotfile):
         x_coord = np.array(ad['x'][srt])
         temp = np.array(ad['Temp'][srt])
         enuc = np.array(ad['enuc'][srt])
+        shock = np.array(ad['Shock'][srt])
 
         self.time = time
         self.x = x_coord
         self.T = temp
         self.enuc = enuc
+        self.shock = shock
 
     def find_x_for_T(self, T_0=1.e9):
         """ given a profile x(T), find the x_0 that corresponds to T_0 """

Exec/science/circular_det/GNUmakefile

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

Exec/science/circular_det/Make.package

@@ -0,0 +1 @@
+

Exec/science/circular_det/README.md

@@ -0,0 +1,5 @@
+# Circular Detonation
+
+This setup initializes a two-dimensional circular (or elliptical)
+detonation to help us understand how shock burning works in a
+multidimensional simulation.

Exec/science/circular_det/_prob_params

@@ -0,0 +1,32 @@
+T_l                  real          1.e9_rt        y
+
+T_r                  real          5.e7_rt        y
+
+dens                 real          1.e8_rt        y
+
+cfrac                real          0.5_rt         y
+
+nfrac                real          0.0_rt         y
+
+ofrac                real          0.0_rt         y
+
+# dimensionless width of transition between hot and cold
+w_T                  real          5.e-4_rt       y
+
+# dimensionless semi-major axis
+a_T                  real          0.3_rt         y
+
+# eccentricity
+ecc_T                real          0.0            y
+
+smallx               real          1.e-12_rt      y
+
+ihe4                 integer       -1
+
+ic12                 integer       -1
+
+in14                 integer       -1
+
+io16                 integer       -1
+
+xn                   real          0.0_rt         n     nspec

Exec/science/circular_det/analysis/contour_compare.py

@@ -0,0 +1,45 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+import yt
+
+ref_pf = "det_x_3km_plt00133"
+compare_pf = "det_x_24km_plt00038"
+
+field = "Temp"
+
+fig, ax = plt.subplots()
+
+legend_elem = []
+legend_labels = []
+
+for pf, label, color in [(ref_pf, "3 km", "C0"),
+                         (compare_pf, "24 km", "C1")]:
+
+    ds = yt.load(pf)
+
+    xmin = ds.domain_left_edge[0]
+    xmax = ds.domain_right_edge[0]
+
+    ymin = ds.domain_left_edge[1]
+    ymax = ds.domain_right_edge[1]
+
+    ref = int(np.prod(ds.ref_factors[0:ds.index.max_level]))
+
+    data = ds.covering_grid(ds.index.max_level,
+                            left_edge=ds.domain_left_edge,
+                            dims=ds.domain_dimensions*ref, fields=field)
+
+    print(data.shape)
+
+    ct = ax.contour(data[field].d[:,:,0].T, levels=4, colors=color,
+                    extent=[xmin, xmax, ymin, ymax])
+
+    legend_elem.append(ct.legend_elements()[0][0])
+    legend_labels.append(label)
+
+ax.set_aspect("equal")
+
+ax.legend(legend_elem, legend_labels)
+fig.savefig("contour.png")
+

Exec/science/circular_det/inputs.24km

@@ -0,0 +1,123 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 200000
+stop_time =  0.024
+
+# 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     = 1.536e8  1.536e8
+amr.n_cell           = 64       64
+
+
+
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+# 0 = Interior           3 = Symmetry
+# 1 = Inflow             4 = SlipWall
+# 2 = Outflow            5 = NoSlipWall
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+castro.lo_bc       =  2   2   4
+castro.hi_bc       =  2   2   4
+
+# WHICH PHYSICS
+castro.do_hydro = 1
+castro.do_react = 1
+
+castro.ppm_type = 1
+castro.ppm_temp_fix = 0
+
+castro.time_integration_method = 3
+
+castro.disable_shock_burning = 1
+castro.shock_detection_threshold = 1
+
+# castro.transverse_reset_density = 1
+
+castro.small_dens   = 1.e-5
+castro.small_temp     = 1.e7
+
+castro.use_flattening = 1
+
+castro.riemann_solver = 1
+
+# TIME STEP CONTROL
+castro.cfl            = 0.4     # cfl number for hyperbolic system
+castro.init_shrink    = 0.1     # 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
+#amr.grid_log        = grdlog  # name of grid logging file
+
+# 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     = 8 8 8 2 2 2 # number of buffer cells in error est
+
+# CHECKPOINT FILES
+amr.check_file      = det_x_24km_chk  # root name of checkpoint file
+amr.check_int       = 1000         # number of timesteps between checkpoints
+
+# PLOTFILES
+amr.plot_file       = det_x_24km_plt  # root name of plotfile
+amr.plot_per = 2.e-3
+amr.derive_plot_vars = ALL
+
+# problem initialization
+
+problem.T_l = 1.1e9
+problem.T_r = 1.75e8
+
+problem.dens = 3.e6
+problem.cfrac = 0.0
+problem.nfrac = 0.01
+
+problem.smallx = 1.e-10
+
+problem.idir = 1
+
+problem.w_T = 1.e-2
+problem.a_T = 0.15
+problem.ecc_T = 0.7
+
+# refinement
+
+amr.refinement_indicators = temperr tempgrad
+
+amr.refine.temperr.max_level = 0
+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
+
+# Microphysics
+
+network.small_x = 1.e-10
+integrator.SMALL_X_SAFE = 1.e-10
+
+integrator.rtol_spec = 1.e-5
+integrator.atol_spec = 1.e-5
+integrator.rtol_enuc = 1.e-5
+integrator.atol_enuc = 1.e-5
+integrator.jacobian = 1
+
+integrator.X_reject_buffer = 4.0
+
+integrator.use_burn_retry = 1
+integrator.retry_swap_jacobian = 1
+
+integrator.retry_rtol_spec = 1.e-5
+integrator.retry_atol_spec = 1.e-5
+integrator.retry_rtol_enuc = 1.e-5
+integrator.retry_atol_enuc = 1.e-5
+
+integrator.ode_max_steps = 10000

Exec/science/circular_det/inputs.3km

@@ -0,0 +1,123 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 200000
+stop_time =  0.024
+
+# 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     = 1.536e8  1.536e8
+amr.n_cell           = 256      256
+
+
+
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+# 0 = Interior           3 = Symmetry
+# 1 = Inflow             4 = SlipWall
+# 2 = Outflow            5 = NoSlipWall
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+castro.lo_bc       =  2   2   4
+castro.hi_bc       =  2   2   4
+
+# WHICH PHYSICS
+castro.do_hydro = 1
+castro.do_react = 1
+
+castro.ppm_type = 1
+castro.ppm_temp_fix = 0
+
+castro.time_integration_method = 3
+
+castro.disable_shock_burning = 1
+castro.shock_detection_threshold = 1
+
+# castro.transverse_reset_density = 1
+
+castro.small_dens   = 1.e-5
+castro.small_temp     = 1.e7
+
+castro.use_flattening = 1
+
+castro.riemann_solver = 1
+
+# TIME STEP CONTROL
+castro.cfl            = 0.4     # cfl number for hyperbolic system
+castro.init_shrink    = 0.1     # 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
+#amr.grid_log        = grdlog  # name of grid logging file
+
+# REFINEMENT / REGRIDDING
+amr.max_level       = 1       # 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     = 8 8 8 2 2 2 # number of buffer cells in error est
+
+# CHECKPOINT FILES
+amr.check_file      = det_x_3km_chk  # root name of checkpoint file
+amr.check_int       = 1000         # number of timesteps between checkpoints
+
+# PLOTFILES
+amr.plot_file       = det_x_3km_plt  # root name of plotfile
+amr.plot_per = 2.e-3
+amr.derive_plot_vars = ALL
+
+# problem initialization
+
+problem.T_l = 1.1e9
+problem.T_r = 1.75e8
+
+problem.dens = 3.e6
+problem.cfrac = 0.0
+problem.nfrac = 0.01
+
+problem.smallx = 1.e-10
+
+problem.idir = 1
+
+problem.w_T = 1.e-2
+problem.a_T = 0.15
+problem.ecc_T = 0.7
+
+# refinement
+
+amr.refinement_indicators = temperr tempgrad
+
+amr.refine.temperr.max_level = 0
+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
+
+# Microphysics
+
+network.small_x = 1.e-10
+integrator.SMALL_X_SAFE = 1.e-10
+
+integrator.rtol_spec = 1.e-5
+integrator.atol_spec = 1.e-5
+integrator.rtol_enuc = 1.e-5
+integrator.atol_enuc = 1.e-5
+integrator.jacobian = 1
+
+integrator.X_reject_buffer = 4.0
+
+integrator.use_burn_retry = 1
+integrator.retry_swap_jacobian = 1
+
+integrator.retry_rtol_spec = 1.e-5
+integrator.retry_atol_spec = 1.e-5
+integrator.retry_rtol_enuc = 1.e-5
+integrator.retry_atol_enuc = 1.e-5
+
+integrator.ode_max_steps = 10000