Merge branch 'next' into more-zstart-zend-loops

Author: ZedThree

examples/dalf3/dalf3.cxx

@@ -71,8 +71,6 @@ class DALF3 : public PhysicsModel {
 
   BoutReal viscosity, hyper_viscosity;
 
-  bool smooth_separatrix;
-
   FieldGroup comms;
 
   std::unique_ptr<Laplacian> phiSolver{nullptr};  // Laplacian solver in X-Z
@@ -132,7 +130,6 @@ class DALF3 : public PhysicsModel {
     viscosity = options["viscosity"].withDefault(-1.0);
     hyper_viscosity = options["hyper_viscosity"].withDefault(-1.0);
     viscosity_par = options["viscosity_par"].withDefault(-1.0);
-    smooth_separatrix = options["smooth_separatrix"].withDefault(false);
 
     filter_z = options["filter_z"].withDefault(false);
 
@@ -466,11 +463,6 @@ class DALF3 : public PhysicsModel {
     ddt(Pe) = -bracket(phi, Pet, bm)
               + Pet * (Kappa(phi - Pe) + B0 * Grad_parP(jpar - Vpar) / B0);
 
-    if (smooth_separatrix) {
-      // Experimental smoothing across separatrix
-      ddt(Vort) += mesh->smoothSeparatrix(Vort);
-    }
-
     if (filter_z) {
       ddt(Pe) = filter(ddt(Pe), 1);
     }

examples/performance/bracket/bracket.cxx

@@ -68,9 +68,6 @@ int main(int argc, char** argv) {
     ITERATOR_TEST_BLOCK("Bracket [2D,3D] ARAKAWA",
                         result = bracket(a, c, BRACKET_ARAKAWA););
 
-    ITERATOR_TEST_BLOCK("Bracket [2D,3D] ARAKAWA_OLD",
-                        result = bracket(a, c, BRACKET_ARAKAWA_OLD););
-
     ITERATOR_TEST_BLOCK("Bracket [2D,3D] SIMPLE",
                         result = bracket(a, c, BRACKET_SIMPLE););
 
@@ -81,21 +78,12 @@ int main(int argc, char** argv) {
     ITERATOR_TEST_BLOCK("Bracket [3D,3D] ARAKAWA",
                         result = bracket(a, b, BRACKET_ARAKAWA););
 
-    ITERATOR_TEST_BLOCK("Bracket [3D,3D] ARAKAWA_OLD",
-                        result = bracket(a, b, BRACKET_ARAKAWA_OLD););
-
     ITERATOR_TEST_BLOCK("Bracket [3D,3D] SIMPLE",
                         result = bracket(a, b, BRACKET_SIMPLE););
 
     ITERATOR_TEST_BLOCK("Bracket [3D,3D] DEFAULT", result = bracket(a, b, BRACKET_STD););
   }
 
-  // Uncomment below for a "correctness" check
-  // Field3D resNew = bracket(a, b, BRACKET_ARAKAWA); mesh->communicate(resNew);
-  // Field3D resOld = bracket(a, b, BRACKET_ARAKAWA_OLD); mesh->communicate(resOld);
-  // time_output << "Max abs diff is
-  // "<<max(abs(resNew-resOld),true)/max(abs(resOld),true)<<std::endl;
-
   if (profileMode) {
     int nthreads = 0;
 #if BOUT_USE_OPENMP

include/bout/difops.hxx

@@ -271,18 +271,16 @@ Field3D b0xGrad_dot_Grad(const Field3D& phi, const Field3D& A,
  * Poisson bracket methods
  */
 enum class BRACKET_METHOD {
-  standard,   ///< Use b0xGrad_dot_Grad
-  simple,     ///< Keep only terms in X-Z
-  arakawa,    ///< Arakawa method in X-Z (optimised)
-  ctu,        ///< Corner Transport Upwind (CTU) method. Explicit method only, needs the
-              ///  timestep from the solver
-  arakawa_old ///< Older version, for regression testing of optimised version.
+  standard, ///< Use b0xGrad_dot_Grad
+  simple,   ///< Keep only terms in X-Z
+  arakawa,  ///< Arakawa method in X-Z
+  ctu,      ///< Corner Transport Upwind (CTU) method. Explicit method only, needs the
+            ///  timestep from the solver
 };
 constexpr BRACKET_METHOD BRACKET_STD = BRACKET_METHOD::standard;
 constexpr BRACKET_METHOD BRACKET_SIMPLE = BRACKET_METHOD::simple;
 constexpr BRACKET_METHOD BRACKET_ARAKAWA = BRACKET_METHOD::arakawa;
 constexpr BRACKET_METHOD BRACKET_CTU = BRACKET_METHOD::ctu;
-constexpr BRACKET_METHOD BRACKET_ARAKAWA_OLD = BRACKET_METHOD::arakawa_old;
 
 /*!
  * Compute advection operator terms, which can be cast as

include/bout/mesh.hxx

@@ -505,9 +505,6 @@ public:
   virtual void addBoundaryPar(std::shared_ptr<BoundaryRegionPar> UNUSED(bndry),
                               BoundaryParType UNUSED(type)) {}
 
-  /// Branch-cut special handling (experimental)
-  virtual Field3D smoothSeparatrix(const Field3D& f) { return f; }
-
   virtual BoutReal GlobalX(int jx) const = 0;      ///< Continuous X index between 0 and 1
   virtual BoutReal GlobalY(int jy) const = 0;      ///< Continuous Y index (0 -> 1)
   virtual BoutReal GlobalX(BoutReal jx) const = 0; ///< Continuous X index between 0 and 1

manual/sphinx/user_docs/bout_options.rst

@@ -950,7 +950,7 @@ automatically set the ``"time_dimension"`` attribute::
   data["field"] = Field3D(3.0);
   
   // Append data to file
-  bout::OptionsIO({{"file", "time.nc"}, {"append", true}})->write(data);
+  bout::OptionsIO::create({{"file", "time.nc"}, {"append", true}})->write(data);
 
 .. note:: By default, `bout::OptionsIO::write` will only write variables
           with a ``"time_dimension"`` of ``"t"``. You can write

manual/sphinx/user_docs/python_boutpp.rst

@@ -160,7 +160,7 @@ A real example - check derivative contributions:
               phi[-1,:,z]=phi_arr
       with open(path+"/equilibrium/phi_eq.dat","rb") as inf:
           phi_arr=np.fromfile(inf,dtype=np.double)
-          bm="BRACKET_ARAKAWA_OLD"
+          bm="BRACKET_ARAKAWA"
 
           for tind in range(len(times)):
               vort     = Field3D.fromCollect("vort"     ,path=path,tind=tind,info=False)

src/mesh/difops.cxx

@@ -771,46 +771,6 @@ Field3D bracket(const Field3D& f, const Field2D& g, BRACKET_METHOD method,
 
     break;
   }
-  case BRACKET_ARAKAWA_OLD: {
-#if not(BOUT_USE_METRIC_3D)
-    const int ncz = mesh->LocalNz;
-    BOUT_OMP_PERF(parallel for)
-    for (int jx = mesh->xstart; jx <= mesh->xend; jx++) {
-      for (int jy = mesh->ystart; jy <= mesh->yend; jy++) {
-        const BoutReal partialFactor = 1.0 / (12 * metric->dz(jx, jy));
-        const BoutReal spacingFactor = partialFactor / metric->dx(jx, jy);
-        for (int jz = mesh->zstart; jz <= mesh->zend; jz++) {
-          const int jzp = jz + 1 < ncz ? jz + 1 : 0;
-          // Above is alternative to const int jzp = (jz + 1) % ncz;
-          const int jzm = jz - 1 >= 0 ? jz - 1 : ncz - 1;
-          // Above is alternative to const int jzmTmp = (jz - 1 + ncz) % ncz;
-
-          // J++ = DDZ(f)*DDX(g) - DDX(f)*DDZ(g)
-          BoutReal Jpp =
-              ((f(jx, jy, jzp) - f(jx, jy, jzm)) * (g(jx + 1, jy) - g(jx - 1, jy))
-               - (f(jx + 1, jy, jz) - f(jx - 1, jy, jz)) * (g(jx, jy) - g(jx, jy)));
-
-          // J+x
-          BoutReal Jpx = (g(jx + 1, jy) * (f(jx + 1, jy, jzp) - f(jx + 1, jy, jzm))
-                          - g(jx - 1, jy) * (f(jx - 1, jy, jzp) - f(jx - 1, jy, jzm))
-                          - g(jx, jy) * (f(jx + 1, jy, jzp) - f(jx - 1, jy, jzp))
-                          + g(jx, jy) * (f(jx + 1, jy, jzm) - f(jx - 1, jy, jzm)));
-
-          // Jx+
-          BoutReal Jxp = (g(jx + 1, jy) * (f(jx, jy, jzp) - f(jx + 1, jy, jz))
-                          - g(jx - 1, jy) * (f(jx - 1, jy, jz) - f(jx, jy, jzm))
-                          - g(jx - 1, jy) * (f(jx, jy, jzp) - f(jx - 1, jy, jz))
-                          + g(jx + 1, jy) * (f(jx + 1, jy, jz) - f(jx, jy, jzm)));
-
-          result(jx, jy, jz) = (Jpp + Jpx + Jxp) * spacingFactor;
-        }
-      }
-    }
-#else
-    throw BoutException("BRACKET_ARAKAWA_OLD not valid with 3D metrics yet.");
-#endif
-    break;
-  }
   case BRACKET_SIMPLE: {
     // Use a subset of terms for comparison to BOUT-06
     result = VDDX(DDZ(f, outloc), g, outloc);
@@ -1090,63 +1050,6 @@ Field3D bracket(const Field3D& f, const Field3D& g, BRACKET_METHOD method,
     }
     break;
   }
-  case BRACKET_ARAKAWA_OLD: {
-    // Arakawa scheme for perpendicular flow
-
-    const int ncz = mesh->LocalNz;
-
-    // We need to discard const qualifier in order to manipulate
-    // storage array directly
-    Field3D f_temp = f;
-    Field3D g_temp = g;
-
-    BOUT_OMP_PERF(parallel for)
-    for (int jx = mesh->xstart; jx <= mesh->xend; jx++) {
-      for (int jy = mesh->ystart; jy <= mesh->yend; jy++) {
-#if not(BOUT_USE_METRIC_3D)
-        const BoutReal spacingFactor =
-            1.0 / (12 * metric->dz(jx, jy) * metric->dx(jx, jy));
-#endif
-        const BoutReal* Fxm = f_temp(jx - 1, jy);
-        const BoutReal* Fx = f_temp(jx, jy);
-        const BoutReal* Fxp = f_temp(jx + 1, jy);
-        const BoutReal* Gxm = g_temp(jx - 1, jy);
-        const BoutReal* Gx = g_temp(jx, jy);
-        const BoutReal* Gxp = g_temp(jx + 1, jy);
-        for (int jz = 0; jz < mesh->LocalNz; jz++) {
-#if BOUT_USE_METRIC_3D
-          const BoutReal spacingFactor =
-              1.0 / (12 * metric->dz(jx, jy, jz) * metric->dx(jx, jy, jz));
-#endif
-          const int jzp = jz + 1 < ncz ? jz + 1 : 0;
-          // Above is alternative to const int jzp = (jz + 1) % ncz;
-          const int jzm = jz - 1 >= 0 ? jz - 1 : ncz - 1;
-          // Above is alternative to const int jzm = (jz - 1 + ncz) % ncz;
-
-          // J++ = DDZ(f)*DDX(g) - DDX(f)*DDZ(g)
-          // NOLINTNEXTLINE
-          BoutReal Jpp = ((Fx[jzp] - Fx[jzm]) * (Gxp[jz] - Gxm[jz])
-                          - (Fxp[jz] - Fxm[jz]) * (Gx[jzp] - Gx[jzm]));
-
-          // J+x
-          // NOLINTNEXTLINE
-          BoutReal Jpx =
-              (Gxp[jz] * (Fxp[jzp] - Fxp[jzm]) - Gxm[jz] * (Fxm[jzp] - Fxm[jzm])
-               - Gx[jzp] * (Fxp[jzp] - Fxm[jzp]) + Gx[jzm] * (Fxp[jzm] - Fxm[jzm]));
-
-          // Jx+
-          // NOLINTNEXTLINE
-          BoutReal Jxp =
-              (Gxp[jzp] * (Fx[jzp] - Fxp[jz]) - Gxm[jzm] * (Fxm[jz] - Fx[jzm])
-               - Gxm[jzp] * (Fx[jzp] - Fxm[jz]) + Gxp[jzm] * (Fxp[jz] - Fx[jzm]));
-
-          result(jx, jy, jz) = (Jpp + Jpx + Jxp) * spacingFactor;
-        }
-      }
-    }
-
-    break;
-  }
   case BRACKET_SIMPLE: {
     // Use a subset of terms for comparison to BOUT-06
     result = VDDX(DDZ(f, outloc), g, outloc) + VDDZ(-DDX(f, outloc), g, outloc);

src/mesh/impls/bout/boutmesh.cxx

@@ -3115,47 +3115,6 @@ void BoutMesh::addBoundaryPar(std::shared_ptr<BoundaryRegionPar> bndry,
   par_boundary[static_cast<int>(BoundaryParType::all)].push_back(bndry);
 }
 
-Field3D BoutMesh::smoothSeparatrix(const Field3D& f) {
-  Field3D result{emptyFrom(f)};
-  if ((ixseps_inner > 0) && (ixseps_inner < nx - 1)) {
-    if (XPROC(ixseps_inner) == PE_XIND) {
-      int x = getLocalXIndex(ixseps_inner);
-      for (int y = 0; y < LocalNy; y++) {
-        for (int z = 0; z < LocalNz; z++) {
-          result(x, y, z) = 0.5 * (f(x, y, z) + f(x - 1, y, z));
-        }
-      }
-    }
-    if (XPROC(ixseps_inner - 1) == PE_XIND) {
-      int x = getLocalXIndex(ixseps_inner - 1);
-      for (int y = 0; y < LocalNy; y++) {
-        for (int z = 0; z < LocalNz; z++) {
-          result(x, y, z) = 0.5 * (f(x, y, z) + f(x + 1, y, z));
-        }
-      }
-    }
-  }
-  if ((ixseps_outer > 0) && (ixseps_outer < nx - 1) && (ixseps_outer != ixseps_inner)) {
-    if (XPROC(ixseps_outer) == PE_XIND) {
-      int x = getLocalXIndex(ixseps_outer);
-      for (int y = 0; y < LocalNy; y++) {
-        for (int z = 0; z < LocalNz; z++) {
-          result(x, y, z) = 0.5 * (f(x, y, z) + f(x - 1, y, z));
-        }
-      }
-    }
-    if (XPROC(ixseps_outer - 1) == PE_XIND) {
-      int x = getLocalXIndex(ixseps_outer - 1);
-      for (int y = 0; y < LocalNy; y++) {
-        for (int z = 0; z < LocalNz; z++) {
-          result(x, y, z) = 0.5 * (f(x, y, z) + f(x + 1, y, z));
-        }
-      }
-    }
-  }
-  return result;
-}
-
 BoutReal BoutMesh::GlobalX(int jx) const {
   if (symmetricGlobalX) {
     // With this definition the boundary sits dx/2 away form the first/last inner points

src/mesh/impls/bout/boutmesh.hxx

@@ -167,8 +167,6 @@ public:
                       BoundaryParType type) override;
   std::set<std::string> getPossibleBoundaries() const override;
 
-  Field3D smoothSeparatrix(const Field3D& f) override;
-
   int getNx() const { return nx; }
   int getNy() const { return ny; }
 

tests/MMS/spatial/advection/runtest

@@ -43,7 +43,6 @@ build_and_log("MMS steady-state advection test")
 # List of options to be passed for each test
 options = [
     ("method=2", "Arakawa", "-^", 2),
-    ("method=4", "Arakawa-old", "-.", 2),
     (
         "method=1 mesh:ddx:upwind=u1 mesh:ddz:upwind=u1",
         "SIMPLE: 1st order upwind",