Updated python output

Grid/algorithms/approx/Chebyshev.h

@@ -269,7 +269,9 @@ public:
     RealD xscale = 2.0/(hi-lo);
     RealD mscale = -(hi+lo)/(hi-lo);
     Linop.HermOp(T0,y);
+    grid->Barrier();
     axpby(T1,xscale,mscale,y,in);
+    grid->Barrier();
 
     // sum = .5 c[0] T0 + c[1] T1
     //    out = ()*T0 + Coeffs[1]*T1;

Grid/communicator/Communicator_mpi3.cc

@@ -260,32 +260,39 @@ CartesianCommunicator::~CartesianCommunicator()
 }
 #ifdef USE_GRID_REDUCTION
 void CartesianCommunicator::GlobalSum(float &f){
+  FlightRecorder::StepLog("GlobalSumP2P");
   CartesianCommunicator::GlobalSumP2P(f);
 }
 void CartesianCommunicator::GlobalSum(double &d)
 {
+  FlightRecorder::StepLog("GlobalSumP2P");
   CartesianCommunicator::GlobalSumP2P(d);
 }
 #else
 void CartesianCommunicator::GlobalSum(float &f){
+  FlightRecorder::StepLog("AllReduce");
   int ierr=MPI_Allreduce(MPI_IN_PLACE,&f,1,MPI_FLOAT,MPI_SUM,communicator);
   assert(ierr==0);
 }
 void CartesianCommunicator::GlobalSum(double &d)
 {
+  FlightRecorder::StepLog("AllReduce");
   int ierr = MPI_Allreduce(MPI_IN_PLACE,&d,1,MPI_DOUBLE,MPI_SUM,communicator);
   assert(ierr==0);
 }
 #endif
 void CartesianCommunicator::GlobalSum(uint32_t &u){
+  FlightRecorder::StepLog("AllReduce");
   int ierr=MPI_Allreduce(MPI_IN_PLACE,&u,1,MPI_UINT32_T,MPI_SUM,communicator);
   assert(ierr==0);
 }
 void CartesianCommunicator::GlobalSum(uint64_t &u){
+  FlightRecorder::StepLog("AllReduce");
   int ierr=MPI_Allreduce(MPI_IN_PLACE,&u,1,MPI_UINT64_T,MPI_SUM,communicator);
   assert(ierr==0);
 }
 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);
   assert(ierr==0);
 }
@@ -794,6 +801,7 @@ void CartesianCommunicator::StencilSendToRecvFromComplete(std::vector<CommsReque
 
 void CartesianCommunicator::StencilBarrier(void)
 {
+  FlightRecorder::StepLog("NodeBarrier");
   MPI_Barrier  (ShmComm);
 }
 //void CartesianCommunicator::SendToRecvFromComplete(std::vector<CommsRequest_t> &list)
@@ -801,11 +809,13 @@ void CartesianCommunicator::StencilBarrier(void)
 //}
 void CartesianCommunicator::Barrier(void)
 {
+  FlightRecorder::StepLog("GridBarrier");
   int ierr = MPI_Barrier(communicator);
   assert(ierr==0);
 }
 void CartesianCommunicator::Broadcast(int root,void* data, int bytes)
 {
+  FlightRecorder::StepLog("Broadcast");
   int ierr=MPI_Bcast(data,
 		     bytes,
 		     MPI_BYTE,
@@ -824,6 +834,7 @@ void CartesianCommunicator::BarrierWorld(void){
 }
 void CartesianCommunicator::BroadcastWorld(int root,void* data, int bytes)
 {
+  FlightRecorder::StepLog("BroadcastWorld");
   int ierr= MPI_Bcast(data,
 		      bytes,
 		      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)
 {
+  FlightRecorder::StepLog("AllToAll");
   // MPI is a pain and uses "int" arguments
   // 64*64*64*128*16 == 500Million elements of data.
   // When 24*4 bytes multiples get 50x 10^9 >>> 2x10^9 Y2K bug.

Grid/communicator/SharedMemoryMPI.cc

@@ -990,7 +990,7 @@ void SharedMemory::SetCommunicator(Grid_MPI_Comm comm)
   }
 #endif
 
-  SharedMemoryTest();
+  //  SharedMemoryTest();
 }
 //////////////////////////////////////////////////////////////////
 // On node barrier

Grid/communicator/SharedMemoryNone.cc

@@ -122,10 +122,10 @@ void GlobalSharedMemory::SharedMemoryZero(void *dest,size_t bytes)
 {
   acceleratorMemSet(dest,0,bytes);
 }
-void GlobalSharedMemory::SharedMemoryCopy(void *dest,void *src,size_t bytes)
-{
-  acceleratorCopyToDevice(src,dest,bytes);
-}
+//void GlobalSharedMemory::SharedMemoryCopy(void *dest,void *src,size_t bytes)
+//{
+//  acceleratorCopyToDevice(src,dest,bytes);
+//}
 ////////////////////////////////////////////////////////
 // Global shared functionality finished
 // Now move to per communicator functionality

Grid/lattice/Lattice_base.h

@@ -236,7 +236,7 @@ public:
   template<class sobj> inline Lattice<vobj> & operator = (const sobj & r){
     vobj vtmp;
     vtmp = r;
-#if 0
+#if 1
     deviceVector<vobj> vvtmp(1);
     acceleratorPut(vvtmp[0],vtmp);
     vobj *vvtmp_p = & vvtmp[0];

Grid/qcd/smearing/WilsonFlow.h

@@ -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

Grid/stencil/Stencil.h

@@ -396,6 +396,7 @@ public:
 					Packets[i].from_rank,Packets[i].do_recv,
 					Packets[i].xbytes,Packets[i].rbytes,i);
     }
+    FlightRecorder::StepLog("Communicate begin has finished");
     // Get comms started then run checksums
     // Having this PRIOR to the dslash seems to make Sunspot work... (!)
     for(int i=0;i<Packets.size();i++){

Grid/threads/Accelerator.h

@@ -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 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) {
-  acceleratorCopyToDevice(to,from,bytes, cudaMemcpyHostToDevice);
+  acceleratorCopyToDevice(from,to,bytes);
   return 0;
 }
 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(							\
 		     sycl::nd_range<3>(global,local),			\
 		     [=] (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 iter2    = item.get_global_id(1);		\

Grid/util/Init.cc

@@ -638,12 +638,11 @@ void Grid_debug_handler_init(void)
   sa.sa_flags    = SA_SIGINFO;
   //  sigaction(SIGSEGV,&sa,NULL);
   sigaction(SIGTRAP,&sa,NULL);
-  sigaction(SIGBUS,&sa,NULL);
+  //  sigaction(SIGBUS,&sa,NULL);
   //  sigaction(SIGUSR2,&sa,NULL);
 
-  feenableexcept( FE_INVALID|FE_OVERFLOW|FE_DIVBYZERO);
-
-  sigaction(SIGFPE,&sa,NULL);
+  //  feenableexcept( FE_INVALID|FE_OVERFLOW|FE_DIVBYZERO);
+  //  sigaction(SIGFPE,&sa,NULL);
   sigaction(SIGKILL,&sa,NULL);
   sigaction(SIGILL,&sa,NULL);
 

HMC/ComputeWilsonFlow.cc

@@ -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
@@ -73,7 +74,7 @@ template <class T> void writeFile(T& in, std::string const fname){
   Grid::emptyUserRecord record;
   Grid::ScidacWriter WR(in.Grid()->IsBoss());
   WR.open(fname);
-  WR.writeScidacFieldRecord(in,record,0);
+  WR.writeScidacFieldRecord(in,record,0); // Lexico
   WR.close();
 #endif
   // 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){
 
+#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

systems/mac-arm/config-command-mpi

@@ -7,8 +7,6 @@ CXX=mpicxx ../../configure \
 	   --enable-unified=yes \
 	   --prefix /Users/peterboyle/QCD/vtk/Grid/install \
 	   --with-lime=$CLIME \
-	   --with-hdf5=$HDF5 \
-	   --with-fftw=$FFTW \
 	   --with-openssl=$OPENSSL \
 	   --with-gmp=$GMP \
 	   --with-mpfr=$MPFR \

systems/sdcc-genoa/config-command

@@ -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"
 
 
 

systems/sdcc-genoa/sourceme.sh

@@ -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

tests/lanczos/Test_dwf_G5R5.cc

@@ -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;
@@ -298,7 +314,7 @@ int main(int argc, char** argv) {
       }
     }
   }
-    for(int i = 0; i < Nconv; i++){
+  for(int i = 0; i < Nconv; i++){
     G5R5Herm.HermOpAndNorm(finalevec[i], G5R5Mevec[i], eMe[i], eMMe[i]);
   }
   cout << "Re<evec, G5R5M(evec)>: " << endl;
@@ -306,6 +322,7 @@ int main(int argc, char** argv) {
   cout << "<G5R5M(evec), G5R5M(evec)>" << endl;
   cout << eMMe << endl;
 
+  
 
 //  vector<LatticeFermion> finalevec(Nconv, FGrid);
 // temporary, until doing rotation
@@ -326,13 +343,41 @@ 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++){
+    Ddwf.M(finalevec[i], G5R5Mevec[i]);
+    for(int j = 0; j < Nconv; j++){
+      std::cout << "<"<<j<<"|Ddwf|"<<i<<"> = "<<innerProduct(finalevec[j],G5R5Mevec[i])<<std::endl;
+    }
+  }
+  for(int i = 0; i < Nconv; i++){
+    RealD t1,t2;
+    G5R5Herm.HermOpAndNorm(finalevec[i], G5R5Mevec[i], t1, t2);
+    for(int j = 0; j < Nconv; j++){
+      std::cout << "<"<<j<<"|G5R5 M|"<<i<<"> = "<<innerProduct(finalevec[j],G5R5Mevec[i])<<std::endl;
+    }
+  }
+  
   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 +386,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");
+  fprintf(fp,"hist, xedges, yedges = np.histogram2d(x, y, bins=%d, range=[[0, %d], [0, %d]])\n",Nconv,Nconv-1,Nconv-1);
+  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();
 }

visualisation/CMakeLists.txt

@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.12 FATAL_ERROR)
 
 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 
   CommonColor

visualisation/FieldDensityAnimate.cxx

@@ -48,15 +48,14 @@ typedef vtkMarchingCubes isosurface;
 
 int mpeg = 0 ;
 int xlate = 0 ;
+int framerate = 10;
 
 template <class T> void readFile(T& out, std::string const fname){
-#ifdef HAVE_LIME
   Grid::emptyUserRecord record;
   Grid::ScidacReader RD;
   RD.open(fname);
   RD.readScidacFieldRecord(out,record);
   RD.close();
-#endif
 }
 using namespace Grid;
 
@@ -208,6 +207,10 @@ int main(int argc, char* argv[])
     xlate = 1;
   }
 
+  if( GridCmdOptionExists(argv,argv+argc,"--fps") ){
+    arg=GridCmdOptionPayload(argv,argv+argc,"--fps");
+    GridCmdOptionInt(arg,framerate);
+  }
   if( GridCmdOptionExists(argv,argv+argc,"--isosurface") ){
     arg=GridCmdOptionPayload(argv,argv+argc,"--isosurface");
     GridCmdOptionFloat(arg,default_contour);
@@ -420,7 +423,7 @@ int main(int argc, char* argv[])
     
     vtkFFMPEGWriter *writer = vtkFFMPEGWriter::New();
     writer->SetFileName("movie.avi");
-    writer->SetRate(1);
+    writer->SetRate(framerate);
     writer->SetInputConnection(imageFilter->GetOutputPort());
     writer->Start();
 
@@ -477,7 +480,7 @@ int main(int argc, char* argv[])
     slidercallback->fu_list = fu_list;
     sliderWidget->AddObserver(vtkCommand::InteractionEvent, slidercallback);
 
-    int timerId = iren->CreateRepeatingTimer(300);
+    int timerId = iren->CreateRepeatingTimer(1000/framerate);
     std::cout << "timerId: " << timerId << std::endl;
 
     // Start the interaction and timer

visualisation/README

@@ -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.
 
+========================================
+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:

visualisation/cmake-command

@@ -1,9 +1,17 @@
-libs=`grid-config --libs`
-ldflags=`grid-config --ldflags`
-cxxflags=`grid-config --cxxflags`
-cxx=`grid-config --cxx`
+export grid_config=/home/paboyle/GPT/install/bin/grid-config
+libs=`$grid_config --libs`
+ldflags=`$grid_config --ldflags`
+cxxflags=`$grid_config --cxxflags`
+cxx=`$grid_config --cxx`
+cc=icx
 
 mkdir 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 "