diff --git a/src/Classic/Algorithms/PartBunch.cpp b/src/Classic/Algorithms/PartBunch.cpp
index b475a680dd1b5f7604acef6ec212cf60b14f0d23..2be8fe9af747162e8cb6e4722b4b1cb65c3e7aa8 100644
--- a/src/Classic/Algorithms/PartBunch.cpp
+++ b/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());
diff --git a/src/Classic/Algorithms/PartBunchBase.h b/src/Classic/Algorithms/PartBunchBase.h
index 76807acfa037f60c17faccd72b4f51a3a7b43870..47b865bd6587ccaacf28de35295d892cb421313a 100644
--- a/src/Classic/Algorithms/PartBunchBase.h
+++ b/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
*/
diff --git a/src/Solvers/EllipticDomain.cpp b/src/Solvers/EllipticDomain.cpp
index 28b6e8f2780ca407ad43de1898c5edfac73781a6..d3365cd53028bd7542eceaadb4ad370342b7e714 100644
--- a/src/Solvers/EllipticDomain.cpp
+++ b/src/Solvers/EllipticDomain.cpp
@@ -1,6 +1,8 @@
//
// 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
@@ -38,52 +40,36 @@ EllipticDomain::EllipticDomain(Vector_t nr, Vector_t hr) {
setHr(hr);
}
-EllipticDomain::EllipticDomain(double semimajor, double semiminor, Vector_t nr, Vector_t hr, std::string interpl) {
- SemiMajor = semimajor;
- SemiMinor = semiminor;
+EllipticDomain::EllipticDomain(double semimajor, double semiminor, Vector_t nr,
+ Vector_t hr, std::string interpl)
+ : semiMajor_m(semimajor)
+ , semiMinor_m(semiminor)
+{
setNr(nr);
setHr(hr);
- if(interpl == "CONSTANT")
- interpolationMethod = CONSTANT;
- else if(interpl == "LINEAR")
- interpolationMethod = LINEAR;
- else if(interpl == "QUADRATIC")
- interpolationMethod = QUADRATIC;
+ if (interpl == "CONSTANT")
+ interpolationMethod_m = CONSTANT;
+ else if (interpl == "LINEAR")
+ interpolationMethod_m = LINEAR;
+ else if (interpl == "QUADRATIC")
+ interpolationMethod_m = QUADRATIC;
}
-EllipticDomain::EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr, std::string interpl) {
- SemiMajor = bgeom->getA();
- SemiMinor = bgeom->getB();
+EllipticDomain::EllipticDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr,
+ std::string interpl)
+ : EllipticDomain(bgeom->getA(), bgeom->getB(), nr, hr, interpl)
+{
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(BoundaryGeometry *bgeom, Vector_t nr, std::string interpl)
+ : EllipticDomain(bgeom, nr,
+ Vector_t(2.0 * bgeom->getA() / nr[0],
+ 2.0 * bgeom->getB() / nr[1],
+ (bgeom->getLength() - bgeom->getS()) / nr[2]),
+ interpl)
+{ }
EllipticDomain::~EllipticDomain() {
//nothing so far
@@ -93,164 +79,115 @@ EllipticDomain::~EllipticDomain() {
// for this geometry we only have to calculate the intersection on
// one x-y-plane
// for the moment we center the ellipse around the center of the grid
-void EllipticDomain::compute(Vector_t hr){
- //there is nothing to be done if the mesh spacings have not changed
- if(hr[0] == getHr()[0] && hr[1] == getHr()[1] && hr[2] == getHr()[2]) {
+void EllipticDomain::compute(Vector_t hr, NDIndex<3> localId){
+ // there is nothing to be done if the mesh spacings in x and y have not changed
+ if (hr[0] == getHr()[0] &&
+ hr[1] == getHr()[1])
+ {
hasGeometryChanged_m = false;
return;
}
+
setHr(hr);
hasGeometryChanged_m = true;
//reset number of points inside domain
nxy_m = 0;
// clear previous coordinate maps
- IdxMap.clear();
- CoordMap.clear();
+ idxMap_m.clear();
+ coordMap_m.clear();
//clear previous intersection points
- IntersectYDir.clear();
- IntersectXDir.clear();
+ intersectYDir_m.clear();
+ intersectXDir_m.clear();
// build a index and coordinate map
int idx = 0;
int x, y;
- for(x = 0; x < nr[0]; x++) {
- for(y = 0; y < nr[1]; y++) {
- if(isInside(x, y, 1)) {
- IdxMap[toCoordIdx(x, y)] = idx;
- CoordMap[idx++] = toCoordIdx(x, y);
+ /* we need to scan the full x-y-plane on all cores
+ * in order to figure out the number of valid
+ * grid points per plane --> otherwise we might
+ * get not unique global indices in the Tpetra::CrsMatrix
+ */
+ for (x = 0; x < nr[0]; ++x) {
+ for (y = 0; y < nr[1]; ++y) {
+ if (isInside(x, y, 1)) {
+ idxMap_m[toCoordIdx(x, y)] = idx;
+ coordMap_m[idx++] = toCoordIdx(x, y);
nxy_m++;
}
-
}
}
- switch(interpolationMethod) {
+ switch (interpolationMethod_m) {
case CONSTANT:
break;
case LINEAR:
case QUADRATIC:
- double smajsq = SemiMajor * SemiMajor;
- double sminsq = SemiMinor * SemiMinor;
+ double smajsq = semiMajor_m * semiMajor_m;
+ double sminsq = semiMinor_m * semiMinor_m;
double yd = 0.0;
double xd = 0.0;
double pos = 0.0;
- double mx = (nr[0] - 1) * hr[0] / 2.0;
- double my = (nr[1] - 1) * hr[1] / 2.0;
- //calculate intersection with the ellipse
- for(x = 0; x < nr[0]; x++) {
- pos = x * hr[0] - mx;
- if (pos <= -SemiMajor || pos >= SemiMajor)
+ // calculate intersection with the ellipse
+ for (x = localId[0].first(); x <= localId[0].last(); x++) {
+ pos = - semiMajor_m + hr[0] * (x + 0.5);
+ if (std::abs(pos) >= semiMajor_m)
{
- IntersectYDir.insert(std::pair<int, double>(x, 0));
- IntersectYDir.insert(std::pair<int, double>(x, 0));
- }else{
- yd = std::abs(sqrt(sminsq - sminsq * pos * pos / smajsq)); // + 0.5*nr[1]*hr[1]);
- IntersectYDir.insert(std::pair<int, double>(x, yd));
- IntersectYDir.insert(std::pair<int, double>(x, -yd));
+ intersectYDir_m.insert(std::pair<int, double>(x, 0));
+ intersectYDir_m.insert(std::pair<int, double>(x, 0));
+ } else {
+ yd = std::abs(std::sqrt(sminsq - sminsq * pos * pos / smajsq));
+ intersectYDir_m.insert(std::pair<int, double>(x, yd));
+ intersectYDir_m.insert(std::pair<int, double>(x, -yd));
}
}
- for(y = 0; y < nr[1]; y++) {
- pos = y * hr[1] - my;
- if (pos <= -SemiMinor || pos >= SemiMinor)
+ for (y = localId[0].first(); y < localId[1].last(); y++) {
+ pos = - semiMinor_m + hr[1] * (y + 0.5);
+ if (std::abs(pos) >= semiMinor_m)
{
- IntersectXDir.insert(std::pair<int, double>(y, 0));
- IntersectXDir.insert(std::pair<int, double>(y, 0));
- }else{
- xd = std::abs(sqrt(smajsq - smajsq * pos * pos / sminsq)); // + 0.5*nr[0]*hr[0]);
- IntersectXDir.insert(std::pair<int, double>(y, xd));
- IntersectXDir.insert(std::pair<int, double>(y, -xd));
+ intersectXDir_m.insert(std::pair<int, double>(y, 0));
+ intersectXDir_m.insert(std::pair<int, double>(y, 0));
+ } else {
+ xd = std::abs(std::sqrt(smajsq - smajsq * pos * pos / sminsq));
+ intersectXDir_m.insert(std::pair<int, double>(y, xd));
+ intersectXDir_m.insert(std::pair<int, double>(y, -xd));
}
}
}
}
-void EllipticDomain::compute(Vector_t hr, NDIndex<3> localId){
- //there is nothing to be done if the mesh spacings have not changed
- if(hr[0] == getHr()[0] && hr[1] == getHr()[1] && hr[2] == getHr()[2]) {
- hasGeometryChanged_m = false;
- return;
- }
- setHr(hr);
- hasGeometryChanged_m = true;
- //reset number of points inside domain
- nxy_m = 0;
-
- // clear previous coordinate maps
- IdxMap.clear();
- CoordMap.clear();
- //clear previous intersection points
- IntersectYDir.clear();
- IntersectXDir.clear();
-
- // build a index and coordinate map
- int idx = 0;
- int x, y;
-
- for(x = localId[0].first(); x<= localId[0].last(); x++) {
- for(y = localId[1].first(); y <= localId[1].last(); y++) {
- if(isInside(x, y, 1)) {
- IdxMap[toCoordIdx(x, y)] = idx;
- CoordMap[idx++] = toCoordIdx(x, y);
- nxy_m++;
- }
- }
- }
-
- switch(interpolationMethod) {
- case CONSTANT:
- break;
- case LINEAR:
- case QUADRATIC:
+void EllipticDomain::resizeMesh(Vector_t& origin, Vector_t& hr, const Vector_t& rmin,
+ const Vector_t& rmax, double dh)
+{
+ Vector_t mymax = Vector_t(0.0, 0.0, 0.0);
- double smajsq = SemiMajor * SemiMajor;
- double sminsq = SemiMinor * SemiMinor;
- double yd = 0.0;
- double xd = 0.0;
- double pos = 0.0;
- double mx = (nr[0] - 1) * hr[0] / 2.0;
- double my = (nr[1] - 1) * hr[1] / 2.0;
+ // apply bounding box increment dh, i.e., "BBOXINCR" input argument
+ double zsize = rmax[2] - rmin[2];
- //calculate intersection with the ellipse
- for(x = localId[0].first(); x <= localId[0].last(); x++) {
- pos = x * hr[0] - mx;
- if (pos <= -SemiMajor || pos >= SemiMajor)
- {
- IntersectYDir.insert(std::pair<int, double>(x, 0));
- IntersectYDir.insert(std::pair<int, double>(x, 0));
- }else{
- yd = std::abs(sqrt(sminsq - sminsq * pos * pos / smajsq)); // + 0.5*nr[1]*hr[1]);
- IntersectYDir.insert(std::pair<int, double>(x, yd));
- IntersectYDir.insert(std::pair<int, double>(x, -yd));
- }
+ setMinMaxZ(rmin[2] - zsize * (1.0 + dh),
+ rmax[2] + zsize * (1.0 + dh));
- }
+ origin = Vector_t(-semiMajor_m, -semiMinor_m, getMinZ());
+ mymax = Vector_t( semiMajor_m, semiMinor_m, getMaxZ());
- for(y = localId[0].first(); y < localId[1].last(); y++) {
- pos = y * hr[1] - my;
- if (pos <= -SemiMinor || pos >= SemiMinor)
- {
- IntersectXDir.insert(std::pair<int, double>(y, 0));
- IntersectXDir.insert(std::pair<int, double>(y, 0));
- }else{
- xd = std::abs(sqrt(smajsq - smajsq * pos * pos / sminsq)); // + 0.5*nr[0]*hr[0]);
- IntersectXDir.insert(std::pair<int, double>(y, xd));
- IntersectXDir.insert(std::pair<int, double>(y, -xd));
- }
- }
- }
+ for (int i = 0; i < 3; ++i)
+ hr[i] = (mymax[i] - origin[i]) / nr[i];
}
-void EllipticDomain::getBoundaryStencil(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
+void EllipticDomain::getBoundaryStencil(int x, int y, int z, double &W,
+ double &E, double &S, double &N,
+ double &F, double &B, double &C,
+ double &scaleFactor)
+{
scaleFactor = 1.0;
- // determine which interpolation method we use for points near the boundary
- switch(interpolationMethod) {
+ // determine which interpolation method we use for points near the boundary
+ switch (interpolationMethod_m) {
case CONSTANT:
constantInterpolation(x, y, z, W, E, S, N, F, B, C, scaleFactor);
break;
@@ -266,9 +203,10 @@ void EllipticDomain::getBoundaryStencil(int x, int y, int z, double &W, double &
assert(C > 0);
}
-void EllipticDomain::getBoundaryStencil(int idx, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
-
-
+void EllipticDomain::getBoundaryStencil(int idx, double &W, double &E,
+ double &S, double &N, double &F,
+ double &B, double &C, double &scaleFactor)
+{
int x = 0, y = 0, z = 0;
getCoord(idx, x, y, z);
getBoundaryStencil(x, y, z, W, E, S, N, F, B, C, scaleFactor);
@@ -283,113 +221,99 @@ void EllipticDomain::getNeighbours(int idx, int &W, int &E, int &S, int &N, int
}
-void EllipticDomain::getNeighbours(int x, int y, int z, int &W, int &E, int &S, int &N, int &F, int &B) {
-
- if(x > 0)
+void EllipticDomain::getNeighbours(int x, int y, int z, int &W, int &E,
+ int &S, int &N, int &F, int &B)
+{
+ if (x > 0)
W = getIdx(x - 1, y, z);
else
W = -1;
- if(x < nr[0] - 1)
+
+ if (x < nr[0] - 1)
E = getIdx(x + 1, y, z);
else
E = -1;
- if(y < nr[1] - 1)
+ if (y < nr[1] - 1)
N = getIdx(x, y + 1, z);
else
N = -1;
- if(y > 0)
+
+ if (y > 0)
S = getIdx(x, y - 1, z);
else
S = -1;
- if(z > 0)
+ if (z > 0)
F = getIdx(x, y, z - 1);
else
F = -1;
- if(z < nr[2] - 1)
+
+ if (z < nr[2] - 1)
B = getIdx(x, y, z + 1);
else
B = -1;
}
-void EllipticDomain::constantInterpolation(int x, int y, int z, double &WV, double &EV, double &SV, double &NV, double &FV, double &BV, double &CV, double &scaleFactor) {
-
+void EllipticDomain::constantInterpolation(int x, int y, int z, double &WV,
+ double &EV, double &SV, double &NV,
+ double &FV, double &BV, double &CV,
+ double &scaleFactor)
+{
scaleFactor = 1.0;
- WV = -1/(hr[0]*hr[0]);
- EV = -1/(hr[0]*hr[0]);
- NV = -1/(hr[1]*hr[1]);
- SV = -1/(hr[1]*hr[1]);
- FV = -1/(hr[2]*hr[2]);
- BV = -1/(hr[2]*hr[2]);
- CV = 2/(hr[0]*hr[0]) + 2/(hr[1]*hr[1]) + 2/(hr[2]*hr[2]);
+ WV = -1.0 / (hr[0] * hr[0]);
+ EV = -1.0 / (hr[0] * hr[0]);
+ NV = -1.0 / (hr[1] * hr[1]);
+ SV = -1.0 / (hr[1] * hr[1]);
+ FV = -1.0 / (hr[2] * hr[2]);
+ BV = -1.0 / (hr[2] * hr[2]);
+
+ CV = 2.0 / (hr[0] * hr[0])
+ + 2.0 / (hr[1] * hr[1])
+ + 2.0 / (hr[2] * hr[2]);
// we are a right boundary point
- if(!isInside(x + 1, y, z))
- EV= 0.0;
+ if (!isInside(x + 1, y, z))
+ EV = 0.0;
// we are a left boundary point
- if(!isInside(x - 1, y, z))
+ if (!isInside(x - 1, y, z))
WV = 0.0;
// we are a upper boundary point
- if(!isInside(x, y + 1, z))
+ if (!isInside(x, y + 1, z))
NV = 0.0;
// we are a lower boundary point
- if(!isInside(x, y - 1, z))
+ if (!isInside(x, y - 1, z))
SV = 0.0;
- if(z == 1 || z == nr[2] - 2) {
-
- // case where we are on the Robin BC in Z-direction
- // where we distinguish two cases
- // IFF: this values should not matter because they
- // never make it into the discretization matrix
- if(z == 1)
- FV = 0.0;
- else
- BV = 0.0;
-
- // add contribution of Robin discretization to center point
- // d the distance between the center of the bunch and the boundary
- //double cx = (x-(nr[0]-1)/2)*hr[0];
- //double cy = (y-(nr[1]-1)/2)*hr[1];
- //double cz = hr[2]*(nr[2]-1);
- //double d = sqrt(cx*cx+cy*cy+cz*cz);
- double d = hr[2] * (nr[2] - 1) / 2;
- CV += 2 / (d * hr[2]);
- //C += 2/((hr[2]*(nr[2]-1)/2.0) * hr[2]);
-
- // scale all stencil-points in z-plane with 0.5 (Robin discretization)
- WV /= 2.0;
- EV /= 2.0;
- NV /= 2.0;
- SV /= 2.0;
- CV /= 2.0;
- scaleFactor *= 0.5;
- }
-
+ // handle boundary condition in z direction
+ robinBoundaryStencil(z, FV, BV, CV);
}
-void EllipticDomain::linearInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
-
+void EllipticDomain::linearInterpolation(int x, int y, int z, double &W,
+ double &E, double &S, double &N,
+ double &F, double &B, double &C,
+ double &scaleFactor)
+{
scaleFactor = 1.0;
- double cx = x * hr[0] - (nr[0] - 1) * hr[0] / 2.0;
- double cy = y * hr[1] - (nr[1] - 1) * hr[1] / 2.0;
+ double cx = - semiMajor_m + hr[0] * (x + 0.5);
+ double cy = - semiMinor_m + hr[1] * (y + 0.5);
double dx = 0.0;
- std::multimap<int, double>::iterator it = IntersectXDir.find(y);
- if(cx < 0)
+ std::multimap<int, double>::iterator it = intersectXDir_m.find(y);
+
+ if (cx < 0)
++it;
dx = it->second;
double dy = 0.0;
- it = IntersectYDir.find(x);
- if(cy < 0)
+ it = intersectYDir_m.find(x);
+ if (cy < 0)
++it;
dy = it->second;
@@ -400,73 +324,55 @@ void EllipticDomain::linearInterpolation(int x, int y, int z, double &W, double
double ds = hr[1];
C = 0.0;
- //we are a right boundary point
- if(!isInside(x + 1, y, z)) {
- C += 1 / ((dx - cx) * de);
+ // we are a right boundary point
+ if (!isInside(x + 1, y, z)) {
+ C += 1.0 / ((dx - cx) * de);
E = 0.0;
} else {
- C += 1 / (de * de);
- E = -1 / (de * de);
+ C += 1.0 / (de * de);
+ E = -1.0 / (de * de);
}
- //we are a left boundary point
- if(!isInside(x - 1, y, z)) {
- C += 1 / ((std::abs(dx) - std::abs(cx)) * dw);
+ // we are a left boundary point
+ if (!isInside(x - 1, y, z)) {
+ C += 1.0 / ((std::abs(dx) - std::abs(cx)) * dw);
W = 0.0;
} else {
- C += 1 / (dw * dw);
- W = -1 / (dw * dw);
+ C += 1.0 / (dw * dw);
+ W = -1.0 / (dw * dw);
}
- //we are a upper boundary point
- if(!isInside(x, y + 1, z)) {
- C += 1 / ((dy - cy) * dn);
+ // we are a upper boundary point
+ if (!isInside(x, y + 1, z)) {
+ C += 1.0 / ((dy - cy) * dn);
N = 0.0;
} else {
- C += 1 / (dn * dn);
- N = -1 / (dn * dn);
+ C += 1.0 / (dn * dn);
+ N = -1.0 / (dn * dn);
}
- //we are a lower boundary point
- if(!isInside(x, y - 1, z)) {
- C += 1 / ((std::abs(dy) - std::abs(cy)) * ds);
+ // we are a lower boundary point
+ if (!isInside(x, y - 1, z)) {
+ C += 1.0 / ((std::abs(dy) - std::abs(cy)) * ds);
S = 0.0;
} else {
- C += 1 / (ds * ds);
- S = -1 / (ds * ds);
+ C += 1.0 / (ds * ds);
+ S = -1.0 / (ds * ds);
}
- F = -1 / (hr[2] * hr[2]);
- B = -1 / (hr[2] * hr[2]);
- C += 2 / (hr[2] * hr[2]);
-
// handle boundary condition in z direction
-/*
- if(z == 0 || z == nr[2] - 1) {
-
- // case where we are on the NEUMAN BC in Z-direction
- // where we distinguish two cases
- if(z == 0)
- F = 0.0;
- else
- B = 0.0;
-
- //hr[2]*(nr2[2]-1)/2 = radius
- double d = hr[2] * (nr[2] - 1) / 2;
- C += 2 / (d * hr[2]);
-
- W /= 2.0;
- E /= 2.0;
- N /= 2.0;
- S /= 2.0;
- C /= 2.0;
- scaleFactor *= 0.5;
-
- }
-*/
+ F = -1.0 / (hr[2] * hr[2]);
+ B = -1.0 / (hr[2] * hr[2]);
+ C += 2.0 / (hr[2] * hr[2]);
+ robinBoundaryStencil(z, F, B, C);
}
-void EllipticDomain::quadraticInterpolation(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
+void EllipticDomain::quadraticInterpolation(int x, int y, int z,
+ double &W, double &E, double &S,
+ double &N, double &F, double &B, double &C,
+ double &scaleFactor)
+{
+ scaleFactor = 1.0;
double cx = (x - (nr[0] - 1) / 2.0) * hr[0];
double cy = (y - (nr[1] - 1) / 2.0) * hr[1];
@@ -474,14 +380,14 @@ void EllipticDomain::quadraticInterpolation(int x, int y, int z, double &W, doub
// since every vector for elliptic domains has ALWAYS size 2 we
// can catch all cases manually
double dx = 0.0;
- std::multimap<int, double>::iterator it = IntersectXDir.find(y);
- if(cx < 0)
+ std::multimap<int, double>::iterator it = intersectXDir_m.find(y);
+ if (cx < 0)
++it;
dx = it->second;
double dy = 0.0;
- it = IntersectYDir.find(x);
- if(cy < 0)
+ it = intersectYDir_m.find(x);
+ if (cy < 0)
++it;
dy = it->second;
@@ -489,197 +395,80 @@ void EllipticDomain::quadraticInterpolation(int x, int y, int z, double &W, doub
double de = hr[0];
double dn = hr[1];
double ds = hr[1];
- W = 1.0;
- E = 1.0;
- N = 1.0;
- S = 1.0;
- F = 1.0;
- B = 1.0;
+
+ W = 0.0;
+ E = 0.0;
+ S = 0.0;
+ N = 0.0;
C = 0.0;
- //TODO: = cx+hr[0] > dx && cx > 0
- //if((x-nr[0]/2.0+1)*hr[0] > dx && cx > 0) {
- ////we are a right boundary point
- ////if(!isInside(x+1,y,z)) {
- //de = dx-cx;
- //}
-
- //if((x-nr[0]/2.0-1)*hr[0] < dx && cx < 0) {
- ////we are a left boundary point
- ////if(!isInside(x-1,y,z)) {
- //dw = std::abs(dx)-std::abs(cx);
- //}
-
- //if((y-nr[1]/2.0+1)*hr[1] > dy && cy > 0) {
- ////we are a upper boundary point
- ////if(!isInside(x,y+1,z)) {
- //dn = dy-cy;
- //}
-
- //if((y-nr[1]/2.0-1)*hr[1] < dy && cy < 0) {
- ////we are a lower boundary point
- ////if(!isInside(x,y-1,z)) {
- //ds = std::abs(dy)-std::abs(cy);
- //}
-
- //TODO: = cx+hr[0] > dx && cx > 0
- //if((x-nr[0]/2.0+1)*hr[0] > dx && cx > 0) {
- //we are a right boundary point
- if(!isInside(x + 1, y, z)) {
- de = dx - cx;
+ // we are a right boundary point
+ if (!isInside(x + 1, y, z)) {
+ double s = dx - cx;
+ C -= 2.0 * (s - 1.0) / (s * de);
E = 0.0;
+ W += (s - 1.0) / ((s + 1.0) * de);
+ } else {
+ C += 1.0 / (de * de);
+ E = -1.0 / (de * de);
}
- //if((x-nr[0]/2.0-1)*hr[0] < dx && cx < 0) {
- //we are a left boundary point
- if(!isInside(x - 1, y, z)) {
- dw = std::abs(dx) - std::abs(cx);
+ // we are a left boundary point
+ if (!isInside(x - 1, y, z)) {
+ double s = std::abs(dx) - std::abs(cx);
+ C -= 2.0 * (s - 1.0) / (s * de);
W = 0.0;
+ E += (s - 1.0) / ((s + 1.0) * de);
+ } else {
+ C += 1.0 / (dw * dw);
+ W = -1.0 / (dw * dw);
}
- //if((y-nr[1]/2.0+1)*hr[1] > dy && cy > 0) {
- //we are a upper boundary point
- if(!isInside(x, y + 1, z)) {
- dn = dy - cy;
+ // we are a upper boundary point
+ if (!isInside(x, y + 1, z)) {
+ double s = dy - cy;
+ C -= 2.0 * (s - 1.0) / (s * dn);
N = 0.0;
+ S += (s - 1.0) / ((s + 1.0) * dn);
+ } else {
+ C += 1.0 / (dn * dn);
+ N = -1.0 / (dn * dn);
}
- //if((y-nr[1]/2.0-1)*hr[1] < dy && cy < 0) {
- //we are a lower boundary point
- if(!isInside(x, y - 1, z)) {
- ds = std::abs(dy) - std::abs(cy);
+ // we are a lower boundary point
+ if (!isInside(x, y - 1, z)) {
+ double s = std::abs(dy) - std::abs(cy);
+ C -= 2.0 * (s - 1.0) / (s * dn);
S = 0.0;
+ N += (s - 1.0) / ((s + 1.0) * dn);
+ } else {
+ C += 1.0 / (ds * ds);
+ S = -1.0 / (ds * ds);
}
- //2/dw*(dw_de)
- W *= -1.0 / (dw * (dw + de));
- E *= -1.0 / (de * (dw + de));
- N *= -1.0 / (dn * (dn + ds));
- S *= -1.0 / (ds * (dn + ds));
- F = -1 / (hr[2] * (hr[2] + hr[2]));
- B = -1 / (hr[2] * (hr[2] + hr[2]));
-
- //TODO: problem when de,dw,dn,ds == 0
- //is NOT a regular BOUND PT
- C += 1 / de * 1 / dw;
- C += 1 / dn * 1 / ds;
- C += 1 / hr[2] * 1 / hr[2];
- scaleFactor = 0.5;
-
-
- //for regular gridpoints no problem with symmetry, just boundary
- //z direction is right
- //implement isLastInside(dir)
- //we have LOCAL x,y coordinates!
-
- /*
- if(dw != 0 && !wIsB)
- W = -1/dw * (dn+ds) * 2*hr[2];
- else
- W = 0;
- if(de != 0 && !eIsB)
- E = -1/de * (dn+ds) * 2*hr[2];
- else
- E = 0;
- if(dn != 0 && !nIsB)
- N = -1/dn * (dw+de) * 2*hr[2];
- else
- N = 0;
- if(ds != 0 && !sIsB)
- S = -1/ds * (dw+de) * 2*hr[2];
- else
- S = 0;
- F = -(dw+de)*(dn+ds)/hr[2];
- B = -(dw+de)*(dn+ds)/hr[2];
- */
-
- //if(dw != 0)
- //W = -2*hr[2]*(dn+ds)/dw;
- //else
- //W = 0;
- //if(de != 0)
- //E = -2*hr[2]*(dn+ds)/de;
- //else
- //E = 0;
- //if(dn != 0)
- //N = -2*hr[2]*(dw+de)/dn;
- //else
- //N = 0;
- //if(ds != 0)
- //S = -2*hr[2]*(dw+de)/ds;
- //else
- //S = 0;
- //F = -(dw+de)*(dn+ds)/hr[2];
- //B = -(dw+de)*(dn+ds)/hr[2];
-
- //// RHS scaleFactor for current 3D index
- //// Factor 0.5 results from the SW/quadratic extrapolation
- //scaleFactor = 0.5*(dw+de)*(dn+ds)*(2*hr[2]);
-
- // catch the case where a point lies on the boundary
- //FIXME: do this more elegant!
- //double m1 = dw*de;
- //double m2 = dn*ds;
- //if(de == 0)
- //m1 = dw;
- //if(dw == 0)
- //m1 = de;
- //if(dn == 0)
- //m2 = ds;
- //if(ds == 0)
- //m2 = dn;
- ////XXX: dn+ds || dw+de can be 0
- ////C = 2*(dn+ds)*(dw+de)/hr[2];
- //C = 2/hr[2];
- //if(dw != 0 || de != 0)
- //C *= (dw+de);
- //if(dn != 0 || ds != 0)
- //C *= (dn+ds);
- //if(dw != 0 || de != 0)
- //C += (2*hr[2])*(dn+ds)*(dw+de)/m1;
- //if(dn != 0 || ds != 0)
- //C += (2*hr[2])*(dw+de)*(dn+ds)/m2;
-
- //handle Neumann case
- //if(z == 0 || z == nr[2]-1) {
-
- //if(z == 0)
- //F = 0.0;
- //else
- //B = 0.0;
-
- ////neumann stuff
- //W = W/2.0;
- //E = E/2.0;
- //N = N/2.0;
- //S = S/2.0;
- //C /= 2.0;
-
- //scaleFactor /= 2.0;
- //}
-
// handle boundary condition in z direction
- if(z == 0 || z == nr[2] - 1) {
+ F = -1.0 / (hr[2] * hr[2]);
+ B = -1.0 / (hr[2] * hr[2]);
+ C += 2.0 / (hr[2] * hr[2]);
+ robinBoundaryStencil(z, F, B, C);
+}
+
+void EllipticDomain::robinBoundaryStencil(int z, double &F, double &B, double &C) {
+ if (z == 0 || z == nr[2] - 1) {
- // case where we are on the NEUMAN BC in Z-direction
+ // case where we are on the Robin BC in Z-direction
// where we distinguish two cases
- if(z == 0)
+ // IFF: this values should not matter because they
+ // never make it into the discretization matrix
+ if (z == 0)
F = 0.0;
else
B = 0.0;
- //C += 2/((hr[2]*(nr[2]-1)/2.0) * hr[2]);
- //hr[2]*(nr2[2]-1)/2 = radius
- double d = hr[2] * (nr[2] - 1) / 2;
- C += 2 / (d * hr[2]);
-
- W /= 2.0;
- E /= 2.0;
- N /= 2.0;
- S /= 2.0;
- C /= 2.0;
- scaleFactor /= 2.0;
-
+ // add contribution of Robin discretization to center point
+ // d the distance between the center of the bunch and the boundary
+ double d = 0.5 * hr[2] * (nr[2] - 1);
+ C += 2.0 / (d * hr[2]);
}
}
diff --git a/src/Solvers/EllipticDomain.h b/src/Solvers/EllipticDomain.h
index 691e82cffb148c0180103eef7fd127c843681fb2..e5478b5ae1fe0cdd90739b84cffa1c9232a5449b 100644
--- a/src/Solvers/EllipticDomain.h
+++ b/src/Solvers/EllipticDomain.h
@@ -1,6 +1,8 @@
//
// 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
diff --git a/src/Solvers/IrregularDomain.h b/src/Solvers/IrregularDomain.h
index 4569f663d9cfbd9a3bdb78ce6875b5b69a438fd0..eb00435e1772f8dca9f019ebc0848650c29e0047 100644
--- a/src/Solvers/IrregularDomain.h
+++ b/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
diff --git a/src/Solvers/MGPoissonSolver.cpp b/src/Solvers/MGPoissonSolver.cpp
index 57836bfd6ea0634b7c4b8ebb34c2525aefff80b8..ddcf129c7ebe3b38e3ae1374d1de41320d266817 100644
--- a/src/Solvers/MGPoissonSolver.cpp
+++ b/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));
diff --git a/src/Solvers/MGPoissonSolver.h b/src/Solvers/MGPoissonSolver.h
index 44bf05aca7449937087cb2a7c415721e9972b1fe..c2e1de9a9fe87dea819218f1b629574b29cd4234 100644
--- a/src/Solvers/MGPoissonSolver.h
+++ b/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;
diff --git a/src/Solvers/PoissonSolver.h b/src/Solvers/PoissonSolver.h
index e2dd5a350b3ec2c1f8a9c2f65e07f8c9592891d6..84aa754cee851194d88d2da0b47501c89b172ed4 100644
--- a/src/Solvers/PoissonSolver.h
+++ b/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*/) {