From 353e709d6f9958faeda4e381efe6d52097bcb571 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 23 Feb 2026 17:22:15 +0000 Subject: [PATCH 01/12] Disable LaplaceXYHypre and test if using 3D metrics --- src/invert/laplacexy/impls/hypre/laplacexy-hypre.cxx | 2 +- src/invert/laplacexy/impls/hypre/laplacexy-hypre.hxx | 6 ++++++ tests/integrated/test-laplacexy2-hypre/CMakeLists.txt | 1 + 3 files changed, 8 insertions(+), 1 deletion(-) diff --git a/src/invert/laplacexy/impls/hypre/laplacexy-hypre.cxx b/src/invert/laplacexy/impls/hypre/laplacexy-hypre.cxx index 9323513f21..439d07a4e5 100644 --- a/src/invert/laplacexy/impls/hypre/laplacexy-hypre.cxx +++ b/src/invert/laplacexy/impls/hypre/laplacexy-hypre.cxx @@ -1,6 +1,6 @@ #include "bout/build_defines.hxx" -#if BOUT_HAS_HYPRE +#if BOUT_HAS_HYPRE and not BOUT_USE_METRIC_3D #include "laplacexy-hypre.hxx" diff --git a/src/invert/laplacexy/impls/hypre/laplacexy-hypre.hxx b/src/invert/laplacexy/impls/hypre/laplacexy-hypre.hxx index f415f9ebf4..f9b1d90fb3 100644 --- a/src/invert/laplacexy/impls/hypre/laplacexy-hypre.hxx +++ b/src/invert/laplacexy/impls/hypre/laplacexy-hypre.hxx @@ -43,6 +43,12 @@ RegisterUnavailableLaplaceXY registerlaplacexyhypre("hypre", "BOUT++ was not configured with HYPRE"); } +#elif BOUT_USE_METRIC_3D +namespace { +RegisterUnavailableLaplaceXY + registerlaplacexyhypre("hypre", "BOUT++ was configured with 3D metrics"); +} + #else // BOUT_HAS_HYPRE class Mesh; diff --git a/tests/integrated/test-laplacexy2-hypre/CMakeLists.txt b/tests/integrated/test-laplacexy2-hypre/CMakeLists.txt index e56829b243..4f16f78713 100644 --- a/tests/integrated/test-laplacexy2-hypre/CMakeLists.txt +++ b/tests/integrated/test-laplacexy2-hypre/CMakeLists.txt @@ -2,5 +2,6 @@ bout_add_integrated_test( test-laplacexy2-hypre SOURCES test-laplacexy.cxx REQUIRES BOUT_HAS_HYPRE + CONFLICTS BOUT_USE_METRIC_3D USE_RUNTEST USE_DATA_BOUT_INP ) From f887aed94e50bcef5b7cf9bc4568f62435b581c0 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 23 Feb 2026 17:24:51 +0000 Subject: [PATCH 02/12] CI: remove (mostly) redundant test This test has significant overlap with other tests, it doesn't really give us any useful information --- .github/workflows/tests.yml | 16 ---------------- 1 file changed, 16 deletions(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index f44ebb87ec..0abc01ef5f 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -121,22 +121,6 @@ jobs: omp_num_threads: 2 on_cron: false - - name: "CMake, new PETSc" - os: ubuntu-latest - cmake_options: "-DBUILD_SHARED_LIBS=ON - -DBOUT_ENABLE_METRIC_3D=ON - -DBOUT_ENABLE_OPENMP=ON - -DBOUT_USE_PETSC=ON - -DBOUT_USE_SLEPC=ON - -DBOUT_USE_SUNDIALS=ON - -DBOUT_USE_HYPRE=OFF - -DBOUT_ENABLE_PYTHON=ON - -DSUNDIALS_ROOT=/home/runner/local - -DPETSC_DIR=/home/runner/local/petsc - -DSLEPC_DIR=/home/runner/local/slepc" - build_petsc: -petsc - on_cron: false - exclude: - is_cron: true config: From 4c96e9b9639e236def795f1efd17289c63e52828 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 23 Feb 2026 17:25:41 +0000 Subject: [PATCH 03/12] CI: rename build to be more accurate --- .github/workflows/tests.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 0abc01ef5f..e45e082232 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -40,7 +40,7 @@ jobs: is_cron: - ${{ github.event_name == 'cron' }} config: - - name: "CMake, PETSc unreleased, ADIOS2" + - name: "PETSc unreleased, ADIOS2, 3D metrics" os: ubuntu-24.04 cmake_options: "-DBUILD_SHARED_LIBS=ON -DBOUT_ENABLE_METRIC_3D=ON From 18ed91d855ebd7ac7d230b59c5ad44aeb818d426 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 23 Feb 2026 18:45:05 +0000 Subject: [PATCH 04/12] Mark unavailable Laplacian solvers that don't work with 3D metrics --- .../laplace/impls/naulin/naulin_laplace.cxx | 6 ++++++ .../laplace/impls/naulin/naulin_laplace.hxx | 15 ++++++++++++++- src/invert/laplace/impls/pcr/pcr.cxx | 6 ++++++ src/invert/laplace/impls/pcr/pcr.hxx | 17 +++++++++++++++-- .../laplace/impls/pcr_thomas/pcr_thomas.cxx | 6 ++++++ .../laplace/impls/pcr_thomas/pcr_thomas.hxx | 15 ++++++++++++++- .../laplace/impls/serial_tri/serial_tri.cxx | 6 ++++++ .../laplace/impls/serial_tri/serial_tri.hxx | 15 ++++++++++++++- src/invert/laplace/impls/spt/spt.cxx | 6 ++++++ src/invert/laplace/impls/spt/spt.hxx | 16 ++++++++++++++-- 10 files changed, 101 insertions(+), 7 deletions(-) diff --git a/src/invert/laplace/impls/naulin/naulin_laplace.cxx b/src/invert/laplace/impls/naulin/naulin_laplace.cxx index f61dc93908..1c90976ce6 100644 --- a/src/invert/laplace/impls/naulin/naulin_laplace.cxx +++ b/src/invert/laplace/impls/naulin/naulin_laplace.cxx @@ -138,6 +138,10 @@ */ // clang-format on +#include "bout/build_defines.hxx" + +#if not BOUT_USE_METRIC_3D + #include #include #include @@ -426,3 +430,5 @@ void LaplaceNaulin::outputVars(Options& output_options, output_options[fmt::format("{}_mean_underrelax_counts", getPerformanceName())] .assignRepeat(naulinsolver_mean_underrelax_counts, time_dimension); } + +#endif diff --git a/src/invert/laplace/impls/naulin/naulin_laplace.hxx b/src/invert/laplace/impls/naulin/naulin_laplace.hxx index 1998a046d7..c4e7d41d9e 100644 --- a/src/invert/laplace/impls/naulin/naulin_laplace.hxx +++ b/src/invert/laplace/impls/naulin/naulin_laplace.hxx @@ -28,11 +28,22 @@ class LaplaceNaulin; #ifndef BOUT_LAP_NAULIN_H #define BOUT_LAP_NAULIN_H +#include #include + +#if BOUT_USE_METRIC_3D + +namespace { +const RegisterUnavailableLaplace + registerlaplacenaulin(LAPLACE_NAULIN, "BOUT++ was configured with 3D metrics"); +} + +#else + #include namespace { -RegisterLaplace registerlaplacenaulin(LAPLACE_NAULIN); +const RegisterLaplace registerlaplacenaulin(LAPLACE_NAULIN); } /// Solves the 2D Laplacian equation @@ -198,4 +209,6 @@ private: void copy_x_boundaries(Field3D& x, const Field3D& x0, Mesh* mesh); }; +#endif // BOUT_USE_METRIC_3D + #endif // BOUT_LAP_NAULIN_H diff --git a/src/invert/laplace/impls/pcr/pcr.cxx b/src/invert/laplace/impls/pcr/pcr.cxx index a33bd7eef2..591edb5080 100644 --- a/src/invert/laplace/impls/pcr/pcr.cxx +++ b/src/invert/laplace/impls/pcr/pcr.cxx @@ -36,6 +36,10 @@ * **************************************************************************/ +#include "bout/build_defines.hxx" + +#if not BOUT_USE_METRIC_3D + #include "pcr.hxx" #include "bout/globals.hxx" @@ -1108,3 +1112,5 @@ void LaplacePCR ::verify_solution(const Matrix& a_ver, output.write("max abs error {}\n", max_error); output.write("max abs error location {} {}\n", max_loc_x, max_loc_z); } + +#endif // BOUT_USE_METRIC_3D diff --git a/src/invert/laplace/impls/pcr/pcr.hxx b/src/invert/laplace/impls/pcr/pcr.hxx index ec4637f56c..2b90e17466 100644 --- a/src/invert/laplace/impls/pcr/pcr.hxx +++ b/src/invert/laplace/impls/pcr/pcr.hxx @@ -29,13 +29,24 @@ class LaplacePCR; #ifndef BOUT_PCR_H #define BOUT_PCR_H -#include +#include #include + +#if BOUT_USE_METRIC_3D + +namespace { +const RegisterUnavailableLaplace + registerlaplacepcr(LAPLACE_PCR, "BOUT++ was configured with 3D metrics"); +} + +#else + +#include #include #include namespace { -RegisterLaplace registerlaplacepcr(LAPLACE_PCR); +const RegisterLaplace registerlaplacepcr(LAPLACE_PCR); } class LaplacePCR : public Laplacian { @@ -175,4 +186,6 @@ private: bool dst{false}; }; +#endif // BOUT_USE_METRIC_3D + #endif // BOUT_PCR_H diff --git a/src/invert/laplace/impls/pcr_thomas/pcr_thomas.cxx b/src/invert/laplace/impls/pcr_thomas/pcr_thomas.cxx index d471775990..1b8c66dfdf 100644 --- a/src/invert/laplace/impls/pcr_thomas/pcr_thomas.cxx +++ b/src/invert/laplace/impls/pcr_thomas/pcr_thomas.cxx @@ -36,6 +36,10 @@ * **************************************************************************/ +#include "bout/build_defines.hxx" + +#if not BOUT_USE_METRIC_3D + #include "pcr_thomas.hxx" #include "bout/globals.hxx" @@ -1190,3 +1194,5 @@ void LaplacePCR_THOMAS ::verify_solution(const Matrix& a_ver, } output.write("max abs error {}\n", max_error); } + +#endif // BOUT_USE_METRIC_3D diff --git a/src/invert/laplace/impls/pcr_thomas/pcr_thomas.hxx b/src/invert/laplace/impls/pcr_thomas/pcr_thomas.hxx index e12a647789..1b5c6b97b9 100644 --- a/src/invert/laplace/impls/pcr_thomas/pcr_thomas.hxx +++ b/src/invert/laplace/impls/pcr_thomas/pcr_thomas.hxx @@ -29,8 +29,19 @@ class LaplacePCR_THOMAS; #ifndef BOUT_PCR_THOMAS_H #define BOUT_PCR_THOMAS_H +#include "bout/build_defines.hxx" +#include "bout/invert_laplace.hxx" + +#if BOUT_USE_METRIC_3D + +namespace { +const RegisterUnavailableLaplace + registerlaplacepcrthomas(LAPLACE_PCR_THOMAS, "BOUT++ was configured with 3D metrics"); +} + +#else + #include -#include #include #include @@ -178,4 +189,6 @@ private: bool dst{false}; }; +#endif // BOUT_USE_METRIC_3D + #endif // BOUT_PCR_THOMAS_H diff --git a/src/invert/laplace/impls/serial_tri/serial_tri.cxx b/src/invert/laplace/impls/serial_tri/serial_tri.cxx index a14e0e4a26..d051ce0c1e 100644 --- a/src/invert/laplace/impls/serial_tri/serial_tri.cxx +++ b/src/invert/laplace/impls/serial_tri/serial_tri.cxx @@ -24,6 +24,10 @@ * **************************************************************************/ +#include "bout/build_defines.hxx" + +#if not BOUT_USE_METRIC_3D + #include "serial_tri.hxx" #include "bout/globals.hxx" @@ -247,3 +251,5 @@ FieldPerp LaplaceSerialTri::solve(const FieldPerp& b, const FieldPerp& x0) { return x; // Result of the inversion } + +#endif // BOUT_USE_METRIC_3D diff --git a/src/invert/laplace/impls/serial_tri/serial_tri.hxx b/src/invert/laplace/impls/serial_tri/serial_tri.hxx index 5b0419fa27..4aed777b7c 100644 --- a/src/invert/laplace/impls/serial_tri/serial_tri.hxx +++ b/src/invert/laplace/impls/serial_tri/serial_tri.hxx @@ -29,8 +29,19 @@ class LaplaceSerialTri; #ifndef BOUT_SERIAL_TRI_H #define BOUT_SERIAL_TRI_H +#include "bout/build_defines.hxx" +#include "bout/invert_laplace.hxx" + +#if BOUT_USE_METRIC_3D + +namespace { +const RegisterUnavailableLaplace + registerlaplacetri(LAPLACE_TRI, "BOUT++ was configured with 3D metrics"); +} + +#else + #include -#include #include namespace { @@ -80,4 +91,6 @@ private: Field2D A, C, D; }; +#endif // BOUT_USE_METRIC_3D + #endif // BOUT_SERIAL_TRI_H diff --git a/src/invert/laplace/impls/spt/spt.cxx b/src/invert/laplace/impls/spt/spt.cxx index cd24ee1acf..34a1a64cae 100644 --- a/src/invert/laplace/impls/spt/spt.cxx +++ b/src/invert/laplace/impls/spt/spt.cxx @@ -31,6 +31,10 @@ * */ +#include "bout/build_defines.hxx" + +#if not BOUT_USE_METRIC_3D + #include #include #include @@ -553,3 +557,5 @@ void LaplaceSPT::SPT_data::allocate(int mm, int nx) { buffer.reallocate(4 * mm); } + +#endif // BOUT_USE_METRIC_3D diff --git a/src/invert/laplace/impls/spt/spt.hxx b/src/invert/laplace/impls/spt/spt.hxx index a9d5b2583f..a3b1acb7c3 100644 --- a/src/invert/laplace/impls/spt/spt.hxx +++ b/src/invert/laplace/impls/spt/spt.hxx @@ -39,10 +39,20 @@ class LaplaceSPT; #ifndef BOUT_SPT_H +#include "bout/build_defines.hxx" +#include "bout/invert_laplace.hxx" + +#if BOUT_USE_METRIC_3D + +namespace { +const RegisterUnavailableLaplace + registerlaplacespt(LAPLACE_SPT, "BOUT++ was configured with 3D metrics"); +} + +#else #define BOUT_SPT_H #include -#include #include #include #include @@ -153,7 +163,9 @@ private: namespace { // Note: After class definition so compiler knows that // registered class is derived from Laplacian -RegisterLaplace registerlaplacespt(LAPLACE_SPT); +const RegisterLaplace registerlaplacespt(LAPLACE_SPT); } // namespace +#endif // BOUT_USE_METRIC_3D + #endif // BOUT_SPT_H From daf63738b55d315bb221bf20bc4d8e28b7a7c10c Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 23 Feb 2026 18:48:08 +0000 Subject: [PATCH 05/12] Get PETSc Laplacian test working with 3D metrics - Remove comparison to default solver - The default is cyclic which doesn't work with 3D metrics - Also we don't actually do any comparison, so it doesn't help us test anything - Use one of PETSc's preconditioners instead of one of our solvers - Prefer MUMPS but fall back to LU, which should always be available --- .../test-petsc_laplace/CMakeLists.txt | 2 - .../test-petsc_laplace/data/BOUT.inp | 9 --- tests/integrated/test-petsc_laplace/runtest | 2 - .../test-petsc_laplace/test_petsc_laplace.cxx | 80 +++++++++++-------- 4 files changed, 46 insertions(+), 47 deletions(-) diff --git a/tests/integrated/test-petsc_laplace/CMakeLists.txt b/tests/integrated/test-petsc_laplace/CMakeLists.txt index 1d2e7896f0..977b33afe2 100644 --- a/tests/integrated/test-petsc_laplace/CMakeLists.txt +++ b/tests/integrated/test-petsc_laplace/CMakeLists.txt @@ -2,8 +2,6 @@ bout_add_integrated_test( test-petsc-laplace SOURCES test_petsc_laplace.cxx REQUIRES BOUT_HAS_PETSC - CONFLICTS - BOUT_USE_METRIC_3D # default preconditioner uses 'cyclic' Laplace solver which is not available with 3d metrics USE_RUNTEST USE_DATA_BOUT_INP PROCESSORS 4 ) diff --git a/tests/integrated/test-petsc_laplace/data/BOUT.inp b/tests/integrated/test-petsc_laplace/data/BOUT.inp index e7c285b54c..15e54ff753 100644 --- a/tests/integrated/test-petsc_laplace/data/BOUT.inp +++ b/tests/integrated/test-petsc_laplace/data/BOUT.inp @@ -67,15 +67,6 @@ outer_boundary_flags = 32 # Identity in boundary ############################################# -[SPT] -#type=spt -all_terms = true -nonuniform = true -#flags=15 -include_yguards = false - -#maxits=10000 - [laplace] all_terms = true nonuniform = true diff --git a/tests/integrated/test-petsc_laplace/runtest b/tests/integrated/test-petsc_laplace/runtest index 83e1006338..5496987128 100755 --- a/tests/integrated/test-petsc_laplace/runtest +++ b/tests/integrated/test-petsc_laplace/runtest @@ -16,11 +16,9 @@ vars = [ ("max_error1", 2.0e-4), ("max_error2", 2.0e-4), ("max_error3", 2.0e-4), - ("max_error4", 1.0e-5), ("max_error5", 2.0e-4), ("max_error6", 2.0e-5), ("max_error7", 2.0e-4), - ("max_error8", 2.0e-5), ] # tol = 1e-4 # Absolute (?) tolerance diff --git a/tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx b/tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx index c42c55d8d6..e869cd1b1d 100644 --- a/tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx +++ b/tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx @@ -26,6 +26,7 @@ #include "bout/bout.hxx" // NOLINT #include "bout/bout_types.hxx" #include "bout/boutexception.hxx" +#include "bout/build_config.hxx" #include "bout/constants.hxx" #include "bout/difops.hxx" #include "bout/field2d.hxx" @@ -39,6 +40,7 @@ #include "fmt/core.h" #include +#include #include #include @@ -112,10 +114,37 @@ int main(int argc, char** argv) { BoutInitialise(argc, argv); { - Options* options = Options::getRoot()->getSection("petsc2nd"); - auto invert = Laplacian::create(options); - options = Options::getRoot()->getSection("petsc4th"); - auto invert_4th = Laplacian::create(options); + // Not be used for 3D metrics + Options::root()["laplace"].setConditionallyUsed(); + + // For 3D metrics, we need to use one of PETSc's preconditioners + // and mumps seems to be the best? Or at least HYPRE doesn't work + // on this test with more than one processor. +#ifdef PETSC_HAVE_MUMPS + constexpr auto petsc_has_mumps = true; +#else + constexpr auto petsc_has_mumps = false; +#endif + + // Preconditioner to use for 3D metrics + constexpr auto petsc_pc = petsc_has_mumps ? "mumps" : "lu"; + + auto& options_2nd = Options::root()["petsc2nd"]; + if constexpr (bout::build::use_metric_3d) { + options_2nd["pctype"] = petsc_pc; + options_2nd["rightprec"].setConditionallyUsed(); + options_2nd["precon"].setConditionallyUsed(); + } + + auto invert = Laplacian::create(&options_2nd); + + auto& options_4th = Options::root()["petsc4th"]; + if constexpr (bout::build::use_metric_3d) { + options_4th["pctype"] = petsc_pc; + options_4th["rightprec"].setConditionallyUsed(); + options_4th["precon"].setConditionallyUsed(); + } + auto invert_4th = Laplacian::create(&options_4th); Options dump; @@ -136,35 +165,27 @@ int main(int argc, char** argv) { const Field3D b_1 = forward_laplace(f_1, a_1, c_1, d_1); - int test_num = 0; - check_laplace(++test_num, "PETSc 2nd order", *invert, INVERT_AC_GRAD, INVERT_AC_GRAD, - a_1, c_1, d_1, b_1, f_1, mesh->ystart, dump); + check_laplace(1, "PETSc 2nd order", *invert, INVERT_AC_GRAD, INVERT_AC_GRAD, a_1, c_1, + d_1, b_1, f_1, mesh->ystart, dump); ///////////////////////////////////////////////// // Test 2: Gaussian x-profiles, 4th order Krylov - check_laplace(++test_num, "PETSc 4th order", *invert_4th, INVERT_AC_GRAD, - INVERT_AC_GRAD, a_1, c_1, d_1, b_1, f_1, mesh->ystart, dump); + check_laplace(2, "PETSc 4th order", *invert_4th, INVERT_AC_GRAD, INVERT_AC_GRAD, a_1, + c_1, d_1, b_1, f_1, mesh->ystart, dump); //////////////////////////////////////////////////////////////////////////////////////// - // Test 3+4: Gaussian x-profiles, z-independent coefficients and compare with SPT method + // Test 3+4: Gaussian x-profiles, z-independent coefficients const Field2D a_3 = DC(a_1); const Field2D c_3 = DC(c_1); const Field2D d_3 = DC(d_1); const Field3D b_3 = forward_laplace(f_1, a_3, c_3, d_3); - check_laplace(++test_num, "with coefficients constant in z, PETSc 2nd order", *invert, + check_laplace(3, "with coefficients constant in z, PETSc 2nd order", *invert, INVERT_AC_GRAD, INVERT_AC_GRAD, a_3, c_3, d_3, b_3, f_1, mesh->ystart, dump); - Options* SPT_options = Options::getRoot()->getSection("SPT"); - auto invert_SPT = Laplacian::create(SPT_options); - - check_laplace(++test_num, "with coefficients constant in z, default solver", - *invert_SPT, INVERT_AC_GRAD, INVERT_AC_GRAD | INVERT_DC_GRAD, a_3, c_3, - d_3, b_3, f_1, mesh->ystart, dump); - ////////////////////////////////////////////// // Test 5: Cosine x-profiles, 2nd order Krylov Field3D f_5 = generate_f5(*mesh); @@ -176,16 +197,14 @@ int main(int argc, char** argv) { const Field3D b_5 = forward_laplace(f_5, a_5, c_5, d_5); - check_laplace(++test_num, "different profiles, PETSc 2nd order", *invert, - INVERT_AC_GRAD, INVERT_AC_GRAD, a_5, c_5, d_5, b_5, f_5, mesh->ystart, - dump); + check_laplace(5, "different profiles, PETSc 2nd order", *invert, INVERT_AC_GRAD, + INVERT_AC_GRAD, a_5, c_5, d_5, b_5, f_5, mesh->ystart, dump); ////////////////////////////////////////////// // Test 6: Cosine x-profiles, 4th order Krylov - check_laplace(++test_num, "different profiles, PETSc 4th order", *invert_4th, - INVERT_AC_GRAD, INVERT_AC_GRAD, a_5, c_5, d_5, b_5, f_5, mesh->ystart, - dump); + check_laplace(6, "different profiles, PETSc 4th order", *invert_4th, INVERT_AC_GRAD, + INVERT_AC_GRAD, a_5, c_5, d_5, b_5, f_5, mesh->ystart, dump); ////////////////////////////////////////////////////////////////////////////////////// // Test 7+8: Cosine x-profiles, z-independent coefficients and compare with SPT method @@ -195,16 +214,11 @@ int main(int argc, char** argv) { const Field2D d_7 = DC(d_5); const Field3D b_7 = forward_laplace(f_5, a_7, c_7, d_7); - check_laplace(++test_num, + check_laplace(7, "different profiles, with coefficients constant in z, PETSc 2nd order", *invert, INVERT_AC_GRAD, INVERT_AC_GRAD, a_7, c_7, d_7, b_7, f_5, mesh->ystart, dump); - check_laplace(++test_num, - "different profiles, with coefficients constant in z, default solver", - *invert_SPT, INVERT_AC_GRAD, INVERT_AC_GRAD | INVERT_DC_GRAD, a_7, c_7, - d_7, b_7, f_5, mesh->ystart, dump); - // Write and close the output file bout::writeDefaultOutputFile(dump); @@ -219,13 +233,11 @@ int main(int argc, char** argv) { BoutReal max_error_at_ystart(const Field3D& error) { const auto* mesh = error.getMesh(); - BoutReal local_max_error = error(mesh->xstart, mesh->ystart, 0); + BoutReal local_max_error = 0.0; for (int jx = mesh->xstart; jx <= mesh->xend; jx++) { for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { - if (local_max_error < error(jx, mesh->ystart, jz)) { - local_max_error = error(jx, mesh->ystart, jz); - } + local_max_error = std::max(local_max_error, error(jx, mesh->ystart, jz)); } } From ec6a97f8ebd6a0909939091b663825c97d396380 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 24 Feb 2026 12:07:38 +0000 Subject: [PATCH 06/12] tests: Refactor petsc-laplace test - use field factory, rather than explicit loops --- tests/integrated/test-petsc_laplace/runtest | 52 +- .../test-petsc_laplace/test_petsc_laplace.cxx | 490 ++++-------------- 2 files changed, 135 insertions(+), 407 deletions(-) diff --git a/tests/integrated/test-petsc_laplace/runtest b/tests/integrated/test-petsc_laplace/runtest index 5496987128..b4379f7754 100755 --- a/tests/integrated/test-petsc_laplace/runtest +++ b/tests/integrated/test-petsc_laplace/runtest @@ -8,25 +8,21 @@ # requires: all_tests # cores: 4 -# Variables to compare -from __future__ import print_function -from builtins import str - -vars = [ - ("max_error1", 2.0e-4), - ("max_error2", 2.0e-4), - ("max_error3", 2.0e-4), - ("max_error5", 2.0e-4), - ("max_error6", 2.0e-5), - ("max_error7", 2.0e-4), -] -# tol = 1e-4 # Absolute (?) tolerance - from boututils.run_wrapper import build_and_log, shell, launch_safe -from boutdata.collect import collect - -# import numpy as np -from sys import stdout, exit +from boutdata.collect import collect, create_cache + +import pathlib +from sys import exit + +errors = [ + "max_error1", + "max_error2", + "max_error3", + "max_error5", + "max_error6", + "max_error7", +] +tol = 2e-4 # Absolute (?) tolerance build_and_log("PETSc Laplacian inversion test") @@ -34,29 +30,25 @@ print("Running PETSc Laplacian inversion test") success = True for nproc in [1, 2, 4]: - # nxpe = 1 - # if nproc > 2: - # nxpe = 2 - cmd = "./test_petsc_laplace" shell("rm data/BOUT.dmp.*.nc") - print(" %d processors...." % nproc) + print(f" {nproc} processors....") s, out = launch_safe(cmd, nproc=nproc, pipe=True, verbose=True) - f = open("run.log." + str(nproc), "w") - f.write(out) - f.close() + + pathlib.Path(f"run.log.{nproc}").write_text(out) + cache = create_cache(path="data", prefix="BOUT.dmp") # Collect output data - for varname, tol in vars: - stdout.write(" Checking " + varname + " ... ") - error = collect(varname, path="data", info=False) + for varname in errors: + print(f" Checking {varname} ... ", end="") + error = collect(varname, path="data", info=False, datafile_cache=cache) if error <= 0: print("Convergence error") success = False elif error > tol: - print("Fail, maximum error is = " + str(error)) + print(f"Fail, maximum error is = {error:e}") success = False else: print("Pass") diff --git a/tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx b/tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx index e869cd1b1d..96e7994637 100644 --- a/tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx +++ b/tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx @@ -31,19 +31,22 @@ #include "bout/difops.hxx" #include "bout/field2d.hxx" #include "bout/field3d.hxx" +#include "bout/field_factory.hxx" #include "bout/invert_laplace.hxx" #include "bout/options.hxx" #include "bout/options_io.hxx" #include "bout/output.hxx" #include "bout/traits.hxx" +#include "bout/vecops.hxx" -#include "fmt/core.h" +#include "fmt/format.h" #include #include #include #include +namespace { BoutReal max_error_at_ystart(const Field3D& error); void apply_flat_boundary(Field3D& bcoef); @@ -101,14 +104,69 @@ Field3D forward_laplace(const Field3D& field, const T& acoef, const T& ccoef, } Field3D generate_f1(const Mesh& mesh); -Field3D generate_a1(const Mesh& mesh); -Field3D generate_c1(const Mesh& mesh); -Field3D generate_d1(const Mesh& mesh); -Field3D generate_f5(const Mesh& mesh); -Field3D generate_a5(const Mesh& mesh); -Field3D generate_c5(const Mesh& mesh); -Field3D generate_d5(const Mesh& mesh); +Field3D generate_d1(Mesh& mesh) { + Options options{{"p_d1", 0.512547}, {"q_d1", 0.30908712}}; + + return FieldFactory::get()->create3D( + "1. + (0.2 * exp(-50. * (x - p_d1)^2 / 4.) * sin((z - 2. * pi * q_d1) * 3.))", + &options, &mesh); +} + +Field3D generate_c1(Mesh& mesh) { + Options options{{"p_c1", 0.18439023}, {"q_c1", 0.401089473}}; + + return FieldFactory::get()->create3D( + "1. + (0.15 * exp(-50. * (x - p_c1)^2 * 2.) * sin((z - 2. * pi * q_c1) * 2.))", + &options, &mesh); +} + +Field3D generate_a1(Mesh& mesh) { + Options options{{"p_a1", 0.612547}, {"q_a1", 0.30908712}}; + + return FieldFactory::get()->create3D( + "-1. + (0.1 * exp(-50. * (x - p_a1)^2 * 2.5) * sin((z - 2. * pi * q_a1) * 7.))", + &options, &mesh); +} + +Field3D generate_f5(Mesh& mesh) { + Options options{{"p_f5", 0.623901}, {"q_f5", 0.01209489}}; + + auto result = FieldFactory::get()->create3D( + "exp(-((50. * (x - p_f5)^2) + 1. - cos((z - 2. * pi * q_f5))))" + "- 50 * (2. * p_f5 * exp(-50. * (-p_f5)^2) * x" + " + (-p_f5 * exp(-50. * (-p_f5)^2) - (1 - p_f5) * exp(-50. * (1 - " + "p_f5)^2)) * x^2)" + " * exp(-(1. - cos(z - 2. * pi * q_f5)))", + &options, &mesh); + result.applyBoundary("neumann"); + checkData(result); + return result; +} + +Field3D generate_d5(Mesh& mesh) { + Options options{{"p_d5", 0.63298589}, {"q_d5", 0.889237890}}; + + return FieldFactory::get()->create3D( + "1. + (p_d5 * cos(2. * pi * x) * sin((z - 2 * pi * q_d5) * 3.))", &options, &mesh); +} + +Field3D generate_c5(Mesh& mesh) { + Options options{{"p_c5", 0.160983834}, {"q_c5", 0.73050121087}}; + + return FieldFactory::get()->create3D( + "1. + (p_c5 * cos(2. * pi * x * 5) * sin((z - 2 * pi * q_c5) * 2.))", &options, + &mesh); +} + +Field3D generate_a5(Mesh& mesh) { + Options options{{"p_a5", 0.5378950}, {"q_a5", 0.2805870}}; + + return FieldFactory::get()->create3D( + "-1. + (p_a5 * cos(2. * pi * x * 2.) * sin((z - 2. * pi * q_a5) * 7.))", &options, + &mesh); +} +} // namespace int main(int argc, char** argv) { @@ -151,6 +209,25 @@ int main(int argc, char** argv) { // Solving equations of the form d*Delp2(f) + 1/c*Grad_perp(c).Grad_perp(f) + a*f = b for various f, a, c, d using bout::globals::mesh; + Field3D loop_x; + Field3D loop_z; + loop_x.allocate(); + loop_z.allocate(); + for (int jx = mesh->xstart; jx <= mesh->xend; jx++) { + const BoutReal x = mesh->GlobalX(jx); + for (int jy = 0; jy < mesh->LocalNy; jy++) { + for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { + const BoutReal z = mesh->GlobalZ(jz); + loop_x(jx, jy, jz) = x; + loop_z(jx, jy, jz) = z; + } + } + } + dump["loop_x"] = loop_x; + dump["loop_z"] = loop_z; + dump["exp_x"] = FieldFactory::get()->create3D("x"); + dump["exp_z"] = FieldFactory::get()->create3D("z"); + // Only Neumann x-boundary conditions are implemented so far, so test functions should be Neumann in x and periodic in z. // Use Field3D's, but solver only works on FieldPerp slices, so only use 1 y-point @@ -231,6 +308,7 @@ int main(int argc, char** argv) { return 0; } +namespace { BoutReal max_error_at_ystart(const Field3D& error) { const auto* mesh = error.getMesh(); BoutReal local_max_error = 0.0; @@ -271,234 +349,23 @@ void apply_flat_boundary(Field3D& bcoef) { } Field3D generate_f1(const Mesh& mesh) { - const BoutReal nx = mesh.GlobalNx - 2 * mesh.xstart - 1; - const BoutReal nz = mesh.GlobalNz; - constexpr BoutReal p = 0.39503274; // NOLINT constexpr BoutReal q = 0.20974396; // NOLINT Field3D result; result.allocate(); for (int jx = mesh.xstart; jx <= mesh.xend; jx++) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; - //make the gradients zero at both x-boundaries - result(jx, jy, jz) = 0. - + exp(-(100. * pow(x - p, 2) + 1. - cos(2. * PI * (z - q)))) - - 50. - * (2. * p * exp(-100. * pow(-p, 2)) * x - + (-p * exp(-100. * pow(-p, 2)) - - (1 - p) * exp(-100. * pow(1 - p, 2))) - * pow(x, 2)) - * exp(-(1. - cos(2. * PI * (z - q)))); - } - } - } - if (mesh.firstX()) { - for (int jx = mesh.xstart - 1; jx >= 0; jx--) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; - - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; - //make the gradients zero at both x-boundaries - result(jx, jy, jz) = 0. - + exp(-(60. * pow(x - p, 2) + 1. - cos(2. * PI * (z - q)))) - - 50. - * (2. * p * exp(-60. * pow(-p, 2)) * x - + (-p * exp(-60. * pow(-p, 2)) - - (1 - p) * exp(-60. * pow(1 - p, 2))) - * pow(x, 2)) - * exp(-(1. - cos(2. * PI * (z - q)))); - } - } - } - } - if (mesh.lastX()) { - for (int jx = mesh.xend + 1; jx < mesh.LocalNx; jx++) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; - //make the gradients zero at both x-boundaries - result(jx, jy, jz) = 0. - + exp(-(60. * pow(x - p, 2) + 1. - cos(2. * PI * (z - q)))) - - 50. - * (2. * p * exp(-60. * pow(-p, 2)) * x - + (-p * exp(-60. * pow(-p, 2)) - - (1 - p) * exp(-60. * pow(1 - p, 2))) - * pow(x, 2)) - * exp(-(1. - cos(2. * PI * (z - q)))); - } - } - } - } - - checkData(result); - result.applyBoundary("neumann"); - return result; -} - -Field3D generate_d1(const Mesh& mesh) { - const BoutReal nx = mesh.GlobalNx - 2 * mesh.xstart - 1; - const BoutReal nz = mesh.GlobalNz; - - constexpr BoutReal p = 0.512547; // NOLINT - constexpr BoutReal q = 0.30908712; // NOLINT - Field3D result; - result.allocate(); - for (int jx = mesh.xstart; jx <= mesh.xend; jx++) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; - result(jx, jy, jz) = - 1. + 0.2 * exp(-50. * pow(x - p, 2) / 4.) * sin(2. * PI * (z - q) * 3.); - } - } - } - if (mesh.firstX()) { - for (int jx = mesh.xstart - 1; jx >= 0; jx--) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; - result(jx, jy, jz) = - 1. + 0.2 * exp(-50. * pow(x - p, 2) / 4.) * sin(2. * PI * (z - q) * 3.); - } - } - } - } - if (mesh.lastX()) { - for (int jx = mesh.xend + 1; jx < mesh.LocalNx; jx++) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; - result(jx, jy, jz) = - 1. + 0.2 * exp(-50. * pow(x - p, 2) / 4.) * sin(2. * PI * (z - q) * 3.); - } - } - } - } - checkData(result); - return result; -} - -Field3D generate_c1(const Mesh& mesh) { - const BoutReal nx = mesh.GlobalNx - 2 * mesh.xstart - 1; - const BoutReal nz = mesh.GlobalNz; - - constexpr BoutReal p = 0.18439023; // NOLINT - constexpr BoutReal q = 0.401089473; // NOLINT - Field3D result; - result.allocate(); - for (int jx = mesh.xstart; jx <= mesh.xend; jx++) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; - result(jx, jy, jz) = - 1. + 0.15 * exp(-50. * pow(x - p, 2) * 2.) * sin(2. * PI * (z - q) * 2.); - } - } - } - if (mesh.firstX()) { - for (int jx = mesh.xstart - 1; jx >= 0; jx--) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; - result(jx, jy, jz) = - 1. + 0.15 * exp(-50. * pow(x - p, 2) * 2.) * sin(2. * PI * (z - q) * 2.); - } - } - } - } - if (mesh.lastX()) { - for (int jx = mesh.xend + 1; jx < mesh.LocalNx; jx++) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; - result(jx, jy, jz) = - 1. + 0.15 * exp(-50. * pow(x - p, 2) * 2.) * sin(2. * PI * (z - q) * 2.); - } - } - } - } - - checkData(result); - return result; -} - -Field3D generate_a1(const Mesh& mesh) { - const BoutReal nx = mesh.GlobalNx - 2 * mesh.xstart - 1; - const BoutReal nz = mesh.GlobalNz; - - constexpr BoutReal p = 0.612547; // NOLINT - constexpr BoutReal q = 0.30908712; // NOLINT - Field3D result; - result.allocate(); - for (int jx = mesh.xstart; jx <= mesh.xend; jx++) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; + const BoutReal x = mesh.GlobalX(jx); for (int jy = 0; jy < mesh.LocalNy; jy++) { for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; - result(jx, jy, jz) = - -1. + 0.1 * exp(-50. * pow(x - p, 2) * 2.5) * sin(2. * PI * (z - q) * 7.); - } - } - } - if (mesh.firstX()) { - for (int jx = mesh.xstart - 1; jx >= 0; jx--) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; - result(jx, jy, jz) = - -1. + 0.1 * exp(-50. * pow(x - p, 2) * 2.5) * sin(2. * PI * (z - q) * 7.); - } - } - } - } - if (mesh.lastX()) { - for (int jx = mesh.xend + 1; jx < mesh.LocalNx; jx++) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; - result(jx, jy, jz) = - -1. + 0.1 * exp(-50. * pow(x - p, 2) * 2.5) * sin(2. * PI * (z - q) * 7.); - } - } - } - } - - checkData(result); - return result; -} - -Field3D generate_f5(const Mesh& mesh) { - const BoutReal nx = mesh.GlobalNx - 2 * mesh.xstart - 1; - const BoutReal nz = mesh.GlobalNz; - constexpr BoutReal p = 0.623901; // NOLINT - constexpr BoutReal q = 0.01209489; // NOLINT - Field3D result; - result.allocate(); - for (int jx = mesh.xstart; jx <= mesh.xend; jx++) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; + const BoutReal z = mesh.GlobalZ(jz); //make the gradients zero at both x-boundaries result(jx, jy, jz) = - 0. + exp(-(50. * pow(x - p, 2) + 1. - cos(2. * PI * (z - q)))) + 0. + exp(-((100. * pow(x - p, 2)) + 1. - cos(2. * PI * (z - q)))) - 50. - * (2. * p * exp(-50. * pow(-p, 2)) * x - + (-p * exp(-50. * pow(-p, 2)) - (1 - p) * exp(-50. * pow(1 - p, 2))) + * (2. * p * exp(-100. * pow(-p, 2)) * x + + (-p * exp(-100. * pow(-p, 2)) + - (1 - p) * exp(-100. * pow(1 - p, 2))) * pow(x, 2)) * exp(-(1. - cos(2. * PI * (z - q)))); } @@ -506,176 +373,45 @@ Field3D generate_f5(const Mesh& mesh) { } if (mesh.firstX()) { for (int jx = mesh.xstart - 1; jx >= 0; jx--) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; - //make the gradients zero at both x-boundaries - result(jx, jy, jz) = 0. - + exp(-(50. * pow(x - p, 2) + 1. - cos(2. * PI * (z - q)))) - - 50. - * (2. * p * exp(-50. * pow(-p, 2)) * x - + (-p * exp(-50. * pow(-p, 2)) - - (1 - p) * exp(-50. * pow(1 - p, 2))) - * pow(x, 2)) - * exp(-(1. - cos(2. * PI * (z - q)))); - } - } - } - } - if (mesh.lastX()) { - for (int jx = mesh.xend + 1; jx < mesh.LocalNx; jx++) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; + const BoutReal x = mesh.GlobalX(jx); for (int jy = 0; jy < mesh.LocalNy; jy++) { for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; + const BoutReal z = mesh.GlobalZ(jz); //make the gradients zero at both x-boundaries - result(jx, jy, jz) = 0. - + exp(-(50. * pow(x - p, 2) + 1. - cos(2. * PI * (z - q)))) - - 50. - * (2. * p * exp(-50. * pow(-p, 2)) * x - + (-p * exp(-50. * pow(-p, 2)) - - (1 - p) * exp(-50. * pow(1 - p, 2))) - * pow(x, 2)) - * exp(-(1. - cos(2. * PI * (z - q)))); - } - } - } - } - result.applyBoundary("neumann"); - checkData(result); - return result; -} - -Field3D generate_d5(const Mesh& mesh) { - const BoutReal nx = mesh.GlobalNx - 2 * mesh.xstart - 1; - const BoutReal nz = mesh.GlobalNz; - constexpr BoutReal p = 0.63298589; // NOLINT - constexpr BoutReal q = 0.889237890; // NOLINT - Field3D result; - result.allocate(); - for (int jx = mesh.xstart; jx <= mesh.xend; jx++) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; - result(jx, jy, jz) = 1. + p * cos(2. * PI * x) * sin(2. * PI * (z - q) * 3.); - } - } - } - if (mesh.firstX()) { - for (int jx = mesh.xstart - 1; jx >= 0; jx--) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; - result(jx, jy, jz) = 1. + p * cos(2. * PI * x) * sin(2. * PI * (z - q) * 3.); - } - } - } - } - if (mesh.lastX()) { - for (int jx = mesh.xend + 1; jx < mesh.LocalNx; jx++) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; - result(jx, jy, jz) = 1. + p * cos(2. * PI * x) * sin(2. * PI * (z - q) * 3.); - } - } - } - } - checkData(result); - return result; -} - -Field3D generate_c5(const Mesh& mesh) { - const BoutReal nx = mesh.GlobalNx - 2 * mesh.xstart - 1; - const BoutReal nz = mesh.GlobalNz; - constexpr BoutReal p = 0.160983834; // NOLINT - constexpr BoutReal q = 0.73050121087; // NOLINT - - Field3D result; - - result.allocate(); - for (int jx = mesh.xstart; jx <= mesh.xend; jx++) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; - result(jx, jy, jz) = 1. + p * cos(2. * PI * x * 5) * sin(2. * PI * (z - q) * 2.); - } - } - } - if (mesh.firstX()) { - for (int jx = mesh.xstart - 1; jx >= 0; jx--) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; result(jx, jy, jz) = - 1. + p * cos(2. * PI * x * 5) * sin(2. * PI * (z - q) * 2.); + 0. + exp(-((60. * pow(x - p, 2)) + 1. - cos(2. * PI * (z - q)))) + - 50. + * (2. * p * exp(-60. * pow(-p, 2)) * x + + (-p * exp(-60. * pow(-p, 2)) + - (1 - p) * exp(-60. * pow(1 - p, 2))) + * pow(x, 2)) + * exp(-(1. - cos(2. * PI * (z - q)))); } } } } if (mesh.lastX()) { for (int jx = mesh.xend + 1; jx < mesh.LocalNx; jx++) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; + const BoutReal x = mesh.GlobalX(jx); for (int jy = 0; jy < mesh.LocalNy; jy++) { for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; + const BoutReal z = mesh.GlobalZ(jz); + //make the gradients zero at both x-boundaries result(jx, jy, jz) = - 1. + p * cos(2. * PI * x * 5) * sin(2. * PI * (z - q) * 2.); + 0. + exp(-((60. * pow(x - p, 2)) + 1. - cos(2. * PI * (z - q)))) + - 50. + * (2. * p * exp(-60. * pow(-p, 2)) * x + + (-p * exp(-60. * pow(-p, 2)) + - (1 - p) * exp(-60. * pow(1 - p, 2))) + * pow(x, 2)) + * exp(-(1. - cos(2. * PI * (z - q)))); } } } } - checkData(result); - return result; -} -Field3D generate_a5(const Mesh& mesh) { - const BoutReal nx = mesh.GlobalNx - 2 * mesh.xstart - 1; - const BoutReal nz = mesh.GlobalNz; - constexpr BoutReal p = 0.5378950; // NOLINT - constexpr BoutReal q = 0.2805870; // NOLINT - Field3D result; - result.allocate(); - for (int jx = mesh.xstart; jx <= mesh.xend; jx++) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; - result(jx, jy, jz) = - -1. + p * cos(2. * PI * x * 2.) * sin(2. * PI * (z - q) * 7.); - } - } - } - if (mesh.firstX()) { - for (int jx = mesh.xstart - 1; jx >= 0; jx--) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; - result(jx, jy, jz) = - -1. + p * cos(2. * PI * x * 2.) * sin(2. * PI * (z - q) * 7.); - } - } - } - } - if (mesh.lastX()) { - for (int jx = mesh.xend + 1; jx < mesh.LocalNx; jx++) { - const BoutReal x = BoutReal(mesh.getGlobalXIndex(jx) - mesh.xstart) / nx; - for (int jy = 0; jy < mesh.LocalNy; jy++) { - for (int jz = mesh.zstart; jz <= mesh.zend; jz++) { - const BoutReal z = BoutReal(jz) / nz; - result(jx, jy, jz) = - -1. + p * cos(2. * PI * x * 2.) * sin(2. * PI * (z - q) * 7.); - } - } - } - } checkData(result); + result.applyBoundary("neumann"); return result; } +} // namespace From 6e0ac334eea7efb1568d0cf20e3eddf0cd9f1098 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 24 Feb 2026 13:18:11 +0000 Subject: [PATCH 07/12] Reduce some duplicated loops in `LaplacePetsc` ctor --- .../laplace/impls/petsc/petsc_laplace.cxx | 142 ++++++------------ 1 file changed, 43 insertions(+), 99 deletions(-) diff --git a/src/invert/laplace/impls/petsc/petsc_laplace.cxx b/src/invert/laplace/impls/petsc/petsc_laplace.cxx index c6a6712f4d..d9076dc6a5 100644 --- a/src/invert/laplace/impls/petsc/petsc_laplace.cxx +++ b/src/invert/laplace/impls/petsc/petsc_laplace.cxx @@ -132,6 +132,9 @@ LaplacePetsc::LaplacePetsc(Options* opt, const CELL_LOC loc, Mesh* mesh_in, MatSetSizes(MatA, localN, localN, size, size); MatSetFromOptions(MatA); + const bool first_x = localmesh->firstX(); + const bool last_x = localmesh->lastX(); + /* Pre allocate memory * nnz denotes an array containing the number of non-zeros in the various rows * for @@ -139,109 +142,49 @@ LaplacePetsc::LaplacePetsc(Options* opt, const CELL_LOC loc, Mesh* mesh_in, * o_nnz - The off-diagonal terms in the matrix (needed when running in * parallel) */ - PetscInt *d_nnz, *o_nnz; + PetscInt *d_nnz = nullptr; + PetscInt *o_nnz = nullptr; PetscMalloc((localN) * sizeof(PetscInt), &d_nnz); PetscMalloc((localN) * sizeof(PetscInt), &o_nnz); if (fourth_order) { - // first and last 2*localmesh-LocalNz entries are the edge x-values that (may) have 'off-diagonal' components (i.e. on another processor) - if (localmesh->firstX() && localmesh->lastX()) { - for (int i = localmesh->zstart; i <= localmesh->zend; i++) { - d_nnz[i] = 15; - d_nnz[localN - 1 - i] = 15; - o_nnz[i] = 0; - o_nnz[localN - 1 - i] = 0; - } - for (int i = (localmesh->LocalNz); i < 2 * (localmesh->LocalNz); i++) { - d_nnz[i] = 20; - d_nnz[localN - 1 - i] = 20; - o_nnz[i] = 0; - o_nnz[localN - 1 - i] = 0; - } - } else if (localmesh->firstX()) { - for (int i = localmesh->zstart; i <= localmesh->zend; i++) { - d_nnz[i] = 15; - d_nnz[localN - 1 - i] = 15; - o_nnz[i] = 0; - o_nnz[localN - 1 - i] = 10; - } - for (int i = (localmesh->LocalNz); i < 2 * (localmesh->LocalNz); i++) { - d_nnz[i] = 20; - d_nnz[localN - 1 - i] = 20; - o_nnz[i] = 0; - o_nnz[localN - 1 - i] = 5; - } - } else if (localmesh->lastX()) { - for (int i = localmesh->zstart; i <= localmesh->zend; i++) { - d_nnz[i] = 15; - d_nnz[localN - 1 - i] = 15; - o_nnz[i] = 10; - o_nnz[localN - 1 - i] = 0; - } - for (int i = (localmesh->LocalNz); i < 2 * (localmesh->LocalNz); i++) { - d_nnz[i] = 20; - d_nnz[localN - 1 - i] = 20; - o_nnz[i] = 5; - o_nnz[localN - 1 - i] = 0; - } - } else { - for (int i = localmesh->zstart; i <= localmesh->zend; i++) { - d_nnz[i] = 15; - d_nnz[localN - 1 - i] = 15; - o_nnz[i] = 10; - o_nnz[localN - 1 - i] = 10; - } - for (int i = (localmesh->LocalNz); i < 2 * (localmesh->LocalNz); i++) { - d_nnz[i] = 20; - d_nnz[localN - 1 - i] = 20; - o_nnz[i] = 5; - o_nnz[localN - 1 - i] = 5; - } + // first and last 2*localmesh-LocalNz entries are the edge x-values that + // (may) have 'off-diagonal' components (i.e. on another processor) + + const int first_first_off = first_x ? 0 : 10; + const int first_last_off = last_x ? 0 : 10; + const int second_first_off = first_x ? 0 : 5; + const int second_last_off = last_x ? 0 : 5; + + for (int i = localmesh->zstart; i <= localmesh->zend; i++) { + d_nnz[i] = 15; + d_nnz[localN - 1 - i] = 15; + o_nnz[i] = first_first_off; + o_nnz[localN - 1 - i] = first_last_off; + } + for (int i = (localmesh->LocalNz); i < 2 * (localmesh->LocalNz); i++) { + d_nnz[i] = 20; + d_nnz[localN - 1 - i] = 20; + o_nnz[i] = second_first_off; + o_nnz[localN - 1 - i] = second_last_off; } - for (int i = 2 * (localmesh->LocalNz); i < localN - 2 * ((localmesh->LocalNz)); i++) { + for (int i = 2 * (localmesh->LocalNz); i < localN - (2 * localmesh->LocalNz); i++) { d_nnz[i] = 25; d_nnz[localN - 1 - i] = 25; o_nnz[i] = 0; o_nnz[localN - 1 - i] = 0; } - - // Use d_nnz and o_nnz for preallocating the matrix - if (localmesh->firstX() && localmesh->lastX()) { - // Only one processor in X - MatSeqAIJSetPreallocation(MatA, 0, d_nnz); - } else { - MatMPIAIJSetPreallocation(MatA, 0, d_nnz, 0, o_nnz); - } } else { - // first and last localmesh->LocalNz entries are the edge x-values that (may) have 'off-diagonal' components (i.e. on another processor) - if (localmesh->firstX() && localmesh->lastX()) { - for (int i = localmesh->zstart; i <= localmesh->zend; i++) { - d_nnz[i] = 6; - d_nnz[localN - 1 - i] = 6; - o_nnz[i] = 0; - o_nnz[localN - 1 - i] = 0; - } - } else if (localmesh->firstX()) { - for (int i = localmesh->zstart; i <= localmesh->zend; i++) { - d_nnz[i] = 6; - d_nnz[localN - 1 - i] = 6; - o_nnz[i] = 0; - o_nnz[localN - 1 - i] = 3; - } - } else if (localmesh->lastX()) { - for (int i = localmesh->zstart; i <= localmesh->zend; i++) { - d_nnz[i] = 6; - d_nnz[localN - 1 - i] = 6; - o_nnz[i] = 3; - o_nnz[localN - 1 - i] = 0; - } - } else { - for (int i = localmesh->zstart; i <= localmesh->zend; i++) { - d_nnz[i] = 6; - d_nnz[localN - 1 - i] = 6; - o_nnz[i] = 3; - o_nnz[localN - 1 - i] = 3; - } + // first and last localmesh->LocalNz entries are the edge x-values that + // (may) have 'off-diagonal' components (i.e. on another processor) + const int first_off = first_x ? 0 : 3; + const int last_off = last_x ? 0 : 3; + + for (int i = localmesh->zstart; i <= localmesh->zend; i++) { + d_nnz[i] = 6; + d_nnz[localN - 1 - i] = 6; + o_nnz[i] = first_off; + o_nnz[localN - 1 - i] = last_off; } for (int i = localmesh->LocalNz; i < localN - (localmesh->LocalNz); i++) { @@ -250,14 +193,15 @@ LaplacePetsc::LaplacePetsc(Options* opt, const CELL_LOC loc, Mesh* mesh_in, o_nnz[i] = 0; o_nnz[localN - 1 - i] = 0; } - - // Use d_nnz and o_nnz for preallocating the matrix - if (localmesh->firstX() && localmesh->lastX()) { - MatSeqAIJSetPreallocation(MatA, 0, d_nnz); - } else { - MatMPIAIJSetPreallocation(MatA, 0, d_nnz, 0, o_nnz); - } } + // Use d_nnz and o_nnz for preallocating the matrix + if (localmesh->firstX() && localmesh->lastX()) { + // Only one processor in X + MatSeqAIJSetPreallocation(MatA, 0, d_nnz); + } else { + MatMPIAIJSetPreallocation(MatA, 0, d_nnz, 0, o_nnz); + } + // Free the d_nnz and o_nnz arrays, as these are will not be used anymore PetscFree(d_nnz); PetscFree(o_nnz); From a5e3f366a79bec2ee1dc3dc8744d818b177c0790 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 24 Feb 2026 13:45:33 +0000 Subject: [PATCH 08/12] Use `std::vector` instead of `PetscMalloc` --- .../laplace/impls/petsc/petsc_laplace.cxx | 31 ++++++++----------- 1 file changed, 13 insertions(+), 18 deletions(-) diff --git a/src/invert/laplace/impls/petsc/petsc_laplace.cxx b/src/invert/laplace/impls/petsc/petsc_laplace.cxx index d9076dc6a5..1f422ad19d 100644 --- a/src/invert/laplace/impls/petsc/petsc_laplace.cxx +++ b/src/invert/laplace/impls/petsc/petsc_laplace.cxx @@ -38,6 +38,8 @@ #include #include +#include + #define KSP_RICHARDSON "richardson" #define KSP_CHEBYSHEV "chebyshev" #define KSP_CG "cg" @@ -134,18 +136,14 @@ LaplacePetsc::LaplacePetsc(Options* opt, const CELL_LOC loc, Mesh* mesh_in, const bool first_x = localmesh->firstX(); const bool last_x = localmesh->lastX(); + // Pre allocate memory. See MatMPIAIJSetPreallocation docs for more info + // Number of non-zero elements in each row of matrix owned by this processor + // ("diagonal submatrices") + std::vector d_nnz(localN, 0); + // Number of non-zero elements in each row of matrix owned by other processors + // ("offdiagonal submatrices") + std::vector o_nnz(localN, 0); - /* Pre allocate memory - * nnz denotes an array containing the number of non-zeros in the various rows - * for - * d_nnz - The diagonal terms in the matrix - * o_nnz - The off-diagonal terms in the matrix (needed when running in - * parallel) - */ - PetscInt *d_nnz = nullptr; - PetscInt *o_nnz = nullptr; - PetscMalloc((localN) * sizeof(PetscInt), &d_nnz); - PetscMalloc((localN) * sizeof(PetscInt), &o_nnz); if (fourth_order) { // first and last 2*localmesh-LocalNz entries are the edge x-values that // (may) have 'off-diagonal' components (i.e. on another processor) @@ -194,17 +192,14 @@ LaplacePetsc::LaplacePetsc(Options* opt, const CELL_LOC loc, Mesh* mesh_in, o_nnz[localN - 1 - i] = 0; } } - // Use d_nnz and o_nnz for preallocating the matrix - if (localmesh->firstX() && localmesh->lastX()) { + + if (first_x and last_x) { // Only one processor in X - MatSeqAIJSetPreallocation(MatA, 0, d_nnz); + MatSeqAIJSetPreallocation(MatA, 0, d_nnz.data()); } else { - MatMPIAIJSetPreallocation(MatA, 0, d_nnz, 0, o_nnz); + MatMPIAIJSetPreallocation(MatA, 0, d_nnz.data(), 0, o_nnz.data()); } - // Free the d_nnz and o_nnz arrays, as these are will not be used anymore - PetscFree(d_nnz); - PetscFree(o_nnz); // Sets up the internal matrix data structures for the later use. MatSetUp(MatA); From 0aa79db9bbe334c1c17f5451260612a0bd6c2612 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 24 Feb 2026 13:51:24 +0000 Subject: [PATCH 09/12] Set member variables in `LaplacePetsc` initialiser --- .clang-tidy | 2 +- .../laplace/impls/petsc/petsc_laplace.cxx | 102 +++++++----------- .../laplace/impls/petsc/petsc_laplace.hxx | 18 ++-- 3 files changed, 52 insertions(+), 70 deletions(-) diff --git a/.clang-tidy b/.clang-tidy index cf29b759cd..121c3f3a60 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -20,7 +20,7 @@ CheckOptions: value: 'MPI_Comm' - key: misc-include-cleaner.IgnoreHeaders - value: 'adios2/.*;bits/.*;cpptrace/.*' + value: 'adios2/.*;bits/.*;cpptrace/.*;petsc.*\.h' --- Disabled checks and reasons: diff --git a/src/invert/laplace/impls/petsc/petsc_laplace.cxx b/src/invert/laplace/impls/petsc/petsc_laplace.cxx index 1f422ad19d..878d18b4f3 100644 --- a/src/invert/laplace/impls/petsc/petsc_laplace.cxx +++ b/src/invert/laplace/impls/petsc/petsc_laplace.cxx @@ -31,7 +31,10 @@ #include "petsc_laplace.hxx" #include +#include #include +#include +#include #include #include #include @@ -53,21 +56,43 @@ #define KSP_BICG "bicg" #define KSP_PREONLY "preonly" -static PetscErrorCode laplacePCapply(PC pc, Vec x, Vec y) { +namespace { +PetscErrorCode laplacePCapply(PC pc, Vec x, Vec y) { PetscFunctionBegin; // NOLINT LaplacePetsc* laplace = nullptr; - const int ierr = PCShellGetContext(pc, reinterpret_cast(&laplace)); // NOLINT - CHKERRQ(ierr); + CHKERRQ(PCShellGetContext(pc, reinterpret_cast(&laplace))); // NOLINT PetscFunctionReturn(laplace->precon(x, y)); // NOLINT } +} // namespace LaplacePetsc::LaplacePetsc(Options* opt, const CELL_LOC loc, Mesh* mesh_in, - Solver* UNUSED(solver)) + [[maybe_unused]] Solver* solver) : Laplacian(opt, loc, mesh_in), A(0.0), C1(1.0), C2(1.0), D(1.0), Ex(0.0), Ez(0.0), - issetD(false), issetC(false), issetE(false), - lib(opt == nullptr ? &(Options::root()["laplace"]) : opt) { + issetD(false), issetC(false), issetE(false), comm(localmesh->getXcomm()), + opts(opt == nullptr ? &(Options::root()["laplace"]) : opt), + // WARNING: only a few of these options actually make sense: see the + // PETSc documentation to work out which they are (possibly + // pbjacobi, sor might be useful choices?) + ksptype((*opts)["ksptype"].doc("KSP solver type").withDefault(KSPGMRES)), + pctype((*opts)["pctype"] + .doc("Preconditioner type. See the PETSc documentation for options") + .withDefault("none")), + richardson_damping_factor((*opts)["richardson_damping_factor"].withDefault(1.0)), + chebyshev_max((*opts)["chebyshev_max"].withDefault(100)), + chebyshev_min((*opts)["chebyshev_min"].withDefault(0.01)), + gmres_max_steps((*opts)["gmres_max_steps"].withDefault(30)), + rtol((*opts)["rtol"].doc("Relative tolerance for KSP solver").withDefault(1e-5)), + atol((*opts)["atol"].doc("Absolute tolerance for KSP solver").withDefault(1e-50)), + dtol((*opts)["dtol"].doc("Divergence tolerance for KSP solver").withDefault(1e5)), + maxits( + (*opts)["maxits"].doc("Maximum number of KSP iterations").withDefault(100000)), + direct((*opts)["direct"].doc("Use direct (LU) solver?").withDefault(false)), + fourth_order( + (*opts)["fourth_order"].doc("Use fourth order stencil").withDefault(false)), + lib(opts) { + A.setLocation(location); C1.setLocation(location); C2.setLocation(location); @@ -75,13 +100,6 @@ LaplacePetsc::LaplacePetsc(Options* opt, const CELL_LOC loc, Mesh* mesh_in, Ex.setLocation(location); Ez.setLocation(location); - // Get Options in Laplace Section - if (!opt) { - opts = Options::getRoot()->getSection("laplace"); - } else { - opts = opt; - } - #if CHECK > 0 // Checking flags are set to something which is not implemented checkFlags(); @@ -93,24 +111,23 @@ LaplacePetsc::LaplacePetsc(Options* opt, const CELL_LOC loc, Mesh* mesh_in, } #endif - // Get communicator for group of processors in X - all points in z-x plane for fixed y. - comm = localmesh->getXcomm(); + const bool first_x = localmesh->firstX(); + const bool last_x = localmesh->lastX(); // Need to determine local size to use based on prior parallelisation // Coefficient values are stored only on local processors. localN = (localmesh->xend - localmesh->xstart + 1) * (localmesh->LocalNz); - if (localmesh->firstX()) { - localN += - localmesh->xstart - * (localmesh->LocalNz); // If on first processor add on width of boundary region + if (first_x) { + // If on first processor add on width of boundary region + localN += localmesh->xstart * (localmesh->LocalNz); } - if (localmesh->lastX()) { - localN += - localmesh->xstart - * (localmesh->LocalNz); // If on last processor add on width of boundary region + if (last_x) { + // If on last processor add on width of boundary region + localN += localmesh->xstart * (localmesh->LocalNz); } // Calculate 'size' (the total number of points in physical grid) + size = localN; if (bout::globals::mpi->MPI_Allreduce(&localN, &size, 1, MPI_INT, MPI_SUM, comm) != MPI_SUCCESS) { throw BoutException("Error in MPI_Allreduce during LaplacePetsc initialisation"); @@ -126,16 +143,11 @@ LaplacePetsc::LaplacePetsc(Options* opt, const CELL_LOC loc, Mesh* mesh_in, VecSetFromOptions(xs); VecDuplicate(xs, &bs); - // Get 4th order solver switch - opts->get("fourth_order", fourth_order, false); - // Set size of (the PETSc) Matrix on each processor to localN x localN MatCreate(comm, &MatA); MatSetSizes(MatA, localN, localN, size, size); MatSetFromOptions(MatA); - const bool first_x = localmesh->firstX(); - const bool last_x = localmesh->lastX(); // Pre allocate memory. See MatMPIAIJSetPreallocation docs for more info // Number of non-zero elements in each row of matrix owned by this processor // ("diagonal submatrices") @@ -206,53 +218,21 @@ LaplacePetsc::LaplacePetsc(Options* opt, const CELL_LOC loc, Mesh* mesh_in, // Declare KSP Context (abstract PETSc object that manages all Krylov methods) KSPCreate(comm, &ksp); - // Get KSP Solver Type (Generalizes Minimal RESidual is the default) - ksptype = (*opts)["ksptype"].doc("KSP solver type").withDefault(KSP_GMRES); - - // Get preconditioner type - // WARNING: only a few of these options actually make sense: see the - // PETSc documentation to work out which they are (possibly - // pbjacobi, sor might be useful choices?) - pctype = (*opts)["pctype"] - .doc("Preconditioner type. See the PETSc documentation for options") - .withDefault("none"); - // Let "user" be a synonym for "shell" if (pctype == "user") { pctype = PCSHELL; } - // Get Options specific to particular solver types - opts->get("richardson_damping_factor", richardson_damping_factor, 1.0, true); - opts->get("chebyshev_max", chebyshev_max, 100, true); - opts->get("chebyshev_min", chebyshev_min, 0.01, true); - opts->get("gmres_max_steps", gmres_max_steps, 30, true); - - // Get Tolerances for KSP solver - rtol = (*opts)["rtol"].doc("Relative tolerance for KSP solver").withDefault(1e-5); - atol = (*opts)["atol"].doc("Absolute tolerance for KSP solver").withDefault(1e-50); - dtol = (*opts)["dtol"].doc("Divergence tolerance for KSP solver").withDefault(1e5); - maxits = (*opts)["maxits"].doc("Maximum number of KSP iterations").withDefault(100000); - - // Get direct solver switch - direct = (*opts)["direct"].doc("Use direct (LU) solver?").withDefault(false); if (direct) { - output << endl - << "Using LU decompostion for direct solution of system" << endl - << endl; + output.write("\nUsing LU decompostion for direct solution of system\n"); } if (pctype == PCSHELL) { - rightprec = (*opts)["rightprec"].doc("Right preconditioning?").withDefault(true); // Options for preconditioner are in a subsection pcsolve = Laplacian::create(opts->getSection("precon")); } - - // Ensure that the matrix is constructed first time - // coefchanged = true; - // lastflag = -1; } FieldPerp LaplacePetsc::solve(const FieldPerp& b) { return solve(b, b); } diff --git a/src/invert/laplace/impls/petsc/petsc_laplace.hxx b/src/invert/laplace/impls/petsc/petsc_laplace.hxx index 1d56abd00b..611bfd6fa1 100644 --- a/src/invert/laplace/impls/petsc/petsc_laplace.hxx +++ b/src/invert/laplace/impls/petsc/petsc_laplace.hxx @@ -220,17 +220,19 @@ private: FieldPerp sol; // solution Field - // Istart is the first row of MatA owned by the process, Iend is 1 greater than the last row. - int Istart, Iend; + /// Istart is the first row of MatA owned by the process, Iend is 1 greater than the last row. + int Istart = -1; + int Iend = -1; - int meshx, meshz, size, - localN; // Mesh sizes, total size, no of points on this processor + /// Mesh sizes, total size, no of points on this processor + int meshx, meshz, size, localN; MPI_Comm comm; - Mat MatA; - Vec xs, bs; // Solution and RHS vectors - KSP ksp; + Mat MatA = nullptr; ///< Stencil matrix + Vec xs = nullptr; ///< Solution vector + Vec bs = nullptr; ///< RHS vector + KSP ksp = nullptr; ///< PETSc solver - Options* opts; // Laplace Section Options Object + Options* opts; ///< Laplace Section Options Object std::string ksptype; ///< KSP solver type std::string pctype; ///< Preconditioner type From 3d34d0cafc5b296734c863237255afa78bab0617 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 24 Feb 2026 14:02:04 +0000 Subject: [PATCH 10/12] Prepare `LaplacePetsc` for Z guards --- .../laplace/impls/petsc/petsc_laplace.cxx | 33 +++++++++++-------- 1 file changed, 19 insertions(+), 14 deletions(-) diff --git a/src/invert/laplace/impls/petsc/petsc_laplace.cxx b/src/invert/laplace/impls/petsc/petsc_laplace.cxx index 878d18b4f3..673f9c0df5 100644 --- a/src/invert/laplace/impls/petsc/petsc_laplace.cxx +++ b/src/invert/laplace/impls/petsc/petsc_laplace.cxx @@ -116,14 +116,15 @@ LaplacePetsc::LaplacePetsc(Options* opt, const CELL_LOC loc, Mesh* mesh_in, // Need to determine local size to use based on prior parallelisation // Coefficient values are stored only on local processors. - localN = (localmesh->xend - localmesh->xstart + 1) * (localmesh->LocalNz); + const auto local_nz = localmesh->zend - localmesh->zstart + 1; + localN = (localmesh->xend - localmesh->xstart + 1) * local_nz; if (first_x) { // If on first processor add on width of boundary region - localN += localmesh->xstart * (localmesh->LocalNz); + localN += localmesh->xstart * local_nz; } if (last_x) { // If on last processor add on width of boundary region - localN += localmesh->xstart * (localmesh->LocalNz); + localN += localmesh->xstart * local_nz; } // Calculate 'size' (the total number of points in physical grid) @@ -134,7 +135,7 @@ LaplacePetsc::LaplacePetsc(Options* opt, const CELL_LOC loc, Mesh* mesh_in, } // Calculate total (physical) grid dimensions - meshz = localmesh->LocalNz; + meshz = local_nz; meshx = size / meshz; // Create PETSc type of vectors for the solution and the RHS vector @@ -165,39 +166,39 @@ LaplacePetsc::LaplacePetsc(Options* opt, const CELL_LOC loc, Mesh* mesh_in, const int second_first_off = first_x ? 0 : 5; const int second_last_off = last_x ? 0 : 5; - for (int i = localmesh->zstart; i <= localmesh->zend; i++) { + for (int i = 0; i <= local_nz; i++) { d_nnz[i] = 15; d_nnz[localN - 1 - i] = 15; o_nnz[i] = first_first_off; o_nnz[localN - 1 - i] = first_last_off; } - for (int i = (localmesh->LocalNz); i < 2 * (localmesh->LocalNz); i++) { + for (int i = local_nz; i < 2 * local_nz; i++) { d_nnz[i] = 20; d_nnz[localN - 1 - i] = 20; o_nnz[i] = second_first_off; o_nnz[localN - 1 - i] = second_last_off; } - for (int i = 2 * (localmesh->LocalNz); i < localN - (2 * localmesh->LocalNz); i++) { + for (int i = 2 * local_nz; i < localN - (2 * local_nz); i++) { d_nnz[i] = 25; d_nnz[localN - 1 - i] = 25; o_nnz[i] = 0; o_nnz[localN - 1 - i] = 0; } } else { - // first and last localmesh->LocalNz entries are the edge x-values that + // first and last local_nz entries are the edge x-values that // (may) have 'off-diagonal' components (i.e. on another processor) const int first_off = first_x ? 0 : 3; const int last_off = last_x ? 0 : 3; - for (int i = localmesh->zstart; i <= localmesh->zend; i++) { + for (int i = 0; i <= local_nz; i++) { d_nnz[i] = 6; d_nnz[localN - 1 - i] = 6; o_nnz[i] = first_off; o_nnz[localN - 1 - i] = last_off; } - for (int i = localmesh->LocalNz; i < localN - (localmesh->LocalNz); i++) { + for (int i = local_nz; i < localN - local_nz; i++) { d_nnz[i] = 9; d_nnz[localN - 1 - i] = 9; o_nnz[i] = 0; @@ -639,7 +640,7 @@ FieldPerp LaplacePetsc::solve(const FieldPerp& b, const FieldPerp& x0) { } if (i != Iend) { - throw BoutException("Petsc index sanity check failed"); + throw BoutException("Petsc index sanity check failed: i={} != Iend={}", i, Iend); } // Assemble Matrix @@ -809,9 +810,13 @@ void LaplacePetsc::Element(int i, int x, int z, int xshift, int zshift, PetscSca // Need to convert LOCAL x to GLOBAL x in order to correctly calculate // PETSC Matrix Index. int xoffset = Istart / meshz; - if (Istart % meshz != 0) { - throw BoutException("Petsc index sanity check 3 failed"); +#if CHECK > 2 + const int rem = Istart % meshz; + if (rem != 0) { + throw BoutException("Petsc index sanity check 3 failed: Istart={} % meshz={} == {}", + Istart, meshz, rem); } +#endif // Calculate the row to be set int row_new = x + xshift; // should never be out of range. @@ -820,7 +825,7 @@ void LaplacePetsc::Element(int i, int x, int z, int xshift, int zshift, PetscSca } // Calculate the column to be set - int col_new = z + zshift; + int col_new = z + zshift - localmesh->zstart; if (col_new < 0) { col_new += meshz; } else if (col_new > meshz - 1) { From df99cdba5963dd293ff0cfb892630e856db938ea Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 24 Feb 2026 14:53:28 +0000 Subject: [PATCH 11/12] Apply clang-tidy fixes --- .../test-petsc_laplace/test_petsc_laplace.cxx | 48 +++++++++---------- 1 file changed, 23 insertions(+), 25 deletions(-) diff --git a/tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx b/tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx index 96e7994637..902156f4c0 100644 --- a/tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx +++ b/tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx @@ -32,6 +32,7 @@ #include "bout/field2d.hxx" #include "bout/field3d.hxx" #include "bout/field_factory.hxx" +#include "bout/globals.hxx" #include "bout/invert_laplace.hxx" #include "bout/options.hxx" #include "bout/options_io.hxx" @@ -106,7 +107,7 @@ Field3D forward_laplace(const Field3D& field, const T& acoef, const T& ccoef, Field3D generate_f1(const Mesh& mesh); Field3D generate_d1(Mesh& mesh) { - Options options{{"p_d1", 0.512547}, {"q_d1", 0.30908712}}; + const Options options{{"p_d1", 0.512547}, {"q_d1", 0.30908712}}; return FieldFactory::get()->create3D( "1. + (0.2 * exp(-50. * (x - p_d1)^2 / 4.) * sin((z - 2. * pi * q_d1) * 3.))", @@ -114,7 +115,7 @@ Field3D generate_d1(Mesh& mesh) { } Field3D generate_c1(Mesh& mesh) { - Options options{{"p_c1", 0.18439023}, {"q_c1", 0.401089473}}; + const Options options{{"p_c1", 0.18439023}, {"q_c1", 0.401089473}}; return FieldFactory::get()->create3D( "1. + (0.15 * exp(-50. * (x - p_c1)^2 * 2.) * sin((z - 2. * pi * q_c1) * 2.))", @@ -122,7 +123,7 @@ Field3D generate_c1(Mesh& mesh) { } Field3D generate_a1(Mesh& mesh) { - Options options{{"p_a1", 0.612547}, {"q_a1", 0.30908712}}; + const Options options{{"p_a1", 0.612547}, {"q_a1", 0.30908712}}; return FieldFactory::get()->create3D( "-1. + (0.1 * exp(-50. * (x - p_a1)^2 * 2.5) * sin((z - 2. * pi * q_a1) * 7.))", @@ -130,7 +131,7 @@ Field3D generate_a1(Mesh& mesh) { } Field3D generate_f5(Mesh& mesh) { - Options options{{"p_f5", 0.623901}, {"q_f5", 0.01209489}}; + const Options options{{"p_f5", 0.623901}, {"q_f5", 0.01209489}}; auto result = FieldFactory::get()->create3D( "exp(-((50. * (x - p_f5)^2) + 1. - cos((z - 2. * pi * q_f5))))" @@ -145,14 +146,14 @@ Field3D generate_f5(Mesh& mesh) { } Field3D generate_d5(Mesh& mesh) { - Options options{{"p_d5", 0.63298589}, {"q_d5", 0.889237890}}; + const Options options{{"p_d5", 0.63298589}, {"q_d5", 0.889237890}}; return FieldFactory::get()->create3D( "1. + (p_d5 * cos(2. * pi * x) * sin((z - 2 * pi * q_d5) * 3.))", &options, &mesh); } Field3D generate_c5(Mesh& mesh) { - Options options{{"p_c5", 0.160983834}, {"q_c5", 0.73050121087}}; + const Options options{{"p_c5", 0.160983834}, {"q_c5", 0.73050121087}}; return FieldFactory::get()->create3D( "1. + (p_c5 * cos(2. * pi * x * 5) * sin((z - 2 * pi * q_c5) * 2.))", &options, @@ -160,7 +161,7 @@ Field3D generate_c5(Mesh& mesh) { } Field3D generate_a5(Mesh& mesh) { - Options options{{"p_a5", 0.5378950}, {"q_a5", 0.2805870}}; + const Options options{{"p_a5", 0.5378950}, {"q_a5", 0.2805870}}; return FieldFactory::get()->create3D( "-1. + (p_a5 * cos(2. * pi * x * 2.) * sin((z - 2. * pi * q_a5) * 7.))", &options, @@ -362,12 +363,11 @@ Field3D generate_f1(const Mesh& mesh) { //make the gradients zero at both x-boundaries result(jx, jy, jz) = 0. + exp(-((100. * pow(x - p, 2)) + 1. - cos(2. * PI * (z - q)))) - - 50. - * (2. * p * exp(-100. * pow(-p, 2)) * x - + (-p * exp(-100. * pow(-p, 2)) - - (1 - p) * exp(-100. * pow(1 - p, 2))) - * pow(x, 2)) - * exp(-(1. - cos(2. * PI * (z - q)))); + - (50. + * (2. * p * exp(-100. * pow(-p, 2)) * x + + (-p * exp(-100. * pow(-p, 2)) - (1 - p) * exp(-100. * pow(1 - p, 2))) + * pow(x, 2)) + * exp(-(1. - cos(2. * PI * (z - q))))); } } } @@ -380,12 +380,11 @@ Field3D generate_f1(const Mesh& mesh) { //make the gradients zero at both x-boundaries result(jx, jy, jz) = 0. + exp(-((60. * pow(x - p, 2)) + 1. - cos(2. * PI * (z - q)))) - - 50. - * (2. * p * exp(-60. * pow(-p, 2)) * x - + (-p * exp(-60. * pow(-p, 2)) - - (1 - p) * exp(-60. * pow(1 - p, 2))) - * pow(x, 2)) - * exp(-(1. - cos(2. * PI * (z - q)))); + - (50. + * (2. * p * exp(-60. * pow(-p, 2)) * x + + (-p * exp(-60. * pow(-p, 2)) - (1 - p) * exp(-60. * pow(1 - p, 2))) + * pow(x, 2)) + * exp(-(1. - cos(2. * PI * (z - q))))); } } } @@ -399,12 +398,11 @@ Field3D generate_f1(const Mesh& mesh) { //make the gradients zero at both x-boundaries result(jx, jy, jz) = 0. + exp(-((60. * pow(x - p, 2)) + 1. - cos(2. * PI * (z - q)))) - - 50. - * (2. * p * exp(-60. * pow(-p, 2)) * x - + (-p * exp(-60. * pow(-p, 2)) - - (1 - p) * exp(-60. * pow(1 - p, 2))) - * pow(x, 2)) - * exp(-(1. - cos(2. * PI * (z - q)))); + - (50. + * (2. * p * exp(-60. * pow(-p, 2)) * x + + (-p * exp(-60. * pow(-p, 2)) - (1 - p) * exp(-60. * pow(1 - p, 2))) + * pow(x, 2)) + * exp(-(1. - cos(2. * PI * (z - q))))); } } } From 3ca3e689f2fedb8b4979998849904890258e6405 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Tue, 24 Feb 2026 15:01:51 +0000 Subject: [PATCH 12/12] Try `sor` preconditioner for 3D metrics in test-petsc-laplace --- tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx b/tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx index 902156f4c0..d4f375e1d7 100644 --- a/tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx +++ b/tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx @@ -186,7 +186,7 @@ int main(int argc, char** argv) { #endif // Preconditioner to use for 3D metrics - constexpr auto petsc_pc = petsc_has_mumps ? "mumps" : "lu"; + constexpr auto petsc_pc = petsc_has_mumps ? "mumps" : "sor"; auto& options_2nd = Options::root()["petsc2nd"]; if constexpr (bout::build::use_metric_3d) {