portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 07:55:35 +00:00
Code Releases Activity

Minor cosmetic changes

Browse Source
This commit is contained in:
Guido Cossu 2017-09-15 11:48:36 +01:00
parent 91eaace19d
commit b542d349b8
1 changed files with 68 additions and 69 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

137
lib/qcd/action/scalar/ScalarImpl.h
View File

@ -91,6 +91,9 @@ class ScalarImplTypes {
#define USE_FFT_ACCELERATION #define USE_FFT_ACCELERATION
#ifdef USE_FFT_ACCELERATION
#define FFT_MASS 0.707
#endif
template <class S, unsigned int N> template <class S, unsigned int N>
@ -113,113 +116,109 @@ class ScalarImplTypes {
typedef Field FermionField; typedef Field FermionField;
typedef Field PropagatorField; typedef Field PropagatorField;
static void MomentaSquare(ComplexField &out)
{
GridBase *grid = out._grid;
const std::vector<int> &l = grid->FullDimensions();
ComplexField kmu(grid);
static void MomentaSquare(ComplexField& out){ for (int mu = 0; mu < grid->Nd(); mu++)
GridBase *grid = out._grid;
const std::vector<int> &l = grid->FullDimensions();
ComplexField kmu(grid);
for(int mu = 0; mu < grid->Nd(); mu++)
{ {
Real twoPiL = M_PI*2.0/l[mu]; Real twoPiL = M_PI * 2.0 / l[mu];
LatticeCoordinate(kmu,mu); LatticeCoordinate(kmu, mu);
kmu = 2.0*sin(0.5*twoPiL*kmu); kmu = 2.0 * sin(0.5 * twoPiL * kmu);
out += kmu*kmu; out += kmu * kmu;
} }
} }
static void MomentumSpacePropagator(ComplexField &out, RealD m) static void MomentumSpacePropagator(ComplexField &out, RealD m)
{ {
GridBase *grid = out._grid; GridBase *grid = out._grid;
ComplexField one(grid); one = Complex(1.0,0.0); ComplexField one(grid);
out = m*m; one = Complex(1.0, 0.0);
out = m * m;
MomentaSquare(out); MomentaSquare(out);
out = one/out; out = one / out;
} }
static inline void generate_momenta(Field &P, GridParallelRNG &pRNG)
static inline void generate_momenta(Field& P, GridParallelRNG& pRNG) { {
#ifndef USE_FFT_ACCELERATION #ifndef USE_FFT_ACCELERATION
Group::GaussianFundamentalLieAlgebraMatrix(pRNG, P); Group::GaussianFundamentalLieAlgebraMatrix(pRNG, P);
#else #else
Field Ptmp(P._grid), Pp(P._grid);
Group::GaussianFundamentalLieAlgebraMatrix(pRNG, Ptmp);
// if we change the mass I need a renormalization here
// transform and multiply by (M*M+p*p)^-1
GridCartesian *Grid = dynamic_cast<GridCartesian*>(P._grid);
FFT theFFT(Grid);
ComplexField p2(Grid);
RealD M = 1.0;
p2= zero;
theFFT.FFT_all_dim(Pp,Ptmp,FFT::forward); Field Pgaussian(P._grid), Pp(P._grid);
ComplexField p2(P._grid); p2 = zero;
RealD M = FFT_MASS;
Group::GaussianFundamentalLieAlgebraMatrix(pRNG, Pgaussian);
FFT theFFT((GridCartesian*)P._grid);
theFFT.FFT_all_dim(Pp, Pgaussian, FFT::forward);
MomentaSquare(p2); MomentaSquare(p2);
p2 += M*M; p2 += M * M;
p2 = sqrt(p2); p2 = sqrt(p2);
Pp *= p2; Pp *= p2;
theFFT.FFT_all_dim(P,Pp,FFT::backward); theFFT.FFT_all_dim(P, Pp, FFT::backward);
#endif //USE_FFT_ACCELERATION #endif //USE_FFT_ACCELERATION
} }
static inline Field projectForce(Field& P) {return P;} static inline Field projectForce(Field& P) {return P;}
static inline void update_field(Field& P, Field& U, double ep) { static inline void update_field(Field &P, Field &U, double ep)
#ifndef USE_FFT_ACCELERATION {
U += P*ep; #ifndef USE_FFT_ACCELERATION
#else U += P * ep;
// Here we can eventually add the Fourier acceleration #else
// FFT transform P(x) -> P(p) // FFT transform P(x) -> P(p)
// divide by (M^2+p^2) M external parameter (how to pass?) // divide by (M^2+p^2) M external parameter (how to pass?)
// P'(p) = P(p)/(M^2+p^2) // P'(p) = P(p)/(M^2+p^2)
// Transform back -> P'(x) // Transform back -> P'(x)
// U += P'(x)*ep // U += P'(x)*ep
// the dynamic cast is safe
GridCartesian *Grid = dynamic_cast<GridCartesian*>(U._grid);
FFT theFFT(Grid);
Field Pp(Grid), Pnew(Grid);
std::vector<int> full_dim = Grid->FullDimensions();
theFFT.FFT_all_dim(Pp,P,FFT::forward); Field Pp(U._grid), P_FFT(U._grid);
RealD M = 1.0; static ComplexField p2(U._grid);
RealD M = FFT_MASS;
FFT theFFT((GridCartesian*)U._grid);
theFFT.FFT_all_dim(Pp, P, FFT::forward);
static bool first_call = true; static bool first_call = true;
static ComplexField p2(Grid); if (first_call)
if (first_call){ {
MomentumSpacePropagator(p2,M); // avoid recomputing
first_call = false; MomentumSpacePropagator(p2, M);
first_call = false;
} }
Pp *= p2; Pp *= p2;
theFFT.FFT_all_dim(Pnew,Pp,FFT::backward); theFFT.FFT_all_dim(P_FFT, Pp, FFT::backward);
U += Pnew * ep; U += P_FFT * ep;
#endif //USE_FFT_ACCELERATION #endif //USE_FFT_ACCELERATION
} }
static inline RealD FieldSquareNorm(Field &U) static inline RealD FieldSquareNorm(Field &U)
{ {
#ifndef USE_FFT_ACCELERATION #ifndef USE_FFT_ACCELERATION
return (TensorRemove(sum(trace(U*U))).real()); return (TensorRemove(sum(trace(U * U))).real());
#else #else
// In case of Fourier acceleration we have to: // In case of Fourier acceleration we have to:
// compute U(p)*U(p)/(M^2+p^2)) Parseval theorem // compute U(p)*U(p)/(M^2+p^2)) Parseval theorem
// 1 FFT needed U(x) -> U(p) // 1 FFT needed U(x) -> U(p)
// M to be passed // M to be passed
GridCartesian *Grid = dynamic_cast<GridCartesian *>(U._grid); FFT theFFT((GridCartesian*)U._grid);
FFT theFFT(Grid); Field Up(U._grid);
Field Up(Grid), Utilde(Grid);
std::vector<int> full_dim = Grid->FullDimensions();
theFFT.FFT_all_dim(Up, U, FFT::forward); theFFT.FFT_all_dim(Up, U, FFT::forward);
RealD M = 1.0; RealD M = FFT_MASS;
ComplexField p2(Grid); ComplexField p2(U._grid);
MomentumSpacePropagator(p2,M); MomentumSpacePropagator(p2, M);
Field Up2 = Up*p2; Field Up2 = Up * p2;
// from the definition of the DFT we need to divide by the volume // from the definition of the DFT we need to divide by the volume
return (-TensorRemove(sum(trace(adj(Up)*Up2))).real()/U._grid->gSites()); return (-TensorRemove(sum(trace(adj(Up) * Up2))).real() / U._grid->gSites());
#endif //USE_FFT_ACCELERATION #endif //USE_FFT_ACCELERATION
} }
static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) { static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) {