diff --git a/Grid/Makefile.am b/Grid/Makefile.am index 8c7cb5c5..45d391ac 100644 --- a/Grid/Makefile.am +++ b/Grid/Makefile.am @@ -54,22 +54,24 @@ Version.h: version-cache include Make.inc include Eigen.inc -extra_sources+=$(WILS_FERMION_FILES) -extra_sources+=$(STAG_FERMION_FILES) +if BUILD_FERMION_INSTANTIATIONS + extra_sources+=$(WILS_FERMION_FILES) + extra_sources+=$(STAG_FERMION_FILES) if BUILD_ZMOBIUS - extra_sources+=$(ZWILS_FERMION_FILES) + extra_sources+=$(ZWILS_FERMION_FILES) endif if BUILD_GPARITY - extra_sources+=$(GP_FERMION_FILES) + extra_sources+=$(GP_FERMION_FILES) endif if BUILD_FERMION_REPS - extra_sources+=$(ADJ_FERMION_FILES) - extra_sources+=$(TWOIND_FERMION_FILES) + extra_sources+=$(ADJ_FERMION_FILES) + extra_sources+=$(TWOIND_FERMION_FILES) endif if BUILD_SP extra_sources+=$(SP_FERMION_FILES) if BUILD_FERMION_REPS - extra_sources+=$(SP_TWOIND_FERMION_FILES) + extra_sources+=$(SP_TWOIND_FERMION_FILES) +endif endif endif diff --git a/Grid/algorithms/blas/BatchedBlas.h b/Grid/algorithms/blas/BatchedBlas.h index bd01ab74..580e8166 100644 --- a/Grid/algorithms/blas/BatchedBlas.h +++ b/Grid/algorithms/blas/BatchedBlas.h @@ -28,6 +28,7 @@ Author: Peter Boyle #pragma once #ifdef GRID_HIP +#include #include #endif #ifdef GRID_CUDA @@ -255,16 +256,29 @@ public: if ( OpB == GridBLAS_OP_N ) hOpB = HIPBLAS_OP_N; if ( OpB == GridBLAS_OP_T ) hOpB = HIPBLAS_OP_T; if ( OpB == GridBLAS_OP_C ) hOpB = HIPBLAS_OP_C; +#if defined(HIP_VERSION_MAJOR) && (HIP_VERSION_MAJOR >=7) auto err = hipblasZgemmBatched(gridblasHandle, hOpA, hOpB, m,n,k, - (hipblasDoubleComplex *) &alpha_p[0], - (hipblasDoubleComplex **)&Amk[0], lda, - (hipblasDoubleComplex **)&Bkn[0], ldb, - (hipblasDoubleComplex *) &beta_p[0], - (hipblasDoubleComplex **)&Cmn[0], ldc, + (hipDoubleComplex *) &alpha_p[0], + (hipDoubleComplex **)&Amk[0], lda, + (hipDoubleComplex **)&Bkn[0], ldb, + (hipDoubleComplex *) &beta_p[0], + (hipDoubleComplex **)&Cmn[0], ldc, batchCount); +#else + auto err = hipblasZgemmBatched(gridblasHandle, + hOpA, + hOpB, + m,n,k, + (hipblasDoubleComplex *) &alpha_p[0], + (hipblasDoubleComplex **)&Amk[0], lda, + (hipblasDoubleComplex **)&Bkn[0], ldb, + (hipblasDoubleComplex *) &beta_p[0], + (hipblasDoubleComplex **)&Cmn[0], ldc, + batchCount); +#endif // std::cout << " hipblas return code " <<(int)err<=7) auto err = hipblasCgemmBatched(gridblasHandle, hOpA, hOpB, m,n,k, - (hipblasComplex *) &alpha_p[0], - (hipblasComplex **)&Amk[0], lda, - (hipblasComplex **)&Bkn[0], ldb, - (hipblasComplex *) &beta_p[0], - (hipblasComplex **)&Cmn[0], ldc, + (hipComplex *) &alpha_p[0], + (hipComplex **)&Amk[0], lda, + (hipComplex **)&Bkn[0], ldb, + (hipComplex *) &beta_p[0], + (hipComplex **)&Cmn[0], ldc, batchCount); +#else + auto err = hipblasCgemmBatched(gridblasHandle, + hOpA, + hOpB, + m,n,k, + (hipblasComplex *) &alpha_p[0], + (hipblasComplex **)&Amk[0], lda, + (hipblasComplex **)&Bkn[0], ldb, + (hipblasComplex *) &beta_p[0], + (hipblasComplex **)&Cmn[0], ldc, + batchCount); +#endif GRID_ASSERT(err==HIPBLAS_STATUS_SUCCESS); #endif #ifdef GRID_CUDA @@ -1094,11 +1121,19 @@ public: GRID_ASSERT(info.size()==batchCount); #ifdef GRID_HIP +#if defined(HIP_VERSION_MAJOR) && (HIP_VERSION_MAJOR >=7) auto err = hipblasZgetrfBatched(gridblasHandle,(int)n, - (hipblasDoubleComplex **)&Ann[0], (int)n, + (hipDoubleComplex **)&Ann[0], (int)n, (int*) &ipiv[0], (int*) &info[0], (int)batchCount); +#else + auto err = hipblasZgetrfBatched(gridblasHandle,(int)n, + (hipblasDoubleComplex **)&Ann[0], (int)n, + (int*) &ipiv[0], + (int*) &info[0], + (int)batchCount); +#endif GRID_ASSERT(err==HIPBLAS_STATUS_SUCCESS); #endif #ifdef GRID_CUDA @@ -1124,11 +1159,20 @@ public: GRID_ASSERT(info.size()==batchCount); #ifdef GRID_HIP +#if defined(HIP_VERSION_MAJOR) && (HIP_VERSION_MAJOR >=7) auto err = hipblasCgetrfBatched(gridblasHandle,(int)n, - (hipblasComplex **)&Ann[0], (int)n, + (hipComplex **)&Ann[0], (int)n, (int*) &ipiv[0], (int*) &info[0], (int)batchCount); +#else + auto err = hipblasCgetrfBatched(gridblasHandle,(int)n, + (hipblasComplex **)&Ann[0], (int)n, + (int*) &ipiv[0], + (int*) &info[0], + (int)batchCount); +#endif + GRID_ASSERT(err==HIPBLAS_STATUS_SUCCESS); #endif #ifdef GRID_CUDA @@ -1201,12 +1245,22 @@ public: GRID_ASSERT(Cnn.size()==batchCount); #ifdef GRID_HIP +#if defined(HIP_VERSION_MAJOR) && (HIP_VERSION_MAJOR >=7) auto err = hipblasZgetriBatched(gridblasHandle,(int)n, - (hipblasDoubleComplex **)&Ann[0], (int)n, + (hipDoubleComplex **)&Ann[0], (int)n, (int*) &ipiv[0], - (hipblasDoubleComplex **)&Cnn[0], (int)n, + (hipDoubleComplex **)&Cnn[0], (int)n, (int*) &info[0], (int)batchCount); +#else + auto err = hipblasZgetriBatched(gridblasHandle,(int)n, + (hipblasDoubleComplex **)&Ann[0], (int)n, + (int*) &ipiv[0], + (hipblasDoubleComplex **)&Cnn[0], (int)n, + (int*) &info[0], + (int)batchCount); + +#endif GRID_ASSERT(err==HIPBLAS_STATUS_SUCCESS); #endif #ifdef GRID_CUDA @@ -1235,12 +1289,21 @@ public: GRID_ASSERT(Cnn.size()==batchCount); #ifdef GRID_HIP +#if defined(HIP_VERSION_MAJOR) && (HIP_VERSION_MAJOR >=7) auto err = hipblasCgetriBatched(gridblasHandle,(int)n, - (hipblasComplex **)&Ann[0], (int)n, + (hipComplex **)&Ann[0], (int)n, (int*) &ipiv[0], - (hipblasComplex **)&Cnn[0], (int)n, + (hipComplex **)&Cnn[0], (int)n, (int*) &info[0], (int)batchCount); +#else + auto err = hipblasCgetriBatched(gridblasHandle,(int)n, + (hipblasComplex **)&Ann[0], (int)n, + (int*) &ipiv[0], + (hipblasComplex **)&Cnn[0], (int)n, + (int*) &info[0], + (int)batchCount); +#endif GRID_ASSERT(err==HIPBLAS_STATUS_SUCCESS); #endif #ifdef GRID_CUDA diff --git a/Grid/algorithms/iterative/AdefMrhs.h b/Grid/algorithms/iterative/AdefMrhs.h index 48034d83..e2090009 100644 --- a/Grid/algorithms/iterative/AdefMrhs.h +++ b/Grid/algorithms/iterative/AdefMrhs.h @@ -92,8 +92,8 @@ class TwoLevelCGmrhs // Vector case virtual void operator() (std::vector &src, std::vector &x) { - // SolveSingleSystem(src,x); - SolvePrecBlockCG(src,x); + SolveSingleSystem(src,x); + // SolvePrecBlockCG(src,x); } //////////////////////////////////////////////////////////////////////////////////////////////////// diff --git a/Grid/algorithms/multigrid/Aggregates.h b/Grid/algorithms/multigrid/Aggregates.h index 2fe7d276..fc527929 100644 --- a/Grid/algorithms/multigrid/Aggregates.h +++ b/Grid/algorithms/multigrid/Aggregates.h @@ -131,7 +131,10 @@ public: RealD scale; TrivialPrecon simple_fine; - PrecGeneralisedConjugateResidualNonHermitian GCR(0.001,10,DiracOp,simple_fine,30,30); + // PrecGeneralisedConjugateResidualNonHermitian GCR(0.001,10,DiracOp,simple_fine,30,30); + // PrecGeneralisedConjugateResidualNonHermitian GCR(0.001,10,DiracOp,simple_fine,12,12); + // PrecGeneralisedConjugateResidualNonHermitian GCR(0.001,30,DiracOp,simple_fine,12,12); + PrecGeneralisedConjugateResidualNonHermitian GCR(0.001,30,DiracOp,simple_fine,10,10); FineField noise(FineGrid); FineField src(FineGrid); FineField guess(FineGrid); diff --git a/Grid/qcd/QCD.h b/Grid/qcd/QCD.h index dbedfa7c..5b35dc64 100644 --- a/Grid/qcd/QCD.h +++ b/Grid/qcd/QCD.h @@ -596,16 +596,32 @@ template inline vobj transposeColour(const vobj &lhs){ ////////////////////////////////////////// // Trace lattice and non-lattice ////////////////////////////////////////// +#define GRID_UNOP(name) name +#define GRID_DEF_UNOP(op, name) \ + template ::value||is_lattice_expr::value,T1>::type * = nullptr> \ + inline auto op(const T1 &arg) ->decltype(LatticeUnaryExpression(GRID_UNOP(name)(), arg)) \ + { \ + return LatticeUnaryExpression(GRID_UNOP(name)(), arg); \ + } + template inline auto traceSpin(const Lattice &lhs) -> Lattice(vobj()))> { return traceIndex(lhs); } + +GridUnopClass(UnaryTraceSpin, traceIndex(a)); +GRID_DEF_UNOP(traceSpin, UnaryTraceSpin); + template inline auto traceColour(const Lattice &lhs) -> Lattice(vobj()))> { return traceIndex(lhs); } + +GridUnopClass(UnaryTraceColour, traceIndex(a)); +GRID_DEF_UNOP(traceColour, UnaryTraceColour); + template inline auto traceSpin(const vobj &lhs) -> Lattice(lhs))> { @@ -617,6 +633,8 @@ inline auto traceColour(const vobj &lhs) -> Lattice(lhs); } +#undef GRID_UNOP +#undef GRID_DEF_UNOP ////////////////////////////////////////// // Current types ////////////////////////////////////////// diff --git a/Grid/qcd/modules/ObservableModules.h b/Grid/qcd/modules/ObservableModules.h index 87fcbb92..cda86185 100644 --- a/Grid/qcd/modules/ObservableModules.h +++ b/Grid/qcd/modules/ObservableModules.h @@ -103,6 +103,18 @@ class PolyakovMod: public ObservableModule, NoParameters>{ PolyakovMod(): ObsBase(NoParameters()){} }; +template < class Impl > +class SpatialPolyakovMod: public ObservableModule, NoParameters>{ + typedef ObservableModule, NoParameters> ObsBase; + using ObsBase::ObsBase; // for constructors + + // acquire resource + virtual void initialize(){ + this->ObservablePtr.reset(new SpatialPolyakovLogger()); + } + public: + SpatialPolyakovMod(): ObsBase(NoParameters()){} +}; template < class Impl > class TopologicalChargeMod: public ObservableModule, TopologyObsParameters>{ diff --git a/Grid/qcd/observables/polyakov_loop.h b/Grid/qcd/observables/polyakov_loop.h index 0b59f549..57812ff6 100644 --- a/Grid/qcd/observables/polyakov_loop.h +++ b/Grid/qcd/observables/polyakov_loop.h @@ -2,11 +2,12 @@ Grid physics library, www.github.com/paboyle/Grid -Source file: ./lib/qcd/modules/polyakov_line.h +Source file: ./Grid/qcd/observables/polyakov_loop.h -Copyright (C) 2017 +Copyright (C) 2025 Author: David Preti +Author: Alexis Verney-Provatas <2414441@swansea.ac.uk> This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -60,4 +61,43 @@ class PolyakovLogger : public HmcObservable { } }; +template +class SpatialPolyakovLogger : public HmcObservable { + public: + // here forces the Impl to be of gauge fields + // if not the compiler will complain + INHERIT_GIMPL_TYPES(Impl); + + // necessary for HmcObservable compatibility + typedef typename Impl::Field Field; + + void TrajectoryComplete(int traj, + Field &U, + GridSerialRNG &sRNG, + GridParallelRNG &pRNG) { + + // Save current numerical output precision + int def_prec = std::cout.precision(); + + // Assume that the dimensions are D=3+1 + int Ndim = 3; + ComplexD polyakov; + + // Iterate over the spatial directions and print the average spatial polyakov loop + // over them + for (int idx=0; idx::avgPolyakovLoop(U, idx); + + std::cout << GridLogMessage + << std::setprecision(std::numeric_limits::digits10 + 1) + << "Polyakov Loop in the " << idx << " spatial direction : [ " << traj << " ] "<< polyakov << std::endl; + + } + + // Return to original output precision + std::cout.precision(def_prec); + + } +}; + NAMESPACE_END(Grid); diff --git a/Grid/qcd/utils/Sp2n.impl.h b/Grid/qcd/utils/Sp2n.impl.h index 48b7c744..5e57c491 100644 --- a/Grid/qcd/utils/Sp2n.impl.h +++ b/Grid/qcd/utils/Sp2n.impl.h @@ -254,9 +254,9 @@ static void testGenerators(GroupName::Sp) { } } -template -static Lattice > > > -ProjectOnGeneralGroup(const Lattice > > > &Umu, GroupName::Sp) { +template +static Lattice > > > +ProjectOnGeneralGroup(const Lattice > > > &Umu, GroupName::Sp) { return ProjectOnSpGroup(Umu); } diff --git a/Grid/qcd/utils/WilsonLoops.h b/Grid/qcd/utils/WilsonLoops.h index 23ad237c..8091cbc8 100644 --- a/Grid/qcd/utils/WilsonLoops.h +++ b/Grid/qcd/utils/WilsonLoops.h @@ -177,25 +177,43 @@ public: } ////////////////////////////////////////////////// - // average over all x,y,z the temporal loop + // average Polyakov loop in mu direction over all directions != mu ////////////////////////////////////////////////// - static ComplexD avgPolyakovLoop(const GaugeField &Umu) { //assume Nd=4 - GaugeMat Ut(Umu.Grid()), P(Umu.Grid()); + static ComplexD avgPolyakovLoop(const GaugeField &Umu, const int mu) { //assume Nd=4 + + // Protect against bad value of mu [0, 3] + if ((mu < 0 ) || (mu > 3)) { + std::cout << GridLogError << "Index is not an integer inclusively between 0 and 3." << std::endl; + exit(1); + } + + // U_loop is U_{mu} + GaugeMat U_loop(Umu.Grid()), P(Umu.Grid()); ComplexD out; int T = Umu.Grid()->GlobalDimensions()[3]; int X = Umu.Grid()->GlobalDimensions()[0]; int Y = Umu.Grid()->GlobalDimensions()[1]; int Z = Umu.Grid()->GlobalDimensions()[2]; - Ut = peekLorentz(Umu,3); //Select temporal direction - P = Ut; - for (int t=1;tGlobalDimensions()[mu]; + + U_loop = peekLorentz(Umu, mu); //Select direction + P = U_loop; + for (int t=1;t + +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + +#include #if Nc == 3 #include @@ -230,3 +234,4 @@ int main(int argc, char **argv) #endif } // main +#endif diff --git a/HMC/FTHMC2p1f_3GeV.cc b/HMC/FTHMC2p1f_3GeV.cc index 36d5caa3..db2b937b 100644 --- a/HMC/FTHMC2p1f_3GeV.cc +++ b/HMC/FTHMC2p1f_3GeV.cc @@ -25,7 +25,11 @@ directory *************************************************************************************/ /* END LEGAL */ -#include + +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + +#include #if Nc == 3 #include @@ -231,5 +235,4 @@ int main(int argc, char **argv) #endif } // main - - +#endif diff --git a/HMC/HMC2p1f_3GeV.cc b/HMC/HMC2p1f_3GeV.cc index 199d4be8..d97b03c6 100644 --- a/HMC/HMC2p1f_3GeV.cc +++ b/HMC/HMC2p1f_3GeV.cc @@ -24,7 +24,11 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#include + +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + +#include #if Nc == 3 #include @@ -230,5 +234,4 @@ int main(int argc, char **argv) #endif } // main - - +#endif diff --git a/HMC/Mobius2p1f.cc b/HMC/Mobius2p1f.cc index 8042d6e6..b6b9da06 100644 --- a/HMC/Mobius2p1f.cc +++ b/HMC/Mobius2p1f.cc @@ -27,7 +27,11 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#include + +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + +#include int main(int argc, char **argv) { using namespace Grid; @@ -195,5 +199,4 @@ int main(int argc, char **argv) { Grid_finalize(); } // main - - +#endif diff --git a/HMC/Mobius2p1fEOFA.cc b/HMC/Mobius2p1fEOFA.cc index e7b0e473..ce7518d9 100644 --- a/HMC/Mobius2p1fEOFA.cc +++ b/HMC/Mobius2p1fEOFA.cc @@ -28,7 +28,11 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#include + +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + +#include #ifdef GRID_DEFAULT_PRECISION_DOUBLE #define MIXED_PRECISION @@ -449,5 +453,4 @@ int main(int argc, char **argv) { Grid_finalize(); } // main - - +#endif diff --git a/HMC/Mobius2p1fEOFA_F1.cc b/HMC/Mobius2p1fEOFA_F1.cc index d9fec4fa..e9f2bb3b 100644 --- a/HMC/Mobius2p1fEOFA_F1.cc +++ b/HMC/Mobius2p1fEOFA_F1.cc @@ -28,7 +28,11 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#include + +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + +#include #ifdef GRID_DEFAULT_PRECISION_DOUBLE #define MIXED_PRECISION @@ -442,5 +446,4 @@ int main(int argc, char **argv) { Grid_finalize(); } // main - - +#endif diff --git a/HMC/Mobius2p1fIDSDRGparityEOFA_40ID.cc b/HMC/Mobius2p1fIDSDRGparityEOFA_40ID.cc index 2825d403..3fe94313 100644 --- a/HMC/Mobius2p1fIDSDRGparityEOFA_40ID.cc +++ b/HMC/Mobius2p1fIDSDRGparityEOFA_40ID.cc @@ -28,7 +28,11 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#include + +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + +#include using namespace Grid; @@ -918,3 +922,5 @@ int main(int argc, char **argv) { return 0; #endif } // main + +#endif diff --git a/HMC/Mobius2p1fIDSDRGparityEOFA_48ID.cc b/HMC/Mobius2p1fIDSDRGparityEOFA_48ID.cc index 4c1b5856..d8a50f6b 100644 --- a/HMC/Mobius2p1fIDSDRGparityEOFA_48ID.cc +++ b/HMC/Mobius2p1fIDSDRGparityEOFA_48ID.cc @@ -28,7 +28,11 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#include + +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + +#include using namespace Grid; @@ -873,3 +877,5 @@ int main(int argc, char **argv) { return 0; #endif } // main + +#endif diff --git a/HMC/Mobius2p1fRHMC.cc b/HMC/Mobius2p1fRHMC.cc index 288a6c54..c49995b8 100644 --- a/HMC/Mobius2p1fRHMC.cc +++ b/HMC/Mobius2p1fRHMC.cc @@ -27,7 +27,11 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#include + +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + +#include int main(int argc, char **argv) { using namespace Grid; @@ -193,5 +197,4 @@ int main(int argc, char **argv) { Grid_finalize(); } // main - - +#endif diff --git a/HMC/Mobius2p1f_DD_EOFA_96I_3level.cc b/HMC/Mobius2p1f_DD_EOFA_96I_3level.cc index c0b5ced1..1aa22332 100644 --- a/HMC/Mobius2p1f_DD_EOFA_96I_3level.cc +++ b/HMC/Mobius2p1f_DD_EOFA_96I_3level.cc @@ -27,7 +27,11 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#include + +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + +#include NAMESPACE_BEGIN(Grid); @@ -512,5 +516,4 @@ int main(int argc, char **argv) { Grid_finalize(); } // main - - +#endif diff --git a/HMC/Mobius2p1f_DD_EOFA_96I_double.cc b/HMC/Mobius2p1f_DD_EOFA_96I_double.cc index 678df981..1c477cb2 100644 --- a/HMC/Mobius2p1f_DD_EOFA_96I_double.cc +++ b/HMC/Mobius2p1f_DD_EOFA_96I_double.cc @@ -27,7 +27,11 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#include + +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + +#include int main(int argc, char **argv) { using namespace Grid; @@ -345,5 +349,4 @@ int main(int argc, char **argv) { Grid_finalize(); } // main - - +#endif diff --git a/HMC/Mobius2p1f_DD_EOFA_96I_mixed.cc b/HMC/Mobius2p1f_DD_EOFA_96I_mixed.cc index becfac53..44a153d2 100644 --- a/HMC/Mobius2p1f_DD_EOFA_96I_mixed.cc +++ b/HMC/Mobius2p1f_DD_EOFA_96I_mixed.cc @@ -27,7 +27,11 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#include + +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + +#include NAMESPACE_BEGIN(Grid); @@ -516,5 +520,4 @@ int main(int argc, char **argv) { Grid_finalize(); } // main - - +#endif diff --git a/HMC/Mobius2p1f_DD_EOFA_96I_mshift.cc b/HMC/Mobius2p1f_DD_EOFA_96I_mshift.cc index 7ae96dcd..d2dd3b9e 100644 --- a/HMC/Mobius2p1f_DD_EOFA_96I_mshift.cc +++ b/HMC/Mobius2p1f_DD_EOFA_96I_mshift.cc @@ -27,7 +27,11 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#include + +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + +#include NAMESPACE_BEGIN(Grid); @@ -567,5 +571,4 @@ int main(int argc, char **argv) { Grid_finalize(); } // main - - +#endif diff --git a/HMC/Mobius2p1f_DD_RHMC.cc b/HMC/Mobius2p1f_DD_RHMC.cc index 39b4c1dd..a95b5b2e 100644 --- a/HMC/Mobius2p1f_DD_RHMC.cc +++ b/HMC/Mobius2p1f_DD_RHMC.cc @@ -27,7 +27,11 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#include + +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + +#include int main(int argc, char **argv) { using namespace Grid; @@ -263,5 +267,4 @@ int main(int argc, char **argv) { Grid_finalize(); } // main - - +#endif diff --git a/HMC/Mobius2p1f_DD_RHMC_96I.cc b/HMC/Mobius2p1f_DD_RHMC_96I.cc index c28296a3..d2d10d2d 100644 --- a/HMC/Mobius2p1f_DD_RHMC_96I.cc +++ b/HMC/Mobius2p1f_DD_RHMC_96I.cc @@ -27,7 +27,11 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#include + +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + +#include int main(int argc, char **argv) { using namespace Grid; @@ -417,5 +421,4 @@ int main(int argc, char **argv) { Grid_finalize(); } // main - - +#endif diff --git a/HMC/Mobius2p1f_DD_RHMC_96I_mixed.cc b/HMC/Mobius2p1f_DD_RHMC_96I_mixed.cc index 864cefcc..5a1bdb15 100644 --- a/HMC/Mobius2p1f_DD_RHMC_96I_mixed.cc +++ b/HMC/Mobius2p1f_DD_RHMC_96I_mixed.cc @@ -27,7 +27,11 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#include + +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + +#include NAMESPACE_BEGIN(Grid); @@ -452,5 +456,4 @@ int main(int argc, char **argv) { Grid_finalize(); } // main - - +#endif diff --git a/HMC/Mobius2p1f_EOFA_96I_hmc.cc b/HMC/Mobius2p1f_EOFA_96I_hmc.cc index d967d55a..b4999fc3 100644 --- a/HMC/Mobius2p1f_EOFA_96I_hmc.cc +++ b/HMC/Mobius2p1f_EOFA_96I_hmc.cc @@ -27,7 +27,11 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#include + +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + +#include NAMESPACE_BEGIN(Grid); @@ -462,5 +466,4 @@ int main(int argc, char **argv) { Grid_finalize(); } // main - - +#endif diff --git a/HMC/Mobius2p1f_EOFA_96I_hmc_double.cc b/HMC/Mobius2p1f_EOFA_96I_hmc_double.cc index 61b32b5c..201866a0 100644 --- a/HMC/Mobius2p1f_EOFA_96I_hmc_double.cc +++ b/HMC/Mobius2p1f_EOFA_96I_hmc_double.cc @@ -27,7 +27,11 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ -#include + +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + +#include @@ -264,5 +268,4 @@ int main(int argc, char **argv) { Grid_finalize(); } // main - - +#endif diff --git a/HMC/disable_examples_without_instantiations.h b/HMC/disable_examples_without_instantiations.h new file mode 100644 index 00000000..308793af --- /dev/null +++ b/HMC/disable_examples_without_instantiations.h @@ -0,0 +1,16 @@ +#include +#pragma once + + +#ifndef ENABLE_FERMION_INSTANTIATIONS +#include + +int main(void) { + std::cout << "This build of Grid was configured to exclude fermion instantiations, " + << "which this example relies on. " + << "Please reconfigure and rebuild Grid with --enable-fermion-instantiations" + << "to run this example." + << std::endl; + return 1; +} +#endif diff --git a/benchmarks/Benchmark_ITT.cc b/benchmarks/Benchmark_ITT.cc index 48980b2c..64255b99 100644 --- a/benchmarks/Benchmark_ITT.cc +++ b/benchmarks/Benchmark_ITT.cc @@ -26,6 +26,9 @@ Author: paboyle See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ +#include "disable_benchmarks_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + #include using namespace Grid; @@ -731,3 +734,5 @@ int main (int argc, char ** argv) Grid_finalize(); } + +#endif diff --git a/benchmarks/Benchmark_dwf.cc b/benchmarks/Benchmark_dwf.cc index d7403705..c28b2686 100644 --- a/benchmarks/Benchmark_dwf.cc +++ b/benchmarks/Benchmark_dwf.cc @@ -20,6 +20,9 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ +#include "disable_benchmarks_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + #include #ifdef GRID_CUDA #define CUDA_PROFILE @@ -439,3 +442,4 @@ void Benchmark(int Ls, Coordinate Dirichlet,bool sloppy) GRID_ASSERT(norm2(src_e)<1.0e-4); GRID_ASSERT(norm2(src_o)<1.0e-4); } +#endif diff --git a/benchmarks/Benchmark_dwf_fp32.cc b/benchmarks/Benchmark_dwf_fp32.cc index 03dce1df..f0927519 100644 --- a/benchmarks/Benchmark_dwf_fp32.cc +++ b/benchmarks/Benchmark_dwf_fp32.cc @@ -20,6 +20,10 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ +#include "disable_benchmarks_without_instantiations.h" + +#ifdef ENABLE_FERMION_INSTANTIATIONS + #include #ifdef GRID_CUDA #define CUDA_PROFILE @@ -439,3 +443,5 @@ void Benchmark(int Ls, Coordinate Dirichlet,bool sloppy) GRID_ASSERT(norm2(src_e)<1.0e-4); GRID_ASSERT(norm2(src_o)<1.0e-4); } + +#endif diff --git a/benchmarks/Benchmark_dwf_fp32_paranoid.cc b/benchmarks/Benchmark_dwf_fp32_paranoid.cc index e0090128..67317313 100644 --- a/benchmarks/Benchmark_dwf_fp32_paranoid.cc +++ b/benchmarks/Benchmark_dwf_fp32_paranoid.cc @@ -20,6 +20,9 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ +#include "disable_benchmarks_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + #include #ifdef GRID_CUDA #define CUDA_PROFILE @@ -385,3 +388,5 @@ int main (int argc, char ** argv) Grid_finalize(); exit(0); } + +#endif diff --git a/benchmarks/Benchmark_dwf_sweep.cc b/benchmarks/Benchmark_dwf_sweep.cc index 2f11eb22..a3ccc9e3 100644 --- a/benchmarks/Benchmark_dwf_sweep.cc +++ b/benchmarks/Benchmark_dwf_sweep.cc @@ -26,6 +26,9 @@ Author: paboyle See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ +#include "disable_benchmarks_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + #include using namespace std; @@ -238,5 +241,4 @@ void benchDw(std::vector & latt4, int Ls, int threads,int report ) } } - - +#endif diff --git a/benchmarks/Benchmark_gparity.cc b/benchmarks/Benchmark_gparity.cc index 421dd3cd..1cfa2127 100644 --- a/benchmarks/Benchmark_gparity.cc +++ b/benchmarks/Benchmark_gparity.cc @@ -1,3 +1,7 @@ +#include "disable_benchmarks_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + + #include #include using namespace std; @@ -155,3 +159,4 @@ int main (int argc, char ** argv) Grid_finalize(); } +#endif diff --git a/benchmarks/Benchmark_halo.cc b/benchmarks/Benchmark_halo.cc index 43138e67..f95c29ad 100644 --- a/benchmarks/Benchmark_halo.cc +++ b/benchmarks/Benchmark_halo.cc @@ -20,6 +20,9 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ +#include "disable_benchmarks_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + #include #ifdef GRID_CUDA #define CUDA_PROFILE @@ -129,3 +132,5 @@ int main (int argc, char ** argv) Grid_finalize(); exit(0); } + +#endif diff --git a/benchmarks/Benchmark_mooee.cc b/benchmarks/Benchmark_mooee.cc index 54235752..d3d39126 100644 --- a/benchmarks/Benchmark_mooee.cc +++ b/benchmarks/Benchmark_mooee.cc @@ -26,6 +26,9 @@ Author: paboyle See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ +#include "disable_benchmarks_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + #include using namespace std; @@ -149,3 +152,5 @@ int main (int argc, char ** argv) Grid_finalize(); } + +#endif diff --git a/benchmarks/Benchmark_schur.cc b/benchmarks/Benchmark_schur.cc index 8171998a..644a158c 100644 --- a/benchmarks/Benchmark_schur.cc +++ b/benchmarks/Benchmark_schur.cc @@ -26,6 +26,9 @@ Author: paboyle See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ +#include "disable_benchmarks_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + #include using namespace std; @@ -172,5 +175,4 @@ void benchDw(std::vector & latt4, int Ls) // Dw.Report(); } - - +#endif diff --git a/benchmarks/Benchmark_staggered.cc b/benchmarks/Benchmark_staggered.cc index a2be7f62..65e04f27 100644 --- a/benchmarks/Benchmark_staggered.cc +++ b/benchmarks/Benchmark_staggered.cc @@ -26,6 +26,9 @@ Author: paboyle See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ +#include "disable_benchmarks_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + #include using namespace std; @@ -110,3 +113,5 @@ int main (int argc, char ** argv) Grid_finalize(); } + +#endif diff --git a/benchmarks/Benchmark_staggeredF.cc b/benchmarks/Benchmark_staggeredF.cc index f7beed2d..e0f4331f 100644 --- a/benchmarks/Benchmark_staggeredF.cc +++ b/benchmarks/Benchmark_staggeredF.cc @@ -26,6 +26,9 @@ Author: paboyle See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ +#include "disable_benchmarks_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + #include using namespace std; @@ -112,3 +115,5 @@ int main (int argc, char ** argv) Grid_finalize(); } + +#endif diff --git a/benchmarks/Benchmark_usqcd.cc b/benchmarks/Benchmark_usqcd.cc index 70489419..34a0d6ec 100644 --- a/benchmarks/Benchmark_usqcd.cc +++ b/benchmarks/Benchmark_usqcd.cc @@ -26,6 +26,10 @@ Author: paboyle See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ +#include "disable_benchmarks_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + + #include #include @@ -978,3 +982,5 @@ int main (int argc, char ** argv) Grid_finalize(); fclose(FP); } + +#endif diff --git a/benchmarks/Benchmark_wilson.cc b/benchmarks/Benchmark_wilson.cc index 3572f6b1..ac2b8ebd 100644 --- a/benchmarks/Benchmark_wilson.cc +++ b/benchmarks/Benchmark_wilson.cc @@ -26,6 +26,9 @@ Author: paboyle See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ +#include "disable_benchmarks_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + #include using namespace std; @@ -258,3 +261,5 @@ int main (int argc, char ** argv) Grid_finalize(); } + +#endif diff --git a/benchmarks/Benchmark_wilson_sweep.cc b/benchmarks/Benchmark_wilson_sweep.cc index 45a10b25..6f840d56 100644 --- a/benchmarks/Benchmark_wilson_sweep.cc +++ b/benchmarks/Benchmark_wilson_sweep.cc @@ -19,6 +19,9 @@ Author: Richard Rollins See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ +#include "disable_benchmarks_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + #include using namespace std; @@ -161,3 +164,5 @@ void bench_wilson_eo ( double flops = (single_site_flops * volume * ncall)/2.0; std::cout << flops/(t1-t0) << "\t\t"; } + +#endif diff --git a/benchmarks/disable_benchmarks_without_instantiations.h b/benchmarks/disable_benchmarks_without_instantiations.h new file mode 100644 index 00000000..0277ccf5 --- /dev/null +++ b/benchmarks/disable_benchmarks_without_instantiations.h @@ -0,0 +1,16 @@ +#include + +#pragma once + +#ifndef ENABLE_FERMION_INSTANTIATIONS +#include + +int main(void) { + std::cout << "This build of Grid was configured to exclude fermion instantiations, " + << "which this benchmark relies on. " + << "Please reconfigure and rebuild Grid with --enable-fermion-instantiations" + << "to run this benchmark." + << std::endl; + return 1; +} +#endif diff --git a/configure.ac b/configure.ac index 55ceed3b..ab018787 100644 --- a/configure.ac +++ b/configure.ac @@ -172,6 +172,12 @@ case ${ac_TRACING} in esac ############### fermions +AC_ARG_ENABLE([fermion-instantiations], + [AS_HELP_STRING([--enable-fermion-instantiations=yes|no],[enable fermion instantiations])], + [ac_FERMION_REPS=${enable_fermion_instantiations}], [ac_FERMION_INSTANTIATIONS=yes]) + +AM_CONDITIONAL(BUILD_FERMION_INSTANTIATIONS, [ test "${ac_FERMION_INSTANTIATIONS}X" == "yesX" ]) + AC_ARG_ENABLE([fermion-reps], [AS_HELP_STRING([--enable-fermion-reps=yes|no],[enable extra fermion representation support])], [ac_FERMION_REPS=${enable_fermion_reps}], [ac_FERMION_REPS=yes]) @@ -194,6 +200,9 @@ AM_CONDITIONAL(BUILD_ZMOBIUS, [ test "${ac_ZMOBIUS}X" == "yesX" ]) case ${ac_FERMION_REPS} in yes) AC_DEFINE([ENABLE_FERMION_REPS],[1],[non QCD fermion reps]);; esac +case ${ac_FERMION_INSTANTIATIONS} in + yes) AC_DEFINE([ENABLE_FERMION_INSTANTIATIONS],[1],[enable fermions]);; +esac case ${ac_GPARITY} in yes) AC_DEFINE([ENABLE_GPARITY],[1],[fermion actions with GPARITY BCs]);; esac diff --git a/examples/Example_Mobius_spectrum.cc b/examples/Example_Mobius_spectrum.cc index 65b7b7de..51609c69 100644 --- a/examples/Example_Mobius_spectrum.cc +++ b/examples/Example_Mobius_spectrum.cc @@ -3,6 +3,9 @@ * without regression / tests being applied */ +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + #include using namespace std; @@ -310,5 +313,4 @@ int main (int argc, char ** argv) Grid_finalize(); } - - +#endif diff --git a/examples/Example_christoph.cc b/examples/Example_christoph.cc index 140491db..eb9a9041 100644 --- a/examples/Example_christoph.cc +++ b/examples/Example_christoph.cc @@ -3,6 +3,9 @@ * without regression / tests being applied */ +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + #include using namespace std; @@ -432,5 +435,4 @@ int main (int argc, char ** argv) Grid_finalize(); } - - +#endif diff --git a/examples/Example_wall_wall_3pt.cc b/examples/Example_wall_wall_3pt.cc index 52f9cf42..c60ed02b 100644 --- a/examples/Example_wall_wall_3pt.cc +++ b/examples/Example_wall_wall_3pt.cc @@ -3,6 +3,9 @@ * without regression / tests being applied */ +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + #include using namespace std; @@ -535,5 +538,4 @@ int main (int argc, char ** argv) Grid_finalize(); } - - +#endif diff --git a/examples/Example_wall_wall_spectrum.cc b/examples/Example_wall_wall_spectrum.cc index 97ae1e06..b8b1a0d7 100644 --- a/examples/Example_wall_wall_spectrum.cc +++ b/examples/Example_wall_wall_spectrum.cc @@ -3,6 +3,9 @@ * without regression / tests being applied */ +#include "disable_examples_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + #include using namespace std; @@ -429,5 +432,4 @@ int main (int argc, char ** argv) Grid_finalize(); } - - +#endif diff --git a/examples/disable_examples_without_instantiations.h b/examples/disable_examples_without_instantiations.h new file mode 100644 index 00000000..8d523557 --- /dev/null +++ b/examples/disable_examples_without_instantiations.h @@ -0,0 +1,15 @@ +#include +#pragma once + +#ifndef ENABLE_FERMION_INSTANTIATIONS +#include + +int main(void) { + std::cout << "This build of Grid was configured to exclude fermion instantiations, " + << "which this example relies on. " + << "Please reconfigure and rebuild Grid with --enable-fermion-instantiations" + << "to run this example." + << std::endl; + return 1; +} +#endif diff --git a/tests/Test_cayley_even_odd_vec.cc b/tests/Test_cayley_even_odd_vec.cc index 243d1f72..f0e1b4de 100644 --- a/tests/Test_cayley_even_odd_vec.cc +++ b/tests/Test_cayley_even_odd_vec.cc @@ -25,6 +25,9 @@ Author: Peter Boyle See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ +#include "disable_tests_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + #include using namespace std; @@ -273,8 +276,6 @@ void TestWhat(What & Ddwf, err = phi-chi; std::cout< * Reimplement the badly named "multigrid" lanczos as compressed Lanczos using the features * in Grid that were intended to be used to support blocked Aggregates, from */ +#include "disable_tests_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + #include #include #include @@ -256,3 +259,4 @@ int main (int argc, char ** argv) { Grid_finalize(); } +#endif diff --git a/tests/Test_dwf_dslash_repro.cc b/tests/Test_dwf_dslash_repro.cc index b0eac64b..57a18b00 100644 --- a/tests/Test_dwf_dslash_repro.cc +++ b/tests/Test_dwf_dslash_repro.cc @@ -25,6 +25,9 @@ Author: Peter Boyle See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ +#include "disable_tests_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + #include using namespace std; @@ -237,3 +240,5 @@ int main (int argc, char ** argv) Grid_finalize(); } + +#endif diff --git a/tests/Test_dwf_mixedcg_prec.cc b/tests/Test_dwf_mixedcg_prec.cc index c87908a3..e5771adc 100644 --- a/tests/Test_dwf_mixedcg_prec.cc +++ b/tests/Test_dwf_mixedcg_prec.cc @@ -25,6 +25,9 @@ Author: Peter Boyle See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ +#include "disable_tests_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + #include using namespace std; @@ -222,3 +225,5 @@ int main (int argc, char ** argv) Grid_finalize(); } + +#endif diff --git a/tests/Test_dwf_mixedcg_prec_halfcomms.cc b/tests/Test_dwf_mixedcg_prec_halfcomms.cc index ff52b0d1..7a59f77d 100644 --- a/tests/Test_dwf_mixedcg_prec_halfcomms.cc +++ b/tests/Test_dwf_mixedcg_prec_halfcomms.cc @@ -25,6 +25,9 @@ Author: Peter Boyle See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ +#include "disable_tests_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS + #include using namespace std; @@ -118,3 +121,4 @@ int main (int argc, char ** argv) Grid_finalize(); } #endif +#endif diff --git a/tests/Test_meson_field.cc b/tests/Test_meson_field.cc index fa428d6a..5f85047f 100644 --- a/tests/Test_meson_field.cc +++ b/tests/Test_meson_field.cc @@ -24,6 +24,8 @@ with this program; if not, write to the Free Software Foundation, Inc., See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ +#include "disable_tests_without_instantiations.h" +#ifdef ENABLE_FERMION_INSTANTIATIONS #include #include @@ -157,3 +159,5 @@ int main(int argc, char *argv[]) return EXIT_SUCCESS; } + +#endif diff --git a/tests/debug/Test_general_coarse_hdcg.cc b/tests/debug/Test_general_coarse_hdcg.cc index 8186b561..9c6be62e 100644 --- a/tests/debug/Test_general_coarse_hdcg.cc +++ b/tests/debug/Test_general_coarse_hdcg.cc @@ -128,6 +128,10 @@ int main (int argc, char ** argv) typedef HermOpAdaptor HermFineMatrix; HermFineMatrix FineHermOp(HermOpEO); + LatticeFermionD src(FrbGrid); + src = ComplexD(1.0); + PowerMethod PM; PM(HermOpEO,src); + //////////////////////////////////////////////////////////// ///////////// Coarse basis and Little Dirac Operator /////// //////////////////////////////////////////////////////////// @@ -150,7 +154,7 @@ int main (int argc, char ** argv) std::cout << "**************************************"< MrhsHermMatrix; - Chebyshev IRLCheby(0.05,40.0,101); // 1 iter + Chebyshev IRLCheby(0.01,16.0,201); // 1 iter MrhsHermMatrix MrhsCoarseOp (mrhs); CoarseVector pm_src(CoarseMrhs); @@ -193,10 +197,10 @@ int main (int argc, char ** argv) PowerMethod cPM; cPM(MrhsCoarseOp,pm_src); - int Nk=nrhs; - int Nm=Nk*3; - // int Nk=36; - // int Nm=144; + // int Nk=16; + // int Nm=Nk*3; + int Nk=32; + int Nm=128; int Nstop=Nk; int Nconv_test_interval=1; @@ -210,7 +214,7 @@ int main (int argc, char ** argv) nrhs, Nk, Nm, - 1e-4,10); + 1e-4,100); int Nconv; std::vector eval(Nm); @@ -231,8 +235,6 @@ int main (int argc, char ** argv) std::cout << "**************************************"< MrhsGuesser; MrhsGuesser.ImportEigenBasis(evec,eval); @@ -252,9 +254,11 @@ int main (int argc, char ** argv) // Extra HDCG parameters ////////////////////////// int maxit=3000; - ConjugateGradient CG(2.0e-1,maxit,false); - RealD lo=2.0; - int ord = 9; + // ConjugateGradient CG(2.0e-1,maxit,false); + // ConjugateGradient CG(1.0e-2,maxit,false); + ConjugateGradient CG(5.0e-2,maxit,false); + RealD lo=0.2; + int ord = 7; DoNothingGuesser DoNothing; HPDSolver HPDSolveMrhs(MrhsCoarseOp,CG,DoNothing); @@ -300,6 +304,19 @@ int main (int argc, char ** argv) ConjugateGradient CGfine(1.0e-8,30000,false); CGfine(HermOpEO, src, result); } + { + std::cout << "**************************************"< HermOp(Ddwf); + ConjugateGradient CGfine(1.0e-8,30000,false); + CGfine(HermOp, src, result); + } #endif Grid_finalize(); return 0; diff --git a/tests/debug/Test_general_coarse_pvdagm.cc b/tests/debug/Test_general_coarse_pvdagm.cc index d34a067b..d0ea894c 100644 --- a/tests/debug/Test_general_coarse_pvdagm.cc +++ b/tests/debug/Test_general_coarse_pvdagm.cc @@ -368,7 +368,10 @@ int main (int argc, char ** argv) TrivialPrecon simple; NonHermitianLinearOperator LinOpCoarse(LittleDiracOpPV); // PrecGeneralisedConjugateResidualNonHermitian L2PGCR(1.0e-4, 100, LinOpCoarse,simple,10,10); - PrecGeneralisedConjugateResidualNonHermitian L2PGCR(1.0e-2, 200, LinOpCoarse,simple,30,30); + // PrecGeneralisedConjugateResidualNonHermitian L2PGCR(3.0e-2, 100, LinOpCoarse,simple,12,12); // 35 outer + // PrecGeneralisedConjugateResidualNonHermitian L2PGCR(5.0e-2, 100, LinOpCoarse,simple,12,12); // 36 outer, 12s + // PrecGeneralisedConjugateResidualNonHermitian L2PGCR(1.0e-1, 100, LinOpCoarse,simple,12,12); // 36 ; 11s + PrecGeneralisedConjugateResidualNonHermitian L2PGCR(3.0e-1, 100, LinOpCoarse,simple,12,12); L2PGCR.Level(3); c_res=Zero(); L2PGCR(c_src,c_res); @@ -400,7 +403,7 @@ int main (int argc, char ** argv) LinOpCoarse, L2PGCR); - PrecGeneralisedConjugateResidualNonHermitian L1PGCR(1.0e-8,1000,PVdagM,TwoLevelPrecon,30,30); + PrecGeneralisedConjugateResidualNonHermitian L1PGCR(1.0e-8,100,PVdagM,TwoLevelPrecon,10,10); L1PGCR.Level(1); f_res=Zero(); diff --git a/tests/disable_tests_without_instantiations.h b/tests/disable_tests_without_instantiations.h new file mode 100644 index 00000000..1fd95db6 --- /dev/null +++ b/tests/disable_tests_without_instantiations.h @@ -0,0 +1,16 @@ +#include + +#pragma once + +#ifndef ENABLE_FERMION_INSTANTIATIONS +#include + +int main(void) { + std::cout << "This build of Grid was configured to exclude fermion instantiations, " + << "which this test relies on. " + << "Please reconfigure and rebuild Grid with --enable-fermion-instantiations" + << "to run this test." + << std::endl; + return 1; +} +#endif diff --git a/visualisation/CLAUDE.md b/visualisation/CLAUDE.md new file mode 100644 index 00000000..c4fb788f --- /dev/null +++ b/visualisation/CLAUDE.md @@ -0,0 +1,171 @@ +# CLAUDE.md + +This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository. + +## What This Is + +VTK-based visualisation and analysis tools for Grid lattice QCD eigenvector density and HMC force data. All programmes link against both Grid (for reading Scidac/ILDG lattice files) and VTK (for rendering). + +## Build + +```bash +cd /Users/peterboyle/QCD/AmSC/Grid/visualisation/build +cmake .. -DVTK_DIR=$HOME/QCD/vtk/VTK-9.4.2-install/lib/cmake/vtk-9.4 +make # e.g. make ControlledVisualise5D +``` + +All executables are built as macOS bundles (`.app`) except `ForceAnalysis`, `FindPeak`, and `DumpField`. + +## Programmes + +### ControlledVisualise5D + +Interactive VTK renderer for 5D DWF eigenvector density (`LatticeComplexD`). Driven via named pipe `/tmp/visualise_cmd`. + +**Launch script**: `/Volumes/X9Pro/visualisation/Grid/visualisation/build/Hdwf_1_long/visualise_controlled.sh` + +**Wire protocol** (one command per line to `/tmp/visualise_cmd`): + +| Command | Effect | +|---------|--------| +| `file ` / `file +N` / `file -N` | Jump to file by index or relative | +| `slice ` / `+N` / `-N` | Set or shift slice coordinate in dimension dim | +| `spin ` | Continuous azimuth rotation at deg/tick (100ms tick); `spin 0` stops | +| `azimuth ` | Single azimuth rotation step | +| `elevation ` | Single elevation rotation step | +| `zoom ` | Camera dolly (>1 = in) | +| `iso ` | Isosurface threshold in RMS units | +| `status` | Print current state | +| `quit` | Exit | + +**Dimension indices for 5D DWF grid** (`--grid 48.32.32.32.32`): + +| dim | axis | size | +|-----|------|------| +| 0 | s (Ls) | 48 | +| 1 | x | 32 | +| 2 | y | 32 | +| 3 | z | 32 | +| 4 | t | 32 | + +**MD time mapping** for trajectory 702 (241 files, τ=3.3–4.0): +- File index N → τ = 3.300000 + N × (1/480) +- τ → file index = round((τ − 3.3) × 480) + +**Display axes**: `--xyz 0.3.4` shows s, z, t. The `--slice` argument sets initial values for all dims; dims not in `--xyz`, `--sum`, or `--loop` are the fixed slice dimensions (x=dim1, y=dim2 with `--xyz 0.3.4`). + +**Spin**: Implemented via `g_spinDeg` global applied on every 100ms poll timer tick inside `CommandHandler::Execute()`. Does not flood the pipe. + +### FindPeak + +Reads a `LatticeComplexD` Scidac file, prints the top-N sites by real value to stderr. + +```bash +./FindPeak --grid 48.32.32.32.32 --mpi 1.1.1.1.1 2>peaks.txt +``` + +Key result: At τ=3.670833 the tunneling hotsite on the s=0 wall is (x=21, y=24, z=2, t=23). + +### ForceAnalysis + +Reads 4D `LatticeComplexD` force snapshot files (Shuhei's snapshots at `/Volumes/X9Pro/visualisation/Shuhei/snapshots/`). Outputs TSV of RMS and hotsite value per file to stderr. + +```bash +./ForceAnalysis --grid 32.32.32.32 --mpi 1.1.1 --hotsite 21.24.2.23 \ + 2>force.tsv 1>/dev/null +``` + +Force components: `Gauge_lat`, `Gauge_smr`, `Jacobian_smr`, `Ferm0047_lat`, `Ferm0047_smr`. + +### DumpField + +Reads a `LatticeComplexD` and dumps via Grid's `<<` operator to stdout for verification. + +### TranscriptToVideo + +Renders a conversation transcript to an MP4 video (1280×720, 10 fps) with a typewriter animation effect, scrolling history, and optional captions. Does **not** link against Grid — pure VTK only. + +#### Transcript format + +``` +[USER] First question text, possibly +continuing on the next line. + +A blank line within a turn creates a paragraph break (visual spacer). + +[ASSISTANT] Response text. +Multiple continuation lines are preserved +as separate display lines, not merged. + +[CAPTION] Caption text shown at bottom of screen in white italic. +[CAPTION] (whitespace-only body clears the caption) +[USER] Next question... +``` + +- Lines beginning `[USER]`, `[ASSISTANT]`, `[CAPTION]` start a new turn. +- Continuation lines (no `[TAG]` prefix) are joined with `\n` — each becomes its own wrapped display line. +- Blank lines within a turn become paragraph-break spacers. +- Markdown emphasis markers (`**`, `*`, `` ` ``) are stripped automatically. +- UTF-8 smart quotes, em-dashes, ellipses, arrows are transliterated to ASCII. + +#### Usage + +```bash +cd /Users/peterboyle/QCD/AmSC/Grid/visualisation/build +# Set runtime library paths first (see Runtime Environment below) +./TranscriptToVideo +``` + +Transcript files live in `/Users/peterboyle/QCD/AmSC/Grid/visualisation/` (e.g. `transcript`, `transcript2`, `transcript3`). + +#### Visual layout + +| Element | Detail | +|---------|--------| +| Background | Near-black navy `(0.04, 0.04, 0.10)` | +| `[USER]` text | Gold `(1.00, 0.84, 0.00)` | +| `[ASSISTANT]` text | Steel blue `(0.68, 0.85, 0.90)` | +| History | Up to 18 lines; brightness fades linearly from 0.85 (newest) to 0.20 (oldest) | +| Caption | Arial italic 20pt white with shadow, centred at bottom | +| Progress bar | Blue, top of frame | +| Typewriter speed | 50 chars/sec (5 chars/frame at 10 fps) | +| Pause between lines | 3 frames (0.3 s) | +| Word-wrap column | 60 chars (body only, after prefix) | + +#### Key implementation notes + +- **Persistent render context**: a single `vtkRenderWindow` is created once and reused for all frames. Creating a new window per frame exhausts the macOS Metal GPU command buffer after ~33 frames (`MTLCommandBufferErrorDomain Code=8`). +- **`SanitiseASCII()`**: replaces multi-byte UTF-8 sequences before passing to VTK's font renderer (which crashes on non-ASCII input). +- Output format is MP4 via `vtkFFMPEGWriter`. `SetOffScreenRendering(1)` is required for headless rendering. + +## Runtime Environment + +All executables in `build/` require Spack-installed HDF5/FFTW/GMP/MPFR on the dynamic linker path: + +```bash +SPACK=/Users/peterboyle/QCD/Spack/spack/opt/spack/darwin-m1 +export DYLD_LIBRARY_PATH=\ +$SPACK/hdf5-1.14.6-2265ms4kymgw6hcnwi6vqehslyfv74t4/lib:\ +$SPACK/fftw-3.3.10-aznn6h3nac5cycidlhrhgjxvntpcbg57/lib:\ +$SPACK/gmp-6.3.0-cwiz4n7ww33fnb3aban2iup4orcr6c7i/lib:\ +$SPACK/mpfr-4.2.1-exgbz4qshmet6tmmuttdewdlunfvtrlb/lib:\ +$DYLD_LIBRARY_PATH +``` + +(These paths are also set by the ControlledVisualise5D launch script.) + +## Key Physics Context + +See `/Volumes/X9Pro/visualisation/analysis_notes_20260407.md` for full analysis. Summary: + +- Near-zero mode of H_DWF localises on the two walls (s=0 and s=47) of the 5D domain wall geometry +- Topology change transfers norm between walls, mediated by a near-zero mode of H_w (Hermitian Wilson at m=−1.8) +- Tunneling hotsite on s=0 wall: (x=21, y=24, z=2, t=23); s=47 wall: (x=4, y=8, z=0, t=20) +- Light fermion pseudofermion force (Ferm0047_smr) peaks at ~20× RMS at the hotsite during tunneling — this is the restoring force that causes topological bounces + +## Grid/VTK interaction notes + +- Grid log messages go to stdout; all data output in analysis programmes uses stderr to avoid interleaving +- `TensorRemove()` is required when extracting a scalar from `peekSite()` result: `real(TensorRemove(peekSite(field, site)))` +- For runtime-determined grid dimensionality use `GridDefaultSimd(latt_size.size(), vComplex::Nsimd())` +- DYLD_LIBRARY_PATH must include Spack HDF5/FFTW/GMP/MPFR paths (see launch script) diff --git a/visualisation/ControlledVisualise5D.cxx b/visualisation/ControlledVisualise5D.cxx new file mode 100644 index 00000000..2ecb04f3 --- /dev/null +++ b/visualisation/ControlledVisualise5D.cxx @@ -0,0 +1,891 @@ +// ControlledVisualise5D.cxx +// Derived from Visualise5D.cxx by Peter Boyle +// +// A minimal-protocol rendering engine for 5D DWF eigenvector-density data. +// Intended to be driven by an external intelligent controller (e.g. Claude) +// that handles all natural-language interpretation and state tracking. +// +// Commands are sent one per line to the named pipe /tmp/visualise_cmd. +// State is reported to stdout after every command. +// +// Wire protocol (all fields whitespace-separated): +// +// slice set Slice[dim] = N (0-based, wraps to lattice size) +// slice + increment Slice[dim] by N +// slice - decrement Slice[dim] by N +// zoom camera Dolly by factor (>1 = in, <1 = out) +// iso set isosurface threshold to x RMS +// file jump to file by absolute index +// file + advance N files +// file - go back N files +// render force a render with current state +// status print current state to stdout +// quit exit cleanly +// +// Dimension indices for 5D DWF grid (e.g. --grid 48.32.32.32.32): +// s=0 (Ls) x=1 y=2 z=3 t=4 +// For a 4D grid (--grid 32.32.32.32): +// x=0 y=1 z=2 t=3 + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#define MPEG +#ifdef MPEG +#include +#endif + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#define USE_FLYING_EDGES +#ifdef USE_FLYING_EDGES +#include +typedef vtkFlyingEdges3D isosurface; +#else +#include +typedef vtkMarchingCubes isosurface; +#endif + +#define CMD_PIPE "/tmp/visualise_cmd" + +static int g_mpeg = 0; +static int g_framerate = 10; + +// ─── Thread-safe command queue ──────────────────────────────────────────────── +static std::queue g_cmdQueue; +static std::mutex g_cmdMutex; +static std::atomic g_running{true}; +static double g_spinDeg = 0.0; // degrees per poll tick; 0 = stopped + +// ─── MPEG recording state ───────────────────────────────────────────────────── +static bool g_recording = false; +static vtkFFMPEGWriter* g_mpegWriter = nullptr; +static vtkWindowToImageFilter* g_imageFilter = nullptr; +static std::string g_recordingFile; // AVI filename for mux step + +// ─── Audio state (PCM audio track, synced to video frames) ─────────────────── +static const int AUDIO_RATE = 44100; +static const double BEEP_FREQ = 800.0; +static const int BEEP_SAMPLES = AUDIO_RATE * 4 / 100; // 40ms beep +static std::vector g_audioBuffer; +static int g_beepRemaining = 0; +static double g_beepPhase = 0.0; +static int g_samplesPerFrame = AUDIO_RATE / 10; // updated at record start + +// Write one video frame worth of audio samples (beep or silence) to the buffer. +static void GenerateAudioFrame() +{ + for (int i = 0; i < g_samplesPerFrame; i++) { + int16_t s = 0; + if (g_beepRemaining > 0) { + int pos = BEEP_SAMPLES - g_beepRemaining; + double env = 1.0; + int fade = AUDIO_RATE / 100; // 10ms fade + if (pos < fade) env = (double)pos / fade; + if (g_beepRemaining < fade) env = (double)g_beepRemaining / fade; + s = (int16_t)(16000.0 * env * std::sin(2.0 * M_PI * BEEP_FREQ * g_beepPhase / AUDIO_RATE)); + g_beepPhase += 1.0; + --g_beepRemaining; + } else { + g_beepPhase = 0.0; + } + g_audioBuffer.push_back(s); + } +} + +static void TriggerBeep() { g_beepRemaining = BEEP_SAMPLES; } + +// Simple mono 16-bit PCM WAV writer. +static void WriteWAV(const std::string& path, const std::vector& buf, int rate) +{ + std::ofstream f(path, std::ios::binary); + int dataBytes = (int)(buf.size() * 2); + int chunkSize = 36 + dataBytes; + int byteRate = rate * 2; + f.write("RIFF", 4); f.write((char*)&chunkSize, 4); + f.write("WAVE", 4); + f.write("fmt ", 4); + int fmtSz = 16; f.write((char*)&fmtSz, 4); + int16_t pcm = 1; f.write((char*)&pcm, 2); + int16_t ch = 1; f.write((char*)&ch, 2); + f.write((char*)&rate, 4); + f.write((char*)&byteRate, 4); + int16_t blk = 2; f.write((char*)&blk, 2); + int16_t bps = 16; f.write((char*)&bps, 2); + f.write("data", 4); f.write((char*)&dataBytes, 4); + f.write((char*)buf.data(), dataBytes); +} + +// Play a short audible beep on the local machine (non-blocking). +static void PlayBeepAudible() +{ + system("afplay /System/Library/Sounds/Tink.aiff -v 0.4 &"); +} + +// ─── Grid I/O ───────────────────────────────────────────────────────────────── +template +void readFile(T& out, const std::string& fname) +{ + Grid::emptyUserRecord record; + Grid::ScidacReader RD; + RD.open(fname); + RD.readScidacFieldRecord(out, record); + RD.close(); +} + +using namespace Grid; + +// ─── Command reader thread ──────────────────────────────────────────────────── +void CommandReaderThread() +{ + mkfifo(CMD_PIPE, 0666); + std::cout << "[cmd] Listening on " << CMD_PIPE << std::endl; + + while (g_running) { + int fd = open(CMD_PIPE, O_RDONLY | O_NONBLOCK); + if (fd < 0) { usleep(200000); continue; } + + int flags = fcntl(fd, F_GETFL); + fcntl(fd, F_SETFL, flags & ~O_NONBLOCK); + + char buf[4096]; + std::string partial; + ssize_t n; + while (g_running && (n = read(fd, buf, sizeof(buf) - 1)) > 0) { + buf[n] = '\0'; + partial += buf; + size_t pos; + while ((pos = partial.find('\n')) != std::string::npos) { + std::string line = partial.substr(0, pos); + if (!line.empty() && line.back() == '\r') line.pop_back(); + if (!line.empty()) { + std::lock_guard lk(g_cmdMutex); + g_cmdQueue.push(line); + } + partial = partial.substr(pos + 1); + } + } + close(fd); + } +} + +// ─── FrameUpdater ───────────────────────────────────────────────────────────── +class FrameUpdater : public vtkCallbackCommand +{ +public: + FrameUpdater() : ffile(0), TimerCount(0), old_file(-1), timerId(-2), maxCount(-1) {} + + static FrameUpdater* New() { return new FrameUpdater; } + + int ffile; + int old_file; + int timerId; + int maxCount; + + Coordinate latt; + Coordinate xyz_dims, xyz_ranges, g_xyz_ranges; + uint64_t xyz_vol; + Coordinate loop_dims, loop_ranges; + uint64_t loop_vol; + Coordinate sum_dims, sum_ranges; + uint64_t sum_vol; + Coordinate slice_dims; + Coordinate Slice; + + std::vector files; + int Nd; + GridBase* grid; + Grid::LatticeComplexD* grid_data; + + double rms; + + vtkImageData* imageData = nullptr; + vtkTextActor* text = nullptr; + isosurface* posExtractor = nullptr; + isosurface* negExtractor = nullptr; + + void SetGrid(GridBase* _grid) + { + grid = _grid; + Nd = grid->Nd(); + latt = grid->GlobalDimensions(); + grid_data = new Grid::LatticeComplexD(grid); + } + + void SetFiles(std::vector list) { files = list; old_file = -1; } + void SetSlice(Coordinate _Slice) { Slice = _Slice; } + + void SetSumDimensions(Coordinate _SumDims) + { + sum_dims = _SumDims; sum_ranges = Coordinate(Nd); sum_vol = 1; + for (int d = 0; d < Nd; d++) { sum_ranges[d] = sum_dims[d] ? latt[d] : 1; sum_vol *= sum_ranges[d]; } + } + + void SetLoopDimensions(Coordinate _LoopDims) + { + loop_dims = _LoopDims; loop_ranges = Coordinate(Nd); loop_vol = 1; + for (int d = 0; d < Nd; d++) { loop_ranges[d] = loop_dims[d] ? latt[d] : 1; loop_vol *= loop_ranges[d]; } + } + + void SetDisplayDimensions(Coordinate _xyz_dims) + { + xyz_dims = _xyz_dims; g_xyz_ranges = Coordinate(Nd); xyz_ranges = Coordinate(3); xyz_vol = 1; + for (int d = 0; d < 3; d++) { xyz_ranges[d] = latt[xyz_dims[d]]; xyz_vol *= xyz_ranges[d]; } + for (int d = 0; d < Nd; d++) { + g_xyz_ranges[d] = 1; + for (int dd = 0; dd < 3; dd++) if (xyz_dims[dd] == d) g_xyz_ranges[d] = latt[d]; + } + } + + void SetSliceDimensions() + { + Coordinate sd; + for (int d = 0; d < Nd; d++) { + if (g_xyz_ranges[d] > 1 || loop_dims[d] || sum_dims[d]) continue; + sd.push_back(d); + } + slice_dims = sd; + std::cout << " Slice dimensions: " << slice_dims << std::endl; + } + + void FillImageData(int loop_idx) + { + Coordinate loop_coor; + Lexicographic::CoorFromIndex(loop_coor, loop_idx, loop_ranges); + Coordinate xyz_coor(3), g_xyz_coor(Nd), sum_coor(Nd); + for (uint64_t xyz = 0; xyz < xyz_vol; xyz++) { + Lexicographic::CoorFromIndex(xyz_coor, xyz, xyz_ranges); + Lexicographic::CoorFromIndex(g_xyz_coor, xyz, g_xyz_ranges); + RealD val = 0.0; + for (uint64_t si = 0; si < sum_vol; si++) { + Lexicographic::CoorFromIndex(sum_coor, si, sum_ranges); + Coordinate site(Nd); + for (int d = 0; d < Nd; d++) + site[d] = (sum_coor[d] + loop_coor[d] + g_xyz_coor[d] + Slice[d]) % latt[d]; + val += real(peekSite(*grid_data, site)); + } + imageData->SetScalarComponentFromDouble(xyz_coor[0], xyz_coor[1], xyz_coor[2], 0, val); + } + imageData->Modified(); + } + + // Reload if needed, fill image, update label, render — no timer advance. + void ForceRender(vtkRenderWindowInteractor* iren) + { + int file = ((TimerCount / (int)loop_vol) + ffile) % (int)files.size(); + if (file != old_file) { + std::cout << "[render] Loading " << files[file] << std::endl; + readFile(*grid_data, files[file]); + old_file = file; + } + FillImageData(TimerCount % (int)loop_vol); + UpdateLabel(file, TimerCount % (int)loop_vol); + iren->GetRenderWindow()->Render(); + } + + virtual void Execute(vtkObject* caller, unsigned long eventId, void* callData) + { + if (vtkCommand::KeyPressEvent == eventId) { + vtkRenderWindowInteractor* iren = static_cast(caller); + std::string key = iren->GetKeySym(); + if (slice_dims.size() > 0) { + int vert = slice_dims[slice_dims.size() - 1]; + int horz = slice_dims[0]; + if (key == "Up") Slice[vert] = (Slice[vert] + 1) % latt[vert]; + if (key == "Down") Slice[vert] = (Slice[vert] + latt[vert] - 1) % latt[vert]; + if (key == "Right") Slice[horz] = (Slice[horz] + 1) % latt[horz]; + if (key == "Left") Slice[horz] = (Slice[horz] + latt[horz] - 1) % latt[horz]; + } + if (key == "greater") ffile = (ffile + 1) % (int)files.size(); + if (key == "less") ffile = (ffile - 1 + (int)files.size()) % (int)files.size(); + ForceRender(iren); + return; + } + + if (vtkCommand::TimerEvent == eventId) { + // timerId == -2: no animation timer (--notime), ignore all timer events + if (timerId < 0) return; + int tid = *(reinterpret_cast(callData)); + if (tid != timerId) return; + int file = ((TimerCount / (int)loop_vol) + ffile) % (int)files.size(); + if (file != old_file) { readFile(*grid_data, files[file]); old_file = file; } + FillImageData(TimerCount % (int)loop_vol); + UpdateLabel(file, TimerCount % (int)loop_vol); + dynamic_cast(caller)->GetRenderWindow()->Render(); + ++TimerCount; + if (TimerCount >= maxCount && timerId > -1) + dynamic_cast(caller)->DestroyTimer(timerId); + } + } + +private: + int TimerCount; + + void UpdateLabel(int file, int loop_idx) + { + Coordinate loop_coor; + Lexicographic::CoorFromIndex(loop_coor, loop_idx, loop_ranges); + // Extract tau value from filename (last '_'-delimited field) + const std::string& path = files[file]; + std::string tau = path.substr(path.rfind('_') + 1); + std::stringstream ss; + ss << "tau = " << tau << "\nSlice " << Slice; + text->SetInput(ss.str().c_str()); + } +}; + +// ─── Typewriter caption state ───────────────────────────────────────────────── +// User caption (gold, upper line) — cleared on new user: instruction +static std::string g_userCaptionFull; +static size_t g_userCaptionPos = 0; +// Claude caption (light blue, lower line) — cleared when user: arrives +static std::string g_claudeCaptionFull; +static size_t g_claudeCaptionPos = 0; +static int g_captionTick = 0; +static const int g_captionRate = 1; // ticks per character (1 x 100ms = 10 chars/sec) + +static std::string WrapText(const std::string& s, int maxCols = 45) { + std::istringstream words(s); + std::string word, line, result; + while (words >> word) { + if (!line.empty() && (int)(line.size() + 1 + word.size()) > maxCols) { + result += line + "\n"; + line = word; + } else { + if (!line.empty()) line += " "; + line += word; + } + } + if (!line.empty()) result += line; + return result; +} + +// ─── CommandHandler ─────────────────────────────────────────────────────────── +// Minimal parser for the wire protocol. Natural-language interpretation +// is handled externally (by Claude) before commands reach this program. +class CommandHandler : public vtkCallbackCommand +{ +public: + static CommandHandler* New() { return new CommandHandler; } + + FrameUpdater* fu; + vtkCamera* camera; + vtkRenderer* renderer; + vtkRenderWindowInteractor* iren; + vtkTextActor* captionActor = nullptr; // claude (light blue, lower) + vtkTextActor* userCaptionActor = nullptr; // user (gold, upper) + int pollTimerId = -1; + double isosurfaceLevel = 1.0; // in RMS units + + void CaptureFrame() { + if (g_recording && g_mpegWriter && g_imageFilter) { + GenerateAudioFrame(); + g_imageFilter->Modified(); + g_mpegWriter->Write(); + } + } + + void SetIsosurface(double level) + { + isosurfaceLevel = std::max(0.0, std::min(10.0, level)); + fu->posExtractor->SetValue(0, isosurfaceLevel * fu->rms); + fu->negExtractor->SetValue(0, -isosurfaceLevel * fu->rms); + fu->posExtractor->Modified(); + fu->negExtractor->Modified(); + } + + void PrintStatus() + { + std::cout << "[status] file = " << fu->ffile + << " : " << fu->files[fu->ffile] << "\n" + << "[status] Slice = " << fu->Slice << "\n" + << "[status] latt = " << fu->latt << "\n" + << "[status] isosurface = " << isosurfaceLevel + << " x RMS (" << isosurfaceLevel * fu->rms << ")" << std::endl; + } + + // Execute one line of the wire protocol. + void RunLine(const std::string& line) + { + std::istringstream iss(line); + std::string verb; + if (!(iss >> verb)) return; + + // ── slice ──────────────────────────────────────────── + if (verb == "slice") { + int dim; std::string valstr; + if (!(iss >> dim >> valstr)) { std::cout << "[cmd] slice: expected dim value" << std::endl; return; } + if (dim < 0 || dim >= fu->Nd) { std::cout << "[cmd] slice: dim out of range" << std::endl; return; } + int n = (int)fu->latt[dim]; + int newval; + if (!valstr.empty() && (valstr[0] == '+' || valstr[0] == '-')) { + int delta = std::stoi(valstr); + newval = ((fu->Slice[dim] + delta) % n + n) % n; + } else { + newval = ((std::stoi(valstr) % n) + n) % n; + } + fu->Slice[dim] = newval; + fu->ForceRender(iren); + PrintStatus(); + } + + // ── zoom ──────────────────────────────────────────────────── + else if (verb == "zoom") { + double factor; + if (!(iss >> factor)) { std::cout << "[cmd] zoom: expected factor" << std::endl; return; } + camera->Dolly(factor); + renderer->ResetCameraClippingRange(); + iren->GetRenderWindow()->Render(); + } + + // ── azimuth ──────────────────────────────────────────────── + else if (verb == "azimuth") { + double deg; + if (!(iss >> deg)) { std::cout << "[cmd] azimuth: expected degrees" << std::endl; return; } + camera->Azimuth(deg); + renderer->ResetCameraClippingRange(); + iren->GetRenderWindow()->Render(); + } + + // ── elevation ────────────────────────────────────────────── + else if (verb == "elevation") { + double deg; + if (!(iss >> deg)) { std::cout << "[cmd] elevation: expected degrees" << std::endl; return; } + camera->Elevation(deg); + renderer->ResetCameraClippingRange(); + iren->GetRenderWindow()->Render(); + } + + // ── spin ────────────────────────────────────────── + // Applies azimuth rotation on every 100ms poll tick. spin 0 stops. + else if (verb == "spin") { + double deg; + if (!(iss >> deg)) { std::cout << "[cmd] spin: expected degrees" << std::endl; return; } + g_spinDeg = deg; + std::cout << "[cmd] spin rate = " << g_spinDeg << " deg/tick" << std::endl; + } + + // ── caption user: / caption claude: / caption ───────── + // user: clears both lines, types user text (gold) on upper line. + // claude: keeps user line, types response (light blue) on lower line. + // caption alone clears both immediately. + else if (verb == "caption") { + std::string rest; + std::getline(iss, rest); + if (!rest.empty() && rest[0] == ' ') rest = rest.substr(1); + if (rest.empty()) { + g_userCaptionFull = ""; g_userCaptionPos = 0; + g_claudeCaptionFull = ""; g_claudeCaptionPos = 0; + g_captionTick = 0; + if (userCaptionActor) userCaptionActor->SetInput(""); + if (captionActor) captionActor->SetInput(""); + iren->GetRenderWindow()->Render(); CaptureFrame(); + } else if (rest.substr(0,5) == "user:") { + // New instruction: clear both, start typing user text + g_claudeCaptionFull = ""; g_claudeCaptionPos = 0; + g_userCaptionFull = WrapText(rest); g_userCaptionPos = 0; + g_captionTick = 0; + if (userCaptionActor) userCaptionActor->SetInput(""); + if (captionActor) captionActor->SetInput(""); + iren->GetRenderWindow()->Render(); CaptureFrame(); + } else { + // claude: or unlabelled — keep user line, type below + g_claudeCaptionFull = WrapText(rest); g_claudeCaptionPos = 0; + g_captionTick = 0; + } + } + + // ── record start / record stop ──────────────────────────── + else if (verb == "record") { +#ifdef MPEG + std::string sub; + if (!(iss >> sub)) { std::cout << "[cmd] record: expected start or stop" << std::endl; return; } + if (sub == "stop") { + if (g_recording && g_mpegWriter) { + g_mpegWriter->End(); + g_mpegWriter->Delete(); g_mpegWriter = nullptr; + g_imageFilter->Delete(); g_imageFilter = nullptr; + g_recording = false; + std::cout << "[cmd] recording stopped: " << g_recordingFile << std::endl; + // Write WAV and mux to MP4 + std::string wavFile = g_recordingFile + ".wav"; + WriteWAV(wavFile, g_audioBuffer, AUDIO_RATE); + g_audioBuffer.clear(); + std::string mp4File = g_recordingFile; + if (mp4File.size() > 4 && mp4File.substr(mp4File.size()-4) == ".avi") + mp4File = mp4File.substr(0, mp4File.size()-4) + ".mp4"; + else + mp4File += ".mp4"; + std::string cmd = "ffmpeg -y -i \"" + g_recordingFile + "\" -i \"" + wavFile + + "\" -c:v copy -c:a aac -shortest \"" + mp4File + "\" 2>/dev/null"; + int ret = system(cmd.c_str()); + if (ret == 0) { + std::cout << "[cmd] muxed output: " << mp4File << std::endl; + unlink(wavFile.c_str()); // clean up intermediate WAV + } else { + std::cout << "[cmd] mux failed (ffmpeg not found?). WAV kept: " << wavFile << std::endl; + } + } else { + std::cout << "[cmd] not recording" << std::endl; + } + } else if (sub == "start") { + std::string fname = "recording.avi"; + iss >> fname; + if (g_recording) { std::cout << "[cmd] already recording" << std::endl; return; } + g_recordingFile = fname; + g_audioBuffer.clear(); + g_samplesPerFrame = AUDIO_RATE / std::max(1, g_framerate); + g_beepRemaining = 0; + g_beepPhase = 0.0; + g_imageFilter = vtkWindowToImageFilter::New(); + g_imageFilter->SetInput(iren->GetRenderWindow()); + g_imageFilter->SetInputBufferTypeToRGB(); + g_mpegWriter = vtkFFMPEGWriter::New(); + g_mpegWriter->SetFileName(fname.c_str()); + g_mpegWriter->SetRate(g_framerate); + g_mpegWriter->SetInputConnection(g_imageFilter->GetOutputPort()); + g_mpegWriter->Start(); + g_recording = true; + std::cout << "[cmd] recording started: " << fname << std::endl; + } else { + std::cout << "[cmd] record: unknown subcommand '" << sub << "'" << std::endl; + } +#else + std::cout << "[cmd] record: MPEG support not compiled" << std::endl; +#endif + } + + // ── iso ────────────────────────────────────────────────────── + else if (verb == "iso") { + double val; + if (!(iss >> val)) { std::cout << "[cmd] iso: expected value" << std::endl; return; } + SetIsosurface(val); + fu->ForceRender(iren); + PrintStatus(); + } + + // ── file ─────────────────────────────────────────────── + else if (verb == "file") { + std::string valstr; + if (!(iss >> valstr)) { std::cout << "[cmd] file: expected index" << std::endl; return; } + int n = (int)fu->files.size(); + int newval; + if (!valstr.empty() && (valstr[0] == '+' || valstr[0] == '-')) { + int delta = std::stoi(valstr); + newval = ((fu->ffile + delta) % n + n) % n; + } else { + newval = ((std::stoi(valstr) % n) + n) % n; + } + fu->ffile = newval; + fu->old_file = -1; + fu->ForceRender(iren); + PrintStatus(); + } + + // ── render ─────────────────────────────────────────────────────────── + else if (verb == "render") { + fu->ForceRender(iren); + } + + // ── status ─────────────────────────────────────────────────────────── + else if (verb == "status") { + PrintStatus(); + } + + // ── quit ───────────────────────────────────────────────────────────── + else if (verb == "quit" || verb == "exit") { + g_running = false; + iren->TerminateApp(); + } + + else { + std::cout << "[cmd] Unknown command: '" << line << "'" << std::endl; + } + } + + virtual void Execute(vtkObject*, unsigned long eventId, void* callData) + { + if (eventId != vtkCommand::TimerEvent) return; + if (pollTimerId >= 0) { + int tid = *(reinterpret_cast(callData)); + if (tid != pollTimerId) return; + } + + std::vector pending; + { + std::lock_guard lk(g_cmdMutex); + while (!g_cmdQueue.empty()) { pending.push_back(g_cmdQueue.front()); g_cmdQueue.pop(); } + } + for (const auto& line : pending) { + std::cout << "[cmd] >> " << line << std::endl; + RunLine(line); + // CaptureFrame() called inside RunLine for caption; for other + // rendering commands capture here (duplicate Modified() is harmless) + CaptureFrame(); + } + + // Typewriter: advance one character every g_captionRate ticks. + // User line types first; claude line starts once user line is complete. + bool typing = (g_userCaptionPos < g_userCaptionFull.size()) || + (g_claudeCaptionPos < g_claudeCaptionFull.size()); + if (typing) { + if (++g_captionTick >= g_captionRate) { + g_captionTick = 0; + bool rendered = false; + if (g_userCaptionPos < g_userCaptionFull.size()) { + ++g_userCaptionPos; + if (userCaptionActor) + userCaptionActor->SetInput(g_userCaptionFull.substr(0, g_userCaptionPos).c_str()); + PlayBeepAudible(); + TriggerBeep(); + rendered = true; + } else if (g_claudeCaptionPos < g_claudeCaptionFull.size()) { + ++g_claudeCaptionPos; + if (captionActor) + captionActor->SetInput(g_claudeCaptionFull.substr(0, g_claudeCaptionPos).c_str()); + rendered = true; + } + if (rendered) { + iren->GetRenderWindow()->Render(); + CaptureFrame(); + } + } + } + + // Apply continuous spin (if active) at poll-timer rate + if (g_spinDeg != 0.0) { + camera->Azimuth(g_spinDeg); + renderer->ResetCameraClippingRange(); + iren->GetRenderWindow()->Render(); + CaptureFrame(); + } + } +}; + +// ─── main ───────────────────────────────────────────────────────────────────── +int main(int argc, char* argv[]) +{ + using namespace Grid; + + Grid_init(&argc, &argv); + GridLogLayout(); + + auto latt_size = GridDefaultLatt(); + auto simd_layout = GridDefaultSimd(latt_size.size(), vComplex::Nsimd()); + auto mpi_layout = GridDefaultMpi(); + GridCartesian Grid(latt_size, simd_layout, mpi_layout); + + double default_contour = 1.0; + std::string arg; + + std::vector file_list({"file1","file2","file3","file4", + "file5","file6","file7","file8"}); + + if (GridCmdOptionExists(argv, argv+argc, "--files")) { + arg = GridCmdOptionPayload(argv, argv+argc, "--files"); + GridCmdOptionCSL(arg, file_list); + } +#ifdef MPEG + if (GridCmdOptionExists(argv, argv+argc, "--mpeg")) g_mpeg = 1; +#endif + if (GridCmdOptionExists(argv, argv+argc, "--fps")) { + arg = GridCmdOptionPayload(argv, argv+argc, "--fps"); + GridCmdOptionInt(arg, g_framerate); + } + if (GridCmdOptionExists(argv, argv+argc, "--isosurface")) { + arg = GridCmdOptionPayload(argv, argv+argc, "--isosurface"); + GridCmdOptionFloat(arg, default_contour); + } + + int NoTime = 0, Nd = Grid.Nd(); + Coordinate Slice(Nd,0), SumDims(Nd,0), LoopDims(Nd,0), XYZDims({0,1,2}); + + if (GridCmdOptionExists(argv, argv+argc, "--slice")) { + arg = GridCmdOptionPayload(argv, argv+argc, "--slice"); + GridCmdOptionIntVector(arg, Slice); + } + if (GridCmdOptionExists(argv, argv+argc, "--sum")) { + arg = GridCmdOptionPayload(argv, argv+argc, "--sum"); + GridCmdOptionIntVector(arg, SumDims); + } + if (GridCmdOptionExists(argv, argv+argc, "--loop")) { + arg = GridCmdOptionPayload(argv, argv+argc, "--loop"); + GridCmdOptionIntVector(arg, LoopDims); + } + if (GridCmdOptionExists(argv, argv+argc, "--xyz")) { + arg = GridCmdOptionPayload(argv, argv+argc, "--xyz"); + GridCmdOptionIntVector(arg, XYZDims); + } + if (GridCmdOptionExists(argv, argv+argc, "--notime")) { NoTime = 1; } + + std::thread cmdThread(CommandReaderThread); + cmdThread.detach(); + + // ── VTK scene ──────────────────────────────────────────────────────────── + vtkNew colors; + std::array posColor{{240,184,160,255}}; colors->SetColor("posColor", posColor.data()); + std::array bkg{{51,77,102,255}}; colors->SetColor("BkgColor", bkg.data()); + + vtkNew renWin; + vtkNew iren; + iren->SetRenderWindow(renWin); + + int frameCount = (int)file_list.size(); + for (int d = 0; d < Nd; d++) if (LoopDims[d]) frameCount *= latt_size[d]; + + vtkNew aCamera; + aCamera->SetViewUp(0,0,-1); aCamera->SetPosition(0,-1000,0); aCamera->SetFocalPoint(0,0,0); + aCamera->ComputeViewPlaneNormal(); aCamera->Azimuth(30.0); aCamera->Elevation(30.0); + + vtkNew aRenderer; + renWin->AddRenderer(aRenderer); + + double nrm, rms, contour; + { LatticeComplexD data(&Grid); readFile(data, file_list[0]); nrm = norm2(data); } + rms = std::sqrt(nrm / Grid.gSites()); + contour = default_contour * rms; + + vtkNew imageData; + imageData->SetDimensions(latt_size[XYZDims[0]], latt_size[XYZDims[1]], latt_size[XYZDims[2]]); + imageData->AllocateScalars(VTK_DOUBLE, 1); + for (int xx=0;xxSetScalarComponentFromDouble(xx,yy,zz,0,0.0); + + vtkNew posExtractor; posExtractor->SetInputData(imageData); posExtractor->SetValue(0, contour); + vtkNew posStripper; posStripper->SetInputConnection(posExtractor->GetOutputPort()); + vtkNew posMapper; posMapper->SetInputConnection(posStripper->GetOutputPort()); posMapper->ScalarVisibilityOff(); + vtkNew pos; pos->SetMapper(posMapper); + pos->GetProperty()->SetDiffuseColor(colors->GetColor3d("posColor").GetData()); + pos->GetProperty()->SetSpecular(0.3); pos->GetProperty()->SetSpecularPower(20); pos->GetProperty()->SetOpacity(0.5); + + vtkNew negExtractor; negExtractor->SetInputData(imageData); negExtractor->SetValue(0, -contour); + vtkNew negStripper; negStripper->SetInputConnection(negExtractor->GetOutputPort()); + vtkNew negMapper; negMapper->SetInputConnection(negStripper->GetOutputPort()); negMapper->ScalarVisibilityOff(); + vtkNew neg; neg->SetMapper(negMapper); + neg->GetProperty()->SetDiffuseColor(colors->GetColor3d("Ivory").GetData()); + + vtkNew outlineData; outlineData->SetInputData(imageData); + vtkNew mapOutline; mapOutline->SetInputConnection(outlineData->GetOutputPort()); + vtkNew outline; outline->SetMapper(mapOutline); + outline->GetProperty()->SetColor(colors->GetColor3d("Black").GetData()); + + vtkNew TextT; + TextT->SetInput("Initialising..."); + TextT->SetPosition(10, 920); + TextT->GetTextProperty()->SetFontSize(24); + TextT->GetTextProperty()->SetColor(colors->GetColor3d("Gold").GetData()); + + // Claude response caption (light blue, lower line) + vtkNew CaptionT; + CaptionT->SetInput(""); + CaptionT->SetPosition(512, 38); + CaptionT->GetTextProperty()->SetFontSize(32); + CaptionT->GetTextProperty()->SetColor(0.6, 0.9, 1.0); + CaptionT->GetTextProperty()->SetJustificationToCentered(); + CaptionT->GetTextProperty()->SetBackgroundColor(0.0, 0.0, 0.0); + CaptionT->GetTextProperty()->SetBackgroundOpacity(0.6); + CaptionT->GetTextProperty()->BoldOn(); + + // User instruction caption (gold, upper line) + vtkNew UserCaptionT; + UserCaptionT->SetInput(""); + UserCaptionT->SetPosition(512, 82); + UserCaptionT->GetTextProperty()->SetFontSize(32); + UserCaptionT->GetTextProperty()->SetColor(1.0, 0.85, 0.0); + UserCaptionT->GetTextProperty()->SetJustificationToCentered(); + UserCaptionT->GetTextProperty()->SetBackgroundColor(0.0, 0.0, 0.0); + UserCaptionT->GetTextProperty()->SetBackgroundOpacity(0.6); + UserCaptionT->GetTextProperty()->BoldOn(); + + aRenderer->AddActor(TextT); aRenderer->AddActor(CaptionT); aRenderer->AddActor(UserCaptionT); aRenderer->AddActor(outline); + aRenderer->AddActor(pos); aRenderer->AddActor(neg); + + vtkNew fu; + fu->SetGrid(&Grid); fu->SetFiles(file_list); fu->SetSlice(Slice); + fu->SetSumDimensions(SumDims); fu->SetLoopDimensions(LoopDims); + fu->SetDisplayDimensions(XYZDims); fu->SetSliceDimensions(); + fu->imageData = imageData; fu->text = TextT; fu->maxCount = frameCount; + fu->posExtractor = posExtractor; fu->negExtractor = negExtractor; fu->rms = rms; + + iren->AddObserver(vtkCommand::TimerEvent, fu); + iren->AddObserver(vtkCommand::KeyPressEvent, fu); + + aRenderer->SetActiveCamera(aCamera); aRenderer->ResetCamera(); + aRenderer->SetBackground(colors->GetColor3d("BkgColor").GetData()); + aCamera->Dolly(1.0); aRenderer->SetViewport(0.0,0.0,1.0,1.0); + aRenderer->ResetCameraClippingRange(); + renWin->SetSize(1024,1024); renWin->SetWindowName("ControlledFieldDensity"); + renWin->Render(); iren->Initialize(); + + // CommandHandler on fast poll timer + vtkNew cmdHandler; + cmdHandler->fu = fu; + cmdHandler->camera = aCamera; + cmdHandler->renderer = aRenderer; + cmdHandler->iren = iren; + cmdHandler->captionActor = CaptionT; + cmdHandler->userCaptionActor = UserCaptionT; + cmdHandler->isosurfaceLevel = default_contour; + iren->AddObserver(vtkCommand::TimerEvent, cmdHandler); + cmdHandler->pollTimerId = iren->CreateRepeatingTimer(100); + + if (g_mpeg == 0 && NoTime == 0) { + fu->timerId = iren->CreateRepeatingTimer(10000 / g_framerate); + } + + if (g_mpeg) { +#ifdef MPEG + vtkWindowToImageFilter* imageFilter = vtkWindowToImageFilter::New(); + imageFilter->SetInput(renWin); imageFilter->SetInputBufferTypeToRGB(); + vtkFFMPEGWriter* writer = vtkFFMPEGWriter::New(); + writer->SetFileName("movie.avi"); writer->SetRate(g_framerate); + writer->SetInputConnection(imageFilter->GetOutputPort()); writer->Start(); + for (int i = 0; i < fu->maxCount; i++) { + fu->Execute(iren, vtkCommand::TimerEvent, &fu->timerId); + imageFilter->Modified(); writer->Write(); + } + writer->End(); writer->Delete(); imageFilter->Delete(); +#else + assert(-1 && "MPEG support not compiled"); +#endif + } else { + iren->Start(); + } + + g_running = false; + Grid_finalize(); + return EXIT_SUCCESS; +} diff --git a/visualisation/ForceAnalysis.cxx b/visualisation/ForceAnalysis.cxx new file mode 100644 index 00000000..e1111be4 --- /dev/null +++ b/visualisation/ForceAnalysis.cxx @@ -0,0 +1,633 @@ +// ForceAnalysis.cxx +// +// Reads a sequence of force snapshot files (LatticeComplexD, real part = force magnitude) +// and produces two outputs: +// +// 1. Tab-separated timeseries to stdout: +// idx Gauge_lat_rms Gauge_lat_hot Gauge_smr_rms ... +// where _rms is the lattice RMS and _hot is the value at --hotsite. +// +// 2. PNG images (one per force component per snapshot) rendered via VTK +// as isosurfaces of the force density, using the same pipeline as +// Visualise5D. Images are written to --pngdir/