ArbitraryDomain.cpp 12.8 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 "Utilities/OpalException.h"
gsell's avatar
gsell committed
38
#include "Index/NDIndex.h"
gsell's avatar
gsell committed
39

40
ArbitraryDomain::ArbitraryDomain( BoundaryGeometry * bgeom,
41
                                  IntVector_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
    }
57 58 59

    throw OpalException("ArbitraryDomain::ArbitraryDomain()",
                        "This domain is currently not available.");
gsell's avatar
gsell committed
60 61
}

62 63
ArbitraryDomain::~ArbitraryDomain() {
    //nothing so far
gsell's avatar
gsell committed
64 65
}

66
void ArbitraryDomain::compute(Vector_t hr, NDIndex<3> localId){
67

68 69
    INFOMSG(level2 << "* Starting the Boundary Intersection Tests..." << endl);

70
    setHr(hr);
71

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

79 80
    hasGeometryChanged_m = true;

81 82 83 84 85 86
    intersectLoX_m.clear();
    intersectHiX_m.clear();
    intersectLoY_m.clear();
    intersectHiY_m.clear();
    intersectLoZ_m.clear();
    intersectHiZ_m.clear();
87

88
    // Calculate intersection
89 90
    Vector_t P, dir, I;
    Vector_t P0 = globalInsideP0_m;
91

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

frey_m's avatar
frey_m committed
97
        P[2] = getZRangeMin() + (idz + 0.5) * hr[2];
98

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

frey_m's avatar
frey_m committed
101
            P[1] = getYRangeMin() + (idy + 0.5) * hr[1];
102

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

frey_m's avatar
frey_m committed
105
                P[0] = getXRangeMin() + (idx + 0.5) * hr[0];
106

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

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
109 110
                    // Fill the map with true or false values for very fast isInside tests
                    // during the rest of the fieldsolve.
frey_m's avatar
frey_m committed
111
                    isInsideMap_m[toCoordIdx(idx, idy, idz)] = true;
112

Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
113 114 115 116 117 118
                    // 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;
119

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

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

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

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

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

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

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

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

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

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

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

194 195 196
    //number of ghost nodes to the left
    int numGhostNodesLeft = 0;
    if(localId[2].first() != 0) {
frey_m's avatar
frey_m committed
197 198
        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
199
                if(isInside(idx, idy, localId[2].first() - zGhostOffsetLeft))
200 201 202 203 204
                    numGhostNodesLeft++;
            }
        }
    }

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

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

    int startIdx = 0;
223
    MPI_Scan(&numtotal, &startIdx, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
224
    startIdx -= numtotal;
225

226
    // Build up index and coord map
227 228
    idxMap_m.clear();
    coordMap_m.clear();
229
    int index = startIdx - numGhostNodesLeft;
230

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

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

    INFOMSG(level2 << "* Done." << endl);
246
}
247

248
void ArbitraryDomain::constantInterpolation(int idx, int idy, int idz,
249
                                            StencilValue_t& value, double& /*scaleFactor*/) const
frey_m's avatar
frey_m committed
250
{
frey_m's avatar
frey_m committed
251 252 253 254 255 256 257
    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]);
258

259
    if(!isInside(idx-1,idy,idz))
260
        value.west = 0.0;
kraus's avatar
kraus committed
261
    if(!isInside(idx+1,idy,idz))
262
        value.east = 0.0;
263

264
    if(!isInside(idx,idy+1,idz))
265
        value.north = 0.0;
kraus's avatar
kraus committed
266
    if(!isInside(idx,idy-1,idz))
267
        value.south = 0.0;
kraus's avatar
kraus committed
268 269

    if(!isInside(idx,idy,idz-1))
270
        value.front = 0.0;
kraus's avatar
kraus committed
271
    if(!isInside(idx,idy,idz+1))
272
        value.back = 0.0;
273
}
274

275
void ArbitraryDomain::linearInterpolation(int idx, int idy, int idz,
276
                                          StencilValue_t& value, double &scaleFactor) const
277
{
278
    scaleFactor = 1;
279

frey_m's avatar
frey_m committed
280 281 282
    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];
283

frey_m's avatar
frey_m committed
284 285 286 287 288 289
    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];
290
    value.center = 0.0;
291

292
    std::tuple<int, int, int> coordxyz(idx, idy, idz);
293

frey_m's avatar
frey_m committed
294
    if (idx == nr_m[0]-1)
295
        dx_e = std::abs(intersectHiX_m.find(coordxyz)->second - cx);
296
    if (idx == 0)
297
        dx_w = std::abs(intersectLoX_m.find(coordxyz)->second - cx);
frey_m's avatar
frey_m committed
298
    if (idy == nr_m[1]-1)
299
        dy_n = std::abs(intersectHiY_m.find(coordxyz)->second - cy);
300
    if (idy == 0)
301
        dy_s = std::abs(intersectLoY_m.find(coordxyz)->second - cy);
frey_m's avatar
frey_m committed
302
    if (idz == nr_m[2]-1)
303
        dz_b = std::abs(intersectHiZ_m.find(coordxyz)->second - cz);
304
    if (idz == 0)
305
        dz_f = std::abs(intersectLoZ_m.find(coordxyz)->second - cz);
306 307

    if(dx_w != 0)
308
        value.west = -(dz_f + dz_b) * (dy_n + dy_s) / dx_w;
309
    else
310
        value.west = 0;
311
    if(dx_e != 0)
312
        value.east = -(dz_f + dz_b) * (dy_n + dy_s) / dx_e;
313
    else
314
        value.east = 0;
315
    if(dy_n != 0)
316
        value.north = -(dz_f + dz_b) * (dx_w + dx_e) / dy_n;
317
    else
318
        value.north = 0;
319
    if(dy_s != 0)
320
        value.south = -(dz_f + dz_b) * (dx_w + dx_e) / dy_s;
321
    else
322
        value.south = 0;
323
    if(dz_f != 0)
324
        value.front = -(dx_w + dx_e) * (dy_n + dy_s) / dz_f;
325
    else
326
        value.front = 0;
327
    if(dz_b != 0)
328
        value.back = -(dx_w + dx_e) * (dy_n + dy_s) / dz_b;
329
    else
330
        value.back = 0;
331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354

    //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
355
    value.center = 2 / hr_m[2];
356
    if(dx_w != 0 || dx_e != 0)
357
        value.center *= (dx_w + dx_e);
358
    if(dy_n != 0 || dy_s != 0)
359
        value.center *= (dy_n + dy_s);
360
    if(dx_w != 0 || dx_e != 0)
361
        value.center += (dz_f + dz_b) * (dy_n + dy_s) * (dx_w + dx_e) / m1;
362
    if(dy_n != 0 || dy_s != 0)
363
        value.center += (dz_f + dz_b) * (dx_w + dx_e) * (dy_n + dy_s) / m2;
364
}