Update for ref

Docs/sphinx_doc/theory/DryEquations.rst

@@ -114,7 +114,7 @@ the velocity field to create a provisional velocity, :math:`\mathbf{u}^*` at the
 then impose the constraint by solving the pressure Poisson equation for :math:`p^\prime`
 
 .. math::
-   \nabla \cdot (\frac{\overline{\rho}}{\rho} \nabla p^\prime = \nabla \cdot \mathbf{u}^*
+   \nabla \cdot ( \frac{\overline{\rho}}{\rho} \nabla p^\prime ) = \nabla \cdot ( \overline{\rho} \mathbf{u}^* )
 
 then setting
 

Source/ERF_make_new_arrays.cpp

@@ -139,8 +139,10 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
     lev_new[Vars::yvel].define(convert(ba, IntVect(0,1,0)), dm, 1, ngrow_vels);
     lev_old[Vars::yvel].define(convert(ba, IntVect(0,1,0)), dm, 1, ngrow_vels);
 
-    lev_new[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, IntVect(ngrow_vels,ngrow_vels,0));
-    lev_old[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, IntVect(ngrow_vels,ngrow_vels,0));
+    // Note that we need the ghost cells in the z-direction if we are doing any
+    // kind of domain decomposition in the vertical (at level 0 or above)
+    lev_new[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow_vels);
+    lev_old[Vars::zvel].define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow_vels);
 
 #ifdef ERF_USE_POISSON_SOLVE
     pp_inc[lev].define(ba, dm, 1, 1);
@@ -463,24 +465,29 @@ ERF::update_terrain_arrays (int lev, Real time)
         }
 
         //
+        // NOTE: we are not currently doing this as described -- we will simply use
+        //       the interpolation from coarse done above
         // Then, if not using real/metgrid data,
         // 1) redefine the terrain at k=0 for every fine box which includes k=0
         // 2) recreate z_phys_nd at every fine node using
         // the data at the bottom of each fine grid
         // which has been either been interpolated from the coarse grid (k>0)
         // or set in init_custom_terrain (k=0)
         //
-        if (init_type != "real" && init_type != "metgrid") {
-            prob->init_custom_terrain(geom[lev],*z_phys_nd[lev],time);
-
-            Vector<Real> zmax(1); // only reduce at k==0
-            reduce_to_max_per_level(zmax, z_phys_nd[lev]);
-            amrex::Print() << "Max terrain elevation = " << zmax[0] << std::endl;
-            AMREX_ALWAYS_ASSERT_WITH_MESSAGE(zlevels_stag[zlevels_stag.size()-1] > zmax[0],
-                "Terrain is taller than domain top!");
-
-            init_terrain_grid(lev,geom[lev],*z_phys_nd[lev],zlevels_stag,phys_bc_type);
-        }
+        if (lev == 0) {
+            if (init_type != "real" && init_type != "metgrid")
+            {
+                prob->init_custom_terrain(geom[lev],*z_phys_nd[lev],time);
+
+                Vector<Real> zmax(1); // only reduce at lev==0
+                reduce_to_max_per_level(zmax, z_phys_nd[lev]);
+                amrex::Print() << "Max terrain elevation = " << zmax[0] << std::endl;
+                AMREX_ALWAYS_ASSERT_WITH_MESSAGE(zlevels_stag[zlevels_stag.size()-1] > zmax[0],
+                    "Terrain is taller than domain top!");
+
+                init_terrain_grid(lev,geom[lev],*z_phys_nd[lev],zlevels_stag,phys_bc_type);
+            } // init_type
+        } // lev == 0
 
         make_J(geom[lev],*z_phys_nd[lev],*detJ_cc[lev]);
         make_areas(geom[lev],*z_phys_nd[lev],*ax[lev],*ay[lev],*az[lev]);

Source/IO/NCPlotFile.cpp

@@ -31,7 +31,7 @@ ERF::writeNCPlotFile (int lev, int which_subdomain, const std::string& dir,
      // total number of cells in this "domain" at this level
      std::vector<int> n_cells;
 
-     // number of points in each blocks at this level
+     // number of points in each block at this level
      std::vector<int> offset;
 
      // set the full IO path for NetCDF output
@@ -58,7 +58,7 @@ ERF::writeNCPlotFile (int lev, int which_subdomain, const std::string& dir,
      }
      for(auto ib=0; ib<nblocks; ib++) {
         auto npts_per_block = grids[lev][ib].length(0)*grids[lev][ib].length(1)*grids[lev][ib].length(2);
-        offset[dm[ib]] = offset[dm[ib]] + npts_per_block;
+        offset[dm[ib]] += npts_per_block;
      }
 
      // We only do single-level writes when using NetCDF format
@@ -232,17 +232,17 @@ ERF::writeNCPlotFile (int lev, int which_subdomain, const std::string& dir,
    for (MFIter fai(*plotMF[lev]); fai.isValid(); ++fai) {
        auto box = fai.validbox();
        if (subdomain.contains(box)) {
+           numpts = box.numPts();
            long unsigned diff = nfai*numpts;
            for(auto ip = 1; ip <= iproc; ++ip) diff += offset[ip-1];
-           numpts = box.numPts();
 
            for (int k(0); k < ncomp; ++k) {
               const auto *data = plotMF[lev]->get(fai).dataPtr(k);
               auto nc_plot_var = ncf.var(plot_var_names[k]);
               nc_plot_var.par_access(NC_INDEPENDENT);
               nc_plot_var.put(data, {diff}, {numpts});
-          }
-          nfai++;
+           }
+           nfai++;
        }
    }
    ncf.close();

Source/TimeIntegration/ERF_make_tau_terms.cpp

@@ -73,7 +73,7 @@ void erf_make_tau_terms (int level, int nrk,
     std::unique_ptr<MultiFab> dflux_z;
 
     if (l_use_diff) {
-        expr    = std::make_unique<MultiFab>(ba  , dm, 1, IntVect(1,1,0));
+        expr    = std::make_unique<MultiFab>(ba  , dm, 1, IntVect(1,1,1));
         dflux_x = std::make_unique<MultiFab>(convert(ba,IntVect(1,0,0)), dm, nvars, 0);
         dflux_y = std::make_unique<MultiFab>(convert(ba,IntVect(0,1,0)), dm, nvars, 0);
         dflux_z = std::make_unique<MultiFab>(convert(ba,IntVect(0,0,1)), dm, nvars, 0);