Note pressure, velocity error norms in refined poiseuille2d

1 files changed, 53 insertions, 14 deletions

diff --git a/apps/adrian/poiseuille2d/poiseuille2d.cpp b/apps/adrian/poiseuille2d/poiseuille2d.cpp
index b2fb48b..034fb08 100644
--- a/apps/adrian/poiseuille2d/poiseuille2d.cpp
+++ b/apps/adrian/poiseuille2d/poiseuille2d.cpp
@@ -91,7 +91,7 @@ void prepareGeometry(Grid2D<T,DESCRIPTOR>& grid)
 }
 
 void prepareLattice(Grid2D<T,DESCRIPTOR>& grid,
-                    Dynamics<T, DESCRIPTOR>& bulkDynamics,
+                    Dynamics<T,DESCRIPTOR>& bulkDynamics,
                     sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryCondition)
 {
   OstreamManager clout(std::cout,"prepareLattice");
@@ -103,7 +103,7 @@ void prepareLattice(Grid2D<T,DESCRIPTOR>& grid,
 
   const T omega = converter.getLatticeRelaxationFrequency();
 
-  sLattice.defineDynamics(sGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>());
+  sLattice.defineDynamics(sGeometry, 0, &instances::getNoDynamics<T,DESCRIPTOR>());
   sLattice.defineDynamics(sGeometry, 1, &bulkDynamics); // bulk
   sLattice.defineDynamics(sGeometry, 2, &bulkDynamics); // walls
   sLattice.defineDynamics(sGeometry, 3, &bulkDynamics); // inflow
@@ -144,9 +144,9 @@ void getResults(const std::string& prefix,
   auto& sGeometry = grid.getSuperGeometry();
 
   SuperVTMwriter2D<T> vtmWriter(prefix + "poiseuille2d");
-  SuperLatticePhysVelocity2D<T, DESCRIPTOR> velocity(sLattice, converter);
-  SuperLatticePhysPressure2D<T, DESCRIPTOR> pressure(sLattice, converter);
-  SuperLatticeGeometry2D<T, DESCRIPTOR> geometry(sLattice, sGeometry);
+  SuperLatticePhysVelocity2D<T,DESCRIPTOR> velocity(sLattice, converter);
+  SuperLatticePhysPressure2D<T,DESCRIPTOR> pressure(sLattice, converter);
+  SuperLatticeGeometry2D<T,DESCRIPTOR> geometry(sLattice, sGeometry);
   vtmWriter.addFunctor(geometry);
   vtmWriter.addFunctor(velocity);
   vtmWriter.addFunctor(pressure);
@@ -169,6 +169,42 @@ void getResults(const std::string& prefix,
   }
 }
 
+void getError(const std::string& prefix,
+              Grid2D<T,DESCRIPTOR>& grid)
+{
+  OstreamManager clout(std::cout, prefix);
+
+  auto& converter = grid.getConverter();
+  auto& sLattice  = grid.getSuperLattice();
+  auto& sGeometry = grid.getSuperGeometry();
+
+  const T maxVelocity = converter.getCharPhysVelocity();
+  const T radius = ly/2;
+  std::vector<T> axisPoint{lx/2, ly/2};
+  std::vector<T> axisDirection{1, 0};
+  Poiseuille2D<T> uSol(axisPoint, axisDirection, maxVelocity, radius);
+
+  const T Lx = converter.getLatticeLength(lx);
+  const T Ly = converter.getLatticeLength(ly);
+  T p0 = 8.*converter.getLatticeViscosity()*converter.getCharLatticeVelocity()*Lx/T(Ly*Ly);
+  AnalyticalLinear2D<T,T> pSol(-converter.getPhysPressure(p0)/lx, 0, converter.getPhysPressure(p0));
+
+  SuperLatticePhysVelocity2D<T,DESCRIPTOR> u(sLattice, converter);
+  SuperLatticePhysPressure2D<T,DESCRIPTOR> p(sLattice, converter);
+  auto fluid = sGeometry.getMaterialIndicator(1);
+
+  int tmp[]= { };
+  T result[2]= { };
+
+  SuperRelativeErrorL2Norm2D<T> velocityError(u, uSol, fluid);
+  velocityError(result, tmp);
+  clout << "velocity-L2-error(rel)=" << result[0] << std::endl;
+
+  SuperRelativeErrorL2Norm2D<T> pressureError(p, pSol, fluid);
+  pressureError(result, tmp);
+  clout << "pressure-L2-error(rel)=" << result[0] << std::endl;
+}
+
 int main(int argc, char* argv[])
 {
   olbInit(&argc, &argv);
@@ -196,26 +232,26 @@ int main(int argc, char* argv[])
     const Vector<T,2> extendY = {0,extend[1]};
 
     coarseGrid.addFineCoupling(fineGrid, origin, extendY);
-    coarseGrid.addCoarseCoupling(fineGrid, origin + Vector<T,2>{coarseDeltaX,0}, extendY);
+    coarseGrid.addCoarseCoupling(fineGrid, origin + Vector<T,2> {coarseDeltaX,0}, extendY);
   }
 
   prepareGeometry(fineGrid);
 
-  BGKdynamics<T, DESCRIPTOR> coarseBulkDynamics(
+  BGKdynamics<T,DESCRIPTOR> coarseBulkDynamics(
     coarseGrid.getConverter().getLatticeRelaxationFrequency(),
-    instances::getBulkMomenta<T, DESCRIPTOR>());
+    instances::getBulkMomenta<T,DESCRIPTOR>());
 
-  sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> coarseSBoundaryCondition(coarseGrid.getSuperLattice());
-  createLocalBoundaryCondition2D<T, DESCRIPTOR>(coarseSBoundaryCondition);
+  sOnLatticeBoundaryCondition2D<T,DESCRIPTOR> coarseSBoundaryCondition(coarseGrid.getSuperLattice());
+  createInterpBoundaryCondition2D<T,DESCRIPTOR>(coarseSBoundaryCondition);
 
   prepareLattice(coarseGrid, coarseBulkDynamics, coarseSBoundaryCondition);
 
-  BGKdynamics<T, DESCRIPTOR> fineBulkDynamics(
+  BGKdynamics<T,DESCRIPTOR> fineBulkDynamics(
     fineGrid.getConverter().getLatticeRelaxationFrequency(),
-    instances::getBulkMomenta<T, DESCRIPTOR>());
+    instances::getBulkMomenta<T,DESCRIPTOR>());
 
-  sOnLatticeBoundaryCondition2D<T, DESCRIPTOR> fineSBoundaryCondition(fineGrid.getSuperLattice());
-  createLocalBoundaryCondition2D<T, DESCRIPTOR>(fineSBoundaryCondition);
+  sOnLatticeBoundaryCondition2D<T,DESCRIPTOR> fineSBoundaryCondition(fineGrid.getSuperLattice());
+  createInterpBoundaryCondition2D<T,DESCRIPTOR>(fineSBoundaryCondition);
 
   prepareLattice(fineGrid, fineBulkDynamics, fineSBoundaryCondition);
 
@@ -252,6 +288,9 @@ int main(int argc, char* argv[])
     converge.takeValue(fineGrid.getSuperLattice().getStatistics().getAverageEnergy(), true);
   }
 
+  getError("coarse", coarseGrid);
+  getError("fine", fineGrid);
+
   timer.stop();
   timer.printSummary();
 }