diff --git a/manual/sphinx/user_docs/variable_init.rst b/manual/sphinx/user_docs/variable_init.rst index 4ac13e0ead..1baaba5e84 100644 --- a/manual/sphinx/user_docs/variable_init.rst +++ b/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`` diff --git a/src/field/field_factory.cxx b/src/field/field_factory.cxx index f65f2e7f55..b93a39eeda 100644 --- a/src/field/field_factory.cxx +++ b/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(nullptr, nullptr, nullptr)); + + // Periodic in the Y direction? + addGenerator("is_periodic_y", std::make_shared()); } Field2D FieldFactory::create2D(const std::string& value, const Options* opt, diff --git a/src/field/fieldgenerators.cxx b/src/field/fieldgenerators.cxx index 1f21f5c262..7e728664f8 100644 --- a/src/field/fieldgenerators.cxx +++ b/src/field/fieldgenerators.cxx @@ -1,6 +1,8 @@ #include "fieldgenerators.hxx" +#include + #include #include diff --git a/src/field/fieldgenerators.hxx b/src/field/fieldgenerators.hxx index 2485b4b82d..050e335448 100644 --- a/src/field/fieldgenerators.hxx +++ b/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 UNUSED(args)) override { + return std::make_shared(); + } + 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 diff --git a/tests/unit/fake_mesh.hxx b/tests/unit/fake_mesh.hxx index 5d0d12de07..2feb43826d 100644 --- a/tests/unit/fake_mesh.hxx +++ b/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 hasBranchCutLower(int UNUSED(jx)) const override { return std::make_pair(false, 0.); diff --git a/tests/unit/field/test_field_factory.cxx b/tests/unit/field/test_field_factory.cxx index 964788b25a..26dc1e2990 100644 --- a/tests/unit/field/test_field_factory.cxx +++ b/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( + [](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( + &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(localmesh)); + + localmesh.ix_separatrix = 0; // All points outside core + + auto output = this->create("is_periodic_y", nullptr, &localmesh); + + auto expected = makeField( + [](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( + &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(localmesh)); + + auto output = this->create("is_periodic_y", nullptr, &localmesh); + + auto expected = makeField( + [&](typename TypeParam::ind_type& index) -> BoutReal { + return index.x() < localmesh.ix_separatrix; + }, + &localmesh); + + EXPECT_TRUE(IsFieldEqual(output, expected)); +}