CSRWakeFunction.cpp 13 KB
Newer Older
gsell's avatar
gsell committed
1
#include "Solvers/CSRWakeFunction.hh"
2 3
#include "Solvers/RootFinderForCSR.h"

4
#include "Algorithms/PartBunchBase.h"
5
#include "Filters/Filter.h"
kraus's avatar
kraus committed
6
#include "Filters/SavitzkyGolay.h"
gsell's avatar
gsell committed
7 8 9
#include "Physics/Physics.h"
#include "AbsBeamline/RBend.h"
#include "AbsBeamline/SBend.h"
10
#include "Utilities/Options.h"
11
#include "Utilities/Util.h"
12

kraus's avatar
kraus committed
13 14
#include <iostream>
#include <fstream>
15
#include <cmath>
16

17

18 19
CSRWakeFunction::CSRWakeFunction(const std::string &name, std::vector<Filter *> filters, const unsigned int &N):
    WakeFunction(name, N),
gsell's avatar
gsell committed
20 21 22
    filters_m(filters.begin(), filters.end()),
    lineDensity_m(),
    dlineDensitydz_m(),
23
    bendRadius_m(0.0),
24
    totalBendAngle_m(0.0)
kraus's avatar
kraus committed
25 26 27 28 29 30 31 32
{
    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
33

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

frey_m's avatar
frey_m committed
37
    const double sPos = bunch->get_sPos();
38 39 40 41
    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
42

43
    if (Ez_m.size() < lineDensity_m.size()) {
44 45 46
        Ez_m.resize(lineDensity_m.size(), 0.0);
        Psi_m.resize(lineDensity_m.size(), 0.0);
    }
gsell's avatar
gsell committed
47 48

    Vector_t smin, smax;
frey_m's avatar
frey_m committed
49
    bunch->get_bounds(smin, smax);
50 51 52
    double minPathLength = smin(2) + sPos - FieldBegin_m;
    if (sPos + smax(2) < FieldBegin_m) return;

53 54
    Ez_m[0] = 0.0;
    // calculate wake field of bunch
55
    for (unsigned int i = 1; i < lineDensity_m.size(); ++i) {
56 57 58 59
        Ez_m[i] = 0.0;

        double angleOfSlice = 0.0;
        double pathLengthOfSlice = minPathLength + i * meshSpacing;
60
        if (pathLengthOfSlice > 0.0)
61 62 63 64 65 66 67 68 69
            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
70
    for (unsigned int i = 0; i < bunch->getLocalNum(); ++i) {
frey_m's avatar
frey_m committed
71
        const Vector_t &R = bunch->R[i];
72 73 74 75
        double distanceToOrigin = (R(2) - meshOrigin) / meshSpacing;

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

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

81
    if (Options::csrDump) {
82
        static std::string oldBendName;
83 84
        static unsigned long counter = 0;

85
        if (oldBendName != bendName_m) counter = 0;
86 87

        const int every = 1;
88
        bool print_criterion = (counter + 1) % every == 0;
89
        if (print_criterion) {
90
            static unsigned int file_number = 0;
91
            if (counter == 0) file_number = 0;
kraus's avatar
kraus committed
92
            if (Ippl::myNode() == 0) {
93
                std::stringstream filename_str;
kraus's avatar
kraus committed
94
                filename_str << "data/" << bendName_m << "-CSRWake" << std::setw(5) << std::setfill('0') << file_number << ".txt";
95 96 97 98

                std::ofstream csr(filename_str.str().c_str());
                csr << std::setprecision(8);
                csr << "# " << sPos + smin(2) - FieldBegin_m << "\t" << sPos + smax(2) - FieldBegin_m << std::endl;
99
                for (unsigned int i = 0; i < lineDensity_m.size(); ++ i) {
100 101 102 103 104 105 106
                  csr << i *meshSpacing << "\t"
                      << Ez_m[i] << "\t"
                      << lineDensity_m[i] << "\t"
                      << dlineDensitydz_m[i] << std::endl;
                }
                csr.close();
                msg << "** wrote " << filename_str.str() << endl;
kraus's avatar
kraus committed
107
            }
108
            ++ file_number;
109
        }
110 111
        ++ counter;
        oldBendName = bendName_m;
112 113
    }
}
gsell's avatar
gsell committed
114

115
void CSRWakeFunction::initialize(ElementBase::ConstSP const& ref) {
116
    if (ref->getType() == ElementBase::RBEND ||
117
       ref->getType() == ElementBase::SBEND) {
118

119
        std::shared_ptr<const Bend2D> bend = std::static_pointer_cast<const Bend2D>(ref);
120 121
        double End;

122
        bendRadius_m = bend->getBendRadius();
123
        bend->getDimensions(Begin_m, End);
124 125 126
        Length_m = bend->getEffectiveLength();
        FieldBegin_m = bend->getEffectiveCenter() - Length_m / 2.0;
        totalBendAngle_m = std::abs(bend->getBendAngle());
127
        bendName_m = bend->getName();
gsell's avatar
gsell committed
128
    }
129 130
}

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

134
    std::vector<Filter *>::const_iterator fit;
135
    for (fit = filters_m.begin(); fit != filters_m.end(); ++ fit) {
gsell's avatar
gsell committed
136 137
        (*fit)->apply(lineDensity_m);
    }
138

gsell's avatar
gsell committed
139
    dlineDensitydz_m.assign(lineDensity_m.begin(), lineDensity_m.end());
kraus's avatar
kraus committed
140
    diffOp_m->calc_derivative(dlineDensitydz_m, meshInfo.second);
141
}
gsell's avatar
gsell committed
142

143
void CSRWakeFunction::calculateContributionInside(size_t sliceNumber, double angleOfSlice, double meshSpacing) {
144
    if (angleOfSlice > totalBendAngle_m || angleOfSlice < 0.0) return;
145

146 147
    const double meshSpacingsup = std::pow(meshSpacing, -1. / 3.);
    double SlippageLength = std::pow(angleOfSlice, 3) * bendRadius_m / 24.;
148
    double relativeSlippageLength = SlippageLength / meshSpacing;
149
    if (relativeSlippageLength > sliceNumber) {
150 151

        /*
152 153
          Break integral into sum of integrals between grid points, then
          use linear interpolation between each grid point.
154 155
        */

156 157 158
        double dx1 = std::pow(sliceNumber, 2. / 3.);
        double dx2 = std::pow(sliceNumber, 5. / 3.);
        double dx3 = std::pow(sliceNumber - 1., 5. / 3.);
159
        Ez_m[sliceNumber] += 0.3 * meshSpacingsup * dlineDensitydz_m[0] * (5. * dx1 - 3. * dx2 + 3. * dx3);
160
        for (unsigned int j = 1; j < sliceNumber; ++ j) {
161 162
            dx1 = dx2;
            dx2 = dx3;
163
            dx3 = std::pow(sliceNumber - j - 1., 5. / 3.);
164 165 166
            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
167

168
    } else if (relativeSlippageLength < 1) {
gsell's avatar
gsell committed
169

170
        // First do transient term.
171
        if (4.0 * relativeSlippageLength <= 1) {
172

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

176
            if (4.0 * relativeSlippageLength < sliceNumber) {
177 178 179

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

            }
183

184
            Ez_m[sliceNumber] += (relativeSlippageLength * lineDensity_m[sliceNumber - 1] + (1. - relativeSlippageLength) * lineDensity_m[sliceNumber]) / std::pow(SlippageLength, 1. / 3.);
185 186 187 188

        }

        // Now do steady state term.
189
        Ez_m[sliceNumber] += (0.3 / meshSpacing) * std::pow(SlippageLength, 2. / 3.) * (5. * dlineDensitydz_m[sliceNumber] - 2. * relativeSlippageLength * (dlineDensitydz_m[sliceNumber] - dlineDensitydz_m[sliceNumber - 1]));
190 191 192

    } else {

193
        if (4. * relativeSlippageLength < sliceNumber) {
194 195 196

            int j = sliceNumber - static_cast<int>(floor(4. * relativeSlippageLength));
            double frac = 4. * relativeSlippageLength - (sliceNumber - j);
197
            Ez_m[sliceNumber] -= (frac * lineDensity_m[j - 1] + (1. - frac) * lineDensity_m[j]) / std::pow(SlippageLength, 1. / 3.);
198 199 200 201 202

        }

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

205 206 207 208
        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.);
209 210 211 212 213 214

        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;
215
        dx3 = std::pow(sliceNumber - j - 1., 5. / 3.);
216
        Ez_m[sliceNumber] += 0.3 * meshSpacingsup * dlineDensitydz_m[j] * (5.*dx1 - 3.*dx2 + 3.*dx3);
217
        for (unsigned int k = j + 1; k < sliceNumber; ++ k) {
218 219
            dx1 = dx2;
            dx2 = dx3;
220
            dx3 = std::pow(sliceNumber - k - 1., 5. / 3.);
221
            Ez_m[sliceNumber] += 0.9 * meshSpacingsup * dlineDensitydz_m[k] * (dx1 - 2.*dx2 + dx3);
gsell's avatar
gsell committed
222
        }
223
        Ez_m[sliceNumber] += 0.9 * meshSpacingsup * dlineDensitydz_m[sliceNumber];
gsell's avatar
gsell committed
224
    }
225
    double prefactor = -2. / std::pow(3. * bendRadius_m * bendRadius_m, 1. / 3.);
226 227
    Ez_m[sliceNumber] *= prefactor;
}
gsell's avatar
gsell committed
228

229
void CSRWakeFunction::calculateContributionAfter(size_t sliceNumber, double angleOfSlice, double meshSpacing) {
230
    if (angleOfSlice <= totalBendAngle_m) return;
231

232
    double Ds_max = bendRadius_m * std::pow(totalBendAngle_m, 3) / 24. * (4. - 3.* totalBendAngle_m / angleOfSlice);
233 234 235

    // First do contribution from particles whose retarded position is
    // prior to the bend.
236
    double Ds_max2 = bendRadius_m * std::pow(totalBendAngle_m, 2) / 6. * (3. * angleOfSlice - 2. * totalBendAngle_m);
237
    int j = 0;
238
    double frac = 0.0;
239
    if (Ds_max2 / meshSpacing < sliceNumber) {
240 241 242
        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
243 244
    }

245 246
    // Now do delta function contribution for particles whose retarded position
    // is in the bend.
247
    if (Ds_max / meshSpacing < sliceNumber) {
248 249 250 251
        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
252

253 254 255 256 257
    // Now do integral contribution for particles whose retarded position is in
    // the bend.

    double angleOverlap = angleOfSlice - totalBendAngle_m;
    int k = sliceNumber;
258
    if (Ds_max / meshSpacing < sliceNumber) {
259 260
        k = j;
        Psi_m[k] = calcPsi(Psi_m[k], angleOverlap, meshSpacing * (k + frac));
261
        if (Psi_m[k] > 0 && Psi_m[k] < totalBendAngle_m)
262 263 264
            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);
265
        if (Psi_m[0] > 0 && Psi_m[0] < totalBendAngle_m)
266
            Ez_m[sliceNumber] += 0.5 * dlineDensitydz_m[0] / (Psi_m[0] + 2.0 * angleOverlap);
gsell's avatar
gsell committed
267
    }
268 269

    // Do rest of integral.
270
    for (unsigned int l = sliceNumber - k + 1; l < sliceNumber; ++ l) {
271
        Psi_m[l] = calcPsi(Psi_m[l], angleOverlap, meshSpacing * (sliceNumber - l));
272
        if (Psi_m[l] > 0 && Psi_m[l] < totalBendAngle_m)
273
            Ez_m[sliceNumber] += dlineDensitydz_m[l] / (Psi_m[l] + 2.0 * angleOverlap);
gsell's avatar
gsell committed
274 275
    }

276 277 278
    // 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);
279
    if (Psi_m[sliceNumber] > 0 && Psi_m[sliceNumber] < totalBendAngle_m)
280 281 282 283
        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
284 285 286 287 288 289 290 291 292 293 294 295
}

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

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

300 301 302
    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
303
            return psi;
304

305
        psi -= residual / (4. * bendRadius_m * psi * psi * psi + 12. * x * bendRadius_m * psi * psi - 24. * Ds);
gsell's avatar
gsell committed
306
    }
307 308 309 310 311 312

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

313
    ERRORMSG("In CSRWakeFunction::calcPsi(): exceed maximum number of iterations!" << endl);
gsell's avatar
gsell committed
314 315 316
    return psi;
}

317
const std::string CSRWakeFunction::getType() const {
gsell's avatar
gsell committed
318
    return "CSRWakeFunction";
kraus's avatar
kraus committed
319
}