/*
* main_plaquette.cpp
*
* Lukas Mazur, 10 Apr 2018 (original Author for Plaquette)
* Vishal Rao , 06 Jan 2024 (generalisation from plaquette to wilson loop)
*
* 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>
#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;
//Constructor to initialize all necessary members.
CalcWilson(Gaugefield<floatT,true,HaloDepth> &gauge, int &length, int &height) : SU3Accessor(gauge.getAccessor()), _length(length),_height(height){
}
__device__ __host__ SU3<floatT> HalfLoop(int &length, int &height,gSite site,int &mu, int &nu)
{
// 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;
typedef GIndexer<All, HaloDepth> GInd;
temp = SU3Accessor.getLink(GInd::getSiteMu(site, mu));
for (int len=1;len<length;len++){
site=GInd::site_up(site, mu);
temp *= SU3Accessor.getLink(GInd::getSiteMu(site, mu));
}
site=GInd::site_up(site,mu);
for (int hei=1;hei<=height;hei++){
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
__device__ __host__ floatT operator()(gSite site) {
/// We need to choose the type of indexer. The first template is the layout of the lattice.
typedef GIndexer<All, HaloDepth> GInd;
int l=_length;
int h=_height;
/// Define a SU(3) matrix
SU3<floatT> temp1,temp2;
floatT result = 0;
for (int nu = 1; nu < 3; nu++)
{
for (int mu = 0; mu < nu; mu++)
{
temp1=HalfLoop(l,h,site,mu,nu);
temp2=HalfLoop(h,l,site,nu,mu);
result+=tr_d(temp1*dagger(temp2));
/* temp1=HalfLoop(h,l,site,mu,nu); */
/* temp2=HalfLoop(l,h,site,nu,mu); */
/* result+=tr_d(temp1*dagger(temp2)); */
//Return the result
//The return value will be stored in the array of the reductionbase at index site.isite.
}
}
return result;
}
};
//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){
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));
//Do the final reduction
floatT Wloop;
redBase.reduce(Wloop, elems);
//Normalize the result
Wloop /= (GInd::getLatData().globalLattice().mult()*9); //3 loops(one averaged loop per plane) times 3 colors.
return Wloop;
}
//
// I am going to paste the hypSmearing.cpp file content here for the moment.
// also call updateAll()
template<class floatT, bool onDevice, size_t HaloDepth, CompressionType comp>
void HypSmearing<floatT, onDevice, HaloDepth, comp>::Su3Unitarize(Gaugefield<floatT, onDevice, HaloDepth, comp> &gauge_out, Gaugefield<floatT, onDevice, HaloDepth, comp> &gauge_base){
HypStaple<floatT, HaloDepth, comp, 4> su_3_unitarize(gauge_out.getAccessor(), gauge_base.getAccessor(), gauge_base.getAccessor(), gauge_base.getAccessor());
gauge_out.iterateOverBulkAllMu(su_3_unitarize);
if(update_all)gauge_out.updateAll();
}
template<class floatT, bool onDevice, size_t HaloDepth, CompressionType comp>
void HypSmearing<floatT, onDevice, HaloDepth, comp>::SmearAll(Gaugefield<floatT, onDevice, HaloDepth, comp> &gauge_out) {
// create level 1 fields
_dummy.iterateOverBulkAllMu(staple3_lvl1_10);
_gauge_lvl1_10 = (1-params.alpha_3) * _gauge_base + params.alpha_3/2 * _dummy;
Su3Unitarize(_gauge_lvl1_10, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl1_20);
_gauge_lvl1_20 = (1-params.alpha_3) * _gauge_base + params.alpha_3/2 * _dummy;
Su3Unitarize(_gauge_lvl1_20, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl1_30);
_gauge_lvl1_30 = (1-params.alpha_3) * _gauge_base + params.alpha_3/2 * _dummy;
Su3Unitarize(_gauge_lvl1_30, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl1_21);
_gauge_lvl1_21 = (1-params.alpha_3) * _gauge_base + params.alpha_3/2 * _dummy;
Su3Unitarize(_gauge_lvl1_21, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl1_31);
_gauge_lvl1_31 = (1-params.alpha_3) * _gauge_base + params.alpha_3/2 * _dummy;
Su3Unitarize(_gauge_lvl1_31, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl1_32);
_gauge_lvl1_32 = (1-params.alpha_3) * _gauge_base + params.alpha_3/2 * _dummy;
Su3Unitarize(_gauge_lvl1_32, _gauge_base);
// now that we have level 1 fields, create level 2 staples
// note: the order of the gauge fields goes in ascending order (10 < 20 < 30, 10 < 21 < 31, 20 < 21 < 32, 30 < 31 < 32)
// this is ASSUMED by HypStaple<floatT, HaloDepth, comp, 2>; DO NOT change this order without also modifying HypStaple<floatT, HaloDepth, comp, 2> and threeLinkStaple_second_level
HypStaple<floatT, HaloDepth, comp, 2> staple3_lvl2_0(_gauge_lvl1_10.getAccessor(), _gauge_lvl1_20.getAccessor(), _gauge_lvl1_30.getAccessor(), _dummy.getAccessor(), 0);
HypStaple<floatT, HaloDepth, comp, 2> staple3_lvl2_1(_gauge_lvl1_10.getAccessor(), _gauge_lvl1_21.getAccessor(), _gauge_lvl1_31.getAccessor(), _dummy.getAccessor(), 1);
HypStaple<floatT, HaloDepth, comp, 2> staple3_lvl2_2(_gauge_lvl1_20.getAccessor(), _gauge_lvl1_21.getAccessor(), _gauge_lvl1_32.getAccessor(), _dummy.getAccessor(), 2);
HypStaple<floatT, HaloDepth, comp, 2> staple3_lvl2_3(_gauge_lvl1_30.getAccessor(), _gauge_lvl1_31.getAccessor(), _gauge_lvl1_32.getAccessor(), _dummy.getAccessor(), 3);
//second level fields
_dummy.iterateOverBulkAllMu(staple3_lvl2_0);
_gauge_lvl2_0 = (1-params.alpha_2) * _gauge_base + params.alpha_2/4 * _dummy;
Su3Unitarize(_gauge_lvl2_0, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl2_1);
_gauge_lvl2_1 = (1-params.alpha_2) * _gauge_base + params.alpha_2/4 * _dummy;
Su3Unitarize(_gauge_lvl2_1, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl2_2);
_gauge_lvl2_2 = (1-params.alpha_2) * _gauge_base + params.alpha_2/4 * _dummy;
Su3Unitarize(_gauge_lvl2_2, _gauge_base);
_dummy.iterateOverBulkAllMu(staple3_lvl2_3);
_gauge_lvl2_3 = (1-params.alpha_2) * _gauge_base + params.alpha_2/4 * _dummy;
Su3Unitarize(_gauge_lvl2_3, _gauge_base);
// now that we have level 2 fields, create level 3 staple
HypStaple<floatT, HaloDepth, comp, 1> staple3_lvl3(_gauge_lvl2_0.getAccessor(), _gauge_lvl2_1.getAccessor(), _gauge_lvl2_2.getAccessor(), _gauge_lvl2_3.getAccessor());
_dummy.iterateOverBulkAllMu(staple3_lvl3);
// OLD VERSION, (MAYBE) DOES NOT WORK FOR SOME REASON
//_gauge_lvl2_0 = (1-params.alpha_1) * _gauge_base + params.alpha_1/6 * _dummy; //reused _gauge_lvl2_0
//Su3Unitarize(_gauge_lvl2_0);
// NEW VERSION, USES EXTRA FIELD RATHER THAN REUSE _gauge_lvl2_0
gauge_out = (1-params.alpha_1) * _gauge_base + params.alpha_1/6 * _dummy; //reused _gauge_lvl2_0
Su3Unitarize(gauge_out, _gauge_base);
}
#define CLASS_INIT(floatT,HALO) \
template class HypSmearing<floatT,true,HALO,R18>;
INIT_PH(CLASS_INIT)
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);
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
for(int no_smearing_steps=0;no_smearing_steps<=100;no_smearing_steps+=5)
{
//we shall make a variable which stores gauge configuration
//doing smearing 5,10,15.... times, (loop will run just 5 times where we left of using gauge=gauge_out)
if (no_smearing_steps!=0){
for(int j=1;j<=5;j++)
{
HypSmearing<floatT, true, HaloDepth, R18> hypsmearing(gauge); //smearing of gauge
hypsmearing.SmearAll(gauge_out); //and storing at gauge_out
gauge=gauge_out;
}
}
const int Ns= LatDim[0];
for (int length=1; length<Ns; 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, redBase, length,height );
printf("Node:%d,Smear_count:%d,Length=%d,Height=%d,Wilson Loop:%1.15e\n",node,no_smearing_steps,length, height,Wloop);
/* rootLogger.info("Reduced Rectangle from rhmc : " , gaugeaction.rectangle()); */
/* rootLogger.info("Reduced plaquette from rhmc : " , gaugeaction.plaquette()); */
}
}
}
}
return 0;
}
//