Code indexing in gitaly is broken and leads to code not being visible to the user. We work on the issue with highest priority.

Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • OPAL/src
  • zheng_d/src
  • ext-rogers_c/src
  • ext-wang_c/src
  • cortes_c/src
  • ext-calvo_p/src
  • ext-edelen_a/src
  • albajacas_a/src
  • kraus/src
  • snuverink_j/OPAL-src
  • adelmann/src
  • muralikrishnan/src
  • wyssling_t/src
  • gsell/src
  • ext-piot_p/src
  • OPAL/opal-src-4-opalx-debug
  • winkle_m/src
17 results
Show changes
Commits on Source (3)
Hide whitespace changes
Inline Side-by-side
src/Solvers/ArbitraryDomain.cpp
View file @ bb04aec0
...@@ -115,74 +115,87 @@ void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){ ...@@ -115,74 +115,87 @@ void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){
// because fastIsInside() creates a number of segments that depends on the // because fastIsInside() creates a number of segments that depends on the
// distance between P and P0. Using the previous P as the new P0 // distance between P and P0. Using the previous P as the new P0
// assures the smallest possible distance in most cases. -DW // assures the smallest possible distance in most cases. -DW
P0 = P; //P0 = P;
std::tuple<int, int, int> pos(idx, idy, idz); std::tuple<int, int, int> pos(idx, idy, idz);
dir = Vector_t(0, 0, 1); dir = Vector_t(0, 0, 1);
if (bgeom_m->intersectRayBoundary(P, dir, I)) { if (bgeom_m->intersectRayBoundary(P, dir, I)) {
intersectHiZ_m.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2])); intersectHiZ_m.insert(
std::pair<std::tuple<int, int, int>, double>(pos, I[2]));
} else { } else {
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY #ifdef DEBUG_INTERSECT_RAY_BOUNDARY
*gmsg << "zdir=+1 " << dir << " x,y,z= " << idx << "," << idy INFOMSG(
<< "," << idz << " P=" << P <<" I=" << I << endl; level2 << "zdir=+1 " << dir << " x,y,z= " << idx << "," << idy << ","
<< idz << " P=" << P << " I=" << I << endl);
#endif #endif
} }
if (bgeom_m->intersectRayBoundary(P, -dir, I)) { if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
intersectLoZ_m.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2])); intersectLoZ_m.insert(
std::pair<std::tuple<int, int, int>, double>(pos, I[2]));
} else { } else {
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY #ifdef DEBUG_INTERSECT_RAY_BOUNDARY
*gmsg << "zdir=-1 " << -dir << " x,y,z= " << idx << "," << idy INFOMSG(
<< "," << idz << " P=" << P <<" I=" << I << endl; level2 << "zdir=-1 " << -dir << " x,y,z= " << idx << "," << idy << ","
<< idz << " P=" << P << " I=" << I << endl);
#endif #endif
} }
dir = Vector_t(0, 1, 0); dir = Vector_t(0, 1, 0);
if (bgeom_m->intersectRayBoundary(P, dir, I)) { if (bgeom_m->intersectRayBoundary(P, dir, I)) {
intersectHiY_m.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1])); intersectHiY_m.insert(
std::pair<std::tuple<int, int, int>, double>(pos, I[1]));
} else { } else {
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY #ifdef DEBUG_INTERSECT_RAY_BOUNDARY
*gmsg << "ydir=+1 " << dir << " x,y,z= " << idx << "," << idy INFOMSG(
<< "," << idz << " P=" << P <<" I=" << I << endl; level2 << "ydir=+1 " << dir << " x,y,z= " << idx << "," << idy << ","
<< idz << " P=" << P << " I=" << I << endl);
#endif #endif
} }
if (bgeom_m->intersectRayBoundary(P, -dir, I)) { if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
intersectLoY_m.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1])); intersectLoY_m.insert(
std::pair<std::tuple<int, int, int>, double>(pos, I[1]));
} else { } else {
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY #ifdef DEBUG_INTERSECT_RAY_BOUNDARY
*gmsg << "ydir=-1" << -dir << " x,y,z= " << idx << "," << idy INFOMSG(
<< "," << idz << " P=" << P <<" I=" << I << endl; level2 << "ydir=-1" << -dir << " x,y,z= " << idx << "," << idy << ","
<< idz << " P=" << P << " I=" << I << endl);
#endif #endif
} }
dir = Vector_t(1, 0, 0); dir = Vector_t(1, 0, 0);
if (bgeom_m->intersectRayBoundary(P, dir, I)) { if (bgeom_m->intersectRayBoundary(P, dir, I)) {
intersectHiX_m.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[0])); intersectHiX_m.insert(
std::pair<std::tuple<int, int, int>, double>(pos, I[0]));
} else { } else {
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY #ifdef DEBUG_INTERSECT_RAY_BOUNDARY
*gmsg << "xdir=+1 " << dir << " x,y,z= " << idx << "," << idy INFOMSG(
<< "," << idz << " P=" << P <<" I=" << I << endl; level2 << "xdir=+1 " << dir << " x,y,z= " << idx << "," << idy << ","
<< idz << " P=" << P << " I=" << I << endl);
#endif #endif
} }
if (bgeom_m->intersectRayBoundary(P, -dir, I)){ if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
intersectLoX_m.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[0])); intersectLoX_m.insert(
std::pair<std::tuple<int, int, int>, double>(pos, I[0]));
} else { } else {
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY #ifdef DEBUG_INTERSECT_RAY_BOUNDARY
*gmsg << "xdir=-1 " << -dir << " x,y,z= " << idx << "," << idy INFOMSG(
<< "," << idz << " P=" << P <<" I=" << I << endl; level2 << "xdir=-1 " << -dir << " x,y,z= " << idx << "," << idy << ","
<< idz << " P=" << P << " I=" << I << endl);
#endif #endif
} }
} else { } else {
isInsideMap_m[toCoordIdx(idx, idy, idz)] = false; isInsideMap_m[toCoordIdx(idx, idy, idz)] = false;
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY #ifdef DEBUG_INTERSECT_RAY_BOUNDARY
*gmsg << "OUTSIDE" << " x,y,z= " << idx << "," << idy INFOMSG(
<< "," << idz << " P=" << P <<" I=" << I << endl; level2 << "OUTSIDE" << " x,y,z= " << idx << "," << idy << "," << idz
<< " P=" << P << " I=" << I << endl);
#endif #endif
} }
} }
...@@ -256,19 +269,19 @@ void ArbitraryDomain::constantInterpolation(int idx, int idy, int idz, ...@@ -256,19 +269,19 @@ void ArbitraryDomain::constantInterpolation(int idx, int idy, int idz,
value.back = -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]); 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)) if (idx == 0 || !isInside(idx - 1, idy, idz))
value.west = 0.0; value.west = 0.0;
if(!isInside(idx+1,idy,idz)) if (idx == (nr_m[0] - 1) || !isInside(idx + 1, idy, idz))
value.east = 0.0; value.east = 0.0;
if(!isInside(idx,idy+1,idz)) if (idy == (nr_m[1] - 1) || !isInside(idx, idy + 1, idz))
value.north = 0.0; value.north = 0.0;
if(!isInside(idx,idy-1,idz)) if (idy == 0 || !isInside(idx, idy - 1, idz))
value.south = 0.0; value.south = 0.0;
if(!isInside(idx,idy,idz-1)) if (idz == 0 || !isInside(idx, idy, idz - 1))
value.front = 0.0; value.front = 0.0;
if(!isInside(idx,idy,idz+1)) if (idz == (nr_m[2] - 1) || !isInside(idx, idy, idz + 1))
value.back = 0.0; value.back = 0.0;
} }
...@@ -277,88 +290,93 @@ void ArbitraryDomain::linearInterpolation(int idx, int idy, int idz, ...@@ -277,88 +290,93 @@ void ArbitraryDomain::linearInterpolation(int idx, int idy, int idz,
{ {
scaleFactor = 1; scaleFactor = 1;
double cx = (idx - (nr_m[0]-1)/2.0)*hr_m[0]; double cx = getXRangeMin() + hr_m[0] * (idx + 0.5);
double cy = (idy - (nr_m[1]-1)/2.0)*hr_m[1]; double cy = getYRangeMin() + hr_m[1] * (idy + 0.5);
double cz = (idz - (nr_m[2]-1)/2.0)*hr_m[2]; double cz = getZRangeMin() + hr_m[2] * (idz + 0.5);
std::tuple<int, int, int> coordxyz(idx, idy, idz);
double dr = 0.0;
double dw = hr_m[0];
double de = hr_m[0];
double dn = hr_m[1];
double ds = hr_m[1];
double df = hr_m[2];
double db = 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; value.center = 0.0;
std::tuple<int, int, int> coordxyz(idx, idy, idz); // we are a left boundary point
if (idx == 0 || !isInside(idx - 1, idy, idz)) {
dr = std::abs(intersectLoX_m.find(coordxyz)->second - cx);
value.center += 1.0 / (dr * dw);
value.west = 0.0;
// ...we are not
} else {
value.center += 1.0 / (dw * dw);
value.west = -1.0 / (dw * dw);
}
// we are a right boundary point
if (idx == (nr_m[0] - 1) || !isInside(idx + 1, idy, idz)) {
dr = std::abs(intersectHiX_m.find(coordxyz)->second - cx);
value.center += 1.0 / (dr * de);
value.east = 0.0;
// ...we are not
} else {
value.center += 1.0 / (de * de);
value.east = -1.0 / (de * de);
}
// we are a lower boundary point
if (idy == 0 || !isInside(idx, idy - 1, idz)) {
dr = std::abs(intersectLoY_m.find(coordxyz)->second - cy);
value.center += 1.0 / (dr * ds);
value.south = 0.0;
// ...we are not
} else {
value.center += 1.0 / (ds * ds);
value.south = -1.0 / (ds * ds);
}
// we are a upper boundary point
if (idy == (nr_m[1] - 1) || !isInside(idx, idy + 1, idz)) {
dr = std::abs(intersectHiY_m.find(coordxyz)->second - cy);
value.center += 1.0 / (dr * dn);
value.north = 0.0;
// ...we are not
} else {
value.center += 1.0 / (dn * dn);
value.north = -1.0 / (dn * dn);
}
// We are a back boundary point
if (idz == 0 || !isInside(idx, idy, idz - 1)) {
dr = std::abs(intersectLoZ_m.find(coordxyz)->second - cz);
value.center += 1.0 / (dr * db);
value.back = 0.0;
// ...we are not
} else {
value.center += 1.0 / (db * db);
value.back = -1.0 / (db * db);
}
// We are a front boundary point
if (idz == (nr_m[2] - 1) || !isInside(idx, idy, idz + 1)) {
dr = std::abs(intersectHiZ_m.find(coordxyz)->second - cz);
value.center += 1.0 / (dr * df);
value.front = 0.0;
// ...we are not
} else {
value.center += 1.0 / (df * df);
value.front = -1.0 / (df * df);
}
if (idx == nr_m[0]-1) if (value.center <= 0.0) {
dx_e = std::abs(intersectHiX_m.find(coordxyz)->second - cx); INFOMSG(
if (idx == 0) level3 << "Got value.center <= 0.0 at idx/idy/idz = " << idx << "/" << idy << "/"
dx_w = std::abs(intersectLoX_m.find(coordxyz)->second - cx); << idz << endl);
if (idy == nr_m[1]-1) INFOMSG(level3 << "This will lead to an exception and quit OPAL." << endl);
dy_n = std::abs(intersectHiY_m.find(coordxyz)->second - cy); }
if (idy == 0)
dy_s = std::abs(intersectLoY_m.find(coordxyz)->second - cy);
if (idz == nr_m[2]-1)
dz_b = std::abs(intersectHiZ_m.find(coordxyz)->second - cz);
if (idz == 0)
dz_f = std::abs(intersectLoZ_m.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;
} }
src/Solvers/IrregularDomain.cpp
View file @ bb04aec0
...@@ -82,8 +82,8 @@ void IrregularDomain::getNeighbours(int id, StencilIndex_t& index) const { ...@@ -82,8 +82,8 @@ void IrregularDomain::getNeighbours(int id, StencilIndex_t& index) const {
void IrregularDomain::getCoord(int idx, int& x, int& y, int& z) const { void IrregularDomain::getCoord(int idx, int& x, int& y, int& z) const {
int xy = coordAccess(idx); int xy = coordAccess(idx);
x = xy % nr_m[0]; x = xy % nr_m[0];
y = (xy - x) / nr_m[0]; y = ((xy - x) / nr_m[0]) % nr_m[1];
z = idx / getNumXY(); z = ((xy - x) / nr_m[0] - y) / nr_m[1];
} }
int IrregularDomain::getIdx(int x, int y, int z) const { int IrregularDomain::getIdx(int x, int y, int z) const {
......
src/Solvers/MGPoissonSolver.cpp
View file @ bb04aec0
...@@ -354,9 +354,11 @@ void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) { ...@@ -354,9 +354,11 @@ void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) {
for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) { for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) { for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
NDIndex<3> l(Index(idx, idx), Index(idy, idy), Index(idz, idz)); NDIndex<3> l(Index(idx, idx), Index(idy, idy), Index(idz, idz));
if (bp_m->isInside(idx, idy, idz)) if (bp_m->isInside(idx, idy, idz)) {
id++;
RHS->replaceGlobalValue(bp_m->getIdx(idx, idy, idz), RHS->replaceGlobalValue(bp_m->getIdx(idx, idy, idz),
4.0 * M_PI * rho.localElement(l) / scaleFactor); 4.0 * M_PI * rho.localElement(l) / scaleFactor);
}
} }
} }
} }
......