more work on 2D

AMReX-Astro · May 30, 2024 · 608d089 · 608d089
1 parent ea687b1
commit 608d089
Show file tree

Hide file tree

Showing 4 changed files with 58 additions and 86 deletions.
diff --git a/Exec/science/circular_det/GNUmakefile b/Exec/science/circular_det/GNUmakefile
@@ -3,7 +3,7 @@ PROFILE    = FALSE
 
 DEBUG      = FALSE
 
-DIM        = 1
+DIM        = 2
 
 COMP	   = gnu
 

diff --git a/Exec/science/circular_det/inputs b/Exec/science/circular_det/inputs
@@ -1,22 +1,22 @@
 # ------------------  INPUTS TO MAIN PROGRAM  -------------------
-max_step = 15000
+max_step = 200000
 stop_time =  0.2
 
 # PROBLEM SIZE & GEOMETRY
 geometry.is_periodic = 0 0 0
 geometry.coord_sys   = 0  # 0 => cart, 1 => RZ  2=>spherical
 geometry.prob_lo     = 0     0     0
-geometry.prob_hi     = 4.e7  2500  2500
-amr.n_cell           = 256   16    16
+geometry.prob_hi     = 5.89824e7  5.89824e7
+amr.n_cell           = 384   384
 
 
 # >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
 # 0 = Interior           3 = Symmetry
 # 1 = Inflow             4 = SlipWall
 # 2 = Outflow            5 = NoSlipWall
 # >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
-castro.lo_bc       =  2   4   4
-castro.hi_bc       =  2   4   4
+castro.lo_bc       =  2   2   4
+castro.hi_bc       =  2   2   4
 
 # WHICH PHYSICS
 castro.do_hydro = 1
@@ -25,20 +25,26 @@ castro.do_react = 1
 castro.ppm_type = 1
 castro.ppm_temp_fix = 0
 
-castro.transverse_reset_density = 1
+castro.time_integration_method = 0
 
-castro.use_flattening = 1
+castro.disable_shock_burning = 0
+castro.shock_detection_threshold = 1
+
+# castro.transverse_reset_density = 1
 
+castro.small_dens   = 1.e-5
 castro.small_temp     = 1.e7
 
-castro.riemann_solver = 0
+castro.use_flattening = 1
+
+castro.riemann_solver = 1
 
 # TIME STEP CONTROL
-castro.cfl            = 0.2     # cfl number for hyperbolic system
-castro.init_shrink    = 0.01     # scale back initial timestep
+castro.cfl            = 0.5     # cfl number for hyperbolic system
+castro.init_shrink    = 0.1     # scale back initial timestep
 castro.change_max     = 1.05    # scale back initial timestep
 
-castro.dtnuc_e = 0.1
+#castro.dtnuc_e = 0.1
 
 # DIAGNOSTICS & VERBOSITY
 castro.sum_interval   = 1       # timesteps between computing mass
@@ -47,42 +53,43 @@ amr.v                 = 1       # verbosity in Amr.cpp
 #amr.grid_log        = grdlog  # name of grid logging file
 
 # REFINEMENT / REGRIDDING
-amr.max_level       = 2       # maximum level number allowed
-amr.ref_ratio       = 2 2 2 2 # refinement ratio
+amr.max_level       = 0       # maximum level number allowed
+amr.ref_ratio       = 4 4 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.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
+amr.n_error_buf     = 8 8 8 2 2 2 # number of buffer cells in error est
 
 # CHECKPOINT FILES
 amr.check_file      = det_x_chk  # root name of checkpoint file
 amr.check_int       = 1000         # number of timesteps between checkpoints
 
 # PLOTFILES
 amr.plot_file       = det_x_plt  # root name of plotfile
-amr.plot_per = 1.e-3
+amr.plot_per = 2.e-3
 amr.derive_plot_vars = ALL
 
 # problem initialization
 
-problem.T_l = 2.e9
-problem.T_r = 2.5e8
+problem.T_l = 1.1e9
+problem.T_r = 1.75e8
 
-problem.dens = 2.e6
+problem.dens = 3.e6
 problem.cfrac = 0.0
+problem.nfrac = 0.01
 
 problem.smallx = 1.e-10
 
 problem.idir = 1
 
 problem.w_T = 1.e-3
-problem.center_T = 3.e-1
+problem.center_T = 0.1
 
 # refinement
 
 amr.refinement_indicators = temperr tempgrad
 
-amr.refine.temperr.max_level = 5
+amr.refine.temperr.max_level = 0
 amr.refine.temperr.value_greater = 4.e9
 amr.refine.temperr.field_name = Temp
 
@@ -92,5 +99,15 @@ amr.refine.tempgrad.field_name = Temp
 
 # Microphysics
 
-network.use_tables = 1
+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
+
 
diff --git a/Exec/science/circular_det/problem_initialize.H b/Exec/science/circular_det/problem_initialize.H
@@ -9,6 +9,15 @@ AMREX_INLINE
 void problem_initialize ()
 {
 
+    const Geometry& dgeom = DefaultGeometry();
+
+    const Real* problo = dgeom.ProbLo();
+    const Real* probhi = dgeom.ProbHi();
+
+    for (int n = 0; n < AMREX_SPACEDIM; ++n) {
+        problem::center[n] = 0.5_rt * (problo[n] + probhi[n]);
+    }
+
     // get the species indices
 
     problem::ihe4 = network_spec_index("helium-4");
@@ -66,36 +75,6 @@ void problem_initialize ()
                                            - (NumSpec - 3) * problem::smallx;
     }
 
-    // Set the ambient material
-
-    problem::ambient_dens = problem::dens;
-    for (int n = 0; n < NumSpec; n++) {
-        problem::ambient_comp[n] = problem::xn[n];
-    }
-
-    eos_t eos_state;
-    eos_state.rho = problem::ambient_dens;
-    for (int n = 0; n < NumSpec; n++) {
-        eos_state.xn[n]  = problem::ambient_comp[n];
-    }
-
-#ifdef AUX_THERMO
-    // set the aux quantities -- we need to do this if we are using the NSE network
-    set_aux_comp_from_X(eos_state);
-#endif
-
-    eos_state.T = problem::T_l;
-
-    eos(eos_input_rt, eos_state);
-
-    problem::ambient_e_l = eos_state.e;
-
-    eos_state.T = problem::T_r;
-
-    eos(eos_input_rt, eos_state);
-
-    problem::ambient_e_r = eos_state.e;
-
 }
 
 #endif
diff --git a/Exec/science/circular_det/problem_initialize_state_data.H b/Exec/science/circular_det/problem_initialize_state_data.H
@@ -16,13 +16,17 @@ void problem_initialize_state_data (int i, int j, int k,
     const Real* probhi = geomdata.ProbHi();
 
     Real width = problem::w_T * (probhi[0] - problo[0]);
-    Real c_T = problo[0] + problem::center_T * (probhi[0] - problo[0]);
+    Real c_T = problem::center_T * (probhi[0] - problo[0]);
 
-    Real xcen = problo[0] + dx[0] * (static_cast<Real>(i) + 0.5_rt);
+    // compute distance from the center
+    Real xx = problo[0] + dx[0] * (static_cast<Real>(i) + 0.5_rt) - problem::center[0];
+    Real yy = problo[1] + dx[1] * (static_cast<Real>(j) + 0.5_rt) - problem::center[1];
+
+    Real r = std::sqrt(xx * xx + yy * yy);
 
     state(i,j,k,URHO) = problem::dens;
 
-    Real sigma = 1.0_rt / (1.0_rt + std::exp(-(c_T - xcen)/ width));
+    Real sigma = 1.0_rt / (1.0_rt + std::exp(-(c_T - r)/ width));
 
     state(i,j,k,UTEMP) = problem::T_l + (problem::T_r - problem::T_l) * (1.0_rt - sigma);
 
@@ -46,41 +50,13 @@ void problem_initialize_state_data (int i, int j, int k,
     }
 #endif
 
-#ifdef NSE_NET
-    // set initial conditions for NSE_NET
-    state(i,j,k,UMUP) = problem::mu_p;
-    state(i,j,k,UMUN) = problem::mu_n;
-
-    burn_t burn_state;
-    burn_state.T = state(i,j,k,UTEMP);
-    burn_state.rho = state(i,j,k,URHO);
-    for (int n = 0; n < NumSpec; ++n) {
-      burn_state.xn[n] = problem::xn[n];
-    }
-    // initialize y_e.
-    burn_state.y_e = 0.5_rt;
-    burn_state.mu_p = problem::mu_p;
-    burn_state.mu_n = problem::mu_n;
-
-    if (burn_state.T > 2.5e9) {
-        auto nse_state = get_actual_nse_state(burn_state, 1.0e-10_rt, true);
-
-      // for (int n = 0; n < NumSpec; ++n) {
-      //   state(i,j,k,UFS+n) = state(i,j,k,URHO) * nse_state.xn[n];
-      // }
-
-      state(i,j,k,UMUP) = nse_state.mu_p;
-      state(i,j,k,UMUN) = nse_state.mu_n;
-    }
-#endif
-
     eos(eos_input_rt, eos_state);
 
-    state(i,j,k,UMX) = state(i,j,k,URHO) * (problem::vel - 2 * problem::vel * (1.0_rt - sigma));
+    state(i,j,k,UMX) = 0.0_rt;
     state(i,j,k,UMY) = 0.0_rt;
     state(i,j,k,UMZ) = 0.0_rt;
     state(i,j,k,UEINT) = state(i,j,k,URHO) * eos_state.e;
-    state(i,j,k,UEDEN) = state(i,j,k,UEINT) + 0.5_rt * state(i,j,k,URHO) * problem::vel * problem::vel;
+    state(i,j,k,UEDEN) = state(i,j,k,UEINT);
 
 }