Commit d9f859dd authored by Tuelin Kaman's avatar Tuelin Kaman
Browse files

Space charge calculation for complicated 3D geometries

parent 42b385ee
......@@ -820,18 +820,13 @@ void PartBunch::computeSelfFields(int binNumber) {
}
void PartBunch::resizeMesh() {
double scaleFactor = 1.0;
//scaleFactor = (Physics::c * dt_m);
//get x, y range and scale to unit-less
double xmin = fs_m->solver_m->getXRangeMin() * scaleFactor ;
double xmax = fs_m->solver_m->getXRangeMax() * scaleFactor;
double ymin = fs_m->solver_m->getYRangeMin() * scaleFactor;
double ymax = fs_m->solver_m->getYRangeMax() * scaleFactor;
// Check if the new domain is larger than bunch max, mins
get_bounds(rmin_m, rmax_m);
//XXX: instead of assert delete oob particles!
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]) {
......@@ -851,21 +846,15 @@ 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);
hr_m[0] = (xmax - xmin) / (nr_m[0]-1);
hr_m[1] = (ymax - ymin) / (nr_m[1]-1);
//hr_m[2] = (rmax_m[2] - rmin_m[2]) / (nr_m[2] - 1);
for(int i = 0; i < 3; i++)
hr_m[i] = (mymax[i] - mymin[i]) / nr_m[i];
// we cannot increase the number of mesh points
// this would require to delete and recreate the
// particle bunch since the FieldLayout is fixed
// in ParticleBase
Vector_t mymin = Vector_t(xmin, ymin, rmin_m[2]);
// rescale mesh
getMesh().set_meshSpacing(&(hr_m[0]));
getMesh().set_origin(mymin);
getMesh().set_origin(mymin);
rho_m.initialize(getMesh(),
getFieldLayout(),
......@@ -877,6 +866,8 @@ void PartBunch::resizeMesh() {
vbc_m);
update();
// setGridIsFixed();
}
void PartBunch::computeSelfFields() {
......@@ -886,8 +877,8 @@ void PartBunch::computeSelfFields() {
if(fs_m->hasValidSolver()) {
if(fs_m->getFieldSolverType() == "MG") // || fs_m->getFieldSolverType() == "FFTBOX") {
resizeMesh();
// if(fs_m->getFieldSolverType() == "SAAMG" && !isGridFixed())
// resizeMesh();
INFOMSG("after resizeMesh" << hr_m << endl);
//scatter charges onto grid
......@@ -911,7 +902,6 @@ void PartBunch::computeSelfFields() {
//divide charge by a 'grid-cube' volume to get [C/m^3]
rho_m *= tmp2;
// #define DBG_SCALARFIELD
#ifdef DBG_SCALARFIELD
INFOMSG("*** START DUMPING SCALAR FIELD ***" << endl);
ofstream fstr1;
......
......@@ -2648,8 +2648,10 @@ void ParallelTTracker::initializeBoundaryGeometry() {
for(unsigned int i = 0; i < itsOpalBeamline_m.sections_m.size(); i++) {
bgf_m = itsOpalBeamline_m.getBoundaryGeometry(i);
if(!bgf_m) continue;
if(!bgf_m)
continue;
else
break;
Distribution *dist = NULL;
Distribution *distrand = NULL;
vector<string> distr_str = bgf_m->getDistributionArray();
......
This diff is collapsed.
/* can handle arbitrary domains
* at the moment implementation for validation against solution with zylinder on [0,1]
*/
#ifndef ARBITRARY_DOMAIN
#define ARBITRARY_DOMAIN
#ifndef ARBITRARY_DOMAIN_H
#define ARBITRARY_DOMAIN_H
#ifdef HAVE_SAAMG_SOLVER
#include "IrregularDomain.h"
#include <mpi.h>
#include <hdf5.h>
#include "H5hut.h"
#include <vector>
#include <map>
#include <string>
#include <math.h>
struct vertex {
h5_float64_t P[3];
};
typedef struct vertex vertex_t;
struct entity {
h5_id_t global_id;
h5_id_t parent_id;
h5_id_t vertexIDs[3];
};
typedef struct entity entity_t;
///PointList maps from an (x,z) resp. (y,z) pair to double values (=intersections with boundary)
// FIXME: replace MultiMap with Vector! (performance)
typedef std::multimap< std::pair<int, int>, double > PointList;
#include <cmath>
#include "IrregularDomain.h"
#include "Structure/BoundaryGeometry.h"
class ArbitraryDomain : public IrregularDomain {
public:
ArbitraryDomain(string fname, Vector_t nr_, Vector_t hr_, NDIndex<3> locidx);
ArbitraryDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr, std::string interpl);
~ArbitraryDomain();
/// load geometry file
void LoadFile();
/// calculates intersection with the elliptic beam pipe
void Compute(Vector_t hr);
/// returns number of nodes in xy plane
int getNumXY(int z);
/// 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);
/// 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 idxyz, 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, double &W, double &E, double &S, double &N, double &F, double &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, double &W, double &E, double &S, double &N, double &F, double &B);
void getNeighbours(int idxyz, int &W, int &E, int &S, int &N, int &F, int &B);
/// returns type of boundary condition
string getType() {return "Geometric";}
std::string getType() {return "Geometric";}
/// queries if a given (x,y,z) coordinate lies inside the domain
inline bool isInside(int x, int y, int z);
inline bool isInside(int x, int y, int z);
/// set a filename where the geometry residues
void setFilename(string fname) {filename = fname;}
/// returns number of nodes in xy plane
int getNumXY(int z);
/// calculates intersection
void Compute(Vector_t hr);
void Compute(Vector_t hr, NDIndex<3> localId);
int getStartIdx() {return startIdx;}
int getStartId() {return startId;}
/// conversion from (x,y,z) to index on the 3D grid
int getIdx(int x, int y, int z);
// int getIdx(int x, int y, int z);
//TODO
double getXRangeMin() { return 1.0; }
double getXRangeMax() { return 1.0; }
double getYRangeMin() { return 1.0; }
double getYRangeMax() { return 1.0; }
double getXRangeMin() { return Geo_mincoords_m(0); }
double getXRangeMax() { return Geo_maxcoords_m(0); }
double getYRangeMin() { return Geo_mincoords_m(1); }
double getYRangeMax() { return Geo_maxcoords_m(1); }
double getZRangeMin() { return Geo_mincoords_m(2); }
double getZRangeMax() { return Geo_maxcoords_m(2); }
bool hasGeometryChanged() { return hasGeometryChanged_m; }
private:
BoundaryGeometry *bgeom_m;
/// PointList maps from an (x,z) resp. (y,z) pair to double values (=intersections with boundary)
typedef std::multimap< std::tuple<int, int, int>, double > PointList;
/// all intersection points with gridlines in Y direction
PointList IntersectYDir;
/// all intersection points with gridlines in X direction
PointList IntersectXDir;
PointList IntersectHiX, IntersectLoX;
/// all intersection points with gridlines in Y direction
PointList IntersectHiY, IntersectLoY;
/// all intersection points with gridlines in Z direction
PointList IntersectZDir;
/// filepath to the mesh
string filename;
/// vector storing the triangles
std::map<h5_id_t, vertex_t> vertices;
std::map<h5_id_t, entity_t> entities;
PointList IntersectHiZ, IntersectLoZ;
int startIdx;
NDIndex<3> localidx;
int startId;
/// here we store the number of nodes in a xy layer for a given z coordinate
std::map<int, int> numXY;
std::map<int, int> numYZ;
std::map<int, int> numXZ;
/// mapping (x,y,z) -> idx
/// number of nodes in the xy plane (for this case: independent of the z coordinate)
int nxy_m;
/// mapping (x,y,z) -> idxyz
std::map<int, int> IdxMap;
/// mapping idx -> (x,y,z)
/// mapping idxyz -> (x,y,z)
std::map<int, int> CoordMap;
/// interpolation type
int interpolationMethod;
/// flag indicating if geometry has changed for the current time-step
bool hasGeometryChanged_m;
Vector_t Geo_nr_m;
Vector_t Geo_hr_m;
Vector_t Geo_mincoords_m;
Vector_t Geo_maxcoords_m;
/// conversion from (x,y,z) to index in xyz plane
inline int toCoordIdx(int x, int y, int z);
// conversion from (x,y,z) to index on the 3D grid
//inline int getIdx(int x, int y, int z);
/// conversion from (x,y,z) to index on the 3D grid
int getIdx(int x, int y, int z);
/// conversion from a 3D index to (x,y,z)
inline void getCoord(int idx, int &x, int &y, int &z);
inline void getCoord(int idxyz, int &x, int &y, int &z);
inline void crossProduct(double A[], double B[], double C[]);
inline double dotProduct(double v1[], double v2[]) { return (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]); }
/// 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);
};
#endif //#ifdef HAVE_SAAMG_SOLVER
......
......@@ -55,7 +55,7 @@ BoxCornerDomain::~BoxCornerDomain() {
// for the moment we center the box corner geometry around the center of the grid
// hr holds the grid-spacings (boundary ellipse embedded in hr-grid)
void BoxCornerDomain::Compute(Vector_t hr) {
void BoxCornerDomain::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]) {
......
......@@ -99,12 +99,15 @@ public:
/// set semi-major
//void setSemiMajor(double sm) {SemiMajor = sm;}
void Compute(Vector_t hr);
void Compute(Vector_t hr, NDIndex<3> localId);
double getXRangeMin() { return -A_m; }
double getXRangeMax() { return A_m; }
double getYRangeMin() { return -B_m;} // actBMin_m; }
double getYRangeMax() { return B_m; } // actBMax_m; }
double getZRangeMin() { return L1_m;}
double getZRangeMax() { return L1_m+L2_m; }
//TODO: ?
int getStartIdx() {return 0;}
......
......@@ -4,6 +4,7 @@ set (_SRCS
MGPoissonSolver.cpp
RectangularDomain.cpp
EllipticDomain.cpp
ArbitraryDomain.cpp
BoxCornerDomain.cpp
TaperDomain.cpp
)
......
#ifdef HAVE_SAAMG_SOLVER
#include <map>
#include <cmath>
#include <iostream>
......@@ -19,11 +20,11 @@ EllipticDomain::EllipticDomain(double semimajor, double semiminor, Vector_t nr,
setNr(nr);
setHr(hr);
if(interpl == "constant")
if(interpl == "CONSTANT")
interpolationMethod = CONSTANT;
else if(interpl == "linear")
else if(interpl == "LINEAR")
interpolationMethod = LINEAR;
else if(interpl == "quadratic")
else if(interpl == "QUADRATIC")
interpolationMethod = QUADRATIC;
}
......@@ -36,21 +37,17 @@ EllipticDomain::~EllipticDomain() {
// one x-y-plane
// for the moment we center the ellipse around the center of the grid
// hr holds the grid-spacings (boundary ellipse embedded in hr-grid)
void EllipticDomain::Compute(Vector_t hr) {
std::cout << " EllipticDomain::Compute" << hr[0] << " " << hr[1] << " " << hr[2] << std::endl;
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;
}
std::cout << " EllipticDomain::Compute" << getHr()[0] << " " << getHr()[1] << " " << getHr()[2] << std::endl;
setHr(hr);
hasGeometryChanged_m = true;
//reset number of points inside domain
nxy_m = 0;
// clear previous coordinate maps
IdxMap.clear();
CoordMap.clear();
......@@ -61,14 +58,12 @@ void EllipticDomain::Compute(Vector_t hr) {
// build a index and coordinate map
register int idx = 0;
register int x, y;
std::cout << " EllipticDomain::Compute" << nr[0] << " " << nr[1] << " " << nr[2] << std::endl;
for(x = 0; x < nr[0]; x++) {
for(y = 0; y < nr[1]; y++) {
if(isInside(x, y, 1)) {
//IdxMap[toCoordIdx(x, y)] = idx++;
std::cout << "x,y " << x << " "<< y << "-> " << toCoordIdx(x,y) << std::endl;
IdxMap[toCoordIdx(x, y)] = idx;
CoordMap[idx++] = toCoordIdx(x, y);
nxy_m++;
......@@ -78,7 +73,6 @@ void EllipticDomain::Compute(Vector_t hr) {
}
switch(interpolationMethod) {
case CONSTANT:
break;
case LINEAR:
......@@ -123,8 +117,9 @@ void EllipticDomain::Compute(Vector_t hr) {
}
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
// determine which interpolation method we use for points near the boundary
switch(interpolationMethod) {
case CONSTANT:
ConstantInterpolation(x, y, z, W, E, S, N, F, B, C, scaleFactor);
......
......@@ -39,12 +39,15 @@ public:
/// set semi-major
void setSemiMajor(double sm) {SemiMajor = sm;}
void Compute(Vector_t hr);
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 getZRangeMin() { return zMin_m; }
double getZRangeMax() { return zMax_m; }
//TODO: ?
int getStartIdx() {return 0;}
......
......@@ -58,6 +58,10 @@ public:
double getXRangeMax() {return a_m;}
double getYRangeMin() {return -a_m;}
double getYRangeMax() {return a_m;}
double getZRangeMin() {return -a_m; }
double getZRangeMax() {return a_m; }
Inform &print(Inform &os) const;
......
......@@ -67,6 +67,9 @@ public:
double getXRangeMax() {return 1.0;}
double getYRangeMin() {return 1.0;}
double getYRangeMax() {return 1.0;}
double getZRangeMin() {return 1.0;}
double getZRangeMax() {return 1.0;}
Inform &print(Inform &os) const;
private:
......
......@@ -5,6 +5,7 @@
#include <vector>
#include <string>
#include "Algorithms/PBunchDefs.h"
#define TK_PRINT(...) do{printf("TK:%s:%d:",__FILE__,__LINE__);printf(__VA_ARGS__);fflush(stdout);}while(0)
/// enumeration corresponding to different interpolation methods at the boundary
enum {
......@@ -21,7 +22,7 @@ public:
/** method to compute the intersection points with the boundary geometry (stored in some appropriate data structure)
* \param hr updated mesh spacings
*/
virtual void Compute(Vector_t hr) = 0;
virtual void Compute(Vector_t hr, NDIndex<3> localId) = 0;
/** method to get the number of gridpoints in a given z plane
* \param z coordinate of the z plane
......@@ -90,6 +91,8 @@ public:
virtual double getXRangeMax() = 0;
virtual double getYRangeMin() = 0;
virtual double getYRangeMax() = 0;
virtual double getZRangeMin() = 0;
virtual double getZRangeMax() = 0;
virtual int getIdx(int x, int y, int z) = 0;
virtual bool hasGeometryChanged() = 0;
......
......@@ -12,6 +12,7 @@
#include "AbstractObjects/OpalData.h"
#include "Utilities/Options.h"
#include <algorithm>
#include "omp.h"
using Physics::c;
......@@ -54,27 +55,27 @@ MGPoissonSolver::MGPoissonSolver(PartBunch &beam,Mesh_t *mesh, FieldLayout_t *fl
// Find CURRENT geometry
currentGeometry = geometries_m[0];
if(currentGeometry->getFilename() == "") {
if(currentGeometry->getTopology() == "ELLIPTIC")
if(currentGeometry->getTopology() == "ELLIPTIC"){
bp = new EllipticDomain(currentGeometry->getA(), currentGeometry->getB(), 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());
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(), localId);
} else {
ERRORMSG("Geometry not known" << endl);
exit(1);
}
} else {
cout << "hr_m = " << hr_m << endl;
localidx_m = layout_m->getLocalNDIndex();
//bp = new ArbitraryDomain(currentGeometry, currentGeometry->getmincoords(), currentGeometry->getmaxcoords(), orig_nr_m, hr_m, localidx_m, interpl);
localId = layout_m->getLocalNDIndex();
bp = new ArbitraryDomain(currentGeometry, orig_nr_m, hr_m, interpl);
}
#if 0
Map = 0;
A = Teuchos::null;
LHS = Teuchos::null;
RHS = Teuchos::null;
#endif
MLPrec = Teuchos::null;
numBlocks_m = Options::numBlocks;
recycleBlocks_m = Options::recycleBlocks;
......@@ -83,25 +84,22 @@ MGPoissonSolver::MGPoissonSolver(PartBunch &beam,Mesh_t *mesh, FieldLayout_t *fl
SetupBelosList();
// setup Belos solver
cout<<"numBlocks_m "<<numBlocks_m<<" " <<"recycleBlocks_m "<<recycleBlocks_m<<endl;
if(numBlocks_m == 0 || recycleBlocks_m == 0)
solver = rcp(new Belos::BlockCGSolMgr<double, MV, OP>());
else
solver = rcp(new Belos::RCGSolMgr<double, MV, OP>());
convStatusTest = rcp(new Belos::StatusTestGenResNorm<ST, MV, OP> (tol));
convStatusTest->defineScaleForm(Belos::NormOfRHS, Belos::TwoNorm);
#if defined(setUserConv)
solver->setUserConvStatusTest(convStatusTest);
#endif
//all timers used here
FunctionTimer1_m = IpplTimings::getTimer("Create Stencil");
FunctionTimer2_m = IpplTimings::getTimer("Prepare RHS");
FunctionTimer3_m = IpplTimings::getTimer("CG");
FunctionTimer4_m = IpplTimings::getTimer("LHS to Rho");
FunctionTimer5_m = IpplTimings::getTimer("recalculate Map");
FunctionTimer6_m = IpplTimings::getTimer("ML");
FunctionTimer1_m = IpplTimings::getTimer("BGF-IndexCoordMap");
FunctionTimer2_m = IpplTimings::getTimer("computeMap");
FunctionTimer3_m = IpplTimings::getTimer("IPPL to RHS");
FunctionTimer4_m = IpplTimings::getTimer("ComputeStencil");
FunctionTimer5_m = IpplTimings::getTimer("ML");
FunctionTimer6_m = IpplTimings::getTimer("Setup");
FunctionTimer7_m = IpplTimings::getTimer("CG");
FunctionTimer8_m = IpplTimings::getTimer("LHS to IPPL");
}
......@@ -119,50 +117,17 @@ void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr, double zshift)
throw OpalException("MGPoissonSolver", "not implemented yet");
}
void MGPoissonSolver::computeMap() {
IpplTimings::startTimer(FunctionTimer5_m);
NDIndex<3> localidx = layout_m->getLocalNDIndex();
void MGPoissonSolver::computeMap(NDIndex<3> localId) {
//FIXME: we have to set this to true after constructor has been called
//(otherwise Map 'size' is 0)
//if(hasParallelDecompositionChanged_m == true) {
if(true) {
if(hasParallelDecompositionChanged_m == true) {
if(Map != 0)
delete Map;
if(useRCB_m)
redistributeWithRCB(localidx);
redistributeWithRCB(localId);
else
IPPLToMap3D(localidx);
#if 0
//we also have to redistribute LHS
if(LHS == Teuchos::null) {
std::cout << "COMPUTEMAP LHS->PutScalar" <<std::endl;
LHS = rcp(new Epetra_Vector(*Map));
LHS->PutScalar(1.0);
} else {
std::cout << "COMPUTEMAP LHS->PutScalar NO" <<std::endl;
RCP<Epetra_Vector> tmplhs = rcp(new Epetra_Vector(*Map));
Epetra_Import importer(*Map, LHS->Map());
tmplhs->Import(*LHS, importer, Add);
LHS = tmplhs;
}
//...and all previously saved LHS
std::deque< Epetra_Vector >::iterator it = OldLHS.begin();
if(OldLHS.size() > 0) {
int n = OldLHS.size();
for(int i = 0; i < n; ++i) {
Epetra_Vector tmplhs = Epetra_Vector(*Map);
Epetra_Import importer(*Map, it->Map());
tmplhs.Import(*it, importer, Add);
*it = tmplhs;
++it;
}
}
#endif
IPPLToMap3D(localId);
}
IpplTimings::stopTimer(FunctionTimer5_m);
extrapolateLHS();
}
void MGPoissonSolver::extrapolateLHS() {
......@@ -171,9 +136,14 @@ void MGPoissonSolver::extrapolateLHS() {
// Pik (x) := (x − xi ) Pi+1,k−1(x) − (x − xi+k ) Pi,k−1(x) /(xi+k − xi )
// k = 1, . . . , n, i = 0, . . . , n − k.
if(LHS == Teuchos::null) {
LHS = rcp(new Epetra_Vector(*Map));
LHS->PutScalar(1.0);
}
std::deque< Epetra_Vector >::iterator it = OldLHS.begin();
if(nLHS == 0)
if(nLHS == 0 || OldLHS.size() == 0)
LHS->PutScalar(1.0);
else if(OldLHS.size() == 1)
*LHS = *it;
......@@ -195,7 +165,7 @@ void MGPoissonSolver::extrapolateLHS() {