/*****************************************************************************
* 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/gaugeFixing/gfix.h"
#include "../modules/observables/polyakovLoop.h"
#include "../modules/observables/polyakovLoop.cpp"
#include "../modules/hyp/hypSmearing.cpp"
#include <algorithm>
#include <iostream>
#include <string>
#include <vector>
#include <cmath>
#define PREC double
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[] = {40, 40, 40, 10}; */
const int LatDim[] = {32, 32, 32, 8};
/* const int LatDim[] = {32, 32, 32, 32}; */
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 );
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);
rootLogger.info("Initialize Gaugefield");
Gaugefield<PREC, true,HaloDepth,R18> gauge(commBase);
/// Initialize gaugefield with unity-matrices.
gauge.one();
typedef double floatT;
int node;
sscanf(argv[1], "%d", &node);
int eqm_point=1;
int jump=1;
int n_gauge_conf_per_node=28;
int x_max=21,y_max=21,z_max=21;
/* bool avg_dir_results=true; //if true, then the average of the results in all directions is taken */
/* bool spatial_wloop=true; //if true, then the spatial wilson loop is calculated */
int no_smearing_steps=4;
int start_conf,end_conf;
int break_conf_no=1; //if we stop the calculation and want to resume; it will begun with reading the configuration number written here.; 1 is default if starting from 1st configuration again
start_conf=eqm_point+(node*n_gauge_conf_per_node*jump)+break_conf_no-1;
end_conf= start_conf +1 -break_conf_no + n_gauge_conf_per_node*jump;
for (int i=start_conf; i<end_conf;i+=jump)
{
rootLogger.info("Read configuration");
try {
gauge.readconf_milc("/root/project1/build_SIMULATeQCD/work/petreczk/l328f21b8400m0004875m009750/conf."+std::to_string(i));
/* gauge.readconf_nersc("/root/project1/build_SIMULATeQCD/applications/try_4_output_dir/after_eqm/links/conf."+std::to_string(i)); */
}
catch(const std::runtime_error& e) {
continue;
}
gauge.updateAll();
Gaugefield<PREC, true,HaloDepth,R18> gauge_out(commBase);
gauge_out.one(); //initialize the gauge field variable
Gaugefield<PREC, true,HaloDepth,R18> gauge_original(commBase);
gauge_original.one(); //initialize the gauge field variable
gauge_original=gauge;
GaugeAction<floatT, true, HaloDepth, R18> gaugeaction(gauge);
printf("plaquette=%.16f\n",gaugeaction.plaquette());
printf("plaquette_ss=%.16f\n",gaugeaction.plaquetteSS());
printf("plaquette_st=%.16f\n",gaugeaction.plaquetteST());
printf("rectangle=%.16f\n",gaugeaction.rectangle());
printf("gauge action density (symanzik_staggered) =%.16f\n",gaugeaction.symanzik_staggered()/(LatDim[0]*LatDim[1]*LatDim[2]*LatDim[3]));
printf("gauge action density =%.16f\n",gaugeaction.symanzik()/(LatDim[0]*LatDim[1]*LatDim[2]*LatDim[3]));
PolyakovLoop<floatT, true, HaloDepth,R18> ploop(gauge);
rootLogger.info("Smear Steps = 0,Polyakov Loop = ",ploop.getPolyakovLoop());
/* for(int ex_dir=0;ex_dir<3;ex_dir++) */
/* { */
/* rootLogger.info("Polyakov Loop (in ex_dir) in before smearing (Real) =",gaugeaction.Planer_Wloop(0,LatDim[ex_dir],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,-1); //and storing at gauge_out
/* HypSmearing_spatial<floatT, true, HaloDepth,R18> HypSmearing(gauge); */
/* HypSmearing.SmearAll_spatial(gauge_out,ex_dir); */
gauge=gauge_out;
/* for(int ex_dir=0;ex_dir<3;ex_dir++) */
/* { */
/* GaugeAction<floatT, true, HaloDepth, R18> gaugeaction1(gauge_out); */
/* rootLogger.info("Polyakov Loop (in ex_dir) in after smearing (Real) =",gaugeaction.Planer_Wloop(0,LatDim[ex_dir]),ex_dir); */
rootLogger.info("Smear Steps = ",j,",Polyakov Loop = ",ploop.getPolyakovLoop());
}
/* TwoGaugeAction<floatT, true, HaloDepth,R18> Twogaugeaction(gauge_out,gauge_original); */
/* int sum=Twogaugeaction.HypSmearExDir_Checker(ex_dir); */
/* printf("Number of links smeared in ex_dir(mu=%d) = %d \n",ex_dir,sum); */
/* printf("finding 2d bresenham results \n"); */
for(int ex_dir=0;ex_dir<3;ex_dir++)
{
for (int a=0;a<x_max;a++)
{
for(int b=a;b<y_max;b++)
{
for(int c=1;c<z_max;c++)
{ if (a==0 and b==0) {continue;} //skip the case where a=0, b=0 and c=1
PREC Wloop=gaugeaction.wLoopBresenham(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);
}
}
}
/* gauge=gauge_original; //reset the gauge field to original after calculating the bresenham results */
}
/* printf(" 2d bresenham results end here \n"); */
/* int &t_dir=ex_dir; */
// calculating planer wilson loops in all directions and averaging the results
/* for (int length=1; length<x_max; length++) */
/* { */
/* for (int height=1; height<t_max; height++) */
/* { */
/* PREC Wloop = gaugeaction.Planer_Wloop(length,height,t_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,t_dir,Wloop); */
/* } */
/* } */
/* int t_plus=spatial_wloop ? t_max:1; */
/* for (int x=0;x<x_max;x+=1) */
/* { */
/* int y_min = avg_dir_results ? x : 0; */
/* for ( int y=y_min; y<y_max;y+=1) */
/* { */
/* int z_min = avg_dir_results ? y : 0; */
/* z_min = spatial_wloop ? 1 : z_min; //if spatial_wloop is true, then z_min should be 1, otherwise it should be 0 */
/* for (int z=z_min; z<z_max;z+=1) */
/* { */
/* for (int t=1; t< t_max; t+=t_plus) */
/* { */
/* if (spatial_wloop and x==0 and y==0){continue;} */
/* else if (!spatial_wloop and x==0 and y==0 and z==0) {continue;} */
/* if (spatial_wloop and avg_dir_results) */
/* { */
/* //CalcWilsonBresenham */
/* //the spatial wilson loop is calculated in 3 spaces only, so if x==y==0 then its just a polyakovLoop in fictitious time direction. */
/* PREC Wloop=gaugeaction.WilsonLoop(x,y,z,t,t_dir,avg_dir_results,spatial_wloop); */
/* printf("Node:%d,Smear_count:%d,Length=%f,Height=%d,ex_dir=%d,Wilson Loop:%1.15e\n",node,no_smearing_steps,pow((x*x+y*y),0.5),z,t_dir,Wloop); */
/* } */
/* else if (spatial_wloop and !avg_dir_results) */
/* { */
/* PREC Wloop=gaugeaction.WilsonLoop(x,y,z,t,t_dir,avg_dir_results,spatial_wloop); */
/* /1* printf("Node:%d,Smear_count:%d,x=%d,y=%d,z=%d,t=%d,ex_dir=%d,Wilson Loop:%1.15e\n",node,no_smearing_steps,x,y,0,z,t_dir,Wloop); *1/ */
/* printf("WILS_LOOP_%d %d %d %d\t%d\t%1.15e\n",no_smearing_steps,x,y,z,t,Wloop); */
/* } */
/* else if (!spatial_wloop and avg_dir_results) */
/* { t_dir=3; //not a spatial wilson loop, so t_dir is always 3 */
/* PREC Wloop=gaugeaction.WilsonLoop(x,y,z,t,t_dir,avg_dir_results,spatial_wloop); */
/* printf("Node:%d,Smear_count:%d,Length=%f,Height=%d,ex_dir=%d,Wilson Loop:%1.15e\n",node,no_smearing_steps,pow((x*x+y*y+z*z),0.5),t,t_dir,Wloop); */
/* } */
/* else if(!spatial_wloop and !avg_dir_results) */
/* { t_dir=3; // not a spatial wilson loop, so t_dir is always 3 */
/* PREC Wloop=gaugeaction.WilsonLoop(x,y,z,t,t_dir,avg_dir_results,spatial_wloop); */
/* /1* printf("Node:%d,Smear_count:%d,x=%d,y=%d,z=%d,t=%d,ex_dir=%d,Wilson Loop:%1.15e\n",node,no_smearing_steps,x,y,z,t,t_dir,Wloop); *1/ */
/* printf("WILS_LOOP_%d %d %d %d\t%d\t%1.15e\n",no_smearing_steps,x,y,z,t,Wloop); */
/* } */
/* } */
/* } */
/* } */
/* } */
/* gauge=gauge_original; */
/* } */
}
return 0;
}