Commit e385f8dc by Thomas Fetzer

### [freeflow][navierstokes] Improve calculation of l2error

parent 649fcea5
 ... ... @@ -21,12 +21,15 @@ * \ingroup NavierStokesTests * \copydoc Dumux::NavierStokesTestL2Error */ #ifndef DUMUX_TEST_L2_ERROR_HH #define DUMUX_TEST_L2_ERROR_HH #ifndef DUMUX_NAVIERSTOKES_STAGGERED_L2_ERROR_HH #define DUMUX_NAVIERSTOKES_STAGGERED_L2_ERROR_HH #include #include // #include #include namespace Dumux { ... ... @@ -34,23 +37,54 @@ namespace Dumux * \ingroup NavierStokesTests * \brief Routines to calculate the discrete L2 error */ template template class NavierStokesTestL2Error { using Indices = typename ModelTraits::Indices; public: /*! * \brief Calculate the L2 error between the analytical solution and the numerical approximation. * * \param problem The problem * \param curSol Vector containing the current solution * \param pOrder Polynomial order of quadrature rule */ template static void printL2Error(const Problem& problem, const SolutionVector& curSol, int pOrder = 1) { const auto l2error = calculateL2Error(problem, curSol, pOrder); std::cout << std::setprecision(8) << "** L2 error (abs/rel) for " << std::setw(6) << problem.fvGridGeometry().numCellCenterDofs() << " cc dofs and " << problem.fvGridGeometry().numFaceDofs() << " face dofs (total: " << problem.fvGridGeometry().numCellCenterDofs() + problem.fvGridGeometry().numFaceDofs() << "): " << std::scientific << "L2(p) = " << l2error.first[Indices::pressureIdx] << " / " << l2error.second[Indices::pressureIdx]; if (verbose) { std::cout << " , L2(vx) = " << l2error.first[Indices::velocityXIdx] << " / " << l2error.second[Indices::velocityXIdx] << " , L2(vy) = " << l2error.first[Indices::velocityYIdx] << " / " << l2error.second[Indices::velocityYIdx]; } else { std::cout << " , L2(v) = " << l2error.first[Indices::velocity(0)] << " / " << l2error.second[Indices::velocity(0)]; } std::cout << std::endl; } /*! * \brief Calculate the L2 error between the analytical solution and the numerical approximation. * * \param problem The problem * \param curSol Vector containing the current solution * \param pOrder Polynomial order of quadrature rule */ template static auto calculateL2Error(const Problem& problem, const SolutionVector& curSol) static auto calculateL2Error(const Problem& problem, const SolutionVector& curSol, int pOrder = 1) { using FVGridGeometry = std::decay_t; PrimaryVariables sumError(0.0), sumReference(0.0), l2NormAbs(0.0), l2NormRel(0.0); PrimaryVariables sumError(0.0), sumReference(0.0), totalVolume(0.0), l2NormAbs(0.0), l2NormRel(0.0); const int numFaceDofs = problem.fvGridGeometry().numFaceDofs(); ... ... @@ -59,8 +93,6 @@ public: std::vector velocityReference(numFaceDofs); std::vector directionIndex(numFaceDofs); Scalar totalVolume = 0.0; for (const auto& element : elements(problem.fvGridGeometry().gridView())) { auto fvGeometry = localView(problem.fvGridGeometry()); ... ... @@ -68,49 +100,66 @@ public: for (auto&& scv : scvs(fvGeometry)) { // treat cell-center dofs const auto dofIdxCellCenter = scv.dofIndex(); const auto& posCellCenter = scv.dofPosition(); const auto analyticalSolutionCellCenter = problem.dirichletAtPos(posCellCenter)[Indices::pressureIdx]; const auto numericalSolutionCellCenter = curSol[FVGridGeometry::cellCenterIdx()][dofIdxCellCenter][Indices::pressureIdx - ModelTraits::dim()]; sumError[Indices::pressureIdx] += squaredDiff_(analyticalSolutionCellCenter, numericalSolutionCellCenter) * scv.volume(); sumReference[Indices::pressureIdx] += analyticalSolutionCellCenter * analyticalSolutionCellCenter * scv.volume(); totalVolume += scv.volume(); // treat face dofs for (auto&& scvf : scvfs(fvGeometry)) // integrate over element using a quadrature rule const auto& rule = Dune::QuadratureRules::rule(element.geometry().type(), pOrder); for (typename Dune::QuadratureRule::const_iterator qIt = rule.begin(); qIt != rule.end(); ++qIt) { const int dofIdxFace = scvf.dofIndex(); const int dirIdx = scvf.directionIndex(); const auto analyticalSolutionFace = problem.dirichletAtPos(scvf.center())[Indices::velocity(dirIdx)]; const auto numericalSolutionFace = curSol[FVGridGeometry::faceIdx()][dofIdxFace][0]; 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; auto localPos = qIt->position(); const auto dofIdxCellCenter = scv.dofIndex(); const auto analyticalSolutionCellCenter = problem.dirichletAtPos(element.geometry().global(localPos))[Indices::pressureIdx]; const auto numericalSolutionCellCenter = curSol[FVGridGeometry::cellCenterIdx()][dofIdxCellCenter][Indices::pressureIdx - ModelTraits::dim()]; Scalar weigth = qIt->weight() * element.geometry().integrationElement(localPos); sumError[Indices::pressureIdx] += squaredDiff_(analyticalSolutionCellCenter, numericalSolutionCellCenter) * weigth; sumReference[Indices::pressureIdx] += analyticalSolutionCellCenter * analyticalSolutionCellCenter * weigth; totalVolume[Indices::pressureIdx] += weigth; // treat face dofs for (auto&& scvf : scvfs(fvGeometry)) { const int dofIdxFace = scvf.dofIndex(); const int dirIdx = scvf.directionIndex(); // only treat the part of the quadrature which belongs to this dof if (((scvf.localFaceIdx() % 2 == 0) && localPos[dirIdx] > 0.5) || ((scvf.localFaceIdx() % 2 == 1) && localPos[dirIdx] < 0.5 + 1e-10)) continue; const auto analyticalSolutionFace = problem.dirichletAtPos(element.geometry().global(localPos))[Indices::velocity(dirIdx)]; const auto numericalSolutionFace = curSol[FVGridGeometry::faceIdx()][dofIdxFace][0]; directionIndex[dofIdxFace] = dirIdx; errorVelocity[dofIdxFace] += squaredDiff_(analyticalSolutionFace, numericalSolutionFace) * weigth; velocityReference[dofIdxFace] += squaredDiff_(analyticalSolutionFace, 0.0) * weigth; staggeredVolume[dofIdxFace] += weigth; } } } } // get the absolute and relative discrete L2-error for cell-center dofs l2NormAbs[Indices::pressureIdx] = std::sqrt(sumError[Indices::pressureIdx] / totalVolume); l2NormAbs[Indices::pressureIdx] = std::sqrt(sumError[Indices::pressureIdx] / totalVolume[Indices::pressureIdx]); l2NormRel[Indices::pressureIdx] = std::sqrt(sumError[Indices::pressureIdx] / sumReference[Indices::pressureIdx]); // 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[Indices::velocity(dirIdx)] += error * volume; sumReference[Indices::velocity(dirIdx)] += ref * volume; if (verbose) { sumError[Indices::velocity(directionIndex[i])] += error; sumReference[Indices::velocity(directionIndex[i])] += ref; totalVolume[Indices::velocity(directionIndex[i])] += staggeredVolume[i]; } else { sumError[Indices::velocity(0)] += error; sumReference[Indices::velocity(0)] += ref; totalVolume[Indices::velocity(0)] += staggeredVolume[i]; } } for(int dirIdx = 0; dirIdx < ModelTraits::dim(); ++dirIdx) { l2NormAbs[Indices::velocity(dirIdx)] = std::sqrt(sumError[Indices::velocity(dirIdx)] / totalVolume); l2NormAbs[Indices::velocity(dirIdx)] = std::sqrt(sumError[Indices::velocity(dirIdx)] / totalVolume[Indices::velocity(dirIdx)]); l2NormRel[Indices::velocity(dirIdx)] = std::sqrt(sumError[Indices::velocity(dirIdx)] / sumReference[Indices::velocity(dirIdx)]); } return std::make_pair(l2NormAbs, l2NormRel); ... ... @@ -125,6 +174,6 @@ private: } }; } //end namespace DUMUX_TEST_L2_ERROR_HH } #endif
 ... ... @@ -30,7 +30,7 @@ #include #include #include #include "l2error.hh" #include namespace Dumux { ... ... @@ -117,17 +117,8 @@ public: { if(printL2Error_) { using L2Error = NavierStokesTestL2Error; const auto l2error = L2Error::calculateL2Error(*this, curSol); const int numCellCenterDofs = this->fvGridGeometry().numCellCenterDofs(); const int numFaceDofs = this->fvGridGeometry().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[Indices::pressureIdx] << " / " << l2error.second[Indices::pressureIdx] << ", L2(vx) = " << l2error.first[Indices::velocityXIdx] << " / " << l2error.second[Indices::velocityXIdx] << ", L2(vy) = " << l2error.first[Indices::velocityYIdx] << " / " << l2error.second[Indices::velocityYIdx] << std::endl; using L2Error = NavierStokesTestL2Error; L2Error::printL2Error(*this, curSol); } } ... ...
 ... ... @@ -28,16 +28,19 @@ done grep "L2 error (abs/rel) for" \$LOGFILE | tee \$L2ERRORFILE echo "reset; \ set log x; \ set xlabel 'refinement'; \ set log y; \ set arrow from graph 0,1 to graph 1,0 nohead lc rgb 'gray'; \ set arrow from graph 0,1 to graph 1,0.5 nohead lc rgb 'gray'" > \$1.gp set ylabel 'l2-error'; \ set arrow from 5,0.04 to 5,0.08 nohead lc 8; \ set arrow from 4,0.08 to 5,0.08 nohead lc 8; \ set arrow from 4,0.08 to 5,0.04 nohead lc 8" > \$1.gp PLOT="plot '\$L2ERRORFILE' u 6:17 w lp t 'pressure', '\$L2ERRORFILE' u 6:23 w lp t 'velocity'" USE_REL_ERROR=2 # set to "2" if the use relative error PLOT="plot '\$L2ERRORFILE' u :`expr 17 + \$USE_REL_ERROR` w lp t 'pressure', '\$L2ERRORFILE' u :`expr 23 + \$USE_REL_ERROR` w lp t 'velocity'" if [ \$2 == 2 ]; then PLOT=\$PLOT", '\$L2ERRORFILE' u 6:29 w lp t 'velocity'" PLOT=\$PLOT", '\$L2ERRORFILE' u :`expr 29 + \$USE_REL_ERROR` w lp t 'velocity'" elif [ \$2 == 3 ]; then PLOT=\$PLOT", '\$L2ERRORFILE' u 6:35 w lp t 'velocity'" PLOT=\$PLOT", '\$L2ERRORFILE' u :`expr 35 + \$USE_REL_ERROR` w lp t 'velocity'" fi echo \$PLOT >> \$1.gp ... ...
 ... ... @@ -30,7 +30,7 @@ #include #include #include #include "l2error.hh" #include namespace Dumux ... ... @@ -97,6 +97,10 @@ public: { printL2Error_ = getParam("Problem.PrintL2Error"); createAnalyticalSolution_(); using CellArray = std::array; const auto numCells = getParam("Grid.Cells"); cellSizeX_ = this->fvGridGeometry().bBoxMax()[0] / numCells[0]; } /*! ... ... @@ -113,17 +117,8 @@ public: { if(printL2Error_) { using L2Error = NavierStokesTestL2Error; const auto l2error = L2Error::calculateL2Error(*this, curSol); const int numCellCenterDofs = this->fvGridGeometry().numCellCenterDofs(); const int numFaceDofs = this->fvGridGeometry().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[Indices::pressureIdx] << " / " << l2error.second[Indices::pressureIdx] << " , L2(vx) = " << l2error.first[Indices::velocityXIdx] << " / " << l2error.second[Indices::velocityXIdx] << " , L2(vy) = " << l2error.first[Indices::velocityYIdx] << " / " << l2error.second[Indices::velocityYIdx] << std::endl; using L2Error = NavierStokesTestL2Error; L2Error::printL2Error(*this, curSol); } } ... ... @@ -174,7 +169,12 @@ public: // set Dirichlet values for the velocity and pressure everywhere values.setDirichlet(Indices::momentumXBalanceIdx); values.setDirichlet(Indices::momentumYBalanceIdx); values.setDirichletCell(Indices::conti0EqIdx); // set a fixed pressure in one cell if (isLowerLeftCell_(globalPos)) values.setDirichletCell(Indices::conti0EqIdx); else values.setOutflow(Indices::conti0EqIdx); return values; } ... ... @@ -291,9 +291,15 @@ private: analyticalVelocity_[ccDofIdx][dirIdx] = analyticalSolutionAtCc[Indices::velocity(dirIdx)]; } } } } bool isLowerLeftCell_(const GlobalPosition& globalPos) const { return globalPos[0] < (this->fvGridGeometry().bBoxMin()[0] + 0.5 * cellSizeX_ + eps_); } Scalar eps_; Scalar cellSizeX_; bool printL2Error_; std::vector analyticalPressure_; std::vector analyticalVelocity_; ... ...
 ... ... @@ -30,7 +30,7 @@ #include #include #include #include "l2error.hh" #include namespace Dumux { ... ... @@ -119,17 +119,8 @@ public: { if(printL2Error_) { using L2Error = NavierStokesTestL2Error; const auto l2error = L2Error::calculateL2Error(*this, curSol); const int numCellCenterDofs = this->fvGridGeometry().numCellCenterDofs(); const int numFaceDofs = this->fvGridGeometry().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[Indices::pressureIdx] << " / " << l2error.second[Indices::pressureIdx] << " , L2(vx) = " << l2error.first[Indices::velocityXIdx] << " / " << l2error.second[Indices::velocityXIdx] << " , L2(vy) = " << l2error.first[Indices::velocityYIdx] << " / " << l2error.second[Indices::velocityYIdx] << std::endl; using L2Error = NavierStokesTestL2Error; L2Error::printL2Error(*this, curSol); } } ... ...
 ... ... @@ -32,7 +32,7 @@ #include #include #include #include "l2error.hh" #include namespace Dumux ... ... @@ -117,16 +117,8 @@ public: { if(printL2Error_) { using L2Error = NavierStokesTestL2Error; const auto l2error = L2Error::calculateL2Error(*this, curSol); const int numCellCenterDofs = this->fvGridGeometry().numCellCenterDofs(); const int numFaceDofs = this->fvGridGeometry().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[Indices::pressureIdx] << " / " << l2error.second[Indices::pressureIdx] << " , L2(vx) = " << l2error.first[Indices::velocityXIdx] << " / " << l2error.second[Indices::velocityXIdx] << std::endl; using L2Error = NavierStokesTestL2Error; L2Error::printL2Error(*this, curSol); } } ... ...
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!