fix for OMP plus some cleanup (#1294)

Docs/sphinx_doc/MeshRefinement.rst

@@ -7,7 +7,8 @@
 Mesh Refinement
 ===============
 
-ERF allows both static and dynamic mesh refinement.
+ERF allows both static and dynamic mesh refinement, as well as three different options for
+how the coarse and fine data.
 
 Note that any tagged region will be covered by one or more boxes.  The user may
 specify the refinement criteria and/or region to be covered, but not the decomposition of the region into
@@ -123,17 +124,18 @@ computed by dividing the variable named rhotheta by the variable named density.
 Coupling Types
 --------------
 
-ERF supports both one-way and two-coupling; this is a run-time input
+ERF supports one-way, two-way, and "mixed" coupling; this is a run-time input
 
 ::
 
-      erf.coupling_type = "OneWay" or "TwoWay"
+      erf.coupling_type = "OneWay" or "TwoWay" or "Mixed"
 
 By one-way coupling, we mean that between each pair of refinement levels,
-only the coarse mesh communicates data to the fine mesh. For cell-centered quantities,
+the coarse mesh communicates data to the fine mesh to serve as boundary conditions
+for the time advance of the fine solution . For cell-centered quantities,
 and face-baced normal momenta on the coarse-fine interface, the coarse data is conservatively
 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
+(outside of the valid fine region) as well as specified and relaxation 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.cf_width = <Int>`` (yellow + blue)
 and analogously the width of the interior Dirichlet region may be specified with
@@ -167,9 +169,9 @@ the RHS (:math:`F`) is given by the following:
    H_{2} &= \frac{1}{50 \Delta t} \frac{{\rm SpecWidth} + {\rm RelaxWidth} - n}{{\rm RelaxWidth} - 1},
    \end{align}
 
-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. The set and
+where :math:`G` is the RHS of the evolution equations, :math:`\psi^{\prime}` is the predicted update without
+relaxation, :math:`\psi^{FP}` is the fine data obtained from spatial and temporal interpolation of the
+coarse data, and :math:`n` is the minimum number of grid points from a lateral boundary. The specified and
 relaxation regions are applied to all dycore variables :math:`\left[\rho \; \rho\Theta \; U\; V\; W \right]`
 on the fine mesh. Finally, we note that time dependent Dirichlet data, provided via an external boundary file,
 may be enforced on the lateral boundary conditions of the domain (coarsest mesh). For such cases,
@@ -178,15 +180,14 @@ the relaxation region width at the domain edges may be specified with ``erf.wrfb
 (yellow). With the boundary file approach, all dycore variables are set and relaxed but
 moisture is only set in the yellow region if it is present within the boundary file.
 
-
 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:
 
 - The fine cell-centered data is conservatively averaged onto the coarse mesh covered by fine mesh.
 
 - A "reflux" operation is performed for all cell-centered data.
 
-Because the normal momentum at the fine level is itself interpolated from the coarse, the
+Because the normal momentum at the fine level is itself interpolated from the coarse level, the
 difference between fluxes -- along the coarse-fine interfaces -- used to update the coarse data and fluxes
 used to update the fine data is due to the difference in the averaging of the advected quantity to the face
 where the flux is defined.
@@ -196,3 +197,7 @@ equation are the momenta themselves, which are interpolated on faces at the coar
 quantities which are advanced in conservation form will lose conservation with one-way coupling.
 Two-way coupling is conservative for these scalars as long as the refluxing operation is included with the
 averaging down.
+
+We define "mixed" coupling as using the two-way coupling algorithm for all cell-centered quantities except for
+:math:`\rho` and :math:`\rho \theta.`
+

Source/Advection/Advection.H

@@ -11,7 +11,7 @@
 
 
 /** Compute advection tendency for density and potential temperature */
-void AdvectionSrcForRho (const amrex::Box& bx, const amrex::Box& valid_bx,
+void AdvectionSrcForRho (const amrex::Box& bx,
                          const amrex::Array4<amrex::Real>& src,
                          const amrex::Array4<const amrex::Real>& rho_u,    // These are being used
                          const amrex::Array4<const amrex::Real>& rho_v,    //  to define the fluxes

Source/Advection/AdvectionSrcForState.cpp

@@ -12,7 +12,6 @@ using namespace amrex;
  * functions when either (or both) spatial order(s) is greater than 2.
  *
  * @param[in] bx box over which the scalars are updated
- * @param[in] valid_bx box that contains only the cells not in the specified or relaxation zones
  * @param[out] advectionSrc tendency for the scalar update equation
  * @param[in] rho_u x-component of momentum
  * @param[in] rho_v y-component of momentum
@@ -32,7 +31,7 @@ using namespace amrex;
  */
 
 void
-AdvectionSrcForRho (const Box& bx, const Box& valid_bx,
+AdvectionSrcForRho (const Box& bx,
                     const Array4<Real>& advectionSrc,
                     const Array4<const Real>& rho_u,
                     const Array4<const Real>& rho_v,
@@ -52,9 +51,6 @@ AdvectionSrcForRho (const Box& bx, const Box& valid_bx,
     BL_PROFILE_VAR("AdvectionSrcForRhoAndTheta", AdvectionSrcForRhoAndTheta);
     auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2];
 
-    // We note that valid_bx is the actual grid, while bx may be a tile within that grid
-    const auto& vbx_hi = ubound(valid_bx);
-
     const Box xbx = surroundingNodes(bx,0);
     const Box ybx = surroundingNodes(bx,1);
     const Box zbx = surroundingNodes(bx,2);

Source/Diffusion/PBLModels.cpp

@@ -58,7 +58,6 @@ ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
 
       const amrex::Box &bx = mfi.growntilebox(1);
       const amrex::Array4<amrex::Real const> &cell_data     = cons_in.array(mfi);
-      const amrex::Array4<amrex::Real const> &cell_data_old = cons_old.array(mfi);
       const amrex::Array4<amrex::Real      > &K_turb = eddyViscosity.array(mfi);
       const amrex::Array4<amrex::Real const> &uvel = xvel.array(mfi);
       const amrex::Array4<amrex::Real const> &vvel = yvel.array(mfi);

Source/TimeIntegration/ERF_slow_rhs_pre.cpp

@@ -464,6 +464,11 @@ void erf_slow_rhs_pre (int level, int finest_level,
     } // l_use_diff
 
 
+    // We must initialize these to zero each RK step
+    S_scratch[IntVar::xmom].setVal(0.);
+    S_scratch[IntVar::ymom].setVal(0.);
+    S_scratch[IntVar::zmom].setVal(0.);
+
     // *************************************************************************
     // Define updates and fluxes in the current RK stage
     // *************************************************************************
@@ -490,11 +495,6 @@ void erf_slow_rhs_pre (int level, int finest_level,
         const Array4<Real>       & cell_rhs   = S_rhs[IntVar::cons].array(mfi);
         const Array4<const Real> & buoyancy_fab = buoyancy.const_array(mfi);
 
-        // We must initialize these to zero each RK step
-        S_scratch[IntVar::xmom][mfi].template setVal<RunOn::Device>(0.,tbx);
-        S_scratch[IntVar::ymom][mfi].template setVal<RunOn::Device>(0.,tby);
-        S_scratch[IntVar::zmom][mfi].template setVal<RunOn::Device>(0.,tbz);
-
         Array4<Real> avg_xmom = S_scratch[IntVar::xmom].array(mfi);
         Array4<Real> avg_ymom = S_scratch[IntVar::ymom].array(mfi);
         Array4<Real> avg_zmom = S_scratch[IntVar::zmom].array(mfi);
@@ -651,7 +651,7 @@ void erf_slow_rhs_pre (int level, int finest_level,
         // **************************************************************************
         // Define updates in the RHS of continuity and potential temperature equations
         // **************************************************************************
-        AdvectionSrcForRho(bx, valid_bx, cell_rhs,
+        AdvectionSrcForRho(bx, cell_rhs,
                            rho_u, rho_v, omega_arr,      // these are being used to build the fluxes
                            avg_xmom, avg_ymom, avg_zmom, // these are being defined from the fluxes
                            z_nd, detJ_arr, dxInv, mf_m, mf_u, mf_v,