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

For ParallelT method: If BoundaryGeometry is allocated, use it for the...

For ParallelT method: If BoundaryGeometry is allocated, use it for the iterative solver in a field free region
parent 1cb64445
......@@ -878,9 +878,9 @@ void PartBunch::computeSelfFields() {
eg_m = Vector_t(0.0);
if(fs_m->hasValidSolver()) {
if (fs_m->getFieldSolverType() == "SAAMG")
resizeMesh();
//use the mesh that is already set
//if (fs_m->getFieldSolverType() == "SAAMG")
// resizeMesh();
INFOMSG("after resizeMesh" << hr_m << endl);
//scatter charges onto grid
......@@ -978,6 +978,20 @@ void PartBunch::computeSelfFields() {
//write out e field
#ifdef DBG_SCALARFIELD
// Immediate debug output:
// Output potential and e-field along the x-, y-, and z-axes
int m1 = (int)nr_m[0]-1;
int m2 = (int)nr_m[0]/2;
for (int i=0; i<m1; i++ )
*gmsg << "Field along x axis E = " << eg_m[i][m2][m2] << " Pot = " << rho_m[i][m2][m2] << endl;
for (int i=0; i<m1; i++ )
*gmsg << "Field along y axis E = " << eg_m[m2][i][m2] << " Pot = " << rho_m[m2][i][m2] << endl;
for (int i=0; i<m1; i++ )
*gmsg << "Field along z axis E = " << eg_m[m2][m2][i] << " Pot = " << rho_m[m2][m2][i] << endl;
INFOMSG("*** START DUMPING E FIELD ***" << endl);
//ostringstream oss;
//MPI_File file;
......
......@@ -252,8 +252,11 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
//calculate intersection
Vector_t P, saveP, dir, I;
Vector_t P0 = Vector_t(0,0,bgeom_m->getmincoords()[2]+hr[2]); //Reference Point inside the boundary
// TODO: Find and set the reference point for any time of geometry
//Reference Point inside the boundary for IsoDar Geo
Vector_t P0 = Vector_t(0,0,bgeom_m->getmincoords()[2]+hr[2]);
//Reference Point where the boundary geometry is cylinder
P0 = Vector_t(0,0,0);
for (int idz = localId[2].first()-zGhostOffsetLeft; idz <= localId[2].last()+zGhostOffsetRight; idz++) {
saveP[2] = (idz - (nr[2]-1)/2.0)*hr[2];
for (int idy = localId[1].first()-yGhostOffsetLeft; idy <= localId[1].last()+yGhostOffsetRight; idy++) {
......@@ -271,7 +274,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
rotateZAxisWithQuaternion(dir, localToGlobalQuaternion_m);
if (bgeom_m->intersectRayBoundary(P, dir, I)) {
setYRangeMax(MIN2(intersectMaxCoords_m(2),I[2]));
// setYRangeMax(MIN2(intersectMaxCoords_m(2),I[2]));
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectHiZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
......@@ -282,7 +285,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
}
if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
setYRangeMin(MAX2(intersectMinCoords_m(2),I[2]));
// setYRangeMin(MAX2(intersectMinCoords_m(2),I[2]));
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectLoZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
......@@ -294,7 +297,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
rotateYAxisWithQuaternion(dir, localToGlobalQuaternion_m);
if (bgeom_m->intersectRayBoundary(P, dir, I)) {
setZRangeMax(MIN2(intersectMaxCoords_m(1),I[1]));
// setZRangeMax(MIN2(intersectMaxCoords_m(1),I[1]));
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectHiY.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1]));
......@@ -305,7 +308,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
}
if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
setZRangeMin(MAX2(intersectMinCoords_m(1),I[1]));
// setZRangeMin(MAX2(intersectMinCoords_m(1),I[1]));
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectLoY.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1]));
......@@ -317,7 +320,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
rotateXAxisWithQuaternion(dir, localToGlobalQuaternion_m);
if (bgeom_m->intersectRayBoundary(P, dir, I)) {
setXRangeMax(MIN2(intersectMaxCoords_m(0),I[0]));
// setXRangeMax(MIN2(intersectMaxCoords_m(0),I[0]));
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectHiX.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[0]));
......@@ -328,7 +331,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
}
if (bgeom_m->intersectRayBoundary(P, -dir, I)){
setXRangeMin(MAX2(intersectMinCoords_m(0),I[0]));
// setXRangeMin(MAX2(intersectMinCoords_m(0),I[0]));
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectLoX.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[0]));
......
......@@ -245,7 +245,7 @@ void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) {
for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
if (bp->isInside(idx, idy, idz))
RHS->Values()[id++] = rho[idx][idy][idz].get()/scaleFactor;
RHS->Values()[id++] = 4*M_PI*rho[idx][idy][idz].get()/scaleFactor;
}
}
}
......
......@@ -262,6 +262,15 @@ void TrackRun::execute() {
*gmsg << " Selected Tracking Method is NOT implemented, good luck ..." << endl;
Beam *beam = Beam::find(Attributes::getString(itsAttr[BEAM]));
if (Attributes::getString(itsAttr[BOUNDARYGEOMETRY]) != "NONE") {
// Ask the dictionary if BoundaryGeometry is allocated.
// If it is allocated use the allocated BoundaryGeometry
if (!OpalData::getInstance()->hasGlobalGeometry()) {
BoundaryGeometry *bg = BoundaryGeometry::find(Attributes::getString(itsAttr[BOUNDARYGEOMETRY]))->
clone(getOpalName() + string("_geometry"));
OpalData::getInstance()->setGlobalGeometry(bg);
}
}
fs = FieldSolver::find(Attributes::getString(itsAttr[FIELDSOLVER]));
fs->initCartesianFields();
Track::block->bunch->setSolver(fs);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment