diff --git a/.clang-tidy b/.clang-tidy index 556b8738f8..4caf530aee 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -9,6 +9,9 @@ CheckOptions: # otherwise this breaks `ASSERT` macros! - key: readability-simplify-boolean-expr.IgnoreMacros value: 'true' + + - key: misc-include-cleaner.IgnoreHeaders + value: 'petsc.*\.h;mpi\.h' --- Disabled checks and reasons: diff --git a/CMakeLists.txt b/CMakeLists.txt index 5c49e434e4..cc368a3b32 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -235,11 +235,13 @@ set(BOUT_SOURCES ./src/invert/laplace/impls/hypre3d/hypre3d_laplace.cxx ./src/invert/laplace/impls/hypre3d/hypre3d_laplace.hxx ./src/invert/laplace/invert_laplace.cxx + ./src/invert/laplacexy/impls/hypre/laplacexy-hypre.cxx + ./src/invert/laplacexy/impls/hypre/laplacexy-hypre.hxx + ./src/invert/laplacexy/impls/petsc/laplacexy-petsc.cxx + ./src/invert/laplacexy/impls/petsc/laplacexy-petsc.hxx + ./src/invert/laplacexy/impls/petsc2/laplacexy-petsc2.cxx + ./src/invert/laplacexy/impls/petsc2/laplacexy-petsc2.hxx ./src/invert/laplacexy/laplacexy.cxx - ./include/bout/invert/laplacexy2.hxx - ./src/invert/laplacexy2/laplacexy2.cxx - ./include/bout/invert/laplacexy2_hypre.hxx - ./src/invert/laplacexy2/laplacexy2_hypre.cxx ./src/invert/laplacexz/impls/cyclic/laplacexz-cyclic.cxx ./src/invert/laplacexz/impls/cyclic/laplacexz-cyclic.hxx ./src/invert/laplacexz/impls/petsc/laplacexz-petsc.cxx diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 849f10e85f..f2444e6bf6 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -30,7 +30,6 @@ add_subdirectory(invertable_operator) add_subdirectory(laplacexy/alfven-wave) add_subdirectory(laplacexy/laplace_perp) add_subdirectory(laplacexy/simple) -add_subdirectory(laplacexy/simple-hypre) add_subdirectory(monitor-newapi) add_subdirectory(orszag-tang) add_subdirectory(preconditioning/wave) diff --git a/examples/dalf3/dalf3.cxx b/examples/dalf3/dalf3.cxx index 493e8f5e45..46a9a4286d 100644 --- a/examples/dalf3/dalf3.cxx +++ b/examples/dalf3/dalf3.cxx @@ -300,7 +300,7 @@ class DALF3 : public PhysicsModel { // LaplaceXY for n=0 solve if (split_n0) { // Create an XY solver for n=0 component - laplacexy = std::make_unique(mesh); + laplacexy = LaplaceXY::create(mesh); phi2D = 0.0; // Starting guess } diff --git a/examples/elm-pb-outerloop/elm_pb_outerloop.cxx b/examples/elm-pb-outerloop/elm_pb_outerloop.cxx index 8e84901806..4ff8e90dcb 100644 --- a/examples/elm-pb-outerloop/elm_pb_outerloop.cxx +++ b/examples/elm-pb-outerloop/elm_pb_outerloop.cxx @@ -28,33 +28,26 @@ /*******************************************************************************/ +#include "bout/build_defines.hxx" +#include "bout/options.hxx" #include #include #include +#include #include #include #include #include #include #include -#include -#include - -#include - -#include -#include #include +#include // Defines BOUT_FOR_RAJA #include #include +#include +#include -#include // Defines BOUT_FOR_RAJA - -#if BOUT_HAS_HYPRE -#include -#endif - -#include +#include CELL_LOC loc = CELL_CENTRE; @@ -91,6 +84,10 @@ BOUT_OVERRIDE_DEFAULT_OPTION("phi:bndry_target", "neumann"); BOUT_OVERRIDE_DEFAULT_OPTION("phi:bndry_xin", "none"); BOUT_OVERRIDE_DEFAULT_OPTION("phi:bndry_xout", "none"); +#if BOUT_HAS_HYPRE +BOUT_OVERRIDE_DEFAULT_OPTION("laplacexy:type", "hypre"); +#endif + /// 3-field ELM simulation class ELMpb : public PhysicsModel { private: @@ -242,11 +239,7 @@ class ELMpb : public PhysicsModel { bool split_n0; // Solve the n=0 component of potential -#if BOUT_HAS_HYPRE - std::unique_ptr laplacexy{nullptr}; // Laplacian solver in X-Y (n=0) -#else std::unique_ptr laplacexy{nullptr}; // Laplacian solver in X-Y (n=0) -#endif Field2D phi2D; // Axisymmetric phi @@ -567,13 +560,11 @@ class ELMpb : public PhysicsModel { split_n0 = options["split_n0"] .doc("Solve zonal (n=0) component of potential using LaplaceXY?") .withDefault(false); + if (split_n0) { // Create an XY solver for n=0 component -#if BOUT_HAS_HYPRE - laplacexy = bout::utils::make_unique(mesh); -#else - laplacexy = bout::utils::make_unique(mesh); -#endif + laplacexy = LaplaceXY::create(mesh); + // Set coefficients for Boussinesq solve laplacexy->setCoefs(1.0, 0.0); phi2D = 0.0; // Starting guess diff --git a/examples/elm-pb/elm_pb.cxx b/examples/elm-pb/elm_pb.cxx index d981796200..fe18a756ae 100644 --- a/examples/elm-pb/elm_pb.cxx +++ b/examples/elm-pb/elm_pb.cxx @@ -22,10 +22,6 @@ #include -#if BOUT_HAS_HYPRE -#include -#endif - #include CELL_LOC loc = CELL_CENTRE; @@ -36,6 +32,10 @@ BOUT_OVERRIDE_DEFAULT_OPTION("phi:bndry_target", "neumann"); BOUT_OVERRIDE_DEFAULT_OPTION("phi:bndry_xin", "none"); BOUT_OVERRIDE_DEFAULT_OPTION("phi:bndry_xout", "none"); +#if BOUT_HAS_HYPRE +BOUT_OVERRIDE_DEFAULT_OPTION("laplacexy:type", "hypre"); +#endif + /// 3-field ELM simulation class ELMpb : public PhysicsModel { private: @@ -188,11 +188,7 @@ class ELMpb : public PhysicsModel { bool split_n0; // Solve the n=0 component of potential -#if BOUT_HAS_HYPRE - std::unique_ptr laplacexy{nullptr}; // Laplacian solver in X-Y (n=0) -#else std::unique_ptr laplacexy{nullptr}; // Laplacian solver in X-Y (n=0) -#endif Field2D phi2D; // Axisymmetric phi @@ -525,11 +521,8 @@ class ELMpb : public PhysicsModel { .withDefault(false); if (split_n0) { // Create an XY solver for n=0 component -#if BOUT_HAS_HYPRE - laplacexy = bout::utils::make_unique(mesh); -#else - laplacexy = bout::utils::make_unique(mesh); -#endif + laplacexy = LaplaceXY::create(mesh); + // Set coefficients for Boussinesq solve laplacexy->setCoefs(1.0, 0.0); phi2D = 0.0; // Starting guess diff --git a/examples/laplacexy/alfven-wave/alfven.cxx b/examples/laplacexy/alfven-wave/alfven.cxx index a4248afcf7..031931c1c4 100644 --- a/examples/laplacexy/alfven-wave/alfven.cxx +++ b/examples/laplacexy/alfven-wave/alfven.cxx @@ -1,10 +1,11 @@ - #include #include #include #include #include +#include + /// Fundamental constants const BoutReal PI = 3.14159265; const BoutReal e0 = 8.854e-12; // Permittivity of free space @@ -32,7 +33,7 @@ class Alfven : public PhysicsModel { bool laplace_perp; // Use Laplace_perp or Delp2? bool split_n0; // Split solve into n=0 and n~=0? - LaplaceXY* laplacexy; // Laplacian solver in X-Y (n=0) + std::unique_ptr laplacexy{nullptr}; // Laplacian solver in X-Y (n=0) bool newXZsolver; std::unique_ptr phiSolver; // Old Laplacian in X-Z @@ -80,7 +81,7 @@ class Alfven : public PhysicsModel { if (split_n0) { // Create an XY solver for n=0 component - laplacexy = new LaplaceXY(mesh); + laplacexy = LaplaceXY::create(mesh); phi2D = 0.0; // Starting guess } diff --git a/examples/laplacexy/laplace_perp/test.cxx b/examples/laplacexy/laplace_perp/test.cxx index 75577bb8b7..e3d9528d45 100644 --- a/examples/laplacexy/laplace_perp/test.cxx +++ b/examples/laplacexy/laplace_perp/test.cxx @@ -56,10 +56,10 @@ int main(int argc, char** argv) { FieldFactory::get()->create2D("input_field", Options::getRoot(), mesh); // Create a LaplaceXY solver - LaplaceXY laplacexy{mesh}; + auto laplacexy = LaplaceXY::create(mesh); // Solve, using 0.0 as starting guess - Field2D solved = laplacexy.solve(input, 0.0); + Field2D solved = laplacexy->solve(input, 0.0); // Need to communicate guard cells mesh->communicate(solved); diff --git a/examples/laplacexy/simple-hypre/.gitignore b/examples/laplacexy/simple-hypre/.gitignore deleted file mode 100644 index eeede075f9..0000000000 --- a/examples/laplacexy/simple-hypre/.gitignore +++ /dev/null @@ -1 +0,0 @@ -test-laplacexy diff --git a/examples/laplacexy/simple-hypre/CMakeLists.txt b/examples/laplacexy/simple-hypre/CMakeLists.txt deleted file mode 100644 index 788ebb1776..0000000000 --- a/examples/laplacexy/simple-hypre/CMakeLists.txt +++ /dev/null @@ -1,15 +0,0 @@ -cmake_minimum_required(VERSION 3.13) - -project(test_laplacexy_hypre LANGUAGES CXX C) - -if (NOT TARGET bout++::bout++) - find_package(bout++ REQUIRED) -endif() - -bout_add_example(test_laplacexy_hypre - SOURCES test-laplacexy-hypre.cxx - REQUIRES BOUT_HAS_HYPRE) - -if(BOUT_HAS_CUDA) - set_source_files_properties(test-laplacexy-hypre.cxx PROPERTIES LANGUAGE CUDA ) -endif() diff --git a/examples/laplacexy/simple-hypre/data/BOUT.inp b/examples/laplacexy/simple-hypre/data/BOUT.inp deleted file mode 100644 index 5da492168d..0000000000 --- a/examples/laplacexy/simple-hypre/data/BOUT.inp +++ /dev/null @@ -1,37 +0,0 @@ -# -# Simple test of the LaplaceXY solver. -# -# Inverts a given function (rhs), writing "rhs" and solution "x" to file -# - -[mesh] - -# Mesh sizes -nx = 20 -ny = 32 -nz = 1 - -# mesh spacing -dx = 1.0 -dy = 1.0 -dz = 1.0 - -[laplacexy] - -ksptype = gmres # Iterative solver type - not used with Hypre interface -pctype = jacobi # Preconditioner. "jacobi", "bjacobi" and "sor" usually good -# On one processor"lu" uses direct solver - -atol = 1e-12 # type: BoutReal, doc: Relative tolerance for Hypre solver -core_bndry_dirichlet = false # type: bool -hypre_print_level = 3 # type: int, doc: Verbosity for Hypre solver. Integer from 0 (silent) to 4 (most verbose). -#hypre_solver_type = gmres # type: HYPRE_SOLVER_TYPE, doc: Type of solver to use when solving Hypre system. Possible values are: gmres, bicgstab -hypre_solver_type = bicgstab # type: HYPRE_SOLVER_TYPE, doc: Type of solver to use when solving Hypre system. Possible values are: gmres, bicgstab -include_y_derivs = true # type: bool, doc: Include Y derivatives in operator to invert? -maxits = 2000 # type: int, doc: Maximum iterations for Hypre solver -print_timing = true # type: bool, doc: Print extra timing information for LaplaceXY2Hypre -rtol = 1e-07 # type: BoutReal, doc: Relative tolerance for Hypre solver -y_bndry_dirichlet = false # type: bool - -# Function to be inverted -rhs = sin(2*pi*x)*sin(y) diff --git a/examples/laplacexy/simple-hypre/test-laplacexy-hypre.cxx b/examples/laplacexy/simple-hypre/test-laplacexy-hypre.cxx deleted file mode 100644 index af2367e1a0..0000000000 --- a/examples/laplacexy/simple-hypre/test-laplacexy-hypre.cxx +++ /dev/null @@ -1,30 +0,0 @@ -#include -#include -#include - -int main(int argc, char** argv) { - BoutInitialise(argc, argv); - { - /// Create a LaplaceXY object - LaplaceXY2Hypre laplacexy(bout::globals::mesh); - - /// Generate rhs function - Field2D rhs = FieldFactory::get()->create2D("laplacexy:rhs", Options::getRoot(), - bout::globals::mesh); - - /// Solution - Field2D solution = 0.0; - - solution = laplacexy.solve(rhs, solution); - - Options dump; - dump["rhs"] = rhs; - dump["x"] = solution; - bout::writeDefaultOutputFile(dump); - } - BoutFinalise(); -#if BOUT_HAS_CUDA - cudaDeviceReset(); -#endif - return 0; -} diff --git a/examples/laplacexy/simple/CMakeLists.txt b/examples/laplacexy/simple/CMakeLists.txt index ec5e7d2492..7859a08259 100644 --- a/examples/laplacexy/simple/CMakeLists.txt +++ b/examples/laplacexy/simple/CMakeLists.txt @@ -6,4 +6,7 @@ if (NOT TARGET bout++::bout++) find_package(bout++ REQUIRED) endif() -bout_add_example(laplacexy-simple SOURCES test-laplacexy.cxx) +bout_add_example(laplacexy-simple + SOURCES test-laplacexy.cxx + DATA_DIRS data hypre +) diff --git a/examples/laplacexy/simple/README.md b/examples/laplacexy/simple/README.md index ce7cbebf00..491d00d4b3 100644 --- a/examples/laplacexy/simple/README.md +++ b/examples/laplacexy/simple/README.md @@ -8,7 +8,7 @@ and preconditioners. See the "ksptype" and "pctype" settings in BOUT.inp Run with - $ ./test-laplacexy -ksp_monitor + $ ./test-laplacexy -laplacexy:petsc:ksp_monitor which should print the KSP norms from PETSc: @@ -31,3 +31,8 @@ which should print the KSP norms from PETSc: 16 KSP Residual norm 4.309526296050e-01 17 KSP Residual norm 1.115269396077e-01 18 KSP Residual norm 4.334487475743e-13 + +HYPRE +----- + +Use the `hypre` directory to use the HYPRE preconditioner instead of PETSc. diff --git a/examples/laplacexy/simple/hypre/BOUT.inp b/examples/laplacexy/simple/hypre/BOUT.inp new file mode 100644 index 0000000000..338e902d4e --- /dev/null +++ b/examples/laplacexy/simple/hypre/BOUT.inp @@ -0,0 +1,23 @@ +# +# Simple test of the LaplaceXY solver using HYPRE +# +# Inverts a given function (rhs), writing "rhs" and solution "x" to file +# + +[mesh] + +# Mesh sizes +nx = 20 +ny = 32 +nz = 1 + +# mesh spacing +dx = 1.0 +dy = 1.0 +dz = 1.0 + +[laplacexy] +type = hypre + +# Function to be inverted +rhs = sin(2*pi*x)*sin(y) diff --git a/examples/laplacexy/simple/test-laplacexy.cxx b/examples/laplacexy/simple/test-laplacexy.cxx index c033eb68e3..49a93aafdb 100644 --- a/examples/laplacexy/simple/test-laplacexy.cxx +++ b/examples/laplacexy/simple/test-laplacexy.cxx @@ -7,20 +7,18 @@ int main(int argc, char** argv) { BoutInitialise(argc, argv); /// Create a LaplaceXY object - LaplaceXY laplacexy(bout::globals::mesh); + auto laplacexy = LaplaceXY::create(bout::globals::mesh); /// Generate rhs function Field2D rhs = FieldFactory::get()->create2D("laplacexy:rhs", Options::getRoot(), bout::globals::mesh); /// Solution - Field2D x = 0.0; - - x = laplacexy.solve(rhs, x); + Field2D result = laplacexy->solve(rhs, 0.0); Options dump; dump["rhs"] = rhs; - dump["x"] = x; + dump["result"] = result; bout::writeDefaultOutputFile(dump); BoutFinalise(); diff --git a/include/bout/invert/laplacexy.hxx b/include/bout/invert/laplacexy.hxx index 5a4dd46512..51605b9db5 100644 --- a/include/bout/invert/laplacexy.hxx +++ b/include/bout/invert/laplacexy.hxx @@ -1,18 +1,18 @@ /************************************************************************** - * Laplacian solver in 2D (X-Y) + * Perpendicular Laplacian inversion in X-Y * - * Equation solved is: + * Equation solved is: * * Div( A * Grad_perp(x) ) + B*x = b * * Intended for use in solving n = 0 component of potential * from inversion of vorticity equation - * + * ************************************************************************** - * Copyright 2015 B.Dudson + * Copyright 2015-2025 BOUT++ Team * * Contact: Ben Dudson, bd512@york.ac.uk - * + * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify @@ -30,196 +30,67 @@ * **************************************************************************/ -#ifndef BOUT_LAPLACE_XY_H -#define BOUT_LAPLACE_XY_H - -#include "bout/build_defines.hxx" +#ifndef BOUT_LAPLACEXY_H +#define BOUT_LAPLACEXY_H -#if !BOUT_HAS_PETSC -// If no PETSc - -#warning LaplaceXY requires PETSc. No LaplaceXY available +#include +#include +#include +#include -#include +#include -class Field2D; -class Mesh; +class LaplaceXY; class Options; class Solver; -/*! - * Create a dummy class so that code will compile - * without PETSc, but will throw an exception if - * LaplaceXY is used. - */ -class LaplaceXY { +class LaplaceXYFactory + : public Factory { public: - LaplaceXY(Mesh* UNUSED(m) = nullptr, Options* UNUSED(opt) = nullptr, - const CELL_LOC UNUSED(loc) = CELL_CENTRE) { - throw BoutException("LaplaceXY requires PETSc. No LaplaceXY available"); - } - void setCoefs(const Field2D& UNUSED(A), const Field2D& UNUSED(B)) {} - const Field2D solve(const Field2D& UNUSED(rhs), const Field2D& UNUSED(x0)) { - throw BoutException("LaplaceXY requires PETSc. No LaplaceXY available"); + static constexpr auto type_name = "LaplaceXY"; + static constexpr auto section_name = "laplacexy"; + static constexpr auto option_name = "type"; + static constexpr auto default_type = "petsc"; + + ReturnType create(Mesh* mesh = nullptr, Options* options = nullptr, + CELL_LOC loc = CELL_CENTRE) const { + return Factory::create(getType(options), mesh, optionsOrDefaultSection(options), loc); } - void savePerformance(Solver&, std::string) { - throw BoutException("LaplaceXY requires PETSc. No LaplaceXY available"); + ReturnType create(const std::string& type, Options* options) const { + return Factory::create(type, nullptr, options, CELL_CENTRE); } + + static void ensureRegistered(); }; -#else // BOUT_HAS_PETSC +template +using RegisterLaplaceXY = LaplaceXYFactory::RegisterInFactory; -#include "bout/solver.hxx" -#include "bout/utils.hxx" -#include -#include -#include - -class Options; -class Solver; +using RegisterUnavailableLaplaceXY = LaplaceXYFactory::RegisterUnavailableInFactory; class LaplaceXY { public: - /*! - * Constructor - */ - LaplaceXY(Mesh* m = nullptr, Options* opt = nullptr, const CELL_LOC loc = CELL_CENTRE); - /*! - * Destructor - */ - ~LaplaceXY(); - - /*! - * Set coefficients (A, B) in equation: - * Div( A * Grad_perp(x) ) + B*x = b - */ - void setCoefs(const Field2D& A, const Field2D& B); - - /*! - * Solve Laplacian in X-Y - * - * Inputs - * ====== - * - * rhs - The field to be inverted. This must be allocated - * and contain valid data. - * x0 - Initial guess at the solution. If this is unallocated - * then an initial guess of zero will be used. - * - * Returns - * ======= - * - * The solution as a Field2D. On failure an exception will be raised - * - */ - const Field2D solve(const Field2D& rhs, const Field2D& x0); - - /*! - * Preconditioner function - * This is called by PETSc via a static function. - * and should not be called by external users - */ - int precon(Vec x, Vec y); - - /*! - * If this method is called, save some performance monitoring information - */ - void savePerformance(Solver& solver, const std::string& name = ""); - -private: - PetscLib lib; ///< Requires PETSc library - Mat MatA; ///< Matrix to be inverted - Vec xs, bs; ///< Solution and RHS vectors - KSP ksp; ///< Krylov Subspace solver - PC pc; ///< Preconditioner + LaplaceXY() = default; + LaplaceXY(const LaplaceXY&) = default; + LaplaceXY(LaplaceXY&&) = delete; + LaplaceXY& operator=(const LaplaceXY&) = default; + LaplaceXY& operator=(LaplaceXY&&) = delete; + virtual ~LaplaceXY() = default; - Mesh* localmesh; ///< The mesh this operates on, provides metrics and communication + virtual void setCoefs(const Field2D& A, const Field2D& B) = 0; - /// default prefix for writing performance logging variables - std::string default_prefix; + virtual Field2D solve(const Field2D& b, const Field2D& x0) = 0; - // Preconditioner - int xstart, xend; - int nloc, nsys; - Matrix acoef, bcoef, ccoef, xvals, bvals; - std::unique_ptr> cr; ///< Tridiagonal solver - - // Use finite volume or finite difference discretization - bool finite_volume{true}; - - // Y derivatives - bool include_y_derivs; // Include Y derivative terms? - - // Boundary conditions - bool x_inner_dirichlet; // Dirichlet on inner X boundary? - bool x_outer_dirichlet; // Dirichlet on outer X boundary? - std::string y_bndry{"neumann"}; // Boundary condition for y-boundary - - // Location of the rhs and solution - CELL_LOC location; - - /*! - * Number of grid points on this processor - */ - int localSize(); - - /*! - * Return the communicator for XY - */ - MPI_Comm communicator(); - - /*! - * Return the global index of a local (x,y) coordinate - * including guard cells. - * Boundary cells have a global index of -1 - * - * To do this, a Field2D (indexXY) is used to store - * the index as a floating point number which is then rounded - * to an integer. Guard cells are filled by communication - * so no additional logic is needed in Mesh. - */ - int globalIndex(int x, int y); - Field2D indexXY; ///< Global index (integer stored as BoutReal) - - // Save performance information? - bool save_performance = false; - - // Running average of number of iterations taken for solve in each output timestep - BoutReal average_iterations = 0.; - - // Variable to store the final result of average_iterations, since output is - // written after all other monitors have been called, and average_iterations - // must be reset in the monitor - BoutReal output_average_iterations = 0.; - - // Running total of number of calls to the solver in each output timestep - int n_calls = 0; - - // Utility methods - void setPreallocationFiniteVolume(PetscInt* d_nnz, PetscInt* o_nnz); - void setPreallocationFiniteDifference(PetscInt* d_nnz, PetscInt* o_nnz); - void setMatrixElementsFiniteVolume(const Field2D& A, const Field2D& B); - void setMatrixElementsFiniteDifference(const Field2D& A, const Field2D& B); - void solveFiniteVolume(const Field2D& x0); - void solveFiniteDifference(const Field2D& x0); - - // Monitor class used to reset performance-monitoring variables for a new - // output timestep - friend class LaplaceXYMonitor; - class LaplaceXYMonitor : public Monitor { - public: - LaplaceXYMonitor(LaplaceXY& owner) : laplacexy(owner) {} - - int call(Solver* /*solver*/, BoutReal /*time*/, int /*iter*/, int /*nout*/) override; - void outputVars(Options& output_options, const std::string& time_dimension) override; - - private: - // LaplaceXY object that this monitor belongs to - LaplaceXY& laplacexy; - }; + static std::unique_ptr create(Mesh* m = nullptr, Options* opt = nullptr, + CELL_LOC loc = CELL_CENTRE) { + return LaplaceXYFactory::getInstance().create(m, opt, loc); + } - LaplaceXYMonitor monitor; +protected: + static const int INVERT_DC_GRAD = 1; + static const int INVERT_AC_GRAD = 2; // Use zero neumann (NOTE: AC is a misnomer) + static const int INVERT_SET = 16; // Set boundary to x0 value + static const int INVERT_RHS = 32; // Set boundary to b value }; -#endif // BOUT_HAS_PETSC -#endif // BOUT_LAPLACE_XY_H +#endif // BOUT_LAPLACEXY_H diff --git a/src/invert/laplacexy2/laplacexy2_hypre.cxx b/src/invert/laplacexy/impls/hypre/laplacexy-hypre.cxx similarity index 95% rename from src/invert/laplacexy2/laplacexy2_hypre.cxx rename to src/invert/laplacexy/impls/hypre/laplacexy-hypre.cxx index a5028169f8..9323513f21 100644 --- a/src/invert/laplacexy2/laplacexy2_hypre.cxx +++ b/src/invert/laplacexy/impls/hypre/laplacexy-hypre.cxx @@ -1,15 +1,26 @@ -#include +#include "bout/build_defines.hxx" #if BOUT_HAS_HYPRE +#include "laplacexy-hypre.hxx" + #include "bout/assert.hxx" +#include "bout/bout_types.hxx" #include "bout/boutcomm.hxx" +#include "bout/coordinates.hxx" +#include "bout/field2d.hxx" +#include "bout/globalindexer.hxx" #include "bout/globals.hxx" #include "bout/mesh.hxx" +#include "bout/operatorstencil.hxx" +#include "bout/options.hxx" #include "bout/output.hxx" +#include "bout/region.hxx" +#include "bout/sys/range.hxx" #include "bout/sys/timer.hxx" -#include +#include +#include #if BOUT_HAS_CUDA && defined(__CUDACC__) #define gpuErrchk(ans) \ @@ -226,7 +237,7 @@ void LaplaceXY2Hypre::setCoefs(const Field2D& A, const Field2D& B) { } } -Field2D LaplaceXY2Hypre::solve(Field2D& rhs, Field2D& x0) { +Field2D LaplaceXY2Hypre::solve(const Field2D& rhs, const Field2D& x0) { Timer timer("invert"); ASSERT1(rhs.getMesh() == localmesh); diff --git a/include/bout/invert/laplacexy2_hypre.hxx b/src/invert/laplacexy/impls/hypre/laplacexy-hypre.hxx similarity index 71% rename from include/bout/invert/laplacexy2_hypre.hxx rename to src/invert/laplacexy/impls/hypre/laplacexy-hypre.hxx index 1df0988d06..f415f9ebf4 100644 --- a/include/bout/invert/laplacexy2_hypre.hxx +++ b/src/invert/laplacexy/impls/hypre/laplacexy-hypre.hxx @@ -33,51 +33,39 @@ #ifndef LAPLACE_XY2_HYPRE_H #define LAPLACE_XY2_HYPRE_H -#include +#include "bout/build_defines.hxx" +#include "bout/invert/laplacexy.hxx" -#if not BOUT_HAS_HYPRE -// If no Hypre +#if !BOUT_HAS_HYPRE -#include "bout/globalindexer.hxx" -#include -#include -#include - -/*! - * Create a dummy class so that code will compile - * without Hypre, but will throw an exception if - * LaplaceXY is used. - */ -class LaplaceXY2Hypre { -public: - LaplaceXY2Hypre(Mesh* UNUSED(m) = nullptr, Options* UNUSED(opt) = nullptr, - const CELL_LOC UNUSED(loc) = CELL_CENTRE) { - throw BoutException("LaplaceXY requires Hypre. No LaplaceXY available"); - } - void setCoefs(const Field2D& UNUSED(A), const Field2D& UNUSED(B)) {} - Field2D solve(const Field2D& UNUSED(rhs), const Field2D& UNUSED(x0)) { - throw BoutException("LaplaceXY requires Hypre. No LaplaceXY available"); - } -}; +namespace { +RegisterUnavailableLaplaceXY + registerlaplacexyhypre("hypre", "BOUT++ was not configured with HYPRE"); +} #else // BOUT_HAS_HYPRE class Mesh; -#include "bout/utils.hxx" -#include +#include "bout/bout_types.hxx" +#include "bout/field2d.hxx" +#include "bout/globalindexer.hxx" +#include "bout/hypre_interface.hxx" -class LaplaceXY2Hypre { +class LaplaceXY2Hypre : public LaplaceXY { public: - LaplaceXY2Hypre(Mesh* m = nullptr, Options* opt = nullptr, - const CELL_LOC loc = CELL_CENTRE); - ~LaplaceXY2Hypre() = default; + LaplaceXY2Hypre(Mesh* m = nullptr, Options* opt = nullptr, CELL_LOC loc = CELL_CENTRE); + LaplaceXY2Hypre(const LaplaceXY2Hypre&) = delete; + LaplaceXY2Hypre(LaplaceXY2Hypre&&) = delete; + LaplaceXY2Hypre& operator=(const LaplaceXY2Hypre&) = delete; + LaplaceXY2Hypre& operator=(LaplaceXY2Hypre&&) = delete; + ~LaplaceXY2Hypre() override = default; /*! * Set coefficients (A, B) in equation: * Div( A * Grad_perp(x) ) + B*x = b */ - void setCoefs(const Field2D& A, const Field2D& B); + void setCoefs(const Field2D& A, const Field2D& B) override; /*! * Solve Laplacian in X-Y @@ -96,7 +84,7 @@ public: * The solution as a Field2D. On failure an exception will be raised * */ - Field2D solve(Field2D& rhs, Field2D& x0); + Field2D solve(const Field2D& rhs, const Field2D& x0) override; /*! * Preconditioner function @@ -131,5 +119,9 @@ private: MPI_Comm communicator(); }; +namespace { +const inline RegisterLaplaceXY registerlaplacexyhypre{"hypre"}; +} + #endif // BOUT_HAS_HYPRE #endif // LAPLACE_XY_HYPRE_H diff --git a/src/invert/laplacexy/impls/petsc/laplacexy-petsc.cxx b/src/invert/laplacexy/impls/petsc/laplacexy-petsc.cxx new file mode 100644 index 0000000000..a116d78475 --- /dev/null +++ b/src/invert/laplacexy/impls/petsc/laplacexy-petsc.cxx @@ -0,0 +1,1899 @@ +#include "bout/build_defines.hxx" + +#if BOUT_HAS_PETSC + +#include "laplacexy-petsc.hxx" + +#include "bout/assert.hxx" +#include "bout/bout_types.hxx" +#include "bout/boutcomm.hxx" +#include "bout/boutexception.hxx" +#include "bout/cyclic_reduction.hxx" +#include "bout/derivs.hxx" +#include "bout/field2d.hxx" +#include "bout/globals.hxx" +#include "bout/options.hxx" +#include "bout/petsclib.hxx" +#include "bout/solver.hxx" +#include "bout/sys/range.hxx" +#include "bout/sys/timer.hxx" +#include "bout/utils.hxx" + +#include +#include + +#include + +namespace { +PetscErrorCode laplacePCapply(PC pc, Vec x, Vec y) { + // Get the context + LaplaceXYpetsc* s = nullptr; + const auto ierr = PCShellGetContext(pc, reinterpret_cast(&s)); + CHKERRQ(ierr); + + PetscFunctionReturn(s->precon(x, y)); +} +} // namespace + +LaplaceXYpetsc::LaplaceXYpetsc(Mesh* m, Options* opt, const CELL_LOC loc) + : lib(opt == nullptr ? &(Options::root()["laplacexy"]) : opt), + localmesh(m == nullptr ? bout::globals::mesh : m), location(loc), monitor(*this) { + Timer timer("invert"); + + if (opt == nullptr) { + // If no options supplied, use default + opt = &(Options::root()["laplacexy"]); + } + + finite_volume = + (*opt)["finite_volume"] + .doc("Use finite volume rather than finite difference discretisation.") + .withDefault(true); + + /////////////////////////////////////////////////// + // Boundary condititions options + if (localmesh->periodicY(localmesh->xstart)) { + // Periodic in Y, so in the core + x_inner_dirichlet = (*opt)["core_bndry_dirichlet"].withDefault(false); + } else { + // Non-periodic, so in the PF region + x_inner_dirichlet = (*opt)["pf_bndry_dirichlet"].withDefault(true); + } + if ((*opt)["y_bndry_dirichlet"].isSet()) { + throw BoutException("y_bndry_dirichlet has been deprecated. Use y_bndry=dirichlet " + "instead."); + } + y_bndry = (*opt)["y_bndry"].withDefault("neumann"); + + // Check value of y_bndry is a supported option + if (not(y_bndry == "dirichlet" or y_bndry == "neumann" or y_bndry == "free_o3")) { + + throw BoutException("Unrecognized option '{}' for laplacexy:ybndry", y_bndry); + } + + if (not finite_volume) { + // Check we can use corner cells + if (not localmesh->include_corner_cells) { + throw BoutException( + "Finite difference form of LaplaceXYpetsc allows non-orthogonal x- and " + "y-directions, so requires mesh:include_corner_cells=true."); + } + } + + // Use name of options section as the default prefix for performance logging variables + default_prefix = opt->name(); + + // Get MPI communicator + auto* comm = BoutComm::get(); + + // Local size + const int localN = localSize(); + + // Create Vectors + VecCreate(comm, &xs); + VecSetSizes(xs, localN, PETSC_DETERMINE); + VecSetFromOptions(xs); + VecDuplicate(xs, &bs); + + // Set size of Matrix on each processor to localN x localN + MatCreate(comm, &MatA); + MatSetSizes(MatA, localN, localN, PETSC_DETERMINE, PETSC_DETERMINE); + MatSetFromOptions(MatA); + + ////////////////////////////////////////////////// + // Specify local indices. This creates a mapping + // from local indices to index, using a Field2D object + + indexXY = -1; // Set all points to -1, indicating out of domain + + int ind = 0; + + // Y boundaries + for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { + // Should not go into corner cells, LaplaceXYpetsc stencil does not include them + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + indexXY(it.ind, localmesh->ystart - 1) = ind++; + } + if ((not finite_volume) and localmesh->hasBndryLowerY()) { + // Corner boundary cells + if (localmesh->firstX()) { + indexXY(localmesh->xstart - 1, localmesh->ystart - 1) = ind++; + } + if (localmesh->lastX()) { + indexXY(localmesh->xend + 1, localmesh->ystart - 1) = ind++; + } + } + for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { + // Should not go into corner cells, LaplaceXYpetsc stencil does not include them + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + indexXY(it.ind, localmesh->yend + 1) = ind++; + } + if ((not finite_volume) and localmesh->hasBndryUpperY()) { + // Corner boundary cells + if (localmesh->firstX()) { + indexXY(localmesh->xstart - 1, localmesh->yend + 1) = ind++; + } + if (localmesh->lastX()) { + indexXY(localmesh->xend + 1, localmesh->yend + 1) = ind++; + } + } + + xstart = localmesh->xstart; + if (localmesh->firstX()) { + xstart -= 1; // Include X guard cells + } + xend = localmesh->xend; + if (localmesh->lastX()) { + xend += 1; + } + for (int x = xstart; x <= xend; x++) { + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + indexXY(x, y) = ind++; + } + } + + ASSERT1(ind == localN); // Reached end of range + + ////////////////////////////////////////////////// + // Allocate storage for preconditioner + + nloc = xend - xstart + 1; // Number of X points on this processor + nsys = localmesh->yend - localmesh->ystart + 1; // Number of separate Y slices + + acoef.reallocate(nsys, nloc); + bcoef.reallocate(nsys, nloc); + ccoef.reallocate(nsys, nloc); + xvals.reallocate(nsys, nloc); + bvals.reallocate(nsys, nloc); + + // Create a cyclic reduction object + cr = bout::utils::make_unique>(localmesh->getXcomm(), nloc); + + ////////////////////////////////////////////////// + // Pre-allocate PETSc storage + + PetscInt* d_nnz = nullptr; + PetscInt* o_nnz = nullptr; + PetscMalloc((localN) * sizeof(PetscInt), &d_nnz); + PetscMalloc((localN) * sizeof(PetscInt), &o_nnz); + + if (finite_volume) { + setPreallocationFiniteVolume(d_nnz, o_nnz); + } else { + setPreallocationFiniteDifference(d_nnz, o_nnz); + } + // Pre-allocate + MatMPIAIJSetPreallocation(MatA, 0, d_nnz, 0, o_nnz); + MatSetUp(MatA); + + PetscFree(d_nnz); + PetscFree(o_nnz); + + // Determine which row/columns of the matrix are locally owned + int Istart = 0; + int Iend = 0; + MatGetOwnershipRange(MatA, &Istart, &Iend); + + // Convert indexXY from local index to global index + indexXY += Istart; + + // Now communicate to fill guard cells + // Note, this includes corner cells if necessary + localmesh->communicate(indexXY); + + ////////////////////////////////////////////////// + // Set up KSP + + // Declare KSP Context + KSPCreate(comm, &ksp); + + // Configure Linear Solver + + const bool direct = (*opt)["direct"].doc("Use a direct LU solver").withDefault(false); + + if (direct) { + KSPGetPC(ksp, &pc); + PCSetType(pc, PCLU); +#if PETSC_VERSION_GE(3, 9, 0) + PCFactorSetMatSolverType(pc, "mumps"); +#else + PCFactorSetMatSolverPackage(pc, "mumps"); +#endif + } else { + + // Convergence Parameters. Solution is considered converged if |r_k| < max( rtol * |b| , atol ) + // where r_k = b - Ax_k. The solution is considered diverged if |r_k| > dtol * |b|. + + const BoutReal rtol = (*opt)["rtol"].doc("Relative tolerance").withDefault(1e-5); + const BoutReal atol = + (*opt)["atol"] + .doc("Absolute tolerance. The solution is considered converged if |Ax-b| " + "< max( rtol * |b| , atol )") + .withDefault(1e-10); + const BoutReal dtol = + (*opt)["dtol"] + .doc("The solution is considered diverged if |Ax-b| > dtol * |b|") + .withDefault(1e3); + const int maxits = (*opt)["maxits"].doc("Maximum iterations").withDefault(100000); + + // Get KSP Solver Type + const std::string ksptype = + (*opt)["ksptype"].doc("KSP solver type").withDefault("gmres"); + + // Get PC type + const std::string pctype = + (*opt)["pctype"].doc("Preconditioner type").withDefault("none"); + + KSPSetType(ksp, ksptype.c_str()); + KSPSetTolerances(ksp, rtol, atol, dtol, maxits); + + KSPSetInitialGuessNonzero(ksp, static_cast(true)); + + KSPGetPC(ksp, &pc); + PCSetType(pc, pctype.c_str()); + + if (pctype == "shell") { + // Using tridiagonal solver as preconditioner + PCShellSetApply(pc, laplacePCapply); + PCShellSetContext(pc, this); + + const bool rightprec = + (*opt)["rightprec"].doc("Use right preconditioning?").withDefault(true); + if (rightprec) { + KSPSetPCSide(ksp, PC_RIGHT); // Right preconditioning + } else { + KSPSetPCSide(ksp, PC_LEFT); // Left preconditioning + } + } + } + + lib.setOptionsFromInputFile(ksp); + + /////////////////////////////////////////////////// + // Including Y derivatives? + + include_y_derivs = (*opt)["include_y_derivs"] + .doc("Include Y derivatives in operator to invert?") + .withDefault(true); + + /////////////////////////////////////////////////// + // Set the default coefficients + Field2D one(1., localmesh); + Field2D zero(0., localmesh); + one.setLocation(location); + zero.setLocation(location); + setCoefs(one, zero); +} + +void LaplaceXYpetsc::setPreallocationFiniteVolume(PetscInt* d_nnz, PetscInt* o_nnz) { + const int localN = localSize(); + + // This discretisation uses a 5-point stencil + for (int i = 0; i < localN; i++) { + // Non-zero elements on this processor + d_nnz[i] = 5; // Star pattern in 2D + // Non-zero elements on neighboring processor + o_nnz[i] = 0; + } + + // X boundaries + if (localmesh->firstX()) { + // Lower X boundary + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + const int localIndex = globalIndex(localmesh->xstart - 1, y); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + + d_nnz[localIndex] = 2; // Diagonal sub-matrix + o_nnz[localIndex] = 0; // Off-diagonal sub-matrix + } + } else { + // On another processor + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + const int localIndex = globalIndex(localmesh->xstart, y); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + d_nnz[localIndex] -= 1; + o_nnz[localIndex] += 1; + } + } + if (localmesh->lastX()) { + // Upper X boundary + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + const int localIndex = globalIndex(localmesh->xend + 1, y); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + d_nnz[localIndex] = 2; // Diagonal sub-matrix + o_nnz[localIndex] = 0; // Off-diagonal sub-matrix + } + } else { + // On another processor + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + const int localIndex = globalIndex(localmesh->xend, y); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + d_nnz[localIndex] -= 1; + o_nnz[localIndex] += 1; + } + } + // Y boundaries + + for (int x = localmesh->xstart; x <= localmesh->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 + { + const int localIndex = globalIndex(x, localmesh->ystart); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + // d_nnz[localIndex] -= 1; // Note: Slightly inefficient + o_nnz[localIndex] += 1; + } + { + const int localIndex = globalIndex(x, localmesh->yend); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + // d_nnz[localIndex] -= 1; // Note: Slightly inefficient + o_nnz[localIndex] += 1; + } + } + + for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { + // Should not go into corner cells, LaplaceXYpetsc stencil does not include them + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + { + const int localIndex = globalIndex(it.ind, localmesh->ystart - 1); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + d_nnz[localIndex] = 2; // Diagonal sub-matrix + o_nnz[localIndex] = 0; // Off-diagonal sub-matrix + } + { + const int localIndex = globalIndex(it.ind, localmesh->ystart); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + d_nnz[localIndex] += 1; + o_nnz[localIndex] -= 1; + } + } + for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { + // Should not go into corner cells, LaplaceXYpetsc stencil does not include them + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + { + const int localIndex = globalIndex(it.ind, localmesh->yend + 1); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + d_nnz[localIndex] = 2; // Diagonal sub-matrix + o_nnz[localIndex] = 0; // Off-diagonal sub-matrix + } + { + const int localIndex = globalIndex(it.ind, localmesh->yend); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + d_nnz[localIndex] += 1; + o_nnz[localIndex] -= 1; + } + } +} + +void LaplaceXYpetsc::setPreallocationFiniteDifference(PetscInt* d_nnz, PetscInt* o_nnz) { + const int localN = localSize(); + + // This discretisation uses a 9-point stencil + for (int i = 0; i < localN; i++) { + // Non-zero elements on this processor + d_nnz[i] = 9; // Square pattern in 2D + // Non-zero elements on neighboring processor + o_nnz[i] = 0; + } + + // X boundaries + if (localmesh->firstX()) { + // Lower X boundary + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + const int localIndex = globalIndex(localmesh->xstart - 1, y); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + + d_nnz[localIndex] = 2; // Diagonal sub-matrix + o_nnz[localIndex] = 0; // Off-diagonal sub-matrix + } + } else { + // On another processor + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + const int localIndex = globalIndex(localmesh->xstart, y); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + d_nnz[localIndex] -= 3; + o_nnz[localIndex] += 3; + } + } + if (localmesh->lastX()) { + // Upper X boundary + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + const int localIndex = globalIndex(localmesh->xend + 1, y); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + d_nnz[localIndex] = 2; // Diagonal sub-matrix + o_nnz[localIndex] = 0; // Off-diagonal sub-matrix + } + } else { + // On another processor + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + const int localIndex = globalIndex(localmesh->xend, y); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + d_nnz[localIndex] -= 3; + o_nnz[localIndex] += 3; + } + } + // Y boundaries + const int y_bndry_stencil_size = (y_bndry == "free_o3") ? 4 : 2; + + for (int x = localmesh->xstart; x <= localmesh->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 + { + const int localIndex = globalIndex(x, localmesh->ystart); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + // d_nnz[localIndex] -= 3; // Note: Slightly inefficient + o_nnz[localIndex] += 3; + } + { + const int localIndex = globalIndex(x, localmesh->yend); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + // d_nnz[localIndex] -= 3; // Note: Slightly inefficient + o_nnz[localIndex] += 3; + } + } + + for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { + // Should not go into corner cells, they are handled specially below + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + { + const int localIndex = globalIndex(it.ind, localmesh->ystart - 1); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + d_nnz[localIndex] = y_bndry_stencil_size; // Diagonal sub-matrix + o_nnz[localIndex] = 0; // Off-diagonal sub-matrix + } + { + const int localIndex = globalIndex(it.ind, localmesh->ystart); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + //d_nnz[localIndex] += 3; + o_nnz[localIndex] -= 3; + } + } + if (localmesh->hasBndryLowerY()) { + if (y_bndry == "dirichlet") { + // special handling for the corners, since we use a free_o3 y-boundary + // condition just in the corners when y_bndry=="dirichlet" + if (localmesh->firstX()) { + const int localIndex = globalIndex(localmesh->xstart - 1, localmesh->ystart - 1); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + d_nnz[localIndex] = 4; + } + if (localmesh->lastX()) { + const int localIndex = globalIndex(localmesh->xend + 1, localmesh->ystart - 1); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + d_nnz[localIndex] = 4; + } + } else { + if (localmesh->firstX()) { + const int localIndex = globalIndex(localmesh->xstart - 1, localmesh->ystart - 1); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + d_nnz[localIndex] = y_bndry_stencil_size; + } + if (localmesh->lastX()) { + const int localIndex = globalIndex(localmesh->xend + 1, localmesh->ystart - 1); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + d_nnz[localIndex] = y_bndry_stencil_size; + } + } + } + for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { + // Should not go into corner cells, they are handled specially below + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + { + const int localIndex = globalIndex(it.ind, localmesh->yend + 1); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + d_nnz[localIndex] = y_bndry_stencil_size; // Diagonal sub-matrix + o_nnz[localIndex] = 0; // Off-diagonal sub-matrix + } + { + const int localIndex = globalIndex(it.ind, localmesh->yend); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + //d_nnz[localIndex] += 3; + o_nnz[localIndex] -= 3; + } + } + if (localmesh->hasBndryUpperY()) { + if (y_bndry == "dirichlet") { + // special handling for the corners, since we use a free_o3 y-boundary + // condition just in the corners when y_bndry=="dirichlet" + if (localmesh->firstX()) { + const int localIndex = globalIndex(localmesh->xstart - 1, localmesh->yend + 1); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + d_nnz[localIndex] = 4; + } + if (localmesh->lastX()) { + const int localIndex = globalIndex(localmesh->xend + 1, localmesh->yend + 1); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + d_nnz[localIndex] = 4; + } + } else { + if (localmesh->firstX()) { + const int localIndex = globalIndex(localmesh->xstart - 1, localmesh->yend + 1); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + d_nnz[localIndex] = y_bndry_stencil_size; + } + if (localmesh->lastX()) { + const int localIndex = globalIndex(localmesh->xend + 1, localmesh->yend + 1); + ASSERT1((localIndex >= 0) && (localIndex < localN)); + d_nnz[localIndex] = y_bndry_stencil_size; + } + } + } +} + +void LaplaceXYpetsc::setCoefs(const Field2D& A, const Field2D& B) { + Timer timer("invert"); + + ASSERT1(A.getMesh() == localmesh); + ASSERT1(B.getMesh() == localmesh); + ASSERT1(A.getLocation() == location); + ASSERT1(B.getLocation() == location); + + if (finite_volume) { + setMatrixElementsFiniteVolume(A, B); + } else { + setMatrixElementsFiniteDifference(A, B); + } + + // X boundaries + if (localmesh->firstX()) { + if (x_inner_dirichlet) { + + // Dirichlet on inner X boundary + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + int row = globalIndex(localmesh->xstart - 1, y); + PetscScalar val = 0.5; + MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); + + int col = globalIndex(localmesh->xstart, y); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + + // Preconditioner + bcoef(y - localmesh->ystart, 0) = 0.5; + ccoef(y - localmesh->ystart, 0) = 0.5; + } + + } else { + + // Neumann on inner X boundary + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + int row = globalIndex(localmesh->xstart - 1, y); + PetscScalar val = 1.0; + MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); + + int col = globalIndex(localmesh->xstart, y); + val = -1.0; + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + + // Preconditioner + bcoef(y - localmesh->ystart, 0) = 1.0; + ccoef(y - localmesh->ystart, 0) = -1.0; + } + } + } + if (localmesh->lastX()) { + // Dirichlet on outer X boundary + + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + int row = globalIndex(localmesh->xend + 1, y); + PetscScalar val = 0.5; + MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); + + int col = globalIndex(localmesh->xend, y); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + + // Preconditioner + acoef(y - localmesh->ystart, localmesh->xend + 1 - xstart) = 0.5; + bcoef(y - localmesh->ystart, localmesh->xend + 1 - xstart) = 0.5; + } + } + + if (y_bndry == "dirichlet") { + // Dirichlet on Y boundaries + for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { + // Should not go into corner cells, LaplaceXYpetsc stencil does not include them + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + int row = globalIndex(it.ind, localmesh->ystart - 1); + PetscScalar val = 0.5; + MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); + + int col = globalIndex(it.ind, localmesh->ystart); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + } + + for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { + // Should not go into corner cells, LaplaceXYpetsc stencil does not include them + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + int row = globalIndex(it.ind, localmesh->yend + 1); + PetscScalar val = 0.5; + MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); + + int col = globalIndex(it.ind, localmesh->yend); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + } + } else if (y_bndry == "neumann") { + // Neumann on Y boundaries + for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { + // Should not go into corner cells, LaplaceXYpetsc stencil does not include them + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + int row = globalIndex(it.ind, localmesh->ystart - 1); + PetscScalar val = 1.0; + MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); + + val = -1.0; + int col = globalIndex(it.ind, localmesh->ystart); + + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + } + + for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { + // Should not go into corner cells, LaplaceXYpetsc stencil does not include them + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + int row = globalIndex(it.ind, localmesh->yend + 1); + PetscScalar val = 1.0; + MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); + + val = -1.0; + int col = globalIndex(it.ind, localmesh->yend); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + } + } else if (y_bndry == "free_o3") { + // 'free_o3' extrapolating boundary condition on Y boundaries + for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { + // Should not go into corner cells, LaplaceXYpetsc stencil does not include them + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + int row = globalIndex(it.ind, localmesh->ystart - 1); + PetscScalar val = 1.0; + MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); + + val = -3.0; + int col = globalIndex(it.ind, localmesh->ystart); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + + val = 3.0; + col = globalIndex(it.ind, localmesh->ystart + 1); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + + val = -1.0; + col = globalIndex(it.ind, localmesh->ystart + 2); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + } + + for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { + // Should not go into corner cells, LaplaceXYpetsc stencil does not include them + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + int row = globalIndex(it.ind, localmesh->yend + 1); + PetscScalar val = 1.0; + MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); + + val = -3.0; + int col = globalIndex(it.ind, localmesh->yend); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + + val = 3.0; + col = globalIndex(it.ind, localmesh->yend - 1); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + + val = -1.0; + col = globalIndex(it.ind, localmesh->yend - 2); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + } + } else { + throw BoutException("Unsupported option for y_bndry"); + } + + if (not finite_volume) { + // Handle corner boundary cells in case we need to include D2DXDY + // Apply the y-boundary-condition to the cells in the x-boundary - this is an + // arbitrary choice, cf. connections around the X-point + + if (localmesh->hasBndryLowerY()) { + if (localmesh->firstX()) { + if (y_bndry == "neumann") { + // Neumann y-bc + // f(xs-1,ys-1) = f(xs-1,ys) + PetscScalar val = 1.0; + int row = globalIndex(localmesh->xstart - 1, localmesh->ystart - 1); + MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); + + val = -1.0; + int col = globalIndex(localmesh->xstart - 1, localmesh->ystart); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + } else if (y_bndry == "free_o3" or y_bndry == "dirichlet") { + // 'free_o3' extrapolating boundary condition on Y boundaries + // f(xs-1,ys-1) = 3*f(xs-1,ys) - 3*f(xs-1,ys+1) + f(xs-1,ys+2) + // + // Use free_o3 at the corners for Dirichlet y-boundaries because we don't know + // what value to pass for the corner + PetscScalar val = 1.0; + int row = globalIndex(localmesh->xstart - 1, localmesh->ystart - 1); + MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); + + val = -3.0; + int col = globalIndex(localmesh->xstart - 1, localmesh->ystart); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + + val = 3.0; + col = globalIndex(localmesh->xstart - 1, localmesh->ystart + 1); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + + val = -1.0; + col = globalIndex(localmesh->xstart - 1, localmesh->ystart + 2); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + } else { + throw BoutException("Unsupported option for y_bndry"); + } + } + if (localmesh->lastX()) { + if (y_bndry == "neumann") { + // Neumann y-bc + // f(xe+1,ys-1) = f(xe+1,ys) + PetscScalar val = 1.0; + int row = globalIndex(localmesh->xend + 1, localmesh->ystart - 1); + MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); + + val = -1.0; + int col = globalIndex(localmesh->xend + 1, localmesh->ystart); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + } else if (y_bndry == "free_o3" or y_bndry == "dirichlet") { + // 'free_o3' extrapolating boundary condition on Y boundaries + // f(xe+1,ys-1) = 3*f(xe+1,ys) - 3*f(xe+1,ys+1) + f(xe+1,ys+2) + // + // Use free_o3 at the corners for Dirichlet y-boundaries because we don't know + // what value to pass for the corner + PetscScalar val = 1.0; + int row = globalIndex(localmesh->xend + 1, localmesh->ystart - 1); + MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); + + val = -3.0; + int col = globalIndex(localmesh->xend + 1, localmesh->ystart); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + + val = 3.0; + col = globalIndex(localmesh->xend + 1, localmesh->ystart + 1); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + + val = -1.0; + col = globalIndex(localmesh->xend + 1, localmesh->ystart + 2); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + } else { + throw BoutException("Unsupported option for y_bndry"); + } + } + } + if (localmesh->hasBndryUpperY()) { + if (localmesh->firstX()) { + if (y_bndry == "neumann") { + // Neumann y-bc + // f(xs-1,ys-1) = f(xs-1,ys) + PetscScalar val = 1.0; + int row = globalIndex(localmesh->xstart - 1, localmesh->yend + 1); + MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); + + val = -1.0; + int col = globalIndex(localmesh->xstart - 1, localmesh->yend); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + } else if (y_bndry == "free_o3" or y_bndry == "dirichlet") { + // 'free_o3' extrapolating boundary condition on Y boundaries + // f(xs-1,ys-1) = 3*f(xs-1,ys) - 3*f(xs-1,ys+1) + f(xs-1,ys+2) + // + // Use free_o3 at the corners for Dirichlet y-boundaries because we don't know + // what value to pass for the corner + PetscScalar val = 1.0; + int row = globalIndex(localmesh->xstart - 1, localmesh->yend + 1); + MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); + + val = -3.0; + int col = globalIndex(localmesh->xstart - 1, localmesh->yend); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + + val = 3.0; + col = globalIndex(localmesh->xstart - 1, localmesh->yend - 1); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + + val = -1.0; + col = globalIndex(localmesh->xstart - 1, localmesh->yend - 2); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + } else { + throw BoutException("Unsupported option for y_bndry"); + } + } + if (localmesh->lastX()) { + if (y_bndry == "neumann") { + // Neumann y-bc + // f(xe+1,ys-1) = f(xe+1,ys) + PetscScalar val = 1.0; + int row = globalIndex(localmesh->xend + 1, localmesh->yend + 1); + MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); + + val = -1.0; + int col = globalIndex(localmesh->xend + 1, localmesh->yend); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + } else if (y_bndry == "free_o3" or y_bndry == "dirichlet") { + // 'free_o3' extrapolating boundary condition on Y boundaries + // f(xe+1,ys-1) = 3*f(xe+1,ys) - 3*f(xe+1,ys+1) + f(xe+1,ys+2) + // + // Use free_o3 at the corners for Dirichlet y-boundaries because we don't know + // what value to pass for the corner + PetscScalar val = 1.0; + int row = globalIndex(localmesh->xend + 1, localmesh->yend + 1); + MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); + + val = -3.0; + int col = globalIndex(localmesh->xend + 1, localmesh->yend); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + + val = 3.0; + col = globalIndex(localmesh->xend + 1, localmesh->yend - 1); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + + val = -1.0; + col = globalIndex(localmesh->xend + 1, localmesh->yend - 2); + MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); + } else { + throw BoutException("Unsupported option for y_bndry"); + } + } + } + } + + // Assemble Matrix + MatAssemblyBegin(MatA, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(MatA, MAT_FINAL_ASSEMBLY); + + // Set the operator +#if PETSC_VERSION_GE(3, 5, 0) + KSPSetOperators(ksp, MatA, MatA); +#else + KSPSetOperators(ksp, MatA, MatA, DIFFERENT_NONZERO_PATTERN); +#endif + + // Set coefficients for preconditioner + cr->setCoefs(acoef, bcoef, ccoef); +} + +void LaplaceXYpetsc::setMatrixElementsFiniteVolume(const Field2D& A, const Field2D& B) { + ////////////////////////////////////////////////// + // Set Matrix elements + // + // (1/J) d/dx ( J * g11 d/dx ) + (1/J) d/dy ( J * g22 d/dy ) + + auto coords = localmesh->getCoordinates(location); + const Field2D J_DC = DC(coords->J); + const Field2D g11_DC = DC(coords->g11); + const Field2D dx_DC = DC(coords->dx); + const Field2D dy_DC = DC(coords->dy); + const Field2D g_22_DC = DC(coords->g_22); + const Field2D g_23_DC = DC(coords->g_23); + const Field2D g23_DC = DC(coords->g23); + + for (int x = localmesh->xstart; x <= localmesh->xend; x++) { + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + // stencil entries + PetscScalar c, xm, xp, ym, yp; + + // XX component + + // Metrics on x+1/2 boundary + BoutReal J = 0.5 * (J_DC(x, y) + J_DC(x + 1, y)); + BoutReal g11 = 0.5 * (g11_DC(x, y) + g11_DC(x + 1, y)); + BoutReal dx = 0.5 * (dx_DC(x, y) + dx_DC(x + 1, y)); + BoutReal Acoef = 0.5 * (A(x, y) + A(x + 1, y)); + + BoutReal val = Acoef * J * g11 / (J_DC(x, y) * dx * dx_DC(x, y)); + xp = val; + c = -val; + + // Metrics on x-1/2 boundary + J = 0.5 * (J_DC(x, y) + J_DC(x - 1, y)); + g11 = 0.5 * (g11_DC(x, y) + g11_DC(x - 1, y)); + dx = 0.5 * (dx_DC(x, y) + dx_DC(x - 1, y)); + Acoef = 0.5 * (A(x, y) + A(x - 1, y)); + + val = Acoef * J * g11 / (J_DC(x, y) * dx * dx_DC(x, y)); + xm = val; + c -= val; + + c += B(x, y); + + // Put values into the preconditioner, X derivatives only + acoef(y - localmesh->ystart, x - xstart) = xm; + bcoef(y - localmesh->ystart, x - xstart) = c; + ccoef(y - localmesh->ystart, x - xstart) = xp; + + if (include_y_derivs) { + // YY component + // Metrics at y+1/2 + J = 0.5 * (J_DC(x, y) + J_DC(x, y + 1)); + BoutReal g_22 = 0.5 * (g_22_DC(x, y) + g_22_DC(x, y + 1)); + BoutReal g23 = 0.5 * (g23_DC(x, y) + g23_DC(x, y + 1)); + BoutReal g_23 = 0.5 * (g_23_DC(x, y) + g_23_DC(x, y + 1)); + BoutReal dy = 0.5 * (dy_DC(x, y) + dy_DC(x, y + 1)); + Acoef = 0.5 * (A(x, y + 1) + A(x, y)); + + val = -Acoef * J * g23 * g_23 / (g_22 * J_DC(x, y) * dy * dy_DC(x, y)); + yp = val; + c -= val; + + // Metrics at y-1/2 + J = 0.5 * (J_DC(x, y) + J_DC(x, y - 1)); + g_22 = 0.5 * (g_22_DC(x, y) + g_22_DC(x, y - 1)); + g23 = 0.5 * (g23_DC(x, y) + g23_DC(x, y - 1)); + g_23 = 0.5 * (g_23_DC(x, y) + g_23_DC(x, y - 1)); + dy = 0.5 * (dy_DC(x, y) + dy_DC(x, y - 1)); + Acoef = 0.5 * (A(x, y - 1) + A(x, y)); + + val = -Acoef * J * g23 * g_23 / (g_22 * J_DC(x, y) * dy * dy_DC(x, y)); + ym = val; + c -= val; + } + + ///////////////////////////////////////////////// + // Now have a 5-point stencil for the Laplacian + + int row = globalIndex(x, y); + + // Set the centre (diagonal) + MatSetValues(MatA, 1, &row, 1, &row, &c, INSERT_VALUES); + + // X + 1 + int col = globalIndex(x + 1, y); + MatSetValues(MatA, 1, &row, 1, &col, &xp, INSERT_VALUES); + + // X - 1 + col = globalIndex(x - 1, y); + MatSetValues(MatA, 1, &row, 1, &col, &xm, INSERT_VALUES); + + if (include_y_derivs) { + // Y + 1 + col = globalIndex(x, y + 1); + MatSetValues(MatA, 1, &row, 1, &col, &yp, INSERT_VALUES); + + // Y - 1 + col = globalIndex(x, y - 1); + MatSetValues(MatA, 1, &row, 1, &col, &ym, INSERT_VALUES); + } + } + } +} + +void LaplaceXYpetsc::setMatrixElementsFiniteDifference(const Field2D& A, + const Field2D& B) { + ////////////////////////////////////////////////// + // Set Matrix elements + // + // Div(A Grad(f)) + B f + // = A Laplace_perp(f) + Grad_perp(A).Grad_perp(f) + B f + // = A*(G1*dfdx + (G2-1/J*d/dy(J/g_22))*dfdy + // + g11*d2fdx2 + (g22-1/g_22)*d2fdy2 + 2*g12*d2fdxdy) + // + g11*dAdx*dfdx + (g22-1/g_22)*dAdy*dfdy + g12*(dAdx*dfdy + dAdy*dfdx) + // + B*f + + auto coords = localmesh->getCoordinates(location); + const Field2D G1_2D = DC(coords->G1); + const Field2D G2_2D = DC(coords->G2); + const Field2D J_2D = DC(coords->J); + const Field2D g11_2D = DC(coords->g11); + const Field2D g_22_2D = DC(coords->g_22); + const Field2D g22_2D = DC(coords->g22); + const Field2D g12_2D = DC(coords->g12); + const Field2D d1_dx_2D = DC(coords->d1_dx); + const Field2D d1_dy_2D = DC(coords->d1_dy); + const Field2D dx_2D = DC(coords->dx); + const Field2D dy_2D = DC(coords->dy); + + const Field2D coef_dfdy = G2_2D - DC(DDY(J_2D / g_22_2D) / J_2D); + + for (int x = localmesh->xstart; x <= localmesh->xend; x++) { + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + // stencil entries + PetscScalar c, xm, xp, ym, yp, xpyp, xpym, xmyp, xmym; + + BoutReal dx = dx_2D(x, y); + + // A*G1*dfdx + BoutReal val = A(x, y) * G1_2D(x, y) / (2. * dx); + xp = val; + xm = -val; + + // A*g11*d2fdx2 + val = A(x, y) * g11_2D(x, y) / SQ(dx); + xp += val; + c = -2. * val; + xm += val; + // Non-uniform grid correction + val = A(x, y) * g11_2D(x, y) * d1_dx_2D(x, y) / (2. * dx); + xp += val; + xm -= val; + + // g11*dAdx*dfdx + val = g11_2D(x, y) * (A(x + 1, y) - A(x - 1, y)) / (4. * SQ(dx)); + xp += val; + xm -= val; + + // B*f + c += B(x, y); + + // Put values into the preconditioner, X derivatives only + acoef(y - localmesh->ystart, x - xstart) = xm; + bcoef(y - localmesh->ystart, x - xstart) = c; + ccoef(y - localmesh->ystart, x - xstart) = xp; + + if (include_y_derivs) { + BoutReal dy = dy_2D(x, y); + BoutReal dAdx = (A(x + 1, y) - A(x - 1, y)) / (2. * dx); + BoutReal dAdy = (A(x, y + 1) - A(x, y - 1)) / (2. * dy); + + // A*(G2-1/J*d/dy(J/g_22))*dfdy + val = A(x, y) * coef_dfdy(x, y) / (2. * dy); + yp = val; + ym = -val; + + // A*(g22-1/g_22)*d2fdy2 + val = A(x, y) * (g22_2D(x, y) - 1. / g_22_2D(x, y)) / SQ(dy); + yp += val; + c -= 2. * val; + ym += val; + // Non-uniform mesh correction + val = A(x, y) * (g22_2D(x, y) - 1. / g_22_2D(x, y)) * d1_dy_2D(x, y) / (2. * dy); + yp += val; + ym -= val; + + // 2*A*g12*d2dfdxdy + val = A(x, y) * g12_2D(x, y) / (2. * dx * dy); + xpyp = val; + xpym = -val; + xmyp = -val; + xmym = val; + + // g22*dAdy*dfdy + val = (g22_2D(x, y) - 1. / g_22_2D(x, y)) * dAdy / (2. * dy); + yp += val; + ym -= val; + + // g12*(dAdx*dfdy + dAdy*dfdx) + val = g12_2D(x, y) * dAdx / (2. * dy); + yp += val; + ym -= val; + val = g12_2D(x, y) * dAdy / (2. * dx); + xp += val; + xm -= val; + } + + ///////////////////////////////////////////////// + // Now have a 9-point stencil for the Laplacian + + int row = globalIndex(x, y); + + // Set the centre (diagonal) + MatSetValues(MatA, 1, &row, 1, &row, &c, INSERT_VALUES); + + // X + 1 + int col = globalIndex(x + 1, y); + MatSetValues(MatA, 1, &row, 1, &col, &xp, INSERT_VALUES); + + // X - 1 + col = globalIndex(x - 1, y); + MatSetValues(MatA, 1, &row, 1, &col, &xm, INSERT_VALUES); + + if (include_y_derivs) { + // Y + 1 + col = globalIndex(x, y + 1); + MatSetValues(MatA, 1, &row, 1, &col, &yp, INSERT_VALUES); + + // Y - 1 + col = globalIndex(x, y - 1); + MatSetValues(MatA, 1, &row, 1, &col, &ym, INSERT_VALUES); + + // X + 1, Y + 1 + col = globalIndex(x + 1, y + 1); + MatSetValues(MatA, 1, &row, 1, &col, &xpyp, INSERT_VALUES); + + // X + 1, Y - 1 + col = globalIndex(x + 1, y - 1); + MatSetValues(MatA, 1, &row, 1, &col, &xpym, INSERT_VALUES); + + // X - 1, Y + 1 + col = globalIndex(x - 1, y + 1); + MatSetValues(MatA, 1, &row, 1, &col, &xmyp, INSERT_VALUES); + + // X - 1, Y - 1 + col = globalIndex(x - 1, y - 1); + MatSetValues(MatA, 1, &row, 1, &col, &xmym, INSERT_VALUES); + } + } + } +} + +LaplaceXYpetsc::~LaplaceXYpetsc() { + PetscBool is_finalised; + PetscFinalized(&is_finalised); + + if (!is_finalised) { + // PetscFinalize may already have destroyed this object + KSPDestroy(&ksp); + } + + VecDestroy(&xs); + VecDestroy(&bs); + MatDestroy(&MatA); +} + +Field2D LaplaceXYpetsc::solve(const Field2D& rhs, const Field2D& x0) { + Timer timer("invert"); + + ASSERT1(rhs.getMesh() == localmesh); + ASSERT1(x0.getMesh() == localmesh); + ASSERT1(rhs.getLocation() == location); + ASSERT1(x0.getLocation() == location); + + // Load initial guess x0 into xs and rhs into bs + + for (int x = localmesh->xstart; x <= localmesh->xend; x++) { + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + int ind = globalIndex(x, y); + + PetscScalar val = x0(x, y); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + val = rhs(x, y); + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + } + + if (finite_volume) { + solveFiniteVolume(x0); + } else { + solveFiniteDifference(x0); + } + + // Assemble RHS Vector + VecAssemblyBegin(bs); + VecAssemblyEnd(bs); + + // Assemble Trial Solution Vector + VecAssemblyBegin(xs); + VecAssemblyEnd(xs); + + // Solve the system + KSPSolve(ksp, bs, xs); + + KSPConvergedReason reason; + KSPGetConvergedReason(ksp, &reason); + + if (reason <= 0) { + throw BoutException("LaplaceXYpetsc failed to converge. Reason {} ({:d})", + KSPConvergedReasons[reason], static_cast(reason)); + } + + if (save_performance) { + // Update performance monitoring information + n_calls++; + + int iterations = 0; + KSPGetIterationNumber(ksp, &iterations); + + average_iterations = BoutReal(n_calls - 1) / BoutReal(n_calls) * average_iterations + + BoutReal(iterations) / BoutReal(n_calls); + } + + ////////////////////////// + // Copy data into result + + Field2D result; + result.allocate(); + result.setLocation(location); + + for (int x = localmesh->xstart; x <= localmesh->xend; x++) { + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + int ind = globalIndex(x, y); + + PetscScalar val; + VecGetValues(xs, 1, &ind, &val); + result(x, y) = val; + } + } + + // Inner X boundary + if (localmesh->firstX()) { + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + int ind = globalIndex(localmesh->xstart - 1, y); + PetscScalar val; + VecGetValues(xs, 1, &ind, &val); + for (int x = localmesh->xstart - 1; x >= 0; x--) { + result(x, y) = val; + } + } + } + + // Outer X boundary + if (localmesh->lastX()) { + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + int ind = globalIndex(localmesh->xend + 1, y); + PetscScalar val; + VecGetValues(xs, 1, &ind, &val); + for (int x = localmesh->xend + 1; x < localmesh->LocalNx; x++) { + result(x, y) = val; + } + } + } + + // Lower Y boundary + for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { + if ( + // Should not go into corner cells, finite-volume LaplaceXYpetsc stencil does not include + // them + (finite_volume and (it.ind < localmesh->xstart or it.ind > localmesh->xend)) + + // Only go into first corner cell for finite-difference + or (it.ind < localmesh->xstart - 1 or it.ind > localmesh->xend + 1)) { + continue; + } + int ind = globalIndex(it.ind, localmesh->ystart - 1); + PetscScalar val; + VecGetValues(xs, 1, &ind, &val); + for (int y = localmesh->ystart - 1; y >= 0; y--) { + result(it.ind, y) = val; + } + } + + // Upper Y boundary + for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { + if ( + // Should not go into corner cells, finite-volume LaplaceXYpetsc stencil does not include + // them + (finite_volume and (it.ind < localmesh->xstart or it.ind > localmesh->xend)) + + // Only go into first corner cell for finite-difference + or (it.ind < localmesh->xstart - 1 or it.ind > localmesh->xend + 1)) { + continue; + } + int ind = globalIndex(it.ind, localmesh->yend + 1); + PetscScalar val; + VecGetValues(xs, 1, &ind, &val); + for (int y = localmesh->yend + 1; y < localmesh->LocalNy; y++) { + result(it.ind, y) = val; + } + } + + return result; +} + +void LaplaceXYpetsc::solveFiniteVolume(const Field2D& x0) { + // Use original LaplaceXYpetsc implementation of passing boundary values for backward + // compatibility + if (localmesh->firstX()) { + if (x_inner_dirichlet) { + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + int ind = globalIndex(localmesh->xstart - 1, y); + + PetscScalar val = x0(localmesh->xstart - 1, y); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + val = 0.5 * (x0(localmesh->xstart - 1, y) + x0(localmesh->xstart, y)); + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + } else { + // Inner X boundary (Neumann) + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + int ind = globalIndex(localmesh->xstart - 1, y); + + PetscScalar val = x0(localmesh->xstart - 1, y); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + val = 0.0; // x0(localmesh->xstart-1,y) - x0(localmesh->xstart,y); + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + } + } + + // Outer X boundary (Dirichlet) + if (localmesh->lastX()) { + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + int ind = globalIndex(localmesh->xend + 1, y); + + PetscScalar val = x0(localmesh->xend + 1, y); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + val = 0.5 * (x0(localmesh->xend, y) + x0(localmesh->xend + 1, y)); + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + } + + if (y_bndry == "dirichlet") { + for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { + // Should not go into corner cells, LaplaceXYpetsc stencil does not include them + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + int ind = globalIndex(it.ind, localmesh->ystart - 1); + + PetscScalar val = x0(it.ind, localmesh->ystart - 1); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + val = 0.5 * (x0(it.ind, localmesh->ystart - 1) + x0(it.ind, localmesh->ystart)); + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + + for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { + // Should not go into corner cells, LaplaceXYpetsc stencil does not include them + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + int ind = globalIndex(it.ind, localmesh->yend + 1); + + PetscScalar val = x0(it.ind, localmesh->yend + 1); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + val = 0.5 * (x0(it.ind, localmesh->yend + 1) + x0(it.ind, localmesh->yend)); + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + } else if (y_bndry == "neumann" or y_bndry == "free_o3") { + // Y boundaries Neumann + for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { + // Should not go into corner cells, LaplaceXYpetsc stencil does not include them + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + int ind = globalIndex(it.ind, localmesh->ystart - 1); + + PetscScalar val = x0(it.ind, localmesh->ystart - 1); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + val = 0.0; + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + + for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { + // Should not go into corner cells, LaplaceXYpetsc stencil does not include them + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + int ind = globalIndex(it.ind, localmesh->yend + 1); + + PetscScalar val = x0(it.ind, localmesh->yend + 1); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + val = 0.0; + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + } else { + throw BoutException("Unsupported option for y_bndry"); + } +} + +void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) { + // For finite-difference implementation pass boundary values in the same way as for + // the 'Laplacian' solvers - the value to use (for Dirichlet boundary conditions) on + // the boundary (which is half way between grid cell and boundary cell) is passed as + // the value in the first boundary cell. + if (localmesh->firstX()) { + if (x_inner_dirichlet) { + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + int ind = globalIndex(localmesh->xstart - 1, y); + + // For the boundary value of the initial guess, use the value that would be set by + // applying the boundary condition to the initial guess + PetscScalar val = 2. * x0(localmesh->xstart - 1, y) - x0(localmesh->xstart, y); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + // Pass the value from boundary cell of x0 as the boundary condition to the rhs + val = x0(localmesh->xstart - 1, y); + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + } else { + // Inner X boundary (Neumann) + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + int ind = globalIndex(localmesh->xstart - 1, y); + + // For the boundary value of the initial guess, use the value that would be set by + // applying the boundary condition to the initial guess + PetscScalar val = x0(localmesh->xstart, y); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + val = 0.0; + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + } + } + + // Outer X boundary (Dirichlet) + if (localmesh->lastX()) { + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + int ind = globalIndex(localmesh->xend + 1, y); + + // For the boundary value of the initial guess, use the value that would be set by + // applying the boundary condition to the initial guess + PetscScalar val = 2. * x0(localmesh->xend + 1, y) - x0(localmesh->xend, y); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + // Pass the value from boundary cell of x0 as the boundary condition to the rhs + val = x0(localmesh->xend + 1, y); + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + } + + if (y_bndry == "dirichlet") { + for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { + // Should not go into corner cells, they are treated specially below + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + int ind = globalIndex(it.ind, localmesh->ystart - 1); + + // For the boundary value of the initial guess, use the value that would be set by + // applying the boundary condition to the initial guess + PetscScalar val = + 2. * x0(it.ind, localmesh->ystart - 1) - x0(it.ind, localmesh->ystart); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + // Pass the value from boundary cell of x0 as the boundary condition to the rhs + val = x0(it.ind, localmesh->ystart - 1); + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + + for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { + // Should not go into corner cells, they are treated specially below + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + int ind = globalIndex(it.ind, localmesh->yend + 1); + + // For the boundary value of the initial guess, use the value that would be set by + // applying the boundary condition to the initial guess + PetscScalar val = + 2. * x0(it.ind, localmesh->yend + 1) - x0(it.ind, localmesh->yend); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + // Pass the value from boundary cell of x0 as the boundary condition to the rhs + val = x0(it.ind, localmesh->yend + 1); + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + + // Use free_o3 for the corner boundary cells + if (localmesh->hasBndryLowerY()) { + if (localmesh->firstX()) { + int ind = globalIndex(localmesh->xstart - 1, localmesh->ystart - 1); + + // For the boundary value of the initial guess, use the value that would be set by + // applying the boundary condition to the initial guess + PetscScalar val = 3. * x0(localmesh->xstart - 1, localmesh->ystart) + - 3. * x0(localmesh->xstart - 1, localmesh->ystart + 1) + + x0(localmesh->xstart - 1, localmesh->ystart + 2); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + // Set the value for the rhs of the boundary condition to zero + val = 0.0; + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + if (localmesh->lastX()) { + int ind = globalIndex(localmesh->xend + 1, localmesh->ystart - 1); + + // For the boundary value of the initial guess, use the value that would be set by + // applying the boundary condition to the initial guess + PetscScalar val = 3. * x0(localmesh->xend + 1, localmesh->ystart) + - 3. * x0(localmesh->xend + 1, localmesh->ystart + 1) + + x0(localmesh->xend + 1, localmesh->ystart + 2); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + // Set the value for the rhs of the boundary condition to zero + val = 0.0; + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + } + if (localmesh->hasBndryUpperY()) { + if (localmesh->firstX()) { + int ind = globalIndex(localmesh->xstart - 1, localmesh->yend + 1); + + // For the boundary value of the initial guess, use the value that would be set by + // applying the boundary condition to the initial guess + PetscScalar val = 3. * x0(localmesh->xstart - 1, localmesh->yend) + - 3. * x0(localmesh->xstart - 1, localmesh->yend - 1) + + x0(localmesh->xstart - 1, localmesh->yend - 2); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + // Set the value for the rhs of the boundary condition to zero + val = 0.0; + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + if (localmesh->lastX()) { + int ind = globalIndex(localmesh->xend + 1, localmesh->yend + 1); + + // For the boundary value of the initial guess, use the value that would be set by + // applying the boundary condition to the initial guess + PetscScalar val = 3. * x0(localmesh->xend + 1, localmesh->yend) + - 3. * x0(localmesh->xend + 1, localmesh->yend - 1) + + x0(localmesh->xend + 1, localmesh->yend - 2); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + // Set the value for the rhs of the boundary condition to zero + val = 0.0; + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + } + } else if (y_bndry == "neumann") { + // Y boundaries Neumann + for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { + // Should not go into corner cells, they are treated specially below + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + int ind = globalIndex(it.ind, localmesh->ystart - 1); + + // For the boundary value of the initial guess, use the value that would be set by + // applying the boundary condition to the initial guess + PetscScalar val = x0(it.ind, localmesh->ystart); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + // Set the value for the rhs of the boundary condition to zero + val = 0.0; + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + if (localmesh->hasBndryLowerY()) { + if (localmesh->firstX()) { + int ind = globalIndex(localmesh->xstart - 1, localmesh->ystart - 1); + + // For the boundary value of the initial guess, use the value that would be set by + // applying the boundary condition to the initial guess + PetscScalar val = x0(localmesh->xstart - 1, localmesh->ystart); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + // Set the value for the rhs of the boundary condition to zero + val = 0.0; + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + if (localmesh->lastX()) { + int ind = globalIndex(localmesh->xend + 1, localmesh->ystart - 1); + + // For the boundary value of the initial guess, use the value that would be set by + // applying the boundary condition to the initial guess + PetscScalar val = x0(localmesh->xend + 1, localmesh->ystart); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + // Set the value for the rhs of the boundary condition to zero + val = 0.0; + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + } + + for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { + // Should not go into corner cells, they are treated specially below + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + int ind = globalIndex(it.ind, localmesh->yend + 1); + + // For the boundary value of the initial guess, use the value that would be set by + // applying the boundary condition to the initial guess + PetscScalar val = x0(it.ind, localmesh->yend); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + // Set the value for the rhs of the boundary condition to zero + val = 0.0; + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + if (localmesh->hasBndryUpperY()) { + if (localmesh->firstX()) { + int ind = globalIndex(localmesh->xstart - 1, localmesh->yend + 1); + + // For the boundary value of the initial guess, use the value that would be set by + // applying the boundary condition to the initial guess + PetscScalar val = x0(localmesh->xstart - 1, localmesh->yend); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + // Set the value for the rhs of the boundary condition to zero + val = 0.0; + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + if (localmesh->lastX()) { + int ind = globalIndex(localmesh->xend + 1, localmesh->yend + 1); + + // For the boundary value of the initial guess, use the value that would be set by + // applying the boundary condition to the initial guess + PetscScalar val = x0(localmesh->xend + 1, localmesh->yend); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + // Set the value for the rhs of the boundary condition to zero + val = 0.0; + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + } + } else if (y_bndry == "free_o3") { + // Y boundaries free_o3 + for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { + // Should not go into corner cells, they are treated specially below + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + int ind = globalIndex(it.ind, localmesh->ystart - 1); + + // For the boundary value of the initial guess, use the value that would be set by + // applying the boundary condition to the initial guess + PetscScalar val = 3. * x0(it.ind, localmesh->ystart) + - 3. * x0(it.ind, localmesh->ystart + 1) + + x0(it.ind, localmesh->ystart + 2); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + // Set the value for the rhs of the boundary condition to zero + val = 0.0; + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + if (localmesh->hasBndryLowerY()) { + if (localmesh->firstX()) { + int ind = globalIndex(localmesh->xstart - 1, localmesh->ystart - 1); + + // For the boundary value of the initial guess, use the value that would be set by + // applying the boundary condition to the initial guess + PetscScalar val = 3. * x0(localmesh->xstart - 1, localmesh->ystart) + - 3. * x0(localmesh->xstart - 1, localmesh->ystart + 1) + + x0(localmesh->xstart - 1, localmesh->ystart + 2); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + // Set the value for the rhs of the boundary condition to zero + val = 0.0; + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + if (localmesh->lastX()) { + int ind = globalIndex(localmesh->xend + 1, localmesh->ystart - 1); + + // For the boundary value of the initial guess, use the value that would be set by + // applying the boundary condition to the initial guess + PetscScalar val = 3. * x0(localmesh->xend + 1, localmesh->ystart) + - 3. * x0(localmesh->xend + 1, localmesh->ystart + 1) + + x0(localmesh->xend + 1, localmesh->ystart + 2); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + // Set the value for the rhs of the boundary condition to zero + val = 0.0; + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + } + + for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { + // Should not go into corner cells, they are treated specially below + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + int ind = globalIndex(it.ind, localmesh->yend + 1); + + // For the boundary value of the initial guess, use the value that would be set by + // applying the boundary condition to the initial guess + PetscScalar val = 3. * x0(it.ind, localmesh->yend) + - 3. * x0(it.ind, localmesh->yend - 1) + + x0(it.ind, localmesh->yend - 2); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + // Set the value for the rhs of the boundary condition to zero + val = 0.0; + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + if (localmesh->hasBndryUpperY()) { + if (localmesh->firstX()) { + int ind = globalIndex(localmesh->xstart - 1, localmesh->yend + 1); + + // For the boundary value of the initial guess, use the value that would be set by + // applying the boundary condition to the initial guess + PetscScalar val = 3. * x0(localmesh->xstart - 1, localmesh->yend) + - 3. * x0(localmesh->xstart - 1, localmesh->yend - 1) + + x0(localmesh->xstart - 1, localmesh->yend - 2); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + // Set the value for the rhs of the boundary condition to zero + val = 0.0; + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + if (localmesh->lastX()) { + int ind = globalIndex(localmesh->xend + 1, localmesh->yend + 1); + + // For the boundary value of the initial guess, use the value that would be set by + // applying the boundary condition to the initial guess + PetscScalar val = 3. * x0(localmesh->xend + 1, localmesh->yend) + - 3. * x0(localmesh->xend + 1, localmesh->yend - 1) + + x0(localmesh->xend + 1, localmesh->yend - 2); + VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); + + // Set the value for the rhs of the boundary condition to zero + val = 0.0; + VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); + } + } + } else { + throw BoutException("Unsupported option for y_bndry"); + } +} + +/*! Preconditioner + * NOTE: For generality, this routine does use globalIndex() in the inner loop, although + * this may be slightly less efficient than incrementing an integer for the global index, + * the finite-volume and finite-difference implementations have slightly different + * indexing patterns, so incrementing an integer would be tricky. + */ +int LaplaceXYpetsc::precon(Vec input, Vec result) { + + for (auto itdwn = localmesh->iterateBndryLowerY(); !itdwn.isDone(); itdwn++) { + // Should not go into corner cells, LaplaceXYpetsc stencil does not include them + if (itdwn.ind < localmesh->xstart or itdwn.ind > localmesh->xend) { + continue; + } + const int ind = globalIndex(itdwn.ind, localmesh->ystart - 1); + PetscScalar val; + VecGetValues(input, 1, &ind, &val); + VecSetValues(result, 1, &ind, &val, INSERT_VALUES); + } + for (auto itup = localmesh->iterateBndryUpperY(); !itup.isDone(); itup++) { + // Should not go into corner cells, LaplaceXYpetsc stencil does not include them + if (itup.ind < localmesh->xstart or itup.ind > localmesh->xend) { + continue; + } + const int ind = globalIndex(itup.ind, localmesh->yend + 1); + PetscScalar val; + VecGetValues(input, 1, &ind, &val); + VecSetValues(result, 1, &ind, &val, INSERT_VALUES); + } + + // Load vector x into bvals array + for (int x = xstart; x <= xend; x++) { + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + const int ind = globalIndex(x, y); + PetscScalar val; + VecGetValues(input, 1, &ind, &val); + bvals(y - localmesh->ystart, x - xstart) = val; + } + } + + // Solve tridiagonal systems using CR solver + cr->solve(bvals, xvals); + + // Save result xvals into y array + for (int x = xstart; x <= xend; x++) { + for (int y = localmesh->ystart; y <= localmesh->yend; y++) { + const int ind = globalIndex(x, y); + PetscScalar val = xvals(y - localmesh->ystart, x - xstart); + VecSetValues(result, 1, &ind, &val, INSERT_VALUES); + } + } + VecAssemblyBegin(result); + VecAssemblyEnd(result); + return 0; +} + +/////////////////////////////////////////////////////////////// + +int LaplaceXYpetsc::localSize() { + + // Bulk of points + const int nx = localmesh->xend - localmesh->xstart + 1; + const int ny = localmesh->yend - localmesh->ystart + 1; + + int n = nx * ny; + + // X boundaries + if (localmesh->firstX()) { + n += ny; + } + if (localmesh->lastX()) { + n += ny; + } + + // Y boundaries + for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { + // Should not go into corner cells, LaplaceXYpetsc stencil does not include them + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + n++; + } + if ((not finite_volume) and localmesh->hasBndryLowerY()) { + if (localmesh->firstX()) { + n++; + } + if (localmesh->lastX()) { + n++; + } + } + for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { + // Should not go into corner cells, LaplaceXYpetsc stencil does not include them + if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { + continue; + } + n++; + } + if ((not finite_volume) and localmesh->hasBndryUpperY()) { + if (localmesh->firstX()) { + n++; + } + if (localmesh->lastX()) { + n++; + } + } + + return n; +} + +int LaplaceXYpetsc::globalIndex(int x, int y) { + if ((x < 0) || (x >= localmesh->LocalNx) || (y < 0) || (y >= localmesh->LocalNy)) { + return -1; // Out of range + } + + // Get the index from a Field2D, round to integer + return static_cast(std::round(indexXY(x, y))); +} + +void LaplaceXYpetsc::savePerformance(Solver& solver, const std::string& name) { + // set flag so that performance monitoring values are calculated + save_performance = true; + + // add values to be saved to the output + if (not name.empty()) { + default_prefix = name; + } + + // add monitor to reset counters/averages for new output timestep + // monitor added to back of queue, so that values are reset after being saved + solver.addMonitor(&monitor, Solver::BACK); +} + +int LaplaceXYpetsc::LaplaceXYMonitor::call(Solver* /*solver*/, BoutReal /*time*/, + int /*iter*/, int /*nout*/) { + laplacexy->output_average_iterations = laplacexy->average_iterations; + + laplacexy->n_calls = 0; + laplacexy->average_iterations = 0.; + + return 0; +} + +void LaplaceXYpetsc::LaplaceXYMonitor::outputVars(Options& output_options, + const std::string& time_dimension) { + output_options[fmt::format("{}_average_iterations", laplacexy->default_prefix)] + .assignRepeat(laplacexy->output_average_iterations, time_dimension); +} + +#endif // BOUT_HAS_PETSC diff --git a/src/invert/laplacexy/impls/petsc/laplacexy-petsc.hxx b/src/invert/laplacexy/impls/petsc/laplacexy-petsc.hxx new file mode 100644 index 0000000000..27d25cf9a7 --- /dev/null +++ b/src/invert/laplacexy/impls/petsc/laplacexy-petsc.hxx @@ -0,0 +1,209 @@ +/************************************************************************** + * Laplacian solver in 2D (X-Y) + * + * Equation solved is: + * + * Div( A * Grad_perp(x) ) + B*x = b + * + * Intended for use in solving n = 0 component of potential + * from inversion of vorticity equation + * + ************************************************************************** + * Copyright 2015 B.Dudson + * + * Contact: Ben Dudson, bd512@york.ac.uk + * + * This file is part of BOUT++. + * + * BOUT++ is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * BOUT++ is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with BOUT++. If not, see . + * + **************************************************************************/ + +#ifndef BOUT_LAPLACE_XY_PETSC_H +#define BOUT_LAPLACE_XY_PETSC_H + +#include "bout/build_defines.hxx" +#include "bout/invert/laplacexy.hxx" + +#if !BOUT_HAS_PETSC + +namespace { +RegisterUnavailableLaplaceXY + registerlaplacexypetsc("petsc", "BOUT++ was not configured with PETSc"); +} + +#else // BOUT_HAS_PETSC + +#include "bout/bout_types.hxx" +#include "bout/cyclic_reduction.hxx" +#include "bout/field2d.hxx" +#include "bout/mesh.hxx" +#include "bout/monitor.hxx" +#include "bout/petsclib.hxx" +#include "bout/utils.hxx" + +#include +#include + +class Options; +class Solver; + +class LaplaceXYpetsc : public LaplaceXY { +public: + LaplaceXYpetsc(Mesh* m = nullptr, Options* opt = nullptr, CELL_LOC loc = CELL_CENTRE); + LaplaceXYpetsc(const LaplaceXYpetsc&) = delete; + LaplaceXYpetsc(LaplaceXYpetsc&&) = delete; + LaplaceXYpetsc& operator=(const LaplaceXYpetsc&) = delete; + LaplaceXYpetsc& operator=(LaplaceXYpetsc&&) = delete; + ~LaplaceXYpetsc() override; + + /*! + * Set coefficients (A, B) in equation: + * Div( A * Grad_perp(x) ) + B*x = b + */ + void setCoefs(const Field2D& A, const Field2D& B) override; + + /*! + * Solve Laplacian in X-Y + * + * Inputs + * ====== + * + * rhs - The field to be inverted. This must be allocated + * and contain valid data. + * x0 - Initial guess at the solution. If this is unallocated + * then an initial guess of zero will be used. + * + * Returns + * ======= + * + * The solution as a Field2D. On failure an exception will be raised + * + */ + Field2D solve(const Field2D& rhs, const Field2D& x0) override; + + /*! + * Preconditioner function + * This is called by PETSc via a static function. + * and should not be called by external users + */ + int precon(Vec x, Vec y); + + /*! + * If this method is called, save some performance monitoring information + */ + void savePerformance(Solver& solver, const std::string& name = ""); + +private: + PetscLib lib; ///< Requires PETSc library + Mat MatA = nullptr; ///< Matrix to be inverted + Vec xs = nullptr; ///< Solution vector + Vec bs = nullptr; ///< RHS vectors + KSP ksp = nullptr; ///< Krylov Subspace solver + PC pc = nullptr; ///< Preconditioner + + Mesh* localmesh; ///< The mesh this operates on, provides metrics and communication + + /// default prefix for writing performance logging variables + std::string default_prefix; + + // Preconditioner + int xstart, xend; + int nloc, nsys; + Matrix acoef, bcoef, ccoef, xvals, bvals; + std::unique_ptr> cr; ///< Tridiagonal solver + + // Use finite volume or finite difference discretization + bool finite_volume{true}; + + // Y derivatives + bool include_y_derivs; // Include Y derivative terms? + + // Boundary conditions + bool x_inner_dirichlet; // Dirichlet on inner X boundary? + bool x_outer_dirichlet{false}; // Dirichlet on outer X boundary? + std::string y_bndry{"neumann"}; // Boundary condition for y-boundary + + // Location of the rhs and solution + CELL_LOC location; + + /*! + * Number of grid points on this processor + */ + int localSize(); + + /*! + * Return the communicator for XY + */ + MPI_Comm communicator(); + + /*! + * Return the global index of a local (x,y) coordinate + * including guard cells. + * Boundary cells have a global index of -1 + * + * To do this, a Field2D (indexXY) is used to store + * the index as a floating point number which is then rounded + * to an integer. Guard cells are filled by communication + * so no additional logic is needed in Mesh. + */ + int globalIndex(int x, int y); + Field2D indexXY; ///< Global index (integer stored as BoutReal) + + // Save performance information? + bool save_performance = false; + + // Running average of number of iterations taken for solve in each output timestep + BoutReal average_iterations = 0.; + + // Variable to store the final result of average_iterations, since output is + // written after all other monitors have been called, and average_iterations + // must be reset in the monitor + BoutReal output_average_iterations = 0.; + + // Running total of number of calls to the solver in each output timestep + int n_calls = 0; + + // Utility methods + void setPreallocationFiniteVolume(PetscInt* d_nnz, PetscInt* o_nnz); + void setPreallocationFiniteDifference(PetscInt* d_nnz, PetscInt* o_nnz); + void setMatrixElementsFiniteVolume(const Field2D& A, const Field2D& B); + void setMatrixElementsFiniteDifference(const Field2D& A, const Field2D& B); + void solveFiniteVolume(const Field2D& x0); + void solveFiniteDifference(const Field2D& x0); + + // Monitor class used to reset performance-monitoring variables for a new + // output timestep + friend class LaplaceXYMonitor; + class LaplaceXYMonitor : public Monitor { + public: + LaplaceXYMonitor(LaplaceXYpetsc& owner) : laplacexy(&owner) {} + + int call(Solver* /*solver*/, BoutReal /*time*/, int /*iter*/, int /*nout*/) override; + void outputVars(Options& output_options, const std::string& time_dimension) override; + + private: + // LaplaceXY object that this monitor belongs to + LaplaceXYpetsc* laplacexy; + }; + + LaplaceXYMonitor monitor; +}; + +namespace { +const inline RegisterLaplaceXY registerlaplacexypetsc{"petsc"}; +} + +#endif // BOUT_HAS_PETSC +#endif // BOUT_LAPLACE_XY_PETSC_H diff --git a/src/invert/laplacexy2/laplacexy2.cxx b/src/invert/laplacexy/impls/petsc2/laplacexy-petsc2.cxx similarity index 90% rename from src/invert/laplacexy2/laplacexy2.cxx rename to src/invert/laplacexy/impls/petsc2/laplacexy-petsc2.cxx index 713fd8e68f..01019466a8 100644 --- a/src/invert/laplacexy2/laplacexy2.cxx +++ b/src/invert/laplacexy/impls/petsc2/laplacexy-petsc2.cxx @@ -2,26 +2,34 @@ #if BOUT_HAS_PETSC and not BOUT_USE_METRIC_3D -#include -#include -#include -#include -#include -#include -#include +#include "laplacexy-petsc2.hxx" + +#include "bout/assert.hxx" +#include "bout/bout_types.hxx" +#include "bout/boutcomm.hxx" +#include "bout/coordinates.hxx" +#include "bout/field2d.hxx" +#include "bout/globalindexer.hxx" +#include "bout/globals.hxx" +#include "bout/operatorstencil.hxx" +#include "bout/options.hxx" +#include "bout/petsc_interface.hxx" +#include "bout/region.hxx" +#include "bout/sys/range.hxx" +#include "bout/sys/timer.hxx" + +#include #include -#include - namespace { Ind2D index2d(Mesh* mesh, int x, int y) { int ny = mesh->LocalNy; - return Ind2D(x * ny + y, ny, 1); + return Ind2D((x * ny) + y, ny, 1); } } // namespace -LaplaceXY2::LaplaceXY2(Mesh* m, Options* opt, const CELL_LOC loc) +LaplaceXYpetsc2::LaplaceXYpetsc2(Mesh* m, Options* opt, const CELL_LOC loc) : localmesh(m == nullptr ? bout::globals::mesh : m), indexConverter(std::make_shared>( localmesh, squareStencil(localmesh))), @@ -100,12 +108,12 @@ LaplaceXY2::LaplaceXY2(Mesh* m, Options* opt, const CELL_LOC loc) // Decide boundary condititions if (localmesh->periodicY(localmesh->xstart)) { // Periodic in Y, so in the core - opt->get("core_bndry_dirichlet", x_inner_dirichlet, false); + x_inner_dirichlet = (*opt)["core_bndry_dirichlet"].withDefault(false); } else { // Non-periodic, so in the PF region - opt->get("pf_bndry_dirichlet", x_inner_dirichlet, true); + x_inner_dirichlet = (*opt)["pf_bndry_dirichlet"].withDefault(true); } - opt->get("y_bndry_dirichlet", y_bndry_dirichlet, false); + y_bndry_dirichlet = (*opt)["y_bndry_dirichlet"].withDefault(false); /////////////////////////////////////////////////// // Including Y derivatives? @@ -120,10 +128,10 @@ LaplaceXY2::LaplaceXY2(Mesh* m, Options* opt, const CELL_LOC loc) Field2D zero(0., localmesh); one.setLocation(location); zero.setLocation(location); - setCoefs(one, zero); + LaplaceXYpetsc2::setCoefs(one, zero); } -void LaplaceXY2::setCoefs(const Field2D& A, const Field2D& B) { +void LaplaceXYpetsc2::setCoefs(const Field2D& A, const Field2D& B) { Timer timer("invert"); ASSERT1(A.getMesh() == localmesh); @@ -290,7 +298,7 @@ void LaplaceXY2::setCoefs(const Field2D& A, const Field2D& B) { #endif } -LaplaceXY2::~LaplaceXY2() { +LaplaceXYpetsc2::~LaplaceXYpetsc2() { PetscBool is_finalised; PetscFinalized(&is_finalised); @@ -300,7 +308,7 @@ LaplaceXY2::~LaplaceXY2() { } } -Field2D LaplaceXY2::solve(const Field2D& rhs, const Field2D& x0) { +Field2D LaplaceXYpetsc2::solve(const Field2D& rhs, const Field2D& x0) { Timer timer("invert"); ASSERT1(rhs.getMesh() == localmesh); @@ -310,7 +318,8 @@ Field2D LaplaceXY2::solve(const Field2D& rhs, const Field2D& x0) { // Load initial guess x0 into xs and rhs into bs - PetscVector xs(x0, indexConverter), bs(rhs, indexConverter); + PetscVector xs(x0, indexConverter); + PetscVector bs(rhs, indexConverter); if (localmesh->firstX()) { if (x_inner_dirichlet) { @@ -381,11 +390,11 @@ Field2D LaplaceXY2::solve(const Field2D& rhs, const Field2D& x0) { // Solve the system KSPSolve(ksp, *bs.get(), *xs.get()); - KSPConvergedReason reason; + KSPConvergedReason reason = KSP_CONVERGED_ITERATING; KSPGetConvergedReason(ksp, &reason); if (reason <= 0) { - throw BoutException("LaplaceXY2 failed to converge. Reason {} ({:d})", + throw BoutException("LaplaceXYpetsc2 failed to converge. Reason {} ({:d})", KSPConvergedReasons[reason], static_cast(reason)); } diff --git a/include/bout/invert/laplacexy2.hxx b/src/invert/laplacexy/impls/petsc2/laplacexy-petsc2.hxx similarity index 61% rename from include/bout/invert/laplacexy2.hxx rename to src/invert/laplacexy/impls/petsc2/laplacexy-petsc2.hxx index 51f75f467d..eb1be636f2 100644 --- a/include/bout/invert/laplacexy2.hxx +++ b/src/invert/laplacexy/impls/petsc2/laplacexy-petsc2.hxx @@ -30,63 +30,47 @@ * **************************************************************************/ -#ifndef BOUT_LAPLACE_XY2_H -#define BOUT_LAPLACE_XY2_H +#ifndef BOUT_LAPLACE_XY_PETSC2_H +#define BOUT_LAPLACE_XY_PETSC2_H #include "bout/build_defines.hxx" - -#if (not BOUT_HAS_PETSC) or BOUT_USE_METRIC_3D -// If no PETSc or 3D metrics - -#warning LaplaceXY2 requires PETSc and 2D metrics. No LaplaceXY2 available - -#include -#include -#include - -/*! - * Create a dummy class so that code will compile - * without PETSc, but will throw an exception if - * LaplaceXY is used. - */ -class LaplaceXY2 { -public: - LaplaceXY2(Mesh* UNUSED(m) = nullptr, Options* UNUSED(opt) = nullptr, - const CELL_LOC UNUSED(loc) = CELL_CENTRE) { - throw BoutException( - "LaplaceXY2 requires PETSc and 2D metrics. No LaplaceXY2 available"); - } - void setCoefs(const Field2D& UNUSED(A), const Field2D& UNUSED(B)) {} - Field2D solve(const Field2D& UNUSED(rhs), const Field2D& UNUSED(x0)) { - throw BoutException( - "LaplaceXY2 requires PETSc and 2D metrics. No LaplaceXY2 available"); - } -}; - -#else // BOUT_HAS_PETSC and 2D metrics - -#include "bout/utils.hxx" -#include -#include -#include -#include - -class LaplaceXY2 { +#include "bout/invert/laplacexy.hxx" + +#if !BOUT_HAS_PETSC +namespace { +RegisterUnavailableLaplaceXY + registerlaplacexypetsc2("petsc2", "BOUT++ was not configured with PETSc"); +} + +#elif BOUT_USE_METRIC_3D +namespace { +RegisterUnavailableLaplaceXY + registerlaplacexypetsc2("petsc2", "BOUT++ was configured with 3D metrics"); +} + +#else // BOUT_HAS_PETSC + +#include "bout/bout_types.hxx" +#include "bout/field2d.hxx" +#include "bout/globalindexer.hxx" +#include "bout/mesh.hxx" +#include "bout/petsc_interface.hxx" +#include "bout/petsclib.hxx" + +class LaplaceXYpetsc2 : public LaplaceXY { public: - /*! - * Constructor - */ - LaplaceXY2(Mesh* m = nullptr, Options* opt = nullptr, const CELL_LOC loc = CELL_CENTRE); - /*! - * Destructor - */ - ~LaplaceXY2(); + LaplaceXYpetsc2(Mesh* m = nullptr, Options* opt = nullptr, CELL_LOC loc = CELL_CENTRE); + LaplaceXYpetsc2(const LaplaceXYpetsc2&) = default; + LaplaceXYpetsc2(LaplaceXYpetsc2&&) = delete; + LaplaceXYpetsc2& operator=(const LaplaceXYpetsc2&) = default; + LaplaceXYpetsc2& operator=(LaplaceXYpetsc2&&) = delete; + ~LaplaceXYpetsc2() override; /*! * Set coefficients (A, B) in equation: * Div( A * Grad_perp(x) ) + B*x = b */ - void setCoefs(const Field2D& A, const Field2D& B); + void setCoefs(const Field2D& A, const Field2D& B) override; /*! * Solve Laplacian in X-Y @@ -105,7 +89,7 @@ public: * The solution as a Field2D. On failure an exception will be raised * */ - Field2D solve(const Field2D& rhs, const Field2D& x0); + Field2D solve(const Field2D& rhs, const Field2D& x0) override; /*! * Preconditioner function @@ -121,8 +105,8 @@ private: IndexerPtr indexConverter; PetscMatrix matrix; ///< Matrix to be inverted - KSP ksp; ///< Krylov Subspace solver - PC pc; ///< Preconditioner + KSP ksp = nullptr; ///< Krylov Subspace solver + PC pc = nullptr; ///< Preconditioner // Y derivatives bool include_y_derivs; // Include Y derivative terms? @@ -140,5 +124,9 @@ private: MPI_Comm communicator(); }; +namespace { +const inline RegisterLaplaceXY registerlaplacexypetsc2{"petsc2"}; +} + #endif // BOUT_HAS_PETSC -#endif // BOUT_LAPLACE_XY2_H +#endif // BOUT_LAPLACE_XY_PETSC2_H diff --git a/src/invert/laplacexy/laplacexy.cxx b/src/invert/laplacexy/laplacexy.cxx index bdd22e29bd..3bac2d9d2b 100644 --- a/src/invert/laplacexy/laplacexy.cxx +++ b/src/invert/laplacexy/laplacexy.cxx @@ -1,1889 +1,10 @@ -#include "bout/build_defines.hxx" - -#if BOUT_HAS_PETSC +// NOLINTBEGIN(misc-include-cleaner, unused-includes) +#include "impls/hypre/laplacexy-hypre.hxx" +#include "impls/petsc/laplacexy-petsc.hxx" +#include "impls/petsc2/laplacexy-petsc2.hxx" +// NOLINTEND(misc-include-cleaner, unused-includes) #include -#include -#include -#include -#include -#include -#include -#include - -#include - -#include - -static PetscErrorCode laplacePCapply(PC pc, Vec x, Vec y) { - int ierr; - - // Get the context - LaplaceXY* s; - ierr = PCShellGetContext(pc, reinterpret_cast(&s)); - CHKERRQ(ierr); - - PetscFunctionReturn(s->precon(x, y)); -} - -LaplaceXY::LaplaceXY(Mesh* m, Options* opt, const CELL_LOC loc) - : lib(opt == nullptr ? &(Options::root()["laplacexy"]) : opt), - localmesh(m == nullptr ? bout::globals::mesh : m), location(loc), monitor(*this) { - Timer timer("invert"); - - if (opt == nullptr) { - // If no options supplied, use default - opt = &(Options::root()["laplacexy"]); - } - - finite_volume = - (*opt)["finite_volume"] - .doc("Use finite volume rather than finite difference discretisation.") - .withDefault(true); - - /////////////////////////////////////////////////// - // Boundary condititions options - if (localmesh->periodicY(localmesh->xstart)) { - // Periodic in Y, so in the core - x_inner_dirichlet = (*opt)["core_bndry_dirichlet"].withDefault(false); - } else { - // Non-periodic, so in the PF region - x_inner_dirichlet = (*opt)["pf_bndry_dirichlet"].withDefault(true); - } - if ((*opt)["y_bndry_dirichlet"].isSet()) { - throw BoutException("y_bndry_dirichlet has been deprecated. Use y_bndry=dirichlet " - "instead."); - } else { - y_bndry = (*opt)["y_bndry"].withDefault("neumann"); - } - - // Check value of y_bndry is a supported option - if (not(y_bndry == "dirichlet" or y_bndry == "neumann" or y_bndry == "free_o3")) { - - throw BoutException("Unrecognized option '{}' for laplacexy:ybndry", y_bndry); - } - - if (not finite_volume) { - // Check we can use corner cells - if (not localmesh->include_corner_cells) { - throw BoutException( - "Finite difference form of LaplaceXY allows non-orthogonal x- and " - "y-directions, so requires mesh:include_corner_cells=true."); - } - } - - // Use name of options section as the default prefix for performance logging variables - default_prefix = opt->name(); - - // Get MPI communicator - auto comm = BoutComm::get(); - - // Local size - const int localN = localSize(); - - // Create Vectors - VecCreate(comm, &xs); - VecSetSizes(xs, localN, PETSC_DETERMINE); - VecSetFromOptions(xs); - VecDuplicate(xs, &bs); - - // Set size of Matrix on each processor to localN x localN - MatCreate(comm, &MatA); - MatSetSizes(MatA, localN, localN, PETSC_DETERMINE, PETSC_DETERMINE); - MatSetFromOptions(MatA); - - ////////////////////////////////////////////////// - // Specify local indices. This creates a mapping - // from local indices to index, using a Field2D object - - indexXY = -1; // Set all points to -1, indicating out of domain - - int ind = 0; - - // Y boundaries - for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { - // Should not go into corner cells, LaplaceXY stencil does not include them - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - indexXY(it.ind, localmesh->ystart - 1) = ind++; - } - if ((not finite_volume) and localmesh->hasBndryLowerY()) { - // Corner boundary cells - if (localmesh->firstX()) { - indexXY(localmesh->xstart - 1, localmesh->ystart - 1) = ind++; - } - if (localmesh->lastX()) { - indexXY(localmesh->xend + 1, localmesh->ystart - 1) = ind++; - } - } - for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { - // Should not go into corner cells, LaplaceXY stencil does not include them - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - indexXY(it.ind, localmesh->yend + 1) = ind++; - } - if ((not finite_volume) and localmesh->hasBndryUpperY()) { - // Corner boundary cells - if (localmesh->firstX()) { - indexXY(localmesh->xstart - 1, localmesh->yend + 1) = ind++; - } - if (localmesh->lastX()) { - indexXY(localmesh->xend + 1, localmesh->yend + 1) = ind++; - } - } - - xstart = localmesh->xstart; - if (localmesh->firstX()) { - xstart -= 1; // Include X guard cells - } - xend = localmesh->xend; - if (localmesh->lastX()) { - xend += 1; - } - for (int x = xstart; x <= xend; x++) { - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - indexXY(x, y) = ind++; - } - } - - ASSERT1(ind == localN); // Reached end of range - - ////////////////////////////////////////////////// - // Allocate storage for preconditioner - - nloc = xend - xstart + 1; // Number of X points on this processor - nsys = localmesh->yend - localmesh->ystart + 1; // Number of separate Y slices - - acoef.reallocate(nsys, nloc); - bcoef.reallocate(nsys, nloc); - ccoef.reallocate(nsys, nloc); - xvals.reallocate(nsys, nloc); - bvals.reallocate(nsys, nloc); - - // Create a cyclic reduction object - cr = bout::utils::make_unique>(localmesh->getXcomm(), nloc); - - ////////////////////////////////////////////////// - // Pre-allocate PETSc storage - - PetscInt *d_nnz, *o_nnz; - PetscMalloc((localN) * sizeof(PetscInt), &d_nnz); - PetscMalloc((localN) * sizeof(PetscInt), &o_nnz); - - if (finite_volume) { - setPreallocationFiniteVolume(d_nnz, o_nnz); - } else { - setPreallocationFiniteDifference(d_nnz, o_nnz); - } - // Pre-allocate - MatMPIAIJSetPreallocation(MatA, 0, d_nnz, 0, o_nnz); - MatSetUp(MatA); - - PetscFree(d_nnz); - PetscFree(o_nnz); - - // Determine which row/columns of the matrix are locally owned - int Istart, Iend; - MatGetOwnershipRange(MatA, &Istart, &Iend); - - // Convert indexXY from local index to global index - indexXY += Istart; - - // Now communicate to fill guard cells - // Note, this includes corner cells if necessary - localmesh->communicate(indexXY); - - ////////////////////////////////////////////////// - // Set up KSP - - // Declare KSP Context - KSPCreate(comm, &ksp); - - // Configure Linear Solver - - const bool direct = (*opt)["direct"].doc("Use a direct LU solver").withDefault(false); - - if (direct) { - KSPGetPC(ksp, &pc); - PCSetType(pc, PCLU); -#if PETSC_VERSION_GE(3, 9, 0) - PCFactorSetMatSolverType(pc, "mumps"); -#else - PCFactorSetMatSolverPackage(pc, "mumps"); -#endif - } else { - - // Convergence Parameters. Solution is considered converged if |r_k| < max( rtol * |b| , atol ) - // where r_k = b - Ax_k. The solution is considered diverged if |r_k| > dtol * |b|. - - const BoutReal rtol = (*opt)["rtol"].doc("Relative tolerance").withDefault(1e-5); - const BoutReal atol = - (*opt)["atol"] - .doc("Absolute tolerance. The solution is considered converged if |Ax-b| " - "< max( rtol * |b| , atol )") - .withDefault(1e-10); - const BoutReal dtol = - (*opt)["dtol"] - .doc("The solution is considered diverged if |Ax-b| > dtol * |b|") - .withDefault(1e3); - const int maxits = (*opt)["maxits"].doc("Maximum iterations").withDefault(100000); - - // Get KSP Solver Type - const std::string ksptype = - (*opt)["ksptype"].doc("KSP solver type").withDefault("gmres"); - - // Get PC type - const std::string pctype = - (*opt)["pctype"].doc("Preconditioner type").withDefault("none"); - - KSPSetType(ksp, ksptype.c_str()); - KSPSetTolerances(ksp, rtol, atol, dtol, maxits); - - KSPSetInitialGuessNonzero(ksp, static_cast(true)); - - KSPGetPC(ksp, &pc); - PCSetType(pc, pctype.c_str()); - - if (pctype == "shell") { - // Using tridiagonal solver as preconditioner - PCShellSetApply(pc, laplacePCapply); - PCShellSetContext(pc, this); - - const bool rightprec = - (*opt)["rightprec"].doc("Use right preconditioning?").withDefault(true); - if (rightprec) { - KSPSetPCSide(ksp, PC_RIGHT); // Right preconditioning - } else { - KSPSetPCSide(ksp, PC_LEFT); // Left preconditioning - } - } - } - - lib.setOptionsFromInputFile(ksp); - - /////////////////////////////////////////////////// - // Including Y derivatives? - - include_y_derivs = (*opt)["include_y_derivs"] - .doc("Include Y derivatives in operator to invert?") - .withDefault(true); - - /////////////////////////////////////////////////// - // Set the default coefficients - Field2D one(1., localmesh); - Field2D zero(0., localmesh); - one.setLocation(location); - zero.setLocation(location); - setCoefs(one, zero); -} - -void LaplaceXY::setPreallocationFiniteVolume(PetscInt* d_nnz, PetscInt* o_nnz) { - const int localN = localSize(); - - // This discretisation uses a 5-point stencil - for (int i = 0; i < localN; i++) { - // Non-zero elements on this processor - d_nnz[i] = 5; // Star pattern in 2D - // Non-zero elements on neighboring processor - o_nnz[i] = 0; - } - - // X boundaries - if (localmesh->firstX()) { - // Lower X boundary - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - const int localIndex = globalIndex(localmesh->xstart - 1, y); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - - d_nnz[localIndex] = 2; // Diagonal sub-matrix - o_nnz[localIndex] = 0; // Off-diagonal sub-matrix - } - } else { - // On another processor - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - const int localIndex = globalIndex(localmesh->xstart, y); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - d_nnz[localIndex] -= 1; - o_nnz[localIndex] += 1; - } - } - if (localmesh->lastX()) { - // Upper X boundary - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - const int localIndex = globalIndex(localmesh->xend + 1, y); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - d_nnz[localIndex] = 2; // Diagonal sub-matrix - o_nnz[localIndex] = 0; // Off-diagonal sub-matrix - } - } else { - // On another processor - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - const int localIndex = globalIndex(localmesh->xend, y); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - d_nnz[localIndex] -= 1; - o_nnz[localIndex] += 1; - } - } - // Y boundaries - - for (int x = localmesh->xstart; x <= localmesh->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 - { - const int localIndex = globalIndex(x, localmesh->ystart); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - // d_nnz[localIndex] -= 1; // Note: Slightly inefficient - o_nnz[localIndex] += 1; - } - { - const int localIndex = globalIndex(x, localmesh->yend); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - // d_nnz[localIndex] -= 1; // Note: Slightly inefficient - o_nnz[localIndex] += 1; - } - } - - for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { - // Should not go into corner cells, LaplaceXY stencil does not include them - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - { - const int localIndex = globalIndex(it.ind, localmesh->ystart - 1); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - d_nnz[localIndex] = 2; // Diagonal sub-matrix - o_nnz[localIndex] = 0; // Off-diagonal sub-matrix - } - { - const int localIndex = globalIndex(it.ind, localmesh->ystart); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - d_nnz[localIndex] += 1; - o_nnz[localIndex] -= 1; - } - } - for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { - // Should not go into corner cells, LaplaceXY stencil does not include them - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - { - const int localIndex = globalIndex(it.ind, localmesh->yend + 1); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - d_nnz[localIndex] = 2; // Diagonal sub-matrix - o_nnz[localIndex] = 0; // Off-diagonal sub-matrix - } - { - const int localIndex = globalIndex(it.ind, localmesh->yend); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - d_nnz[localIndex] += 1; - o_nnz[localIndex] -= 1; - } - } -} - -void LaplaceXY::setPreallocationFiniteDifference(PetscInt* d_nnz, PetscInt* o_nnz) { - const int localN = localSize(); - - // This discretisation uses a 9-point stencil - for (int i = 0; i < localN; i++) { - // Non-zero elements on this processor - d_nnz[i] = 9; // Square pattern in 2D - // Non-zero elements on neighboring processor - o_nnz[i] = 0; - } - - // X boundaries - if (localmesh->firstX()) { - // Lower X boundary - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - const int localIndex = globalIndex(localmesh->xstart - 1, y); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - - d_nnz[localIndex] = 2; // Diagonal sub-matrix - o_nnz[localIndex] = 0; // Off-diagonal sub-matrix - } - } else { - // On another processor - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - const int localIndex = globalIndex(localmesh->xstart, y); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - d_nnz[localIndex] -= 3; - o_nnz[localIndex] += 3; - } - } - if (localmesh->lastX()) { - // Upper X boundary - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - const int localIndex = globalIndex(localmesh->xend + 1, y); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - d_nnz[localIndex] = 2; // Diagonal sub-matrix - o_nnz[localIndex] = 0; // Off-diagonal sub-matrix - } - } else { - // On another processor - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - const int localIndex = globalIndex(localmesh->xend, y); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - d_nnz[localIndex] -= 3; - o_nnz[localIndex] += 3; - } - } - // Y boundaries - const int y_bndry_stencil_size = (y_bndry == "free_o3") ? 4 : 2; - - for (int x = localmesh->xstart; x <= localmesh->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 - { - const int localIndex = globalIndex(x, localmesh->ystart); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - // d_nnz[localIndex] -= 3; // Note: Slightly inefficient - o_nnz[localIndex] += 3; - } - { - const int localIndex = globalIndex(x, localmesh->yend); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - // d_nnz[localIndex] -= 3; // Note: Slightly inefficient - o_nnz[localIndex] += 3; - } - } - - for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { - // Should not go into corner cells, they are handled specially below - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - { - const int localIndex = globalIndex(it.ind, localmesh->ystart - 1); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - d_nnz[localIndex] = y_bndry_stencil_size; // Diagonal sub-matrix - o_nnz[localIndex] = 0; // Off-diagonal sub-matrix - } - { - const int localIndex = globalIndex(it.ind, localmesh->ystart); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - //d_nnz[localIndex] += 3; - o_nnz[localIndex] -= 3; - } - } - if (localmesh->hasBndryLowerY()) { - if (y_bndry == "dirichlet") { - // special handling for the corners, since we use a free_o3 y-boundary - // condition just in the corners when y_bndry=="dirichlet" - if (localmesh->firstX()) { - const int localIndex = globalIndex(localmesh->xstart - 1, localmesh->ystart - 1); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - d_nnz[localIndex] = 4; - } - if (localmesh->lastX()) { - const int localIndex = globalIndex(localmesh->xend + 1, localmesh->ystart - 1); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - d_nnz[localIndex] = 4; - } - } else { - if (localmesh->firstX()) { - const int localIndex = globalIndex(localmesh->xstart - 1, localmesh->ystart - 1); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - d_nnz[localIndex] = y_bndry_stencil_size; - } - if (localmesh->lastX()) { - const int localIndex = globalIndex(localmesh->xend + 1, localmesh->ystart - 1); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - d_nnz[localIndex] = y_bndry_stencil_size; - } - } - } - for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { - // Should not go into corner cells, they are handled specially below - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - { - const int localIndex = globalIndex(it.ind, localmesh->yend + 1); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - d_nnz[localIndex] = y_bndry_stencil_size; // Diagonal sub-matrix - o_nnz[localIndex] = 0; // Off-diagonal sub-matrix - } - { - const int localIndex = globalIndex(it.ind, localmesh->yend); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - //d_nnz[localIndex] += 3; - o_nnz[localIndex] -= 3; - } - } - if (localmesh->hasBndryUpperY()) { - if (y_bndry == "dirichlet") { - // special handling for the corners, since we use a free_o3 y-boundary - // condition just in the corners when y_bndry=="dirichlet" - if (localmesh->firstX()) { - const int localIndex = globalIndex(localmesh->xstart - 1, localmesh->yend + 1); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - d_nnz[localIndex] = 4; - } - if (localmesh->lastX()) { - const int localIndex = globalIndex(localmesh->xend + 1, localmesh->yend + 1); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - d_nnz[localIndex] = 4; - } - } else { - if (localmesh->firstX()) { - const int localIndex = globalIndex(localmesh->xstart - 1, localmesh->yend + 1); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - d_nnz[localIndex] = y_bndry_stencil_size; - } - if (localmesh->lastX()) { - const int localIndex = globalIndex(localmesh->xend + 1, localmesh->yend + 1); - ASSERT1((localIndex >= 0) && (localIndex < localN)); - d_nnz[localIndex] = y_bndry_stencil_size; - } - } - } -} - -void LaplaceXY::setCoefs(const Field2D& A, const Field2D& B) { - Timer timer("invert"); - - ASSERT1(A.getMesh() == localmesh); - ASSERT1(B.getMesh() == localmesh); - ASSERT1(A.getLocation() == location); - ASSERT1(B.getLocation() == location); - - if (finite_volume) { - setMatrixElementsFiniteVolume(A, B); - } else { - setMatrixElementsFiniteDifference(A, B); - } - - // X boundaries - if (localmesh->firstX()) { - if (x_inner_dirichlet) { - - // Dirichlet on inner X boundary - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - int row = globalIndex(localmesh->xstart - 1, y); - PetscScalar val = 0.5; - MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); - - int col = globalIndex(localmesh->xstart, y); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - - // Preconditioner - bcoef(y - localmesh->ystart, 0) = 0.5; - ccoef(y - localmesh->ystart, 0) = 0.5; - } - - } else { - - // Neumann on inner X boundary - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - int row = globalIndex(localmesh->xstart - 1, y); - PetscScalar val = 1.0; - MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); - - int col = globalIndex(localmesh->xstart, y); - val = -1.0; - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - - // Preconditioner - bcoef(y - localmesh->ystart, 0) = 1.0; - ccoef(y - localmesh->ystart, 0) = -1.0; - } - } - } - if (localmesh->lastX()) { - // Dirichlet on outer X boundary - - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - int row = globalIndex(localmesh->xend + 1, y); - PetscScalar val = 0.5; - MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); - - int col = globalIndex(localmesh->xend, y); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - - // Preconditioner - acoef(y - localmesh->ystart, localmesh->xend + 1 - xstart) = 0.5; - bcoef(y - localmesh->ystart, localmesh->xend + 1 - xstart) = 0.5; - } - } - - if (y_bndry == "dirichlet") { - // Dirichlet on Y boundaries - for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { - // Should not go into corner cells, LaplaceXY stencil does not include them - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - int row = globalIndex(it.ind, localmesh->ystart - 1); - PetscScalar val = 0.5; - MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); - - int col = globalIndex(it.ind, localmesh->ystart); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - } - - for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { - // Should not go into corner cells, LaplaceXY stencil does not include them - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - int row = globalIndex(it.ind, localmesh->yend + 1); - PetscScalar val = 0.5; - MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); - - int col = globalIndex(it.ind, localmesh->yend); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - } - } else if (y_bndry == "neumann") { - // Neumann on Y boundaries - for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { - // Should not go into corner cells, LaplaceXY stencil does not include them - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - int row = globalIndex(it.ind, localmesh->ystart - 1); - PetscScalar val = 1.0; - MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); - - val = -1.0; - int col = globalIndex(it.ind, localmesh->ystart); - - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - } - - for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { - // Should not go into corner cells, LaplaceXY stencil does not include them - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - int row = globalIndex(it.ind, localmesh->yend + 1); - PetscScalar val = 1.0; - MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); - - val = -1.0; - int col = globalIndex(it.ind, localmesh->yend); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - } - } else if (y_bndry == "free_o3") { - // 'free_o3' extrapolating boundary condition on Y boundaries - for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { - // Should not go into corner cells, LaplaceXY stencil does not include them - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - int row = globalIndex(it.ind, localmesh->ystart - 1); - PetscScalar val = 1.0; - MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); - - val = -3.0; - int col = globalIndex(it.ind, localmesh->ystart); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - - val = 3.0; - col = globalIndex(it.ind, localmesh->ystart + 1); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - - val = -1.0; - col = globalIndex(it.ind, localmesh->ystart + 2); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - } - - for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { - // Should not go into corner cells, LaplaceXY stencil does not include them - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - int row = globalIndex(it.ind, localmesh->yend + 1); - PetscScalar val = 1.0; - MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); - - val = -3.0; - int col = globalIndex(it.ind, localmesh->yend); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - - val = 3.0; - col = globalIndex(it.ind, localmesh->yend - 1); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - - val = -1.0; - col = globalIndex(it.ind, localmesh->yend - 2); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - } - } else { - throw BoutException("Unsupported option for y_bndry"); - } - - if (not finite_volume) { - // Handle corner boundary cells in case we need to include D2DXDY - // Apply the y-boundary-condition to the cells in the x-boundary - this is an - // arbitrary choice, cf. connections around the X-point - - if (localmesh->hasBndryLowerY()) { - if (localmesh->firstX()) { - if (y_bndry == "neumann") { - // Neumann y-bc - // f(xs-1,ys-1) = f(xs-1,ys) - PetscScalar val = 1.0; - int row = globalIndex(localmesh->xstart - 1, localmesh->ystart - 1); - MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); - - val = -1.0; - int col = globalIndex(localmesh->xstart - 1, localmesh->ystart); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - } else if (y_bndry == "free_o3" or y_bndry == "dirichlet") { - // 'free_o3' extrapolating boundary condition on Y boundaries - // f(xs-1,ys-1) = 3*f(xs-1,ys) - 3*f(xs-1,ys+1) + f(xs-1,ys+2) - // - // Use free_o3 at the corners for Dirichlet y-boundaries because we don't know - // what value to pass for the corner - PetscScalar val = 1.0; - int row = globalIndex(localmesh->xstart - 1, localmesh->ystart - 1); - MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); - - val = -3.0; - int col = globalIndex(localmesh->xstart - 1, localmesh->ystart); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - - val = 3.0; - col = globalIndex(localmesh->xstart - 1, localmesh->ystart + 1); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - - val = -1.0; - col = globalIndex(localmesh->xstart - 1, localmesh->ystart + 2); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - } else { - throw BoutException("Unsupported option for y_bndry"); - } - } - if (localmesh->lastX()) { - if (y_bndry == "neumann") { - // Neumann y-bc - // f(xe+1,ys-1) = f(xe+1,ys) - PetscScalar val = 1.0; - int row = globalIndex(localmesh->xend + 1, localmesh->ystart - 1); - MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); - - val = -1.0; - int col = globalIndex(localmesh->xend + 1, localmesh->ystart); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - } else if (y_bndry == "free_o3" or y_bndry == "dirichlet") { - // 'free_o3' extrapolating boundary condition on Y boundaries - // f(xe+1,ys-1) = 3*f(xe+1,ys) - 3*f(xe+1,ys+1) + f(xe+1,ys+2) - // - // Use free_o3 at the corners for Dirichlet y-boundaries because we don't know - // what value to pass for the corner - PetscScalar val = 1.0; - int row = globalIndex(localmesh->xend + 1, localmesh->ystart - 1); - MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); - - val = -3.0; - int col = globalIndex(localmesh->xend + 1, localmesh->ystart); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - - val = 3.0; - col = globalIndex(localmesh->xend + 1, localmesh->ystart + 1); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - - val = -1.0; - col = globalIndex(localmesh->xend + 1, localmesh->ystart + 2); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - } else { - throw BoutException("Unsupported option for y_bndry"); - } - } - } - if (localmesh->hasBndryUpperY()) { - if (localmesh->firstX()) { - if (y_bndry == "neumann") { - // Neumann y-bc - // f(xs-1,ys-1) = f(xs-1,ys) - PetscScalar val = 1.0; - int row = globalIndex(localmesh->xstart - 1, localmesh->yend + 1); - MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); - - val = -1.0; - int col = globalIndex(localmesh->xstart - 1, localmesh->yend); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - } else if (y_bndry == "free_o3" or y_bndry == "dirichlet") { - // 'free_o3' extrapolating boundary condition on Y boundaries - // f(xs-1,ys-1) = 3*f(xs-1,ys) - 3*f(xs-1,ys+1) + f(xs-1,ys+2) - // - // Use free_o3 at the corners for Dirichlet y-boundaries because we don't know - // what value to pass for the corner - PetscScalar val = 1.0; - int row = globalIndex(localmesh->xstart - 1, localmesh->yend + 1); - MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); - - val = -3.0; - int col = globalIndex(localmesh->xstart - 1, localmesh->yend); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - - val = 3.0; - col = globalIndex(localmesh->xstart - 1, localmesh->yend - 1); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - - val = -1.0; - col = globalIndex(localmesh->xstart - 1, localmesh->yend - 2); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - } else { - throw BoutException("Unsupported option for y_bndry"); - } - } - if (localmesh->lastX()) { - if (y_bndry == "neumann") { - // Neumann y-bc - // f(xe+1,ys-1) = f(xe+1,ys) - PetscScalar val = 1.0; - int row = globalIndex(localmesh->xend + 1, localmesh->yend + 1); - MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); - - val = -1.0; - int col = globalIndex(localmesh->xend + 1, localmesh->yend); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - } else if (y_bndry == "free_o3" or y_bndry == "dirichlet") { - // 'free_o3' extrapolating boundary condition on Y boundaries - // f(xe+1,ys-1) = 3*f(xe+1,ys) - 3*f(xe+1,ys+1) + f(xe+1,ys+2) - // - // Use free_o3 at the corners for Dirichlet y-boundaries because we don't know - // what value to pass for the corner - PetscScalar val = 1.0; - int row = globalIndex(localmesh->xend + 1, localmesh->yend + 1); - MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES); - - val = -3.0; - int col = globalIndex(localmesh->xend + 1, localmesh->yend); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - - val = 3.0; - col = globalIndex(localmesh->xend + 1, localmesh->yend - 1); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - - val = -1.0; - col = globalIndex(localmesh->xend + 1, localmesh->yend - 2); - MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES); - } else { - throw BoutException("Unsupported option for y_bndry"); - } - } - } - } - - // Assemble Matrix - MatAssemblyBegin(MatA, MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(MatA, MAT_FINAL_ASSEMBLY); - - // Set the operator -#if PETSC_VERSION_GE(3, 5, 0) - KSPSetOperators(ksp, MatA, MatA); -#else - KSPSetOperators(ksp, MatA, MatA, DIFFERENT_NONZERO_PATTERN); -#endif - - // Set coefficients for preconditioner - cr->setCoefs(acoef, bcoef, ccoef); -} - -void LaplaceXY::setMatrixElementsFiniteVolume(const Field2D& A, const Field2D& B) { - ////////////////////////////////////////////////// - // Set Matrix elements - // - // (1/J) d/dx ( J * g11 d/dx ) + (1/J) d/dy ( J * g22 d/dy ) - - auto coords = localmesh->getCoordinates(location); - const Field2D J_DC = DC(coords->J); - const Field2D g11_DC = DC(coords->g11); - const Field2D dx_DC = DC(coords->dx); - const Field2D dy_DC = DC(coords->dy); - const Field2D g_22_DC = DC(coords->g_22); - const Field2D g_23_DC = DC(coords->g_23); - const Field2D g23_DC = DC(coords->g23); - - for (int x = localmesh->xstart; x <= localmesh->xend; x++) { - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - // stencil entries - PetscScalar c, xm, xp, ym, yp; - - // XX component - - // Metrics on x+1/2 boundary - BoutReal J = 0.5 * (J_DC(x, y) + J_DC(x + 1, y)); - BoutReal g11 = 0.5 * (g11_DC(x, y) + g11_DC(x + 1, y)); - BoutReal dx = 0.5 * (dx_DC(x, y) + dx_DC(x + 1, y)); - BoutReal Acoef = 0.5 * (A(x, y) + A(x + 1, y)); - - BoutReal val = Acoef * J * g11 / (J_DC(x, y) * dx * dx_DC(x, y)); - xp = val; - c = -val; - - // Metrics on x-1/2 boundary - J = 0.5 * (J_DC(x, y) + J_DC(x - 1, y)); - g11 = 0.5 * (g11_DC(x, y) + g11_DC(x - 1, y)); - dx = 0.5 * (dx_DC(x, y) + dx_DC(x - 1, y)); - Acoef = 0.5 * (A(x, y) + A(x - 1, y)); - - val = Acoef * J * g11 / (J_DC(x, y) * dx * dx_DC(x, y)); - xm = val; - c -= val; - - c += B(x, y); - - // Put values into the preconditioner, X derivatives only - acoef(y - localmesh->ystart, x - xstart) = xm; - bcoef(y - localmesh->ystart, x - xstart) = c; - ccoef(y - localmesh->ystart, x - xstart) = xp; - - if (include_y_derivs) { - // YY component - // Metrics at y+1/2 - J = 0.5 * (J_DC(x, y) + J_DC(x, y + 1)); - BoutReal g_22 = 0.5 * (g_22_DC(x, y) + g_22_DC(x, y + 1)); - BoutReal g23 = 0.5 * (g23_DC(x, y) + g23_DC(x, y + 1)); - BoutReal g_23 = 0.5 * (g_23_DC(x, y) + g_23_DC(x, y + 1)); - BoutReal dy = 0.5 * (dy_DC(x, y) + dy_DC(x, y + 1)); - Acoef = 0.5 * (A(x, y + 1) + A(x, y)); - - val = -Acoef * J * g23 * g_23 / (g_22 * J_DC(x, y) * dy * dy_DC(x, y)); - yp = val; - c -= val; - - // Metrics at y-1/2 - J = 0.5 * (J_DC(x, y) + J_DC(x, y - 1)); - g_22 = 0.5 * (g_22_DC(x, y) + g_22_DC(x, y - 1)); - g23 = 0.5 * (g23_DC(x, y) + g23_DC(x, y - 1)); - g_23 = 0.5 * (g_23_DC(x, y) + g_23_DC(x, y - 1)); - dy = 0.5 * (dy_DC(x, y) + dy_DC(x, y - 1)); - Acoef = 0.5 * (A(x, y - 1) + A(x, y)); - - val = -Acoef * J * g23 * g_23 / (g_22 * J_DC(x, y) * dy * dy_DC(x, y)); - ym = val; - c -= val; - } - - ///////////////////////////////////////////////// - // Now have a 5-point stencil for the Laplacian - - int row = globalIndex(x, y); - - // Set the centre (diagonal) - MatSetValues(MatA, 1, &row, 1, &row, &c, INSERT_VALUES); - - // X + 1 - int col = globalIndex(x + 1, y); - MatSetValues(MatA, 1, &row, 1, &col, &xp, INSERT_VALUES); - - // X - 1 - col = globalIndex(x - 1, y); - MatSetValues(MatA, 1, &row, 1, &col, &xm, INSERT_VALUES); - - if (include_y_derivs) { - // Y + 1 - col = globalIndex(x, y + 1); - MatSetValues(MatA, 1, &row, 1, &col, &yp, INSERT_VALUES); - - // Y - 1 - col = globalIndex(x, y - 1); - MatSetValues(MatA, 1, &row, 1, &col, &ym, INSERT_VALUES); - } - } - } -} - -void LaplaceXY::setMatrixElementsFiniteDifference(const Field2D& A, const Field2D& B) { - ////////////////////////////////////////////////// - // Set Matrix elements - // - // Div(A Grad(f)) + B f - // = A Laplace_perp(f) + Grad_perp(A).Grad_perp(f) + B f - // = A*(G1*dfdx + (G2-1/J*d/dy(J/g_22))*dfdy - // + g11*d2fdx2 + (g22-1/g_22)*d2fdy2 + 2*g12*d2fdxdy) - // + g11*dAdx*dfdx + (g22-1/g_22)*dAdy*dfdy + g12*(dAdx*dfdy + dAdy*dfdx) - // + B*f - - auto coords = localmesh->getCoordinates(location); - const Field2D G1_2D = DC(coords->G1); - const Field2D G2_2D = DC(coords->G2); - const Field2D J_2D = DC(coords->J); - const Field2D g11_2D = DC(coords->g11); - const Field2D g_22_2D = DC(coords->g_22); - const Field2D g22_2D = DC(coords->g22); - const Field2D g12_2D = DC(coords->g12); - const Field2D d1_dx_2D = DC(coords->d1_dx); - const Field2D d1_dy_2D = DC(coords->d1_dy); - const Field2D dx_2D = DC(coords->dx); - const Field2D dy_2D = DC(coords->dy); - - const Field2D coef_dfdy = G2_2D - DC(DDY(J_2D / g_22_2D) / J_2D); - - for (int x = localmesh->xstart; x <= localmesh->xend; x++) { - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - // stencil entries - PetscScalar c, xm, xp, ym, yp, xpyp, xpym, xmyp, xmym; - - BoutReal dx = dx_2D(x, y); - - // A*G1*dfdx - BoutReal val = A(x, y) * G1_2D(x, y) / (2. * dx); - xp = val; - xm = -val; - - // A*g11*d2fdx2 - val = A(x, y) * g11_2D(x, y) / SQ(dx); - xp += val; - c = -2. * val; - xm += val; - // Non-uniform grid correction - val = A(x, y) * g11_2D(x, y) * d1_dx_2D(x, y) / (2. * dx); - xp += val; - xm -= val; - - // g11*dAdx*dfdx - val = g11_2D(x, y) * (A(x + 1, y) - A(x - 1, y)) / (4. * SQ(dx)); - xp += val; - xm -= val; - - // B*f - c += B(x, y); - - // Put values into the preconditioner, X derivatives only - acoef(y - localmesh->ystart, x - xstart) = xm; - bcoef(y - localmesh->ystart, x - xstart) = c; - ccoef(y - localmesh->ystart, x - xstart) = xp; - - if (include_y_derivs) { - BoutReal dy = dy_2D(x, y); - BoutReal dAdx = (A(x + 1, y) - A(x - 1, y)) / (2. * dx); - BoutReal dAdy = (A(x, y + 1) - A(x, y - 1)) / (2. * dy); - - // A*(G2-1/J*d/dy(J/g_22))*dfdy - val = A(x, y) * coef_dfdy(x, y) / (2. * dy); - yp = val; - ym = -val; - - // A*(g22-1/g_22)*d2fdy2 - val = A(x, y) * (g22_2D(x, y) - 1. / g_22_2D(x, y)) / SQ(dy); - yp += val; - c -= 2. * val; - ym += val; - // Non-uniform mesh correction - val = A(x, y) * (g22_2D(x, y) - 1. / g_22_2D(x, y)) * d1_dy_2D(x, y) / (2. * dy); - yp += val; - ym -= val; - - // 2*A*g12*d2dfdxdy - val = A(x, y) * g12_2D(x, y) / (2. * dx * dy); - xpyp = val; - xpym = -val; - xmyp = -val; - xmym = val; - - // g22*dAdy*dfdy - val = (g22_2D(x, y) - 1. / g_22_2D(x, y)) * dAdy / (2. * dy); - yp += val; - ym -= val; - - // g12*(dAdx*dfdy + dAdy*dfdx) - val = g12_2D(x, y) * dAdx / (2. * dy); - yp += val; - ym -= val; - val = g12_2D(x, y) * dAdy / (2. * dx); - xp += val; - xm -= val; - } - - ///////////////////////////////////////////////// - // Now have a 9-point stencil for the Laplacian - - int row = globalIndex(x, y); - - // Set the centre (diagonal) - MatSetValues(MatA, 1, &row, 1, &row, &c, INSERT_VALUES); - - // X + 1 - int col = globalIndex(x + 1, y); - MatSetValues(MatA, 1, &row, 1, &col, &xp, INSERT_VALUES); - - // X - 1 - col = globalIndex(x - 1, y); - MatSetValues(MatA, 1, &row, 1, &col, &xm, INSERT_VALUES); - - if (include_y_derivs) { - // Y + 1 - col = globalIndex(x, y + 1); - MatSetValues(MatA, 1, &row, 1, &col, &yp, INSERT_VALUES); - - // Y - 1 - col = globalIndex(x, y - 1); - MatSetValues(MatA, 1, &row, 1, &col, &ym, INSERT_VALUES); - - // X + 1, Y + 1 - col = globalIndex(x + 1, y + 1); - MatSetValues(MatA, 1, &row, 1, &col, &xpyp, INSERT_VALUES); - - // X + 1, Y - 1 - col = globalIndex(x + 1, y - 1); - MatSetValues(MatA, 1, &row, 1, &col, &xpym, INSERT_VALUES); - - // X - 1, Y + 1 - col = globalIndex(x - 1, y + 1); - MatSetValues(MatA, 1, &row, 1, &col, &xmyp, INSERT_VALUES); - - // X - 1, Y - 1 - col = globalIndex(x - 1, y - 1); - MatSetValues(MatA, 1, &row, 1, &col, &xmym, INSERT_VALUES); - } - } - } -} - -LaplaceXY::~LaplaceXY() { - PetscBool is_finalised; - PetscFinalized(&is_finalised); - - if (!is_finalised) { - // PetscFinalize may already have destroyed this object - KSPDestroy(&ksp); - } - - VecDestroy(&xs); - VecDestroy(&bs); - MatDestroy(&MatA); -} - -const Field2D LaplaceXY::solve(const Field2D& rhs, const Field2D& x0) { - Timer timer("invert"); - - ASSERT1(rhs.getMesh() == localmesh); - ASSERT1(x0.getMesh() == localmesh); - ASSERT1(rhs.getLocation() == location); - ASSERT1(x0.getLocation() == location); - - // Load initial guess x0 into xs and rhs into bs - - for (int x = localmesh->xstart; x <= localmesh->xend; x++) { - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - int ind = globalIndex(x, y); - - PetscScalar val = x0(x, y); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - val = rhs(x, y); - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - } - - if (finite_volume) { - solveFiniteVolume(x0); - } else { - solveFiniteDifference(x0); - } - - // Assemble RHS Vector - VecAssemblyBegin(bs); - VecAssemblyEnd(bs); - - // Assemble Trial Solution Vector - VecAssemblyBegin(xs); - VecAssemblyEnd(xs); - - // Solve the system - KSPSolve(ksp, bs, xs); - - KSPConvergedReason reason; - KSPGetConvergedReason(ksp, &reason); - - if (reason <= 0) { - throw BoutException("LaplaceXY failed to converge. Reason {} ({:d})", - KSPConvergedReasons[reason], static_cast(reason)); - } - - if (save_performance) { - // Update performance monitoring information - n_calls++; - - int iterations = 0; - KSPGetIterationNumber(ksp, &iterations); - - average_iterations = BoutReal(n_calls - 1) / BoutReal(n_calls) * average_iterations - + BoutReal(iterations) / BoutReal(n_calls); - } - - ////////////////////////// - // Copy data into result - - Field2D result; - result.allocate(); - result.setLocation(location); - - for (int x = localmesh->xstart; x <= localmesh->xend; x++) { - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - int ind = globalIndex(x, y); - - PetscScalar val; - VecGetValues(xs, 1, &ind, &val); - result(x, y) = val; - } - } - - // Inner X boundary - if (localmesh->firstX()) { - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - int ind = globalIndex(localmesh->xstart - 1, y); - PetscScalar val; - VecGetValues(xs, 1, &ind, &val); - for (int x = localmesh->xstart - 1; x >= 0; x--) { - result(x, y) = val; - } - } - } - - // Outer X boundary - if (localmesh->lastX()) { - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - int ind = globalIndex(localmesh->xend + 1, y); - PetscScalar val; - VecGetValues(xs, 1, &ind, &val); - for (int x = localmesh->xend + 1; x < localmesh->LocalNx; x++) { - result(x, y) = val; - } - } - } - - // Lower Y boundary - for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { - if ( - // Should not go into corner cells, finite-volume LaplaceXY stencil does not include - // them - (finite_volume and (it.ind < localmesh->xstart or it.ind > localmesh->xend)) - - // Only go into first corner cell for finite-difference - or (it.ind < localmesh->xstart - 1 or it.ind > localmesh->xend + 1)) { - continue; - } - int ind = globalIndex(it.ind, localmesh->ystart - 1); - PetscScalar val; - VecGetValues(xs, 1, &ind, &val); - for (int y = localmesh->ystart - 1; y >= 0; y--) { - result(it.ind, y) = val; - } - } - - // Upper Y boundary - for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { - if ( - // Should not go into corner cells, finite-volume LaplaceXY stencil does not include - // them - (finite_volume and (it.ind < localmesh->xstart or it.ind > localmesh->xend)) - - // Only go into first corner cell for finite-difference - or (it.ind < localmesh->xstart - 1 or it.ind > localmesh->xend + 1)) { - continue; - } - int ind = globalIndex(it.ind, localmesh->yend + 1); - PetscScalar val; - VecGetValues(xs, 1, &ind, &val); - for (int y = localmesh->yend + 1; y < localmesh->LocalNy; y++) { - result(it.ind, y) = val; - } - } - - return result; -} - -void LaplaceXY::solveFiniteVolume(const Field2D& x0) { - // Use original LaplaceXY implementation of passing boundary values for backward - // compatibility - if (localmesh->firstX()) { - if (x_inner_dirichlet) { - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - int ind = globalIndex(localmesh->xstart - 1, y); - - PetscScalar val = x0(localmesh->xstart - 1, y); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - val = 0.5 * (x0(localmesh->xstart - 1, y) + x0(localmesh->xstart, y)); - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - } else { - // Inner X boundary (Neumann) - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - int ind = globalIndex(localmesh->xstart - 1, y); - - PetscScalar val = x0(localmesh->xstart - 1, y); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - val = 0.0; // x0(localmesh->xstart-1,y) - x0(localmesh->xstart,y); - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - } - } - - // Outer X boundary (Dirichlet) - if (localmesh->lastX()) { - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - int ind = globalIndex(localmesh->xend + 1, y); - - PetscScalar val = x0(localmesh->xend + 1, y); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - val = 0.5 * (x0(localmesh->xend, y) + x0(localmesh->xend + 1, y)); - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - } - - if (y_bndry == "dirichlet") { - for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { - // Should not go into corner cells, LaplaceXY stencil does not include them - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - int ind = globalIndex(it.ind, localmesh->ystart - 1); - - PetscScalar val = x0(it.ind, localmesh->ystart - 1); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - val = 0.5 * (x0(it.ind, localmesh->ystart - 1) + x0(it.ind, localmesh->ystart)); - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - - for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { - // Should not go into corner cells, LaplaceXY stencil does not include them - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - int ind = globalIndex(it.ind, localmesh->yend + 1); - - PetscScalar val = x0(it.ind, localmesh->yend + 1); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - val = 0.5 * (x0(it.ind, localmesh->yend + 1) + x0(it.ind, localmesh->yend)); - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - } else if (y_bndry == "neumann" or y_bndry == "free_o3") { - // Y boundaries Neumann - for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { - // Should not go into corner cells, LaplaceXY stencil does not include them - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - int ind = globalIndex(it.ind, localmesh->ystart - 1); - - PetscScalar val = x0(it.ind, localmesh->ystart - 1); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - val = 0.0; - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - - for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { - // Should not go into corner cells, LaplaceXY stencil does not include them - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - int ind = globalIndex(it.ind, localmesh->yend + 1); - - PetscScalar val = x0(it.ind, localmesh->yend + 1); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - val = 0.0; - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - } else { - throw BoutException("Unsupported option for y_bndry"); - } -} - -void LaplaceXY::solveFiniteDifference(const Field2D& x0) { - // For finite-difference implementation pass boundary values in the same way as for - // the 'Laplacian' solvers - the value to use (for Dirichlet boundary conditions) on - // the boundary (which is half way between grid cell and boundary cell) is passed as - // the value in the first boundary cell. - if (localmesh->firstX()) { - if (x_inner_dirichlet) { - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - int ind = globalIndex(localmesh->xstart - 1, y); - - // For the boundary value of the initial guess, use the value that would be set by - // applying the boundary condition to the initial guess - PetscScalar val = 2. * x0(localmesh->xstart - 1, y) - x0(localmesh->xstart, y); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - // Pass the value from boundary cell of x0 as the boundary condition to the rhs - val = x0(localmesh->xstart - 1, y); - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - } else { - // Inner X boundary (Neumann) - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - int ind = globalIndex(localmesh->xstart - 1, y); - - // For the boundary value of the initial guess, use the value that would be set by - // applying the boundary condition to the initial guess - PetscScalar val = x0(localmesh->xstart, y); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - val = 0.0; - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - } - } - - // Outer X boundary (Dirichlet) - if (localmesh->lastX()) { - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - int ind = globalIndex(localmesh->xend + 1, y); - - // For the boundary value of the initial guess, use the value that would be set by - // applying the boundary condition to the initial guess - PetscScalar val = 2. * x0(localmesh->xend + 1, y) - x0(localmesh->xend, y); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - // Pass the value from boundary cell of x0 as the boundary condition to the rhs - val = x0(localmesh->xend + 1, y); - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - } - - if (y_bndry == "dirichlet") { - for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { - // Should not go into corner cells, they are treated specially below - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - int ind = globalIndex(it.ind, localmesh->ystart - 1); - - // For the boundary value of the initial guess, use the value that would be set by - // applying the boundary condition to the initial guess - PetscScalar val = - 2. * x0(it.ind, localmesh->ystart - 1) - x0(it.ind, localmesh->ystart); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - // Pass the value from boundary cell of x0 as the boundary condition to the rhs - val = x0(it.ind, localmesh->ystart - 1); - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - - for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { - // Should not go into corner cells, they are treated specially below - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - int ind = globalIndex(it.ind, localmesh->yend + 1); - - // For the boundary value of the initial guess, use the value that would be set by - // applying the boundary condition to the initial guess - PetscScalar val = - 2. * x0(it.ind, localmesh->yend + 1) - x0(it.ind, localmesh->yend); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - // Pass the value from boundary cell of x0 as the boundary condition to the rhs - val = x0(it.ind, localmesh->yend + 1); - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - - // Use free_o3 for the corner boundary cells - if (localmesh->hasBndryLowerY()) { - if (localmesh->firstX()) { - int ind = globalIndex(localmesh->xstart - 1, localmesh->ystart - 1); - - // For the boundary value of the initial guess, use the value that would be set by - // applying the boundary condition to the initial guess - PetscScalar val = 3. * x0(localmesh->xstart - 1, localmesh->ystart) - - 3. * x0(localmesh->xstart - 1, localmesh->ystart + 1) - + x0(localmesh->xstart - 1, localmesh->ystart + 2); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - // Set the value for the rhs of the boundary condition to zero - val = 0.0; - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - if (localmesh->lastX()) { - int ind = globalIndex(localmesh->xend + 1, localmesh->ystart - 1); - - // For the boundary value of the initial guess, use the value that would be set by - // applying the boundary condition to the initial guess - PetscScalar val = 3. * x0(localmesh->xend + 1, localmesh->ystart) - - 3. * x0(localmesh->xend + 1, localmesh->ystart + 1) - + x0(localmesh->xend + 1, localmesh->ystart + 2); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - // Set the value for the rhs of the boundary condition to zero - val = 0.0; - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - } - if (localmesh->hasBndryUpperY()) { - if (localmesh->firstX()) { - int ind = globalIndex(localmesh->xstart - 1, localmesh->yend + 1); - - // For the boundary value of the initial guess, use the value that would be set by - // applying the boundary condition to the initial guess - PetscScalar val = 3. * x0(localmesh->xstart - 1, localmesh->yend) - - 3. * x0(localmesh->xstart - 1, localmesh->yend - 1) - + x0(localmesh->xstart - 1, localmesh->yend - 2); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - // Set the value for the rhs of the boundary condition to zero - val = 0.0; - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - if (localmesh->lastX()) { - int ind = globalIndex(localmesh->xend + 1, localmesh->yend + 1); - - // For the boundary value of the initial guess, use the value that would be set by - // applying the boundary condition to the initial guess - PetscScalar val = 3. * x0(localmesh->xend + 1, localmesh->yend) - - 3. * x0(localmesh->xend + 1, localmesh->yend - 1) - + x0(localmesh->xend + 1, localmesh->yend - 2); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - // Set the value for the rhs of the boundary condition to zero - val = 0.0; - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - } - } else if (y_bndry == "neumann") { - // Y boundaries Neumann - for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { - // Should not go into corner cells, they are treated specially below - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - int ind = globalIndex(it.ind, localmesh->ystart - 1); - - // For the boundary value of the initial guess, use the value that would be set by - // applying the boundary condition to the initial guess - PetscScalar val = x0(it.ind, localmesh->ystart); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - // Set the value for the rhs of the boundary condition to zero - val = 0.0; - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - if (localmesh->hasBndryLowerY()) { - if (localmesh->firstX()) { - int ind = globalIndex(localmesh->xstart - 1, localmesh->ystart - 1); - - // For the boundary value of the initial guess, use the value that would be set by - // applying the boundary condition to the initial guess - PetscScalar val = x0(localmesh->xstart - 1, localmesh->ystart); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - // Set the value for the rhs of the boundary condition to zero - val = 0.0; - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - if (localmesh->lastX()) { - int ind = globalIndex(localmesh->xend + 1, localmesh->ystart - 1); - - // For the boundary value of the initial guess, use the value that would be set by - // applying the boundary condition to the initial guess - PetscScalar val = x0(localmesh->xend + 1, localmesh->ystart); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - // Set the value for the rhs of the boundary condition to zero - val = 0.0; - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - } - - for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { - // Should not go into corner cells, they are treated specially below - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - int ind = globalIndex(it.ind, localmesh->yend + 1); - - // For the boundary value of the initial guess, use the value that would be set by - // applying the boundary condition to the initial guess - PetscScalar val = x0(it.ind, localmesh->yend); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - // Set the value for the rhs of the boundary condition to zero - val = 0.0; - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - if (localmesh->hasBndryUpperY()) { - if (localmesh->firstX()) { - int ind = globalIndex(localmesh->xstart - 1, localmesh->yend + 1); - - // For the boundary value of the initial guess, use the value that would be set by - // applying the boundary condition to the initial guess - PetscScalar val = x0(localmesh->xstart - 1, localmesh->yend); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - // Set the value for the rhs of the boundary condition to zero - val = 0.0; - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - if (localmesh->lastX()) { - int ind = globalIndex(localmesh->xend + 1, localmesh->yend + 1); - - // For the boundary value of the initial guess, use the value that would be set by - // applying the boundary condition to the initial guess - PetscScalar val = x0(localmesh->xend + 1, localmesh->yend); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - // Set the value for the rhs of the boundary condition to zero - val = 0.0; - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - } - } else if (y_bndry == "free_o3") { - // Y boundaries free_o3 - for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { - // Should not go into corner cells, they are treated specially below - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - int ind = globalIndex(it.ind, localmesh->ystart - 1); - - // For the boundary value of the initial guess, use the value that would be set by - // applying the boundary condition to the initial guess - PetscScalar val = 3. * x0(it.ind, localmesh->ystart) - - 3. * x0(it.ind, localmesh->ystart + 1) - + x0(it.ind, localmesh->ystart + 2); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - // Set the value for the rhs of the boundary condition to zero - val = 0.0; - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - if (localmesh->hasBndryLowerY()) { - if (localmesh->firstX()) { - int ind = globalIndex(localmesh->xstart - 1, localmesh->ystart - 1); - - // For the boundary value of the initial guess, use the value that would be set by - // applying the boundary condition to the initial guess - PetscScalar val = 3. * x0(localmesh->xstart - 1, localmesh->ystart) - - 3. * x0(localmesh->xstart - 1, localmesh->ystart + 1) - + x0(localmesh->xstart - 1, localmesh->ystart + 2); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - // Set the value for the rhs of the boundary condition to zero - val = 0.0; - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - if (localmesh->lastX()) { - int ind = globalIndex(localmesh->xend + 1, localmesh->ystart - 1); - - // For the boundary value of the initial guess, use the value that would be set by - // applying the boundary condition to the initial guess - PetscScalar val = 3. * x0(localmesh->xend + 1, localmesh->ystart) - - 3. * x0(localmesh->xend + 1, localmesh->ystart + 1) - + x0(localmesh->xend + 1, localmesh->ystart + 2); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - // Set the value for the rhs of the boundary condition to zero - val = 0.0; - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - } - - for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { - // Should not go into corner cells, they are treated specially below - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - int ind = globalIndex(it.ind, localmesh->yend + 1); - - // For the boundary value of the initial guess, use the value that would be set by - // applying the boundary condition to the initial guess - PetscScalar val = 3. * x0(it.ind, localmesh->yend) - - 3. * x0(it.ind, localmesh->yend - 1) - + x0(it.ind, localmesh->yend - 2); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - // Set the value for the rhs of the boundary condition to zero - val = 0.0; - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - if (localmesh->hasBndryUpperY()) { - if (localmesh->firstX()) { - int ind = globalIndex(localmesh->xstart - 1, localmesh->yend + 1); - - // For the boundary value of the initial guess, use the value that would be set by - // applying the boundary condition to the initial guess - PetscScalar val = 3. * x0(localmesh->xstart - 1, localmesh->yend) - - 3. * x0(localmesh->xstart - 1, localmesh->yend - 1) - + x0(localmesh->xstart - 1, localmesh->yend - 2); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - // Set the value for the rhs of the boundary condition to zero - val = 0.0; - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - if (localmesh->lastX()) { - int ind = globalIndex(localmesh->xend + 1, localmesh->yend + 1); - - // For the boundary value of the initial guess, use the value that would be set by - // applying the boundary condition to the initial guess - PetscScalar val = 3. * x0(localmesh->xend + 1, localmesh->yend) - - 3. * x0(localmesh->xend + 1, localmesh->yend - 1) - + x0(localmesh->xend + 1, localmesh->yend - 2); - VecSetValues(xs, 1, &ind, &val, INSERT_VALUES); - - // Set the value for the rhs of the boundary condition to zero - val = 0.0; - VecSetValues(bs, 1, &ind, &val, INSERT_VALUES); - } - } - } else { - throw BoutException("Unsupported option for y_bndry"); - } -} - -/*! Preconditioner - * NOTE: For generality, this routine does use globalIndex() in the inner loop, although - * this may be slightly less efficient than incrementing an integer for the global index, - * the finite-volume and finite-difference implementations have slightly different - * indexing patterns, so incrementing an integer would be tricky. - */ -int LaplaceXY::precon(Vec input, Vec result) { - - for (auto itdwn = localmesh->iterateBndryLowerY(); !itdwn.isDone(); itdwn++) { - // Should not go into corner cells, LaplaceXY stencil does not include them - if (itdwn.ind < localmesh->xstart or itdwn.ind > localmesh->xend) { - continue; - } - const int ind = globalIndex(itdwn.ind, localmesh->ystart - 1); - PetscScalar val; - VecGetValues(input, 1, &ind, &val); - VecSetValues(result, 1, &ind, &val, INSERT_VALUES); - } - for (auto itup = localmesh->iterateBndryUpperY(); !itup.isDone(); itup++) { - // Should not go into corner cells, LaplaceXY stencil does not include them - if (itup.ind < localmesh->xstart or itup.ind > localmesh->xend) { - continue; - } - const int ind = globalIndex(itup.ind, localmesh->yend + 1); - PetscScalar val; - VecGetValues(input, 1, &ind, &val); - VecSetValues(result, 1, &ind, &val, INSERT_VALUES); - } - - // Load vector x into bvals array - for (int x = xstart; x <= xend; x++) { - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - const int ind = globalIndex(x, y); - PetscScalar val; - VecGetValues(input, 1, &ind, &val); - bvals(y - localmesh->ystart, x - xstart) = val; - } - } - - // Solve tridiagonal systems using CR solver - cr->solve(bvals, xvals); - - // Save result xvals into y array - for (int x = xstart; x <= xend; x++) { - for (int y = localmesh->ystart; y <= localmesh->yend; y++) { - const int ind = globalIndex(x, y); - PetscScalar val = xvals(y - localmesh->ystart, x - xstart); - VecSetValues(result, 1, &ind, &val, INSERT_VALUES); - } - } - VecAssemblyBegin(result); - VecAssemblyEnd(result); - return 0; -} - -/////////////////////////////////////////////////////////////// - -int LaplaceXY::localSize() { - - // Bulk of points - const int nx = localmesh->xend - localmesh->xstart + 1; - const int ny = localmesh->yend - localmesh->ystart + 1; - - int n = nx * ny; - - // X boundaries - if (localmesh->firstX()) { - n += ny; - } - if (localmesh->lastX()) { - n += ny; - } - - // Y boundaries - for (RangeIterator it = localmesh->iterateBndryLowerY(); !it.isDone(); it++) { - // Should not go into corner cells, LaplaceXY stencil does not include them - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - n++; - } - if ((not finite_volume) and localmesh->hasBndryLowerY()) { - if (localmesh->firstX()) { - n++; - } - if (localmesh->lastX()) { - n++; - } - } - for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) { - // Should not go into corner cells, LaplaceXY stencil does not include them - if (it.ind < localmesh->xstart or it.ind > localmesh->xend) { - continue; - } - n++; - } - if ((not finite_volume) and localmesh->hasBndryUpperY()) { - if (localmesh->firstX()) { - n++; - } - if (localmesh->lastX()) { - n++; - } - } - - return n; -} - -int LaplaceXY::globalIndex(int x, int y) { - if ((x < 0) || (x >= localmesh->LocalNx) || (y < 0) || (y >= localmesh->LocalNy)) { - return -1; // Out of range - } - - // Get the index from a Field2D, round to integer - return static_cast(std::round(indexXY(x, y))); -} - -void LaplaceXY::savePerformance(Solver& solver, const std::string& name) { - // set flag so that performance monitoring values are calculated - save_performance = true; - - // add values to be saved to the output - if (not name.empty()) { - default_prefix = name; - } - - // add monitor to reset counters/averages for new output timestep - // monitor added to back of queue, so that values are reset after being saved - solver.addMonitor(&monitor, Solver::BACK); -} - -int LaplaceXY::LaplaceXYMonitor::call(Solver* /*solver*/, BoutReal /*time*/, int /*iter*/, - int /*nout*/) { - laplacexy.output_average_iterations = laplacexy.average_iterations; - - laplacexy.n_calls = 0; - laplacexy.average_iterations = 0.; - - return 0; -} - -void LaplaceXY::LaplaceXYMonitor::outputVars(Options& output_options, - const std::string& time_dimension) { - output_options[fmt::format("{}_average_iterations", laplacexy.default_prefix)] - .assignRepeat(laplacexy.output_average_iterations, time_dimension); -} - -#endif // BOUT_HAS_PETSC +// DO NOT REMOVE: ensures linker keeps all symbols in this TU +void LaplaceXYFactory::ensureRegistered() {} diff --git a/tests/integrated/test-laplacexy-fv/test-laplacexy.cxx b/tests/integrated/test-laplacexy-fv/test-laplacexy.cxx index 8408f08c00..be2f993b31 100644 --- a/tests/integrated/test-laplacexy-fv/test-laplacexy.cxx +++ b/tests/integrated/test-laplacexy-fv/test-laplacexy.cxx @@ -23,18 +23,23 @@ * **************************************************************************/ -#include -#include -#include -#include -#include -#include +#include "bout/bout.hxx" +#include "bout/bout_types.hxx" +#include "bout/difops.hxx" +#include "bout/field2d.hxx" +#include "bout/initialprofiles.hxx" +#include "bout/invert/laplacexy.hxx" +#include "bout/options.hxx" +#include "bout/options_io.hxx" +#include "bout/output.hxx" + +#include int main(int argc, char** argv) { BoutInitialise(argc, argv); - LaplaceXY laplacexy; + auto laplacexy = LaplaceXY::create(); // Solving equations of the form // Div(A Grad_perp(f)) + B*f = rhs @@ -57,9 +62,9 @@ int main(int argc, char** argv) { Field2D rhs, rhs_check; rhs = Laplace_perpXY(a, f) + b * f; - laplacexy.setCoefs(a, b); + laplacexy->setCoefs(a, b); - sol = laplacexy.solve(rhs, 0.); + sol = laplacexy->solve(rhs, 0.); error = (f - sol) / f; absolute_error = f - sol; max_error = max(abs(absolute_error), true); diff --git a/tests/integrated/test-laplacexy-short/test-laplacexy.cxx b/tests/integrated/test-laplacexy-short/test-laplacexy.cxx index 3486760810..7c284290b3 100644 --- a/tests/integrated/test-laplacexy-short/test-laplacexy.cxx +++ b/tests/integrated/test-laplacexy-short/test-laplacexy.cxx @@ -24,31 +24,39 @@ **************************************************************************/ #include -#include +#include #include +#include +#include +#include #include #include #include +#include +#include +#include + +#include int main(int argc, char** argv) { BoutInitialise(argc, argv); using bout::globals::mesh; - auto coords = mesh->getCoordinates(); + auto* coords = mesh->getCoordinates(); auto& opt = Options::root(); - LaplaceXY laplacexy; + auto laplacexy = LaplaceXY::create(); bool include_y_derivs = opt["laplacexy"]["include_y_derivs"]; // 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); @@ -61,7 +69,8 @@ int main(int argc, char** argv) { //////////////////////////////////////////////////////////////////////////////////////// - Field2D rhs, rhs_check; + Field2D rhs; + Field2D rhs_check; if (include_y_derivs) { rhs = a * DC(Laplace_perp(f)) + DC(Grad_perp(a) * Grad_perp(f)) + b * f; } else { @@ -69,14 +78,14 @@ int main(int argc, char** argv) { a * DC(Delp2(f, CELL_DEFAULT, false)) + DC(coords->g11 * DDX(a) * DDX(f)) + b * f; } - laplacexy.setCoefs(a, b); + laplacexy->setCoefs(a, b); - sol = laplacexy.solve(rhs, 0.); - error = (f - sol) / f; - absolute_error = f - sol; - max_error = max(abs(absolute_error), true); + auto sol = laplacexy->solve(rhs, 0.); + const auto error = (f - sol) / f; + const auto absolute_error = f - sol; + const auto 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); if (include_y_derivs) { diff --git a/tests/integrated/test-laplacexy/test-laplacexy.cxx b/tests/integrated/test-laplacexy/test-laplacexy.cxx index 3142d359b1..543274d37d 100644 --- a/tests/integrated/test-laplacexy/test-laplacexy.cxx +++ b/tests/integrated/test-laplacexy/test-laplacexy.cxx @@ -24,11 +24,19 @@ **************************************************************************/ #include -#include +#include #include +#include +#include +#include #include #include #include +#include +#include +#include + +#include using bout::globals::mesh; @@ -64,10 +72,10 @@ int main(int argc, char** argv) { rhs = a * Delp2(f, CELL_DEFAULT, false) + coords->g11 * DDX(a) * DDX(f) + b * f; } - LaplaceXY laplacexy; - laplacexy.setCoefs(a, b); + auto laplacexy = LaplaceXY::create(); + laplacexy->setCoefs(a, b); - Field2D solution = laplacexy.solve(rhs, 0.); + 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); diff --git a/tests/integrated/test-laplacexy2-hypre/test-laplacexy.cxx b/tests/integrated/test-laplacexy2-hypre/test-laplacexy.cxx index 2fb8a10777..e04de9c675 100644 --- a/tests/integrated/test-laplacexy2-hypre/test-laplacexy.cxx +++ b/tests/integrated/test-laplacexy2-hypre/test-laplacexy.cxx @@ -24,17 +24,22 @@ **************************************************************************/ #include -#include -#include +#include +#include +#include #include -#include +#include #include +#include +#include + +#include int main(int argc, char** argv) { BoutInitialise(argc, argv); - LaplaceXY2Hypre laplacexy; + auto laplacexy = LaplaceXYFactory::getInstance().create("hypre", nullptr); // Solving equations of the form // Div(A Grad_perp(f)) + B*f = rhs @@ -54,10 +59,10 @@ int main(int argc, char** argv) { Field2D rhs = Laplace_perpXY(a, f) + b * f; - laplacexy.setCoefs(a, b); + laplacexy->setCoefs(a, b); Field2D guess = 0.0; - Field2D sol = laplacexy.solve(rhs, guess); + Field2D sol = laplacexy->solve(rhs, guess); Field2D error = (f - sol) / f; // Absolute value of relative error: abs((f - sol)/f) Field2D absolute_error = abs(f - sol); diff --git a/tools/pylib/boutconfig/__init__.py.cin b/tools/pylib/boutconfig/__init__.py.cin index d031d599b5..b65b25a75a 100644 --- a/tools/pylib/boutconfig/__init__.py.cin +++ b/tools/pylib/boutconfig/__init__.py.cin @@ -31,6 +31,7 @@ config = { "has_fftw": "@BOUT_HAS_FFTW@", "petsc_has_sundials": "@PETSC_HAS_SUNDIALS@", "metric_type": "@BOUT_METRIC_TYPE@", + "has_hypre": "@BOUT_HAS_HYPRE@", } for k, v in config.items():