From c9955c379f99951f19a5d9d1edb859b69d848a1d Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sun, 21 Jun 2026 13:55:08 -0700 Subject: [PATCH 01/18] Stencil expressions: DDZ_C2 and DDZ_C4 operators Return BinaryExpr types, perform derivatives in Z using CoordinateAccessor to get dz. Appear to inline in Hasegawa-Wakatani example. --- examples/hasegawa-wakatani/hw.cxx | 3 +- include/bout/stencil_expr.hxx | 80 +++++++++++++++++++++++++++++++ 2 files changed, 82 insertions(+), 1 deletion(-) create mode 100644 include/bout/stencil_expr.hxx diff --git a/examples/hasegawa-wakatani/hw.cxx b/examples/hasegawa-wakatani/hw.cxx index c3f8717597..7bea5938d8 100644 --- a/examples/hasegawa-wakatani/hw.cxx +++ b/examples/hasegawa-wakatani/hw.cxx @@ -3,6 +3,7 @@ #include #include #include +#include class HW : public PhysicsModel { private: @@ -110,7 +111,7 @@ class HW : public PhysicsModel { } ddt(n) = - -bracket(phi, n, bm) + alpha * (nonzonal_phi - nonzonal_n) - kappa * DDZ(phi); + -bracket(phi, n, bm) + alpha * (nonzonal_phi - nonzonal_n) - kappa * DDZ_C2(phi); ddt(vort) = -bracket(phi, vort, bm) + alpha * (nonzonal_phi - nonzonal_n); diff --git a/include/bout/stencil_expr.hxx b/include/bout/stencil_expr.hxx new file mode 100644 index 0000000000..1e67c4779f --- /dev/null +++ b/include/bout/stencil_expr.hxx @@ -0,0 +1,80 @@ +#pragma once +#ifndef BOUT_STENCIL_EXPR_HXX +#define BOUT_STENCIL_EXPR_HXX + +#include "bout/coordinates_accessor.hxx" +#include "bout/field3d.hxx" +#include "bout/fieldops.hxx" +#include "bout/single_index_ops.hxx" + +namespace bout::stencil { + +struct DDZ_C2_Op { + CoordinatesAccessor coords; + int nz{0}; + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal operator()(int idx, const LView& lhs, + const RView&) const { + const int izp = i_zp(idx, nz); + const int izm = i_zm(idx, nz); + + return 0.5 * (lhs(izp) - lhs(izm)) / coords.dz(idx); + } +}; + +struct DDZ_C4_Op { + CoordinatesAccessor coords; + int nz{0}; + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal operator()(int idx, const LView& lhs, + const RView&) const { + const int izp = i_zp(idx, nz); + const int izm = i_zm(idx, nz); + const int izp2 = i_zp(izp, nz); + const int izm2 = i_zm(izm, nz); + + return (-lhs(izp2) + 8.0 * lhs(izp) - 8.0 * lhs(izm) + lhs(izm2)) + / (12.0 * coords.dz(idx)); + } +}; + +using DDZExprC2 = BinaryExpr; +using DDZExprC4 = BinaryExpr; + +} // namespace bout::stencil + +inline bout::stencil::DDZExprC2 DDZ_C2(const Field3D& f) { + checkData(f); + + const auto region_id = f.getMesh()->getRegionID("RGN_NOBNDRY"); + + return bout::stencil::DDZExprC2{ + static_cast(f), + static_cast(f), + bout::stencil::DDZ_C2_Op{CoordinatesAccessor{f.getCoordinates()}, f.getNz()}, + f.getMesh(), + f.getLocation(), + f.getDirections(), + region_id, + f.getMesh()->getRegion("RGN_NOBNDRY")}; +} + +inline bout::stencil::DDZExprC4 DDZ_C4(const Field3D& f) { + checkData(f); + + const auto region_id = f.getMesh()->getRegionID("RGN_NOBNDRY"); + + return bout::stencil::DDZExprC4{ + static_cast(f), + static_cast(f), + bout::stencil::DDZ_C4_Op{CoordinatesAccessor{f.getCoordinates()}, f.getNz()}, + f.getMesh(), + f.getLocation(), + f.getDirections(), + region_id, + f.getMesh()->getRegion("RGN_NOBNDRY")}; +} + +#endif From acb30aba4e76e239c54ddbf9d5c73b882a76cbae Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Mon, 22 Jun 2026 15:02:03 -0700 Subject: [PATCH 02/18] Stencil expressions: bracket_arakawa Implements Arakawa bracket in X-Z as a BinaryExpr. --- examples/hasegawa-wakatani/hw.cxx | 6 ++-- include/bout/stencil_expr.hxx | 57 +++++++++++++++++++++++++++++++ 2 files changed, 60 insertions(+), 3 deletions(-) diff --git a/examples/hasegawa-wakatani/hw.cxx b/examples/hasegawa-wakatani/hw.cxx index 7bea5938d8..e72d844a00 100644 --- a/examples/hasegawa-wakatani/hw.cxx +++ b/examples/hasegawa-wakatani/hw.cxx @@ -110,10 +110,10 @@ class HW : public PhysicsModel { nonzonal_phi -= averageY(DC(phi)); } - ddt(n) = - -bracket(phi, n, bm) + alpha * (nonzonal_phi - nonzonal_n) - kappa * DDZ_C2(phi); + ddt(n) = -bracket_arakawa(phi, n) + alpha * (nonzonal_phi - nonzonal_n) + - kappa * DDZ_C2(phi); - ddt(vort) = -bracket(phi, vort, bm) + alpha * (nonzonal_phi - nonzonal_n); + ddt(vort) = -bracket_arakawa(phi, vort) + alpha * (nonzonal_phi - nonzonal_n); return 0; } diff --git a/include/bout/stencil_expr.hxx b/include/bout/stencil_expr.hxx index 1e67c4779f..807a1c3852 100644 --- a/include/bout/stencil_expr.hxx +++ b/include/bout/stencil_expr.hxx @@ -5,6 +5,7 @@ #include "bout/coordinates_accessor.hxx" #include "bout/field3d.hxx" #include "bout/fieldops.hxx" +#include "bout/mesh.hxx" #include "bout/single_index_ops.hxx" namespace bout::stencil { @@ -40,8 +41,44 @@ struct DDZ_C4_Op { } }; +struct BracketArakawaOp { + CoordinatesAccessor coords; + int ny{0}; + int nz{0}; + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal operator()(int idx, const LView& f, + const RView& g) const { + const int ixp = i_xp(idx, ny, nz); + const int ixm = i_xm(idx, ny, nz); + const int izp = i_zp(idx, nz); + const int izm = i_zm(idx, nz); + + const int izpxp = i_xp(izp, ny, nz); + const int izpxm = i_xm(izp, ny, nz); + const int izmxp = i_xp(izm, ny, nz); + const int izmxm = i_xm(izm, ny, nz); + + // J++ = DDZ(f)*DDX(g) - DDX(f)*DDZ(g) + const BoutReal Jpp = + ((f(izp) - f(izm)) * (g(ixp) - g(ixm)) - (f(ixp) - f(ixm)) * (g(izp) - g(izm))); + + // J+x + const BoutReal Jpx = + (g(ixp) * (f(izpxp) - f(izmxp)) - g(ixm) * (f(izpxm) - f(izmxm)) + - g(izp) * (f(izpxp) - f(izpxm)) + g(izm) * (f(izmxp) - f(izmxm))); + + // Jx+ + const BoutReal Jxp = (g(izpxp) * (f(izp) - f(ixp)) - g(izmxm) * (f(ixm) - f(izm)) + - g(izpxm) * (f(izp) - f(ixm)) + g(izmxp) * (f(ixp) - f(izm))); + + return (Jpp + Jpx + Jxp) / (12.0 * coords.dx(idx) * coords.dz(idx)); + } +}; + using DDZExprC2 = BinaryExpr; using DDZExprC4 = BinaryExpr; +using BracketArakawaExpr = BinaryExpr; } // namespace bout::stencil @@ -77,4 +114,24 @@ inline bout::stencil::DDZExprC4 DDZ_C4(const Field3D& f) { f.getMesh()->getRegion("RGN_NOBNDRY")}; } +inline bout::stencil::BracketArakawaExpr bracket_arakawa(const Field3D& f, + const Field3D& g) { + checkData(f); + checkData(g); + ASSERT1_FIELDS_COMPATIBLE(f, g); + + const auto region_id = f.getMesh()->getRegionID("RGN_NOBNDRY"); + + return bout::stencil::BracketArakawaExpr{ + static_cast(f), + static_cast(g), + bout::stencil::BracketArakawaOp{CoordinatesAccessor{f.getCoordinates()}, f.getNy(), + f.getNz()}, + f.getMesh(), + f.getLocation(), + f.getDirections(), + region_id, + f.getMesh()->getRegion("RGN_NOBNDRY")}; +} + #endif From 4fd34afd448135c0af48761face2e6cda99fb9bf Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Mon, 22 Jun 2026 22:54:25 -0700 Subject: [PATCH 03/18] stencil expr: DDX_Dispatch Testing the performance impact of runtime dispatch. The hope is to preserve the ability to switch method at runtime. The conditional that selects the method is the same for all iterations, so hopefully good branch prediction and no warp divergence. DDZ_Dispatch(f, method) selects between DDZ_C2 and DDZ_C4 at runtime. --- examples/hasegawa-wakatani/hw.cxx | 2 +- include/bout/stencil_expr.hxx | 61 +++++++++++++++++++++++++++++++ 2 files changed, 62 insertions(+), 1 deletion(-) diff --git a/examples/hasegawa-wakatani/hw.cxx b/examples/hasegawa-wakatani/hw.cxx index e72d844a00..39ffdd709c 100644 --- a/examples/hasegawa-wakatani/hw.cxx +++ b/examples/hasegawa-wakatani/hw.cxx @@ -111,7 +111,7 @@ class HW : public PhysicsModel { } ddt(n) = -bracket_arakawa(phi, n) + alpha * (nonzonal_phi - nonzonal_n) - - kappa * DDZ_C2(phi); + - kappa * DDZ_Dispatch(phi, DIFF_C2); ddt(vort) = -bracket_arakawa(phi, vort) + alpha * (nonzonal_phi - nonzonal_n); diff --git a/include/bout/stencil_expr.hxx b/include/bout/stencil_expr.hxx index 807a1c3852..9820012fa2 100644 --- a/include/bout/stencil_expr.hxx +++ b/include/bout/stencil_expr.hxx @@ -41,6 +41,44 @@ struct DDZ_C4_Op { } }; +struct DDZ_Dispatch_Op { + CoordinatesAccessor coords; + int nz{0}; + DIFF_METHOD method{DIFF_DEFAULT}; + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c2(int idx, const LView& lhs) const { + const int izp = i_zp(idx, nz); + const int izm = i_zm(idx, nz); + + return 0.5 * (lhs(izp) - lhs(izm)) / coords.dz(idx); + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c4(int idx, const LView& lhs) const { + const int izp = i_zp(idx, nz); + const int izm = i_zm(idx, nz); + const int izp2 = i_zp(izp, nz); + const int izm2 = i_zm(izm, nz); + + return (-lhs(izp2) + 8.0 * lhs(izp) - 8.0 * lhs(izm) + lhs(izm2)) + / (12.0 * coords.dz(idx)); + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal operator()(int idx, const LView& lhs, + const RView&) const { + switch (method) { + case DIFF_C2: + return apply_c2(idx, lhs); + case DIFF_C4: + return apply_c4(idx, lhs); + default: + return 0.0; + } + } +}; + struct BracketArakawaOp { CoordinatesAccessor coords; int ny{0}; @@ -78,6 +116,7 @@ struct BracketArakawaOp { using DDZExprC2 = BinaryExpr; using DDZExprC4 = BinaryExpr; +using DDZDispatchExpr = BinaryExpr; using BracketArakawaExpr = BinaryExpr; } // namespace bout::stencil @@ -114,6 +153,28 @@ inline bout::stencil::DDZExprC4 DDZ_C4(const Field3D& f) { f.getMesh()->getRegion("RGN_NOBNDRY")}; } +inline bout::stencil::DDZDispatchExpr DDZ_Dispatch(const Field3D& f, DIFF_METHOD method) { + checkData(f); + + if ((method != DIFF_C2) && (method != DIFF_C4)) { + throw BoutException("DDZ_Dispatch only supports DIFF_C2 and DIFF_C4, got {:s}", + toString(method)); + } + + const auto region_id = f.getMesh()->getRegionID("RGN_NOBNDRY"); + + return bout::stencil::DDZDispatchExpr{ + static_cast(f), + static_cast(f), + bout::stencil::DDZ_Dispatch_Op{CoordinatesAccessor{f.getCoordinates()}, f.getNz(), + method}, + f.getMesh(), + f.getLocation(), + f.getDirections(), + region_id, + f.getMesh()->getRegion("RGN_NOBNDRY")}; +} + inline bout::stencil::BracketArakawaExpr bracket_arakawa(const Field3D& f, const Field3D& g) { checkData(f); From 1238f4b908a2f20449ae83b2cea5ee0bf6854e03 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 23 Jun 2026 12:36:42 -0700 Subject: [PATCH 04/18] pow: Convert to BinaryExpr lazy function --- include/bout/field.hxx | 163 ++++++++++++++++++---------- include/bout/field3d.hxx | 19 +++- include/bout/fieldops.hxx | 10 ++ src/field/field3d.cxx | 16 --- tests/unit/field/test_field2d.cxx | 25 +++++ tests/unit/field/test_field3d.cxx | 26 +++++ tests/unit/field/test_fieldperp.cxx | 27 +++++ 7 files changed, 213 insertions(+), 73 deletions(-) diff --git a/include/bout/field.hxx b/include/bout/field.hxx index a8e1952546..51281d08de 100644 --- a/include/bout/field.hxx +++ b/include/bout/field.hxx @@ -539,53 +539,127 @@ inline BoutReal mean(const BinaryExpr& f, bool allpe = false, return bout::reduce::Mean::finalize(state); } -/// 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") { - - ASSERT1(areFieldsCompatible(lhs, rhs)); +class Field3DParallel; +class FieldPerp; - T result{emptyFrom(lhs)}; +namespace bout::detail { +template +struct expression_result { + using type = std::decay_t; +}; - BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs[i], rhs[i]); } +template +struct expression_result> { + using type = ResT; +}; - checkData(result); - return result; -} +template +using expression_result_t = typename expression_result>::type; -template > -T pow(const T& lhs, BoutReal rhs, const std::string& rgn = "RGN_ALL") { +template +inline constexpr bool is_expression_field_v = + is_expr_field2d_v || is_expr_field3d_v || is_expr_fieldperp_v; - // Check if the inputs are allocated - checkData(lhs); - checkData(rhs); +template +inline constexpr bool is_same_expression_rank_v = + (is_expr_field2d_v && is_expr_field2d_v) + || (is_expr_field3d_v && is_expr_field3d_v) + || (is_expr_fieldperp_v && is_expr_fieldperp_v); - T result{emptyFrom(lhs)}; +template +std::optional getPerpYIndex(const T& value) { + if constexpr (std::is_same_v, ::FieldPerp>) { + return value.getIndex(); + } else { + return std::nullopt; + } +} - BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs[i], rhs); } +template +std::optional getPerpYIndex(const BinaryExpr& expr) { + if constexpr (std::is_same_v) { + return expr.getIndex(); + } else { + return std::nullopt; + } +} - checkData(result); - return result; +template +std::optional getExpressionRegionID(const Expr& expr, + const std::string& region_name) { + std::optional region_id{}; + if (expr.getMesh()->hasRegion3D(region_name)) { + region_id = expr.getMesh()->getRegionID(region_name); + } + return expr.getMesh()->getCommonRegion(region_id, expr.getRegionID()); } -template > -T pow(BoutReal lhs, const T& rhs, const std::string& rgn = "RGN_ALL") { +template +std::optional getExpressionRegionID(const L& lhs, const R& rhs, + const std::string& region_name) { + return lhs.getMesh()->getCommonRegion(getExpressionRegionID(lhs, region_name), + rhs.getRegionID()); +} +} // namespace bout::detail - // Check if the inputs are allocated - checkData(lhs); - checkData(rhs); +/// 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 , + typename = std::enable_if_t>> +auto pow(const L& lhs, const R& rhs, const std::string& rgn = "RGN_ALL") { - // Define and allocate the output result - T result{emptyFrom(rhs)}; + if constexpr (bout::utils::is_Field_v> + && bout::utils::is_Field_v>) { + ASSERT1(areFieldsCompatible(lhs, rhs)); + } else { + ASSERT1_EXPR_COMPATIBLE(lhs, rhs); + } + + return BinaryExpr{ + static_cast(lhs), + static_cast(rhs), + bout::op::Pow{}, + lhs.getMesh(), + lhs.getLocation(), + lhs.getDirections(), + bout::detail::getExpressionRegionID(lhs, rhs, rgn), + lhs.getMesh()->template getRegion(rgn), + bout::detail::getPerpYIndex(lhs)}; +} - BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs, rhs[i]); } +template , + typename = std::enable_if_t>> +auto pow(const T& lhs, BoutReal rhs, const std::string& rgn = "RGN_ALL") { + + return BinaryExpr, bout::op::Pow>{ + static_cast(lhs), + static_cast::View>(rhs), + bout::op::Pow{}, + lhs.getMesh(), + lhs.getLocation(), + lhs.getDirections(), + bout::detail::getExpressionRegionID(lhs, rgn), + lhs.getMesh()->template getRegion(rgn), + bout::detail::getPerpYIndex(lhs)}; +} - checkData(result); - return result; +template , + typename = std::enable_if_t>> +auto pow(BoutReal lhs, const T& rhs, const std::string& rgn = "RGN_ALL") { + + return BinaryExpr, T, bout::op::Pow>{ + static_cast::View>(lhs), + static_cast(rhs), + bout::op::Pow{}, + rhs.getMesh(), + rhs.getLocation(), + rhs.getDirections(), + bout::detail::getExpressionRegionID(rhs, rgn), + rhs.getMesh()->template getRegion(rgn), + bout::detail::getPerpYIndex(rhs)}; } /*! @@ -604,29 +678,6 @@ T pow(BoutReal lhs, const T& rhs, const std::string& rgn = "RGN_ALL") { * result for non-finite numbers * */ -class Field3DParallel; -class FieldPerp; - -namespace bout::detail { -template -std::optional getPerpYIndex(const T& value) { - if constexpr (std::is_same_v, ::FieldPerp>) { - return value.getIndex(); - } else { - return std::nullopt; - } -} - -template -std::optional getPerpYIndex(const BinaryExpr& expr) { - if constexpr (std::is_same_v) { - return expr.getIndex(); - } else { - return std::nullopt; - } -} -} // namespace bout::detail - #ifdef FIELD_FUNC #error This macro has already been defined #else diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index 09c7f1ce55..abcb52f892 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -900,7 +900,24 @@ 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"); +template +std::enable_if_t && is_expr_field2d_v, + BinaryExpr> +pow(const L& lhs, const R& rhs, const std::string& rgn = "RGN_ALL") { + ASSERT1_EXPR_COMPATIBLE(lhs, rhs); + auto regionID = bout::detail::getExpressionRegionID(lhs, rhs, rgn); + int mesh_nz = lhs.getMesh()->LocalNz; + return BinaryExpr{ + static_cast(lhs), + static_cast(rhs).setScale(1, mesh_nz), + bout::op::Pow{}, + lhs.getMesh(), + lhs.getLocation(), + lhs.getDirections(), + regionID, + (regionID.has_value() ? lhs.getMesh()->getRegion(regionID.value()) + : lhs.getMesh()->getRegion(rgn))}; +} FieldPerp pow(const Field3D& lhs, const FieldPerp& rhs, const std::string& rgn = "RGN_ALL"); diff --git a/include/bout/fieldops.hxx b/include/bout/fieldops.hxx index 53e28042c9..1d3193f638 100644 --- a/include/bout/fieldops.hxx +++ b/include/bout/fieldops.hxx @@ -118,6 +118,16 @@ struct Div { return a / b; } }; +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); + } +}; struct IfElse { bool condition; diff --git a/src/field/field3d.cxx b/src/field/field3d.cxx index 7915579440..0d7687d0ca 100644 --- a/src/field/field3d.cxx +++ b/src/field/field3d.cxx @@ -688,22 +688,6 @@ void Field3D::swapData(Field3D& other) { std::swap(data, other.data); } //////////////// NON-MEMBER FUNCTIONS ////////////////// -Field3D pow(const Field3D& lhs, const Field2D& rhs, const std::string& rgn) { - - // Check if the inputs are allocated - checkData(lhs); - checkData(rhs); - ASSERT1_FIELDS_COMPATIBLE(lhs, rhs); - - // Define and allocate the output result - Field3D result{emptyFrom(lhs)}; - - BOUT_FOR(i, result.getRegion(rgn)) { result[i] = ::pow(lhs[i], rhs[i]); } - - checkData(result); - return result; -} - FieldPerp pow(const Field3D& lhs, const FieldPerp& rhs, const std::string& rgn) { checkData(lhs); diff --git a/tests/unit/field/test_field2d.cxx b/tests/unit/field/test_field2d.cxx index a91acd4a40..a5920d5d19 100644 --- a/tests/unit/field/test_field2d.cxx +++ b/tests/unit/field/test_field2d.cxx @@ -1168,6 +1168,31 @@ TEST_F(Field2DTest, PowField2DField2D) { EXPECT_TRUE(IsFieldEqual(c, 64.0)); } +TEST_F(Field2DTest, PowExpressionUsesBinaryExpr) { + Field2D field; + + field = 2.0; + const auto expr = field + 1.0; + + EXPECT_TRUE( + (std::is_same_v, + BinaryExpr, bout::op::Pow>>)); + EXPECT_TRUE((std::is_same_v, + BinaryExpr, + Constant, bout::op::Pow>>)); + EXPECT_TRUE(IsFieldEqual(pow(expr, 2.0), 9.0)); +} + +TEST_F(Field2DTest, PowRegionLimitedExpressionConstructsField2D) { + Field2D field; + + field = 2.0; + + Field2D result = pow(field, 2.0, "RGN_NOBNDRY"); + + EXPECT_TRUE(IsFieldEqual(result, 4.0, "RGN_NOBNDRY")); +} + TEST_F(Field2DTest, Sqrt) { Field2D field; diff --git a/tests/unit/field/test_field3d.cxx b/tests/unit/field/test_field3d.cxx index 905b182018..8502837642 100644 --- a/tests/unit/field/test_field3d.cxx +++ b/tests/unit/field/test_field3d.cxx @@ -1942,6 +1942,32 @@ TEST_F(Field3DTest, PowField3DField3D) { EXPECT_TRUE(IsFieldEqual(c, 64.0)); } +TEST_F(Field3DTest, PowExpressionUsesBinaryExpr) { + Field3D field; + + field = 2.0; + const auto expr = field + 1.0; + + EXPECT_TRUE( + (std::is_same_v, + BinaryExpr, bout::op::Pow>>)); + EXPECT_TRUE((std::is_same_v, + BinaryExpr, + Constant, bout::op::Pow>>)); + EXPECT_TRUE(IsFieldEqual(pow(expr, 2.0), 9.0)); +} + +TEST_F(Field3DTest, PowRegionArgumentSetsRegionID) { + Field3D field; + + field = 2.0; + const auto expr = pow(field, 2.0, "RGN_NOBNDRY"); + + ASSERT_TRUE(expr.getRegionID().has_value()); + EXPECT_EQ(expr.getRegionID().value(), field.getMesh()->getRegionID("RGN_NOBNDRY")); + EXPECT_TRUE(IsFieldEqual(expr, 4.0, "RGN_NOBNDRY")); +} + TEST_F(Field3DTest, Sqrt) { Field3D field; diff --git a/tests/unit/field/test_fieldperp.cxx b/tests/unit/field/test_fieldperp.cxx index 46f07d589f..60a15ea693 100644 --- a/tests/unit/field/test_fieldperp.cxx +++ b/tests/unit/field/test_fieldperp.cxx @@ -1569,6 +1569,33 @@ TEST_F(FieldPerpTest, PowFieldPerpFieldPerp) { EXPECT_TRUE(IsFieldEqual(c, 64.0)); } +TEST_F(FieldPerpTest, PowExpressionUsesBinaryExpr) { + FieldPerp field; + field.setIndex(0); + + field = 2.0; + const auto expr = field + 1.0; + + EXPECT_TRUE((std::is_same_v< + std::decay_t, + BinaryExpr, bout::op::Pow>>)); + EXPECT_TRUE((std::is_same_v, + BinaryExpr, + Constant, bout::op::Pow>>)); + EXPECT_TRUE(IsFieldEqual(pow(expr, 2.0), 9.0)); +} + +TEST_F(FieldPerpTest, PowRegionLimitedExpressionConstructsFieldPerp) { + FieldPerp field; + field.setIndex(0); + + field = 2.0; + + FieldPerp result = pow(field, 2.0, "RGN_NOBNDRY"); + + EXPECT_TRUE(IsFieldEqual(result, 4.0, "RGN_NOBNDRY")); +} + TEST_F(FieldPerpTest, Sqrt) { FieldPerp field; field.setIndex(0); From 19502b404423558f876b8a48e3141191f879daf0 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 23 Jun 2026 13:15:32 -0700 Subject: [PATCH 05/18] DDZ_Dispatch: Extend to handle staggered locations Retains the runtime dispatch of DDZ, while building a BinaryExpr lazy expression. --- examples/hasegawa-wakatani/hw.cxx | 2 +- include/bout/stencil_expr.hxx | 96 +++++++++++++++++-- tests/unit/CMakeLists.txt | 1 + tests/unit/include/bout/test_stencil_expr.cxx | 78 +++++++++++++++ 4 files changed, 168 insertions(+), 9 deletions(-) create mode 100644 tests/unit/include/bout/test_stencil_expr.cxx diff --git a/examples/hasegawa-wakatani/hw.cxx b/examples/hasegawa-wakatani/hw.cxx index 39ffdd709c..42dd923c59 100644 --- a/examples/hasegawa-wakatani/hw.cxx +++ b/examples/hasegawa-wakatani/hw.cxx @@ -111,7 +111,7 @@ class HW : public PhysicsModel { } ddt(n) = -bracket_arakawa(phi, n) + alpha * (nonzonal_phi - nonzonal_n) - - kappa * DDZ_Dispatch(phi, DIFF_C2); + - kappa * DDZ_Dispatch(phi, CELL_DEFAULT, DIFF_C2); ddt(vort) = -bracket_arakawa(phi, vort) + alpha * (nonzonal_phi - nonzonal_n); diff --git a/include/bout/stencil_expr.hxx b/include/bout/stencil_expr.hxx index 9820012fa2..7e031305f4 100644 --- a/include/bout/stencil_expr.hxx +++ b/include/bout/stencil_expr.hxx @@ -2,12 +2,18 @@ #ifndef BOUT_STENCIL_EXPR_HXX #define BOUT_STENCIL_EXPR_HXX +#include "bout/bout_types.hxx" +#include "bout/boutexception.hxx" +#include "bout/build_config.hxx" #include "bout/coordinates_accessor.hxx" +#include "bout/field.hxx" #include "bout/field3d.hxx" #include "bout/fieldops.hxx" #include "bout/mesh.hxx" #include "bout/single_index_ops.hxx" +#include + namespace bout::stencil { struct DDZ_C2_Op { @@ -44,10 +50,12 @@ struct DDZ_C4_Op { struct DDZ_Dispatch_Op { CoordinatesAccessor coords; int nz{0}; + STAGGER stagger{STAGGER::None}; DIFF_METHOD method{DIFF_DEFAULT}; template - BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c2(int idx, const LView& lhs) const { + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c2_none(int idx, + const LView& lhs) const { const int izp = i_zp(idx, nz); const int izm = i_zm(idx, nz); @@ -55,7 +63,8 @@ struct DDZ_Dispatch_Op { } template - BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c4(int idx, const LView& lhs) const { + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c4_none(int idx, + const LView& lhs) const { const int izp = i_zp(idx, nz); const int izm = i_zm(idx, nz); const int izp2 = i_zp(izp, nz); @@ -65,6 +74,70 @@ struct DDZ_Dispatch_Op { / (12.0 * coords.dz(idx)); } + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c2_c2l(int idx, + const LView& lhs) const { + const int izm = i_zm(idx, nz); + + return (lhs(idx) - lhs(izm)) / coords.dz(idx); + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c4_c2l(int idx, + const LView& lhs) const { + const int izm = i_zm(idx, nz); + const int izp = i_zp(idx, nz); + const int izm2 = i_zm(izm, nz); + + return (27.0 * (lhs(idx) - lhs(izm)) - (lhs(izp) - lhs(izm2))) + / (24.0 * coords.dz(idx)); + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c2_l2c(int idx, + const LView& lhs) const { + const int izp = i_zp(idx, nz); + + return (lhs(izp) - lhs(idx)) / coords.dz(idx); + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c4_l2c(int idx, + const LView& lhs) const { + const int izp = i_zp(idx, nz); + const int izm = i_zm(idx, nz); + const int izp2 = i_zp(izp, nz); + + return (27.0 * (lhs(izp) - lhs(idx)) - (lhs(izp2) - lhs(izm))) + / (24.0 * coords.dz(idx)); + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c2(int idx, const LView& lhs) const { + switch (stagger) { + case STAGGER::None: + return apply_c2_none(idx, lhs); + case STAGGER::C2L: + return apply_c2_c2l(idx, lhs); + case STAGGER::L2C: + return apply_c2_l2c(idx, lhs); + } + return 0.0; + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c4(int idx, const LView& lhs) const { + switch (stagger) { + case STAGGER::None: + return apply_c4_none(idx, lhs); + case STAGGER::C2L: + return apply_c4_c2l(idx, lhs); + case STAGGER::L2C: + return apply_c4_l2c(idx, lhs); + } + return 0.0; + } + template BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal operator()(int idx, const LView& lhs, const RView&) const { @@ -153,7 +226,10 @@ inline bout::stencil::DDZExprC4 DDZ_C4(const Field3D& f) { f.getMesh()->getRegion("RGN_NOBNDRY")}; } -inline bout::stencil::DDZDispatchExpr DDZ_Dispatch(const Field3D& f, DIFF_METHOD method) { +inline bout::stencil::DDZDispatchExpr +DDZ_Dispatch(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT, + DIFF_METHOD method = DIFF_DEFAULT, + const std::string& region = "RGN_NOBNDRY") { checkData(f); if ((method != DIFF_C2) && (method != DIFF_C4)) { @@ -161,18 +237,22 @@ inline bout::stencil::DDZDispatchExpr DDZ_Dispatch(const Field3D& f, DIFF_METHOD toString(method)); } - const auto region_id = f.getMesh()->getRegionID("RGN_NOBNDRY"); + const auto resolved_outloc = (outloc == CELL_DEFAULT) ? f.getLocation() : outloc; + const auto stagger = + f.getMesh()->getStagger(f.getLocation(), resolved_outloc, CELL_ZLOW); + const auto region_id = f.getMesh()->getRegionID(region); return bout::stencil::DDZDispatchExpr{ static_cast(f), static_cast(f), - bout::stencil::DDZ_Dispatch_Op{CoordinatesAccessor{f.getCoordinates()}, f.getNz(), - method}, + bout::stencil::DDZ_Dispatch_Op{ + CoordinatesAccessor{f.getCoordinates(resolved_outloc)}, f.getNz(), stagger, + method}, f.getMesh(), - f.getLocation(), + resolved_outloc, f.getDirections(), region_id, - f.getMesh()->getRegion("RGN_NOBNDRY")}; + f.getMesh()->getRegion(region)}; } inline bout::stencil::BracketArakawaExpr bracket_arakawa(const Field3D& f, diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index 4963aaf1f0..0e87a4d406 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -80,6 +80,7 @@ set(serial_tests_source ./include/bout/test_petsc_vector.cxx ./include/bout/test_region.cxx ./include/bout/test_single_index_ops.cxx + ./include/bout/test_stencil_expr.cxx ./include/bout/test_stencil.cxx ./include/bout/test_template_combinations.cxx ./include/bout/test_traits.cxx diff --git a/tests/unit/include/bout/test_stencil_expr.cxx b/tests/unit/include/bout/test_stencil_expr.cxx new file mode 100644 index 0000000000..898d9d5694 --- /dev/null +++ b/tests/unit/include/bout/test_stencil_expr.cxx @@ -0,0 +1,78 @@ +#include "gtest/gtest.h" + +#include "test_extras.hxx" + +#include "bout/derivs.hxx" +#include "bout/stencil_expr.hxx" + +#include + +#include "fake_mesh_fixture.hxx" + +namespace { + +Field3D makeTestField(Mesh* mesh, CELL_LOC location) { + Field3D result(mesh, location); + result.allocate(); + + for (auto i : result.getRegion("RGN_ALL")) { + result[i] = 0.1 * i.x() + 0.2 * i.y() + 0.3 * i.z() + 0.05 * i.x() * i.z(); + } + + return result; +} + +class DDZDispatchExprTest : public FakeMeshFixture {}; + +class DDZDispatchExprParamTest + : public DDZDispatchExprTest, + public ::testing::WithParamInterface> { +}; + +std::string paramToString( + const ::testing::TestParamInfo>& param) { + const auto [inloc, outloc, method] = param.param; + return toString(inloc) + "_to_" + toString(outloc) + "_" + toString(method); +} + +} // namespace + +INSTANTIATE_TEST_SUITE_P( + SupportedLocations, DDZDispatchExprParamTest, + ::testing::Values(std::make_tuple(CELL_CENTRE, CELL_CENTRE, DIFF_C2), + std::make_tuple(CELL_CENTRE, CELL_CENTRE, DIFF_C4), + std::make_tuple(CELL_CENTRE, CELL_ZLOW, DIFF_C2), + std::make_tuple(CELL_CENTRE, CELL_ZLOW, DIFF_C4), + std::make_tuple(CELL_ZLOW, CELL_CENTRE, DIFF_C2), + std::make_tuple(CELL_ZLOW, CELL_CENTRE, DIFF_C4), + std::make_tuple(CELL_ZLOW, CELL_ZLOW, DIFF_C2), + std::make_tuple(CELL_ZLOW, CELL_ZLOW, DIFF_C4)), + paramToString); + +TEST_P(DDZDispatchExprParamTest, MatchesDDZ) { + const auto [inloc, outloc, method] = GetParam(); + auto input = makeTestField(mesh_staggered, inloc); + + const auto actual = Field3D{DDZ_Dispatch(input, outloc, method)}; + const auto expected = DDZ(input, outloc, toString(method)); + + EXPECT_EQ(actual.getLocation(), outloc); + EXPECT_TRUE(IsFieldEqual(actual, expected, "RGN_NOBNDRY")); +} + +TEST_F(DDZDispatchExprTest, UsesRequestedRegion) { + auto input = makeTestField(mesh_staggered, CELL_CENTRE); + + const auto actual = Field3D{DDZ_Dispatch(input, CELL_ZLOW, DIFF_C2, "RGN_ALL")}; + const auto expected = DDZ(input, CELL_ZLOW, toString(DIFF_C2), "RGN_ALL"); + + EXPECT_EQ(actual.getLocation(), CELL_ZLOW); + EXPECT_TRUE(IsFieldEqual(actual, expected, "RGN_ALL")); +} + +TEST_F(DDZDispatchExprTest, RejectsUnsupportedMethods) { + auto input = makeTestField(mesh_staggered, CELL_CENTRE); + + EXPECT_THROW((void)DDZ_Dispatch(input, CELL_DEFAULT, DIFF_DEFAULT), BoutException); + EXPECT_THROW((void)DDZ_Dispatch(input, CELL_DEFAULT, DIFF_FFT), BoutException); +} From 6c15db49d40b4661eb62012ae17925df31454ab0 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 23 Jun 2026 14:24:08 -0700 Subject: [PATCH 06/18] Mesh::DerivativeDefaults default method Stores the default numerical methods as DIFF_METHOD enums. This enables runtime dispatch with minimal overhead in expressions. --- examples/hasegawa-wakatani/hw.cxx | 2 +- include/bout/mesh.hxx | 97 +++++++++++++++++++ include/bout/stencil_expr.hxx | 11 ++- src/mesh/index_derivs.cxx | 94 ++++++++++++++++++ tests/unit/include/bout/test_stencil_expr.cxx | 17 +++- tests/unit/mesh/test_mesh.cxx | 38 ++++++++ 6 files changed, 254 insertions(+), 5 deletions(-) diff --git a/examples/hasegawa-wakatani/hw.cxx b/examples/hasegawa-wakatani/hw.cxx index 42dd923c59..f4e546908f 100644 --- a/examples/hasegawa-wakatani/hw.cxx +++ b/examples/hasegawa-wakatani/hw.cxx @@ -111,7 +111,7 @@ class HW : public PhysicsModel { } ddt(n) = -bracket_arakawa(phi, n) + alpha * (nonzonal_phi - nonzonal_n) - - kappa * DDZ_Dispatch(phi, CELL_DEFAULT, DIFF_C2); + - kappa * DDZ_Dispatch(phi); ddt(vort) = -bracket_arakawa(phi, vort) + alpha * (nonzonal_phi - nonzonal_n); diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index 929e9d5aa7..2cef181157 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -59,6 +59,7 @@ class Mesh; #include "bout/sys/range.hxx" // RangeIterator #include "bout/unused.hxx" +#include #include #include #include @@ -95,6 +96,94 @@ using comm_handle = void*; class Mesh { public: + struct DerivativeDefaults { + static constexpr int num_directions = 3; + static constexpr int num_staggers = 3; + static constexpr int num_deriv_kinds = 5; + + std::array values{}; + + DerivativeDefaults() { + for (auto direction : {DIRECTION::X, DIRECTION::Y, DIRECTION::Z}) { + for (auto stagger : {STAGGER::None, STAGGER::C2L, STAGGER::L2C}) { + for (auto deriv : {DERIV::Standard, DERIV::StandardSecond, + DERIV::StandardFourth, DERIV::Upwind, DERIV::Flux}) { + set(direction, deriv, builtinDefaultMethod(deriv), stagger); + } + } + } + } + + static DIFF_METHOD builtinDefaultMethod(DERIV deriv) { + switch (deriv) { + case DERIV::Standard: + case DERIV::StandardSecond: + case DERIV::StandardFourth: + return DIFF_C2; + case DERIV::Upwind: + case DERIV::Flux: + return DIFF_U1; + } + throw BoutException("Unhandled derivative kind in builtinDefaultMethod"); + } + + static int directionIndex(DIRECTION direction) { + switch (direction) { + case DIRECTION::X: + return 0; + case DIRECTION::Y: + case DIRECTION::YOrthogonal: + case DIRECTION::YAligned: + return 1; + case DIRECTION::Z: + return 2; + } + throw BoutException("Unhandled direction in DerivativeDefaults"); + } + + static int staggerIndex(STAGGER stagger) { + switch (stagger) { + case STAGGER::None: + return 0; + case STAGGER::C2L: + return 1; + case STAGGER::L2C: + return 2; + } + throw BoutException("Unhandled stagger in DerivativeDefaults"); + } + + static int derivIndex(DERIV deriv) { + switch (deriv) { + case DERIV::Standard: + return 0; + case DERIV::StandardSecond: + return 1; + case DERIV::StandardFourth: + return 2; + case DERIV::Upwind: + return 3; + case DERIV::Flux: + return 4; + } + throw BoutException("Unhandled derivative kind in DerivativeDefaults"); + } + + DIFF_METHOD get(DIRECTION direction, DERIV deriv, + STAGGER stagger = STAGGER::None) const { + return values[(directionIndex(direction) * num_staggers + staggerIndex(stagger)) + * num_deriv_kinds + + derivIndex(deriv)]; + } + + void set(DIRECTION direction, DERIV deriv, DIFF_METHOD method, + STAGGER stagger = STAGGER::None) { + values[(directionIndex(direction) * num_staggers + staggerIndex(stagger)) + * num_deriv_kinds + + derivIndex(deriv)] = method; + } + }; + /// Constructor for a "bare", uninitialised Mesh /// Only useful for testing Mesh() : source(nullptr), options(nullptr), include_corner_cells(true) {} @@ -719,6 +808,14 @@ public: /// Fraction of modes to filter. This is set in derivs_init from option "ddz:fft_filter" BoutReal fft_derivs_filter{0.0}; + /// Concrete default methods initialised from options in derivs_init + DerivativeDefaults derivative_defaults{}; + + DIFF_METHOD getDefaultMethod(DIRECTION direction, DERIV deriv, + STAGGER stagger = STAGGER::None) const { + return derivative_defaults.get(direction, deriv, stagger); + } + /// Determines the resultant output stagger location in derivatives /// given the input and output location. Also checks that the /// combination of locations is allowed diff --git a/include/bout/stencil_expr.hxx b/include/bout/stencil_expr.hxx index 7e031305f4..bbccc267bf 100644 --- a/include/bout/stencil_expr.hxx +++ b/include/bout/stencil_expr.hxx @@ -232,14 +232,19 @@ DDZ_Dispatch(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& region = "RGN_NOBNDRY") { checkData(f); + const auto resolved_outloc = (outloc == CELL_DEFAULT) ? f.getLocation() : outloc; + const auto stagger = + f.getMesh()->getStagger(f.getLocation(), resolved_outloc, CELL_ZLOW); + + if (method == DIFF_DEFAULT) { + method = f.getMesh()->getDefaultMethod(DIRECTION::Z, DERIV::Standard, stagger); + } + if ((method != DIFF_C2) && (method != DIFF_C4)) { throw BoutException("DDZ_Dispatch only supports DIFF_C2 and DIFF_C4, got {:s}", toString(method)); } - const auto resolved_outloc = (outloc == CELL_DEFAULT) ? f.getLocation() : outloc; - const auto stagger = - f.getMesh()->getStagger(f.getLocation(), resolved_outloc, CELL_ZLOW); const auto region_id = f.getMesh()->getRegionID(region); return bout::stencil::DDZDispatchExpr{ diff --git a/src/mesh/index_derivs.cxx b/src/mesh/index_derivs.cxx index 70fc47b538..f01df89b68 100644 --- a/src/mesh/index_derivs.cxx +++ b/src/mesh/index_derivs.cxx @@ -23,10 +23,58 @@ #include "bout/build_defines.hxx" #include "bout/traits.hxx" +#include "bout/utils.hxx" #include #include #include +#include +#include + +namespace { + +DIFF_METHOD parseConcreteDiffMethod(const std::string& method_name) { + const auto method = uppercase(method_name); + + if (method == "U1") { + return DIFF_U1; + } + if (method == "U2") { + return DIFF_U2; + } + if (method == "C2") { + return DIFF_C2; + } + if (method == "W2") { + return DIFF_W2; + } + if (method == "W3") { + return DIFF_W3; + } + if (method == "C4") { + return DIFF_C4; + } + if (method == "U3") { + return DIFF_U3; + } + if (method == "FFT") { + return DIFF_FFT; + } + if (method == "SPLIT") { + return DIFF_SPLIT; + } + if (method == "S2") { + return DIFF_S2; + } + if (method == "DEFAULT") { + throw BoutException("Default derivative options must resolve to a concrete method"); + } + + throw BoutException("Unknown differential method '{:s}'", method_name); +} + +} // namespace + /******************************************************************************* * Helper routines *******************************************************************************/ @@ -37,6 +85,52 @@ void Mesh::derivs_init(Options* options) { // of derivative. DerivativeStore::getInstance().initialise(options); DerivativeStore::getInstance().initialise(options); + + auto backup_section = options->getSection("diff"); + const std::array, 3> directions{{ + {DIRECTION::X, "ddx"}, + {DIRECTION::Y, "ddy"}, + {DIRECTION::Z, "ddz"}, + }}; + const std::array, 5> deriv_types{{ + {DERIV::Standard, "first"}, + {DERIV::StandardSecond, "second"}, + {DERIV::StandardFourth, "fourth"}, + {DERIV::Upwind, "upwind"}, + {DERIV::Flux, "flux"}, + }}; + + derivative_defaults = DerivativeDefaults{}; + + for (const auto& [direction, section_name] : directions) { + auto specific_section = options->getSection(section_name); + auto staggered_section = options->getSection(section_name + "stag"); + + for (const auto& [deriv, option_name] : deriv_types) { + auto default_method = DerivativeDefaults::builtinDefaultMethod(deriv); + auto default_name = toString(default_method); + + if (specific_section->isSet(option_name)) { + specific_section->get(option_name, default_name, default_name); + } else if (backup_section->isSet(option_name)) { + backup_section->get(option_name, default_name, default_name); + } + + default_method = parseConcreteDiffMethod(default_name); + derivative_defaults.set(direction, deriv, default_method, STAGGER::None); + + auto staggered_name = toString(default_method); + auto staggered_method = default_method; + if (staggered_section->isSet(option_name)) { + staggered_section->get(option_name, staggered_name, staggered_name); + staggered_method = parseConcreteDiffMethod(staggered_name); + } + + derivative_defaults.set(direction, deriv, staggered_method, STAGGER::C2L); + derivative_defaults.set(direction, deriv, staggered_method, STAGGER::L2C); + } + } + // Get the fraction of modes filtered out in FFT derivatives options->getSection("ddz")->get("fft_filter", fft_derivs_filter, 0.0); } diff --git a/tests/unit/include/bout/test_stencil_expr.cxx b/tests/unit/include/bout/test_stencil_expr.cxx index 898d9d5694..9d3154c29a 100644 --- a/tests/unit/include/bout/test_stencil_expr.cxx +++ b/tests/unit/include/bout/test_stencil_expr.cxx @@ -7,6 +7,7 @@ #include +#include "fake_mesh.hxx" #include "fake_mesh_fixture.hxx" namespace { @@ -73,6 +74,20 @@ TEST_F(DDZDispatchExprTest, UsesRequestedRegion) { TEST_F(DDZDispatchExprTest, RejectsUnsupportedMethods) { auto input = makeTestField(mesh_staggered, CELL_CENTRE); - EXPECT_THROW((void)DDZ_Dispatch(input, CELL_DEFAULT, DIFF_DEFAULT), BoutException); EXPECT_THROW((void)DDZ_Dispatch(input, CELL_DEFAULT, DIFF_FFT), BoutException); } + +TEST_F(DDZDispatchExprTest, ResolvesDefaultMethodFromMesh) { + auto input = makeTestField(mesh_staggered, CELL_CENTRE); + + Options diff_options{{"ddz", {{"first", "C4"}}}, {"ddzstag", {{"first", "C2"}}}}; + static_cast(mesh_staggered)->initDerivs(&diff_options); + + const auto centre_actual = Field3D{DDZ_Dispatch(input, CELL_CENTRE, DIFF_DEFAULT)}; + const auto centre_expected = DDZ(input, CELL_CENTRE, toString(DIFF_C4)); + EXPECT_TRUE(IsFieldEqual(centre_actual, centre_expected, "RGN_NOBNDRY")); + + const auto staggered_actual = Field3D{DDZ_Dispatch(input, CELL_ZLOW, DIFF_DEFAULT)}; + const auto staggered_expected = DDZ(input, CELL_ZLOW, toString(DIFF_C2)); + EXPECT_TRUE(IsFieldEqual(staggered_actual, staggered_expected, "RGN_NOBNDRY")); +} diff --git a/tests/unit/mesh/test_mesh.cxx b/tests/unit/mesh/test_mesh.cxx index 8d64f1d0a3..0b30176704 100644 --- a/tests/unit/mesh/test_mesh.cxx +++ b/tests/unit/mesh/test_mesh.cxx @@ -201,6 +201,44 @@ TEST_F(MeshTest, GetStringNoSource) { EXPECT_EQ(string_value, ""); } +TEST_F(MeshTest, GetDefaultMethodUsesBuiltinDefaults) { + Options options; + + localmesh.initDerivs(&options); + + EXPECT_EQ(localmesh.getDefaultMethod(DIRECTION::X, DERIV::Standard), DIFF_C2); + EXPECT_EQ(localmesh.getDefaultMethod(DIRECTION::Y, DERIV::StandardSecond), DIFF_C2); + EXPECT_EQ(localmesh.getDefaultMethod(DIRECTION::Z, DERIV::StandardFourth), DIFF_C2); + EXPECT_EQ(localmesh.getDefaultMethod(DIRECTION::Z, DERIV::Upwind), DIFF_U1); + EXPECT_EQ(localmesh.getDefaultMethod(DIRECTION::Z, DERIV::Flux), DIFF_U1); + EXPECT_EQ(localmesh.getDefaultMethod(DIRECTION::Z, DERIV::Standard, STAGGER::C2L), + DIFF_C2); + EXPECT_EQ(localmesh.getDefaultMethod(DIRECTION::Z, DERIV::Standard, STAGGER::L2C), + DIFF_C2); +} + +TEST_F(MeshTest, GetDefaultMethodReadsDirectionAndStaggeredOptions) { + Options options{{"diff", {{"first", "W2"}, {"upwind", "U2"}}}, + {"ddz", {{"first", "C4"}}}, + {"ddzstag", {{"first", "S2"}}}}; + + localmesh.initDerivs(&options); + + EXPECT_EQ(localmesh.getDefaultMethod(DIRECTION::X, DERIV::Standard), DIFF_W2); + EXPECT_EQ(localmesh.getDefaultMethod(DIRECTION::Y, DERIV::Upwind), DIFF_U2); + EXPECT_EQ(localmesh.getDefaultMethod(DIRECTION::Z, DERIV::Standard), DIFF_C4); + EXPECT_EQ(localmesh.getDefaultMethod(DIRECTION::Z, DERIV::Standard, STAGGER::C2L), + DIFF_S2); + EXPECT_EQ(localmesh.getDefaultMethod(DIRECTION::Z, DERIV::Standard, STAGGER::L2C), + DIFF_S2); +} + +TEST_F(MeshTest, GetDefaultMethodRejectsDefaultOptionValue) { + Options options{{"ddz", {{"first", "DEFAULT"}}}}; + + EXPECT_THROW(localmesh.initDerivs(&options), BoutException); +} + TEST_F(MeshTest, GetStringNoSourceWithDefault) { std::string string_value; const std::string default_value = "some default"; From 38fb3c29b5fad7faa0e169ef932d95bf4b1a3871 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 23 Jun 2026 15:27:49 -0700 Subject: [PATCH 07/18] Replace DDZ(Field3D) with BinaryExpr implementation Retains runtime choice of method, and handling of staggered inputs or outputs. Only C2 and C4 methods are supported. --- examples/hasegawa-wakatani/hw.cxx | 4 +- examples/performance/ddz/ddz.cxx | 12 +--- include/bout/derivs.hxx | 29 ++++----- include/bout/stencil_expr.hxx | 62 +++++++++++++++++-- .../laplace/impls/naulin/naulin_laplace.cxx | 4 +- src/sys/derivs.cxx | 46 +++++++++++--- tests/unit/include/bout/test_stencil_expr.cxx | 16 ++--- 7 files changed, 120 insertions(+), 53 deletions(-) diff --git a/examples/hasegawa-wakatani/hw.cxx b/examples/hasegawa-wakatani/hw.cxx index f4e546908f..45864fa7eb 100644 --- a/examples/hasegawa-wakatani/hw.cxx +++ b/examples/hasegawa-wakatani/hw.cxx @@ -110,8 +110,8 @@ class HW : public PhysicsModel { nonzonal_phi -= averageY(DC(phi)); } - ddt(n) = -bracket_arakawa(phi, n) + alpha * (nonzonal_phi - nonzonal_n) - - kappa * DDZ_Dispatch(phi); + ddt(n) = + -bracket_arakawa(phi, n) + alpha * (nonzonal_phi - nonzonal_n) - kappa * DDZ(phi); ddt(vort) = -bracket_arakawa(phi, vort) + alpha * (nonzonal_phi - nonzonal_n); diff --git a/examples/performance/ddz/ddz.cxx b/examples/performance/ddz/ddz.cxx index 5d3bf2c961..0b1be902f2 100644 --- a/examples/performance/ddz/ddz.cxx +++ b/examples/performance/ddz/ddz.cxx @@ -65,17 +65,9 @@ int main(int argc, char** argv) { // Nested loops over block data ITERATOR_TEST_BLOCK("DDZ Default", result = DDZ(a);); - ITERATOR_TEST_BLOCK("DDZ C2", result = DDZ(a, CELL_DEFAULT, "DIFF_C2");); + ITERATOR_TEST_BLOCK("DDZ C2", result = DDZ(a, CELL_DEFAULT, DIFF_C2);); - ITERATOR_TEST_BLOCK("DDZ C4", result = DDZ(a, CELL_DEFAULT, "DIFF_C4");); - - ITERATOR_TEST_BLOCK("DDZ S2", result = DDZ(a, CELL_DEFAULT, "DIFF_S2");); - - ITERATOR_TEST_BLOCK("DDZ W2", result = DDZ(a, CELL_DEFAULT, "DIFF_W2");); - - ITERATOR_TEST_BLOCK("DDZ W3", result = DDZ(a, CELL_DEFAULT, "DIFF_W3");); - - ITERATOR_TEST_BLOCK("DDZ FFT", result = DDZ(a, CELL_DEFAULT, "DIFF_FFT");); + ITERATOR_TEST_BLOCK("DDZ C4", result = DDZ(a, CELL_DEFAULT, DIFF_C4);); if (profileMode) { int nthreads = 0; diff --git a/include/bout/derivs.hxx b/include/bout/derivs.hxx index a8d9279378..5adb58b544 100644 --- a/include/bout/derivs.hxx +++ b/include/bout/derivs.hxx @@ -31,6 +31,7 @@ #include "bout/field2d.hxx" #include "bout/field3d.hxx" +#include "bout/stencil_expr.hxx" #include "bout/vector2d.hxx" #include "bout/vector3d.hxx" @@ -102,22 +103,6 @@ Coordinates::FieldMetric DDY(const Field2D& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", const std::string& region = "RGN_NOBNDRY"); -/// Calculate first partial derivative in Z -/// -/// \f$\partial / \partial z\f$ -/// -/// @param[in] f The field to be differentiated -/// @param[in] outloc The cell location where the result is desired. If -/// staggered grids is not enabled then this has no effect -/// If not given, defaults to CELL_DEFAULT -/// @param[in] method Differencing method to use. This overrides the default -/// If not given, defaults to DIFF_DEFAULT -/// @param[in] region What region is expected to be calculated -/// If not given, defaults to RGN_NOBNDRY -Field3D DDZ(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT, - const std::string& method = "DEFAULT", - const std::string& region = "RGN_NOBNDRY"); - /// Calculate first partial derivative in Z /// /// \f$\partial / \partial z\f$ @@ -147,7 +132,11 @@ Coordinates::FieldMetric DDZ(const Field2D& f, CELL_LOC outloc = CELL_DEFAULT, /// @param[in] region What region is expected to be calculated /// If not given, defaults to RGN_NOBNDRY Vector3D DDZ(const Vector3D& f, CELL_LOC outloc = CELL_DEFAULT, - const std::string& method = "DEFAULT", + DIFF_METHOD method = DIFF_DEFAULT, + const std::string& region = "RGN_NOBNDRY"); + +/// Compatibility overload for string-based callers. +Vector3D DDZ(const Vector3D& f, CELL_LOC outloc, const std::string& method, const std::string& region = "RGN_NOBNDRY"); /// Calculate first partial derivative in Z @@ -163,7 +152,11 @@ Vector3D DDZ(const Vector3D& f, CELL_LOC outloc = CELL_DEFAULT, /// @param[in] region What region is expected to be calculated /// If not given, defaults to RGN_NOBNDRY Vector2D DDZ(const Vector2D& f, CELL_LOC outloc = CELL_DEFAULT, - const std::string& method = "DEFAULT", + DIFF_METHOD method = DIFF_DEFAULT, + const std::string& region = "RGN_NOBNDRY"); + +/// Compatibility overload for string-based callers. +Vector2D DDZ(const Vector2D& f, CELL_LOC outloc, const std::string& method, const std::string& region = "RGN_NOBNDRY"); ////////// SECOND DERIVATIVES ////////// diff --git a/include/bout/stencil_expr.hxx b/include/bout/stencil_expr.hxx index bbccc267bf..8542fc06dc 100644 --- a/include/bout/stencil_expr.hxx +++ b/include/bout/stencil_expr.hxx @@ -11,6 +11,7 @@ #include "bout/fieldops.hxx" #include "bout/mesh.hxx" #include "bout/single_index_ops.hxx" +#include "bout/utils.hxx" #include @@ -194,6 +195,53 @@ using BracketArakawaExpr = BinaryExprgetRegion("RGN_NOBNDRY")}; } -inline bout::stencil::DDZDispatchExpr -DDZ_Dispatch(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT, - DIFF_METHOD method = DIFF_DEFAULT, - const std::string& region = "RGN_NOBNDRY") { +inline bout::stencil::DDZDispatchExpr DDZ(const Field3D& f, + CELL_LOC outloc = CELL_DEFAULT, + DIFF_METHOD method = DIFF_DEFAULT, + const std::string& region = "RGN_NOBNDRY") { checkData(f); const auto resolved_outloc = (outloc == CELL_DEFAULT) ? f.getLocation() : outloc; @@ -260,6 +308,12 @@ DDZ_Dispatch(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT, f.getMesh()->getRegion(region)}; } +inline bout::stencil::DDZDispatchExpr DDZ(const Field3D& f, CELL_LOC outloc, + const std::string& method, + const std::string& region = "RGN_NOBNDRY") { + return DDZ(f, outloc, parseDDZMethodString(method), region); +} + inline bout::stencil::BracketArakawaExpr bracket_arakawa(const Field3D& f, const Field3D& g) { checkData(f); diff --git a/src/invert/laplace/impls/naulin/naulin_laplace.cxx b/src/invert/laplace/impls/naulin/naulin_laplace.cxx index ae8e78d1ff..a323bc6c08 100644 --- a/src/invert/laplace/impls/naulin/naulin_laplace.cxx +++ b/src/invert/laplace/impls/naulin/naulin_laplace.cxx @@ -247,7 +247,7 @@ Field3D LaplaceNaulin::solve(const Field3D& rhs, const Field3D& x0) { Field3D coef_y = DDY(C2coef, location, "C2") / C1TimesD; // z-component of 1./(C1*D) * Grad_perp(C2) - Field3D coef_z = DDZ(C2coef, location, "FFT") / C1TimesD; + Field3D coef_z = DDZ(C2coef, location, DIFF_C4) / C1TimesD; Field3D AOverD = Acoef / Dcoef; @@ -286,7 +286,7 @@ Field3D LaplaceNaulin::solve(const Field3D& rhs, const Field3D& x0) { auto calc_b_guess = [&](const Field3D& x_in) { // Derivatives of x Field3D ddx_x = DDX(x_in, location, "C2"); - Field3D ddz_x = DDZ(x_in, location, "FFT"); + Field3D ddz_x = DDZ(x_in, location, DIFF_C4); return rhsOverD - (coords->g11 * coef_x_AC * ddx_x + coords->g33 * coef_z * ddz_x + coords->g13 * (coef_x_AC * ddz_x + coef_z * ddx_x)) diff --git a/src/sys/derivs.cxx b/src/sys/derivs.cxx index e449dbcd30..b0a5b14388 100644 --- a/src/sys/derivs.cxx +++ b/src/sys/derivs.cxx @@ -86,12 +86,6 @@ Coordinates::FieldMetric DDY(const Field2D& f, CELL_LOC outloc, const std::strin ////////////// Z DERIVATIVE ///////////////// -Field3D DDZ(const Field3D& f, CELL_LOC outloc, const std::string& method, - const std::string& region) { - return bout::derivatives::index::DDZ(f, outloc, method, region) - / f.getCoordinates(outloc)->dz; -} - Coordinates::FieldMetric DDZ(const Field2D& f, CELL_LOC UNUSED(outloc), const std::string& UNUSED(method), const std::string& UNUSED(region)) { @@ -100,7 +94,7 @@ Coordinates::FieldMetric DDZ(const Field2D& f, CELL_LOC UNUSED(outloc), return tmp; } -Vector3D DDZ(const Vector3D& v, CELL_LOC outloc, const std::string& method, +Vector3D DDZ(const Vector3D& v, CELL_LOC outloc, DIFF_METHOD method, const std::string& region) { Vector3D result(v.getMesh()); const Coordinates* metric = v.x.getCoordinates(outloc); @@ -131,8 +125,37 @@ Vector3D DDZ(const Vector3D& v, CELL_LOC outloc, const std::string& method, return result; } -Vector2D DDZ(const Vector2D& v, CELL_LOC UNUSED(outloc), - const std::string& UNUSED(method), const std::string& UNUSED(region)) { +Vector3D DDZ(const Vector3D& v, CELL_LOC outloc, const std::string& method, + const std::string& region) { + Vector3D result(v.getMesh()); + const Coordinates* metric = v.x.getCoordinates(outloc); + + if (v.covariant) { + result.x = DDZ(v.x, outloc, method, region) - v.x * metric->G1_13 + - v.y * metric->G2_13 - v.z * metric->G3_13; + result.y = DDZ(v.y, outloc, method, region) - v.x * metric->G1_23 + - v.y * metric->G2_23 - v.z * metric->G3_23; + result.z = DDZ(v.z, outloc, method, region) - v.x * metric->G1_33 + - v.y * metric->G2_33 - v.z * metric->G3_33; + result.covariant = true; + } else { + result.x = DDZ(v.x, outloc, method, region) + v.x * metric->G1_13 + + v.y * metric->G1_23 + v.z * metric->G1_33; + result.y = DDZ(v.y, outloc, method, region) + v.x * metric->G2_13 + + v.y * metric->G2_23 + v.z * metric->G2_33; + result.z = DDZ(v.z, outloc, method, region) + v.x * metric->G3_13 + + v.y * metric->G3_23 + v.z * metric->G3_33; + result.covariant = false; + } + + ASSERT2(((outloc == CELL_DEFAULT) && (result.getLocation() == v.getLocation())) + || (result.getLocation() == outloc)); + + return result; +} + +Vector2D DDZ(const Vector2D& v, CELL_LOC UNUSED(outloc), DIFF_METHOD UNUSED(method), + const std::string& UNUSED(region)) { Vector2D result(v.getMesh()); result.covariant = v.covariant; @@ -147,6 +170,11 @@ Vector2D DDZ(const Vector2D& v, CELL_LOC UNUSED(outloc), return result; } +Vector2D DDZ(const Vector2D& v, CELL_LOC outloc, const std::string& method, + const std::string& region) { + return DDZ(v, outloc, parseDDZMethodString(method), region); +} + /******************************************************************************* * 2nd derivative *******************************************************************************/ diff --git a/tests/unit/include/bout/test_stencil_expr.cxx b/tests/unit/include/bout/test_stencil_expr.cxx index 9d3154c29a..e2e36e3b00 100644 --- a/tests/unit/include/bout/test_stencil_expr.cxx +++ b/tests/unit/include/bout/test_stencil_expr.cxx @@ -54,7 +54,7 @@ TEST_P(DDZDispatchExprParamTest, MatchesDDZ) { const auto [inloc, outloc, method] = GetParam(); auto input = makeTestField(mesh_staggered, inloc); - const auto actual = Field3D{DDZ_Dispatch(input, outloc, method)}; + const auto actual = Field3D{DDZ(input, outloc, method)}; const auto expected = DDZ(input, outloc, toString(method)); EXPECT_EQ(actual.getLocation(), outloc); @@ -64,8 +64,8 @@ TEST_P(DDZDispatchExprParamTest, MatchesDDZ) { TEST_F(DDZDispatchExprTest, UsesRequestedRegion) { auto input = makeTestField(mesh_staggered, CELL_CENTRE); - const auto actual = Field3D{DDZ_Dispatch(input, CELL_ZLOW, DIFF_C2, "RGN_ALL")}; - const auto expected = DDZ(input, CELL_ZLOW, toString(DIFF_C2), "RGN_ALL"); + const auto actual = Field3D{DDZ(input, CELL_ZLOW, DIFF_C2, "RGN_ALL")}; + const auto expected = DDZ(input, CELL_ZLOW, "C2", "RGN_ALL"); EXPECT_EQ(actual.getLocation(), CELL_ZLOW); EXPECT_TRUE(IsFieldEqual(actual, expected, "RGN_ALL")); @@ -74,7 +74,7 @@ TEST_F(DDZDispatchExprTest, UsesRequestedRegion) { TEST_F(DDZDispatchExprTest, RejectsUnsupportedMethods) { auto input = makeTestField(mesh_staggered, CELL_CENTRE); - EXPECT_THROW((void)DDZ_Dispatch(input, CELL_DEFAULT, DIFF_FFT), BoutException); + EXPECT_THROW((void)DDZ(input, CELL_DEFAULT, DIFF_FFT), BoutException); } TEST_F(DDZDispatchExprTest, ResolvesDefaultMethodFromMesh) { @@ -83,11 +83,11 @@ TEST_F(DDZDispatchExprTest, ResolvesDefaultMethodFromMesh) { Options diff_options{{"ddz", {{"first", "C4"}}}, {"ddzstag", {{"first", "C2"}}}}; static_cast(mesh_staggered)->initDerivs(&diff_options); - const auto centre_actual = Field3D{DDZ_Dispatch(input, CELL_CENTRE, DIFF_DEFAULT)}; - const auto centre_expected = DDZ(input, CELL_CENTRE, toString(DIFF_C4)); + const auto centre_actual = Field3D{DDZ(input, CELL_CENTRE, DIFF_DEFAULT)}; + const auto centre_expected = DDZ(input, CELL_CENTRE, "C4"); EXPECT_TRUE(IsFieldEqual(centre_actual, centre_expected, "RGN_NOBNDRY")); - const auto staggered_actual = Field3D{DDZ_Dispatch(input, CELL_ZLOW, DIFF_DEFAULT)}; - const auto staggered_expected = DDZ(input, CELL_ZLOW, toString(DIFF_C2)); + const auto staggered_actual = Field3D{DDZ(input, CELL_ZLOW, DIFF_DEFAULT)}; + const auto staggered_expected = DDZ(input, CELL_ZLOW, "C2"); EXPECT_TRUE(IsFieldEqual(staggered_actual, staggered_expected, "RGN_NOBNDRY")); } From ccf8f692bb6366362d6ff0727e01d65800ad4228 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 24 Jun 2026 15:41:00 -0700 Subject: [PATCH 08/18] DDZ: Correct exception message --- include/bout/stencil_expr.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/bout/stencil_expr.hxx b/include/bout/stencil_expr.hxx index 8542fc06dc..9a02faf36b 100644 --- a/include/bout/stencil_expr.hxx +++ b/include/bout/stencil_expr.hxx @@ -289,7 +289,7 @@ inline bout::stencil::DDZDispatchExpr DDZ(const Field3D& f, } if ((method != DIFF_C2) && (method != DIFF_C4)) { - throw BoutException("DDZ_Dispatch only supports DIFF_C2 and DIFF_C4, got {:s}", + throw BoutException("DDZ only supports DIFF_C2 and DIFF_C4, got {:s}", toString(method)); } From 8c85aa0252dab7406bbde0dff56758eab08877c1 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 29 Aug 2026 10:55:02 -0700 Subject: [PATCH 09/18] Derivative default resets Fixes failing serial tests --- include/bout/deriv_store.hxx | 12 ++++++++++++ include/bout/fieldperp.hxx | 5 ----- src/field/fieldperp.cxx | 35 ----------------------------------- src/sys/options.cxx | 11 ++++++++++- tests/unit/mesh/test_mesh.cxx | 1 + 5 files changed, 23 insertions(+), 41 deletions(-) diff --git a/include/bout/deriv_store.hxx b/include/bout/deriv_store.hxx index 7c08802cae..0a995fd7a3 100644 --- a/include/bout/deriv_store.hxx +++ b/include/bout/deriv_store.hxx @@ -308,6 +308,8 @@ struct DerivativeStore { }; void initialise(Options* options) { + defaultMethods.clear(); + setDefaults(); // To replicate the existing behaviour we first search for a section called //"dd?" and if the option isn't in there we search a section called "diff" @@ -355,6 +357,11 @@ struct DerivativeStore { backupSection->get(derivName, theDefault, ""); } + if (uppercase(theDefault) == toString(DIFF_DEFAULT)) { + throw BoutException("Default derivative options must resolve to a concrete " + "method"); + } + // Now we have the default method we should store it in defaultMethods theDefault = uppercase(theDefault); defaultMethods[getKey(theDirection, STAGGER::None, theDerivTypeString)] = @@ -377,6 +384,11 @@ struct DerivativeStore { specificSection->get(derivName, theDefault, ""); } + if (uppercase(theDefault) == toString(DIFF_DEFAULT)) { + throw BoutException("Default derivative options must resolve to a concrete " + "method"); + } + // Now we have the default method we should store it in defaultMethods theDefault = uppercase(theDefault); defaultMethods[getKey(theDirection, STAGGER::L2C, theDerivTypeString)] = diff --git a/include/bout/fieldperp.hxx b/include/bout/fieldperp.hxx index 2f8408d82b..36b79a3a85 100644 --- a/include/bout/fieldperp.hxx +++ b/include/bout/fieldperp.hxx @@ -409,11 +409,6 @@ 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/src/field/fieldperp.cxx b/src/field/fieldperp.cxx index bca551c29d..b7b2d9d731 100644 --- a/src/field/fieldperp.cxx +++ b/src/field/fieldperp.cxx @@ -153,41 +153,6 @@ 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/sys/options.cxx b/src/sys/options.cxx index eee8018f5d..9580a20fa0 100644 --- a/src/sys/options.cxx +++ b/src/sys/options.cxx @@ -3,6 +3,7 @@ #include "bout/array.hxx" #include "bout/bout_types.hxx" #include "bout/boutexception.hxx" +#include "bout/deriv_store.hxx" #include "bout/field2d.hxx" #include "bout/field3d.hxx" #include "bout/field_factory.hxx" // Used for parsing expressions @@ -45,7 +46,15 @@ Options& Options::root() { return root_instance; } -void Options::cleanup() { root() = Options{}; } +void Options::cleanup() { + root() = Options{}; + + // Derivative defaults are configured from options and stored in global singletons. + // Restore their option-derived defaults without clearing the registered kernels. + Options defaults; + DerivativeStore::getInstance().initialise(&defaults); + DerivativeStore::getInstance().initialise(&defaults); +} Options Options::copy() const { Options result; diff --git a/tests/unit/mesh/test_mesh.cxx b/tests/unit/mesh/test_mesh.cxx index 0b30176704..d8ae034a73 100644 --- a/tests/unit/mesh/test_mesh.cxx +++ b/tests/unit/mesh/test_mesh.cxx @@ -14,6 +14,7 @@ class MeshTest : public ::testing::Test { public: MeshTest() : localmesh(nx, ny, nz) {} + ~MeshTest() override { Options::cleanup(); } static const int nx = 3; static const int ny = 5; static const int nz = 7; From d05c2274cf37def29d7d2d9692bc06b52faf11bc Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sun, 30 Aug 2026 11:29:37 -0700 Subject: [PATCH 10/18] Interchange instability: Use C4 DDZ method --- tests/integrated/test-interchange-instability/data_1/BOUT.inp | 2 +- tests/integrated/test-interchange-instability/data_10/BOUT.inp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/integrated/test-interchange-instability/data_1/BOUT.inp b/tests/integrated/test-interchange-instability/data_1/BOUT.inp index d664def68f..ecc022bdd0 100644 --- a/tests/integrated/test-interchange-instability/data_1/BOUT.inp +++ b/tests/integrated/test-interchange-instability/data_1/BOUT.inp @@ -30,7 +30,7 @@ upwind = W3 [mesh:ddz] -first = FFT +first = C4 second = FFT upwind = W3 diff --git a/tests/integrated/test-interchange-instability/data_10/BOUT.inp b/tests/integrated/test-interchange-instability/data_10/BOUT.inp index acbc66c839..6a9799e182 100644 --- a/tests/integrated/test-interchange-instability/data_10/BOUT.inp +++ b/tests/integrated/test-interchange-instability/data_10/BOUT.inp @@ -30,7 +30,7 @@ upwind = W3 [mesh:ddz] -first = FFT +first = C4 second = FFT upwind = W3 From 85648fdde17f7e9d93c4121fac8e2c075d0bad72 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sun, 30 Aug 2026 20:40:30 -0700 Subject: [PATCH 11/18] DDX implement lazy version --- include/bout/coordinates.hxx | 6 + include/bout/coordinates_accessor.hxx | 4 +- include/bout/derivs.hxx | 4 - include/bout/stencil_expr.hxx | 424 +++++++++++++++++- src/mesh/coordinates.cxx | 2 + src/mesh/coordinates_accessor.cxx | 18 +- src/mesh/mesh.cxx | 1 + src/sys/derivs.cxx | 18 +- tests/unit/fake_mesh.hxx | 3 + tests/unit/include/bout/test_stencil_expr.cxx | 156 +++++++ 10 files changed, 610 insertions(+), 26 deletions(-) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index a090dcd959..c6641a0664 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -80,6 +80,12 @@ public: FieldMetric g_22, FieldMetric g_33, FieldMetric g_12, FieldMetric g_13, FieldMetric g_23, FieldMetric ShiftTorsion, FieldMetric IntShiftTorsion); + ~Coordinates(); + Coordinates(const Coordinates&) = delete; + Coordinates& operator=(const Coordinates&) = delete; + Coordinates(Coordinates&&) noexcept = default; + Coordinates& operator=(Coordinates&&) noexcept = default; + /// Add variables to \p output_options, for post-processing void outputVars(Options& output_options); diff --git a/include/bout/coordinates_accessor.hxx b/include/bout/coordinates_accessor.hxx index 2376ab5039..a117718274 100644 --- a/include/bout/coordinates_accessor.hxx +++ b/include/bout/coordinates_accessor.hxx @@ -54,7 +54,8 @@ struct CoordinatesAccessor { d1_dx, d1_dy, d1_dz, // Grid spacing non-uniformity - J, // Jacobian + IntShiftTorsion, + J, // Jacobian B, Byup, Bydown, // Magnetic field magnitude @@ -119,6 +120,7 @@ struct CoordinatesAccessor { COORD_FN(dx, dy, dz); COORD_FN(d1_dx, d1_dy, d1_dz); + COORD_FN(IntShiftTorsion); COORD_FN(J); COORD_FN(B, Byup, Bydown); COORD_FN(G1, G3); diff --git a/include/bout/derivs.hxx b/include/bout/derivs.hxx index b004d4a8fa..4684a95a3e 100644 --- a/include/bout/derivs.hxx +++ b/include/bout/derivs.hxx @@ -54,10 +54,6 @@ /// If not given, defaults to DIFF_DEFAULT /// @param[in] region What region is expected to be calculated /// If not given, defaults to RGN_NOBNDRY -Field3D DDX(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT, - const std::string& method = "DEFAULT", - const std::string& region = "RGN_NOBNDRY"); - /// Calculate first partial derivative in X /// /// \f$\partial / \partial x\f$ diff --git a/include/bout/stencil_expr.hxx b/include/bout/stencil_expr.hxx index 9a02faf36b..36f8d87c72 100644 --- a/include/bout/stencil_expr.hxx +++ b/include/bout/stencil_expr.hxx @@ -9,14 +9,313 @@ #include "bout/field.hxx" #include "bout/field3d.hxx" #include "bout/fieldops.hxx" +#include "bout/index_derivs.hxx" #include "bout/mesh.hxx" #include "bout/single_index_ops.hxx" #include "bout/utils.hxx" +#include #include namespace bout::stencil { +struct DDX_C2_Op { + CoordinatesAccessor coords; + int ny{0}; + int nz{0}; + bool inc_int_shear{false}; + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_ddz_c2(int idx, + const LView& lhs) const { + const int izp = i_zp(idx, nz); + const int izm = i_zm(idx, nz); + + return 0.5 * (lhs(izp) - lhs(izm)) / coords.dz(idx); + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal operator()(int idx, const LView& lhs, + const RView&) const { + const int ixp = i_xp(idx, ny, nz); + const int ixm = i_xm(idx, ny, nz); + + BoutReal result = 0.5 * (lhs(ixp) - lhs(ixm)) / coords.dx(idx); + if (inc_int_shear) { + result += coords.IntShiftTorsion(idx) * apply_ddz_c2(idx, lhs); + } + return result; + } +}; + +struct DDX_C4_Op { + CoordinatesAccessor coords; + int ny{0}; + int nz{0}; + bool inc_int_shear{false}; + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_ddz_c4(int idx, + const LView& lhs) const { + const int izp = i_zp(idx, nz); + const int izm = i_zm(idx, nz); + const int izp2 = i_zp(izp, nz); + const int izm2 = i_zm(izm, nz); + + return (-lhs(izp2) + 8.0 * lhs(izp) - 8.0 * lhs(izm) + lhs(izm2)) + / (12.0 * coords.dz(idx)); + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal operator()(int idx, const LView& lhs, + const RView&) const { + const int ixp = i_xp(idx, ny, nz); + const int ixm = i_xm(idx, ny, nz); + const int ixp2 = i_xp(ixp, ny, nz); + const int ixm2 = i_xm(ixm, ny, nz); + + BoutReal result = (-lhs(ixp2) + 8.0 * lhs(ixp) - 8.0 * lhs(ixm) + lhs(ixm2)) + / (12.0 * coords.dx(idx)); + if (inc_int_shear) { + result += coords.IntShiftTorsion(idx) * apply_ddz_c4(idx, lhs); + } + return result; + } +}; + +struct DDX_Dispatch_Op { + CoordinatesAccessor coords; + int ny{0}; + int nz{0}; + STAGGER stagger{STAGGER::None}; + DIFF_METHOD method{DIFF_DEFAULT}; + bool inc_int_shear{false}; + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_ddz_c2(int idx, + const LView& lhs) const { + const int izp = i_zp(idx, nz); + const int izm = i_zm(idx, nz); + + return 0.5 * (lhs(izp) - lhs(izm)) / coords.dz(idx); + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_ddz_c4(int idx, + const LView& lhs) const { + const int izp = i_zp(idx, nz); + const int izm = i_zm(idx, nz); + const int izp2 = i_zp(izp, nz); + const int izm2 = i_zm(izm, nz); + + return (-lhs(izp2) + 8.0 * lhs(izp) - 8.0 * lhs(izm) + lhs(izm2)) + / (12.0 * coords.dz(idx)); + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c2_none(int idx, + const LView& lhs) const { + const int ixp = i_xp(idx, ny, nz); + const int ixm = i_xm(idx, ny, nz); + + return 0.5 * (lhs(ixp) - lhs(ixm)) / coords.dx(idx); + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c4_none(int idx, + const LView& lhs) const { + const int ixp = i_xp(idx, ny, nz); + const int ixm = i_xm(idx, ny, nz); + const int ixp2 = i_xp(ixp, ny, nz); + const int ixm2 = i_xm(ixm, ny, nz); + + return (-lhs(ixp2) + 8.0 * lhs(ixp) - 8.0 * lhs(ixm) + lhs(ixm2)) + / (12.0 * coords.dx(idx)); + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_w2_none(int idx, + const LView& lhs) const { + const int ixp = i_xp(idx, ny, nz); + const int ixm = i_xm(idx, ny, nz); + + const BoutReal dc = 0.5 * (lhs(ixp) - lhs(ixm)); + const BoutReal dl = lhs(idx) - lhs(ixm); + const BoutReal dr = lhs(ixp) - lhs(idx); + + const BoutReal isl = SQ(dl); + const BoutReal isr = SQ(dr); + const BoutReal isc = (13.0 / 3.0) * SQ(lhs(ixp) - 2.0 * lhs(idx) + lhs(ixm)) + + 0.25 * SQ(lhs(ixp) - lhs(ixm)); + + const BoutReal al = 0.25 / SQ(WENO_SMALL + isl); + const BoutReal ar = 0.25 / SQ(WENO_SMALL + isr); + const BoutReal ac = 0.5 / SQ(WENO_SMALL + isc); + const BoutReal sa = al + ar + ac; + + return (al * dl + ar * dr + ac * dc) / (sa * coords.dx(idx)); + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_w3_none(int idx, + const LView& lhs) const { + const int ixp = i_xp(idx, ny, nz); + const int ixm = i_xm(idx, ny, nz); + const int ixp2 = i_xp(ixp, ny, nz); + const int ixm2 = i_xm(ixm, ny, nz); + + BoutReal ma = fabs(lhs(idx)); + ma = BOUTMAX(ma, fabs(lhs(ixm))); + ma = BOUTMAX(ma, fabs(lhs(ixp))); + ma = BOUTMAX(ma, fabs(lhs(ixm2))); + ma = BOUTMAX(ma, fabs(lhs(ixp2))); + + const BoutReal sp_mm = lhs(ixm2) + ma; + const BoutReal sp_m = lhs(ixm) + ma; + const BoutReal sp_c = lhs(idx) + ma; + const BoutReal sp_p = lhs(ixp) + ma; + const BoutReal sm_m = ma - lhs(ixm); + const BoutReal sm_c = ma - lhs(idx); + const BoutReal sm_p = ma - lhs(ixp); + const BoutReal sm_pp = ma - lhs(ixp2); + + BoutReal r = (WENO_SMALL + SQ(sp_c - 2.0 * sp_m + sp_mm)) + / (WENO_SMALL + SQ(sp_p - 2.0 * sp_c + sp_m)); + BoutReal deriv = -sp_mm + 3.0 * sp_m - 3.0 * sp_c + sp_p; + const BoutReal w_pos = 1.0 / (1.0 + 2.0 * r * r); + const BoutReal pos = 0.25 * ((sp_p - sp_m) - w_pos * deriv); + + r = (WENO_SMALL + SQ(sm_pp - 2.0 * sm_p + sm_c)) + / (WENO_SMALL + SQ(sm_p - 2.0 * sm_c + sm_m)); + deriv = -sm_m + 3.0 * sm_c - 3.0 * sm_p + sm_pp; + const BoutReal w_neg = 1.0 / (1.0 + 2.0 * r * r); + const BoutReal neg = -0.25 * ((sm_p - sm_m) - w_neg * deriv); + + return (pos + neg) / coords.dx(idx); + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_s2_none(int idx, + const LView& lhs) const { + const int ixp = i_xp(idx, ny, nz); + const int ixm = i_xm(idx, ny, nz); + const int ixp2 = i_xp(ixp, ny, nz); + const int ixm2 = i_xm(ixm, ny, nz); + + BoutReal result = (-lhs(ixp2) + 8.0 * lhs(ixp) - 8.0 * lhs(ixm) + lhs(ixm2)) / 12.0; + result += SIGN(lhs(idx)) + * (lhs(ixp2) - 4.0 * lhs(ixp) + 6.0 * lhs(idx) - 4.0 * lhs(ixm) + lhs(ixm2)) + / 12.0; + + return result / coords.dx(idx); + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c2_c2l(int idx, + const LView& lhs) const { + const int ixm = i_xm(idx, ny, nz); + + return (lhs(idx) - lhs(ixm)) / coords.dx(idx); + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c4_c2l(int idx, + const LView& lhs) const { + const int ixm = i_xm(idx, ny, nz); + const int ixp = i_xp(idx, ny, nz); + const int ixm2 = i_xm(ixm, ny, nz); + + return (27.0 * (lhs(idx) - lhs(ixm)) - (lhs(ixp) - lhs(ixm2))) + / (24.0 * coords.dx(idx)); + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c2_l2c(int idx, + const LView& lhs) const { + const int ixp = i_xp(idx, ny, nz); + + return (lhs(ixp) - lhs(idx)) / coords.dx(idx); + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c4_l2c(int idx, + const LView& lhs) const { + const int ixp = i_xp(idx, ny, nz); + const int ixm = i_xm(idx, ny, nz); + const int ixp2 = i_xp(ixp, ny, nz); + + return (27.0 * (lhs(ixp) - lhs(idx)) - (lhs(ixp2) - lhs(ixm))) + / (24.0 * coords.dx(idx)); + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c2(int idx, const LView& lhs) const { + switch (stagger) { + case STAGGER::None: + return apply_c2_none(idx, lhs); + case STAGGER::C2L: + return apply_c2_c2l(idx, lhs); + case STAGGER::L2C: + return apply_c2_l2c(idx, lhs); + } + return 0.0; + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c4(int idx, const LView& lhs) const { + switch (stagger) { + case STAGGER::None: + return apply_c4_none(idx, lhs); + case STAGGER::C2L: + return apply_c4_c2l(idx, lhs); + case STAGGER::L2C: + return apply_c4_l2c(idx, lhs); + } + return 0.0; + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_ddx(int idx, const LView& lhs) const { + switch (method) { + case DIFF_C2: + return apply_c2(idx, lhs); + case DIFF_C4: + return apply_c4(idx, lhs); + case DIFF_W2: + return apply_w2_none(idx, lhs); + case DIFF_W3: + return apply_w3_none(idx, lhs); + case DIFF_S2: + return apply_s2_none(idx, lhs); + default: + return 0.0; + } + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_torsion(int idx, + const LView& lhs) const { + if (!inc_int_shear) { + return 0.0; + } + + switch (method) { + case DIFF_C2: + return coords.IntShiftTorsion(idx) * apply_ddz_c2(idx, lhs); + case DIFF_C4: + return coords.IntShiftTorsion(idx) * apply_ddz_c4(idx, lhs); + default: + return 0.0; + } + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal operator()(int idx, const LView& lhs, + const RView&) const { + return apply_ddx(idx, lhs) + apply_torsion(idx, lhs); + } +}; + struct DDZ_C2_Op { CoordinatesAccessor coords; int nz{0}; @@ -191,13 +490,41 @@ struct BracketArakawaOp { using DDZExprC2 = BinaryExpr; using DDZExprC4 = BinaryExpr; using DDZDispatchExpr = BinaryExpr; +using DDXExprC2 = BinaryExpr; +using DDXExprC4 = BinaryExpr; +using DDXDispatchExpr = BinaryExpr; using BracketArakawaExpr = BinaryExpr; } // namespace bout::stencil namespace { -inline DIFF_METHOD parseDDZMethodString(const std::string& method) { +inline int requiredDDXGuards(DIFF_METHOD method, STAGGER stagger) { + if (stagger != STAGGER::None) { + switch (method) { + case DIFF_C2: + return 1; + case DIFF_C4: + return 2; + default: + return -1; + } + } + + switch (method) { + case DIFF_C2: + case DIFF_W2: + return 1; + case DIFF_C4: + case DIFF_W3: + case DIFF_S2: + return 2; + default: + return -1; + } +} + +inline DIFF_METHOD parseField3DMethodString(const std::string& method) { auto normalized = uppercase(method); if (normalized.rfind("DIFF_", 0) == 0) { normalized = normalized.substr(5); @@ -237,11 +564,102 @@ inline DIFF_METHOD parseDDZMethodString(const std::string& method) { return DIFF_SPLIT; } - throw BoutException("Unknown DDZ method '{:s}'", method); + throw BoutException("Unknown field derivative method '{:s}'", method); } } // namespace +inline bout::stencil::DDXExprC2 DDX_C2(const Field3D& f) { + checkData(f); + + const auto region_id = f.getMesh()->getRegionID("RGN_NOBNDRY"); + + return bout::stencil::DDXExprC2{ + static_cast(f), + static_cast(f), + bout::stencil::DDX_C2_Op{CoordinatesAccessor{f.getCoordinates()}, f.getNy(), + f.getNz(), f.getMesh()->IncIntShear}, + f.getMesh(), + f.getLocation(), + f.getDirections(), + region_id, + f.getMesh()->getRegion("RGN_NOBNDRY")}; +} + +inline bout::stencil::DDXExprC4 DDX_C4(const Field3D& f) { + checkData(f); + + const auto region_id = f.getMesh()->getRegionID("RGN_NOBNDRY"); + + return bout::stencil::DDXExprC4{ + static_cast(f), + static_cast(f), + bout::stencil::DDX_C4_Op{CoordinatesAccessor{f.getCoordinates()}, f.getNy(), + f.getNz(), f.getMesh()->IncIntShear}, + f.getMesh(), + f.getLocation(), + f.getDirections(), + region_id, + f.getMesh()->getRegion("RGN_NOBNDRY")}; +} + +inline bout::stencil::DDXDispatchExpr DDX(const Field3D& f, + CELL_LOC outloc = CELL_DEFAULT, + DIFF_METHOD method = DIFF_DEFAULT, + const std::string& region = "RGN_NOBNDRY") { + checkData(f); + + const auto resolved_outloc = (outloc == CELL_DEFAULT) ? f.getLocation() : outloc; + const auto stagger = + f.getMesh()->getStagger(f.getLocation(), resolved_outloc, CELL_XLOW); + + if (method == DIFF_DEFAULT) { + method = f.getMesh()->getDefaultMethod(DIRECTION::X, DERIV::Standard, stagger); + } + + const bool supported_none = (method == DIFF_C2) || (method == DIFF_C4) + || (method == DIFF_W2) || (method == DIFF_W3) + || (method == DIFF_S2); + const bool supported_staggered = (method == DIFF_C2) || (method == DIFF_C4); + + if ((stagger == STAGGER::None) ? !supported_none : !supported_staggered) { + throw BoutException("DDX only supports DIFF_C2, DIFF_C4, DIFF_W2, DIFF_W3, " + "and DIFF_S2 for unstaggered grids, and DIFF_C2/DIFF_C4 " + "for staggered grids; got {:s}", + toString(method)); + } + + if (f.getMesh()->IncIntShear && (method != DIFF_C2) && (method != DIFF_C4)) { + throw BoutException("DDX with integrated shear only supports DIFF_C2 and " + "DIFF_C4, got {:s}", + toString(method)); + } + + const auto required_guards = requiredDDXGuards(method, stagger); + ASSERT2(required_guards >= 0); + ASSERT2(f.getMesh()->getNguard(DIRECTION::X) >= required_guards); + + const auto region_id = f.getMesh()->getRegionID(region); + + return bout::stencil::DDXDispatchExpr{ + static_cast(f), + static_cast(f), + bout::stencil::DDX_Dispatch_Op{ + CoordinatesAccessor{f.getCoordinates(resolved_outloc)}, f.getNy(), f.getNz(), + stagger, method, f.getMesh()->IncIntShear}, + f.getMesh(), + resolved_outloc, + f.getDirections(), + region_id, + f.getMesh()->getRegion(region)}; +} + +inline bout::stencil::DDXDispatchExpr DDX(const Field3D& f, CELL_LOC outloc, + const std::string& method, + const std::string& region = "RGN_NOBNDRY") { + return DDX(f, outloc, parseField3DMethodString(method), region); +} + inline bout::stencil::DDZExprC2 DDZ_C2(const Field3D& f) { checkData(f); @@ -311,7 +729,7 @@ inline bout::stencil::DDZDispatchExpr DDZ(const Field3D& f, inline bout::stencil::DDZDispatchExpr DDZ(const Field3D& f, CELL_LOC outloc, const std::string& method, const std::string& region = "RGN_NOBNDRY") { - return DDZ(f, outloc, parseDDZMethodString(method), region); + return DDZ(f, outloc, parseField3DMethodString(method), region); } inline bout::stencil::BracketArakawaExpr bracket_arakawa(const Field3D& f, diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index bad8ebbd84..3df74a5aae 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -358,6 +358,8 @@ Coordinates::Coordinates(Mesh* mesh, Options* options, const CELL_LOC loc, } } +Coordinates::~Coordinates() { CoordinatesAccessor::clear(this); } + void Coordinates::readFromMesh(Options* options, const std::string& suffix) { if (options == nullptr) { options = Options::getRoot()->getSection("mesh"); diff --git a/src/mesh/coordinates_accessor.cxx b/src/mesh/coordinates_accessor.cxx index a8a1f84f5d..272c67a22e 100644 --- a/src/mesh/coordinates_accessor.cxx +++ b/src/mesh/coordinates_accessor.cxx @@ -57,6 +57,7 @@ CoordinatesAccessor::CoordinatesAccessor(const Coordinates* coords) { for (const auto& ind : coords->dx().getRegion("RGN_ALL")) { COPY_STRIPE(dx, dy, dz); COPY_STRIPE(d1_dx, d1_dy, d1_dz); + COPY_STRIPE(IntShiftTorsion); COPY_STRIPE(J); if (coords->Bxy().isAllocated()) { @@ -70,11 +71,24 @@ CoordinatesAccessor::CoordinatesAccessor(const Coordinates* coords) { coords->Bxy().ydown()[ind]; } } - - COPY_STRIPE(G1, G3); COPY_STRIPE(g11, g12, g13, g22, g23, g33); COPY_STRIPE(g_11, g_12, g_13, g_22, g_23, g_33); } + + // G1/G3 may be computed from derivatives of the metric coefficients. Populate the + // base coordinate data first so any recursive accessor construction during that work + // reuses this seeded cache rather than trying to build a second incomplete copy. + const auto& G1 = coords->G1(); + const auto& G3 = coords->G3(); + + for (const auto& ind : coords->dx().getRegion("RGN_ALL")) { + if (G1.isAllocated()) { + COPY_STRIPE(G1); + } + if (G3.isAllocated()) { + COPY_STRIPE(G3); + } + } } std::size_t CoordinatesAccessor::clear(const Coordinates* coords) { diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx index de17645e30..c8a18e8768 100644 --- a/src/mesh/mesh.cxx +++ b/src/mesh/mesh.cxx @@ -586,6 +586,7 @@ std::shared_ptr Mesh::getCoordinatesSmart(CELL_LOC location) { inserted.first->second->communicateMetricTensor(); inserted.first->second->communicateDz(); + inserted.first->second->g_values(); return inserted.first->second; } diff --git a/src/sys/derivs.cxx b/src/sys/derivs.cxx index 5f88933fe0..fa61b73b71 100644 --- a/src/sys/derivs.cxx +++ b/src/sys/derivs.cxx @@ -60,20 +60,6 @@ ////////////// X DERIVATIVE ///////////////// -Field3D DDX(const Field3D& f, CELL_LOC outloc, const std::string& method, - const std::string& region) { - const auto& coords = *f.getCoordinates(outloc); - - Field3D result = bout::derivatives::index::DDX(f, outloc, method, region) / coords.dx(); - - if (f.getMesh()->IncIntShear) { - // Using BOUT-06 style shifting - result += coords.IntShiftTorsion() * DDZ(f, outloc, method, region); - } - - return result; -} - bout::FieldMetric DDX(const Field2D& f, CELL_LOC outloc, const std::string& method, const std::string& region) { ASSERT1(f.getLocation() == outloc || outloc == CELL_DEFAULT); @@ -144,7 +130,7 @@ Vector3D DDZ(const Vector3D& v, CELL_LOC outloc, DIFF_METHOD method, Vector3D DDZ(const Vector3D& v, CELL_LOC outloc, const std::string& method, const std::string& region) { - return DDZ(v, outloc, parseDDZMethodString(method), region); + return DDZ(v, outloc, parseField3DMethodString(method), region); } Vector2D DDZ(const Vector2D& v, CELL_LOC UNUSED(outloc), DIFF_METHOD UNUSED(method), @@ -165,7 +151,7 @@ Vector2D DDZ(const Vector2D& v, CELL_LOC UNUSED(outloc), DIFF_METHOD UNUSED(meth Vector2D DDZ(const Vector2D& v, CELL_LOC outloc, const std::string& method, const std::string& region) { - return DDZ(v, outloc, parseDDZMethodString(method), region); + return DDZ(v, outloc, parseField3DMethodString(method), region); } /******************************************************************************* diff --git a/tests/unit/fake_mesh.hxx b/tests/unit/fake_mesh.hxx index 216d04e18a..6a98c8528b 100644 --- a/tests/unit/fake_mesh.hxx +++ b/tests/unit/fake_mesh.hxx @@ -92,6 +92,9 @@ public: void setCoordinates(std::shared_ptr coords, CELL_LOC location = CELL_CENTRE) { coords_map[location] = std::move(coords); + if (coords_map[location] != nullptr) { + coords_map[location]->g_values(); + } } void setGridDataSource(GridDataSource* source_in) { source = source_in; } diff --git a/tests/unit/include/bout/test_stencil_expr.cxx b/tests/unit/include/bout/test_stencil_expr.cxx index e2e36e3b00..bf9bf67f6d 100644 --- a/tests/unit/include/bout/test_stencil_expr.cxx +++ b/tests/unit/include/bout/test_stencil_expr.cxx @@ -3,6 +3,7 @@ #include "test_extras.hxx" #include "bout/derivs.hxx" +#include "bout/index_derivs_interface.hxx" #include "bout/stencil_expr.hxx" #include @@ -23,6 +24,52 @@ Field3D makeTestField(Mesh* mesh, CELL_LOC location) { return result; } +Field3D expectedDDX(const Field3D& input, CELL_LOC outloc, const std::string& method, + const std::string& region = "RGN_NOBNDRY") { + const auto resolved_outloc = (outloc == CELL_DEFAULT) ? input.getLocation() : outloc; + const auto* coords = input.getCoordinates(resolved_outloc); + auto dx = Field3D{coords->dx()}.setLocation(resolved_outloc); + Field3D result = bout::derivatives::index::DDX(input, resolved_outloc, method, region); + result /= dx; + + if (input.getMesh()->IncIntShear) { + auto dz = Field3D{coords->dz()}.setLocation(resolved_outloc); + auto torsion = Field3D{coords->IntShiftTorsion()}.setLocation(resolved_outloc); + auto torsion_term = + bout::derivatives::index::DDZ(input, resolved_outloc, method, region); + torsion_term /= dz; + torsion_term *= torsion; + result += torsion_term; + } + + return result; +} + +class DDXDispatchExprTest : public FakeMeshFixture_tmpl<7, 5, 7> { +public: + DDXDispatchExprTest() { + for (auto* current_mesh : {bout::globals::mesh, mesh_staggered}) { + current_mesh->xstart = 2; + current_mesh->xend = current_mesh->LocalNx - 3; + current_mesh->addRegion3D(x_safe_region, Region(2, current_mesh->LocalNx - 3, + 0, current_mesh->LocalNy - 1, + 0, current_mesh->LocalNz - 1, + current_mesh->LocalNy, + current_mesh->LocalNz)); + } + } + + static constexpr auto x_safe_region = "RGN_XSAFE"; +}; + +class DDXDispatchExprParamTest + : public DDXDispatchExprTest, + public ::testing::WithParamInterface> { +}; + +class DDXDispatchExprMethodTest : public DDXDispatchExprTest, + public ::testing::WithParamInterface {}; + class DDZDispatchExprTest : public FakeMeshFixture {}; class DDZDispatchExprParamTest @@ -36,8 +83,27 @@ std::string paramToString( return toString(inloc) + "_to_" + toString(outloc) + "_" + toString(method); } +std::string methodToString(const ::testing::TestParamInfo& param) { + return toString(param.param); +} + } // namespace +INSTANTIATE_TEST_SUITE_P( + SupportedLocations, DDXDispatchExprParamTest, + ::testing::Values(std::make_tuple(CELL_CENTRE, CELL_CENTRE, DIFF_C2), + std::make_tuple(CELL_CENTRE, CELL_CENTRE, DIFF_C4), + std::make_tuple(CELL_CENTRE, CELL_XLOW, DIFF_C2), + std::make_tuple(CELL_CENTRE, CELL_XLOW, DIFF_C4), + std::make_tuple(CELL_XLOW, CELL_CENTRE, DIFF_C2), + std::make_tuple(CELL_XLOW, CELL_CENTRE, DIFF_C4), + std::make_tuple(CELL_XLOW, CELL_XLOW, DIFF_C2), + std::make_tuple(CELL_XLOW, CELL_XLOW, DIFF_C4)), + paramToString); + +INSTANTIATE_TEST_SUITE_P(SupportedMethods, DDXDispatchExprMethodTest, + ::testing::Values(DIFF_W2, DIFF_W3, DIFF_S2), methodToString); + INSTANTIATE_TEST_SUITE_P( SupportedLocations, DDZDispatchExprParamTest, ::testing::Values(std::make_tuple(CELL_CENTRE, CELL_CENTRE, DIFF_C2), @@ -50,6 +116,31 @@ INSTANTIATE_TEST_SUITE_P( std::make_tuple(CELL_ZLOW, CELL_ZLOW, DIFF_C4)), paramToString); +TEST_P(DDXDispatchExprParamTest, MatchesDDX) { + const auto [inloc, outloc, method] = GetParam(); + auto input = makeTestField(mesh_staggered, inloc); + + const auto actual = + Field3D{DDX(input, outloc, method, DDXDispatchExprTest::x_safe_region)}; + const auto expected = + expectedDDX(input, outloc, toString(method), DDXDispatchExprTest::x_safe_region); + + EXPECT_EQ(actual.getLocation(), outloc); + EXPECT_TRUE(IsFieldEqual(actual, expected, DDXDispatchExprTest::x_safe_region)); +} + +TEST_P(DDXDispatchExprMethodTest, MatchesDDX) { + auto input = makeTestField(mesh_staggered, CELL_CENTRE); + + const auto actual = + Field3D{DDX(input, CELL_CENTRE, GetParam(), DDXDispatchExprTest::x_safe_region)}; + const auto expected = expectedDDX(input, CELL_CENTRE, toString(GetParam()), + DDXDispatchExprTest::x_safe_region); + + EXPECT_EQ(actual.getLocation(), CELL_CENTRE); + EXPECT_TRUE(IsFieldEqual(actual, expected, DDXDispatchExprTest::x_safe_region)); +} + TEST_P(DDZDispatchExprParamTest, MatchesDDZ) { const auto [inloc, outloc, method] = GetParam(); auto input = makeTestField(mesh_staggered, inloc); @@ -61,6 +152,71 @@ TEST_P(DDZDispatchExprParamTest, MatchesDDZ) { EXPECT_TRUE(IsFieldEqual(actual, expected, "RGN_NOBNDRY")); } +TEST_F(DDXDispatchExprTest, UsesRequestedRegion) { + auto input = makeTestField(mesh_staggered, CELL_CENTRE); + + const auto actual = Field3D{DDX(input, CELL_XLOW, DIFF_C2, x_safe_region)}; + const auto expected = expectedDDX(input, CELL_XLOW, "C2", x_safe_region); + + EXPECT_EQ(actual.getLocation(), CELL_XLOW); + EXPECT_TRUE(IsFieldEqual(actual, expected, x_safe_region)); +} + +TEST_F(DDXDispatchExprTest, RejectsUnsupportedMethods) { + auto input = makeTestField(mesh_staggered, CELL_CENTRE); + + EXPECT_THROW((void)DDX(input, CELL_DEFAULT, DIFF_U1), BoutException); +} + +TEST_F(DDXDispatchExprTest, RejectsUnsupportedShearMethods) { + auto input = makeTestField(mesh_staggered, CELL_CENTRE); + + mesh_staggered->IncIntShear = true; + + EXPECT_THROW((void)DDX(input, CELL_DEFAULT, DIFF_W2), BoutException); +} + +TEST_F(DDXDispatchExprTest, ResolvesCentredDefaultMethodFromMesh) { + auto input = makeTestField(mesh_staggered, CELL_CENTRE); + + Options diff_options{{"ddx", {{"first", "C4"}}}, {"ddxstag", {{"first", "C2"}}}}; + static_cast(mesh_staggered)->initDerivs(&diff_options); + + EXPECT_EQ(mesh_staggered->getDefaultMethod(DIRECTION::X, DERIV::Standard), DIFF_C4); + const auto centre_actual = + Field3D{DDX(input, CELL_CENTRE, DIFF_DEFAULT, x_safe_region)}; + const auto centre_expected = expectedDDX(input, CELL_CENTRE, "C4", x_safe_region); + EXPECT_TRUE(IsFieldEqual(centre_actual, centre_expected, x_safe_region)); +} + +TEST_F(DDXDispatchExprTest, ResolvesStaggeredDefaultMethodFromMesh) { + auto input = makeTestField(mesh_staggered, CELL_CENTRE); + + Options diff_options{{"ddx", {{"first", "C4"}}}, {"ddxstag", {{"first", "C2"}}}}; + static_cast(mesh_staggered)->initDerivs(&diff_options); + + EXPECT_EQ(mesh_staggered->getDefaultMethod(DIRECTION::X, DERIV::Standard, STAGGER::C2L), + DIFF_C2); + const auto staggered_actual = + Field3D{DDX(input, CELL_XLOW, DIFF_DEFAULT, x_safe_region)}; + const auto staggered_expected = expectedDDX(input, CELL_XLOW, "C2", x_safe_region); + EXPECT_TRUE(IsFieldEqual(staggered_actual, staggered_expected, x_safe_region)); +} + +TEST_F(DDXDispatchExprTest, IncludesIntegratedShearCorrection) { + auto input = makeTestField(mesh_staggered, CELL_CENTRE); + + mesh_staggered->IncIntShear = true; + auto torsion = bout::FieldMetric(0.125, mesh_staggered); + mesh_staggered->getCoordinates()->setIntShiftTorsion(torsion); + CoordinatesAccessor::clear(mesh_staggered->getCoordinates()); + + const auto actual = Field3D{DDX(input, CELL_CENTRE, DIFF_C2, x_safe_region)}; + const auto expected = expectedDDX(input, CELL_CENTRE, "C2", x_safe_region); + + EXPECT_TRUE(IsFieldEqual(actual, expected, x_safe_region)); +} + TEST_F(DDZDispatchExprTest, UsesRequestedRegion) { auto input = makeTestField(mesh_staggered, CELL_CENTRE); From 23d308b0db34250de16cb92f6b63fcf7fe9b0288 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Mon, 31 Aug 2026 12:41:10 -0700 Subject: [PATCH 12/18] DDZ: Restore eager DDZ with FFT method DDZ is restored to the previous (eager) implementation that can dispatch to the FFT method. The new lazy implementation is `DDZ_stencil` because it is limited to stencil-based methods. --- include/bout/derivs.hxx | 24 ++++++++++++++++-- include/bout/stencil_expr.hxx | 18 ++++++------- src/sys/derivs.cxx | 25 ++++++++++++++++--- tests/unit/include/bout/test_stencil_expr.cxx | 10 ++++---- 4 files changed, 57 insertions(+), 20 deletions(-) diff --git a/include/bout/derivs.hxx b/include/bout/derivs.hxx index 4684a95a3e..4dd53f7f5e 100644 --- a/include/bout/derivs.hxx +++ b/include/bout/derivs.hxx @@ -125,10 +125,30 @@ bout::FieldMetric DDZ(const Field2D& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", const std::string& region = "RGN_NOBNDRY"); +/// Calculate first partial derivative in Z +/// +/// \f$\partial / \partial z\f$ +/// +/// @param[in] f The field to be differentiated +/// @param[in] outloc The cell location where the result is desired. If +/// staggered grids is not enabled then this has no effect +/// If not given, defaults to CELL_DEFAULT +/// @param[in] method Differencing method to use. This overrides the default +/// If not given, defaults to DIFF_DEFAULT +/// @param[in] region What region is expected to be calculated +/// If not given, defaults to RGN_NOBNDRY +Field3D DDZ(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT, + DIFF_METHOD method = DIFF_DEFAULT, const std::string& region = "RGN_NOBNDRY"); + +/// Compatibility overload for string-based callers. +Field3D DDZ(const Field3D& f, CELL_LOC outloc, const std::string& method, + const std::string& region = "RGN_NOBNDRY"); + /// Calculate the Z derivative using the FFT-based implementation. /// -/// This is kept separate from the lazy-expression DDZ overload because FFT -/// differentiation cannot be represented by the stencil-expression machinery. +/// This is kept separate from the stencil-expression `DDZ_stencil` overload +/// because FFT differentiation cannot be represented by the stencil-expression +/// machinery. Field3D DDZ_FFT(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& region = "RGN_NOBNDRY"); diff --git a/include/bout/stencil_expr.hxx b/include/bout/stencil_expr.hxx index 36f8d87c72..ed27a890dc 100644 --- a/include/bout/stencil_expr.hxx +++ b/include/bout/stencil_expr.hxx @@ -692,10 +692,10 @@ inline bout::stencil::DDZExprC4 DDZ_C4(const Field3D& f) { f.getMesh()->getRegion("RGN_NOBNDRY")}; } -inline bout::stencil::DDZDispatchExpr DDZ(const Field3D& f, - CELL_LOC outloc = CELL_DEFAULT, - DIFF_METHOD method = DIFF_DEFAULT, - const std::string& region = "RGN_NOBNDRY") { +inline bout::stencil::DDZDispatchExpr +DDZ_stencil(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT, + DIFF_METHOD method = DIFF_DEFAULT, + const std::string& region = "RGN_NOBNDRY") { checkData(f); const auto resolved_outloc = (outloc == CELL_DEFAULT) ? f.getLocation() : outloc; @@ -707,7 +707,7 @@ inline bout::stencil::DDZDispatchExpr DDZ(const Field3D& f, } if ((method != DIFF_C2) && (method != DIFF_C4)) { - throw BoutException("DDZ only supports DIFF_C2 and DIFF_C4, got {:s}", + throw BoutException("DDZ_stencil only supports DIFF_C2 and DIFF_C4, got {:s}", toString(method)); } @@ -726,10 +726,10 @@ inline bout::stencil::DDZDispatchExpr DDZ(const Field3D& f, f.getMesh()->getRegion(region)}; } -inline bout::stencil::DDZDispatchExpr DDZ(const Field3D& f, CELL_LOC outloc, - const std::string& method, - const std::string& region = "RGN_NOBNDRY") { - return DDZ(f, outloc, parseField3DMethodString(method), region); +inline bout::stencil::DDZDispatchExpr +DDZ_stencil(const Field3D& f, CELL_LOC outloc, const std::string& method, + const std::string& region = "RGN_NOBNDRY") { + return DDZ_stencil(f, outloc, parseField3DMethodString(method), region); } inline bout::stencil::BracketArakawaExpr bracket_arakawa(const Field3D& f, diff --git a/src/sys/derivs.cxx b/src/sys/derivs.cxx index fa61b73b71..0e2ca41345 100644 --- a/src/sys/derivs.cxx +++ b/src/sys/derivs.cxx @@ -84,9 +84,24 @@ bout::FieldMetric DDY(const Field2D& f, CELL_LOC outloc, const std::string& meth ////////////// Z DERIVATIVE ///////////////// +Field3D DDZ(const Field3D& f, CELL_LOC outloc, DIFF_METHOD method, + const std::string& region) { + const auto resolved_outloc = (outloc == CELL_DEFAULT) ? f.getLocation() : outloc; + return bout::derivatives::index::DDZ(f, outloc, toString(method), region) + / interp_to(f.getCoordinates(resolved_outloc)->dz(), resolved_outloc, region); +} + +Field3D DDZ(const Field3D& f, CELL_LOC outloc, const std::string& method, + const std::string& region) { + const auto resolved_outloc = (outloc == CELL_DEFAULT) ? f.getLocation() : outloc; + return bout::derivatives::index::DDZ(f, outloc, method, region) + / interp_to(f.getCoordinates(resolved_outloc)->dz(), resolved_outloc, region); +} + Field3D DDZ_FFT(const Field3D& f, CELL_LOC outloc, const std::string& region) { + const auto resolved_outloc = (outloc == CELL_DEFAULT) ? f.getLocation() : outloc; return bout::derivatives::index::DDZ(f, outloc, "FFT", region) - / f.getCoordinates(outloc)->dz(); + / f.getCoordinates(resolved_outloc)->dz(); } bout::FieldMetric DDZ(const Field2D& f, CELL_LOC UNUSED(outloc), @@ -383,13 +398,15 @@ bout::FieldMetric D2DYDZ(const Field2D& f, CELL_LOC outloc, #endif } -Field3D D2DYDZ(const Field3D& f, CELL_LOC outloc, - [[maybe_unused]] const std::string& method, const std::string& region) { +Field3D D2DYDZ(const Field3D& f, CELL_LOC outloc, const std::string& method, + const std::string& region) { // If staggering in z, take y-derivative at f's location. const auto y_location = (outloc == CELL_ZLOW or f.getLocation() == CELL_ZLOW) ? CELL_DEFAULT : outloc; + const auto y_method = + (parseField3DMethodString(method) == DIFF_FFT) ? std::string{"DEFAULT"} : method; - return DDZ(DDY(f, y_location, method, region), outloc, method, region); + return DDZ(DDY(f, y_location, y_method, region), outloc, method, region); } /******************************************************************************* diff --git a/tests/unit/include/bout/test_stencil_expr.cxx b/tests/unit/include/bout/test_stencil_expr.cxx index bf9bf67f6d..9cfe7b5f97 100644 --- a/tests/unit/include/bout/test_stencil_expr.cxx +++ b/tests/unit/include/bout/test_stencil_expr.cxx @@ -145,7 +145,7 @@ TEST_P(DDZDispatchExprParamTest, MatchesDDZ) { const auto [inloc, outloc, method] = GetParam(); auto input = makeTestField(mesh_staggered, inloc); - const auto actual = Field3D{DDZ(input, outloc, method)}; + const auto actual = Field3D{DDZ_stencil(input, outloc, method)}; const auto expected = DDZ(input, outloc, toString(method)); EXPECT_EQ(actual.getLocation(), outloc); @@ -220,7 +220,7 @@ TEST_F(DDXDispatchExprTest, IncludesIntegratedShearCorrection) { TEST_F(DDZDispatchExprTest, UsesRequestedRegion) { auto input = makeTestField(mesh_staggered, CELL_CENTRE); - const auto actual = Field3D{DDZ(input, CELL_ZLOW, DIFF_C2, "RGN_ALL")}; + const auto actual = Field3D{DDZ_stencil(input, CELL_ZLOW, DIFF_C2, "RGN_ALL")}; const auto expected = DDZ(input, CELL_ZLOW, "C2", "RGN_ALL"); EXPECT_EQ(actual.getLocation(), CELL_ZLOW); @@ -230,7 +230,7 @@ TEST_F(DDZDispatchExprTest, UsesRequestedRegion) { TEST_F(DDZDispatchExprTest, RejectsUnsupportedMethods) { auto input = makeTestField(mesh_staggered, CELL_CENTRE); - EXPECT_THROW((void)DDZ(input, CELL_DEFAULT, DIFF_FFT), BoutException); + EXPECT_THROW((void)DDZ_stencil(input, CELL_DEFAULT, DIFF_FFT), BoutException); } TEST_F(DDZDispatchExprTest, ResolvesDefaultMethodFromMesh) { @@ -239,11 +239,11 @@ TEST_F(DDZDispatchExprTest, ResolvesDefaultMethodFromMesh) { Options diff_options{{"ddz", {{"first", "C4"}}}, {"ddzstag", {{"first", "C2"}}}}; static_cast(mesh_staggered)->initDerivs(&diff_options); - const auto centre_actual = Field3D{DDZ(input, CELL_CENTRE, DIFF_DEFAULT)}; + const auto centre_actual = Field3D{DDZ_stencil(input, CELL_CENTRE, DIFF_DEFAULT)}; const auto centre_expected = DDZ(input, CELL_CENTRE, "C4"); EXPECT_TRUE(IsFieldEqual(centre_actual, centre_expected, "RGN_NOBNDRY")); - const auto staggered_actual = Field3D{DDZ(input, CELL_ZLOW, DIFF_DEFAULT)}; + const auto staggered_actual = Field3D{DDZ_stencil(input, CELL_ZLOW, DIFF_DEFAULT)}; const auto staggered_expected = DDZ(input, CELL_ZLOW, "C2"); EXPECT_TRUE(IsFieldEqual(staggered_actual, staggered_expected, "RGN_NOBNDRY")); } From e00fcd19ee41cc8248c5c8feac84dda8a9075190 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Mon, 31 Aug 2026 13:39:35 -0700 Subject: [PATCH 13/18] Naulin Laplace: Revert to eager FFT method --- .../laplace/impls/naulin/naulin_laplace.cxx | 23 ++++++++++++------- .../test-naulin-laplace/data/BOUT.inp | 4 ++-- .../test_naulin_laplace.cxx | 4 ++-- 3 files changed, 19 insertions(+), 12 deletions(-) diff --git a/src/invert/laplace/impls/naulin/naulin_laplace.cxx b/src/invert/laplace/impls/naulin/naulin_laplace.cxx index f41b89bb36..3e66ff55ee 100644 --- a/src/invert/laplace/impls/naulin/naulin_laplace.cxx +++ b/src/invert/laplace/impls/naulin/naulin_laplace.cxx @@ -14,9 +14,6 @@ * and so that all Neumann boundary conditions can be used at least when * DC(A/D)!=0. * - * CHANGELOG - * ========= - * ************************************************************************** * Copyright 2018 - 2026 BOUT++ contributors * @@ -43,17 +40,27 @@ #if not BOUT_USE_METRIC_3D +#include +#include #include #include #include #include +#include +#include #include #include #include +#include #include #include "naulin_laplace.hxx" +#include + +#include +#include + LaplaceNaulin::LaplaceNaulin(Options* opt, const CELL_LOC loc, Mesh* mesh_in, Solver* UNUSED(solver)) : Laplacian(opt, loc, mesh_in), Acoef(0.0), C1coef(1.0), C2coef(0.0), Dcoef(1.0), @@ -142,13 +149,13 @@ Field3D LaplaceNaulin::solve(const Field3D& rhs, const Field3D& x0) { Field3D C1TimesD = C1coef * Dcoef; // This is needed several times // x-component of 1./(C1*D) * Grad_perp(C2) - Field3D coef_x = DDX(C2coef, location, "C2") / C1TimesD; + Field3D coef_x = DDX(C2coef, location, DIFF_C2) / C1TimesD; // y-component of 1./(C1*D) * Grad_perp(C2) Field3D coef_y = DDY(C2coef, location, "C2") / C1TimesD; // z-component of 1./(C1*D) * Grad_perp(C2) - Field3D coef_z = DDZ(C2coef, location, DIFF_C4) / C1TimesD; + Field3D coef_z = DDZ(C2coef, location, DIFF_FFT) / C1TimesD; Field3D AOverD = Acoef / Dcoef; @@ -186,8 +193,8 @@ Field3D LaplaceNaulin::solve(const Field3D& rhs, const Field3D& x0) { auto calc_b_guess = [&](const Field3D& x_in) { // Derivatives of x - Field3D ddx_x = DDX(x_in, location, "C2"); - Field3D ddz_x = DDZ(x_in, location, DIFF_C4); + const Field3D ddx_x = DDX(x_in, location, DIFF_C2); + const Field3D ddz_x = DDZ(x_in, location, DIFF_FFT); return rhsOverD - (coords->g11() * coef_x_AC * ddx_x + coords->g33() * coef_z * ddz_x + coords->g13() * (coef_x_AC * ddz_x + coef_z * ddx_x)) @@ -211,7 +218,7 @@ Field3D LaplaceNaulin::solve(const Field3D& rhs, const Field3D& x0) { return std::make_pair(b, x); }; - Field3D b = calc_b_guess(x0); + const Field3D b = calc_b_guess(x0); // Need to make a copy of x0 here to make sure we don't change x0 auto b_x_pair = calc_b_x_pair(b, x0); auto b_x_pair_old = b_x_pair; diff --git a/tests/integrated/test-naulin-laplace/data/BOUT.inp b/tests/integrated/test-naulin-laplace/data/BOUT.inp index b171d1755d..eb4ea02b72 100644 --- a/tests/integrated/test-naulin-laplace/data/BOUT.inp +++ b/tests/integrated/test-naulin-laplace/data/BOUT.inp @@ -23,8 +23,8 @@ g_13 = 2.38908100128/(1.+.01*y*dy) g_23 = 2.14121198654175/(1.+.01*y*dy) [mesh:ddz] -first = C4 -second = C4 +first = fft +second = fft ############################################# diff --git a/tests/integrated/test-naulin-laplace/test_naulin_laplace.cxx b/tests/integrated/test-naulin-laplace/test_naulin_laplace.cxx index 456f5fd72d..a4af049d14 100644 --- a/tests/integrated/test-naulin-laplace/test_naulin_laplace.cxx +++ b/tests/integrated/test-naulin-laplace/test_naulin_laplace.cxx @@ -305,8 +305,8 @@ int main(int argc, char** argv) { Field3D this_Grad_perp_dot_Grad_perp(const Field3D& f, const Field3D& g) { const auto* coords = f.getCoordinates(); Field3D result = coords->g11() * ::DDX(f) * ::DDX(g) - + coords->g33() * ::DDZ_FFT(f) * ::DDZ_FFT(g) - + coords->g13() * (DDX(f) * DDZ_FFT(g) + DDZ_FFT(f) * DDX(g)); + + coords->g33() * ::DDZ(f) * ::DDZ(g) + + coords->g13() * (DDX(f) * DDZ(g) + DDZ(f) * DDX(g)); return result; } From f0173ddfacdfc734871efd45719e348609c69f3a Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Mon, 31 Aug 2026 16:55:36 -0700 Subject: [PATCH 14/18] Revert test-invpar input and address clang-tidy comments --- examples/hasegawa-wakatani/hw.cxx | 9 +++- include/bout/deriv_store.hxx | 42 +++++++++---------- include/bout/physicsmodel.hxx | 2 +- .../laplace/impls/naulin/naulin_laplace.cxx | 13 +++--- src/mesh/coordinates_accessor.cxx | 12 ++++-- src/sys/derivs.cxx | 1 + tests/integrated/test-invpar/data/BOUT.inp | 4 +- 7 files changed, 47 insertions(+), 36 deletions(-) diff --git a/examples/hasegawa-wakatani/hw.cxx b/examples/hasegawa-wakatani/hw.cxx index 45864fa7eb..1eefd54fda 100644 --- a/examples/hasegawa-wakatani/hw.cxx +++ b/examples/hasegawa-wakatani/hw.cxx @@ -1,9 +1,14 @@ +#include #include +#include #include #include #include #include +#include + +#include class HW : public PhysicsModel { private: @@ -110,8 +115,8 @@ class HW : public PhysicsModel { nonzonal_phi -= averageY(DC(phi)); } - ddt(n) = - -bracket_arakawa(phi, n) + alpha * (nonzonal_phi - nonzonal_n) - kappa * DDZ(phi); + ddt(n) = -bracket_arakawa(phi, n) + alpha * (nonzonal_phi - nonzonal_n) + - kappa * DDZ_stencil(phi); ddt(vort) = -bracket_arakawa(phi, vort) + alpha * (nonzonal_phi - nonzonal_n); diff --git a/include/bout/deriv_store.hxx b/include/bout/deriv_store.hxx index 0a995fd7a3..078e59f599 100644 --- a/include/bout/deriv_store.hxx +++ b/include/bout/deriv_store.hxx @@ -4,10 +4,9 @@ * Definition of derivative methods storage class * ************************************************************************** - * Copyright 2018 - * D.Dickinson, P.Hill, B.Dudson + * Copyright 2018 - 2026 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -26,22 +25,24 @@ * **************************************************************************/ -#ifndef __DERIV_STORE_HXX__ -#define __DERIV_STORE_HXX__ +#ifndef DERIV_STORE_HXX +#define DERIV_STORE_HXX #include #include #include #include +#include #include #include -#include "bout/field3d.hxx" -#include - #include #include +#include #include +#include +#include +#include /// Here we have a templated singleton that is used to store DerivativeFunctions /// for all types of derivatives. It is templated on the FieldType (2D or 3D) as @@ -100,9 +101,8 @@ struct DerivativeStore { auto key = getKey(direction, stagger, toString(derivType)); if (isEmpty(key)) { return std::set{}; - } else { - return registeredMethods.at(key); } + return registeredMethods.at(key); }; /// Outputs a list of all registered method names for the @@ -313,18 +313,18 @@ struct DerivativeStore { // To replicate the existing behaviour we first search for a section called //"dd?" and if the option isn't in there we search a section called "diff" - auto backupSection = options->getSection("diff"); + auto* backupSection = options->getSection("diff"); - std::map directions = {{DIRECTION::X, "ddx"}, - {DIRECTION::Y, "ddy"}, - {DIRECTION::YOrthogonal, "ddy"}, - {DIRECTION::Z, "ddz"}}; + const std::map directions = {{DIRECTION::X, "ddx"}, + {DIRECTION::Y, "ddy"}, + {DIRECTION::YOrthogonal, "ddy"}, + {DIRECTION::Z, "ddz"}}; - std::map derivTypes = {{DERIV::Standard, "first"}, - {DERIV::StandardSecond, "second"}, - {DERIV::StandardFourth, "fourth"}, - {DERIV::Upwind, "upwind"}, - {DERIV::Flux, "flux"}}; + const std::map derivTypes = {{DERIV::Standard, "first"}, + {DERIV::StandardSecond, "second"}, + {DERIV::StandardFourth, "fourth"}, + {DERIV::Upwind, "upwind"}, + {DERIV::Flux, "flux"}}; for (const auto& direction : directions) { for (const auto& deriv : derivTypes) { @@ -344,7 +344,7 @@ struct DerivativeStore { //------------------------------------------------------------- // The direction specific section to consider - auto specificSection = options->getSection(direction.second); + auto* specificSection = options->getSection(direction.second); // Find the appropriate value for theDefault either from // the input file or if not found then use the value in diff --git a/include/bout/physicsmodel.hxx b/include/bout/physicsmodel.hxx index f9160f5709..7091493947 100644 --- a/include/bout/physicsmodel.hxx +++ b/include/bout/physicsmodel.hxx @@ -428,7 +428,7 @@ private: #define BOUTMAIN(ModelClass) \ int main(int argc, char** argv) { \ try { \ - int init_err = BoutInitialise(argc, argv); \ + const int init_err = BoutInitialise(argc, argv); \ if (init_err < 0) { \ return 0; \ } \ diff --git a/src/invert/laplace/impls/naulin/naulin_laplace.cxx b/src/invert/laplace/impls/naulin/naulin_laplace.cxx index 3e66ff55ee..1e9021c06d 100644 --- a/src/invert/laplace/impls/naulin/naulin_laplace.cxx +++ b/src/invert/laplace/impls/naulin/naulin_laplace.cxx @@ -47,6 +47,7 @@ #include #include #include +#include #include #include #include @@ -146,24 +147,24 @@ Field3D LaplaceNaulin::solve(const Field3D& rhs, const Field3D& x0) { Field3D rhsOverD = rhs / Dcoef; - Field3D C1TimesD = C1coef * Dcoef; // This is needed several times + const Field3D C1TimesD = C1coef * Dcoef; // This is needed several times // x-component of 1./(C1*D) * Grad_perp(C2) - Field3D coef_x = DDX(C2coef, location, DIFF_C2) / C1TimesD; + const Field3D coef_x = DDX(C2coef, location, DIFF_C2) / C1TimesD; // y-component of 1./(C1*D) * Grad_perp(C2) - Field3D coef_y = DDY(C2coef, location, "C2") / C1TimesD; + const Field3D coef_y = DDY(C2coef, location, "C2") / C1TimesD; // z-component of 1./(C1*D) * Grad_perp(C2) Field3D coef_z = DDZ(C2coef, location, DIFF_FFT) / C1TimesD; - Field3D AOverD = Acoef / Dcoef; + const Field3D AOverD = Acoef / Dcoef; // Split coefficients into DC and AC parts so that delp2solver can use DC part. // This allows all-Neumann boundary conditions as long as AOverD_DC is non-zero - Field2D C1coefTimesD_DC = DC(C1TimesD); - Field2D C2coef_DC = DC(C2coef); + const Field2D C1coefTimesD_DC = DC(C1TimesD); + const Field2D C2coef_DC = DC(C2coef); // Our naming is slightly misleading here, as coef_x_AC may actually have a // DC component, as the AC components of C2coef and C1coefTimesD are not diff --git a/src/mesh/coordinates_accessor.cxx b/src/mesh/coordinates_accessor.cxx index 272c67a22e..a7c8c71023 100644 --- a/src/mesh/coordinates_accessor.cxx +++ b/src/mesh/coordinates_accessor.cxx @@ -1,5 +1,9 @@ #include "bout/coordinates_accessor.hxx" +#include "bout/array.hxx" +#include "bout/assert.hxx" +#include "bout/bout_types.hxx" #include "bout/build_defines.hxx" +#include "bout/coordinates.hxx" #include "bout/macro_for_each.hxx" #include "bout/mesh.hxx" @@ -40,10 +44,10 @@ CoordinatesAccessor::CoordinatesAccessor(const Coordinates* coords) { // Copy data from Coordinates variable into data array // Uses the symbol to look up the corresponding Offset -#define COPY_STRIPE1(symbol) \ - if (coords->symbol().isAllocated()) { \ - data[stripe_size * ind.ind + static_cast(Offset::symbol)] = \ - coords->symbol()[ind]; \ +#define COPY_STRIPE1(symbol) \ + if (coords->symbol().isAllocated()) { \ + data[(stripe_size * ind.ind) + static_cast(Offset::symbol)] = \ + coords->symbol()[ind]; \ } // Implement copy for each argument diff --git a/src/sys/derivs.cxx b/src/sys/derivs.cxx index 0e2ca41345..f25787af61 100644 --- a/src/sys/derivs.cxx +++ b/src/sys/derivs.cxx @@ -50,6 +50,7 @@ #include #include #include +#include #include #include diff --git a/tests/integrated/test-invpar/data/BOUT.inp b/tests/integrated/test-invpar/data/BOUT.inp index b7935fc298..dfe708b688 100644 --- a/tests/integrated/test-invpar/data/BOUT.inp +++ b/tests/integrated/test-invpar/data/BOUT.inp @@ -19,8 +19,8 @@ g_22 = 1.1 + 0.2*cos(y) g22 = 1 / g_22 [mesh:ddz] -first = c4 -second = c4 +first = fft +second = fft [solver] # Test of parallel tridiagonal solver From 3fc221a31209328e3a6146a77f35ecf4ef1001ff Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Mon, 31 Aug 2026 17:32:58 -0700 Subject: [PATCH 15/18] Bxy yup/down guards in CoordinateAccessor Ensure that regions are propagated through index derivatives, and ensure that Bxy has yup/down slices before trying to access them. --- include/bout/index_derivs_interface.hxx | 10 ++++++++-- src/mesh/coordinates_accessor.cxx | 4 ++-- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/include/bout/index_derivs_interface.hxx b/include/bout/index_derivs_interface.hxx index cbb1e5ac53..3f20df857b 100644 --- a/include/bout/index_derivs_interface.hxx +++ b/include/bout/index_derivs_interface.hxx @@ -85,7 +85,9 @@ T flowDerivative(const T& vel, const T& f, CELL_LOC outloc, const std::string& m const int nPoint = localmesh->getNpoints(direction); if (nPoint == 1) { - return zeroFrom(f).setLocation(outloc); + T result{zeroFrom(f).setLocation(outloc)}; + result.setRegion(region); + return result; } // Lookup the method @@ -94,6 +96,7 @@ T flowDerivative(const T& vel, const T& f, CELL_LOC outloc, const std::string& m // Create the result field T result{emptyFrom(f).setLocation(outloc)}; + result.setRegion(region); // Apply method derivativeMethod(vel, f, result, region); @@ -146,7 +149,9 @@ T standardDerivative(const T& f, CELL_LOC outloc, const std::string& method, const int nPoint = localmesh->getNpoints(direction); if (nPoint == 1) { - return zeroFrom(f).setLocation(outloc); + T result{zeroFrom(f).setLocation(outloc)}; + result.setRegion(region); + return result; } // Lookup the method @@ -155,6 +160,7 @@ T standardDerivative(const T& f, CELL_LOC outloc, const std::string& method, // Create the result field T result{emptyFrom(f).setLocation(outloc)}; + result.setRegion(region); // Apply method derivativeMethod(f, result, region); diff --git a/src/mesh/coordinates_accessor.cxx b/src/mesh/coordinates_accessor.cxx index a7c8c71023..c3a63637a3 100644 --- a/src/mesh/coordinates_accessor.cxx +++ b/src/mesh/coordinates_accessor.cxx @@ -66,11 +66,11 @@ CoordinatesAccessor::CoordinatesAccessor(const Coordinates* coords) { if (coords->Bxy().isAllocated()) { data[(stripe_size * ind.ind) + static_cast(Offset::B)] = coords->Bxy()[ind]; - if (coords->Bxy().yup().isAllocated()) { + if (coords->Bxy().hasParallelSlices() && coords->Bxy().yup().isAllocated()) { data[stripe_size * ind.ind + static_cast(Offset::Byup)] = coords->Bxy().yup()[ind]; } - if (coords->Bxy().ydown().isAllocated()) { + if (coords->Bxy().hasParallelSlices() && coords->Bxy().ydown().isAllocated()) { data[stripe_size * ind.ind + static_cast(Offset::Bydown)] = coords->Bxy().ydown()[ind]; } From ff7f003e6bec81b3983e15789a5c68753da391c9 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Mon, 31 Aug 2026 21:23:28 -0700 Subject: [PATCH 16/18] DDY_stencil: Lazy evaluation of Y derivative Requires yup/ydown fields so that the to/fromFieldAligned transformations are not needed. --- include/bout/stencil_expr.hxx | 145 ++++++++++++++++++ tests/unit/include/bout/test_stencil_expr.cxx | 108 +++++++++++++ 2 files changed, 253 insertions(+) diff --git a/include/bout/stencil_expr.hxx b/include/bout/stencil_expr.hxx index ed27a890dc..61ac66a5d8 100644 --- a/include/bout/stencil_expr.hxx +++ b/include/bout/stencil_expr.hxx @@ -316,6 +316,77 @@ struct DDX_Dispatch_Op { } }; +struct DDY_C2_Op { + CoordinatesAccessor coords; + int nz{0}; + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal operator()(int idx, const LView& lhs, + const RView&) const { + const int iyp = i_yp(idx, nz); + const int iym = i_ym(idx, nz); + + return 0.5 * (lhs.yup(0)(iyp) - lhs.ydown(0)(iym)) / coords.dy(idx); + } +}; + +struct DDY_C4_Op { + CoordinatesAccessor coords; + int nz{0}; + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal operator()(int idx, const LView& lhs, + const RView&) const { + const int iyp = i_yp(idx, nz); + const int iym = i_ym(idx, nz); + const int iyp2 = i_yp(iyp, nz); + const int iym2 = i_ym(iym, nz); + + return (-lhs.yup(1)(iyp2) + 8.0 * lhs.yup(0)(iyp) - 8.0 * lhs.ydown(0)(iym) + + lhs.ydown(1)(iym2)) + / (12.0 * coords.dy(idx)); + } +}; + +struct DDY_Dispatch_Op { + CoordinatesAccessor coords; + int nz{0}; + DIFF_METHOD method{DIFF_DEFAULT}; + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c2(int idx, const LView& lhs) const { + const int iyp = i_yp(idx, nz); + const int iym = i_ym(idx, nz); + + return 0.5 * (lhs.yup(0)(iyp) - lhs.ydown(0)(iym)) / coords.dy(idx); + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal apply_c4(int idx, const LView& lhs) const { + const int iyp = i_yp(idx, nz); + const int iym = i_ym(idx, nz); + const int iyp2 = i_yp(iyp, nz); + const int iym2 = i_ym(iym, nz); + + return (-lhs.yup(1)(iyp2) + 8.0 * lhs.yup(0)(iyp) - 8.0 * lhs.ydown(0)(iym) + + lhs.ydown(1)(iym2)) + / (12.0 * coords.dy(idx)); + } + + template + BOUT_HOST_DEVICE BOUT_FORCEINLINE BoutReal operator()(int idx, const LView& lhs, + const RView&) const { + switch (method) { + case DIFF_C2: + return apply_c2(idx, lhs); + case DIFF_C4: + return apply_c4(idx, lhs); + default: + return 0.0; + } + } +}; + struct DDZ_C2_Op { CoordinatesAccessor coords; int nz{0}; @@ -490,6 +561,9 @@ struct BracketArakawaOp { using DDZExprC2 = BinaryExpr; using DDZExprC4 = BinaryExpr; using DDZDispatchExpr = BinaryExpr; +using DDYExprC2 = BinaryExpr; +using DDYExprC4 = BinaryExpr; +using DDYDispatchExpr = BinaryExpr; using DDXExprC2 = BinaryExpr; using DDXExprC4 = BinaryExpr; using DDXDispatchExpr = BinaryExpr; @@ -524,6 +598,17 @@ inline int requiredDDXGuards(DIFF_METHOD method, STAGGER stagger) { } } +inline int requiredDDYParallelSlices(DIFF_METHOD method) { + switch (method) { + case DIFF_C2: + return 1; + case DIFF_C4: + return 2; + default: + return -1; + } +} + inline DIFF_METHOD parseField3DMethodString(const std::string& method) { auto normalized = uppercase(method); if (normalized.rfind("DIFF_", 0) == 0) { @@ -660,6 +745,66 @@ inline bout::stencil::DDXDispatchExpr DDX(const Field3D& f, CELL_LOC outloc, return DDX(f, outloc, parseField3DMethodString(method), region); } +inline bout::stencil::DDYDispatchExpr +DDY_stencil(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT, + DIFF_METHOD method = DIFF_DEFAULT, + const std::string& region = "RGN_NOBNDRY") { + checkData(f); + + if (!f.hasParallelSlices()) { + throw BoutException("DDY_stencil requires parallel slices. Use eager DDY for " + "field-aligned transforms or communicate/apply parallel " + "boundaries before calling DDY_stencil."); + } + + ASSERT1(f.getDirectionY() == YDirectionType::Standard); + + const auto resolved_outloc = (outloc == CELL_DEFAULT) ? f.getLocation() : outloc; + const auto stagger = + f.getMesh()->getStagger(f.getLocation(), resolved_outloc, CELL_YLOW); + + if (stagger != STAGGER::None) { + throw BoutException("DDY_stencil currently only supports unstaggered outputs, got " + "{:s} to {:s}", + toString(f.getLocation()), toString(resolved_outloc)); + } + + if (method == DIFF_DEFAULT) { + method = f.getMesh()->getDefaultMethod(DIRECTION::Y, DERIV::Standard, stagger); + } + + const auto required_slices = requiredDDYParallelSlices(method); + if (required_slices < 0) { + throw BoutException("DDY_stencil only supports DIFF_C2 and DIFF_C4, got {:s}", + toString(method)); + } + + if (static_cast(f.numberParallelSlices()) < required_slices) { + throw BoutException("DDY_stencil with {:s} requires {:d} parallel slice pair(s), " + "but field only has {:d}", + toString(method), required_slices, f.numberParallelSlices()); + } + + const auto region_id = f.getMesh()->getRegionID(region); + + return bout::stencil::DDYDispatchExpr{ + static_cast(f), + static_cast(f), + bout::stencil::DDY_Dispatch_Op{ + CoordinatesAccessor{f.getCoordinates(resolved_outloc)}, f.getNz(), method}, + f.getMesh(), + resolved_outloc, + f.getDirections(), + region_id, + f.getMesh()->getRegion(region)}; +} + +inline bout::stencil::DDYDispatchExpr +DDY_stencil(const Field3D& f, CELL_LOC outloc, const std::string& method, + const std::string& region = "RGN_NOBNDRY") { + return DDY_stencil(f, outloc, parseField3DMethodString(method), region); +} + inline bout::stencil::DDZExprC2 DDZ_C2(const Field3D& f) { checkData(f); diff --git a/tests/unit/include/bout/test_stencil_expr.cxx b/tests/unit/include/bout/test_stencil_expr.cxx index 9cfe7b5f97..a58a64bb23 100644 --- a/tests/unit/include/bout/test_stencil_expr.cxx +++ b/tests/unit/include/bout/test_stencil_expr.cxx @@ -24,6 +24,28 @@ Field3D makeTestField(Mesh* mesh, CELL_LOC location) { return result; } +void fillTestField(Field3D& result) { + for (auto i : result.getRegion("RGN_ALL")) { + result[i] = 0.1 * i.x() + 0.2 * i.y() + 0.3 * i.z() + 0.05 * i.x() * i.z(); + } +} + +Field3DParallel makeParallelTestField(Mesh* mesh, CELL_LOC location) { + Field3DParallel result(mesh, location); + result.allocate(); + fillTestField(result); + result.splitParallelSlices(); + + for (size_t slice = 0; slice < result.numberParallelSlices(); ++slice) { + result.yup(slice).allocate(); + result.ydown(slice).allocate(); + fillTestField(result.yup(slice)); + fillTestField(result.ydown(slice)); + } + + return result; +} + Field3D expectedDDX(const Field3D& input, CELL_LOC outloc, const std::string& method, const std::string& region = "RGN_NOBNDRY") { const auto resolved_outloc = (outloc == CELL_DEFAULT) ? input.getLocation() : outloc; @@ -77,6 +99,25 @@ class DDZDispatchExprParamTest public ::testing::WithParamInterface> { }; +class DDYDispatchExprTest : public FakeMeshFixture {}; + +class DDYDispatchExprTwoSliceTest : public FakeMeshFixture_tmpl<5, 7, 7> { +public: + DDYDispatchExprTwoSliceTest() { + for (auto* current_mesh : {bout::globals::mesh, mesh_staggered}) { + current_mesh->ystart = 2; + current_mesh->yend = current_mesh->LocalNy - 3; + current_mesh->addRegion3D(y_safe_region, Region(0, current_mesh->LocalNx - 1, + 2, current_mesh->LocalNy - 3, + 0, current_mesh->LocalNz - 1, + current_mesh->LocalNy, + current_mesh->LocalNz)); + } + } + + static constexpr auto y_safe_region = "RGN_YSAFE"; +}; + std::string paramToString( const ::testing::TestParamInfo>& param) { const auto [inloc, outloc, method] = param.param; @@ -247,3 +288,70 @@ TEST_F(DDZDispatchExprTest, ResolvesDefaultMethodFromMesh) { const auto staggered_expected = DDZ(input, CELL_ZLOW, "C2"); EXPECT_TRUE(IsFieldEqual(staggered_actual, staggered_expected, "RGN_NOBNDRY")); } + +TEST_F(DDYDispatchExprTest, MatchesDDYC2) { + auto input = makeParallelTestField(mesh_staggered, CELL_CENTRE); + + const auto actual = Field3D{DDY_stencil(input, CELL_CENTRE, DIFF_C2)}; + const auto expected = DDY(input, CELL_CENTRE, "C2"); + + EXPECT_EQ(actual.getLocation(), CELL_CENTRE); + EXPECT_TRUE(IsFieldEqual(actual, expected, "RGN_NOBNDRY")); +} + +TEST_F(DDYDispatchExprTest, RejectsWithoutParallelSlices) { + auto input = makeTestField(mesh_staggered, CELL_CENTRE); + + EXPECT_THROW((void)DDY_stencil(input, CELL_CENTRE, DIFF_C2), BoutException); +} + +TEST_F(DDYDispatchExprTest, RejectsUnsupportedMethods) { + auto input = makeParallelTestField(mesh_staggered, CELL_CENTRE); + + EXPECT_THROW((void)DDY_stencil(input, CELL_CENTRE, DIFF_W3), BoutException); +} + +TEST_F(DDYDispatchExprTest, RejectsStaggeredOutputs) { + auto input = makeParallelTestField(mesh_staggered, CELL_CENTRE); + + EXPECT_THROW((void)DDY_stencil(input, CELL_YLOW, DIFF_C2), BoutException); +} + +TEST_F(DDYDispatchExprTest, RejectsC4WithOneSlicePair) { + auto input = makeParallelTestField(mesh_staggered, CELL_CENTRE); + + EXPECT_THROW((void)DDY_stencil(input, CELL_CENTRE, DIFF_C4), BoutException); +} + +TEST_F(DDYDispatchExprTest, ResolvesDefaultMethodFromMesh) { + auto input = makeParallelTestField(mesh_staggered, CELL_CENTRE); + + Options diff_options{{"ddy", {{"first", "C2"}}}}; + static_cast(mesh_staggered)->initDerivs(&diff_options); + + const auto actual = Field3D{DDY_stencil(input, CELL_CENTRE, DIFF_DEFAULT)}; + const auto expected = DDY(input, CELL_CENTRE, "C2"); + EXPECT_TRUE(IsFieldEqual(actual, expected, "RGN_NOBNDRY")); +} + +TEST_F(DDYDispatchExprTwoSliceTest, MatchesDDYC4WithTwoSlicePairs) { + auto input = makeParallelTestField(mesh_staggered, CELL_CENTRE); + + const auto actual = Field3D{DDY_stencil(input, CELL_CENTRE, DIFF_C4, y_safe_region)}; + const auto expected = DDY(input, CELL_CENTRE, "C4", y_safe_region); + + EXPECT_EQ(actual.getLocation(), CELL_CENTRE); + EXPECT_TRUE(IsFieldEqual(actual, expected, y_safe_region)); +} + +TEST_F(DDYDispatchExprTwoSliceTest, ResolvesDefaultMethodFromMesh) { + auto input = makeParallelTestField(mesh_staggered, CELL_CENTRE); + + Options diff_options{{"ddy", {{"first", "C4"}}}}; + static_cast(mesh_staggered)->initDerivs(&diff_options); + + const auto actual = + Field3D{DDY_stencil(input, CELL_CENTRE, DIFF_DEFAULT, y_safe_region)}; + const auto expected = DDY(input, CELL_CENTRE, "C4", y_safe_region); + EXPECT_TRUE(IsFieldEqual(actual, expected, y_safe_region)); +} From 15ab3c1ee10b3e3266e8daec2a0290e794b6a548 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 1 Sep 2026 07:38:23 -0700 Subject: [PATCH 17/18] DDY_stencil: Take Field3DParallel Alternative that takes Field3D and checks hasParallelSlices(). --- include/bout/stencil_expr.hxx | 39 +++++++++++++++++++++++++++++++---- 1 file changed, 35 insertions(+), 4 deletions(-) diff --git a/include/bout/stencil_expr.hxx b/include/bout/stencil_expr.hxx index 61ac66a5d8..b8dc6065cf 100644 --- a/include/bout/stencil_expr.hxx +++ b/include/bout/stencil_expr.hxx @@ -745,10 +745,12 @@ inline bout::stencil::DDXDispatchExpr DDX(const Field3D& f, CELL_LOC outloc, return DDX(f, outloc, parseField3DMethodString(method), region); } +namespace detail { + +template inline bout::stencil::DDYDispatchExpr -DDY_stencil(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT, - DIFF_METHOD method = DIFF_DEFAULT, - const std::string& region = "RGN_NOBNDRY") { +makeDDYStencilExpr(const FieldType& f, CELL_LOC outloc, DIFF_METHOD method, + const std::string& region) { checkData(f); if (!f.hasParallelSlices()) { @@ -799,8 +801,37 @@ DDY_stencil(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT, f.getMesh()->getRegion(region)}; } +} // namespace detail + +inline bout::stencil::DDYDispatchExpr +DDY_stencil(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT, + DIFF_METHOD method = DIFF_DEFAULT, + const std::string& region = "RGN_NOBNDRY") { + return detail::makeDDYStencilExpr(f, outloc, method, region); +} + +inline bout::stencil::DDYDispatchExpr +DDY_stencil(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT, + DIFF_METHOD method = DIFF_DEFAULT, + const std::string& region = "RGN_NOBNDRY") { + if (!f.hasParallelSlices()) { + throw BoutException("DDY_stencil requires parallel slices. Use eager DDY for " + "field-aligned transforms or communicate/apply parallel " + "boundaries before calling DDY_stencil."); + } + + return detail::makeDDYStencilExpr(f, outloc, method, region); +} + inline bout::stencil::DDYDispatchExpr -DDY_stencil(const Field3D& f, CELL_LOC outloc, const std::string& method, +DDY_stencil(const Field3DParallel& f, CELL_LOC outloc, const std::string& method, + const std::string& region = "RGN_NOBNDRY") { + return DDY_stencil(f, outloc, parseField3DMethodString(method), region); +} + +inline bout::stencil::DDYDispatchExpr +DDY_stencil(const Field3D& f, CELL_LOC outloc = CELL_DEFAULT, + const std::string& method = "DEFAULT", const std::string& region = "RGN_NOBNDRY") { return DDY_stencil(f, outloc, parseField3DMethodString(method), region); } From c43ebf8d88e2d6003029a71ca4396a902e17b780 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 1 Sep 2026 08:56:30 -0700 Subject: [PATCH 18/18] Petsc3DAMG: Set guards to zero if not finite Issue arises because only one guard cell is set in toField. --- src/invert/laplace/impls/petsc3damg/petsc3damg.cxx | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/invert/laplace/impls/petsc3damg/petsc3damg.cxx b/src/invert/laplace/impls/petsc3damg/petsc3damg.cxx index d254ea4a68..5ee479a89e 100644 --- a/src/invert/laplace/impls/petsc3damg/petsc3damg.cxx +++ b/src/invert/laplace/impls/petsc3damg/petsc3damg.cxx @@ -213,6 +213,7 @@ Field3D LaplacePetsc3dAmg::solve(const Field3D& b_in, const Field3D& x0) { if (updateRequired) { updateMatrix3D(); } + PetscVector rhs(b_in, indexer); PetscVector guess(x0, indexer); @@ -245,8 +246,14 @@ Field3D LaplacePetsc3dAmg::solve(const Field3D& b_in, const Field3D& x0) { KSPConvergedReasons[reason], static_cast(reason)); } - // Create field from result + // Reconstruct the PETSc result, then patch any deeper guards the vector does not + // populate so shifted-metric slice generation never sees NaNs. Field3D solution = guess.toField(); + BOUT_FOR_SERIAL(i, solution.getRegion("RGN_GUARDS")) { + if (!std::isfinite(solution[i])) { + solution[i] = 0.0; + } + } localmesh->communicate(solution); if (solution.hasParallelSlices()) { BOUT_FOR(i, indexer->getRegionLowerY()) { solution.ydown()[i] = solution[i]; }