From fd11e55793c44b13319c80542df29d1407c64455 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 27 Oct 2026 09:35:39 -0700 Subject: [PATCH 1/5] Div_par_fvv_heating: conserve kinetic and thermal energy Implementation changes: - Use the transport operators' split particle flux for numerical heating, correcting the doubled wave-speed contribution. - Match open-boundary momentum fluxes to Div_par_fvv in both flux modes. - Use the receiving cell's velocity for left-going kinetic energy flux. - Allow configurable Y guards in the fake mesh and add energy-balance regressions for MC, MinMod and Upwind reconstruction. Documentation changes: - Explain the discrete kinetic-energy identity and flux orientation in operator comments, and document the two-guard test fixture. Tests and validation: - Nine energy-balance cases failed before the fix; all 12 new cases pass. - Check local and integrated energy balance with nonuniform geometry, periodic flow, both open-boundary flux modes and uniform velocity. - Full BOUT serial suite: 2017 passed, four skipped and three failed. - Changed code passes clang-format checks and git diff --check. Notes: - Changes apply to the non-FCI implementation. - Two full-suite failures require PVODE, disabled in the isolated build; the remaining failure is the signal-handler death test. - Hermes' Python reconstruction and tracked heating-defect regression still need updating for the corrected operator. --- include/bout/fv_ops_impl.hxx | 28 +++---- tests/unit/fake_mesh.hxx | 14 ++-- tests/unit/fake_mesh_fixture.hxx | 8 +- tests/unit/mesh/test_fv_ops.cxx | 123 +++++++++++++++++++++++++++++++ 4 files changed, 151 insertions(+), 22 deletions(-) diff --git a/include/bout/fv_ops_impl.hxx b/include/bout/fv_ops_impl.hxx index 67b08664c8..5ac81d4e35 100644 --- a/include/bout/fv_ops_impl.hxx +++ b/include/bout/fv_ops_impl.hxx @@ -1170,7 +1170,7 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in, } else { flux_mom = (s.R * sv.R * sv.R) + (BOUTMAX(wave_speed(i, j, k), std::abs(sv.c), std::abs(sv.p)) - * (s.R * sv.R - n_mid * v_mid)); + * n_mid * (sv.R - v_mid)); } // Assume that particle flux is fixed to boundary value @@ -1191,14 +1191,16 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in, const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j + 1, k), std::abs(sv.c), std::abs(sv.p)); - // Viscous heating due to relaxation of velocity towards midpoint - result(i, j, k) += - (amax + 0.5 * sv.R) * s.R * (sv.c - sv.p) * (sv.R - v_mid) * flux_factor_rc; + // Use the same split particle flux as Div_par_mod and Div_par_fvv. + const BoutReal flux_part = s.R * 0.5 * (sv.R + amax); + + // Difference between kinetic energy leaving this cell and entering + // the next cell, using d(n*v^2/2)/dt = v*d(n*v)/dt - v^2*dn/dt/2. + result(i, j, k) += flux_part * (sv.c - sv.p) * (sv.R - v_mid) * flux_factor_rc; // Kinetic energy flow into next cell. // Note: Different from flow out of this cell; the difference // is in the viscous heating. - const BoutReal flux_part = s.R * 0.5 * (sv.R + amax); const BoutReal flux_mom = flux_part * sv.R; flow_ylow(i, j + 1, k) += @@ -1224,7 +1226,7 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in, // Add flux due to difference in boundary values flux_mom = (s.L * sv.L * sv.L) - (BOUTMAX(wave_speed(i, j, k), std::abs(sv.c), std::abs(sv.m)) - * (s.L * sv.L - n_mid * v_mid)); + * n_mid * (sv.L - v_mid)); } // Assume that density flux is fixed to boundary value @@ -1241,17 +1243,17 @@ Field3D Div_par_fvv_heating(const Field3D& f_in, const Field3D& v_in, const BoutReal amax = BOUTMAX(wave_speed(i, j, k), wave_speed(i, j - 1, k), std::abs(sv.c), std::abs(sv.m)); - // Viscous heating due to relaxation - result(i, j, k) += - (amax - 0.5 * sv.L) * s.L * (sv.c - sv.m) * (sv.L - v_mid) * flux_factor_lc; + const BoutReal flux_part = s.L * 0.5 * (sv.L - amax); - // Kinetic energy flow into this cell. - // Note: Different from flow out of left cell; the difference + // The left-going split flux has the opposite orientation. + result(i, j, k) -= flux_part * (sv.c - sv.m) * (sv.L - v_mid) * flux_factor_lc; + + // Kinetic energy flow into the previous cell. + // Note: Different from flow out of this cell; the difference // is in the viscous heating. - const BoutReal flux_part = s.L * 0.5 * (sv.L - amax); const BoutReal flux_mom = flux_part * sv.L; - flow_ylow(i, j, k) += (sv.c * flux_mom - 0.5 * SQ(sv.c) * flux_part) * area_lc; + flow_ylow(i, j, k) += (sv.m * flux_mom - 0.5 * SQ(sv.m) * flux_part) * area_lc; } } } diff --git a/tests/unit/fake_mesh.hxx b/tests/unit/fake_mesh.hxx index 216d04e18a..81e116f4aa 100644 --- a/tests/unit/fake_mesh.hxx +++ b/tests/unit/fake_mesh.hxx @@ -35,20 +35,20 @@ class Options; /// - This is a mesh for a single process, so the global and local /// indices are the same. /// -/// - There is a single guard cell at each of the start/end x/y grids. +/// - There is one X guard cell and, by default, one Y guard cell at each end. /// /// - Only the **grid** information is assumed to be used -- anything /// else will likely **not** work! class FakeMesh : public Mesh { public: - FakeMesh(int nx, int ny, int nz, MpiWrapper& mpi_in) { + FakeMesh(int nx, int ny, int nz, MpiWrapper& mpi_in, int myg = 1) { // Mesh only on one process, so global and local indices are the // same GlobalNx = nx; GlobalNy = ny; GlobalNz = nz; GlobalNxNoBoundaries = nx - 2; - GlobalNyNoBoundaries = ny - 2; + GlobalNyNoBoundaries = ny - 2 * myg; GlobalNzNoBoundaries = nz; LocalNx = nx; LocalNy = ny; @@ -59,20 +59,20 @@ public: // These bits only for ADIOS2, also boring due to single process MapCountX = nx - 2; - MapCountY = ny - 2; + MapCountY = ny - 2 * myg; MapCountZ = nz; MapGlobalX = nx; MapGlobalY = ny; MapGlobalZ = nz; MapLocalX = nx - 2; - MapLocalY = ny - 2; + MapLocalY = ny - 2 * myg; MapLocalZ = nz; // Small "inner" region xstart = 1; xend = nx - 2; - ystart = 1; - yend = ny - 2; + ystart = myg; + yend = ny - myg - 1; zstart = 0; // no guards zend = nz - 1; diff --git a/tests/unit/fake_mesh_fixture.hxx b/tests/unit/fake_mesh_fixture.hxx index 021f37a24d..a2b8417418 100644 --- a/tests/unit/fake_mesh_fixture.hxx +++ b/tests/unit/fake_mesh_fixture.hxx @@ -34,11 +34,15 @@ /// Use this template class directly to use different sized grid: /// /// using MyTest = FakeMeshFixture_tmpl<7, 9, 11>; -template +/// +/// Use the final template argument to retain two Y guards for reconstruction: +/// +/// using ParallelTest = FakeMeshFixture_tmpl<5, 20, 3, false, 2>; +template class FakeMeshFixture_tmpl : public ::testing::Test { public: FakeMeshFixture_tmpl() - : mesh_m(NX, NY, NZ, mpi), mesh_staggered_m(NX, NY, NZ, mpi), + : mesh_m(NX, NY, NZ, mpi, MYG), mesh_staggered_m(NX, NY, NZ, mpi, MYG), mesh_staggered(&mesh_staggered_m) { bout::globals::mpi = &mpi; diff --git a/tests/unit/mesh/test_fv_ops.cxx b/tests/unit/mesh/test_fv_ops.cxx index 2aecd707c4..92a3340413 100644 --- a/tests/unit/mesh/test_fv_ops.cxx +++ b/tests/unit/mesh/test_fv_ops.cxx @@ -5,6 +5,7 @@ #include "bout/constants.hxx" #include "bout/coordinates.hxx" #include "bout/fv_ops.hxx" +#include "bout/fv_ops_impl.hxx" #include "bout/mesh.hxx" #include "bout/output.hxx" @@ -33,3 +34,125 @@ TEST_F(FVOpsTest, Div_a_Grad_perp_FCI_coeffs) { (*op)(a, f, true); (*op)(a, f, false); } + +namespace { + +// Two Y guards let the parallel fluxes reconstruct their neighboring guard cell. +// Keep the existing one-guard fixture as the default for other unit tests. +template +class FVVHeatingTest : public FakeMeshFixture_tmpl<5, 20, 3, false, 2> { +public: + struct Fields { + Field3D density{1.0}, velocity{0.0}, wave{1.0}; + }; + + Fields fields(BoutReal flow_sign, bool uniform_velocity = false) { + Fields result; + const int ny_physical = mesh->yend - mesh->ystart + 1; + Coordinates::FieldMetric jacobian{1.0}, dy{1.0}; + for (int x = 0; x < mesh->LocalNx; ++x) { + for (int y = 0; y < mesh->LocalNy; ++y) { + const BoutReal theta = TWOPI * (y - mesh->ystart) / ny_physical; + for (int z = 0; z < mesh->LocalNz; ++z) { + const BoutReal phase = theta + 0.2 * x + 0.3 * z; + result.density(x, y, z) = 1.0 + 0.3 * cos(phase); + result.velocity(x, y, z) = + uniform_velocity ? 0.6 : flow_sign * (0.4 + 0.6 * sin(2.0 * phase)); + result.wave(x, y, z) = 0.8 + 0.1 * cos(phase + 0.2); + jacobian(x, y, z) = 1.0 + 0.2 * cos(theta); + dy(x, y, z) = 0.7 + 0.1 * sin(theta); + } + } + } + mesh->getCoordinates()->setJ(jacobian, false); + mesh->getCoordinates()->setDy(dy, false); + return result; + } + + void checkEnergyBalance(bool periodic, bool fix_momentum_flux, BoutReal flow_sign) { + static_cast(mesh)->ix_separatrix = periodic ? mesh->LocalNx : 0; + const auto f = fields(flow_sign); + const auto* coord = mesh->getCoordinates(); + Field3D particle_flow, kinetic_flow; + // The density operator fixes its physical boundary flux to midpoint n*v. + const auto div_particle = + FV::Div_par_mod(f.density, f.velocity, f.wave, particle_flow, true); + const auto div_momentum = + FV::Div_par_fvv(f.density, f.velocity, f.wave, fix_momentum_flux); + const auto heating = FV::Div_par_fvv_heating( + f.density, f.velocity, f.wave, kinetic_flow, fix_momentum_flux); + + BoutReal integrated_change = 0.0, outward = 0.0, heating_magnitude = 0.0; + for (int x = mesh->xstart; x <= mesh->xend; ++x) { + for (int z = 0; z < mesh->LocalNz; ++z) { + outward += kinetic_flow(x, mesh->yend + 1, z) - kinetic_flow(x, mesh->ystart, z); + if (!periodic) { + // Both momentum boundary modes should give the prescribed midpoint + // kinetic energy flux, with heating compensating their different work. + for (int face : {mesh->ystart, mesh->yend + 1}) { + const BoutReal n_mid = + 0.5 * (f.density(x, face - 1, z) + f.density(x, face, z)); + const BoutReal v_mid = + 0.5 * (f.velocity(x, face - 1, z) + f.velocity(x, face, z)); + const BoutReal area = face == mesh->ystart + ? coord->cell_area_ylow()(x, face, z) + : coord->cell_area_yhigh()(x, face - 1, z); + EXPECT_NEAR(kinetic_flow(x, face, z), + 0.5 * n_mid * v_mid * v_mid * v_mid * area, 1e-13); + } + } + for (int y = mesh->ystart; y <= mesh->yend; ++y) { + SCOPED_TRACE(::testing::Message() + << "x=" << x << " y=" << y << " z=" << z << " sign=" << flow_sign); + const BoutReal v = f.velocity(x, y, z); + // d(n*v^2/2)/dt = v*d(n*v)/dt - v^2*dn/dt/2, for unit mass. + const BoutReal kinetic_rhs = + -v * div_momentum(x, y, z) + 0.5 * v * v * div_particle(x, y, z); + const BoutReal volume = coord->cell_volume()(x, y, z); + // kinetic_flow is already area integrated, and is indexed at y-low. + const BoutReal divergence = + (kinetic_flow(x, y + 1, z) - kinetic_flow(x, y, z)) / volume; + EXPECT_NEAR(kinetic_rhs + heating(x, y, z) + divergence, 0.0, 1e-13); + integrated_change += (kinetic_rhs + heating(x, y, z)) * volume; + heating_magnitude += std::abs(heating(x, y, z)) * volume; + } + } + } + EXPECT_GT(heating_magnitude, 1e-4); + EXPECT_NEAR(integrated_change + outward, 0.0, 2e-12); + if (periodic) { + EXPECT_NEAR(outward, 0.0, 1e-13); + } + } +}; + +using HeatingLimiters = ::testing::Types; +TYPED_TEST_SUITE(FVVHeatingTest, HeatingLimiters); + +TYPED_TEST(FVVHeatingTest, PeriodicKineticAndThermalEnergyBalance) { + this->checkEnergyBalance(true, false, 1.0); +} + +TYPED_TEST(FVVHeatingTest, OpenMidpointMomentumEnergyBalance) { + for (BoutReal sign : {1.0, -1.0}) { + this->checkEnergyBalance(false, true, sign); + } +} + +TYPED_TEST(FVVHeatingTest, OpenReconstructedMomentumEnergyBalance) { + for (BoutReal sign : {1.0, -1.0}) { + this->checkEnergyBalance(false, false, sign); + } +} + +TYPED_TEST(FVVHeatingTest, UniformVelocityHasNoPeriodicHeating) { + const auto f = this->fields(1.0, true); + Field3D flow; + const auto heating = + FV::Div_par_fvv_heating(f.density, f.velocity, f.wave, flow); + BOUT_FOR_SERIAL(i, heating.getRegion("RGN_NOBNDRY")) { + EXPECT_DOUBLE_EQ(heating[i], 0.0); + } +} + +} // namespace From 9f2e89856fa62600f28cc1b22b5e1e43ac78967f Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 27 Oct 2026 09:53:18 -0700 Subject: [PATCH 2/5] Div_a_Grad_perp: retain Z diffusion without Y guard cells Implementation changes: - Calculate Z flux independently of the Y flux loop so MYG=0 retains perpendicular diffusion in Z. - Omit Y fluxes and the Y derivative in the Z flux when the global mesh has one physical Y cell, including when Y guards are present. - Preserve Y derivatives for resolved global Y meshes and require guards when those derivatives are needed. Documentation changes: - Explain the independent Z flux and singleton Y behavior in comments. Tests and validation: - Add seven regressions covering Z Fourier modes, independence from Y guard values, variable-coefficient conservation and integration by parts, and retained Y diffusion on a resolved mesh. - Three initial regression cases fail before the fix; all seven new tests and all 27 finite-volume tests pass after it. - Full BOUT serial suite: 2024 passed, four skipped and three failed. - Changed code passes clang-format checks and git diff --check. Notes: - Full-suite failures match the previous build: two default-solver tests require disabled PVODE, and the signal-handler death test fails. --- src/mesh/fv_ops.cxx | 48 +++++++++-------- tests/unit/mesh/test_fv_ops.cxx | 96 +++++++++++++++++++++++++++++++++ 2 files changed, 121 insertions(+), 23 deletions(-) diff --git a/src/mesh/fv_ops.cxx b/src/mesh/fv_ops.cxx index ff854b6cfa..bdbcc63901 100644 --- a/src/mesh/fv_ops.cxx +++ b/src/mesh/fv_ops.cxx @@ -84,7 +84,11 @@ Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f) { } } - // Y and Z fluxes require Y derivatives + // A singleton global Y mesh has no Y derivatives, even if guards exist. + const bool y_derivatives = mesh->GlobalNyNoBoundaries > 1; + ASSERT1(!y_derivatives || mesh->ystart > 0); + + // Y fluxes and the mixed part of the Z flux require Y derivatives. // Fields containing values along the magnetic field @@ -109,7 +113,7 @@ Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f) { } // Y flux - if (mesh->ystart > 0) { + if (y_derivatives) { BOUT_FOR(i, yzresult.getRegion("RGN_NOBNDRY")) { auto ikp = i.zp(); auto ikm = i.zm(); @@ -155,29 +159,27 @@ Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f) { yzresult[i] -= fout / (dy.c[i] * J.c[i]); } - // Z flux - { + } + } - // Coefficient in front of df/dy term - const BoutReal coef = g_23.c[i] / (dy.up[i.yp()] + 2. * dy.c[i] + dy.down[i.ym()]) - / SQ(J.c[i] * Bxy.c[i]); - - const BoutReal fout = - 0.25 * (a_slice.c[i] + a_slice.c[ikp]) - * (J.c[i] * coord->g33()[i] + J.c[ikp] * coord->g33()[ikp]) - * ( // df/dz - (f_slice.c[ikp] - f_slice.c[i]) / dz.c[i] - // - g_yz * df/dy / SQ(J*B) - - coef - * (f_slice.up[i.yp()] + f_slice.up[ikp.yp()] - f_slice.down[i.ym()] - - f_slice.down[ikp.ym()])); - - yzresult[i] += fout / (J.c[i] * dz.c[i]); - yzresult[ikp] -= fout / (J.c[ikp] * dz.c[ikp]); - } + // Z flux exists independently of Y guards or Y derivatives. + BOUT_FOR(i, yzresult.getRegion("RGN_NOBNDRY")) { + const auto ikp = i.zp(); + BoutReal y_correction = 0.0; + if (y_derivatives) { + const BoutReal coef = g_23.c[i] / (dy.up[i.yp()] + 2. * dy.c[i] + dy.down[i.ym()]) + / SQ(J.c[i] * Bxy.c[i]); + y_correction = coef + * (f_slice.up[i.yp()] + f_slice.up[ikp.yp()] - f_slice.down[i.ym()] + - f_slice.down[ikp.ym()]); } - } else { - ASSERT1(mesh->yend == 0); + + const BoutReal fout = 0.25 * (a_slice.c[i] + a_slice.c[ikp]) + * (J.c[i] * coord->g33()[i] + J.c[ikp] * coord->g33()[ikp]) + * ((f_slice.c[ikp] - f_slice.c[i]) / dz.c[i] - y_correction); + + yzresult[i] += fout / (J.c[i] * dz.c[i]); + yzresult[ikp] -= fout / (J.c[ikp] * dz.c[ikp]); } // Check if we need to transform back diff --git a/tests/unit/mesh/test_fv_ops.cxx b/tests/unit/mesh/test_fv_ops.cxx index 92a3340413..1955d24327 100644 --- a/tests/unit/mesh/test_fv_ops.cxx +++ b/tests/unit/mesh/test_fv_ops.cxx @@ -37,6 +37,102 @@ TEST_F(FVOpsTest, Div_a_Grad_perp_FCI_coeffs) { namespace { +template +struct PerpGrid { + static constexpr int ny = NY; + static constexpr int myg = MYG; +}; + +// Exercise a singleton Y with and without guards, and a resolved Y mesh. +template +class PerpDiffusionTest : public FakeMeshFixture_tmpl<5, Grid::ny, 16, false, Grid::myg> { +public: + void SetUp() override { + auto* coord = mesh->getCoordinates(); + // Nonzero g23 ensures the tests also exercise the mixed Y-Z terms. + const BoutReal cross = 0.3; + const BoutReal determinant = 1.0 - cross * cross; + coord->setMetricTensorJB( + ContravariantMetricTensor(1.0, 1.0, 1.0, 0.0, 0.0, cross), + CovariantMetricTensor(1.0, 1.0 / determinant, 1.0 / determinant, 0.0, 0.0, + -cross / determinant), + Coordinates::FieldMetric{1.0 / sqrt(determinant)}, Coordinates::FieldMetric{1.0}); + coord->setDy(Coordinates::FieldMetric{0.8}, false); + coord->setDz(Coordinates::FieldMetric{0.7}, false); + } + + Field3D zMode(bool vary_y_guards = false) { + Field3D f{0.0}; + BOUT_FOR_SERIAL(i, f.getRegion("RGN_ALL")) { + const BoutReal amplitude = + vary_y_guards ? 1.0 + 0.25 * SQ(i.y() - mesh->ystart) : 1.0; + f[i] = amplitude * cos(TWOPI * i.z() / mesh->LocalNz); + } + return f; + } +}; + +using PerpGrids = ::testing::Types, PerpGrid<3, 1>, PerpGrid<6, 1>>; +TYPED_TEST_SUITE(PerpDiffusionTest, PerpGrids); + +TYPED_TEST(PerpDiffusionTest, ZFourierMode) { + const auto f = this->zMode(mesh->GlobalNyNoBoundaries == 1); + const BoutReal diffusivity = 2.5; + const auto result = FV::Div_a_Grad_perp(Field3D{diffusivity}, f); + const auto* coord = mesh->getCoordinates(); + BOUT_FOR_SERIAL(i, result.getRegion("RGN_NOBNDRY")) { + // The centered periodic Z Laplacian has this exact discrete eigenvalue. + const BoutReal eigenvalue = -4.0 * diffusivity * coord->g33()[i] + * SQ(sin(PI / mesh->LocalNz)) / SQ(coord->dz()[i]); + EXPECT_NEAR(result[i], eigenvalue * f[i], 1e-13); + } +} + +TYPED_TEST(PerpDiffusionTest, VariableCoefficientConservesAndDissipates) { + const auto f = this->zMode(); + Field3D a{0.0}; + Coordinates::FieldMetric jacobian{1.0}; + BOUT_FOR_SERIAL(i, a.getRegion("RGN_ALL")) { + a[i] = 2.0 + 0.2 * i.x() + 0.4 * cos(TWOPI * i.z() / mesh->LocalNz); + jacobian(i.x(), i.y(), i.z()) = 1.0 + 0.1 * i.x(); + } + auto* coord = mesh->getCoordinates(); + coord->setJ(jacobian, false); + const auto result = FV::Div_a_Grad_perp(a, f); + + BoutReal integral = 0.0, dissipation = 0.0, expected_dissipation = 0.0; + BOUT_FOR_SERIAL(i, result.getRegion("RGN_NOBNDRY")) { + const BoutReal volume = coord->cell_volume()[i]; + integral += result[i] * volume; + dissipation += f[i] * result[i] * volume; + const BoutReal gradient = (f[i.zp()] - f[i]) / coord->dz()[i]; + expected_dissipation -= + 0.5 * (a[i] + a[i.zp()]) * coord->g33()[i] * SQ(gradient) * volume; + } + EXPECT_NEAR(integral, 0.0, 1e-12); + EXPECT_LT(dissipation, -1.0); + EXPECT_NEAR(dissipation, expected_dissipation, 1e-12); +} + +using PerpResolvedYTest = PerpDiffusionTest>; + +TEST_F(PerpResolvedYTest, YFourierMode) { + Field3D f{0.0}; + BOUT_FOR_SERIAL(i, f.getRegion("RGN_ALL")) { + f[i] = cos(TWOPI * (i.y() - mesh->ystart) / mesh->GlobalNyNoBoundaries); + } + const auto result = FV::Div_a_Grad_perp(Field3D{1.0}, f); + const auto* coord = mesh->getCoordinates(); + BOUT_FOR_SERIAL(i, result.getRegion("RGN_NOBNDRY")) { + // The perpendicular Y coefficient in this nonorthogonal metric is nonzero. + const BoutReal perp_yy = + -coord->g23()[i] * coord->g_23()[i] / SQ(coord->J()[i] * coord->Bxy()[i]); + const BoutReal eigenvalue = + -4.0 * perp_yy * SQ(sin(PI / mesh->GlobalNyNoBoundaries)) / SQ(coord->dy()[i]); + EXPECT_NEAR(result[i], eigenvalue * f[i], 1e-13); + } +} + // Two Y guards let the parallel fluxes reconstruct their neighboring guard cell. // Keep the existing one-guard fixture as the default for other unit tests. template From c6916cd9e69c0446415d3d7c3cf74d1a2be71da7 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 30 Oct 2026 11:58:09 -0700 Subject: [PATCH 3/5] FV::Div_a_Grad_perp unit tests: Fix MYG=0, increase coverage The fixture now handles cases without Y guards. Added a test case that exercises the mixed Y-Z derivatives. --- tests/unit/fake_mesh_fixture.hxx | 5 +++-- tests/unit/mesh/test_fv_ops.cxx | 32 ++++++++++++++++++++++++++++++++ 2 files changed, 35 insertions(+), 2 deletions(-) diff --git a/tests/unit/fake_mesh_fixture.hxx b/tests/unit/fake_mesh_fixture.hxx index a2b8417418..f029edd91d 100644 --- a/tests/unit/fake_mesh_fixture.hxx +++ b/tests/unit/fake_mesh_fixture.hxx @@ -62,7 +62,8 @@ public: test_coords->setD1_dy(0.2); test_coords->setD1_dz(0.0); - if (bout::build::use_metric_3d) { + // Parallel slices are only allocated when Y guard cells exist. + if (bout::build::use_metric_3d && MYG > 0) { bout::FieldMetric mutable_Bxy = test_coords->Bxy(); mutable_Bxy.splitParallelSlices(); mutable_Bxy.yup() = test_coords->Bxy(); @@ -104,7 +105,7 @@ public: test_coords_staggered->setD1_dy(0.2); test_coords_staggered->setD1_dz(0.0); - if (bout::build::use_metric_3d) { + if (bout::build::use_metric_3d && MYG > 0) { bout::FieldMetric mutable_Bxy = test_coords_staggered->Bxy(); mutable_Bxy.splitParallelSlices(); test_coords_staggered->setBxy(mutable_Bxy); diff --git a/tests/unit/mesh/test_fv_ops.cxx b/tests/unit/mesh/test_fv_ops.cxx index 1955d24327..f7ea4a6e33 100644 --- a/tests/unit/mesh/test_fv_ops.cxx +++ b/tests/unit/mesh/test_fv_ops.cxx @@ -133,6 +133,38 @@ TEST_F(PerpResolvedYTest, YFourierMode) { } } +TEST_F(PerpResolvedYTest, MixedYZFourierMode) { + const BoutReal ky = TWOPI / mesh->GlobalNyNoBoundaries; + const BoutReal diffusivity = 2.5; + const auto* coord = mesh->getCoordinates(); + // Reversing kz reverses the mixed contribution, leaving the pure terms unchanged. + for (BoutReal sign : {1.0, -1.0}) { + SCOPED_TRACE(::testing::Message() << "Z mode sign=" << sign); + const BoutReal kz = sign * TWOPI / mesh->LocalNz; + Field3D f{0.0}; + BOUT_FOR_SERIAL(i, f.getRegion("RGN_ALL")) { + // Include periodic Y guards as well as the physical cells. + f[i] = cos(ky * (i.y() - mesh->ystart) + kz * i.z()); + } + const auto result = FV::Div_a_Grad_perp(Field3D{diffusivity}, f); + BOUT_FOR_SERIAL(i, result.getRegion("RGN_NOBNDRY")) { + const BoutReal inverse_JB_squared = 1.0 / SQ(coord->J()[i] * coord->Bxy()[i]); + const BoutReal perp_yy = -coord->g23()[i] * coord->g_23()[i] * inverse_JB_squared; + const BoutReal perp_yz = + 0.5 + * (coord->g23()[i] - coord->g33()[i] * coord->g_23()[i] * inverse_JB_squared); + // The pure terms use three-point second derivatives; the mixed term + // uses centered first derivatives in both Y and Z. + const BoutReal eigenvalue = + -diffusivity + * (4.0 * perp_yy * SQ(sin(0.5 * ky)) / SQ(coord->dy()[i]) + + 4.0 * coord->g33()[i] * SQ(sin(0.5 * kz)) / SQ(coord->dz()[i]) + + 2.0 * perp_yz * sin(ky) * sin(kz) / (coord->dy()[i] * coord->dz()[i])); + EXPECT_NEAR(result[i], eigenvalue * f[i], 1e-13); + } + } +} + // Two Y guards let the parallel fluxes reconstruct their neighboring guard cell. // Keep the existing one-guard fixture as the default for other unit tests. template From 78c0962426295ab0aba66a973bff9648e590b16d Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sun, 8 Nov 2026 09:59:12 -0800 Subject: [PATCH 4/5] Apply batched suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> Co-authored-by: David Bold --- src/mesh/fv_ops.cxx | 2 +- tests/unit/fake_mesh.hxx | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/mesh/fv_ops.cxx b/src/mesh/fv_ops.cxx index bdbcc63901..52553f4720 100644 --- a/src/mesh/fv_ops.cxx +++ b/src/mesh/fv_ops.cxx @@ -85,7 +85,7 @@ Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f) { } // A singleton global Y mesh has no Y derivatives, even if guards exist. - const bool y_derivatives = mesh->GlobalNyNoBoundaries > 1; + const bool y_derivatives = mesh->GlobalNy > 1; ASSERT1(!y_derivatives || mesh->ystart > 0); // Y fluxes and the mixed part of the Z flux require Y derivatives. diff --git a/tests/unit/fake_mesh.hxx b/tests/unit/fake_mesh.hxx index 81e116f4aa..062bd5fea5 100644 --- a/tests/unit/fake_mesh.hxx +++ b/tests/unit/fake_mesh.hxx @@ -48,7 +48,7 @@ public: GlobalNy = ny; GlobalNz = nz; GlobalNxNoBoundaries = nx - 2; - GlobalNyNoBoundaries = ny - 2 * myg; + GlobalNyNoBoundaries = ny - (2 * myg); GlobalNzNoBoundaries = nz; LocalNx = nx; LocalNy = ny; @@ -59,13 +59,13 @@ public: // These bits only for ADIOS2, also boring due to single process MapCountX = nx - 2; - MapCountY = ny - 2 * myg; + MapCountY = ny - (2 * myg); MapCountZ = nz; MapGlobalX = nx; MapGlobalY = ny; MapGlobalZ = nz; MapLocalX = nx - 2; - MapLocalY = ny - 2 * myg; + MapLocalY = ny - (2 * myg); MapLocalZ = nz; // Small "inner" region From a086cd67bf7df69aaf42d15bc647f7a2ea8a86dd Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Sun, 8 Nov 2026 15:57:38 -0800 Subject: [PATCH 5/5] FV ops unit test: Include Y derivs when GlobalNy>1 Now consistent with changed FV operator. Many small Clang tidy fixes. --- tests/unit/mesh/test_fv_ops.cxx | 34 +++++++++++++++++++-------------- 1 file changed, 20 insertions(+), 14 deletions(-) diff --git a/tests/unit/mesh/test_fv_ops.cxx b/tests/unit/mesh/test_fv_ops.cxx index f7ea4a6e33..4f627280e4 100644 --- a/tests/unit/mesh/test_fv_ops.cxx +++ b/tests/unit/mesh/test_fv_ops.cxx @@ -1,13 +1,19 @@ #include "gtest/gtest.h" #include +#include "bout/bout_types.hxx" +#include "bout/boutexception.hxx" #include "bout/build_defines.hxx" #include "bout/constants.hxx" #include "bout/coordinates.hxx" +#include "bout/field.hxx" +#include "bout/field3d.hxx" #include "bout/fv_ops.hxx" #include "bout/fv_ops_impl.hxx" +#include "bout/globals.hxx" #include "bout/mesh.hxx" -#include "bout/output.hxx" +#include "bout/region.hxx" +#include "bout/utils.hxx" #include "fake_mesh.hxx" #include "fake_mesh_fixture.hxx" @@ -30,7 +36,7 @@ TEST_F(FVOpsTest, Div_a_Grad_perp_FCI_coeffs) { {"dagp_fv_volume", 1.0}})); auto op = FV::Div_a_Grad_perp_FCI::create(mesh, 1); - Field3D a{1.0}, f{1.0}; + const Field3D a{1.0}, f{1.0}; (*op)(a, f, true); (*op)(a, f, false); } @@ -51,7 +57,7 @@ class PerpDiffusionTest : public FakeMeshFixture_tmpl<5, Grid::ny, 16, false, Gr auto* coord = mesh->getCoordinates(); // Nonzero g23 ensures the tests also exercise the mixed Y-Z terms. const BoutReal cross = 0.3; - const BoutReal determinant = 1.0 - cross * cross; + const BoutReal determinant = 1.0 - (cross * cross); coord->setMetricTensorJB( ContravariantMetricTensor(1.0, 1.0, 1.0, 0.0, 0.0, cross), CovariantMetricTensor(1.0, 1.0 / determinant, 1.0 / determinant, 0.0, 0.0, @@ -65,7 +71,7 @@ class PerpDiffusionTest : public FakeMeshFixture_tmpl<5, Grid::ny, 16, false, Gr Field3D f{0.0}; BOUT_FOR_SERIAL(i, f.getRegion("RGN_ALL")) { const BoutReal amplitude = - vary_y_guards ? 1.0 + 0.25 * SQ(i.y() - mesh->ystart) : 1.0; + vary_y_guards ? 1.0 + (0.25 * SQ(i.y() - mesh->ystart)) : 1.0; f[i] = amplitude * cos(TWOPI * i.z() / mesh->LocalNz); } return f; @@ -76,7 +82,7 @@ using PerpGrids = ::testing::Types, PerpGrid<3, 1>, PerpGrid<6, 1 TYPED_TEST_SUITE(PerpDiffusionTest, PerpGrids); TYPED_TEST(PerpDiffusionTest, ZFourierMode) { - const auto f = this->zMode(mesh->GlobalNyNoBoundaries == 1); + const auto f = this->zMode(mesh->GlobalNy == 1); const BoutReal diffusivity = 2.5; const auto result = FV::Div_a_Grad_perp(Field3D{diffusivity}, f); const auto* coord = mesh->getCoordinates(); @@ -93,8 +99,8 @@ TYPED_TEST(PerpDiffusionTest, VariableCoefficientConservesAndDissipates) { Field3D a{0.0}; Coordinates::FieldMetric jacobian{1.0}; BOUT_FOR_SERIAL(i, a.getRegion("RGN_ALL")) { - a[i] = 2.0 + 0.2 * i.x() + 0.4 * cos(TWOPI * i.z() / mesh->LocalNz); - jacobian(i.x(), i.y(), i.z()) = 1.0 + 0.1 * i.x(); + a[i] = 2.0 + (0.2 * i.x()) + (0.4 * cos(TWOPI * i.z() / mesh->LocalNz)); + jacobian(i.x(), i.y(), i.z()) = 1.0 + (0.1 * i.x()); } auto* coord = mesh->getCoordinates(); coord->setJ(jacobian, false); @@ -138,13 +144,13 @@ TEST_F(PerpResolvedYTest, MixedYZFourierMode) { const BoutReal diffusivity = 2.5; const auto* coord = mesh->getCoordinates(); // Reversing kz reverses the mixed contribution, leaving the pure terms unchanged. - for (BoutReal sign : {1.0, -1.0}) { + for (const BoutReal sign : {1.0, -1.0}) { SCOPED_TRACE(::testing::Message() << "Z mode sign=" << sign); const BoutReal kz = sign * TWOPI / mesh->LocalNz; Field3D f{0.0}; BOUT_FOR_SERIAL(i, f.getRegion("RGN_ALL")) { // Include periodic Y guards as well as the physical cells. - f[i] = cos(ky * (i.y() - mesh->ystart) + kz * i.z()); + f[i] = cos(ky * (i.y() - mesh->ystart) + (kz * i.z())); } const auto result = FV::Div_a_Grad_perp(Field3D{diffusivity}, f); BOUT_FOR_SERIAL(i, result.getRegion("RGN_NOBNDRY")) { @@ -182,13 +188,13 @@ class FVVHeatingTest : public FakeMeshFixture_tmpl<5, 20, 3, false, 2> { for (int y = 0; y < mesh->LocalNy; ++y) { const BoutReal theta = TWOPI * (y - mesh->ystart) / ny_physical; for (int z = 0; z < mesh->LocalNz; ++z) { - const BoutReal phase = theta + 0.2 * x + 0.3 * z; - result.density(x, y, z) = 1.0 + 0.3 * cos(phase); + const BoutReal phase = theta + (0.2 * x) + (0.3 * z); + result.density(x, y, z) = 1.0 + (0.3 * cos(phase)); result.velocity(x, y, z) = uniform_velocity ? 0.6 : flow_sign * (0.4 + 0.6 * sin(2.0 * phase)); - result.wave(x, y, z) = 0.8 + 0.1 * cos(phase + 0.2); - jacobian(x, y, z) = 1.0 + 0.2 * cos(theta); - dy(x, y, z) = 0.7 + 0.1 * sin(theta); + result.wave(x, y, z) = 0.8 + (0.1 * cos(phase + 0.2)); + jacobian(x, y, z) = 1.0 + (0.2 * cos(theta)); + dy(x, y, z) = 0.7 + (0.1 * sin(theta)); } } }