From aed0328358bb018bb16aaf5dd6a79cfa250bfc89 Mon Sep 17 00:00:00 2001 From: Zhi Chen <62574124+zhichen3@users.noreply.github.com> Date: Thu, 14 Nov 2024 08:04:34 -0500 Subject: [PATCH] xrb spherical problem setup (#2972) Initial problem setup for full-star xrb flame in spherical shell based on the flame_wave problem. --- Exec/science/xrb_spherical/GNUmakefile | 40 ++++ Exec/science/xrb_spherical/Make.package | 2 + Exec/science/xrb_spherical/README.md | 6 + Exec/science/xrb_spherical/_prob_params | 68 ++++++ .../xrb_spherical/analysis/r_profile.py | 56 +++++ Exec/science/xrb_spherical/analysis/slice.py | 94 ++++++++ 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 +++++ 11 files changed, 770 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/README.md create mode 100644 Exec/science/xrb_spherical/_prob_params create mode 100755 Exec/science/xrb_spherical/analysis/r_profile.py create mode 100755 Exec/science/xrb_spherical/analysis/slice.py 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..d53475be8d --- /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 = FALSE + +USE_ROTATION = FALSE +USE_DIFFUSION = FALSE + +# 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 := subch_base + +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/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 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/analysis/r_profile.py b/Exec/science/xrb_spherical/analysis/r_profile.py new file mode 100755 index 0000000000..468f2ba945 --- /dev/null +++ b/Exec/science/xrb_spherical/analysis/r_profile.py @@ -0,0 +1,56 @@ +#!/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] +print(ds.domain_left_edge[1]) +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 +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, 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(float(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, loc="lower left") + +#fig.set_size_inches(10.0, 9.0) +plt.tight_layout() +plt.savefig("{}_profiles.png".format(os.path.basename(plotfile))) diff --git a/Exec/science/xrb_spherical/analysis/slice.py b/Exec/science/xrb_spherical/analysis/slice.py new file mode 100755 index 0000000000..7e55d3848f --- /dev/null +++ b/Exec/science/xrb_spherical/analysis/slice.py @@ -0,0 +1,94 @@ +#!/usr/bin/env python3 + +import sys +import os +import yt +import numpy as np +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). +""" + +def slice(fname:str, field:str, + loc: str = "top", width_factor: float = 3.0) -> 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 = CastroDataset(fname) + 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 = width_factor * dr + box_widths = (width, width) + + 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.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)+0.5*width, r_center*np.cos(thetar))} + + # 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") + 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}") + +if __name__ == "__main__": + + if len(sys.argv) < 3: + raise Exception("Please enter parameters in order of: fname field_name width_factor[optional] loc[optional]") + + fname = sys.argv[1] + field = sys.argv[2] + loc = "top" + width_factor = 3.0 + + if len(sys.argv) == 4: + width_factor = float(sys.argv[3]) + elif len(sys.argv) > 4: + loc = sys.argv[4] + + slice(fname, field, loc=loc, width_factor=width_factor) 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..17864ed3f3 --- /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 = 1.1e6 0 +geometry.prob_hi = 1.13072e6 3.1415926 +amr.n_cell = 768 2304 #192 1152 + +# >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< +# 0 = Interior 3 = Symmetry +# 1 = Inflow 4 = SlipWall +# 2 = Outflow 5 = NoSlipWall +# >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< +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 +castro.allow_non_unit_aspect_zones = 1 + +# 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 +castro.use_pslope = 1 +castro.ppm_well_balanced = 1 +castro.pslope_cutoff_density = 1.e4 + +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 accommodate 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 = 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 +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-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 + +# problem initialization + +problem.dtemp = 1.2e9 +problem.theta_half_max = 1.745e-2 +problem.theta_half_width = 5.279e-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 = 400.0 #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.r_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..0431916a8e --- /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 (theta < problem::theta_half_max) { + f = 1.0_rt; + + } else if (theta > problem::theta_half_max + problem::theta_half_width) { + f = 0.0_rt; + + } else { + f = -(theta - 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..48b2f254c6 --- /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 < geomdata.ProbLo(0) + 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] < geomdata.ProbLo(0) + problem::refine_height) { + if (level < problem::max_base_tagging_level) { + tag(i,j,k) = TagBox::SET; + } + } + } + +} +#endif