Generalized multi-level relaxation (#1233)

* This has a bug on the second fill patch with two levels. * This ran. * hitting uninitialized data. * This runs but is not perfectly verified for correctness and has more of an affect on the boundaries than expected. * Hack update in set region for the fast integrator. * This is a mess but it's close. Need to fix the sponge factors near corners. The mask seems OK but the logic is hard to follow. * THIS WORKS! * fix plotfile and general clean up. * Remove G2E from slow_pre. Plan to remove update_interior_ghost; this will be overwritten with FIP. Need to remove G2E from fast_rhs. * Remove grids to evolve. Start style changes in microphysics and radiation. * General clean up. * More clean up. * Cleanup and fix for intermediate level with dynamic refinement. * Blocking factors * Use complement method to build mask. * Clean up blocking factors, modify inputs for dynamic ref, and Weiqun edits to DM/BA. * Update Source/BoundaryConditions/ERF_FillPatcher.cpp Co-authored-by: Weiqun Zhang <[email protected]> * Update Source/BoundaryConditions/ERF_FillPatcher.cpp Co-authored-by: Weiqun Zhang <[email protected]> * Update Source/BoundaryConditions/ERF_FillPatcher.cpp Co-authored-by: Weiqun Zhang <[email protected]> * Const and new fill function to be called by wrappers. * Reduce number of mask setvals. --------- Co-authored-by: Ann Almgren <[email protected]> Co-authored-by: Weiqun Zhang <[email protected]>
erf-model · Oct 4, 2023 · f0d8d23 · f0d8d23
1 parent 0598c89
commit f0d8d23
Show file tree

Hide file tree

Showing 29 changed files with 825 additions and 880 deletions.
diff --git a/...egTests/DynamicRefinement/inputs_onelevel → ...Tests/DynamicRefinement/inputs_threelevel b/...egTests/DynamicRefinement/inputs_onelevel → ...Tests/DynamicRefinement/inputs_threelevel
@@ -6,15 +6,13 @@ amrex.fpe_trap_invalid = 1
 fabarray.mfiter_tile_size = 1024 1024 1024
 
 # PROBLEM SIZE & GEOMETRY
-geometry.prob_lo     = -12  -12  -1
-geometry.prob_hi     =  12   12   1
-amr.n_cell           =  96   96   4
-
 geometry.prob_lo     =  -6  -6  -1
 geometry.prob_hi     =   6   6   1
 amr.n_cell           =  96  96   4
 
-amr.blocking_factor = 32 32 4
+amr.blocking_factor_x = 16  16  16
+amr.blocking_factor_y = 16  16  16
+amr.blocking_factor_z = 1 1 1
 
 geometry.is_periodic = 1 1 0
 
@@ -31,25 +29,29 @@ erf.v               = 1       # verbosity in ERF.cpp
 amr.v               = 1       # verbosity in Amr.cpp
 
 # REFINEMENT / REGRIDDING
-amr.max_level      = 1       # maximum level number allowed
-amr.ref_ratio_vect = 2 2 1
-erf.refinement_indicators  = hi_scal1
+amr.max_level      = 2       # maximum level number allowed
+amr.ref_ratio_vect = 2 2 1 4 4 1
+erf.refinement_indicators  = hi_scal1 hi_scal2
 erf.hi_scal1.max_level     = 1
 erf.hi_scal1.field_name    = scalar
 erf.hi_scal1.value_greater = 1. 
+erf.hi_scal2.max_level     = 2
+erf.hi_scal2.field_name    = scalar
+erf.hi_scal2.value_greater = 2.
 amr.n_error_buf  = 4
 erf.regrid_int   = 2 
-erf.cf_width     = 0 # Internal relaxation
-erf.cf_set_width = 0 # Internal Dirichlet
+erf.cf_width     = 4      # Internal relaxation (not ready w/ multiple boxes at a level)
+erf.cf_set_width = 1      # Internal Dirichlet  (not ready w/ multiple boxes at a level)
+amr.iterate_grids = true # Prevents level 1 BA decomposition (use true to decompose)
 
 # CHECKPOINT FILES
 erf.check_file      = chk        # root name of checkpoint file
 erf.check_int       = -100        # number of timesteps between checkpoints
 
 # PLOTFILES
 erf.plot_file_1     = plt        # prefix of plotfile name
-erf.plot_int_1      = 100        # number of timesteps between plotfiles
-erf.plot_vars_1     = density x_velocity y_velocity z_velocity pressure theta temp scalar pres_hse dens_hse pert_pres pert_dens
+erf.plot_int_1      = 1        # number of timesteps between plotfiles
+erf.plot_vars_1     = density rhotheta x_velocity y_velocity z_velocity pressure theta temp scalar pres_hse dens_hse pert_pres pert_dens mapfac
 
 # SOLVER CHOICE
 erf.alpha_T = 0.0

diff --git a/Exec/RegTests/DynamicRefinement/inputs_twolevel b/Exec/RegTests/DynamicRefinement/inputs_twolevel
@@ -6,15 +6,13 @@ amrex.fpe_trap_invalid = 1
 fabarray.mfiter_tile_size = 1024 1024 1024
 
 # PROBLEM SIZE & GEOMETRY
-geometry.prob_lo     = -12  -12  -1
-geometry.prob_hi     =  12   12   1
-amr.n_cell           =  96   96   4
-
 geometry.prob_lo     =  -6  -6  -1
 geometry.prob_hi     =   6   6   1
 amr.n_cell           =  96  96   4
 
-amr.blocking_factor = 32 32 4
+amr.blocking_factor_x = 32 32
+amr.blocking_factor_y = 32 32
+amr.blocking_factor_z = 1 1    
 
 geometry.is_periodic = 1 1 0
 
@@ -31,28 +29,24 @@ erf.v               = 1       # verbosity in ERF.cpp
 amr.v               = 1       # verbosity in Amr.cpp
 
 # REFINEMENT / REGRIDDING
-amr.max_level      = 2       # maximum level number allowed
-amr.ref_ratio_vect = 2 2 1 4 4 1
-erf.refinement_indicators  = hi_scal1 hi_scal2
+amr.max_level      = 1       # maximum level number allowed
+amr.ref_ratio_vect = 2 2 1
+erf.refinement_indicators  = hi_scal1
 erf.hi_scal1.max_level     = 1
 erf.hi_scal1.field_name    = scalar
 erf.hi_scal1.value_greater = 1. 
-erf.hi_scal2.max_level     = 2
-erf.hi_scal2.field_name    = scalar
-erf.hi_scal2.value_greater = 2.
 amr.n_error_buf  = 4
 erf.regrid_int   = 2 
-erf.cf_width     = 0      # Internal relaxation (not ready w/ multiple boxes at a level)
-erf.cf_set_width = 0      # Internal Dirichlet  (not ready w/ multiple boxes at a level)
-amr.iterate_grids = false # Prevents level 1 BA decomposition (use true to decompose)
+erf.cf_width     = 0 # Internal relaxation
+erf.cf_set_width = 0 # Internal Dirichlet
 
 # CHECKPOINT FILES
 erf.check_file      = chk        # root name of checkpoint file
 erf.check_int       = -100        # number of timesteps between checkpoints
 
 # PLOTFILES
 erf.plot_file_1     = plt        # prefix of plotfile name
-erf.plot_int_1      = 100        # number of timesteps between plotfiles
+erf.plot_int_1      = 50        # number of timesteps between plotfiles
 erf.plot_vars_1     = density x_velocity y_velocity z_velocity pressure theta temp scalar pres_hse dens_hse pert_pres pert_dens
 
 # SOLVER CHOICE

diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp
@@ -17,7 +17,7 @@ PhysBCFunctNoOp null_bc;
  */
 
 void
-ERF::FillPatch (int lev, Real time, const Vector<MultiFab*>& mfs)
+ERF::FillPatch (int lev, Real time, const Vector<MultiFab*>& mfs, bool fillset)
 {
     BL_PROFILE_VAR("ERF::FillPatch()",ERF_FillPatch);
     int bccomp;
@@ -73,11 +73,11 @@ ERF::FillPatch (int lev, Real time, const Vector<MultiFab*>& mfs)
     } // var_idx
 
     // Coarse-Fine set region
-    if (lev>0 && coupling_type=="OneWay" && cf_set_width>0) {
-        FPr_c[lev-1].fill(*mfs[Vars::cons], time, null_bc, domain_bcs_type, true);
-        FPr_u[lev-1].fill(*mfs[Vars::xvel], time, null_bc, domain_bcs_type, true);
-        FPr_v[lev-1].fill(*mfs[Vars::yvel], time, null_bc, domain_bcs_type, true);
-        FPr_w[lev-1].fill(*mfs[Vars::zvel], time, null_bc, domain_bcs_type, true);
+    if (lev>0 && coupling_type=="OneWay" && cf_set_width>0 && fillset) {
+        FPr_c[lev-1].FillSet(*mfs[Vars::cons], time, null_bc, domain_bcs_type);
+        FPr_u[lev-1].FillSet(*mfs[Vars::xvel], time, null_bc, domain_bcs_type);
+        FPr_v[lev-1].FillSet(*mfs[Vars::yvel], time, null_bc, domain_bcs_type);
+        FPr_w[lev-1].FillSet(*mfs[Vars::zvel], time, null_bc, domain_bcs_type);
     }
 
     // ***************************************************************************
@@ -228,10 +228,14 @@ ERF::FillIntermediatePatch (int lev, Real time,
 
     // Coarse-Fine set region
     if (lev>0 && coupling_type=="OneWay" && cf_set_width>0) {
-        FPr_c[lev-1].fill(*mfs[Vars::cons], time, null_bc, domain_bcs_type, true);
-        FPr_u[lev-1].fill(*mfs[Vars::xvel], time, null_bc, domain_bcs_type, true);
-        FPr_v[lev-1].fill(*mfs[Vars::yvel], time, null_bc, domain_bcs_type, true);
-        FPr_w[lev-1].fill(*mfs[Vars::zvel], time, null_bc, domain_bcs_type, true);
+        if (cons_only) {
+            FPr_c[lev-1].FillSet(*mfs[Vars::cons], time, null_bc, domain_bcs_type);
+        } else {
+            FPr_c[lev-1].FillSet(*mfs[Vars::cons], time, null_bc, domain_bcs_type);
+            FPr_u[lev-1].FillSet(*mfs[Vars::xvel], time, null_bc, domain_bcs_type);
+            FPr_v[lev-1].FillSet(*mfs[Vars::yvel], time, null_bc, domain_bcs_type);
+            FPr_w[lev-1].FillSet(*mfs[Vars::zvel], time, null_bc, domain_bcs_type);
+        }
     }
 
     // ***************************************************************************
@@ -251,10 +255,7 @@ ERF::FillIntermediatePatch (int lev, Real time,
     (*physbcs[lev])(mfs,icomp_cons,ncomp_cons,ngvect_cons,ngvect_vels,init_type,cons_only,BCVars::cons_bc,time);
     // ***************************************************************************
 
-    //
-    // It is important that we apply the MOST bcs after we have imposed all the others
-    //    so that we have enough information in the ghost cells to calculate the viscosity
-    //
+    // MOST boundary conditions
     if (!(cons_only && ncomp_cons == 1) && m_most && allow_most_bcs)
         m_most->impose_most_bcs(lev,mfs,eddyDiffs);
 }

diff --git a/Source/BoundaryConditions/ERF_FillPatcher.H b/Source/BoundaryConditions/ERF_FillPatcher.H
@@ -2,37 +2,57 @@
 #define ERF_FILLPATCHER_H_
 
 #include <AMReX_FillPatchUtil.H>
+#include <AMReX_Interp_C.H>
+#include <AMReX_MFInterp_C.H>
 
 class ERFFillPatcher
 {
 public:
 
-    ERFFillPatcher (amrex::BoxArray const& fba, amrex::DistributionMapping  fdm,
+    ERFFillPatcher (amrex::BoxArray const& fba, amrex::DistributionMapping const& fdm,
                     amrex::Geometry const& fgeom,
-                    amrex::BoxArray  cba, amrex::DistributionMapping  cdm,
+                    amrex::BoxArray const& cba, amrex::DistributionMapping const& cdm,
                     amrex::Geometry const& cgeom,
-                    int nghost, int nghost_subset, int ncomp, amrex::InterpBase* interp);
-
-    ~ERFFillPatcher ( )
-    {
-        delete m_cf_fine_data;
-        delete m_cf_fine_subset_data;
-        delete m_cf_crse_data[0];
-        delete m_cf_crse_data[1];
-    }
+                    int nghost, int nghost_set, int ncomp, amrex::InterpBase* interp);
 
-    void Define (amrex::BoxArray const& fba, amrex::DistributionMapping  fdm,
+    void Define (amrex::BoxArray const& fba, amrex::DistributionMapping const& fdm,
                  amrex::Geometry const& fgeom,
-                 amrex::BoxArray  cba, amrex::DistributionMapping  cdm,
+                 amrex::BoxArray const& cba, amrex::DistributionMapping const& cdm,
                  amrex::Geometry const& cgeom,
-                 int nghost, int nghost_subset, int ncomp, amrex::InterpBase* interp);
+                 int nghost, int nghost_sset, int ncomp,
+                 amrex::InterpBase* interp);
+
+    void BuildMask (amrex::BoxArray const& fba, int nghost, int nghost_set);
 
-    void registerCoarseData (amrex::Vector<amrex::MultiFab const*> const& crse_data,
+    void RegisterCoarseData (amrex::Vector<amrex::MultiFab const*> const& crse_data,
                              amrex::Vector<amrex::Real> const& crse_time);
 
+    void InterpFace (amrex::MultiFab& fine,
+                     amrex::MultiFab const& crse,
+                     int mask_val);
+
+    void InterpCell (amrex::MultiFab& fine,
+                     amrex::MultiFab const& crse,
+                     amrex::Vector<amrex::BCRec> const& bcr,
+                     int mask_val);
+
+    int GetSetMaskVal () { return m_set_mask; }
+
+    int GetRelaxMaskVal () { return m_relax_mask; }
+
+    amrex::iMultiFab* GetMask () { return m_cf_mask.get(); }
+
     template <typename BC>
-    void fill (amrex::MultiFab& mf, amrex::Real time,
-              BC& cbc, amrex::Vector<amrex::BCRec> const& bcs, bool fill_subset=false);
+    void FillSet (amrex::MultiFab& mf, amrex::Real time,
+                  BC& cbc, amrex::Vector<amrex::BCRec> const& bcs);
+
+    template <typename BC>
+    void FillRelax (amrex::MultiFab& mf, amrex::Real time,
+                    BC& cbc, amrex::Vector<amrex::BCRec> const& bcs);
+
+    template <typename BC>
+    void Fill (amrex::MultiFab& mf, amrex::Real time,
+               BC& cbc, amrex::Vector<amrex::BCRec> const& bcs, int mask_val);
 
 private:
 
@@ -47,15 +67,17 @@ private:
     int m_ncomp;
     amrex::InterpBase* m_interp;
     amrex::IntVect m_ratio;
-    amrex::Vector<amrex::MultiFab*> m_cf_crse_data;
-    amrex::MultiFab* m_cf_fine_data;
-    amrex::MultiFab* m_cf_fine_subset_data;
+    std::unique_ptr<amrex::MultiFab> m_cf_crse_data_old;
+    std::unique_ptr<amrex::MultiFab> m_cf_crse_data_new;
+    std::unique_ptr<amrex::iMultiFab> m_cf_mask;
     amrex::Vector<amrex::Real> m_crse_times;
     amrex::Real m_dt_crse;
+    int m_set_mask{2};
+    int m_relax_mask{1};
 };
 
 /*
- * Fill data at the coarse-fine boundary
+ * Fill fine data in the set region
  *
  * @param[out] mf    MultiFab to be filled
  * @param[in]  time  Time at which to fill data
@@ -64,8 +86,41 @@ private:
  */
 template <typename BC>
 void
-ERFFillPatcher::fill (amrex::MultiFab& mf, amrex::Real time,
-                      BC& cbc, amrex::Vector<amrex::BCRec> const& bcs, bool fill_subset)
+ERFFillPatcher::FillSet (amrex::MultiFab& mf, amrex::Real time,
+                         BC& cbc, amrex::Vector<amrex::BCRec> const& bcs)
+{
+    Fill(mf,time,cbc,bcs,m_set_mask);
+}
+
+/*
+ * Fill fine data in the relax region
+ *
+ * @param[out] mf    MultiFab to be filled
+ * @param[in]  time  Time at which to fill data
+ * @param[in]  cbc   Coarse boundary condition
+ * @param[in]  bcs   Vector of boundary conditions
+ */
+template <typename BC>
+void
+ERFFillPatcher::FillRelax (amrex::MultiFab& mf, amrex::Real time,
+                           BC& cbc, amrex::Vector<amrex::BCRec> const& bcs)
+{
+    Fill(mf,time,cbc,bcs,m_relax_mask);
+}
+
+/*
+ * Fill fine data in the relax region
+ *
+ * @param[out] mf    MultiFab to be filled
+ * @param[in]  time  Time at which to fill data
+ * @param[in]  cbc   Coarse boundary condition
+ * @param[in]  bcs   Vector of boundary conditions
+ * @param[in]  mask_val Value to assign mask array
+ */
+template <typename BC>
+void
+ERFFillPatcher::Fill (amrex::MultiFab& mf, amrex::Real time,
+                      BC& cbc, amrex::Vector<amrex::BCRec> const& bcs, int mask_val)
 {
     constexpr amrex::Real eps = std::numeric_limits<float>::epsilon();
     AMREX_ALWAYS_ASSERT((time >= m_crse_times[0]-eps) && (time <= m_crse_times[1]+eps));
@@ -75,32 +130,28 @@ ERFFillPatcher::fill (amrex::MultiFab& mf, amrex::Real time,
     amrex::Real fac_old = 1.0 - fac_new;
 
     // Boundary condition operator
-    cbc(*(m_cf_crse_data[0]), 0, m_ncomp, amrex::IntVect(0), time, 0);
+    cbc(*(m_cf_crse_data_old), 0, m_ncomp, amrex::IntVect(0), time, 0);
 
     // Coarse MF to hold time interpolated data
-    amrex::MultiFab crse_data_time_interp(m_cf_crse_data[0]->boxArray(), m_cf_crse_data[0]->DistributionMap(),
+    amrex::MultiFab crse_data_time_interp(m_cf_crse_data_old->boxArray(),
+                                          m_cf_crse_data_old->DistributionMap(),
                                           m_ncomp, amrex::IntVect{0});
 
     // Time interpolate the coarse data
     amrex::MultiFab::LinComb(crse_data_time_interp,
-                             fac_old, *(m_cf_crse_data[0]), 0,
-                             fac_new, *(m_cf_crse_data[1]), 0,
+                             fac_old, *(m_cf_crse_data_old), 0,
+                             fac_new, *(m_cf_crse_data_new), 0,
                              0, m_ncomp, amrex::IntVect{0});
 
-    // Ensure fine domain box is correct index type
-    amrex::Box fdest_dom = amrex::convert(m_fgeom.Domain(),m_cf_fine_data->boxArray().ixType());
-
-    // Spatially interpolate the time-interpolated coarse data
-    amrex::FillPatchInterp(*m_cf_fine_data, 0, crse_data_time_interp, 0, m_ncomp, amrex::IntVect(0),
-                           m_cgeom, m_fgeom, fdest_dom, m_ratio, m_interp, bcs, 0);
-
-    // Fill whole region or subset?
-    if (fill_subset) {
-        //amrex::MultiFab::Copy(m_cf_fine_subset_data, m_cf_fine_data, 0, 0, m_ncomp, amrex::IntVect{0});
-        m_cf_fine_subset_data->ParallelCopy(*m_cf_fine_data, 0, 0, m_ncomp, amrex::IntVect{0}, amrex::IntVect{0});
-        mf.ParallelCopy(*m_cf_fine_subset_data, 0, 0, m_ncomp, amrex::IntVect{0}, amrex::IntVect{0});
+    // Call correct spatial interpolation type
+    amrex::IndexType m_ixt = mf.boxArray().ixType();
+    int ixt_sum = m_ixt[0]+m_ixt[1]+m_ixt[2];
+    if (ixt_sum == 0) {
+        InterpCell(mf,crse_data_time_interp,bcs,mask_val);
+    } else if (ixt_sum == 1) {
+        InterpFace(mf,crse_data_time_interp,mask_val);
     } else {
-        mf.ParallelCopy(*m_cf_fine_data, 0, 0, m_ncomp, amrex::IntVect{0}, amrex::IntVect{0});
+        amrex::Abort("ERF_FillPatcher only supports face linear and cell cons linear interp!");
     }
 }
 #endif