| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // vi: set et ts=4 sw=4 sts=4: | ||
| 3 | // | ||
| 4 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup NavierStokesTests | ||
| 10 | * \copydoc Dumux::NavierStokesTestError | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_TEST_FREEFLOW_NAVIERSTOKES_ERRORS_HH | ||
| 13 | #define DUMUX_TEST_FREEFLOW_NAVIERSTOKES_ERRORS_HH | ||
| 14 | |||
| 15 | #include <vector> | ||
| 16 | #include <cmath> | ||
| 17 | #include <fstream> | ||
| 18 | |||
| 19 | #include <dune/common/fvector.hh> | ||
| 20 | #include <dune/common/indices.hh> | ||
| 21 | |||
| 22 | #include <dumux/discretization/method.hh> | ||
| 23 | #include <dumux/discretization/extrusion.hh> | ||
| 24 | #include <dumux/io/format.hh> | ||
| 25 | #include <dumux/geometry/diameter.hh> | ||
| 26 | |||
| 27 | namespace Dumux { | ||
| 28 | |||
| 29 | /*! | ||
| 30 | * \brief Compute errors between an analytical solution and the numerical approximation | ||
| 31 | */ | ||
| 32 | template<class Problem, class Scalar = double> | ||
| 33 |
1/4✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
12 | class NavierStokesErrors |
| 34 | { | ||
| 35 | using GridGeometry = std::decay_t<decltype(std::declval<Problem>().gridGeometry())>; | ||
| 36 | using SubControlVolume = typename GridGeometry::SubControlVolume; | ||
| 37 | using Extrusion = Extrusion_t<GridGeometry>; | ||
| 38 | using Indices = typename Problem::Indices; | ||
| 39 | using Element = typename GridGeometry::LocalView::Element; | ||
| 40 | using PrimaryVariables = std::decay_t<decltype( | ||
| 41 | std::declval<Problem>().dirichlet(std::declval<Element>(), std::declval<SubControlVolume>()) | ||
| 42 | )>; | ||
| 43 | |||
| 44 | static constexpr int dim = GridGeometry::GridView::dimension; | ||
| 45 | |||
| 46 | public: | ||
| 47 | template<class SolutionVector> | ||
| 48 | 6 | NavierStokesErrors(std::shared_ptr<const Problem> problem, | |
| 49 | const SolutionVector& curSol, | ||
| 50 | Scalar time = 0.0) | ||
| 51 |
4/10✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 6 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 6 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 6 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
|
6 | : problem_(problem) |
| 52 |
1/2✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
|
6 | { calculateErrors_(curSol, time); } |
| 53 | |||
| 54 | /*! | ||
| 55 | * \brief Computes errors between an analytical solution and the numerical approximation | ||
| 56 | * | ||
| 57 | * \param curSol The current solution vector | ||
| 58 | * \param time The current time | ||
| 59 | */ | ||
| 60 | template<class SolutionVector> | ||
| 61 | void update(const SolutionVector& curSol, Scalar time = 0.0) | ||
| 62 | { calculateErrors_(curSol, time); } | ||
| 63 | |||
| 64 | //! The (absolute) discrete l2 error | ||
| 65 | ✗ | const PrimaryVariables& l2Absolute() const { return l2Absolute_; } | |
| 66 | //! The relative discrete l2 error (relative to the discrete l2 norm of the reference solution) | ||
| 67 | ✗ | const PrimaryVariables& l2Relative() const { return l2Relative_; } | |
| 68 | //! The (absolute) discrete l-infinity error | ||
| 69 | ✗ | const PrimaryVariables& lInfAbsolute() const { return lInfAbsolute_; } | |
| 70 | //! The relative discrete l-infinity error (relative to the discrete loo norm of the reference solution) | ||
| 71 | ✗ | const PrimaryVariables& lInfRelative() const { return lInfRelative_; } | |
| 72 | |||
| 73 | //! Time corresponding to the error (returns 0 per default) | ||
| 74 | ✗ | Scalar time() const { return time_; } | |
| 75 | |||
| 76 | private: | ||
| 77 | template<class SolutionVector> | ||
| 78 | 6 | void calculateErrors_(const SolutionVector& curSol, Scalar time) | |
| 79 | { | ||
| 80 | // store time information | ||
| 81 | 6 | time_ = time; | |
| 82 | |||
| 83 | // calculate helping variables | ||
| 84 | 6 | Scalar totalVolume = 0.0; | |
| 85 | |||
| 86 | 6 | PrimaryVariables sumReference(0.0); | |
| 87 | 6 | PrimaryVariables sumError(0.0); | |
| 88 | 6 | PrimaryVariables maxReference(0.0); | |
| 89 | 6 | PrimaryVariables maxError(0.0); | |
| 90 | |||
| 91 | 18 | auto fvGeometry = localView(problem_->gridGeometry()); | |
| 92 |
1/2✓ Branch 4 taken 67206 times.
✗ Branch 5 not taken.
|
134424 | for (const auto& element : elements(problem_->gridGeometry().gridView())) |
| 93 | { | ||
| 94 | 67200 | fvGeometry.bindElement(element); | |
| 95 | |||
| 96 |
4/4✓ Branch 2 taken 67200 times.
✓ Branch 3 taken 67200 times.
✓ Branch 4 taken 67200 times.
✓ Branch 5 taken 67200 times.
|
201600 | for (const auto& scv : scvs(fvGeometry)) |
| 97 | { | ||
| 98 | 67200 | totalVolume += Extrusion::volume(fvGeometry, scv); | |
| 99 | |||
| 100 | // compute the pressure errors | ||
| 101 | 67200 | const auto analyticalSolutionCellCenter | |
| 102 |
2/2✓ Branch 3 taken 1826 times.
✓ Branch 4 taken 65374 times.
|
201600 | = problem_->analyticalSolution(scv.dofPosition(), time)[Indices::pressureIdx]; |
| 103 | 67200 | const auto numericalSolutionCellCenter | |
| 104 |
6/6✓ Branch 0 taken 1826 times.
✓ Branch 1 taken 65374 times.
✓ Branch 2 taken 1826 times.
✓ Branch 3 taken 65374 times.
✓ Branch 4 taken 1826 times.
✓ Branch 5 taken 65374 times.
|
201600 | = curSol[GridGeometry::cellCenterIdx()][scv.dofIndex()][Indices::pressureIdx - dim]; |
| 105 | |||
| 106 | 134400 | const Scalar pError = absDiff_(analyticalSolutionCellCenter, numericalSolutionCellCenter); | |
| 107 | 134400 | const Scalar pReference = absDiff_(analyticalSolutionCellCenter, 0.0); | |
| 108 | |||
| 109 |
4/4✓ Branch 0 taken 1826 times.
✓ Branch 1 taken 65374 times.
✓ Branch 2 taken 1799 times.
✓ Branch 3 taken 65401 times.
|
69026 | maxError[Indices::pressureIdx] = std::max(maxError[Indices::pressureIdx], pError); |
| 110 |
4/4✓ Branch 0 taken 1799 times.
✓ Branch 1 taken 65401 times.
✓ Branch 2 taken 1799 times.
✓ Branch 3 taken 65401 times.
|
136199 | maxReference[Indices::pressureIdx] = std::max(maxReference[Indices::pressureIdx], pReference); |
| 111 | 134400 | sumError[Indices::pressureIdx] += pError * pError * Extrusion::volume(fvGeometry, scv); | |
| 112 | 134400 | sumReference[Indices::pressureIdx] += pReference * pReference * Extrusion::volume(fvGeometry, scv); | |
| 113 | |||
| 114 |
4/4✓ Branch 0 taken 268800 times.
✓ Branch 1 taken 67200 times.
✓ Branch 2 taken 268800 times.
✓ Branch 3 taken 67200 times.
|
672000 | for (const auto& scvf : scvfs(fvGeometry)) |
| 115 | { | ||
| 116 | // compute the velocity errors | ||
| 117 | 268800 | const auto velocityIndex = Indices::velocity(scvf.directionIndex()); | |
| 118 | |||
| 119 | 268800 | const Scalar staggeredHalfVolume = Extrusion::volume(fvGeometry, | |
| 120 | typename GridGeometry::Traits::FaceSubControlVolume( | ||
| 121 | 1612800 | 0.5*(scv.center() + scvf.center()), 0.5*Extrusion::volume(fvGeometry, scv) | |
| 122 | ) | ||
| 123 | ); | ||
| 124 | |||
| 125 |
2/2✓ Branch 3 taken 1512 times.
✓ Branch 4 taken 267288 times.
|
806400 | const auto analyticalSolutionFace = problem_->analyticalSolution(scvf.center(), time)[velocityIndex]; |
| 126 |
6/6✓ Branch 0 taken 1512 times.
✓ Branch 1 taken 267288 times.
✓ Branch 2 taken 1512 times.
✓ Branch 3 taken 267288 times.
✓ Branch 4 taken 1512 times.
✓ Branch 5 taken 267288 times.
|
806400 | const auto numericalSolutionFace = curSol[GridGeometry::faceIdx()][scvf.dofIndex()][0]; |
| 127 | |||
| 128 | 537600 | const Scalar velError = absDiff_(analyticalSolutionFace, numericalSolutionFace); | |
| 129 | 537600 | const Scalar velReference = absDiff_(analyticalSolutionFace, 0.0); | |
| 130 | |||
| 131 |
4/4✓ Branch 0 taken 1512 times.
✓ Branch 1 taken 267288 times.
✓ Branch 2 taken 7701 times.
✓ Branch 3 taken 261099 times.
|
270312 | maxError[velocityIndex] = std::max(maxError[velocityIndex], velError); |
| 132 |
4/4✓ Branch 0 taken 7701 times.
✓ Branch 1 taken 261099 times.
✓ Branch 2 taken 7701 times.
✓ Branch 3 taken 261099 times.
|
545301 | maxReference[velocityIndex] = std::max(maxReference[velocityIndex], velReference); |
| 133 | 268800 | sumError[velocityIndex] += velError * velError * staggeredHalfVolume; | |
| 134 | 537600 | sumReference[velocityIndex] += velReference * velReference * staggeredHalfVolume; | |
| 135 | } | ||
| 136 | } | ||
| 137 | } | ||
| 138 | |||
| 139 | // calculate errors | ||
| 140 |
2/2✓ Branch 0 taken 18 times.
✓ Branch 1 taken 6 times.
|
24 | for (int i = 0; i < PrimaryVariables::size(); ++i) |
| 141 | { | ||
| 142 | 36 | l2Absolute_[i] = std::sqrt(sumError[i] / totalVolume); | |
| 143 | 54 | l2Relative_[i] = std::sqrt(sumError[i] / sumReference[i]); | |
| 144 | 36 | lInfAbsolute_[i] = maxError[i]; | |
| 145 | 72 | lInfRelative_[i] = maxError[i] / maxReference[i]; | |
| 146 | } | ||
| 147 | 6 | } | |
| 148 | |||
| 149 | template<class T> | ||
| 150 | ✗ | T absDiff_(const T& a, const T& b) const | |
| 151 |
12/12✓ Branch 0 taken 1826 times.
✓ Branch 1 taken 65374 times.
✓ Branch 2 taken 1826 times.
✓ Branch 3 taken 65374 times.
✓ Branch 4 taken 1826 times.
✓ Branch 5 taken 65374 times.
✓ Branch 6 taken 1512 times.
✓ Branch 7 taken 267288 times.
✓ Branch 8 taken 1512 times.
✓ Branch 9 taken 267288 times.
✓ Branch 10 taken 1512 times.
✓ Branch 11 taken 267288 times.
|
1008000 | { using std::abs; return abs(a-b); } |
| 152 | |||
| 153 | std::shared_ptr<const Problem> problem_; | ||
| 154 | |||
| 155 | PrimaryVariables l2Absolute_; | ||
| 156 | PrimaryVariables l2Relative_; | ||
| 157 | PrimaryVariables lInfAbsolute_; | ||
| 158 | PrimaryVariables lInfRelative_; | ||
| 159 | Scalar time_; | ||
| 160 | }; | ||
| 161 | |||
| 162 | template<class Problem, class SolutionVector> | ||
| 163 | NavierStokesErrors(std::shared_ptr<Problem>, SolutionVector&&) | ||
| 164 | -> NavierStokesErrors<Problem>; | ||
| 165 | |||
| 166 | template<class Problem, class SolutionVector, class Scalar> | ||
| 167 | NavierStokesErrors(std::shared_ptr<Problem>, SolutionVector&&, Scalar) | ||
| 168 | -> NavierStokesErrors<Problem, Scalar>; | ||
| 169 | |||
| 170 | |||
| 171 | /*! | ||
| 172 | * \brief An error CSV file writer | ||
| 173 | */ | ||
| 174 | template<class Problem, class Scalar = double> | ||
| 175 | ✗ | class NavierStokesErrorCSVWriter | |
| 176 | { | ||
| 177 | using GridGeometry = std::decay_t<decltype(std::declval<Problem>().gridGeometry())>; | ||
| 178 | using Indices = typename Problem::Indices; | ||
| 179 | |||
| 180 | static constexpr int dim = GridGeometry::GridView::dimension; | ||
| 181 | |||
| 182 | public: | ||
| 183 | ✗ | NavierStokesErrorCSVWriter(std::shared_ptr<const Problem> problem, const std::string& suffix = "") | |
| 184 | ✗ | : name_(problem->name() + (suffix.empty() ? "_error" : "_error_" + suffix)) | |
| 185 | { | ||
| 186 | ✗ | const int numCCDofs = problem->gridGeometry().numCellCenterDofs(); | |
| 187 | ✗ | const int numFaceDofs = problem->gridGeometry().numFaceDofs(); | |
| 188 | |||
| 189 | // print auxiliary file with the number of dofs | ||
| 190 | ✗ | std::ofstream logFileDofs(name_ + "_dofs.csv", std::ios::trunc); | |
| 191 | ✗ | logFileDofs << "cc dofs, face dofs, all dofs\n" | |
| 192 | ✗ | << numCCDofs << ", " << numFaceDofs << ", " << numCCDofs + numFaceDofs << "\n"; | |
| 193 | |||
| 194 | // clear error file | ||
| 195 | ✗ | std::ofstream logFile(name_ + ".csv", std::ios::trunc); | |
| 196 | |||
| 197 | // write header | ||
| 198 | ✗ | logFile << "time"; | |
| 199 | ✗ | for (const std::string e : { "L2Abs", "L2Rel", "LinfAbs", "LinfRel" }) | |
| 200 | ✗ | printError_(logFile, errorNames_(e), "{:s}"); | |
| 201 | ✗ | logFile << "\n"; | |
| 202 | ✗ | } | |
| 203 | |||
| 204 | ✗ | void printErrors(const NavierStokesErrors<Problem>& errors) const | |
| 205 | { | ||
| 206 | // append to error file | ||
| 207 | ✗ | std::ofstream logFile(name_ + ".csv", std::ios::app); | |
| 208 | |||
| 209 | ✗ | logFile << Fmt::format("{:.5e}", errors.time()); | |
| 210 | ✗ | printError_(logFile, errors.l2Absolute()); | |
| 211 | ✗ | printError_(logFile, errors.l2Relative()); | |
| 212 | ✗ | printError_(logFile, errors.lInfAbsolute()); | |
| 213 | ✗ | printError_(logFile, errors.lInfRelative()); | |
| 214 | |||
| 215 | ✗ | logFile << "\n"; | |
| 216 | ✗ | } | |
| 217 | |||
| 218 | private: | ||
| 219 | template<class Error> | ||
| 220 | ✗ | void printError_(std::ofstream& logFile, const Error& error, const std::string& format = "{:.5e}") const | |
| 221 | { | ||
| 222 | ✗ | logFile << Fmt::vformat(", " + format, Fmt::make_format_args(error[Indices::pressureIdx])); | |
| 223 | ✗ | for (int dirIdx = 0; dirIdx < dim; ++dirIdx) | |
| 224 | ✗ | logFile << Fmt::vformat(", " + format, Fmt::make_format_args(error[Indices::velocity(dirIdx)])); | |
| 225 | ✗ | } | |
| 226 | |||
| 227 | ✗ | std::vector<std::string> errorNames_(const std::string& e) const | |
| 228 | { | ||
| 229 | if constexpr (dim == 1) | ||
| 230 | return { e + "(u)", e + "(p)" }; | ||
| 231 | else if constexpr (dim == 2) | ||
| 232 | ✗ | return { e + "(u)", e + "(v)", e + "(p)" }; | |
| 233 | else | ||
| 234 | return { e + "(u)", e + "(v)", e + "(w)", e + "(p)" }; | ||
| 235 | } | ||
| 236 | |||
| 237 | std::string name_; | ||
| 238 | }; | ||
| 239 | |||
| 240 | template<class Problem> | ||
| 241 | NavierStokesErrorCSVWriter(std::shared_ptr<Problem>) | ||
| 242 | -> NavierStokesErrorCSVWriter<Problem>; | ||
| 243 | |||
| 244 | template<class Problem> | ||
| 245 | NavierStokesErrorCSVWriter(std::shared_ptr<Problem>, const std::string&) | ||
| 246 | -> NavierStokesErrorCSVWriter<Problem>; | ||
| 247 | |||
| 248 | |||
| 249 | /*! | ||
| 250 | * \brief Append errors to the log file during a convergence test | ||
| 251 | */ | ||
| 252 | template<class Problem> | ||
| 253 | ✗ | void convergenceTestAppendErrors(std::ofstream& logFile, | |
| 254 | std::shared_ptr<Problem> problem, | ||
| 255 | const NavierStokesErrors<Problem>& errors) | ||
| 256 | { | ||
| 257 | ✗ | const auto numCCDofs = problem->gridGeometry().numCellCenterDofs(); | |
| 258 | ✗ | const auto numFaceDofs = problem->gridGeometry().numFaceDofs(); | |
| 259 | |||
| 260 | ✗ | logFile << Fmt::format( | |
| 261 | "[ConvergenceTest] numCCDofs = {} numFaceDofs = {}", | ||
| 262 | numCCDofs, numFaceDofs | ||
| 263 | ); | ||
| 264 | |||
| 265 | ✗ | const auto print = [&](const auto& e, const std::string& name){ | |
| 266 | using Indices = typename Problem::Indices; | ||
| 267 | ✗ | logFile << Fmt::format( | |
| 268 | " {0}(p) = {1} {0}(vx) = {2} {0}(vy) = {3}", | ||
| 269 | name, e[Indices::pressureIdx], e[Indices::velocityXIdx], e[Indices::velocityYIdx] | ||
| 270 | ); | ||
| 271 | }; | ||
| 272 | |||
| 273 | ✗ | print(errors.l2Absolute(), "L2Abs"); | |
| 274 | ✗ | print(errors.l2Relative(), "L2Rel"); | |
| 275 | ✗ | print(errors.lInfAbsolute(), "LinfAbs"); | |
| 276 | ✗ | print(errors.lInfRelative(), "LinfRel"); | |
| 277 | |||
| 278 | ✗ | logFile << "\n"; | |
| 279 | ✗ | } | |
| 280 | |||
| 281 | /*! | ||
| 282 | * \brief Append errors to the log file during a convergence test | ||
| 283 | */ | ||
| 284 | template<class Problem> | ||
| 285 | ✗ | void convergenceTestAppendErrors(std::shared_ptr<Problem> problem, | |
| 286 | const NavierStokesErrors<Problem>& errors) | ||
| 287 | { | ||
| 288 | ✗ | const auto logFileName = problem->name() + ".log"; | |
| 289 | ✗ | std::ofstream logFile(logFileName, std::ios::app); | |
| 290 | ✗ | convergenceTestAppendErrors(logFile, problem, errors); | |
| 291 | ✗ | } | |
| 292 | |||
| 293 | namespace NavierStokesTest { | ||
| 294 | |||
| 295 | /*! | ||
| 296 | * \brief Compute errors between an analytical solution and the numerical approximation for momentum eq. | ||
| 297 | */ | ||
| 298 | template<class Problem, class Scalar = double> | ||
| 299 |
1/5✓ Branch 0 taken 23 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
46 | class ErrorsSubProblem |
| 300 | { | ||
| 301 | static constexpr bool isMomentumProblem = Problem::isMomentumProblem(); | ||
| 302 | static constexpr int dim | ||
| 303 | = std::decay_t<decltype(std::declval<Problem>().gridGeometry())>::GridView::dimension; | ||
| 304 | |||
| 305 | using ErrorVector = typename std::conditional_t< isMomentumProblem, | ||
| 306 | Dune::FieldVector<Scalar, dim>, | ||
| 307 | Dune::FieldVector<Scalar, 1> >; | ||
| 308 | |||
| 309 | public: | ||
| 310 | template<class SolutionVector> | ||
| 311 | 151 | ErrorsSubProblem(std::shared_ptr<const Problem> problem, | |
| 312 | const SolutionVector& curSol, | ||
| 313 | Scalar time = 0.0) | ||
| 314 |
4/10✓ Branch 2 taken 87 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 87 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 87 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 87 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
|
151 | : problem_(problem) |
| 315 |
1/2✓ Branch 1 taken 87 times.
✗ Branch 2 not taken.
|
151 | { calculateErrors_(curSol, time); } |
| 316 | |||
| 317 | /*! | ||
| 318 | * \brief Computes errors between an analytical solution and the numerical approximation | ||
| 319 | * | ||
| 320 | * \param curSol The current solution vector | ||
| 321 | * \param time The current time | ||
| 322 | */ | ||
| 323 | template<class SolutionVector> | ||
| 324 | void update(const SolutionVector& curSol, Scalar time = 0.0) | ||
| 325 | 138 | { calculateErrors_(curSol, time); } | |
| 326 | |||
| 327 | //! The (absolute) discrete l2 error | ||
| 328 | 144 | const ErrorVector& l2Absolute() const { return l2Absolute_; } | |
| 329 | //! The relative discrete l2 error (relative to the discrete l2 norm of the reference solution) | ||
| 330 | 144 | const ErrorVector& l2Relative() const { return l2Relative_; } | |
| 331 | //! The (absolute) discrete l-infinity error | ||
| 332 | 144 | const ErrorVector& lInfAbsolute() const { return lInfAbsolute_; } | |
| 333 | //! The relative discrete l-infinity error (relative to the discrete loo norm of the reference solution) | ||
| 334 | 144 | const ErrorVector& lInfRelative() const { return lInfRelative_; } | |
| 335 | |||
| 336 | //! Time corresponding to the error (returns 0 per default) | ||
| 337 | ✗ | Scalar time() const { return time_; } | |
| 338 | |||
| 339 | //! Maximum diameter of primal grid elements | ||
| 340 | ✗ | Scalar hMax() const { return hMax_; } | |
| 341 | |||
| 342 | //! Volume of domain | ||
| 343 | ✗ | Scalar totalVolume() const { return totalVolume_; } | |
| 344 | |||
| 345 | //! Number of scvs that have been considered in error calculation | ||
| 346 | 155 | Scalar numDofs() const { return numDofs_; } | |
| 347 | |||
| 348 | private: | ||
| 349 | template<class SolutionVector> | ||
| 350 | 427 | void calculateErrors_(const SolutionVector& curSol, Scalar time) | |
| 351 | { | ||
| 352 | // store time information | ||
| 353 | 427 | time_ = time; | |
| 354 | |||
| 355 | // calculate helping variables | ||
| 356 | 427 | totalVolume_ = 0.0; | |
| 357 | 427 | hMax_ = 0.0; | |
| 358 |
3/6✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
|
1281 | numDofs_ = problem_->gridGeometry().numDofs(); |
| 359 | |||
| 360 | using GridGeometry = std::decay_t<decltype(std::declval<Problem>().gridGeometry())>; | ||
| 361 | // We do not consider the overlapping Dofs, i.e. elements, in the errors | ||
| 362 | if constexpr (GridGeometry::discMethod == DiscretizationMethods::pq1bubble) | ||
| 363 | 20 | numDofs_ -= problem_->gridGeometry().gridView().size(0); | |
| 364 | |||
| 365 |
1/2✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
|
427 | ErrorVector sumReference(0.0); |
| 366 |
1/2✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
|
427 | ErrorVector sumError(0.0); |
| 367 |
1/2✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
|
427 | ErrorVector maxReference(0.0); |
| 368 |
1/2✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
|
427 | ErrorVector maxError(0.0); |
| 369 | |||
| 370 | using namespace Dune::Indices; | ||
| 371 | |||
| 372 |
4/10✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 28 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 22 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
|
1325 | auto fvGeometry = localView(problem_->gridGeometry()); |
| 373 |
17/24✓ Branch 1 taken 28 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1636369 times.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 28 times.
✓ Branch 8 taken 12 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 7100 times.
✓ Branch 11 taken 12 times.
✓ Branch 12 taken 7072 times.
✓ Branch 13 taken 40 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 7072 times.
✓ Branch 16 taken 28 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 42022 times.
✓ Branch 19 taken 28 times.
✓ Branch 20 taken 34950 times.
✓ Branch 21 taken 28 times.
✓ Branch 23 taken 34950 times.
✗ Branch 24 not taken.
✓ Branch 26 taken 34950 times.
✗ Branch 27 not taken.
|
6662638 | for (const auto& element : elements(problem_->gridGeometry().gridView())) |
| 374 | { | ||
| 375 |
5/5✓ Branch 1 taken 42022 times.
✓ Branch 2 taken 1045 times.
✓ Branch 3 taken 1635111 times.
✓ Branch 4 taken 135 times.
✓ Branch 5 taken 41887 times.
|
3296145 | hMax_ = std::max(hMax_, diameter(element.geometry())); |
| 376 |
1/2✓ Branch 1 taken 42022 times.
✗ Branch 2 not taken.
|
3295914 | fvGeometry.bindElement(element); |
| 377 |
8/8✓ Branch 0 taken 303690 times.
✓ Branch 1 taken 3289478 times.
✓ Branch 2 taken 1883702 times.
✓ Branch 3 taken 4077802 times.
✓ Branch 4 taken 1597730 times.
✓ Branch 5 taken 798894 times.
✓ Branch 6 taken 2177 times.
✓ Branch 7 taken 514419 times.
|
21375874 | for (const auto& scv : scvs(fvGeometry)) |
| 378 | { | ||
| 379 | using Extrusion = Extrusion_t<GridGeometry>; | ||
| 380 | |||
| 381 | // velocity errors | ||
| 382 | if constexpr (isMomentumProblem) | ||
| 383 | { | ||
| 384 | if constexpr (GridGeometry::discMethod == DiscretizationMethods::fcstaggered) | ||
| 385 | { | ||
| 386 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 189 times.
|
6432816 | totalVolume_ += Extrusion::volume(fvGeometry, scv); |
| 387 | // compute the velocity errors | ||
| 388 | using Indices = typename Problem::Indices; | ||
| 389 | 6432816 | const auto velIdx = Indices::velocity(scv.dofAxis()); | |
| 390 | 6432816 | const auto analyticalSolution | |
| 391 |
6/8✓ Branch 0 taken 3 times.
✓ Branch 1 taken 189 times.
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 19323 times.
✓ Branch 4 taken 3203813 times.
✓ Branch 5 taken 189 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
19298448 | = problem_->analyticalSolution(scv.dofPosition(), time)[velIdx]; |
| 392 | 6432816 | const auto numericalSolution | |
| 393 |
2/2✓ Branch 0 taken 19137 times.
✓ Branch 1 taken 3203999 times.
|
6432816 | = curSol[scv.dofIndex()][0]; |
| 394 | |||
| 395 | 12865632 | const Scalar vError = absDiff_(analyticalSolution, numericalSolution); | |
| 396 | 12865632 | const Scalar vReference = absDiff_(analyticalSolution, 0.0); | |
| 397 | |||
| 398 |
4/4✓ Branch 0 taken 19134 times.
✓ Branch 1 taken 3203810 times.
✓ Branch 2 taken 25080 times.
✓ Branch 3 taken 3197864 times.
|
6471176 | maxError[velIdx] = std::max(maxError[velIdx], vError); |
| 399 |
4/4✓ Branch 0 taken 25176 times.
✓ Branch 1 taken 3197960 times.
✓ Branch 2 taken 25080 times.
✓ Branch 3 taken 3197864 times.
|
12914836 | maxReference[velIdx] = std::max(maxReference[velIdx], vReference); |
| 400 | 12865248 | sumError[velIdx] += vError * vError * Extrusion::volume(fvGeometry, scv); | |
| 401 | 19298064 | sumReference[velIdx] += vReference * vReference * Extrusion::volume(fvGeometry, scv); | |
| 402 | } | ||
| 403 | else if constexpr (GridGeometry::discMethod == DiscretizationMethods::fcdiamond | ||
| 404 | || GridGeometry::discMethod == DiscretizationMethods::box) | ||
| 405 | { | ||
| 406 | 156674 | totalVolume_ += Extrusion::volume(fvGeometry, scv); | |
| 407 |
2/2✓ Branch 0 taken 274948 times.
✓ Branch 1 taken 137474 times.
|
470022 | for (int dirIdx = 0; dirIdx < dim; ++dirIdx) |
| 408 | { | ||
| 409 | 313348 | const auto analyticalSolution | |
| 410 |
2/2✓ Branch 3 taken 4438 times.
✓ Branch 4 taken 270510 times.
|
940044 | = problem_->analyticalSolution(scv.dofPosition(), time)[dirIdx]; |
| 411 | 313348 | const auto numericalSolution | |
| 412 |
4/4✓ Branch 0 taken 4438 times.
✓ Branch 1 taken 270510 times.
✓ Branch 2 taken 4438 times.
✓ Branch 3 taken 270510 times.
|
626696 | = curSol[scv.dofIndex()][dirIdx]; |
| 413 | |||
| 414 | 626696 | const Scalar vError = absDiff_(analyticalSolution, numericalSolution); | |
| 415 | 626696 | const Scalar vReference = absDiff_(analyticalSolution, 0.0); | |
| 416 | |||
| 417 |
4/4✓ Branch 0 taken 4438 times.
✓ Branch 1 taken 270510 times.
✓ Branch 2 taken 1092 times.
✓ Branch 3 taken 273856 times.
|
318230 | maxError[dirIdx] = std::max(maxError[dirIdx], vError); |
| 418 |
4/4✓ Branch 0 taken 1092 times.
✓ Branch 1 taken 273856 times.
✓ Branch 2 taken 1092 times.
✓ Branch 3 taken 273856 times.
|
627970 | maxReference[dirIdx] = std::max(maxReference[dirIdx], vReference); |
| 419 | 626696 | sumError[dirIdx] += vError * vError * Extrusion::volume(fvGeometry, scv); | |
| 420 | 940044 | sumReference[dirIdx] += vReference * vReference * Extrusion::volume(fvGeometry, scv); | |
| 421 | } | ||
| 422 | } | ||
| 423 | else if constexpr (GridGeometry::discMethod == DiscretizationMethods::pq1bubble) | ||
| 424 | { | ||
| 425 |
2/2✓ Branch 0 taken 100000 times.
✓ Branch 1 taken 25800 times.
|
125800 | if(!scv.isOverlapping()) |
| 426 | { | ||
| 427 | 100000 | totalVolume_ += Extrusion::volume(fvGeometry, scv); | |
| 428 |
2/2✓ Branch 0 taken 200000 times.
✓ Branch 1 taken 100000 times.
|
300000 | for (int dirIdx = 0; dirIdx < dim; ++dirIdx) |
| 429 | { | ||
| 430 | 200000 | const auto analyticalSolution | |
| 431 |
2/2✓ Branch 3 taken 2371 times.
✓ Branch 4 taken 197629 times.
|
600000 | = problem_->analyticalSolution(scv.dofPosition(), time)[dirIdx]; |
| 432 | 200000 | const auto numericalSolution | |
| 433 |
4/4✓ Branch 0 taken 2371 times.
✓ Branch 1 taken 197629 times.
✓ Branch 2 taken 2371 times.
✓ Branch 3 taken 197629 times.
|
400000 | = curSol[scv.dofIndex()][dirIdx]; |
| 434 | |||
| 435 | 400000 | const Scalar vError = absDiff_(analyticalSolution, numericalSolution); | |
| 436 | 400000 | const Scalar vReference = absDiff_(analyticalSolution, 0.0); | |
| 437 | |||
| 438 |
4/4✓ Branch 0 taken 2371 times.
✓ Branch 1 taken 197629 times.
✓ Branch 2 taken 191 times.
✓ Branch 3 taken 199809 times.
|
202371 | maxError[dirIdx] = std::max(maxError[dirIdx], vError); |
| 439 |
4/4✓ Branch 0 taken 191 times.
✓ Branch 1 taken 199809 times.
✓ Branch 2 taken 191 times.
✓ Branch 3 taken 199809 times.
|
400191 | maxReference[dirIdx] = std::max(maxReference[dirIdx], vReference); |
| 440 | 400000 | sumError[dirIdx] += vError * vError * Extrusion::volume(fvGeometry, scv); | |
| 441 | 600000 | sumReference[dirIdx] += vReference * vReference * Extrusion::volume(fvGeometry, scv); | |
| 442 | } | ||
| 443 | } | ||
| 444 | } | ||
| 445 | else | ||
| 446 | DUNE_THROW(Dune::InvalidStateException, | ||
| 447 | "Unknown momentum discretization scheme in Navier-Stokes error calculation"); | ||
| 448 | } | ||
| 449 | // pressure errors | ||
| 450 | else | ||
| 451 | { | ||
| 452 |
2/2✓ Branch 0 taken 2271 times.
✓ Branch 1 taken 540597 times.
|
1643336 | totalVolume_ += Extrusion::volume(fvGeometry, scv); |
| 453 | // compute the pressure errors | ||
| 454 | using Indices = typename Problem::Indices; | ||
| 455 | 1643336 | const auto analyticalSolution | |
| 456 |
6/6✓ Branch 0 taken 2271 times.
✓ Branch 1 taken 540597 times.
✓ Branch 2 taken 2271 times.
✓ Branch 3 taken 540597 times.
✓ Branch 4 taken 2271 times.
✓ Branch 5 taken 540597 times.
|
4930008 | = problem_->analyticalSolution(scv.dofPosition(), time)[Indices::pressureIdx]; |
| 457 | 1643336 | const auto numericalSolution | |
| 458 |
4/6✓ Branch 0 taken 5625 times.
✓ Branch 1 taken 816043 times.
✓ Branch 2 taken 5607 times.
✓ Branch 3 taken 796861 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
3248272 | = curSol[scv.dofIndex()][Indices::pressureIdx]; |
| 459 | |||
| 460 | 3286672 | const Scalar pError = absDiff_(analyticalSolution, numericalSolution); | |
| 461 | 3286672 | const Scalar pReference = absDiff_(analyticalSolution, 0.0); | |
| 462 | |||
| 463 | 3286672 | maxError[0] = std::max(maxError[0], pError); | |
| 464 |
2/2✓ Branch 0 taken 6878 times.
✓ Branch 1 taken 814790 times.
|
1643336 | maxReference[0] = std::max(maxReference[0], pReference); |
| 465 | 1643336 | sumError[0] += pError * pError * Extrusion::volume(fvGeometry, scv); | |
| 466 | 3286672 | sumReference[0] += pReference * pReference * Extrusion::volume(fvGeometry, scv); | |
| 467 | } | ||
| 468 | } | ||
| 469 | } | ||
| 470 | |||
| 471 | // calculate errors | ||
| 472 |
2/2✓ Branch 0 taken 346 times.
✓ Branch 1 taken 225 times.
|
1073 | for (int i = 0; i < ErrorVector::size(); ++i) |
| 473 | { | ||
| 474 | 1084 | l2Absolute_[i] = std::sqrt(sumError[i] / totalVolume_); | |
| 475 | 1522 | l2Relative_[i] = std::sqrt(sumError[i] / sumReference[i]); | |
| 476 | 1084 | lInfAbsolute_[i] = maxError[i]; | |
| 477 | 1960 | lInfRelative_[i] = maxError[i] / maxReference[i]; | |
| 478 | } | ||
| 479 | 427 | } | |
| 480 | |||
| 481 | template<class T> | ||
| 482 | ✗ | T absDiff_(const T& a, const T& b) const | |
| 483 |
12/12✓ Branch 0 taken 12282 times.
✓ Branch 1 taken 1259390 times.
✓ Branch 2 taken 12282 times.
✓ Branch 3 taken 1259390 times.
✓ Branch 4 taken 12282 times.
✓ Branch 5 taken 1259390 times.
✓ Branch 6 taken 19289 times.
✓ Branch 7 taken 3228791 times.
✓ Branch 8 taken 19289 times.
✓ Branch 9 taken 3228791 times.
✓ Branch 10 taken 19289 times.
✓ Branch 11 taken 3228791 times.
|
13559256 | { using std::abs; return abs(a-b); } |
| 484 | |||
| 485 | std::shared_ptr<const Problem> problem_; | ||
| 486 | |||
| 487 | ErrorVector l2Absolute_; | ||
| 488 | ErrorVector l2Relative_; | ||
| 489 | ErrorVector lInfAbsolute_; | ||
| 490 | ErrorVector lInfRelative_; | ||
| 491 | Scalar time_; | ||
| 492 | Scalar hMax_; | ||
| 493 | Scalar totalVolume_; | ||
| 494 | std::size_t numDofs_; | ||
| 495 | }; | ||
| 496 | |||
| 497 | /*! | ||
| 498 | * \brief Compute errors between an analytical solution and the numerical approximation | ||
| 499 | */ | ||
| 500 | template<class MomentumProblem, class MassProblem, class Scalar = double> | ||
| 501 | class Errors | ||
| 502 | { | ||
| 503 | static constexpr int dim | ||
| 504 | = std::decay_t<decltype(std::declval<MomentumProblem>().gridGeometry())>::GridView::dimension; | ||
| 505 | |||
| 506 | using ErrorVector = Dune::FieldVector<Scalar, dim+1>; | ||
| 507 | using NumSubProblemVector = Dune::FieldVector<Scalar, 2>; | ||
| 508 | |||
| 509 | public: | ||
| 510 | template<class SolutionVector> | ||
| 511 | 32 | Errors(std::shared_ptr<const MomentumProblem> momentumProblem, | |
| 512 | std::shared_ptr<const MassProblem> massProblem, | ||
| 513 | const SolutionVector& curSol, | ||
| 514 | Scalar time = 0.0) | ||
| 515 | 64 | : momentumErrors_(momentumProblem, curSol[Dune::Indices::_0], time) | |
| 516 |
11/30✓ Branch 3 taken 32 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 32 times.
✗ Branch 6 not taken.
✓ Branch 10 taken 32 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 32 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 32 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 32 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 32 times.
✗ Branch 22 not taken.
✓ Branch 24 taken 32 times.
✗ Branch 25 not taken.
✓ Branch 27 taken 32 times.
✗ Branch 28 not taken.
✓ Branch 30 taken 32 times.
✗ Branch 31 not taken.
✓ Branch 33 taken 32 times.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
|
352 | , massErrors_(massProblem, curSol[Dune::Indices::_1], time) |
| 517 |
1/2✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
|
32 | { update(curSol, time); } |
| 518 | |||
| 519 | /*! | ||
| 520 | * \brief Computes errors between an analytical solution and the numerical approximation | ||
| 521 | * | ||
| 522 | * \param curSol The current solution vector | ||
| 523 | * \param time The current time | ||
| 524 | */ | ||
| 525 | template<class SolutionVector> | ||
| 526 | 69 | void update(const SolutionVector& curSol, Scalar time = 0.0) | |
| 527 | { | ||
| 528 | 138 | momentumErrors_.update(curSol[Dune::Indices::_0], time); | |
| 529 | 138 | massErrors_.update(curSol[Dune::Indices::_1], time); | |
| 530 | |||
| 531 | 69 | time_ = time; | |
| 532 | |||
| 533 | 138 | const auto& l2AbsoluteMomentum = momentumErrors_.l2Absolute(); | |
| 534 | 138 | const auto& l2RelativeMomentum = momentumErrors_.l2Relative(); | |
| 535 | 138 | const auto& lInfAbsoluteMomentum = momentumErrors_.lInfAbsolute(); | |
| 536 | 138 | const auto& lInfRelativeMomentum = momentumErrors_.lInfRelative(); | |
| 537 | |||
| 538 | 138 | const auto& l2AbsoluteMass = massErrors_.l2Absolute(); | |
| 539 | 138 | const auto& l2RelativeMass = massErrors_.l2Relative(); | |
| 540 | 138 | const auto& lInfAbsoluteMass = massErrors_.lInfAbsolute(); | |
| 541 | 138 | const auto& lInfRelativeMass = massErrors_.lInfRelative(); | |
| 542 | |||
| 543 | 69 | l2Absolute_[0] = l2AbsoluteMass[0]; | |
| 544 | 69 | l2Relative_[0] = l2RelativeMass[0]; | |
| 545 | 69 | lInfAbsolute_[0] = lInfAbsoluteMass[0]; | |
| 546 | 69 | lInfRelative_[0] = lInfRelativeMass[0]; | |
| 547 | |||
| 548 | 345 | std::copy( l2AbsoluteMomentum.begin(), l2AbsoluteMomentum.end(), l2Absolute_.begin() + 1 ); | |
| 549 | 345 | std::copy( l2RelativeMomentum.begin(), l2RelativeMomentum.end(), l2Relative_.begin() + 1 ); | |
| 550 | 345 | std::copy( lInfAbsoluteMomentum.begin(), lInfAbsoluteMomentum.end(), lInfAbsolute_.begin() + 1 ); | |
| 551 | 345 | std::copy( lInfRelativeMomentum.begin(), lInfRelativeMomentum.end(), lInfRelative_.begin() + 1 ); | |
| 552 | |||
| 553 | 138 | hMax_[0] = massErrors_.hMax(); | |
| 554 | 138 | totalVolume_[0] = massErrors_.totalVolume(); | |
| 555 | 207 | numDofs_[0] = massErrors_.numDofs(); | |
| 556 | 138 | hMax_[1] = momentumErrors_.hMax(); | |
| 557 | 138 | totalVolume_[1] = momentumErrors_.totalVolume(); | |
| 558 | 207 | numDofs_[1] = momentumErrors_.numDofs(); | |
| 559 | 69 | } | |
| 560 | |||
| 561 | //! The (absolute) discrete l2 error | ||
| 562 | 52 | const ErrorVector& l2Absolute() const { return l2Absolute_; } | |
| 563 | //! The relative discrete l2 error (relative to the discrete l2 norm of the reference solution) | ||
| 564 | 52 | const ErrorVector& l2Relative() const { return l2Relative_; } | |
| 565 | //! The (absolute) discrete l-infinity error | ||
| 566 | 52 | const ErrorVector& lInfAbsolute() const { return lInfAbsolute_; } | |
| 567 | //! The relative discrete l-infinity error (relative to the discrete loo norm of the reference solution) | ||
| 568 | 52 | const ErrorVector& lInfRelative() const { return lInfRelative_; } | |
| 569 | |||
| 570 | //! Volume of scvs considered in error calculation | ||
| 571 | const NumSubProblemVector& totalVolume() const { return totalVolume_; } | ||
| 572 | //! Number of scvs considered in error calculation | ||
| 573 | const NumSubProblemVector& numDofs() const { return numDofs_; } | ||
| 574 | //! Maximum diameter of primal grid elements | ||
| 575 | const NumSubProblemVector& hMax() const { return hMax_; } | ||
| 576 | |||
| 577 | //! Time corresponding to the error (returns 0 per default) | ||
| 578 | ✗ | Scalar time() const { return time_; } | |
| 579 | |||
| 580 | private: | ||
| 581 | ErrorsSubProblem<MomentumProblem, Scalar> momentumErrors_; | ||
| 582 | ErrorsSubProblem<MassProblem, Scalar> massErrors_; | ||
| 583 | |||
| 584 | ErrorVector l2Absolute_; | ||
| 585 | ErrorVector l2Relative_; | ||
| 586 | ErrorVector lInfAbsolute_; | ||
| 587 | ErrorVector lInfRelative_; | ||
| 588 | |||
| 589 | NumSubProblemVector hMax_; | ||
| 590 | NumSubProblemVector totalVolume_; | ||
| 591 | NumSubProblemVector numDofs_; | ||
| 592 | |||
| 593 | Scalar time_; | ||
| 594 | }; | ||
| 595 | |||
| 596 | template<class Problem, class SolutionVector> | ||
| 597 | ErrorsSubProblem(std::shared_ptr<Problem>, SolutionVector&&) | ||
| 598 | -> ErrorsSubProblem<Problem>; | ||
| 599 | |||
| 600 | template<class Problem, class SolutionVector, class Scalar> | ||
| 601 | ErrorsSubProblem(std::shared_ptr<Problem>, SolutionVector&&, Scalar) | ||
| 602 | -> ErrorsSubProblem<Problem, Scalar>; | ||
| 603 | |||
| 604 | template<class MomentumProblem, class MassProblem, class SolutionVector> | ||
| 605 | Errors(std::shared_ptr<MomentumProblem>, std::shared_ptr<MassProblem>, SolutionVector&&) | ||
| 606 | -> Errors<MomentumProblem, MassProblem>; | ||
| 607 | |||
| 608 | template<class MomentumProblem, class MassProblem, class SolutionVector, class Scalar> | ||
| 609 | Errors(std::shared_ptr<MomentumProblem>, std::shared_ptr<MassProblem>, SolutionVector&&, Scalar) | ||
| 610 | -> Errors<MomentumProblem, MassProblem, Scalar>; | ||
| 611 | |||
| 612 | |||
| 613 | /*! | ||
| 614 | * \brief An error CSV file writer | ||
| 615 | */ | ||
| 616 | template<class MomentumProblem, class MassProblem, class Scalar = double> | ||
| 617 |
3/10✓ Branch 0 taken 17 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 7 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
40 | class ErrorCSVWriter |
| 618 | { | ||
| 619 | static constexpr int dim | ||
| 620 | = std::decay_t<decltype(std::declval<MomentumProblem>().gridGeometry())>::GridView::dimension; | ||
| 621 | |||
| 622 | public: | ||
| 623 | 20 | ErrorCSVWriter(std::shared_ptr<const MomentumProblem> momentumProblem, | |
| 624 | std::shared_ptr<const MassProblem> massProblem, | ||
| 625 | const std::string& suffix = "") | ||
| 626 |
11/20✓ Branch 0 taken 13 times.
✓ Branch 1 taken 7 times.
✓ Branch 2 taken 13 times.
✓ Branch 3 taken 7 times.
✓ Branch 5 taken 13 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 13 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 7 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 20 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 20 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 20 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✓ Branch 23 taken 20 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
|
40 | : name_(massProblem->name() + (suffix.empty() ? "_error" : "_error_" + suffix)) |
| 627 | { | ||
| 628 |
3/6✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 20 times.
✗ Branch 8 not taken.
|
60 | const int numCCDofs = massProblem->gridGeometry().numDofs(); |
| 629 |
3/6✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 20 times.
✗ Branch 8 not taken.
|
60 | const int numFaceDofs = momentumProblem->gridGeometry().numDofs(); |
| 630 | |||
| 631 | // print auxiliary file with the number of dofs | ||
| 632 |
2/6✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
40 | std::ofstream logFileDofs(name_ + "_dofs.csv", std::ios::trunc); |
| 633 | 20 | logFileDofs << "cc dofs, face dofs, all dofs\n" | |
| 634 |
4/8✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 20 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 20 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 20 times.
✗ Branch 13 not taken.
|
60 | << numCCDofs << ", " << numFaceDofs << ", " << numCCDofs + numFaceDofs << "\n"; |
| 635 | |||
| 636 | // clear error file | ||
| 637 |
2/4✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
|
40 | std::ofstream logFile(name_ + ".csv", std::ios::trunc); |
| 638 | |||
| 639 | // write header | ||
| 640 |
1/2✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
|
20 | logFile << "time"; |
| 641 | using ErrorNames = std::vector<std::string>; | ||
| 642 |
4/8✓ Branch 0 taken 80 times.
✓ Branch 1 taken 20 times.
✓ Branch 3 taken 80 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 80 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
100 | for (const std::string e : { "L2Abs", "L2Rel", "LinfAbs", "LinfRel" }) |
| 643 |
16/36✓ Branch 1 taken 80 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 80 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 80 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 80 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 80 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 80 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 80 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 80 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 80 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 320 times.
✓ Branch 29 taken 80 times.
✗ Branch 30 not taken.
✓ Branch 31 taken 320 times.
✗ Branch 32 not taken.
✓ Branch 33 taken 80 times.
✗ Branch 34 not taken.
✓ Branch 35 taken 80 times.
✗ Branch 36 not taken.
✓ Branch 37 taken 80 times.
✗ Branch 38 not taken.
✓ Branch 39 taken 80 times.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
|
480 | printError_(logFile, ErrorNames({ e + "(p)", e + "(u)", e + "(v)", e + "(w)" }), "{:s}"); |
| 644 |
1/2✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
|
20 | logFile << "\n"; |
| 645 | 20 | } | |
| 646 | |||
| 647 | 47 | void printErrors(const Errors<MomentumProblem, MassProblem>& errors) const | |
| 648 | { | ||
| 649 | // append to error file | ||
| 650 |
1/2✓ Branch 2 taken 47 times.
✗ Branch 3 not taken.
|
94 | std::ofstream logFile(name_ + ".csv", std::ios::app); |
| 651 | |||
| 652 |
3/8✓ Branch 2 taken 47 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 47 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 47 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
47 | logFile << Fmt::format("{:.5e}", errors.time()); |
| 653 |
5/12✓ Branch 1 taken 47 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 47 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 47 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 47 times.
✓ Branch 12 taken 47 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
94 | printError_(logFile, errors.l2Absolute()); |
| 654 |
5/12✓ Branch 1 taken 47 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 47 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 47 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 47 times.
✓ Branch 12 taken 47 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
94 | printError_(logFile, errors.l2Relative()); |
| 655 |
5/12✓ Branch 1 taken 47 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 47 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 47 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 47 times.
✓ Branch 12 taken 47 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
94 | printError_(logFile, errors.lInfAbsolute()); |
| 656 |
5/12✓ Branch 1 taken 47 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 47 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 47 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 47 times.
✓ Branch 12 taken 47 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
94 | printError_(logFile, errors.lInfRelative()); |
| 657 | |||
| 658 |
1/2✓ Branch 1 taken 47 times.
✗ Branch 2 not taken.
|
47 | logFile << "\n"; |
| 659 | 47 | } | |
| 660 | |||
| 661 | private: | ||
| 662 | template<class Error> | ||
| 663 | 536 | void printError_(std::ofstream& logFile, const Error& error, const std::string& format = "{:.5e}") const | |
| 664 | { | ||
| 665 | using MassIndices = typename MassProblem::Indices; | ||
| 666 | using MomIndices = typename MomentumProblem::Indices; | ||
| 667 |
10/23✓ Branch 1 taken 268 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 268 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 268 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 268 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 268 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 268 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✓ Branch 19 taken 268 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 80 times.
✓ Branch 22 taken 188 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 188 times.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
|
1984 | logFile << Fmt::vformat(", " + format, Fmt::make_format_args(error[MassIndices::pressureIdx])); |
| 668 |
2/2✓ Branch 0 taken 528 times.
✓ Branch 1 taken 268 times.
|
1592 | for (int dirIdx = 0; dirIdx < dim; ++dirIdx) |
| 669 |
9/22✓ Branch 1 taken 528 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 528 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 528 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 528 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 528 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 528 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 528 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 528 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 528 times.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
|
4224 | logFile << Fmt::vformat(", " + format, Fmt::make_format_args(error[MomIndices::velocity(dirIdx)+1])); |
| 670 | 536 | } | |
| 671 | |||
| 672 | std::string name_; | ||
| 673 | }; | ||
| 674 | |||
| 675 | template<class MomentumProblem, class MassProblem> | ||
| 676 | ErrorCSVWriter(std::shared_ptr<MomentumProblem>, std::shared_ptr<MassProblem>) | ||
| 677 | -> ErrorCSVWriter<MomentumProblem, MassProblem>; | ||
| 678 | |||
| 679 | template<class MomentumProblem, class MassProblem> | ||
| 680 | ErrorCSVWriter(std::shared_ptr<MomentumProblem>, std::shared_ptr<MassProblem>, const std::string&) | ||
| 681 | -> ErrorCSVWriter<MomentumProblem, MassProblem>; | ||
| 682 | |||
| 683 | |||
| 684 | /*! | ||
| 685 | * \brief Append errors to the log file during a convergence test | ||
| 686 | */ | ||
| 687 | template<class MomentumProblem, class MassProblem> | ||
| 688 | 5 | void convergenceTestAppendErrors(std::ofstream& logFile, | |
| 689 | std::shared_ptr<MomentumProblem> momentumProblem, | ||
| 690 | std::shared_ptr<MassProblem> massProblem, | ||
| 691 | const Errors<MomentumProblem, MassProblem>& errors) | ||
| 692 | { | ||
| 693 | 15 | const auto numCCDofs = massProblem->gridGeometry().numDofs(); | |
| 694 | 15 | const auto numFaceDofs = momentumProblem->gridGeometry().numDofs(); | |
| 695 | |||
| 696 |
2/6✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 5 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
5 | logFile << Fmt::format( |
| 697 | "[ConvergenceTest] numCCDofs = {} numFaceDofs = {}", | ||
| 698 | numCCDofs, numFaceDofs | ||
| 699 | ); | ||
| 700 | |||
| 701 | 45 | const auto print = [&](const auto& e, const std::string& name){ | |
| 702 | using MassIndices = typename MassProblem::Indices; | ||
| 703 | using MomIndices = typename MomentumProblem::Indices; | ||
| 704 |
2/6✓ Branch 6 taken 20 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 20 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
100 | logFile << Fmt::format( |
| 705 | " {0}(p) = {1} {0}(vx) = {2} {0}(vy) = {3}", | ||
| 706 | name, e[MassIndices::pressureIdx], e[MomIndices::velocityXIdx], e[MomIndices::velocityYIdx] | ||
| 707 | ); | ||
| 708 | }; | ||
| 709 | |||
| 710 |
5/12✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 5 times.
✓ Branch 12 taken 5 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
10 | print(errors.l2Absolute(), "L2Abs"); |
| 711 |
5/12✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 5 times.
✓ Branch 12 taken 5 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
10 | print(errors.l2Relative(), "L2Rel"); |
| 712 |
5/12✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 5 times.
✓ Branch 12 taken 5 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
10 | print(errors.lInfAbsolute(), "LinfAbs"); |
| 713 |
4/10✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 5 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
10 | print(errors.lInfRelative(), "LinfRel"); |
| 714 | |||
| 715 | 5 | logFile << "\n"; | |
| 716 | 5 | } | |
| 717 | |||
| 718 | |||
| 719 | /*! | ||
| 720 | * \brief Append errors to the log file during a convergence test | ||
| 721 | */ | ||
| 722 | template<class MomentumProblem, class MassProblem> | ||
| 723 | 5 | void convergenceTestAppendErrors(std::shared_ptr<MomentumProblem> momentumProblem, | |
| 724 | std::shared_ptr<MassProblem> massProblem, | ||
| 725 | const Errors<MomentumProblem, MassProblem>& errors) | ||
| 726 | { | ||
| 727 |
0/2✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
15 | const auto logFileName = massProblem->name() + ".log"; |
| 728 |
2/4✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
|
10 | std::ofstream logFile(logFileName, std::ios::app); |
| 729 |
3/10✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 5 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
15 | convergenceTestAppendErrors(logFile, momentumProblem, massProblem, errors); |
| 730 | 5 | } | |
| 731 | |||
| 732 | /*! | ||
| 733 | * \brief Append errors to the log file during a convergence test | ||
| 734 | */ | ||
| 735 | template<class MomentumProblem> | ||
| 736 | 6 | void convergenceTestAppendErrorsMomentum(std::ofstream& logFile, | |
| 737 | std::shared_ptr<MomentumProblem> problem, | ||
| 738 | const ErrorsSubProblem<MomentumProblem>& errors) | ||
| 739 | { | ||
| 740 | 18 | const auto numFaceDofs = problem->gridGeometry().numDofs(); | |
| 741 | |||
| 742 |
2/6✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 6 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
6 | logFile << Fmt::format( |
| 743 | "[ConvergenceTest] numFaceDofs = {}", | ||
| 744 | numFaceDofs | ||
| 745 | ); | ||
| 746 | |||
| 747 | 54 | const auto print = [&](const auto& e, const std::string& name){ | |
| 748 | using MomIndices = typename MomentumProblem::Indices; | ||
| 749 |
2/6✓ Branch 5 taken 24 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 24 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
96 | logFile << Fmt::format( |
| 750 | " {0}(vx) = {1} {0}(vy) = {2}", | ||
| 751 | name, e[MomIndices::velocityXIdx], e[MomIndices::velocityYIdx] | ||
| 752 | ); | ||
| 753 | }; | ||
| 754 | |||
| 755 |
5/12✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 6 times.
✓ Branch 12 taken 6 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
12 | print(errors.l2Absolute(), "L2Abs"); |
| 756 |
5/12✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 6 times.
✓ Branch 12 taken 6 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
12 | print(errors.l2Relative(), "L2Rel"); |
| 757 |
5/12✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 6 times.
✓ Branch 12 taken 6 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
12 | print(errors.lInfAbsolute(), "LinfAbs"); |
| 758 |
4/10✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 6 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
12 | print(errors.lInfRelative(), "LinfRel"); |
| 759 | |||
| 760 | 6 | logFile << "\n"; | |
| 761 | 6 | } | |
| 762 | |||
| 763 | |||
| 764 | /*! | ||
| 765 | * \brief Append errors to the log file during a convergence test | ||
| 766 | */ | ||
| 767 | template<class MomentumProblem> | ||
| 768 | 6 | void convergenceTestAppendErrorsMomentum(std::shared_ptr<MomentumProblem> problem, | |
| 769 | const ErrorsSubProblem<MomentumProblem>& errors) | ||
| 770 | { | ||
| 771 |
0/2✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
18 | const auto logFileName = problem->name() + ".log"; |
| 772 |
2/4✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
|
12 | std::ofstream logFile(logFileName, std::ios::app); |
| 773 |
2/6✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
12 | convergenceTestAppendErrorsMomentum(logFile, problem, errors); |
| 774 | 6 | } | |
| 775 | |||
| 776 | } // end namespace NavierStokesTest | ||
| 777 | |||
| 778 | } // end namespace Dumux | ||
| 779 | |||
| 780 | #endif | ||
| 781 |