diff --git a/Exec/science/Detonation/inputs-det-x.sdc2 b/Exec/science/Detonation/inputs-det-x.sdc2 new file mode 100644 index 0000000000..286674d7c4 --- /dev/null +++ b/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 diff --git a/Exec/science/Detonation/sdc_tests/inputs.1d.sdc b/Exec/science/Detonation/sdc_tests/inputs.1d.sdc index cef4192c7e..69f2ab321e 100644 --- a/Exec/science/Detonation/sdc_tests/inputs.1d.sdc +++ b/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 diff --git a/Exec/science/Detonation/sdc_tests/make_plots.py b/Exec/science/Detonation/sdc_tests/make_plots.py index 4482f499ab..cb841ac204 100755 --- a/Exec/science/Detonation/sdc_tests/make_plots.py +++ b/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() diff --git a/Exec/science/Detonation/sdc_tests/probin.sdc b/Exec/science/Detonation/sdc_tests/probin.sdc index 2a3dad2295..af3c14fbdb 100644 --- a/Exec/science/Detonation/sdc_tests/probin.sdc +++ b/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 / diff --git a/Exec/science/Detonation/sdc_tests/setup_runs.sh b/Exec/science/Detonation/sdc_tests/setup_runs.sh index d8b01a8720..a101579f95 100755 --- a/Exec/science/Detonation/sdc_tests/setup_runs.sh +++ b/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,6 +59,46 @@ 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 #============================================================================= @@ -66,7 +106,7 @@ done 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 diff --git a/Source/driver/Castro_advance.cpp b/Source/driver/Castro_advance.cpp index b8ce22aff4..4f23c9eefd 100644 --- a/Source/driver/Castro_advance.cpp +++ b/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)); diff --git a/Source/driver/Castro_advance_sdc.cpp b/Source/driver/Castro_advance_sdc.cpp index a5951a22dc..8f5542eec1 100644 --- a/Source/driver/Castro_advance_sdc.cpp +++ b/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 diff --git a/Source/driver/Castro_sdc.cpp b/Source/driver/Castro_sdc.cpp index dd4b8c8174..c3a7e3d843 100644 --- a/Source/driver/Castro_sdc.cpp +++ b/Source/driver/Castro_sdc.cpp @@ -28,6 +28,9 @@ Castro::do_sdc_update(int m_start, int m_end, Real dt) { const int* domain_lo = geom.Domain().loVect(); const int* domain_hi = geom.Domain().hiVect(); + // the timestep from m to m+1 + Real dt_m = (dt_sdc[m_end] - dt_sdc[m_start]) * dt; + #ifdef REACTIONS // SDC_Source_Type is only defined for 4th order MultiFab tmp; @@ -42,16 +45,34 @@ Castro::do_sdc_update(int m_start, int m_end, Real dt) { const Box& bx = mfi.tilebox(); - ca_sdc_compute_C4(BL_TO_FORTRAN_BOX(bx), - BL_TO_FORTRAN_3D((*A_new[m_start])[mfi]), - BL_TO_FORTRAN_3D((*A_old[0])[mfi]), - BL_TO_FORTRAN_3D((*A_old[1])[mfi]), - BL_TO_FORTRAN_3D((*A_old[2])[mfi]), - BL_TO_FORTRAN_3D((*R_old[0])[mfi]), - BL_TO_FORTRAN_3D((*R_old[1])[mfi]), - BL_TO_FORTRAN_3D((*R_old[2])[mfi]), - BL_TO_FORTRAN_3D(C_source[mfi]), - &m_start); + if (sdc_quadrature == 0) { + + ca_sdc_compute_C4_lobatto(BL_TO_FORTRAN_BOX(bx), + &dt_m, &dt, + BL_TO_FORTRAN_3D((*A_new[m_start])[mfi]), + BL_TO_FORTRAN_3D((*A_old[0])[mfi]), + BL_TO_FORTRAN_3D((*A_old[1])[mfi]), + BL_TO_FORTRAN_3D((*A_old[2])[mfi]), + BL_TO_FORTRAN_3D((*R_old[0])[mfi]), + BL_TO_FORTRAN_3D((*R_old[1])[mfi]), + BL_TO_FORTRAN_3D((*R_old[2])[mfi]), + BL_TO_FORTRAN_3D(C_source[mfi]), + &m_start); + } else { + ca_sdc_compute_C4_radau(BL_TO_FORTRAN_BOX(bx), + &dt_m, &dt, + BL_TO_FORTRAN_3D((*A_new[m_start])[mfi]), + BL_TO_FORTRAN_3D((*A_old[0])[mfi]), + BL_TO_FORTRAN_3D((*A_old[1])[mfi]), + BL_TO_FORTRAN_3D((*A_old[2])[mfi]), + BL_TO_FORTRAN_3D((*A_old[3])[mfi]), + BL_TO_FORTRAN_3D((*R_old[0])[mfi]), + BL_TO_FORTRAN_3D((*R_old[1])[mfi]), + BL_TO_FORTRAN_3D((*R_old[2])[mfi]), + BL_TO_FORTRAN_3D((*R_old[3])[mfi]), + BL_TO_FORTRAN_3D(C_source[mfi]), + &m_start); + } } // need to construct the time for this stage -- but it is not really @@ -66,8 +87,6 @@ Castro::do_sdc_update(int m_start, int m_end, Real dt) { // staging place so we can do a FillPatch MultiFab& S_new = get_new_data(State_Type); - Real dt_m = dt/2.0; - for (MFIter mfi(S_new); mfi.isValid(); ++mfi) { const Box& bx = mfi.tilebox(); @@ -97,6 +116,8 @@ Castro::do_sdc_update(int m_start, int m_end, Real dt) { FArrayBox U_new_center; FArrayBox R_new; + FArrayBox C2; + for (MFIter mfi(*k_new[0]); mfi.isValid(); ++mfi) { const Box& bx = mfi.tilebox(); @@ -109,21 +130,48 @@ Castro::do_sdc_update(int m_start, int m_end, Real dt) { // second order SDC reaction update -- we don't care about // the difference between cell-centers and averages - ca_sdc_update_o2(BL_TO_FORTRAN_BOX(bx), &dt, + // first compute the source term, C -- this differs depending + // on whether we are Lobatto or Radau + C2.resize(bx, NUM_STATE); + + if (sdc_quadrature == 0) { + + ca_sdc_compute_C2_lobatto(BL_TO_FORTRAN_BOX(bx), + &dt_m, &dt, + BL_TO_FORTRAN_3D((*A_new[m_start])[mfi]), + BL_TO_FORTRAN_3D((*A_old[0])[mfi]), + BL_TO_FORTRAN_3D((*A_old[1])[mfi]), + BL_TO_FORTRAN_3D((*R_old[0])[mfi]), + BL_TO_FORTRAN_3D((*R_old[1])[mfi]), + BL_TO_FORTRAN_3D(C2), + &m_start); + + } else { + + ca_sdc_compute_C2_radau(BL_TO_FORTRAN_BOX(bx), + &dt_m, &dt, + BL_TO_FORTRAN_3D((*A_new[m_start])[mfi]), + BL_TO_FORTRAN_3D((*A_old[0])[mfi]), + BL_TO_FORTRAN_3D((*A_old[1])[mfi]), + BL_TO_FORTRAN_3D((*A_old[2])[mfi]), + BL_TO_FORTRAN_3D((*R_old[0])[mfi]), + BL_TO_FORTRAN_3D((*R_old[1])[mfi]), + BL_TO_FORTRAN_3D((*R_old[2])[mfi]), + BL_TO_FORTRAN_3D(C2), + &m_start); + + } + + ca_sdc_update_o2(BL_TO_FORTRAN_BOX(bx), &dt_m, BL_TO_FORTRAN_3D((*k_new[m_start])[mfi]), BL_TO_FORTRAN_3D((*k_new[m_end])[mfi]), BL_TO_FORTRAN_3D((*A_new[m_start])[mfi]), - BL_TO_FORTRAN_3D((*A_old[0])[mfi]), - BL_TO_FORTRAN_3D((*A_old[1])[mfi]), - BL_TO_FORTRAN_3D((*R_old[0])[mfi]), - BL_TO_FORTRAN_3D((*R_old[1])[mfi]), + BL_TO_FORTRAN_3D((*R_old[m_start])[mfi]), + BL_TO_FORTRAN_3D(C2), &sdc_iteration, &m_start); } else { - // the timestep from m to m+1 - Real dt_m = dt/2.0; - // fourth order SDC reaction update -- we need to respect the // difference between cell-centers and averages @@ -182,22 +230,50 @@ Castro::do_sdc_update(int m_start, int m_end, Real dt) { #else // pure advection if (sdc_order == 2) { - ca_sdc_update_advection_o2(BL_TO_FORTRAN_BOX(bx), &dt, - BL_TO_FORTRAN_3D((*k_new[m_start])[mfi]), - BL_TO_FORTRAN_3D((*k_new[m_end])[mfi]), - BL_TO_FORTRAN_3D((*A_new[m_start])[mfi]), - BL_TO_FORTRAN_3D((*A_old[0])[mfi]), - BL_TO_FORTRAN_3D((*A_old[1])[mfi]), - &m_start); + + if (sdc_quadrature == 0) { + ca_sdc_update_advection_o2_lobatto(BL_TO_FORTRAN_BOX(bx), &dt_m, &dt, + BL_TO_FORTRAN_3D((*k_new[m_start])[mfi]), + BL_TO_FORTRAN_3D((*k_new[m_end])[mfi]), + BL_TO_FORTRAN_3D((*A_new[m_start])[mfi]), + BL_TO_FORTRAN_3D((*A_old[0])[mfi]), + BL_TO_FORTRAN_3D((*A_old[1])[mfi]), + &m_start); + + } else { + ca_sdc_update_advection_o2_radau(BL_TO_FORTRAN_BOX(bx), &dt_m, &dt, + BL_TO_FORTRAN_3D((*k_new[m_start])[mfi]), + BL_TO_FORTRAN_3D((*k_new[m_end])[mfi]), + BL_TO_FORTRAN_3D((*A_new[m_start])[mfi]), + BL_TO_FORTRAN_3D((*A_old[0])[mfi]), + BL_TO_FORTRAN_3D((*A_old[1])[mfi]), + BL_TO_FORTRAN_3D((*A_old[2])[mfi]), + &m_start); + + } + } else { - ca_sdc_update_advection_o4(BL_TO_FORTRAN_BOX(bx), &dt, - BL_TO_FORTRAN_3D((*k_new[m_start])[mfi]), - BL_TO_FORTRAN_3D((*k_new[m_end])[mfi]), - BL_TO_FORTRAN_3D((*A_new[m_start])[mfi]), - BL_TO_FORTRAN_3D((*A_old[0])[mfi]), - BL_TO_FORTRAN_3D((*A_old[1])[mfi]), - BL_TO_FORTRAN_3D((*A_old[2])[mfi]), - &m_start); + if (sdc_quadrature == 0) { + ca_sdc_update_advection_o4_lobatto(BL_TO_FORTRAN_BOX(bx), &dt_m, &dt, + BL_TO_FORTRAN_3D((*k_new[m_start])[mfi]), + BL_TO_FORTRAN_3D((*k_new[m_end])[mfi]), + BL_TO_FORTRAN_3D((*A_new[m_start])[mfi]), + BL_TO_FORTRAN_3D((*A_old[0])[mfi]), + BL_TO_FORTRAN_3D((*A_old[1])[mfi]), + BL_TO_FORTRAN_3D((*A_old[2])[mfi]), + &m_start); + } else { + ca_sdc_update_advection_o4_radau(BL_TO_FORTRAN_BOX(bx), &dt_m, &dt, + BL_TO_FORTRAN_3D((*k_new[m_start])[mfi]), + BL_TO_FORTRAN_3D((*k_new[m_end])[mfi]), + BL_TO_FORTRAN_3D((*A_new[m_start])[mfi]), + BL_TO_FORTRAN_3D((*A_old[0])[mfi]), + BL_TO_FORTRAN_3D((*A_old[1])[mfi]), + BL_TO_FORTRAN_3D((*A_old[2])[mfi]), + BL_TO_FORTRAN_3D((*A_old[3])[mfi]), + &m_start); + } + } #endif diff --git a/Source/driver/Castro_setup.cpp b/Source/driver/Castro_setup.cpp index b7843e8893..5850405680 100644 --- a/Source/driver/Castro_setup.cpp +++ b/Source/driver/Castro_setup.cpp @@ -1090,28 +1090,63 @@ Castro::variableSetUp () #endif - if (sdc_order == 2) { - // Gauss-Lobatto (trapezoid) - SDC_NODES = 2; - dt_sdc.resize(SDC_NODES); - dt_sdc = {0.0, 1.0}; + if (sdc_quadrature == 0) { + // Gauss-Lobatto - node_weights.resize(SDC_NODES); - node_weights = {0.5, 0.5}; + if (sdc_order == 2) { + // trapezoid + SDC_NODES = 2; + + dt_sdc.resize(SDC_NODES); + dt_sdc = {0.0, 1.0}; + + node_weights.resize(SDC_NODES); + node_weights = {0.5, 0.5}; + + } else if (sdc_order == 4) { + // Simpsons + SDC_NODES = 3; + + dt_sdc.resize(SDC_NODES); + dt_sdc = {0.0, 0.5, 1.0}; + + node_weights.resize(SDC_NODES); + node_weights = {1.0/6.0, 4.0/6.0, 1.0/6.0}; + + } else { + amrex::Error("invalid value of sdc_order"); + } - } else if (sdc_order == 4) { - // Gauss-Lobatto (Simpsons) - SDC_NODES = 3; + } else if (sdc_quadrature == 1) { + // Radau - dt_sdc.resize(SDC_NODES); - dt_sdc = {0.0, 0.5, 1.0}; + if (sdc_order == 2) { + SDC_NODES = 3; + + dt_sdc.resize(SDC_NODES); + dt_sdc = {0.0, 1.0/3.0, 1.0}; + + node_weights.resize(SDC_NODES); + node_weights = {0.0, 3.0/4.0, 1.0/4.0}; + + } else if (sdc_order == 4) { + SDC_NODES = 4; + + dt_sdc.resize(SDC_NODES); + dt_sdc = {0.0, (4.0 - std::sqrt(6.0))/10.0, (4.0 + std::sqrt(6.0))/10.0, 1.0}; + + node_weights.resize(SDC_NODES); + node_weights = {0.0, (16.0 - std::sqrt(6.0))/36.0, (16.0 + std::sqrt(6.0))/36.0, 1.0/9.0}; + + } else { + amrex::Error("invalid value of sdc_order"); + } node_weights.resize(SDC_NODES); node_weights = {1.0/6.0, 4.0/6.0, 1.0/6.0}; } else { - amrex::Error("invalid value of sdc_order"); + amrex::Error("invalid value of sdc_quadrature"); } - } diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters index 313e6e42eb..1083dec51b 100644 --- a/Source/driver/_cpp_parameters +++ b/Source/driver/_cpp_parameters @@ -239,6 +239,9 @@ hse_reflect_vels int 0 y # valid options are 2 and 4 sdc_order int 2 y +# which quadrature type to use with SDC? 0 = Gauss-Lobatto, 1 = Radau +sdc_quadrature int 0 y + # number of extra SDC iterations to take beyond the order sdc_extra int 0 y diff --git a/Source/driver/meth_params.F90 b/Source/driver/meth_params.F90 index f2bdda0203..d82ff0535e 100644 --- a/Source/driver/meth_params.F90 +++ b/Source/driver/meth_params.F90 @@ -164,6 +164,7 @@ module meth_params_module integer, allocatable, save :: hse_interp_temp integer, allocatable, save :: hse_reflect_vels integer, allocatable, save :: sdc_order + integer, allocatable, save :: sdc_quadrature integer, allocatable, save :: sdc_extra integer, allocatable, save :: sdc_solver real(rt), allocatable, save :: sdc_solver_tol_dens @@ -258,6 +259,7 @@ module meth_params_module attributes(managed) :: hse_interp_temp attributes(managed) :: hse_reflect_vels attributes(managed) :: sdc_order +attributes(managed) :: sdc_quadrature attributes(managed) :: sdc_extra attributes(managed) :: sdc_solver attributes(managed) :: sdc_solver_tol_dens @@ -379,6 +381,7 @@ module meth_params_module !$acc create(hse_interp_temp) & !$acc create(hse_reflect_vels) & !$acc create(sdc_order) & + !$acc create(sdc_quadrature) & !$acc create(sdc_extra) & !$acc create(sdc_solver) & !$acc create(sdc_solver_tol_dens) & @@ -619,6 +622,8 @@ subroutine ca_set_castro_method_params() bind(C, name="ca_set_castro_method_para hse_reflect_vels = 0; allocate(sdc_order) sdc_order = 2; + allocate(sdc_quadrature) + sdc_quadrature = 0; allocate(sdc_extra) sdc_extra = 0; allocate(sdc_solver) @@ -742,6 +747,7 @@ subroutine ca_set_castro_method_params() bind(C, name="ca_set_castro_method_para call pp%query("hse_interp_temp", hse_interp_temp) call pp%query("hse_reflect_vels", hse_reflect_vels) call pp%query("sdc_order", sdc_order) + call pp%query("sdc_quadrature", sdc_quadrature) call pp%query("sdc_extra", sdc_extra) call pp%query("sdc_solver", sdc_solver) call pp%query("sdc_solver_tol_dens", sdc_solver_tol_dens) @@ -800,21 +806,22 @@ subroutine ca_set_castro_method_params() bind(C, name="ca_set_castro_method_para !$acc device(use_pslope, limit_fluxes_on_small_dens, density_reset_method) & !$acc device(allow_small_energy, do_sponge, sponge_implicit) & !$acc device(first_order_hydro, hse_zero_vels, hse_interp_temp) & - !$acc device(hse_reflect_vels, sdc_order, sdc_extra) & - !$acc device(sdc_solver, sdc_solver_tol_dens, sdc_solver_tol_spec) & - !$acc device(sdc_solver_tol_ener, sdc_solver_atol, sdc_solver_relax_factor) & - !$acc device(sdc_solve_for_rhoe, sdc_use_analytic_jac, cfl) & - !$acc device(dtnuc_e, dtnuc_X, dtnuc_X_threshold) & - !$acc device(do_react, react_T_min, react_T_max) & - !$acc device(react_rho_min, react_rho_max, disable_shock_burning) & - !$acc device(T_guess, diffuse_temp, diffuse_cutoff_density) & - !$acc device(diffuse_cutoff_density_hi, diffuse_cond_scale_fac, do_grav) & - !$acc device(grav_source_type, do_rotation, rot_period) & - !$acc device(rot_period_dot, rotation_include_centrifugal, rotation_include_coriolis) & - !$acc device(rotation_include_domegadt, state_in_rotating_frame, rot_source_type) & - !$acc device(implicit_rotation_update, rot_axis, use_point_mass) & - !$acc device(point_mass, point_mass_fix_solution, do_acc) & - !$acc device(grown_factor, track_grid_losses, const_grav, get_g_from_phi) + !$acc device(hse_reflect_vels, sdc_order, sdc_quadrature) & + !$acc device(sdc_extra, sdc_solver, sdc_solver_tol_dens) & + !$acc device(sdc_solver_tol_spec, sdc_solver_tol_ener, sdc_solver_atol) & + !$acc device(sdc_solver_relax_factor, sdc_solve_for_rhoe, sdc_use_analytic_jac) & + !$acc device(cfl, dtnuc_e, dtnuc_X) & + !$acc device(dtnuc_X_threshold, do_react, react_T_min) & + !$acc device(react_T_max, react_rho_min, react_rho_max) & + !$acc device(disable_shock_burning, T_guess, diffuse_temp) & + !$acc device(diffuse_cutoff_density, diffuse_cutoff_density_hi, diffuse_cond_scale_fac) & + !$acc device(do_grav, grav_source_type, do_rotation) & + !$acc device(rot_period, rot_period_dot, rotation_include_centrifugal) & + !$acc device(rotation_include_coriolis, rotation_include_domegadt, state_in_rotating_frame) & + !$acc device(rot_source_type, implicit_rotation_update, rot_axis) & + !$acc device(use_point_mass, point_mass, point_mass_fix_solution) & + !$acc device(do_acc, grown_factor, track_grid_losses) & + !$acc device(const_grav, get_g_from_phi) #ifdef GRAVITY @@ -1055,6 +1062,9 @@ subroutine ca_finalize_meth_params() bind(C, name="ca_finalize_meth_params") if (allocated(sdc_order)) then deallocate(sdc_order) end if + if (allocated(sdc_quadrature)) then + deallocate(sdc_quadrature) + end if if (allocated(sdc_extra)) then deallocate(sdc_extra) end if diff --git a/Source/driver/param_includes/castro_defaults.H b/Source/driver/param_includes/castro_defaults.H index d79fe94e87..8bb0e64a20 100644 --- a/Source/driver/param_includes/castro_defaults.H +++ b/Source/driver/param_includes/castro_defaults.H @@ -70,6 +70,7 @@ int Castro::hse_zero_vels = 0; int Castro::hse_interp_temp = 0; int Castro::hse_reflect_vels = 0; int Castro::sdc_order = 2; +int Castro::sdc_quadrature = 0; int Castro::sdc_extra = 0; int Castro::sdc_solver = 1; amrex::Real Castro::sdc_solver_tol_dens = 1.e-6; diff --git a/Source/driver/param_includes/castro_job_info_tests.H b/Source/driver/param_includes/castro_job_info_tests.H index de2b2f8889..75646752ff 100644 --- a/Source/driver/param_includes/castro_job_info_tests.H +++ b/Source/driver/param_includes/castro_job_info_tests.H @@ -65,6 +65,7 @@ jobInfoFile << (Castro::hse_zero_vels == 0 ? " " : "[*] ") << "castro.hse_zer jobInfoFile << (Castro::hse_interp_temp == 0 ? " " : "[*] ") << "castro.hse_interp_temp = " << Castro::hse_interp_temp << std::endl; jobInfoFile << (Castro::hse_reflect_vels == 0 ? " " : "[*] ") << "castro.hse_reflect_vels = " << Castro::hse_reflect_vels << std::endl; jobInfoFile << (Castro::sdc_order == 2 ? " " : "[*] ") << "castro.sdc_order = " << Castro::sdc_order << std::endl; +jobInfoFile << (Castro::sdc_quadrature == 0 ? " " : "[*] ") << "castro.sdc_quadrature = " << Castro::sdc_quadrature << std::endl; jobInfoFile << (Castro::sdc_extra == 0 ? " " : "[*] ") << "castro.sdc_extra = " << Castro::sdc_extra << std::endl; jobInfoFile << (Castro::sdc_solver == 1 ? " " : "[*] ") << "castro.sdc_solver = " << Castro::sdc_solver << std::endl; jobInfoFile << (Castro::sdc_solver_tol_dens == 1.e-6 ? " " : "[*] ") << "castro.sdc_solver_tol_dens = " << Castro::sdc_solver_tol_dens << std::endl; diff --git a/Source/driver/param_includes/castro_params.H b/Source/driver/param_includes/castro_params.H index db5bf1610a..9fb74b9c73 100644 --- a/Source/driver/param_includes/castro_params.H +++ b/Source/driver/param_includes/castro_params.H @@ -70,6 +70,7 @@ static int hse_zero_vels; static int hse_interp_temp; static int hse_reflect_vels; static int sdc_order; +static int sdc_quadrature; static int sdc_extra; static int sdc_solver; static amrex::Real sdc_solver_tol_dens; diff --git a/Source/driver/param_includes/castro_queries.H b/Source/driver/param_includes/castro_queries.H index ecbe137974..3d7080dca1 100644 --- a/Source/driver/param_includes/castro_queries.H +++ b/Source/driver/param_includes/castro_queries.H @@ -70,6 +70,7 @@ pp.query("hse_zero_vels", hse_zero_vels); pp.query("hse_interp_temp", hse_interp_temp); pp.query("hse_reflect_vels", hse_reflect_vels); pp.query("sdc_order", sdc_order); +pp.query("sdc_quadrature", sdc_quadrature); pp.query("sdc_extra", sdc_extra); pp.query("sdc_solver", sdc_solver); pp.query("sdc_solver_tol_dens", sdc_solver_tol_dens); diff --git a/Source/driver/sdc_util.F90 b/Source/driver/sdc_util.F90 index 2dec54fbe4..21489e4943 100644 --- a/Source/driver/sdc_util.F90 +++ b/Source/driver/sdc_util.F90 @@ -55,7 +55,7 @@ subroutine sdc_solve(dt_m, U_old, U_new, C, sdc_iteration) ! reaction update. It either directly calls the Newton method or first ! tries VODE and then does the Newton update. - use meth_params_module, only : NVAR, sdc_solver + use meth_params_module, only : NVAR, sdc_solver, URHO, UTEMP, UEINT, UFS use amrex_constants_module, only : ZERO, HALF, ONE use burn_type_module, only : burn_t use react_util_module @@ -72,6 +72,10 @@ subroutine sdc_solve(dt_m, U_old, U_new, C, sdc_iteration) integer :: ierr + ! for debugging + real(rt) :: U_orig(NVAR) + + U_orig(:) = U_old(:) if (sdc_solver == NEWTON_SOLVE) then ! we are going to assume we already have a good guess for the @@ -81,6 +85,13 @@ subroutine sdc_solve(dt_m, U_old, U_new, C, sdc_iteration) ! failing? if (ierr /= NEWTON_SUCCESS) then + print *, "Newton convergence failure" + print *, " input state:" + print *, " density: ", U_orig(URHO) + print *, " temperature : ", U_orig(UTEMP) + print *, " (rho e): ", U_orig(UEINT) + print *, " mass fractions: ", U_orig(UFS:UFS-1+nspec)/U_orig(URHO) + print *, " " call castro_error("Newton subcycling failed in sdc_solve") end if @@ -117,6 +128,7 @@ subroutine sdc_newton_subdivide(dt_m, U_old, U_new, C, sdc_iteration, ierr) use amrex_constants_module, only : ZERO, HALF, ONE use network, only : nspec, nspec_evolve use rpar_sdc_module + use extern_probin_module, only : small_x implicit none @@ -142,16 +154,18 @@ subroutine sdc_newton_subdivide(dt_m, U_old, U_new, C, sdc_iteration, ierr) do while (nsub < MAX_NSUB .and. ierr /= NEWTON_SUCCESS) dt_sub = dt_m / nsub do isub = 1, nsub - call sdc_newton_solve(dt_sub, U_begin, U_new, C, sdc_iteration, ierr) - U_begin(:) = U_new(:) ! normalize species do n = 1, nspec - U_begin(UFS-1+n) = max(ZERO, U_begin(UFS-1+n)) + U_begin(UFS-1+n) = max(small_x, U_begin(UFS-1+n)) end do + sum_rhoX = sum(U_begin(UFS:UFS-1+nspec)) U_begin(UFS:UFS-1+nspec) = U_begin(UFS:UFS-1+nspec) * U_begin(URHO)/sum_rhoX + call sdc_newton_solve(dt_sub, U_begin, U_new, C, sdc_iteration, ierr) + U_begin(:) = U_new(:) + end do nsub = nsub * 2 end do @@ -176,6 +190,7 @@ subroutine sdc_newton_solve(dt_m, U_old, U_new, C, sdc_iteration, ierr) use react_util_module use network, only : nspec, nspec_evolve use rpar_sdc_module + use extern_probin_module, only : small_x implicit none @@ -215,6 +230,9 @@ subroutine sdc_newton_solve(dt_m, U_old, U_new, C, sdc_iteration, ierr) integer :: max_newton_iter + real(rt) :: xn(nspec) + integer :: k + ierr = NEWTON_SUCCESS ! the tolerance we are solving to may depend on the iteration @@ -307,6 +325,18 @@ subroutine sdc_newton_solve(dt_m, U_old, U_new, C, sdc_iteration, ierr) U_react(:) = U_react(:) + dU_react(:) + ! we still need to normalize here + xn(1:nspec_evolve) = U_react(1:nspec_evolve)/U_react(0) + xn(nspec_evolve+1:nspec) = rpar(irp_spec:irp_spec-1+(nspec-nspec_evolve))/U_react(0) + + do k = 1, nspec + xn(k) = max(small_x, xn(k)) + end do + xn(:) = xn(:)/sum(xn) + + U_react(1:nspec_evolve) = U_react(0) * xn(1:nspec_evolve) + rpar(irp_spec:irp_spec-1+(nspec-nspec_evolve)) = U_react(0) * xn(nspec_evolve+1:nspec) + eps_tot(0) = tol_dens * abs(U_react(0)) + sdc_solver_atol ! for species, atol is the mass fraction limit, so we multiply by density to get a partial density limit eps_tot(1:nspec_evolve) = tol_spec * abs(U_react(1:nspec_evolve)) + sdc_solver_atol * abs(U_react(0)) @@ -567,6 +597,7 @@ subroutine jac_ode(n, time, U, ml, mu, Jac, nrowpd, rpar, ipar) use eos_module, only : eos use rpar_sdc_module use react_util_module + use extern_probin_module, only : small_x implicit none @@ -641,70 +672,8 @@ subroutine jac_ode(n, time, U, ml, mu, Jac, nrowpd, rpar, ipar) end subroutine jac_ode - subroutine f_sdc(n, U, f, iflag, rpar) - ! this is used by the Newton solve to compute the Jacobian via differencing - - use rpar_sdc_module - use meth_params_module, only : nvar, URHO, UFS, UEDEN, UMX, UMZ, UEINT, UTEMP, sdc_solve_for_rhoe - use network, only : nspec, nspec_evolve - use burn_type_module - use react_util_module - - ! this computes the function we need to zero for the SDC update - implicit none - - integer,intent(in) :: n - real(rt), intent(in) :: U(0:n-1) - real(rt), intent(out) :: f(0:n-1) - integer, intent(inout) :: iflag !! leave this untouched - real(rt), intent(inout) :: rpar(0:n_rpar-1) - - real(rt) :: U_full(nvar), R_full(nvar) - real(rt) :: R_react(0:n-1), f_source(0:n-1) - type(burn_t) :: burn_state - - real(rt) :: dt_m - - ! we are not solving the momentum equations - ! create a full state -- we need this for some interfaces - U_full(URHO) = U(0) - U_full(UFS:UFS-1+nspec_evolve) = U(1:nspec_evolve) - if (sdc_solve_for_rhoe == 1) then - U_full(UEINT) = U(nspec_evolve+1) - U_full(UEDEN) = rpar(irp_evar) - else - U_full(UEDEN) = U(nspec_evolve+1) - U_full(UEINT) = rpar(irp_evar) - endif - - U_full(UMX:UMZ) = rpar(irp_mom:irp_mom+2) - U_full(UFS+nspec_evolve:UFS-1+nspec) = rpar(irp_spec:irp_spec-1+(nspec-nspec_evolve)) - - ! unpack rpar - dt_m = rpar(irp_dt) - f_source(:) = rpar(irp_f_source:irp_f_source-1+nspec_evolve+2) - ! initial guess for T - U_full(UTEMP) = rpar(irp_temp) - - call single_zone_react_source(U_full, R_full, 0,0,0, burn_state) - - ! update guess for next time - rpar(irp_temp) = U_full(UTEMP) - - R_react(0) = R_full(URHO) - R_react(1:nspec_evolve) = R_full(UFS:UFS-1+nspec_evolve) - if (sdc_solve_for_rhoe == 1) then - R_react(nspec_evolve+1) = R_full(UEINT) - else - R_react(nspec_evolve+1) = R_full(UEDEN) - endif - - f(:) = U(:) - dt_m * R_react(:) - f_source(:) - - end subroutine f_sdc - - subroutine f_sdc_jac(n, U, f, Jac, ldjac, iflag, rpar) + subroutine f_sdc_jac(neq, U, f, Jac, ldjac, iflag, rpar) ! this is used with the Newton solve and returns f and the Jacobian use rpar_sdc_module @@ -716,26 +685,28 @@ subroutine f_sdc_jac(n, U, f, Jac, ldjac, iflag, rpar) use eos_type_module, only : eos_t, eos_input_re use eos_module, only : eos use amrex_constants_module, only : ZERO, HALF, ONE + use extern_probin_module, only : small_x ! this computes the function we need to zero for the SDC update implicit none - integer,intent(in) :: n, ldjac - real(rt), intent(in) :: U(0:n-1) - real(rt), intent(out) :: f(0:n-1) - real(rt), intent(out) :: Jac(0:ldjac-1,0:n-1) + integer,intent(in) :: neq, ldjac + real(rt), intent(in) :: U(0:neq-1) + real(rt), intent(out) :: f(0:neq-1) + real(rt), intent(out) :: Jac(0:ldjac-1,0:neq-1) integer, intent(inout) :: iflag !! leave this untouched real(rt), intent(inout) :: rpar(0:n_rpar-1) real(rt) :: U_full(nvar), R_full(nvar) - real(rt) :: R_react(0:n-1), f_source(0:n-1) + real(rt) :: R_react(0:neq-1), f_source(0:neq-1) type(burn_t) :: burn_state type(eos_t) :: eos_state real(rt) :: dt_m real(rt) :: denom real(rt) :: dRdw(0:nspec_evolve+1, 0:nspec_evolve+1), dwdU(0:nspec_evolve+1, 0:nspec_evolve+1) - integer :: m + integer :: m, k + real(rt) :: sum_rhoX ! we are not solving the momentum equations ! create a full state -- we need this for some interfaces @@ -752,6 +723,14 @@ subroutine f_sdc_jac(n, U, f, Jac, ldjac, iflag, rpar) U_full(UMX:UMZ) = rpar(irp_mom:irp_mom+2) U_full(UFS+nspec_evolve:UFS-1+nspec) = rpar(irp_spec:irp_spec-1+(nspec-nspec_evolve)) + ! normalize the species + do k = 1, nspec + U_full(UFS-1+k) = max(small_x, U_full(UFS-1+k)) + end do + + sum_rhoX = sum(U_full(UFS:UFS-1+nspec)) + U_full(UFS:UFS-1+nspec) = U_full(UFS:UFS-1+nspec) * U_full(URHO)/sum_rhoX + ! unpack rpar dt_m = rpar(irp_dt) f_source(:) = rpar(irp_f_source:irp_f_source-1+nspec_evolve+2) @@ -828,14 +807,18 @@ subroutine f_sdc_jac(n, U, f, Jac, ldjac, iflag, rpar) end subroutine f_sdc_jac #endif - subroutine ca_sdc_update_advection_o2(lo, hi, dt_m, & - k_m, kmlo, kmhi, & - k_n, knlo, knhi, & - A_m, Amlo, Amhi, & - A_0_old, A0lo, A0hi, & - A_1_old, A1lo, A1hi, & - m_start) bind(C, name="ca_sdc_update_advection_o2") + subroutine ca_sdc_update_advection_o2_lobatto(lo, hi, dt_m, dt, & + k_m, kmlo, kmhi, & + k_n, knlo, knhi, & + A_m, Amlo, Amhi, & + A_0_old, A0lo, A0hi, & + A_1_old, A1lo, A1hi, & + m_start) bind(C, name="ca_sdc_update_advection_o2_lobatto") ! update k_m to k_n via advection -- this is a second-order accurate update + ! for the Gauss-Lobatto discretization of the time nodes + + ! here, dt_m is the update for this stage, from one time node to the next + ! dt is the update over the whole timestep, n to n+1 use meth_params_module, only : NVAR use amrex_constants_module, only : HALF @@ -843,7 +826,7 @@ subroutine ca_sdc_update_advection_o2(lo, hi, dt_m, & implicit none integer, intent(in) :: lo(3), hi(3) - real(rt), intent(in) :: dt_m + real(rt), intent(in) :: dt_m, dt integer, intent(in) :: kmlo(3), kmhi(3) integer, intent(in) :: knlo(3), knhi(3) integer, intent(in) :: Amlo(3), Amhi(3) @@ -861,35 +844,40 @@ subroutine ca_sdc_update_advection_o2(lo, hi, dt_m, & integer :: i, j, k + ! Gauss-Lobatto / trapezoid + do k = lo(3), hi(3) do j = lo(2), hi(2) do i = lo(1), hi(1) - k_n(i,j,k,:) = k_m(i,j,k,:) + HALF * dt_m * (A_0_old(i,j,k,:) + A_1_old(i,j,k,:)) - enddo - enddo - enddo + k_n(i,j,k,:) = k_m(i,j,k,:) + HALF * dt * (A_0_old(i,j,k,:) + A_1_old(i,j,k,:)) + end do + end do + end do - end subroutine ca_sdc_update_advection_o2 + end subroutine ca_sdc_update_advection_o2_lobatto - subroutine ca_sdc_update_advection_o4(lo, hi, dt, & - k_m, kmlo, kmhi, & - k_n, knlo, knhi, & - A_m, Amlo, Amhi, & - A_0_old, A0lo, A0hi, & - A_1_old, A1lo, A1hi, & - A_2_old, A2lo, A2hi, & - m_start) bind(C, name="ca_sdc_update_advection_o4") + subroutine ca_sdc_update_advection_o2_radau(lo, hi, dt_m, dt, & + k_m, kmlo, kmhi, & + k_n, knlo, knhi, & + A_m, Amlo, Amhi, & + A_0_old, A0lo, A0hi, & + A_1_old, A1lo, A1hi, & + A_2_old, A2lo, A2hi, & + m_start) bind(C, name="ca_sdc_update_advection_o2_radau") ! update k_m to k_n via advection -- this is a second-order accurate update - ! dt is the total timestep from n to n+1 + ! for the Radau discretization of the time nodes + + ! here, dt_m is the update for this stage, from one time node to the next + ! dt is the update over the whole timestep, n to n+1 use meth_params_module, only : NVAR - use amrex_constants_module, only : HALF, TWO, FIVE, EIGHT + use amrex_constants_module, only : HALF, FIVE implicit none integer, intent(in) :: lo(3), hi(3) - real(rt), intent(in) :: dt + real(rt), intent(in) :: dt_m, dt integer, intent(in) :: kmlo(3), kmhi(3) integer, intent(in) :: knlo(3), knhi(3) integer, intent(in) :: Amlo(3), Amhi(3) @@ -908,10 +896,82 @@ subroutine ca_sdc_update_advection_o4(lo, hi, dt, & real(rt), intent(in) :: A_2_old(A2lo(1):A2hi(1), A2lo(2):A2hi(2), A2lo(3):A2hi(3), NVAR) integer :: i, j, k - real(rt) :: dt_m + ! Radau + + if (m_start == 0) then + + do k = lo(3), hi(3) + do j = lo(2), hi(2) + do i = lo(1), hi(1) + + k_n(i,j,k,:) = k_m(i,j,k,:) + & + dt_m * (A_m(i,j,k,:) - A_0_old(i,j,k,:)) + & + dt/12.0_rt * (FIVE*A_1_old(i,j,k,:) - A_2_old(i,j,k,:)) + + end do + end do + end do + + else if (m_start == 1) then + + do k = lo(3), hi(3) + do j = lo(2), hi(2) + do i = lo(1), hi(1) + + k_n(i,j,k,:) = k_m(i,j,k,:) + & + dt_m * (A_m(i,j,k,:) - A_1_old(i,j,k,:)) + & + dt/3.0_rt * (A_1_old(i,j,k,:) + A_2_old(i,j,k,:)) + + end do + end do + end do + + end if + + end subroutine ca_sdc_update_advection_o2_radau + + + subroutine ca_sdc_update_advection_o4_lobatto(lo, hi, dt_m, dt, & + k_m, kmlo, kmhi, & + k_n, knlo, knhi, & + A_m, Amlo, Amhi, & + A_0_old, A0lo, A0hi, & + A_1_old, A1lo, A1hi, & + A_2_old, A2lo, A2hi, & + m_start) bind(C, name="ca_sdc_update_advection_o4_lobatto") + ! update k_m to k_n via advection -- this is a fourth order accurate update + + ! here, dt_m is the update for this stage, from one time node to the next + ! dt is the update over the whole timestep, n to n+1 + + use meth_params_module, only : NVAR + use amrex_constants_module, only : HALF, TWO, FIVE, EIGHT - dt_m = HALF * dt + implicit none + + integer, intent(in) :: lo(3), hi(3) + real(rt), intent(in) :: dt_m, dt + integer, intent(in) :: kmlo(3), kmhi(3) + integer, intent(in) :: knlo(3), knhi(3) + integer, intent(in) :: Amlo(3), Amhi(3) + integer, intent(in) :: A0lo(3), A0hi(3) + integer, intent(in) :: A1lo(3), A1hi(3) + integer, intent(in) :: A2lo(3), A2hi(3) + integer, intent(in) :: m_start + + + real(rt), intent(in) :: k_m(kmlo(1):kmhi(1), kmlo(2):kmhi(2), kmlo(3):kmhi(3), NVAR) + real(rt), intent(inout) :: k_n(knlo(1):knhi(1), knlo(2):knhi(2), knlo(3):knhi(3), NVAR) + + real(rt), intent(in) :: A_m(Amlo(1):Amhi(1), Amlo(2):Amhi(2), Amlo(3):Amhi(3), NVAR) + real(rt), intent(in) :: A_0_old(A0lo(1):A0hi(1), A0lo(2):A0hi(2), A0lo(3):A0hi(3), NVAR) + real(rt), intent(in) :: A_1_old(A1lo(1):A1hi(1), A1lo(2):A1hi(2), A1lo(3):A1hi(3), NVAR) + real(rt), intent(in) :: A_2_old(A2lo(1):A2hi(1), A2lo(2):A2hi(2), A2lo(3):A2hi(3), NVAR) + + integer :: i, j, k + + ! Gauss-Lobatto (Simpsons) if (m_start == 0) then @@ -938,23 +998,118 @@ subroutine ca_sdc_update_advection_o4(lo, hi, dt, & enddo else - call castro_error("error in ca_sdc_update_advection_o4 -- should not be here") + call castro_error("error in ca_sdc_update_advection_o4_lobatto -- should not be here") endif - end subroutine ca_sdc_update_advection_o4 + end subroutine ca_sdc_update_advection_o4_lobatto + + + subroutine ca_sdc_update_advection_o4_radau(lo, hi, dt_m, dt, & + k_m, kmlo, kmhi, & + k_n, knlo, knhi, & + A_m, Amlo, Amhi, & + A_0_old, A0lo, A0hi, & + A_1_old, A1lo, A1hi, & + A_2_old, A2lo, A2hi, & + A_3_old, A3lo, A3hi, & + m_start) bind(C, name="ca_sdc_update_advection_o4_radau") + ! update k_m to k_n via advection -- this is a fourth-order accurate update + + ! here, dt_m is the update for this stage, from one time node to the next + ! dt is the update over the whole timestep, n to n+1 + + use meth_params_module, only : NVAR + use amrex_constants_module, only : HALF, TWO, FIVE, EIGHT + + implicit none + + integer, intent(in) :: lo(3), hi(3) + real(rt), intent(in) :: dt_m, dt + integer, intent(in) :: kmlo(3), kmhi(3) + integer, intent(in) :: knlo(3), knhi(3) + integer, intent(in) :: Amlo(3), Amhi(3) + integer, intent(in) :: A0lo(3), A0hi(3) + integer, intent(in) :: A1lo(3), A1hi(3) + integer, intent(in) :: A2lo(3), A2hi(3) + integer, intent(in) :: A3lo(3), A3hi(3) + integer, intent(in) :: m_start + + + real(rt), intent(in) :: k_m(kmlo(1):kmhi(1), kmlo(2):kmhi(2), kmlo(3):kmhi(3), NVAR) + real(rt), intent(inout) :: k_n(knlo(1):knhi(1), knlo(2):knhi(2), knlo(3):knhi(3), NVAR) + + real(rt), intent(in) :: A_m(Amlo(1):Amhi(1), Amlo(2):Amhi(2), Amlo(3):Amhi(3), NVAR) + real(rt), intent(in) :: A_0_old(A0lo(1):A0hi(1), A0lo(2):A0hi(2), A0lo(3):A0hi(3), NVAR) + real(rt), intent(in) :: A_1_old(A1lo(1):A1hi(1), A1lo(2):A1hi(2), A1lo(3):A1hi(3), NVAR) + real(rt), intent(in) :: A_2_old(A2lo(1):A2hi(1), A2lo(2):A2hi(2), A2lo(3):A2hi(3), NVAR) + real(rt), intent(in) :: A_3_old(A3lo(1):A3hi(1), A3lo(2):A3hi(2), A3lo(3):A3hi(3), NVAR) + + integer :: i, j, k + + ! Gauss-Lobatto (Simpsons) + + if (m_start == 0) then + + do k = lo(3), hi(3) + do j = lo(2), hi(2) + do i = lo(1), hi(1) + k_n(i,j,k,:) = k_m(i,j,k,:) + & + dt_m * (A_m(i,j,k,:) - A_0_old(i,j,k,:)) + & + dt/1800.0_rt * ((-35.0_rt*sqrt(6.0_rt) + 440.0_rt)*A_1_old(i,j,k,:) + & + (-169.0_rt*sqrt(6.0_rt) + 296.0_rt)*A_2_old(i,j,k,:) + & + (-16.0_rt + 24.0_rt*sqrt(6.0_rt))*A_3_old(i,j,k,:)) + enddo + enddo + enddo + + else if (m_start == 1) then + + do k = lo(3), hi(3) + do j = lo(2), hi(2) + do i = lo(1), hi(1) + k_n(i,j,k,:) = k_m(i,j,k,:) + & + dt_m * (A_m(i,j,k,:) - A_1_old(i,j,k,:)) + & + dt/150.0_rt * ((-12.0_rt + 17.0_rt*sqrt(6.0_rt))*A_1_old(i,j,k,:) + & + (12.0_rt + 17.0_rt*sqrt(6.0_rt))*A_2_old(i,j,k,:) + & + (-4.0_rt*sqrt(6.0_rt))*A_3_old(i,j,k,:)) + enddo + enddo + enddo + + else if (m_start == 2) then + + do k = lo(3), hi(3) + do j = lo(2), hi(2) + do i = lo(1), hi(1) + k_n(i,j,k,:) = k_m(i,j,k,:) + & + dt_m * (A_m(i,j,k,:) - A_2_old(i,j,k,:)) + & + dt/600.0_rt * ((168.0_rt - 73.0_rt*sqrt(6.0_rt))*A_1_old(i,j,k,:) + & + (120.0_rt + 5.0_rt*sqrt(6.0_rt))*A_2_old(i,j,k,:) + & + (72.0_rt + 8.0_rt*sqrt(6.0_rt))*A_3_old(i,j,k,:)) + enddo + enddo + enddo + + else + call castro_error("error in ca_sdc_update_advection_o4_radau -- should not be here") + endif + + end subroutine ca_sdc_update_advection_o4_radau #ifdef REACTIONS - subroutine ca_sdc_compute_C4(lo, hi, & - A_m, Amlo, Amhi, & - A_0_old, A0lo, A0hi, & - A_1_old, A1lo, A1hi, & - A_2_old, A2lo, A2hi, & - R_0_old, R0lo, R0hi, & - R_1_old, R1lo, R1hi, & - R_2_old, R2lo, R2hi, & - C, Clo, Chi, & - m_start) bind(C, name="ca_sdc_compute_C4") + subroutine ca_sdc_compute_C4_lobatto(lo, hi, & + dt_m, dt, & + A_m, Amlo, Amhi, & + A_0_old, A0lo, A0hi, & + A_1_old, A1lo, A1hi, & + A_2_old, A2lo, A2hi, & + R_0_old, R0lo, R0hi, & + R_1_old, R1lo, R1hi, & + R_2_old, R2lo, R2hi, & + C, Clo, Chi, & + m_start) bind(C, name="ca_sdc_compute_C4_lobatto") + ! compute the 'C' term for the 4th-order solve with reactions ! note: this 'C' is cell-averages @@ -974,6 +1129,7 @@ subroutine ca_sdc_compute_C4(lo, hi, & integer, intent(in) :: Clo(3), Chi(3) integer, intent(in) :: m_start + real(rt), intent(in) :: dt_m, dt real(rt), intent(in) :: A_m(Amlo(1):Amhi(1), Amlo(2):Amhi(2), Amlo(3):Amhi(3), NVAR) real(rt), intent(in) :: A_0_old(A0lo(1):A0hi(1), A0lo(2):A0hi(2), A0lo(3):A0hi(3), NVAR) real(rt), intent(in) :: A_1_old(A1lo(1):A1hi(1), A1lo(2):A1hi(2), A1lo(3):A1hi(3), NVAR) @@ -986,13 +1142,14 @@ subroutine ca_sdc_compute_C4(lo, hi, & integer :: i, j, k real(rt) :: integral(NVAR) + ! Gauss-Lobatto (Simpsons) + do k = lo(3), hi(3) do j = lo(2), hi(2) do i = lo(1), hi(1) - ! compute the integral (without the dt). Note that each of these is over - ! dt/2 if (m_start == 0) then + ! compute the integral from [t_m, t_{m+1}], normalized by dt_m integral(:) = ONE/12.0_rt * (FIVE*(A_0_old(i,j,k,:) + R_0_old(i,j,k,:)) + & EIGHT*(A_1_old(i,j,k,:) + R_1_old(i,j,k,:)) - & (A_2_old(i,j,k,:) + R_2_old(i,j,k,:))) @@ -1000,6 +1157,7 @@ subroutine ca_sdc_compute_C4(lo, hi, & C(i,j,k,:) = (A_m(i,j,k,:) - A_0_old(i,j,k,:)) - R_1_old(i,j,k,:) + integral else if (m_start == 1) then + ! compute the integral from [t_m, t_{m+1}], normalized by dt_m integral(:) = ONE/12.0_rt * (-(A_0_old(i,j,k,:) + R_0_old(i,j,k,:)) + & EIGHT*(A_1_old(i,j,k,:) + R_1_old(i,j,k,:)) + & FIVE*(A_2_old(i,j,k,:) + R_2_old(i,j,k,:))) @@ -1014,7 +1172,102 @@ subroutine ca_sdc_compute_C4(lo, hi, & enddo enddo - end subroutine ca_sdc_compute_C4 + end subroutine ca_sdc_compute_C4_lobatto + + + subroutine ca_sdc_compute_C4_radau(lo, hi, & + dt_m, dt, & + A_m, Amlo, Amhi, & + A_0_old, A0lo, A0hi, & + A_1_old, A1lo, A1hi, & + A_2_old, A2lo, A2hi, & + A_3_old, A3lo, A3hi, & + R_0_old, R0lo, R0hi, & + R_1_old, R1lo, R1hi, & + R_2_old, R2lo, R2hi, & + R_3_old, R3lo, R3hi, & + C, Clo, Chi, & + m_start) bind(C, name="ca_sdc_compute_C4_radau") + ! compute the 'C' term for the 4th-order solve with reactions + ! note: this 'C' is cell-averages + + use meth_params_module, only : NVAR + use amrex_constants_module, only : ONE, HALF, TWO, FIVE, EIGHT + + implicit none + + integer, intent(in) :: lo(3), hi(3) + integer, intent(in) :: Amlo(3), Amhi(3) + integer, intent(in) :: A0lo(3), A0hi(3) + integer, intent(in) :: A1lo(3), A1hi(3) + integer, intent(in) :: A2lo(3), A2hi(3) + integer, intent(in) :: A3lo(3), A3hi(3) + integer, intent(in) :: R0lo(3), R0hi(3) + integer, intent(in) :: R1lo(3), R1hi(3) + integer, intent(in) :: R2lo(3), R2hi(3) + integer, intent(in) :: R3lo(3), R3hi(3) + integer, intent(in) :: Clo(3), Chi(3) + integer, intent(in) :: m_start + + real(rt), intent(in) :: dt_m, dt + real(rt), intent(in) :: A_m(Amlo(1):Amhi(1), Amlo(2):Amhi(2), Amlo(3):Amhi(3), NVAR) + real(rt), intent(in) :: A_0_old(A0lo(1):A0hi(1), A0lo(2):A0hi(2), A0lo(3):A0hi(3), NVAR) + real(rt), intent(in) :: A_1_old(A1lo(1):A1hi(1), A1lo(2):A1hi(2), A1lo(3):A1hi(3), NVAR) + real(rt), intent(in) :: A_2_old(A2lo(1):A2hi(1), A2lo(2):A2hi(2), A2lo(3):A2hi(3), NVAR) + real(rt), intent(in) :: A_3_old(A3lo(1):A3hi(1), A3lo(2):A3hi(2), A3lo(3):A3hi(3), NVAR) + real(rt), intent(in) :: R_0_old(R0lo(1):R0hi(1), R0lo(2):R0hi(2), R0lo(3):R0hi(3), NVAR) + real(rt), intent(in) :: R_1_old(R1lo(1):R1hi(1), R1lo(2):R1hi(2), R1lo(3):R1hi(3), NVAR) + real(rt), intent(in) :: R_2_old(R2lo(1):R2hi(1), R2lo(2):R2hi(2), R2lo(3):R2hi(3), NVAR) + real(rt), intent(in) :: R_3_old(R3lo(1):R3hi(1), R3lo(2):R3hi(2), R3lo(3):R3hi(3), NVAR) + real(rt), intent(out) :: C(Clo(1):Chi(1), Clo(2):Chi(2), Clo(3):Chi(3), NVAR) + + integer :: i, j, k + real(rt) :: integral(NVAR) + + ! Gauss-Lobatto (Simpsons) + + do k = lo(3), hi(3) + do j = lo(2), hi(2) + do i = lo(1), hi(1) + + if (m_start == 0) then + ! compute the integral from [t_m, t_{m+1}], normalized by dt_m + integral(:) = (dt/dt_m) * (ONE/1800.0_rt) * & + ((-35.0_rt*sqrt(6.0_rt) + 440.0_rt)*(A_1_old(i,j,k,:) + R_1_old(i,j,k,:)) + & + (-169.0_rt*sqrt(6.0_rt) + 296.0_rt)*(A_2_old(i,j,k,:) + R_2_old(i,j,k,:)) + & + (-16.0_rt + 24.0_rt*sqrt(6.0_rt))*(A_3_old(i,j,k,:) + R_3_old(i,j,k,:))) + + C(i,j,k,:) = (A_m(i,j,k,:) - A_0_old(i,j,k,:)) - R_1_old(i,j,k,:) + integral + + else if (m_start == 1) then + ! compute the integral from [t_m, t_{m+1}], normalized by dt_m + integral(:) = (dt/dt_m) * (ONE/150.0_rt) * & + ((-12.0_rt + 17.0_rt*sqrt(6.0_rt))*(A_1_old(i,j,k,:) + R_1_old(i,j,k,:)) + & + (12.0_rt + 17.0_rt*sqrt(6.0_rt))*(A_2_old(i,j,k,:) + R_2_old(i,j,k,:)) + & + (-4.0_rt*sqrt(6.0_rt))*(A_3_old(i,j,k,:) + R_3_old(i,j,k,:))) + + C(i,j,k,:) = (A_m(i,j,k,:) - A_1_old(i,j,k,:)) - R_2_old(i,j,k,:) + integral + + + else if (m_start == 2) then + ! compute the integral from [t_m, t_{m+1}], normalized by dt_m + integral(:) = (dt/dt_m) * (ONE/600.0_rt) * & + ((168.0_rt - 73.0_rt*sqrt(6.0_rt))*(A_1_old(i,j,k,:) + R_1_old(i,j,k,:)) + & + (120.0_rt + 5.0_rt*sqrt(6.0_rt))*(A_2_old(i,j,k,:) + R_2_old(i,j,k,:)) + & + (72.0_rt + 8.0_rt*sqrt(6.0_rt))*(A_3_old(i,j,k,:) + R_3_old(i,j,k,:))) + + C(i,j,k,:) = (A_m(i,j,k,:) - A_2_old(i,j,k,:)) - R_3_old(i,j,k,:) + integral + + else + call castro_error("error in ca_sdc_compute_C4 -- should not be here") + endif + + end do + end do + end do + + end subroutine ca_sdc_compute_C4_radau + subroutine ca_sdc_compute_initial_guess(lo, hi, & U_old, Uo_lo, Uo_hi, & @@ -1046,6 +1299,8 @@ subroutine ca_sdc_compute_initial_guess(lo, hi, & integer, intent(in) :: sdc_iteration integer :: i, j, k + ! Here dt_m is the timestep to update from time node m to m+1 + do k = lo(3), hi(3) do j = lo(2), hi(2) do i = lo(1), hi(1) @@ -1063,18 +1318,158 @@ subroutine ca_sdc_compute_initial_guess(lo, hi, & end subroutine ca_sdc_compute_initial_guess + subroutine ca_sdc_compute_C2_lobatto(lo, hi, dt_m, dt, & + A_m, Amlo, Amhi, & + A_0_old, A0lo, A0hi, & + A_1_old, A1lo, A1hi, & + R_0_old, R0lo, R0hi, & + R_1_old, R1lo, R1hi, & + C, Clo, Chi, & + m_start) bind(C, name="ca_sdc_compute_C2_lobatto") + ! compute the source term C for the 2nd order Lobatto update + + ! Here, dt_m is the timestep between time-nodes m and m+1 + + use meth_params_module, only : NVAR + use amrex_constants_module, only : ZERO, HALF + use burn_type_module, only : burn_t + use network, only : nspec, nspec_evolve + use react_util_module + + implicit none + + integer, intent(in) :: lo(3), hi(3) + real(rt), intent(in) :: dt_m, dt + integer, intent(in) :: Amlo(3), Amhi(3) + integer, intent(in) :: A0lo(3), A0hi(3) + integer, intent(in) :: A1lo(3), A1hi(3) + integer, intent(in) :: R0lo(3), R0hi(3) + integer, intent(in) :: R1lo(3), R1hi(3) + integer, intent(in) :: Clo(3), Chi(3) + integer, intent(in) :: m_start + + real(rt), intent(in) :: A_m(Amlo(1):Amhi(1), Amlo(2):Amhi(2), Amlo(3):Amhi(3), NVAR) + real(rt), intent(in) :: A_0_old(A0lo(1):A0hi(1), A0lo(2):A0hi(2), A0lo(3):A0hi(3), NVAR) + real(rt), intent(in) :: A_1_old(A1lo(1):A1hi(1), A1lo(2):A1hi(2), A1lo(3):A1hi(3), NVAR) + + real(rt), intent(in) :: R_0_old(R0lo(1):R0hi(1), R0lo(2):R0hi(2), R0lo(3):R0hi(3), NVAR) + real(rt), intent(in) :: R_1_old(R1lo(1):R1hi(1), R1lo(2):R1hi(2), R1lo(3):R1hi(3), NVAR) + + real(rt), intent(out) :: C(Clo(1):Chi(1), Clo(2):Chi(2), Clo(3):Chi(3), NVAR) + + integer :: i, j, k + + do k = lo(3), hi(3) + do j = lo(2), hi(2) + do i = lo(1), hi(1) + + ! construct the source term to the update for 2nd order + ! Lobatto, there is no advective correction, and we have + ! C = - R(U^{m+1,k}) + I_m^{m+1}/dt + C(i,j,k,:) = -R_1_old(i,j,k,:) + & + HALF * (A_0_old(i,j,k,:) + A_1_old(i,j,k,:)) + & + HALF * (R_0_old(i,j,k,:) + R_1_old(i,j,k,:)) + end do + end do + end do + + end subroutine ca_sdc_compute_C2_lobatto + + + subroutine ca_sdc_compute_C2_radau(lo, hi, dt_m, dt, & + A_m, Amlo, Amhi, & + A_0_old, A0lo, A0hi, & + A_1_old, A1lo, A1hi, & + A_2_old, A2lo, A2hi, & + R_0_old, R0lo, R0hi, & + R_1_old, R1lo, R1hi, & + R_2_old, R2lo, R2hi, & + C, Clo, Chi, & + m_start) bind(C, name="ca_sdc_compute_C2_radau") + ! compute the source term C for the 2nd order Radau update + + ! Here, dt_m is the timestep between time-nodes m and m+1 + + use meth_params_module, only : NVAR + use amrex_constants_module, only : ZERO, HALF, ONE, FIVE + use burn_type_module, only : burn_t + use network, only : nspec, nspec_evolve + use react_util_module + + implicit none + + integer, intent(in) :: lo(3), hi(3) + real(rt), intent(in) :: dt_m, dt + integer, intent(in) :: Amlo(3), Amhi(3) + integer, intent(in) :: A0lo(3), A0hi(3) + integer, intent(in) :: A1lo(3), A1hi(3) + integer, intent(in) :: A2lo(3), A2hi(3) + integer, intent(in) :: R0lo(3), R0hi(3) + integer, intent(in) :: R1lo(3), R1hi(3) + integer, intent(in) :: R2lo(3), R2hi(3) + integer, intent(in) :: Clo(3), Chi(3) + integer, intent(in) :: m_start + + real(rt), intent(in) :: A_m(Amlo(1):Amhi(1), Amlo(2):Amhi(2), Amlo(3):Amhi(3), NVAR) + real(rt), intent(in) :: A_0_old(A0lo(1):A0hi(1), A0lo(2):A0hi(2), A0lo(3):A0hi(3), NVAR) + real(rt), intent(in) :: A_1_old(A1lo(1):A1hi(1), A1lo(2):A1hi(2), A1lo(3):A1hi(3), NVAR) + real(rt), intent(in) :: A_2_old(A2lo(1):A2hi(1), A2lo(2):A2hi(2), A2lo(3):A2hi(3), NVAR) + + real(rt), intent(in) :: R_0_old(R0lo(1):R0hi(1), R0lo(2):R0hi(2), R0lo(3):R0hi(3), NVAR) + real(rt), intent(in) :: R_1_old(R1lo(1):R1hi(1), R1lo(2):R1hi(2), R1lo(3):R1hi(3), NVAR) + real(rt), intent(in) :: R_2_old(R2lo(1):R2hi(1), R2lo(2):R2hi(2), R2lo(3):R2hi(3), NVAR) + + real(rt), intent(out) :: C(Clo(1):Chi(1), Clo(2):Chi(2), Clo(3):Chi(3), NVAR) + + integer :: i, j, k + + ! construct the source term to the update for 2nd order + ! Radau + + if (m_start == 0) then + do k = lo(3), hi(3) + do j = lo(2), hi(2) + do i = lo(1), hi(1) + + C(i,j,k,:) = -R_1_old(i,j,k,:) + & + (A_m(i,j,k,:) - A_0_old(i,j,k,:)) + & + (dt/dt_m) * (ONE/12.0_rt) * & + FIVE*(A_1_old(i,j,k,:) + R_1_old(i,j,k,:)) - & + (A_2_old(i,j,k,:) + R_2_old(i,j,k,:)) + end do + end do + end do + + else if (m_start == 1) then + do k = lo(3), hi(3) + do j = lo(2), hi(2) + do i = lo(1), hi(1) + + C(i,j,k,:) = -R_2_old(i,j,k,:) + & + (A_m(i,j,k,:) - A_1_old(i,j,k,:)) + & + (dt/dt_m) * (ONE/3.0_rt) * & + (A_1_old(i,j,k,:) + R_1_old(i,j,k,:)) + & + (A_2_old(i,j,k,:) + R_2_old(i,j,k,:)) + end do + end do + end do + end if + + end subroutine ca_sdc_compute_C2_radau + + subroutine ca_sdc_update_o2(lo, hi, dt_m, & k_m, kmlo, kmhi, & k_n, knlo, knhi, & A_m, Amlo, Amhi, & - A_0_old, A0lo, A0hi, & - A_1_old, A1lo, A1hi, & - R_0_old, R0lo, R0hi, & - R_1_old, R1lo, R1hi, & + R_m_old, Rmlo, Rmhi, & + C, Clo, Chi, & sdc_iteration, & m_start) bind(C, name="ca_sdc_update_o2") ! update k_m to k_n via advection -- this is a second-order accurate update + ! Here, dt_m is the timestep between time-nodes m and m+1 + use meth_params_module, only : NVAR use amrex_constants_module, only : ZERO, HALF use burn_type_module, only : burn_t @@ -1088,10 +1483,8 @@ subroutine ca_sdc_update_o2(lo, hi, dt_m, & integer, intent(in) :: kmlo(3), kmhi(3) integer, intent(in) :: knlo(3), knhi(3) integer, intent(in) :: Amlo(3), Amhi(3) - integer, intent(in) :: A0lo(3), A0hi(3) - integer, intent(in) :: A1lo(3), A1hi(3) - integer, intent(in) :: R0lo(3), R0hi(3) - integer, intent(in) :: R1lo(3), R1hi(3) + integer, intent(in) :: Rmlo(3), Rmhi(3) + integer, intent(in) :: Clo(3), Chi(3) integer, intent(in) :: sdc_iteration, m_start @@ -1099,32 +1492,21 @@ subroutine ca_sdc_update_o2(lo, hi, dt_m, & real(rt), intent(inout) :: k_n(knlo(1):knhi(1), knlo(2):knhi(2), knlo(3):knhi(3), NVAR) real(rt), intent(in) :: A_m(Amlo(1):Amhi(1), Amlo(2):Amhi(2), Amlo(3):Amhi(3), NVAR) - real(rt), intent(in) :: A_0_old(A0lo(1):A0hi(1), A0lo(2):A0hi(2), A0lo(3):A0hi(3), NVAR) - real(rt), intent(in) :: A_1_old(A1lo(1):A1hi(1), A1lo(2):A1hi(2), A1lo(3):A1hi(3), NVAR) - - real(rt), intent(in) :: R_0_old(R0lo(1):R0hi(1), R0lo(2):R0hi(2), R0lo(3):R0hi(3), NVAR) - real(rt), intent(in) :: R_1_old(R1lo(1):R1hi(1), R1lo(2):R1hi(2), R1lo(3):R1hi(3), NVAR) + real(rt), intent(in) :: R_m_old(Rmlo(1):Rmhi(1), Rmlo(2):Rmhi(2), Rmlo(3):Rmhi(3), NVAR) + real(rt), intent(in) :: C(Clo(1):Chi(1), Clo(2):Chi(2), Clo(3):Chi(3), NVAR) integer :: i, j, k type(burn_t) :: burn_state - real(rt) :: U_old(NVAR), U_new(NVAR), C(NVAR), R_full(NVAR) - + real(rt) :: U_old(NVAR), U_new(NVAR), R_full(NVAR), C_zone(NVAR) - ! now consider the reacting system do k = lo(3), hi(3) do j = lo(2), hi(2) do i = lo(1), hi(1) U_old(:) = k_m(i,j,k,:) - - ! construct the source term to the update - ! for 2nd order, there is no advective correction, and we have - ! C = - R(U^{m+1,k}) + I_m^{m+1}/dt - C(:) = -R_1_old(i,j,k,:) + & - HALF * (A_0_old(i,j,k,:) + A_1_old(i,j,k,:)) + & - HALF * (R_0_old(i,j,k,:) + R_1_old(i,j,k,:)) + C_zone(:) = C(i,j,k,:) ! only burn if we are within the temperature and density ! limits for burning @@ -1140,7 +1522,7 @@ subroutine ca_sdc_update_o2(lo, hi, dt_m, & ! the first iteration, let's try to extrapolate forward ! in time. if (sdc_iteration == 0) then - U_new(:) = U_old(:) + dt_m * A_m(i,j,k,:) + dt_m * R_0_old(i,j,k,:) + U_new(:) = U_old(:) + dt_m * A_m(i,j,k,:) + dt_m * R_m_old(i,j,k,:) else U_new(:) = k_n(i,j,k,:) endif @@ -1153,7 +1535,7 @@ subroutine ca_sdc_update_o2(lo, hi, dt_m, & end if - U_new(:) = U_old(:) + dt_m * R_full(:) + dt_m * C(:) + U_new(:) = U_old(:) + dt_m * R_full(:) + dt_m * C_zone(:) ! copy back to k_n k_n(i,j,k,:) = U_new(:) @@ -1174,7 +1556,7 @@ subroutine ca_sdc_update_centers_o4(lo, hi, dt_m, & ! update U_old to U_new on cell-centers. This is an implicit ! solve because of reactions. Here U_old corresponds to time node ! m and U_new is node m+1. dt_m is the timestep between m and - ! m+1, which is dt_m = dt/2. + ! m+1 use meth_params_module, only : NVAR use react_util_module, only : okay_to_burn diff --git a/Source/hydro/Castro_hydro_F.H b/Source/hydro/Castro_hydro_F.H index 17546fc893..bcc16859c0 100644 --- a/Source/hydro/Castro_hydro_F.H +++ b/Source/hydro/Castro_hydro_F.H @@ -503,9 +503,10 @@ extern "C" const amrex::Real* dx); #ifndef AMREX_USE_CUDA - void ca_sdc_update_advection_o2 + void ca_sdc_update_advection_o2_lobatto (const int* lo, const int* hi, - const amrex::Real* time, + const amrex::Real* dt_m, + const amrex::Real* dt, const BL_FORT_FAB_ARG_3D(k_m), BL_FORT_FAB_ARG_3D(k_n), BL_FORT_FAB_ARG_3D(A_m), @@ -513,9 +514,10 @@ extern "C" BL_FORT_FAB_ARG_3D(A_1_old), const int* m_start); - void ca_sdc_update_advection_o4 + void ca_sdc_update_advection_o2_radau (const int* lo, const int* hi, - const amrex::Real* time, + const amrex::Real* dt_m, + const amrex::Real* dt, const BL_FORT_FAB_ARG_3D(k_m), BL_FORT_FAB_ARG_3D(k_n), BL_FORT_FAB_ARG_3D(A_m), @@ -524,17 +526,66 @@ extern "C" BL_FORT_FAB_ARG_3D(A_2_old), const int* m_start); -#ifdef REACTIONS - void ca_sdc_update_o2 + void ca_sdc_update_advection_o4_lobatto (const int* lo, const int* hi, - const amrex::Real* time, + const amrex::Real* dt_m, + const amrex::Real* dt, + const BL_FORT_FAB_ARG_3D(k_m), + BL_FORT_FAB_ARG_3D(k_n), + BL_FORT_FAB_ARG_3D(A_m), + BL_FORT_FAB_ARG_3D(A_0_old), + BL_FORT_FAB_ARG_3D(A_1_old), + BL_FORT_FAB_ARG_3D(A_2_old), + const int* m_start); + + void ca_sdc_update_advection_o4_radau + (const int* lo, const int* hi, + const amrex::Real* dt_m, + const amrex::Real* dt, const BL_FORT_FAB_ARG_3D(k_m), BL_FORT_FAB_ARG_3D(k_n), BL_FORT_FAB_ARG_3D(A_m), BL_FORT_FAB_ARG_3D(A_0_old), BL_FORT_FAB_ARG_3D(A_1_old), - BL_FORT_FAB_ARG_3D(R_0_old), - BL_FORT_FAB_ARG_3D(R_1_old), + BL_FORT_FAB_ARG_3D(A_2_old), + BL_FORT_FAB_ARG_3D(A_3_old), + const int* m_start); + +#ifdef REACTIONS + void ca_sdc_compute_C2_lobatto + (const int* lo, const int* hi, + const amrex::Real* dt_m, + const amrex::Real* dt, + const BL_FORT_FAB_ARG_3D(A_m), + const BL_FORT_FAB_ARG_3D(A_0_old), + const BL_FORT_FAB_ARG_3D(A_1_old), + const BL_FORT_FAB_ARG_3D(R_0_old), + const BL_FORT_FAB_ARG_3D(R_1_old), + BL_FORT_FAB_ARG_3D(C2), + const int* m_start); + + void ca_sdc_compute_C2_radau + (const int* lo, const int* hi, + const amrex::Real* dt_m, + const amrex::Real* dt, + const BL_FORT_FAB_ARG_3D(A_m), + const BL_FORT_FAB_ARG_3D(A_0_old), + const BL_FORT_FAB_ARG_3D(A_1_old), + const BL_FORT_FAB_ARG_3D(A_2_old), + const BL_FORT_FAB_ARG_3D(R_0_old), + const BL_FORT_FAB_ARG_3D(R_1_old), + const BL_FORT_FAB_ARG_3D(R_2_old), + BL_FORT_FAB_ARG_3D(C2), + const int* m_start); + + void ca_sdc_update_o2 + (const int* lo, const int* hi, + const amrex::Real* dt_m, + const BL_FORT_FAB_ARG_3D(k_m), + BL_FORT_FAB_ARG_3D(k_n), + const BL_FORT_FAB_ARG_3D(A_m), + const BL_FORT_FAB_ARG_3D(R_m_old), + const BL_FORT_FAB_ARG_3D(C2), const int* sdc_iteration, const int* m_start); @@ -554,16 +605,33 @@ extern "C" const BL_FORT_FAB_ARG_3D(C), const BL_FORT_FAB_ARG_3D(R_new)); - void ca_sdc_compute_C4(const int* lo, const int* hi, - const BL_FORT_FAB_ARG_3D(A_m), - const BL_FORT_FAB_ARG_3D(A_0), - const BL_FORT_FAB_ARG_3D(A_1), - const BL_FORT_FAB_ARG_3D(A_2), - const BL_FORT_FAB_ARG_3D(R_0), - const BL_FORT_FAB_ARG_3D(R_1), - const BL_FORT_FAB_ARG_3D(R_2), - BL_FORT_FAB_ARG_3D(C), - const int* m_start); + void ca_sdc_compute_C4_lobatto(const int* lo, const int* hi, + const amrex::Real* dt_m, + const amrex::Real* dt, + const BL_FORT_FAB_ARG_3D(A_m), + const BL_FORT_FAB_ARG_3D(A_0), + const BL_FORT_FAB_ARG_3D(A_1), + const BL_FORT_FAB_ARG_3D(A_2), + const BL_FORT_FAB_ARG_3D(R_0), + const BL_FORT_FAB_ARG_3D(R_1), + const BL_FORT_FAB_ARG_3D(R_2), + BL_FORT_FAB_ARG_3D(C), + const int* m_start); + + void ca_sdc_compute_C4_radau(const int* lo, const int* hi, + const amrex::Real* dt_m, + const amrex::Real* dt, + const BL_FORT_FAB_ARG_3D(A_m), + const BL_FORT_FAB_ARG_3D(A_0), + const BL_FORT_FAB_ARG_3D(A_1), + const BL_FORT_FAB_ARG_3D(A_2), + const BL_FORT_FAB_ARG_3D(A_3), + const BL_FORT_FAB_ARG_3D(R_0), + const BL_FORT_FAB_ARG_3D(R_1), + const BL_FORT_FAB_ARG_3D(R_2), + const BL_FORT_FAB_ARG_3D(R_3), + BL_FORT_FAB_ARG_3D(C), + const int* m_start); void ca_sdc_compute_initial_guess(const int* lo, const int* hi, const BL_FORT_FAB_ARG_3D(U_old), diff --git a/Source/hydro/Castro_mol_hydro.cpp b/Source/hydro/Castro_mol_hydro.cpp index 65ba80b856..0bc85ce0fe 100644 --- a/Source/hydro/Castro_mol_hydro.cpp +++ b/Source/hydro/Castro_mol_hydro.cpp @@ -24,7 +24,7 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) const Real strt_time = ParallelDescriptor::second(); if (verbose && ParallelDescriptor::IOProcessor()) { - std::cout << "... SDC iteration: " << sdc_iteration << "; current node: " << current_sdc_node << std::endl; + std::cout << "... construct advection term, SDC iteration: " << sdc_iteration << "; current node: " << current_sdc_node << std::endl; }