From 1447e98a8a3abd43c404b5444c3b2daca4059bf6 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 14 Nov 2018 10:03:53 +0000 Subject: [PATCH 1/4] Use Field::getCoordinates() where possible in Laplacian solvers Neater than previous mesh->getCoordinates(location) calls. --- src/invert/laplace/impls/cyclic/cyclic_laplace.cxx | 2 +- src/invert/laplace/impls/naulin/naulin_laplace.cxx | 2 +- src/invert/parderiv/impls/cyclic/cyclic.cxx | 2 +- src/invert/parderiv/impls/serial/serial.cxx | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx b/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx index a802b4beef..3661b95dda 100644 --- a/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx +++ b/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx @@ -258,7 +258,7 @@ const Field3D LaplaceCyclic::solve(const Field3D &rhs, const Field3D &x0) { x.allocate(); x.setLocation(location); - Coordinates *coord = mesh->getCoordinates(location); + Coordinates *coord = rhs.getCoordinates(); // Get the width of the boundary diff --git a/src/invert/laplace/impls/naulin/naulin_laplace.cxx b/src/invert/laplace/impls/naulin/naulin_laplace.cxx index fa065b0f12..9aad382d78 100644 --- a/src/invert/laplace/impls/naulin/naulin_laplace.cxx +++ b/src/invert/laplace/impls/naulin/naulin_laplace.cxx @@ -172,7 +172,7 @@ const Field3D LaplaceNaulin::solve(const Field3D &rhs, const Field3D &x0) { ASSERT1(Acoef.getLocation() == location); Mesh *mesh = rhs.getMesh(); - Coordinates *coords = mesh->getCoordinates(location); + Coordinates *coords = rhs.getCoordinates(); Field3D x(x0); // Result Field3D rhsOverD = rhs/Dcoef; diff --git a/src/invert/parderiv/impls/cyclic/cyclic.cxx b/src/invert/parderiv/impls/cyclic/cyclic.cxx index 29d4bc8dde..f6d0a3eb2f 100644 --- a/src/invert/parderiv/impls/cyclic/cyclic.cxx +++ b/src/invert/parderiv/impls/cyclic/cyclic.cxx @@ -89,7 +89,7 @@ const Field3D InvertParCR::solve(const Field3D &f) { result.allocate(); result.setLocation(f.getLocation()); - Coordinates *coord = mesh->getCoordinates(f.getLocation()); + Coordinates *coord = f.getCoordinates(); // Create cyclic reduction object CyclicReduce *cr = diff --git a/src/invert/parderiv/impls/serial/serial.cxx b/src/invert/parderiv/impls/serial/serial.cxx index 7cc0fd128e..200650efd2 100644 --- a/src/invert/parderiv/impls/serial/serial.cxx +++ b/src/invert/parderiv/impls/serial/serial.cxx @@ -67,7 +67,7 @@ const Field3D InvertParSerial::solve(const Field3D &f) { result.allocate(); result.setLocation(f.getLocation()); - Coordinates *coord = mesh->getCoordinates(f.getLocation()); + Coordinates *coord = f.getCoordinates(); // Loop over flux-surfaces SurfaceIter surf(mesh); From 7965a6208e6b6f15fff6643037df474d8cb70e06 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 14 Nov 2018 10:35:41 +0000 Subject: [PATCH 2/4] Fix location setting in more Laplacian solvers Pass the location to delp2solver in LaplaceNaulin. Set location of result in solve() method of LaplaceXY, LaplaceXZcyclic and LaplaceXZpetsc. --- src/invert/laplace/impls/naulin/naulin_laplace.cxx | 2 +- src/invert/laplacexy/laplacexy.cxx | 1 + src/invert/laplacexz/impls/cyclic/laplacexz-cyclic.cxx | 1 + src/invert/laplacexz/impls/petsc/laplacexz-petsc.cxx | 1 + 4 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/invert/laplace/impls/naulin/naulin_laplace.cxx b/src/invert/laplace/impls/naulin/naulin_laplace.cxx index 9aad382d78..81e1297899 100644 --- a/src/invert/laplace/impls/naulin/naulin_laplace.cxx +++ b/src/invert/laplace/impls/naulin/naulin_laplace.cxx @@ -135,7 +135,7 @@ LaplaceNaulin::LaplaceNaulin(Options *opt, const CELL_LOC loc) OPTION(opt, rtol, 1.e-7); OPTION(opt, atol, 1.e-20); OPTION(opt, maxits, 100); - delp2solver = create(opt->getSection("delp2solver")); + delp2solver = create(opt->getSection("delp2solver"), location); std::string delp2type; opt->getSection("delp2solver")->get("type", delp2type, "cyclic"); // Check delp2solver is using an FFT scheme, otherwise it will not exactly diff --git a/src/invert/laplacexy/laplacexy.cxx b/src/invert/laplacexy/laplacexy.cxx index 6793880b8d..aa3b2d47aa 100644 --- a/src/invert/laplacexy/laplacexy.cxx +++ b/src/invert/laplacexy/laplacexy.cxx @@ -634,6 +634,7 @@ const Field2D LaplaceXY::solve(const Field2D &rhs, const Field2D &x0) { Field2D result; result.allocate(); + result.setLocation(rhs.getLocation()); for(int x=mesh->xstart;x<= mesh->xend;x++) { for(int y=mesh->ystart;y<=mesh->yend;y++) { diff --git a/src/invert/laplacexz/impls/cyclic/laplacexz-cyclic.cxx b/src/invert/laplacexz/impls/cyclic/laplacexz-cyclic.cxx index 8567b0c740..115c325267 100644 --- a/src/invert/laplacexz/impls/cyclic/laplacexz-cyclic.cxx +++ b/src/invert/laplacexz/impls/cyclic/laplacexz-cyclic.cxx @@ -244,6 +244,7 @@ Field3D LaplaceXZcyclic::solve(const Field3D &rhs, const Field3D &x0) { Field3D result(mesh); result.allocate(); + result.setLocation(rhs.getLocation()); ind = 0; for(int y=mesh->ystart; y <= mesh->yend; y++) { diff --git a/src/invert/laplacexz/impls/petsc/laplacexz-petsc.cxx b/src/invert/laplacexz/impls/petsc/laplacexz-petsc.cxx index 6095919064..db433b6b54 100644 --- a/src/invert/laplacexz/impls/petsc/laplacexz-petsc.cxx +++ b/src/invert/laplacexz/impls/petsc/laplacexz-petsc.cxx @@ -589,6 +589,7 @@ Field3D LaplaceXZpetsc::solve(const Field3D &bin, const Field3D &x0in) { Field3D result; result.allocate(); + result.setLocation(bin.getLocation()); for (auto &it : slice) { /// Get y index From 447d543e91ed53f48141c2b5fcb10a07c8cba52b Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 14 Nov 2018 12:12:01 +0000 Subject: [PATCH 3/4] Remove usage of global 'mesh' pointer in Laplacian solvers Allow Laplacian solvers to be initialized with a local Mesh object, stored as Laplacian::localmesh. Use localmesh everywhere instead of global mesh, or mesh from Field*D.getMesh() to ensure consistency between initialization and solve methods. Includes ASSERT1 checks that 'localmesh == rhs.getMesh()', etc. --- include/invert_laplace.hxx | 21 +- include/invert_parderiv.hxx | 18 +- .../laplace/impls/cyclic/cyclic_laplace.cxx | 119 ++++----- .../laplace/impls/cyclic/cyclic_laplace.hxx | 5 +- .../impls/multigrid/multigrid_laplace.cxx | 75 +++--- .../impls/multigrid/multigrid_laplace.hxx | 16 +- .../laplace/impls/mumps/mumps_laplace.cxx | 226 +++++++++--------- .../laplace/impls/mumps/mumps_laplace.hxx | 18 +- .../laplace/impls/naulin/naulin_laplace.cxx | 30 +-- .../laplace/impls/naulin/naulin_laplace.hxx | 12 +- src/invert/laplace/impls/pdd/pdd.cxx | 142 +++++------ src/invert/laplace/impls/pdd/pdd.hxx | 7 +- .../laplace/impls/petsc/petsc_laplace.cxx | 159 ++++++------ .../laplace/impls/petsc/petsc_laplace.hxx | 18 +- .../laplace/impls/serial_band/serial_band.cxx | 54 +++-- .../laplace/impls/serial_band/serial_band.hxx | 5 +- .../laplace/impls/serial_tri/serial_tri.cxx | 32 +-- .../laplace/impls/serial_tri/serial_tri.hxx | 5 +- .../laplace/impls/shoot/shoot_laplace.cxx | 53 ++-- .../laplace/impls/shoot/shoot_laplace.hxx | 5 +- src/invert/laplace/impls/spt/spt.cxx | 171 ++++++------- src/invert/laplace/impls/spt/spt.hxx | 11 +- src/invert/laplace/invert_laplace.cxx | 74 +++--- src/invert/laplace/laplacefactory.cxx | 38 +-- src/invert/laplace/laplacefactory.hxx | 2 +- src/invert/parderiv/impls/cyclic/cyclic.cxx | 32 +-- src/invert/parderiv/impls/cyclic/cyclic.hxx | 27 ++- src/invert/parderiv/impls/serial/serial.cxx | 41 ++-- src/invert/parderiv/impls/serial/serial.hxx | 27 ++- src/invert/parderiv/invert_parderiv.cxx | 4 +- src/invert/parderiv/parderiv_factory.cxx | 14 +- src/invert/parderiv/parderiv_factory.hxx | 6 +- 32 files changed, 799 insertions(+), 668 deletions(-) diff --git a/include/invert_laplace.hxx b/include/invert_laplace.hxx index 77d43e5565..005cc61a78 100644 --- a/include/invert_laplace.hxx +++ b/include/invert_laplace.hxx @@ -105,14 +105,14 @@ const int INVERT_OUT_RHS = 32768; ///< Use input value in RHS at outer boundary /// Base class for Laplacian inversion class Laplacian { public: - Laplacian(Options *options = nullptr, const CELL_LOC loc = CELL_CENTRE); + Laplacian(Options *options = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh* mesh_in = mesh); virtual ~Laplacian() {} /// Set coefficients for inversion. Re-builds matrices if necessary virtual void setCoefA(const Field2D &val) = 0; virtual void setCoefA(const Field3D &val) { setCoefA(DC(val)); } virtual void setCoefA(BoutReal r) { - Field2D f(r); + Field2D f(r, localmesh); f.setLocation(location); setCoefA(f); } @@ -120,7 +120,7 @@ public: virtual void setCoefC(const Field2D &val) = 0; virtual void setCoefC(const Field3D &val) { setCoefC(DC(val)); } virtual void setCoefC(BoutReal r) { - Field2D f(r); + Field2D f(r, localmesh); f.setLocation(location); setCoefC(f); } @@ -130,7 +130,7 @@ public: } virtual void setCoefC1(const Field3D &val) { setCoefC1(DC(val)); } virtual void setCoefC1(BoutReal r) { - Field2D f(r); + Field2D f(r, localmesh); f.setLocation(location); setCoefC1(f); } @@ -140,7 +140,7 @@ public: } virtual void setCoefC2(const Field3D &val) { setCoefC2(DC(val)); } virtual void setCoefC2(BoutReal r) { - Field2D f(r); + Field2D f(r, localmesh); f.setLocation(location); setCoefC2(f); } @@ -148,7 +148,7 @@ public: virtual void setCoefD(const Field2D &val) = 0; virtual void setCoefD(const Field3D &val) { setCoefD(DC(val)); } virtual void setCoefD(BoutReal r) { - Field2D f(r); + Field2D f(r, localmesh); f.setLocation(location); setCoefD(f); } @@ -156,7 +156,7 @@ public: virtual void setCoefEx(const Field2D &val) = 0; virtual void setCoefEx(const Field3D &val) { setCoefEx(DC(val)); } virtual void setCoefEx(BoutReal r) { - Field2D f(r); + Field2D f(r, localmesh); f.setLocation(location); setCoefEx(f); } @@ -164,7 +164,7 @@ public: virtual void setCoefEz(const Field2D &val) = 0; virtual void setCoefEz(const Field3D &val) { setCoefEz(DC(val)); } virtual void setCoefEz(BoutReal r) { - Field2D f(r); + Field2D f(r, localmesh); f.setLocation(location); setCoefEz(f); } @@ -192,7 +192,7 @@ public: * * @param[in] opt The options section to use. By default "laplace" will be used */ - static Laplacian *create(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE); + static Laplacian *create(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = mesh); static Laplacian* defaultInstance(); ///< Return pointer to global singleton static void cleanup(); ///< Frees all memory @@ -227,7 +227,8 @@ protected: const Field2D *a, const Field2D *ccoef, const Field2D *d, bool includeguards=true); - CELL_LOC location; + CELL_LOC location; ///< staggered grid location of this solver + Mesh* localmesh; ///< Mesh object for this solver private: /// Singleton instance static Laplacian *instance; diff --git a/include/invert_parderiv.hxx b/include/invert_parderiv.hxx index 1ab6a9c4bb..0f562a2f79 100644 --- a/include/invert_parderiv.hxx +++ b/include/invert_parderiv.hxx @@ -64,7 +64,8 @@ public: * with pure virtual members, so can't be created directly. * To create an InvertPar object call the create() static function. */ - InvertPar(Options *UNUSED(opt)) {} + InvertPar(Options *UNUSED(opt), Mesh *mesh_in = mesh) + : localmesh(mesh_in) {} virtual ~InvertPar() {} /*! @@ -72,7 +73,7 @@ public: * * Note: For consistency this should be renamed "create" and take an Options* argument */ - static InvertPar* Create(); + static InvertPar* Create(Mesh *mesh_in = mesh); /*! * Solve the system of equations @@ -100,35 +101,38 @@ public: */ virtual void setCoefA(const Field2D &f) = 0; virtual void setCoefA(const Field3D &f) {setCoefA(DC(f));} - virtual void setCoefA(BoutReal f) {setCoefA(Field2D(f));} + virtual void setCoefA(BoutReal f) {setCoefA(Field2D(f, localmesh));} /*! * Set the Grad2_par2 coefficient B */ virtual void setCoefB(const Field2D &f) = 0; virtual void setCoefB(const Field3D &f) {setCoefB(DC(f));} - virtual void setCoefB(BoutReal f) {setCoefB(Field2D(f));} + virtual void setCoefB(BoutReal f) {setCoefB(Field2D(f, localmesh));} /*! * Set the D2DYDZ coefficient C */ virtual void setCoefC(const Field2D &f) = 0; virtual void setCoefC(const Field3D &f) {setCoefB(DC(f));} - virtual void setCoefC(BoutReal f) {setCoefB(Field2D(f));} + virtual void setCoefC(BoutReal f) {setCoefB(Field2D(f, localmesh));} /*! * Set the D2DZ2 coefficient D */ virtual void setCoefD(const Field2D &f) = 0; virtual void setCoefD(const Field3D &f) {setCoefB(DC(f));} - virtual void setCoefD(BoutReal f) {setCoefB(Field2D(f));} + virtual void setCoefD(BoutReal f) {setCoefB(Field2D(f, localmesh));} /*! * Set the DDY coefficient E */ virtual void setCoefE(const Field2D &f) = 0; virtual void setCoefE(const Field3D &f) {setCoefB(DC(f));} - virtual void setCoefE(BoutReal f) {setCoefB(Field2D(f));} + virtual void setCoefE(BoutReal f) {setCoefB(Field2D(f, localmesh));} + +protected: + Mesh* localmesh; ///< Mesh object for this solver private: }; diff --git a/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx b/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx index 3661b95dda..4d0df4e7c4 100644 --- a/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx +++ b/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx @@ -43,8 +43,8 @@ #include "cyclic_laplace.hxx" -LaplaceCyclic::LaplaceCyclic(Options *opt, const CELL_LOC loc) - : Laplacian(opt, loc), Acoef(0.0), Ccoef(1.0), Dcoef(1.0) { +LaplaceCyclic::LaplaceCyclic(Options *opt, const CELL_LOC loc, Mesh *mesh_in) + : Laplacian(opt, loc, mesh_in), Acoef(0.0), Ccoef(1.0), Dcoef(1.0) { Acoef.setLocation(location); Ccoef.setLocation(location); Dcoef.setLocation(location); @@ -54,7 +54,7 @@ LaplaceCyclic::LaplaceCyclic(Options *opt, const CELL_LOC loc) OPTION(opt, dst, false); if(dst) { - nmode = mesh->LocalNz-2; + nmode = localmesh->LocalNz-2; }else nmode = maxmode+1; // Number of Z modes. maxmode set in invert_laplace.cxx from options @@ -62,13 +62,13 @@ LaplaceCyclic::LaplaceCyclic(Options *opt, const CELL_LOC loc) // Allocate arrays - xs = mesh->xstart; // Starting X index - if(mesh->firstX() && !mesh->periodicX){ // Only want to include guard cells at boundaries (unless periodic in x) + xs = localmesh->xstart; // Starting X index + if(localmesh->firstX() && !localmesh->periodicX){ // Only want to include guard cells at boundaries (unless periodic in x) xs = 0; } - xe = mesh->xend; // Last X index - if(mesh->lastX() && !mesh->periodicX){ // Only want to include guard cells at boundaries (unless periodic in x) - xe = mesh->LocalNx-1; + xe = localmesh->xend; // Last X index + if(localmesh->lastX() && !localmesh->periodicX){ // Only want to include guard cells at boundaries (unless periodic in x) + xe = localmesh->LocalNx-1; } int n = xe - xs + 1; // Number of X points on this processor, // including boundaries but not guard cells @@ -80,8 +80,8 @@ LaplaceCyclic::LaplaceCyclic(Options *opt, const CELL_LOC loc) bcmplx = Matrix(nmode, n); // Create a cyclic reduction object, operating on dcomplex values - cr = new CyclicReduce(mesh->getXcomm(), n); - cr->setPeriodic(mesh->periodicX); + cr = new CyclicReduce(localmesh->getXcomm(), n); + cr->setPeriodic(localmesh->periodicX); } LaplaceCyclic::~LaplaceCyclic() { @@ -90,11 +90,12 @@ LaplaceCyclic::~LaplaceCyclic() { } const FieldPerp LaplaceCyclic::solve(const FieldPerp &rhs, const FieldPerp &x0) { - Mesh *mesh = rhs.getMesh(); - FieldPerp x(mesh); // Result + ASSERT1(localmesh == rhs.getMesh() && localmesh == x0.getMesh()); + + FieldPerp x(localmesh); // Result x.allocate(); - Coordinates *coord = mesh->getCoordinates(location); + Coordinates *coord = localmesh->getCoordinates(location); int jy = rhs.getIndex(); // Get the Y index x.setIndex(jy); @@ -102,8 +103,8 @@ const FieldPerp LaplaceCyclic::solve(const FieldPerp &rhs, const FieldPerp &x0) // Get the width of the boundary // If the flags to assign that only one guard cell should be used is set - int inbndry = mesh->xstart, outbndry=mesh->xstart; - if((global_flags & INVERT_BOTH_BNDRY_ONE) || (mesh->xstart < 2)) { + int inbndry = localmesh->xstart, outbndry=localmesh->xstart; + if((global_flags & INVERT_BOTH_BNDRY_ONE) || (localmesh->xstart < 2)) { inbndry = outbndry = 1; } if(inner_boundary_flags & INVERT_BNDRY_ONE) @@ -115,7 +116,7 @@ const FieldPerp LaplaceCyclic::solve(const FieldPerp &rhs, const FieldPerp &x0) BOUT_OMP(parallel) { /// Create a local thread-scope working array auto k1d = - Array(mesh->LocalNz); // ZFFT routine expects input of this length + Array(localmesh->LocalNz); // ZFFT routine expects input of this length // Loop over X indices, including boundaries but not guard cells. (unless periodic // in x) @@ -123,13 +124,13 @@ const FieldPerp LaplaceCyclic::solve(const FieldPerp &rhs, const FieldPerp &x0) for (int ix = xs; ix <= xe; ix++) { // Take DST in Z direction and put result in k1d - if (((ix < inbndry) && (inner_boundary_flags & INVERT_SET) && mesh->firstX()) || + if (((ix < inbndry) && (inner_boundary_flags & INVERT_SET) && localmesh->firstX()) || ((xe - ix < outbndry) && (outer_boundary_flags & INVERT_SET) && - mesh->lastX())) { + localmesh->lastX())) { // Use the values in x0 in the boundary - DST(x0[ix] + 1, mesh->LocalNz - 2, std::begin(k1d)); + DST(x0[ix] + 1, localmesh->LocalNz - 2, std::begin(k1d)); } else { - DST(rhs[ix] + 1, mesh->LocalNz - 2, std::begin(k1d)); + DST(rhs[ix] + 1, localmesh->LocalNz - 2, std::begin(k1d)); } // Copy into array, transposing so kz is first index @@ -141,7 +142,7 @@ const FieldPerp LaplaceCyclic::solve(const FieldPerp &rhs, const FieldPerp &x0) // including boundary conditions BOUT_OMP(for nowait) for (int kz = 0; kz < nmode; kz++) { - BoutReal zlen = coord->dz * (mesh->LocalNz - 3); + BoutReal zlen = coord->dz * (localmesh->LocalNz - 3); BoutReal kwave = kz * 2.0 * PI / (2. * zlen); // wave number is 1/[rad]; DST has extra 2. @@ -162,27 +163,27 @@ const FieldPerp LaplaceCyclic::solve(const FieldPerp &rhs, const FieldPerp &x0) BOUT_OMP(parallel) { /// Create a local thread-scope working array auto k1d = - Array(mesh->LocalNz); // ZFFT routine expects input of this length + Array(localmesh->LocalNz); // ZFFT routine expects input of this length BOUT_OMP(for nowait) for (int ix = xs; ix <= xe; ix++) { for (int kz = 0; kz < nmode; kz++) k1d[kz] = xcmplx(kz, ix - xs); - for (int kz = nmode; kz < (mesh->LocalNz); kz++) + for (int kz = nmode; kz < (localmesh->LocalNz); kz++) k1d[kz] = 0.0; // Filtering out all higher harmonics - DST_rev(std::begin(k1d), mesh->LocalNz - 2, x[ix] + 1); + DST_rev(std::begin(k1d), localmesh->LocalNz - 2, x[ix] + 1); x(ix, 0) = -x(ix, 2); - x(ix, mesh->LocalNz - 1) = -x(ix, mesh->LocalNz - 3); + x(ix, localmesh->LocalNz - 1) = -x(ix, localmesh->LocalNz - 3); } } }else { BOUT_OMP(parallel) { /// Create a local thread-scope working array - auto k1d = Array((mesh->LocalNz) / 2 + + auto k1d = Array((localmesh->LocalNz) / 2 + 1); // ZFFT routine expects input of this length // Loop over X indices, including boundaries but not guard cells (unless periodic in @@ -191,13 +192,13 @@ const FieldPerp LaplaceCyclic::solve(const FieldPerp &rhs, const FieldPerp &x0) for (int ix = xs; ix <= xe; ix++) { // Take FFT in Z direction, apply shift, and put result in k1d - if (((ix < inbndry) && (inner_boundary_flags & INVERT_SET) && mesh->firstX()) || + if (((ix < inbndry) && (inner_boundary_flags & INVERT_SET) && localmesh->firstX()) || ((xe - ix < outbndry) && (outer_boundary_flags & INVERT_SET) && - mesh->lastX())) { + localmesh->lastX())) { // Use the values in x0 in the boundary - rfft(x0[ix], mesh->LocalNz, std::begin(k1d)); + rfft(x0[ix], localmesh->LocalNz, std::begin(k1d)); } else { - rfft(rhs[ix], mesh->LocalNz, std::begin(k1d)); + rfft(rhs[ix], localmesh->LocalNz, std::begin(k1d)); } // Copy into array, transposing so kz is first index @@ -227,7 +228,7 @@ const FieldPerp LaplaceCyclic::solve(const FieldPerp &rhs, const FieldPerp &x0) BOUT_OMP(parallel) { /// Create a local thread-scope working array - auto k1d = Array((mesh->LocalNz) / 2 + + auto k1d = Array((localmesh->LocalNz) / 2 + 1); // ZFFT routine expects input of this length BOUT_OMP(for nowait) @@ -235,10 +236,10 @@ const FieldPerp LaplaceCyclic::solve(const FieldPerp &rhs, const FieldPerp &x0) for (int kz = 0; kz < nmode; kz++) k1d[kz] = xcmplx(kz, ix - xs); - for (int kz = nmode; kz < (mesh->LocalNz) / 2 + 1; kz++) + for (int kz = nmode; kz < (localmesh->LocalNz) / 2 + 1; kz++) k1d[kz] = 0.0; // Filtering out all higher harmonics - irfft(std::begin(k1d), mesh->LocalNz, x[ix]); + irfft(std::begin(k1d), localmesh->LocalNz, x[ix]); } } } @@ -250,11 +251,11 @@ const Field3D LaplaceCyclic::solve(const Field3D &rhs, const Field3D &x0) { ASSERT1(rhs.getLocation() == location); ASSERT1(x0.getLocation() == location); + ASSERT1(localmesh == rhs.getMesh() && localmesh == x0.getMesh()); Timer timer("invert"); - Mesh *mesh = rhs.getMesh(); - Field3D x(mesh); // Result + Field3D x(localmesh); // Result x.allocate(); x.setLocation(location); @@ -263,8 +264,8 @@ const Field3D LaplaceCyclic::solve(const Field3D &rhs, const Field3D &x0) { // Get the width of the boundary // If the flags to assign that only one guard cell should be used is set - int inbndry = mesh->xstart, outbndry = mesh->xstart; - if ((global_flags & INVERT_BOTH_BNDRY_ONE) || (mesh->xstart < 2)) { + int inbndry = localmesh->xstart, outbndry = localmesh->xstart; + if ((global_flags & INVERT_BOTH_BNDRY_ONE) || (localmesh->xstart < 2)) { inbndry = outbndry = 1; } if (inner_boundary_flags & INVERT_BNDRY_ONE) @@ -275,17 +276,17 @@ const Field3D LaplaceCyclic::solve(const Field3D &rhs, const Field3D &x0) { int nx = xe - xs + 1; // Number of X points on this processor // Get range of Y indices - int ys = mesh->ystart, ye = mesh->yend; + int ys = localmesh->ystart, ye = localmesh->yend; - if (mesh->hasBndryLowerY()) { + if (localmesh->hasBndryLowerY()) { if (include_yguards) ys = 0; // Mesh contains a lower boundary and we are solving in the guard cells ys += extra_yguards_lower; } - if (mesh->hasBndryUpperY()) { + if (localmesh->hasBndryUpperY()) { if (include_yguards) - ye = mesh->LocalNy - + ye = localmesh->LocalNy - 1; // Contains upper boundary and we are solving in the guard cells ye -= extra_yguards_upper; @@ -306,7 +307,7 @@ const Field3D LaplaceCyclic::solve(const Field3D &rhs, const Field3D &x0) { BOUT_OMP(parallel) { /// Create a local thread-scope working array auto k1d = - Array(mesh->LocalNz); // ZFFT routine expects input of this length + Array(localmesh->LocalNz); // ZFFT routine expects input of this length // Loop over X and Y indices, including boundaries but not guard cells. // (unless periodic in x) @@ -318,13 +319,13 @@ const Field3D LaplaceCyclic::solve(const Field3D &rhs, const Field3D &x0) { // Take DST in Z direction and put result in k1d - if (((ix < inbndry) && (inner_boundary_flags & INVERT_SET) && mesh->firstX()) || + if (((ix < inbndry) && (inner_boundary_flags & INVERT_SET) && localmesh->firstX()) || ((xe - ix < outbndry) && (outer_boundary_flags & INVERT_SET) && - mesh->lastX())) { + localmesh->lastX())) { // Use the values in x0 in the boundary - DST(x0(ix, iy) + 1, mesh->LocalNz - 2, std::begin(k1d)); + DST(x0(ix, iy) + 1, localmesh->LocalNz - 2, std::begin(k1d)); } else { - DST(rhs(ix, iy) + 1, mesh->LocalNz - 2, std::begin(k1d)); + DST(rhs(ix, iy) + 1, localmesh->LocalNz - 2, std::begin(k1d)); } // Copy into array, transposing so kz is first index @@ -340,7 +341,7 @@ const Field3D LaplaceCyclic::solve(const Field3D &rhs, const Field3D &x0) { int iy = ys + ind / nmode; int kz = ind % nmode; - BoutReal zlen = coord->dz * (mesh->LocalNz - 3); + BoutReal zlen = coord->dz * (localmesh->LocalNz - 3); BoutReal kwave = kz * 2.0 * PI / (2. * zlen); // wave number is 1/[rad]; DST has extra 2. @@ -361,7 +362,7 @@ const Field3D LaplaceCyclic::solve(const Field3D &rhs, const Field3D &x0) { BOUT_OMP(parallel) { /// Create a local thread-scope working array auto k1d = - Array(mesh->LocalNz); // ZFFT routine expects input of this length + Array(localmesh->LocalNz); // ZFFT routine expects input of this length BOUT_OMP(for nowait) for (int ind = 0; ind < nxny; ++ind) { // Loop over X and Y @@ -372,19 +373,19 @@ const Field3D LaplaceCyclic::solve(const Field3D &rhs, const Field3D &x0) { for (int kz = 0; kz < nmode; kz++) k1d[kz] = xcmplx3D((iy - ys) * nmode + kz, ix - xs); - for (int kz = nmode; kz < mesh->LocalNz; kz++) + for (int kz = nmode; kz < localmesh->LocalNz; kz++) k1d[kz] = 0.0; // Filtering out all higher harmonics - DST_rev(std::begin(k1d), mesh->LocalNz - 2, &x(ix, iy, 1)); + DST_rev(std::begin(k1d), localmesh->LocalNz - 2, &x(ix, iy, 1)); x(ix, iy, 0) = -x(ix, iy, 2); - x(ix, iy, mesh->LocalNz - 1) = -x(ix, iy, mesh->LocalNz - 3); + x(ix, iy, localmesh->LocalNz - 1) = -x(ix, iy, localmesh->LocalNz - 3); } } } else { BOUT_OMP(parallel) { /// Create a local thread-scope working array - auto k1d = Array(mesh->LocalNz / 2 + + auto k1d = Array(localmesh->LocalNz / 2 + 1); // ZFFT routine expects input of this length // Loop over X and Y indices, including boundaries but not guard cells @@ -398,13 +399,13 @@ const Field3D LaplaceCyclic::solve(const Field3D &rhs, const Field3D &x0) { // Take FFT in Z direction, apply shift, and put result in k1d - if (((ix < inbndry) && (inner_boundary_flags & INVERT_SET) && mesh->firstX()) || + if (((ix < inbndry) && (inner_boundary_flags & INVERT_SET) && localmesh->firstX()) || ((xe - ix < outbndry) && (outer_boundary_flags & INVERT_SET) && - mesh->lastX())) { + localmesh->lastX())) { // Use the values in x0 in the boundary - rfft(x0(ix, iy), mesh->LocalNz, std::begin(k1d)); + rfft(x0(ix, iy), localmesh->LocalNz, std::begin(k1d)); } else { - rfft(rhs(ix, iy), mesh->LocalNz, std::begin(k1d)); + rfft(rhs(ix, iy), localmesh->LocalNz, std::begin(k1d)); } // Copy into array, transposing so kz is first index @@ -437,7 +438,7 @@ const Field3D LaplaceCyclic::solve(const Field3D &rhs, const Field3D &x0) { // FFT back to real space BOUT_OMP(parallel) { /// Create a local thread-scope working array - auto k1d = Array((mesh->LocalNz) / 2 + + auto k1d = Array((localmesh->LocalNz) / 2 + 1); // ZFFT routine expects input of this length BOUT_OMP(for nowait) @@ -449,10 +450,10 @@ const Field3D LaplaceCyclic::solve(const Field3D &rhs, const Field3D &x0) { for (int kz = 0; kz < nmode; kz++) k1d[kz] = xcmplx3D((iy - ys) * nmode + kz, ix - xs); - for (int kz = nmode; kz < mesh->LocalNz / 2 + 1; kz++) + for (int kz = nmode; kz < localmesh->LocalNz / 2 + 1; kz++) k1d[kz] = 0.0; // Filtering out all higher harmonics - irfft(std::begin(k1d), mesh->LocalNz, x(ix, iy)); + irfft(std::begin(k1d), localmesh->LocalNz, x(ix, iy)); } } } diff --git a/src/invert/laplace/impls/cyclic/cyclic_laplace.hxx b/src/invert/laplace/impls/cyclic/cyclic_laplace.hxx index bb1d9f31c2..b007542569 100644 --- a/src/invert/laplace/impls/cyclic/cyclic_laplace.hxx +++ b/src/invert/laplace/impls/cyclic/cyclic_laplace.hxx @@ -44,22 +44,25 @@ class LaplaceCyclic; */ class LaplaceCyclic : public Laplacian { public: - LaplaceCyclic(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE); + LaplaceCyclic(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = mesh); ~LaplaceCyclic(); using Laplacian::setCoefA; void setCoefA(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh = val.getMesh()); Acoef = val; } using Laplacian::setCoefC; void setCoefC(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh = val.getMesh()); Ccoef = val; } using Laplacian::setCoefD; void setCoefD(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh = val.getMesh()); Dcoef = val; } using Laplacian::setCoefEx; diff --git a/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx b/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx index 17e0e6b127..567ec0e583 100644 --- a/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx +++ b/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx @@ -37,8 +37,8 @@ BoutReal soltime=0.0,settime=0.0; -LaplaceMultigrid::LaplaceMultigrid(Options *opt, const CELL_LOC loc) : - Laplacian(opt, loc), +LaplaceMultigrid::LaplaceMultigrid(Options *opt, const CELL_LOC loc, Mesh *mesh_in) : + Laplacian(opt, loc, mesh_in), A(0.0), C1(1.0), C2(1.0), D(1.0) { TRACE("LaplaceMultigrid::LaplaceMultigrid(Options *opt)"); @@ -84,20 +84,20 @@ LaplaceMultigrid::LaplaceMultigrid(Options *opt, const CELL_LOC loc) : throw BoutException("Attempted to set Laplacian outer boundary inversion flag that is not implemented in LaplaceMultigrid."); } - commX = mesh->getXcomm(); + commX = localmesh->getXcomm(); - Nx_local = mesh->xend - mesh->xstart + 1; // excluding guard cells - Nx_global = mesh->GlobalNx - 2*mesh->xstart; // excluding guard cells + 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="<GlobalNz; + Nz_global = localmesh->GlobalNz; Nz_local = Nz_global; // No parallelization in z-direction (for now) // //else { - // Nz_local = mesh->zend - mesh->zstart + 1; // excluding guard cells - // Nz_global = mesh->GlobalNz - 2*mesh->zstart; // excluding guard cells + // Nz_local = localmesh->zend - localmesh->zstart + 1; // excluding guard cells + // Nz_global = localmesh->GlobalNz - 2*localmesh->zstart; // excluding guard cells // } if (mgcount==0) { output <<"Nz="<getCoordinates(location); + Coordinates *coords = localmesh->getCoordinates(location); yindex = b_in.getIndex(); int level = kMG->mglevel-1; @@ -231,7 +232,7 @@ BOUT_OMP(parallel default(shared) ) BOUT_OMP(for collapse(2)) for (int i=1; ixstart; + int i2 = i-1+localmesh->xstart; int k2 = k-1; x[i*lz2+k] = x0[i2][k2]; } @@ -243,13 +244,13 @@ BOUT_OMP(parallel default(shared) ) BOUT_OMP(for collapse(2)) for (int i=1; ixstart; + int i2 = i-1+localmesh->xstart; int k2 = k-1; b[i*lz2+k] = b_in(i2, k2); } } - if (mesh->firstX()) { + if (localmesh->firstX()) { if ( inner_boundary_flags & INVERT_AC_GRAD ) { // Neumann boundary condition if ( inner_boundary_flags & INVERT_SET ) { @@ -258,7 +259,7 @@ BOUT_OMP(parallel default(shared) ) BOUT_OMP(for) for (int k=1; kxstart-1, k2)*sqrt(coords->g_11(mesh->xstart, yindex))*coords->dx(mesh->xstart, yindex); + 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 @@ -277,7 +278,7 @@ BOUT_OMP(parallel default(shared) ) BOUT_OMP(for) for (int k=1; kxstart-1, k2); + x[k] = 2.*x0(localmesh->xstart-1, k2); // this is the value to set at the inner boundary } } @@ -292,7 +293,7 @@ BOUT_OMP(for) } } } - if (mesh->lastX()) { + if (localmesh->lastX()) { if ( outer_boundary_flags & INVERT_AC_GRAD ) { // Neumann boundary condition if ( inner_boundary_flags & INVERT_SET ) { @@ -301,7 +302,7 @@ BOUT_OMP(parallel default(shared) ) BOUT_OMP(for) for (int k=1; kxend+1, k2)*sqrt(coords->g_11(mesh->xend, yindex))*coords->dx(mesh->xend, yindex); + 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 } } @@ -323,7 +324,7 @@ BOUT_OMP(parallel default(shared) ) BOUT_OMP(for) for (int k=1; kxend+1, k2); + x[(lxx+1)*lz2+k]=2.*x0(localmesh->xend+1, k2); // this is the value to set at the outer boundary } } @@ -423,13 +424,13 @@ BOUT_OMP(for) } } - FieldPerp result(mesh); + FieldPerp result(localmesh); result.allocate(); result.setIndex(yindex); #if CHECK>2 // Make any unused elements NaN so that user does not try to do calculations with them - const auto ®ion = mesh->getRegionPerp("RGN_ALL"); + const auto ®ion = localmesh->getRegionPerp("RGN_ALL"); BOUT_FOR(i, region) { result[i] = BoutNaN; } @@ -439,27 +440,27 @@ BOUT_OMP(parallel default(shared) ) BOUT_OMP(for collapse(2)) for (int i=1; ixstart; + int i2 = i-1+localmesh->xstart; int k2 = k-1; result(i2, k2) = x[i*lz2+k]; } } - if (mesh->firstX()) { + if (localmesh->firstX()) { if ( inner_boundary_flags & INVERT_AC_GRAD ) { // Neumann boundary condition if ( inner_boundary_flags & INVERT_SET ) { // guard cells of x0 specify gradient to set at inner boundary - int i2 = -1+mesh->xstart; + int i2 = -1+localmesh->xstart; BOUT_OMP(parallel default(shared) ) BOUT_OMP(for) for (int k=1; kxstart-1, k2)*sqrt(coords->g_11(mesh->xstart, yindex))*coords->dx(mesh->xstart, yindex); + 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+mesh->xstart; + int i2 = -1+localmesh->xstart; BOUT_OMP(parallel default(shared) ) BOUT_OMP(for) for (int k=1; kxstart; + int i2 = -1+localmesh->xstart; BOUT_OMP(parallel default(shared) ) BOUT_OMP(for) for (int k=1; kxstart-1,k2) - x[lz2+k]; + result(i2, k2) = 2.*x0(localmesh->xstart-1,k2) - x[lz2+k]; } } else { // zero value inner boundary condition - int i2 = -1+mesh->xstart; + int i2 = -1+localmesh->xstart; BOUT_OMP(parallel default(shared) ) BOUT_OMP(for) for (int k=1; klastX()) { + if (localmesh->lastX()) { if ( outer_boundary_flags & INVERT_AC_GRAD ) { // Neumann boundary condition if ( inner_boundary_flags & INVERT_SET ) { // guard cells of x0 specify gradient to set at outer boundary - int i2 = lxx+mesh->xstart; + int i2 = lxx+localmesh->xstart; BOUT_OMP(parallel default(shared) ) BOUT_OMP(for) for (int k=1; kxend+1, k2)*sqrt(coords->g_11(mesh->xend, yindex))*coords->dx(mesh->xend, yindex); + 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+mesh->xstart; + int i2 = lxx+localmesh->xstart; BOUT_OMP(parallel default(shared) ) BOUT_OMP(for) for (int k=1; kxstart; + int i2 = lxx+localmesh->xstart; BOUT_OMP(parallel default(shared) ) BOUT_OMP(for) for (int k=1; kxend+1,k2) - x[lxx*lz2+k]; + result(i2, k2) = 2.*x0(localmesh->xend+1,k2) - x[lxx*lz2+k]; } } else { // zero value inner boundary condition - int i2 = lxx+mesh->xstart; + int i2 = lxx+localmesh->xstart; BOUT_OMP(parallel default(shared) ) BOUT_OMP(for) for (int k=1; kgetCoordinates(location); + Coordinates *coords = localmesh->getCoordinates(location); BoutReal *mat; mat = kMG->matmg[level]; int llx = kMG->lnx[level]; @@ -564,7 +565,7 @@ BOUT_OMP(parallel default(shared)) BOUT_OMP(for collapse(2)) for (int i=1; ixstart; + int i2 = i-1+localmesh->xstart; int k2 = k-1; int k2p = (k2+1)%Nz_global; int k2m = (k2+Nz_global-1)%Nz_global; diff --git a/src/invert/laplace/impls/multigrid/multigrid_laplace.hxx b/src/invert/laplace/impls/multigrid/multigrid_laplace.hxx index 6c22d8105a..00eaf96a6b 100644 --- a/src/invert/laplace/impls/multigrid/multigrid_laplace.hxx +++ b/src/invert/laplace/impls/multigrid/multigrid_laplace.hxx @@ -132,28 +132,33 @@ private: class LaplaceMultigrid : public Laplacian { public: - LaplaceMultigrid(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE); + LaplaceMultigrid(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = mesh); ~LaplaceMultigrid() {}; 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"); } @@ -161,28 +166,35 @@ public: 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; } const FieldPerp solve(const FieldPerp &b) override { - FieldPerp zero(b.getMesh()); + ASSERT1(localmesh == b.getMesh()); + + FieldPerp zero(localmesh); zero = 0.; zero.setIndex(b.getIndex()); return solve(b, zero); diff --git a/src/invert/laplace/impls/mumps/mumps_laplace.cxx b/src/invert/laplace/impls/mumps/mumps_laplace.cxx index df8967babc..dd30e4b15e 100644 --- a/src/invert/laplace/impls/mumps/mumps_laplace.cxx +++ b/src/invert/laplace/impls/mumps/mumps_laplace.cxx @@ -34,8 +34,8 @@ #include #include -LaplaceMumps::LaplaceMumps(Options *opt, const CELL_LOC loc) : - Laplacian(opt, loc), +LaplaceMumps::LaplaceMumps(Options *opt, const CELL_LOC loc, Mesh *mesh_in = mesh) : + Laplacian(opt, loc, mesh_in), A(0.0), C1(1.0), C2(1.0), D(1.0), Ex(0.0), Ez(0.0), issetD(false), issetC(false), issetE(false) { @@ -65,21 +65,21 @@ LaplaceMumps::LaplaceMumps(Options *opt, const CELL_LOC loc) : if (outer_boundary_flags & ~implemented_boundary_flags) { throw BoutException("Attempted to set Laplacian inversion boundary condition flag that is not implemented in mumps_laplace.cxx"); } - if(mesh->periodicX) { - throw BoutException("LaplaceMumps does not work with periodicity in the x direction (mesh->PeriodicX == true). Change boundary conditions or use serial-tri or cyclic solver instead"); + if(localmesh->periodicX) { + throw BoutException("LaplaceMumps does not work with periodicity in the x direction (localmesh->PeriodicX == true). Change boundary conditions or use serial-tri or cyclic solver instead"); } #endif // Get communicator for group of processors in X - all points in z-x plane for fixed y. - comm = mesh->getXcomm(); + comm = localmesh->getXcomm(); // Need to determine local size to use based on prior parallelisation // Coefficient values are stored only on local processors. - localN = (mesh->xend - mesh->xstart + 1) * (mesh->LocalNz); - if(mesh->firstX()) - localN += mesh->xstart * (mesh->LocalNz); // If on first processor add on width of boundary region - if(mesh->lastX()) - localN += mesh->xstart * (mesh->LocalNz); // If on last processor add on width of boundary region + localN = (localmesh->xend - localmesh->xstart + 1) * (localmesh->LocalNz); + if(localmesh->firstX()) + localN += localmesh->xstart * (localmesh->LocalNz); // If on first processor add on width of boundary region + if(localmesh->lastX()) + localN += localmesh->xstart * (localmesh->LocalNz); // If on last processor add on width of boundary region // Calculate total number of points in physical grid @@ -87,11 +87,11 @@ LaplaceMumps::LaplaceMumps(Options *opt, const CELL_LOC loc) : throw BoutException("Error in MPI_Allreduce during LaplacePetsc initialisation"); // Calculate total (physical) grid dimensions - meshz = mesh->GlobalNz-1; + meshz = localmesh->GlobalNz-1; meshx = size / meshz; // Calculate number of guard cells in x-direction - nxguards = mesh->LocalNx - (mesh->xend-mesh->xstart+1); + nxguards = localmesh->LocalNx - (localmesh->xend-localmesh->xstart+1); // Get implementation specific options opts->get("fourth_order", fourth_order, false); @@ -99,7 +99,7 @@ LaplaceMumps::LaplaceMumps(Options *opt, const CELL_LOC loc) : sol.allocate(); - mumps_struc.comm_fortran = (MUMPS_INT) MPI_Comm_c2f(mesh->getXcomm()); // MPI communicator for MUMPS, in fortran format + mumps_struc.comm_fortran = (MUMPS_INT) MPI_Comm_c2f(localmesh->getXcomm()); // MPI communicator for MUMPS, in fortran format mumps_struc.sym = 0; // Solve using unsymmetric matrix mumps_struc.par = 1; // Use the host processor (rank 0) to do work for the solution @@ -112,39 +112,39 @@ LaplaceMumps::LaplaceMumps(Options *opt, const CELL_LOC loc) : // nz is the total number of non-zero elements in the matrix, nz_loc is the number of non-zero elements on this processor if (fourth_order) { mumps_struc.nz = 25*(meshx-nxguards)*meshz; - mumps_struc.nz_loc = 25*(mesh->xend-mesh->xstart+1)*(mesh->LocalNz); + mumps_struc.nz_loc = 25*(localmesh->xend-localmesh->xstart+1)*(localmesh->LocalNz); } else { mumps_struc.nz = 9*(meshx-nxguards)*meshz; - mumps_struc.nz_loc = 9*(mesh->xend-mesh->xstart+1)*(mesh->LocalNz); + mumps_struc.nz_loc = 9*(localmesh->xend-localmesh->xstart+1)*(localmesh->LocalNz); } if (inner_boundary_flags & INVERT_AC_GRAD) { if (fourth_order) { - mumps_struc.nz += 5*meshz*mesh->xstart; - if (mesh->firstX()) mumps_struc.nz_loc += 5*mesh->xstart*(mesh->LocalNz); + mumps_struc.nz += 5*meshz*localmesh->xstart; + if (localmesh->firstX()) mumps_struc.nz_loc += 5*localmesh->xstart*(localmesh->LocalNz); } else { - mumps_struc.nz += 3*meshz*(mesh->xstart); - if (mesh->firstX()) mumps_struc.nz_loc += 3*mesh->xstart*(mesh->LocalNz); + mumps_struc.nz += 3*meshz*(localmesh->xstart); + if (localmesh->firstX()) mumps_struc.nz_loc += 3*localmesh->xstart*(localmesh->LocalNz); } } else { - mumps_struc.nz += mesh->xstart*meshz; - if (mesh->firstX()) mumps_struc.nz_loc += mesh->xstart*(mesh->LocalNz); + mumps_struc.nz += localmesh->xstart*meshz; + if (localmesh->firstX()) mumps_struc.nz_loc += localmesh->xstart*(localmesh->LocalNz); } if (outer_boundary_flags & INVERT_AC_GRAD) { if (fourth_order) { - mumps_struc.nz += 5*(mesh->LocalNx-mesh->xend-1)*meshz; - if (mesh->lastX()) mumps_struc.nz_loc += 5*(mesh->LocalNx-mesh->xend-1)*(mesh->LocalNz); + mumps_struc.nz += 5*(localmesh->LocalNx-localmesh->xend-1)*meshz; + if (localmesh->lastX()) mumps_struc.nz_loc += 5*(localmesh->LocalNx-localmesh->xend-1)*(localmesh->LocalNz); } else { - mumps_struc.nz += 3*(mesh->LocalNx-mesh->xend-1)*meshz; - if (mesh->lastX()) mumps_struc.nz_loc += 3*(mesh->LocalNx-mesh->xend-1)*(mesh->LocalNz); + mumps_struc.nz += 3*(localmesh->LocalNx-localmesh->xend-1)*meshz; + if (localmesh->lastX()) mumps_struc.nz_loc += 3*(localmesh->LocalNx-localmesh->xend-1)*(localmesh->LocalNz); } } else { - mumps_struc.nz += (mesh->LocalNx-mesh->xend-1)*meshz; - if (mesh->lastX()) mumps_struc.nz_loc += (mesh->LocalNx-mesh->xend-1)*(mesh->LocalNz); + mumps_struc.nz += (localmesh->LocalNx-localmesh->xend-1)*meshz; + if (localmesh->lastX()) mumps_struc.nz_loc += (localmesh->LocalNx-localmesh->xend-1)*(localmesh->LocalNz); } // // These would be needed if giving the matrix only on the host processor, or possibly if providing the structure on the host processor for analysis // mumps_struc.irn = new MUMPS_INT[mumps_struc.nz]; @@ -154,13 +154,13 @@ LaplaceMumps::LaplaceMumps(Options *opt, const CELL_LOC loc) : mumps_struc.jcn_loc = new MUMPS_INT[mumps_struc.nz_loc]; // list of GLOBAL column indices of local matrix entries mumps_struc.a_loc = new BoutReal[mumps_struc.nz_loc]; // the matrix entries - if (mesh->firstX()) { + if (localmesh->firstX()) { mumps_struc.nrhs = 1; // number of right hand side vectors mumps_struc.lrhs = mumps_struc.n; // leading dimension of rhs (i.e. length of vector) // mumps_struc.rhs = new BoutReal[mumps_struc.lrhs*mumps_struc.nrhs]; // rhs, to be provided on the rank-0 processor only } -// if (mesh->firstX()) mumps_struc.sol_loc = *sol.getData(); // pointer to the array to put the solution in, starts at 0 on first processor -// else mumps_struc.sol_loc = *sol.getData() + mesh->xstart*meshz; // pointer to the array to put the solution in, starts at mesh->xstart +// if (localmesh->firstX()) mumps_struc.sol_loc = *sol.getData(); // pointer to the array to put the solution in, starts at 0 on first processor +// else mumps_struc.sol_loc = *sol.getData() + localmesh->xstart*meshz; // pointer to the array to put the solution in, starts at localmesh->xstart // mumps_struc.lsol_loc = localN; // size of the (local) solution array // mumps_struc.isol_loc = new MUMPS_INT[localN]; // list of indices of the solution array (though this is all local points) mumps_struc.icntl[2] = 0; // Suppress output of global information @@ -180,43 +180,43 @@ LaplaceMumps::LaplaceMumps(Options *opt, const CELL_LOC loc) : // mumps_struc.job = MUMPS_JOB_ALL; // iteration_count = repeat_analysis; -// localrhssize = (mesh->xend-mesh->xstart+1)*mesh->LocalNy*mesh->LocalNz; -// if (mesh->lastX()) { -// localrhssize += (mesh->LocalNx-mesh->xend-1)*mesh->LocalNy*mesh->LocalNz; +// localrhssize = (localmesh->xend-localmesh->xstart+1)*localmesh->LocalNy*localmesh->LocalNz; +// if (localmesh->lastX()) { +// localrhssize += (localmesh->LocalNx-localmesh->xend-1)*localmesh->LocalNy*localmesh->LocalNz; // } -// if (mesh->firstX()) { -// localrhssize += mesh->xstart*mesh->LocalNy*mesh->LocalNz; +// if (localmesh->firstX()) { +// localrhssize += localmesh->xstart*localmesh->LocalNy*localmesh->LocalNz; // -// int nxpe = mesh->NXPE; +// int nxpe = localmesh->NXPE; // localrhs_size_array = new int[nxpe]; // localrhs_size_array[0] = localrhssize; // if (nxpe>1) { // for (int i=1; ixend-mesh->xstart+1)*mesh->LocalNy*mesh->LocalNz; -// localrhs_size_array[nxpe-1] = (mesh->LocalNx-mesh->xstart)*mesh->LocalNy*mesh->LocalNz; +// localrhs_size_array[i] = (localmesh->xend-localmesh->xstart+1)*localmesh->LocalNy*localmesh->LocalNz; +// localrhs_size_array[nxpe-1] = (localmesh->LocalNx-localmesh->xstart)*localmesh->LocalNy*localmesh->LocalNz; // } // rhs_positions = new int[nxpe]; // rhs_positions[0] = 0; // for (int i=1; iLocalNy*mesh->LocalNz]; -// rhs_slice = new BoutReal[meshx*mesh->LocalNz]; +// rhs = new BoutReal[meshx*localmesh->LocalNy*localmesh->LocalNz]; +// rhs_slice = new BoutReal[meshx*localmesh->LocalNz]; // } - localrhssize = (mesh->xend-mesh->xstart+1)*(mesh->LocalNz); - if (mesh->lastX()) { - localrhssize += (mesh->LocalNx-mesh->xend-1)*(mesh->LocalNz); + localrhssize = (localmesh->xend-localmesh->xstart+1)*(localmesh->LocalNz); + if (localmesh->lastX()) { + localrhssize += (localmesh->LocalNx-localmesh->xend-1)*(localmesh->LocalNz); } - if (mesh->firstX()) { - localrhssize += mesh->xstart*(mesh->LocalNz); + if (localmesh->firstX()) { + localrhssize += localmesh->xstart*(localmesh->LocalNz); - int nxpe = mesh->NXPE; + int nxpe = localmesh->NXPE; localrhs_size_array = Array(nxpe); localrhs_size_array[0] = localrhssize; if (nxpe>1) { for (int i=1; ixend-mesh->xstart+1)*(mesh->LocalNz); - localrhs_size_array[nxpe-1] = (mesh->LocalNx-mesh->xstart)*(mesh->LocalNz); + localrhs_size_array[i] = (localmesh->xend-localmesh->xstart+1)*(localmesh->LocalNz); + localrhs_size_array[nxpe-1] = (localmesh->LocalNx-localmesh->xstart)*(localmesh->LocalNz); } rhs_positions = Array(nxpe); rhs_positions[0] = 0; @@ -229,9 +229,9 @@ LaplaceMumps::LaplaceMumps(Options *opt, const CELL_LOC loc) : // Set Arrays of matrix indices, using i (0<=ifirstX()) - for (int x=0; xxstart; x++) - for (int z=0; zLocalNz; z++) { + if (localmesh->firstX()) + for (int x=0; xxstart; x++) + for (int z=0; zLocalNz; z++) { int x0 = x; int xp = x+1; int xpp = x+2; @@ -263,13 +263,13 @@ LaplaceMumps::LaplaceMumps(Options *opt, const CELL_LOC loc) : // mumps_struc.isol_loc[j] = x0*meshz + z0 + 1; // Indices for fortran arrays that start at 1 // j++; } - for (int x=mesh->xstart; x<=mesh->xend; x++) - for (int z=0; zLocalNz; z++) { - int xmm = mesh->XGLOBAL(x)-2; - int xm = mesh->XGLOBAL(x)-1; - int x0 = mesh->XGLOBAL(x); - int xp = mesh->XGLOBAL(x)+1; - int xpp = mesh->XGLOBAL(x)+2; + for (int x=localmesh->xstart; x<=localmesh->xend; x++) + for (int z=0; zLocalNz; z++) { + int xmm = localmesh->XGLOBAL(x)-2; + int xm = localmesh->XGLOBAL(x)-1; + int x0 = localmesh->XGLOBAL(x); + int xp = localmesh->XGLOBAL(x)+1; + int xpp = localmesh->XGLOBAL(x)+2; int zmm = (z-2<0) ? (z-2+meshz) : (z-2); int zm = (z-1<0) ? (z-1+meshz) : (z-1); int z0 = z; @@ -384,12 +384,12 @@ LaplaceMumps::LaplaceMumps(Options *opt, const CELL_LOC loc) : // mumps_struc.isol_loc[j] = x0*meshz + z0 + 1; // Indices for fortran arrays that start at 1 // j++; } - if (mesh->lastX()) - for (int x=mesh->xend+1; xLocalNx; x++) - for (int z=0; zLocalNz; z++) { - int xmm = mesh->XGLOBAL(mesh->xend)+x-mesh->xend-2; - int xm = mesh->XGLOBAL(mesh->xend)+x-mesh->xend-1; - int x0 = mesh->XGLOBAL(mesh->xend)+x-mesh->xend; + if (localmesh->lastX()) + for (int x=localmesh->xend+1; xLocalNx; x++) + for (int z=0; zLocalNz; z++) { + int xmm = localmesh->XGLOBAL(localmesh->xend)+x-localmesh->xend-2; + int xm = localmesh->XGLOBAL(localmesh->xend)+x-localmesh->xend-1; + int x0 = localmesh->XGLOBAL(localmesh->xend)+x-localmesh->xend; int z0 = z; if(outer_boundary_flags & INVERT_AC_GRAD) { mumps_struc.irn_loc[i] = x0*meshz + z0 + 1; // Indices for fortran arrays that start at 1 @@ -442,17 +442,17 @@ LaplaceMumps::LaplaceMumps(Options *opt, const CELL_LOC loc) : // #if CHECK > 0 // msg_stack.push("Laplacian::solve(Field3D)"); // #endif -// int ys = mesh->ystart, ye = mesh->yend; +// int ys = localmesh->ystart, ye = localmesh->yend; // -// if(mesh->hasBndryLowerY()) { +// if(localmesh->hasBndryLowerY()) { // if (include_yguards) // ys = 0; // Mesh contains a lower boundary and we are solving in the guard cells // // ys += extra_yguards_lower; // } -// if(mesh->hasBndryUpperY()) { +// if(localmesh->hasBndryUpperY()) { // if (include_yguards) -// ye = mesh->LocalNy-1; // Contains upper boundary and we are solving in the guard cells +// ye = localmesh->LocalNy-1; // Contains upper boundary and we are solving in the guard cells // // ye -= extra_yguards_upper; // } @@ -460,43 +460,43 @@ LaplaceMumps::LaplaceMumps(Options *opt, const CELL_LOC loc) : // Field3D x = copy(b); // Force new memory allocation as we will mess around with x's data via pointers (i.e. 'unsafely') // // BoutReal* localrhs = **x.getData(); // Input the rhs in the solution field as solution will be returned in place by MUMPS -// if (!mesh->firstX()) localrhs += mesh->xstart*mesh->LocalNy*mesh->LocalNz; -// MPI_Gatherv(localrhs,localrhssize,MPI_DOUBLE,rhs,localrhs_size_array,rhs_positions,MPI_DOUBLE,0,mesh->getXcomm()); +// if (!localmesh->firstX()) localrhs += localmesh->xstart*localmesh->LocalNy*localmesh->LocalNz; +// MPI_Gatherv(localrhs,localrhssize,MPI_DOUBLE,rhs,localrhs_size_array,rhs_positions,MPI_DOUBLE,0,localmesh->getXcomm()); // // if ( ++iteration_count > repeat_analysis ) { // mumps_struc.job = MUMPS_JOB_ALL; // for(int jy=ys; jy <= ye; jy++) { -// if (mesh->firstX()) +// if (localmesh->firstX()) // for(int jx=0; jxLocalNy*mesh->LocalNz + jy*mesh->LocalNz + jz]; +// rhs_slice[jx*meshz+jz] = rhs[jx*localmesh->LocalNy*localmesh->LocalNz + jy*localmesh->LocalNz + jz]; // // solve(rhs_slice,jy); // -// if (mesh->firstX()) +// if (localmesh->firstX()) // for(int jx=0; jxLocalNy*mesh->LocalNz + jy*mesh->LocalNz + jz] = rhs_slice[jx*meshz+jz]; +// rhs[jx*localmesh->LocalNy*localmesh->LocalNz + jy*localmesh->LocalNz + jz] = rhs_slice[jx*meshz+jz]; // } // mumps_struc.job = MUMPS_JOB_BOTH; // } // else { // for(int jy=ys; jy <= ye; jy++) { -// if (mesh->firstX()) +// if (localmesh->firstX()) // for(int jx=0; jxLocalNy*mesh->LocalNz + jy*mesh->LocalNz + jz]; +// rhs_slice[jx*meshz+jz] = rhs[jx*localmesh->LocalNy*localmesh->LocalNz + jy*localmesh->LocalNz + jz]; // // solve(rhs_slice,jy); // -// if (mesh->firstX()) +// if (localmesh->firstX()) // for(int jx=0; jxLocalNy*mesh->LocalNz + jy*mesh->LocalNz + jz] = rhs_slice[jx*meshz+jz]; +// rhs[jx*localmesh->LocalNy*localmesh->LocalNz + jy*localmesh->LocalNz + jz] = rhs_slice[jx*meshz+jz]; // } // } // -// MPI_Scatterv(rhs,localrhs_size_array,rhs_positions,MPI_DOUBLE,localrhs,localrhssize,MPI_DOUBLE,0,mesh->getXcomm()); // Scatters solution from host back to localrhs (which points to x's data) on all processors +// MPI_Scatterv(rhs,localrhs_size_array,rhs_positions,MPI_DOUBLE,localrhs,localrhssize,MPI_DOUBLE,0,localmesh->getXcomm()); // Scatters solution from host back to localrhs (which points to x's data) on all processors // // #if CHECK > 0 // msg_stack.pop(); @@ -516,17 +516,17 @@ LaplaceMumps::LaplaceMumps(Options *opt, const CELL_LOC loc) : // #if CHECK > 0 // msg_stack.push("Laplacian::solve(Field3D)"); // #endif -// int ys = mesh->ystart, ye = mesh->yend; +// int ys = localmesh->ystart, ye = localmesh->yend; // -// if(mesh->hasBndryLowerY()) { +// if(localmesh->hasBndryLowerY()) { // if (include_yguards) // ys = 0; // Mesh contains a lower boundary and we are solving in the guard cells // // ys += extra_yguards_lower; // } -// if(mesh->hasBndryUpperY()) { +// if(localmesh->hasBndryUpperY()) { // if (include_yguards) -// ye = mesh->LocalNy-1; // Contains upper boundary and we are solving in the guard cells +// ye = localmesh->LocalNy-1; // Contains upper boundary and we are solving in the guard cells // // ye -= extra_yguards_upper; // } @@ -563,36 +563,38 @@ const FieldPerp LaplaceMumps::solve(const FieldPerp &b, const FieldPerp &x0) { } const FieldPerp LaplaceMumps::solve(const FieldPerp &b) { + ASSERT1(localmesh == b.getMesh()); + int y = b.getIndex(); sol = 0.; sol.setIndex(y); // Set boundary conditions through rhs if needed if (!(inner_boundary_flags & INVERT_RHS)) { - if (mesh->firstX()) - for (int z=0; zLocalNz; z++) - for (int x=mesh->xstart-1; x>=0; x--) { + if (localmesh->firstX()) + for (int z=0; zLocalNz; z++) + for (int x=localmesh->xstart-1; x>=0; x--) { b[x][z]=0.; } } if (!(outer_boundary_flags & INVERT_RHS)) { - if (mesh->lastX()) - for (int z=0; zLocalNz; z++) - for (int x=mesh->xend+1; xLocalNx; x++) { + if (localmesh->lastX()) + for (int z=0; zLocalNz; z++) + for (int x=localmesh->xend+1; xLocalNx; x++) { b[x][z]=0.; } } BoutReal* bdata = *b.getData(); int xs,xe; - if (mesh->firstX()) xs=0; - else xs=mesh->xstart; - if (mesh->lastX()) xe=mesh->LocalNx-1; - else xe=mesh->xend; + if (localmesh->firstX()) xs=0; + else xs=localmesh->xstart; + if (localmesh->lastX()) xe=localmesh->LocalNx-1; + else xe=localmesh->xend; for (int x=xs; x<=xe; x++) - for (int z=0; zLocalNz; z++) - localrhs[(x-xs)*(mesh->LocalNz)+z] = bdata[x*mesh->LocalNz+z]; + for (int z=0; zLocalNz; z++) + localrhs[(x-xs)*(localmesh->LocalNz)+z] = bdata[x*localmesh->LocalNz+z]; MPI_Gatherv(localrhs,localrhssize,MPI_DOUBLE,rhs,localrhs_size_array,rhs_positions,MPI_DOUBLE,0,comm); @@ -602,8 +604,8 @@ const FieldPerp LaplaceMumps::solve(const FieldPerp &b) { BoutReal* soldata = *sol.getData(); for (int x=xs; x<=xe; x++) - for (int z=0; zLocalNz; z++) - soldata[x*mesh->LocalNz+z] = localrhs[(x-xs)*(mesh->LocalNz)+z]; + for (int z=0; zLocalNz; z++) + soldata[x*localmesh->LocalNz+z] = localrhs[(x-xs)*(localmesh->LocalNz)+z]; return sol; } @@ -613,14 +615,14 @@ void LaplaceMumps::solve(BoutReal* rhs, int y) { { Timer timer("mumpssetup"); int i = 0; - Coordinates *coord = mesh->coordinates(location); + Coordinates *coord = localmesh->coordinates(location); // Set Matrix Elements corresponding to index lists created in constructor (x,z) loop over rows - // X=0 to mesh->xstart-1 defines the boundary region of the domain. - if( mesh->firstX() ) - for(int x=0; xxstart; x++) - for(int z=0; zLocalNz; z++) { + // X=0 to localmesh->xstart-1 defines the boundary region of the domain. + if( localmesh->firstX() ) + for(int x=0; xxstart; x++) + for(int z=0; zLocalNz; z++) { // Set values corresponding to nodes adjacent in x if Neumann Boundary Conditions are required. if(inner_boundary_flags & INVERT_AC_GRAD) if( fourth_order ) { @@ -653,8 +655,8 @@ void LaplaceMumps::solve(BoutReal* rhs, int y) { } // Main domain with Laplacian operator - for(int x=mesh->xstart; x <= mesh->xend; x++) - for(int z=0; zLocalNz; z++) { + for(int x=localmesh->xstart; x <= localmesh->xend; x++) + for(int z=0; zLocalNz; z++) { BoutReal A0, A1, A2, A3, A4, A5; A0 = A[x][y][z]; Coeffs( x, y, z, A1, A2, A3, A4, A5 ); @@ -806,10 +808,10 @@ void LaplaceMumps::solve(BoutReal* rhs, int y) { } } - // X=mesh->xend+1 to mesh->LocalNx-1 defines the upper boundary region of the domain. - if( mesh->lastX() ) - for(int x=mesh->xend+1; xLocalNx; x++) - for(int z=0; zLocalNz; z++) { + // X=localmesh->xend+1 to localmesh->LocalNx-1 defines the upper boundary region of the domain. + if( localmesh->lastX() ) + for(int x=localmesh->xend+1; xLocalNx; x++) + for(int z=0; zLocalNz; z++) { // Set values corresponding to nodes adjacent in x if Neumann Boundary Conditions are required. if(outer_boundary_flags & INVERT_AC_GRAD) { @@ -856,7 +858,7 @@ void LaplaceMumps::solve(BoutReal* rhs, int y) { void LaplaceMumps::Coeffs( int x, int y, int z, BoutReal &coef1, BoutReal &coef2, BoutReal &coef3, BoutReal &coef4, BoutReal &coef5 ) { - Coordinates *coord = mesh->coordinates(location); + Coordinates *coord = localmesh->coordinates(location); coef1 = coord->g11[x][y]; // X 2nd derivative coefficient coef2 = coord->g33[x][y]; // Z 2nd derivative coefficient @@ -872,7 +874,7 @@ void LaplaceMumps::Coeffs( int x, int y, int z, BoutReal &coef1, BoutReal &coef2 if(nonuniform) { // non-uniform mesh correction - if((x != 0) && (x != (mesh->LocalNx-1))) + if((x != 0) && (x != (localmesh->LocalNx-1))) { //coef4 += coord->g11[jx][jy]*0.25*( (1.0/dx[jx+1][jy]) - (1.0/dx[jx-1][jy]) )/dx[jx][jy]; // SHOULD BE THIS (?) //coef4 -= 0.5 * ( ( coord->dx[x+1][y] - coord->dx[x-1][y] ) / SQ ( coord->dx[x][y] ) ) * coef1; // BOUT-06 term @@ -893,8 +895,8 @@ void LaplaceMumps::Coeffs( int x, int y, int z, BoutReal &coef1, BoutReal &coef2 // A second/fourth order derivative term if (issetC) { -// if( (x > 0) && (x < (mesh->LocalNx-1)) ) //Valid if doing second order derivative, not if fourth: should only be called for xstart<=x<=xend anyway - if( (x > 1) && (x < (mesh->LocalNx-2)) ) { +// if( (x > 0) && (x < (localmesh->LocalNx-1)) ) //Valid if doing second order derivative, not if fourth: should only be called for xstart<=x<=xend anyway + if( (x > 1) && (x < (localmesh->LocalNx-2)) ) { int zp = z+1; if (zp > meshz-1) zp -= meshz; int zm = z-1; diff --git a/src/invert/laplace/impls/mumps/mumps_laplace.hxx b/src/invert/laplace/impls/mumps/mumps_laplace.hxx index bac607afa0..6c37e63fa9 100644 --- a/src/invert/laplace/impls/mumps/mumps_laplace.hxx +++ b/src/invert/laplace/impls/mumps/mumps_laplace.hxx @@ -36,7 +36,7 @@ class LaplaceMumps; class LaplaceMumps : public Laplacian { public: - LaplaceMumps(Options *UNUSED(opt) = nullptr, const CELL_LOC UNUSED(loc) = CELL_CENTRE) { + LaplaceMumps(Options *UNUSED(opt) = nullptr, const CELL_LOC UNUSED(loc) = CELL_CENTRE, Mesh *UNUSED(mesh_in) = mesh) { throw BoutException("Mumps library not available"); } @@ -76,7 +76,7 @@ public: class LaplaceMumps : public Laplacian { public: - LaplaceMumps(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE); + LaplaceMumps(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = mesh); ~LaplaceMumps() { mumps_struc.job = -2; dmumps_c(&mumps_struc); @@ -88,72 +88,86 @@ public: 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; issetC = true; } void setCoefC1(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); C1 = val; issetC = true; } void setCoefC2(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); C2 = val; issetC = true; } void setCoefD(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); D = val; issetD = true; } void setCoefEx(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); Ex = val; issetE = true; } void setCoefEz(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); Ez = val; issetE = true; } 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; issetC = true; } void setCoefC1(const Field3D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); C1 = val; issetC = true; } void setCoefC2(const Field3D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); C2 = val; issetC = true; } void setCoefD(const Field3D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); D = val; issetD = true; } void setCoefEx(const Field3D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); Ex = val; issetE = true; } void setCoefEz(const Field3D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); Ez = val; issetE = true; } diff --git a/src/invert/laplace/impls/naulin/naulin_laplace.cxx b/src/invert/laplace/impls/naulin/naulin_laplace.cxx index 81e1297899..87e1093a61 100644 --- a/src/invert/laplace/impls/naulin/naulin_laplace.cxx +++ b/src/invert/laplace/impls/naulin/naulin_laplace.cxx @@ -120,8 +120,8 @@ #include "naulin_laplace.hxx" -LaplaceNaulin::LaplaceNaulin(Options *opt, const CELL_LOC loc) - : Laplacian(opt, loc), Acoef(0.0), C1coef(1.0), C2coef(0.0), Dcoef(1.0), +LaplaceNaulin::LaplaceNaulin(Options *opt, const CELL_LOC loc, Mesh *mesh_in) + : Laplacian(opt, loc, mesh_in), Acoef(0.0), C1coef(1.0), C2coef(0.0), Dcoef(1.0), delp2solver(nullptr), naulinsolver_mean_its(0.), ncalls(0) { ASSERT1(opt != nullptr); // An Options pointer should always be passed in by LaplaceFactory @@ -135,7 +135,7 @@ LaplaceNaulin::LaplaceNaulin(Options *opt, const CELL_LOC loc) OPTION(opt, rtol, 1.e-7); OPTION(opt, atol, 1.e-20); OPTION(opt, maxits, 100); - delp2solver = create(opt->getSection("delp2solver"), location); + delp2solver = create(opt->getSection("delp2solver"), location, localmesh); std::string delp2type; opt->getSection("delp2solver")->get("type", delp2type, "cyclic"); // Check delp2solver is using an FFT scheme, otherwise it will not exactly @@ -170,8 +170,8 @@ const Field3D LaplaceNaulin::solve(const Field3D &rhs, const Field3D &x0) { ASSERT1(C1coef.getLocation() == location); ASSERT1(C2coef.getLocation() == location); ASSERT1(Acoef.getLocation() == location); + ASSERT1(localmesh == rhs.getMesh() && localmesh == x0.getMesh()); - Mesh *mesh = rhs.getMesh(); Coordinates *coords = rhs.getCoordinates(); Field3D x(x0); // Result @@ -203,12 +203,12 @@ const Field3D LaplaceNaulin::solve(const Field3D &rhs, const Field3D &x0) { if ( (inner_boundary_flags & INVERT_SET) || (outer_boundary_flags & INVERT_SET) ) // This passes in the boundary conditions from x0's guard cells - copy_x_boundaries(x, x0, mesh); + copy_x_boundaries(x, x0, localmesh); // NB need to pass x in case boundary flags require 'x0', even if // delp2solver is not iterative and does not use an initial guess x = delp2solver->solve(b, x); - mesh->communicate(x); + localmesh->communicate(x); // re-calculate the rhs from the new solution // Use here to calculate an error, can also use for the next iteration @@ -233,17 +233,17 @@ const Field3D LaplaceNaulin::solve(const Field3D &rhs, const Field3D &x0) { return x; } -void LaplaceNaulin::copy_x_boundaries(Field3D &x, const Field3D &x0, Mesh *mesh) { - if (mesh->firstX()) { - for (int i=mesh->xstart-1; i>=0; i--) - for (int j=mesh->ystart; j<=mesh->yend; j++) - for (int k=0; kLocalNz; k++) +void LaplaceNaulin::copy_x_boundaries(Field3D &x, const Field3D &x0, Mesh *localmesh) { + if (localmesh->firstX()) { + for (int i=localmesh->xstart-1; i>=0; i--) + for (int j=localmesh->ystart; j<=localmesh->yend; j++) + for (int k=0; kLocalNz; k++) x(i, j, k) = x0(i, j, k); } - if (mesh->lastX()) { - for (int i=mesh->xend+1; iLocalNx; i++) - for (int j=mesh->ystart; j<=mesh->yend; j++) - for (int k=0; kLocalNz; k++) + if (localmesh->lastX()) { + for (int i=localmesh->xend+1; iLocalNx; i++) + for (int j=localmesh->ystart; j<=localmesh->yend; j++) + for (int k=0; kLocalNz; k++) x(i, j, k) = x0(i, j, k); } } diff --git a/src/invert/laplace/impls/naulin/naulin_laplace.hxx b/src/invert/laplace/impls/naulin/naulin_laplace.hxx index c8d2e8c359..b58d63f678 100644 --- a/src/invert/laplace/impls/naulin/naulin_laplace.hxx +++ b/src/invert/laplace/impls/naulin/naulin_laplace.hxx @@ -37,7 +37,7 @@ class LaplaceNaulin; */ class LaplaceNaulin : public Laplacian { public: - LaplaceNaulin(Options *opt = NULL, const CELL_LOC loc = CELL_CENTRE); + LaplaceNaulin(Options *opt = NULL, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = mesh); ~LaplaceNaulin(); // ACoef is not implemented because the delp2solver that we use can probably @@ -45,44 +45,54 @@ public: // where we allow Dcoef to be a Field3D void setCoefA(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); Acoef = val; } void setCoefA(const Field3D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); Acoef = val; } void setCoefC(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); setCoefC1(val); setCoefC2(val); } void setCoefC(const Field3D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); setCoefC1(val); setCoefC2(val); } void setCoefC1(const Field3D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); C1coef = val; } void setCoefC1(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); C1coef = val; } void setCoefC2(const Field3D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); C2coef = val; } void setCoefC2(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); C2coef = val; } void setCoefD(const Field3D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); Dcoef = val; } void setCoefD(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); Dcoef = val; } void setCoefEx(const Field2D &UNUSED(val)) override { diff --git a/src/invert/laplace/impls/pdd/pdd.cxx b/src/invert/laplace/impls/pdd/pdd.cxx index f71b551ccb..4a2fbbf46a 100644 --- a/src/invert/laplace/impls/pdd/pdd.cxx +++ b/src/invert/laplace/impls/pdd/pdd.cxx @@ -39,9 +39,11 @@ #include "pdd.hxx" const FieldPerp LaplacePDD::solve(const FieldPerp &b) { + ASSERT1(localmesh == b.getMesh()); + PDD_data data; - FieldPerp x(b.getMesh()); + FieldPerp x(localmesh); x.allocate(); start(b, data); @@ -52,17 +54,18 @@ const FieldPerp LaplacePDD::solve(const FieldPerp &b) { } const Field3D LaplacePDD::solve(const Field3D &b) { - Mesh *mesh = b.getMesh(); - Field3D x(mesh); + ASSERT1(localmesh == b.getMesh()); + + Field3D x(localmesh); x.allocate(); - FieldPerp xperp(mesh); + FieldPerp xperp(localmesh); xperp.allocate(); - int ys = mesh->ystart, ye = mesh->yend; - if(mesh->hasBndryLowerY()) + int ys = localmesh->ystart, ye = localmesh->yend; + if(localmesh->hasBndryLowerY()) ys = 0; // Mesh contains a lower boundary - if(mesh->hasBndryUpperY()) - ye = mesh->LocalNy-1; // Contains upper boundary + if(localmesh->hasBndryUpperY()) + ye = localmesh->LocalNy-1; // Contains upper boundary if(low_mem) { // Solve one slice at a time @@ -110,37 +113,38 @@ const Field3D LaplacePDD::solve(const Field3D &b) { /// @param[in] b RHS values (Ax = b) /// @param[in] data Internal data used for multiple calls in parallel mode void LaplacePDD::start(const FieldPerp &b, PDD_data &data) { + ASSERT1(localmesh == b.getMesh()); + int ix, kz; - Mesh *mesh = b.getMesh(); - int ncz = mesh->LocalNz; + int ncz = localmesh->LocalNz; data.jy = b.getIndex(); - if(mesh->firstX() && mesh->lastX()) + if(localmesh->firstX() && localmesh->lastX()) throw BoutException("Error: PDD method only works for NXPE > 1\n"); - if(mesh->periodicX) { - throw BoutException("LaplacePDD does not work with periodicity in the x direction (mesh->PeriodicX == true). Change boundary conditions or use serial-tri or cyclic solver instead"); + if(localmesh->periodicX) { + throw BoutException("LaplacePDD does not work with periodicity in the x direction (localmesh->PeriodicX == true). Change boundary conditions or use serial-tri or cyclic solver instead"); } if (data.bk.empty()) { // Need to allocate working memory // RHS vector - data.bk = Matrix(maxmode + 1, mesh->LocalNx); + data.bk = Matrix(maxmode + 1, localmesh->LocalNx); // Matrix to be solved - data.avec = Matrix(maxmode + 1, mesh->LocalNx); - data.bvec = Matrix(maxmode + 1, mesh->LocalNx); - data.cvec = Matrix(maxmode + 1, mesh->LocalNx); + data.avec = Matrix(maxmode + 1, localmesh->LocalNx); + data.bvec = Matrix(maxmode + 1, localmesh->LocalNx); + data.cvec = Matrix(maxmode + 1, localmesh->LocalNx); // Working vectors - data.v = Matrix(maxmode + 1, mesh->LocalNx); - data.w = Matrix(maxmode + 1, mesh->LocalNx); + data.v = Matrix(maxmode + 1, localmesh->LocalNx); + data.w = Matrix(maxmode + 1, localmesh->LocalNx); // Result - data.xk = Matrix(maxmode + 1, mesh->LocalNx); + data.xk = Matrix(maxmode + 1, localmesh->LocalNx); // Communication buffers. Space for 2 complex values for each kz data.snd = Array(4 * (maxmode + 1)); @@ -152,7 +156,7 @@ void LaplacePDD::start(const FieldPerp &b, PDD_data &data) { /// Take FFTs of data Array bk1d(ncz / 2 + 1); ///< 1D in Z for taking FFTs - for(ix=0; ix < mesh->LocalNx; ix++) { + for(ix=0; ix < localmesh->LocalNx; ix++) { rfft(b[ix], ncz, std::begin(bk1d)); for(kz = 0; kz <= maxmode; kz++) data.bk(kz, ix) = bk1d[kz]; @@ -160,7 +164,7 @@ void LaplacePDD::start(const FieldPerp &b, PDD_data &data) { /// Create the matrices to be inverted (one for each z point) - BoutReal kwaveFactor = 2.0 * PI / mesh->getCoordinates(location)->zlength(); + BoutReal kwaveFactor = 2.0 * PI / localmesh->getCoordinates(location)->zlength(); /// Set matrix elements for (int kz = 0; kz <= maxmode; kz++) { @@ -169,8 +173,8 @@ void LaplacePDD::start(const FieldPerp &b, PDD_data &data) { outer_boundary_flags, &Acoef, &Ccoef, &Dcoef); } - Array e(mesh->LocalNx); - for (ix = 0; ix < mesh->LocalNx; ix++) + Array e(localmesh->LocalNx); + for (ix = 0; ix < localmesh->LocalNx; ix++) e[ix] = 0.0; // Do we need this? for(kz = 0; kz <= maxmode; kz++) { @@ -180,51 +184,51 @@ void LaplacePDD::start(const FieldPerp &b, PDD_data &data) { dcomplex v0, x0; // Values to be sent to processor i-1 - if(mesh->firstX()) { + if(localmesh->firstX()) { // Domain includes inner boundary tridag(&data.avec(kz, 0), &data.bvec(kz, 0), &data.cvec(kz, 0), &data.bk(kz, 0), - &data.xk(kz, 0), mesh->xend + 1); + &data.xk(kz, 0), localmesh->xend + 1); // Add C (row m-1) from next processor - e[mesh->xend] = data.cvec(kz, mesh->xend); + e[localmesh->xend] = data.cvec(kz, localmesh->xend); tridag(&data.avec(kz, 0), &data.bvec(kz, 0), &data.cvec(kz, 0), std::begin(e), - &data.w(kz, 0), mesh->xend + 1); + &data.w(kz, 0), localmesh->xend + 1); - }else if(mesh->lastX()) { + }else if(localmesh->lastX()) { // Domain includes outer boundary - tridag(&data.avec(kz, mesh->xstart), &data.bvec(kz, mesh->xstart), - &data.cvec(kz, mesh->xstart), &data.bk(kz, mesh->xstart), - &data.xk(kz, mesh->xstart), mesh->xend - mesh->xend + 1); + tridag(&data.avec(kz, localmesh->xstart), &data.bvec(kz, localmesh->xstart), + &data.cvec(kz, localmesh->xstart), &data.bk(kz, localmesh->xstart), + &data.xk(kz, localmesh->xstart), localmesh->xend - localmesh->xend + 1); // Add A (row 0) from previous processor - e[0] = data.avec(kz, mesh->xstart); - tridag(&data.avec(kz, mesh->xstart), &data.bvec(kz, mesh->xstart), - &data.cvec(kz, mesh->xstart), std::begin(e), &data.v(kz, mesh->xstart), - mesh->xend + 1); + e[0] = data.avec(kz, localmesh->xstart); + tridag(&data.avec(kz, localmesh->xstart), &data.bvec(kz, localmesh->xstart), + &data.cvec(kz, localmesh->xstart), std::begin(e), &data.v(kz, localmesh->xstart), + localmesh->xend + 1); - x0 = data.xk(kz, mesh->xstart); - v0 = data.v(kz, mesh->xstart); + x0 = data.xk(kz, localmesh->xstart); + v0 = data.v(kz, localmesh->xstart); }else { // No boundaries - tridag(&data.avec(kz, mesh->xstart), &data.bvec(kz, mesh->xstart), - &data.cvec(kz, mesh->xstart), &data.bk(kz, mesh->xstart), - &data.xk(kz, mesh->xstart), mesh->xend - mesh->xstart + 1); + tridag(&data.avec(kz, localmesh->xstart), &data.bvec(kz, localmesh->xstart), + &data.cvec(kz, localmesh->xstart), &data.bk(kz, localmesh->xstart), + &data.xk(kz, localmesh->xstart), localmesh->xend - localmesh->xstart + 1); // Add A (row 0) from previous processor - e[0] = data.avec(kz, mesh->xstart); - tridag(&data.avec(kz, mesh->xstart), &data.bvec(kz, mesh->xstart), - &data.cvec(kz, mesh->xstart), &e[mesh->xstart], &data.v(kz, mesh->xstart), - mesh->xend - mesh->xstart + 1); + e[0] = data.avec(kz, localmesh->xstart); + tridag(&data.avec(kz, localmesh->xstart), &data.bvec(kz, localmesh->xstart), + &data.cvec(kz, localmesh->xstart), &e[localmesh->xstart], &data.v(kz, localmesh->xstart), + localmesh->xend - localmesh->xstart + 1); e[0] = 0.0; // Add C (row m-1) from next processor - e[mesh->xend] = data.cvec(kz, mesh->xend); - tridag(&data.avec(kz, mesh->xstart), &data.bvec(kz, mesh->xstart), - &data.cvec(kz, mesh->xstart), &e[mesh->xstart], &data.v(kz, mesh->xstart), - mesh->xend - mesh->xstart + 1); - e[mesh->xend] = 0.0; + e[localmesh->xend] = data.cvec(kz, localmesh->xend); + tridag(&data.avec(kz, localmesh->xstart), &data.bvec(kz, localmesh->xstart), + &data.cvec(kz, localmesh->xstart), &e[localmesh->xstart], &data.v(kz, localmesh->xstart), + localmesh->xend - localmesh->xstart + 1); + e[localmesh->xend] = 0.0; } // Put values into communication buffers @@ -236,17 +240,17 @@ void LaplacePDD::start(const FieldPerp &b, PDD_data &data) { // Stage 3: Communicate x0, v0 from node i to i-1 - if(!mesh->lastX()) { + if(!localmesh->lastX()) { // All except the last processor expect to receive data // Post async receive data.recv_handle = - mesh->irecvXOut(std::begin(data.rcv), 4 * (maxmode + 1), PDD_COMM_XV); + localmesh->irecvXOut(std::begin(data.rcv), 4 * (maxmode + 1), PDD_COMM_XV); } - if(!mesh->firstX()) { + if(!localmesh->firstX()) { // Send the data - mesh->sendXIn(std::begin(data.snd), 4 * (maxmode + 1), PDD_COMM_XV); + localmesh->sendXIn(std::begin(data.snd), 4 * (maxmode + 1), PDD_COMM_XV); } } @@ -255,8 +259,8 @@ void LaplacePDD::start(const FieldPerp &b, PDD_data &data) { void LaplacePDD::next(PDD_data &data) { // Wait for x0 and v0 to arrive from processor i+1 - if(!mesh->lastX()) { - mesh->wait(data.recv_handle); + if(!localmesh->lastX()) { + localmesh->wait(data.recv_handle); /*! Now solving on all except the last processor * @@ -273,18 +277,18 @@ void LaplacePDD::next(PDD_data &data) { x0 = dcomplex(data.rcv[4*kz], data.rcv[4*kz+1]); v0 = dcomplex(data.rcv[4*kz+2], data.rcv[4*kz+3]); - data.y2i[kz] = (data.xk(kz, mesh->xend) - data.w(kz, mesh->xend) * x0) / - (1. - data.w(kz, mesh->xend) * v0); + data.y2i[kz] = (data.xk(kz, localmesh->xend) - data.w(kz, localmesh->xend) * x0) / + (1. - data.w(kz, localmesh->xend) * v0); } } - if(!mesh->firstX()) { + if(!localmesh->firstX()) { // All except pe=0 receive values from i-1. Posting async receive data.recv_handle = - mesh->irecvXIn(std::begin(data.rcv), 2 * (maxmode + 1), PDD_COMM_Y); + localmesh->irecvXIn(std::begin(data.rcv), 2 * (maxmode + 1), PDD_COMM_Y); } - if(!mesh->lastX()) { + if(!localmesh->lastX()) { // Send value to the (i+1)th processor for(int kz = 0; kz <= maxmode; kz++) { @@ -292,7 +296,7 @@ void LaplacePDD::next(PDD_data &data) { data.snd[2*kz+1] = data.y2i[kz].imag(); } - mesh->sendXOut(std::begin(data.snd), 2 * (maxmode + 1), PDD_COMM_Y); + localmesh->sendXOut(std::begin(data.snd), 2 * (maxmode + 1), PDD_COMM_Y); } } @@ -303,32 +307,32 @@ void LaplacePDD::finish(PDD_data &data, FieldPerp &x) { x.allocate(); x.setIndex(data.jy); - if(!mesh->lastX()) { + if(!localmesh->lastX()) { for(kz = 0; kz <= maxmode; kz++) { - for(ix=0; ix < mesh->LocalNx; ix++) + for(ix=0; ix < localmesh->LocalNx; ix++) data.xk(kz, ix) -= data.w(kz, ix) * data.y2i[kz]; } } - if(!mesh->firstX()) { - mesh->wait(data.recv_handle); + if(!localmesh->firstX()) { + localmesh->wait(data.recv_handle); for(kz = 0; kz <= maxmode; kz++) { dcomplex y2m = dcomplex(data.rcv[2*kz], data.rcv[2*kz+1]); - for(ix=0; ix < mesh->LocalNx; ix++) + for(ix=0; ix < localmesh->LocalNx; ix++) data.xk(kz, ix) -= data.v(kz, ix) * y2m; } } // Have result in Fourier space. Convert back to BoutReal space - int ncz = mesh->LocalNz; + int ncz = localmesh->LocalNz; Array xk1d(ncz / 2 + 1); ///< 1D in Z for taking FFTs for (kz = maxmode; kz <= ncz / 2; kz++) xk1d[kz] = 0.0; - for(ix=0; ixLocalNx; ix++){ + for(ix=0; ixLocalNx; ix++){ for(kz = 0; kz <= maxmode; kz++) { xk1d[kz] = data.xk(kz, ix); diff --git a/src/invert/laplace/impls/pdd/pdd.hxx b/src/invert/laplace/impls/pdd/pdd.hxx index e9fdccf342..0dff1dd6e0 100644 --- a/src/invert/laplace/impls/pdd/pdd.hxx +++ b/src/invert/laplace/impls/pdd/pdd.hxx @@ -40,8 +40,8 @@ class LaplacePDD; class LaplacePDD : public Laplacian { public: - LaplacePDD(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE) - : Laplacian(opt, loc), Acoef(0.0), Ccoef(1.0), Dcoef(1.0), PDD_COMM_XV(123), + LaplacePDD(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = mesh) + : Laplacian(opt, loc, mesh_in), Acoef(0.0), Ccoef(1.0), Dcoef(1.0), PDD_COMM_XV(123), PDD_COMM_Y(456) { Acoef.setLocation(location); Ccoef.setLocation(location); @@ -52,16 +52,19 @@ public: using Laplacian::setCoefA; void setCoefA(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); Acoef = val; } using Laplacian::setCoefC; void setCoefC(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); Ccoef = val; } using Laplacian::setCoefD; void setCoefD(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); Dcoef = val; } using Laplacian::setCoefEx; diff --git a/src/invert/laplace/impls/petsc/petsc_laplace.cxx b/src/invert/laplace/impls/petsc/petsc_laplace.cxx index 27f01a3161..1555345456 100644 --- a/src/invert/laplace/impls/petsc/petsc_laplace.cxx +++ b/src/invert/laplace/impls/petsc/petsc_laplace.cxx @@ -57,8 +57,8 @@ static PetscErrorCode laplacePCapply(PC pc,Vec x,Vec y) { PetscFunctionReturn(s->precon(x, y)); } -LaplacePetsc::LaplacePetsc(Options *opt, const CELL_LOC loc) : - Laplacian(opt, loc), +LaplacePetsc::LaplacePetsc(Options *opt, const CELL_LOC loc, Mesh *mesh_in) : + Laplacian(opt, loc, mesh_in), A(0.0), C1(1.0), C2(1.0), D(1.0), Ex(0.0), Ez(0.0), issetD(false), issetC(false), issetE(false) { @@ -92,28 +92,28 @@ LaplacePetsc::LaplacePetsc(Options *opt, const CELL_LOC loc) : if ( outer_boundary_flags & ~implemented_boundary_flags ) { throw BoutException("Attempted to set Laplacian inversion boundary flag that is not implemented in petsc_laplace.cxx"); } - if(mesh->periodicX) { - throw BoutException("LaplacePetsc does not work with periodicity in the x direction (mesh->PeriodicX == true). Change boundary conditions or use serial-tri or cyclic solver instead"); + if(localmesh->periodicX) { + throw BoutException("LaplacePetsc does not work with periodicity in the x direction (localmesh->PeriodicX == true). Change boundary conditions or use serial-tri or cyclic solver instead"); } #endif // Get communicator for group of processors in X - all points in z-x plane for fixed y. - comm = mesh->getXcomm(); + comm = localmesh->getXcomm(); // Need to determine local size to use based on prior parallelisation // Coefficient values are stored only on local processors. - localN = (mesh->xend - mesh->xstart + 1) * (mesh->LocalNz); - if(mesh->firstX()) - localN += mesh->xstart * (mesh->LocalNz); // If on first processor add on width of boundary region - if(mesh->lastX()) - localN += mesh->xstart * (mesh->LocalNz); // If on last processor add on width of boundary region + localN = (localmesh->xend - localmesh->xstart + 1) * (localmesh->LocalNz); + if(localmesh->firstX()) + localN += localmesh->xstart * (localmesh->LocalNz); // If on first processor add on width of boundary region + if(localmesh->lastX()) + localN += localmesh->xstart * (localmesh->LocalNz); // If on last processor add on width of boundary region // Calculate 'size' (the total number of points in physical grid) if(MPI_Allreduce(&localN, &size, 1, MPI_INT, MPI_SUM, comm) != MPI_SUCCESS) throw BoutException("Error in MPI_Allreduce during LaplacePetsc initialisation"); // Calculate total (physical) grid dimensions - meshz = mesh->LocalNz; + meshz = localmesh->LocalNz; meshx = size / meshz; // Create PETSc type of vectors for the solution and the RHS vector @@ -143,43 +143,43 @@ LaplacePetsc::LaplacePetsc(Options *opt, const CELL_LOC loc) : PetscMalloc( (localN)*sizeof(PetscInt), &d_nnz ); PetscMalloc( (localN)*sizeof(PetscInt), &o_nnz ); if (fourth_order) { - // first and last 2*mesh-LocalNz entries are the edge x-values that (may) have 'off-diagonal' components (i.e. on another processor) - if ( mesh->firstX() && mesh->lastX() ) { - for (int i=0; iLocalNz; i++) { + // first and last 2*localmesh-LocalNz entries are the edge x-values that (may) have 'off-diagonal' components (i.e. on another processor) + if ( localmesh->firstX() && localmesh->lastX() ) { + for (int i=0; iLocalNz; i++) { d_nnz[i]=15; d_nnz[localN-1-i]=15; o_nnz[i]=0; o_nnz[localN-1-i]=0; } - for (int i=(mesh->LocalNz); i<2*(mesh->LocalNz); i++) { + for (int i=(localmesh->LocalNz); i<2*(localmesh->LocalNz); i++) { d_nnz[i]=20; d_nnz[localN-1-i]=20; o_nnz[i]=0; o_nnz[localN-1-i]=0; } } - else if ( mesh->firstX() ) { - for (int i=0; iLocalNz; i++) { + else if ( localmesh->firstX() ) { + for (int i=0; iLocalNz; i++) { d_nnz[i]=15; d_nnz[localN-1-i]=15; o_nnz[i]=0; o_nnz[localN-1-i]=10; } - for (int i=(mesh->LocalNz); i<2*(mesh->LocalNz); i++) { + for (int i=(localmesh->LocalNz); i<2*(localmesh->LocalNz); i++) { d_nnz[i]=20; d_nnz[localN-1-i]=20; o_nnz[i]=0; o_nnz[localN-1-i]=5; } } - else if ( mesh->lastX() ) { - for (int i=0; iLocalNz; i++) { + else if ( localmesh->lastX() ) { + for (int i=0; iLocalNz; i++) { d_nnz[i]=15; d_nnz[localN-1-i]=15; o_nnz[i]=10; o_nnz[localN-1-i]=0; } - for (int i=(mesh->LocalNz); i<2*(mesh->LocalNz); i++) { + for (int i=(localmesh->LocalNz); i<2*(localmesh->LocalNz); i++) { d_nnz[i]=20; d_nnz[localN-1-i]=20; o_nnz[i]=5; @@ -187,13 +187,13 @@ LaplacePetsc::LaplacePetsc(Options *opt, const CELL_LOC loc) : } } else { - for (int i=0; iLocalNz; i++) { + for (int i=0; iLocalNz; i++) { d_nnz[i]=15; d_nnz[localN-1-i]=15; o_nnz[i]=10; o_nnz[localN-1-i]=10; } - for (int i=(mesh->LocalNz); i<2*(mesh->LocalNz); i++) { + for (int i=(localmesh->LocalNz); i<2*(localmesh->LocalNz); i++) { d_nnz[i]=20; d_nnz[localN-1-i]=20; o_nnz[i]=5; @@ -201,7 +201,7 @@ LaplacePetsc::LaplacePetsc(Options *opt, const CELL_LOC loc) : } } - for (int i=2*(mesh->LocalNz); iLocalNz));i++) { + for (int i=2*(localmesh->LocalNz); iLocalNz));i++) { d_nnz[i]=25; d_nnz[localN-1-i]=25; o_nnz[i]=0; @@ -209,7 +209,7 @@ LaplacePetsc::LaplacePetsc(Options *opt, const CELL_LOC loc) : } // Use d_nnz and o_nnz for preallocating the matrix - if (mesh->firstX() && mesh->lastX()) { + if (localmesh->firstX() && localmesh->lastX()) { // Only one processor in X MatSeqAIJSetPreallocation( MatA, 0, d_nnz ); }else { @@ -217,25 +217,25 @@ LaplacePetsc::LaplacePetsc(Options *opt, const CELL_LOC loc) : } } else { - // first and last mesh-LocalNz entries are the edge x-values that (may) have 'off-diagonal' components (i.e. on another processor) - if ( mesh->firstX() && mesh->lastX() ) { - for (int i=0; iLocalNz; i++) { + // first and last localmesh->LocalNz entries are the edge x-values that (may) have 'off-diagonal' components (i.e. on another processor) + if ( localmesh->firstX() && localmesh->lastX() ) { + for (int i=0; iLocalNz; i++) { d_nnz[i]=6; d_nnz[localN-1-i]=6; o_nnz[i]=0; o_nnz[localN-1-i]=0; } } - else if ( mesh->firstX() ) { - for (int i=0; iLocalNz; i++) { + else if ( localmesh->firstX() ) { + for (int i=0; iLocalNz; i++) { d_nnz[i]=6; d_nnz[localN-1-i]=6; o_nnz[i]=0; o_nnz[localN-1-i]=3; } } - else if ( mesh->lastX() ) { - for (int i=0; iLocalNz; i++) { + else if ( localmesh->lastX() ) { + for (int i=0; iLocalNz; i++) { d_nnz[i]=6; d_nnz[localN-1-i]=6; o_nnz[i]=3; @@ -243,7 +243,7 @@ LaplacePetsc::LaplacePetsc(Options *opt, const CELL_LOC loc) : } } else { - for (int i=0; iLocalNz; i++) { + for (int i=0; iLocalNz; i++) { d_nnz[i]=6; d_nnz[localN-1-i]=6; o_nnz[i]=3; @@ -251,7 +251,7 @@ LaplacePetsc::LaplacePetsc(Options *opt, const CELL_LOC loc) : } } - for (int i=mesh->LocalNz; iLocalNz);i++) { + for (int i=localmesh->LocalNz; iLocalNz);i++) { d_nnz[i]=9; d_nnz[localN-1-i]=9; o_nnz[i]=0; @@ -259,7 +259,7 @@ LaplacePetsc::LaplacePetsc(Options *opt, const CELL_LOC loc) : } // Use d_nnz and o_nnz for preallocating the matrix - if (mesh->firstX() && mesh->lastX()) { + if (localmesh->firstX() && localmesh->lastX()) { MatSeqAIJSetPreallocation( MatA, 0, d_nnz ); } else { MatMPIAIJSetPreallocation( MatA, 0, d_nnz, 0, o_nnz ); @@ -344,6 +344,8 @@ const FieldPerp LaplacePetsc::solve(const FieldPerp &b) { */ const FieldPerp LaplacePetsc::solve(const FieldPerp &b, const FieldPerp &x0) { TRACE("LaplacePetsc::solve"); + + ASSERT1(localmesh == b.getMesh() && localmesh == x0.getMesh()); #if CHECK > 0 // Checking flags are set to something which is not implemented (see @@ -361,7 +363,7 @@ const FieldPerp LaplacePetsc::solve(const FieldPerp &b, const FieldPerp &x0) { #endif // Get the metric tensor - Coordinates* coord = mesh->getCoordinates(location); + Coordinates* coord = localmesh->getCoordinates(location); int y = b.getIndex(); // Get the Y index sol.setIndex(y); // Initialize the solution field. @@ -389,11 +391,11 @@ const FieldPerp LaplacePetsc::solve(const FieldPerp &b, const FieldPerp &x0) { * In other word the indexing is done in a row-major order, but starting at * bottom left rather than top left */ - // X=0 to mesh->xstart-1 defines the boundary region of the domain. + // X=0 to localmesh->xstart-1 defines the boundary region of the domain. // Set the values for the inner boundary region - if( mesh->firstX() ) { - for(int x=0; xxstart; x++) { - for(int z=0; zLocalNz; z++) { + if( localmesh->firstX() ) { + for(int x=0; xxstart; x++) { + for(int z=0; zLocalNz; z++) { PetscScalar val; // Value of element to be set in the matrix // If Neumann Boundary Conditions are set. if(inner_boundary_flags & INVERT_AC_GRAD) { @@ -460,8 +462,8 @@ const FieldPerp LaplacePetsc::solve(const FieldPerp &b, const FieldPerp &x0) { } // Set the values for the main domain - for(int x=mesh->xstart; x <= mesh->xend; x++) { - for(int z=0; zLocalNz; z++) { + for(int x=localmesh->xstart; x <= localmesh->xend; x++) { + for(int z=0; zLocalNz; z++) { // NOTE: Only A0 is the A from setCoefA () BoutReal A0, A1, A2, A3, A4, A5; A0 = A(x,y,z); @@ -635,11 +637,11 @@ const FieldPerp LaplacePetsc::solve(const FieldPerp &b, const FieldPerp &x0) { } } - // X=mesh->xend+1 to mesh->LocalNx-1 defines the upper boundary region of the domain. + // X=localmesh->xend+1 to localmesh->LocalNx-1 defines the upper boundary region of the domain. // Set the values for the outer boundary region - if( mesh->lastX() ) { - for(int x=mesh->xend+1; xLocalNx; x++) { - for(int z=0; zLocalNz; z++) { + if( localmesh->lastX() ) { + for(int x=localmesh->xend+1; xLocalNx; x++) { + for(int z=0; zLocalNz; z++) { // Set Diagonal Values to 1 PetscScalar val = 1; Element(i,x,z, 0, 0, val, MatA ); @@ -803,9 +805,9 @@ const FieldPerp LaplacePetsc::solve(const FieldPerp &b, const FieldPerp &x0) { // Add data to FieldPerp Object i = Istart; // Set the inner boundary values - if(mesh->firstX()) { - for(int x=0; xxstart; x++) { - for(int z=0; zLocalNz; z++) { + if(localmesh->firstX()) { + for(int x=0; xxstart; x++) { + for(int z=0; zLocalNz; z++) { PetscScalar val = 0; VecGetValues(xs, 1, &i, &val ); sol[x][z] = val; @@ -815,8 +817,8 @@ const FieldPerp LaplacePetsc::solve(const FieldPerp &b, const FieldPerp &x0) { } // Set the main domain values - for(int x=mesh->xstart; x <= mesh->xend; x++) { - for(int z=0; zLocalNz; z++) { + for(int x=localmesh->xstart; x <= localmesh->xend; x++) { + for(int z=0; zLocalNz; z++) { PetscScalar val = 0; VecGetValues(xs, 1, &i, &val ); sol[x][z] = val; @@ -825,9 +827,9 @@ const FieldPerp LaplacePetsc::solve(const FieldPerp &b, const FieldPerp &x0) { } // Set the outer boundary values - if(mesh->lastX()) { - for(int x=mesh->xend+1; xLocalNx; x++) { - for(int z=0;z < mesh->LocalNz; z++) { + if(localmesh->lastX()) { + for(int x=localmesh->xend+1; xLocalNx; x++) { + for(int z=0;z < localmesh->LocalNz; z++) { PetscScalar val = 0; VecGetValues(xs, 1, &i, &val ); sol[x][z] = val; @@ -871,7 +873,7 @@ void LaplacePetsc::Element(int i, int x, int z, // Calculate the row to be set int row_new = x + xshift; // should never be out of range. - if( !mesh->firstX() ) row_new += (xoffset - mesh->xstart); + if( !localmesh->firstX() ) row_new += (xoffset - localmesh->xstart); // Calculate the column to be set int col_new = z + zshift; @@ -935,7 +937,7 @@ void LaplacePetsc::Element(int i, int x, int z, */ void LaplacePetsc::Coeffs( int x, int y, int z, BoutReal &coef1, BoutReal &coef2, BoutReal &coef3, BoutReal &coef4, BoutReal &coef5 ) { - Coordinates *coord = mesh->getCoordinates(location); // Get metric tensor + Coordinates *coord = localmesh->getCoordinates(location); // Get metric tensor coef1 = coord->g11(x,y); // X 2nd derivative coefficient coef2 = coord->g33(x,y); // Z 2nd derivative coefficient @@ -954,12 +956,12 @@ void LaplacePetsc::Coeffs( int x, int y, int z, BoutReal &coef1, BoutReal &coef2 if(nonuniform) { // non-uniform mesh correction - if((x != 0) && (x != (mesh->LocalNx-1))) { + if((x != 0) && (x != (localmesh->LocalNx-1))) { coef4 -= 0.5 * ( ( coord->dx(x+1,y) - coord->dx(x-1,y) ) / SQ(coord->dx(x,y)) ) * coef1; // BOUT-06 term } } - if(mesh->IncIntShear) { + if(localmesh->IncIntShear) { // d2dz2 term coef2 += coord->g11(x,y) * coord->IntShiftTorsion(x,y) * coord->IntShiftTorsion(x,y); // Mixed derivative @@ -976,8 +978,8 @@ void LaplacePetsc::Coeffs( int x, int y, int z, BoutReal &coef1, BoutReal &coef2 // A second/fourth order derivative term if (issetC) { -// if( (x > 0) && (x < (mesh->LocalNx-1)) ) //Valid if doing second order derivative, not if fourth: should only be called for xstart<=x<=xend anyway - if( (x > 1) && (x < (mesh->LocalNx-2)) ) { +// if( (x > 0) && (x < (localmesh->LocalNx-1)) ) //Valid if doing second order derivative, not if fourth: should only be called for xstart<=x<=xend anyway + if( (x > 1) && (x < (localmesh->LocalNx-2)) ) { int zp = z+1; // z plus 1 if (zp > meshz-1) zp -= meshz; int zm = z-1; // z minus 1 @@ -1024,13 +1026,16 @@ void LaplacePetsc::Coeffs( int x, int y, int z, BoutReal &coef1, BoutReal &coef2 void LaplacePetsc::vecToField(Vec xs, FieldPerp &f) { + + ASSERT1(localmesh == f.getMesh()); + f.allocate(); int i = Istart; - if(mesh->firstX()) + if(localmesh->firstX()) { - for(int x=0; xxstart; x++) + for(int x=0; xxstart; x++) { - for(int z=0; zLocalNz; z++) + for(int z=0; zLocalNz; z++) { PetscScalar val; VecGetValues(xs, 1, &i, &val ); @@ -1040,9 +1045,9 @@ void LaplacePetsc::vecToField(Vec xs, FieldPerp &f) { } } - for(int x=mesh->xstart; x <= mesh->xend; x++) + for(int x=localmesh->xstart; x <= localmesh->xend; x++) { - for(int z=0; zLocalNz; z++) + for(int z=0; zLocalNz; z++) { PetscScalar val; VecGetValues(xs, 1, &i, &val ); @@ -1051,11 +1056,11 @@ void LaplacePetsc::vecToField(Vec xs, FieldPerp &f) { } } - if(mesh->lastX()) + if(localmesh->lastX()) { - for(int x=mesh->xend+1; xLocalNx; x++) + for(int x=localmesh->xend+1; xLocalNx; x++) { - for(int z=0;z < mesh->LocalNz; z++) + for(int z=0;z < localmesh->LocalNz; z++) { PetscScalar val; VecGetValues(xs, 1, &i, &val ); @@ -1068,10 +1073,12 @@ void LaplacePetsc::vecToField(Vec xs, FieldPerp &f) { } void LaplacePetsc::fieldToVec(const FieldPerp &f, Vec bs) { + ASSERT1(localmesh == f.getMesh()); + int i = Istart; - if(mesh->firstX()) { - for(int x=0; xxstart; x++) { - for(int z=0; zLocalNz; z++) { + if(localmesh->firstX()) { + for(int x=0; xxstart; x++) { + for(int z=0; zLocalNz; z++) { PetscScalar val = f[x][z]; VecSetValues( bs, 1, &i, &val, INSERT_VALUES ); i++; // Increment row in Petsc matrix @@ -1079,17 +1086,17 @@ void LaplacePetsc::fieldToVec(const FieldPerp &f, Vec bs) { } } - for(int x=mesh->xstart; x <= mesh->xend; x++) { - for(int z=0; zLocalNz; z++) { + for(int x=localmesh->xstart; x <= localmesh->xend; x++) { + for(int z=0; zLocalNz; z++) { PetscScalar val = f[x][z]; VecSetValues( bs, 1, &i, &val, INSERT_VALUES ); i++; // Increment row in Petsc matrix } } - if(mesh->lastX()) { - for(int x=mesh->xend+1; xLocalNx; x++) { - for(int z=0;z < mesh->LocalNz; z++) { + if(localmesh->lastX()) { + for(int x=localmesh->xend+1; xLocalNx; x++) { + for(int z=0;z < localmesh->LocalNz; z++) { PetscScalar val = f[x][z]; VecSetValues( bs, 1, &i, &val, INSERT_VALUES ); i++; // Increment row in Petsc matrix diff --git a/src/invert/laplace/impls/petsc/petsc_laplace.hxx b/src/invert/laplace/impls/petsc/petsc_laplace.hxx index 919be782bc..61b094ace2 100644 --- a/src/invert/laplace/impls/petsc/petsc_laplace.hxx +++ b/src/invert/laplace/impls/petsc/petsc_laplace.hxx @@ -37,7 +37,7 @@ class LaplacePetsc; class LaplacePetsc : public Laplacian { public: - LaplacePetsc(Options *UNUSED(opt) = nullptr, const CELL_LOC UNUSED(loc) = CELL_CENTRE) { + LaplacePetsc(Options *UNUSED(opt) = nullptr, const CELL_LOC UNUSED(loc) = CELL_CENTRE, Mesh *UNUSED(mesh_in) = mesh) { throw BoutException("No PETSc solver available"); } @@ -68,7 +68,7 @@ public: class LaplacePetsc : public Laplacian { public: - LaplacePetsc(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE); + LaplacePetsc(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = mesh); ~LaplacePetsc() { KSPDestroy( &ksp ); VecDestroy( &xs ); @@ -78,12 +78,14 @@ public: void setCoefA(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); A = val; /*Acoefchanged = true;*/ if(pcsolve) pcsolve->setCoefA(val); } void setCoefC(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); C1 = val; C2 = val; issetC = true; /*coefchanged = true;*/ @@ -91,28 +93,33 @@ public: } void setCoefC1(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); C1 = val; issetC = true; } void setCoefC2(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); C2 = val; issetC = true; } void setCoefD(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); D = val; issetD = true; /*coefchanged = true;*/ if(pcsolve) pcsolve->setCoefD(val); } void setCoefEx(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); Ex = val; issetE = true; /*coefchanged = true;*/ if(pcsolve) pcsolve->setCoefEx(val); } void setCoefEz(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); Ez = val; issetE = true; /*coefchanged = true;*/ if(pcsolve) pcsolve->setCoefEz(val); @@ -120,12 +127,14 @@ public: void setCoefA(const Field3D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); A = val; /*Acoefchanged = true;*/ if(pcsolve) pcsolve->setCoefA(val); } void setCoefC(const Field3D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); C1 = val; C2 = val; issetC = true; /*coefchanged = true;*/ @@ -133,28 +142,33 @@ public: } void setCoefC1(const Field3D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); C1 = val; issetC = true; } void setCoefC2(const Field3D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); C2 = val; issetC = true; } void setCoefD(const Field3D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); D = val; issetD = true; /*coefchanged = true;*/ if(pcsolve) pcsolve->setCoefD(val); } void setCoefEx(const Field3D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); Ex = val; issetE = true; /*coefchanged = true;*/ if(pcsolve) pcsolve->setCoefEx(val); } void setCoefEz(const Field3D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh == val.getMesh()); Ez = val; issetE = true; /*coefchanged = true;*/ if(pcsolve) pcsolve->setCoefEz(val); diff --git a/src/invert/laplace/impls/serial_band/serial_band.cxx b/src/invert/laplace/impls/serial_band/serial_band.cxx index f2dfa04a66..b9541bf867 100644 --- a/src/invert/laplace/impls/serial_band/serial_band.cxx +++ b/src/invert/laplace/impls/serial_band/serial_band.cxx @@ -38,40 +38,41 @@ //#define SECONDORDER // Define to use 2nd order differencing -LaplaceSerialBand::LaplaceSerialBand(Options *opt, const CELL_LOC loc) : Laplacian(opt, loc), Acoef(0.0), Ccoef(1.0), Dcoef(1.0) { +LaplaceSerialBand::LaplaceSerialBand(Options *opt, const CELL_LOC loc, Mesh *mesh_in) + : Laplacian(opt, loc, mesh_in), Acoef(0.0), Ccoef(1.0), Dcoef(1.0) { Acoef.setLocation(location); Ccoef.setLocation(location); Dcoef.setLocation(location); - if(!mesh->firstX() || !mesh->lastX()) - throw BoutException("LaplaceSerialBand only works for mesh->NXPE = 1"); - if(mesh->periodicX) { - throw BoutException("LaplaceSerialBand does not work with periodicity in the x direction (mesh->PeriodicX == true). Change boundary conditions or use serial-tri or cyclic solver instead"); + if(!localmesh->firstX() || !localmesh->lastX()) + throw BoutException("LaplaceSerialBand only works for localmesh->NXPE = 1"); + if(localmesh->periodicX) { + throw BoutException("LaplaceSerialBand does not work with periodicity in the x direction (localmesh->PeriodicX == true). Change boundary conditions or use serial-tri or cyclic solver instead"); } // Allocate memory - int ncz = mesh->LocalNz; - bk = Matrix(mesh->LocalNx, ncz / 2 + 1); - bk1d = Array(mesh->LocalNx); + int ncz = localmesh->LocalNz; + bk = Matrix(localmesh->LocalNx, ncz / 2 + 1); + bk1d = Array(localmesh->LocalNx); //Initialise bk to 0 as we only visit 0<= kz <= maxmode in solve for(int kz=maxmode+1; kz < ncz/2 + 1; kz++){ - for (int ix=0; ixLocalNx; ix++){ + for (int ix=0; ixLocalNx; ix++){ bk(ix, kz) = 0.0; } } - xk = Matrix(mesh->LocalNx, ncz / 2 + 1); - xk1d = Array(mesh->LocalNx); + xk = Matrix(localmesh->LocalNx, ncz / 2 + 1); + xk1d = Array(localmesh->LocalNx); //Initialise xk to 0 as we only visit 0<= kz <= maxmode in solve for(int kz=maxmode+1; kz < ncz/2 + 1; kz++){ - for (int ix=0; ixLocalNx; ix++){ + for (int ix=0; ixLocalNx; ix++){ xk(ix, kz) = 0.0; } } - A = Matrix(mesh->LocalNx, 5); + A = Matrix(localmesh->LocalNx, 5); } const FieldPerp LaplaceSerialBand::solve(const FieldPerp &b) { @@ -79,25 +80,26 @@ const FieldPerp LaplaceSerialBand::solve(const FieldPerp &b) { } const FieldPerp LaplaceSerialBand::solve(const FieldPerp &b, const FieldPerp &x0) { - Mesh *mesh = b.getMesh(); - FieldPerp x(mesh); + ASSERT1(localmesh == b.getMesh() && localmesh == x0.getMesh()); + + FieldPerp x(localmesh); x.allocate(); int jy = b.getIndex(); x.setIndex(jy); - Coordinates *coord = mesh->getCoordinates(location); + Coordinates *coord = localmesh->getCoordinates(location); - int ncz = mesh->LocalNz; - int ncx = mesh->LocalNx-1; + int ncz = localmesh->LocalNz; + int ncx = localmesh->LocalNx-1; - int xbndry = mesh->xstart; // Width of the x boundary + int xbndry = localmesh->xstart; // Width of the x boundary // If the flags to assign that only one guard cell should be used is set - if((global_flags & INVERT_BOTH_BNDRY_ONE) || (mesh->xstart < 2)) + if((global_flags & INVERT_BOTH_BNDRY_ONE) || (localmesh->xstart < 2)) xbndry = 1; BOUT_OMP(parallel for) - for(int ix=0;ixLocalNx;ix++) { + for(int ix=0;ixLocalNx;ix++) { // for fixed ix,jy set a complex vector rho(z) if(((ix < xbndry) && (inner_boundary_flags & INVERT_SET)) || @@ -115,7 +117,7 @@ const FieldPerp LaplaceSerialBand::solve(const FieldPerp &b, const FieldPerp &x0 xend = ncx-xbndry; }else { xstart = 2; - xend = mesh->LocalNx-2; + xend = localmesh->LocalNx-2; } for(int iz=0;iz<=maxmode;iz++) { @@ -129,7 +131,7 @@ const FieldPerp LaplaceSerialBand::solve(const FieldPerp &b, const FieldPerp &x0 kwave=iz*2.0*PI/coord->zlength(); // wave number is 1/[rad] // set bk1d - for(int ix=0;ixLocalNx;ix++) + for(int ix=0;ixLocalNx;ix++) bk1d[ix] = bk(ix, iz); // Fill in interior points @@ -166,14 +168,14 @@ const FieldPerp LaplaceSerialBand::solve(const FieldPerp &b, const FieldPerp &x0 } if(nonuniform) { - // non-uniform mesh correction + // non-uniform localmesh correction if((ix != 0) && (ix != ncx)) coef4 += coord->g11(ix,jy)*( (1.0/coord->dx(ix+1,jy)) - (1.0/coord->dx(ix-1,jy)) )/(2.0*coord->dx(ix,jy)); } // A first order derivative term (1/c)\nabla_perp c\cdot\nabla_\perp x - if((ix > 1) && (ix < (mesh->LocalNx-2))) + if((ix > 1) && (ix < (localmesh->LocalNx-2))) coef4 += coord->g11(ix,jy) * (Ccoef(ix-2,jy) - 8.*Ccoef(ix-1,jy) + 8.*Ccoef(ix+1,jy) - Ccoef(ix+2,jy)) / (12.*coord->dx(ix,jy)*(Ccoef(ix,jy))); // Put into matrix @@ -392,7 +394,7 @@ const FieldPerp LaplaceSerialBand::solve(const FieldPerp &b, const FieldPerp &x0 } // Perform inversion - cband_solve(A, mesh->LocalNx, 2, 2, bk1d); + cband_solve(A, localmesh->LocalNx, 2, 2, bk1d); if((global_flags & INVERT_KX_ZERO) && (iz == 0)) { // Set the Kx = 0, n = 0 component to zero. For now just subtract diff --git a/src/invert/laplace/impls/serial_band/serial_band.hxx b/src/invert/laplace/impls/serial_band/serial_band.hxx index 7051b5925d..05a7649e32 100644 --- a/src/invert/laplace/impls/serial_band/serial_band.hxx +++ b/src/invert/laplace/impls/serial_band/serial_band.hxx @@ -36,22 +36,25 @@ class LaplaceSerialBand; class LaplaceSerialBand : public Laplacian { public: - LaplaceSerialBand(Options *opt = nullptr, const CELL_LOC = CELL_CENTRE); + LaplaceSerialBand(Options *opt = nullptr, const CELL_LOC = CELL_CENTRE, Mesh *mesh_in = mesh); ~LaplaceSerialBand(){}; using Laplacian::setCoefA; void setCoefA(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh = val.getMesh()); Acoef = val; } using Laplacian::setCoefC; void setCoefC(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh = val.getMesh()); Ccoef = val; } using Laplacian::setCoefD; void setCoefD(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh = val.getMesh()); Dcoef = val; } using Laplacian::setCoefEx; diff --git a/src/invert/laplace/impls/serial_tri/serial_tri.cxx b/src/invert/laplace/impls/serial_tri/serial_tri.cxx index 7080f18b31..0a2cbe02b1 100644 --- a/src/invert/laplace/impls/serial_tri/serial_tri.cxx +++ b/src/invert/laplace/impls/serial_tri/serial_tri.cxx @@ -37,13 +37,14 @@ #include -LaplaceSerialTri::LaplaceSerialTri(Options *opt, CELL_LOC loc) : Laplacian(opt, loc), A(0.0), C(1.0), D(1.0) { +LaplaceSerialTri::LaplaceSerialTri(Options *opt, CELL_LOC loc, Mesh *mesh_in) + : Laplacian(opt, loc, mesh_in), A(0.0), C(1.0), D(1.0) { A.setLocation(location); C.setLocation(location); D.setLocation(location); - if(!mesh->firstX() || !mesh->lastX()) { - throw BoutException("LaplaceSerialTri only works for mesh->NXPE = 1"); + if(!localmesh->firstX() || !localmesh->lastX()) { + throw BoutException("LaplaceSerialTri only works for localmesh->NXPE = 1"); } } @@ -72,24 +73,25 @@ const FieldPerp LaplaceSerialTri::solve(const FieldPerp &b) { * \return The inverted variable. */ const FieldPerp LaplaceSerialTri::solve(const FieldPerp &b, const FieldPerp &x0) { - Mesh *mesh = b.getMesh(); - FieldPerp x(mesh); + ASSERT1(localmesh == b.getMesh() && localmesh == x0.getMesh()); + + FieldPerp x(localmesh); x.allocate(); int jy = b.getIndex(); x.setIndex(jy); - int ncz = mesh->LocalNz; // No of z pnts - int ncx = mesh->LocalNx; // No of x pnts + int ncz = localmesh->LocalNz; // No of z pnts + int ncx = localmesh->LocalNx; // No of x pnts - BoutReal kwaveFactor = 2.0 * PI / mesh->getCoordinates(location)->zlength(); + BoutReal kwaveFactor = 2.0 * PI / localmesh->getCoordinates(location)->zlength(); // Setting the width of the boundary. // NOTE: The default is a width of 2 guard cells - int inbndry = mesh->xstart, outbndry=mesh->xstart; + int inbndry = localmesh->xstart, outbndry=localmesh->xstart; // If the flags to assign that only one guard cell should be used is set - if((global_flags & INVERT_BOTH_BNDRY_ONE) || (mesh->xstart < 2)) { + if((global_flags & INVERT_BOTH_BNDRY_ONE) || (localmesh->xstart < 2)) { inbndry = outbndry = 1; } if (inner_boundary_flags & INVERT_BNDRY_ONE) @@ -185,7 +187,7 @@ const FieldPerp LaplaceSerialTri::solve(const FieldPerp &b, const FieldPerp &x0) outer_boundary_flags, &A, &C, &D); ///////// PERFORM INVERSION ///////// - if (!mesh->periodicX) { + if (!localmesh->periodicX) { // Call tridiagonal solver tridag(std::begin(avec), std::begin(bvec), std::begin(cvec), std::begin(bk1d), std::begin(xk1d), ncx); @@ -193,7 +195,7 @@ const FieldPerp LaplaceSerialTri::solve(const FieldPerp &b, const FieldPerp &x0) } else { // Periodic in X, so cyclic tridiagonal - int xs = mesh->xstart; + int xs = localmesh->xstart; cyclic_tridag(&avec[xs], &bvec[xs], &cvec[xs], &bk1d[xs], &xk1d[xs], ncx - 2 * xs); // Copy boundary regions @@ -206,11 +208,11 @@ const FieldPerp LaplaceSerialTri::solve(const FieldPerp &b, const FieldPerp &x0) // If the global flag is set to INVERT_KX_ZERO if ((global_flags & INVERT_KX_ZERO) && (kz == 0)) { dcomplex offset(0.0); - for (int ix = mesh->xstart; ix <= mesh->xend; ix++) { + for (int ix = localmesh->xstart; ix <= localmesh->xend; ix++) { offset += xk1d[ix]; } - offset /= static_cast(mesh->xend - mesh->xstart + 1); - for (int ix = mesh->xstart; ix <= mesh->xend; ix++) { + offset /= static_cast(localmesh->xend - localmesh->xstart + 1); + for (int ix = localmesh->xstart; ix <= localmesh->xend; ix++) { xk1d[ix] -= offset; } } diff --git a/src/invert/laplace/impls/serial_tri/serial_tri.hxx b/src/invert/laplace/impls/serial_tri/serial_tri.hxx index 86e98e49f1..519dc9a361 100644 --- a/src/invert/laplace/impls/serial_tri/serial_tri.hxx +++ b/src/invert/laplace/impls/serial_tri/serial_tri.hxx @@ -35,22 +35,25 @@ class LaplaceSerialTri; class LaplaceSerialTri : public Laplacian { public: - LaplaceSerialTri(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE); + LaplaceSerialTri(Options *opt = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = mesh); ~LaplaceSerialTri(){}; using Laplacian::setCoefA; void setCoefA(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh = val.getMesh()); A = val; } using Laplacian::setCoefC; void setCoefC(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh = val.getMesh()); C = val; } using Laplacian::setCoefD; void setCoefD(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh = val.getMesh()); D = val; } using Laplacian::setCoefEx; diff --git a/src/invert/laplace/impls/shoot/shoot_laplace.cxx b/src/invert/laplace/impls/shoot/shoot_laplace.cxx index 0db36a3c65..86430f928c 100644 --- a/src/invert/laplace/impls/shoot/shoot_laplace.cxx +++ b/src/invert/laplace/impls/shoot/shoot_laplace.cxx @@ -36,23 +36,23 @@ #include #include -LaplaceShoot::LaplaceShoot(Options *opt, const CELL_LOC loc) - : Laplacian(opt, loc), Acoef(0.0), Ccoef(1.0), Dcoef(1.0) { +LaplaceShoot::LaplaceShoot(Options *opt, const CELL_LOC loc, Mesh *mesh_in) + : Laplacian(opt, loc, mesh_in), Acoef(0.0), Ccoef(1.0), Dcoef(1.0) { throw BoutException("LaplaceShoot is a test implementation and does not currently work. Please select a different implementation."); Acoef.setLocation(location); Ccoef.setLocation(location); Dcoef.setLocation(location); - if(mesh->periodicX) { - throw BoutException("LaplaceShoot does not work with periodicity in the x direction (mesh->PeriodicX == true). Change boundary conditions or use serial-tri or cyclic solver instead"); + if(localmesh->periodicX) { + throw BoutException("LaplaceShoot does not work with periodicity in the x direction (localmesh->PeriodicX == true). Change boundary conditions or use serial-tri or cyclic solver instead"); } nmode = maxmode + 1; // Number of Z modes. maxmode set in invert_laplace.cxx from options // Allocate memory - int size = (mesh->LocalNz)/2 + 1; + int size = (localmesh->LocalNz)/2 + 1; km = Array(size); kc = Array(size); kp = Array(size); @@ -69,20 +69,21 @@ LaplaceShoot::LaplaceShoot(Options *opt, const CELL_LOC loc) } const FieldPerp LaplaceShoot::solve(const FieldPerp &rhs) { - Mesh *mesh = rhs.getMesh(); - FieldPerp x(mesh); // Result + ASSERT1(localmesh = rhs.getMesh()); + + FieldPerp x(localmesh); // Result x.allocate(); int jy = rhs.getIndex(); // Get the Y index x.setIndex(jy); - Coordinates *coord = mesh->getCoordinates(location); + Coordinates *coord = localmesh->getCoordinates(location); // Get the width of the boundary - int inbndry = mesh->xstart, outbndry=mesh->xstart; + int inbndry = localmesh->xstart, outbndry=localmesh->xstart; // If the flags to assign that only one guard cell should be used is set - if((global_flags & INVERT_BOTH_BNDRY_ONE) || (mesh->xstart < 2)) { + if((global_flags & INVERT_BOTH_BNDRY_ONE) || (localmesh->xstart < 2)) { inbndry = outbndry = 1; } if(inner_boundary_flags & INVERT_BNDRY_ONE) @@ -91,14 +92,14 @@ const FieldPerp LaplaceShoot::solve(const FieldPerp &rhs) { outbndry = 1; int xs, xe; - xs = mesh->xstart; // Starting X index - if(mesh->firstX()) + xs = localmesh->xstart; // Starting X index + if(localmesh->firstX()) xs = inbndry; - xe = mesh->xend; // Last X index - if(mesh->lastX()) - xe = mesh->LocalNx-outbndry-1; + xe = localmesh->xend; // Last X index + if(localmesh->lastX()) + xe = localmesh->LocalNx-outbndry-1; - if(mesh->lastX()) { + if(localmesh->lastX()) { // Set initial value and gradient to zero // by setting kc and kp @@ -107,15 +108,15 @@ const FieldPerp LaplaceShoot::solve(const FieldPerp &rhs) { kp[i] = 0.0; } - for(int ix=xe;ixLocalNx;ix++) - for(int iz=0;izLocalNz;iz++) { + for(int ix=xe;ixLocalNx;ix++) + for(int iz=0;izLocalNz;iz++) { x(ix, iz) = 0.0; } }else { // Wait for processor outer X - comm_handle handle = mesh->irecvXOut(std::begin(buffer), 4 * maxmode, jy); - mesh->wait(handle); + comm_handle handle = localmesh->irecvXOut(std::begin(buffer), 4 * maxmode, jy); + localmesh->wait(handle); // Copy into kc, kp for(int i=0;iLocalNz, x[xe]); + irfft(std::begin(kc), localmesh->LocalNz, x[xe]); } // kc and kp now set to result at x and x+1 respectively // Use b at x to get km at x-1 // Loop inwards from edge for(int ix=xe; ix >= xs; ix--) { - rfft(rhs[ix], mesh->LocalNz, std::begin(rhsk)); + rfft(rhs[ix], localmesh->LocalNz, std::begin(rhsk)); for(int kz=0; kzzlength()); // wave number is 1/[rad] @@ -147,7 +148,7 @@ const FieldPerp LaplaceShoot::solve(const FieldPerp &rhs) { } // Inverse FFT to get x[ix-1] - irfft(std::begin(km), mesh->LocalNz, x[ix - 1]); + irfft(std::begin(km), localmesh->LocalNz, x[ix - 1]); // Cycle km->kc->kp std::swap(kp, kc); @@ -155,7 +156,7 @@ const FieldPerp LaplaceShoot::solve(const FieldPerp &rhs) { } // Finished on this processor. Send data to next inner processor - if(!mesh->firstX()) { + if(!localmesh->firstX()) { // Should be able to send dcomplex buffers. For now copy into BoutReal buffer for(int i=0;isendXIn(std::begin(buffer), 4 * maxmode, jy); + localmesh->sendXIn(std::begin(buffer), 4 * maxmode, jy); }else { // Set inner boundary for(int ix=xs-2;ix>=0;ix--) { - for(int iz=0;izLocalNz;iz++) { + for(int iz=0;izLocalNz;iz++) { x(ix, iz) = x(xs - 1, iz); } } diff --git a/src/invert/laplace/impls/shoot/shoot_laplace.hxx b/src/invert/laplace/impls/shoot/shoot_laplace.hxx index 999fb20ba2..db6e1bd617 100644 --- a/src/invert/laplace/impls/shoot/shoot_laplace.hxx +++ b/src/invert/laplace/impls/shoot/shoot_laplace.hxx @@ -37,22 +37,25 @@ class LaplaceShoot; class LaplaceShoot : public Laplacian { public: - LaplaceShoot(Options *opt = nullptr, const CELL_LOC = CELL_CENTRE); + LaplaceShoot(Options *opt = nullptr, const CELL_LOC = CELL_CENTRE, Mesh *mesh_in = mesh); ~LaplaceShoot(){}; using Laplacian::setCoefA; void setCoefA(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh = val.getMesh()); Acoef = val; } using Laplacian::setCoefC; void setCoefC(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh = val.getMesh()); Ccoef = val; } using Laplacian::setCoefD; void setCoefD(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh = val.getMesh()); Dcoef = val; } using Laplacian::setCoefEx; diff --git a/src/invert/laplace/impls/spt/spt.cxx b/src/invert/laplace/impls/spt/spt.cxx index b7867a87d5..e75cad6cae 100644 --- a/src/invert/laplace/impls/spt/spt.cxx +++ b/src/invert/laplace/impls/spt/spt.cxx @@ -41,23 +41,23 @@ #include "spt.hxx" -LaplaceSPT::LaplaceSPT(Options *opt, const CELL_LOC loc) - : Laplacian(opt, loc), Acoef(0.0), Ccoef(1.0), Dcoef(1.0) { +LaplaceSPT::LaplaceSPT(Options *opt, const CELL_LOC loc, Mesh *mesh_in) + : Laplacian(opt, loc, mesh_in), Acoef(0.0), Ccoef(1.0), Dcoef(1.0) { Acoef.setLocation(location); Ccoef.setLocation(location); Dcoef.setLocation(location); - if(mesh->periodicX) { - throw BoutException("LaplaceSPT does not work with periodicity in the x direction (mesh->PeriodicX == true). Change boundary conditions or use serial-tri or cyclic solver instead"); + if(localmesh->periodicX) { + throw BoutException("LaplaceSPT does not work with periodicity in the x direction (localmesh->PeriodicX == true). Change boundary conditions or use serial-tri or cyclic solver instead"); } // Get start and end indices - ys = mesh->ystart; - ye = mesh->yend; - if(mesh->hasBndryLowerY() && include_yguards) + ys = localmesh->ystart; + ye = localmesh->yend; + if(localmesh->hasBndryLowerY() && include_yguards) ys = 0; // Mesh contains a lower boundary - if(mesh->hasBndryUpperY() && include_yguards) - ye = mesh->LocalNy-1; // Contains upper boundary + if(localmesh->hasBndryUpperY() && include_yguards) + ye = localmesh->LocalNy-1; // Contains upper boundary alldata = new SPT_data[ye - ys + 1]; alldata -= ys; // Re-number indices to start at ys @@ -66,7 +66,7 @@ LaplaceSPT::LaplaceSPT(Options *opt, const CELL_LOC loc) } // Temporary array for taking FFTs - int ncz = mesh->LocalNz; + int ncz = localmesh->LocalNz; dc1d = Array(ncz / 2 + 1); } @@ -80,27 +80,28 @@ const FieldPerp LaplaceSPT::solve(const FieldPerp &b) { } const FieldPerp LaplaceSPT::solve(const FieldPerp &b, const FieldPerp &x0) { - Mesh *mesh = b.getMesh(); - FieldPerp x(mesh); + ASSERT1(localmesh == b.getMesh() && localmesh == x0.getMesh()); + + FieldPerp x(localmesh); x.allocate(); if( (inner_boundary_flags & INVERT_SET) || (outer_boundary_flags & INVERT_SET) ) { FieldPerp bs = copy(b); - int xbndry = mesh->xstart; + int xbndry = localmesh->xstart; // If the flags to assign that only one guard cell should be used is set - if((global_flags & INVERT_BOTH_BNDRY_ONE) || (mesh->xstart < 2)) + if((global_flags & INVERT_BOTH_BNDRY_ONE) || (localmesh->xstart < 2)) xbndry = 1; - if((inner_boundary_flags & INVERT_SET) && mesh->firstX()) { + if((inner_boundary_flags & INVERT_SET) && localmesh->firstX()) { // Copy x0 inner boundary into bs for(int ix=0;ixLocalNz;iz++) + for(int iz=0;izLocalNz;iz++) bs[ix][iz] = x0[ix][iz]; } - if((outer_boundary_flags & INVERT_SET) && mesh->lastX()) { + if((outer_boundary_flags & INVERT_SET) && localmesh->lastX()) { // Copy x0 outer boundary into bs - for(int ix=mesh->LocalNx-1;ix>=mesh->LocalNx-xbndry;ix--) - for(int iz=0;izLocalNz;iz++) + for(int ix=localmesh->LocalNx-1;ix>=localmesh->LocalNx-xbndry;ix--) + for(int iz=0;izLocalNz;iz++) bs[ix][iz] = x0[ix][iz]; } start(bs, slicedata); @@ -113,17 +114,17 @@ const FieldPerp LaplaceSPT::solve(const FieldPerp &b, const FieldPerp &x0) { /// Extracts perpendicular slices from 3D fields and inverts separately /*! - * In parallel (mesh->NXPE > 1) this tries to overlap computation and communication. + * In parallel (localmesh->NXPE > 1) this tries to overlap computation and communication. * This is done at the expense of more memory useage. Setting low_mem * in the config file uses less memory, and less communication overlap */ const Field3D LaplaceSPT::solve(const Field3D &b) { ASSERT1(b.getLocation() == location); + ASSERT1(localmesh = b.getMesh()); Timer timer("invert"); - Mesh *mesh = b.getMesh(); - Field3D x(mesh); + Field3D x(localmesh); x.allocate(); for(int jy=ys; jy <= ye; jy++) { @@ -142,7 +143,7 @@ const Field3D LaplaceSPT::solve(const Field3D &b) { running = next(alldata[jy]) == 0; }while(running); - FieldPerp xperp(mesh); + FieldPerp xperp(localmesh); xperp.allocate(); // All calculations finished. Get result @@ -157,27 +158,29 @@ const Field3D LaplaceSPT::solve(const Field3D &b) { } const Field3D LaplaceSPT::solve(const Field3D &b, const Field3D &x0) { - if( ((inner_boundary_flags & INVERT_SET) && mesh->firstX()) || - ((outer_boundary_flags & INVERT_SET) && mesh->lastX()) ) { + ASSERT1(localmesh == b.getMesh() && localmesh == x0.getMesh()); + + if( ((inner_boundary_flags & INVERT_SET) && localmesh->firstX()) || + ((outer_boundary_flags & INVERT_SET) && localmesh->lastX()) ) { Field3D bs = copy(b); - int xbndry = mesh->xstart; + int xbndry = localmesh->xstart; // If the flags to assign that only one guard cell should be used is set - if((global_flags & INVERT_BOTH_BNDRY_ONE) || (mesh->xstart < 2)) + if((global_flags & INVERT_BOTH_BNDRY_ONE) || (localmesh->xstart < 2)) xbndry = 1; - if((inner_boundary_flags & INVERT_SET) && mesh->firstX()) { + if((inner_boundary_flags & INVERT_SET) && localmesh->firstX()) { // Copy x0 inner boundary into bs for(int ix=0;ixLocalNy;iy++) - for(int iz=0;izLocalNz;iz++) + for(int iy=0;iyLocalNy;iy++) + for(int iz=0;izLocalNz;iz++) bs(ix,iy,iz) = x0(ix,iy,iz); } - if((outer_boundary_flags & INVERT_SET) && mesh->lastX()) { + if((outer_boundary_flags & INVERT_SET) && localmesh->lastX()) { // Copy x0 outer boundary into bs - for(int ix=mesh->LocalNx-1;ix>=mesh->LocalNx-xbndry;ix--) - for(int iy=0;iyLocalNy;iy++) - for(int iz=0;izLocalNz;iz++) + for(int ix=localmesh->LocalNx-1;ix>=localmesh->LocalNx-xbndry;ix--) + for(int iy=0;iyLocalNy;iy++) + for(int iz=0;izLocalNz;iz++) bs(ix,iy,iz) = x0(ix,iy,iz); } return solve(bs); @@ -235,8 +238,8 @@ void LaplaceSPT::tridagForward(dcomplex *a, dcomplex *b, dcomplex *c, * @param[inout] u Result to be solved (Au = r) * @param[in] n Size of the problem * @param[in] gam Intermediate values produced by the forward part - * @param[inout] gp gam from the processor mesh->PE_XIND + 1, and returned to mesh->PE_XIND - 1 - * @param[inout] up u from processor mesh->PE_XIND + 1, and returned to mesh->PE_XIND - 1 + * @param[inout] gp gam from the processor localmesh->PE_XIND + 1, and returned to localmesh->PE_XIND - 1 + * @param[inout] up u from processor localmesh->PE_XIND + 1, and returned to localmesh->PE_XIND - 1 */ void LaplaceSPT::tridagBack(dcomplex *u, int n, dcomplex *gam, dcomplex &gp, dcomplex &up) { @@ -256,39 +259,39 @@ void LaplaceSPT::tridagBack(dcomplex *u, int n, /// /// This is a reference code which performs the same operations as the /// serial code. To invert a single XZ slice (FieldPerp object), data -/// must pass from the innermost processor (mesh->PE_XIND = 0) to the -/// outermost (mesh->PE_XIND = mesh->NXPE-1) and back again. +/// must pass from the innermost processor (localmesh->PE_XIND = 0) to the +/// outermost (localmesh->PE_XIND = localmesh->NXPE-1) and back again. /// /// Some parallelism is achieved by running several inversions /// simultaneously, so while processor #1 is inverting Y=0, processor /// #0 is starting on Y=1. This works ok as long as the number of /// slices to be inverted is greater than the number of X processors -/// (MYSUB > mesh->NXPE). If MYSUB < mesh->NXPE then not all +/// (MYSUB > localmesh->NXPE). If MYSUB < localmesh->NXPE then not all /// processors can be busy at once, and so efficiency will fall /// sharply. /// /// @param[in] b RHS values (Ax = b) /// @param[out] data Structure containing data needed for second half of inversion int LaplaceSPT::start(const FieldPerp &b, SPT_data &data) { - if(mesh->firstX() && mesh->lastX()) - throw BoutException("Error: SPT method only works for mesh->NXPE > 1\n"); + if(localmesh->firstX() && localmesh->lastX()) + throw BoutException("Error: SPT method only works for localmesh->NXPE > 1\n"); data.jy = b.getIndex(); - int mm = mesh->LocalNz/2 + 1; - data.allocate(mm, mesh->LocalNx); // Make sure data is allocated. Already allocated -> does nothing + int mm = localmesh->LocalNz/2 + 1; + data.allocate(mm, localmesh->LocalNx); // Make sure data is allocated. Already allocated -> does nothing /// Take FFTs of data - int ncz = mesh->LocalNz; + int ncz = localmesh->LocalNz; - for(int ix=0; ix < mesh->LocalNx; ix++) { + for(int ix=0; ix < localmesh->LocalNx; ix++) { rfft(b[ix], ncz, std::begin(dc1d)); for(int kz = 0; kz <= maxmode; kz++) data.bk(kz, ix) = dc1d[kz]; } - BoutReal kwaveFactor = 2.0 * PI / mesh->getCoordinates(location)->zlength(); + BoutReal kwaveFactor = 2.0 * PI / localmesh->getCoordinates(location)->zlength(); /// Set matrix elements for (int kz = 0; kz <= maxmode; kz++) { @@ -300,13 +303,13 @@ int LaplaceSPT::start(const FieldPerp &b, SPT_data &data) { data.proc = 0; //< Starts at processor 0 data.dir = 1; - if(mesh->firstX()) { + if(localmesh->firstX()) { BOUT_OMP(parallel for) for(int kz = 0; kz <= maxmode; kz++) { dcomplex bet, u0; // Start tridiagonal solve tridagForward(&data.avec(kz, 0), &data.bvec(kz, 0), &data.cvec(kz, 0), - &data.bk(kz, 0), &data.xk(kz, 0), mesh->xend + 1, &data.gam(kz, 0), + &data.bk(kz, 0), &data.xk(kz, 0), localmesh->xend + 1, &data.gam(kz, 0), bet, u0, true); // Load intermediate values into buffers data.buffer[4*kz] = bet.real(); @@ -316,16 +319,16 @@ int LaplaceSPT::start(const FieldPerp &b, SPT_data &data) { } // Send data - mesh->sendXOut(std::begin(data.buffer), 4 * (maxmode + 1), data.comm_tag); + localmesh->sendXOut(std::begin(data.buffer), 4 * (maxmode + 1), data.comm_tag); - }else if(mesh->PE_XIND == 1) { + }else if(localmesh->PE_XIND == 1) { // Post a receive data.recv_handle = - mesh->irecvXIn(std::begin(data.buffer), 4 * (maxmode + 1), data.comm_tag); + localmesh->irecvXIn(std::begin(data.buffer), 4 * (maxmode + 1), data.comm_tag); } data.proc++; // Now moved onto the next processor - if(mesh->NXPE == 2) + if(localmesh->NXPE == 2) data.dir = -1; // Special case. Otherwise reversal handled in spt_continue return 0; @@ -341,13 +344,13 @@ int LaplaceSPT::next(SPT_data &data) { if(data.proc < 0) // Already finished return 1; - if(mesh->PE_XIND == data.proc) { + if(localmesh->PE_XIND == data.proc) { /// This processor's turn to do inversion // Wait for data to arrive - mesh->wait(data.recv_handle); + localmesh->wait(data.recv_handle); - if(mesh->lastX()) { + if(localmesh->lastX()) { // Last processor, turn-around BOUT_OMP(parallel for) @@ -356,16 +359,16 @@ int LaplaceSPT::next(SPT_data &data) { dcomplex gp, up; bet = dcomplex(data.buffer[4*kz], data.buffer[4*kz + 1]); u0 = dcomplex(data.buffer[4*kz + 2], data.buffer[4*kz + 3]); - tridagForward(&data.avec(kz, mesh->xstart), &data.bvec(kz, mesh->xstart), - &data.cvec(kz, mesh->xstart), &data.bk(kz, mesh->xstart), - &data.xk(kz, mesh->xstart), mesh->xend + 1, - &data.gam(kz, mesh->xstart), bet, u0); + tridagForward(&data.avec(kz, localmesh->xstart), &data.bvec(kz, localmesh->xstart), + &data.cvec(kz, localmesh->xstart), &data.bk(kz, localmesh->xstart), + &data.xk(kz, localmesh->xstart), localmesh->xend + 1, + &data.gam(kz, localmesh->xstart), bet, u0); // Back-substitute gp = 0.0; up = 0.0; - tridagBack(&data.xk(kz, mesh->xstart), mesh->LocalNx - mesh->xstart, - &data.gam(kz, mesh->xstart), gp, up); + tridagBack(&data.xk(kz, localmesh->xstart), localmesh->LocalNx - localmesh->xstart, + &data.gam(kz, localmesh->xstart), gp, up); data.buffer[4*kz] = gp.real(); data.buffer[4*kz + 1] = gp.imag(); data.buffer[4*kz + 2] = up.real(); @@ -380,10 +383,10 @@ int LaplaceSPT::next(SPT_data &data) { dcomplex bet, u0; bet = dcomplex(data.buffer[4*kz], data.buffer[4*kz + 1]); u0 = dcomplex(data.buffer[4*kz + 2], data.buffer[4*kz + 3]); - tridagForward(&data.avec(kz, mesh->xstart), &data.bvec(kz, mesh->xstart), - &data.cvec(kz, mesh->xstart), &data.bk(kz, mesh->xstart), - &data.xk(kz, mesh->xstart), mesh->xend - mesh->xstart + 1, - &data.gam(kz, mesh->xstart), bet, u0); + tridagForward(&data.avec(kz, localmesh->xstart), &data.bvec(kz, localmesh->xstart), + &data.cvec(kz, localmesh->xstart), &data.bk(kz, localmesh->xstart), + &data.xk(kz, localmesh->xstart), localmesh->xend - localmesh->xstart + 1, + &data.gam(kz, localmesh->xstart), bet, u0); // Load intermediate values into buffers data.buffer[4*kz] = bet.real(); data.buffer[4*kz + 1] = bet.imag(); @@ -391,7 +394,7 @@ int LaplaceSPT::next(SPT_data &data) { data.buffer[4*kz + 3] = u0.imag(); } - }else if(mesh->firstX()) { + }else if(localmesh->firstX()) { // Back to the start BOUT_OMP(parallel for) @@ -400,7 +403,7 @@ BOUT_OMP(parallel for) gp = dcomplex(data.buffer[4*kz], data.buffer[4*kz + 1]); up = dcomplex(data.buffer[4*kz + 2], data.buffer[4*kz + 3]); - tridagBack(&data.xk(kz, 0), mesh->xend + 1, &data.gam(kz, 0), gp, up); + tridagBack(&data.xk(kz, 0), localmesh->xend + 1, &data.gam(kz, 0), gp, up); } }else { @@ -411,8 +414,8 @@ BOUT_OMP(parallel for) dcomplex gp = dcomplex(data.buffer[4*kz], data.buffer[4*kz + 1]); dcomplex up = dcomplex(data.buffer[4*kz + 2], data.buffer[4*kz + 3]); - tridagBack(&data.xk(kz, mesh->xstart), mesh->xend - mesh->xstart + 1, - &data.gam(kz, mesh->xstart), gp, up); + tridagBack(&data.xk(kz, localmesh->xstart), localmesh->xend - localmesh->xstart + 1, + &data.gam(kz, localmesh->xstart), gp, up); data.buffer[4*kz] = gp.real(); data.buffer[4*kz + 1] = gp.imag(); @@ -421,29 +424,29 @@ BOUT_OMP(parallel for) } } - if(mesh->PE_XIND != 0) { // If not finished yet + if(localmesh->PE_XIND != 0) { // If not finished yet /// Send data if(data.dir > 0) { - mesh->sendXOut(std::begin(data.buffer), 4 * (maxmode + 1), data.comm_tag); + localmesh->sendXOut(std::begin(data.buffer), 4 * (maxmode + 1), data.comm_tag); }else - mesh->sendXIn(std::begin(data.buffer), 4 * (maxmode + 1), data.comm_tag); + localmesh->sendXIn(std::begin(data.buffer), 4 * (maxmode + 1), data.comm_tag); } - }else if(mesh->PE_XIND == data.proc + data.dir) { + }else if(localmesh->PE_XIND == data.proc + data.dir) { // This processor is next, post receive if(data.dir > 0) { data.recv_handle = - mesh->irecvXIn(std::begin(data.buffer), 4 * (maxmode + 1), data.comm_tag); + localmesh->irecvXIn(std::begin(data.buffer), 4 * (maxmode + 1), data.comm_tag); }else data.recv_handle = - mesh->irecvXOut(std::begin(data.buffer), 4 * (maxmode + 1), data.comm_tag); + localmesh->irecvXOut(std::begin(data.buffer), 4 * (maxmode + 1), data.comm_tag); } data.proc += data.dir; - if(data.proc == mesh->NXPE-1) + if(data.proc == localmesh->NXPE-1) data.dir = -1; // Reverses direction at the end return 0; @@ -454,8 +457,8 @@ BOUT_OMP(parallel for) /// @param[inout] data Structure keeping track of calculation /// @param[out] x The result void LaplaceSPT::finish(SPT_data &data, FieldPerp &x) { - int ncx = mesh->LocalNx-1; - int ncz = mesh->LocalNz; + int ncx = localmesh->LocalNx-1; + int ncz = localmesh->LocalNz; x.allocate(); x.setIndex(data.jy); @@ -479,17 +482,17 @@ void LaplaceSPT::finish(SPT_data &data, FieldPerp &x) { irfft(std::begin(dc1d), ncz, x[ix]); } - if(!mesh->firstX()) { + if(!localmesh->firstX()) { // Set left boundary to zero (Prevent unassigned values in corners) - for(int ix=0; ixxstart; ix++){ - for(int kz=0;kzLocalNz;kz++) + for(int ix=0; ixxstart; ix++){ + for(int kz=0;kzLocalNz;kz++) x(ix,kz) = 0.0; } } - if(!mesh->lastX()) { + if(!localmesh->lastX()) { // Same for right boundary - for(int ix=mesh->xend+1; ixLocalNx; ix++){ - for(int kz=0;kzLocalNz;kz++) + for(int ix=localmesh->xend+1; ixLocalNx; ix++){ + for(int kz=0;kzLocalNz;kz++) x(ix,kz) = 0.0; } } diff --git a/src/invert/laplace/impls/spt/spt.hxx b/src/invert/laplace/impls/spt/spt.hxx index 564f9c844c..44b7d98eb8 100644 --- a/src/invert/laplace/impls/spt/spt.hxx +++ b/src/invert/laplace/impls/spt/spt.hxx @@ -50,12 +50,12 @@ class LaplaceSPT; /*! * This is a reference code which performs the same operations as the serial code. * To invert a single XZ slice (FieldPerp object), data must pass from the innermost - * processor (mesh->PE_XIND = 0) to the outermost (mesh->PE_XIND = mesh->NXPE-1) and back again. + * processor (localmesh->PE_XIND = 0) to the outermost (localmesh->PE_XIND = localmesh->NXPE-1) and back again. * * Some parallelism is achieved by running several inversions simultaneously, so while * processor #1 is inverting Y=0, processor #0 is starting on Y=1. This works ok as long - * as the number of slices to be inverted is greater than the number of X processors (MYSUB > mesh->NXPE). - * If MYSUB < mesh->NXPE then not all processors can be busy at once, and so efficiency will fall sharply. + * as the number of slices to be inverted is greater than the number of X processors (MYSUB > localmesh->NXPE). + * If MYSUB < localmesh->NXPE then not all processors can be busy at once, and so efficiency will fall sharply. * * @param[in] b RHS values (Ax = b) * @param[in] flags Inversion settings (see boundary.h for values) @@ -66,22 +66,25 @@ class LaplaceSPT; */ class LaplaceSPT : public Laplacian { public: - LaplaceSPT(Options *opt = nullptr, const CELL_LOC = CELL_CENTRE); + LaplaceSPT(Options *opt = nullptr, const CELL_LOC = CELL_CENTRE, Mesh *mesh_in = mesh); ~LaplaceSPT(); using Laplacian::setCoefA; void setCoefA(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh = val.getMesh()); Acoef = val; } using Laplacian::setCoefC; void setCoefC(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh = val.getMesh()); Ccoef = val; } using Laplacian::setCoefD; void setCoefD(const Field2D &val) override { ASSERT1(val.getLocation() == location); + ASSERT1(localmesh = val.getMesh()); Dcoef = val; } using Laplacian::setCoefEx; diff --git a/src/invert/laplace/invert_laplace.cxx b/src/invert/laplace/invert_laplace.cxx index c9478eb20d..2e66fa3767 100644 --- a/src/invert/laplace/invert_laplace.cxx +++ b/src/invert/laplace/invert_laplace.cxx @@ -51,7 +51,8 @@ **********************************************************************************/ /// Laplacian inversion initialisation. Called once at the start to get settings -Laplacian::Laplacian(Options *options, const CELL_LOC loc) : location(loc) { +Laplacian::Laplacian(Options *options, const CELL_LOC loc, Mesh *mesh_in) + : location(loc), localmesh(mesh_in) { if (options == nullptr) { // Use the default options @@ -69,7 +70,7 @@ Laplacian::Laplacian(Options *options, const CELL_LOC loc) : location(loc) { BoutReal filter; ///< Fraction of Z modes to filter out. Between 0 and 1 OPTION(options, filter, 0.0); - int ncz = mesh->LocalNz; + int ncz = localmesh->LocalNz; // convert filtering into an integer number of modes maxmode = ROUND((1.0 - filter) * static_cast(ncz / 2)); // Can be overriden by max_mode option @@ -80,7 +81,7 @@ Laplacian::Laplacian(Options *options, const CELL_LOC loc) : location(loc) { OPTION(options, low_mem, false); OPTION(options, nonuniform, - mesh->getCoordinates(location)->non_uniform); // Default is the mesh setting + localmesh->getCoordinates(location)->non_uniform); // Default is the mesh setting OPTION(options, all_terms, true); // Include first derivative terms @@ -103,13 +104,13 @@ Laplacian::Laplacian(Options *options, const CELL_LOC loc) : location(loc) { OPTION2(options, extra_yguards_lower, extra_yguards_upper, 0); } -Laplacian* Laplacian::create(Options *opts, const CELL_LOC location) { +Laplacian* Laplacian::create(Options *opts, const CELL_LOC location, Mesh *mesh_in) { // Factory pattern: // 1. getInstance() is making an instance of LaplacianFactory // 2. createLaplacian() is accessing this instance and returning a Laplacian // form one of the child classes of the Laplacian (the laplace solver // implementations) - return LaplaceFactory::getInstance()->createLaplacian(opts, location); + return LaplaceFactory::getInstance()->createLaplacian(opts, location, mesh_in); } Laplacian *Laplacian::instance = nullptr; @@ -135,26 +136,25 @@ const Field3D Laplacian::solve(const Field3D &b) { TRACE("Laplacian::solve(Field3D)"); ASSERT1(b.getLocation() == location); - - Mesh *mesh = b.getMesh(); + ASSERT1(localmesh = b.getMesh()); Timer timer("invert"); - int ys = mesh->ystart, ye = mesh->yend; + int ys = localmesh->ystart, ye = localmesh->yend; - if(mesh->hasBndryLowerY()) { + if(localmesh->hasBndryLowerY()) { if (include_yguards) ys = 0; // Mesh contains a lower boundary and we are solving in the guard cells ys += extra_yguards_lower; } - if(mesh->hasBndryUpperY()) { + if(localmesh->hasBndryUpperY()) { if (include_yguards) - ye = mesh->LocalNy-1; // Contains upper boundary and we are solving in the guard cells + ye = localmesh->LocalNy-1; // Contains upper boundary and we are solving in the guard cells ye -= extra_yguards_upper; } - Field3D x(mesh); + Field3D x(localmesh); x.allocate(); int status = 0; @@ -198,18 +198,18 @@ const Field3D Laplacian::solve(const Field3D &b, const Field3D &x0) { ASSERT1(b.getLocation() == location); ASSERT1(x0.getLocation() == location); + ASSERT1(localmesh == b.getMesh() && localmesh == x0.getMesh()); Timer timer("invert"); - Mesh *mesh = b.getMesh(); // Setting the start and end range of the y-slices - int ys = mesh->ystart, ye = mesh->yend; - if(mesh->hasBndryLowerY() && include_yguards) + int ys = localmesh->ystart, ye = localmesh->yend; + if(localmesh->hasBndryLowerY() && include_yguards) ys = 0; // Mesh contains a lower boundary - if(mesh->hasBndryUpperY() && include_yguards) - ye = mesh->LocalNy-1; // Contains upper boundary + if(localmesh->hasBndryUpperY() && include_yguards) + ye = localmesh->LocalNy-1; // Contains upper boundary - Field3D x(mesh); + Field3D x(localmesh); x.allocate(); int status = 0; @@ -249,7 +249,7 @@ void Laplacian::tridagCoefs(int jx, int jy, int jz, ASSERT1(ccoef == nullptr || ccoef->getLocation() == loc); ASSERT1(d == nullptr || d->getLocation() == loc); - Coordinates *coord = mesh->getCoordinates(loc); + Coordinates *coord = localmesh->getCoordinates(loc); BoutReal kwave=jz*2.0*PI/coord->zlength(); // wave number is 1/[rad] @@ -275,7 +275,7 @@ void Laplacian::tridagCoefs(int jx, int jy, BoutReal kwave, * Input: * jx - The current x index * jy - The current y index - * kwave - The mode number multiplied with (2*pi)/mesh->zlength(), where + * kwave - The mode number multiplied with (2*pi)/localmesh->zlength(), where * zlength() is the length of the full z domain (usually 2*pi) * a - Lower diagonal of the tridiagonal matrix. DO NOT CONFUSE WITH A * b - The main diagonal @@ -296,7 +296,7 @@ void Laplacian::tridagCoefs(int jx, int jy, BoutReal kwave, BoutReal coef1, coef2, coef3, coef4, coef5; - Coordinates *coord = mesh->getCoordinates(loc); + Coordinates *coord = localmesh->getCoordinates(loc); coef1=coord->g11(jx,jy); ///< X 2nd derivative coefficient coef2=coord->g33(jx,jy); ///< Z 2nd derivative coefficient @@ -321,18 +321,18 @@ void Laplacian::tridagCoefs(int jx, int jy, BoutReal kwave, if(nonuniform) { // non-uniform mesh correction - if((jx != 0) && (jx != (mesh->LocalNx-1))) { + if((jx != 0) && (jx != (localmesh->LocalNx-1))) { coef4 -= 0.5*((coord->dx(jx+1,jy) - coord->dx(jx-1,jy))/SQ(coord->dx(jx,jy)))*coef1; } } if (ccoef != nullptr) { // A first order derivative term - if((jx > 0) && (jx < (mesh->LocalNx-1))) + if((jx > 0) && (jx < (localmesh->LocalNx-1))) coef4 += coord->g11(jx,jy) * ((*ccoef)(jx+1,jy) - (*ccoef)(jx-1,jy)) / (2.*coord->dx(jx,jy)*((*ccoef)(jx,jy))); } - if(mesh->IncIntShear) { + if(localmesh->IncIntShear) { // d2dz2 term coef2 += coord->g11(jx,jy) * coord->IntShiftTorsion(jx,jy) * coord->IntShiftTorsion(jx,jy); // Mixed derivative @@ -361,7 +361,7 @@ void Laplacian::tridagMatrix(dcomplex **avec, dcomplex **bvec, dcomplex **cvec, ASSERT1(ccoef->getLocation() == location); ASSERT1(d->getLocation() == location); - Coordinates *coord = mesh->getCoordinates(location); + Coordinates *coord = localmesh->getCoordinates(location); BOUT_OMP(parallel for) for(int kz = 0; kz <= maxmode; kz++) { @@ -424,27 +424,27 @@ void Laplacian::tridagMatrix(dcomplex *avec, dcomplex *bvec, dcomplex *cvec, ASSERT1(d->getLocation() == location); int xs = 0; // xstart set to the start of x on this processor (including ghost points) - int xe = mesh->LocalNx-1; // xend set to the end of x on this processor (including ghost points) + int xe = localmesh->LocalNx-1; // xend set to the end of x on this processor (including ghost points) - Coordinates *coord = mesh->getCoordinates(location); + Coordinates *coord = localmesh->getCoordinates(location); // Do not want boundary cells if x is periodic for cyclic solver. Only other solver which // works with periodicX is serial_tri, which uses includeguards==true, so the below isn't called. if(!includeguards) { - if(!mesh->firstX() || mesh->periodicX) - xs = mesh->xstart; // Inner edge is a guard cell - if(!mesh->lastX() || mesh->periodicX) - xe = mesh->xend; // Outer edge is a guard cell + if(!localmesh->firstX() || localmesh->periodicX) + xs = localmesh->xstart; // Inner edge is a guard cell + if(!localmesh->lastX() || localmesh->periodicX) + xe = localmesh->xend; // Outer edge is a guard cell } int ncx = xe - xs; // Total number of points in x to be used // Setting the width of the boundary. - // NOTE: The default is a width of (mesh->xstart) guard cells - int inbndry = mesh->xstart, outbndry=mesh->xstart; + // NOTE: The default is a width of (localmesh->xstart) guard cells + int inbndry = localmesh->xstart, outbndry=localmesh->xstart; // If the flags to assign that only one guard cell should be used is set - if((global_flags & INVERT_BOTH_BNDRY_ONE) || (mesh->xstart < 2)) { + if((global_flags & INVERT_BOTH_BNDRY_ONE) || (localmesh->xstart < 2)) { inbndry = outbndry = 1; } if(inner_boundary_flags & INVERT_BNDRY_ONE) @@ -463,8 +463,8 @@ void Laplacian::tridagMatrix(dcomplex *avec, dcomplex *bvec, dcomplex *cvec, } // Set the boundary conditions if x is not periodic - if(!mesh->periodicX) { - if(mesh->firstX()) { + if(!localmesh->periodicX) { + if(localmesh->firstX()) { // INNER BOUNDARY ON THIS PROCESSOR // If no user specified value is set on inner boundary, set the first @@ -627,7 +627,7 @@ void Laplacian::tridagMatrix(dcomplex *avec, dcomplex *bvec, dcomplex *cvec, } } } - if(mesh->lastX()) { + if(localmesh->lastX()) { // OUTER BOUNDARY ON THIS PROCESSOR // If no user specified value is set on outer boundary, set the last diff --git a/src/invert/laplace/laplacefactory.cxx b/src/invert/laplace/laplacefactory.cxx index 6a04e6e68e..5446426508 100644 --- a/src/invert/laplace/laplacefactory.cxx +++ b/src/invert/laplace/laplacefactory.cxx @@ -37,35 +37,35 @@ LaplaceFactory* LaplaceFactory::getInstance() { return instance; } -Laplacian* LaplaceFactory::createLaplacian(Options *options, const CELL_LOC loc) { +Laplacian* LaplaceFactory::createLaplacian(Options *options, const CELL_LOC loc, Mesh *mesh_in) { if (options == nullptr) options = Options::getRoot()->getSection("laplace"); string type; - if(mesh->firstX() && mesh->lastX()) { + if(mesh_in->firstX() && mesh_in->lastX()) { // Can use serial algorithm options->get("type", type, LAPLACE_CYCLIC); if(strcasecmp(type.c_str(), LAPLACE_TRI) == 0) { - return new LaplaceSerialTri(options, loc); + return new LaplaceSerialTri(options, loc, mesh_in); }else if(strcasecmp(type.c_str(), LAPLACE_BAND) == 0) { - return new LaplaceSerialBand(options, loc); + return new LaplaceSerialBand(options, loc, mesh_in); }else if(strcasecmp(type.c_str(), LAPLACE_SPT) == 0) { - return new LaplaceSPT(options, loc); + return new LaplaceSPT(options, loc, mesh_in); }else if(strcasecmp(type.c_str(), LAPLACE_PETSC) == 0) { - return new LaplacePetsc(options, loc); + return new LaplacePetsc(options, loc, mesh_in); }else if(strcasecmp(type.c_str(), LAPLACE_MUMPS) == 0) { - return new LaplaceMumps(options, loc); + return new LaplaceMumps(options, loc, mesh_in); }else if(strcasecmp(type.c_str(), LAPLACE_CYCLIC) == 0) { - return new LaplaceCyclic(options, loc); + return new LaplaceCyclic(options, loc, mesh_in); }else if(strcasecmp(type.c_str(), LAPLACE_SHOOT) == 0) { - return new LaplaceShoot(options, loc); + return new LaplaceShoot(options, loc, mesh_in); }else if(strcasecmp(type.c_str(), LAPLACE_MULTIGRID) == 0) { - return new LaplaceMultigrid(options, loc); + return new LaplaceMultigrid(options, loc, mesh_in); }else if(strcasecmp(type.c_str(), LAPLACE_NAULIN) == 0) { - return new LaplaceNaulin(options, loc); + return new LaplaceNaulin(options, loc, mesh_in); }else { throw BoutException("Unknown serial Laplacian solver type '%s'", type.c_str()); } @@ -75,21 +75,21 @@ Laplacian* LaplaceFactory::createLaplacian(Options *options, const CELL_LOC loc) // Parallel algorithm if(strcasecmp(type.c_str(), LAPLACE_PDD) == 0) { - return new LaplacePDD(options, loc); + return new LaplacePDD(options, loc, mesh_in); }else if(strcasecmp(type.c_str(), LAPLACE_SPT) == 0) { - return new LaplaceSPT(options, loc); + return new LaplaceSPT(options, loc, mesh_in); }else if(strcasecmp(type.c_str(), LAPLACE_PETSC) == 0) { - return new LaplacePetsc(options, loc); + return new LaplacePetsc(options, loc, mesh_in); }else if(strcasecmp(type.c_str(), LAPLACE_MUMPS) == 0) { - return new LaplaceMumps(options, loc); + return new LaplaceMumps(options, loc, mesh_in); }else if(strcasecmp(type.c_str(), LAPLACE_CYCLIC) == 0) { - return new LaplaceCyclic(options, loc); + return new LaplaceCyclic(options, loc, mesh_in); }else if(strcasecmp(type.c_str(), LAPLACE_SHOOT) == 0) { - return new LaplaceShoot(options, loc); + return new LaplaceShoot(options, loc, mesh_in); }else if(strcasecmp(type.c_str(), LAPLACE_MULTIGRID) == 0) { - return new LaplaceMultigrid(options, loc); + return new LaplaceMultigrid(options, loc, mesh_in); }else if(strcasecmp(type.c_str(), LAPLACE_NAULIN) == 0) { - return new LaplaceNaulin(options, loc); + return new LaplaceNaulin(options, loc, mesh_in); }else { throw BoutException("Unknown parallel Laplacian solver type '%s'", type.c_str()); } diff --git a/src/invert/laplace/laplacefactory.hxx b/src/invert/laplace/laplacefactory.hxx index e71543a073..215ee8790d 100644 --- a/src/invert/laplace/laplacefactory.hxx +++ b/src/invert/laplace/laplacefactory.hxx @@ -11,7 +11,7 @@ class LaplaceFactory { /// Return a pointer to the only instance static LaplaceFactory* getInstance(); - Laplacian *createLaplacian(Options *options = nullptr, const CELL_LOC loc = CELL_CENTRE); + Laplacian *createLaplacian(Options *options = nullptr, const CELL_LOC loc = CELL_CENTRE, Mesh *mesh_in = mesh); private: LaplaceFactory() {} // Prevent instantiation of this class diff --git a/src/invert/parderiv/impls/cyclic/cyclic.cxx b/src/invert/parderiv/impls/cyclic/cyclic.cxx index f6d0a3eb2f..2dfc7c2de3 100644 --- a/src/invert/parderiv/impls/cyclic/cyclic.cxx +++ b/src/invert/parderiv/impls/cyclic/cyclic.cxx @@ -48,18 +48,19 @@ #include -InvertParCR::InvertParCR(Options *opt) : InvertPar(opt), A(1.0), B(0.0), C(0.0), D(0.0), E(0.0) { +InvertParCR::InvertParCR(Options *opt, Mesh *mesh_in) + : InvertPar(opt, mesh_in), A(1.0), B(0.0), C(0.0), D(0.0), E(0.0) { // Number of k equations to solve for each x location - nsys = 1 + (mesh->LocalNz)/2; + nsys = 1 + (localmesh->LocalNz)/2; - rhs = Matrix(mesh->LocalNy, nsys); + rhs = Matrix(localmesh->LocalNy, nsys); // Find out if we are on a boundary - int size = mesh->LocalNy-4; - SurfaceIter surf(mesh); + int size = localmesh->LocalNy-4; + SurfaceIter surf(localmesh); for(surf.first(); !surf.isDone(); surf.next()) { BoutReal ts; - int n = mesh->LocalNy-4; + int n = localmesh->LocalNy-4; if(!surf.closed(ts)) { // Open field line if(surf.firstY()) @@ -84,8 +85,9 @@ InvertParCR::~InvertParCR() { const Field3D InvertParCR::solve(const Field3D &f) { TRACE("InvertParCR::solve(Field3D)"); - Mesh *mesh = f.getMesh(); - Field3D result(mesh); + ASSERT1(localmesh == f.getMesh()); + + Field3D result(localmesh); result.allocate(); result.setLocation(f.getLocation()); @@ -96,7 +98,7 @@ const Field3D InvertParCR::solve(const Field3D &f) { new CyclicReduce(); // Loop over flux-surfaces - SurfaceIter surf(mesh); + SurfaceIter surf(localmesh); for(surf.first(); !surf.isDone(); surf.next()) { int x = surf.xpos; @@ -105,7 +107,7 @@ const Field3D InvertParCR::solve(const Field3D &f) { bool closed = surf.closed(ts); // Number of rows - int y0 = 0, size = mesh->LocalNy-4; // If no boundaries + int y0 = 0, size = localmesh->LocalNy-4; // If no boundaries if(!closed) { if(surf.firstY()) { y0 += 2; @@ -120,13 +122,13 @@ const Field3D InvertParCR::solve(const Field3D &f) { cr->setPeriodic(closed); // Take Fourier transform - for(int y=0;yLocalNy-4;y++) - rfft(f(x, y + 2), mesh->LocalNz, &rhs(y + y0, 0)); + for(int y=0;yLocalNy-4;y++) + rfft(f(x, y + 2), localmesh->LocalNz, &rhs(y + y0, 0)); // Set up tridiagonal system for(int k=0; kzlength(); // wave number is 1/[rad] - for(int y=0;yLocalNy-4;y++) { + for(int y=0;yLocalNy-4;y++) { BoutReal acoef = A(x, y+2); // Constant BoutReal bcoef = B(x, y+2) / coord->g_22(x,y+2); // d2dy2 @@ -165,7 +167,7 @@ const Field3D InvertParCR::solve(const Field3D &f) { for(int k=0; kzlength(); // wave number is 1/[rad] dcomplex phase(cos(kwave*ts) , sin(kwave*ts)); - c(k, mesh->LocalNy - 5) *= phase; + c(k, localmesh->LocalNy - 5) *= phase; } } }else { @@ -206,7 +208,7 @@ const Field3D InvertParCR::solve(const Field3D &f) { // Inverse Fourier transform for(int y=0;yLocalNz, result(x, y + 2 - y0)); + irfft(&rhs(y, 0), localmesh->LocalNz, result(x, y + 2 - y0)); } // Delete cyclic reduction object diff --git a/src/invert/parderiv/impls/cyclic/cyclic.hxx b/src/invert/parderiv/impls/cyclic/cyclic.hxx index 2de04d0575..de4758315c 100644 --- a/src/invert/parderiv/impls/cyclic/cyclic.hxx +++ b/src/invert/parderiv/impls/cyclic/cyclic.hxx @@ -45,22 +45,37 @@ class InvertParCR : public InvertPar { public: - InvertParCR(Options *opt); + InvertParCR(Options *opt, Mesh *mesh_in = mesh); ~InvertParCR(); using InvertPar::solve; const Field3D solve(const Field3D &f) override; using InvertPar::setCoefA; - void setCoefA(const Field2D &f) override { A = f; } + void setCoefA(const Field2D &f) override { + ASSERT1(localmesh == f.getMesh()); + A = f; + } using InvertPar::setCoefB; - void setCoefB(const Field2D &f) override { B = f; } + void setCoefB(const Field2D &f) override { + ASSERT1(localmesh == f.getMesh()); + B = f; + } using InvertPar::setCoefC; - void setCoefC(const Field2D &f) override { C = f; } + void setCoefC(const Field2D &f) override { + ASSERT1(localmesh == f.getMesh()); + C = f; + } using InvertPar::setCoefD; - void setCoefD(const Field2D &f) override { D = f; } + void setCoefD(const Field2D &f) override { + ASSERT1(localmesh == f.getMesh()); + D = f; + } using InvertPar::setCoefE; - void setCoefE(const Field2D &f) override { E = f; } + void setCoefE(const Field2D &f) override { + ASSERT1(localmesh == f.getMesh()); + E = f; + } private: Field2D A, B, C, D, E; diff --git a/src/invert/parderiv/impls/serial/serial.cxx b/src/invert/parderiv/impls/serial/serial.cxx index 200650efd2..65b56ed8e8 100644 --- a/src/invert/parderiv/impls/serial/serial.cxx +++ b/src/invert/parderiv/impls/serial/serial.cxx @@ -51,26 +51,29 @@ #include -InvertParSerial::InvertParSerial(Options *opt) : InvertPar(opt), A(1.0), B(0.0), C(0.0), D(0.0), E(0.0) { - rhs = Matrix(mesh->LocalNy, (mesh->LocalNz)/2 + 1); - rhsk = Array(mesh->LocalNy-4); - xk = Array(mesh->LocalNy-4); - a = Array(mesh->LocalNy-4); - b = Array(mesh->LocalNy-4); - c = Array(mesh->LocalNy-4); +InvertParSerial::InvertParSerial(Options *opt, Mesh *mesh_in) + : InvertPar(opt, mesh_in), A(1.0), B(0.0), C(0.0), D(0.0), E(0.0) { + rhs = Matrix(localmesh->LocalNy, (localmesh->LocalNz)/2 + 1); + rhsk = Array(localmesh->LocalNy-4); + xk = Array(localmesh->LocalNy-4); + a = Array(localmesh->LocalNy-4); + b = Array(localmesh->LocalNy-4); + c = Array(localmesh->LocalNy-4); } const Field3D InvertParSerial::solve(const Field3D &f) { TRACE("InvertParSerial::solve(Field3D)"); - Field3D result(f.getMesh()); + ASSERT1(localmesh == f.getMesh()); + + Field3D result(localmesh); result.allocate(); result.setLocation(f.getLocation()); Coordinates *coord = f.getCoordinates(); // Loop over flux-surfaces - SurfaceIter surf(mesh); + SurfaceIter surf(localmesh); for(surf.first(); !surf.isDone(); surf.next()) { int x = surf.xpos; BoutReal ts; // Twist-shift angle @@ -78,20 +81,20 @@ const Field3D InvertParSerial::solve(const Field3D &f) { throw BoutException("InvertParSerial doesn't handle open surfaces"); // Take Fourier transform - for(int y=0;yLocalNy-4;y++) - rfft(f(x,y+2), mesh->LocalNz, &rhs(y, 0)); + for(int y=0;yLocalNy-4;y++) + rfft(f(x,y+2), localmesh->LocalNz, &rhs(y, 0)); // Solve cyclic tridiagonal system for each k - int nyq = (mesh->LocalNz)/2; + int nyq = (localmesh->LocalNz)/2; for(int k=0;k<=nyq;k++) { // Copy component of rhs into 1D array - for(int y=0;yLocalNy-4;y++) + for(int y=0;yLocalNy-4;y++) rhsk[y] = rhs(y, k); BoutReal kwave=k*2.0*PI/coord->zlength(); // wave number is 1/[rad] // Set up tridiagonal system - for(int y=0;yLocalNy-4;y++) { + for(int y=0;yLocalNy-4;y++) { BoutReal acoef = A(x, y+2); // Constant BoutReal bcoef = B(x, y+2) / coord->g_22(x,y+2); // d2dy2 BoutReal ccoef = C(x, y+2); // d2dydz @@ -113,19 +116,19 @@ const Field3D InvertParSerial::solve(const Field3D &f) { // Modify coefficients across twist-shift dcomplex phase(cos(kwave*ts) , -sin(kwave*ts)); a[0] *= phase; - c[mesh->LocalNy-5] /= phase; + c[localmesh->LocalNy-5] /= phase; // Solve cyclic tridiagonal system - cyclic_tridag(std::begin(a), std::begin(b), std::begin(c), std::begin(rhsk), std::begin(xk), mesh->LocalNy-4); + cyclic_tridag(std::begin(a), std::begin(b), std::begin(c), std::begin(rhsk), std::begin(xk), localmesh->LocalNy-4); // Put back into rhs array - for(int y=0;yLocalNy-4;y++) + for(int y=0;yLocalNy-4;y++) rhs(y, k) = xk[y]; } // Inverse Fourier transform - for(int y=0;yLocalNy-4;y++) - irfft(&rhs(y, 0), mesh->LocalNz, result(x,y+2)); + for(int y=0;yLocalNy-4;y++) + irfft(&rhs(y, 0), localmesh->LocalNz, result(x,y+2)); } return result; diff --git a/src/invert/parderiv/impls/serial/serial.hxx b/src/invert/parderiv/impls/serial/serial.hxx index a009c966c3..547c5be5d4 100644 --- a/src/invert/parderiv/impls/serial/serial.hxx +++ b/src/invert/parderiv/impls/serial/serial.hxx @@ -47,22 +47,37 @@ class InvertParSerial : public InvertPar { public: - InvertParSerial(Options* opt); + InvertParSerial(Options* opt, Mesh* mesh_in = mesh); ~InvertParSerial(){}; using InvertPar::solve; const Field3D solve(const Field3D &f) override; using InvertPar::setCoefA; - void setCoefA(const Field2D &f) override { A = f; } + void setCoefA(const Field2D &f) override { + ASSERT1(localmesh == f.getMesh()); + A = f; + } using InvertPar::setCoefB; - void setCoefB(const Field2D &f) override { B = f; } + void setCoefB(const Field2D &f) override { + ASSERT1(localmesh == f.getMesh()); + B = f; + } using InvertPar::setCoefC; - void setCoefC(const Field2D &f) override { C = f; } + void setCoefC(const Field2D &f) override { + ASSERT1(localmesh == f.getMesh()); + C = f; + } using InvertPar::setCoefD; - void setCoefD(const Field2D &f) override { D = f; } + void setCoefD(const Field2D &f) override { + ASSERT1(localmesh == f.getMesh()); + D = f; + } using InvertPar::setCoefE; - void setCoefE(const Field2D &f) override { E = f; } + void setCoefE(const Field2D &f) override { + ASSERT1(localmesh == f.getMesh()); + E = f; + } private: Field2D A, B, C, D, E; diff --git a/src/invert/parderiv/invert_parderiv.cxx b/src/invert/parderiv/invert_parderiv.cxx index 1938a383b6..9a06c9449e 100644 --- a/src/invert/parderiv/invert_parderiv.cxx +++ b/src/invert/parderiv/invert_parderiv.cxx @@ -30,8 +30,8 @@ #include #include "parderiv_factory.hxx" -InvertPar* InvertPar::Create() { - return ParDerivFactory::getInstance()->createInvertPar(); +InvertPar* InvertPar::Create(Mesh* mesh_in) { + return ParDerivFactory::getInstance()->createInvertPar(mesh_in); } const Field2D InvertPar::solve(const Field2D &f) { diff --git a/src/invert/parderiv/parderiv_factory.cxx b/src/invert/parderiv/parderiv_factory.cxx index 1408dd8960..2fa681e08f 100644 --- a/src/invert/parderiv/parderiv_factory.cxx +++ b/src/invert/parderiv/parderiv_factory.cxx @@ -23,14 +23,14 @@ ParDerivFactory* ParDerivFactory::getInstance() { return instance; } -InvertPar* ParDerivFactory::createInvertPar() { +InvertPar* ParDerivFactory::createInvertPar(Mesh *mesh_in) { // Get the default options section Options *opt = Options::getRoot()->getSection(default_section); - return createInvertPar( opt ); + return createInvertPar(opt, mesh_in); } -InvertPar* ParDerivFactory::createInvertPar(const char* type, Options *opt) { +InvertPar* ParDerivFactory::createInvertPar(const char* type, Options *opt, Mesh *mesh_in) { int NPES; MPI_Comm_size(BoutComm::get(), &NPES); @@ -38,17 +38,17 @@ InvertPar* ParDerivFactory::createInvertPar(const char* type, Options *opt) { opt = Options::getRoot()->getSection(default_section); if(!strcasecmp(type, PARDERIVSERIAL)) { - return new InvertParSerial(opt); + return new InvertParSerial(opt, mesh_in); }else if(!strcasecmp(type, PARDERIVCYCLIC)) { - return new InvertParCR(opt); + return new InvertParCR(opt, mesh_in); } throw BoutException("No such ParDeriv solver exists in this build, type: %s", type); } -InvertPar* ParDerivFactory::createInvertPar(Options *opts) { +InvertPar* ParDerivFactory::createInvertPar(Options *opts, Mesh *mesh_in) { string type; opts->get("type", type, "cyclic"); - return createInvertPar(type.c_str(), opts); + return createInvertPar(type.c_str(), opts, mesh_in); } diff --git a/src/invert/parderiv/parderiv_factory.hxx b/src/invert/parderiv/parderiv_factory.hxx index 9bdc314419..d5a3a5cb17 100644 --- a/src/invert/parderiv/parderiv_factory.hxx +++ b/src/invert/parderiv/parderiv_factory.hxx @@ -11,9 +11,9 @@ class ParDerivFactory { /// Return a pointer to the only instance static ParDerivFactory* getInstance(); - InvertPar* createInvertPar(); - InvertPar *createInvertPar(const char *type, Options *opt = nullptr); - InvertPar* createInvertPar(Options *opts); + InvertPar* createInvertPar(Mesh* mesh_in = mesh); + InvertPar *createInvertPar(const char *type, Options *opt = nullptr, Mesh* mesh_in = mesh); + InvertPar* createInvertPar(Options *opts, Mesh* mesh_in = mesh); private: ParDerivFactory() {} // Prevent instantiation of this class static ParDerivFactory* instance; ///< The only instance of this class (Singleton) From 8c887838723e1337248864c59750e9fea70c4f0b Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 14 Nov 2018 13:43:15 +0000 Subject: [PATCH 4/4] Save Coordinates* pointer in Laplacian Always need the same Coordinates*, so call getCoordinates(location) once in constructor instead of in each call to solve(). --- include/invert_laplace.hxx | 2 + .../laplace/impls/cyclic/cyclic_laplace.cxx | 12 +-- .../impls/multigrid/multigrid_laplace.cxx | 3 - .../laplace/impls/naulin/naulin_laplace.cxx | 1 - src/invert/laplace/impls/pdd/pdd.cxx | 2 +- .../laplace/impls/petsc/petsc_laplace.cxx | 81 +++++++++---------- .../laplace/impls/serial_band/serial_band.cxx | 80 +++++++++--------- .../laplace/impls/serial_tri/serial_tri.cxx | 2 +- .../laplace/impls/shoot/shoot_laplace.cxx | 4 +- src/invert/laplace/impls/spt/spt.cxx | 2 +- src/invert/laplace/invert_laplace.cxx | 78 +++++++++--------- 11 files changed, 123 insertions(+), 144 deletions(-) diff --git a/include/invert_laplace.hxx b/include/invert_laplace.hxx index 005cc61a78..17d1e48507 100644 --- a/include/invert_laplace.hxx +++ b/include/invert_laplace.hxx @@ -229,6 +229,8 @@ protected: bool includeguards=true); CELL_LOC location; ///< staggered grid location of this solver Mesh* localmesh; ///< Mesh object for this solver + Coordinates* coords; ///< Coordinates object, so we only have to call + /// localmesh->getCoordinates(location) once private: /// Singleton instance static Laplacian *instance; diff --git a/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx b/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx index 4d0df4e7c4..ff31e5a77b 100644 --- a/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx +++ b/src/invert/laplace/impls/cyclic/cyclic_laplace.cxx @@ -95,8 +95,6 @@ const FieldPerp LaplaceCyclic::solve(const FieldPerp &rhs, const FieldPerp &x0) FieldPerp x(localmesh); // Result x.allocate(); - Coordinates *coord = localmesh->getCoordinates(location); - int jy = rhs.getIndex(); // Get the Y index x.setIndex(jy); @@ -142,7 +140,7 @@ const FieldPerp LaplaceCyclic::solve(const FieldPerp &rhs, const FieldPerp &x0) // including boundary conditions BOUT_OMP(for nowait) for (int kz = 0; kz < nmode; kz++) { - BoutReal zlen = coord->dz * (localmesh->LocalNz - 3); + BoutReal zlen = coords->dz * (localmesh->LocalNz - 3); BoutReal kwave = kz * 2.0 * PI / (2. * zlen); // wave number is 1/[rad]; DST has extra 2. @@ -210,7 +208,7 @@ const FieldPerp LaplaceCyclic::solve(const FieldPerp &rhs, const FieldPerp &x0) // including boundary conditions BOUT_OMP(for nowait) for (int kz = 0; kz < nmode; kz++) { - BoutReal kwave = kz * 2.0 * PI / (coord->zlength()); // wave number is 1/[rad] + BoutReal kwave = kz * 2.0 * PI / (coords->zlength()); // wave number is 1/[rad] tridagMatrix(&a(kz, 0), &b(kz, 0), &c(kz, 0), &bcmplx(kz, 0), jy, kz, // True for the component constant (DC) in Z kwave, // Z wave number @@ -259,8 +257,6 @@ const Field3D LaplaceCyclic::solve(const Field3D &rhs, const Field3D &x0) { x.allocate(); x.setLocation(location); - Coordinates *coord = rhs.getCoordinates(); - // Get the width of the boundary // If the flags to assign that only one guard cell should be used is set @@ -341,7 +337,7 @@ const Field3D LaplaceCyclic::solve(const Field3D &rhs, const Field3D &x0) { int iy = ys + ind / nmode; int kz = ind % nmode; - BoutReal zlen = coord->dz * (localmesh->LocalNz - 3); + BoutReal zlen = coords->dz * (localmesh->LocalNz - 3); BoutReal kwave = kz * 2.0 * PI / (2. * zlen); // wave number is 1/[rad]; DST has extra 2. @@ -421,7 +417,7 @@ const Field3D LaplaceCyclic::solve(const Field3D &rhs, const Field3D &x0) { int iy = ys + ind / nmode; int kz = ind % nmode; - BoutReal kwave = kz * 2.0 * PI / (coord->zlength()); // wave number is 1/[rad] + BoutReal kwave = kz * 2.0 * PI / (coords->zlength()); // wave number is 1/[rad] tridagMatrix(&a3D(ind, 0), &b3D(ind, 0), &c3D(ind, 0), &bcmplx3D(ind, 0), iy, kz, // True for the component constant (DC) in Z kwave, // Z wave number diff --git a/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx b/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx index 567ec0e583..d3e419e54f 100644 --- a/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx +++ b/src/invert/laplace/impls/multigrid/multigrid_laplace.cxx @@ -209,8 +209,6 @@ const FieldPerp LaplaceMultigrid::solve(const FieldPerp &b_in, const FieldPerp & BoutReal t0,t1; - Coordinates *coords = localmesh->getCoordinates(location); - yindex = b_in.getIndex(); int level = kMG->mglevel-1; int lzz = kMG->lnz[level]; @@ -555,7 +553,6 @@ void LaplaceMultigrid::generateMatrixF(int level) { // Set (fine-level) matrix entries - Coordinates *coords = localmesh->getCoordinates(location); BoutReal *mat; mat = kMG->matmg[level]; int llx = kMG->lnx[level]; diff --git a/src/invert/laplace/impls/naulin/naulin_laplace.cxx b/src/invert/laplace/impls/naulin/naulin_laplace.cxx index 87e1093a61..a216d9de70 100644 --- a/src/invert/laplace/impls/naulin/naulin_laplace.cxx +++ b/src/invert/laplace/impls/naulin/naulin_laplace.cxx @@ -172,7 +172,6 @@ const Field3D LaplaceNaulin::solve(const Field3D &rhs, const Field3D &x0) { ASSERT1(Acoef.getLocation() == location); ASSERT1(localmesh == rhs.getMesh() && localmesh == x0.getMesh()); - Coordinates *coords = rhs.getCoordinates(); Field3D x(x0); // Result Field3D rhsOverD = rhs/Dcoef; diff --git a/src/invert/laplace/impls/pdd/pdd.cxx b/src/invert/laplace/impls/pdd/pdd.cxx index 4a2fbbf46a..23c09145c7 100644 --- a/src/invert/laplace/impls/pdd/pdd.cxx +++ b/src/invert/laplace/impls/pdd/pdd.cxx @@ -164,7 +164,7 @@ void LaplacePDD::start(const FieldPerp &b, PDD_data &data) { /// Create the matrices to be inverted (one for each z point) - BoutReal kwaveFactor = 2.0 * PI / localmesh->getCoordinates(location)->zlength(); + BoutReal kwaveFactor = 2.0 * PI / coords->zlength(); /// Set matrix elements for (int kz = 0; kz <= maxmode; kz++) { diff --git a/src/invert/laplace/impls/petsc/petsc_laplace.cxx b/src/invert/laplace/impls/petsc/petsc_laplace.cxx index 1555345456..d8cf5a0ff6 100644 --- a/src/invert/laplace/impls/petsc/petsc_laplace.cxx +++ b/src/invert/laplace/impls/petsc/petsc_laplace.cxx @@ -362,9 +362,6 @@ const FieldPerp LaplacePetsc::solve(const FieldPerp &b, const FieldPerp &x0) { } #endif - // Get the metric tensor - Coordinates* coord = localmesh->getCoordinates(location); - int y = b.getIndex(); // Get the Y index sol.setIndex(y); // Initialize the solution field. sol = 0.; @@ -402,21 +399,21 @@ const FieldPerp LaplacePetsc::solve(const FieldPerp &b, const FieldPerp &x0) { // Set values corresponding to nodes adjacent in x if( fourth_order ) { // Fourth Order Accuracy on Boundary - Element(i,x,z, 0, 0, -25.0 / (12.0*coord->dx(x,y)) / sqrt(coord->g_11(x,y)), MatA ); - Element(i,x,z, 1, 0, 4.0 / coord->dx(x,y) / sqrt(coord->g_11(x,y)), MatA ); - Element(i,x,z, 2, 0, -3.0 / coord->dx(x,y) / sqrt(coord->g_11(x,y)), MatA ); - Element(i,x,z, 3, 0, 4.0 / (3.0*coord->dx(x,y)) / sqrt(coord->g_11(x,y)), MatA ); - Element(i,x,z, 4, 0, -1.0 / (4.0*coord->dx(x,y)) / sqrt(coord->g_11(x,y)), MatA ); + Element(i,x,z, 0, 0, -25.0 / (12.0*coords->dx(x,y)) / sqrt(coords->g_11(x,y)), MatA ); + Element(i,x,z, 1, 0, 4.0 / coords->dx(x,y) / sqrt(coords->g_11(x,y)), MatA ); + Element(i,x,z, 2, 0, -3.0 / coords->dx(x,y) / sqrt(coords->g_11(x,y)), MatA ); + Element(i,x,z, 3, 0, 4.0 / (3.0*coords->dx(x,y)) / sqrt(coords->g_11(x,y)), MatA ); + Element(i,x,z, 4, 0, -1.0 / (4.0*coords->dx(x,y)) / sqrt(coords->g_11(x,y)), MatA ); } else { // // Second Order Accuracy on Boundary -// Element(i,x,z, 0, 0, -3.0 / (2.0*coord->dx(x,y)), MatA ); -// Element(i,x,z, 1, 0, 2.0 / coord->dx(x,y), MatA ); -// Element(i,x,z, 2, 0, -1.0 / (2.0*coord->dx(x,y)), MatA ); +// Element(i,x,z, 0, 0, -3.0 / (2.0*coords->dx(x,y)), MatA ); +// Element(i,x,z, 1, 0, 2.0 / coords->dx(x,y), MatA ); +// Element(i,x,z, 2, 0, -1.0 / (2.0*coords->dx(x,y)), MatA ); // // Element(i,x,z, 3, 0, 0.0, MatA ); // Reset these elements to 0 in case 4th order flag was used previously: not allowed now // // Element(i,x,z, 4, 0, 0.0, MatA ); // Second Order Accuracy on Boundary, set half-way between grid points - Element(i,x,z, 0, 0, -1.0 / coord->dx(x,y) / sqrt(coord->g_11(x,y)), MatA ); - Element(i,x,z, 1, 0, 1.0 / coord->dx(x,y) / sqrt(coord->g_11(x,y)), MatA ); + Element(i,x,z, 0, 0, -1.0 / coords->dx(x,y) / sqrt(coords->g_11(x,y)), MatA ); + Element(i,x,z, 1, 0, 1.0 / coords->dx(x,y) / sqrt(coords->g_11(x,y)), MatA ); Element(i,x,z, 2, 0, 0.0, MatA ); // Element(i,x,z, 3, 0, 0.0, MatA ); // Reset these elements to 0 in case 4th order flag was used previously: not allowed now // Element(i,x,z, 4, 0, 0.0, MatA ); @@ -473,11 +470,11 @@ const FieldPerp LaplacePetsc::solve(const FieldPerp &b, const FieldPerp &x0) { // Set the matrix coefficients Coeffs( x, y, z, A1, A2, A3, A4, A5 ); - BoutReal dx = coord->dx(x,y); - BoutReal dx2 = SQ(coord->dx(x,y)); - BoutReal dz = coord->dz; - BoutReal dz2 = SQ(coord->dz); - BoutReal dxdz = coord->dx(x,y) * coord->dz; + BoutReal dx = coords->dx(x,y); + BoutReal dx2 = SQ(coords->dx(x,y)); + BoutReal dz = coords->dz; + BoutReal dz2 = SQ(coords->dz); + BoutReal dxdz = coords->dx(x,y) * coords->dz; ASSERT3(finite(A1)); ASSERT3(finite(A2)); @@ -651,22 +648,22 @@ const FieldPerp LaplacePetsc::solve(const FieldPerp &b, const FieldPerp &x0) { // Set values corresponding to nodes adjacent in x if( fourth_order ) { // Fourth Order Accuracy on Boundary - Element(i,x,z, 0, 0, 25.0 / (12.0*coord->dx(x,y)) / sqrt(coord->g_11(x,y)), MatA ); - Element(i,x,z, -1, 0, -4.0 / coord->dx(x,y) / sqrt(coord->g_11(x,y)), MatA ); - Element(i,x,z, -2, 0, 3.0 / coord->dx(x,y) / sqrt(coord->g_11(x,y)), MatA ); - Element(i,x,z, -3, 0, -4.0 / (3.0*coord->dx(x,y)) / sqrt(coord->g_11(x,y)), MatA ); - Element(i,x,z, -4, 0, 1.0 / (4.0*coord->dx(x,y)) / sqrt(coord->g_11(x,y)), MatA ); + Element(i,x,z, 0, 0, 25.0 / (12.0*coords->dx(x,y)) / sqrt(coords->g_11(x,y)), MatA ); + Element(i,x,z, -1, 0, -4.0 / coords->dx(x,y) / sqrt(coords->g_11(x,y)), MatA ); + Element(i,x,z, -2, 0, 3.0 / coords->dx(x,y) / sqrt(coords->g_11(x,y)), MatA ); + Element(i,x,z, -3, 0, -4.0 / (3.0*coords->dx(x,y)) / sqrt(coords->g_11(x,y)), MatA ); + Element(i,x,z, -4, 0, 1.0 / (4.0*coords->dx(x,y)) / sqrt(coords->g_11(x,y)), MatA ); } else { // // Second Order Accuracy on Boundary -// Element(i,x,z, 0, 0, 3.0 / (2.0*coord->dx(x,y)), MatA ); -// Element(i,x,z, -1, 0, -2.0 / coord->dx(x,y), MatA ); -// Element(i,x,z, -2, 0, 1.0 / (2.0*coord->dx(x,y)), MatA ); +// Element(i,x,z, 0, 0, 3.0 / (2.0*coords->dx(x,y)), MatA ); +// Element(i,x,z, -1, 0, -2.0 / coords->dx(x,y), MatA ); +// Element(i,x,z, -2, 0, 1.0 / (2.0*coords->dx(x,y)), MatA ); // // Element(i,x,z, -3, 0, 0.0, MatA ); // Reset these elements to 0 in case 4th order flag was used previously: not allowed now // // Element(i,x,z, -4, 0, 0.0, MatA ); // Second Order Accuracy on Boundary, set half-way between grid points - Element(i,x,z, 0, 0, 1.0 / coord->dx(x,y) / sqrt(coord->g_11(x,y)), MatA ); - Element(i,x,z, -1, 0, -1.0 / coord->dx(x,y) / sqrt(coord->g_11(x,y)), MatA ); + Element(i,x,z, 0, 0, 1.0 / coords->dx(x,y) / sqrt(coords->g_11(x,y)), MatA ); + Element(i,x,z, -1, 0, -1.0 / coords->dx(x,y) / sqrt(coords->g_11(x,y)), MatA ); Element(i,x,z, -2, 0, 0.0, MatA ); // Element(i,x,z, -3, 0, 0.0, MatA ); // Reset these elements to 0 in case 4th order flag was used previously: not allowed now // Element(i,x,z, -4, 0, 0.0, MatA ); @@ -937,18 +934,16 @@ void LaplacePetsc::Element(int i, int x, int z, */ void LaplacePetsc::Coeffs( int x, int y, int z, BoutReal &coef1, BoutReal &coef2, BoutReal &coef3, BoutReal &coef4, BoutReal &coef5 ) { - Coordinates *coord = localmesh->getCoordinates(location); // Get metric tensor - - coef1 = coord->g11(x,y); // X 2nd derivative coefficient - coef2 = coord->g33(x,y); // Z 2nd derivative coefficient - coef3 = 2.*coord->g13(x,y); // X-Z mixed derivative coefficient + coef1 = coords->g11(x,y); // X 2nd derivative coefficient + coef2 = coords->g33(x,y); // Z 2nd derivative coefficient + coef3 = 2.*coords->g13(x,y); // X-Z mixed derivative coefficient coef4 = 0.0; coef5 = 0.0; // If global flag all_terms are set (true by default) if (all_terms) { - coef4 = coord->G1(x,y); // X 1st derivative - coef5 = coord->G3(x,y); // Z 1st derivative + coef4 = coords->G1(x,y); // X 1st derivative + coef5 = coords->G3(x,y); // Z 1st derivative ASSERT3(finite(coef4)); ASSERT3(finite(coef5)); @@ -957,13 +952,13 @@ void LaplacePetsc::Coeffs( int x, int y, int z, BoutReal &coef1, BoutReal &coef2 if(nonuniform) { // non-uniform mesh correction if((x != 0) && (x != (localmesh->LocalNx-1))) { - coef4 -= 0.5 * ( ( coord->dx(x+1,y) - coord->dx(x-1,y) ) / SQ(coord->dx(x,y)) ) * coef1; // BOUT-06 term + coef4 -= 0.5 * ( ( coords->dx(x+1,y) - coords->dx(x-1,y) ) / SQ(coords->dx(x,y)) ) * coef1; // BOUT-06 term } } if(localmesh->IncIntShear) { // d2dz2 term - coef2 += coord->g11(x,y) * coord->IntShiftTorsion(x,y) * coord->IntShiftTorsion(x,y); + coef2 += coords->g11(x,y) * coords->IntShiftTorsion(x,y) * coords->IntShiftTorsion(x,y); // Mixed derivative coef3 = 0.0; // This cancels out } @@ -993,19 +988,19 @@ void LaplacePetsc::Coeffs( int x, int y, int z, BoutReal &coef1, BoutReal &coef2 int zmm = z-2; // z minus 1 minus 1 if (zmm<0) zmm += meshz; // Fourth order discretization of C in x - ddx_C = (-C2(x+2,y,z) + 8.*C2(x+1,y,z) - 8.*C2(x-1,y,z) + C2(x-2,y,z)) / (12.*coord->dx(x,y)*(C1(x,y,z))); + ddx_C = (-C2(x+2,y,z) + 8.*C2(x+1,y,z) - 8.*C2(x-1,y,z) + C2(x-2,y,z)) / (12.*coords->dx(x,y)*(C1(x,y,z))); // Fourth order discretization of C in z - ddz_C = (-C2(x,y,zpp) + 8.*C2(x,y,zp) - 8.*C2(x,y,zm) + C2(x,y,zmm)) / (12.*coord->dz*(C1(x,y,z))); + ddz_C = (-C2(x,y,zpp) + 8.*C2(x,y,zp) - 8.*C2(x,y,zm) + C2(x,y,zmm)) / (12.*coords->dz*(C1(x,y,z))); } else { // Second order discretization of C in x - ddx_C = (C2(x+1,y,z) - C2(x-1,y,z)) / (2.*coord->dx(x,y)*(C1(x,y,z))); + ddx_C = (C2(x+1,y,z) - C2(x-1,y,z)) / (2.*coords->dx(x,y)*(C1(x,y,z))); // Second order discretization of C in z - ddz_C = (C2(x,y,zp) - C2(x,y,zm)) / (2.*coord->dz*(C1(x,y,z))); + ddz_C = (C2(x,y,zp) - C2(x,y,zm)) / (2.*coords->dz*(C1(x,y,z))); } - coef4 += coord->g11(x,y) * ddx_C + coord->g13(x,y) * ddz_C; - coef5 += coord->g13(x,y) * ddx_C + coord->g33(x,y) * ddz_C; + coef4 += coords->g11(x,y) * ddx_C + coords->g13(x,y) * ddz_C; + coef5 += coords->g13(x,y) * ddx_C + coords->g33(x,y) * ddz_C; } } diff --git a/src/invert/laplace/impls/serial_band/serial_band.cxx b/src/invert/laplace/impls/serial_band/serial_band.cxx index b9541bf867..46835d8174 100644 --- a/src/invert/laplace/impls/serial_band/serial_band.cxx +++ b/src/invert/laplace/impls/serial_band/serial_band.cxx @@ -88,8 +88,6 @@ const FieldPerp LaplaceSerialBand::solve(const FieldPerp &b, const FieldPerp &x0 int jy = b.getIndex(); x.setIndex(jy); - Coordinates *coord = localmesh->getCoordinates(location); - int ncz = localmesh->LocalNz; int ncx = localmesh->LocalNx-1; @@ -128,7 +126,7 @@ const FieldPerp LaplaceSerialBand::solve(const FieldPerp &b, const FieldPerp &x0 ///////// PERFORM INVERSION ///////// // shift freqs according to FFT convention - kwave=iz*2.0*PI/coord->zlength(); // wave number is 1/[rad] + kwave=iz*2.0*PI/coords->zlength(); // wave number is 1/[rad] // set bk1d for(int ix=0;ixLocalNx;ix++) @@ -150,9 +148,9 @@ const FieldPerp LaplaceSerialBand::solve(const FieldPerp &b, const FieldPerp &x0 A(ix, 4) = 0.; #else // Set coefficients - coef1 = coord->g11(ix,jy); // X 2nd derivative - coef2 = coord->g33(ix,jy); // Z 2nd derivative - coef3 = coord->g13(ix,jy); // X-Z mixed derivatives + coef1 = coords->g11(ix,jy); // X 2nd derivative + coef2 = coords->g33(ix,jy); // Z 2nd derivative + coef3 = coords->g13(ix,jy); // X-Z mixed derivatives coef4 = 0.0; // X 1st derivative coef5 = 0.0; // Z 1st derivative coef6 = Acoef(ix,jy); // Constant @@ -163,26 +161,26 @@ const FieldPerp LaplaceSerialBand::solve(const FieldPerp &b, const FieldPerp &x0 coef3 *= Dcoef(ix,jy); if(all_terms) { - coef4 = coord->G1(ix,jy); - coef5 = coord->G3(ix,jy); + coef4 = coords->G1(ix,jy); + coef5 = coords->G3(ix,jy); } if(nonuniform) { // non-uniform localmesh correction if((ix != 0) && (ix != ncx)) - coef4 += coord->g11(ix,jy)*( (1.0/coord->dx(ix+1,jy)) - (1.0/coord->dx(ix-1,jy)) )/(2.0*coord->dx(ix,jy)); + coef4 += coords->g11(ix,jy)*( (1.0/coords->dx(ix+1,jy)) - (1.0/coords->dx(ix-1,jy)) )/(2.0*coords->dx(ix,jy)); } // A first order derivative term (1/c)\nabla_perp c\cdot\nabla_\perp x if((ix > 1) && (ix < (localmesh->LocalNx-2))) - coef4 += coord->g11(ix,jy) * (Ccoef(ix-2,jy) - 8.*Ccoef(ix-1,jy) + 8.*Ccoef(ix+1,jy) - Ccoef(ix+2,jy)) / (12.*coord->dx(ix,jy)*(Ccoef(ix,jy))); + coef4 += coords->g11(ix,jy) * (Ccoef(ix-2,jy) - 8.*Ccoef(ix-1,jy) + 8.*Ccoef(ix+1,jy) - Ccoef(ix+2,jy)) / (12.*coords->dx(ix,jy)*(Ccoef(ix,jy))); // Put into matrix - coef1 /= 12.* SQ(coord->dx(ix,jy)); + coef1 /= 12.* SQ(coords->dx(ix,jy)); coef2 *= SQ(kwave); - coef3 *= kwave / (12. * coord->dx(ix,jy)); - coef4 /= 12. * coord->dx(ix,jy); + coef3 *= kwave / (12. * coords->dx(ix,jy)); + coef4 /= 12. * coords->dx(ix,jy); coef5 *= kwave; A(ix, 0) = dcomplex(-coef1 + coef4, coef3); @@ -198,9 +196,9 @@ const FieldPerp LaplaceSerialBand::solve(const FieldPerp &b, const FieldPerp &x0 int ix = 1; - coef1=coord->g11(ix,jy)/(SQ(coord->dx(ix,jy))); - coef2=coord->g33(ix,jy); - coef3= kwave * coord->g13(ix,jy)/(2. * coord->dx(ix,jy)); + coef1=coords->g11(ix,jy)/(SQ(coords->dx(ix,jy))); + coef2=coords->g33(ix,jy); + coef3= kwave * coords->g13(ix,jy)/(2. * coords->dx(ix,jy)); // Multiply Delp2 component by a factor coef1 *= Dcoef(ix,jy); @@ -215,9 +213,9 @@ const FieldPerp LaplaceSerialBand::solve(const FieldPerp &b, const FieldPerp &x0 ix = ncx-1; - coef1=coord->g11(ix,jy)/(SQ(coord->dx(ix,jy))); - coef2=coord->g33(ix,jy); - coef3= kwave * coord->g13(ix,jy)/(2. * coord->dx(ix,jy)); + coef1=coords->g11(ix,jy)/(SQ(coords->dx(ix,jy))); + coef2=coords->g33(ix,jy); + coef3= kwave * coords->g13(ix,jy)/(2. * coords->dx(ix,jy)); A(ix, 0) = 0.0; A(ix, 1) = dcomplex(coef1, -coef3); @@ -253,8 +251,8 @@ const FieldPerp LaplaceSerialBand::solve(const FieldPerp &b, const FieldPerp &x0 for (int ix=0;ixg_11(ix, jy)) / coord->dx(ix, jy); - A(ix, 3) = .5 / sqrt(coord->g_11(ix, jy)) / coord->dx(ix, jy); + A(ix, 2) = -.5 / sqrt(coords->g_11(ix, jy)) / coords->dx(ix, jy); + A(ix, 3) = .5 / sqrt(coords->g_11(ix, jy)) / coords->dx(ix, jy); A(ix, 4) = 0.; } @@ -275,18 +273,18 @@ const FieldPerp LaplaceSerialBand::solve(const FieldPerp &b, const FieldPerp &x0 for (int ix=0;ixg_22(ix, jy)); - A(ix, 3) = 4. / sqrt(coord->g_22(ix + 1, jy)); - A(ix, 4) = -1. / sqrt(coord->g_22(ix + 2, jy)); + A(ix, 2) = -3. / sqrt(coords->g_22(ix, jy)); + A(ix, 3) = 4. / sqrt(coords->g_22(ix + 1, jy)); + A(ix, 4) = -1. / sqrt(coords->g_22(ix + 2, jy)); } } else if(inner_boundary_flags & INVERT_DC_GRADPARINV) { for (int ix=0;ixg_22(ix, jy)); - A(ix, 3) = 4. * sqrt(coord->g_22(ix + 1, jy)); - A(ix, 4) = -sqrt(coord->g_22(ix + 2, jy)); + A(ix, 2) = -3. * sqrt(coords->g_22(ix, jy)); + A(ix, 3) = 4. * sqrt(coords->g_22(ix + 1, jy)); + A(ix, 4) = -sqrt(coords->g_22(ix + 2, jy)); } } else if (inner_boundary_flags & INVERT_DC_LAP) { @@ -319,27 +317,27 @@ const FieldPerp LaplaceSerialBand::solve(const FieldPerp &b, const FieldPerp &x0 int ix = 1; - coef1=coord->g11(ix,jy)/(12.* SQ(coord->dx(ix,jy))); + coef1=coords->g11(ix,jy)/(12.* SQ(coords->dx(ix,jy))); - coef2=coord->g33(ix,jy); + coef2=coords->g33(ix,jy); - coef3= kwave * coord->g13(ix,jy)/(2. * coord->dx(ix,jy)); + coef3= kwave * coords->g13(ix,jy)/(2. * coords->dx(ix,jy)); coef4 = Acoef(ix,jy); // Combine 4th order at 1 with 2nd order at 0 A(1, 0) = 0.0; // Not used A(1, 1) = dcomplex( - (14. - SQ(coord->dx(0, jy) * kwave) * coord->g33(0, jy) / coord->g11(0, jy)) * + (14. - SQ(coords->dx(0, jy) * kwave) * coords->g33(0, jy) / coords->g11(0, jy)) * coef1, -coef3); A(1, 2) = dcomplex(-29. * coef1 - SQ(kwave) * coef2 + coef4, 0.0); A(1, 3) = dcomplex(16. * coef1, coef3); A(1, 4) = dcomplex(-coef1, 0.0); - coef1=coord->g11(ix,jy)/(SQ(coord->dx(ix,jy))); - coef2=coord->g33(ix,jy); - coef3= kwave * coord->g13(ix,jy)/(2. * coord->dx(ix,jy)); + coef1=coords->g11(ix,jy)/(SQ(coords->dx(ix,jy))); + coef2=coords->g33(ix,jy); + coef3= kwave * coords->g13(ix,jy)/(2. * coords->dx(ix,jy)); // Use 2nd order at 1 A(0, 0) = 0.0; // Should never be used @@ -361,11 +359,11 @@ const FieldPerp LaplaceSerialBand::solve(const FieldPerp &b, const FieldPerp &x0 int ix = ncx-1; - coef1=coord->g11(ix,jy)/(12.* SQ(coord->dx(ix,jy))); + coef1=coords->g11(ix,jy)/(12.* SQ(coords->dx(ix,jy))); - coef2=coord->g33(ix,jy); + coef2=coords->g33(ix,jy); - coef3= kwave * coord->g13(ix,jy)/(2. * coord->dx(ix,jy)); + coef3= kwave * coords->g13(ix,jy)/(2. * coords->dx(ix,jy)); coef4 = Acoef(ix,jy); @@ -375,14 +373,14 @@ const FieldPerp LaplaceSerialBand::solve(const FieldPerp &b, const FieldPerp &x0 A(ix, 2) = dcomplex(-29. * coef1 - SQ(kwave) * coef2 + coef4, 0.0); A(ix, 3) = dcomplex( (14. - - SQ(coord->dx(ncx, jy) * kwave) * coord->g33(ncx, jy) / coord->g11(ncx, jy)) * + SQ(coords->dx(ncx, jy) * kwave) * coords->g33(ncx, jy) / coords->g11(ncx, jy)) * coef1, coef3); A(ix, 4) = 0.0; // Not used - coef1=coord->g11(ix,jy)/(SQ(coord->dx(ix,jy))); - coef2=coord->g33(ix,jy); - coef3= kwave * coord->g13(ix,jy)/(2. * coord->dx(ix,jy)); + coef1=coords->g11(ix,jy)/(SQ(coords->dx(ix,jy))); + coef2=coords->g33(ix,jy); + coef3= kwave * coords->g13(ix,jy)/(2. * coords->dx(ix,jy)); // Use 2nd order at ncx - 1 A(ncx, 0) = dcomplex(coef1, -coef3); diff --git a/src/invert/laplace/impls/serial_tri/serial_tri.cxx b/src/invert/laplace/impls/serial_tri/serial_tri.cxx index 0a2cbe02b1..fc4c226ac7 100644 --- a/src/invert/laplace/impls/serial_tri/serial_tri.cxx +++ b/src/invert/laplace/impls/serial_tri/serial_tri.cxx @@ -84,7 +84,7 @@ const FieldPerp LaplaceSerialTri::solve(const FieldPerp &b, const FieldPerp &x0) int ncz = localmesh->LocalNz; // No of z pnts int ncx = localmesh->LocalNx; // No of x pnts - BoutReal kwaveFactor = 2.0 * PI / localmesh->getCoordinates(location)->zlength(); + BoutReal kwaveFactor = 2.0 * PI / coords->zlength(); // Setting the width of the boundary. // NOTE: The default is a width of 2 guard cells diff --git a/src/invert/laplace/impls/shoot/shoot_laplace.cxx b/src/invert/laplace/impls/shoot/shoot_laplace.cxx index 86430f928c..3b515e6ed7 100644 --- a/src/invert/laplace/impls/shoot/shoot_laplace.cxx +++ b/src/invert/laplace/impls/shoot/shoot_laplace.cxx @@ -77,8 +77,6 @@ const FieldPerp LaplaceShoot::solve(const FieldPerp &rhs) { int jy = rhs.getIndex(); // Get the Y index x.setIndex(jy); - Coordinates *coord = localmesh->getCoordinates(location); - // Get the width of the boundary int inbndry = localmesh->xstart, outbndry=localmesh->xstart; @@ -135,7 +133,7 @@ const FieldPerp LaplaceShoot::solve(const FieldPerp &rhs) { rfft(rhs[ix], localmesh->LocalNz, std::begin(rhsk)); for(int kz=0; kzzlength()); // wave number is 1/[rad] + BoutReal kwave=kz*2.0*PI/(coords->zlength()); // wave number is 1/[rad] // Get the coefficients dcomplex a,b,c; diff --git a/src/invert/laplace/impls/spt/spt.cxx b/src/invert/laplace/impls/spt/spt.cxx index e75cad6cae..42b1f30935 100644 --- a/src/invert/laplace/impls/spt/spt.cxx +++ b/src/invert/laplace/impls/spt/spt.cxx @@ -291,7 +291,7 @@ int LaplaceSPT::start(const FieldPerp &b, SPT_data &data) { data.bk(kz, ix) = dc1d[kz]; } - BoutReal kwaveFactor = 2.0 * PI / localmesh->getCoordinates(location)->zlength(); + BoutReal kwaveFactor = 2.0 * PI / coords->zlength(); /// Set matrix elements for (int kz = 0; kz <= maxmode; kz++) { diff --git a/src/invert/laplace/invert_laplace.cxx b/src/invert/laplace/invert_laplace.cxx index 2e66fa3767..e120f86244 100644 --- a/src/invert/laplace/invert_laplace.cxx +++ b/src/invert/laplace/invert_laplace.cxx @@ -65,6 +65,8 @@ Laplacian::Laplacian(Options *options, const CELL_LOC loc, Mesh *mesh_in) location = CELL_CENTRE; } + coords = localmesh->getCoordinates(location); + // Communication option. Controls if asyncronous sends are used options->get("async", async_send, true); @@ -81,7 +83,7 @@ Laplacian::Laplacian(Options *options, const CELL_LOC loc, Mesh *mesh_in) OPTION(options, low_mem, false); OPTION(options, nonuniform, - localmesh->getCoordinates(location)->non_uniform); // Default is the mesh setting + coords->non_uniform); // Default is the mesh setting OPTION(options, all_terms, true); // Include first derivative terms @@ -249,9 +251,7 @@ void Laplacian::tridagCoefs(int jx, int jy, int jz, ASSERT1(ccoef == nullptr || ccoef->getLocation() == loc); ASSERT1(d == nullptr || d->getLocation() == loc); - Coordinates *coord = localmesh->getCoordinates(loc); - - BoutReal kwave=jz*2.0*PI/coord->zlength(); // wave number is 1/[rad] + BoutReal kwave=jz*2.0*PI/coords->zlength(); // wave number is 1/[rad] tridagCoefs(jx, jy, kwave, a, b, c, @@ -296,18 +296,16 @@ void Laplacian::tridagCoefs(int jx, int jy, BoutReal kwave, BoutReal coef1, coef2, coef3, coef4, coef5; - Coordinates *coord = localmesh->getCoordinates(loc); - - coef1=coord->g11(jx,jy); ///< X 2nd derivative coefficient - coef2=coord->g33(jx,jy); ///< Z 2nd derivative coefficient - coef3=2.*coord->g13(jx,jy); ///< X-Z mixed derivative coefficient + coef1=coords->g11(jx,jy); ///< X 2nd derivative coefficient + coef2=coords->g33(jx,jy); ///< Z 2nd derivative coefficient + coef3=2.*coords->g13(jx,jy); ///< X-Z mixed derivative coefficient coef4 = 0.0; coef5 = 0.0; // If global flag all_terms are set (true by default) if(all_terms) { - coef4 = coord->G1(jx,jy); // X 1st derivative - coef5 = coord->G3(jx,jy); // Z 1st derivative + coef4 = coords->G1(jx,jy); // X 1st derivative + coef5 = coords->G3(jx,jy); // Z 1st derivative } if (d != nullptr) { @@ -322,26 +320,26 @@ void Laplacian::tridagCoefs(int jx, int jy, BoutReal kwave, if(nonuniform) { // non-uniform mesh correction if((jx != 0) && (jx != (localmesh->LocalNx-1))) { - coef4 -= 0.5*((coord->dx(jx+1,jy) - coord->dx(jx-1,jy))/SQ(coord->dx(jx,jy)))*coef1; + coef4 -= 0.5*((coords->dx(jx+1,jy) - coords->dx(jx-1,jy))/SQ(coords->dx(jx,jy)))*coef1; } } if (ccoef != nullptr) { // A first order derivative term if((jx > 0) && (jx < (localmesh->LocalNx-1))) - coef4 += coord->g11(jx,jy) * ((*ccoef)(jx+1,jy) - (*ccoef)(jx-1,jy)) / (2.*coord->dx(jx,jy)*((*ccoef)(jx,jy))); + coef4 += coords->g11(jx,jy) * ((*ccoef)(jx+1,jy) - (*ccoef)(jx-1,jy)) / (2.*coords->dx(jx,jy)*((*ccoef)(jx,jy))); } if(localmesh->IncIntShear) { // d2dz2 term - coef2 += coord->g11(jx,jy) * coord->IntShiftTorsion(jx,jy) * coord->IntShiftTorsion(jx,jy); + coef2 += coords->g11(jx,jy) * coords->IntShiftTorsion(jx,jy) * coords->IntShiftTorsion(jx,jy); // Mixed derivative coef3 = 0.0; // This cancels out } - coef1 /= SQ(coord->dx(jx,jy)); - coef3 /= 2.*coord->dx(jx,jy); - coef4 /= 2.*coord->dx(jx,jy); + coef1 /= SQ(coords->dx(jx,jy)); + coef3 /= 2.*coords->dx(jx,jy); + coef4 /= 2.*coords->dx(jx,jy); a = dcomplex(coef1 - coef4,-kwave*coef3); b = dcomplex(-2.0*coef1 - SQ(kwave)*coef2,kwave*coef5); @@ -361,11 +359,9 @@ void Laplacian::tridagMatrix(dcomplex **avec, dcomplex **bvec, dcomplex **cvec, ASSERT1(ccoef->getLocation() == location); ASSERT1(d->getLocation() == location); - Coordinates *coord = localmesh->getCoordinates(location); - BOUT_OMP(parallel for) for(int kz = 0; kz <= maxmode; kz++) { - BoutReal kwave=kz*2.0*PI/coord->zlength(); // wave number is 1/[rad] + BoutReal kwave=kz*2.0*PI/coords->zlength(); // wave number is 1/[rad] tridagMatrix(avec[kz], bvec[kz], cvec[kz], bk[kz], @@ -426,8 +422,6 @@ void Laplacian::tridagMatrix(dcomplex *avec, dcomplex *bvec, dcomplex *cvec, int xs = 0; // xstart set to the start of x on this processor (including ghost points) int xe = localmesh->LocalNx-1; // xend set to the end of x on this processor (including ghost points) - Coordinates *coord = localmesh->getCoordinates(location); - // Do not want boundary cells if x is periodic for cyclic solver. Only other solver which // works with periodicX is serial_tri, which uses includeguards==true, so the below isn't called. if(!includeguards) { @@ -481,8 +475,8 @@ void Laplacian::tridagMatrix(dcomplex *avec, dcomplex *bvec, dcomplex *cvec, // Zero gradient at inner boundary for (int ix=0;ixg_11(ix,jy))/coord->dx(ix,jy); - cvec[ix] = 1./sqrt(coord->g_11(ix,jy))/coord->dx(ix,jy); + bvec[ix] = -1./sqrt(coords->g_11(ix,jy))/coords->dx(ix,jy); + cvec[ix] = 1./sqrt(coords->g_11(ix,jy))/coords->dx(ix,jy); } } else if(inner_boundary_flags & INVERT_DC_GRAD) { @@ -496,15 +490,15 @@ void Laplacian::tridagMatrix(dcomplex *avec, dcomplex *bvec, dcomplex *cvec, else if(inner_boundary_flags & INVERT_DC_GRADPAR) { for (int ix=0;ixg_22(ix,jy)); - cvec[ix] = -1.0/sqrt(coord->g_22(ix+1,jy)); + bvec[ix] = 1.0/sqrt(coords->g_22(ix,jy)); + cvec[ix] = -1.0/sqrt(coords->g_22(ix+1,jy)); } } else if(inner_boundary_flags & INVERT_DC_GRADPARINV) { for (int ix=0;ixg_22(ix,jy)); - cvec[ix] = -sqrt(coord->g_22(ix+1,jy)); + bvec[ix] = sqrt(coords->g_22(ix,jy)); + cvec[ix] = -sqrt(coords->g_22(ix+1,jy)); } } else if (inner_boundary_flags & INVERT_DC_LAP) { @@ -519,7 +513,7 @@ void Laplacian::tridagMatrix(dcomplex *avec, dcomplex *bvec, dcomplex *cvec, for (int ix=0;ixdx(ix,jy)/sqrt(coord->g11(ix,jy))); + cvec[ix] = -exp(-k*coords->dx(ix,jy)/sqrt(coords->g11(ix,jy))); } } else if (inner_boundary_flags & INVERT_IN_CYLINDER){ @@ -581,8 +575,8 @@ void Laplacian::tridagMatrix(dcomplex *avec, dcomplex *bvec, dcomplex *cvec, // Zero gradient at inner boundary for (int ix=0;ixg_11(ix,jy))/coord->dx(ix,jy); - cvec[ix] = dcomplex(1.,0.)/sqrt(coord->g_11(ix,jy))/coord->dx(ix,jy); + bvec[ix] = dcomplex(-1.,0.)/sqrt(coords->g_11(ix,jy))/coords->dx(ix,jy); + cvec[ix] = dcomplex(1.,0.)/sqrt(coords->g_11(ix,jy))/coords->dx(ix,jy); } } else if(inner_boundary_flags & INVERT_AC_GRAD) { @@ -598,7 +592,7 @@ void Laplacian::tridagMatrix(dcomplex *avec, dcomplex *bvec, dcomplex *cvec, for (int ix=0;ixg33(ix,jy)/coord->g11(ix,jy))*kwave*coord->dx(ix,jy)); + cvec[ix] = -exp(-1.0*sqrt(coords->g33(ix,jy)/coords->g11(ix,jy))*kwave*coords->dx(ix,jy)); } } else if (inner_boundary_flags & INVERT_IN_CYLINDER) { @@ -644,8 +638,8 @@ void Laplacian::tridagMatrix(dcomplex *avec, dcomplex *bvec, dcomplex *cvec, if(outer_boundary_flags & INVERT_DC_GRAD && ( outer_boundary_flags & INVERT_SET || outer_boundary_flags & INVERT_RHS)) { // Zero gradient at outer boundary for (int ix=0;ixg_11(ncx-ix,jy))/coord->dx(ncx-ix,jy); - bvec[ncx-ix]=dcomplex(1.,0.)/sqrt(coord->g_11(ncx-ix,jy))/coord->dx(ncx-ix,jy); + avec[ncx-ix]=dcomplex(-1.,0.)/sqrt(coords->g_11(ncx-ix,jy))/coords->dx(ncx-ix,jy); + bvec[ncx-ix]=dcomplex(1.,0.)/sqrt(coords->g_11(ncx-ix,jy))/coords->dx(ncx-ix,jy); cvec[ncx-ix]=dcomplex(0.,0.); } } @@ -659,15 +653,15 @@ void Laplacian::tridagMatrix(dcomplex *avec, dcomplex *bvec, dcomplex *cvec, } else if(inner_boundary_flags & INVERT_DC_GRADPAR) { for (int ix=0;ixg_22(ncx-ix+1,jy)); - bvec[ncx-ix] = -1.0/sqrt(coord->g_22(ncx-ix,jy)); + avec[ncx-ix] = 1.0/sqrt(coords->g_22(ncx-ix+1,jy)); + bvec[ncx-ix] = -1.0/sqrt(coords->g_22(ncx-ix,jy)); cvec[ncx-ix] = 0.0; } } else if(inner_boundary_flags & INVERT_DC_GRADPARINV) { for (int ix=0;ixg_22(ncx-ix-1,jy)); - bvec[ncx-ix] = -sqrt(coord->g_22(ncx-ix,jy)); + avec[ncx-ix] = sqrt(coords->g_22(ncx-ix-1,jy)); + bvec[ncx-ix] = -sqrt(coords->g_22(ncx-ix,jy)); cvec[ncx-ix] = 0.0; } } @@ -683,7 +677,7 @@ void Laplacian::tridagMatrix(dcomplex *avec, dcomplex *bvec, dcomplex *cvec, for (int ix=0;ixdx(ncx-ix,jy)/sqrt(coord->g11(ncx-ix,jy))); + avec[ncx-ix] = -exp(-k*coords->dx(ncx-ix,jy)/sqrt(coords->g11(ncx-ix,jy))); } } else { @@ -702,8 +696,8 @@ void Laplacian::tridagMatrix(dcomplex *avec, dcomplex *bvec, dcomplex *cvec, if(outer_boundary_flags & INVERT_AC_GRAD && ( outer_boundary_flags & INVERT_SET || outer_boundary_flags & INVERT_RHS)) { // Zero gradient at outer boundary for (int ix=0;ixg_11(ncx-ix,jy))/coord->dx(ncx-ix,jy); - bvec[ncx-ix]=dcomplex(1.,0.)/sqrt(coord->g_11(ncx-ix,jy))/coord->dx(ncx-ix,jy); + avec[ncx-ix]=dcomplex(-1.,0.)/sqrt(coords->g_11(ncx-ix,jy))/coords->dx(ncx-ix,jy); + bvec[ncx-ix]=dcomplex(1.,0.)/sqrt(coords->g_11(ncx-ix,jy))/coords->dx(ncx-ix,jy); cvec[ncx-ix]=dcomplex(0.,0.); } } @@ -718,7 +712,7 @@ void Laplacian::tridagMatrix(dcomplex *avec, dcomplex *bvec, dcomplex *cvec, else if(outer_boundary_flags & INVERT_AC_LAP) { // Use decaying zero-Laplacian solution in the boundary for (int ix=0;ixg33(xe-ix,jy)/coord->g11(xe-ix,jy))*kwave*coord->dx(xe-ix,jy)); + avec[ncx-ix] = -exp(-1.0*sqrt(coords->g33(xe-ix,jy)/coords->g11(xe-ix,jy))*kwave*coords->dx(xe-ix,jy)); bvec[ncx-ix] = 1.0; cvec[ncx-ix] = 0.0; }