From 3ff0599108f46815b222324323041f0710a8831e Mon Sep 17 00:00:00 2001 From: Zhi Date: Mon, 7 Oct 2024 14:42:27 -0400 Subject: [PATCH] initial xrb spherical problem setup --- Exec/science/xrb_spherical/GNUmakefile | 40 ++++ Exec/science/xrb_spherical/Make.package | 2 + Exec/science/xrb_spherical/_prob_params | 68 ++++++ Exec/science/xrb_spherical/initial_model.H | 1 + Exec/science/xrb_spherical/inputs.He.1000Hz | 147 ++++++++++++ .../xrb_spherical/problem_initialize.H | 222 ++++++++++++++++++ .../problem_initialize_state_data.H | 78 ++++++ Exec/science/xrb_spherical/problem_tagging.H | 56 +++++ 8 files changed, 614 insertions(+) create mode 100644 Exec/science/xrb_spherical/GNUmakefile create mode 100644 Exec/science/xrb_spherical/Make.package create mode 100644 Exec/science/xrb_spherical/_prob_params create mode 120000 Exec/science/xrb_spherical/initial_model.H create mode 100644 Exec/science/xrb_spherical/inputs.He.1000Hz create mode 100644 Exec/science/xrb_spherical/problem_initialize.H create mode 100644 Exec/science/xrb_spherical/problem_initialize_state_data.H create mode 100644 Exec/science/xrb_spherical/problem_tagging.H diff --git a/Exec/science/xrb_spherical/GNUmakefile b/Exec/science/xrb_spherical/GNUmakefile new file mode 100644 index 0000000000..b6e0d4bf4a --- /dev/null +++ b/Exec/science/xrb_spherical/GNUmakefile @@ -0,0 +1,40 @@ +PRECISION = DOUBLE +PROFILE = FALSE + +DEBUG = FALSE + +DIM = 2 + +COMP = gnu + +USE_MPI = TRUE + +USE_GRAV = TRUE +USE_REACT = TRUE + +USE_ROTATION = TRUE +USE_DIFFUSION = TRUE + +# define the location of the CASTRO top directory +CASTRO_HOME ?= ../../.. + +USE_JACOBIAN_CACHING = TRUE +USE_MODEL_PARSER = TRUE +NUM_MODELS := 2 + +# This sets the EOS directory in $(MICROPHYSICS_HOME)/eos +EOS_DIR := helmholtz + +# This sets the network directory in $(MICROPHYSICS_HOME)/networks +NETWORK_DIR := aprox13 + +INTEGRATOR_DIR := VODE + +CONDUCTIVITY_DIR := stellar + +PROBLEM_DIR ?= ./ + +Bpack := $(PROBLEM_DIR)/Make.package +Blocs := $(PROBLEM_DIR) + +include $(CASTRO_HOME)/Exec/Make.Castro diff --git a/Exec/science/xrb_spherical/Make.package b/Exec/science/xrb_spherical/Make.package new file mode 100644 index 0000000000..e5cc052427 --- /dev/null +++ b/Exec/science/xrb_spherical/Make.package @@ -0,0 +1,2 @@ +CEXE_headers += initial_model.H + diff --git a/Exec/science/xrb_spherical/_prob_params b/Exec/science/xrb_spherical/_prob_params new file mode 100644 index 0000000000..89c1775615 --- /dev/null +++ b/Exec/science/xrb_spherical/_prob_params @@ -0,0 +1,68 @@ + +dtemp real 3.81e8_rt y + +theta_half_max real 1.745e-2_rt y + +theta_half_width real 4.9e-3_rt y + +# cutoff mass fraction of the first species for refinement +X_min real 1.e-4_rt y + +# do we dynamically refine based on density? or based on height? +tag_by_density integer 1 y + +# used for tagging if tag_by_density = 1 +cutoff_density real 500.e0_rt y + +# used if we are refining based on height rather than density +refine_height real 3600 y + +T_hi real 5.e8_rt y + +T_star real 1.e8_rt y + +T_lo real 5.e7_rt y + +dens_base real 2.e6_rt y + +H_star real 500.e0_rt y + +atm_delta real 25.e0_rt y + +fuel1_name string "helium-4" y + +fuel2_name string "" y + +fuel3_name string "" y + +fuel4_name string "" y + +ash1_name string "iron-56" y + +ash2_name string "" y + +ash3_name string "" y + +fuel1_frac real 1.0_rt y + +fuel2_frac real 0.0_rt y + +fuel3_frac real 0.0_rt y + +fuel4_frac real 0.0_rt y + +ash1_frac real 1.0_rt y + +ash2_frac real 0.0_rt y + +ash3_frac real 0.0_rt y + +low_density_cutoff real 1.e-4_rt y + +smallx real 1.e-10_rt y + +r_refine_distance real 1.e30_rt y + +max_hse_tagging_level integer 2 y + +max_base_tagging_level integer 1 y diff --git a/Exec/science/xrb_spherical/initial_model.H b/Exec/science/xrb_spherical/initial_model.H new file mode 120000 index 0000000000..9c923c3113 --- /dev/null +++ b/Exec/science/xrb_spherical/initial_model.H @@ -0,0 +1 @@ +../flame_wave/initial_model.H \ No newline at end of file diff --git a/Exec/science/xrb_spherical/inputs.He.1000Hz b/Exec/science/xrb_spherical/inputs.He.1000Hz new file mode 100644 index 0000000000..4f465e61d0 --- /dev/null +++ b/Exec/science/xrb_spherical/inputs.He.1000Hz @@ -0,0 +1,147 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 9900000 +stop_time = 3.0 + +# PROBLEM SIZE & GEOMETRY +geometry.is_periodic = 0 0 +geometry.coord_sys = 2 # 0 => cart, 1 => RZ 2=>spherical +geometry.prob_lo = 0 0 +geometry.prob_hi = 3.072e4 3.14159 +amr.n_cell = 192 1152 + +# >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< +# 0 = Interior 3 = Symmetry +# 1 = Inflow 4 = SlipWall +# 2 = Outflow 5 = NoSlipWall +# >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< +castro.lo_bc = 1 3 # Inflow in lower R and Symmetry across Theta +castro.hi_bc = 2 2 # Outflow in upper boundaries + +# Hydrostatic condition along R-direction +castro.yl_ext_bc_type = 0 +castro.hse_interp_temp = 0 +castro.hse_fixed_temp = 1.e8 # this should match problem.T_star +castro.hse_reflect_vels = 0 + +# Fill ambient states with outflow velocity in R-direction +castro.fill_ambient_bc = 1 +castro.ambient_fill_dir = 0 +castro.ambient_outflow_vel = 1 + +# WHICH PHYSICS +castro.do_hydro = 1 +castro.do_react = 1 +castro.do_rotation = 1 +castro.do_grav = 1 +castro.do_sponge = 1 + +castro.small_temp = 1.e6 +castro.small_dens = 1.e-5 + +castro.ppm_type = 1 +castro.grav_source_type = 2 + +gravity.gravity_type = ConstantGrav + +# 1.4 Solar Mass NS with radius ~11 km +gravity.const_grav = -1.5e14 + +# 1000Hz Spinning Frequency +# Might want to use a more realistic spinning frequency like 500Hz +castro.rotational_period = 0.001 +# Centrifugal is not important since NS would simply deform to accomondate for it +castro.rotation_include_centrifugal = 0 + +castro.diffuse_temp = 1 +castro.diffuse_cutoff_density_hi = 5.e4 +castro.diffuse_cutoff_density = 2.e4 + +castro.diffuse_cond_scale_fac = 1.0 + +castro.react_rho_min = 1.e2 +castro.react_rho_max = 1.5e7 + +castro.react_T_min = 6.e7 + +castro.sponge_upper_density = 1.e2 +castro.sponge_lower_density = 1.e0 +castro.sponge_timescale = 1.e-7 + +castro.abundance_failure_tolerance = 0.1 +castro.abundance_failure_rho_cutoff = 1.0 + +# TIME STEP CONTROL +castro.cfl = 0.8 # cfl number for hyperbolic system +castro.init_shrink = 0.1 # scale back initial timestep +castro.change_max = 1.1 # max time step growth + +castro.use_retry = 1 +castro.max_subcycles = 16 + +castro.retry_small_density_cutoff = 1.0 + +# DIAGNOSTICS & VERBOSITY +castro.sum_interval = 0 # timesteps between computing mass +castro.v = 1 # verbosity in Castro.cpp +amr.v = 1 # verbosity in Amr.cpp + +# REFINEMENT / REGRIDDING +amr.max_level = 2 # maximum level number allowed +amr.ref_ratio = 4 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 = 128 +amr.n_error_buf = 2 2 2 2 # number of buffer cells in error est + +# CHECKPOINT FILES +amr.check_file = flame_wave_1000Hz_chk # root name of checkpoint file +amr.check_int = 250 # number of timesteps between checkpoints + +# PLOTFILES +amr.plot_file = flame_wave_1000Hz_plt # root name of plotfile +amr.plot_per = 5.e-3 # number of seconds between plotfiles +amr.derive_plot_vars = ALL + +amr.small_plot_file = flame_wave_1000Hz_smallplt # root name of plotfile +amr.small_plot_per = 1.e-4 # number of seconds between plotfiles +amr.small_plot_vars = density Temp +amr.derive_small_plot_vars = abar x_velocity y_velocity z_velocity enuc + +# problem initialization + +problem.dtemp = 1.2e9 +problem.theta_half_max = 1.745e-2 +problem.theta_half_width = 4.9e-3 + +problem.dens_base = 3.43e6 + +problem.T_star = 1.e8 +problem.T_hi = 2.e8 +problem.T_lo = 8.e6 + +problem.H_star = 2000.e0 +problem.atm_delta = 50.0 + +problem.fuel1_name = "helium-4" +problem.fuel1_frac = 1.0e0 + +problem.ash1_name = "nickel-56" +problem.ash1_frac = 1.0e0 + +problem.low_density_cutoff = 1.e-4 + +problem.cutoff_density = 2.5e4 +problem.max_hse_tagging_level = 3 +problem.max_base_tagging_level = 2 + +problem.X_min = 1.e-2 + +problem.x_refine_distance = 9.216e4 + + +# Microphysics + +integrator.rtol_spec = 1.e-6 +integrator.atol_spec = 1.e-6 + +network.use_tables = 1 diff --git a/Exec/science/xrb_spherical/problem_initialize.H b/Exec/science/xrb_spherical/problem_initialize.H new file mode 100644 index 0000000000..cfb11741ff --- /dev/null +++ b/Exec/science/xrb_spherical/problem_initialize.H @@ -0,0 +1,222 @@ +#ifndef problem_initialize_H +#define problem_initialize_H + +#include +#include +#include +#include +#include +#include +#include + +AMREX_INLINE +void problem_initialize () +{ + + const Geometry& dgeom = DefaultGeometry(); + + const Real* problo = dgeom.ProbLo(); + const Real* probhi = dgeom.ProbHi(); + + // check to make sure that small_dens is less than low_density_cutoff + // if not, funny things can happen above the atmosphere + + if (small_dens >= 0.99_rt * problem::low_density_cutoff) { + amrex::Error("ERROR: small_dens should be set lower than low_density_cutoff"); + } + + // make sure hse_fixed_temp is the same as T_star, if it's specified + + if (hse_fixed_temp > 0.0_rt && hse_fixed_temp != problem::T_star) { + amrex::Error("ERROR: hse_fixed_temp should be the same as T_star"); + } + + // get the species indices + + bool species_defined = true; + + int ifuel1 = network_spec_index(problem::fuel1_name); + if (ifuel1 < 0) { + species_defined = false; + } + + int ifuel2; + if (!problem::fuel2_name.empty()) { + ifuel2 = network_spec_index(problem::fuel2_name); + if (ifuel2 < 0) { + species_defined = false; + } + } + + int ifuel3; + if (!problem::fuel3_name.empty()) { + ifuel3 = network_spec_index(problem::fuel3_name); + if (ifuel3 < 0) { + species_defined = false; + } + } + + int ifuel4; + if (!problem::fuel4_name.empty()) { + ifuel4 = network_spec_index(problem::fuel4_name); + if (ifuel4 < 0) { + species_defined = false; + } + } + + int iash1 = network_spec_index(problem::ash1_name); + if (iash1 < 0) { + species_defined = false; + } + + int iash2; + if (!problem::ash2_name.empty()) { + iash2 = network_spec_index(problem::ash2_name); + if (iash2 < 0) { + species_defined = false; + } + } + + int iash3; + if (!problem::ash3_name.empty()) { + iash3 = network_spec_index(problem::ash3_name); + if (iash3 < 0) { + species_defined = false; + } + } + + if (! species_defined) { + std::cout << ifuel1 << " " << ifuel2 << " " << ifuel3 << " " << ifuel4 << std::endl; + std::cout << iash1 << " " << iash2 << " "<< iash3 << std::endl; + amrex::Error("ERROR: species not defined"); + } + + model_t model_params; + + // set the composition of the underlying star + + + for (Real &X : model_params.xn_star) { + X = problem::smallx; + } + model_params.xn_star[iash1] = problem::ash1_frac; + if (!problem::ash2_name.empty()) { + model_params.xn_star[iash2] = problem::ash2_frac; + } + if (!problem::ash3_name.empty()) { + model_params.xn_star[iash3] = problem::ash3_frac; + } + + // and the composition of the accreted layer + + for (Real &X : model_params.xn_base) { + X = problem::smallx; + } + model_params.xn_base[ifuel1] = problem::fuel1_frac; + if (!problem::fuel2_name.empty()) { + model_params.xn_base[ifuel2] = problem::fuel2_frac; + } + if (!problem::fuel3_name.empty()) { + model_params.xn_base[ifuel3] = problem::fuel3_frac; + } + if (!problem::fuel4_name.empty()) { + model_params.xn_base[ifuel4] = problem::fuel4_frac; + } + + // check if they sum to 1 + + Real sumX = 0.0_rt; + for (Real X : model_params.xn_star) { + sumX += X; + } + if (std::abs(sumX - 1.0_rt) > NumSpec * problem::smallx) { + amrex::Error("ERROR: ash mass fractions don't sum to 1"); + } + + sumX = 0.0_rt; + for (Real X : model_params.xn_base) { + sumX += X; + } + if (std::abs(sumX - 1.0_rt) > NumSpec * problem::smallx) { + amrex::Error("ERROR: fuel mass fractions don't sum to 1"); + } + + // we are going to generate an initial model from probl(1), i.e. r_min, to + // probhi(1), i.e., r_max, with nx_model zones. But to allow for a interpolated + // lower boundary, we'll add 4 ghostcells to this, so we need to + // compute dx + + // we use the fine grid dx for the model resolution + auto fine_geom = global::the_amr_ptr->Geom(global::the_amr_ptr->maxLevel()); + + auto dx = fine_geom.CellSizeArray(); + auto dx_model = dx[0]; + + int nx_model = static_cast((probhi[0] - problo[0]) / + dx_model); + + int ng = 4; + + // now generate the initial models + + model_params.dens_base = problem::dens_base; + model_params.T_star = problem::T_star; + model_params.T_hi = problem::T_hi; + model_params.T_lo = problem::T_lo; + + model_params.H_star = problem::H_star; + model_params.atm_delta = problem::atm_delta; + + model_params.low_density_cutoff = problem::low_density_cutoff; + + generate_initial_model(nx_model + ng, + problo[0] - ng * dx_model, + probhi[0], + model_params, 0); + + // now create a perturbed model -- we want the same base conditions + // a hotter temperature + + model_params.T_hi = model_params.T_hi + problem::dtemp; + + generate_initial_model(nx_model + ng, + problo[0] - ng * dx_model, + probhi[0], + model_params, 1); + + // set center + + for (int d = 0; d < AMREX_SPACEDIM; d++) { + // problem::center[d] = 0.5_rt * (problo[d] + probhi[d]); + problem::center[d] = 0.0_rt; + } + + // 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; +} + +#endif diff --git a/Exec/science/xrb_spherical/problem_initialize_state_data.H b/Exec/science/xrb_spherical/problem_initialize_state_data.H new file mode 100644 index 0000000000..088b1b42a9 --- /dev/null +++ b/Exec/science/xrb_spherical/problem_initialize_state_data.H @@ -0,0 +1,78 @@ +#ifndef problem_initialize_state_data_H +#define problem_initialize_state_data_H + +#include +#include +#include +#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(); + + Real r = problo[0] + dx[0] * (static_cast(i) + 0.5_rt) - problem::center[0]; + Real theta = problo[1] + dx[1] * (static_cast(j) + 0.5_rt) - problem::center[1]; + + // blending factor + + Real f; + + if (r < problem::theta_half_max) { + f = 1.0_rt; + + } else if (r > problem::theta_half_max + problem::theta_half_width) { + f = 0.0_rt; + + } else { + f = -(r - problem::theta_half_max) / problem::theta_half_width + 1.0_rt; + } + + state(i,j,k,URHO) = f * interpolate(r, model::idens, 1) + + (1.0_rt - f) * interpolate(r, model::idens, 0); + + state(i,j,k,UTEMP) = f * interpolate(r, model::itemp, 1) + + (1.0_rt - f) * interpolate(r, model::itemp, 0); + + Real temppres = f * interpolate(r, model::ipres, 1) + + (1.0_rt - f) * interpolate(r, model::ipres, 0); + + for (int n = 0; n < NumSpec; n++) { + state(i,j,k,UFS+n) = f * interpolate(r, model::ispec+n, 1) + + (1.0_rt - f) * interpolate(r, model::ispec+n, 0); + } + + eos_t eos_state; + eos_state.rho = state(i,j,k,URHO); + eos_state.T = state(i,j,k,UTEMP); + eos_state.p = temppres; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = state(i,j,k,UFS+n); + } + + eos(eos_input_rp, eos_state); + + state(i,j,k,UTEMP) = eos_state.T; + state(i,j,k,UEINT) = eos_state.rho * eos_state.e; + state(i,j,k,UEDEN) = state(i,j,k,UEINT); + + // Initial velocities = 0 + + state(i,j,k,UMX) = 0.e0_rt; + state(i,j,k,UMY) = 0.e0_rt; + state(i,j,k,UMZ) = 0.e0_rt; + + // convert to partial densities + + for (int n = 0; n < NumSpec; n++) { + state(i,j,k,UFS+n) = state(i,j,k,URHO) * state(i,j,k,UFS+n); + } +} + +#endif diff --git a/Exec/science/xrb_spherical/problem_tagging.H b/Exec/science/xrb_spherical/problem_tagging.H new file mode 100644 index 0000000000..877bf1c8cc --- /dev/null +++ b/Exec/science/xrb_spherical/problem_tagging.H @@ -0,0 +1,56 @@ +#ifndef problem_tagging_H +#define problem_tagging_H +#include +#include +#include + +/// +/// Define problem-specific tagging criteria +/// +/// @param i x-index +/// @param j y-index +/// @param k z-index +/// @param tag tag array (TagBox) +/// @param state simulation state (Fab) +/// @param level AMR level +/// @param geomdata geometry data +/// +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void problem_tagging(int i, int j, int k, + Array4 const& tag, + Array4 const& state, + int level, const GeometryData& geomdata) +{ + + GpuArray loc; + position(i, j, k, geomdata, loc); + + if (problem::tag_by_density) { + if (state(i,j,k,URHO) > problem::cutoff_density && + state(i,j,k,UFS) / state(i,j,k,URHO) > problem::X_min) { + + Real dist = std::abs(loc[0]); + + if (level < problem::max_hse_tagging_level && dist < problem::r_refine_distance) { + tag(i,j,k) = TagBox::SET; + } + } + + if (state(i,j,k,URHO) > problem::cutoff_density) { + if (level < problem::max_base_tagging_level) { + tag(i,j,k) = TagBox::SET; + } + } + + } else { + + // tag everything below a certain height + if (loc[0] < problem::refine_height) { + if (level < problem::max_base_tagging_level) { + tag(i,j,k) = TagBox::SET; + } + } + } + +} +#endif