Towards Dynamic Mesh Refinement (#1172)

* Allow only filling of external ghost cells and enable dynamic mesh refinement with relaxation/set regions.

* Delete backtrace file.

* Clean up.

* Relocate crse register after it has been fillpatched.

---------

Co-authored-by: Aaron Lattanzi <[email protected]>

Docs/sphinx_doc/MeshRefinement.rst

@@ -135,9 +135,9 @@ and face-baced normal momenta on the coarse-fine interface, the coarse data is c
 interpolated to the fine mesh. The interpolated data is utilized to specify ghost cell data
 (outside of the valid fine region) as well as specify and relax data inside the lateral boundaries
 of the fine region. More specifically, a user may specify the total width of the interior
-Dirichlet and relaxation region with ``erf.wrfbdy_width = <Int>`` (yellow + blue)
+Dirichlet and relaxation region with ``erf.cf_width = <Int>`` (yellow + blue)
 and analogously the width of the interior Dirichlet region may be specified with
-``erf.wrfbdy_set_width = <Int>`` (yellow).
+``erf.cf_set_width = <Int>`` (yellow).
 
 .. |wrfbdy| image:: figures/wrfbdy_BCs.png
  :width: 600
@@ -169,7 +169,12 @@ the RHS (:math:`F`) is given by the following:
 
 where :math:`G` is the RHS of the NS equations, :math:`\psi^{\prime}` is the predicted update without
 relaxation, :math:`\psi^{FP}` is the fine patch data obtained from space-time interpolation of the
-coarse mesh, and :math:`n` is the minimum number of grid point from a lateral boundary.
+coarse mesh, and :math:`n` is the minimum number of grid point from a lateral boundary. Finally, we
+note that time dependent Dirichlet data, provided via an external file, may be enforced on the
+lateral boundary conditions of the domain (coarsest mesh). For such cases, the relaxation region width
+at the domain edges may be specified with ``erf.wrfbdy_width = <Int>`` (yellow + blue) while the
+interior Dirichlet region may be specified with ``erf.wrfbdy_set_width = <Int>`` (yellow).
+
 
 By two-way coupling, we mean that in additional to specifying ghost cell data (outside of the valid fine region),
 the fine mesh communicates data back to the coarse mesh in two ways:

Exec/RegTests/DynamicRefinement/inputs

@@ -78,9 +78,8 @@ erf.hi_scal2.max_level = 2
 erf.hi_scal2.field_name = scalar
 erf.hi_scal2.value_greater = 10.
 
-amr.n_error_buf = 4
-erf.regrid_int = 2 
-erf.wrfbdy_width = 4
-erf.wrfbdy_set_width = 1
+amr.n_error_buf  = 4
+erf.regrid_int  = 2 
+erf.cf_width  = 4
+erf.cf_set_width = 1
 
-erf.use_NumDiff = 0.9

Source/BoundaryConditions/ERF_FillPatch.cpp

@@ -73,7 +73,7 @@ ERF::FillPatch (int lev, Real time, const Vector<MultiFab*>& mfs)
  } // var_idx
 
  // Coarse-Fine set region
- if (lev>0 && coupling_type=="OneWay") {
+ 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);
@@ -227,7 +227,7 @@ ERF::FillIntermediatePatch (int lev, Real time,
  } // var_idx
 
  // Coarse-Fine set region
- if (lev>0 && coupling_type=="OneWay") {
+ 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);

Source/BoundaryConditions/ERF_FillPatcher.H

@@ -13,6 +13,20 @@ public:
  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];
+ }
+
+ void Define (amrex::BoxArray const& fba, amrex::DistributionMapping fdm,
+ amrex::Geometry const& fgeom,
+ amrex::BoxArray cba, amrex::DistributionMapping cdm,
+ amrex::Geometry const& cgeom,
+ int nghost, int nghost_subset, int ncomp, amrex::InterpBase* interp);
+
  void registerCoarseData (amrex::Vector<amrex::MultiFab const*> const& crse_data,
  amrex::Vector<amrex::Real> const& crse_time);
 
@@ -33,9 +47,9 @@ 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;
+ amrex::Vector<amrex::MultiFab*> m_cf_crse_data;
+ amrex::MultiFab* m_cf_fine_data;
+ amrex::MultiFab* m_cf_fine_subset_data;
  amrex::Vector<amrex::Real> m_crse_times;
  amrex::Real m_dt_crse;
 };
@@ -61,34 +75,32 @@ 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[0]), 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[0]->boxArray(), m_cf_crse_data[0]->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[0]), 0,
+ fac_new, *(m_cf_crse_data[1]), 0,
  0, m_ncomp, 0);
 
  // Ensure fine domain box is correct index type
- amrex::Box fdest_dom = amrex::convert(m_fgeom.Domain(),m_cf_fine_data.boxArray().ixType());
+ 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),
+ 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}, mf.nGrowVect());
+ 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}, mf.nGrowVect());
  } else {
- mf.ParallelCopy(m_cf_fine_data, 0, 0, m_ncomp, amrex::IntVect{0}, mf.nGrowVect());
+ mf.ParallelCopy(*m_cf_fine_data, 0, 0, m_ncomp, amrex::IntVect{0}, mf.nGrowVect());
  }
-
-
 }
 #endif

Source/BoundaryConditions/ERF_FillPatcher.cpp

@@ -23,23 +23,64 @@ ERFFillPatcher::ERFFillPatcher (BoxArray const& fba, DistributionMapping fdm,
  int nghost, int nghost_subset,
  int ncomp, InterpBase* interp)
  : m_fba(fba),
- m_cba(std::move(cba)),
- m_fdm(std::move(fdm)),
- m_cdm(std::move(cdm)),
+ m_cba(cba),
+ m_fdm(fdm),
+ m_cdm(cdm),
  m_fgeom(fgeom),
- m_cgeom(cgeom),
- m_nghost(nghost),
- m_nghost_subset(nghost_subset),
- m_ncomp(ncomp),
- m_interp(interp)
+ m_cgeom(cgeom)
 {
- AMREX_ALWAYS_ASSERT(nghost < 0);
- AMREX_ALWAYS_ASSERT(nghost_subset < 0);
- AMREX_ALWAYS_ASSERT(nghost < nghost_subset);
  AMREX_ALWAYS_ASSERT(fba.ixType() == cba.ixType());
 
  // Vector to hold times for coarse data
  m_crse_times.resize(2);
+ m_cf_crse_data.resize(2);
+
+ // Init MF patches
+ m_cf_fine_data = nullptr; m_cf_fine_subset_data = nullptr;
+ m_cf_crse_data[0] = nullptr; m_cf_crse_data[1] = nullptr;
+
+ // Define the coarse and fine MFs
+ Define(fba, fdm, fgeom, cba, cdm, cgeom,
+ nghost, nghost_subset, ncomp, interp);
+}
+
+
+/*
+ * Redefine the coarse and fine patch MultiFabs.
+ *
+ * @param[in] fba BoxArray of data to be filled at fine level
+ * @param[in] fdm DistributionMapping of data to be filled at fine level
+ * @param[in] fgeom container of geometry infomation at fine level
+ * @param[in] cba BoxArray of data to be filled at coarse level
+ * @param[in] cdm DistributionMapping of data to be filled at coarse level
+ * @param[in] cgeom container of geometry infomation at coarse level
+ * @param[in] nghost number of ghost cells to be filled
+ * @param[in] ncomp number of components to be filled
+ * @param[in] interp interpolation operator to be used
+ */
+void ERFFillPatcher::Define (BoxArray const& fba, DistributionMapping fdm,
+ Geometry const& fgeom,
+ BoxArray cba, DistributionMapping cdm,
+ Geometry const& cgeom,
+ int nghost, int nghost_subset,
+ int ncomp, InterpBase* interp)
+{
+ AMREX_ALWAYS_ASSERT(nghost < 0);
+ AMREX_ALWAYS_ASSERT(nghost_subset <= 0);
+ AMREX_ALWAYS_ASSERT(nghost <= nghost_subset);
+
+ // Set data memebers
+ m_fba = fba; m_cba = cba;
+ m_fdm = fdm; m_cdm = cdm;
+ m_fgeom = fgeom; m_cgeom = cgeom;
+ m_nghost = nghost; m_nghost_subset = nghost_subset;
+ m_ncomp = ncomp; m_interp = interp;
+
+ // Delete old MFs if they exist
+ if (m_cf_fine_data) {
+ delete m_cf_fine_data; delete m_cf_fine_subset_data;
+ delete m_cf_crse_data[0]; delete m_cf_crse_data[1];
+ }
 
  // Index type for the BL/BA
  IndexType m_ixt = fba.ixType();
@@ -108,16 +149,16 @@ ERFFillPatcher::ERFFillPatcher (BoxArray const& fba, DistributionMapping fdm,
  DistributionMapping gcf_dm(gcf_fba);
  DistributionMapping cf_dm_s(cf_fba_s);
 
+
  // Fine patch to hold the time-interpolated state
- m_cf_fine_data.define(gcf_fba, gcf_dm, m_ncomp, 0);
+ m_cf_fine_data = new MultiFab (gcf_fba, gcf_dm, m_ncomp, 0);
 
  // Fine subset patch to hold the time-interpolated state
- //m_cf_fine_subset_data.define(cf_fba_s, cf_dm_s, m_ncomp, 0);
- m_cf_fine_subset_data.define(cf_fba_s, gcf_dm, m_ncomp, 0);
+ m_cf_fine_subset_data = new MultiFab (cf_fba_s, gcf_dm, m_ncomp, 0);
 
  // Two coarse patches to hold the data to be interpolated
- m_cf_crse_data.emplace_back(gcf_cba, gcf_dm, m_ncomp, 0);
- m_cf_crse_data.emplace_back(gcf_cba, gcf_dm, m_ncomp, 0);
+ m_cf_crse_data[0] = new MultiFab (gcf_cba, gcf_dm, m_ncomp, 0);
+ m_cf_crse_data[1] = new MultiFab (gcf_cba, gcf_dm, m_ncomp, 0);
 }
 
 
@@ -139,10 +180,10 @@ void ERFFillPatcher::registerCoarseData (Vector<MultiFab const*> const& crse_dat
  // m_cf_crse_data into ghost cells in the z-dir. So we need
  // to include ghost cells for crse_data when doing the copy
  IntVect src_ng = crse_data[0]->nGrowVect();
- IntVect dst_ng = m_cf_crse_data[0].nGrowVect();
- m_cf_crse_data[0].ParallelCopy(*crse_data[0], 0, 0, m_ncomp,
+ IntVect dst_ng = m_cf_crse_data[0]->nGrowVect();
+ m_cf_crse_data[0]->ParallelCopy(*crse_data[0], 0, 0, m_ncomp,
  src_ng, dst_ng, m_cgeom.periodicity()); // old data
- m_cf_crse_data[1].ParallelCopy(*crse_data[1], 0, 0, m_ncomp,
+ m_cf_crse_data[1]->ParallelCopy(*crse_data[1], 0, 0, m_ncomp,
  src_ng, dst_ng, m_cgeom.periodicity()); // new data
 
  m_crse_times[0] = crse_time[0]; // time of "old" coarse data

Source/ERF.H

@@ -306,7 +306,11 @@ private:
 
  void define_grids_to_evolve (int lev, const amrex::BoxArray& ba); // NOLINT
 
- void Define_ERFFillPatchers ();
+ void Construct_ERFFillPatchers (int lev);
+
+ void Define_ERFFillPatchers (int lev);
+
+ void Register_ERFFillPatchers (int lev);
 
  void init1DArrays ();
 
@@ -508,6 +512,8 @@ private:
 #endif
 
  // Fillpatcher classes for coarse-fine boundaries
+ int cf_width{0};
+ int cf_set_width{0};
  amrex::Vector<ERFFillPatcher> FPr_c;
  amrex::Vector<ERFFillPatcher> FPr_u;
  amrex::Vector<ERFFillPatcher> FPr_v;