EllipticDomain.cpp 9.25 KB
Newer Older
1
//
frey_m's avatar
frey_m committed
2
// Class EllipticDomain
3 4 5
//   This class provides an elliptic beam pipe. The mesh adapts to the bunch size
//   in the longitudinal direction. At the intersection of the mesh with the
//   beam pipe, three stencil interpolation methods are available.
6
//
7
// Copyright (c) 2008,        Yves Ineichen, ETH Zürich,
frey_m's avatar
frey_m committed
8
//               2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland
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
//
kraus's avatar
kraus committed
28 29
#include "Solvers/EllipticDomain.h"

gsell's avatar
gsell committed
30 31 32 33 34 35 36
#include <map>
#include <cmath>
#include <iostream>

//FIXME: ORDER HOW TO TRAVERSE NODES IS FIXED, THIS SHOULD BE MORE GENERIC! (PLACES MARKED)


37
EllipticDomain::EllipticDomain(BoundaryGeometry *bgeom, IntVector_t nr, Vector_t hr,
38
                               std::string interpl)
39
    : RegularDomain(nr, hr, interpl)
40
{
41 42 43 44 45
    Vector_t min(-bgeom->getA(), -bgeom->getB(), bgeom->getS());
    Vector_t max( bgeom->getA(),  bgeom->getB(), bgeom->getS() + bgeom->getLength());
    setRangeMin(min);
    setRangeMax(max);
    setMinMaxZ(min[2], max[2]);
gsell's avatar
gsell committed
46 47 48 49 50 51 52 53 54
}

EllipticDomain::~EllipticDomain() {
    //nothing so far
}

// for this geometry we only have to calculate the intersection on
// one x-y-plane
// for the moment we center the ellipse around the center of the grid
55
void EllipticDomain::compute(Vector_t hr, NDIndex<3> localId) {
56
    // there is nothing to be done if the mesh spacings in x and y have not changed
57
    if (hr[0] == getHr()[0] &&
58
        hr[1] == getHr()[1])
59
    {
60 61 62
        hasGeometryChanged_m = false;
        return;
    }
63

64 65 66 67 68
    setHr(hr);
    hasGeometryChanged_m = true;
    //reset number of points inside domain

    // clear previous coordinate maps
69 70
    idxMap_m.clear();
    coordMap_m.clear();
71
    //clear previous intersection points
72 73
    intersectYDir_m.clear();
    intersectXDir_m.clear();
74 75

    // build a index and coordinate map
76 77
    int idx = 0;
    int x, y;
78

79 80
    int nxy = 0;

81 82 83 84 85
    /* we need to scan the full x-y-plane on all cores
     * in order to figure out the number of valid
     * grid points per plane --> otherwise we might
     * get not unique global indices in the Tpetra::CrsMatrix
     */
frey_m's avatar
frey_m committed
86 87
    for (x = 0; x < nr_m[0]; ++x) {
        for (y = 0; y < nr_m[1]; ++y) {
88 89 90
            if (isInside(x, y, 1)) {
                idxMap_m[toCoordIdx(x, y)] = idx;
                coordMap_m[idx++] = toCoordIdx(x, y);
91
                nxy++;
92 93 94 95
            }
        }
    }

96 97
    setNumXY(nxy);

98
    switch (interpolationMethod_m) {
99 100 101 102 103
        case CONSTANT:
            break;
        case LINEAR:
        case QUADRATIC:

104 105
            double smajsq = getXRangeMax() * getXRangeMax();
            double sminsq = getYRangeMax() * getYRangeMax();
106 107 108 109
            double yd = 0.0;
            double xd = 0.0;
            double pos = 0.0;

110
            // calculate intersection with the ellipse
frey_m's avatar
frey_m committed
111
            for (x = localId[0].first(); x <= localId[0].last(); x++) {
frey_m's avatar
frey_m committed
112
                pos = getXRangeMin() + hr_m[0] * (x + 0.5);
113
                if (std::abs(pos) >= getXRangeMax())
114
                {
115 116
                    intersectYDir_m.insert(std::pair<int, double>(x, 0));
                    intersectYDir_m.insert(std::pair<int, double>(x, 0));
frey_m's avatar
frey_m committed
117
                } else {
118
                    yd = std::abs(std::sqrt(sminsq - sminsq * pos * pos / smajsq));
119 120
                    intersectYDir_m.insert(std::pair<int, double>(x, yd));
                    intersectYDir_m.insert(std::pair<int, double>(x, -yd));
kraus's avatar
kraus committed
121
                }
122 123 124

            }

frey_m's avatar
frey_m committed
125
            for (y = localId[0].first(); y < localId[1].last(); y++) {
frey_m's avatar
frey_m committed
126
                pos = getYRangeMin() + hr_m[1] * (y + 0.5);
127
                if (std::abs(pos) >= getYRangeMax())
128
                {
129 130
                    intersectXDir_m.insert(std::pair<int, double>(y, 0));
                    intersectXDir_m.insert(std::pair<int, double>(y, 0));
frey_m's avatar
frey_m committed
131
                } else {
132
                    xd = std::abs(std::sqrt(smajsq - smajsq * pos * pos / sminsq));
133 134
                    intersectXDir_m.insert(std::pair<int, double>(y, xd));
                    intersectXDir_m.insert(std::pair<int, double>(y, -xd));
kraus's avatar
kraus committed
135
                }
136 137
            }
    }
138 139
}

140
void EllipticDomain::linearInterpolation(int x, int y, int z, StencilValue_t& value,
141
                                         double &scaleFactor) const
142
{
gsell's avatar
gsell committed
143 144
    scaleFactor = 1.0;

frey_m's avatar
frey_m committed
145 146
    double cx = getXRangeMin() + hr_m[0] * (x + 0.5);
    double cy = getYRangeMin() + hr_m[1] * (y + 0.5);
gsell's avatar
gsell committed
147 148

    double dx = 0.0;
149
    std::multimap<int, double>::const_iterator it = intersectXDir_m.find(y);
150 151

    if (cx < 0)
152
        ++it;
gsell's avatar
gsell committed
153 154 155
    dx = it->second;

    double dy = 0.0;
156 157
    it = intersectYDir_m.find(x);
    if (cy < 0)
158
        ++it;
gsell's avatar
gsell committed
159 160 161
    dy = it->second;


frey_m's avatar
frey_m committed
162 163 164 165
    double dw = hr_m[0];
    double de = hr_m[0];
    double dn = hr_m[1];
    double ds = hr_m[1];
166
    value.center = 0.0;
gsell's avatar
gsell committed
167

168
    // we are a right boundary point
169
    if (!isInside(x + 1, y, z)) {
170 171
        value.center += 1.0 / ((dx - cx) * de);
        value.east = 0.0;
gsell's avatar
gsell committed
172
    } else {
173 174
        value.center += 1.0 / (de * de);
        value.east = -1.0 / (de * de);
gsell's avatar
gsell committed
175 176
    }

177
    // we are a left boundary point
178
    if (!isInside(x - 1, y, z)) {
179 180
        value.center += 1.0 / ((std::abs(dx) - std::abs(cx)) * dw);
        value.west = 0.0;
gsell's avatar
gsell committed
181
    } else {
182 183
        value.center += 1.0 / (dw * dw);
        value.west = -1.0 / (dw * dw);
gsell's avatar
gsell committed
184 185
    }

186
    // we are a upper boundary point
187
    if (!isInside(x, y + 1, z)) {
188 189
        value.center += 1.0 / ((dy - cy) * dn);
        value.north = 0.0;
gsell's avatar
gsell committed
190
    } else {
191 192
        value.center += 1.0 / (dn * dn);
        value.north = -1.0 / (dn * dn);
gsell's avatar
gsell committed
193 194
    }

195
    // we are a lower boundary point
196
    if (!isInside(x, y - 1, z)) {
197 198
        value.center += 1.0 / ((std::abs(dy) - std::abs(cy)) * ds);
        value.south = 0.0;
gsell's avatar
gsell committed
199
    } else {
200 201
        value.center += 1.0 / (ds * ds);
        value.south = -1.0 / (ds * ds);
gsell's avatar
gsell committed
202 203 204
    }

    // handle boundary condition in z direction
frey_m's avatar
frey_m committed
205 206 207
    value.front = -1.0 / (hr_m[2] * hr_m[2]);
    value.back  = -1.0 / (hr_m[2] * hr_m[2]);
    value.center += 2.0 / (hr_m[2] * hr_m[2]);
208
    robinBoundaryStencil(z, value.front, value.back, value.center);
gsell's avatar
gsell committed
209 210
}

211
void EllipticDomain::quadraticInterpolation(int x, int y, int z,
212
                                            StencilValue_t& value,
213
                                            double &scaleFactor) const
214
{
215 216
    scaleFactor = 1.0;

frey_m's avatar
frey_m committed
217 218
    double cx = (x - (nr_m[0] - 1) / 2.0) * hr_m[0];
    double cy = (y - (nr_m[1] - 1) / 2.0) * hr_m[1];
gsell's avatar
gsell committed
219 220 221 222

    // since every vector for elliptic domains has ALWAYS size 2 we
    // can catch all cases manually
    double dx = 0.0;
223
    std::multimap<int, double>::const_iterator it = intersectXDir_m.find(y);
224
    if (cx < 0)
225
        ++it;
gsell's avatar
gsell committed
226 227 228
    dx = it->second;

    double dy = 0.0;
229 230
    it = intersectYDir_m.find(x);
    if (cy < 0)
231
        ++it;
gsell's avatar
gsell committed
232 233
    dy = it->second;

frey_m's avatar
frey_m committed
234 235 236 237
    double dw = hr_m[0];
    double de = hr_m[0];
    double dn = hr_m[1];
    double ds = hr_m[1];
238

239 240 241 242 243
    value.west = 0.0;
    value.east = 0.0;
    value.south = 0.0;
    value.north = 0.0;
    value.center = 0.0;
gsell's avatar
gsell committed
244

245
    // we are a right boundary point
246
    if (!isInside(x + 1, y, z)) {
247
        double s = dx - cx;
248 249 250
        value.center -= 2.0 * (s - 1.0) / (s * de);
        value.east = 0.0;
        value.west += (s - 1.0) / ((s + 1.0) * de);
251
    } else {
252 253
        value.center += 1.0 / (de * de);
        value.east = -1.0 / (de * de);
gsell's avatar
gsell committed
254 255
    }

256
    // we are a left boundary point
257
    if (!isInside(x - 1, y, z)) {
258
        double s = std::abs(dx) - std::abs(cx);
259 260 261
        value.center -= 2.0 * (s - 1.0) / (s * de);
        value.west = 0.0;
        value.east += (s - 1.0) / ((s + 1.0) * de);
262
    } else {
263 264
        value.center += 1.0 / (dw * dw);
        value.west = -1.0 / (dw * dw);
gsell's avatar
gsell committed
265 266
    }

267
    // we are a upper boundary point
268
    if (!isInside(x, y + 1, z)) {
269
        double s = dy - cy;
270 271 272
        value.center -= 2.0 * (s - 1.0) / (s * dn);
        value.north = 0.0;
        value.south += (s - 1.0) / ((s + 1.0) * dn);
273
    } else {
274 275
        value.center += 1.0 / (dn * dn);
        value.north = -1.0 / (dn * dn);
gsell's avatar
gsell committed
276 277
    }

278
    // we are a lower boundary point
279
    if (!isInside(x, y - 1, z)) {
280
        double s = std::abs(dy) - std::abs(cy);
281 282 283
        value.center -= 2.0 * (s - 1.0) / (s * dn);
        value.south = 0.0;
        value.north += (s - 1.0) / ((s + 1.0) * dn);
284
    } else {
285 286
        value.center += 1.0 / (ds * ds);
        value.south = -1.0 / (ds * ds);
gsell's avatar
gsell committed
287 288 289
    }

    // handle boundary condition in z direction
frey_m's avatar
frey_m committed
290 291 292
    value.front = -1.0 / (hr_m[2] * hr_m[2]);
    value.back  = -1.0 / (hr_m[2] * hr_m[2]);
    value.center += 2.0 / (hr_m[2] * hr_m[2]);
293
    robinBoundaryStencil(z, value.front, value.back, value.center);
294 295
}

296 297 298 299 300 301 302
// 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: