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/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/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 diff --git a/include/bout/field.hxx b/include/bout/field.hxx index b39a82eb0b..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" @@ -539,53 +541,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")); +} + +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)); +} - ASSERT1(areFieldsCompatible(lhs, rhs)); +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)); } /*! @@ -608,6 +790,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 +837,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 +866,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 @@ -694,40 +889,41 @@ 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 > 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 +938,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 @@ -874,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/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..5562c58b49 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,33 @@ 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; + } + /// 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()); + 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], yup_fields.data(), ydown_fields.data(), + static_cast(numberParallelSlices())}; + } + operator View() const { + return View{const_cast(&data[0]), yup_fields.data(), ydown_fields.data(), + static_cast(numberParallelSlices())}; + } //operator View() const { return View{&data[0]}; } ///////////////////////////////////////////////////////// @@ -856,31 +884,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. @@ -907,7 +910,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"); @@ -1024,6 +1026,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 +1061,44 @@ public: Field3D& asField3D() { return *this; } const Field3D& asField3D() const { return *this; } + struct View { + Field3D::View base; + + 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 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 < 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 < 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)}; } + operator View() const { return View{static_cast(*this)}; } + Field3DParallel& operator*=(const Field3D&); Field3DParallel& operator/=(const Field3D&); Field3DParallel& operator+=(const Field3D&); @@ -1074,6 +1121,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(); @@ -1119,16 +1201,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/include/bout/fieldops.hxx b/include/bout/fieldops.hxx index f36061ceaa..6a809b7893 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 @@ -381,6 +383,58 @@ 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, + 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, + mesh->getCommonRegion(lhs.ydown(slice).getRegionID(), + rhs.ydown(slice).getRegionID()), + indices, + yindex}; + } + //operator ResT() { return ResT{*this}; } struct View { typename L::View lhs; @@ -396,13 +450,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..834c82be4d 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])}; } @@ -405,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/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 ------------------------- diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 8700096e24..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); @@ -936,6 +920,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/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/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 f957652a62..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; } @@ -366,3 +367,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/fake_mesh_fixture.hxx b/tests/unit/fake_mesh_fixture.hxx index 2758dbe416..92d8ee2a13 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; @@ -148,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_field2d.cxx b/tests/unit/field/test_field2d.cxx index a91acd4a40..91721714ec 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; @@ -1317,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 905b182018..6a6ecb9aea 100644 --- a/tests/unit/field/test_field3d.cxx +++ b/tests/unit/field/test_field3d.cxx @@ -7,13 +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 @@ -26,6 +29,7 @@ using namespace bout::globals; // Reuse the "standard" fixture for FakeMesh using Field3DTest = FakeMeshFixture; +using Field3DTestFCI = FakeMeshFixtureFCI; TEST_F(Field3DTest, Is3D) { Field3D field; @@ -74,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); @@ -91,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); @@ -166,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); @@ -293,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); @@ -305,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); } @@ -316,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); } @@ -433,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 @@ -772,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(); @@ -883,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) { @@ -901,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) { @@ -919,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) { @@ -951,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); } @@ -971,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); } @@ -982,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); @@ -994,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); @@ -1006,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); @@ -1027,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; @@ -1056,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]; @@ -1229,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)); } @@ -1259,7 +1264,7 @@ TEST_F(Field3DTest, AssignFromInvalid) { TEST_F(Field3DTest, AssignFromField2D) { Field3D field; - Field2D field2(2.0); + const Field2D field2(2.0); field = field2; @@ -1942,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; @@ -1962,7 +1998,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; @@ -1972,11 +2008,410 @@ 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}; + + // Not FCI so parallel slices are not calculated + EXPECT_FALSE(result.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(result, 4.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(Field3DTestFCI, Field3DParallelAssignmentFromLazyExprRequiresField3DSlicesInFci) { + 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(Field3DTestFCI, Field3DParallelAssignmentFromLazyExprUsesField3DSlicesInFci) { + 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(Field3DTestFCI, Field3DParallelAssignmentFromField3DExprUsesSlicesInFci) { + 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(Field3DTestFCI, Field3DParallelAssignmentFromScalarExprUsesSlicesInFci) { + 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(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, 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; + + 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()); + + field = 2.0; + field.splitParallelSlices(); + field.yup() = 3.0; + field.ydown() = 4.0; + field.resetRegionParallel(); + + Field3D rhs; + EXPECT_TRUE(rhs.isFci()); + rhs = 3.0; + rhs.splitParallelSlices(); + rhs.yup() = 4.0; + rhs.ydown() = 5.0; + rhs.resetRegionParallel(); + + const Field3D prod = field * rhs; + + EXPECT_FALSE(prod.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(prod, 6.0)); + + const Field3DParallel prodpar = field.asField3DParallel() * rhs; + EXPECT_TRUE((std::is_same_v, Field3DParallel>)); + EXPECT_TRUE(prodpar.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(prodpar, 6.0)); + EXPECT_TRUE(IsFieldEqual(prodpar.yup(), 12.0, "RGN_YPAR_+1")); + EXPECT_TRUE(IsFieldEqual(prodpar.ydown(), 20.0, "RGN_YPAR_-1")); +} + +TEST_F(Field3DTestFCI, DivField3DParallelPreservesParallelSlices) { + Field3D field; + + field = 2.0; + field.splitParallelSlices(); + field.yup() = 3.0; + field.ydown() = 4.0; + field.resetRegionParallel(); + + Field3DParallel rhs{1.0}; + EXPECT_TRUE(IsFieldEqual(rhs.ydown(), 1, "RGN_YPAR_-1")); + rhs *= 3.0; + EXPECT_TRUE(IsFieldEqual(rhs.ydown(), 3, "RGN_YPAR_-1")); + rhs.yup() *= 4; + EXPECT_TRUE(IsFieldEqual(rhs.ydown(), 3, "RGN_YPAR_-1")); + rhs.ydown() *= 20. / 3; + EXPECT_TRUE(IsFieldEqual(rhs.ydown(), 20, "RGN_YPAR_-1")); + + const Field3D prod = field / rhs.asField3D(); + + EXPECT_FALSE(prod.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(prod, 2. / 3)); + + const Field3DParallel prodpar = field / rhs; + EXPECT_TRUE(prodpar.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(prodpar, 2. / 3)); + EXPECT_TRUE(IsFieldEqual(prodpar.yup(), .25, "RGN_YPAR_+1")); + EXPECT_TRUE(IsFieldEqual(prodpar.ydown(), 4. / 20, "RGN_YPAR_-1")); +} + +TEST_F(Field3DTestFCI, AddField3DParallelPreservesParallelSlices) { + Field3D field; + field = 2.0; + field.splitParallelSlices(); + field.yup() = 3.0; + field.ydown() = 4.0; + field.resetRegionParallel(); + + Field3D rhs; + rhs = 3.0; + rhs.splitParallelSlices(); + rhs.yup() = 4.0; + rhs.ydown() = 5.0; + rhs.resetRegionParallel(); + + const Field3D res = field + rhs; + + EXPECT_FALSE(res.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(res, 5.0)); + + const Field3DParallel respar = field.asField3DParallel() + rhs; + EXPECT_TRUE((std::is_same_v, Field3DParallel>)); + EXPECT_TRUE(respar.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(respar, 5.0)); + EXPECT_TRUE(IsFieldEqual(respar.yup(), 7.0, "RGN_YPAR_+1")); + EXPECT_TRUE(IsFieldEqual(respar.ydown(), 9.0, "RGN_YPAR_-1")); +} + +TEST_F(Field3DTestFCI, SubField3DParallelPreservesParallelSlices) { + Field3D field; + field = 2.0; + field.splitParallelSlices(); + field.yup() = 3.0; + field.ydown() = 4.0; + field.resetRegionParallel(); + + Field3D rhs; + rhs = 3.0; + rhs.splitParallelSlices(); + rhs.yup() = 5.0; + rhs.ydown() = 7.0; + rhs.resetRegionParallel(); + + const Field3D res = field - rhs; + + EXPECT_FALSE(res.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(res, -1.0)); + + const Field3DParallel respar = field.asField3DParallel() - rhs; + EXPECT_TRUE((std::is_same_v, Field3DParallel>)); + EXPECT_TRUE(respar.hasParallelSlices()); + EXPECT_TRUE(IsFieldEqual(respar, -1.0)); + EXPECT_TRUE(IsFieldEqual(respar.yup(), -2.0, "RGN_YPAR_+1")); + EXPECT_TRUE(IsFieldEqual(respar.ydown(), -3.0, "RGN_YPAR_-1")); } TEST_F(Field3DTest, Abs) { @@ -2108,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_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}; 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);