CSRWakeFunction.cpp 13.7 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
//
// Class CSRWakeFunction
//
// Copyright (c) 2008 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// 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/>.
//
gsell's avatar
gsell committed
17
#include "Solvers/CSRWakeFunction.hh"
18 19
#include "Solvers/RootFinderForCSR.h"

20
#include "AbstractObjects/OpalData.h"
21
#include "Algorithms/PartBunchBase.h"
22
#include "Filters/Filter.h"
kraus's avatar
kraus committed
23
#include "Filters/SavitzkyGolay.h"
gsell's avatar
gsell committed
24 25 26
#include "Physics/Physics.h"
#include "AbsBeamline/RBend.h"
#include "AbsBeamline/SBend.h"
27
#include "Utilities/Options.h"
28
#include "Utilities/Util.h"
29

kraus's avatar
kraus committed
30 31
#include <iostream>
#include <fstream>
32
#include <cmath>
33

34

35 36
CSRWakeFunction::CSRWakeFunction(const std::string &name, std::vector<Filter *> filters, const unsigned int &N):
    WakeFunction(name, N),
gsell's avatar
gsell committed
37 38 39
    filters_m(filters.begin(), filters.end()),
    lineDensity_m(),
    dlineDensitydz_m(),
40
    bendRadius_m(0.0),
41
    totalBendAngle_m(0.0)
kraus's avatar
kraus committed
42 43 44 45 46 47 48 49
{
    if (filters_m.size() == 0) {
        defaultFilter_m.reset(new SavitzkyGolayFilter(7, 3, 3, 3));
        filters_m.push_back(defaultFilter_m.get());
    }

    diffOp_m = filters_m.back();
}
gsell's avatar
gsell committed
50

frey_m's avatar
frey_m committed
51
void CSRWakeFunction::apply(PartBunchBase<double, 3> *bunch) {
52
    Inform msg("CSRWake ");
gsell's avatar
gsell committed
53

frey_m's avatar
frey_m committed
54
    const double sPos = bunch->get_sPos();
55 56 57 58
    std::pair<double, double> meshInfo;
    calculateLineDensity(bunch, meshInfo);
    const double &meshSpacing = meshInfo.second;
    const double &meshOrigin = meshInfo.first + 0.5 * meshSpacing;
gsell's avatar
gsell committed
59

60
    if (Ez_m.size() < lineDensity_m.size()) {
61 62 63
        Ez_m.resize(lineDensity_m.size(), 0.0);
        Psi_m.resize(lineDensity_m.size(), 0.0);
    }
gsell's avatar
gsell committed
64 65

    Vector_t smin, smax;
frey_m's avatar
frey_m committed
66
    bunch->get_bounds(smin, smax);
67 68 69
    double minPathLength = smin(2) + sPos - FieldBegin_m;
    if (sPos + smax(2) < FieldBegin_m) return;

70 71
    Ez_m[0] = 0.0;
    // calculate wake field of bunch
72
    for (unsigned int i = 1; i < lineDensity_m.size(); ++i) {
73 74 75 76
        Ez_m[i] = 0.0;

        double angleOfSlice = 0.0;
        double pathLengthOfSlice = minPathLength + i * meshSpacing;
77
        if (pathLengthOfSlice > 0.0)
78 79 80 81 82 83 84 85 86
            angleOfSlice = pathLengthOfSlice / bendRadius_m;

        calculateContributionInside(i, angleOfSlice, meshSpacing);
        calculateContributionAfter(i, angleOfSlice, meshSpacing);

        Ez_m[i] /= (4. * Physics::pi * Physics::epsilon_0);
    }

    // calculate the wake field seen by the particles
87
    for (unsigned int i = 0; i < bunch->getLocalNum(); ++i) {
frey_m's avatar
frey_m committed
88
        const Vector_t &R = bunch->R[i];
89 90 91 92
        double distanceToOrigin = (R(2) - meshOrigin) / meshSpacing;

        unsigned int indexz = (unsigned int)floor(distanceToOrigin);
        double leverz = distanceToOrigin - indexz;
93
        PAssert_LT(indexz, lineDensity_m.size() - 1);
94

frey_m's avatar
frey_m committed
95
        bunch->Ef[i](2) += (1. - leverz) * Ez_m[indexz] + leverz * Ez_m[indexz + 1];
96 97
    }

98
    if (Options::csrDump) {
99
        static std::string oldBendName;
100 101
        static unsigned long counter = 0;

102
        if (oldBendName != bendName_m) counter = 0;
103 104

        const int every = 1;
105
        bool print_criterion = (counter + 1) % every == 0;
106
        if (print_criterion) {
107
            static unsigned int file_number = 0;
108
            if (counter == 0) file_number = 0;
kraus's avatar
kraus committed
109
            if (Ippl::myNode() == 0) {
110
                std::stringstream filename_str;
111 112 113 114 115
                filename_str << bendName_m << "-CSRWake" << std::setw(5) << std::setfill('0') << file_number << ".txt";
                std::string fname = Util::combineFilePath({
                    OpalData::getInstance()->getAuxiliaryOutputDirectory(),
                    filename_str.str()
                });
116

117
                std::ofstream csr(fname);
118 119
                csr << std::setprecision(8);
                csr << "# " << sPos + smin(2) - FieldBegin_m << "\t" << sPos + smax(2) - FieldBegin_m << std::endl;
120
                for (unsigned int i = 0; i < lineDensity_m.size(); ++ i) {
121 122 123 124 125 126
                  csr << i *meshSpacing << "\t"
                      << Ez_m[i] << "\t"
                      << lineDensity_m[i] << "\t"
                      << dlineDensitydz_m[i] << std::endl;
                }
                csr.close();
127
                msg << "** wrote " << fname << endl;
kraus's avatar
kraus committed
128
            }
129
            ++ file_number;
130
        }
131 132
        ++ counter;
        oldBendName = bendName_m;
133 134
    }
}
gsell's avatar
gsell committed
135

136
void CSRWakeFunction::initialize(const ElementBase *ref) {
137
    if (ref->getType() == ElementBase::RBEND ||
138
       ref->getType() == ElementBase::SBEND) {
139 140 141 142

        const Bend2D *bend = static_cast<const Bend2D *>(ref);
        double End;

143
        bendRadius_m = bend->getBendRadius();
144
        bend->getDimensions(Begin_m, End);
145 146 147
        Length_m = bend->getEffectiveLength();
        FieldBegin_m = bend->getEffectiveCenter() - Length_m / 2.0;
        totalBendAngle_m = std::abs(bend->getBendAngle());
148
        bendName_m = bend->getName();
gsell's avatar
gsell committed
149
    }
150 151
}

frey_m's avatar
frey_m committed
152 153
void CSRWakeFunction::calculateLineDensity(PartBunchBase<double, 3> *bunch, std::pair<double, double> &meshInfo) {
    bunch->calcLineDensity(nBins_m, lineDensity_m, meshInfo);
gsell's avatar
gsell committed
154

155
    std::vector<Filter *>::const_iterator fit;
156
    for (fit = filters_m.begin(); fit != filters_m.end(); ++ fit) {
gsell's avatar
gsell committed
157 158
        (*fit)->apply(lineDensity_m);
    }
159

gsell's avatar
gsell committed
160
    dlineDensitydz_m.assign(lineDensity_m.begin(), lineDensity_m.end());
kraus's avatar
kraus committed
161
    diffOp_m->calc_derivative(dlineDensitydz_m, meshInfo.second);
162
}
gsell's avatar
gsell committed
163

164
void CSRWakeFunction::calculateContributionInside(size_t sliceNumber, double angleOfSlice, double meshSpacing) {
165
    if (angleOfSlice > totalBendAngle_m || angleOfSlice < 0.0) return;
166

167 168
    const double meshSpacingsup = std::pow(meshSpacing, -1. / 3.);
    double SlippageLength = std::pow(angleOfSlice, 3) * bendRadius_m / 24.;
169
    double relativeSlippageLength = SlippageLength / meshSpacing;
170
    if (relativeSlippageLength > sliceNumber) {
171 172

        /*
173 174
          Break integral into sum of integrals between grid points, then
          use linear interpolation between each grid point.
175 176
        */

177 178 179
        double dx1 = std::pow(sliceNumber, 2. / 3.);
        double dx2 = std::pow(sliceNumber, 5. / 3.);
        double dx3 = std::pow(sliceNumber - 1., 5. / 3.);
180
        Ez_m[sliceNumber] += 0.3 * meshSpacingsup * dlineDensitydz_m[0] * (5. * dx1 - 3. * dx2 + 3. * dx3);
181
        for (unsigned int j = 1; j < sliceNumber; ++ j) {
182 183
            dx1 = dx2;
            dx2 = dx3;
184
            dx3 = std::pow(sliceNumber - j - 1., 5. / 3.);
185 186 187
            Ez_m[sliceNumber] += 0.9 * meshSpacingsup * dlineDensitydz_m[j] * (dx1 - 2.* dx2 + dx3);
        }
        Ez_m[sliceNumber] += 0.9 * meshSpacingsup * dlineDensitydz_m[sliceNumber];
gsell's avatar
gsell committed
188

189
    } else if (relativeSlippageLength < 1) {
gsell's avatar
gsell committed
190

191
        // First do transient term.
192
        if (4.0 * relativeSlippageLength <= 1) {
193

194
            Ez_m[sliceNumber] += 3.0 * std::pow(SlippageLength, 2.0 / 3.0) * (lineDensity_m[sliceNumber] - lineDensity_m[sliceNumber - 1]) / meshSpacing;
195 196
        } else {

197
            if (4.0 * relativeSlippageLength < sliceNumber) {
198 199 200

                int j = sliceNumber - static_cast<int>(floor(4.0 * relativeSlippageLength));
                double frac = 4.0 * relativeSlippageLength - (sliceNumber - j);
201
                Ez_m[sliceNumber] -= (frac * lineDensity_m[j - 1] + (1. - frac) * lineDensity_m[j]) / std::pow(SlippageLength, 1. / 3.);
gsell's avatar
gsell committed
202 203

            }
204

205
            Ez_m[sliceNumber] += (relativeSlippageLength * lineDensity_m[sliceNumber - 1] + (1. - relativeSlippageLength) * lineDensity_m[sliceNumber]) / std::pow(SlippageLength, 1. / 3.);
206 207 208 209

        }

        // Now do steady state term.
210
        Ez_m[sliceNumber] += (0.3 / meshSpacing) * std::pow(SlippageLength, 2. / 3.) * (5. * dlineDensitydz_m[sliceNumber] - 2. * relativeSlippageLength * (dlineDensitydz_m[sliceNumber] - dlineDensitydz_m[sliceNumber - 1]));
211 212 213

    } else {

214
        if (4. * relativeSlippageLength < sliceNumber) {
215 216 217

            int j = sliceNumber - static_cast<int>(floor(4. * relativeSlippageLength));
            double frac = 4. * relativeSlippageLength - (sliceNumber - j);
218
            Ez_m[sliceNumber] -= (frac * lineDensity_m[j - 1] + (1. - frac) * lineDensity_m[j]) / std::pow(SlippageLength, 1. / 3.);
219 220 221 222 223

        }

        int j = sliceNumber - static_cast<int>(floor(SlippageLength / meshSpacing));
        double frac = relativeSlippageLength - (sliceNumber - j);
224
        Ez_m[sliceNumber] += (frac * lineDensity_m[j - 1] + (1. - frac) * lineDensity_m[j]) / std::pow(SlippageLength, 1. / 3.);
225

226 227 228 229
        double dx1 = std::pow(sliceNumber - j + frac, 2. / 3.);
        double dx2 = std::pow(sliceNumber - j, 2. / 3.);
        double dx3 = std::pow(sliceNumber - j + frac, 5. / 3.);
        double dx4 = std::pow(sliceNumber - j, 5. / 3.);
230 231 232 233 234 235

        Ez_m[sliceNumber] += 1.5 * meshSpacingsup * dlineDensitydz_m[j - 1] * (dx1 - dx2);
        Ez_m[sliceNumber] += 0.3 * meshSpacingsup * (dlineDensitydz_m[j] - dlineDensitydz_m[j - 1]) * (5.*(dx1 - dx2) + 3.*(dx3 - dx4));

        dx1 = dx2;
        dx2 = dx4;
236
        dx3 = std::pow(sliceNumber - j - 1., 5. / 3.);
237
        Ez_m[sliceNumber] += 0.3 * meshSpacingsup * dlineDensitydz_m[j] * (5.*dx1 - 3.*dx2 + 3.*dx3);
238
        for (unsigned int k = j + 1; k < sliceNumber; ++ k) {
239 240
            dx1 = dx2;
            dx2 = dx3;
241
            dx3 = std::pow(sliceNumber - k - 1., 5. / 3.);
242
            Ez_m[sliceNumber] += 0.9 * meshSpacingsup * dlineDensitydz_m[k] * (dx1 - 2.*dx2 + dx3);
gsell's avatar
gsell committed
243
        }
244
        Ez_m[sliceNumber] += 0.9 * meshSpacingsup * dlineDensitydz_m[sliceNumber];
gsell's avatar
gsell committed
245
    }
246
    double prefactor = -2. / std::pow(3. * bendRadius_m * bendRadius_m, 1. / 3.);
247 248
    Ez_m[sliceNumber] *= prefactor;
}
gsell's avatar
gsell committed
249

250
void CSRWakeFunction::calculateContributionAfter(size_t sliceNumber, double angleOfSlice, double meshSpacing) {
251
    if (angleOfSlice <= totalBendAngle_m) return;
252

253
    double Ds_max = bendRadius_m * std::pow(totalBendAngle_m, 3) / 24. * (4. - 3.* totalBendAngle_m / angleOfSlice);
254 255 256

    // First do contribution from particles whose retarded position is
    // prior to the bend.
257
    double Ds_max2 = bendRadius_m * std::pow(totalBendAngle_m, 2) / 6. * (3. * angleOfSlice - 2. * totalBendAngle_m);
258
    int j = 0;
259
    double frac = 0.0;
260
    if (Ds_max2 / meshSpacing < sliceNumber) {
261 262 263
        j = sliceNumber - static_cast<int>(floor(Ds_max2 / meshSpacing));
        frac = Ds_max2 / meshSpacing - (sliceNumber - j);
        Ez_m[sliceNumber] -= (frac * lineDensity_m[j - 1] + (1. - frac) * lineDensity_m[j]) / (2. * angleOfSlice - totalBendAngle_m);
gsell's avatar
gsell committed
264 265
    }

266 267
    // Now do delta function contribution for particles whose retarded position
    // is in the bend.
268
    if (Ds_max / meshSpacing < sliceNumber) {
269 270 271 272
        j = sliceNumber - static_cast<int>(floor(Ds_max / meshSpacing));
        frac = Ds_max / meshSpacing - (sliceNumber - j);
        Ez_m[sliceNumber] += (frac * lineDensity_m[j - 1] + (1.0 - frac) * lineDensity_m[j]) / (2. * angleOfSlice - totalBendAngle_m);
    }
gsell's avatar
gsell committed
273

274 275 276 277 278
    // Now do integral contribution for particles whose retarded position is in
    // the bend.

    double angleOverlap = angleOfSlice - totalBendAngle_m;
    int k = sliceNumber;
279
    if (Ds_max / meshSpacing < sliceNumber) {
280 281
        k = j;
        Psi_m[k] = calcPsi(Psi_m[k], angleOverlap, meshSpacing * (k + frac));
282
        if (Psi_m[k] > 0 && Psi_m[k] < totalBendAngle_m)
283 284 285
            Ez_m[sliceNumber] += 0.5 * (frac * dlineDensitydz_m[sliceNumber - k - 1] + (1.0 - frac) * dlineDensitydz_m[sliceNumber - k]) / (Psi_m[k] + 2.0 * angleOverlap);
    } else {
        Psi_m[0] = calcPsi(Psi_m[0], angleOverlap, meshSpacing * sliceNumber);
286
        if (Psi_m[0] > 0 && Psi_m[0] < totalBendAngle_m)
287
            Ez_m[sliceNumber] += 0.5 * dlineDensitydz_m[0] / (Psi_m[0] + 2.0 * angleOverlap);
gsell's avatar
gsell committed
288
    }
289 290

    // Do rest of integral.
291
    for (unsigned int l = sliceNumber - k + 1; l < sliceNumber; ++ l) {
292
        Psi_m[l] = calcPsi(Psi_m[l], angleOverlap, meshSpacing * (sliceNumber - l));
293
        if (Psi_m[l] > 0 && Psi_m[l] < totalBendAngle_m)
294
            Ez_m[sliceNumber] += dlineDensitydz_m[l] / (Psi_m[l] + 2.0 * angleOverlap);
gsell's avatar
gsell committed
295 296
    }

297 298 299
    // We don't go right to the end as there is a singularity in the numerical integral that we don't quite know
    // how to deal with properly yet. This introduces a very slight error in the calculation (fractions of a percent).
    Psi_m[sliceNumber] = calcPsi(Psi_m[sliceNumber], angleOverlap, meshSpacing / 4.0);
300
    if (Psi_m[sliceNumber] > 0 && Psi_m[sliceNumber] < totalBendAngle_m)
301 302 303 304
        Ez_m[sliceNumber] += 0.5 * dlineDensitydz_m[sliceNumber] / (Psi_m[sliceNumber] + 2.0 * angleOverlap);

    double prefactor = -4 / bendRadius_m;
    Ez_m[sliceNumber] *= prefactor;
gsell's avatar
gsell committed
305 306 307 308 309 310 311 312 313 314 315 316
}

double CSRWakeFunction::calcPsi(const double &psiInitial, const double &x, const double &Ds) const {
    /** solve the equation
     *  \f[
     *  \Delta s = \frac{R \Psi^3}{24} \frac{\Psi + 4x}{\Psi + x}
     *  \f]
     *  for \f$\Psi\f$ using Newtons method.
     */

    const int Nmax = 100;
    const double eps = 1e-10;
317

318 319
    double psi = std::pow(24. * Ds / bendRadius_m, 1. / 3.);
    if (psiInitial != 0.0) psi = psiInitial;
gsell's avatar
gsell committed
320

321 322 323
    for (int i = 0; i < Nmax; ++i) {
    	double residual = bendRadius_m * psi * psi * psi * (psi + 4. * x) - 24. * Ds * psi - 24. * Ds * x;
        if (std::abs(residual) < eps)
gsell's avatar
gsell committed
324
            return psi;
325

326
        psi -= residual / (4. * bendRadius_m * psi * psi * psi + 12. * x * bendRadius_m * psi * psi - 24. * Ds);
gsell's avatar
gsell committed
327
    }
328 329 330 331 332 333

    RootFinderForCSR rootFinder(bendRadius_m, 4 * x * bendRadius_m, -24 * Ds, -24 * Ds * x);
    if (rootFinder.hasPositiveRealRoots()) {
        return rootFinder.searchRoot(eps);
    }

334
    ERRORMSG("In CSRWakeFunction::calcPsi(): exceed maximum number of iterations!" << endl);
gsell's avatar
gsell committed
335 336 337
    return psi;
}

338
const std::string CSRWakeFunction::getType() const {
gsell's avatar
gsell committed
339
    return "CSRWakeFunction";
kraus's avatar
kraus committed
340
}