GreenWakeFunction.cpp 17 KB
Newer Older
gsell's avatar
gsell committed
1
#include "Solvers/GreenWakeFunction.hh"
2
#include "Algorithms/PartBunchBase.h"
kraus's avatar
kraus committed
3
#include "Utilities/GeneralClassicException.h"
4
#include "AbsBeamline/ElementBase.h"
5

gsell's avatar
gsell committed
6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
#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


/**
 * @brief   just a Testfunction!  Calculate the energy of the Wakefunction with the lambda
 *
 *
 * @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]   ref
 * @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,
41
                                     std::vector<Filter *> filters,
gsell's avatar
gsell committed
42 43 44 45 46 47 48 49
                                     int NBIN,
                                     double Z0,
                                     double radius,
                                     double sigma,
                                     int acMode,
                                     double tau,
                                     int direction,
                                     bool constLength,
50
                                     std::string fname):
51
    WakeFunction(name, NBIN),
gsell's avatar
gsell committed
52
    lineDensity_m(),
53
    //~ FftWField_m(0),
gsell's avatar
gsell committed
54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
    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() {
72 73 74
    //~ if(FftWField_m != 0) {
        //~ delete[] FftWField_m;
    //~ }
gsell's avatar
gsell committed
75 76 77 78 79 80 81 82 83 84 85 86 87
}


/**
 * @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
 */
88
std::pair<int, int> GreenWakeFunction::distrIndices(int vectLen) {
gsell's avatar
gsell committed
89

90
    std::pair<int, int> dist;
gsell's avatar
gsell committed
91 92 93 94 95 96 97 98 99 100 101 102

    //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
103
void GreenWakeFunction::apply(PartBunchBase<double, 3> *bunch) {
gsell's avatar
gsell committed
104 105

    Vector_t rmin, rmax;
frey_m's avatar
frey_m committed
106 107
    double charge = bunch->getChargePerParticle();
    // CKR: was bunch->Q[1] changed it;
gsell's avatar
gsell committed
108
    //FIXME: why 1? bunch,getTotalCharge()
frey_m's avatar
frey_m committed
109
    // or bunch->getChargePerParticle()?
gsell's avatar
gsell committed
110 111
    double K = 0; // constant to normalize the lineDensity_m to 1
    double spacing, mindist;
112
    std::vector<double> OutEnergy(NBin_m);
gsell's avatar
gsell committed
113

frey_m's avatar
frey_m committed
114 115
    bunch->calcBeamParameters();
    bunch->get_bounds(rmin, rmax);
gsell's avatar
gsell committed
116 117 118 119 120
    //FIXME IFF: do we have unitless r's here? is that what we want?

    mindist = rmin(2);
    switch(direction_m) {
        case LONGITUDINAL:
121
            spacing = std::abs(rmax(2) - rmin(2));
gsell's avatar
gsell committed
122 123 124 125 126
            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
127
            throw GeneralClassicException("GreenWakeFunction", "invalid direction specified");
gsell's avatar
gsell committed
128 129 130 131 132
    }
    assert(NBin_m > 0);
    spacing /= (NBin_m - 1); //FIXME: why -1? CKR: because grid spacings = grid points - 1

    // Calculate the Wakefield if needed
133 134
    if(FftWField_m.empty()) {
        FftWField_m.resize(2*NBin_m-1);
gsell's avatar
gsell committed
135 136 137 138 139 140 141 142 143 144
        if(filename_m != "") {
            setWakeFromFile(NBin_m, spacing);
        } else {
            CalcWakeFFT(spacing);
        }
    } else if(!constLength_m) {
        CalcWakeFFT(spacing);
    }

    // Calculate the line density of the particle bunch
145
    std::pair<double, double> meshInfo;
frey_m's avatar
frey_m committed
146
    bunch->calcLineDensity(nBins_m, lineDensity_m, meshInfo);
gsell's avatar
gsell committed
147 148 149 150 151 152 153 154

#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
155
    for(std::vector<Filter *>::const_iterator fit = filters_m.begin(); fit != filters_m.end(); ++fit) {
gsell's avatar
gsell committed
156 157 158 159 160 161 162 163 164
        (*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
165
    compEnergy(K, charge, lineDensity_m, OutEnergy.data());
gsell's avatar
gsell committed
166 167 168 169 170

    // 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
171
            for(unsigned int i = 0; i < bunch->getLocalNum(); i++) {
gsell's avatar
gsell committed
172 173 174

                //FIXME: Stimmt das????????? (von den einheiten)
                // calculate bin containing particle
frey_m's avatar
frey_m committed
175
                int idx = (int)(floor((bunch->R[i](2) - mindist) / spacing));
gsell's avatar
gsell committed
176 177 178 179
                //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
180
                bunch->Ef[i](2) += dE;
gsell's avatar
gsell committed
181 182 183 184 185

            }
            break;

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

                // calculate bin containing particle
frey_m's avatar
frey_m committed
189
                int idx = (int)(floor((bunch->R[i](2) - mindist) / spacing));
gsell's avatar
gsell committed
190 191 192 193 194 195
                //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
196
                double dist = sqrt(bunch->R[i](0) * bunch->R[i](0) + bunch->R[i](1) * bunch->R[i](1));
gsell's avatar
gsell committed
197
                assert(dist > 0);
frey_m's avatar
frey_m committed
198 199
                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
200 201 202 203 204

            }
            break;

        default:
kraus's avatar
kraus committed
205
            throw GeneralClassicException("GreenWakeFunction", "invalid direction specified");
gsell's avatar
gsell committed
206 207 208 209 210 211 212 213 214 215 216 217 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
    }

#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
243
    std::vector<double> pLambda(N);
gsell's avatar
gsell committed
244 245 246 247 248 249 250 251 252 253 254 255

    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
256
    gsl_fft_real_transform(pLambda.data(), 1, N, real, work);
gsell's avatar
gsell committed
257 258 259 260 261 262 263 264 265 266 267 268 269
    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);

270
    gsl_fft_halfcomplex_inverse(pLambda.data(), 1, N, hc, work);
gsell's avatar
gsell committed
271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298

    // 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,
299
                                   std::vector<double> lambda,
gsell's avatar
gsell committed
300 301 302
                                   double *OutEnergy) {
    int N = 2 * NBin_m - 1;
    // Allocate Space for the zero padded lambda and its Fourier Transformed
303
    std::vector<double> pLambda(N);
gsell's avatar
gsell committed
304 305 306 307 308 309 310 311 312 313 314 315

    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
316
    gsl_fft_real_transform(pLambda.data(), 1, N, real, work);
gsell's avatar
gsell committed
317 318 319 320 321 322 323 324 325 326 327 328 329 330
    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);

331
    gsl_fft_halfcomplex_inverse(pLambda.data(), 1, N, hc, work);
gsell's avatar
gsell committed
332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362

    // 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);

363
    const std::pair<int, int> myDist = distrIndices(NBin_m);
gsell's avatar
gsell committed
364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392
    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());

393
#ifdef ENABLE_WAKE_TESTS_FFT_OUT
394
    std::vector<double> wf(2*NBin_m-1);
gsell's avatar
gsell committed
395 396 397
    for(int i = 0; i < 2 * NBin_m - 1; ++ i) {
        wf[i] = FftWField_m[i];
    }
398
#endif
gsell's avatar
gsell committed
399
    // calculate the FFT of the Wakefield
400
    gsl_fft_real_transform(FftWField_m.data(), 1, M, real, work);
gsell's avatar
gsell committed
401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435


#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 ");
436
    std::string name;
gsell's avatar
gsell committed
437
    char temp[256];
snuverink_j's avatar
snuverink_j committed
438
    int Np;
gsell's avatar
gsell committed
439 440 441 442 443 444 445 446
    double dummy;
    gsl_fft_real_wavetable *real;
    gsl_fft_real_workspace *work;
    std::ifstream fs;

    fs.open(filename_m.c_str());

    if(fs.fail()) {
447
        throw GeneralClassicException("GreenWakeFunction::setWakeFromFile",
gsell's avatar
gsell committed
448 449 450 451 452 453 454
                            "Open file operation failed, please check if \""
                            + filename_m +  "\" really exists.");
    }

    fs >> name;
    msg << " SSDS1 read = " << name << endl;
    if(name.compare("SDDS1") != 0) {
455
        throw GeneralClassicException("GreenWakeFunction::setWakeFromFile",
gsell's avatar
gsell committed
456 457 458 459 460 461 462 463 464 465 466 467
                            " 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) {
468
        throw GeneralClassicException("GreenWakeFunction::setWakeFromFile",
gsell's avatar
gsell committed
469 470 471 472 473
                            " The particle number should be bigger than zero! Please check the first line of file \""
                            + filename_m +  "\".");
    }

    msg  << " Np = " << Np << endl;
474 475
    std::vector<double> wake(Np);
    std::vector<double> dist(Np);
476

gsell's avatar
gsell committed
477 478 479 480 481 482
    // read the wakefunction
    for(int i = 0; i < Np; i ++) {
        if(!fs.eof()) {
            fs >> dist[i] >> wake[i] >> dummy;
        }
        if(fs.eof()) {
483
            throw GeneralClassicException("GreenWakeFunction::setWakeFromFile",
gsell's avatar
gsell committed
484 485 486 487 488 489
                                " 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 (??)

490
    FftWField_m.resize(NBin_m);
gsell's avatar
gsell committed
491 492

    for(int i = 0; i < NBin_m; i ++) {
snuverink_j's avatar
snuverink_j committed
493
        int j = 0;
gsell's avatar
gsell committed
494 495 496 497 498 499 500 501 502 503 504
        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);

505
    gsl_fft_real_transform(FftWField_m.data(), 1, NBin_m, real, work);
gsell's avatar
gsell committed
506 507 508 509 510

    gsl_fft_real_wavetable_free(real);
    gsl_fft_real_workspace_free(work);
}

511
const std::string GreenWakeFunction::getType() const {
gsell's avatar
gsell committed
512
    return "GreenWakeFunction";
513
}