From 2d0c580c99cb6deef7a436ca8391be14b4af8168 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 5 Jul 2023 10:20:48 -0400 Subject: [PATCH] remove the SDC code paths (#394) this is all untested and doesn't work with the current Microphysics and we have no one pushing to use it at the moment. We can revive it from git in the future if needed --- CHANGES.md | 5 + Source/Maestro.H | 63 -- Source/MaestroAdvanceSdc.cpp | 1232 ----------------------------- Source/MaestroBurner.cpp | 179 ----- Source/MaestroCheckpoint.cpp | 8 - Source/MaestroDensityAdvance.cpp | 249 ------ Source/MaestroDiag.cpp | 14 - Source/MaestroEnthalpyAdvance.cpp | 370 --------- Source/MaestroEvolve.cpp | 4 - Source/MaestroInit.cpp | 147 ---- Source/MaestroPlot.cpp | 14 - Source/MaestroReact.cpp | 336 -------- Source/MaestroRegrid.cpp | 19 - Source/MaestroSetup.cpp | 1 - Source/MaestroThermal.cpp | 198 ----- Source/Make.package | 1 - Source/param/_cpp_parameters | 12 - 17 files changed, 5 insertions(+), 2847 deletions(-) delete mode 100644 Source/MaestroAdvanceSdc.cpp 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 #-----------------------------------------------------------------------------