Options is_periodic_y function

manual/sphinx/user_docs/variable_init.rst

@@ -81,17 +81,19 @@ following values are also already defined:
 .. _tab-initexprvals:
 .. table:: Initialisation expression values
 
-   +--------+------------------------------------------------------------------------------------+
-   | Name   | Description                                                                        |
-   +========+====================================================================================+
-   | x      | :math:`x` position between :math:`0` and :math:`1`                                 |
-   +--------+------------------------------------------------------------------------------------+
-   | y      | :math:`y` angle-like position, definition depends on topology of grid              |
-   +--------+------------------------------------------------------------------------------------+
-   | z      | :math:`z` position between :math:`0` and :math:`2\pi` (excluding the last point)   |
-   +--------+------------------------------------------------------------------------------------+
-   | pi  π  | :math:`3.1415\ldots`                                                               |
-   +--------+------------------------------------------------------------------------------------+
+   +----------------+------------------------------------------------------------------------------------+
+   | Name           | Description                                                                        |
+   +================+====================================================================================+
+   | x              | :math:`x` position between :math:`0` and :math:`1`                                 |
+   +----------------+------------------------------------------------------------------------------------+
+   | y              | :math:`y` angle-like position, definition depends on topology of grid              |
+   +----------------+------------------------------------------------------------------------------------+
+   | z              | :math:`z` position between :math:`0` and :math:`2\pi` (excluding the last point)   |
+   +----------------+------------------------------------------------------------------------------------+
+   | pi  π          | :math:`3.1415\ldots`                                                               |
+   +----------------+------------------------------------------------------------------------------------+
+   | is_periodic_y  | :math:`1` in core region where Y is periodic. :math:`0` otherwise                  |
+   +----------------+------------------------------------------------------------------------------------+
 
 
 By default, :math:`x` is defined as ``(i+0.5) / (nx - 2*MXG)``, where ``MXG``

src/field/field_factory.cxx

@@ -1,7 +1,7 @@
 /**************************************************************************
- * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu
+ * Copyright 2010-2025 BOUT++ contributors
  *
- * Contact: Ben Dudson, bd512@york.ac.uk
+ * Contact: Ben Dudson, dudson2@llnl.gov
  *
  * This file is part of BOUT++.
  *
@@ -176,6 +176,9 @@ FieldFactory::FieldFactory(Mesh* localmesh, Options* opt)
 
   // Where switch function
   addGenerator("where", std::make_shared<FieldWhere>(nullptr, nullptr, nullptr));
+
+  // Periodic in the Y direction?
+  addGenerator("is_periodic_y", std::make_shared<FieldPeriodicY>());
 }
 
 Field2D FieldFactory::create2D(const std::string& value, const Options* opt,

src/field/fieldgenerators.cxx

@@ -1,6 +1,8 @@
 
 #include "fieldgenerators.hxx"
 
+#include <memory>
+
 #include <bout/constants.hxx>
 #include <bout/utils.hxx>
 

src/field/fieldgenerators.hxx

@@ -352,4 +352,29 @@ private:
   FieldGeneratorPtr test, gt0, lt0;
 };
 
+/// Function that evaluates to 1 when Y is periodic (i.e. in the core), 0 otherwise
+/// Note: Assumes symmetricGlobalX
+class FieldPeriodicY : public FieldGenerator {
+public:
+  FieldPeriodicY() = default;
+  FieldGeneratorPtr clone(const std::list<FieldGeneratorPtr> UNUSED(args)) override {
+    return std::make_shared<FieldPeriodicY>();
+  }
+  BoutReal generate(const bout::generator::Context& ctx) override {
+    const Mesh* mesh = ctx.getMesh();
+    const BoutReal local_inner_boundary =
+        0.5 * (mesh->GlobalX(mesh->xstart - 1) + mesh->GlobalX(mesh->xstart));
+    const BoutReal local_outer_boundary =
+        0.5 * (mesh->GlobalX(mesh->xend + 1) + mesh->GlobalX(mesh->xend));
+    const int local_index = mesh->xstart
+                            + int(((ctx.x() - local_inner_boundary)
+                                   / (local_outer_boundary - local_inner_boundary))
+                                  * (mesh->xend - mesh->xstart + 1));
+    if (mesh->periodicY(local_index)) {
+      return 1.0;
+    }
+    return 0.0;
+  }
+};
+
 #endif // BOUT_FIELDGENERATORS_H

tests/unit/fake_mesh.hxx

@@ -131,8 +131,14 @@ public:
   }
   MPI_Comm getXcomm(int UNUSED(jy)) const override { return BoutComm::get(); }
   MPI_Comm getYcomm(int UNUSED(jx)) const override { return BoutComm::get(); }
-  bool periodicY(int UNUSED(jx)) const override { return true; }
-  bool periodicY(int UNUSED(jx), BoutReal& UNUSED(ts)) const override { return true; }
+
+  // Periodic Y
+  int ix_separatrix{1000000}; // separatrix index
+
+  bool periodicY(int jx) const override { return jx < ix_separatrix; }
+  bool periodicY(int jx, BoutReal& UNUSED(ts)) const override {
+    return jx < ix_separatrix;
+  }
   int numberOfYBoundaries() const override { return 1; }
   std::pair<bool, BoutReal> hasBranchCutLower(int UNUSED(jx)) const override {
     return std::make_pair(false, 0.);

tests/unit/field/test_field_factory.cxx

@@ -969,3 +969,63 @@ TEST_F(FieldFactoryCreateAndTransformTest, Create3DCantTransform) {
 
   EXPECT_TRUE(IsFieldEqual(output, expected));
 }
+
+TYPED_TEST(FieldFactoryCreationTest, CreatePeriodicY) {
+  auto output = this->create("is_periodic_y");
+
+  auto expected = makeField<TypeParam>(
+      [](typename TypeParam::ind_type& index) -> BoutReal {
+        return mesh->periodicY(index.x());
+      },
+      mesh);
+
+  EXPECT_TRUE(IsFieldEqual(output, expected));
+}
+
+TYPED_TEST(FieldFactoryCreationTest, CreatePeriodicYoutsideCore) {
+  FakeMesh localmesh{5, 1, 1};
+  localmesh.createDefaultRegions();
+  localmesh.setCoordinates(std::make_shared<Coordinates>(
+      &localmesh, Field2D{1.0}, Field2D{1.0}, BoutReal{1.0}, Field2D{1.0}, Field2D{0.0},
+      Field2D{1.0}, Field2D{1.0}, Field2D{1.0}, Field2D{0.0}, Field2D{0.0}, Field2D{0.0},
+      Field2D{1.0}, Field2D{1.0}, Field2D{1.0}, Field2D{0.0}, Field2D{0.0}, Field2D{0.0},
+      Field2D{0.0}, Field2D{0.0}));
+  // No call to Coordinates::geometry() needed here
+
+  localmesh.getCoordinates()->setParallelTransform(
+      bout::utils::make_unique<ParallelTransformIdentity>(localmesh));
+
+  localmesh.ix_separatrix = 0; // All points outside core
+
+  auto output = this->create("is_periodic_y", nullptr, &localmesh);
+
+  auto expected = makeField<TypeParam>(
+      [](typename TypeParam::ind_type& index) -> BoutReal { return 0.0; }, &localmesh);
+
+  EXPECT_TRUE(IsFieldEqual(output, expected));
+}
+
+TYPED_TEST(FieldFactoryCreationTest, CreatePeriodicYacrossSeparatrix) {
+  FakeMesh localmesh{5, 1, 1};
+  localmesh.createDefaultRegions();
+  localmesh.ix_separatrix = 2; // All points in core
+  localmesh.setCoordinates(std::make_shared<Coordinates>(
+      &localmesh, Field2D{1.0}, Field2D{1.0}, BoutReal{1.0}, Field2D{1.0}, Field2D{0.0},
+      Field2D{1.0}, Field2D{1.0}, Field2D{1.0}, Field2D{0.0}, Field2D{0.0}, Field2D{0.0},
+      Field2D{1.0}, Field2D{1.0}, Field2D{1.0}, Field2D{0.0}, Field2D{0.0}, Field2D{0.0},
+      Field2D{0.0}, Field2D{0.0}));
+  // No call to Coordinates::geometry() needed here
+
+  localmesh.getCoordinates()->setParallelTransform(
+      bout::utils::make_unique<ParallelTransformIdentity>(localmesh));
+
+  auto output = this->create("is_periodic_y", nullptr, &localmesh);
+
+  auto expected = makeField<TypeParam>(
+      [&](typename TypeParam::ind_type& index) -> BoutReal {
+        return index.x() < localmesh.ix_separatrix;
+      },
+      &localmesh);
+
+  EXPECT_TRUE(IsFieldEqual(output, expected));
+}