From b07ccf3d5e7801490526d82284c61827a005961f Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Fri, 2 Oct 2026 12:02:54 -0400
Subject: [PATCH] Force test reorg

---
 tests/forces/Test_bdy.cc                  | 84 +----------------------
 tests/forces/Test_double_ratio.cc         | 81 ++--------------------
 tests/forces/Test_dwf_ratio_4dpf_force.cc | 66 +-----------------
 tests/forces/Test_dwf_ratio_leftprec.cc   | 64 +----------------
 4 files changed, 11 insertions(+), 284 deletions(-)

diff --git a/tests/forces/Test_bdy.cc b/tests/forces/Test_bdy.cc
index d328b49fb..106a73f31 100644
--- a/tests/forces/Test_bdy.cc
+++ b/tests/forces/Test_bdy.cc
@@ -99,91 +99,12 @@ Author: Peter Boyle <pboyle@bnl.gov>
 using namespace std;
 using namespace Grid;
 
+#include "ForceTest.h"
+
 typedef MobiusFermionD FermionAction;
 typedef WilsonImplD FimplD;
 typedef WilsonImplD FermionImplPolicy;
 
-template<class Gimpl>
-void ForceTest(Action<LatticeGaugeField> &action,LatticeGaugeField & U,MomentumFilterBase<LatticeGaugeField> &Filter)
-{
-  GridBase *UGrid = U.Grid();
-
-  std::vector<int> seeds({1,2,3,5});
-  GridSerialRNG            sRNG;         sRNG.SeedFixedIntegers(seeds);
-  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds);
-
-  LatticeColourMatrix Pmu(UGrid); 
-  LatticeGaugeField P(UGrid); 
-  LatticeGaugeField UdSdU(UGrid); 
-
-  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
-  std::cout << GridLogMessage << " Force test for "<<action.action_name()<<std::endl;
-  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
-  
-  RealD eps=0.005;
-
-  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
-  std::cout << GridLogMessage << " Refresh "<<action.action_name()<<std::endl;
-  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
-  
-  Gimpl::generate_momenta(P,sRNG,RNG4);
-  Filter.applyFilter(P);
-
-#if  0
-  FieldMetaData header;
-  std::string file("./ckpoint_lat.2000");
-  NerscIO::readConfiguration(U,header,file);
-#else
-  U = 1.0;
-#endif
-  action.refresh(U,sRNG,RNG4);
-
-  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
-  std::cout << GridLogMessage << " Action "<<action.action_name()<<std::endl;
-  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
-
-  RealD S1 = action.S(U);
-
-  Gimpl::update_field(P,U,eps);
-
-  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
-  std::cout << GridLogMessage << " Derivative "<<action.action_name()<<std::endl;
-  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
-  action.deriv(U,UdSdU);
-  UdSdU = Ta(UdSdU);
-  Filter.applyFilter(UdSdU);
-
-  DumpSliceNorm("Force",UdSdU,Nd-1);
-  
-  Gimpl::update_field(P,U,eps);
-  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
-  std::cout << GridLogMessage << " Action "<<action.action_name()<<std::endl;
-  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
-  
-  RealD S2 = action.S(U);
-
-  // Use the derivative
-  LatticeComplex dS(UGrid); dS = Zero();
-  for(int mu=0;mu<Nd;mu++){
-    auto UdSdUmu = PeekIndex<LorentzIndex>(UdSdU,mu);
-    Pmu= PeekIndex<LorentzIndex>(P,mu);
-    dS = dS - trace(Pmu*UdSdUmu)*eps*2.0*2.0;
-  }
-  ComplexD dSpred    = sum(dS);
-  RealD diff =  S2-S1-dSpred.real();
-
-  std::cout<< GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
-  std::cout<< GridLogMessage << "S1 : "<< S1    <<std::endl;
-  std::cout<< GridLogMessage << "S2 : "<< S2    <<std::endl;
-  std::cout<< GridLogMessage << "dS : "<< S2-S1 <<std::endl;
-  std::cout<< GridLogMessage << "dSpred : "<< dSpred.real() <<std::endl;
-  std::cout<< GridLogMessage << "diff : "<< diff<<std::endl;
-  std::cout<< GridLogMessage << "*********************************************************"<<std::endl;
-  //  GRID_ASSERT(diff<1.0);
-  std::cout<< GridLogMessage << "Done" <<std::endl;
-  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
-}
-
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
@@ -230,6 +151,7 @@ int main (int argc, char ** argv)
 
   ///////////////////// Gauge Field and Gauge Forces ////////////////////////////
   LatticeGaugeField U(UGrid);
+  U = 1.0;
 
   RealD beta=6.0;
   WilsonGaugeActionR PlaqAction(beta);
diff --git a/tests/forces/Test_double_ratio.cc b/tests/forces/Test_double_ratio.cc
index 418ba87b5..1306ee7d2 100644
--- a/tests/forces/Test_double_ratio.cc
+++ b/tests/forces/Test_double_ratio.cc
@@ -303,83 +303,7 @@ typedef MobiusFermionD FermionAction;
 typedef WilsonImplD FimplD;
 typedef WilsonImplD FermionImplPolicy;
 
-template<class Gimpl>
-void ForceTest(Action<LatticeGaugeField> &action,LatticeGaugeField & U,MomentumFilterBase<LatticeGaugeField> &Filter)
-{
-  GridBase *UGrid = U.Grid();
-
-  std::vector<int> seeds({1,2,3,5});
-  GridSerialRNG            sRNG;         sRNG.SeedFixedIntegers(seeds);
-  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds);
-
-  LatticeColourMatrix Pmu(UGrid); 
-  LatticeGaugeField P(UGrid); 
-  LatticeGaugeField UdSdU(UGrid); 
-
-  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
-  std::cout << GridLogMessage << " Force test for "<<action.action_name()<<std::endl;
-  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
-  
-  RealD eps=0.005;
-
-  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
-  std::cout << GridLogMessage << " Refresh "<<action.action_name()<<std::endl;
-  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
-  
-  Gimpl::generate_momenta(P,sRNG,RNG4);
-  Filter.applyFilter(P);
-
-  FieldMetaData header;
-  std::string file("./ckpoint_lat.2000");
-  NerscIO::readConfiguration(U,header,file);
- 
-  action.refresh(U,sRNG,RNG4);
-
-  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
-  std::cout << GridLogMessage << " Action "<<action.action_name()<<std::endl;
-  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
-
-  RealD S1 = action.S(U);
-
-  Gimpl::update_field(P,U,eps);
-
-  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
-  std::cout << GridLogMessage << " Derivative "<<action.action_name()<<std::endl;
-  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
-  action.deriv(U,UdSdU);
-  UdSdU = Ta(UdSdU);
-  Filter.applyFilter(UdSdU);
-
-  DumpSliceNorm("Force",UdSdU,Nd-1);
-  
-  Gimpl::update_field(P,U,eps);
-  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
-  std::cout << GridLogMessage << " Action "<<action.action_name()<<std::endl;
-  std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
-  
-  RealD S2 = action.S(U);
-
-  // Use the derivative
-  LatticeComplex dS(UGrid); dS = Zero();
-  for(int mu=0;mu<Nd;mu++){
-    auto UdSdUmu = PeekIndex<LorentzIndex>(UdSdU,mu);
-    Pmu= PeekIndex<LorentzIndex>(P,mu);
-    dS = dS - trace(Pmu*UdSdUmu)*eps*2.0*2.0;
-  }
-  ComplexD dSpred    = sum(dS);
-  RealD diff =  S2-S1-dSpred.real();
-
-  std::cout<< GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
-  std::cout<< GridLogMessage << "S1 : "<< S1    <<std::endl;
-  std::cout<< GridLogMessage << "S2 : "<< S2    <<std::endl;
-  std::cout<< GridLogMessage << "dS : "<< S2-S1 <<std::endl;
-  std::cout<< GridLogMessage << "dSpred : "<< dSpred.real() <<std::endl;
-  std::cout<< GridLogMessage << "diff : "<< diff<<std::endl;
-  std::cout<< GridLogMessage << "*********************************************************"<<std::endl;
-  //  GRID_ASSERT(diff<1.0);
-  std::cout<< GridLogMessage << "Done" <<std::endl;
-  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
-}
+#include "ForceTest.h"
 
 int main (int argc, char ** argv)
 {
@@ -427,6 +351,9 @@ int main (int argc, char ** argv)
 
   ///////////////////// Gauge Field and Gauge Forces ////////////////////////////
   LatticeGaugeField U(UGrid);
+  FieldMetaData header;
+  std::string file("./ckpoint_lat.2000");
+  NerscIO::readConfiguration(U,header,file);
 
   RealD beta=6.0;
   WilsonGaugeActionR PlaqAction(beta);
diff --git a/tests/forces/Test_dwf_ratio_4dpf_force.cc b/tests/forces/Test_dwf_ratio_4dpf_force.cc
index 41f1f6290..c0411f705 100644
--- a/tests/forces/Test_dwf_ratio_4dpf_force.cc
+++ b/tests/forces/Test_dwf_ratio_4dpf_force.cc
@@ -35,7 +35,7 @@ Author: Peter Boyle <pboyle@bnl.gov>
 //     composition identity [P M^-1 V Pdag][P V^-1 M Pdag] = 1 for the chosen
 //     (P,Pdag) wall pair.  NO PREDICTION is made about which convention
 //     passes -- that is what the test decides.
-//  T2 ForceTest (idiom from Test_double_ratio.cc) : midpoint-derivative
+//  T2 ForceTest (tests/forces/ForceTest.h) : midpoint-derivative
 //     check of deriv against S.  Should PASS for BOTH conventions (S and
 //     deriv use the same literal-adjoint pair by construction).
 //  T3 Trivial-ratio control (V == M) : T1 with NumOp = DenOp.  The solve
@@ -124,69 +124,7 @@ RealD HeatbathIdentityTest(TwoFlavourRatio4DPseudoFermionAction<Impl> &action,
   return defect;
 }
 
-////////////////////////////////////////////////////////////////////
-// T2 : ForceTest idiom from Test_double_ratio.cc (midpoint derivative)
-////////////////////////////////////////////////////////////////////
-template<class Gimpl>
-void ForceTest(Action<LatticeGaugeField> &action,LatticeGaugeField & U,MomentumFilterBase<LatticeGaugeField> &Filter)
-{
-  GridBase *UGrid = U.Grid();
-
-  std::vector<int> seeds({1,2,3,5});
-  GridSerialRNG            sRNG;         sRNG.SeedFixedIntegers(seeds);
-  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds);
-
-  LatticeColourMatrix Pmu(UGrid);
-  LatticeGaugeField P(UGrid);
-  LatticeGaugeField UdSdU(UGrid);
-
-  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
-  std::cout << GridLogMessage << " Force test for "<<action.action_name()<<std::endl;
-  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
-
-  RealD eps=0.005;
-
-  Gimpl::generate_momenta(P,sRNG,RNG4);
-  Filter.applyFilter(P);
-
-  action.refresh(U,sRNG,RNG4);
-
-  RealD S1 = action.S(U);
-
-  Gimpl::update_field(P,U,eps);
-
-  action.deriv(U,UdSdU);
-  UdSdU = Ta(UdSdU);
-  Filter.applyFilter(UdSdU);
-
-  DumpSliceNorm("Force",UdSdU,Nd-1);
-
-  Gimpl::update_field(P,U,eps);
-
-  RealD S2 = action.S(U);
-
-  // Use the derivative
-  LatticeComplex dS(UGrid); dS = Zero();
-  for(int mu=0;mu<Nd;mu++){
-    auto UdSdUmu = PeekIndex<LorentzIndex>(UdSdU,mu);
-    Pmu= PeekIndex<LorentzIndex>(P,mu);
-    dS = dS - trace(Pmu*UdSdUmu)*eps*2.0*2.0;
-  }
-  ComplexD dSpred    = sum(dS);
-  RealD diff =  S2-S1-dSpred.real();
-
-  std::cout<< GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
-  std::cout<< GridLogMessage << "S1 : "<< S1    <<std::endl;
-  std::cout<< GridLogMessage << "S2 : "<< S2    <<std::endl;
-  std::cout<< GridLogMessage << "dS : "<< S2-S1 <<std::endl;
-  std::cout<< GridLogMessage << "dSpred : "<< dSpred.real() <<std::endl;
-  std::cout<< GridLogMessage << "diff : "<< diff<<std::endl;
-  std::cout<< GridLogMessage << "diff/dS : "<< diff/(S2-S1)<<std::endl;
-  std::cout<< GridLogMessage << "*********************************************************"<<std::endl;
-  //  GRID_ASSERT(diff<1.0);
-  std::cout<< GridLogMessage << "Done" <<std::endl;
-  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
-}
+#include "ForceTest.h"
 
 int main (int argc, char ** argv)
 {
diff --git a/tests/forces/Test_dwf_ratio_leftprec.cc b/tests/forces/Test_dwf_ratio_leftprec.cc
index 113e262b3..659012364 100644
--- a/tests/forces/Test_dwf_ratio_leftprec.cc
+++ b/tests/forces/Test_dwf_ratio_leftprec.cc
@@ -38,7 +38,7 @@ Author: Peter Boyle <pboyle@bnl.gov>
 //             BOTH classes (E0a also validates the twin-eta capture).
 //   E1      : S_classic == S_leftprec       (relative, ~1e-8)
 //   E2      : deriv_classic == deriv_leftprec (pointwise field norm, ~1e-8)
-//   F1      : ForceTest (Test_double_ratio.cc idiom) on the LeftPrec class.
+//   F1      : ForceTest (tests/forces/ForceTest.h) on the LeftPrec class.
 //
 // All asserts are hard: this is the regression gate for the new class.
 // Run small, e.g.:  ./Test_dwf_ratio_leftprec --grid 8.8.8.8
@@ -138,67 +138,7 @@ public:
   }
 };
 
-////////////////////////////////////////////////////////////////////
-// ForceTest idiom from Test_double_ratio.cc (midpoint derivative)
-////////////////////////////////////////////////////////////////////
-template<class Gimpl>
-void ForceTest(Action<LatticeGaugeField> &action,LatticeGaugeField & U,MomentumFilterBase<LatticeGaugeField> &Filter)
-{
-  GridBase *UGrid = U.Grid();
-
-  std::vector<int> seeds({1,2,3,5});
-  GridSerialRNG            sRNG;         sRNG.SeedFixedIntegers(seeds);
-  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds);
-
-  LatticeColourMatrix Pmu(UGrid);
-  LatticeGaugeField P(UGrid);
-  LatticeGaugeField UdSdU(UGrid);
-
-  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
-  std::cout << GridLogMessage << " Force test for "<<action.action_name()<<std::endl;
-  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
-
-  RealD eps=0.005;
-
-  Gimpl::generate_momenta(P,sRNG,RNG4);
-  Filter.applyFilter(P);
-
-  action.refresh(U,sRNG,RNG4);
-
-  RealD S1 = action.S(U);
-
-  Gimpl::update_field(P,U,eps);
-
-  action.deriv(U,UdSdU);
-  UdSdU = Ta(UdSdU);
-  Filter.applyFilter(UdSdU);
-
-  DumpSliceNorm("Force",UdSdU,Nd-1);
-
-  Gimpl::update_field(P,U,eps);
-
-  RealD S2 = action.S(U);
-
-  LatticeComplex dS(UGrid); dS = Zero();
-  for(int mu=0;mu<Nd;mu++){
-    auto UdSdUmu = PeekIndex<LorentzIndex>(UdSdU,mu);
-    Pmu= PeekIndex<LorentzIndex>(P,mu);
-    dS = dS - trace(Pmu*UdSdUmu)*eps*2.0*2.0;
-  }
-  ComplexD dSpred    = sum(dS);
-  RealD diff =  S2-S1-dSpred.real();
-
-  std::cout<< GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
-  std::cout<< GridLogMessage << "S1 : "<< S1    <<std::endl;
-  std::cout<< GridLogMessage << "S2 : "<< S2    <<std::endl;
-  std::cout<< GridLogMessage << "dS : "<< S2-S1 <<std::endl;
-  std::cout<< GridLogMessage << "dSpred : "<< dSpred.real() <<std::endl;
-  std::cout<< GridLogMessage << "diff : "<< diff<<std::endl;
-  std::cout<< GridLogMessage << "diff/dS : "<< diff/(S2-S1)<<std::endl;
-  std::cout<< GridLogMessage << "*********************************************************"<<std::endl;
-  std::cout<< GridLogMessage << "Done" <<std::endl;
-  std::cout << GridLogMessage << "*********************************************************"<<std::endl;
-}
+#include "ForceTest.h"
 
 int main (int argc, char ** argv)
 {