Radau quadrature for SDC

Exec/science/Detonation/inputs-det-x.sdc2

@@ -0,0 +1,77 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 15000
+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.e4  2500  2500
+amr.n_cell           = 512   16    16
+
+
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+# 0 = Interior           3 = Symmetry
+# 1 = Inflow             4 = SlipWall
+# 2 = Outflow            5 = NoSlipWall
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+castro.lo_bc       =  3   4   4
+castro.hi_bc       =  2   4   4
+
+# WHICH PHYSICS
+castro.do_hydro = 1
+castro.do_react = 1
+
+castro.time_integration_method = 2
+castro.sdc_order = 2
+castro.sdc_quadrature = 0
+#castro.sdc_extra = 1
+
+castro.sdc_solve_for_rhoe  = 1
+
+castro.sdc_solver = 1
+castro.sdc_solver_tol_dens = 1.e-5
+castro.sdc_solver_tol_spec = 1.e-5
+castro.sdc_solver_tol_ener = 1.e-5
+castro.sdc_use_analytic_jac = 1
+
+castro.ppm_type = 0
+
+castro.use_flattening = 1
+
+castro.small_temp     = 1.e7
+
+castro.riemann_solver = 0
+
+# TIME STEP CONTROL
+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.dt_cutoff      = 5.e-20  # level 0 timestep below which we halt
+castro.dtnuc_e = 0.25
+
+# DIAGNOSTICS & VERBOSITY
+castro.sum_interval   = 1       # timesteps between computing mass
+castro.v              = 1       # verbosity in Castro.cpp
+amr.v                 = 1       # verbosity in Amr.cpp
+#amr.grid_log        = grdlog  # name of grid logging file
+
+# 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 = 4       # 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      = det_x_chk  # root name of checkpoint file
+amr.check_int       = 300         # number of timesteps between checkpoints
+
+# PLOTFILES
+amr.plot_file       = det_x_plt  # root name of plotfile
+amr.plot_per = 1.e-7
+amr.derive_plot_vars = ALL
+
+#PROBIN FILENAME
+amr.probin_file = probin-det-x

Exec/science/Detonation/sdc_tests/inputs.1d.sdc

@@ -1,6 +1,6 @@
 # ------------------  INPUTS TO MAIN PROGRAM  -------------------
-max_step = 40000
-stop_time =  0.2
+max_step = 400000
+stop_time =  3.e-6
 
 # PROBLEM SIZE & GEOMETRY
 geometry.is_periodic = 0 0 0

Exec/science/Detonation/sdc_tests/make_plots.py

@@ -133,9 +133,10 @@ def get_velocity(self):
     def get_data(self, n=-1):
         """get the temperature and energy generation rate from the nth
         plotfile (starting to count at 0).""" 
-
-        return Profile(os.path.join(self.name, self.files[n]))
-
+        try:
+            return Profile(os.path.join(self.name, self.files[n]))
+        except IndexError:
+            return None
 
 if __name__ == "__main__":
 
@@ -153,9 +154,15 @@ def get_data(self, n=-1):
     runs.sort()
     print(runs)
 
-    # make a plot of the profiles for different CFLs and different resolutions
-    nzones = set([q.nzones for q in runs])
-    integrators = set([q.integrator for q in runs])
+    nzones = list(set([q.nzones for q in runs]))
+    nzones.sort()
+    integrators = list(set([q.integrator for q in runs]))
+    integrators.sort()
+    cfls = list(set([q.cfl for q in runs]))
+    cfls.sort()
+
+    # make a plot of the profiles for different CFLs with the same
+    # integrator and resolutions
 
     idx = 5
     for nz in nzones:
@@ -165,17 +172,19 @@ def get_data(self, n=-1):
             if len(dset) == 0:
                 continue
 
-            fig = plt.figure(1)
-            fig.clear()
+            fig, axes = plt.subplots(ncols=1, nrows=3, constrained_layout=True, squeeze=True)
 
             fig.set_size_inches(7.0, 9.0)
 
-            ax_T = fig.add_subplot(311)
-            ax_e = fig.add_subplot(312)
-            ax_he = fig.add_subplot(313)
+            ax_T = axes[0]
+            ax_e = axes[1]
+            ax_he = axes[2]
 
             for q in dset:
                 pf = q.get_data(idx)
+                if not pf:
+                    continue
+
                 print("working on {}, time = {}".format(q.name, pf.time))
 
                 ax_T.plot(pf.x, pf.T, label="CFL = {}".format(q.cfl))
@@ -198,5 +207,52 @@ def get_data(self, n=-1):
 
             fig.suptitle("{}, number of zones = {}".format(intg, nz))
 
-            fig.tight_layout()
             fig.savefig("det_cfl_compare_{}_nz{}.png".format(intg.replace(" ", "_"), nz))
+            plt.close()
+
+    # make a plot of the profiles for different integrators with the same CFLs and resolutions
+
+    idx = 5
+    for nz in nzones:
+        for c in cfls:
+            dset = [q for q in runs if q.cfl == c and q.nzones == nz]
+            dset.sort(key=lambda q: q.integrator)
+            if len(dset) == 0:
+                continue
+
+            fig, axes = plt.subplots(ncols=1, nrows=3, constrained_layout=True, squeeze=True)
+
+            fig.set_size_inches(7.0, 9.0)
+
+            ax_T = axes[0]
+            ax_e = axes[1]
+            ax_he = axes[2]
+
+            for q in dset:
+                pf = q.get_data(idx)
+                if not pf:
+                    continue
+                print("working on {}, time = {}".format(q.name, pf.time))
+
+                ax_T.plot(pf.x, pf.T, label="{}".format(q.integrator))
+                ax_e.plot(pf.x, pf.enuc)
+                ax_he.plot(pf.x, pf.rho_he4)
+
+            ax_T.legend(frameon=False)
+
+            ax_T.set_ylabel(r"$T$ (K)")
+            ax_e.set_ylabel(r"$H_\mathrm{nuc}$ (erg/g/s)")
+            ax_he.set_ylabel(r"$\rho X({}^4\mathrm{He})$ (g/cm${}^3$)")
+            ax_he.set_xlabel("x (cm)")
+
+            ax_T.set_xlim(5000, 15000)
+            ax_e.set_xlim(5000, 15000)
+            ax_he.set_xlim(5000, 15000)
+
+            ax_e.set_yscale("log")
+            ax_he.set_yscale("log")
+
+            fig.suptitle("number of zones = {}, CFL = {}".format(nz, c))
+
+            fig.savefig("det_integrator_compare_cfl{}_nz{}.png".format(c, nz))
+            plt.close()

Exec/science/Detonation/sdc_tests/probin.sdc

@@ -36,4 +36,6 @@
 !   rtol_enuc = 1.e-8
    call_eos_in_rhs = T
    do_constant_volume_burn = T
+
+   small_x = 1.e-10
 /

Exec/science/Detonation/sdc_tests/setup_runs.sh

@@ -26,7 +26,7 @@ stop_time=3.e-6"
 RUNPARAMS="
 castro.time_integration_method=2
 castro.sdc_order=4
-castro.sdc_quadrature=1
+castro.sdc_quadrature=0
 castro.limit_fourth_order=1
 castro.use_reconstructed_gamma1=1
 castro.sdc_solve_for_rhoe=1
@@ -59,14 +59,54 @@ do
 done
 
 
+#=============================================================================
+# Lobatto SDC-2
+#=============================================================================
+
+RUNPARAMS="
+castro.time_integration_method=2
+castro.sdc_order=2
+castro.sdc_quadrature=0
+castro.limit_fourth_order=1
+castro.use_reconstructed_gamma1=1
+castro.sdc_solve_for_rhoe=1
+castro.sdc_solver_tol_dens=1.e-8
+castro.sdc_solver_tol_spec=1.e-8
+castro.sdc_solver_tol_ener=1.e-8
+castro.sdc_solver=2"
+
+for c in ${CFL}
+do
+
+    for nz in ${NZONES}
+    do
+        rdir=det_z${nz}_c${c}_lobatto_sdc2
+        if [ ! -d ${rdir} ]; then
+            mkdir ${rdir}
+        fi
+
+        cd ${rdir}
+        for nf in ${NEEDED_FILES}
+        do
+            if [ ! -f ${nf} ]; then
+                cp ../${nf} .
+            fi
+        done
+
+        nohup ${CASTRO_EXEC} inputs.1d.sdc amr.plot_file=${rdir}_plt ${GLOBAL_RUNPARAMS} ${RUNPARAMS} castro.cfl=${c} amr.n_cell=${nz} >& out &
+        cd ..
+    done
+done
+
+
 #=============================================================================
 # Radau SDC-4
 #=============================================================================
 
 RUNPARAMS="
 castro.time_integration_method=2
 castro.sdc_order=4
-castro.sdc_quadrature=2
+castro.sdc_quadrature=1
 castro.limit_fourth_order=1
 castro.use_reconstructed_gamma1=1
 castro.sdc_solve_for_rhoe=1
@@ -99,6 +139,46 @@ do
 done
 
 
+#=============================================================================
+# Radau SDC-2
+#=============================================================================
+
+RUNPARAMS="
+castro.time_integration_method=2
+castro.sdc_order=2
+castro.sdc_quadrature=1
+castro.limit_fourth_order=1
+castro.use_reconstructed_gamma1=1
+castro.sdc_solve_for_rhoe=1
+castro.sdc_solver_tol_dens=1.e-8
+castro.sdc_solver_tol_spec=1.e-8
+castro.sdc_solver_tol_ener=1.e-8
+castro.sdc_solver=2"
+
+for c in ${CFL}
+do
+
+    for nz in ${NZONES}
+    do
+        rdir=det_z${nz}_c${c}_radau_sdc2
+        if [ ! -d ${rdir} ]; then
+            mkdir ${rdir}
+        fi
+
+        cd ${rdir}
+        for nf in ${NEEDED_FILES}
+        do
+            if [ ! -f ${nf} ]; then
+                cp ../${nf} .
+            fi
+        done
+
+        nohup ${CASTRO_EXEC} inputs.1d.sdc amr.plot_file=${rdir}_plt ${GLOBAL_RUNPARAMS} ${RUNPARAMS} castro.cfl=${c} amr.n_cell=${nz} >& out &
+        cd ..
+    done
+done
+
+
 
 #=============================================================================
 # Strang CTU

Source/driver/Castro_advance.cpp

@@ -540,6 +540,7 @@ Castro::initialize_advance(Real time, Real dt, int amr_iteration, int amr_ncycle
 
       MultiFab& S_old = get_old_data(State_Type);
       k_new.resize(SDC_NODES);
+
       k_new[0].reset(new MultiFab(S_old, amrex::make_alias, 0, NUM_STATE));
       for (int n = 1; n < SDC_NODES; ++n) {
 	k_new[n].reset(new MultiFab(grids, dmap, NUM_STATE, 0));

Source/driver/Castro_advance_sdc.cpp

@@ -77,10 +77,10 @@ Castro::do_advance_sdc (Real time,
 
 
     // the next chunk of code constructs the advective term for the
-    // current node, m.  First we get the sources, then convert to
-    // primitive, then get the source term from the MOL driver.  Note,
-    // for m = 0, we only have to do all of this the first iteration,
-    // since that state never changes
+    // current node, m.  First we get the sources, then full hydro
+    // source term from the MOL driver.  Note, for m = 0, we only have
+    // to do all of this the first iteration, since that state never
+    // changes
     if (!(sdc_iteration > 0 && m == 0) &&
         !(sdc_iteration == sdc_order+sdc_extra-1 && m == SDC_NODES-1)) {
 
@@ -205,6 +205,9 @@ Castro::do_advance_sdc (Real time,
     // this if we are on the final node (since there is nothing to
     // update to
     if (m < SDC_NODES-1) {
+
+      amrex::Print() << "... doing the SDC update, iteration = " << sdc_iteration << " from node " << m << " to " << m+1 << std::endl;
+
       do_sdc_update(m, m+1, dt); //(dt_sdc[m+1] - dt_sdc[m])*dt);
 
       // we now have a new value of k_new[m+1], do a clean_state on it