p3m3dHeating.cpp 27.6 KB
Newer Older
frey_m's avatar
frey_m committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
//
// Application p3m3dHeating
//   mpirun -np 4 ./p3m3dHeating Nx Ny Nz l_beam l_box particleDensity r_cut alpha dt
//                               eps iterations charge_per_part m_part printEvery
//
//   alpha is the splitting parameter for pm and pp,
//   eps is the smoothing factor and Si are the coordinates of the charged sphere center
//
// Copyright (c) 2016, Benjamin Ulmer, ETH Zürich
// All rights reserved
//
// Implemented as part of the Master thesis
// "The P3M Model on Emerging Computer Architectures With Application to Microbunching"
// (http://amas.web.psi.ch/people/aadelmann/ETH-Accel-Lecture-1/projectscompleted/cse/thesisBUlmer.pdf)
//
// 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/>.
//
Andreas Adelmann's avatar
Andreas Adelmann committed
26 27 28 29 30 31 32 33 34 35 36 37 38
#include "Ippl.h"
#include <string>
#include <vector>
#include <iostream>
#include <cfloat>
#include <fstream>
#include <iomanip>
#include <complex>
#include "H5hut.h"
#include "Particle/BoxParticleCachingPolicy.h"
#include "Particle/PairBuilder/HashPairBuilderPeriodic.h"
#include "Particle/PairBuilder/HashPairBuilderPeriodicParallel.h"
#include "Particle/PairBuilder/PairConditions.h"
39
#include "Utility/PAssert.h"
Andreas Adelmann's avatar
Andreas Adelmann committed
40 41 42 43 44 45 46 47
#include "math.h"

#include <random>

#include "VTKFieldWriterParallel.hpp"
#include "ChargedParticleFactory.hpp"


Andreas Adelmann's avatar
Andreas Adelmann committed
48
// dimension of our positions
Andreas Adelmann's avatar
Andreas Adelmann committed
49 50 51 52 53 54 55 56 57 58 59 60 61
const unsigned Dim = 3;
//const double ke=1./(4.*M_PI*8.8e-14);
const double ke=2.532638e8;
// some typedefs
typedef UniformCartesian<Dim, double>                                 Mesh_t;
typedef BoxParticleCachingPolicy<double, Dim, Mesh_t>                 CachingPolicy_t;
typedef ParticleSpatialLayout<double, Dim, Mesh_t, CachingPolicy_t>   playout_t;
typedef playout_t::SingleParticlePos_t                                Vector_t;
typedef Cell                                                          Center_t;
typedef CenteredFieldLayout<Dim, Mesh_t, Center_t>                    FieldLayout_t;
typedef Field<double, Dim, Mesh_t, Center_t>                          Field_t;
typedef Field<int, Dim, Mesh_t, Center_t>                             IField_t;
typedef Field<Vector_t, Dim, Mesh_t, Center_t>                        VField_t;
gsell's avatar
gsell committed
62
typedef Field<std::complex<double>, Dim, Mesh_t, Center_t>                        CxField_t;
Andreas Adelmann's avatar
Andreas Adelmann committed
63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81
typedef FFT<CCTransform, Dim, double>                                 FFT_t;

typedef IntCIC                                                        IntrplCIC_t;
typedef IntNGP                                                        IntrplNGP_t;
typedef IntTSC                                                        IntrplTSC_t;

typedef UniformCartesian<2, double>                                   Mesh2d_t;
typedef CenteredFieldLayout<2, Mesh2d_t, Center_t>                    FieldLayout2d_t;
typedef Field<double, 2, Mesh2d_t, Center_t>                          Field2d_t;

template<class T>
struct ApplyField;

//This is the periodic Greens function with regularization parameter epsilon.
template<unsigned int Dim>
struct SpecializedGreensFunction { };

template<>
struct SpecializedGreensFunction<3> {
82 83 84 85 86 87 88 89 90 91 92 93 94 95 96
    template<class T, class FT, class FT2>
    static void calculate(Vektor<T, 3> &hrsq, FT &grn, FT2 *grnI, double alpha,double eps) {
        double r;
        NDIndex<3> elem0=NDIndex<3>(Index(0,0), Index(0,0),Index(0,0));
        grn = grnI[0] * hrsq[0] + grnI[1] * hrsq[1] + grnI[2] * hrsq[2];
        NDIndex<3> lDomain_m = grn.getLayout().getLocalNDIndex();
        NDIndex<3> elem;
        for (int i=lDomain_m[0].min(); i<=lDomain_m[0].max(); ++i) {
            elem[0]=Index(i,i);
            for (int j=lDomain_m[1].min(); j<=lDomain_m[1].max(); ++j) {
                elem[1]=Index(j,j);
                for (int k=lDomain_m[2].min(); k<=lDomain_m[2].max(); ++k) {
                    elem[2]=Index(k,k);
                    r = real(sqrt(grn.localElement(elem)));
                    if(elem==elem0) {
gsell's avatar
gsell committed
97
                        //grn.localElement(elem) = ke*std::complex<double>(2*alpha/sqrt(M_PI)) ;
98 99 100
                        grn.localElement(elem) = 0 ;
                    }
                    else
gsell's avatar
gsell committed
101
                        grn.localElement(elem) = ke*std::complex<double>(erf(alpha*r)/(r+eps));
102 103 104 105
                }
            }
        }
    }
Andreas Adelmann's avatar
Andreas Adelmann committed
106 107 108 109 110
};



template<class PL>
111
class ChargedParticles : public IpplParticleBase<PL> {
112
public:
113 114 115 116 117 118 119 120
    ParticleAttrib<double>      Q;
    ParticleAttrib<double>      m;
    ParticleAttrib<double>      Phi; //electrostatic potential
    ParticleAttrib<Vector_t>    EF;
    ParticleAttrib<Vector_t>    v; //velocity of the particles
    ParticleAttrib<int> ID; //velocity of the particles

    ChargedParticles(PL* pl, Vektor<double,3> nr, e_dim_tag /*decomp*/[Dim],Vektor<double,3> extend_l_, Vektor<double,3> extend_r_) :
121
    IpplParticleBase<PL>(pl),
122 123 124 125 126 127 128 129 130 131
    nr_m(nr),
    extend_l(extend_l_),
    extend_r(extend_r_)
    {
        this->addAttribute(Q);
        this->addAttribute(m);
        this->addAttribute(Phi);
        this->addAttribute(EF);
        this->addAttribute(v);
        this->addAttribute(ID);
132

133 134 135 136
        for (unsigned int i = 0; i < 2 * Dim; ++i) {
            //use periodic boundary conditions for the particles
            this->getBConds()[i] = ParticlePeriodicBCond;
            //boundary conditions used for interpolation kernels allow writes to ghost cells
137

138 139 140 141 142 143 144 145 146 147 148 149 150
            if (Ippl::getNodes()>1) {
                bc_m[i] = new ParallelInterpolationFace<double, Dim, Mesh_t, Center_t>(i);
                //std periodic boundary conditions for gradient computations etc.
                vbc_m[i] = new ParallelPeriodicFace<Vector_t, Dim, Mesh_t, Center_t>(i);
                bcp_m[i] = new ParallelPeriodicFace<double, Dim, Mesh_t, Center_t>(i);
            }
            else {
                bc_m[i] = new InterpolationFace<double, Dim, Mesh_t, Center_t>(i);
                //std periodic boundary conditions for gradient computations etc.
                vbc_m[i] = new PeriodicFace<Vector_t, Dim, Mesh_t, Center_t>(i);
                bcp_m[i] = new PeriodicFace<double, Dim, Mesh_t, Center_t>(i);
            }
        }
151

152 153 154 155
        for (unsigned int d = 0;d<Dim;++d) {
            rmax_m[d] = extend_r[d];
            rmin_m[d] = extend_l[d];
        }
156

157 158
        domain_m = this->getFieldLayout().getDomain();
        lDomain_m = this->getFieldLayout().getLocalNDIndex(); // local domain
159

160 161 162
        //initialize the FFT
        bool compressTemps = true;
        fft_m = new FFT_t(domain_m,compressTemps);
163

164 165 166 167
        fft_m->setDirectionName(+1, "forward");
        fft_m->setDirectionName(-1, "inverse");
        INFOMSG("INIT FFT DONE"<<endl);
    }
168

169
    inline const Mesh_t& getMesh() const { return this->getLayout().getLayout().getMesh(); }
170

171
    inline Mesh_t& getMesh() { return this->getLayout().getLayout().getMesh(); }
172

173 174 175
    inline const FieldLayout_t& getFieldLayout() const {
        return dynamic_cast<FieldLayout_t&>( this->getLayout().getLayout().getFieldLayout());
    }
176

177 178 179
    inline FieldLayout_t& getFieldLayout() {
        return dynamic_cast<FieldLayout_t&>(this->getLayout().getLayout().getFieldLayout());
    }
180

181 182 183 184 185 186 187 188 189
    void update()
    {
        //should only be needed if meshspacing changes -----------
        for (unsigned int d = 0;d<Dim;++d) {
            hr_m[d] = (extend_r[d] - extend_l[d]) / (nr_m[d]);
        }
        this->getMesh().set_meshSpacing(&(hr_m[0]));
        this->getMesh().set_origin(extend_l);
        //--------------------------------------------------------
190

191 192 193 194 195 196
        //init resets the meshes to 0 ?!
        rhocmpl_m.initialize(getMesh(), getFieldLayout(), GuardCellSizes<Dim>(1));
        grncmpl_m.initialize(getMesh(), getFieldLayout(), GuardCellSizes<Dim>(1));
        rho_m.initialize(getMesh(), getFieldLayout(), GuardCellSizes<Dim>(1),bc_m);
        phi_m.initialize(getMesh(), getFieldLayout(), GuardCellSizes<Dim>(1),bcp_m);
        eg_m.initialize(getMesh(), getFieldLayout(), GuardCellSizes<Dim>(1), vbc_m);
197

198 199
        domain_m = this->getFieldLayout().getDomain();
        lDomain_m = this->getFieldLayout().getLocalNDIndex();
200

201
        IpplParticleBase<PL>::update();
202
    }
203

204
    void compute_temperature() {
Andreas Adelmann's avatar
Andreas Adelmann committed
205 206 207 208
        Inform m("compute_temperature ");
        double loc_temp[Dim]={0.0,0.0,0.0};
        double loc_avg_vel[Dim]={0.0,0.0,0.0};

209
        for(unsigned long k = 0; k < this->getLocalNum(); ++k) {
210 211 212
          for(unsigned i = 0; i < Dim; i++) {
            loc_avg_vel[i]   += this->v[k](i);
          }
213 214 215
        }
        reduce(&(loc_avg_vel[0]), &(loc_avg_vel[0]) + Dim,
               &(avg_vel[0]), OpAddAssign());
Andreas Adelmann's avatar
Andreas Adelmann committed
216

217 218 219 220
        const double N =  static_cast<double>(this->getTotalNum());
        avg_vel[0]=avg_vel[0]/N;
        avg_vel[1]=avg_vel[1]/N;
        avg_vel[2]=avg_vel[2]/N;
221

Andreas Adelmann's avatar
Andreas Adelmann committed
222
        m << "avg_vel[0]= " << avg_vel[0] << " avg_vel[1]= " << avg_vel[1] << " avg_vel[2]= " << avg_vel[2] <<  endl;
223

224
        for(unsigned long k = 0; k < this->getLocalNum(); ++k) {
225 226 227
          for(unsigned i = 0; i < Dim; i++) {
            loc_temp[i]   += (this->v[k](i)-avg_vel[i])*(this->v[k](i)-avg_vel[i]);
          }
228 229 230 231 232 233 234
        }
        reduce(&(loc_temp[0]), &(loc_temp[0]) + Dim,
               &(temperature[0]), OpAddAssign());
        temperature[0]=temperature[0]/N;
        temperature[1]=temperature[1]/N;
        temperature[2]=temperature[2]/N;
    }
235

236 237
    void calcMoments() {
        double part[2 * Dim];
238

239 240 241
        double loc_centroid[2 * Dim]={};
        double loc_moment[2 * Dim][2 * Dim]={};
        double moments[2 * Dim][2 * Dim]={};
242

243 244 245 246 247 248 249
        for(unsigned i = 0; i < 2 * Dim; i++) {
            loc_centroid[i] = 0.0;
            for(unsigned j = 0; j <= i; j++) {
                loc_moment[i][j] = 0.0;
                loc_moment[j][i] = loc_moment[i][j];
            }
        }
250

251 252 253 254 255 256 257 258
        //double p0=m0*gamma*beta0;
        for(unsigned long k = 0; k < this->getLocalNum(); ++k) {
            part[1] = this->v[k](0);
            part[3] = this->v[k](1);
            part[5] = this->v[k](2);
            part[0] = this->R[k](0);
            part[2] = this->R[k](1);
            part[4] = this->R[k](2);
259

260 261 262 263 264 265 266
            for(unsigned i = 0; i < 2 * Dim; i++) {
                loc_centroid[i]   += part[i];
                for(unsigned j = 0; j <= i; j++) {
                    loc_moment[i][j]   += part[i] * part[j];
                }
            }
        }
267

268 269 270 271 272
        for(unsigned i = 0; i < 2 * Dim; i++) {
            for(unsigned j = 0; j < i; j++) {
                loc_moment[j][i] = loc_moment[i][j];
            }
        }
273

274 275
        reduce(&(loc_moment[0][0]), &(loc_moment[0][0]) + 2 * Dim * 2 * Dim,
               &(moments[0][0]), OpAddAssign());
276

277 278
        reduce(&(loc_centroid[0]), &(loc_centroid[0]) + 2 * Dim,
               &(centroid_m[0]), OpAddAssign());
279

280 281 282 283 284 285 286
        for(unsigned i = 0; i < 2 * Dim; i++) {
            for(unsigned j = 0; j <= i; j++) {
                moments_m[i][j] = moments[i][j];
                moments_m[j][i] = moments[i][j];
            }
        }
    }
287 288


289 290 291 292
    void computeBeamStatistics() {
        //const size_t locNp = this->getLocalNum();
        const double N =  static_cast<double>(this->getTotalNum());
        const double zero = 0.0;
293

294 295 296 297 298 299 300 301 302 303
        Vector_t eps2, fac, rsqsum, vsqsum, rvsum;
        for(unsigned int i = 0 ; i < Dim; i++) {
            rmean_m(i) = centroid_m[2 * i] / N;
            vmean_m(i) = centroid_m[(2 * i) + 1] / N;
            rsqsum(i) = moments_m[2 * i][2 * i] - N * rmean_m(i) * rmean_m(i);
            vsqsum(i) = moments_m[(2 * i) + 1][(2 * i) + 1] - N * vmean_m(i) * vmean_m(i);
            if(vsqsum(i) < 0)
                vsqsum(i) = 0;
            rvsum(i) = moments_m[(2 * i)][(2 * i) + 1] - N * rmean_m(i) * vmean_m(i);
        }
304

305 306
        eps2 = (rsqsum * vsqsum - rvsum * rvsum) / (N * N);
        rvsum /= N;
307

308 309 310 311 312 313 314 315
        for(unsigned int i = 0 ; i < Dim; i++) {
            rrms_m(i) = sqrt(rsqsum(i) / N);
            vrms_m(i) = sqrt(vsqsum(i) / N);
            eps_m(i)  =  std::sqrt(std::max(eps2(i), zero));
            double tmp = rrms_m(i) * vrms_m(i);
            fac(i) = (tmp == 0) ? zero : 1.0 / tmp;
        }
        rvrms_m = rvsum * fac;
316

317
    }
318

319 320 321 322 323 324 325 326 327
    void calc_kinetic_energy() {
        double loc_kinetic_energy=0;
        double v2;
        for (unsigned i=0; i<this->getLocalNum(); ++i) {
            v2 = (v[i][0]*v[i][0]+v[i][1]*v[i][1]+v[i][2]*v[i][2]);
            loc_kinetic_energy+=0.5*m[i]*v2;
        }
        reduce(loc_kinetic_energy,kinetic_energy, OpAddAssign());
    }
328

329 330 331 332 333
    void  calc_field_energy() {
        NDIndex<3> elem;
        double cell_volume = hr_m[0]*hr_m[1]*hr_m[2];
        field_energy=0;
        field_energy=0.5*cell_volume*sum(dot(eg_m,eg_m));
334

335 336 337 338
        rhomax=max(abs(rho_m))/(hr_m[0]*hr_m[1]*hr_m[2]);
        //rhomax=max(rho_m);
        integral_phi_m=0.5*sum(rho_m*phi_m);
    }
339

340 341 342 343 344 345
    void  calc_potential_energy() {
        potential_energy = 0;
        for (unsigned i=0; i<this->getLocalNum(); ++i) {
            potential_energy+=0.5*(Q[i])*Phi[i];
        }
    }
346

347 348 349 350 351 352
    void calc_Amplitude_E(){
        //computes the maximum amplitude in the electric field
        AmplitudeEfield=max(sqrt(dot(eg_m,eg_m)));
        eg_m=eg_m*Vektor<double,3>(0,0,1);
        AmplitudeEFz=max(sqrt(dot(eg_m,eg_m)));
    }
353

354
    void computeAvgSpaceChargeForces() {
Andreas Adelmann's avatar
Andreas Adelmann committed
355
        Inform m("computeAvgSpaceChargeForces ");
356

357 358 359 360 361 362 363
        const double N =  static_cast<double>(this->getTotalNum());
        double locAvgEF[Dim]={};
        for (unsigned i=0; i<this->getLocalNum(); ++i) {
            locAvgEF[0]+=fabs(EF[i](0));
            locAvgEF[1]+=fabs(EF[i](1));
            locAvgEF[2]+=fabs(EF[i](2));
        }
364

365 366
        reduce(&(locAvgEF[0]), &(locAvgEF[0]) + Dim,
               &(globSumEF[0]), OpAddAssign());
367 368


Andreas Adelmann's avatar
Andreas Adelmann committed
369
        m << "globSumEF = " << globSumEF[0] << "\t" << globSumEF[1] << "\t" << globSumEF[2] << endl;
370

371 372 373
        avgEF[0]=globSumEF[0]/N;
        avgEF[1]=globSumEF[1]/N;
        avgEF[2]=globSumEF[2]/N;
Andreas Adelmann's avatar
Andreas Adelmann committed
374

375
    }
376

377 378 379
    void applyConstantFocusing(double f,double beam_radius) {
        double focusingForce=sqrt(dot(avgEF,avgEF));
        for (unsigned i=0; i<this->getLocalNum(); ++i) {
Andreas Adelmann's avatar
Andreas Adelmann committed
380
            EF[i]+=this->R[i]/beam_radius*f*focusingForce;
381 382
        }
    }
383 384


385
    void calculatePairForces(double interaction_radius, double eps, double alpha);
386

387
    void calculateGridForces(double /*interaction_radius*/, double alpha, double eps, int /*it*/=0, bool /*normalizeSphere*/=0) {
388 389 390 391 392 393
        // (1) scatter charge to charge density grid and transform to fourier space
        //this->Q.scatter(this->rho_m, this->R, IntrplTSC_t());
        rho_m[domain_m]=0; //!!!!!! there has to be a better way than setting rho to 0 every time
        this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());
        //this->Q.scatter(this->rho_m, this->R, IntrplNGP_t());
        //dumpVTKScalar(rho_m,this,it,"RhoInterpol");
394

395 396 397
        //rhocmpl_m[domain_m] = rho_m[domain_m];
        rhocmpl_m[domain_m] = rho_m[domain_m]/(hr_m[0]*hr_m[1]*hr_m[2]);
        RhoSum=sum(real(rhocmpl_m));
398

399 400 401 402 403 404 405 406 407
        //std::cout << "total charge in densitty field before ion subtraction is" << sum(real(rhocmpl_m))<< std::endl;
        //std::cout << "max total charge in densitty field before ion subtraction is" << max(real(rhocmpl_m)) << std::endl;
        //subtract the background charge of the ions
        /*
         if (normalizeSphere)
         rhocmpl_m[domain_m]=1./(hr_m[0]*hr_m[1]*hr_m[2]*(nr_m[0]*nr_m[1]*nr_m[2]))+rhocmpl_m[domain_m];
         else
         rhocmpl_m[domain_m]=1.+rhocmpl_m[domain_m];
         */
408

409
        //std::cout << "total charge in densitty field after ion subtraction is" << sum(real(rhocmpl_m)) << std::endl;
410

411 412
        //compute rhoHat and store in rhocmpl_m
        fft_m->transform("inverse", rhocmpl_m);
413

414 415
        // (2) compute Greens function in real space and transform to fourier space
        //calcGrealSpace(alpha,eps);
416

417 418 419
        /////////compute G with Index Magic///////////////////
        // Fields used to eliminate excess calculation in greensFunction()
        IField_t grnIField_m[3];
420

421 422 423
        // mesh and layout objects for rho_m
        Mesh_t *mesh_m = &(getMesh());
        FieldLayout_t *layout_m = &(getFieldLayout());
424

425 426 427 428 429 430 431 432 433 434 435
        //This loop stores in grnIField_m[i] the index of the ith dimension mirrored at the central axis. e.g. grnIField_m[0]=[(0 1 2 3 ... 3 2 1) ; (0 1 2 3 ... 3 2 1; ...)]
        for (int i = 0; i < 3; ++i) {
            grnIField_m[i].initialize(*mesh_m, *layout_m);
            grnIField_m[i][domain_m] = where(lt(domain_m[i], nr_m[i]/2),
                                             domain_m[i] * domain_m[i],
                                             (nr_m[i]-domain_m[i]) *
                                             (nr_m[i]-domain_m[i]));
        }
        Vector_t hrsq(hr_m * hr_m);
        SpecializedGreensFunction<3>::calculate(hrsq, grncmpl_m, grnIField_m, alpha,eps);
        /////////////////////////////////////////////////
436

437 438 439 440
        //transform G -> Ghat and store in grncmpl_m
        fft_m->transform("inverse", grncmpl_m);
        //multiply in fourier space and obtain PhiHat in rhocmpl_m
        rhocmpl_m *= grncmpl_m;
441

442 443 444
        // (3) Backtransformation: compute potential field in real space and E=-Grad Phi
        //compute electrostatic potential Phi in real space by FFT PhiHat -> Phi and store it in rhocmpl_m
        fft_m->transform("forward", rhocmpl_m);
445

446 447 448
        //take only the real part and store in phi_m (has periodic bc instead of interpolation bc)
        phi_m = real(rhocmpl_m)*hr_m[0]*hr_m[1]*hr_m[2];
        //dumpVTKScalar( phi_m, this,it, "Phi_m") ;
449

450 451
        //compute Electric field on the grid by -Grad(Phi) store in eg_m
        eg_m = -Grad1Ord(phi_m, eg_m);
452

453 454 455 456 457
        //interpolate the electric field to the particle positions
        EF.gather(eg_m, this->R,  IntrplCIC_t());
        //interpolate electrostatic potenital to the particle positions
        Phi.gather(phi_m, this->R, IntrplCIC_t());
    }
458 459


460 461 462 463 464 465
    Vector_t getRmin() {
        return this->rmin_m;
    }
    Vector_t getRmax() {
        return this->rmax_m;
    }
466

Andreas Adelmann's avatar
Andreas Adelmann committed
467
    Vector_t get_hr() { return hr_m;}
468

469 470 471
    void closeH5(){
        H5CloseFile(H5f_m);
    }
472

473
    void openH5(std::string fn){
474 475 476
        h5_prop_t props = H5CreateFileProp ();
        MPI_Comm comm = Ippl::getComm();
        h5_err_t h5err = H5SetPropFileMPIOCollective (props, &comm);
frey_m's avatar
frey_m committed
477 478 479
#if defined (NDEBUG)
        (void)h5err;
#endif
480
        PAssert (h5err != H5_ERR);
481
        H5f_m = H5OpenFile (fn.c_str(), H5_O_RDONLY, props);
482
        PAssert (H5f_m != (h5_file_t)H5_ERR);
frey_m's avatar
frey_m committed
483
        H5CloseProp (props);
484
    }
485 486 487



488 489 490 491
    //private:
    BConds<double, Dim, Mesh_t, Center_t> bc_m;
    BConds<double, Dim, Mesh_t, Center_t> bcp_m;
    BConds<Vector_t, Dim, Mesh_t, Center_t> vbc_m;
492

493 494
    CxField_t rhocmpl_m;
    CxField_t grncmpl_m;
495

496 497
    Field_t rho_m;
    Field_t phi_m;
498

499
    VField_t eg_m;
500

501 502 503 504 505 506 507 508 509 510
    Vektor<int,Dim> nr_m;
    Vector_t hr_m;
    Vector_t rmax_m;
    Vector_t rmin_m;
    Vektor<double,Dim> extend_l;
    Vektor<double,Dim> extend_r;
    Mesh_t *mesh_m;
    FieldLayout_t *layout_m;
    NDIndex<Dim> domain_m;
    NDIndex<Dim> lDomain_m;
511

512 513 514 515 516 517 518 519 520
    double kinetic_energy;
    double field_energy;
    double field_energy_gather;
    double integral_phi_m;
    double potential_energy;
    double AmplitudeEfield;
    double AmplitudeEFz;
    double total_charge;
    double rhomax;
521

522
    FFT_t *fft_m;
523

524
    e_dim_tag decomp_m[Dim];
525

526 527
    Vektor<int,Dim> Nx;
    Vektor<int,Dim> Nv;
Andreas Adelmann's avatar
Andreas Adelmann committed
528
    Vektor<double,Dim> Vmax;
529 530 531 532 533 534 535 536
    //Fields for tracking distribution function
    Field2d_t f_m;
    Mesh2d_t mesh2d_m;
    NDIndex<2> domain2d_m;
    FieldLayout2d_t *layout2d_m;
    //TEMP debug variable
    double RhoSum=0;
    
537
    h5_file_t H5f_m;
538 539
    double temperature[Dim];
    double avg_vel[Dim];
540 541


542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559
    //Moment calculations:
    /// 6x6 matrix of the moments of the beam
    //FMatrix<double, 2 * Dim, 2 * Dim> moments_m;
    double moments_m[2*Dim][2*Dim];
    /// holds the centroid of the beam
    double centroid_m[2 * Dim];
    /// rms beam size (m)
    Vector_t rrms_m;
    /// rms momenta
    Vector_t vrms_m;
    /// mean position (m)
    Vector_t rmean_m;
    /// mean momenta
    Vector_t vmean_m;
    /// rms emittance (not normalized)
    Vector_t eps_m;
    /// rms correlation
    Vector_t rvrms_m;
560 561


562 563
    Vektor<double,Dim> avgEF;
    double globSumEF[Dim];
564

Andreas Adelmann's avatar
Andreas Adelmann committed
565 566 567 568
};

template<class T>
struct ApplyField {
569 570 571 572 573
    ApplyField(T c, double r, double epsilon, double alpha) : C(c), R(r), eps(epsilon), a(alpha) {}
    void operator()(std::size_t i, std::size_t j, ChargedParticles<playout_t> &P,Vektor<double,3> &shift) const
    {
        Vector_t diff = P.R[i] - (P.R[j]+shift);
        double sqr = 0;
574

575 576
        for (unsigned d = 0; d<Dim; ++d)
            sqr += diff[d]*diff[d];
577

578 579 580 581 582 583 584
        //compute r with softening parameter, unsoftened r is obtained by sqrt(sqr)
        if(sqr!=0) {
            double r = std::sqrt(sqr+eps*eps);
            //for order two transition
            if (P.Q[i]!=0 && P.Q[j]!=0) {
                //compute potential energy
                double phi =ke*(1.-erf(a*sqrt(sqr)))/r;
585

586 587
                //compute force
                Vector_t Fij = ke*C*(diff/sqrt(sqr))*((2.*a*exp(-a*a*sqr))/(sqrt(M_PI)*r)+(1.-erf(a*sqrt(sqr)))/(r*r));
588

589 590 591 592 593 594 595 596 597 598 599 600 601 602
                //Actual Force is F_ij multiplied by Qi*Qj
                //The electrical field on particle i is E=F/q_i and hence:
                P.EF[i] -= P.Q[j]*Fij;
                P.EF[j] += P.Q[i]*Fij;
                //update potential per particle
                P.Phi[i] += P.Q[j]*phi;
                P.Phi[j] += P.Q[i]*phi;
            }
        }
    }
    T C;
    double R;
    double eps;
    double a;
Andreas Adelmann's avatar
Andreas Adelmann committed
603 604
};

605 606 607 608 609 610 611 612 613 614 615 616 617 618 619
template<class PL>
void ChargedParticles<PL>::calculatePairForces(double interaction_radius, double eps, double alpha)
{
    if (interaction_radius>0){
        if (Ippl::getNodes() > 1) {
            HashPairBuilderPeriodicParallel< ChargedParticles<playout_t> > HPB(*this);
            HPB.for_each(RadiusCondition<double, Dim>(interaction_radius), ApplyField<double>(-1,interaction_radius,eps,alpha),extend_l, extend_r);
        }
        else {
            HashPairBuilderPeriodic< ChargedParticles<playout_t> > HPB(*this);
            HPB.for_each(RadiusCondition<double, Dim>(interaction_radius), ApplyField<double>(-1,interaction_radius,eps,alpha),extend_l, extend_r);
        }
    }
}

Andreas Adelmann's avatar
Andreas Adelmann committed
620
int main(int argc, char *argv[]){
621 622 623
    Ippl ippl(argc, argv);
    Inform msg(argv[0]);
    Inform msg2all(argv[0],INFORM_ALL_NODES);
624

625 626
    IpplTimings::TimerRef allTimer = IpplTimings::getTimer("AllTimer");
    IpplTimings::startTimer(allTimer);
627

628
    Vektor<int,Dim> nr;
629

630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645
    nr = Vektor<int,Dim>(atoi(argv[1]),atoi(argv[2]),atoi(argv[3]));
    int param = 4;
    double beam_radius =atof(argv[param++]);
    double box_length =atof(argv[param++]);
    //double part_density =atof(argv[param++]);
    int Nparticle =atoi(argv[param++]);
    double interaction_radius = atof(argv[param++]);
    //read the remaining sim params
    double alpha =atof(argv[param++]);
    double dt = atof(argv[param++]);
    double eps = atof(argv[param++]);
    int iterations =  atoi(argv[param++]);
    double charge_per_part =  atof(argv[param++]);
    double mass_per_part =  atof(argv[param++]);
    double focusingForce =  atof(argv[param++]);
    int print_every =  atof(argv[param++]);
646 647


648 649 650 651 652
    ///////// setup the initial layout ///////////////////////////////////////
    e_dim_tag decomp[Dim];
    Mesh_t *mesh;
    FieldLayout_t *FL;
    ChargedParticles<playout_t>  *P;
653

654 655 656
    NDIndex<Dim> domain;
    for (unsigned i=0; i<Dim; i++)
        domain[i] = domain[i] = Index(nr[i]+1);
657

658 659
    for (unsigned d=0; d < Dim; ++d)
        decomp[d] = PARALLEL;
660

661 662 663 664
    // create mesh and layout objects for this problem domain
    mesh          = new Mesh_t(domain);
    FL            = new FieldLayout_t(*mesh, decomp);
    playout_t* PL = new playout_t(*FL, *mesh);
665

Andreas Adelmann's avatar
Andreas Adelmann committed
666 667
    PL->setAllCacheDimensions(interaction_radius);
    PL->enableCaching();
668

669 670 671 672
    /////////// Create the particle distribution /////////////////////////////////////////////////////
    double L = box_length/2.;
    Vektor<double,Dim> extend_l(-L,-L,-L);
    Vektor<double,Dim> extend_r(L,L,L);
673

Andreas Adelmann's avatar
Andreas Adelmann committed
674
    Vektor<double,Dim> Vmax(6,6,6);
675 676
    P = new ChargedParticles<playout_t>(PL, nr, decomp, extend_l, extend_r);
    createParticleDistributionHeating(P,extend_l,extend_r,beam_radius, Nparticle,charge_per_part,mass_per_part);
677

678 679 680 681
    //COmpute and write temperature
    P->compute_temperature();
    writeTemperature(P,0);
    /////////////////////////////////////////////////////////////////////////////////////////////
682

683 684 685 686 687
    /////// Print mesh informations ////////////////////////////////////////////////////////////
    INFOMSG(P->getMesh() << endl);
    INFOMSG(P->getFieldLayout() << endl);
    msg << endl << endl;
    Ippl::Comm->barrier();
688

Andreas Adelmann's avatar
Andreas Adelmann committed
689
    //dumpParticlesCSV(P,0);
690

691 692 693
    INFOMSG(P->getMesh() << endl);
    INFOMSG(P->getFieldLayout() << endl);
    msg << endl << endl;
694

Andreas Adelmann's avatar
Andreas Adelmann committed
695 696 697
    msg<<"number of particles = " << P->getTotalNum() << endl;
    msg<<"Total charge Q      = " << P->total_charge << endl;

698 699 700 701
    ////////////////////////////////////////////////////////////////////////////////////////////
    std::string fname;
    fname = "data/particleData";
    fname += ".h5part";
702

703 704 705
    P->openH5(fname);
    dumpH5partVelocity(P,0);
    unsigned printid=1;
706

Andreas Adelmann's avatar
Andreas Adelmann committed
707 708
    msg << "Starting iterations ..." << endl;
    P->compute_temperature();
709
    // calculate initial space charge forces
710 711
    P->calculateGridForces(interaction_radius,alpha,0,0,0);
    P->calculatePairForces(interaction_radius,eps,alpha);
712

713
    //avg space charge forces for constant focusing
714
    P->computeAvgSpaceChargeForces();
Andreas Adelmann's avatar
Andreas Adelmann committed
715

716
    //dumpVTKVector(P->eg_m, P,0,"EFieldAfterPMandPP");
717

718 719 720 721 722 723 724
    //compute quantities to check correctness:
    /*
     P->calc_field_energy();
     P->calc_potential_energy();
     P->calc_kinetic_energy();
     writeEnergy(P,0);
     */
725

Andreas Adelmann's avatar
Andreas Adelmann committed
726 727
    IpplTimings::TimerRef gridTimer = IpplTimings::getTimer("GridTimer");
    IpplTimings::TimerRef particleTimer = IpplTimings::getTimer("ParticleTimer");
728

729
    for (int it=0; it<iterations; it++) {
Andreas Adelmann's avatar
Andreas Adelmann committed
730
      /*
731 732 733
        P->calcMoments();
        P->computeBeamStatistics();
        writeBeamStatisticsVelocity(P,it);
734

735 736 737
        P->calc_kinetic_energy();
        P->calc_field_energy();
        writeEnergy(P,it);
738
      */
739 740 741
        // advance the particle positions
        // basic leapfrogging timestep scheme.  velocities are offset
        // by half a timestep from the positions.
742

743 744 745
        assign(P->R, P->R + dt * P->v);
        // update particle distribution across processors
        P->update();
746

747
        // compute the electric field
748 749


Andreas Adelmann's avatar
Andreas Adelmann committed
750 751
        IpplTimings::startTimer(gridTimer);
        P->calculateGridForces(interaction_radius,alpha,0,it+1,0);
752
        IpplTimings::stopTimer(gridTimer);
753

754 755 756
        IpplTimings::startTimer(particleTimer);
        P->calculatePairForces(interaction_radius,eps,alpha);
        IpplTimings::stopTimer(particleTimer);
757

758
        //P->update();
759

760 761 762 763
        //second part of leapfrog: advance velocitites
        //P->computeAvgSpaceChargeForces();
        //if (Ippl::myNode()==0)
        //std::cout <<"avg E-Field = " << P->avgEF << std::endl;
Andreas Adelmann's avatar
Andreas Adelmann committed
764

765
        P->applyConstantFocusing(focusingForce,beam_radius);
Andreas Adelmann's avatar
Andreas Adelmann committed
766

767
        assign(P->v, P->v + dt * P->Q/P->m * (P->EF));
768

Andreas Adelmann's avatar
Andreas Adelmann committed
769 770
        P->compute_temperature();

771 772 773 774 775 776 777 778
        if (it%print_every==0){
            //dumpConservedQuantities(P,printid);
            //compute quantities
            /*
             P->calc_field_energy();
             P->calc_kinetic_energy();
             P->calc_potential_energy();
             writeEnergy(P,printid);
779
             */
780 781
            P->compute_temperature();
            writeTemperature(P,it+1);
782

783 784
            dumpH5partVelocity(P,printid++);
        }
785

786 787 788
        msg << "Finished iteration " << it << endl;
    }
    Ippl::Comm->barrier();
789

790 791
    P->closeH5();
    Ippl::Comm->barrier();
792

793
    IpplTimings::stopTimer(allTimer);
794

795
    IpplTimings::print();
796 797


798 799 800
    delete P;
    delete FL;
    delete mesh;
801

802
    return 0;
Andreas Adelmann's avatar
Andreas Adelmann committed
803
}