| 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 | * \brief Test for the instationary staggered grid Navier-Stokes model | ||
| 11 | * with analytical solution (Angeli et al. 2017, \cite Angeli2017). | ||
| 12 | */ | ||
| 13 | |||
| 14 | #ifndef DUMUX_ANGELI_TEST_PROBLEM_HH | ||
| 15 | #define DUMUX_ANGELI_TEST_PROBLEM_HH | ||
| 16 | |||
| 17 | #include <dune/geometry/quadraturerules.hh> | ||
| 18 | |||
| 19 | #include <dumux/common/parameters.hh> | ||
| 20 | #include <dumux/common/properties.hh> | ||
| 21 | |||
| 22 | #include <dumux/freeflow/navierstokes/boundarytypes.hh> | ||
| 23 | |||
| 24 | namespace Dumux { | ||
| 25 | |||
| 26 | /*! | ||
| 27 | * \ingroup NavierStokesTests | ||
| 28 | * \brief Test problem for the staggered grid (Angeli et al. 2017, \cite Angeli2017). | ||
| 29 | * | ||
| 30 | * The unsteady, 2D, incompressible Navier-Stokes equations for a zero source and a Newtonian | ||
| 31 | * flow is solved and compared to an analytical solution (sums/products of trigonometric functions). | ||
| 32 | * The velocities and pressures decay exponentially. The Dirichlet boundary conditions are | ||
| 33 | * time-dependent and consistent with the analytical solution. | ||
| 34 | */ | ||
| 35 | template <class TypeTag, class BaseProblem> | ||
| 36 | 4 | class AngeliTestProblem : public BaseProblem | |
| 37 | { | ||
| 38 | using ParentType = BaseProblem; | ||
| 39 | |||
| 40 | using BoundaryTypes = typename ParentType::BoundaryTypes; | ||
| 41 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 42 | using GridView = typename GridGeometry::GridView; | ||
| 43 | using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>; | ||
| 44 | using InitialValues = typename ParentType::InitialValues; | ||
| 45 | using Sources = typename ParentType::Sources; | ||
| 46 | using DirichletValues = typename ParentType::DirichletValues; | ||
| 47 | using BoundaryFluxes = typename ParentType::BoundaryFluxes; | ||
| 48 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 49 | using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>; | ||
| 50 | using SubControlVolume = typename GridGeometry::SubControlVolume; | ||
| 51 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 52 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 53 | |||
| 54 | static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld; | ||
| 55 | using Element = typename FVElementGeometry::Element; | ||
| 56 | |||
| 57 | using GlobalPosition = typename Element::Geometry::GlobalCoordinate; | ||
| 58 | using VelocityVector = Dune::FieldVector<Scalar, dimWorld>; | ||
| 59 | |||
| 60 | using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>; | ||
| 61 | |||
| 62 | |||
| 63 | public: | ||
| 64 | using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices; | ||
| 65 | |||
| 66 | 8 | AngeliTestProblem(std::shared_ptr<const GridGeometry> gridGeometry, std::shared_ptr<CouplingManager> couplingManager) | |
| 67 |
7/20✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 4 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 4 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 4 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 4 times.
✓ Branch 18 taken 4 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
|
32 | : ParentType(gridGeometry, couplingManager) |
| 68 | { | ||
| 69 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
8 | kinematicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0); |
| 70 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
8 | useNeumann_ = getParam<bool>("Problem.UseNeumann", false); |
| 71 |
1/2✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
|
8 | interpolateExactVelocity_ = getParam<bool>("Problem.InterpolateExactVelocity", false); |
| 72 | 8 | } | |
| 73 | |||
| 74 | /*! | ||
| 75 | * \name Boundary conditions | ||
| 76 | */ | ||
| 77 | // \{ | ||
| 78 | |||
| 79 | /*! | ||
| 80 | * \brief Specifies which kind of boundary condition should be | ||
| 81 | * used for which equation on a given boundary control volume. | ||
| 82 | * | ||
| 83 | * \param globalPos The position of the center of the finite volume | ||
| 84 | */ | ||
| 85 | 55600 | BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const | |
| 86 | { | ||
| 87 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 54000 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 13600 times.
|
123200 | BoundaryTypes values; |
| 88 | |||
| 89 | if constexpr (ParentType::isMomentumProblem()) | ||
| 90 | { | ||
| 91 | static constexpr Scalar eps = 1e-8; | ||
| 92 |
1/12✗ Branch 0 not taken.
✓ Branch 1 taken 27800 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
55600 | if (useNeumann_ && (globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps || |
| 93 | ✗ | globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps)) | |
| 94 | values.setAllNeumann(); | ||
| 95 | else | ||
| 96 | values.setAllDirichlet(); | ||
| 97 | } | ||
| 98 | else | ||
| 99 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 54000 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 13600 times.
|
67600 | values.setAllNeumann(); |
| 100 | |||
| 101 | 55600 | return values; | |
| 102 | |||
| 103 | } | ||
| 104 | |||
| 105 | //! Enable internal Dirichlet constraints | ||
| 106 | static constexpr bool enableInternalDirichletConstraints() | ||
| 107 | { return !ParentType::isMomentumProblem(); } | ||
| 108 | |||
| 109 | |||
| 110 | /*! | ||
| 111 | * \brief Tag a degree of freedom to carry internal Dirichlet constraints. | ||
| 112 | * If true is returned for a dof, the equation for this dof is replaced | ||
| 113 | * by the constraint that its primary variable values must match the | ||
| 114 | * user-defined values obtained from the function internalDirichlet(), | ||
| 115 | * which must be defined in the problem. | ||
| 116 | * | ||
| 117 | * \param element The finite element | ||
| 118 | * \param scv The sub-control volume | ||
| 119 | */ | ||
| 120 | 602100 | std::bitset<DirichletValues::dimension> hasInternalDirichletConstraint(const Element& element, const SubControlVolume& scv) const | |
| 121 | { | ||
| 122 | 602100 | std::bitset<DirichletValues::dimension> values; | |
| 123 | |||
| 124 | // We don't need internal Dirichlet conditions if a Neumann BC is set for the momentum balance (which accounts for the pressure). | ||
| 125 | // If only Dirichlet BCs are set for the momentum balance, fix the pressure at some cells such that the solution is fully defined. | ||
| 126 |
1/2✓ Branch 0 taken 602100 times.
✗ Branch 1 not taken.
|
602100 | if (!useNeumann_) |
| 127 | { | ||
| 128 | 1806300 | const auto fvGeometry = localView(this->gridGeometry()).bindElement(element); | |
| 129 |
4/4✓ Branch 0 taken 42228 times.
✓ Branch 1 taken 559872 times.
✓ Branch 2 taken 42228 times.
✓ Branch 3 taken 559872 times.
|
1204200 | if (fvGeometry.hasBoundaryScvf()) |
| 130 | 42228 | values.set(Indices::pressureIdx); | |
| 131 | } | ||
| 132 | |||
| 133 | 602100 | return values; | |
| 134 | } | ||
| 135 | |||
| 136 | /*! | ||
| 137 | * \brief Define the values of internal Dirichlet constraints for a degree of freedom. | ||
| 138 | * \param element The finite element | ||
| 139 | * \param scv The sub-control volume | ||
| 140 | */ | ||
| 141 | ✗ | DirichletValues internalDirichlet(const Element& element, const SubControlVolume& scv) const | |
| 142 | 135000 | { return analyticalSolution(scv.center(), time_); } | |
| 143 | |||
| 144 | |||
| 145 | |||
| 146 | /*! | ||
| 147 | * \brief Evaluate the boundary conditions for a dirichlet | ||
| 148 | * control volume face (velocities) | ||
| 149 | * | ||
| 150 | * \param element The finite element | ||
| 151 | * \param scvf the sub control volume face | ||
| 152 | * | ||
| 153 | * Concerning the usage of averagedVelocity_, see the explanation of the initial function. | ||
| 154 | */ | ||
| 155 | 245680 | DirichletValues dirichlet(const Element& element, const SubControlVolumeFace& scvf) const | |
| 156 | { | ||
| 157 |
2/2✓ Branch 0 taken 9440 times.
✓ Branch 1 taken 113400 times.
|
245680 | if (ParentType::isMomentumProblem() && interpolateExactVelocity_) |
| 158 | { | ||
| 159 | 56640 | const auto fvGeometry = localView(this->gridGeometry()).bindElement(element); | |
| 160 | 18880 | return velocityDirichlet_(fvGeometry.geometry(scvf)); | |
| 161 | } | ||
| 162 | else | ||
| 163 | 453600 | return analyticalSolution(scvf.center(), time_); | |
| 164 | } | ||
| 165 | |||
| 166 | /*! | ||
| 167 | * \brief Evaluates the boundary conditions for a Neumann control volume. | ||
| 168 | * | ||
| 169 | * \param element The element for which the Neumann boundary condition is set | ||
| 170 | * \param fvGeometry The fvGeometry | ||
| 171 | * \param elemVolVars The element volume variables | ||
| 172 | * \param elemFaceVars The element face variables | ||
| 173 | * \param scvf The boundary sub control volume face | ||
| 174 | */ | ||
| 175 | template<class ElementVolumeVariables, class ElementFluxVariablesCache> | ||
| 176 | 108000 | BoundaryFluxes neumann(const Element& element, | |
| 177 | const FVElementGeometry& fvGeometry, | ||
| 178 | const ElementVolumeVariables& elemVolVars, | ||
| 179 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 180 | const SubControlVolumeFace& scvf) const | ||
| 181 | { | ||
| 182 | 108000 | BoundaryFluxes values(0.0); | |
| 183 | |||
| 184 | if constexpr (!ParentType::isMomentumProblem()) | ||
| 185 | { | ||
| 186 | 216000 | const auto insideDensity = elemVolVars[scvf.insideScvIdx()].density(); | |
| 187 | 324000 | values[Indices::conti0EqIdx] = this->faceVelocity(element, fvGeometry, scvf) * insideDensity * scvf.unitOuterNormal(); | |
| 188 | } | ||
| 189 | else | ||
| 190 | { | ||
| 191 | using std::sin; | ||
| 192 | using std::cos; | ||
| 193 | ✗ | Dune::FieldMatrix<Scalar, 2, 2> momentumFlux(0.0); | |
| 194 | ✗ | const auto x = scvf.ipGlobal()[0]; | |
| 195 | ✗ | const auto y = scvf.ipGlobal()[1]; | |
| 196 | ✗ | const Scalar mu = kinematicViscosity_; | |
| 197 | ✗ | const Scalar t = time_; | |
| 198 | ✗ | momentumFlux[0][0] = M_PI*M_PI *(-4.0*mu*exp(10.0*M_PI*M_PI*mu*t)*sin(M_PI*x)*sin(2*M_PI*y) + (4.0*sin(2*M_PI*y)*sin(2*M_PI*y)*cos(M_PI*x)*cos(M_PI*x) - 1.0*cos(2*M_PI*x) - 0.25*cos(4*M_PI*y))*exp(5.0*M_PI*M_PI*mu*t))*exp(-15.0*M_PI*M_PI*mu*t); | |
| 199 | ✗ | momentumFlux[0][1] = M_PI*M_PI *( 3.0*mu*exp(10.0*M_PI*M_PI*mu*t) - 2.0*exp(5.0*M_PI*M_PI*mu*t)*sin(M_PI*x)*sin(2*M_PI*y))*exp(-15.0*M_PI*M_PI*mu*t)*cos(M_PI*x)*cos(2*M_PI*y); | |
| 200 | ✗ | momentumFlux[1][0] = momentumFlux[0][1]; | |
| 201 | ✗ | momentumFlux[1][1] = M_PI*M_PI *( 4.0*mu*exp(10.0*M_PI*M_PI*mu*t)*sin(M_PI*x)*sin(2*M_PI*y) + (sin(M_PI*x)*sin(M_PI*x)*cos(2*M_PI*y)*cos(2*M_PI*y) - 1.0*cos(2*M_PI*x) - 0.25*cos(4*M_PI*y))*exp(5.0*M_PI*M_PI*mu*t))*exp(-15.0*M_PI*M_PI*mu*t); | |
| 202 | |||
| 203 | ✗ | const auto normal = scvf.unitOuterNormal(); | |
| 204 | ✗ | momentumFlux.mv(normal, values); | |
| 205 | } | ||
| 206 | |||
| 207 | 108000 | return values; | |
| 208 | } | ||
| 209 | |||
| 210 | |||
| 211 | |||
| 212 | /*! | ||
| 213 | * \brief Returns the analytical solution of the problem at a given time and position. | ||
| 214 | * \param globalPos The global position | ||
| 215 | * \param time The current simulation time | ||
| 216 | */ | ||
| 217 | 1310160 | DirichletValues analyticalSolution(const GlobalPosition& globalPos, Scalar time) const | |
| 218 | { | ||
| 219 | 2620320 | const Scalar x = globalPos[0]; | |
| 220 | 2620320 | const Scalar y = globalPos[1]; | |
| 221 | |||
| 222 | 1310160 | DirichletValues values; | |
| 223 | using std::exp; | ||
| 224 | using std::sin; | ||
| 225 | using std::cos; | ||
| 226 | |||
| 227 | if constexpr (ParentType::isMomentumProblem()) | ||
| 228 | { | ||
| 229 | 2070320 | values[Indices::velocityXIdx] = - 2.0 * M_PI * exp(- 5.0 * kinematicViscosity_ * M_PI * M_PI * time) * cos(M_PI * x) * sin(2.0 * M_PI * y); | |
| 230 | 2070320 | values[Indices::velocityYIdx] = M_PI * exp(- 5.0 * kinematicViscosity_ * M_PI * M_PI * time) * sin(M_PI * x) * cos(2.0 * M_PI * y); | |
| 231 | } | ||
| 232 | else | ||
| 233 | 275000 | values[Indices::pressureIdx] = - 0.25 * exp(-10.0 * kinematicViscosity_ * M_PI * M_PI * time) * M_PI * M_PI * (4.0 * cos(2.0 * M_PI * x) + cos(4.0 * M_PI * y)); | |
| 234 | |||
| 235 | 1310160 | return values; | |
| 236 | } | ||
| 237 | |||
| 238 | // \} | ||
| 239 | |||
| 240 | /*! | ||
| 241 | * \name Volume terms | ||
| 242 | */ | ||
| 243 | // \{ | ||
| 244 | |||
| 245 | /*! | ||
| 246 | * \brief Evaluates the initial value for a control volume (velocity) | ||
| 247 | */ | ||
| 248 | 10200 | InitialValues initial(const SubControlVolume& scv) const | |
| 249 | { | ||
| 250 | // We can either evaluate the analytical solution point-wise at the | ||
| 251 | // momentum balance integration points or use an averaged velocity. | ||
| 252 | |||
| 253 | // Not using the quadrature to integrate and average the velocity | ||
| 254 | // will yield a spurious pressure solution in the first time step | ||
| 255 | // because the discrete mass balance equation is not fulfilled when just taking | ||
| 256 | // point-wise values of the analytical solution. | ||
| 257 | |||
| 258 |
2/2✓ Branch 0 taken 5100 times.
✓ Branch 1 taken 5100 times.
|
10200 | if (!interpolateExactVelocity_) |
| 259 | 10200 | return analyticalSolution(scv.dofPosition(), 0.0); | |
| 260 | |||
| 261 | // Get the element intersection/facet corresponding corresponding to the dual grid scv | ||
| 262 | // (where the velocity DOFs are located) and use a quadrature to get the average velocity | ||
| 263 | // on that facet. | ||
| 264 | 10200 | const auto& element = this->gridGeometry().element(scv.elementIndex()); | |
| 265 | |||
| 266 |
1/2✗ Branch 3 not taken.
✓ Branch 4 taken 15200 times.
|
25400 | for (const auto& intersection : intersections(this->gridGeometry().gridView(), element)) |
| 267 | { | ||
| 268 |
4/4✓ Branch 0 taken 5100 times.
✓ Branch 1 taken 10100 times.
✓ Branch 2 taken 5100 times.
✓ Branch 3 taken 10100 times.
|
30400 | if (intersection.indexInInside() == scv.indexInElement()) |
| 269 | 5100 | return velocityDirichlet_(intersection.geometry()); | |
| 270 | |||
| 271 | } | ||
| 272 | ✗ | DUNE_THROW(Dune::InvalidStateException, "No intersection found"); | |
| 273 | } | ||
| 274 | |||
| 275 | |||
| 276 | /*! | ||
| 277 | * \brief Evaluates the initial value for an element (pressure) | ||
| 278 | */ | ||
| 279 | 5000 | InitialValues initial(const Element& element) const | |
| 280 | { | ||
| 281 | 10000 | return analyticalSolution(element.geometry().center(), 0.0); | |
| 282 | } | ||
| 283 | |||
| 284 | // \} | ||
| 285 | |||
| 286 | /*! | ||
| 287 | * \brief Updates the time information | ||
| 288 | */ | ||
| 289 | ✗ | void updateTime(Scalar t) | |
| 290 | { | ||
| 291 |
2/4✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
|
20 | time_ = t; |
| 292 | ✗ | } | |
| 293 | |||
| 294 | private: | ||
| 295 | template<class Geometry> | ||
| 296 | 14540 | DirichletValues velocityDirichlet_(const Geometry& geo) const | |
| 297 | { | ||
| 298 | 14540 | DirichletValues priVars(0.0); | |
| 299 | 19640 | const auto& quad = Dune::QuadratureRules<Scalar, Geometry::mydimension>::rule(geo.type(), 3); | |
| 300 | 101780 | for (auto&& qp : quad) | |
| 301 | { | ||
| 302 | 58160 | const auto w = qp.weight()*geo.integrationElement(qp.position()); | |
| 303 | 47960 | const auto globalPos = geo.global(qp.position()); | |
| 304 | 29080 | const auto sol = analyticalSolution(globalPos, time_); | |
| 305 | 87240 | priVars += sol*w; | |
| 306 | } | ||
| 307 | 14540 | priVars /= geo.volume(); | |
| 308 | 14540 | return priVars; | |
| 309 | } | ||
| 310 | |||
| 311 | Scalar kinematicViscosity_; | ||
| 312 | Scalar time_ = 0.0; | ||
| 313 | bool useNeumann_; | ||
| 314 | bool interpolateExactVelocity_; | ||
| 315 | }; | ||
| 316 | } // end namespace Dumux | ||
| 317 | |||
| 318 | #endif | ||
| 319 |