35 const double ke=2.532638e8;
37 template<
unsigned int Dim>
42 template<
class T,
class FT,
class FT2>
46 grn = grnI[0] * hrsq[0] + grnI[1] * hrsq[1] + grnI[2] * hrsq[2];
47 NDIndex<3> lDomain_m = grn.getLayout().getLocalNDIndex();
49 for (
int i=lDomain_m[0].
min(); i<=lDomain_m[0].max(); ++i) {
51 for (
int j=lDomain_m[1].
min(); j<=lDomain_m[1].max(); ++j) {
53 for (
int k=lDomain_m[2].
min(); k<=lDomain_m[2].max(); ++k) {
55 r =
real(
sqrt(grn.localElement(elem)));
58 grn.localElement(elem) = 0 ;
76 for (
unsigned d = 0; d<
Dim; ++d)
77 sqr += diff[d]*diff[d];
83 if (P.
Q[i]!=0 && P.
Q[j]!=0) {
92 P.
Ef[i] -= P.
Q[j]*Fij;
93 P.
Ef[j] += P.
Q[i]*Fij;
95 P.
Phi[i] += P.
Q[j]*phi;
96 P.
Phi[j] += P.
Q[i]*phi;
115 interaction_radius_m(interaction_radius),
119 Inform msg(
"P3MPoissonSolver::P3MPoissonSolver ");
123 for (
unsigned int i = 0; i < 2 *
Dim; ++i) {
138 for (
unsigned int i = 0; i < 2 *
Dim; ++i) {
167 for(
unsigned int i = 0; i < 3; i++) {
179 bool compressTemps =
true;
182 fft_m->setDirectionName(+1,
"forward");
183 fft_m->setDirectionName(-1,
"inverse");
190 for (
unsigned int i = 0; i < 2 *
Dim; ++i) {
205 for (
unsigned int i = 0; i < 2 *
Dim; ++i) {
220 for(
unsigned int i = 0; i < 3; i++) {
232 bool compressTemps =
true;
238 fft_m->setDirectionName(+1,
"forward");
239 fft_m->setDirectionName(-1,
"inverse");
251 if (interaction_radius>0){
268 Inform msg (
"calculateGridForces ");
284 for (
int i = 0; i < 3; ++i) {
330 Inform m(
"computeAvgSpaceChargeForces ");
332 const double N =
static_cast<double>(bunch->
getTotalNum());
333 double locAvgEf[
Dim]={};
335 locAvgEf[0]+=
fabs(bunch->
Ef[i](0));
336 locAvgEf[1]+=
fabs(bunch->
Ef[i](1));
337 locAvgEf[2]+=
fabs(bunch->
Ef[i](2));
340 reduce(&(locAvgEf[0]), &(locAvgEf[0]) +
Dim,
357 Rrel=bunch->
R[i] - beam_center;
358 bunch->
Ef[i] += Rrel/r*f*scFoc;
372 Inform m(
"compute_temperature ");
374 double loc_avg_vel[
Dim]={0.0,0.0,0.0};
375 double avg_vel[
Dim]={0.0,0.0,0.0};
377 for(
unsigned long k = 0; k < bunch->
getLocalNum(); ++k) {
378 for(
unsigned i = 0; i <
Dim; i++) {
379 loc_avg_vel[i] += bunch->
P[k](i);
382 reduce(&(loc_avg_vel[0]), &(loc_avg_vel[0]) +
Dim,
385 const double N =
static_cast<double>(bunch->
getTotalNum());
386 avg_vel[0]=avg_vel[0]/N;
387 avg_vel[1]=avg_vel[1]/N;
388 avg_vel[2]=avg_vel[2]/N;
390 m <<
"avg_vel[0]= " << avg_vel[0] <<
" avg_vel[1]= " << avg_vel[1] <<
" avg_vel[2]= " << avg_vel[2] <<
endl;
407 Inform msg(
"P3MPoissonSolver::test ");
410 const double mi = 1.0;
411 const double qi = -1.0;
412 const double qom = qi/mi;
413 const double beam_radius = 0.001774;
414 const double f = 1.5;
415 const double dt = bunch->
getdT();
430 for (
int i=0; i<3; i++) {
435 msg << *
this <<
endl;
444 for (
int it=0; it<1000; it++) {
450 assign(bunch->
R, bunch->
R + dt * bunch->
P);
458 assign(bunch->
P, bunch->
P + dt * qom * bunch->
Ef);
464 msg <<
"Finished iteration " << it <<
endl;
470 os <<
"* ************* P 3 M - P o i s s o n S o l v e r *************** " <<
endl;
471 os <<
"* h " <<
hr_m <<
'\n';
473 os <<
"* ALPHA " <<
alpha_m <<
'\n';
474 os <<
"* EPSILON " <<
eps_m <<
'\n';
475 os <<
"* Extend L " <<
extend_l <<
'\n';
476 os <<
"* Extend R " <<
extend_r <<
'\n';
477 os <<
"* hr " << hr_m <<
'\n';
478 os <<
"* nr " <<
nr_m <<
'\n';
479 os <<
"* *************************************************************** " <<
endl;
ParticleAttrib< Vector_t > P
double interaction_radius_m
The global OPAL structure.
ParticleAttrib< Vector_t > Ef
void applyConstantFocusing(PartBunchBase< double, 3 > *bunch, double f, double r)
void scatter(Field< T, Dim, M, C > &f, const ParticleAttrib< Vektor< PT, Dim > > &pp, const IntOp &intop) const
FTps< T, N > erf(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Error function.
static void calculate(Vektor< T, 3 > &hrsq, FT &grn, FT2 *grnI, double alpha, double eps)
PETE_TUTree< FnFabs, typename T::PETE_Expr_t > fabs(const PETE_Expr< T > &l)
ParticleAttrib< double > Q
FFT< CCTransform, 3, double > FFTC_t
Tps< T > exp(const Tps< T > &x)
Exponential.
std::complex< double > dcomplex
void for_each(const Pred &pred, const OP &op, Vektor< double, 3 > extend_l, Vektor< double, 3 > extend_r)
P3MPoissonSolver(Mesh_t *mesh, FieldLayout_t *fl, double interaction_radius, double alpha, double eps)
void test(PartBunchBase< double, 3 > *bunch)
BConds< Vector_t, Dim, Mesh_t, Center_t > vbc_m
void assign(const BareField< T, Dim > &a, RHS b, OP op, ExprTag< true >)
size_t getTotalNum() const
double dot(const Vector3D &lhs, const Vector3D &rhs)
Vector dot product.
FLieGenerator< T, N > real(const FLieGenerator< std::complex< T >, N > &)
Take real part of a complex generator.
constexpr double alpha
The fine structure constant, no dimension.
bool reduce(Communicate &, InputIterator, InputIterator, OutputIterator, const ReduceOp &, bool *IncludeVal=0)
static OpalData * getInstance()
void compute_temperature(PartBunchBase< double, 3 > *bunch)
void calcBeamParameters()
constexpr double c
The velocity of light in m/s.
PETE_TBTree< OpLT, Index::PETE_Expr_t, PETE_Scalar< double > > lt(const Index &idx, double x)
void for_each(const Pred &pred, const OP &op, Vektor< double, 3 > extend_l, Vektor< double, 3 > extend_r)
PETE_TTTree< OpWhere, typename Cond_t::PETE_Expr_t, typename True_t::PETE_Expr_t, PETE_Scalar< Vektor< T, Dim > > > where(const PETE_Expr< Cond_t > &c, const PETE_Expr< True_t > &t, const Vektor< T, Dim > &f)
void computeAvgSpaceChargeForces(PartBunchBase< double, 3 > *bunch)
std::unique_ptr< FFTC_t > fft_m
void computePotential(Field_t &rho, Vector_t hr, double zshift)
BConds< double, Dim, Mesh_t, Center_t > bcp_m
size_t getLocalNum() const
void operator()(std::size_t i, std::size_t j, PartBunch &P, Vektor< double, 3 > &shift) const
Tps< T > sqrt(const Tps< T > &x)
Square root.
Vektor< double, 3 > extend_l
Inform & print(Inform &os) const
ParticleAttrib< double > Phi
BConds< double, Dim, Mesh_t, Center_t > bc_m
Field< Vektor< T, 3U >, 3U, UniformCartesian< 3U, MFLOAT >, Cell > & Grad1Ord(Field< T, 3U, UniformCartesian< 3U, MFLOAT >, Vert > &x, Field< Vektor< T, 3U >, 3U, UniformCartesian< 3U, MFLOAT >, Cell > &r)
Old Grad operator.
MFLOAT get_meshSpacing(unsigned d) const
void calculatePairForces(PartBunchBase< double, 3 > *bunch, double interaction_radius, double alpha, double eps)
ApplyField(T c, double r, double epsilon, double alpha)
void gather(const Field< T, Dim, M, C > &f, const ParticleAttrib< Vektor< PT, Dim > > &pp, const IntOp &)
void calculateGridForces(PartBunchBase< double, 3 > *bunch, double interaction_radius, double alpha, double eps)
void dumpSDDS(PartBunchBase< double, 3 > *beam, Vector_t FDext[], const double &azimuth=-1) const
void initialize(Layout_t &)
ParticleAttrib< double > M
Vektor< double, Dim > avgEF_m
Vektor< double, 3 > extend_r
T::PETE_Expr_t::PETE_Return_t min(const PETE_Expr< T > &expr, NDIndex< D > &loc)
const NDIndex< Dim > & getDomain() const
Inform & endl(Inform &inf)