From 23ddf98e76471e9b28be2be94ae4367c57d039c6 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 9 Jul 2026 13:48:40 -0700 Subject: [PATCH 01/14] Delete unneededed expr.hxx This was an experiment in expression templates, that has now been replaced by the BinaryExpr machinery. --- CMakeLists.txt | 1 - examples/performance/arithmetic/.gitignore | 1 - .../performance/arithmetic/arithmetic.cxx | 106 --------- examples/performance/arithmetic/data/BOUT.inp | 5 - examples/performance/arithmetic/run.sh | 5 - .../performance/arithmetic_3d2d/.gitignore | 1 - .../arithmetic_3d2d/arithmetic_3d2d.cxx | 116 ---------- .../performance/arithmetic_3d2d/data/BOUT.inp | 5 - examples/performance/arithmetic_3d2d/run.sh | 5 - include/bout/expr.hxx | 208 ------------------ 10 files changed, 453 deletions(-) delete mode 100644 examples/performance/arithmetic/.gitignore delete mode 100644 examples/performance/arithmetic/arithmetic.cxx delete mode 100644 examples/performance/arithmetic/data/BOUT.inp delete mode 100755 examples/performance/arithmetic/run.sh delete mode 100644 examples/performance/arithmetic_3d2d/.gitignore delete mode 100644 examples/performance/arithmetic_3d2d/arithmetic_3d2d.cxx delete mode 100644 examples/performance/arithmetic_3d2d/data/BOUT.inp delete mode 100644 examples/performance/arithmetic_3d2d/run.sh delete mode 100644 include/bout/expr.hxx diff --git a/CMakeLists.txt b/CMakeLists.txt index e6e424b55c..23919726cc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -128,7 +128,6 @@ set(BOUT_SOURCES ./include/bout/deriv_store.hxx ./include/bout/derivs.hxx ./include/bout/difops.hxx - ./include/bout/expr.hxx ./include/bout/fft.hxx ./include/bout/field.hxx ./include/bout/field2d.hxx diff --git a/examples/performance/arithmetic/.gitignore b/examples/performance/arithmetic/.gitignore deleted file mode 100644 index 077be4cbd0..0000000000 --- a/examples/performance/arithmetic/.gitignore +++ /dev/null @@ -1 +0,0 @@ -arithmetic \ No newline at end of file diff --git a/examples/performance/arithmetic/arithmetic.cxx b/examples/performance/arithmetic/arithmetic.cxx deleted file mode 100644 index fc2357978a..0000000000 --- a/examples/performance/arithmetic/arithmetic.cxx +++ /dev/null @@ -1,106 +0,0 @@ -/* - * Timing of arithmetic operations - * - */ - -#include - -#include - -#include - -using SteadyClock = std::chrono::time_point; -using Duration = std::chrono::duration; -using namespace std::chrono; - -#define TIMEIT(elapsed, ...) \ - { \ - SteadyClock start = steady_clock::now(); \ - { __VA_ARGS__; } \ - Duration diff = steady_clock::now() - start; \ - diff *= 1000 * 1000; \ - elapsed.min = diff > elapsed.min ? elapsed.min : diff; \ - elapsed.max = diff < elapsed.max ? elapsed.max : diff; \ - elapsed.count++; \ - elapsed.avg = elapsed.avg * (1 - 1. / elapsed.count) + diff / elapsed.count; \ - } - -struct Durations { - Duration max; - Duration min; - Duration avg; - int count; -}; - -class Arithmetic : public PhysicsModel { -protected: - int init(bool) { - - Field3D a = 1.0; - Field3D b = 2.0; - Field3D c = 3.0; - a.setRegion("RGN_ALL"); - b.setRegion("RGN_NOBNDRY"); - - Field3D result1, result2, result3, result4; - - // Using Field methods (classic operator overloading) - - result1 = 2. * a + b * c; -#define dur_init {Duration::min(), Duration::max(), Duration::zero(), 0} - Durations elapsed1 = dur_init, elapsed2 = dur_init, elapsed3 = dur_init, - elapsed4 = dur_init; - - for (int ik = 0; ik < 1e3; ++ik) { - TIMEIT(elapsed1, result1 = 2. * a + b * c;); - - // Using C loops - result2.allocate(); - BoutReal* rd = &result2(0, 0, 0); - BoutReal* ad = &a(0, 0, 0); - BoutReal* bd = &b(0, 0, 0); - BoutReal* cd = &c(0, 0, 0); - TIMEIT( - elapsed2, - for (int i = 0, iend = (mesh->LocalNx * mesh->LocalNy * mesh->LocalNz) - 1; - i != iend; i++) { - *rd = 2. * (*ad) + (*bd) * (*cd); - rd++; - ad++; - bd++; - cd++; - }); - - // Template expressions - TIMEIT(elapsed3, result3 = eval3D(add(mul(2, a), mul(b, c)));); - - // Range iterator - result4.allocate(); - TIMEIT(elapsed4, for (auto i : result4) result4[i] = 2. * a[i] + b[i] * c[i];); - } - - output.enable(); - output << "TIMING | minimum | mean | maximum\n" - << "----------- | ---------- | ---------- | ----------\n"; - //#define PRINT(str,elapsed) output << str << elapsed.min.count()<< - //elapsed.avg.count()<< elapsed.max.count() << endl; -#define PRINT(str, elapsed) \ - output.write("{:s} | {:7.3f} us | {:7.3f} us | {:7.3f} us\n", str, \ - elapsed.min.count(), elapsed.avg.count(), elapsed.max.count()) - PRINT("Fields: ", elapsed1); - PRINT("C loop: ", elapsed2); - PRINT("Templates: ", elapsed3); - PRINT("Range For: ", elapsed4); - output.disable(); - SOLVE_FOR(n); - return 0; - } - - int rhs(BoutReal) { - ddt(n) = 0; - return 0; - } - Field3D n; -}; - -BOUTMAIN(Arithmetic); diff --git a/examples/performance/arithmetic/data/BOUT.inp b/examples/performance/arithmetic/data/BOUT.inp deleted file mode 100644 index 0deb623c4b..0000000000 --- a/examples/performance/arithmetic/data/BOUT.inp +++ /dev/null @@ -1,5 +0,0 @@ -MZ = 1024 - -[mesh] -nx = 50 -ny = 2 diff --git a/examples/performance/arithmetic/run.sh b/examples/performance/arithmetic/run.sh deleted file mode 100755 index 3a1cc844a6..0000000000 --- a/examples/performance/arithmetic/run.sh +++ /dev/null @@ -1,5 +0,0 @@ -#!/usr/bin/env bash - -make || exit - -./arithmetic -q -q -q diff --git a/examples/performance/arithmetic_3d2d/.gitignore b/examples/performance/arithmetic_3d2d/.gitignore deleted file mode 100644 index 14968af063..0000000000 --- a/examples/performance/arithmetic_3d2d/.gitignore +++ /dev/null @@ -1 +0,0 @@ -arithmetic_3d2d \ No newline at end of file diff --git a/examples/performance/arithmetic_3d2d/arithmetic_3d2d.cxx b/examples/performance/arithmetic_3d2d/arithmetic_3d2d.cxx deleted file mode 100644 index 83167a5b42..0000000000 --- a/examples/performance/arithmetic_3d2d/arithmetic_3d2d.cxx +++ /dev/null @@ -1,116 +0,0 @@ -/* - * Timing of arithmetic operations (Field3D/Field2D mixed) - * - */ - -#include - -#include - -#include -#include -#include - -using SteadyClock = std::chrono::time_point; -using Duration = std::chrono::duration; -using namespace std::chrono; - -#define TIMEIT(NAME, ...) \ - { \ - SteadyClock start = steady_clock::now(); \ - __VA_ARGS__ \ - Duration diff = steady_clock::now() - start; \ - auto elapsed = elapsedMap[NAME]; \ - elapsed.min = std::min(diff, elapsed.min); \ - elapsed.max = std::max(diff, elapsed.max); \ - elapsed.count++; \ - elapsed.avg = elapsed.avg * (1 - 1. / elapsed.count) + diff / elapsed.count; \ - elapsedMap[NAME] = elapsed; \ - } - -struct Durations { - Duration max; - Duration min; - Duration avg; - int count; - Durations() - : max(Duration::min()), min(Duration::max()), avg(Duration::zero()), count(0){}; -}; - -class Arithmetic : public PhysicsModel { -protected: - std::map elapsedMap; - - int init(bool) { - Field3D a = 1.0; - Field3D b = 2.0; - Field2D c = 3.0; - - Field3D result1, result2, result3, result4; - - // Using Field methods (classic operator overloading) - result1 = 2. * a + b * c; - - for (int ik = 0; ik < 1e2; ++ik) { - result1.allocate(); - TIMEIT("Fields", result1 = 2. * a + b * c;); - - // Using C loops - result2.allocate(); - BoutReal* rd = &result2(0, 0, 0); - BoutReal* ad = &a(0, 0, 0); - BoutReal* bd = &b(0, 0, 0); - BoutReal* cd = &c(0, 0, 0); - TIMEIT( - "C loop", - for (int i = 0, iend = (mesh->LocalNx * mesh->LocalNy) - 1; i != iend; i++) { - for (int j = 0, jend = mesh->LocalNz - 1; j != jend; j++) { - *rd = 2. * (*ad) + (*bd) * (*cd); - rd++; - ad++; - bd++; - } - cd++; - }); - - // Template expressions - result3.allocate(); - TIMEIT("Templates", result3 = eval3D(add(mul(2, a), mul(b, c)));); - - // Range iterator - result4.allocate(); - TIMEIT("Range For", for (auto i : result4) result4[i] = 2. * a[i] + b[i] * c[i];); - } - - output.enable(); - constexpr int width = 15; - output << std::setw(width) << "TIMING"; - output << std::setw(width) << "min"; - output << std::setw(width) << "avg"; - output << std::setw(width) << "max"; - output << "\n======"; - for (int i = 0; i < 4 * width; ++i) { - output << "="; - }; - output << "\n"; - - for (const auto& approach : elapsedMap) { - output << std::setw(width) << approach.first; - output << std::setw(width) << approach.second.min.count(); - output << std::setw(width) << approach.second.avg.count(); - output << std::setw(width) << approach.second.max.count(); - output << "\n"; - } - output.disable(); - SOLVE_FOR(n); - return 0; - } - - int rhs(BoutReal) { - ddt(n) = 0; - return 0; - } - Field3D n; -}; - -BOUTMAIN(Arithmetic); diff --git a/examples/performance/arithmetic_3d2d/data/BOUT.inp b/examples/performance/arithmetic_3d2d/data/BOUT.inp deleted file mode 100644 index 0deb623c4b..0000000000 --- a/examples/performance/arithmetic_3d2d/data/BOUT.inp +++ /dev/null @@ -1,5 +0,0 @@ -MZ = 1024 - -[mesh] -nx = 50 -ny = 2 diff --git a/examples/performance/arithmetic_3d2d/run.sh b/examples/performance/arithmetic_3d2d/run.sh deleted file mode 100644 index ee36808c21..0000000000 --- a/examples/performance/arithmetic_3d2d/run.sh +++ /dev/null @@ -1,5 +0,0 @@ -#!/usr/bin/env bash - -make || exit - -./arithmetic_3d2d -q -q -q diff --git a/include/bout/expr.hxx b/include/bout/expr.hxx deleted file mode 100644 index 267af202ed..0000000000 --- a/include/bout/expr.hxx +++ /dev/null @@ -1,208 +0,0 @@ -/************************************************************************** - * - * Operators, and support for template expressions - * - * Originally based on article by Klaus Kreft & Angelika Langer - * http://www.angelikalanger.com/Articles/Cuj/ExpressionTemplates/ExpressionTemplates.htm - * - * Parts adapted from Blitz++ library - * - **************************************************************************/ - -#ifndef BOUT_EXPR_H -#define BOUT_EXPR_H - -#warning expr.hxx is deprecated. Do not use! - -#include -#include -#include - -/// Literal class to capture BoutReal values in expressions -class Literal { -public: - /// Type of this expression - using type = Literal; - - Literal(BoutReal v) : val(v) {} - ~Literal() {} - BoutReal operator()(int x, int y, int z) const { return val; } - -private: - const BoutReal val; -}; - -class Field3DExpr { -public: - using type = Field3D; - - Field3DExpr(const Field3D& f) : data(&f(0, 0, 0)) {} - const BoutReal& operator()(int x, int y, int z) const { - return data[(x * bout::globals::mesh->LocalNy + y) * bout::globals::mesh->LocalNz - + z]; - } - -private: - const BoutReal* data; -}; - -class Field2DExpr { -public: - using type = Field2D; - - Field2DExpr(const Field2D& f) : data(&f(0, 0)) {} - const BoutReal& operator()(int x, int y, int z) const { - return data[x * bout::globals::mesh->LocalNy + y]; - } - -private: - const BoutReal* data; -}; - -/// Expression traits, to convert doubles etc. to Literal - -template -struct exprTraits { - using expr_type = ExprT; -}; - -template <> -struct exprTraits { - using expr_type = Literal; -}; - -template <> -struct exprTraits { - using expr_type = Literal; -}; - -template <> -struct exprTraits { - using expr_type = Literal; -}; - -/////////////////////////////////////////////// -// asExpr: convert objects to expressions - -template -struct asExpr { - using type = T; - static const T& getExpr(const T& x) { return x; } -}; - -template <> -struct asExpr { - using type = Literal; - static const Literal getExpr(const int& x) { return Literal(x); } -}; - -template <> -struct asExpr { - using type = Literal; - static const Literal getExpr(const double& x) { return Literal(x); } -}; - -template <> -struct asExpr { - using type = Literal; - static const Literal getExpr(const float& x) { return Literal(x); } -}; - -template <> -struct asExpr { - using type = Field3DExpr; - static const Field3DExpr getExpr(const Field3D& x) { return Field3DExpr(x); } -}; - -///////////////////////////////////////////////////////////// -// Type promotion. Work out the type of a calculation, -// based on the type of the arguments - -template // If in doubt, convert to Field3D -struct PromoteType { - using type = Field3D; -}; - -///////////////////////////////////////////////////////////// -// Binary expressions - -template -class BinaryExpr { -public: - BinaryExpr(const ExprT1& e1, const ExprT2& e2) : _expr1(e1), _expr2(e2) {} - - // Work out the type of the inputs - using ltype = typename exprTraits::expr_type; - using rtype = typename exprTraits::expr_type; - - /// Type of the resulting expression - using type = typename PromoteType::type; - - BoutReal operator()(int x, int y, int z) const { - return BinOp::apply((_expr1)(x, y, z), (_expr2)(x, y, z)); - } - -private: - ltype const _expr1; - rtype const _expr2; -}; - -template -struct BinaryResult { - using arg1 = typename asExpr::type; - using arg2 = typename asExpr::type; - using type = BinaryExpr; -}; - -/// Binary operator classes - -#define DEFINE_BINARY_OP(name, op) \ - struct name { \ - template \ - static inline T apply(T a, T b) { \ - return a op b; \ - } \ - }; - -DEFINE_BINARY_OP(Add, +) -DEFINE_BINARY_OP(Subtract, -) -DEFINE_BINARY_OP(Multiply, *) -DEFINE_BINARY_OP(Divide, /) - -struct Power { - template - static inline T apply(T a, T b) { - return pow(a, b); - } -}; - -/// Define functions add, mul which use operator structs -#define DEFINE_OVERLOAD_FUNC(name, func) \ - template \ - typename BinaryResult::type func(const ExprT1& e1, \ - const ExprT2& e2) { \ - using type = typename BinaryResult::type; \ - return type(asExpr::getExpr(e1), asExpr::getExpr(e2)); \ - } - -/// Addition of two Expressions -DEFINE_OVERLOAD_FUNC(Add, add); -/// Multiplication of two Expressions -DEFINE_OVERLOAD_FUNC(Multiply, mul); - -/// A function to evaluate expressions -template -const Field3D eval3D(Expr e) { - Field3D result; - result.allocate(); - for (int i = 0; i < bout::globals::mesh->LocalNx; i++) { - for (int j = 0; j < bout::globals::mesh->LocalNy; j++) { - for (int k = 0; k < bout::globals::mesh->LocalNz; k++) { - result(i, j, k) = e(i, j, k); - } - } - } - return result; -} - -#endif // BOUT_EXPR_H From 49774de6ffb6701bd26e8ef67434e7cb25435b1a Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 9 Jul 2026 20:47:10 -0700 Subject: [PATCH 02/14] Field3DParallel BinaryExpr templates Views handle yup/ydown slices, evaluating lazy expressions in parallel slices if assigned to a Field3DParallel. --- include/bout/bout_types.hxx | 4 + include/bout/field2d.hxx | 4 + include/bout/field3d.hxx | 147 ++++++++++++++++++++++------ include/bout/fieldops.hxx | 76 ++++++++++++++- include/bout/fieldperp.hxx | 4 + tests/unit/field/test_field3d.cxx | 155 ++++++++++++++++++++++++++++++ 6 files changed, 362 insertions(+), 28 deletions(-) diff --git a/include/bout/bout_types.hxx b/include/bout/bout_types.hxx index 7747b937b3..dee2ea22e7 100644 --- a/include/bout/bout_types.hxx +++ b/include/bout/bout_types.hxx @@ -149,6 +149,10 @@ struct Constant { T v; View(T v) : v(v) {} BOUT_HOST_DEVICE T operator()(int) const { return v; } + BOUT_HOST_DEVICE bool hasParallelSlices() const { return false; } + BOUT_HOST_DEVICE int numberParallelSlices() const { return 0; } + BOUT_HOST_DEVICE View yup(int = 0) const { return *this; } + BOUT_HOST_DEVICE View ydown(int = 0) const { return *this; } }; operator View() const { return {val}; } }; diff --git a/include/bout/field2d.hxx b/include/bout/field2d.hxx index 540680bb73..ce3009f4ad 100644 --- a/include/bout/field2d.hxx +++ b/include/bout/field2d.hxx @@ -345,6 +345,10 @@ public: this->div = div; return *this; } + BOUT_HOST_DEVICE BOUT_FORCEINLINE bool hasParallelSlices() const { return false; } + BOUT_HOST_DEVICE BOUT_FORCEINLINE int numberParallelSlices() const { return 0; } + BOUT_HOST_DEVICE BOUT_FORCEINLINE View yup(int = 0) const { return *this; } + BOUT_HOST_DEVICE BOUT_FORCEINLINE View ydown(int = 0) const { return *this; } }; operator View() { return View{&data[0]}; } operator View() const { return View{const_cast(&data[0])}; } diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index 905e736999..02e582c20a 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -465,6 +465,10 @@ public: struct View { BoutReal* data; + const Field3D* yup_fields{nullptr}; + const Field3D* ydown_fields{nullptr}; + int num_parallel_slices{0}; + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal operator()(int idx) const { return data[idx]; } @@ -478,9 +482,34 @@ public: "Field3D::View does not support setScale()"); return *this; } + BOUT_HOST_DEVICE BOUT_FORCEINLINE bool hasParallelSlices() const { + return num_parallel_slices > 0; + } + BOUT_HOST_DEVICE BOUT_FORCEINLINE int numberParallelSlices() const { + return num_parallel_slices; + } + BOUT_FORCEINLINE View yup(int slice = 0) const { + ASSERT2(slice < num_parallel_slices); + ASSERT2(yup_fields[slice].isAllocated()); + return static_cast(yup_fields[slice]); + } + BOUT_FORCEINLINE View ydown(int slice = 0) const { + ASSERT2(slice < num_parallel_slices); + ASSERT2(ydown_fields[slice].isAllocated()); + return static_cast(ydown_fields[slice]); + } }; - operator View() { return View{&data[0]}; } - operator View() const { return View{const_cast(&data[0])}; } + operator View() { + return View{&data[0], hasParallelSlices() ? yup_fields.data() : nullptr, + hasParallelSlices() ? ydown_fields.data() : nullptr, + static_cast(numberParallelSlices())}; + } + operator View() const { + return View{const_cast(&data[0]), + hasParallelSlices() ? yup_fields.data() : nullptr, + hasParallelSlices() ? ydown_fields.data() : nullptr, + static_cast(numberParallelSlices())}; + } //operator View() const { return View{&data[0]}; } ///////////////////////////////////////////////////////// @@ -856,31 +885,6 @@ if_else(bool condition, const L& lhs, const R& rhs) { rhs.getMesh()->getRegion("RGN_ALL")}; } -Field3DParallel operator+(const Field3D& lhs, const Field3DParallel& rhs); -Field3DParallel operator-(const Field3D& lhs, const Field3DParallel& rhs); -Field3DParallel operator*(const Field3D& lhs, const Field3DParallel& rhs); -Field3DParallel operator/(const Field3D& lhs, const Field3DParallel& rhs); - -Field3DParallel operator+(const Field3DParallel& lhs, const Field3D& rhs); -Field3DParallel operator-(const Field3DParallel& lhs, const Field3D& rhs); -Field3DParallel operator*(const Field3DParallel& lhs, const Field3D& rhs); -Field3DParallel operator/(const Field3DParallel& lhs, const Field3D& rhs); - -Field3DParallel operator+(const Field3DParallel& lhs, const Field3DParallel& rhs); -Field3DParallel operator-(const Field3DParallel& lhs, const Field3DParallel& rhs); -Field3DParallel operator*(const Field3DParallel& lhs, const Field3DParallel& rhs); -Field3DParallel operator/(const Field3DParallel& lhs, const Field3DParallel& rhs); - -Field3DParallel operator+(BoutReal lhs, const Field3DParallel& rhs); -Field3DParallel operator-(BoutReal lhs, const Field3DParallel& rhs); -Field3DParallel operator*(BoutReal lhs, const Field3DParallel& rhs); -Field3DParallel operator/(BoutReal lhs, const Field3DParallel& rhs); - -Field3DParallel operator+(const Field3DParallel& lhs, BoutReal rhs); -Field3DParallel operator-(const Field3DParallel& lhs, BoutReal rhs); -Field3DParallel operator*(const Field3DParallel& lhs, BoutReal rhs); -Field3DParallel operator/(const Field3DParallel& lhs, BoutReal rhs); - /*! * Unary minus. Returns the negative of given field, * iterates over whole domain including guard/boundary cells. @@ -1024,6 +1028,13 @@ public: explicit Field3DParallel(Types... args) : Field3D(std::move(args)...) { ensureFieldAligned(); } + template || is_expr_field3d_v>> + Field3DParallel(const BinaryExpr& expr) + : Field3DParallel(expr.getMesh(), expr.getLocation(), expr.getDirections(), + expr.getRegionID()) { + *this = expr; + } Field3DParallel(const Field3D& f) : Field3D(f) { ensureFieldAligned(); } Field3DParallel(const Field3D& f, bool isRef) : Field3D(f), isRef(isRef) { ensureFieldAligned(); @@ -1052,6 +1063,53 @@ public: Field3D& asField3D() { return *this; } const Field3D& asField3D() const { return *this; } + struct View { + Field3D::View base; + const Field3D* yup_fields{nullptr}; + const Field3D* ydown_fields{nullptr}; + int num_parallel_slices{0}; + + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal operator()(int idx) const { + return base(idx); + } + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal& operator[](int idx) const { + return base[idx]; + } + + template + View& setScale(Mul /*unused*/, Div /*unused*/) { + static_assert(!std::is_same_v, + "Field3DParallel::View does not support setScale()"); + return *this; + } + + BOUT_FORCEINLINE bool hasParallelSlices() const { return num_parallel_slices > 0; } + BOUT_FORCEINLINE int numberParallelSlices() const { return num_parallel_slices; } + BOUT_FORCEINLINE View yup(int slice = 0) const { + ASSERT2(slice < num_parallel_slices); + ASSERT2(yup_fields[slice].isAllocated()); + return View{static_cast(yup_fields[slice]), nullptr, nullptr, 0}; + } + BOUT_FORCEINLINE View ydown(int slice = 0) const { + ASSERT2(slice < num_parallel_slices); + ASSERT2(ydown_fields[slice].isAllocated()); + return View{static_cast(ydown_fields[slice]), nullptr, nullptr, 0}; + } + }; + + operator View() { + return View{static_cast(*this), + hasParallelSlices() ? yup_fields.data() : nullptr, + hasParallelSlices() ? ydown_fields.data() : nullptr, + static_cast(numberParallelSlices())}; + } + operator View() const { + return View{static_cast(*this), + hasParallelSlices() ? yup_fields.data() : nullptr, + hasParallelSlices() ? ydown_fields.data() : nullptr, + static_cast(numberParallelSlices())}; + } + Field3DParallel& operator*=(const Field3D&); Field3DParallel& operator/=(const Field3D&); Field3DParallel& operator+=(const Field3D&); @@ -1074,6 +1132,41 @@ public: ensureFieldAligned(); return *this; } + template + std::enable_if_t || is_expr_field3d_v, Field3DParallel&> + operator=(const BinaryExpr& expr) { + if (getMesh() != expr.getMesh()) { + clearParallelSlices(); + fieldmesh = expr.getMesh(); + data = Array{}; + } + if (isFci()) { + if (!hasParallelSlices()) { + splitParallelSlices(); + } + } else if (hasParallelSlices()) { + clearParallelSlices(); + } + + setRegion(expr.getRegionID()); + setLocation(expr.getLocation()); + setDirections(expr.getDirections()); + allocate(); + expr.evaluate(static_cast(*this).data); + + if (isFci()) { + ASSERT2(expr.hasParallelSlices()); + ASSERT2(expr.numberParallelSlices() == static_cast(numberParallelSlices())); + for (int i = 0; i < expr.numberParallelSlices(); ++i) { + yup(i).allocate(); + ydown(i).allocate(); + expr.yup(i).evaluate(static_cast(yup(i)).data); + expr.ydown(i).evaluate(static_cast(ydown(i)).data); + } + } + + return *this; + } Field3DParallel& operator=(BoutReal); Field3DParallel& allocate(); diff --git a/include/bout/fieldops.hxx b/include/bout/fieldops.hxx index f36061ceaa..a09734ad21 100644 --- a/include/bout/fieldops.hxx +++ b/include/bout/fieldops.hxx @@ -381,6 +381,43 @@ struct BinaryExpr { } BOUT_HOST_DEVICE BOUT_FORCEINLINE int regionIdx(int idx) const { return indices[idx]; } + bool hasParallelSlices() const { + if constexpr (is_expr_constant_v && is_expr_constant_v) { + return false; + } else if constexpr (is_expr_constant_v) { + return rhs.hasParallelSlices(); + } else if constexpr (is_expr_constant_v) { + return lhs.hasParallelSlices(); + } else { + return lhs.hasParallelSlices() && rhs.hasParallelSlices(); + } + } + int numberParallelSlices() const { + if (!hasParallelSlices()) { + return 0; + } + if constexpr (is_expr_constant_v && is_expr_constant_v) { + return 0; + } else if constexpr (is_expr_constant_v) { + return rhs.numberParallelSlices(); + } else if constexpr (is_expr_constant_v) { + return lhs.numberParallelSlices(); + } else { + ASSERT2(lhs.numberParallelSlices() == rhs.numberParallelSlices()); + return lhs.numberParallelSlices(); + } + } + auto yup(int slice = 0) const { + return BinaryExpr{lhs.yup(slice), rhs.yup(slice), f, + mesh, location, directions, + regionID, indices, yindex}; + } + auto ydown(int slice = 0) const { + return BinaryExpr{ + lhs.ydown(slice), rhs.ydown(slice), f, mesh, location, + directions, regionID, indices, yindex}; + } + //operator ResT() { return ResT{*this}; } struct View { typename L::View lhs; @@ -396,13 +433,50 @@ struct BinaryExpr { this->div = div; return *this; } + BOUT_HOST_DEVICE BOUT_FORCEINLINE bool hasParallelSlices() const { + if constexpr (is_expr_constant_v && is_expr_constant_v) { + return false; + } else if constexpr (is_expr_constant_v) { + return rhs.hasParallelSlices(); + } else if constexpr (is_expr_constant_v) { + return lhs.hasParallelSlices(); + } else { + return lhs.hasParallelSlices() && rhs.hasParallelSlices(); + } + } + BOUT_HOST_DEVICE BOUT_FORCEINLINE int numberParallelSlices() const { + if (!hasParallelSlices()) { + return 0; + } + if constexpr (is_expr_constant_v && is_expr_constant_v) { + return 0; + } else if constexpr (is_expr_constant_v) { + return rhs.numberParallelSlices(); + } else if constexpr (is_expr_constant_v) { + return lhs.numberParallelSlices(); + } else { + ASSERT2(lhs.numberParallelSlices() == rhs.numberParallelSlices()); + return lhs.numberParallelSlices(); + } + } + BOUT_HOST_DEVICE BOUT_FORCEINLINE auto yup(int slice = 0) const { + auto result = *this; + result.lhs = lhs.yup(slice); + result.rhs = rhs.yup(slice); + return result; + } + BOUT_HOST_DEVICE BOUT_FORCEINLINE auto ydown(int slice = 0) const { + auto result = *this; + result.lhs = lhs.ydown(slice); + result.rhs = rhs.ydown(slice); + return result; + } BOUT_HOST_DEVICE BOUT_FORCEINLINE int size() const { return num_indices; } BOUT_HOST_DEVICE BOUT_FORCEINLINE int regionIdx(int idx) const { return indices[idx]; } BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal operator()(int idx) const { return f((idx * mul) / div, lhs, rhs); // single‐pass fusion - //return f(lhs((idx * mul) / div), rhs((idx * mul) / div)); // single‐pass fusion } }; diff --git a/include/bout/fieldperp.hxx b/include/bout/fieldperp.hxx index 36a116e1b5..83bb3b79b0 100644 --- a/include/bout/fieldperp.hxx +++ b/include/bout/fieldperp.hxx @@ -352,6 +352,10 @@ public: this->div = div; return *this; } + BOUT_HOST_DEVICE BOUT_FORCEINLINE bool hasParallelSlices() const { return false; } + BOUT_HOST_DEVICE BOUT_FORCEINLINE int numberParallelSlices() const { return 0; } + BOUT_HOST_DEVICE BOUT_FORCEINLINE View yup(int = 0) const { return *this; } + BOUT_HOST_DEVICE BOUT_FORCEINLINE View ydown(int = 0) const { return *this; } }; operator View() { return View{&data[0]}; } operator View() const { return View{const_cast(&data[0])}; } diff --git a/tests/unit/field/test_field3d.cxx b/tests/unit/field/test_field3d.cxx index 905b182018..f0b672f324 100644 --- a/tests/unit/field/test_field3d.cxx +++ b/tests/unit/field/test_field3d.cxx @@ -12,6 +12,7 @@ #include "bout/field3d.hxx" #include "bout/mesh.hxx" #include "bout/output.hxx" +#include "bout/paralleltransform.hxx" #include "bout/unused.hxx" #include "bout/utils.hxx" @@ -27,6 +28,34 @@ using namespace bout::globals; // Reuse the "standard" fixture for FakeMesh using Field3DTest = FakeMeshFixture; +class MockFciParallelTransform : public ParallelTransform { +public: + explicit MockFciParallelTransform(Mesh& mesh_in) : ParallelTransform(mesh_in) {} + + void calcParallelSlices(Field3D& f) override { + if (!f.hasParallelSlices()) { + f.splitParallelSlices(); + } + for (size_t i = 0; i < f.numberParallelSlices(); ++i) { + f.yup(i) = f; + f.ydown(i) = f; + } + } + + Field3D toFieldAligned(const Field3D& f, const std::string&) override { return f; } + FieldPerp toFieldAligned(const FieldPerp& f, const std::string&) override { return f; } + Field3D fromFieldAligned(const Field3D& f, const std::string&) override { return f; } + FieldPerp fromFieldAligned(const FieldPerp& f, const std::string&) override { + return f; + } + + bool canToFromFieldAligned() const override { return false; } + bool requiresTwistShift(bool, YDirectionType) override { return false; } + +protected: + void checkInputGrid() override {} +}; + TEST_F(Field3DTest, Is3D) { Field3D field; @@ -1979,6 +2008,132 @@ TEST_F(Field3DTest, SQField3DParallelPreservesParallelSlices) { EXPECT_TRUE(IsFieldEqual(squared.ydown(), 16.0)); } +TEST_F(Field3DTest, Field3DParallelArithmeticReturnsLazyExpr) { + Field3DParallel parallel; + Field3D field; + + parallel = 2.0; + parallel.splitParallelSlices(); + parallel.yup() = 3.0; + parallel.ydown() = 4.0; + field = 5.0; + + const auto expr = parallel + field; + + EXPECT_TRUE( + (std::is_same_v, + BinaryExpr>)); + + Field3DParallel result{expr}; + + EXPECT_FALSE(result.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(result, 7.0)); +} + +TEST_F(Field3DTest, Field3DParallelAssignmentFromLazyExprDiscardsSlicesWhenNotFci) { + Field3DParallel parallel; + Field3DParallel result; + + parallel = 2.0; + parallel.splitParallelSlices(); + parallel.yup() = 3.0; + parallel.ydown() = 4.0; + + const auto expr = 10.0 - parallel; + + EXPECT_TRUE((std::is_same_v< + std::decay_t, + BinaryExpr, Field3DParallel, bout::op::Sub>>)); + + result = expr; + + EXPECT_FALSE(result.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(result, 8.0)); +} + +#if CHECK >= 2 +TEST_F(Field3DTest, Field3DParallelAssignmentFromLazyExprRequiresField3DSlicesInFci) { + mesh->getCoordinates()->setParallelTransform( + bout::utils::make_unique(*mesh)); + + Field3DParallel parallel; + Field3D field; + + parallel = 2.0; + parallel.yup() = 3.0; + parallel.ydown() = 4.0; + field = 5.0; + + EXPECT_THROW(Field3DParallel result{parallel + field}, BoutException); +} +#endif + +TEST_F(Field3DTest, Field3DParallelAssignmentFromLazyExprUsesField3DSlicesInFci) { + mesh->getCoordinates()->setParallelTransform( + bout::utils::make_unique(*mesh)); + + Field3DParallel parallel; + Field3D field; + + parallel = 2.0; + parallel.yup() = 3.0; + parallel.ydown() = 4.0; + + field = 5.0; + field.splitParallelSlices(); + field.yup() = 7.0; + field.ydown() = 11.0; + + Field3DParallel result{parallel + field}; + + EXPECT_TRUE(result.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(result, 7.0)); + EXPECT_TRUE(IsFieldEqual(result.yup(), 10.0)); + EXPECT_TRUE(IsFieldEqual(result.ydown(), 15.0)); +} + +TEST_F(Field3DTest, Field3DParallelAssignmentFromField3DExprUsesSlicesInFci) { + mesh->getCoordinates()->setParallelTransform( + bout::utils::make_unique(*mesh)); + + Field3D lhs; + Field3D rhs; + + lhs = 2.0; + lhs.splitParallelSlices(); + lhs.yup() = 3.0; + lhs.ydown() = 4.0; + + rhs = 5.0; + rhs.splitParallelSlices(); + rhs.yup() = 7.0; + rhs.ydown() = 11.0; + + Field3DParallel result{lhs * rhs}; + + EXPECT_TRUE(result.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(result, 10.0)); + EXPECT_TRUE(IsFieldEqual(result.yup(), 21.0)); + EXPECT_TRUE(IsFieldEqual(result.ydown(), 44.0)); +} + +TEST_F(Field3DTest, Field3DParallelAssignmentFromScalarExprUsesSlicesInFci) { + mesh->getCoordinates()->setParallelTransform( + bout::utils::make_unique(*mesh)); + + Field3DParallel parallel; + parallel = 2.0; + parallel.yup() = 3.0; + parallel.ydown() = 4.0; + + Field3DParallel result{parallel + 1.0}; + + EXPECT_TRUE(result.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(result, 3.0)); + EXPECT_TRUE(IsFieldEqual(result.yup(), 4.0)); + EXPECT_TRUE(IsFieldEqual(result.ydown(), 5.0)); +} + TEST_F(Field3DTest, Abs) { Field3D field; From 470abe40140ca1e1a65634676fabfc6e8fdd27f6 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 8 Jul 2026 14:21:53 +0200 Subject: [PATCH 03/14] Add unit test to ensure Field3DParallel works for products with FCI --- tests/unit/field/test_field3d.cxx | 936 +++++++++++++++--------------- 1 file changed, 484 insertions(+), 452 deletions(-) diff --git a/tests/unit/field/test_field3d.cxx b/tests/unit/field/test_field3d.cxx index f0b672f324..89286a8f80 100644 --- a/tests/unit/field/test_field3d.cxx +++ b/tests/unit/field/test_field3d.cxx @@ -27,6 +27,7 @@ using namespace bout::globals; // Reuse the "standard" fixture for FakeMesh using Field3DTest = FakeMeshFixture; +using Field3DTestFCI = FakeMeshFixture_tmpl<3, 5, 7, true>; class MockFciParallelTransform : public ParallelTransform { public: @@ -2132,618 +2133,649 @@ TEST_F(Field3DTest, Field3DParallelAssignmentFromScalarExprUsesSlicesInFci) { EXPECT_TRUE(IsFieldEqual(result, 3.0)); EXPECT_TRUE(IsFieldEqual(result.yup(), 4.0)); EXPECT_TRUE(IsFieldEqual(result.ydown(), 5.0)); -} -TEST_F(Field3DTest, Abs) { - Field3D field; + 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")); + } - field = -31.0; - EXPECT_TRUE(IsFieldEqual(abs(field), 31.0)); -} + TEST_F(Field3DTest, Abs) { + Field3D field; -TEST_F(Field3DTest, AbsExpressionUsesAbsOp) { - Field3D field; + field = -31.0; + EXPECT_TRUE(IsFieldEqual(abs(field), 31.0)); + } - field = -2.0; - const auto expr = field + 1.0; + TEST_F(Field3DTest, AbsExpressionUsesAbsOp) { + Field3D field; - EXPECT_TRUE((std::is_same_v, - BinaryExpr, - std::decay_t, bout::op::abs>>)); - EXPECT_TRUE(IsFieldEqual(abs(expr), 1.0)); - EXPECT_TRUE(IsFieldEqual(abs(expr, "RGN_ALL"), 1.0)); -} + field = -2.0; + const auto expr = field + 1.0; -TEST_F(Field3DTest, RegionLimitedExpressionConstructsField3D) { - Field3D field; + EXPECT_TRUE( + (std::is_same_v, + BinaryExpr, + std::decay_t, bout::op::abs>>)); + EXPECT_TRUE(IsFieldEqual(abs(expr), 1.0)); + EXPECT_TRUE(IsFieldEqual(abs(expr, "RGN_ALL"), 1.0)); + } - field = -31.0; + TEST_F(Field3DTest, RegionLimitedExpressionConstructsField3D) { + Field3D field; - Field3D result = abs(field, "RGN_NOBNDRY"); + field = -31.0; - EXPECT_TRUE(IsFieldEqual(result, 31.0, "RGN_NOBNDRY")); -} + Field3D result = abs(field, "RGN_NOBNDRY"); -TEST_F(Field3DTest, Exp) { - Field3D field; + EXPECT_TRUE(IsFieldEqual(result, 31.0, "RGN_NOBNDRY")); + } - field = 2.5; - const BoutReal expected = 12.182493960703473; - EXPECT_TRUE(IsFieldEqual(exp(field), expected)); -} + TEST_F(Field3DTest, Exp) { + Field3D field; -TEST_F(Field3DTest, Log) { - Field3D field; + field = 2.5; + const BoutReal expected = 12.182493960703473; + EXPECT_TRUE(IsFieldEqual(exp(field), expected)); + } - field = 12.182493960703473; - const BoutReal expected = 2.5; - EXPECT_TRUE(IsFieldEqual(log(field), expected)); -} + TEST_F(Field3DTest, Log) { + Field3D field; -TEST_F(Field3DTest, LogExp) { - Field3D field; + field = 12.182493960703473; + const BoutReal expected = 2.5; + EXPECT_TRUE(IsFieldEqual(log(field), expected)); + } - field = 2.5; - const BoutReal expected = 2.5; - EXPECT_TRUE(IsFieldEqual(log(exp(field)), expected)); -} + TEST_F(Field3DTest, LogExp) { + Field3D field; -TEST_F(Field3DTest, Sin) { - Field3D field; + field = 2.5; + const BoutReal expected = 2.5; + EXPECT_TRUE(IsFieldEqual(log(exp(field)), expected)); + } - field = PI / 2.0; - EXPECT_TRUE(IsFieldEqual(sin(field), 1.0)); + TEST_F(Field3DTest, Sin) { + Field3D field; - field = PI; - EXPECT_TRUE(IsFieldEqual(sin(field), 0.0)); -} + field = PI / 2.0; + EXPECT_TRUE(IsFieldEqual(sin(field), 1.0)); -TEST_F(Field3DTest, Cos) { - Field3D field; + field = PI; + EXPECT_TRUE(IsFieldEqual(sin(field), 0.0)); + } - field = PI / 2.0; - EXPECT_TRUE(IsFieldEqual(cos(field), 0.0)); + TEST_F(Field3DTest, Cos) { + Field3D field; - field = PI; - EXPECT_TRUE(IsFieldEqual(cos(field), -1.0)); -} + field = PI / 2.0; + EXPECT_TRUE(IsFieldEqual(cos(field), 0.0)); -TEST_F(Field3DTest, Tan) { - Field3D field; + field = PI; + EXPECT_TRUE(IsFieldEqual(cos(field), -1.0)); + } - field = PI / 4.0; - EXPECT_TRUE(IsFieldEqual(tan(field), 1.0)); + TEST_F(Field3DTest, Tan) { + Field3D field; - field = PI; - EXPECT_TRUE(IsFieldEqual(tan(field), 0.0)); -} + field = PI / 4.0; + EXPECT_TRUE(IsFieldEqual(tan(field), 1.0)); -TEST_F(Field3DTest, Sinh) { - Field3D field; + field = PI; + EXPECT_TRUE(IsFieldEqual(tan(field), 0.0)); + } - field = 1.0; - const BoutReal expected = 1.1752011936438014; - EXPECT_TRUE(IsFieldEqual(sinh(field), expected)); + TEST_F(Field3DTest, Sinh) { + Field3D field; - field = -1.0; - EXPECT_TRUE(IsFieldEqual(sinh(field), -expected)); -} + field = 1.0; + const BoutReal expected = 1.1752011936438014; + EXPECT_TRUE(IsFieldEqual(sinh(field), expected)); -TEST_F(Field3DTest, Cosh) { - Field3D field; + field = -1.0; + EXPECT_TRUE(IsFieldEqual(sinh(field), -expected)); + } - field = 1.0; - const BoutReal expected = 1.5430806348152437; - EXPECT_TRUE(IsFieldEqual(cosh(field), expected)); + TEST_F(Field3DTest, Cosh) { + Field3D field; - field = -1.0; - EXPECT_TRUE(IsFieldEqual(cosh(field), expected)); -} + field = 1.0; + const BoutReal expected = 1.5430806348152437; + EXPECT_TRUE(IsFieldEqual(cosh(field), expected)); -TEST_F(Field3DTest, Tanh) { - Field3D field; + field = -1.0; + EXPECT_TRUE(IsFieldEqual(cosh(field), expected)); + } - field = 1.0; - const BoutReal expected = 0.761594155955764; - EXPECT_TRUE(IsFieldEqual(tanh(field), expected)); + TEST_F(Field3DTest, Tanh) { + Field3D field; - field = -1.0; - EXPECT_TRUE(IsFieldEqual(tanh(field), -expected)); -} + field = 1.0; + const BoutReal expected = 0.761594155955764; + EXPECT_TRUE(IsFieldEqual(tanh(field), expected)); -TEST_F(Field3DTest, Floor) { - Field3D field; + field = -1.0; + EXPECT_TRUE(IsFieldEqual(tanh(field), -expected)); + } - field = 50.0; - field(1, 1, 1) = 49.9; - field(2, 3, 4) = -20; + TEST_F(Field3DTest, Floor) { + Field3D field; - const BoutReal floor_value = 50.0; + field = 50.0; + field(1, 1, 1) = 49.9; + field(2, 3, 4) = -20; - EXPECT_TRUE(IsFieldEqual(floor(field, floor_value), floor_value)); -} + const BoutReal floor_value = 50.0; -TEST_F(Field3DTest, Min) { - Field3D field; + EXPECT_TRUE(IsFieldEqual(floor(field, floor_value), floor_value)); + } - field = 50.0; - field(0, 0, 0) = -99.0; - field(1, 1, 1) = 60.0; - field(1, 2, 2) = 40.0; - field(2, 4, 3) = 99.0; + TEST_F(Field3DTest, Min) { + Field3D field; - // min doesn't include guard cells - const BoutReal min_value = 40.0; + field = 50.0; + field(0, 0, 0) = -99.0; + field(1, 1, 1) = 60.0; + field(1, 2, 2) = 40.0; + field(2, 4, 3) = 99.0; - EXPECT_EQ(min(field, false), min_value); - EXPECT_EQ(min(field, false, "RGN_ALL"), -99.0); - EXPECT_EQ(min(field, true, "RGN_ALL"), -99.0); -} + // min doesn't include guard cells + const BoutReal min_value = 40.0; -TEST_F(Field3DTest, MinBinaryExpr) { - Field3D field; + EXPECT_EQ(min(field, false), min_value); + EXPECT_EQ(min(field, false, "RGN_ALL"), -99.0); + EXPECT_EQ(min(field, true, "RGN_ALL"), -99.0); + } - field = 50.0; - field(0, 0, 0) = -99.0; - field(1, 1, 1) = 60.0; - field(1, 2, 2) = 40.0; - field(2, 4, 3) = 99.0; + TEST_F(Field3DTest, MinBinaryExpr) { + Field3D field; - const auto expr = field / 2.0 - 5.0; + field = 50.0; + field(0, 0, 0) = -99.0; + field(1, 1, 1) = 60.0; + field(1, 2, 2) = 40.0; + field(2, 4, 3) = 99.0; - EXPECT_EQ(min(expr, false), 15.0); - EXPECT_EQ(min(expr, false, "RGN_ALL"), -54.5); -} + const auto expr = field / 2.0 - 5.0; -TEST_F(Field3DTest, Max) { - Field3D field; + EXPECT_EQ(min(expr, false), 15.0); + EXPECT_EQ(min(expr, false, "RGN_ALL"), -54.5); + } - field = 50.0; - field(0, 0, 0) = -99.0; - field(1, 1, 1) = 40.0; - field(1, 2, 2) = 60.0; - field(2, 4, 3) = 99.0; + TEST_F(Field3DTest, Max) { + Field3D field; - // max doesn't include guard cells - const BoutReal max_value = 60.0; + field = 50.0; + field(0, 0, 0) = -99.0; + field(1, 1, 1) = 40.0; + field(1, 2, 2) = 60.0; + field(2, 4, 3) = 99.0; - EXPECT_EQ(max(field, false), max_value); - EXPECT_EQ(max(field, false, "RGN_ALL"), 99.0); - EXPECT_EQ(max(field, true, "RGN_ALL"), 99.0); -} + // max doesn't include guard cells + const BoutReal max_value = 60.0; -TEST_F(Field3DTest, MaxBinaryExpr) { - Field3D field; + EXPECT_EQ(max(field, false), max_value); + EXPECT_EQ(max(field, false, "RGN_ALL"), 99.0); + EXPECT_EQ(max(field, true, "RGN_ALL"), 99.0); + } - field = 50.0; - field(0, 0, 0) = -99.0; - field(1, 1, 1) = 40.0; - field(1, 2, 2) = 60.0; - field(2, 4, 3) = 99.0; + TEST_F(Field3DTest, MaxBinaryExpr) { + Field3D field; - const auto expr = field / 2.0 - 5.0; + field = 50.0; + field(0, 0, 0) = -99.0; + field(1, 1, 1) = 40.0; + field(1, 2, 2) = 60.0; + field(2, 4, 3) = 99.0; - EXPECT_EQ(max(expr, false), 25.0); - EXPECT_EQ(max(expr, false, "RGN_ALL"), 44.5); -} + const auto expr = field / 2.0 - 5.0; -TEST_F(Field3DTest, Mean) { - Field3D field; + EXPECT_EQ(max(expr, false), 25.0); + EXPECT_EQ(max(expr, false, "RGN_ALL"), 44.5); + } - field = 50.0; - field(0, 0, 0) = 1.0; - field(1, 1, 1) = 40.0; - field(1, 2, 2) = 60.0; - field(2, 4, 3) = 109.0; + TEST_F(Field3DTest, Mean) { + Field3D field; - // mean doesn't include guard cells by default - const int npoints_all = nx * ny * nz; - const BoutReal mean_value_nobndry = 50.0; - const BoutReal mean_value_all = 50.0 + 10.0 / npoints_all; + field = 50.0; + field(0, 0, 0) = 1.0; + field(1, 1, 1) = 40.0; + field(1, 2, 2) = 60.0; + field(2, 4, 3) = 109.0; - EXPECT_EQ(mean(field, false), mean_value_nobndry); - EXPECT_EQ(mean(field, false, "RGN_ALL"), mean_value_all); - EXPECT_EQ(mean(field, true, "RGN_ALL"), mean_value_all); -} + // mean doesn't include guard cells by default + const int npoints_all = nx * ny * nz; + const BoutReal mean_value_nobndry = 50.0; + const BoutReal mean_value_all = 50.0 + 10.0 / npoints_all; -TEST_F(Field3DTest, MeanBinaryExpr) { - Field3D field; - - field = 50.0; - field(0, 0, 0) = 1.0; - field(1, 1, 1) = 40.0; - field(1, 2, 2) = 60.0; - field(2, 4, 3) = 109.0; + EXPECT_EQ(mean(field, false), mean_value_nobndry); + EXPECT_EQ(mean(field, false, "RGN_ALL"), mean_value_all); + EXPECT_EQ(mean(field, true, "RGN_ALL"), mean_value_all); + } - const int npoints_all = nx * ny * nz; - const BoutReal mean_value_nobndry = 103.0; - const BoutReal mean_value_all = 103.0 + 20.0 / npoints_all; - const auto expr = field * 2.0 + 3.0; + TEST_F(Field3DTest, MeanBinaryExpr) { + Field3D field; - EXPECT_EQ(mean(expr, false), mean_value_nobndry); - EXPECT_EQ(mean(expr, false, "RGN_ALL"), mean_value_all); -} + field = 50.0; + field(0, 0, 0) = 1.0; + field(1, 1, 1) = 40.0; + field(1, 2, 2) = 60.0; + field(2, 4, 3) = 109.0; -TEST_F(Field3DTest, DC) { - Field3D field; + const int npoints_all = nx * ny * nz; + const BoutReal mean_value_nobndry = 103.0; + const BoutReal mean_value_all = 103.0 + 20.0 / npoints_all; + const auto expr = field * 2.0 + 3.0; - field = 1.0; - for (const auto& i : field) { - field[i] = i.z(); + EXPECT_EQ(mean(expr, false), mean_value_nobndry); + EXPECT_EQ(mean(expr, false, "RGN_ALL"), mean_value_all); } - EXPECT_TRUE(IsFieldEqual(DC(field), 3.0)); -} + TEST_F(Field3DTest, DC) { + Field3D field; -TEST_F(Field3DTest, Swap) { - WithQuietOutput quiet{output_info}; + field = 1.0; + for (const auto& i : field) { + field[i] = i.z(); + } - // First field - Field3D first(1., mesh_staggered); + EXPECT_TRUE(IsFieldEqual(DC(field), 3.0)); + } - first.setLocation(CELL_XLOW); + TEST_F(Field3DTest, Swap) { + WithQuietOutput quiet{output_info}; - first.splitParallelSlices(); - first.yup() = 1.5; - first.ydown() = 0.5; + // First field + Field3D first(1., mesh_staggered); - ddt(first) = 1.1; + first.setLocation(CELL_XLOW); - // Mesh for second field - constexpr int second_nx = Field3DTest::nx + 2; - constexpr int second_ny = Field3DTest::ny + 2; - constexpr int second_nz = Field3DTest::nz + 2; + first.splitParallelSlices(); + first.yup() = 1.5; + first.ydown() = 0.5; - FakeMesh second_mesh{second_nx, second_ny, second_nz}; - second_mesh.setCoordinates(nullptr); - second_mesh.StaggerGrids = false; - second_mesh.createDefaultRegions(); + ddt(first) = 1.1; - // Second field - Field3D second(2., &second_mesh); + // Mesh for second field + constexpr int second_nx = Field3DTest::nx + 2; + constexpr int second_ny = Field3DTest::ny + 2; + constexpr int second_nz = Field3DTest::nz + 2; - second.splitParallelSlices(); - second.yup() = 2.2; - second.ydown() = 1.2; + FakeMesh second_mesh{second_nx, second_ny, second_nz}; + second_mesh.setCoordinates(nullptr); + second_mesh.StaggerGrids = false; + second_mesh.createDefaultRegions(); - ddt(second) = 2.4; + // Second field + Field3D second(2., &second_mesh); - // Basic sanity check - EXPECT_TRUE(IsFieldEqual(first, 1.0)); - EXPECT_TRUE(IsFieldEqual(second, 2.0)); + second.splitParallelSlices(); + second.yup() = 2.2; + second.ydown() = 1.2; - // swap is marked noexcept, so absolutely should not throw! - ASSERT_NO_THROW(swap(first, second)); + ddt(second) = 2.4; - // Values - EXPECT_TRUE(IsFieldEqual(first, 2.0)); - EXPECT_TRUE(IsFieldEqual(second, 1.0)); + // Basic sanity check + EXPECT_TRUE(IsFieldEqual(first, 1.0)); + EXPECT_TRUE(IsFieldEqual(second, 2.0)); - EXPECT_TRUE(IsFieldEqual(first.yup(), 2.2)); - EXPECT_TRUE(IsFieldEqual(first.ydown(), 1.2)); + // swap is marked noexcept, so absolutely should not throw! + ASSERT_NO_THROW(swap(first, second)); - EXPECT_TRUE(IsFieldEqual(second.yup(), 1.5)); - EXPECT_TRUE(IsFieldEqual(second.ydown(), 0.5)); + // Values + EXPECT_TRUE(IsFieldEqual(first, 2.0)); + EXPECT_TRUE(IsFieldEqual(second, 1.0)); - EXPECT_TRUE(IsFieldEqual(ddt(first), 2.4)); - EXPECT_TRUE(IsFieldEqual(ddt(second), 1.1)); + EXPECT_TRUE(IsFieldEqual(first.yup(), 2.2)); + EXPECT_TRUE(IsFieldEqual(first.ydown(), 1.2)); - // Mesh properties - EXPECT_EQ(first.getMesh(), &second_mesh); - EXPECT_EQ(second.getMesh(), mesh_staggered); + EXPECT_TRUE(IsFieldEqual(second.yup(), 1.5)); + EXPECT_TRUE(IsFieldEqual(second.ydown(), 0.5)); - EXPECT_EQ(first.getNx(), second_nx); - EXPECT_EQ(first.getNy(), second_ny); - EXPECT_EQ(first.getNz(), second_nz); + EXPECT_TRUE(IsFieldEqual(ddt(first), 2.4)); + EXPECT_TRUE(IsFieldEqual(ddt(second), 1.1)); - EXPECT_EQ(second.getNx(), Field3DTest::nx); - EXPECT_EQ(second.getNy(), Field3DTest::ny); - EXPECT_EQ(second.getNz(), Field3DTest::nz); + // Mesh properties + EXPECT_EQ(first.getMesh(), &second_mesh); + EXPECT_EQ(second.getMesh(), mesh_staggered); - EXPECT_EQ(first.getLocation(), CELL_CENTRE); - EXPECT_EQ(second.getLocation(), CELL_XLOW); + EXPECT_EQ(first.getNx(), second_nx); + EXPECT_EQ(first.getNy(), second_ny); + EXPECT_EQ(first.getNz(), second_nz); - // We don't check the boundaries, but the data is protected and - // there are no inquiry functions -} + EXPECT_EQ(second.getNx(), Field3DTest::nx); + EXPECT_EQ(second.getNy(), Field3DTest::ny); + EXPECT_EQ(second.getNz(), Field3DTest::nz); -TEST_F(Field3DTest, MoveCtor) { - // First field - Field3D first(1., mesh_staggered); + EXPECT_EQ(first.getLocation(), CELL_CENTRE); + EXPECT_EQ(second.getLocation(), CELL_XLOW); - first.setLocation(CELL_XLOW); + // We don't check the boundaries, but the data is protected and + // there are no inquiry functions + } - first.splitParallelSlices(); - first.yup() = 1.5; - first.ydown() = 0.5; + TEST_F(Field3DTest, MoveCtor) { + // First field + Field3D first(1., mesh_staggered); - ddt(first) = 1.1; + first.setLocation(CELL_XLOW); - // Second field - Field3D second{std::move(first)}; + first.splitParallelSlices(); + first.yup() = 1.5; + first.ydown() = 0.5; - // Values - EXPECT_TRUE(IsFieldEqual(second, 1.0)); + ddt(first) = 1.1; - EXPECT_TRUE(IsFieldEqual(second.yup(), 1.5)); - EXPECT_TRUE(IsFieldEqual(second.ydown(), 0.5)); + // Second field + Field3D second{std::move(first)}; - EXPECT_TRUE(IsFieldEqual(ddt(second), 1.1)); + // Values + EXPECT_TRUE(IsFieldEqual(second, 1.0)); - // Mesh properties - EXPECT_EQ(second.getMesh(), mesh_staggered); + EXPECT_TRUE(IsFieldEqual(second.yup(), 1.5)); + EXPECT_TRUE(IsFieldEqual(second.ydown(), 0.5)); - EXPECT_EQ(second.getNx(), Field3DTest::nx); - EXPECT_EQ(second.getNy(), Field3DTest::ny); - EXPECT_EQ(second.getNz(), Field3DTest::nz); + EXPECT_TRUE(IsFieldEqual(ddt(second), 1.1)); - EXPECT_EQ(second.getLocation(), CELL_XLOW); + // Mesh properties + EXPECT_EQ(second.getMesh(), mesh_staggered); - // We don't check the boundaries, but the data is protected and - // there are no inquiry functions -} + EXPECT_EQ(second.getNx(), Field3DTest::nx); + EXPECT_EQ(second.getNy(), Field3DTest::ny); + EXPECT_EQ(second.getNz(), Field3DTest::nz); -TEST_F(Field3DTest, FillField) { - Field3D f{mesh}; + EXPECT_EQ(second.getLocation(), CELL_XLOW); - fillField(f, {{{1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}}, + // We don't check the boundaries, but the data is protected and + // there are no inquiry functions + } - {{1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}}, + TEST_F(Field3DTest, FillField) { + Field3D f{mesh}; + + fillField(f, {{{1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}}, + + {{1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}}, + + {{1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}}}); + + EXPECT_TRUE(IsFieldEqual(f, 1.)); + + fillField(f, {{{0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}}, + + {{0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}}, + + {{0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}}}); + + Field3D g{mesh}; + g.allocate(); + BOUT_FOR_SERIAL(i, g.getRegion("RGN_ALL")) { g[i] = i.z(); } + + EXPECT_TRUE(IsFieldEqual(f, g)); + } - {{1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}}}); +#if BOUT_HAS_FFTW + namespace bout { + namespace testing { - EXPECT_TRUE(IsFieldEqual(f, 1.)); + // Amplitudes for the nth wavenumber + constexpr int k0{1}; + constexpr int k1{2}; + constexpr int k2{3}; - fillField(f, {{{0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}}, + const BoutReal box_size{TWOPI / Field3DTest::nz}; - {{0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}}, + // Helper function for the filter and lowpass tests + BoutReal zWaves(Field3D::ind_type& i) { + return 1.0 + std::sin(k0 * i.z() * box_size) + std::cos(k1 * i.z() * box_size) + + std::sin(k2 * i.z() * box_size); + } + } // namespace testing + } // namespace bout - {{0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}}}); + TEST_F(Field3DTest, Filter) { - Field3D g{mesh}; - g.allocate(); - BOUT_FOR_SERIAL(i, g.getRegion("RGN_ALL")) { g[i] = i.z(); } + using namespace bout::testing; - EXPECT_TRUE(IsFieldEqual(f, g)); -} + auto input = makeField(zWaves, bout::globals::mesh); -#if BOUT_HAS_FFTW -namespace bout { -namespace testing { + auto expected = makeField( + [&](Field3D::ind_type& i) { return std::cos(k1 * i.z() * box_size); }, + bout::globals::mesh); -// Amplitudes for the nth wavenumber -constexpr int k0{1}; -constexpr int k1{2}; -constexpr int k2{3}; + auto output = filter(input, 2); -const BoutReal box_size{TWOPI / Field3DTest::nz}; + EXPECT_TRUE(IsFieldEqual(output, expected)); + } -// Helper function for the filter and lowpass tests -BoutReal zWaves(Field3D::ind_type& i) { - return 1.0 + std::sin(k0 * i.z() * box_size) + std::cos(k1 * i.z() * box_size) - + std::sin(k2 * i.z() * box_size); -} -} // namespace testing -} // namespace bout + TEST_F(Field3DTest, LowPassOneArg) { -TEST_F(Field3DTest, Filter) { + using namespace bout::testing; - using namespace bout::testing; + auto input = makeField(zWaves, bout::globals::mesh); - auto input = makeField(zWaves, bout::globals::mesh); + auto expected = makeField( + [&](Field3D::ind_type& i) { + return 1.0 + std::sin(k0 * i.z() * box_size) + std::cos(k1 * i.z() * box_size); + }, + bout::globals::mesh); - auto expected = makeField( - [&](Field3D::ind_type& i) { return std::cos(k1 * i.z() * box_size); }, - bout::globals::mesh); + auto output = lowPass(input, 2); - auto output = filter(input, 2); + EXPECT_TRUE(IsFieldEqual(output, expected)); + } - EXPECT_TRUE(IsFieldEqual(output, expected)); -} + TEST_F(Field3DTest, LowPassOneArgNothing) { -TEST_F(Field3DTest, LowPassOneArg) { + using namespace bout::testing; - using namespace bout::testing; + auto input = makeField(zWaves, bout::globals::mesh); - auto input = makeField(zWaves, bout::globals::mesh); + auto output = lowPass(input, 20); - auto expected = makeField( - [&](Field3D::ind_type& i) { - return 1.0 + std::sin(k0 * i.z() * box_size) + std::cos(k1 * i.z() * box_size); - }, - bout::globals::mesh); + EXPECT_TRUE(IsFieldEqual(output, input)); + } - auto output = lowPass(input, 2); + TEST_F(Field3DTest, LowPassTwoArg) { - EXPECT_TRUE(IsFieldEqual(output, expected)); -} + using namespace bout::testing; -TEST_F(Field3DTest, LowPassOneArgNothing) { + auto input = makeField(zWaves, bout::globals::mesh); - using namespace bout::testing; + auto expected = makeField( + [&](Field3D::ind_type& i) { + return std::sin(k0 * i.z() * box_size) + std::cos(k1 * i.z() * box_size); + }, + bout::globals::mesh); - auto input = makeField(zWaves, bout::globals::mesh); + auto output = lowPass(input, 2, false); - auto output = lowPass(input, 20); + EXPECT_TRUE(IsFieldEqual(output, expected)); - EXPECT_TRUE(IsFieldEqual(output, input)); -} + // Check passing int still works + auto output2 = lowPass(input, 2, 0); -TEST_F(Field3DTest, LowPassTwoArg) { + EXPECT_TRUE(IsFieldEqual(output2, expected)); - using namespace bout::testing; + // Calling lowPass with an int that is not 0 or 1 is an error + EXPECT_THROW(lowPass(input, 2, -1), BoutException); + EXPECT_THROW(lowPass(input, 2, 2), BoutException); + } - auto input = makeField(zWaves, bout::globals::mesh); + TEST_F(Field3DTest, LowPassTwoArgKeepZonal) { - auto expected = makeField( - [&](Field3D::ind_type& i) { - return std::sin(k0 * i.z() * box_size) + std::cos(k1 * i.z() * box_size); - }, - bout::globals::mesh); + using namespace bout::testing; - auto output = lowPass(input, 2, false); + auto input = makeField(zWaves, bout::globals::mesh); - EXPECT_TRUE(IsFieldEqual(output, expected)); + auto expected = makeField( + [&](Field3D::ind_type& i) { + return 1.0 + std::sin(k0 * i.z() * box_size) + std::cos(k1 * i.z() * box_size); + }, + bout::globals::mesh); - // Check passing int still works - auto output2 = lowPass(input, 2, 0); + auto output = lowPass(input, 2, true); - EXPECT_TRUE(IsFieldEqual(output2, expected)); + EXPECT_TRUE(IsFieldEqual(output, expected)); - // Calling lowPass with an int that is not 0 or 1 is an error - EXPECT_THROW(lowPass(input, 2, -1), BoutException); - EXPECT_THROW(lowPass(input, 2, 2), BoutException); -} + // Check passing int still works + auto output2 = lowPass(input, 2, 1); -TEST_F(Field3DTest, LowPassTwoArgKeepZonal) { + EXPECT_TRUE(IsFieldEqual(output2, expected)); + } - using namespace bout::testing; + TEST_F(Field3DTest, LowPassTwoArgNothing) { - auto input = makeField(zWaves, bout::globals::mesh); + using namespace bout::testing; - auto expected = makeField( - [&](Field3D::ind_type& i) { - return 1.0 + std::sin(k0 * i.z() * box_size) + std::cos(k1 * i.z() * box_size); - }, - bout::globals::mesh); + auto input = makeField(zWaves, bout::globals::mesh); - auto output = lowPass(input, 2, true); + auto output = lowPass(input, 20, true); - EXPECT_TRUE(IsFieldEqual(output, expected)); + EXPECT_TRUE(IsFieldEqual(output, input)); + } +#endif - // Check passing int still works - auto output2 = lowPass(input, 2, 1); + TEST_F(Field3DTest, OperatorEqualsField3D) { + Field3D field; - EXPECT_TRUE(IsFieldEqual(output2, expected)); -} + // Create field with non-default arguments so we can check they get copied + // to 'field'. + // Note that Average z-direction type is not really allowed for Field3D, but + // we don't check anywhere at the moment. + Field3D field2{ + mesh_staggered, CELL_XLOW, {YDirectionType::Aligned, ZDirectionType::Average}}; -TEST_F(Field3DTest, LowPassTwoArgNothing) { + field = field2; - using namespace bout::testing; + EXPECT_TRUE(areFieldsCompatible(field, field2)); + EXPECT_EQ(field.getMesh(), field2.getMesh()); + EXPECT_EQ(field.getLocation(), field2.getLocation()); + EXPECT_EQ(field.getDirectionY(), field2.getDirectionY()); + EXPECT_EQ(field.getDirectionZ(), field2.getDirectionZ()); + } - auto input = makeField(zWaves, bout::globals::mesh); + TEST_F(Field3DTest, OperatorEqualsBinaryExprCopiesMetadata) { + Field3D source{ + mesh_staggered, CELL_XLOW, {YDirectionType::Aligned, ZDirectionType::Average}}; + source = 9.; - auto output = lowPass(input, 20, true); + Field3D target(mesh_staggered); + target = 0.; + target.splitParallelSlices(); - EXPECT_TRUE(IsFieldEqual(output, input)); -} -#endif + target = sqrt(source); -TEST_F(Field3DTest, OperatorEqualsField3D) { - Field3D field; + EXPECT_EQ(target.getMesh(), source.getMesh()); + EXPECT_EQ(target.getLocation(), source.getLocation()); + EXPECT_EQ(target.getDirectionY(), source.getDirectionY()); + EXPECT_EQ(target.getDirectionZ(), source.getDirectionZ()); + EXPECT_FALSE(target.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(target, 3.)); + } - // Create field with non-default arguments so we can check they get copied - // to 'field'. - // Note that Average z-direction type is not really allowed for Field3D, but - // we don't check anywhere at the moment. - Field3D field2{ - mesh_staggered, CELL_XLOW, {YDirectionType::Aligned, ZDirectionType::Average}}; + TEST_F(Field3DTest, EmptyFrom) { + // Create field with non-default arguments so we can check they get copied + // to 'field2'. + // Note that Average z-direction type is not really allowed for Field3D, but + // we don't check anywhere at the moment. + Field3D field{ + mesh_staggered, CELL_XLOW, {YDirectionType::Aligned, ZDirectionType::Average}}; + field = 5.; + + Field3D field2{emptyFrom(field)}; + EXPECT_EQ(field2.getMesh(), mesh_staggered); + EXPECT_EQ(field2.getLocation(), CELL_XLOW); + EXPECT_EQ(field2.getDirectionY(), YDirectionType::Aligned); + EXPECT_EQ(field2.getDirectionZ(), ZDirectionType::Average); + EXPECT_TRUE(field2.isAllocated()); + } - field = field2; + TEST_F(Field3DTest, ZeroFrom) { + // Create field with non-default arguments so we can check they get copied + // to 'field2'. + // Note that Average z-direction type is not really allowed for Field3D, but + // we don't check anywhere at the moment. + Field3D field{ + mesh_staggered, CELL_XLOW, {YDirectionType::Aligned, ZDirectionType::Average}}; + field = 5.; + + Field3D field2{zeroFrom(field)}; + EXPECT_EQ(field2.getMesh(), mesh_staggered); + EXPECT_EQ(field2.getLocation(), CELL_XLOW); + EXPECT_EQ(field2.getDirectionY(), YDirectionType::Aligned); + EXPECT_EQ(field2.getDirectionZ(), ZDirectionType::Average); + EXPECT_TRUE(field2.isAllocated()); + EXPECT_TRUE(IsFieldEqual(field2, 0.)); + } - EXPECT_TRUE(areFieldsCompatible(field, field2)); - EXPECT_EQ(field.getMesh(), field2.getMesh()); - EXPECT_EQ(field.getLocation(), field2.getLocation()); - EXPECT_EQ(field.getDirectionY(), field2.getDirectionY()); - EXPECT_EQ(field.getDirectionZ(), field2.getDirectionZ()); -} - -TEST_F(Field3DTest, OperatorEqualsBinaryExprCopiesMetadata) { - Field3D source{ - mesh_staggered, CELL_XLOW, {YDirectionType::Aligned, ZDirectionType::Average}}; - source = 9.; - - Field3D target(mesh_staggered); - target = 0.; - target.splitParallelSlices(); - - target = sqrt(source); - - EXPECT_EQ(target.getMesh(), source.getMesh()); - EXPECT_EQ(target.getLocation(), source.getLocation()); - EXPECT_EQ(target.getDirectionY(), source.getDirectionY()); - EXPECT_EQ(target.getDirectionZ(), source.getDirectionZ()); - EXPECT_FALSE(target.hasParallelSlices()); - EXPECT_TRUE(IsFieldEqual(target, 3.)); -} - -TEST_F(Field3DTest, EmptyFrom) { - // Create field with non-default arguments so we can check they get copied - // to 'field2'. - // Note that Average z-direction type is not really allowed for Field3D, but - // we don't check anywhere at the moment. - Field3D field{ - mesh_staggered, CELL_XLOW, {YDirectionType::Aligned, ZDirectionType::Average}}; - field = 5.; - - Field3D field2{emptyFrom(field)}; - EXPECT_EQ(field2.getMesh(), mesh_staggered); - EXPECT_EQ(field2.getLocation(), CELL_XLOW); - EXPECT_EQ(field2.getDirectionY(), YDirectionType::Aligned); - EXPECT_EQ(field2.getDirectionZ(), ZDirectionType::Average); - EXPECT_TRUE(field2.isAllocated()); -} - -TEST_F(Field3DTest, ZeroFrom) { - // Create field with non-default arguments so we can check they get copied - // to 'field2'. - // Note that Average z-direction type is not really allowed for Field3D, but - // we don't check anywhere at the moment. - Field3D field{ - mesh_staggered, CELL_XLOW, {YDirectionType::Aligned, ZDirectionType::Average}}; - field = 5.; - - Field3D field2{zeroFrom(field)}; - EXPECT_EQ(field2.getMesh(), mesh_staggered); - EXPECT_EQ(field2.getLocation(), CELL_XLOW); - EXPECT_EQ(field2.getDirectionY(), YDirectionType::Aligned); - EXPECT_EQ(field2.getDirectionZ(), ZDirectionType::Average); - EXPECT_TRUE(field2.isAllocated()); - EXPECT_TRUE(IsFieldEqual(field2, 0.)); -} - -TEST_F(Field3DTest, Field3DParallel) { - Field3DParallel field(1.0); - field = 1.0; + TEST_F(Field3DTest, Field3DParallel) { + Field3DParallel field(1.0); + field = 1.0; - Field3D field2 = field; + Field3D field2 = field; - auto& field3 = field.asField3D(); + auto& field3 = field.asField3D(); - field *= 2; + field *= 2; - EXPECT_TRUE(IsFieldEqual(field, 2.0)); - EXPECT_TRUE(IsFieldEqual(field2, 1.0)); - EXPECT_TRUE(IsFieldEqual(field3, 2.0)); + EXPECT_TRUE(IsFieldEqual(field, 2.0)); + EXPECT_TRUE(IsFieldEqual(field2, 1.0)); + EXPECT_TRUE(IsFieldEqual(field3, 2.0)); - field3.asField3DParallel() *= 3; + field3.asField3DParallel() *= 3; - EXPECT_TRUE(IsFieldEqual(field3, 6.0)); -} + EXPECT_TRUE(IsFieldEqual(field3, 6.0)); + } // Restore compiler warnings #pragma GCC diagnostic pop From c23f2d532b8093beab2fc37868ae420f1b7ac124 Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 9 Jul 2026 11:44:46 +0200 Subject: [PATCH 04/14] Allow FCI fixture --- tests/unit/fake_mesh_fixture.hxx | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) 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; From 356fd0c3e48c6a5b905e9565576bea11378cfed8 Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 9 Jul 2026 11:44:33 +0200 Subject: [PATCH 05/14] Move MockParallelTransform To make it reusable --- tests/unit/fake_mesh.hxx | 51 +++++++++++++++++++++++++ tests/unit/field/test_field_factory.cxx | 51 ------------------------- 2 files changed, 51 insertions(+), 51 deletions(-) diff --git a/tests/unit/fake_mesh.hxx b/tests/unit/fake_mesh.hxx index f957652a62..6559e554f2 100644 --- a/tests/unit/fake_mesh.hxx +++ b/tests/unit/fake_mesh.hxx @@ -366,3 +366,54 @@ 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() = 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; +}; 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}; From 8f66c008d013104f684a2c4bf7dd9f404e5e73a1 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 9 Jul 2026 22:12:40 -0700 Subject: [PATCH 06/14] test_field3d: Tidying, define FakeMeshFixtureFCI Implement some clang-tidy suggestions, use FakeMeshFixtureFCI for tests with FCI parallel transform. --- tests/unit/fake_mesh_fixture.hxx | 1 + tests/unit/field/test_field3d.cxx | 1091 ++++++++++++++--------------- 2 files changed, 528 insertions(+), 564 deletions(-) diff --git a/tests/unit/fake_mesh_fixture.hxx b/tests/unit/fake_mesh_fixture.hxx index 3c261a6121..92d8ee2a13 100644 --- a/tests/unit/fake_mesh_fixture.hxx +++ b/tests/unit/fake_mesh_fixture.hxx @@ -153,3 +153,4 @@ public: }; using FakeMeshFixture = FakeMeshFixture_tmpl<3, 5, 7>; +using FakeMeshFixtureFCI = FakeMeshFixture_tmpl<3, 5, 7, true>; diff --git a/tests/unit/field/test_field3d.cxx b/tests/unit/field/test_field3d.cxx index 89286a8f80..173991a31d 100644 --- a/tests/unit/field/test_field3d.cxx +++ b/tests/unit/field/test_field3d.cxx @@ -7,14 +7,16 @@ #include "gtest/gtest.h" #include "test_extras.hxx" +#include "bout/bout_types.hxx" #include "bout/boutexception.hxx" #include "bout/constants.hxx" +#include "bout/field2d.hxx" #include "bout/field3d.hxx" #include "bout/mesh.hxx" #include "bout/output.hxx" #include "bout/paralleltransform.hxx" +#include "bout/region.hxx" #include "bout/unused.hxx" -#include "bout/utils.hxx" #include #include @@ -27,35 +29,7 @@ using namespace bout::globals; // Reuse the "standard" fixture for FakeMesh using Field3DTest = FakeMeshFixture; -using Field3DTestFCI = FakeMeshFixture_tmpl<3, 5, 7, true>; - -class MockFciParallelTransform : public ParallelTransform { -public: - explicit MockFciParallelTransform(Mesh& mesh_in) : ParallelTransform(mesh_in) {} - - void calcParallelSlices(Field3D& f) override { - if (!f.hasParallelSlices()) { - f.splitParallelSlices(); - } - for (size_t i = 0; i < f.numberParallelSlices(); ++i) { - f.yup(i) = f; - f.ydown(i) = f; - } - } - - Field3D toFieldAligned(const Field3D& f, const std::string&) override { return f; } - FieldPerp toFieldAligned(const FieldPerp& f, const std::string&) override { return f; } - Field3D fromFieldAligned(const Field3D& f, const std::string&) override { return f; } - FieldPerp fromFieldAligned(const FieldPerp& f, const std::string&) override { - return f; - } - - bool canToFromFieldAligned() const override { return false; } - bool requiresTwistShift(bool, YDirectionType) override { return false; } - -protected: - void checkInputGrid() override {} -}; +using Field3DTestFCI = FakeMeshFixtureFCI; TEST_F(Field3DTest, Is3D) { Field3D field; @@ -104,14 +78,14 @@ TEST_F(Field3DTest, GetGridSizes) { } TEST_F(Field3DTest, CreateOnGivenMesh) { - int test_nx = Field3DTest::nx + 2; - int test_ny = Field3DTest::ny + 2; - int test_nz = Field3DTest::nz + 2; + const int test_nx = Field3DTest::nx + 2; + const int test_ny = Field3DTest::ny + 2; + const int test_nz = Field3DTest::nz + 2; FakeMesh fieldmesh{test_nx, test_ny, test_nz}; fieldmesh.setCoordinates(nullptr); - Field3D field{&fieldmesh}; + const Field3D field{&fieldmesh}; EXPECT_EQ(field.getNx(), test_nx); EXPECT_EQ(field.getNy(), test_ny); @@ -121,17 +95,17 @@ TEST_F(Field3DTest, CreateOnGivenMesh) { TEST_F(Field3DTest, CopyCheckFieldmesh) { WithQuietOutput quiet{output_info}; - int test_nx = Field3DTest::nx + 2; - int test_ny = Field3DTest::ny + 2; - int test_nz = Field3DTest::nz + 2; + const int test_nx = Field3DTest::nx + 2; + const int test_ny = Field3DTest::ny + 2; + const int test_nz = Field3DTest::nz + 2; FakeMesh fieldmesh{test_nx, test_ny, test_nz}; fieldmesh.setCoordinates(nullptr); fieldmesh.createDefaultRegions(); - Field3D field{0.0, &fieldmesh}; + const Field3D field{0.0, &fieldmesh}; - Field3D field2{field}; + const Field3D field2{field}; EXPECT_EQ(field2.getNx(), test_nx); EXPECT_EQ(field2.getNy(), test_ny); @@ -196,10 +170,10 @@ TEST_F(Field3DTest, CreateCopyOnNullMesh) { TEST_F(Field3DTest, TimeDeriv) { Field3D field; - auto deriv = field.timeDeriv(); + auto* deriv = field.timeDeriv(); EXPECT_NE(&field, deriv); - auto deriv2 = field.timeDeriv(); + auto* deriv2 = field.timeDeriv(); EXPECT_EQ(deriv, deriv2); EXPECT_EQ(&(ddt(field)), deriv); @@ -323,9 +297,9 @@ TEST_F(Field3DTest, ConstYnext) { const Field3D& field2 = field; - auto& yup = field2.ynext(1); + const auto& yup = field2.ynext(1); EXPECT_NE(&field2, &yup); - auto& ydown = field2.ynext(-1); + const auto& ydown = field2.ynext(-1); EXPECT_NE(&field2, &ydown); EXPECT_NE(&yup, &ydown); @@ -335,9 +309,9 @@ TEST_F(Field3DTest, ConstYnext) { } TEST_F(Field3DTest, GetGlobalMesh) { - Field3D field; + const Field3D field; - auto localmesh = field.getMesh(); + auto* localmesh = field.getMesh(); EXPECT_EQ(localmesh, mesh); } @@ -346,9 +320,9 @@ TEST_F(Field3DTest, GetLocalMesh) { FakeMesh myMesh{nx + 1, ny + 2, nz + 3}; myMesh.setCoordinates(nullptr); - Field3D field(&myMesh); + const Field3D field(&myMesh); - auto localmesh = field.getMesh(); + auto* localmesh = field.getMesh(); EXPECT_EQ(localmesh, &myMesh); } @@ -463,8 +437,9 @@ TEST_F(Field3DTest, IterateOverRegionInd3D_RGN_ALL) { // We use a set in case for some reason the iterator doesn't visit // each point in the order we expect - std::set> test_indices{{0, 0, 0}, {0, 0, 1}, {0, 1, 0}, {1, 0, 0}, - {0, 1, 1}, {1, 0, 1}, {1, 1, 0}, {1, 1, 1}}; + const std::set> test_indices{{0, 0, 0}, {0, 0, 1}, {0, 1, 0}, + {1, 0, 0}, {0, 1, 1}, {1, 0, 1}, + {1, 1, 0}, {1, 1, 1}}; const int num_sentinels = test_indices.size(); // Assign sentinel value to watch out for to our chosen points @@ -802,7 +777,7 @@ TEST_F(Field3DTest, IterateOverRGN_ZGUARDS) { test_indices.insert({1, 1, 1}); // This is the set of indices actually inside the region we want - std::set> region_indices; + const std::set> region_indices; const int num_sentinels = region_indices.size(); @@ -913,7 +888,7 @@ TEST_F(Field3DTest, IterateOver2DRGN_NOBNDRY) { EXPECT_EQ(i.z(), 0); } - EXPECT_EQ(sum, nx * ny - 2 * nx - 2 * (ny - 2)); + EXPECT_EQ(sum, (nx * ny) - (2 * nx) - (2 * (ny - 2))); } TEST_F(Field3DTest, IterateOver2DRGN_NOX) { @@ -931,7 +906,7 @@ TEST_F(Field3DTest, IterateOver2DRGN_NOX) { EXPECT_EQ(i.z(), 0); } - EXPECT_EQ(sum, nx * ny - 2 * ny); + EXPECT_EQ(sum, (nx * ny) - (2 * ny)); } TEST_F(Field3DTest, IterateOver2DRGN_NOY) { @@ -949,7 +924,7 @@ TEST_F(Field3DTest, IterateOver2DRGN_NOY) { EXPECT_EQ(i.z(), 0); } - EXPECT_EQ(sum, nx * ny - 2 * nx); + EXPECT_EQ(sum, (nx * ny) - (2 * nx)); } TEST_F(Field3DTest, Indexing) { @@ -981,7 +956,7 @@ TEST_F(Field3DTest, IndexingInd3D) { } } - Ind3D ind{(2 * ny + 2) * nz + 2}; + const Ind3D ind{(((2 * ny) + 2) * nz) + 2}; EXPECT_DOUBLE_EQ(field[ind], 6); } @@ -1001,7 +976,7 @@ TEST_F(Field3DTest, ConstIndexingInd3D) { const Field3D field2{field1}; - Ind3D ind{(2 * ny + 2) * nz + 2}; + const Ind3D ind{(((2 * ny) + 2) * nz) + 2}; EXPECT_DOUBLE_EQ(field2[ind], 6); } @@ -1012,7 +987,7 @@ TEST_F(Field3DTest, IndexingInd2D) { int ix = 1, iy = 2, iz = 3; field(ix, iy, iz) = sentinel; - Ind2D ind{iy + ny * ix, ny, 1}; + const Ind2D ind{iy + (ny * ix), ny, 1}; EXPECT_DOUBLE_EQ(field(ind, iz), sentinel); field(ind, iz) = -sentinel; EXPECT_DOUBLE_EQ(field(ix, iy, iz), -sentinel); @@ -1024,7 +999,7 @@ TEST_F(Field3DTest, ConstIndexingInd2D) { int ix = 1, iy = 2, iz = 3; field(ix, iy, iz) = sentinel; - Ind2D ind{iy + ny * ix, ny, 1}; + const Ind2D ind{iy + (ny * ix), ny, 1}; const Field3D field2{field}; EXPECT_DOUBLE_EQ(field2(ind, iz), sentinel); @@ -1036,7 +1011,7 @@ TEST_F(Field3DTest, IndexingIndPerp) { int ix = 1, iy = 2, iz = 3; field(ix, iy, iz) = sentinel; - IndPerp ind{iz + nz * ix, 1, nz}; + const IndPerp ind{iz + (nz * ix), 1, nz}; EXPECT_DOUBLE_EQ(field(ind, iy), sentinel); field(ind, iy) = -sentinel; EXPECT_DOUBLE_EQ(field(ix, iy, iz), -sentinel); @@ -1057,7 +1032,7 @@ TEST_F(Field3DTest, IndexingToZPointer) { for (int i = 0; i < nx; ++i) { for (int j = 0; j < ny; ++j) { - auto tmp = field(i, j); + auto* tmp = field(i, j); for (int k = 0; k < nz; ++k) { EXPECT_EQ(tmp[k], i + j + k); tmp[k] = -1.0; @@ -1086,7 +1061,7 @@ TEST_F(Field3DTest, ConstIndexingToZPointer) { for (int i = 0; i < nx; ++i) { for (int j = 0; j < ny; ++j) { - auto tmp = field(i, j); + const auto* tmp = field(i, j); for (int k = 0; k < nz; ++k) { EXPECT_EQ(tmp[k], 1.0); field2(i, j, k) = tmp[k]; @@ -1259,15 +1234,15 @@ TEST_F(Field3DTest, CreateFromBoutReal) { } TEST_F(Field3DTest, CreateFromField3D) { - Field3D field(99.0); - Field3D result(field); + const Field3D field(99.0); + const Field3D result(field); EXPECT_TRUE(IsFieldEqual(result, 99.0)); } TEST_F(Field3DTest, CreateFromField2D) { - Field2D field(99.0); - Field3D result(field); + const Field2D field(99.0); + const Field3D result(field); EXPECT_TRUE(IsFieldEqual(result, 99.0)); } @@ -1289,7 +1264,7 @@ TEST_F(Field3DTest, AssignFromInvalid) { TEST_F(Field3DTest, AssignFromField2D) { Field3D field; - Field2D field2(2.0); + const Field2D field2(2.0); field = field2; @@ -2053,10 +2028,7 @@ TEST_F(Field3DTest, Field3DParallelAssignmentFromLazyExprDiscardsSlicesWhenNotFc } #if CHECK >= 2 -TEST_F(Field3DTest, Field3DParallelAssignmentFromLazyExprRequiresField3DSlicesInFci) { - mesh->getCoordinates()->setParallelTransform( - bout::utils::make_unique(*mesh)); - +TEST_F(Field3DTestFCI, Field3DParallelAssignmentFromLazyExprRequiresField3DSlicesInFci) { Field3DParallel parallel; Field3D field; @@ -2069,10 +2041,7 @@ TEST_F(Field3DTest, Field3DParallelAssignmentFromLazyExprRequiresField3DSlicesIn } #endif -TEST_F(Field3DTest, Field3DParallelAssignmentFromLazyExprUsesField3DSlicesInFci) { - mesh->getCoordinates()->setParallelTransform( - bout::utils::make_unique(*mesh)); - +TEST_F(Field3DTestFCI, Field3DParallelAssignmentFromLazyExprUsesField3DSlicesInFci) { Field3DParallel parallel; Field3D field; @@ -2093,10 +2062,7 @@ TEST_F(Field3DTest, Field3DParallelAssignmentFromLazyExprUsesField3DSlicesInFci) EXPECT_TRUE(IsFieldEqual(result.ydown(), 15.0)); } -TEST_F(Field3DTest, Field3DParallelAssignmentFromField3DExprUsesSlicesInFci) { - mesh->getCoordinates()->setParallelTransform( - bout::utils::make_unique(*mesh)); - +TEST_F(Field3DTestFCI, Field3DParallelAssignmentFromField3DExprUsesSlicesInFci) { Field3D lhs; Field3D rhs; @@ -2118,10 +2084,7 @@ TEST_F(Field3DTest, Field3DParallelAssignmentFromField3DExprUsesSlicesInFci) { EXPECT_TRUE(IsFieldEqual(result.ydown(), 44.0)); } -TEST_F(Field3DTest, Field3DParallelAssignmentFromScalarExprUsesSlicesInFci) { - mesh->getCoordinates()->setParallelTransform( - bout::utils::make_unique(*mesh)); - +TEST_F(Field3DTestFCI, Field3DParallelAssignmentFromScalarExprUsesSlicesInFci) { Field3DParallel parallel; parallel = 2.0; parallel.yup() = 3.0; @@ -2133,649 +2096,649 @@ TEST_F(Field3DTest, Field3DParallelAssignmentFromScalarExprUsesSlicesInFci) { EXPECT_TRUE(IsFieldEqual(result, 3.0)); EXPECT_TRUE(IsFieldEqual(result.yup(), 4.0)); EXPECT_TRUE(IsFieldEqual(result.ydown(), 5.0)); +} - 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, MulField3DParallelPreservesParallelSlices) { + Field3D field; + EXPECT_TRUE(field.isFci()); - TEST_F(Field3DTest, Abs) { - Field3D field; + field = 2.0; + field.splitParallelSlices(); + field.yup() = 3.0; + field.ydown() = 4.0; + field.resetRegionParallel(); - field = -31.0; - EXPECT_TRUE(IsFieldEqual(abs(field), 31.0)); - } + Field3D rhs; + EXPECT_TRUE(rhs.isFci()); + rhs = 3.0; + rhs.splitParallelSlices(); + rhs.yup() = 4.0; + rhs.ydown() = 5.0; + rhs.resetRegionParallel(); - TEST_F(Field3DTest, AbsExpressionUsesAbsOp) { - Field3D field; + const Field3D prod = field * rhs; - field = -2.0; - const auto expr = field + 1.0; + EXPECT_FALSE(prod.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(prod, 6.0)); - EXPECT_TRUE( - (std::is_same_v, - BinaryExpr, - std::decay_t, bout::op::abs>>)); - EXPECT_TRUE(IsFieldEqual(abs(expr), 1.0)); - EXPECT_TRUE(IsFieldEqual(abs(expr, "RGN_ALL"), 1.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(Field3DTest, RegionLimitedExpressionConstructsField3D) { - Field3D field; +TEST_F(Field3DTest, Abs) { + Field3D field; - field = -31.0; + field = -31.0; + EXPECT_TRUE(IsFieldEqual(abs(field), 31.0)); +} - Field3D result = abs(field, "RGN_NOBNDRY"); +TEST_F(Field3DTest, AbsExpressionUsesAbsOp) { + Field3D field; - EXPECT_TRUE(IsFieldEqual(result, 31.0, "RGN_NOBNDRY")); - } + field = -2.0; + const auto expr = field + 1.0; - TEST_F(Field3DTest, Exp) { - Field3D field; + EXPECT_TRUE((std::is_same_v, + BinaryExpr, + std::decay_t, bout::op::abs>>)); + EXPECT_TRUE(IsFieldEqual(abs(expr), 1.0)); + EXPECT_TRUE(IsFieldEqual(abs(expr, "RGN_ALL"), 1.0)); +} - field = 2.5; - const BoutReal expected = 12.182493960703473; - EXPECT_TRUE(IsFieldEqual(exp(field), expected)); - } +TEST_F(Field3DTest, RegionLimitedExpressionConstructsField3D) { + Field3D field; - TEST_F(Field3DTest, Log) { - Field3D field; + field = -31.0; - field = 12.182493960703473; - const BoutReal expected = 2.5; - EXPECT_TRUE(IsFieldEqual(log(field), expected)); - } + Field3D result = abs(field, "RGN_NOBNDRY"); - TEST_F(Field3DTest, LogExp) { - Field3D field; + EXPECT_TRUE(IsFieldEqual(result, 31.0, "RGN_NOBNDRY")); +} - field = 2.5; - const BoutReal expected = 2.5; - EXPECT_TRUE(IsFieldEqual(log(exp(field)), expected)); - } +TEST_F(Field3DTest, Exp) { + Field3D field; - TEST_F(Field3DTest, Sin) { - Field3D field; + field = 2.5; + const BoutReal expected = 12.182493960703473; + EXPECT_TRUE(IsFieldEqual(exp(field), expected)); +} - field = PI / 2.0; - EXPECT_TRUE(IsFieldEqual(sin(field), 1.0)); +TEST_F(Field3DTest, Log) { + Field3D field; - field = PI; - EXPECT_TRUE(IsFieldEqual(sin(field), 0.0)); - } + field = 12.182493960703473; + const BoutReal expected = 2.5; + EXPECT_TRUE(IsFieldEqual(log(field), expected)); +} - TEST_F(Field3DTest, Cos) { - Field3D field; +TEST_F(Field3DTest, LogExp) { + Field3D field; - field = PI / 2.0; - EXPECT_TRUE(IsFieldEqual(cos(field), 0.0)); + field = 2.5; + const BoutReal expected = 2.5; + EXPECT_TRUE(IsFieldEqual(log(exp(field)), expected)); +} - field = PI; - EXPECT_TRUE(IsFieldEqual(cos(field), -1.0)); - } +TEST_F(Field3DTest, Sin) { + Field3D field; - TEST_F(Field3DTest, Tan) { - Field3D field; + field = PI / 2.0; + EXPECT_TRUE(IsFieldEqual(sin(field), 1.0)); - field = PI / 4.0; - EXPECT_TRUE(IsFieldEqual(tan(field), 1.0)); + field = PI; + EXPECT_TRUE(IsFieldEqual(sin(field), 0.0)); +} - field = PI; - EXPECT_TRUE(IsFieldEqual(tan(field), 0.0)); - } +TEST_F(Field3DTest, Cos) { + Field3D field; - TEST_F(Field3DTest, Sinh) { - Field3D field; + field = PI / 2.0; + EXPECT_TRUE(IsFieldEqual(cos(field), 0.0)); - field = 1.0; - const BoutReal expected = 1.1752011936438014; - EXPECT_TRUE(IsFieldEqual(sinh(field), expected)); + field = PI; + EXPECT_TRUE(IsFieldEqual(cos(field), -1.0)); +} - field = -1.0; - EXPECT_TRUE(IsFieldEqual(sinh(field), -expected)); - } +TEST_F(Field3DTest, Tan) { + Field3D field; - TEST_F(Field3DTest, Cosh) { - Field3D field; + field = PI / 4.0; + EXPECT_TRUE(IsFieldEqual(tan(field), 1.0)); - field = 1.0; - const BoutReal expected = 1.5430806348152437; - EXPECT_TRUE(IsFieldEqual(cosh(field), expected)); + field = PI; + EXPECT_TRUE(IsFieldEqual(tan(field), 0.0)); +} - field = -1.0; - EXPECT_TRUE(IsFieldEqual(cosh(field), expected)); - } +TEST_F(Field3DTest, Sinh) { + Field3D field; - TEST_F(Field3DTest, Tanh) { - Field3D field; + field = 1.0; + const BoutReal expected = 1.1752011936438014; + EXPECT_TRUE(IsFieldEqual(sinh(field), expected)); - field = 1.0; - const BoutReal expected = 0.761594155955764; - EXPECT_TRUE(IsFieldEqual(tanh(field), expected)); + field = -1.0; + EXPECT_TRUE(IsFieldEqual(sinh(field), -expected)); +} - field = -1.0; - EXPECT_TRUE(IsFieldEqual(tanh(field), -expected)); - } +TEST_F(Field3DTest, Cosh) { + Field3D field; - TEST_F(Field3DTest, Floor) { - Field3D field; + field = 1.0; + const BoutReal expected = 1.5430806348152437; + EXPECT_TRUE(IsFieldEqual(cosh(field), expected)); - field = 50.0; - field(1, 1, 1) = 49.9; - field(2, 3, 4) = -20; + field = -1.0; + EXPECT_TRUE(IsFieldEqual(cosh(field), expected)); +} - const BoutReal floor_value = 50.0; +TEST_F(Field3DTest, Tanh) { + Field3D field; - EXPECT_TRUE(IsFieldEqual(floor(field, floor_value), floor_value)); - } + field = 1.0; + const BoutReal expected = 0.761594155955764; + EXPECT_TRUE(IsFieldEqual(tanh(field), expected)); - TEST_F(Field3DTest, Min) { - Field3D field; + field = -1.0; + EXPECT_TRUE(IsFieldEqual(tanh(field), -expected)); +} - field = 50.0; - field(0, 0, 0) = -99.0; - field(1, 1, 1) = 60.0; - field(1, 2, 2) = 40.0; - field(2, 4, 3) = 99.0; +TEST_F(Field3DTest, Floor) { + Field3D field; - // min doesn't include guard cells - const BoutReal min_value = 40.0; + field = 50.0; + field(1, 1, 1) = 49.9; + field(2, 3, 4) = -20; - EXPECT_EQ(min(field, false), min_value); - EXPECT_EQ(min(field, false, "RGN_ALL"), -99.0); - EXPECT_EQ(min(field, true, "RGN_ALL"), -99.0); - } + const BoutReal floor_value = 50.0; - TEST_F(Field3DTest, MinBinaryExpr) { - Field3D field; + EXPECT_TRUE(IsFieldEqual(floor(field, floor_value), floor_value)); +} - field = 50.0; - field(0, 0, 0) = -99.0; - field(1, 1, 1) = 60.0; - field(1, 2, 2) = 40.0; - field(2, 4, 3) = 99.0; +TEST_F(Field3DTest, Min) { + Field3D field; - const auto expr = field / 2.0 - 5.0; + field = 50.0; + field(0, 0, 0) = -99.0; + field(1, 1, 1) = 60.0; + field(1, 2, 2) = 40.0; + field(2, 4, 3) = 99.0; - EXPECT_EQ(min(expr, false), 15.0); - EXPECT_EQ(min(expr, false, "RGN_ALL"), -54.5); - } + // min doesn't include guard cells + const BoutReal min_value = 40.0; - TEST_F(Field3DTest, Max) { - Field3D field; + EXPECT_EQ(min(field, false), min_value); + EXPECT_EQ(min(field, false, "RGN_ALL"), -99.0); + EXPECT_EQ(min(field, true, "RGN_ALL"), -99.0); +} - field = 50.0; - field(0, 0, 0) = -99.0; - field(1, 1, 1) = 40.0; - field(1, 2, 2) = 60.0; - field(2, 4, 3) = 99.0; +TEST_F(Field3DTest, MinBinaryExpr) { + Field3D field; - // max doesn't include guard cells - const BoutReal max_value = 60.0; + field = 50.0; + field(0, 0, 0) = -99.0; + field(1, 1, 1) = 60.0; + field(1, 2, 2) = 40.0; + field(2, 4, 3) = 99.0; - EXPECT_EQ(max(field, false), max_value); - EXPECT_EQ(max(field, false, "RGN_ALL"), 99.0); - EXPECT_EQ(max(field, true, "RGN_ALL"), 99.0); - } + const auto expr = field / 2.0 - 5.0; - TEST_F(Field3DTest, MaxBinaryExpr) { - Field3D field; + EXPECT_EQ(min(expr, false), 15.0); + EXPECT_EQ(min(expr, false, "RGN_ALL"), -54.5); +} - field = 50.0; - field(0, 0, 0) = -99.0; - field(1, 1, 1) = 40.0; - field(1, 2, 2) = 60.0; - field(2, 4, 3) = 99.0; +TEST_F(Field3DTest, Max) { + Field3D field; - const auto expr = field / 2.0 - 5.0; + field = 50.0; + field(0, 0, 0) = -99.0; + field(1, 1, 1) = 40.0; + field(1, 2, 2) = 60.0; + field(2, 4, 3) = 99.0; - EXPECT_EQ(max(expr, false), 25.0); - EXPECT_EQ(max(expr, false, "RGN_ALL"), 44.5); - } + // max doesn't include guard cells + const BoutReal max_value = 60.0; - TEST_F(Field3DTest, Mean) { - Field3D field; + EXPECT_EQ(max(field, false), max_value); + EXPECT_EQ(max(field, false, "RGN_ALL"), 99.0); + EXPECT_EQ(max(field, true, "RGN_ALL"), 99.0); +} - field = 50.0; - field(0, 0, 0) = 1.0; - field(1, 1, 1) = 40.0; - field(1, 2, 2) = 60.0; - field(2, 4, 3) = 109.0; +TEST_F(Field3DTest, MaxBinaryExpr) { + Field3D field; - // mean doesn't include guard cells by default - const int npoints_all = nx * ny * nz; - const BoutReal mean_value_nobndry = 50.0; - const BoutReal mean_value_all = 50.0 + 10.0 / npoints_all; + field = 50.0; + field(0, 0, 0) = -99.0; + field(1, 1, 1) = 40.0; + field(1, 2, 2) = 60.0; + field(2, 4, 3) = 99.0; - EXPECT_EQ(mean(field, false), mean_value_nobndry); - EXPECT_EQ(mean(field, false, "RGN_ALL"), mean_value_all); - EXPECT_EQ(mean(field, true, "RGN_ALL"), mean_value_all); - } + const auto expr = field / 2.0 - 5.0; - TEST_F(Field3DTest, MeanBinaryExpr) { - Field3D field; + EXPECT_EQ(max(expr, false), 25.0); + EXPECT_EQ(max(expr, false, "RGN_ALL"), 44.5); +} - field = 50.0; - field(0, 0, 0) = 1.0; - field(1, 1, 1) = 40.0; - field(1, 2, 2) = 60.0; - field(2, 4, 3) = 109.0; +TEST_F(Field3DTest, Mean) { + Field3D field; - const int npoints_all = nx * ny * nz; - const BoutReal mean_value_nobndry = 103.0; - const BoutReal mean_value_all = 103.0 + 20.0 / npoints_all; - const auto expr = field * 2.0 + 3.0; + field = 50.0; + field(0, 0, 0) = 1.0; + field(1, 1, 1) = 40.0; + field(1, 2, 2) = 60.0; + field(2, 4, 3) = 109.0; - EXPECT_EQ(mean(expr, false), mean_value_nobndry); - EXPECT_EQ(mean(expr, false, "RGN_ALL"), mean_value_all); - } + // mean doesn't include guard cells by default + const int npoints_all = nx * ny * nz; + const BoutReal mean_value_nobndry = 50.0; + const BoutReal mean_value_all = 50.0 + 10.0 / npoints_all; - TEST_F(Field3DTest, DC) { - Field3D field; + EXPECT_EQ(mean(field, false), mean_value_nobndry); + EXPECT_EQ(mean(field, false, "RGN_ALL"), mean_value_all); + EXPECT_EQ(mean(field, true, "RGN_ALL"), mean_value_all); +} - field = 1.0; - for (const auto& i : field) { - field[i] = i.z(); - } +TEST_F(Field3DTest, MeanBinaryExpr) { + Field3D field; + + field = 50.0; + field(0, 0, 0) = 1.0; + field(1, 1, 1) = 40.0; + field(1, 2, 2) = 60.0; + field(2, 4, 3) = 109.0; + + const int npoints_all = nx * ny * nz; + const BoutReal mean_value_nobndry = 103.0; + const BoutReal mean_value_all = 103.0 + 20.0 / npoints_all; + const auto expr = field * 2.0 + 3.0; + + EXPECT_EQ(mean(expr, false), mean_value_nobndry); + EXPECT_EQ(mean(expr, false, "RGN_ALL"), mean_value_all); +} + +TEST_F(Field3DTest, DC) { + Field3D field; - EXPECT_TRUE(IsFieldEqual(DC(field), 3.0)); + field = 1.0; + for (const auto& i : field) { + field[i] = i.z(); } - TEST_F(Field3DTest, Swap) { - WithQuietOutput quiet{output_info}; + EXPECT_TRUE(IsFieldEqual(DC(field), 3.0)); +} - // First field - Field3D first(1., mesh_staggered); +TEST_F(Field3DTest, Swap) { + WithQuietOutput quiet{output_info}; - first.setLocation(CELL_XLOW); + // First field + Field3D first(1., mesh_staggered); - first.splitParallelSlices(); - first.yup() = 1.5; - first.ydown() = 0.5; + first.setLocation(CELL_XLOW); - ddt(first) = 1.1; + first.splitParallelSlices(); + first.yup() = 1.5; + first.ydown() = 0.5; - // Mesh for second field - constexpr int second_nx = Field3DTest::nx + 2; - constexpr int second_ny = Field3DTest::ny + 2; - constexpr int second_nz = Field3DTest::nz + 2; + ddt(first) = 1.1; - FakeMesh second_mesh{second_nx, second_ny, second_nz}; - second_mesh.setCoordinates(nullptr); - second_mesh.StaggerGrids = false; - second_mesh.createDefaultRegions(); + // Mesh for second field + constexpr int second_nx = Field3DTest::nx + 2; + constexpr int second_ny = Field3DTest::ny + 2; + constexpr int second_nz = Field3DTest::nz + 2; - // Second field - Field3D second(2., &second_mesh); + FakeMesh second_mesh{second_nx, second_ny, second_nz}; + second_mesh.setCoordinates(nullptr); + second_mesh.StaggerGrids = false; + second_mesh.createDefaultRegions(); - second.splitParallelSlices(); - second.yup() = 2.2; - second.ydown() = 1.2; + // Second field + Field3D second(2., &second_mesh); - ddt(second) = 2.4; + second.splitParallelSlices(); + second.yup() = 2.2; + second.ydown() = 1.2; - // Basic sanity check - EXPECT_TRUE(IsFieldEqual(first, 1.0)); - EXPECT_TRUE(IsFieldEqual(second, 2.0)); + ddt(second) = 2.4; - // swap is marked noexcept, so absolutely should not throw! - ASSERT_NO_THROW(swap(first, second)); + // Basic sanity check + EXPECT_TRUE(IsFieldEqual(first, 1.0)); + EXPECT_TRUE(IsFieldEqual(second, 2.0)); - // Values - EXPECT_TRUE(IsFieldEqual(first, 2.0)); - EXPECT_TRUE(IsFieldEqual(second, 1.0)); + // swap is marked noexcept, so absolutely should not throw! + ASSERT_NO_THROW(swap(first, second)); - EXPECT_TRUE(IsFieldEqual(first.yup(), 2.2)); - EXPECT_TRUE(IsFieldEqual(first.ydown(), 1.2)); + // Values + EXPECT_TRUE(IsFieldEqual(first, 2.0)); + EXPECT_TRUE(IsFieldEqual(second, 1.0)); - EXPECT_TRUE(IsFieldEqual(second.yup(), 1.5)); - EXPECT_TRUE(IsFieldEqual(second.ydown(), 0.5)); + EXPECT_TRUE(IsFieldEqual(first.yup(), 2.2)); + EXPECT_TRUE(IsFieldEqual(first.ydown(), 1.2)); - EXPECT_TRUE(IsFieldEqual(ddt(first), 2.4)); - EXPECT_TRUE(IsFieldEqual(ddt(second), 1.1)); + EXPECT_TRUE(IsFieldEqual(second.yup(), 1.5)); + EXPECT_TRUE(IsFieldEqual(second.ydown(), 0.5)); - // Mesh properties - EXPECT_EQ(first.getMesh(), &second_mesh); - EXPECT_EQ(second.getMesh(), mesh_staggered); + EXPECT_TRUE(IsFieldEqual(ddt(first), 2.4)); + EXPECT_TRUE(IsFieldEqual(ddt(second), 1.1)); - EXPECT_EQ(first.getNx(), second_nx); - EXPECT_EQ(first.getNy(), second_ny); - EXPECT_EQ(first.getNz(), second_nz); + // Mesh properties + EXPECT_EQ(first.getMesh(), &second_mesh); + EXPECT_EQ(second.getMesh(), mesh_staggered); - EXPECT_EQ(second.getNx(), Field3DTest::nx); - EXPECT_EQ(second.getNy(), Field3DTest::ny); - EXPECT_EQ(second.getNz(), Field3DTest::nz); + EXPECT_EQ(first.getNx(), second_nx); + EXPECT_EQ(first.getNy(), second_ny); + EXPECT_EQ(first.getNz(), second_nz); - EXPECT_EQ(first.getLocation(), CELL_CENTRE); - EXPECT_EQ(second.getLocation(), CELL_XLOW); + EXPECT_EQ(second.getNx(), Field3DTest::nx); + EXPECT_EQ(second.getNy(), Field3DTest::ny); + EXPECT_EQ(second.getNz(), Field3DTest::nz); - // We don't check the boundaries, but the data is protected and - // there are no inquiry functions - } + EXPECT_EQ(first.getLocation(), CELL_CENTRE); + EXPECT_EQ(second.getLocation(), CELL_XLOW); - TEST_F(Field3DTest, MoveCtor) { - // First field - Field3D first(1., mesh_staggered); + // We don't check the boundaries, but the data is protected and + // there are no inquiry functions +} - first.setLocation(CELL_XLOW); +TEST_F(Field3DTest, MoveCtor) { + // First field + Field3D first(1., mesh_staggered); - first.splitParallelSlices(); - first.yup() = 1.5; - first.ydown() = 0.5; + first.setLocation(CELL_XLOW); - ddt(first) = 1.1; + first.splitParallelSlices(); + first.yup() = 1.5; + first.ydown() = 0.5; - // Second field - Field3D second{std::move(first)}; + ddt(first) = 1.1; - // Values - EXPECT_TRUE(IsFieldEqual(second, 1.0)); + // Second field + Field3D second{std::move(first)}; - EXPECT_TRUE(IsFieldEqual(second.yup(), 1.5)); - EXPECT_TRUE(IsFieldEqual(second.ydown(), 0.5)); + // Values + EXPECT_TRUE(IsFieldEqual(second, 1.0)); - EXPECT_TRUE(IsFieldEqual(ddt(second), 1.1)); + EXPECT_TRUE(IsFieldEqual(second.yup(), 1.5)); + EXPECT_TRUE(IsFieldEqual(second.ydown(), 0.5)); - // Mesh properties - EXPECT_EQ(second.getMesh(), mesh_staggered); + EXPECT_TRUE(IsFieldEqual(ddt(second), 1.1)); - EXPECT_EQ(second.getNx(), Field3DTest::nx); - EXPECT_EQ(second.getNy(), Field3DTest::ny); - EXPECT_EQ(second.getNz(), Field3DTest::nz); + // Mesh properties + EXPECT_EQ(second.getMesh(), mesh_staggered); - EXPECT_EQ(second.getLocation(), CELL_XLOW); + EXPECT_EQ(second.getNx(), Field3DTest::nx); + EXPECT_EQ(second.getNy(), Field3DTest::ny); + EXPECT_EQ(second.getNz(), Field3DTest::nz); - // We don't check the boundaries, but the data is protected and - // there are no inquiry functions - } + EXPECT_EQ(second.getLocation(), CELL_XLOW); - TEST_F(Field3DTest, FillField) { - Field3D f{mesh}; - - fillField(f, {{{1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}}, - - {{1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}}, - - {{1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}, - {1., 1., 1., 1., 1., 1., 1.}}}); - - EXPECT_TRUE(IsFieldEqual(f, 1.)); - - fillField(f, {{{0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}}, - - {{0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}}, - - {{0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}, - {0., 1., 2., 3., 4., 5., 6.}}}); - - Field3D g{mesh}; - g.allocate(); - BOUT_FOR_SERIAL(i, g.getRegion("RGN_ALL")) { g[i] = i.z(); } - - EXPECT_TRUE(IsFieldEqual(f, g)); - } + // We don't check the boundaries, but the data is protected and + // there are no inquiry functions +} -#if BOUT_HAS_FFTW - namespace bout { - namespace testing { +TEST_F(Field3DTest, FillField) { + Field3D f{mesh}; - // Amplitudes for the nth wavenumber - constexpr int k0{1}; - constexpr int k1{2}; - constexpr int k2{3}; + fillField(f, {{{1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}}, - const BoutReal box_size{TWOPI / Field3DTest::nz}; + {{1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}}, - // Helper function for the filter and lowpass tests - BoutReal zWaves(Field3D::ind_type& i) { - return 1.0 + std::sin(k0 * i.z() * box_size) + std::cos(k1 * i.z() * box_size) - + std::sin(k2 * i.z() * box_size); - } - } // namespace testing - } // namespace bout + {{1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}, + {1., 1., 1., 1., 1., 1., 1.}}}); - TEST_F(Field3DTest, Filter) { + EXPECT_TRUE(IsFieldEqual(f, 1.)); - using namespace bout::testing; + fillField(f, {{{0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}}, - auto input = makeField(zWaves, bout::globals::mesh); + {{0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}}, - auto expected = makeField( - [&](Field3D::ind_type& i) { return std::cos(k1 * i.z() * box_size); }, - bout::globals::mesh); + {{0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}, + {0., 1., 2., 3., 4., 5., 6.}}}); - auto output = filter(input, 2); + Field3D g{mesh}; + g.allocate(); + BOUT_FOR_SERIAL(i, g.getRegion("RGN_ALL")) { g[i] = i.z(); } - EXPECT_TRUE(IsFieldEqual(output, expected)); - } + EXPECT_TRUE(IsFieldEqual(f, g)); +} - TEST_F(Field3DTest, LowPassOneArg) { +#if BOUT_HAS_FFTW +namespace bout { +namespace testing { - using namespace bout::testing; +// Amplitudes for the nth wavenumber +constexpr int k0{1}; +constexpr int k1{2}; +constexpr int k2{3}; - auto input = makeField(zWaves, bout::globals::mesh); +const BoutReal box_size{TWOPI / Field3DTest::nz}; - auto expected = makeField( - [&](Field3D::ind_type& i) { - return 1.0 + std::sin(k0 * i.z() * box_size) + std::cos(k1 * i.z() * box_size); - }, - bout::globals::mesh); +// Helper function for the filter and lowpass tests +BoutReal zWaves(Field3D::ind_type& i) { + return 1.0 + std::sin(k0 * i.z() * box_size) + std::cos(k1 * i.z() * box_size) + + std::sin(k2 * i.z() * box_size); +} +} // namespace testing +} // namespace bout - auto output = lowPass(input, 2); +TEST_F(Field3DTest, Filter) { - EXPECT_TRUE(IsFieldEqual(output, expected)); - } + using namespace bout::testing; - TEST_F(Field3DTest, LowPassOneArgNothing) { + auto input = makeField(zWaves, bout::globals::mesh); - using namespace bout::testing; + auto expected = makeField( + [&](Field3D::ind_type& i) { return std::cos(k1 * i.z() * box_size); }, + bout::globals::mesh); - auto input = makeField(zWaves, bout::globals::mesh); + auto output = filter(input, 2); - auto output = lowPass(input, 20); + EXPECT_TRUE(IsFieldEqual(output, expected)); +} - EXPECT_TRUE(IsFieldEqual(output, input)); - } +TEST_F(Field3DTest, LowPassOneArg) { - TEST_F(Field3DTest, LowPassTwoArg) { + using namespace bout::testing; - using namespace bout::testing; + auto input = makeField(zWaves, bout::globals::mesh); - auto input = makeField(zWaves, bout::globals::mesh); + auto expected = makeField( + [&](Field3D::ind_type& i) { + return 1.0 + std::sin(k0 * i.z() * box_size) + std::cos(k1 * i.z() * box_size); + }, + bout::globals::mesh); - auto expected = makeField( - [&](Field3D::ind_type& i) { - return std::sin(k0 * i.z() * box_size) + std::cos(k1 * i.z() * box_size); - }, - bout::globals::mesh); + auto output = lowPass(input, 2); - auto output = lowPass(input, 2, false); + EXPECT_TRUE(IsFieldEqual(output, expected)); +} - EXPECT_TRUE(IsFieldEqual(output, expected)); +TEST_F(Field3DTest, LowPassOneArgNothing) { - // Check passing int still works - auto output2 = lowPass(input, 2, 0); + using namespace bout::testing; - EXPECT_TRUE(IsFieldEqual(output2, expected)); + auto input = makeField(zWaves, bout::globals::mesh); - // Calling lowPass with an int that is not 0 or 1 is an error - EXPECT_THROW(lowPass(input, 2, -1), BoutException); - EXPECT_THROW(lowPass(input, 2, 2), BoutException); - } + auto output = lowPass(input, 20); - TEST_F(Field3DTest, LowPassTwoArgKeepZonal) { + EXPECT_TRUE(IsFieldEqual(output, input)); +} - using namespace bout::testing; +TEST_F(Field3DTest, LowPassTwoArg) { - auto input = makeField(zWaves, bout::globals::mesh); + using namespace bout::testing; - auto expected = makeField( - [&](Field3D::ind_type& i) { - return 1.0 + std::sin(k0 * i.z() * box_size) + std::cos(k1 * i.z() * box_size); - }, - bout::globals::mesh); + auto input = makeField(zWaves, bout::globals::mesh); - auto output = lowPass(input, 2, true); + auto expected = makeField( + [&](Field3D::ind_type& i) { + return std::sin(k0 * i.z() * box_size) + std::cos(k1 * i.z() * box_size); + }, + bout::globals::mesh); - EXPECT_TRUE(IsFieldEqual(output, expected)); + auto output = lowPass(input, 2, false); - // Check passing int still works - auto output2 = lowPass(input, 2, 1); + EXPECT_TRUE(IsFieldEqual(output, expected)); - EXPECT_TRUE(IsFieldEqual(output2, expected)); - } + // Check passing int still works + auto output2 = lowPass(input, 2, 0); - TEST_F(Field3DTest, LowPassTwoArgNothing) { + EXPECT_TRUE(IsFieldEqual(output2, expected)); - using namespace bout::testing; + // Calling lowPass with an int that is not 0 or 1 is an error + EXPECT_THROW(lowPass(input, 2, -1), BoutException); + EXPECT_THROW(lowPass(input, 2, 2), BoutException); +} - auto input = makeField(zWaves, bout::globals::mesh); +TEST_F(Field3DTest, LowPassTwoArgKeepZonal) { - auto output = lowPass(input, 20, true); + using namespace bout::testing; - EXPECT_TRUE(IsFieldEqual(output, input)); - } -#endif + auto input = makeField(zWaves, bout::globals::mesh); - TEST_F(Field3DTest, OperatorEqualsField3D) { - Field3D field; + auto expected = makeField( + [&](Field3D::ind_type& i) { + return 1.0 + std::sin(k0 * i.z() * box_size) + std::cos(k1 * i.z() * box_size); + }, + bout::globals::mesh); - // Create field with non-default arguments so we can check they get copied - // to 'field'. - // Note that Average z-direction type is not really allowed for Field3D, but - // we don't check anywhere at the moment. - Field3D field2{ - mesh_staggered, CELL_XLOW, {YDirectionType::Aligned, ZDirectionType::Average}}; + auto output = lowPass(input, 2, true); - field = field2; + EXPECT_TRUE(IsFieldEqual(output, expected)); - EXPECT_TRUE(areFieldsCompatible(field, field2)); - EXPECT_EQ(field.getMesh(), field2.getMesh()); - EXPECT_EQ(field.getLocation(), field2.getLocation()); - EXPECT_EQ(field.getDirectionY(), field2.getDirectionY()); - EXPECT_EQ(field.getDirectionZ(), field2.getDirectionZ()); - } + // Check passing int still works + auto output2 = lowPass(input, 2, 1); - TEST_F(Field3DTest, OperatorEqualsBinaryExprCopiesMetadata) { - Field3D source{ - mesh_staggered, CELL_XLOW, {YDirectionType::Aligned, ZDirectionType::Average}}; - source = 9.; + EXPECT_TRUE(IsFieldEqual(output2, expected)); +} - Field3D target(mesh_staggered); - target = 0.; - target.splitParallelSlices(); +TEST_F(Field3DTest, LowPassTwoArgNothing) { - target = sqrt(source); + using namespace bout::testing; - EXPECT_EQ(target.getMesh(), source.getMesh()); - EXPECT_EQ(target.getLocation(), source.getLocation()); - EXPECT_EQ(target.getDirectionY(), source.getDirectionY()); - EXPECT_EQ(target.getDirectionZ(), source.getDirectionZ()); - EXPECT_FALSE(target.hasParallelSlices()); - EXPECT_TRUE(IsFieldEqual(target, 3.)); - } + auto input = makeField(zWaves, bout::globals::mesh); - TEST_F(Field3DTest, EmptyFrom) { - // Create field with non-default arguments so we can check they get copied - // to 'field2'. - // Note that Average z-direction type is not really allowed for Field3D, but - // we don't check anywhere at the moment. - Field3D field{ - mesh_staggered, CELL_XLOW, {YDirectionType::Aligned, ZDirectionType::Average}}; - field = 5.; - - Field3D field2{emptyFrom(field)}; - EXPECT_EQ(field2.getMesh(), mesh_staggered); - EXPECT_EQ(field2.getLocation(), CELL_XLOW); - EXPECT_EQ(field2.getDirectionY(), YDirectionType::Aligned); - EXPECT_EQ(field2.getDirectionZ(), ZDirectionType::Average); - EXPECT_TRUE(field2.isAllocated()); - } + auto output = lowPass(input, 20, true); - TEST_F(Field3DTest, ZeroFrom) { - // Create field with non-default arguments so we can check they get copied - // to 'field2'. - // Note that Average z-direction type is not really allowed for Field3D, but - // we don't check anywhere at the moment. - Field3D field{ - mesh_staggered, CELL_XLOW, {YDirectionType::Aligned, ZDirectionType::Average}}; - field = 5.; - - Field3D field2{zeroFrom(field)}; - EXPECT_EQ(field2.getMesh(), mesh_staggered); - EXPECT_EQ(field2.getLocation(), CELL_XLOW); - EXPECT_EQ(field2.getDirectionY(), YDirectionType::Aligned); - EXPECT_EQ(field2.getDirectionZ(), ZDirectionType::Average); - EXPECT_TRUE(field2.isAllocated()); - EXPECT_TRUE(IsFieldEqual(field2, 0.)); - } + EXPECT_TRUE(IsFieldEqual(output, input)); +} +#endif - TEST_F(Field3DTest, Field3DParallel) { - Field3DParallel field(1.0); - field = 1.0; +TEST_F(Field3DTest, OperatorEqualsField3D) { + Field3D field; - Field3D field2 = field; + // Create field with non-default arguments so we can check they get copied + // to 'field'. + // Note that Average z-direction type is not really allowed for Field3D, but + // we don't check anywhere at the moment. + Field3D field2{ + mesh_staggered, CELL_XLOW, {YDirectionType::Aligned, ZDirectionType::Average}}; - auto& field3 = field.asField3D(); + field = field2; - field *= 2; + EXPECT_TRUE(areFieldsCompatible(field, field2)); + EXPECT_EQ(field.getMesh(), field2.getMesh()); + EXPECT_EQ(field.getLocation(), field2.getLocation()); + EXPECT_EQ(field.getDirectionY(), field2.getDirectionY()); + EXPECT_EQ(field.getDirectionZ(), field2.getDirectionZ()); +} + +TEST_F(Field3DTest, OperatorEqualsBinaryExprCopiesMetadata) { + Field3D source{ + mesh_staggered, CELL_XLOW, {YDirectionType::Aligned, ZDirectionType::Average}}; + source = 9.; + + Field3D target(mesh_staggered); + target = 0.; + target.splitParallelSlices(); + + target = sqrt(source); + + EXPECT_EQ(target.getMesh(), source.getMesh()); + EXPECT_EQ(target.getLocation(), source.getLocation()); + EXPECT_EQ(target.getDirectionY(), source.getDirectionY()); + EXPECT_EQ(target.getDirectionZ(), source.getDirectionZ()); + EXPECT_FALSE(target.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(target, 3.)); +} + +TEST_F(Field3DTest, EmptyFrom) { + // Create field with non-default arguments so we can check they get copied + // to 'field2'. + // Note that Average z-direction type is not really allowed for Field3D, but + // we don't check anywhere at the moment. + Field3D field{ + mesh_staggered, CELL_XLOW, {YDirectionType::Aligned, ZDirectionType::Average}}; + field = 5.; + + Field3D field2{emptyFrom(field)}; + EXPECT_EQ(field2.getMesh(), mesh_staggered); + EXPECT_EQ(field2.getLocation(), CELL_XLOW); + EXPECT_EQ(field2.getDirectionY(), YDirectionType::Aligned); + EXPECT_EQ(field2.getDirectionZ(), ZDirectionType::Average); + EXPECT_TRUE(field2.isAllocated()); +} + +TEST_F(Field3DTest, ZeroFrom) { + // Create field with non-default arguments so we can check they get copied + // to 'field2'. + // Note that Average z-direction type is not really allowed for Field3D, but + // we don't check anywhere at the moment. + Field3D field{ + mesh_staggered, CELL_XLOW, {YDirectionType::Aligned, ZDirectionType::Average}}; + field = 5.; + + Field3D field2{zeroFrom(field)}; + EXPECT_EQ(field2.getMesh(), mesh_staggered); + EXPECT_EQ(field2.getLocation(), CELL_XLOW); + EXPECT_EQ(field2.getDirectionY(), YDirectionType::Aligned); + EXPECT_EQ(field2.getDirectionZ(), ZDirectionType::Average); + EXPECT_TRUE(field2.isAllocated()); + EXPECT_TRUE(IsFieldEqual(field2, 0.)); +} + +TEST_F(Field3DTest, Field3DParallel) { + Field3DParallel field(1.0); + field = 1.0; - EXPECT_TRUE(IsFieldEqual(field, 2.0)); - EXPECT_TRUE(IsFieldEqual(field2, 1.0)); - EXPECT_TRUE(IsFieldEqual(field3, 2.0)); + Field3D field2 = field; - field3.asField3DParallel() *= 3; + auto& field3 = field.asField3D(); - EXPECT_TRUE(IsFieldEqual(field3, 6.0)); - } + field *= 2; + + EXPECT_TRUE(IsFieldEqual(field, 2.0)); + EXPECT_TRUE(IsFieldEqual(field2, 1.0)); + EXPECT_TRUE(IsFieldEqual(field3, 2.0)); + + field3.asField3DParallel() *= 3; + + EXPECT_TRUE(IsFieldEqual(field3, 6.0)); +} // Restore compiler warnings #pragma GCC diagnostic pop From ceaaf980e6079f0370924c3ea91769b45df3c97d Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 9 Jul 2026 11:54:06 +0200 Subject: [PATCH 07/14] Functions of Field3DParallel preserve slices More unit tests --- include/bout/field.hxx | 119 ++++++++-------- include/bout/fieldops.hxx | 2 + src/field/field3d.cxx | 5 + tests/unit/field/test_field3d.cxx | 222 +++++++++++++++++++++++++++++- 4 files changed, 289 insertions(+), 59 deletions(-) diff --git a/include/bout/field.hxx b/include/bout/field.hxx index b39a82eb0b..b20e8915e9 100644 --- a/include/bout/field.hxx +++ b/include/bout/field.hxx @@ -608,6 +608,20 @@ class Field3DParallel; class FieldPerp; namespace bout::detail { +template +using UnaryFieldResult_t = + std::conditional_t, ::Field3DParallel>, ::Field3D, + std::decay_t>; + +template +std::optional getUnaryRegionID(const Mesh* mesh, const std::string& region_name) { + if constexpr (std::is_same_v, ::Field3D>) { + return bout::detail::getField3DRegionID(mesh, region_name); + } else { + return std::nullopt; + } +} + template std::optional getPerpYIndex(const T& value) { if constexpr (std::is_same_v, ::FieldPerp>) { @@ -641,35 +655,23 @@ std::optional getPerpYIndex(const BinaryExpr& expr) { }; \ template > \ inline auto name(const T& f, const std::string& rgn = "RGN_ALL") { \ - if constexpr (std::is_same_v) { \ - /* Check if the input is allocated */ \ - checkData(f); \ - /* Define and allocate the output result */ \ - T result{emptyFrom(f)}; \ - BOUT_FOR(d, result.getRegion(rgn)) { result[d] = func(f[d]); } \ - for (int i = 0; i < f.numberParallelSlices(); ++i) { \ - result.yup(i) = func(f.yup(i)); \ - result.ydown(i) = func(f.ydown(i)); \ - } \ - result.name = std::string(#name "(") + f.name + std::string(")"); \ - 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)}; \ - } \ + using ResT = bout::detail::UnaryFieldResult_t; \ + return BinaryExpr{ \ + static_cast(f), \ + static_cast(f), \ + bout::op::name{}, \ + f.getMesh(), \ + f.getLocation(), \ + f.getDirections(), \ + bout::detail::getUnaryRegionID(f.getMesh(), rgn), \ + f.getMesh()->template getRegion(rgn), \ + bout::detail::getPerpYIndex(f)}; \ } \ template \ inline auto name(const BinaryExpr& f) { \ - return BinaryExpr, BinaryExpr, \ - bout::op::name>{ \ + using UnaryResT = bout::detail::UnaryFieldResult_t; \ + return BinaryExpr, \ + BinaryExpr, bout::op::name>{ \ static_cast::View>(f), \ static_cast::View>(f), \ bout::op::name{}, \ @@ -682,7 +684,18 @@ std::optional getPerpYIndex(const BinaryExpr& expr) { } \ template \ inline auto name(const BinaryExpr& f, const std::string& rgn) { \ - return name(ResT{f}, rgn); \ + using UnaryResT = bout::detail::UnaryFieldResult_t; \ + return BinaryExpr, \ + BinaryExpr, bout::op::name>{ \ + static_cast::View>(f), \ + static_cast::View>(f), \ + bout::op::name{}, \ + f.getMesh(), \ + f.getLocation(), \ + f.getDirections(), \ + bout::detail::getUnaryRegionID(f.getMesh(), rgn), \ + f.getMesh()->template getRegion(rgn), \ + bout::detail::getPerpYIndex(f)}; \ } #endif @@ -698,36 +711,23 @@ struct Square { template > inline auto SQ(const T& f, const std::string& rgn = "RGN_ALL") { - if constexpr (std::is_same_v) { - checkData(f); - T result{emptyFrom(f)}; - if (f.hasParallelSlices() and !result.hasParallelSlices()) { - result.splitParallelSlices(); - } - BOUT_FOR(d, result.getRegion(rgn)) { result[d] = ::SQ(f[d]); } - for (size_t i = 0; i < f.numberParallelSlices(); ++i) { - result.yup(i) = SQ(f.yup(i), rgn); - result.ydown(i) = SQ(f.ydown(i), rgn); - } - result.name = std::string("SQ(") + f.name + std::string(")"); - checkData(result); - return result; - } else { - return BinaryExpr{static_cast(f), - static_cast(f), - bout::op::Square{}, - f.getMesh(), - f.getLocation(), - f.getDirections(), - std::nullopt, - f.getRegion(rgn), - bout::detail::getPerpYIndex(f)}; - } + using ResT = bout::detail::UnaryFieldResult_t; + return BinaryExpr{ + static_cast(f), + static_cast(f), + bout::op::Square{}, + f.getMesh(), + f.getLocation(), + f.getDirections(), + bout::detail::getUnaryRegionID(f.getMesh(), rgn), + f.getMesh()->template getRegion(rgn), + bout::detail::getPerpYIndex(f)}; } template inline auto SQ(const BinaryExpr& f) { - return BinaryExpr, BinaryExpr, + using UnaryResT = bout::detail::UnaryFieldResult_t; + return BinaryExpr, BinaryExpr, bout::op::Square>{ static_cast::View>(f), static_cast::View>(f), @@ -742,7 +742,18 @@ inline auto SQ(const BinaryExpr& f) { template inline auto SQ(const BinaryExpr& f, const std::string& rgn) { - return SQ(ResT{f}, rgn); + using UnaryResT = bout::detail::UnaryFieldResult_t; + return BinaryExpr, BinaryExpr, + bout::op::Square>{ + static_cast::View>(f), + static_cast::View>(f), + bout::op::Square{}, + f.getMesh(), + f.getLocation(), + f.getDirections(), + bout::detail::getUnaryRegionID(f.getMesh(), rgn), + f.getMesh()->template getRegion(rgn), + bout::detail::getPerpYIndex(f)}; } /// Square root of \p f over region \p rgn diff --git a/include/bout/fieldops.hxx b/include/bout/fieldops.hxx index a09734ad21..52d2c82a13 100644 --- a/include/bout/fieldops.hxx +++ b/include/bout/fieldops.hxx @@ -12,6 +12,7 @@ #include #include #include +#include #include #if BOUT_HAS_CUDA @@ -29,6 +30,7 @@ namespace bout::detail { // It is used because Mesh is an incomplete type so methods cannot be called // in the template functions in this header file. const Region& getField3DRegion(const Mesh* mesh, std::optional regionID); +size_t getField3DRegionID(const Mesh* mesh, const std::string& region_name); } // namespace bout::detail template diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 8700096e24..5da438dcdf 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -936,6 +936,11 @@ const Region& getField3DRegion(const Mesh* mesh, std::optional re return mesh->getRegion("RGN_ALL"); } +size_t getField3DRegionID(const Mesh* mesh, const std::string& region_name) { + ASSERT1(mesh != nullptr); + return mesh->getRegionID(region_name); +} + } // namespace bout::detail void swap(Field3D& first, Field3D& second) noexcept { diff --git a/tests/unit/field/test_field3d.cxx b/tests/unit/field/test_field3d.cxx index 173991a31d..38f51a8164 100644 --- a/tests/unit/field/test_field3d.cxx +++ b/tests/unit/field/test_field3d.cxx @@ -1977,11 +1977,16 @@ TEST_F(Field3DTest, SQField3DParallelPreservesParallelSlices) { const auto squared = SQ(field); - EXPECT_TRUE((std::is_same_v, Field3DParallel>)); - EXPECT_TRUE(squared.hasParallelSlices()); - EXPECT_TRUE(IsFieldEqual(squared, 4.0)); - EXPECT_TRUE(IsFieldEqual(squared.yup(), 9.0)); - EXPECT_TRUE(IsFieldEqual(squared.ydown(), 16.0)); + EXPECT_TRUE((std::is_same_v< + std::decay_t, + BinaryExpr>)); + + Field3DParallel result{squared}; + + EXPECT_FALSE(result.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(result, 4.0)); + EXPECT_TRUE(IsFieldEqual(result.yup(), 9.0)); + EXPECT_TRUE(IsFieldEqual(result.ydown(), 16.0)); } TEST_F(Field3DTest, Field3DParallelArithmeticReturnsLazyExpr) { @@ -2098,6 +2103,127 @@ TEST_F(Field3DTestFCI, Field3DParallelAssignmentFromScalarExprUsesSlicesInFci) { EXPECT_TRUE(IsFieldEqual(result.ydown(), 5.0)); } +TEST_F(Field3DTestFCI, SqrtField3DParallelPreservesParallelSlices) { + Field3DParallel field; + + field = 4.0; + field.splitParallelSlices(); + field.yup() = 9.0; + field.ydown() = 16.0; + field.resetRegionParallel(); + + { + const auto expr = sqrt(field); + + EXPECT_TRUE((std::is_same_v< + std::decay_t, + BinaryExpr>)); + + const Field3DParallel res{expr}; + + 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, SqrtField3DParallelExprPreservesParallelSlices) { + Field3DParallel lhs; + Field3D rhs; + + lhs = 3.0; + lhs.yup() = 8.0; + lhs.ydown() = 15.0; + + rhs = 1.0; + rhs.splitParallelSlices(); + rhs.yup() = 1.0; + rhs.ydown() = 1.0; + + Field3DParallel result{sqrt(lhs + rhs)}; + + EXPECT_TRUE(result.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(result, 2.0)); + EXPECT_TRUE(IsFieldEqual(result.yup(), 3.0)); + EXPECT_TRUE(IsFieldEqual(result.ydown(), 4.0)); +} + +TEST_F(Field3DTestFCI, SqrtField3DParallelRegionExprPreservesParallelSlices) { + Field3DParallel lhs; + Field3D rhs; + + lhs = 3.0; + lhs.yup() = 8.0; + lhs.ydown() = 15.0; + + rhs = 1.0; + rhs.splitParallelSlices(); + rhs.yup() = 1.0; + rhs.ydown() = 1.0; + + const auto base = lhs + rhs; + const auto expr = sqrt(base, "RGN_ALL"); + + EXPECT_TRUE((std::is_same_v, + BinaryExpr, + std::decay_t, bout::op::sqrt>>)); + + Field3DParallel result{expr}; + + EXPECT_TRUE(result.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(result, 2.0)); + EXPECT_TRUE(IsFieldEqual(result.yup(), 3.0)); + EXPECT_TRUE(IsFieldEqual(result.ydown(), 4.0)); +} + +TEST_F(Field3DTestFCI, SQField3DParallelPreservesParallelSlices) { + Field3DParallel field; + + field = 2.0; + field.yup() = 3.0; + field.ydown() = 4.0; + + const auto expr = SQ(field); + + EXPECT_TRUE((std::is_same_v< + std::decay_t, + BinaryExpr>)); + + Field3DParallel result{expr}; + + EXPECT_TRUE(result.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(result, 4.0)); + EXPECT_TRUE(IsFieldEqual(result.yup(), 9.0)); + EXPECT_TRUE(IsFieldEqual(result.ydown(), 16.0)); +} + +TEST_F(Field3DTestFCI, SQField3DParallelExprPreservesParallelSlices) { + Field3DParallel lhs; + Field3D rhs; + + lhs = 2.0; + lhs.yup() = 3.0; + lhs.ydown() = 4.0; + + rhs = 1.0; + rhs.splitParallelSlices(); + rhs.yup() = 2.0; + rhs.ydown() = 3.0; + + Field3DParallel result{SQ(lhs + rhs)}; + + EXPECT_TRUE(result.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(result, 9.0)); + EXPECT_TRUE(IsFieldEqual(result.yup(), 25.0)); + EXPECT_TRUE(IsFieldEqual(result.ydown(), 49.0)); +} + TEST_F(Field3DTestFCI, MulField3DParallelPreservesParallelSlices) { Field3D field; EXPECT_TRUE(field.isFci()); @@ -2129,6 +2255,92 @@ TEST_F(Field3DTestFCI, MulField3DParallelPreservesParallelSlices) { 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; From b2f52cde587b5322e38213285ccefa88e9ea150e Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 10 Jul 2026 08:24:28 -0700 Subject: [PATCH 08/14] Field3DParallel unit test If the parallel transform is not FCI then parallel slices are not calculated on assignment to Field3DParallel. --- tests/unit/field/test_field3d.cxx | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/tests/unit/field/test_field3d.cxx b/tests/unit/field/test_field3d.cxx index 38f51a8164..eee92c806e 100644 --- a/tests/unit/field/test_field3d.cxx +++ b/tests/unit/field/test_field3d.cxx @@ -1967,7 +1967,7 @@ TEST_F(Field3DTest, SQExpressionUsesSquareOp) { EXPECT_TRUE(IsFieldEqual(SQ(expr), 9.0)); } -TEST_F(Field3DTest, SQField3DParallelPreservesParallelSlices) { +TEST_F(Field3DTest, SQField3DParallelNoFCIDropsParallelSlices) { Field3DParallel field; field = 2.0; @@ -1983,10 +1983,9 @@ TEST_F(Field3DTest, SQField3DParallelPreservesParallelSlices) { Field3DParallel result{squared}; + // Not FCI so parallel slices are not calculated EXPECT_FALSE(result.hasParallelSlices()); EXPECT_TRUE(IsFieldEqual(result, 4.0)); - EXPECT_TRUE(IsFieldEqual(result.yup(), 9.0)); - EXPECT_TRUE(IsFieldEqual(result.ydown(), 16.0)); } TEST_F(Field3DTest, Field3DParallelArithmeticReturnsLazyExpr) { From 858c5c1961a8970c1b3a4e8443ddc693b1d82084 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 10 Jul 2026 09:59:49 -0700 Subject: [PATCH 09/14] Lazy pow function on fields --- include/bout/field.hxx | 232 ++++++++++++++++++++++++++---- include/bout/field3d.hxx | 1 - include/bout/fieldperp.hxx | 5 + include/bout/fv_ops.hxx | 14 +- src/field/field3d.cxx | 16 --- src/field/fieldperp.cxx | 35 +++++ tests/unit/field/test_field2d.cxx | 13 ++ tests/unit/field/test_field3d.cxx | 74 ++++++++++ 8 files changed, 340 insertions(+), 50 deletions(-) diff --git a/include/bout/field.hxx b/include/bout/field.hxx index b20e8915e9..05bccf250c 100644 --- a/include/bout/field.hxx +++ b/include/bout/field.hxx @@ -539,53 +539,233 @@ inline BoutReal mean(const BinaryExpr& f, bool allpe = false, return bout::reduce::Mean::finalize(state); } +namespace bout::op { +struct Pow { + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal operator()(int idx, const LView& L, + const RView& R) const { + return ::pow(L(idx), R(idx)); + } + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal operator()(BoutReal a, BoutReal b) const { + return ::pow(a, b); + } +}; +}; // namespace bout::op + +namespace bout::detail { +template +std::optional getPerpYIndex(const T& value); + +template +std::optional getPerpYIndex(const BinaryExpr& expr); + +template +std::optional getPowRegionID(const Mesh* mesh, const std::string& region_name) { + if constexpr (std::is_same_v) { + return bout::detail::getField3DRegionID(mesh, region_name); + } else { + return std::nullopt; + } +} + +template +auto makePowExpr(const LView& lhs_view, const RView& rhs_view, Mesh* mesh, + CELL_LOC location, DirectionTypes directions, + std::optional regionID, const Region& region, + std::optional yindex = std::nullopt) { + return BinaryExpr{lhs_view, rhs_view, bout::op::Pow{}, + mesh, location, directions, + regionID, region, yindex}; +} +} // namespace bout::detail + /// Exponent: pow(lhs, lhs) is \p lhs raised to the power of \p rhs /// /// This loops over the entire domain, including guard/boundary cells by /// default (can be changed using the \p rgn argument) /// If CHECK >= 3 then the result will be checked for non-finite numbers -template > -T pow(const T& lhs, const T& rhs, const std::string& rgn = "RGN_ALL") { +template +std::enable_if_t && is_expr_field2d_v, + BinaryExpr> +pow(const L& lhs, const R& rhs) { + ASSERT1_EXPR_COMPATIBLE(lhs, rhs); + return bout::detail::makePowExpr( + static_cast(lhs), static_cast(rhs), + lhs.getMesh(), lhs.getLocation(), lhs.getDirections(), std::nullopt, + lhs.getMesh()->getRegion2D("RGN_ALL")); +} - ASSERT1(areFieldsCompatible(lhs, rhs)); +template +std::enable_if_t && is_expr_field2d_v, + BinaryExpr> +pow(const L& lhs, const R& rhs, const std::string& rgn) { + ASSERT1_EXPR_COMPATIBLE(lhs, rhs); + return bout::detail::makePowExpr( + static_cast(lhs), static_cast(rhs), + lhs.getMesh(), lhs.getLocation(), lhs.getDirections(), std::nullopt, + lhs.getMesh()->getRegion2D(rgn)); +} + +template +std::enable_if_t && is_expr_field3d_v, + BinaryExpr> +pow(const L& lhs, const R& rhs) { + ASSERT1_EXPR_COMPATIBLE(lhs, rhs); + auto regionID = lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID()); + return bout::detail::makePowExpr( + static_cast(lhs), static_cast(rhs), + lhs.getMesh(), lhs.getLocation(), lhs.getDirections(), regionID, + (regionID.has_value() ? lhs.getMesh()->getRegion(regionID.value()) + : lhs.getMesh()->getRegion("RGN_ALL")), + bout::detail::getPerpYIndex(lhs)); +} - T result{emptyFrom(lhs)}; +template +std::enable_if_t && is_expr_field3d_v, + BinaryExpr> +pow(const L& lhs, const R& rhs, const std::string& rgn) { + ASSERT1_EXPR_COMPATIBLE(lhs, rhs); + return bout::detail::makePowExpr( + static_cast(lhs), static_cast(rhs), + lhs.getMesh(), lhs.getLocation(), lhs.getDirections(), + bout::detail::getPowRegionID(lhs.getMesh(), rgn), + lhs.getMesh()->getRegion(rgn), bout::detail::getPerpYIndex(lhs)); +} - BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs[i], rhs[i]); } +template +std::enable_if_t && is_expr_field2d_v, + BinaryExpr> +pow(const L& lhs, const R& rhs) { + ASSERT1_EXPR_COMPATIBLE(lhs, rhs); + int mesh_nz = lhs.getMesh()->LocalNz; + return bout::detail::makePowExpr( + static_cast(lhs), + static_cast(rhs).setScale(1, mesh_nz), lhs.getMesh(), + lhs.getLocation(), lhs.getDirections(), lhs.getRegionID(), + lhs.getMesh()->getRegion("RGN_ALL"), bout::detail::getPerpYIndex(lhs)); +} - checkData(result); - return result; +template +std::enable_if_t && is_expr_field2d_v, + BinaryExpr> +pow(const L& lhs, const R& rhs, const std::string& rgn) { + ASSERT1_EXPR_COMPATIBLE(lhs, rhs); + int mesh_nz = lhs.getMesh()->LocalNz; + return bout::detail::makePowExpr( + static_cast(lhs), + static_cast(rhs).setScale(1, mesh_nz), lhs.getMesh(), + lhs.getLocation(), lhs.getDirections(), + bout::detail::getPowRegionID(lhs.getMesh(), rgn), + lhs.getMesh()->getRegion(rgn), bout::detail::getPerpYIndex(lhs)); } -template > -T pow(const T& lhs, BoutReal rhs, const std::string& rgn = "RGN_ALL") { +template +std::enable_if_t && is_expr_field3d_v, + BinaryExpr> +pow(const L& lhs, const R& rhs) { + ASSERT1_EXPR_COMPATIBLE(lhs, rhs); + int mesh_nz = rhs.getMesh()->LocalNz; + return bout::detail::makePowExpr( + static_cast(lhs).setScale(1, mesh_nz), + static_cast(rhs), rhs.getMesh(), rhs.getLocation(), + rhs.getDirections(), rhs.getRegionID(), rhs.getMesh()->getRegion("RGN_ALL"), + bout::detail::getPerpYIndex(rhs)); +} - // Check if the inputs are allocated - checkData(lhs); - checkData(rhs); +template +std::enable_if_t && is_expr_field3d_v, + BinaryExpr> +pow(const L& lhs, const R& rhs, const std::string& rgn) { + ASSERT1_EXPR_COMPATIBLE(lhs, rhs); + int mesh_nz = rhs.getMesh()->LocalNz; + return bout::detail::makePowExpr( + static_cast(lhs).setScale(1, mesh_nz), + static_cast(rhs), rhs.getMesh(), rhs.getLocation(), + rhs.getDirections(), bout::detail::getPowRegionID(rhs.getMesh(), rgn), + rhs.getMesh()->getRegion(rgn), bout::detail::getPerpYIndex(rhs)); +} - T result{emptyFrom(lhs)}; +template +std::enable_if_t && is_expr_constant_v, + BinaryExpr, bout::op::Pow>> +pow(const L& lhs, R rhs) { + return bout::detail::makePowExpr>( + static_cast(lhs), static_cast::View>(rhs), + lhs.getMesh(), lhs.getLocation(), lhs.getDirections(), std::nullopt, + lhs.getMesh()->getRegion2D("RGN_ALL")); +} - BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs[i], rhs); } +template +std::enable_if_t && is_expr_constant_v, + BinaryExpr, bout::op::Pow>> +pow(const L& lhs, R rhs, const std::string& rgn) { + return bout::detail::makePowExpr>( + static_cast(lhs), static_cast::View>(rhs), + lhs.getMesh(), lhs.getLocation(), lhs.getDirections(), std::nullopt, + lhs.getMesh()->getRegion2D(rgn)); +} - checkData(result); - return result; +template +std::enable_if_t && is_expr_field2d_v, + BinaryExpr, R, bout::op::Pow>> +pow(L lhs, const R& rhs) { + return bout::detail::makePowExpr, R>( + static_cast::View>(lhs), static_cast(rhs), + rhs.getMesh(), rhs.getLocation(), rhs.getDirections(), std::nullopt, + rhs.getMesh()->getRegion2D("RGN_ALL")); } -template > -T pow(BoutReal lhs, const T& rhs, const std::string& rgn = "RGN_ALL") { +template +std::enable_if_t && is_expr_field2d_v, + BinaryExpr, R, bout::op::Pow>> +pow(L lhs, const R& rhs, const std::string& rgn) { + return bout::detail::makePowExpr, R>( + static_cast::View>(lhs), static_cast(rhs), + rhs.getMesh(), rhs.getLocation(), rhs.getDirections(), std::nullopt, + rhs.getMesh()->getRegion2D(rgn)); +} - // Check if the inputs are allocated - checkData(lhs); - checkData(rhs); +template +std::enable_if_t && is_expr_constant_v, + BinaryExpr, bout::op::Pow>> +pow(const L& lhs, R rhs) { + return bout::detail::makePowExpr>( + static_cast(lhs), static_cast::View>(rhs), + lhs.getMesh(), lhs.getLocation(), lhs.getDirections(), lhs.getRegionID(), + lhs.getMesh()->getRegion("RGN_ALL"), bout::detail::getPerpYIndex(lhs)); +} - // Define and allocate the output result - T result{emptyFrom(rhs)}; +template +std::enable_if_t && is_expr_constant_v, + BinaryExpr, bout::op::Pow>> +pow(const L& lhs, R rhs, const std::string& rgn) { + return bout::detail::makePowExpr>( + static_cast(lhs), static_cast::View>(rhs), + lhs.getMesh(), lhs.getLocation(), lhs.getDirections(), + bout::detail::getPowRegionID(lhs.getMesh(), rgn), + lhs.getMesh()->getRegion(rgn), bout::detail::getPerpYIndex(lhs)); +} - BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs, rhs[i]); } +template +std::enable_if_t && is_expr_field3d_v, + BinaryExpr, R, bout::op::Pow>> +pow(L lhs, const R& rhs) { + return bout::detail::makePowExpr, R>( + static_cast::View>(lhs), static_cast(rhs), + rhs.getMesh(), rhs.getLocation(), rhs.getDirections(), rhs.getRegionID(), + rhs.getMesh()->getRegion("RGN_ALL"), bout::detail::getPerpYIndex(rhs)); +} - checkData(result); - return result; +template +std::enable_if_t && is_expr_field3d_v, + BinaryExpr, R, bout::op::Pow>> +pow(L lhs, const R& rhs, const std::string& rgn) { + return bout::detail::makePowExpr, R>( + static_cast::View>(lhs), static_cast(rhs), + rhs.getMesh(), rhs.getLocation(), rhs.getDirections(), + bout::detail::getPowRegionID(rhs.getMesh(), rgn), + rhs.getMesh()->getRegion(rgn), bout::detail::getPerpYIndex(rhs)); } /*! diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index 02e582c20a..63e20c349e 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -911,7 +911,6 @@ inline auto operator-(const Field3D& f) { /// This loops over the entire domain, including guard/boundary cells by /// default (can be changed using the \p rgn argument). /// If CHECK >= 3 then the result will be checked for non-finite numbers -Field3D pow(const Field3D& lhs, const Field2D& rhs, const std::string& rgn = "RGN_ALL"); FieldPerp pow(const Field3D& lhs, const FieldPerp& rhs, const std::string& rgn = "RGN_ALL"); diff --git a/include/bout/fieldperp.hxx b/include/bout/fieldperp.hxx index 83bb3b79b0..834c82be4d 100644 --- a/include/bout/fieldperp.hxx +++ b/include/bout/fieldperp.hxx @@ -409,6 +409,11 @@ FieldPerp operator/(const FieldPerp& lhs, const Field2D& rhs); FieldPerp operator/(const FieldPerp& lhs, BoutReal rhs); FieldPerp operator/(BoutReal lhs, const FieldPerp& rhs); +FieldPerp pow(const FieldPerp& lhs, const FieldPerp& rhs, + const std::string& rgn = "RGN_ALL"); +FieldPerp pow(const FieldPerp& lhs, BoutReal rhs, const std::string& rgn = "RGN_ALL"); +FieldPerp pow(BoutReal lhs, const FieldPerp& rhs, const std::string& rgn = "RGN_ALL"); + /*! * Unary minus. Returns the negative of given field, * iterates over whole domain including guard/boundary cells. diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 306ff3301f..67db42cf0b 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -79,13 +79,13 @@ Field3D D4DY4_Index(const Field3D& f, bool bndry_flux = true); // Forward declarations of flux limiters // If you want to use your own flux limiter, you need to // #include to instantiate the templates. -class Upwind; -class Fromm; -class MinMod; -class MC; -class Superbee; -class VanAlbada; -class WENO3; +struct Upwind; +struct Fromm; +struct MinMod; +struct MC; +struct Superbee; +struct VanAlbada; +struct WENO3; /*! * Communicate fluxes between processors diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 5da438dcdf..95ea1de05a 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -686,22 +686,6 @@ void Field3D::swapData(Field3D& other) { std::swap(data, other.data); } //////////////// NON-MEMBER FUNCTIONS ////////////////// -Field3D pow(const Field3D& lhs, const Field2D& rhs, const std::string& rgn) { - - // Check if the inputs are allocated - checkData(lhs); - checkData(rhs); - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - // Define and allocate the output result - Field3D result{emptyFrom(lhs)}; - - BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs[i], rhs[i]); } - - checkData(result); - return result; -} - FieldPerp pow(const Field3D& lhs, const FieldPerp& rhs, const std::string& rgn) { checkData(lhs); diff --git a/src/field/fieldperp.cxx b/src/field/fieldperp.cxx index b7b2d9d731..bca551c29d 100644 --- a/src/field/fieldperp.cxx +++ b/src/field/fieldperp.cxx @@ -153,6 +153,41 @@ FieldPerp fromFieldAligned(const FieldPerp& f, const std::string& region) { ///////////////////////////////////////////////// // functions +FieldPerp pow(const FieldPerp& lhs, const FieldPerp& rhs, const std::string& rgn) { + checkData(lhs); + checkData(rhs); + ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); + + FieldPerp result{emptyFrom(lhs)}; + + BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs[i], rhs[i]); } + + checkData(result); + return result; +} + +FieldPerp pow(const FieldPerp& lhs, BoutReal rhs, const std::string& rgn) { + checkData(lhs); + + FieldPerp result{emptyFrom(lhs)}; + + BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs[i], rhs); } + + checkData(result); + return result; +} + +FieldPerp pow(BoutReal lhs, const FieldPerp& rhs, const std::string& rgn) { + checkData(rhs); + + FieldPerp result{emptyFrom(rhs)}; + + BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs, rhs[i]); } + + checkData(result); + return result; +} + const FieldPerp sliceXZ(const Field3D& f, int y) { // Source field should be valid checkData(f); diff --git a/tests/unit/field/test_field2d.cxx b/tests/unit/field/test_field2d.cxx index a91acd4a40..1420c79bc2 100644 --- a/tests/unit/field/test_field2d.cxx +++ b/tests/unit/field/test_field2d.cxx @@ -1168,6 +1168,19 @@ TEST_F(Field2DTest, PowField2DField2D) { EXPECT_TRUE(IsFieldEqual(c, 64.0)); } +TEST_F(Field2DTest, PowExpressionUsesPowOp) { + Field2D field; + + field = 2.0; + const auto expr = field + 1.0; + + EXPECT_TRUE((std::is_same_v, + BinaryExpr, + Constant, bout::op::Pow>>)); + EXPECT_TRUE(IsFieldEqual(pow(expr, 2.0), 9.0)); + EXPECT_TRUE(IsFieldEqual(pow(expr, 2.0, "RGN_ALL"), 9.0)); +} + TEST_F(Field2DTest, Sqrt) { Field2D field; diff --git a/tests/unit/field/test_field3d.cxx b/tests/unit/field/test_field3d.cxx index eee92c806e..4ddc7dde5c 100644 --- a/tests/unit/field/test_field3d.cxx +++ b/tests/unit/field/test_field3d.cxx @@ -1947,6 +1947,37 @@ TEST_F(Field3DTest, PowField3DField3D) { EXPECT_TRUE(IsFieldEqual(c, 64.0)); } +TEST_F(Field3DTest, PowExpressionUsesPowOp) { + Field3D field; + + field = 2.0; + const auto expr = field + 1.0; + + EXPECT_TRUE((std::is_same_v, + BinaryExpr, + Constant, bout::op::Pow>>)); + EXPECT_TRUE(IsFieldEqual(pow(expr, 2.0), 9.0)); +} + +TEST_F(Field3DTest, PowMixedField2DField3DExpressionUsesPowOp) { + Field2D lhs; + Field3D rhs; + + lhs = 2.0; + rhs = 3.0; + + const auto lhs_expr = lhs + 1.0; + const auto rhs_expr = rhs + 1.0; + const auto expr = pow(lhs_expr, rhs_expr); + + EXPECT_TRUE( + (std::is_same_v, + BinaryExpr, + std::decay_t, bout::op::Pow>>)); + EXPECT_TRUE(IsFieldEqual(expr, 81.0)); + EXPECT_TRUE(IsFieldEqual(pow(lhs_expr, rhs_expr, "RGN_ALL"), 81.0)); +} + TEST_F(Field3DTest, Sqrt) { Field3D field; @@ -2132,6 +2163,49 @@ TEST_F(Field3DTestFCI, SqrtField3DParallelPreservesParallelSlices) { } } +TEST_F(Field3DTestFCI, PowField3DParallelExprPreservesParallelSlices) { + Field3DParallel field; + + field = 2.0; + field.yup() = 3.0; + field.ydown() = 4.0; + + const auto expr = pow(field + 1.0, 2.0); + + EXPECT_TRUE((std::is_same_v, + BinaryExpr, + Constant, bout::op::Pow>>)); + + Field3DParallel result{expr}; + + EXPECT_TRUE(result.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(result, 9.0)); + EXPECT_TRUE(IsFieldEqual(result.yup(), 16.0)); + EXPECT_TRUE(IsFieldEqual(result.ydown(), 25.0)); +} + +TEST_F(Field3DTestFCI, PowField3DParallelRegionExprPreservesParallelSlices) { + Field3DParallel field; + + field = 2.0; + field.yup() = 3.0; + field.ydown() = 4.0; + + const auto base = field + 1.0; + const auto expr = pow(base, 2.0, "RGN_ALL"); + + EXPECT_TRUE((std::is_same_v, + BinaryExpr, + Constant, bout::op::Pow>>)); + + Field3DParallel result{expr}; + + EXPECT_TRUE(result.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(result, 9.0)); + EXPECT_TRUE(IsFieldEqual(result.yup(), 16.0)); + EXPECT_TRUE(IsFieldEqual(result.ydown(), 25.0)); +} + TEST_F(Field3DTestFCI, SqrtField3DParallelExprPreservesParallelSlices) { Field3DParallel lhs; Field3D rhs; From c808fa7bf87d5f2be63d790058efa4ae1939e05b Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 10 Jul 2026 10:26:40 -0700 Subject: [PATCH 10/14] floor: Return BinaryExpr Makes the `floor` function lazy, so it can be inlined into template expressions. --- include/bout/field.hxx | 99 +++++++++++++++++------------ tests/unit/field/test_field2d.cxx | 13 ++++ tests/unit/field/test_field3d.cxx | 55 ++++++++++++++++ tests/unit/field/test_fieldperp.cxx | 14 ++++ 4 files changed, 142 insertions(+), 39 deletions(-) diff --git a/include/bout/field.hxx b/include/bout/field.hxx index 05bccf250c..67ce45dc4b 100644 --- a/include/bout/field.hxx +++ b/include/bout/field.hxx @@ -33,10 +33,12 @@ class Field; #include #include #include +#include #include "bout/bout_types.hxx" #include "bout/boutcomm.hxx" #include "bout/boutexception.hxx" +#include "bout/build_config.hxx" #include "bout/field_data.hxx" #include "bout/region.hxx" #include "bout/traits.hxx" @@ -887,6 +889,20 @@ struct Square { return ::SQ(value); } }; + +struct Floor { + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal operator()(int idx, const LView& L, + const RView& R) const { + const BoutReal value = L(idx); + const BoutReal floor_value = R(idx); + return value < floor_value ? floor_value : value; + } + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal operator()(BoutReal value, + BoutReal floor_value) const { + return value < floor_value ? floor_value : value; + } +}; }; // namespace bout::op template > @@ -1065,46 +1081,51 @@ class Field3DParallel; /// @param[in] f The floor value /// @param[in] rgn The region to calculate the result over template > -inline T floor(const T& var, BoutReal f, const std::string& rgn = "RGN_ALL") { - checkData(var); - T result = copy(var); +inline auto floor(const T& var, BoutReal f, const std::string& rgn = "RGN_ALL") { + using ResT = bout::detail::UnaryFieldResult_t; + return BinaryExpr, bout::op::Floor>{ + static_cast(var), + static_cast::View>(f), + bout::op::Floor{}, + var.getMesh(), + var.getLocation(), + var.getDirections(), + bout::detail::getUnaryRegionID(var.getMesh(), rgn), + var.getMesh()->template getRegion(rgn), + bout::detail::getPerpYIndex(var)}; +} - BOUT_FOR(d, var.getRegion(rgn)) { - if (result[d] < f) { - result[d] = f; - } - } - if constexpr (std::is_same_v) { - if (var.hasParallelSlices()) { - for (size_t i = 0; i < result.numberParallelSlices(); ++i) { - if (result.yup(i).isAllocated()) { - BOUT_FOR(d, result.yup(i).getRegion(rgn)) { - if (result.yup(i)[d] < f) { - result.yup(i)[d] = f; - } - } - } else { - if (result.isFci()) { - throw BoutException("Expected parallel slice to be allocated"); - } - } - if (result.ydown(i).isAllocated()) { - BOUT_FOR(d, result.ydown(i).getRegion(rgn)) { - if (result.ydown(i)[d] < f) { - result.ydown(i)[d] = f; - } - } - } else { - if (result.isFci()) { - throw BoutException("Expected parallel slice to be allocated"); - } - } - } - } - } else { - result.clearParallelSlices(); - } - return result; +template +inline auto floor(const BinaryExpr& var, BoutReal f) { + using UnaryResT = bout::detail::UnaryFieldResult_t; + return BinaryExpr, Constant, + bout::op::Floor>{ + static_cast::View>(var), + static_cast::View>(f), + bout::op::Floor{}, + var.getMesh(), + var.getLocation(), + var.getDirections(), + var.getRegionID(), + var.indices, + bout::detail::getPerpYIndex(var)}; +} + +template +inline auto floor(const BinaryExpr& var, BoutReal f, + const std::string& rgn) { + using UnaryResT = bout::detail::UnaryFieldResult_t; + return BinaryExpr, Constant, + bout::op::Floor>{ + static_cast::View>(var), + static_cast::View>(f), + bout::op::Floor{}, + var.getMesh(), + var.getLocation(), + var.getDirections(), + bout::detail::getUnaryRegionID(var.getMesh(), rgn), + var.getMesh()->template getRegion(rgn), + bout::detail::getPerpYIndex(var)}; } #undef FIELD_FUNC diff --git a/tests/unit/field/test_field2d.cxx b/tests/unit/field/test_field2d.cxx index 1420c79bc2..91721714ec 100644 --- a/tests/unit/field/test_field2d.cxx +++ b/tests/unit/field/test_field2d.cxx @@ -1330,6 +1330,19 @@ TEST_F(Field2DTest, Floor) { EXPECT_TRUE(IsFieldEqual(floor(field, floor_value), floor_value)); } +TEST_F(Field2DTest, FloorExpressionUsesFloorOp) { + Field2D field; + + field = 2.0; + const auto expr = field + 1.0; + + EXPECT_TRUE((std::is_same_v, + BinaryExpr, + Constant, bout::op::Floor>>)); + EXPECT_TRUE(IsFieldEqual(floor(expr, 5.0), 5.0)); + EXPECT_TRUE(IsFieldEqual(floor(expr, 5.0, "RGN_ALL"), 5.0)); +} + TEST_F(Field2DTest, Min) { Field2D field; diff --git a/tests/unit/field/test_field3d.cxx b/tests/unit/field/test_field3d.cxx index 4ddc7dde5c..6a6ecb9aea 100644 --- a/tests/unit/field/test_field3d.cxx +++ b/tests/unit/field/test_field3d.cxx @@ -2543,6 +2543,61 @@ TEST_F(Field3DTest, Floor) { EXPECT_TRUE(IsFieldEqual(floor(field, floor_value), floor_value)); } +TEST_F(Field3DTest, FloorExpressionUsesFloorOp) { + Field3D field; + + field = 2.0; + const auto expr = field + 1.0; + + EXPECT_TRUE((std::is_same_v, + BinaryExpr, + Constant, bout::op::Floor>>)); + EXPECT_TRUE(IsFieldEqual(floor(expr, 5.0), 5.0)); +} + +TEST_F(Field3DTestFCI, FloorField3DParallelExprPreservesParallelSlices) { + Field3DParallel field; + + field = 2.0; + field.yup() = 3.0; + field.ydown() = 4.0; + + const auto expr = floor(field + 1.0, 4.5); + + EXPECT_TRUE((std::is_same_v, + BinaryExpr, + Constant, bout::op::Floor>>)); + + Field3DParallel result{expr}; + + EXPECT_TRUE(result.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(result, 4.5)); + EXPECT_TRUE(IsFieldEqual(result.yup(), 4.5)); + EXPECT_TRUE(IsFieldEqual(result.ydown(), 5.0)); +} + +TEST_F(Field3DTestFCI, FloorField3DParallelRegionExprPreservesParallelSlices) { + Field3DParallel field; + + field = 2.0; + field.yup() = 3.0; + field.ydown() = 4.0; + + const auto base = field + 1.0; + const auto expr = floor(base, 4.5, "RGN_ALL"); + + EXPECT_TRUE((std::is_same_v, + BinaryExpr, + Constant, bout::op::Floor>>)); + + Field3DParallel result{expr}; + + EXPECT_TRUE(result.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(result, 4.5)); + EXPECT_TRUE(IsFieldEqual(result.yup(), 4.5)); + EXPECT_TRUE(IsFieldEqual(result.ydown(), 5.0)); +} + TEST_F(Field3DTest, Min) { Field3D field; diff --git a/tests/unit/field/test_fieldperp.cxx b/tests/unit/field/test_fieldperp.cxx index 46f07d589f..754316e077 100644 --- a/tests/unit/field/test_fieldperp.cxx +++ b/tests/unit/field/test_fieldperp.cxx @@ -1713,6 +1713,20 @@ TEST_F(FieldPerpTest, Floor) { EXPECT_TRUE(IsFieldEqual(floor(field, floor_value), floor_value)); } +TEST_F(FieldPerpTest, FloorExpressionUsesFloorOp) { + FieldPerp field; + field.setIndex(0); + + field = 2.0; + const auto expr = field + 1.0; + + EXPECT_TRUE((std::is_same_v, + BinaryExpr, + Constant, bout::op::Floor>>)); + EXPECT_TRUE(IsFieldEqual(floor(expr, 5.0), 5.0)); + EXPECT_TRUE(IsFieldEqual(floor(expr, 5.0, "RGN_ALL"), 5.0)); +} + TEST_F(FieldPerpTest, Min) { FieldPerp field; field.setIndex(0); From c7493e49edb7d6300e2cd0d13fac677b3b4e93bc Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 10 Jul 2026 13:07:07 -0700 Subject: [PATCH 11/14] generated_fieldops: Remove eager functions Arithmetic operators now use the lazy `BinaryExpr` expressions. --- include/bout/field3d.hxx | 10 - src/field/gen_fieldops.jinja | 6 +- src/field/gen_fieldops.py | 18 + src/field/generated_fieldops.cxx | 1336 ------------------------------ tests/unit/fake_mesh.hxx | 3 +- 5 files changed, 24 insertions(+), 1349 deletions(-) diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index 63e20c349e..ade329bb6c 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -1211,16 +1211,6 @@ struct is_expr_field3d> : std::integral_constant>::value || is_expr_field3d_v>> {}; -Field3D operator+(const Field2D& lhs, const Field3DParallel& rhs); -Field3D operator-(const Field2D& lhs, const Field3DParallel& rhs); -Field3D operator*(const Field2D& lhs, const Field3DParallel& rhs); -Field3D operator/(const Field2D& lhs, const Field3DParallel& rhs); - -Field3D operator+(const Field3DParallel& lhs, const Field2D& rhs); -Field3D operator-(const Field3DParallel& lhs, const Field2D& rhs); -Field3D operator*(const Field3DParallel& lhs, const Field2D& rhs); -Field3D operator/(const Field3DParallel& lhs, const Field2D& rhs); - inline Field3DParallel filledFrom(const Field3DParallel& f, const std::function& func) { diff --git a/src/field/gen_fieldops.jinja b/src/field/gen_fieldops.jinja index 913acadf7b..89f46469bd 100644 --- a/src/field/gen_fieldops.jinja +++ b/src/field/gen_fieldops.jinja @@ -1,6 +1,7 @@ {% set use_parallel_arg = lhs.field_type == "Field3DParallel" or rhs.field_type == "Field3DParallel" %} {% set use_raja_path = region_loop == "BOUT_FOR_RAJA" and not use_parallel_arg %} +{% if emit_free_binary %} // Provide the C++ wrapper for {{operator_name}} of {{lhs}} and {{rhs}} {{out}} operator{{operator}}(const {{lhs.passByReference}}, const {{rhs.passByReference}}) { {% if lhs != "BoutReal" and rhs != "BoutReal" %} @@ -140,8 +141,9 @@ checkData({{out.name}}); return {{out.name}}; } +{% endif %} -{% if out.field_type == lhs.field_type and lhs == "Field3D" %} +{% if emit_update_inplace %} // Provide the C++ operator to update {{lhs}} by {{operator_name}} with {{rhs}} {{lhs}} &{{lhs}}::update_{{operator_name}}_inplace(const {{rhs.passByReference}}) { // only if data is unique we update the field @@ -234,7 +236,7 @@ {% endif %} -{% if out.field_type == lhs.field_type %} +{% if emit_member_op_equals %} // Provide the C++ operator to update {{lhs}} by {{operator_name}} with {{rhs}} {{lhs}} &{{lhs}}::operator{{operator}}=(const {{rhs.passByReference}}) { // only if data is unique we update the field diff --git a/src/field/gen_fieldops.py b/src/field/gen_fieldops.py index 6610286af5..e77530ff97 100755 --- a/src/field/gen_fieldops.py +++ b/src/field/gen_fieldops.py @@ -254,6 +254,21 @@ def returnType(f1, f2): return copy(field3D) +def emit_free_binary_wrapper(out): + """Return True if this operator still needs an eager non-member wrapper.""" + return out.field_type == "FieldPerp" + + +def emit_update_inplace(lhs, out): + """Return True if this operator needs a Field3D update_*_inplace definition.""" + return out.field_type == lhs.field_type and lhs.field_type == "Field3D" + + +def emit_member_operator_equals(lhs, out): + """Return True if this operator needs an in-place operator definition.""" + return out.field_type == lhs.field_type + + if __name__ == "__main__": parser = argparse.ArgumentParser( description="Generate code for the Field arithmetic operators" @@ -365,6 +380,9 @@ def returnType(f1, f2): "out": out, "lhs": lhs, "rhs": rhs, + "emit_free_binary": emit_free_binary_wrapper(out), + "emit_update_inplace": emit_update_inplace(lhs, out), + "emit_member_op_equals": emit_member_operator_equals(lhs, out), # "region_loop": region_loop, "region_name": region_name, diff --git a/src/field/generated_fieldops.cxx b/src/field/generated_fieldops.cxx index d47b2a7a89..1e9ade4ba4 100644 --- a/src/field/generated_fieldops.cxx +++ b/src/field/generated_fieldops.cxx @@ -9,23 +9,6 @@ #include #include -// Provide the C++ wrapper for multiplication of Field3D and Field3D -Field3D operator*(const Field3D& lhs, const Field3D& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3D result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] * rhs[index]; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3D by multiplication with Field3D Field3D& Field3D::update_multiplication_inplace(const Field3D& rhs) { // only if data is unique we update the field @@ -78,23 +61,6 @@ Field3D& Field3D::operator*=(const Field3D& rhs) { return *this; } -// Provide the C++ wrapper for division of Field3D and Field3D -Field3D operator/(const Field3D& lhs, const Field3D& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3D result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] / rhs[index]; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3D by division with Field3D Field3D& Field3D::update_division_inplace(const Field3D& rhs) { // only if data is unique we update the field @@ -147,23 +113,6 @@ Field3D& Field3D::operator/=(const Field3D& rhs) { return *this; } -// Provide the C++ wrapper for addition of Field3D and Field3D -Field3D operator+(const Field3D& lhs, const Field3D& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3D result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] + rhs[index]; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3D by addition with Field3D Field3D& Field3D::update_addition_inplace(const Field3D& rhs) { // only if data is unique we update the field @@ -216,23 +165,6 @@ Field3D& Field3D::operator+=(const Field3D& rhs) { return *this; } -// Provide the C++ wrapper for subtraction of Field3D and Field3D -Field3D operator-(const Field3D& lhs, const Field3D& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3D result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] - rhs[index]; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3D by subtraction with Field3D Field3D& Field3D::update_subtraction_inplace(const Field3D& rhs) { // only if data is unique we update the field @@ -285,28 +217,6 @@ Field3D& Field3D::operator-=(const Field3D& rhs) { return *this; } -// Provide the C++ wrapper for multiplication of Field3D and Field2D -Field3D operator*(const Field3D& lhs, const Field2D& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3D result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getRegionID()); - - Mesh* localmesh = lhs.getMesh(); - - BOUT_FOR_SERIAL(index, rhs.getRegion("RGN_ALL")) { - const auto base_ind = localmesh->ind2Dto3D(index); - for (int jz = 0; jz < localmesh->LocalNz; ++jz) { - result[base_ind + jz] = lhs[base_ind + jz] * rhs[index]; - } - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3D by multiplication with Field2D Field3D& Field3D::update_multiplication_inplace(const Field2D& rhs) { // only if data is unique we update the field @@ -365,29 +275,6 @@ Field3D& Field3D::operator*=(const Field2D& rhs) { return *this; } -// Provide the C++ wrapper for division of Field3D and Field2D -Field3D operator/(const Field3D& lhs, const Field2D& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3D result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getRegionID()); - - Mesh* localmesh = lhs.getMesh(); - - BOUT_FOR_SERIAL(index, rhs.getRegion("RGN_ALL")) { - const auto base_ind = localmesh->ind2Dto3D(index); - const auto tmp = 1.0 / rhs[index]; - for (int jz = 0; jz < localmesh->LocalNz; ++jz) { - result[base_ind + jz] = lhs[base_ind + jz] * tmp; - } - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3D by division with Field2D Field3D& Field3D::update_division_inplace(const Field2D& rhs) { // only if data is unique we update the field @@ -448,28 +335,6 @@ Field3D& Field3D::operator/=(const Field2D& rhs) { return *this; } -// Provide the C++ wrapper for addition of Field3D and Field2D -Field3D operator+(const Field3D& lhs, const Field2D& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3D result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getRegionID()); - - Mesh* localmesh = lhs.getMesh(); - - BOUT_FOR_SERIAL(index, rhs.getRegion("RGN_ALL")) { - const auto base_ind = localmesh->ind2Dto3D(index); - for (int jz = 0; jz < localmesh->LocalNz; ++jz) { - result[base_ind + jz] = lhs[base_ind + jz] + rhs[index]; - } - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3D by addition with Field2D Field3D& Field3D::update_addition_inplace(const Field2D& rhs) { // only if data is unique we update the field @@ -528,28 +393,6 @@ Field3D& Field3D::operator+=(const Field2D& rhs) { return *this; } -// Provide the C++ wrapper for subtraction of Field3D and Field2D -Field3D operator-(const Field3D& lhs, const Field2D& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3D result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getRegionID()); - - Mesh* localmesh = lhs.getMesh(); - - BOUT_FOR_SERIAL(index, rhs.getRegion("RGN_ALL")) { - const auto base_ind = localmesh->ind2Dto3D(index); - for (int jz = 0; jz < localmesh->LocalNz; ++jz) { - result[base_ind + jz] = lhs[base_ind + jz] - rhs[index]; - } - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3D by subtraction with Field2D Field3D& Field3D::update_subtraction_inplace(const Field2D& rhs) { // only if data is unique we update the field @@ -684,22 +527,6 @@ FieldPerp operator-(const Field3D& lhs, const FieldPerp& rhs) { return result; } -// Provide the C++ wrapper for multiplication of Field3D and BoutReal -Field3D operator*(const Field3D& lhs, const BoutReal rhs) { - - Field3D result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getRegionID()); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] * rhs; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3D by multiplication with BoutReal Field3D& Field3D::update_multiplication_inplace(const BoutReal rhs) { // only if data is unique we update the field @@ -746,23 +573,6 @@ Field3D& Field3D::operator*=(const BoutReal rhs) { return *this; } -// Provide the C++ wrapper for division of Field3D and BoutReal -Field3D operator/(const Field3D& lhs, const BoutReal rhs) { - - Field3D result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getRegionID()); - - const auto tmp = 1.0 / rhs; - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] * tmp; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3D by division with BoutReal Field3D& Field3D::update_division_inplace(const BoutReal rhs) { // only if data is unique we update the field @@ -811,22 +621,6 @@ Field3D& Field3D::operator/=(const BoutReal rhs) { return *this; } -// Provide the C++ wrapper for addition of Field3D and BoutReal -Field3D operator+(const Field3D& lhs, const BoutReal rhs) { - - Field3D result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getRegionID()); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] + rhs; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3D by addition with BoutReal Field3D& Field3D::update_addition_inplace(const BoutReal rhs) { // only if data is unique we update the field @@ -873,22 +667,6 @@ Field3D& Field3D::operator+=(const BoutReal rhs) { return *this; } -// Provide the C++ wrapper for subtraction of Field3D and BoutReal -Field3D operator-(const Field3D& lhs, const BoutReal rhs) { - - Field3D result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getRegionID()); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] - rhs; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3D by subtraction with BoutReal Field3D& Field3D::update_subtraction_inplace(const BoutReal rhs) { // only if data is unique we update the field @@ -935,109 +713,6 @@ Field3D& Field3D::operator-=(const BoutReal rhs) { return *this; } -// Provide the C++ wrapper for multiplication of Field2D and Field3D -Field3D operator*(const Field2D& lhs, const Field3D& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3D result{emptyFrom(rhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(rhs.getRegionID()); - - Mesh* localmesh = lhs.getMesh(); - - BOUT_FOR_SERIAL(index, lhs.getRegion("RGN_ALL")) { - const auto base_ind = localmesh->ind2Dto3D(index); - for (int jz = 0; jz < localmesh->LocalNz; ++jz) { - result[base_ind + jz] = lhs[index] * rhs[base_ind + jz]; - } - } - checkData(result); - return result; -} - -// Provide the C++ wrapper for division of Field2D and Field3D -Field3D operator/(const Field2D& lhs, const Field3D& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3D result{emptyFrom(rhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(rhs.getRegionID()); - - Mesh* localmesh = lhs.getMesh(); - - BOUT_FOR_SERIAL(index, lhs.getRegion("RGN_ALL")) { - const auto base_ind = localmesh->ind2Dto3D(index); - for (int jz = 0; jz < localmesh->LocalNz; ++jz) { - result[base_ind + jz] = lhs[index] / rhs[base_ind + jz]; - } - } - checkData(result); - return result; -} - -// Provide the C++ wrapper for addition of Field2D and Field3D -Field3D operator+(const Field2D& lhs, const Field3D& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3D result{emptyFrom(rhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(rhs.getRegionID()); - - Mesh* localmesh = lhs.getMesh(); - - BOUT_FOR_SERIAL(index, lhs.getRegion("RGN_ALL")) { - const auto base_ind = localmesh->ind2Dto3D(index); - for (int jz = 0; jz < localmesh->LocalNz; ++jz) { - result[base_ind + jz] = lhs[index] + rhs[base_ind + jz]; - } - } - checkData(result); - return result; -} - -// Provide the C++ wrapper for subtraction of Field2D and Field3D -Field3D operator-(const Field2D& lhs, const Field3D& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3D result{emptyFrom(rhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(rhs.getRegionID()); - - Mesh* localmesh = lhs.getMesh(); - - BOUT_FOR_SERIAL(index, lhs.getRegion("RGN_ALL")) { - const auto base_ind = localmesh->ind2Dto3D(index); - for (int jz = 0; jz < localmesh->LocalNz; ++jz) { - result[base_ind + jz] = lhs[index] - rhs[base_ind + jz]; - } - } - checkData(result); - return result; -} - -// Provide the C++ wrapper for multiplication of Field2D and Field2D -Field2D operator*(const Field2D& lhs, const Field2D& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field2D result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] * rhs[index]; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field2D by multiplication with Field2D Field2D& Field2D::operator*=(const Field2D& rhs) { // only if data is unique we update the field @@ -1058,21 +733,6 @@ Field2D& Field2D::operator*=(const Field2D& rhs) { return *this; } -// Provide the C++ wrapper for division of Field2D and Field2D -Field2D operator/(const Field2D& lhs, const Field2D& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field2D result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] / rhs[index]; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field2D by division with Field2D Field2D& Field2D::operator/=(const Field2D& rhs) { // only if data is unique we update the field @@ -1093,21 +753,6 @@ Field2D& Field2D::operator/=(const Field2D& rhs) { return *this; } -// Provide the C++ wrapper for addition of Field2D and Field2D -Field2D operator+(const Field2D& lhs, const Field2D& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field2D result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] + rhs[index]; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field2D by addition with Field2D Field2D& Field2D::operator+=(const Field2D& rhs) { // only if data is unique we update the field @@ -1128,21 +773,6 @@ Field2D& Field2D::operator+=(const Field2D& rhs) { return *this; } -// Provide the C++ wrapper for subtraction of Field2D and Field2D -Field2D operator-(const Field2D& lhs, const Field2D& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field2D result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] - rhs[index]; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field2D by subtraction with Field2D Field2D& Field2D::operator-=(const Field2D& rhs) { // only if data is unique we update the field @@ -1239,20 +869,6 @@ FieldPerp operator-(const Field2D& lhs, const FieldPerp& rhs) { return result; } -// Provide the C++ wrapper for multiplication of Field2D and BoutReal -Field2D operator*(const Field2D& lhs, const BoutReal rhs) { - - Field2D result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] * rhs; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field2D by multiplication with BoutReal Field2D& Field2D::operator*=(const BoutReal rhs) { // only if data is unique we update the field @@ -1272,21 +888,6 @@ Field2D& Field2D::operator*=(const BoutReal rhs) { return *this; } -// Provide the C++ wrapper for division of Field2D and BoutReal -Field2D operator/(const Field2D& lhs, const BoutReal rhs) { - - Field2D result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - const auto tmp = 1.0 / rhs; - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] * tmp; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field2D by division with BoutReal Field2D& Field2D::operator/=(const BoutReal rhs) { // only if data is unique we update the field @@ -1307,20 +908,6 @@ Field2D& Field2D::operator/=(const BoutReal rhs) { return *this; } -// Provide the C++ wrapper for addition of Field2D and BoutReal -Field2D operator+(const Field2D& lhs, const BoutReal rhs) { - - Field2D result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] + rhs; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field2D by addition with BoutReal Field2D& Field2D::operator+=(const BoutReal rhs) { // only if data is unique we update the field @@ -1340,20 +927,6 @@ Field2D& Field2D::operator+=(const BoutReal rhs) { return *this; } -// Provide the C++ wrapper for subtraction of Field2D and BoutReal -Field2D operator-(const Field2D& lhs, const BoutReal rhs) { - - Field2D result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] - rhs; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field2D by subtraction with BoutReal Field2D& Field2D::operator-=(const BoutReal rhs) { // only if data is unique we update the field @@ -2006,126 +1579,6 @@ FieldPerp& FieldPerp::operator-=(const BoutReal rhs) { return *this; } -// Provide the C++ wrapper for multiplication of BoutReal and Field3D -Field3D operator*(const BoutReal lhs, const Field3D& rhs) { - - Field3D result{emptyFrom(rhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(rhs.getRegionID()); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs * rhs[index]; - } - checkData(result); - return result; -} - -// Provide the C++ wrapper for division of BoutReal and Field3D -Field3D operator/(const BoutReal lhs, const Field3D& rhs) { - - Field3D result{emptyFrom(rhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(rhs.getRegionID()); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs / rhs[index]; - } - checkData(result); - return result; -} - -// Provide the C++ wrapper for addition of BoutReal and Field3D -Field3D operator+(const BoutReal lhs, const Field3D& rhs) { - - Field3D result{emptyFrom(rhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(rhs.getRegionID()); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs + rhs[index]; - } - checkData(result); - return result; -} - -// Provide the C++ wrapper for subtraction of BoutReal and Field3D -Field3D operator-(const BoutReal lhs, const Field3D& rhs) { - - Field3D result{emptyFrom(rhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(rhs.getRegionID()); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs - rhs[index]; - } - checkData(result); - return result; -} - -// Provide the C++ wrapper for multiplication of BoutReal and Field2D -Field2D operator*(const BoutReal lhs, const Field2D& rhs) { - - Field2D result{emptyFrom(rhs)}; - checkData(lhs); - checkData(rhs); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs * rhs[index]; - } - checkData(result); - return result; -} - -// Provide the C++ wrapper for division of BoutReal and Field2D -Field2D operator/(const BoutReal lhs, const Field2D& rhs) { - - Field2D result{emptyFrom(rhs)}; - checkData(lhs); - checkData(rhs); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs / rhs[index]; - } - checkData(result); - return result; -} - -// Provide the C++ wrapper for addition of BoutReal and Field2D -Field2D operator+(const BoutReal lhs, const Field2D& rhs) { - - Field2D result{emptyFrom(rhs)}; - checkData(lhs); - checkData(rhs); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs + rhs[index]; - } - checkData(result); - return result; -} - -// Provide the C++ wrapper for subtraction of BoutReal and Field2D -Field2D operator-(const BoutReal lhs, const Field2D& rhs) { - - Field2D result{emptyFrom(rhs)}; - checkData(lhs); - checkData(rhs); - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs - rhs[index]; - } - checkData(result); - return result; -} - // Provide the C++ wrapper for multiplication of BoutReal and FieldPerp FieldPerp operator*(const BoutReal lhs, const FieldPerp& rhs) { @@ -2182,216 +1635,6 @@ FieldPerp operator-(const BoutReal lhs, const FieldPerp& rhs) { return result; } -// Provide the C++ wrapper for multiplication of Field3D and Field3DParallel -Field3DParallel operator*(const Field3D& lhs, const Field3DParallel& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3DParallel result{emptyFrom(rhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); - if (result.isFci()) { - - ASSERT2(lhs.hasParallelSlices()); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - ASSERT2(lhs.ydown(i).isAllocated()); - ASSERT2(lhs.yup(i).isAllocated()); - } - - ASSERT2(rhs.hasParallelSlices()); - for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { - ASSERT2(rhs.ydown(i).isAllocated()); - ASSERT2(rhs.yup(i).isAllocated()); - } - result.splitParallelSlices(); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - result.yup(i) = lhs.yup(i) * rhs.yup(i); - result.ydown(i) = lhs.ydown(i) * rhs.ydown(i); - } - - ASSERT2(result.hasParallelSlices()); - for (size_t i{0}; i < result.numberParallelSlices(); ++i) { - ASSERT2(result.ydown(i).isAllocated()); - ASSERT2(result.yup(i).isAllocated()); - } - } - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] * rhs[index]; - } - checkData(result); - return result; -} - -// Provide the C++ wrapper for division of Field3D and Field3DParallel -Field3DParallel operator/(const Field3D& lhs, const Field3DParallel& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3DParallel result{emptyFrom(rhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); - if (result.isFci()) { - - ASSERT2(lhs.hasParallelSlices()); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - ASSERT2(lhs.ydown(i).isAllocated()); - ASSERT2(lhs.yup(i).isAllocated()); - } - - ASSERT2(rhs.hasParallelSlices()); - for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { - ASSERT2(rhs.ydown(i).isAllocated()); - ASSERT2(rhs.yup(i).isAllocated()); - } - result.splitParallelSlices(); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - result.yup(i) = lhs.yup(i) / rhs.yup(i); - result.ydown(i) = lhs.ydown(i) / rhs.ydown(i); - } - - ASSERT2(result.hasParallelSlices()); - for (size_t i{0}; i < result.numberParallelSlices(); ++i) { - ASSERT2(result.ydown(i).isAllocated()); - ASSERT2(result.yup(i).isAllocated()); - } - } - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] / rhs[index]; - } - checkData(result); - return result; -} - -// Provide the C++ wrapper for addition of Field3D and Field3DParallel -Field3DParallel operator+(const Field3D& lhs, const Field3DParallel& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3DParallel result{emptyFrom(rhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); - if (result.isFci()) { - - ASSERT2(lhs.hasParallelSlices()); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - ASSERT2(lhs.ydown(i).isAllocated()); - ASSERT2(lhs.yup(i).isAllocated()); - } - - ASSERT2(rhs.hasParallelSlices()); - for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { - ASSERT2(rhs.ydown(i).isAllocated()); - ASSERT2(rhs.yup(i).isAllocated()); - } - result.splitParallelSlices(); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - result.yup(i) = lhs.yup(i) + rhs.yup(i); - result.ydown(i) = lhs.ydown(i) + rhs.ydown(i); - } - - ASSERT2(result.hasParallelSlices()); - for (size_t i{0}; i < result.numberParallelSlices(); ++i) { - ASSERT2(result.ydown(i).isAllocated()); - ASSERT2(result.yup(i).isAllocated()); - } - } - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] + rhs[index]; - } - checkData(result); - return result; -} - -// Provide the C++ wrapper for subtraction of Field3D and Field3DParallel -Field3DParallel operator-(const Field3D& lhs, const Field3DParallel& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3DParallel result{emptyFrom(rhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); - if (result.isFci()) { - - ASSERT2(lhs.hasParallelSlices()); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - ASSERT2(lhs.ydown(i).isAllocated()); - ASSERT2(lhs.yup(i).isAllocated()); - } - - ASSERT2(rhs.hasParallelSlices()); - for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { - ASSERT2(rhs.ydown(i).isAllocated()); - ASSERT2(rhs.yup(i).isAllocated()); - } - result.splitParallelSlices(); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - result.yup(i) = lhs.yup(i) - rhs.yup(i); - result.ydown(i) = lhs.ydown(i) - rhs.ydown(i); - } - - ASSERT2(result.hasParallelSlices()); - for (size_t i{0}; i < result.numberParallelSlices(); ++i) { - ASSERT2(result.ydown(i).isAllocated()); - ASSERT2(result.yup(i).isAllocated()); - } - } - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] - rhs[index]; - } - checkData(result); - return result; -} - -// Provide the C++ wrapper for multiplication of Field3DParallel and Field3D -Field3DParallel operator*(const Field3DParallel& lhs, const Field3D& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3DParallel result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); - if (result.isFci()) { - - ASSERT2(lhs.hasParallelSlices()); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - ASSERT2(lhs.ydown(i).isAllocated()); - ASSERT2(lhs.yup(i).isAllocated()); - } - - ASSERT2(rhs.hasParallelSlices()); - for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { - ASSERT2(rhs.ydown(i).isAllocated()); - ASSERT2(rhs.yup(i).isAllocated()); - } - result.splitParallelSlices(); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - result.yup(i) = lhs.yup(i) * rhs.yup(i); - result.ydown(i) = lhs.ydown(i) * rhs.ydown(i); - } - - ASSERT2(result.hasParallelSlices()); - for (size_t i{0}; i < result.numberParallelSlices(); ++i) { - ASSERT2(result.ydown(i).isAllocated()); - ASSERT2(result.yup(i).isAllocated()); - } - } - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] * rhs[index]; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3DParallel by multiplication with Field3D Field3DParallel& Field3DParallel::operator*=(const Field3D& rhs) { // only if data is unique we update the field @@ -2431,48 +1674,6 @@ Field3DParallel& Field3DParallel::operator*=(const Field3D& rhs) { return *this; } -// Provide the C++ wrapper for division of Field3DParallel and Field3D -Field3DParallel operator/(const Field3DParallel& lhs, const Field3D& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3DParallel result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); - if (result.isFci()) { - - ASSERT2(lhs.hasParallelSlices()); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - ASSERT2(lhs.ydown(i).isAllocated()); - ASSERT2(lhs.yup(i).isAllocated()); - } - - ASSERT2(rhs.hasParallelSlices()); - for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { - ASSERT2(rhs.ydown(i).isAllocated()); - ASSERT2(rhs.yup(i).isAllocated()); - } - result.splitParallelSlices(); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - result.yup(i) = lhs.yup(i) / rhs.yup(i); - result.ydown(i) = lhs.ydown(i) / rhs.ydown(i); - } - - ASSERT2(result.hasParallelSlices()); - for (size_t i{0}; i < result.numberParallelSlices(); ++i) { - ASSERT2(result.ydown(i).isAllocated()); - ASSERT2(result.yup(i).isAllocated()); - } - } - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] / rhs[index]; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3DParallel by division with Field3D Field3DParallel& Field3DParallel::operator/=(const Field3D& rhs) { // only if data is unique we update the field @@ -2512,48 +1713,6 @@ Field3DParallel& Field3DParallel::operator/=(const Field3D& rhs) { return *this; } -// Provide the C++ wrapper for addition of Field3DParallel and Field3D -Field3DParallel operator+(const Field3DParallel& lhs, const Field3D& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3DParallel result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); - if (result.isFci()) { - - ASSERT2(lhs.hasParallelSlices()); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - ASSERT2(lhs.ydown(i).isAllocated()); - ASSERT2(lhs.yup(i).isAllocated()); - } - - ASSERT2(rhs.hasParallelSlices()); - for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { - ASSERT2(rhs.ydown(i).isAllocated()); - ASSERT2(rhs.yup(i).isAllocated()); - } - result.splitParallelSlices(); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - result.yup(i) = lhs.yup(i) + rhs.yup(i); - result.ydown(i) = lhs.ydown(i) + rhs.ydown(i); - } - - ASSERT2(result.hasParallelSlices()); - for (size_t i{0}; i < result.numberParallelSlices(); ++i) { - ASSERT2(result.ydown(i).isAllocated()); - ASSERT2(result.yup(i).isAllocated()); - } - } - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] + rhs[index]; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3DParallel by addition with Field3D Field3DParallel& Field3DParallel::operator+=(const Field3D& rhs) { // only if data is unique we update the field @@ -2593,48 +1752,6 @@ Field3DParallel& Field3DParallel::operator+=(const Field3D& rhs) { return *this; } -// Provide the C++ wrapper for subtraction of Field3DParallel and Field3D -Field3DParallel operator-(const Field3DParallel& lhs, const Field3D& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3DParallel result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); - if (result.isFci()) { - - ASSERT2(lhs.hasParallelSlices()); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - ASSERT2(lhs.ydown(i).isAllocated()); - ASSERT2(lhs.yup(i).isAllocated()); - } - - ASSERT2(rhs.hasParallelSlices()); - for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { - ASSERT2(rhs.ydown(i).isAllocated()); - ASSERT2(rhs.yup(i).isAllocated()); - } - result.splitParallelSlices(); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - result.yup(i) = lhs.yup(i) - rhs.yup(i); - result.ydown(i) = lhs.ydown(i) - rhs.ydown(i); - } - - ASSERT2(result.hasParallelSlices()); - for (size_t i{0}; i < result.numberParallelSlices(); ++i) { - ASSERT2(result.ydown(i).isAllocated()); - ASSERT2(result.yup(i).isAllocated()); - } - } - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] - rhs[index]; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3DParallel by subtraction with Field3D Field3DParallel& Field3DParallel::operator-=(const Field3D& rhs) { // only if data is unique we update the field @@ -2674,48 +1791,6 @@ Field3DParallel& Field3DParallel::operator-=(const Field3D& rhs) { return *this; } -// Provide the C++ wrapper for multiplication of Field3DParallel and Field3DParallel -Field3DParallel operator*(const Field3DParallel& lhs, const Field3DParallel& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3DParallel result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); - if (result.isFci()) { - - ASSERT2(lhs.hasParallelSlices()); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - ASSERT2(lhs.ydown(i).isAllocated()); - ASSERT2(lhs.yup(i).isAllocated()); - } - - ASSERT2(rhs.hasParallelSlices()); - for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { - ASSERT2(rhs.ydown(i).isAllocated()); - ASSERT2(rhs.yup(i).isAllocated()); - } - result.splitParallelSlices(); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - result.yup(i) = lhs.yup(i) * rhs.yup(i); - result.ydown(i) = lhs.ydown(i) * rhs.ydown(i); - } - - ASSERT2(result.hasParallelSlices()); - for (size_t i{0}; i < result.numberParallelSlices(); ++i) { - ASSERT2(result.ydown(i).isAllocated()); - ASSERT2(result.yup(i).isAllocated()); - } - } - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] * rhs[index]; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3DParallel by multiplication with Field3DParallel Field3DParallel& Field3DParallel::operator*=(const Field3DParallel& rhs) { // only if data is unique we update the field @@ -2755,48 +1830,6 @@ Field3DParallel& Field3DParallel::operator*=(const Field3DParallel& rhs) { return *this; } -// Provide the C++ wrapper for division of Field3DParallel and Field3DParallel -Field3DParallel operator/(const Field3DParallel& lhs, const Field3DParallel& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3DParallel result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); - if (result.isFci()) { - - ASSERT2(lhs.hasParallelSlices()); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - ASSERT2(lhs.ydown(i).isAllocated()); - ASSERT2(lhs.yup(i).isAllocated()); - } - - ASSERT2(rhs.hasParallelSlices()); - for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { - ASSERT2(rhs.ydown(i).isAllocated()); - ASSERT2(rhs.yup(i).isAllocated()); - } - result.splitParallelSlices(); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - result.yup(i) = lhs.yup(i) / rhs.yup(i); - result.ydown(i) = lhs.ydown(i) / rhs.ydown(i); - } - - ASSERT2(result.hasParallelSlices()); - for (size_t i{0}; i < result.numberParallelSlices(); ++i) { - ASSERT2(result.ydown(i).isAllocated()); - ASSERT2(result.yup(i).isAllocated()); - } - } - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] / rhs[index]; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3DParallel by division with Field3DParallel Field3DParallel& Field3DParallel::operator/=(const Field3DParallel& rhs) { // only if data is unique we update the field @@ -2836,48 +1869,6 @@ Field3DParallel& Field3DParallel::operator/=(const Field3DParallel& rhs) { return *this; } -// Provide the C++ wrapper for addition of Field3DParallel and Field3DParallel -Field3DParallel operator+(const Field3DParallel& lhs, const Field3DParallel& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3DParallel result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); - if (result.isFci()) { - - ASSERT2(lhs.hasParallelSlices()); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - ASSERT2(lhs.ydown(i).isAllocated()); - ASSERT2(lhs.yup(i).isAllocated()); - } - - ASSERT2(rhs.hasParallelSlices()); - for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { - ASSERT2(rhs.ydown(i).isAllocated()); - ASSERT2(rhs.yup(i).isAllocated()); - } - result.splitParallelSlices(); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - result.yup(i) = lhs.yup(i) + rhs.yup(i); - result.ydown(i) = lhs.ydown(i) + rhs.ydown(i); - } - - ASSERT2(result.hasParallelSlices()); - for (size_t i{0}; i < result.numberParallelSlices(); ++i) { - ASSERT2(result.ydown(i).isAllocated()); - ASSERT2(result.yup(i).isAllocated()); - } - } - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] + rhs[index]; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3DParallel by addition with Field3DParallel Field3DParallel& Field3DParallel::operator+=(const Field3DParallel& rhs) { // only if data is unique we update the field @@ -2917,48 +1908,6 @@ Field3DParallel& Field3DParallel::operator+=(const Field3DParallel& rhs) { return *this; } -// Provide the C++ wrapper for subtraction of Field3DParallel and Field3DParallel -Field3DParallel operator-(const Field3DParallel& lhs, const Field3DParallel& rhs) { - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - Field3DParallel result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getMesh()->getCommonRegion(lhs.getRegionID(), rhs.getRegionID())); - if (result.isFci()) { - - ASSERT2(lhs.hasParallelSlices()); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - ASSERT2(lhs.ydown(i).isAllocated()); - ASSERT2(lhs.yup(i).isAllocated()); - } - - ASSERT2(rhs.hasParallelSlices()); - for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { - ASSERT2(rhs.ydown(i).isAllocated()); - ASSERT2(rhs.yup(i).isAllocated()); - } - result.splitParallelSlices(); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - result.yup(i) = lhs.yup(i) - rhs.yup(i); - result.ydown(i) = lhs.ydown(i) - rhs.ydown(i); - } - - ASSERT2(result.hasParallelSlices()); - for (size_t i{0}; i < result.numberParallelSlices(); ++i) { - ASSERT2(result.ydown(i).isAllocated()); - ASSERT2(result.yup(i).isAllocated()); - } - } - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] - rhs[index]; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3DParallel by subtraction with Field3DParallel Field3DParallel& Field3DParallel::operator-=(const Field3DParallel& rhs) { // only if data is unique we update the field @@ -2998,42 +1947,6 @@ Field3DParallel& Field3DParallel::operator-=(const Field3DParallel& rhs) { return *this; } -// Provide the C++ wrapper for multiplication of Field3DParallel and BoutReal -Field3DParallel operator*(const Field3DParallel& lhs, const BoutReal rhs) { - - Field3DParallel result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getRegionID()); - if (result.isFci()) { - - ASSERT2(lhs.hasParallelSlices()); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - ASSERT2(lhs.ydown(i).isAllocated()); - ASSERT2(lhs.yup(i).isAllocated()); - } - - result.splitParallelSlices(); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - result.yup(i) = lhs.yup(i) * rhs; - result.ydown(i) = lhs.ydown(i) * rhs; - } - - ASSERT2(result.hasParallelSlices()); - for (size_t i{0}; i < result.numberParallelSlices(); ++i) { - ASSERT2(result.ydown(i).isAllocated()); - ASSERT2(result.yup(i).isAllocated()); - } - } - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] * rhs; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3DParallel by multiplication with BoutReal Field3DParallel& Field3DParallel::operator*=(const BoutReal rhs) { // only if data is unique we update the field @@ -3070,43 +1983,6 @@ Field3DParallel& Field3DParallel::operator*=(const BoutReal rhs) { return *this; } -// Provide the C++ wrapper for division of Field3DParallel and BoutReal -Field3DParallel operator/(const Field3DParallel& lhs, const BoutReal rhs) { - - Field3DParallel result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getRegionID()); - if (result.isFci()) { - - ASSERT2(lhs.hasParallelSlices()); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - ASSERT2(lhs.ydown(i).isAllocated()); - ASSERT2(lhs.yup(i).isAllocated()); - } - - result.splitParallelSlices(); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - result.yup(i) = lhs.yup(i) / rhs; - result.ydown(i) = lhs.ydown(i) / rhs; - } - - ASSERT2(result.hasParallelSlices()); - for (size_t i{0}; i < result.numberParallelSlices(); ++i) { - ASSERT2(result.ydown(i).isAllocated()); - ASSERT2(result.yup(i).isAllocated()); - } - } - - const auto tmp = 1.0 / rhs; - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] * tmp; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3DParallel by division with BoutReal Field3DParallel& Field3DParallel::operator/=(const BoutReal rhs) { // only if data is unique we update the field @@ -3143,42 +2019,6 @@ Field3DParallel& Field3DParallel::operator/=(const BoutReal rhs) { return *this; } -// Provide the C++ wrapper for addition of Field3DParallel and BoutReal -Field3DParallel operator+(const Field3DParallel& lhs, const BoutReal rhs) { - - Field3DParallel result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getRegionID()); - if (result.isFci()) { - - ASSERT2(lhs.hasParallelSlices()); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - ASSERT2(lhs.ydown(i).isAllocated()); - ASSERT2(lhs.yup(i).isAllocated()); - } - - result.splitParallelSlices(); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - result.yup(i) = lhs.yup(i) + rhs; - result.ydown(i) = lhs.ydown(i) + rhs; - } - - ASSERT2(result.hasParallelSlices()); - for (size_t i{0}; i < result.numberParallelSlices(); ++i) { - ASSERT2(result.ydown(i).isAllocated()); - ASSERT2(result.yup(i).isAllocated()); - } - } - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] + rhs; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3DParallel by addition with BoutReal Field3DParallel& Field3DParallel::operator+=(const BoutReal rhs) { // only if data is unique we update the field @@ -3215,42 +2055,6 @@ Field3DParallel& Field3DParallel::operator+=(const BoutReal rhs) { return *this; } -// Provide the C++ wrapper for subtraction of Field3DParallel and BoutReal -Field3DParallel operator-(const Field3DParallel& lhs, const BoutReal rhs) { - - Field3DParallel result{emptyFrom(lhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(lhs.getRegionID()); - if (result.isFci()) { - - ASSERT2(lhs.hasParallelSlices()); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - ASSERT2(lhs.ydown(i).isAllocated()); - ASSERT2(lhs.yup(i).isAllocated()); - } - - result.splitParallelSlices(); - for (size_t i{0}; i < lhs.numberParallelSlices(); ++i) { - result.yup(i) = lhs.yup(i) - rhs; - result.ydown(i) = lhs.ydown(i) - rhs; - } - - ASSERT2(result.hasParallelSlices()); - for (size_t i{0}; i < result.numberParallelSlices(); ++i) { - ASSERT2(result.ydown(i).isAllocated()); - ASSERT2(result.yup(i).isAllocated()); - } - } - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs[index] - rhs; - } - checkData(result); - return result; -} - // Provide the C++ operator to update Field3DParallel by subtraction with BoutReal Field3DParallel& Field3DParallel::operator-=(const BoutReal rhs) { // only if data is unique we update the field @@ -3286,143 +2090,3 @@ Field3DParallel& Field3DParallel::operator-=(const BoutReal rhs) { } return *this; } - -// Provide the C++ wrapper for multiplication of BoutReal and Field3DParallel -Field3DParallel operator*(const BoutReal lhs, const Field3DParallel& rhs) { - - Field3DParallel result{emptyFrom(rhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(rhs.getRegionID()); - if (result.isFci()) { - - ASSERT2(rhs.hasParallelSlices()); - for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { - ASSERT2(rhs.ydown(i).isAllocated()); - ASSERT2(rhs.yup(i).isAllocated()); - } - result.splitParallelSlices(); - for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { - result.yup(i) = lhs * rhs.yup(i); - result.ydown(i) = lhs * rhs.ydown(i); - } - - ASSERT2(result.hasParallelSlices()); - for (size_t i{0}; i < result.numberParallelSlices(); ++i) { - ASSERT2(result.ydown(i).isAllocated()); - ASSERT2(result.yup(i).isAllocated()); - } - } - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs * rhs[index]; - } - checkData(result); - return result; -} - -// Provide the C++ wrapper for division of BoutReal and Field3DParallel -Field3DParallel operator/(const BoutReal lhs, const Field3DParallel& rhs) { - - Field3DParallel result{emptyFrom(rhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(rhs.getRegionID()); - if (result.isFci()) { - - ASSERT2(rhs.hasParallelSlices()); - for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { - ASSERT2(rhs.ydown(i).isAllocated()); - ASSERT2(rhs.yup(i).isAllocated()); - } - result.splitParallelSlices(); - for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { - result.yup(i) = lhs / rhs.yup(i); - result.ydown(i) = lhs / rhs.ydown(i); - } - - ASSERT2(result.hasParallelSlices()); - for (size_t i{0}; i < result.numberParallelSlices(); ++i) { - ASSERT2(result.ydown(i).isAllocated()); - ASSERT2(result.yup(i).isAllocated()); - } - } - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs / rhs[index]; - } - checkData(result); - return result; -} - -// Provide the C++ wrapper for addition of BoutReal and Field3DParallel -Field3DParallel operator+(const BoutReal lhs, const Field3DParallel& rhs) { - - Field3DParallel result{emptyFrom(rhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(rhs.getRegionID()); - if (result.isFci()) { - - ASSERT2(rhs.hasParallelSlices()); - for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { - ASSERT2(rhs.ydown(i).isAllocated()); - ASSERT2(rhs.yup(i).isAllocated()); - } - result.splitParallelSlices(); - for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { - result.yup(i) = lhs + rhs.yup(i); - result.ydown(i) = lhs + rhs.ydown(i); - } - - ASSERT2(result.hasParallelSlices()); - for (size_t i{0}; i < result.numberParallelSlices(); ++i) { - ASSERT2(result.ydown(i).isAllocated()); - ASSERT2(result.yup(i).isAllocated()); - } - } - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs + rhs[index]; - } - checkData(result); - return result; -} - -// Provide the C++ wrapper for subtraction of BoutReal and Field3DParallel -Field3DParallel operator-(const BoutReal lhs, const Field3DParallel& rhs) { - - Field3DParallel result{emptyFrom(rhs)}; - checkData(lhs); - checkData(rhs); - - result.setRegion(rhs.getRegionID()); - if (result.isFci()) { - - ASSERT2(rhs.hasParallelSlices()); - for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { - ASSERT2(rhs.ydown(i).isAllocated()); - ASSERT2(rhs.yup(i).isAllocated()); - } - result.splitParallelSlices(); - for (size_t i{0}; i < rhs.numberParallelSlices(); ++i) { - result.yup(i) = lhs - rhs.yup(i); - result.ydown(i) = lhs - rhs.ydown(i); - } - - ASSERT2(result.hasParallelSlices()); - for (size_t i{0}; i < result.numberParallelSlices(); ++i) { - ASSERT2(result.ydown(i).isAllocated()); - ASSERT2(result.yup(i).isAllocated()); - } - } - - BOUT_FOR_SERIAL(index, result.getValidRegionWithDefault("RGN_ALL")) { - result[index] = lhs - rhs[index]; - } - checkData(result); - return result; -} diff --git a/tests/unit/fake_mesh.hxx b/tests/unit/fake_mesh.hxx index 6559e554f2..216d04e18a 100644 --- a/tests/unit/fake_mesh.hxx +++ b/tests/unit/fake_mesh.hxx @@ -13,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -90,7 +91,7 @@ public: void setCoordinates(std::shared_ptr coords, CELL_LOC location = CELL_CENTRE) { - coords_map[location] = coords; + coords_map[location] = std::move(coords); } void setGridDataSource(GridDataSource* source_in) { source = source_in; } From d41d93c85d37993a5109815c66d714728a9c4f95 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 10 Jul 2026 16:30:57 -0700 Subject: [PATCH 12/14] Field3DParallel::View simplify Remove unnecessary data members that dublicate the Field3D::View base. --- include/bout/field3d.hxx | 46 ++++++++++++++++------------------------ 1 file changed, 18 insertions(+), 28 deletions(-) diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index ade329bb6c..5562c58b49 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -488,11 +488,13 @@ public: BOUT_HOST_DEVICE BOUT_FORCEINLINE int numberParallelSlices() const { return num_parallel_slices; } + /// Not a DEVICE function because it dereferences a Field3D pointer BOUT_FORCEINLINE View yup(int slice = 0) const { ASSERT2(slice < num_parallel_slices); ASSERT2(yup_fields[slice].isAllocated()); return static_cast(yup_fields[slice]); } + /// Not a DEVICE function because it dereferences a Field3D pointer BOUT_FORCEINLINE View ydown(int slice = 0) const { ASSERT2(slice < num_parallel_slices); ASSERT2(ydown_fields[slice].isAllocated()); @@ -500,14 +502,11 @@ public: } }; operator View() { - return View{&data[0], hasParallelSlices() ? yup_fields.data() : nullptr, - hasParallelSlices() ? ydown_fields.data() : nullptr, + return View{&data[0], yup_fields.data(), ydown_fields.data(), static_cast(numberParallelSlices())}; } operator View() const { - return View{const_cast(&data[0]), - hasParallelSlices() ? yup_fields.data() : nullptr, - hasParallelSlices() ? ydown_fields.data() : nullptr, + return View{const_cast(&data[0]), yup_fields.data(), ydown_fields.data(), static_cast(numberParallelSlices())}; } //operator View() const { return View{&data[0]}; } @@ -1064,9 +1063,6 @@ public: struct View { Field3D::View base; - const Field3D* yup_fields{nullptr}; - const Field3D* ydown_fields{nullptr}; - int num_parallel_slices{0}; BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal operator()(int idx) const { return base(idx); @@ -1082,32 +1078,26 @@ public: return *this; } - BOUT_FORCEINLINE bool hasParallelSlices() const { return num_parallel_slices > 0; } - BOUT_FORCEINLINE int numberParallelSlices() const { return num_parallel_slices; } + BOUT_FORCEINLINE bool hasParallelSlices() const { return base.hasParallelSlices(); } + BOUT_FORCEINLINE int numberParallelSlices() const { + return base.numberParallelSlices(); + } + /// Not a DEVICE function because it dereferences a Field3D pointer BOUT_FORCEINLINE View yup(int slice = 0) const { - ASSERT2(slice < num_parallel_slices); - ASSERT2(yup_fields[slice].isAllocated()); - return View{static_cast(yup_fields[slice]), nullptr, nullptr, 0}; + ASSERT2(slice < base.num_parallel_slices); + ASSERT2(base.yup_fields[slice].isAllocated()); + return View{static_cast(base.yup_fields[slice])}; } + /// Not a DEVICE function because it dereferences a Field3D pointer BOUT_FORCEINLINE View ydown(int slice = 0) const { - ASSERT2(slice < num_parallel_slices); - ASSERT2(ydown_fields[slice].isAllocated()); - return View{static_cast(ydown_fields[slice]), nullptr, nullptr, 0}; + ASSERT2(slice < base.num_parallel_slices); + ASSERT2(base.ydown_fields[slice].isAllocated()); + return View{static_cast(base.ydown_fields[slice])}; } }; - operator View() { - return View{static_cast(*this), - hasParallelSlices() ? yup_fields.data() : nullptr, - hasParallelSlices() ? ydown_fields.data() : nullptr, - static_cast(numberParallelSlices())}; - } - operator View() const { - return View{static_cast(*this), - hasParallelSlices() ? yup_fields.data() : nullptr, - hasParallelSlices() ? ydown_fields.data() : nullptr, - static_cast(numberParallelSlices())}; - } + operator View() { return View{static_cast(*this)}; } + operator View() const { return View{static_cast(*this)}; } Field3DParallel& operator*=(const Field3D&); Field3DParallel& operator/=(const Field3D&); From 3c6565df8e66cfeb7c70ea1ccaad90001c7da664 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 10 Jul 2026 16:49:40 -0700 Subject: [PATCH 13/14] manual: Updates for Field3DParallel expressions --- manual/sphinx/developer_docs/data_types.rst | 33 ++++++++++------- .../sphinx/user_docs/algebraic_operators.rst | 17 +++++++-- manual/sphinx/user_docs/field_expressions.rst | 35 +++++++++++++++---- 3 files changed, 64 insertions(+), 21 deletions(-) diff --git a/manual/sphinx/developer_docs/data_types.rst b/manual/sphinx/developer_docs/data_types.rst index 7feb3945aa..ba4faa6e42 100644 --- a/manual/sphinx/developer_docs/data_types.rst +++ b/manual/sphinx/developer_docs/data_types.rst @@ -515,37 +515,44 @@ central type is ``BinaryExpr``, which stores: result - a cached list of linear region indices describing where the expression is valid +- for 3D expressions, any parallel-slice structure needed to evaluate + matching ``yup`` and ``ydown`` expressions -`Field2D`, `Field3D`, and `FieldPerp` act as expression leaves by -providing lightweight ``View`` types. Those views are the device- and -backend-friendly objects used by the expression evaluator. +`Field2D`, `Field3D`, `Field3DParallel`, and `FieldPerp` act as +expression leaves by providing lightweight ``View`` types. Constants are +also wrapped in small view objects so they can participate in the same +expression machinery. These views are the device- and backend-friendly +objects used by the expression evaluator. Materialization happens when a field is constructed or assigned from an expression, when an expression is stored in `Options`, or when a scalar reduction such as ``min`` or ``mean`` is requested. The same mechanism is also used to propagate metadata such as mesh, staggered location, -directions, and `FieldPerp` y-index. +directions, `FieldPerp` y-index, and for `Field3DParallel` the +forward/backward parallel slices on FCI meshes. The unary algebraic helpers in ``include/bout/field.hxx`` build on the same mechanism. Functions such as ``sqrt``, ``abs``, ``SQ``, -``if_else``, ``if_else_zero``, ``min``, ``max``, and ``mean`` can all -operate directly on lazy expressions. +``pow``, ``floor``, ``if_else``, ``if_else_zero``, ``min``, ``max``, +and ``mean`` can all operate directly on lazy expressions. Generated eager operators ~~~~~~~~~~~~~~~~~~~~~~~~~ -The eager arithmetic operators and in-place update paths are still -generated automatically using the `Jinja`_ templating system. The main -files are: +The remaining eager arithmetic operators and in-place update paths are +still generated automatically using the `Jinja`_ templating system. The +main files are: - ``src/field/gen_fieldops.jinja`` - ``src/field/gen_fieldops.py`` - ``src/field/generated_fieldops.cxx`` -The generated code handles the broad matrix of combinations between -`BoutReal`, `Field2D`, `Field3D`, `Field3DParallel`, and `FieldPerp`, -including several mixed-rank and in-place cases where hand-maintaining -all overloads would be error-prone. +The generated code still handles combinations between `BoutReal`, +`Field2D`, `Field3D`, `Field3DParallel`, and `FieldPerp`, especially the +eager `FieldPerp` wrappers and in-place update operators where +hand-maintaining all overloads would be error-prone. More free +non-member arithmetic now lives in the header-defined lazy-expression +path, so the generator covers a smaller subset than it used to. The generated loops now also depend on the configured execution backend. At configure time, the generator is told whether to emit RAJA-based, diff --git a/manual/sphinx/user_docs/algebraic_operators.rst b/manual/sphinx/user_docs/algebraic_operators.rst index b8c40d4dc5..6955ab5580 100644 --- a/manual/sphinx/user_docs/algebraic_operators.rst +++ b/manual/sphinx/user_docs/algebraic_operators.rst @@ -30,7 +30,7 @@ Common operators +------------------------------------------+------------------------------------------------------+ | ``mean(f, allpe=true, region)`` | Mean (optionally over all processes) | +------------------------------------------+------------------------------------------------------+ - | ``pow(lhs, rhs, region)`` | :math:`\mathtt{lhs}^\mathtt{rhs}` | + | ``pow(lhs, rhs[, region])`` | :math:`\mathtt{lhs}^\mathtt{rhs}` | +------------------------------------------+------------------------------------------------------+ | ``SQ(f, region)`` | Square of ``f`` | +------------------------------------------+------------------------------------------------------+ @@ -54,7 +54,8 @@ Common operators +------------------------------------------+------------------------------------------------------+ | ``tanh(f, region)`` | :math:`\tanh(f)` | +------------------------------------------+------------------------------------------------------+ - | ``floor(f, region)`` | Returns a field with the floor of `f` at each point | + | ``floor(f, floor_value[, region])`` | Returns a field where values below ``floor_value`` | + | | are replaced by ``floor_value`` | +------------------------------------------+------------------------------------------------------+ | ``filter(f, n, region)`` | Calculate the amplitude of the Fourier mode in the | | | z-direction with mode number `n` | @@ -84,12 +85,18 @@ Common operators These operators can usually be combined directly in expressions:: Field3D rhs = sqrt(SQ(n) + SQ(T)); + Field3D profile = pow(n + n0, 1.5); Field3D masked = if_else(use_drive, source * profile, sink * profile); BoutReal max_error = max(abs(lhs - rhs), true); Reductions such as ``min``, ``max``, and ``mean`` can operate directly on an expression, so an intermediate field is often unnecessary. +``pow`` also participates in lazy field expressions, including mixed +`Field2D`/`Field3D` cases where the `Field2D` operand is broadcast in +``z``. `FieldPerp` also has ``pow`` overloads for `FieldPerp` with +`FieldPerp` or a scalar. + Region arguments ---------------- @@ -108,6 +115,12 @@ When a region-limited expression is materialized into a field, only the selected region is guaranteed to contain valid values. This is the same performance-oriented convention used by other field operators. +`Field3DParallel` follows the same rule for the main field data. On FCI +meshes, materializing a lazy expression into `Field3DParallel` also +preserves the forward and backward parallel slices when those slices are +available on the expression operands. Materializing the same expression +into plain `Field3D` does not preserve those slices. + Further reading --------------- diff --git a/manual/sphinx/user_docs/field_expressions.rst b/manual/sphinx/user_docs/field_expressions.rst index 62a50988d6..7c82d48884 100644 --- a/manual/sphinx/user_docs/field_expressions.rst +++ b/manual/sphinx/user_docs/field_expressions.rst @@ -5,9 +5,9 @@ Field Expressions BOUT++ field algebra now supports *lazy expressions* for many common operations. Instead of creating a temporary field for every ``+``, -``-``, ``*``, ``/``, ``sqrt`` or ``abs``, BOUT++ can keep the expression -symbolic and evaluate it only when a concrete field or scalar result is -needed. +``-``, ``*``, ``/``, ``pow``, ``sqrt`` or ``abs``, BOUT++ can keep the +expression symbolic and evaluate it only when a concrete field or scalar +result is needed. This keeps ordinary model code readable while reducing temporary allocations and extra loops over the mesh. It is especially helpful for @@ -21,9 +21,9 @@ The following operations can form lazy expressions over `Field2D`, sense: - Arithmetic operators: ``+``, ``-``, ``*``, ``/`` +- Binary algebraic helpers such as ``pow`` and ``floor`` - Unary algebraic operators such as ``sqrt``, ``abs``, ``exp``, ``log``, - ``sin``, ``cos``, ``tan``, ``sinh``, ``cosh``, ``tanh``, ``floor``, - and ``SQ`` + ``sin``, ``cos``, ``tan``, ``sinh``, ``cosh``, ``tanh``, and ``SQ`` - Simple conditionals with ``if_else`` and ``if_else_zero`` - Reductions such as ``min``, ``max``, and ``mean`` @@ -32,7 +32,7 @@ For example:: Field3D n, T; Field3D result; - result = sqrt(SQ(n) + SQ(T)); + result = sqrt(SQ(n) + pow(T + 1.0, 2.0)); The right-hand side can stay lazy until the assignment to ``result``. @@ -93,6 +93,8 @@ Several mixed-type combinations are supported directly: - `Field2D` with `Field3D`: the 2D quantity is broadcast in ``z`` - `FieldPerp` with matching perpendicular data: the operation uses the `FieldPerp` y-index +- `pow` follows the same mixed-rank pattern, so either operand can be + `Field2D` or `Field3D` and the result is a `Field3D` - expressions involving metric components may return `Coordinates::FieldMetric`, which is `Field2D` or `Field3D` depending on how BOUT++ was built @@ -121,6 +123,27 @@ expression or zero:: This is particularly convenient when optional source terms are enabled or disabled by compile-time or run-time logic. +`Field3DParallel` and FCI +------------------------- + +Lazy expressions can also carry parallel-slice information. This matters +when working with `Field3DParallel` on FCI meshes: + +- materializing into `Field3DParallel` preserves ``yup`` and ``ydown`` + slices +- materializing the same expression into plain `Field3D` keeps only the + main field values +- on FCI meshes, operands contributing to a `Field3DParallel` + expression must have compatible parallel slices available + +For example:: + + Field3DParallel f_par, g_par; + Field3DParallel result = sqrt(f_par + g_par); + +On a non-FCI mesh, assigning a lazy expression to `Field3DParallel` +still evaluates the main field, but no parallel slices are retained. + Reductions on expressions ------------------------- From a972f3f349e1a0e5f7075fca45078c8c10eff6b2 Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 30 Jul 2026 07:21:41 +0200 Subject: [PATCH 14/14] Use regionID of parallel slices for parallel slices --- include/bout/fieldops.hxx | 25 ++++++++++++++++++++----- 1 file changed, 20 insertions(+), 5 deletions(-) diff --git a/include/bout/fieldops.hxx b/include/bout/fieldops.hxx index 52d2c82a13..6a809b7893 100644 --- a/include/bout/fieldops.hxx +++ b/include/bout/fieldops.hxx @@ -410,14 +410,29 @@ struct BinaryExpr { } } auto yup(int slice = 0) const { - return BinaryExpr{lhs.yup(slice), rhs.yup(slice), f, - mesh, location, directions, - regionID, indices, yindex}; + return BinaryExpr{ + lhs.yup(slice), + rhs.yup(slice), + f, + mesh, + location, + directions, + mesh->getCommonRegion(lhs.yup(slice).getRegionID(), rhs.yup(slice).getRegionID()), + indices, + yindex}; } auto ydown(int slice = 0) const { return BinaryExpr{ - lhs.ydown(slice), rhs.ydown(slice), f, mesh, location, - directions, regionID, indices, yindex}; + lhs.ydown(slice), + rhs.ydown(slice), + f, + mesh, + location, + directions, + mesh->getCommonRegion(lhs.ydown(slice).getRegionID(), + rhs.ydown(slice).getRegionID()), + indices, + yindex}; } //operator ResT() { return ResT{*this}; }