Merge branch 'development' into rm_cfl_violation

CHANGES.md

@@ -1,3 +1,18 @@
+# 23.08
+
+  * Time evolution without subcycling on the fine levels, which is enabled via
+    the runtime parameter amr.subcycling_mode = "None", has been significantly
+    rewritten to improve computational performance for certain cases. When the
+    fine levels do not subcycle, the evolution can be done by processing all
+    of the hydro updates on the fine level together and then immediately doing
+    the flux correction to sync the coarse and fine level fluxes at the
+    boundary between levels. This is how many AMR codes that do not subcycle
+    are written. Castro now does this when the user chooses not to subcycle.
+    The benefit of this approach is most evidence for problems with Poisson
+    gravity that use a multigrid solve, as we can significantly reduce the
+    number of Poisson solves per step, performing only a single composite
+    (multi-level) solve at the new-time simultaneously for all levels. (#2505)
+
 # 23.07
 
   * The parameter castro.state_nghost, which allowed State_Type to have ghost

Source/diffusion/Castro_diffusion.cpp

@@ -31,7 +31,7 @@ Castro::construct_old_diff_source(MultiFab& source, MultiFab& state_in, Real tim
 #endif
         ParallelDescriptor::ReduceRealMax(run_time,IOProc);
 
-        amrex::Print() << "Castro::construct_old_diff_source() time = " << run_time << "\n" << "\n";
+        amrex::Print() << "Castro::construct_old_diff_source() time = " << run_time << " on level " << level << "\n" << "\n";
 #ifdef BL_LAZY
         });
 #endif
@@ -68,7 +68,7 @@ Castro::construct_new_diff_source(MultiFab& source, MultiFab& state_old, MultiFa
 #endif
         ParallelDescriptor::ReduceRealMax(run_time,IOProc);
 
-        amrex::Print() << "Castro::construct_new_diff_source() time = " << run_time << "\n" << "\n";
+        amrex::Print() << "Castro::construct_new_diff_source() time = " << run_time << " on level " << level << "\n" << "\n";
 #ifdef BL_LAZY
         });
 #endif

Source/driver/Castro.H

@@ -1036,8 +1036,9 @@ public:
 ///
 /// @param crse_level
 /// @param fine_level
+/// @param in_post_timestep
 ///
-    void reflux (int crse_level, int fine_level);
+    void reflux (int crse_level, int fine_level, bool in_post_timestep);
 
 
 ///

Source/driver/Castro.cpp

@@ -1951,7 +1951,9 @@ Castro::post_timestep (int iteration_local)
     // will also do the sync solve associated with the reflux.
 
     if (do_reflux && level < parent->finestLevel()) {
-      reflux(level, level+1);
+        if (parent->subcyclingMode() != "None") {
+            reflux(level, level+1, true);
+        }
     }
 
     // Ensure consistency with finer grids.
@@ -2614,9 +2616,31 @@ Castro::FluxRegFineAdd() {
 
 }
 
+// reflux() synchronizes fluxes between levels and has two modes of operation.
+//
+// When in_post_timestep = true, we are performing the reflux in AmrLevel's
+// post_timestep() routine. This takes place after this level has completed
+// its advance as well as all fine levels above this level have completed
+// their advance. The main activities are:
+//
+// (1) Limit the flux correction so that it won't cause unphysical densities.
+// (2) Update any copies of the full flux arrays with the flux correction.
+// (3) Perform the flux correction on StateData using FluxRegister's Reflux().
+// (4) Do a sync solve for Poisson gravity now that the density has changed.
+// (5) Recompute the new-time sources now that StateData has changed.
+//
+// When in_post_timestep = false, we are performing the reflux immediately
+// after the hydro step has taken place on all levels between crse_level and
+// fine_level. This happens when amr.subcycling_mode = None, or in general when
+// the fine level does not subcycle relative to the coarse level. In this mode
+// we can skip steps (4) and (5), as those steps are only necessary to fix the
+// updates on these levels that took place without full information about the
+// hydro solve being synchronized. The advance on crse_level and fine_level will
+// then be followed by a normal computation of the new-time gravity and new-time
+// sources, which will not need to be corrected later.
 
 void
-Castro::reflux(int crse_level, int fine_level)
+Castro::reflux (int crse_level, int fine_level, bool in_post_timestep)
 {
     BL_PROFILE("Castro::reflux()");
 
@@ -2630,7 +2654,7 @@ Castro::reflux(int crse_level, int fine_level)
     Vector<std::unique_ptr<MultiFab> > drho(nlevs);
     Vector<std::unique_ptr<MultiFab> > dphi(nlevs);
 
-    if (do_grav && gravity->get_gravity_type() == "PoissonGrav" && gravity->NoSync() == 0)  {
+    if (do_grav && gravity->get_gravity_type() == "PoissonGrav" && gravity->NoSync() == 0 && in_post_timestep)  {
 
         for (int lev = crse_level; lev <= fine_level; ++lev) {
 
@@ -2783,7 +2807,7 @@ Castro::reflux(int crse_level, int fine_level)
             // Update the coarse fluxes MultiFabs using the reflux data. This should only make
             // a difference if we re-evaluate the source terms later.
 
-            if (update_sources_after_reflux) {
+            if (update_sources_after_reflux || !in_post_timestep) {
 
                 MultiFab::Add(*crse_lev.fluxes[idir], temp_fluxes[idir], 0, 0, crse_lev.fluxes[idir]->nComp(), 0);
 
@@ -2803,7 +2827,7 @@ Castro::reflux(int crse_level, int fine_level)
 #ifdef GRAVITY
         int ilev = lev - crse_level - 1;
 
-        if (do_grav && gravity->get_gravity_type() == "PoissonGrav" && gravity->NoSync() == 0) {
+        if (do_grav && gravity->get_gravity_type() == "PoissonGrav" && gravity->NoSync() == 0 && in_post_timestep) {
             reg->Reflux(*drho[ilev], crse_lev.volume, 1.0, 0, URHO, 1, crse_lev.geom);
             amrex::average_down(*drho[ilev + 1], *drho[ilev], 0, 1, getLevel(lev).crse_ratio);
         }
@@ -2826,7 +2850,7 @@ Castro::reflux(int crse_level, int fine_level)
 
             reg->Reflux(crse_state, dr, 1.0, 0, UMX, 1, crse_lev.geom);
 
-            if (update_sources_after_reflux) {
+            if (update_sources_after_reflux || !in_post_timestep) {
 
                 MultiFab tmp_fluxes(crse_lev.P_radial.boxArray(),
                                     crse_lev.P_radial.DistributionMap(),
@@ -2863,7 +2887,7 @@ Castro::reflux(int crse_level, int fine_level)
 
             reg->Reflux(crse_lev.get_new_data(Rad_Type), crse_lev.volume, 1.0, 0, 0, Radiation::nGroups, crse_lev.geom);
 
-            if (update_sources_after_reflux) {
+            if (update_sources_after_reflux || !in_post_timestep) {
 
                 for (int idir = 0; idir < AMREX_SPACEDIM; ++idir) {
 
@@ -2893,7 +2917,7 @@ Castro::reflux(int crse_level, int fine_level)
 #endif
 
 #ifdef GRAVITY
-        if (do_grav && gravity->get_gravity_type() == "PoissonGrav" && gravity->NoSync() == 0)  {
+        if (do_grav && gravity->get_gravity_type() == "PoissonGrav" && gravity->NoSync() == 0 && in_post_timestep)  {
 
             reg = &getLevel(lev).phi_reg;
 
@@ -2922,7 +2946,7 @@ Castro::reflux(int crse_level, int fine_level)
     // Do the sync solve across all levels.
 
 #ifdef GRAVITY
-    if (do_grav && gravity->get_gravity_type() == "PoissonGrav" && gravity->NoSync() == 0) {
+    if (do_grav && gravity->get_gravity_type() == "PoissonGrav" && gravity->NoSync() == 0 && in_post_timestep) {
       gravity->gravity_sync(crse_level, fine_level, amrex::GetVecOfPtrs(drho), amrex::GetVecOfPtrs(dphi));
     }
 #endif
@@ -2938,7 +2962,7 @@ Castro::reflux(int crse_level, int fine_level)
     // ghost zone fills like diffusion depend on the data in the
     // coarser levels.
 
-    if (update_sources_after_reflux &&
+    if (update_sources_after_reflux && in_post_timestep &&
         (time_integration_method == CornerTransportUpwind ||
          time_integration_method == SimplifiedSpectralDeferredCorrections)) {
 

Source/driver/Castro_advance.cpp

@@ -38,6 +38,11 @@ Castro::advance (Real time,
   //    amr_ncycle    : the number of subcycles at this level
 
 {
+    if (parent->subcyclingMode() == "None" && level > 0) {
+        amrex::Print() << "\n  Advance at this level has already been completed.\n\n";
+        return dt;
+    }
+
     BL_PROFILE("Castro::advance()");
 
     // Save the wall time when we started the step.
@@ -48,7 +53,15 @@ Castro::advance (Real time,
 
     Real dt_new = dt;
 
-    initialize_advance(time, dt, amr_iteration);
+    int max_level_to_advance = level;
+
+    if (parent->subcyclingMode() == "None" && level == 0) {
+        max_level_to_advance = parent->finestLevel();
+    }
+
+    for (int lev = level; lev <= max_level_to_advance; ++lev) {
+        getLevel(lev).initialize_advance(time, dt, amr_iteration);
+    }
 
     // Do the advance.
 
@@ -77,27 +90,29 @@ Castro::advance (Real time,
         post_step_regrid = 1;
     }
 
+    for (int lev = level; lev <= max_level_to_advance; ++lev) {
 #ifdef AUX_UPDATE
-    advance_aux(time, dt);
+        getLevel(lev).advance_aux(time, dt);
 #endif
 
 #ifdef GRAVITY
-    // Update the point mass.
-    if (use_point_mass == 1) {
-        pointmass_update(time, dt);
-    }
+        // Update the point mass.
+        if (use_point_mass == 1) {
+            getLevel(lev).pointmass_update(time, dt);
+        }
 #endif
 
 #ifdef RADIATION
-    MultiFab& S_new = get_new_data(State_Type);
-    final_radiation_call(S_new, amr_iteration, amr_ncycle);
+        MultiFab& S_new = getLevel(lev).get_new_data(State_Type);
+        getLevel(lev).final_radiation_call(S_new, amr_iteration, amr_ncycle);
 #endif
 
 #ifdef AMREX_PARTICLES
-    advance_particles(amr_iteration, time, dt);
+        getLevel(lev).advance_particles(amr_iteration, time, dt);
 #endif
 
-    finalize_advance();
+        getLevel(lev).finalize_advance();
+    }
 
     return dt_new;
 }
@@ -560,7 +575,7 @@ Castro::finalize_advance()
 {
     BL_PROFILE("Castro::finalize_advance()");
 
-    if (do_reflux == 1) {
+    if (do_reflux == 1 && parent->subcyclingMode() != "None") {
         FluxRegCrseInit();
         FluxRegFineAdd();
     }
@@ -570,7 +585,7 @@ Castro::finalize_advance()
     // the fluxes from the full timestep (this will be used
     // later during the reflux operation).
 
-    if (do_reflux == 1 && update_sources_after_reflux == 1) {
+    if (do_reflux == 1 && update_sources_after_reflux == 1 && parent->subcyclingMode() != "None") {
         for (int idir = 0; idir < AMREX_SPACEDIM; ++idir) {
             MultiFab::Copy(*mass_fluxes[idir], *fluxes[idir], URHO, 0, 1, 0);
         }
@@ -614,14 +629,34 @@ Castro::finalize_advance()
 
     // Record how many zones we have advanced.
 
-    num_zones_advanced += static_cast<Real>(grids.numPts()) / static_cast<Real>(getLevel(0).grids.numPts());
+    int max_level_to_advance = level;
+
+    if (parent->subcyclingMode() == "None") {
+        max_level_to_advance = parent->finestLevel();
+    }
+
+    long num_pts_advanced = 0;
+
+    for (int lev = level; lev <= max_level_to_advance; ++lev) {
+        num_pts_advanced += getLevel(lev).grids.numPts();
+    }
+
+    num_zones_advanced += static_cast<Real>(num_pts_advanced) / static_cast<Real>(getLevel(0).grids.numPts());
 
     Real wall_time = ParallelDescriptor::second() - wall_time_start;
-    Real fom_advance = static_cast<Real>(grids.numPts()) / wall_time / 1.e6;
+
+    Real fom_advance = static_cast<Real>(num_pts_advanced) / wall_time / 1.e6;
 
     if (verbose >= 1) {
-        amrex::Print() << "  Zones advanced per microsecond at this level: "
-                       << fom_advance << std::endl << std::endl;
+        if (max_level_to_advance > 0) {
+            if (level == 0) {
+                amrex::Print() << "  Zones advanced per microsecond from level " << level << " to level "
+                               << max_level_to_advance << ": " << fom_advance << std::endl << std::endl;
+            }
+        }
+        else {
+            amrex::Print() << "  Zones advanced per microsecond at this level: "
+                           << fom_advance << std::endl << std::endl;
+        }
     }
-
 }