1
0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 15:55:37 +00:00

Scalar action

This commit is contained in:
Antonin Portelli 2016-12-29 22:42:58 +01:00
parent 4c3fd9fa3f
commit 8c3cc32364
2 changed files with 216 additions and 0 deletions

View File

@ -292,4 +292,9 @@ typedef MobiusFermion<GparityWilsonImplD> GparityMobiusFermionD;
#include <Grid/qcd/action/pseudofermion/OneFlavourEvenOddRational.h>
#include <Grid/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h>
////////////////////
// Scalar actions
////////////////////
#include <Grid/qcd/action/scalar/Scalar.h>
#endif

View File

@ -0,0 +1,211 @@
#ifndef QCD_SCALAR_ACTION_H
#define QCD_SCALAR_ACTION_H
#define INHERIT_SIMPL_TYPES(Impl)\
typedef typename Impl::SiteScalar SiteScalar; \
typedef typename Impl::SiteSpinor SiteSpinor; \
typedef typename Impl::SitePropagator SitePropagator; \
typedef typename Impl::ScalarField ScalarField; \
typedef typename Impl::FermionField FermionField; \
typedef typename Impl::PropagatorField PropagatorField; \
typedef typename Impl::StencilImpl StencilImpl;
namespace Grid{
namespace QCD{
// Scalar implementation class ///////////////////////////////////////////////
// FIXME: it is not very nice to have the FImpl aliases
template <class S,
class Representation = FundamentalRep<1>,
class _Coeff_t = RealD>
class ScalarImpl:
public PeriodicGaugeImpl<GaugeImplTypes<S, Representation::Dimension>>
{
public:
static constexpr unsigned int rDim = Representation::Dimension;
public:
// gauge types
typedef PeriodicGaugeImpl<GaugeImplTypes<S, rDim>> Gimpl;
INHERIT_GIMPL_TYPES(Gimpl);
// site types
// (using classes instead of aliases to allow for partial specialisation)
template <typename vtype, unsigned int d>
class iImplScalar
{
public:
typedef iScalar<iScalar<iVector<vtype, d>>> type;
};
template <typename vtype>
class iImplScalar<vtype, 1>
{
public:
typedef iScalar<iScalar<iScalar<vtype>>> type;
};
template <typename vtype, unsigned int d>
class iImplPropagator
{
public:
typedef iScalar<iScalar<iMatrix<vtype, d>>> type;
};
template <typename vtype>
class iImplPropagator<vtype, 1>
{
public:
typedef iScalar<iScalar<iScalar<vtype>>> type;
};
// type aliases
typedef typename iImplScalar<S, rDim>::type SiteScalar;
typedef SiteScalar SiteSpinor;
typedef typename iImplPropagator<S, rDim>::type SitePropagator;
typedef Lattice<SiteScalar> ScalarField;
typedef ScalarField FermionField;
typedef Lattice<SitePropagator> PropagatorField;
typedef CartesianStencil<SiteScalar, SiteScalar> StencilImpl;
};
// single scalar implementation
typedef ScalarImpl<vComplex> ScalarImplR;
// Scalar action /////////////////////////////////////////////////////////////
template <typename SImpl>
class Scalar:
public CheckerBoardedSparseMatrixBase<typename SImpl::ScalarField>,
public SImpl
{
public:
INHERIT_GIMPL_TYPES(SImpl);
INHERIT_SIMPL_TYPES(SImpl);
public:
// constructor
Scalar(GaugeField &_Umu, GridCartesian &Sgrid, GridRedBlackCartesian &Hgrid,
RealD _mass)
: _grid(&Sgrid)
, _cbgrid(&Hgrid)
, mass(_mass)
, Lebesgue(_grid)
, LebesgueEvenOdd(_cbgrid)
, Umu(&Sgrid)
, UmuEven(&Hgrid)
, UmuOdd(&Hgrid)
{
Umu = _Umu;
pickCheckerboard(Even, UmuEven, Umu);
pickCheckerboard(Odd, UmuOdd, Umu);
}
// grid access
virtual GridBase *RedBlackGrid(void) {return _grid;}
// half checkerboard operations
// FIXME: do implementation
virtual void Meooe(const ScalarField &in, ScalarField &out)
{
assert(0);
}
virtual void Mooee(const ScalarField &in, ScalarField &out)
{
assert(0);
}
virtual void MooeeInv(const ScalarField &in, ScalarField &out)
{
assert(0);
}
virtual void MeooeDag(const ScalarField &in, ScalarField &out)
{
assert(0);
}
virtual void MooeeDag(const ScalarField &in, ScalarField &out)
{
assert(0);
}
virtual void MooeeInvDag(const ScalarField &in, ScalarField &out)
{
assert(0);
}
// free propagators
static void MomentumSpacePropagator(ScalarField &out, RealD m);
static void FreePropagator(const ScalarField &in, ScalarField &out,
const ScalarField &momKernel);
static void FreePropagator(const ScalarField &in, ScalarField &out, RealD m);
public:
RealD mass;
GridBase *_grid;
GridBase *_cbgrid;
// Defines the stencils for even and odd
StencilImpl Stencil;
StencilImpl StencilEven;
StencilImpl StencilOdd;
// Copy of the gauge field, with even and odd subsets
GaugeField Umu;
GaugeField UmuEven;
GaugeField UmuOdd;
LebesgueOrder Lebesgue;
LebesgueOrder LebesgueEvenOdd;
};
template <typename SImpl>
void Scalar<SImpl>::MomentumSpacePropagator(ScalarField &out, RealD m)
{
GridBase *grid = out._grid;
ScalarField kmu(grid);
const unsigned int nd = grid->_ndimension;
std::vector<int> &l = grid->_fdimensions;
out = m*m;
for(int mu = 0; mu < nd; mu++)
{
Real twoPiL = M_PI*2./l[mu];
LatticeCoordinate(kmu,mu);
kmu = 2.*sin(.5*twoPiL*kmu);
out = out + kmu*kmu;
}
}
template <typename SImpl>
void Scalar<SImpl>::FreePropagator(const ScalarField &in, ScalarField &out,
const ScalarField &FTKernel)
{
FFT fft((GridCartesian *)in._grid);
ScalarField inFT(in._grid);
fft.FFT_all_dim(inFT, in, FFT::forward);
inFT = inFT*FTKernel;
fft.FFT_all_dim(out, inFT, FFT::backward);
}
template <typename SImpl>
void Scalar<SImpl>::FreePropagator(const ScalarField &in, ScalarField &out,
RealD m)
{
ScalarField FTKernel(in._grid);
MomentumSpacePropagator(FTKernel, m);
FreePropagator(in, out, FTKernel);
}
template <class SImpl>
void ScalarToProp(typename SImpl::PropagatorField &p,
const typename SImpl::ScalarField &s,
const int c)
{
for(int i = 0; i < SImpl::rDim; ++i)
{
pokeColour(p, peekColour(s, i), i);
}
}
template <class SImpl>
void PropToScalar(typename SImpl::ScalarField &s,
const typename SImpl::PropagatorField &p,
const int c)
{
for(int i = 0; i < SImpl::rDim; ++i)
{
pokeColour(s, peekColour(p, i), i);
}
}
}}
#endif