Make mixed second-derivative operators more correct

include/derivs.hxx

@@ -694,11 +694,15 @@ inline const Field2D FDDZ(const Field2D& v, const Field2D& f, CELL_LOC outloc,
 ///                    If not given, defaults to DIFF_DEFAULT
 /// @param[in] region  What region is expected to be calculated
 ///                    If not given, defaults to RGN_NOBNDRY
+/// @param[in] dfdy_boundary_condition Boundary condition to use to set the guard cells of
+///                                    df/dy, before calculating the x-derivative.
 const Field3D D2DXDY(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT,
-                     const std::string& method = "DEFAULT", REGION region = RGN_NOBNDRY);
+                     const std::string& method = "DEFAULT", REGION region = RGN_NOBNDRY,
+                     const std::string& dfdy_boundary_condition = "free_o3");
 inline const Field3D D2DXDY(const Field3D& f, CELL_LOC outloc, DIFF_METHOD method,
-                            REGION region = RGN_NOBNDRY) {
-  return D2DXDY(f, outloc, toString(method), region);
+                            REGION region = RGN_NOBNDRY,
+                            const std::string& dfdy_boundary_condition = "free_o3") {
+  return D2DXDY(f, outloc, toString(method), region, dfdy_boundary_condition);
 };
 
 /// Calculate mixed partial derivative in x and y
@@ -713,11 +717,15 @@ inline const Field3D D2DXDY(const Field3D& f, CELL_LOC outloc, DIFF_METHOD metho
 ///                    If not given, defaults to DIFF_DEFAULT
 /// @param[in] region  What region is expected to be calculated
 ///                    If not given, defaults to RGN_NOBNDRY
+/// @param[in] dfdy_boundary_condition Boundary condition to use to set the guard cells of
+///                                    df/dy, before calculating the x-derivative.
 const Field2D D2DXDY(const Field2D& f, CELL_LOC outloc = CELL_DEFAULT,
-                     const std::string& method = "DEFAULT", REGION region = RGN_NOBNDRY);
+                     const std::string& method = "DEFAULT", REGION region = RGN_NOBNDRY,
+                     const std::string& dfdy_boundary_condition = "free_o3");
 inline const Field2D D2DXDY(const Field2D& f, CELL_LOC outloc, DIFF_METHOD method,
-                            REGION region = RGN_NOBNDRY) {
-  return D2DXDY(f, outloc, toString(method), region);
+                            REGION region = RGN_NOBNDRY,
+                            const std::string& dfdy_boundary_condition = "free_o3") {
+  return D2DXDY(f, outloc, toString(method), region, dfdy_boundary_condition);
 };
 
 /// Calculate mixed partial derivative in x and z

manual/sphinx/user_docs/differential_operators.rst

@@ -240,6 +240,84 @@ store using key `"C2"` for all three directions and both fields with
 no staggering.
 
 
+.. _sec-diffmethod-mixedsecond:
+
+Mixed second-derivative operators
+---------------------------------
+
+Coordinate derivatives commute, as long as the coordinates are globally well-defined, i.e.
+
+.. math::
+
+    \frac{\partial}{\partial x} \left(\frac{\partial}{\partial y} f \right)
+    = \frac{\partial}{\partial y} \left(\frac{\partial}{\partial x} f \right) \\
+    \frac{\partial}{\partial y} \left(\frac{\partial}{\partial z} f \right)
+    = \frac{\partial}{\partial z} \left(\frac{\partial}{\partial y} f \right) \\
+    \frac{\partial}{\partial z} \left(\frac{\partial}{\partial x} f \right)
+    = \frac{\partial}{\partial x} \left(\frac{\partial}{\partial z} f \right)
+
+When using ``paralleltransform = shifted`` or ``paralleltransform = fci`` (see
+:ref:`sec-parallel-transforms`) we do not have globally well-defined coordinates. In those
+cases the coordinate systems are field-aligned, but the grid points are at constant
+toroidal angle. The field-aligned coordinates are defined locally, on planes of constant
+:math:`y`. There are different coordinate systems for each plane. However, within each
+local coordinate system the derivatives do commute. :math:`y`-derivatives are taken in the
+local field-aligned coordinate system, so mixed derivatives are calculated as
+
+::
+
+    D2DXDY(f) = DDX(DDY(f))
+    D2DYDZ(f) = DDZ(DDY(f))
+
+This order is simpler -- the alternative is possible. Using second-order central
+difference operators for the y-derivatives we could calculate (not worring about
+communications or boundary conditions here)
+
+::
+
+    Field3D D2DXDY(Field3D f) {
+      auto result{emptyFrom(f)};
+      auto& coords = \*f.getCoordinates()
+
+      auto dfdx_yup = DDX(f.yup());
+      auto dfdx_ydown = DDX(f.ydown());
+
+      BOUT_FOR(i, f.getRegion()) {
+        result[i] = (dfdx_yup[i.yp()] - dfdx_ydown[i.ym()]) / (2. * coords.dy[i])
+      }
+
+      return result;
+    }
+
+This would give equivalent results to the previous form [#]_ as ``yup`` and ``ydown`` give
+the values of ``f`` one grid point along the magnetic field *in the local field-aligned
+coordinate system*.
+
+The :math:`x\mathrm{-}z` derivative is unaffected as it is taken entirely on a plane of
+constant :math:`y` anyway. It is evaluated as
+
+::
+
+    D2DXDZ(f) = DDZ(DDX(f))
+
+As the ``z``-direction is periodic and the ``z``-grid is not split across processors,
+``DDZ`` does not require any guard cells. By taking ``DDZ`` second, we do not have to
+communicate or set boundary conditions on the result of ``DDX`` or ``DDY`` before taking
+``DDZ``.
+
+The derivatives in ``D2DXDY(f)`` are applied in two steps. First ``dfdy = DDY(f)`` is
+calculated; ``dfdy`` is communicated and has a boundary condition applied so that all the
+x-guard cells are filled. The boundary condition is ``free_o3`` by default (3rd order
+extrapolation into the boundary cells), but can be specified with the fifth argument to
+``D2DXDY`` (see :ref:`sec-bndryopts` for possible options). Second ``DDX(dfdy)`` is
+calculated, and returned from the function.
+
+.. [#] Equivalent but not exactly the same numerically. Expanding out the derivatives in
+       second-order central-difference form shows that the two differ in the grid points
+       at which they evaluate ``dx`` and ``dy``. As long as the grid spacings are smooth
+       this should not affect the order of accuracy of the scheme (?).
+
+
 .. _sec-diffmethod-nonuniform:
 
 Non-uniform meshes
@@ -693,7 +771,7 @@ well defined Fourier transform. This means that
 In our case, we are dealing with periodic boundary conditions. Strictly
 speaking, the Fourier transform does not exist in such cases, but it is
 possible to define a Fourier transform in the limit which in the end
-lead to the Fourier series [1]_ By discretising the spatial domain, it
+lead to the Fourier series [#]_ By discretising the spatial domain, it
 is no longer possible to represent the infinite amount of Fourier modes,
 but only :math:`N+1` number of modes, where :math:`N` is the number of
 points (this includes the modes with negative frequencies, and the
@@ -724,5 +802,5 @@ The discrete version of equation (:eq:`f_derivative`) thus gives
 
    \partial_z^n F(x,y)_k = (i k)^n F(x,y)_k
 
-.. [1] For more detail see Bracewell, R. N. - The Fourier Transform
+.. [#] For more detail see Bracewell, R. N. - The Fourier Transform
        and Its Applications 3rd Edition chapter 10

src/sys/derivs.cxx

@@ -276,26 +276,46 @@ const Field2D D4DZ4(const Field2D &f, CELL_LOC outloc, const std::string &method
 /*!
  * Mixed derivative in X and Y
  *
- * This first takes derivatives in X, then in Y.
+ * This first takes derivatives in Y, then in X.
  *
- * ** Applies Neumann boundary in Y, communicates
+ * ** Communicates and applies boundary in X.
  */
-const Field2D D2DXDY(const Field2D &f, CELL_LOC outloc, const std::string &method, REGION region) {
-  Field2D dfdy = DDY(f, outloc, method, RGN_NOY);
+const Field2D D2DXDY(const Field2D& f, CELL_LOC outloc, const std::string& method,
+    REGION region, const std::string& dfdy_boundary_condition) {
+
+  // If staggering in x, take y-derivative at f's location.
+  const auto y_location =
+    (outloc == CELL_XLOW or f.getLocation() == CELL_XLOW) ? CELL_DEFAULT : outloc;
+
+  Field2D dfdy = DDY(f, y_location, method, region);
+
+  // Set x-guard cells and x-boundary cells before calculating DDX
   f.getMesh()->communicate(dfdy);
+  dfdy.applyBoundary(dfdy_boundary_condition);
+
   return DDX(dfdy, outloc, method, region);
 }
 
 /*!
  * Mixed derivative in X and Y
  *
- * This first takes derivatives in X, then in Y.
+ * This first takes derivatives in Y, then in X.
  *
- * ** Applies Neumann boundary in Y, communicates
+ * ** Communicates and applies boundary in X.
  */
-const Field3D D2DXDY(const Field3D &f, CELL_LOC outloc, const std::string &method, REGION region) {
-  Field3D dfdy = DDY(f, outloc, method, RGN_NOY);
+const Field3D D2DXDY(const Field3D& f, CELL_LOC outloc, const std::string& method,
+    REGION region, const std::string& dfdy_boundary_condition) {
+
+  // If staggering in x, take y-derivative at f's location.
+  const auto y_location =
+    (outloc == CELL_XLOW or f.getLocation() == CELL_XLOW) ? CELL_DEFAULT : outloc;
+
+  Field3D dfdy = DDY(f, y_location, method, region);
+
+  // Set x-guard cells and x-boundary cells before calculating DDX
   f.getMesh()->communicate(dfdy);
+  dfdy.applyBoundary(dfdy_boundary_condition);
+
   return DDX(dfdy, outloc, method, region);
 }
 
@@ -309,26 +329,11 @@ const Field2D D2DXDZ(const Field2D &f, CELL_LOC outloc,
 
 /// X-Z mixed derivative
 const Field3D D2DXDZ(const Field3D &f, CELL_LOC outloc, const std::string &method, REGION region) {
-  // Take derivative in Z, including in X boundaries. Then take derivative in X
-  // Maybe should average results of DDX(DDZ) and DDZ(DDX)?
-  ASSERT1(outloc == CELL_DEFAULT || outloc == f.getLocation());
-  // region specifies what the combined derivative should return
-  // Therefore we need to add the X boundary to the inner derivative
-  // RGN_NOY and RGN_NOZ include the X boundary, therefore we need to
-  // throw - or add communication code.
-  REGION region_inner;
-  switch (region){
-  case RGN_NOBNDRY:
-    region_inner = RGN_NOY;
-    break;
-  case RGN_NOX:
-    region_inner = RGN_ALL;
-    break;
-  default:
-    throw BoutException("Unhandled region case in D2DXDZ");
-  }
+  // If staggering in z, take x-derivative at f's location.
+  const auto x_location =
+    (outloc == CELL_ZLOW or f.getLocation() == CELL_ZLOW) ? CELL_DEFAULT : outloc;
 
-  return DDX(DDZ(f, outloc,method, region_inner),outloc,method,region);;
+  return DDZ(DDX(f, x_location, method, region), outloc, method, region);
 }
 
 const Field2D D2DYDZ(const Field2D &f, CELL_LOC outloc,
@@ -339,32 +344,13 @@ const Field2D D2DYDZ(const Field2D &f, CELL_LOC outloc,
   return zeroFrom(f).setLocation(outloc);
 }
 
-const Field3D D2DYDZ(const Field3D& f, CELL_LOC outloc,
-                     MAYBE_UNUSED(const std::string& method), REGION UNUSED(region)) {
-  Coordinates *coords = f.getCoordinates(outloc);
+const Field3D D2DYDZ(const Field3D& f, CELL_LOC outloc, const std::string& method,
+                     REGION region) {
+  // If staggering in z, take y-derivative at f's location.
+  const auto y_location =
+    (outloc == CELL_ZLOW or f.getLocation() == CELL_ZLOW) ? CELL_DEFAULT : outloc;
 
-  Field3D result{emptyFrom(f)};
-  ASSERT1(outloc == CELL_DEFAULT || outloc == f.getLocation());
-  result.allocate();
-  result.setLocation(f.getLocation());
-  ASSERT1(method == "DEFAULT");
-  for(int i=f.getMesh()->xstart;i<=f.getMesh()->xend;i++)
-    for(int j=f.getMesh()->ystart;j<=f.getMesh()->yend;j++)
-      for(int k=0;k<f.getMesh()->LocalNz;k++) {
-        int kp = (k+1) % (f.getMesh()->LocalNz);
-        int km = (k-1+f.getMesh()->LocalNz) % (f.getMesh()->LocalNz);
-        result(i,j,k) = 0.25*( +(f(i,j+1,kp) - f(i,j-1,kp))
-                               -(f(i,j+1,km) - f(i,j-1,km)) )
-                    / (coords->dy(i,j) * coords->dz);
-      }
-  // TODO: use region aware implementation
-  // BOUT_FOR(i, f.getRegion(region)) {
-  // result[i] = 0.25*( +(f[i.offset(0,1, 1)] - f[i.offset(0,-1, 1)])
-  //                              / (coords->dy[i.yp()])
-  //                    -(f[i.offset(0,1,-1)] - f[i.offset(0,-1,-1)])
-  //                              / (coords->dy[i.ym()]))
-  //   / coords->dz; }
-  return result;
+  return DDZ(DDY(f, y_location, method, region), outloc, method, region);
 }
 
 /*******************************************************************************