change subchandra model generation to allow different core X (AMReX-Astro#3166)

this now allows Ne20 and also the C/O ratio can be something other
than 1/2.  Each core X is specified separately now.

Author: zingale

Exec/science/subchandra/_prob_params

@@ -54,7 +54,14 @@ temp_core            real         1.e7                   y
 
 temp_base            real         4.e8                   y
 
-mixed_co_wd          int          1                      y
+# mass fraction of C12 to put in the core
+core_XC12            real         0.5                    y
+
+# mass fraction of O16 to put in the core
+core_XO16            real         0.5                    y
+
+# mass fraction of Ne20 to put in the core
+core_XNe20           real         0.0                    y
 
 low_density_cutoff   real         1.e-4                  y
 

Exec/science/subchandra/initial_model.H

@@ -34,10 +34,12 @@ struct model_t {
     int ic12;
     int in14;
     int io16;
+    int ine20;
 
     amrex::Real X_N14;
     amrex::Real X_C12;
     amrex::Real X_O16;
+    amrex::Real X_Ne20;
 };
 
 
@@ -428,7 +430,8 @@ build_star(const model_t& model_params, const Real dx,
             (model::profile(0).state(i, model::ispec+model_params.ihe4) <=
              network_rp::small_x)) {
             core_X += model::profile(0).state(i, model::ispec+model_params.ic12) +
-                model::profile(0).state(i, model::ispec+model_params.io16);
+                      model::profile(0).state(i, model::ispec+model_params.io16) +
+                      model::profile(0).state(i, model::ispec+model_params.ine20);
         }
 
         mass_he += vol * model::profile(0).state(i, model::idens) * layer_X;

Exec/science/subchandra/inputs_2d.N14.coarse

@@ -118,7 +118,8 @@ problem.delta = 5.e6
 problem.temp_core = 1.e7
 problem.temp_base = 1.75e8
 
-problem.mixed_co_wd = 1
+problem.core_XC12 = 0.5
+problem.core_XO16 = 0.5
 
 problem.X_N14 = 0.01
 

Exec/science/subchandra/inputs_2d.NSE

@@ -113,7 +113,8 @@ problem.delta = 5.e6
 problem.temp_core = 1.e7
 problem.temp_base = 1.75e8
 
-problem.mixed_co_wd = 1
+problem.core_XC12 = 0.5
+problem.core_XO16 = 0.5
 
 problem.X_N14 = 0.0
 

Exec/science/subchandra/inputs_2d.ONe

@@ -0,0 +1,182 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+
+amr.plot_files_output = 1
+amr.checkpoint_files_output = 1
+
+max_step = 1000000
+stop_time = 10.0
+
+geometry.is_periodic = 0 0
+geometry.coord_sys = 1         # r-z coordinates
+
+geometry.prob_lo   =  0.    0.
+geometry.prob_hi   =  5.12e9 1.024e10
+
+amr.n_cell         = 640 1280
+
+castro.lo_bc       =  3 2
+castro.hi_bc       =  2 2
+
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+# 0 = Interior           3 = Symmetry
+# 1 = Inflow             4 = SlipWall
+# 2 = Outflow            5 = NoSlipWall
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+
+castro.do_hydro = 1
+castro.do_grav  = 1
+castro.do_react = 1
+castro.do_sponge = 1
+
+castro.react_rho_min = 100.0
+castro.react_T_min = 5.e7
+
+castro.disable_shock_burning = 1
+castro.ppm_type = 1
+castro.ppm_temp_fix = 0
+castro.use_pslope = 1
+
+castro.use_flattening = 1
+
+castro.riemann_solver = 1
+
+# Full self-gravity with the Poisson equation
+gravity.gravity_type = PoissonGrav
+
+# Multipole expansion includes terms up to r**(-max_multipole_order)
+gravity.max_multipole_order = 6
+
+# Tolerance for multigrid solver for phi solves
+gravity.abs_tol = 1.e-10
+
+# Use sync solve for gravity after refluxing
+#gravity.no_sync = 0
+
+# Disable the use of the lagged composite correction for the potential
+gravity.do_composite_phi_correction = 0
+
+castro.sponge_upper_density = 1.e4
+castro.sponge_lower_density = 1.e2
+castro.sponge_timescale     = 1.e-3
+
+castro.cfl            = 0.2     # cfl number for hyperbolic system
+castro.init_shrink    = 0.05     # scale back initial timestep by this factor
+castro.change_max     = 1.025    # factor by which dt is allowed to change each timestep
+castro.sum_interval   = 5       # timesteps between computing and printing volume averages
+castro.update_sources_after_reflux = 0
+castro.time_integration_method = 3
+
+castro.use_retry = 1
+castro.retry_subcycle_factor = 0.5
+castro.max_subcycles = 32
+
+castro.abundance_failure_rho_cutoff = 1.0
+
+amr.v               = 1       # control verbosity in Amr.cpp
+castro.v            = 1       # control verbosity in Castro.cpp
+
+castro.small_dens   = 1.e-5
+castro.small_temp   = 1.e5
+
+#castro.dtnuc_e = 0.25
+#castro.dtnuc_X = 0.25
+
+# gridding
+
+amr.max_level      = 2       # maximum level number allowed
+
+amr.ref_ratio       = 2 2 2 2 # refinement ratio
+amr.regrid_int      = 2   # how often to regrid
+amr.n_error_buf     = 2 2 2 2 # number of buffer cells in error est
+amr.grid_eff        = 0.7     # what constitutes an efficient grid
+
+amr.max_grid_size   = 256      # maximum grid size allowed -- used to control parallelism
+amr.blocking_factor = 32       # block factor in grid generation
+
+# I/O
+
+amr.check_file      = subch_chk     # root name of checkpoint file
+amr.check_int       = 50     # number of timesteps between checkpoints
+amr.plot_file       = subch_plt     # root name of plot file
+amr.plot_int        = -1     # number of timesteps between plotfiles
+amr.plot_per        = 2.e-3
+
+amr.derive_plot_vars = ALL
+castro.store_burn_weights = 1
+
+# problem parameters
+
+problem.pert_temp_factor = 20.0
+problem.pert_rad_factor = 0.5
+problem.R_pert = 1.e7
+
+problem.M_WD = 1.1
+problem.M_He  = 0.05
+
+problem.delta = 5.e6
+
+problem.temp_core = 1.e7
+problem.temp_base = 1.75e8
+
+problem.core_XC12 = 0.0
+problem.core_XO16 = 0.7
+problem.core_XNe20 = 0.3
+
+problem.X_N14 = 0.0
+
+problem.low_density_cutoff = 1.e-4
+problem.temp_fluff = 7.5e7
+
+problem.model_min_res = 5.e5
+
+# tagging
+
+amr.refinement_indicators = tempgrad denerr dencore temperr dencutoff
+
+amr.refine.tempgrad.field_name = Temp
+amr.refine.tempgrad.relative_gradient = 2.0
+amr.refine.tempgrad.max_level = 2
+
+amr.refine.denerr.field_name = density
+amr.refine.denerr.value_greater = 1.0
+amr.refine.denerr.max_level = 2
+
+# this just refines the very center to capture the convergence of the
+# compression wave
+
+amr.refine.dencore.field_name = density
+amr.refine.dencore.value_greater = 7.5e7
+amr.refine.dencore.max_level = 3
+
+# now we want to refine regions where T > 2.e8 up to the maximum
+# level, but only if rho > 1.e3 -- otherwise, we don't care about
+# following things
+
+amr.refine.temperr.field_name = Temp
+amr.refine.temperr.value_greater = 1.e8
+amr.refine.temperr.max_level = 3
+
+amr.refine.dencutoff.derefine = 1
+amr.refine.dencutoff.field_name = density
+amr.refine.dencutoff.value_less = 1.e4
+
+# Microphysics
+
+network.small_x = 1.e-10
+integrator.SMALL_X_SAFE = 1.e-10
+
+integrator.rtol_spec = 1.e-5
+integrator.atol_spec = 1.e-5
+integrator.rtol_enuc = 1.e-5
+integrator.atol_enuc = 1.e-5
+integrator.jacobian = 1
+
+integrator.X_reject_buffer = 4.0
+
+# disable jacobian caching in VODE
+integrator.use_jacobian_caching = 0
+
+integrator.ode_max_steps = 20000
+
+integrator.use_burn_retry = 1
+integrator.retry_swap_jacobian = 1

Exec/science/subchandra/inputs_3d.N14.coarse

@@ -128,7 +128,8 @@ problem.delta = 5.e6
 problem.temp_core = 1.e7
 problem.temp_base = 1.75e8
 
-problem.mixed_co_wd = 1
+problem.core_XC12 = 0.5
+problem.core_XO16 = 0.5
 
 problem.X_N14 = 0.01
 

Exec/science/subchandra/problem_initialize.H

@@ -61,31 +61,55 @@ void problem_initialize ()
     int ihe4 = network_spec_index("helium-4");
     int ic12 = network_spec_index("carbon-12");
     int io16 = network_spec_index("oxygen-16");
+    int ine20 = network_spec_index("neon-20");
 
     int in14 = network_spec_index("nitrogen-14");
 
-    if (ihe4 < 0 || ic12 < 0 || io16 < 0) {
+    if (ihe4 < 0 || ic12 < 0 || io16 < 0 || ine20 < 0) {
         amrex::Error("ERROR: species not defined");
     }
 
     if (in14 < 0 && problem::X_N14 > 0.0_rt) {
         amrex::Error("ERROR: N14 not defined");
     }
 
+    // core mass fractions should sum to 1
+    if (std::abs(problem::core_XC12 +
+                 problem::core_XO16 +
+                 problem::core_XNe20 - 1.0_rt) > 1.e-10) {
+        std::cout << problem::core_XC12 << " " << problem::core_XO16 << " " << problem::core_XNe20 << std::endl;
+        amrex::Error("ERROR: core mass fractions should sum to 1");
+    }
+
     for (auto& X : model_params.xn_core) {
         X = network_rp::small_x;
     }
     for (auto& X : model_params.xn_he) {
         X = network_rp::small_x;
     }
 
-    if (problem::mixed_co_wd) {
-        model_params.xn_core[ic12] = 0.5_rt - 0.5_rt * (NumSpec - 1) * network_rp::small_x;
-        model_params.xn_core[io16] = 0.5_rt - 0.5_rt * (NumSpec - 1) * network_rp::small_x;
-    } else {
-        model_params.xn_core[ic12] = 1.0_rt - (NumSpec - 1) * network_rp::small_x;
+    int ncore = 0;
+    if (problem::core_XC12 > 0.0) {
+        ncore++;
     }
 
+    if (problem::core_XO16 > 0.0) {
+        ncore++;
+    }
+
+    if (problem::core_XNe20 > 0.0) {
+        ncore++;
+    }
+
+    // we need to initialize the core to be 1 - other_xn * small_x.
+    // We'll distribute this remainder over the core nuclei
+
+    int other_xn = NumSpec - ncore;
+
+    model_params.xn_core[ic12] = problem::core_XC12 - other_xn * network_rp::small_x / ncore;
+    model_params.xn_core[io16] = problem::core_XO16 - other_xn * network_rp::small_x / ncore;;
+    model_params.xn_core[ine20] = problem::core_XNe20 - other_xn * network_rp::small_x / ncore;;
+
     if (problem::X_N14 > 0.0_rt) {
         model_params.xn_he[in14] = problem::X_N14;
     }
@@ -108,6 +132,7 @@ void problem_initialize ()
     model_params.ic12 = ic12;
     model_params.in14 = in14;
     model_params.io16 = io16;
+    model_params.ine20 = ine20;
 
     model_params.X_N14 = problem::X_N14;
     model_params.X_C12 = problem::X_C12;