From 3ff0599108f46815b222324323041f0710a8831e Mon Sep 17 00:00:00 2001 From: Zhi Date: Mon, 7 Oct 2024 14:42:27 -0400 Subject: [PATCH 01/25] 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 From a18761db4a5a6eaa8f395ea837f86ba51bf5d876 Mon Sep 17 00:00:00 2001 From: Zhi Date: Mon, 7 Oct 2024 16:43:35 -0400 Subject: [PATCH 02/25] update inputfile --- Exec/science/xrb_spherical/inputs.He.1000Hz | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/Exec/science/xrb_spherical/inputs.He.1000Hz b/Exec/science/xrb_spherical/inputs.He.1000Hz index 4f465e61d0..843f9f40ba 100644 --- a/Exec/science/xrb_spherical/inputs.He.1000Hz +++ b/Exec/science/xrb_spherical/inputs.He.1000Hz @@ -5,9 +5,9 @@ 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 +geometry.prob_lo = 1.1e6 0 +geometry.prob_hi = 1.13072e6 3.14159 +amr.n_cell = 192 1152 # >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< # 0 = Interior 3 = Symmetry @@ -17,11 +17,14 @@ amr.n_cell = 192 1152 castro.lo_bc = 1 3 # Inflow in lower R and Symmetry across Theta castro.hi_bc = 2 2 # Outflow in upper boundaries +# Allow non-square zones +castro.allow_non_unit_aspect_zones = 1 + # Hydrostatic condition along R-direction -castro.yl_ext_bc_type = 0 +castro.xl_ext_bc_type = 1 castro.hse_interp_temp = 0 castro.hse_fixed_temp = 1.e8 # this should match problem.T_star -castro.hse_reflect_vels = 0 +castro.hse_reflect_vels = 1 # Fill ambient states with outflow velocity in R-direction castro.fill_ambient_bc = 1 @@ -49,6 +52,7 @@ 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 @@ -136,7 +140,7 @@ problem.max_base_tagging_level = 2 problem.X_min = 1.e-2 -problem.x_refine_distance = 9.216e4 +problem.r_refine_distance = 9.216e4 # Microphysics From a7ce238629736a7bc16c0b70da97cf2849441e8e Mon Sep 17 00:00:00 2001 From: Zhi Date: Mon, 7 Oct 2024 17:12:41 -0400 Subject: [PATCH 03/25] update --- Exec/science/xrb_spherical/inputs.He.1000Hz | 2 +- Exec/science/xrb_spherical/problem_tagging.H | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Exec/science/xrb_spherical/inputs.He.1000Hz b/Exec/science/xrb_spherical/inputs.He.1000Hz index 843f9f40ba..860fb39959 100644 --- a/Exec/science/xrb_spherical/inputs.He.1000Hz +++ b/Exec/science/xrb_spherical/inputs.He.1000Hz @@ -90,7 +90,7 @@ 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.max_level = 0 #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 diff --git a/Exec/science/xrb_spherical/problem_tagging.H b/Exec/science/xrb_spherical/problem_tagging.H index 877bf1c8cc..48b2f254c6 100644 --- a/Exec/science/xrb_spherical/problem_tagging.H +++ b/Exec/science/xrb_spherical/problem_tagging.H @@ -31,7 +31,7 @@ void problem_tagging(int i, int j, int k, Real dist = std::abs(loc[0]); - if (level < problem::max_hse_tagging_level && dist < problem::r_refine_distance) { + if (level < problem::max_hse_tagging_level && dist < geomdata.ProbLo(0) + problem::r_refine_distance) { tag(i,j,k) = TagBox::SET; } } @@ -45,7 +45,7 @@ void problem_tagging(int i, int j, int k, } else { // tag everything below a certain height - if (loc[0] < problem::refine_height) { + if (loc[0] < geomdata.ProbLo(0) + problem::refine_height) { if (level < problem::max_base_tagging_level) { tag(i,j,k) = TagBox::SET; } From 6842f31fb35075666a2ac10c4fd52b64c42bef18 Mon Sep 17 00:00:00 2001 From: Zhi Date: Mon, 7 Oct 2024 17:13:31 -0400 Subject: [PATCH 04/25] update gnumake --- Exec/science/xrb_spherical/GNUmakefile | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Exec/science/xrb_spherical/GNUmakefile b/Exec/science/xrb_spherical/GNUmakefile index b6e0d4bf4a..f6cea4a9ad 100644 --- a/Exec/science/xrb_spherical/GNUmakefile +++ b/Exec/science/xrb_spherical/GNUmakefile @@ -10,10 +10,10 @@ COMP = gnu USE_MPI = TRUE USE_GRAV = TRUE -USE_REACT = TRUE +USE_REACT = FALSE -USE_ROTATION = TRUE -USE_DIFFUSION = TRUE +USE_ROTATION = FALSE +USE_DIFFUSION = FALSE # define the location of the CASTRO top directory CASTRO_HOME ?= ../../.. From 457688c3f898ff933647c283d1f9b820ab0c3242 Mon Sep 17 00:00:00 2001 From: Zhi Date: Mon, 7 Oct 2024 17:14:08 -0400 Subject: [PATCH 05/25] update --- Exec/science/xrb_spherical/inputs.He.1000Hz | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Exec/science/xrb_spherical/inputs.He.1000Hz b/Exec/science/xrb_spherical/inputs.He.1000Hz index 860fb39959..b8db7efb5f 100644 --- a/Exec/science/xrb_spherical/inputs.He.1000Hz +++ b/Exec/science/xrb_spherical/inputs.He.1000Hz @@ -113,7 +113,7 @@ amr.derive_small_plot_vars = abar x_velocity y_velocity z_velocity enuc # problem initialization -problem.dtemp = 1.2e9 +problem.dtemp = 0.0 #1.2e9 problem.theta_half_max = 1.745e-2 problem.theta_half_width = 4.9e-3 From 4169303dd80300a1c7e68289f3f8458eb33907bd Mon Sep 17 00:00:00 2001 From: Zhi Date: Tue, 8 Oct 2024 13:04:18 -0400 Subject: [PATCH 06/25] chnage input file --- Exec/science/xrb_spherical/inputs.He.1000Hz | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Exec/science/xrb_spherical/inputs.He.1000Hz b/Exec/science/xrb_spherical/inputs.He.1000Hz index b8db7efb5f..0c1569a5ab 100644 --- a/Exec/science/xrb_spherical/inputs.He.1000Hz +++ b/Exec/science/xrb_spherical/inputs.He.1000Hz @@ -124,7 +124,7 @@ problem.T_hi = 2.e8 problem.T_lo = 8.e6 problem.H_star = 2000.e0 -problem.atm_delta = 50.0 +problem.atm_delta = 400.0 #50.0 problem.fuel1_name = "helium-4" problem.fuel1_frac = 1.0e0 From d67290417c07a4ae9767da7c9afee897f802f866 Mon Sep 17 00:00:00 2001 From: Zhi Date: Wed, 9 Oct 2024 14:00:09 -0400 Subject: [PATCH 07/25] update input file --- Exec/science/xrb_spherical/inputs.He.1000Hz | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Exec/science/xrb_spherical/inputs.He.1000Hz b/Exec/science/xrb_spherical/inputs.He.1000Hz index 0c1569a5ab..72a4800ad5 100644 --- a/Exec/science/xrb_spherical/inputs.He.1000Hz +++ b/Exec/science/xrb_spherical/inputs.He.1000Hz @@ -5,8 +5,8 @@ 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 = 1.1e6 0 -geometry.prob_hi = 1.13072e6 3.14159 +geometry.prob_lo = 1.1e6 0.392699 #0.0 +geometry.prob_hi = 1.13072e6 2.748893 #3.14159 amr.n_cell = 192 1152 # >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< @@ -107,7 +107,7 @@ 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_per = 1.e-6 #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 From 7f2176649b79d80516bb118ab713d8dab67f5d78 Mon Sep 17 00:00:00 2001 From: Zhi Date: Thu, 10 Oct 2024 14:59:00 -0400 Subject: [PATCH 08/25] add a slice plot script --- Exec/science/xrb_spherical/analysis/slice.py | 75 ++++++++++++++++++++ 1 file changed, 75 insertions(+) create mode 100755 Exec/science/xrb_spherical/analysis/slice.py diff --git a/Exec/science/xrb_spherical/analysis/slice.py b/Exec/science/xrb_spherical/analysis/slice.py new file mode 100755 index 0000000000..6c38d76822 --- /dev/null +++ b/Exec/science/xrb_spherical/analysis/slice.py @@ -0,0 +1,75 @@ +#!/usr/bin/env python3 + +import sys +import os +import yt +import numpy as np +import matplotlib.pyplot as plt +from yt.units import cm + +""" +Given a plot file and field name, it gives slice plots at the top, +middle, and bottom of the domain (shell). +""" + +def slice(fname:str, field:str, loc:str="top") -> None: + """ + A slice plot of the dataset for Spherical2D geometry. + + Parameter + ======================= + fname: plot file name + field: field parameter + loc: location on the domain. {top, mid, bot} + """ + + ds = yt.load(fname, hint="castro") + currentTime = ds.current_time.in_units("s") + print(f"Current time of this plot file is {currentTime} s") + + # Some geometry properties + rr = ds.domain_right_edge[0].in_units("km") + rl = ds.domain_left_edge[0].in_units("km") + dr = rr - rl + r_center = 0.5 * (rr + rl) + + thetar = ds.domain_right_edge[1] + thetal = ds.domain_left_edge[1] + dtheta = thetar - thetal + theta_center = 0.5 * (thetar + thetal) + + # Domain width of the slice plot + width = (2.0*dr, 2.0*dr) + + loc = loc.lower() + loc_options = ["top", "mid", "bot"] + + if loc not in loc_options: + raise Exception("loc parameter must be top, mid or bot") + + # Centers for the Top, Mid and Bot panels + centers = {"top":(r_center*np.sin(thetal)+0.8*dr, r_center*np.cos(thetal)), + "mid":(r_center*np.sin(theta_center), r_center*np.cos(theta_center)), + "bot":(r_center*np.sin(thetar)+0.8*dr, r_center*np.cos(thetar))} + + sp = yt.SlicePlot(ds, 'phi', field, width=width) + + sp.set_center(centers[loc]) + + sp.set_cmap(field, "viridis") + if field in ["x_velocity", "y_velocity", "z_velocity"]: + sp.set_cmap(field, "coolwarm") + + # sp.annotate_text((0.05, 0.05), f"{currentTime.in_cgs():8.5f} s") + sp.save(f"{ds}_{loc}") + +if __name__ == "__main__": + + if len(sys.argv) < 3: + raise Exception("Please enter parameters in order of: fname field_name loc[optional]") + fname = sys.argv[1] + field = sys.argv[2] + if len(sys.argv) > 3: + loc = sys.argv[3] + + slice(fname, field, loc=loc) From 934042ec63e899b5bc41c4412a74cdd9ec983712 Mon Sep 17 00:00:00 2001 From: Zhi Date: Fri, 11 Oct 2024 15:52:50 -0400 Subject: [PATCH 09/25] use pslope --- Exec/science/xrb_spherical/inputs.He.1000Hz | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Exec/science/xrb_spherical/inputs.He.1000Hz b/Exec/science/xrb_spherical/inputs.He.1000Hz index 72a4800ad5..d55d64e14d 100644 --- a/Exec/science/xrb_spherical/inputs.He.1000Hz +++ b/Exec/science/xrb_spherical/inputs.He.1000Hz @@ -15,7 +15,7 @@ amr.n_cell = 192 1152 # 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 +castro.hi_bc = 2 3 # Outflow in upper boundaries # Allow non-square zones castro.allow_non_unit_aspect_zones = 1 @@ -43,6 +43,8 @@ castro.small_dens = 1.e-5 castro.ppm_type = 1 castro.grav_source_type = 2 +castro.use_pslope = 1 +castro.pslope_cutoff_density = 1.e4 gravity.gravity_type = ConstantGrav From 7dca9c5640cfa4536d686593d01276ae622b305d Mon Sep 17 00:00:00 2001 From: Zhi Date: Fri, 11 Oct 2024 17:27:40 -0400 Subject: [PATCH 10/25] fix --- Exec/science/xrb_spherical/problem_initialize_state_data.H | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Exec/science/xrb_spherical/problem_initialize_state_data.H b/Exec/science/xrb_spherical/problem_initialize_state_data.H index 088b1b42a9..0431916a8e 100644 --- a/Exec/science/xrb_spherical/problem_initialize_state_data.H +++ b/Exec/science/xrb_spherical/problem_initialize_state_data.H @@ -24,14 +24,14 @@ void problem_initialize_state_data (int i, int j, int k, Real f; - if (r < problem::theta_half_max) { + if (theta < problem::theta_half_max) { f = 1.0_rt; - } else if (r > problem::theta_half_max + problem::theta_half_width) { + } else if (theta > 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; + f = -(theta - problem::theta_half_max) / problem::theta_half_width + 1.0_rt; } state(i,j,k,URHO) = f * interpolate(r, model::idens, 1) + From 0a1cfb0a8a3f0272ad55d90bf80bb515e19d731e Mon Sep 17 00:00:00 2001 From: Zhi Date: Fri, 11 Oct 2024 21:20:30 -0400 Subject: [PATCH 11/25] update script --- Exec/science/xrb_spherical/analysis/slice.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Exec/science/xrb_spherical/analysis/slice.py b/Exec/science/xrb_spherical/analysis/slice.py index 6c38d76822..42e5b12a14 100755 --- a/Exec/science/xrb_spherical/analysis/slice.py +++ b/Exec/science/xrb_spherical/analysis/slice.py @@ -67,8 +67,11 @@ def slice(fname:str, field:str, loc:str="top") -> None: if len(sys.argv) < 3: raise Exception("Please enter parameters in order of: fname field_name loc[optional]") + fname = sys.argv[1] field = sys.argv[2] + loc = "top" + if len(sys.argv) > 3: loc = sys.argv[3] From d68885a1c190a10f2310b9888a33eb676336d0c7 Mon Sep 17 00:00:00 2001 From: Zhi Date: Mon, 14 Oct 2024 19:17:26 -0400 Subject: [PATCH 12/25] update script --- Exec/science/xrb_spherical/analysis/slice.py | 10 +++++----- Exec/science/xrb_spherical/inputs.He.1000Hz | 4 ++-- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/Exec/science/xrb_spherical/analysis/slice.py b/Exec/science/xrb_spherical/analysis/slice.py index 42e5b12a14..e459771a8e 100755 --- a/Exec/science/xrb_spherical/analysis/slice.py +++ b/Exec/science/xrb_spherical/analysis/slice.py @@ -5,7 +5,7 @@ import yt import numpy as np import matplotlib.pyplot as plt -from yt.units import cm +from yt.frontends.boxlib.api import CastroDataset """ Given a plot file and field name, it gives slice plots at the top, @@ -23,7 +23,7 @@ def slice(fname:str, field:str, loc:str="top") -> None: loc: location on the domain. {top, mid, bot} """ - ds = yt.load(fname, hint="castro") + ds = CastroDataset(fname) currentTime = ds.current_time.in_units("s") print(f"Current time of this plot file is {currentTime} s") @@ -39,7 +39,7 @@ def slice(fname:str, field:str, loc:str="top") -> None: theta_center = 0.5 * (thetar + thetal) # Domain width of the slice plot - width = (2.0*dr, 2.0*dr) + width = (3.0*dr, 3.0*dr) loc = loc.lower() loc_options = ["top", "mid", "bot"] @@ -48,9 +48,9 @@ def slice(fname:str, field:str, loc:str="top") -> None: raise Exception("loc parameter must be top, mid or bot") # Centers for the Top, Mid and Bot panels - centers = {"top":(r_center*np.sin(thetal)+0.8*dr, r_center*np.cos(thetal)), + centers = {"top":(r_center*np.sin(thetal)+1.48*dr, r_center*np.cos(thetal)), "mid":(r_center*np.sin(theta_center), r_center*np.cos(theta_center)), - "bot":(r_center*np.sin(thetar)+0.8*dr, r_center*np.cos(thetar))} + "bot":(r_center*np.sin(thetar)+1.48*dr, r_center*np.cos(thetar))} sp = yt.SlicePlot(ds, 'phi', field, width=width) diff --git a/Exec/science/xrb_spherical/inputs.He.1000Hz b/Exec/science/xrb_spherical/inputs.He.1000Hz index d55d64e14d..db5c4ab4a5 100644 --- a/Exec/science/xrb_spherical/inputs.He.1000Hz +++ b/Exec/science/xrb_spherical/inputs.He.1000Hz @@ -5,8 +5,8 @@ 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 = 1.1e6 0.392699 #0.0 -geometry.prob_hi = 1.13072e6 2.748893 #3.14159 +geometry.prob_lo = 1.1e6 0 +geometry.prob_hi = 1.13072e6 3.1415926 amr.n_cell = 192 1152 # >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< From c99bde87d0277b831046abe190d4bc5dc94bcc67 Mon Sep 17 00:00:00 2001 From: Zhi Date: Mon, 14 Oct 2024 21:46:30 -0400 Subject: [PATCH 13/25] update --- Exec/science/xrb_spherical/GNUmakefile | 4 ++-- Exec/science/xrb_spherical/analysis/slice.py | 19 ++++++++++++------- Exec/science/xrb_spherical/inputs.He.1000Hz | 2 +- 3 files changed, 15 insertions(+), 10 deletions(-) diff --git a/Exec/science/xrb_spherical/GNUmakefile b/Exec/science/xrb_spherical/GNUmakefile index f6cea4a9ad..f6c6800f13 100644 --- a/Exec/science/xrb_spherical/GNUmakefile +++ b/Exec/science/xrb_spherical/GNUmakefile @@ -10,7 +10,7 @@ COMP = gnu USE_MPI = TRUE USE_GRAV = TRUE -USE_REACT = FALSE +USE_REACT = TRUE USE_ROTATION = FALSE USE_DIFFUSION = FALSE @@ -26,7 +26,7 @@ NUM_MODELS := 2 EOS_DIR := helmholtz # This sets the network directory in $(MICROPHYSICS_HOME)/networks -NETWORK_DIR := aprox13 +NETWORK_DIR := subch_base INTEGRATOR_DIR := VODE diff --git a/Exec/science/xrb_spherical/analysis/slice.py b/Exec/science/xrb_spherical/analysis/slice.py index e459771a8e..cf539a4ee8 100755 --- a/Exec/science/xrb_spherical/analysis/slice.py +++ b/Exec/science/xrb_spherical/analysis/slice.py @@ -12,7 +12,8 @@ middle, and bottom of the domain (shell). """ -def slice(fname:str, field:str, loc:str="top") -> None: +def slice(fname:str, field:str, + loc: str = "top", width_factor: float = 3.0) -> None: """ A slice plot of the dataset for Spherical2D geometry. @@ -39,7 +40,8 @@ def slice(fname:str, field:str, loc:str="top") -> None: theta_center = 0.5 * (thetar + thetal) # Domain width of the slice plot - width = (3.0*dr, 3.0*dr) + width = width_factor * dr + box_widths = (width, width) loc = loc.lower() loc_options = ["top", "mid", "bot"] @@ -48,11 +50,11 @@ def slice(fname:str, field:str, loc:str="top") -> None: raise Exception("loc parameter must be top, mid or bot") # Centers for the Top, Mid and Bot panels - centers = {"top":(r_center*np.sin(thetal)+1.48*dr, r_center*np.cos(thetal)), + centers = {"top":(r_center*np.sin(thetal)+0.5*width, r_center*np.cos(thetal)), "mid":(r_center*np.sin(theta_center), r_center*np.cos(theta_center)), - "bot":(r_center*np.sin(thetar)+1.48*dr, r_center*np.cos(thetar))} + "bot":(r_center*np.sin(thetar)+0.5*width, r_center*np.cos(thetar))} - sp = yt.SlicePlot(ds, 'phi', field, width=width) + sp = yt.SlicePlot(ds, 'phi', field, width=box_widths) sp.set_center(centers[loc]) @@ -71,8 +73,11 @@ def slice(fname:str, field:str, loc:str="top") -> None: fname = sys.argv[1] field = sys.argv[2] loc = "top" + width_factor = 3.0 - if len(sys.argv) > 3: + if len(sys.argv) == 4: loc = sys.argv[3] + elif len(sys.argv) > 4: + width_factor = float(sys.argv[4]) - slice(fname, field, loc=loc) + slice(fname, field, loc=loc, width_factor=width_factor) diff --git a/Exec/science/xrb_spherical/inputs.He.1000Hz b/Exec/science/xrb_spherical/inputs.He.1000Hz index db5c4ab4a5..ca42ec3bbd 100644 --- a/Exec/science/xrb_spherical/inputs.He.1000Hz +++ b/Exec/science/xrb_spherical/inputs.He.1000Hz @@ -117,7 +117,7 @@ amr.derive_small_plot_vars = abar x_velocity y_velocity z_velocity enuc problem.dtemp = 0.0 #1.2e9 problem.theta_half_max = 1.745e-2 -problem.theta_half_width = 4.9e-3 +problem.theta_half_width = 5.279e-3 problem.dens_base = 3.43e6 From eecbdf57efb686f2cab9fff698653234459c8410 Mon Sep 17 00:00:00 2001 From: Zhi Date: Mon, 14 Oct 2024 22:41:31 -0400 Subject: [PATCH 14/25] fix spelling --- Exec/science/xrb_spherical/inputs.He.1000Hz | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Exec/science/xrb_spherical/inputs.He.1000Hz b/Exec/science/xrb_spherical/inputs.He.1000Hz index ca42ec3bbd..3f9661ccf7 100644 --- a/Exec/science/xrb_spherical/inputs.He.1000Hz +++ b/Exec/science/xrb_spherical/inputs.He.1000Hz @@ -55,7 +55,7 @@ gravity.const_grav = -1.5e14 # 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 +# Centrifugal is not important since NS would simply deform to accommodate for it castro.rotation_include_centrifugal = 0 castro.diffuse_temp = 1 @@ -115,7 +115,7 @@ amr.derive_small_plot_vars = abar x_velocity y_velocity z_velocity enuc # problem initialization -problem.dtemp = 0.0 #1.2e9 +problem.dtemp = 1.2e9 problem.theta_half_max = 1.745e-2 problem.theta_half_width = 5.279e-3 From 3c63920d958681897c794076e8108aee8e3a0330 Mon Sep 17 00:00:00 2001 From: Zhi Date: Tue, 15 Oct 2024 13:50:09 -0400 Subject: [PATCH 15/25] update colorbar limit --- Exec/science/xrb_spherical/analysis/slice.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/Exec/science/xrb_spherical/analysis/slice.py b/Exec/science/xrb_spherical/analysis/slice.py index cf539a4ee8..9c7b3676c6 100755 --- a/Exec/science/xrb_spherical/analysis/slice.py +++ b/Exec/science/xrb_spherical/analysis/slice.py @@ -61,6 +61,13 @@ def slice(fname:str, field:str, sp.set_cmap(field, "viridis") if field in ["x_velocity", "y_velocity", "z_velocity"]: sp.set_cmap(field, "coolwarm") + elif field == "Temp": + sp.set_zlim(f, 5.e7, 2.5e9) + sp.set_cmap(f, "magma_r") + elif field == "enuc": + sp.set_zlim(f, 1.e18, 1.e20) + elif field == "density": + sp.set_zlim(f, 1.e-3, 5.e8) # sp.annotate_text((0.05, 0.05), f"{currentTime.in_cgs():8.5f} s") sp.save(f"{ds}_{loc}") From 15ed991bfa48c46478026447cf724e46c1a301d8 Mon Sep 17 00:00:00 2001 From: Zhi Date: Wed, 16 Oct 2024 12:00:47 -0400 Subject: [PATCH 16/25] add README --- Exec/science/xrb_spherical/README.md | 6 ++++++ 1 file changed, 6 insertions(+) create mode 100644 Exec/science/xrb_spherical/README.md diff --git a/Exec/science/xrb_spherical/README.md b/Exec/science/xrb_spherical/README.md new file mode 100644 index 0000000000..325667c9c1 --- /dev/null +++ b/Exec/science/xrb_spherical/README.md @@ -0,0 +1,6 @@ +# xrb_spherical + +This is the full-star XRB flame setup based on flame_wave. +This setup uses a spherical 2D geometry to model XRB flame +on a spherical shell with initial temperature perturbation +on the north pole. \ No newline at end of file From 321fd689f25ed40b1c1b55c8f1df62b789fa44b3 Mon Sep 17 00:00:00 2001 From: Zhi Date: Wed, 16 Oct 2024 12:06:40 -0400 Subject: [PATCH 17/25] update script --- Exec/science/xrb_spherical/analysis/slice.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Exec/science/xrb_spherical/analysis/slice.py b/Exec/science/xrb_spherical/analysis/slice.py index 9c7b3676c6..5ac7905f39 100755 --- a/Exec/science/xrb_spherical/analysis/slice.py +++ b/Exec/science/xrb_spherical/analysis/slice.py @@ -75,7 +75,7 @@ def slice(fname:str, field:str, if __name__ == "__main__": if len(sys.argv) < 3: - raise Exception("Please enter parameters in order of: fname field_name loc[optional]") + raise Exception("Please enter parameters in order of: fname field_name width_factor[optional] loc[optional]") fname = sys.argv[1] field = sys.argv[2] @@ -83,8 +83,8 @@ def slice(fname:str, field:str, width_factor = 3.0 if len(sys.argv) == 4: - loc = sys.argv[3] + width_factor = float(sys.argv[3]) elif len(sys.argv) > 4: - width_factor = float(sys.argv[4]) + loc = sys.argv[4] slice(fname, field, loc=loc, width_factor=width_factor) From 8eb00d666330cc4d18d1d9c03a6591847b4c5ba2 Mon Sep 17 00:00:00 2001 From: zhi Date: Thu, 24 Oct 2024 14:42:37 -0400 Subject: [PATCH 18/25] update input --- Exec/science/xrb_spherical/inputs.He.1000Hz | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Exec/science/xrb_spherical/inputs.He.1000Hz b/Exec/science/xrb_spherical/inputs.He.1000Hz index 3f9661ccf7..356cc6299d 100644 --- a/Exec/science/xrb_spherical/inputs.He.1000Hz +++ b/Exec/science/xrb_spherical/inputs.He.1000Hz @@ -14,7 +14,7 @@ amr.n_cell = 192 1152 # 1 = Inflow 4 = SlipWall # 2 = Outflow 5 = NoSlipWall # >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< -castro.lo_bc = 1 3 # Inflow in lower R and Symmetry across Theta +castro.lo_bc = 3 3 # Inflow in lower R and Symmetry across Theta castro.hi_bc = 2 3 # Outflow in upper boundaries # Allow non-square zones @@ -44,6 +44,7 @@ castro.small_dens = 1.e-5 castro.ppm_type = 1 castro.grav_source_type = 2 castro.use_pslope = 1 +castro.ppm_well_balanced = 1 castro.pslope_cutoff_density = 1.e4 gravity.gravity_type = ConstantGrav From 4337b0daf26d895368a252289aab124d49b0fbb6 Mon Sep 17 00:00:00 2001 From: zhi Date: Thu, 24 Oct 2024 15:33:54 -0400 Subject: [PATCH 19/25] fix boundary condition in input --- Exec/science/xrb_spherical/inputs.He.1000Hz | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/Exec/science/xrb_spherical/inputs.He.1000Hz b/Exec/science/xrb_spherical/inputs.He.1000Hz index 356cc6299d..b6cc876ee5 100644 --- a/Exec/science/xrb_spherical/inputs.He.1000Hz +++ b/Exec/science/xrb_spherical/inputs.He.1000Hz @@ -7,7 +7,7 @@ geometry.is_periodic = 0 0 geometry.coord_sys = 2 # 0 => cart, 1 => RZ 2=>spherical geometry.prob_lo = 1.1e6 0 geometry.prob_hi = 1.13072e6 3.1415926 -amr.n_cell = 192 1152 +amr.n_cell = 1536 2304 #192 1152 # >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< # 0 = Interior 3 = Symmetry @@ -20,12 +20,6 @@ castro.hi_bc = 2 3 # Outflow in upper boundaries # Allow non-square zones castro.allow_non_unit_aspect_zones = 1 -# Hydrostatic condition along R-direction -castro.xl_ext_bc_type = 1 -castro.hse_interp_temp = 0 -castro.hse_fixed_temp = 1.e8 # this should match problem.T_star -castro.hse_reflect_vels = 1 - # Fill ambient states with outflow velocity in R-direction castro.fill_ambient_bc = 1 castro.ambient_fill_dir = 0 @@ -145,7 +139,6 @@ problem.X_min = 1.e-2 problem.r_refine_distance = 9.216e4 - # Microphysics integrator.rtol_spec = 1.e-6 From e155b99e2b4a964372abe8a3ff001cf74f3fa288 Mon Sep 17 00:00:00 2001 From: zhi Date: Thu, 24 Oct 2024 16:59:21 -0400 Subject: [PATCH 20/25] add 1 r-profile plot at given theta --- .../xrb_spherical/analysis/r_profile.py | 54 +++++++++++++++++++ 1 file changed, 54 insertions(+) create mode 100755 Exec/science/xrb_spherical/analysis/r_profile.py diff --git a/Exec/science/xrb_spherical/analysis/r_profile.py b/Exec/science/xrb_spherical/analysis/r_profile.py new file mode 100755 index 0000000000..14fbfe1058 --- /dev/null +++ b/Exec/science/xrb_spherical/analysis/r_profile.py @@ -0,0 +1,54 @@ +#!/usr/bin/env python3 + +# Spherical R profile at different theta + +import os +import sys +import yt +import matplotlib.pyplot as plt +import numpy as np +from functools import reduce +import itertools + +import matplotlib.ticker as ptick +from yt.frontends.boxlib.api import CastroDataset +from yt.units import cm + + +plotfile = sys.argv[1] +ds = CastroDataset(plotfile) + +rmin = ds.domain_left_edge[0] +rmax = rmin + 5000.0*cm +#rmax = ds.domain_right_edge[0] + +fig, _ax = plt.subplots(2,2) + +axes = list(itertools.chain(*_ax)) + +fig.set_size_inches(7.0, 8.0) + +fields = ["Temp", "density", "x_velocity", "y_velocity"] +nice_names = [r"$T$ (K)", r"$\rho$ (g/${cm}^3$)", r"$u$ (cm/s)", r"$v$ (cm/s)"] + +# 4 rays at different theta values +thetas = [0, 0.125*np.pi, 0.25*np.pi, 0.5*np.pi] + +for i, f in enumerate(fields): + + for theta in thetas: + ray = ds.ray((rmin*np.sin(theta), rmin*np.cos(theta), 0*cm), + (rmax*np.sin(theta), rmax*np.cos(theta), 0*cm)) + isrt = np.argsort(ray["t"]) + axes[i].plot(ray['r'][isrt], ray[f][isrt], label="theta = {}".format(theta)) + + axes[i].set_xlabel(r"$r$ (cm)") + axes[i].set_ylabel(nice_names[i]) + axes[i].set_yscale("symlog") + + if i == 0: + axes[0].legend(frameon=False) + +#fig.set_size_inches(10.0, 9.0) +plt.tight_layout() +plt.savefig("{}_profiles.png".format(os.path.basename(plotfile))) From 4ddef2cff9b3792542dcc635280a763dc8fd45d8 Mon Sep 17 00:00:00 2001 From: zhi Date: Thu, 24 Oct 2024 17:24:59 -0400 Subject: [PATCH 21/25] update script --- Exec/science/xrb_spherical/analysis/r_profile.py | 4 ++-- Exec/science/xrb_spherical/inputs.He.1000Hz | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Exec/science/xrb_spherical/analysis/r_profile.py b/Exec/science/xrb_spherical/analysis/r_profile.py index 14fbfe1058..549fe77b45 100755 --- a/Exec/science/xrb_spherical/analysis/r_profile.py +++ b/Exec/science/xrb_spherical/analysis/r_profile.py @@ -40,14 +40,14 @@ ray = ds.ray((rmin*np.sin(theta), rmin*np.cos(theta), 0*cm), (rmax*np.sin(theta), rmax*np.cos(theta), 0*cm)) isrt = np.argsort(ray["t"]) - axes[i].plot(ray['r'][isrt], ray[f][isrt], label="theta = {}".format(theta)) + axes[i].plot(ray['r'][isrt], ray[f][isrt], label=r"$\theta$ = {:.4f}".format(theta)) axes[i].set_xlabel(r"$r$ (cm)") axes[i].set_ylabel(nice_names[i]) axes[i].set_yscale("symlog") if i == 0: - axes[0].legend(frameon=False) + axes[0].legend(frameon=False, loc="lower left") #fig.set_size_inches(10.0, 9.0) plt.tight_layout() diff --git a/Exec/science/xrb_spherical/inputs.He.1000Hz b/Exec/science/xrb_spherical/inputs.He.1000Hz index b6cc876ee5..17864ed3f3 100644 --- a/Exec/science/xrb_spherical/inputs.He.1000Hz +++ b/Exec/science/xrb_spherical/inputs.He.1000Hz @@ -7,7 +7,7 @@ geometry.is_periodic = 0 0 geometry.coord_sys = 2 # 0 => cart, 1 => RZ 2=>spherical geometry.prob_lo = 1.1e6 0 geometry.prob_hi = 1.13072e6 3.1415926 -amr.n_cell = 1536 2304 #192 1152 +amr.n_cell = 768 2304 #192 1152 # >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< # 0 = Interior 3 = Symmetry From 1d5ff06982afce16bb7ca1188643cfd42c9aada8 Mon Sep 17 00:00:00 2001 From: zhi Date: Tue, 29 Oct 2024 18:21:27 -0400 Subject: [PATCH 22/25] update script --- Exec/science/xrb_spherical/analysis/r_profile.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/Exec/science/xrb_spherical/analysis/r_profile.py b/Exec/science/xrb_spherical/analysis/r_profile.py index 549fe77b45..b03ff8bb85 100755 --- a/Exec/science/xrb_spherical/analysis/r_profile.py +++ b/Exec/science/xrb_spherical/analysis/r_profile.py @@ -21,7 +21,7 @@ rmin = ds.domain_left_edge[0] rmax = rmin + 5000.0*cm #rmax = ds.domain_right_edge[0] - +print(ds.domain_left_edge[1]) fig, _ax = plt.subplots(2,2) axes = list(itertools.chain(*_ax)) @@ -32,13 +32,14 @@ nice_names = [r"$T$ (K)", r"$\rho$ (g/${cm}^3$)", r"$u$ (cm/s)", r"$v$ (cm/s)"] # 4 rays at different theta values -thetas = [0, 0.125*np.pi, 0.25*np.pi, 0.5*np.pi] +thetas = [0, 0.0155625*np.pi, 0.03125*np.pi] for i, f in enumerate(fields): for theta in thetas: - ray = ds.ray((rmin*np.sin(theta), rmin*np.cos(theta), 0*cm), - (rmax*np.sin(theta), rmax*np.cos(theta), 0*cm)) + # simply go from (rmin, theta) -> (rmax, theta). Doesn't need to convert to physical R-Z + ray = ds.ray((rmin, theta*cm, 0*cm), + (rmax, theta*cm, 0*cm)) isrt = np.argsort(ray["t"]) axes[i].plot(ray['r'][isrt], ray[f][isrt], label=r"$\theta$ = {:.4f}".format(theta)) From 5a13fc5178766ec966562228e07794003d3dae5a Mon Sep 17 00:00:00 2001 From: zhi Date: Tue, 29 Oct 2024 18:33:20 -0400 Subject: [PATCH 23/25] add comment about making sliceplots --- Exec/science/xrb_spherical/analysis/slice.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/Exec/science/xrb_spherical/analysis/slice.py b/Exec/science/xrb_spherical/analysis/slice.py index 5ac7905f39..7e55d3848f 100755 --- a/Exec/science/xrb_spherical/analysis/slice.py +++ b/Exec/science/xrb_spherical/analysis/slice.py @@ -7,6 +7,8 @@ import matplotlib.pyplot as plt from yt.frontends.boxlib.api import CastroDataset +from yt.units import cm + """ Given a plot file and field name, it gives slice plots at the top, middle, and bottom of the domain (shell). @@ -54,8 +56,10 @@ def slice(fname:str, field:str, "mid":(r_center*np.sin(theta_center), r_center*np.cos(theta_center)), "bot":(r_center*np.sin(thetar)+0.5*width, r_center*np.cos(thetar))} - sp = yt.SlicePlot(ds, 'phi', field, width=box_widths) + # Note we can also set center during SlicePlot, however then we would enter in [r_center, theta_center, 0] + # rather than the physical R-Z coordinate if we do it via sp.set_center + sp = yt.SlicePlot(ds, 'phi', field, width=box_widths) sp.set_center(centers[loc]) sp.set_cmap(field, "viridis") From 4f53b47bf3bff26c2798304825673c3f4934d269 Mon Sep 17 00:00:00 2001 From: zhi Date: Wed, 30 Oct 2024 15:09:30 -0400 Subject: [PATCH 24/25] update script --- Exec/science/xrb_spherical/analysis/r_profile.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/Exec/science/xrb_spherical/analysis/r_profile.py b/Exec/science/xrb_spherical/analysis/r_profile.py index b03ff8bb85..468f2ba945 100755 --- a/Exec/science/xrb_spherical/analysis/r_profile.py +++ b/Exec/science/xrb_spherical/analysis/r_profile.py @@ -32,16 +32,17 @@ nice_names = [r"$T$ (K)", r"$\rho$ (g/${cm}^3$)", r"$u$ (cm/s)", r"$v$ (cm/s)"] # 4 rays at different theta values -thetas = [0, 0.0155625*np.pi, 0.03125*np.pi] +thetal = ds.domain_left_edge[1] +thetar = ds.domain_right_edge[1] +thetas = [thetal, 0.25*thetar, 0.5*thetar, 0.75*thetar] for i, f in enumerate(fields): for theta in thetas: # simply go from (rmin, theta) -> (rmax, theta). Doesn't need to convert to physical R-Z - ray = ds.ray((rmin, theta*cm, 0*cm), - (rmax, theta*cm, 0*cm)) + ray = ds.ray((rmin, theta, 0*cm), (rmax, theta, 0*cm)) isrt = np.argsort(ray["t"]) - axes[i].plot(ray['r'][isrt], ray[f][isrt], label=r"$\theta$ = {:.4f}".format(theta)) + axes[i].plot(ray['r'][isrt], ray[f][isrt], label=r"$\theta$ = {:.4f}".format(float(theta))) axes[i].set_xlabel(r"$r$ (cm)") axes[i].set_ylabel(nice_names[i]) From 96a07fa1c445d26537121c7d994f7c7f36f11f64 Mon Sep 17 00:00:00 2001 From: Zhi Date: Wed, 13 Nov 2024 13:53:11 -0500 Subject: [PATCH 25/25] update makefile --- Exec/science/xrb_spherical/GNUmakefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Exec/science/xrb_spherical/GNUmakefile b/Exec/science/xrb_spherical/GNUmakefile index f6c6800f13..d53475be8d 100644 --- a/Exec/science/xrb_spherical/GNUmakefile +++ b/Exec/science/xrb_spherical/GNUmakefile @@ -10,7 +10,7 @@ COMP = gnu USE_MPI = TRUE USE_GRAV = TRUE -USE_REACT = TRUE +USE_REACT = FALSE USE_ROTATION = FALSE USE_DIFFUSION = FALSE