diff --git a/Exec/science/Detonation/nse_runs/detonation.py b/Exec/science/Detonation/nse_runs/detonation.py index d8b58cee97..2fc2c57375 100644 --- a/Exec/science/Detonation/nse_runs/detonation.py +++ b/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 """ diff --git a/Exec/science/circular_det/GNUmakefile b/Exec/science/circular_det/GNUmakefile new file mode 100644 index 0000000000..46f051fb7e --- /dev/null +++ b/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 diff --git a/Exec/science/circular_det/Make.package b/Exec/science/circular_det/Make.package new file mode 100644 index 0000000000..8b13789179 --- /dev/null +++ b/Exec/science/circular_det/Make.package @@ -0,0 +1 @@ + diff --git a/Exec/science/circular_det/README.md b/Exec/science/circular_det/README.md new file mode 100644 index 0000000000..08a27b90f2 --- /dev/null +++ b/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. diff --git a/Exec/science/circular_det/_prob_params b/Exec/science/circular_det/_prob_params new file mode 100644 index 0000000000..44e0db3a3a --- /dev/null +++ b/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 diff --git a/Exec/science/circular_det/analysis/contour_compare.py b/Exec/science/circular_det/analysis/contour_compare.py new file mode 100644 index 0000000000..53554710c6 --- /dev/null +++ b/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") + diff --git a/Exec/science/circular_det/inputs.24km b/Exec/science/circular_det/inputs.24km new file mode 100644 index 0000000000..190c37d2ec --- /dev/null +++ b/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 diff --git a/Exec/science/circular_det/inputs.3km b/Exec/science/circular_det/inputs.3km new file mode 100644 index 0000000000..6460fa3420 --- /dev/null +++ b/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 diff --git a/Exec/science/circular_det/problem_initialize.H b/Exec/science/circular_det/problem_initialize.H new file mode 100644 index 0000000000..eb2327ded1 --- /dev/null +++ b/Exec/science/circular_det/problem_initialize.H @@ -0,0 +1,80 @@ +#ifndef problem_initialize_H +#define problem_initialize_H + +#include +#include +#include + +AMREX_INLINE +void problem_initialize () +{ + + const Geometry& dgeom = DefaultGeometry(); + + const Real* problo = dgeom.ProbLo(); + const Real* probhi = dgeom.ProbHi(); + + for (int n = 0; n < AMREX_SPACEDIM; ++n) { + problem::center[n] = 0.5_rt * (problo[n] + probhi[n]); + } + + // 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 + + if (problem::cfrac > 1.0_rt || problem::cfrac < 0.0_rt) { + amrex::Error("ERROR: cfrac must fall between 0 and 1"); + } + + // make sure that the nitrogen fraction falls between 0 and 1 + + if (problem::nfrac > 1.0_rt || problem::nfrac < 0.0_rt) { + amrex::Error("ERROR: nfrac must fall between 0 and 1"); + } + + // make sure that the oxygen fraction falls between 0 and 1 + + if (problem::ofrac > 1.0_rt || problem::ofrac < 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 + + if (problem::cfrac + problem::nfrac + problem::ofrac > 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[n] = problem::smallx; + } + + problem::xn[problem::ic12] = amrex::max(problem::cfrac, problem::smallx); + problem::xn[problem::io16] = amrex::max(problem::ofrac, problem::smallx); + + if (problem::in14 >= 0) { + problem::xn[problem::in14] = amrex::max(problem::nfrac, problem::smallx); + problem::xn[problem::ihe4] = 1.0_rt - problem::xn[problem::ic12] + - problem::xn[problem::in14] + - problem::xn[problem::io16] + - (NumSpec - 4) * problem::smallx; + } + else { + problem::xn[problem::ihe4] = 1.0_rt - problem::xn[problem::ic12] + - problem::xn[problem::io16] + - (NumSpec - 3) * problem::smallx; + } + +} + +#endif diff --git a/Exec/science/circular_det/problem_initialize_state_data.H b/Exec/science/circular_det/problem_initialize_state_data.H new file mode 100644 index 0000000000..579e3346c8 --- /dev/null +++ b/Exec/science/circular_det/problem_initialize_state_data.H @@ -0,0 +1,67 @@ +#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(); + + // semi-major axis + Real a = problem::a_T * (probhi[0] - problo[0]); + + // semi-minor axis + Real b = a * std::sqrt(1.0 - problem::ecc_T * problem::ecc_T); + + // compute distance from the center + Real xx = problo[0] + dx[0] * (static_cast(i) + 0.5_rt) - problem::center[0]; + Real yy = problo[1] + dx[1] * (static_cast(j) + 0.5_rt) - problem::center[1]; + + Real r = std::sqrt((xx / a) * (xx / a) + + (yy / b) * (yy / b)); + + state(i,j,k,URHO) = problem::dens; + + Real sigma = 1.0_rt / (1.0_rt + std::exp(-(1.0_rt - r)/ problem::w_T)); + + 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[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[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); + + for (int n = 0; n < NumAux; n++) { + state(i,j,k,UFX+n) = state(i,j,k,URHO) * eos_state.aux[n]; + } +#endif + + eos(eos_input_rt, eos_state); + + state(i,j,k,UMX) = 0.0_rt; + 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); + +} + +#endif diff --git a/Exec/science/subchandra/GNUmakefile b/Exec/science/subchandra/GNUmakefile index c7fe376fd3..f667161d2f 100644 --- a/Exec/science/subchandra/GNUmakefile +++ b/Exec/science/subchandra/GNUmakefile @@ -13,6 +13,7 @@ USE_MPI = TRUE USE_GRAV = TRUE USE_REACT = TRUE +USE_SHOCK_VAR = TRUE USE_MODEL_PARSER = TRUE