From a28f7ad5314cb5fd4b4582ebfe262516d04e9a74 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Thu, 3 Sep 2026 23:18:06 -0400
Subject: [PATCH] Grids in FermionOperator

---
 .../fermion/ContinuedFractionFermion5D.h      |  3 +--
 Grid/qcd/action/fermion/DomainWallFermion.h   |  3 +--
 Grid/qcd/action/fermion/FermionOperator.h     | 25 +++++++++++--------
 .../action/fermion/PartialFractionFermion5D.h |  3 +--
 4 files changed, 18 insertions(+), 16 deletions(-)

diff --git a/Grid/qcd/action/fermion/ContinuedFractionFermion5D.h b/Grid/qcd/action/fermion/ContinuedFractionFermion5D.h
index 24886ee28..0e4f3769e 100644
--- a/Grid/qcd/action/fermion/ContinuedFractionFermion5D.h
+++ b/Grid/qcd/action/fermion/ContinuedFractionFermion5D.h
@@ -66,8 +66,7 @@ public:
     FermionField in_k(in.Grid());
     FermionField prop_k(in.Grid());
 
-    GRID_ASSERT(in.Grid() == this->FermionGrid());
-    PlannedFFT<typename FermionField::vector_object> &theFFT = this->FermionGridFFT();
+    PlannedFFT<typename FermionField::vector_object> &theFFT = this->ThePlannedFFT(in.Grid());
 
     //phase for boundary condition
     ComplexField coor(in.Grid());
diff --git a/Grid/qcd/action/fermion/DomainWallFermion.h b/Grid/qcd/action/fermion/DomainWallFermion.h
index 6180824dd..680c203eb 100644
--- a/Grid/qcd/action/fermion/DomainWallFermion.h
+++ b/Grid/qcd/action/fermion/DomainWallFermion.h
@@ -45,8 +45,7 @@ public:
 	FermionField in_k(in.Grid());
 	FermionField prop_k(in.Grid());
 
-	GRID_ASSERT(in.Grid() == this->FermionGrid());
-	PlannedFFT<typename FermionField::vector_object> &theFFT = this->FermionGridFFT();
+	PlannedFFT<typename FermionField::vector_object> &theFFT = this->ThePlannedFFT(in.Grid());
 
 	//phase for boundary condition
 	ComplexField coor(in.Grid());
diff --git a/Grid/qcd/action/fermion/FermionOperator.h b/Grid/qcd/action/fermion/FermionOperator.h
index 867f1183d..a2abaa781 100644
--- a/Grid/qcd/action/fermion/FermionOperator.h
+++ b/Grid/qcd/action/fermion/FermionOperator.h
@@ -45,7 +45,7 @@ public:
   INHERIT_IMPL_TYPES(Impl);
 
   FermionOperator(const ImplParams &p= ImplParams()) : Impl(p) {};
-  virtual ~FermionOperator(void) { if ( _fermionGridFFT ) delete _fermionGridFFT; }
+  virtual ~FermionOperator(void) { if ( _thePlannedFFT ) delete _thePlannedFFT; }
 
   virtual FermionField &tmp(void) = 0;
 
@@ -96,22 +96,27 @@ public:
   virtual void  MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _m,std::vector<double> twist) { GRID_ASSERT(0);};
 
 protected:
-  // Cached planned FFT on the fermion grid -- a general utility (FreePropagator
+  // Cached planned FFT for the ACTIVE field grid -- a general utility (FreePropagator
   // today; smoother, smearings and other users anticipated).  Frontier FFTW plan
   // create+destroy is ~22 ms/call (measured, Test_fft_prop PLANCOST) -- ~4x the
   // transform itself and ~80% of an unplanned call -- so a per-call `FFT theFFT(grid)`
-  // dominates.  Lazily built once on FermionGrid() and reused; owned, deleted in the dtor.
-  PlannedFFT<typename FermionField::vector_object> *_fermionGridFFT{nullptr};
-  PlannedFFT<typename FermionField::vector_object> & FermionGridFFT(void) {
-    if ( _fermionGridFFT == nullptr )
-      _fermionGridFFT = new PlannedFFT<typename FermionField::vector_object>((GridCartesian *)this->FermionGrid());
-    return *_fermionGridFFT;
+  // dominates.  We cache BOTH the FFT and the grid it was built on; when called with a
+  // different grid (e.g. the 5D propagator vs the 4D physical propagator) we rebuild.
+  // Lazily built, owned, deleted in the dtor.  Raw pointer, no smart pointers.
+  PlannedFFT<typename FermionField::vector_object> *_thePlannedFFT{nullptr};
+  GridBase                                         *_thePlannedFFTGrid{nullptr};
+  PlannedFFT<typename FermionField::vector_object> & ThePlannedFFT(GridBase *grid) {
+    if ( _thePlannedFFT == nullptr || _thePlannedFFTGrid != grid ) {
+      if ( _thePlannedFFT != nullptr ) delete _thePlannedFFT;   // active grid changed: rebuild
+      _thePlannedFFT     = new PlannedFFT<typename FermionField::vector_object>((GridCartesian *)grid);
+      _thePlannedFFTGrid = grid;
+    }
+    return *_thePlannedFFT;
   }
 public:
   virtual void  FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector<Complex> boundary,std::vector<double> twist)
       {
-	GRID_ASSERT(in.Grid() == this->FermionGrid());
-	PlannedFFT<typename FermionField::vector_object> &theFFT = this->FermionGridFFT();
+	PlannedFFT<typename FermionField::vector_object> &theFFT = this->ThePlannedFFT(in.Grid());
 
 	typedef typename Simd::scalar_type Scalar;
 
diff --git a/Grid/qcd/action/fermion/PartialFractionFermion5D.h b/Grid/qcd/action/fermion/PartialFractionFermion5D.h
index 68c9680c3..2933d5b12 100644
--- a/Grid/qcd/action/fermion/PartialFractionFermion5D.h
+++ b/Grid/qcd/action/fermion/PartialFractionFermion5D.h
@@ -96,8 +96,7 @@ public:
     FermionField in_k(in.Grid());
     FermionField prop_k(in.Grid());
 
-    GRID_ASSERT(in.Grid() == this->FermionGrid());
-    PlannedFFT<typename FermionField::vector_object> &theFFT = this->FermionGridFFT();
+    PlannedFFT<typename FermionField::vector_object> &theFFT = this->ThePlannedFFT(in.Grid());
 
     //phase for boundary condition
     ComplexField coor(in.Grid());