Merge branch '542-saamg-wrong-results-with-elliptic-geometry' into 'master'

Resolve "SAAMG: Wrong results with elliptic geometry"

Closes #542

See merge request !371

src/Classic/Algorithms/PartBunch.cpp

@@ -123,8 +123,7 @@ void PartBunch::computeSelfFields(int binNumber) {
 
     if(fs_m->hasValidSolver()) {
         /// Mesh the whole domain
-        if(fs_m->getFieldSolverType() == "SAAMG")
-            resizeMesh();
+        resizeMesh();
 
         /// Scatter charge onto space charge grid.
         this->Q *= this->dt;
@@ -328,12 +327,14 @@ void PartBunch::computeSelfFields(int binNumber) {
 }
 
 void PartBunch::resizeMesh() {
+    if (fs_m->getFieldSolverType() != "SAAMG") {
+        return;
+    }
+
     double xmin = fs_m->solver_m->getXRangeMin();
     double xmax = fs_m->solver_m->getXRangeMax();
     double ymin = fs_m->solver_m->getYRangeMin();
     double ymax = fs_m->solver_m->getYRangeMax();
-    double zmin = fs_m->solver_m->getZRangeMin();
-    double zmax = fs_m->solver_m->getZRangeMax();
 
     if(xmin > rmin_m[0] || xmax < rmax_m[0] ||
        ymin > rmin_m[1] || ymax < rmax_m[1]) {
@@ -354,14 +355,14 @@ void PartBunch::resizeMesh() {
         boundp();
         get_bounds(rmin_m, rmax_m);
     }
-    Vector_t mymin = Vector_t(xmin, ymin , zmin);
-    Vector_t mymax = Vector_t(xmax, ymax , zmax);
 
-    for (int i = 0; i < 3; i++)
-        hr_m[i]   = (mymax[i] - mymin[i])/nr_m[i];
+    Vector_t origin = Vector_t(0.0, 0.0, 0.0);
+
+    // update the mesh origin and mesh spacing hr_m
+    fs_m->solver_m->resizeMesh(origin, hr_m, rmin_m, rmax_m, dh_m);
 
     getMesh().set_meshSpacing(&(hr_m[0]));
-    getMesh().set_origin(mymin);
+    getMesh().set_origin(origin);
 
     rho_m.initialize(getMesh(),
                      getFieldLayout(),
@@ -384,8 +385,7 @@ void PartBunch::computeSelfFields() {
 
     if(fs_m->hasValidSolver()) {
         //mesh the whole domain
-        if(fs_m->getFieldSolverType() == "SAAMG")
-            resizeMesh();
+        resizeMesh();
 
         //scatter charges onto grid
         this->Q *= this->dt;
@@ -510,8 +510,7 @@ void PartBunch::computeSelfFields_cycl(double gamma) {
 
     if(fs_m->hasValidSolver()) {
         /// mesh the whole domain
-        if(fs_m->getFieldSolverType() == "SAAMG")
-            resizeMesh();
+        resizeMesh();
 
         /// scatter particles charge onto grid.
         this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());
@@ -639,8 +638,7 @@ void PartBunch::computeSelfFields_cycl(int bin) {
 
     if(fs_m->hasValidSolver()) {
         /// mesh the whole domain
-        if(fs_m->getFieldSolverType() == "SAAMG")
-            resizeMesh();
+        resizeMesh();
 
         /// scatter particles charge onto grid.
         this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());

src/Classic/Algorithms/PartBunchBase.h

@@ -411,6 +411,8 @@ public:
 //     virtual void setFieldLayout(FieldLayout_t* fLayout) = 0;
     virtual FieldLayout_t &getFieldLayout() = 0;
 
+    virtual void resizeMesh() { };
+
     /*
      * Wrapped member functions of IpplParticleBase
      */

src/Solvers/EllipticDomain.h

 //
 // Class EllipticDomain
-//   :FIXME: add brief description
+//   This class provides an elliptic beam pipe. The mesh adapts to the bunch size
+//   in the longitudinal direction. At the intersection of the mesh with the
+//   beam pipe, three stencil interpolation methods are available.
 //
 // Copyright (c) 2008,        Yves Ineichen, ETH Zürich,
 //               2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland
@@ -32,100 +34,127 @@
 #include <cmath>
 #include "IrregularDomain.h"
 #include "Structure/BoundaryGeometry.h"
+#include "Utilities/OpalException.h"
 
 class EllipticDomain : public IrregularDomain {
 
 public:
 
     EllipticDomain(Vector_t nr, Vector_t hr);
-    EllipticDomain(double semimajor, double semiminor, Vector_t nr, Vector_t hr, std::string interpl);
-    EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr, std::string interpl);
+
+    EllipticDomain(double semimajor, double semiminor, Vector_t nr,
+                   Vector_t hr, std::string interpl);
+
+    EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr,
+                   Vector_t hr, std::string interpl);
+
     EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, std::string interpl);
 
     ~EllipticDomain();
 
     /// returns discretization at (x,y,z)
-    void getBoundaryStencil(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
+    void getBoundaryStencil(int x, int y, int z,
+                            double &W, double &E, double &S,
+                            double &N, double &F, double &B,
+                            double &C, double &scaleFactor);
+
     /// returns discretization at 3D index
-    void getBoundaryStencil(int idx, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
+    void getBoundaryStencil(int idx, double &W, double &E,
+                            double &S, double &N, double &F,
+                            double &B, double &C, double &scaleFactor);
+
     /// returns index of neighbours at (x,y,z)
-    void getNeighbours(int x, int y, int z, int &W, int &E, int &S, int &N, int &F, int &B);
+    void getNeighbours(int x, int y, int z, int &W, int &E,
+                       int &S, int &N, int &F, int &B);
+
     /// returns index of neighbours at 3D index
     void getNeighbours(int idx, int &W, int &E, int &S, int &N, int &F, int &B);
+
     /// returns type of boundary condition
     std::string getType() {return "Elliptic";}
+
     /// queries if a given (x,y,z) coordinate lies inside the domain
     inline bool isInside(int x, int y, int z) {
-        double xx = (x - (nr[0] - 1) / 2.0) * hr[0];
-        double yy = (y - (nr[1] - 1) / 2.0) * hr[1];
-        return ((xx * xx / (SemiMajor * SemiMajor) + yy * yy / (SemiMinor * SemiMinor) < 1) && z != 0 && z != nr[2] - 1);
+        double xx = - semiMajor_m + hr[0] * (x + 0.5);
+        double yy = - semiMinor_m + hr[1] * (y + 0.5);
+
+        bool isInsideEllipse = (xx * xx / (semiMajor_m * semiMajor_m) +
+                                yy * yy / (semiMinor_m * semiMinor_m) < 1);
+
+        return (isInsideEllipse && z > 0 && z < nr[2] - 1);
     }
 
     int getNumXY(int /*z*/) { return nxy_m; }
-    /// set semi-minor
-    void setSemiMinor(double sm) {SemiMinor = sm;}
-    /// set semi-major
-    void setSemiMajor(double sm) {SemiMajor = sm;}
 
-    /// calculates intersection 
-    void compute(Vector_t hr);
+
+    /// calculates intersection
+    void compute(Vector_t /*hr*/) {
+        throw OpalException("EllipticDomain::compute()", "This function is not available.");
+    }
+
     void compute(Vector_t hr, NDIndex<3> localId);
 
-    double getXRangeMin() { return -SemiMajor; }
-    double getXRangeMax() { return SemiMajor;  }
-    double getYRangeMin() { return -SemiMinor; }
-    double getYRangeMax() { return SemiMinor;  }
+    double getXRangeMin() { return -semiMajor_m; }
+    double getXRangeMax() { return semiMajor_m;  }
+    double getYRangeMin() { return -semiMinor_m; }
+    double getYRangeMax() { return semiMinor_m;  }
     double getZRangeMin() { return zMin_m; }
     double getZRangeMax() { return zMax_m; }
 
+    bool hasGeometryChanged() { return hasGeometryChanged_m; }
 
-    //TODO: ?
-    int getStartIdx() {return 0;}
 
-    bool hasGeometryChanged() { return hasGeometryChanged_m; }
+    void resizeMesh(Vector_t& origin, Vector_t& hr, const Vector_t& rmin,
+                    const Vector_t& rmax, double dh);
 
 private:
 
     /// Map from a single coordinate (x or y) to a list of intersection values with
     /// boundary.
-    typedef std::multimap<int, double> EllipticPointList;
+    typedef std::multimap<int, double> EllipticPointList_t;
 
     /// all intersection points with grid lines in X direction
-    EllipticPointList IntersectXDir;
+    EllipticPointList_t intersectXDir_m;
 
     /// all intersection points with grid lines in Y direction
-    EllipticPointList IntersectYDir;
+    EllipticPointList_t intersectYDir_m;
 
     /// mapping (x,y,z) -> idx
-    std::map<int, int> IdxMap;
+    std::map<int, int> idxMap_m;
 
     /// mapping idx -> (x,y,z)
-    std::map<int, int> CoordMap;
+    std::map<int, int> coordMap_m;
 
     /// semi-major of the ellipse
-    double SemiMajor;
+    double semiMajor_m;
+
     /// semi-minor of the ellipse
-    double SemiMinor;
+    double semiMinor_m;
+
     /// number of nodes in the xy plane (for this case: independent of the z coordinate)
     int nxy_m;
+
     /// interpolation type
-    int interpolationMethod;
+    int interpolationMethod_m;
+
     /// flag indicating if geometry has changed for the current time-step
     bool hasGeometryChanged_m;
 
     /// conversion from (x,y) to index in xy plane
     inline int toCoordIdx(int x, int y) { return y * nr[0] + x; }
+
     /// conversion from (x,y,z) to index on the 3D grid
     inline int getIdx(int x, int y, int z) {
-        if(isInside(x, y, z) && x >= 0 && y >= 0 && z >= 0)
-            return IdxMap[toCoordIdx(x, y)] + (z - 1) * nxy_m;
+        if (isInside(x, y, z))
+            return idxMap_m[toCoordIdx(x, y)] + (z - 1) * nxy_m;
         else
             return -1;
     }
+
     /// conversion from a 3D index to (x,y,z)
     inline void getCoord(int idx, int &x, int &y, int &z) {
         int ixy = idx % nxy_m;
-        int xy = CoordMap[ixy];
+        int xy = coordMap_m[ixy];
         int inr = (int)nr[0];
         x = xy % inr;
         y = (xy - x) / nr[0];
@@ -133,10 +162,23 @@ private:
     }
 
     /// different interpolation methods for boundary points
-    void constantInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
-    void linearInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
-    void quadraticInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
-
+    void constantInterpolation(int x, int y, int z,
+                               double &W, double &E, double &S,
+                               double &N, double &F, double &B,
+                               double &C, double &scaleFactor);
+
+    void linearInterpolation(int x, int y, int z, double &W,
+                             double &E, double &S, double &N,
+                             double &F, double &B, double &C,
+                             double &scaleFactor);
+
+    void quadraticInterpolation(int x, int y, int z, double &W,
+                                double &E, double &S, double &N,
+                                double &F, double &B, double &C,
+                                double &scaleFactor);
+
+    /// function to handle the open boundary condition in longitudinal direction
+    void robinBoundaryStencil(int z, double &F, double &B, double &C);
 };
 
 #endif //#ifdef ELLIPTICAL_DOMAIN_H

src/Solvers/IrregularDomain.h

@@ -129,6 +129,22 @@ public:
     virtual bool hasGeometryChanged() = 0;
     virtual ~IrregularDomain() {};
 
+    virtual void resizeMesh(Vector_t& origin, Vector_t& hr,
+                            const Vector_t& /*rmin*/, const Vector_t& /*rmax*/, double /*dh*/) {
+        double xmin = getXRangeMin();
+        double xmax = getXRangeMax();
+        double ymin = getYRangeMin();
+        double ymax = getYRangeMax();
+        double zmin = getZRangeMin();
+        double zmax = getZRangeMax();
+
+        origin = Vector_t(xmin, ymin , zmin);
+        Vector_t mymax = Vector_t(xmax, ymax , zmax);
+
+        for (int i = 0; i < 3; i++)
+            hr[i] = (mymax[i] - origin[i]) / nr[i];
+    };
+
 protected:
 
     // a irregular domain is always defined on a grid

src/Solvers/MGPoissonSolver.cpp

@@ -353,8 +353,10 @@ void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) {
     for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
         for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
             for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
+                NDIndex<3> l(Index(idx, idx), Index(idy, idy), Index(idz, idz));
                 if (bp_m->isInside(idx, idy, idz))
-                        RHS->getDataNonConst()[id++] = 4*M_PI*rho[idx][idy][idz].get()/scaleFactor;
+                        RHS->replaceGlobalValue(bp_m->getIdx(idx, idy, idz),
+                                                4.0 * M_PI * rho.localElement(l) / scaleFactor);
             }
         }
     }
@@ -508,6 +510,7 @@ void MGPoissonSolver::IPPLToMap3D(NDIndex<3> localId) {
             }
         }
     }
+
     int indexbase = 0;
     map_p = Teuchos::rcp(new TpetraMap_t(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(),
                                          &MyGlobalElements[0], NumMyElements, indexbase, comm_mp));

src/Solvers/MGPoissonSolver.h

@@ -125,8 +125,16 @@ public:
 
     void extrapolateLHS();
 
+    void resizeMesh(Vector_t& origin, Vector_t& hr, const Vector_t& rmin,
+                    const Vector_t& rmax, double dh)
+    {
+        bp_m->resizeMesh(origin, hr, rmin, rmax, dh);
+    }
+
     Inform &print(Inform &os) const;
 
+
+
 private:
 
     bool isMatrixfilled_m;
@@ -171,6 +179,7 @@ private:
 
     /// Map holding the processor distribution of data
     Teuchos::RCP<TpetraMap_t> map_p;
+
     /// communicator used by Trilinos
     Teuchos::RCP<const Comm_t> comm_mp;
 

src/Solvers/PoissonSolver.h

@@ -58,6 +58,11 @@ public:
     virtual void test(PartBunchBase<double, 3> *bunch) = 0 ;
     virtual ~PoissonSolver(){};
 
+    virtual void resizeMesh(Vector_t& /*origin*/, Vector_t& /*hr*/,
+                            const Vector_t& /*rmin*/, const Vector_t& /*rmax*/,
+                            double /*dh*/)
+    { };
+
 };
 
 inline Inform &operator<<(Inform &os, const PoissonSolver &/*fs*/) {