diff options
| -rw-r--r-- | apps/adrian/poiseuille2d/poiseuille2d.cpp | 67 | ||||
| -rw-r--r-- | 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<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(); } 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<T,DESCRIPTOR>::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<T>( *_domainF, _converter->getConversionFactorLength(), @@ -82,8 +81,7 @@ Grid2D<T,DESCRIPTOR>::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<T>( *_domainF, _converter->getConversionFactorLength(), |