diff --git a/include/bout/field.hxx b/include/bout/field.hxx index a8e1952546..ad3ca5cf6b 100644 --- a/include/bout/field.hxx +++ b/include/bout/field.hxx @@ -630,9 +630,9 @@ std::optional getPerpYIndex(const BinaryExpr& expr) { #ifdef FIELD_FUNC #error This macro has already been defined #else -#define FIELD_FUNC(name, func) \ +#define FIELD_FUNC(opname, func) \ namespace bout::op { \ - struct name { \ + struct opname { \ template \ BOUT_HOST_DEVICE BoutReal operator()(int idx, const LView& L, const RView&) const { \ return func(L(idx)); \ @@ -640,7 +640,7 @@ std::optional getPerpYIndex(const BinaryExpr& expr) { }; \ }; \ template > \ - inline auto name(const T& f, const std::string& rgn = "RGN_ALL") { \ + inline auto opname(const T& f, const std::string& rgn = "RGN_ALL") { \ if constexpr (std::is_same_v) { \ /* Check if the input is allocated */ \ checkData(f); \ @@ -651,28 +651,28 @@ std::optional getPerpYIndex(const BinaryExpr& expr) { result.yup(i) = func(f.yup(i)); \ result.ydown(i) = func(f.ydown(i)); \ } \ - result.name = std::string(#name "(") + f.name + std::string(")"); \ + result.name = fmt::format(#opname "({})", f.name); \ checkData(result); \ return result; \ } else { \ - return BinaryExpr{static_cast(f), \ - static_cast(f), \ - bout::op::name{}, \ - f.getMesh(), \ - f.getLocation(), \ - f.getDirections(), \ - std::nullopt, \ - f.getRegion(rgn), \ - bout::detail::getPerpYIndex(f)}; \ + return BinaryExpr{static_cast(f), \ + static_cast(f), \ + bout::op::opname{}, \ + f.getMesh(), \ + f.getLocation(), \ + f.getDirections(), \ + std::nullopt, \ + f.getRegion(rgn), \ + bout::detail::getPerpYIndex(f)}; \ } \ } \ template \ - inline auto name(const BinaryExpr& f) { \ + inline auto opname(const BinaryExpr& f) { \ return BinaryExpr, BinaryExpr, \ - bout::op::name>{ \ + bout::op::opname>{ \ static_cast::View>(f), \ static_cast::View>(f), \ - bout::op::name{}, \ + bout::op::opname{}, \ f.getMesh(), \ f.getLocation(), \ f.getDirections(), \ @@ -681,8 +681,8 @@ std::optional getPerpYIndex(const BinaryExpr& expr) { bout::detail::getPerpYIndex(f)}; \ } \ template \ - inline auto name(const BinaryExpr& f, const std::string& rgn) { \ - return name(ResT{f}, rgn); \ + inline auto opname(const BinaryExpr& f, const std::string& rgn) { \ + return opname(ResT{f}, rgn); \ } #endif @@ -709,7 +709,9 @@ inline auto SQ(const T& f, const std::string& rgn = "RGN_ALL") { result.yup(i) = SQ(f.yup(i), rgn); result.ydown(i) = SQ(f.ydown(i), rgn); } +#if TRACK result.name = std::string("SQ(") + f.name + std::string(")"); +#endif checkData(result); return result; } else { diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index 905e736999..77413f292e 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -684,6 +684,88 @@ protected: } }; +/// Field3DParallel is intended to behave like Field3D, but preserve parallel +/// Fields. +/// Operations on Field3D, like multiplication, exp and floor only work on the +/// "main" field, Field3DParallel will retain the parallel slices. +class Field3DParallel : public Field3D { +public: + template + explicit Field3DParallel(Types... args) : Field3D(std::move(args)...) { + ensureFieldAligned(); + } + Field3DParallel(const Field3D& f) : Field3D(f) { ensureFieldAligned(); } + Field3DParallel(const Field3D& f, bool isRef) : Field3D(f), isRef(isRef) { + ensureFieldAligned(); + } + Field3DParallel(const Field2D& f) : Field3D(f) { ensureFieldAligned(); } + // Explicitly needed, as DirectionTypes is sometimes constructed from a + // brace enclosed list + explicit Field3DParallel(Mesh* localmesh = nullptr, CELL_LOC location_in = CELL_CENTRE, + DirectionTypes directions_in = {YDirectionType::Standard, + ZDirectionType::Standard}, + std::optional regionID = {}) + : Field3D(localmesh, location_in, directions_in, regionID) { + if (isFci()) { + splitParallelSlices(); + } + ensureFieldAligned(); + } + explicit Field3DParallel(Array data, Mesh* localmesh, + CELL_LOC location = CELL_CENTRE, + DirectionTypes directions_in = {YDirectionType::Standard, + ZDirectionType::Standard}) + : Field3D(std::move(data), localmesh, location, directions_in) { + ensureFieldAligned(); + } + explicit Field3DParallel(BoutReal, Mesh* mesh = nullptr); + Field3D& asField3D() { return *this; } + const Field3D& asField3D() const { return *this; } + + Field3DParallel& operator*=(const Field3D&); + Field3DParallel& operator/=(const Field3D&); + Field3DParallel& operator+=(const Field3D&); + Field3DParallel& operator-=(const Field3D&); + Field3DParallel& operator*=(const Field3DParallel&); + Field3DParallel& operator/=(const Field3DParallel&); + Field3DParallel& operator+=(const Field3DParallel&); + Field3DParallel& operator-=(const Field3DParallel&); + Field3DParallel& operator*=(BoutReal); + Field3DParallel& operator/=(BoutReal); + Field3DParallel& operator+=(BoutReal); + Field3DParallel& operator-=(BoutReal); + Field3DParallel& operator=(const Field3D& rhs) { + Field3D::operator=(rhs); + ensureFieldAligned(); + return *this; + } + Field3DParallel& operator=(Field3D&& rhs) { + Field3D::operator=(std::move(rhs)); + ensureFieldAligned(); + return *this; + } + Field3DParallel& operator=(BoutReal); + Field3DParallel& allocate(); + +private: + void ensureFieldAligned(); + bool isRef{false}; +}; + +// We need checkData for templates +#if CHECK > 0 +/// Throw an exception if \p f is not allocated or if any +/// elements are non-finite (for CHECK > 2). +/// Loops over all points including the boundaries by +/// default (can be changed using the \p rgn argument +void checkData(const Field3D& f, const std::string& region = "RGN_NOBNDRY"); +#else +/// Ignored with disabled CHECK; Throw an exception if \p f is not +/// allocated or if any elements are non-finite (for CHECK > 2) +inline void checkData([[maybe_unused]] const Field3D& f, + [[maybe_unused]] const std::string& region = "RGN_NOBNDRY") {}; +#endif + // Non-member overloaded operators template @@ -698,9 +780,10 @@ FieldPerp operator/(const Field3D& lhs, const FieldPerp& rhs); #define FIELD3D_FIELD3D_FIELD3D_OP(OP_SYM, OP_TYPE) \ template && is_expr_field3d_v>> \ - BinaryExpr operator OP_SYM(const L& lhs, \ - const R& rhs) { \ + auto operator OP_SYM(const L& lhs, const R& rhs) { \ ASSERT1_EXPR_COMPATIBLE(lhs, rhs); \ + static_assert(!std::is_same_v); \ + static_assert(!std::is_same_v); \ auto regionID = \ lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID()); \ return BinaryExpr{ \ @@ -911,19 +994,6 @@ Field3D pow(const Field3D& lhs, const Field2D& rhs, const std::string& rgn = "RG FieldPerp pow(const Field3D& lhs, const FieldPerp& rhs, const std::string& rgn = "RGN_ALL"); -#if CHECK > 0 -/// Throw an exception if \p f is not allocated or if any -/// elements are non-finite (for CHECK > 2). -/// Loops over all points including the boundaries by -/// default (can be changed using the \p rgn argument -void checkData(const Field3D& f, const std::string& region = "RGN_NOBNDRY"); -#else -/// Ignored with disabled CHECK; Throw an exception if \p f is not -/// allocated or if any elements are non-finite (for CHECK > 2) -inline void checkData([[maybe_unused]] const Field3D& f, - [[maybe_unused]] const std::string& region = "RGN_NOBNDRY") {}; -#endif - /// Fourier filtering, removes all except one mode /// /// @param[in] var Variable to apply filter to @@ -1014,74 +1084,6 @@ inline Field3D copy(const Field3D& f) { return result; } -/// Field3DParallel is intended to behave like Field3D, but preserve parallel -/// Fields. -/// Operations on Field3D, like multiplication, exp and floor only work on the -/// "main" field, Field3DParallel will retain the parallel slices. -class Field3DParallel : public Field3D { -public: - template - explicit Field3DParallel(Types... args) : Field3D(std::move(args)...) { - ensureFieldAligned(); - } - Field3DParallel(const Field3D& f) : Field3D(f) { ensureFieldAligned(); } - Field3DParallel(const Field3D& f, bool isRef) : Field3D(f), isRef(isRef) { - ensureFieldAligned(); - } - Field3DParallel(const Field2D& f) : Field3D(f) { ensureFieldAligned(); } - // Explicitly needed, as DirectionTypes is sometimes constructed from a - // brace enclosed list - explicit Field3DParallel(Mesh* localmesh = nullptr, CELL_LOC location_in = CELL_CENTRE, - DirectionTypes directions_in = {YDirectionType::Standard, - ZDirectionType::Standard}, - std::optional regionID = {}) - : Field3D(localmesh, location_in, directions_in, regionID) { - if (isFci()) { - splitParallelSlices(); - } - ensureFieldAligned(); - } - explicit Field3DParallel(Array data, Mesh* localmesh, - CELL_LOC location = CELL_CENTRE, - DirectionTypes directions_in = {YDirectionType::Standard, - ZDirectionType::Standard}) - : Field3D(std::move(data), localmesh, location, directions_in) { - ensureFieldAligned(); - } - explicit Field3DParallel(BoutReal, Mesh* mesh = nullptr); - Field3D& asField3D() { return *this; } - const Field3D& asField3D() const { return *this; } - - Field3DParallel& operator*=(const Field3D&); - Field3DParallel& operator/=(const Field3D&); - Field3DParallel& operator+=(const Field3D&); - Field3DParallel& operator-=(const Field3D&); - Field3DParallel& operator*=(const Field3DParallel&); - Field3DParallel& operator/=(const Field3DParallel&); - Field3DParallel& operator+=(const Field3DParallel&); - Field3DParallel& operator-=(const Field3DParallel&); - Field3DParallel& operator*=(BoutReal); - Field3DParallel& operator/=(BoutReal); - Field3DParallel& operator+=(BoutReal); - Field3DParallel& operator-=(BoutReal); - Field3DParallel& operator=(const Field3D& rhs) { - Field3D::operator=(rhs); - ensureFieldAligned(); - return *this; - } - Field3DParallel& operator=(Field3D&& rhs) { - Field3D::operator=(std::move(rhs)); - ensureFieldAligned(); - return *this; - } - Field3DParallel& operator=(BoutReal); - Field3DParallel& allocate(); - -private: - void ensureFieldAligned(); - bool isRef{false}; -}; - Field3DParallel Field3D::asField3DParallel() { if (isAllocated()) { allocate(); diff --git a/include/bout/fieldops.hxx b/include/bout/fieldops.hxx index f36061ceaa..b26c0c07f7 100644 --- a/include/bout/fieldops.hxx +++ b/include/bout/fieldops.hxx @@ -56,10 +56,10 @@ struct is_expr_fieldperp : std::false_type {}; template <> struct is_expr_fieldperp : std::true_type {}; +// Helper variable template template inline constexpr bool is_expr_fieldperp_v = is_expr_fieldperp>::value; -// Helper variable template template inline constexpr bool is_expr_field3d_v = is_expr_field3d>::value; diff --git a/tests/unit/fake_mesh.hxx b/tests/unit/fake_mesh.hxx index f957652a62..2ed085de1c 100644 --- a/tests/unit/fake_mesh.hxx +++ b/tests/unit/fake_mesh.hxx @@ -1,5 +1,6 @@ #pragma once +#include "bout/paralleltransform.hxx" #include #include #include @@ -366,3 +367,61 @@ public: private: Options values; ///< Store values to be returned by get() }; + +// A mock ParallelTransform to test transform_from_field_aligned +// property of FieldFactory. For now, the transform just returns the +// negative of the input. Ideally, this will get moved to GoogleMock +// when we start using it. +// +// Can turn off the ability to do the transform. Should still be valid +class MockParallelTransform : public ParallelTransform { +public: + MockParallelTransform(Mesh& mesh, bool allow_transform_) + : ParallelTransform(mesh), allow_transform(allow_transform_) {} + MockParallelTransform(const MockParallelTransform&) = delete; + MockParallelTransform& operator=(const MockParallelTransform&) = delete; + MockParallelTransform(MockParallelTransform&& other) = delete; + MockParallelTransform& operator=(MockParallelTransform&& other) = delete; + ~MockParallelTransform() override = default; + + void calcParallelSlices(Field3D& /*f*/) override {} + + bool canToFromFieldAligned() const override { return allow_transform; } + + bool requiresTwistShift(bool /*twist_shift_enabled*/, + YDirectionType /*ytype*/) override { + return false; + } + + void checkInputGrid() override {} + + Field3D fromFieldAligned(const Field3D& f, const std::string& /*region*/) override { + if (f.getDirectionY() != YDirectionType::Aligned) { + throw BoutException("Unaligned field passed to fromFieldAligned"); + } + return -f; + } + + FieldPerp fromFieldAligned(const FieldPerp& f, const std::string& /*region*/) override { + if (f.getDirectionY() != YDirectionType::Aligned) { + throw BoutException("Unaligned field passed to fromFieldAligned"); + } + return -f; + } + + Field3D toFieldAligned(const Field3D& f, const std::string& /*region*/) override { + if (f.getDirectionY() != YDirectionType::Standard) { + throw BoutException("Aligned field passed to toFieldAligned"); + } + return -f; + } + FieldPerp toFieldAligned(const FieldPerp& f, const std::string& /*region*/) override { + if (f.getDirectionY() != YDirectionType::Standard) { + throw BoutException("Aligned field passed to toFieldAligned"); + } + return -f; + } + +private: + bool allow_transform; +}; diff --git a/tests/unit/fake_mesh_fixture.hxx b/tests/unit/fake_mesh_fixture.hxx index 2758dbe416..3c261a6121 100644 --- a/tests/unit/fake_mesh_fixture.hxx +++ b/tests/unit/fake_mesh_fixture.hxx @@ -31,7 +31,7 @@ /// Use this template class directly to use different sized grid: /// /// using MyTest = FakeMeshFixture_tmpl<7, 9, 11>; -template +template class FakeMeshFixture_tmpl : public ::testing::Test { public: FakeMeshFixture_tmpl() @@ -113,6 +113,11 @@ public: mesh_staggered_m.setCoordinates(test_coords_staggered, CELL_XLOW); mesh_staggered_m.setCoordinates(test_coords_staggered, CELL_YLOW); mesh_staggered_m.setCoordinates(test_coords_staggered, CELL_ZLOW); + + if constexpr (FCI) { + mesh_m.getCoordinates()->setParallelTransform( + bout::utils::make_unique(mesh_m, false)); + } } FakeMeshFixture_tmpl(const FakeMeshFixture_tmpl&) = delete; diff --git a/tests/unit/field/test_field3d.cxx b/tests/unit/field/test_field3d.cxx index 905b182018..7f73f6d85b 100644 --- a/tests/unit/field/test_field3d.cxx +++ b/tests/unit/field/test_field3d.cxx @@ -26,6 +26,7 @@ using namespace bout::globals; // Reuse the "standard" fixture for FakeMesh using Field3DTest = FakeMeshFixture; +using Field3DTestFCI = FakeMeshFixture_tmpl<3, 5, 7, true>; TEST_F(Field3DTest, Is3D) { Field3D field; @@ -1979,6 +1980,147 @@ TEST_F(Field3DTest, SQField3DParallelPreservesParallelSlices) { EXPECT_TRUE(IsFieldEqual(squared.ydown(), 16.0)); } +TEST_F(Field3DTestFCI, SqrtField3DParallelPreservesParallelSlices) { + Field3DParallel field; + + field = 4.0; + field.splitParallelSlices(); + field.yup() = 9.0; + field.ydown() = 16.0; + field.resetRegionParallel(); + + { + const Field3D res = sqrt(field); + + EXPECT_TRUE(res.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(res, 2.0)); + EXPECT_TRUE(IsFieldEqual(res.yup(), 3.0)); + EXPECT_TRUE(IsFieldEqual(res.ydown(), 4.0)); + } + { + const Field3D res = sqrt(field.asField3D()); + + EXPECT_FALSE(res.hasParallelSlices()); + } +} + +TEST_F(Field3DTestFCI, MulField3DParallelPreservesParallelSlices) { + Field3D field; + EXPECT_TRUE(field.isFci()); + + field = 2.0; + field.splitParallelSlices(); + field.yup() = 3.0; + field.ydown() = 4.0; + field.resetRegionParallel(); + + Field3D rhs; + EXPECT_TRUE(rhs.isFci()); + rhs = 3.0; + rhs.splitParallelSlices(); + rhs.yup() = 4.0; + rhs.ydown() = 5.0; + rhs.resetRegionParallel(); + + const Field3D prod = field * rhs; + + EXPECT_FALSE(prod.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(prod, 6.0)); + + const Field3DParallel prodpar = field.asField3DParallel() * rhs; + EXPECT_TRUE((std::is_same_v, Field3DParallel>)); + EXPECT_TRUE(prodpar.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(prodpar, 6.0)); + EXPECT_TRUE(IsFieldEqual(prodpar.yup(), 12.0, "RGN_YPAR_+1")); + EXPECT_TRUE(IsFieldEqual(prodpar.ydown(), 20.0, "RGN_YPAR_-1")); +} + +TEST_F(Field3DTestFCI, DivField3DParallelPreservesParallelSlices) { + Field3D field; + + field = 2.0; + field.splitParallelSlices(); + field.yup() = 3.0; + field.ydown() = 4.0; + field.resetRegionParallel(); + + Field3DParallel rhs{1.0}; + EXPECT_TRUE(IsFieldEqual(rhs.ydown(), 1, "RGN_YPAR_-1")); + rhs *= 3.0; + EXPECT_TRUE(IsFieldEqual(rhs.ydown(), 3, "RGN_YPAR_-1")); + rhs.yup() *= 4; + EXPECT_TRUE(IsFieldEqual(rhs.ydown(), 3, "RGN_YPAR_-1")); + rhs.ydown() *= 20. / 3; + EXPECT_TRUE(IsFieldEqual(rhs.ydown(), 20, "RGN_YPAR_-1")); + + const Field3D prod = field / rhs.asField3D(); + + EXPECT_FALSE(prod.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(prod, 2. / 3)); + + const Field3DParallel prodpar = field / rhs; + EXPECT_TRUE(prodpar.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(prodpar, 2. / 3)); + EXPECT_TRUE(IsFieldEqual(prodpar.yup(), .25, "RGN_YPAR_+1")); + EXPECT_TRUE(IsFieldEqual(prodpar.ydown(), 4. / 20, "RGN_YPAR_-1")); +} + +TEST_F(Field3DTestFCI, AddField3DParallelPreservesParallelSlices) { + Field3D field; + field = 2.0; + field.splitParallelSlices(); + field.yup() = 3.0; + field.ydown() = 4.0; + field.resetRegionParallel(); + + Field3D rhs; + rhs = 3.0; + rhs.splitParallelSlices(); + rhs.yup() = 4.0; + rhs.ydown() = 5.0; + rhs.resetRegionParallel(); + + const Field3D res = field + rhs; + + EXPECT_FALSE(res.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(res, 5.0)); + + const Field3DParallel respar = field.asField3DParallel() + rhs; + EXPECT_TRUE((std::is_same_v, Field3DParallel>)); + EXPECT_TRUE(respar.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(respar, 5.0)); + EXPECT_TRUE(IsFieldEqual(respar.yup(), 7.0, "RGN_YPAR_+1")); + EXPECT_TRUE(IsFieldEqual(respar.ydown(), 9.0, "RGN_YPAR_-1")); +} + +TEST_F(Field3DTestFCI, SubField3DParallelPreservesParallelSlices) { + Field3D field; + field = 2.0; + field.splitParallelSlices(); + field.yup() = 3.0; + field.ydown() = 4.0; + field.resetRegionParallel(); + + Field3D rhs; + rhs = 3.0; + rhs.splitParallelSlices(); + rhs.yup() = 5.0; + rhs.ydown() = 7.0; + rhs.resetRegionParallel(); + + const Field3D res = field - rhs; + + EXPECT_FALSE(res.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(res, -1.0)); + + const Field3DParallel respar = field.asField3DParallel() - rhs; + EXPECT_TRUE((std::is_same_v, Field3DParallel>)); + EXPECT_TRUE(respar.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(respar, -1.0)); + EXPECT_TRUE(IsFieldEqual(respar.yup(), -2.0, "RGN_YPAR_+1")); + EXPECT_TRUE(IsFieldEqual(respar.ydown(), -3.0, "RGN_YPAR_-1")); +} + TEST_F(Field3DTest, Abs) { Field3D field; diff --git a/tests/unit/field/test_field_factory.cxx b/tests/unit/field/test_field_factory.cxx index 9db9fcef10..8cb08ed230 100644 --- a/tests/unit/field/test_field_factory.cxx +++ b/tests/unit/field/test_field_factory.cxx @@ -832,57 +832,6 @@ TEST_F(FieldFactoryTest, FuzzyFind) { EXPECT_EQ(CAPS_matches.size(), 1); } -// A mock ParallelTransform to test transform_from_field_aligned -// property of FieldFactory. For now, the transform just returns the -// negative of the input. Ideally, this will get moved to GoogleMock -// when we start using it. -// -// Can turn off the ability to do the transform. Should still be valid -class MockParallelTransform : public ParallelTransform { -public: - MockParallelTransform(Mesh& mesh, bool allow_transform_) - : ParallelTransform(mesh), allow_transform(allow_transform_) {} - ~MockParallelTransform() = default; - - void calcParallelSlices(Field3D&) override {} - - bool canToFromFieldAligned() const override { return allow_transform; } - - bool requiresTwistShift(bool, YDirectionType) override { return false; } - - void checkInputGrid() override {} - - Field3D fromFieldAligned(const Field3D& f, const std::string&) override { - if (f.getDirectionY() != YDirectionType::Aligned) { - throw BoutException("Unaligned field passed to fromFieldAligned"); - } - return -f; - } - - FieldPerp fromFieldAligned(const FieldPerp& f, const std::string&) override { - if (f.getDirectionY() != YDirectionType::Aligned) { - throw BoutException("Unaligned field passed to fromFieldAligned"); - } - return -f; - } - - Field3D toFieldAligned(const Field3D& f, const std::string&) override { - if (f.getDirectionY() != YDirectionType::Standard) { - throw BoutException("Aligned field passed to toFieldAligned"); - } - return -f; - } - FieldPerp toFieldAligned(const FieldPerp& f, const std::string&) override { - if (f.getDirectionY() != YDirectionType::Standard) { - throw BoutException("Aligned field passed to toFieldAligned"); - } - return -f; - } - -private: - const bool allow_transform; -}; - class FieldFactoryCreateAndTransformTest : public FakeMeshFixture { public: WithQuietOutput quiet_info{output_info};