1
0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-07 09:57:07 +01:00

Compare commits

...

4 Commits

Author SHA1 Message Date
Peter Boyle
dc546aaa4b Updated config options for BNL cluster 2025-05-13 18:44:47 -04:00
Peter Boyle
5364d580c9 Output chirality, eigenvector density files and python source lego plot 2025-05-13 18:44:47 -04:00
Peter Boyle
2a9a6347e3 Do not require Grid format RNGs and also to the 5Li reporting 2025-05-13 18:44:47 -04:00
Peter Boyle
cfdb56f314 Run measurements at t=0 too 2025-05-13 18:44:46 -04:00
5 changed files with 139 additions and 21 deletions

View File

@ -252,6 +252,11 @@ void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const{
out = in;
RealD taus = 0.;
// Perform initial t=0 measurements
for(auto const &meas : this->functions)
meas.second(0,taus,out);
for (unsigned int step = 1; step <= Nstep; step++) { //step indicates the number of smearing steps applied at the time of measurement
auto start = std::chrono::high_resolution_clock::now();
evolve_step(out, taus);
@ -336,6 +341,11 @@ void WilsonFlowAdaptive<Gimpl>::smear(GaugeField& out, const GaugeField& in) con
RealD taus = 0.;
RealD eps = init_epsilon;
unsigned int step = 0;
// Perform initial t=0 measurements
for(auto const &meas : this->functions)
meas.second(step,taus,out);
do{
int step_success = evolve_step_adaptive(out, taus, eps);
step += step_success; //step will not be incremented if the integration step fails

View File

@ -66,6 +66,7 @@ namespace Grid{
};
}
template <class T> void writeFile(T& in, std::string const fname){
#ifdef HAVE_LIME
// Ref: https://github.com/paboyle/Grid/blob/feature/scidac-wp1/tests/debug/Test_general_coarse_hdcg_phys48.cc#L111
@ -107,8 +108,18 @@ int main(int argc, char **argv) {
for (int conf = CPar.StartConfiguration; conf <= CPar.EndConfiguration; conf+= CPar.Skip){
#if 0
CPNersc.CheckpointRestore(conf, Umu, sRNG, pRNG);
#else
// Don't require Grid format RNGs
FieldMetaData header;
std::string file, filesmr;
file = CPar.conf_path + "/" + CPar.conf_prefix + "." + std::to_string(conf);
filesmr = CPar.conf_path + "/" + CPar.conf_smr_prefix + "." + std::to_string(conf);
NerscIO::readConfiguration(Umu,header,file);
#endif
std::cout << std::setprecision(15);
std::cout << GridLogMessage << "Initial plaquette: "<< WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu) << std::endl;
@ -116,6 +127,7 @@ int main(int argc, char **argv) {
std::string file_post = CPar.conf_prefix + "." + std::to_string(conf);
WilsonFlow<PeriodicGimplR> WF(WFPar.step_size,WFPar.steps,WFPar.meas_interval);
WF.addMeasurement(WFPar.meas_interval_density, [&file_pre,&file_post,&conf](int step, RealD t, const typename PeriodicGimplR::GaugeField &U){
typedef typename PeriodicGimplR::GaugeLinkField GaugeMat;
@ -165,33 +177,48 @@ int main(int argc, char **argv) {
//double coeff = 2.0 / (1.0 * Nd * (Nd - 1)) / 3.0;
//Plq = coeff * Plq;
int tau = std::round(t);
std::string efile = file_pre + "E_dnsty_" + std::to_string(tau) + "_" + file_post;
writeFile(R,efile);
std::string tfile = file_pre + "Top_dnsty_" + std::to_string(tau) + "_" + file_post;
writeFile(qfield,tfile);
RealD WFlow_TC5Li = WilsonLoops<PeriodicGimplR>::TopologicalCharge5Li(U);
int tau = std::round(t);
std::string efile = file_pre + "E_dnsty_" + std::to_string(tau) + "_" + file_post;
// writeFile(R,efile);
std::string tfile = file_pre + "Top_dnsty_" + std::to_string(tau) + "_" + file_post;
// writeFile(qfield,tfile);
std::string ufile = file_pre + "U_" + std::to_string(tau) + "_" + file_post;
{
PeriodicGimplR::GaugeField Ucopy = U;
NerscIO::writeConfiguration(Ucopy,ufile);
}
RealD E = real(sum(R))/ RealD(U.Grid()->gSites());
RealD T = real( sum(qfield) );
Coordinate scoor; for (int mu=0; mu < Nd; mu++) scoor[mu] = 0;
RealD E0 = real(peekSite(R,scoor));
RealD T0 = real(peekSite(qfield,scoor));
std::cout << GridLogMessage << "[WilsonFlow] Saved energy density (clover) & topo. charge density: " << conf << " " << step << " " << tau << " "
<< "(E_avg,T_sum) " << E << " " << T << " (E, T at origin) " << E0 << " " << T0 << std::endl;
<< "(E_avg,T_sum) " << E << " " << T << " (E, T at origin) " << E0 << " " << T0 << " Q5Li "<< WFlow_TC5Li << std::endl;
});
int t=WFPar.maxTau;
WF.smear(Uflow, Umu);
NerscIO::writeConfiguration(Uflow,filesmr);
RealD WFlow_plaq = WilsonLoops<PeriodicGimplR>::avgPlaquette(Uflow);
RealD WFlow_TC = WilsonLoops<PeriodicGimplR>::TopologicalCharge(Uflow);
RealD WFlow_TC5Li = WilsonLoops<PeriodicGimplR>::TopologicalCharge5Li(Uflow);
RealD WFlow_T0 = WF.energyDensityPlaquette(t,Uflow); // t
RealD WFlow_EC = WF.energyDensityCloverleaf(t,Uflow);
std::cout << GridLogMessage << "Plaquette "<< conf << " " << WFlow_plaq << std::endl;
std::cout << GridLogMessage << "T0 "<< conf << " " << WFlow_T0 << std::endl;
std::cout << GridLogMessage << "TC0 "<< conf << " " << WFlow_EC << std::endl;
std::cout << GridLogMessage << "TopologicalCharge "<< conf << " " << WFlow_TC << std::endl;
std::cout << GridLogMessage << "Plaquette "<< conf << " " << WFlow_plaq << std::endl;
std::cout << GridLogMessage << "T0 "<< conf << " " << WFlow_T0 << std::endl;
std::cout << GridLogMessage << "TC0 "<< conf << " " << WFlow_EC << std::endl;
std::cout << GridLogMessage << "TopologicalCharge "<< conf << " " << WFlow_TC << std::endl;
std::cout << GridLogMessage << "TopologicalCharge5Li "<< conf << " " << WFlow_TC5Li<< std::endl;
std::cout<< GridLogMessage << " Admissibility check:\n";
const double sp_adm = 0.067; // admissible threshold

View File

@ -1,3 +1,12 @@
spack load c-lime
spack load fftw
spack load hdf5+cxx
export FFTW=`spack find --paths fftw | grep ^fftw | awk '{print $2}' `
export HDF5=`spack find --paths hdf5+cxx | grep ^hdf5 | awk '{print $2}' `
export CLIME=`spack find --paths c-lime | grep ^c-lime | awk '{print $2}' `
../../configure \
--enable-comms=mpi-auto \
--enable-unified=yes \
@ -5,12 +14,16 @@
--enable-shm-fast-path=shmopen \
--enable-accelerator=none \
--enable-simd=AVX512 \
--disable-accelerator-cshift \
--with-lime=$CLIME \
--with-hdf5=$HDF5 \
--with-fftw=$FFTW \
--disable-fermion-reps \
--disable-gparity \
CXX=clang++ \
MPICXX=mpicxx \
CXXFLAGS="-std=c++17"
LIBS=-llime \
LDFLAGS=-L$CLIME/lib/ \
CXXFLAGS="-std=c++17 -fPIE"

View File

@ -1,4 +1,5 @@
source $HOME/spack/share/spack/setup-env.sh
spack load llvm@17.0.4
export LD_LIBRARY_PATH=/direct/sdcc+u/paboyle/spack/opt/spack/linux-almalinux8-icelake/gcc-8.5.0/llvm-17.0.4-laufdrcip63ivkadmtgoepwmj3dtztdu/lib:$LD_LIBRARY_PATH
module load openmpi
module load openmpi/4.1.8
spack load c-lime

View File

@ -31,12 +31,23 @@ directory
using namespace std;
using namespace Grid;
;
//typedef WilsonFermionD FermionOp;
typedef DomainWallFermionD FermionOp;
typedef typename DomainWallFermionD::FermionField FermionField;
template <class T> void writeFile(T& in, std::string const fname){
#ifdef HAVE_LIME
// Ref: https://github.com/paboyle/Grid/blob/feature/scidac-wp1/tests/debug/Test_general_coarse_hdcg_phys48.cc#L111
std::cout << Grid::GridLogMessage << "Writes to: " << fname << std::endl;
Grid::emptyUserRecord record;
Grid::ScidacWriter WR(in.Grid()->IsBoss());
WR.open(fname);
WR.writeScidacFieldRecord(in,record,0);
WR.close();
#endif
// What is the appropriate way to throw error?
}
RealD AllZero(RealD x) { return 0.; }
@ -121,7 +132,7 @@ int main(int argc, char** argv) {
int Ls=16;
RealD M5=1.8;
RealD mass = -1.0;
RealD mass = 0.01;
mass=LanParams.mass;
Ls=LanParams.Ls;
@ -159,10 +170,10 @@ int main(int argc, char** argv) {
U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
}
*/
int Nstop = 10;
int Nk = 20;
int Nstop = Nk;
int Np = 80;
Nstop=LanParams.Nstop;
Nk=LanParams.Nk;
Np=LanParams.Np;
@ -201,8 +212,12 @@ int main(int argc, char** argv) {
int Nconv;
IRL.calc(eval, evec, src, Nconv);
std::cout << mass <<" : " << eval << std::endl;
std::cout << mass <<" : " << eval << std::endl;
std::cout << " #evecs " << evec.size() << std::endl;
std::cout << " Nconv " << Nconv << std::endl;
std::cout << " Nm " << Nm << std::endl;
if ( Nconv > evec.size() ) Nconv = evec.size();
#if 0
Gamma g5(Gamma::Algebra::Gamma5) ;
ComplexD dot;
@ -232,6 +247,7 @@ int main(int argc, char** argv) {
vector<LatticeFermion> finalevec(Nconv, FGrid);
vector<RealD> eMe(Nconv), eMMe(Nconv);
for(int i = 0; i < Nconv; i++){
cout << "calculate the matrix element["<<i<<"]" << endl;
G5R5Herm.HermOpAndNorm(evec[i], G5R5Mevec[i], eMe[i], eMMe[i]);
}
cout << "Re<evec, G5R5M(evec)>: " << endl;
@ -326,13 +342,27 @@ int main(int argc, char** argv) {
axpby_ssp(G5evec[i], -1., finalevec[i], 0., G5evec[i], j, j);
}
}
for(int i = 0; i < Nconv; i++){
chiral_matrix_real[i].resize(Nconv);
chiral_matrix[i].resize(Nconv);
std::string evfile("./evec_density");
evfile = evfile+"_"+std::to_string(i);
auto evdensity = localInnerProduct(finalevec[i],finalevec[i] );
writeFile(evdensity,evfile);
for(int j = 0; j < Nconv; j++){
chiral_matrix[i][j] = innerProduct(finalevec[i], G5evec[j]);
std::cout <<" chiral_matrix_real signed "<<i<<" "<<j<<" "<< chiral_matrix_real[i][j] << std::endl;
chiral_matrix_real[i][j] = abs(chiral_matrix[i][j]);
std::cout <<" chiral_matrix_real "<<i<<" "<<j<<" "<< chiral_matrix_real[i][j] << std::endl;
if ( chiral_matrix_real[i][j] > 0.8 ) {
auto g5density = localInnerProduct(finalevec[i], G5evec[j]);
std::string chfile("./chiral_density_");
chfile = chfile +std::to_string(i)+"_"+std::to_string(j);
writeFile(g5density,chfile);
}
}
}
for(int i = 0; i < Nconv; i++){
@ -341,6 +371,43 @@ int main(int argc, char** argv) {
}
}
FILE *fp = fopen("lego-plot.py","w"); assert(fp!=NULL);
#define PYTHON_LINE(A) fprintf(fp,A"\n");
PYTHON_LINE("import matplotlib.pyplot as plt");
PYTHON_LINE("import numpy as np");
PYTHON_LINE("");
PYTHON_LINE("fig = plt.figure()");
PYTHON_LINE("ax = fig.add_subplot(projection='3d')");
PYTHON_LINE("");
PYTHON_LINE("x, y = np.random.rand(2, 100) * 4");
PYTHON_LINE("hist, xedges, yedges = np.histogram2d(x, y, bins=10, range=[[0, 9], [0, 9]])");
PYTHON_LINE("");
PYTHON_LINE("# Construct arrays for the anchor positions of the 16 bars");
PYTHON_LINE("xpos, ypos = np.meshgrid(xedges[:-1] + 0.25, yedges[:-1] + 0.25, indexing=\"ij\")");
PYTHON_LINE("xpos = xpos.ravel()");
PYTHON_LINE("ypos = ypos.ravel()");
PYTHON_LINE("zpos = 0");
PYTHON_LINE("");
PYTHON_LINE("# Construct arrays with the dimensions for the 16 bars.");
PYTHON_LINE("dx = dy = 0.5 * np.ones_like(zpos)");
PYTHON_LINE("dz = np.array([");
for(int i = 0; i < Nconv; i++){
fprintf(fp,"\t[ ");
for(int j = 0; j < Nconv; j++){
fprintf(fp,"%lf ",chiral_matrix_real[i][j]);
if(j<Nconv-1) fprintf(fp,",");
else fprintf(fp," ");
}
fprintf(fp,"]");
if(i<Nconv-1) fprintf(fp,",\n");
else fprintf(fp,"\n");
}
PYTHON_LINE("\t])");
PYTHON_LINE("dz = dz.ravel()");
PYTHON_LINE("ax.bar3d(xpos, ypos, zpos, dx, dy, dz, zsort='average')");
PYTHON_LINE("plt.show()");
fclose(fp);
Grid_finalize();
}