/******************************************************************************
* File: main_plaquette.cpp
*
* Author: Vishal Rao
* Created: 02/08/24
* Description: This program contains the function to calculate the wilson loop with hyp smearing and generation of more points using bresenham algorithm described in :https://arxiv.org/pdf/hep-lat/0005018.pdf
*****************************************************************************
* main_plaquette.cpp
* This is just an example how a very very basic program works. Look at src/testing/main_GeneralOperatorTest.cpp to see
* how to write more advanced GPU code.
*
*/
#include "../simulateqcd.h"
/* #include "../modules/hyp/hypParameters.h" was already included in ../modules/hypSmearing.h */
/* #include "../modules/hyp/smearParameters.h" */
#include "../modules/hyp/hypSmearing.h"
#include "../modules/hyp/hypSmearing.cpp"
#include <algorithm>
#include <iostream>
#include <string>
#include <vector>
#include <cmath>
#define PREC double
template<class floatT,size_t HaloDepth>
struct CalcWilson{
//Gauge accessor to access the gauge field
SU3Accessor<floatT> SU3Accessor;
int _length=1;
int _height=1;
int _nu_in=0; // nu_in is the direction of fictitious time , in this direction we don't do smearing.
//Constructor to initialize all necessary members.
CalcWilson(Gaugefield<floatT,true,HaloDepth> &gauge, int &length, int &height, int &nu_in) : SU3Accessor(gauge.getAccessor()), _length(length),_height(height), _nu_in(nu_in){
}
__device__ __host__ SU3<floatT> HalfLoop(int &length, int &height,gSite site,int &mu, int &nu)
{
// nu_in is the direction of fictitious time , in this direction we don't do smearing. // we have not taken site variable as reference as if it were take as reference, it will update the values at same mEmory location
SU3<floatT> temp;
temp=SU3<floatT>(1, 0, 0, 0, 1, 0, 0, 0, 1); //initialize to identity matrix
typedef GIndexer<All, HaloDepth> GInd;
for (int il=0;il<length;il++){
temp *= SU3Accessor.getLink(GInd::getSiteMu(site, mu));
site=GInd::site_up(site, mu);
}
for (int ih=0;ih<height;ih++){
temp *= SU3Accessor.getLink(GInd::getSiteMu(site, nu));
site=GInd::site_up(site,nu);
}
return temp;
}
//This is the operator that is called inside the Kernel // nu_in is the direction of fictitious time , in this direction we don't do smearing.
__device__ __host__ floatT operator()(gSite site) {
/// We need to choose the type of indexer. The first template is the layout of the lattice. // nu_in is the direction of fictitious time , in this direction we don't do smearing.
typedef GIndexer<All, HaloDepth> GInd;
int l=_length;
int h=_height;
/// Define a SU(3) matrix
SU3<floatT> temp1,temp2;
floatT result = 0;
int mu1,nu1,mu2,nu2;
switch (_nu_in) {
case 0: {nu1=0;nu2=0;mu1=1;mu2=2;}
break;
case 1: {nu1=1;nu2=1;mu1=0;mu2=2;}
break;
case 2: {nu1=2;nu2=2;mu1=0;mu2=1;}
break;
}
temp1=HalfLoop(l,h,site,mu1,nu1);
temp2=HalfLoop(h,l,site,nu1,mu1);
result+=tr_d(temp2*dagger(temp1));
temp1=HalfLoop(l,h,site,mu2,nu2);
temp2=HalfLoop(h,l,site,nu2,mu2);
result+=tr_d(temp2*dagger(temp1));
return result/2.0;
}
};
template<class floatT,size_t HaloDepth>
struct CalcWilsonBresenham{
int _x,_y,_z, _ex_dir;
SU3Accessor<floatT> SU3Accessor;
CalcWilsonBresenham(Gaugefield<floatT,true,HaloDepth> &gauge, int &x,int &y, int &z, int &ex_dir) : SU3Accessor(gauge.getAccessor()), _x(x),_y(y),_z(z), _ex_dir(ex_dir){}
//This is the operator that is called inside the Kernel
__device__ __host__ SU3<floatT> BresenhamLoopXYL(int &x,int &y, int &z, gSite site, int &mu, int &nu)
{
SU3<floatT> temp;
temp=SU3<floatT>(1, 0, 0, 0, 1, 0, 0, 0, 1); //initialize to identity matrix
int cmax,cmin,cmax_dir,cmin_dir;
if (x>=y) {cmax=x; cmin=y; cmax_dir=mu; cmin_dir=nu;}
else {cmax=y; cmin=x; cmax_dir=nu; cmin_dir=mu;}
typedef GIndexer<All, HaloDepth> GInd;
int cmax2=2*cmax;
int cmin2=2*cmin;
int chi=cmin2-cmax;
for (int ix=0;ix<cmax;ix++)
{
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,cmax_dir));//step in cmax
site=GInd::site_up(site,cmax_dir);
if(chi>=0)
{
chi-=cmax2;
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,cmin_dir)); //step in cmin
site=GInd::site_up(site,cmin_dir);
}
chi+=cmin2;
}
for(int iz=0;iz<z;iz++)
{
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,3-mu-nu)); //step in cmin
site=GInd::site_up(site,3-mu-nu);
}
return temp;
}
__device__ __host__ SU3<floatT> BresenhamLoopXYU(int &x,int &y, int &z, gSite site, int &mu, int &nu)
{
SU3<floatT> temp;
temp=SU3<floatT>(1, 0, 0, 0, 1, 0, 0, 0, 1); //initialize to identity matrix
int cmax,cmin,cmax_dir,cmin_dir;
if (x>=y) {cmax=x; cmin=y; cmax_dir=mu; cmin_dir=nu;}
else {cmax=y; cmin=x; cmax_dir=nu; cmin_dir=mu;}
typedef GIndexer<All, HaloDepth> GInd;
int cmax2=2*cmax;
int cmin2=2*cmin;
int chi=cmin2-cmax;
for(int iz=0;iz<z;iz++)
{
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,3-mu-nu)); //step in cmin
site=GInd::site_up(site,3-mu-nu);
}
for (int ix=0;ix<cmax;ix++)
{
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,cmax_dir));//step in cmax
site=GInd::site_up(site,cmax_dir);
if(chi>=0)
{
chi-=cmax2;
temp*=SU3Accessor.getLink(GInd::getSiteMu(site,cmin_dir)); //step in cmin
site=GInd::site_up(site,cmin_dir);
}
chi+=cmin2;
}
return temp;
}
__device__ __host__ floatT operator()(gSite site) {
typedef GIndexer<All, HaloDepth> GInd;
/// Define a SU(3) matrix
floatT result = 0;
SU3<floatT> temp1,temp2;
int mu,nu;
switch (_ex_dir) {
case 0: mu=1;nu=2;break;
case 1: mu=0;nu=2;break;
case 2: mu=0;nu=1;break;
}
temp1=BresenhamLoopXYL(_x,_y,_z,site,mu,nu);
temp2=BresenhamLoopXYU(_x,_y,_z,site,mu,nu);
result+=tr_d(temp1*dagger(temp2));
//replacing x<->y
temp1=BresenhamLoopXYL(_y,_x,_z,site,mu,nu);
temp2=BresenhamLoopXYU(_y,_x,_z,site,mu,nu);
result+=tr_d(temp1*dagger(temp2));
return result/2.0;
}
};
//Function to compute the wilson loop using the above struct CalcWilson.
template<class floatT, size_t HaloDepth>
floatT WilsonLoop(Gaugefield<floatT,true, HaloDepth> &gauge, LatticeContainer<true,floatT> &redBase, int &l, int &h, int &ex_dir){
typedef GIndexer<All,HaloDepth> GInd;
const size_t elems = GInd::getLatData().vol4;
//Make sure, redBase is large enough
redBase.adjustSize(elems);
// we will iterate the process of finding wilson loop over whole bulk.
redBase.template iterateOverBulk<All, HaloDepth>(CalcWilson<floatT, HaloDepth>(gauge, l,h,ex_dir)); //ex_dir is excluded direction for smearing
//Do the final reduction
floatT Wloop;
redBase.reduce(Wloop, elems);
//Normalize the result
Wloop /= (GInd::getLatData().globalLattice().mult()*3); //3 colors
return Wloop;
}
//defining function for bresenham
template<class floatT, size_t HaloDepth>
floatT WilsonLoopBresenham(Gaugefield<floatT,true, HaloDepth> &gauge, LatticeContainer<true,floatT> &redBase, int &x, int &y,int &z, int &ex_dir){
typedef GIndexer<All,HaloDepth> GInd;
const size_t elems = GInd::getLatData().vol4;
//Make sure, redBase is large enough
redBase.adjustSize(elems);
// we will iterate the process of finding wilson loop over whole bulk.
redBase.template iterateOverBulk<All, HaloDepth>(CalcWilsonBresenham<floatT, HaloDepth>(gauge,x,y,z, ex_dir));
//Do the final reduction
floatT Wloop;
redBase.reduce(Wloop, elems);
//Normalize the result
const int n_colors=3;
Wloop /= (GInd::getLatData().globalLattice().mult()*n_colors); //3 loops(one averaged loop per plane) times 3 colors=9.
return Wloop;
}
int main(int argc, char *argv[]) {
stdLogger.setVerbosity(DEBUG);
/// Initialize parameter class. This class can also read parameter from textfiles!
LatticeParameters param;
/// Initialize the Lattice dimension
const int LatDim[] = {32, 32, 32, 8};
const int NodeDim[] = {1, 1, 1, 1};
/// Just pass these dimensions to the parameter class
param.latDim.set(LatDim);
param.nodeDim.set(NodeDim);
/// Initialize a timer
StopWatch<true> timer;
/// Initialize the CommunicationBase. This class handles the communitation between different Cores/GPU's.
CommunicationBase commBase(&argc, &argv, true);
commBase.init(param.nodeDim());
const size_t HaloDepth = 1; //since it is not multi-gpu code
/// highlight the output differently.
rootLogger.info("Initialize Lattice");
/// Initialize the Indexer on GPU and CPU.
initIndexer(HaloDepth,param,commBase);
typedef GIndexer<All,HaloDepth> GInd;
rootLogger.info("Initialize Gaugefield");
Gaugefield<PREC, true,HaloDepth> gauge(commBase);
typedef double floatT;
GaugeAction<floatT, true, HaloDepth, R18> gaugeaction(gauge);
/// Initialize gaugefield with unity-matrices.
gauge.one();
/// Initialize LatticeContainer. This is in principle the "array", where the values of the plaquette are
/// stored which are summed up in the end
LatticeContainer<true,PREC> redBase(commBase);
/// We need to tell the Reductionbase how large our Array will be
redBase.adjustSize(GInd::getLatData().vol4);
int eqm_point=155;
int last_point=1150;
int node;
sscanf(argv[1], "%d", &node);
//sscanf(argv[2], "%d", &height);
int no_smearing_steps=1;
for (int i=eqm_point; i<=last_point;i+=5)
{
rootLogger.info("Read configuration");
gauge.readconf_nersc("/root/project1/build_SIMULATeQCD/applications/try_20_output_dir/after_eqm/node"+std::to_string(node)+"/l328f21b6285m0039185m0783706a_"+std::to_string(node)+"."+std::to_string(i));
gauge.updateAll();
Gaugefield<PREC, true,HaloDepth> gauge_out(commBase);
gauge_out.one(); //initialize the gauge field variable
Gaugefield<PREC, true,HaloDepth> gauge_copy(commBase);
gauge_copy.one(); //initialize the gauge field variable
gauge_copy=gauge;
for(int ex_dir=0;ex_dir<3;ex_dir++)
{
for(int j=1;j<=no_smearing_steps;j++)
{
/* HypSmearing<floatT, true, HaloDepth, R18> hypsmearing(gauge); //smearing of gauge */
/* hypsmearing.SmearAll(gauge_out); //and storing at gauge_out */
/* hypsmearing.SmearAll_spatial(gauge_out, ex_dir); //and storing at gauge_out */
HypSmearing_spatial<floatT, true, HaloDepth,R18> HypSmearing_spatial(gauge);
HypSmearing_spatial.SmearAll_spatial(gauge_out);
gauge=gauge_out;
}
gauge=gauge_copy;
const int Ns= LatDim[0];
for (int length=1; length<=1; length++)
{
for (int height=1; height<Ns; height++)
{
PREC Wloop = 0;
/* timer.start(); */
/// compute wilson loop with smeared gauge_out
Wloop = WilsonLoop<PREC,HaloDepth>(gauge_out, redBase, length,height, ex_dir );
printf("Node:%d,Smear_count:%d,Length=%d,Height=%d,ex_dir=%d,Wilson Loop:%1.15e\n",node,no_smearing_steps,length, height,ex_dir,Wloop);
/* rootLogger.info("Reduced Rectangle from rhmc : " , gaugeaction.rectangle()); */
/* rootLogger.info("Reduced plaquette from rhmc : " , gaugeaction.plaquette()); */
}
}
// we shall now implement bresenham algorithm to calculate the new wilson loops.int x_coord=length;
/* for (int a=1;a<Ns;a++) */
/* { */
/* for(int b=a;b<Ns;b++) */
/* { */
/* for(int c=1;c<Ns;c++) */
/* { */
/* PREC Wloop=WilsonLoopBresenham<PREC,HaloDepth>(gauge_out, redBase,a,b,c, ex_dir ); */
/* printf("Node:%d,Smear_count:%d,Length=%f,Height=%d,ex_dir=%d,Wilson Loop:%1.15e\n",node,no_smearing_steps,pow((a*a+b*b),0.5), c,ex_dir,Wloop); */
/* } */
/* } */
/* } */
}
}
return 0;
}