Commit c668aaf9 authored by Kilian Weishaupt's avatar Kilian Weishaupt Committed by Timo Koch
Browse files

[staggeredGrid][tests] Improve Kovasznay test

* Add ref solution
* Add L2 error norm function
parent a2e60866
......@@ -38,9 +38,11 @@ add_dumux_test(test_donea test_donea test_donea.cc
--command "${CMAKE_CURRENT_BINARY_DIR}/test_donea")
add_dumux_test(test_kovasznay test_kovasznay test_kovasznay.cc
${CMAKE_CURRENT_BINARY_DIR}/test_kovasznay)
python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
--script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_kovasznay-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_kovasznay-00002.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_kovasznay")
#install sources
install(FILES
......
......@@ -143,6 +143,11 @@ public:
Problem,
Name);
printL2Error_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag,
bool,
Problem,
PrintL2Error);
kinematicViscosity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, LiquidKinematicViscosity);
Scalar reynoldsNumber = 1.0 / kinematicViscosity_;
lambda_ = 0.5 * reynoldsNumber
......@@ -176,6 +181,23 @@ public:
return false;
}
void postTimeStep() const
{
if(printL2Error_)
{
const auto l2error = calculateL2Error();
const int numCellCenterDofs = this->model().numCellCenterDofs();
const int numFaceDofs = this->model().numFaceDofs();
std::cout << std::setprecision(8) << "** L2 error (abs/rel) for "
<< std::setw(6) << numCellCenterDofs << " cc dofs and " << numFaceDofs << " face dofs (total: " << numCellCenterDofs + numFaceDofs << "): "
<< std::scientific
<< "L2(p) = " << l2error.first[pressureIdx] << " / " << l2error.second[pressureIdx]
<< ", L2(vx) = " << l2error.first[velocityXIdx] << " / " << l2error.second[velocityXIdx]
<< ", L2(vy) = " << l2error.first[velocityYIdx] << " / " << l2error.second[velocityYIdx]
<< std::endl;
}
}
/*!
* \brief Return the temperature within the domain in [K].
*
......@@ -329,12 +351,18 @@ public:
* \brief Calculate the L2 error between the analytical solution and the numerical approximation.
*
*/
BoundaryValues calculateL2Error()
auto calculateL2Error() const
{
BoundaryValues sumError(0.0), sumReference(0.0), l2NormAbs(0.0), l2NormRel(0.0);
BoundaryValues analyticalSolution(0.0);
BoundaryValues numericalSolution(0.0); // TODO type
int numElements(0);
const int numFaceDofs = this->model().numFaceDofs();
std::vector<Scalar> staggeredVolume(numFaceDofs);
std::vector<Scalar> errorVelocity(numFaceDofs);
std::vector<Scalar> velocityReference(numFaceDofs);
std::vector<int> directionIndex(numFaceDofs);
Scalar totalVolume = 0.0;
for (const auto& element : elements(this->gridView()))
{
......@@ -343,60 +371,60 @@ public:
for (auto&& scv : scvs(fvGeometry))
{
auto dofIdxGlobal = scv.dofIndex();
const auto& globalPos = scv.dofPosition();
analyticalSolution = dirichletAtPos(globalPos);
numericalSolution[pressureIdx] = this->model().curSol()[cellCenterIdx][pressureIdx][dofIdxGlobal]; // Multitypeblockvector
numericalSolution[velocityXIdx] = this->model().curSol()[faceIdx][velocityXIdx][dofIdxGlobal];
numericalSolution[velocityYIdx] = this->model().curSol()[faceIdx][velocityYIdx][dofIdxGlobal];
sumError[pressureIdx] += (analyticalSolution[pressureIdx] - numericalSolution[pressureIdx]) * (analyticalSolution[pressureIdx] - numericalSolution[pressureIdx]);
sumError[velocityXIdx] += (analyticalSolution[velocityXIdx] - numericalSolution[velocityXIdx]) * (analyticalSolution[velocityXIdx] - numericalSolution[velocityXIdx]);
sumError[velocityYIdx] += (analyticalSolution[velocityYIdx] - numericalSolution[velocityYIdx]) * (analyticalSolution[velocityYIdx] - numericalSolution[velocityYIdx]);
sumReference[pressureIdx] += (analyticalSolution[pressureIdx]) * (analyticalSolution[pressureIdx]);
sumReference[velocityXIdx] += (analyticalSolution[velocityXIdx]) * (analyticalSolution[velocityXIdx]);
sumReference[velocityYIdx] += (analyticalSolution[velocityYIdx]) * (analyticalSolution[velocityYIdx]);
++numElements;
// treat cell-center dofs
const auto dofIdxCellCenter = scv.dofIndex();
const auto& posCellCenter = scv.dofPosition();
const auto analyticalSolutionCellCenter = dirichletAtPos(posCellCenter)[cellCenterIdx];
const auto numericalSolutionCellCenter = this->model().curSol()[cellCenterIdx][dofIdxCellCenter];
sumError[cellCenterIdx] += squaredDiff(analyticalSolutionCellCenter, numericalSolutionCellCenter) * scv.volume();
sumReference[cellCenterIdx] += analyticalSolutionCellCenter * analyticalSolutionCellCenter * scv.volume();
totalVolume += scv.volume();
// treat face dofs
for (auto&& scvf : scvfs(fvGeometry))
{
const int dofIdxFace = scvf.dofIndexSelf();
const int dirIdx = scvf.directionIndex();
const auto analyticalSolutionFace = dirichletAtPos(scvf.center())[faceIdx][dirIdx];
const auto numericalSolutionFace = this->model().curSol()[faceIdx][dofIdxFace][momentumBalanceIdx];
directionIndex[dofIdxFace] = dirIdx;
errorVelocity[dofIdxFace] = squaredDiff(analyticalSolutionFace, numericalSolutionFace);
velocityReference[dofIdxFace] = squaredDiff(analyticalSolutionFace, 0.0);
const Scalar staggeredHalfVolume = 0.5 * scv.volume();
staggeredVolume[dofIdxFace] = staggeredVolume[dofIdxFace] + staggeredHalfVolume;
}
}
}
// l2NormAbs = std::sqrt(sumError / numElements);
l2NormRel[pressureIdx] = std::sqrt(sumError[pressureIdx] / numElements / sumReference[pressureIdx]);
l2NormRel[velocityXIdx] = std::sqrt(sumError[velocityXIdx] / numElements / sumReference[velocityXIdx]);
l2NormRel[velocityYIdx] = std::sqrt(sumError[velocityYIdx] / numElements / sumReference[velocityYIdx]);
// get the absolute and relative discrete L2-error for cell-center dofs
l2NormAbs[cellCenterIdx] = std::sqrt(sumError[cellCenterIdx] / totalVolume);
l2NormRel[cellCenterIdx] = std::sqrt(sumError[cellCenterIdx] / sumReference[cellCenterIdx]);
return l2NormRel; // TODO return both errors??
}
// get the absolute and relative discrete L2-error for face dofs
for(int i = 0; i < numFaceDofs; ++i)
{
const int dirIdx = directionIndex[i];
const auto error = errorVelocity[i];
const auto ref = velocityReference[i];
const auto volume = staggeredVolume[i];
sumError[faceIdx][dirIdx] += error * volume;
sumReference[faceIdx][dirIdx] += ref * volume;
}
/*!
* \brief Write the L2 error into an output file
*
*/
void writeOutput(const bool verbose = true)
{
ParentType::writeOutput(verbose);
BoundaryValues l2error = calculateL2Error();
std::cout.precision(8);
std::cout << "** L2 error for "
<< std::setw(6) << this->gridView().size(0)
" elements: "
std::scientific
"L2(p) = "
l2error[pressureIdx]
", L2(vx) = "
l2error[velocityXIdx]
<< ", L2(vy) = "
l2error[velocityYIdx]
std::endl;
for(int dirIdx = 0; dirIdx < dimWorld; ++dirIdx)
{
l2NormAbs[faceIdx][dirIdx] = std::sqrt(sumError[faceIdx][dirIdx] / totalVolume);
l2NormRel[faceIdx][dirIdx] = std::sqrt(sumError[faceIdx][dirIdx] / sumReference[faceIdx][dirIdx]);
}
return std::make_pair(l2NormAbs, l2NormRel);
}
private:
template<class T>
T squaredDiff(const T& a, const T& b) const
{
return (a-b)*(a-b);
}
bool isLowerLeftCell_(const GlobalPosition& globalPos) const
{
......@@ -409,6 +437,7 @@ private:
Scalar kinematicViscosity_;
Scalar lambda_;
bool printL2Error_;
};
} //end namespace
......
[TimeManager]
DtInitial = 1e-1 # [s]
DtInitial = 1 # [s]
TEnd = 2 # [s]
[Grid]
LowerLeft = -0.5 -0.5
UpperRight = 2 1.5
Cells = 100 100
Cells = 50 50
[Problem]
Name = test_kovasznay # name passed to the output routines
InletVelocity = 0
LiquidDensity = 1
EnableGravity = false
LiquidKinematicViscosity = 0.025
PrintL2Error = false
[ Newton ]
MaxSteps = 10
......
This diff is collapsed.
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment