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