ArbitraryDomain.cpp 13.2 KB
Newer Older
1
//
2
// Class ArbitraryDomain
kraus's avatar
kraus committed
3
//   Interface to iterative solver and boundary geometry
4 5
//   for space charge calculation
//
6
// Copyright (c) 2008,        Yves Ineichen, ETH Zürich,
frey_m's avatar
frey_m committed
7 8
//               2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland
//                      2016, Daniel Winklehner, Massachusetts Institute of Technology
9
//               2017 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
frey_m's avatar
frey_m committed
10
// All rights reserved
11
//
12 13 14 15 16
// 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)
frey_m's avatar
frey_m committed
17 18 19 20 21 22 23 24 25 26
//
// 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 <https://www.gnu.org/licenses/>.
27 28
//

29
//#define DEBUG_INTERSECT_RAY_BOUNDARY
30

kraus's avatar
kraus committed
31
#include "Solvers/ArbitraryDomain.h"
snuverink_j's avatar
snuverink_j committed
32
#include "Structure/BoundaryGeometry.h"
kraus's avatar
kraus committed
33

34 35
#include <cmath>
#include <iostream>
snuverink_j's avatar
snuverink_j committed
36
#include <tuple>
37
#include <cassert>
38
#include "Utilities/OpalException.h"
gsell's avatar
gsell committed
39

40
ArbitraryDomain::ArbitraryDomain( BoundaryGeometry * bgeom,
kraus's avatar
kraus committed
41
                                  Vector_t nr,
frey_m's avatar
frey_m committed
42 43
                                  Vector_t hr,
                                  std::string interpl)
44
    : IrregularDomain(nr, hr, interpl)
frey_m's avatar
frey_m committed
45
{
46
    bgeom_m  = bgeom;
47 48 49 50

    setRangeMin(bgeom->getmincoords());
    setRangeMax(bgeom->getmaxcoords());

51 52
    bool have_inside_pt = bgeom->getInsidePoint(globalInsideP0_m);
    if (have_inside_pt == false) {
53 54 55
        throw OpalException(
            "ArbitraryDomain::ArbitraryDomain()",
            "No point inside geometry found/set!");
56
    }
gsell's avatar
gsell committed
57 58
}

59 60
ArbitraryDomain::~ArbitraryDomain() {
    //nothing so far
gsell's avatar
gsell committed
61 62
}

63
void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){
64

65 66
    INFOMSG(level2 << "* Starting the Boundary Intersection Tests..." << endl);

67
    setHr(hr);
68

69
    int zGhostOffsetLeft  = (localId[2].first()== 0) ? 0 : 1;
frey_m's avatar
frey_m committed
70
    int zGhostOffsetRight = (localId[2].last() == nr_m[2] - 1) ? 0 : 1;
71
    int yGhostOffsetLeft  = (localId[1].first()== 0) ? 0 : 1;
frey_m's avatar
frey_m committed
72
    int yGhostOffsetRight = (localId[1].last() == nr_m[1] - 1) ? 0 : 1;
73
    int xGhostOffsetLeft  = (localId[0].first()== 0) ? 0 : 1;
frey_m's avatar
frey_m committed
74
    int xGhostOffsetRight = (localId[0].last() == nr_m[0] - 1) ? 0 : 1;
kraus's avatar
kraus committed
75

76 77
    hasGeometryChanged_m = true;

78 79 80 81 82 83
    intersectLoX_m.clear();
    intersectHiX_m.clear();
    intersectLoY_m.clear();
    intersectHiY_m.clear();
    intersectLoZ_m.clear();
    intersectHiZ_m.clear();
84

85
    // Calculate intersection
86 87
    Vector_t P, dir, I;
    Vector_t P0 = globalInsideP0_m;
88

89
    // We cannot assume that the geometry is symmetric about the xy, xz, and yz planes!
90 91
    // In my spiral inflector simulation, this is not the case for z direction for
    // example (-0.13 to +0.025). -DW
92
    for (int idz = localId[2].first()-zGhostOffsetLeft; idz <= localId[2].last()+zGhostOffsetRight; idz++) {
93

94
        P[2] = min_m[2] + (idz + 0.5) * hr[2];
95

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
96
        for (int idy = localId[1].first()-yGhostOffsetLeft; idy <= localId[1].last()+yGhostOffsetRight; idy++) {
97

98
            P[1] = min_m[1] + (idy + 0.5) * hr[1];
99

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
100
            for (int idx = localId[0].first()-xGhostOffsetLeft; idx <= localId[0].last()+xGhostOffsetRight; idx++) {
101

102
                P[0] = min_m[0] + (idx + 0.5) * hr[0];
103

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
104
                if (bgeom_m->fastIsInside(P0, P) % 2 == 0) {
kraus's avatar
kraus committed
105

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
106 107 108
                    // 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;
109

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
110 111 112 113 114 115
                    // 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;
116

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
117
                    std::tuple<int, int, int> pos(idx, idy, idz);
118

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
119
                    dir = Vector_t(0, 0, 1);
120

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
121
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
122
                        intersectHiZ_m.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
123
                    } else {
124
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
125 126
                        *gmsg << "zdir=+1 " << dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
127
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
128
                    }
129

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
130
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
131
                        intersectLoZ_m.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
132
                    } else {
133
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
134 135
                        *gmsg << "zdir=-1 " << -dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
136
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
137
                    }
kraus's avatar
kraus committed
138

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
139
                    dir = Vector_t(0, 1, 0);
140

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
141
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
142
                        intersectHiY_m.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
143
                    } else {
144
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
145 146
                        *gmsg << "ydir=+1 " << dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
147
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
148
                    }
149

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
150
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
151
                        intersectLoY_m.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
152
                    } else {
153
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
154 155
                        *gmsg << "ydir=-1" << -dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
156
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
157
                    }
158

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
159
                    dir = Vector_t(1, 0, 0);
160

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
161
                    if (bgeom_m->intersectRayBoundary(P, dir, I)) {
162
                        intersectHiX_m.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[0]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
163
                    } else {
164
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
165 166
                        *gmsg << "xdir=+1 " << dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
167
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
168
                    }
169

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
170
                    if (bgeom_m->intersectRayBoundary(P, -dir, I)){
171
                        intersectLoX_m.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[0]));
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
172
                    } else {
173
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
174 175
                        *gmsg << "xdir=-1 " << -dir << " x,y,z= " << idx << "," << idy
                              << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
176
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
177 178 179
                    }
                } else {
                    IsInsideMap[toCoordIdx(idx, idy, idz)] = false;
180
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
frey_m's avatar
frey_m committed
181 182
                    *gmsg << "OUTSIDE" << " x,y,z= " << idx << "," << idy
                          << "," << idz << " P=" << P <<" I=" << I << endl;
kraus's avatar
kraus committed
183
#endif
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
184 185 186 187
                }
            }
        }
    }
188

189 190
    INFOMSG(level2 << "* Finding number of ghost nodes to the left..." << endl);

191 192 193
    //number of ghost nodes to the left
    int numGhostNodesLeft = 0;
    if(localId[2].first() != 0) {
frey_m's avatar
frey_m committed
194 195
        for(int idx = 0; idx < nr_m[0]; idx++) {
            for(int idy = 0; idy < nr_m[1]; idy++) {
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
196
                if(isInside(idx, idy, localId[2].first() - zGhostOffsetLeft))
197 198 199 200 201
                    numGhostNodesLeft++;
            }
        }
    }

202 203
    INFOMSG(level2 << "* Finding number of xy points in each plane along z..." << endl);

204 205 206 207
    //xy points in z plane
    int numtotal = 0;
    numXY.clear();
    for(int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
snuverink_j's avatar
snuverink_j committed
208
        int numxy = 0;
frey_m's avatar
frey_m committed
209 210
        for(int idx = 0; idx < nr_m[0]; idx++) {
            for(int idy = 0; idy < nr_m[1]; idy++) {
211 212 213 214 215 216 217 218 219
                if(isInside(idx, idy, idz))
                    numxy++;
            }
        }
        numXY[idz-localId[2].first()] = numxy;
        numtotal += numxy;
    }

    int startIdx = 0;
220
    MPI_Scan(&numtotal, &startIdx, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
221
    startIdx -= numtotal;
222

223
    // Build up index and coord map
224 225
    idxMap_m.clear();
    coordMap_m.clear();
226
    int index = startIdx - numGhostNodesLeft;
227

228 229
    INFOMSG(level2 << "* Building up index and coordinate map..." << endl);

230
    for(int idz = localId[2].first() - zGhostOffsetLeft; idz <= localId[2].last() + zGhostOffsetRight; idz++) {
frey_m's avatar
frey_m committed
231 232
        for(int idy = 0; idy < nr_m[1]; idy++) {
            for(int idx = 0; idx < nr_m[0]; idx++) {
kraus's avatar
kraus committed
233
                if(isInside(idx, idy, idz)) {
234 235
                    idxMap_m[toCoordIdx(idx, idy, idz)] = index;
                    coordMap_m[index] = toCoordIdx(idx, idy, idz);
236 237
                    index++;
                }
238
            }
239 240
        }
    }
241 242

    INFOMSG(level2 << "* Done." << endl);
243
}
244

245 246
// Conversion from (x,y,z) to index in xyz plane
inline int ArbitraryDomain::toCoordIdx(int idx, int idy, int idz) {
frey_m's avatar
frey_m committed
247
    return (idz * nr_m[1] + idy) * nr_m[0]  + idx;
gsell's avatar
gsell committed
248 249
}

250

frey_m's avatar
frey_m committed
251 252
int ArbitraryDomain::indexAccess(int x, int y, int z) {
    return idxMap_m[toCoordIdx(x, y, z)];
gsell's avatar
gsell committed
253 254
}

255
inline bool ArbitraryDomain::isInside(int idx, int idy, int idz) {
256

257 258 259
    return IsInsideMap[toCoordIdx(idx, idy, idz)];
}

gsell's avatar
gsell committed
260
int ArbitraryDomain::getNumXY(int z) {
kraus's avatar
kraus committed
261

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
262
    return numXY[z];
263
}
gsell's avatar
gsell committed
264

265 266
void ArbitraryDomain::constantInterpolation(int idx, int idy, int idz,
                                            StencilValue_t& value, double& /*scaleFactor*/)
frey_m's avatar
frey_m committed
267
{
frey_m's avatar
frey_m committed
268 269 270 271 272 273 274
    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]);
275

276
    if(!isInside(idx-1,idy,idz))
277
        value.west = 0.0;
kraus's avatar
kraus committed
278
    if(!isInside(idx+1,idy,idz))
279
        value.east = 0.0;
280

281
    if(!isInside(idx,idy+1,idz))
282
        value.north = 0.0;
kraus's avatar
kraus committed
283
    if(!isInside(idx,idy-1,idz))
284
        value.south = 0.0;
kraus's avatar
kraus committed
285 286

    if(!isInside(idx,idy,idz-1))
287
        value.front = 0.0;
kraus's avatar
kraus committed
288
    if(!isInside(idx,idy,idz+1))
289
        value.back = 0.0;
290
}
291

292 293
void ArbitraryDomain::linearInterpolation(int idx, int idy, int idz,
                                          StencilValue_t& value, double &scaleFactor)
294
{
295
    scaleFactor = 1;
296

frey_m's avatar
frey_m committed
297 298 299
    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];
300

frey_m's avatar
frey_m committed
301 302 303 304 305 306
    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];
307
    value.center = 0.0;
308

309
    std::tuple<int, int, int> coordxyz(idx, idy, idz);
310

frey_m's avatar
frey_m committed
311
    if (idx == nr_m[0]-1)
312
        dx_e = std::abs(intersectHiX_m.find(coordxyz)->second - cx);
313
    if (idx == 0)
314
        dx_w = std::abs(intersectLoX_m.find(coordxyz)->second - cx);
frey_m's avatar
frey_m committed
315
    if (idy == nr_m[1]-1)
316
        dy_n = std::abs(intersectHiY_m.find(coordxyz)->second - cy);
317
    if (idy == 0)
318
        dy_s = std::abs(intersectLoY_m.find(coordxyz)->second - cy);
frey_m's avatar
frey_m committed
319
    if (idz == nr_m[2]-1)
320
        dz_b = std::abs(intersectHiZ_m.find(coordxyz)->second - cz);
321
    if (idz == 0)
322
        dz_f = std::abs(intersectLoZ_m.find(coordxyz)->second - cz);
323 324

    if(dx_w != 0)
325
        value.west = -(dz_f + dz_b) * (dy_n + dy_s) / dx_w;
326
    else
327
        value.west = 0;
328
    if(dx_e != 0)
329
        value.east = -(dz_f + dz_b) * (dy_n + dy_s) / dx_e;
330
    else
331
        value.east = 0;
332
    if(dy_n != 0)
333
        value.north = -(dz_f + dz_b) * (dx_w + dx_e) / dy_n;
334
    else
335
        value.north = 0;
336
    if(dy_s != 0)
337
        value.south = -(dz_f + dz_b) * (dx_w + dx_e) / dy_s;
338
    else
339
        value.south = 0;
340
    if(dz_f != 0)
341
        value.front = -(dx_w + dx_e) * (dy_n + dy_s) / dz_f;
342
    else
343
        value.front = 0;
344
    if(dz_b != 0)
345
        value.back = -(dx_w + dx_e) * (dy_n + dy_s) / dz_b;
346
    else
347
        value.back = 0;
348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371

    //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;

frey_m's avatar
frey_m committed
372
    value.center = 2 / hr_m[2];
373
    if(dx_w != 0 || dx_e != 0)
374
        value.center *= (dx_w + dx_e);
375
    if(dy_n != 0 || dy_s != 0)
376
        value.center *= (dy_n + dy_s);
377
    if(dx_w != 0 || dx_e != 0)
378
        value.center += (dz_f + dz_b) * (dy_n + dy_s) * (dx_w + dx_e) / m1;
379
    if(dy_n != 0 || dy_s != 0)
380
        value.center += (dz_f + dz_b) * (dx_w + dx_e) * (dy_n + dy_s) / m2;
381
}
gsell's avatar
gsell committed
382

383 384 385 386 387 388 389
// 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: