From 5c6e3a76ac9f0bd3ce34fef53fa83624ed90ba85 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Thu, 6 Nov 2025 15:09:29 +0000 Subject: [PATCH 01/38] Port `FV::Div_par_mod` from Hermes-3 --- include/bout/fv_ops.hxx | 251 +++++++++++++++++++++++++++- tests/MMS/spatial/fci/data/BOUT.inp | 1 + tests/MMS/spatial/fci/fci_mms.cxx | 5 + tests/MMS/spatial/fci/mms.py | 2 + tests/MMS/spatial/fci/runtest | 9 +- 5 files changed, 265 insertions(+), 3 deletions(-) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 0ec1fbe3ad..adba5f21d7 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -5,12 +5,16 @@ #ifndef BOUT_FV_OPS_H #define BOUT_FV_OPS_H -#include "bout/build_defines.hxx" +#include + +#include "bout/bout_types.hxx" #include "bout/field3d.hxx" #include "bout/globals.hxx" +#include "bout/mesh.hxx" +#include "bout/output_bout_types.hxx" // NOLINT(unused-includes) +#include "bout/region.hxx" #include "bout/utils.hxx" #include "bout/vector2d.hxx" -#include namespace FV { /*! @@ -524,5 +528,248 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { * X-Z Finite Volume diffusion operator */ Field3D Div_Perp_Lap(const Field3D& a, const Field3D& f, CELL_LOC outloc = CELL_DEFAULT); + +/// Finite volume parallel divergence +/// +/// NOTE: Modified version, applies limiter to velocity and field +/// Performs better (smaller overshoots) than Div_par +/// +/// Preserves the sum of f*J*dx*dy*dz over the domain +/// +/// @param[in] f_in The field being advected. +/// This will be reconstructed at cell faces +/// using the given CellEdges method +/// @param[in] v_in The advection velocity. +/// This will be interpolated to cell boundaries +/// using linear interpolation +/// @param[in] wave_speed_in Local maximum speed of all waves in the system at each +// point in space +/// @param[in] fixflux Fix the flux at the boundary to be the value at the +/// midpoint (for boundary conditions) +/// +/// @param[out] flow_ylow Flow at the lower Y cell boundary +/// Already includes area factor * flux +template +Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, Field3D& flow_ylow, + bool fixflux = true) { + + Coordinates* coord = f_in.getCoordinates(); + + if (f_in.isFci()) { + // Use mid-point (cell boundary) averages + if (flow_ylow.isAllocated()) { + flow_ylow = emptyFrom(flow_ylow); + } + + ASSERT1(f_in.hasParallelSlices()); + ASSERT1(v_in.hasParallelSlices()); + + const auto& f_up = f_in.yup(); + const auto& f_down = f_in.ydown(); + + const auto& v_up = v_in.yup(); + const auto& v_down = v_in.ydown(); + + Field3D result{emptyFrom(f_in)}; + BOUT_FOR(i, f_in.getRegion("RGN_NOBNDRY")) { + const auto iyp = i.yp(); + const auto iym = i.ym(); + + result[i] = (0.25 * (f_in[i] + f_up[iyp]) * (v_in[i] + v_up[iyp]) + * (coord->J[i] + coord->J.yup()[iyp]) + / (sqrt(coord->g_22[i]) + sqrt(coord->g_22.yup()[iyp])) + - 0.25 * (f_in[i] + f_down[iym]) * (v_in[i] + v_down[iym]) + * (coord->J[i] + coord->J.ydown()[iym]) + / (sqrt(coord->g_22[i]) + sqrt(coord->g_22.ydown()[iym]))) + / (coord->dy[i] * coord->J[i]); + } + return result; + } + ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); + ASSERT1_FIELDS_COMPATIBLE(f_in, wave_speed_in); + + const Mesh* mesh = f_in.getMesh(); + + CellEdges cellboundary; + + ASSERT2(f_in.getDirectionY() == v_in.getDirectionY()); + ASSERT2(f_in.getDirectionY() == wave_speed_in.getDirectionY()); + const bool are_unaligned = + ((f_in.getDirectionY() == YDirectionType::Standard) + and (v_in.getDirectionY() == YDirectionType::Standard) + and (wave_speed_in.getDirectionY() == YDirectionType::Standard)); + + const Field3D f = are_unaligned ? toFieldAligned(f_in, "RGN_NOX") : f_in; + const Field3D v = are_unaligned ? toFieldAligned(v_in, "RGN_NOX") : v_in; + const Field3D wave_speed = + are_unaligned ? toFieldAligned(wave_speed_in, "RGN_NOX") : wave_speed_in; + + Field3D result{zeroFrom(f)}; + flow_ylow = zeroFrom(f); + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + const bool is_periodic_y = mesh->periodicY(i); + const bool is_first_y = mesh->firstY(i); + const bool is_last_y = mesh->lastY(i); + + // Only need one guard cell, so no need to communicate fluxes Instead + // calculate in guard cells to get fluxes consistent between processors, but + // don't include the boundary cell. Note that this implies special handling + // of boundaries later + const int ys = (!is_first_y || is_periodic_y) ? mesh->ystart - 1 : mesh->ystart; + const int ye = (!is_last_y || is_periodic_y) ? mesh->yend + 1 : mesh->yend; + + for (int j = ys; j <= ye; j++) { + // Pre-calculate factors which multiply fluxes +#if not(BOUT_USE_METRIC_3D) + // For right cell boundaries + const BoutReal common_factor_r = + (coord->J(i, j) + coord->J(i, j + 1)) + / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1))); + + const BoutReal flux_factor_rc = + common_factor_r / (coord->dy(i, j) * coord->J(i, j)); + const BoutReal flux_factor_rp = + common_factor_r / (coord->dy(i, j + 1) * coord->J(i, j + 1)); + + const BoutReal area_rp = + common_factor_r * coord->dx(i, j + 1) * coord->dz(i, j + 1); + + // For left cell boundaries + const BoutReal common_factor_l = + (coord->J(i, j) + coord->J(i, j - 1)) + / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1))); + + const BoutReal flux_factor_lc = + common_factor_l / (coord->dy(i, j) * coord->J(i, j)); + const BoutReal flux_factor_lm = + common_factor_l / (coord->dy(i, j - 1) * coord->J(i, j - 1)); + + const BoutReal area_lc = common_factor_l * coord->dx(i, j) * coord->dz(i, j); +#endif + for (int k = 0; k < mesh->LocalNz; k++) { +#if BOUT_USE_METRIC_3D + // For right cell boundaries + const BoutReal common_factor_r = + (coord->J(i, j, k) + coord->J(i, j + 1, k)) + / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); + + const BoutReal flux_factor_rc = + common_factor_r / (coord->dy(i, j, k) * coord->J(i, j, k)); + const BoutReal flux_factor_rp = + common_factor_r / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); + + const BoutReal area_rp = + common_factor_r * coord->dx(i, j + 1, k) * coord->dz(i, j + 1, k); + + // For left cell boundaries + const BoutReal common_factor_l = + (coord->J(i, j, k) + coord->J(i, j - 1, k)) + / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); + + const BoutReal flux_factor_lc = + common_factor_l / (coord->dy(i, j, k) * coord->J(i, j, k)); + const BoutReal flux_factor_lm = + common_factor_l / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); + + const BoutReal area_lc = + common_factor_l * coord->dx(i, j, k) * coord->dz(i, j, k); +#endif + + //////////////////////////////////////////// + // Reconstruct f at the cell faces + // This calculates s.R and s.L for the Right and Left + // face values on this cell + + // Reconstruct f at the cell faces + // TODO(peter): We can remove this #ifdef guard after switching to C++20 +#if __cpp_designated_initializers >= 201707L + Stencil1D s{.c = f(i, j, k), .m = f(i, j - 1, k), .p = f(i, j + 1, k)}; +#else + Stencil1D s{f(i, j, k), f(i, j - 1, k), f(i, j + 1, k), BoutNaN, + BoutNaN, BoutNaN, BoutNaN}; +#endif + cellboundary(s); // Calculate s.R and s.L + + //////////////////////////////////////////// + // Reconstruct v at the cell faces + // TODO(peter): We can remove this #ifdef guard after switching to C++20 +#if __cpp_designated_initializers >= 201707L + Stencil1D sv{.c = v(i, j, k), .m = v(i, j - 1, k), .p = v(i, j + 1, k)}; +#else + Stencil1D sv{v(i, j, k), v(i, j - 1, k), v(i, j + 1, k), BoutNaN, + BoutNaN, BoutNaN, BoutNaN}; +#endif + cellboundary(sv); // Calculate sv.R and sv.L + + //////////////////////////////////////////// + // Right boundary + + BoutReal flux = BoutNaN; + + if (is_last_y && (j == mesh->yend) && !is_periodic_y) { + // Last point in domain + + // Calculate velocity at right boundary (y+1/2) + const BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j + 1, k)); + + const BoutReal bndryval = 0.5 * (s.c + s.p); + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = bndryval * vpar; + } else { + // Add flux due to difference in boundary values + flux = (s.R * vpar) + (wave_speed(i, j, k) * (s.R - bndryval)); + } + + } else { + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k), + fabs(v(i, j, k)), fabs(v(i, j + 1, k))); + + flux = s.R * 0.5 * (sv.R + amax); + } + + result(i, j, k) += flux * flux_factor_rc; + result(i, j + 1, k) -= flux * flux_factor_rp; + + flow_ylow(i, j + 1, k) += flux * area_rp; + + //////////////////////////////////////////// + // Calculate at left boundary + + if (is_first_y && (j == mesh->ystart) && !is_periodic_y) { + // First point in domain + const BoutReal bndryval = 0.5 * (s.c + s.m); + const BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j - 1, k)); + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = bndryval * vpar; + } else { + // Add flux due to difference in boundary values + flux = (s.L * vpar) - (wave_speed(i, j, k) * (s.L - bndryval)); + } + } else { + + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k), + fabs(v(i, j, k)), fabs(v(i, j - 1, k))); + + flux = s.L * 0.5 * (sv.L - amax); + } + + result(i, j, k) -= flux * flux_factor_lc; + result(i, j - 1, k) += flux * flux_factor_lm; + + flow_ylow(i, j, k) += flux * area_lc; + } + } + } + if (are_unaligned) { + flow_ylow = fromFieldAligned(flow_ylow, "RGN_NOBNDRY"); + } + return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result; +} } // namespace FV #endif // BOUT_FV_OPS_H diff --git a/tests/MMS/spatial/fci/data/BOUT.inp b/tests/MMS/spatial/fci/data/BOUT.inp index 93e2101473..9171178d5d 100644 --- a/tests/MMS/spatial/fci/data/BOUT.inp +++ b/tests/MMS/spatial/fci/data/BOUT.inp @@ -5,6 +5,7 @@ div_par_solution = (0.01*x + 0.045)*(-12.5663706143592*cos(y - 2*z) - 6.28318530 div_par_K_grad_par_solution = (0.01*x + 0.045)*(6.28318530717959*sin(y - z) - 0.628318530717959*sin(y - z)/(0.01*x + 0.045))*(6.28318530717959*(0.01*x + 0.045)*(-2*cos(y - 2*z) - cos(y - z)) + 0.628318530717959*cos(y - 2*z) + 0.628318530717959*cos(y - z))/((0.01*x + 0.045)^2 + 1.0) + (6.28318530717959*(0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0))*cos(y - z)/sqrt((0.01*x + 0.045)^2 + 1.0) K = cos(y - z) laplace_par_solution = (0.01*x + 0.045)*(6.28318530717959*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/((0.01*x + 0.045)*sqrt((0.01*x + 0.045)^2 + 1.0)))/sqrt((0.01*x + 0.045)^2 + 1.0) +FV_div_par_mod_solution = (0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*((sin(y - 2*z) + sin(y - z))*sin(y - z)/(0.01*x + 0.045) + (-2*cos(y - 2*z) - cos(y - z))*cos(y - z)/(0.01*x + 0.045)) - 0.628318530717959*(sin(y - 2*z) + sin(y - z))*sin(y - z)/(0.01*x + 0.045) + 0.628318530717959*(cos(y - 2*z) + cos(y - z))*cos(y - z)/(0.01*x + 0.045))/sqrt((0.01*x + 0.045)^2 + 1.0) [mesh] symmetricglobalx = true diff --git a/tests/MMS/spatial/fci/fci_mms.cxx b/tests/MMS/spatial/fci/fci_mms.cxx index 7967452f3d..ca8ee8cd9f 100644 --- a/tests/MMS/spatial/fci/fci_mms.cxx +++ b/tests/MMS/spatial/fci/fci_mms.cxx @@ -4,6 +4,7 @@ #include "bout/field.hxx" #include "bout/field3d.hxx" #include "bout/field_factory.hxx" +#include "bout/fv_ops.hxx" #include "bout/globals.hxx" #include "bout/options.hxx" #include "bout/options_io.hxx" @@ -64,6 +65,10 @@ int main(int argc, char** argv) { fci_op_test("div_par_K_grad_par", dump, input, Div_par_K_Grad_par(K, input)); fci_op_test("laplace_par", dump, input, Laplace_par(input)); + // Finite volume methods + Field3D flow_ylow; + fci_op_test("FV_div_par_mod", dump, input, FV::Div_par_mod(input, K, K, flow_ylow)); + bout::writeDefaultOutputFile(dump); BoutFinalise(); diff --git a/tests/MMS/spatial/fci/mms.py b/tests/MMS/spatial/fci/mms.py index b28e337ac0..2fb2bd6aa6 100755 --- a/tests/MMS/spatial/fci/mms.py +++ b/tests/MMS/spatial/fci/mms.py @@ -16,6 +16,7 @@ f = sin(y - z) + sin(y - 2 * z) K = cos(z - y) + Lx = 0.1 Ly = 10.0 Lz = 1.0 @@ -61,6 +62,7 @@ def FCI_Laplace_par(f: Expr) -> Expr: ("div_par_solution", FCI_div_par(f)), ("div_par_K_grad_par_solution", FCI_div_par_K_grad_par(f, K)), ("laplace_par_solution", FCI_Laplace_par(f)), + ("FV_div_par_mod_solution", FCI_div_par(f * K)), ): expr_str = exprToStr(expr) print(f"{name} = {expr_str}") diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index 34340e53f4..7e960cb024 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -23,7 +23,14 @@ from scipy.interpolate import RectBivariateSpline as RBS DIRECTORY = "data" NPROC = 2 MTHREAD = 2 -OPERATORS = ("grad_par", "grad2_par2", "div_par", "div_par_K_grad_par", "laplace_par") +OPERATORS = ( + "grad_par", + "grad2_par2", + "div_par", + "div_par_K_grad_par", + "laplace_par", + "FV_div_par_mod", +) # Note that we need at least _2_ interior points for hermite spline # interpolation due to an awkwardness with the boundaries NX = 4 From 3059d128d44aacca266587a460a5b010afd3ba00 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Thu, 6 Nov 2025 15:43:38 +0000 Subject: [PATCH 02/38] Port `FV::Div_par_fvv` from Hermes-3 --- include/bout/fv_ops.hxx | 208 ++++++++++++++++++++++++++++ tests/MMS/spatial/fci/data/BOUT.inp | 1 + tests/MMS/spatial/fci/fci_mms.cxx | 1 + tests/MMS/spatial/fci/mms.py | 1 + tests/MMS/spatial/fci/runtest | 1 + 5 files changed, 212 insertions(+) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index adba5f21d7..bd0fa812d8 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -771,5 +771,213 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, } return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result; } + +/// This operator calculates Div_par(f v v) +/// It is used primarily (only?) in the parallel momentum equation. +/// +/// This operator is used rather than Div(f fv) so that the values of +/// f and v are consistent with other advection equations: The product +/// fv is not interpolated to cell boundaries. +template +Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true) { + ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); + const Mesh* mesh = f_in.getMesh(); + const Coordinates* coord = f_in.getCoordinates(); + CellEdges cellboundary; + + if (f_in.isFci()) { + // FCI version, using yup/down fields + ASSERT1(f_in.hasParallelSlices()); + ASSERT1(v_in.hasParallelSlices()); + + const auto& B = coord->Bxy; + const auto& B_up = coord->Bxy.yup(); + const auto& B_down = coord->Bxy.ydown(); + + const auto& f_up = f_in.yup(); + const auto& f_down = f_in.ydown(); + + const auto& v_up = v_in.yup(); + const auto& v_down = v_in.ydown(); + + const auto& g_22 = coord->g_22; + const auto& dy = coord->dy; + + Field3D result{emptyFrom(f_in)}; + BOUT_FOR(i, f_in.getRegion("RGN_NOBNDRY")) { + const auto iyp = i.yp(); + const auto iym = i.ym(); + + // Maximum local wave speed + const BoutReal amax = + BOUTMAX(wave_speed_in[i], fabs(v_in[i]), fabs(v_up[iyp]), fabs(v_down[iym])); + + const BoutReal term = (f_up[iyp] * v_up[iyp] * v_up[iyp] / B_up[iyp]) + - (f_down[iym] * v_down[iym] * v_down[iym] / B_down[iym]); + + // Penalty terms. This implementation is very dissipative. + BoutReal penalty = + (amax * (f_in[i] * v_in[i] - f_up[iyp] * v_up[iyp]) / (B[i] + B_up[iyp])) + + (amax * (f_in[i] * v_in[i] - f_down[iym] * v_down[iym]) + / (B[i] + B_down[iym])); + + if (fabs(penalty) > fabs(term) and penalty * v_in[i] > 0) { + if (term * penalty > 0) { + penalty = term; + } else { + penalty = -term; + } + } + + result[i] = B[i] * (term + penalty) / (2 * dy[i] * sqrt(g_22[i])); + +#if CHECK > 0 + if (!std::isfinite(result[i])) { + throw BoutException("Non-finite value in Div_par_fvv at {}\n" + "fup {} vup {} fdown {} vdown {} amax {}\n", + "B {} Bup {} Bdown {} dy {} sqrt(g_22} {}", i, f_up[i], + v_up[i], f_down[i], v_down[i], amax, B[i], B_up[i], B_down[i], + dy[i], sqrt(g_22[i])); + } +#endif + } + return result; + } + + ASSERT1(areFieldsCompatible(f_in, wave_speed_in)); + + /// Ensure that f, v and wave_speed are field aligned + Field3D f = toFieldAligned(f_in, "RGN_NOX"); + Field3D v = toFieldAligned(v_in, "RGN_NOX"); + Field3D wave_speed = toFieldAligned(wave_speed_in, "RGN_NOX"); + + Field3D result{zeroFrom(f)}; + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + const bool is_periodic_y = mesh->periodicY(i); + const bool is_first_y = mesh->firstY(i); + const bool is_last_y = mesh->lastY(i); + + // Only need one guard cell, so no need to communicate fluxes Instead + // calculate in guard cells to get fluxes consistent between processors, but + // don't include the boundary cell. Note that this implies special handling + // of boundaries later + const int ys = (!is_first_y || is_periodic_y) ? mesh->ystart - 1 : mesh->ystart; + const int ye = (!is_last_y || is_periodic_y) ? mesh->yend + 1 : mesh->yend; + + for (int j = ys; j <= ye; j++) { + // Pre-calculate factors which multiply fluxes + + for (int k = 0; k < mesh->LocalNz; k++) { + // For right cell boundaries + const BoutReal common_factor_r = + (coord->J(i, j, k) + coord->J(i, j + 1, k)) + / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); + + const BoutReal flux_factor_rc = + common_factor_r / (coord->dy(i, j, k) * coord->J(i, j, k)); + const BoutReal flux_factor_rp = + common_factor_r / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); + + // For left cell boundaries + const BoutReal common_factor_l = + (coord->J(i, j, k) + coord->J(i, j - 1, k)) + / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); + + const BoutReal flux_factor_lc = + common_factor_l / (coord->dy(i, j, k) * coord->J(i, j, k)); + const BoutReal flux_factor_lm = + common_factor_l / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); + + //////////////////////////////////////////// + // Reconstruct f at the cell faces + // This calculates s.R and s.L for the Right and Left + // face values on this cell + + // Reconstruct f at the cell faces +#if __cpp_designated_initializers >= 201707L + Stencil1D s{.c = f(i, j, k), .m = f(i, j - 1, k), .p = f(i, j + 1, k)}; +#else + Stencil1D s{f(i, j, k), f(i, j - 1, k), f(i, j + 1, k), BoutNaN, + BoutNaN, BoutNaN, BoutNaN}; +#endif + cellboundary(s); // Calculate s.R and s.L + + //////////////////////////////////////////// + // Reconstruct v at the cell faces + // TODO(peter): We can remove this #ifdef guard after switching to C++20 +#if __cpp_designated_initializers >= 201707L + Stencil1D sv{.c = v(i, j, k), .m = v(i, j - 1, k), .p = v(i, j + 1, k)}; +#else + Stencil1D sv{v(i, j, k), v(i, j - 1, k), v(i, j + 1, k), BoutNaN, + BoutNaN, BoutNaN, BoutNaN}; +#endif + cellboundary(sv); + + //////////////////////////////////////////// + // Right boundary + + // Calculate velocity at right boundary (y+1/2) + const BoutReal v_mid_r = 0.5 * (sv.c + sv.p); + // And mid-point density at right boundary + const BoutReal n_mid_r = 0.5 * (s.c + s.p); + BoutReal flux = NAN; + + if (mesh->lastY(i) && (j == mesh->yend) && !mesh->periodicY(i)) { + // Last point in domain + + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = n_mid_r * v_mid_r * v_mid_r; + } else { + // Add flux due to difference in boundary values + flux = (s.R * sv.R * sv.R) // Use right cell edge values + + (BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.p)) * n_mid_r + * (sv.R - v_mid_r)); // Damp differences in velocity, not flux + } + } else { + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k), + fabs(sv.c), fabs(sv.p)); + + flux = s.R * 0.5 * (sv.R + amax) * sv.R; + } + + result(i, j, k) += flux * flux_factor_rc; + result(i, j + 1, k) -= flux * flux_factor_rp; + + //////////////////////////////////////////// + // Calculate at left boundary + + const BoutReal v_mid_l = 0.5 * (sv.c + sv.m); + const BoutReal n_mid_l = 0.5 * (s.c + s.m); + + if (mesh->firstY(i) && (j == mesh->ystart) && !mesh->periodicY(i)) { + // First point in domain + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = n_mid_l * v_mid_l * v_mid_l; + } else { + // Add flux due to difference in boundary values + flux = (s.L * sv.L * sv.L) + - (BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.m)) * n_mid_l + * (sv.L - v_mid_l)); + } + } else { + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k), + fabs(sv.c), fabs(sv.m)); + + flux = s.L * 0.5 * (sv.L - amax) * sv.L; + } + + result(i, j, k) -= flux * flux_factor_lc; + result(i, j - 1, k) += flux * flux_factor_lm; + } + } + } + return fromFieldAligned(result, "RGN_NOBNDRY"); +} } // namespace FV #endif // BOUT_FV_OPS_H diff --git a/tests/MMS/spatial/fci/data/BOUT.inp b/tests/MMS/spatial/fci/data/BOUT.inp index 9171178d5d..99edee6ecc 100644 --- a/tests/MMS/spatial/fci/data/BOUT.inp +++ b/tests/MMS/spatial/fci/data/BOUT.inp @@ -6,6 +6,7 @@ div_par_K_grad_par_solution = (0.01*x + 0.045)*(6.28318530717959*sin(y - z) - 0. K = cos(y - z) laplace_par_solution = (0.01*x + 0.045)*(6.28318530717959*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/((0.01*x + 0.045)*sqrt((0.01*x + 0.045)^2 + 1.0)))/sqrt((0.01*x + 0.045)^2 + 1.0) FV_div_par_mod_solution = (0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*((sin(y - 2*z) + sin(y - z))*sin(y - z)/(0.01*x + 0.045) + (-2*cos(y - 2*z) - cos(y - z))*cos(y - z)/(0.01*x + 0.045)) - 0.628318530717959*(sin(y - 2*z) + sin(y - z))*sin(y - z)/(0.01*x + 0.045) + 0.628318530717959*(cos(y - 2*z) + cos(y - z))*cos(y - z)/(0.01*x + 0.045))/sqrt((0.01*x + 0.045)^2 + 1.0) +FV_div_par_fvv_solution = (0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(2*(sin(y - 2*z) + sin(y - z))*sin(y - z)*cos(y - z)/(0.01*x + 0.045) + (-2*cos(y - 2*z) - cos(y - z))*cos(y - z)^2/(0.01*x + 0.045)) - 1.25663706143592*(sin(y - 2*z) + sin(y - z))*sin(y - z)*cos(y - z)/(0.01*x + 0.045) + 0.628318530717959*(cos(y - 2*z) + cos(y - z))*cos(y - z)^2/(0.01*x + 0.045))/sqrt((0.01*x + 0.045)^2 + 1.0) [mesh] symmetricglobalx = true diff --git a/tests/MMS/spatial/fci/fci_mms.cxx b/tests/MMS/spatial/fci/fci_mms.cxx index ca8ee8cd9f..13744c4965 100644 --- a/tests/MMS/spatial/fci/fci_mms.cxx +++ b/tests/MMS/spatial/fci/fci_mms.cxx @@ -68,6 +68,7 @@ int main(int argc, char** argv) { // Finite volume methods Field3D flow_ylow; fci_op_test("FV_div_par_mod", dump, input, FV::Div_par_mod(input, K, K, flow_ylow)); + fci_op_test("FV_div_par_fvv", dump, input, FV::Div_par_fvv(input, K, K)); bout::writeDefaultOutputFile(dump); diff --git a/tests/MMS/spatial/fci/mms.py b/tests/MMS/spatial/fci/mms.py index 2fb2bd6aa6..62089c7e21 100755 --- a/tests/MMS/spatial/fci/mms.py +++ b/tests/MMS/spatial/fci/mms.py @@ -63,6 +63,7 @@ def FCI_Laplace_par(f: Expr) -> Expr: ("div_par_K_grad_par_solution", FCI_div_par_K_grad_par(f, K)), ("laplace_par_solution", FCI_Laplace_par(f)), ("FV_div_par_mod_solution", FCI_div_par(f * K)), + ("FV_div_par_fvv_solution", FCI_div_par(f * K * K)), ): expr_str = exprToStr(expr) print(f"{name} = {expr_str}") diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index 7e960cb024..c0f0a45132 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -30,6 +30,7 @@ OPERATORS = ( "div_par_K_grad_par", "laplace_par", "FV_div_par_mod", + "FV_div_par_fvv", ) # Note that we need at least _2_ interior points for hermite spline # interpolation due to an awkwardness with the boundaries From 60f194643a0819e0c7a9ac0dff4013968dac80e2 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Thu, 6 Nov 2025 16:12:41 +0000 Subject: [PATCH 03/38] Port `Div_par_K_Grad_par_mod` from Hermes-3 --- include/bout/difops.hxx | 4 ++ src/mesh/difops.cxx | 102 ++++++++++++++++++++++++++++ tests/MMS/spatial/fci/data/BOUT.inp | 1 + tests/MMS/spatial/fci/fci_mms.cxx | 6 +- tests/MMS/spatial/fci/mms.py | 1 + tests/MMS/spatial/fci/runtest | 1 + 6 files changed, 114 insertions(+), 1 deletion(-) diff --git a/include/bout/difops.hxx b/include/bout/difops.hxx index 18220b63ad..2cd99f8d33 100644 --- a/include/bout/difops.hxx +++ b/include/bout/difops.hxx @@ -195,6 +195,10 @@ Field3D Div_par_K_Grad_par(const Field3D& kY, const Field2D& f, Field3D Div_par_K_Grad_par(const Field3D& kY, const Field3D& f, CELL_LOC outloc = CELL_DEFAULT); +/// Version with energy flow diagnostic +Field3D Div_par_K_Grad_par_mod(const Field3D& k, const Field3D& f, Field3D& flow_ylow, + bool bndry_flux = true); + /*! * Perpendicular Laplacian operator * diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 09433b0685..09d2441951 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -366,6 +366,108 @@ Field3D Div_par_K_Grad_par(const Field3D& kY, const Field3D& f, CELL_LOC outloc) + Div_par(kY, outloc) * Grad_par(f, outloc); } +Field3D Div_par_K_Grad_par_mod(const Field3D& Kin, const Field3D& fin, + Field3D& flow_ylow, bool bndry_flux) { + TRACE("FV::Div_par_K_Grad_par_mod"); + + ASSERT2(Kin.getLocation() == fin.getLocation()); + + Mesh* mesh = Kin.getMesh(); + Coordinates* coord = fin.getCoordinates(); + + if (Kin.hasParallelSlices() && fin.hasParallelSlices()) { + // Using parallel slices. + // Note: Y slices may use different coordinate systems + // -> Only B, dy and g_22 can be used in yup/ydown + // Others (e.g J) may not be averaged between y planes. + + const auto& K_up = Kin.yup(); + const auto& K_down = Kin.ydown(); + + const auto& f_up = fin.yup(); + const auto& f_down = fin.ydown(); + + Field3D result{zeroFrom(fin)}; + flow_ylow = zeroFrom(fin); + + BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { + const auto iyp = i.yp(); + const auto iym = i.ym(); + + // Upper cell edge + const BoutReal c_up = 0.5 * (Kin[i] + K_up[iyp]); // K at the upper boundary + const BoutReal J_up = 0.5 * (coord->J[i] + coord->J.yup()[iyp]); // Jacobian at boundary + const BoutReal g_22_up = 0.5 * (coord->g_22[i] + coord->g_22.yup()[iyp]); + const BoutReal gradient_up = 2. * (f_up[iyp] - fin[i]) / (coord->dy[i] + coord->dy.yup()[iyp]); + + const BoutReal flux_up = c_up * J_up * gradient_up / g_22_up; + + // Lower cell edge + const BoutReal c_down = 0.5 * (Kin[i] + K_down[iym]); // K at the lower boundary + const BoutReal J_down = 0.5 * (coord->J[i] + coord->J.ydown()[iym]); // Jacobian at boundary + const BoutReal g_22_down = 0.5 * (coord->g_22[i] + coord->g_22.ydown()[iym]); + const BoutReal gradient_down = 2. * (fin[i] - f_down[iym]) / (coord->dy[i] + coord->dy.ydown()[iym]); + + const BoutReal flux_down = c_down * J_down * gradient_down / g_22_down; + + result[i] = (flux_up - flux_down) / (coord->dy[i] * coord->J[i]); + } + + return result; + } + + // Calculate in field-aligned coordinates + const auto& K = toFieldAligned(Kin, "RGN_NOX"); + const auto& f = toFieldAligned(fin, "RGN_NOX"); + + Field3D result{zeroFrom(f)}; + flow_ylow = zeroFrom(f); + + BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { + // Calculate flux at upper surface + + const auto iyp = i.yp(); + const auto iym = i.ym(); + + if (bndry_flux || mesh->periodicY(i.x()) || !mesh->lastY(i.x()) + || (i.y() != mesh->yend)) { + + BoutReal c = 0.5 * (K[i] + K[iyp]); // K at the upper boundary + BoutReal J = 0.5 * (coord->J[i] + coord->J[iyp]); // Jacobian at boundary + BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iyp]); + + BoutReal gradient = 2. * (f[iyp] - f[i]) / (coord->dy[i] + coord->dy[iyp]); + + BoutReal flux = c * J * gradient / g_22; + + result[i] += flux / (coord->dy[i] * coord->J[i]); + } + + // Calculate flux at lower surface + if (bndry_flux || mesh->periodicY(i.x()) || !mesh->firstY(i.x()) + || (i.y() != mesh->ystart)) { + BoutReal c = 0.5 * (K[i] + K[iym]); // K at the lower boundary + BoutReal J = 0.5 * (coord->J[i] + coord->J[iym]); // Jacobian at boundary + + BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iym]); + + BoutReal gradient = 2. * (f[i] - f[iym]) / (coord->dy[i] + coord->dy[iym]); + + BoutReal flux = c * J * gradient / g_22; + + result[i] -= flux / (coord->dy[i] * coord->J[i]); + flow_ylow[i] = -flux * coord->dx[i] * coord->dz[i]; + } + } + + // Shifted to field aligned coordinates, so need to shift back + result = fromFieldAligned(result, "RGN_NOBNDRY"); + flow_ylow = fromFieldAligned(flow_ylow); + + return result; +} + + /******************************************************************************* * Delp2 * perpendicular Laplacian operator diff --git a/tests/MMS/spatial/fci/data/BOUT.inp b/tests/MMS/spatial/fci/data/BOUT.inp index 99edee6ecc..76ac3035c9 100644 --- a/tests/MMS/spatial/fci/data/BOUT.inp +++ b/tests/MMS/spatial/fci/data/BOUT.inp @@ -7,6 +7,7 @@ K = cos(y - z) laplace_par_solution = (0.01*x + 0.045)*(6.28318530717959*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/((0.01*x + 0.045)*sqrt((0.01*x + 0.045)^2 + 1.0)))/sqrt((0.01*x + 0.045)^2 + 1.0) FV_div_par_mod_solution = (0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*((sin(y - 2*z) + sin(y - z))*sin(y - z)/(0.01*x + 0.045) + (-2*cos(y - 2*z) - cos(y - z))*cos(y - z)/(0.01*x + 0.045)) - 0.628318530717959*(sin(y - 2*z) + sin(y - z))*sin(y - z)/(0.01*x + 0.045) + 0.628318530717959*(cos(y - 2*z) + cos(y - z))*cos(y - z)/(0.01*x + 0.045))/sqrt((0.01*x + 0.045)^2 + 1.0) FV_div_par_fvv_solution = (0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(2*(sin(y - 2*z) + sin(y - z))*sin(y - z)*cos(y - z)/(0.01*x + 0.045) + (-2*cos(y - 2*z) - cos(y - z))*cos(y - z)^2/(0.01*x + 0.045)) - 1.25663706143592*(sin(y - 2*z) + sin(y - z))*sin(y - z)*cos(y - z)/(0.01*x + 0.045) + 0.628318530717959*(cos(y - 2*z) + cos(y - z))*cos(y - z)^2/(0.01*x + 0.045))/sqrt((0.01*x + 0.045)^2 + 1.0) +div_par_K_grad_par_mod_solution = (0.01*x + 0.045)*(6.28318530717959*sin(y - z) - 0.628318530717959*sin(y - z)/(0.01*x + 0.045))*(6.28318530717959*(0.01*x + 0.045)*(-2*cos(y - 2*z) - cos(y - z)) + 0.628318530717959*cos(y - 2*z) + 0.628318530717959*cos(y - z))/((0.01*x + 0.045)^2 + 1.0) + (6.28318530717959*(0.01*x + 0.045)*(6.28318530717959*(0.01*x + 0.045)*(-4*sin(y - 2*z) - sin(y - z)) + 1.25663706143592*sin(y - 2*z) + 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0) + 0.628318530717959*(6.28318530717959*(0.01*x + 0.045)*(2*sin(y - 2*z) + sin(y - z)) - 0.628318530717959*sin(y - 2*z) - 0.628318530717959*sin(y - z))/sqrt((0.01*x + 0.045)^2 + 1.0))*cos(y - z)/sqrt((0.01*x + 0.045)^2 + 1.0) [mesh] symmetricglobalx = true diff --git a/tests/MMS/spatial/fci/fci_mms.cxx b/tests/MMS/spatial/fci/fci_mms.cxx index 13744c4965..17408aeeab 100644 --- a/tests/MMS/spatial/fci/fci_mms.cxx +++ b/tests/MMS/spatial/fci/fci_mms.cxx @@ -59,14 +59,18 @@ int main(int argc, char** argv) { // Add mesh geometry variables mesh->outputVars(dump); + // Dummy variable for *_mod overloads + Field3D flow_ylow; + fci_op_test("grad_par", dump, input, Grad_par(input)); fci_op_test("grad2_par2", dump, input, Grad2_par2(input)); fci_op_test("div_par", dump, input, Div_par(input)); fci_op_test("div_par_K_grad_par", dump, input, Div_par_K_Grad_par(K, input)); + fci_op_test("div_par_K_grad_par_mod", dump, input, + Div_par_K_Grad_par_mod(K, input, flow_ylow)); fci_op_test("laplace_par", dump, input, Laplace_par(input)); // Finite volume methods - Field3D flow_ylow; fci_op_test("FV_div_par_mod", dump, input, FV::Div_par_mod(input, K, K, flow_ylow)); fci_op_test("FV_div_par_fvv", dump, input, FV::Div_par_fvv(input, K, K)); diff --git a/tests/MMS/spatial/fci/mms.py b/tests/MMS/spatial/fci/mms.py index 62089c7e21..801a8d3f26 100755 --- a/tests/MMS/spatial/fci/mms.py +++ b/tests/MMS/spatial/fci/mms.py @@ -61,6 +61,7 @@ def FCI_Laplace_par(f: Expr) -> Expr: ("grad2_par2_solution", FCI_grad2_par2(f)), ("div_par_solution", FCI_div_par(f)), ("div_par_K_grad_par_solution", FCI_div_par_K_grad_par(f, K)), + ("div_par_K_grad_par_mod_solution", FCI_div_par_K_grad_par(f, K)), ("laplace_par_solution", FCI_Laplace_par(f)), ("FV_div_par_mod_solution", FCI_div_par(f * K)), ("FV_div_par_fvv_solution", FCI_div_par(f * K * K)), diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index c0f0a45132..73babc9691 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -28,6 +28,7 @@ OPERATORS = ( "grad2_par2", "div_par", "div_par_K_grad_par", + "div_par_K_grad_par_mod", "laplace_par", "FV_div_par_mod", "FV_div_par_fvv", From f53b8f3c5ec7d6af148200e8a2af4ab1ec80aac0 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Thu, 6 Nov 2025 16:59:34 +0000 Subject: [PATCH 04/38] Add MMS tests for finite volume operators --- tests/MMS/CMakeLists.txt | 1 + .../MMS/spatial/finite-volume/CMakeLists.txt | 6 + tests/MMS/spatial/finite-volume/data/BOUT.inp | 23 ++ tests/MMS/spatial/finite-volume/fv_mms.cxx | 61 +++++ tests/MMS/spatial/finite-volume/makefile | 6 + tests/MMS/spatial/finite-volume/mms.py | 62 ++++++ tests/MMS/spatial/finite-volume/runtest | 209 ++++++++++++++++++ 7 files changed, 368 insertions(+) create mode 100644 tests/MMS/spatial/finite-volume/CMakeLists.txt create mode 100644 tests/MMS/spatial/finite-volume/data/BOUT.inp create mode 100644 tests/MMS/spatial/finite-volume/fv_mms.cxx create mode 100644 tests/MMS/spatial/finite-volume/makefile create mode 100755 tests/MMS/spatial/finite-volume/mms.py create mode 100755 tests/MMS/spatial/finite-volume/runtest diff --git a/tests/MMS/CMakeLists.txt b/tests/MMS/CMakeLists.txt index a6667bcfa5..510385d54c 100644 --- a/tests/MMS/CMakeLists.txt +++ b/tests/MMS/CMakeLists.txt @@ -8,6 +8,7 @@ add_subdirectory(spatial/d2dx2) add_subdirectory(spatial/d2dz2) add_subdirectory(spatial/diffusion) add_subdirectory(spatial/fci) +add_subdirectory(spatial/finite-volume) add_subdirectory(time) add_subdirectory(time-petsc) add_subdirectory(wave-1d) diff --git a/tests/MMS/spatial/finite-volume/CMakeLists.txt b/tests/MMS/spatial/finite-volume/CMakeLists.txt new file mode 100644 index 0000000000..6d9c839a05 --- /dev/null +++ b/tests/MMS/spatial/finite-volume/CMakeLists.txt @@ -0,0 +1,6 @@ +bout_add_mms_test(MMS-spatial-finite-volume + SOURCES fv_mms.cxx + USE_RUNTEST + USE_DATA_BOUT_INP + PROCESSORS 2 +) diff --git a/tests/MMS/spatial/finite-volume/data/BOUT.inp b/tests/MMS/spatial/finite-volume/data/BOUT.inp new file mode 100644 index 0000000000..afab4a34d5 --- /dev/null +++ b/tests/MMS/spatial/finite-volume/data/BOUT.inp @@ -0,0 +1,23 @@ +input_field = 0.1*sin(2.0*y) + 1 +K = 0.1*cos(3.0*y) + 1 +FV_Div_par_mod_solution = -0.188495559215388*(0.1*sin(2.0*y) + 1)*sin(3.0*y) + 0.125663706143592*(0.1*cos(3.0*y) + 1)*cos(2.0*y) +FV_Div_par_fvv_solution = -0.376991118430775*(0.1*sin(2.0*y) + 1)*(0.1*cos(3.0*y) + 1)*sin(3.0*y) + 0.125663706143592*(0.1*cos(3.0*y) + 1)^2*cos(2.0*y) +FV_Div_par_solution = -0.188495559215388*(0.1*sin(2.0*y) + 1)*sin(3.0*y) + 0.125663706143592*(0.1*cos(3.0*y) + 1)*cos(2.0*y) +FV_Div_par_K_Grad_par_solution = -0.15791367041743*(0.1*cos(3.0*y) + 1)*sin(2.0*y) - 0.0236870505626145*sin(3.0*y)*cos(2.0*y) +FV_Div_par_K_Grad_par_mod_solution = -0.15791367041743*(0.1*cos(3.0*y) + 1)*sin(2.0*y) - 0.0236870505626145*sin(3.0*y)*cos(2.0*y) + +[mesh] +MXG = 0 + +nx = 1 +ny = 128 +nz = 1 + +Ly = 10 + +dy = Ly / ny +J = 1 # Identity metric + +[mesh:ddy] +first = C2 +second = C2 diff --git a/tests/MMS/spatial/finite-volume/fv_mms.cxx b/tests/MMS/spatial/finite-volume/fv_mms.cxx new file mode 100644 index 0000000000..6b45ef3259 --- /dev/null +++ b/tests/MMS/spatial/finite-volume/fv_mms.cxx @@ -0,0 +1,61 @@ +#include "bout/bout.hxx" +#include "bout/field.hxx" +#include "bout/field3d.hxx" +#include "bout/field_factory.hxx" +#include "bout/fv_ops.hxx" +#include "bout/globals.hxx" +#include "bout/options.hxx" +#include "bout/options_io.hxx" +#include "bout/utils.hxx" + +#include + +#include +#include + +namespace { +auto fv_op_test(const std::string& name, Options& dump, const Field3D& input, + const Field3D& result) { + auto* mesh = input.getMesh(); + const Field3D solution{FieldFactory::get()->create3D(fmt::format("{}_solution", name), + Options::getRoot(), mesh)}; + const Field3D error{result - solution}; + + dump[fmt::format("{}_l_2", name)] = sqrt(mean(SQ(error), true, "RGN_NOBNDRY")); + dump[fmt::format("{}_l_inf", name)] = max(abs(error), true, "RGN_NOBNDRY"); + + dump[fmt::format("{}_result", name)] = result; + dump[fmt::format("{}_error", name)] = error; + dump[fmt::format("{}_input", name)] = input; + dump[fmt::format("{}_solution", name)] = solution; +} +} // namespace + +int main(int argc, char** argv) { + BoutInitialise(argc, argv); + + using bout::globals::mesh; + + Field3D input{FieldFactory::get()->create3D("input_field", Options::getRoot(), mesh)}; + Field3D K{FieldFactory::get()->create3D("K", Options::getRoot(), mesh)}; + + // Communicate to calculate parallel transform. + mesh->communicate(input, K); + + Options dump; + // Add mesh geometry variables + mesh->outputVars(dump); + + // Dummy variable for *_mod overloads + Field3D flow_ylow; + + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, K, K)); + fv_op_test("FV_Div_par_mod", dump, input, FV::Div_par_mod(input, K, K, flow_ylow)); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, K, K)); + fv_op_test("FV_Div_par_K_Grad_par", dump, input, FV::Div_par_K_Grad_par(K, input)); + fv_op_test("FV_Div_par_K_Grad_par_mod", dump, input, FV::Div_par_K_Grad_par(K, input)); + + bout::writeDefaultOutputFile(dump); + + BoutFinalise(); +} diff --git a/tests/MMS/spatial/finite-volume/makefile b/tests/MMS/spatial/finite-volume/makefile new file mode 100644 index 0000000000..88ba6c77e7 --- /dev/null +++ b/tests/MMS/spatial/finite-volume/makefile @@ -0,0 +1,6 @@ + +BOUT_TOP = ../../../.. + +SOURCEC = fci_mms.cxx + +include $(BOUT_TOP)/make.config diff --git a/tests/MMS/spatial/finite-volume/mms.py b/tests/MMS/spatial/finite-volume/mms.py new file mode 100755 index 0000000000..c8b473138d --- /dev/null +++ b/tests/MMS/spatial/finite-volume/mms.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python3 +# +# Generate manufactured solution and sources for FCI test +# + +from math import pi +import warnings + +from boututils.boutwarnings import AlwaysWarning +from boutdata.data import BoutOptionsFile +from boutdata.mms import exprToStr, y, Grad_par, Div_par, Metric +from sympy import sin, cos, Expr + +warnings.simplefilter("ignore", AlwaysWarning) + +# Length of the y domain +Ly = 10.0 + +# Identity +metric = Metric() + +# Define solution in terms of input x,y,z +f = 1 + 0.1 * sin(2 * y) +K = 1 + 0.1 * cos(3 * y) + +# Turn solution into real x and z coordinates +replace = [(y, metric.y * 2 * pi / Ly)] + +f = f.subs(replace) +K = K.subs(replace) + +# Substitute back to get input y coordinates +replace = [ (metric.y, y*Ly/(2*pi) ) ] + + +def Grad2_par2(f: Expr) -> Expr: + return Grad_par(Grad_par(f)) + + +def Div_par_K_Grad_par(f: Expr, K: Expr) -> Expr: + return (K * Grad2_par2(f)) + (Div_par(K) * Grad_par(f)) + + +############################################ +# Equations solved + +options = BoutOptionsFile("data/BOUT.inp") + +for name, expr in ( + ("input_field", f), + ("K", K), + ("FV_Div_par_solution", Div_par(f * K)), + ("FV_Div_par_K_Grad_par_solution", Div_par_K_Grad_par(f, K)), + ("FV_Div_par_K_Grad_par_mod_solution", Div_par_K_Grad_par(f, K)), + ("FV_Div_par_mod_solution", Div_par(f * K)), + ("FV_Div_par_fvv_solution", Div_par(f * K * K)), +): + expr_str = exprToStr(expr.subs(replace)) + print(f"{name} = {expr_str}") + options[name] = expr_str + +options.write("data/BOUT.inp", overwrite=True) diff --git a/tests/MMS/spatial/finite-volume/runtest b/tests/MMS/spatial/finite-volume/runtest new file mode 100755 index 0000000000..a836f5e735 --- /dev/null +++ b/tests/MMS/spatial/finite-volume/runtest @@ -0,0 +1,209 @@ +#!/usr/bin/env python3 +# +# Python script to run and analyse MMS test +# + +import argparse +import json +import pathlib +import sys +from time import time + +from boutdata.collect import collect +from boututils.run_wrapper import build_and_log, launch_safe +from numpy import array, log, polyfit + +# Global parameters +DIRECTORY = "data" +NPROC = 2 +MTHREAD = 2 +OPERATORS = ( + "FV_Div_par", + "FV_Div_par_K_Grad_par", + "FV_Div_par_K_Grad_par_mod", + "FV_Div_par_mod", + "FV_Div_par_fvv", +) +# Resolution in y and z +NLIST = [8, 16, 32, 64] +dx = 1.0 / array(NLIST) + + +def quiet_collect(name: str) -> float: + # Index to return a plain (numpy) float rather than `BoutArray` + return collect( + name, + tind=[1, 1], + info=False, + path=DIRECTORY, + xguards=False, + yguards=False, + )[()] + + +def assert_convergence(error, dx, name, expected) -> bool: + fit = polyfit(log(dx), log(error), 1) + order = fit[0] + print(f"{name} convergence order = {order:f} (fit)", end="") + + order = log(error[-2] / error[-1]) / log(dx[-2] / dx[-1]) + print(f", {order:f} (small spacing)", end="") + + # Should be close to the expected order + success = order > expected * 0.95 + print(f"\t............ {'PASS' if success else 'FAIL'}") + + return success + + +def run_fv_operators(nz: int, name: str) -> dict[str, float]: + # Command to run + args = f"MZ={nz} mesh:ny={nz} {name}" + cmd = f"./fv_mms {args}" + print(f"Running command: {cmd}", end="") + + # Launch using MPI + start = time() + status, out = launch_safe(cmd, nproc=NPROC, mthread=MTHREAD, pipe=True) + print(f" ... done in {time() - start:.3}s") + + # Save output to log file + pathlib.Path(f"run.log.{nz}").write_text(out) + + if status: + print(f"Run failed!\nOutput was:\n{out}") + sys.exit(status) + + return { + operator: { + "l_2": quiet_collect(f"{operator}_l_2"), + "l_inf": quiet_collect(f"{operator}_l_inf"), + } + for operator in OPERATORS + } + + +def transpose( + errors: list[dict[str, dict[str, float]]], +) -> dict[str, dict[str, list[float]]]: + """Turn a list of {operator: error} into a dict of {operator: [errors]}""" + + kinds = ("l_2", "l_inf") + result = {operator: {kind: [] for kind in kinds} for operator in OPERATORS} + for error in errors: + for k, v in error.items(): + for kind in kinds: + result[k][kind].append(v[kind]) + return result + + +def check_fv_operators(name: str, case: dict) -> bool: + failures = [] + + order = case["order"] + args = case["args"] + + all_errors = [] + + for n in NLIST: + errors = run_fv_operators(n, args) + all_errors.append(errors) + + for operator in OPERATORS: + l_2 = errors[operator]["l_2"] + l_inf = errors[operator]["l_inf"] + + print(f"{operator} errors: l-2 {l_2:f} l-inf {l_inf:f}") + + final_errors = transpose(all_errors) + for operator in OPERATORS: + test_name = f"{operator} {name}" + success = assert_convergence( + final_errors[operator]["l_2"], dx, test_name, order + ) + if not success: + failures.append(test_name) + + return final_errors, failures + + +def make_plots(cases: dict[str, dict]): + try: + import matplotlib.pyplot as plt + except ImportError: + print("No matplotlib") + return + + num_operators = len(OPERATORS) + fig, axes = plt.subplots(1, num_operators, figsize=(num_operators * 4, 4)) + + for ax, operator in zip(axes, OPERATORS): + for name, case in cases.items(): + ax.loglog(dx, case[operator]["l_2"], "-", label=f"{name} $l_2$") + ax.loglog(dx, case[operator]["l_inf"], "--", label=f"{name} $l_\\inf$") + ax.legend(loc="upper left") + ax.grid() + ax.set_title(f"Error scaling for {operator}") + ax.set_xlabel(r"Mesh spacing $\delta x$") + ax.set_ylabel("Error norm") + + fig.tight_layout() + fig.savefig("fv_mms.pdf") + print("Plot saved to fv_mms.pdf") + + if args.show_plots: + plt.show() + plt.close() + + +if __name__ == "__main__": + build_and_log("Finite volume MMS test") + + parser = argparse.ArgumentParser("Error scaling test for finite volume operators") + parser.add_argument( + "--make-plots", action="store_true", help="Create plots of error scaling" + ) + parser.add_argument( + "--show-plots", + action="store_true", + help="Stop and show plots, implies --make-plots", + ) + parser.add_argument( + "--dump-errors", + type=str, + help="Output file to dump errors as JSON", + default="fv_operator_errors.json", + ) + + args = parser.parse_args() + + success = True + failures = [] + + cases = { + "default": { + "order": 2, + "args": "", + }, + } + + for name, case in cases.items(): + error2, failures_ = check_fv_operators(name, case) + case.update(error2) + failures.extend(failures_) + success &= len(failures) == 0 + + if args.dump_errors: + pathlib.Path(args.dump_errors).write_text(json.dumps(cases)) + + if args.make_plots or args.show_plots: + make_plots(cases) + + if success: + print("\nAll tests passed") + else: + print("\nSome tests failed:") + for failure in failures: + print(f"\t{failure}") + + sys.exit(0 if success else 1) From 207802fdf4e9c2d328bef8d86f3929b5f0f9ce06 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Fri, 7 Nov 2025 14:20:14 +0000 Subject: [PATCH 05/38] Add tests for all FV slope limiters --- tests/MMS/spatial/finite-volume/data/BOUT.inp | 12 +-- tests/MMS/spatial/finite-volume/fv_mms.cxx | 51 +++++++++---- tests/MMS/spatial/finite-volume/mms.py | 18 +++-- tests/MMS/spatial/finite-volume/runtest | 73 +++++++++---------- 4 files changed, 87 insertions(+), 67 deletions(-) diff --git a/tests/MMS/spatial/finite-volume/data/BOUT.inp b/tests/MMS/spatial/finite-volume/data/BOUT.inp index afab4a34d5..029011e437 100644 --- a/tests/MMS/spatial/finite-volume/data/BOUT.inp +++ b/tests/MMS/spatial/finite-volume/data/BOUT.inp @@ -1,10 +1,10 @@ input_field = 0.1*sin(2.0*y) + 1 -K = 0.1*cos(3.0*y) + 1 -FV_Div_par_mod_solution = -0.188495559215388*(0.1*sin(2.0*y) + 1)*sin(3.0*y) + 0.125663706143592*(0.1*cos(3.0*y) + 1)*cos(2.0*y) -FV_Div_par_fvv_solution = -0.376991118430775*(0.1*sin(2.0*y) + 1)*(0.1*cos(3.0*y) + 1)*sin(3.0*y) + 0.125663706143592*(0.1*cos(3.0*y) + 1)^2*cos(2.0*y) -FV_Div_par_solution = -0.188495559215388*(0.1*sin(2.0*y) + 1)*sin(3.0*y) + 0.125663706143592*(0.1*cos(3.0*y) + 1)*cos(2.0*y) -FV_Div_par_K_Grad_par_solution = -0.15791367041743*(0.1*cos(3.0*y) + 1)*sin(2.0*y) - 0.0236870505626145*sin(3.0*y)*cos(2.0*y) -FV_Div_par_K_Grad_par_mod_solution = -0.15791367041743*(0.1*cos(3.0*y) + 1)*sin(2.0*y) - 0.0236870505626145*sin(3.0*y)*cos(2.0*y) +FV_Div_par_mod_solution = -0.188495559215388*sin(3.0*y) +FV_Div_par_fvv_solution = -0.376991118430775*(0.1*cos(3.0*y) + 1)*sin(3.0*y)/(0.1*sin(2.0*y) + 1) - 0.125663706143592*(0.1*cos(3.0*y) + 1)^2*cos(2.0*y)/(0.1*sin(2.0*y) + 1)^2 +FV_Div_par_solution = -0.188495559215388*sin(3.0*y) +FV_Div_par_K_Grad_par_solution = 0.125663706143592*(-0.188495559215388*sin(3.0*y)/(0.1*sin(2.0*y) + 1) - 0.125663706143592*(0.1*cos(3.0*y) + 1)*cos(2.0*y)/(0.1*sin(2.0*y) + 1)^2)*cos(2.0*y) - 0.15791367041743*(0.1*cos(3.0*y) + 1)*sin(2.0*y)/(0.1*sin(2.0*y) + 1) +FV_Div_par_K_Grad_par_mod_solution = 0.125663706143592*(-0.188495559215388*sin(3.0*y)/(0.1*sin(2.0*y) + 1) - 0.125663706143592*(0.1*cos(3.0*y) + 1)*cos(2.0*y)/(0.1*sin(2.0*y) + 1)^2)*cos(2.0*y) - 0.15791367041743*(0.1*cos(3.0*y) + 1)*sin(2.0*y)/(0.1*sin(2.0*y) + 1) +v = (0.1*cos(3.0*y) + 1)/(0.1*sin(2.0*y) + 1) [mesh] MXG = 0 diff --git a/tests/MMS/spatial/finite-volume/fv_mms.cxx b/tests/MMS/spatial/finite-volume/fv_mms.cxx index 6b45ef3259..09f9986e72 100644 --- a/tests/MMS/spatial/finite-volume/fv_mms.cxx +++ b/tests/MMS/spatial/finite-volume/fv_mms.cxx @@ -15,19 +15,20 @@ namespace { auto fv_op_test(const std::string& name, Options& dump, const Field3D& input, - const Field3D& result) { + const Field3D& result, std::string suffix = "") { auto* mesh = input.getMesh(); const Field3D solution{FieldFactory::get()->create3D(fmt::format("{}_solution", name), Options::getRoot(), mesh)}; const Field3D error{result - solution}; - dump[fmt::format("{}_l_2", name)] = sqrt(mean(SQ(error), true, "RGN_NOBNDRY")); - dump[fmt::format("{}_l_inf", name)] = max(abs(error), true, "RGN_NOBNDRY"); + dump[fmt::format("{}{}_l_2", name, suffix)] = + sqrt(mean(SQ(error), true, "RGN_NOBNDRY")); + dump[fmt::format("{}{}_l_inf", name, suffix)] = max(abs(error), true, "RGN_NOBNDRY"); - dump[fmt::format("{}_result", name)] = result; - dump[fmt::format("{}_error", name)] = error; - dump[fmt::format("{}_input", name)] = input; - dump[fmt::format("{}_solution", name)] = solution; + dump[fmt::format("{}{}_result", name, suffix)] = result; + dump[fmt::format("{}{}_error", name, suffix)] = error; + dump[fmt::format("{}{}_input", name, suffix)] = input; + dump[fmt::format("{}{}_solution", name, suffix)] = solution; } } // namespace @@ -37,23 +38,45 @@ int main(int argc, char** argv) { using bout::globals::mesh; Field3D input{FieldFactory::get()->create3D("input_field", Options::getRoot(), mesh)}; - Field3D K{FieldFactory::get()->create3D("K", Options::getRoot(), mesh)}; + Field3D v{FieldFactory::get()->create3D("v", Options::getRoot(), mesh)}; // Communicate to calculate parallel transform. - mesh->communicate(input, K); + mesh->communicate(input, v); Options dump; // Add mesh geometry variables mesh->outputVars(dump); + dump["v"] = v; // Dummy variable for *_mod overloads Field3D flow_ylow; - fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, K, K)); - fv_op_test("FV_Div_par_mod", dump, input, FV::Div_par_mod(input, K, K, flow_ylow)); - fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, K, K)); - fv_op_test("FV_Div_par_K_Grad_par", dump, input, FV::Div_par_K_Grad_par(K, input)); - fv_op_test("FV_Div_par_K_Grad_par_mod", dump, input, FV::Div_par_K_Grad_par(K, input)); + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), "_MC"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), "_MC"); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), "_MC"); + + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), "_Upwind"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), "_Upwind"); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), + "_Upwind"); + + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), "_Fromm"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), "_Fromm"); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), + "_Fromm"); + + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), "_MinMod"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), "_MinMod"); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), + "_MinMod"); + + fv_op_test("FV_Div_par_K_Grad_par", dump, input, FV::Div_par_K_Grad_par(v, input)); + fv_op_test("FV_Div_par_K_Grad_par_mod", dump, input, + Div_par_K_Grad_par_mod(v, input, flow_ylow)); bout::writeDefaultOutputFile(dump); diff --git a/tests/MMS/spatial/finite-volume/mms.py b/tests/MMS/spatial/finite-volume/mms.py index c8b473138d..dfcfce9a09 100755 --- a/tests/MMS/spatial/finite-volume/mms.py +++ b/tests/MMS/spatial/finite-volume/mms.py @@ -21,13 +21,15 @@ # Define solution in terms of input x,y,z f = 1 + 0.1 * sin(2 * y) -K = 1 + 0.1 * cos(3 * y) +fv = 1 + 0.1 * cos(3 * y) + # Turn solution into real x and z coordinates replace = [(y, metric.y * 2 * pi / Ly)] f = f.subs(replace) -K = K.subs(replace) +fv = fv.subs(replace) +v = fv / f # Substitute back to get input y coordinates replace = [ (metric.y, y*Ly/(2*pi) ) ] @@ -48,12 +50,12 @@ def Div_par_K_Grad_par(f: Expr, K: Expr) -> Expr: for name, expr in ( ("input_field", f), - ("K", K), - ("FV_Div_par_solution", Div_par(f * K)), - ("FV_Div_par_K_Grad_par_solution", Div_par_K_Grad_par(f, K)), - ("FV_Div_par_K_Grad_par_mod_solution", Div_par_K_Grad_par(f, K)), - ("FV_Div_par_mod_solution", Div_par(f * K)), - ("FV_Div_par_fvv_solution", Div_par(f * K * K)), + ("v", v), + ("FV_Div_par_solution", Div_par(f * v)), + ("FV_Div_par_K_Grad_par_solution", Div_par_K_Grad_par(f, v)), + ("FV_Div_par_K_Grad_par_mod_solution", Div_par_K_Grad_par(f, v)), + ("FV_Div_par_mod_solution", Div_par(f * v)), + ("FV_Div_par_fvv_solution", Div_par(f * v * v)), ): expr_str = exprToStr(expr.subs(replace)) print(f"{name} = {expr_str}") diff --git a/tests/MMS/spatial/finite-volume/runtest b/tests/MMS/spatial/finite-volume/runtest index a836f5e735..5a02be1e50 100755 --- a/tests/MMS/spatial/finite-volume/runtest +++ b/tests/MMS/spatial/finite-volume/runtest @@ -17,13 +17,27 @@ from numpy import array, log, polyfit DIRECTORY = "data" NPROC = 2 MTHREAD = 2 -OPERATORS = ( - "FV_Div_par", - "FV_Div_par_K_Grad_par", - "FV_Div_par_K_Grad_par_mod", - "FV_Div_par_mod", - "FV_Div_par_fvv", -) +OPERATORS = { + # Slope-limiters necessarily reduce the accuracy in places + "FV_Div_par_MC": 1.5, + "FV_Div_par_mod_MC": 1.5, + "FV_Div_par_fvv_MC": 1.5, + + "FV_Div_par_Upwind": 1, + "FV_Div_par_mod_Upwind": 1, + "FV_Div_par_fvv_Upwind": 1, + + "FV_Div_par_Fromm": 1.5, + "FV_Div_par_mod_Fromm": 1.5, + "FV_Div_par_fvv_Fromm": 1.5, + + "FV_Div_par_MinMod": 1.5, + "FV_Div_par_mod_MinMod": 1.5, + "FV_Div_par_fvv_MinMod": 1.5, + + "FV_Div_par_K_Grad_par": 2, + "FV_Div_par_K_Grad_par_mod": 2, +} # Resolution in y and z NLIST = [8, 16, 32, 64] dx = 1.0 / array(NLIST) @@ -56,10 +70,9 @@ def assert_convergence(error, dx, name, expected) -> bool: return success -def run_fv_operators(nz: int, name: str) -> dict[str, float]: +def run_fv_operators(nz: int) -> dict[str, float]: # Command to run - args = f"MZ={nz} mesh:ny={nz} {name}" - cmd = f"./fv_mms {args}" + cmd = f"./fv_mms MZ={nz} mesh:ny={nz}" print(f"Running command: {cmd}", end="") # Launch using MPI @@ -97,16 +110,13 @@ def transpose( return result -def check_fv_operators(name: str, case: dict) -> bool: +def test_fv_operators() -> bool: failures = [] - order = case["order"] - args = case["args"] - all_errors = [] for n in NLIST: - errors = run_fv_operators(n, args) + errors = run_fv_operators(n) all_errors.append(errors) for operator in OPERATORS: @@ -116,13 +126,12 @@ def check_fv_operators(name: str, case: dict) -> bool: print(f"{operator} errors: l-2 {l_2:f} l-inf {l_inf:f}") final_errors = transpose(all_errors) - for operator in OPERATORS: - test_name = f"{operator} {name}" + for operator, order in OPERATORS.items(): success = assert_convergence( - final_errors[operator]["l_2"], dx, test_name, order + final_errors[operator]["l_2"], dx, operator, order ) if not success: - failures.append(test_name) + failures.append(operator) return final_errors, failures @@ -138,9 +147,8 @@ def make_plots(cases: dict[str, dict]): fig, axes = plt.subplots(1, num_operators, figsize=(num_operators * 4, 4)) for ax, operator in zip(axes, OPERATORS): - for name, case in cases.items(): - ax.loglog(dx, case[operator]["l_2"], "-", label=f"{name} $l_2$") - ax.loglog(dx, case[operator]["l_inf"], "--", label=f"{name} $l_\\inf$") + ax.loglog(dx, cases[operator]["l_2"], "-", label="$l_2$") + ax.loglog(dx, cases[operator]["l_inf"], "--", label="$l_\\inf$") ax.legend(loc="upper left") ax.grid() ax.set_title(f"Error scaling for {operator}") @@ -177,27 +185,14 @@ if __name__ == "__main__": args = parser.parse_args() - success = True - failures = [] - - cases = { - "default": { - "order": 2, - "args": "", - }, - } - - for name, case in cases.items(): - error2, failures_ = check_fv_operators(name, case) - case.update(error2) - failures.extend(failures_) - success &= len(failures) == 0 + error2, failures = test_fv_operators() + success = len(failures) == 0 if args.dump_errors: - pathlib.Path(args.dump_errors).write_text(json.dumps(cases)) + pathlib.Path(args.dump_errors).write_text(json.dumps(error2)) if args.make_plots or args.show_plots: - make_plots(cases) + make_plots(error2) if success: print("\nAll tests passed") From a230642c788a491ce0b8264b3c2ce24d84e326b0 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Fri, 7 Nov 2025 14:25:53 +0000 Subject: [PATCH 06/38] Port `Superbee` finite volume limiter from Hermes-3 Includes bug fix: ```diff - BoutReal gL = n.c - n.L; - BoutReal gR = n.R - n.c; + BoutReal gL = n.c - n.m; + BoutReal gR = n.p - n.c; ``` --- include/bout/fv_ops.hxx | 45 ++++++++++++++++++++++ tests/MMS/spatial/finite-volume/fv_mms.cxx | 7 ++++ tests/MMS/spatial/finite-volume/runtest | 4 ++ 3 files changed, 56 insertions(+) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index bd0fa812d8..37d7b3b6fe 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -170,6 +170,51 @@ private: } }; +/// Superbee limiter +/// +/// This corresponds to the limiter function +/// φ(r) = max(0, min(2r, 1), min(r,2) +/// +/// The value at cell right (i.e. i + 1/2) is: +/// +/// n.R = n.c - φ(r) (n.c - (n.p + n.c)/2) +/// = n.c + φ(r) (n.p - n.c)/2 +/// +/// Four regimes: +/// a) r < 1/2 -> φ(r) = 2r +/// n.R = n.c + gL +/// b) 1/2 < r < 1 -> φ(r) = 1 +/// n.R = n.c + gR/2 +/// c) 1 < r < 2 -> φ(r) = r +/// n.R = n.c + gL/2 +/// d) 2 < r -> φ(r) = 2 +/// n.R = n.c + gR +/// +/// where the left and right gradients are: +/// gL = n.c - n.m +/// gR = n.p - n.c +/// +struct Superbee { + void operator()(Stencil1D& n) { + BoutReal gL = n.c - n.m; + BoutReal gR = n.p - n.c; + + // r = gL / gR + // Limiter is φ(r) + if (gL * gR < 0) { + // Different signs => Zero gradient + n.L = n.R = n.c; + } else { + BoutReal sign = SIGN(gL); + gL = fabs(gL); + gR = fabs(gR); + BoutReal half_slope = sign * BOUTMAX(BOUTMIN(gL, 0.5 * gR), BOUTMIN(gR, 0.5 * gL)); + n.L = n.c - half_slope; + n.R = n.c + half_slope; + } + } +}; + /*! * Communicate fluxes between processors * Takes values in guard cells, and adds them to cells diff --git a/tests/MMS/spatial/finite-volume/fv_mms.cxx b/tests/MMS/spatial/finite-volume/fv_mms.cxx index 09f9986e72..edf4bbc16a 100644 --- a/tests/MMS/spatial/finite-volume/fv_mms.cxx +++ b/tests/MMS/spatial/finite-volume/fv_mms.cxx @@ -74,6 +74,13 @@ int main(int argc, char** argv) { fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), "_MinMod"); + fv_op_test("FV_Div_par", dump, input, FV::Div_par(input, v, v), + "_Superbee"); + fv_op_test("FV_Div_par_mod", dump, input, + FV::Div_par_mod(input, v, v, flow_ylow), "_Superbee"); + fv_op_test("FV_Div_par_fvv", dump, input, FV::Div_par_fvv(input, v, v), + "_Superbee"); + fv_op_test("FV_Div_par_K_Grad_par", dump, input, FV::Div_par_K_Grad_par(v, input)); fv_op_test("FV_Div_par_K_Grad_par_mod", dump, input, Div_par_K_Grad_par_mod(v, input, flow_ylow)); diff --git a/tests/MMS/spatial/finite-volume/runtest b/tests/MMS/spatial/finite-volume/runtest index 5a02be1e50..b38a6359ac 100755 --- a/tests/MMS/spatial/finite-volume/runtest +++ b/tests/MMS/spatial/finite-volume/runtest @@ -35,6 +35,10 @@ OPERATORS = { "FV_Div_par_mod_MinMod": 1.5, "FV_Div_par_fvv_MinMod": 1.5, + "FV_Div_par_Superbee": 1.5, + "FV_Div_par_mod_Superbee": 1.5, + "FV_Div_par_fvv_Superbee": 1.5, + "FV_Div_par_K_Grad_par": 2, "FV_Div_par_K_Grad_par_mod": 2, } From 30dc67002eff2e2f216a937d8f10e3496e8281c8 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Fri, 7 Nov 2025 14:55:39 +0000 Subject: [PATCH 07/38] Fix clang-tidy warnings for fv_ops --- include/bout/fv_ops.hxx | 172 +++++++++---------- src/mesh/difops.cxx | 85 +++++----- src/mesh/fv_ops.cxx | 185 +++++++++++---------- tests/MMS/spatial/finite-volume/fv_mms.cxx | 1 + 4 files changed, 223 insertions(+), 220 deletions(-) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 37d7b3b6fe..cd9a3536c1 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -5,33 +5,38 @@ #ifndef BOUT_FV_OPS_H #define BOUT_FV_OPS_H -#include - +#include "bout/assert.hxx" #include "bout/bout_types.hxx" +#include "bout/boutexception.hxx" +#include "bout/build_defines.hxx" +#include "bout/coordinates.hxx" +#include "bout/field.hxx" +#include "bout/field2d.hxx" #include "bout/field3d.hxx" #include "bout/globals.hxx" #include "bout/mesh.hxx" -#include "bout/output_bout_types.hxx" // NOLINT(unused-includes) +#include "bout/output_bout_types.hxx" // NOLINT(unused-includes, misc-include-cleaner) #include "bout/region.hxx" #include "bout/utils.hxx" #include "bout/vector2d.hxx" +#include + namespace FV { /*! * Div ( a Grad_perp(f) ) -- ∇⊥ ( a ⋅ ∇⊥ f) -- Vorticity */ -Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& x); +Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f); [[deprecated("Please use Div_a_Grad_perp instead")]] inline Field3D -Div_a_Laplace_perp(const Field3D& a, const Field3D& x) { - return Div_a_Grad_perp(a, x); +Div_a_Laplace_perp(const Field3D& a, const Field3D& f) { + return Div_a_Grad_perp(a, f); } /*! * Divergence of a parallel diffusion Div( k * Grad_par(f) ) */ -const Field3D Div_par_K_Grad_par(const Field3D& k, const Field3D& f, - bool bndry_flux = true); +Field3D Div_par_K_Grad_par(const Field3D& k, const Field3D& f, bool bndry_flux = true); /*! * 4th-order derivative in Y, using derivatives @@ -53,7 +58,7 @@ const Field3D Div_par_K_Grad_par(const Field3D& k, const Field3D& f, * * No fluxes through domain boundaries */ -const Field3D D4DY4(const Field3D& d, const Field3D& f); +Field3D D4DY4(const Field3D& d, const Field3D& f); /*! * 4th-order dissipation term @@ -71,18 +76,24 @@ const Field3D D4DY4(const Field3D& d, const Field3D& f); * f_2 | f_1 | f_0 | * f_b */ -const Field3D D4DY4_Index(const Field3D& f, bool bndry_flux = true); +Field3D D4DY4_Index(const Field3D& f, bool bndry_flux = true); /*! * Stencil used for Finite Volume calculations * which includes cell face values L and R */ struct Stencil1D { - // Cell centre values - BoutReal c, m, p, mm, pp; - - // Left and right cell face values - BoutReal L, R; + /// Cell centre values + BoutReal c; + BoutReal m; + BoutReal p; + BoutReal mm = BoutNaN; + BoutReal pp = BoutNaN; + + /// Left cell face value + BoutReal L = BoutNaN; + /// Right cell face value + BoutReal R = BoutNaN; }; /*! @@ -97,8 +108,8 @@ struct Upwind { */ struct Fromm { void operator()(Stencil1D& n) { - n.L = n.c - 0.25 * (n.p - n.m); - n.R = n.c + 0.25 * (n.p - n.m); + n.L = n.c - (0.25 * (n.p - n.m)); + n.R = n.c + (0.25 * (n.p - n.m)); } }; @@ -114,9 +125,9 @@ struct MinMod { void operator()(Stencil1D& n) { // Choose the gradient within the cell // as the minimum (smoothest) solution - BoutReal slope = _minmod(n.p - n.c, n.c - n.m); - n.L = n.c - 0.5 * slope; - n.R = n.c + 0.5 * slope; + const BoutReal slope = _minmod(n.p - n.c, n.c - n.m); + n.L = n.c - (0.5 * slope); + n.R = n.c + (0.5 * slope); } private: @@ -127,7 +138,7 @@ private: * returns zero, otherwise chooses the value * with the minimum magnitude. */ - BoutReal _minmod(BoutReal a, BoutReal b) { + static BoutReal _minmod(BoutReal a, BoutReal b) { if (a * b <= 0.0) { return 0.0; } @@ -149,17 +160,17 @@ private: */ struct MC { void operator()(Stencil1D& n) { - BoutReal slope = minmod(2. * (n.p - n.c), // 2*right difference - 0.5 * (n.p - n.m), // Central difference - 2. * (n.c - n.m)); // 2*left difference - n.L = n.c - 0.5 * slope; - n.R = n.c + 0.5 * slope; + const BoutReal slope = minmod(2. * (n.p - n.c), // 2*right difference + 0.5 * (n.p - n.m), // Central difference + 2. * (n.c - n.m)); // 2*left difference + n.L = n.c - (0.5 * slope); + n.R = n.c + (0.5 * slope); } private: // Return zero if any signs are different // otherwise return the value with the minimum magnitude - BoutReal minmod(BoutReal a, BoutReal b, BoutReal c) { + static BoutReal minmod(BoutReal a, BoutReal b, BoutReal c) { // if any of the signs are different, return zero gradient if ((a * b <= 0.0) || (a * c <= 0.0)) { return 0.0; @@ -196,8 +207,8 @@ private: /// struct Superbee { void operator()(Stencil1D& n) { - BoutReal gL = n.c - n.m; - BoutReal gR = n.p - n.c; + const BoutReal gL = n.c - n.m; + const BoutReal gR = n.p - n.c; // r = gL / gR // Limiter is φ(r) @@ -205,10 +216,11 @@ struct Superbee { // Different signs => Zero gradient n.L = n.R = n.c; } else { - BoutReal sign = SIGN(gL); - gL = fabs(gL); - gR = fabs(gR); - BoutReal half_slope = sign * BOUTMAX(BOUTMIN(gL, 0.5 * gR), BOUTMIN(gR, 0.5 * gL)); + const BoutReal sign = SIGN(gL); + const BoutReal abs_gL = fabs(gL); + const BoutReal abs_gR = fabs(gR); + const BoutReal half_slope = + sign * BOUTMAX(BOUTMIN(abs_gL, 0.5 * abs_gR), BOUTMIN(abs_gR, 0.5 * abs_gL)); n.L = n.c - half_slope; n.R = n.c + half_slope; } @@ -238,13 +250,13 @@ void communicateFluxes(Field3D& f); /// /// NB: Uses to/from FieldAligned coordinates template -const Field3D Div_par(const Field3D& f_in, const Field3D& v_in, - const Field3D& wave_speed_in, bool fixflux = true) { +Field3D Div_par(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, + bool fixflux = true) { ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); ASSERT1_FIELDS_COMPATIBLE(f_in, wave_speed_in); - Mesh* mesh = f_in.getMesh(); + Mesh const* mesh = f_in.getMesh(); CellEdges cellboundary; @@ -264,29 +276,17 @@ const Field3D Div_par(const Field3D& f_in, const Field3D& v_in, Field3D result{zeroFrom(f)}; - // Only need one guard cell, so no need to communicate fluxes - // Instead calculate in guard cells to preserve fluxes - int ys = mesh->ystart - 1; - int ye = mesh->yend + 1; - for (int i = mesh->xstart; i <= mesh->xend; i++) { + const bool is_periodic_y = mesh->periodicY(i); + const bool is_first_y = mesh->firstY(i); + const bool is_last_y = mesh->lastY(i); - if (!mesh->firstY(i) || mesh->periodicY(i)) { - // Calculate in guard cell to get fluxes consistent between processors - ys = mesh->ystart - 1; - } else { - // Don't include the boundary cell. Note that this implies special - // handling of boundaries later - ys = mesh->ystart; - } - - if (!mesh->lastY(i) || mesh->periodicY(i)) { - // Calculate in guard cells - ye = mesh->yend + 1; - } else { - // Not in boundary cells - ye = mesh->yend; - } + // Only need one guard cell, so no need to communicate fluxes Instead + // calculate in guard cells to get fluxes consistent between processors, but + // don't include the boundary cell. Note that this implies special handling + // of boundaries later + const int ys = (!is_first_y || is_periodic_y) ? mesh->ystart - 1 : mesh->ystart; + const int ye = (!is_last_y || is_periodic_y) ? mesh->yend + 1 : mesh->yend; for (int j = ys; j <= ye; j++) { // Pre-calculate factors which multiply fluxes @@ -295,16 +295,16 @@ const Field3D Div_par(const Field3D& f_in, const Field3D& v_in, BoutReal common_factor = (coord->J(i, j) + coord->J(i, j + 1)) / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1))); - BoutReal flux_factor_rc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - BoutReal flux_factor_rp = + const BoutReal flux_factor_rc = common_factor / (coord->dy(i, j) * coord->J(i, j)); + const BoutReal flux_factor_rp = common_factor / (coord->dy(i, j + 1) * coord->J(i, j + 1)); // For left cell boundaries common_factor = (coord->J(i, j) + coord->J(i, j - 1)) / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1))); - BoutReal flux_factor_lc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - BoutReal flux_factor_lm = + const BoutReal flux_factor_lc = common_factor / (coord->dy(i, j) * coord->J(i, j)); + const BoutReal flux_factor_lm = common_factor / (coord->dy(i, j - 1) * coord->J(i, j - 1)); #endif for (int k = mesh->zstart; k <= mesh->zend; k++) { @@ -347,23 +347,23 @@ const Field3D Div_par(const Field3D& f_in, const Field3D& v_in, // Calculate velocity at right boundary (y+1/2) BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j + 1, k)); - BoutReal flux; + BoutReal flux = NAN; - if (mesh->lastY(i) && (j == mesh->yend) && !mesh->periodicY(i)) { + if (is_last_y && (j == mesh->yend) && !is_periodic_y) { // Last point in domain - BoutReal bndryval = 0.5 * (s.c + s.p); + const BoutReal bndryval = 0.5 * (s.c + s.p); if (fixflux) { // Use mid-point to be consistent with boundary conditions flux = bndryval * vpar; } else { // Add flux due to difference in boundary values - flux = s.R * vpar + wave_speed(i, j, k) * (s.R - bndryval); + flux = (s.R * vpar) + (wave_speed(i, j, k) * (s.R - bndryval)); } } else { // Maximum wave speed in the two cells - BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k)); + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k)); if (vpar > amax) { // Supersonic flow out of this cell @@ -385,20 +385,20 @@ const Field3D Div_par(const Field3D& f_in, const Field3D& v_in, vpar = 0.5 * (v(i, j, k) + v(i, j - 1, k)); - if (mesh->firstY(i) && (j == mesh->ystart) && !mesh->periodicY(i)) { + if (is_first_y && (j == mesh->ystart) && !is_periodic_y) { // First point in domain - BoutReal bndryval = 0.5 * (s.c + s.m); + const BoutReal bndryval = 0.5 * (s.c + s.m); if (fixflux) { // Use mid-point to be consistent with boundary conditions flux = bndryval * vpar; } else { // Add flux due to difference in boundary values - flux = s.L * vpar - wave_speed(i, j, k) * (s.L - bndryval); + flux = (s.L * vpar) - (wave_speed(i, j, k) * (s.L - bndryval)); } } else { // Maximum wave speed in the two cells - BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k)); + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k)); if (vpar < -amax) { // Supersonic out of this cell @@ -432,11 +432,11 @@ const Field3D Div_par(const Field3D& f_in, const Field3D& v_in, * */ template -const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { +Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { ASSERT1(n_in.getLocation() == v.getLocation()); ASSERT1_FIELDS_COMPATIBLE(n_in, v.x); - Mesh* mesh = n_in.getMesh(); + const Mesh* mesh = n_in.getMesh(); CellEdges cellboundary; @@ -455,10 +455,10 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { // Calculate velocities - BoutReal vU = 0.25 * (vz[i.zp()] + vz[i]) * (coord->J[i.zp()] + coord->J[i]); - BoutReal vD = 0.25 * (vz[i.zm()] + vz[i]) * (coord->J[i.zm()] + coord->J[i]); - BoutReal vL = 0.25 * (vx[i.xm()] + vx[i]) * (coord->J[i.xm()] + coord->J[i]); - BoutReal vR = 0.25 * (vx[i.xp()] + vx[i]) * (coord->J[i.xp()] + coord->J[i]); + const BoutReal vU = 0.25 * (vz[i.zp()] + vz[i]) * (coord->J[i.zp()] + coord->J[i]); + const BoutReal vD = 0.25 * (vz[i.zm()] + vz[i]) * (coord->J[i.zm()] + coord->J[i]); + const BoutReal vL = 0.25 * (vx[i.xm()] + vx[i]) * (coord->J[i.xm()] + coord->J[i]); + const BoutReal vR = 0.25 * (vx[i.xp()] + vx[i]) * (coord->J[i.xp()] + coord->J[i]); // X direction Stencil1D s; @@ -473,7 +473,7 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { if ((i.x() == mesh->xend) && (mesh->lastX())) { // At right boundary in X if (bndry_flux) { - BoutReal flux; + BoutReal flux = NAN; if (vR > 0.0) { // Flux to boundary flux = vR * s.R; @@ -488,7 +488,7 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { // Not at a boundary if (vR > 0.0) { // Flux out into next cell - BoutReal flux = vR * s.R; + const BoutReal flux = vR * s.R; result[i] += flux / (coord->dx[i] * coord->J[i]); result[i.xp()] -= flux / (coord->dx[i.xp()] * coord->J[i.xp()]); } @@ -500,7 +500,7 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { // At left boundary in X if (bndry_flux) { - BoutReal flux; + BoutReal flux = NAN; if (vL < 0.0) { // Flux to boundary flux = vL * s.L; @@ -514,7 +514,7 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { } else { // Not at a boundary if (vL < 0.0) { - BoutReal flux = vL * s.L; + const BoutReal flux = vL * s.L; result[i] -= flux / (coord->dx[i] * coord->J[i]); result[i.xm()] += flux / (coord->dx[i.xm()] * coord->J[i.xm()]); } @@ -531,12 +531,12 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { cellboundary(s); if (vU > 0.0) { - BoutReal flux = vU * s.R; + const BoutReal flux = vU * s.R; result[i] += flux / (coord->J[i] * coord->dz[i]); result[i.zp()] -= flux / (coord->J[i.zp()] * coord->dz[i.zp()]); } if (vD < 0.0) { - BoutReal flux = vD * s.L; + const BoutReal flux = vD * s.L; result[i] -= flux / (coord->J[i] * coord->dz[i]); result[i.zm()] += flux / (coord->J[i.zm()] * coord->dz[i.zm()]); } @@ -556,13 +556,13 @@ const Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { // Y velocities on y boundaries - BoutReal vU = 0.25 * (vy[i] + vy[i.yp()]) * (coord->J[i] + coord->J[i.yp()]); - BoutReal vD = 0.25 * (vy[i] + vy[i.ym()]) * (coord->J[i] + coord->J[i.ym()]); + const BoutReal vU = 0.25 * (vy[i] + vy[i.yp()]) * (coord->J[i] + coord->J[i.yp()]); + const BoutReal vD = 0.25 * (vy[i] + vy[i.ym()]) * (coord->J[i] + coord->J[i.ym()]); // n (advected quantity) on y boundaries // Note: Use unshifted n_in variable - BoutReal nU = 0.5 * (n[i] + n[i.yp()]); - BoutReal nD = 0.5 * (n[i] + n[i.ym()]); + const BoutReal nU = 0.5 * (n[i] + n[i.yp()]); + const BoutReal nD = 0.5 * (n[i] + n[i.ym()]); yresult[i] = (nU * vU - nD * vD) / (coord->J[i] * coord->dy[i]); } diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 09d2441951..56773f3c4c 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -25,19 +25,20 @@ #include "bout/build_defines.hxx" -#include -#include -#include -#include -#include -#include -#include -#include - -#include // Delp2 uses same coefficients as inversion code - -#include -#include +#include "bout/assert.hxx" +#include "bout/derivs.hxx" +#include "bout/difops.hxx" +#include "bout/fft.hxx" +#include "bout/field2d.hxx" +#include "bout/globals.hxx" +#include "bout/interpolation.hxx" +#include "bout/invert_laplace.hxx" // Delp2 uses same coefficients as inversion code +#include "bout/msg_stack.hxx" +#include "bout/region.hxx" +#include "bout/solver.hxx" +#include "bout/unused.hxx" +#include "bout/utils.hxx" +#include "bout/vecops.hxx" #include @@ -366,14 +367,14 @@ Field3D Div_par_K_Grad_par(const Field3D& kY, const Field3D& f, CELL_LOC outloc) + Div_par(kY, outloc) * Grad_par(f, outloc); } -Field3D Div_par_K_Grad_par_mod(const Field3D& Kin, const Field3D& fin, - Field3D& flow_ylow, bool bndry_flux) { +Field3D Div_par_K_Grad_par_mod(const Field3D& Kin, const Field3D& fin, Field3D& flow_ylow, + bool bndry_flux) { TRACE("FV::Div_par_K_Grad_par_mod"); ASSERT2(Kin.getLocation() == fin.getLocation()); - Mesh* mesh = Kin.getMesh(); - Coordinates* coord = fin.getCoordinates(); + const Mesh* mesh = Kin.getMesh(); + const Coordinates* coord = fin.getCoordinates(); if (Kin.hasParallelSlices() && fin.hasParallelSlices()) { // Using parallel slices. @@ -395,18 +396,22 @@ Field3D Div_par_K_Grad_par_mod(const Field3D& Kin, const Field3D& fin, const auto iym = i.ym(); // Upper cell edge - const BoutReal c_up = 0.5 * (Kin[i] + K_up[iyp]); // K at the upper boundary - const BoutReal J_up = 0.5 * (coord->J[i] + coord->J.yup()[iyp]); // Jacobian at boundary + const BoutReal c_up = 0.5 * (Kin[i] + K_up[iyp]); // K at the upper boundary + const BoutReal J_up = + 0.5 * (coord->J[i] + coord->J.yup()[iyp]); // Jacobian at boundary const BoutReal g_22_up = 0.5 * (coord->g_22[i] + coord->g_22.yup()[iyp]); - const BoutReal gradient_up = 2. * (f_up[iyp] - fin[i]) / (coord->dy[i] + coord->dy.yup()[iyp]); + const BoutReal gradient_up = + 2. * (f_up[iyp] - fin[i]) / (coord->dy[i] + coord->dy.yup()[iyp]); const BoutReal flux_up = c_up * J_up * gradient_up / g_22_up; // Lower cell edge - const BoutReal c_down = 0.5 * (Kin[i] + K_down[iym]); // K at the lower boundary - const BoutReal J_down = 0.5 * (coord->J[i] + coord->J.ydown()[iym]); // Jacobian at boundary + const BoutReal c_down = 0.5 * (Kin[i] + K_down[iym]); // K at the lower boundary + const BoutReal J_down = + 0.5 * (coord->J[i] + coord->J.ydown()[iym]); // Jacobian at boundary const BoutReal g_22_down = 0.5 * (coord->g_22[i] + coord->g_22.ydown()[iym]); - const BoutReal gradient_down = 2. * (fin[i] - f_down[iym]) / (coord->dy[i] + coord->dy.ydown()[iym]); + const BoutReal gradient_down = + 2. * (fin[i] - f_down[iym]) / (coord->dy[i] + coord->dy.ydown()[iym]); const BoutReal flux_down = c_down * J_down * gradient_down / g_22_down; @@ -425,35 +430,32 @@ Field3D Div_par_K_Grad_par_mod(const Field3D& Kin, const Field3D& fin, BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { // Calculate flux at upper surface - + const auto ix = i.x(); + const auto iy = i.y(); const auto iyp = i.yp(); const auto iym = i.ym(); - if (bndry_flux || mesh->periodicY(i.x()) || !mesh->lastY(i.x()) - || (i.y() != mesh->yend)) { - - BoutReal c = 0.5 * (K[i] + K[iyp]); // K at the upper boundary - BoutReal J = 0.5 * (coord->J[i] + coord->J[iyp]); // Jacobian at boundary - BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iyp]); + const bool is_periodic_y = mesh->periodicY(ix); - BoutReal gradient = 2. * (f[iyp] - f[i]) / (coord->dy[i] + coord->dy[iyp]); + if (bndry_flux || is_periodic_y || !mesh->lastY(ix) || (iy != mesh->yend)) { + const BoutReal c = 0.5 * (K[i] + K[iyp]); // K at the upper boundary + const BoutReal J = 0.5 * (coord->J[i] + coord->J[iyp]); // Jacobian at boundary + const BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iyp]); + const BoutReal gradient = 2. * (f[iyp] - f[i]) / (coord->dy[i] + coord->dy[iyp]); - BoutReal flux = c * J * gradient / g_22; + const BoutReal flux = c * J * gradient / g_22; result[i] += flux / (coord->dy[i] * coord->J[i]); } // Calculate flux at lower surface - if (bndry_flux || mesh->periodicY(i.x()) || !mesh->firstY(i.x()) - || (i.y() != mesh->ystart)) { - BoutReal c = 0.5 * (K[i] + K[iym]); // K at the lower boundary - BoutReal J = 0.5 * (coord->J[i] + coord->J[iym]); // Jacobian at boundary + if (bndry_flux || is_periodic_y || !mesh->firstY(ix) || (iy != mesh->ystart)) { + const BoutReal c = 0.5 * (K[i] + K[iym]); // K at the lower boundary + const BoutReal J = 0.5 * (coord->J[i] + coord->J[iym]); // Jacobian at boundary + const BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iym]); + const BoutReal gradient = 2. * (f[i] - f[iym]) / (coord->dy[i] + coord->dy[iym]); - BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iym]); - - BoutReal gradient = 2. * (f[i] - f[iym]) / (coord->dy[i] + coord->dy[iym]); - - BoutReal flux = c * J * gradient / g_22; + const BoutReal flux = c * J * gradient / g_22; result[i] -= flux / (coord->dy[i] * coord->J[i]); flow_ylow[i] = -flux * coord->dx[i] * coord->dz[i]; @@ -467,7 +469,6 @@ Field3D Div_par_K_Grad_par_mod(const Field3D& Kin, const Field3D& fin, return result; } - /******************************************************************************* * Delp2 * perpendicular Laplacian operator diff --git a/src/mesh/fv_ops.cxx b/src/mesh/fv_ops.cxx index fab8beb794..6b8d8a6f21 100644 --- a/src/mesh/fv_ops.cxx +++ b/src/mesh/fv_ops.cxx @@ -1,7 +1,16 @@ -#include -#include -#include -#include +#include "bout/fv_ops.hxx" + +#include "bout/assert.hxx" +#include "bout/bout_types.hxx" +#include "bout/boutexception.hxx" +#include "bout/build_config.hxx" +#include "bout/coordinates.hxx" +#include "bout/field2d.hxx" +#include "bout/field3d.hxx" +#include "bout/globals.hxx" +#include "bout/msg_stack.hxx" +#include "bout/region.hxx" +#include "bout/utils.hxx" namespace { template @@ -33,28 +42,19 @@ Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f) { // Flux in x - int xs = mesh->xstart - 1; - int xe = mesh->xend; - - /* - if(mesh->firstX()) - xs += 1; - */ - /* - if(mesh->lastX()) - xe -= 1; - */ + const int xs = mesh->xstart - 1; + const int xe = mesh->xend; for (int i = xs; i <= xe; i++) { for (int j = mesh->ystart; j <= mesh->yend; j++) { for (int k = mesh->zstart; k <= mesh->zend; k++) { // Calculate flux from i to i+1 - BoutReal fout = 0.5 * (a(i, j, k) + a(i + 1, j, k)) - * (coord->J(i, j, k) * coord->g11(i, j, k) - + coord->J(i + 1, j, k) * coord->g11(i + 1, j, k)) - * (f(i + 1, j, k) - f(i, j, k)) - / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)); + const BoutReal fout = 0.5 * (a(i, j, k) + a(i + 1, j, k)) + * (coord->J(i, j, k) * coord->g11(i, j, k) + + coord->J(i + 1, j, k) * coord->g11(i + 1, j, k)) + * (f(i + 1, j, k) - f(i, j, k)) + / (coord->dx(i, j, k) + coord->dx(i + 1, j, k)); result(i, j, k) += fout / (coord->dx(i, j, k) * coord->J(i, j, k)); result(i + 1, j, k) -= fout / (coord->dx(i + 1, j, k) * coord->J(i + 1, j, k)); @@ -178,14 +178,13 @@ Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f) { return result; } -const Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, - bool bndry_flux) { +Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, bool bndry_flux) { ASSERT2(Kin.getLocation() == fin.getLocation()); - Mesh* mesh = Kin.getMesh(); + const Mesh* mesh = Kin.getMesh(); - bool use_parallel_slices = (Kin.hasParallelSlices() && fin.hasParallelSlices()); + const bool use_parallel_slices = (Kin.hasParallelSlices() && fin.hasParallelSlices()); const auto& K = use_parallel_slices ? Kin : toFieldAligned(Kin, "RGN_NOX"); const auto& f = use_parallel_slices ? fin : toFieldAligned(fin, "RGN_NOX"); @@ -209,13 +208,13 @@ const Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, if (bndry_flux || mesh->periodicY(i.x()) || !mesh->lastY(i.x()) || (i.y() != mesh->yend)) { - BoutReal c = 0.5 * (K[i] + Kup[iyp]); // K at the upper boundary - BoutReal J = 0.5 * (coord->J[i] + coord->J[iyp]); // Jacobian at boundary - BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iyp]); + const BoutReal c = 0.5 * (K[i] + Kup[iyp]); // K at the upper boundary + const BoutReal J = 0.5 * (coord->J[i] + coord->J[iyp]); // Jacobian at boundary + const BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iyp]); - BoutReal gradient = 2. * (fup[iyp] - f[i]) / (coord->dy[i] + coord->dy[iyp]); + const BoutReal gradient = 2. * (fup[iyp] - f[i]) / (coord->dy[i] + coord->dy[iyp]); - BoutReal flux = c * J * gradient / g_22; + const BoutReal flux = c * J * gradient / g_22; result[i] += flux / (coord->dy[i] * coord->J[i]); } @@ -223,14 +222,15 @@ const Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, // Calculate flux at lower surface if (bndry_flux || mesh->periodicY(i.x()) || !mesh->firstY(i.x()) || (i.y() != mesh->ystart)) { - BoutReal c = 0.5 * (K[i] + Kdown[iym]); // K at the lower boundary - BoutReal J = 0.5 * (coord->J[i] + coord->J[iym]); // Jacobian at boundary + const BoutReal c = 0.5 * (K[i] + Kdown[iym]); // K at the lower boundary + const BoutReal J = 0.5 * (coord->J[i] + coord->J[iym]); // Jacobian at boundary - BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iym]); + const BoutReal g_22 = 0.5 * (coord->g_22[i] + coord->g_22[iym]); - BoutReal gradient = 2. * (f[i] - fdown[iym]) / (coord->dy[i] + coord->dy[iym]); + const BoutReal gradient = + 2. * (f[i] - fdown[iym]) / (coord->dy[i] + coord->dy[iym]); - BoutReal flux = c * J * gradient / g_22; + const BoutReal flux = c * J * gradient / g_22; result[i] -= flux / (coord->dy[i] * coord->J[i]); } @@ -244,10 +244,10 @@ const Field3D Div_par_K_Grad_par(const Field3D& Kin, const Field3D& fin, return result; } -const Field3D D4DY4(const Field3D& d_in, const Field3D& f_in) { +Field3D D4DY4(const Field3D& d_in, const Field3D& f_in) { ASSERT1_FIELDS_COMPATIBLE(d_in, f_in); - Mesh* mesh = d_in.getMesh(); + const Mesh* mesh = d_in.getMesh(); Coordinates* coord = f_in.getCoordinates(); @@ -263,9 +263,9 @@ const Field3D D4DY4(const Field3D& d_in, const Field3D& f_in) { for (int i = mesh->xstart; i <= mesh->xend; i++) { // Check for boundaries - bool yperiodic = mesh->periodicY(i); - bool has_upper_boundary = !yperiodic && mesh->lastY(i); - bool has_lower_boundary = !yperiodic && mesh->firstY(i); + const bool yperiodic = mesh->periodicY(i); + const bool has_upper_boundary = !yperiodic && mesh->lastY(i); + const bool has_lower_boundary = !yperiodic && mesh->firstY(i); // Always calculate fluxes at upper Y cell boundary const int ystart = @@ -281,15 +281,15 @@ const Field3D D4DY4(const Field3D& d_in, const Field3D& f_in) { for (int j = ystart; j <= yend; j++) { for (int k = mesh->zstart; k <= mesh->zend; k++) { - BoutReal dy3 = SQ(coord->dy(i, j, k)) * coord->dy(i, j, k); + const BoutReal dy3 = SQ(coord->dy(i, j, k)) * coord->dy(i, j, k); // 3rd derivative at upper boundary - BoutReal d3fdy3 = + const BoutReal d3fdy3 = (f(i, j + 2, k) - 3. * f(i, j + 1, k) + 3. * f(i, j, k) - f(i, j - 1, k)) / dy3; - BoutReal flux = 0.5 * (d(i, j, k) + d(i, j + 1, k)) - * (coord->J(i, j, k) + coord->J(i, j + 1, k)) * d3fdy3; + const BoutReal flux = 0.5 * (d(i, j, k) + d(i, j + 1, k)) + * (coord->J(i, j, k) + coord->J(i, j + 1, k)) * d3fdy3; result(i, j, k) += flux / (coord->J(i, j, k) * coord->dy(i, j, k)); result(i, j + 1, k) -= flux / (coord->J(i, j + 1, k) * coord->dy(i, j + 1, k)); @@ -301,8 +301,8 @@ const Field3D D4DY4(const Field3D& d_in, const Field3D& f_in) { return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result; } -const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { - Mesh* mesh = f_in.getMesh(); +Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { + const Mesh* mesh = f_in.getMesh(); // Convert to field aligned coordinates const bool is_unaligned = (f_in.getDirectionY() == YDirectionType::Standard); @@ -313,10 +313,10 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { Coordinates* coord = f_in.getCoordinates(); for (int i = mesh->xstart; i <= mesh->xend; i++) { - bool yperiodic = mesh->periodicY(i); + const bool yperiodic = mesh->periodicY(i); - bool has_upper_boundary = !yperiodic && mesh->lastY(i); - bool has_lower_boundary = !yperiodic && mesh->firstY(i); + const bool has_upper_boundary = !yperiodic && mesh->lastY(i); + const bool has_lower_boundary = !yperiodic && mesh->firstY(i); for (int j = mesh->ystart; j <= mesh->yend; j++) { @@ -341,8 +341,8 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { // Not on domain boundary // 3rd derivative at right cell boundary - const BoutReal d3fdx3 = - (f(i, j + 2, k) - 3. * f(i, j + 1, k) + 3. * f(i, j, k) - f(i, j - 1, k)); + const BoutReal d3fdx3 = (f(i, j + 2, k) - (3. * f(i, j + 1, k)) + + (3. * f(i, j, k)) - f(i, j - 1, k)); result(i, j, k) += d3fdx3 * factor_rc; result(i, j + 1, k) -= d3fdx3 * factor_rp; @@ -363,10 +363,10 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { common_factor / (coord->J(i, j + 1, k) * coord->dy(i, j + 1, k)); const BoutReal d3fdx3 = - -((16. / 5) * 0.5 * (f(i, j + 1, k) + f(i, j, k)) // Boundary value f_b - - 6. * f(i, j, k) // f_0 - + 4. * f(i, j - 1, k) // f_1 - - (6. / 5) * f(i, j - 2, k) // f_2 + -(((16. / 5) * 0.5 * (f(i, j + 1, k) + f(i, j, k))) // Boundary value f_b + - (6. * f(i, j, k)) // f_0 + + (4. * f(i, j - 1, k)) // f_1 + - ((6. / 5) * f(i, j - 2, k)) // f_2 ); result(i, j, k) += d3fdx3 * factor_rc; @@ -392,8 +392,8 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { common_factor / (coord->J(i, j - 1, k) * coord->dy(i, j - 1, k)); // Not on a domain boundary - const BoutReal d3fdx3 = - (f(i, j + 1, k) - 3. * f(i, j, k) + 3. * f(i, j - 1, k) - f(i, j - 2, k)); + const BoutReal d3fdx3 = (f(i, j + 1, k) - (3. * f(i, j, k)) + + (3. * f(i, j - 1, k)) - f(i, j - 2, k)); result(i, j, k) -= d3fdx3 * factor_lc; result(i, j - 1, k) += d3fdx3 * factor_lm; @@ -410,10 +410,10 @@ const Field3D D4DY4_Index(const Field3D& f_in, bool bndry_flux) { const BoutReal factor_lm = common_factor / (coord->J(i, j - 1, k) * coord->dy(i, j - 1, k)); const BoutReal d3fdx3 = - -(-(16. / 5) * 0.5 * (f(i, j - 1, k) + f(i, j, k)) // Boundary value f_b - + 6. * f(i, j, k) // f_0 - - 4. * f(i, j + 1, k) // f_1 - + (6. / 5) * f(i, j + 2, k) // f_2 + -((-(16. / 5) * 0.5 * (f(i, j - 1, k) + f(i, j, k))) // Boundary value f_b + + (6. * f(i, j, k)) // f_0 + - (4. * f(i, j + 1, k)) // f_1 + + ((6. / 5) * f(i, j + 2, k)) // f_2 ); result(i, j, k) -= d3fdx3 * factor_lc; @@ -436,8 +436,9 @@ void communicateFluxes(Field3D& f) { throw BoutException("communicateFluxes: Sorry!"); } - int size = mesh->LocalNy * mesh->LocalNz; - comm_handle xin, xout; + const int size = mesh->LocalNy * mesh->LocalNz; + comm_handle xin = nullptr; + comm_handle xout = nullptr; // Cache results to silence spurious compiler warning about xin, // xout possibly being uninitialised when used const bool not_first = mesh->periodicX || !mesh->firstX(); @@ -496,45 +497,45 @@ Field3D Div_Perp_Lap(const Field3D& a, const Field3D& f, CELL_LOC outloc) { // o --- gD --- o // Coordinates* coords = a.getCoordinates(outloc); - Mesh* mesh = f.getMesh(); + const Mesh* mesh = f.getMesh(); for (int i = mesh->xstart; i <= mesh->xend; i++) { for (int j = mesh->ystart; j <= mesh->yend; j++) { for (int k = 0; k < mesh->LocalNz; k++) { // wrap k-index around as Z is (currently) periodic. - int kp = (k + 1) % (mesh->LocalNz); - int km = (k - 1 + mesh->LocalNz) % (mesh->LocalNz); + const int kp = (k + 1) % (mesh->LocalNz); + const int km = (k - 1 + mesh->LocalNz) % (mesh->LocalNz); // Calculate gradients on cell faces -- assumes constant grid spacing - BoutReal gR = - (coords->g11(i, j, k) + coords->g11(i + 1, j, k)) - * (f(i + 1, j, k) - f(i, j, k)) - / (coords->dx(i + 1, j, k) + coords->dx(i, j, k)) - + 0.5 * (coords->g13(i, j, k) + coords->g13(i + 1, j, k)) - * (f(i + 1, j, kp) - f(i + 1, j, km) + f(i, j, kp) - f(i, j, km)) - / (4. * coords->dz(i, j, k)); - - BoutReal gL = - (coords->g11(i - 1, j, k) + coords->g11(i, j, k)) - * (f(i, j, k) - f(i - 1, j, k)) - / (coords->dx(i - 1, j, k) + coords->dx(i, j, k)) - + 0.5 * (coords->g13(i - 1, j, k) + coords->g13(i, j, k)) - * (f(i - 1, j, kp) - f(i - 1, j, km) + f(i, j, kp) - f(i, j, km)) - / (4 * coords->dz(i, j, k)); - - BoutReal gD = - coords->g13(i, j, k) - * (f(i + 1, j, km) - f(i - 1, j, km) + f(i + 1, j, k) - f(i - 1, j, k)) - / (4. * coords->dx(i, j, k)) - + coords->g33(i, j, k) * (f(i, j, k) - f(i, j, km)) / coords->dz(i, j, k); - - BoutReal gU = - coords->g13(i, j, k) - * (f(i + 1, j, kp) - f(i - 1, j, kp) + f(i + 1, j, k) - f(i - 1, j, k)) - / (4. * coords->dx(i, j, k)) - + coords->g33(i, j, k) * (f(i, j, kp) - f(i, j, k)) / coords->dz(i, j, k); + const BoutReal gR = + ((coords->g11(i, j, k) + coords->g11(i + 1, j, k)) + * (f(i + 1, j, k) - f(i, j, k)) + / (coords->dx(i + 1, j, k) + coords->dx(i, j, k))) + + (0.5 * (coords->g13(i, j, k) + coords->g13(i + 1, j, k)) + * (f(i + 1, j, kp) - f(i + 1, j, km) + f(i, j, kp) - f(i, j, km)) + / (4. * coords->dz(i, j, k))); + + const BoutReal gL = + ((coords->g11(i - 1, j, k) + coords->g11(i, j, k)) + * (f(i, j, k) - f(i - 1, j, k)) + / (coords->dx(i - 1, j, k) + coords->dx(i, j, k))) + + (0.5 * (coords->g13(i - 1, j, k) + coords->g13(i, j, k)) + * (f(i - 1, j, kp) - f(i - 1, j, km) + f(i, j, kp) - f(i, j, km)) + / (4 * coords->dz(i, j, k))); + + const BoutReal gD = + (coords->g13(i, j, k) + * (f(i + 1, j, km) - f(i - 1, j, km) + f(i + 1, j, k) - f(i - 1, j, k)) + / (4. * coords->dx(i, j, k))) + + (coords->g33(i, j, k) * (f(i, j, k) - f(i, j, km)) / coords->dz(i, j, k)); + + const BoutReal gU = + (coords->g13(i, j, k) + * (f(i + 1, j, kp) - f(i - 1, j, kp) + f(i + 1, j, k) - f(i - 1, j, k)) + / (4. * coords->dx(i, j, k))) + + (coords->g33(i, j, k) * (f(i, j, kp) - f(i, j, k)) / coords->dz(i, j, k)); // Flow right BoutReal flux = gR * 0.25 * (coords->J(i + 1, j, k) + coords->J(i, j, k)) diff --git a/tests/MMS/spatial/finite-volume/fv_mms.cxx b/tests/MMS/spatial/finite-volume/fv_mms.cxx index edf4bbc16a..19f3f14610 100644 --- a/tests/MMS/spatial/finite-volume/fv_mms.cxx +++ b/tests/MMS/spatial/finite-volume/fv_mms.cxx @@ -1,4 +1,5 @@ #include "bout/bout.hxx" +#include "bout/difops.hxx" #include "bout/field.hxx" #include "bout/field3d.hxx" #include "bout/field_factory.hxx" From fcfe963eece7eb1a4c7afe3b0996acf1767b204e Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 12 Nov 2025 09:53:30 +0100 Subject: [PATCH 08/38] Move check up --- include/bout/fv_ops.hxx | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index cd9a3536c1..07a2c2976b 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -600,12 +600,10 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, bool fixflux = true) { Coordinates* coord = f_in.getCoordinates(); + ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); if (f_in.isFci()) { // Use mid-point (cell boundary) averages - if (flow_ylow.isAllocated()) { - flow_ylow = emptyFrom(flow_ylow); - } ASSERT1(f_in.hasParallelSlices()); ASSERT1(v_in.hasParallelSlices()); @@ -631,7 +629,6 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, } return result; } - ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); ASSERT1_FIELDS_COMPATIBLE(f_in, wave_speed_in); const Mesh* mesh = f_in.getMesh(); From b209a5ad52f302f54983ed2b96037faab4b4a2c1 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 20 Jan 2026 09:27:22 +0100 Subject: [PATCH 09/38] Ensure FCI path is always used for FCI --- src/mesh/difops.cxx | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 56773f3c4c..b0b5cc60d5 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -376,7 +376,9 @@ Field3D Div_par_K_Grad_par_mod(const Field3D& Kin, const Field3D& fin, Field3D& const Mesh* mesh = Kin.getMesh(); const Coordinates* coord = fin.getCoordinates(); - if (Kin.hasParallelSlices() && fin.hasParallelSlices()) { + if (Kin.isFci()) { + ASSERT1(Kin.hasParallelSlices()); + ASSERT1(fin.hasParallelSlices()); // Using parallel slices. // Note: Y slices may use different coordinate systems // -> Only B, dy and g_22 can be used in yup/ydown From e0480fc00c367e3c4df08e46326f1b79dd033898 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 4 Feb 2026 10:47:41 +0100 Subject: [PATCH 10/38] FV_div_par_fvv seems to be reduced order due to slope limiting character --- tests/MMS/spatial/fci/runtest | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/tests/MMS/spatial/fci/runtest b/tests/MMS/spatial/fci/runtest index 73babc9691..a08f36fbb8 100755 --- a/tests/MMS/spatial/fci/runtest +++ b/tests/MMS/spatial/fci/runtest @@ -33,6 +33,11 @@ OPERATORS = ( "FV_div_par_mod", "FV_div_par_fvv", ) + +# div_par_fvv is also tested in ../finite-volume, where 1.5th order is achieved +operator_order = { + "FV_div_par_fvv": 1, +} # Note that we need at least _2_ interior points for hermite spline # interpolation due to an awkwardness with the boundaries NX = 4 @@ -172,7 +177,10 @@ def check_fci_operators(name: str, case: dict) -> bool: for operator in OPERATORS: test_name = f"{operator} {name}" success = assert_convergence( - final_errors[operator]["l_2"], dx, test_name, order + final_errors[operator]["l_2"], + dx, + test_name, + operator_order.get(operator, order), ) if not success: failures.append(test_name) From 7ffc7facd5ebef7259c0560220d390f7a7d17af2 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 4 Feb 2026 11:34:25 +0100 Subject: [PATCH 11/38] Add doc to option Co-authored-by: Peter Hill --- src/mesh/interpolation/hermite_spline_xz.cxx | 31 ++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx index 27a4f1d614..1c140b0de2 100644 --- a/src/mesh/interpolation/hermite_spline_xz.cxx +++ b/src/mesh/interpolation/hermite_spline_xz.cxx @@ -107,6 +107,18 @@ XZHermiteSpline::XZHermiteSpline(int y_offset, Mesh* meshin) h10_x(localmesh), h11_x(localmesh), h00_z(localmesh), h01_z(localmesh), h10_z(localmesh), h11_z(localmesh) { + if constexpr (monotonic) { + if (options == nullptr) { + options = &Options::root()["mesh:paralleltransform:xzinterpolation"]; + } + abs_fac_monotonic = (*options)["atol"] + .doc("Absolute tolerance for clipping overshoot") + .withDefault(abs_fac_monotonic); + rel_fac_monotonic = (*options)["rtol"] + .doc("Relative tolerance for clipping overshoot") + .withDefault(rel_fac_monotonic); + } + // Index arrays contain guard cells in order to get subscripts right i_corner.reallocate(localmesh->LocalNx, localmesh->LocalNy, localmesh->LocalNz); k_corner.reallocate(localmesh->LocalNx, localmesh->LocalNy, localmesh->LocalNz); @@ -422,6 +434,25 @@ Field3D XZHermiteSpline::interpolate(const Field3D& f, const std::string& region f_interp[iyp] = +f_z * h00_z[i] + f_zp1 * h01_z[i] + fz_z * h10_z[i] + fz_zp1 * h11_z[i]; + if constexpr (monotonic) { +#endif + const auto corners = {(*gf)[IndG3D(g3dinds[i][0])], (*gf)[IndG3D(g3dinds[i][1])], + (*gf)[IndG3D(g3dinds[i][2])], (*gf)[IndG3D(g3dinds[i][3])]}; + const auto minmax = std::minmax(corners); + + const auto diff = + ((minmax.second - minmax.first) * rel_fac_monotonic) + abs_fac_monotonic; + f_interp[iyp] = std::max(f_interp[iyp], minmax.first - diff); + f_interp[iyp] = std::min(f_interp[iyp], minmax.second + diff); + } +#if USE_NEW_WEIGHTS and defined(HS_USE_PETSC) + ASSERT2(std::isfinite(cptr[int(i)])); + } + VecRestoreArrayRead(result, &cptr); +#elif USE_NEW_WEIGHTS + ASSERT2(std::isfinite(f_interp[iyp])); + } +#else ASSERT2(std::isfinite(f_interp[iyp]) || i.x() < localmesh->xstart || i.x() > localmesh->xend); } From b2a73839911b3f247b1800d7563307667eebfb34 Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 20 Mar 2026 21:01:44 +0100 Subject: [PATCH 12/38] prek fixes --- include/bout/difops.hxx | 28 +++++++++---------- include/bout/fv_ops.hxx | 12 ++++---- src/mesh/interpolation/hermite_spline_xz.cxx | 1 - .../MMS/spatial/finite-volume/CMakeLists.txt | 6 ++-- tests/MMS/spatial/finite-volume/mms.py | 2 +- tests/MMS/spatial/finite-volume/runtest | 9 +----- 6 files changed, 25 insertions(+), 33 deletions(-) diff --git a/include/bout/difops.hxx b/include/bout/difops.hxx index 2cd99f8d33..2070cc30d3 100644 --- a/include/bout/difops.hxx +++ b/include/bout/difops.hxx @@ -1,11 +1,11 @@ /*!****************************************************************************** * \file difops.hxx - * + * * Differential operators * * Changelog: * - * 2009-01 Ben Dudson + * 2009-01 Ben Dudson * * Added two optional parameters which can be put in any order * These determine the method to use (DIFF_METHOD) * and CELL_LOC location of the result. @@ -15,7 +15,7 @@ * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu * * Contact: Ben Dudson, bd512@york.ac.uk - * + * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify @@ -30,7 +30,7 @@ * * You should have received a copy of the GNU Lesser General Public License * along with BOUT++. If not, see . - * + * *******************************************************************************/ #ifndef BOUT_DIFOPS_H @@ -81,7 +81,7 @@ Field3D Grad_parP(const Field3D& apar, const Field3D& f); * \f[ * v\mathbf{b}_0 \cdot \nabla f * \f] - * + * * * @param[in] v The velocity in y direction * @param[in] f The scalar field to be differentiated @@ -175,7 +175,7 @@ inline Field3D Grad2_par2(const Field3D& f, CELL_LOC outloc, DIFF_METHOD method) /*! * Parallel divergence of diffusive flux, K*Grad_par - * + * * \f[ * \nabla \cdot ( \mathbf{b}_0 kY (\mathbf{b}_0 \cdot \nabla) f ) * \f] @@ -203,8 +203,8 @@ Field3D Div_par_K_Grad_par_mod(const Field3D& k, const Field3D& f, Field3D& flow * Perpendicular Laplacian operator * * This version only includes terms in X and Z, dropping - * derivatives in Y. This is the inverse operation to - * the Laplacian inversion class. + * derivatives in Y. This is the inverse operation to + * the Laplacian inversion class. * * For the full perpendicular Laplacian, use Laplace_perp */ @@ -216,7 +216,7 @@ FieldPerp Delp2(const FieldPerp& f, CELL_LOC outloc = CELL_DEFAULT, bool useFFT /*! * Perpendicular Laplacian, keeping y derivatives * - * + * */ Coordinates::FieldMetric Laplace_perp(const Field2D& f, CELL_LOC outloc = CELL_DEFAULT, @@ -250,18 +250,18 @@ Field2D Laplace_perpXY(const Field2D& A, const Field2D& f); /*! * Terms of form b0 x Grad(phi) dot Grad(A) - * + * */ Coordinates::FieldMetric b0xGrad_dot_Grad(const Field2D& phi, const Field2D& A, CELL_LOC outloc = CELL_DEFAULT); /*! - * Terms of form + * Terms of form * * \f[ * \mathbf{b}_0 \times \nabla \phi \cdot \nabla A * \f] - * + * * @param[in] phi The scalar potential * @param[in] A The field being advected * @param[in] outloc The cell location where the result is defined. By default the same as A. @@ -295,13 +295,13 @@ constexpr BRACKET_METHOD BRACKET_CTU = BRACKET_METHOD::ctu; * \f[ * [f, g] = (1/B) \mathbf{b}_0 \times \nabla f \cdot \nabla g * \f] - * + * * @param[in] f The potential * @param[in] g The field being advected * @param[in] method The method to use * @param[in] outloc The cell location where the result is defined. Default is the same as g * @param[in] solver Pointer to the time integration solver - * + * */ Coordinates::FieldMetric bracket(const Field2D& f, const Field2D& g, BRACKET_METHOD method = BRACKET_STD, diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 07a2c2976b..027a1ef77b 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -115,7 +115,7 @@ struct Fromm { /*! * Second order slope limiter method - * + * * Limits slope to minimum absolute value * of left and right gradients. If at a maximum * or minimum slope set to zero, i.e. reverts @@ -152,8 +152,8 @@ private: /*! * Monotonised Central (MC) second order slope limiter (Van Leer) - * - * Limits the slope based on taking the slope with + * + * Limits the slope based on taking the slope with * the minimum absolute value from central, 2*left and * 2*right. If any of these slopes have different signs * then the slope reverts to zero (i.e. 1st-order upwinding). @@ -256,7 +256,7 @@ Field3D Div_par(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_sp ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); ASSERT1_FIELDS_COMPATIBLE(f_in, wave_speed_in); - Mesh const* mesh = f_in.getMesh(); + const Mesh* mesh = f_in.getMesh(); CellEdges cellboundary; @@ -423,9 +423,9 @@ Field3D Div_par(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_sp * Div ( n * v ) -- Magnetic drifts * * This uses the expression - * + * * Div( A ) = 1/J * d/di ( J * A^i ) - * + * * Hence the input vector should be contravariant * * Note: Uses to/from FieldAligned diff --git a/src/mesh/interpolation/hermite_spline_xz.cxx b/src/mesh/interpolation/hermite_spline_xz.cxx index e2760ac0e3..5b9271c93f 100644 --- a/src/mesh/interpolation/hermite_spline_xz.cxx +++ b/src/mesh/interpolation/hermite_spline_xz.cxx @@ -509,7 +509,6 @@ Field3D XZHermiteSplineBase::interpolate( f_interp[iyp] = +f_z * h00_z[i] + f_zp1 * h01_z[i] + fz_z * h10_z[i] + fz_zp1 * h11_z[i]; - if constexpr (monotonic) { const auto corners = {(*gf)[IndG3D(g3dinds[i][0])], (*gf)[IndG3D(g3dinds[i][1])], (*gf)[IndG3D(g3dinds[i][2])], (*gf)[IndG3D(g3dinds[i][3])]}; diff --git a/tests/MMS/spatial/finite-volume/CMakeLists.txt b/tests/MMS/spatial/finite-volume/CMakeLists.txt index 6d9c839a05..7180eb7d98 100644 --- a/tests/MMS/spatial/finite-volume/CMakeLists.txt +++ b/tests/MMS/spatial/finite-volume/CMakeLists.txt @@ -1,6 +1,6 @@ -bout_add_mms_test(MMS-spatial-finite-volume +bout_add_mms_test( + MMS-spatial-finite-volume SOURCES fv_mms.cxx - USE_RUNTEST - USE_DATA_BOUT_INP + USE_RUNTEST USE_DATA_BOUT_INP PROCESSORS 2 ) diff --git a/tests/MMS/spatial/finite-volume/mms.py b/tests/MMS/spatial/finite-volume/mms.py index dfcfce9a09..a95ecc1328 100755 --- a/tests/MMS/spatial/finite-volume/mms.py +++ b/tests/MMS/spatial/finite-volume/mms.py @@ -32,7 +32,7 @@ v = fv / f # Substitute back to get input y coordinates -replace = [ (metric.y, y*Ly/(2*pi) ) ] +replace = [(metric.y, y * Ly / (2 * pi))] def Grad2_par2(f: Expr) -> Expr: diff --git a/tests/MMS/spatial/finite-volume/runtest b/tests/MMS/spatial/finite-volume/runtest index b38a6359ac..bcd4672545 100755 --- a/tests/MMS/spatial/finite-volume/runtest +++ b/tests/MMS/spatial/finite-volume/runtest @@ -22,23 +22,18 @@ OPERATORS = { "FV_Div_par_MC": 1.5, "FV_Div_par_mod_MC": 1.5, "FV_Div_par_fvv_MC": 1.5, - "FV_Div_par_Upwind": 1, "FV_Div_par_mod_Upwind": 1, "FV_Div_par_fvv_Upwind": 1, - "FV_Div_par_Fromm": 1.5, "FV_Div_par_mod_Fromm": 1.5, "FV_Div_par_fvv_Fromm": 1.5, - "FV_Div_par_MinMod": 1.5, "FV_Div_par_mod_MinMod": 1.5, "FV_Div_par_fvv_MinMod": 1.5, - "FV_Div_par_Superbee": 1.5, "FV_Div_par_mod_Superbee": 1.5, "FV_Div_par_fvv_Superbee": 1.5, - "FV_Div_par_K_Grad_par": 2, "FV_Div_par_K_Grad_par_mod": 2, } @@ -131,9 +126,7 @@ def test_fv_operators() -> bool: final_errors = transpose(all_errors) for operator, order in OPERATORS.items(): - success = assert_convergence( - final_errors[operator]["l_2"], dx, operator, order - ) + success = assert_convergence(final_errors[operator]["l_2"], dx, operator, order) if not success: failures.append(operator) From 24c99925c6a368134ba25fd1c31dcd913ca66d10 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 15 Jun 2026 14:22:24 +0200 Subject: [PATCH 13/38] Remove unused header --- src/mesh/difops.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index b0b5cc60d5..f23e1ed01b 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -38,7 +38,6 @@ #include "bout/solver.hxx" #include "bout/unused.hxx" #include "bout/utils.hxx" -#include "bout/vecops.hxx" #include From 639808c795879e6948a305c6b99631f00c2ba2ff Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 16 Jun 2026 14:08:32 +0200 Subject: [PATCH 14/38] Use cell area and volumes for div_par_mod --- include/bout/fv_ops.hxx | 59 +++++++++++------------------------------ 1 file changed, 15 insertions(+), 44 deletions(-) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 3a1a4eaee8..9ffc235918 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -710,12 +710,10 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, const auto iym = i.ym(); result[i] = (0.25 * (f_in[i] + f_up[iyp]) * (v_in[i] + v_up[iyp]) - * (coord->J[i] + coord->J.yup()[iyp]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22.yup()[iyp])) + * coord->cell_area_yhigh()[i] - 0.25 * (f_in[i] + f_down[iym]) * (v_in[i] + v_down[iym]) - * (coord->J[i] + coord->J.ydown()[iym]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22.ydown()[iym]))) - / (coord->dy[i] * coord->J[i]); + * coord->cell_area_ylow()[i]) + / coord->cell_volume()[i]; } return result; } @@ -756,57 +754,30 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, // Pre-calculate factors which multiply fluxes #if not(BOUT_USE_METRIC_3D) // For right cell boundaries - const BoutReal common_factor_r = - (coord->J(i, j) + coord->J(i, j + 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1))); + const BoutReal area_rp = coord->cell_area_yhigh()(i, j); - const BoutReal flux_factor_rc = - common_factor_r / (coord->dy(i, j) * coord->J(i, j)); - const BoutReal flux_factor_rp = - common_factor_r / (coord->dy(i, j + 1) * coord->J(i, j + 1)); - - const BoutReal area_rp = - common_factor_r * coord->dx(i, j + 1) * coord->dz(i, j + 1); + const BoutReal flux_factor_rc = area_rp / coord->cell_volume()(i, j); + const BoutReal flux_factor_rp = area_rp / coord->cell_volume()(i, j + 1); // For left cell boundaries - const BoutReal common_factor_l = - (coord->J(i, j) + coord->J(i, j - 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1))); + const BoutReal area_lc = coord->cell_area_ylow()(i, j); - const BoutReal flux_factor_lc = - common_factor_l / (coord->dy(i, j) * coord->J(i, j)); - const BoutReal flux_factor_lm = - common_factor_l / (coord->dy(i, j - 1) * coord->J(i, j - 1)); - - const BoutReal area_lc = common_factor_l * coord->dx(i, j) * coord->dz(i, j); + const BoutReal flux_factor_lc = area_lc / coord->cell_volume()(i, j); + const BoutReal flux_factor_lm = area_lc / coord->cell_volume()(i, j - 1); #endif for (int k = 0; k < mesh->LocalNz; k++) { #if BOUT_USE_METRIC_3D // For right cell boundaries - const BoutReal common_factor_r = - (coord->J(i, j, k) + coord->J(i, j + 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); - - const BoutReal flux_factor_rc = - common_factor_r / (coord->dy(i, j, k) * coord->J(i, j, k)); - const BoutReal flux_factor_rp = - common_factor_r / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); + const BoutReal area_rp = coord->cell_area_yhigh()(i, j, k); - const BoutReal area_rp = - common_factor_r * coord->dx(i, j + 1, k) * coord->dz(i, j + 1, k); + const BoutReal flux_factor_rc = area_rp / coord->cell_volume()(i, j, k); + const BoutReal flux_factor_rp = area_rp / coord->cell_volume()(i, j + 1, k); // For left cell boundaries - const BoutReal common_factor_l = - (coord->J(i, j, k) + coord->J(i, j - 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); - - const BoutReal flux_factor_lc = - common_factor_l / (coord->dy(i, j, k) * coord->J(i, j, k)); - const BoutReal flux_factor_lm = - common_factor_l / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); + const BoutReal area_lc = coord->cell_area_ylow()(i, j, k); - const BoutReal area_lc = - common_factor_l * coord->dx(i, j, k) * coord->dz(i, j, k); + const BoutReal flux_factor_lc = area_lc / coord->cell_volume()(i, j, k); + const BoutReal flux_factor_lm = area_lc / coord->cell_volume()(i, j - 1, k); #endif //////////////////////////////////////////// From 5a2a3dc0e5e19a16c3e3e652717cb917f09a13d7 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 16 Jun 2026 14:11:43 +0200 Subject: [PATCH 15/38] Use cell area and volumes for Div_par_fvv --- include/bout/fv_ops.hxx | 22 +++++++--------------- 1 file changed, 7 insertions(+), 15 deletions(-) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 9ffc235918..6198ca81d6 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -974,24 +974,16 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, for (int k = 0; k < mesh->LocalNz; k++) { // For right cell boundaries - const BoutReal common_factor_r = - (coord->J(i, j, k) + coord->J(i, j + 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); + const BoutReal area_r = coord->cell_area_yhigh()(i, j, k); - const BoutReal flux_factor_rc = - common_factor_r / (coord->dy(i, j, k) * coord->J(i, j, k)); - const BoutReal flux_factor_rp = - common_factor_r / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); + const BoutReal flux_factor_rc = area_r / coord->cell_volume()(i, j, k); + const BoutReal flux_factor_rp = area_r / coord->cell_volume()(i, j + 1, k); // For left cell boundaries - const BoutReal common_factor_l = - (coord->J(i, j, k) + coord->J(i, j - 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); - - const BoutReal flux_factor_lc = - common_factor_l / (coord->dy(i, j, k) * coord->J(i, j, k)); - const BoutReal flux_factor_lm = - common_factor_l / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); + const BoutReal area_l = coord->cell_area_ylow()(i, j, k); + + const BoutReal flux_factor_lc = area_l / coord->cell_volume()(i, j, k); + const BoutReal flux_factor_lm = area_l / coord->cell_volume()(i, j - 1, k); //////////////////////////////////////////// // Reconstruct f at the cell faces From 820e2219e365ee8250eda82d5ae403e712e6abc6 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 16 Jun 2026 14:17:06 +0200 Subject: [PATCH 16/38] Use cell area and volumes for Div_par --- include/bout/fv_ops.hxx | 41 ++++++++++++----------------------------- 1 file changed, 12 insertions(+), 29 deletions(-) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 6198ca81d6..e96dfca1d7 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -382,41 +382,24 @@ Field3D Div_par(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_sp // Pre-calculate factors which multiply fluxes #if not(BOUT_USE_METRIC_3D) // For right cell boundaries - BoutReal common_factor = (coord->J(i, j) + coord->J(i, j + 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1))); - - const BoutReal flux_factor_rc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - const BoutReal flux_factor_rp = - common_factor / (coord->dy(i, j + 1) * coord->J(i, j + 1)); - + const BoutReal area_r = coord->cell_area_yhigh()(i, j); + const BoutReal flux_factor_rc = area_r / coord->cell_volume()(i, j); + const BoutReal flux_factor_rp = area_r / coord->cell_volume()(i, j + 1); // For left cell boundaries - common_factor = (coord->J(i, j) + coord->J(i, j - 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1))); - - const BoutReal flux_factor_lc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - const BoutReal flux_factor_lm = - common_factor / (coord->dy(i, j - 1) * coord->J(i, j - 1)); + const BoutReal area_l = coord->cell_area_ylow()(i, j); + const BoutReal flux_factor_lc = area_l / coord->cell_volume()(i, j); + const BoutReal flux_factor_lm = area_l / coord->cell_volume()(i, j - 1); #endif for (int k = mesh->zstart; k <= mesh->zend; k++) { #if BOUT_USE_METRIC_3D // For right cell boundaries - BoutReal common_factor = - (coord->J(i, j, k) + coord->J(i, j + 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); - - BoutReal flux_factor_rc = - common_factor / (coord->dy(i, j, k) * coord->J(i, j, k)); - BoutReal flux_factor_rp = - common_factor / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); - + const BoutReal area_r = coord->cell_area_yhigh()(i, j, k); + const BoutReal flux_factor_rc = area_r / coord->cell_volume()(i, j, k); + const BoutReal flux_factor_rp = area_r / coord->cell_volume()(i, j + 1, k); // For left cell boundaries - common_factor = (coord->J(i, j, k) + coord->J(i, j - 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); - - BoutReal flux_factor_lc = - common_factor / (coord->dy(i, j, k) * coord->J(i, j, k)); - BoutReal flux_factor_lm = - common_factor / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); + const BoutReal area_l = coord->cell_area_ylow()(i, j, k); + const BoutReal flux_factor_lc = area_l / coord->cell_volume()(i, j, k); + const BoutReal flux_factor_lm = area_l / coord->cell_volume()(i, j - 1, k); #endif //////////////////////////////////////////// From 3b0117ea11652b576fbcf8e2dcb2be9163227ce3 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 17 Jun 2026 10:46:58 +0200 Subject: [PATCH 17/38] Also fill the boundaries to the extend possible --- src/mesh/coordinates.cxx | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 39bcf16936..2e1c4ad389 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -2099,10 +2099,17 @@ void Coordinates::_compute_cell_area_y() const { _cell_area_yhigh.emplace(emptyFrom(area_centre)); // We cannot setLocation, as that would trigger the computation of staggered // metrics. - ASSERT0(mesh->ystart > 0); - BOUT_FOR(i, mesh->getRegion("RGN_NOY")) { - (*_cell_area_ylow)[i] = 0.5 * (area_centre[i] + area_centre[i.ym()]); - (*_cell_area_yhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.yp()]); + BOUT_FOR(i, mesh->getRegion("RGN_ALL")) { + if (i.y() > 0) { + (*_cell_area_ylow)[i] = 0.5 * (area_centre[i] + area_centre[i.ym()]); + } else { + (*_cell_area_ylow)[i] = BoutNaN; + } + if (i.y() < mesh->LocalNy - 1) { + (*_cell_area_yhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.yp()]); + } else { + (*_cell_area_yhigh)[i] = BoutNaN; + } } } } From 2e5367a60f20d820b5b758854a03cea28168feb4 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 23 Jun 2026 14:30:10 +0200 Subject: [PATCH 18/38] Move implementations to different file --- include/bout/fv_ops.hxx | 971 +----------------------- src/mesh/fv_ops.cxx | 69 ++ tests/unit/include/bout/test_fv_ops.cxx | 2 +- 3 files changed, 84 insertions(+), 958 deletions(-) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 3a1a4eaee8..d256bad32a 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -77,245 +77,14 @@ Field3D D4DY4(const Field3D& d, const Field3D& f); */ Field3D D4DY4_Index(const Field3D& f, bool bndry_flux = true); -/*! - * Stencil used for Finite Volume calculations - * which includes cell face values L and R - */ -struct Stencil1D { - /// Cell centre values - BoutReal c; - BoutReal m; - BoutReal p; - BoutReal mm = BoutNaN; - BoutReal pp = BoutNaN; - - /// Left cell face value - BoutReal L = BoutNaN; - /// Right cell face value - BoutReal R = BoutNaN; -}; - -/*! - * First order upwind for testing - */ -struct Upwind { - void operator()(Stencil1D& n) { n.L = n.R = n.c; } -}; - -/*! - * Fromm method - */ -struct Fromm { - void operator()(Stencil1D& n) { - n.L = n.c - (0.25 * (n.p - n.m)); - n.R = n.c + (0.25 * (n.p - n.m)); - } -}; - -/*! - * Second order slope limiter method - * - * Limits slope to minimum absolute value - * of left and right gradients. If at a maximum - * or minimum slope set to zero, i.e. reverts - * to first order upwinding - */ -struct MinMod { - void operator()(Stencil1D& n) { - // Choose the gradient within the cell - // as the minimum (smoothest) solution - const BoutReal slope = _minmod(n.p - n.c, n.c - n.m); - n.L = n.c - (0.5 * slope); - n.R = n.c + (0.5 * slope); - } - -private: - /*! - * Internal helper function for minmod slope limiter - * - * If the inputs have different signs then - * returns zero, otherwise chooses the value - * with the minimum magnitude. - */ - static BoutReal _minmod(BoutReal a, BoutReal b) { - if (a * b <= 0.0) { - return 0.0; - } - - if (fabs(a) < fabs(b)) { - return a; - } - return b; - } -}; - -/*! - * Monotonised Central (MC) second order slope limiter (Van Leer) - * - * Limits the slope based on taking the slope with - * the minimum absolute value from central, 2*left and - * 2*right. If any of these slopes have different signs - * then the slope reverts to zero (i.e. 1st-order upwinding). - */ -struct MC { - void operator()(Stencil1D& n) { - const BoutReal slope = minmod(2. * (n.p - n.c), // 2*right difference - 0.5 * (n.p - n.m), // Central difference - 2. * (n.c - n.m)); // 2*left difference - n.L = n.c - (0.5 * slope); - n.R = n.c + (0.5 * slope); - } - -private: - // Return zero if any signs are different - // otherwise return the value with the minimum magnitude - static BoutReal minmod(BoutReal a, BoutReal b, BoutReal c) { - // if any of the signs are different, return zero gradient - if ((a * b <= 0.0) || (a * c <= 0.0)) { - return 0.0; - } - - // Return the minimum absolute value - return SIGN(a) * BOUTMIN(fabs(a), fabs(b), fabs(c)); - } -}; - -/// Superbee limiter -/// -/// This corresponds to the limiter function -/// φ(r) = max(0, min(2r, 1), min(r,2) -/// -/// The value at cell right (i.e. i + 1/2) is: -/// -/// n.R = n.c - φ(r) (n.c - (n.p + n.c)/2) -/// = n.c + φ(r) (n.p - n.c)/2 -/// -/// Four regimes: -/// a) r < 1/2 -> φ(r) = 2r -/// n.R = n.c + gL -/// b) 1/2 < r < 1 -> φ(r) = 1 -/// n.R = n.c + gR/2 -/// c) 1 < r < 2 -> φ(r) = r -/// n.R = n.c + gL/2 -/// d) 2 < r -> φ(r) = 2 -/// n.R = n.c + gR -/// -/// where the left and right gradients are: -/// gL = n.c - n.m -/// gR = n.p - n.c -/// -struct Superbee { - void operator()(Stencil1D& n) { - const BoutReal gL = n.c - n.m; - const BoutReal gR = n.p - n.c; - - // r = gL / gR - // Limiter is φ(r) - if (gL * gR < 0) { - // Different signs => Zero gradient - n.L = n.R = n.c; - } else { - const BoutReal sign = SIGN(gL); - const BoutReal abs_gL = fabs(gL); - const BoutReal abs_gR = fabs(gR); - const BoutReal half_slope = - sign * BOUTMAX(BOUTMIN(abs_gL, 0.5 * abs_gR), BOUTMIN(abs_gR, 0.5 * abs_gL)); - n.L = n.c - half_slope; - n.R = n.c + half_slope; - } - } -}; - -/*! - * Symmetric Van Albada second order slope limiter - * - * Uses a smooth (differentiable) approximation to `max(a*b, 0)` to avoid - * introducing a kink at extrema, which can be helpful for nonlinear solvers - * and finite-difference Jacobian calculations. - * - * The limited slope is calculated from the left and right differences - * `dl = c - m` and `dr = p - c` as - * - * slope = (pos(dl*dr) * (dl + dr)) / (dl^2 + dr^2) - * - * where `pos(x)` is a smooth approximation to `max(x, 0)`. - */ -struct VanAlbada { - void operator()(Stencil1D& n) { - const BoutReal dl = n.c - n.m; - const BoutReal dr = n.p - n.c; - - const BoutReal denom = dl * dl + dr * dr; - - // Smoothness parameters: - // - keep division well-defined when dl=dr=0 - // - provide a differentiable approximation to max(dl*dr, 0) - const BoutReal eps = 1e-12 * denom + 1e-30; - - const BoutReal ab = dl * dr; - const BoutReal ab_pos = 0.5 * (ab + sqrt(ab * ab + eps * eps)); - - const BoutReal slope = (ab_pos * (dl + dr)) / (denom + eps); - - n.L = n.c - 0.5 * slope; - n.R = n.c + 0.5 * slope; - } -}; - -/*! - * WENO3-JS (Jiang-Shu) reconstruction to cell faces - * - * This is a third-order essentially non-oscillatory reconstruction using two - * candidate second-order polynomials and smoothness-weighted blending. - * - * Unlike TVD slope limiters (e.g. ``MC``), WENO reconstruction is generally - * smooth (differentiable) for all inputs, but it does not enforce strict - * monotonicity. - * - * Uses only the three-point stencil (`m`, `c`, `p`), so it is a drop-in - * replacement anywhere `Stencil1D` is populated with those values. - */ -struct WENO3 { - void operator()(Stencil1D& n) { - // Right face (between c and p): value from cell c (left state at i+1/2) - const BoutReal p0_r = 0.5 * (-n.m + 3.0 * n.c); - const BoutReal p1_r = 0.5 * (n.c + n.p); - - const BoutReal beta0_r = SQ(n.c - n.m); - const BoutReal beta1_r = SQ(n.p - n.c); - - // Left face (between m and c): value from cell c (right state at i-1/2) - const BoutReal p0_l = 0.5 * (-n.p + 3.0 * n.c); - const BoutReal p1_l = 0.5 * (n.m + n.c); - - const BoutReal beta0_l = beta1_r; - const BoutReal beta1_l = beta0_r; - - // Smoothness parameter (scaled to local variation) - const BoutReal eps = 1e-12 * (beta0_r + beta1_r) + 1e-30; - - // Linear weights for WENO3-JS - constexpr BoutReal d0 = 1.0 / 3.0; - constexpr BoutReal d1 = 2.0 / 3.0; - - // Right face weights - const BoutReal a0_r = d0 / SQ(eps + beta0_r); - const BoutReal a1_r = d1 / SQ(eps + beta1_r); - const BoutReal wsum_r = a0_r + a1_r; - const BoutReal w0_r = a0_r / wsum_r; - const BoutReal w1_r = a1_r / wsum_r; - - // Left face weights (mirrored) - const BoutReal a0_l = d0 / SQ(eps + beta0_l); - const BoutReal a1_l = d1 / SQ(eps + beta1_l); - const BoutReal wsum_l = a0_l + a1_l; - const BoutReal w0_l = a0_l / wsum_l; - const BoutReal w1_l = a1_l / wsum_l; - - n.R = w0_r * p0_r + w1_r * p1_r; - n.L = w0_l * p0_l + w1_l * p1_l; - } -}; +// FluxLimiter +class Upwind; +class Fromm; +class MinMod; +class MC; +class Superbee; +class VanAlbada; +class WENO3; /*! * Communicate fluxes between processors @@ -341,173 +110,7 @@ void communicateFluxes(Field3D& f); /// NB: Uses to/from FieldAligned coordinates template Field3D Div_par(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, - bool fixflux = true) { - - ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); - ASSERT1_FIELDS_COMPATIBLE(f_in, wave_speed_in); - - const Mesh* mesh = f_in.getMesh(); - - CellEdges cellboundary; - - ASSERT2(f_in.getDirectionY() == v_in.getDirectionY()); - ASSERT2(f_in.getDirectionY() == wave_speed_in.getDirectionY()); - const bool are_unaligned = - ((f_in.getDirectionY() == YDirectionType::Standard) - and (v_in.getDirectionY() == YDirectionType::Standard) - and (wave_speed_in.getDirectionY() == YDirectionType::Standard)); - - Field3D f = are_unaligned ? toFieldAligned(f_in, "RGN_NOX") : f_in; - Field3D v = are_unaligned ? toFieldAligned(v_in, "RGN_NOX") : v_in; - Field3D wave_speed = - are_unaligned ? toFieldAligned(wave_speed_in, "RGN_NOX") : wave_speed_in; - - Coordinates* coord = f_in.getCoordinates(); - - Field3D result{zeroFrom(f)}; - - for (int i = mesh->xstart; i <= mesh->xend; i++) { - const bool is_periodic_y = mesh->periodicY(i); - const bool is_first_y = mesh->firstY(i); - const bool is_last_y = mesh->lastY(i); - - // Only need one guard cell, so no need to communicate fluxes Instead - // calculate in guard cells to get fluxes consistent between processors, but - // don't include the boundary cell. Note that this implies special handling - // of boundaries later - const int ys = (!is_first_y || is_periodic_y) ? mesh->ystart - 1 : mesh->ystart; - const int ye = (!is_last_y || is_periodic_y) ? mesh->yend + 1 : mesh->yend; - - for (int j = ys; j <= ye; j++) { - // Pre-calculate factors which multiply fluxes -#if not(BOUT_USE_METRIC_3D) - // For right cell boundaries - BoutReal common_factor = (coord->J(i, j) + coord->J(i, j + 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1))); - - const BoutReal flux_factor_rc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - const BoutReal flux_factor_rp = - common_factor / (coord->dy(i, j + 1) * coord->J(i, j + 1)); - - // For left cell boundaries - common_factor = (coord->J(i, j) + coord->J(i, j - 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1))); - - const BoutReal flux_factor_lc = common_factor / (coord->dy(i, j) * coord->J(i, j)); - const BoutReal flux_factor_lm = - common_factor / (coord->dy(i, j - 1) * coord->J(i, j - 1)); -#endif - for (int k = mesh->zstart; k <= mesh->zend; k++) { -#if BOUT_USE_METRIC_3D - // For right cell boundaries - BoutReal common_factor = - (coord->J(i, j, k) + coord->J(i, j + 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); - - BoutReal flux_factor_rc = - common_factor / (coord->dy(i, j, k) * coord->J(i, j, k)); - BoutReal flux_factor_rp = - common_factor / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); - - // For left cell boundaries - common_factor = (coord->J(i, j, k) + coord->J(i, j - 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); - - BoutReal flux_factor_lc = - common_factor / (coord->dy(i, j, k) * coord->J(i, j, k)); - BoutReal flux_factor_lm = - common_factor / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); -#endif - - //////////////////////////////////////////// - // Reconstruct f at the cell faces - // This calculates s.R and s.L for the Right and Left - // face values on this cell - - // Reconstruct f at the cell faces - Stencil1D s; - s.c = f(i, j, k); - s.m = f(i, j - 1, k); - s.p = f(i, j + 1, k); - - cellboundary(s); // Calculate s.R and s.L - - //////////////////////////////////////////// - // Right boundary - - // Calculate velocity at right boundary (y+1/2) - BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j + 1, k)); - BoutReal flux = NAN; - - if (is_last_y && (j == mesh->yend) && !is_periodic_y) { - // Last point in domain - - const BoutReal bndryval = 0.5 * (s.c + s.p); - if (fixflux) { - // Use mid-point to be consistent with boundary conditions - flux = bndryval * vpar; - } else { - // Add flux due to difference in boundary values - flux = (s.R * vpar) + (wave_speed(i, j, k) * (s.R - bndryval)); - } - } else { - - // Maximum wave speed in the two cells - const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k)); - - if (vpar > amax) { - // Supersonic flow out of this cell - flux = s.R * vpar; - } else if (vpar < -amax) { - // Supersonic flow into this cell - flux = 0.0; - } else { - // Subsonic flow, so a mix of right and left fluxes - flux = s.R * 0.5 * (vpar + amax); - } - } - - result(i, j, k) += flux * flux_factor_rc; - result(i, j + 1, k) -= flux * flux_factor_rp; - - //////////////////////////////////////////// - // Calculate at left boundary - - vpar = 0.5 * (v(i, j, k) + v(i, j - 1, k)); - - if (is_first_y && (j == mesh->ystart) && !is_periodic_y) { - // First point in domain - const BoutReal bndryval = 0.5 * (s.c + s.m); - if (fixflux) { - // Use mid-point to be consistent with boundary conditions - flux = bndryval * vpar; - } else { - // Add flux due to difference in boundary values - flux = (s.L * vpar) - (wave_speed(i, j, k) * (s.L - bndryval)); - } - } else { - - // Maximum wave speed in the two cells - const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k)); - - if (vpar < -amax) { - // Supersonic out of this cell - flux = s.L * vpar; - } else if (vpar > amax) { - // Supersonic into this cell - flux = 0.0; - } else { - flux = s.L * 0.5 * (vpar - amax); - } - } - - result(i, j, k) -= flux * flux_factor_lc; - result(i, j - 1, k) += flux * flux_factor_lm; - } - } - } - return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result; -} + bool fixflux = true); /*! * Div ( n * v ) -- Magnetic drifts @@ -522,142 +125,7 @@ Field3D Div_par(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_sp * */ template -Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { - ASSERT1(n_in.getLocation() == v.getLocation()); - ASSERT1_FIELDS_COMPATIBLE(n_in, v.x); - - const Mesh* mesh = n_in.getMesh(); - - CellEdges cellboundary; - - Coordinates* coord = n_in.getCoordinates(); - - if (v.covariant) { - // Got a covariant vector instead - throw BoutException("Div_f_v passed a covariant v"); - } - - Field3D result{zeroFrom(n_in)}; - - Field3D vx = v.x; - Field3D vz = v.z; - Field3D n = n_in; - - BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { - // Calculate velocities - const BoutReal vU = 0.25 * (vz[i.zp()] + vz[i]) * (coord->J[i.zp()] + coord->J[i]); - const BoutReal vD = 0.25 * (vz[i.zm()] + vz[i]) * (coord->J[i.zm()] + coord->J[i]); - const BoutReal vL = 0.25 * (vx[i.xm()] + vx[i]) * (coord->J[i.xm()] + coord->J[i]); - const BoutReal vR = 0.25 * (vx[i.xp()] + vx[i]) * (coord->J[i.xp()] + coord->J[i]); - - // X direction - Stencil1D s; - s.c = n[i]; - s.m = n[i.xm()]; - s.mm = n[i.xmm()]; - s.p = n[i.xp()]; - s.pp = n[i.xpp()]; - - cellboundary(s); - - if ((i.x() == mesh->xend) && (mesh->lastX())) { - // At right boundary in X - if (bndry_flux) { - BoutReal flux = NAN; - if (vR > 0.0) { - // Flux to boundary - flux = vR * s.R; - } else { - // Flux in from boundary - flux = vR * 0.5 * (n[i.xp()] + n[i]); - } - result[i] += flux / (coord->dx[i] * coord->J[i]); - result[i.xp()] -= flux / (coord->dx[i.xp()] * coord->J[i.xp()]); - } - } else { - // Not at a boundary - if (vR > 0.0) { - // Flux out into next cell - const BoutReal flux = vR * s.R; - result[i] += flux / (coord->dx[i] * coord->J[i]); - result[i.xp()] -= flux / (coord->dx[i.xp()] * coord->J[i.xp()]); - } - } - - // Left side - - if ((i.x() == mesh->xstart) && (mesh->firstX())) { - // At left boundary in X - - if (bndry_flux) { - BoutReal flux = NAN; - if (vL < 0.0) { - // Flux to boundary - flux = vL * s.L; - } else { - // Flux in from boundary - flux = vL * 0.5 * (n[i.xm()] + n[i]); - } - result[i] -= flux / (coord->dx[i] * coord->J[i]); - result[i.xm()] += flux / (coord->dx[i.xm()] * coord->J[i.xm()]); - } - } else { - // Not at a boundary - if (vL < 0.0) { - const BoutReal flux = vL * s.L; - result[i] -= flux / (coord->dx[i] * coord->J[i]); - result[i.xm()] += flux / (coord->dx[i.xm()] * coord->J[i.xm()]); - } - } - - /// NOTE: Need to communicate fluxes - - // Z direction - s.m = n[i.zm()]; - s.mm = n[i.zmm()]; - s.p = n[i.zp()]; - s.pp = n[i.zpp()]; - - cellboundary(s); - - if (vU > 0.0) { - const BoutReal flux = vU * s.R; - result[i] += flux / (coord->J[i] * coord->dz[i]); - result[i.zp()] -= flux / (coord->J[i.zp()] * coord->dz[i.zp()]); - } - if (vD < 0.0) { - const BoutReal flux = vD * s.L; - result[i] -= flux / (coord->J[i] * coord->dz[i]); - result[i.zm()] += flux / (coord->J[i.zm()] * coord->dz[i.zm()]); - } - } - - communicateFluxes(result); - - // Y advection - // Currently just using simple centered differences - // so no fluxes need to be exchanged - - n = toFieldAligned(n_in, "RGN_NOX"); - Field3D vy = toFieldAligned(v.y, "RGN_NOX"); - - Field3D yresult = 0.0; - yresult.setDirectionY(YDirectionType::Aligned); - - BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { - // Y velocities on y boundaries - const BoutReal vU = 0.25 * (vy[i] + vy[i.yp()]) * (coord->J[i] + coord->J[i.yp()]); - const BoutReal vD = 0.25 * (vy[i] + vy[i.ym()]) * (coord->J[i] + coord->J[i.ym()]); - - // n (advected quantity) on y boundaries - // Note: Use unshifted n_in variable - const BoutReal nU = 0.5 * (n[i] + n[i.yp()]); - const BoutReal nD = 0.5 * (n[i] + n[i.ym()]); - - yresult[i] = (nU * vU - nD * vD) / (coord->J[i] * coord->dy[i]); - } - return result + fromFieldAligned(yresult, "RGN_NOBNDRY"); -} +Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux); /*! * X-Z Finite Volume diffusion operator @@ -687,222 +155,7 @@ Field3D Div_Perp_Lap(const Field3D& a, const Field3D& f, CELL_LOC outloc = CELL_ template Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, Field3D& flow_ylow, - bool fixflux = true) { - - Coordinates* coord = f_in.getCoordinates(); - ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); - - if (f_in.isFci()) { - // Use mid-point (cell boundary) averages - - ASSERT1(f_in.hasParallelSlices()); - ASSERT1(v_in.hasParallelSlices()); - - const auto& f_up = f_in.yup(); - const auto& f_down = f_in.ydown(); - - const auto& v_up = v_in.yup(); - const auto& v_down = v_in.ydown(); - - Field3D result{emptyFrom(f_in)}; - BOUT_FOR(i, f_in.getRegion("RGN_NOBNDRY")) { - const auto iyp = i.yp(); - const auto iym = i.ym(); - - result[i] = (0.25 * (f_in[i] + f_up[iyp]) * (v_in[i] + v_up[iyp]) - * (coord->J[i] + coord->J.yup()[iyp]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22.yup()[iyp])) - - 0.25 * (f_in[i] + f_down[iym]) * (v_in[i] + v_down[iym]) - * (coord->J[i] + coord->J.ydown()[iym]) - / (sqrt(coord->g_22[i]) + sqrt(coord->g_22.ydown()[iym]))) - / (coord->dy[i] * coord->J[i]); - } - return result; - } - ASSERT1_FIELDS_COMPATIBLE(f_in, wave_speed_in); - - const Mesh* mesh = f_in.getMesh(); - - CellEdges cellboundary; - - ASSERT2(f_in.getDirectionY() == v_in.getDirectionY()); - ASSERT2(f_in.getDirectionY() == wave_speed_in.getDirectionY()); - const bool are_unaligned = - ((f_in.getDirectionY() == YDirectionType::Standard) - and (v_in.getDirectionY() == YDirectionType::Standard) - and (wave_speed_in.getDirectionY() == YDirectionType::Standard)); - - const Field3D f = are_unaligned ? toFieldAligned(f_in, "RGN_NOX") : f_in; - const Field3D v = are_unaligned ? toFieldAligned(v_in, "RGN_NOX") : v_in; - const Field3D wave_speed = - are_unaligned ? toFieldAligned(wave_speed_in, "RGN_NOX") : wave_speed_in; - - Field3D result{zeroFrom(f)}; - flow_ylow = zeroFrom(f); - - for (int i = mesh->xstart; i <= mesh->xend; i++) { - const bool is_periodic_y = mesh->periodicY(i); - const bool is_first_y = mesh->firstY(i); - const bool is_last_y = mesh->lastY(i); - - // Only need one guard cell, so no need to communicate fluxes Instead - // calculate in guard cells to get fluxes consistent between processors, but - // don't include the boundary cell. Note that this implies special handling - // of boundaries later - const int ys = (!is_first_y || is_periodic_y) ? mesh->ystart - 1 : mesh->ystart; - const int ye = (!is_last_y || is_periodic_y) ? mesh->yend + 1 : mesh->yend; - - for (int j = ys; j <= ye; j++) { - // Pre-calculate factors which multiply fluxes -#if not(BOUT_USE_METRIC_3D) - // For right cell boundaries - const BoutReal common_factor_r = - (coord->J(i, j) + coord->J(i, j + 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1))); - - const BoutReal flux_factor_rc = - common_factor_r / (coord->dy(i, j) * coord->J(i, j)); - const BoutReal flux_factor_rp = - common_factor_r / (coord->dy(i, j + 1) * coord->J(i, j + 1)); - - const BoutReal area_rp = - common_factor_r * coord->dx(i, j + 1) * coord->dz(i, j + 1); - - // For left cell boundaries - const BoutReal common_factor_l = - (coord->J(i, j) + coord->J(i, j - 1)) - / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1))); - - const BoutReal flux_factor_lc = - common_factor_l / (coord->dy(i, j) * coord->J(i, j)); - const BoutReal flux_factor_lm = - common_factor_l / (coord->dy(i, j - 1) * coord->J(i, j - 1)); - - const BoutReal area_lc = common_factor_l * coord->dx(i, j) * coord->dz(i, j); -#endif - for (int k = 0; k < mesh->LocalNz; k++) { -#if BOUT_USE_METRIC_3D - // For right cell boundaries - const BoutReal common_factor_r = - (coord->J(i, j, k) + coord->J(i, j + 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); - - const BoutReal flux_factor_rc = - common_factor_r / (coord->dy(i, j, k) * coord->J(i, j, k)); - const BoutReal flux_factor_rp = - common_factor_r / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); - - const BoutReal area_rp = - common_factor_r * coord->dx(i, j + 1, k) * coord->dz(i, j + 1, k); - - // For left cell boundaries - const BoutReal common_factor_l = - (coord->J(i, j, k) + coord->J(i, j - 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); - - const BoutReal flux_factor_lc = - common_factor_l / (coord->dy(i, j, k) * coord->J(i, j, k)); - const BoutReal flux_factor_lm = - common_factor_l / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); - - const BoutReal area_lc = - common_factor_l * coord->dx(i, j, k) * coord->dz(i, j, k); -#endif - - //////////////////////////////////////////// - // Reconstruct f at the cell faces - // This calculates s.R and s.L for the Right and Left - // face values on this cell - - // Reconstruct f at the cell faces - // TODO(peter): We can remove this #ifdef guard after switching to C++20 -#if __cpp_designated_initializers >= 201707L - Stencil1D s{.c = f(i, j, k), .m = f(i, j - 1, k), .p = f(i, j + 1, k)}; -#else - Stencil1D s{f(i, j, k), f(i, j - 1, k), f(i, j + 1, k), BoutNaN, - BoutNaN, BoutNaN, BoutNaN}; -#endif - cellboundary(s); // Calculate s.R and s.L - - //////////////////////////////////////////// - // Reconstruct v at the cell faces - // TODO(peter): We can remove this #ifdef guard after switching to C++20 -#if __cpp_designated_initializers >= 201707L - Stencil1D sv{.c = v(i, j, k), .m = v(i, j - 1, k), .p = v(i, j + 1, k)}; -#else - Stencil1D sv{v(i, j, k), v(i, j - 1, k), v(i, j + 1, k), BoutNaN, - BoutNaN, BoutNaN, BoutNaN}; -#endif - cellboundary(sv); // Calculate sv.R and sv.L - - //////////////////////////////////////////// - // Right boundary - - BoutReal flux = BoutNaN; - - if (is_last_y && (j == mesh->yend) && !is_periodic_y) { - // Last point in domain - - // Calculate velocity at right boundary (y+1/2) - const BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j + 1, k)); - - const BoutReal bndryval = 0.5 * (s.c + s.p); - if (fixflux) { - // Use mid-point to be consistent with boundary conditions - flux = bndryval * vpar; - } else { - // Add flux due to difference in boundary values - flux = (s.R * vpar) + (wave_speed(i, j, k) * (s.R - bndryval)); - } - - } else { - // Maximum wave speed in the two cells - const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k), - fabs(v(i, j, k)), fabs(v(i, j + 1, k))); - - flux = s.R * 0.5 * (sv.R + amax); - } - - result(i, j, k) += flux * flux_factor_rc; - result(i, j + 1, k) -= flux * flux_factor_rp; - - flow_ylow(i, j + 1, k) += flux * area_rp; - - //////////////////////////////////////////// - // Calculate at left boundary - - if (is_first_y && (j == mesh->ystart) && !is_periodic_y) { - // First point in domain - const BoutReal bndryval = 0.5 * (s.c + s.m); - const BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j - 1, k)); - if (fixflux) { - // Use mid-point to be consistent with boundary conditions - flux = bndryval * vpar; - } else { - // Add flux due to difference in boundary values - flux = (s.L * vpar) - (wave_speed(i, j, k) * (s.L - bndryval)); - } - } else { - - // Maximum wave speed in the two cells - const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k), - fabs(v(i, j, k)), fabs(v(i, j - 1, k))); - - flux = s.L * 0.5 * (sv.L - amax); - } - - result(i, j, k) -= flux * flux_factor_lc; - result(i, j - 1, k) += flux * flux_factor_lm; - - flow_ylow(i, j, k) += flux * area_lc; - } - } - } - if (are_unaligned) { - flow_ylow = fromFieldAligned(flow_ylow, "RGN_NOBNDRY"); - } - return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result; -} + bool fixflux = true); /// This operator calculates Div_par(f v v) /// It is used primarily (only?) in the parallel momentum equation. @@ -912,204 +165,8 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, /// fv is not interpolated to cell boundaries. template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, - const Field3D& wave_speed_in, bool fixflux = true) { - ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); - const Mesh* mesh = f_in.getMesh(); - const Coordinates* coord = f_in.getCoordinates(); - CellEdges cellboundary; - - if (f_in.isFci()) { - // FCI version, using yup/down fields - ASSERT1(f_in.hasParallelSlices()); - ASSERT1(v_in.hasParallelSlices()); - - const auto& B = coord->Bxy; - const auto& B_up = coord->Bxy.yup(); - const auto& B_down = coord->Bxy.ydown(); - - const auto& f_up = f_in.yup(); - const auto& f_down = f_in.ydown(); - - const auto& v_up = v_in.yup(); - const auto& v_down = v_in.ydown(); - - const auto& g_22 = coord->g_22; - const auto& dy = coord->dy; - - Field3D result{emptyFrom(f_in)}; - BOUT_FOR(i, f_in.getRegion("RGN_NOBNDRY")) { - const auto iyp = i.yp(); - const auto iym = i.ym(); - - // Maximum local wave speed - const BoutReal amax = - BOUTMAX(wave_speed_in[i], fabs(v_in[i]), fabs(v_up[iyp]), fabs(v_down[iym])); - - const BoutReal term = (f_up[iyp] * v_up[iyp] * v_up[iyp] / B_up[iyp]) - - (f_down[iym] * v_down[iym] * v_down[iym] / B_down[iym]); - - // Penalty terms. This implementation is very dissipative. - BoutReal penalty = - (amax * (f_in[i] * v_in[i] - f_up[iyp] * v_up[iyp]) / (B[i] + B_up[iyp])) - + (amax * (f_in[i] * v_in[i] - f_down[iym] * v_down[iym]) - / (B[i] + B_down[iym])); - - if (fabs(penalty) > fabs(term) and penalty * v_in[i] > 0) { - if (term * penalty > 0) { - penalty = term; - } else { - penalty = -term; - } - } - - result[i] = B[i] * (term + penalty) / (2 * dy[i] * sqrt(g_22[i])); - -#if CHECK > 0 - if (!std::isfinite(result[i])) { - throw BoutException("Non-finite value in Div_par_fvv at {}\n" - "fup {} vup {} fdown {} vdown {} amax {}\n", - "B {} Bup {} Bdown {} dy {} sqrt(g_22} {}", i, f_up[i], - v_up[i], f_down[i], v_down[i], amax, B[i], B_up[i], B_down[i], - dy[i], sqrt(g_22[i])); - } -#endif - } - return result; - } - - ASSERT1(areFieldsCompatible(f_in, wave_speed_in)); - - /// Ensure that f, v and wave_speed are field aligned - Field3D f = toFieldAligned(f_in, "RGN_NOX"); - Field3D v = toFieldAligned(v_in, "RGN_NOX"); - Field3D wave_speed = toFieldAligned(wave_speed_in, "RGN_NOX"); - - Field3D result{zeroFrom(f)}; - - for (int i = mesh->xstart; i <= mesh->xend; i++) { - const bool is_periodic_y = mesh->periodicY(i); - const bool is_first_y = mesh->firstY(i); - const bool is_last_y = mesh->lastY(i); - - // Only need one guard cell, so no need to communicate fluxes Instead - // calculate in guard cells to get fluxes consistent between processors, but - // don't include the boundary cell. Note that this implies special handling - // of boundaries later - const int ys = (!is_first_y || is_periodic_y) ? mesh->ystart - 1 : mesh->ystart; - const int ye = (!is_last_y || is_periodic_y) ? mesh->yend + 1 : mesh->yend; - - for (int j = ys; j <= ye; j++) { - // Pre-calculate factors which multiply fluxes - - for (int k = 0; k < mesh->LocalNz; k++) { - // For right cell boundaries - const BoutReal common_factor_r = - (coord->J(i, j, k) + coord->J(i, j + 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); - - const BoutReal flux_factor_rc = - common_factor_r / (coord->dy(i, j, k) * coord->J(i, j, k)); - const BoutReal flux_factor_rp = - common_factor_r / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); - - // For left cell boundaries - const BoutReal common_factor_l = - (coord->J(i, j, k) + coord->J(i, j - 1, k)) - / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); - - const BoutReal flux_factor_lc = - common_factor_l / (coord->dy(i, j, k) * coord->J(i, j, k)); - const BoutReal flux_factor_lm = - common_factor_l / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); - - //////////////////////////////////////////// - // Reconstruct f at the cell faces - // This calculates s.R and s.L for the Right and Left - // face values on this cell - - // Reconstruct f at the cell faces -#if __cpp_designated_initializers >= 201707L - Stencil1D s{.c = f(i, j, k), .m = f(i, j - 1, k), .p = f(i, j + 1, k)}; -#else - Stencil1D s{f(i, j, k), f(i, j - 1, k), f(i, j + 1, k), BoutNaN, - BoutNaN, BoutNaN, BoutNaN}; -#endif - cellboundary(s); // Calculate s.R and s.L - - //////////////////////////////////////////// - // Reconstruct v at the cell faces - // TODO(peter): We can remove this #ifdef guard after switching to C++20 -#if __cpp_designated_initializers >= 201707L - Stencil1D sv{.c = v(i, j, k), .m = v(i, j - 1, k), .p = v(i, j + 1, k)}; -#else - Stencil1D sv{v(i, j, k), v(i, j - 1, k), v(i, j + 1, k), BoutNaN, - BoutNaN, BoutNaN, BoutNaN}; -#endif - cellboundary(sv); - - //////////////////////////////////////////// - // Right boundary - - // Calculate velocity at right boundary (y+1/2) - const BoutReal v_mid_r = 0.5 * (sv.c + sv.p); - // And mid-point density at right boundary - const BoutReal n_mid_r = 0.5 * (s.c + s.p); - BoutReal flux = NAN; - - if (mesh->lastY(i) && (j == mesh->yend) && !mesh->periodicY(i)) { - // Last point in domain - - if (fixflux) { - // Use mid-point to be consistent with boundary conditions - flux = n_mid_r * v_mid_r * v_mid_r; - } else { - // Add flux due to difference in boundary values - flux = (s.R * sv.R * sv.R) // Use right cell edge values - + (BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.p)) * n_mid_r - * (sv.R - v_mid_r)); // Damp differences in velocity, not flux - } - } else { - // Maximum wave speed in the two cells - const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k), - fabs(sv.c), fabs(sv.p)); - - flux = s.R * 0.5 * (sv.R + amax) * sv.R; - } - - result(i, j, k) += flux * flux_factor_rc; - result(i, j + 1, k) -= flux * flux_factor_rp; - - //////////////////////////////////////////// - // Calculate at left boundary - - const BoutReal v_mid_l = 0.5 * (sv.c + sv.m); - const BoutReal n_mid_l = 0.5 * (s.c + s.m); - - if (mesh->firstY(i) && (j == mesh->ystart) && !mesh->periodicY(i)) { - // First point in domain - if (fixflux) { - // Use mid-point to be consistent with boundary conditions - flux = n_mid_l * v_mid_l * v_mid_l; - } else { - // Add flux due to difference in boundary values - flux = (s.L * sv.L * sv.L) - - (BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.m)) * n_mid_l - * (sv.L - v_mid_l)); - } - } else { - // Maximum wave speed in the two cells - const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k), - fabs(sv.c), fabs(sv.m)); - - flux = s.L * 0.5 * (sv.L - amax) * sv.L; - } - - result(i, j, k) -= flux * flux_factor_lc; - result(i, j - 1, k) += flux * flux_factor_lm; - } - } - } - return fromFieldAligned(result, "RGN_NOBNDRY"); -} + const Field3D& wave_speed_in, bool fixflux = true); +// extern template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, +// const Field3D& wave_speed_in, bool fixflux = true); } // namespace FV #endif // BOUT_FV_OPS_H diff --git a/src/mesh/fv_ops.cxx b/src/mesh/fv_ops.cxx index 6b8d8a6f21..a1333b1406 100644 --- a/src/mesh/fv_ops.cxx +++ b/src/mesh/fv_ops.cxx @@ -1,4 +1,5 @@ #include "bout/fv_ops.hxx" +#include "bout/fv_ops_impl.hxx" #include "bout/assert.hxx" #include "bout/bout_types.hxx" @@ -563,4 +564,72 @@ Field3D Div_Perp_Lap(const Field3D& a, const Field3D& f, CELL_LOC outloc) { return result; } +// BOUT_ENUM_CLASS(FluxLimiter, Upwind, Fromm, MinMod, MC, Superbee, VanAlbada, WENO3); + +template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_par(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux); +template Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, Field3D& flow_ylow, + bool fixflux = true); + +template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_par(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux); +template Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, Field3D& flow_ylow, + bool fixflux = true); + +template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_par(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux); +template Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, Field3D& flow_ylow, + bool fixflux = true); + +template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_par(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux); +template Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, Field3D& flow_ylow, + bool fixflux = true); + +template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_par(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, + bool bndry_flux); +template Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, Field3D& flow_ylow, + bool fixflux = true); + +template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, + bool fixflux = true); +template Field3D Div_par(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, + bool bndry_flux); +template Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, Field3D& flow_ylow, + bool fixflux = true); + +template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_par(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux = true); +template Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux); +template Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, Field3D& flow_ylow, + bool fixflux = true); + } // Namespace FV diff --git a/tests/unit/include/bout/test_fv_ops.cxx b/tests/unit/include/bout/test_fv_ops.cxx index e78c0cca41..5064b076b9 100644 --- a/tests/unit/include/bout/test_fv_ops.cxx +++ b/tests/unit/include/bout/test_fv_ops.cxx @@ -1,6 +1,6 @@ #include "gtest/gtest.h" -#include +#include TEST(FVOpsLimiterTest, VanAlbadaConstant) { FV::Stencil1D s{}; From 7a1a902ce46b4d119e673c9cd41e247c07859cf4 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 23 Jun 2026 14:30:42 +0200 Subject: [PATCH 19/38] Add to cmake file --- CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 0c42fbdeb2..239274c453 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -140,6 +140,7 @@ set(BOUT_SOURCES ./include/bout/fieldgroup.hxx ./include/bout/fieldperp.hxx ./include/bout/fv_ops.hxx + ./include/bout/fv_ops_impl.hxx ./include/bout/generic_factory.hxx ./include/bout/globalfield.hxx ./include/bout/globalindexer.hxx From 392f84217ec9dceb75ff4f72319fe19d09f3f075 Mon Sep 17 00:00:00 2001 From: David Bold Date: Tue, 23 Jun 2026 14:37:31 +0200 Subject: [PATCH 20/38] Add implementation file --- include/bout/fv_ops_impl.hxx | 1054 ++++++++++++++++++++++++++++++++++ 1 file changed, 1054 insertions(+) create mode 100644 include/bout/fv_ops_impl.hxx diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx new file mode 100644 index 0000000000..8326f73d28 --- /dev/null +++ b/include/bout/fv_ops_impl.hxx @@ -0,0 +1,1054 @@ +/* + Finite-volume discretisation methods. Flux-conservative form + */ + +#ifndef BOUT_FV_OPS_IMPL_H +#define BOUT_FV_OPS_IMPL_H + +#include "bout/assert.hxx" +#include "bout/bout_types.hxx" +#include "bout/boutexception.hxx" +#include "bout/build_defines.hxx" +#include "bout/coordinates.hxx" +#include "bout/field.hxx" +#include "bout/field3d.hxx" +#include "bout/globals.hxx" +#include "bout/mesh.hxx" +#include "bout/output_bout_types.hxx" // NOLINT(unused-includes, misc-include-cleaner) +#include "bout/region.hxx" +#include "bout/utils.hxx" +#include "bout/vector2d.hxx" + +#include + +namespace FV { +/*! + * Stencil used for Finite Volume calculations + * which includes cell face values L and R + */ +struct Stencil1D { + /// Cell centre values + BoutReal c; + BoutReal m; + BoutReal p; + BoutReal mm = BoutNaN; + BoutReal pp = BoutNaN; + + /// Left cell face value + BoutReal L = BoutNaN; + /// Right cell face value + BoutReal R = BoutNaN; +}; + +/*! + * First order upwind for testing + */ +struct Upwind { + void operator()(Stencil1D& n) { n.L = n.R = n.c; } +}; + +/*! + * Fromm method + */ +struct Fromm { + void operator()(Stencil1D& n) { + n.L = n.c - (0.25 * (n.p - n.m)); + n.R = n.c + (0.25 * (n.p - n.m)); + } +}; + +/*! + * Second order slope limiter method + * + * Limits slope to minimum absolute value + * of left and right gradients. If at a maximum + * or minimum slope set to zero, i.e. reverts + * to first order upwinding + */ +struct MinMod { + void operator()(Stencil1D& n) { + // Choose the gradient within the cell + // as the minimum (smoothest) solution + const BoutReal slope = _minmod(n.p - n.c, n.c - n.m); + n.L = n.c - (0.5 * slope); + n.R = n.c + (0.5 * slope); + } + +private: + /*! + * Internal helper function for minmod slope limiter + * + * If the inputs have different signs then + * returns zero, otherwise chooses the value + * with the minimum magnitude. + */ + static BoutReal _minmod(BoutReal a, BoutReal b) { + if (a * b <= 0.0) { + return 0.0; + } + + if (fabs(a) < fabs(b)) { + return a; + } + return b; + } +}; + +/*! + * Monotonised Central (MC) second order slope limiter (Van Leer) + * + * Limits the slope based on taking the slope with + * the minimum absolute value from central, 2*left and + * 2*right. If any of these slopes have different signs + * then the slope reverts to zero (i.e. 1st-order upwinding). + */ +struct MC { + void operator()(Stencil1D& n) { + const BoutReal slope = minmod(2. * (n.p - n.c), // 2*right difference + 0.5 * (n.p - n.m), // Central difference + 2. * (n.c - n.m)); // 2*left difference + n.L = n.c - (0.5 * slope); + n.R = n.c + (0.5 * slope); + } + +private: + // Return zero if any signs are different + // otherwise return the value with the minimum magnitude + static BoutReal minmod(BoutReal a, BoutReal b, BoutReal c) { + // if any of the signs are different, return zero gradient + if ((a * b <= 0.0) || (a * c <= 0.0)) { + return 0.0; + } + + // Return the minimum absolute value + return SIGN(a) * BOUTMIN(fabs(a), fabs(b), fabs(c)); + } +}; + +/// Superbee limiter +/// +/// This corresponds to the limiter function +/// φ(r) = max(0, min(2r, 1), min(r,2) +/// +/// The value at cell right (i.e. i + 1/2) is: +/// +/// n.R = n.c - φ(r) (n.c - (n.p + n.c)/2) +/// = n.c + φ(r) (n.p - n.c)/2 +/// +/// Four regimes: +/// a) r < 1/2 -> φ(r) = 2r +/// n.R = n.c + gL +/// b) 1/2 < r < 1 -> φ(r) = 1 +/// n.R = n.c + gR/2 +/// c) 1 < r < 2 -> φ(r) = r +/// n.R = n.c + gL/2 +/// d) 2 < r -> φ(r) = 2 +/// n.R = n.c + gR +/// +/// where the left and right gradients are: +/// gL = n.c - n.m +/// gR = n.p - n.c +/// +struct Superbee { + void operator()(Stencil1D& n) { + const BoutReal gL = n.c - n.m; + const BoutReal gR = n.p - n.c; + + // r = gL / gR + // Limiter is φ(r) + if (gL * gR < 0) { + // Different signs => Zero gradient + n.L = n.R = n.c; + } else { + const BoutReal sign = SIGN(gL); + const BoutReal abs_gL = fabs(gL); + const BoutReal abs_gR = fabs(gR); + const BoutReal half_slope = + sign * BOUTMAX(BOUTMIN(abs_gL, 0.5 * abs_gR), BOUTMIN(abs_gR, 0.5 * abs_gL)); + n.L = n.c - half_slope; + n.R = n.c + half_slope; + } + } +}; + +/*! + * Symmetric Van Albada second order slope limiter + * + * Uses a smooth (differentiable) approximation to `max(a*b, 0)` to avoid + * introducing a kink at extrema, which can be helpful for nonlinear solvers + * and finite-difference Jacobian calculations. + * + * The limited slope is calculated from the left and right differences + * `dl = c - m` and `dr = p - c` as + * + * slope = (pos(dl*dr) * (dl + dr)) / (dl^2 + dr^2) + * + * where `pos(x)` is a smooth approximation to `max(x, 0)`. + */ +struct VanAlbada { + void operator()(Stencil1D& n) { + const BoutReal dl = n.c - n.m; + const BoutReal dr = n.p - n.c; + + const BoutReal denom = dl * dl + dr * dr; + + // Smoothness parameters: + // - keep division well-defined when dl=dr=0 + // - provide a differentiable approximation to max(dl*dr, 0) + const BoutReal eps = 1e-12 * denom + 1e-30; + + const BoutReal ab = dl * dr; + const BoutReal ab_pos = 0.5 * (ab + sqrt(ab * ab + eps * eps)); + + const BoutReal slope = (ab_pos * (dl + dr)) / (denom + eps); + + n.L = n.c - 0.5 * slope; + n.R = n.c + 0.5 * slope; + } +}; + +/*! + * WENO3-JS (Jiang-Shu) reconstruction to cell faces + * + * This is a third-order essentially non-oscillatory reconstruction using two + * candidate second-order polynomials and smoothness-weighted blending. + * + * Unlike TVD slope limiters (e.g. ``MC``), WENO reconstruction is generally + * smooth (differentiable) for all inputs, but it does not enforce strict + * monotonicity. + * + * Uses only the three-point stencil (`m`, `c`, `p`), so it is a drop-in + * replacement anywhere `Stencil1D` is populated with those values. + */ +struct WENO3 { + void operator()(Stencil1D& n) { + // Right face (between c and p): value from cell c (left state at i+1/2) + const BoutReal p0_r = 0.5 * (-n.m + 3.0 * n.c); + const BoutReal p1_r = 0.5 * (n.c + n.p); + + const BoutReal beta0_r = SQ(n.c - n.m); + const BoutReal beta1_r = SQ(n.p - n.c); + + // Left face (between m and c): value from cell c (right state at i-1/2) + const BoutReal p0_l = 0.5 * (-n.p + 3.0 * n.c); + const BoutReal p1_l = 0.5 * (n.m + n.c); + + const BoutReal beta0_l = beta1_r; + const BoutReal beta1_l = beta0_r; + + // Smoothness parameter (scaled to local variation) + const BoutReal eps = 1e-12 * (beta0_r + beta1_r) + 1e-30; + + // Linear weights for WENO3-JS + constexpr BoutReal d0 = 1.0 / 3.0; + constexpr BoutReal d1 = 2.0 / 3.0; + + // Right face weights + const BoutReal a0_r = d0 / SQ(eps + beta0_r); + const BoutReal a1_r = d1 / SQ(eps + beta1_r); + const BoutReal wsum_r = a0_r + a1_r; + const BoutReal w0_r = a0_r / wsum_r; + const BoutReal w1_r = a1_r / wsum_r; + + // Left face weights (mirrored) + const BoutReal a0_l = d0 / SQ(eps + beta0_l); + const BoutReal a1_l = d1 / SQ(eps + beta1_l); + const BoutReal wsum_l = a0_l + a1_l; + const BoutReal w0_l = a0_l / wsum_l; + const BoutReal w1_l = a1_l / wsum_l; + + n.R = w0_r * p0_r + w1_r * p1_r; + n.L = w0_l * p0_l + w1_l * p1_l; + } +}; + +/*! + * Communicate fluxes between processors + * Takes values in guard cells, and adds them to cells + */ +void communicateFluxes(Field3D& f); + +/// Finite volume parallel divergence +/// +/// Preserves the sum of f*J*dx*dy*dz over the domain +/// +/// @param[in] f_in The field being advected. +/// This will be reconstructed at cell faces +/// using the given CellEdges method +/// @param[in] v_in The advection velocity. +/// This will be interpolated to cell boundaries +/// using linear interpolation +/// @param[in] wave_speed_in Local maximum speed of all waves in the system at each +// point in space +/// @param[in] fixflux Fix the flux at the boundary to be the value at the +/// midpoint (for boundary conditions) +/// +/// NB: Uses to/from FieldAligned coordinates +template +Field3D Div_par(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, + bool fixflux) { + + ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); + ASSERT1_FIELDS_COMPATIBLE(f_in, wave_speed_in); + + const Mesh* mesh = f_in.getMesh(); + + CellEdges cellboundary; + + ASSERT2(f_in.getDirectionY() == v_in.getDirectionY()); + ASSERT2(f_in.getDirectionY() == wave_speed_in.getDirectionY()); + const bool are_unaligned = + ((f_in.getDirectionY() == YDirectionType::Standard) + and (v_in.getDirectionY() == YDirectionType::Standard) + and (wave_speed_in.getDirectionY() == YDirectionType::Standard)); + + Field3D f = are_unaligned ? toFieldAligned(f_in, "RGN_NOX") : f_in; + Field3D v = are_unaligned ? toFieldAligned(v_in, "RGN_NOX") : v_in; + Field3D wave_speed = + are_unaligned ? toFieldAligned(wave_speed_in, "RGN_NOX") : wave_speed_in; + + Coordinates* coord = f_in.getCoordinates(); + + Field3D result{zeroFrom(f)}; + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + const bool is_periodic_y = mesh->periodicY(i); + const bool is_first_y = mesh->firstY(i); + const bool is_last_y = mesh->lastY(i); + + // Only need one guard cell, so no need to communicate fluxes Instead + // calculate in guard cells to get fluxes consistent between processors, but + // don't include the boundary cell. Note that this implies special handling + // of boundaries later + const int ys = (!is_first_y || is_periodic_y) ? mesh->ystart - 1 : mesh->ystart; + const int ye = (!is_last_y || is_periodic_y) ? mesh->yend + 1 : mesh->yend; + + for (int j = ys; j <= ye; j++) { + // Pre-calculate factors which multiply fluxes +#if not(BOUT_USE_METRIC_3D) + // For right cell boundaries + BoutReal common_factor = (coord->J(i, j) + coord->J(i, j + 1)) + / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1))); + + const BoutReal flux_factor_rc = common_factor / (coord->dy(i, j) * coord->J(i, j)); + const BoutReal flux_factor_rp = + common_factor / (coord->dy(i, j + 1) * coord->J(i, j + 1)); + + // For left cell boundaries + common_factor = (coord->J(i, j) + coord->J(i, j - 1)) + / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1))); + + const BoutReal flux_factor_lc = common_factor / (coord->dy(i, j) * coord->J(i, j)); + const BoutReal flux_factor_lm = + common_factor / (coord->dy(i, j - 1) * coord->J(i, j - 1)); +#endif + for (int k = mesh->zstart; k <= mesh->zend; k++) { +#if BOUT_USE_METRIC_3D + // For right cell boundaries + BoutReal common_factor = + (coord->J(i, j, k) + coord->J(i, j + 1, k)) + / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); + + BoutReal flux_factor_rc = + common_factor / (coord->dy(i, j, k) * coord->J(i, j, k)); + BoutReal flux_factor_rp = + common_factor / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); + + // For left cell boundaries + common_factor = (coord->J(i, j, k) + coord->J(i, j - 1, k)) + / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); + + BoutReal flux_factor_lc = + common_factor / (coord->dy(i, j, k) * coord->J(i, j, k)); + BoutReal flux_factor_lm = + common_factor / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); +#endif + + //////////////////////////////////////////// + // Reconstruct f at the cell faces + // This calculates s.R and s.L for the Right and Left + // face values on this cell + + // Reconstruct f at the cell faces + Stencil1D s; + s.c = f(i, j, k); + s.m = f(i, j - 1, k); + s.p = f(i, j + 1, k); + + cellboundary(s); // Calculate s.R and s.L + + //////////////////////////////////////////// + // Right boundary + + // Calculate velocity at right boundary (y+1/2) + BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j + 1, k)); + BoutReal flux = NAN; + + if (is_last_y && (j == mesh->yend) && !is_periodic_y) { + // Last point in domain + + const BoutReal bndryval = 0.5 * (s.c + s.p); + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = bndryval * vpar; + } else { + // Add flux due to difference in boundary values + flux = (s.R * vpar) + (wave_speed(i, j, k) * (s.R - bndryval)); + } + } else { + + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k)); + + if (vpar > amax) { + // Supersonic flow out of this cell + flux = s.R * vpar; + } else if (vpar < -amax) { + // Supersonic flow into this cell + flux = 0.0; + } else { + // Subsonic flow, so a mix of right and left fluxes + flux = s.R * 0.5 * (vpar + amax); + } + } + + result(i, j, k) += flux * flux_factor_rc; + result(i, j + 1, k) -= flux * flux_factor_rp; + + //////////////////////////////////////////// + // Calculate at left boundary + + vpar = 0.5 * (v(i, j, k) + v(i, j - 1, k)); + + if (is_first_y && (j == mesh->ystart) && !is_periodic_y) { + // First point in domain + const BoutReal bndryval = 0.5 * (s.c + s.m); + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = bndryval * vpar; + } else { + // Add flux due to difference in boundary values + flux = (s.L * vpar) - (wave_speed(i, j, k) * (s.L - bndryval)); + } + } else { + + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k)); + + if (vpar < -amax) { + // Supersonic out of this cell + flux = s.L * vpar; + } else if (vpar > amax) { + // Supersonic into this cell + flux = 0.0; + } else { + flux = s.L * 0.5 * (vpar - amax); + } + } + + result(i, j, k) -= flux * flux_factor_lc; + result(i, j - 1, k) += flux * flux_factor_lm; + } + } + } + return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result; +} + +/*! + * Div ( n * v ) -- Magnetic drifts + * + * This uses the expression + * + * Div( A ) = 1/J * d/di ( J * A^i ) + * + * Hence the input vector should be contravariant + * + * Note: Uses to/from FieldAligned + * + */ +template +Field3D Div_f_v(const Field3D& n_in, const Vector3D& v, bool bndry_flux) { + ASSERT1(n_in.getLocation() == v.getLocation()); + ASSERT1_FIELDS_COMPATIBLE(n_in, v.x); + + const Mesh* mesh = n_in.getMesh(); + + CellEdges cellboundary; + + Coordinates* coord = n_in.getCoordinates(); + + if (v.covariant) { + // Got a covariant vector instead + throw BoutException("Div_f_v passed a covariant v"); + } + + Field3D result{zeroFrom(n_in)}; + + Field3D vx = v.x; + Field3D vz = v.z; + Field3D n = n_in; + + BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { + // Calculate velocities + const BoutReal vU = 0.25 * (vz[i.zp()] + vz[i]) * (coord->J[i.zp()] + coord->J[i]); + const BoutReal vD = 0.25 * (vz[i.zm()] + vz[i]) * (coord->J[i.zm()] + coord->J[i]); + const BoutReal vL = 0.25 * (vx[i.xm()] + vx[i]) * (coord->J[i.xm()] + coord->J[i]); + const BoutReal vR = 0.25 * (vx[i.xp()] + vx[i]) * (coord->J[i.xp()] + coord->J[i]); + + // X direction + Stencil1D s; + s.c = n[i]; + s.m = n[i.xm()]; + s.mm = n[i.xmm()]; + s.p = n[i.xp()]; + s.pp = n[i.xpp()]; + + cellboundary(s); + + if ((i.x() == mesh->xend) && (mesh->lastX())) { + // At right boundary in X + if (bndry_flux) { + BoutReal flux = NAN; + if (vR > 0.0) { + // Flux to boundary + flux = vR * s.R; + } else { + // Flux in from boundary + flux = vR * 0.5 * (n[i.xp()] + n[i]); + } + result[i] += flux / (coord->dx[i] * coord->J[i]); + result[i.xp()] -= flux / (coord->dx[i.xp()] * coord->J[i.xp()]); + } + } else { + // Not at a boundary + if (vR > 0.0) { + // Flux out into next cell + const BoutReal flux = vR * s.R; + result[i] += flux / (coord->dx[i] * coord->J[i]); + result[i.xp()] -= flux / (coord->dx[i.xp()] * coord->J[i.xp()]); + } + } + + // Left side + + if ((i.x() == mesh->xstart) && (mesh->firstX())) { + // At left boundary in X + + if (bndry_flux) { + BoutReal flux = NAN; + if (vL < 0.0) { + // Flux to boundary + flux = vL * s.L; + } else { + // Flux in from boundary + flux = vL * 0.5 * (n[i.xm()] + n[i]); + } + result[i] -= flux / (coord->dx[i] * coord->J[i]); + result[i.xm()] += flux / (coord->dx[i.xm()] * coord->J[i.xm()]); + } + } else { + // Not at a boundary + if (vL < 0.0) { + const BoutReal flux = vL * s.L; + result[i] -= flux / (coord->dx[i] * coord->J[i]); + result[i.xm()] += flux / (coord->dx[i.xm()] * coord->J[i.xm()]); + } + } + + /// NOTE: Need to communicate fluxes + + // Z direction + s.m = n[i.zm()]; + s.mm = n[i.zmm()]; + s.p = n[i.zp()]; + s.pp = n[i.zpp()]; + + cellboundary(s); + + if (vU > 0.0) { + const BoutReal flux = vU * s.R; + result[i] += flux / (coord->J[i] * coord->dz[i]); + result[i.zp()] -= flux / (coord->J[i.zp()] * coord->dz[i.zp()]); + } + if (vD < 0.0) { + const BoutReal flux = vD * s.L; + result[i] -= flux / (coord->J[i] * coord->dz[i]); + result[i.zm()] += flux / (coord->J[i.zm()] * coord->dz[i.zm()]); + } + } + + communicateFluxes(result); + + // Y advection + // Currently just using simple centered differences + // so no fluxes need to be exchanged + + n = toFieldAligned(n_in, "RGN_NOX"); + Field3D vy = toFieldAligned(v.y, "RGN_NOX"); + + Field3D yresult = 0.0; + yresult.setDirectionY(YDirectionType::Aligned); + + BOUT_FOR(i, result.getRegion("RGN_NOBNDRY")) { + // Y velocities on y boundaries + const BoutReal vU = 0.25 * (vy[i] + vy[i.yp()]) * (coord->J[i] + coord->J[i.yp()]); + const BoutReal vD = 0.25 * (vy[i] + vy[i.ym()]) * (coord->J[i] + coord->J[i.ym()]); + + // n (advected quantity) on y boundaries + // Note: Use unshifted n_in variable + const BoutReal nU = 0.5 * (n[i] + n[i.yp()]); + const BoutReal nD = 0.5 * (n[i] + n[i.ym()]); + + yresult[i] = (nU * vU - nD * vD) / (coord->J[i] * coord->dy[i]); + } + return result + fromFieldAligned(yresult, "RGN_NOBNDRY"); +} + +/// Finite volume parallel divergence +/// +/// NOTE: Modified version, applies limiter to velocity and field +/// Performs better (smaller overshoots) than Div_par +/// +/// Preserves the sum of f*J*dx*dy*dz over the domain +/// +/// @param[in] f_in The field being advected. +/// This will be reconstructed at cell faces +/// using the given CellEdges method +/// @param[in] v_in The advection velocity. +/// This will be interpolated to cell boundaries +/// using linear interpolation +/// @param[in] wave_speed_in Local maximum speed of all waves in the system at each +// point in space +/// @param[in] fixflux Fix the flux at the boundary to be the value at the +/// midpoint (for boundary conditions) +/// +/// @param[out] flow_ylow Flow at the lower Y cell boundary +/// Already includes area factor * flux +template +Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, Field3D& flow_ylow, bool fixflux) { + + Coordinates* coord = f_in.getCoordinates(); + ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); + + if (f_in.isFci()) { + // Use mid-point (cell boundary) averages + + ASSERT1(f_in.hasParallelSlices()); + ASSERT1(v_in.hasParallelSlices()); + + const auto& f_up = f_in.yup(); + const auto& f_down = f_in.ydown(); + + const auto& v_up = v_in.yup(); + const auto& v_down = v_in.ydown(); + + Field3D result{emptyFrom(f_in)}; + BOUT_FOR(i, f_in.getRegion("RGN_NOBNDRY")) { + const auto iyp = i.yp(); + const auto iym = i.ym(); + + result[i] = (0.25 * (f_in[i] + f_up[iyp]) * (v_in[i] + v_up[iyp]) + * (coord->J[i] + coord->J.yup()[iyp]) + / (sqrt(coord->g_22[i]) + sqrt(coord->g_22.yup()[iyp])) + - 0.25 * (f_in[i] + f_down[iym]) * (v_in[i] + v_down[iym]) + * (coord->J[i] + coord->J.ydown()[iym]) + / (sqrt(coord->g_22[i]) + sqrt(coord->g_22.ydown()[iym]))) + / (coord->dy[i] * coord->J[i]); + } + return result; + } + ASSERT1_FIELDS_COMPATIBLE(f_in, wave_speed_in); + + const Mesh* mesh = f_in.getMesh(); + + CellEdges cellboundary; + + ASSERT2(f_in.getDirectionY() == v_in.getDirectionY()); + ASSERT2(f_in.getDirectionY() == wave_speed_in.getDirectionY()); + const bool are_unaligned = + ((f_in.getDirectionY() == YDirectionType::Standard) + and (v_in.getDirectionY() == YDirectionType::Standard) + and (wave_speed_in.getDirectionY() == YDirectionType::Standard)); + + const Field3D f = are_unaligned ? toFieldAligned(f_in, "RGN_NOX") : f_in; + const Field3D v = are_unaligned ? toFieldAligned(v_in, "RGN_NOX") : v_in; + const Field3D wave_speed = + are_unaligned ? toFieldAligned(wave_speed_in, "RGN_NOX") : wave_speed_in; + + Field3D result{zeroFrom(f)}; + flow_ylow = zeroFrom(f); + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + const bool is_periodic_y = mesh->periodicY(i); + const bool is_first_y = mesh->firstY(i); + const bool is_last_y = mesh->lastY(i); + + // Only need one guard cell, so no need to communicate fluxes Instead + // calculate in guard cells to get fluxes consistent between processors, but + // don't include the boundary cell. Note that this implies special handling + // of boundaries later + const int ys = (!is_first_y || is_periodic_y) ? mesh->ystart - 1 : mesh->ystart; + const int ye = (!is_last_y || is_periodic_y) ? mesh->yend + 1 : mesh->yend; + + for (int j = ys; j <= ye; j++) { + // Pre-calculate factors which multiply fluxes +#if not(BOUT_USE_METRIC_3D) + // For right cell boundaries + const BoutReal common_factor_r = + (coord->J(i, j) + coord->J(i, j + 1)) + / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j + 1))); + + const BoutReal flux_factor_rc = + common_factor_r / (coord->dy(i, j) * coord->J(i, j)); + const BoutReal flux_factor_rp = + common_factor_r / (coord->dy(i, j + 1) * coord->J(i, j + 1)); + + const BoutReal area_rp = + common_factor_r * coord->dx(i, j + 1) * coord->dz(i, j + 1); + + // For left cell boundaries + const BoutReal common_factor_l = + (coord->J(i, j) + coord->J(i, j - 1)) + / (sqrt(coord->g_22(i, j)) + sqrt(coord->g_22(i, j - 1))); + + const BoutReal flux_factor_lc = + common_factor_l / (coord->dy(i, j) * coord->J(i, j)); + const BoutReal flux_factor_lm = + common_factor_l / (coord->dy(i, j - 1) * coord->J(i, j - 1)); + + const BoutReal area_lc = common_factor_l * coord->dx(i, j) * coord->dz(i, j); +#endif + for (int k = 0; k < mesh->LocalNz; k++) { +#if BOUT_USE_METRIC_3D + // For right cell boundaries + const BoutReal common_factor_r = + (coord->J(i, j, k) + coord->J(i, j + 1, k)) + / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); + + const BoutReal flux_factor_rc = + common_factor_r / (coord->dy(i, j, k) * coord->J(i, j, k)); + const BoutReal flux_factor_rp = + common_factor_r / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); + + const BoutReal area_rp = + common_factor_r * coord->dx(i, j + 1, k) * coord->dz(i, j + 1, k); + + // For left cell boundaries + const BoutReal common_factor_l = + (coord->J(i, j, k) + coord->J(i, j - 1, k)) + / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); + + const BoutReal flux_factor_lc = + common_factor_l / (coord->dy(i, j, k) * coord->J(i, j, k)); + const BoutReal flux_factor_lm = + common_factor_l / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); + + const BoutReal area_lc = + common_factor_l * coord->dx(i, j, k) * coord->dz(i, j, k); +#endif + + //////////////////////////////////////////// + // Reconstruct f at the cell faces + // This calculates s.R and s.L for the Right and Left + // face values on this cell + + // Reconstruct f at the cell faces + // TODO(peter): We can remove this #ifdef guard after switching to C++20 +#if __cpp_designated_initializers >= 201707L + Stencil1D s{.c = f(i, j, k), .m = f(i, j - 1, k), .p = f(i, j + 1, k)}; +#else + Stencil1D s{f(i, j, k), f(i, j - 1, k), f(i, j + 1, k), BoutNaN, + BoutNaN, BoutNaN, BoutNaN}; +#endif + cellboundary(s); // Calculate s.R and s.L + + //////////////////////////////////////////// + // Reconstruct v at the cell faces + // TODO(peter): We can remove this #ifdef guard after switching to C++20 +#if __cpp_designated_initializers >= 201707L + Stencil1D sv{.c = v(i, j, k), .m = v(i, j - 1, k), .p = v(i, j + 1, k)}; +#else + Stencil1D sv{v(i, j, k), v(i, j - 1, k), v(i, j + 1, k), BoutNaN, + BoutNaN, BoutNaN, BoutNaN}; +#endif + cellboundary(sv); // Calculate sv.R and sv.L + + //////////////////////////////////////////// + // Right boundary + + BoutReal flux = BoutNaN; + + if (is_last_y && (j == mesh->yend) && !is_periodic_y) { + // Last point in domain + + // Calculate velocity at right boundary (y+1/2) + const BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j + 1, k)); + + const BoutReal bndryval = 0.5 * (s.c + s.p); + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = bndryval * vpar; + } else { + // Add flux due to difference in boundary values + flux = (s.R * vpar) + (wave_speed(i, j, k) * (s.R - bndryval)); + } + + } else { + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k), + fabs(v(i, j, k)), fabs(v(i, j + 1, k))); + + flux = s.R * 0.5 * (sv.R + amax); + } + + result(i, j, k) += flux * flux_factor_rc; + result(i, j + 1, k) -= flux * flux_factor_rp; + + flow_ylow(i, j + 1, k) += flux * area_rp; + + //////////////////////////////////////////// + // Calculate at left boundary + + if (is_first_y && (j == mesh->ystart) && !is_periodic_y) { + // First point in domain + const BoutReal bndryval = 0.5 * (s.c + s.m); + const BoutReal vpar = 0.5 * (v(i, j, k) + v(i, j - 1, k)); + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = bndryval * vpar; + } else { + // Add flux due to difference in boundary values + flux = (s.L * vpar) - (wave_speed(i, j, k) * (s.L - bndryval)); + } + } else { + + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k), + fabs(v(i, j, k)), fabs(v(i, j - 1, k))); + + flux = s.L * 0.5 * (sv.L - amax); + } + + result(i, j, k) -= flux * flux_factor_lc; + result(i, j - 1, k) += flux * flux_factor_lm; + + flow_ylow(i, j, k) += flux * area_lc; + } + } + } + if (are_unaligned) { + flow_ylow = fromFieldAligned(flow_ylow, "RGN_NOBNDRY"); + } + return are_unaligned ? fromFieldAligned(result, "RGN_NOBNDRY") : result; +} + +/// This operator calculates Div_par(f v v) +/// It is used primarily (only?) in the parallel momentum equation. +/// +/// This operator is used rather than Div(f fv) so that the values of +/// f and v are consistent with other advection equations: The product +/// fv is not interpolated to cell boundaries. +template +Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, + const Field3D& wave_speed_in, bool fixflux) { + ASSERT1_FIELDS_COMPATIBLE(f_in, v_in); + const Mesh* mesh = f_in.getMesh(); + const Coordinates* coord = f_in.getCoordinates(); + CellEdges cellboundary; + + if (f_in.isFci()) { + // FCI version, using yup/down fields + ASSERT1(f_in.hasParallelSlices()); + ASSERT1(v_in.hasParallelSlices()); + + const auto& B = coord->Bxy; + const auto& B_up = coord->Bxy.yup(); + const auto& B_down = coord->Bxy.ydown(); + + const auto& f_up = f_in.yup(); + const auto& f_down = f_in.ydown(); + + const auto& v_up = v_in.yup(); + const auto& v_down = v_in.ydown(); + + const auto& g_22 = coord->g_22; + const auto& dy = coord->dy; + + Field3D result{emptyFrom(f_in)}; + BOUT_FOR(i, f_in.getRegion("RGN_NOBNDRY")) { + const auto iyp = i.yp(); + const auto iym = i.ym(); + + // Maximum local wave speed + const BoutReal amax = + BOUTMAX(wave_speed_in[i], fabs(v_in[i]), fabs(v_up[iyp]), fabs(v_down[iym])); + + const BoutReal term = (f_up[iyp] * v_up[iyp] * v_up[iyp] / B_up[iyp]) + - (f_down[iym] * v_down[iym] * v_down[iym] / B_down[iym]); + + // Penalty terms. This implementation is very dissipative. + BoutReal penalty = + (amax * (f_in[i] * v_in[i] - f_up[iyp] * v_up[iyp]) / (B[i] + B_up[iyp])) + + (amax * (f_in[i] * v_in[i] - f_down[iym] * v_down[iym]) + / (B[i] + B_down[iym])); + + if (fabs(penalty) > fabs(term) and penalty * v_in[i] > 0) { + if (term * penalty > 0) { + penalty = term; + } else { + penalty = -term; + } + } + + result[i] = B[i] * (term + penalty) / (2 * dy[i] * sqrt(g_22[i])); + +#if CHECK > 0 + if (!std::isfinite(result[i])) { + throw BoutException("Non-finite value in Div_par_fvv at {}\n" + "fup {} vup {} fdown {} vdown {} amax {}\n", + "B {} Bup {} Bdown {} dy {} sqrt(g_22} {}", i, f_up[i], + v_up[i], f_down[i], v_down[i], amax, B[i], B_up[i], B_down[i], + dy[i], sqrt(g_22[i])); + } +#endif + } + return result; + } + + ASSERT1(areFieldsCompatible(f_in, wave_speed_in)); + + /// Ensure that f, v and wave_speed are field aligned + Field3D f = toFieldAligned(f_in, "RGN_NOX"); + Field3D v = toFieldAligned(v_in, "RGN_NOX"); + Field3D wave_speed = toFieldAligned(wave_speed_in, "RGN_NOX"); + + Field3D result{zeroFrom(f)}; + + for (int i = mesh->xstart; i <= mesh->xend; i++) { + const bool is_periodic_y = mesh->periodicY(i); + const bool is_first_y = mesh->firstY(i); + const bool is_last_y = mesh->lastY(i); + + // Only need one guard cell, so no need to communicate fluxes Instead + // calculate in guard cells to get fluxes consistent between processors, but + // don't include the boundary cell. Note that this implies special handling + // of boundaries later + const int ys = (!is_first_y || is_periodic_y) ? mesh->ystart - 1 : mesh->ystart; + const int ye = (!is_last_y || is_periodic_y) ? mesh->yend + 1 : mesh->yend; + + for (int j = ys; j <= ye; j++) { + // Pre-calculate factors which multiply fluxes + + for (int k = 0; k < mesh->LocalNz; k++) { + // For right cell boundaries + const BoutReal common_factor_r = + (coord->J(i, j, k) + coord->J(i, j + 1, k)) + / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j + 1, k))); + + const BoutReal flux_factor_rc = + common_factor_r / (coord->dy(i, j, k) * coord->J(i, j, k)); + const BoutReal flux_factor_rp = + common_factor_r / (coord->dy(i, j + 1, k) * coord->J(i, j + 1, k)); + + // For left cell boundaries + const BoutReal common_factor_l = + (coord->J(i, j, k) + coord->J(i, j - 1, k)) + / (sqrt(coord->g_22(i, j, k)) + sqrt(coord->g_22(i, j - 1, k))); + + const BoutReal flux_factor_lc = + common_factor_l / (coord->dy(i, j, k) * coord->J(i, j, k)); + const BoutReal flux_factor_lm = + common_factor_l / (coord->dy(i, j - 1, k) * coord->J(i, j - 1, k)); + + //////////////////////////////////////////// + // Reconstruct f at the cell faces + // This calculates s.R and s.L for the Right and Left + // face values on this cell + + // Reconstruct f at the cell faces +#if __cpp_designated_initializers >= 201707L + Stencil1D s{.c = f(i, j, k), .m = f(i, j - 1, k), .p = f(i, j + 1, k)}; +#else + Stencil1D s{f(i, j, k), f(i, j - 1, k), f(i, j + 1, k), BoutNaN, + BoutNaN, BoutNaN, BoutNaN}; +#endif + cellboundary(s); // Calculate s.R and s.L + + //////////////////////////////////////////// + // Reconstruct v at the cell faces + // TODO(peter): We can remove this #ifdef guard after switching to C++20 +#if __cpp_designated_initializers >= 201707L + Stencil1D sv{.c = v(i, j, k), .m = v(i, j - 1, k), .p = v(i, j + 1, k)}; +#else + Stencil1D sv{v(i, j, k), v(i, j - 1, k), v(i, j + 1, k), BoutNaN, + BoutNaN, BoutNaN, BoutNaN}; +#endif + cellboundary(sv); + + //////////////////////////////////////////// + // Right boundary + + // Calculate velocity at right boundary (y+1/2) + const BoutReal v_mid_r = 0.5 * (sv.c + sv.p); + // And mid-point density at right boundary + const BoutReal n_mid_r = 0.5 * (s.c + s.p); + BoutReal flux = NAN; + + if (mesh->lastY(i) && (j == mesh->yend) && !mesh->periodicY(i)) { + // Last point in domain + + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = n_mid_r * v_mid_r * v_mid_r; + } else { + // Add flux due to difference in boundary values + flux = (s.R * sv.R * sv.R) // Use right cell edge values + + (BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.p)) * n_mid_r + * (sv.R - v_mid_r)); // Damp differences in velocity, not flux + } + } else { + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k), + fabs(sv.c), fabs(sv.p)); + + flux = s.R * 0.5 * (sv.R + amax) * sv.R; + } + + result(i, j, k) += flux * flux_factor_rc; + result(i, j + 1, k) -= flux * flux_factor_rp; + + //////////////////////////////////////////// + // Calculate at left boundary + + const BoutReal v_mid_l = 0.5 * (sv.c + sv.m); + const BoutReal n_mid_l = 0.5 * (s.c + s.m); + + if (mesh->firstY(i) && (j == mesh->ystart) && !mesh->periodicY(i)) { + // First point in domain + if (fixflux) { + // Use mid-point to be consistent with boundary conditions + flux = n_mid_l * v_mid_l * v_mid_l; + } else { + // Add flux due to difference in boundary values + flux = (s.L * sv.L * sv.L) + - (BOUTMAX(wave_speed(i, j, k), fabs(sv.c), fabs(sv.m)) * n_mid_l + * (sv.L - v_mid_l)); + } + } else { + // Maximum wave speed in the two cells + const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k), + fabs(sv.c), fabs(sv.m)); + + flux = s.L * 0.5 * (sv.L - amax) * sv.L; + } + + result(i, j, k) -= flux * flux_factor_lc; + result(i, j - 1, k) += flux * flux_factor_lm; + } + } + } + return fromFieldAligned(result, "RGN_NOBNDRY"); +} +} // namespace FV +#endif // BOUT_FV_OPS_H From e5ba73d8a978308d9934a36efba4cfab2e2fde68 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 24 Jun 2026 11:32:51 +0200 Subject: [PATCH 21/38] Improve code comments Co-authored-by: Peter Hill --- include/bout/fv_ops.hxx | 4 +--- src/mesh/fv_ops.cxx | 3 +-- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index d256bad32a..89caaadf97 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -77,7 +77,7 @@ Field3D D4DY4(const Field3D& d, const Field3D& f); */ Field3D D4DY4_Index(const Field3D& f, bool bndry_flux = true); -// FluxLimiter +// Forward declarations of flux limiters class Upwind; class Fromm; class MinMod; @@ -166,7 +166,5 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, bool fixflux = true); -// extern template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, -// const Field3D& wave_speed_in, bool fixflux = true); } // namespace FV #endif // BOUT_FV_OPS_H diff --git a/src/mesh/fv_ops.cxx b/src/mesh/fv_ops.cxx index a1333b1406..1bd00d34a9 100644 --- a/src/mesh/fv_ops.cxx +++ b/src/mesh/fv_ops.cxx @@ -564,8 +564,7 @@ Field3D Div_Perp_Lap(const Field3D& a, const Field3D& f, CELL_LOC outloc) { return result; } -// BOUT_ENUM_CLASS(FluxLimiter, Upwind, Fromm, MinMod, MC, Superbee, VanAlbada, WENO3); - +// Explicit instantiations of flux-limited finite volume methods template Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, bool fixflux = true); template Field3D Div_par(const Field3D& f_in, const Field3D& v_in, From f8a246996d733c52e6183797b02123e20f14b5e4 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 24 Jun 2026 11:35:51 +0200 Subject: [PATCH 22/38] Remove pre-C++20 fall back path --- include/bout/fv_ops_impl.hxx | 23 ----------------------- 1 file changed, 23 deletions(-) diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx index 8326f73d28..e27fd32828 100644 --- a/include/bout/fv_ops_impl.hxx +++ b/include/bout/fv_ops_impl.hxx @@ -754,24 +754,12 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, // face values on this cell // Reconstruct f at the cell faces - // TODO(peter): We can remove this #ifdef guard after switching to C++20 -#if __cpp_designated_initializers >= 201707L Stencil1D s{.c = f(i, j, k), .m = f(i, j - 1, k), .p = f(i, j + 1, k)}; -#else - Stencil1D s{f(i, j, k), f(i, j - 1, k), f(i, j + 1, k), BoutNaN, - BoutNaN, BoutNaN, BoutNaN}; -#endif cellboundary(s); // Calculate s.R and s.L //////////////////////////////////////////// // Reconstruct v at the cell faces - // TODO(peter): We can remove this #ifdef guard after switching to C++20 -#if __cpp_designated_initializers >= 201707L Stencil1D sv{.c = v(i, j, k), .m = v(i, j - 1, k), .p = v(i, j + 1, k)}; -#else - Stencil1D sv{v(i, j, k), v(i, j - 1, k), v(i, j + 1, k), BoutNaN, - BoutNaN, BoutNaN, BoutNaN}; -#endif cellboundary(sv); // Calculate sv.R and sv.L //////////////////////////////////////////// @@ -967,23 +955,12 @@ Field3D Div_par_fvv(const Field3D& f_in, const Field3D& v_in, // face values on this cell // Reconstruct f at the cell faces -#if __cpp_designated_initializers >= 201707L Stencil1D s{.c = f(i, j, k), .m = f(i, j - 1, k), .p = f(i, j + 1, k)}; -#else - Stencil1D s{f(i, j, k), f(i, j - 1, k), f(i, j + 1, k), BoutNaN, - BoutNaN, BoutNaN, BoutNaN}; -#endif cellboundary(s); // Calculate s.R and s.L //////////////////////////////////////////// // Reconstruct v at the cell faces - // TODO(peter): We can remove this #ifdef guard after switching to C++20 -#if __cpp_designated_initializers >= 201707L Stencil1D sv{.c = v(i, j, k), .m = v(i, j - 1, k), .p = v(i, j + 1, k)}; -#else - Stencil1D sv{v(i, j, k), v(i, j - 1, k), v(i, j + 1, k), BoutNaN, - BoutNaN, BoutNaN, BoutNaN}; -#endif cellboundary(sv); //////////////////////////////////////////// From 8f8e0359bcd2285dae281a3051caad334ef470ed Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 24 Jun 2026 11:45:56 +0200 Subject: [PATCH 23/38] Apply clang-tidy fixes --- include/bout/fv_ops_impl.hxx | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx index e27fd32828..91117a6961 100644 --- a/include/bout/fv_ops_impl.hxx +++ b/include/bout/fv_ops_impl.hxx @@ -28,9 +28,9 @@ namespace FV { */ struct Stencil1D { /// Cell centre values - BoutReal c; - BoutReal m; - BoutReal p; + BoutReal c{}; + BoutReal m{}; + BoutReal p{}; BoutReal mm = BoutNaN; BoutReal pp = BoutNaN; @@ -190,20 +190,20 @@ struct VanAlbada { const BoutReal dl = n.c - n.m; const BoutReal dr = n.p - n.c; - const BoutReal denom = dl * dl + dr * dr; + const BoutReal denom = (dl * dl) + (dr * dr); // Smoothness parameters: // - keep division well-defined when dl=dr=0 // - provide a differentiable approximation to max(dl*dr, 0) - const BoutReal eps = 1e-12 * denom + 1e-30; + const BoutReal eps = (1e-12 * denom) + 1e-30; const BoutReal ab = dl * dr; - const BoutReal ab_pos = 0.5 * (ab + sqrt(ab * ab + eps * eps)); + const BoutReal ab_pos = 0.5 * (ab + sqrt((ab * ab) + (eps * eps))); const BoutReal slope = (ab_pos * (dl + dr)) / (denom + eps); - n.L = n.c - 0.5 * slope; - n.R = n.c + 0.5 * slope; + n.L = n.c - (0.5 * slope); + n.R = n.c + (0.5 * slope); } }; @@ -237,7 +237,7 @@ struct WENO3 { const BoutReal beta1_l = beta0_r; // Smoothness parameter (scaled to local variation) - const BoutReal eps = 1e-12 * (beta0_r + beta1_r) + 1e-30; + const BoutReal eps = (1e-12 * (beta0_r + beta1_r)) + 1e-30; // Linear weights for WENO3-JS constexpr BoutReal d0 = 1.0 / 3.0; @@ -257,8 +257,8 @@ struct WENO3 { const BoutReal w0_l = a0_l / wsum_l; const BoutReal w1_l = a1_l / wsum_l; - n.R = w0_r * p0_r + w1_r * p1_r; - n.L = w0_l * p0_l + w1_l * p1_l; + n.R = (w0_r * p0_r) + (w1_r * p1_r); + n.L = (w0_l * p0_l) + (w1_l * p1_l); } }; From bc980b1366fc6d3feec0c033e21f0c9464de6334 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 24 Jun 2026 12:07:02 +0200 Subject: [PATCH 24/38] Ensure calls to _compute_... are guarded --- include/bout/coordinates.hxx | 140 ++++++++++++++++++++++++----------- 1 file changed, 98 insertions(+), 42 deletions(-) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index 6d803ab28a..a8f4c09c5c 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -116,103 +116,159 @@ public: FieldMetric& g_22_yhigh(); // Cell Areas const FieldMetric& cell_area_xlow() const { - if (!_cell_area_xlow.has_value()) { + if (_cell_area_xlow.has_value()) { + return *_cell_area_xlow; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_x(); + ASSERT2(_cell_area_xlow.has_value()); + return *_cell_area_xlow; } - ASSERT2(_cell_area_xlow.has_value()); - return *_cell_area_xlow; } const FieldMetric& cell_area_xhigh() const { - if (!_cell_area_xhigh.has_value()) { + if (_cell_area_xhigh.has_value()) { + return *_cell_area_xhigh; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_x(); + ASSERT2(_cell_area_xhigh.has_value()); + return *_cell_area_xhigh; } - ASSERT2(_cell_area_xhigh.has_value()); - return *_cell_area_xhigh; } const FieldMetric& cell_area_ylow() const { - if (!_cell_area_ylow.has_value()) { + if (_cell_area_ylow.has_value()) { + return *_cell_area_ylow; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_y(); + ASSERT2(_cell_area_ylow.has_value()); + return *_cell_area_ylow; } - ASSERT2(_cell_area_ylow.has_value()); - return *_cell_area_ylow; } const FieldMetric& cell_area_yhigh() const { - if (!_cell_area_yhigh.has_value()) { + if (_cell_area_yhigh.has_value()) { + return *_cell_area_yhigh; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_y(); + ASSERT2(_cell_area_yhigh.has_value()); + return *_cell_area_yhigh; } - ASSERT2(_cell_area_yhigh.has_value()); - return *_cell_area_yhigh; } const FieldMetric& cell_area_zlow() const { - if (!_cell_area_zlow.has_value()) { + if (_cell_area_zlow.has_value()) { + return *_cell_area_zlow; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_z(); + ASSERT2(_cell_area_zlow.has_value()); + return *_cell_area_zlow; } - ASSERT2(_cell_area_zlow.has_value()); - return *_cell_area_zlow; } const FieldMetric& cell_area_zhigh() const { - if (!_cell_area_zhigh.has_value()) { + if (_cell_area_zhigh.has_value()) { + return *_cell_area_zhigh; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_z(); + ASSERT2(_cell_area_zhigh.has_value()); + return *_cell_area_zhigh; } - ASSERT2(_cell_area_zhigh.has_value()); - return *_cell_area_zhigh; } FieldMetric& cell_area_xlow() { - if (!_cell_area_xlow.has_value()) { + if (_cell_area_xlow.has_value()) { + return *_cell_area_xlow; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_x(); + ASSERT2(_cell_area_xlow.has_value()); + return *_cell_area_xlow; } - ASSERT2(_cell_area_xlow.has_value()); - return *_cell_area_xlow; } FieldMetric& cell_area_xhigh() { - if (!_cell_area_xhigh.has_value()) { + if (_cell_area_xhigh.has_value()) { + return *_cell_area_xhigh; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_x(); + ASSERT2(_cell_area_xhigh.has_value()); + return *_cell_area_xhigh; } - ASSERT2(_cell_area_xhigh.has_value()); - return *_cell_area_xhigh; } FieldMetric& cell_area_ylow() { - if (!_cell_area_ylow.has_value()) { + if (_cell_area_ylow.has_value()) { + return *_cell_area_ylow; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_y(); + ASSERT2(_cell_area_ylow.has_value()); + return *_cell_area_ylow; } - ASSERT2(_cell_area_ylow.has_value()); - return *_cell_area_ylow; } FieldMetric& cell_area_yhigh() { - if (!_cell_area_yhigh.has_value()) { + if (_cell_area_yhigh.has_value()) { + return *_cell_area_yhigh; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_y(); + ASSERT2(_cell_area_yhigh.has_value()); + return *_cell_area_yhigh; } - ASSERT2(_cell_area_yhigh.has_value()); - return *_cell_area_yhigh; } FieldMetric& cell_area_zlow() { - if (!_cell_area_zlow.has_value()) { + if (_cell_area_zlow.has_value()) { + return *_cell_area_zlow; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_z(); + ASSERT2(_cell_area_zlow.has_value()); + return *_cell_area_zlow; } - ASSERT2(_cell_area_zlow.has_value()); - return *_cell_area_zlow; } FieldMetric& cell_area_zhigh() { - if (!_cell_area_zhigh.has_value()) { + if (_cell_area_zhigh.has_value()) { + return *_cell_area_zhigh; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_area_z(); + ASSERT2(_cell_area_zhigh.has_value()); + return *_cell_area_zhigh; } - ASSERT2(_cell_area_zhigh.has_value()); - return *_cell_area_zhigh; } // Cell Volume const FieldMetric& cell_volume() const { - if (!_cell_volume.has_value()) { + if (_cell_volume.has_value()) { + return *_cell_volume; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_volume(); + ASSERT2(_cell_volume.has_value()); + return *_cell_volume; } - ASSERT2(_cell_volume.has_value()); - return *_cell_volume; } FieldMetric& cell_volume() { - if (!_cell_volume.has_value()) { + if (_cell_volume.has_value()) { + return *_cell_volume; + } + BOUT_OMP_SAFE(critical) + { _compute_cell_volume(); + ASSERT2(_cell_volume.has_value()); + return *_cell_volume; } - ASSERT2(_cell_volume.has_value()); - return *_cell_volume; } private: From 9ff3c73ee5ab0f61d9a493b3513558834cce9ab7 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 24 Jun 2026 12:46:20 +0200 Subject: [PATCH 25/38] Add comment about impl header --- include/bout/fv_ops.hxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index 89caaadf97..ab791efb41 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -78,6 +78,8 @@ Field3D D4DY4(const Field3D& d, const Field3D& f); Field3D D4DY4_Index(const Field3D& f, bool bndry_flux = true); // Forward declarations of flux limiters +// If you want to use your own flux limiter, you need to +// #include to instantiate the templates. class Upwind; class Fromm; class MinMod; From b98954f95021b9ad5bbd60795eaab84e571ee2f5 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 24 Jun 2026 12:46:36 +0200 Subject: [PATCH 26/38] Add comments about explicit instantiation --- include/bout/fv_ops_impl.hxx | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx index 91117a6961..81d4c67c25 100644 --- a/include/bout/fv_ops_impl.hxx +++ b/include/bout/fv_ops_impl.hxx @@ -43,6 +43,7 @@ struct Stencil1D { /*! * First order upwind for testing */ +// NB: The templates need to be explicitly instantiated in fv_ops.cxx struct Upwind { void operator()(Stencil1D& n) { n.L = n.R = n.c; } }; @@ -50,6 +51,7 @@ struct Upwind { /*! * Fromm method */ +// NB: The templates need to be explicitly instantiated in fv_ops.cxx struct Fromm { void operator()(Stencil1D& n) { n.L = n.c - (0.25 * (n.p - n.m)); @@ -65,6 +67,7 @@ struct Fromm { * or minimum slope set to zero, i.e. reverts * to first order upwinding */ +// NB: The templates need to be explicitly instantiated in fv_ops.cxx struct MinMod { void operator()(Stencil1D& n) { // Choose the gradient within the cell @@ -102,6 +105,7 @@ private: * 2*right. If any of these slopes have different signs * then the slope reverts to zero (i.e. 1st-order upwinding). */ +// NB: The templates need to be explicitly instantiated in fv_ops.cxx struct MC { void operator()(Stencil1D& n) { const BoutReal slope = minmod(2. * (n.p - n.c), // 2*right difference @@ -149,6 +153,7 @@ private: /// gL = n.c - n.m /// gR = n.p - n.c /// +// NB: The templates need to be explicitly instantiated in fv_ops.cxx struct Superbee { void operator()(Stencil1D& n) { const BoutReal gL = n.c - n.m; @@ -185,6 +190,7 @@ struct Superbee { * * where `pos(x)` is a smooth approximation to `max(x, 0)`. */ +// NB: The templates need to be explicitly instantiated in fv_ops.cxx struct VanAlbada { void operator()(Stencil1D& n) { const BoutReal dl = n.c - n.m; @@ -220,6 +226,7 @@ struct VanAlbada { * Uses only the three-point stencil (`m`, `c`, `p`), so it is a drop-in * replacement anywhere `Stencil1D` is populated with those values. */ +// NB: The templates need to be explicitly instantiated in fv_ops.cxx struct WENO3 { void operator()(Stencil1D& n) { // Right face (between c and p): value from cell c (left state at i+1/2) From e68d1659aca232c6e01278325817dd9c01dff56a Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 24 Jun 2026 12:53:49 +0200 Subject: [PATCH 27/38] Do not return from critical areas --- include/bout/coordinates.hxx | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index a8f4c09c5c..2917e39388 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -123,8 +123,8 @@ public: { _compute_cell_area_x(); ASSERT2(_cell_area_xlow.has_value()); - return *_cell_area_xlow; } + return *_cell_area_xlow; } const FieldMetric& cell_area_xhigh() const { if (_cell_area_xhigh.has_value()) { @@ -134,8 +134,8 @@ public: { _compute_cell_area_x(); ASSERT2(_cell_area_xhigh.has_value()); - return *_cell_area_xhigh; } + return *_cell_area_xhigh; } const FieldMetric& cell_area_ylow() const { if (_cell_area_ylow.has_value()) { @@ -145,8 +145,8 @@ public: { _compute_cell_area_y(); ASSERT2(_cell_area_ylow.has_value()); - return *_cell_area_ylow; } + return *_cell_area_ylow; } const FieldMetric& cell_area_yhigh() const { if (_cell_area_yhigh.has_value()) { @@ -156,8 +156,8 @@ public: { _compute_cell_area_y(); ASSERT2(_cell_area_yhigh.has_value()); - return *_cell_area_yhigh; } + return *_cell_area_yhigh; } const FieldMetric& cell_area_zlow() const { if (_cell_area_zlow.has_value()) { @@ -167,8 +167,8 @@ public: { _compute_cell_area_z(); ASSERT2(_cell_area_zlow.has_value()); - return *_cell_area_zlow; } + return *_cell_area_zlow; } const FieldMetric& cell_area_zhigh() const { if (_cell_area_zhigh.has_value()) { @@ -178,8 +178,8 @@ public: { _compute_cell_area_z(); ASSERT2(_cell_area_zhigh.has_value()); - return *_cell_area_zhigh; } + return *_cell_area_zhigh; } FieldMetric& cell_area_xlow() { if (_cell_area_xlow.has_value()) { @@ -189,8 +189,8 @@ public: { _compute_cell_area_x(); ASSERT2(_cell_area_xlow.has_value()); - return *_cell_area_xlow; } + return *_cell_area_xlow; } FieldMetric& cell_area_xhigh() { if (_cell_area_xhigh.has_value()) { @@ -200,8 +200,8 @@ public: { _compute_cell_area_x(); ASSERT2(_cell_area_xhigh.has_value()); - return *_cell_area_xhigh; } + return *_cell_area_xhigh; } FieldMetric& cell_area_ylow() { if (_cell_area_ylow.has_value()) { @@ -211,8 +211,8 @@ public: { _compute_cell_area_y(); ASSERT2(_cell_area_ylow.has_value()); - return *_cell_area_ylow; } + return *_cell_area_ylow; } FieldMetric& cell_area_yhigh() { if (_cell_area_yhigh.has_value()) { @@ -222,8 +222,8 @@ public: { _compute_cell_area_y(); ASSERT2(_cell_area_yhigh.has_value()); - return *_cell_area_yhigh; } + return *_cell_area_yhigh; } FieldMetric& cell_area_zlow() { if (_cell_area_zlow.has_value()) { @@ -233,8 +233,8 @@ public: { _compute_cell_area_z(); ASSERT2(_cell_area_zlow.has_value()); - return *_cell_area_zlow; } + return *_cell_area_zlow; } FieldMetric& cell_area_zhigh() { if (_cell_area_zhigh.has_value()) { @@ -244,8 +244,8 @@ public: { _compute_cell_area_z(); ASSERT2(_cell_area_zhigh.has_value()); - return *_cell_area_zhigh; } + return *_cell_area_zhigh; } // Cell Volume const FieldMetric& cell_volume() const { @@ -256,8 +256,8 @@ public: { _compute_cell_volume(); ASSERT2(_cell_volume.has_value()); - return *_cell_volume; } + return *_cell_volume; } FieldMetric& cell_volume() { if (_cell_volume.has_value()) { @@ -267,8 +267,8 @@ public: { _compute_cell_volume(); ASSERT2(_cell_volume.has_value()); - return *_cell_volume; } + return *_cell_volume; } private: From a99536026e6683285a2f1dd3c78029cf18fb1243 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 24 Jun 2026 13:01:08 +0200 Subject: [PATCH 28/38] Add missing header --- src/mesh/difops.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 1eef4d7524..79ffbcc44b 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -35,7 +35,9 @@ #include #include #include +#include #include +#include #include #include #include From 5284cc185075cda283c2424d82e9649a633d9e2d Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 24 Jun 2026 13:12:22 +0200 Subject: [PATCH 29/38] Remove TRACE --- src/mesh/difops.cxx | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/mesh/difops.cxx b/src/mesh/difops.cxx index 79ffbcc44b..8ecd9d64ff 100644 --- a/src/mesh/difops.cxx +++ b/src/mesh/difops.cxx @@ -35,7 +35,6 @@ #include #include #include -#include #include #include #include @@ -377,8 +376,6 @@ Field3D Div_par_K_Grad_par(const Field3D& kY, const Field3D& f, CELL_LOC outloc) Field3D Div_par_K_Grad_par_mod(const Field3D& Kin, const Field3D& fin, Field3D& flow_ylow, bool bndry_flux) { - TRACE("FV::Div_par_K_Grad_par_mod"); - ASSERT2(Kin.getLocation() == fin.getLocation()); const Mesh* mesh = Kin.getMesh(); From 1672208242c2bf2c6428129c6945efa49cb24f3d Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 24 Jun 2026 14:10:57 +0200 Subject: [PATCH 30/38] Move critical section to compute functions --- include/bout/coordinates.hxx | 98 +++++------------ src/mesh/coordinates.cxx | 202 ++++++++++++++++++++--------------- 2 files changed, 146 insertions(+), 154 deletions(-) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index 2917e39388..a21f62cd4d 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -119,132 +119,96 @@ public: if (_cell_area_xlow.has_value()) { return *_cell_area_xlow; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_x(); - ASSERT2(_cell_area_xlow.has_value()); - } + _compute_cell_area_x(); + ASSERT2(_cell_area_xlow.has_value()); return *_cell_area_xlow; } const FieldMetric& cell_area_xhigh() const { if (_cell_area_xhigh.has_value()) { return *_cell_area_xhigh; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_x(); - ASSERT2(_cell_area_xhigh.has_value()); - } + _compute_cell_area_x(); + ASSERT2(_cell_area_xhigh.has_value()); return *_cell_area_xhigh; } const FieldMetric& cell_area_ylow() const { if (_cell_area_ylow.has_value()) { return *_cell_area_ylow; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_y(); - ASSERT2(_cell_area_ylow.has_value()); - } + _compute_cell_area_y(); + ASSERT2(_cell_area_ylow.has_value()); return *_cell_area_ylow; } const FieldMetric& cell_area_yhigh() const { if (_cell_area_yhigh.has_value()) { return *_cell_area_yhigh; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_y(); - ASSERT2(_cell_area_yhigh.has_value()); - } + _compute_cell_area_y(); + ASSERT2(_cell_area_yhigh.has_value()); return *_cell_area_yhigh; } const FieldMetric& cell_area_zlow() const { if (_cell_area_zlow.has_value()) { return *_cell_area_zlow; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_z(); - ASSERT2(_cell_area_zlow.has_value()); - } + _compute_cell_area_z(); + ASSERT2(_cell_area_zlow.has_value()); return *_cell_area_zlow; } const FieldMetric& cell_area_zhigh() const { if (_cell_area_zhigh.has_value()) { return *_cell_area_zhigh; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_z(); - ASSERT2(_cell_area_zhigh.has_value()); - } + _compute_cell_area_z(); + ASSERT2(_cell_area_zhigh.has_value()); return *_cell_area_zhigh; } FieldMetric& cell_area_xlow() { if (_cell_area_xlow.has_value()) { return *_cell_area_xlow; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_x(); - ASSERT2(_cell_area_xlow.has_value()); - } + _compute_cell_area_x(); + ASSERT2(_cell_area_xlow.has_value()); return *_cell_area_xlow; } FieldMetric& cell_area_xhigh() { if (_cell_area_xhigh.has_value()) { return *_cell_area_xhigh; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_x(); - ASSERT2(_cell_area_xhigh.has_value()); - } + _compute_cell_area_x(); + ASSERT2(_cell_area_xhigh.has_value()); return *_cell_area_xhigh; } FieldMetric& cell_area_ylow() { if (_cell_area_ylow.has_value()) { return *_cell_area_ylow; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_y(); - ASSERT2(_cell_area_ylow.has_value()); - } + _compute_cell_area_y(); + ASSERT2(_cell_area_ylow.has_value()); return *_cell_area_ylow; } FieldMetric& cell_area_yhigh() { if (_cell_area_yhigh.has_value()) { return *_cell_area_yhigh; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_y(); - ASSERT2(_cell_area_yhigh.has_value()); - } + _compute_cell_area_y(); + ASSERT2(_cell_area_yhigh.has_value()); return *_cell_area_yhigh; } FieldMetric& cell_area_zlow() { if (_cell_area_zlow.has_value()) { return *_cell_area_zlow; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_z(); - ASSERT2(_cell_area_zlow.has_value()); - } + _compute_cell_area_z(); + ASSERT2(_cell_area_zlow.has_value()); return *_cell_area_zlow; } FieldMetric& cell_area_zhigh() { if (_cell_area_zhigh.has_value()) { return *_cell_area_zhigh; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_area_z(); - ASSERT2(_cell_area_zhigh.has_value()); - } + _compute_cell_area_z(); + ASSERT2(_cell_area_zhigh.has_value()); return *_cell_area_zhigh; } // Cell Volume @@ -252,22 +216,16 @@ public: if (_cell_volume.has_value()) { return *_cell_volume; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_volume(); - ASSERT2(_cell_volume.has_value()); - } + _compute_cell_volume(); + ASSERT2(_cell_volume.has_value()); return *_cell_volume; } FieldMetric& cell_volume() { if (_cell_volume.has_value()) { return *_cell_volume; } - BOUT_OMP_SAFE(critical) - { - _compute_cell_volume(); - ASSERT2(_cell_volume.has_value()); - } + _compute_cell_volume(); + ASSERT2(_cell_volume.has_value()); return *_cell_volume; } diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 057ffa65b5..afb63f981d 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -2015,17 +2015,23 @@ const Coordinates::FieldMetric& Coordinates::g_22_ylow() const { if (_g_22_ylow.has_value()) { return *_g_22_ylow; } - _g_22_ylow.emplace(emptyFrom(g_22)); - //_g_22_ylow->setLocation(CELL_YLOW); - auto* mesh = Bxy.getMesh(); - if (Bxy.isFci()) { - if (mesh->get(_g_22_ylow.value(), "g_22_cell_ylow", 0.0, false) != 0) { - throw BoutException("The grid file does not contain `g_22_cell_ylow`."); - } - } else { - ASSERT0(mesh->ystart > 0); - BOUT_FOR(i, g_22.getRegion("RGN_NOY")) { - _g_22_ylow.value()[i] = SQ(0.5 * (std::sqrt(g_22[i]) + std::sqrt(g_22[i.ym()]))); + BOUT_OMP_SAFE(critical) + { + if (!_g_22_ylow.has_value()) { + _g_22_ylow.emplace(emptyFrom(g_22)); + //_g_22_ylow->setLocation(CELL_YLOW); + auto* mesh = Bxy.getMesh(); + if (Bxy.isFci()) { + if (mesh->get(_g_22_ylow.value(), "g_22_cell_ylow", 0.0, false) != 0) { + throw BoutException("The grid file does not contain `g_22_cell_ylow`."); + } + } else { + ASSERT0(mesh->ystart > 0); + BOUT_FOR(i, g_22.getRegion("RGN_NOY")) { + _g_22_ylow.value()[i] = + SQ(0.5 * (std::sqrt(g_22[i]) + std::sqrt(g_22[i.ym()]))); + } + } } } return g_22_ylow(); @@ -2035,99 +2041,127 @@ const Coordinates::FieldMetric& Coordinates::g_22_yhigh() const { if (_g_22_yhigh.has_value()) { return *_g_22_yhigh; } - _g_22_yhigh.emplace(emptyFrom(g_22)); - auto* mesh = Bxy.getMesh(); - if (Bxy.isFci()) { - if (mesh->get(_g_22_yhigh.value(), "g_22_cell_yhigh", 0.0, false) != 0) { - throw BoutException("The grid file does not contain `g_22_cell_yhigh`."); - } - } else { - ASSERT0(mesh->ystart > 0); - BOUT_FOR(i, g_22.getRegion("RGN_NOY")) { - _g_22_yhigh.value()[i] = SQ(0.5 * (std::sqrt(g_22[i]) + std::sqrt(g_22[i.yp()]))); + BOUT_OMP_SAFE(critical) + { + if (!_g_22_yhigh.has_value()) { + _g_22_yhigh.emplace(emptyFrom(g_22)); + auto* mesh = Bxy.getMesh(); + if (Bxy.isFci()) { + if (mesh->get(_g_22_yhigh.value(), "g_22_cell_yhigh", 0.0, false) != 0) { + throw BoutException("The grid file does not contain `g_22_cell_yhigh`."); + } + } else { + ASSERT0(mesh->ystart > 0); + BOUT_FOR(i, g_22.getRegion("RGN_NOY")) { + _g_22_yhigh.value()[i] = + SQ(0.5 * (std::sqrt(g_22[i]) + std::sqrt(g_22[i.yp()]))); + } + } } } return g_22_yhigh(); } void Coordinates::_compute_cell_area_x() const { - const FieldMetric area_centre = sqrt(g_22 * g_33 - SQ(g_23)) * dy * dz; - _cell_area_xlow.emplace(emptyFrom(area_centre)); - _cell_area_xhigh.emplace(emptyFrom(area_centre)); - // We cannot setLocation, as that would trigger the computation of staggered - // metrics. - auto* mesh = Bxy.getMesh(); - ASSERT0(mesh->xstart > 0); - BOUT_FOR(i, area_centre.getRegion("RGN_NOX")) { - (*_cell_area_xlow)[i] = 0.5 * (area_centre[i] + area_centre[i.xm()]); - (*_cell_area_xhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.xp()]); + BOUT_OMP_SAFE(critical) + { + if (!_cell_area_xlow.has_value()) { + const FieldMetric area_centre = sqrt(g_22 * g_33 - SQ(g_23)) * dy * dz; + _cell_area_xlow.emplace(emptyFrom(area_centre)); + _cell_area_xhigh.emplace(emptyFrom(area_centre)); + // We cannot setLocation, as that would trigger the computation of staggered + // metrics. + auto* mesh = Bxy.getMesh(); + ASSERT0(mesh->xstart > 0); + BOUT_FOR(i, area_centre.getRegion("RGN_NOX")) { + (*_cell_area_xlow)[i] = 0.5 * (area_centre[i] + area_centre[i.xm()]); + (*_cell_area_xhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.xp()]); + } + } } } void Coordinates::_compute_cell_area_y() const { - auto* mesh = Bxy.getMesh(); - if (g_11.isFci()) { - const FieldMetric jxz_centre = sqrt(g_11 * g_33 - SQ(g_13)); - auto jxz_ylow = emptyFrom(jxz_centre); - auto jxz_yhigh = emptyFrom(jxz_centre); - - auto By_c = emptyFrom(jxz_centre); - auto By_h = emptyFrom(jxz_yhigh); - auto By_l = emptyFrom(jxz_ylow); - if (mesh->get(By_c, "By", 0.0, false, CELL_CENTRE) != 0) { - throw BoutException("The grid file does not contain `By`."); - } - if (mesh->get(By_l, "By_cell_ylow", 0.0, false) != 0) { - throw BoutException("The grid file does not contain `By_cell_ylow`."); - } - if (mesh->get(By_h, "By_cell_yhigh", 0.0, false) != 0) { - throw BoutException("The grid file does not contain `By_cell_yhigh`."); - } - BOUT_FOR(i, By_c.getRegion("RGN_NOY")) { - jxz_ylow[i] = By_c[i] / By_l[i] * jxz_centre[i]; - jxz_yhigh[i] = By_c[i] / By_h[i] * jxz_centre[i]; - } - ASSERT3(isUniform(dx, true, "RGN_ALL")); - ASSERT2(isUniform(dx, false, "RGN_ALL")); - ASSERT3(isUniform(dz, true, "RGN_ALL")); - ASSERT2(isUniform(dz, false, "RGN_ALL")); - _cell_area_ylow.emplace(jxz_ylow * dx * dz); - _cell_area_yhigh.emplace(jxz_yhigh * dx * dz); - } else { - // Field aligned - const FieldMetric area_centre = sqrt(g_11 * g_33 - SQ(g_13)) * dx * dz; - _cell_area_ylow.emplace(emptyFrom(area_centre)); - _cell_area_yhigh.emplace(emptyFrom(area_centre)); - // We cannot setLocation, as that would trigger the computation of staggered - // metrics. - BOUT_FOR(i, mesh->getRegion("RGN_ALL")) { - if (i.y() > 0) { - (*_cell_area_ylow)[i] = 0.5 * (area_centre[i] + area_centre[i.ym()]); - } else { - (*_cell_area_ylow)[i] = BoutNaN; - } - if (i.y() < mesh->LocalNy - 1) { - (*_cell_area_yhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.yp()]); + BOUT_OMP_SAFE(critical) + { + if (!_cell_area_ylow.has_value()) { + auto* mesh = Bxy.getMesh(); + if (g_11.isFci()) { + const FieldMetric jxz_centre = sqrt(g_11 * g_33 - SQ(g_13)); + auto jxz_ylow = emptyFrom(jxz_centre); + auto jxz_yhigh = emptyFrom(jxz_centre); + + auto By_c = emptyFrom(jxz_centre); + auto By_h = emptyFrom(jxz_yhigh); + auto By_l = emptyFrom(jxz_ylow); + if (mesh->get(By_c, "By", 0.0, false, CELL_CENTRE) != 0) { + throw BoutException("The grid file does not contain `By`."); + } + if (mesh->get(By_l, "By_cell_ylow", 0.0, false) != 0) { + throw BoutException("The grid file does not contain `By_cell_ylow`."); + } + if (mesh->get(By_h, "By_cell_yhigh", 0.0, false) != 0) { + throw BoutException("The grid file does not contain `By_cell_yhigh`."); + } + BOUT_FOR(i, By_c.getRegion("RGN_NOY")) { + jxz_ylow[i] = By_c[i] / By_l[i] * jxz_centre[i]; + jxz_yhigh[i] = By_c[i] / By_h[i] * jxz_centre[i]; + } + ASSERT3(isUniform(dx, true, "RGN_ALL")); + ASSERT2(isUniform(dx, false, "RGN_ALL")); + ASSERT3(isUniform(dz, true, "RGN_ALL")); + ASSERT2(isUniform(dz, false, "RGN_ALL")); + _cell_area_ylow.emplace(jxz_ylow * dx * dz); + _cell_area_yhigh.emplace(jxz_yhigh * dx * dz); } else { - (*_cell_area_yhigh)[i] = BoutNaN; + // Field aligned + const FieldMetric area_centre = sqrt(g_11 * g_33 - SQ(g_13)) * dx * dz; + _cell_area_ylow.emplace(emptyFrom(area_centre)); + _cell_area_yhigh.emplace(emptyFrom(area_centre)); + // We cannot setLocation, as that would trigger the computation of staggered + // metrics. + BOUT_FOR(i, mesh->getRegion("RGN_ALL")) { + if (i.y() > 0) { + (*_cell_area_ylow)[i] = 0.5 * (area_centre[i] + area_centre[i.ym()]); + } else { + (*_cell_area_ylow)[i] = BoutNaN; + } + if (i.y() < mesh->LocalNy - 1) { + (*_cell_area_yhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.yp()]); + } else { + (*_cell_area_yhigh)[i] = BoutNaN; + } + } } } } } void Coordinates::_compute_cell_area_z() const { - const FieldMetric area_centre = sqrt(g_11 * g_22 - SQ(g_12)) * dx * dy; - _cell_area_zlow.emplace(emptyFrom(area_centre)); - _cell_area_zhigh.emplace(emptyFrom(area_centre)); - // We cannot setLocation, as that would trigger the computation of staggered - // metrics. - BOUT_FOR(i, area_centre.getRegion("RGN_NOZ")) { - (*_cell_area_zlow)[i] = 0.5 * (area_centre[i] + area_centre[i.zm()]); - (*_cell_area_zhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.zp()]); + BOUT_OMP_SAFE(critical) + { + if (!_cell_volume.has_value()) { + const FieldMetric area_centre = sqrt(g_11 * g_22 - SQ(g_12)) * dx * dy; + _cell_area_zlow.emplace(emptyFrom(area_centre)); + _cell_area_zhigh.emplace(emptyFrom(area_centre)); + // We cannot setLocation, as that would trigger the computation of staggered + // metrics. + BOUT_FOR(i, area_centre.getRegion("RGN_NOZ")) { + (*_cell_area_zlow)[i] = 0.5 * (area_centre[i] + area_centre[i.zm()]); + (*_cell_area_zhigh)[i] = 0.5 * (area_centre[i] + area_centre[i.zp()]); + } + } } } -void Coordinates::_compute_cell_volume() const { _cell_volume.emplace(J * dx * dy * dz); } +void Coordinates::_compute_cell_volume() const { + BOUT_OMP_SAFE(critical) + { + if (!_cell_volume.has_value()) { + _cell_volume.emplace(J * dx * dy * dz); + } + } +} std::shared_ptr Coordinates::makeYBoundary(YBndryType type) const { return std::make_shared(type, localoptions, *localmesh); From b89867dd3a54d17d1c7b4174a79baabebec44c54 Mon Sep 17 00:00:00 2001 From: David Bold Date: Wed, 24 Jun 2026 14:11:38 +0200 Subject: [PATCH 31/38] Add DDY template for expression --- include/bout/derivs.hxx | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/include/bout/derivs.hxx b/include/bout/derivs.hxx index a8d9279378..92f094f550 100644 --- a/include/bout/derivs.hxx +++ b/include/bout/derivs.hxx @@ -86,6 +86,13 @@ Field3D DDY(const Field3DParallel& f, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", const std::string& region = "RGN_NOBNDRY"); +template +std::enable_if_t && !bout::utils::is_Field3D_v, const Field3D> +DDY(const E& expr, CELL_LOC outloc = CELL_DEFAULT, const std::string& method = "DEFAULT", + const std::string& region = "RGN_NOBNDRY") { + return DDY(Field3D{expr}, outloc, method, region); +} + /// Calculate first partial derivative in Y /// /// \f$\partial / \partial y\f$ From 26876b7d6bf786b6c47dcc9ebd06b009939db8fc Mon Sep 17 00:00:00 2001 From: David Bold Date: Fri, 26 Jun 2026 14:09:49 +0200 Subject: [PATCH 32/38] Add missing used header --- include/bout/derivs.hxx | 1 + 1 file changed, 1 insertion(+) diff --git a/include/bout/derivs.hxx b/include/bout/derivs.hxx index 92f094f550..14bc5c2824 100644 --- a/include/bout/derivs.hxx +++ b/include/bout/derivs.hxx @@ -35,6 +35,7 @@ #include "bout/vector3d.hxx" #include "bout/bout_types.hxx" +#include ////////// FIRST DERIVATIVES ////////// From bf5f01df35c6c57b4af507d1b22b8e0a9b7a5da6 Mon Sep 17 00:00:00 2001 From: David Bold Date: Mon, 29 Jun 2026 09:07:12 +0200 Subject: [PATCH 33/38] Remove unused headers --- include/bout/fv_ops.hxx | 3 --- tests/MMS/spatial/finite-volume/fv_mms.cxx | 1 - 2 files changed, 4 deletions(-) diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index ab791efb41..b2b0e49120 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -7,7 +7,6 @@ #include "bout/assert.hxx" #include "bout/bout_types.hxx" -#include "bout/boutexception.hxx" #include "bout/build_defines.hxx" #include "bout/coordinates.hxx" #include "bout/field.hxx" @@ -19,8 +18,6 @@ #include "bout/utils.hxx" #include "bout/vector2d.hxx" -#include - namespace FV { /*! * Div ( a Grad_perp(f) ) -- ∇⊥ ( a ⋅ ∇⊥ f) -- Vorticity diff --git a/tests/MMS/spatial/finite-volume/fv_mms.cxx b/tests/MMS/spatial/finite-volume/fv_mms.cxx index 19f3f14610..1de6052e89 100644 --- a/tests/MMS/spatial/finite-volume/fv_mms.cxx +++ b/tests/MMS/spatial/finite-volume/fv_mms.cxx @@ -7,7 +7,6 @@ #include "bout/globals.hxx" #include "bout/options.hxx" #include "bout/options_io.hxx" -#include "bout/utils.hxx" #include From f893f3da21bbe19bdaba4f5bb17b547d8b2dd326 Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 2 Jul 2026 11:18:56 +0200 Subject: [PATCH 34/38] Fix nonorthogonal MMS test for hermes-3 --- src/mesh/coordinates.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index afb63f981d..10b560f71c 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -2066,7 +2066,7 @@ void Coordinates::_compute_cell_area_x() const { BOUT_OMP_SAFE(critical) { if (!_cell_area_xlow.has_value()) { - const FieldMetric area_centre = sqrt(g_22 * g_33 - SQ(g_23)) * dy * dz; + const FieldMetric area_centre = J / sqrt(g_11) * dy * dz; _cell_area_xlow.emplace(emptyFrom(area_centre)); _cell_area_xhigh.emplace(emptyFrom(area_centre)); // We cannot setLocation, as that would trigger the computation of staggered @@ -2087,7 +2087,7 @@ void Coordinates::_compute_cell_area_y() const { if (!_cell_area_ylow.has_value()) { auto* mesh = Bxy.getMesh(); if (g_11.isFci()) { - const FieldMetric jxz_centre = sqrt(g_11 * g_33 - SQ(g_13)); + const FieldMetric jxz_centre = J / sqrt(g_22); auto jxz_ylow = emptyFrom(jxz_centre); auto jxz_yhigh = emptyFrom(jxz_centre); @@ -2115,7 +2115,7 @@ void Coordinates::_compute_cell_area_y() const { _cell_area_yhigh.emplace(jxz_yhigh * dx * dz); } else { // Field aligned - const FieldMetric area_centre = sqrt(g_11 * g_33 - SQ(g_13)) * dx * dz; + const FieldMetric area_centre = J / sqrt(g_22) * dx * dz; _cell_area_ylow.emplace(emptyFrom(area_centre)); _cell_area_yhigh.emplace(emptyFrom(area_centre)); // We cannot setLocation, as that would trigger the computation of staggered @@ -2141,7 +2141,7 @@ void Coordinates::_compute_cell_area_z() const { BOUT_OMP_SAFE(critical) { if (!_cell_volume.has_value()) { - const FieldMetric area_centre = sqrt(g_11 * g_22 - SQ(g_12)) * dx * dy; + const FieldMetric area_centre = J / sqrt(g_33) * dx * dy; _cell_area_zlow.emplace(emptyFrom(area_centre)); _cell_area_zhigh.emplace(emptyFrom(area_centre)); // We cannot setLocation, as that would trigger the computation of staggered From 6a0e276028fbdd2e817f5a67813e64c9a1417850 Mon Sep 17 00:00:00 2001 From: David Bold Date: Thu, 2 Jul 2026 11:52:10 +0200 Subject: [PATCH 35/38] Ensure fv_ops.hxx is parsed before fv_ops_impl.hxx We need to ensure that fv_ops_impl.hxx is never parsed before fv_ops.hxx, as fv_ops.hxx provides default arguments, which is forbiddef for a redeclaration. To ensure this, fv_ops.hxx is included in fv_ops_impl.hxx, so that fv_ops_impl.hxx can be included and fv_ops.hxx can appear afterwards without causing an error. --- include/bout/fv_ops_impl.hxx | 1 + 1 file changed, 1 insertion(+) diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx index 67efdd6476..3755467027 100644 --- a/include/bout/fv_ops_impl.hxx +++ b/include/bout/fv_ops_impl.hxx @@ -12,6 +12,7 @@ #include "bout/coordinates.hxx" #include "bout/field.hxx" #include "bout/field3d.hxx" +#include "bout/fv_ops.hxx" // NOLINT(unused-includes, misc-include-cleaner) #include "bout/globals.hxx" #include "bout/mesh.hxx" #include "bout/output_bout_types.hxx" // NOLINT(unused-includes, misc-include-cleaner) From 0b187d94572086b6fae3a376e43e375a72f1ecdc Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 4 Jul 2026 08:54:01 -0700 Subject: [PATCH 36/38] Coordinates::_compute_cell_area_z check _compute_cell_area_z Was testing _cell_volume.has_value() rather than _cell_area_zlow.has_value(). Added regression test. --- src/mesh/coordinates.cxx | 2 +- tests/unit/mesh/test_coordinates.cxx | 27 +++++++++++++++++++++++++++ 2 files changed, 28 insertions(+), 1 deletion(-) diff --git a/src/mesh/coordinates.cxx b/src/mesh/coordinates.cxx index 10b560f71c..e7a7579d45 100644 --- a/src/mesh/coordinates.cxx +++ b/src/mesh/coordinates.cxx @@ -2140,7 +2140,7 @@ void Coordinates::_compute_cell_area_y() const { void Coordinates::_compute_cell_area_z() const { BOUT_OMP_SAFE(critical) { - if (!_cell_volume.has_value()) { + if (!_cell_area_zlow.has_value()) { const FieldMetric area_centre = J / sqrt(g_33) * dx * dy; _cell_area_zlow.emplace(emptyFrom(area_centre)); _cell_area_zhigh.emplace(emptyFrom(area_centre)); diff --git a/tests/unit/mesh/test_coordinates.cxx b/tests/unit/mesh/test_coordinates.cxx index ea233cd9ca..8e1aa2e13f 100644 --- a/tests/unit/mesh/test_coordinates.cxx +++ b/tests/unit/mesh/test_coordinates.cxx @@ -374,3 +374,30 @@ TEST_F(CoordinatesTest, CellAreasUpdate) { EXPECT_TRUE(IsFieldEqual(coords.cell_volume(), 8.0)); } + +TEST_F(CoordinatesTest, CellAreaZComputedAfterCellVolume) { + Coordinates coords{mesh, + FieldMetric{1.0}, // dx + FieldMetric{1.0}, // dy + FieldMetric{1.0}, // dz + FieldMetric{6.0}, // J + FieldMetric{1.0}, // Bxy + FieldMetric{1.0}, // g11 + FieldMetric{1.0}, // g22 + FieldMetric{9.0}, // g33 + FieldMetric{0.0}, // g12 + FieldMetric{0.0}, // g13 + FieldMetric{0.0}, // g23 + FieldMetric{4.0}, // g_11 + FieldMetric{1.0}, // g_22 + FieldMetric{9.0}, // g_33 + FieldMetric{0.0}, // g_12 + FieldMetric{0.0}, // g_13 + FieldMetric{0.0}, // g_23 + FieldMetric{0.0}, // ShiftTorsion + FieldMetric{0.0}}; // IntShiftTorsion + + EXPECT_TRUE(IsFieldEqual(coords.cell_volume(), 6.0)); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_zlow(), 2.0, "RGN_NOZ")); + EXPECT_TRUE(IsFieldEqual(coords.cell_area_zhigh(), 2.0, "RGN_NOZ")); +} From 2c41fcabfcffe7e5197f2609aaa2bba2dcb9823d Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 4 Jul 2026 08:55:48 -0700 Subject: [PATCH 37/38] FV ops: Document diagnostics for FCI args Field-aligned operation returns a flow diagnostic, but this is not currently done for FCI. Set to zero and document behavior. --- include/bout/coordinates.hxx | 9 +-------- include/bout/difops.hxx | 1 + include/bout/fv_ops.hxx | 1 + include/bout/fv_ops_impl.hxx | 1 + tests/MMS/spatial/finite-volume/makefile | 2 +- 5 files changed, 5 insertions(+), 9 deletions(-) diff --git a/include/bout/coordinates.hxx b/include/bout/coordinates.hxx index a21f62cd4d..14019728d0 100644 --- a/include/bout/coordinates.hxx +++ b/include/bout/coordinates.hxx @@ -1,15 +1,8 @@ /************************************************************************** * Describes coordinate systems * - * ChangeLog - * ========= - * - * 2014-11-10 Ben Dudson - * * Created by separating metric from Mesh - * - * ************************************************************************** - * Copyright 2014-2025 BOUT++ contributors + * Copyright 2014-2026 BOUT++ contributors * * Contact: Ben Dudson, dudson2@llnl.gov * diff --git a/include/bout/difops.hxx b/include/bout/difops.hxx index 2070cc30d3..0754c1feb6 100644 --- a/include/bout/difops.hxx +++ b/include/bout/difops.hxx @@ -196,6 +196,7 @@ Field3D Div_par_K_Grad_par(const Field3D& kY, const Field3D& f, CELL_LOC outloc = CELL_DEFAULT); /// Version with energy flow diagnostic +/// For FCI fields, `flow_ylow` is currently returned as zero. Field3D Div_par_K_Grad_par_mod(const Field3D& k, const Field3D& f, Field3D& flow_ylow, bool bndry_flux = true); diff --git a/include/bout/fv_ops.hxx b/include/bout/fv_ops.hxx index b2b0e49120..94c599c58c 100644 --- a/include/bout/fv_ops.hxx +++ b/include/bout/fv_ops.hxx @@ -151,6 +151,7 @@ Field3D Div_Perp_Lap(const Field3D& a, const Field3D& f, CELL_LOC outloc = CELL_ /// /// @param[out] flow_ylow Flow at the lower Y cell boundary /// Already includes area factor * flux +/// For FCI fields this diagnostic is currently set to zero. template Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, const Field3D& wave_speed_in, Field3D& flow_ylow, diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx index 3755467027..7da980e584 100644 --- a/include/bout/fv_ops_impl.hxx +++ b/include/bout/fv_ops_impl.hxx @@ -635,6 +635,7 @@ Field3D Div_par_mod(const Field3D& f_in, const Field3D& v_in, const auto& v_down = v_in.ydown(); Field3D result{emptyFrom(f_in)}; + flow_ylow = zeroFrom(f_in); BOUT_FOR(i, f_in.getRegion("RGN_NOBNDRY")) { const auto iyp = i.yp(); const auto iym = i.ym(); diff --git a/tests/MMS/spatial/finite-volume/makefile b/tests/MMS/spatial/finite-volume/makefile index 88ba6c77e7..7111f5ed8e 100644 --- a/tests/MMS/spatial/finite-volume/makefile +++ b/tests/MMS/spatial/finite-volume/makefile @@ -1,6 +1,6 @@ BOUT_TOP = ../../../.. -SOURCEC = fci_mms.cxx +SOURCEC = fv_mms.cxx include $(BOUT_TOP)/make.config From 8fb81ad62617baeb304b2a6478d277d2717f4e81 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sat, 4 Jul 2026 10:29:14 -0700 Subject: [PATCH 38/38] Documentation of FV operators and metric caching Short descriptions of the ported operators, SuperBee slope limiter, and metric derived quantity caching. --- manual/sphinx/developer_docs/mesh.rst | 23 +++++++ .../user_docs/differential_operators.rst | 65 +++++++++++++++++++ 2 files changed, 88 insertions(+) diff --git a/manual/sphinx/developer_docs/mesh.rst b/manual/sphinx/developer_docs/mesh.rst index ff3e40d4b0..b1861fc416 100644 --- a/manual/sphinx/developer_docs/mesh.rst +++ b/manual/sphinx/developer_docs/mesh.rst @@ -307,3 +307,26 @@ because they are needed in a lot of the code. They shouldn’t change after initialisation, unless the physics model starts doing fancy things with deforming meshes. In that case it is up to the user to ensure they are updated. + +.. _sec-derived-geometric-quantities: + +Derived geometric quantities +~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +`Coordinates` also provides several quantities derived from the metric tensor +and Jacobian, rather than requiring each operator to reconstruct them locally. +Examples include: + +* `g_22` evaluated at the lower and upper `y` cell faces +* cell-face areas in the `x`, `y`, and `z` directions +* cell volumes + +These are used by conservative operators, especially the finite-volume +operators documented in :ref:`sec-finite-volume-operators`, where fluxes are +naturally expressed as a face area multiplied by a face flux and divided by a +cell volume. + +In the current implementation these quantities are computed lazily from the +current metric data and then cached for reuse. This keeps the operator code +closer to the discrete flux expressions while centralising the geometry +construction in `Coordinates`. diff --git a/manual/sphinx/user_docs/differential_operators.rst b/manual/sphinx/user_docs/differential_operators.rst index 71dbeb4f9e..3de13001ab 100644 --- a/manual/sphinx/user_docs/differential_operators.rst +++ b/manual/sphinx/user_docs/differential_operators.rst @@ -503,6 +503,8 @@ neglected if :math:`g_{xy}` and :math:`g_{yz}` are non-zero. An example of usage of the brackets can be found in for example ``examples/MMS/advection`` or ``examples/blob2d``. +.. _sec-finite-volume-operators: + Finite volume, conservative finite difference methods ----------------------------------------------------- @@ -537,6 +539,8 @@ The methods can be used by including the Some methods (those with templates) are defined in the header, but others are defined in :doc:`src/mesh/fv_ops.cxx<../_breathe_autogen/file/fv__ops_8cxx>`. +These operators use derived geometric quantities such as cell-face areas and +cell volumes; see :ref:`sec-derived-geometric-quantities`. Parallel divergence ``Div_par`` @@ -570,6 +574,45 @@ be changed at compile time e.g:: A list of available limiters is given in section :ref:`sec-slope-limiters` below. +Modified parallel divergence ``Div_par_mod`` +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This is a modified version of ``FV::Div_par`` which applies the limiter to both +the advected field and the velocity. This typically gives smaller overshoots +than ``FV::Div_par`` for the same limiter choice. + +:: + + template + Field3D Div_par_mod(const Field3D &f_in, const Field3D &v_in, + const Field3D &a, Field3D &flow_ylow, + bool fixflux=true); + + +The extra output argument ``flow_ylow`` stores the flow through the lower +:math:`y` cell boundary, including the area factor. This can be useful as a +diagnostic in energy or flux budgets. For FCI fields this diagnostic is +currently returned as zero. + + +Parallel momentum flux ``Div_par_fvv`` +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This operator calculates the divergence of :math:`f v v`, and is mainly useful +for the parallel momentum equation: + +:: + + template + Field3D Div_par_fvv(const Field3D &f_in, const Field3D &v_in, + const Field3D &a, bool fixflux=true); + + +This is provided separately rather than forming :math:`fv` first, so that the +reconstructed values of :math:`f` and :math:`v` remain consistent with the +other finite-volume advection operators. + + Example and convergence test ++++++++++++++++++++++++++++ @@ -620,6 +663,24 @@ This is done by calculating the flux :math:`k\partial_{||}\left(f\right)` on cel using central differencing. +Parallel diffusion with flow diagnostic ``Div_par_K_Grad_par_mod`` +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This variant of the parallel diffusion operator also returns a lower-face flow +diagnostic: + +:: + + Field3D Div_par_K_Grad_par_mod(const Field3D &k, const Field3D &f, + Field3D &flow_ylow, + bool bndry_flux=true); + + +As for ``Div_par_mod``, ``flow_ylow`` stores the lower :math:`y` boundary flow +including the area factor. For FCI fields this diagnostic is currently returned +as zero. + + Advection in 3D ~~~~~~~~~~~~~~~ @@ -667,6 +728,10 @@ values. Several slope limiters are defined in ``fv_ops.hxx``: to ``MinMod``. It has smaller dissipation than ``MinMod`` so is the default. +* ``Superbee`` - A more compressive TVD limiter than ``MC`` or ``MinMod``. + It tends to sharpen steep gradients and contacts more aggressively, at the + cost of being less smooth. + * ``VanAlbada`` - A smooth (differentiable) symmetric slope limiter which avoids piecewise branches at extrema. This can be useful for nonlinear solvers and finite-difference Jacobian calculations.