P3MPoissonSolver.h 4.29 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
//
// Class P3MPoissonSolver
//   This class contains methods for solving Poisson's equation for the
//   space charge portion of the calculation.
//
// 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
23 24
#ifndef P3M_POISSON_SOLVER_H_
#define P3M_POISSON_SOLVER_H_
25
const unsigned Dim = 3;
Andreas Adelmann's avatar
Andreas Adelmann committed
26 27 28 29 30 31 32 33 34 35

#ifdef dontOPTIMIZE_FIELD_ASSIGNMENT
#define FIELDASSIGNOPTIMIZATION __attribute__((optimize(0)))
#else
#define FIELDASSIGNOPTIMIZATION
#endif

#include <memory>
//////////////////////////////////////////////////////////////
#include "PoissonSolver.h"
frey_m's avatar
frey_m committed
36 37 38 39
#include "Algorithms/PartBunchBase.h"

template <class T, unsigned Dim>
class PartBunchBase;
Andreas Adelmann's avatar
Andreas Adelmann committed
40 41 42 43 44 45

//////////////////////////////////////////////////////////////

class P3MPoissonSolver : public PoissonSolver {
public:
    // constructor and destructor
46
    P3MPoissonSolver(Mesh_t *mesh, FieldLayout_t *fl, double interaction_radius, double alpha, double eps);
Andreas Adelmann's avatar
Andreas Adelmann committed
47 48 49

    ~P3MPoissonSolver();

50
    void initFields();
51

frey_m's avatar
frey_m committed
52
    void calculateGridForces(PartBunchBase<double, 3> *bunch, double interaction_radius, double alpha, double eps);
53

frey_m's avatar
frey_m committed
54
    void calculatePairForces(PartBunchBase<double, 3> *bunch, double interaction_radius, double alpha, double eps);
Andreas Adelmann's avatar
Andreas Adelmann committed
55 56 57 58 59 60 61 62 63

    // given a charge-density field rho and a set of mesh spacings hr,
    // compute the scalar potential with image charges at  -z
    void computePotential(Field_t &rho, Vector_t hr, double zshift);

    // given a charge-density field rho and a set of mesh spacings hr,
    // compute the scalar potential in open space
    void computePotential(Field_t &rho, Vector_t hr);

frey_m's avatar
frey_m committed
64 65
    void applyConstantFocusing(PartBunchBase<double, 3> *bunch, double f, double r);
    void test(PartBunchBase<double, 3> *bunch);
66

67 68 69 70 71 72
    double getXRangeMin(unsigned short /*level*/) {return 1.0;}
    double getXRangeMax(unsigned short /*level*/) {return 1.0;}
    double getYRangeMin(unsigned short /*level*/) {return 1.0;}
    double getYRangeMax(unsigned short /*level*/) {return 1.0;}
    double getZRangeMin(unsigned short /*level*/) {return 1.0;}
    double getZRangeMax(unsigned short /*level*/) {return 1.0;}
Andreas Adelmann's avatar
Andreas Adelmann committed
73

frey_m's avatar
frey_m committed
74 75
    void computeAvgSpaceChargeForces(PartBunchBase<double, 3> *bunch);
    void compute_temperature(PartBunchBase<double, 3> *bunch);
Andreas Adelmann's avatar
Andreas Adelmann committed
76 77
    Inform &print(Inform &os) const;
private:
78

79 80 81
    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;
82

83 84 85
    // rho_m is the charge-density field with mesh doubled in each dimension
    Field_t rho_m;
    Field_t phi_m;
86

87
    VField_t eg_m;
88

Andreas Adelmann's avatar
Andreas Adelmann committed
89 90 91
    // real field with layout of complex field: domain3_m
    Field_t greentr_m;

92 93
    CxField_t rhocmpl_m;
    CxField_t grncmpl_m;
94

Andreas Adelmann's avatar
Andreas Adelmann committed
95 96 97 98
    // grntr_m is the Fourier transformed Green's function
    // domain3_m and mesh3_ are used
    CxField_t grntr_m;

99 100
    // the FFT object
    std::unique_ptr<FFTC_t> fft_m;
101

Andreas Adelmann's avatar
Andreas Adelmann committed
102 103 104 105 106

    // Fields used to eliminate excess calculation in greensFunction()
    // mesh2_m and layout2_m are used
    IField_t grnIField_m[3];

107

Andreas Adelmann's avatar
Andreas Adelmann committed
108 109 110 111
    // mesh and layout objects for rho_m
    Mesh_t *mesh_m;
    FieldLayout_t *layout_m;

112

Andreas Adelmann's avatar
Andreas Adelmann committed
113 114 115 116 117 118
    // tmp
    Field_t tmpgreen;

    // domains for the various fields
    NDIndex<3> domain_m;             // original domain, gridsize
    // mesh and gridsize defined outside of P3M class, given as
119 120


Andreas Adelmann's avatar
Andreas Adelmann committed
121
    NDIndex<3> domainP3MConstruct_m;
122 123


124 125 126
    double interaction_radius_m;
    double alpha_m;
    double eps_m;
127

Andreas Adelmann's avatar
Andreas Adelmann committed
128 129
    Vector_t hr_m;
    Vektor<int, 3> nr_m;
130

131 132 133
    // for tests
    Vektor<double,Dim> avgEF_m;
    double globSumEf_m[Dim];
Andreas Adelmann's avatar
Andreas Adelmann committed
134 135


136 137 138
public:
    Vektor<double,3> extend_l;
    Vektor<double,3> extend_r;
Andreas Adelmann's avatar
Andreas Adelmann committed
139

140

Andreas Adelmann's avatar
Andreas Adelmann committed
141 142 143 144 145 146 147 148
};

inline Inform &operator<<(Inform &os, const P3MPoissonSolver &fs) {
    return fs.print(os);
}



frey_m's avatar
frey_m committed
149
#endif