From eb92f732e9f41e86a1f0843981aa7dcb9702e939 Mon Sep 17 00:00:00 2001 From: Peter Hill Date: Mon, 9 Nov 2026 15:43:05 +0000 Subject: [PATCH] Remove geometric multigrid Laplacian inversion solver Unmaintained and doesn't seem to be used in practice --- CHANGELOG.md | 1 + CMakeLists.txt | 4 - include/bout/invert_laplace.hxx | 1 - manual/sphinx/user_docs/laplacian.rst | 19 +- .../laplace/impls/multigrid/multigrid_alg.cxx | 873 ------------------ .../impls/multigrid/multigrid_laplace.cxx | 715 -------------- .../impls/multigrid/multigrid_laplace.hxx | 249 ----- .../impls/multigrid/multigrid_solver.cxx | 771 ---------------- src/invert/laplace/invert_laplace.cxx | 1 - tests/integrated/CMakeLists.txt | 1 - .../test-multigrid_laplace/CMakeLists.txt | 9 - .../test-multigrid_laplace/data/BOUT.inp | 121 --- .../data/BOUT_jy127.inp | 123 --- .../test-multigrid_laplace/data/BOUT_jy4.inp | 123 --- .../test-multigrid_laplace/data/BOUT_jy63.inp | 123 --- .../data/BOUT_unsheared.inp | 114 --- .../test-multigrid_laplace/plotresults.py | 53 -- .../test_multigrid_laplace.cxx | 326 ------- .../test_multigrid_laplace.py | 51 - .../test_multigrid_laplace_unsheared.py | 52 -- 20 files changed, 5 insertions(+), 3725 deletions(-) delete mode 100644 src/invert/laplace/impls/multigrid/multigrid_alg.cxx delete mode 100644 src/invert/laplace/impls/multigrid/multigrid_laplace.cxx delete mode 100644 src/invert/laplace/impls/multigrid/multigrid_laplace.hxx delete mode 100644 src/invert/laplace/impls/multigrid/multigrid_solver.cxx delete mode 100644 tests/integrated/test-multigrid_laplace/CMakeLists.txt delete mode 100644 tests/integrated/test-multigrid_laplace/data/BOUT.inp delete mode 100644 tests/integrated/test-multigrid_laplace/data/BOUT_jy127.inp delete mode 100644 tests/integrated/test-multigrid_laplace/data/BOUT_jy4.inp delete mode 100644 tests/integrated/test-multigrid_laplace/data/BOUT_jy63.inp delete mode 100644 tests/integrated/test-multigrid_laplace/data/BOUT_unsheared.inp delete mode 100755 tests/integrated/test-multigrid_laplace/plotresults.py delete mode 100644 tests/integrated/test-multigrid_laplace/test_multigrid_laplace.cxx delete mode 100755 tests/integrated/test-multigrid_laplace/test_multigrid_laplace.py delete mode 100755 tests/integrated/test-multigrid_laplace/test_multigrid_laplace_unsheared.py diff --git a/CHANGELOG.md b/CHANGELOG.md index a0cbfd5915..b557d3d61d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,7 @@ ```console $ bout-upgrader v6 v6_upgrader file/to/update.cxx ```` +- Removed `multigrid` Laplacian inversion solver ## [v5.2.0](https://github.com/boutproject/BOUT-dev/tree/v5.2.0 diff --git a/CMakeLists.txt b/CMakeLists.txt index 90fe822aa3..028a98f109 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -248,10 +248,6 @@ set(BOUT_SOURCES ./src/invert/laplace/impls/hypre3d/hypre3d_laplace.hxx ./src/invert/laplace/impls/iterative_parallel_tri/iterative_parallel_tri.cxx ./src/invert/laplace/impls/iterative_parallel_tri/iterative_parallel_tri.hxx - ./src/invert/laplace/impls/multigrid/multigrid_alg.cxx - ./src/invert/laplace/impls/multigrid/multigrid_laplace.cxx - ./src/invert/laplace/impls/multigrid/multigrid_laplace.hxx - ./src/invert/laplace/impls/multigrid/multigrid_solver.cxx ./src/invert/laplace/impls/naulin/naulin_laplace.cxx ./src/invert/laplace/impls/naulin/naulin_laplace.hxx ./src/invert/laplace/impls/pcr/pcr.cxx diff --git a/include/bout/invert_laplace.hxx b/include/bout/invert_laplace.hxx index 76212dcd76..b7f7f916b0 100644 --- a/include/bout/invert_laplace.hxx +++ b/include/bout/invert_laplace.hxx @@ -66,7 +66,6 @@ constexpr auto LAPLACE_PETSCAMG = "petscamg"; constexpr auto LAPLACE_PETSC3DAMG = "petsc3damg"; constexpr auto LAPLACE_HYPRE3D = "hypre3d"; constexpr auto LAPLACE_CYCLIC = "cyclic"; -constexpr auto LAPLACE_MULTIGRID = "multigrid"; constexpr auto LAPLACE_NAULIN = "naulin"; constexpr auto LAPLACE_IPT = "ipt"; constexpr auto LAPLACE_PCR = "pcr"; diff --git a/manual/sphinx/user_docs/laplacian.rst b/manual/sphinx/user_docs/laplacian.rst index 175abfe082..631fb03fde 100644 --- a/manual/sphinx/user_docs/laplacian.rst +++ b/manual/sphinx/user_docs/laplacian.rst @@ -41,8 +41,6 @@ implementations are listed in table :numref:`tab-laplacetypes`. | petsc3damg | Serial/parallel. Solves full 3D operator (with | PETSc (section :ref:`sec-PETSc-install`) | | | y-derivatives) with algebraic multigrid. | | +------------------------+--------------------------------------------------------------+------------------------------------------+ - | multigrid | Serial/parallel. Geometric multigrid, no Boussinesq | | - +------------------------+--------------------------------------------------------------+------------------------------------------+ | `naulin | Serial/parallel. Iterative treatment of non-Boussinesq terms | | | `__ | | | +------------------------+--------------------------------------------------------------+------------------------------------------+ @@ -87,11 +85,10 @@ and :math:`n_z` are the number of grid-points in the :math:`x` and :math:`z` directions respectively. -In the second approach, the full :math:`2`\ -D system is solved. The -available solvers for this approach are 'multigrid' using a multigrid -algorithm; 'naulin' using an iterative scheme to correct the FFT-based -approach; or 'petsc' using KSP linear solvers from the PETSc library -(this requires PETSc to be built with BOUT++). +In the second approach, the full :math:`2`\ -D system is solved. The available +solvers for this approach are ``naulin`` using an iterative scheme to correct +the FFT-based approach; or ``petsc`` using KSP linear solvers from the PETSc +library (this requires PETSc to be built with BOUT++). The `Laplacian` class is defined in ``invert_laplace.hxx`` and solves @@ -780,14 +777,6 @@ Cyclic algorithm This is now the default solver in both serial and parallel. It is an FFT-based solver using a cyclic reduction algorithm. -.. _sec-multigrid: - -Multigrid solver -~~~~~~~~~~~~~~~~ - -A solver using a geometric multigrid algorithm was introduced by projects in -2015 and 2016 of CCFE and the EUROfusion HLST. - .. _sec-naulin: Naulin solver diff --git a/src/invert/laplace/impls/multigrid/multigrid_alg.cxx b/src/invert/laplace/impls/multigrid/multigrid_alg.cxx deleted file mode 100644 index f741f53d55..0000000000 --- a/src/invert/laplace/impls/multigrid/multigrid_alg.cxx +++ /dev/null @@ -1,873 +0,0 @@ -/************************************************************************** - * Perpendicular Laplacian inversion. - * Using Geometrical Multigrid Solver - * - * Equation solved is: - d*\nabla^2_\perp x + (1/c1)\nabla_perp c2\cdot\nabla_\perp x + a x = b - * - ************************************************************************** - * Copyright 2016 K.S. Kang - * - * 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 . - * - **************************************************************************/ - -#include "bout/build_defines.hxx" - -#if not BOUT_USE_METRIC_3D - -#include "multigrid_laplace.hxx" -#include "bout/unused.hxx" -#include - -// Define basic multigrid algorithm - -MultigridAlg::MultigridAlg(int level, int lx, int lz, int gx, int gz, MPI_Comm comm, - int check) - : mglevel(level), pcheck(check), commMG(comm) { - - if (pcheck > 0) { - output << "Construct MG " << level << endl; - } - - // Memory allocate for Multigrid - gnx.reallocate(mglevel); - gnz.reallocate(mglevel); - lnx.reallocate(mglevel); - lnz.reallocate(mglevel); - - gnx[mglevel - 1] = gx; - gnz[mglevel - 1] = gz; - lnx[mglevel - 1] = lx; - lnz[mglevel - 1] = lz; - if (mglevel > 1) { - for (int i = mglevel - 1; i > 0; i--) { - gnx[i - 1] = gnx[i] / 2; - gnz[i - 1] = gnz[i] / 2; - lnx[i - 1] = lnx[i] / 2; - lnz[i - 1] = lnz[i] / 2; - } - } - - // Could be replaced with a Matrix - matmg = new BoutReal*[mglevel]; - for (int i = 0; i < mglevel; i++) { - matmg[i] = new BoutReal[(lnx[i] + 2) * (lnz[i] + 2) * 9]; - } -} - -MultigridAlg::~MultigridAlg() { - output << "End deconstruction Malg AAAA " << numP << endl; - for (int i = 0; i < mglevel; i++) { - delete[] matmg[i]; - } - delete[] matmg; -} - -void MultigridAlg::getSolution(BoutReal* x, BoutReal* b, int flag) { - - // swap ghost cells of initial guess - communications(x, mglevel - 1); - // don't think the ghost cells of the rhs are used, so don't need to be communicated - // /JTO 17-4-2019 - - if (flag == 0) { - //Solve exaclty - if (mglevel == 1) { - pGMRES(x, b, mglevel - 1, 1); - } else if (mgplag == 1) { - pGMRES(x, b, mglevel - 1, 1); - } else { - solveMG(x, b, mglevel - 1); - } - } else { - cycleMG(mglevel - 1, x, b); - if (flag > 1) { - int level = mglevel - 1; - int ldim = (lnx[level] + 2) * (lnz[level] + 2); - Array y(ldim); - Array r(ldim); - for (int n = 1; n < flag; n++) { - residualVec(level, x, b, std::begin(r)); - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int i = 0; i < ldim; i++) { - y[i] = 0.0; - } - cycleMG(level, std::begin(y), std::begin(r)); - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int i = 0; i < ldim; i++) { - x[i] = x[i] + y[i]; - } - } - } - } -} - -void MultigridAlg::cycleMG(int level, BoutReal* sol, BoutReal* rhs) { - if (level == 0) { - lowestSolver(sol, rhs, 0); - } else { - Array r((lnx[level] + 2) * (lnz[level] + 2)); - Array pr((lnx[level - 1] + 2) * (lnz[level - 1] + 2)); - Array y((lnx[level - 1] + 2) * (lnz[level - 1] + 2)); - Array iy((lnx[level] + 2) * (lnz[level] + 2)); - - smoothings(level, sol, rhs); - - residualVec(level, sol, rhs, std::begin(r)); - - projection(level, std::begin(r), std::begin(pr)); - - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int i = 0; i < (lnx[level - 1] + 2) * (lnz[level - 1] + 2); i++) { - y[i] = 0.0; - } - - cycleMG(level - 1, std::begin(y), std::begin(pr)); - - prolongation(level - 1, std::begin(y), std::begin(iy)); - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int i = 0; i < (lnx[level] + 2) * (lnz[level] + 2); i++) { - sol[i] += iy[i]; - } - - smoothings(level, sol, rhs); - } -} - -void MultigridAlg::projection(int level, BoutReal* r, BoutReal* pr) { - - BOUT_OMP_PERF(parallel default(shared)) - { - BOUT_OMP_PERF(for) - for (int i = 0; i < (lnx[level - 1] + 2) * (lnz[level - 1] + 2); i++) { - pr[i] = 0.; - } - int xend = lnx[level - 1] + 1; - int zend = lnz[level - 1] + 1; - BOUT_OMP_PERF(for collapse(2)) - for (int i = 1; i < xend; i++) { - for (int k = 1; k < zend; k++) { - int i2 = 2 * i - 1; - int k2 = 2 * k - 1; - int nn = i * (lnz[level - 1] + 2) + k; - int n0 = i2 * (lnz[level] + 2) + k2; - int n1 = n0 + 1; - int n2 = n0 + lnz[level] + 2; - int n3 = n2 + 1; - pr[nn] = (r[n0] + r[n1] + r[n2] + r[n3]) / 4.0; - } - } - } - communications(pr, level - 1); -} - -void MultigridAlg::prolongation(int level, BoutReal* x, BoutReal* ix) { - - BOUT_OMP_PERF(parallel default(shared)) - { - BOUT_OMP_PERF(for) - for (int i = 0; i < (lnx[level + 1] + 2) * (lnz[level + 1] + 2); i++) { - ix[i] = 0.; - } - - int xend = lnx[level] + 1; - int zend = lnz[level] + 1; - BOUT_OMP_PERF(for collapse(2)) - for (int i = 1; i < xend; i++) { - for (int k = 1; k < zend; k++) { - int i2 = 2 * i - 1; - int k2 = 2 * k - 1; - int nn = i * (lnz[level] + 2) + k; - int n0 = i2 * (lnz[level + 1] + 2) + k2; - int n1 = n0 + 1; - int n2 = n0 + lnz[level + 1] + 2; - int n3 = n2 + 1; - ix[n0] = x[nn]; - ix[n1] = x[nn]; - ix[n2] = x[nn]; - ix[n3] = x[nn]; - } - } - } - communications(ix, level + 1); -} - -void MultigridAlg::smoothings(int level, BoutReal* x, BoutReal* b) { - - int dim; - int mm = lnz[level] + 2; - dim = mm * (lnx[level] + 2); - if (mgsm == 0) { - Array x0(dim); - BOUT_OMP_PERF(parallel default(shared)) - for (int num = 0; num < 2; num++) { - BOUT_OMP_PERF(for) - for (int i = 0; i < dim; i++) { - x0[i] = x[i]; - } - - int xend = lnx[level] + 1; - int zend = lnz[level] + 1; - BOUT_OMP_PERF(for collapse(2)) - for (int i = 1; i < xend; i++) { - for (int k = 1; k < zend; k++) { - int nn = i * mm + k; - BoutReal val = b[nn] - matmg[level][nn * 9 + 3] * x0[nn - 1] - - matmg[level][nn * 9 + 5] * x0[nn + 1] - - matmg[level][nn * 9 + 1] * x0[nn - mm] - - matmg[level][nn * 9 + 7] * x0[nn + mm] - - matmg[level][nn * 9] * x0[nn - mm - 1] - - matmg[level][nn * 9 + 2] * x0[nn - mm + 1] - - matmg[level][nn * 9 + 6] * x0[nn + mm - 1] - - matmg[level][nn * 9 + 8] * x0[nn + mm + 1]; - if (fabs(matmg[level][nn * 9 + 4]) < atol) { - throw BoutException("Error at matmg({:d}-{:d})", level, nn); - } - - x[nn] = (1.0 - omega) * x[nn] + omega * val / matmg[level][nn * 9 + 4]; - } - } - communications(x, level); - } - } else { - for (int i = 1; i < lnx[level] + 1; i++) { - for (int k = 1; k < lnz[level] + 1; k++) { - int nn = i * mm + k; - BoutReal val = b[nn] - matmg[level][nn * 9 + 3] * x[nn - 1] - - matmg[level][nn * 9 + 5] * x[nn + 1] - - matmg[level][nn * 9 + 1] * x[nn - mm] - - matmg[level][nn * 9 + 7] * x[nn + mm] - - matmg[level][nn * 9] * x[nn - mm - 1] - - matmg[level][nn * 9 + 2] * x[nn - mm + 1] - - matmg[level][nn * 9 + 6] * x[nn + mm - 1] - - matmg[level][nn * 9 + 8] * x[nn + mm + 1]; - if (fabs(matmg[level][nn * 9 + 4]) < atol) { - throw BoutException("Error at matmg({:d}-{:d})", level, nn); - } - x[nn] = val / matmg[level][nn * 9 + 4]; - } - } - communications(x, level); - for (int i = lnx[level]; i > 0; i--) { - for (int k = lnz[level]; k > 0; k--) { - int nn = i * mm + k; - BoutReal val = b[nn] - matmg[level][nn * 9 + 3] * x[nn - 1] - - matmg[level][nn * 9 + 5] * x[nn + 1] - - matmg[level][nn * 9 + 1] * x[nn - mm] - - matmg[level][nn * 9 + 7] * x[nn + mm] - - matmg[level][nn * 9] * x[nn - mm - 1] - - matmg[level][nn * 9 + 2] * x[nn - mm + 1] - - matmg[level][nn * 9 + 6] * x[nn + mm - 1] - - matmg[level][nn * 9 + 8] * x[nn + mm + 1]; - if (fabs(matmg[level][nn * 9 + 4]) < atol) { - throw BoutException("Error at matmg({:d}-{:d})", level, nn); - } - x[nn] = val / matmg[level][nn * 9 + 4]; - } - } - communications(x, level); - } -} - -void MultigridAlg::pGMRES(BoutReal* sol, BoutReal* rhs, int level, int iplag) { - int it, etest = 1, MAXIT; - BoutReal ini_e, error, a0, a1, rederr, perror; - BoutReal** v; - BoutReal c[MAXGM + 1], s[MAXGM + 1], y[MAXGM + 1], g[MAXGM + 1], - h[MAXGM + 1][MAXGM + 1]; - - if ((level == 0) || (iplag == 0)) { - MAXIT = 40000; - } else { - MAXIT = 500; - } - - int ldim = (lnx[level] + 2) * (lnz[level] + 2); - // Could we use a Matrix here? - v = new BoutReal*[MAXGM + 1]; - for (int i = 0; i < MAXGM + 1; i++) { - v[i] = new BoutReal[ldim]; - } - - Array p(ldim); - Array q(ldim); - Array r(ldim); - - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int i = 0; i < ldim; i++) { - sol[i] = 0.0; - } - int num = 0; - - ini_e = sqrt(vectorProd(level, rhs, rhs)); - if ((pcheck == 1) && (rProcI == 0)) { - output << numP << "--In GMRES ini " << ini_e << endl; - } - if (ini_e < atol * rtol) { - if ((pcheck == 1) && (rProcI == 0)) { - output << numP << "Don't need to solve. E= " << ini_e << endl; - } - // Clean up memory before returning from method - for (int i = 0; i < MAXGM + 1; i++) { - delete[] v[i]; - } - delete[] v; - return; - } - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int i = 0; i < ldim; i++) { - r[i] = 0.0; - } - if (iplag == 0) { - smoothings(level, std::begin(r), rhs); - } else { - cycleMG(level, std::begin(r), rhs); - } - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int i = 0; i < ldim; i++) { - v[0][i] = r[i]; - } - perror = ini_e; - do { - a1 = vectorProd(level, v[0], v[0]); - a1 = sqrt(a1); - if (fabs(a1) < atol * rtol) { - output << num << " First a1 in GMRES is wrong at level " << level << ": " << a1 - << endl; - } - a0 = 1.0 / a1; - g[0] = a1; - BOUT_OMP_PERF(parallel default(shared)) - { - BOUT_OMP_PERF(for) - for (int i = 0; i < ldim; i++) { - v[0][i] *= a0; - } - BOUT_OMP_PERF(for) - for (int i = 1; i < MAXGM + 1; i++) { - g[i] = 0.0; - } - } - for (it = 0; it < MAXGM; it++) { - multiAVec(level, v[it], std::begin(q)); - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int i = 0; i < ldim; i++) { - v[it + 1][i] = 0.0; - } - - if (iplag == 0) { - smoothings(level, v[it + 1], std::begin(q)); - } else { - cycleMG(level, v[it + 1], std::begin(q)); - } - - for (int i = 0; i < it + 1; i++) { - h[i][it] = vectorProd(level, v[it + 1], v[i]); - } - for (int i = 0; i < it + 1; i++) { - a0 = -h[i][it]; - for (int k = 0; k < ldim; k++) { - v[it + 1][k] += a0 * v[i][k]; - } - } - a1 = vectorProd(level, v[it + 1], v[it + 1]); - a1 = sqrt(a1); - - // if ldim==9 then there is only one grid point at this level, so the - // solution will be exact, the residual will vanish and we will exit this - // loop on the first iteration, so the value of a0=1/a1=infinity will - // never be used. Therefore this check is not needed in that case - if (fabs(a1) < atol * rtol && ldim > 9) { - output << num << " Second a1 in GMRES is wrong at level " << level << ": " << a1 - << endl; - } - a0 = 1.0 / a1; - h[it + 1][it] = a1; - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int i = 0; i < ldim; i++) { - v[it + 1][i] *= a0; - } - - for (int i = 0; i < it; i++) { - a0 = c[i] * h[i][it] - s[i] * h[i + 1][it]; - a1 = s[i] * h[i][it] + c[i] * h[i + 1][it]; - h[i][it] = a0; - h[i + 1][it] = a1; - } - a0 = h[it][it] * h[it][it] + h[it + 1][it] * h[it + 1][it]; - if (a0 < 0.0) { - throw BoutException("a0 in GMRES is negative \n"); - } - a0 = sqrt(a0); - if (fabs(a0) < atol * rtol) { - throw BoutException("a0 in GMRES is wrong \n"); - } - c[it] = h[it][it] / a0; - s[it] = -h[it + 1][it] / a0; - h[it][it] = a0; - h[it + 1][it] = 0.0; - a0 = c[it] * g[it] - s[it] * g[it + 1]; - a1 = s[it] * g[it] + c[it] * g[it + 1]; - g[it] = a0; - g[it + 1] = a1; - - /* Get solution y and x_m*/ - for (int i = it; i >= 0; i--) { - y[i] = g[i]; - for (int j = i + 1; j <= it; j++) { - y[i] -= h[i][j] * y[j]; - } - y[i] = y[i] / h[i][i]; - } - BOUT_OMP_PERF(parallel default(shared)) - { - BOUT_OMP_PERF(for) - for (int i = 0; i < ldim; i++) { - p[i] = sol[i]; - } - BOUT_OMP_PERF(for) - for (int k = 0; k < ldim; k++) { - for (int i = 0; i <= it; i++) { - p[k] += y[i] * v[i][k]; - } - } - } - - /* Get r_m and test convergence.*/ - residualVec(level, std::begin(p), rhs, std::begin(r)); - error = sqrt(vectorProd(level, std::begin(r), std::begin(r))); - num += 1; - if (error > dtol) { - throw BoutException("GMRES reached dtol with error {:16.10f} at iteration {:d}\n", - error, num); - } - if (num > MAXIT) { - throw BoutException( - "GMRES reached MAXIT with error {:16.10f} at iteration {:d}\n", error, num); - } - if (error <= rtol * ini_e + atol) { - etest = 0; - break; - } - // J. Omotani, 27/2/2018: I think this test is intended to check for slow - // convergence of the GMRES solve, and 'abort' if it is converging - // slowly. This is OK on a coarse level solver, because at worst it means - // the top-level iteration will have to continue but the top level - // iteration should only be stopped by the previous test against the - // tolerance. - // Therefore, check that this is not the top-level solver before applying - // this test. - if ((level < mglevel - 1) && (fabs(perror - error) / error < rtol)) { - if (it == 0) { - etest = 0; - } - num -= 1; - break; - } - perror = error; - } - /* Restart with new initial */ - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int i = 0; i < ldim; i++) { - v[0][i] = 0.0; - } - if (iplag == 0) { - smoothings(level, v[0], std::begin(r)); - } else { - cycleMG(level, v[0], std::begin(r)); - } - - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int i = 0; i < ldim; i++) { - sol[i] = p[i]; - } - if (num > MAXIT) { - throw BoutException(" GMRES Iteration limit.\n"); - } - // if((etest == 1) & (xProcI == 0)) - // printf("Restart GMRES %d | %20.14f\n",num,error/ini_e); - } while (etest == 1); - - if (level == 0) { - if ((rProcI == 0) && (pcheck == 1)) { - rederr = log(error / ini_e) / static_cast(num); - rederr = exp(rederr); - printf("The average error reduction of GMRES %d: %14.8f(%18.10f)\n", num, rederr, - error); - fflush(stdout); - } - } else { - if ((rProcI == 0) && (pcheck == 1)) { - rederr = log(error / ini_e) / static_cast(num); - rederr = exp(rederr); - printf("The average error reduction of PGMRES %d: %14.8f(%18.10f)\n", num, rederr, - error); - fflush(stdout); - } - } - for (int i = 0; i < MAXGM + 1; i++) { - delete[] v[i]; - } - delete[] v; -} - -void MultigridAlg::setMultigridC(int UNUSED(plag)) { - - int level = mglevel - 1; - for (int n = level; n > 0; n--) { - setMatrixC(n); - if (pcheck == 2) { - output << n << "matrix in Alg = " << lnx[n] << "," << lnz[n] << endl; - output << gnx[n] << "," << gnz[n] << endl; - } - } -} - -void MultigridAlg::lowestSolver(BoutReal* x, BoutReal* b, int UNUSED(plag)) { - pGMRES(x, b, 0, 0); -} - -BoutReal MultigridAlg::vectorProd(int level, BoutReal* x, BoutReal* y) { - // nx does not include guard cells - - BoutReal val; - BoutReal ini_e = 0.0; - BOUT_OMP_PERF(parallel default(shared)) - { - int xend = lnx[level] + 1; - int zend = lnz[level] + 1; - BOUT_OMP_PERF(for reduction(+:ini_e) collapse(2)) - for (int i = 1; i < xend; i++) { - for (int k = 1; k < zend; k++) { - int ii = i * (lnz[level] + 2) + k; - ini_e += x[ii] * y[ii]; - } - } - } - if (numP > 1) { - bout::globals::mpi->MPI_Allreduce(&ini_e, &val, 1, MPI_DOUBLE, MPI_SUM, commMG); - } else { - val = ini_e; - } - - return (val); -} - -void MultigridAlg::multiAVec(int level, BoutReal* x, BoutReal* b) { - - int mm = lnz[level] + 2; - BOUT_OMP_PERF(parallel default(shared)) - { - BOUT_OMP_PERF(for) - for (int i = 0; i < mm * (lnx[level] + 2); i++) { - b[i] = 0.0; - } - - int xend = lnx[level] + 1; - int zend = lnz[level] + 1; - BOUT_OMP_PERF(for collapse(2)) - for (int i = 1; i < xend; i++) { - for (int k = 1; k < zend; k++) { - int nn = i * mm + k; - b[nn] = matmg[level][nn * 9 + 4] * x[nn] + matmg[level][nn * 9 + 3] * x[nn - 1] - + matmg[level][nn * 9 + 5] * x[nn + 1] - + matmg[level][nn * 9 + 1] * x[nn - mm] - + matmg[level][nn * 9 + 7] * x[nn + mm] - + matmg[level][nn * 9] * x[nn - mm - 1] - + matmg[level][nn * 9 + 2] * x[nn - mm + 1] - + matmg[level][nn * 9 + 6] * x[nn + mm - 1] - + matmg[level][nn * 9 + 8] * x[nn + mm + 1]; - } - } - } - communications(b, level); -} - -void MultigridAlg::residualVec(int level, BoutReal* x, BoutReal* b, BoutReal* r) { - - int mm; - mm = lnz[level] + 2; - BOUT_OMP_PERF(parallel default(shared)) - { - BOUT_OMP_PERF(for) - for (int i = 0; i < mm * (lnx[level] + 2); i++) { - r[i] = 0.0; - } - - int xend = lnx[level] + 1; - int zend = lnz[level] + 1; - BOUT_OMP_PERF(for collapse(2)) - for (int i = 1; i < xend; i++) { - for (int k = 1; k < zend; k++) { - int nn = i * mm + k; - BoutReal val = - matmg[level][nn * 9 + 4] * x[nn] + matmg[level][nn * 9 + 3] * x[nn - 1] - + matmg[level][nn * 9 + 5] * x[nn + 1] + matmg[level][nn * 9 + 1] * x[nn - mm] - + matmg[level][nn * 9 + 7] * x[nn + mm] - + matmg[level][nn * 9] * x[nn - mm - 1] - + matmg[level][nn * 9 + 2] * x[nn - mm + 1] - + matmg[level][nn * 9 + 6] * x[nn + mm - 1] - + matmg[level][nn * 9 + 8] * x[nn + mm + 1]; - r[nn] = b[nn] - val; - } - } - } - communications(r, level); -} - -void MultigridAlg::setMatrixC(int level) { - - BoutReal ratio = 8.0; - - BOUT_OMP_PERF(parallel default(shared)) - { - BOUT_OMP_PERF(for) - for (int i = 0; i < (lnx[level - 1] + 2) * (lnz[level - 1] + 2) * 9; i++) { - matmg[level - 1][i] = 0.0; - } - - int xend = lnx[level - 1] + 1; - int zend = lnz[level - 1] + 1; - BOUT_OMP_PERF(for collapse(2)) - for (int i = 1; i < xend; i++) { - for (int k = 1; k < zend; k++) { - int i2 = 2 * i - 1; - int k2 = 2 * k - 1; - int mm = i * (lnz[level - 1] + 2) + k; - int m0 = i2 * (lnz[level] + 2) + k2; - int m1 = i2 * (lnz[level] + 2) + k2 + 1; - int m2 = (i2 + 1) * (lnz[level] + 2) + k2; - int m3 = (i2 + 1) * (lnz[level] + 2) + k2 + 1; - BoutReal val = matmg[level][m0 * 9 + 4] + matmg[level][m1 * 9 + 4]; - val += matmg[level][m2 * 9 + 4] + matmg[level][m3 * 9 + 4]; - val += matmg[level][m0 * 9 + 5] + matmg[level][m1 * 9 + 3]; - val += matmg[level][m2 * 9 + 5] + matmg[level][m3 * 9 + 3]; - val += matmg[level][m0 * 9 + 7] + matmg[level][m2 * 9 + 1]; - val += matmg[level][m1 * 9 + 7] + matmg[level][m3 * 9 + 1]; - val += matmg[level][m0 * 9 + 8] + matmg[level][m3 * 9]; - val += matmg[level][m1 * 9 + 6] + matmg[level][m2 * 9 + 2]; - matmg[level - 1][mm * 9 + 4] = val / ratio; - val = matmg[level][m0 * 9 + 1] + matmg[level][m1 * 9 + 1]; - val += matmg[level][m0 * 9 + 2] + matmg[level][m1 * 9]; - matmg[level - 1][mm * 9 + 1] = val / ratio; - val = matmg[level][m0 * 9 + 3] + matmg[level][m2 * 9 + 3]; - val += matmg[level][m0 * 9 + 6] + matmg[level][m2 * 9]; - matmg[level - 1][mm * 9 + 3] = val / ratio; - val = matmg[level][m1 * 9 + 5] + matmg[level][m3 * 9 + 5]; - val += matmg[level][m1 * 9 + 8] + matmg[level][m3 * 9 + 2]; - matmg[level - 1][mm * 9 + 5] = val / ratio; - val = matmg[level][m2 * 9 + 7] + matmg[level][m3 * 9 + 7]; - val += matmg[level][m2 * 9 + 8] + matmg[level][m3 * 9 + 6]; - matmg[level - 1][mm * 9 + 7] = val / ratio; - matmg[level - 1][mm * 9] = matmg[level][m0 * 9] / ratio; - matmg[level - 1][mm * 9 + 2] = matmg[level][m1 * 9 + 2] / ratio; - matmg[level - 1][mm * 9 + 6] = matmg[level][m2 * 9 + 6] / ratio; - matmg[level - 1][mm * 9 + 8] = matmg[level][m3 * 9 + 8] / ratio; - } - } - } -} - -void MultigridAlg::communications(BoutReal* x, int level) { - // Note: currently z-direction communications and x-direction communications are done - // independently. They should really be done together if zNP>1 and xNP>1 (with a single - // MPI_Waitall after all the MPI_Isend and MPI_Irecv calls, instead of the two - // MPI_Waitalls there are now. As there are never any z-communications at the moment, it - // is not worth implementing for now. - - MPI_Status status[4]; - int stag, rtag; - [[maybe_unused]] int ierr; - - if (zNP > 1) { - MPI_Request requests[] = {MPI_REQUEST_NULL, MPI_REQUEST_NULL, MPI_REQUEST_NULL, - MPI_REQUEST_NULL}; - MPI_Datatype xvector; - - // Create datatype to hold z-direction guard cells, which are not continuous in memory - ierr = MPI_Type_vector(lnx[level], 1, lnz[level] + 2, MPI_DOUBLE, &xvector); - ASSERT1(ierr == MPI_SUCCESS); - - ierr = MPI_Type_commit(&xvector); - ASSERT1(ierr == MPI_SUCCESS); - - // Receive from z- - rtag = zProcM; - ierr = bout::globals::mpi->MPI_Irecv(&x[lnz[level] + 2], 1, xvector, zProcM, rtag, - commMG, &requests[2]); - ASSERT1(ierr == MPI_SUCCESS); - - // Receive from z+ - rtag = zProcP + numP; - ierr = bout::globals::mpi->MPI_Irecv(&x[2 * (lnz[level] + 2) - 1], 1, xvector, zProcP, - rtag, commMG, &requests[3]); - ASSERT1(ierr == MPI_SUCCESS); - - // Send to z+ - stag = rProcI; - ierr = bout::globals::mpi->MPI_Isend(&x[2 * (lnz[level] + 2) - 2], 1, xvector, zProcP, - stag, commMG, &requests[0]); - ASSERT1(ierr == MPI_SUCCESS); - - // Send to z- - stag = rProcI + numP; - ierr = bout::globals::mpi->MPI_Isend(&x[lnz[level] + 3], 1, xvector, zProcM, stag, - commMG, &requests[1]); - ASSERT1(ierr == MPI_SUCCESS); - - // Wait for communications to complete - ierr = bout::globals::mpi->MPI_Waitall(4, requests, status); - ASSERT1(ierr == MPI_SUCCESS); - - ierr = MPI_Type_free(&xvector); - ASSERT1(ierr == MPI_SUCCESS); - } else { - for (int i = 1; i < lnx[level] + 1; i++) { - x[i * (lnz[level] + 2)] = x[(i + 1) * (lnz[level] + 2) - 2]; - x[(i + 1) * (lnz[level] + 2) - 1] = x[i * (lnz[level] + 2) + 1]; - } - } - if (xNP > 1) { - // Note: periodic x-direction not handled here - - MPI_Request requests[] = {MPI_REQUEST_NULL, MPI_REQUEST_NULL, MPI_REQUEST_NULL, - MPI_REQUEST_NULL}; - - if (xProcI > 0) { - // Receive from x- - rtag = xProcM; - ierr = bout::globals::mpi->MPI_Irecv(&x[0], lnz[level] + 2, MPI_DOUBLE, xProcM, - rtag, commMG, &requests[2]); - ASSERT1(ierr == MPI_SUCCESS); - } - - if (xProcI < xNP - 1) { - // Receive from x+ - rtag = xProcP + xNP; - ; - ierr = bout::globals::mpi->MPI_Irecv(&x[(lnx[level] + 1) * (lnz[level] + 2)], - lnz[level] + 2, MPI_DOUBLE, xProcP, rtag, - commMG, &requests[3]); - ASSERT1(ierr == MPI_SUCCESS); - - // Send to x+ - stag = rProcI; - ierr = - bout::globals::mpi->MPI_Isend(&x[lnx[level] * (lnz[level] + 2)], lnz[level] + 2, - MPI_DOUBLE, xProcP, stag, commMG, &requests[0]); - ASSERT1(ierr == MPI_SUCCESS); - } - - if (xProcI > 0) { - // Send to x- - stag = rProcI + xNP; - ierr = bout::globals::mpi->MPI_Isend(&x[lnz[level] + 2], lnz[level] + 2, MPI_DOUBLE, - xProcM, stag, commMG, &requests[1]); - ASSERT1(ierr == MPI_SUCCESS); - } - - // Wait for communications to complete - ierr = bout::globals::mpi->MPI_Waitall(4, requests, status); - ASSERT1(ierr == MPI_SUCCESS); - } else { - for (int i = 0; i < lnz[level] + 2; i++) { - x[i] = x[lnx[level] * (lnz[level] + 2) + i]; - x[(lnx[level] + 1) * (lnz[level] + 2) + i] = x[(lnz[level] + 2) + i]; - } - } -} - -void MultigridAlg::solveMG(BoutReal* sol, BoutReal* rhs, int level) { - int m, MAXIT = 150; - BoutReal ini_e, perror, error, rederr; - int ldim = (lnx[level] + 2) * (lnz[level] + 2); - - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int i = 0; i < ldim; i++) { - sol[i] = 0.0; - } - - ini_e = vectorProd(level, rhs, rhs); - if (ini_e < 0.0) { - throw BoutException("In MG Initial Error {:10.4e} \n", ini_e); - } - ini_e = sqrt(ini_e); - if ((pcheck == 1) && (rProcI == 0)) { - printf("%d \n In MGsolve ini = %24.18f\n", numP, ini_e); - } - Array y(ldim); - Array r(ldim); - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int i = 0; i < ldim; i++) { - r[i] = rhs[i]; - } - - perror = ini_e; - for (m = 0; m < MAXIT; m++) { - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int i = 0; i < ldim; i++) { - y[i] = 0.0; - } - cycleMG(level, std::begin(y), std::begin(r)); - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int i = 0; i < ldim; i++) { - sol[i] = sol[i] + y[i]; - } - residualVec(level, sol, rhs, std::begin(r)); - error = sqrt(vectorProd(level, std::begin(r), std::begin(r))); - if ((pcheck == 1) && (rProcI == 0)) { - printf("%d \n In MGsolve error = %24.18f\n", m, error); - } - if (error < rtol * ini_e + atol) { - break; - } - if ((fabs(perror - error) / error < rtol) || (error > dtol)) { - throw BoutException("In MG Limited Error {:10.4e} \n", error); - } - perror = error; - } - - if ((rProcI == 0) && (pcheck == 1)) { - rederr = log(error / ini_e) / (static_cast(m) + 1.0); - rederr = exp(rederr); - if (m == MAXIT) { - printf("Reached maximum iteration: %14.8f\n", error); - } - printf("The average error reduction of MG %d: %14.8f(%18.12f)\n", m + 1, rederr, - error); - fflush(stdout); - } -} - -#endif diff --git a/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx b/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx deleted file mode 100644 index c06a83585b..0000000000 --- a/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx +++ /dev/null @@ -1,715 +0,0 @@ -/************************************************************************** - * Perpendicular Laplacian inversion. - * Using Geometrical Multigrid Solver - * - * Equation solved is: - * d*\nabla^2_\perp x + (1/c1)\nabla_perp c2\cdot\nabla_\perp x + a x = b - * - ************************************************************************** - * Copyright 2016 K.S. Kang - * - * 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 . - * - **************************************************************************/ - -#include "multigrid_laplace.hxx" -#include "bout/build_defines.hxx" - -#if not BOUT_USE_METRIC_3D - -#include -#include -#include -#include -#include - -#if BOUT_USE_OPENMP -#include -#endif - -BoutReal soltime = 0.0, settime = 0.0; - -LaplaceMultigrid::LaplaceMultigrid(Options* opt, const CELL_LOC loc, Mesh* mesh_in, - Solver* UNUSED(solver)) - : Laplacian(opt, loc, mesh_in), A(0.0), C1(1.0), C2(1.0), D(1.0) { - - // periodic x-direction not handled: see MultigridAlg::communications - ASSERT1(!localmesh->periodicX); - - A.setLocation(location); - C1.setLocation(location); - C2.setLocation(location); - D.setLocation(location); - - // Get Options in Laplace Section - if (!opt) { - opts = Options::getRoot()->getSection("laplace"); - } else { - opts = opt; - } - opts->get("multigridlevel", mglevel, 100, true); - opts->get("rtol", rtol, pow(10.0, -8), true); - opts->get("atol", atol, pow(10.0, -20), true); - opts->get("dtol", dtol, pow(10.0, 5), true); - opts->get("smtype", mgsm, 1, true); - if (mgsm != 0 && omp_get_max_threads() > 1) { - output_warn << "WARNING: in multigrid Laplace solver, for smtype!=0 the smoothing " - "cannot be parallelised with OpenMP threads." - << endl - << " Consider using smtype=0 instead when using OpenMP threads." - << endl; - } - opts->get("jacomega", omega, 0.8, true); - opts->get("solvertype", mgplag, 1, true); - opts->get("cftype", cftype, 0, true); - opts->get("mergempi", mgmpi, 63, true); - opts->get("checking", pcheck, 0, true); - mgcount = 0; - - // Initialize, allocate memory, etc. - comms_tagbase = 385; // Some random number - - constexpr int implemented_global_flags = INVERT_START_NEW; - if (isGlobalFlagSet(~implemented_global_flags)) { - throw BoutException("Attempted to set Laplacian inversion flag that is not " - "implemented in LaplaceMultigrid."); - } - // INVERT_DC_GRAD does not actually do anything, but harmless to set while comparing to Fourier solver with Neumann boundary conditions - constexpr int implemented_boundary_flags = INVERT_AC_GRAD + INVERT_SET + INVERT_DC_GRAD; - if (isInnerBoundaryFlagSet(~implemented_boundary_flags)) { - throw BoutException("Attempted to set Laplacian inner boundary inversion flag that " - "is not implemented in LaplaceMultigrid."); - } - if (isOuterBoundaryFlagSet(~implemented_boundary_flags)) { - throw BoutException("Attempted to set Laplacian outer boundary inversion flag that " - "is not implemented in LaplaceMultigrid."); - } - - commX = localmesh->getXcomm(); - - Nx_local = localmesh->xend - localmesh->xstart + 1; // excluding guard cells - Nx_global = localmesh->GlobalNx - 2 * localmesh->xstart; // excluding guard cells - - if (mgcount == 0) { - output << "Nx=" << Nx_global << "(" << Nx_local << ")" << endl; - } - Nz_global = localmesh->GlobalNz; - Nz_local = Nz_global; // No parallelization in z-direction (for now) - - if (mgcount == 0) { - output << "Nz=" << Nz_global << "(" << Nz_local << ")" << endl; - } - - // Compute available levels along x-direction - if (mglevel > 1) { - int nn = Nx_global; - for (int n = mglevel; n > 1; n--) { - if (nn % 2 != 0) { - output << "Size of global x-domain is not a multiple of 2^" << mglevel - 1 - << " mglevel is changed to " << mglevel - n + 1 << endl; - mglevel = mglevel - n + 1; - break; - } - nn = nn / 2; - } - // ... and check the same for z-direction - nn = Nz_global; - for (int n = mglevel; n > 1; n--) { - if (nn % 2 != 0) { - output << "Size of global z-domain is not a multiple of 2^ " << mglevel - 1 - << " mglevel is changed to " << mglevel - n + 1 << endl; - mglevel = mglevel - n + 1; - break; - } - nn = nn / 2; - } - } else { - mglevel = 1; - } - - // Compute available levels on each processor along x-direction - // aclevel is the number of levels that can be used in parallel, i.e. set by - // the grid size on a single processor - // If the number of levels is higher than aclevel, then the grid is collected - // to a single processor, and a new multigrid solver (called sMG) is created - // to run in serial to compute the coarsest (mglevel-aclevel) levels - int aclevel, adlevel; - if (mglevel > 1) { - int nn = Nx_local; - aclevel = mglevel; - for (int n = aclevel; n > 1; n--) { - if (nn % 2 != 0) { - output << "Size of local x-domain is not a multiple of 2^" << aclevel - << " aclevel is changed to " << aclevel - n + 1 << endl; - aclevel = aclevel - n + 1; - break; - } - nn = nn / 2; - } - // ... and check the same for z-direction - nn = Nz_local; - for (int n = aclevel; n > 1; n--) { - if (nn % 2 != 0) { - output << "Size of local z-domain is not a multiple of 2^ " << aclevel - << " aclevel is changed to " << aclevel - n << endl; - aclevel = aclevel - n + 1; - break; - } - nn = nn / 2; - } - } else { - aclevel = 1; - } - adlevel = mglevel - aclevel; - - kMG = bout::utils::make_unique(aclevel, Nx_local, Nz_local, Nx_global, - adlevel, mgmpi, commX, pcheck); - kMG->mgplag = mgplag; - kMG->mgsm = mgsm; - kMG->cftype = cftype; - kMG->rtol = rtol; - kMG->atol = atol; - kMG->dtol = dtol; - kMG->omega = omega; - kMG->setValueS(); - - // Set up Multigrid Cycle - - x.reallocate((Nx_local + 2) * (Nz_local + 2)); - b.reallocate((Nx_local + 2) * (Nz_local + 2)); - - if (mgcount == 0) { - output << " Smoothing type is "; - if (mgsm == 0) { - output << "Jacobi smoother"; - output << "with omega = " << omega << endl; - } else if (mgsm == 1) { - output << " Gauss-Seidel smoother" << endl; - } else { - throw BoutException("Undefined smoother"); - } - output << "Solver type is "; - if (mglevel == 1) { - output << "PGMRES with simple Preconditioner" << endl; - } else if (mgplag == 1) { - output << "PGMRES with multigrid Preconditioner" << endl; - } else { - output << "Multigrid solver with merging " << mgmpi << endl; - } - BOUT_OMP_SAFE(parallel) - BOUT_OMP_SAFE(master) - { - output << "Num threads = " << omp_get_num_threads() << endl; - } - } -} - -FieldPerp LaplaceMultigrid::solve(const FieldPerp& b_in, const FieldPerp& x0) { - - ASSERT1(localmesh == b_in.getMesh() && localmesh == x0.getMesh()); - ASSERT1(b_in.getLocation() == location); - ASSERT1(x0.getLocation() == location); - - checkData(b_in); - checkData(x0); - ASSERT3(b_in.getIndex() == x0.getIndex()); - - BoutReal t0, t1; - - yindex = b_in.getIndex(); - int level = kMG->mglevel - 1; - int lzz = kMG->lnz[level]; - int lz2 = lzz + 2; - int lxx = kMG->lnx[level]; - - if (isGlobalFlagSet(INVERT_START_NEW)) { - // set initial guess to zero - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for collapse(2)) - for (int i = 1; i < lxx + 1; i++) { - for (int k = 1; k < lzz + 1; k++) { - x[i * lz2 + k] = 0.; - } - } - } else { - // Read initial guess into local array, ignoring guard cells - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for collapse(2)) - for (int i = 1; i < lxx + 1; i++) { - for (int k = 1; k < lzz + 1; k++) { - int i2 = i - 1 + localmesh->xstart; - int k2 = k - 1; - x[i * lz2 + k] = x0[i2][k2]; - } - } - } - - // Read RHS into local array - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for collapse(2)) - for (int i = 1; i < lxx + 1; i++) { - for (int k = 1; k < lzz + 1; k++) { - int i2 = i - 1 + localmesh->xstart; - int k2 = k - 1; - b[i * lz2 + k] = b_in(i2, k2); - } - } - - if (localmesh->firstX()) { - if (isInnerBoundaryFlagSet(INVERT_AC_GRAD)) { - // Neumann boundary condition - if (isInnerBoundaryFlagSet(INVERT_SET)) { - // guard cells of x0 specify gradient to set at inner boundary - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int k = 1; k < lzz + 1; k++) { - int k2 = k - 1; - x[k] = -x0(localmesh->xstart - 1, k2) - * sqrt(coords->g_11(localmesh->xstart, yindex)) - * coords->dx(localmesh->xstart, yindex); - } - } else { - // zero gradient inner boundary condition - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int k = 1; k < lzz + 1; k++) { - // set inner guard cells - x[k] = 0.0; - } - } - } else { - // Dirichlet boundary condition - if (isInnerBoundaryFlagSet(INVERT_SET)) { - // guard cells of x0 specify value to set at inner boundary - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int k = 1; k < lzz + 1; k++) { - int k2 = k - 1; - x[k] = 2. * x0(localmesh->xstart - 1, k2); - // this is the value to set at the inner boundary - } - } else { - // zero value inner boundary condition - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int k = 1; k < lzz + 1; k++) { - // set inner guard cells - x[k] = 0.; - } - } - } - } - if (localmesh->lastX()) { - if (isOuterBoundaryFlagSet(INVERT_AC_GRAD)) { - // Neumann boundary condition - if (isInnerBoundaryFlagSet(INVERT_SET)) { - // guard cells of x0 specify gradient to set at outer boundary - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int k = 1; k < lzz + 1; k++) { - int k2 = k - 1; - x[(lxx + 1) * lz2 + k] = x0(localmesh->xend + 1, k2) - * sqrt(coords->g_11(localmesh->xend, yindex)) - * coords->dx(localmesh->xend, yindex); - // this is the value to set the gradient to at the outer boundary - } - } else { - // zero gradient outer boundary condition - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int k = 1; k < lzz + 1; k++) { - // set outer guard cells - x[(lxx + 1) * lz2 + k] = 0.; - } - } - } else { - // Dirichlet boundary condition - if (isOuterBoundaryFlagSet(INVERT_SET)) { - // guard cells of x0 specify value to set at outer boundary - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int k = 1; k < lzz + 1; k++) { - int k2 = k - 1; - x[(lxx + 1) * lz2 + k] = 2. * x0(localmesh->xend + 1, k2); - // this is the value to set at the outer boundary - } - } else { - // zero value inner boundary condition - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int k = 1; k < lzz + 1; k++) { - // set outer guard cells - x[(lxx + 1) * lz2 + k] = 0.; - } - } - } - } - - // Exchange ghost cells of initial guess - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int i = 0; i < lxx + 2; i++) { - x[i * lz2] = x[(i + 1) * lz2 - 2]; - x[(i + 1) * lz2 - 1] = x[i * lz2 + 1]; - } - - if (mgcount == 0) { - soltime = 0.0; - settime = 0.0; - } else { - if (kMG->pcheck > 0) { - kMG->setPcheck(0); - } - } - - t0 = bout::globals::mpi->MPI_Wtime(); - generateMatrixF(level); - - if (kMG->xNP > 1) { - bout::globals::mpi->MPI_Barrier(commX); - } - - if ((pcheck == 3) && (mgcount == 0)) { - std::string outfile = fmt::format("test_matF_{:d}.mat", kMG->rProcI); - output << "Out file= " << outfile << endl; - FILE* outf = fopen(outfile.c_str(), "w"); - int dim = (lxx + 2) * (lzz + 2); - fprintf(outf, "dim = %d (%d, %d)\n", dim, lxx, lzz); - - for (int i = 0; i < dim; i++) { - fprintf(outf, "%d ==", i); - for (int j = 0; j < 9; j++) { - fprintf(outf, "%12.6f,", kMG->matmg[level][i * 9 + j]); - } - fprintf(outf, "\n"); - } - fclose(outf); - } - - if (level > 0) { - kMG->setMultigridC(0); - } - - if ((pcheck == 3) && (mgcount == 0)) { - for (int i = level; i > 0; i--) { - output << i << "dimension= " << kMG->lnx[i - 1] << "(" << kMG->gnx[i - 1] << ")," - << kMG->lnz[i - 1] << endl; - - std::string outfile = fmt::format("test_matC{:1d}_{:d}.mat", i, kMG->rProcI); - output << "Out file= " << outfile << endl; - FILE* outf = fopen(outfile.c_str(), "w"); - int dim = (kMG->lnx[i - 1] + 2) * (kMG->lnz[i - 1] + 2); - fprintf(outf, "dim = %d (%d,%d)\n", dim, kMG->lnx[i - 1], kMG->lnz[i - 1]); - - for (int ii = 0; ii < dim; ii++) { - fprintf(outf, "%d ==", ii); - for (int j = 0; j < 9; j++) { - fprintf(outf, "%12.6f,", kMG->matmg[i - 1][ii * 9 + j]); - } - fprintf(outf, "\n"); - } - fclose(outf); - } - } - - t1 = bout::globals::mpi->MPI_Wtime(); - settime += t1 - t0; - - // Compute solution. - - mgcount++; - if (pcheck > 0) { - t0 = bout::globals::mpi->MPI_Wtime(); - } - - kMG->getSolution(std::begin(x), std::begin(b), 0); - - if (pcheck > 0) { - t1 = bout::globals::mpi->MPI_Wtime(); - soltime += t1 - t0; - if (mgcount % 300 == 0) { - output << "Accumulated execution time at " << mgcount << " Sol " << soltime << " ( " - << settime << " )" << endl; - settime = 0.; - soltime = 0.; - } - } - - FieldPerp result{emptyFrom(b_in)}; - -#if CHECK > 2 - // Make any unused elements NaN so that user does not try to do calculations with them - BOUT_FOR(i, result.getRegion("RGN_ALL")) { result[i] = BoutNaN; } -#endif - - // Copy solution into a FieldPerp to return - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for collapse(2)) - for (int i = 1; i < lxx + 1; i++) { - for (int k = 1; k < lzz + 1; k++) { - int i2 = i - 1 + localmesh->xstart; - int k2 = k - 1; - result(i2, k2) = x[i * lz2 + k]; - } - } - if (localmesh->firstX()) { - if (isInnerBoundaryFlagSet(INVERT_AC_GRAD)) { - // Neumann boundary condition - if (isInnerBoundaryFlagSet(INVERT_SET)) { - // guard cells of x0 specify gradient to set at inner boundary - int i2 = -1 + localmesh->xstart; - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int k = 1; k < lzz + 1; k++) { - int k2 = k - 1; - result(i2, k2) = x[lz2 + k] - - x0(localmesh->xstart - 1, k2) - * sqrt(coords->g_11(localmesh->xstart, yindex)) - * coords->dx(localmesh->xstart, yindex); - } - } else { - // zero gradient inner boundary condition - int i2 = -1 + localmesh->xstart; - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int k = 1; k < lzz + 1; k++) { - int k2 = k - 1; - result(i2, k2) = x[lz2 + k]; - } - } - } else { - // Dirichlet boundary condition - if (isInnerBoundaryFlagSet(INVERT_SET)) { - // guard cells of x0 specify value to set at inner boundary - int i2 = -1 + localmesh->xstart; - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int k = 1; k < lzz + 1; k++) { - int k2 = k - 1; - result(i2, k2) = 2. * x0(localmesh->xstart - 1, k2) - x[lz2 + k]; - } - } else { - // zero value inner boundary condition - int i2 = -1 + localmesh->xstart; - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int k = 1; k < lzz + 1; k++) { - int k2 = k - 1; - result(i2, k2) = -x[lz2 + k]; - } - } - } - } - if (localmesh->lastX()) { - if (isOuterBoundaryFlagSet(INVERT_AC_GRAD)) { - // Neumann boundary condition - if (isInnerBoundaryFlagSet(INVERT_SET)) { - // guard cells of x0 specify gradient to set at outer boundary - int i2 = lxx + localmesh->xstart; - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int k = 1; k < lzz + 1; k++) { - int k2 = k - 1; - result(i2, k2) = x[lxx * lz2 + k] - + x0(localmesh->xend + 1, k2) - * sqrt(coords->g_11(localmesh->xend, yindex)) - * coords->dx(localmesh->xend, yindex); - } - } else { - // zero gradient outer boundary condition - int i2 = lxx + localmesh->xstart; - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int k = 1; k < lzz + 1; k++) { - int k2 = k - 1; - result(i2, k2) = x[lxx * lz2 + k]; - } - } - } else { - // Dirichlet boundary condition - if (isOuterBoundaryFlagSet(INVERT_SET)) { - // guard cells of x0 specify value to set at outer boundary - int i2 = lxx + localmesh->xstart; - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int k = 1; k < lzz + 1; k++) { - int k2 = k - 1; - result(i2, k2) = 2. * x0(localmesh->xend + 1, k2) - x[lxx * lz2 + k]; - } - } else { - // zero value inner boundary condition - int i2 = lxx + localmesh->xstart; - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int k = 1; k < lzz + 1; k++) { - int k2 = k - 1; - result(i2, k2) = -x[lxx * lz2 + k]; - } - } - } - } - - checkData(result); - - return result; -} - -void LaplaceMultigrid::generateMatrixF(int level) { - - // Set (fine-level) matrix entries - - BoutReal* mat; - mat = kMG->matmg[level]; - int llx = kMG->lnx[level]; - int llz = kMG->lnz[level]; - - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for collapse(2)) - for (int i = 1; i < llx + 1; i++) { - for (int k = 1; k < llz + 1; k++) { - int i2 = i - 1 + localmesh->xstart; - int k2 = k - 1; - int k2p = (k2 + 1) % Nz_global; - int k2m = (k2 + Nz_global - 1) % Nz_global; - - BoutReal dz = coords->dz(i2, yindex); - BoutReal ddx_C = (C2(i2 + 1, yindex, k2) - C2(i2 - 1, yindex, k2)) / 2. - / coords->dx(i2, yindex) / C1(i2, yindex, k2); - BoutReal ddz_C = - (C2(i2, yindex, k2p) - C2(i2, yindex, k2m)) / 2. / dz / C1(i2, yindex, k2); - - BoutReal ddx = D(i2, yindex, k2) * coords->g11(i2, yindex) / coords->dx(i2, yindex) - / coords->dx(i2, yindex); - // coefficient of 2nd derivative stencil (x-direction) - - BoutReal ddz = D(i2, yindex, k2) * coords->g33(i2, yindex) / SQ(dz); - // coefficient of 2nd derivative stencil (z-direction) - - BoutReal dxdz = - D(i2, yindex, k2) * 2. * coords->g13(i2, yindex) / coords->dx(i2, yindex) / dz; - // coefficient of mixed derivative stencil (could assume zero, at least initially, - // if easier; then check this is true in constructor) - - BoutReal dxd = - (D(i2, yindex, k2) * coords->G1(i2, yindex) + coords->g11(i2, yindex) * ddx_C - + coords->g13(i2, yindex) - * ddz_C // (could assume zero, at least initially, if easier; then check this is true in constructor) - ) - / coords->dx(i2, yindex); // coefficient of 1st derivative stencil (x-direction) - if (nonuniform) { - // add correction for non-uniform dx - dxd += D(i2, yindex, k2) * coords->d1_dx(i2, yindex); - } - - BoutReal dzd = - (D(i2, yindex, k2) * coords->G3(i2, yindex) + coords->g33(i2, yindex) * ddz_C - + coords->g13(i2, yindex) - * ddx_C // (could assume zero, at least initially, if easier; then check - // this is true in constructor) - ) - / dz; // coefficient of 1st derivative stencil (z-direction) - - int ic = i * (llz + 2) + k; - mat[ic * 9] = dxdz / 4.; - mat[ic * 9 + 1] = ddx - dxd / 2.; - mat[ic * 9 + 2] = -dxdz / 4.; - mat[ic * 9 + 3] = ddz - dzd / 2.; - mat[ic * 9 + 4] = - A(i2, yindex, k2) - 2. * (ddx + ddz); // coefficient of no-derivative component - mat[ic * 9 + 5] = ddz + dzd / 2.; - mat[ic * 9 + 6] = -dxdz / 4.; - mat[ic * 9 + 7] = ddx + dxd / 2.; - mat[ic * 9 + 8] = dxdz / 4.; - } - } - - // Here put boundary conditions - - if (kMG->rProcI == 0) { - if (isInnerBoundaryFlagSet(INVERT_AC_GRAD)) { - // Neumann boundary condition - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int k = 1; k < llz + 1; k++) { - int ic = llz + 2 + k; - mat[ic * 9 + 3] += mat[ic * 9]; - mat[ic * 9 + 4] += mat[ic * 9 + 1]; - mat[ic * 9 + 5] += mat[ic * 9 + 2]; - b[ic] -= mat[ic * 9] * x[k - 1]; - b[ic] -= mat[ic * 9 + 1] * x[k]; - b[ic] -= mat[ic * 9 + 2] * x[k + 1]; - mat[ic * 9] = 0.; - mat[ic * 9 + 1] = 0.; - mat[ic * 9 + 2] = 0.; - } - } else { - // Dirichlet boundary condition - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int k = 1; k < llz + 1; k++) { - int ic = llz + 2 + k; - mat[ic * 9 + 3] -= mat[ic * 9]; - mat[ic * 9 + 4] -= mat[ic * 9 + 1]; - mat[ic * 9 + 5] -= mat[ic * 9 + 2]; - b[ic] -= mat[ic * 9] * x[k - 1]; - b[ic] -= mat[ic * 9 + 1] * x[k]; - b[ic] -= mat[ic * 9 + 2] * x[k + 1]; - mat[ic * 9] = 0.; - mat[ic * 9 + 1] = 0.; - mat[ic * 9 + 2] = 0.; - } - } - } - if (kMG->rProcI == kMG->xNP - 1) { - if (isOuterBoundaryFlagSet(INVERT_AC_GRAD)) { - // Neumann boundary condition - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int k = 1; k < llz + 1; k++) { - int ic = llx * (llz + 2) + k; - mat[ic * 9 + 3] += mat[ic * 9 + 6]; - mat[ic * 9 + 4] += mat[ic * 9 + 7]; - mat[ic * 9 + 5] += mat[ic * 9 + 8]; - b[ic] -= mat[ic * 9 + 6] * x[(llx + 1) * (llz + 2) + k - 1]; - b[ic] -= mat[ic * 9 + 7] * x[(llx + 1) * (llz + 2) + k]; - b[ic] -= mat[ic * 9 + 8] * x[(llx + 1) * (llz + 2) + k + 1]; - mat[ic * 9 + 6] = 0.; - mat[ic * 9 + 7] = 0.; - mat[ic * 9 + 8] = 0.; - } - } else { - // Dirichlet boundary condition - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int k = 1; k < llz + 1; k++) { - int ic = llx * (llz + 2) + k; - mat[ic * 9 + 3] -= mat[ic * 9 + 6]; - mat[ic * 9 + 4] -= mat[ic * 9 + 7]; - mat[ic * 9 + 5] -= mat[ic * 9 + 8]; - b[ic] -= mat[ic * 9 + 6] * x[(llx + 1) * (llz + 2) + k - 1]; - b[ic] -= mat[ic * 9 + 7] * x[(llx + 1) * (llz + 2) + k]; - b[ic] -= mat[ic * 9 + 8] * x[(llx + 1) * (llz + 2) + k + 1]; - mat[ic * 9 + 6] = 0.; - mat[ic * 9 + 7] = 0.; - mat[ic * 9 + 8] = 0.; - } - } - } -} - -#endif // BOUT_USE_METRIC_3D diff --git a/src/invert/laplace/impls/multigrid/multigrid_laplace.hxx b/src/invert/laplace/impls/multigrid/multigrid_laplace.hxx deleted file mode 100644 index c728902815..0000000000 --- a/src/invert/laplace/impls/multigrid/multigrid_laplace.hxx +++ /dev/null @@ -1,249 +0,0 @@ -/************************************************************************** - * Perpendicular Laplacian inversion. - * Using Geometrical Multigrid Solver - * - * Equation solved is: - d*\nabla^2_\perp x + (1/c1)\nabla_perp c2\cdot\nabla_\perp x + a x = b - * - ************************************************************************** - * Copyright 2015 K.S. Kang - * Modified version Aeptember 2015 - * - * 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_MULTIGRID_LAPLACE_H -#define BOUT_MULTIGRID_LAPLACE_H - -#include "bout/build_defines.hxx" -#include "bout/invert_laplace.hxx" - -#if BOUT_USE_METRIC_3D - -namespace { -RegisterUnavailableLaplace - registerlaplacemultigrid(LAPLACE_MULTIGRID, "BOUT++ was configured with 3D metrics"); -} - -#else - -#include - -#include -#include -#include -#include -#include - -#define MAXGM 15 - -// In multigrid_alg.cxx - -class MultigridAlg { -public: - MultigridAlg(int, int, int, int, int, MPI_Comm, int); - virtual ~MultigridAlg(); - - void setMultigridC(int); - void getSolution(BoutReal*, BoutReal*, int); - - int mglevel, mgplag, cftype, mgsm, pcheck, xNP, zNP, rProcI; - BoutReal rtol, atol, dtol, omega; - Array gnx, gnz, lnx, lnz; - BoutReal** matmg; - -protected: - /******* Start implementation ********/ - int numP, xProcI, zProcI, xProcP, xProcM, zProcP, zProcM; - - MPI_Comm commMG; - - void communications(BoutReal*, int); - void setMatrixC(int); - - void cycleMG(int, BoutReal*, BoutReal*); - void smoothings(int, BoutReal*, BoutReal*); - void projection(int, BoutReal*, BoutReal*); - void prolongation(int, BoutReal*, BoutReal*); - void pGMRES(BoutReal*, BoutReal*, int, int); - void solveMG(BoutReal*, BoutReal*, int); - void multiAVec(int, BoutReal*, BoutReal*); - void residualVec(int, BoutReal*, BoutReal*, BoutReal*); - BoutReal vectorProd(int, BoutReal*, BoutReal*); - - virtual void lowestSolver(BoutReal*, BoutReal*, int); -}; - -// Define three different type of multigrid solver -// in multigrid_solver.cxx - -class MultigridSerial : public MultigridAlg { -public: - MultigridSerial(int level, int gx, int gz, MPI_Comm comm, int check); - ~MultigridSerial(){}; - - void convertMatrixF(BoutReal*); -}; - -class Multigrid2DPf1D : public MultigridAlg { -public: - Multigrid2DPf1D(int, int, int, int, int, int, int, int, MPI_Comm, int); - ~Multigrid2DPf1D(){}; - - void setMultigridC(int); - void setPcheck(int); - void setValueS(); - int kflag; - -private: - std::unique_ptr sMG; - void convertMatrixFS(int); - void lowestSolver(BoutReal*, BoutReal*, int); -}; - -class Multigrid1DP : public MultigridAlg { -public: - Multigrid1DP(int, int, int, int, int, int, MPI_Comm, int); - ~Multigrid1DP(){}; - void setMultigridC(int); - void setPcheck(int); - void setValueS(); - int kflag; - -private: - MPI_Comm comm2D; - std::unique_ptr sMG; - std::unique_ptr rMG; - void convertMatrixF2D(int); - void convertMatrixFS(int); - void lowestSolver(BoutReal*, BoutReal*, int); -}; - -class LaplaceMultigrid : public Laplacian { -public: - LaplaceMultigrid(Options* opt = nullptr, const CELL_LOC loc = CELL_CENTRE, - Mesh* mesh_in = nullptr, Solver* solver = nullptr); - ~LaplaceMultigrid(){}; - - using Laplacian::setCoefA; - using Laplacian::setCoefC; - using Laplacian::setCoefC1; - using Laplacian::setCoefC2; - using Laplacian::setCoefD; - using Laplacian::setCoefEx; - using Laplacian::setCoefEz; - - void setCoefA(const Field2D& val) override { - ASSERT1(val.getLocation() == location); - ASSERT1(localmesh == val.getMesh()); - A = val; - } - void setCoefC(const Field2D& val) override { - ASSERT1(val.getLocation() == location); - ASSERT1(localmesh == val.getMesh()); - C1 = val; - C2 = val; - } - void setCoefC1(const Field2D& val) override { - ASSERT1(val.getLocation() == location); - ASSERT1(localmesh == val.getMesh()); - C1 = val; - } - void setCoefC2(const Field2D& val) override { - ASSERT1(val.getLocation() == location); - ASSERT1(localmesh == val.getMesh()); - C2 = val; - } - void setCoefD(const Field2D& val) override { - ASSERT1(val.getLocation() == location); - ASSERT1(localmesh == val.getMesh()); - D = val; - } - void setCoefEx(const Field2D& UNUSED(val)) override { - throw BoutException("setCoefEx is not implemented in LaplaceMultigrid"); - } - void setCoefEz(const Field2D& UNUSED(val)) override { - throw BoutException("setCoefEz is not implemented in LaplaceMultigrid"); - } - - void setCoefA(const Field3D& val) override { - ASSERT1(val.getLocation() == location); - ASSERT1(localmesh == val.getMesh()); - A = val; - } - void setCoefC(const Field3D& val) override { - ASSERT1(val.getLocation() == location); - ASSERT1(localmesh == val.getMesh()); - C1 = val; - C2 = val; - } - void setCoefC1(const Field3D& val) override { - ASSERT1(val.getLocation() == location); - ASSERT1(localmesh == val.getMesh()); - C1 = val; - } - void setCoefC2(const Field3D& val) override { - ASSERT1(val.getLocation() == location); - ASSERT1(localmesh == val.getMesh()); - C2 = val; - } - void setCoefD(const Field3D& val) override { - ASSERT1(val.getLocation() == location); - ASSERT1(localmesh == val.getMesh()); - D = val; - } - - bool uses3DCoefs() const override { return true; } - - using Laplacian::solve; - FieldPerp solve(const FieldPerp& b) override { - ASSERT1(localmesh == b.getMesh()); - - return solve(b, zeroFrom(b)); - } - FieldPerp solve(const FieldPerp& b_in, const FieldPerp& x0) override; - -private: - Field3D A, C1, C2, D; // ODE Coefficients - int Nx_local, Nx_global, Nz_local, Nz_global; // Local and global grid sizes - int yindex; // y-position of the current solution phase - Array x; // solution vector - Array b; // RHS vector - std::unique_ptr kMG; - - /******* Start implementation ********/ - int mglevel, mgplag, cftype, mgsm, pcheck; - int mgcount, mgmpi; - - Options* opts; - BoutReal rtol, atol, dtol, omega; - MPI_Comm commX; - - int comms_tagbase; - - void generateMatrixF(int); -}; - -namespace { -RegisterLaplace registerlaplacemultigrid(LAPLACE_MULTIGRID); -} - -#endif // BOUT_USE_METRIC_3D - -#endif // BOUT_MULTIGRID_LAPLACE_H diff --git a/src/invert/laplace/impls/multigrid/multigrid_solver.cxx b/src/invert/laplace/impls/multigrid/multigrid_solver.cxx deleted file mode 100644 index 074e482320..0000000000 --- a/src/invert/laplace/impls/multigrid/multigrid_solver.cxx +++ /dev/null @@ -1,771 +0,0 @@ -/************************************************************************** - * Perpendicular Laplacian inversion. - * Using Geometrical Multigrid Solver - * - * Equation solved is: - d*\nabla^2_\perp x + (1/c1)\nabla_perp c2\cdot\nabla_\perp x + a x = b - * - ************************************************************************** - * Copyright 2015 K.S. Kang - * - * 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 . - * - **************************************************************************/ - -#include "bout/build_defines.hxx" - -#if not BOUT_USE_METRIC_3D - -#include "multigrid_laplace.hxx" -#include "bout/unused.hxx" -#include - -#include -#include - -Multigrid1DP::Multigrid1DP(int level, int lx, int lz, int gx, int dl, int merge, - MPI_Comm comm, int check) - : MultigridAlg(level, lx, lz, gx, lz, comm, check) { - - mglevel = level; - - if (pcheck > 0) { - output << "Construct MG1DP " << endl; - } - commMG = comm; - MPI_Comm_size(commMG, &xNP); - MPI_Comm_rank(commMG, &rProcI); - xProcI = rProcI; - if (xNP > 1) { - if (xProcI == 0) { - xProcM = xNP - 1; - } else { - xProcM = xProcI - 1; - } - if (xProcI == xNP - 1) { - xProcP = 0; - } else { - xProcP = xProcI + 1; - } - } else { - xProcI = 0; - xProcM = 0; - xProcP = 0; - } - zNP = 1; - numP = xNP; - zProcI = xProcI; - zProcP = xProcI; - zProcM = xProcI; - - if (pcheck == 1) { - output << "In MG1DP level " << mglevel << " xNP=" << xNP << "(" << xProcI << ")" - << endl; - output << "lest level is " << dl << "(" << numP << ")" << endl; - } - - int nz, kk, nx; - if (dl > 0) { - // Find levels for more coarser spaces - if (numP > merge) { - int nn = numP; - int mm = static_cast(sqrt(numP)); - kk = 1; - for (int n = nn; n > 1; n--) { - if (nn % 2 != 0) { - n = 1; - } else { - kk += 1; - nn = nn / 2; - } - } - nn = lnz[0]; - int kz = 1; - for (int n = nn; n > 1; n--) { - if (nn % 2 != 0) { - n = 1; - } else { - kz += 1; - nn = nn / 2; - } - } - - if (kz < kk) { - kk = kz; - } - nz = 1; - nx = xNP; - for (int n = 0; n < kk; n++) { - if (nz * 2 <= mm) { - nz = 2 * nz; - nx = nx / 2; - } else { - n = kk; - } - } - - lx = gnx[0] / nx; - lz = lnz[0] / nz; - // 2D parallelization - kk = 1; - int llx = lx; - int llz = lz; - for (int n = dl; n > 0; n--) { - if ((llx % 2 == 0) && (llz % 2 == 0)) { - kk += 1; - llx = llx / 2; - llz = llz / 2; - } else { - n = 1; - } - } - if (kk > 0) { - kflag = 1; - } else { - kflag = 2; - } - } else { - kflag = 2; - } - if (kflag == 1) { - if (pcheck == 1) { - output << "To MG2DP " << kk << "xNP=" << nx << "(" << nz << ")" << endl; - output << "lest level is " << dl - kk + 1 << "(" << lx << ", " << lz << ")" - << endl; - } - int colors = rProcI / nz; - int keys = rProcI / nz; - MPI_Comm_split(commMG, colors, keys, &comm2D); - rMG = bout::utils::make_unique( - kk, lx, lz, gnx[0], lnz[0], dl - kk + 1, nx, nz, commMG, pcheck); - } else { - int nn = gnx[0]; - int mm = gnz[0]; - int kk = 1; - for (int n = dl; n > 0; n--) { - if ((nn % 2 == 0) && (mm % 2 == 0)) { - kk += 1; - nn = nn / 2; - mm = mm / 2; - } else { - n = 1; - } - } - if (pcheck == 1) { - output << "To Ser " << kk << " xNP=" << xNP << "(" << zNP << ")" << endl; - output << kflag << " total dim " << gnx[0] << "(" << lnz[0] << ")" << endl; - } - sMG = bout::utils::make_unique(kk, gnx[0], lnz[0], commMG, pcheck); - } - } else { - kflag = 0; - } -} - -void Multigrid1DP::setMultigridC(int UNUSED(plag)) { - - int level = mglevel - 1; - for (int n = level; n > 0; n--) { - if (pcheck == 2) { - output << n << "matrix in 1DP = " << lnx[n - 1] << "," << lnz[n - 1] << endl; - output << gnx[n - 1] << "," << gnz[n - 1] << endl; - } - setMatrixC(n); - } - - if (kflag == 1) { - level = rMG->mglevel - 1; - convertMatrixF2D(level); - - if (level > 0) { - rMG->setMultigridC(0); - } - - if (pcheck == 2) { - for (int i = level; i >= 0; i--) { - std::string outfile = fmt::format("2DP_matC{:1d}_{:d}.mat", i, rMG->rProcI); - output << "Out file= " << outfile << endl; - FILE* outf = fopen(outfile.c_str(), "w"); - int dim = (rMG->lnx[i] + 2) * (rMG->lnz[i] + 2); - fprintf(outf, "dim = %d (%d, %d)\n", dim, rMG->lnx[i], rMG->lnz[i]); - - for (int ii = 0; ii < dim; ii++) { - fprintf(outf, "%d ==", ii); - for (int j = 0; j < 9; j++) { - fprintf(outf, "%12.6f,", rMG->matmg[i][ii * 9 + j]); - } - fprintf(outf, "\n"); - } - fclose(outf); - } - } - } else if (kflag == 2) { - level = sMG->mglevel - 1; - convertMatrixFS(level); - - if (level > 0) { - sMG->setMultigridC(0); - } - if (pcheck == 3) { - for (int i = level; i >= 0; i--) { - std::string outfile = fmt::format("S1D_matC{:1d}_{:d}.mat", i, sMG->rProcI); - output << "Out file= " << outfile << endl; - FILE* outf = fopen(outfile.c_str(), "w"); - int dim = (sMG->lnx[i] + 2) * (sMG->lnz[i] + 2); - fprintf(outf, "dim = %d\n", dim); - - for (int ii = 0; ii < dim; ii++) { - fprintf(outf, "%d ==", ii); - for (int j = 0; j < 9; j++) { - fprintf(outf, "%12.6f,", sMG->matmg[i][ii * 9 + j]); - } - fprintf(outf, "\n"); - } - fclose(outf); - } - } - } -} - -void Multigrid1DP::setValueS() { - - if (kflag == 1) { - rMG->mgplag = mgplag; - rMG->mgsm = mgsm; - rMG->cftype = cftype; - rMG->rtol = rtol; - rMG->atol = atol; - rMG->dtol = dtol; - rMG->omega = omega; - rMG->setValueS(); - } else if (kflag == 2) { - sMG->mgplag = mgplag; - sMG->mgsm = mgsm; - sMG->cftype = cftype; - sMG->rtol = rtol; - sMG->atol = atol; - sMG->dtol = dtol; - sMG->omega = omega; - } -} - -void Multigrid1DP::setPcheck(int check) { - - pcheck = check; - if (kflag == 1) { - rMG->setPcheck(check); - } else if (kflag == 2) { - sMG->pcheck = check; - } -} - -void Multigrid1DP::lowestSolver(BoutReal* x, BoutReal* b, int UNUSED(plag)) { - - if (kflag == 1) { - int level = rMG->mglevel - 1; - int dim = (rMG->lnx[level] + 2) * (rMG->lnz[level] + 2); - Array y(dim); - Array r(dim); - - int ggx = rMG->lnx[level]; - int dimg = (ggx + 2) * (gnz[0] + 2); - Array yl(dimg); - Array yg(dimg); - - int nx = (xProcI % rMG->zNP) * lnx[0]; - - BOUT_OMP_PERF(parallel default(shared)) - { - BOUT_OMP_PERF(for) - for (int i = 0; i < dim; i++) { - y[i] = 0.0; - r[i] = 0.0; - } - - BOUT_OMP_PERF(for) - for (int i = 0; i < dimg; i++) { - yl[i] = 0.0; - yg[i] = 0.0; - } - - int xend = lnx[0] + 1; - int zend = lnz[0] + 1; - BOUT_OMP_PERF(for collapse(2)) - for (int ix = 1; ix < xend; ix++) { - for (int iz = 1; iz < zend; iz++) { - int nn = (nx + ix) * (lnz[0] + 2) + iz; - int mm = ix * (lnz[0] + 2) + iz; - yl[nn] = b[mm]; - } - } - } - bout::globals::mpi->MPI_Allreduce(std::begin(yl), std::begin(yg), dimg, MPI_DOUBLE, - MPI_SUM, comm2D); - - int nz = (xProcI % rMG->zNP) * (rMG->lnz[level]); - BOUT_OMP_PERF(parallel default(shared)) - { - int xend = rMG->lnx[level] + 1; - int zend = rMG->lnz[level] + 1; - BOUT_OMP_PERF(for collapse(2)) - for (int ix = 1; ix < xend; ix++) { - for (int iz = 1; iz < zend; iz++) { - int nn = ix * (lnz[0] + 2) + nz + iz; - int mm = ix * (rMG->lnz[level] + 2) + iz; - r[mm] = yg[nn]; - } - } - } - - rMG->getSolution(std::begin(y), std::begin(r), 1); - - BOUT_OMP_PERF(parallel default(shared)) - { - BOUT_OMP_PERF(for) - for (int i = 0; i < dimg; i++) { - yl[i] = 0.0; - yg[i] = 0.0; - } - - int xend = rMG->lnx[level] + 1; - int zend = rMG->lnz[level] + 1; - BOUT_OMP_PERF(for collapse(2)) - for (int ix = 1; ix < xend; ix++) { - for (int iz = 1; iz < zend; iz++) { - int nn = ix * (lnz[0] + 2) + nz + iz; - int mm = ix * (rMG->lnz[level] + 2) + iz; - yl[nn] = y[mm]; - } - } - } - bout::globals::mpi->MPI_Allreduce(std::begin(yl), std::begin(yg), dimg, MPI_DOUBLE, - MPI_SUM, comm2D); - - BOUT_OMP_PERF(parallel default(shared)) - { - int xend = lnx[0] + 1; - int zend = lnz[0] + 1; - BOUT_OMP_PERF(for collapse(2)) - for (int ix = 1; ix < xend; ix++) { - for (int iz = 1; iz < zend; iz++) { - int nn = (nx + ix) * (lnz[0] + 2) + iz; - int mm = ix * (lnz[0] + 2) + iz; - x[mm] = yg[nn]; - } - } - } - communications(x, 0); - } else if (kflag == 2) { - int level = sMG->mglevel - 1; - int dim = (sMG->lnx[level] + 2) * (sMG->lnz[level] + 2); - Array y(dim); - Array r(dim); - int nx = xProcI * lnx[0]; - BOUT_OMP_PERF(parallel default(shared)) - { - BOUT_OMP_PERF(for) - for (int i = 0; i < dim; i++) { - y[i] = 0.0; - r[i] = 0.0; - } - int xend = lnx[0] + 1; - int zend = lnz[0] + 1; - BOUT_OMP_PERF(for collapse(2)) - for (int ix = 1; ix < xend; ix++) { - for (int iz = 1; iz < zend; iz++) { - int nn = (nx + ix) * (lnz[0] + 2) + iz; - int mm = ix * (lnz[0] + 2) + iz; - y[nn] = b[mm]; - } - } - } - bout::globals::mpi->MPI_Allreduce(std::begin(y), std::begin(r), dim, MPI_DOUBLE, - MPI_SUM, commMG); - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int i = 0; i < dim; i++) { - y[i] = 0.0; - } - sMG->getSolution(std::begin(y), std::begin(r), 1); - - BOUT_OMP_PERF(parallel default(shared)) - { - int xend = lnx[0] + 1; - int zend = lnz[0] + 1; - BOUT_OMP_PERF(for collapse(2)) - for (int ix = 1; ix < xend; ix++) { - for (int iz = 1; iz < zend; iz++) { - int nn = (nx + ix) * (lnz[0] + 2) + iz; - int mm = ix * (lnz[0] + 2) + iz; - x[mm] = y[nn]; - } - } - } - communications(x, 0); - } else { - pGMRES(x, b, 0, 0); - } -} - -void Multigrid1DP::convertMatrixF2D(int level) { - - int ggx = rMG->lnx[level]; - int dim = (ggx + 2) * (gnz[0] + 2); - Array yl(dim * 9); - Array yg(dim * 9); - int nx = (xProcI % rMG->zNP) * lnx[0]; - BOUT_OMP_PERF(parallel default(shared)) - { - BOUT_OMP_PERF(for) - for (int i = 0; i < dim * 9; i++) { - yl[i] = 0.0; - yg[i] = 0.0; - } - BOUT_OMP_PERF(for) - for (int i = 0; i < (rMG->lnx[level] + 2) * (rMG->lnz[level] + 2) * 9; i++) { - rMG->matmg[level][i] = 0.0; - } - - int xend = lnx[0] + 1; - int zend = lnz[0] + 1; - BOUT_OMP_PERF(for collapse(2)) - for (int ix = 1; ix < xend; ix++) { - for (int iz = 1; iz < zend; iz++) { - int nn = (nx + ix) * (lnz[0] + 2) + iz; - int mm = ix * (lnz[0] + 2) + iz; - for (int k = 0; k < 9; k++) { - yl[nn * 9 + k] = matmg[0][mm * 9 + k]; - } - } - } - } - if (pcheck == 3) { - std::string outfile = fmt::format("2DP_CP_{}.mat", rProcI); - output << "Out file= " << outfile << endl; - FILE* outf = fopen(outfile.c_str(), "w"); - fprintf(outf, "dim = (%d, %d)\n", ggx, gnz[0]); - - for (int ii = 0; ii < dim; ii++) { - fprintf(outf, "%d ==", ii); - for (int j = 0; j < 9; j++) { - fprintf(outf, "%12.6f,", yl[ii * 9 + j]); - } - fprintf(outf, "\n"); - } - fclose(outf); - } - bout::globals::mpi->MPI_Allreduce(std::begin(yl), std::begin(yg), dim * 9, MPI_DOUBLE, - MPI_SUM, comm2D); - - if (pcheck == 3) { - - std::string outfile = fmt::format("2DP_Conv_{}.mat", rProcI); - output << "Out file= " << outfile << endl; - FILE* outf = fopen(outfile.c_str(), "w"); - fprintf(outf, "dim = (%d, %d)\n", ggx, gnz[0]); - - for (int ii = 0; ii < dim; ii++) { - fprintf(outf, "%d ==", ii); - for (int j = 0; j < 9; j++) { - fprintf(outf, "%12.6f,", yg[ii * 9 + j]); - } - fprintf(outf, "\n"); - } - fclose(outf); - } - int nz = (xProcI % rMG->zNP) * (rMG->lnz[level]); - - BOUT_OMP_PERF(parallel default(shared)) - { - int xend = rMG->lnx[level] + 1; - int zend = rMG->lnz[level] + 1; - BOUT_OMP_PERF(for collapse(2)) - for (int ix = 1; ix < xend; ix++) { - for (int iz = 1; iz < zend; iz++) { - int nn = ix * (lnz[0] + 2) + nz + iz; - int mm = ix * (rMG->lnz[level] + 2) + iz; - for (int k = 0; k < 9; k++) { - rMG->matmg[level][mm * 9 + k] = yg[nn * 9 + k]; - } - } - } - } -} - -void Multigrid1DP::convertMatrixFS(int level) { - - int dim = (gnx[0] + 2) * (gnz[0] + 2); - Array yl(dim * 9); - BoutReal* yg = sMG->matmg[level]; - int nx = xProcI * lnx[0]; - BOUT_OMP_PERF(parallel default(shared)) - { - BOUT_OMP_PERF(for) - for (int i = 0; i < dim * 9; i++) { - yl[i] = 0.0; - yg[i] = 0.0; - } - int xend = lnx[0] + 1; - int zend = lnz[0] + 1; - BOUT_OMP_PERF(for collapse(2)) - for (int ix = 1; ix < xend; ix++) { - for (int iz = 1; iz < zend; iz++) { - int nn = (nx + ix) * (lnz[0] + 2) + iz; - int mm = ix * (lnz[0] + 2) + iz; - for (int k = 0; k < 9; k++) { - yl[nn * 9 + k] = matmg[0][mm * 9 + k]; - } - } - } - } - bout::globals::mpi->MPI_Allreduce(std::begin(yl), yg, dim * 9, MPI_DOUBLE, MPI_SUM, - commMG); -} - -Multigrid2DPf1D::Multigrid2DPf1D(int level, int lx, int lz, int gx, int gz, int dl, - int px, int pz, MPI_Comm comm, int check) - : MultigridAlg(level, lx, lz, gx, gz, comm, check) { - - mglevel = level; - - /* Momory allocate for Multigrid */ - xNP = px; - zNP = pz; - numP = px * pz; - commMG = comm; - MPI_Comm_rank(comm, &rProcI); - xProcI = rProcI / zNP; - zProcI = rProcI % zNP; - if (xProcI == 0) { - xProcM = numP - zNP + zProcI; - } else { - xProcM = rProcI - zNP; - } - if (xProcI == xNP - 1) { - xProcP = zProcI; - } else { - xProcP = rProcI + zNP; - } - if (zProcI == 0) { - zProcM = rProcI + zNP - 1; - } else { - zProcM = rProcI - 1; - } - if (zProcI == zNP - 1) { - zProcP = xProcI * zNP; - } else { - zProcP = rProcI + 1; - } - if (pcheck == 2) { - output << "In 2DP " << mglevel << "xNP=" << xNP << "(" << zNP << ")" << dl << endl; - for (int i = mglevel - 1; i >= 0; i--) { - output << i << " loc dim " << lnx[i] << "," << lnz[i] << endl; - output << i << " glo dim " << gnx[i] << "," << gnz[i] << endl; - } - } - - if (dl > 0) { - int nn = gnx[0]; - int mm = gnz[0]; - int kk = 1; - for (int n = dl; n > 0; n--) { - if ((nn % 2 == 0) && (mm % 2 == 0)) { - kk += 1; - nn = nn / 2; - mm = mm / 2; - } else { - n = 0; - } - } - if (pcheck == 2) { - output << "In 2DP To Ser" << kk << "xNP=" << xNP << "(" << zNP << ")" << endl; - output << "total dim" << gnx[0] << "(" << gnz[0] << ")" << endl; - } - kflag = 2; - sMG = bout::utils::make_unique(kk, gnx[0], gnz[0], commMG, pcheck); - } else { - kflag = 0; - } -} - -void Multigrid2DPf1D::setMultigridC(int UNUSED(plag)) { - - int level = mglevel - 1; - for (int n = level; n > 0; n--) { - setMatrixC(n); - if (pcheck == 2) { - output << n << "Cmatrix in 2DP = " << lnx[n - 1] << "," << lnz[n - 1] << endl; - output << gnx[n - 1] << "," << gnz[n - 1] << endl; - } - } - if (kflag == 2) { - level = sMG->mglevel - 1; - convertMatrixFS(level); - if (level > 0) { - sMG->setMultigridC(0); - } - if (pcheck == 2) { - for (int i = level; i >= 0; i--) { - std::string outfile = fmt::format("S2D_matC{:1d}_{:d}.mat", i, sMG->rProcI); - FILE* outf = fopen(outfile.c_str(), "w"); - output << "Out file= " << outfile << endl; - int dim = (sMG->lnx[i] + 2) * (sMG->lnz[i] + 2); - fprintf(outf, "dim = %d\n", dim); - - for (int ii = 0; ii < dim; ii++) { - fprintf(outf, "%d ==", ii); - for (int j = 0; j < 9; j++) { - fprintf(outf, "%12.6f,", sMG->matmg[i][ii * 9 + j]); - } - fprintf(outf, "\n"); - } - fclose(outf); - } - } - } -} - -void Multigrid2DPf1D::setValueS() { - if (kflag == 2) { - sMG->mgplag = mgplag; - sMG->mgsm = mgsm; - sMG->cftype = cftype; - sMG->rtol = rtol; - sMG->atol = atol; - sMG->dtol = dtol; - sMG->omega = omega; - } -} - -void Multigrid2DPf1D::setPcheck(int check) { - - pcheck = check; - if (kflag == 2) { - sMG->pcheck = check; - } -} - -void Multigrid2DPf1D::lowestSolver(BoutReal* x, BoutReal* b, int UNUSED(plag)) { - - if (kflag == 2) { - int level = sMG->mglevel - 1; - int dim = (sMG->lnx[level] + 2) * (sMG->lnz[level] + 2); - Array y(dim); - Array r(dim); - int nx = xProcI * lnx[0]; - int nz = zProcI * lnz[0]; - BOUT_OMP_PERF(parallel default(shared)) - { - BOUT_OMP_PERF(for) - for (int i = 0; i < dim; i++) { - y[i] = 0.0; - r[i] = 0.0; - } - - int xend = lnx[0] + 1; - int zend = lnz[0] + 1; - BOUT_OMP_PERF(for collapse(2)) - for (int ix = 1; ix < xend; ix++) { - for (int iz = 1; iz < zend; iz++) { - int nn = (nx + ix) * (gnz[0] + 2) + nz + iz; - int mm = ix * (lnz[0] + 2) + iz; - y[nn] = b[mm]; - } - } - } - bout::globals::mpi->MPI_Allreduce(std::begin(y), std::begin(r), dim, MPI_DOUBLE, - MPI_SUM, commMG); - BOUT_OMP_PERF(parallel default(shared)) - BOUT_OMP_PERF(for) - for (int i = 0; i < dim; i++) { - y[i] = 0.0; - } - sMG->getSolution(std::begin(y), std::begin(r), 1); - BOUT_OMP_PERF(parallel default(shared)) - { - int xend = lnx[0] + 1; - int zend = lnz[0] + 1; - BOUT_OMP_PERF(for collapse(2)) - for (int ix = 1; ix < xend; ix++) { - for (int iz = 1; iz < zend; iz++) { - int nn = (nx + ix) * (gnz[0] + 2) + nz + iz; - int mm = ix * (lnz[0] + 2) + iz; - x[mm] = y[nn]; - } - } - } - communications(x, 0); - } else { - pGMRES(x, b, 0, 0); - } -} - -void Multigrid2DPf1D::convertMatrixFS(int level) { - - int dim = (gnx[0] + 2) * (gnz[0] + 2); - Array yl(dim * 9); - BoutReal* yg = sMG->matmg[level]; - int nx = xProcI * lnx[0]; - int nz = zProcI * lnz[0]; - BOUT_OMP_PERF(parallel default(shared)) - { - BOUT_OMP_PERF(for) - for (int i = 0; i < dim * 9; i++) { - yl[i] = 0.0; - yg[i] = 0.0; - } - int xend = lnx[0] + 1; - int zend = lnz[0] + 1; - BOUT_OMP_PERF(for collapse(2)) - for (int ix = 1; ix < xend; ix++) { - for (int iz = 1; iz < zend; iz++) { - int nn = (nx + ix) * (gnz[0] + 2) + nz + iz; - int mm = ix * (lnz[0] + 2) + iz; - for (int k = 0; k < 9; k++) { - yl[nn * 9 + k] = matmg[0][mm * 9 + k]; - } - } - } - } - bout::globals::mpi->MPI_Allreduce(std::begin(yl), yg, dim * 9, MPI_DOUBLE, MPI_SUM, - commMG); -} - -MultigridSerial::MultigridSerial(int level, int gx, int gz, MPI_Comm comm, int check) - : MultigridAlg(level, gx, gz, gx, gz, comm, check) { - - xNP = 1; - zNP = 1; - numP = 1; - MPI_Comm_rank(commMG, &rProcI); - xProcI = rProcI; - zProcI = rProcI; - xProcM = rProcI; - xProcP = rProcI; - zProcM = rProcI; - zProcP = rProcI; - if (pcheck == 2) { - output << "In SerMG " << mglevel << "xNP=" << xNP << "(" << zNP << ")" << endl; - for (int i = mglevel - 1; i >= 0; i--) { - output << i << " Ser loc dim " << lnx[i] << "," << lnz[i] << endl; - output << i << " Ser glo dim " << gnx[i] << "," << gnz[i] << endl; - } - } -} - -#endif diff --git a/src/invert/laplace/invert_laplace.cxx b/src/invert/laplace/invert_laplace.cxx index 0ec5f471da..b4cb8fe064 100644 --- a/src/invert/laplace/invert_laplace.cxx +++ b/src/invert/laplace/invert_laplace.cxx @@ -51,7 +51,6 @@ #include "impls/cyclic/cyclic_laplace.hxx" #include "impls/hypre3d/hypre3d_laplace.hxx" #include "impls/iterative_parallel_tri/iterative_parallel_tri.hxx" -#include "impls/multigrid/multigrid_laplace.hxx" #include "impls/naulin/naulin_laplace.hxx" #include "impls/pcr/pcr.hxx" #include "impls/pcr_thomas/pcr_thomas.hxx" diff --git a/tests/integrated/CMakeLists.txt b/tests/integrated/CMakeLists.txt index 59a3091971..e158d480be 100644 --- a/tests/integrated/CMakeLists.txt +++ b/tests/integrated/CMakeLists.txt @@ -30,7 +30,6 @@ add_subdirectory(test-laplacexy2-hypre) add_subdirectory(test-laplacexy-fv) add_subdirectory(test-laplacexy-short) add_subdirectory(test-laplacexz) -add_subdirectory(test-multigrid_laplace) add_subdirectory(test-naulin-laplace) add_subdirectory(test-options-netcdf) add_subdirectory(test-petsc_laplace) diff --git a/tests/integrated/test-multigrid_laplace/CMakeLists.txt b/tests/integrated/test-multigrid_laplace/CMakeLists.txt deleted file mode 100644 index b64b0e166c..0000000000 --- a/tests/integrated/test-multigrid_laplace/CMakeLists.txt +++ /dev/null @@ -1,9 +0,0 @@ -bout_add_integrated_test( - test-multigrid-laplace - SOURCES test_multigrid_laplace.cxx - CONFLICTS BOUT_USE_METRIC_3D - USE_RUNTEST USE_DATA_BOUT_INP - EXTRA_FILES data/BOUT_jy4.inp data/BOUT_jy63.inp data/BOUT_jy127.inp - data/BOUT_unsheared.inp - PROCESSORS 3 -) diff --git a/tests/integrated/test-multigrid_laplace/data/BOUT.inp b/tests/integrated/test-multigrid_laplace/data/BOUT.inp deleted file mode 100644 index 1f3133b781..0000000000 --- a/tests/integrated/test-multigrid_laplace/data/BOUT.inp +++ /dev/null @@ -1,121 +0,0 @@ - -[mesh] -nx = 100 -nz = 128 -ny = 2 -ixseps1 = -1 -ixseps2 = 0 -dx = 2.66390416073e-05 -dy = 0.0490873865783 -dz = 0.001636246219277382 -g11 = 0.0350116938353*(1.+.01*y*dy) -g22 = 2.25850105285645*(1.+.01*y*dy) -g33 = 2.65228796005*(1.+.01*y*dy) -g12 = 0.0 -g13 = -0.0381551384926*(1.+.01*y*dy) -g23 = -2.20591616630554*(1.+.01*y*dy) -g_11 = 31.1654663086/(1.+.01*y*dy) -g_22 = 2.5341296196/(1.+.01*y*dy) -g_33 = 2.19225430489/(1.+.01*y*dy) -g_12 = 2.33345580101013/(1.+.01*y*dy) -g_13 = 2.38908100128/(1.+.01*y*dy) -g_23 = 2.14121198654175/(1.+.01*y*dy) - -############################################# - -[laplace] -type = multigrid -all_terms = true -rtol = 1e-10 -# Note, set a fairly high atol here because the RHS in this test is large, O(10^6) -atol = 1e-10 -dtol = 1000000000000.0 -include_yguards = false - -############################################# -# input variables for tests -############################################# - -[f1] -# zero value x boundary conditions -p1 = 0.39503274 -q1 = 0.20974396 -function = sin(2*pi*3*x)*cos(z-p1) + sin(2*pi*x)*cos(3*z-q1) - -[d1] -r1 = 0.512547 -s1 = 0.30908712 -function = 1. + .1*cos(4.*x-r1)*sin(3*z-s1) - -[c1] -t1 = 0.18439023 -u1 = 0.401089473 -function = 1. + .1*sin(3.2*x-t1)*sin(2*z-u1) - -[a1] -w1 = 0.612547 -x1 = 0.30908712 -function = -1. + .1*sin(1.43*x-w1)*sin(5*z-x1) - -[f2] -# zero gradient x boundary conditions -p2 = 0.39503274 -q2 = 0.20974396 -function = cos(2*pi*3*x)*cos(z-p2) + cos(2*pi*x)*cos(3*z-q2) - -[d2] -r2 = 0.512547 -s2 = 0.30908712 -function = 1. + .1*cos(4.*x-r2)*sin(3*z-s2) - -[c2] -t2 = 0.18439023 -u2 = 0.401089473 -function = 1. + .1*sin(3.2*x-t2)*sin(2*z-u2) - -[a2] -w2 = 0.612547 -x2 = 0.30908712 -function = -1. + .1*sin(1.43*x-w2)*sin(5*z-x2) - -[f3] -# set value x boundary conditions -p3 = 0.39503274 -q3 = 0.20974396 -function = sin(2*pi*3*x-q3)*cos(z-p3) + sin(2*pi*x-p3)*cos(3*z-q3) - -[d3] -r3 = 0.512547 -s3 = 0.30908712 -function = 1. + .1*cos(4.*x-r3)*sin(3*z-s3) - -[c3] -t3 = 0.18439023 -u3 = 0.401089473 -function = 1. + .1*sin(3.2*x-t3)*sin(2*z-u3) - -[a3] -w3 = 0.612547 -x3 = 0.30908712 -function = -1. + .1*sin(1.43*x-w3)*sin(5*z-x3) - -[f4] -# set gradient x boundary conditions -p4 = 0.39503274 -q4 = 0.20974396 -function = sin(2*pi*3*x+q4)*cos(z-p4) + sin(2*pi*x+p4)*cos(3*z-q4) - -[d4] -r4 = 0.512547 -s4 = 0.30908712 -function = 1. + .1*cos(4.*x-r4)*sin(3*z-s4) - -[c4] -t4 = 0.18439023 -u4 = 0.401089473 -function = 1. + .1*sin(3.2*x-t4)*sin(2*z-u4) - -[a4] -w4 = 0.612547 -x4 = 0.30908712 -function = -1. + .1*sin(1.43*x-w4)*sin(5*z-x4) diff --git a/tests/integrated/test-multigrid_laplace/data/BOUT_jy127.inp b/tests/integrated/test-multigrid_laplace/data/BOUT_jy127.inp deleted file mode 100644 index ccd807f51e..0000000000 --- a/tests/integrated/test-multigrid_laplace/data/BOUT_jy127.inp +++ /dev/null @@ -1,123 +0,0 @@ -dump_format = "nc" - -[mesh] -nx = 132 -nz = 128 -ny = 2 -ixseps1 = -1 -ixseps2 = 0 -dx = 2.66390416073e-05 -dy = 0.0490873865783 -dz = 0.001636246219277382 -g11 = 0.00443281093612313*(1.+.01*y*dy) -g22 = 2.25850105285645*(1.+.01*y*dy) -g33 = 110.63729095459*(1.+.01*y*dy) -g12 = 0.0 -g13 = 0.616832256317139*(1.+.01*y*dy) -g23 = -7.38632392883301*(1.+.01*y*dy) -g_11 = 30129.087890625/(1.+.01*y*dy) -g_22 = 16.9609527587891/(1.+.01*y*dy) -g_33 = 1.54435157775879/(1.+.01*y*dy) -g_12 = -702.816772460938/(1.+.01*y*dy) -g_13 = -214.898834228516/(1.+.01*y*dy) -g_23 = 5.05073070526123/(1.+.01*y*dy) - -############################################# - -[laplace] -type = multigrid -all_terms = true -rtol = 1e-09 -# Note, set a fairly high atol here because the RHS in this test is large, O(10^6) -atol = 1e-10 -dtol = 1000000000000.0 -include_yguards = false -maxits = 1000 - -############################################# -# input variables for tests -############################################# - -[f1] -# zero value x boundary conditions -p1 = 0.39503274 -q1 = 0.20974396 -function = sin(2*pi*3*x)*cos(z-p1) + sin(2*pi*x)*cos(3*z-q1) - -[d1] -r1 = 0.512547 -s1 = 0.30908712 -function = 1. + .1*cos(4.*x-r1)*sin(3*z-s1) - -[c1] -t1 = 0.18439023 -u1 = 0.401089473 -function = 1. + .1*sin(3.2*x-t1)*sin(2*z-u1) - -[a1] -w1 = 0.612547 -x1 = 0.30908712 -function = -1. + .1*sin(1.43*x-w1)*sin(5*z-x1) - -[f2] -# zero gradient x boundary conditions -p2 = 0.39503274 -q2 = 0.20974396 -function = cos(2*pi*3*x)*cos(z-p2) + cos(2*pi*x)*cos(3*z-q2) - -[d2] -r2 = 0.512547 -s2 = 0.30908712 -function = 1. + .1*cos(4.*x-r2)*sin(3*z-s2) - -[c2] -t2 = 0.18439023 -u2 = 0.401089473 -function = 1. + .1*sin(3.2*x-t2)*sin(2*z-u2) - -[a2] -w2 = 0.612547 -x2 = 0.30908712 -function = -1. + .1*sin(1.43*x-w2)*sin(5*z-x2) - -[f3] -# set value x boundary conditions -p3 = 0.39503274 -q3 = 0.20974396 -function = sin(2*pi*3*x-q3)*cos(z-p3) + sin(2*pi*x-p3)*cos(3*z-q3) - -[d3] -r3 = 0.512547 -s3 = 0.30908712 -function = 1. + .1*cos(4.*x-r3)*sin(3*z-s3) - -[c3] -t3 = 0.18439023 -u3 = 0.401089473 -function = 1. + .1*sin(3.2*x-t3)*sin(2*z-u3) - -[a3] -w3 = 0.612547 -x3 = 0.30908712 -function = -1. + .1*sin(1.43*x-w3)*sin(5*z-x3) - -[f4] -# set gradient x boundary conditions -p4 = 0.39503274 -q4 = 0.20974396 -function = sin(2*pi*3*x+q4)*cos(z-p4) + sin(2*pi*x+p4)*cos(3*z-q4) - -[d4] -r4 = 0.512547 -s4 = 0.30908712 -function = 1. + .1*cos(4.*x-r4)*sin(3*z-s4) - -[c4] -t4 = 0.18439023 -u4 = 0.401089473 -function = 1. + .1*sin(3.2*x-t4)*sin(2*z-u4) - -[a4] -w4 = 0.612547 -x4 = 0.30908712 -function = -1. + .1*sin(1.43*x-w4)*sin(5*z-x4) diff --git a/tests/integrated/test-multigrid_laplace/data/BOUT_jy4.inp b/tests/integrated/test-multigrid_laplace/data/BOUT_jy4.inp deleted file mode 100644 index 6660237cc9..0000000000 --- a/tests/integrated/test-multigrid_laplace/data/BOUT_jy4.inp +++ /dev/null @@ -1,123 +0,0 @@ -dump_format = "nc" - -[mesh] -nx = 132 -nz = 128 -ny = 2 -ixseps1 = -1 -ixseps2 = 0 -dx = 2.66390416072682e-05 -dy = 0.0490873865783 -dz = 0.001636246219277382 -g11 = 0.00432878965511918*(1.+.01*y*dy) -g22 = 2.25850105285645*(1.+.01*y*dy) -g33 = 108.512680053711*(1.+.01*y*dy) -g12 = 0.0 -g13 = -0.602102994918823*(1.+.01*y*dy) -g23 = -7.38268566131592*(1.+.01*y*dy) -g_11 = 30857.376953125/(1.+.01*y*dy) -g_22 = 17.357889175415/(1.+.01*y*dy) -g_33 = 1.5830215215683/(1.+.01*y*dy) -g_12 = 719.755859375/(1.+.01*y*dy) -g_13 = 220.186721801758/(1.+.01*y*dy) -g_23 = 5.17464876174927/(1.+.01*y*dy) - -############################################# - -[laplace] -type = multigrid -all_terms = true -rtol = 1e-09 -# Note, set a fairly high atol here because the RHS in this test is large, O(10^6) -atol = 1e-10 -dtol = 1000000000000.0 -include_yguards = false -maxits = 1000 - -############################################# -# input variables for tests -############################################# - -[f1] -# zero value x boundary conditions -p1 = 0.39503274 -q1 = 0.20974396 -function = sin(2*pi*3*x)*cos(z-p1) + sin(2*pi*x)*cos(3*z-q1) - -[d1] -r1 = 0.512547 -s1 = 0.30908712 -function = 1. + .1*cos(4.*x-r1)*sin(3*z-s1) - -[c1] -t1 = 0.18439023 -u1 = 0.401089473 -function = 1. + .1*sin(3.2*x-t1)*sin(2*z-u1) - -[a1] -w1 = 0.612547 -x1 = 0.30908712 -function = -1. + .1*sin(1.43*x-w1)*sin(5*z-x1) - -[f2] -# zero gradient x boundary conditions -p2 = 0.39503274 -q2 = 0.20974396 -function = cos(2*pi*3*x)*cos(z-p2) + cos(2*pi*x)*cos(3*z-q2) - -[d2] -r2 = 0.512547 -s2 = 0.30908712 -function = 1. + .1*cos(4.*x-r2)*sin(3*z-s2) - -[c2] -t2 = 0.18439023 -u2 = 0.401089473 -function = 1. + .1*sin(3.2*x-t2)*sin(2*z-u2) - -[a2] -w2 = 0.612547 -x2 = 0.30908712 -function = -1. + .1*sin(1.43*x-w2)*sin(5*z-x2) - -[f3] -# set value x boundary conditions -p3 = 0.39503274 -q3 = 0.20974396 -function = sin(2*pi*3*x-q3)*cos(z-p3) + sin(2*pi*x-p3)*cos(3*z-q3) - -[d3] -r3 = 0.512547 -s3 = 0.30908712 -function = 1. + .1*cos(4.*x-r3)*sin(3*z-s3) - -[c3] -t3 = 0.18439023 -u3 = 0.401089473 -function = 1. + .1*sin(3.2*x-t3)*sin(2*z-u3) - -[a3] -w3 = 0.612547 -x3 = 0.30908712 -function = -1. + .1*sin(1.43*x-w3)*sin(5*z-x3) - -[f4] -# set gradient x boundary conditions -p4 = 0.39503274 -q4 = 0.20974396 -function = sin(2*pi*3*x+q4)*cos(z-p4) + sin(2*pi*x+p4)*cos(3*z-q4) - -[d4] -r4 = 0.512547 -s4 = 0.30908712 -function = 1. + .1*cos(4.*x-r4)*sin(3*z-s4) - -[c4] -t4 = 0.18439023 -u4 = 0.401089473 -function = 1. + .1*sin(3.2*x-t4)*sin(2*z-u4) - -[a4] -w4 = 0.612547 -x4 = 0.30908712 -function = -1. + .1*sin(1.43*x-w4)*sin(5*z-x4) diff --git a/tests/integrated/test-multigrid_laplace/data/BOUT_jy63.inp b/tests/integrated/test-multigrid_laplace/data/BOUT_jy63.inp deleted file mode 100644 index 3dc4d09ef5..0000000000 --- a/tests/integrated/test-multigrid_laplace/data/BOUT_jy63.inp +++ /dev/null @@ -1,123 +0,0 @@ -dump_format = "nc" - -[mesh] -nx = 132 -nz = 128 -ny = 2 -ixseps1 = -1 -ixseps2 = 0 -dx = 2.66390416073e-05 -dy = 0.0490873865783 -dz = 0.001636246219277382 -g11 = 0.0350116938353*(1.+.01*y*dy) -g22 = 2.25850105285645*(1.+.01*y*dy) -g33 = 2.65228796005*(1.+.01*y*dy) -g12 = 0.0 -g13 = -0.0381551384926*(1.+.01*y*dy) -g23 = -2.20591616630554*(1.+.01*y*dy) -g_11 = 31.1654663086/(1.+.01*y*dy) -g_22 = 2.5341296196/(1.+.01*y*dy) -g_33 = 2.19225430489/(1.+.01*y*dy) -g_12 = 2.33345580101013/(1.+.01*y*dy) -g_13 = 2.38908100128/(1.+.01*y*dy) -g_23 = 2.14121198654175/(1.+.01*y*dy) - -############################################# - -[laplace] -type = multigrid -all_terms = true -rtol = 1e-09 -# Note, set a fairly high atol here because the RHS in this test is large, O(10^6) -atol = 1e-10 -dtol = 1000000000000.0 -include_yguards = false -maxits = 1000 - -############################################# -# input variables for tests -############################################# - -[f1] -# zero value x boundary conditions -p1 = 0.39503274 -q1 = 0.20974396 -function = sin(2*pi*3*x)*cos(z-p1) + sin(2*pi*x)*cos(3*z-q1) - -[d1] -r1 = 0.512547 -s1 = 0.30908712 -function = 1. + .1*cos(4.*x-r1)*sin(3*z-s1) - -[c1] -t1 = 0.18439023 -u1 = 0.401089473 -function = 1. + .1*sin(3.2*x-t1)*sin(2*z-u1) - -[a1] -w1 = 0.612547 -x1 = 0.30908712 -function = -1. + .1*sin(1.43*x-w1)*sin(5*z-x1) - -[f2] -# zero gradient x boundary conditions -p2 = 0.39503274 -q2 = 0.20974396 -function = cos(2*pi*3*x)*cos(z-p2) + cos(2*pi*x)*cos(3*z-q2) - -[d2] -r2 = 0.512547 -s2 = 0.30908712 -function = 1. + .1*cos(4.*x-r2)*sin(3*z-s2) - -[c2] -t2 = 0.18439023 -u2 = 0.401089473 -function = 1. + .1*sin(3.2*x-t2)*sin(2*z-u2) - -[a2] -w2 = 0.612547 -x2 = 0.30908712 -function = -1. + .1*sin(1.43*x-w2)*sin(5*z-x2) - -[f3] -# set value x boundary conditions -p3 = 0.39503274 -q3 = 0.20974396 -function = sin(2*pi*3*x-q3)*cos(z-p3) + sin(2*pi*x-p3)*cos(3*z-q3) - -[d3] -r3 = 0.512547 -s3 = 0.30908712 -function = 1. + .1*cos(4.*x-r3)*sin(3*z-s3) - -[c3] -t3 = 0.18439023 -u3 = 0.401089473 -function = 1. + .1*sin(3.2*x-t3)*sin(2*z-u3) - -[a3] -w3 = 0.612547 -x3 = 0.30908712 -function = -1. + .1*sin(1.43*x-w3)*sin(5*z-x3) - -[f4] -# set gradient x boundary conditions -p4 = 0.39503274 -q4 = 0.20974396 -function = sin(2*pi*3*x+q4)*cos(z-p4) + sin(2*pi*x+p4)*cos(3*z-q4) - -[d4] -r4 = 0.512547 -s4 = 0.30908712 -function = 1. + .1*cos(4.*x-r4)*sin(3*z-s4) - -[c4] -t4 = 0.18439023 -u4 = 0.401089473 -function = 1. + .1*sin(3.2*x-t4)*sin(2*z-u4) - -[a4] -w4 = 0.612547 -x4 = 0.30908712 -function = -1. + .1*sin(1.43*x-w4)*sin(5*z-x4) diff --git a/tests/integrated/test-multigrid_laplace/data/BOUT_unsheared.inp b/tests/integrated/test-multigrid_laplace/data/BOUT_unsheared.inp deleted file mode 100644 index 222020256c..0000000000 --- a/tests/integrated/test-multigrid_laplace/data/BOUT_unsheared.inp +++ /dev/null @@ -1,114 +0,0 @@ -#grid = "grids/flat_grid.nc" -dump_format = "nc" - -[mesh] -nx = 100 -ny = 1 -mz = 128 - -myg = 0 - -############################################# - -[laplace] -type = multigrid -all_terms = true -rtol = 1e-10 -include_yguards = false -maxits = 1000 - -############################################# - -[SPT] -#type=spt -all_terms = true -#flags=15 -include_yguards = false - -############################################# -# input variables for tests -############################################# - -[f1] -# zero value x boundary conditions -p1 = 0.39503274 -q1 = 0.20974396 -function = sin(2*pi*3*x)*cos(z-p1) + sin(2*pi*x)*cos(3*z-q1) - -[d1] -r1 = 0.512547 -s1 = 0.30908712 -function = 1. + .1*cos(4.*x-r1)*sin(3*z-s1) - -[c1] -t1 = 0.18439023 -u1 = 0.401089473 -function = 1. + .1*sin(3.2*x-t1)*sin(2*z-u1) - -[a1] -w1 = 0.612547 -x1 = 0.30908712 -function = -1. + .1*sin(1.43*x-w1)*sin(5*z-x1) - -[f2] -# zero gradient x boundary conditions -p2 = 0.39503274 -q2 = 0.20974396 -function = cos(2*pi*3*x)*cos(z-p2) + cos(2*pi*x)*cos(3*z-q2) - -[d2] -r2 = 0.512547 -s2 = 0.30908712 -function = 1. + .1*cos(4.*x-r2)*sin(3*z-s2) - -[c2] -t2 = 0.18439023 -u2 = 0.401089473 -function = 1. + .1*sin(3.2*x-t2)*sin(2*z-u2) - -[a2] -w2 = 0.612547 -x2 = 0.30908712 -function = -1. + .1*sin(1.43*x-w2)*sin(5*z-x2) - -[f3] -# set value x boundary conditions -p3 = 0.39503274 -q3 = 0.20974396 -function = sin(2*pi*3*x-q3)*cos(z-p3) + sin(2*pi*x-p3)*cos(3*z-q3) - -[d3] -r3 = 0.512547 -s3 = 0.30908712 -function = 1. + .1*cos(4.*x-r3)*sin(3*z-s3) - -[c3] -t3 = 0.18439023 -u3 = 0.401089473 -function = 1. + .1*sin(3.2*x-t3)*sin(2*z-u3) - -[a3] -w3 = 0.612547 -x3 = 0.30908712 -function = -1. + .1*sin(1.43*x-w3)*sin(5*z-x3) - -[f4] -# set gradient x boundary conditions -p4 = 0.39503274 -q4 = 0.20974396 -function = sin(2*pi*3*x+q4)*cos(z-p4) + sin(2*pi*x+p4)*cos(3*z-q4) - -[d4] -r4 = 0.512547 -s4 = 0.30908712 -function = 1. + .1*cos(4.*x-r4)*sin(3*z-s4) - -[c4] -t4 = 0.18439023 -u4 = 0.401089473 -function = 1. + .1*sin(3.2*x-t4)*sin(2*z-u4) - -[a4] -w4 = 0.612547 -x4 = 0.30908712 -function = -1. + .1*sin(1.43*x-w4)*sin(5*z-x4) diff --git a/tests/integrated/test-multigrid_laplace/plotresults.py b/tests/integrated/test-multigrid_laplace/plotresults.py deleted file mode 100755 index 37d594e85d..0000000000 --- a/tests/integrated/test-multigrid_laplace/plotresults.py +++ /dev/null @@ -1,53 +0,0 @@ -#!/usr/bin/env python3 - -from boutdata.data import BoutOutputs -from matplotlib import pyplot -from sys import exit - -Ntests = 4 - -d = BoutOutputs("data") -# d = BoutOutputs("data", caching=True, info=False) - -for i in range(1, Ntests + 1): - pyplot.figure(i) - pyplot.subplot(231) - pyplot.pcolor(d["f" + str(i)][1:-1, 0, :].T) - pyplot.colorbar() - pyplot.title("f" + str(i)) - pyplot.xlabel("X") - pyplot.ylabel("Z") - pyplot.subplot(232) - pyplot.pcolor(d["sol" + str(i)][1:-1, 0, :].T) - pyplot.colorbar() - pyplot.title("sol" + str(i)) - pyplot.xlabel("X") - pyplot.ylabel("Z") - pyplot.subplot(233) - pyplot.pcolor(d["absolute_error" + str(i)][1:-1, 0, :].T) - pyplot.colorbar() - pyplot.title("absolute_error" + str(i)) - pyplot.xlabel("X") - pyplot.ylabel("Z") - pyplot.subplot(234) - pyplot.pcolor(d["b" + str(i)][2:-2, 0, :].T) - pyplot.colorbar() - pyplot.title("b" + str(i)) - pyplot.xlabel("X") - pyplot.ylabel("Z") - pyplot.subplot(235) - pyplot.pcolor(d["bcheck" + str(i)][2:-2, 0, :].T) - pyplot.colorbar() - pyplot.title("bcheck" + str(i)) - pyplot.xlabel("X") - pyplot.ylabel("Z") - pyplot.subplot(236) - pyplot.pcolor(d["b" + str(i)][2:-2, 0, :].T - d["bcheck" + str(i)][2:-2, 0, :].T) - pyplot.colorbar() - pyplot.title("b" + str(i) + "-bcheck" + str(i)) - pyplot.xlabel("X") - pyplot.ylabel("Z") - -pyplot.show() - -exit(0) diff --git a/tests/integrated/test-multigrid_laplace/test_multigrid_laplace.cxx b/tests/integrated/test-multigrid_laplace/test_multigrid_laplace.cxx deleted file mode 100644 index b9f80a6ae0..0000000000 --- a/tests/integrated/test-multigrid_laplace/test_multigrid_laplace.cxx +++ /dev/null @@ -1,326 +0,0 @@ -/************************************************************************** - * Testing Perpendicular Laplacian inversion using multigrid solver - * - ************************************************************************** - * Copyright 2013 J.T. Omotani - * - * 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 . - * - **************************************************************************/ - -#include -#include -#include -#include -#include -#include -#include -#include - -BoutReal max_error_at_ystart(const Field3D& error); -Field3D this_Grad_perp_dot_Grad_perp(const Field3D& f, const Field3D& g); - -int main(int argc, char** argv) { - - // Initialise BOUT++, setting up mesh - BoutInitialise(argc, argv); - - auto invert = Laplacian::create(); - - // Solving equations of the form d*Grad_perp2(f) + 1/c*Grad_perp(c).Grad_perp(f) + a*f = b for various boundary conditions - Field3D f1, a1, b1, c1, d1, sol1, bcheck1; - Field3D absolute_error1; - BoutReal max_error1; //Output of test - - using bout::globals::mesh; - - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - // Test 1: zero-value Dirichlet boundaries - f1 = FieldFactory::get()->create3D("f1:function", Options::getRoot(), mesh); - d1 = FieldFactory::get()->create3D("d1:function", Options::getRoot(), mesh); - c1 = FieldFactory::get()->create3D("c1:function", Options::getRoot(), mesh); - a1 = FieldFactory::get()->create3D("a1:function", Options::getRoot(), mesh); - - b1 = d1 * Delp2(f1, CELL_DEFAULT, false) + this_Grad_perp_dot_Grad_perp(c1, f1) / c1 - + a1 * f1; - sol1 = 0.; - - invert->setInnerBoundaryFlags(0); - invert->setOuterBoundaryFlags(0); - invert->setCoefA(a1); - invert->setCoefC(c1); - invert->setCoefD(d1); - - try { - sol1 = invert->solve(sliceXZ(b1, mesh->ystart)); - mesh->communicate(sol1); - checkData(sol1); - bcheck1 = d1 * Delp2(sol1, CELL_DEFAULT, false) - + this_Grad_perp_dot_Grad_perp(c1, sol1) / c1 + a1 * sol1; - absolute_error1 = f1 - sol1; - max_error1 = max_error_at_ystart(abs(absolute_error1, "RGN_NOBNDRY")); - } catch (BoutException& err) { - output << "BoutException occured in invert->solve(b1): " << err.what() << endl - << "Laplacian inversion failed to converge (probably)" << endl; - max_error1 = -1; - sol1 = -1.; - bcheck1 = -1.; - absolute_error1 = -1.; - } - - d1 = Delp2(f1, CELL_DEFAULT, false); - c1 = this_Grad_perp_dot_Grad_perp(c1, f1) / c1; - a1 = a1 * f1; - - output << endl << "Test 1: zero Dirichlet" << endl; - output << "Magnitude of maximum absolute error is " << max_error1 << endl; - - Options dump; - dump["a1"] = a1; - dump["b1"] = b1; - dump["c1"] = c1; - dump["d1"] = d1; - dump["f1"] = f1; - dump["sol1"] = sol1; - dump["bcheck1"] = bcheck1; - dump["absolute_error1"] = absolute_error1; - dump["max_error1"] = max_error1; - - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - Field3D f2, a2, b2, c2, d2, sol2, bcheck2; - Field3D absolute_error2; - BoutReal max_error2; //Output of test - // Test 2: zero-value Neumann boundaries - f2 = FieldFactory::get()->create3D("f2:function", Options::getRoot(), mesh); - d2 = FieldFactory::get()->create3D("d2:function", Options::getRoot(), mesh); - c2 = FieldFactory::get()->create3D("c2:function", Options::getRoot(), mesh); - a2 = FieldFactory::get()->create3D("a2:function", Options::getRoot(), mesh); - - b2 = d2 * Delp2(f2, CELL_DEFAULT, false) + this_Grad_perp_dot_Grad_perp(c2, f2) / c2 - + a2 * f2; - sol2 = 0.; - - invert->setInnerBoundaryFlags(INVERT_AC_GRAD); - invert->setOuterBoundaryFlags(INVERT_AC_GRAD); - invert->setCoefA(a2); - invert->setCoefC(c2); - invert->setCoefD(d2); - - try { - sol2 = invert->solve(sliceXZ(b2, mesh->ystart)); - mesh->communicate(sol2); - bcheck2 = d2 * Delp2(sol2, CELL_DEFAULT, false) - + this_Grad_perp_dot_Grad_perp(c2, sol2) / c2 + a2 * sol2; - absolute_error2 = f2 - sol2; - max_error2 = max_error_at_ystart(abs(absolute_error2, "RGN_NOBNDRY")); - } catch (BoutException& err) { - output << "BoutException occured in invert->solve(b2): " << err.what() << endl - << "Laplacian inversion failed to converge (probably)" << endl; - max_error2 = -1; - sol2 = -1.; - bcheck2 = -1.; - absolute_error2 = -1.; - } - - output << endl << "Test 2: zero Neumann" << endl; - output << "Magnitude of maximum absolute error is " << max_error2 << endl; - - dump["a2"] = a2; - dump["b2"] = b2; - dump["c2"] = c2; - dump["d2"] = d2; - dump["f2"] = f2; - dump["sol2"] = sol2; - dump["bcheck2"] = bcheck2; - dump["absolute_error2"] = absolute_error2; - dump["max_error2"] = max_error2; - - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - Field3D f3, a3, b3, c3, d3, sol3, bcheck3; - Field3D absolute_error3; - BoutReal max_error3; //Output of test - // Test 3: set-value Dirichlet boundaries - f3 = FieldFactory::get()->create3D("f3:function", Options::getRoot(), mesh); - d3 = FieldFactory::get()->create3D("d3:function", Options::getRoot(), mesh); - c3 = FieldFactory::get()->create3D("c3:function", Options::getRoot(), mesh); - a3 = FieldFactory::get()->create3D("a3:function", Options::getRoot(), mesh); - - b3 = d3 * Delp2(f3, CELL_DEFAULT, false) + this_Grad_perp_dot_Grad_perp(c3, f3) / c3 - + a3 * f3; - sol3 = 0.; - - invert->setInnerBoundaryFlags(INVERT_SET); - invert->setOuterBoundaryFlags(INVERT_SET); - invert->setCoefA(a3); - invert->setCoefC(c3); - invert->setCoefD(d3); - - // make field to pass in boundary conditions - Field3D x0 = 0.; - if (mesh->firstX()) { - for (int k = mesh->zstart; k <= mesh->zend; k++) { - x0(mesh->xstart - 1, mesh->ystart, k) = - 0.5 - * (f3(mesh->xstart - 1, mesh->ystart, k) + f3(mesh->xstart, mesh->ystart, k)); - } - } - if (mesh->lastX()) { - for (int k = mesh->zstart; k <= mesh->zend; k++) { - x0(mesh->xend + 1, mesh->ystart, k) = - 0.5 * (f3(mesh->xend + 1, mesh->ystart, k) + f3(mesh->xend, mesh->ystart, k)); - } - } - - try { - sol3 = invert->solve(sliceXZ(b3, mesh->ystart), sliceXZ(x0, mesh->ystart)); - mesh->communicate(sol3); - bcheck3 = d3 * Delp2(sol3, CELL_DEFAULT, false) - + this_Grad_perp_dot_Grad_perp(c3, f3) / c3 + a3 * sol3; - absolute_error3 = f3 - sol3; - max_error3 = max_error_at_ystart(abs(absolute_error3, "RGN_NOBNDRY")); - } catch (BoutException& err) { - output << "BoutException occured in invert->solve(b3): " << err.what() << endl - << "Laplacian inversion failed to converge (probably)" << endl; - max_error3 = -1; - sol3 = -1.; - bcheck3 = -1.; - absolute_error3 = -1.; - } - - output << endl << "Test 3: set Dirichlet" << endl; - output << "Magnitude of maximum absolute error is " << max_error3 << endl; - - dump["a3"] = a3; - dump["b3"] = b3; - dump["c3"] = c3; - dump["d3"] = d3; - dump["f3"] = f3; - dump["sol3"] = sol3; - dump["bcheck3"] = bcheck3; - dump["absolute_error3"] = absolute_error3; - dump["max_error3"] = max_error3; - - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - Field3D f4, a4, b4, c4, d4, sol4, bcheck4; - Field3D absolute_error4; - BoutReal max_error4; //Output of test - // Test 4: set-value Neumann boundaries - f4 = FieldFactory::get()->create3D("f4:function", Options::getRoot(), mesh); - d4 = FieldFactory::get()->create3D("d4:function", Options::getRoot(), mesh); - c4 = FieldFactory::get()->create3D("c4:function", Options::getRoot(), mesh); - a4 = FieldFactory::get()->create3D("a4:function", Options::getRoot(), mesh); - - b4 = d4 * Delp2(f4, CELL_DEFAULT, false) + this_Grad_perp_dot_Grad_perp(c4, f4) / c4 - + a4 * f4; - sol4 = 0.; - - invert->setInnerBoundaryFlags(INVERT_AC_GRAD + INVERT_SET); - invert->setOuterBoundaryFlags(INVERT_AC_GRAD + INVERT_SET); - invert->setCoefA(a4); - invert->setCoefC(c4); - invert->setCoefD(d4); - - // make field to pass in boundary conditions - x0 = 0.; - if (mesh->firstX()) { - for (int k = mesh->zstart; k <= mesh->zend; k++) { - x0(mesh->xstart - 1, mesh->ystart, k) = - (f4(mesh->xstart, mesh->ystart, k) - f4(mesh->xstart - 1, mesh->ystart, k)) - / mesh->getCoordinates()->dx(mesh->xstart, mesh->ystart, k) - / sqrt(mesh->getCoordinates()->g_11(mesh->xstart, mesh->ystart, k)); - } - } - if (mesh->lastX()) { - for (int k = mesh->zstart; k <= mesh->zend; k++) { - x0(mesh->xend + 1, mesh->ystart, k) = - (f4(mesh->xend + 1, mesh->ystart, k) - f4(mesh->xend, mesh->ystart, k)) - / mesh->getCoordinates()->dx(mesh->xend, mesh->ystart, k) - / sqrt(mesh->getCoordinates()->g_11(mesh->xend, mesh->ystart, k)); - } - } - - try { - sol4 = invert->solve(sliceXZ(b4, mesh->ystart), sliceXZ(x0, mesh->ystart)); - mesh->communicate(sol4); - bcheck4 = d4 * Delp2(sol4, CELL_DEFAULT, false) - + this_Grad_perp_dot_Grad_perp(c4, sol4) / c4 + a4 * sol4; - absolute_error4 = f4 - sol4; - max_error4 = max_error_at_ystart(abs(absolute_error4, "RGN_NOBNDRY")); - } catch (BoutException& err) { - output << "BoutException occured in invert->solve(b4): " << err.what() << endl - << "Laplacian inversion failed to converge (probably)" << endl; - max_error4 = -1; - sol4 = -1.; - bcheck4 = -1.; - absolute_error4 = -1.; - } - - output << endl << "Test 4: set Neumann" << endl; - output << "Magnitude of maximum absolute error is " << max_error4 << endl; - - dump["a4"] = a4; - dump["b4"] = b4; - dump["c4"] = c4; - dump["d4"] = d4; - dump["f4"] = f4; - dump["sol4"] = sol4; - dump["bcheck4"] = bcheck4; - dump["absolute_error4"] = absolute_error4; - dump["max_error4"] = max_error4; - //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - - output << "\nFinished running test.\n\n"; - - bout::writeDefaultOutputFile(dump); - - MPI_Barrier(BoutComm::get()); // Wait for all processors to write data - - bout::checkForUnusedOptions(); - BoutFinalise(); - return 0; -} - -Field3D this_Grad_perp_dot_Grad_perp(const Field3D& f, const Field3D& g) { - auto* mesh = f.getMesh(); - - Field3D result = mesh->getCoordinates()->g11() * ::DDX(f) * ::DDX(g) - + mesh->getCoordinates()->g33() * ::DDZ(f) * ::DDZ(g) - + mesh->getCoordinates()->g13() * (DDX(f) * DDZ(g) + DDZ(f) * DDX(g)); - - return result; -} - -BoutReal max_error_at_ystart(const Field3D& error) { - auto* mesh = error.getMesh(); - - BoutReal local_max_error = error(mesh->xstart, mesh->ystart, 0); - - for (int jx = mesh->xstart; jx <= mesh->xend; jx++) { - for (int jz = mesh->zstart; jz <= mesh->zend; jz++) { - if (local_max_error < error(jx, mesh->ystart, jz)) { - local_max_error = error(jx, mesh->ystart, jz); - } - } - } - - BoutReal max_error; - - MPI_Allreduce(&local_max_error, &max_error, 1, MPI_DOUBLE, MPI_MAX, BoutComm::get()); - - return max_error; -} diff --git a/tests/integrated/test-multigrid_laplace/test_multigrid_laplace.py b/tests/integrated/test-multigrid_laplace/test_multigrid_laplace.py deleted file mode 100755 index d29a7596f3..0000000000 --- a/tests/integrated/test-multigrid_laplace/test_multigrid_laplace.py +++ /dev/null @@ -1,51 +0,0 @@ -#!/usr/bin/env python3 - -# -# Run the test, check the error -# - -import pytest -from boututils.run_wrapper import shell, launch_safe -from boutdata.collect import collect - -tol = 2e-6 # Absolute tolerance -numTests = 4 # We test 4 different boundary conditions (with slightly different inputs for each) - - -@pytest.mark.parametrize("nproc", [1, 2, 4]) -@pytest.mark.parametrize( - "inputfile", ["BOUT_jy4.inp", "BOUT_jy63.inp", "BOUT_jy127.inp"] -) -def test_multigrid_laplace(nproc, inputfile): - # set nxpe on the command line as we only use solution from one point in y, - # so splitting in y-direction is redundant (and also doesn't help test the multigrid solver) - cmd = f"./test_multigrid_laplace -f {inputfile} NXPE={nproc} input:error_on_unused_options=false" - - shell(["rm -f data/BOUT.dmp.*.nc"]) - - print("Running multigrid Laplacian inversion test") - print(f"{nproc} processors, input file is: {inputfile}") - - s, out = launch_safe(cmd, nproc=nproc, pipe=True) - - # Save a log file uniquely named for this specific combination - input_file_name = inputfile.replace(".inp", "") - with open(f"run.log.{nproc}.{input_file_name}", "w") as f: - f.write(out) - - # Collect errors - failures = [] - for i in range(1, numTests + 1): - var_name = f"max_error{i}" - e = collect(var_name, path="data", info=False) - - if e < 0.0: - failures.append( - f"Sub-test {i} ({var_name}): Solver did not converge (error = {e})" - ) - elif e > tol: - failures.append( - f"Sub-test {i} ({var_name}): Absolute error {e} exceeds tolerance {tol}" - ) - - assert not failures, "Some failed tests:\n" + "\n".join(failures) diff --git a/tests/integrated/test-multigrid_laplace/test_multigrid_laplace_unsheared.py b/tests/integrated/test-multigrid_laplace/test_multigrid_laplace_unsheared.py deleted file mode 100755 index 627db84371..0000000000 --- a/tests/integrated/test-multigrid_laplace/test_multigrid_laplace_unsheared.py +++ /dev/null @@ -1,52 +0,0 @@ -#!/usr/bin/env python3 - -# -# Run the test, check the error -# - -import pytest -from boututils.run_wrapper import shell, launch_safe -from boutdata.collect import collect - -tol = 1e-9 # Absolute tolerance -numTests = 4 # We test 4 different boundary conditions (with slightly different inputs for each) - - -@pytest.mark.parametrize("nproc", [1, 3]) -def test_multigrid_laplace_unsheared(nproc): - # Make sure we don't use too many cores: - # Reduce number of OpenMP threads when using multiple MPI processes - mthread = 1 if nproc > 1 else 2 - - # set nxpe on the command line as we only use solution from one point in y, - # so splitting in y-direction is redundant (and also doesn't help test the multigrid solver) - inputfile = "BOUT_unsheared.inp" - cmd = f"./test_multigrid_laplace -f {inputfile} NXPE={nproc}" - - shell(["rm -f data/BOUT.dmp.*.nc"]) - - print( - f"Running unsheared multigrid Laplacian test: {nproc} procs, {mthread} thread(s)" - ) - - s, out = launch_safe(cmd, nproc=nproc, mthread=mthread, pipe=True) - - with open(f"run.log.{nproc}", "w") as f: - f.write(out) - - # Collect errors - failures = [] - for i in range(1, numTests + 1): - var_name = f"max_error{i}" - e = collect(var_name, path="data", info=False) - - if e < 0.0: - failures.append( - f"Sub-test {i} ({var_name}): Solver did not converge (error = {e})" - ) - elif e > tol: - failures.append( - f"Sub-test {i} ({var_name}): Absolute error {e:.4e} exceeds tolerance {tol:.4e}" - ) - - assert not failures, "Some failed tests:\n" + "\n".join(failures)