diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index d865c72999..90c63e0233 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -47,9 +47,10 @@ jobs: -DBOUT_USE_PETSC=ON -DBOUT_USE_SLEPC=ON -DBOUT_USE_SUNDIALS=ON + -DBOUT_USE_HYPRE=OFF -DBOUT_USE_ADIOS2=ON -DBOUT_ENABLE_PYTHON=ON - -DADIOS2_ROOT=/home/runner/local/adios2 + -DADIOS2_ROOT=/home/runner/local -DSUNDIALS_ROOT=/home/runner/local -DPETSC_DIR=/home/runner/local/petsc -DSLEPC_DIR=/home/runner/local/slepc" @@ -74,6 +75,7 @@ jobs: -DBOUT_USE_PETSC=ON -DBOUT_USE_SLEPC=ON -DBOUT_USE_SUNDIALS=ON + -DBOUT_USE_HYPRE=ON -DSUNDIALS_ROOT=/home/runner/local" on_cron: false @@ -86,6 +88,7 @@ jobs: -DBOUT_USE_PETSC=ON -DBOUT_USE_SLEPC=ON -DBOUT_USE_SUNDIALS=ON + -DBOUT_USE_HYPRE=ON -DSUNDIALS_ROOT=/home/runner/local" on_cron: false @@ -97,6 +100,7 @@ jobs: -DBOUT_USE_PETSC=ON -DBOUT_USE_SLEPC=ON -DBOUT_USE_SUNDIALS=ON + -DBOUT_USE_HYPRE=ON -DBOUT_BUILD_DOCS=OFF -DSUNDIALS_ROOT=/home/runner/local" omp_num_threads: 2 @@ -110,6 +114,7 @@ jobs: -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" omp_num_threads: 2 @@ -123,6 +128,7 @@ jobs: -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 @@ -130,25 +136,7 @@ jobs: build_petsc: -petsc on_cron: false - - name: "Coverage" - os: ubuntu-latest - cmake_options: "-DBUILD_SHARED_LIBS=ON - -DCMAKE_BUILD_TYPE=Debug - -DCHECK=3 - -DENABLE_COVERAGE=ON - -DBOUT_USE_PETSC=ON - -DBOUT_USE_SLEPC=ON - -DBOUT_USE_HDF5=ON - -DBOUT_USE_SUNDIALS=ON - -DBOUT_ENABLE_PYTHON=ON - -DSUNDIALS_ROOT=/home/runner/local" - unit_only: YES - on_cron: false exclude: - # Don't run the coverage tests if the branch isn't master or next - - is_master_or_next: false - config: - name: "Coverage" - is_cron: true config: on_cron: false @@ -178,6 +166,7 @@ jobs: slepc-dev liblapack-dev libparpack2-dev + libhypre-dev - uses: actions/checkout@v4 with: @@ -210,22 +199,6 @@ jobs: - name: Build BOUT++ run: UNIT_ONLY=${{ matrix.config.unit_only }} ./.ci_with_cmake.sh ${{ matrix.config.cmake_options }} - - name: Capture coverage - if: ${{ matrix.config.name == 'Coverage' }} - # Explicitly run the coverage capture target, because - # ci_script.sh also does the upload, and we're going to do - # that ourselves in the next step - run: | - # Ensure that there is a corresponding .gcda file for every .gcno file - # This is to try and make the coverage report slightly more accurate - # It still won't include, e.g. any solvers we don't build with though - find . -name "*.gcno" -exec sh -c 'touch -a "${1%.gcno}.gcda"' _ {} \; - make -C build code-coverage-capture - - - name: Upload coverage - if: ${{ matrix.config.name == 'Coverage' }} - uses: codecov/codecov-action@v4 - Fedora: # This is its own job as it doesn't use most of the steps of the # standard_tests diff --git a/bout++Config.cmake.in b/bout++Config.cmake.in index 5af0dc43ea..d461e58e4b 100644 --- a/bout++Config.cmake.in +++ b/bout++Config.cmake.in @@ -91,7 +91,7 @@ elseif(EXISTS "@NCXX4_CONFIG@") set(NCXX4_CONFIG "@NCXX4_CONFIG@") elseif(EXISTS "@NCXX_BINARY_DIR@") # If we downloaded netCDF-cxx4, then we need to add its build directory to our search paths - list(APPEND CMAKE_PREFIX_PATH "@NCXX_BINARY_DIR@") + list(APPEND CMAKE_PREFIX_PATH "@NCXX_BINARY_DIR@") endif() if(EXISTS "@PVODE_ROOT@") set(PVODE_ROOT "@PVODE_ROOT@") @@ -103,6 +103,12 @@ endif() if(EXISTS "@Libuuid_ROOT@") set(Libuuid_ROOT "@Libuuid_ROOT@") endif() +if(EXISTS "@ADIOS2_ROOT@") + set(ADIOS2_ROOT "@ADIOS2_ROOT@") +elseif(EXISTS "@ADIOS2_BINARY_DIR@") + # If we downloaded ADIOS2, then we need to add its build directory to our search paths + list(APPEND CMAKE_PREFIX_PATH "@ADIOS2_BINARY_DIR@") +endif() if(@BOUT_USE_SYSTEM_MPARK_VARIANT@) set(mpark_variant_ROOT "@mpark_variant_ROOT@") @@ -117,6 +123,11 @@ endif() set(MPIEXEC_EXECUTABLE @MPIEXEC_EXECUTABLE@) find_dependency(MPI @MPI_CXX_VERSION@ EXACT) +if (BOUT_HAS_ADIOS2 OR BOUT_HAS_HYPRE OR BOUT_HAS_SUNDIALS) + # These libraries may also require MPI_C + enable_language(C) + find_dependency(MPI COMPONENTS C) +endif() if (BOUT_USE_OPENMP) find_dependency(OpenMP) @@ -158,5 +169,8 @@ endif() if (BOUT_USE_UUID_SYSTEM_GENERATOR) find_dependency(Libuuid) endif() +if (BOUT_HAS_ADIOS2) + find_dependency(ADIOS2 @ADIOS2_VERSION@ EXACT) +endif() include("${CMAKE_CURRENT_LIST_DIR}/bout++Targets.cmake") diff --git a/cmake/FindHYPRE.cmake b/cmake/FindHYPRE.cmake index 2412981c80..1b9a5ca6f9 100644 --- a/cmake/FindHYPRE.cmake +++ b/cmake/FindHYPRE.cmake @@ -7,6 +7,7 @@ include(FindPackageHandleStandardArgs) find_package(HYPRE CONFIG QUIET) if (HYPRE_FOUND) + message(STATUS "Found HYPRE: ${HYPRE_VERSION}") return() endif() @@ -14,7 +15,7 @@ find_path(HYPRE_INCLUDE_DIR NAMES HYPRE.h DOC "HYPRE include directories" REQUIRED - PATH_SUFFIXES include + PATH_SUFFIXES include include/hypre ) find_library(HYPRE_LIBRARY diff --git a/cmake/SetupBOUTThirdParty.cmake b/cmake/SetupBOUTThirdParty.cmake index 9c49fe6fdc..b7c222eb72 100644 --- a/cmake/SetupBOUTThirdParty.cmake +++ b/cmake/SetupBOUTThirdParty.cmake @@ -190,9 +190,12 @@ endif() message(STATUS "NetCDF support: ${BOUT_USE_NETCDF}") set(BOUT_HAS_NETCDF ${BOUT_USE_NETCDF}) -option(BOUT_USE_ADIOS2 "Enable support for ADIOS output" ON) +option(BOUT_USE_ADIOS2 "Enable support for ADIOS output" OFF) option(BOUT_DOWNLOAD_ADIOS2 "Download and build ADIOS2" OFF) if (BOUT_USE_ADIOS2) + enable_language(C) + find_package(MPI REQUIRED COMPONENTS C) + if (BOUT_DOWNLOAD_ADIOS2) message(STATUS "Downloading and configuring ADIOS2") include(FetchContent) @@ -211,18 +214,11 @@ if (BOUT_USE_ADIOS2) # Note: SST requires but doesn't check at configure time set(ADIOS2_USE_SST OFF CACHE BOOL "" FORCE) FetchContent_MakeAvailable(adios2) - target_link_libraries(bout++ PUBLIC adios2::cxx11_mpi) message(STATUS "ADIOS2 done configuring") else() - find_package(ADIOS2) - if (ADIOS2_FOUND) - ENABLE_LANGUAGE(C) - find_package(MPI REQUIRED COMPONENTS C) - target_link_libraries(bout++ PUBLIC adios2::cxx11_mpi MPI::MPI_C) - else() - set(BOUT_USE_ADIOS2 OFF) - endif() + find_package(ADIOS2 REQUIRED) endif() + target_link_libraries(bout++ PUBLIC adios2::cxx11_mpi MPI::MPI_C) endif() message(STATUS "ADIOS2 support: ${BOUT_USE_ADIOS2}") set(BOUT_HAS_ADIOS2 ${BOUT_USE_ADIOS2}) diff --git a/examples/elm-pb/data-hypre/BOUT.inp b/examples/elm-pb/data-hypre/BOUT.inp index 4d9aa9481f..45dd59b536 100644 --- a/examples/elm-pb/data-hypre/BOUT.inp +++ b/examples/elm-pb/data-hypre/BOUT.inp @@ -1,26 +1,18 @@ -wall_limit = 1.55 # wall time limit (in hours) - zperiod = 15 # Fraction of a torus to simulate MZ = 16 # Number of points in Z -dump_format = "nc" # Dump file format. "nc" = NetCDF, "pdb" = PDB -restart_format = "nc" # Restart file format - [mesh] - staggergrids = false # Use staggered grids file = "cbm18_dens8.grid_nx68ny64.nc" # Grid file [mesh:paralleltransform] - -type = shifted # Use shifted metric method +type = shiftedinterp ################################################## # derivative methods [mesh:ddx] - first = C4 # order of first x derivatives second = C4 # order of second x derivatives upwind = W3 # order of upwinding method W3 = Weno3 @@ -36,20 +28,6 @@ first = C4 # Z derivatives can be done using FFT second = C4 upwind = W3 -[output] -shiftoutput = true # Put the output into field-aligned coordinates - -################################################## -# Laplacian inversion routines - -[laplace] -type = hypre3d - -flags = 0 # Flags for Laplacian inversion - -rtol = 1e-09 -atol = 1e-14 - ################################################## # FFTs @@ -63,23 +41,12 @@ fft_measurement_flag = measure # If using FFTW, perform tests to determine fast [solver] # mudq, mldq, mukeep, mlkeep preconditioner options -atol = 1e-08 # absolute tolerance -rtol = 1e-05 # relative tolerance +atol = 1.0e-8 # absolute tolerance +rtol = 1.0e-5 # relative tolerance use_precon = false # Use preconditioner: User-supplied or BBD -use_jacobian = false # Use user-supplied Jacobian mxstep = 5000 # Number of internal steps between outputs -adams_moulton = false # Use Adams-Moulton method (default is BDF) -func_iter = false # Functional iteration (default is Newton) -output_step = 1 # time between outputs -# settings file for BOUT++ -# High-Beta reduced MHD case - -################################################## -# Global settings used by the core code - -nout = 40 # number of time-steps ################################################## # settings for high-beta reduced MHD @@ -116,7 +83,7 @@ diamag_phi0 = true # Balance ExB against Vd for stationary equilibrium bm_exb_flag = 0 bm_mag_flag = 0 ################################################################## -withflow = false # With flow or not +withflow = false # With flow or not D_0 = 1.3e5 # differential potential D_s = 20 # shear parameter K_H_term = false # Contain K-H term @@ -125,8 +92,8 @@ x0 = 0.855 # peak location D_min = 3000 # constant ################################################################## -eHall = false # Include electron pressue effects in Ohm's law? -AA = 2.0 # ion mass in units of proton mass +eHall = false # Include electron pressure effects in Ohm's law? +AA = 2.0 # ion mass in units of proton mass noshear = false # zero all shear @@ -174,14 +141,14 @@ damp_t_const = 0.01 # Damping time constant ## Parallel pressure diffusion -diffusion_par = -1.0 # Parallel pressure diffusion (< 0 = none) +diffusion_par = -1.0 # Parallel pressure diffusion (< 0 = none) diffusion_p4 = -1e-05 # parallel hyper-viscous diffusion for pressure (< 0 = none) -diffusion_u4 = 1e-05 # parallel hyper-viscous diffusion for vorticity (< 0 = none) +diffusion_u4 = -1e-05 # parallel hyper-viscous diffusion for vorticity (< 0 = none) diffusion_a4 = -1e-05 # parallel hyper-viscous diffusion for vector potential (< 0 = none) ## heat source in pressure in watts -heating_P = -1 # heat power in watts (< 0 = none) +heating_P = -1 # heat power in watts (< 0 = none) hp_width = 0.1 # heat width, in percentage of nx (< 0 = none) hp_length = 0.3 # heat length in percentage of nx (< 0 = none) @@ -205,7 +172,7 @@ su_lengthr = 0.1 # right edge sink length in percentage of nx (< 0 = none) ## Viscosity and Hyper-viscosity viscos_par = -0.1 # Parallel viscosity (< 0 = none) -viscos_perp = -1.0 # Perpendicular +viscos_perp = -1.0 # Perpendicular viscosity (< 0 = none) hyperviscos = -1.0 # Radial hyper viscosity ## Compressional terms (only when compress = true) @@ -213,14 +180,11 @@ phi_curv = true # Include curvature*Grad(phi) in P equation # gamma = 1.6666 [phiSolver] -#inner_boundary_flags = 1 + 2 + 128 # INVERT_DC_GRAD + INVERT_AC_GRAD + INVERT_BNDRY_ONE -#outer_boundary_flags = 1 + 2 + 128 # INVERT_DC_GRAD + INVERT_AC_GRAD + INVERT_BNDRY_ONE -inner_boundary_flags = 1 + 4 # INVERT_DC_GRAD + INVERT_AC_LAP -outer_boundary_flags = 1 + 4 # INVERT_DC_GRAD + INVERT_AC_LAP +type = hypre3d +inner_boundary_flags = 0 # Dirichlet +outer_boundary_flags = 0 # Dirichlet [aparSolver] -#inner_boundary_flags = 1 + 2 + 128 # INVERT_DC_GRAD + INVERT_AC_GRAD + INVERT_BNDRY_ONE -#outer_boundary_flags = 1 + 2 + 128 # INVERT_DC_GRAD + INVERT_AC_GRAD + INVERT_BNDRY_ONE inner_boundary_flags = 1 + 4 # INVERT_DC_GRAD + INVERT_AC_LAP outer_boundary_flags = 1 + 4 # INVERT_DC_GRAD + INVERT_AC_LAP @@ -271,12 +235,9 @@ bndry_ydown = free_o3 # Zero gradient in the core bndry_core = neumann -[Vpar] - -bndry_core = neumann - [phi] bndry_xin = none bndry_xout = none -bndry_target = neumann +#bndry_target = neumann + diff --git a/examples/elm-pb/data/BOUT.inp b/examples/elm-pb/data/BOUT.inp index 8d1fedff25..a37d69a43a 100644 --- a/examples/elm-pb/data/BOUT.inp +++ b/examples/elm-pb/data/BOUT.inp @@ -94,7 +94,7 @@ diamag_phi0 = true # Balance ExB against Vd for stationary equilibrium bm_exb_flag = 0 bm_mag_flag = 0 ################################################################## -withflow = false # With flow or not +withflow = false # With flow or not D_0 = 1.3e5 # differential potential D_s = 20 # shear parameter K_H_term = false # Contain K-H term @@ -104,7 +104,7 @@ D_min = 3000 # constant ################################################################## eHall = false # Include electron pressure effects in Ohm's law? -AA = 2.0 # ion mass in units of proton mass +AA = 2.0 # ion mass in units of proton mass noshear = false # zero all shear @@ -152,14 +152,14 @@ damp_t_const = 0.01 # Damping time constant ## Parallel pressure diffusion -diffusion_par = -1.0 # Parallel pressure diffusion (< 0 = none) +diffusion_par = -1.0 # Parallel pressure diffusion (< 0 = none) diffusion_p4 = -1e-05 # parallel hyper-viscous diffusion for pressure (< 0 = none) -diffusion_u4 = -1e-05 # parallel hyper-viscous diffusion for vorticity (< 0 = none) +diffusion_u4 = -1e-05 # parallel hyper-viscous diffusion for vorticity (< 0 = none) diffusion_a4 = -1e-05 # parallel hyper-viscous diffusion for vector potential (< 0 = none) ## heat source in pressure in watts -heating_P = -1 # heat power in watts (< 0 = none) +heating_P = -1 # heat power in watts (< 0 = none) hp_width = 0.1 # heat width, in percentage of nx (< 0 = none) hp_length = 0.3 # heat length in percentage of nx (< 0 = none) diff --git a/examples/elm-pb/elm_pb.cxx b/examples/elm-pb/elm_pb.cxx index f108e58e2f..f830f3d98a 100644 --- a/examples/elm-pb/elm_pb.cxx +++ b/examples/elm-pb/elm_pb.cxx @@ -1142,6 +1142,9 @@ class ELMpb : public PhysicsModel { // Create a solver for the Laplacian phiSolver = Laplacian::create(&globalOptions["phiSolver"]); + // Save performance metrics to output, using the + // given name as the prefix. + phiSolver->savePerformance(*solver, "phiSolver"); aparSolver = Laplacian::create(&globalOptions["aparSolver"], loc); diff --git a/include/bout/hypre_interface.hxx b/include/bout/hypre_interface.hxx index 86791fbd3d..ae392de4f3 100644 --- a/include/bout/hypre_interface.hxx +++ b/include/bout/hypre_interface.hxx @@ -710,6 +710,15 @@ public: bool isAssembled() const { return assembled; } + /// Call HYPRE_IJMatrixPrint + void print() const { + if (!assembled) { + output_warn.write(""); + return; + } + HYPRE_IJMatrixPrint(*hypre_matrix, "hypreIJ.matrix"); + } + HypreMatrix yup(int index = 0) { return ynext(index + 1); } HypreMatrix ydown(int index = 0) { return ynext(-index - 1); } HypreMatrix ynext(int dir) { @@ -894,20 +903,30 @@ public: void setMaxIter(int max_iter) { checkHypreError(solverSetMaxIter(solver, max_iter)); } - double getFinalRelResNorm() { + double getFinalRelResNorm() const { HYPRE_Real resnorm{}; checkHypreError(solverGetFinalRelativeResidualNorm(solver, &resnorm)); return resnorm; } - int getNumItersTaken() { + /// Return the number of solver iterations taken + int getNumItersTaken() const { HYPRE_Int iters{}; checkHypreError(solverGetNumIterations(solver, &iters)); return iters; } + /// Return the number of BoomerAMG preconditioner iterations taken + int getNumItersTakenAMG() const { + HYPRE_Int iters{}; + checkHypreError(HYPRE_BoomerAMGGetNumIterations(precon, &iters)); + return iters; + } + + /// Set the HypreMatrix to be used in the solver void setMatrix(HypreMatrix* A_) { A = A_; } + /// Enable BoomerAMG preconditioner with the given HypreMatrix int setupAMG(HypreMatrix* P_) { CALI_CXX_MARK_FUNCTION; diff --git a/include/bout/invert/laplacexy2_hypre.hxx b/include/bout/invert/laplacexy2_hypre.hxx index 9a488ed6ff..1df0988d06 100644 --- a/include/bout/invert/laplacexy2_hypre.hxx +++ b/include/bout/invert/laplacexy2_hypre.hxx @@ -118,7 +118,6 @@ private: // Boundary conditions bool x_inner_dirichlet; // Dirichlet on inner X boundary? - bool x_outer_dirichlet; // Dirichlet on outer X boundary? bool y_bndry_dirichlet; // Dirichlet on Y boundary? bool print_timing; diff --git a/include/bout/paralleltransform.hxx b/include/bout/paralleltransform.hxx index 231fc5220c..49dea67743 100644 --- a/include/bout/paralleltransform.hxx +++ b/include/bout/paralleltransform.hxx @@ -216,7 +216,7 @@ public: std::vector getWeightsForYApproximation(int UNUSED(i), int UNUSED(j), int UNUSED(k), int UNUSED(yoffset)) override { - throw BoutException("ParallelTransform::getWeightsForYApproximation not implemented" + throw BoutException("ParallelTransform::getWeightsForYApproximation not implemented " "for `type = shifted`. Try `type = shiftedinterp`"); } diff --git a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx index 8de1a54099..c50be1db85 100644 --- a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx +++ b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx @@ -3,9 +3,9 @@ * Using Hypre Solvers * ************************************************************************** - * Copyright 2021 J.Omotani, C.MacMackin + * Copyright 2021 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -42,15 +42,14 @@ #include -LaplaceHypre3d::LaplaceHypre3d(Options* opt, const CELL_LOC loc, Mesh* mesh_in, - Solver* solver) +LaplaceHypre3d::LaplaceHypre3d(Options* opt, const CELL_LOC loc, Mesh* mesh_in, Solver*) : Laplacian(opt, loc, mesh_in), A(0.0), C1(1.0), C2(1.0), D(1.0), Ex(0.0), Ez(0.0), opts(opt == nullptr ? Options::getRoot()->getSection("laplace") : opt), lowerY(localmesh->iterateBndryLowerY()), upperY(localmesh->iterateBndryUpperY()), indexer(std::make_shared>( localmesh, getStencil(localmesh, lowerY, upperY))), operator3D(indexer), solution(indexer), rhs(indexer), - linearSystem(*localmesh, *opts), monitor(*this) { + linearSystem(*localmesh, *opts) { // Provide basic initialisation of field coefficients, etc. // Get relevent options from user input A.setLocation(location); @@ -60,7 +59,13 @@ LaplaceHypre3d::LaplaceHypre3d(Options* opt, const CELL_LOC loc, Mesh* mesh_in, Ex.setLocation(location); Ez.setLocation(location); + // Print matrix coefficients? + print_matrix = (*opts)["print_matrix"] + .doc("Print matrix coefficients when the operator is updated?") + .withDefault(false); + // Initialise Hypre objects + // Note: linearSystem is a bout::HypreSystem object linearSystem.setMatrix(&operator3D); linearSystem.setRHSVector(&rhs); linearSystem.setSolutionVector(&solution); @@ -71,23 +76,23 @@ LaplaceHypre3d::LaplaceHypre3d(Options* opt, const CELL_LOC loc, Mesh* mesh_in, // Checking flags are set to something which is not implemented // This is done binary (which is possible as each flag is a power of 2) - if (global_flags & ~implemented_flags) { + if (isGlobalFlagSet(~implemented_flags)) { throw BoutException("Attempted to set Laplacian inversion flag that is not " "implemented in LaplaceHypre3d"); } - if (inner_boundary_flags & ~implemented_boundary_flags) { + if (isInnerBoundaryFlagSet(~implemented_boundary_flags)) { throw BoutException("Attempted to set Laplacian inversion boundary flag that is not " "implemented in LaplaceHypre3d"); } - if (outer_boundary_flags & ~implemented_boundary_flags) { + if (isOuterBoundaryFlagSet(~implemented_boundary_flags)) { throw BoutException("Attempted to set Laplacian inversion boundary flag that is not " "implemented in LaplaceHypre3d"); } - if (lower_boundary_flags & ~implemented_boundary_flags) { + if ((lower_boundary_flags & ~implemented_boundary_flags) != 0) { throw BoutException("Attempted to set Laplacian inversion boundary flag that is not " "implemented in LaplaceHypre3d"); } - if (upper_boundary_flags & ~implemented_boundary_flags) { + if ((upper_boundary_flags & ~implemented_boundary_flags) != 0) { throw BoutException("Attempted to set Laplacian inversion boundary flag that is not " "implemented in LaplaceHypre3d"); } @@ -123,7 +128,7 @@ LaplaceHypre3d::LaplaceHypre3d(Options* opt, const CELL_LOC loc, Mesh* mesh_in, } BOUT_FOR_SERIAL(i, indexer->getRegionLowerY()) { - if (lower_boundary_flags & INVERT_AC_GRAD) { + if ((lower_boundary_flags & INVERT_AC_GRAD) != 0) { // Neumann on lower Y boundary operator3D(i, i) = -1. / coords->dy[i] / sqrt(coords->g_22[i]); operator3D(i, i.yp()) = 1. / coords->dy[i] / sqrt(coords->g_22[i]); @@ -135,7 +140,7 @@ LaplaceHypre3d::LaplaceHypre3d(Options* opt, const CELL_LOC loc, Mesh* mesh_in, } BOUT_FOR_SERIAL(i, indexer->getRegionUpperY()) { - if (upper_boundary_flags & INVERT_AC_GRAD) { + if ((upper_boundary_flags & INVERT_AC_GRAD) != 0) { // Neumann on upper Y boundary operator3D(i, i) = 1. / coords->dy[i] / sqrt(coords->g_22[i]); operator3D(i, i.ym()) = -1. / coords->dy[i] / sqrt(coords->g_22[i]); @@ -145,15 +150,6 @@ LaplaceHypre3d::LaplaceHypre3d(Options* opt, const CELL_LOC loc, Mesh* mesh_in, operator3D(i, i.ym()) = 0.5; } } - - // FIXME: This needs to be converted to outputVars - if (solver == nullptr) { - output_warn << "Warning: Need to pass a Solver to " - "Laplacian::create() to get iteration counts in the output." - << endl; - } else { - solver->addMonitor(&monitor); - } } LaplaceHypre3d::~LaplaceHypre3d() {} @@ -200,9 +196,9 @@ Field3D LaplaceHypre3d::solve(const Field3D& b_in, const Field3D& x0) { } BOUT_FOR_SERIAL(i, indexer->getRegionLowerY()) { - const BoutReal val = (lower_boundary_flags & INVERT_SET) ? x0[i] : 0.; + const BoutReal val = ((lower_boundary_flags & INVERT_SET) != 0) ? x0[i] : 0.; ASSERT1(std::isfinite(val)); - if (!(lower_boundary_flags & INVERT_RHS)) { + if ((lower_boundary_flags & INVERT_RHS) == 0) { b[i] = val; } else { ASSERT1(std::isfinite(b[i])); @@ -210,9 +206,9 @@ Field3D LaplaceHypre3d::solve(const Field3D& b_in, const Field3D& x0) { } BOUT_FOR_SERIAL(i, indexer->getRegionUpperY()) { - const BoutReal val = (upper_boundary_flags & INVERT_SET) ? x0[i] : 0.; + const BoutReal val = ((upper_boundary_flags & INVERT_SET) != 0) ? x0[i] : 0.; ASSERT1(std::isfinite(val)); - if (!(upper_boundary_flags & INVERT_RHS)) { + if ((upper_boundary_flags & INVERT_RHS) == 0) { b[i] = val; } else { ASSERT1(std::isfinite(b[i])); @@ -239,6 +235,7 @@ Field3D LaplaceHypre3d::solve(const Field3D& b_in, const Field3D& x0) { // Increment counters n_solves++; cumulative_iterations += linearSystem.getNumItersTaken(); + cumulative_amg_iterations += linearSystem.getNumItersTakenAMG(); CALI_MARK_END("LaplaceHypre3d_solve:solve"); @@ -246,6 +243,8 @@ Field3D LaplaceHypre3d::solve(const Field3D& b_in, const Field3D& x0) { // Create field from solution Field3D result = solution.toField(); + + // Fill guard and boundary cells localmesh->communicate(result); if (result.hasParallelSlices()) { BOUT_FOR(i, indexer->getRegionLowerY()) { result.ydown()[i] = result[i]; } @@ -277,7 +276,7 @@ void LaplaceHypre3d::updateMatrix3D() { const Field3D dc_dx = issetC ? DDX(C2) : Field3D(); const Field3D dc_dy = issetC ? DDY(C2) : Field3D(); const Field3D dc_dz = issetC ? DDZ(C2) : Field3D(); - const Field2D dJ_dy = DDY(coords->J / coords->g_22); + const auto dJ_dy = DDY(coords->J / coords->g_22); // Set up the matrix for the internal points on the grid. // Boundary conditions were set in the constructor. @@ -286,7 +285,8 @@ void LaplaceHypre3d::updateMatrix3D() { // avoid confusing it with the x-index. // Calculate coefficients for the terms in the differential operator - BoutReal C_df_dx = coords->G1[l], C_df_dz = coords->G3[l]; + BoutReal C_df_dx = coords->G1[l]; + BoutReal C_df_dz = coords->G3[l]; if (issetD) { C_df_dx *= D[l]; C_df_dz *= D[l]; @@ -304,9 +304,9 @@ void LaplaceHypre3d::updateMatrix3D() { C_df_dz += Ez[l]; } - BoutReal C_d2f_dx2 = coords->g11[l], - C_d2f_dy2 = (coords->g22[l] - 1.0 / coords->g_22[l]), - C_d2f_dz2 = coords->g33[l]; + BoutReal C_d2f_dx2 = coords->g11[l]; + BoutReal C_d2f_dy2 = (coords->g22[l] - 1.0 / coords->g_22[l]); + BoutReal C_d2f_dz2 = coords->g33[l]; if (issetD) { C_d2f_dx2 *= D[l]; C_d2f_dy2 *= D[l]; @@ -343,8 +343,8 @@ void LaplaceHypre3d::updateMatrix3D() { // The values stored in the y-boundary are already interpolated // up/down, so we don't want the matrix to do any such // interpolation there. - const int yup = (l.y() == localmesh->yend && upperY.intersects(l.x())) ? -1 : 0, - ydown = (l.y() == localmesh->ystart && lowerY.intersects(l.x())) ? -1 : 0; + const int yup = (l.y() == localmesh->yend && upperY.intersects(l.x())) ? -1 : 0; + const int ydown = (l.y() == localmesh->ystart && lowerY.intersects(l.x())) ? -1 : 0; operator3D.yup(yup)(l, l.yp()) = 0.0; operator3D.ydown(ydown)(l, l.ym()) = 0.0; operator3D.yup(yup)(l, l.xp().yp()) = 0.0; @@ -376,7 +376,8 @@ void LaplaceHypre3d::updateMatrix3D() { C_d2f_dy2 *= D[l]; } - BoutReal C_d2f_dxdy = 2 * coords->g12[l], C_d2f_dydz = 2 * coords->g23[l]; + BoutReal C_d2f_dxdy = 2 * coords->g12[l]; + BoutReal C_d2f_dydz = 2 * coords->g23[l]; if (issetD) { C_d2f_dxdy *= D[l]; C_d2f_dydz *= D[l]; @@ -396,8 +397,8 @@ void LaplaceHypre3d::updateMatrix3D() { // The values stored in the y-boundary are already interpolated // up/down, so we don't want the matrix to do any such // interpolation there. - const int yup = (l.y() == localmesh->yend && upperY.intersects(l.x())) ? -1 : 0, - ydown = (l.y() == localmesh->ystart && lowerY.intersects(l.x())) ? -1 : 0; + const int yup = (l.y() == localmesh->yend && upperY.intersects(l.x())) ? -1 : 0; + const int ydown = (l.y() == localmesh->ystart && lowerY.intersects(l.x())) ? -1 : 0; operator3D.yup(yup)(l, l.yp()) += C_df_dy + C_d2f_dy2; operator3D.ydown(ydown)(l, l.ym()) += -C_df_dy + C_d2f_dy2; @@ -411,6 +412,11 @@ void LaplaceHypre3d::updateMatrix3D() { operator3D.ydown(ydown)(l, l.ym().zm()) += C_d2f_dydz; } operator3D.assemble(); + + if (print_matrix) { + operator3D.print(); + } + linearSystem.setupAMG(&operator3D); updateRequired = false; @@ -437,19 +443,15 @@ OperatorStencil LaplaceHypre3d::getStencil(Mesh* localmesh, OffsetInd3D zero; // Add interior cells - const std::vector interpolatedUpElements = {zero.yp(), zero.xp().yp(), - zero.xm().yp(), zero.yp().zp(), - zero.yp().zm()}, - interpolatedDownElements = { - zero.ym(), zero.xp().ym(), zero.xm().ym(), - zero.ym().zp(), zero.ym().zm()}; - std::set interiorStencil = {zero, zero.xp(), - zero.xm(), zero.zp(), - zero.zm(), zero.xp().zp(), - zero.xp().zm(), zero.xm().zp(), - zero.xm().zm()}, - lowerEdgeStencil = interiorStencil, - upperEdgeStencil = interiorStencil; + const std::vector interpolatedUpElements = { + zero.yp(), zero.xp().yp(), zero.xm().yp(), zero.yp().zp(), zero.yp().zm()}; + const std::vector interpolatedDownElements = { + zero.ym(), zero.xp().ym(), zero.xm().ym(), zero.ym().zp(), zero.ym().zm()}; + std::set interiorStencil = { + zero, zero.xp(), zero.xm(), zero.zp(), zero.zm(), + zero.xp().zp(), zero.xp().zm(), zero.xm().zp(), zero.xm().zm()}; + std::set lowerEdgeStencil = interiorStencil; + std::set upperEdgeStencil = interiorStencil; for (const auto& i : interpolatedDownElements) { for (auto& j : interpPattern) { @@ -466,9 +468,11 @@ OperatorStencil LaplaceHypre3d::getStencil(Mesh* localmesh, upperEdgeStencil.insert(i); } const std::vector interiorStencilVector(interiorStencil.begin(), - interiorStencil.end()), - lowerEdgeStencilVector(lowerEdgeStencil.begin(), lowerEdgeStencil.end()), - upperEdgeStencilVector(upperEdgeStencil.begin(), upperEdgeStencil.end()); + interiorStencil.end()); + const std::vector lowerEdgeStencilVector(lowerEdgeStencil.begin(), + lowerEdgeStencil.end()); + const std::vector upperEdgeStencilVector(upperEdgeStencil.begin(), + upperEdgeStencil.end()); // If there is a lower y-boundary then create a part of the stencil // for cells immediately adjacent to it. @@ -536,19 +540,35 @@ OperatorStencil LaplaceHypre3d::getStencil(Mesh* localmesh, return stencil; } -int LaplaceHypre3d::Hypre3dMonitor::call(Solver*, BoutReal, int, int) { - if (laplace.n_solves == 0) { - laplace.average_iterations = 0.0; - return 0; - } +void LaplaceHypre3d::outputVars(Options& output_options, + const std::string& time_dimension) const { + BoutReal mean_iterations = 0.0; + BoutReal mean_amg_iterations = 0.0; + BoutReal rel_res_norm = linearSystem.getFinalRelResNorm(); - // Calculate average and reset counters - laplace.average_iterations = static_cast(laplace.cumulative_iterations) - / static_cast(laplace.n_solves); + if (n_solves > 0) { + // Calculate average + mean_iterations = + static_cast(cumulative_iterations) / static_cast(n_solves); - output_info.write("\nHypre3d average iterations: {}\n", laplace.average_iterations); + mean_amg_iterations = static_cast(cumulative_amg_iterations) + / static_cast(n_solves); + } - return 0; + std::string name = getPerformanceName(); + + output_options[fmt::format("{}_mean_its", name)] + .assignRepeat(mean_iterations, time_dimension) + .setAttributes({{"source", "hypre3d_laplace"}, + {"description", "Mean number of solver iterations"}}); + output_options[fmt::format("{}_mean_amg_its", name)] + .assignRepeat(mean_amg_iterations, time_dimension) + .setAttributes({{"source", "hypre3d_laplace"}, + {"description", "Mean number of BoomerAMG iterations"}}); + output_options[fmt::format("{}_rel_res_norm", name)] + .assignRepeat(rel_res_norm, time_dimension) + .setAttributes({{"source", "hypre3d_laplace"}, + {"description", "Final relative residual norm"}}); } #endif // BOUT_HAS_HYPRE diff --git a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.hxx b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.hxx index f3ca44df46..d58ef6f688 100644 --- a/src/invert/laplace/impls/hypre3d/hypre3d_laplace.hxx +++ b/src/invert/laplace/impls/hypre3d/hypre3d_laplace.hxx @@ -6,9 +6,9 @@ *ex\nabla_x x + ez\nabla_z x + a x = b\f$ * ************************************************************************** - * Copyright 2021 J. Omotani, C. MacMackin + * Copyright 2021 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -201,28 +201,36 @@ public: bout::HypreVector rhs; bout::HypreSystem linearSystem; - // used to save some statistics + // used to save Hypre solver statistics int n_solves = 0; int cumulative_iterations = 0; - BoutReal average_iterations = 0.0; - class Hypre3dMonitor : public Monitor { - public: - Hypre3dMonitor(LaplaceHypre3d& laplace_in) : laplace(laplace_in) {} - int call(Solver*, BoutReal, int, int) override; - - private: - LaplaceHypre3d& laplace; - }; - Hypre3dMonitor monitor; - - bool use_precon; // Switch for preconditioning - bool rightprec; // Right preconditioning + int cumulative_amg_iterations = 0; + + /// Save performance metrics + /// This is called by LaplacianMonitor, that is enabled + /// if Laplacian::savePerformance(Solver&, string&) is called. + /// + /// # Example + /// + /// // Create a Laplacian solver using "phiSolver" options + /// phiSolver = Laplacian::create(&Options::root()["phiSolver"]); + /// // Enable output diagnostics with "phiSolver" prefix + /// phiSolver->savePerformance(solver, "phiSolver"); + /// + /// Saves the following variables to the output: + /// - phiSolver_mean_its Mean number of solver iterations taken + /// - phiSolver_mean_amg_its Mean number of BoomerAMG preconditioner iterations + /// - phiSolver_rel_res_norm The final solver relative residual norm + void outputVars(Options& output_options, + const std::string& time_dimension) const override; // These are the implemented flags - static constexpr int implemented_flags = INVERT_START_NEW, - implemented_boundary_flags = - INVERT_AC_GRAD + INVERT_SET + INVERT_RHS; + static constexpr int implemented_flags = INVERT_START_NEW; + static constexpr int implemented_boundary_flags = + INVERT_AC_GRAD + INVERT_SET + INVERT_RHS; + + bool print_matrix; ///< Print matrix coefficients }; #endif // BOUT_HAS_HYPRE diff --git a/src/invert/laplacexy2/laplacexy2.cxx b/src/invert/laplacexy2/laplacexy2.cxx index 8b5e98d747..713fd8e68f 100644 --- a/src/invert/laplacexy2/laplacexy2.cxx +++ b/src/invert/laplacexy2/laplacexy2.cxx @@ -14,10 +14,12 @@ #include +namespace { Ind2D index2d(Mesh* mesh, int x, int y) { int ny = mesh->LocalNy; return Ind2D(x * ny + y, ny, 1); } +} // namespace LaplaceXY2::LaplaceXY2(Mesh* m, Options* opt, const CELL_LOC loc) : localmesh(m == nullptr ? bout::globals::mesh : m), diff --git a/src/invert/laplacexy2/laplacexy2_hypre.cxx b/src/invert/laplacexy2/laplacexy2_hypre.cxx index b18903be3a..a5028169f8 100644 --- a/src/invert/laplacexy2/laplacexy2_hypre.cxx +++ b/src/invert/laplacexy2/laplacexy2_hypre.cxx @@ -8,7 +8,6 @@ #include "bout/mesh.hxx" #include "bout/output.hxx" #include "bout/sys/timer.hxx" -#include "bout/utils.hxx" #include @@ -25,11 +24,6 @@ inline void gpuAssert(cudaError_t code, const char* file, int line, bool abort = } #endif -Ind2D index2d(Mesh* mesh, int x, int y) { - int ny = mesh->LocalNy; - return Ind2D(x * ny + y, ny, 1); -} - LaplaceXY2Hypre::LaplaceXY2Hypre(Mesh* m, Options* opt, const CELL_LOC loc) : localmesh(m == nullptr ? bout::globals::mesh : m), indexConverter(std::make_shared>( diff --git a/src/mesh/boundary_factory.cxx b/src/mesh/boundary_factory.cxx index 00282566a9..d112a216ad 100644 --- a/src/mesh/boundary_factory.cxx +++ b/src/mesh/boundary_factory.cxx @@ -318,7 +318,7 @@ BoundaryOpBase* BoundaryFactory::createFromOptions(const string& varname, /// Then (all, all) if (region->isParallel) { // Different default for parallel boundary regions - varOpts->get(prefix + "par_all", set, "parallel_dirichlet"); + varOpts->get(prefix + "par_all", set, "parallel_dirichlet_o2"); } else { varOpts->get(prefix + "all", set, "dirichlet"); } diff --git a/src/sys/options/options_adios.cxx b/src/sys/options/options_adios.cxx index 529ad5b7f5..b3acbaada6 100644 --- a/src/sys/options/options_adios.cxx +++ b/src/sys/options/options_adios.cxx @@ -245,7 +245,7 @@ bool readAttribute(adios2::IO& io, const std::string& name, const std::string& t return false; } -Options OptionsADIOS::read() { +Options OptionsADIOS::read([[maybe_unused]] bool lazy) { Timer timer("io"); // Open file diff --git a/src/sys/options/options_adios.hxx b/src/sys/options/options_adios.hxx index 73170c9501..6e81508ca3 100644 --- a/src/sys/options/options_adios.hxx +++ b/src/sys/options/options_adios.hxx @@ -50,7 +50,7 @@ public: OptionsADIOS& operator=(OptionsADIOS&&) noexcept = default; /// Read options from file - Options read() override; + Options read(bool lazy = true) override; /// Write options to file void write(const Options& options, const std::string& time_dim) override; diff --git a/src/sys/options/options_netcdf.hxx b/src/sys/options/options_netcdf.hxx index 657068f5a8..09a2da3220 100644 --- a/src/sys/options/options_netcdf.hxx +++ b/src/sys/options/options_netcdf.hxx @@ -51,12 +51,12 @@ public: /// Write options to file void write(const Options& options) { write(options, "t"); } - void write(const Options& options, const std::string& time_dim); + void write(const Options& options, const std::string& time_dim) override; /// Check that all variables with the same time dimension have the /// same size in that dimension. Throws BoutException if there are /// any differences, otherwise is silent - void verifyTimesteps() const; + void verifyTimesteps() const override; private: enum class FileMode { diff --git a/tests/integrated/test-laplace-hypre3d/data_circular_core-sol/BOUT.inp b/tests/integrated/test-laplace-hypre3d/data_circular_core-sol/BOUT.inp index e16eb9752a..dd7f1ba9ce 100644 --- a/tests/integrated/test-laplace-hypre3d/data_circular_core-sol/BOUT.inp +++ b/tests/integrated/test-laplace-hypre3d/data_circular_core-sol/BOUT.inp @@ -62,8 +62,8 @@ g_23 = Bt*hthe*Rxy/Bp # zShift = \int_{theta0}^{theta} {nu dtheta} arctanarg = (r - 1)*tan(theta/2.)/sqrt(1. - r^2) -zshift = r^2*(r*sin(theta)/((r^2 - 1)*(r*cos(theta) - 1.)) - 2.*atan(arctanarg)/(1. - r^2)^1.5) -shiftangle = r^2 * 2.*pi/(1. - r^2)^1.5 +zShift = r^2*(r*sin(theta)/((r^2 - 1)*(r*cos(theta) - 1.)) - 2.*atan(arctanarg)/(1. - r^2)^1.5) +ShiftAngle = r^2 * 2.*pi/(1. - r^2)^1.5 [mesh:paralleltransform] type = shiftedinterp diff --git a/tests/integrated/test-laplace-hypre3d/data_circular_core/BOUT.inp b/tests/integrated/test-laplace-hypre3d/data_circular_core/BOUT.inp index 30fe3c0738..956718054f 100644 --- a/tests/integrated/test-laplace-hypre3d/data_circular_core/BOUT.inp +++ b/tests/integrated/test-laplace-hypre3d/data_circular_core/BOUT.inp @@ -59,8 +59,8 @@ g_23 = Bt*hthe*Rxy/Bp # zShift = \int_{theta0}^{theta} {nu dtheta} arctanarg = (r - 1)*tan(theta/2.)/sqrt(1. - r^2) -zshift = r^2*(r*sin(theta)/((r^2 - 1)*(r*cos(theta) - 1.)) - 2.*atan(arctanarg)/(1. - r^2)^1.5) -shiftangle = r^2 * 2.*pi/(1. - r^2)^1.5 +zShift = r^2*(r*sin(theta)/((r^2 - 1)*(r*cos(theta) - 1.)) - 2.*atan(arctanarg)/(1. - r^2)^1.5) +ShiftAngle = r^2 * 2.*pi/(1. - r^2)^1.5 [mesh:paralleltransform] type = shiftedinterp diff --git a/tests/integrated/test-laplacexy2-hypre/data/BOUT.inp b/tests/integrated/test-laplacexy2-hypre/data/BOUT.inp index f79045e4d4..3f99c56d0c 100644 --- a/tests/integrated/test-laplacexy2-hypre/data/BOUT.inp +++ b/tests/integrated/test-laplacexy2-hypre/data/BOUT.inp @@ -29,6 +29,7 @@ rtol = 1e-14 core_bndry_dirichlet = true pf_bndry_dirichlet = true y_bndry = dirichlet +hypre_solver_type = gmres [f] # make an input: diff --git a/tests/integrated/test-laplacexy2-hypre/runtest b/tests/integrated/test-laplacexy2-hypre/runtest index c35c3609da..55f5982a0c 100755 --- a/tests/integrated/test-laplacexy2-hypre/runtest +++ b/tests/integrated/test-laplacexy2-hypre/runtest @@ -7,10 +7,11 @@ # requires: hypre # cores: 8 -from boututils.run_wrapper import build_and_log, shell, shell_safe, launch, launch_safe -from boutdata.collect import collect from sys import exit +from boutdata.collect import collect +from boututils.run_wrapper import build_and_log, launch_safe, shell + tol = 5.0e-8 argslist = [ @@ -44,15 +45,14 @@ success = True for nproc in [8]: print(" %d processors...." % nproc) for args in argslist: - cmd = "test-laplacexy " + args + cmd = "./test-laplacexy " + args shell("rm data/BOUT.dmp.*.nc > /dev/null 2>&1") - s, out = launch(cmd, nproc=nproc, pipe=True, verbose=True) + s, out = launch_safe(cmd, nproc=nproc, pipe=True, verbose=True) - f = open("run.log." + str(nproc), "w") - f.write(out) - f.close() + with open(f"run.log.{nproc}", "w") as f: + f.write(out) # Collect output data error = collect("max_error", path="data", info=False) diff --git a/tests/integrated/test-laplacexy2-hypre/test-laplacexy.cxx b/tests/integrated/test-laplacexy2-hypre/test-laplacexy.cxx index 0ae2cbf312..2fb8a10777 100644 --- a/tests/integrated/test-laplacexy2-hypre/test-laplacexy.cxx +++ b/tests/integrated/test-laplacexy2-hypre/test-laplacexy.cxx @@ -63,7 +63,7 @@ int main(int argc, char** argv) { Field2D absolute_error = abs(f - sol); BoutReal max_error = max(absolute_error, true); - output << "Magnitude of maximum absolute error is " << max_error << endl; + output.write("Magnitude of maximum absolute error is {}\n", max_error); sol.getMesh()->communicate(sol); Field2D rhs_check = Laplace_perpXY(a, sol);