mirror of
https://github.com/paboyle/Grid.git
synced 2024-11-09 23:45:36 +00:00
Extra examples / solutions
This commit is contained in:
Peter Boyle
parent
40098424c7
commit
6135ad530e
122
examples/Example_Laplacian_inverter.cc
Normal file
122
examples/Example_Laplacian_inverter.cc
Normal file
@ -0,0 +1,122 @@
|
||||
#include <Grid/Grid.h>
|
||||
using namespace Grid;
|
||||
|
||||
// Function used for Chebyshev smearing
|
||||
//
|
||||
Real MomentumSmearing(Real p2)
|
||||
{
|
||||
return (1 - 4.0*p2) * exp(-p2/4);
|
||||
}
|
||||
Real DistillationSmearing(Real p2)
|
||||
{
|
||||
if ( p2 > 0.5 ) return 0.0;
|
||||
else return 1.0;
|
||||
}
|
||||
|
||||
// Flip sign to make prop to p^2, not -p^2 relative to last example
|
||||
template<class Gimpl,class Field> class CovariantLaplacianCshift : public SparseMatrixBase<Field>
|
||||
{
|
||||
public:
|
||||
INHERIT_GIMPL_TYPES(Gimpl);
|
||||
|
||||
GridBase *grid;
|
||||
GaugeField U;
|
||||
|
||||
CovariantLaplacianCshift(GaugeField &_U) :
|
||||
grid(_U.Grid()),
|
||||
U(_U) { };
|
||||
|
||||
virtual GridBase *Grid(void) { return grid; };
|
||||
|
||||
virtual void M (const Field &in, Field &out)
|
||||
{
|
||||
out=Zero();
|
||||
for(int mu=0;mu<Nd-1;mu++) {
|
||||
GaugeLinkField Umu = PeekIndex<LorentzIndex>(U, mu); // NB: Inefficent
|
||||
out = out - Gimpl::CovShiftForward(Umu,mu,in);
|
||||
out = out - Gimpl::CovShiftBackward(Umu,mu,in);
|
||||
out = out + 2.0*in;
|
||||
}
|
||||
};
|
||||
virtual void Mdag (const Field &in, Field &out) { M(in,out);}; // Laplacian is hermitian
|
||||
virtual void Mdiag (const Field &in, Field &out) {assert(0);}; // Unimplemented need only for multigrid
|
||||
virtual void Mdir (const Field &in, Field &out,int dir, int disp){assert(0);}; // Unimplemented need only for multigrid
|
||||
virtual void MdirAll (const Field &in, std::vector<Field> &out) {assert(0);}; // Unimplemented need only for multigrid
|
||||
};
|
||||
|
||||
|
||||
|
||||
int main(int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc, &argv);
|
||||
|
||||
typedef LatticeColourVector Field;
|
||||
|
||||
auto latt_size = GridDefaultLatt();
|
||||
auto simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
|
||||
auto mpi_layout = GridDefaultMpi();
|
||||
|
||||
GridCartesian Grid(latt_size,simd_layout,mpi_layout);
|
||||
GridParallelRNG RNG(&Grid); RNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
|
||||
|
||||
|
||||
LatticeGaugeField U(&Grid);
|
||||
|
||||
SU<Nc>::ColdConfiguration(RNG,U);
|
||||
|
||||
typedef CovariantLaplacianCshift <PeriodicGimplR,Field> Laplacian_t;
|
||||
Laplacian_t Laplacian(U);
|
||||
|
||||
ColourVector ColourKronecker;
|
||||
ColourKronecker = Zero();
|
||||
ColourKronecker()()(0) = 1.0;
|
||||
|
||||
Coordinate site({latt_size[0]/2,
|
||||
latt_size[1]/2,
|
||||
latt_size[2]/2,
|
||||
0});
|
||||
|
||||
Field kronecker(&Grid);
|
||||
kronecker = Zero();
|
||||
pokeSite(ColourKronecker,kronecker,site);
|
||||
|
||||
|
||||
Field psi(&Grid), chi(&Grid);
|
||||
|
||||
//////////////////////////////////////
|
||||
// Classic Wuppertal smearing
|
||||
//////////////////////////////////////
|
||||
|
||||
Integer Iterations = 80;
|
||||
Real width = 2.0;
|
||||
Real coeff = (width*width) / Real(4*Iterations);
|
||||
|
||||
chi=kronecker;
|
||||
// chi = (1-p^2/2N)^N kronecker
|
||||
for(int n = 0; n < Iterations; ++n) {
|
||||
Laplacian.M(chi,psi);
|
||||
chi = chi - coeff*psi;
|
||||
}
|
||||
|
||||
std::cout << " Wuppertal smeared operator is chi = \n" << chi <<std::endl;
|
||||
|
||||
/////////////////////////////////////
|
||||
// Chebyshev smearing
|
||||
/////////////////////////////////////
|
||||
RealD lo = 0.0;
|
||||
RealD hi = 128;
|
||||
HermitianLinearOperator<Laplacian_t,Field> HermOp(Laplacian);
|
||||
|
||||
Chebyshev<Field> ChebySmear(lo,hi,20,DistillationSmearing);
|
||||
// Chebyshev<Field> ChebySmear(lo,hi,20,MomentumSmearing);
|
||||
{
|
||||
std::ofstream of("chebysmear");
|
||||
ChebySmear.csv(of);
|
||||
}
|
||||
|
||||
ChebySmear(HermOp,kronecker,chi);
|
||||
|
||||
std::cout << " Chebyshev smeared operator is chi = \n" << chi <<std::endl;
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
127
examples/Example_Laplacian_smearing.cc
Normal file
127
examples/Example_Laplacian_smearing.cc
Normal file
@ -0,0 +1,127 @@
|
||||
#include <Grid/Grid.h>
|
||||
using namespace Grid;
|
||||
|
||||
// Function used for Chebyshev smearing
|
||||
//
|
||||
Real MomentumSmearing(Real p2)
|
||||
{
|
||||
return (1 - 4.0*p2) * exp(-p2/4);
|
||||
}
|
||||
Real DistillationSmearing(Real p2)
|
||||
{
|
||||
if ( p2 > 0.5 ) return 0.0;
|
||||
else return 1.0;
|
||||
}
|
||||
|
||||
// Flip sign to make prop to p^2, not -p^2 relative to last example
|
||||
template<class Gimpl,class Field> class CovariantLaplacianCshift : public SparseMatrixBase<Field>
|
||||
{
|
||||
public:
|
||||
INHERIT_GIMPL_TYPES(Gimpl);
|
||||
|
||||
GridBase *grid;
|
||||
GaugeField U;
|
||||
|
||||
CovariantLaplacianCshift(GaugeField &_U) :
|
||||
grid(_U.Grid()),
|
||||
U(_U) { };
|
||||
|
||||
virtual GridBase *Grid(void) { return grid; };
|
||||
|
||||
virtual void M (const Field &in, Field &out)
|
||||
{
|
||||
out=Zero();
|
||||
for(int mu=0;mu<Nd-1;mu++) {
|
||||
GaugeLinkField Umu = PeekIndex<LorentzIndex>(U, mu); // NB: Inefficent
|
||||
out = out - Gimpl::CovShiftForward(Umu,mu,in);
|
||||
out = out - Gimpl::CovShiftBackward(Umu,mu,in);
|
||||
out = out + 2.0*in;
|
||||
}
|
||||
};
|
||||
virtual void Mdag (const Field &in, Field &out) { M(in,out);}; // Laplacian is hermitian
|
||||
virtual void Mdiag (const Field &in, Field &out) {assert(0);}; // Unimplemented need only for multigrid
|
||||
virtual void Mdir (const Field &in, Field &out,int dir, int disp){assert(0);}; // Unimplemented need only for multigrid
|
||||
virtual void MdirAll (const Field &in, std::vector<Field> &out) {assert(0);}; // Unimplemented need only for multigrid
|
||||
};
|
||||
|
||||
|
||||
|
||||
int main(int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc, &argv);
|
||||
|
||||
typedef LatticeColourVector Field;
|
||||
|
||||
auto latt_size = GridDefaultLatt();
|
||||
auto simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
|
||||
auto mpi_layout = GridDefaultMpi();
|
||||
|
||||
GridCartesian Grid(latt_size,simd_layout,mpi_layout);
|
||||
GridParallelRNG RNG(&Grid); RNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
|
||||
|
||||
|
||||
LatticeGaugeField U(&Grid);
|
||||
|
||||
SU<Nc>::ColdConfiguration(RNG,U);
|
||||
|
||||
typedef CovariantLaplacianCshift <PeriodicGimplR,Field> Laplacian_t;
|
||||
Laplacian_t Laplacian(U);
|
||||
|
||||
|
||||
ColourVector ColourKronecker;
|
||||
ColourKronecker = Zero();
|
||||
ColourKronecker()()(0) = 1.0;
|
||||
|
||||
Coordinate site({latt_size[0]/2,
|
||||
latt_size[1]/2,
|
||||
latt_size[2]/2,
|
||||
0});
|
||||
|
||||
Field kronecker(&Grid);
|
||||
kronecker = Zero();
|
||||
pokeSite(ColourKronecker,kronecker,site);
|
||||
|
||||
|
||||
Field psi(&Grid), chi(&Grid);
|
||||
|
||||
//////////////////////////////////////
|
||||
// Classic Wuppertal smearing
|
||||
//////////////////////////////////////
|
||||
|
||||
Integer Iterations = 80;
|
||||
Real width = 2.0;
|
||||
Real coeff = (width*width) / Real(4*Iterations);
|
||||
|
||||
chi=kronecker;
|
||||
// chi = (1-p^2/2N)^N kronecker
|
||||
for(int n = 0; n < Iterations; ++n) {
|
||||
Laplacian.M(chi,psi);
|
||||
chi = chi - coeff*psi;
|
||||
}
|
||||
|
||||
std::cout << " Wuppertal smeared operator is chi = \n" << chi <<std::endl;
|
||||
|
||||
/////////////////////////////////////
|
||||
// Chebyshev smearing
|
||||
/////////////////////////////////////
|
||||
RealD lo = 0.0;
|
||||
RealD hi = 12.0; // Analytic free field bound
|
||||
HermitianLinearOperator<Laplacian_t,Field> HermOp(Laplacian);
|
||||
|
||||
std::cout << " Checking spectral range of our POSITIVE definite operator \n";
|
||||
PowerMethod<Field> PM;
|
||||
PM(HermOp,kronecker);
|
||||
|
||||
Chebyshev<Field> ChebySmear(lo,hi,20,DistillationSmearing);
|
||||
// Chebyshev<Field> ChebySmear(lo,hi,20,MomentumSmearing);
|
||||
{
|
||||
std::ofstream of("chebysmear");
|
||||
ChebySmear.csv(of);
|
||||
}
|
||||
|
||||
ChebySmear(HermOp,kronecker,chi);
|
||||
|
||||
std::cout << " Chebyshev smeared operator is chi = \n" << chi <<std::endl;
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
127
examples/Example_Laplacian_solver.cc
Normal file
127
examples/Example_Laplacian_solver.cc
Normal file
@ -0,0 +1,127 @@
|
||||
#include <Grid/Grid.h>
|
||||
using namespace Grid;
|
||||
|
||||
template<class Field>
|
||||
void SimpleConjugateGradient(LinearOperatorBase<Field> &HPDop,const Field &b, Field &x)
|
||||
{
|
||||
RealD cp, c, alpha, d, beta, ssq, qq;
|
||||
RealD Tolerance=1.0e-10;
|
||||
int MaxIterations=10000;
|
||||
|
||||
Field p(b), mmp(b), r(b);
|
||||
|
||||
HPDop.HermOpAndNorm(x, mmp, d, beta);
|
||||
|
||||
r = b - mmp;
|
||||
p = r;
|
||||
|
||||
alpha = norm2(p);
|
||||
cp = alpha;
|
||||
ssq = norm2(b);
|
||||
|
||||
RealD rsq = Tolerance * Tolerance * ssq;
|
||||
|
||||
for (int k = 1; k <= MaxIterations; k++) {
|
||||
c = cp;
|
||||
|
||||
HPDop.HermOp(p, mmp);
|
||||
|
||||
ComplexD dc = innerProduct(p,mmp);
|
||||
d = dc.real();
|
||||
alpha = c / d;
|
||||
|
||||
cp = axpy_norm(r, -alpha, mmp, r);
|
||||
beta = cp / c;
|
||||
|
||||
x = x + alpha* p ;
|
||||
p = r + beta* p ;
|
||||
|
||||
std::cout << "iteration "<<k<<" cp " <<std::sqrt(cp/ssq) << std::endl;
|
||||
if (cp <= rsq) {
|
||||
return;
|
||||
}
|
||||
}
|
||||
assert(0);
|
||||
}
|
||||
|
||||
|
||||
|
||||
// Flip sign to make prop to p^2, not -p^2 relative to last example
|
||||
template<class Gimpl,class Field> class CovariantLaplacianCshift : public SparseMatrixBase<Field>
|
||||
{
|
||||
public:
|
||||
INHERIT_GIMPL_TYPES(Gimpl);
|
||||
|
||||
GridBase *grid;
|
||||
GaugeField U;
|
||||
RealD m2=1.0e-2;
|
||||
CovariantLaplacianCshift(GaugeField &_U) :
|
||||
grid(_U.Grid()),
|
||||
U(_U) { };
|
||||
|
||||
virtual GridBase *Grid(void) { return grid; };
|
||||
|
||||
virtual void M (const Field &in, Field &out)
|
||||
{
|
||||
out=Zero();
|
||||
for(int mu=0;mu<Nd-1;mu++) {
|
||||
GaugeLinkField Umu = PeekIndex<LorentzIndex>(U, mu); // NB: Inefficent
|
||||
out = out - Gimpl::CovShiftForward(Umu,mu,in);
|
||||
out = out - Gimpl::CovShiftBackward(Umu,mu,in);
|
||||
out = out + 2.0*in + m2*in;
|
||||
}
|
||||
};
|
||||
virtual void Mdag (const Field &in, Field &out) { M(in,out);}; // Laplacian is hermitian
|
||||
virtual void Mdiag (const Field &in, Field &out) {assert(0);}; // Unimplemented need only for multigrid
|
||||
virtual void Mdir (const Field &in, Field &out,int dir, int disp){assert(0);}; // Unimplemented need only for multigrid
|
||||
virtual void MdirAll (const Field &in, std::vector<Field> &out) {assert(0);}; // Unimplemented need only for multigrid
|
||||
};
|
||||
|
||||
|
||||
|
||||
int main(int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc, &argv);
|
||||
|
||||
typedef LatticeColourVector Field;
|
||||
|
||||
auto latt_size = GridDefaultLatt();
|
||||
auto simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
|
||||
auto mpi_layout = GridDefaultMpi();
|
||||
|
||||
GridCartesian Grid(latt_size,simd_layout,mpi_layout);
|
||||
GridParallelRNG RNG(&Grid); RNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
|
||||
|
||||
|
||||
LatticeGaugeField U(&Grid);
|
||||
|
||||
SU<Nc>::ColdConfiguration(RNG,U);
|
||||
|
||||
typedef CovariantLaplacianCshift <PeriodicGimplR,Field> Laplacian_t;
|
||||
Laplacian_t Laplacian(U);
|
||||
|
||||
|
||||
ColourVector ColourKronecker;
|
||||
ColourKronecker = Zero();
|
||||
ColourKronecker()()(0) = 1.0;
|
||||
|
||||
Coordinate site({0,0,0,0}); // Point source at origin
|
||||
|
||||
Field kronecker(&Grid);
|
||||
kronecker = Zero();
|
||||
pokeSite(ColourKronecker,kronecker,site);
|
||||
|
||||
Field psi(&Grid); psi=Zero();
|
||||
|
||||
HermitianLinearOperator<Laplacian_t,Field> HermOp(Laplacian);
|
||||
SimpleConjugateGradient(HermOp, kronecker,psi);
|
||||
|
||||
Field r(&Grid);
|
||||
Laplacian.M(psi,r);
|
||||
r=kronecker-r;
|
||||
|
||||
std::cout << "True residual "<< norm2(r) <<std::endl;
|
||||
std::cout << psi<<std::endl;
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user