Create problem to study He flame movement in a double det

Exec/science/subch_planar/GNUmakefile

@@ -0,0 +1,31 @@
+PRECISION        = DOUBLE
+PROFILE          = FALSE
+DEBUG            = FALSE
+DIM              = 2
+
+COMP	         = gnu
+
+USE_MPI          = TRUE
+USE_OMP          = FALSE
+
+USE_GRAV         = TRUE
+USE_REACT        = TRUE
+
+USE_SHOCK_VAR = TRUE
+
+USE_MODEL_PARSER = TRUE
+
+CASTRO_HOME ?= ../../..
+
+# This sets the EOS directory in $(MICROPHYSICS_HOME)/eos
+EOS_DIR     := helmholtz
+
+# This sets the EOS directory in $(MICROPHYSICS_HOME)/networks
+NETWORK_DIR := subch_simple
+
+PROBLEM_DIR ?= ./
+
+Bpack   := $(PROBLEM_DIR)/Make.package
+Blocs   := $(PROBLEM_DIR)
+
+include $(CASTRO_HOME)/Exec/Make.Castro

Exec/science/subch_planar/Make.package

@@ -0,0 +1 @@
+

Exec/science/subch_planar/Problem.H

@@ -0,0 +1,2 @@
+amrex::Real maxVal (const std::string& name, amrex::Real time);
+

Exec/science/subch_planar/README

@@ -0,0 +1 @@
+We want to study the impact of shock detection and burning treatment in the context of a double detonation by simulating a lateral flame in the accreted He layer.

Exec/science/subch_planar/_prob_params

@@ -0,0 +1,19 @@
+H_min                 real           1.e-4_rt      y
+
+cutoff_density        real           50.0_rt       y
+
+model_name            string         ""            y
+
+apply_vel_field       integer        0             y
+
+velpert_scale         real           1.0e2_rt      y
+
+velpert_amplitude     real           1.0e2_rt      y
+
+velpert_height_loc    real           6.5e3_rt      y
+
+num_vortices          integer        16            y
+
+max_num_vortices      integer        16            n
+
+xloc_vortices         real           0.0_rt        n        16

Exec/science/subch_planar/inputs_2d

@@ -0,0 +1,85 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 1000000
+stop_time =  1000.0
+
+# PROBLEM SIZE & GEOMETRY
+geometry.is_periodic = 0       0
+geometry.coord_sys   = 1                  # 0 => cart, 1 => RZ  2=>spherical
+geometry.prob_lo     = 0.0     0.0
+geometry.prob_hi     = 1.6e9   3.2e9
+amr.n_cell           = 1024    2048
+
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+# 0 = Interior           3 = Symmetry
+# 1 = Inflow             4 = SlipWall
+# 2 = Outflow            5 = NoSlipWall
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+castro.lo_bc       =  0   3
+castro.hi_bc       =  0   2
+
+castro.fill_ambient_bc = 1
+castro.ambient_fill_dir = 1
+castro.ambient_outflow_vel = 1
+
+
+# WHICH PHYSICS
+castro.do_hydro = 1
+castro.do_react = 1
+castro.add_ext_src   = 0
+castro.do_grav = 1
+castro.do_sponge = 1
+
+castro.ppm_type = 1
+castro.grav_source_type = 2
+castro.use_pslope = 1
+castro.pslope_cutoff_density = 1.e4
+
+gravity.gravity_type = ConstantGrav
+gravity.const_grav   = -1.466e9
+
+# TIME STEP CONTROL
+castro.cfl            = 0.7     # cfl number for hyperbolic system
+castro.init_shrink    = 0.1     # scale back initial timestep
+castro.change_max     = 1.1     # max time step growth
+
+# SPONGE
+castro.sponge_upper_density = 5.0e4
+castro.sponge_lower_density = 8.0e1
+castro.sponge_timescale     = 1.0e-3
+
+# DIAGNOSTICS & VERBOSITY
+castro.sum_interval   = 1       # timesteps between computing mass
+castro.v              = 1       # verbosity in Castro.cpp
+amr.v                 = 1       # verbosity in Amr.cpp
+
+# REFINEMENT / REGRIDDING 
+amr.max_level       = 0       # maximum level number allowed
+amr.ref_ratio       = 2 2 2 2 # refinement ratio
+amr.regrid_int      = 2 2 2 2 # how often to regrid
+amr.blocking_factor = 8       # block factor in grid generation
+amr.max_grid_size   = 64 
+amr.n_error_buf     = 2 2 2 2 # number of buffer cells in error est
+
+# CHECKPOINT FILES
+amr.check_file      = planar_chk    # root name of checkpoint file
+amr.check_int       = 100       # number of timesteps between checkpoints
+
+# PLOTFILES
+amr.plot_file       = planar_plt    # root name of plotfile
+amr.plot_int        = 50       # number of timesteps between plotfiles
+amr.derive_plot_vars = ALL
+
+# PROBLEM PARAMETERS
+problem.model_name =  "toy_subch.hse.tanh.delta_50.00km.dx_73.24km"
+
+problem.apply_vel_field = 1
+problem.velpert_height_loc = 1475.0
+problem.velpert_scale = 1.e2
+problem.velpert_amplitude = 1.e5
+problem.num_vortices = 16
+
+# MICROPHYSICS
+integrator.jacobian = 1
+
+integrator.atol_spec = 1.e-6
+integrator.rtol_spec = 1.e-6

Exec/science/subch_planar/problem_initialize.H

@@ -0,0 +1,82 @@
+#ifndef problem_initialize_H
+#define problem_initialize_H
+
+#include <prob_parameters.H>
+#include <eos.H>
+#include <model_parser.H>
+#include <ambient.H>
+
+AMREX_INLINE
+void problem_initialize ()
+{
+
+    const Geometry& dgeom = DefaultGeometry();
+
+    const Real* problo = dgeom.ProbLo();
+    const Real* probhi = dgeom.ProbHi();
+
+    if (problem::num_vortices > problem::max_num_vortices) {
+        amrex::Error("num_vortices too large, please increase max_num_vortices and the size of xloc_vortices");
+    }
+
+    // Read initial model
+    read_model_file(problem::model_name);
+
+    if (ParallelDescriptor::IOProcessor()) {
+        for (int i = 0; i < model::npts; i++) {
+            std::cout << i << " " << model::profile(0).r(i) << " " << model::profile(0).state(i,model::idens) << std::endl;
+        }
+    }
+
+    // velocity perturbation stuff
+
+    Real offset = (probhi[0] - problo[0]) / problem::num_vortices;
+
+    for (int i = 0; i < problem::num_vortices; i++) {
+        problem::xloc_vortices[i] = (static_cast<Real>(i) + 0.5_rt) * offset + problo[0];
+    }
+
+    // Set up Castro data logs for this problem
+
+    if (castro::sum_interval > 0 && amrex::ParallelDescriptor::IOProcessor()) {
+
+        Castro::problem_data_logs.resize(1);
+
+        Castro::problem_data_logs[0].reset(new std::fstream);
+        Castro::problem_data_logs[0]->open("xrb_diag.out", std::ios::out | std::ios::app);
+        if (!Castro::problem_data_logs[0]->good()) {
+            amrex::FileOpenFailed("xrb_diag.out");
+        }
+
+    }
+
+    // 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
+

Exec/science/subch_planar/problem_initialize_state_data.H

@@ -0,0 +1,116 @@
+#ifndef problem_initialize_state_data_H
+#define problem_initialize_state_data_H
+
+#include <prob_parameters.H>
+#include <eos.H>
+#include <model_parser.H>
+
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void problem_initialize_state_data (int i, int j, int k,
+                                    Array4<Real> const& state,
+                                    const GeometryData& geomdata)
+{
+
+    const Real* dx = geomdata.CellSize();
+    const Real* problo = geomdata.ProbLo();
+
+    Real x = problo[0] + dx[0] * (static_cast<Real>(i) + 0.5_rt);
+
+    Real y = 0.0;
+#if AMREX_SPACEDIM >= 2
+    y = problo[1] + dx[1] * (static_cast<Real>(j) + 0.5_rt);
+#endif
+
+    Real z = 0.0;
+#if AMREX_SPACEDIM == 3
+    z = problo[2] + dx[2] * (static_cast<Real>(k) + 0.5_rt);
+#endif
+
+#if AMREX_SPACEDIM == 2
+    Real height = y;
+#else
+    Real height = z;
+#endif
+
+    state(i,j,k,URHO) = interpolate(height, model::idens);
+    state(i,j,k,UTEMP) = interpolate(height, model::itemp);
+    for (int n = 0; n < NumSpec; n++) {
+        state(i,j,k,UFS+n) = interpolate(height, model::ispec+n);
+    }
+
+    eos_t eos_state;
+    eos_state.rho = state(i,j,k,URHO);
+    eos_state.T = state(i,j,k,UTEMP);
+    for (int n = 0; n < NumSpec; n++) {
+        eos_state.xn[n] = state(i,j,k,UFS+n);
+    }
+
+    eos(eos_input_rt, eos_state);
+
+    state(i,j,k,UEINT) = state(i,j,k,URHO) * eos_state.e;
+    state(i,j,k,UEDEN) = state(i,j,k,UEINT);
+
+    for (int n = 0; n < NumSpec; n++) {
+        state(i,j,k,UFS+n) = state(i,j,k,URHO) * state(i,j,k,UFS+n);
+    }
+
+    // Initial velocities = 0
+    state(i,j,k,UMX) = 0.0_rt;
+    state(i,j,k,UMY) = 0.0_rt;
+    state(i,j,k,UMZ) = 0.0_rt;
+
+    //  Now add the velocity perturbation
+    if (problem::apply_vel_field) {
+
+        Real zdist = z - problem::velpert_height_loc;
+        Real ydist = y - problem::velpert_height_loc;
+
+        Real upert[3] = {0.0_rt};
+
+        // at the moment this is really just a 2-d perturbation
+
+        // loop over each vortex
+
+        for (int vortex = 0; vortex < problem::num_vortices; vortex++) {
+
+            Real xdist = x - problem::xloc_vortices[vortex];
+
+#if AMREX_SPACEDIM == 2
+            Real r = std::sqrt(xdist * xdist + ydist * ydist);
+
+            upert[0] += -ydist * problem::velpert_amplitude *
+                std::exp(-r * r / (2.0_rt * problem::velpert_scale)) *
+                std::pow(-1.0_rt, vortex+1);
+
+            upert[1] += xdist * problem::velpert_amplitude *
+                std::exp(-r * r / (2.0_rt * problem::velpert_scale)) *
+                std::pow(-1.0_rt, vortex+1);
+
+            upert[2] = 0.0_rt;
+#else
+            Real r = std::sqrt(xdist * xdist + zdist * zdist);
+
+            upert[0] += -zdist * problem::velpert_amplitude *
+                std::exp(-r * r / (2.0_rt * problem::velpert_scale)) *
+                std::pow(-1.0_rt, vortex+1);
+
+            upert[2] += xdist * problem::velpert_amplitude *
+                std::exp(-r * r / (2.0_rt * problem::velpert_scale)) *
+                std::pow(-1.0_rt, vortex+1);
+
+            upert[1] = 0.0_rt;
+#endif
+
+        }
+
+        state(i,j,k,UMX) = state(i,j,k,URHO) * upert[0];
+        state(i,j,k,UMY) = state(i,j,k,URHO) * upert[1];
+        state(i,j,k,UMZ) = state(i,j,k,URHO) * upert[2];
+
+        state(i,j,k,UEDEN) += 0.5_rt * state(i,j,k,URHO) * (upert[0] * upert[0] +
+                                                            upert[1] * upert[1] +
+                                                            upert[2] * upert[2]);
+    }
+}
+#endif
+