diff --git a/CHANGES.md b/CHANGES.md index 107dd7b9a..459b30fe9 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,3 +1,8 @@ +# 23.08 + + * remove the SDC code paths. This is no longer supported by + Microphysics and is not being tested. + # 23.06 * the test_diffusion unit test was cleaned up and now gives the diff --git a/Source/Maestro.H b/Source/Maestro.H index c72dc433f..69ffa51c8 100644 --- a/Source/Maestro.H +++ b/Source/Maestro.H @@ -141,10 +141,6 @@ class Maestro : public amrex::AmrCore { /// @param is_initIter is it the initial iteration? void AdvanceTimeStepAverage(bool is_initIter); - // advance solution at all levels for a single time step using SDC - // instead of Strang splitting - void AdvanceTimeStepSDC(bool is_initIter); - // end MaestroAdvance.cpp functions //////////// @@ -498,19 +494,6 @@ class Maestro : public amrex::AmrCore { const amrex::Vector>& w0mac, const BaseState& rho0_predicted_edge); - // SDC - void DensityAdvanceSDC( - int which_step, amrex::Vector& scalold, - amrex::Vector& scalnew, - amrex::Vector>& sedge, - amrex::Vector>& sflux, - amrex::Vector& scal_force, - amrex::Vector& etarhoflux, - amrex::Vector>& umac, - const amrex::Vector>& w0mac, - const BaseState& rho0_predicted_edge); - //////////////////////// - //////////// // MaestroDiag.cpp functions @@ -577,18 +560,6 @@ class Maestro : public amrex::AmrCore { const amrex::Vector>& w0mac, const amrex::Vector& thermal); - // SDC - void EnthalpyAdvanceSDC( - int which_step, amrex::Vector& scalold, - amrex::Vector& scalnew, - amrex::Vector>& sedge, - amrex::Vector>& sflux, - amrex::Vector& scal_force, - amrex::Vector>& umac, - const amrex::Vector>& w0mac, - const amrex::Vector& thermal); - //////////////////////// - //////////////////////// // MaestroFillData.cpp functions @@ -809,8 +780,6 @@ class Maestro : public amrex::AmrCore { /// Performs the divu iteration void DivuIter(int istep_divu_iter); - // SDC - void DivuIterSDC(int istep_divu_iter); /// Performs the initial iteration to initialize `gradpi` void InitIter(); @@ -1436,14 +1405,6 @@ class Maestro : public amrex::AmrCore { const BaseState& p0, const amrex::Real dt_in, const amrex::Real time_in); - void ReactSDC(const amrex::Vector& s_in, - amrex::Vector& s_out, - amrex::Vector& rho_Hext, - const BaseState& p0, const amrex::Real dt_in, - const amrex::Real time_in, - amrex::Vector& source); - -#ifndef SDC void Burner(const amrex::Vector& s_in, amrex::Vector& s_out, const amrex::Vector& rho_Hext, @@ -1452,19 +1413,7 @@ class Maestro : public amrex::AmrCore { const BaseState& p0, const amrex::Real dt_in, const amrex::Real time_in); -#else - void Burner(const amrex::Vector& s_in, - amrex::Vector& s_out, - const BaseState& p0, const amrex::Real dt_in, - const amrex::Real time_in, - const amrex::Vector& source); -#endif - // compute heating terms, rho_omegadot and rho_Hnuc - void MakeReactionRates(amrex::Vector& rho_omegadot, - amrex::Vector& rho_Hnuc, - const amrex::Vector& scal); - void MakeIntraCoeffs(const amrex::Vector& scal1, const amrex::Vector& scal2, amrex::Vector& cp, @@ -1714,17 +1663,6 @@ class Maestro : public amrex::AmrCore { const amrex::Vector& Xkcoeff2, const amrex::Vector& pcoeff2); - void ThermalConductSDC(int which_step, - const amrex::Vector& s_old, - amrex::Vector& s_hat, - const amrex::Vector& s_new, - const amrex::Vector& hcoeff1, - const amrex::Vector& Xkcoeff1, - const amrex::Vector& pcoeff1, - const amrex::Vector& hcoeff2, - const amrex::Vector& Xkcoeff2, - const amrex::Vector& pcoeff2); - void MakeExplicitThermalHterm(amrex::Vector& thermal, const amrex::Vector& scal, const amrex::Vector& hcoeff); @@ -1942,7 +1880,6 @@ class Maestro : public amrex::AmrCore { amrex::Vector gpi; amrex::Vector dSdt; amrex::Vector pi; // nodal - amrex::Vector intra; // for sdc /// this doesn't have to be persistent, but we make it so that we avoid /// continually creating and filling temporaries diff --git a/Source/MaestroAdvanceSdc.cpp b/Source/MaestroAdvanceSdc.cpp deleted file mode 100644 index 5920274a2..000000000 --- a/Source/MaestroAdvanceSdc.cpp +++ /dev/null @@ -1,1232 +0,0 @@ - -#include - -using namespace amrex; - -// advance a single level for a single time step, updates flux registers -void Maestro::AdvanceTimeStepSDC(bool is_initIter) { - // timer for profiling - BL_PROFILE_VAR("Maestro::AdvanceTimeStepSDC()", AdvanceTimeStepSDC); - - // cell-centered MultiFabs needed within the AdvanceTimeStep routine - Vector shat(finest_level + 1); - Vector rhohalf(finest_level + 1); - Vector cphalf(finest_level + 1); - Vector xihalf(finest_level + 1); - Vector macrhs(finest_level + 1); - Vector macphi(finest_level + 1); - Vector S_cc_nph(finest_level + 1); - Vector rho_omegadot(finest_level + 1); - Vector diff_old(finest_level + 1); - Vector diff_new(finest_level + 1); - Vector diff_hat(finest_level + 1); - Vector diff_hterm_new(finest_level + 1); - Vector diff_hterm_hat(finest_level + 1); - Vector rho_Hnuc(finest_level + 1); - Vector rho_Hext(finest_level + 1); - Vector sdc_source(finest_level + 1); - Vector aofs(finest_level + 1); - Vector intra_rhoh0(finest_level + 1); - - Vector delta_gamma1_term(finest_level + 1); - Vector delta_gamma1(finest_level + 1); - Vector p0_cart(finest_level + 1); - Vector delta_p_term(finest_level + 1); - - Vector Tcoeff1(finest_level + 1); - Vector hcoeff1(finest_level + 1); - Vector Xkcoeff1(finest_level + 1); - Vector pcoeff1(finest_level + 1); - Vector Tcoeff2(finest_level + 1); - Vector hcoeff2(finest_level + 1); - Vector Xkcoeff2(finest_level + 1); - Vector pcoeff2(finest_level + 1); - Vector scal_force(finest_level + 1); - Vector delta_chi(finest_level + 1); - Vector sponge(finest_level + 1); - - // face-centered in the dm-direction (planar only) - Vector etarhoflux_dummy(finest_level + 1); - - // face-centered - Vector > umac(finest_level + 1); - Vector > sedge(finest_level + 1); - Vector > sflux(finest_level + 1); - - //////////////////////// - // needed for spherical routines only - - // cell-centered - Vector w0_force_cart_dummy(finest_level + 1); - - // face-centered - Vector > w0mac(finest_level + 1); - Vector > w0mac_dummy(finest_level + 1); - - // end spherical-only MultiFabs - //////////////////////// - - // vectors store the multilevel 1D states as one very long array - // these are cell-centered - BaseState grav_cell_nph(base_geom.max_radial_level + 1, - base_geom.nr_fine); - BaseState rho0_nph(base_geom.max_radial_level + 1, base_geom.nr_fine); - BaseState p0_nph(base_geom.max_radial_level + 1, base_geom.nr_fine); - BaseState p0_minus_peosbar(base_geom.max_radial_level + 1, - base_geom.nr_fine); - BaseState peosbar(base_geom.max_radial_level + 1, base_geom.nr_fine); - BaseState w0_force_dummy(base_geom.max_radial_level + 1, - base_geom.nr_fine); - BaseState Sbar(base_geom.max_radial_level + 1, base_geom.nr_fine); - BaseState beta0_nph(base_geom.max_radial_level + 1, - base_geom.nr_fine); - BaseState gamma1bar_nph(base_geom.max_radial_level + 1, - base_geom.nr_fine); - BaseState delta_gamma1_termbar(base_geom.max_radial_level + 1, - base_geom.nr_fine); - BaseState delta_rhoh0(base_geom.max_radial_level + 1, - base_geom.nr_fine); - - // vectors store the multilevel 1D states as one very long array - // these are edge-centered - BaseState w0_old(base_geom.max_radial_level + 1, - base_geom.nr_fine + 1); - BaseState rho0_pred_edge_dummy(base_geom.max_radial_level + 1, - base_geom.nr_fine + 1); - - bool is_predictor; - bool split_projection = true; - - Print() << "\nTimestep " << istep << " starts with TIME = " << t_old - << " DT = " << dt << std::endl - << std::endl; - - if (maestro_verbose > 0) { - Print() << "Cell Count:" << std::endl; - for (int lev = 0; lev <= finest_level; ++lev) { - Print() << "Level " << lev << ", " << CountCells(lev) << " cells" - << std::endl; - } - } - - for (int lev = 0; lev <= finest_level; ++lev) { - // cell-centered MultiFabs - shat[lev].define(grids[lev], dmap[lev], Nscal, ng_s); - rhohalf[lev].define(grids[lev], dmap[lev], 1, 1); - cphalf[lev].define(grids[lev], dmap[lev], 1, 0); - xihalf[lev].define(grids[lev], dmap[lev], NumSpec, 0); - macrhs[lev].define(grids[lev], dmap[lev], 1, 0); - macphi[lev].define(grids[lev], dmap[lev], 1, 1); - S_cc_nph[lev].define(grids[lev], dmap[lev], 1, 0); - rho_omegadot[lev].define(grids[lev], dmap[lev], NumSpec, 0); - diff_old[lev].define(grids[lev], dmap[lev], 1, 0); - diff_new[lev].define(grids[lev], dmap[lev], 1, 0); - diff_hat[lev].define(grids[lev], dmap[lev], 1, 0); - diff_hterm_new[lev].define(grids[lev], dmap[lev], 1, 0); - diff_hterm_hat[lev].define(grids[lev], dmap[lev], 1, 0); - rho_Hnuc[lev].define(grids[lev], dmap[lev], 1, 0); - rho_Hext[lev].define(grids[lev], dmap[lev], 1, 0); - sdc_source[lev].define(grids[lev], dmap[lev], Nscal, 0); - aofs[lev].define(grids[lev], dmap[lev], Nscal, 0); - intra_rhoh0[lev].define(grids[lev], dmap[lev], 1, 0); - delta_gamma1_term[lev].define(grids[lev], dmap[lev], 1, 0); - delta_gamma1[lev].define(grids[lev], dmap[lev], 1, 0); - p0_cart[lev].define(grids[lev], dmap[lev], 1, 0); - delta_p_term[lev].define(grids[lev], dmap[lev], 1, 0); - Tcoeff1[lev].define(grids[lev], dmap[lev], 1, 1); - hcoeff1[lev].define(grids[lev], dmap[lev], 1, 1); - Xkcoeff1[lev].define(grids[lev], dmap[lev], NumSpec, 1); - pcoeff1[lev].define(grids[lev], dmap[lev], 1, 1); - Tcoeff2[lev].define(grids[lev], dmap[lev], 1, 1); - hcoeff2[lev].define(grids[lev], dmap[lev], 1, 1); - Xkcoeff2[lev].define(grids[lev], dmap[lev], NumSpec, 1); - pcoeff2[lev].define(grids[lev], dmap[lev], 1, 1); - - if (ppm_trace_forces == 0) { - scal_force[lev].define(grids[lev], dmap[lev], Nscal, 1); - } else { - // we need more ghostcells if we are tracing the forces - scal_force[lev].define(grids[lev], dmap[lev], Nscal, ng_s); - } - delta_chi[lev].define(grids[lev], dmap[lev], 1, 0); - sponge[lev].define(grids[lev], dmap[lev], 1, 0); - - // face-centered in the dm-direction (planar only) - AMREX_D_TERM( - etarhoflux_dummy[lev].define(convert(grids[lev], nodal_flag_x), - dmap[lev], 1, 1); - , etarhoflux_dummy[lev].define(convert(grids[lev], nodal_flag_y), - dmap[lev], 1, 1); - , etarhoflux_dummy[lev].define(convert(grids[lev], nodal_flag_z), - dmap[lev], 1, 1);); - - // face-centered arrays of MultiFabs - AMREX_D_TERM(umac[lev][0].define(convert(grids[lev], nodal_flag_x), - dmap[lev], 1, 1); - , umac[lev][1].define(convert(grids[lev], nodal_flag_y), - dmap[lev], 1, 1); - , umac[lev][2].define(convert(grids[lev], nodal_flag_z), - dmap[lev], 1, 1);); - AMREX_D_TERM(sedge[lev][0].define(convert(grids[lev], nodal_flag_x), - dmap[lev], Nscal, 0); - , sedge[lev][1].define(convert(grids[lev], nodal_flag_y), - dmap[lev], Nscal, 0); - , sedge[lev][2].define(convert(grids[lev], nodal_flag_z), - dmap[lev], Nscal, 0);); - AMREX_D_TERM(sflux[lev][0].define(convert(grids[lev], nodal_flag_x), - dmap[lev], Nscal, 0); - , sflux[lev][1].define(convert(grids[lev], nodal_flag_y), - dmap[lev], Nscal, 0); - , sflux[lev][2].define(convert(grids[lev], nodal_flag_z), - dmap[lev], Nscal, 0);); - - // initialize umac - for (int d = 0; d < AMREX_SPACEDIM; ++d) { - umac[lev][d].setVal(0.); - } - - // initialize intra_rhoh0 - intra_rhoh0[lev].setVal(0.); - - rho_Hext[lev].setVal(0.); - } - -#if (AMREX_SPACEDIM == 3) - for (int lev = 0; lev <= finest_level; ++lev) { - w0mac[lev][0].define(convert(grids[lev], nodal_flag_x), dmap[lev], 1, - 1); - w0mac[lev][1].define(convert(grids[lev], nodal_flag_y), dmap[lev], 1, - 1); - w0mac[lev][2].define(convert(grids[lev], nodal_flag_z), dmap[lev], 1, - 1); - w0mac_dummy[lev][0].define(convert(grids[lev], nodal_flag_x), dmap[lev], - 1, 1); - w0mac_dummy[lev][1].define(convert(grids[lev], nodal_flag_y), dmap[lev], - 1, 1); - w0mac_dummy[lev][2].define(convert(grids[lev], nodal_flag_z), dmap[lev], - 1, 1); - } -#endif - - for (int lev = 0; lev <= finest_level; ++lev) { - w0_force_cart_dummy[lev].define(grids[lev], dmap[lev], AMREX_SPACEDIM, - 1); - w0_force_cart_dummy[lev].setVal(0.); - } - - // set etarhoflux_dummy to zero - for (int lev = 0; lev <= finest_level; ++lev) { - etarhoflux_dummy[lev].setVal(0.); - } - -#if (AMREX_SPACEDIM == 3) - // initialize MultiFabs and Vectors to ZERO - for (int lev = 0; lev <= finest_level; ++lev) { - for (int d = 0; d < AMREX_SPACEDIM; ++d) { - w0mac[lev][d].setVal(0.); - w0mac_dummy[lev][d].setVal(0.); - } - } -#endif - - // initialize to zero - Sbar.setVal(0.); - w0.setVal(0.0); - - // set dummy variables to zero - w0_force_dummy.setVal(0.0); - rho0_pred_edge_dummy.setVal(0.0); - - // make the sponge for all levels - if (do_sponge) { - SpongeInit(rho0_old); - MakeSponge(sponge); - } - - ////////////////////////////////////////////////////////////////////////////// - // STEP 1 -- Compute advection velocities - ////////////////////////////////////////////////////////////////////////////// - - if (maestro_verbose >= 1) { - Print() << "<<< STEP 1 : Compute advection velocities >>>" << std::endl; - } - - if (t_old == 0.) { - // this is either a pressure iteration or the first time step - // set S_cc_nph = (1/2) (S_cc_old + S_cc_new) - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::LinComb(S_cc_nph[lev], 0.5, S_cc_old[lev], 0, 0.5, - S_cc_new[lev], 0, 0, 1, 0); - } - } else { - // set S_cc_nph = S_cc_old + (dt/2) * dSdt - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::LinComb(S_cc_nph[lev], 1.0, S_cc_old[lev], 0, 0.5 * dt, - dSdt[lev], 0, 0, 1, 0); - } - } - // no ghost cells for S_cc_nph - AverageDown(S_cc_nph, 0, 1); - - // compute p0_minus_peosbar = p0_old - peosbar_old (for making w0) and - // compute delta_p_term = peos_old - p0_old (for RHS of projections) - if (dpdt_factor > 0.0) { - // peos_old (delta_p_term) now holds the thermodynamic p computed from sold(rho,h,X) - PfromRhoH(sold, sold, delta_p_term); - - // compute peosbar = Avg(peos_old) - Average(delta_p_term, peosbar, 0); - - // compute p0_minus_peosbar = p0_old - peosbar - p0_minus_peosbar.copy(p0_old - peosbar); - - // put p0_old on cart - Put1dArrayOnCart(p0_old, p0_cart, false, false, bcs_f, 0); - - // compute delta_p_term = peos_old - p0_old - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Subtract(delta_p_term[lev], p0_cart[lev], 0, 0, 1, 0); - } - - } else { - // these should have no effect if dpdt_factor <= 0 - p0_minus_peosbar.setVal(0.0); - for (int lev = 0; lev <= finest_level; ++lev) { - delta_p_term[lev].setVal(0.); - } - } - - if (evolve_base_state) { - if (split_projection) { - // compute Sbar = average(S_cc_nph) - Average(S_cc_nph, Sbar, 0); - - // save old-time value - w0_old.copy(w0); - - // compute w0, w0_force - is_predictor = true; - Makew0(w0_old, w0_force_dummy, Sbar, rho0_old, rho0_old, p0_old, - p0_old, gamma1bar_old, gamma1bar_old, p0_minus_peosbar, dt, - dtold, is_predictor); - - Put1dArrayOnCart(w0, w0_cart, true, true, bcs_u, 0, 1); -#if (AMREX_SPACEDIM == 3) - if (spherical) { - // put w0 on Cartesian edges - MakeW0mac(w0mac); - } -#endif - } - } - - // compute unprojected MAC velocities - is_predictor = true; - AdvancePremac(umac, w0mac_dummy, w0_force_cart_dummy); - - for (int lev = 0; lev <= finest_level; ++lev) { - delta_chi[lev].setVal(0.); - macphi[lev].setVal(0.); - delta_gamma1_term[lev].setVal(0.); - } - - if (evolve_base_state && !split_projection) { - Sbar.copy(1.0 / (gamma1bar_old * p0_old) * (p0_old - p0_nm1) / dtold); - } - - // compute RHS for MAC projection, beta0*(S_cc-Sbar) + beta0*delta_chi - MakeRHCCforMacProj(macrhs, rho0_old, S_cc_nph, Sbar, beta0_old, - delta_gamma1_term, gamma1bar_old, p0_old, delta_p_term, - delta_chi, is_predictor); - - if (spherical && evolve_base_state && split_projection) { - // subtract w0mac from umac - Addw0(umac, w0mac, -1.); - } - - // wallclock time - Real start_total_macproj = ParallelDescriptor::second(); - - // MAC projection - // includes spherical option in C++ function - MacProj(umac, macphi, macrhs, beta0_old, is_predictor); - - // wallclock time - Real end_total_macproj = ParallelDescriptor::second() - start_total_macproj; - ParallelDescriptor::ReduceRealMax(end_total_macproj, - ParallelDescriptor::IOProcessorNumber()); - - if (spherical && evolve_base_state && split_projection) { - // add w0mac back to umac - Addw0(umac, w0mac, 1.); - } - - ////////////////////////////////////////////////////////////////////////////// - // STEP 2 -- Predictor - ////////////////////////////////////////////////////////////////////////////// - - if (maestro_verbose >= 1) { - Print() << "<<< STEP 2 : Predictor >>>" << std::endl; - } - - ////////////////////////////////////////////////////////////////////////////// - // STEP 2A -- compute advective flux divergences - ////////////////////////////////////////////////////////////////////////////// - - if (maestro_verbose >= 1) { - Print() << "<<< STEP 2A : compute advective flux div >>>" << std::endl; - } - - // no need to advect the base state density - rho0_new.copy(rho0_old); - - // set diff to zero - for (int lev = 0; lev <= finest_level; ++lev) { - diff_old[lev].setVal(0.); - diff_new[lev].setVal(0.); - diff_hat[lev].setVal(0.); - } - - // diff is the forcing for rhoh or temperature - if (use_thermal_diffusion) { - MakeThermalCoeffs(sold, Tcoeff1, hcoeff1, Xkcoeff1, pcoeff1); - - MakeExplicitThermal(diff_old, sold, Tcoeff1, hcoeff1, Xkcoeff1, pcoeff1, - p0_old, temp_diffusion_formulation); - } - - // copy sold into shat - // temperature will be overwritten later after enthalpy advance - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Copy(shat[lev], sold[lev], 0, 0, Nscal, 0); - } - - if (maestro_verbose >= 1) { - Print() << " : density_advance >>>" << std::endl; - } - - // set sedge and sflux to zero - for (int lev = 0; lev <= finest_level; ++lev) { - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - sedge[lev][idim].setVal(0.); - sflux[lev][idim].setVal(0.); - } - } - - // advect rhoX and rho - DensityAdvanceSDC(1, sold, shat, sedge, sflux, scal_force, etarhoflux_dummy, - umac, w0mac, rho0_pred_edge_dummy); - - if (evolve_base_state) { - // correct the base state density by "averaging" - Average(shat, rho0_new, Rho); - ComputeCutoffCoords(rho0_new); - } - - // update grav_cell_new - if (evolve_base_state) { - MakeGravCell(grav_cell_new, rho0_new); - } else { - grav_cell_new.copy(grav_cell_old); - } - - // base state pressure update - if (evolve_base_state) { - // set new p0 through HSE - p0_new.copy(p0_old); - - EnforceHSE(rho0_new, p0_new, grav_cell_new); - - // compute p0_nph - p0_nph.copy(0.5 * (p0_old + p0_new)); - - // hold dp0/dt in psi for enthalpy advance - psi.copy((p0_new - p0_old) / dt); - } else { - p0_new.copy(p0_old); - p0_nph.copy(p0_old); - } - - // base state enthalpy update - if (evolve_base_state) { - // compute rhoh0_old by "averaging" - Average(sold, rhoh0_old, RhoH); - } else { - rhoh0_new.copy(rhoh0_old); - } - - if (maestro_verbose >= 1) { - Print() << " : enthalpy_advance >>>" << std::endl; - } - - EnthalpyAdvanceSDC(1, sold, shat, sedge, sflux, scal_force, umac, w0mac, - diff_old); - - // base state enthalpy update - if (evolve_base_state) { - // compute rhoh0_new by "averaging" - Average(shat, rhoh0_new, RhoH); - - // store (rhoh0_hat - rhoh0_old)/dt in delta_rhoh0 - delta_rhoh0.copy((rhoh0_new - rhoh0_old) / dt); - } - - // extract aofs = (shat - sold) / dt - intra - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::LinComb(aofs[lev], 1.0 / dt, shat[lev], 0, -1.0 / dt, - sold[lev], 0, 0, Nscal, 0); - MultiFab::Subtract(aofs[lev], intra[lev], 0, 0, Nscal, 0); - } - - ////////////////////////////////////////////////////////////////////////////// - // STEP 2B (optional) -- compute diffusive flux divergence - ////////////////////////////////////////////////////////////////////////////// - - if (maestro_verbose >= 1) { - Print() << "<<< STEP 2B : compute diffusive flux div >>>" << std::endl; - } - - if (use_thermal_diffusion) { - // 1 = predictor, 2 = corrector - ThermalConductSDC(1, sold, shat, snew, hcoeff1, - Xkcoeff1, pcoeff1, hcoeff2, Xkcoeff2, pcoeff2); - - // note p0_new => p0_hat if evolve_base_state = T - MakeExplicitThermal(diff_hat, shat, Tcoeff1, hcoeff1, Xkcoeff1, pcoeff1, - p0_new, temp_diffusion_formulation); - } - - ////////////////////////////////////////////////////////////////////////////// - // STEP 2C -- advance thermodynamic variables - ////////////////////////////////////////////////////////////////////////////// - - if (maestro_verbose >= 1) { - Print() << "<<< STEP 2C : advance thermo variables >>>" << std::endl; - } - - // build sdc_source - for (int lev = 0; lev <= finest_level; ++lev) { - sdc_source[lev].setVal(0.); - MultiFab::Add(sdc_source[lev], aofs[lev], FirstSpec, FirstSpec, NumSpec, - 0); - MultiFab::LinComb(sdc_source[lev], 0.5, diff_old[lev], 0, 0.5, - diff_hat[lev], 0, RhoH, 1, 0); - MultiFab::Add(sdc_source[lev], aofs[lev], RhoH, RhoH, 1, 0); - } - - // wallclock time - Real start_total_react = ParallelDescriptor::second(); - - ReactSDC(sold, snew, rho_Hext, p0_old, dt, t_old, sdc_source); - - // wallclock time - Real end_total_react = ParallelDescriptor::second() - start_total_react; - ParallelDescriptor::ReduceRealMax(end_total_react, - ParallelDescriptor::IOProcessorNumber()); - - // extract IR = [ (snew - sold)/dt - sdc_source ] - - for (int lev = 0; lev <= finest_level; ++lev) { - intra[lev].setVal(0.); - // species source term - MultiFab::LinComb(intra[lev], 1.0 / dt, snew[lev], FirstSpec, -1.0 / dt, - sold[lev], FirstSpec, FirstSpec, NumSpec, 0); - MultiFab::Subtract(intra[lev], sdc_source[lev], FirstSpec, FirstSpec, - NumSpec, 0); - // enthalpy source term - MultiFab::LinComb(intra[lev], 1.0 / dt, snew[lev], RhoH, -1.0 / dt, - sold[lev], RhoH, RhoH, 1, 0); - MultiFab::Subtract(intra[lev], sdc_source[lev], RhoH, RhoH, 1, 0); - } - - // massage the rhoh intra term into the proper form, depending on - // what we are predicting. Note: we do this before we deal with - // the species terms, since some enthalpy types need this default - // species intra. - - // first create rhohalf -- a lot of forms need this. - FillPatch(0.5 * (t_old + t_new), rhohalf, sold, snew, Rho, 0, 1, Rho, - bcs_s); - - if (evolve_base_state) { - // update base state density and pressure - Average(snew, rho0_new, Rho); - ComputeCutoffCoords(rho0_new); - - if (use_etarho) { - // compute the new etarho - if (!spherical) { - MakeEtarho(etarhoflux_dummy); -#if AMREX_SPACEDIM == 3 - } else { - MakeEtarhoSphr(sold, snew, umac, w0mac_dummy); -#endif - } - } - - MakeGravCell(grav_cell_new, rho0_new); - - EnforceHSE(rho0_new, p0_new, grav_cell_new); - - // hold dp0/dt in psi for Make_S_cc - psi.copy((p0_new - p0_old) / dt); - - // update base state enthalpy - Average(snew, rhoh0_new, RhoH); - - // compute intra_rhoh0 = (rhoh0_new - rhoh0_old)/dt - // - (rhoh0_hat - rhoh0_old)/dt - delta_rhoh0.copy((rhoh0_new - rhoh0_old) / dt - delta_rhoh0); - Put1dArrayOnCart(delta_rhoh0, intra_rhoh0, false, false, bcs_f, 0); - } - - // now update temperature - if (use_tfromp) { - TfromRhoP(snew, p0_new); - } else { - TfromRhoH(snew, p0_new); - } - - if (enthalpy_pred_type == predict_rhohprime) { - // intra is only different from predict_rhoh if rhoh0 is not constant - if (evolve_base_state) { - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Subtract(intra[lev], intra_rhoh0[lev], 0, RhoH, 1, 0); - } - } - - } else if (enthalpy_pred_type == predict_h) { - // we want this in terms of h, not (rho h) - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Divide(intra[lev], rhohalf[lev], 0, RhoH, 1, 0); - } - - } else if ((enthalpy_pred_type == predict_T_then_rhohprime) || - (enthalpy_pred_type == predict_T_then_h)) { - // for predict_T_*, the intra force needs to be in the temp_comp - // slot, since temperature is what is predicted. - - // first make the thermodynamic coefficients at the half-time - MakeIntraCoeffs(sold, snew, cphalf, xihalf); - - // overwrite intra(temp_comp). We want to create - // I_T = (1 / (rho c_p)) [ (rhoh_new - rhoh_old)/dt - A_rhoh - - // sum_k xi_k ( (rhoX_new - rhoX_old)/dt - A_rhoX ) ] - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Copy(intra[lev], intra[lev], RhoH, Temp, 1, 0); - - for (int comp = 0; comp < NumSpec; ++comp) { - // multiple xi by intra and store in xi - MultiFab::Multiply(xihalf[lev], intra[lev], FirstSpec + comp, - comp, 1, 0); - - // subtract from intra temp - MultiFab::Subtract(intra[lev], xihalf[lev], comp, Temp, 1, 0); - } - - MultiFab::Divide(intra[lev], rhohalf[lev], 0, Temp, 1, 0); - MultiFab::Divide(intra[lev], cphalf[lev], 0, Temp, 1, 0); - } - } - - // for some species_pred_types, we need to make intra in terms of - // X, NOT rhoX - if ((species_pred_type == predict_rhoprime_and_X) || - (species_pred_type == predict_rho_and_X)) { - for (int lev = 0; lev <= finest_level; ++lev) { - for (int comp = FirstSpec; comp < FirstSpec + NumSpec; ++comp) { - MultiFab::Divide(intra[lev], rhohalf[lev], 0, comp, 1, 0); - } - } - } - - // compute new-time coefficients and diffusion term - if (use_thermal_diffusion) { - MakeThermalCoeffs(snew, Tcoeff2, hcoeff2, Xkcoeff2, pcoeff2); - - MakeExplicitThermal(diff_new, snew, Tcoeff2, hcoeff2, Xkcoeff2, pcoeff2, - p0_new, temp_diffusion_formulation); - } - - if (evolve_base_state) { - // compute beta0 and gamma1bar - MakeGamma1bar(snew, gamma1bar_new, p0_new); - - MakeBeta0(beta0_new, rho0_new, p0_new, gamma1bar_new, grav_cell_new); - } else { - // Just pass beta0 and gamma1bar through if not evolving base state - beta0_new.copy(beta0_old); - gamma1bar_new.copy(gamma1bar_old); - } - - gamma1bar_nph.copy(0.5 * (gamma1bar_old + gamma1bar_new)); - beta0_nph.copy(0.5 * (beta0_old + beta0_new)); - - ////////////////////////////////////////////////////////////////////////////// - // Corrector loop - ////////////////////////////////////////////////////////////////////////////// - - for (int misdc = 0; misdc < sdc_iters; ++misdc) { - ////////////////////////////////////////////////////////////////////////////// - // STEP 3 -- Update advection velocities - ////////////////////////////////////////////////////////////////////////////// - - if (sdc_couple_mac_velocity) { - if (maestro_verbose >= 1) { - Print() - << "<<< STEP 3 : Update advection velocities (MISDC iter = " - << misdc << ") >>>" << std::endl; - } - - MakeReactionRates(rho_omegadot, rho_Hnuc, snew); - - // compute S at cell-centers - Make_S_cc(S_cc_new, delta_gamma1_term, delta_gamma1, snew, uold, - rho_omegadot, rho_Hnuc, rho_Hext, diff_new, p0_new, - gamma1bar_new, delta_gamma1_termbar); - - // set S_cc_nph = (1/2) (S_cc_old + S_cc_new) - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::LinComb(S_cc_nph[lev], 0.5, S_cc_old[lev], 0, 0.5, - S_cc_new[lev], 0, 0, 1, 0); - } - AverageDown(S_cc_nph, 0, 1); - - // compute p0_minus_peosbar = p0_new - peosbar_new (for making w0) and - // set delta_p_term = peos_new - p0_new (for RHS of projection) - if (dpdt_factor > 0.) { - // peos now holds "peos_new", the thermodynamic p computed from snew(rho,h,X) - PfromRhoH(snew, snew, delta_p_term); - - // compute peosbar = Avg(peos_new) - Average(delta_p_term, peosbar, 0); - - // compute p0_minus_peosbar = p0_new - peosbar - p0_minus_peosbar.copy(p0_new - peosbar); - - // put p0_new on cart - Put1dArrayOnCart(p0_new, p0_cart, false, false, bcs_f, 0); - - // set delta_p_term = peos_new - p0_new - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Subtract(delta_p_term[lev], p0_cart[lev], 0, 0, 1, - 0); - } - - } else { - // these should have no effect if dpdt_factor <= 0 - p0_minus_peosbar.setVal(0.); - for (int lev = 0; lev <= finest_level; ++lev) { - delta_p_term[lev].setVal(0.); - } - } - - if (evolve_base_state) { - if (split_projection) { - // compute Sbar = average(S_cc_nph) - Average(S_cc_nph, Sbar, 0); - - // compute Sbar = Sbar + delta_gamma1_termbar - if (use_delta_gamma1_term) { - Sbar += delta_gamma1_termbar; - } - - // compute w0, w0_force - is_predictor = false; - Makew0(w0_old, w0_force_dummy, Sbar, rho0_old, rho0_new, - p0_old, p0_new, gamma1bar_old, gamma1bar_new, - p0_minus_peosbar, dt, dtold, is_predictor); - - Put1dArrayOnCart(w0, w0_cart, true, true, bcs_u, 0, 1); -#if (AMREX_SPACEDIM == 3) - if (spherical) { - // put w0 on Cartesian edges - MakeW0mac(w0mac); - } -#endif - } - } - - // compute unprojected MAC velocities - is_predictor = false; - AdvancePremac(umac, w0mac_dummy, w0_force_cart_dummy); - - if (evolve_base_state && !split_projection) { - Sbar.copy(1.0 / (gamma1bar_nph * p0_nph) * (p0_nph - p0_old) / - dt); - } - - // compute RHS for MAC projection, beta0*(S_cc-Sbar) + beta0*delta_chi - MakeRHCCforMacProj(macrhs, rho0_new, S_cc_nph, Sbar, beta0_nph, - delta_gamma1_term, gamma1bar_new, p0_new, - delta_p_term, delta_chi, is_predictor); - - if (spherical && evolve_base_state && split_projection) { - // subtract w0mac from umac - Addw0(umac, w0mac, -1.); - } - - // wallclock time - Real start_total_macproj_corrector = ParallelDescriptor::second(); - - // MAC projection - // includes spherical option in C++ function - MacProj(umac, macphi, macrhs, beta0_nph, is_predictor); - - // wallclock time - Real end_total_macproj_corrector = - ParallelDescriptor::second() - start_total_macproj_corrector; - ParallelDescriptor::ReduceRealMax( - end_total_macproj_corrector, - ParallelDescriptor::IOProcessorNumber()); - - if (spherical && evolve_base_state && split_projection) { - // add w0mac back to umac - Addw0(umac, w0mac, 1.); - } - } // end sdc_couple_mac_velocity - - ////////////////////////////////////////////////////////////////////////////// - // STEP 4A -- compute advective flux divergences - ////////////////////////////////////////////////////////////////////////////// - - if (maestro_verbose >= 1) { - Print() << "<<< STEP 4 : Corrector loop (MISDC iter = " << misdc - << ") >>>" << std::endl; - Print() << "<<< STEP 4A : compute advective flux div (SDC iter = " - << misdc << ") >>>" << std::endl; - } - - // no need to advect the base state density - rho0_new.copy(rho0_old); - - if (maestro_verbose >= 1) { - Print() << " : density_advance >>>" << std::endl; - } - - // advect rhoX, rho, and tracers - DensityAdvanceSDC(2, sold, shat, sedge, sflux, scal_force, - etarhoflux_dummy, umac, w0mac, rho0_pred_edge_dummy); - - if (evolve_base_state) { - // correct the base state density by "averaging" - Average(shat, rho0_new, Rho); - ComputeCutoffCoords(rho0_new); - } - - // update grav_cell_new, rho0_nph, grav_cell_nph - if (evolve_base_state) { - MakeGravCell(grav_cell_new, rho0_new); - - rho0_nph.copy(0.5 * (rho0_old + rho0_new)); - - MakeGravCell(grav_cell_nph, rho0_nph); - } else { - rho0_nph.copy(rho0_old); - grav_cell_nph.copy(grav_cell_old); - } - - // base state pressure update - if (evolve_base_state) { - // set new p0 through HSE - p0_new.copy(p0_old); - - EnforceHSE(rho0_new, p0_new, grav_cell_new); - - p0_nph.copy(0.5 * (p0_old + p0_new)); - - // hold dp0/dt in psi for enthalpy advance - psi.copy((p0_new - p0_old) / dt); - } - - // enthalpy update - if (maestro_verbose >= 1) { - Print() << " : enthalpy_advance >>>" << std::endl; - } - - EnthalpyAdvanceSDC(2, sold, shat, sedge, sflux, scal_force, umac, w0mac, - diff_old); - - // base state enthalpy update - if (evolve_base_state) { - Average(shat, rhoh0_new, RhoH); - - // store (rhoh0_hat - rhoh0_old)/dt in delta_rhoh0 - delta_rhoh0.copy((rhoh0_new - rhoh0_old) / dt); - } - - // extract aofs = (shat - sold) / dt - intra - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::LinComb(aofs[lev], 1.0 / dt, shat[lev], 0, -1.0 / dt, - sold[lev], 0, 0, Nscal, 0); - MultiFab::Subtract(aofs[lev], intra[lev], 0, 0, Nscal, 0); - } - - ////////////////////////////////////////////////////////////////////////////// - // STEP 4B (optional) -- compute diffusive flux divergences - ////////////////////////////////////////////////////////////////////////////// - - if (maestro_verbose >= 1) { - Print() << "<<< STEP 4B : compute diffusive flux div (SDC iter = " - << misdc << ") >>>" << std::endl; - } - - // set diff_hterm to zero - for (int lev = 0; lev <= finest_level; ++lev) { - diff_hterm_new[lev].setVal(0.); - diff_hterm_hat[lev].setVal(0.); - } - - if (use_thermal_diffusion) { - // 1 = predictor, 2 = corrector - ThermalConductSDC(2, sold, shat, snew, hcoeff1, - Xkcoeff1, pcoeff1, hcoeff2, Xkcoeff2, pcoeff2); - - // compute diff_hat using shat, p0_new, and new coefficients from previous iteration - // note p0_new = p0_hat if evolve_base_state = T - MakeExplicitThermal(diff_hat, shat, Tcoeff2, hcoeff2, Xkcoeff2, - pcoeff2, p0_new, temp_diffusion_formulation); - - // compute only the h term in diff_hat and diff_new - MakeExplicitThermalHterm(diff_hterm_hat, shat, hcoeff2); - MakeExplicitThermalHterm(diff_hterm_new, snew, hcoeff2); - } - - ////////////////////////////////////////////////////////////////////////////// - // STEP 4C -- advance thermodynamic variables - ////////////////////////////////////////////////////////////////////////////// - - if (maestro_verbose >= 1) { - Print() << "<<< STEP 4C: advance thermo variables (MISDC iter = " - << misdc << ") >>>" << std::endl; - } - - // build sdc_source - for (int lev = 0; lev <= finest_level; ++lev) { - sdc_source[lev].setVal(0.); - MultiFab::Copy(sdc_source[lev], diff_old[lev], 0, RhoH, 1, 0); - MultiFab::Add(sdc_source[lev], diff_hat[lev], 0, RhoH, 1, 0); - MultiFab::Add(sdc_source[lev], diff_hterm_hat[lev], 0, RhoH, 1, 0); - MultiFab::Subtract(sdc_source[lev], diff_hterm_new[lev], 0, RhoH, 1, - 0); - sdc_source[lev].mult(0.5, RhoH, 1, 0); - MultiFab::Add(sdc_source[lev], aofs[lev], 0, 0, Nscal, 0); - } - - // wallclock time - Real start_total_react_corrector = ParallelDescriptor::second(); - - ReactSDC(sold, snew, rho_Hext, p0_new, dt, t_old, sdc_source); - - // wallclock time - Real end_total_react_corrector = - ParallelDescriptor::second() - start_total_react_corrector; - ParallelDescriptor::ReduceRealMax( - end_total_react_corrector, ParallelDescriptor::IOProcessorNumber()); - - // extract IR = [ (snew - sold)/dt - sdc_source ] - for (int lev = 0; lev <= finest_level; ++lev) { - intra[lev].setVal(0.); - // species source term - MultiFab::LinComb(intra[lev], 1.0 / dt, snew[lev], FirstSpec, - -1.0 / dt, sold[lev], FirstSpec, FirstSpec, - NumSpec, 0); - MultiFab::Subtract(intra[lev], sdc_source[lev], FirstSpec, - FirstSpec, NumSpec, 0); - // enthalpy source term - MultiFab::LinComb(intra[lev], 1.0 / dt, snew[lev], RhoH, -1.0 / dt, - sold[lev], RhoH, RhoH, 1, 0); - MultiFab::Subtract(intra[lev], sdc_source[lev], RhoH, RhoH, 1, 0); - } - - // massage the rhoh intra term into the proper form, depending on - // what we are predicting. Note: we do this before we deal with - // the species terms, since some enthalpy types need this default - // species intra. - - // create rhohalf -- a lot of forms need this. - FillPatch(0.5 * (t_old + t_new), rhohalf, sold, snew, Rho, 0, 1, Rho, - bcs_s); - - if (evolve_base_state) { - // update base state density and pressure - Average(snew, rho0_new, Rho); - ComputeCutoffCoords(rho0_new); - - if (use_etarho) { - // compute the new etarho - if (!spherical) { - MakeEtarho(etarhoflux_dummy); -#if AMREX_SPACEDIM == 3 - } else { - MakeEtarhoSphr(sold, snew, umac, w0mac_dummy); -#endif - } - } - - MakeGravCell(grav_cell_new, rho0_new); - - EnforceHSE(rho0_new, p0_new, grav_cell_new); - - // hold dp0/dt in psi for Make_S_cc - psi.copy((p0_new - p0_old) / dt); - - // also update base state enthalpy - Average(snew, rhoh0_new, RhoH); - - // compute intra_rhoh0 = (rhoh0_new - rhoh0_old)/dt - // - (rhoh0_hat - rhoh0_old)/dt - delta_rhoh0.copy((rhoh0_new - rhoh0_old) / dt - delta_rhoh0); - Put1dArrayOnCart(delta_rhoh0, intra_rhoh0, false, false, bcs_f, 0); - } - - // now update temperature - if (use_tfromp) { - TfromRhoP(snew, p0_new); - } else { - TfromRhoH(snew, p0_new); - } - - if (is_initIter) { - for (int lev = 0; lev <= finest_level; ++lev) { - intra[lev].setVal(0.); - intra_rhoh0[lev].setVal(0.); - } - } - - if (enthalpy_pred_type == predict_rhohprime) { - // intra is only different from predict_rhoh if rhoh0 is not constant - if (evolve_base_state) { - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Subtract(intra[lev], intra_rhoh0[lev], 0, RhoH, 1, - 0); - } - } - - } else if (enthalpy_pred_type == predict_h) { - // we want this in terms of h, not (rho h) - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Divide(intra[lev], rhohalf[lev], 0, RhoH, 1, 0); - } - - } else if ((enthalpy_pred_type == predict_T_then_rhohprime) || - (enthalpy_pred_type == predict_T_then_h)) { - // for predict_T_*, the intra force needs to be in the temp_comp - // slot, since temperature is what is predicted. - - // first make the thermodynamic coefficients at the half-time - MakeIntraCoeffs(sold, snew, cphalf, xihalf); - - // overwrite intra(temp_comp). We want to create - // I_T = (1 / (rho c_p)) [ (rhoh_new - rhoh_old)/dt - A_rhoh - - // sum_k xi_k ( (rhoX_new - rhoX_old)/dt - A_rhoX ) ] - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Copy(intra[lev], intra[lev], RhoH, Temp, 1, 0); - - for (int comp = 0; comp < NumSpec; ++comp) { - // multiple xi by intra and store in xi - MultiFab::Multiply(xihalf[lev], intra[lev], - FirstSpec + comp, comp, 1, 0); - - // subtract from intra temp - MultiFab::Subtract(intra[lev], xihalf[lev], comp, Temp, 1, - 0); - } - - MultiFab::Divide(intra[lev], rhohalf[lev], 0, Temp, 1, 0); - MultiFab::Divide(intra[lev], cphalf[lev], 0, Temp, 1, 0); - } - } - - // for some species_pred_types, we need to make intra in terms of - // X, NOT rhoX - if ((species_pred_type == predict_rhoprime_and_X) || - (species_pred_type == predict_rho_and_X)) { - for (int lev = 0; lev <= finest_level; ++lev) { - for (int comp = FirstSpec; comp < FirstSpec + NumSpec; ++comp) { - MultiFab::Divide(intra[lev], rhohalf[lev], 0, comp, 1, 0); - } - } - } - - // compute new-time coefficients and diffusion term - if (use_thermal_diffusion) { - MakeThermalCoeffs(snew, Tcoeff2, hcoeff2, Xkcoeff2, pcoeff2); - - MakeExplicitThermal(diff_new, snew, Tcoeff2, hcoeff2, Xkcoeff2, - pcoeff2, p0_new, temp_diffusion_formulation); - } - - if (evolve_base_state) { - // compute beta0 and gamma1bar - MakeGamma1bar(snew, gamma1bar_new, p0_new); - - MakeBeta0(beta0_new, rho0_new, p0_new, gamma1bar_new, - grav_cell_new); - } else { - // Just pass beta0 and gamma1bar through if not evolving base state - beta0_new.copy(beta0_old); - gamma1bar_new.copy(gamma1bar_old); - } - - gamma1bar_nph.copy(0.5 * (gamma1bar_old + gamma1bar_new)); - beta0_nph.copy(0.5 * (beta0_old + beta0_new)); - - } // end loop over misdc iterations - - ////////////////////////////////////////////////////////////////////////////// - // STEP 5 -- Advance velocity and dynamic pressure - ////////////////////////////////////////////////////////////////////////////// - - if (maestro_verbose >= 1) { - Print() << "<<< STEP 5 : Advance velocity and dynamic pressure >>>" - << std::endl; - } - - MakeReactionRates(rho_omegadot, rho_Hnuc, snew); - - Make_S_cc(S_cc_new, delta_gamma1_term, delta_gamma1, snew, uold, - rho_omegadot, rho_Hnuc, rho_Hext, diff_new, p0_new, gamma1bar_new, - delta_gamma1_termbar); - - if (evolve_base_state) { - if (split_projection) { - // compute Sbar = average(S_cc_new) - Average(S_cc_new, Sbar, 0); - - // compute Sbar = Sbar + delta_gamma1_termbar - if (use_delta_gamma1_term) { - Sbar += delta_gamma1_termbar; - } - - // compute w0, w0_force - is_predictor = false; - Makew0(w0_old, w0_force_dummy, Sbar, rho0_new, rho0_new, p0_new, - p0_new, gamma1bar_new, gamma1bar_new, p0_minus_peosbar, dt, - dtold, is_predictor); - - Put1dArrayOnCart(w0, w0_cart, true, true, bcs_u, 0, 1); - } - } - - // define dSdt = (S_cc_new - S_cc_old) / dt - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::LinComb(dSdt[lev], -1. / dt, S_cc_old[lev], 0, 1. / dt, - S_cc_new[lev], 0, 0, 1, 0); - } - - // Define rho at half time using the new rho from Step 4 - FillPatch(0.5 * (t_old + t_new), rhohalf, sold, snew, Rho, 0, 1, Rho, - bcs_s); - - VelocityAdvance(rhohalf, umac, w0mac_dummy, w0_force_cart_dummy, rho0_nph, - grav_cell_nph, sponge); - - if (evolve_base_state && is_initIter) { - // throw away w0 by setting w0 = w0_old - w0.copy(w0_old); - } - - if (spherical && evolve_base_state && split_projection) { - // subtract w0 from uold and unew for nodal projection - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Subtract(uold[lev], w0_cart[lev], 0, 0, AMREX_SPACEDIM, - 0); - MultiFab::Subtract(unew[lev], w0_cart[lev], 0, 0, AMREX_SPACEDIM, - 0); - } - } - - if (evolve_base_state && !split_projection) { - Sbar.copy((p0_new - p0_old) / (dt * gamma1bar_new * p0_new)); - } - - int proj_type; - - // Project the new velocity field - if (is_initIter) { - proj_type = pressure_iters_comp; - - // rhcc_for_nodalproj needs to contain - // (beta0^nph S^1 - beta0^n S^0 ) / dt - - Vector rhcc_for_nodalproj_old(finest_level + 1); - for (int lev = 0; lev <= finest_level; ++lev) { - rhcc_for_nodalproj_old[lev].define(grids[lev], dmap[lev], 1, 1); - MultiFab::Copy(rhcc_for_nodalproj_old[lev], rhcc_for_nodalproj[lev], - 0, 0, 1, 1); - } - - MakeRHCCforNodalProj(rhcc_for_nodalproj, S_cc_new, Sbar, beta0_nph, - delta_gamma1_term); - - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Subtract(rhcc_for_nodalproj[lev], - rhcc_for_nodalproj_old[lev], 0, 0, 1, 1); - rhcc_for_nodalproj[lev].mult(1. / dt, 0, 1, 1); - } - - } else { - proj_type = regular_timestep_comp; - - MakeRHCCforNodalProj(rhcc_for_nodalproj, S_cc_new, Sbar, beta0_nph, - delta_gamma1_term); - - // compute delta_p_term = peos_new - p0_new (for RHS of projection) - if (dpdt_factor > 0.) { - // peos now holds "peos_new", the thermodynamic p computed from snew(rho,h,X) - PfromRhoH(snew, snew, delta_p_term); - - // put p0_new on cart - Put1dArrayOnCart(p0_new, p0_cart, false, false, bcs_f, 0); - - // compute delta_p_term = peos_new - p0_new - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Subtract(delta_p_term[lev], p0_cart[lev], 0, 0, 1, 0); - } - - CorrectRHCCforNodalProj(rhcc_for_nodalproj, rho0_new, beta0_nph, - gamma1bar_new, p0_new, delta_p_term); - } - } - - // wallclock time - const Real start_total_nodalproj = ParallelDescriptor::second(); - - // call nodal projection - NodalProj(proj_type, rhcc_for_nodalproj, 0, false); - - // wallclock time - Real end_total_nodalproj = - ParallelDescriptor::second() - start_total_nodalproj; - ParallelDescriptor::ReduceRealMax(end_total_nodalproj, - ParallelDescriptor::IOProcessorNumber()); - - if (spherical && evolve_base_state && split_projection) { - // add w0 back to unew - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Add(unew[lev], w0_cart[lev], 0, 0, AMREX_SPACEDIM, 0); - } - AverageDown(unew, 0, AMREX_SPACEDIM); - FillPatch(t_new, unew, unew, unew, 0, 0, AMREX_SPACEDIM, 0, bcs_u, 1); - } - - beta0_nm1.copy(0.5 * (beta0_old + beta0_new)); - - if (!is_initIter) { - if (!fix_base_state) { - // compute tempbar by "averaging" - Average(snew, tempbar, Temp); - } - } - - Print() << "\nTimestep " << istep << " ends with TIME = " << t_new - << " DT = " << dt << std::endl; - - // print wallclock time - if (maestro_verbose > 0) { - Print() << "Time to solve mac proj : " << end_total_macproj << '\n'; - Print() << "Time to solve nodal proj : " << end_total_nodalproj << '\n'; - Print() << "Time to solve reactions : " << end_total_react << '\n'; - } -} diff --git a/Source/MaestroBurner.cpp b/Source/MaestroBurner.cpp index 31063e894..e085c079e 100644 --- a/Source/MaestroBurner.cpp +++ b/Source/MaestroBurner.cpp @@ -3,7 +3,6 @@ using namespace amrex; -#ifndef SDC void Maestro::Burner(const Vector& s_in, Vector& s_out, const Vector& rho_Hext, Vector& rho_omegadot, Vector& rho_Hnuc, @@ -235,181 +234,3 @@ void Maestro::Burner(const Vector& s_in, Vector& s_out, } } } - -#else -// SDC burner -void Maestro::Burner(const Vector& s_in, Vector& s_out, - const BaseState& p0, const Real dt_in, - const Real time_in, const Vector& source) { - // timer for profiling - BL_PROFILE_VAR("Maestro::BurnerSDC()", BurnerSDC); - - // Put tempbar_init on cart - Vector p0_cart(finest_level + 1); - const auto ispec_threshold = network_spec_index(burner_threshold_species); - - for (int lev = 0; lev <= finest_level; ++lev) { - p0_cart[lev].define(grids[lev], dmap[lev], 1, 0); - p0_cart[lev].setVal(0.); - } - if (spherical) { - Put1dArrayOnCart(p0, p0_cart, false, false, bcs_f, 0); - } - - for (int lev = 0; lev <= finest_level; ++lev) { - // create mask assuming refinement ratio = 2 - int finelev = lev + 1; - if (lev == finest_level) finelev = finest_level; - - const BoxArray& fba = s_in[finelev].boxArray(); - const iMultiFab& mask = makeFineMask(s_in[lev], fba, IntVect(2)); - - // loop over boxes (make sure mfi takes a cell-centered multifab as an argument) -#ifdef _OPENMP -#pragma omp parallel -#endif - for (MFIter mfi(s_in[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) { - // Get the index space of the valid region - const Box& tileBox = mfi.tilebox(); - - const bool use_mask = !(lev == finest_level); - - const Array4 s_in_arr = s_in[lev].array(mfi); - const Array4 s_out_arr = s_out[lev].array(mfi); - const Array4 p0_cart_arr = p0_cart[lev].array(mfi); - const Array4 source_arr = source[lev].array(mfi); - const Array4 mask_arr = mask.array(mfi); - - const auto p0_arr = p0.const_array(); - - ParallelFor(tileBox, [=] AMREX_GPU_DEVICE(int i, int j, int k) { - if (use_mask && mask_arr(i, j, k)) - return; // cell is covered by finer cells - - auto r = (AMREX_SPACEDIM == 2) ? j : k; - - Real sdc_rhoX[NumSpec]; - for (int n = 0; n < NumSpec; ++n) { - sdc_rhoX[n] = source_arr(i, j, k, FirstSpec + n); - } - auto sdc_rhoh = source_arr(i, j, k, RhoH); - auto sdc_p0 = spherical ? p0_cart_arr(i, j, k) : p0_arr(lev, r); - - auto rho_in = s_in_arr(i, j, k, Rho); - Real rhoX_in[NumSpec]; - for (int n = 0; n < NumSpec; ++n) { - rhoX_in[n] = s_in_arr(i, j, k, FirstSpec + n); - } -#if NAUX_NET > 0 - Real aux_in[NumAux]; - for (int n = 0; n < NumAux; ++n) { - aux_in[n] = s_in_arr(i, j, k, FirstAux + n); - } -#endif - auto rhoh_in = s_in_arr(i, j, k, RhoH); - - Real x_test = (ispec_threshold > 0) - ? rhoX_in[ispec_threshold] / rho_in - : 0.0; - - burn_t state_in; - burn_t state_out; - - Real rhoX_out[NumSpec]; -#if NAUX_NET > 0 - Real aux_out[NumAux]; -#endif - Real rho_out = 0.0; - Real rhoh_out = 0.0; - - // if the threshold species is not in the network, then we burn - // normally. if it is in the network, make sure the mass - // fraction is above the cutoff. - if ((rho_in > burning_cutoff_density_lo && - rho_in < burning_cutoff_density_hi) && - (ispec_threshold < 0 || - (ispec_threshold > 0 && - x_test > burner_threshold_cutoff))) { - state_in.p0 = sdc_p0; - state_in.rho = rho_in; - for (int n = 0; n < NumSpec; ++n) { - state_in.y[n] = rhoX_in[n]; - } - state_in.y[SENTH] = rhoh_in; - for (int n = 0; n < NumSpec; ++n) { - state_in.ydot_a[n] = sdc_rhoX[n]; - } - state_in.ydot_a[SENTH] = sdc_rhoh; -#if NAUX_NET > 0 - for (int n = 0; n < NumAux; ++n) { - state_in.aux[n] = rhoaux_in[n] / rho_in; - } -#endif - state_in.success = true; - - // initialize state_out the same - state_out.p0 = sdc_p0; - state_out.rho = rho_in; - for (int n = 0; n < NumSpec; ++n) { - state_out.y[n] = rhoX_in[n]; - } - state_out.y[NumSpec] = rhoh_in; - for (int n = 0; n < NumSpec; ++n) { - state_out.ydot_a[n] = sdc_rhoX[n]; - } - state_out.ydot_a[NumSpec] = sdc_rhoh; -#if NAUX_NET > 0 - for (int n = 0; n < NumAux; ++n) { - state_out.aux[n] = rhoaux_in[n] / rho_in; - } -#endif - state_out.success = true; - - integrator(state_out, dt_in); - - for (int n = 0; n < NumSpec; ++n) { - rho_out += state_out.y[n]; - rhoX_out[n] = state_out.y[n]; - } - rhoh_out = state_out.y[NumSpec]; -#if NAUX_NET > 0 - for (int n = 0; n < NumAux; ++n) { - rhoaux_out[n] = rho_out * state_out.aux[n]; - } -#endif - } else { - rho_out = rho_in; - for (int n = 0; n < NumSpec; ++n) { - rho_out += sdc_rhoX[n] * dt_in; - rhoX_out[n] = rhoX_in[n] + sdc_rhoX[n] * dt_in; - } - rhoh_out = rhoh_in + sdc_rhoh * dt_in; -#if NAUX_NET > 0 - for (int n = 0; n < NumAux; ++n) { - rhoaux_out[n] = rho_out / rho_in * rhoaux_in[n]; - } -#endif - } - - // update the density - s_out_arr(i, j, k, Rho) = rho_out; - - // update the species - for (int n = 0; n < NumSpec; ++n) { - s_out_arr(i, j, k, FirstSpec + n) = rhoX_out[n]; - } -#if NAUX_NET > 0 - for (int n = 0; n < NumAux; ++n) { - s_out_arr(i, j, k, FirstAux + n) = rhoaux_out[n]; - } -#endif - - // update the enthalpy -- include the change due to external heating - s_out_arr(i, j, k, RhoH) = rhoh_out; - - // pass the tracers through (currently not implemented) - }); - } - } -} -#endif diff --git a/Source/MaestroCheckpoint.cpp b/Source/MaestroCheckpoint.cpp index 366c0c875..e6381e8d0 100644 --- a/Source/MaestroCheckpoint.cpp +++ b/Source/MaestroCheckpoint.cpp @@ -111,10 +111,6 @@ void Maestro::WriteCheckPoint(int step) { VisMF::Write(S_cc_new[lev], amrex::MultiFabFileFullPrefix(lev, checkpointname, "Level_", "S_cc_new")); -#ifdef SDC - VisMF::Write(intra[lev], amrex::MultiFabFileFullPrefix( - lev, checkpointname, "Level_", "intra")); -#endif } // Restore the previous FAB format. @@ -267,10 +263,6 @@ int Maestro::ReadCheckPoint() { amrex::MultiFabFileFullPrefix(lev, restart_file, "Level_", "S_cc_new")); -#ifdef SDC - VisMF::Read(intra[lev], amrex::MultiFabFileFullPrefix( - lev, restart_file, "Level_", "intra")); -#endif } // get the elapsed CPU time to now; diff --git a/Source/MaestroDensityAdvance.cpp b/Source/MaestroDensityAdvance.cpp index b7f3a41db..59502a297 100644 --- a/Source/MaestroDensityAdvance.cpp +++ b/Source/MaestroDensityAdvance.cpp @@ -235,252 +235,3 @@ void Maestro::DensityAdvance( UpdateScal(scalold, scalnew, sflux, scal_force, FirstSpec, NumSpec, p0_new_cart); } - -// Density advance for SDC using intra(global var) -void Maestro::DensityAdvanceSDC( - int which_step, Vector& scalold, Vector& scalnew, - Vector >& sedge, - Vector >& sflux, - Vector& scal_force, Vector& etarhoflux, - Vector >& umac, - const Vector >& w0mac, - const BaseState& rho0_predicted_edge) { - // timer for profiling - BL_PROFILE_VAR("Maestro::DensityAdvanceSDC()", DensityAdvanceSDC); - - BaseState rho0_edge_old(base_geom.max_radial_level + 1, - base_geom.nr_fine + 1); - BaseState rho0_edge_new(base_geom.max_radial_level + 1, - base_geom.nr_fine + 1); - - if (!spherical) { - // create edge-centered base state quantities. - // Note: rho0_edge_{old,new} - // contains edge-centered quantities created via spatial interpolation. - // This is to be contrasted to rho0_predicted_edge which is the half-time - // edge state created in advect_base. - CelltoEdge(rho0_old, rho0_edge_old); - CelltoEdge(rho0_new, rho0_edge_new); - } - - ////////////////////////////////// - // Create source terms at time n - ////////////////////////////////// - - // source terms for X and for tracers include reaction forcing terms - for (int lev = 0; lev <= finest_level; ++lev) { - scal_force[lev].setVal(0.); - MultiFab::Add(scal_force[lev], intra[lev], FirstSpec, FirstSpec, - NumSpec, 0); - } - - if (finest_level == 0) { - // fill periodic ghost cells - for (int lev = 0; lev <= finest_level; ++lev) { - scal_force[lev].FillBoundary(geom[lev].periodicity()); - } - } - // fill ghost cells behind physical boundaries - // !!!!!! uncertain about this - FillPatch(t_old, scal_force, scal_force, scal_force, FirstSpec, FirstSpec, - NumSpec, FirstSpec, bcs_f); - - Vector rho0_old_cart(finest_level + 1); - for (int lev = 0; lev <= finest_level; ++lev) { - rho0_old_cart[lev].define(grids[lev], dmap[lev], 1, 1); - } - - Put1dArrayOnCart(rho0_old, rho0_old_cart, false, false, bcs_s, Rho); - - ///////////////////////////////////////////////////////////////// - // Subtract w0 from MAC velocities (MAC velocities has w0 already). - ///////////////////////////////////////////////////////////////// - - Addw0(umac, w0mac, -1.); - - ///////////////////////////////////////////////////////////////// - // Compute source terms - ///////////////////////////////////////////////////////////////// - - // ** density source term ** - - // Make source term for rho or rho' - if (species_pred_type == predict_rhoprime_and_X) { - // rho' source term - // this is needed for pred_rhoprime_and_X - ModifyScalForce(scal_force, scalold, umac, rho0_edge_old, rho0_old_cart, - Rho, bcs_s, 0); - } else if (species_pred_type == predict_rho_and_X) { - // rho source term - ModifyScalForce(scal_force, scalold, umac, rho0_edge_old, rho0_old_cart, - Rho, bcs_s, 1); - } - - // ** species source term ** - - // for species_pred_types predict_rhoprime_and_X and - // predict_rho_and_X, there is no force for X. - - // for predict_rhoX, we are predicting (rho X) - // as a conservative equation, and there is no force. - - ///////////////////////////////////////////////////////////////// - // Add w0 to MAC velocities (trans velocities already have w0). - ///////////////////////////////////////////////////////////////// - - Addw0(umac, w0mac, 1.); - - ///////////////////////////////////////////////////////////////// - // Create the edge states of (rho X)' or X and rho' - ///////////////////////////////////////////////////////////////// - - if ((species_pred_type == predict_rhoprime_and_X) || - (species_pred_type == predict_rho_and_X)) { - // we are predicting X to the edges, so convert the scalar - // data to those quantities - - // convert (rho X) --> X in scalold - ConvertRhoXToX(scalold, true); - } - - if (species_pred_type == predict_rhoprime_and_X) { - // convert rho -> rho' in scalold - // . this is needed for predict_rhoprime_and_X - PutInPertForm(scalold, rho0_old, Rho, 0, bcs_f, true); - } - - // predict species at the edges -- note, either X or (rho X) will be - // predicted here, depending on species_pred_type - - const auto is_vel = false; // false - if (species_pred_type == predict_rhoprime_and_X || - species_pred_type == predict_rho_and_X) { - // we are predicting X to the edges, using the advective form of - // the prediction - MakeEdgeScal(scalold, sedge, umac, scal_force, is_vel, bcs_s, Nscal, - FirstSpec, FirstSpec, NumSpec, false); - - } else if (species_pred_type == predict_rhoX) { - MakeEdgeScal(scalold, sedge, umac, scal_force, is_vel, bcs_s, Nscal, - FirstSpec, FirstSpec, NumSpec, true); - } - - // predict rho or rho' at the edges (depending on species_pred_type) - if (species_pred_type == predict_rhoprime_and_X || - species_pred_type == predict_rho_and_X) { - MakeEdgeScal(scalold, sedge, umac, scal_force, is_vel, bcs_s, Nscal, - Rho, Rho, 1, false); - - } else if (species_pred_type == predict_rhoX) { - for (int lev = 0; lev <= finest_level; ++lev) { - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - MultiFab::Copy(sedge[lev][idim], sedge[lev][idim], FirstSpec, - Rho, 1, 0); - for (int ispec = 1; ispec < NumSpec; ++ispec) { - MultiFab::Add(sedge[lev][idim], sedge[lev][idim], - FirstSpec + ispec, Rho, 1, 0); - } - } - } - } - - if (species_pred_type == predict_rhoprime_and_X) { - // convert rho' -> rho in scalold - PutInPertForm(scalold, rho0_old, Rho, Rho, bcs_s, false); - } - - if ((species_pred_type == predict_rhoprime_and_X) || - (species_pred_type == predict_rho_and_X)) { - // convert X --> (rho X) in scalold - ConvertRhoXToX(scalold, false); - } - - ///////////////////////////////////////////////////////////////// - // Compute fluxes - ///////////////////////////////////////////////////////////////// - - if (which_step == 1) { - Vector > rho0mac_old(finest_level + - 1); - -#if (AMREX_SPACEDIM == 3) - if (spherical) { - for (int lev = 0; lev <= finest_level; ++lev) { - AMREX_D_TERM( - rho0mac_old[lev][0].define( - convert(grids[lev], nodal_flag_x), dmap[lev], 1, 1); - , rho0mac_old[lev][1].define( - convert(grids[lev], nodal_flag_y), dmap[lev], 1, 1); - , rho0mac_old[lev][2].define( - convert(grids[lev], nodal_flag_z), dmap[lev], 1, 1);); - } - MakeS0mac(rho0_old, rho0mac_old); - } -#endif - - // compute species fluxes - MakeRhoXFlux(scalold, sflux, etarhoflux, sedge, umac, rho0_old, - rho0_edge_old, rho0mac_old, rho0_old, rho0_edge_old, - rho0mac_old, rho0_predicted_edge, FirstSpec, NumSpec); - - } else if (which_step == 2) { - Vector > rho0mac_old(finest_level + - 1); - Vector > rho0mac_new(finest_level + - 1); - -#if (AMREX_SPACEDIM == 3) - if (spherical) { - for (int lev = 0; lev <= finest_level; ++lev) { - AMREX_D_TERM( - rho0mac_old[lev][0].define( - convert(grids[lev], nodal_flag_x), dmap[lev], 1, 1); - , rho0mac_old[lev][1].define( - convert(grids[lev], nodal_flag_y), dmap[lev], 1, 1); - , rho0mac_old[lev][2].define( - convert(grids[lev], nodal_flag_z), dmap[lev], 1, 1);); - - AMREX_D_TERM( - rho0mac_new[lev][0].define( - convert(grids[lev], nodal_flag_x), dmap[lev], 1, 1); - , rho0mac_new[lev][1].define( - convert(grids[lev], nodal_flag_y), dmap[lev], 1, 1); - , rho0mac_new[lev][2].define( - convert(grids[lev], nodal_flag_z), dmap[lev], 1, 1);); - } - MakeS0mac(rho0_old, rho0mac_old); - MakeS0mac(rho0_new, rho0mac_new); - } -#endif - - // compute species fluxes - MakeRhoXFlux(scalold, sflux, etarhoflux, sedge, umac, rho0_old, - rho0_edge_old, rho0mac_old, rho0_new, rho0_edge_new, - rho0mac_new, rho0_predicted_edge, FirstSpec, NumSpec); - } - - //************************************************************************** - // 1) Set force for (rho X)_i at time n+1/2 = 0. - // 2) Update (rho X)_i with conservative differencing. - // 3) Define density as the sum of the (rho X)_i - // 4) Update tracer with conservative differencing as well. - //************************************************************************** - - // reaction forcing terms - for (int lev = 0; lev <= finest_level; ++lev) { - scal_force[lev].setVal(0.); - MultiFab::Add(scal_force[lev], intra[lev], FirstSpec, FirstSpec, - NumSpec, 0); - } - - Vector p0_new_cart(finest_level + 1); - for (int lev = 0; lev <= finest_level; ++lev) { - p0_new_cart[lev].define(grids[lev], dmap[lev], 1, 1); - } - - Put1dArrayOnCart(p0_new, p0_new_cart, false, false, bcs_f, 0); - - // p0 only used in rhoh update so it's an optional parameter - UpdateScal(scalold, scalnew, sflux, scal_force, FirstSpec, NumSpec, - p0_new_cart); -} diff --git a/Source/MaestroDiag.cpp b/Source/MaestroDiag.cpp index a857d061e..b93d9f635 100644 --- a/Source/MaestroDiag.cpp +++ b/Source/MaestroDiag.cpp @@ -30,7 +30,6 @@ void Maestro::DiagFile(const int step, const Real t_in, Vector rho_Hext(finest_level + 1); Vector rho_omegadot(finest_level + 1); Vector rho_Hnuc(finest_level + 1); - Vector sdc_source(finest_level + 1); #if (AMREX_SPACEDIM == 3) if (spherical) { @@ -61,12 +60,8 @@ void Maestro::DiagFile(const int step, const Real t_in, rho_Hext[lev].define(grids[lev], dmap[lev], 1, 0); rho_omegadot[lev].define(grids[lev], dmap[lev], NumSpec, 0); rho_Hnuc[lev].define(grids[lev], dmap[lev], 1, 0); - sdc_source[lev].define(grids[lev], dmap[lev], Nscal, 0); - - sdc_source[lev].setVal(0.0); } -#ifndef SDC if (dt < small_dt) { React(s_in, stemp, rho_Hext, rho_omegadot, rho_Hnuc, p0_in, small_dt, t_in); @@ -74,15 +69,6 @@ void Maestro::DiagFile(const int step, const Real t_in, React(s_in, stemp, rho_Hext, rho_omegadot, rho_Hnuc, p0_in, dt * 0.5, t_in); } -#else - if (dt < small_dt) { - ReactSDC(s_in, stemp, rho_Hext, p0_in, small_dt, t_in, sdc_source); - } else { - ReactSDC(s_in, stemp, rho_Hext, p0_in, dt * 0.5, t_in, sdc_source); - } - - MakeReactionRates(rho_omegadot, rho_Hnuc, s_in); -#endif // initialize diagnosis variables // diag_temp.out diff --git a/Source/MaestroEnthalpyAdvance.cpp b/Source/MaestroEnthalpyAdvance.cpp index 9a92b44b6..5dd33e15f 100644 --- a/Source/MaestroEnthalpyAdvance.cpp +++ b/Source/MaestroEnthalpyAdvance.cpp @@ -350,373 +350,3 @@ void Maestro::EnthalpyAdvance( UpdateScal(scalold, scalnew, sflux, scal_force, RhoH, 1, p0_new_cart); } - -// Enthalpy advance for SDC using intra(global var) -void Maestro::EnthalpyAdvanceSDC( - int which_step, Vector& scalold, Vector& scalnew, - Vector >& sedge, - Vector >& sflux, - Vector& scal_force, - Vector >& umac, - const Vector >& w0mac, - const Vector& thermal) { - // timer for profiling - BL_PROFILE_VAR("Maestro::EnthalpyAdvanceSDC()", EnthalpyAdvanceSDC); - - // Create cell-centered base state quantity - BaseState h0_old(base_geom.max_radial_level + 1, base_geom.nr_fine); - BaseState h0_new(base_geom.max_radial_level + 1, base_geom.nr_fine); - - // Create edge-centered base state quantities. - // Note: rho0_edge_{old,new} and rhoh0_edge_{old,new} - // contain edge-centered quantities created via spatial interpolation. - BaseState rho0_edge_old(base_geom.max_radial_level + 1, - base_geom.nr_fine + 1); - BaseState rho0_edge_new(base_geom.max_radial_level + 1, - base_geom.nr_fine + 1); - BaseState rhoh0_edge_old(base_geom.max_radial_level + 1, - base_geom.nr_fine + 1); - BaseState rhoh0_edge_new(base_geom.max_radial_level + 1, - base_geom.nr_fine + 1); - - if (!spherical) { - CelltoEdge(rho0_old, rho0_edge_old); - CelltoEdge(rho0_new, rho0_edge_new); - CelltoEdge(rhoh0_old, rhoh0_edge_old); - CelltoEdge(rhoh0_new, rhoh0_edge_new); - } - - if (enthalpy_pred_type == predict_h || - enthalpy_pred_type == predict_hprime) { - // convert (rho h) -> h - ConvertRhoHToH(scalold, true); - } - - ////////////////////////////////// - // Create scalar source term at time n - ////////////////////////////////// - - for (int lev = 0; lev <= finest_level; ++lev) { - scal_force[lev].setVal(0.); - } - - Vector rhoh0_old_cart(finest_level + 1); - for (int lev = 0; lev <= finest_level; ++lev) { - rhoh0_old_cart[lev].define(grids[lev], dmap[lev], 1, 1); - // needed to avoid NaNs in filling corner ghost cells with 2 physical boundaries - rhoh0_old_cart[lev].setVal(0.); - } - - ///////////////////////////////////////////////////////////////// - // Subtract w0 from MAC velocities (MAC velocities has w0 already). - ///////////////////////////////////////////////////////////////// - - Addw0(umac, w0mac, -1.); - - ///////////////////////////////////////////////////////////////// - // Compute forcing terms - ///////////////////////////////////////////////////////////////// - - if (enthalpy_pred_type == predict_rhohprime) { - // make force for (rho h)' - MakeRhoHForce(scal_force, 1, thermal, umac, 1, 1); - - Put1dArrayOnCart(rhoh0_old, rhoh0_old_cart, false, false, bcs_s, RhoH); - - ModifyScalForce(scal_force, scalold, umac, rhoh0_edge_old, - rhoh0_old_cart, RhoH, bcs_s, 0); - } else if (enthalpy_pred_type == predict_h || - enthalpy_pred_type == predict_rhoh) { - // make force for (rho h) - MakeRhoHForce(scal_force, 1, thermal, umac, 1, 1); - - // make force for h by calling mkrhohforce then dividing by rho - if (enthalpy_pred_type == predict_h) { - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Divide(scal_force[lev], scalold[lev], Rho, RhoH, 1, - 1); - } - } - } else if (enthalpy_pred_type == predict_hprime) { - // first compute h0_old - // make force for hprime - Abort( - "MaestroEnthalpyAdvance does not support enthalpy_pred_type == " - "predict_hprime"); - } else if (enthalpy_pred_type == predict_T_then_rhohprime || - enthalpy_pred_type == predict_T_then_h || - enthalpy_pred_type == predict_Tprime_then_h) { - // make force for temperature - MakeTempForce(scal_force, scalold, thermal, umac); - } - - // source terms for X and for tracers include reaction forcing terms - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Add(scal_force[lev], intra[lev], RhoH, RhoH, 1, 0); - } - - if (finest_level == 0) { - // fill periodic ghost cells - for (int lev = 0; lev <= finest_level; ++lev) { - scal_force[lev].FillBoundary(geom[lev].periodicity()); - } - } - // fill ghost cells behind physical boundaries - // !!!!!! uncertain about this - FillPatch(t_old, scal_force, scal_force, scal_force, RhoH, RhoH, 1, RhoH, - bcs_f); - - ////////////////////////////////// - // Add w0 to MAC velocities (trans velocities already have w0). - ////////////////////////////////// - - Addw0(umac, w0mac, 1.); - - ////////////////////////////////// - // Create the edge states of (rho h)' or h or T - ////////////////////////////////// - - if (enthalpy_pred_type == predict_rhohprime) { - // convert (rho h) -> (rho h)' - PutInPertForm(scalold, rhoh0_old, RhoH, 0, bcs_f, true); - } - - if (enthalpy_pred_type == predict_hprime) { - // convert h -> h' - Abort("MaestroEnthalpyAdvance predict_hprime"); - } - - if (enthalpy_pred_type == predict_Tprime_then_h) { - // convert T -> T' - PutInPertForm(scalold, tempbar, Temp, 0, bcs_f, true); - } - - // predict either T, h, or (rho h)' at the edges - int pred_comp = 0; - if (enthalpy_pred_type == predict_T_then_rhohprime || - enthalpy_pred_type == predict_T_then_h || - enthalpy_pred_type == predict_Tprime_then_h) { - pred_comp = Temp; - } else { - pred_comp = RhoH; - } - - if (enthalpy_pred_type == predict_rhoh) { - // use the conservative form of the prediction - MakeEdgeScal(scalold, sedge, umac, scal_force, false, bcs_s, Nscal, - pred_comp, pred_comp, 1, true); - } else { - // use the advective form of the prediction - MakeEdgeScal(scalold, sedge, umac, scal_force, false, bcs_s, Nscal, - pred_comp, pred_comp, 1, true); - } - - if (enthalpy_pred_type == predict_rhohprime) { - // convert (rho h)' -> (rho h) - PutInPertForm(scalold, rhoh0_old, RhoH, RhoH, bcs_s, false); - } - - if (enthalpy_pred_type == predict_hprime) { - // convert h' -> h - Abort("MaestroEnthalpyAdvance predict_hprime"); - } - - if (enthalpy_pred_type == predict_Tprime_then_h) { - // convert T' -> T - PutInPertForm(scalold, tempbar, Temp, Temp, bcs_s, false); - } - - if (enthalpy_pred_type == predict_h || - enthalpy_pred_type == predict_hprime) { - // convert (rho h) -> h - ConvertRhoHToH(scalold, false); - } - - // Compute enthalpy edge states if we were predicting temperature. This - // needs to be done after the state was returned to the full state. - if ((enthalpy_pred_type == predict_T_then_rhohprime) || - (enthalpy_pred_type == predict_T_then_h) || - (enthalpy_pred_type == predict_Tprime_then_h)) { - HfromRhoTedge(sedge, rho0_edge_old, rhoh0_edge_old, rho0_edge_new, - rhoh0_edge_new); - } - - ////////////////////////////////// - // Compute fluxes - ////////////////////////////////// - - // for which_step .eq. 1, we pass in only the old base state quantities - // for which_step .eq. 2, we pass in the old and new for averaging within mkflux - if (which_step == 1) { - Vector > rho0mac_old(finest_level + - 1); - Vector > rhoh0mac_old( - finest_level + 1); - Vector > h0mac_old(finest_level + - 1); - -#if (AMREX_SPACEDIM == 3) - if (spherical) { - h0_old.copy(rhoh0_old / rho0_old); - - for (int lev = 0; lev <= finest_level; ++lev) { - AMREX_D_TERM( - rho0mac_old[lev][0].define( - convert(grids[lev], nodal_flag_x), dmap[lev], 1, 1); - , rho0mac_old[lev][1].define( - convert(grids[lev], nodal_flag_y), dmap[lev], 1, 1); - , rho0mac_old[lev][2].define( - convert(grids[lev], nodal_flag_z), dmap[lev], 1, 1);); - - AMREX_D_TERM( - rhoh0mac_old[lev][0].define( - convert(grids[lev], nodal_flag_x), dmap[lev], 1, 1); - , rhoh0mac_old[lev][1].define( - convert(grids[lev], nodal_flag_y), dmap[lev], 1, 1); - , rhoh0mac_old[lev][2].define( - convert(grids[lev], nodal_flag_z), dmap[lev], 1, 1);); - - AMREX_D_TERM( - h0mac_old[lev][0].define(convert(grids[lev], nodal_flag_x), - dmap[lev], 1, 1); - , h0mac_old[lev][1].define( - convert(grids[lev], nodal_flag_y), dmap[lev], 1, 1); - , h0mac_old[lev][2].define( - convert(grids[lev], nodal_flag_z), dmap[lev], 1, 1);); - } - MakeS0mac(rho0_old, rho0mac_old); - MakeS0mac(rhoh0_old, rhoh0mac_old); - MakeS0mac(h0_old, h0mac_old); - } -#endif - - // compute enthalpy fluxes - MakeRhoHFlux(scalold, sflux, sedge, umac, w0mac, rho0_old, - rho0_edge_old, rho0mac_old, rho0_old, rho0_edge_old, - rho0mac_old, rhoh0_old, rhoh0_edge_old, rhoh0mac_old, - rhoh0_old, rhoh0_edge_old, rhoh0mac_old, h0mac_old, - h0mac_old); - } else if (which_step == 2) { - Vector > rho0mac_old(finest_level + - 1); - Vector > rhoh0mac_old( - finest_level + 1); - Vector > h0mac_old(finest_level + - 1); - Vector > rho0mac_new(finest_level + - 1); - Vector > rhoh0mac_new( - finest_level + 1); - Vector > h0mac_new(finest_level + - 1); - -#if (AMREX_SPACEDIM == 3) - if (spherical) { - h0_old.copy(rhoh0_old / rho0_old); - h0_new.copy(rhoh0_new / rho0_new); - - for (int lev = 0; lev <= finest_level; ++lev) { - AMREX_D_TERM( - rho0mac_old[lev][0].define( - convert(grids[lev], nodal_flag_x), dmap[lev], 1, 1); - , rho0mac_old[lev][1].define( - convert(grids[lev], nodal_flag_y), dmap[lev], 1, 1); - , rho0mac_old[lev][2].define( - convert(grids[lev], nodal_flag_z), dmap[lev], 1, 1);); - - AMREX_D_TERM( - rhoh0mac_old[lev][0].define( - convert(grids[lev], nodal_flag_x), dmap[lev], 1, 1); - , rhoh0mac_old[lev][1].define( - convert(grids[lev], nodal_flag_y), dmap[lev], 1, 1); - , rhoh0mac_old[lev][2].define( - convert(grids[lev], nodal_flag_z), dmap[lev], 1, 1);); - - AMREX_D_TERM( - h0mac_old[lev][0].define(convert(grids[lev], nodal_flag_x), - dmap[lev], 1, 1); - , h0mac_old[lev][1].define( - convert(grids[lev], nodal_flag_y), dmap[lev], 1, 1); - , h0mac_old[lev][2].define( - convert(grids[lev], nodal_flag_z), dmap[lev], 1, 1);); - - AMREX_D_TERM( - rho0mac_new[lev][0].define( - convert(grids[lev], nodal_flag_x), dmap[lev], 1, 1); - , rho0mac_new[lev][1].define( - convert(grids[lev], nodal_flag_y), dmap[lev], 1, 1); - , rho0mac_new[lev][2].define( - convert(grids[lev], nodal_flag_z), dmap[lev], 1, 1);); - - AMREX_D_TERM( - rhoh0mac_new[lev][0].define( - convert(grids[lev], nodal_flag_x), dmap[lev], 1, 1); - , rhoh0mac_new[lev][1].define( - convert(grids[lev], nodal_flag_y), dmap[lev], 1, 1); - , rhoh0mac_new[lev][2].define( - convert(grids[lev], nodal_flag_z), dmap[lev], 1, 1);); - - AMREX_D_TERM( - h0mac_new[lev][0].define(convert(grids[lev], nodal_flag_x), - dmap[lev], 1, 1); - , h0mac_new[lev][1].define( - convert(grids[lev], nodal_flag_y), dmap[lev], 1, 1); - , h0mac_new[lev][2].define( - convert(grids[lev], nodal_flag_z), dmap[lev], 1, 1);); - } - MakeS0mac(rho0_old, rho0mac_old); - MakeS0mac(rhoh0_old, rhoh0mac_old); - MakeS0mac(h0_old, h0mac_old); - MakeS0mac(rho0_new, rho0mac_new); - MakeS0mac(rhoh0_new, rhoh0mac_new); - MakeS0mac(h0_new, h0mac_new); - } -#endif - - // compute enthalpy fluxes - MakeRhoHFlux(scalold, sflux, sedge, umac, w0mac, rho0_old, - rho0_edge_old, rho0mac_old, rho0_new, rho0_edge_new, - rho0mac_new, rhoh0_old, rhoh0_edge_old, rhoh0mac_old, - rhoh0_new, rhoh0_edge_new, rhoh0mac_new, h0mac_old, - h0mac_new); - } - - for (int lev = 0; lev <= finest_level; ++lev) { - scal_force[lev].setVal(0.); - } - - ////////////////////////////////// - // Subtract w0 from MAC velocities - ////////////////////////////////// - - Addw0(umac, w0mac, -1.); - - //************************************************************************** - // 1) Create (rho h)' force at time n+1/2. - // (NOTE: we don't worry about filling ghost cells of the scal_force - // because we only need them in valid regions...) - // 2) Update (rho h) with conservative differencing. - //************************************************************************** - - MakeRhoHForce(scal_force, 0, thermal, umac, 0, which_step); - - // reaction forcing terms - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Add(scal_force[lev], intra[lev], RhoH, RhoH, 1, 0); - } - - ////////////////////////////////// - // Add w0 to MAC velocities - ////////////////////////////////// - - Addw0(umac, w0mac, 1.); - - Vector p0_new_cart(finest_level + 1); - for (int lev = 0; lev <= finest_level; ++lev) { - p0_new_cart[lev].define(grids[lev], dmap[lev], 1, 1); - } - - Put1dArrayOnCart(p0_new, p0_new_cart, false, false, bcs_f, 0); - - UpdateScal(scalold, scalnew, sflux, scal_force, RhoH, 1, p0_new_cart); -} diff --git a/Source/MaestroEvolve.cpp b/Source/MaestroEvolve.cpp index 4c552dc5f..6aa86858e 100644 --- a/Source/MaestroEvolve.cpp +++ b/Source/MaestroEvolve.cpp @@ -79,9 +79,6 @@ void Maestro::Evolve() { Real start_total = ParallelDescriptor::second(); // advance the solution by dt -#ifdef SDC - AdvanceTimeStepSDC(false); -#else if (use_exact_base_state || average_base_state) { // new temporal algorithm AdvanceTimeStepAverage(false); @@ -89,7 +86,6 @@ void Maestro::Evolve() { // original temporal algorithm AdvanceTimeStep(false); } -#endif t_old = t_new; diff --git a/Source/MaestroInit.cpp b/Source/MaestroInit.cpp index ab81fbe59..8871a33c3 100644 --- a/Source/MaestroInit.cpp +++ b/Source/MaestroInit.cpp @@ -73,10 +73,6 @@ void Maestro::Init() { } pi[lev].define(convert(grids[lev], nodal_flag), dmap[lev], 1, 0); // nodal -#ifdef SDC - intra[lev].define(grids[lev], dmap[lev], Nscal, 0); // for sdc - intra[lev].setVal(0.); -#endif } for (int lev = 0; lev <= finest_level; ++lev) { @@ -171,11 +167,7 @@ void Maestro::Init() { if (init_divu_iter > 0) { for (int i = 1; i <= init_divu_iter; ++i) { Print() << "Doing initial divu iteration #" << i << std::endl; -#ifdef SDC - DivuIterSDC(i); -#else DivuIter(i); -#endif } if (plot_int > 0 || plot_deltat > 0) { @@ -346,7 +338,6 @@ void Maestro::MakeNewLevelFromScratch(int lev, [[maybe_unused]] Real time, const rhcc_for_nodalproj[lev].define(ba, dm, 1, 1); pi[lev].define(convert(ba, nodal_flag), dm, 1, 0); // nodal - intra[lev].define(ba, dm, Nscal, 0); // for sdc sold[lev].setVal(0.); snew[lev].setVal(0.); @@ -359,7 +350,6 @@ void Maestro::MakeNewLevelFromScratch(int lev, [[maybe_unused]] Real time, const w0_cart[lev].setVal(0.); rhcc_for_nodalproj[lev].setVal(0.); pi[lev].setVal(0.); - intra[lev].setVal(0.); if (spherical) { normal[lev].define(ba, dm, 3, 1); @@ -473,11 +463,7 @@ void Maestro::InitProj() { delta_gamma1_term); // perform a nodal projection -#ifndef SDC NodalProj(initial_projection_comp, rhcc_for_nodalproj); -#else - NodalProj(initial_projection_comp, rhcc_for_nodalproj, false); -#endif } void Maestro::DivuIter(int istep_divu_iter) { @@ -615,139 +601,6 @@ void Maestro::DivuIter(int istep_divu_iter) { } } -// SDC -void Maestro::DivuIterSDC(int istep_divu_iter) { - // timer for profiling - BL_PROFILE_VAR("Maestro::DivuIterSDC()", DivuIterSDC); - - Vector stemp(finest_level + 1); - Vector rho_Hext(finest_level + 1); - Vector rho_omegadot(finest_level + 1); - Vector rho_Hnuc(finest_level + 1); - Vector thermal(finest_level + 1); - Vector rhohalf(finest_level + 1); - Vector Tcoeff(finest_level + 1); - Vector hcoeff(finest_level + 1); - Vector Xkcoeff(finest_level + 1); - Vector pcoeff(finest_level + 1); - Vector delta_gamma1(finest_level + 1); - Vector delta_gamma1_term(finest_level + 1); - Vector sdc_source(finest_level + 1); - - BaseState Sbar(base_geom.max_radial_level + 1, base_geom.nr_fine); - BaseState w0_force(base_geom.max_radial_level + 1, base_geom.nr_fine); - BaseState p0_minus_pthermbar(base_geom.max_radial_level + 1, - base_geom.nr_fine); - BaseState delta_gamma1_termbar(base_geom.max_radial_level + 1, - base_geom.nr_fine); - - Sbar.setVal(0.); - etarho_ec.setVal(0.0); - w0_force.setVal(0.0); - psi.setVal(0.0); - etarho_cc.setVal(0.0); - p0_minus_pthermbar.setVal(0.); - delta_gamma1_termbar.setVal(0.); - - for (int lev = 0; lev <= finest_level; ++lev) { - stemp[lev].define(grids[lev], dmap[lev], Nscal, 0); - rho_Hext[lev].define(grids[lev], dmap[lev], 1, 0); - rho_omegadot[lev].define(grids[lev], dmap[lev], NumSpec, 0); - rho_Hnuc[lev].define(grids[lev], dmap[lev], 1, 0); - thermal[lev].define(grids[lev], dmap[lev], 1, 0); - rhohalf[lev].define(grids[lev], dmap[lev], 1, 1); - Tcoeff[lev].define(grids[lev], dmap[lev], 1, 1); - hcoeff[lev].define(grids[lev], dmap[lev], 1, 1); - Xkcoeff[lev].define(grids[lev], dmap[lev], NumSpec, 1); - pcoeff[lev].define(grids[lev], dmap[lev], 1, 1); - delta_gamma1[lev].define(grids[lev], dmap[lev], 1, 1); - delta_gamma1_term[lev].define(grids[lev], dmap[lev], 1, 1); - sdc_source[lev].define(grids[lev], dmap[lev], Nscal, 0); - - // divu_iters do not use density weighting - rhohalf[lev].setVal(1.); - sdc_source[lev].setVal(0.); - } - - ReactSDC(sold, stemp, rho_Hext, p0_old, 0.5 * dt, t_old, sdc_source); - - if (use_thermal_diffusion) { - MakeThermalCoeffs(sold, Tcoeff, hcoeff, Xkcoeff, pcoeff); - - MakeExplicitThermal(thermal, sold, Tcoeff, hcoeff, Xkcoeff, pcoeff, - p0_old, temp_diffusion_formulation); - } else { - for (int lev = 0; lev <= finest_level; ++lev) { - thermal[lev].setVal(0.); - } - } - - MakeReactionRates(rho_omegadot, rho_Hnuc, sold); - - // compute S at cell-centers - Make_S_cc(S_cc_old, delta_gamma1_term, delta_gamma1, sold, uold, - rho_omegadot, rho_Hnuc, rho_Hext, thermal, p0_old, gamma1bar_old, - delta_gamma1_termbar); - - // NOTE: not sure if valid for use_exact_base_state - if (evolve_base_state) { - if ((use_exact_base_state || average_base_state) && - use_delta_gamma1_term) { - Sbar += delta_gamma1_termbar; - } else { - Average(S_cc_old, Sbar, 0); - - // compute Sbar = Sbar + delta_gamma1_termbar - if (use_delta_gamma1_term) { - Sbar += delta_gamma1_termbar; - } - - const auto is_predictor = true; - Makew0(w0, w0_force, Sbar, rho0_old, rho0_old, p0_old, p0_old, - gamma1bar_old, gamma1bar_old, p0_minus_pthermbar, dt, dt, - is_predictor); - } - } - - // make the nodal rhs for projection beta0*(S_cc-Sbar) + beta0*delta_chi - MakeRHCCforNodalProj(rhcc_for_nodalproj, S_cc_old, Sbar, beta0_old, - delta_gamma1_term); - - // perform a nodal projection - NodalProj(divu_iters_comp, rhcc_for_nodalproj, istep_divu_iter, false); - - Real dt_hold = dt; - - // compute new time step - EstDt(); - - if (maestro_verbose > 0) { - Print() << "Call to estdt at end of istep_divu_iter = " - << istep_divu_iter << " gives dt = " << dt << std::endl; - } - - dt *= init_shrink; - if (maestro_verbose > 0) { - Print() << "Multiplying dt by init_shrink; dt = " << dt << std::endl; - } - - if (dt > dt_hold) { - if (maestro_verbose > 0) { - Print() << "Ignoring this new dt since it's larger than the " - "previous dt = " - << dt_hold << std::endl; - } - dt = amrex::min(dt_hold, dt); - } - - if (fixed_dt != -1.0) { - // fixed dt - dt = fixed_dt; - if (maestro_verbose > 0) { - Print() << "Setting fixed dt = " << dt << std::endl; - } - } -} void Maestro::InitIter() { // timer for profiling diff --git a/Source/MaestroPlot.cpp b/Source/MaestroPlot.cpp index cfb669684..cfad73b74 100644 --- a/Source/MaestroPlot.cpp +++ b/Source/MaestroPlot.cpp @@ -349,19 +349,14 @@ Vector Maestro::PlotFileMF( Vector rho_Hext(finest_level + 1); Vector rho_omegadot(finest_level + 1); Vector rho_Hnuc(finest_level + 1); - Vector sdc_source(finest_level + 1); for (int lev = 0; lev <= finest_level; ++lev) { stemp[lev].define(grids[lev], dmap[lev], Nscal, 0); rho_Hext[lev].define(grids[lev], dmap[lev], 1, 0); rho_omegadot[lev].define(grids[lev], dmap[lev], NumSpec, 0); rho_Hnuc[lev].define(grids[lev], dmap[lev], 1, 0); - sdc_source[lev].define(grids[lev], dmap[lev], Nscal, 0); - - sdc_source[lev].setVal(0.); } -#ifndef SDC if (dt_in < small_dt) { React(s_in, stemp, rho_Hext, rho_omegadot, rho_Hnuc, p0_in, small_dt, t_in); @@ -369,15 +364,6 @@ Vector Maestro::PlotFileMF( React(s_in, stemp, rho_Hext, rho_omegadot, rho_Hnuc, p0_in, dt_in * 0.5, t_in); } -#else - if (dt_in < small_dt) { - ReactSDC(s_in, stemp, rho_Hext, p0_in, small_dt, t_in, sdc_source); - } else { - ReactSDC(s_in, stemp, rho_Hext, p0_in, dt_in * 0.5, t_in, sdc_source); - } - - MakeReactionRates(rho_omegadot, rho_Hnuc, s_in); -#endif if (plot_spec || plot_omegadot) { // omegadot diff --git a/Source/MaestroReact.cpp b/Source/MaestroReact.cpp index 19d09bd2f..c48546ab8 100644 --- a/Source/MaestroReact.cpp +++ b/Source/MaestroReact.cpp @@ -38,12 +38,10 @@ void Maestro::React(const Vector& s_in, Vector& s_out, // apply burning term if (do_burning) { -#ifndef SDC // do the burning, update rho_omegadot and rho_Hnuc // we pass in rho_Hext so that we can add it to rhoh in case we applied heating Burner(s_in, s_out, rho_Hext, rho_omegadot, rho_Hnuc, p0, dt_in, time_in); -#endif // pass temperature through for seeding the temperature update eos call for (int lev = 0; lev <= finest_level; ++lev) { MultiFab::Copy(s_out[lev], s_in[lev], Temp, Temp, 1, 0); @@ -79,337 +77,3 @@ void Maestro::React(const Vector& s_in, Vector& s_out, TfromRhoH(s_out, p0); } } - -// SDC subroutines -// compute heating term, rho_Hext, then -// react the state over dt_in and update s_out -void Maestro::ReactSDC(const Vector& s_in, Vector& s_out, - Vector& rho_Hext, const BaseState& p0, - const Real dt_in, const Real time_in, - Vector& source) { - - amrex::ignore_unused(time_in); - - // timer for profiling - BL_PROFILE_VAR("Maestro::ReactSDC()", ReactSDC); - - // external heating - if (do_heating) { - // computing heating term - MakeHeating(rho_Hext, s_in); - - if (do_burning) { - // add to source for enthalpy - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Add(source[lev], rho_Hext[lev], 0, RhoH, 1, 0); - } - } else { - // if we aren't burning, then we should just copy the old state to the - // new and only update the rhoh component with the heating term - // s_out = s_in + dt_in * rho_Hext - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Copy(s_out[lev], s_in[lev], 0, 0, Nscal, 0); - MultiFab::Saxpy(s_out[lev], dt_in, rho_Hext[lev], 0, RhoH, 1, - 0); - } - } - } else { - // not heating, so we zero rho_Hext - for (int lev = 0; lev <= finest_level; ++lev) { - rho_Hext[lev].setVal(0.); - } - } - - // apply burning term - if (do_burning) { - // copy s_in into s_out to fill coarse grids that are masked - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Copy(s_out[lev], s_in[lev], 0, 0, Nscal, 0); - } -#ifdef SDC - // do the burning, update s_out - Burner(s_in, s_out, p0, dt_in, time_in, source); -#endif - } - - // if we aren't doing any heating/burning, then just copy the old to the new - if (!do_heating && !do_burning) { - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Copy(s_out[lev], s_in[lev], 0, 0, Nscal, 0); - } - } - - // average down and fill ghost cells - AverageDown(s_out, 0, Nscal); - FillPatch(t_old, s_out, s_out, s_out, 0, 0, Nscal, 0, bcs_s); - - // average down (no ghost cells) - if (do_heating) { - AverageDown(rho_Hext, 0, 1); - } - - // now update temperature - if (use_tfromp) { - TfromRhoP(s_out, p0); - } else { - TfromRhoH(s_out, p0); - } -} - -// #ifndef SDC -// void Maestro::Burner(const Vector& s_in, Vector& s_out, -// const Vector& rho_Hext, -// Vector& rho_omegadot, Vector& rho_Hnuc, -// const BaseState& p0, const Real dt_in, -// const Real time_in) { -// // timer for profiling -// BL_PROFILE_VAR("Maestro::Burner()", Burner); - -// // Put tempbar_init on cart -// Vector tempbar_init_cart(finest_level + 1); - -// if (spherical) { -// for (int lev = 0; lev <= finest_level; ++lev) { -// tempbar_init_cart[lev].define(grids[lev], dmap[lev], 1, 0); -// tempbar_init_cart[lev].setVal(0.); -// } - -// if (drive_initial_convection) { -// Put1dArrayOnCart(tempbar_init, tempbar_init_cart, false, false, -// bcs_f, 0); -// } -// } - -// RealVector tempbar_init_vec((base_geom.max_radial_level + 1) * -// base_geom.nr_fine); -// if (!spherical) { -// tempbar_init.toVector(tempbar_init_vec); -// } - -// for (int lev = 0; lev <= finest_level; ++lev) { -// // get references to the MultiFabs at level lev -// const MultiFab& s_in_mf = s_in[lev]; -// MultiFab& s_out_mf = s_out[lev]; -// const MultiFab& rho_Hext_mf = rho_Hext[lev]; -// MultiFab& rho_omegadot_mf = rho_omegadot[lev]; -// MultiFab& rho_Hnuc_mf = rho_Hnuc[lev]; -// const MultiFab& tempbar_cart_mf = tempbar_init_cart[lev]; - -// // create mask assuming refinement ratio = 2 -// int finelev = lev + 1; -// if (lev == finest_level) { -// finelev = finest_level; -// } - -// const BoxArray& fba = s_in[finelev].boxArray(); -// const iMultiFab& mask = makeFineMask(s_in_mf, fba, IntVect(2)); - -// // loop over boxes (make sure mfi takes a cell-centered multifab as an argument) -// #ifdef _OPENMP -// #pragma omp parallel -// #endif -// for (MFIter mfi(s_in_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { -// // Get the index space of the valid region -// const Box& tileBox = mfi.tilebox(); - -// const bool use_mask = !(lev == finest_level); - -// // call fortran subroutine -// // use macros in AMReX_ArrayLim.H to pass in each FAB's data, -// // lo/hi coordinates (including ghost cells), and/or the # of components -// // We will also pass "validBox", which specifies the "valid" region. -// if (spherical) { -// #pragma gpu box(tileBox) -// burner_loop_sphr(AMREX_INT_ANYD(tileBox.loVect()), -// AMREX_INT_ANYD(tileBox.hiVect()), -// BL_TO_FORTRAN_ANYD(s_in_mf[mfi]), -// BL_TO_FORTRAN_ANYD(s_out_mf[mfi]), -// BL_TO_FORTRAN_ANYD(rho_Hext_mf[mfi]), -// BL_TO_FORTRAN_ANYD(rho_omegadot_mf[mfi]), -// BL_TO_FORTRAN_ANYD(rho_Hnuc_mf[mfi]), -// BL_TO_FORTRAN_ANYD(tempbar_cart_mf[mfi]), -// dt_in, time_in, BL_TO_FORTRAN_ANYD(mask[mfi]), -// (int)use_mask); -// } else { -// #pragma gpu box(tileBox) -// burner_loop(AMREX_INT_ANYD(tileBox.loVect()), -// AMREX_INT_ANYD(tileBox.hiVect()), lev, -// BL_TO_FORTRAN_ANYD(s_in_mf[mfi]), -// BL_TO_FORTRAN_ANYD(s_out_mf[mfi]), -// BL_TO_FORTRAN_ANYD(rho_Hext_mf[mfi]), -// BL_TO_FORTRAN_ANYD(rho_omegadot_mf[mfi]), -// BL_TO_FORTRAN_ANYD(rho_Hnuc_mf[mfi]), -// tempbar_init_vec.dataPtr(), dt_in, time_in, -// BL_TO_FORTRAN_ANYD(mask[mfi]), (int)use_mask); -// } -// } -// } -// } - -// #else -// // SDC burner -// void Maestro::Burner(const Vector& s_in, Vector& s_out, -// const BaseState& p0, const Real dt_in, -// const Real time_in, const Vector& source) { -// // timer for profiling -// BL_PROFILE_VAR("Maestro::BurnerSDC()", BurnerSDC); - -// Vector p0_cart(finest_level + 1); - -// // make a Fortran-friendly RealVector of p0 -// RealVector p0_vec((base_geom.max_radial_level + 1) * base_geom.nr_fine, -// 0.0); - -// if (spherical) { -// for (int lev = 0; lev <= finest_level; ++lev) { -// p0_cart[lev].define(grids[lev], dmap[lev], 1, 0); -// p0_cart[lev].setVal(0.); -// } - -// Put1dArrayOnCart(p0, p0_cart, false, false, bcs_f, 0); -// } else { -// p0.toVector(p0_vec); -// } - -// for (int lev = 0; lev <= finest_level; ++lev) { -// // get references to the MultiFabs at level lev -// const MultiFab& s_in_mf = s_in[lev]; -// MultiFab& s_out_mf = s_out[lev]; -// const MultiFab& p0_cart_mf = p0_cart[lev]; -// const MultiFab& source_mf = source[lev]; - -// // create mask assuming refinement ratio = 2 -// int finelev = lev + 1; -// if (lev == finest_level) finelev = finest_level; - -// const BoxArray& fba = s_in[finelev].boxArray(); -// const iMultiFab& mask = makeFineMask(s_in_mf, fba, IntVect(2)); - -// // loop over boxes (make sure mfi takes a cell-centered multifab as an argument) -// #ifdef _OPENMP -// #pragma omp parallel -// #endif -// for (MFIter mfi(s_in_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) { -// // Get the index space of the valid region -// const Box& tileBox = mfi.tilebox(); - -// int use_mask = !(lev == finest_level); - -// // call fortran subroutine - -// if (spherical) { -// #pragma gpu box(tileBox) -// burner_loop_sphr(AMREX_INT_ANYD(tileBox.loVect()), -// AMREX_INT_ANYD(tileBox.hiVect()), -// BL_TO_FORTRAN_ANYD(s_in_mf[mfi]), -// BL_TO_FORTRAN_ANYD(s_out_mf[mfi]), -// BL_TO_FORTRAN_ANYD(source_mf[mfi]), -// BL_TO_FORTRAN_ANYD(p0_cart_mf[mfi]), dt_in, -// time_in, BL_TO_FORTRAN_ANYD(mask[mfi]), -// use_mask); -// } else { -// #pragma gpu box(tileBox) -// burner_loop(AMREX_INT_ANYD(tileBox.loVect()), -// AMREX_INT_ANYD(tileBox.hiVect()), lev, -// BL_TO_FORTRAN_ANYD(s_in_mf[mfi]), -// BL_TO_FORTRAN_ANYD(s_out_mf[mfi]), -// BL_TO_FORTRAN_ANYD(source_mf[mfi]), -// p0_vec.dataPtr(), dt_in, time_in, -// BL_TO_FORTRAN_ANYD(mask[mfi]), use_mask); -// } -// } -// } -// } -// #endif - -// compute heating terms, rho_omegadot and rho_Hnuc -void Maestro::MakeReactionRates(Vector& rho_omegadot, - Vector& rho_Hnuc, - const Vector& scal) { - // timer for profiling - BL_PROFILE_VAR("Maestro::MakeReactionRates()", MakeReactionRates); - - const auto ispec_threshold = network_spec_index(burner_threshold_species); - - const Real temp_min = EOSData::mintemp; - const Real temp_max = EOSData::maxtemp; - - for (int lev = 0; lev <= finest_level; ++lev) { - if (!do_burning) { - rho_omegadot[lev].setVal(0.); - rho_Hnuc[lev].setVal(0.); - } else { - // loop over boxes (make sure mfi takes a cell-centered multifab as an argument) -#ifdef _OPENMP -#pragma omp parallel -#endif - for (MFIter mfi(scal[lev], TilingIfNotGPU()); mfi.isValid(); - ++mfi) { - // Get the index space of the valid region - const Box& tilebox = mfi.tilebox(); - - const Array4 rho_omegadot_arr = - rho_omegadot[lev].array(mfi); - const Array4 rho_Hnuc_arr = rho_Hnuc[lev].array(mfi); - const Array4 scal_arr = scal[lev].array(mfi); - - ParallelFor(tilebox, [=] AMREX_GPU_DEVICE(int i, int j, int k) { - auto rho = scal_arr(i, j, k, Rho); - Real x_in[NumSpec]; - for (auto n = 0; n < NumSpec; ++n) { - x_in[n] = scal_arr(i, j, k, FirstSpec + n) / rho; - } - Real x_test = - (ispec_threshold > 0) ? x_in[ispec_threshold] : 0.0; - - eos_t eos_state; - burn_t state; - - // if the threshold species is not in the network, then we burn - // normally. if it is in the network, make sure the mass - // fraction is above the cutoff. - if ((rho > burning_cutoff_density_lo && - rho < burning_cutoff_density_hi) && - (ispec_threshold < 0 || - (ispec_threshold > 0 && - x_test > burner_threshold_cutoff))) { - eos_state.rho = scal_arr(i, j, k, Rho); - for (auto n = 0; n < NumSpec; ++n) { - eos_state.xn[n] = scal_arr(i, j, k, FirstSpec + n) / - eos_state.rho; - } - eos_state.h = scal_arr(i, j, k, RhoH) / eos_state.rho; - eos_state.T = std::sqrt(temp_min * temp_max); - - // call the EOS with input rh to set T for rate evaluation - eos(eos_input_rh, eos_state); - eos_to_burn(eos_state, state); - - Array1D ydot; - actual_rhs(state, ydot); - // Note ydot is 1-based - - for (auto n = 0; n < NumSpec; ++n) { - rho_omegadot_arr(i, j, k, n) = - state.rho * aion[n] * ydot(1 + n); - } - // only necessary if nspec_evolve < nspec - // rho_omegadot_arr(i, j, k, NumSpec) = 0.0; - rho_Hnuc_arr(i, j, k) = state.rho * ydot(net_ienuc); - } else { - for (auto n = 0; n < NumSpec; ++n) { - rho_omegadot_arr(i, j, k, n) = 0.0; - } - rho_Hnuc_arr(i, j, k) = 0.0; - } - }); - } - } - } - - if (do_burning) { - // average down (no ghost cells) - AverageDown(rho_omegadot, 0, NumSpec); - AverageDown(rho_Hnuc, 0, 1); - } -} diff --git a/Source/MaestroRegrid.cpp b/Source/MaestroRegrid.cpp index 0c4bdf555..192a4e90f 100644 --- a/Source/MaestroRegrid.cpp +++ b/Source/MaestroRegrid.cpp @@ -223,9 +223,6 @@ void Maestro::RemakeLevel(int lev, Real time, const BoxArray& ba, const int ng_w = w0_cart[lev].nGrow(); const int ng_r = rhcc_for_nodalproj[lev].nGrow(); const int ng_p = pi[lev].nGrow(); -#ifdef SDC - const int ng_i = intra[lev].nGrow(); -#endif MultiFab snew_state(ba, dm, Nscal, ng_snew); MultiFab sold_state(ba, dm, Nscal, ng_snew); @@ -238,9 +235,6 @@ void Maestro::RemakeLevel(int lev, Real time, const BoxArray& ba, MultiFab w0_cart_state(ba, dm, AMREX_SPACEDIM, ng_w); MultiFab rhcc_for_nodalproj_state(ba, dm, 1, ng_r); MultiFab pi_state(convert(ba, nodal_flag), dm, 1, ng_p); -#ifdef SDC - MultiFab intra_state(ba, dm, Nscal, ng_i); -#endif FillPatch(lev, time, sold_state, sold, sold, 0, 0, Nscal, 0, bcs_s); std::swap(sold_state, sold[lev]); @@ -264,10 +258,6 @@ void Maestro::RemakeLevel(int lev, Real time, const BoxArray& ba, std::swap(w0_cart_state, w0_cart[lev]); std::swap(rhcc_for_nodalproj_state, rhcc_for_nodalproj[lev]); std::swap(pi_state, pi[lev]); -#ifdef SDC - FillPatch(lev, time, intra_state, intra, intra, 0, 0, Nscal, 0, bcs_f); - std::swap(intra_state, intra[lev]); -#endif if (spherical) { const int ng_n = normal[lev].nGrow(); @@ -304,9 +294,6 @@ void Maestro::MakeNewLevelFromCoarse(int lev, Real time, const BoxArray& ba, rhcc_for_nodalproj[lev].define(ba, dm, 1, 1); pi[lev].define(convert(ba, nodal_flag), dm, 1, 0); // nodal -#ifdef SDC - intra[lev].define(ba, dm, Nscal, 0); -#endif if (spherical) { normal[lev].define(ba, dm, 3, 1); @@ -325,9 +312,6 @@ void Maestro::MakeNewLevelFromCoarse(int lev, Real time, const BoxArray& ba, bcs_f); FillCoarsePatch(lev, time, gpi[lev], gpi, gpi, 0, 0, AMREX_SPACEDIM, bcs_f); FillCoarsePatch(lev, time, dSdt[lev], dSdt, dSdt, 0, 0, 1, bcs_f); -#ifdef SDC - FillCoarsePatch(lev, time, intra[lev], intra, intra, 0, 0, Nscal, bcs_f); -#endif } // within a call to AmrCore::regrid, this function deletes all data @@ -348,9 +332,6 @@ void Maestro::ClearLevel(int lev) { w0_cart[lev].clear(); rhcc_for_nodalproj[lev].clear(); pi[lev].clear(); -#ifdef SDC - intra[lev].clear(); -#endif if (spherical) { normal[lev].clear(); cell_cc_to_r[lev].clear(); diff --git a/Source/MaestroSetup.cpp b/Source/MaestroSetup.cpp index a53f6e2df..e609da6cd 100644 --- a/Source/MaestroSetup.cpp +++ b/Source/MaestroSetup.cpp @@ -180,7 +180,6 @@ void Maestro::Setup() { gpi.resize(max_level + 1); dSdt.resize(max_level + 1); pi.resize(max_level + 1); - intra.resize(max_level + 1); w0_cart.resize(max_level + 1); rhcc_for_nodalproj.resize(max_level + 1); normal.resize(max_level + 1); diff --git a/Source/MaestroThermal.cpp b/Source/MaestroThermal.cpp index af4cbcc5d..26d91f648 100644 --- a/Source/MaestroThermal.cpp +++ b/Source/MaestroThermal.cpp @@ -521,201 +521,3 @@ void Maestro::ThermalConduct(const Vector& s1, Vector& s2, // fill ghost cells FillPatch(t_old, s2, s2, s2, RhoH, RhoH, 1, RhoH, bcs_s); } - -// SDC version -void Maestro::ThermalConductSDC( - int which_step, const Vector& s_old, Vector& s_hat, - const Vector& s_new, const Vector& hcoeff1, - const Vector& Xkcoeff1, const Vector& pcoeff1, - const Vector& hcoeff2, const Vector& Xkcoeff2, - const Vector& pcoeff2) { - // timer for profiling - BL_PROFILE_VAR("Maestro::ThermalConductSDC()", ThermalConductSDC); - - // Dummy coefficient matrix, holds all zeros - Vector Dcoeff(finest_level + 1); - for (int lev = 0; lev <= finest_level; ++lev) { - Dcoeff[lev].define(grids[lev], dmap[lev], 1, 1); - Dcoeff[lev].setVal(0.); - } - - // solverrhs will hold solver RHS = (rho h)^hat - // dt/2 div . ( hcoeff1 grad h^1) - - // dt/2 sum_k div . (Xkcoeff2 grad X_k^hat + Xkcoeff1 grad X_k^1) - - // dt/2 div . ( pcoeff2 grad p_0^new + pcoeff1 grad p_0^old) - Vector solverrhs(finest_level + 1); - Vector resid(finest_level + 1); - for (int lev = 0; lev <= finest_level; ++lev) { - solverrhs[lev].define(grids[lev], dmap[lev], 1, 0); - resid[lev].define(grids[lev], dmap[lev], 1, 0); - } - - // compute RHS = (rho h)^hat - // = (rho h)^1 + dt . (aofs + intra) - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Copy(solverrhs[lev], s_hat[lev], RhoH, 0, 1, 0); - } - - // compute resid = div(hcoeff1 grad h^1) - sum_k div(Xkcoeff1 grad Xk^1) - div(pcoeff1 grad p0_old) - MakeExplicitThermal(resid, s_old, Dcoeff, hcoeff1, Xkcoeff1, pcoeff1, - p0_old, 2); - - // RHS = solverrhs + dt/2 * resid1 - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::LinComb(solverrhs[lev], 1.0, solverrhs[lev], 0, dt / 2.0, - resid[lev], 0, 0, 1, 0); - } - - // predictor - if (which_step == 1) { - // compute resid = 0 - sum_k div(Xkcoeff1 grad Xk^hat) - div(pcoeff1 grad p0_new) - MakeExplicitThermal(resid, s_hat, Dcoeff, Dcoeff, Xkcoeff1, pcoeff1, - p0_new, 2); - - } - // corrector - else if (which_step == 2) { - // compute resid = div(hcoeff2 grad h2) - MakeExplicitThermalHterm(resid, s_new, hcoeff2); - - // RHS = solverrhs + dt/2 * resid - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Saxpy(solverrhs[lev], -dt / 2.0, resid[lev], 0, 0, 1, 0); - } - - // compute resid = - sum_k div(Xkcoeff2 grad Xk^hat) - div(pcoeff2 grad p0_new) - MakeExplicitThermal(resid, s_hat, Dcoeff, Dcoeff, Xkcoeff2, pcoeff2, - p0_new, 2); - } - - // RHS = solverrhs + dt/2 * resid - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::LinComb(solverrhs[lev], 1.0, solverrhs[lev], 0, dt / 2.0, - resid[lev], 0, 0, 1, 0); - } - - // LHS coefficients for solver - Vector acoef(finest_level + 1); - Vector > face_bcoef(finest_level + 1); - Vector phi(finest_level + 1); - for (int lev = 0; lev <= finest_level; ++lev) { - acoef[lev].define(grids[lev], dmap[lev], 1, 1); - AMREX_D_TERM( - face_bcoef[lev][0].define(convert(grids[lev], nodal_flag_x), - dmap[lev], 1, 0); - , face_bcoef[lev][1].define(convert(grids[lev], nodal_flag_y), - dmap[lev], 1, 0); - , face_bcoef[lev][2].define(convert(grids[lev], nodal_flag_z), - dmap[lev], 1, 0);); - phi[lev].define(grids[lev], dmap[lev], 1, 1); - } - - // copy rho^hat=rho^new into acoef - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Copy(acoef[lev], s_hat[lev], Rho, 0, 1, 1); - } - - // average face-centered Bcoefficients - if (which_step == 1) { - PutDataOnFaces(hcoeff1, face_bcoef, true); - } else if (which_step == 2) { - PutDataOnFaces(hcoeff2, face_bcoef, true); - } - - // initialize value of phi to h^(2) as a guess - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Copy(phi[lev], s_hat[lev], RhoH, 0, 1, 1); - MultiFab::Divide(phi[lev], s_hat[lev], Rho, 0, 1, 1); - } - - // - // Set up implicit solve using MLABecLaplacian class - // - LPInfo info; - - // Only pass up to defined level to prevent looping over undefined grids. - MLABecLaplacian mlabec(Geom(0, finest_level), grids, dmap, info); - - // order of stencil - int linop_maxorder = 2; - mlabec.setMaxOrder(linop_maxorder); - - // set boundaries for mlabec using enthalpy bc's - std::array mlmg_lobc; - std::array mlmg_hibc; - - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - if (Geom(0).isPeriodic(idim)) { - mlmg_lobc[idim] = mlmg_hibc[idim] = LinOpBCType::Periodic; - } else { - // lo-side BCs - if (bcs_s[RhoH].lo(idim) == BCType::foextrap) { - mlmg_lobc[idim] = LinOpBCType::Neumann; - } else if (bcs_s[RhoH].lo(idim) == BCType::ext_dir) { - mlmg_lobc[idim] = LinOpBCType::Dirichlet; - } else { - mlmg_lobc[idim] = LinOpBCType::Neumann; - } - - // hi-side BCs - if (bcs_s[RhoH].hi(idim) == BCType::foextrap) { - mlmg_hibc[idim] = LinOpBCType::Neumann; - } else if (bcs_s[RhoH].hi(idim) == BCType::ext_dir) { - mlmg_hibc[idim] = LinOpBCType::Dirichlet; - } else { - mlmg_hibc[idim] = LinOpBCType::Neumann; - } - } - } - - mlabec.setDomainBC(mlmg_lobc, mlmg_hibc); - - for (int lev = 0; lev <= finest_level; ++lev) { - mlabec.setLevelBC(lev, &phi[lev]); - } - - // set timestep - if (which_step == 1) { - mlabec.setScalars(1.0, -dt / 2.0); - } else if (which_step == 2) { - mlabec.setScalars(1.0, -dt); - } - - for (int lev = 0; lev <= finest_level; ++lev) { - mlabec.setACoeffs(lev, acoef[lev]); - mlabec.setBCoeffs(lev, amrex::GetArrOfConstPtrs(face_bcoef[lev])); - } - - // solve A - timestep * div B grad phi = RHS - - // build an MLMG solver - MLMG thermal_mlmg(mlabec); - - // set solver parameters - thermal_mlmg.setVerbose(mg_verbose); - thermal_mlmg.setBottomVerbose(cg_verbose); - - // tolerance parameters taken from original MAESTRO fortran code - Real thermal_tol_abs = -1.e0; - for (int lev = 0; lev <= finest_level; ++lev) { - thermal_tol_abs = amrex::max(thermal_tol_abs, phi[lev].norm0()); - } - const Real solver_tol_abs = eps_mac * thermal_tol_abs; - const Real solver_tol_rel = eps_mac; - - // solve for phi - thermal_mlmg.solve(GetVecOfPtrs(phi), GetVecOfConstPtrs(solverrhs), - solver_tol_rel, solver_tol_abs); - - // load new rho*h into s2 - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Copy(s_hat[lev], phi[lev], 0, RhoH, 1, 1); - MultiFab::Multiply(s_hat[lev], s_hat[lev], Rho, RhoH, 1, 1); - } - - // average fine data onto coarser cells - AverageDown(s_hat, RhoH, 1); - - // fill ghost cells - FillPatch(t_old, s_hat, s_hat, s_hat, RhoH, RhoH, 1, RhoH, bcs_s); -} diff --git a/Source/Make.package b/Source/Make.package index b385aa4a1..50d4f4c3c 100644 --- a/Source/Make.package +++ b/Source/Make.package @@ -3,7 +3,6 @@ CEXE_sources += BaseStateGeometry.cpp CEXE_sources += Maestro.cpp CEXE_sources += MaestroAdvance.cpp CEXE_sources += MaestroAdvanceAvg.cpp -CEXE_sources += MaestroAdvanceSdc.cpp CEXE_sources += MaestroAdvectBase.cpp CEXE_sources += MaestroAdvection.cpp CEXE_sources += MaestroAverage.cpp diff --git a/Source/param/_cpp_parameters b/Source/param/_cpp_parameters index 0a8ffea2e..70c96f4c6 100644 --- a/Source/param/_cpp_parameters +++ b/Source/param/_cpp_parameters @@ -600,18 +600,6 @@ store_particle_vels bool false do_heating bool false -#----------------------------------------------------------------------------- -# category: SDC -#----------------------------------------------------------------------------- - -# how many SDC iterations to take (requires the code be compiled with -# {\tt SDC := t} -sdc_iters int 1 - -# recompute MAC velocity at the beginning of each SDC iter -sdc_couple_mac_velocity bool false - - #----------------------------------------------------------------------------- # category: GPU #-----------------------------------------------------------------------------