// // Class ArbitraryDomain // Interface to iterative solver and boundary geometry // for space charge calculation // // Copyright (c) 2008, Yves Ineichen, ETH Zürich, // 2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland // 2016, Daniel Winklehner, Massachusetts Institute of Technology // 2017 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland // All rights reserved // // Implemented as part of the master thesis // "A Parallel Multigrid Solver for Beam Dynamics" // and the paper // "A fast parallel Poisson solver on irregular domains applied to beam dynamics simulations" // (https://doi.org/10.1016/j.jcp.2010.02.022) // // This file is part of OPAL. // // OPAL is free software: you can redistribute it and/or modify // it under the terms of the GNU General Public License as published by // the Free Software Foundation, either version 3 of the License, or // (at your option) any later version. // // You should have received a copy of the GNU General Public License // along with OPAL. If not, see . // //#define DEBUG_INTERSECT_RAY_BOUNDARY #include "Solvers/ArbitraryDomain.h" #include "Structure/BoundaryGeometry.h" #include #include #include #include #include "Utilities/OpalException.h" ArbitraryDomain::ArbitraryDomain( BoundaryGeometry * bgeom, Vector_t nr, Vector_t hr, std::string interpl) : IrregularDomain(nr, hr, interpl) { bgeom_m = bgeom; setRangeMin(bgeom->getmincoords()); setRangeMax(bgeom->getmaxcoords()); bool have_inside_pt = bgeom->getInsidePoint(globalInsideP0_m); if (have_inside_pt == false) { throw OpalException( "ArbitraryDomain::ArbitraryDomain()", "No point inside geometry found/set!"); } startId = 0; } ArbitraryDomain::~ArbitraryDomain() { //nothing so far } void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){ INFOMSG(level2 << "* Starting the Boundary Intersection Tests..." << endl); setHr(hr); int zGhostOffsetLeft = (localId[2].first()== 0) ? 0 : 1; int zGhostOffsetRight = (localId[2].last() == nr_m[2] - 1) ? 0 : 1; int yGhostOffsetLeft = (localId[1].first()== 0) ? 0 : 1; int yGhostOffsetRight = (localId[1].last() == nr_m[1] - 1) ? 0 : 1; int xGhostOffsetLeft = (localId[0].first()== 0) ? 0 : 1; int xGhostOffsetRight = (localId[0].last() == nr_m[0] - 1) ? 0 : 1; hasGeometryChanged_m = true; IntersectLoX.clear(); IntersectHiX.clear(); IntersectLoY.clear(); IntersectHiY.clear(); IntersectLoZ.clear(); IntersectHiZ.clear(); // Calculate intersection Vector_t P, dir, I; // Vector_t saveP, saveP_old; Vector_t P0 = globalInsideP0_m; // We cannot assume that the geometry is symmetric about the xy, xz, and yz planes! // In my spiral inflector simulation, this is not the case for z direction for // example (-0.13 to +0.025). -DW for (int idz = localId[2].first()-zGhostOffsetLeft; idz <= localId[2].last()+zGhostOffsetRight; idz++) { //saveP_old[2] = (idz - (nr_m[2]-1)/2.0)*hr[2]; P[2] = min_m[2] + (idz + 0.5) * hr[2]; for (int idy = localId[1].first()-yGhostOffsetLeft; idy <= localId[1].last()+yGhostOffsetRight; idy++) { //saveP_old[1] = (idy - (nr_m[1]-1)/2.0)*hr[1]; P[1] = min_m[1] + (idy + 0.5) * hr[1]; for (int idx = localId[0].first()-xGhostOffsetLeft; idx <= localId[0].last()+xGhostOffsetRight; idx++) { //saveP_old[0] = (idx - (nr_m[0]-1)/2.0)*hr[0]; P[0] = min_m[0] + (idx + 0.5) * hr[0]; // *gmsg << "Now working on point " << saveP << " (original was " << saveP_old << ")" << endl; //P = saveP; if (bgeom_m->fastIsInside(P0, P) % 2 == 0) { // Fill the map with true or false values for very fast isInside tests // during the rest of the fieldsolve. IsInsideMap[toCoordIdx(idx, idy, idz)] = true; // Replace the old reference point with the new point (which we know is // inside because we just tested for it. This makes the algorithm faster // because fastIsInside() creates a number of segments that depends on the // distance between P and P0. Using the previous P as the new P0 // assures the smallest possible distance in most cases. -DW P0 = P; std::tuple pos(idx, idy, idz); dir = Vector_t(0, 0, 1); if (bgeom_m->intersectRayBoundary(P, dir, I)) { IntersectHiZ.insert(std::pair< std::tuple, double >(pos, I[2])); } else { #ifdef DEBUG_INTERSECT_RAY_BOUNDARY *gmsg << "zdir=+1 " << dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl; #endif } if (bgeom_m->intersectRayBoundary(P, -dir, I)) { IntersectLoZ.insert(std::pair< std::tuple, double >(pos, I[2])); } else { #ifdef DEBUG_INTERSECT_RAY_BOUNDARY *gmsg << "zdir=-1 " << -dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl; #endif } dir = Vector_t(0, 1, 0); if (bgeom_m->intersectRayBoundary(P, dir, I)) { IntersectHiY.insert(std::pair< std::tuple, double >(pos, I[1])); } else { #ifdef DEBUG_INTERSECT_RAY_BOUNDARY *gmsg << "ydir=+1 " << dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl; #endif } if (bgeom_m->intersectRayBoundary(P, -dir, I)) { IntersectLoY.insert(std::pair< std::tuple, double >(pos, I[1])); } else { #ifdef DEBUG_INTERSECT_RAY_BOUNDARY *gmsg << "ydir=-1" << -dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl; #endif } dir = Vector_t(1, 0, 0); if (bgeom_m->intersectRayBoundary(P, dir, I)) { IntersectHiX.insert(std::pair< std::tuple, double >(pos, I[0])); } else { #ifdef DEBUG_INTERSECT_RAY_BOUNDARY *gmsg << "xdir=+1 " << dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl; #endif } if (bgeom_m->intersectRayBoundary(P, -dir, I)){ IntersectLoX.insert(std::pair< std::tuple, double >(pos, I[0])); } else { #ifdef DEBUG_INTERSECT_RAY_BOUNDARY *gmsg << "xdir=-1 " << -dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl; #endif } } else { IsInsideMap[toCoordIdx(idx, idy, idz)] = false; #ifdef DEBUG_INTERSECT_RAY_BOUNDARY *gmsg << "OUTSIDE" << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl; #endif } } } } INFOMSG(level2 << "* Finding number of ghost nodes to the left..." << endl); //number of ghost nodes to the left int numGhostNodesLeft = 0; if(localId[2].first() != 0) { for(int idx = 0; idx < nr_m[0]; idx++) { for(int idy = 0; idy < nr_m[1]; idy++) { if(isInside(idx, idy, localId[2].first() - zGhostOffsetLeft)) numGhostNodesLeft++; } } } INFOMSG(level2 << "* Finding number of xy points in each plane along z..." << endl); //xy points in z plane int numtotal = 0; numXY.clear(); for(int idz = localId[2].first(); idz <= localId[2].last(); idz++) { int numxy = 0; for(int idx = 0; idx < nr_m[0]; idx++) { for(int idy = 0; idy < nr_m[1]; idy++) { if(isInside(idx, idy, idz)) numxy++; } } numXY[idz-localId[2].first()] = numxy; numtotal += numxy; } int startIdx = 0; MPI_Scan(&numtotal, &startIdx, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); startIdx -= numtotal; // Build up index and coord map idxMap_m.clear(); coordMap_m.clear(); int index = startIdx - numGhostNodesLeft; INFOMSG(level2 << "* Building up index and coordinate map..." << endl); for(int idz = localId[2].first() - zGhostOffsetLeft; idz <= localId[2].last() + zGhostOffsetRight; idz++) { for(int idy = 0; idy < nr_m[1]; idy++) { for(int idx = 0; idx < nr_m[0]; idx++) { if(isInside(idx, idy, idz)) { idxMap_m[toCoordIdx(idx, idy, idz)] = index; coordMap_m[index] = toCoordIdx(idx, idy, idz); index++; } } } } INFOMSG(level2 << "* Done." << endl); } // Conversion from (x,y,z) to index in xyz plane inline int ArbitraryDomain::toCoordIdx(int idx, int idy, int idz) { return (idz * nr_m[1] + idy) * nr_m[0] + idx; } // Conversion from (x,y,z) to index on the 3D grid int ArbitraryDomain::getIdx(int idx, int idy, int idz) { if(isInside(idx, idy, idz) && idx >= 0 && idy >= 0 && idz >= 0) return idxMap_m[toCoordIdx(idx, idy, idz)]; else return -1; } // Conversion from a 3D index to (x,y,z) inline void ArbitraryDomain::getCoord(int idxyz, int &idx, int &idy, int &idz) { int id = coordMap_m[idxyz]; idx = id % (int)nr_m[0]; id /= nr_m[0]; idy = id % (int)nr_m[1]; id /= nr_m[1]; idz = id; } inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) { return IsInsideMap[toCoordIdx(idx, idy, idz)]; } /* inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) { Vector_t P; P[0] = min_m[0] + (idx + 0.5) * hr_m[0]; P[1] = min_m[1] + (idy + 0.5) * hr_m[1]; P[2] = min_m[2] + (idz + 0.5) * hr_m[2]; return (bgeom_m->fastIsInside(globalInsideP0_m, P) % 2 == 0); } */ /* inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) { Vector_t P; P[0] = (idx - (nr_m[0]-1)/2.0) * hr_m[0]; P[1] = (idy - (nr_m[1]-1)/2.0) * hr_m[1]; P[2] = (idz - (nr_m[2]-1)/2.0) * hr_m[2]; bool ret = false; int countH, countL; std::multimap < std::tuple, double >::iterator itrH, itrL; std::tuple coordxyz(idx, idy, idz); //check if z is inside with x,y coords itrH = IntersectHiZ.find(coordxyz); itrL = IntersectLoZ.find(coordxyz); countH = IntersectHiZ.count(coordxyz); countL = IntersectLoZ.count(coordxyz); if(countH == 1 && countL == 1) ret = (P[2] <= itrH->second) && (P[2] >= itrL->second); //check if y is inside with x,z coords itrH = IntersectHiY.find(coordxyz); itrL = IntersectLoY.find(coordxyz); countH = IntersectHiY.count(coordxyz); countL = IntersectLoY.count(coordxyz); if(countH == 1 && countL == 1) ret = ret && (P[1] <= itrH->second) && (P[1] >= itrL->second); //check if x is inside with y,z coords itrH = IntersectHiX.find(coordxyz); itrL = IntersectLoX.find(coordxyz); countH = IntersectHiX.count(coordxyz); countL = IntersectLoX.count(coordxyz); if(countH == 1 && countL == 1) ret = ret && (P[0] <= itrH->second) && (P[0] >= itrL->second); return ret; } */ int ArbitraryDomain::getNumXY(int z) { return numXY[z]; } void ArbitraryDomain::constantInterpolation(int idx, int idy, int idz, StencilValue_t& value, double& /*scaleFactor*/) { value.west = -1/(hr_m[0]*hr_m[0]); value.east = -1/(hr_m[0]*hr_m[0]); value.north = -1/(hr_m[1]*hr_m[1]); value.south = -1/(hr_m[1]*hr_m[1]); value.front = -1/(hr_m[2]*hr_m[2]); value.back = -1/(hr_m[2]*hr_m[2]); value.center = 2/(hr_m[0]*hr_m[0]) + 2/(hr_m[1]*hr_m[1]) + 2/(hr_m[2]*hr_m[2]); if(!isInside(idx-1,idy,idz)) value.west = 0.0; if(!isInside(idx+1,idy,idz)) value.east = 0.0; if(!isInside(idx,idy+1,idz)) value.north = 0.0; if(!isInside(idx,idy-1,idz)) value.south = 0.0; if(!isInside(idx,idy,idz-1)) value.front = 0.0; if(!isInside(idx,idy,idz+1)) value.back = 0.0; } void ArbitraryDomain::linearInterpolation(int idx, int idy, int idz, StencilValue_t& value, double &scaleFactor) { scaleFactor = 1; double cx = (idx - (nr_m[0]-1)/2.0)*hr_m[0]; double cy = (idy - (nr_m[1]-1)/2.0)*hr_m[1]; double cz = (idz - (nr_m[2]-1)/2.0)*hr_m[2]; double dx_w=hr_m[0]; double dx_e=hr_m[0]; double dy_n=hr_m[1]; double dy_s=hr_m[1]; double dz_f=hr_m[2]; double dz_b=hr_m[2]; value.center = 0.0; std::tuple coordxyz(idx, idy, idz); if (idx == nr_m[0]-1) dx_e = std::abs(IntersectHiX.find(coordxyz)->second - cx); if (idx == 0) dx_w = std::abs(IntersectLoX.find(coordxyz)->second - cx); if (idy == nr_m[1]-1) dy_n = std::abs(IntersectHiY.find(coordxyz)->second - cy); if (idy == 0) dy_s = std::abs(IntersectLoY.find(coordxyz)->second - cy); if (idz == nr_m[2]-1) dz_b = std::abs(IntersectHiZ.find(coordxyz)->second - cz); if (idz == 0) dz_f = std::abs(IntersectLoZ.find(coordxyz)->second - cz); if(dx_w != 0) value.west = -(dz_f + dz_b) * (dy_n + dy_s) / dx_w; else value.west = 0; if(dx_e != 0) value.east = -(dz_f + dz_b) * (dy_n + dy_s) / dx_e; else value.east = 0; if(dy_n != 0) value.north = -(dz_f + dz_b) * (dx_w + dx_e) / dy_n; else value.north = 0; if(dy_s != 0) value.south = -(dz_f + dz_b) * (dx_w + dx_e) / dy_s; else value.south = 0; if(dz_f != 0) value.front = -(dx_w + dx_e) * (dy_n + dy_s) / dz_f; else value.front = 0; if(dz_b != 0) value.back = -(dx_w + dx_e) * (dy_n + dy_s) / dz_b; else value.back = 0; //RHS scaleFactor for current 3D index //0.5* comes from discretiztaion //scaleFactor = 0.5*(dw+de)*(dn+ds)*(df+db); scaleFactor = 0.5; if(dx_w + dx_e != 0) scaleFactor *= (dx_w + dx_e); if(dy_n + dy_s != 0) scaleFactor *= (dy_n + dy_s); if(dz_f + dz_b != 0) scaleFactor *= (dz_f + dz_b); //catch the case where a point lies on the boundary double m1 = dx_w * dx_e; double m2 = dy_n * dy_s; if(dx_e == 0) m1 = dx_w; if(dx_w == 0) m1 = dx_e; if(dy_n == 0) m2 = dy_s; if(dy_s == 0) m2 = dy_n; value.center = 2 / hr_m[2]; if(dx_w != 0 || dx_e != 0) value.center *= (dx_w + dx_e); if(dy_n != 0 || dy_s != 0) value.center *= (dy_n + dy_s); if(dx_w != 0 || dx_e != 0) value.center += (dz_f + dz_b) * (dy_n + dy_s) * (dx_w + dx_e) / m1; if(dy_n != 0 || dy_s != 0) value.center += (dz_f + dz_b) * (dx_w + dx_e) * (dy_n + dy_s) / m2; } // vi: set et ts=4 sw=4 sts=4: // Local Variables: // mode:c // c-basic-offset: 4 // indent-tabs-mode: nil // require-final-newline: nil // End: