Cache lattice positions for grid coupling

2 files changed, 43 insertions, 45 deletions

diff --git a/src/refinement/coupler2D.h b/src/refinement/coupler2D.h
index b218eff..e911575 100644
--- a/src/refinement/coupler2D.h
+++ b/src/refinement/coupler2D.h
@@ -39,12 +39,14 @@ protected:
   const int  _fineSize;
   const bool _vertical;
 
-  Vector<T,2>   _physOrigin;
-  Vector<int,3> _coarseOrigin;
-  Vector<int,3> _fineOrigin;
+  Vector<T,2> _physOrigin;
 
-  Vector<int,3> getFineLatticeR(int y) const;
-  Vector<int,3> getCoarseLatticeR(int y) const;
+  const Vector<int,3>& getFineLatticeR(int y) const;
+  const Vector<int,3>& getCoarseLatticeR(int y) const;
+
+private:
+  std::vector<Vector<int,3>> _coarseLatticeR;
+  std::vector<Vector<int,3>> _fineLatticeR;
 
 public:
   Coupler2D(Grid2D<T,DESCRIPTOR>& coarse, Grid2D<T,DESCRIPTOR>& fine,
diff --git a/src/refinement/coupler2D.hh b/src/refinement/coupler2D.hh
index ab082f4..5614d16 100644
--- a/src/refinement/coupler2D.hh
+++ b/src/refinement/coupler2D.hh
@@ -32,37 +32,15 @@ namespace olb {
 
 
 template <typename T, template<typename> class DESCRIPTOR>
-Vector<int,3> Coupler2D<T,DESCRIPTOR>::getFineLatticeR(int y) const
+const Vector<int,3>& Coupler2D<T,DESCRIPTOR>::getFineLatticeR(int y) const
 {
-  if (_vertical) {
-    const auto physR = _physOrigin + Vector<T,2>{0, y*_fine.getConverter().getPhysDeltaX()};
-    Vector<int,3> latticeR{};
-    _fine.getCuboidGeometry().getLatticeR(physR, latticeR);
-    return latticeR;
-  }
-  else {
-    const auto physR = _physOrigin + Vector<T,2>{y*_fine.getConverter().getPhysDeltaX(), 0};
-    Vector<int,3> latticeR{};
-    _fine.getCuboidGeometry().getLatticeR(physR, latticeR);
-    return latticeR;
-  }
+  return _fineLatticeR[y];
 }
 
 template <typename T, template<typename> class DESCRIPTOR>
-Vector<int,3> Coupler2D<T,DESCRIPTOR>::getCoarseLatticeR(int y) const
+const Vector<int,3>& Coupler2D<T,DESCRIPTOR>::getCoarseLatticeR(int y) const
 {
-  if (_vertical) {
-    const auto physR = _physOrigin + Vector<T,2>{0, y*_coarse.getConverter().getPhysDeltaX()};
-    Vector<int,3> latticeR{};
-    _coarse.getCuboidGeometry().getLatticeR(physR, latticeR);
-    return latticeR;
-  }
-  else {
-    const auto physR = _physOrigin + Vector<T,2>{y*_coarse.getConverter().getPhysDeltaX(), 0};
-    Vector<int,3> latticeR{};
-    _coarse.getCuboidGeometry().getLatticeR(physR, latticeR);
-    return latticeR;
-  }
+  return _coarseLatticeR[y];
 }
 
 template <typename T, template<typename> class DESCRIPTOR>
@@ -72,7 +50,9 @@ Coupler2D<T,DESCRIPTOR>::Coupler2D(Grid2D<T,DESCRIPTOR>& coarse, Grid2D<T,DESCRI
   _fine(fine),
   _coarseSize(floor(extend.norm() / coarse.getConverter().getPhysDeltaX() + 0.5)+1),
   _fineSize(2*_coarseSize-1),
-  _vertical(util::nearZero(extend[0]))
+  _vertical(util::nearZero(extend[0])),
+  _coarseLatticeR(_coarseSize),
+  _fineLatticeR(_fineSize)
 {
   OLB_ASSERT(util::nearZero(extend[0]) || util::nearZero(extend[1]), "Coupling is only implemented alongside unit vectors");
 
@@ -83,8 +63,16 @@ Coupler2D<T,DESCRIPTOR>::Coupler2D(Grid2D<T,DESCRIPTOR>& coarse, Grid2D<T,DESCRI
   coarseGeometry.getLatticeR(origin, tmpLatticeOrigin);
   _physOrigin = coarseGeometry.getPhysR(tmpLatticeOrigin.toStdVector());
 
-  coarseGeometry.getLatticeR(_physOrigin, _coarseOrigin);
-  fineGeometry.getLatticeR(_physOrigin, _fineOrigin);
+  const T deltaX = _fine.getConverter().getPhysDeltaX();
+  const Vector<T,2> stepPhysR = _vertical ? Vector<T,2>{0, deltaX} : Vector<T,2>{deltaX, 0};
+
+  for (int i=0; i < _fineSize; ++i) {
+    if (i % 2 == 0) {
+      coarseGeometry.getLatticeR(_physOrigin + i*stepPhysR, _coarseLatticeR[i/2]);
+    }
+
+    fineGeometry.getLatticeR(_physOrigin + i*stepPhysR, _fineLatticeR[i]);
+  }
 }
 
 
@@ -97,9 +85,13 @@ FineCoupler2D<T,DESCRIPTOR>::FineCoupler2D(Grid2D<T,DESCRIPTOR>& coarse, Grid2D<
   _c2f_fneq(this->_coarseSize, Vector<T,DESCRIPTOR<T>::q>(T{}))
 {
   OstreamManager clout(std::cout,"C2F");
-  clout << "coarse origin: " << this->_coarseOrigin[0] << " " << this->_coarseOrigin[1] << " " << this->_coarseOrigin[2] << std::endl;
-  clout << "fine origin:   " << this->_fineOrigin[0] << " " << this->_fineOrigin[1] << " " << this->_fineOrigin[2] << std::endl;
-  clout << "fine size:   " << this->_fineSize << std::endl;
+
+  const auto& coarseOrigin = this->getCoarseLatticeR(0);
+  const auto& fineOrigin   = this->getFineLatticeR(0);
+
+  clout << "coarse origin: " << coarseOrigin[0] << " " << coarseOrigin[1] << " " << coarseOrigin[2] << std::endl;
+  clout << "fine origin:   " << fineOrigin[0]   << " " << fineOrigin[1]   << " " << fineOrigin[2]   << std::endl;
+  clout << "fine size:     " << this->_fineSize << std::endl;
 }
 
 template <typename T, template<typename> class DESCRIPTOR>
@@ -172,8 +164,8 @@ void FineCoupler2D<T,DESCRIPTOR>::couple()
   auto& fineLattice   = this->_fine.getSuperLattice();
 
   for (int y=0; y < this->_coarseSize; ++y) {
-    const auto coarsePos = this->getCoarseLatticeR(y);
-    const auto finePos   = this->getFineLatticeR(2*y);
+    const auto& coarsePos = this->getCoarseLatticeR(y);
+    const auto& finePos   = this->getFineLatticeR(2*y);
 
     T fEq[DESCRIPTOR<T>::q] {};
     Cell<T,DESCRIPTOR> coarseCell;
@@ -247,7 +239,7 @@ void FineCoupler2D<T,DESCRIPTOR>::couple()
            );
     const T uSqr = u[0]*u[0] + u[1]*u[1];
 
-    const auto finePos = this->getFineLatticeR(1);
+    const auto& finePos = this->getFineLatticeR(1);
     Cell<T,DESCRIPTOR> fineCell;
     fineLattice.get(finePos, fineCell);
 
@@ -282,7 +274,7 @@ void FineCoupler2D<T,DESCRIPTOR>::couple()
            );
     const T uSqr = u[0]*u[0] + u[1]*u[1];
 
-    const auto finePos = this->getFineLatticeR(this->_fineSize-2);
+    const auto& finePos = this->getFineLatticeR(this->_fineSize-2);
     Cell<T,DESCRIPTOR> fineCell;
     fineLattice.get(finePos, fineCell);
 
@@ -329,8 +321,12 @@ CoarseCoupler2D<T,DESCRIPTOR>::CoarseCoupler2D(
   Coupler2D<T,DESCRIPTOR>(coarse, fine, origin, extend)
 {
   OstreamManager clout(std::cout,"F2C");
-  clout << "coarse origin: " << this->_coarseOrigin[0] << " " << this->_coarseOrigin[1] << " " << this->_coarseOrigin[2] << std::endl;
-  clout << "fine origin:   " << this->_fineOrigin[0] << " " << this->_fineOrigin[1] << " " << this->_fineOrigin[2] << std::endl;
+
+  const auto& coarseOrigin = this->getCoarseLatticeR(0);
+  const auto& fineOrigin   = this->getFineLatticeR(0);
+
+  clout << "coarse origin: " << coarseOrigin[0] << " " << coarseOrigin[1] << " " << coarseOrigin[2] << std::endl;
+  clout << "fine origin:   " << fineOrigin[0]   << " " << fineOrigin[1]   << " " << fineOrigin[2]   << std::endl;
   clout << "coarse size:   " << this->_coarseSize << std::endl;
 }
 
@@ -341,8 +337,8 @@ void CoarseCoupler2D<T,DESCRIPTOR>::couple()
   auto& coarseLattice = this->_coarse.getSuperLattice();
 
   for (int y=0; y < this->_coarseSize; ++y) {
-    const auto finePos   = this->getFineLatticeR(2*y);
-    const auto coarsePos = this->getCoarseLatticeR(y);
+    const auto& finePos   = this->getFineLatticeR(2*y);
+    const auto& coarsePos = this->getCoarseLatticeR(y);
 
     T fEq[DESCRIPTOR<T>::q] {};
     Cell<T,DESCRIPTOR> fineCell;