diff --git a/.build_petsc_for_ci.sh b/.build_petsc_for_ci.sh new file mode 100755 index 0000000000..ce69944c9c --- /dev/null +++ b/.build_petsc_for_ci.sh @@ -0,0 +1,60 @@ +#!/bin/bash + +set -e + +if test $BUILD_PETSC ; then + if [[ ! -d $HOME/local/petsc/include/petsc ]]; then + echo "****************************************" + echo "Building PETSc" + echo "****************************************" + + git clone -b release https://gitlab.com/petsc/petsc.git petsc --depth=1 + + unset PETSC_DIR + unset PETSC_ARCH + + pushd petsc + ./configure \ + --with-mpi=yes \ + --with-precision=double \ + --with-scalar-type=real \ + --with-shared-libraries=1 \ + --with-debugging=0 \ + {C,CXX,F}OPTFLAGS="-O3" \ + --prefix=$HOME/local/petsc + + make && make install + popd + + echo "****************************************" + echo " Finished building PETSc" + echo "****************************************" + + echo "****************************************" + echo "Building SLEPc" + echo "****************************************" + + git clone -b release https://gitlab.com/slepc/slepc.git slepc --depth=1 + + pushd slepc + unset SLEPC_DIR + unset SLEPC_ARCH + PETSC_DIR=$HOME/local/petsc ./configure --prefix=$HOME/local/slepc + + make SLEPC_DIR=$(pwd) PETSC_DIR=$HOME/local/petsc + make SLEPC_DIR=$(pwd) PETSC_DIR=$HOME/local/petsc install + popd + + echo "****************************************" + echo " Finished building SLEPc" + echo "****************************************" + else + echo "****************************************" + echo " PETSc already installed" + echo "****************************************" + fi +else + echo "****************************************" + echo " PETSc not requested" + echo "****************************************" +fi diff --git a/.ci_fedora.sh b/.ci_fedora.sh index 5106480dc5..1a3d2d92b2 100755 --- a/.ci_fedora.sh +++ b/.ci_fedora.sh @@ -37,7 +37,7 @@ then cat /etc/os-release # Ignore weak depencies echo "install_weak_deps=False" >> /etc/dnf/dnf.conf - time dnf -y install dnf-plugins-core python3-pip cmake + time dnf -y install dnf5-plugins python3-pip cmake # Allow to override packages - see #2073 time dnf copr enable -y davidsch/fixes4bout || : time dnf -y upgrade @@ -49,7 +49,7 @@ then sudo -u test ${0/\/tmp/\/home\/test} $mpi ## If we are called as normal user, run test else - pip install --user zoidberg + pip install --user zoidberg natsort . /etc/profile.d/modules.sh module load mpi/${1}-x86_64 export OMPI_MCA_rmaps_base_oversubscribe=yes diff --git a/.clang-tidy b/.clang-tidy index fbd5d2f27c..7a1c58af80 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -1,5 +1,5 @@ --- -Checks: 'clang-diagnostic-*,clang-analyzer-*,-*,performance-*,readability-*,bugprone-*,clang-analyzer-*,cppcoreguidelines-*,mpi-*,misc-*,-readability-magic-numbers,-cppcoreguidelines-avoid-magic-numbers,-misc-non-private-member-variables-in-classes,-cppcoreguidelines-pro-bounds-array-to-pointer-decay,-cppcoreguidelines-pro-type-vararg,-clang-analyzer-optin.mpi*,-bugprone-exception-escape,-cppcoreguidelines-pro-bounds-pointer-arithmetic' +Checks: 'clang-diagnostic-*,clang-analyzer-*,-*,performance-*,readability-*,bugprone-*,clang-analyzer-*,cppcoreguidelines-*,mpi-*,misc-*,-readability-magic-numbers,-cppcoreguidelines-avoid-magic-numbers,-misc-non-private-member-variables-in-classes,-cppcoreguidelines-pro-bounds-array-to-pointer-decay,-cppcoreguidelines-pro-type-vararg,-clang-analyzer-optin.mpi*,-bugprone-exception-escape,-cppcoreguidelines-pro-bounds-pointer-arithmetic,-readability-function-cognitive-complexity' WarningsAsErrors: '' HeaderFilterRegex: '' AnalyzeTemporaryDtors: false diff --git a/.docker/build.sh b/.docker/build.sh new file mode 100644 index 0000000000..12b56c767f --- /dev/null +++ b/.docker/build.sh @@ -0,0 +1,29 @@ +set -ex + +if which docker &> /dev/null ; then + cmd="time sudo docker" +else + cmd="time podman" +fi + +file=$1 +(test $file && test -f $file) || file=.docker/fedora/Dockerfile +test $# -gt 0 && shift + +if test -x $file +then + COMMIT=$(git rev-parse HEAD) $file $@ > Dockerfile +else + cp $file Dockerfile +fi + +$cmd login -p $DOCKER_TOKEN -u $DOCKER_USER + +cat Dockerfile + +$cmd build -t mobydick . --build-arg=COMMIT=$(git rev-parse HEAD) --build-arg=MPI=$MPI --build-arg=CMAKE_OPTIONS="$CMAKE_OPTIONS" --build-arg=BASE=$BASE +for tag in $TAGS +do + $cmd tag mobydick $tag + $cmd push $tag +done diff --git a/.docker/fedora/Dockerfile b/.docker/fedora/Dockerfile new file mode 100644 index 0000000000..56445f8cea --- /dev/null +++ b/.docker/fedora/Dockerfile @@ -0,0 +1,39 @@ +ARG BASE +FROM ghcr.io/dschwoerer/bout-container-base:$BASE + +# ---------------------------------------------------------------- +# Build and install BOUT++ +# ---------------------------------------------------------------- +# user: boutuser +# password: boutforever + +ARG MPI +ARG URL +ARG COMMIT +ARG CMAKE_OPTIONS + +RUN sudo ls + +# Checkout submodules now so configure later is fast, and iterating on +# it less painful +RUN git clone $URL \ + && chown -R boutuser /home/boutuser/BOUT-dev \ + && cd BOUT-dev \ + && git checkout $COMMIT \ + && git submodule update --init --recursive + + +WORKDIR /home/boutuser/BOUT-dev + +RUN cmake -S . -B build -DCMAKE_INSTALL_PREFIX=/opt/bout++/ \ + -DBOUT_GENERATE_FIELDOPS=OFF \ + -DBOUT_USE_PETSC=ON -DPETSc_ROOT=/opt/petsc \ + -DBOUT_ENABLE_PYTHON=ON \ + -DBOUT_USE_SUNDIALS=ON -DSUNDIALS_ROOT=/usr/lib64/$MPI/ -DSUNDIALS_INCLUDE_DIR=/usr/include/$MPI-x86_64/sundials/ \ + $CMAKE_OPTIONS || (cat /home/boutuser/BOUT-dev/build/CMakeFiles/CMake{Output,Error}.log ; exit 1) + + +RUN make -C build -j 2 +RUN sudo make -C build install + +RUN find /opt/bout++ diff --git a/.github/workflows/clang-tidy-review.yml b/.github/workflows/clang-tidy-review.yml index 86161f583d..24e87cd447 100644 --- a/.github/workflows/clang-tidy-review.yml +++ b/.github/workflows/clang-tidy-review.yml @@ -22,7 +22,7 @@ jobs: submodules: true - name: Run clang-tidy - uses: ZedThree/clang-tidy-review@v0.8.4 + uses: ZedThree/clang-tidy-review@v0.13.1 id: review with: build_dir: build @@ -44,3 +44,6 @@ jobs: -DBOUT_ENABLE_PYTHON=OFF \ -DCMAKE_EXPORT_COMPILE_COMMANDS=On \ -DBOUT_UPDATE_GIT_SUBMODULE=OFF + + - name: Upload clang-tidy fixes + uses: ZedThree/clang-tidy-review/upload@v0.13.1 diff --git a/.github/workflows/docker.yml b/.github/workflows/docker.yml new file mode 100644 index 0000000000..dda81cd8fc --- /dev/null +++ b/.github/workflows/docker.yml @@ -0,0 +1,84 @@ +name: Docker + +on: + push: + branches: + - master + - next + # Add your branch here if you want containers for it + - db-WIP + - docker-ci + +env: + REGISTRY: ghcr.io + IMAGE_NAME: ${{ github.repository }} + +jobs: + build: + name: Build ${{ matrix.config.name }} Container (${{ matrix.metric3d.name }} ; ${{matrix.mpi}}) + runs-on: ubuntu-latest + permissions: + contents: read + packages: write + + strategy: + fail-fast: false + matrix: + mpi: [mpich] + metric3d: + - name: "With OpenMP" + cmake: ON + base_prefix: "openmp-" + tag_prefix: "3d-" + - name: "Without OpenMP" + cmake: OFF + base_prefix: "" + tag_prefix: "" + config: + - name: "Debug" + tag_postfix: "debug" + options: "-DCHECK=3" + base_postfix: "mini" + + - name: "Optimised" + tag_postfix: "opt" + options: "-DCHECK=0" + base_postfix: "mini" + + - name: "Full" + tag_postfix: "full" + options: "-DCHECK=3" + base_postfix: "full" + + steps: + - name: Checkout repository + uses: actions/checkout@v3 + + - name: Log in to the Container registry + uses: docker/login-action@65b78e6e13532edd9afa3aa52ac7964289d1a9c1 + with: + registry: ${{ env.REGISTRY }} + username: ${{ github.actor }} + password: ${{ secrets.GITHUB_TOKEN }} + + - name: Extract metadata (tags, labels) for Docker + id: meta + uses: docker/metadata-action@9ec57ed1fcdbf14dcef7dfbe97b2010124a938b7 + with: + images: ${{ env.REGISTRY }}/${{ env.IMAGE_NAME }} + flavor: | + prefix=${{ matrix.mpi }}-${{matrix.metric3d.tag_prefix}}${{ matrix.config.tag_postfix }}- + + - name: Build and push Docker image + uses: docker/build-push-action@f2a1d5e99d037542a71f64918e516c093c6f3fc4 + with: + build-args: | + BASE=${{ matrix.mpi }}-${{ matrix.metric3d.base_prefix }}${{ matrix.config.base_postfix }}-main + MPI=${{ matrix.mpi }} + CMAKE_OPTIONS=${{ matrix.config.options }} -DBOUT_ENABLE_METRIC_3D=${{ matrix.metric3d.cmake }} -DBOUT_ENABLE_OPENMP=${{ matrix.metric3d.cmake }} + COMMIT=${{ github.sha }} + URL=${{ github.server_url }}/${{ github.repository }} + context: .docker/fedora/ + push: true + tags: ${{ steps.meta.outputs.tags }} + labels: ${{ steps.meta.outputs.labels }} diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index a4d293b802..da209d6595 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -37,9 +37,9 @@ jobs: - name: "Optimised, static" os: ubuntu-20.04 - cmake_options: "-DBOUT_ENABLE_CHECKS=OFF + cmake_options: "-DCHECK=0 -DBUILD_SHARED_LIBS=OFF - -DCXXFLAGS=-Ofast + -DCMAKE_CXX_FLAGS=-Ofast -DBOUT_ENABLE_SIGNAL=OFF -DBOUT_ENABLE_TRACK=OFF -DBOUT_ENABLE_BACKTRACE=OFF @@ -50,7 +50,7 @@ jobs: - name: "Debug, shared" os: ubuntu-20.04 - cmake_options: "-DBOUT_ENABLE_CHECKS=3 + cmake_options: "-DCHECK=3 -DCMAKE_BUILD_TYPE=Debug -DBOUT_ENABLE_SIGNAL=ON -DBOUT_ENABLE_TRACK=ON @@ -83,6 +83,21 @@ jobs: -DSUNDIALS_ROOT=/home/runner/local" omp_num_threads: 2 + - name: "CMake, new PETSc" + os: ubuntu-20.04 + 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_ENABLE_PYTHON=ON + -DSUNDIALS_ROOT=/home/runner/local + -DPETSC_DIR=/home/runner/local/petsc + -DSLEPC_DIR=/home/runner/local/slepc" + + build_petsc: -petsc + - name: "Coverage" os: ubuntu-20.04 cmake_options: "-DBUILD_SHARED_LIBS=ON @@ -134,7 +149,7 @@ jobs: - name: Install pip packages run: | - ./.pip_install_for_ci.sh 'cython~=0.29' 'netcdf4~=1.5' 'sympy~=1.5' 'gcovr' 'cmake' 'h5py' zoidberg fastcov + ./.pip_install_for_ci.sh 'cython~=0.29' 'netcdf4~=1.5' 'sympy~=1.5' 'gcovr' 'cmake' zoidberg fastcov # Add the pip install location to the runner's PATH echo ~/.local/bin >> $GITHUB_PATH @@ -142,12 +157,15 @@ jobs: uses: actions/cache@v2 with: path: /home/runner/local - key: bout-sundials-${{ matrix.config.os }} + key: bout-sundials-${{ matrix.config.os }}${{ matrix.config.build_petsc }} - name: Build SUNDIALS run: ./.build_sundials_for_ci.sh - - name: Build (CMake) + - name: Build PETSc + run: BUILD_PETSC=${{ matrix.config.build_petsc }} ./.build_petsc_for_ci.sh + + - name: Build BOUT++ run: UNIT_ONLY=${{ matrix.config.unit_only }} ./.ci_with_cmake.sh ${{ matrix.config.cmake_options }} - name: Capture coverage @@ -176,7 +194,7 @@ jobs: with: submodules: true - name: Build Fedora rawhide - run: ./.ci_fedora.sh setup mpich rawhide + run: ./.ci_fedora.sh setup openmpi rawhide shell: bash env: TRAVIS_BUILD_DIR: ${{ github.workspace }} diff --git a/.gitignore b/.gitignore index 7fee4eaae7..7ddf9526ab 100644 --- a/.gitignore +++ b/.gitignore @@ -83,3 +83,5 @@ coverage/ /dist/ /manual/sphinx/_build /_version.txt +/BOUT++-v*.tar.gz +/BOUT++-v*.tar.xz diff --git a/.pip_install_for_ci.sh b/.pip_install_for_ci.sh index 26ca80481a..4a5258cc2d 100755 --- a/.pip_install_for_ci.sh +++ b/.pip_install_for_ci.sh @@ -4,7 +4,7 @@ set -e export PATH=${HOME}/.local/bin:${PATH} pip3 install --user --upgrade pip~=20.0 setuptools~=46.1 -pip3 install --user --upgrade scipy~=1.4 numpy~=1.18 +pip3 install --user --upgrade scipy~=1.4 numpy~=1.18 natsort~=8.1.0 for package in $@ do if test $package == "cython" diff --git a/CMakeLists.txt b/CMakeLists.txt index 8864c9d113..6bb9c963cd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -15,6 +15,13 @@ if(${CMAKE_VERSION} VERSION_GREATER_EQUAL 3.22) cmake_policy(SET CMP0127 NEW) endif() +if ("${CMAKE_CURRENT_BINARY_DIR}" STREQUAL "${CMAKE_CURRENT_SOURCE_DIR}") + option(BOUT_ALLOW_INSOURCE_BUILD "Whether BOUT++ should really allow to build in source." OFF) + if (NOT ${BOUT_ALLOW_INSOURCE_BUILD}) + message(FATAL_ERROR "BOUT++ does not recommend in source builds. Try building out of source, e.g. with `cmake -S . -B build` or set -DBOUT_ALLOW_INSOURCE_BUILD=ON - but things may break!") + endif() +endif() + # CMake currently doesn't support proper semver # Set the version here, strip any extra tags to use in `project` # We try to use git to get a full description, inspired by setuptools_scm @@ -410,8 +417,13 @@ execute_process(COMMAND ${Python3_EXECUTABLE} -c "import site ; print('/'.join(s ) set(CMAKE_INSTALL_PYTHON_SITEARCH lib/${PYTHON_SITEPATH_SUFFIX} CACHE STRING "Location to install python arch-specific modules") -set(BOUT_ENABLE_PYTHON MAYBE CACHE STRING "Build the Python interface") -if (BOUT_ENABLE_PYTHON OR BOUT_ENABLE_PYTHON STREQUAL "MAYBE") +set(ON_OFF_AUTO ON OFF AUTO) +set(BOUT_ENABLE_PYTHON AUTO CACHE STRING "Build the Python interface") +set_property(CACHE BOUT_ENABLE_PYTHON PROPERTY STRINGS ${ON_OFF_AUTO}) +if (NOT BOUT_ENABLE_PYTHON IN_LIST ON_OFF_AUTO) + message(FATAL_ERROR "BOUT_ENABLE_PYTHON must be one of ${ON_OFF_AUTO}") +endif() +if (BOUT_ENABLE_PYTHON OR BOUT_ENABLE_PYTHON STREQUAL "AUTO") add_subdirectory(tools/pylib/_boutpp_build) else() set(BOUT_ENABLE_PYTHON OFF) @@ -723,7 +735,7 @@ endif() add_custom_target(dist - COMMAND ${CMAKE_SOURCE_DIR}/bin/bout-archive-helper.sh v${BOUT_FULL_VERSION} + COMMAND ${Python3_EXECUTABLE} ${CMAKE_SOURCE_DIR}/tools/pylib/_boutpp_build/backend.py dist # there is no cmake equivalent to `mv` - so only works on systems that are not inentionally non-POSIX complient COMMAND mv BOUT++-v${BOUT_FULL_VERSION}.tar.gz ${CMAKE_CURRENT_BINARY_DIR}/ WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} @@ -761,6 +773,26 @@ endif() set(ISINSTALLED "FALSE") +set(_CONFIG_LDFLAGS) +string(REPLACE " " ";" CONFIG_LDFLAGS_LIST ${CONFIG_LDFLAGS}) +foreach (flag ${CONFIG_LDFLAGS_LIST}) + string(REGEX MATCH "^-.*$" isopt "${flag}") + if (isopt) + # message("${flag} is an option") + set(_CONFIG_LDFLAGS "${_CONFIG_LDFLAGS} ${flag}") + # All good + else() + if(EXISTS ${flag}) + # message("${flag} is a file") + set(_CONFIG_LDFLAGS "${_CONFIG_LDFLAGS} ${flag}") + else() + message("Fixing ${flag} to -l${flag}") + set(_CONFIG_LDFLAGS "${_CONFIG_LDFLAGS} -l${flag}") + endif() + endif() +endforeach() +set( CONFIG_LDFLAGS ${_CONFIG_LDFLAGS}) + # This version of the file allows the build directory to be used directly configure_file(bin/bout-config.in bin/bout-config @ONLY) configure_file(tools/pylib/boutconfig/__init__.py.cin tools/pylib/boutconfig/__init__.py @ONLY) diff --git a/cmake/FindPETSc.cmake b/cmake/FindPETSc.cmake index 1ddf95d2d0..c93d464673 100644 --- a/cmake/FindPETSc.cmake +++ b/cmake/FindPETSc.cmake @@ -264,7 +264,7 @@ show : include(Check${PETSC_LANGUAGE_BINDINGS}SourceRuns) - macro (PETSC_TEST_RUNS includes libraries runs) + macro (petsc_test_compiles includes libraries runs) message(STATUS "PETSc test with : ${includes} ${libraries}" ) if (PETSC_VERSION VERSION_GREATER 3.1) set (_PETSC_TSDestroy "TSDestroy(&ts)") @@ -287,12 +287,12 @@ int main(int argc,char *argv[]) { return 0; } ") - multipass_source_runs ("${includes}" "${libraries}" "${_PETSC_TEST_SOURCE}" ${runs} "${PETSC_LANGUAGE_BINDINGS}") + multipass_source_compiles ("${includes}" "${libraries}" "${_PETSC_TEST_SOURCE}" ${runs} "${PETSC_LANGUAGE_BINDINGS}") if (${${runs}}) - set (PETSC_EXECUTABLE_RUNS "YES" CACHE BOOL + set (PETSC_EXECUTABLE_COMPILES "YES" CACHE BOOL "Can the system successfully run a PETSc executable? This variable can be manually set to \"YES\" to force CMake to accept a given PETSc configuration, but this will almost always result in a broken build. If you change PETSC_DIR, PETSC_ARCH, or PETSC_CURRENT you would have to reset this variable." FORCE) endif (${${runs}}) - endmacro (PETSC_TEST_RUNS) + endmacro () find_path (PETSC_INCLUDE_DIR petscts.h @@ -314,14 +314,14 @@ int main(int argc,char *argv[]) { set (petsc_mpi_include_dirs "${MPI_${PETSC_LANGUAGE_BINDINGS}_INCLUDE_DIRS}") #set (petsc_additional_libraries "MPI::MPI_${PETSC_LANGUAGE_BINDINGS}${petsc_openmp_library}") - petsc_test_runs ("${petsc_includes_minimal};${petsc_mpi_include_dirs}" + petsc_test_compiles ("${petsc_includes_minimal};${petsc_mpi_include_dirs}" "${PETSC_LIBRARIES_TS};${petsc_additional_libraries}" petsc_works_minimal) if (petsc_works_minimal) message (STATUS "Minimal PETSc includes and libraries work. This probably means we are building with shared libs.") set (petsc_includes_needed "${petsc_includes_minimal}") else (petsc_works_minimal) # Minimal includes fail, see if just adding full includes fixes it - petsc_test_runs ("${petsc_includes_all};${petsc_mpi_include_dirs}" + petsc_test_compiles ("${petsc_includes_all};${petsc_mpi_include_dirs}" "${PETSC_LIBRARIES_TS};${petsc_additional_libraries}" petsc_works_allincludes) if (petsc_works_allincludes) # It does, we just need all the includes ( @@ -332,7 +332,7 @@ int main(int argc,char *argv[]) { foreach (pkg SYS VEC MAT DM KSP SNES TS ALL) list (APPEND PETSC_LIBRARIES_${pkg} ${petsc_libraries_external}) endforeach (pkg) - petsc_test_runs ("${petsc_includes_minimal};${petsc_mpi_include_dirs}" + petsc_test_compiles ("${petsc_includes_minimal};${petsc_mpi_include_dirs}" "${PETSC_LIBRARIES_TS};${petsc_additional_libraries}" petsc_works_alllibraries) if (petsc_works_alllibraries) @@ -341,7 +341,7 @@ int main(int argc,char *argv[]) { else (petsc_works_alllibraries) # It looks like we really need everything, should have listened to Matt set (petsc_includes_needed ${petsc_includes_all}) - petsc_test_runs ("${petsc_includes_all};${petsc_mpi_include_dirs}" + petsc_test_compiles ("${petsc_includes_all};${petsc_mpi_include_dirs}" "${PETSC_LIBRARIES_TS};${petsc_additional_libraries}" petsc_works_all) if (petsc_works_all) # We fail anyways @@ -372,7 +372,7 @@ if (NOT PETSC_INCLUDES AND NOT TARGET PETSc::PETSc) pkg_search_module(PkgPETSC PETSc>3.4.0 petsc>3.4.0) set (PETSC_LIBRARIES ${PkgPETSC_LINK_LIBRARIES} CACHE STRING "PETSc libraries" FORCE) set (PETSC_INCLUDES ${PkgPETSC_INCLUDE_DIRS} CACHE STRING "PETSc include path" FORCE) - set (PETSC_EXECUTABLE_RUNS "YES" CACHE BOOL + set (PETSC_EXECUTABLE_COMPILES "YES" CACHE BOOL "Can the system successfully run a PETSc executable? This variable can be manually set to \"YES\" to force CMake to accept a given PETSc configuration, but this will almost always result in a broken build. If you change PETSC_DIR, PETSC_ARCH, or PETSC_CURRENT you would have to reset this variable." FORCE) endif() endif() @@ -380,7 +380,7 @@ endif() # Note that we have forced values for all these choices. If you # change these, you are telling the system to trust you that they # work. It is likely that you will end up with a broken build. -mark_as_advanced (PETSC_INCLUDES PETSC_LIBRARIES PETSC_COMPILER PETSC_DEFINITIONS PETSC_MPIEXEC PETSC_EXECUTABLE_RUNS) +mark_as_advanced (PETSC_INCLUDES PETSC_LIBRARIES PETSC_COMPILER PETSC_DEFINITIONS PETSC_MPIEXEC PETSC_EXECUTABLE_COMPILES) include (FindPackageHandleStandardArgs) find_package_handle_standard_args (PETSc diff --git a/cmake/FindPackageMultipass.cmake b/cmake/FindPackageMultipass.cmake index fbf06a7f0f..2452096b56 100644 --- a/cmake/FindPackageMultipass.cmake +++ b/cmake/FindPackageMultipass.cmake @@ -32,6 +32,8 @@ # Always runs the given test, use this when you need to re-run tests # because parent variables have made old cache entries stale. +include(CheckCXXSourceCompiles) + macro (FIND_PACKAGE_MULTIPASS _name _current) string (TOUPPER ${_name} _NAME) set (_args ${ARGV}) @@ -104,3 +106,27 @@ endmacro (MULTIPASS_SOURCE_RUNS) macro (MULTIPASS_C_SOURCE_RUNS includes libraries source runs) multipass_source_runs("${includes}" "${libraries}" "${source}" ${runs} "C") endmacro (MULTIPASS_C_SOURCE_RUNS) + +macro (MULTIPASS_SOURCE_COMPILES includes libraries source runs language) + include (Check${language}SourceRuns) + # This is a ridiculous hack. CHECK_${language}_SOURCE_* thinks that if the + # *name* of the return variable doesn't change, then the test does + # not need to be re-run. We keep an internal count which we + # increment to guarantee that every test name is unique. If we've + # gotten here, then the configuration has changed enough that the + # test *needs* to be rerun. + if (NOT MULTIPASS_TEST_COUNT) + set (MULTIPASS_TEST_COUNT 00) + endif (NOT MULTIPASS_TEST_COUNT) + math (EXPR _tmp "${MULTIPASS_TEST_COUNT} + 1") # Why can't I add to a cache variable? + set (MULTIPASS_TEST_COUNT ${_tmp} CACHE INTERNAL "Unique test ID") + set (testname MULTIPASS_TEST_${MULTIPASS_TEST_COUNT}_${runs}) + set (CMAKE_REQUIRED_INCLUDES ${includes}) + set (CMAKE_REQUIRED_LIBRARIES ${libraries}) + if(${language} STREQUAL "C") + check_c_source_compiles ("${source}" ${testname}) + elseif(${language} STREQUAL "CXX") + check_cxx_source_compiles ("${source}" ${testname}) + endif() + set (${runs} "${${testname}}") +endmacro () diff --git a/cmake/FindSLEPc.cmake b/cmake/FindSLEPc.cmake index 2db1a25a1f..3add8eba8b 100644 --- a/cmake/FindSLEPc.cmake +++ b/cmake/FindSLEPc.cmake @@ -137,10 +137,9 @@ show : endif() if (SLEPC_SKIP_BUILD_TESTS) - set(SLEPC_TEST_RUNS TRUE) set(SLEPC_VERSION "UNKNOWN") set(SLEPC_VERSION_OK TRUE) -elseif (SLEPC_LIBRARIES AND SLEPC_INCLUDE_DIRS AND NOT SLEPC_TEST_RUNS) +elseif (SLEPC_LIBRARIES AND SLEPC_INCLUDE_DIRS) # Set flags for building test program set(CMAKE_REQUIRED_INCLUDES ${SLEPC_INCLUDE_DIRS}) @@ -192,72 +191,6 @@ int main() { set(SLEPC_VERSION_OK TRUE CACHE BOOL "") endif() mark_as_advanced(SLEPC_VERSION_OK) - - # Run SLEPc test program - set(SLEPC_TEST_LIB_CPP - "${CMAKE_CURRENT_BINARY_DIR}/CMakeFiles/slepc_test_lib.cpp") - file(WRITE ${SLEPC_TEST_LIB_CPP} " -#include \"petsc.h\" -#include \"slepceps.h\" -int main() -{ - PetscErrorCode ierr; - int argc = 0; - char** argv = NULL; - ierr = SlepcInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL); - EPS eps; - ierr = EPSCreate(PETSC_COMM_SELF, &eps); CHKERRQ(ierr); - //ierr = EPSSetFromOptions(eps); CHKERRQ(ierr); -#if PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR <= 1 - ierr = EPSDestroy(eps); CHKERRQ(ierr); -#else - ierr = EPSDestroy(&eps); CHKERRQ(ierr); -#endif - ierr = SlepcFinalize(); CHKERRQ(ierr); - return 0; -} -") - - try_run( - SLEPC_TEST_LIB_EXITCODE - SLEPC_TEST_LIB_COMPILED - ${CMAKE_CURRENT_BINARY_DIR} - ${SLEPC_TEST_LIB_CPP} - CMAKE_FLAGS "-DINCLUDE_DIRECTORIES:STRING=${CMAKE_REQUIRED_INCLUDES}" - LINK_LIBRARIES ${CMAKE_REQUIRED_LIBRARIES} - COMPILE_OUTPUT_VARIABLE SLEPC_TEST_LIB_COMPILE_OUTPUT - RUN_OUTPUT_VARIABLE SLEPC_TEST_LIB_OUTPUT - ) - - if (SLEPC_TEST_LIB_COMPILED AND SLEPC_TEST_LIB_EXITCODE EQUAL 0) - message(STATUS "Performing test SLEPC_TEST_RUNS - Success") - set(SLEPC_TEST_RUNS TRUE CACHE BOOL "SLEPc test program can run") - else() - message(STATUS "Performing test SLEPC_TEST_RUNS - Failed") - - # Test program does not run - try adding SLEPc 3rd party libs and test again - list(APPEND CMAKE_REQUIRED_LIBRARIES ${SLEPC_EXTERNAL_LIBRARIES}) - - try_run( - SLEPC_TEST_3RD_PARTY_LIBS_EXITCODE - SLEPC_TEST_3RD_PARTY_LIBS_COMPILED - ${CMAKE_CURRENT_BINARY_DIR} - ${SLEPC_TEST_LIB_CPP} - CMAKE_FLAGS "-DINCLUDE_DIRECTORIES:STRING=${CMAKE_REQUIRED_INCLUDES}" - LINK_LIBRARIES ${CMAKE_REQUIRED_LIBRARIES} - COMPILE_OUTPUT_VARIABLE SLEPC_TEST_3RD_PARTY_LIBS_COMPILE_OUTPUT - RUN_OUTPUT_VARIABLE SLEPC_TEST_3RD_PARTY_LIBS_OUTPUT - ) - - if (SLEPC_TEST_3RD_PARTY_LIBS_COMPILED AND SLEPC_TEST_3RD_PARTY_LIBS_EXITCODE EQUAL 0) - message(STATUS "Performing test SLEPC_TEST_3RD_PARTY_LIBS_RUNS - Success") - set(SLEPC_LIBRARIES ${SLEPC_LIBRARIES} ${SLEPC_EXTERNAL_LIBRARIES} - CACHE STRING "SLEPc libraries." FORCE) - set(SLEPC_TEST_RUNS TRUE CACHE BOOL "SLEPc test program can run") - else() - message(STATUS "Performing test SLEPC_TEST_3RD_PARTY_LIBS_RUNS - Failed") - endif() - endif() endif() # Standard package handling @@ -266,8 +199,7 @@ find_package_handle_standard_args(SLEPc FOUND_VAR SLEPC_FOUND FAIL_MESSAGE "SLEPc could not be found. Be sure to set SLEPC_DIR, PETSC_DIR, and PETSC_ARCH." VERSION_VAR SLEPC_VERSION - REQUIRED_VARS SLEPC_LIBRARIES SLEPC_DIR SLEPC_INCLUDE_DIRS SLEPC_TEST_RUNS - SLEPC_VERSION_OK) + REQUIRED_VARS SLEPC_LIBRARIES SLEPC_DIR SLEPC_INCLUDE_DIRS SLEPC_VERSION_OK) if (SLEPC_FOUND) if (NOT TARGET SLEPc::SLEPc) diff --git a/cmake/FindSUNDIALS.cmake b/cmake/FindSUNDIALS.cmake index 41f97890f7..1ecb5db429 100644 --- a/cmake/FindSUNDIALS.cmake +++ b/cmake/FindSUNDIALS.cmake @@ -31,6 +31,15 @@ include(FindPackageHandleStandardArgs) +find_package(SUNDIALS CONFIG QUIET) + +if (SUNDIALS_FOUND) + if (TARGET SUNDIALS::nvecparallel) + return() + else() + message(STATUS "SUNDIALS found but not SUNDIALS::nvecparallel") + endif() +endif() find_path(SUNDIALS_INCLUDE_DIR sundials_config.h diff --git a/cmake/SetupBOUTThirdParty.cmake b/cmake/SetupBOUTThirdParty.cmake index 00d796fd5d..78dfc6b56e 100644 --- a/cmake/SetupBOUTThirdParty.cmake +++ b/cmake/SetupBOUTThirdParty.cmake @@ -167,16 +167,20 @@ if (BOUT_USE_NETCDF) # Use our own FindnetCDF module which uses nc-config find_package(netCDF REQUIRED) FetchContent_MakeAvailable(netcdf-cxx4) + target_link_libraries(bout++ PUBLIC netCDF::netcdf-cxx4) else() - find_package(netCDFCxx REQUIRED) - foreach(FLAG ${netCDF_CXX_LIBRARY} ${netCDF_LIBRARIES}) - if("${FLAG}" STREQUAL "netcdf") - set(FLAG -lnetcdf) - endif() - set(CONFIG_LDFLAGS "${CONFIG_LDFLAGS} ${FLAG}") - endforeach() + find_package(netCDFCxx) + if (netCDFCxx_FOUND) + set(CONFIG_LDFLAGS "${CONFIG_LDFLAGS} ${netCDF_CXX_LIBRARY} ${netCDF_LIBRARIES}") + target_link_libraries(bout++ PUBLIC netCDF::netcdf-cxx4) + else() + find_package(PkgConfig REQUIRED) + pkg_check_modules(NETCDF REQUIRED IMPORTED_TARGET netcdf-cxx4) + target_link_libraries(bout++ PUBLIC PkgConfig::NETCDF) + list(JOIN NETCDF_LDFLAGS " " NETCDF_LDFLAGS_STRING) + set(CONFIG_LDFLAGS "${CONFIG_LDFLAGS} ${NETCDF_LDFLAGS_STRING}") + endif() endif() - target_link_libraries(bout++ PUBLIC netCDF::netcdf-cxx4) endif() message(STATUS "NetCDF support: ${BOUT_USE_NETCDF}") set(BOUT_HAS_NETCDF ${BOUT_USE_NETCDF}) @@ -190,18 +194,29 @@ endif() message(STATUS "FFTW support: ${BOUT_USE_FFTW}") set(BOUT_HAS_FFTW ${BOUT_USE_FFTW}) -option(BOUT_USE_LAPACK "Enable support for LAPACK" ON) +set(ON_OFF_AUTO ON OFF AUTO) +option(BOUT_USE_LAPACK "Enable support for LAPACK" AUTO) +set_property(CACHE BOUT_USE_LAPACK PROPERTY STRINGS ${ON_OFF_AUTO}) + +set(LAPACK_FOUND OFF) if (BOUT_USE_LAPACK) if (NOT CMAKE_SYSTEM_NAME STREQUAL "CrayLinuxEnvironment") # Cray wrappers sort this out for us - find_package(LAPACK REQUIRED) - target_link_libraries(bout++ PUBLIC "${LAPACK_LIBRARIES}") - string(JOIN " " CONFIG_LAPACK_LIBRARIES ${LAPACK_LIBRARIES}) - set(CONFIG_LDFLAGS "${CONFIG_LDFLAGS} ${CONFIG_LAPACK_LIBRARIES}") + if (BOUT_USE_LAPACK STREQUAL ON) + find_package(LAPACK REQUIRED) + else() + find_package(LAPACK) + endif() + if (LAPACK_FOUND) + target_link_libraries(bout++ PUBLIC "${LAPACK_LIBRARIES}") + string(JOIN " " CONFIG_LAPACK_LIBRARIES ${LAPACK_LIBRARIES}) + set(CONFIG_LDFLAGS "${CONFIG_LDFLAGS} ${CONFIG_LAPACK_LIBRARIES}") + endif() endif() endif() -message(STATUS "LAPACK support: ${BOUT_USE_LAPACK}") -set(BOUT_HAS_LAPACK ${BOUT_USE_LAPACK}) + +message(STATUS "LAPACK support: ${LAPACK_FOUND}") +set(BOUT_HAS_LAPACK ${LAPACK_FOUND}) option(BOUT_USE_SLEPC "Enable support for SLEPc eigen solver" OFF) if (BOUT_USE_SLEPC) diff --git a/examples/laplacexy/simple-hypre/test-laplacexy-hypre.cxx b/examples/laplacexy/simple-hypre/test-laplacexy-hypre.cxx index 154870ec9e..af2367e1a0 100644 --- a/examples/laplacexy/simple-hypre/test-laplacexy-hypre.cxx +++ b/examples/laplacexy/simple-hypre/test-laplacexy-hypre.cxx @@ -13,12 +13,14 @@ int main(int argc, char** argv) { bout::globals::mesh); /// Solution - Field2D x = 0.0; + Field2D solution = 0.0; - x = laplacexy.solve(rhs, x); + solution = laplacexy.solve(rhs, solution); - SAVE_ONCE2(rhs, x); - bout::globals::dump.write(); // Save output file + Options dump; + dump["rhs"] = rhs; + dump["x"] = solution; + bout::writeDefaultOutputFile(dump); } BoutFinalise(); #if BOUT_HAS_CUDA diff --git a/examples/tokamak-2fluid/2fluid.cxx b/examples/tokamak-2fluid/2fluid.cxx index 9108edb5bd..03458d2d4e 100644 --- a/examples/tokamak-2fluid/2fluid.cxx +++ b/examples/tokamak-2fluid/2fluid.cxx @@ -4,6 +4,7 @@ * for cross-benchmarking etc. *******************************************************************************/ +#include "bout/bout_types.hxx" #include #include @@ -98,6 +99,9 @@ class TwoFluid : public PhysicsModel { std::unique_ptr aparSolver; Field2D acoef; // Coefficient in the Helmholtz equation + /// Location of possibly staggered fields + CELL_LOC stagger_loc = CELL_LOC::deflt; + int init(bool UNUSED(restarting)) override { TRACE("int init(bool) "); @@ -154,15 +158,15 @@ class TwoFluid : public PhysicsModel { // READ OPTIONS // Read some parameters - auto globalOptions = Options::root(); - auto options = globalOptions["2fluid"]; + auto& globalOptions = Options::root(); + auto& options = globalOptions["2fluid"]; AA = options["AA"].withDefault(2.0); ZZ = options["ZZ"].withDefault(1.0); estatic = options["estatic"].withDefault(false); ZeroElMass = options["ZeroElMass"].withDefault(false); - zeff = options["zeff"].withDefault(1.0); + zeff = options["Zeff"].withDefault(1.0); nu_perp = options["nu_perp"].withDefault(0.0); ShearFactor = options["ShearFactor"].withDefault(1.0); OhmPe = options["OhmPe"].withDefault(true); @@ -170,8 +174,7 @@ class TwoFluid : public PhysicsModel { curv_upwind = options["curv_upwind"].withDefault(false); // Choose method to use for Poisson bracket advection terms - int bracket_method; - bracket_method = options["bracket_method"].withDefault(0); + int bracket_method = options["bracket_method"].withDefault(0); switch (bracket_method) { case 0: { bm = BRACKET_STD; @@ -206,79 +209,92 @@ class TwoFluid : public PhysicsModel { laplace_extra_rho_term = options["laplace_extra_rho_term"].withDefault(false); vort_include_pi = options["vort_include_pi"].withDefault(false); - lowPass_z = options["lowPass_z"].withDefault(-1); + lowPass_z = options["low_pass_z"].withDefault(-1); + + //////////////////////////////////////////////////////// + // Equation terms + + auto& ni_options = globalOptions["Ni"]; + ni_options.setConditionallyUsed(); + evolve_ni = ni_options["evolve"].withDefault(true); + + auto& rho_options = globalOptions["rho"]; + rho_options.setConditionallyUsed(); + evolve_rho = rho_options["evolve"].withDefault(true); + + auto& vi_options = globalOptions["Vi"]; + vi_options.setConditionallyUsed(); + evolve_vi = vi_options["evolve"].withDefault(true); + + auto& te_options = globalOptions["Te"]; + te_options.setConditionallyUsed(); + evolve_ti = te_options["evolve"].withDefault(true); - evolve_ni = globalOptions["Ni"]["evolve"].withDefault(true); - evolve_rho = globalOptions["rho"]["evolve"].withDefault(true); - evolve_vi = globalOptions["vi"]["evolve"].withDefault(true); - evolve_ti = globalOptions["ti"]["evolve"].withDefault(true); - evolve_ajpar = globalOptions["Ajpar"]["evolve"].withDefault(true); + auto& ti_options = globalOptions["Ti"]; + ti_options.setConditionallyUsed(); + evolve_ti = ti_options["evolve"].withDefault(true); + + auto& ajpar_options = globalOptions["Ajpar"]; + ajpar_options.setConditionallyUsed(); + evolve_ajpar = ajpar_options["evolve"].withDefault(true); if (ZeroElMass) { - evolve_ajpar = 0; // Don't need ajpar - calculated from ohm's law + evolve_ajpar = false; // Don't need ajpar - calculated from ohm's law } - //////////////////////////////////////////////////////// - // Equation terms - if (evolve_ni) { - Options* options = &globalOptions["Ni"]; - options->get("ni1_phi0", ni_ni1_phi0, false); - options->get("ni0_phi1", ni_ni0_phi1, false); - options->get("ni1_phi1", ni_ni1_phi1, false); - options->get("nit_phit", ni_nit_phit, false); - options->get("vi1_ni0", ni_vi1_ni0, false); - options->get("vi0_ni1", ni_vi0_ni1, false); - options->get("vi1_ni1", ni_vi1_ni1, false); - options->get("vit_nit", ni_vit_nit, false); - options->get("jpar1", ni_jpar1, false); - options->get("pe1", ni_pe1, false); - options->get("ni0_curv_phi1", ni_ni0_curv_phi1, false); - options->get("ni1_curv_phi0", ni_ni1_curv_phi0, false); - options->get("ni1_curv_phi1", ni_ni1_curv_phi1, false); - options->get("nit_curv_phit", ni_nit_curv_phit, false); + ni_ni1_phi0 = ni_options["ni1_phi0"].withDefault(false); + ni_ni0_phi1 = ni_options["ni0_phi1"].withDefault(false); + ni_ni1_phi1 = ni_options["ni1_phi1"].withDefault(false); + ni_nit_phit = ni_options["nit_phit"].withDefault(false); + ni_vi1_ni0 = ni_options["vi1_ni0"].withDefault(false); + ni_vi0_ni1 = ni_options["vi0_ni1"].withDefault(false); + ni_vi1_ni1 = ni_options["vi1_ni1"].withDefault(false); + ni_vit_nit = ni_options["vit_nit"].withDefault(false); + ni_jpar1 = ni_options["jpar1"].withDefault(false); + ni_pe1 = ni_options["pe1"].withDefault(false); + ni_ni0_curv_phi1 = ni_options["ni0_curv_phi1"].withDefault(false); + ni_ni1_curv_phi0 = ni_options["ni1_curv_phi0"].withDefault(false); + ni_ni1_curv_phi1 = ni_options["ni1_curv_phi1"].withDefault(false); + ni_nit_curv_phit = ni_options["nit_curv_phit"].withDefault(false); } if (evolve_rho) { - Options* options = &globalOptions["rho"]; - options->get("rho0_phi1", rho_rho0_phi1, false); - options->get("rho1_phi0", rho_rho1_phi0, false); - options->get("rho1_phi1", rho_rho1_phi1, false); - options->get("vi1_rho0", rho_vi1_rho0, false); - options->get("vi0_rho1", rho_vi0_rho1, false); - options->get("vi1_rho1", rho_vi1_rho1, false); - options->get("pei1", rho_pei1, false); - options->get("jpar1", rho_jpar1, false); - options->get("rho1", rho_rho1, false); + rho_rho0_phi1 = rho_options["rho0_phi1"].withDefault(false); + rho_rho1_phi0 = rho_options["rho1_phi0"].withDefault(false); + rho_rho1_phi1 = rho_options["rho1_phi1"].withDefault(false); + rho_vi1_rho0 = rho_options["vi1_rho0"].withDefault(false); + rho_vi0_rho1 = rho_options["vi0_rho1"].withDefault(false); + rho_vi1_rho1 = rho_options["vi1_rho1"].withDefault(false); + rho_pei1 = rho_options["pei1"].withDefault(false); + rho_jpar1 = rho_options["jpar1"].withDefault(false); + rho_rho1 = rho_options["rho1"].withDefault(false); } if (evolve_vi) { - Options* options = &globalOptions["vi"]; - options->get("vi0_phi1", vi_vi0_phi1, false); - options->get("vi1_phi0", vi_vi1_phi0, false); - options->get("vi1_phi1", vi_vi1_phi1, false); - options->get("vit_phit", vi_vit_phit, false); - options->get("vi1_vi0", vi_vi1_vi0, false); - options->get("vi0_vi1", vi_vi0_vi1, false); - options->get("vi1_vi1", vi_vi1_vi1, false); - options->get("vit_vit", vi_vit_vit, false); - options->get("pei1", vi_pei1, false); - options->get("peit", vi_peit, false); - options->get("vi1", vi_vi1, false); + vi_vi0_phi1 = vi_options["vi0_phi1"].withDefault(false); + vi_vi1_phi0 = vi_options["vi1_phi0"].withDefault(false); + vi_vi1_phi1 = vi_options["vi1_phi1"].withDefault(false); + vi_vit_phit = vi_options["vit_phit"].withDefault(false); + vi_vi1_vi0 = vi_options["vi1_vi0"].withDefault(false); + vi_vi0_vi1 = vi_options["vi0_vi1"].withDefault(false); + vi_vi1_vi1 = vi_options["vi1_vi1"].withDefault(false); + vi_vit_vit = vi_options["vit_vit"].withDefault(false); + vi_pei1 = vi_options["pei1"].withDefault(false); + vi_peit = vi_options["peit"].withDefault(false); + vi_vi1 = vi_options["vi1"].withDefault(false); } if (evolve_te) { - Options* options = &globalOptions["te"]; - options->get("te1_phi0", te_te1_phi0, false); - options->get("te0_phi1", te_te0_phi1, false); - options->get("te1_phi1", te_te1_phi1, false); + te_te1_phi0 = te_options["te1_phi0"].withDefault(false); + te_te0_phi1 = te_options["te0_phi1"].withDefault(false); + te_te1_phi1 = te_options["te1_phi1"].withDefault(false); } if (evolve_ti) { - Options* options = &globalOptions["ti"]; - options->get("ti1_phi0", ti_ti1_phi0, false); - options->get("ti0_phi1", ti_ti0_phi1, false); - options->get("ti1_phi1", ti_ti1_phi1, false); + ti_ti1_phi0 = ti_options["ti1_phi0"].withDefault(false); + ti_ti0_phi1 = ti_options["ti0_phi1"].withDefault(false); + ti_ti1_phi1 = ti_options["ti1_phi1"].withDefault(false); } //////////////////////////////////////////////////////// @@ -337,13 +353,15 @@ class TwoFluid : public PhysicsModel { //////////////////////////////////////////////////////// // SHIFTED GRIDS LOCATION + stagger_loc = CELL_LOC::ylow; + // Velocities defined on cell boundaries - Vi.setLocation(CELL_YLOW); - Ajpar.setLocation(CELL_YLOW); + Vi.setLocation(stagger_loc); + Ajpar.setLocation(stagger_loc); // Apar and jpar too - Apar.setLocation(CELL_YLOW); - jpar.setLocation(CELL_YLOW); + Apar.setLocation(stagger_loc); + jpar.setLocation(stagger_loc); } //////////////////////////////////////////////////////// @@ -490,16 +508,17 @@ class TwoFluid : public PhysicsModel { dump.addOnce(wci, "wci"); // Create a solver for the Laplacian - phiSolver = Laplacian::create(&options["phiSolver"]); + phiSolver = Laplacian::create(&globalOptions["phiSolver"]); if (laplace_extra_rho_term) { // Include the first order term Grad_perp Ni dot Grad_perp phi phiSolver->setCoefC(Ni0); } + globalOptions["aparSolver"].setConditionallyUsed(); + if (!(estatic || ZeroElMass)) { // Create a solver for the electromagnetic potential - aparSolver = Laplacian::create(&options["aparSolver"], - mesh->StaggerGrids ? CELL_YLOW : CELL_CENTRE); + aparSolver = Laplacian::create(&globalOptions["aparSolver"], stagger_loc); if (mesh->StaggerGrids) { acoef = (-0.5 * beta_p / fmei) * interp_to(Ni0, CELL_YLOW); } else { @@ -599,10 +618,13 @@ class TwoFluid : public PhysicsModel { if (ZeroElMass) { // Set jpar,Ve,Ajpar neglecting the electron inertia term // Calculate Jpar, communicating across processors - jpar = -(Ni0 * Grad_par(phi, CELL_YLOW)) / (fmei * 0.51 * nu); + + jpar = + -interp_to(Ni0, stagger_loc) * Grad_par(phi, stagger_loc) / (fmei * 0.51 * nu); if (OhmPe) { - jpar += (Te0 * Grad_par(Ni, CELL_YLOW)) / (fmei * 0.51 * nu); + jpar += + interp_to(Te0, stagger_loc) * Grad_par(Ni, stagger_loc) / (fmei * 0.51 * nu); } // Need to communicate jpar diff --git a/examples/tokamak-2fluid/data/BOUT.inp b/examples/tokamak-2fluid/data/BOUT.inp index 563af66cd6..214bafc19a 100644 --- a/examples/tokamak-2fluid/data/BOUT.inp +++ b/examples/tokamak-2fluid/data/BOUT.inp @@ -45,14 +45,6 @@ first = C4 second = C4 upwind = U1 -################################################## -# Laplacian inversion settings - -[laplace] -all_terms = false -nonuniform = true -filter = 0.2 # Remove the top 20% of modes (BOUT-06 zwindow=0.4) - ################################################## # Solver settings @@ -81,8 +73,6 @@ estatic = true # if true, electrostatic (Apar = 0). (BOUT-06 = esop) ZeroElMass = true # Use Ohms law without electron inertia bout_jpar = true # Use BOUT-06 method to calculate ZeroElMass jpar -bout_exb = true # Use the BOUT-06 subset of ExB terms - curv_upwind = false # Use upwinding for b0xkappa_dot_Grad terms laplace_extra_rho_term = false # include Grad_perp(Ni) dot Grad_perp(phi) term diff --git a/externalpackages/PVODE/include/pvode/nvector.h b/externalpackages/PVODE/include/pvode/nvector.h index 8232995847..1095bdb474 100644 --- a/externalpackages/PVODE/include/pvode/nvector.h +++ b/externalpackages/PVODE/include/pvode/nvector.h @@ -88,12 +88,13 @@ namespace pvode { /* Set types real and integer for MPI calls. */ - +#ifndef PVEC_REAL_MPI_TYPE #if (LLNL_DOUBLE == 1) #define PVEC_REAL_MPI_TYPE MPI_DOUBLE #else #define PVEC_REAL_MPI_TYPE MPI_FLOAT #endif +#endif // Ugh, potentionanlly bad but need to do this for now when compiling with both SUNDIALS and PVODE #ifndef PVEC_INTEGER_MPI_TYPE diff --git a/externalpackages/PVODE/precon/band.cpp b/externalpackages/PVODE/precon/band.cpp index c1d0f6f1e2..43151320d0 100644 --- a/externalpackages/PVODE/precon/band.cpp +++ b/externalpackages/PVODE/precon/band.cpp @@ -185,8 +185,8 @@ integer gbfa(real **a, integer n, integer mu, integer ml, integer smu, if (col_k[storage_l] == ZERO) return(k+1); /* swap a(l,k) and a(k,k) if necessary */ - - if (swap = (l != k)) { + swap = (l != k); + if (swap) { temp = col_k[storage_l]; col_k[storage_l] = *diag_k; *diag_k = temp; diff --git a/externalpackages/PVODE/source/cvode.cpp b/externalpackages/PVODE/source/cvode.cpp index e28a3b1ae6..135662a40f 100644 --- a/externalpackages/PVODE/source/cvode.cpp +++ b/externalpackages/PVODE/source/cvode.cpp @@ -179,7 +179,7 @@ namespace pvode { #define MSG_Y0_NULL CVM "y0=NULL illegal.\n\n" -#define MSG_BAD_N CVM "N=%ld < 1 illegal.\n\n" +#define MSG_BAD_N CVM "N=%d < 1 illegal.\n\n" #define MSG_BAD_LMM_1 CVM "lmm=%d illegal.\n" #define MSG_BAD_LMM_2 "The legal values are ADAMS=%d and BDF=%d.\n\n" @@ -1826,6 +1826,8 @@ static int CVnls(CVodeMem cv_mem, int nflag) case FUNCTIONAL : return(CVnlsFunctional(cv_mem)); case NEWTON : return(CVnlsNewton(cv_mem, nflag)); } + fprintf(errfp, "Should be unreachable ..."); + return (ERR_FAILURE); } /***************** CVnlsFunctional ******************************** @@ -2392,6 +2394,8 @@ static int CVHandleFailure(CVodeMem cv_mem, int kflag) case SOLVE_FAILED: fprintf(errfp, MSG_SOLVE_FAILED, tn); return(SOLVE_FAILURE); } + fprintf(errfp, "Should be unreachable ..."); + return (ERR_FAILURE); } /*******************************************************************/ diff --git a/externalpackages/boutdata b/externalpackages/boutdata index 2d4c5bb4b2..ab59ef9138 160000 --- a/externalpackages/boutdata +++ b/externalpackages/boutdata @@ -1 +1 @@ -Subproject commit 2d4c5bb4b2bdc0b74d3d267559a5624aa599b95f +Subproject commit ab59ef913884918a5bc5ee84101e6d6b833dcd6c diff --git a/externalpackages/boututils b/externalpackages/boututils index a79a00a54f..c2e97a226a 160000 --- a/externalpackages/boututils +++ b/externalpackages/boututils @@ -1 +1 @@ -Subproject commit a79a00a54f69663117a93dd42ffa6004783432c8 +Subproject commit c2e97a226a8a7728b8f5085792f7744484943512 diff --git a/include/bout/boutcomm.hxx b/include/bout/boutcomm.hxx index 77038844a8..fea401af02 100644 --- a/include/bout/boutcomm.hxx +++ b/include/bout/boutcomm.hxx @@ -40,7 +40,7 @@ public: static BoutComm* getInstance(); /// Shortcut method - static MPI_Comm get(); + static MPI_Comm& get(); static void setArgs(int& c, char**& v); @@ -53,7 +53,7 @@ public: void setComm(MPI_Comm c); // Getters - MPI_Comm getComm(); + MPI_Comm& getComm(); bool isSet(); private: diff --git a/include/bout/field.hxx b/include/bout/field.hxx index 51719de82e..be347f60d0 100644 --- a/include/bout/field.hxx +++ b/include/bout/field.hxx @@ -345,7 +345,8 @@ inline bool isUniform(const T& f, bool allpe = false, auto element = f[*f.getRegion(region).begin()]; // TODO: maybe parallise this loop, as the early return is unlikely BOUT_FOR_SERIAL(i, f.getRegion(region)) { - if (f[i] != element) { + // by default we only check for exact equality, as that should cover most cases + if (f[i] != element and (not almost_equal(f[i], element, 10))) { result = false; break; } @@ -370,8 +371,12 @@ inline BoutReal getUniform(const T& f, MAYBE_UNUSED(bool allpe) = false, const std::string& region = "RGN_ALL") { #if CHECK > 1 if (not isUniform(f, allpe, region)) { - throw BoutException("Requested getUniform({}, {}, {}) but Field is not const", f.name, - allpe, region); + const BoutReal f1 = min(f, allpe, region); + const BoutReal f2 = max(f, allpe, region); + throw BoutException("Requested getUniform({}, {}, {}) but Field is not const " + "([{:.15f}...{:.15f}] Δ={:e} {:e}ε)", + f.name, allpe, region, f1, f2, f2 - f1, + (f2 - f1) / (f1 + f2) / std::numeric_limits::epsilon()); } #endif return f[*f.getRegion(region).begin()]; diff --git a/include/bout/invertable_operator.hxx b/include/bout/invertable_operator.hxx index 49b254187a..7168324b75 100644 --- a/include/bout/invertable_operator.hxx +++ b/include/bout/invertable_operator.hxx @@ -430,7 +430,8 @@ public: KSPConvergedReason reason; ierr = KSPGetConvergedReason(ksp, &reason); if (reason <= 0) { - throw BoutException("KSPSolve failed with reason {:d}.", reason); + throw BoutException("KSPSolve failed. Reason {} ({:d})", + KSPConvergedReasons[reason], static_cast(reason)); } // Probably want to remove the following in the long run diff --git a/include/bout/output.hxx b/include/bout/output.hxx index bb202cc8df..ef09aa7ee5 100644 --- a/include/bout/output.hxx +++ b/include/bout/output.hxx @@ -137,11 +137,11 @@ private: class DummyOutput : public Output { public: template - void write(const S&, const Args&...){}; + void write([[maybe_unused]] const S& format, [[maybe_unused]] const Args&... args){}; template - void print(const S&, const Args&...){}; - void write(const std::string& s) override{}; - void print(const std::string& s) override{}; + void print([[maybe_unused]] const S& format, [[maybe_unused]] const Args&... args){}; + void write([[maybe_unused]] const std::string& message) override{}; + void print([[maybe_unused]] const std::string& message) override{}; void enable() override{}; void disable() override{}; void enable(MAYBE_UNUSED(bool enable)){}; diff --git a/include/bout/petsclib.hxx b/include/bout/petsclib.hxx index e731405e1c..dd3daa28aa 100644 --- a/include/bout/petsclib.hxx +++ b/include/bout/petsclib.hxx @@ -64,6 +64,10 @@ class Options; #include #include +#include "bout/boutexception.hxx" + +#define BOUT_DO_PETSC(cmd) PetscLib::assertIerr(cmd, #cmd) + /*! * Handles initialisation and finalisation of PETSc library. * The first instance which is created initialises PETSc @@ -111,6 +115,14 @@ public: */ static void cleanup(); + static inline void assertIerr(PetscErrorCode ierr, std::string op = "PETSc operation") { + if (ierr) { + throw BoutException("{:s} failed with {:d}", op, ierr); + } + } + + static BoutException SNESFailure(SNES& snes); + private: static int count; ///< How many instances? static char help[]; ///< Help string diff --git a/include/bout/sundials_backports.hxx b/include/bout/sundials_backports.hxx index 9c96de1935..c5e0f3ab15 100644 --- a/include/bout/sundials_backports.hxx +++ b/include/bout/sundials_backports.hxx @@ -39,7 +39,7 @@ inline void SUNNonlinSolFree(MAYBE_UNUSED(SUNNonlinearSolver solver)) {} #if SUNDIALS_VERSION_MAJOR < 6 namespace sundials { struct Context { - Context(MPI_Comm comm MAYBE_UNUSED() = MPI_COMM_NULL) {} + Context(void* comm MAYBE_UNUSED()) {} }; } // namespace sundials diff --git a/include/bout/utils.hxx b/include/bout/utils.hxx index d221abcbad..450d3eaf87 100644 --- a/include/bout/utils.hxx +++ b/include/bout/utils.hxx @@ -661,30 +661,16 @@ std::string trimComments(const std::string& s, const std::string& c = "#;"); /// https://en.wikipedia.org/wiki/Damerau%E2%80%93Levenshtein_distance#Optimal_string_alignment_distance std::string::size_type editDistance(const std::string& str1, const std::string& str2); -/// the bout_vsnprintf macro: -/// The first argument is an char * buffer of length len. -/// It needs to have been allocated with new[], as it may be -/// reallocated. -/// len: the length of said buffer. May be changed, mussn't be const. -/// fmt: the const char * descriping the format. -/// note that fmt should be the first argument of the function of type -/// const char * and has to be directly followed by the variable arguments. -#define bout_vsnprintf(buf, len, fmt) \ - { \ - va_list va; \ - va_start(va, fmt); \ - int _vsnprintflen = vsnprintf(buf, len, fmt, va); \ - va_end(va); \ - if (_vsnprintflen + 1 > int(len)) { \ - _vsnprintflen += 1; \ - delete[] buf; \ - buf = new char[_vsnprintflen]; \ - len = _vsnprintflen; \ - va_start(va, fmt); \ - vsnprintf(buf, len, fmt, va); \ - va_end(va); \ - } \ - } +// from https://en.cppreference.com/w/cpp/types/numeric_limits/epsilon +template +typename std::enable_if::is_integer, bool>::type +almost_equal(T x, T y, int ulp = 2) { + // the machine epsilon has to be scaled to the magnitude of the values used + // and multiplied by the desired precision in ULPs (units in the last place) + return std::fabs(x - y) <= std::numeric_limits::epsilon() * std::fabs(x + y) * ulp + // unless the result is subnormal + || std::fabs(x - y) < std::numeric_limits::min(); +} /// Convert pointer or reference to pointer /// This allows consistent handling of both in macros, templates diff --git a/manual/Makefile b/manual/Makefile index 1e1d674a3f..6da41e4a71 100644 --- a/manual/Makefile +++ b/manual/Makefile @@ -1,5 +1,5 @@ # Makefile for the reference and user manuals -.PHONY:all clean sphinx html pdf doxygen breathe-autogen +.PHONY:all clean sphinx html pdf doxygen breathe-autogen maybe-doxygen BOUT_TOP?=.. @@ -14,26 +14,30 @@ manual: all pdf: sphinx-pdf html: sphinx-html man: sphinx-man -sphinx: sphinx-pdf +sphinx: sphinx-html -sphinx-pdf: doxygen +sphinx-pdf: maybe-doxygen PYTHONPATH=$(BOUT_TOP)/tools/pylib:$$PYTHONPATH $(sphinx-build) -b latex sphinx/ build/ cd build && latexmk -pdf BOUT -interaction=batchmode test -e BOUT.pdf || ln -s build/BOUT.pdf . @echo "Documentation is available in $$(pwd)/BOUT.pdf" -sphinx-html: doxygen +sphinx-html: maybe-doxygen PYTHONPATH=$(BOUT_TOP)/tools/pylib:$$PYTHONPATH $(sphinx-build) -b html sphinx/ html/ @echo "Documentation available in file://$$(pwd)/html/index.html" -sphinx-man: doxygen +sphinx-man: maybe-doxygen PYTHONPATH=$(BOUT_TOP)/tools/pylib:$$PYTHONPATH $(sphinx-build) -b man sphinx/ man/ @echo "Documentation available in $$(pwd)/man/bout.1" # Run doxygen, ignore if it fails (leading '-') -doxygen: +maybe-doxygen: -cd doxygen && doxygen Doxyfile + +doxygen: + cd doxygen && doxygen Doxyfile + # Run breathe-apidoc, ignore if it fails (leading '-') breathe-autogen: doxygen -breathe-apidoc -f -o sphinx/_breathe_autogen doxygen/bout/xml diff --git a/manual/README.md b/manual/README.md index 1003f0ed19..9f69c00548 100644 --- a/manual/README.md +++ b/manual/README.md @@ -12,30 +12,28 @@ To build the manual locally, you need at least "sphinx" and "recommonmark", which you can install using pip (or pip3): ```bash -$ pip install --user sphinx -$ pip install --user recommonmark +$ pip install -r sphinx/requirements.txt ``` -These documents can be built into a PDF using "sphinx-build": +To get a local html version, run ```bash -$ make +$ make html ``` +This should create a file "index.html" in the "manual/html" directory. + To use e.g. "sphinx-build-3" instead of "sphinx-build", run ```bash $ make sphinx-build=sphinx-build-3 ``` -This should create a file "BOUT.pdf" in the "manual" directory. - -To get a local html version, run +These documents can be built into a PDF using "sphinx-build": ```bash -$ make html +$ make ``` - -This should create a file "index.html" in the "manual/html" directory. +This should create a file "BOUT.pdf" in the "manual" directory. ### API documentation @@ -53,14 +51,7 @@ It is possible to build the API documentation into the main manual using "breathe". Install breathe: ```bash -$ pip install --user breathe -``` - -and comment out the following line in `sphinx/conf.py`: - -```python -# Disable breathe -has_breathe = False +$ pip install breathe ``` You can then build the sphinx documentation as normal. diff --git a/manual/RELEASE_HOWTO.md b/manual/RELEASE_HOWTO.md index 98eb013c69..4103e6f197 100644 --- a/manual/RELEASE_HOWTO.md +++ b/manual/RELEASE_HOWTO.md @@ -20,6 +20,7 @@ releases - Raise issues for any tests that fail - Possibly run `clang-tidy`, `clang-check`, `coverity`, etc. - [ ] Review pinned pip package versions for CI +- [ ] Review bundled libraries Before merging PR: @@ -54,6 +55,7 @@ Before merging PR: After PR is merged: +- [ ] Make tarball: `make dist` from build directory. Ensure you are on a tag and correct version is used for archive and folder within. - [ ] Try to summarise the changes! - [ ] Make [GitHub Release][gh_release], include change summary **NB:** tag should have leading `v` diff --git a/manual/sphinx/_apidoc/boutpp.rst b/manual/sphinx/_apidoc/boutpp.rst index d28d7c7d18..71a092ca17 100644 --- a/manual/sphinx/_apidoc/boutpp.rst +++ b/manual/sphinx/_apidoc/boutpp.rst @@ -13,16 +13,4 @@ Module contents .. automodule:: boutpp :members: - -Undocumented functions ----------------------- - -.. automodule:: boutpp - :undoc-members: - -Special and inherited functions -------------------------------- - -.. automodule:: boutpp - :special-members: - :inherited-members: + :imported-members: diff --git a/manual/sphinx/conf.py b/manual/sphinx/conf.py index a643774553..1e73d9febc 100755 --- a/manual/sphinx/conf.py +++ b/manual/sphinx/conf.py @@ -49,7 +49,6 @@ def __getattr__(cls, name): return MagicMock() MOCK_MODULES = [ - "h5py", "netCDF4", "mayavi2", "enthought", @@ -67,6 +66,8 @@ def __getattr__(cls, name): "scipy.ndimage.filters", "scipy.ndimage.morphology", "scipy.spatial", + "past", + "crosslines", ] sys.modules.update((mod_name, Mock()) for mod_name in MOCK_MODULES) print(os.environ) @@ -86,6 +87,7 @@ def __getattr__(cls, name): + " -DBOUT_ENABLE_PYTHON=ON" + " -DBOUT_UPDATE_GIT_SUBMODULE=OFF" + " -DBOUT_TESTS=OFF" + + " -DBOUT_ALLOW_INSOURCE_BUILD=ON" + f" -DPython_ROOT_DIR={pydir}" + f" -Dmpark_variant_DIR={pwd}/externalpackages/mpark.variant/" + f" -Dfmt_DIR={pwd}/externalpackages/fmt/" @@ -174,8 +176,8 @@ def __getattr__(cls, name): # General information about the project. project = "BOUT++" -copyright = "2017, B. Dudson" -author = "The BOUT++ team" +copyright = "2017-2023" +author = "B. Dudson and The BOUT++ team" # The version info for the project you're documenting, acts as replacement for # |version| and |release|, also used in various other places throughout the @@ -231,6 +233,7 @@ def __getattr__(cls, name): # html_theme_options = dict( repository_url="https://github.com/boutproject/BOUT-dev", + repository_branch="master", path_to_docs="manual/sphinx", use_edit_page_button=True, use_repository_button=True, diff --git a/manual/sphinx/developer_docs/file_io.rst b/manual/sphinx/developer_docs/file_io.rst index ad44c97277..aece8cfb1e 100644 --- a/manual/sphinx/developer_docs/file_io.rst +++ b/manual/sphinx/developer_docs/file_io.rst @@ -92,7 +92,7 @@ similar interface to ``Datafile``, and just like ``Datafile`` it also stores pointers to variables. This means that it still suffers from all of the downsides of ``Datafile``, and developers are encouraged to move to the ``outputVars`` approach. The `SAVE_ONCE`/`SAVE_REPEAT` -macros also work through `DataFileFacade`` -- this means that they +macros also work through ``DataFileFacade`` -- this means that they cannot be used outside of ``PhysicsModel`` methods! diff --git a/manual/sphinx/user_docs/advanced_install.rst b/manual/sphinx/user_docs/advanced_install.rst index fddbbcec8b..b2b6b49c80 100644 --- a/manual/sphinx/user_docs/advanced_install.rst +++ b/manual/sphinx/user_docs/advanced_install.rst @@ -113,32 +113,38 @@ Cori First set up the environment by loading the correct modules. For Bash shell use: .. code-block:: bash + source config/cori/setup-env-cgpu.sh and for C shell: .. code-block:: csh + source config/cori/setup-env-cgpu.sh Then configure BOUT++ by running a script which calls CMake. Under bash: .. code-block:: bash + ./config/cori/config-bout-cgpu.sh and C shell: .. code-block:: csh + ./config/cori/config-bout-cgpu.csh At the time of writing, Hypre linking is not working with CUDA. If you come across errors with the above configuration, try turning off Hypre support: .. code-block:: bash + ./config/cori/config-bout-cgpu-nohypre.sh or .. code-block:: csh + ./config/cori/config-bout-cgpu-nohypre.csh See section :ref:`sec-gpusupport` for details of compiling and running diff --git a/manual/sphinx/user_docs/installing.rst b/manual/sphinx/user_docs/installing.rst index b09d219543..7f8c26ffbe 100644 --- a/manual/sphinx/user_docs/installing.rst +++ b/manual/sphinx/user_docs/installing.rst @@ -443,27 +443,27 @@ configuration:: If not, see :ref:`sec-advancedinstall` for some things you can try to resolve common problems. -Working with an active `conda` environment +Working with an active ``conda`` environment ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -When `conda` is used, it installs separate versions of several libraries. These +When ``conda`` is used, it installs separate versions of several libraries. These can cause warnings or even failures when linking BOUT++ executables. There are several alternatives to deal with this problem: -* The simplest but least convenient option is to use `conda deactivate` before +* The simplest but least convenient option is to use ``conda deactivate`` before configuring, compiling, or running any BOUT++ program. * You might sometimes want to link to the conda-installed libraries. This is probably not ideal for production runs on an HPC system (as conda downloads binary packages that will not be optimized for specific hardware), but can be a simple way to get packages for testing or on a personal computer. In this - case just keep your `conda` environment active, and with luck the libraries + case just keep your ``conda`` environment active, and with luck the libraries should be picked up by the standard search mechanisms. * In case you do want a fully optimized and as-stable-as-possible build for production runs, it is probably best not to depend on any conda packages for - compiling or running BOUT++ executables (restrict `conda` to providing Python + compiling or running BOUT++ executables (restrict ``conda`` to providing Python packages for post-processing, and their dependencies). Passing - `-DBOUT_IGNORE_CONDA_ENV=ON` (default `OFF`) excludes anything in the conda + ``-DBOUT_IGNORE_CONDA_ENV=ON`` (default ``OFF``) excludes anything in the conda environment from CMake search paths. This should totally separate BOUT++ from - the `conda` environment. + the ``conda`` environment. .. _sec-config-nls: @@ -520,11 +520,11 @@ You can also install all the packages directly (see the documentation in the `bo `__ repos for the most up to date list) using pip:: - $ pip install --user numpy scipy matplotlib sympy netCDF4 h5py future importlib-metadata + $ pip install --user numpy scipy matplotlib sympy netCDF4 future importlib-metadata or conda:: - $ conda install numpy scipy matplotlib sympy netcdf4 h5py future importlib-metadata + $ conda install numpy scipy matplotlib sympy netcdf4 future importlib-metadata They may also be available from your Linux system's package manager. diff --git a/manual/sphinx/user_docs/output_and_post.rst b/manual/sphinx/user_docs/output_and_post.rst index ff32da04d0..71f4098113 100644 --- a/manual/sphinx/user_docs/output_and_post.rst +++ b/manual/sphinx/user_docs/output_and_post.rst @@ -30,7 +30,7 @@ Requirements The Python tools provided with BOUT++ make heavy use of numpy_ and scipy_, as well as matplotlib_ for the plotting routines. In order -to read BOUT++ output in Python, you will need either netcdf4_ or h5py_. +to read BOUT++ output in Python, you will need either netcdf4_. While we try to ensure that the Python tools are compatible with both Python 2 and 3, we officially only support Python 3. @@ -55,7 +55,6 @@ supported versions of numpy, scipy, netcdf4, matplotlib and jinja2. .. _scipy: http://www.scipy.org/ .. _matplotlib: https://www.matplotlib.org .. _netcdf4: http://unidata.github.io/netcdf4-python/ -.. _h5py: http://www.h5py.org .. _Jinja2: http://jinja.pocoo.org/ .. _installation instructions: https://www.scipy.org/install.html diff --git a/manual/sphinx/user_docs/parallel-transforms.rst b/manual/sphinx/user_docs/parallel-transforms.rst index ad7295bbee..3ee3eccfb8 100644 --- a/manual/sphinx/user_docs/parallel-transforms.rst +++ b/manual/sphinx/user_docs/parallel-transforms.rst @@ -70,8 +70,8 @@ setting is .. code-block:: cfg - [mesh] - paralleltransform = identity + [mesh:paralleltransform] + type = identity This then uses the `ParallelTransformIdentity` class to calculate the @@ -100,8 +100,8 @@ The shifted metric method is selected using: .. code-block:: cfg - [mesh] - paralleltransform = shifted + [mesh:paralleltransform] + type = shifted so that mesh uses the `ShiftedMetric` class to calculate parallel transforms. During initialisation, this class reads a quantity zShift @@ -141,8 +141,8 @@ calculation of parallel slices. Select it by using: .. code-block:: cfg - [mesh] - paralleltransform = shifted + [mesh:paralleltransform] + type = shifted calcParallelSlices_on_communicate = false With these settings, inputs to parallel derivative or interpolation @@ -177,8 +177,8 @@ To use the FCI method for parallel transforms, set .. code-block:: cfg - [mesh] - paralleltransform = fci + [mesh:paralleltransform] + type = fci which causes the `FCITransform` class to be used for parallel transforms. This reads four variables (3D fields) from the input diff --git a/manual/sphinx/user_docs/time_integration.rst b/manual/sphinx/user_docs/time_integration.rst index 9d5a34ea89..7d0b56abcc 100644 --- a/manual/sphinx/user_docs/time_integration.rst +++ b/manual/sphinx/user_docs/time_integration.rst @@ -151,6 +151,20 @@ for a certain variable (e.g. ``[n]``), setting the option ``positivity_constraint`` to one of ``positive``, ``non_negative``, ``negative``, or ``non_positive``. +Additional options can be used to modify the behaviour of the linear and +nonlinear solvers: + +- ``cvode_nonlinear_convergence_coef`` specifies the safety factor + used in the nonlinear convergence test. Passed as a parameter to + `CVodeSetNonlinConvCoef + `_. + +- ``cvode_linear_convergence_coef`` specifies the factor by which the + Krylov linear solver’s convergence test constant is reduced from the + nonlinear solver test constant. Passed as a parameter to + `CVodeSetEpsLin + `_. + IMEX-BDF2 --------- @@ -358,7 +372,10 @@ iterations within a given range. +---------------------------+---------------+----------------------------------------------------+ | ksp_type | gmres | PETSc KSP linear solver | +---------------------------+---------------+----------------------------------------------------+ -| pc_type | ilu / bjacobi | PETSc PC preconditioner | +| pc_type | ilu / bjacobi | PETSc PC preconditioner (try hypre in parallel) | ++---------------------------+---------------+----------------------------------------------------+ +| pc_hypre_type | pilut | If ``pc_type = hypre``. | +| | | Hypre preconditioner type: euclid, boomeramg | +---------------------------+---------------+----------------------------------------------------+ | max_nonlinear_iterations | 20 | If exceeded, solve restarts with timestep / 2 | +---------------------------+---------------+----------------------------------------------------+ diff --git a/src/bout++.cxx b/src/bout++.cxx index 79a6738f83..7d7b38b947 100644 --- a/src/bout++.cxx +++ b/src/bout++.cxx @@ -822,7 +822,7 @@ BoutMonitor::BoutMonitor(BoutReal timestep, Options& options) .doc(_("Name of file whose existence triggers a stop")) .withDefault("BOUT.stop"))) {} -int BoutMonitor::call(Solver* solver, BoutReal t, int iter, int NOUT) { +int BoutMonitor::call(Solver* solver, BoutReal t, MAYBE_UNUSED(int iter), int NOUT) { TRACE("BoutMonitor::call({:e}, {:d}, {:d})", t, iter, NOUT); // Increment Solver's iteration counter, and set the global `iteration` diff --git a/src/invert/laplace/impls/petsc/petsc_laplace.cxx b/src/invert/laplace/impls/petsc/petsc_laplace.cxx index 107f3912d9..1df560f288 100644 --- a/src/invert/laplace/impls/petsc/petsc_laplace.cxx +++ b/src/invert/laplace/impls/petsc/petsc_laplace.cxx @@ -148,8 +148,6 @@ LaplacePetsc::LaplacePetsc(Options* opt, const CELL_LOC loc, Mesh* mesh_in, MatCreate(comm, &MatA); MatSetSizes(MatA, localN, localN, size, size); MatSetFromOptions(MatA); - // if (fourth_order) MatMPIAIJSetPreallocation( MatA, 25, PETSC_NULL, 10, PETSC_NULL ); - // else MatMPIAIJSetPreallocation( MatA, 9, PETSC_NULL, 3, PETSC_NULL ); /* Pre allocate memory * nnz denotes an array containing the number of non-zeros in the various rows @@ -854,9 +852,11 @@ FieldPerp LaplacePetsc::solve(const FieldPerp& b, const FieldPerp& x0) { KSPGetConvergedReason(ksp, &reason); if (reason == -3) { // Too many iterations, might be fixed by taking smaller timestep throw BoutIterationFail("petsc_laplace: too many iterations"); - } else if (reason <= 0) { - output << "KSPConvergedReason is " << reason << endl; - throw BoutException("petsc_laplace: inversion failed to converge."); + } + if (reason <= 0) { + throw BoutException( + "petsc_laplace: inversion failed to converge. KSPConvergedReason: {} ({})", + KSPConvergedReasons[reason], reason); } // Add data to FieldPerp Object diff --git a/src/invert/laplace/impls/petsc3damg/petsc3damg.cxx b/src/invert/laplace/impls/petsc3damg/petsc3damg.cxx index 460bec865b..24d61af783 100644 --- a/src/invert/laplace/impls/petsc3damg/petsc3damg.cxx +++ b/src/invert/laplace/impls/petsc3damg/petsc3damg.cxx @@ -265,8 +265,9 @@ Field3D LaplacePetsc3dAmg::solve(const Field3D& b_in, const Field3D& x0) { throw BoutIterationFail("Petsc3dAmg: too many iterations"); } if (reason <= 0) { - output << "KSPConvergedReason is " << reason << "\n"; - throw BoutException("Petsc3dAmg: inversion failed to converge."); + throw BoutException( + "Petsc3dAmg: inversion failed to converge. KSPConvergedReason: {} ({})", + KSPConvergedReasons[reason], reason); } // Create field from result @@ -511,7 +512,9 @@ void LaplacePetsc3dAmg::updateMatrix3D() { // Set the relative and absolute tolerances PCSetType(pc, pctype.c_str()); +#if PETSC_VERSION_LT(3, 18, 0) PCGAMGSetSymGraph(pc, PETSC_TRUE); +#endif } lib.setOptionsFromInputFile(ksp); diff --git a/src/invert/laplacexy/laplacexy.cxx b/src/invert/laplacexy/laplacexy.cxx index 439022e9bb..c1bc616bb1 100644 --- a/src/invert/laplacexy/laplacexy.cxx +++ b/src/invert/laplacexy/laplacexy.cxx @@ -1199,7 +1199,8 @@ const Field2D LaplaceXY::solve(const Field2D& rhs, const Field2D& x0) { KSPGetConvergedReason(ksp, &reason); if (reason <= 0) { - throw BoutException("LaplaceXY failed to converge. Reason {:d}", reason); + throw BoutException("LaplaceXY failed to converge. Reason {} ({:d})", + KSPConvergedReasons[reason], static_cast(reason)); } if (save_performance) { diff --git a/src/invert/laplacexy2/laplacexy2.cxx b/src/invert/laplacexy2/laplacexy2.cxx index dc7a4afa7f..8b5e98d747 100644 --- a/src/invert/laplacexy2/laplacexy2.cxx +++ b/src/invert/laplacexy2/laplacexy2.cxx @@ -383,7 +383,8 @@ Field2D LaplaceXY2::solve(const Field2D& rhs, const Field2D& x0) { KSPGetConvergedReason(ksp, &reason); if (reason <= 0) { - throw BoutException("LaplaceXY2 failed to converge. Reason {:d}", reason); + throw BoutException("LaplaceXY2 failed to converge. Reason {} ({:d})", + KSPConvergedReasons[reason], static_cast(reason)); } // Convert result into a Field2D diff --git a/src/invert/laplacexz/impls/petsc/laplacexz-petsc.cxx b/src/invert/laplacexz/impls/petsc/laplacexz-petsc.cxx index ff28f6d51a..8ddb058ab2 100644 --- a/src/invert/laplacexz/impls/petsc/laplacexz-petsc.cxx +++ b/src/invert/laplacexz/impls/petsc/laplacexz-petsc.cxx @@ -942,7 +942,8 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) { KSPGetConvergedReason(it.ksp, &reason); if (reason <= 0) { - throw BoutException("LaplaceXZ failed to converge. Reason {:d}", reason); + throw BoutException("LaplaceXZ failed to converge. Reason {} ({:d})", + KSPConvergedReasons[reason], static_cast(reason)); } ////////////////////////// diff --git a/src/invert/pardiv/impls/cyclic/pardiv_cyclic.cxx b/src/invert/pardiv/impls/cyclic/pardiv_cyclic.cxx index 32cfea8306..c17e5af64d 100644 --- a/src/invert/pardiv/impls/cyclic/pardiv_cyclic.cxx +++ b/src/invert/pardiv/impls/cyclic/pardiv_cyclic.cxx @@ -40,16 +40,15 @@ #if not BOUT_USE_METRIC_3D +#include #include -#include -#include -#include -#include -#include -#include -#include - +#include +#include +#include +#include +#include #include +#include #include diff --git a/src/invert/pardiv/impls/cyclic/pardiv_cyclic.hxx b/src/invert/pardiv/impls/cyclic/pardiv_cyclic.hxx index 2ce4f1acc9..fdcdab055e 100644 --- a/src/invert/pardiv/impls/cyclic/pardiv_cyclic.hxx +++ b/src/invert/pardiv/impls/cyclic/pardiv_cyclic.hxx @@ -52,9 +52,9 @@ RegisterUnavailableInvertParDiv registerinvertpardivcyclic{ #else -#include "dcomplex.hxx" -#include "utils.hxx" -#include +#include +#include +#include class InvertParDivCR : public InvertParDiv { public: diff --git a/src/mesh/fv_ops.cxx b/src/mesh/fv_ops.cxx index c5fe2fc0ee..0a5d5f9624 100644 --- a/src/mesh/fv_ops.cxx +++ b/src/mesh/fv_ops.cxx @@ -96,7 +96,9 @@ Field3D Div_a_Grad_perp(const Field3D& a, const Field3D& f) { // Result of the Y and Z fluxes Field3D yzresult(0.0, mesh); - // yzresult.allocate(); + if (!fci) { + yzresult.setDirectionY(YDirectionType::Aligned); + } // Y flux diff --git a/src/physics/physicsmodel.cxx b/src/physics/physicsmodel.cxx index f0fbffb0fb..cac4bda5cc 100644 --- a/src/physics/physicsmodel.cxx +++ b/src/physics/physicsmodel.cxx @@ -99,7 +99,6 @@ void PhysicsModel::initialise(Solver* s) { solver = s; bout::experimental::addBuildFlagsToOptions(output_options); - mesh->outputVars(output_options); // Restart option const bool restarting = Options::root()["restart"].withDefault(false); @@ -113,6 +112,8 @@ void PhysicsModel::initialise(Solver* s) { throw BoutException("Couldn't initialise physics model"); } + mesh->outputVars(output_options); + // Post-initialise, which reads restart files // This function can be overridden by the user if (postInit(restarting) != 0) { diff --git a/src/solver/impls/arkode/arkode.cxx b/src/solver/impls/arkode/arkode.cxx index 09aa5b9ea2..13e0dd817e 100644 --- a/src/solver/impls/arkode/arkode.cxx +++ b/src/solver/impls/arkode/arkode.cxx @@ -226,7 +226,7 @@ ArkodeSolver::ArkodeSolver(Options* opts) .withDefault(false)), optimize( (*options)["optimize"].doc("Use ARKode optimal parameters").withDefault(false)), - suncontext(MPI_COMM_WORLD) { + suncontext(static_cast(&BoutComm::get())) { has_constraints = false; // This solver doesn't have constraints // Add diagnostics to output diff --git a/src/solver/impls/cvode/cvode.cxx b/src/solver/impls/cvode/cvode.cxx index 4701d02f69..70eb3b8841 100644 --- a/src/solver/impls/cvode/cvode.cxx +++ b/src/solver/impls/cvode/cvode.cxx @@ -174,7 +174,16 @@ CvodeSolver::CvodeSolver(Options* opts) .doc("Use right preconditioner? Otherwise use left.") .withDefault(false)), use_jacobian((*options)["use_jacobian"].withDefault(false)), - suncontext(MPI_COMM_WORLD) { + cvode_nonlinear_convergence_coef( + (*options)["cvode_nonlinear_convergence_coef"] + .doc("Safety factor used in the nonlinear convergence test") + .withDefault(0.1)), + cvode_linear_convergence_coef( + (*options)["cvode_linear_convergence_coef"] + .doc("Factor by which the Krylov linear solver’s convergence test constant " + "is reduced from the nonlinear solver test constant.") + .withDefault(0.05)), + suncontext(static_cast(&BoutComm::get())) { has_constraints = false; // This solver doesn't have constraints canReset = true; @@ -458,6 +467,10 @@ int CvodeSolver::init() { #endif } + // Set internal tolerance factors + CVodeSetNonlinConvCoef(cvode_mem, cvode_nonlinear_convergence_coef); + CVodeSetEpsLin(cvode_mem, cvode_linear_convergence_coef); + cvode_initialised = true; return 0; diff --git a/src/solver/impls/cvode/cvode.hxx b/src/solver/impls/cvode/cvode.hxx index b2ed83568b..fa8b972bca 100644 --- a/src/solver/impls/cvode/cvode.hxx +++ b/src/solver/impls/cvode/cvode.hxx @@ -123,6 +123,8 @@ private: /// Use right preconditioner? Otherwise use left. bool rightprec; bool use_jacobian; + BoutReal cvode_nonlinear_convergence_coef; + BoutReal cvode_linear_convergence_coef; // Diagnostics from CVODE int nsteps{0}; diff --git a/src/solver/impls/ida/ida.cxx b/src/solver/impls/ida/ida.cxx index cf057ac74c..b008ebf903 100644 --- a/src/solver/impls/ida/ida.cxx +++ b/src/solver/impls/ida/ida.cxx @@ -101,7 +101,7 @@ IdaSolver::IdaSolver(Options* opts) correct_start((*options)["correct_start"] .doc("Correct the initial values") .withDefault(true)), - suncontext(MPI_COMM_WORLD) { + suncontext(static_cast(&BoutComm::get())) { has_constraints = true; // This solver has constraints } diff --git a/src/solver/impls/imex-bdf2/imex-bdf2.cxx b/src/solver/impls/imex-bdf2/imex-bdf2.cxx index aae92ad27d..11954d331a 100644 --- a/src/solver/impls/imex-bdf2/imex-bdf2.cxx +++ b/src/solver/impls/imex-bdf2/imex-bdf2.cxx @@ -676,9 +676,9 @@ void IMEXBDF2::constructSNES(SNES* snesIn) { BoutComm::get(), nlocal, nlocal, // Local sizes PETSC_DETERMINE, PETSC_DETERMINE, // Global sizes 3, // Number of nonzero entries in diagonal portion of local submatrix - PETSC_NULL, + nullptr, 0, // Number of nonzeros per row in off-diagonal portion of local submatrix - PETSC_NULL, &Jmf); + nullptr, &Jmf); #if PETSC_VERSION_GE(3, 4, 0) SNESSetJacobian(*snesIn, Jmf, Jmf, SNESComputeJacobianDefault, this); @@ -1240,7 +1240,7 @@ PetscErrorCode IMEXBDF2::solve_implicit(BoutReal curtime, BoutReal gamma) { if (verbose) { output << "SNES failed to converge with reason " << reason << endl; } - throw BoutException("SNES failed to converge. Reason: {:d}\n", reason); + throw PetscLib::SNESFailure(snesUse); } int its; diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index ad8e99dd38..1b81ca36b6 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -261,11 +261,10 @@ int PetscSolver::init() { CHKERRQ(ierr); #if PETSC_VERSION_GE(3, 7, 0) - ierr = PetscOptionsGetBool(PETSC_NULL, PETSC_NULL, "-interpolate", &interpolate, - PETSC_NULL); + ierr = PetscOptionsGetBool(nullptr, nullptr, "-interpolate", &interpolate, nullptr); CHKERRQ(ierr); #else - ierr = PetscOptionsGetBool(PETSC_NULL, "-interpolate", &interpolate, PETSC_NULL); + ierr = PetscOptionsGetBool(nullptr, "-interpolate", &interpolate, nullptr); CHKERRQ(ierr); #endif @@ -273,21 +272,21 @@ int PetscSolver::init() { // run, if they didn't then use the standard monitor function. TODO: // use PetscFList #if PETSC_VERSION_GE(3, 7, 0) - ierr = PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-output_name", this->output_name, + ierr = PetscOptionsGetString(nullptr, nullptr, "-output_name", this->output_name, sizeof this->output_name, &output_flag); CHKERRQ(ierr); #else - ierr = PetscOptionsGetString(PETSC_NULL, "-output_name", this->output_name, + ierr = PetscOptionsGetString(nullptr, "-output_name", this->output_name, sizeof this->output_name, &output_flag); CHKERRQ(ierr); #endif // If the output_name is not specified then use the standard monitor function - if (output_flag) { - ierr = SNESMonitorSet(snes, PetscSNESMonitor, this, PETSC_NULL); + if (output_flag != 0U) { + ierr = SNESMonitorSet(snes, PetscSNESMonitor, this, nullptr); CHKERRQ(ierr); } else { - ierr = TSMonitorSet(ts, PetscMonitor, this, PETSC_NULL); + ierr = TSMonitorSet(ts, PetscMonitor, this, nullptr); CHKERRQ(ierr); } @@ -399,11 +398,11 @@ int PetscSolver::init() { // Create Jacobian matrix to be used by preconditioner output_info << " Get Jacobian matrix at simtime " << simtime << "\n"; #if PETSC_VERSION_GE(3, 7, 0) - ierr = PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-J_load", load_file, + ierr = PetscOptionsGetString(nullptr, nullptr, "-J_load", load_file, PETSC_MAX_PATH_LEN - 1, &J_load); CHKERRQ(ierr); #else - ierr = PetscOptionsGetString(PETSC_NULL, "-J_load", load_file, PETSC_MAX_PATH_LEN - 1, + ierr = PetscOptionsGetString(nullptr, "-J_load", load_file, PETSC_MAX_PATH_LEN - 1, &J_load); CHKERRQ(ierr); #endif @@ -451,26 +450,26 @@ int PetscSolver::init() { // Get nonzero pattern of J - color_none !!! prealloc = cols * dof * dof; - ierr = MatSeqAIJSetPreallocation(J, prealloc, PETSC_NULL); + ierr = MatSeqAIJSetPreallocation(J, prealloc, nullptr); CHKERRQ(ierr); - ierr = MatMPIAIJSetPreallocation(J, prealloc, PETSC_NULL, prealloc, PETSC_NULL); + ierr = MatMPIAIJSetPreallocation(J, prealloc, nullptr, prealloc, nullptr); CHKERRQ(ierr); - prealloc = - cols; // why nonzeros=295900, allocated nonzeros=2816000/12800000 (*dof*dof), number of mallocs used during MatSetValues calls =256? - ierr = MatSeqBAIJSetPreallocation(J, dof, prealloc, PETSC_NULL); + // why nonzeros=295900, allocated nonzeros=2816000/12800000 (*dof*dof), number of mallocs used during MatSetValues calls =256? + prealloc = cols; + ierr = MatSeqBAIJSetPreallocation(J, dof, prealloc, nullptr); CHKERRQ(ierr); - ierr = MatMPIBAIJSetPreallocation(J, dof, prealloc, PETSC_NULL, prealloc, PETSC_NULL); + ierr = MatMPIBAIJSetPreallocation(J, dof, prealloc, nullptr, prealloc, nullptr); CHKERRQ(ierr); #if PETSC_VERSION_GE(3, 7, 0) - ierr = PetscOptionsHasName(PETSC_NULL, PETSC_NULL, "-J_slowfd", &J_slowfd); + ierr = PetscOptionsHasName(nullptr, nullptr, "-J_slowfd", &J_slowfd); CHKERRQ(ierr); #else - ierr = PetscOptionsHasName(PETSC_NULL, "-J_slowfd", &J_slowfd); + ierr = PetscOptionsHasName(nullptr, "-J_slowfd", &J_slowfd); CHKERRQ(ierr); #endif - if (J_slowfd) { // create Jacobian matrix by slow fd - ierr = SNESSetJacobian(snes, J, J, SNESComputeJacobianDefault, PETSC_NULL); + if (J_slowfd != 0U) { // create Jacobian matrix by slow fd + ierr = SNESSetJacobian(snes, J, J, SNESComputeJacobianDefault, nullptr); CHKERRQ(ierr); output_info << "SNESComputeJacobian J by slow fd...\n"; @@ -524,15 +523,15 @@ int PetscSolver::init() { // Write J in binary for study - see ~petsc/src/mat/examples/tests/ex124.c #if PETSC_VERSION_GE(3, 7, 0) - ierr = PetscOptionsHasName(PETSC_NULL, PETSC_NULL, "-J_write", &J_write); + ierr = PetscOptionsHasName(nullptr, nullptr, "-J_write", &J_write); CHKERRQ(ierr); #else - ierr = PetscOptionsHasName(PETSC_NULL, "-J_write", &J_write); + ierr = PetscOptionsHasName(nullptr, "-J_write", &J_write); CHKERRQ(ierr); #endif - if (J_write) { - PetscViewer viewer; - output_info << "\n[" << rank << "] Test TSComputeRHSJacobian() ...\n"; + if (J_write != 0U) { + PetscViewer viewer = nullptr; + output_info.write("\n[{:d}] Test TSComputeRHSJacobian() ...\n", rank); #if PETSC_VERSION_GE(3, 5, 0) ierr = TSComputeRHSJacobian(ts, simtime, u, J, J); CHKERRQ(ierr); @@ -542,22 +541,14 @@ int PetscSolver::init() { CHKERRQ(ierr); #endif - ierr = PetscSynchronizedPrintf(comm, "[{:d}] TSComputeRHSJacobian is done\n", rank); - CHKERRQ(ierr); + output.write("[{:d}] TSComputeRHSJacobian is done\n", rank); -#if PETSC_VERSION_GE(3, 5, 0) - ierr = PetscSynchronizedFlush(comm, PETSC_STDOUT); - CHKERRQ(ierr); -#else - ierr = PetscSynchronizedFlush(comm); - CHKERRQ(ierr); -#endif - if (J_slowfd) { - output_info << "[" << rank << "] writing J in binary to data/Jrhs_dense.dat...\n"; + if (J_slowfd != 0U) { + output_info.write("[{:d}] writing J in binary to data/Jrhs_dense.dat...\n", rank); ierr = PetscViewerBinaryOpen(comm, "data/Jrhs_dense.dat", FILE_MODE_WRITE, &viewer); CHKERRQ(ierr); } else { - output_info << "[" << rank << "] writing J in binary to data/Jrhs_sparse.dat...\n"; + output_info.write("[{:d}] writing J in binary to data/Jrhs_sparse.dat...\n", rank); ierr = PetscViewerBinaryOpen(comm, "data/Jrhs_sparse.dat", FILE_MODE_WRITE, &viewer); CHKERRQ(ierr); @@ -579,8 +570,6 @@ int PetscSolver::init() { **************************************************************************/ PetscErrorCode PetscSolver::run() { - PetscErrorCode ierr; - FILE* fp = nullptr; // Set when the next call to monitor is desired next_output = simtime + getOutputTimestep(); @@ -592,23 +581,18 @@ PetscErrorCode PetscSolver::run() { bout_snes_time = bout::globals::mpi->MPI_Wtime(); } - ierr = TSSolve(ts, u); - CHKERRQ(ierr); + CHKERRQ(TSSolve(ts, u)); // Gawd, everything is a hack - if (this->output_flag) { - ierr = PetscFOpen(PETSC_COMM_WORLD, this->output_name, "w", &fp); - CHKERRQ(ierr); - ierr = PetscFPrintf(PETSC_COMM_WORLD, fp, - "SNES Iteration, KSP Iterations, Wall Time, Norm\n"); - CHKERRQ(ierr); + if ((this->output_flag != 0U) and (BoutComm::rank() == 0)) { + Output petsc_info(output_name); + // Don't write to stdout + petsc_info.disable(); + petsc_info.write("SNES Iteration, KSP Iterations, Wall Time, Norm\n"); for (const auto& info : snes_list) { - ierr = PetscFPrintf(PETSC_COMM_WORLD, fp, "{:d}, {:d}, {:e}, {:e}\n", info.it, - info.linear_its, info.time, info.norm); - CHKERRQ(ierr); + petsc_info.write("{:d}, {:d}, {:e}, {:e}\n", info.it, info.linear_its, info.time, + info.norm); } - ierr = PetscFClose(PETSC_COMM_WORLD, fp); - CHKERRQ(ierr); } PetscFunctionReturn(0); @@ -876,7 +860,7 @@ PetscErrorCode PhysicsSNESApply(SNES snes, Vec x) { Mat A, B; PetscFunctionBegin; - ierr = SNESGetJacobian(snes, &A, &B, PETSC_NULL, PETSC_NULL); + ierr = SNESGetJacobian(snes, &A, &B, nullptr, nullptr); CHKERRQ(ierr); #if PETSC_VERSION_GE(3, 5, 0) ierr = SNESComputeJacobian(snes, x, A, B); @@ -890,7 +874,7 @@ PetscErrorCode PhysicsSNESApply(SNES snes, Vec x) { CHKERRQ(ierr); ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr); - ierr = SNESGetFunction(snes, &F, PETSC_NULL, PETSC_NULL); + ierr = SNESGetFunction(snes, &F, nullptr, nullptr); CHKERRQ(ierr); ierr = SNESComputeFunction(snes, x, F); CHKERRQ(ierr); @@ -913,7 +897,7 @@ PetscErrorCode PhysicsSNESApply(SNES snes, Vec x) { << ", F \\cdot Fout " << dot << " "; #if PETSC_VERSION_GE(3, 5, 0) Vec func; - ierr = SNESGetFunction(snes, &func, PETSC_NULL, PETSC_NULL); + ierr = SNESGetFunction(snes, &func, nullptr, nullptr); CHKERRQ(ierr); ierr = VecNorm(func, NORM_2, &fnorm); CHKERRQ(ierr); @@ -962,7 +946,7 @@ PetscErrorCode PetscMonitor(TS ts, PetscInt UNUSED(step), PetscReal t, Vec X, vo ierr = TSGetMaxTime(ts, &tfinal); CHKERRQ(ierr); #else - ierr = TSGetDuration(ts, PETSC_NULL, &tfinal); + ierr = TSGetDuration(ts, nullptr, &tfinal); CHKERRQ(ierr); #endif diff --git a/src/solver/impls/slepc/slepc.cxx b/src/solver/impls/slepc/slepc.cxx index 3f6b7a91bf..f90afe717e 100644 --- a/src/solver/impls/slepc/slepc.cxx +++ b/src/solver/impls/slepc/slepc.cxx @@ -162,11 +162,13 @@ SlepcSolver::SlepcSolver(Options* options) { .withDefault(0); tol = options_ref["tol"].doc("SLEPc tolerance").withDefault(1.0e-6); - maxIt = options_ref["maxIt"].doc("Maximum iterations").withDefault(PETSC_DEFAULT); + maxIt = options_ref["maxIt"] + .doc("Maximum iterations") + .withDefault(static_cast(PETSC_DEFAULT)); mpd = options_ref["mpd"] .doc("Maximum dimension allowed for the projected problem") - .withDefault(PETSC_DEFAULT); + .withDefault(static_cast(PETSC_DEFAULT)); ddtMode = options_ref["ddtMode"].withDefault(true); diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 69280bd6c2..2143aad96a 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -108,6 +108,10 @@ SNESSolver::SNESSolver(Options* opts) (*options)["pc_type"] .doc("Preconditioner type. By default lets PETSc decide (ilu or bjacobi)") .withDefault("default")), + pc_hypre_type((*options)["pc_hypre_type"] + .doc("hypre preconditioner type: euclid, pilut, parasails, " + "boomeramg, ams, ads") + .withDefault("pilut")), line_search_type((*options)["line_search_type"] .doc("Line search type: basic, bt, l2, cp, nleqerr") .withDefault("default")), @@ -146,6 +150,7 @@ int SNESSolver::init() { int ierr; // Vectors + output_info.write("Creating vector\n"); ierr = VecCreate(BoutComm::get(), &snes_x); CHKERRQ(ierr); ierr = VecSetSizes(snes_x, nlocal, PETSC_DECIDE); @@ -167,6 +172,7 @@ int SNESSolver::init() { } // Nonlinear solver interface (SNES) + output_info.write("Create SNES\n"); SNESCreate(BoutComm::get(), &snes); // Set the callback function @@ -234,10 +240,26 @@ int SNESSolver::init() { std::vector o_nnz(localN); // Set values for most points + const int ncells_x = (mesh->LocalNx > 1) ? 2 : 0; + const int ncells_y = (mesh->LocalNy > 1) ? 2 : 0; + const int ncells_z = (mesh->LocalNz > 1) ? 2 : 0; + + const auto star_pattern = (1 + ncells_x + ncells_y) * (n3d + n2d) + ncells_z * n3d; + + // Offsets. Start with the central cell + std::vector> xyoffsets{{0, 0}}; + if (ncells_x != 0) { + // Stencil includes points in X + xyoffsets.push_back({-1, 0}); + xyoffsets.push_back({1, 0}); + } + if (ncells_y != 0) { + // Stencil includes points in Y + xyoffsets.push_back({0, -1}); + xyoffsets.push_back({0, 1}); + } - const auto star_pattern_2d = 5 * (n3d + n2d); - const auto star_pattern_3d = 7 * n3d + 5 * n2d; - const auto star_pattern = (mesh->LocalNz > 1) ? star_pattern_3d : star_pattern_2d; + output_info.write("Star pattern: {} non-zero entries\n", star_pattern); for (int i = 0; i < localN; i++) { // Non-zero elements on this processor d_nnz[i] = star_pattern; @@ -246,148 +268,154 @@ int SNESSolver::init() { } // X boundaries - if (mesh->firstX()) { - // Lower X boundary - for (int y = mesh->ystart; y <= mesh->yend; y++) { - for (int z = 0; z < mesh->LocalNz; z++) { - int localIndex = ROUND(index(mesh->xstart, y, z)); - ASSERT2((localIndex >= 0) && (localIndex < localN)); - const int num_fields = (z == 0) ? n2d + n3d : n3d; - for (int i = 0; i < num_fields; i++) { - d_nnz[localIndex + i] -= (n3d + n2d); + if (ncells_x != 0) { + if (mesh->firstX()) { + // Lower X boundary + for (int y = mesh->ystart; y <= mesh->yend; y++) { + for (int z = 0; z < mesh->LocalNz; z++) { + int localIndex = ROUND(index(mesh->xstart, y, z)); + ASSERT2((localIndex >= 0) && (localIndex < localN)); + const int num_fields = (z == 0) ? n2d + n3d : n3d; + for (int i = 0; i < num_fields; i++) { + d_nnz[localIndex + i] -= (n3d + n2d); + } } } - } - } else { - // On another processor - for (int y = mesh->ystart; y <= mesh->yend; y++) { - for (int z = 0; z < mesh->LocalNz; z++) { - int localIndex = ROUND(index(mesh->xstart, y, z)); - ASSERT2((localIndex >= 0) && (localIndex < localN)); - const int num_fields = (z == 0) ? n2d + n3d : n3d; - for (int i = 0; i < num_fields; i++) { - d_nnz[localIndex + i] -= (n3d + n2d); - o_nnz[localIndex + i] += (n3d + n2d); + } else { + // On another processor + for (int y = mesh->ystart; y <= mesh->yend; y++) { + for (int z = 0; z < mesh->LocalNz; z++) { + int localIndex = ROUND(index(mesh->xstart, y, z)); + ASSERT2((localIndex >= 0) && (localIndex < localN)); + const int num_fields = (z == 0) ? n2d + n3d : n3d; + for (int i = 0; i < num_fields; i++) { + d_nnz[localIndex + i] -= (n3d + n2d); + o_nnz[localIndex + i] += (n3d + n2d); + } } } } - } - - if (mesh->lastX()) { - // Upper X boundary - for (int y = mesh->ystart; y <= mesh->yend; y++) { - for (int z = 0; z < mesh->LocalNz; z++) { - int localIndex = ROUND(index(mesh->xend, y, z)); - ASSERT2((localIndex >= 0) && (localIndex < localN)); - const int num_fields = (z == 0) ? n2d + n3d : n3d; - for (int i = 0; i < num_fields; i++) { - d_nnz[localIndex + i] -= (n3d + n2d); + if (mesh->lastX()) { + // Upper X boundary + for (int y = mesh->ystart; y <= mesh->yend; y++) { + for (int z = 0; z < mesh->LocalNz; z++) { + int localIndex = ROUND(index(mesh->xend, y, z)); + ASSERT2((localIndex >= 0) && (localIndex < localN)); + const int num_fields = (z == 0) ? n2d + n3d : n3d; + for (int i = 0; i < num_fields; i++) { + d_nnz[localIndex + i] -= (n3d + n2d); + } } } - } - } else { - // On another processor - for (int y = mesh->ystart; y <= mesh->yend; y++) { - for (int z = 0; z < mesh->LocalNz; z++) { - int localIndex = ROUND(index(mesh->xend, y, z)); - ASSERT2((localIndex >= 0) && (localIndex < localN)); - const int num_fields = (z == 0) ? n2d + n3d : n3d; - for (int i = 0; i < num_fields; i++) { - d_nnz[localIndex + i] -= (n3d + n2d); - o_nnz[localIndex + i] += (n3d + n2d); + } else { + // On another processor + for (int y = mesh->ystart; y <= mesh->yend; y++) { + for (int z = 0; z < mesh->LocalNz; z++) { + int localIndex = ROUND(index(mesh->xend, y, z)); + ASSERT2((localIndex >= 0) && (localIndex < localN)); + const int num_fields = (z == 0) ? n2d + n3d : n3d; + for (int i = 0; i < num_fields; i++) { + d_nnz[localIndex + i] -= (n3d + n2d); + o_nnz[localIndex + i] += (n3d + n2d); + } } } } } // Y boundaries - - for (int x = mesh->xstart; x <= mesh->xend; x++) { - // Default to no boundary - // NOTE: This assumes that communications in Y are to other - // processors. If Y is communicated with this processor (e.g. NYPE=1) - // then this will result in PETSc warnings about out of range allocations - - // z = 0 case - int localIndex = ROUND(index(x, mesh->ystart, 0)); - ASSERT2(localIndex >= 0); - - // All 2D and 3D fields - for (int i = 0; i < n2d + n3d; i++) { - o_nnz[localIndex + i] += (n3d + n2d); - } - - for (int z = 1; z < mesh->LocalNz; z++) { - localIndex = ROUND(index(x, mesh->ystart, z)); - - // Only 3D fields - for (int i = 0; i < n3d; i++) { + if (ncells_y != 0) { + for (int x = mesh->xstart; x <= mesh->xend; x++) { + // Default to no boundary + // NOTE: This assumes that communications in Y are to other + // processors. If Y is communicated with this processor (e.g. NYPE=1) + // then this will result in PETSc warnings about out of range allocations + + // z = 0 case + int localIndex = ROUND(index(x, mesh->ystart, 0)); + ASSERT2(localIndex >= 0); + + // All 2D and 3D fields + for (int i = 0; i < n2d + n3d; i++) { o_nnz[localIndex + i] += (n3d + n2d); + d_nnz[localIndex + i] -= (n3d + n2d); } - } - // z = 0 case - localIndex = ROUND(index(x, mesh->yend, 0)); - // All 2D and 3D fields - for (int i = 0; i < n2d + n3d; i++) { - o_nnz[localIndex + i] += (n3d + n2d); - } + for (int z = 1; z < mesh->LocalNz; z++) { + localIndex = ROUND(index(x, mesh->ystart, z)); - for (int z = 1; z < mesh->LocalNz; z++) { - localIndex = ROUND(index(x, mesh->yend, z)); + // Only 3D fields + for (int i = 0; i < n3d; i++) { + o_nnz[localIndex + i] += (n3d + n2d); + d_nnz[localIndex + i] -= (n3d + n2d); + } + } - // Only 3D fields - for (int i = 0; i < n3d; i++) { + // z = 0 case + localIndex = ROUND(index(x, mesh->yend, 0)); + // All 2D and 3D fields + for (int i = 0; i < n2d + n3d; i++) { o_nnz[localIndex + i] += (n3d + n2d); + d_nnz[localIndex + i] -= (n3d + n2d); } - } - } - for (RangeIterator it = mesh->iterateBndryLowerY(); !it.isDone(); it++) { - // A boundary, so no communication + for (int z = 1; z < mesh->LocalNz; z++) { + localIndex = ROUND(index(x, mesh->yend, z)); - // z = 0 case - int localIndex = ROUND(index(it.ind, mesh->ystart, 0)); - if (localIndex < 0) { - // This can occur because it.ind includes values in x boundary e.g. x=0 - continue; - } - // All 2D and 3D fields - for (int i = 0; i < n2d + n3d; i++) { - o_nnz[localIndex + i] -= (n3d + n2d); + // Only 3D fields + for (int i = 0; i < n3d; i++) { + o_nnz[localIndex + i] += (n3d + n2d); + d_nnz[localIndex + i] -= (n3d + n2d); + } + } } - for (int z = 1; z < mesh->LocalNz; z++) { - int localIndex = ROUND(index(it.ind, mesh->ystart, z)); + for (RangeIterator it = mesh->iterateBndryLowerY(); !it.isDone(); it++) { + // A boundary, so no communication - // Only 3D fields - for (int i = 0; i < n3d; i++) { + // z = 0 case + int localIndex = ROUND(index(it.ind, mesh->ystart, 0)); + if (localIndex < 0) { + // This can occur because it.ind includes values in x boundary e.g. x=0 + continue; + } + // All 2D and 3D fields + for (int i = 0; i < n2d + n3d; i++) { o_nnz[localIndex + i] -= (n3d + n2d); } - } - } - for (RangeIterator it = mesh->iterateBndryUpperY(); !it.isDone(); it++) { - // A boundary, so no communication + for (int z = 1; z < mesh->LocalNz; z++) { + int localIndex = ROUND(index(it.ind, mesh->ystart, z)); - // z = 0 case - int localIndex = ROUND(index(it.ind, mesh->yend, 0)); - if (localIndex < 0) { - continue; // Out of domain + // Only 3D fields + for (int i = 0; i < n3d; i++) { + o_nnz[localIndex + i] -= (n3d + n2d); + } + } } - // All 2D and 3D fields - for (int i = 0; i < n2d + n3d; i++) { - o_nnz[localIndex + i] -= (n3d + n2d); - } + for (RangeIterator it = mesh->iterateBndryUpperY(); !it.isDone(); it++) { + // A boundary, so no communication - for (int z = 1; z < mesh->LocalNz; z++) { - int localIndex = ROUND(index(it.ind, mesh->yend, z)); + // z = 0 case + int localIndex = ROUND(index(it.ind, mesh->yend, 0)); + if (localIndex < 0) { + continue; // Out of domain + } - // Only 3D fields - for (int i = 0; i < n3d; i++) { + // All 2D and 3D fields + for (int i = 0; i < n2d + n3d; i++) { o_nnz[localIndex + i] -= (n3d + n2d); } + + for (int z = 1; z < mesh->LocalNz; z++) { + int localIndex = ROUND(index(it.ind, mesh->yend, z)); + + // Only 3D fields + for (int i = 0; i < n3d; i++) { + o_nnz[localIndex + i] -= (n3d + n2d); + } + } } } @@ -395,15 +423,19 @@ int SNESSolver::init() { // Pre-allocate MatMPIAIJSetPreallocation(Jmf, 0, d_nnz.data(), 0, o_nnz.data()); + MatSeqAIJSetPreallocation(Jmf, 0, d_nnz.data()); MatSetUp(Jmf); - MatSetOption(Jmf, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); + MatSetOption(Jmf, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE); // Determine which row/columns of the matrix are locally owned int Istart, Iend; MatGetOwnershipRange(Jmf, &Istart, &Iend); // Convert local into global indices - index += Istart; + // Note: Not in the boundary cells, to keep -1 values + for (const auto& i : mesh->getRegion3D("RGN_NOBNDRY")) { + index[i] += Istart; + } // Now communicate to fill guard cells mesh->communicate(index); @@ -413,11 +445,6 @@ int SNESSolver::init() { output_progress.write("Marking non-zero Jacobian entries\n"); - // Offsets for a 5-point pattern - constexpr std::size_t stencil_size = 5; - const std::array xoffset = {0, -1, 1, 0, 0}; - const std::array yoffset = {0, 0, 0, -1, 1}; - PetscScalar val = 1.0; for (int x = mesh->xstart; x <= mesh->xend; x++) { @@ -430,9 +457,9 @@ int SNESSolver::init() { PetscInt row = ind0 + i; // Loop through each point in the 5-point stencil - for (std::size_t c = 0; c < stencil_size; c++) { - int xi = x + xoffset[c]; - int yi = y + yoffset[c]; + for (const auto& xyoffset : xyoffsets) { + int xi = x + xyoffset.first; + int yi = y + xyoffset.second; if ((xi < 0) || (yi < 0) || (xi >= mesh->LocalNx) || (yi >= mesh->LocalNy)) { @@ -448,8 +475,8 @@ int SNESSolver::init() { // Depends on all variables on this cell for (int j = 0; j < n2d; j++) { PetscInt col = ind2 + j; - - MatSetValues(Jmf, 1, &row, 1, &col, &val, INSERT_VALUES); + ierr = MatSetValues(Jmf, 1, &row, 1, &col, &val, INSERT_VALUES); + CHKERRQ(ierr); } } } @@ -468,13 +495,14 @@ int SNESSolver::init() { // Depends on 2D fields for (int j = 0; j < n2d; j++) { PetscInt col = ind0 + j; - MatSetValues(Jmf, 1, &row, 1, &col, &val, INSERT_VALUES); + ierr = MatSetValues(Jmf, 1, &row, 1, &col, &val, INSERT_VALUES); + CHKERRQ(ierr); } - // 5 point star pattern - for (std::size_t c = 0; c < stencil_size; c++) { - int xi = x + xoffset[c]; - int yi = y + yoffset[c]; + // Star pattern + for (const auto& xyoffset : xyoffsets) { + int xi = x + xyoffset.first; + int yi = y + xyoffset.second; if ((xi < 0) || (yi < 0) || (xi >= mesh->LocalNx) || (yi >= mesh->LocalNy)) { @@ -493,7 +521,12 @@ int SNESSolver::init() { // 3D fields on this cell for (int j = 0; j < n3d; j++) { PetscInt col = ind2 + j; - MatSetValues(Jmf, 1, &row, 1, &col, &val, INSERT_VALUES); + ierr = MatSetValues(Jmf, 1, &row, 1, &col, &val, INSERT_VALUES); + if (ierr != 0) { + output.write("ERROR: {} : ({}, {}) -> ({}, {}) : {} -> {}\n", row, x, + y, xi, yi, ind2, ind2 + n3d - 1); + } + CHKERRQ(ierr); } } @@ -509,7 +542,8 @@ int SNESSolver::init() { } for (int j = 0; j < n3d; j++) { PetscInt col = ind2 + j; - MatSetValues(Jmf, 1, &row, 1, &col, &val, INSERT_VALUES); + ierr = MatSetValues(Jmf, 1, &row, 1, &col, &val, INSERT_VALUES); + CHKERRQ(ierr); } int zm = (z - 1 + nz) % nz; @@ -519,7 +553,8 @@ int SNESSolver::init() { } for (int j = 0; j < n3d; j++) { PetscInt col = ind2 + j; - MatSetValues(Jmf, 1, &row, 1, &col, &val, INSERT_VALUES); + ierr = MatSetValues(Jmf, 1, &row, 1, &col, &val, INSERT_VALUES); + CHKERRQ(ierr); } } } @@ -564,9 +599,9 @@ int SNESSolver::init() { BoutComm::get(), nlocal, nlocal, // Local sizes PETSC_DETERMINE, PETSC_DETERMINE, // Global sizes 3, // Number of nonzero entries in diagonal portion of local submatrix - PETSC_NULL, + nullptr, 0, // Number of nonzeros per row in off-diagonal portion of local submatrix - PETSC_NULL, &Jmf); + nullptr, &Jmf); #if PETSC_VERSION_GE(3, 4, 0) SNESSetJacobian(snes, Jmf, Jmf, SNESComputeJacobianDefault, this); #else @@ -633,6 +668,15 @@ int SNESSolver::init() { // Set PC type from input if (pc_type != "default") { PCSetType(pc, pc_type.c_str()); + + if (pc_type == "hypre") { +#if PETSC_HAVE_HYPRE + // Set the type of hypre preconditioner + PCHYPRESetType(pc, pc_hypre_type.c_str()); +#else + throw BoutException("PETSc was not configured with Hypre."); +#endif + } } } @@ -726,7 +770,7 @@ int SNESSolver::run() { // Print diagnostics to help identify source of the problem output.write("\n======== SNES failed =========\n"); - output.write("\nReturn code: {}, reason: {}\n", ierr, reason); + output.write("\nReturn code: {}, reason: {}\n", ierr, static_cast(reason)); for (const auto& f : f2d) { output.write( "{} : ({} -> {}), ddt: ({} -> {})\n", f.name, @@ -748,6 +792,12 @@ int SNESSolver::run() { // Restore state VecCopy(x0, snes_x); + // Recalculate the Jacobian + if (saved_jacobian_lag == 0) { + SNESGetLagJacobian(snes, &saved_jacobian_lag); + SNESSetLagJacobian(snes, 1); + } + // Check lock state PetscInt lock_state; VecLockGet(snes_x, &lock_state); @@ -816,7 +866,7 @@ int SNESSolver::run() { output.print("\r"); // Carriage return for printing to screen output.write("Time: {}, timestep: {}, nl iter: {}, lin iter: {}, reason: {}", - simtime, timestep, nl_its, lin_its, reason); + simtime, timestep, nl_its, lin_its, static_cast(reason)); if (snes_failures > 0) { output.write(", SNES failures: {}", snes_failures); snes_failures = 0; diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx index c00f4b2581..2021402cd7 100644 --- a/src/solver/impls/snes/snes.hxx +++ b/src/solver/impls/snes/snes.hxx @@ -125,6 +125,7 @@ private: bool kspsetinitialguessnonzero; ///< Set initial guess to non-zero int maxl; ///< Maximum linear iterations std::string pc_type; ///< Preconditioner type + std::string pc_hypre_type; ///< Hypre preconditioner type std::string line_search_type; ///< Line search type bool matrix_free; ///< Use matrix free Jacobian int lag_jacobian; ///< Re-use Jacobian diff --git a/src/sys/boutcomm.cxx b/src/sys/boutcomm.cxx index ecaeb866bb..2d78a24076 100644 --- a/src/sys/boutcomm.cxx +++ b/src/sys/boutcomm.cxx @@ -25,7 +25,7 @@ void BoutComm::setComm(MPI_Comm c) { hasBeenSet = true; } -MPI_Comm BoutComm::getComm() { +MPI_Comm& BoutComm::getComm() { if (comm == MPI_COMM_NULL) { // No communicator set. Initialise MPI MPI_Init(pargc, pargv); @@ -39,7 +39,7 @@ MPI_Comm BoutComm::getComm() { bool BoutComm::isSet() { return hasBeenSet; } // Static functions below. Must use getInstance() -MPI_Comm BoutComm::get() { return getInstance()->getComm(); } +MPI_Comm& BoutComm::get() { return getInstance()->getComm(); } void BoutComm::setArgs(int& c, char**& v) { getInstance()->pargc = &c; diff --git a/src/sys/output.cxx b/src/sys/output.cxx index 0ab1d2e0c8..c07cd53474 100644 --- a/src/sys/output.cxx +++ b/src/sys/output.cxx @@ -109,5 +109,3 @@ ConditionalOutput output_progress(Output::getInstance()); ConditionalOutput output_error(Output::getInstance()); ConditionalOutput output_verbose(Output::getInstance(), false); ConditionalOutput output(Output::getInstance()); - -#undef bout_vsnprint_pre diff --git a/src/sys/petsclib.cxx b/src/sys/petsclib.cxx index 91db0c44ea..b99c51320c 100644 --- a/src/sys/petsclib.cxx +++ b/src/sys/petsclib.cxx @@ -5,9 +5,10 @@ #include "bout/boutcomm.hxx" #include "bout/openmpwrap.hxx" #include "bout/options.hxx" +#include #include -#include +#include "petscsnes.h" // Define all the static member variables int PetscLib::count = 0; @@ -23,16 +24,18 @@ PetscLib::PetscLib(Options* opt) { if (count == 0) { // Initialise PETSc + // Load global PETSc options from the [petsc] section of the input + // Note: This should be before PetscInitialize so that some options + // can modify initialization e.g. -log_view. + setPetscOptions(Options::root()["petsc"], ""); + output << "Initialising PETSc\n"; PETSC_COMM_WORLD = BoutComm::getInstance()->getComm(); - PetscInitialize(pargc, pargv, PETSC_NULL, help); + PetscInitialize(pargc, pargv, nullptr, help); PetscPopSignalHandler(); PetscLogEventRegister("Total BOUT++", 0, &USER_EVENT); PetscLogEventBegin(USER_EVENT, 0, 0, 0, 0); - - // Load global PETSc options from the [petsc] section of the input - setPetscOptions(Options::root()["petsc"], ""); } if (opt != nullptr and opt->isSection()) { @@ -66,27 +69,15 @@ PetscLib::~PetscLib() { } void PetscLib::setOptionsFromInputFile(KSP& ksp) { - auto ierr = KSPSetOptionsPrefix(ksp, options_prefix.c_str()); - if (ierr) { - throw BoutException("KSPSetOptionsPrefix failed with error {}", ierr); - } + assertIerr(KSPSetOptionsPrefix(ksp, options_prefix.c_str()), "KSPSetOptionsPrefix"); - ierr = KSPSetFromOptions(ksp); - if (ierr) { - throw BoutException("KSPSetFromOptions failed with error {}", ierr); - } + assertIerr(KSPSetFromOptions(ksp), "KSPSetFromOptions"); } void PetscLib::setOptionsFromInputFile(SNES& snes) { - auto ierr = SNESSetOptionsPrefix(snes, options_prefix.c_str()); - if (ierr) { - throw BoutException("SNESSetOptionsPrefix failed with error %i", ierr); - } + BOUT_DO_PETSC(SNESSetOptionsPrefix(snes, options_prefix.c_str())); - ierr = SNESSetFromOptions(snes); - if (ierr) { - throw BoutException("SNESSetFromOptions failed with error %i", ierr); - } + BOUT_DO_PETSC(SNESSetFromOptions(snes)); } void PetscLib::cleanup() { @@ -133,4 +124,17 @@ void PetscLib::setPetscOptions(Options& options, const std::string& prefix) { } } } + +BoutException PetscLib::SNESFailure(SNES& snes) { + SNESConvergedReason reason; + BOUT_DO_PETSC(SNESGetConvergedReason(snes, &reason)); +#if PETSC_VERSION_GE(3, 15, 0) + const char* message; + BOUT_DO_PETSC(SNESGetConvergedReasonString(snes, &message)); +#else + const char* message{""}; +#endif + return BoutException("SNES failed to converge. Reason: {} ({:d})", message, + static_cast(reason)); +} #endif // BOUT_HAS_PETSC diff --git a/src/sys/slepclib.cxx b/src/sys/slepclib.cxx index 1dc89bc047..79fdfa9bbc 100644 --- a/src/sys/slepclib.cxx +++ b/src/sys/slepclib.cxx @@ -18,7 +18,7 @@ SlepcLib::SlepcLib() { // Initialise SLEPc output << "Initialising SLEPc\n"; - SlepcInitialize(pargc, pargv, PETSC_NULL, help); + SlepcInitialize(pargc, pargv, nullptr, help); PetscLogEventRegister("Total BOUT++", 0, &USER_EVENT); PetscLogEventBegin(USER_EVENT, 0, 0, 0, 0); } diff --git a/tests/integrated/test-bout-override-default-option/test-bout-override-default-option.cxx b/tests/integrated/test-bout-override-default-option/test-bout-override-default-option.cxx index 24cb6bc0b2..4c50579c47 100644 --- a/tests/integrated/test-bout-override-default-option/test-bout-override-default-option.cxx +++ b/tests/integrated/test-bout-override-default-option/test-bout-override-default-option.cxx @@ -1,4 +1,4 @@ -#include +#include // Use an integrated test for what is effectively a unit test of the // BOUT_OVERRIDE_DEFAULT_OPTION() macro because the functionality relies on the state of diff --git a/tests/integrated/test-laplace-hypre3d/test-laplace3d.cxx b/tests/integrated/test-laplace-hypre3d/test-laplace3d.cxx index 69524d44da..64b5dec1d4 100644 --- a/tests/integrated/test-laplace-hypre3d/test-laplace3d.cxx +++ b/tests/integrated/test-laplace-hypre3d/test-laplace3d.cxx @@ -125,8 +125,13 @@ int main(int argc, char** argv) { Field3D error = rhs_check - rhs; BoutReal error_max = max(abs(error), true); - SAVE_ONCE(f, rhs, rhs_check, error, error_max); - bout::globals::dump.write(); + Options dump; + dump["f"] = f; + dump["rhs"] = rhs; + dump["rhs_check"] = rhs_check; + dump["error"] = error; + dump["error_max"] = error_max; + bout::writeDefaultOutputFile(dump); laplace_solver.reset(nullptr); BoutFinalise(); diff --git a/tests/integrated/test-laplace-petsc3d/data_circular_core/BOUT.inp b/tests/integrated/test-laplace-petsc3d/data_circular_core/BOUT.inp index ae7cfaa5b4..6474b2604b 100644 --- a/tests/integrated/test-laplace-petsc3d/data_circular_core/BOUT.inp +++ b/tests/integrated/test-laplace-petsc3d/data_circular_core/BOUT.inp @@ -76,6 +76,7 @@ atol = 1e-13 [laplace:petsc] mg_levels_ksp_max_it = 3 +mg_levels_pc_type = sor [input] transform_from_field_aligned = false diff --git a/tests/integrated/test-laplacexy/test-laplacexy.cxx b/tests/integrated/test-laplacexy/test-laplacexy.cxx index 62449a76ee..3142d359b1 100644 --- a/tests/integrated/test-laplacexy/test-laplacexy.cxx +++ b/tests/integrated/test-laplacexy/test-laplacexy.cxx @@ -30,27 +30,20 @@ #include #include -using bout::globals::dump; using bout::globals::mesh; int main(int argc, char** argv) { BoutInitialise(argc, argv); - auto coords = mesh->getCoordinates(); - - auto& opt = Options::root(); - - LaplaceXY laplacexy; - - bool include_y_derivs = opt["laplacexy"]["include_y_derivs"]; + auto* coords = mesh->getCoordinates(); // Solving equations of the form // Div(A Grad_perp(f)) + B*f = rhs // A*Laplace_perp(f) + Grad_perp(A).Grad_perp(f) + B*f = rhs - Field2D f, a, b, sol; - Field2D error, absolute_error; //Absolute value of relative error: abs((f - sol)/f) - BoutReal max_error; //Output of test + Field2D f; + Field2D a; + Field2D b; initial_profile("f", f); initial_profile("a", a); @@ -63,42 +56,45 @@ int main(int argc, char** argv) { //////////////////////////////////////////////////////////////////////////////////////// - Field2D rhs, rhs_check; + Field2D rhs; + const bool include_y_derivs = Options::root()["laplacexy"]["include_y_derivs"]; if (include_y_derivs) { rhs = a * Laplace_perp(f) + Grad_perp(a) * Grad_perp(f) + b * f; } else { rhs = a * Delp2(f, CELL_DEFAULT, false) + coords->g11 * DDX(a) * DDX(f) + b * f; } + LaplaceXY laplacexy; laplacexy.setCoefs(a, b); - sol = laplacexy.solve(rhs, 0.); - error = (f - sol) / f; - absolute_error = f - sol; - max_error = max(abs(absolute_error), true); + Field2D solution = laplacexy.solve(rhs, 0.); + Field2D relative_error = (f - solution) / f; + Field2D absolute_error = f - solution; + BoutReal max_error = max(abs(absolute_error), true); - output << "Magnitude of maximum absolute error is " << max_error << endl; + output.write("Magnitude of maximum absolute error is {}\n", max_error); - mesh->communicate(sol); + mesh->communicate(solution); + Field2D rhs_check; if (include_y_derivs) { - rhs_check = a * Laplace_perp(sol) + Grad_perp(a) * Grad_perp(sol) + b * sol; - } else { rhs_check = - a * Delp2(sol, CELL_DEFAULT, false) + coords->g11 * DDX(a) * DDX(sol) + b * sol; + a * Laplace_perp(solution) + Grad_perp(a) * Grad_perp(solution) + b * solution; + } else { + rhs_check = a * Delp2(solution, CELL_DEFAULT, false) + + coords->g11 * DDX(a) * DDX(solution) + b * solution; } - dump.add(a, "a"); - dump.add(b, "b"); - dump.add(f, "f"); - dump.add(sol, "sol"); - dump.add(error, "error"); - dump.add(absolute_error, "absolute_error"); - dump.add(max_error, "max_error"); - dump.add(rhs, "rhs"); - dump.add(rhs_check, "rhs_check"); - - dump.write(); - dump.close(); + Options dump; + dump["a"] = a; + dump["b"] = b; + dump["f"] = f; + dump["sol"] = solution; + dump["relative_error"] = relative_error; + dump["absolute_error"] = absolute_error; + dump["max_error"] = max_error; + dump["rhs"] = rhs; + dump["rhs_check"] = rhs_check; + bout::writeDefaultOutputFile(dump); MPI_Barrier(BoutComm::get()); // Wait for all processors to write data diff --git a/tests/integrated/test-laplacexy2-hypre/test-laplacexy.cxx b/tests/integrated/test-laplacexy2-hypre/test-laplacexy.cxx index 77fbdb197c..0ae2cbf312 100644 --- a/tests/integrated/test-laplacexy2-hypre/test-laplacexy.cxx +++ b/tests/integrated/test-laplacexy2-hypre/test-laplacexy.cxx @@ -59,8 +59,8 @@ int main(int argc, char** argv) { Field2D guess = 0.0; Field2D sol = laplacexy.solve(rhs, guess); Field2D error = (f - sol) / f; - Field2D absolute_error = - abs(f - sol); // Absolute value of relative error: abs((f - sol)/f) + // Absolute value of relative error: abs((f - sol)/f) + Field2D absolute_error = abs(f - sol); BoutReal max_error = max(absolute_error, true); output << "Magnitude of maximum absolute error is " << max_error << endl; @@ -68,20 +68,17 @@ int main(int argc, char** argv) { sol.getMesh()->communicate(sol); Field2D rhs_check = Laplace_perpXY(a, sol); - using bout::globals::dump; - - dump.add(a, "a"); - dump.add(b, "b"); - dump.add(f, "f"); - dump.add(sol, "sol"); - dump.add(error, "error"); - dump.add(absolute_error, "absolute_error"); - dump.add(max_error, "max_error"); - dump.add(rhs, "rhs"); - dump.add(rhs_check, "rhs_check"); - - dump.write(); - dump.close(); + Options dump; + dump["a"] = a; + dump["b"] = b; + dump["f"] = f; + dump["sol"] = sol; + dump["error"] = error; + dump["absolute_error"] = absolute_error; + dump["max_error"] = max_error; + dump["rhs"] = rhs; + dump["rhs_check"] = rhs_check; + bout::writeDefaultOutputFile(dump); MPI_Barrier(BoutComm::get()); // Wait for all processors to write data diff --git a/tests/integrated/test-squash/runtest b/tests/integrated/test-squash/runtest index 09bf2f1d14..692d561c59 100755 --- a/tests/integrated/test-squash/runtest +++ b/tests/integrated/test-squash/runtest @@ -7,6 +7,7 @@ import numpy as np from boututils.run_wrapper import launch_safe, shell_safe, build_and_log import argparse import re +import os.path # requires: all_tests @@ -81,7 +82,7 @@ def verify(f1, f2): parser = argparse.ArgumentParser(description="Test the bout-squashoutput wrapper") parser.add_argument( - "executable", help="Path to bout-squashoutput", default="../../../bin" + "executable", help="Path to bout-squashoutput", default="../../../bin", nargs="?" ) args = parser.parse_args() @@ -89,6 +90,9 @@ build_and_log("Squash test") bout_squashoutput = args.executable + "/bout-squashoutput" +if not os.path.exists(bout_squashoutput): + bout_squashoutput = "bout-squashoutput" + print("Run once to get normal data") timed_shell_safe("./squash -q -q -q nout=2") timed_shell_safe("mv data/BOUT.dmp.0.nc f1.nc") diff --git a/tests/unit/invert/laplace/test_laplace_petsc3damg.cxx b/tests/unit/invert/laplace/test_laplace_petsc3damg.cxx index 832ccc8e6b..b2831808a4 100644 --- a/tests/unit/invert/laplace/test_laplace_petsc3damg.cxx +++ b/tests/unit/invert/laplace/test_laplace_petsc3damg.cxx @@ -99,7 +99,7 @@ class Petsc3dAmgTest public: WithQuietOutput info{output_info}, warn{output_warn}, progress{output_progress}, all{output}; - Petsc3dAmgTest() : FakeMeshFixture(), solver(getOptions(GetParam())) { + Petsc3dAmgTest() : solver(&getOptions(GetParam())) { PetscErrorPrintf = PetscErrorPrintfNone; int nx = mesh->GlobalNx, ny = mesh->GlobalNy, nz = mesh->GlobalNz; static_cast(bout::globals::mesh) @@ -139,20 +139,23 @@ class Petsc3dAmgTest ForwardOperator forward; private: - static Options* getOptions(std::tuple param) { - Options* options = Options::getRoot()->getSection("laplace"); - (*options)["type"] = "petsc3damg"; - (*options)["inner_boundary_flags"] = + static Options& getOptions(std::tuple param) { + auto& options = Options::root()["laplace"]; + options["type"] = "petsc3damg"; + options["inner_boundary_flags"] = (std::get<0>(param) ? INVERT_AC_GRAD : 0) + INVERT_RHS; - (*options)["outer_boundary_flags"] = + options["outer_boundary_flags"] = (std::get<1>(param) ? INVERT_AC_GRAD : 0) + INVERT_RHS; - (*options)["lower_boundary_flags"] = + options["lower_boundary_flags"] = (std::get<2>(param) ? INVERT_AC_GRAD : 0) + INVERT_RHS; - (*options)["upper_boundary_flags"] = + options["upper_boundary_flags"] = (std::get<3>(param) ? INVERT_AC_GRAD : 0) + INVERT_RHS; - (*options)["fourth_order"] = false; - (*options)["atol"] = tol / 30; // Need to specify smaller than desired tolerance to - (*options)["rtol"] = tol / 30; // ensure it is satisfied for every element. + options["fourth_order"] = false; + options["atol"] = tol / 30; // Need to specify smaller than desired tolerance to + options["rtol"] = tol / 30; // ensure it is satisfied for every element. + auto& petsc_options = options["petsc"]; + petsc_options["mg_levels_ksp_max_it"] = 4; + petsc_options["mg_levels_pc_type"] = "sor"; return options; } }; diff --git a/tools/pylib/_boutpp_build/CMakeLists.txt b/tools/pylib/_boutpp_build/CMakeLists.txt index 7ceae04f0e..6b88986a28 100644 --- a/tools/pylib/_boutpp_build/CMakeLists.txt +++ b/tools/pylib/_boutpp_build/CMakeLists.txt @@ -3,7 +3,7 @@ macro(bout_python_maybe_error VAR NAME) if (NOT ${VAR}) set(_error_msg "${NAME} is required for the Python interface") - if (NOT "${BOUT_ENABLE_PYTHON}" STREQUAL "MAYBE") + if (NOT "${BOUT_ENABLE_PYTHON}" STREQUAL "AUTO") message(FATAL_ERROR ${_error_msg}) else() message(WARNING ${_error_msg}) @@ -49,7 +49,7 @@ if (NOT ${PYTHON_WORKING} EQUAL 0) endif() # No errors? We can build the interface! -if ("${BOUT_ENABLE_PYTHON}" STREQUAL "MAYBE") +if ("${BOUT_ENABLE_PYTHON}" STREQUAL "AUTO") set(BOUT_ENABLE_PYTHON ON PARENT_SCOPE) endif() diff --git a/tools/pylib/_boutpp_build/backend.py b/tools/pylib/_boutpp_build/backend.py index 65a15c913c..b4c9b306b8 100644 --- a/tools/pylib/_boutpp_build/backend.py +++ b/tools/pylib/_boutpp_build/backend.py @@ -6,6 +6,7 @@ import subprocess # corelib import re # corelib import pathlib # corelib +import contextlib # corelib try: import packaging.tags # packaging @@ -27,17 +28,28 @@ def getversion(): """ global version if version is None: + with contextlib.suppress(KeyError): + version = os.environ["BOUT_PRETEND_VERSION"] + return version + _bout_previous_version = "v5.0.0" - _bout_next_version = "5.1.0" + _bout_next_version = "v5.1.0" try: - tmp = run2(f"git describe --tags --match={_bout_previous_version}").strip() - tmp = re.sub(f"{_bout_previous_version}-", f"{_bout_next_version}.dev", tmp) - if useLocalVersion: - tmp = re.sub("-", "+", tmp) - else: - tmp = re.sub("-.*", "", tmp) - version = tmp + try: + version = run2("git describe --exact-match --tags HEAD").strip() + except subprocess.CalledProcessError: + tmp = run2( + f"git describe --tags --match={_bout_previous_version}" + ).strip() + tmp = re.sub( + f"{_bout_previous_version}-", f"{_bout_next_version}.dev", tmp + ) + if useLocalVersion: + tmp = re.sub("-", "+", tmp) + else: + tmp = re.sub("-.*", "", tmp) + version = tmp with open("_version.txt", "w") as f: f.write(version + "\n") except subprocess.CalledProcessError: @@ -145,6 +157,7 @@ def build_sdist(sdist_directory, config_settings=None): print(config_settings, sdist_directory) enable_gz = True enable_xz = False + external = {"fmt", "mpark.variant"} if config_settings is not None: global useLocalVersion, pkgname for k, v in config_settings.items(): @@ -159,6 +172,10 @@ def build_sdist(sdist_directory, config_settings=None): enable_xz = True else: raise ValueError(f"unknown option {v} for {k}") + if k == "dist": + enable_xz = True + pkgname = "BOUT++" + external.add("googletest") if k == "useLocalVersion": useLocalVersion = False if k == "nightly": @@ -168,7 +185,7 @@ def build_sdist(sdist_directory, config_settings=None): fname = f"{prefix}.tar" run(f"git archive HEAD --prefix {prefix}/ -o {sdist_directory}/{fname}") _, tmp = tempfile.mkstemp(suffix=".tar") - for ext in "fmt", "mpark.variant": + for ext in sorted(external): run( f"git archive --remote=externalpackages/{ext} HEAD --prefix {prefix}/externalpackages/{ext}/ --format=tar > {tmp}" ) @@ -268,22 +285,70 @@ def parse(fn): def nightly(): + """ + Build the python sdist for upload to boutpp-nightly. + """ return build_sdist(os.getcwd() + "/dist/", dict(nightly=True)) def sdist(): + """ + Build the python sdist + """ return build_sdist(os.getcwd() + "/dist/") def wheel(): + """ + Build the python binary wheel + """ return build_wheel(os.getcwd() + "/dist/") +def dist(): + """ + Build an archive for BOUT++ release + """ + return build_sdist(os.getcwd(), config_settings=dict(dist=True)) + + +def help(): + """ + Print this help + """ + table = [] + for k, v in todos.items(): + try: + doc = v.__doc__.strip() + doc = " : " + doc + except: + doc = "" + table.append((k, doc)) + maxkey = max([len(k) for k, _ in table]) + fmt = f" %-{maxkey}s%s" + print(f"{sys.argv[0]} [command] [command]") + for row in table: + print(fmt % row) + + +todos = dict( + nightly=nightly, + sdist=sdist, + wheel=wheel, + dist=dist, + version=lambda: print(getversion()), + help=help, +) +todos.update( + { + "--help": help, + "-?": help, + "-h": help, + } +) + if __name__ == "__main__": import sys - todos = dict( - nightly=nightly, sdist=sdist, wheel=wheel, version=lambda: print(getversion()) - ) for todo in sys.argv[1:]: todos[todo]() diff --git a/tools/pylib/_boutpp_build/boutpp.pyx.jinja b/tools/pylib/_boutpp_build/boutpp.pyx.jinja index 0d8cef173a..657e2f28c1 100644 --- a/tools/pylib/_boutpp_build/boutpp.pyx.jinja +++ b/tools/pylib/_boutpp_build/boutpp.pyx.jinja @@ -1500,7 +1500,7 @@ def interp_to(Field3D f3d,location): def setOption(name, value, source="PyInterface", force=False): """ Set an option in the global Options tree. Prefer - `Options.set` to avoid unexpected results if several Option + ``Options.set`` to avoid unexpected results if several Option roots are avalaible. Parameters