Commit 908438c0 authored by frey_m's avatar frey_m
Browse files

yMerge branch '543-saamg-enable-rectangulardomain' of gitlab.psi.ch:OPAL/src...

yMerge branch '543-saamg-enable-rectangulardomain' of gitlab.psi.ch:OPAL/src into 565-follow-up-from-resolve-saamg-enable-rectangulardomain

+ add base class for regular domains
parents 8451cf66 f194d79b
......@@ -19,6 +19,7 @@ if (ENABLE_SAAMG_SOLVER)
IrregularDomain.cpp
MGPoissonSolver.cpp
RectangularDomain.cpp
RegularDomain.cpp
)
list (APPEND HDRS
......@@ -28,6 +29,7 @@ if (ENABLE_SAAMG_SOLVER)
IrregularDomain.h
MGPoissonSolver.h
RectangularDomain.h
RegularDomain.h
)
endif ()
......
......@@ -37,7 +37,7 @@
EllipticDomain::EllipticDomain(BoundaryGeometry *bgeom, IntVector_t nr, Vector_t hr,
std::string interpl)
: IrregularDomain(nr, hr, interpl)
: RegularDomain(nr, hr, interpl)
{
Vector_t min(-bgeom->getA(), -bgeom->getB(), bgeom->getS());
Vector_t max( bgeom->getA(), bgeom->getB(), bgeom->getS() + bgeom->getLength());
......@@ -50,13 +50,6 @@ EllipticDomain::~EllipticDomain() {
//nothing so far
}
int EllipticDomain::getNumXY() const {
return nxy_m;
}
// 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
......@@ -72,7 +65,6 @@ void EllipticDomain::compute(Vector_t hr, NDIndex<3> localId) {
setHr(hr);
hasGeometryChanged_m = true;
//reset number of points inside domain
nxy_m = 0;
// clear previous coordinate maps
idxMap_m.clear();
......@@ -85,6 +77,8 @@ void EllipticDomain::compute(Vector_t hr, NDIndex<3> localId) {
int idx = 0;
int x, y;
int nxy = 0;
/* 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
......@@ -95,11 +89,13 @@ void EllipticDomain::compute(Vector_t hr, NDIndex<3> localId) {
if (isInside(x, y, 1)) {
idxMap_m[toCoordIdx(x, y)] = idx;
coordMap_m[idx++] = toCoordIdx(x, y);
nxy_m++;
nxy++;
}
}
}
setNumXY(nxy);
switch (interpolationMethod_m) {
case CONSTANT:
break;
......@@ -142,24 +138,6 @@ void EllipticDomain::compute(Vector_t hr, NDIndex<3> localId) {
}
}
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);
// apply bounding box increment dh, i.e., "BBOXINCR" input argument
double zsize = rmax[2] - rmin[2];
setMinMaxZ(rmin[2] - zsize * (1.0 + dh),
rmax[2] + zsize * (1.0 + dh));
origin = Vector_t(getXRangeMin(), getYRangeMin(), getMinZ());
mymax = Vector_t(getXRangeMax(), getYRangeMax(), getMaxZ());
for (int i = 0; i < 3; ++i)
hr[i] = (mymax[i] - origin[i]) / nr_m[i];
}
void EllipticDomain::constantInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor) const
{
......
......@@ -32,11 +32,11 @@
#include <map>
#include <string>
#include <cmath>
#include "IrregularDomain.h"
#include "Solvers/RegularDomain.h"
#include "Structure/BoundaryGeometry.h"
#include "Utilities/OpalException.h"
class EllipticDomain : public IrregularDomain {
class EllipticDomain : public RegularDomain {
public:
EllipticDomain(BoundaryGeometry *bgeom, IntVector_t nr,
......@@ -44,8 +44,6 @@ public:
~EllipticDomain();
int getNumXY() const override;
/// queries if a given (x,y,z) coordinate lies inside the domain
bool isInside(int x, int y, int z) const {
double xx = getXRangeMin() + hr_m[0] * (x + 0.5);
......@@ -60,9 +58,6 @@ public:
/// calculates intersection
void compute(Vector_t hr, NDIndex<3> localId);
void resizeMesh(Vector_t& origin, Vector_t& hr, const Vector_t& rmin,
const Vector_t& rmax, double dh) override;
private:
/// Map from a single coordinate (x or y) to a list of intersection values with
......@@ -75,19 +70,16 @@ private:
/// all intersection points with grid lines in Y direction
EllipticPointList_t intersectYDir_m;
/// number of nodes in the xy plane (for this case: independent of the z coordinate)
int nxy_m;
/// conversion from (x,y) to index in xy plane
int toCoordIdx(int x, int y) const { return y * nr_m[0] + x; }
/// conversion from (x,y,z) to index on the 3D grid
int indexAccess(int x, int y, int z) const {
return idxMap_m.at(toCoordIdx(x, y)) + z * nxy_m;
return idxMap_m.at(toCoordIdx(x, y)) + z * getNumXY();
}
int coordAccess(int idx) const {
int ixy = idx % nxy_m;
int ixy = idx % getNumXY();
return coordMap_m.at(ixy);
}
......
......@@ -40,6 +40,7 @@
#include "ArbitraryDomain.h"
#include "EllipticDomain.h"
#include "BoxCornerDomain.h"
#include "RectangularDomain.h"
#include "Track/Track.h"
#include "Physics/Physics.h"
......@@ -127,6 +128,11 @@ MGPoissonSolver::MGPoissonSolver ( PartBunch *beam,
currentGeometry->getL2(),
orig_nr_m, hr_m, interpl));
bp_m->compute(itsBunch_m->get_hr(), layout_m->getLocalNDIndex());
} else if (currentGeometry->getTopology() == "RECTANGULAR") {
bp_m = std::unique_ptr<IrregularDomain>(
new RectangularDomain(currentGeometry->getA(),
currentGeometry->getB(),
orig_nr_m, hr_m));
} else {
throw OpalException("MGPoissonSolver::MGPoissonSolver",
"Geometry not known");
......
//
// Class RectangularDomain
// :FIXME: add brief description
// This class provides a rectangular beam pipe. The mesh adapts to the bunch
// in longitudinal direction.
//
// Copyright (c) 2008, Yves Ineichen, ETH Zürich,
// 2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland
......@@ -29,8 +30,7 @@
#include "Utilities/OpalException.h"
RectangularDomain::RectangularDomain(double a, double b, IntVector_t nr, Vector_t hr)
: IrregularDomain(nr, hr, "CONSTANT")
, nxy_m(nr[0] * nr[1])
: RegularDomain(nr, hr, "CONSTANT")
{
setRangeMin(Vector_t(-a, -b, getMinZ()));
setRangeMax(Vector_t( a, b, getMaxZ()));
......@@ -43,669 +43,36 @@ void RectangularDomain::compute(Vector_t hr, NDIndex<3> /*localId*/){
void RectangularDomain::constantInterpolation(int x, int y, int z, StencilValue_t& value,
double &scaleFactor) const
{
//scaleFactor = 1.0;
//W = -1/hr_m_m[0]*1/hr_m[0];
//E = -1/hr_m[0]*1/hr_m[0];
//N = -1/hr_m[1]*1/hr_m[1];
//S = -1/hr_m[1]*1/hr_m[1];
//F = -1/hr_m[2]*1/hr_m[2];
//B = -1/hr_m[2]*1/hr_m[2];
//C = 2/hr_m[0]*1/hr_m[0] + 2/hr_m[1]*1/hr_m[1] + 2/hr_m[2]*1/hr_m[2];
scaleFactor = hr_m[0] * hr_m[1] * hr_m[2];
value.west = -hr_m[1] * hr_m[2] / hr_m[0];
value.east = -hr_m[1] * hr_m[2] / hr_m[0];
value.north = -hr_m[0] * hr_m[2] / hr_m[1];
value.south = -hr_m[0] * hr_m[2] / hr_m[1];
value.front = -hr_m[0] * hr_m[1] / hr_m[2];
value.back = -hr_m[0] * hr_m[1] / hr_m[2];
value.center = 2 * hr_m[1] * hr_m[2] / hr_m[0] + 2 * hr_m[0] * hr_m[2] / hr_m[1] + 2 * hr_m[0] * hr_m[1] / hr_m[2];
if(!isInside(x + 1, y, z))
scaleFactor = 1.0;
value.west = -1.0 / (hr[0] * hr[0]);
value.east = -1.0 / (hr[0] * hr[0]);
value.north = -1.0 / (hr[1] * hr[1]);
value.south = -1.0 / (hr[1] * hr[1]);
value.front = -1.0 / (hr[2] * hr[2]);
value.back = -1.0 / (hr[2] * hr[2]);
value.center = 2.0 / (hr[0] * hr[0])
+ 2.0 / (hr[1] * hr[1])
+ 2.0 / (hr[2] * hr[2]);
if (!isInside(x + 1, y, z))
value.east = 0.0; //we are a right boundary point
if(!isInside(x - 1, y, z))
if (!isInside(x - 1, y, z))
value.west = 0.0; //we are a left boundary point
if(!isInside(x, y + 1, z))
if (!isInside(x, y + 1, z))
value.north = 0.0; //we are a upper boundary point
if(!isInside(x, y - 1, z))
if (!isInside(x, y - 1, z))
value.south = 0.0; //we are a lower boundary point
//dirichlet
if(z == 0)
if (z == 0)
value.front = 0.0;
if(z == nr_m[2] - 1)
value.back = 0.0;
if(false /*z == 0 || z == nr_m[2]-1*/) {
//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 == 0)
value.front = 0.0;
else
value.back = 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_m[0]-1)/2)*hr_m[0];
//double cy = (y-(nr_m[1]-1)/2)*hr_m[1];
//double cz = hr_m[2]*(nr_m[2]-1);
//double d = sqrt(cx*cx+cy*cy+cz*cz);
double d = hr_m[2] * (nr_m[2] - 1) / 2;
//cout << cx << " " << cy << " " << cz << " " << 2/(d*hr_m[2]) << endl;
value.center += 2 / (d * hr_m[2]);
//C += 2/((hr_m[2]*(nr_m[2]-1)/2.0) * hr_m[2]);
//scale all stencil-points in z-plane with 0.5 (Neumann discretization)
value.west /= 2.0;
value.east /= 2.0;
value.north /= 2.0;
value.south /= 2.0;
value.center /= 2.0;
scaleFactor *= 0.5;
}
//simple check if center value of stencil is positive
#ifdef DEBUG
if(value.center <= 0)
throw OpalException("RectangularDomain::getBoundaryStencil",
"Stencil C is <= 0! This case should never occure!");
#endif
}
/*
void MGPoissonSolver::getNeighbours(const int idx, int& W, int& E, int& S, int& N, int& F, int& B, int& numOutOfDomain)
{
int ixy, iz, iy, ix;
int nx = nr_m[0];
int ny = nr_m[1];
int nz = nr_m[2];
int nxy = nx*ny;
ixy = idx % nxy;
ix = ixy % nx;
iy = (ixy - ix) / nx;
iz = (idx - ixy) / nxy;
numOutOfDomain = 0;
if (iz == 0)
{F = -1; numOutOfDomain++;}
else
F = idx - nxy;
if (iz == nz - 1)
{B = -1; numOutOfDomain++;}
else
B = idx + nxy;
if (ix == 0)
{W = -1; numOutOfDomain++;}
else
W = ixy - 1;
if (ix == nx - 1)
{E = -1; numOutOfDomain++;}
else
E = ixy + 1;
if (iy == 0)
{S = -1; numOutOfDomain++;}
else
S = ixy - nx;
if (iy == ny - 1)
{N = -1; numOutOfDomain++;}
else
N = ixy + nx;
int cor = iz * nxy;
switch(W) {
case -1: break;
default: W += cor;
}
switch(E) {
case -1: break;
default: E += cor;
}
switch(S) {
case -1: break;
default: S += cor;
}
switch(N) {
case -1: break;
default: N += cor;
}
}
//at the moment this stencil has neumann everywhere except on the z=0 plane (Dirichlet)
//this code is experimental/not-optimized/not-final at the moment
inline Epetra_CrsMatrix* MGPoissonSolver::stencil3DOneSidedDirichlet(Vector_t hr)
{
Epetra_CrsMatrix* Matrix = new Epetra_CrsMatrix(Copy, *Map, 7);
Vector_t hr2 = hr*hr;
double a = 2*(1/hr2[0]+1/hr2[1]+1/hr2[2]);
double b = -1/hr2[0];
double d = -1/hr2[1];
double f = -1/hr2[2];
int NumMyElements = Map->NumMyElements();
int* MyGlobalElements = Map->MyGlobalElements();
int W, E, S, N, F, B, numout;
std::vector<double> Values(6);
std::vector<int> Indices(6);
for (int i = 0 ; i < NumMyElements ; ++i)
{
getNeighbours(MyGlobalElements[i], W, E, S, N, F, B, numout);
int NumEntries = 0;
double diag = a;
if(F != -1) {
switch(numout) {
case 3: {
if(W == -1) {
Indices[NumEntries] = E;
Values[NumEntries++] = 0.25*b;
}
if(E == -1) {
Indices[NumEntries] = W;
Values[NumEntries++] = 0.25*b;
}
if(N == -1) {
Indices[NumEntries] = S;
Values[NumEntries++] = 0.25*d;
}
if(S == -1) {
Indices[NumEntries] = N;
Values[NumEntries++] = 0.25*d;
}
if(B == -1) {
Indices[NumEntries] = F;
Values[NumEntries++] = 0.25*f;
}
if(F == -1) {
Indices[NumEntries] = B;
Values[NumEntries++] = 0.25*f;
}
diag *= 0.125;
if(NumEntries != 3)
cout << "ERROR: while calc corners in stencil" << endl;
break;
}
case 2: {
if(W != -1 && E != -1) {
Indices[NumEntries] = W;
Values[NumEntries++] = 0.25*b;
Indices[NumEntries] = E;
Values[NumEntries++] = 0.25*b;
}
if(E == -1 && W != -1) {
Indices[NumEntries] = W;
Values[NumEntries++] = 0.5*b;
}
if(W == -1 && E != -1) {
Indices[NumEntries] = E;
Values[NumEntries++] = 0.5*b;
}
if(N != -1 && S != -1) {
Indices[NumEntries] = N;
Values[NumEntries++] = 0.25*d;
Indices[NumEntries] = S;
Values[NumEntries++] = 0.25*d;
}
if(S == -1 && N != -1) {
Indices[NumEntries] = N;
Values[NumEntries++] = 0.5*d;
}
if(N == -1 && S != -1) {
Indices[NumEntries] = S;
Values[NumEntries++] = 0.5*d;
}
if(B != -1 && F != -1) {
Indices[NumEntries] = B;
Values[NumEntries++] = 0.25*f;
Indices[NumEntries] = F;
Values[NumEntries++] = 0.25*f;
}
if(B == -1 && F != -1) {
Indices[NumEntries] = F;
Values[NumEntries++] = 0.5*f;
}
if(F == -1 && B != -1) {
Indices[NumEntries] = B;
Values[NumEntries++] = 0.5*f;
}
diag *= 0.25;
if(NumEntries != 4)
cout << "ERROR: while calc edge stencil elements" << endl;
break;
}
case 1: {
if(W != -1 && E != -1) {
Indices[NumEntries] = W;
Values[NumEntries++] = 0.5*b;
Indices[NumEntries] = E;
Values[NumEntries++] = 0.5*b;
}
if(E == -1 && W != -1) {
Indices[NumEntries] = W;
Values[NumEntries++] = b;
}
if(W == -1 && E != -1) {
Indices[NumEntries] = E;
Values[NumEntries++] = b;
}
if(N != -1 && S != -1) {
Indices[NumEntries] = S;
Values[NumEntries++] = 0.5*d;
Indices[NumEntries] = N;
Values[NumEntries++] = 0.5*d;
}
if(S == -1 && N != -1) {
Indices[NumEntries] = N;
Values[NumEntries++] = d;
}
if(N == -1 && S != -1) {
Indices[NumEntries] = S;
Values[NumEntries++] = d;
}
if(B != -1 && F != -1) {
Indices[NumEntries] = F;
Values[NumEntries++] = 0.5*f;
Indices[NumEntries] = B;
Values[NumEntries++] = 0.5*f;
}
if(B == -1 && F != -1) {
Indices[NumEntries] = F;
Values[NumEntries++] = f;
}
if(F == -1 && B != -1) {
Indices[NumEntries] = B;
Values[NumEntries++] = f;
}
diag *= 0.5;
if(NumEntries != 5)
cout << "ERROR: calc surface elements of stencil" << endl;
break;
}
case 0: {
if(W != -1) {
Indices[NumEntries] = W;
Values[NumEntries++] = b;
} if(E != -1) {
Indices[NumEntries] = E;
Values[NumEntries++] = b;
} if(S != -1) {
Indices[NumEntries] = S;
Values[NumEntries++] = d;
} if(N != -1) {
Indices[NumEntries] = N;
Values[NumEntries++] = d;
} if(F != -1) {
Indices[NumEntries] = F;
Values[NumEntries++] = f;
} if(B != -1) {
Indices[NumEntries] = B;
Values[NumEntries++] = f;
}
break;
}
default: {
cout << "ERROR CREATING THE STENCIL: points out of domain (" << numout << ") is >3 OR <0" << endl;
}
}
} else {
//F == -1 --> this is our dirichlet boundary surface
switch(numout) {
case 3: {
if(W == -1) {
Indices[NumEntries] = E;
Values[NumEntries++] = 0.5*b;
}
if(E == -1) {
Indices[NumEntries] = W;
Values[NumEntries++] = 0.5*b;
}
if(N == -1) {
Indices[NumEntries] = S;
Values[NumEntries++] = 0.5*d;
}
if(S == -1) {
Indices[NumEntries] = N;
Values[NumEntries++] = 0.5*d;
}
Indices[NumEntries] = B;
Values[NumEntries++] = 0.25*f;
if(NumEntries != 3)
cout << "ERROR: calc corner (below) elements of stencil" << endl;
diag *= 0.25;
break;
}
case 2:{
if(W != -1 && E != -1) {
Indices[NumEntries] = W;
Values[NumEntries++] = 0.5*b;
Indices[NumEntries] = E;
Values[NumEntries++] = 0.5*b;
}
if(W == -1 && E != -1) {
Indices[NumEntries] = E;
Values[NumEntries++] = b;
}
if(E == -1 && W != -1) {
Indices[NumEntries] = W;
Values[NumEntries++] = b;
}
if(N != -1 && S != -1) {
Indices[NumEntries] = N;
Values[NumEntries++] = 0.5*d;
Indices[NumEntries] = S;
Values[NumEntries++] = 0.5*d;
}
if(N == -1 && S != -1) {
Indices[NumEntries] = S;
Values[NumEntries++] = d;
}
if(S == -1 && N != -1) {
Indices[NumEntries] = N;
Values[NumEntries++] = d;
}
Indices[NumEntries] = B;