Change F2C restriction, some cleanup

1 files changed, 150 insertions, 137 deletions

diff --git a/apps/adrian/poiseuille2d/poiseuille2d.cpp b/apps/adrian/poiseuille2d/poiseuille2d.cpp
index 8bb982d..22c35cb 100644
--- a/apps/adrian/poiseuille2d/poiseuille2d.cpp
+++ b/apps/adrian/poiseuille2d/poiseuille2d.cpp
@@ -27,25 +27,19 @@
 #include "olb2D.hh"
 
 #include <vector>
-#include <cmath>
-#include <iostream>
-#include <iomanip>
-#include <fstream>
 
 using namespace olb;
 using namespace olb::descriptors;
-using namespace olb::graphics;
-using namespace std;
 
 typedef double T;
 
 #define DESCRIPTOR D2Q9Descriptor
 
-const T lx  = 2.;             // length of the channel
+const T lx  = 4.;             // length of the channel
 const T ly  = 1.;             // height of the channel
-const int N = 16;             // resolution of the model
-const T Re  = 10.;            // Reynolds number
-const T maxPhysT = 40.;       // max. simulation time in s, SI unit
+const int N = 30;             // resolution of the model
+const T Re  = 100.;           // Reynolds number
+const T maxPhysT = 60.;       // max. simulation time in s, SI unit
 const T physInterval = 0.25;  // interval for the convergence check in s
 const T residuum = 1e-5;      // residuum for the convergence check
 
@@ -61,14 +55,11 @@ void prepareGeometry(UnitConverter<T,DESCRIPTOR> const& converter,
   // Set material number for bounce back boundaries
   superGeometry.rename(2,1,0,1);
 
-  Vector<T,2> extend;
-  Vector<T,2> origin;
   T physSpacing = converter.getPhysDeltaX();
+  Vector<T,2> extend{  physSpacing / 2, ly };
+  Vector<T,2> origin{ -physSpacing / 4, 0  };
 
   // Set material number for inflow
-  extend[1] = ly;
-  extend[0] = physSpacing / 2;
-  origin[0] -= physSpacing / 4;
   IndicatorCuboid2D<T> inflow(extend, origin);
   superGeometry.rename(1,3,1,inflow);
 
@@ -77,6 +68,9 @@ void prepareGeometry(UnitConverter<T,DESCRIPTOR> const& converter,
   IndicatorCuboid2D<T> outflow(extend, origin);
   superGeometry.rename(1,4,1,outflow);
 
+  //IndicatorCuboid2D<T> obstacle(Vector<T,2>{0.5,0.2}, Vector<T,2>{1.75,0.4});
+  //superGeometry.rename(1,2,obstacle);
+
   // Removes all not needed boundary voxels outside the surface
   superGeometry.clean();
   // Removes all not needed boundary voxels inside the surface
@@ -106,30 +100,28 @@ void prepareCoarseLattice(UnitConverter<T,DESCRIPTOR> const& converter,
 
   sBoundaryCondition.addVelocityBoundary(superGeometry, 3, omega);
 
-  T Lx = converter.getLatticeLength(lx);
-  T Ly = converter.getLatticeLength(ly);
-
-  T p0 =8.*converter.getLatticeViscosity()*converter.getCharLatticeVelocity()*Lx/(Ly*Ly);
+  const T Lx = converter.getLatticeLength(lx);
+  const T Ly = converter.getLatticeLength(ly);
+  const T p0 = 8.*converter.getLatticeViscosity()*converter.getCharLatticeVelocity()*Lx/(Ly*Ly);
   AnalyticalLinear2D<T,T> rho(-p0/lx*DESCRIPTOR<T>::invCs2, 0, p0*DESCRIPTOR<T>::invCs2+1);
 
-  AnalyticalConst2D<T,T> iniRho(1);
-  AnalyticalConst2D<T,T> iniU(std::vector<T> {0.,0.});
-
-  sLattice.defineRhoU(superGeometry, 1, iniRho, iniU);
-  sLattice.iniEquilibrium(superGeometry, 1, iniRho, iniU);
-  sLattice.defineRhoU(superGeometry, 2, iniRho, iniU);
-  sLattice.iniEquilibrium(superGeometry, 2, iniRho, iniU);
+  const T maxVelocity = converter.getCharLatticeVelocity();
+  const T radius = ly/2;
+  std::vector<T> axisPoint{0, ly/2};
+  std::vector<T> axisDirection{1, 0};
+  Poiseuille2D<T> u(axisPoint, axisDirection, maxVelocity, radius);
 
-  T maxVelocity = converter.getCharLatticeVelocity();
-  T distance2Wall = converter.getConversionFactorLength();
-  Poiseuille2D<T> u(superGeometry, 3, maxVelocity, distance2Wall);
+  sLattice.defineRhoU(superGeometry, 1, rho, u);
+  sLattice.iniEquilibrium(superGeometry, 1, rho, u);
+  sLattice.defineRhoU(superGeometry, 2, rho, u);
+  sLattice.iniEquilibrium(superGeometry, 2, rho, u);
 
   sLattice.defineRhoU(superGeometry, 3, rho, u);
   sLattice.iniEquilibrium(superGeometry, 3, rho, u);
 
   sLattice.initialize();
 
-  clout << "Prepare Lattice ... OK" << std::endl;
+  clout << "Prepare coarse lattice ... OK" << std::endl;
 }
 
 void prepareFineLattice(UnitConverter<T,DESCRIPTOR> const& converter,
@@ -139,7 +131,7 @@ void prepareFineLattice(UnitConverter<T,DESCRIPTOR> const& converter,
                         SuperGeometry2D<T>& superGeometry)
 {
   OstreamManager clout(std::cout,"prepareLattice");
-  clout << "Prepare Lattice ..." << std::endl;
+  clout << "Prepare fine lattice ..." << std::endl;
 
   const T omega = converter.getLatticeRelaxationFrequency();
 
@@ -150,25 +142,28 @@ void prepareFineLattice(UnitConverter<T,DESCRIPTOR> const& converter,
 
   sBoundaryCondition.addPressureBoundary(superGeometry, 4, omega);
 
-  T Lx = converter.getLatticeLength(lx);
-  T Ly = converter.getLatticeLength(ly);
-
-  T p0 =8.*converter.getLatticeViscosity()*converter.getCharLatticeVelocity()*Lx/(Ly*Ly);
+  const T Lx = converter.getLatticeLength(lx);
+  const T Ly = converter.getLatticeLength(ly);
+  const T p0 = 8.*converter.getLatticeViscosity()*converter.getCharLatticeVelocity()*Lx/(Ly*Ly);
   AnalyticalLinear2D<T,T> rho(-p0/lx*DESCRIPTOR<T>::invCs2, 0, p0*DESCRIPTOR<T>::invCs2+1);
 
-  AnalyticalConst2D<T,T> iniRho(1);
-  AnalyticalConst2D<T,T> iniU(std::vector<T> {0.,0.});
+  const T maxVelocity = converter.getCharLatticeVelocity();
+  const T radius = ly/2;
+  std::vector<T> axisPoint{lx, ly/2};
+  std::vector<T> axisDirection{1, 0};
+  Poiseuille2D<T> u(axisPoint, axisDirection, maxVelocity, radius);
+
+  sLattice.defineRhoU(superGeometry, 1, rho, u);
+  sLattice.iniEquilibrium(superGeometry, 1, rho, u);
+  sLattice.defineRhoU(superGeometry, 2, rho, u);
+  sLattice.iniEquilibrium(superGeometry, 2, rho, u);
 
-  sLattice.defineRhoU(superGeometry, 1, iniRho, iniU);
-  sLattice.iniEquilibrium(superGeometry, 1, iniRho, iniU);
-  sLattice.defineRhoU(superGeometry, 2, iniRho, iniU);
-  sLattice.iniEquilibrium(superGeometry, 2, iniRho, iniU);
-  sLattice.defineRhoU(superGeometry, 4, iniRho, iniU);
-  sLattice.iniEquilibrium(superGeometry, 4, iniRho, iniU);
+  sLattice.defineRhoU(superGeometry, 4, rho, u);
+  sLattice.iniEquilibrium(superGeometry, 4, rho, u);
 
   sLattice.initialize();
 
-  clout << "Prepare Lattice ... OK" << std::endl;
+  clout << "Prepare fine lattice ... OK" << std::endl;
 }
 
 void getResults(const std::string& prefix,
@@ -214,6 +209,13 @@ private:
   std::unique_ptr<SuperLattice2D<T,DESCRIPTOR>> _lattice;
 
 public:
+  static std::unique_ptr<Grid<T,DESCRIPTOR>> make(IndicatorF2D<T>& domainF, int resolution, T tau, int re)
+  {
+    return std::unique_ptr<Grid<T,DESCRIPTOR>>(
+             new Grid<T,DESCRIPTOR>(domainF, resolution, tau, re)
+           );
+  }
+
   Grid(IndicatorF2D<T>& domainF, int resolution, T tau, int re):
     _converter(new UnitConverterFromResolutionAndRelaxationTime<T,DESCRIPTOR>(
                  resolution, // resolution: number of voxels per charPhysL
@@ -221,7 +223,8 @@ public:
                  T{1},       // charPhysLength: reference length of simulation geometry
                  T{1},       // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
                  T{1./re},   // physViscosity: physical kinematic viscosity in __m^2 / s__
-                 T{1})),     // physDensity: physical density in __kg / m^3__
+                 T{1},       // physDensity: physical density in __kg / m^3__
+                 T{1})),
     _cuboids(new CuboidGeometry2D<T>(
                domainF,
                _converter->getConversionFactorLength(),
@@ -272,6 +275,17 @@ public:
   {
     return 1./getScalingFactor();
   }
+
+  std::unique_ptr<Grid<T,DESCRIPTOR>> refine(IndicatorF2D<T>& domainF)
+  {
+    return std::unique_ptr<Grid<T,DESCRIPTOR>>(
+             new Grid<T,DESCRIPTOR>(
+               domainF,
+               2*getConverter().getResolution(),
+               2.0*getConverter().getLatticeRelaxationTime() - 0.5,
+               getConverter().getReynoldsNumber()
+             ));
+  }
 };
 
 template <typename T, template<typename> class DESCRIPTOR>
@@ -281,8 +295,8 @@ void computeFeq(const Cell<T,DESCRIPTOR>& cell, T fEq[DESCRIPTOR<T>::q])
   T u[2] {};
   cell.computeRhoU(rho, u);
   const T uSqr = u[0]*u[0] + u[1]*u[1];
-  for (int i=0; i < DESCRIPTOR<T>::q; ++i) {
-    fEq[i] = lbHelpers<T,DESCRIPTOR>::equilibrium(i, rho, u, uSqr);
+  for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
+    fEq[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr);
   }
 }
 
@@ -331,8 +345,8 @@ void storeC2F(
 
     T fNeq[DESCRIPTOR<T>::q] {};
     computeFneq(coarseLattice.get(x,y), fNeq);
-    for (int i=0; i < DESCRIPTOR<T>::q; ++i) {
-      c2f_fneq[y][i] = fNeq[i];
+    for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
+      c2f_fneq[y][iPop] = fNeq[iPop];
     }
   }
 }
@@ -358,8 +372,8 @@ void interpolateStoredC2F(
 
     T fNeq[DESCRIPTOR<T>::q] {};
     computeFneq(coarseLattice.get(x,y), fNeq);
-    for (int i=0; i < DESCRIPTOR<T>::q; ++i) {
-      c2f_fneq[y][i] = order2interpolation(fNeq[i], c2f_fneq[y][i]);
+    for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
+      c2f_fneq[y][iPop] = order2interpolation(fNeq[iPop], c2f_fneq[y][iPop]);
     }
   }
 }
@@ -382,12 +396,12 @@ void coupleC2F(
     T fEq[DESCRIPTOR<T>::q] {};
     computeFeq(coarseLattice.get(x,y), fEq);
 
-    for (int i=0; i < DESCRIPTOR<T>::q; ++i) {
-      fineLattice.get(0,2*y)[i] = fEq[i] + coarse.getScalingFactor() * c2f_fneq[y][i];
+    for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
+      fineLattice.get(0,2*y)[iPop] = fEq[iPop] + coarse.getScalingFactor() * c2f_fneq[y][iPop];
     }
   }
 
-  for (int y=2; y < coarseLattice.getNy()-2; ++y) {
+  for (int y=1; y < coarseLattice.getNy()-2; ++y) {
     const T rho = order4interpolation(
                     c2f_rho[y-1],
                     c2f_rho[y],
@@ -409,15 +423,15 @@ void coupleC2F(
            );
     const T uSqr = u[0]*u[0] + u[1]*u[1];
 
-    for (int i=0; i < DESCRIPTOR<T>::q; ++i) {
-      const T fneqi = order4interpolation(
-                        c2f_fneq[y-1][i],
-                        c2f_fneq[y][i],
-                        c2f_fneq[y+1][i],
-                        c2f_fneq[y+2][i]
-                      );
+    for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
+      const T fneq = order4interpolation(
+                       c2f_fneq[y-1][iPop],
+                       c2f_fneq[y][iPop],
+                       c2f_fneq[y+1][iPop],
+                       c2f_fneq[y+2][iPop]
+                     );
 
-      fineLattice.get(0,1+2*y)[i] = lbHelpers<T,DESCRIPTOR>::equilibrium(i, rho, u, uSqr) + coarse.getScalingFactor() * fneqi;
+      fineLattice.get(0,1+2*y)[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + coarse.getScalingFactor() * fneq;
     }
   }
 
@@ -440,56 +454,62 @@ void coupleC2F(
            );
     const T uSqr = u[0]*u[0] + u[1]*u[1];
 
-    for (int i=0; i < DESCRIPTOR<T>::q; ++i) {
-      const T fneqi = order3interpolation(
-                        c2f_fneq[0][i],
-                        c2f_fneq[1][i],
-                        c2f_fneq[2][i]
-                      );
-      fineLattice.get(0,1)[i] = lbHelpers<T,DESCRIPTOR>::equilibrium(i, rho, u, uSqr) + coarse.getScalingFactor() * fneqi;
+    for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
+      const T fneq = order3interpolation(
+                       c2f_fneq[0][iPop],
+                       c2f_fneq[1][iPop],
+                       c2f_fneq[2][iPop]
+                     );
+      fineLattice.get(0,1)[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + coarse.getScalingFactor() * fneq;
     }
   }
 
   {
+    const int maxY = coarseLattice.getNy()-1;
     const T rho = order3interpolation(
-                    c2f_rho[coarseLattice.getNy()-1],
-                    c2f_rho[coarseLattice.getNy()-2],
-                    c2f_rho[coarseLattice.getNy()-3]
+                    c2f_rho[maxY  ],
+                    c2f_rho[maxY-1],
+                    c2f_rho[maxY-2]
                   );
     T u[2] {};
     u[0] = order3interpolation(
-             c2f_u[coarseLattice.getNy()-1][0],
-             c2f_u[coarseLattice.getNy()-2][0],
-             c2f_u[coarseLattice.getNy()-3][0]
+             c2f_u[maxY  ][0],
+             c2f_u[maxY-1][0],
+             c2f_u[maxY-2][0]
            );
     u[1] = order3interpolation(
-             c2f_u[coarseLattice.getNy()-1][1],
-             c2f_u[coarseLattice.getNy()-2][1],
-             c2f_u[coarseLattice.getNy()-3][1]
+             c2f_u[maxY  ][1],
+             c2f_u[maxY-1][1],
+             c2f_u[maxY-2][1]
            );
     const T uSqr = u[0]*u[0] + u[1]*u[1];
 
-    for (int i=0; i < DESCRIPTOR<T>::q; ++i) {
-      const T fneqi = order3interpolation(
-                        c2f_fneq[coarseLattice.getNy()-1][i],
-                        c2f_fneq[coarseLattice.getNy()-2][i],
-                        c2f_fneq[coarseLattice.getNy()-3][i]
-                      );
-      fineLattice.get(0,fineLattice.getNy()-2)[i] = lbHelpers<T,DESCRIPTOR>::equilibrium(i, rho, u, uSqr) + coarse.getScalingFactor() * fneqi;
+    for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
+      const T fneq = order3interpolation(
+                       c2f_fneq[maxY  ][iPop],
+                       c2f_fneq[maxY-1][iPop],
+                       c2f_fneq[maxY-2][iPop]
+                     );
+      fineLattice.get(0,fineLattice.getNy()-2)[iPop] = lbHelpers<T,DESCRIPTOR>::equilibrium(iPop, rho, u, uSqr) + coarse.getScalingFactor() * fneq;
     }
   }
 }
 
 template <typename T, template<typename> class DESCRIPTOR>
-T computeRestrictedFneq(const BlockLatticeView2D<T,DESCRIPTOR>& lattice, int x, int y)
+void computeRestrictedFneq(const BlockLatticeView2D<T,DESCRIPTOR>& lattice, int x, int y, T restrictedFneq[DESCRIPTOR<T>::q])
 {
-  T sum = T{0};
-  for (int i=0; i < DESCRIPTOR<T>::q; ++i) {
+  for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
     T fNeq[DESCRIPTOR<T>::q] {};
-    computeFneq(lattice.get(x + DESCRIPTOR<T>::c[i][0], y + DESCRIPTOR<T>::c[i][1]), fNeq);
-    sum += fNeq[i];
+    computeFneq(lattice.get(x + DESCRIPTOR<T>::c[iPop][0], y + DESCRIPTOR<T>::c[iPop][1]), fNeq);
+
+    for (int jPop=0; jPop < DESCRIPTOR<T>::q; ++jPop) {
+      restrictedFneq[jPop] += fNeq[jPop];
+    }
+  }
+
+  for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
+    restrictedFneq[iPop] /= DESCRIPTOR<T>::q;
   }
-  return sum / DESCRIPTOR<T>::q;
 }
 
 template <typename T, template<typename> class DESCRIPTOR>
@@ -503,10 +523,11 @@ void coupleF2C(Grid<T,DESCRIPTOR>& coarse, Grid<T,DESCRIPTOR>& fine)
     T fEq[DESCRIPTOR<T>::q] {};
     computeFeq(fineLattice.get(2,2*y), fEq);
 
-    const T restrictedFneq = computeRestrictedFneq(fineLattice, 2, 2*y);
+    T fNeq[DESCRIPTOR<T>::q] {};
+    computeRestrictedFneq(fineLattice, 2, 2*y, fNeq);
 
-    for (int i=0; i < DESCRIPTOR<T>::q; ++i) {
-      coarseLattice.get(x,y)[i] = fEq[i] + coarse.getInvScalingFactor() * restrictedFneq;
+    for (int iPop=0; iPop < DESCRIPTOR<T>::q; ++iPop) {
+      coarseLattice.get(x,y)[iPop] = fEq[iPop] + coarse.getInvScalingFactor() * fNeq[iPop];
     }
   }
 }
@@ -526,101 +547,93 @@ int main(int argc, char* argv[])
   Vector<T,2> fineExtend   { 0.5 * lx + coarseDeltaX, ly  };
   IndicatorCuboid2D<T> fineCuboid(fineExtend, fineOrigin);
 
-  Grid<T,DESCRIPTOR> coarseGrid(
-    coarseCuboid,
-    N,
-    0.8,
-    Re);
-  Grid<T,DESCRIPTOR> fineGrid(
-    fineCuboid,
-    2*coarseGrid.getConverter().getResolution(),
-    2.0*coarseGrid.getConverter().getLatticeRelaxationTime() - 0.5,
-    Re);
+  auto coarseGrid = Grid<T,DESCRIPTOR>::make(coarseCuboid, N, 0.8, Re);
+  auto fineGrid   = coarseGrid->refine(fineCuboid);
 
-  prepareGeometry(coarseGrid.getConverter(), coarseGrid.getSuperGeometry());
-  prepareGeometry(fineGrid.getConverter(), fineGrid.getSuperGeometry());
+  prepareGeometry(coarseGrid->getConverter(), coarseGrid->getSuperGeometry());
+  prepareGeometry(fineGrid->getConverter(), fineGrid->getSuperGeometry());
 
   Dynamics<T, DESCRIPTOR>* coarseBulkDynamics;
   coarseBulkDynamics = new BGKdynamics<T, DESCRIPTOR>(
-    coarseGrid.getConverter().getLatticeRelaxationFrequency(),
+    coarseGrid->getConverter().getLatticeRelaxationFrequency(),
     instances::getBulkMomenta<T, DESCRIPTOR>());
 
-  sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> coarseSBoundaryCondition(coarseGrid.getSuperLattice());
+  sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> coarseSBoundaryCondition(coarseGrid->getSuperLattice());
   createLocalBoundaryCondition2D<T, DESCRIPTOR>(coarseSBoundaryCondition);
 
   prepareCoarseLattice(
-    coarseGrid.getConverter(),
-    coarseGrid.getSuperLattice(),
+    coarseGrid->getConverter(),
+    coarseGrid->getSuperLattice(),
     *coarseBulkDynamics,
     coarseSBoundaryCondition,
-    coarseGrid.getSuperGeometry());
+    coarseGrid->getSuperGeometry());
 
   Dynamics<T, DESCRIPTOR>* fineBulkDynamics;
   fineBulkDynamics = new BGKdynamics<T, DESCRIPTOR>(
-    fineGrid.getConverter().getLatticeRelaxationFrequency(),
+    fineGrid->getConverter().getLatticeRelaxationFrequency(),
     instances::getBulkMomenta<T, DESCRIPTOR>());
 
-  sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> fineSBoundaryCondition(fineGrid.getSuperLattice());
+  sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> fineSBoundaryCondition(fineGrid->getSuperLattice());
   createLocalBoundaryCondition2D<T, DESCRIPTOR>(fineSBoundaryCondition);
 
   prepareFineLattice(
-    fineGrid.getConverter(),
-    fineGrid.getSuperLattice(),
+    fineGrid->getConverter(),
+    fineGrid->getSuperLattice(),
     *fineBulkDynamics,
     fineSBoundaryCondition,
-    fineGrid.getSuperGeometry());
+    fineGrid->getSuperGeometry());
 
   clout << "starting simulation..." << endl;
   Timer<T> timer(
-    coarseGrid.getConverter().getLatticeTime(maxPhysT),
-    coarseGrid.getSuperGeometry().getStatistics().getNvoxel());
+    coarseGrid->getConverter().getLatticeTime(maxPhysT),
+    coarseGrid->getSuperGeometry().getStatistics().getNvoxel());
   util::ValueTracer<T> converge(
-    fineGrid.getConverter().getLatticeTime(physInterval),
+    fineGrid->getConverter().getLatticeTime(physInterval),
     residuum);
   timer.start();
 
-  const auto sizeC2F = coarseGrid.getSuperLattice().getBlockLattice(0).getNy();
+  const auto sizeC2F = coarseGrid->getSuperLattice().getBlockLattice(0).getNy();
   T c2f_rho[sizeC2F] {};
   T c2f_u[sizeC2F][2] {};
   T c2f_fneq[sizeC2F][DESCRIPTOR<T>::q] {};
 
-  for (int iT = 0; iT < coarseGrid.getConverter().getLatticeTime(maxPhysT); ++iT) {
+  for (int iT = 0; iT < coarseGrid->getConverter().getLatticeTime(maxPhysT); ++iT) {
     if (converge.hasConverged()) {
       clout << "Simulation converged." << endl;
       break;
     }
 
-    storeC2F(coarseGrid, c2f_rho, c2f_u, c2f_fneq);
-    coarseGrid.getSuperLattice().collideAndStream();
-    interpolateStoredC2F(coarseGrid, c2f_rho, c2f_u, c2f_fneq);
+    storeC2F(*coarseGrid, c2f_rho, c2f_u, c2f_fneq);
+    coarseGrid->getSuperLattice().collideAndStream();
 
-    fineGrid.getSuperLattice().collideAndStream();
-    coupleC2F(coarseGrid, fineGrid, c2f_rho, c2f_u, c2f_fneq);
+    interpolateStoredC2F(*coarseGrid, c2f_rho, c2f_u, c2f_fneq);
+    fineGrid->getSuperLattice().collideAndStream();
+    coupleC2F(*coarseGrid, *fineGrid, c2f_rho, c2f_u, c2f_fneq);
 
-    fineGrid.getSuperLattice().collideAndStream();
-    storeC2F(coarseGrid, c2f_rho, c2f_u, c2f_fneq);
-    coupleC2F(coarseGrid, fineGrid, c2f_rho, c2f_u, c2f_fneq);
+    storeC2F(*coarseGrid, c2f_rho, c2f_u, c2f_fneq);
+    fineGrid->getSuperLattice().collideAndStream();
+    coupleC2F(*coarseGrid, *fineGrid, c2f_rho, c2f_u, c2f_fneq);
 
-    coupleF2C(coarseGrid, fineGrid);
+    coupleF2C(*coarseGrid, *fineGrid);
 
     getResults(
       "coarse_",
-      coarseGrid.getSuperLattice(),
-      coarseGrid.getConverter(),
+      coarseGrid->getSuperLattice(),
+      coarseGrid->getConverter(),
       iT,
-      coarseGrid.getSuperGeometry(),
+      coarseGrid->getSuperGeometry(),
       timer,
       converge.hasConverged());
     getResults(
       "fine_",
-      fineGrid.getSuperLattice(),
-      fineGrid.getConverter(),
+      fineGrid->getSuperLattice(),
+      fineGrid->getConverter(),
       iT,
-      fineGrid.getSuperGeometry(),
+      fineGrid->getSuperGeometry(),
       timer,
       converge.hasConverged());
 
-    converge.takeValue(fineGrid.getSuperLattice().getStatistics().getAverageEnergy(), true);
+    converge.takeValue(fineGrid->getSuperLattice().getStatistics().getAverageEnergy(), true);
   }
 
   timer.stop();