summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--apps/adrian/poiseuille2d/poiseuille2d.cpp67
-rw-r--r--src/refinement/grid2D.hh6
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(),