1
0
mirror of https://github.com/paboyle/Grid.git synced 2026-04-23 03:56:13 +01:00

Compare commits

..

1 Commits

Author SHA1 Message Date
portelli 6e40e22004 aarch64 libunwind compatibility fix 2025-10-30 16:45:43 +01:00
81 changed files with 151 additions and 6488 deletions
+7 -9
View File
@@ -54,24 +54,22 @@ Version.h: version-cache
include Make.inc
include Eigen.inc
if BUILD_FERMION_INSTANTIATIONS
extra_sources+=$(WILS_FERMION_FILES)
extra_sources+=$(STAG_FERMION_FILES)
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)
endif
extra_sources+=$(SP_TWOIND_FERMION_FILES)
endif
endif
+16 -79
View File
@@ -28,7 +28,6 @@ Author: Peter Boyle <pboyle@bnl.gov>
#pragma once
#ifdef GRID_HIP
#include <hip/hip_version.h>
#include <hipblas/hipblas.h>
#endif
#ifdef GRID_CUDA
@@ -256,29 +255,16 @@ 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,
(hipDoubleComplex *) &alpha_p[0],
(hipDoubleComplex **)&Amk[0], lda,
(hipDoubleComplex **)&Bkn[0], ldb,
(hipDoubleComplex *) &beta_p[0],
(hipDoubleComplex **)&Cmn[0], ldc,
(hipblasDoubleComplex *) &alpha_p[0],
(hipblasDoubleComplex **)&Amk[0], lda,
(hipblasDoubleComplex **)&Bkn[0], ldb,
(hipblasDoubleComplex *) &beta_p[0],
(hipblasDoubleComplex **)&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<<std::endl;
GRID_ASSERT(err==HIPBLAS_STATUS_SUCCESS);
#endif
@@ -517,30 +503,17 @@ 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 = hipblasCgemmBatched(gridblasHandle,
hOpA,
hOpB,
m,n,k,
(hipComplex *) &alpha_p[0],
(hipComplex **)&Amk[0], lda,
(hipComplex **)&Bkn[0], ldb,
(hipComplex *) &beta_p[0],
(hipComplex **)&Cmn[0], ldc,
(hipblasComplex *) &alpha_p[0],
(hipblasComplex **)&Amk[0], lda,
(hipblasComplex **)&Bkn[0], ldb,
(hipblasComplex *) &beta_p[0],
(hipblasComplex **)&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
@@ -1121,19 +1094,11 @@ public:
GRID_ASSERT(info.size()==batchCount);
#ifdef GRID_HIP
#if defined(HIP_VERSION_MAJOR) && (HIP_VERSION_MAJOR >=7)
auto err = hipblasZgetrfBatched(gridblasHandle,(int)n,
(hipDoubleComplex **)&Ann[0], (int)n,
(hipblasDoubleComplex **)&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
@@ -1159,20 +1124,11 @@ public:
GRID_ASSERT(info.size()==batchCount);
#ifdef GRID_HIP
#if defined(HIP_VERSION_MAJOR) && (HIP_VERSION_MAJOR >=7)
auto err = hipblasCgetrfBatched(gridblasHandle,(int)n,
(hipComplex **)&Ann[0], (int)n,
(hipblasComplex **)&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
@@ -1245,22 +1201,12 @@ public:
GRID_ASSERT(Cnn.size()==batchCount);
#ifdef GRID_HIP
#if defined(HIP_VERSION_MAJOR) && (HIP_VERSION_MAJOR >=7)
auto err = hipblasZgetriBatched(gridblasHandle,(int)n,
(hipDoubleComplex **)&Ann[0], (int)n,
(hipblasDoubleComplex **)&Ann[0], (int)n,
(int*) &ipiv[0],
(hipDoubleComplex **)&Cnn[0], (int)n,
(hipblasDoubleComplex **)&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
@@ -1289,21 +1235,12 @@ public:
GRID_ASSERT(Cnn.size()==batchCount);
#ifdef GRID_HIP
#if defined(HIP_VERSION_MAJOR) && (HIP_VERSION_MAJOR >=7)
auto err = hipblasCgetriBatched(gridblasHandle,(int)n,
(hipComplex **)&Ann[0], (int)n,
(hipblasComplex **)&Ann[0], (int)n,
(int*) &ipiv[0],
(hipComplex **)&Cnn[0], (int)n,
(hipblasComplex **)&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
+2 -2
View File
@@ -92,8 +92,8 @@ class TwoLevelCGmrhs
// Vector case
virtual void operator() (std::vector<Field> &src, std::vector<Field> &x)
{
SolveSingleSystem(src,x);
// SolvePrecBlockCG(src,x);
// SolveSingleSystem(src,x);
SolvePrecBlockCG(src,x);
}
////////////////////////////////////////////////////////////////////////////////////////////////////
+6 -9
View File
@@ -97,7 +97,7 @@ public:
RealD scale;
ConjugateGradient<FineField> CG(1.0e-4,2000,false);
ConjugateGradient<FineField> CG(1.0e-3,400,false);
FineField noise(FineGrid);
FineField Mn(FineGrid);
@@ -131,10 +131,7 @@ public:
RealD scale;
TrivialPrecon<FineField> simple_fine;
// PrecGeneralisedConjugateResidualNonHermitian<FineField> GCR(0.001,10,DiracOp,simple_fine,30,30);
// PrecGeneralisedConjugateResidualNonHermitian<FineField> GCR(0.001,10,DiracOp,simple_fine,12,12);
// PrecGeneralisedConjugateResidualNonHermitian<FineField> GCR(0.001,30,DiracOp,simple_fine,12,12);
PrecGeneralisedConjugateResidualNonHermitian<FineField> GCR(0.001,30,DiracOp,simple_fine,10,10);
PrecGeneralisedConjugateResidualNonHermitian<FineField> GCR(0.001,30,DiracOp,simple_fine,12,12);
FineField noise(FineGrid);
FineField src(FineGrid);
FineField guess(FineGrid);
@@ -149,16 +146,16 @@ public:
DiracOp.Op(noise,Mn); std::cout<<GridLogMessage << "noise ["<<b<<"] <n|Op|n> "<<innerProduct(noise,Mn)<<std::endl;
for(int i=0;i<3;i++){
for(int i=0;i<2;i++){
// void operator() (const Field &src, Field &psi){
#if 1
if (i==0)std::cout << GridLogMessage << " inverting on noise "<<std::endl;
std::cout << GridLogMessage << " inverting on noise "<<std::endl;
src = noise;
guess=Zero();
GCR(src,guess);
subspace[b] = guess;
#else
if (i==0)std::cout << GridLogMessage << " inverting on zero "<<std::endl;
std::cout << GridLogMessage << " inverting on zero "<<std::endl;
src=Zero();
guess = noise;
GCR(src,guess);
@@ -170,7 +167,7 @@ public:
}
DiracOp.Op(noise,Mn); std::cout<<GridLogMessage << "filtered["<<b<<"] <f|Op|f> "<<innerProduct(noise,Mn)<<" <f|OpDagOp|f>"<<norm2(Mn)<<std::endl;
DiracOp.Op(noise,Mn); std::cout<<GridLogMessage << "filtered["<<b<<"] <f|Op|f> "<<innerProduct(noise,Mn)<<std::endl;
subspace[b] = noise;
}
+1 -2
View File
@@ -27,8 +27,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
/* END LEGAL */
#include <Grid/GridCore.h>
void GridAbort(void) { abort(); }
NAMESPACE_BEGIN(Grid);
///////////////////////////////////////////////////////////////////////////////////////////////////
@@ -36,6 +34,7 @@ NAMESPACE_BEGIN(Grid);
///////////////////////////////////////////////////////////////////////////////////////////////////
Grid_MPI_Comm CartesianCommunicator::communicator_world;
void GridAbort(void) { abort(); }
void CartesianCommunicator::Init(int *argc, char *** arv)
{
+2 -2
View File
@@ -124,7 +124,7 @@ template<class vobj> void Cshift_simple(Lattice<vobj>& ret,const Lattice<vobj> &
void *hsend_buf = (void *)&hrhs[0];
void *hrecv_buf = (void *)&hret[0];
acceleratorCopyFromDevice(send_buf,hsend_buf,bytes);
acceleratorCopyFromDevice(&send_buf[0],&hsend_buf[0],bytes);
grid->SendToRecvFrom(hsend_buf,
xmit_to_rank,
@@ -132,7 +132,7 @@ template<class vobj> void Cshift_simple(Lattice<vobj>& ret,const Lattice<vobj> &
recv_from_rank,
bytes);
acceleratorCopyToDevice(hrecv_buf,recv_buf,bytes);
acceleratorCopyToDevice(&hrecv_buf[0],&recv_buf[0],bytes);
#endif
}
}
-18
View File
@@ -596,32 +596,16 @@ template<int Index,class vobj> inline vobj transposeColour(const vobj &lhs){
//////////////////////////////////////////
// Trace lattice and non-lattice
//////////////////////////////////////////
#define GRID_UNOP(name) name
#define GRID_DEF_UNOP(op, name) \
template <typename T1, typename std::enable_if<is_lattice<T1>::value||is_lattice_expr<T1>::value,T1>::type * = nullptr> \
inline auto op(const T1 &arg) ->decltype(LatticeUnaryExpression<GRID_UNOP(name),T1>(GRID_UNOP(name)(), arg)) \
{ \
return LatticeUnaryExpression<GRID_UNOP(name),T1>(GRID_UNOP(name)(), arg); \
}
template<int Index,class vobj>
inline auto traceSpin(const Lattice<vobj> &lhs) -> Lattice<decltype(traceIndex<SpinIndex>(vobj()))>
{
return traceIndex<SpinIndex>(lhs);
}
GridUnopClass(UnaryTraceSpin, traceIndex<SpinIndex>(a));
GRID_DEF_UNOP(traceSpin, UnaryTraceSpin);
template<int Index,class vobj>
inline auto traceColour(const Lattice<vobj> &lhs) -> Lattice<decltype(traceIndex<ColourIndex>(vobj()))>
{
return traceIndex<ColourIndex>(lhs);
}
GridUnopClass(UnaryTraceColour, traceIndex<ColourIndex>(a));
GRID_DEF_UNOP(traceColour, UnaryTraceColour);
template<int Index,class vobj>
inline auto traceSpin(const vobj &lhs) -> Lattice<decltype(traceIndex<SpinIndex>(lhs))>
{
@@ -633,8 +617,6 @@ inline auto traceColour(const vobj &lhs) -> Lattice<decltype(traceIndex<ColourIn
return traceIndex<ColourIndex>(lhs);
}
#undef GRID_UNOP
#undef GRID_DEF_UNOP
//////////////////////////////////////////
// Current types
//////////////////////////////////////////
@@ -411,7 +411,7 @@ void WilsonKernels<Impl>::DhopDirKernel( StencilImpl &st, DoubledGaugeField &U,S
#undef LoopBody
}
#if 0
#ifdef GRID_SYCL
extern "C" {
ulong SYCL_EXTERNAL __attribute__((overloadable)) intel_get_cycle_counter( void );
uint SYCL_EXTERNAL __attribute__((overloadable)) intel_get_active_channel_mask( void );
+3 -6
View File
@@ -138,13 +138,10 @@ public:
//auto start = std::chrono::high_resolution_clock::now();
autoView(U_v,U,AcceleratorWrite);
autoView(P_v,P,AcceleratorRead);
typedef typename Field::vector_object vobj;
const int Nsimd = vobj::Nsimd();
accelerator_for(ss, P.Grid()->oSites(),Nsimd,{
accelerator_for(ss, P.Grid()->oSites(),1,{
for (int mu = 0; mu < Nd; mu++) {
auto tmp = Exponentiate(P_v(ss)(mu), ep, Nexp) * U_v(ss)(mu);
tmp = Group::ProjectOnGeneralGroup(tmp);
coalescedWrite(U_v[ss](mu),tmp);
U_v[ss](mu) = Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu);
U_v[ss](mu) = Group::ProjectOnGeneralGroup(U_v[ss](mu));
}
});
//auto end = std::chrono::high_resolution_clock::now();
-12
View File
@@ -103,18 +103,6 @@ class PolyakovMod: public ObservableModule<PolyakovLogger<Impl>, NoParameters>{
PolyakovMod(): ObsBase(NoParameters()){}
};
template < class Impl >
class SpatialPolyakovMod: public ObservableModule<SpatialPolyakovLogger<Impl>, NoParameters>{
typedef ObservableModule<SpatialPolyakovLogger<Impl>, NoParameters> ObsBase;
using ObsBase::ObsBase; // for constructors
// acquire resource
virtual void initialize(){
this->ObservablePtr.reset(new SpatialPolyakovLogger<Impl>());
}
public:
SpatialPolyakovMod(): ObsBase(NoParameters()){}
};
template < class Impl >
class TopologicalChargeMod: public ObservableModule<TopologicalCharge<Impl>, TopologyObsParameters>{
+2 -42
View File
@@ -2,12 +2,11 @@
Grid physics library, www.github.com/paboyle/Grid
Source file: ./Grid/qcd/observables/polyakov_loop.h
Source file: ./lib/qcd/modules/polyakov_line.h
Copyright (C) 2025
Copyright (C) 2017
Author: David Preti <david.preti@csic.es>
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
@@ -61,43 +60,4 @@ class PolyakovLogger : public HmcObservable<typename Impl::Field> {
}
};
template <class Impl>
class SpatialPolyakovLogger : public HmcObservable<typename Impl::Field> {
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<Ndim; idx++) {
polyakov = WilsonLoops<Impl>::avgPolyakovLoop(U, idx);
std::cout << GridLogMessage
<< std::setprecision(std::numeric_limits<Real>::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);
+4 -7
View File
@@ -291,8 +291,8 @@ public:
int idx=0;
for(int mu=0;mu<4;mu++){
for(int nu=0;nu<4;nu++){
if ( mu!=nu) assert(this->StoutSmearing->SmearRho[idx]==rho);
else assert(this->StoutSmearing->SmearRho[idx]==0.0);
if ( mu!=nu) GRID_ASSERT(this->StoutSmearing->SmearRho[idx]==rho);
else GRID_ASSERT(this->StoutSmearing->SmearRho[idx]==0.0);
idx++;
}}
//////////////////////////////////////////////////////////////////
@@ -825,7 +825,6 @@ public:
virtual void fill_smearedSet(GaugeField &U)
{
this->ThinLinks = &U; // attach the smearing routine to the field U
std::cout << GridLogMessage << " fill_smearedSet " << WilsonLoops<PeriodicGimplR>::avgPlaquette(U) << std::endl;
// check the pointer is not null
if (this->ThinLinks == NULL)
@@ -847,8 +846,6 @@ public:
ApplyMask(smeared_A,smearLvl);
smeared_B = previous_u;
ApplyMask(smeared_B,smearLvl);
std::cout << GridLogMessage << " smeared_A " << norm2(smeared_A) << std::endl;
std::cout << GridLogMessage << " smeared_B " << norm2(smeared_B) << std::endl;
// Replace only the masked portion
this->SmearedSet[smearLvl] = previous_u-smeared_B + smeared_A;
previous_u = this->SmearedSet[smearLvl];
@@ -937,10 +934,10 @@ public:
SmearedConfigurationMasked(GridCartesian* _UGrid, unsigned int Nsmear, Smear_Stout<Gimpl>& Stout)
: SmearedConfiguration<Gimpl>(_UGrid, Nsmear,Stout)
{
assert(Nsmear%(2*Nd)==0); // Or multiply by 8??
GRID_ASSERT(Nsmear%(2*Nd)==0); // Or multiply by 8??
// was resized in base class
assert(this->SmearedSet.size()==Nsmear);
GRID_ASSERT(this->SmearedSet.size()==Nsmear);
GridRedBlackCartesian * UrbGrid;
UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(_UGrid);
+1 -1
View File
@@ -54,7 +54,7 @@ public:
// Usual cases are not used
//////////////////////////////////
virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG &pRNG){ GRID_ASSERT(0);};
virtual RealD S(const GaugeField &U) { GRID_ASSERT(0); return 0; }
virtual RealD S(const GaugeField &U) { GRID_ASSERT(0); }
virtual void deriv(const GaugeField &U, GaugeField &dSdU) { GRID_ASSERT(0); }
//////////////////////////////////
+3 -3
View File
@@ -254,9 +254,9 @@ static void testGenerators(GroupName::Sp) {
}
}
template <class vtype, int N>
static Lattice<iScalar<iScalar<iMatrix<vtype, N> > > >
ProjectOnGeneralGroup(const Lattice<iScalar<iScalar<iMatrix<vtype, N> > > > &Umu, GroupName::Sp) {
template <int N>
static Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > >
ProjectOnGeneralGroup(const Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu, GroupName::Sp) {
return ProjectOnSpGroup(Umu);
}
+8 -26
View File
@@ -177,43 +177,25 @@ public:
}
//////////////////////////////////////////////////
// average Polyakov loop in mu direction over all directions != mu
// average over all x,y,z the temporal loop
//////////////////////////////////////////////////
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());
static ComplexD avgPolyakovLoop(const GaugeField &Umu) { //assume Nd=4
GaugeMat Ut(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];
// Number of sites in mu direction
int N_mu = Umu.Grid()->GlobalDimensions()[mu];
U_loop = peekLorentz(Umu, mu); //Select direction
P = U_loop;
for (int t=1;t<N_mu;t++){
P = Gimpl::CovShiftForward(U_loop,mu,P);
Ut = peekLorentz(Umu,3); //Select temporal direction
P = Ut;
for (int t=1;t<T;t++){
P = Gimpl::CovShiftForward(Ut,3,P);
}
RealD norm = 1.0/(Nc*X*Y*Z*T);
out = sum(trace(P))*norm;
return out;
}
/////////////////////////////////////////////////
// overload for temporal Polyakov loop
/////////////////////////////////////////////////
static ComplexD avgPolyakovLoop(const GaugeField &Umu) {
return avgPolyakovLoop(Umu, 3);
}
}
//////////////////////////////////////////////////
// average over traced single links
+4 -6
View File
@@ -751,7 +751,7 @@ public:
obj.xbytes = xbytes;
obj.rbytes = rbytes;
obj.cb = cb;
for(int i=0;i<CachedTransfers.size();i++){
if ( (CachedTransfers[i].direction ==direction)
&&(CachedTransfers[i].OrthogPlane==OrthogPlane)
@@ -763,13 +763,11 @@ public:
){
// FIXME worry about duplicate with partial compression
// Wont happen as DWF has no duplicates, but...
// AddCopy(CachedTransfers[i].recv_buf,recv_buf,rbytes);
// std::cout << "Duplicate dir " <<direction<<" "<<" OrthogPlane "<<OrthogPlane<<" Dest"<<DestProc <<" xbytes " <<xbytes<<" lane "<< lane<<" cb "<<cb<<std::endl;
return 0;
// return 1;
AddCopy(CachedTransfers[i].recv_buf,recv_buf,rbytes);
return 1;
}
}
CachedTransfers.push_back(obj);
return 0;
}
+1 -1
View File
@@ -755,7 +755,7 @@ void Grid_generic_handler(int sig,siginfo_t *si,void * ptr)
sig_print_uint(si->si_code);
SIGLOG("\n");
ucontext_t *uc= (ucontext_t *)ptr;
unw_context_t *uc= (unw_context_t *)ptr;
SIGLOG("Backtrace:\n");
#ifdef HAVE_UNWIND
+1 -6
View File
@@ -24,11 +24,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include<Grid/Grid.h>
#include <Grid/Grid.h>
#if Nc == 3
#include <Grid/qcd/smearing/GaugeConfigurationMasked.h>
@@ -234,4 +230,3 @@ int main(int argc, char **argv)
#endif
} // main
#endif
+3 -6
View File
@@ -25,11 +25,7 @@ directory
*************************************************************************************/
/* END LEGAL */
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include<Grid/Grid.h>
#include <Grid/Grid.h>
#if Nc == 3
#include <Grid/qcd/smearing/GaugeConfigurationMasked.h>
@@ -235,4 +231,5 @@ int main(int argc, char **argv)
#endif
} // main
#endif
+3 -6
View File
@@ -24,11 +24,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include<Grid/Grid.h>
#include <Grid/Grid.h>
#if Nc == 3
#include <Grid/qcd/smearing/GaugeConfigurationMasked.h>
@@ -234,4 +230,5 @@ int main(int argc, char **argv)
#endif
} // main
#endif
+3 -6
View File
@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include<Grid/Grid.h>
#include <Grid/Grid.h>
int main(int argc, char **argv) {
using namespace Grid;
@@ -199,4 +195,5 @@ int main(int argc, char **argv) {
Grid_finalize();
} // main
#endif
+3 -6
View File
@@ -28,11 +28,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include<Grid/Grid.h>
#include <Grid/Grid.h>
#ifdef GRID_DEFAULT_PRECISION_DOUBLE
#define MIXED_PRECISION
@@ -453,4 +449,5 @@ int main(int argc, char **argv) {
Grid_finalize();
} // main
#endif
+3 -6
View File
@@ -28,11 +28,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include<Grid/Grid.h>
#include <Grid/Grid.h>
#ifdef GRID_DEFAULT_PRECISION_DOUBLE
#define MIXED_PRECISION
@@ -446,4 +442,5 @@ int main(int argc, char **argv) {
Grid_finalize();
} // main
#endif
+1 -7
View File
@@ -28,11 +28,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include<Grid/Grid.h>
#include <Grid/Grid.h>
using namespace Grid;
@@ -922,5 +918,3 @@ int main(int argc, char **argv) {
return 0;
#endif
} // main
#endif
+1 -7
View File
@@ -28,11 +28,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include<Grid/Grid.h>
#include <Grid/Grid.h>
using namespace Grid;
@@ -877,5 +873,3 @@ int main(int argc, char **argv) {
return 0;
#endif
} // main
#endif
+3 -6
View File
@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include<Grid/Grid.h>
#include <Grid/Grid.h>
int main(int argc, char **argv) {
using namespace Grid;
@@ -197,4 +193,5 @@ int main(int argc, char **argv) {
Grid_finalize();
} // main
#endif
+3 -6
View File
@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include<Grid/Grid.h>
#include <Grid/Grid.h>
NAMESPACE_BEGIN(Grid);
@@ -516,4 +512,5 @@ int main(int argc, char **argv) {
Grid_finalize();
} // main
#endif
+3 -6
View File
@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include<Grid/Grid.h>
#include <Grid/Grid.h>
int main(int argc, char **argv) {
using namespace Grid;
@@ -349,4 +345,5 @@ int main(int argc, char **argv) {
Grid_finalize();
} // main
#endif
+3 -6
View File
@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include<Grid/Grid.h>
#include <Grid/Grid.h>
NAMESPACE_BEGIN(Grid);
@@ -520,4 +516,5 @@ int main(int argc, char **argv) {
Grid_finalize();
} // main
#endif
+3 -6
View File
@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include<Grid/Grid.h>
#include <Grid/Grid.h>
NAMESPACE_BEGIN(Grid);
@@ -571,4 +567,5 @@ int main(int argc, char **argv) {
Grid_finalize();
} // main
#endif
+3 -6
View File
@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include<Grid/Grid.h>
#include <Grid/Grid.h>
int main(int argc, char **argv) {
using namespace Grid;
@@ -267,4 +263,5 @@ int main(int argc, char **argv) {
Grid_finalize();
} // main
#endif
+3 -6
View File
@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include<Grid/Grid.h>
#include <Grid/Grid.h>
int main(int argc, char **argv) {
using namespace Grid;
@@ -421,4 +417,5 @@ int main(int argc, char **argv) {
Grid_finalize();
} // main
#endif
+3 -6
View File
@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include<Grid/Grid.h>
#include <Grid/Grid.h>
NAMESPACE_BEGIN(Grid);
@@ -456,4 +452,5 @@ int main(int argc, char **argv) {
Grid_finalize();
} // main
#endif
+3 -6
View File
@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include<Grid/Grid.h>
#include <Grid/Grid.h>
NAMESPACE_BEGIN(Grid);
@@ -466,4 +462,5 @@ int main(int argc, char **argv) {
Grid_finalize();
} // main
#endif
+3 -6
View File
@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include<Grid/Grid.h>
#include <Grid/Grid.h>
@@ -268,4 +264,5 @@ int main(int argc, char **argv) {
Grid_finalize();
} // main
#endif
@@ -1,16 +0,0 @@
#include <Grid/Grid.h>
#pragma once
#ifndef ENABLE_FERMION_INSTANTIATIONS
#include <iostream>
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
-5
View File
@@ -26,9 +26,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
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 <Grid/Grid.h>
using namespace Grid;
@@ -734,5 +731,3 @@ int main (int argc, char ** argv)
Grid_finalize();
}
#endif
-4
View File
@@ -20,9 +20,6 @@
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 <Grid/Grid.h>
#ifdef GRID_CUDA
#define CUDA_PROFILE
@@ -442,4 +439,3 @@ 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
-6
View File
@@ -20,10 +20,6 @@
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 <Grid/Grid.h>
#ifdef GRID_CUDA
#define CUDA_PROFILE
@@ -443,5 +439,3 @@ 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
@@ -20,9 +20,6 @@
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 <Grid/Grid.h>
#ifdef GRID_CUDA
#define CUDA_PROFILE
@@ -388,5 +385,3 @@ int main (int argc, char ** argv)
Grid_finalize();
exit(0);
}
#endif
+2 -4
View File
@@ -26,9 +26,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
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 <Grid/Grid.h>
using namespace std;
@@ -241,4 +238,5 @@ void benchDw(std::vector<int> & latt4, int Ls, int threads,int report )
}
}
#endif
-5
View File
@@ -1,7 +1,3 @@
#include "disable_benchmarks_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include <Grid/Grid.h>
#include <sstream>
using namespace std;
@@ -159,4 +155,3 @@ int main (int argc, char ** argv)
Grid_finalize();
}
#endif
-5
View File
@@ -20,9 +20,6 @@
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 <Grid/Grid.h>
#ifdef GRID_CUDA
#define CUDA_PROFILE
@@ -132,5 +129,3 @@ int main (int argc, char ** argv)
Grid_finalize();
exit(0);
}
#endif
-5
View File
@@ -26,9 +26,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
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 <Grid/Grid.h>
using namespace std;
@@ -152,5 +149,3 @@ int main (int argc, char ** argv)
Grid_finalize();
}
#endif
+2 -4
View File
@@ -26,9 +26,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
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 <Grid/Grid.h>
using namespace std;
@@ -175,4 +172,5 @@ void benchDw(std::vector<int> & latt4, int Ls)
// Dw.Report();
}
#endif
-5
View File
@@ -26,9 +26,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
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 <Grid/Grid.h>
using namespace std;
@@ -113,5 +110,3 @@ int main (int argc, char ** argv)
Grid_finalize();
}
#endif
-5
View File
@@ -26,9 +26,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
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 <Grid/Grid.h>
using namespace std;
@@ -115,5 +112,3 @@ int main (int argc, char ** argv)
Grid_finalize();
}
#endif
-6
View File
@@ -26,10 +26,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
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 <Grid/Grid.h>
#include <Grid/algorithms/blas/BatchedBlas.h>
@@ -982,5 +978,3 @@ int main (int argc, char ** argv)
Grid_finalize();
fclose(FP);
}
#endif
-5
View File
@@ -26,9 +26,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
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 <Grid/Grid.h>
using namespace std;
@@ -261,5 +258,3 @@ int main (int argc, char ** argv)
Grid_finalize();
}
#endif
-5
View File
@@ -19,9 +19,6 @@ Author: Richard Rollins <rprollins@users.noreply.github.com>
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 <Grid/Grid.h>
using namespace std;
@@ -164,5 +161,3 @@ void bench_wilson_eo (
double flops = (single_site_flops * volume * ncall)/2.0;
std::cout << flops/(t1-t0) << "\t\t";
}
#endif
@@ -1,16 +0,0 @@
#include <Grid/Grid.h>
#pragma once
#ifndef ENABLE_FERMION_INSTANTIATIONS
#include <iostream>
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
+10 -13
View File
@@ -172,12 +172,6 @@ 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])
@@ -200,9 +194,6 @@ 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
@@ -418,10 +409,16 @@ AC_SEARCH_LIBS([unw_backtrace], [unwind],
[have_unwind=true],
[AC_MSG_WARN(libunwind library was not found in your system.)])
AC_SEARCH_LIBS([_Ux86_64_step], [unwind-x86_64],
[AC_DEFINE([HAVE_UNWIND_X86_64], [1], [Define to 1 if you have the `libunwind-x86_64' library])]
[have_unwind_x86_64=true],
[AC_MSG_WARN(libunwind library was not found in your system.)])
AS_CASE([$host_cpu], [x86_64],
[AC_SEARCH_LIBS([_Ux86_64_step], [unwind-x86_64],
[AC_DEFINE([HAVE_UNWIND_X86_64], [1], [Define to 1 if you have the `libunwind-x86_64' library])]
[have_unwind_x86_64=true],
[AC_MSG_WARN(libunwind library was not found in your system.)])],
[aarch64],
[AC_SEARCH_LIBS([_Uaarch64_step], [unwind-aarch64],
[AC_DEFINE([HAVE_UNWIND_AARCH64], [1], [Define to 1 if you have the `libunwind-aarch64' library])]
[have_unwind_aarch64=true],
[AC_MSG_WARN(libunwind library was not found in your system.)])])
AC_SEARCH_LIBS([SHA256_Init], [crypto],
[AC_DEFINE([HAVE_CRYPTO], [1], [Define to 1 if you have the `OpenSSL' library])]
+2 -4
View File
@@ -3,9 +3,6 @@
* without regression / tests being applied
*/
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include <Grid/Grid.h>
using namespace std;
@@ -313,4 +310,5 @@ int main (int argc, char ** argv)
Grid_finalize();
}
#endif
+2 -4
View File
@@ -3,9 +3,6 @@
* without regression / tests being applied
*/
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include <Grid/Grid.h>
using namespace std;
@@ -435,4 +432,5 @@ int main (int argc, char ** argv)
Grid_finalize();
}
#endif
+2 -4
View File
@@ -3,9 +3,6 @@
* without regression / tests being applied
*/
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include <Grid/Grid.h>
using namespace std;
@@ -538,4 +535,5 @@ int main (int argc, char ** argv)
Grid_finalize();
}
#endif
+2 -4
View File
@@ -3,9 +3,6 @@
* without regression / tests being applied
*/
#include "disable_examples_without_instantiations.h"
#ifdef ENABLE_FERMION_INSTANTIATIONS
#include <Grid/Grid.h>
using namespace std;
@@ -432,4 +429,5 @@ int main (int argc, char ** argv)
Grid_finalize();
}
#endif
@@ -1,15 +0,0 @@
#include <Grid/Grid.h>
#pragma once
#ifndef ENABLE_FERMION_INSTANTIATIONS
#include <iostream>
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
+1 -1
View File
@@ -10,7 +10,7 @@ export HDF5=/opt/cray/pe/hdf5/1.12.2.3/gnu/9.1
--disable-gparity \
--disable-fermion-reps \
--enable-shm=nvlink \
--enable-checksum-comms=no \
--enable-checksum-comms=yes \
--enable-log-views=yes \
--enable-accelerator=sycl \
--enable-accelerator-aware-mpi=no \
+3 -4
View File
@@ -25,9 +25,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
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 <Grid/Grid.h>
using namespace std;
@@ -276,6 +273,8 @@ void TestWhat(What & Ddwf,
err = phi-chi;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
}
#endif
@@ -30,9 +30,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
* 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 <Grid/Grid.h>
#include <Grid/algorithms/iterative/ImplicitlyRestartedLanczos.h>
#include <Grid/algorithms/iterative/LocalCoherenceLanczos.h>
@@ -259,4 +256,3 @@ int main (int argc, char ** argv) {
Grid_finalize();
}
#endif
-5
View File
@@ -25,9 +25,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
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 <Grid/Grid.h>
using namespace std;
@@ -240,5 +237,3 @@ int main (int argc, char ** argv)
Grid_finalize();
}
#endif
-5
View File
@@ -25,9 +25,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
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 <Grid/Grid.h>
using namespace std;
@@ -225,5 +222,3 @@ int main (int argc, char ** argv)
Grid_finalize();
}
#endif
-4
View File
@@ -25,9 +25,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
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 <Grid/Grid.h>
using namespace std;
@@ -121,4 +118,3 @@ int main (int argc, char ** argv)
Grid_finalize();
}
#endif
#endif
-4
View File
@@ -24,8 +24,6 @@ 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 <Grid/Grid.h>
#include <Grid/qcd/utils/A2Autils.h>
@@ -159,5 +157,3 @@ int main(int argc, char *argv[])
return EXIT_SUCCESS;
}
#endif
+13 -30
View File
@@ -128,10 +128,6 @@ int main (int argc, char ** argv)
typedef HermOpAdaptor<LatticeFermionD> HermFineMatrix;
HermFineMatrix FineHermOp(HermOpEO);
LatticeFermionD src(FrbGrid);
src = ComplexD(1.0);
PowerMethod<LatticeFermionD> PM; PM(HermOpEO,src);
////////////////////////////////////////////////////////////
///////////// Coarse basis and Little Dirac Operator ///////
////////////////////////////////////////////////////////////
@@ -154,7 +150,7 @@ int main (int argc, char ** argv)
std::cout << "**************************************"<<std::endl;
std::cout << "Create Subspace"<<std::endl;
std::cout << "**************************************"<<std::endl;
Aggregates.CreateSubspaceChebyshev(RNG5,HermOpEO,nbasis,35.,0.01,500);// <== last run
Aggregates.CreateSubspaceChebyshevNew(RNG5,HermOpEO,95.);
std::cout << "**************************************"<<std::endl;
std::cout << "Refine Subspace"<<std::endl;
@@ -189,7 +185,7 @@ int main (int argc, char ** argv)
std::cout << "**************************************"<<std::endl;
typedef HermitianLinearOperator<MultiGeneralCoarsenedMatrix_t,CoarseVector> MrhsHermMatrix;
Chebyshev<CoarseVector> IRLCheby(0.01,16.0,201); // 1 iter
Chebyshev<CoarseVector> IRLCheby(0.05,40.0,101); // 1 iter
MrhsHermMatrix MrhsCoarseOp (mrhs);
CoarseVector pm_src(CoarseMrhs);
@@ -197,10 +193,10 @@ int main (int argc, char ** argv)
PowerMethod<CoarseVector> cPM;
cPM(MrhsCoarseOp,pm_src);
// int Nk=16;
// int Nm=Nk*3;
int Nk=32;
int Nm=128;
int Nk=nrhs;
int Nm=Nk*3;
// int Nk=36;
// int Nm=144;
int Nstop=Nk;
int Nconv_test_interval=1;
@@ -214,7 +210,7 @@ int main (int argc, char ** argv)
nrhs,
Nk,
Nm,
1e-4,100);
1e-4,10);
int Nconv;
std::vector<RealD> eval(Nm);
@@ -235,6 +231,8 @@ int main (int argc, char ** argv)
std::cout << "**************************************"<<std::endl;
std::cout << " Recompute coarse evecs "<<std::endl;
std::cout << "**************************************"<<std::endl;
evec.resize(Nm,Coarse5d);
eval.resize(Nm);
for(int r=0;r<nrhs;r++){
random(CRNG,c_src[r]);
}
@@ -245,7 +243,7 @@ int main (int argc, char ** argv)
// Deflation guesser object
///////////////////////
std::cout << "**************************************"<<std::endl;
std::cout << " Reimport coarse evecs "<<evec.size()<<" "<<eval.size()<<std::endl;
std::cout << " Reimport coarse evecs "<<std::endl;
std::cout << "**************************************"<<std::endl;
MultiRHSDeflation<CoarseVector> MrhsGuesser;
MrhsGuesser.ImportEigenBasis(evec,eval);
@@ -254,11 +252,9 @@ int main (int argc, char ** argv)
// Extra HDCG parameters
//////////////////////////
int maxit=3000;
// ConjugateGradient<CoarseVector> CG(2.0e-1,maxit,false);
// ConjugateGradient<CoarseVector> CG(1.0e-2,maxit,false);
ConjugateGradient<CoarseVector> CG(5.0e-2,maxit,false);
RealD lo=0.2;
int ord = 7;
ConjugateGradient<CoarseVector> CG(2.0e-1,maxit,false);
RealD lo=2.0;
int ord = 9;
DoNothingGuesser<CoarseVector> DoNothing;
HPDSolver<CoarseVector> HPDSolveMrhs(MrhsCoarseOp,CG,DoNothing);
@@ -304,19 +300,6 @@ int main (int argc, char ** argv)
ConjugateGradient<LatticeFermionD> CGfine(1.0e-8,30000,false);
CGfine(HermOpEO, src, result);
}
{
std::cout << "**************************************"<<std::endl;
std::cout << "Calling MdagM CG"<<std::endl;
std::cout << "**************************************"<<std::endl;
LatticeFermion result(FGrid); result=Zero();
LatticeFermion src(FGrid); random(RNG5,src);
result=Zero();
MdagMLinearOperator<MobiusFermionD, LatticeFermionD> HermOp(Ddwf);
ConjugateGradient<LatticeFermionD> CGfine(1.0e-8,30000,false);
CGfine(HermOp, src, result);
}
#endif
Grid_finalize();
return 0;
+2 -5
View File
@@ -368,10 +368,7 @@ int main (int argc, char ** argv)
TrivialPrecon<CoarseVector> simple;
NonHermitianLinearOperator<LittleDiracOperator,CoarseVector> LinOpCoarse(LittleDiracOpPV);
// PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(1.0e-4, 100, LinOpCoarse,simple,10,10);
// PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(3.0e-2, 100, LinOpCoarse,simple,12,12); // 35 outer
// PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(5.0e-2, 100, LinOpCoarse,simple,12,12); // 36 outer, 12s
// PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(1.0e-1, 100, LinOpCoarse,simple,12,12); // 36 ; 11s
PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(3.0e-1, 100, LinOpCoarse,simple,12,12);
PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(3.0e-2, 100, LinOpCoarse,simple,10,10);
L2PGCR.Level(3);
c_res=Zero();
L2PGCR(c_src,c_res);
@@ -403,7 +400,7 @@ int main (int argc, char ** argv)
LinOpCoarse,
L2PGCR);
PrecGeneralisedConjugateResidualNonHermitian<LatticeFermion> L1PGCR(1.0e-8,100,PVdagM,TwoLevelPrecon,10,10);
PrecGeneralisedConjugateResidualNonHermitian<LatticeFermion> L1PGCR(1.0e-8,1000,PVdagM,TwoLevelPrecon,16,16);
L1PGCR.Level(1);
f_res=Zero();
@@ -1,493 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_padded_cell.cc
Copyright (C) 2023
Author: Peter Boyle <paboyle@ph.ed.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
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
#include <Grid/lattice/PaddedCell.h>
#include <Grid/stencil/GeneralLocalStencil.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
#include <Grid/algorithms/iterative/BiCGSTAB.h>
using namespace std;
using namespace Grid;
template<class Matrix,class Field>
class PVdagMLinearOperator : public LinearOperatorBase<Field> {
Matrix &_Mat;
Matrix &_PV;
public:
PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV){};
void OpDiag (const Field &in, Field &out) { assert(0); }
void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); }
void OpDirAll (const Field &in, std::vector<Field> &out){ assert(0); };
void Op (const Field &in, Field &out){
// std::cout << GridLogMessage<< "Op: PVdag M "<<std::endl;
Field tmp(in.Grid());
_Mat.M(in,tmp);
_PV.Mdag(tmp,out);
}
void AdjOp (const Field &in, Field &out){
// std::cout << GridLogMessage<<"AdjOp: Mdag PV "<<std::endl;
Field tmp(in.Grid());
_PV.M(in,tmp);
_Mat.Mdag(tmp,out);
}
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
assert(0);
}
void HermOp(const Field &in, Field &out){
// std::cout <<GridLogMessage<< "HermOp: Mdag PV PVdag M"<<std::endl;
Field tmp(in.Grid());
Op(in,tmp);
AdjOp(tmp,out);
// std::cout << "HermOp done "<<norm2(out)<<std::endl;
}
};
template<class Matrix,class Field>
class MdagPVLinearOperator : public LinearOperatorBase<Field> {
Matrix &_Mat;
Matrix &_PV;
public:
MdagPVLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV){};
void OpDiag (const Field &in, Field &out) { assert(0); }
void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); }
void OpDirAll (const Field &in, std::vector<Field> &out){ assert(0); };
void Op (const Field &in, Field &out){
Field tmp(in.Grid());
// std::cout <<GridLogMessage<< "Op: PVdag M "<<std::endl;
_PV.M(in,tmp);
_Mat.Mdag(tmp,out);
}
void AdjOp (const Field &in, Field &out){
// std::cout <<GridLogMessage<< "AdjOp: Mdag PV "<<std::endl;
Field tmp(in.Grid());
_Mat.M(in,tmp);
_PV.Mdag(tmp,out);
}
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
assert(0);
}
void HermOp(const Field &in, Field &out){
// std::cout << GridLogMessage<<"HermOp: PVdag M Mdag PV "<<std::endl;
Field tmp(in.Grid());
Op(in,tmp);
AdjOp(tmp,out);
// std::cout << "HermOp done "<<norm2(out)<<std::endl;
}
};
template<class Matrix,class Field>
class ShiftedPVdagMLinearOperator : public LinearOperatorBase<Field> {
Matrix &_Mat;
Matrix &_PV;
RealD shift;
public:
ShiftedPVdagMLinearOperator(RealD _shift,Matrix &Mat,Matrix &PV): shift(_shift),_Mat(Mat),_PV(PV){};
void OpDiag (const Field &in, Field &out) { assert(0); }
void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); }
void OpDirAll (const Field &in, std::vector<Field> &out){ assert(0); };
void Op (const Field &in, Field &out){
// std::cout << "Op: PVdag M "<<std::endl;
Field tmp(in.Grid());
_Mat.M(in,tmp);
_PV.Mdag(tmp,out);
out = out + shift * in;
}
void AdjOp (const Field &in, Field &out){
// std::cout << "AdjOp: Mdag PV "<<std::endl;
Field tmp(in.Grid());
_PV.M(tmp,out);
_Mat.Mdag(in,tmp);
out = out + shift * in;
}
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
void HermOp(const Field &in, Field &out){
// std::cout << "HermOp: Mdag PV PVdag M"<<std::endl;
Field tmp(in.Grid());
Op(in,tmp);
AdjOp(tmp,out);
}
};
template<class Fobj,class CComplex,int nbasis>
class MGPreconditionerSVD : public LinearFunction< Lattice<Fobj> > {
public:
using LinearFunction<Lattice<Fobj> >::operator();
typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField FineField;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
typedef LinearOperatorBase<FineField> FineOperator;
typedef LinearFunction <FineField> FineSmoother;
typedef LinearOperatorBase<CoarseVector> CoarseOperator;
typedef LinearFunction <CoarseVector> CoarseSolver;
Aggregates & _FineToCoarse;
Aggregates & _CoarseToFine;
FineOperator & _FineOperator;
FineSmoother & _PreSmoother;
FineSmoother & _PostSmoother;
CoarseOperator & _CoarseOperator;
CoarseSolver & _CoarseSolve;
int level; void Level(int lv) {level = lv; };
MGPreconditionerSVD(Aggregates &FtoC,
Aggregates &CtoF,
FineOperator &Fine,
FineSmoother &PreSmoother,
FineSmoother &PostSmoother,
CoarseOperator &CoarseOperator_,
CoarseSolver &CoarseSolve_)
: _FineToCoarse(FtoC),
_CoarseToFine(CtoF),
_FineOperator(Fine),
_PreSmoother(PreSmoother),
_PostSmoother(PostSmoother),
_CoarseOperator(CoarseOperator_),
_CoarseSolve(CoarseSolve_),
level(1) { }
virtual void operator()(const FineField &in, FineField & out)
{
GridBase *CoarseGrid = _FineToCoarse.CoarseGrid;
// auto CoarseGrid = _CoarseOperator.Grid();
CoarseVector Csrc(CoarseGrid);
CoarseVector Csol(CoarseGrid);
FineField vec1(in.Grid());
FineField vec2(in.Grid());
std::cout<<GridLogMessage << "Calling PreSmoother " <<std::endl;
// std::cout<<GridLogMessage << "Calling PreSmoother input residual "<<norm2(in) <<std::endl;
double t;
// Fine Smoother
// out = in;
out = Zero();
t=-usecond();
_PreSmoother(in,out);
t+=usecond();
std::cout<<GridLogMessage << "PreSmoother took "<< t/1000.0<< "ms" <<std::endl;
// Update the residual
_FineOperator.Op(out,vec1); sub(vec1, in ,vec1);
// std::cout<<GridLogMessage <<"Residual-1 now " <<norm2(vec1)<<std::endl;
// Fine to Coarse
t=-usecond();
_FineToCoarse.ProjectToSubspace (Csrc,vec1);
t+=usecond();
std::cout<<GridLogMessage << "Project to coarse took "<< t/1000.0<< "ms" <<std::endl;
// Coarse correction
t=-usecond();
Csol = Zero();
_CoarseSolve(Csrc,Csol);
//Csol=Zero();
t+=usecond();
std::cout<<GridLogMessage << "Coarse solve took "<< t/1000.0<< "ms" <<std::endl;
// Coarse to Fine
t=-usecond();
// _CoarseOperator.PromoteFromSubspace(_Aggregates,Csol,vec1);
_CoarseToFine.PromoteFromSubspace(Csol,vec1);
add(out,out,vec1);
t+=usecond();
std::cout<<GridLogMessage << "Promote to this level took "<< t/1000.0<< "ms" <<std::endl;
// Residual
_FineOperator.Op(out,vec1); sub(vec1 ,in , vec1);
// std::cout<<GridLogMessage <<"Residual-2 now " <<norm2(vec1)<<std::endl;
// Fine Smoother
t=-usecond();
// vec2=vec1;
vec2=Zero();
_PostSmoother(vec1,vec2);
t+=usecond();
std::cout<<GridLogMessage << "PostSmoother took "<< t/1000.0<< "ms" <<std::endl;
add( out,out,vec2);
std::cout<<GridLogMessage << "Done " <<std::endl;
}
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
const int Ls=16;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
// Construct a coarsened grid
Coordinate clatt = GridDefaultLatt();
for(int d=0;d<clatt.size();d++){
clatt[d] = clatt[d]/2;
// clatt[d] = clatt[d]/4;
}
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
std::vector<int> cseeds({5,6,7,8});
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
LatticeFermion src(FGrid); random(RNG5,src);
LatticeFermion result(FGrid); result=Zero();
LatticeFermion ref(FGrid); ref=Zero();
LatticeFermion tmp(FGrid);
LatticeFermion err(FGrid);
LatticeGaugeField Umu(UGrid);
FieldMetaData header;
std::string file("ckpoint_lat.4000");
NerscIO::readConfiguration(Umu,header,file);
RealD mass=0.01;
RealD M5=1.8;
DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
DomainWallFermionD Dpv(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0,M5);
const int nbasis = 30;
const int cb = 0 ;
NextToNearestStencilGeometry5D geom(Coarse5d);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
typedef PVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> PVdagM_t;
typedef MdagPVLinearOperator<DomainWallFermionD,LatticeFermionD> MdagPV_t;
typedef ShiftedPVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> ShiftedPVdagM_t;
PVdagM_t PVdagM(Ddwf,Dpv);
MdagPV_t MdagPV(Ddwf,Dpv);
// ShiftedPVdagM_t ShiftedPVdagM(2.0,Ddwf,Dpv); // 355
// ShiftedPVdagM_t ShiftedPVdagM(1.0,Ddwf,Dpv); // 246
// ShiftedPVdagM_t ShiftedPVdagM(0.5,Ddwf,Dpv); // 183
// ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // 145
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 134
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 127 -- NULL space via inverse iteration
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 57 -- NULL space via inverse iteration; 3 iterations
// ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // 57 , tighter inversion
// ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // nbasis 20 -- 49 iters
// ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // nbasis 20 -- 70 iters; asymmetric
// ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // 58; Loosen coarse, tighten fine
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 56 ...
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 51 ... with 24 vecs
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 31 ... with 24 vecs and 2^4 blocking
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 43 ... with 16 vecs and 2^4 blocking, sloppier
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 35 ... with 20 vecs and 2^4 blocking
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 35 ... with 20 vecs and 2^4 blocking, looser coarse
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 64 ... with 20 vecs, Christoph setup, and 2^4 blocking, looser coarse
ShiftedPVdagM_t ShiftedPVdagM(0.01,Ddwf,Dpv); //
// Run power method on HOA??
PowerMethod<LatticeFermion> PM;
// PM(PVdagM,src);
// PM(MdagPV,src);
// Warning: This routine calls PVdagM.Op, not PVdagM.HermOp
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
Subspace V(Coarse5d,FGrid,cb);
Subspace U(Coarse5d,FGrid,cb);
// Breeds right singular vectors with call to HermOp (V)
V.CreateSubspaceChebyshev(RNG5,PVdagM,
nbasis,
4000.0,0.003,
500);
// Breeds left singular vectors with call to HermOp (U)
// U.CreateSubspaceChebyshev(RNG5,PVdagM,
U.CreateSubspaceChebyshev(RNG5,MdagPV,
nbasis,
4000.0,0.003,
500);
typedef Aggregation<vSpinColourVector,vTComplex,2*nbasis> CombinedSubspace;
CombinedSubspace CombinedUV(Coarse5d,FGrid,cb);
for(int b=0;b<nbasis;b++){
CombinedUV.subspace[b] = V.subspace[b];
CombinedUV.subspace[b+nbasis] = U.subspace[b];
}
int bl, br;
std::cout <<" <V| PVdagM| V> " <<std::endl;
for(bl=0;bl<nbasis;bl++){
for(br=0;br<nbasis;br++){
PVdagM.Op(V.subspace[br],src);
std::cout <<bl<<" "<<br<<"\t"<<innerProduct(V.subspace[bl],src)<<std::endl;
}}
std::cout <<" <V| PVdagM| U> " <<std::endl;
for(bl=0;bl<nbasis;bl++){
for(br=0;br<nbasis;br++){
PVdagM.Op(U.subspace[br],src);
std::cout <<bl<<" "<<br<<"\t"<<innerProduct(V.subspace[bl],src)<<std::endl;
}}
std::cout <<" <U| PVdagM| V> " <<std::endl;
for(bl=0;bl<nbasis;bl++){
for(br=0;br<nbasis;br++){
PVdagM.Op(V.subspace[br],src);
std::cout <<bl<<" "<<br<<"\t"<<innerProduct(U.subspace[bl],src)<<std::endl;
}}
std::cout <<" <U| PVdagM| U> " <<std::endl;
for(bl=0;bl<nbasis;bl++){
for(br=0;br<nbasis;br++){
PVdagM.Op(U.subspace[br],src);
std::cout <<bl<<" "<<br<<"\t"<<innerProduct(U.subspace[bl],src)<<std::endl;
}}
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperatorV;
typedef LittleDiracOperatorV::CoarseVector CoarseVectorV;
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,2*nbasis> LittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
V.Orthogonalise();
for(int b =0 ; b<nbasis;b++){
CoarseVectorV c_src (Coarse5d);
V.ProjectToSubspace (c_src,U.subspace[b]);
V.PromoteFromSubspace(c_src,src);
std::cout << " Completeness of U in V ["<< b<<"] "<< std::sqrt(norm2(src)/norm2(U.subspace[b]))<<std::endl;
}
CoarseVector c_src (Coarse5d);
CoarseVector c_res (Coarse5d);
CoarseVector c_proj(Coarse5d);
LittleDiracOperator LittleDiracOpPV(geom,FGrid,Coarse5d);
LittleDiracOpPV.CoarsenOperator(PVdagM,CombinedUV,CombinedUV);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"Testing coarsened operator "<<std::endl;
Complex one(1.0);
c_src = one; // 1 in every element for vector 1.
blockPromote(c_src,err,CombinedUV.subspace);
LatticeFermion prom(FGrid);
prom=Zero();
for(int b=0;b<nbasis*2;b++){
prom=prom+CombinedUV.subspace[b];
}
std::cout<<GridLogMessage<<"c_src "<<norm2(c_src)<<std::endl;
std::cout<<GridLogMessage<<"prom "<<norm2(prom)<<std::endl;
PVdagM.Op(prom,tmp);
blockProject(c_proj,tmp,CombinedUV.subspace);
std::cout<<GridLogMessage<<" Called Big Dirac Op "<<norm2(tmp)<<std::endl;
LittleDiracOpPV.M(c_src,c_res);
std::cout<<GridLogMessage<<" Called Little Dirac Op c_src "<< norm2(c_src) << " c_res "<< norm2(c_res) <<std::endl;
std::cout<<GridLogMessage<<"Little dop : "<<norm2(c_res)<<std::endl;
std::cout<<GridLogMessage<<"Big dop in subspace : "<<norm2(c_proj)<<std::endl;
c_proj = c_proj - c_res;
std::cout<<GridLogMessage<<" ldop error: "<<norm2(c_proj)<<std::endl;
/**********
* Some solvers
**********
*/
///////////////////////////////////////
// Coarse grid solver test
///////////////////////////////////////
std::cout<<GridLogMessage<<"******************* "<<std::endl;
std::cout<<GridLogMessage<<" Coarse Grid Solve -- Level 3 "<<std::endl;
std::cout<<GridLogMessage<<"******************* "<<std::endl;
TrivialPrecon<CoarseVector> simple;
NonHermitianLinearOperator<LittleDiracOperator,CoarseVector> LinOpCoarse(LittleDiracOpPV);
// PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(1.0e-4, 100, LinOpCoarse,simple,10,10);
PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(1.0e-2, 10, LinOpCoarse,simple,20,20);
L2PGCR.Level(3);
c_res=Zero();
L2PGCR(c_src,c_res);
////////////////////////////////////////
// Fine grid smoother
////////////////////////////////////////
std::cout<<GridLogMessage<<"******************* "<<std::endl;
std::cout<<GridLogMessage<<" Fine Grid Smoother -- Level 2 "<<std::endl;
std::cout<<GridLogMessage<<"******************* "<<std::endl;
TrivialPrecon<LatticeFermionD> simple_fine;
// NonHermitianLinearOperator<PVdagM_t,LatticeFermionD> LinOpSmooth(PVdagM);
PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> SmootherGCR(0.01,1,ShiftedPVdagM,simple_fine,16,16);
SmootherGCR.Level(2);
LatticeFermionD f_src(FGrid);
LatticeFermionD f_res(FGrid);
f_src = one; // 1 in every element for vector 1.
f_res=Zero();
SmootherGCR(f_src,f_res);
typedef MGPreconditionerSVD<vSpinColourVector, vTComplex,nbasis*2> TwoLevelMG;
TwoLevelMG TwoLevelPrecon(CombinedUV,CombinedUV,
PVdagM,
simple_fine,
SmootherGCR,
LinOpCoarse,
L2PGCR);
PrecGeneralisedConjugateResidualNonHermitian<LatticeFermion> L1PGCR(1.0e-8,1000,PVdagM,TwoLevelPrecon,20,20);
L1PGCR.Level(1);
f_res=Zero();
L1PGCR(f_src,f_res);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage << "Done "<< std::endl;
Grid_finalize();
return 0;
}
@@ -1,492 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_padded_cell.cc
Copyright (C) 2023
Author: Peter Boyle <paboyle@ph.ed.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
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
#include <Grid/lattice/PaddedCell.h>
#include <Grid/stencil/GeneralLocalStencil.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
#include <Grid/algorithms/iterative/BiCGSTAB.h>
using namespace std;
using namespace Grid;
template<class Matrix,class Field>
class PVdagMLinearOperator : public LinearOperatorBase<Field> {
Matrix &_Mat;
Matrix &_PV;
public:
PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV){};
void OpDiag (const Field &in, Field &out) { assert(0); }
void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); }
void OpDirAll (const Field &in, std::vector<Field> &out){ assert(0); };
void Op (const Field &in, Field &out){
// std::cout << GridLogMessage<< "Op: PVdag M "<<std::endl;
Field tmp(in.Grid());
_Mat.M(in,tmp);
_PV.Mdag(tmp,out);
}
void AdjOp (const Field &in, Field &out){
// std::cout << GridLogMessage<<"AdjOp: Mdag PV "<<std::endl;
Field tmp(in.Grid());
_PV.M(in,tmp);
_Mat.Mdag(tmp,out);
}
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
HermOp(in,out);
ComplexD dot = innerProduct(in,out);
n1=real(dot);
n2=norm2(out);
}
void HermOp(const Field &in, Field &out){
// std::cout <<GridLogMessage<< "HermOp: Mdag PV PVdag M"<<std::endl;
Field tmp(in.Grid());
Op(in,tmp);
AdjOp(tmp,out);
// std::cout << "HermOp done "<<norm2(out)<<std::endl;
}
};
template<class Matrix,class Field>
class MdagPVLinearOperator : public LinearOperatorBase<Field> {
Matrix &_Mat;
Matrix &_PV;
public:
MdagPVLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV){};
void OpDiag (const Field &in, Field &out) { assert(0); }
void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); }
void OpDirAll (const Field &in, std::vector<Field> &out){ assert(0); };
void Op (const Field &in, Field &out){
Field tmp(in.Grid());
// std::cout <<GridLogMessage<< "Op: PVdag M "<<std::endl;
_PV.M(in,tmp);
_Mat.Mdag(tmp,out);
}
void AdjOp (const Field &in, Field &out){
// std::cout <<GridLogMessage<< "AdjOp: Mdag PV "<<std::endl;
Field tmp(in.Grid());
_Mat.M(in,tmp);
_PV.Mdag(tmp,out);
}
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
ComplexD dot = innerProduct(in,out);
n1=real(dot);
n2=norm2(out);
}
void HermOp(const Field &in, Field &out){
// std::cout << GridLogMessage<<"HermOp: PVdag M Mdag PV "<<std::endl;
Field tmp(in.Grid());
Op(in,tmp);
AdjOp(tmp,out);
// std::cout << "HermOp done "<<norm2(out)<<std::endl;
}
};
template<class Matrix,class Field>
class ShiftedPVdagMLinearOperator : public LinearOperatorBase<Field> {
Matrix &_Mat;
Matrix &_PV;
RealD shift;
public:
ShiftedPVdagMLinearOperator(RealD _shift,Matrix &Mat,Matrix &PV): shift(_shift),_Mat(Mat),_PV(PV){};
void OpDiag (const Field &in, Field &out) { assert(0); }
void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); }
void OpDirAll (const Field &in, std::vector<Field> &out){ assert(0); };
void Op (const Field &in, Field &out){
// std::cout << "Op: PVdag M "<<std::endl;
Field tmp(in.Grid());
_Mat.M(in,tmp);
_PV.Mdag(tmp,out);
out = out + shift * in;
}
void AdjOp (const Field &in, Field &out){
// std::cout << "AdjOp: Mdag PV "<<std::endl;
Field tmp(in.Grid());
_PV.M(tmp,out);
_Mat.Mdag(in,tmp);
out = out + shift * in;
}
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
void HermOp(const Field &in, Field &out){
// std::cout << "HermOp: Mdag PV PVdag M"<<std::endl;
Field tmp(in.Grid());
Op(in,tmp);
AdjOp(tmp,out);
}
};
template<class Fobj,class CComplex,int nbasis>
class MGPreconditionerSVD : public LinearFunction< Lattice<Fobj> > {
public:
using LinearFunction<Lattice<Fobj> >::operator();
typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField FineField;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
typedef LinearOperatorBase<FineField> FineOperator;
typedef LinearFunction <FineField> FineSmoother;
typedef LinearOperatorBase<CoarseVector> CoarseOperator;
typedef LinearFunction <CoarseVector> CoarseSolver;
Aggregates & _FineToCoarse;
Aggregates & _CoarseToFine;
FineOperator & _FineOperator;
FineSmoother & _PreSmoother;
FineSmoother & _PostSmoother;
CoarseOperator & _CoarseOperator;
CoarseSolver & _CoarseSolve;
int level; void Level(int lv) {level = lv; };
MGPreconditionerSVD(Aggregates &FtoC,
Aggregates &CtoF,
FineOperator &Fine,
FineSmoother &PreSmoother,
FineSmoother &PostSmoother,
CoarseOperator &CoarseOperator_,
CoarseSolver &CoarseSolve_)
: _FineToCoarse(FtoC),
_CoarseToFine(CtoF),
_FineOperator(Fine),
_PreSmoother(PreSmoother),
_PostSmoother(PostSmoother),
_CoarseOperator(CoarseOperator_),
_CoarseSolve(CoarseSolve_),
level(1) { }
virtual void operator()(const FineField &in, FineField & out)
{
GridBase *CoarseGrid = _FineToCoarse.CoarseGrid;
// auto CoarseGrid = _CoarseOperator.Grid();
CoarseVector Csrc(CoarseGrid);
CoarseVector Csol(CoarseGrid);
FineField vec1(in.Grid());
FineField vec2(in.Grid());
std::cout<<GridLogMessage << "Calling PreSmoother " <<std::endl;
// std::cout<<GridLogMessage << "Calling PreSmoother input residual "<<norm2(in) <<std::endl;
double t;
// Fine Smoother
// out = in;
out = Zero();
t=-usecond();
_PreSmoother(in,out);
t+=usecond();
std::cout<<GridLogMessage << "PreSmoother took "<< t/1000.0<< "ms" <<std::endl;
// Update the residual
_FineOperator.Op(out,vec1); sub(vec1, in ,vec1);
// std::cout<<GridLogMessage <<"Residual-1 now " <<norm2(vec1)<<std::endl;
// Fine to Coarse
t=-usecond();
_FineToCoarse.ProjectToSubspace (Csrc,vec1);
t+=usecond();
std::cout<<GridLogMessage << "Project to coarse took "<< t/1000.0<< "ms" <<std::endl;
// Coarse correction
t=-usecond();
Csol = Zero();
_CoarseSolve(Csrc,Csol);
//Csol=Zero();
t+=usecond();
std::cout<<GridLogMessage << "Coarse solve took "<< t/1000.0<< "ms" <<std::endl;
// Coarse to Fine
t=-usecond();
// _CoarseOperator.PromoteFromSubspace(_Aggregates,Csol,vec1);
_CoarseToFine.PromoteFromSubspace(Csol,vec1);
add(out,out,vec1);
t+=usecond();
std::cout<<GridLogMessage << "Promote to this level took "<< t/1000.0<< "ms" <<std::endl;
// Residual
_FineOperator.Op(out,vec1); sub(vec1 ,in , vec1);
// std::cout<<GridLogMessage <<"Residual-2 now " <<norm2(vec1)<<std::endl;
// Fine Smoother
t=-usecond();
// vec2=vec1;
vec2=Zero();
_PostSmoother(vec1,vec2);
t+=usecond();
std::cout<<GridLogMessage << "PostSmoother took "<< t/1000.0<< "ms" <<std::endl;
add( out,out,vec2);
std::cout<<GridLogMessage << "Done " <<std::endl;
}
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
const int Ls=16;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
// Construct a coarsened grid
Coordinate clatt = GridDefaultLatt();
for(int d=0;d<clatt.size();d++){
clatt[d] = clatt[d]/2;
// clatt[d] = clatt[d]/4;
}
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
std::vector<int> cseeds({5,6,7,8});
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
LatticeFermion src(FGrid); random(RNG5,src);
LatticeFermion result(FGrid); result=Zero();
LatticeFermion ref(FGrid); ref=Zero();
LatticeFermion tmp(FGrid);
LatticeFermion err(FGrid);
LatticeGaugeField Umu(UGrid);
FieldMetaData header;
std::string file("ckpoint_lat.4000");
NerscIO::readConfiguration(Umu,header,file);
RealD mass=0.01;
RealD M5=1.8;
DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
DomainWallFermionD Dpv(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0,M5);
const int nbasis = 20;
const int cb = 0 ;
NextToNearestStencilGeometry5D geom(Coarse5d);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
typedef PVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> PVdagM_t;
typedef MdagPVLinearOperator<DomainWallFermionD,LatticeFermionD> MdagPV_t;
typedef ShiftedPVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> ShiftedPVdagM_t;
PVdagM_t PVdagM(Ddwf,Dpv);
MdagPV_t MdagPV(Ddwf,Dpv);
// ShiftedPVdagM_t ShiftedPVdagM(2.0,Ddwf,Dpv); // 355
// ShiftedPVdagM_t ShiftedPVdagM(1.0,Ddwf,Dpv); // 246
// ShiftedPVdagM_t ShiftedPVdagM(0.5,Ddwf,Dpv); // 183
// ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // 145
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 134
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 127 -- NULL space via inverse iteration
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 57 -- NULL space via inverse iteration; 3 iterations
// ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // 57 , tighter inversion
// ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // nbasis 20 -- 49 iters
// ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // nbasis 20 -- 70 iters; asymmetric
// ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // 58; Loosen coarse, tighten fine
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 56 ...
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 51 ... with 24 vecs
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 31 ... with 24 vecs and 2^4 blocking
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 43 ... with 16 vecs and 2^4 blocking, sloppier
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 35 ... with 20 vecs and 2^4 blocking
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 35 ... with 20 vecs and 2^4 blocking, looser coarse
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 64 ... with 20 vecs, Christoph setup, and 2^4 blocking, looser coarse
ShiftedPVdagM_t ShiftedPVdagM(0.01,Ddwf,Dpv); //
// Run power method on HOA??
PowerMethod<LatticeFermion> PM;
// PM(PVdagM,src);
// PM(MdagPV,src);
// Warning: This routine calls PVdagM.Op, not PVdagM.HermOp
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
Subspace V(Coarse5d,FGrid,cb);
Subspace U(Coarse5d,FGrid,cb);
// Breeds right singular vectors with call to HermOp (V)
V.CreateSubspace(RNG5,PVdagM,nbasis);
// Breeds left singular vectors with call to HermOp (U)
// U.CreateSubspaceChebyshev(RNG5,MdagPV,
U.CreateSubspace(RNG5,PVdagM,nbasis);
typedef Aggregation<vSpinColourVector,vTComplex,2*nbasis> CombinedSubspace;
CombinedSubspace CombinedUV(Coarse5d,FGrid,cb);
for(int b=0;b<nbasis;b++){
CombinedUV.subspace[b] = V.subspace[b];
CombinedUV.subspace[b+nbasis] = U.subspace[b];
}
int bl, br;
std::cout <<" <V| PVdagM| V> " <<std::endl;
for(bl=0;bl<nbasis;bl++){
for(br=0;br<nbasis;br++){
PVdagM.Op(V.subspace[br],src);
std::cout <<bl<<" "<<br<<"\t"<<innerProduct(V.subspace[bl],src)<<std::endl;
}}
std::cout <<" <V| PVdagM| U> " <<std::endl;
for(bl=0;bl<nbasis;bl++){
for(br=0;br<nbasis;br++){
PVdagM.Op(U.subspace[br],src);
std::cout <<bl<<" "<<br<<"\t"<<innerProduct(V.subspace[bl],src)<<std::endl;
}}
std::cout <<" <U| PVdagM| V> " <<std::endl;
for(bl=0;bl<nbasis;bl++){
for(br=0;br<nbasis;br++){
PVdagM.Op(V.subspace[br],src);
std::cout <<bl<<" "<<br<<"\t"<<innerProduct(U.subspace[bl],src)<<std::endl;
}}
std::cout <<" <U| PVdagM| U> " <<std::endl;
for(bl=0;bl<nbasis;bl++){
for(br=0;br<nbasis;br++){
PVdagM.Op(U.subspace[br],src);
std::cout <<bl<<" "<<br<<"\t"<<innerProduct(U.subspace[bl],src)<<std::endl;
}}
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperatorV;
typedef LittleDiracOperatorV::CoarseVector CoarseVectorV;
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,2*nbasis> LittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
V.Orthogonalise();
for(int b =0 ; b<nbasis;b++){
CoarseVectorV c_src (Coarse5d);
V.ProjectToSubspace (c_src,U.subspace[b]);
V.PromoteFromSubspace(c_src,src);
std::cout << " Completeness of U in V ["<< b<<"] "<< std::sqrt(norm2(src)/norm2(U.subspace[b]))<<std::endl;
}
CoarseVector c_src (Coarse5d);
CoarseVector c_res (Coarse5d);
CoarseVector c_proj(Coarse5d);
LittleDiracOperator LittleDiracOpPV(geom,FGrid,Coarse5d);
LittleDiracOpPV.CoarsenOperator(PVdagM,CombinedUV,CombinedUV);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"Testing coarsened operator "<<std::endl;
Complex one(1.0);
c_src = one; // 1 in every element for vector 1.
blockPromote(c_src,err,CombinedUV.subspace);
LatticeFermion prom(FGrid);
prom=Zero();
for(int b=0;b<nbasis*2;b++){
prom=prom+CombinedUV.subspace[b];
}
std::cout<<GridLogMessage<<"c_src "<<norm2(c_src)<<std::endl;
std::cout<<GridLogMessage<<"prom "<<norm2(prom)<<std::endl;
PVdagM.Op(prom,tmp);
blockProject(c_proj,tmp,CombinedUV.subspace);
std::cout<<GridLogMessage<<" Called Big Dirac Op "<<norm2(tmp)<<std::endl;
LittleDiracOpPV.M(c_src,c_res);
std::cout<<GridLogMessage<<" Called Little Dirac Op c_src "<< norm2(c_src) << " c_res "<< norm2(c_res) <<std::endl;
std::cout<<GridLogMessage<<"Little dop : "<<norm2(c_res)<<std::endl;
std::cout<<GridLogMessage<<"Big dop in subspace : "<<norm2(c_proj)<<std::endl;
c_proj = c_proj - c_res;
std::cout<<GridLogMessage<<" ldop error: "<<norm2(c_proj)<<std::endl;
/**********
* Some solvers
**********
*/
///////////////////////////////////////
// Coarse grid solver test
///////////////////////////////////////
std::cout<<GridLogMessage<<"******************* "<<std::endl;
std::cout<<GridLogMessage<<" Coarse Grid Solve -- Level 3 "<<std::endl;
std::cout<<GridLogMessage<<"******************* "<<std::endl;
TrivialPrecon<CoarseVector> simple;
NonHermitianLinearOperator<LittleDiracOperator,CoarseVector> LinOpCoarse(LittleDiracOpPV);
// PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(1.0e-4, 100, LinOpCoarse,simple,10,10);
PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(1.0e-2, 10, LinOpCoarse,simple,20,20);
L2PGCR.Level(3);
c_res=Zero();
L2PGCR(c_src,c_res);
////////////////////////////////////////
// Fine grid smoother
////////////////////////////////////////
std::cout<<GridLogMessage<<"******************* "<<std::endl;
std::cout<<GridLogMessage<<" Fine Grid Smoother -- Level 2 "<<std::endl;
std::cout<<GridLogMessage<<"******************* "<<std::endl;
TrivialPrecon<LatticeFermionD> simple_fine;
// NonHermitianLinearOperator<PVdagM_t,LatticeFermionD> LinOpSmooth(PVdagM);
PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> SmootherGCR(0.01,1,ShiftedPVdagM,simple_fine,16,16);
SmootherGCR.Level(2);
LatticeFermionD f_src(FGrid);
LatticeFermionD f_res(FGrid);
f_src = one; // 1 in every element for vector 1.
f_res=Zero();
SmootherGCR(f_src,f_res);
typedef MGPreconditionerSVD<vSpinColourVector, vTComplex,nbasis*2> TwoLevelMG;
TwoLevelMG TwoLevelPrecon(CombinedUV,CombinedUV,
PVdagM,
simple_fine,
SmootherGCR,
LinOpCoarse,
L2PGCR);
PrecGeneralisedConjugateResidualNonHermitian<LatticeFermion> L1PGCR(1.0e-8,1000,PVdagM,TwoLevelPrecon,20,20);
L1PGCR.Level(1);
f_res=Zero();
L1PGCR(f_src,f_res);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage << "Done "<< std::endl;
Grid_finalize();
return 0;
}
@@ -1,479 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_padded_cell.cc
Copyright (C) 2023
Author: Peter Boyle <paboyle@ph.ed.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
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
#include <Grid/lattice/PaddedCell.h>
#include <Grid/stencil/GeneralLocalStencil.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
#include <Grid/algorithms/iterative/BiCGSTAB.h>
using namespace std;
using namespace Grid;
template<class Matrix,class Field>
class PVdagMLinearOperator : public LinearOperatorBase<Field> {
Matrix &_Mat;
Matrix &_PV;
public:
PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV){};
void OpDiag (const Field &in, Field &out) { assert(0); }
void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); }
void OpDirAll (const Field &in, std::vector<Field> &out){ assert(0); };
void Op (const Field &in, Field &out){
// std::cout << GridLogMessage<< "Op: PVdag M "<<std::endl;
Field tmp(in.Grid());
_Mat.M(in,tmp);
_PV.Mdag(tmp,out);
}
void AdjOp (const Field &in, Field &out){
// std::cout << GridLogMessage<<"AdjOp: Mdag PV "<<std::endl;
Field tmp(in.Grid());
_PV.M(in,tmp);
_Mat.Mdag(tmp,out);
}
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
assert(0);
}
void HermOp(const Field &in, Field &out){
// std::cout <<GridLogMessage<< "HermOp: Mdag PV PVdag M"<<std::endl;
Field tmp(in.Grid());
Op(in,tmp);
AdjOp(tmp,out);
// std::cout << "HermOp done "<<norm2(out)<<std::endl;
}
};
template<class Matrix,class Field>
class MdagPVLinearOperator : public LinearOperatorBase<Field> {
Matrix &_Mat;
Matrix &_PV;
public:
MdagPVLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV){};
void OpDiag (const Field &in, Field &out) { assert(0); }
void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); }
void OpDirAll (const Field &in, std::vector<Field> &out){ assert(0); };
void Op (const Field &in, Field &out){
Field tmp(in.Grid());
// std::cout <<GridLogMessage<< "Op: PVdag M "<<std::endl;
_PV.M(in,tmp);
_Mat.Mdag(tmp,out);
}
void AdjOp (const Field &in, Field &out){
// std::cout <<GridLogMessage<< "AdjOp: Mdag PV "<<std::endl;
Field tmp(in.Grid());
_Mat.M(in,tmp);
_PV.Mdag(tmp,out);
}
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
assert(0);
}
void HermOp(const Field &in, Field &out){
// std::cout << GridLogMessage<<"HermOp: PVdag M Mdag PV "<<std::endl;
Field tmp(in.Grid());
Op(in,tmp);
AdjOp(tmp,out);
// std::cout << "HermOp done "<<norm2(out)<<std::endl;
}
};
template<class Matrix,class Field>
class ShiftedPVdagMLinearOperator : public LinearOperatorBase<Field> {
Matrix &_Mat;
Matrix &_PV;
RealD shift;
public:
ShiftedPVdagMLinearOperator(RealD _shift,Matrix &Mat,Matrix &PV): shift(_shift),_Mat(Mat),_PV(PV){};
void OpDiag (const Field &in, Field &out) { assert(0); }
void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); }
void OpDirAll (const Field &in, std::vector<Field> &out){ assert(0); };
void Op (const Field &in, Field &out){
// std::cout << "Op: PVdag M "<<std::endl;
Field tmp(in.Grid());
_Mat.M(in,tmp);
_PV.Mdag(tmp,out);
out = out + shift * in;
}
void AdjOp (const Field &in, Field &out){
// std::cout << "AdjOp: Mdag PV "<<std::endl;
Field tmp(in.Grid());
_PV.M(tmp,out);
_Mat.Mdag(in,tmp);
out = out + shift * in;
}
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
void HermOp(const Field &in, Field &out){
// std::cout << "HermOp: Mdag PV PVdag M"<<std::endl;
Field tmp(in.Grid());
Op(in,tmp);
AdjOp(tmp,out);
}
};
template<class Fobj,class CComplex,int nbasis>
class MGPreconditionerSVD : public LinearFunction< Lattice<Fobj> > {
public:
using LinearFunction<Lattice<Fobj> >::operator();
typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField FineField;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
typedef LinearOperatorBase<FineField> FineOperator;
typedef LinearFunction <FineField> FineSmoother;
typedef LinearOperatorBase<CoarseVector> CoarseOperator;
typedef LinearFunction <CoarseVector> CoarseSolver;
///////////////////////////////
// SVD is M = U S Vdag
//
// Define a subset of Vc and Uc in Complex_f,c matrix
// - these are the coarsening, non-square matrices
//
// Solve a coarse approx to
//
// M psi = eta
//
// via
//
// Uc^dag U S Vdag Vc Vc^dag psi = Uc^dag eta
//
// M_coarse Vc^dag psi = M_coarse psi_c = eta_c
//
///////////////////////////////
Aggregates & _U;
Aggregates & _V;
FineOperator & _FineOperator;
FineSmoother & _PreSmoother;
FineSmoother & _PostSmoother;
CoarseOperator & _CoarseOperator;
CoarseSolver & _CoarseSolve;
int level; void Level(int lv) {level = lv; };
MGPreconditionerSVD(Aggregates &U,
Aggregates &V,
FineOperator &Fine,
FineSmoother &PreSmoother,
FineSmoother &PostSmoother,
CoarseOperator &CoarseOperator_,
CoarseSolver &CoarseSolve_)
: _U(U),
_V(V),
_FineOperator(Fine),
_PreSmoother(PreSmoother),
_PostSmoother(PostSmoother),
_CoarseOperator(CoarseOperator_),
_CoarseSolve(CoarseSolve_),
level(1) { }
virtual void operator()(const FineField &in, FineField & out)
{
GridBase *CoarseGrid = _U.CoarseGrid;
// auto CoarseGrid = _CoarseOperator.Grid();
CoarseVector Csrc(CoarseGrid);
CoarseVector Csol(CoarseGrid);
FineField vec1(in.Grid());
FineField vec2(in.Grid());
std::cout<<GridLogMessage << "Calling PreSmoother " <<std::endl;
// std::cout<<GridLogMessage << "Calling PreSmoother input residual "<<norm2(in) <<std::endl;
double t;
// Fine Smoother
// out = in;
out = Zero();
t=-usecond();
_PreSmoother(in,out);
t+=usecond();
std::cout<<GridLogMessage << "PreSmoother took "<< t/1000.0<< "ms" <<std::endl;
// Update the residual
_FineOperator.Op(out,vec1); sub(vec1, in ,vec1);
// std::cout<<GridLogMessage <<"Residual-1 now " <<norm2(vec1)<<std::endl;
// Uc^dag U S Vdag Vc Vc^dag psi = Uc^dag eta
// Fine to Coarse
t=-usecond();
_U.ProjectToSubspace (Csrc,vec1);
t+=usecond();
std::cout<<GridLogMessage << "Project to coarse took "<< t/1000.0<< "ms" <<std::endl;
// Coarse correction
t=-usecond();
Csol = Zero();
_CoarseSolve(Csrc,Csol);
//Csol=Zero();
t+=usecond();
std::cout<<GridLogMessage << "Coarse solve took "<< t/1000.0<< "ms" <<std::endl;
// Coarse to Fine
t=-usecond();
// _CoarseOperator.PromoteFromSubspace(_Aggregates,Csol,vec1);
_V.PromoteFromSubspace(Csol,vec1);
add(out,out,vec1);
t+=usecond();
std::cout<<GridLogMessage << "Promote to this level took "<< t/1000.0<< "ms" <<std::endl;
// Residual
_FineOperator.Op(out,vec1); sub(vec1 ,in , vec1);
// std::cout<<GridLogMessage <<"Residual-2 now " <<norm2(vec1)<<std::endl;
// Fine Smoother
t=-usecond();
// vec2=vec1;
vec2=Zero();
_PostSmoother(vec1,vec2);
t+=usecond();
std::cout<<GridLogMessage << "PostSmoother took "<< t/1000.0<< "ms" <<std::endl;
add( out,out,vec2);
std::cout<<GridLogMessage << "Done " <<std::endl;
}
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
const int Ls=16;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
// Construct a coarsened grid
Coordinate clatt = GridDefaultLatt();
for(int d=0;d<clatt.size();d++){
clatt[d] = clatt[d]/2;
// clatt[d] = clatt[d]/4;
}
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
std::vector<int> cseeds({5,6,7,8});
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
LatticeFermion src(FGrid); random(RNG5,src);
LatticeFermion result(FGrid); result=Zero();
LatticeFermion ref(FGrid); ref=Zero();
LatticeFermion tmp(FGrid);
LatticeFermion err(FGrid);
LatticeGaugeField Umu(UGrid);
FieldMetaData header;
std::string file("ckpoint_lat.4000");
NerscIO::readConfiguration(Umu,header,file);
RealD mass=0.01;
RealD M5=1.8;
DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
DomainWallFermionD Dpv(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0,M5);
const int nbasis = 60;
const int cb = 0 ;
NextToNearestStencilGeometry5D geom(Coarse5d);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
typedef PVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> PVdagM_t;
typedef MdagPVLinearOperator<DomainWallFermionD,LatticeFermionD> MdagPV_t;
typedef ShiftedPVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> ShiftedPVdagM_t;
PVdagM_t PVdagM(Ddwf,Dpv);
MdagPV_t MdagPV(Ddwf,Dpv);
// ShiftedPVdagM_t ShiftedPVdagM(2.0,Ddwf,Dpv); // 355
// ShiftedPVdagM_t ShiftedPVdagM(1.0,Ddwf,Dpv); // 246
// ShiftedPVdagM_t ShiftedPVdagM(0.5,Ddwf,Dpv); // 183
// ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // 145
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 134
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 127 -- NULL space via inverse iteration
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 57 -- NULL space via inverse iteration; 3 iterations
// ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // 57 , tighter inversion
// ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // nbasis 20 -- 49 iters
// ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // nbasis 20 -- 70 iters; asymmetric
// ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // 58; Loosen coarse, tighten fine
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 56 ...
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 51 ... with 24 vecs
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 31 ... with 24 vecs and 2^4 blocking
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 43 ... with 16 vecs and 2^4 blocking, sloppier
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 35 ... with 20 vecs and 2^4 blocking
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 35 ... with 20 vecs and 2^4 blocking, looser coarse
// ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 64 ... with 20 vecs, Christoph setup, and 2^4 blocking, looser coarse
ShiftedPVdagM_t ShiftedPVdagM(0.01,Ddwf,Dpv); //
// Run power method on HOA??
PowerMethod<LatticeFermion> PM;
PM(PVdagM,src);
PM(MdagPV,src);
// Warning: This routine calls PVdagM.Op, not PVdagM.HermOp
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
Subspace V(Coarse5d,FGrid,cb);
// Subspace U(Coarse5d,FGrid,cb);
// Breeds right singular vectors with call to HermOp
V.CreateSubspaceChebyshev(RNG5,PVdagM,
nbasis,
4000.0,0.003,
300);
// Breeds left singular vectors with call to HermOp
// U.CreateSubspaceChebyshev(RNG5,MdagPV,
// nbasis,
// 4000.0,0.003,
// 300);
// U.subspace=V.subspace;
// typedef Aggregation<vSpinColourVector,vTComplex,2*nbasis> CombinedSubspace;
// CombinedSubspace CombinedUV(Coarse5d,FGrid,cb);
// for(int b=0;b<nbasis;b++){
// CombinedUV.subspace[b] = V.subspace[b];
// CombinedUV.subspace[b+nbasis] = U.subspace[b];
// }
// typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,2*nbasis> LittleDiracOperator;
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
LittleDiracOperator LittleDiracOpPV(geom,FGrid,Coarse5d);
LittleDiracOpPV.CoarsenOperator(PVdagM,V,V);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"Testing coarsened operator "<<std::endl;
CoarseVector c_src (Coarse5d);
CoarseVector c_res (Coarse5d);
CoarseVector c_proj(Coarse5d);
Complex one(1.0);
c_src = one; // 1 in every element for vector 1.
// blockPromote(c_src,err,CoarseToFine.subspace);
LatticeFermion prom(FGrid);
prom=Zero();
for(int b=0;b<nbasis;b++){
prom=prom+V.subspace[b];
}
std::cout<<GridLogMessage<<"c_src "<<norm2(c_src)<<std::endl;
std::cout<<GridLogMessage<<"prom "<<norm2(prom)<<std::endl;
PVdagM.Op(prom,tmp);
blockProject(c_proj,tmp,V.subspace);
std::cout<<GridLogMessage<<" Called Big Dirac Op "<<norm2(tmp)<<std::endl;
LittleDiracOpPV.M(c_src,c_res);
std::cout<<GridLogMessage<<" Called Little Dirac Op c_src "<< norm2(c_src) << " c_res "<< norm2(c_res) <<std::endl;
std::cout<<GridLogMessage<<"Little dop : "<<norm2(c_res)<<std::endl;
std::cout<<GridLogMessage<<"Big dop in subspace : "<<norm2(c_proj)<<std::endl;
c_proj = c_proj - c_res;
std::cout<<GridLogMessage<<" ldop error: "<<norm2(c_proj)<<std::endl;
/**********
* Some solvers
**********
*/
///////////////////////////////////////
// Coarse grid solver test
///////////////////////////////////////
std::cout<<GridLogMessage<<"******************* "<<std::endl;
std::cout<<GridLogMessage<<" Coarse Grid Solve -- Level 3 "<<std::endl;
std::cout<<GridLogMessage<<"******************* "<<std::endl;
TrivialPrecon<CoarseVector> simple;
NonHermitianLinearOperator<LittleDiracOperator,CoarseVector> LinOpCoarse(LittleDiracOpPV);
// PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(1.0e-4, 100, LinOpCoarse,simple,10,10);
PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L3PGCR(1.0e-4, 10, LinOpCoarse,simple,20,20);
L3PGCR.Level(3);
c_res=Zero();
L3PGCR(c_src,c_res);
////////////////////////////////////////
// Fine grid smoother
////////////////////////////////////////
std::cout<<GridLogMessage<<"******************* "<<std::endl;
std::cout<<GridLogMessage<<" Fine Grid Smoother -- Level 2 "<<std::endl;
std::cout<<GridLogMessage<<"******************* "<<std::endl;
TrivialPrecon<LatticeFermionD> simple_fine;
// NonHermitianLinearOperator<PVdagM_t,LatticeFermionD> LinOpSmooth(PVdagM);
PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> SmootherGCR(0.01,1,ShiftedPVdagM,simple_fine,16,16);
SmootherGCR.Level(2);
LatticeFermionD f_src(FGrid);
LatticeFermionD f_res(FGrid);
f_src = one; // 1 in every element for vector 1.
f_res=Zero();
SmootherGCR(f_src,f_res);
// typedef MGPreconditionerSVD<vSpinColourVector, vTComplex,nbasis*2> TwoLevelMG;
typedef MGPreconditionerSVD<vSpinColourVector, vTComplex,nbasis> TwoLevelMG;
// TwoLevelMG TwoLevelPrecon(CombinedUV,CombinedUV,
TwoLevelMG TwoLevelPrecon(V,V,
PVdagM,
simple_fine,
SmootherGCR,
LinOpCoarse,
L3PGCR);
PrecGeneralisedConjugateResidualNonHermitian<LatticeFermion> L1PGCR(1.0e-8,1000,PVdagM,TwoLevelPrecon,16,16);
L1PGCR.Level(1);
f_res=Zero();
L1PGCR(f_src,f_res);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage << "Done "<< std::endl;
Grid_finalize();
return 0;
}
-333
View File
@@ -1,333 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_padded_cell.cc
Copyright (C) 2023
Author: Peter Boyle <paboyle@ph.ed.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
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
#include <Grid/lattice/PaddedCell.h>
#include <Grid/stencil/GeneralLocalStencil.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
#include <Grid/algorithms/iterative/BiCGSTAB.h>
using namespace std;
using namespace Grid;
template<class Fobj,class CComplex,int nbasis>
class MGPreconditioner : public LinearFunction< Lattice<Fobj> > {
public:
using LinearFunction<Lattice<Fobj> >::operator();
typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField FineField;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
typedef LinearOperatorBase<FineField> FineOperator;
typedef LinearFunction <FineField> FineSmoother;
typedef LinearOperatorBase<CoarseVector> CoarseOperator;
typedef LinearFunction <CoarseVector> CoarseSolver;
Aggregates & _Aggregates;
FineOperator & _FineOperator;
FineSmoother & _PreSmoother;
FineSmoother & _PostSmoother;
CoarseOperator & _CoarseOperator;
CoarseSolver & _CoarseSolve;
int level; void Level(int lv) {level = lv; };
MGPreconditioner(Aggregates &Agg,
FineOperator &Fine,
FineSmoother &PreSmoother,
FineSmoother &PostSmoother,
CoarseOperator &CoarseOperator_,
CoarseSolver &CoarseSolve_)
: _Aggregates(Agg),
_FineOperator(Fine),
_PreSmoother(PreSmoother),
_PostSmoother(PostSmoother),
_CoarseOperator(CoarseOperator_),
_CoarseSolve(CoarseSolve_),
level(1) { }
virtual void operator()(const FineField &in, FineField & out)
{
GridBase *CoarseGrid = _Aggregates.CoarseGrid;
// auto CoarseGrid = _CoarseOperator.Grid();
CoarseVector Csrc(CoarseGrid);
CoarseVector Csol(CoarseGrid);
FineField vec1(in.Grid());
FineField vec2(in.Grid());
std::cout<<GridLogMessage << "Calling PreSmoother " <<std::endl;
// std::cout<<GridLogMessage << "Calling PreSmoother input residual "<<norm2(in) <<std::endl;
double t;
// Fine Smoother
// out = in;
out = Zero();
t=-usecond();
_PreSmoother(in,out);
t+=usecond();
std::cout<<GridLogMessage << "PreSmoother took "<< t/1000.0<< "ms" <<std::endl;
// Update the residual
_FineOperator.Op(out,vec1); sub(vec1, in ,vec1);
// std::cout<<GridLogMessage <<"Residual-1 now " <<norm2(vec1)<<std::endl;
// Fine to Coarse
t=-usecond();
_Aggregates.ProjectToSubspace (Csrc,vec1);
t+=usecond();
std::cout<<GridLogMessage << "Project to coarse took "<< t/1000.0<< "ms" <<std::endl;
// Coarse correction
t=-usecond();
Csol = Zero();
_CoarseSolve(Csrc,Csol);
//Csol=Zero();
t+=usecond();
std::cout<<GridLogMessage << "Coarse solve took "<< t/1000.0<< "ms" <<std::endl;
// Coarse to Fine
t=-usecond();
// _CoarseOperator.PromoteFromSubspace(_Aggregates,Csol,vec1);
_Aggregates.PromoteFromSubspace(Csol,vec1);
add(out,out,vec1);
t+=usecond();
std::cout<<GridLogMessage << "Promote to this level took "<< t/1000.0<< "ms" <<std::endl;
// Residual
_FineOperator.Op(out,vec1); sub(vec1 ,in , vec1);
// std::cout<<GridLogMessage <<"Residual-2 now " <<norm2(vec1)<<std::endl;
// Fine Smoother
t=-usecond();
// vec2=vec1;
vec2=Zero();
_PostSmoother(vec1,vec2);
t+=usecond();
std::cout<<GridLogMessage << "PostSmoother took "<< t/1000.0<< "ms" <<std::endl;
add( out,out,vec2);
std::cout<<GridLogMessage << "Done " <<std::endl;
}
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
const int Ls=16;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = UGrid;
GridRedBlackCartesian * FrbGrid = UrbGrid;
// Construct a coarsened grid
Coordinate clatt = GridDefaultLatt();
for(int d=0;d<clatt.size();d++){
clatt[d] = clatt[d]/2;
//clatt[d] = clatt[d]/4;
}
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
std::vector<int> seeds4({1,2,3,4});
std::vector<int> cseeds({5,6,7,8});
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
GridParallelRNG CRNG(Coarse4d);CRNG.SeedFixedIntegers(cseeds);
Complex one(1.0);
LatticeFermion src(FGrid); src=one;
LatticeFermion result(FGrid); result=Zero();
LatticeFermion ref(FGrid); ref=Zero();
LatticeFermion tmp(FGrid);
LatticeFermion err(FGrid);
LatticeFermion precsrc(FGrid);
LatticeGaugeField Umu(UGrid);
FieldMetaData header;
std::string file("ckpoint_lat");
NerscIO::readConfiguration(Umu,header,file);
RealD csw =0.0;
RealD mass=-0.92;
WilsonCloverFermionD Dw(Umu,*UGrid,*UrbGrid,mass,csw,csw);
const int nbasis = 20;
const int cb = 0 ;
LatticeFermion prom(FGrid);
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,2*nbasis> LittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
NearestStencilGeometry4D geom(Coarse4d);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
// Warning: This routine calls Linop.Op, not LinOpo.HermOp
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
Subspace Aggregates(Coarse4d,FGrid,cb);
NonHermitianLinearOperator<WilsonCloverFermionD,LatticeFermion> LinOpDw(Dw);
ShiftedNonHermitianLinearOperator<WilsonCloverFermionD,LatticeFermion> ShiftedLinOpDw(Dw,0.01);
Aggregates.CreateSubspaceGCR(RNG4,
LinOpDw,
nbasis);
typedef Aggregation<vSpinColourVector,vTComplex,2*nbasis> CombinedSubspace;
CombinedSubspace CombinedUV(Coarse4d,UGrid,cb);
for(int b=0;b<nbasis;b++){
Gamma G5(Gamma::Algebra::Gamma5);
CombinedUV.subspace[b] = Aggregates.subspace[b];
CombinedUV.subspace[b+nbasis] = G5*Aggregates.subspace[b];
}
LittleDiracOperator LittleDiracOp(geom,FGrid,Coarse4d);
LittleDiracOp.CoarsenOperator(LinOpDw,CombinedUV);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"Testing coarsened operator "<<std::endl;
CoarseVector c_src (Coarse4d);
CoarseVector c_res (Coarse4d);
CoarseVector c_proj(Coarse4d);
std::vector<LatticeFermion> subspace(2*nbasis,FGrid);
subspace=CombinedUV.subspace;
c_src = one; // 1 in every element for vector 1.
blockPromote(c_src,err,subspace);
prom=Zero();
for(int b=0;b<2*nbasis;b++){
prom=prom+subspace[b];
}
err=err-prom;
std::cout<<GridLogMessage<<"Promoted back from subspace: err "<<norm2(err)<<std::endl;
std::cout<<GridLogMessage<<"c_src "<<norm2(c_src)<<std::endl;
std::cout<<GridLogMessage<<"prom "<<norm2(prom)<<std::endl;
LinOpDw.Op(prom,tmp);
blockProject(c_proj,tmp,subspace);
std::cout<<GridLogMessage<<" Called Big Dirac Op "<<norm2(tmp)<<std::endl;
LittleDiracOp.M(c_src,c_res);
std::cout<<GridLogMessage<<" Called Little Dirac Op c_src "<< norm2(c_src) << " c_res "<< norm2(c_res) <<std::endl;
std::cout<<GridLogMessage<<"Little dop : "<<norm2(c_res)<<std::endl;
// std::cout<<GridLogMessage<<" Little "<< c_res<<std::endl;
std::cout<<GridLogMessage<<"Big dop in subspace : "<<norm2(c_proj)<<std::endl;
// std::cout<<GridLogMessage<<" Big "<< c_proj<<std::endl;
c_proj = c_proj - c_res;
std::cout<<GridLogMessage<<" ldop error: "<<norm2(c_proj)<<std::endl;
// std::cout<<GridLogMessage<<" error "<< c_proj<<std::endl;
/**********
* Some solvers
**********
*/
// CG
{
MdagMLinearOperator<WilsonFermionD,LatticeFermion> HermOp(Dw);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
Dw.Mdag(src,precsrc);
CG(HermOp,precsrc,result);
result=Zero();
}
///////////////////////////////////////
// Coarse grid solver test
///////////////////////////////////////
std::cout<<GridLogMessage<<"******************* "<<std::endl;
std::cout<<GridLogMessage<<" Coarse Grid Solve -- Level 3 "<<std::endl;
std::cout<<GridLogMessage<<"******************* "<<std::endl;
TrivialPrecon<CoarseVector> simple;
NonHermitianLinearOperator<LittleDiracOperator,CoarseVector> LinOpCoarse(LittleDiracOp);
ShiftedNonHermitianLinearOperator<LittleDiracOperator,CoarseVector> ShiftedLinOpCoarse(LittleDiracOp,0.001);
// ShiftedNonHermitianLinearOperator<LittleDiracOperator,CoarseVector> ShiftedLinOpCoarse(LittleDiracOp,0.01);
// ShiftedNonHermitianLinearOperator<LittleDiracOperator,CoarseVector> ShiftedLinOpCoarse(LinOpCoarse,0.001);
// PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(1.0e-4, 100, LinOpCoarse,simple,10,10);
// PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(1.0e-1, 100, LinOpCoarse,simple,30,30);
PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(2.0e-1, 50, ShiftedLinOpCoarse,simple,50,50);
L2PGCR.Level(3);
c_res=Zero();
L2PGCR(c_src,c_res);
////////////////////////////////////////
// Fine grid smoother
////////////////////////////////////////
std::cout<<GridLogMessage<<"******************* "<<std::endl;
std::cout<<GridLogMessage<<" Fine Grid Smoother -- Level 2 "<<std::endl;
std::cout<<GridLogMessage<<"******************* "<<std::endl;
TrivialPrecon<LatticeFermionD> simple_fine;
PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> SmootherGCR(0.1,1,ShiftedLinOpDw,simple_fine,4,4);
SmootherGCR.Level(2);
LatticeFermionD f_src(FGrid);
LatticeFermionD f_res(FGrid);
f_src = one; // 1 in every element for vector 1.
f_res=Zero();
SmootherGCR(f_src,f_res);
typedef MGPreconditioner<vSpinColourVector, vTComplex,2*nbasis> TwoLevelMG;
TwoLevelMG TwoLevelPrecon(CombinedUV,
LinOpDw,
simple_fine,
SmootherGCR,
LinOpCoarse,
L2PGCR);
PrecGeneralisedConjugateResidualNonHermitian<LatticeFermion> L1PGCR(1.0e-8,1000,LinOpDw,TwoLevelPrecon,16,16);
L1PGCR.Level(1);
f_res=Zero();
L1PGCR(f_src,f_res);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage << "Done "<< std::endl;
Grid_finalize();
return 0;
}
@@ -1,326 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_padded_cell.cc
Copyright (C) 2023
Author: Peter Boyle <paboyle@ph.ed.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
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
#include <Grid/lattice/PaddedCell.h>
#include <Grid/stencil/GeneralLocalStencil.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
#include <Grid/algorithms/iterative/BiCGSTAB.h>
using namespace std;
using namespace Grid;
template<class Fobj,class CComplex,int nbasis>
class MGPreconditioner : public LinearFunction< Lattice<Fobj> > {
public:
using LinearFunction<Lattice<Fobj> >::operator();
typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField FineField;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
typedef LinearOperatorBase<FineField> FineOperator;
typedef LinearFunction <FineField> FineSmoother;
typedef LinearOperatorBase<CoarseVector> CoarseOperator;
typedef LinearFunction <CoarseVector> CoarseSolver;
Aggregates & _Aggregates;
FineOperator & _FineOperator;
FineSmoother & _PreSmoother;
FineSmoother & _PostSmoother;
CoarseOperator & _CoarseOperator;
CoarseSolver & _CoarseSolve;
int level; void Level(int lv) {level = lv; };
MGPreconditioner(Aggregates &Agg,
FineOperator &Fine,
FineSmoother &PreSmoother,
FineSmoother &PostSmoother,
CoarseOperator &CoarseOperator_,
CoarseSolver &CoarseSolve_)
: _Aggregates(Agg),
_FineOperator(Fine),
_PreSmoother(PreSmoother),
_PostSmoother(PostSmoother),
_CoarseOperator(CoarseOperator_),
_CoarseSolve(CoarseSolve_),
level(1) { }
virtual void operator()(const FineField &in, FineField & out)
{
GridBase *CoarseGrid = _Aggregates.CoarseGrid;
// auto CoarseGrid = _CoarseOperator.Grid();
CoarseVector Csrc(CoarseGrid);
CoarseVector Csol(CoarseGrid);
FineField vec1(in.Grid());
FineField vec2(in.Grid());
std::cout<<GridLogMessage << "Calling PreSmoother " <<std::endl;
// std::cout<<GridLogMessage << "Calling PreSmoother input residual "<<norm2(in) <<std::endl;
double t;
// Fine Smoother
// out = in;
out = Zero();
t=-usecond();
_PreSmoother(in,out);
t+=usecond();
std::cout<<GridLogMessage << "PreSmoother took "<< t/1000.0<< "ms" <<std::endl;
// Update the residual
_FineOperator.Op(out,vec1); sub(vec1, in ,vec1);
// std::cout<<GridLogMessage <<"Residual-1 now " <<norm2(vec1)<<std::endl;
// Fine to Coarse
t=-usecond();
_Aggregates.ProjectToSubspace (Csrc,vec1);
t+=usecond();
std::cout<<GridLogMessage << "Project to coarse took "<< t/1000.0<< "ms" <<std::endl;
// Coarse correction
t=-usecond();
Csol = Zero();
_CoarseSolve(Csrc,Csol);
//Csol=Zero();
t+=usecond();
std::cout<<GridLogMessage << "Coarse solve took "<< t/1000.0<< "ms" <<std::endl;
// Coarse to Fine
t=-usecond();
// _CoarseOperator.PromoteFromSubspace(_Aggregates,Csol,vec1);
_Aggregates.PromoteFromSubspace(Csol,vec1);
add(out,out,vec1);
t+=usecond();
std::cout<<GridLogMessage << "Promote to this level took "<< t/1000.0<< "ms" <<std::endl;
// Residual
_FineOperator.Op(out,vec1); sub(vec1 ,in , vec1);
// std::cout<<GridLogMessage <<"Residual-2 now " <<norm2(vec1)<<std::endl;
// Fine Smoother
t=-usecond();
// vec2=vec1;
vec2=Zero();
_PostSmoother(vec1,vec2);
t+=usecond();
std::cout<<GridLogMessage << "PostSmoother took "<< t/1000.0<< "ms" <<std::endl;
add( out,out,vec2);
std::cout<<GridLogMessage << "Done " <<std::endl;
}
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
const int Ls=16;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = UGrid;
GridRedBlackCartesian * FrbGrid = UrbGrid;
// Construct a coarsened grid
Coordinate clatt = GridDefaultLatt();
for(int d=0;d<clatt.size();d++){
clatt[d] = clatt[d]/2;
// clatt[d] = clatt[d]/4;
}
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
std::vector<int> seeds4({1,2,3,4});
std::vector<int> cseeds({5,6,7,8});
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
GridParallelRNG CRNG(Coarse4d);CRNG.SeedFixedIntegers(cseeds);
Complex one(1.0);
LatticeFermion src(FGrid); src=one;
LatticeFermion result(FGrid); result=Zero();
LatticeFermion ref(FGrid); ref=Zero();
LatticeFermion tmp(FGrid);
LatticeFermion err(FGrid);
LatticeFermion precsrc(FGrid);
LatticeGaugeField Umu(UGrid);
FieldMetaData header;
std::string file("ckpoint_lat");
NerscIO::readConfiguration(Umu,header,file);
RealD csw =0.0;
RealD mass=-0.92;
WilsonCloverFermionD Dw(Umu,*UGrid,*UrbGrid,mass,csw,csw);
const int nbasis = 40;
const int cb = 0 ;
LatticeFermion prom(FGrid);
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
NearestStencilGeometry4D geom(Coarse4d);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
// Warning: This routine calls Linop.Op, not LinOpo.HermOp
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
Subspace Aggregates(Coarse4d,FGrid,cb);
NonHermitianLinearOperator<WilsonCloverFermionD,LatticeFermion> LinOpDw(Dw);
ShiftedNonHermitianLinearOperator<WilsonCloverFermionD,LatticeFermion> ShiftedLinOpDw(Dw,0.01);
Aggregates.CreateSubspaceGCR(RNG4,
LinOpDw,
nbasis);
LittleDiracOperator LittleDiracOp(geom,FGrid,Coarse4d);
LittleDiracOp.CoarsenOperator(LinOpDw,Aggregates);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"Testing coarsened operator "<<std::endl;
CoarseVector c_src (Coarse4d);
CoarseVector c_res (Coarse4d);
CoarseVector c_proj(Coarse4d);
std::vector<LatticeFermion> subspace(nbasis,FGrid);
subspace=Aggregates.subspace;
c_src = one; // 1 in every element for vector 1.
blockPromote(c_src,err,subspace);
prom=Zero();
for(int b=0;b<nbasis;b++){
prom=prom+subspace[b];
}
err=err-prom;
std::cout<<GridLogMessage<<"Promoted back from subspace: err "<<norm2(err)<<std::endl;
std::cout<<GridLogMessage<<"c_src "<<norm2(c_src)<<std::endl;
std::cout<<GridLogMessage<<"prom "<<norm2(prom)<<std::endl;
LinOpDw.Op(prom,tmp);
blockProject(c_proj,tmp,subspace);
std::cout<<GridLogMessage<<" Called Big Dirac Op "<<norm2(tmp)<<std::endl;
LittleDiracOp.M(c_src,c_res);
std::cout<<GridLogMessage<<" Called Little Dirac Op c_src "<< norm2(c_src) << " c_res "<< norm2(c_res) <<std::endl;
std::cout<<GridLogMessage<<"Little dop : "<<norm2(c_res)<<std::endl;
// std::cout<<GridLogMessage<<" Little "<< c_res<<std::endl;
std::cout<<GridLogMessage<<"Big dop in subspace : "<<norm2(c_proj)<<std::endl;
// std::cout<<GridLogMessage<<" Big "<< c_proj<<std::endl;
c_proj = c_proj - c_res;
std::cout<<GridLogMessage<<" ldop error: "<<norm2(c_proj)<<std::endl;
// std::cout<<GridLogMessage<<" error "<< c_proj<<std::endl;
/**********
* Some solvers
**********
*/
// CG
{
MdagMLinearOperator<WilsonFermionD,LatticeFermion> HermOp(Dw);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
Dw.Mdag(src,precsrc);
CG(HermOp,precsrc,result);
result=Zero();
}
///////////////////////////////////////
// Coarse grid solver test
///////////////////////////////////////
std::cout<<GridLogMessage<<"******************* "<<std::endl;
std::cout<<GridLogMessage<<" Coarse Grid Solve -- Level 3 "<<std::endl;
std::cout<<GridLogMessage<<"******************* "<<std::endl;
TrivialPrecon<CoarseVector> simple;
NonHermitianLinearOperator<LittleDiracOperator,CoarseVector> LinOpCoarse(LittleDiracOp);
ShiftedNonHermitianLinearOperator<LittleDiracOperator,CoarseVector> ShiftedLinOpCoarse(LittleDiracOp,0.001);
// ShiftedNonHermitianLinearOperator<LittleDiracOperator,CoarseVector> ShiftedLinOpCoarse(LittleDiracOp,0.01);
// ShiftedNonHermitianLinearOperator<LittleDiracOperator,CoarseVector> ShiftedLinOpCoarse(LinOpCoarse,0.001);
// PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(1.0e-4, 100, LinOpCoarse,simple,10,10);
PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(1.0e-1, 100, LinOpCoarse,simple,30,30);
// PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(2.0e-1, 50, ShiftedLinOpCoarse,simple,50,50);
L2PGCR.Level(3);
c_res=Zero();
L2PGCR(c_src,c_res);
////////////////////////////////////////
// Fine grid smoother
////////////////////////////////////////
std::cout<<GridLogMessage<<"******************* "<<std::endl;
std::cout<<GridLogMessage<<" Fine Grid Smoother -- Level 2 "<<std::endl;
std::cout<<GridLogMessage<<"******************* "<<std::endl;
TrivialPrecon<LatticeFermionD> simple_fine;
PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> SmootherGCR(0.1,1,ShiftedLinOpDw,simple_fine,6,6);
SmootherGCR.Level(2);
LatticeFermionD f_src(FGrid);
LatticeFermionD f_res(FGrid);
f_src = one; // 1 in every element for vector 1.
f_res=Zero();
SmootherGCR(f_src,f_res);
typedef MGPreconditioner<vSpinColourVector, vTComplex,nbasis> TwoLevelMG;
TwoLevelMG TwoLevelPrecon(Aggregates,
LinOpDw,
simple_fine,
SmootherGCR,
LinOpCoarse,
L2PGCR);
PrecGeneralisedConjugateResidualNonHermitian<LatticeFermion> L1PGCR(1.0e-8,1000,LinOpDw,TwoLevelPrecon,16,16);
L1PGCR.Level(1);
f_res=Zero();
L1PGCR(f_src,f_res);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage << "Done "<< std::endl;
Grid_finalize();
return 0;
}
@@ -1,320 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_padded_cell.cc
Copyright (C) 2023
Author: Peter Boyle <paboyle@ph.ed.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
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
#include <Grid/lattice/PaddedCell.h>
#include <Grid/stencil/GeneralLocalStencil.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
#include <Grid/algorithms/iterative/BiCGSTAB.h>
using namespace std;
using namespace Grid;
template<class Fobj,class CComplex,int nbasis>
class MGPreconditioner : public LinearFunction< Lattice<Fobj> > {
public:
using LinearFunction<Lattice<Fobj> >::operator();
typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField FineField;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
typedef LinearOperatorBase<FineField> FineOperator;
typedef LinearFunction <FineField> FineSmoother;
typedef LinearOperatorBase<CoarseVector> CoarseOperator;
typedef LinearFunction <CoarseVector> CoarseSolver;
Aggregates & _Aggregates;
FineOperator & _FineOperator;
FineSmoother & _PreSmoother;
FineSmoother & _PostSmoother;
CoarseOperator & _CoarseOperator;
CoarseSolver & _CoarseSolve;
int level; void Level(int lv) {level = lv; };
MGPreconditioner(Aggregates &Agg,
FineOperator &Fine,
FineSmoother &PreSmoother,
FineSmoother &PostSmoother,
CoarseOperator &CoarseOperator_,
CoarseSolver &CoarseSolve_)
: _Aggregates(Agg),
_FineOperator(Fine),
_PreSmoother(PreSmoother),
_PostSmoother(PostSmoother),
_CoarseOperator(CoarseOperator_),
_CoarseSolve(CoarseSolve_),
level(1) { }
virtual void operator()(const FineField &in, FineField & out)
{
GridBase *CoarseGrid = _Aggregates.CoarseGrid;
// auto CoarseGrid = _CoarseOperator.Grid();
CoarseVector Csrc(CoarseGrid);
CoarseVector Csol(CoarseGrid);
FineField vec1(in.Grid());
FineField vec2(in.Grid());
std::cout<<GridLogMessage << "Calling PreSmoother " <<std::endl;
// std::cout<<GridLogMessage << "Calling PreSmoother input residual "<<norm2(in) <<std::endl;
double t;
// Fine Smoother
// out = in;
out = Zero();
t=-usecond();
_PreSmoother(in,out);
t+=usecond();
std::cout<<GridLogMessage << "PreSmoother took "<< t/1000.0<< "ms" <<std::endl;
// Update the residual
_FineOperator.Op(out,vec1); sub(vec1, in ,vec1);
// std::cout<<GridLogMessage <<"Residual-1 now " <<norm2(vec1)<<std::endl;
// Fine to Coarse
t=-usecond();
_Aggregates.ProjectToSubspace (Csrc,vec1);
t+=usecond();
std::cout<<GridLogMessage << "Project to coarse took "<< t/1000.0<< "ms" <<std::endl;
// Coarse correction
t=-usecond();
Csol = Zero();
_CoarseSolve(Csrc,Csol);
//Csol=Zero();
t+=usecond();
std::cout<<GridLogMessage << "Coarse solve took "<< t/1000.0<< "ms" <<std::endl;
// Coarse to Fine
t=-usecond();
// _CoarseOperator.PromoteFromSubspace(_Aggregates,Csol,vec1);
_Aggregates.PromoteFromSubspace(Csol,vec1);
add(out,out,vec1);
t+=usecond();
std::cout<<GridLogMessage << "Promote to this level took "<< t/1000.0<< "ms" <<std::endl;
// Residual
_FineOperator.Op(out,vec1); sub(vec1 ,in , vec1);
// std::cout<<GridLogMessage <<"Residual-2 now " <<norm2(vec1)<<std::endl;
// Fine Smoother
t=-usecond();
// vec2=vec1;
vec2=Zero();
_PostSmoother(vec1,vec2);
t+=usecond();
std::cout<<GridLogMessage << "PostSmoother took "<< t/1000.0<< "ms" <<std::endl;
add( out,out,vec2);
std::cout<<GridLogMessage << "Done " <<std::endl;
}
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
const int Ls=16;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = UGrid;
GridRedBlackCartesian * FrbGrid = UrbGrid;
// Construct a coarsened grid
Coordinate clatt = GridDefaultLatt();
for(int d=0;d<clatt.size();d++){
clatt[d] = clatt[d]/2;
// clatt[d] = clatt[d]/4;
}
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
std::vector<int> seeds4({1,2,3,4});
std::vector<int> cseeds({5,6,7,8});
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
GridParallelRNG CRNG(Coarse4d);CRNG.SeedFixedIntegers(cseeds);
LatticeFermion src(FGrid); random(RNG4,src);
LatticeFermion result(FGrid); result=Zero();
LatticeFermion ref(FGrid); ref=Zero();
LatticeFermion tmp(FGrid);
LatticeFermion err(FGrid);
LatticeGaugeField Umu(UGrid);
FieldMetaData header;
std::string file("ckpoint_lat");
NerscIO::readConfiguration(Umu,header,file);
RealD csw =0.0;
RealD mass=-0.92;
WilsonCloverFermionD Dw(Umu,*UGrid,*UrbGrid,mass,csw,csw);
const int nbasis = 20;
const int cb = 0 ;
LatticeFermion prom(FGrid);
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,2*nbasis> LittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
NearestStencilGeometry4D geom(Coarse4d);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
// Warning: This routine calls Linop.Op, not LinOpo.HermOp
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
Subspace Aggregates(Coarse4d,FGrid,cb);
MdagMLinearOperator<WilsonCloverFermionD,LatticeFermion> MdagMOpDw(Dw);
NonHermitianLinearOperator<WilsonCloverFermionD,LatticeFermion> LinOpDw(Dw);
ShiftedNonHermitianLinearOperator<WilsonCloverFermionD,LatticeFermion> ShiftedLinOpDw(Dw,0.5);
// Aggregates.CreateSubspaceGCR(RNG4,
// LinOpDw,
// nbasis);
Aggregates.CreateSubspace(RNG4,MdagMOpDw,nbasis);
typedef Aggregation<vSpinColourVector,vTComplex,2*nbasis> CombinedSubspace;
CombinedSubspace CombinedUV(Coarse4d,UGrid,cb);
for(int b=0;b<nbasis;b++){
Gamma G5(Gamma::Algebra::Gamma5);
CombinedUV.subspace[b] = Aggregates.subspace[b];
CombinedUV.subspace[b+nbasis] = G5*Aggregates.subspace[b];
}
LittleDiracOperator LittleDiracOp(geom,FGrid,Coarse4d);
LittleDiracOp.CoarsenOperator(LinOpDw,CombinedUV);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"Testing coarsened operator "<<std::endl;
CoarseVector c_src (Coarse4d);
CoarseVector c_res (Coarse4d);
CoarseVector c_proj(Coarse4d);
std::vector<LatticeFermion> subspace(2*nbasis,FGrid);
subspace=CombinedUV.subspace;
Complex one(1.0);
c_src = one; // 1 in every element for vector 1.
blockPromote(c_src,err,subspace);
prom=Zero();
for(int b=0;b<2*nbasis;b++){
prom=prom+subspace[b];
}
err=err-prom;
std::cout<<GridLogMessage<<"Promoted back from subspace: err "<<norm2(err)<<std::endl;
std::cout<<GridLogMessage<<"c_src "<<norm2(c_src)<<std::endl;
std::cout<<GridLogMessage<<"prom "<<norm2(prom)<<std::endl;
LinOpDw.Op(prom,tmp);
blockProject(c_proj,tmp,subspace);
std::cout<<GridLogMessage<<" Called Big Dirac Op "<<norm2(tmp)<<std::endl;
LittleDiracOp.M(c_src,c_res);
std::cout<<GridLogMessage<<" Called Little Dirac Op c_src "<< norm2(c_src) << " c_res "<< norm2(c_res) <<std::endl;
std::cout<<GridLogMessage<<"Little dop : "<<norm2(c_res)<<std::endl;
// std::cout<<GridLogMessage<<" Little "<< c_res<<std::endl;
std::cout<<GridLogMessage<<"Big dop in subspace : "<<norm2(c_proj)<<std::endl;
// std::cout<<GridLogMessage<<" Big "<< c_proj<<std::endl;
c_proj = c_proj - c_res;
std::cout<<GridLogMessage<<" ldop error: "<<norm2(c_proj)<<std::endl;
// std::cout<<GridLogMessage<<" error "<< c_proj<<std::endl;
/**********
* Some solvers
**********
*/
///////////////////////////////////////
// Coarse grid solver test
///////////////////////////////////////
std::cout<<GridLogMessage<<"******************* "<<std::endl;
std::cout<<GridLogMessage<<" Coarse Grid Solve -- Level 3 "<<std::endl;
std::cout<<GridLogMessage<<"******************* "<<std::endl;
TrivialPrecon<CoarseVector> simple;
NonHermitianLinearOperator<LittleDiracOperator,CoarseVector> LinOpCoarse(LittleDiracOp);
// PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(1.0e-4, 100, LinOpCoarse,simple,10,10);
PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(1.0e-2, 100, LinOpCoarse,simple,30,30);
L2PGCR.Level(3);
c_res=Zero();
L2PGCR(c_src,c_res);
////////////////////////////////////////
// Fine grid smoother
////////////////////////////////////////
std::cout<<GridLogMessage<<"******************* "<<std::endl;
std::cout<<GridLogMessage<<" Fine Grid Smoother -- Level 2 "<<std::endl;
std::cout<<GridLogMessage<<"******************* "<<std::endl;
TrivialPrecon<LatticeFermionD> simple_fine;
PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> SmootherGCR(0.01,1,ShiftedLinOpDw,simple_fine,4,4);
SmootherGCR.Level(2);
LatticeFermionD f_src(FGrid);
LatticeFermionD f_res(FGrid);
f_src = one; // 1 in every element for vector 1.
f_res=Zero();
SmootherGCR(f_src,f_res);
typedef MGPreconditioner<vSpinColourVector, vTComplex,2*nbasis> TwoLevelMG;
TwoLevelMG TwoLevelPrecon(CombinedUV,
LinOpDw,
simple_fine,
SmootherGCR,
LinOpCoarse,
L2PGCR);
PrecGeneralisedConjugateResidualNonHermitian<LatticeFermion> L1PGCR(1.0e-8,1000,LinOpDw,TwoLevelPrecon,32,32);
L1PGCR.Level(1);
f_res=Zero();
L1PGCR(f_src,f_res);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage << "Done "<< std::endl;
Grid_finalize();
return 0;
}
@@ -1,312 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_padded_cell.cc
Copyright (C) 2023
Author: Peter Boyle <paboyle@ph.ed.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
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
#include <Grid/lattice/PaddedCell.h>
#include <Grid/stencil/GeneralLocalStencil.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
#include <Grid/algorithms/iterative/BiCGSTAB.h>
using namespace std;
using namespace Grid;
template<class Fobj,class CComplex,int nbasis>
class MGPreconditioner : public LinearFunction< Lattice<Fobj> > {
public:
using LinearFunction<Lattice<Fobj> >::operator();
typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField FineField;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
typedef LinearOperatorBase<FineField> FineOperator;
typedef LinearFunction <FineField> FineSmoother;
typedef LinearOperatorBase<CoarseVector> CoarseOperator;
typedef LinearFunction <CoarseVector> CoarseSolver;
Aggregates & _Aggregates;
FineOperator & _FineOperator;
FineSmoother & _PreSmoother;
FineSmoother & _PostSmoother;
CoarseOperator & _CoarseOperator;
CoarseSolver & _CoarseSolve;
int level; void Level(int lv) {level = lv; };
MGPreconditioner(Aggregates &Agg,
FineOperator &Fine,
FineSmoother &PreSmoother,
FineSmoother &PostSmoother,
CoarseOperator &CoarseOperator_,
CoarseSolver &CoarseSolve_)
: _Aggregates(Agg),
_FineOperator(Fine),
_PreSmoother(PreSmoother),
_PostSmoother(PostSmoother),
_CoarseOperator(CoarseOperator_),
_CoarseSolve(CoarseSolve_),
level(1) { }
virtual void operator()(const FineField &in, FineField & out)
{
GridBase *CoarseGrid = _Aggregates.CoarseGrid;
// auto CoarseGrid = _CoarseOperator.Grid();
CoarseVector Csrc(CoarseGrid);
CoarseVector Csol(CoarseGrid);
FineField vec1(in.Grid());
FineField vec2(in.Grid());
std::cout<<GridLogMessage << "Calling PreSmoother " <<std::endl;
// std::cout<<GridLogMessage << "Calling PreSmoother input residual "<<norm2(in) <<std::endl;
double t;
// Fine Smoother
// out = in;
out = Zero();
t=-usecond();
_PreSmoother(in,out);
t+=usecond();
std::cout<<GridLogMessage << "PreSmoother took "<< t/1000.0<< "ms" <<std::endl;
// Update the residual
_FineOperator.Op(out,vec1); sub(vec1, in ,vec1);
// std::cout<<GridLogMessage <<"Residual-1 now " <<norm2(vec1)<<std::endl;
// Fine to Coarse
t=-usecond();
_Aggregates.ProjectToSubspace (Csrc,vec1);
t+=usecond();
std::cout<<GridLogMessage << "Project to coarse took "<< t/1000.0<< "ms" <<std::endl;
// Coarse correction
t=-usecond();
Csol = Zero();
_CoarseSolve(Csrc,Csol);
//Csol=Zero();
t+=usecond();
std::cout<<GridLogMessage << "Coarse solve took "<< t/1000.0<< "ms" <<std::endl;
// Coarse to Fine
t=-usecond();
// _CoarseOperator.PromoteFromSubspace(_Aggregates,Csol,vec1);
_Aggregates.PromoteFromSubspace(Csol,vec1);
add(out,out,vec1);
t+=usecond();
std::cout<<GridLogMessage << "Promote to this level took "<< t/1000.0<< "ms" <<std::endl;
// Residual
_FineOperator.Op(out,vec1); sub(vec1 ,in , vec1);
// std::cout<<GridLogMessage <<"Residual-2 now " <<norm2(vec1)<<std::endl;
// Fine Smoother
t=-usecond();
// vec2=vec1;
vec2=Zero();
_PostSmoother(vec1,vec2);
t+=usecond();
std::cout<<GridLogMessage << "PostSmoother took "<< t/1000.0<< "ms" <<std::endl;
add( out,out,vec2);
std::cout<<GridLogMessage << "Done " <<std::endl;
}
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
const int Ls=16;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = UGrid;
GridRedBlackCartesian * FrbGrid = UrbGrid;
// Construct a coarsened grid
Coordinate clatt = GridDefaultLatt();
for(int d=0;d<clatt.size();d++){
clatt[d] = clatt[d]/2;
// clatt[d] = clatt[d]/4;
}
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
std::vector<int> seeds4({1,2,3,4});
std::vector<int> cseeds({5,6,7,8});
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
GridParallelRNG CRNG(Coarse4d);CRNG.SeedFixedIntegers(cseeds);
LatticeFermion src(FGrid); random(RNG4,src);
LatticeFermion result(FGrid); result=Zero();
LatticeFermion ref(FGrid); ref=Zero();
LatticeFermion tmp(FGrid);
LatticeFermion err(FGrid);
LatticeGaugeField Umu(UGrid);
FieldMetaData header;
std::string file("ckpoint_lat");
NerscIO::readConfiguration(Umu,header,file);
RealD csw =0.0;
RealD mass=-0.92;
WilsonCloverFermionD Dw(Umu,*UGrid,*UrbGrid,mass,csw,csw);
const int nbasis = 40;
const int cb = 0 ;
LatticeFermion prom(FGrid);
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
NearestStencilGeometry4D geom(Coarse4d);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
// Warning: This routine calls Linop.Op, not LinOpo.HermOp
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
Subspace Aggregates(Coarse4d,FGrid,cb);
MdagMLinearOperator<WilsonCloverFermionD,LatticeFermion> MdagMOpDw(Dw);
NonHermitianLinearOperator<WilsonCloverFermionD,LatticeFermion> LinOpDw(Dw);
ShiftedNonHermitianLinearOperator<WilsonCloverFermionD,LatticeFermion> ShiftedLinOpDw(Dw,0.5);
// Aggregates.CreateSubspaceGCR(RNG4,
// LinOpDw,
// nbasis);
Aggregates.CreateSubspace(RNG4,MdagMOpDw,nbasis);
LittleDiracOperator LittleDiracOp(geom,FGrid,Coarse4d);
LittleDiracOp.CoarsenOperator(LinOpDw,Aggregates);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"Testing coarsened operator "<<std::endl;
CoarseVector c_src (Coarse4d);
CoarseVector c_res (Coarse4d);
CoarseVector c_proj(Coarse4d);
std::vector<LatticeFermion> subspace(nbasis,FGrid);
subspace=Aggregates.subspace;
Complex one(1.0);
c_src = one; // 1 in every element for vector 1.
blockPromote(c_src,err,subspace);
prom=Zero();
for(int b=0;b<nbasis;b++){
prom=prom+subspace[b];
}
err=err-prom;
std::cout<<GridLogMessage<<"Promoted back from subspace: err "<<norm2(err)<<std::endl;
std::cout<<GridLogMessage<<"c_src "<<norm2(c_src)<<std::endl;
std::cout<<GridLogMessage<<"prom "<<norm2(prom)<<std::endl;
LinOpDw.Op(prom,tmp);
blockProject(c_proj,tmp,subspace);
std::cout<<GridLogMessage<<" Called Big Dirac Op "<<norm2(tmp)<<std::endl;
LittleDiracOp.M(c_src,c_res);
std::cout<<GridLogMessage<<" Called Little Dirac Op c_src "<< norm2(c_src) << " c_res "<< norm2(c_res) <<std::endl;
std::cout<<GridLogMessage<<"Little dop : "<<norm2(c_res)<<std::endl;
// std::cout<<GridLogMessage<<" Little "<< c_res<<std::endl;
std::cout<<GridLogMessage<<"Big dop in subspace : "<<norm2(c_proj)<<std::endl;
// std::cout<<GridLogMessage<<" Big "<< c_proj<<std::endl;
c_proj = c_proj - c_res;
std::cout<<GridLogMessage<<" ldop error: "<<norm2(c_proj)<<std::endl;
// std::cout<<GridLogMessage<<" error "<< c_proj<<std::endl;
/**********
* Some solvers
**********
*/
///////////////////////////////////////
// Coarse grid solver test
///////////////////////////////////////
std::cout<<GridLogMessage<<"******************* "<<std::endl;
std::cout<<GridLogMessage<<" Coarse Grid Solve -- Level 3 "<<std::endl;
std::cout<<GridLogMessage<<"******************* "<<std::endl;
TrivialPrecon<CoarseVector> simple;
NonHermitianLinearOperator<LittleDiracOperator,CoarseVector> LinOpCoarse(LittleDiracOp);
// PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(1.0e-4, 100, LinOpCoarse,simple,10,10);
PrecGeneralisedConjugateResidualNonHermitian<CoarseVector> L2PGCR(1.0e-2, 100, LinOpCoarse,simple,30,30);
L2PGCR.Level(3);
c_res=Zero();
L2PGCR(c_src,c_res);
////////////////////////////////////////
// Fine grid smoother
////////////////////////////////////////
std::cout<<GridLogMessage<<"******************* "<<std::endl;
std::cout<<GridLogMessage<<" Fine Grid Smoother -- Level 2 "<<std::endl;
std::cout<<GridLogMessage<<"******************* "<<std::endl;
TrivialPrecon<LatticeFermionD> simple_fine;
PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> SmootherGCR(0.01,1,ShiftedLinOpDw,simple_fine,6,6);
SmootherGCR.Level(2);
LatticeFermionD f_src(FGrid);
LatticeFermionD f_res(FGrid);
f_src = one; // 1 in every element for vector 1.
f_res=Zero();
SmootherGCR(f_src,f_res);
typedef MGPreconditioner<vSpinColourVector, vTComplex,nbasis> TwoLevelMG;
TwoLevelMG TwoLevelPrecon(Aggregates,
LinOpDw,
simple_fine,
SmootherGCR,
LinOpCoarse,
L2PGCR);
PrecGeneralisedConjugateResidualNonHermitian<LatticeFermion> L1PGCR(1.0e-8,1000,LinOpDw,TwoLevelPrecon,32,32);
L1PGCR.Level(1);
f_res=Zero();
L1PGCR(f_src,f_res);
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
std::cout<<GridLogMessage<<std::endl;
std::cout<<GridLogMessage << "Done "<< std::endl;
Grid_finalize();
return 0;
}
+1 -1
View File
@@ -490,7 +490,7 @@ public:
}
}
assert(s==nshift);
GRID_ASSERT(s==nshift);
coalescedWrite(gStaple_v[ss],stencil_ss);
}
);
@@ -1,16 +0,0 @@
#include <Grid/Grid.h>
#pragma once
#ifndef ENABLE_FERMION_INSTANTIATIONS
#include <iostream>
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
-7
View File
@@ -62,7 +62,6 @@ void ForceTest(Action<LatticeGaugeField> &action,ConfigurationBase<LatticeGaugeF
Gimpl::generate_momenta(P,sRNG,RNG4);
// Filter.applyFilter(P);
std::cout << GridLogMessage << "Initial momenta " << norm2(P) << std::endl;
action.refresh(smU,sRNG,RNG4);
@@ -71,8 +70,6 @@ void ForceTest(Action<LatticeGaugeField> &action,ConfigurationBase<LatticeGaugeF
std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
RealD S1 = action.S(smU);
std::cout << GridLogMessage << "Initial action " << S1 << std::endl;
Gimpl::update_field(P,U,eps);
smU.set_Field(U);
@@ -83,7 +80,6 @@ void ForceTest(Action<LatticeGaugeField> &action,ConfigurationBase<LatticeGaugeF
action.deriv(smU,UdSdU);
UdSdU = Ta(UdSdU);
// Filter.applyFilter(UdSdU);
std::cout << GridLogMessage << "Derivative " << norm2(UdSdU) << std::endl;
DumpSliceNorm("Force",UdSdU,Nd-1);
@@ -95,7 +91,6 @@ void ForceTest(Action<LatticeGaugeField> &action,ConfigurationBase<LatticeGaugeF
std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
RealD S2 = action.S(smU);
std::cout << GridLogMessage << "Final action " << S1 << std::endl;
// Use the derivative
LatticeComplex dS(UGrid); dS = Zero();
@@ -150,8 +145,6 @@ int main (int argc, char ** argv)
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds);
SU<Nc>::HotConfiguration(RNG4,U);
#endif
std::cout << GridLogMessage << "Initial plaquette: " << WilsonLoops<PeriodicGimplR>::avgPlaquette(U) << std::endl;
WilsonGaugeActionR PlaqAction(6.0);
-171
View File
@@ -1,171 +0,0 @@
# 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 <target> # 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 <N>` / `file +N` / `file -N` | Jump to file by index or relative |
| `slice <dim> <N>` / `+N` / `-N` | Set or shift slice coordinate in dimension dim |
| `spin <deg>` | Continuous azimuth rotation at deg/tick (100ms tick); `spin 0` stops |
| `azimuth <deg>` | Single azimuth rotation step |
| `elevation <deg>` | Single elevation rotation step |
| `zoom <factor>` | Camera dolly (>1 = in) |
| `iso <value>` | 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.34.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 <file> 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 \
<files...> 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_file> <output.mp4>
```
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)
-891
View File
@@ -1,891 +0,0 @@
// 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 <dim> <N> set Slice[dim] = N (0-based, wraps to lattice size)
// slice <dim> +<N> increment Slice[dim] by N
// slice <dim> -<N> decrement Slice[dim] by N
// zoom <factor> camera Dolly by factor (>1 = in, <1 = out)
// iso <value> set isosurface threshold to <value> x RMS
// file <index> jump to file by absolute index
// file +<N> advance N files
// file -<N> 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 <vtkActor.h>
#include <vtkCamera.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkOutlineFilter.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkStripper.h>
#include <vtkImageData.h>
#include <vtkCallbackCommand.h>
#include <vtkTextActor.h>
#include <vtkTextProperty.h>
#include <vtkProperty2D.h>
#include <vtkWindowToImageFilter.h>
#define MPEG
#ifdef MPEG
#include <vtkFFMPEGWriter.h>
#endif
#include <array>
#include <string>
#include <vector>
#include <queue>
#include <mutex>
#include <thread>
#include <atomic>
#include <sstream>
#include <iostream>
#include <fstream>
#include <cmath>
#include <cstdlib>
#include <sys/stat.h>
#include <fcntl.h>
#include <unistd.h>
#include <Grid/Grid.h>
#define USE_FLYING_EDGES
#ifdef USE_FLYING_EDGES
#include <vtkFlyingEdges3D.h>
typedef vtkFlyingEdges3D isosurface;
#else
#include <vtkMarchingCubes.h>
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<std::string> g_cmdQueue;
static std::mutex g_cmdMutex;
static std::atomic<bool> 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<int16_t> 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<int16_t>& 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 <class T>
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<std::mutex> 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<std::string> 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<std::string> 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<vtkRenderWindowInteractor*>(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<int*>(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<vtkRenderWindowInteractor*>(caller)->GetRenderWindow()->Render();
++TimerCount;
if (TimerCount >= maxCount && timerId > -1)
dynamic_cast<vtkRenderWindowInteractor*>(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 <dim> <N|+N|-N> ────────────────────────────────────────────
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 <factor> ────────────────────────────────────────────────────
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 <degrees> ────────────────────────────────────────────────
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 <degrees> ──────────────────────────────────────────────
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 <degrees_per_tick> ──────────────────────────────────────────
// 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: <text> / caption claude: <text> / 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 <filename> / record stop ────────────────────────────
else if (verb == "record") {
#ifdef MPEG
std::string sub;
if (!(iss >> sub)) { std::cout << "[cmd] record: expected start <file> 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 <value> ──────────────────────────────────────────────────────
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 <index|+N|-N> ───────────────────────────────────────────────
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<int*>(callData));
if (tid != pollTimerId) return;
}
std::vector<std::string> pending;
{
std::lock_guard<std::mutex> 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<std::string> 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<vtkNamedColors> colors;
std::array<unsigned char,4> posColor{{240,184,160,255}}; colors->SetColor("posColor", posColor.data());
std::array<unsigned char,4> bkg{{51,77,102,255}}; colors->SetColor("BkgColor", bkg.data());
vtkNew<vtkRenderWindow> renWin;
vtkNew<vtkRenderWindowInteractor> 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<vtkCamera> 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<vtkRenderer> 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<vtkImageData> imageData;
imageData->SetDimensions(latt_size[XYZDims[0]], latt_size[XYZDims[1]], latt_size[XYZDims[2]]);
imageData->AllocateScalars(VTK_DOUBLE, 1);
for (int xx=0;xx<latt_size[XYZDims[0]];xx++)
for (int yy=0;yy<latt_size[XYZDims[1]];yy++)
for (int zz=0;zz<latt_size[XYZDims[2]];zz++)
imageData->SetScalarComponentFromDouble(xx,yy,zz,0,0.0);
vtkNew<isosurface> posExtractor; posExtractor->SetInputData(imageData); posExtractor->SetValue(0, contour);
vtkNew<vtkStripper> posStripper; posStripper->SetInputConnection(posExtractor->GetOutputPort());
vtkNew<vtkPolyDataMapper> posMapper; posMapper->SetInputConnection(posStripper->GetOutputPort()); posMapper->ScalarVisibilityOff();
vtkNew<vtkActor> 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<isosurface> negExtractor; negExtractor->SetInputData(imageData); negExtractor->SetValue(0, -contour);
vtkNew<vtkStripper> negStripper; negStripper->SetInputConnection(negExtractor->GetOutputPort());
vtkNew<vtkPolyDataMapper> negMapper; negMapper->SetInputConnection(negStripper->GetOutputPort()); negMapper->ScalarVisibilityOff();
vtkNew<vtkActor> neg; neg->SetMapper(negMapper);
neg->GetProperty()->SetDiffuseColor(colors->GetColor3d("Ivory").GetData());
vtkNew<vtkOutlineFilter> outlineData; outlineData->SetInputData(imageData);
vtkNew<vtkPolyDataMapper> mapOutline; mapOutline->SetInputConnection(outlineData->GetOutputPort());
vtkNew<vtkActor> outline; outline->SetMapper(mapOutline);
outline->GetProperty()->SetColor(colors->GetColor3d("Black").GetData());
vtkNew<vtkTextActor> 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<vtkTextActor> 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<vtkTextActor> 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<FrameUpdater> 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<CommandHandler> 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;
}
-633
View File
@@ -1,633 +0,0 @@
// 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/<label>_<idx>.png.
// These can be read back by Claude to interpret spatial structure.
//
// Usage:
// ForceAnalysis --grid 32.32.32.32 --mpi 1.1.1.1
// --snapdir /path/to/snapshots
// --first 0 --last 1920 --step 10
// --hotsite x.y.z.t
// --pngdir /path/to/output/pngs
// --isosurface 1.0 (contour in units of field RMS)
// --fixediso 0.05 (fixed absolute contour, overrides --isosurface)
// --slice t (which dimension to fix for 3D display, default: t)
// --sliceval 2 (value of that dimension, default: 0)
//
// Dimension order on the 32^4 lattice: x=0 y=1 z=2 t=3
#include <vtkActor.h>
#include <vtkActor2D.h>
#include <vtkCamera.h>
#include <vtkImageActor.h>
#include <vtkImageMapper3D.h>
#include <vtkImageData.h>
#include <vtkImageMapToColors.h>
#include <vtkLookupTable.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkOutlineFilter.h>
#include <vtkPolyData.h>
#include <vtkPolyDataMapper.h>
#include <vtkPolyDataMapper2D.h>
#include <vtkProperty.h>
#include <vtkProperty2D.h>
#include <vtkPoints.h>
#include <vtkCellArray.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkStripper.h>
#include <vtkCallbackCommand.h>
#include <vtkTextActor.h>
#include <vtkTextProperty.h>
#include <vtkWindowToImageFilter.h>
#include <vtkPNGWriter.h>
#define USE_FLYING_EDGES
#ifdef USE_FLYING_EDGES
#include <vtkFlyingEdges3D.h>
typedef vtkFlyingEdges3D isosurface;
#else
#include <vtkMarchingCubes.h>
typedef vtkMarchingCubes isosurface;
#endif
#include <Grid/Grid.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <string>
#include <vector>
#include <cmath>
#include <memory>
#include <sys/stat.h>
using namespace Grid;
// ─── I/O ─────────────────────────────────────────────────────────────────────
template <class T>
bool tryReadFile(T& out, const std::string& fname)
{
std::ifstream test(fname);
if (!test.good()) return false;
test.close();
emptyUserRecord record;
ScidacReader RD;
RD.open(fname);
RD.readScidacFieldRecord(out, record);
RD.close();
return true;
}
// ─── Fill a 3D vtkImageData slice from a 4D lattice field ────────────────────
// Sums over the sliced dimension at sliceval, displays the remaining 3 dims.
void fillImageData(vtkImageData* img,
LatticeComplexD& field,
const Coordinate& latt_size,
int slice_dim, int sliceval)
{
// Display dims = all dims except slice_dim, in order
std::vector<int> disp;
for (int d = 0; d < 4; d++) if (d != slice_dim) disp.push_back(d);
int Nx = latt_size[disp[0]];
int Ny = latt_size[disp[1]];
int Nz = latt_size[disp[2]];
for (int ix = 0; ix < Nx; ix++)
for (int iy = 0; iy < Ny; iy++)
for (int iz = 0; iz < Nz; iz++) {
Coordinate site(4);
site[disp[0]] = ix;
site[disp[1]] = iy;
site[disp[2]] = iz;
site[slice_dim] = sliceval;
RealD val = real(peekSite(field, site));
img->SetScalarComponentFromDouble(ix, iy, iz, 0, val);
}
img->Modified();
}
// ─── 2D heatmap: persistent context ───────────────────────────────────────────
// Renders a fixed (dim1=v1, dim2=v2) slice of the 4D force field as a
// diverging blue→white→red colour map, with a fixed symmetric colour scale
// so brightness directly encodes force magnitude across all frames.
// A white cross-hair marks the hotsite projection onto the slice.
struct HeatmapCtx {
// image pipeline
vtkNew<vtkImageData> img;
vtkNew<vtkLookupTable> lut;
vtkNew<vtkImageMapToColors> colorMap;
vtkNew<vtkImageActor> imgActor;
// colour scale legend (text, avoids needing RenderingAnnotation module)
vtkNew<vtkTextActor> cbar;
// hotsite cross-hair (2D overlay actors)
vtkNew<vtkPolyData> crossPD;
vtkNew<vtkPoints> crossPts;
vtkNew<vtkCellArray> crossLines;
vtkNew<vtkActor2D> crossActor;
// title
vtkNew<vtkTextActor> titleAct;
// renderer / window
vtkNew<vtkRenderer> ren;
vtkNew<vtkRenderWindow> renWin;
vtkNew<vtkWindowToImageFilter> w2i;
vtkNew<vtkPNGWriter> writer;
int Nx = 0, Ny = 0; // display dimensions of the slice
double scale = 0.07; // colour range: [-scale, +scale]
int hotX = -1, hotY = -1; // hotsite projection onto (Nx,Ny) plane
// pixel coords of the image origin in the render window
int imgOffX = 60, imgOffY = 40;
int imgW = 0, imgH = 0; // rendered pixel size of each lattice cell
void init(int nx, int ny, double sc, int hx, int hy)
{
Nx = nx; Ny = ny; scale = sc;
hotX = hx; hotY = hy;
const int WIN_W = 900, WIN_H = 700;
// Make cells square and as large as possible within the central area
int cellW = (WIN_W - 160) / Nx;
int cellH = (WIN_H - 120) / Ny;
imgW = std::min(cellW, cellH);
imgH = imgW;
imgOffX = (WIN_W - Nx * imgW) / 2;
imgOffY = 60;
// --- Image data (scalar field, one component) ---
img->SetDimensions(Nx, Ny, 1);
img->SetSpacing(imgW, imgH, 1);
img->SetOrigin(imgOffX, imgOffY, 0);
img->AllocateScalars(VTK_DOUBLE, 1);
// --- Diverging LUT: blue(-scale) → white(0) → red(+scale) ---
lut->SetNumberOfTableValues(512);
lut->SetRange(-scale, scale);
lut->SetNanColor(0.2, 0.2, 0.2, 1.0);
for (int i = 0; i < 512; ++i) {
double t = i / 511.0; // 0=blue, 0.5=white, 1=red
double r = (t > 0.5) ? 1.0 : 2.0 * t;
double g = (t < 0.5) ? 2.0 * t : 2.0 * (1.0 - t);
double b = (t < 0.5) ? 1.0 : 2.0 * (1.0 - t);
lut->SetTableValue(i, r, g, b, 1.0);
}
lut->Build();
// --- Colour map pipeline ---
colorMap->SetInputData(img);
colorMap->SetLookupTable(lut);
colorMap->Update();
imgActor->GetMapper()->SetInputConnection(colorMap->GetOutputPort());
// --- Colour scale legend (text) ---
{
std::ostringstream ss;
ss << std::scientific << std::setprecision(2)
<< "blue=-" << sc << " white=0 red=+" << sc;
cbar->SetInput(ss.str().c_str());
}
cbar->GetTextProperty()->SetFontFamilyToCourier();
cbar->GetTextProperty()->SetFontSize(13);
cbar->GetTextProperty()->SetColor(0.9, 0.9, 0.9);
cbar->SetDisplayPosition(10, 10);
// --- Cross-hair at hotsite (2D display coords) ---
if (hotX >= 0 && hotY >= 0) {
double cx = imgOffX + (hotX + 0.5) * imgW;
double cy = imgOffY + (hotY + 0.5) * imgH;
double arm = imgW * 0.8;
crossPts->InsertNextPoint(cx - arm, cy, 0);
crossPts->InsertNextPoint(cx + arm, cy, 0);
crossPts->InsertNextPoint(cx, cy - arm, 0);
crossPts->InsertNextPoint(cx, cy + arm, 0);
vtkIdType seg0[2] = {0, 1};
vtkIdType seg1[2] = {2, 3};
crossLines->InsertNextCell(2, seg0);
crossLines->InsertNextCell(2, seg1);
crossPD->SetPoints(crossPts);
crossPD->SetLines(crossLines);
vtkNew<vtkPolyDataMapper2D> crossMap;
crossMap->SetInputData(crossPD);
crossActor->SetMapper(crossMap);
crossActor->GetProperty()->SetColor(1, 1, 1);
crossActor->GetProperty()->SetLineWidth(2.0);
}
// --- Title ---
titleAct->GetTextProperty()->SetFontFamilyToCourier();
titleAct->GetTextProperty()->SetFontSize(16);
titleAct->GetTextProperty()->SetColor(1, 1, 0);
titleAct->SetDisplayPosition(10, WIN_H - 30);
// --- Renderer (2D parallel projection so image fills correctly) ---
ren->SetBackground(0.08, 0.08, 0.12);
ren->AddActor(imgActor);
ren->AddActor2D(cbar);
ren->AddActor2D(crossActor);
ren->AddActor2D(titleAct);
ren->GetActiveCamera()->ParallelProjectionOn();
// Set up camera to look straight down at the image plane
ren->GetActiveCamera()->SetPosition(WIN_W/2.0, WIN_H/2.0, 1000);
ren->GetActiveCamera()->SetFocalPoint(WIN_W/2.0, WIN_H/2.0, 0);
ren->GetActiveCamera()->SetViewUp(0, 1, 0);
ren->GetActiveCamera()->SetParallelScale(WIN_H / 2.0);
ren->ResetCameraClippingRange();
renWin->AddRenderer(ren);
renWin->SetSize(WIN_W, WIN_H);
renWin->SetOffScreenRendering(1);
renWin->SetMultiSamples(0);
w2i->SetInput(renWin);
w2i->SetInputBufferTypeToRGB();
w2i->ReadFrontBufferOff();
}
};
void renderHeatmap(HeatmapCtx& ctx,
LatticeComplexD& field,
const Coordinate& latt_size,
int dim1, int val1, // first fixed dimension
int dim2, int val2, // second fixed dimension
const std::string& title,
const std::string& outpath)
{
// Display dimensions: the two dims that are NOT fixed
std::vector<int> disp;
for (int d = 0; d < 4; d++)
if (d != dim1 && d != dim2) disp.push_back(d);
int Nx = latt_size[disp[0]];
int Ny = latt_size[disp[1]];
// Fill image data
for (int ix = 0; ix < Nx; ix++) {
for (int iy = 0; iy < Ny; iy++) {
Coordinate site(4);
site[disp[0]] = ix;
site[disp[1]] = iy;
site[dim1] = val1;
site[dim2] = val2;
RealD val = real(TensorRemove(peekSite(field, site)));
ctx.img->SetScalarComponentFromDouble(ix, iy, 0, 0, val);
}
}
ctx.img->Modified();
ctx.colorMap->Update();
ctx.titleAct->SetInput(title.c_str());
ctx.renWin->Render();
ctx.w2i->Modified();
ctx.w2i->Update();
ctx.writer->SetFileName(outpath.c_str());
ctx.writer->SetInputConnection(ctx.w2i->GetOutputPort());
ctx.writer->Write();
}
// ─── Persistent rendering context (created once, reused every frame) ──────────
// Avoids Metal GPU context exhaustion on macOS when rendering hundreds of frames.
struct RenderCtx {
vtkNew<vtkNamedColors> colors;
vtkNew<vtkImageData> imageData;
vtkNew<isosurface> posEx, negEx;
vtkNew<vtkStripper> posSt, negSt;
vtkNew<vtkPolyDataMapper> posMap, negMap, outMap;
vtkNew<vtkActor> posAct, negAct, outAct;
vtkNew<vtkOutlineFilter> outF;
vtkNew<vtkTextActor> label;
vtkNew<vtkRenderer> ren;
vtkNew<vtkCamera> cam;
vtkNew<vtkRenderWindow> renWin;
vtkNew<vtkWindowToImageFilter> w2i;
vtkNew<vtkPNGWriter> writer;
void init(int Nx, int Ny, int Nz)
{
std::array<unsigned char,4> posColor{{240,184,160,255}};
colors->SetColor("posColor", posColor.data());
std::array<unsigned char,4> bkg{{51,77,102,255}};
colors->SetColor("BkgColor", bkg.data());
imageData->SetDimensions(Nx, Ny, Nz);
imageData->AllocateScalars(VTK_DOUBLE, 1);
posEx->SetInputData(imageData); posEx->SetValue(0, 1.0);
posSt->SetInputConnection(posEx->GetOutputPort());
posMap->SetInputConnection(posSt->GetOutputPort());
posMap->ScalarVisibilityOff();
posAct->SetMapper(posMap);
posAct->GetProperty()->SetDiffuseColor(colors->GetColor3d("posColor").GetData());
posAct->GetProperty()->SetSpecular(0.3);
posAct->GetProperty()->SetSpecularPower(20);
posAct->GetProperty()->SetOpacity(0.6);
negEx->SetInputData(imageData); negEx->SetValue(0, -1.0);
negSt->SetInputConnection(negEx->GetOutputPort());
negMap->SetInputConnection(negSt->GetOutputPort());
negMap->ScalarVisibilityOff();
negAct->SetMapper(negMap);
negAct->GetProperty()->SetDiffuseColor(colors->GetColor3d("Ivory").GetData());
negAct->GetProperty()->SetOpacity(0.6);
outF->SetInputData(imageData);
outMap->SetInputConnection(outF->GetOutputPort());
outAct->SetMapper(outMap);
outAct->GetProperty()->SetColor(colors->GetColor3d("Black").GetData());
label->SetPosition(10, 10);
label->GetTextProperty()->SetFontFamilyToCourier();
label->GetTextProperty()->SetFontSize(18);
label->GetTextProperty()->SetColor(colors->GetColor3d("Gold").GetData());
ren->AddActor(posAct);
ren->AddActor(negAct);
ren->AddActor(outAct);
ren->AddActor2D(label);
ren->SetBackground(colors->GetColor3d("BkgColor").GetData());
cam->SetViewUp(0,0,-1);
cam->SetPosition(0,-1000,0);
cam->SetFocalPoint(0,0,0);
cam->ComputeViewPlaneNormal();
cam->Azimuth(30.0);
cam->Elevation(30.0);
ren->SetActiveCamera(cam);
renWin->AddRenderer(ren);
renWin->SetSize(800, 600);
renWin->SetOffScreenRendering(1);
renWin->SetMultiSamples(0);
w2i->SetInput(renWin);
w2i->SetInputBufferTypeToRGB();
w2i->ReadFrontBufferOff();
}
};
// ─── Render one force field snapshot to a PNG (reuses existing RenderCtx) ─────
void renderPNG(RenderCtx& ctx,
LatticeComplexD& field,
const Coordinate& latt_size,
int slice_dim, int sliceval,
double contour,
const std::string& title,
const std::string& outpath)
{
// Update image data
fillImageData(ctx.imageData, field, latt_size, slice_dim, sliceval);
// Update isosurface levels
ctx.posEx->SetValue(0, contour);
ctx.negEx->SetValue(0, -contour);
// Update label
ctx.label->SetInput(title.c_str());
// Reset camera to fit the (possibly new) data bounds
ctx.ren->ResetCamera();
ctx.cam->Dolly(1.2);
ctx.ren->ResetCameraClippingRange();
ctx.renWin->Render();
ctx.w2i->Modified();
ctx.w2i->Update();
ctx.writer->SetFileName(outpath.c_str());
ctx.writer->SetInputConnection(ctx.w2i->GetOutputPort());
ctx.writer->Write();
}
// ─── main ─────────────────────────────────────────────────────────────────────
int main(int argc, char* argv[])
{
Grid_init(&argc, &argv);
GridLogMessage.Active(0);
GridLogIterative.Active(0);
GridLogDebug.Active(0);
GridLogPerformance.Active(0);
GridLogComms.Active(0);
GridLogDslash.Active(0);
GridLogMemory.Active(0);
// ── CLI ──────────────────────────────────────────────────────────────────
std::string snapdir = ".";
std::string pngdir = "";
int first = 0, last = 1920, step = 1;
int slice_dim = 3, sliceval = 0; // default: fix t=0, display xyz
double iso_rms = 1.0;
double fixed_iso = -1.0; // if >0, use this absolute contour
double tau_start = -1.0; // if >=0, display MD time tau = tau_start + idx*tau_step
double tau_step = 0.0;
// Heatmap mode: fix two dimensions, show 2D colour map
bool do_heatmap = false;
int slice_dim2 = -1, sliceval2 = 0;
double heat_scale = -1.0; // if >0, fixed symmetric colour scale; else auto
Coordinate hotsite({0,0,0,0});
bool has_hotsite = false;
std::string arg;
if (GridCmdOptionExists(argv, argv+argc, "--snapdir"))
snapdir = GridCmdOptionPayload(argv, argv+argc, "--snapdir");
if (GridCmdOptionExists(argv, argv+argc, "--pngdir"))
pngdir = GridCmdOptionPayload(argv, argv+argc, "--pngdir");
if (GridCmdOptionExists(argv, argv+argc, "--first")) {
arg = GridCmdOptionPayload(argv, argv+argc, "--first");
GridCmdOptionInt(arg, first);
}
if (GridCmdOptionExists(argv, argv+argc, "--last")) {
arg = GridCmdOptionPayload(argv, argv+argc, "--last");
GridCmdOptionInt(arg, last);
}
if (GridCmdOptionExists(argv, argv+argc, "--step")) {
arg = GridCmdOptionPayload(argv, argv+argc, "--step");
GridCmdOptionInt(arg, step);
}
if (GridCmdOptionExists(argv, argv+argc, "--slicedim")) {
arg = GridCmdOptionPayload(argv, argv+argc, "--slicedim");
GridCmdOptionInt(arg, slice_dim);
}
if (GridCmdOptionExists(argv, argv+argc, "--sliceval")) {
arg = GridCmdOptionPayload(argv, argv+argc, "--sliceval");
GridCmdOptionInt(arg, sliceval);
}
if (GridCmdOptionExists(argv, argv+argc, "--isosurface")) {
arg = GridCmdOptionPayload(argv, argv+argc, "--isosurface");
GridCmdOptionFloat(arg, iso_rms);
}
if (GridCmdOptionExists(argv, argv+argc, "--fixediso")) {
arg = GridCmdOptionPayload(argv, argv+argc, "--fixediso");
GridCmdOptionFloat(arg, fixed_iso);
}
if (GridCmdOptionExists(argv, argv+argc, "--taustart")) {
arg = GridCmdOptionPayload(argv, argv+argc, "--taustart");
GridCmdOptionFloat(arg, tau_start);
}
if (GridCmdOptionExists(argv, argv+argc, "--taustep")) {
arg = GridCmdOptionPayload(argv, argv+argc, "--taustep");
GridCmdOptionFloat(arg, tau_step);
}
if (GridCmdOptionExists(argv, argv+argc, "--hotsite")) {
arg = GridCmdOptionPayload(argv, argv+argc, "--hotsite");
GridCmdOptionIntVector(arg, hotsite);
has_hotsite = true;
}
if (GridCmdOptionExists(argv, argv+argc, "--heatmap"))
do_heatmap = true;
if (GridCmdOptionExists(argv, argv+argc, "--slicedim2")) {
arg = GridCmdOptionPayload(argv, argv+argc, "--slicedim2");
GridCmdOptionInt(arg, slice_dim2);
}
if (GridCmdOptionExists(argv, argv+argc, "--sliceval2")) {
arg = GridCmdOptionPayload(argv, argv+argc, "--sliceval2");
GridCmdOptionInt(arg, sliceval2);
}
if (GridCmdOptionExists(argv, argv+argc, "--heatscale")) {
arg = GridCmdOptionPayload(argv, argv+argc, "--heatscale");
GridCmdOptionFloat(arg, heat_scale);
}
bool do_png = !pngdir.empty();
if (do_png) mkdir(pngdir.c_str(), 0755);
// ── Grid setup ───────────────────────────────────────────────────────────
auto latt_size = GridDefaultLatt();
auto simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
auto mpi_layout = GridDefaultMpi();
GridCartesian grid(latt_size, simd_layout, mpi_layout);
LatticeComplexD field(&grid);
// Force components
struct ForceSpec { std::string prefix; std::string label; };
std::vector<ForceSpec> forces = {
{ "F_IwasakiGaugeAction_lat", "Gauge_lat" },
{ "F_IwasakiGaugeAction_smr", "Gauge_smr" },
{ "F_JacobianAction_smr", "Jacobian" },
{ "F_TwoFlavourEvenOddRatioPseudoFermionActiondet_0.0047_det_0.05_lat", "Ferm0047_lat" },
{ "F_TwoFlavourEvenOddRatioPseudoFermionActiondet_0.0047_det_0.05_smr", "Ferm0047_smr" },
{ "F_TwoFlavourEvenOddRatioPseudoFermionActiondet_0.05_det_0.1_lat", "Ferm005_lat" },
{ "F_TwoFlavourEvenOddRatioPseudoFermionActiondet_0.1_det_0.25_lat", "Ferm01_lat" },
{ "F_TwoFlavourEvenOddRatioPseudoFermionActiondet_0.25_det_0.5_lat", "Ferm025_lat" },
{ "F_TwoFlavourEvenOddRatioPseudoFermionActiondet_0.5_det_1_lat", "Ferm05_lat" },
};
// ── Stdout header ─────────────────────────────────────────────────────────
std::cerr << "idx";
for (auto& fs : forces) {
std::cerr << "\t" << fs.label << "_rms";
if (has_hotsite) std::cerr << "\t" << fs.label << "_hot";
}
std::cerr << "\n";
// ── Persistent render contexts (one GPU context for all frames) ──────────
std::unique_ptr<RenderCtx> ctx; // isosurface mode
std::unique_ptr<HeatmapCtx> hctx; // heatmap mode
// ── Main loop ─────────────────────────────────────────────────────────────
for (int idx = first; idx <= last; idx += step) {
std::cerr << idx;
for (auto& fs : forces) {
std::string fname = snapdir + "/" + fs.prefix + "." + std::to_string(idx);
if (!tryReadFile(field, fname)) {
std::cerr << "\t-";
if (has_hotsite) std::cerr << "\t-";
continue;
}
// RMS (real part)
RealD sumsq = 0.0;
for (int i = 0; i < grid.gSites(); i++) {
Coordinate site;
Lexicographic::CoorFromIndex(site, i, latt_size);
RealD v = real(peekSite(field, site));
sumsq += v * v;
}
RealD rms = std::sqrt(sumsq / grid.gSites());
std::cerr << "\t" << rms;
if (has_hotsite) {
RealD hval = real(TensorRemove(peekSite(field, hotsite)));
std::cerr << "\t" << hval;
}
// PNG output (isosurface or heatmap)
if (do_png) {
// Build title string
std::ostringstream title;
title << fs.label << " ";
if (tau_start >= 0.0 && tau_step > 0.0) {
double tau = tau_start + idx * tau_step;
title << std::fixed << std::setprecision(6) << "tau=" << tau;
} else {
title << "idx=" << idx;
}
title << " rms=" << std::scientific << std::setprecision(3) << rms;
std::ostringstream outpath;
outpath << pngdir << "/" << fs.label
<< "_" << std::setfill('0') << std::setw(6) << idx << ".png";
if (do_heatmap && slice_dim2 >= 0) {
// ── Heatmap mode ────────────────────────────────────────
// Display dims = the two that are NOT fixed
std::vector<int> disp;
for (int d = 0; d < 4; d++)
if (d != slice_dim && d != slice_dim2) disp.push_back(d);
if (!hctx) {
double sc = (heat_scale > 0) ? heat_scale : rms * 20.0;
// Hotsite projection onto display plane
int hx = -1, hy = -1;
if (has_hotsite) {
hx = hotsite[disp[0]];
hy = hotsite[disp[1]];
}
hctx = std::make_unique<HeatmapCtx>();
hctx->init(latt_size[disp[0]], latt_size[disp[1]], sc, hx, hy);
}
title << " scale=+-" << std::fixed << std::setprecision(4) << hctx->scale;
renderHeatmap(*hctx, field, latt_size,
slice_dim, sliceval,
slice_dim2, sliceval2,
title.str(), outpath.str());
} else {
// ── Isosurface mode ─────────────────────────────────────
double contour = (fixed_iso > 0) ? fixed_iso : iso_rms * rms;
title << " iso=" << contour;
if (!ctx) {
std::vector<int> disp;
for (int d = 0; d < 4; d++) if (d != slice_dim) disp.push_back(d);
ctx = std::make_unique<RenderCtx>();
ctx->init(latt_size[disp[0]], latt_size[disp[1]], latt_size[disp[2]]);
}
renderPNG(*ctx, field, latt_size, slice_dim, sliceval,
contour, title.str(), outpath.str());
}
}
}
std::cerr << "\n";
std::cerr.flush();
}
Grid_finalize();
return EXIT_SUCCESS;
}
-742
View File
@@ -1,742 +0,0 @@
// TranscriptToVideo.cxx
//
// Reads a conversation transcript file with [User] / [Claude] turns and
// renders it to an AVI using vtkFFMPEGWriter at 1280x720, 10 fps.
//
// Transcript format:
// [USER] Some question or command, possibly spanning
// multiple continuation lines.
// [ASSISTANT] A response, also possibly
// spanning multiple lines.
// ...
//
// Rules:
// - A line beginning "[User]" or "[Claude]" starts a new turn.
// - Any subsequent non-blank line that does NOT begin with "[" is a
// continuation of the previous turn (joined with a single space).
// - Blank lines are ignored.
//
// Usage:
// ./TranscriptToVideo <transcript.txt> <output.avi>
//
// Typewriter speed : 10 chars/sec → 1 frame/char at 10 fps
// Pause after turn : 0.5 s → 5 frames
// Word-wrap column : 62
#include <vtkActor.h>
#include <vtkActor2D.h>
#include <vtkCellArray.h>
#include <vtkFFMPEGWriter.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkPoints.h>
#include <vtkPolyData.h>
#include <vtkPolyDataMapper2D.h>
#include <vtkProperty2D.h>
#include <vtkRenderWindow.h>
#include <vtkRenderer.h>
#include <vtkTextActor.h>
#include <vtkTextProperty.h>
#include <vtkWindowToImageFilter.h>
#include <algorithm>
#include <fstream>
#include <iostream>
#include <sstream>
#include <string>
#include <vector>
// ---------------------------------------------------------------------------
// Constants
// ---------------------------------------------------------------------------
static const int WIDTH = 1280;
static const int HEIGHT = 720;
static const int FPS = 10;
static const int CHARS_PER_FRAME = 5; // 50 chars/sec at 10 fps
static const int PAUSE_FRAMES = 3; // 0.3 s
static const int WRAP_COLS = 60;
static const int MAX_HISTORY = 18; // visible completed lines
static const int FONT_SIZE = 18;
static const int TITLE_SIZE = 23;
static const int LINE_HEIGHT = 28; // pixels between lines
static const int MARGIN_LEFT = 48;
static const int MARGIN_TOP = 58; // below title bar
// Colours (R, G, B in 01)
static const double COL_BG[3] = { 0.04, 0.04, 0.10 }; // near-black navy
static const double COL_USER[3] = { 1.00, 0.84, 0.00 }; // gold
static const double COL_CLAUDE[3] = { 0.68, 0.85, 0.90 }; // light steel blue
static const double COL_TITLE[3] = { 1.00, 1.00, 1.00 }; // white
static const double COL_BAR[3] = { 0.30, 0.55, 0.80 }; // progress bar blue
static const double COL_LABEL[3] = { 0.65, 0.65, 0.65 }; // dim grey for speaker tag
// ---------------------------------------------------------------------------
// Structs
// ---------------------------------------------------------------------------
enum class Speaker { User, Claude, Caption };
struct Turn {
Speaker speaker;
std::string text; // full unwrapped text
};
// A rendered display line (already word-wrapped, tagged with speaker)
struct DisplayLine {
Speaker speaker;
std::string prefix; // "[USER] " or "[ASSISTANT] "
std::string body; // wrapped segment
bool isFirst; // first line of this turn (prefix printed)
bool isBlank; // spacer between turns (no text rendered)
bool isCaption; // caption update — body holds text (empty = clear)
DisplayLine() : speaker(Speaker::User), isFirst(false), isBlank(false), isCaption(false) {}
};
// ---------------------------------------------------------------------------
// Replace common multi-byte UTF-8 sequences with ASCII equivalents so that
// VTK's font renderer (which only handles ASCII reliably) does not crash.
// Any remaining non-ASCII byte is replaced with '?'.
// ---------------------------------------------------------------------------
static std::string SanitiseASCII(const std::string& s)
{
std::string out;
out.reserve(s.size());
const unsigned char* p = (const unsigned char*)s.data();
const unsigned char* end = p + s.size();
while (p < end) {
unsigned char c = *p;
if (c < 0x80) {
out += (char)c;
++p;
} else if (c == 0xE2 && (p + 2) < end) {
// 3-byte sequence starting 0xE2
unsigned char b1 = *(p+1), b2 = *(p+2);
// U+2018 ' U+2019 ' (0xE2 0x80 0x98 / 0x99)
if (b1 == 0x80 && (b2 == 0x98 || b2 == 0x99)) { out += '\''; p += 3; }
// U+201C " U+201D " (0xE2 0x80 0x9C / 0x9D)
else if (b1 == 0x80 && (b2 == 0x9C || b2 == 0x9D)) { out += '"'; p += 3; }
// U+2013 en-dash (0xE2 0x80 0x93)
else if (b1 == 0x80 && b2 == 0x93) { out += '-'; p += 3; }
// U+2014 em-dash (0xE2 0x80 0x94)
else if (b1 == 0x80 && b2 == 0x94) { out += '-'; p += 3; }
// U+2026 ellipsis (0xE2 0x80 0xA6)
else if (b1 == 0x80 && b2 == 0xA6) { out += "..."; p += 3; }
// U+2192 arrow (0xE2 0x86 0x92)
else if (b1 == 0x86 && b2 == 0x92) { out += "->"; p += 3; }
// U+00D7 multiplication sign (0xC3 0x97) — caught below, but
// U+22xx math operators: replace with '~'
else { out += '~'; p += 3; }
} else if (c == 0xC3 && (p + 1) < end) {
// 2-byte latin-1 supplement
unsigned char b1 = *(p+1);
if (b1 == 0x97) { out += 'x'; p += 2; } // U+00D7 ×
else if (b1 == 0xB7) { out += '/'; p += 2; } // U+00F7 ÷
else { out += '?'; p += 2; }
} else {
// Unknown multi-byte: skip the whole sequence
out += '?';
++p;
while (p < end && (*p & 0xC0) == 0x80) ++p; // skip continuation bytes
}
}
return out;
}
// ---------------------------------------------------------------------------
// Strip markdown emphasis markers (* and `) so they don't appear in the video.
// ---------------------------------------------------------------------------
static std::string StripMarkdown(const std::string& s)
{
std::string out;
out.reserve(s.size());
for (char c : s) {
if (c == '*' || c == '`') continue;
out += c;
}
return out;
}
// ---------------------------------------------------------------------------
// Word-wrap: splits `text` into lines of at most maxCols characters.
// ---------------------------------------------------------------------------
static std::vector<std::string> WrapText(const std::string& text, int maxCols)
{
std::vector<std::string> lines;
std::istringstream words(text);
std::string word, line;
while (words >> word) {
if (!line.empty() && (int)(line.size() + 1 + word.size()) > maxCols) {
lines.push_back(line);
line = word;
} else {
if (!line.empty()) line += ' ';
line += word;
}
}
if (!line.empty()) lines.push_back(line);
if (lines.empty()) lines.push_back("");
return lines;
}
// ---------------------------------------------------------------------------
// Parse transcript file into list of Turns.
// ---------------------------------------------------------------------------
static std::vector<Turn> ParseTranscript(const std::string& path)
{
std::ifstream f(path);
if (!f) {
std::cerr << "Cannot open transcript: " << path << "\n";
std::exit(1);
}
std::vector<Turn> turns;
std::string line;
while (std::getline(f, line)) {
// Trim trailing CR (Windows files)
if (!line.empty() && line.back() == '\r') line.pop_back();
if (line.empty()) {
// Blank line within a dialogue turn = paragraph break.
// Caption turns don't support paragraph breaks.
if (!turns.empty() && turns.back().speaker != Speaker::Caption) {
turns.back().text += '\n'; // '\n' sentinel: expands to blank spacer
}
continue;
}
if (line.size() >= 6 && line.substr(0, 6) == "[USER]") {
Turn t;
t.speaker = Speaker::User;
t.text = line.substr(6);
while (!t.text.empty() && t.text.front() == ' ') t.text.erase(t.text.begin());
t.text = SanitiseASCII(t.text);
turns.push_back(std::move(t));
} else if (line.size() >= 11 && line.substr(0, 11) == "[ASSISTANT]") {
Turn t;
t.speaker = Speaker::Claude;
t.text = line.substr(11);
while (!t.text.empty() && t.text.front() == ' ') t.text.erase(t.text.begin());
t.text = SanitiseASCII(t.text);
turns.push_back(std::move(t));
} else if (line.size() >= 9 && line.substr(0, 9) == "[CAPTION]") {
Turn t;
t.speaker = Speaker::Caption;
t.text = line.substr(9);
while (!t.text.empty() && t.text.front() == ' ') t.text.erase(t.text.begin());
t.text = SanitiseASCII(t.text);
turns.push_back(std::move(t));
} else if (!turns.empty()) {
// Continuation line — strip leading whitespace, preserve as separate line
size_t start = line.find_first_not_of(" \t");
if (start != std::string::npos) {
if (!turns.back().text.empty()) turns.back().text += '\n';
turns.back().text += SanitiseASCII(line.substr(start));
}
}
}
return turns;
}
// ---------------------------------------------------------------------------
// Expand all Turns into DisplayLines (word-wrapped).
// ---------------------------------------------------------------------------
static std::vector<DisplayLine> ExpandToDisplayLines(const std::vector<Turn>& turns)
{
std::vector<DisplayLine> out;
// Prefix widths kept equal for alignment
const std::string userPfx = "[USER] ";
const std::string claudePfx = "[ASSISTANT] ";
// Track previous non-caption speaker to know when to insert blank spacers
Speaker prevSpeaker = Speaker::Caption; // sentinel: no spacer before first real turn
for (size_t ti = 0; ti < turns.size(); ++ti) {
const Turn& t = turns[ti];
// Caption turn: emit one special DisplayLine, no spacer, no history entry
if (t.speaker == Speaker::Caption) {
DisplayLine dl;
dl.isCaption = true;
dl.speaker = Speaker::Caption;
// body: whitespace-only → clear; otherwise wrap lines joined with \n
std::string trimmed = t.text;
size_t first = trimmed.find_first_not_of(" \t\r\n");
if (first == std::string::npos) {
dl.body = ""; // signal to clear caption
} else {
// Wrap to ~90 cols for the wider caption zone
auto lines = WrapText(trimmed, 90);
for (size_t i = 0; i < lines.size(); ++i) {
if (i > 0) dl.body += "\n";
dl.body += lines[i];
}
}
out.push_back(dl);
continue;
}
// Insert a blank spacer before each new dialogue turn (not before the first)
if (prevSpeaker != Speaker::Caption) {
DisplayLine blank;
blank.isBlank = true;
blank.speaker = t.speaker;
out.push_back(blank);
}
prevSpeaker = t.speaker;
const std::string& pfx = (t.speaker == Speaker::User) ? userPfx : claudePfx;
int bodyWidth = WRAP_COLS - (int)pfx.size();
if (bodyWidth < 20) bodyWidth = 20;
// Split the turn text on '\n' to get individual source lines.
// Empty source lines become blank spacers (paragraph breaks within a turn).
// Non-empty source lines are stripped of markdown markers and word-wrapped.
std::vector<std::string> srcLines;
{
std::string seg;
for (char ch : t.text) {
if (ch == '\n') { srcLines.push_back(seg); seg.clear(); }
else seg += ch;
}
srcLines.push_back(seg);
}
bool firstOfTurn = true;
for (const auto& srcLine : srcLines) {
// Strip markdown emphasis markers
std::string stripped = StripMarkdown(srcLine);
// Trim leading/trailing whitespace
size_t f = stripped.find_first_not_of(" \t");
if (f == std::string::npos) {
// Blank — paragraph break spacer within the turn
DisplayLine blank;
blank.isBlank = true;
blank.speaker = t.speaker;
out.push_back(blank);
continue;
}
size_t l = stripped.find_last_not_of(" \t");
stripped = stripped.substr(f, l - f + 1);
auto wrapped = WrapText(stripped, bodyWidth);
for (size_t i = 0; i < wrapped.size(); ++i) {
DisplayLine dl;
dl.speaker = t.speaker;
dl.prefix = pfx;
dl.body = wrapped[i];
dl.isFirst = firstOfTurn && (i == 0);
out.push_back(dl);
}
firstOfTurn = false;
}
}
return out;
}
// ---------------------------------------------------------------------------
// Helper: set actor text to `s`, configure font/colour, position.
// ---------------------------------------------------------------------------
static void ConfigureTextActor(vtkTextActor* a, int fontSize,
double r, double g, double b)
{
a->GetTextProperty()->SetFontFamilyToCourier();
a->GetTextProperty()->SetFontSize(fontSize);
a->GetTextProperty()->SetColor(r, g, b);
a->GetTextProperty()->SetBold(0);
a->GetTextProperty()->SetItalic(0);
a->GetTextProperty()->ShadowOff();
a->GetTextProperty()->SetJustificationToLeft();
a->GetTextProperty()->SetVerticalJustificationToBottom();
}
// ---------------------------------------------------------------------------
// Create a thin horizontal progress bar actor (2D polygon).
// Returns the actor; caller adds to renderer.
// ---------------------------------------------------------------------------
static vtkActor2D* MakeProgressBar(vtkPolyData*& pd, vtkPoints*& pts)
{
pts = vtkPoints::New();
vtkCellArray* cells = vtkCellArray::New();
// 4 points, updated every frame
pts->SetNumberOfPoints(4);
pts->SetPoint(0, 0, HEIGHT - 8, 0);
pts->SetPoint(1, 0, HEIGHT - 2, 0);
pts->SetPoint(2, 100, HEIGHT - 2, 0);
pts->SetPoint(3, 100, HEIGHT - 8, 0);
vtkIdType quad[4] = { 0, 1, 2, 3 };
cells->InsertNextCell(4, quad);
pd = vtkPolyData::New();
pd->SetPoints(pts);
pd->SetPolys(cells);
cells->Delete();
vtkPolyDataMapper2D* mapper = vtkPolyDataMapper2D::New();
mapper->SetInputData(pd);
vtkActor2D* actor = vtkActor2D::New();
actor->SetMapper(mapper);
actor->GetProperty()->SetColor(COL_BAR[0], COL_BAR[1], COL_BAR[2]);
mapper->Delete();
return actor;
}
// ---------------------------------------------------------------------------
// Main
// ---------------------------------------------------------------------------
int main(int argc, char* argv[])
{
if (argc < 3) {
std::cerr << "Usage: TranscriptToVideo <transcript.txt> <output.avi>\n";
return 1;
}
const std::string transcriptPath = argv[1];
const std::string outputPath = argv[2];
// -----------------------------------------------------------------------
// Parse & expand
// -----------------------------------------------------------------------
auto turns = ParseTranscript(transcriptPath);
auto displayLines = ExpandToDisplayLines(turns);
const int totalLines = (int)displayLines.size();
if (totalLines == 0) {
std::cerr << "Transcript has no parseable turns.\n";
return 1;
}
// Total frames (rough estimate for progress bar denominator)
// Each display line: on average ~40 chars + PAUSE_FRAMES
// We'll compute exact total below after knowing char counts.
int totalFrames = 0;
for (const auto& dl : displayLines) {
int chars = (int)dl.body.size();
totalFrames += (chars + CHARS_PER_FRAME - 1) / CHARS_PER_FRAME + 1 + PAUSE_FRAMES;
}
// -----------------------------------------------------------------------
// VTK Setup
// -----------------------------------------------------------------------
vtkNew<vtkRenderer> ren;
vtkNew<vtkRenderWindow> renWin;
ren->SetBackground(COL_BG[0], COL_BG[1], COL_BG[2]);
renWin->AddRenderer(ren);
renWin->SetSize(WIDTH, HEIGHT);
renWin->SetOffScreenRendering(1);
renWin->SetMultiSamples(0);
// -----------------------------------------------------------------------
// Title actor
// -----------------------------------------------------------------------
vtkNew<vtkTextActor> titleActor;
titleActor->SetInput("Transcript");
ConfigureTextActor(titleActor, TITLE_SIZE,
COL_TITLE[0], COL_TITLE[1], COL_TITLE[2]);
titleActor->GetTextProperty()->SetBold(1);
titleActor->GetTextProperty()->SetJustificationToCentered();
titleActor->SetDisplayPosition(WIDTH / 2, HEIGHT - 36);
ren->AddActor2D(titleActor);
// Thin separator line below title (drawn as a narrow quad)
{
vtkPoints* lpts = vtkPoints::New();
vtkCellArray* lcell = vtkCellArray::New();
lpts->InsertNextPoint(0, HEIGHT - 44, 0);
lpts->InsertNextPoint(WIDTH, HEIGHT - 44, 0);
lpts->InsertNextPoint(WIDTH, HEIGHT - 42, 0);
lpts->InsertNextPoint(0, HEIGHT - 42, 0);
vtkIdType q[4] = {0,1,2,3};
lcell->InsertNextCell(4, q);
vtkPolyData* lpd = vtkPolyData::New();
lpd->SetPoints(lpts);
lpd->SetPolys(lcell);
vtkPolyDataMapper2D* lmap = vtkPolyDataMapper2D::New();
lmap->SetInputData(lpd);
vtkActor2D* lineActor = vtkActor2D::New();
lineActor->SetMapper(lmap);
lineActor->GetProperty()->SetColor(0.3, 0.3, 0.5);
ren->AddActor2D(lineActor);
lpts->Delete(); lcell->Delete(); lpd->Delete();
lmap->Delete(); lineActor->Delete();
}
// -----------------------------------------------------------------------
// History text actors (MAX_HISTORY lines, reused with shifting content)
// -----------------------------------------------------------------------
std::vector<vtkTextActor*> histActors(MAX_HISTORY);
for (int i = 0; i < MAX_HISTORY; ++i) {
histActors[i] = vtkTextActor::New();
histActors[i]->SetInput("");
ConfigureTextActor(histActors[i], FONT_SIZE, 0.5, 0.5, 0.5);
// Position: bottom of history area = just above current-line area
int y = MARGIN_TOP + (MAX_HISTORY - 1 - i) * LINE_HEIGHT;
histActors[i]->SetDisplayPosition(MARGIN_LEFT, HEIGHT - y);
ren->AddActor2D(histActors[i]);
}
// Current (actively typing) line actor — two actors: prefix + body
vtkNew<vtkTextActor> curPfxActor; // "[User] " in dim colour
vtkNew<vtkTextActor> curBodyActor; // body text in vivid colour
vtkNew<vtkTextActor> cursorActor; // blinking block
ConfigureTextActor(curPfxActor, FONT_SIZE, COL_LABEL[0], COL_LABEL[1], COL_LABEL[2]);
ConfigureTextActor(curBodyActor, FONT_SIZE, 1, 1, 1); // will be overridden per turn
ConfigureTextActor(cursorActor, FONT_SIZE, 1, 1, 1);
cursorActor->SetInput("|");
int curLineY = HEIGHT - (MARGIN_TOP + MAX_HISTORY * LINE_HEIGHT);
// Place current line below history
curPfxActor->SetDisplayPosition(MARGIN_LEFT, curLineY);
ren->AddActor2D(curPfxActor);
ren->AddActor2D(curBodyActor);
ren->AddActor2D(cursorActor);
// -----------------------------------------------------------------------
// Progress bar
// -----------------------------------------------------------------------
vtkPolyData* barPD = nullptr;
vtkPoints* barPts = nullptr;
vtkActor2D* barActor = MakeProgressBar(barPD, barPts);
ren->AddActor2D(barActor);
// Progress label
vtkNew<vtkTextActor> progLabelActor;
progLabelActor->SetInput("");
ConfigureTextActor(progLabelActor, 13,
COL_LABEL[0], COL_LABEL[1], COL_LABEL[2]);
progLabelActor->SetDisplayPosition(MARGIN_LEFT, 12);
ren->AddActor2D(progLabelActor);
// -----------------------------------------------------------------------
// Caption actor — bottom-centre, Arial, white, initially hidden
// -----------------------------------------------------------------------
vtkNew<vtkTextActor> captionActor;
captionActor->SetInput("");
captionActor->GetTextProperty()->SetFontFamilyToArial();
captionActor->GetTextProperty()->SetFontSize(20);
captionActor->GetTextProperty()->SetColor(1.0, 1.0, 1.0);
captionActor->GetTextProperty()->SetBold(0);
captionActor->GetTextProperty()->SetItalic(1);
captionActor->GetTextProperty()->ShadowOn();
captionActor->GetTextProperty()->SetShadowOffset(1, -1);
captionActor->GetTextProperty()->SetJustificationToCentered();
captionActor->GetTextProperty()->SetVerticalJustificationToBottom();
// Position: centred horizontally, in the gap between typing line and progress label
captionActor->SetDisplayPosition(WIDTH / 2, 32);
ren->AddActor2D(captionActor);
// -----------------------------------------------------------------------
// FFMPEG writer
// -----------------------------------------------------------------------
vtkNew<vtkWindowToImageFilter> w2i;
w2i->SetInput(renWin);
w2i->SetScale(1);
w2i->ReadFrontBufferOff();
vtkNew<vtkFFMPEGWriter> writer;
writer->SetInputConnection(w2i->GetOutputPort());
writer->SetFileName(outputPath.c_str());
writer->SetRate(FPS);
writer->SetBitRate(4000);
writer->SetBitRateTolerance(400);
writer->Start();
// -----------------------------------------------------------------------
// Helper: measure pixel width of a string in the current font
// (approximate: Courier is monospace, so width ≈ chars × charWidth)
// At font size 17 in Courier, one character ≈ 10.2px wide.
// -----------------------------------------------------------------------
const double CHAR_PX = 10.2;
auto bodyX = [&](const std::string& pfx) -> int {
return MARGIN_LEFT + (int)(pfx.size() * CHAR_PX);
};
// -----------------------------------------------------------------------
// Render one frame
// -----------------------------------------------------------------------
int frameCount = 0;
auto renderFrame = [&]() {
renWin->Render();
w2i->Modified();
writer->Write();
++frameCount;
};
// -----------------------------------------------------------------------
// History ring — holds completed display lines (most-recent last)
// -----------------------------------------------------------------------
std::vector<DisplayLine> history;
auto refreshHistory = [&]() {
// slot 0 = bottom row (newest), slot MAX_HISTORY-1 = top row (oldest).
// Brightness fades linearly from 0.85 (slot 0) to 0.20 (slot MAX-1).
// Blank spacer lines show as empty strings.
// The [USER]/[ASSISTANT] prefix stays at a fixed bright level so the
// speaker is always identifiable even in dim history.
int n = (int)history.size();
for (int slot = 0; slot < MAX_HISTORY; ++slot) {
int idx = n - 1 - slot; // slot 0 → newest, slot MAX-1 → oldest
if (idx < 0) {
histActors[slot]->SetInput("");
continue;
}
const auto& hl = history[idx];
if (hl.isBlank) {
histActors[slot]->SetInput("");
continue;
}
// Graduated brightness: bright near bottom, dim near top
double bodyBright = 0.20 + 0.65 * (1.0 - (double)slot / (MAX_HISTORY - 1));
const double* col = (hl.speaker == Speaker::User) ? COL_USER : COL_CLAUDE;
if (hl.isFirst) {
// Prefix stays vivid; body fades
// We render prefix + body as one string but colour the whole line
// at body brightness — a compromise since VTK TextActor is single-colour.
// Use a slightly higher floor for the prefix line so the tag is readable.
double pfxBright = std::min(1.0, bodyBright + 0.25);
std::string txt = hl.prefix + hl.body;
histActors[slot]->SetInput(txt.c_str());
histActors[slot]->GetTextProperty()->SetColor(
col[0] * pfxBright, col[1] * pfxBright, col[2] * pfxBright);
} else {
std::string txt = std::string(hl.prefix.size(), ' ') + hl.body;
histActors[slot]->SetInput(txt.c_str());
histActors[slot]->GetTextProperty()->SetColor(
col[0] * bodyBright, col[1] * bodyBright, col[2] * bodyBright);
}
}
};
// -----------------------------------------------------------------------
// Main animation loop
// -----------------------------------------------------------------------
int dlIdx = 0; // index into displayLines
// Count total turns to feed into progress label
// (display line index → turn index: count isFirst lines)
int totalTurns = 0;
for (const auto& dl : displayLines) if (dl.isFirst) ++totalTurns;
int turnsSeen = 0;
for (int li = 0; li < totalLines; ++li) {
const DisplayLine& dl = displayLines[li];
// Caption update: swap the bottom caption, no history, no typewriter
if (dl.isCaption) {
captionActor->SetInput(dl.body.c_str());
renderFrame();
continue;
}
// Blank spacer: push to history silently with no typewriter frames
if (dl.isBlank) {
history.push_back(dl);
refreshHistory();
continue;
}
// Vivid colour for current speaker
double cr, cg, cb;
if (dl.speaker == Speaker::User) {
cr = COL_USER[0]; cg = COL_USER[1]; cb = COL_USER[2];
} else {
cr = COL_CLAUDE[0]; cg = COL_CLAUDE[1]; cb = COL_CLAUDE[2];
}
curBodyActor->GetTextProperty()->SetColor(cr, cg, cb);
cursorActor->GetTextProperty()->SetColor(cr, cg, cb);
// Prefix: vivid speaker colour for first line, indent for continuation
if (dl.isFirst) {
curPfxActor->SetInput(dl.prefix.c_str());
curPfxActor->GetTextProperty()->SetColor(cr, cg, cb);
} else {
curPfxActor->SetInput(std::string(dl.prefix.size(), ' ').c_str());
curPfxActor->GetTextProperty()->SetColor(0, 0, 0);
}
int bx = bodyX(dl.prefix);
curBodyActor->SetDisplayPosition(bx, curLineY);
cursorActor->SetDisplayPosition(bx, curLineY);
if (dl.isFirst) ++turnsSeen;
// Update progress label
{
char buf[64];
std::snprintf(buf, sizeof(buf), "Turn %d / %d", turnsSeen, totalTurns);
progLabelActor->SetInput(buf);
}
// Update progress bar width
{
double frac = (totalFrames > 0) ? (double)frameCount / totalFrames : 0.0;
double barW = frac * (WIDTH - 2 * MARGIN_LEFT);
barPts->SetPoint(0, MARGIN_LEFT, HEIGHT - 8, 0);
barPts->SetPoint(1, MARGIN_LEFT, HEIGHT - 2, 0);
barPts->SetPoint(2, MARGIN_LEFT + barW, HEIGHT - 2, 0);
barPts->SetPoint(3, MARGIN_LEFT + barW, HEIGHT - 8, 0);
barPts->Modified();
barPD->Modified();
}
// Typewriter: reveal CHARS_PER_FRAME characters per rendered frame.
// Always render the fully-complete state last.
const std::string& body = dl.body;
int ci = 0;
while (true) {
curBodyActor->SetInput(body.substr(0, ci).c_str());
double cx = bx + ci * CHAR_PX;
cursorActor->SetDisplayPosition((int)cx, curLineY);
cursorActor->SetVisibility((frameCount % 2 == 0) ? 1 : 0);
renderFrame();
if (ci >= (int)body.size()) break;
ci = std::min(ci + CHARS_PER_FRAME, (int)body.size());
}
// Line is complete — move it to history
history.push_back(dl);
refreshHistory();
// Clear current line display
curPfxActor->SetInput("");
curBodyActor->SetInput("");
cursorActor->SetVisibility(0);
// Pause frames
for (int p = 0; p < PAUSE_FRAMES; ++p) {
// Update progress bar
{
double frac = (totalFrames > 0) ? (double)frameCount / totalFrames : 0.0;
double barW = frac * (WIDTH - 2 * MARGIN_LEFT);
barPts->SetPoint(0, MARGIN_LEFT, HEIGHT - 8, 0);
barPts->SetPoint(1, MARGIN_LEFT, HEIGHT - 2, 0);
barPts->SetPoint(2, MARGIN_LEFT + barW, HEIGHT - 2, 0);
barPts->SetPoint(3, MARGIN_LEFT + barW, HEIGHT - 8, 0);
barPts->Modified();
barPD->Modified();
}
renderFrame();
}
}
// -----------------------------------------------------------------------
// Finish
// -----------------------------------------------------------------------
writer->End();
// Tidy up manual ref-counted objects
barActor->Delete();
barPD->Delete();
barPts->Delete();
for (auto* a : histActors) a->Delete();
std::cerr << "Wrote " << frameCount << " frames ("
<< frameCount / FPS << " s) to " << outputPath << "\n";
return 0;
}
-729
View File
@@ -1,729 +0,0 @@
// Derived from VTK/Examples/Cxx/Medical2.cxx
// The example reads a volume dataset, extracts two isosurfaces that
// represent the skin and bone, and then displays them.
//
// Modified heavily by Peter Boyle to display lattice field theory data as movies and compare multiple files
#include <vtkActor.h>
#include <vtkCamera.h>
#include <vtkMetaImageReader.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkOutlineFilter.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkStripper.h>
#include <vtkImageData.h>
#include <vtkVersion.h>
#include <vtkCallbackCommand.h>
#include <vtkTextActor.h>
#include <vtkTextProperty.h>
#define MPEG
#ifdef MPEG
#include <vtkFFMPEGWriter.h>
#endif
#include <vtkProperty2D.h>
#include <vtkSliderWidget.h>
#include <vtkSliderRepresentation2D.h>
#include <vtkWindowToImageFilter.h>
#include <array>
#include <string>
#include <Grid/Grid.h>
#define USE_FLYING_EDGES
#ifdef USE_FLYING_EDGES
#include <vtkFlyingEdges3D.h>
typedef vtkFlyingEdges3D isosurface;
#else
#include <vtkMarchingCubes.h>
typedef vtkMarchingCubes isosurface;
#endif
int mpeg = 0 ;
int framerate = 10;
template <class T> void readFile(T& out, std::string const fname){
Grid::emptyUserRecord record;
Grid::ScidacReader RD;
RD.open(fname);
RD.readScidacFieldRecord(out,record);
RD.close();
}
using namespace Grid;
class FrameUpdater : public vtkCallbackCommand
{
public:
FrameUpdater() {
ffile=0;
TimerCount = 0;
xoff = 0;
t = 0;
imageData = nullptr;
timerId = 0;
maxCount = -1;
old_file=-1;
}
static FrameUpdater* New()
{
FrameUpdater* cb = new FrameUpdater;
cb->TimerCount = 0;
return cb;
}
//
// Must map a x,y,z + frame index into
// i) a d-dimensional site Coordinate
// ii) a file name
// Need a:
// loop_ranges
// sum_ranges
// loop_vol -- map loop_idx -> loop_coor
// sum_vol -- map sum_idx -> sum_coor with Lexicographic
//
/*
* Just set this up
*/
int old_file ; // Cache, avoid reread
Coordinate latt;
Coordinate xyz_dims ; // List lattice dimensions corresponding to xyz_dims displayed
Coordinate xyz_ranges ; // 3-vector
Coordinate g_xyz_ranges; // Nd-vector
uint64_t xyz_vol ;
Coordinate loop_dims; // List lattice dimensions put into movie time
Coordinate loop_ranges; // movie time ranges
uint64_t loop_vol;
Coordinate sum_dims; // List lattice dimensions summed
Coordinate sum_ranges; // summation ranges
uint64_t sum_vol;
Coordinate slice_dims; // List slice dimensions
Coordinate Slice;
std::vector<std::string> files; // file list that is looped over
int Nd;
GridBase *grid;
Grid::LatticeComplexD *grid_data;
void SetGrid(GridBase *_grid)
{
grid = _grid;
Nd=grid->Nd();
latt = grid->GlobalDimensions();
grid_data = new Grid::LatticeComplexD(grid);
}
void SetFiles(std::vector<std::string> list) { files = list; old_file = -1; }
void SetSlice(Coordinate _Slice) { Slice = _Slice;} // Offset / skew for lattice coords
void SetSumDimensions (Coordinate _SumDims ) {
sum_ranges=Coordinate(Nd);
sum_dims = _SumDims; // 1 hot for dimensions summed
sum_vol = 1;
for(int d=0;d<sum_dims.size();d++){
if ( sum_dims[d] == 1 ) sum_ranges[d] = latt[d];
else sum_ranges[d] = 1;
sum_vol*=sum_ranges[d];
}
}
void SetLoopDimensions(Coordinate _LoopDims) {
loop_ranges=Coordinate(Nd);
loop_dims= _LoopDims;
loop_vol = 1;
for(int d=0;d<loop_dims.size();d++){
if ( loop_dims[d] == 1 ) loop_ranges[d] = latt[d];
else loop_ranges[d] = 1;
loop_vol*=loop_ranges[d];
}
} //
void SetDisplayDimensions(Coordinate _xyz_dims ) {
g_xyz_ranges=Coordinate(Nd);
xyz_ranges=Coordinate(3);
xyz_dims = _xyz_dims;
xyz_vol = 1;
for(int d=0;d<3;d++){
xyz_ranges[d] = latt[xyz_dims[d]];
xyz_vol *= xyz_ranges[d];
}
// Find dim extents for grid
int dd=0;
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(void) {
Coordinate _slice_dims;
for ( int d=0;d<Nd;d++){
int is_slice = 1;
if(g_xyz_ranges[d]>1) is_slice = 0;
if(loop_dims[d]) is_slice = 0;
if(sum_dims[d] ) is_slice = 0;
if(is_slice) _slice_dims.push_back(d);
}
slice_dims = _slice_dims;
std::cout << " Setting Slice Dimensions to "<<slice_dims<<std::endl;
}
virtual void Execute(vtkObject* caller, unsigned long eventId,void* vtkNotUsed(callData))
{
const int max=256;
char text_string[max];
auto latt_size = grid->GlobalDimensions();
if ( vtkCommand::KeyPressEvent == eventId ) {
vtkRenderWindowInteractor* iren = static_cast<vtkRenderWindowInteractor*>(caller);
std::string key = iren->GetKeySym();
std::cout << "Pressed: " << key << std::endl;
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) % files.size();
}
if ( key == "less" ) {
ffile = (ffile - 1 + files.size()) % files.size();
}
std::cout <<"Slice " <<Slice <<std::endl;
std::cout <<"File " <<ffile <<std::endl;
}
// Make a new frame for frame index TimerCount
if ( vtkCommand::TimerEvent == eventId || vtkCommand::KeyPressEvent == eventId)
{
int file = ((this->TimerCount / loop_vol) + ffile )%files.size();
if ( file != old_file ) {
readFile(*grid_data,files[file]);
old_file = file;
}
RealD max, min, max_abs,min_abs;
Coordinate max_site;
Coordinate min_site;
Coordinate max_abs_site;
Coordinate min_abs_site;
for(int idx=0;idx<grid->gSites();idx++){
Coordinate site;
Lexicographic::CoorFromIndex (site,idx,latt);
RealD val=real(peekSite(*grid_data,site));
if (idx==0){
max = min = val;
max_abs = min_abs = fabs(val);
max_site = site;
min_site = site;
min_abs_site = site;
max_abs_site = site;
} else {
if ( val > max ) {
max=val;
max_site = site;
}
if ( fabs(val) > max_abs ) {
max_abs=fabs(val);
max_abs_site = site;
}
if ( val < min ) {
min=val;
min_site = site;
}
if ( fabs(val) < min_abs ) {
min_abs=fabs(val);
min_abs_site = site;
}
}
}
std::cout << " abs_max "<<max_abs<<" at " << max_abs_site<<std::endl;
std::cout << " abs_min "<<min_abs<<" at " << min_abs_site<<std::endl;
std::cout << " max "<<max<<" at " << max_site<<std::endl;
std::cout << " min "<<min<<" at " << min_site<<std::endl;
// Looped dimensions, map index to coordinate
int loop_idx = this->TimerCount % loop_vol;
Coordinate loop_coor;
Lexicographic::CoorFromIndex (loop_coor,loop_idx,loop_ranges);
// Loop over xyz sites
Coordinate xyz_coor(3);
Coordinate g_xyz_coor(Nd);
Coordinate 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 sum_value = 0.0;
for(uint64_t sum_idx = 0 ; sum_idx< sum_vol; sum_idx++){
Lexicographic::CoorFromIndex (sum_coor,sum_idx,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];
}
sum_value+= real(peekSite(*grid_data,site));
if(xyz==0) std::cout << "sum "<<sum_idx<<" "<<sum_value<<std::endl;
}
imageData->SetScalarComponentFromDouble(xyz_coor[0],xyz_coor[1],xyz_coor[2],0,sum_value);
}
imageData->Modified();
std::stringstream ss;
ss<< files[file] <<"\nSlice "<<Slice << "\nLoop " <<loop_coor<<"\nSummed "<<sum_dims;
text->SetInput(ss.str().c_str());
vtkRenderWindowInteractor* iren = dynamic_cast<vtkRenderWindowInteractor*>(caller);
iren->GetRenderWindow()->Render();
}
if ( vtkCommand::TimerEvent == eventId ) {
++this->TimerCount;
std::cout << " This was a timer event count "<<this->TimerCount << std::endl;
}
if (this->TimerCount >= this->maxCount) {
vtkRenderWindowInteractor* iren = dynamic_cast<vtkRenderWindowInteractor*>(caller);
if (this->timerId > -1)
{
iren->DestroyTimer(this->timerId);
}
}
}
private:
int TimerCount;
int ffile;
int xoff;
int t;
public:
vtkImageData* imageData = nullptr;
vtkTextActor* text = nullptr;
vtkFFMPEGWriter *writer = nullptr;
int timerId ;
int maxCount ;
double rms;
isosurface * posExtractor;
isosurface * negExtractor;
};
class SliderCallback : public vtkCommand
{
public:
static SliderCallback* New()
{
return new SliderCallback;
}
virtual void Execute(vtkObject* caller, unsigned long eventId, void* callData)
{
vtkSliderWidget *sliderWidget = vtkSliderWidget::SafeDownCast(caller);
if (sliderWidget)
{
contour = ((vtkSliderRepresentation *)sliderWidget->GetRepresentation())->GetValue();
}
fu->posExtractor->SetValue(0, SliderCallback::contour*fu->rms);
fu->negExtractor->SetValue(0, -SliderCallback::contour*fu->rms);
fu->posExtractor->Modified();
fu->negExtractor->Modified();
}
public:
static double contour;
FrameUpdater * fu;
};
FrameUpdater * KBfu;
void KeypressCallbackFunction(vtkObject* caller, long unsigned int eventId,
void* clientData, void* callData)
{
std::cout << "Keypress callback" << std::endl;
vtkRenderWindowInteractor* iren = static_cast<vtkRenderWindowInteractor*>(caller);
std::cout << "Pressed: " << iren->GetKeySym() << std::endl;
// imageData->Modified();
}
double SliderCallback::contour;
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::cout << argc << " command Line arguments "<<std::endl;
for(int c=0;c<argc;c++) {
std::cout << " - "<<argv[c]<<std::endl;
}
std::vector<std::string> 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") ){
mpeg = 1;
}
#endif
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);
}
for(int c=0;c<file_list.size();c++) {
std::cout << " file: "<<file_list[c]<<std::endl;
}
int NoTime = 0;
int Nd; Nd = Grid.Nd();
Coordinate Slice(Nd,0);
Coordinate SumDims(Nd,0);
Coordinate LoopDims(Nd,0);
Coordinate 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);
std::cout << "xyz : "<<XYZDims<<std::endl;
}
if( GridCmdOptionExists(argv,argv+argc,"--notime") ){
NoTime = 1;
std::cout << "Suppressing time loop"<<std::endl;
}
// Common things:
vtkNew<vtkNamedColors> colors;
std::array<unsigned char, 4> posColor{{240, 184, 160, 255}}; colors->SetColor("posColor", posColor.data());
std::array<unsigned char, 4> bkg{{51, 77, 102, 255}}; colors->SetColor("BkgColor", bkg.data());
// Create the renderer, the render window, and the interactor. The renderer
// draws into the render window, the interactor enables mouse- and
// keyboard-based interaction with the data within the render window.
//
vtkNew<vtkRenderWindow> renWin;
vtkNew<vtkRenderWindowInteractor> iren;
iren->SetRenderWindow(renWin);
// std::vector<LatticeComplexD> data(file_list.size(),&Grid);
// FieldMetaData header;
int frameCount = file_list.size();
for(int d=0;d<Grid.Nd();d++) {
if ( LoopDims[d] ) frameCount*= latt_size[d];
}
// It is convenient to create an initial view of the data. The FocalPoint
// and Position form a vector direction. Later on (ResetCamera() method)
// this vector is used to position the camera to look at the data in
// this direction.
vtkNew<vtkCamera> 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<vtkRenderer> aRenderer;
renWin->AddRenderer(aRenderer);
double vol = Grid.gSites();
std::cout << "Reading "<<file_list[0]<<std::endl;
double nrm, nrmbar,rms, contour;
{
LatticeComplexD data(&Grid);
readFile(data,file_list[0]);
nrm = norm2(data);
}
nrmbar = nrm/vol;
rms = sqrt(nrmbar);
contour = default_contour * rms; // default to 1 x RMS
// The following reader is used to read a series of 2D slices (images)
// that compose the volume. The slice dimensions are set, and the
// pixel spacing. The data Endianness must also be specified. The reader
// uses the FilePrefix in combination with the slice number to construct
// filenames using the format FilePrefix.%d. (In this case the FilePrefix
// is the root name of the file: quarter.)
vtkNew<vtkImageData> imageData;
imageData->SetDimensions(latt_size[0],latt_size[1],latt_size[2]);
imageData->AllocateScalars(VTK_DOUBLE, 1);
for(int xx=0;xx<latt_size[0];xx++){
for(int yy=0;yy<latt_size[1];yy++){
for(int zz=0;zz<latt_size[2];zz++){
Coordinate site({xx,yy,zz,0});
RealD value = 0;
imageData->SetScalarComponentFromDouble(xx,yy,zz,0,value);
}}}
vtkNew<isosurface> posExtractor;
posExtractor->SetInputData(imageData);
posExtractor->SetValue(0, contour);
vtkNew<vtkStripper> posStripper;
posStripper->SetInputConnection(posExtractor->GetOutputPort());
vtkNew<vtkPolyDataMapper> posMapper;
posMapper->SetInputConnection(posStripper->GetOutputPort());
posMapper->ScalarVisibilityOff();
vtkNew<vtkActor> 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);
// An isosurface, or contour value is set
// The triangle stripper is used to create triangle strips from the
// isosurface; these render much faster on may systems.
vtkNew<isosurface> negExtractor;
negExtractor->SetInputData(imageData);
negExtractor->SetValue(0, -contour);
vtkNew<vtkStripper> negStripper;
negStripper->SetInputConnection(negExtractor->GetOutputPort());
vtkNew<vtkPolyDataMapper> negMapper;
negMapper->SetInputConnection(negStripper->GetOutputPort());
negMapper->ScalarVisibilityOff();
vtkNew<vtkActor> neg;
neg->SetMapper(negMapper);
neg->GetProperty()->SetDiffuseColor(colors->GetColor3d("Ivory").GetData());
// An outline provides context around the data.
vtkNew<vtkOutlineFilter> outlineData;
outlineData->SetInputData(imageData);
vtkNew<vtkPolyDataMapper> mapOutline;
mapOutline->SetInputConnection(outlineData->GetOutputPort());
vtkNew<vtkActor> outline;
outline->SetMapper(mapOutline);
outline->GetProperty()->SetColor(colors->GetColor3d("Black").GetData());
vtkNew<vtkTextActor> Text;
// Text->SetInput(file_list[f].c_str());
Text->SetPosition2(0,0);
Text->GetTextProperty()->SetFontSize(24);
Text->GetTextProperty()->SetColor(colors->GetColor3d("Gold").GetData());
vtkNew<vtkTextActor> TextT;
TextT->SetInput("T=0");
TextT->SetPosition(0,.7*1025);
TextT->GetTextProperty()->SetFontSize(24);
TextT->GetTextProperty()->SetColor(colors->GetColor3d("Gold").GetData());
// Actors are added to the renderer. An initial camera view is created.
// The Dolly() method moves the camera towards the FocalPoint,
// thereby enlarging the image.
// aRenderer->AddActor(Text);
aRenderer->AddActor(TextT);
aRenderer->AddActor(outline);
aRenderer->AddActor(pos);
aRenderer->AddActor(neg);
// Sign up to receive TimerEvent
vtkNew<FrameUpdater> 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->grid_data = &data[f];
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);
// double nf = file_list.size();
// std::cout << " Adding renderer " <<f<<" of "<<nf<<std::endl;
aRenderer->SetViewport(0.0, 0.0,1.0 , 1.0);
// Note that when camera movement occurs (as it does in the Dolly()
// method), the clipping planes often need adjusting. Clipping planes
// consist of two planes: near and far along the view direction. The
// near plane clips out objects in front of the plane; the far plane
// clips out objects behind the plane. This way only what is drawn
// between the planes is actually rendered.
aRenderer->ResetCameraClippingRange();
// Set a background color for the renderer and set the size of the
// render window (expressed in pixels).
// Initialize the event loop and then start it.
renWin->SetSize(1024, 1024);
renWin->SetWindowName("FieldDensity");
renWin->Render();
// Take a pointer to the FrameUpdater for keypress mgt.
// KBfu = fu;
// vtkNew<vtkCallbackCommand> keypressCallback;
// keypressCallback->SetCallback(KeypressCallbackFunction);
// iren->AddObserver(vtkCommand::KeyPressEvent,keypressCallback);
iren->Initialize();
if ( mpeg ) {
#ifdef MPEG
vtkWindowToImageFilter *imageFilter = vtkWindowToImageFilter::New();
imageFilter->SetInput( renWin );
imageFilter->SetInputBufferTypeToRGB();
vtkFFMPEGWriter *writer = vtkFFMPEGWriter::New();
writer->SetFileName("movie.avi");
writer->SetRate(framerate);
writer->SetInputConnection(imageFilter->GetOutputPort());
writer->Start();
for(int i=0;i<fu->maxCount;i++){
fu->Execute(iren,vtkCommand::TimerEvent,nullptr);
imageFilter->Modified();
writer->Write();
}
writer->End();
writer->Delete();
#else
assert(-1 && "MPEG support not compiled");
#endif
} else {
// Add control of contour threshold
// Create a slider widget
vtkSmartPointer<vtkSliderRepresentation2D> sliderRep = vtkSmartPointer<vtkSliderRepresentation2D>::New();
sliderRep->SetMinimumValue(0.0);
sliderRep->SetMaximumValue(10.0);
sliderRep->SetValue(1.0);
sliderRep->SetTitleText("Fraction RMS");
// Set color properties:
// Change the color of the knob that slides
// sliderRep->GetSliderProperty()->SetColor(colors->GetColor3d("Green").GetData());
sliderRep->GetTitleProperty()->SetColor(colors->GetColor3d("AliceBlue").GetData());
sliderRep->GetLabelProperty()->SetColor(colors->GetColor3d("AliceBlue").GetData());
sliderRep->GetSelectedProperty()->SetColor(colors->GetColor3d("DeepPink").GetData());
// Change the color of the bar
sliderRep->GetTubeProperty()->SetColor(colors->GetColor3d("MistyRose").GetData());
sliderRep->GetCapProperty()->SetColor(colors->GetColor3d("Yellow").GetData());
sliderRep->SetSliderLength(0.05);
sliderRep->SetSliderWidth(0.025);
sliderRep->SetEndCapLength(0.02);
sliderRep->GetPoint1Coordinate()->SetCoordinateSystemToNormalizedDisplay();
sliderRep->GetPoint1Coordinate()->SetValue(0.1, 0.1);
sliderRep->GetPoint2Coordinate()->SetCoordinateSystemToNormalizedDisplay();
sliderRep->GetPoint2Coordinate()->SetValue(0.9, 0.1);
vtkSmartPointer<vtkSliderWidget> sliderWidget = vtkSmartPointer<vtkSliderWidget>::New();
sliderWidget->SetInteractor(iren);
sliderWidget->SetRepresentation(sliderRep);
sliderWidget->SetAnimationModeToAnimate();
sliderWidget->EnabledOn();
// Create the slider callback
vtkSmartPointer<SliderCallback> slidercallback = vtkSmartPointer<SliderCallback>::New();
slidercallback->fu = fu;
sliderWidget->AddObserver(vtkCommand::InteractionEvent, slidercallback);
if ( NoTime==0 ) {
int timerId = iren->CreateRepeatingTimer(10000/framerate);
std::cout << "timerId "<<timerId<<std::endl;
}
// Start the interaction and timer
iren->Start();
}
Grid_finalize();
return EXIT_SUCCESS;
}