Field3DParallel + more FCI changes

CMakeLists.txt

@@ -271,11 +271,12 @@ set(BOUT_SOURCES
   ./src/mesh/interpolation/hermite_spline_z.cxx
   ./src/mesh/interpolation/interpolation_z.cxx
   ./src/mesh/interpolation/lagrange_4pt_xz.cxx
-  ./src/mesh/interpolation/monotonic_hermite_spline_xz.cxx
   ./src/mesh/invert3x3.hxx
   ./src/mesh/mesh.cxx
   ./src/mesh/parallel/fci.cxx
   ./src/mesh/parallel/fci.hxx
+  ./src/mesh/parallel/fci_comm.cxx
+  ./src/mesh/parallel/fci_comm.hxx
   ./src/mesh/parallel/identity.cxx
   ./src/mesh/parallel/shiftedmetric.cxx
   ./src/mesh/parallel/shiftedmetricinterp.cxx
@@ -376,6 +377,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, \
@@ -594,6 +596,9 @@ else()
 endif()
 set(BOUT_USE_METRIC_3D ${BOUT_ENABLE_METRIC_3D})
 
+option(BOUT_ENABLE_FCI_AUTOMAGIC "Enable (slow?) automatic features for FCI" ON)
+set(BOUT_USE_FCI_AUTOMAGIC ${BOUT_ENABLE_FCI_AUTOMAGIC})
+
 include(CheckCXXSourceCompiles)
 check_cxx_source_compiles("int main() { const char* name = __PRETTY_FUNCTION__; }"
   HAS_PRETTY_FUNCTION)

bin/bout-build-deps.sh

@@ -10,7 +10,7 @@ NCVER=${NCVER:-4.7.4}
 NCCXXVER=${NCCXXVER:-4.3.1}
 FFTWVER=${FFTWVER:-3.3.9}
 SUNVER=${SUNVER:-5.7.0}
-PETSCVER=${PETSCVER:-3.15.0}
+PETSCVER=${PETSCVER:-3.21.4}
 
 
 HDF5FLAGS=${HDF5FLAGS:-}
@@ -147,7 +147,7 @@ petsc() {
     test -z $PETSC_DIR || error "\$PETSC_DIR is set ($PETSC_DIR) - please unset"
     test -z $PETSC_ARCH || error "\$PETSC_ARCH is set ($PETSC_ARCH) - please unset"
     cd $BUILD
-    wget -c https://ftp.mcs.anl.gov/pub/petsc/release-snapshots/petsc-$PETSCVER.tar.gz || :
+    wget -c https://web.cels.anl.gov/projects/petsc/download/release-snapshots/petsc-$PETSCVER.tar.gz || :
     tar -xf petsc-$PETSCVER.tar.gz
     cd petsc-$PETSCVER
     unset PETSC_DIR

cmake/SetupBOUTThirdParty.cmake

@@ -285,7 +285,7 @@ if (BOUT_USE_SUNDIALS)
     set(EXAMPLES_ENABLE_C OFF CACHE BOOL "" FORCE)
     set(EXAMPLES_INSTALL OFF CACHE BOOL "" FORCE)
     set(ENABLE_MPI ${BOUT_USE_MPI} CACHE BOOL "" FORCE)
-    set(ENABLE_OPENMP OFF CACHE BOOL "" FORCE)
+    set(ENABLE_OPENMP ${BOUT_USE_OPENMP} CACHE BOOL "" FORCE)
     if (BUILD_SHARED_LIBS)
       set(BUILD_STATIC_LIBS OFF CACHE BOOL "" FORCE)
     else()

cmake_build_defines.hxx.in

@@ -35,6 +35,7 @@
 #cmakedefine BOUT_METRIC_TYPE @BOUT_METRIC_TYPE@
 #cmakedefine01 BOUT_USE_METRIC_3D
 #cmakedefine01 BOUT_USE_MSGSTACK
+#cmakedefine01 BOUT_USE_FCI_AUTOMAGIC
 
 // CMake build does not support legacy interface
 #define BOUT_HAS_LEGACY_NETCDF 0

examples/fci-wave/div-integrate/BOUT.inp

@@ -38,7 +38,7 @@ bndry_par_ydown = parallel_neumann_o2
 
 [v]
 
-bndry_par_all = parallel_dirichlet
+bndry_par_all = parallel_dirichlet_o2
 
 [solver]
 output_step = 0.1

examples/fci-wave/div/BOUT.inp

@@ -38,7 +38,7 @@ bndry_par_ydown = parallel_neumann_o2
 
 [v]
 
-bndry_par_all = parallel_dirichlet
+bndry_par_all = parallel_dirichlet_o2
 
 [solver]
 output_step = 0.1

examples/fci-wave/logn/BOUT.inp

@@ -38,7 +38,7 @@ bndry_par_ydown = parallel_neumann_o2
 
 [nv]
 
-bndry_par_all = parallel_dirichlet
+bndry_par_all = parallel_dirichlet_o2
 
 [solver]
 output_step = 0.1

externalpackages/PVODE/include/pvode/band.h

@@ -138,7 +138,7 @@ typedef struct bandmat_type {
  *                                                                *
  ******************************************************************/
 
-#define PVODE_BAND_ELEM(A,i,j) ((A->data)[j][i-j+(A->smu)])
+#define PVODE_BAND_ELEM(A, i, j) (((A)->data)[j][(i) - (j) + ((A)->smu)])
 
 /******************************************************************
  *                                                                *
@@ -153,7 +153,7 @@ typedef struct bandmat_type {
  *                                                                *
  ******************************************************************/
 
-#define PVODE_BAND_COL(A,j) (((A->data)[j])+(A->smu))
+#define PVODE_BAND_COL(A, j) ((((A)->data)[j]) + ((A)->smu))
 
 /******************************************************************
  *                                                                *

externalpackages/PVODE/precon/band.h

@@ -153,7 +153,7 @@ typedef struct bandmat_type {
  *                                                                *
  ******************************************************************/
 
-#define BAND_COL(A, j) (((A->data)[j]) + (A->smu))
+#define PVODE_BAND_COL(A, j) ((((A)->data)[j]) + ((A)->smu))
 
 /******************************************************************
  *                                                                *

include/bout/boundary_iterator.hxx

@@ -0,0 +1,217 @@
+#pragma once
+
+#include "bout/assert.hxx"
+#include "bout/bout_types.hxx"
+#include "bout/mesh.hxx"
+#include "bout/parallel_boundary_region.hxx"
+#include "bout/region.hxx"
+#include "bout/sys/parallel_stencils.hxx"
+#include "bout/sys/range.hxx"
+
+#include <algorithm>
+#include <functional>
+
+class BoundaryRegionIter {
+public:
+  BoundaryRegionIter(int x, int y, int bx, int by, Mesh* mesh)
+      : dir(bx + by), x(x), y(y), bx(bx), by(by), localmesh(mesh) {
+    ASSERT3(bx * by == 0);
+  }
+  bool operator!=(const BoundaryRegionIter& rhs) const { return ind() != rhs.ind(); }
+
+  Ind3D ind() const { return xyz2ind(x, y, z); }
+  BoundaryRegionIter& operator++() {
+    ASSERT3(z < nz());
+    z++;
+    if (z == nz()) {
+      z = 0;
+      _next();
+    }
+    return *this;
+  }
+  virtual void _next() = 0;
+  BoundaryRegionIter& operator*() { return *this; }
+
+  void dirichlet_o2(Field3D& f, BoutReal value) const {
+    ynext(f) = parallel_stencil::dirichlet_o2(1, f[ind()], 0.5, value);
+  }
+
+  BoutReal extrapolate_grad_o2(const Field3D& f) const { return f[ind()] - yprev(f); }
+
+  BoutReal extrapolate_sheath_o2(const Field3D& f) const {
+    return (f[ind()] * 3 - yprev(f)) * 0.5;
+  }
+
+  BoutReal extrapolate_next_o2(const Field3D& f) const {
+    return (2 * f[ind()]) - yprev(f);
+  }
+
+  BoutReal
+  extrapolate_next_o2(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
+    return (2 * f(0, ind())) - f(0, ind().yp(-by).xp(-bx));
+  }
+
+  BoutReal interpolate_sheath_o2(const Field3D& f) const {
+    return (f[ind()] + ynext(f)) * 0.5;
+  }
+
+  BoutReal
+  interpolate_sheath_o2(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
+    return (f(0, ind()) + f(0, ind().yp(-by).xp(-bx))) * 0.5;
+  }
+
+  BoutReal
+  extrapolate_sheath_o2(const std::function<BoutReal(int yoffset, Ind3D ind)>& f) const {
+    return 0.5 * (3 * f(0, ind()) - f(0, ind().yp(-by).xp(-bx)));
+  }
+
+  BoutReal extrapolate_sheath_free(const Field3D& f, SheathLimitMode mode) const {
+    const BoutReal fac =
+        bout::parallel_boundary_region::limitFreeScale(yprev(f), ythis(f), mode);
+    const BoutReal val = ythis(f);
+    const BoutReal next = mode == SheathLimitMode::linear_free ? val + fac : val * fac;
+    return 0.5 * (val + next);
+  }
+
+  void set_free(Field3D& f, SheathLimitMode mode) const {
+    const BoutReal fac =
+        bout::parallel_boundary_region::limitFreeScale(yprev(f), ythis(f), mode);
+    BoutReal val = ythis(f);
+    if (mode == SheathLimitMode::linear_free) {
+      for (int i = 1; i <= localmesh->ystart; ++i) {
+        val += fac;
+        f[ind().yp(by * i).xp(bx * i)] = val;
+      }
+    } else {
+      for (int i = 1; i <= localmesh->ystart; ++i) {
+        val *= fac;
+        f[ind().yp(by * i).xp(bx * i)] = val;
+      }
+    }
+  }
+
+  void limitFree(Field3D& f) const {
+    const BoutReal fac =
+        bout::parallel_boundary_region::limitFreeScale(yprev(f), ythis(f));
+    BoutReal val = ythis(f);
+    for (int i = 1; i <= localmesh->ystart; ++i) {
+      val *= fac;
+      f[ind().yp(by * i).xp(bx * i)] = val;
+    }
+  }
+
+  bool is_lower() const {
+    ASSERT2(bx == 0);
+    return by == -1;
+  }
+
+  void neumann_o1(Field3D& f, BoutReal grad) const {
+    BoutReal val = ythis(f);
+    for (int i = 1; i <= localmesh->ystart; ++i) {
+      val += grad;
+      f[ind().yp(by * i).xp(bx * i)] = val;
+    }
+  }
+
+  void neumann_o2(Field3D& f, BoutReal grad) const {
+    BoutReal val = yprev(f) + grad;
+    for (int i = 1; i <= localmesh->ystart; ++i) {
+      val += grad;
+      f[ind().yp(by * i).xp(bx * i)] = val;
+    }
+  }
+
+  void limit_at_least(Field3D& f, BoutReal value) const {
+    if (ynext(f) < value) {
+      ynext(f) = value;
+    }
+  }
+
+  BoutReal& ynext(Field3D& f) const { return f[ind().yp(by).xp(bx)]; }
+  const BoutReal& ynext(const Field3D& f) const { return f[ind().yp(by).xp(bx)]; }
+  BoutReal& yprev(Field3D& f) const { return f[ind().yp(-by).xp(-bx)]; }
+  const BoutReal& yprev(const Field3D& f) const { return f[ind().yp(-by).xp(-bx)]; }
+  BoutReal& ythis(Field3D& f) const { return f[ind()]; }
+  const BoutReal& ythis(const Field3D& f) const { return f[ind()]; }
+
+  void setYPrevIfValid(Field3D& f, BoutReal val) const { yprev(f) = val; }
+  void setAll(Field3D& f, const BoutReal val) const {
+    for (int i = -localmesh->ystart; i <= localmesh->ystart; ++i) {
+      f[ind().yp(by * i).xp(bx * i)] = val;
+    }
+  }
+
+  static int abs_offset() { return 1; }
+
+#if BOUT_USE_METRIC_3D == 0
+  BoutReal& ynext(Field2D& f) const { return f[ind().yp(by).xp(bx)]; }
+  const BoutReal& ynext(const Field2D& f) const { return f[ind().yp(by).xp(bx)]; }
+  BoutReal& yprev(Field2D& f) const { return f[ind().yp(-by).xp(-bx)]; }
+  const BoutReal& yprev(const Field2D& f) const { return f[ind().yp(-by).xp(-bx)]; }
+#endif
+
+  const int dir;
+  virtual ~BoundaryRegionIter() = default;
+
+protected:
+  int z{0};
+  int x;
+  int y;
+  const int bx;
+  const int by;
+
+private:
+  Mesh* localmesh;
+  int nx() const { return localmesh->LocalNx; }
+  int ny() const { return localmesh->LocalNy; }
+  int nz() const { return localmesh->LocalNz; }
+
+  Ind3D xyz2ind(int x, int y, int z) const {
+    return Ind3D{(x * ny() + y) * nz() + z, ny(), nz()};
+  }
+};
+
+class BoundaryRegionIterY : public BoundaryRegionIter {
+public:
+  BoundaryRegionIterY(RangeIterator r, int y, int dir, bool is_end, Mesh* mesh)
+      : BoundaryRegionIter(r.ind, y, 0, dir, mesh), r(r), is_end(is_end) {}
+
+  bool operator!=(const BoundaryRegionIterY& rhs) {
+    ASSERT2(y == rhs.y);
+    if (is_end) {
+      if (rhs.is_end) {
+        return false;
+      }
+      return !rhs.r.isDone();
+    }
+    if (rhs.is_end) {
+      return !r.isDone();
+    }
+    return x != rhs.x;
+  }
+
+  virtual void _next() override {
+    ++r;
+    x = r.ind;
+  }
+
+private:
+  RangeIterator r;
+  bool is_end;
+};
+
+class NewBoundaryRegionY {
+public:
+  NewBoundaryRegionY(Mesh* mesh, bool lower, RangeIterator r)
+      : mesh(mesh), lower(lower), r(std::move(r)) {}
+  BoundaryRegionIterY begin(bool begin = true) {
+    return BoundaryRegionIterY(r, lower ? mesh->ystart : mesh->yend, lower ? -1 : +1,
+                               !begin, mesh);
+  }
+  BoundaryRegionIterY end() { return begin(false); }
+
+private:
+  Mesh* mesh;
+  bool lower;
+  RangeIterator r;
+};

include/bout/boundary_region.hxx

@@ -4,6 +4,9 @@ class BoundaryRegion;
 #ifndef BOUT_BNDRY_REGION_H
 #define BOUT_BNDRY_REGION_H
 
+#include "bout/mesh.hxx"
+#include "bout/region.hxx"
+#include "bout/sys/parallel_stencils.hxx"
 #include <string>
 #include <utility>