Mesh the whole domain; take into account the bunch offset; ippl to 3D mapping...

Mesh the whole domain; take into account the bunch offset; ippl to 3D mapping for multicore simulations are handled

classic/5.0/src/Algorithms/PartBunch.cpp

@@ -701,6 +701,10 @@ void PartBunch::computeSelfFields(int binNumber) {
     eg_m = Vector_t(0.0);
 
     if(fs_m->hasValidSolver()) {
+         /// Mesh the whole domain
+         if(fs_m->getFieldSolverType() == "SAAMG")
+             resizeMesh();
+
         /// Scatter charge onto space charge grid.
         this->Q *= this->dt;
         if(!interpolationCacheSet_m) {
@@ -893,8 +897,8 @@ void PartBunch::resizeMesh() {
     double xmax = fs_m->solver_m->getXRangeMax();
     double ymin = fs_m->solver_m->getYRangeMin();
     double ymax = fs_m->solver_m->getYRangeMax();
-    double zmin = rmin_m[2]; //fs_m->solver_m->getZRangeMin();
-    double zmax = rmax_m[2]; //fs_m->solver_m->getZRangeMax();
+    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]) {
@@ -915,12 +919,11 @@ void PartBunch::resizeMesh() {
         boundp();
         get_bounds(rmin_m, rmax_m);
     }
-    // extend domain with extra "ghost" point
-    Vector_t mymin = Vector_t(xmin, ymin , zmin-hr_m[2]);
-    Vector_t mymax = Vector_t(xmax, ymax , zmax+hr_m[2]);
+    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];
+        hr_m[i]   = (mymax[i] - mymin[i])/nr_m[i];
 
     getMesh().set_meshSpacing(&(hr_m[0]));
     getMesh().set_origin(mymin);
@@ -945,10 +948,9 @@ void PartBunch::computeSelfFields() {
     eg_m = Vector_t(0.0);
 
     if(fs_m->hasValidSolver()) {
-        //use the mesh that is already set
-        //if (fs_m->getFieldSolverType() == "SAAMG")
-        //   resizeMesh();
-	// INFOMSG("after resizeMesh" << hr_m << endl);
+        //mesh the whole domain
+        if(fs_m->getFieldSolverType() == "SAAMG")
+            resizeMesh();
 
         INFOMSG("mesh size" << hr_m << endl);
         //scatter charges onto grid
@@ -1220,6 +1222,9 @@ void PartBunch::computeSelfFields_cycl(double gamma) {
     eg_m = Vector_t(0.0);
 
     if(fs_m->hasValidSolver()) {
+        /// mesh the whole domain
+        if(fs_m->getFieldSolverType() == "SAAMG")
+            resizeMesh();
         /// scatter particles charge onto grid.
         this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());
 
@@ -1412,6 +1417,9 @@ void PartBunch::computeSelfFields_cycl(int bin) {
     double gamma = getBinGamma(bin);
 
     if(fs_m->hasValidSolver()) {
+        /// mesh the whole domain
+        if(fs_m->getFieldSolverType() == "SAAMG")
+            resizeMesh();
 
         /// scatter particles charge onto grid.
         this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());
@@ -2950,4 +2958,4 @@ bool PartBunch::WeHaveEnergyBins() {
         return dist_m->GetNumberOfEnergyBins() > 0;
     else
         return false;
-}
\ No newline at end of file
+}

src/Algorithms/ParallelTTracker.cpp

@@ -3252,4 +3252,4 @@ Vector_t ParallelTTracker::calcMeanP() const {
     }
     reduce(meanP, meanP, OpAddAssign());
     return meanP / Vector_t(itsBunch->getTotalNum());
-}
\ No newline at end of file
+}

src/Solvers/ArbitraryDomain.h

@@ -39,7 +39,7 @@ public:
 
     /// returns number of nodes in xy plane
     int getNumXY(int idz);
-    /// calculates intersection 
+    // calculates intersection  
     void Compute(Vector_t hr);
     // calculates intersection with rotated and shifted geometry 
     void Compute(Vector_t hr, NDIndex<3> localId);
@@ -93,7 +93,7 @@ private:
     std::map<int, int> numXZ;
 
     // Number of nodes in the xy plane (for this case: independent of the z coordinate)
-    int nxy_m;
+    int nxy_m[1000];
     // mapping (x,y,z) -> idxyz
     std::map<int, int> IdxMap;
     // Mapping idxyz -> (x,y,z)
@@ -106,6 +106,7 @@ private:
 
     Vector_t Geo_nr_m;
     Vector_t Geo_hr_m;
+    Vector_t geomCentroid_m;
     Vector_t minCoords_m;
     Vector_t maxCoords_m;
 

src/Solvers/EllipticDomain.cpp

@@ -29,6 +29,39 @@ EllipticDomain::EllipticDomain(double semimajor, double semiminor, Vector_t nr,
         interpolationMethod = QUADRATIC;
 }
 
+EllipticDomain::EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr, std::string interpl) {
+    SemiMajor = bgeom->getA();
+    SemiMinor = bgeom->getB();
+    setMinMaxZ(bgeom->getS(), bgeom->getLength());
+    setNr(nr);
+    setHr(hr);
+
+    if(interpl == "CONSTANT")
+        interpolationMethod = CONSTANT;
+    else if(interpl == "LINEAR")
+        interpolationMethod = LINEAR;
+    else if(interpl == "QUADRATIC")
+        interpolationMethod = QUADRATIC;
+}
+
+EllipticDomain::EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, std::string interpl) {
+    SemiMajor = bgeom->getA();
+    SemiMinor = bgeom->getB();
+    setMinMaxZ(bgeom->getS(), bgeom->getLength());
+    Vector_t hr_m; 
+    hr_m[0] = (getXRangeMax()-getXRangeMin())/nr[0];
+    hr_m[1] = (getYRangeMax()-getYRangeMin())/nr[1];
+    hr_m[2] = (getZRangeMax()-getZRangeMin())/nr[2];
+    setHr(hr_m);
+
+    if(interpl == "CONSTANT")
+        interpolationMethod = CONSTANT;
+    else if(interpl == "LINEAR")
+        interpolationMethod = LINEAR;
+    else if(interpl == "QUADRATIC")
+        interpolationMethod = QUADRATIC;
+}
+
 EllipticDomain::~EllipticDomain() {
     //nothing so far
 }
@@ -143,7 +176,6 @@ void EllipticDomain::Compute(Vector_t hr, NDIndex<3> localId){
                 CoordMap[idx++] = toCoordIdx(x, y);
                 nxy_m++;
             }
-
         }
     }
 

src/Solvers/EllipticDomain.h

@@ -7,6 +7,7 @@
 #include <string>
 #include <cmath>
 #include "IrregularDomain.h"
+#include "Structure/BoundaryGeometry.h"
 
 class EllipticDomain : public IrregularDomain {
 
@@ -14,6 +15,9 @@ 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(BoundaryGeometry *bgeom, Vector_t nr, std::string interpl);
+
     ~EllipticDomain();
 
     /// returns discretization at (x,y,z)

src/Solvers/MGPoissonSolver.cpp

 #ifdef HAVE_SAAMG_SOLVER
-#define DBG_STENCIL
+//#define DBG_STENCIL
 
 #include "Solvers/MGPoissonSolver.h"
 
@@ -10,6 +10,7 @@
 //#include "RectangularDomain.h"
 
 #include "Algorithms/PartBunch.h"
+#include "Track/Track.h"
 #include "Physics/Physics.h"
 #include "Attributes/Attributes.h"
 #include "ValueDefinitions/RealVariable.h"
@@ -94,6 +95,7 @@ MGPoissonSolver::MGPoissonSolver ( PartBunch &beam,
     forcePreconditionerRecomputation_m = false;
 
     hasParallelDecompositionChanged_m = true;
+    repartFreq_m = 1000; 
     useRCB_m = false;
     if(Ippl::Info->getOutputLevel() > 1)
         verbose_m = true;
@@ -104,7 +106,7 @@ MGPoissonSolver::MGPoissonSolver ( PartBunch &beam,
     currentGeometry = geometries_m[0];
     if(currentGeometry->getFilename() == "") {
         if(currentGeometry->getTopology() == "ELLIPTIC"){
-            bp = new EllipticDomain(currentGeometry->getA(), currentGeometry->getB(), orig_nr_m, hr_m, interpl);
+            bp = new EllipticDomain(currentGeometry, orig_nr_m, hr_m, interpl);
 	} else if (currentGeometry->getTopology() == "BOXCORNER") {
             bp = new BoxCornerDomain(currentGeometry->getA(), currentGeometry->getB(), currentGeometry->getC(), currentGeometry->getLength(),currentGeometry->getL1(), currentGeometry->getL2(), orig_nr_m, hr_m, interpl);
             bp->Compute(itsBunch_m->get_hr());
@@ -112,7 +114,7 @@ MGPoissonSolver::MGPoissonSolver ( PartBunch &beam,
             throw OpalException("MGPoissonSolver::MGPoissonSolver",
                                 "Geometry not known");
         }
-    } else
+    } else 
 	bp = new ArbitraryDomain(currentGeometry, orig_nr_m, hr_m, interpl);
 
     Map = 0;
@@ -120,6 +122,7 @@ MGPoissonSolver::MGPoissonSolver ( PartBunch &beam,
     LHS = Teuchos::null;
     RHS = Teuchos::null;
     MLPrec = 0;
+    prec_m = Teuchos::null;
 
     numBlocks_m = Options::numBlocks;
     recycleBlocks_m = Options::recycleBlocks;
@@ -169,6 +172,8 @@ void MGPoissonSolver::deletePtr() {
     A = Teuchos::null;
     LHS = Teuchos::null;
     RHS = Teuchos::null;
+    if(Map) delete Map;
+    if(MLPrec) delete MLPrec; MLPrec=0;
     prec_m = Teuchos::null;
 }
 
@@ -177,13 +182,11 @@ void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr, double zshift)
 }
 
 void MGPoissonSolver::computeMap(NDIndex<3> localId) {
-    if (hasParallelDecompositionChanged_m){
-        if (Map !=0 ) delete Map;
+    if (itsBunch_m->getLocalTrackStep()%repartFreq_m == 0){
+        deletePtr(); 
         if(useRCB_m)
             redistributeWithRCB(localId);
-        else if ( bp->getType() == "Geometric" )
-            IPPLToMap3DGeo(localId);
-        else if ( bp->getType() == "Elliptic" )
+        else
             IPPLToMap3D(localId);
 
         extrapolateLHS();
@@ -255,17 +258,14 @@ void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) {
     nr_m[1] = orig_nr_m[1];
     nr_m[2] = orig_nr_m[2];
 
-    //bp->setMinMaxZ(itsBunch_m->get_origin()[2], itsBunch_m->get_maxExtend()[2]);
     bp->setGlobalMeanR(itsBunch_m->getGlobalMeanR());
     bp->setGlobalToLocalQuaternion(itsBunch_m->getGlobalToLocalQuaternion());
     bp->setNr(nr_m);
     NDIndex<3> localId = layout_m->getLocalNDIndex();
 
     IpplTimings::startTimer(FunctionTimer1_m);
-    if ( bp->getType() == "Geometric" ) {
+    if(!itsBunch_m->getLocalTrackStep())
         bp->Compute(hr, localId);
-    } else if ( bp->getType() == "Elliptic"  )
-        bp->Compute(hr);
     IpplTimings::stopTimer(FunctionTimer1_m);
 
     // Define the Map
@@ -274,8 +274,8 @@ void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) {
     IpplTimings::stopTimer(FunctionTimer2_m);
 
     // Allocate the RHS with the new Epetra Map
-    if (Teuchos::is_null(RHS))
-      RHS = rcp(new Epetra_Vector(*Map));
+    if (Teuchos::is_null(RHS)) 
+        RHS = rcp(new Epetra_Vector(*Map));
     RHS->PutScalar(0.0);
     // get charge densities from IPPL field and store in Epetra vector (RHS)
     IpplTimings::startTimer(FunctionTimer3_m);
@@ -293,23 +293,21 @@ void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) {
 
     // build discretization matrix
     IpplTimings::startTimer(FunctionTimer4_m);
-    if (Teuchos::is_null(A))
+    if(Teuchos::is_null(A))
         A = rcp(new Epetra_CrsMatrix(Copy, *Map,  7, true));
     ComputeStencil(hr, RHS);
     IpplTimings::stopTimer(FunctionTimer4_m);
+
 #ifdef DBG_STENCIL
     EpetraExt::RowMatrixToMatlabFile("DiscrStencil.dat", *A);
 #endif
 
     IpplTimings::startTimer(FunctionTimer5_m);
-    if ( MLPrec == 0 ){
+    if(!MLPrec) {
         MLPrec = new ML_Epetra::MultiLevelPreconditioner(*A, MLList_m);
     } else if (precmode_m == REUSE_HIERARCHY) {
         MLPrec->ReComputePreconditioner();
     } else if (precmode_m == REUSE_PREC){
-    } else if (precmode_m == STD_PREC) {
-        delete MLPrec; MLPrec = 0;
-        MLPrec = new ML_Epetra::MultiLevelPreconditioner(*A, MLList_m);
     }
     IpplTimings::stopTimer(FunctionTimer5_m);
 
@@ -320,7 +318,8 @@ void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) {
     problem_ptr->setOperator(A);
     problem_ptr->setLHS(LHS);
     problem_ptr->setRHS(RHS);
-    prec_m = Teuchos::rcp ( new Belos::EpetraPrecOp ( rcp(MLPrec,false)));
+    if(Teuchos::is_null(prec_m))
+        prec_m = Teuchos::rcp ( new Belos::EpetraPrecOp ( rcp(MLPrec,false)));
     problem_ptr->setLeftPrec(prec_m);
     solver_ptr->setProblem( problem_ptr);
     if (!problem_ptr->isProblemSet()) {
@@ -376,8 +375,13 @@ void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) {
         }
     }
     IpplTimings::stopTimer(FunctionTimer8_m);
-    if ( hasParallelDecompositionChanged_m )
-        deletePtr();
+
+    if(itsBunch_m->getLocalTrackStep()+1 == (long long)Track::block->localTimeSteps.front()) {
+        A = Teuchos::null;
+        LHS = Teuchos::null;
+        RHS = Teuchos::null;
+        prec_m = Teuchos::null;
+    }	
 }
 
 
@@ -456,21 +460,6 @@ void MGPoissonSolver::IPPLToMap3D(NDIndex<3> localId) {
     Map = new Epetra_Map(-1, NumMyElements, &MyGlobalElements[0], 0, Comm);
 }
 
-void MGPoissonSolver::IPPLToMap3DGeo(NDIndex<3> localId) {
-    int NumMyElements = 0;
-    std::vector<int> MyGlobalElements;
-
-    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++) {
-                    MyGlobalElements.push_back(bp->getIdx(idx, idy, idz));
-                    NumMyElements++;
-            }
-	}
-    }
-    Map = new Epetra_Map(-1, NumMyElements, &MyGlobalElements[0], 0, Comm);
-}
-
 void MGPoissonSolver::ComputeStencil(Vector_t hr, Teuchos::RCP<Epetra_Vector> RHS) {
 
     A->PutScalar(0.0);
@@ -569,4 +558,4 @@ Inform &MGPoissonSolver::print(Inform &os) const {
     return os;
 }
 
-#endif /* HAVE_SAAMG_SOLVER */
\ No newline at end of file
+#endif /* HAVE_SAAMG_SOLVER */

src/Solvers/MGPoissonSolver.h

@@ -67,6 +67,10 @@ namespace ML_Epetra {
 
 #pragma GCC diagnostic pop
 
+// using Teuchos::RCP;
+// using Teuchos::rcp;
+// using namespace ML_Epetra;
+// using namespace Isorropia;
 //////////////////////////////////////////////////////////////
 
 typedef UniformCartesian<3, double> Mesh_t;
@@ -141,6 +145,7 @@ private:
     bool hasGeometryChanged_m;
     /// flag is set when OPAL changed decomposition of mesh
     bool hasParallelDecompositionChanged_m;
+    int repartFreq_m;
     /// flag specifying if problem is redistributed with RCB
     bool useRCB_m;
     /// flag specifying if we are verbose
@@ -248,7 +253,6 @@ private:
     /// converts IPPL grid to a 3D Epetra_Map
     /// \param localId local IPPL grid node indices
     void IPPLToMap3D(NDIndex<3> localId);
-    void IPPLToMap3DGeo(NDIndex<3> localId);
 
     /** returns a discretized stencil that has Neumann BC in z direction and
      * Dirichlet BC on the surface of a specified geometry
@@ -360,4 +364,4 @@ int main(int argc, char *argv[]) {
 
 #endif /* #ifndef MG_POISSON_SOLVER_H_ */
 
-#endif /* #ifdef HAVE_SAAMG_SOLVER */
\ No newline at end of file
+#endif /* #ifdef HAVE_SAAMG_SOLVER */