From 5aeafca5883f7af387363ede78596ff665b36d0c Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Thu, 24 Jan 2019 13:42:57 +0100 Subject: Note pressure, velocity error norms in refined poiseuille2d --- apps/adrian/poiseuille2d/poiseuille2d.cpp | 67 ++++++++++++++++++++++++------- src/refinement/grid2D.hh | 6 +-- 2 files changed, 55 insertions(+), 18 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& grid) } void prepareLattice(Grid2D& grid, - Dynamics& bulkDynamics, + Dynamics& bulkDynamics, sOnLatticeBoundaryCondition2D& sBoundaryCondition) { OstreamManager clout(std::cout,"prepareLattice"); @@ -103,7 +103,7 @@ void prepareLattice(Grid2D& grid, const T omega = converter.getLatticeRelaxationFrequency(); - sLattice.defineDynamics(sGeometry, 0, &instances::getNoDynamics()); + sLattice.defineDynamics(sGeometry, 0, &instances::getNoDynamics()); 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 vtmWriter(prefix + "poiseuille2d"); - SuperLatticePhysVelocity2D velocity(sLattice, converter); - SuperLatticePhysPressure2D pressure(sLattice, converter); - SuperLatticeGeometry2D geometry(sLattice, sGeometry); + SuperLatticePhysVelocity2D velocity(sLattice, converter); + SuperLatticePhysPressure2D pressure(sLattice, converter); + SuperLatticeGeometry2D 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& 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 axisPoint{lx/2, ly/2}; + std::vector axisDirection{1, 0}; + Poiseuille2D 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 pSol(-converter.getPhysPressure(p0)/lx, 0, converter.getPhysPressure(p0)); + + SuperLatticePhysVelocity2D u(sLattice, converter); + SuperLatticePhysPressure2D p(sLattice, converter); + auto fluid = sGeometry.getMaterialIndicator(1); + + int tmp[]= { }; + T result[2]= { }; + + SuperRelativeErrorL2Norm2D velocityError(u, uSol, fluid); + velocityError(result, tmp); + clout << "velocity-L2-error(rel)=" << result[0] << std::endl; + + SuperRelativeErrorL2Norm2D 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 extendY = {0,extend[1]}; coarseGrid.addFineCoupling(fineGrid, origin, extendY); - coarseGrid.addCoarseCoupling(fineGrid, origin + Vector{coarseDeltaX,0}, extendY); + coarseGrid.addCoarseCoupling(fineGrid, origin + Vector {coarseDeltaX,0}, extendY); } prepareGeometry(fineGrid); - BGKdynamics coarseBulkDynamics( + BGKdynamics coarseBulkDynamics( coarseGrid.getConverter().getLatticeRelaxationFrequency(), - instances::getBulkMomenta()); + instances::getBulkMomenta()); - sOnLatticeBoundaryCondition2D coarseSBoundaryCondition(coarseGrid.getSuperLattice()); - createLocalBoundaryCondition2D(coarseSBoundaryCondition); + sOnLatticeBoundaryCondition2D coarseSBoundaryCondition(coarseGrid.getSuperLattice()); + createInterpBoundaryCondition2D(coarseSBoundaryCondition); prepareLattice(coarseGrid, coarseBulkDynamics, coarseSBoundaryCondition); - BGKdynamics fineBulkDynamics( + BGKdynamics fineBulkDynamics( fineGrid.getConverter().getLatticeRelaxationFrequency(), - instances::getBulkMomenta()); + instances::getBulkMomenta()); - sOnLatticeBoundaryCondition2D fineSBoundaryCondition(fineGrid.getSuperLattice()); - createLocalBoundaryCondition2D(fineSBoundaryCondition); + sOnLatticeBoundaryCondition2D fineSBoundaryCondition(fineGrid.getSuperLattice()); + createInterpBoundaryCondition2D(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(); } diff --git a/src/refinement/grid2D.hh b/src/refinement/grid2D.hh index 2bf2f7f..f1caf84 100644 --- a/src/refinement/grid2D.hh +++ b/src/refinement/grid2D.hh @@ -46,8 +46,7 @@ Grid2D::Grid2D( 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})), + T{1})), // physDensity: physical density in __kg / m^3__ _cuboids(new CuboidGeometry2D( *_domainF, _converter->getConversionFactorLength(), @@ -82,8 +81,7 @@ Grid2D::Grid2D( 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})), + T{1})), // physDensity: physical density in __kg / m^3__ _cuboids(new CuboidGeometry2D( *_domainF, _converter->getConversionFactorLength(), -- cgit v1.2.3