Compare commits

...
12 Commits
18 changed files with 202 additions and 46 deletions

No files matched your search

+2
View File
@@ -269,7 +269,9 @@ public:
RealD xscale = 2.0/(hi-lo); RealD xscale = 2.0/(hi-lo);
RealD mscale = -(hi+lo)/(hi-lo); RealD mscale = -(hi+lo)/(hi-lo);
Linop.HermOp(T0,y); Linop.HermOp(T0,y);
grid->Barrier();
axpby(T1,xscale,mscale,y,in); axpby(T1,xscale,mscale,y,in);
grid->Barrier();
// sum = .5 c[0] T0 + c[1] T1 // sum = .5 c[0] T0 + c[1] T1
// out = ()*T0 + Coeffs[1]*T1; // out = ()*T0 + Coeffs[1]*T1;
+12
View File
@@ -260,32 +260,39 @@ CartesianCommunicator::~CartesianCommunicator()
} }
#ifdef USE_GRID_REDUCTION #ifdef USE_GRID_REDUCTION
void CartesianCommunicator::GlobalSum(float &f){ void CartesianCommunicator::GlobalSum(float &f){
FlightRecorder::StepLog("GlobalSumP2P");
CartesianCommunicator::GlobalSumP2P(f); CartesianCommunicator::GlobalSumP2P(f);
} }
void CartesianCommunicator::GlobalSum(double &d) void CartesianCommunicator::GlobalSum(double &d)
{ {
FlightRecorder::StepLog("GlobalSumP2P");
CartesianCommunicator::GlobalSumP2P(d); CartesianCommunicator::GlobalSumP2P(d);
} }
#else #else
void CartesianCommunicator::GlobalSum(float &f){ void CartesianCommunicator::GlobalSum(float &f){
FlightRecorder::StepLog("AllReduce");
int ierr=MPI_Allreduce(MPI_IN_PLACE,&f,1,MPI_FLOAT,MPI_SUM,communicator); int ierr=MPI_Allreduce(MPI_IN_PLACE,&f,1,MPI_FLOAT,MPI_SUM,communicator);
assert(ierr==0); assert(ierr==0);
} }
void CartesianCommunicator::GlobalSum(double &d) void CartesianCommunicator::GlobalSum(double &d)
{ {
FlightRecorder::StepLog("AllReduce");
int ierr = MPI_Allreduce(MPI_IN_PLACE,&d,1,MPI_DOUBLE,MPI_SUM,communicator); int ierr = MPI_Allreduce(MPI_IN_PLACE,&d,1,MPI_DOUBLE,MPI_SUM,communicator);
assert(ierr==0); assert(ierr==0);
} }
#endif #endif
void CartesianCommunicator::GlobalSum(uint32_t &u){ void CartesianCommunicator::GlobalSum(uint32_t &u){
FlightRecorder::StepLog("AllReduce");
int ierr=MPI_Allreduce(MPI_IN_PLACE,&u,1,MPI_UINT32_T,MPI_SUM,communicator); int ierr=MPI_Allreduce(MPI_IN_PLACE,&u,1,MPI_UINT32_T,MPI_SUM,communicator);
assert(ierr==0); assert(ierr==0);
} }
void CartesianCommunicator::GlobalSum(uint64_t &u){ void CartesianCommunicator::GlobalSum(uint64_t &u){
FlightRecorder::StepLog("AllReduce");
int ierr=MPI_Allreduce(MPI_IN_PLACE,&u,1,MPI_UINT64_T,MPI_SUM,communicator); int ierr=MPI_Allreduce(MPI_IN_PLACE,&u,1,MPI_UINT64_T,MPI_SUM,communicator);
assert(ierr==0); assert(ierr==0);
} }
void CartesianCommunicator::GlobalSumVector(uint64_t* u,int N){ void CartesianCommunicator::GlobalSumVector(uint64_t* u,int N){
FlightRecorder::StepLog("AllReduceVector");
int ierr=MPI_Allreduce(MPI_IN_PLACE,u,N,MPI_UINT64_T,MPI_SUM,communicator); int ierr=MPI_Allreduce(MPI_IN_PLACE,u,N,MPI_UINT64_T,MPI_SUM,communicator);
assert(ierr==0); assert(ierr==0);
} }
@@ -794,6 +801,7 @@ void CartesianCommunicator::StencilSendToRecvFromComplete(std::vector<CommsReque
void CartesianCommunicator::StencilBarrier(void) void CartesianCommunicator::StencilBarrier(void)
{ {
FlightRecorder::StepLog("NodeBarrier");
MPI_Barrier (ShmComm); MPI_Barrier (ShmComm);
} }
//void CartesianCommunicator::SendToRecvFromComplete(std::vector<CommsRequest_t> &list) //void CartesianCommunicator::SendToRecvFromComplete(std::vector<CommsRequest_t> &list)
@@ -801,11 +809,13 @@ void CartesianCommunicator::StencilBarrier(void)
//} //}
void CartesianCommunicator::Barrier(void) void CartesianCommunicator::Barrier(void)
{ {
FlightRecorder::StepLog("GridBarrier");
int ierr = MPI_Barrier(communicator); int ierr = MPI_Barrier(communicator);
assert(ierr==0); assert(ierr==0);
} }
void CartesianCommunicator::Broadcast(int root,void* data, int bytes) void CartesianCommunicator::Broadcast(int root,void* data, int bytes)
{ {
FlightRecorder::StepLog("Broadcast");
int ierr=MPI_Bcast(data, int ierr=MPI_Bcast(data,
bytes, bytes,
MPI_BYTE, MPI_BYTE,
@@ -824,6 +834,7 @@ void CartesianCommunicator::BarrierWorld(void){
} }
void CartesianCommunicator::BroadcastWorld(int root,void* data, int bytes) void CartesianCommunicator::BroadcastWorld(int root,void* data, int bytes)
{ {
FlightRecorder::StepLog("BroadcastWorld");
int ierr= MPI_Bcast(data, int ierr= MPI_Bcast(data,
bytes, bytes,
MPI_BYTE, MPI_BYTE,
@@ -846,6 +857,7 @@ void CartesianCommunicator::AllToAll(int dim,void *in,void *out,uint64_t words,
} }
void CartesianCommunicator::AllToAll(void *in,void *out,uint64_t words,uint64_t bytes) void CartesianCommunicator::AllToAll(void *in,void *out,uint64_t words,uint64_t bytes)
{ {
FlightRecorder::StepLog("AllToAll");
// MPI is a pain and uses "int" arguments // MPI is a pain and uses "int" arguments
// 64*64*64*128*16 == 500Million elements of data. // 64*64*64*128*16 == 500Million elements of data.
// When 24*4 bytes multiples get 50x 10^9 >>> 2x10^9 Y2K bug. // When 24*4 bytes multiples get 50x 10^9 >>> 2x10^9 Y2K bug.
+1 -1
View File
@@ -990,7 +990,7 @@ void SharedMemory::SetCommunicator(Grid_MPI_Comm comm)
} }
#endif #endif
SharedMemoryTest(); // SharedMemoryTest();
} }
////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////
// On node barrier // On node barrier
+4 -4
View File
@@ -122,10 +122,10 @@ void GlobalSharedMemory::SharedMemoryZero(void *dest,size_t bytes)
{ {
acceleratorMemSet(dest,0,bytes); acceleratorMemSet(dest,0,bytes);
} }
void GlobalSharedMemory::SharedMemoryCopy(void *dest,void *src,size_t bytes) //void GlobalSharedMemory::SharedMemoryCopy(void *dest,void *src,size_t bytes)
{ //{
acceleratorCopyToDevice(src,dest,bytes); // acceleratorCopyToDevice(src,dest,bytes);
} //}
//////////////////////////////////////////////////////// ////////////////////////////////////////////////////////
// Global shared functionality finished // Global shared functionality finished
// Now move to per communicator functionality // Now move to per communicator functionality
+1 -1
View File
@@ -236,7 +236,7 @@ public:
template<class sobj> inline Lattice<vobj> & operator = (const sobj & r){ template<class sobj> inline Lattice<vobj> & operator = (const sobj & r){
vobj vtmp; vobj vtmp;
vtmp = r; vtmp = r;
#if 0 #if 1
deviceVector<vobj> vvtmp(1); deviceVector<vobj> vvtmp(1);
acceleratorPut(vvtmp[0],vtmp); acceleratorPut(vvtmp[0],vtmp);
vobj *vvtmp_p = & vvtmp[0]; vobj *vvtmp_p = & vvtmp[0];
+10
View File
@@ -252,6 +252,11 @@ void WilsonFlow<Gimpl>::smear(GaugeField& out, const GaugeField& in) const{
out = in; out = in;
RealD taus = 0.; 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 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(); auto start = std::chrono::high_resolution_clock::now();
evolve_step(out, taus); evolve_step(out, taus);
@@ -336,6 +341,11 @@ void WilsonFlowAdaptive<Gimpl>::smear(GaugeField& out, const GaugeField& in) con
RealD taus = 0.; RealD taus = 0.;
RealD eps = init_epsilon; RealD eps = init_epsilon;
unsigned int step = 0; unsigned int step = 0;
// Perform initial t=0 measurements
for(auto const &meas : this->functions)
meas.second(step,taus,out);
do{ do{
int step_success = evolve_step_adaptive(out, taus, eps); int step_success = evolve_step_adaptive(out, taus, eps);
step += step_success; //step will not be incremented if the integration step fails step += step_success; //step will not be incremented if the integration step fails
+1
View File
@@ -396,6 +396,7 @@ public:
Packets[i].from_rank,Packets[i].do_recv, Packets[i].from_rank,Packets[i].do_recv,
Packets[i].xbytes,Packets[i].rbytes,i); Packets[i].xbytes,Packets[i].rbytes,i);
} }
FlightRecorder::StepLog("Communicate begin has finished");
// Get comms started then run checksums // Get comms started then run checksums
// Having this PRIOR to the dslash seems to make Sunspot work... (!) // Having this PRIOR to the dslash seems to make Sunspot work... (!)
for(int i=0;i<Packets.size();i++){ for(int i=0;i<Packets.size();i++){
+2 -2
View File
@@ -251,7 +251,7 @@ inline void acceleratorCopyToDevice(const void *from,void *to,size_t bytes) { c
inline void acceleratorCopyFromDevice(const void *from,void *to,size_t bytes){ cudaMemcpy(to,from,bytes, cudaMemcpyDeviceToHost);} inline void acceleratorCopyFromDevice(const void *from,void *to,size_t bytes){ cudaMemcpy(to,from,bytes, cudaMemcpyDeviceToHost);}
inline void acceleratorMemSet(void *base,int value,size_t bytes) { cudaMemset(base,value,bytes);} inline void acceleratorMemSet(void *base,int value,size_t bytes) { cudaMemset(base,value,bytes);}
inline acceleratorEvent_t acceleratorCopyToDeviceAsynch(void *from, void *to, size_t bytes, cudaStream_t stream = copyStream) { inline acceleratorEvent_t acceleratorCopyToDeviceAsynch(void *from, void *to, size_t bytes, cudaStream_t stream = copyStream) {
acceleratorCopyToDevice(to,from,bytes, cudaMemcpyHostToDevice); acceleratorCopyToDevice(from,to,bytes);
return 0; return 0;
} }
inline acceleratorEvent_t acceleratorCopyFromDeviceAsynch(void *from, void *to, size_t bytes, cudaStream_t stream = copyStream) { inline acceleratorEvent_t acceleratorCopyFromDeviceAsynch(void *from, void *to, size_t bytes, cudaStream_t stream = copyStream) {
@@ -337,7 +337,7 @@ accelerator_inline int acceleratorSIMTlane(int Nsimd) {
cgh.parallel_for( \ cgh.parallel_for( \
sycl::nd_range<3>(global,local), \ sycl::nd_range<3>(global,local), \
[=] (sycl::nd_item<3> item) /*mutable*/ \ [=] (sycl::nd_item<3> item) /*mutable*/ \
[[intel::reqd_sub_group_size(16)]] \ [[sycl::reqd_sub_group_size(16)]] \
{ \ { \
auto iter1 = item.get_global_id(0); \ auto iter1 = item.get_global_id(0); \
auto iter2 = item.get_global_id(1); \ auto iter2 = item.get_global_id(1); \
+3 -4
View File
@@ -638,12 +638,11 @@ void Grid_debug_handler_init(void)
sa.sa_flags = SA_SIGINFO; sa.sa_flags = SA_SIGINFO;
// sigaction(SIGSEGV,&sa,NULL); // sigaction(SIGSEGV,&sa,NULL);
sigaction(SIGTRAP,&sa,NULL); sigaction(SIGTRAP,&sa,NULL);
sigaction(SIGBUS,&sa,NULL); // sigaction(SIGBUS,&sa,NULL);
// sigaction(SIGUSR2,&sa,NULL); // sigaction(SIGUSR2,&sa,NULL);
feenableexcept( FE_INVALID|FE_OVERFLOW|FE_DIVBYZERO); // feenableexcept( FE_INVALID|FE_OVERFLOW|FE_DIVBYZERO);
// sigaction(SIGFPE,&sa,NULL);
sigaction(SIGFPE,&sa,NULL);
sigaction(SIGKILL,&sa,NULL); sigaction(SIGKILL,&sa,NULL);
sigaction(SIGILL,&sa,NULL); sigaction(SIGILL,&sa,NULL);
+39 -12
View File
@@ -66,6 +66,7 @@ namespace Grid{
}; };
} }
template <class T> void writeFile(T& in, std::string const fname){ template <class T> void writeFile(T& in, std::string const fname){
#ifdef HAVE_LIME #ifdef HAVE_LIME
// Ref: https://github.com/paboyle/Grid/blob/feature/scidac-wp1/tests/debug/Test_general_coarse_hdcg_phys48.cc#L111 // Ref: https://github.com/paboyle/Grid/blob/feature/scidac-wp1/tests/debug/Test_general_coarse_hdcg_phys48.cc#L111
@@ -73,7 +74,7 @@ template <class T> void writeFile(T& in, std::string const fname){
Grid::emptyUserRecord record; Grid::emptyUserRecord record;
Grid::ScidacWriter WR(in.Grid()->IsBoss()); Grid::ScidacWriter WR(in.Grid()->IsBoss());
WR.open(fname); WR.open(fname);
WR.writeScidacFieldRecord(in,record,0); WR.writeScidacFieldRecord(in,record,0); // Lexico
WR.close(); WR.close();
#endif #endif
// What is the appropriate way to throw error? // What is the appropriate way to throw error?
@@ -107,8 +108,18 @@ int main(int argc, char **argv) {
for (int conf = CPar.StartConfiguration; conf <= CPar.EndConfiguration; conf+= CPar.Skip){ for (int conf = CPar.StartConfiguration; conf <= CPar.EndConfiguration; conf+= CPar.Skip){
#if 0
CPNersc.CheckpointRestore(conf, Umu, sRNG, pRNG); 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 << std::setprecision(15);
std::cout << GridLogMessage << "Initial plaquette: "<< WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu) << std::endl; 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); std::string file_post = CPar.conf_prefix + "." + std::to_string(conf);
WilsonFlow<PeriodicGimplR> WF(WFPar.step_size,WFPar.steps,WFPar.meas_interval); 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){ 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; 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; //double coeff = 2.0 / (1.0 * Nd * (Nd - 1)) / 3.0;
//Plq = coeff * Plq; //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 E = real(sum(R))/ RealD(U.Grid()->gSites());
RealD T = real( sum(qfield) ); RealD T = real( sum(qfield) );
Coordinate scoor; for (int mu=0; mu < Nd; mu++) scoor[mu] = 0; Coordinate scoor; for (int mu=0; mu < Nd; mu++) scoor[mu] = 0;
RealD E0 = real(peekSite(R,scoor)); RealD E0 = real(peekSite(R,scoor));
RealD T0 = real(peekSite(qfield,scoor)); RealD T0 = real(peekSite(qfield,scoor));
std::cout << GridLogMessage << "[WilsonFlow] Saved energy density (clover) & topo. charge density: " << conf << " " << step << " " << tau << " " 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; int t=WFPar.maxTau;
WF.smear(Uflow, Umu); WF.smear(Uflow, Umu);
// NerscIO::writeConfiguration(Uflow,filesmr);
RealD WFlow_plaq = WilsonLoops<PeriodicGimplR>::avgPlaquette(Uflow); RealD WFlow_plaq = WilsonLoops<PeriodicGimplR>::avgPlaquette(Uflow);
RealD WFlow_TC = WilsonLoops<PeriodicGimplR>::TopologicalCharge(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_T0 = WF.energyDensityPlaquette(t,Uflow); // t
RealD WFlow_EC = WF.energyDensityCloverleaf(t,Uflow); RealD WFlow_EC = WF.energyDensityCloverleaf(t,Uflow);
std::cout << GridLogMessage << "Plaquette "<< conf << " " << WFlow_plaq << std::endl; std::cout << GridLogMessage << "Plaquette "<< conf << " " << WFlow_plaq << std::endl;
std::cout << GridLogMessage << "T0 "<< conf << " " << WFlow_T0 << std::endl; std::cout << GridLogMessage << "T0 "<< conf << " " << WFlow_T0 << std::endl;
std::cout << GridLogMessage << "TC0 "<< conf << " " << WFlow_EC << std::endl; std::cout << GridLogMessage << "TC0 "<< conf << " " << WFlow_EC << std::endl;
std::cout << GridLogMessage << "TopologicalCharge "<< conf << " " << WFlow_TC << 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"; std::cout<< GridLogMessage << " Admissibility check:\n";
const double sp_adm = 0.067; // admissible threshold const double sp_adm = 0.067; // admissible threshold
-2
View File
@@ -7,8 +7,6 @@ CXX=mpicxx ../../configure \
--enable-unified=yes \ --enable-unified=yes \
--prefix /Users/peterboyle/QCD/vtk/Grid/install \ --prefix /Users/peterboyle/QCD/vtk/Grid/install \
--with-lime=$CLIME \ --with-lime=$CLIME \
--with-hdf5=$HDF5 \
--with-fftw=$FFTW \
--with-openssl=$OPENSSL \ --with-openssl=$OPENSSL \
--with-gmp=$GMP \ --with-gmp=$GMP \
--with-mpfr=$MPFR \ --with-mpfr=$MPFR \
+15 -2
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 \ ../../configure \
--enable-comms=mpi-auto \ --enable-comms=mpi-auto \
--enable-unified=yes \ --enable-unified=yes \
@@ -5,12 +14,16 @@
--enable-shm-fast-path=shmopen \ --enable-shm-fast-path=shmopen \
--enable-accelerator=none \ --enable-accelerator=none \
--enable-simd=AVX512 \ --enable-simd=AVX512 \
--disable-accelerator-cshift \ --with-lime=$CLIME \
--with-hdf5=$HDF5 \
--with-fftw=$FFTW \
--disable-fermion-reps \ --disable-fermion-reps \
--disable-gparity \ --disable-gparity \
CXX=clang++ \ CXX=clang++ \
MPICXX=mpicxx \ MPICXX=mpicxx \
CXXFLAGS="-std=c++17" LIBS=-llime \
LDFLAGS=-L$CLIME/lib/ \
CXXFLAGS="-std=c++17 -fPIE"
+2 -1
View File
@@ -1,4 +1,5 @@
source $HOME/spack/share/spack/setup-env.sh source $HOME/spack/share/spack/setup-env.sh
spack load llvm@17.0.4 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 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
+74 -7
View File
@@ -31,12 +31,23 @@ directory
using namespace std; using namespace std;
using namespace Grid; using namespace Grid;
;
//typedef WilsonFermionD FermionOp; //typedef WilsonFermionD FermionOp;
typedef DomainWallFermionD FermionOp; typedef DomainWallFermionD FermionOp;
typedef typename DomainWallFermionD::FermionField FermionField; 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.; } RealD AllZero(RealD x) { return 0.; }
@@ -121,7 +132,7 @@ int main(int argc, char** argv) {
int Ls=16; int Ls=16;
RealD M5=1.8; RealD M5=1.8;
RealD mass = -1.0; RealD mass = 0.01;
mass=LanParams.mass; mass=LanParams.mass;
Ls=LanParams.Ls; Ls=LanParams.Ls;
@@ -159,10 +170,10 @@ int main(int argc, char** argv) {
U[mu] = PeekIndex<LorentzIndex>(Umu, mu); U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
} }
*/ */
int Nstop = 10;
int Nk = 20; int Nk = 20;
int Nstop = Nk;
int Np = 80; int Np = 80;
Nstop=LanParams.Nstop; Nstop=LanParams.Nstop;
Nk=LanParams.Nk; Nk=LanParams.Nk;
Np=LanParams.Np; Np=LanParams.Np;
@@ -201,8 +212,12 @@ int main(int argc, char** argv) {
int Nconv; int Nconv;
IRL.calc(eval, evec, src, 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 #if 0
Gamma g5(Gamma::Algebra::Gamma5) ; Gamma g5(Gamma::Algebra::Gamma5) ;
ComplexD dot; ComplexD dot;
@@ -232,6 +247,7 @@ int main(int argc, char** argv) {
vector<LatticeFermion> finalevec(Nconv, FGrid); vector<LatticeFermion> finalevec(Nconv, FGrid);
vector<RealD> eMe(Nconv), eMMe(Nconv); vector<RealD> eMe(Nconv), eMMe(Nconv);
for(int i = 0; i < Nconv; i++){ for(int i = 0; i < Nconv; i++){
cout << "calculate the matrix element["<<i<<"]" << endl;
G5R5Herm.HermOpAndNorm(evec[i], G5R5Mevec[i], eMe[i], eMMe[i]); G5R5Herm.HermOpAndNorm(evec[i], G5R5Mevec[i], eMe[i], eMMe[i]);
} }
cout << "Re<evec, G5R5M(evec)>: " << endl; 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); axpby_ssp(G5evec[i], -1., finalevec[i], 0., G5evec[i], j, j);
} }
} }
for(int i = 0; i < Nconv; i++){ for(int i = 0; i < Nconv; i++){
chiral_matrix_real[i].resize(Nconv); chiral_matrix_real[i].resize(Nconv);
chiral_matrix[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++){ for(int j = 0; j < Nconv; j++){
chiral_matrix[i][j] = innerProduct(finalevec[i], G5evec[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]); chiral_matrix_real[i][j] = abs(chiral_matrix[i][j]);
std::cout <<" chiral_matrix_real "<<i<<" "<<j<<" "<< chiral_matrix_real[i][j] << std::endl; 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++){ 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(); Grid_finalize();
} }
+1 -1
View File
@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.12 FATAL_ERROR)
project(GridViewer) project(GridViewer)
list(APPEND CMAKE_PREFIX_PATH "/Users/peterboyle/QCD/vtk/VTK-9.4.2-install/") list(APPEND CMAKE_PREFIX_PATH "/home/paboyle/Visualisation/install/")
find_package(VTK COMPONENTS find_package(VTK COMPONENTS
CommonColor CommonColor
+7 -4
View File
@@ -48,15 +48,14 @@ typedef vtkMarchingCubes isosurface;
int mpeg = 0 ; int mpeg = 0 ;
int xlate = 0 ; int xlate = 0 ;
int framerate = 10;
template <class T> void readFile(T& out, std::string const fname){ template <class T> void readFile(T& out, std::string const fname){
#ifdef HAVE_LIME
Grid::emptyUserRecord record; Grid::emptyUserRecord record;
Grid::ScidacReader RD; Grid::ScidacReader RD;
RD.open(fname); RD.open(fname);
RD.readScidacFieldRecord(out,record); RD.readScidacFieldRecord(out,record);
RD.close(); RD.close();
#endif
} }
using namespace Grid; using namespace Grid;
@@ -208,6 +207,10 @@ int main(int argc, char* argv[])
xlate = 1; xlate = 1;
} }
if( GridCmdOptionExists(argv,argv+argc,"--fps") ){
arg=GridCmdOptionPayload(argv,argv+argc,"--fps");
GridCmdOptionInt(arg,framerate);
}
if( GridCmdOptionExists(argv,argv+argc,"--isosurface") ){ if( GridCmdOptionExists(argv,argv+argc,"--isosurface") ){
arg=GridCmdOptionPayload(argv,argv+argc,"--isosurface"); arg=GridCmdOptionPayload(argv,argv+argc,"--isosurface");
GridCmdOptionFloat(arg,default_contour); GridCmdOptionFloat(arg,default_contour);
@@ -420,7 +423,7 @@ int main(int argc, char* argv[])
vtkFFMPEGWriter *writer = vtkFFMPEGWriter::New(); vtkFFMPEGWriter *writer = vtkFFMPEGWriter::New();
writer->SetFileName("movie.avi"); writer->SetFileName("movie.avi");
writer->SetRate(1); writer->SetRate(framerate);
writer->SetInputConnection(imageFilter->GetOutputPort()); writer->SetInputConnection(imageFilter->GetOutputPort());
writer->Start(); writer->Start();
@@ -477,7 +480,7 @@ int main(int argc, char* argv[])
slidercallback->fu_list = fu_list; slidercallback->fu_list = fu_list;
sliderWidget->AddObserver(vtkCommand::InteractionEvent, slidercallback); sliderWidget->AddObserver(vtkCommand::InteractionEvent, slidercallback);
int timerId = iren->CreateRepeatingTimer(300); int timerId = iren->CreateRepeatingTimer(1000/framerate);
std::cout << "timerId: " << timerId << std::endl; std::cout << "timerId: " << timerId << std::endl;
// Start the interaction and timer // Start the interaction and timer
+15
View File
@@ -73,6 +73,21 @@ each to:
VTK really should make it easier to pick up the flags required for FFMPEG linkage, especially as they are very quirky on MacOS. VTK really should make it easier to pick up the flags required for FFMPEG linkage, especially as they are very quirky on MacOS.
========================================
Aurora compilation:
========================================
module load ffmpeg
download & untar: VTK-7.0.2
mkdir build
cd build
ccmake ../
"t"
Enable: VTK_MODULE_ENABLE_VTK_IOFFMPEG YES
"configure" ; should "discover" the installed ffmpeg module
Still need an "X" connection to make the MPEG files.
======================================== ========================================
Grid: Grid:
+13 -5
View File
@@ -1,9 +1,17 @@
libs=`grid-config --libs` export grid_config=/home/paboyle/GPT/install/bin/grid-config
ldflags=`grid-config --ldflags` libs=`$grid_config --libs`
cxxflags=`grid-config --cxxflags` ldflags=`$grid_config --ldflags`
cxx=`grid-config --cxx` cxxflags=`$grid_config --cxxflags`
cxx=`$grid_config --cxx`
cc=icx
mkdir build mkdir build
cd build cd build
LDFLAGS="$ldflags $libs " cmake .. -DCMAKE_CXX_COMPILER=$cxx -DCMAKE_CXX_FLAGS=$cxxflags echo CC $cc
echo CXX $cxx
echo CXXFLAGS $cxxflags
echo LDFLAGS $ldflags
echo LIBS $libs
LDFLAGS="$ldflags $libs " cmake .. -DCMAKE_C_COMPILER=$cc -DCMAKE_CXX_COMPILER="$cxx" -DCMAKE_CXX_FLAGS="$cxxflags "