GreenWakeFunction.cpp 17.4 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
//
// Class GreenWakeFunction
//
// 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/GreenWakeFunction.hh"
18
#include "Algorithms/PartBunchBase.h"
kraus's avatar
kraus committed
19
#include "Utilities/GeneralClassicException.h"
20

gsell's avatar
gsell committed
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
#include <fstream>
#include <string>
#include <vector>
#include <istream>
#include <iostream>  // Needed for stream I/O
#include <iomanip>   // Needed for I/O manipulators
#include "gsl/gsl_fft_real.h"
#include "gsl/gsl_fft_halfcomplex.h"


//IFF: TEST
//#define ENABLE_WAKE_DEBUG
//#define ENABLE_WAKE_DUMP
//#define ENABLE_WAKE_TESTS_FFT_OUT


/**
 *
 * @todo        In this code one can only apply either the longitudinal wakefield or the transversal wakefield. One should implement that both wakefields can be applied to the particle beam
 * @todo        NBins must be set equal to MT of the fieldsolver. This should be changed. (the length of lineDensity_m must be NBins and not to MT)
 *
 * @param[in]   NBIN number of Bins
 * @param[in]   Z0 impedance of the tube
 * @param[in]   radius radius of the tube
 * @param[in]   sigma material constant
 * @param[in]   acMode 1 for AC and 2 for DC
 * @param[in]   tau material constant
 * @param[in]   direction 0 for transversal and 1 for Longitudinal
 * @param[in]   constLength  true if the length of the particle bunch is considered as constant
 * @param[in]   fname read wake from file
 */
GreenWakeFunction::GreenWakeFunction(const std::string &name,
53
                                     std::vector<Filter *> filters,
gsell's avatar
gsell committed
54 55 56 57 58 59 60 61
                                     int NBIN,
                                     double Z0,
                                     double radius,
                                     double sigma,
                                     int acMode,
                                     double tau,
                                     int direction,
                                     bool constLength,
62
                                     std::string fname):
63
    WakeFunction(name, NBIN),
gsell's avatar
gsell committed
64
    lineDensity_m(),
65
    //~ FftWField_m(0),
gsell's avatar
gsell committed
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83
    NBin_m(NBIN),
    Z0_m(Z0),
    radius_m(radius),
    sigma_m(sigma),
    acMode_m(acMode),
    tau_m(tau),
    direction_m(direction),
    constLength_m(constLength),
    filename_m(fname),
    filters_m(filters.begin(), filters.end()) {
#ifdef ENABLE_WAKE_DEBUG
    *gmsg << "* ************* W A K E ************************************************************ " << endl;
    *gmsg << "* Entered GreenWakeFunction::GreenWakeFunction " << '\n';
    *gmsg << "* ********************************************************************************** " << endl;
#endif
}

GreenWakeFunction::~GreenWakeFunction() {
84 85 86
    //~ if(FftWField_m != 0) {
        //~ delete[] FftWField_m;
    //~ }
gsell's avatar
gsell committed
87 88 89 90 91 92 93 94 95 96 97 98 99
}


/**
 * @brief   given a vector of length N, distribute the indexes among the available processors
 *
 *
 * @todo        make this function general available
 *
 *
 * @param[in]   length of vector
 * @param[out]  first: lowIndex, second: hiIndex
 */
100
std::pair<int, int> GreenWakeFunction::distrIndices(int vectLen) {
gsell's avatar
gsell committed
101

102
    std::pair<int, int> dist;
gsell's avatar
gsell committed
103 104 105 106 107 108 109 110 111 112 113 114

    //IFF: properly distribute remainder
    int rem = vectLen - (vectLen / Ippl::getNodes()) * Ippl::getNodes();
    int tmp = (rem > Ippl::myNode()) ? 1 : 0;
    int locBunchRange = vectLen / Ippl::getNodes() + tmp;

    dist.first = locBunchRange * Ippl::myNode() + (1 - tmp) * rem;
    dist.second = dist.first + locBunchRange - 1;

    return dist;
}

frey_m's avatar
frey_m committed
115
void GreenWakeFunction::apply(PartBunchBase<double, 3> *bunch) {
gsell's avatar
gsell committed
116 117

    Vector_t rmin, rmax;
frey_m's avatar
frey_m committed
118 119
    double charge = bunch->getChargePerParticle();
    // CKR: was bunch->Q[1] changed it;
gsell's avatar
gsell committed
120
    //FIXME: why 1? bunch,getTotalCharge()
frey_m's avatar
frey_m committed
121
    // or bunch->getChargePerParticle()?
gsell's avatar
gsell committed
122 123
    double K = 0; // constant to normalize the lineDensity_m to 1
    double spacing, mindist;
124
    std::vector<double> OutEnergy(NBin_m);
gsell's avatar
gsell committed
125

frey_m's avatar
frey_m committed
126 127
    bunch->calcBeamParameters();
    bunch->get_bounds(rmin, rmax);
gsell's avatar
gsell committed
128 129 130 131 132
    //FIXME IFF: do we have unitless r's here? is that what we want?

    mindist = rmin(2);
    switch(direction_m) {
        case LONGITUDINAL:
133
            spacing = std::abs(rmax(2) - rmin(2));
gsell's avatar
gsell committed
134 135 136 137 138
            break; //FIXME: Kann mann das Spacing immer ändern?
        case TRANSVERSAL:
            spacing = rmax(0) * rmax(0) + rmax(1) * rmax(1);
            break;
        default:
kraus's avatar
kraus committed
139
            throw GeneralClassicException("GreenWakeFunction", "invalid direction specified");
gsell's avatar
gsell committed
140 141 142 143 144
    }
    assert(NBin_m > 0);
    spacing /= (NBin_m - 1); //FIXME: why -1? CKR: because grid spacings = grid points - 1

    // Calculate the Wakefield if needed
145 146
    if(FftWField_m.empty()) {
        FftWField_m.resize(2*NBin_m-1);
gsell's avatar
gsell committed
147 148 149 150 151 152 153 154 155 156
        if(filename_m != "") {
            setWakeFromFile(NBin_m, spacing);
        } else {
            CalcWakeFFT(spacing);
        }
    } else if(!constLength_m) {
        CalcWakeFFT(spacing);
    }

    // Calculate the line density of the particle bunch
157
    std::pair<double, double> meshInfo;
frey_m's avatar
frey_m committed
158
    bunch->calcLineDensity(nBins_m, lineDensity_m, meshInfo);
gsell's avatar
gsell committed
159 160 161 162 163 164 165 166

#ifdef ENABLE_WAKE_DEBUG
    *gmsg << "* ************* W A K E ************************************************************ " << endl;
    *gmsg << "* GreenWakeFunction::apply  lineDensity_m.size() = " << lineDensity_m.size() << endl;
    *gmsg << "* ********************************************************************************** " << endl;
#endif

    // smooth the line density of the particle bunch
167
    for(std::vector<Filter *>::const_iterator fit = filters_m.begin(); fit != filters_m.end(); ++fit) {
gsell's avatar
gsell committed
168 169 170 171 172 173 174 175 176
        (*fit)->apply(lineDensity_m);
    }

    for(unsigned int i = 0; i < lineDensity_m.size(); i++) {
        K += lineDensity_m[i];
    }
    K = 1 / K;

    // compute the kick due to the wakefield
177
    compEnergy(K, charge, lineDensity_m, OutEnergy.data());
gsell's avatar
gsell committed
178 179 180 181 182

    // Add the right OutEnergy[i] to all the particles
    //FIXME: can we specify LONG AND TRANS?
    switch(direction_m) {
        case LONGITUDINAL:
frey_m's avatar
frey_m committed
183
            for(unsigned int i = 0; i < bunch->getLocalNum(); i++) {
gsell's avatar
gsell committed
184 185 186

                //FIXME: Stimmt das????????? (von den einheiten)
                // calculate bin containing particle
frey_m's avatar
frey_m committed
187
                int idx = (int)(floor((bunch->R[i](2) - mindist) / spacing));
gsell's avatar
gsell committed
188 189 190 191
                //IFF: should be ok
                if(idx == NBin_m) idx--;
                assert(idx >= 0 && idx < NBin_m);
                double dE = OutEnergy[idx];
frey_m's avatar
frey_m committed
192
                bunch->Ef[i](2) += dE;
gsell's avatar
gsell committed
193 194 195 196 197

            }
            break;

        case TRANSVERSAL:
frey_m's avatar
frey_m committed
198
            for(unsigned int i = 0; i < bunch->getLocalNum(); i++) {
gsell's avatar
gsell committed
199 200

                // calculate bin containing particle
frey_m's avatar
frey_m committed
201
                int idx = (int)(floor((bunch->R[i](2) - mindist) / spacing));
gsell's avatar
gsell committed
202 203 204 205 206 207
                //IFF: should be ok
                if(idx == NBin_m) idx--;
                assert(idx >= 0 && idx < NBin_m);
                double dE = OutEnergy[idx];

                // ACHTUNG spacing auch in transversal richtung
frey_m's avatar
frey_m committed
208
                double dist = sqrt(bunch->R[i](0) * bunch->R[i](0) + bunch->R[i](1) * bunch->R[i](1));
gsell's avatar
gsell committed
209
                assert(dist > 0);
frey_m's avatar
frey_m committed
210 211
                bunch->Ef[i](0) += dE * bunch->R[i](0) / dist;
                bunch->Ef[i](1) += dE * bunch->R[i](1) / dist;
gsell's avatar
gsell committed
212 213 214 215 216

            }
            break;

        default:
kraus's avatar
kraus committed
217
            throw GeneralClassicException("GreenWakeFunction", "invalid direction specified");
gsell's avatar
gsell committed
218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254
    }

#ifdef ENABLE_WAKE_DUMP
    ofstream  f2("OutEnergy.dat");
    f2 << "# Energy of the Wake calculated in Opal\n"
       << "# Z0 = " << Z0_m << "\n"
       << "# radius = " << radius << "\n"
       << "# sigma = " << sigma << "\n"
       << "# c = " << c << "\n"
       << "# acMode = " << acMode << "\n"
       << "# tau = " << tau << "\n"
       << "# direction = " << direction << "\n"
       << "# spacing = " << spacing << "\n"
       << "# Lbunch = " << NBin_m << "\n";
    for(int i = 0; i < NBin_m; i++) {
        f2 << i + 1 << " " << OutEnergy[i] << "\n";
    }
    f2.flush();
    f2.close();
#endif
}

/**
 * @brief   just a Testfunction!  Calculate the energy of the Wakefunction with the lambda
 *
 *
 * @param[in]   K a constant
 * @param[in]   charge a constant
 * @param[in]   lambda the distribution of the Particles
 * @param[out]  OutEnergy this is the Output
 */
void GreenWakeFunction::compEnergy(const double K,
                                   const double charge,
                                   const double *lambda,
                                   double *OutEnergy) {
    int N = 2 * NBin_m - 1;
    // Allocate Space for the zero padded lambda and its Fourier Transformed
255
    std::vector<double> pLambda(N);
gsell's avatar
gsell committed
256 257 258 259 260 261 262 263 264 265 266 267

    gsl_fft_halfcomplex_wavetable *hc;
    gsl_fft_real_wavetable *real = gsl_fft_real_wavetable_alloc(N);
    gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(N);

    // fill the arrays with data
    for(int i = 0; i < NBin_m; i ++) {
        pLambda[N-i-1] = 0.0;
        pLambda[i] = lambda[i];
    }

    //FFT of the lambda
268
    gsl_fft_real_transform(pLambda.data(), 1, N, real, work);
gsell's avatar
gsell committed
269 270 271 272 273 274 275 276 277 278 279 280 281
    gsl_fft_real_wavetable_free(real);

    // convolution -> multiplication in Fourier space
    pLambda[0] *= FftWField_m[0];
    for(int i = 1; i < N; i += 2) {
        double temp = pLambda[i];
        pLambda[i] = FftWField_m[i] * pLambda[i] - FftWField_m[i+1] * pLambda[i+1];
        pLambda[i+1] = FftWField_m[i] * pLambda[i+1] + FftWField_m[i+1] * temp;
    }

    // inverse transform to get c, the convolution of a and b;
    hc = gsl_fft_halfcomplex_wavetable_alloc(N);

282
    gsl_fft_halfcomplex_inverse(pLambda.data(), 1, N, hc, work);
gsell's avatar
gsell committed
283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310

    // Write the result to the output:
    for(int i = 0; i < NBin_m; i ++) {
        OutEnergy[i] = -charge * K * pLambda[i] / (2.0 * NBin_m) * N; // CKR: I doubt that the multiplication with N is correct,
        //      put it here to get the same result as with FFTW
        //      My suspicion: S. Pauli has forgotten that if you
        //      do an fft followed by an inverse fft you'll get
        //      N times your original data

    }


    gsl_fft_halfcomplex_wavetable_free(hc);
    gsl_fft_real_workspace_free(work);
}

/**
 * @brief   Calculate the energy of the Wakefunction with the lambda
 *
 *
 * @param[in]   K a constant
 * @param[in]   charge a constant
 * @param[in]   lambda the distribution of the Particles
 * @param[out]  OutEnergy this is the Output
 *
 */
void GreenWakeFunction::compEnergy(const double K,
                                   const double charge,
311
                                   std::vector<double> lambda,
gsell's avatar
gsell committed
312 313 314
                                   double *OutEnergy) {
    int N = 2 * NBin_m - 1;
    // Allocate Space for the zero padded lambda and its Fourier Transformed
315
    std::vector<double> pLambda(N);
gsell's avatar
gsell committed
316 317 318 319 320 321 322 323 324 325 326 327

    gsl_fft_halfcomplex_wavetable *hc;
    gsl_fft_real_wavetable *real = gsl_fft_real_wavetable_alloc(N);
    gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(N);

    // fill the arrays with data
    for(int i = 0; i < NBin_m; i ++) {
        pLambda[N-i-1] = 0.0;
        pLambda[i] = lambda[i];
    }

    //FFT of the lambda
328
    gsl_fft_real_transform(pLambda.data(), 1, N, real, work);
gsell's avatar
gsell committed
329 330 331 332 333 334 335 336 337 338 339 340 341 342
    gsl_fft_real_wavetable_free(real);


    // Convolution -> just a multiplication in Fourier space
    pLambda[0] *= FftWField_m[0];
    for(int i = 1; i < N; i += 2) {
        double temp = pLambda[i];
        pLambda[i] = FftWField_m[i] * pLambda[i] - FftWField_m[i+1] * pLambda[i+1];
        pLambda[i+1] = FftWField_m[i] * pLambda[i+1] + FftWField_m[i+1] * temp;
    }

    // IFFT
    hc = gsl_fft_halfcomplex_wavetable_alloc(N);

343
    gsl_fft_halfcomplex_inverse(pLambda.data(), 1, N, hc, work);
gsell's avatar
gsell committed
344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374

    // Write the result to the output:
    for(int i = 0; i < NBin_m; i ++) {
        OutEnergy[i] = -charge * K * pLambda[i] / (2.0 * NBin_m) * N; // CKR: I doubt that the multiplication with N is correct,
        //      put it here to get the same result as with FFTW
        //      My suspicion: S. Pauli has forgotten that if you
        //      do an fft followed by an inverse fft you'll get
        //      N times your original data

    }


    gsl_fft_halfcomplex_wavetable_free(hc);
    gsl_fft_real_workspace_free(work);
}

/**
 * @brief   Calculate the FFT of the Wakefunction
 *
 * @param[in]   spacing distance between 2 slice in the line distribution
 *
 */
void GreenWakeFunction::CalcWakeFFT(double spacing) {
    // Set integration properties
    double  a = 1, b = 1000000;
    unsigned int N = 1000000;
    int M = 2 * NBin_m - 1;

    gsl_fft_real_wavetable *real = gsl_fft_real_wavetable_alloc(M);
    gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(M);

375
    const std::pair<int, int> myDist = distrIndices(NBin_m);
gsell's avatar
gsell committed
376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404
    const int lowIndex = myDist.first;
    const int hiIndex  = myDist.second;

    for(int i = 0; i < M; i ++) {
        FftWField_m[i] = 0.0;
    }

    /**
      Calculate the Wakefield on all processors
      */

    //     if (Ippl::myNode() != Ippl::getNodes()-1) {
    for(int i = lowIndex; i <= hiIndex; i ++) {
        Wake w(i * spacing, Z0_m, radius_m, sigma_m, acMode_m, tau_m, direction_m);
        FftWField_m[i] = simpson(w, a, b, N);
    }
    //     } else {
    //         //IFF: changed  to <= with new distr
    //         for (int i = lowIndex; i <= hiIndex; i ++) {
    //             Wake w(i*spacing, Z0_m, radius_m, sigma_m, acMode_m, tau_m, direction_m);
    //             FftWField[i] = simpson(w,a,b,N);
    //         }
    //     }

    /**
      Reduce the results
      */
    reduce(&(FftWField_m[0]), &(FftWField_m[0]) + NBin_m, &(FftWField_m[0]), OpAddAssign());

405
#ifdef ENABLE_WAKE_TESTS_FFT_OUT
406
    std::vector<double> wf(2*NBin_m-1);
gsell's avatar
gsell committed
407 408 409
    for(int i = 0; i < 2 * NBin_m - 1; ++ i) {
        wf[i] = FftWField_m[i];
    }
410
#endif
gsell's avatar
gsell committed
411
    // calculate the FFT of the Wakefield
412
    gsl_fft_real_transform(FftWField_m.data(), 1, M, real, work);
gsell's avatar
gsell committed
413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447


#ifdef ENABLE_WAKE_TESTS_FFT_OUT
    ofstream  f2("FFTwake.dat");
    f2 << "# FFT of the Wake calculated in Opal" << "\n"
       << "# Z0 = " << Z0_m << "\n"
       << "# radius = " << radius_m << "\n"
       << "# sigma = " << sigma_m << "\n"
       << "# mode = " << acMode_m << "\n"
       << "# tau = " << tau_m << "\n"
       << "# direction = " << direction_m << "\n"
       << "# spacing = " << spacing << "\n"
       << "# Lbunch = " << NBin_m << "\n";

    f2 << "0\t" << FftWField_m[0] << "\t0.0\t" << wf[0] << "\n";
    for(int i = 1; i < M; i += 2) {
        f2 << (i + 1) / 2 << "\t"
           << FftWField_m[i] << "\t"
           << FftWField_m[i + 1] << "\t"
           << wf[(i+1)/2] << "\n";
    }


    f2.flush();
    f2.close();
#endif
    gsl_fft_real_wavetable_free(real);
    gsl_fft_real_workspace_free(work);
}

/**
 * @brief   reads in the wakefield from file
 */
void GreenWakeFunction::setWakeFromFile(int NBin_m, double spacing) {
    Inform msg("readSDDS ");
448
    std::string name;
gsell's avatar
gsell committed
449
    char temp[256];
snuverink_j's avatar
snuverink_j committed
450
    int Np;
gsell's avatar
gsell committed
451 452 453 454 455 456 457 458
    double dummy;
    gsl_fft_real_wavetable *real;
    gsl_fft_real_workspace *work;
    std::ifstream fs;

    fs.open(filename_m.c_str());

    if(fs.fail()) {
459
        throw GeneralClassicException("GreenWakeFunction::setWakeFromFile",
gsell's avatar
gsell committed
460 461 462 463 464 465 466
                            "Open file operation failed, please check if \""
                            + filename_m +  "\" really exists.");
    }

    fs >> name;
    msg << " SSDS1 read = " << name << endl;
    if(name.compare("SDDS1") != 0) {
467
        throw GeneralClassicException("GreenWakeFunction::setWakeFromFile",
gsell's avatar
gsell committed
468 469 470 471 472 473 474 475 476 477 478 479
                            " No SDDS1 File. A SDDS1 file should start with a SDDS1 String. Check file \""
                            + filename_m +  "\" ");
    }

    for(int i = 0; i < 6; i++) {
        fs.getline(temp, 256);
        msg << "line " << i << " :   " << temp << endl;
    }

    fs >> Np;
    msg << " header read" << endl;
    if(Np <= 0) {
480
        throw GeneralClassicException("GreenWakeFunction::setWakeFromFile",
gsell's avatar
gsell committed
481 482 483 484 485
                            " The particle number should be bigger than zero! Please check the first line of file \""
                            + filename_m +  "\".");
    }

    msg  << " Np = " << Np << endl;
486 487
    std::vector<double> wake(Np);
    std::vector<double> dist(Np);
488

gsell's avatar
gsell committed
489 490 491 492 493 494
    // read the wakefunction
    for(int i = 0; i < Np; i ++) {
        if(!fs.eof()) {
            fs >> dist[i] >> wake[i] >> dummy;
        }
        if(fs.eof()) {
495
            throw GeneralClassicException("GreenWakeFunction::setWakeFromFile",
gsell's avatar
gsell committed
496 497 498 499 500 501
                                " End of file reached before the whole wakefield is imported, please check file \""
                                + filename_m +  "\".");
        }
    }
    // if needed interpolate the wake in a way that the wake form the file fits to the wake needs in the code (??)

502
    FftWField_m.resize(NBin_m);
gsell's avatar
gsell committed
503 504

    for(int i = 0; i < NBin_m; i ++) {
snuverink_j's avatar
snuverink_j committed
505
        int j = 0;
gsell's avatar
gsell committed
506 507 508 509 510 511 512 513 514 515 516
        while(dist[j] < i * spacing) {
            j++;
        }
        -- j; // i * spacing should probably between dist[j] and dist[j+1]
        // linear interpolation
        FftWField_m[i] = wake[j] + ((wake[j+1] - wake[j]) / (dist[j+1] - dist[j]) * (i * spacing - dist[j]));
    }

    real = gsl_fft_real_wavetable_alloc(NBin_m);
    work = gsl_fft_real_workspace_alloc(NBin_m);

517
    gsl_fft_real_transform(FftWField_m.data(), 1, NBin_m, real, work);
gsell's avatar
gsell committed
518 519 520 521 522

    gsl_fft_real_wavetable_free(real);
    gsl_fft_real_workspace_free(work);
}

523
const std::string GreenWakeFunction::getType() const {
gsell's avatar
gsell committed
524
    return "GreenWakeFunction";
525
}