Add more tests for FCI operators

CMakeLists.txt

@@ -376,6 +376,7 @@ if (zoidberg_FOUND EQUAL 0)
 else()
   set(zoidberg_FOUND OFF)
 endif()
+message(STATUS "Found Zoidberg for FCI tests: ${zoidberg_FOUND}")
 
 option(BOUT_GENERATE_FIELDOPS "Automatically re-generate the Field arithmetic operators from the Python templates. \
 Requires Python3, clang-format, and Jinja2. Turn this OFF to skip generating them if, for example, \

include/bout/difops.hxx

@@ -40,7 +40,9 @@
 #include "bout/field3d.hxx"
 
 #include "bout/bout_types.hxx"
-#include "bout/solver.hxx"
+#include "bout/coordinates.hxx"
+
+class Solver;
 
 /*!
  * Parallel derivative (central differencing) in Y

src/mesh/coordinates.cxx

@@ -1588,7 +1588,7 @@ Field3D Coordinates::Div_par(const Field3D& f, CELL_LOC outloc,
 
   // Need Bxy at location of f, which might be different from location of this
   // Coordinates object
-  auto Bxy_floc = f.getCoordinates()->Bxy;
+  const auto& Bxy_floc = f.getCoordinates()->Bxy;
 
   if (!f.hasParallelSlices()) {
     // No yup/ydown fields. The Grad_par operator will

src/mesh/interpolation/hermite_spline_xz.cxx

@@ -157,7 +157,7 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z
 
   const int ny = localmesh->LocalNy;
   const int nz = localmesh->LocalNz;
-  const int xend = (localmesh->xend - localmesh->xstart + 1) * localmesh->getNXPE()
+  const int xend = ((localmesh->xend - localmesh->xstart + 1) * localmesh->getNXPE())
                    + localmesh->xstart - 1;
 #ifdef HS_USE_PETSC
   IndConverter conv{localmesh};
@@ -177,7 +177,19 @@ void XZHermiteSpline::calcWeights(const Field3D& delta_x, const Field3D& delta_z
     BoutReal t_x = delta_x(x, y, z) - static_cast<BoutReal>(i_corn);
     BoutReal t_z = delta_z(x, y, z) - static_cast<BoutReal>(k_corner(x, y, z));
 
-    // NOTE: A (small) hack to avoid one-sided differences
+    // NOTE: A (small) hack to avoid one-sided differences. We need at
+    // least 2 interior points due to an awkwardness with the
+    // boundaries. The splines need derivatives in x, but we don't
+    // know the value in the boundaries, so _any_ interpolation in the
+    // last interior cell can't be done. Instead, we fudge the
+    // interpolation in the last cell to be at the extreme right-hand
+    // edge of the previous cell (that is, exactly on the last
+    // interior point). However, this doesn't work with only one
+    // interior point, because we have to do something similar to the
+    // _first_ cell, and these two fudges cancel out and we end up
+    // indexing into the boundary anyway.
+    // TODO(peter): Can we remove this if we apply (dirichlet?) BCs to
+    // the X derivatives? Note that we need at least _2_
     if (i_corn >= xend) {
       i_corn = xend - 1;
       t_x = 1.0;
@@ -346,6 +358,8 @@ Field3D XZHermiteSpline::interpolate(const Field3D& f, const std::string& region
   ASSERT1(f.getMesh() == localmesh);
   Field3D f_interp{emptyFrom(f)};
 
+  // TODO(peter): Should we apply dirichlet BCs to derivatives?
+
 #if USE_NEW_WEIGHTS
 #ifdef HS_USE_PETSC
   BoutReal* ptr;

src/mesh/parallel/fci.cxx

@@ -37,79 +37,99 @@
  **************************************************************************/
 
 #include "fci.hxx"
+
+#include "bout/assert.hxx"
+#include "bout/bout_types.hxx"
+#include "bout/boutexception.hxx"
+#include "bout/field2d.hxx"
+#include "bout/field3d.hxx"
+#include "bout/field_data.hxx"
+#include "bout/mesh.hxx"
+#include "bout/msg_stack.hxx"
+#include "bout/options.hxx"
 #include "bout/parallel_boundary_op.hxx"
 #include "bout/parallel_boundary_region.hxx"
-#include <bout/bout_types.hxx>
-#include <bout/constants.hxx>
-#include <bout/mesh.hxx>
-#include <bout/msg_stack.hxx>
-#include <bout/utils.hxx>
+#include "bout/paralleltransform.hxx"
+#include "bout/region.hxx"
+
+#include <fmt/format.h>
 
+#include <array>
+#include <cmath>
+#include <cstddef>
+#include <cstdlib>
+#include <memory>
 #include <string>
+#include <string_view>
+
+using namespace std::string_view_literals;
 
-FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& UNUSED(dy), Options& options,
-               int offset_, const std::shared_ptr<BoundaryRegionPar>& inner_boundary,
+FCIMap::FCIMap(Mesh& mesh, [[maybe_unused]] const Coordinates::FieldMetric& dy,
+               Options& options, int offset,
+               const std::shared_ptr<BoundaryRegionPar>& inner_boundary,
                const std::shared_ptr<BoundaryRegionPar>& outer_boundary, bool zperiodic)
-    : map_mesh(mesh), offset(offset_),
-      region_no_boundary(map_mesh.getRegion("RGN_NOBNDRY")),
+    : map_mesh(&mesh), offset_(offset),
+      region_no_boundary(map_mesh->getRegion("RGN_NOBNDRY")),
       corner_boundary_mask(map_mesh) {
 
-  TRACE("Creating FCIMAP for direction {:d}", offset);
+  TRACE("Creating FCIMAP for direction {:d}", offset_);
 
-  if (offset == 0) {
+  if (offset_ == 0) {
     throw BoutException(
         "FCIMap called with offset = 0; You probably didn't mean to do that");
   }
 
   auto& interpolation_options = options["xzinterpolation"];
-  interp =
-      XZInterpolationFactory::getInstance().create(&interpolation_options, &map_mesh);
-  interp->setYOffset(offset);
+  interp = XZInterpolationFactory::getInstance().create(&interpolation_options, map_mesh);
+  interp->setYOffset(offset_);
 
   interp_corner =
-      XZInterpolationFactory::getInstance().create(&interpolation_options, &map_mesh);
-  interp_corner->setYOffset(offset);
+      XZInterpolationFactory::getInstance().create(&interpolation_options, map_mesh);
+  interp_corner->setYOffset(offset_);
 
   // Index-space coordinates of forward/backward points
-  Field3D xt_prime{&map_mesh}, zt_prime{&map_mesh};
+  Field3D xt_prime{map_mesh};
+  Field3D zt_prime{map_mesh};
 
   // Real-space coordinates of grid points
-  Field3D R{&map_mesh}, Z{&map_mesh};
+  Field3D R{map_mesh};
+  Field3D Z{map_mesh};
 
   // Real-space coordinates of forward/backward points
-  Field3D R_prime{&map_mesh}, Z_prime{&map_mesh};
+  Field3D R_prime{map_mesh};
+  Field3D Z_prime{map_mesh};
 
-  map_mesh.get(R, "R", 0.0, false);
-  map_mesh.get(Z, "Z", 0.0, false);
+  map_mesh->get(R, "R", 0.0, false);
+  map_mesh->get(Z, "Z", 0.0, false);
 
   // Get a unique name for a field based on the sign/magnitude of the offset
-  const auto parallel_slice_field_name = [&](std::string field) -> std::string {
-    const std::string direction = (offset > 0) ? "forward" : "backward";
+  const auto parallel_slice_field_name = [&](std::string_view field) -> std::string {
+    const auto direction = (offset_ > 0) ? "forward"sv : "backward"sv;
     // We only have a suffix for parallel slices beyond the first
     // This is for backwards compatibility
-    const std::string slice_suffix =
-        (std::abs(offset) > 1) ? "_" + std::to_string(std::abs(offset)) : "";
-    return direction + "_" + field + slice_suffix;
+    const auto slice_suffix =
+        (std::abs(offset_) > 1) ? fmt::format("_{}", std::abs(offset_)) : "";
+    return fmt::format("{}_{}{}", direction, field, slice_suffix);
   };
 
   // If we can't read in any of these fields, things will silently not
   // work, so best throw
-  if (map_mesh.get(xt_prime, parallel_slice_field_name("xt_prime"), 0.0, false) != 0) {
+  if (map_mesh->get(xt_prime, parallel_slice_field_name("xt_prime"), 0.0, false) != 0) {
     throw BoutException("Could not read {:s} from grid file!\n"
                         "  Either add it to the grid file, or reduce MYG",
                         parallel_slice_field_name("xt_prime"));
   }
-  if (map_mesh.get(zt_prime, parallel_slice_field_name("zt_prime"), 0.0, false) != 0) {
+  if (map_mesh->get(zt_prime, parallel_slice_field_name("zt_prime"), 0.0, false) != 0) {
     throw BoutException("Could not read {:s} from grid file!\n"
                         "  Either add it to the grid file, or reduce MYG",
                         parallel_slice_field_name("zt_prime"));
   }
-  if (map_mesh.get(R_prime, parallel_slice_field_name("R"), 0.0, false) != 0) {
+  if (map_mesh->get(R_prime, parallel_slice_field_name("R"), 0.0, false) != 0) {
     throw BoutException("Could not read {:s} from grid file!\n"
                         "  Either add it to the grid file, or reduce MYG",
                         parallel_slice_field_name("R"));
   }
-  if (map_mesh.get(Z_prime, parallel_slice_field_name("Z"), 0.0, false) != 0) {
+  if (map_mesh->get(Z_prime, parallel_slice_field_name("Z"), 0.0, false) != 0) {
     throw BoutException("Could not read {:s} from grid file!\n"
                         "  Either add it to the grid file, or reduce MYG",
                         parallel_slice_field_name("Z"));
@@ -157,25 +177,26 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& UNUSED(dy), Options&
     interp->calcWeights(xt_prime, zt_prime);
   }
 
-  const int ncz = map_mesh.LocalNz;
+  const int ncz = map_mesh->LocalNz;
 
   BoutMask to_remove(map_mesh);
-  const int xend =
-      map_mesh.xstart + (map_mesh.xend - map_mesh.xstart + 1) * map_mesh.getNXPE() - 1;
+  const int xend = map_mesh->xstart
+                   + ((map_mesh->xend - map_mesh->xstart + 1) * map_mesh->getNXPE()) - 1;
   // Serial loop because call to BoundaryRegionPar::addPoint
   // (probably?) can't be done in parallel
   BOUT_FOR_SERIAL(i, xt_prime.getRegion("RGN_NOBNDRY")) {
     // z is periodic, so make sure the z-index wraps around
     if (zperiodic) {
-      zt_prime[i] = zt_prime[i]
-                    - ncz * (static_cast<int>(zt_prime[i] / static_cast<BoutReal>(ncz)));
+      zt_prime[i] =
+          zt_prime[i]
+          - (ncz * (static_cast<int>(zt_prime[i] / static_cast<BoutReal>(ncz))));
 
       if (zt_prime[i] < 0.0) {
         zt_prime[i] += ncz;
       }
     }
 
-    if ((xt_prime[i] >= map_mesh.xstart) and (xt_prime[i] <= xend)) {
+    if ((xt_prime[i] >= map_mesh->xstart) and (xt_prime[i] <= xend)) {
       // Not a boundary
       continue;
     }
@@ -215,7 +236,7 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& UNUSED(dy), Options&
     const BoutReal dR_dz = 0.5 * (R[i_zp] - R[i_zm]);
     const BoutReal dZ_dz = 0.5 * (Z[i_zp] - Z[i_zm]);
 
-    const BoutReal det = dR_dx * dZ_dz - dR_dz * dZ_dx; // Determinant of 2x2 matrix
+    const BoutReal det = (dR_dx * dZ_dz) - (dR_dz * dZ_dx); // Determinant of 2x2 matrix
 
     const BoutReal dR = R_prime[i] - R[i];
     const BoutReal dZ = Z_prime[i] - Z[i];
@@ -228,9 +249,9 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& UNUSED(dy), Options&
     // outer boundary. However, if any of the surrounding points are negative,
     // that also means inner. So to differentiate between inner and outer we
     // need at least 2 points in the domain.
-    ASSERT2(map_mesh.xend - map_mesh.xstart >= 2);
-    auto boundary = (xt_prime[i] < map_mesh.xstart) ? inner_boundary : outer_boundary;
-    boundary->add_point(x, y, z, x + dx, y + 0.5 * offset,
+    ASSERT2(map_mesh->xend - map_mesh->xstart >= 2);
+    auto boundary = (xt_prime[i] < map_mesh->xstart) ? inner_boundary : outer_boundary;
+    boundary->add_point(x, y, z, x + dx, y + (0.5 * offset_),
                         z + dz, // Intersection point in local index space
                         0.5,    // Distance to intersection
                         1       // Default to that there is a point in the other direction
@@ -240,21 +261,22 @@ FCIMap::FCIMap(Mesh& mesh, const Coordinates::FieldMetric& UNUSED(dy), Options&
 
   interp->setRegion(region_no_boundary);
 
-  const auto region = fmt::format("RGN_YPAR_{:+d}", offset);
-  if (not map_mesh.hasRegion3D(region)) {
+  const auto region = fmt::format("RGN_YPAR_{:+d}", offset_);
+  if (not map_mesh->hasRegion3D(region)) {
     // The valid region for this slice
-    map_mesh.addRegion3D(
-        region, Region<Ind3D>(map_mesh.xstart, map_mesh.xend, map_mesh.ystart + offset,
-                              map_mesh.yend + offset, 0, map_mesh.LocalNz - 1,
-                              map_mesh.LocalNy, map_mesh.LocalNz));
+    map_mesh->addRegion3D(region, Region<Ind3D>(map_mesh->xstart, map_mesh->xend,
+                                                map_mesh->ystart + offset_,
+                                                map_mesh->yend + offset_, 0,
+                                                map_mesh->LocalNz - 1, map_mesh->LocalNy,
+                                                map_mesh->LocalNz));
   }
 }
 
 Field3D FCIMap::integrate(Field3D& f) const {
   TRACE("FCIMap::integrate");
 
   ASSERT1(f.getDirectionY() == YDirectionType::Standard);
-  ASSERT1(&map_mesh == f.getMesh());
+  ASSERT1(map_mesh == f.getMesh());
 
   // Cell centre values
   Field3D centre = interp->interpolate(f);
@@ -269,7 +291,7 @@ Field3D FCIMap::integrate(Field3D& f) const {
 #endif
 
   BOUT_FOR(i, region_no_boundary) {
-    const auto inext = i.yp(offset);
+    const auto inext = i.yp(offset_);
     const BoutReal f_c = centre[inext];
     const auto izm = i.zm();
     const int x = i.x();
@@ -278,7 +300,7 @@ Field3D FCIMap::integrate(Field3D& f) const {
     const int zm = izm.z();
     if (corner_boundary_mask(x, y, z) || corner_boundary_mask(x - 1, y, z)
         || corner_boundary_mask(x, y, zm) || corner_boundary_mask(x - 1, y, zm)
-        || (x == map_mesh.xstart)) {
+        || (x == map_mesh->xstart)) {
       // One of the corners leaves the domain.
       // Use the cell centre value, since boundary conditions are not
       // currently applied to corners.
@@ -299,19 +321,69 @@ Field3D FCIMap::integrate(Field3D& f) const {
   return result;
 }
 
+FCITransform::FCITransform(Mesh& mesh, const Coordinates::FieldMetric& dy, bool zperiodic,
+                           Options* opt)
+    : ParallelTransform(mesh, opt), R{&mesh}, Z{&mesh} {
+
+  // check the coordinate system used for the grid data source
+  FCITransform::checkInputGrid();
+
+  // Real-space coordinates of grid cells
+  mesh.get(R, "R", 0.0, false);
+  mesh.get(Z, "Z", 0.0, false);
+
+  auto forward_boundary_xin =
+      std::make_shared<BoundaryRegionPar>("FCI_forward", BNDRY_PAR_FWD_XIN, +1, &mesh);
+  auto backward_boundary_xin =
+      std::make_shared<BoundaryRegionPar>("FCI_backward", BNDRY_PAR_BKWD_XIN, -1, &mesh);
+  auto forward_boundary_xout =
+      std::make_shared<BoundaryRegionPar>("FCI_forward", BNDRY_PAR_FWD_XOUT, +1, &mesh);
+  auto backward_boundary_xout =
+      std::make_shared<BoundaryRegionPar>("FCI_backward", BNDRY_PAR_BKWD_XOUT, -1, &mesh);
+
+  // Add the boundary region to the mesh's vector of parallel boundaries
+  mesh.addBoundaryPar(forward_boundary_xin, BoundaryParType::xin_fwd);
+  mesh.addBoundaryPar(backward_boundary_xin, BoundaryParType::xin_bwd);
+  mesh.addBoundaryPar(forward_boundary_xout, BoundaryParType::xout_fwd);
+  mesh.addBoundaryPar(backward_boundary_xout, BoundaryParType::xout_bwd);
+
+  field_line_maps.reserve(static_cast<std::size_t>(mesh.ystart) * 2);
+  for (int offset = 1; offset < mesh.ystart + 1; ++offset) {
+    field_line_maps.emplace_back(mesh, dy, options, offset, forward_boundary_xin,
+                                 forward_boundary_xout, zperiodic);
+    field_line_maps.emplace_back(mesh, dy, options, -offset, backward_boundary_xin,
+                                 backward_boundary_xout, zperiodic);
+  }
+  ASSERT0(mesh.ystart == 1);
+  const std::array bndries = {forward_boundary_xin, forward_boundary_xout,
+                              backward_boundary_xin, backward_boundary_xout};
+  for (const auto& bndry : bndries) {
+    for (const auto& bndry2 : bndries) {
+      if (bndry->dir == bndry2->dir) {
+        continue;
+      }
+      for (bndry->first(); !bndry->isDone(); bndry->next()) {
+        if (bndry2->contains(*bndry)) {
+          bndry->setValid(0);
+        }
+      }
+    }
+  }
+}
+
 void FCITransform::checkInputGrid() {
   std::string parallel_transform;
   if (mesh.isDataSourceGridFile()
-      && !mesh.get(parallel_transform, "parallel_transform")) {
+      && (mesh.get(parallel_transform, "parallel_transform") == 0)) {
     if (parallel_transform != "fci") {
       throw BoutException(
           "Incorrect parallel transform type '" + parallel_transform
           + "' used "
             "to generate metric components for FCITransform. Should be 'fci'.");
     }
   } // else: parallel_transform variable not found in grid input, indicates older input
-    //       file or grid from options so must rely on the user having ensured the type is
-    //       correct
+  //       file or grid from options so must rely on the user having ensured the type is
+  //       correct
 }
 
 void FCITransform::calcParallelSlices(Field3D& f) {
@@ -327,8 +399,8 @@ void FCITransform::calcParallelSlices(Field3D& f) {
 
   // Interpolate f onto yup and ydown fields
   for (const auto& map : field_line_maps) {
-    f.ynext(map.offset) = map.interpolate(f);
-    f.ynext(map.offset).setRegion(fmt::format("RGN_YPAR_{:+d}", map.offset));
+    f.ynext(map.offset()) = map.interpolate(f);
+    f.ynext(map.offset()).setRegion(fmt::format("RGN_YPAR_{:+d}", map.offset()));
   }
 }
 
@@ -345,7 +417,7 @@ void FCITransform::integrateParallelSlices(Field3D& f) {
 
   // Integrate f onto yup and ydown fields
   for (const auto& map : field_line_maps) {
-    f.ynext(map.offset) = map.integrate(f);
+    f.ynext(map.offset()) = map.integrate(f);
   }
 }