| 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 NavierStokesModel | ||
| 10 | * \copydoc Dumux::NavierStokesMomentumProblem | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_NAVIERSTOKES_MOMENTUM_PROBLEM_HH | ||
| 13 | #define DUMUX_NAVIERSTOKES_MOMENTUM_PROBLEM_HH | ||
| 14 | |||
| 15 | #include <dune/common/exceptions.hh> | ||
| 16 | #include <dune/common/fvector.hh> | ||
| 17 | #include <dune/common/typetraits.hh> | ||
| 18 | #include <dumux/common/numeqvector.hh> | ||
| 19 | #include <dumux/common/properties.hh> | ||
| 20 | #include <dumux/common/fvproblemwithspatialparams.hh> | ||
| 21 | #include <dumux/discretization/method.hh> | ||
| 22 | #include <dumux/freeflow/navierstokes/momentum/boundarytypes.hh> | ||
| 23 | |||
| 24 | namespace Dumux { | ||
| 25 | |||
| 26 | template<class TypeTag, class DiscretizationMethod> | ||
| 27 | class NavierStokesMomentumProblemImpl; | ||
| 28 | |||
| 29 | template<class TypeTag> | ||
| 30 | class NavierStokesMomentumProblemImpl<TypeTag, DiscretizationMethods::FCStaggered> | ||
| 31 | : public FVProblemWithSpatialParams<TypeTag> | ||
| 32 | { | ||
| 33 | using ParentType = FVProblemWithSpatialParams<TypeTag>; | ||
| 34 | using Implementation = GetPropType<TypeTag, Properties::Problem>; | ||
| 35 | |||
| 36 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 37 | using GridView = typename GridGeometry::GridView; | ||
| 38 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 39 | |||
| 40 | using GridVariables = GetPropType<TypeTag, Properties::GridVariables>; | ||
| 41 | using GridVolumeVariables = typename GridVariables::GridVolumeVariables; | ||
| 42 | using ElementVolumeVariables = typename GridVolumeVariables::LocalView; | ||
| 43 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 44 | |||
| 45 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 46 | using SubControlVolume = typename FVElementGeometry::SubControlVolume; | ||
| 47 | using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; | ||
| 48 | using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices; | ||
| 49 | |||
| 50 | enum { | ||
| 51 | dim = GridView::dimension, | ||
| 52 | dimWorld = GridView::dimensionworld | ||
| 53 | }; | ||
| 54 | |||
| 55 | using GlobalPosition = typename SubControlVolumeFace::GlobalPosition; | ||
| 56 | using VelocityVector = Dune::FieldVector<Scalar, dimWorld>; | ||
| 57 | using GravityVector = Dune::FieldVector<Scalar, dimWorld>; | ||
| 58 | using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>; | ||
| 59 | using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>; | ||
| 60 | |||
| 61 | static constexpr bool isCoupled_ = !std::is_empty_v<CouplingManager>; | ||
| 62 | |||
| 63 | |||
| 64 | public: | ||
| 65 | //! These types are used in place of the typical NumEqVector type. | ||
| 66 | //! In the numeqvector assembly type, only one equation per DOF (face) is considered | ||
| 67 | //! while the type here provides one entry for each world dimension. | ||
| 68 | using InitialValues = Dune::FieldVector<Scalar, dimWorld>; | ||
| 69 | using Sources = Dune::FieldVector<Scalar, dimWorld>; | ||
| 70 | using DirichletValues = Dune::FieldVector<Scalar, dimWorld>; | ||
| 71 | using BoundaryFluxes = Dune::FieldVector<Scalar, dimWorld>; | ||
| 72 | |||
| 73 | using MomentumFluxType = Dune::FieldVector<Scalar, dimWorld>; | ||
| 74 | |||
| 75 | //! Export the boundary types. | ||
| 76 | using BoundaryTypes = NavierStokesMomentumBoundaryTypes<ModelTraits::dim()>; | ||
| 77 | |||
| 78 | //! This problem is used for the momentum balance model. | ||
| 79 | static constexpr bool isMomentumProblem() { return true; } | ||
| 80 | |||
| 81 | /*! | ||
| 82 | * \brief The constructor | ||
| 83 | * \param gridGeometry The finite volume grid geometry | ||
| 84 | * \param couplingManager The coupling manager (couples mass and momentum equations) | ||
| 85 | * \param paramGroup The parameter group in which to look for runtime parameters first (default is "") | ||
| 86 | */ | ||
| 87 | 71 | NavierStokesMomentumProblemImpl(std::shared_ptr<const GridGeometry> gridGeometry, | |
| 88 | std::shared_ptr<CouplingManager> couplingManager, | ||
| 89 | const std::string& paramGroup = "") | ||
| 90 | : ParentType(gridGeometry, paramGroup) | ||
| 91 | 141 | , gravity_(0.0) | |
| 92 |
1/2✓ Branch 2 taken 71 times.
✗ Branch 3 not taken.
|
212 | , couplingManager_(couplingManager) |
| 93 | { | ||
| 94 |
3/4✓ Branch 1 taken 71 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 5 times.
✓ Branch 4 taken 66 times.
|
71 | if (getParamFromGroup<bool>(paramGroup, "Problem.EnableGravity")) |
| 95 | 5 | gravity_[dim-1] = -9.81; | |
| 96 | |||
| 97 |
1/2✓ Branch 1 taken 71 times.
✗ Branch 2 not taken.
|
71 | enableInertiaTerms_ = getParamFromGroup<bool>(paramGroup, "Problem.EnableInertiaTerms"); |
| 98 | 71 | } | |
| 99 | |||
| 100 | /*! | ||
| 101 | * \brief The constructor for usage without a coupling manager | ||
| 102 | * \param gridGeometry The finite volume grid geometry | ||
| 103 | * \param paramGroup The parameter group in which to look for runtime parameters first (default is "") | ||
| 104 | */ | ||
| 105 | 5 | NavierStokesMomentumProblemImpl(std::shared_ptr<const GridGeometry> gridGeometry, | |
| 106 | const std::string& paramGroup = "") | ||
| 107 |
1/2✓ Branch 2 taken 5 times.
✗ Branch 3 not taken.
|
10 | : NavierStokesMomentumProblemImpl(gridGeometry, {}, paramGroup) |
| 108 | 5 | {} | |
| 109 | |||
| 110 | /*! | ||
| 111 | * \brief Evaluate the source term for all phases within a given | ||
| 112 | * sub-control-volume. | ||
| 113 | * | ||
| 114 | * This is the method for the case where the source term is | ||
| 115 | * potentially solution dependent and requires some quantities that | ||
| 116 | * are specific to the fully-implicit method. | ||
| 117 | * | ||
| 118 | * \param element The finite element | ||
| 119 | * \param fvGeometry The finite-volume geometry | ||
| 120 | * \param elemVolVars All volume variables for the element | ||
| 121 | * \param scv The sub control volume | ||
| 122 | * | ||
| 123 | * For this method, the return parameter stores the conserved quantity rate | ||
| 124 | * generated or annihilate per volume unit. Positive values mean | ||
| 125 | * that the conserved quantity is created, negative ones mean that it vanishes. | ||
| 126 | * E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$. | ||
| 127 | */ | ||
| 128 | 96362830 | Sources source(const Element& element, | |
| 129 | const FVElementGeometry& fvGeometry, | ||
| 130 | const ElementVolumeVariables& elemVolVars, | ||
| 131 | const SubControlVolume& scv) const | ||
| 132 | { | ||
| 133 | // forward to solution independent, fully-implicit specific interface | ||
| 134 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 34273600 times.
|
96362830 | return asImp_().sourceAtPos(scv.dofPosition()); |
| 135 | } | ||
| 136 | |||
| 137 | /*! | ||
| 138 | * \brief Evaluate the source term for all phases within a given | ||
| 139 | * sub-control-volume. | ||
| 140 | * | ||
| 141 | * \param globalPos The position of the center of the finite volume | ||
| 142 | * for which the source term ought to be | ||
| 143 | * specified in global coordinates | ||
| 144 | * | ||
| 145 | * For this method, the values parameter stores the conserved quantity rate | ||
| 146 | * generated or annihilate per volume unit. Positive values mean | ||
| 147 | * that the conserved quantity is created, negative ones mean that it vanishes. | ||
| 148 | * E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$. | ||
| 149 | */ | ||
| 150 | 45849870 | Sources sourceAtPos(const GlobalPosition& globalPos) const | |
| 151 | { | ||
| 152 | //! As a default, i.e. if the user's problem does not overload any source method | ||
| 153 | //! return 0.0 (no source terms) | ||
| 154 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 45849870 times.
|
91699740 | return Sources(0.0); |
| 155 | } | ||
| 156 | |||
| 157 | /*! | ||
| 158 | * \brief Adds contribution of point sources for a specific sub control volume | ||
| 159 | * to the values. | ||
| 160 | * Caution: Only overload this method in the implementation if you know | ||
| 161 | * what you are doing. | ||
| 162 | */ | ||
| 163 | ✗ | Sources scvPointSources(const Element& element, | |
| 164 | const FVElementGeometry& fvGeometry, | ||
| 165 | const ElementVolumeVariables& elemVolVars, | ||
| 166 | const SubControlVolume& scv) const | ||
| 167 | { | ||
| 168 | ✗ | if (!this->pointSourceMap().empty()) | |
| 169 | ✗ | DUNE_THROW(Dune::NotImplemented, "scvPointSources not implemented"); | |
| 170 | |||
| 171 | ✗ | return Sources(0.0); | |
| 172 | } | ||
| 173 | |||
| 174 | /*! | ||
| 175 | * \brief Specifies which kind of boundary condition should be | ||
| 176 | * used for which equation on a given boundary segment. | ||
| 177 | * | ||
| 178 | * \param element The finite element | ||
| 179 | * \param scvf The sub control volume face | ||
| 180 | */ | ||
| 181 |
1/4✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
1455352 | auto boundaryTypes(const Element& element, |
| 182 | const SubControlVolumeFace& scvf) const | ||
| 183 | { | ||
| 184 | // Forward it to the method which only takes the global coordinate. | ||
| 185 | // We evaluate the boundary type at the center of the sub control volume face | ||
| 186 | // in order to avoid ambiguities at domain corners. | ||
| 187 |
16/16✓ Branch 0 taken 180806 times.
✓ Branch 1 taken 253344 times.
✓ Branch 2 taken 197992 times.
✓ Branch 3 taken 268588 times.
✓ Branch 4 taken 26900 times.
✓ Branch 5 taken 165396 times.
✓ Branch 6 taken 59296 times.
✓ Branch 7 taken 82880 times.
✓ Branch 8 taken 108992 times.
✓ Branch 9 taken 45632 times.
✓ Branch 10 taken 36172 times.
✓ Branch 11 taken 23876 times.
✓ Branch 12 taken 2436 times.
✓ Branch 13 taken 1500 times.
✓ Branch 14 taken 240 times.
✓ Branch 15 taken 1302 times.
|
1455352 | return asImp_().boundaryTypesAtPos(scvf.center()); |
| 188 | } | ||
| 189 | |||
| 190 | /*! | ||
| 191 | * \brief Evaluate the boundary conditions for a dirichlet | ||
| 192 | * control volume face. | ||
| 193 | * | ||
| 194 | * \param element The finite element | ||
| 195 | * \param scvf the sub control volume face | ||
| 196 | * \note used for cell-centered discretization schemes | ||
| 197 | */ | ||
| 198 |
6/14✓ Branch 0 taken 61740 times.
✓ Branch 1 taken 185220 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 taken 3136 times.
✓ Branch 9 taken 9408 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 6630 times.
✓ Branch 13 taken 19482 times.
|
1326644 | DirichletValues dirichlet(const Element& element, const SubControlVolumeFace& scvf) const |
| 199 | { | ||
| 200 |
9/15✓ Branch 0 taken 61740 times.
✓ Branch 1 taken 191266 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 5600 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 920 times.
✓ Branch 9 taken 3136 times.
✓ Branch 10 taken 9408 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 6630 times.
✓ Branch 14 taken 19482 times.
✓ Branch 2 taken 6046 times.
|
1600646 | return asImp_().dirichletAtPos(scvf.ipGlobal()); |
| 201 | } | ||
| 202 | |||
| 203 | /*! | ||
| 204 | * \brief Evaluates the boundary conditions for a Neumann control volume. | ||
| 205 | * | ||
| 206 | * \param element The element for which the Neumann boundary condition is set | ||
| 207 | * \param fvGeometry The fvGeometry | ||
| 208 | * \param elemVolVars The element volume variables | ||
| 209 | * \param elemFluxVarsCache The element flux variables cache | ||
| 210 | * \param scvf The boundary sub control volume face | ||
| 211 | */ | ||
| 212 | template<class ElementFluxVariablesCache> | ||
| 213 | ✗ | BoundaryFluxes neumann(const Element& element, | |
| 214 | const FVElementGeometry& fvGeometry, | ||
| 215 | const ElementVolumeVariables& elemVolVars, | ||
| 216 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 217 | const SubControlVolumeFace& scvf) const | ||
| 218 | ✗ | { return asImp_().neumannAtPos(scvf.ipGlobal()); } | |
| 219 | |||
| 220 | /*! | ||
| 221 | * \brief Returns the neumann flux at a given position. | ||
| 222 | */ | ||
| 223 | ✗ | BoundaryFluxes neumannAtPos(const GlobalPosition& globalPos) const | |
| 224 | ✗ | { return BoundaryFluxes(0.0); } //! A default, i.e. if the user's does not overload any neumann method | |
| 225 | |||
| 226 | /*! | ||
| 227 | * \brief Returns the acceleration due to gravity. | ||
| 228 | * | ||
| 229 | * If the <tt>Problem.EnableGravity</tt> parameter is true, this means | ||
| 230 | * \f$\boldsymbol{g} = ( 0,\dots,\ -9.81)^T \f$, else \f$\boldsymbol{g} = ( 0,\dots, 0)^T \f$ | ||
| 231 | */ | ||
| 232 | const GravityVector& gravity() const | ||
| 233 | 104220774 | { return gravity_; } | |
| 234 | |||
| 235 | /*! | ||
| 236 | * \brief Returns whether inertia terms should be considered. | ||
| 237 | */ | ||
| 238 | 343077434 | bool enableInertiaTerms() const | |
| 239 |
8/10✓ Branch 0 taken 106164576 times.
✓ Branch 1 taken 201463822 times.
✓ Branch 2 taken 118980 times.
✓ Branch 3 taken 192034 times.
✓ Branch 4 taken 34371500 times.
✓ Branch 5 taken 634630 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 65984 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 65908 times.
|
343077434 | { return enableInertiaTerms_; } |
| 240 | |||
| 241 | /*! | ||
| 242 | * \brief Returns the pressure at a given sub control volume face. | ||
| 243 | * \note Overload this if a fixed pressure shall be prescribed (e.g., given by an analytical solution). | ||
| 244 | */ | ||
| 245 | 104054148 | Scalar pressure(const Element& element, | |
| 246 | const FVElementGeometry& fvGeometry, | ||
| 247 | const SubControlVolumeFace& scvf) const | ||
| 248 | { | ||
| 249 | if constexpr (isCoupled_) | ||
| 250 | 103819484 | return couplingManager_->pressure(element, fvGeometry, scvf); | |
| 251 | else | ||
| 252 | 234664 | return asImp_().pressureAtPos(scvf.ipGlobal()); | |
| 253 | } | ||
| 254 | |||
| 255 | /*! | ||
| 256 | * \brief Returns the pressure at a given position. | ||
| 257 | */ | ||
| 258 | Scalar pressureAtPos(const GlobalPosition&) const | ||
| 259 | { | ||
| 260 | DUNE_THROW(Dune::NotImplemented, "pressureAtPos not implemented"); | ||
| 261 | } | ||
| 262 | |||
| 263 | /*! | ||
| 264 | * \brief Returns a reference pressure at a given sub control volume face. | ||
| 265 | * This pressure is subtracted from the actual pressure for the momentum balance | ||
| 266 | * which potentially helps to improve numerical accuracy by avoiding issues related do floating point arithmetic. | ||
| 267 | * \note Overload this for reference pressures other than zero. | ||
| 268 | */ | ||
| 269 | Scalar referencePressure(const Element& element, | ||
| 270 | const FVElementGeometry& fvGeometry, | ||
| 271 | const SubControlVolumeFace& scvf) const | ||
| 272 | { return 0.0; } | ||
| 273 | |||
| 274 | /*! | ||
| 275 | * \brief Returns the density at a given sub control volume face. | ||
| 276 | * \note Overload this if a fixed density shall be prescribed. | ||
| 277 | */ | ||
| 278 | 67194798 | Scalar density(const Element& element, | |
| 279 | const FVElementGeometry& fvGeometry, | ||
| 280 | const SubControlVolumeFace& scvf) const | ||
| 281 | { | ||
| 282 | if constexpr (isCoupled_) | ||
| 283 |
2/2✓ Branch 1 taken 3 times.
✓ Branch 2 taken 7355627 times.
|
67194798 | return couplingManager_->density(element, fvGeometry, scvf); |
| 284 | else | ||
| 285 | return asImp_().densityAtPos(scvf.ipGlobal()); | ||
| 286 | } | ||
| 287 | |||
| 288 | /*! | ||
| 289 | * \brief Returns the density at a given sub control volume. | ||
| 290 | * \note Overload this if a fixed density shall be prescribed. | ||
| 291 | */ | ||
| 292 | 226036228 | Scalar density(const Element& element, | |
| 293 | const SubControlVolume& scv, | ||
| 294 | const bool isPreviousTimeStep = false) const | ||
| 295 | { | ||
| 296 | if constexpr (isCoupled_) | ||
| 297 |
0/2✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
226036228 | return couplingManager_->density(element, scv, isPreviousTimeStep); |
| 298 | else | ||
| 299 | return asImp_().densityAtPos(scv.dofPosition()); | ||
| 300 | } | ||
| 301 | |||
| 302 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
132768324 | auto insideAndOutsideDensity(const Element& element, |
| 303 | const FVElementGeometry& fvGeometry, | ||
| 304 | const SubControlVolumeFace& scvf, | ||
| 305 | const bool isPreviousTimeStep = false) const | ||
| 306 | { | ||
| 307 | if constexpr (isCoupled_) | ||
| 308 |
4/7✓ Branch 1 taken 3 times.
✓ Branch 2 taken 14538917 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 3 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 6 taken 30437 times.
|
132768324 | return couplingManager_->insideAndOutsideDensity(element, fvGeometry, scvf, isPreviousTimeStep); |
| 309 | else | ||
| 310 | { | ||
| 311 | ✗ | const auto rho = asImp_().densityAtPos(scvf.ipGlobal()); | |
| 312 | ✗ | return std::make_pair(rho, rho); | |
| 313 | } | ||
| 314 | } | ||
| 315 | |||
| 316 | /*! | ||
| 317 | * \brief Returns the density at a given position. | ||
| 318 | */ | ||
| 319 | Scalar densityAtPos(const GlobalPosition&) const | ||
| 320 | { | ||
| 321 | DUNE_THROW(Dune::NotImplemented, "densityAtPos not implemented"); | ||
| 322 | } | ||
| 323 | |||
| 324 | /*! | ||
| 325 | * \brief Returns the effective dynamic viscosity at a given sub control volume face. | ||
| 326 | * \note Overload this if a fixed viscosity shall be prescribed. | ||
| 327 | */ | ||
| 328 | 309689520 | Scalar effectiveViscosity(const Element& element, | |
| 329 | const FVElementGeometry& fvGeometry, | ||
| 330 | const SubControlVolumeFace& scvf) const | ||
| 331 | { | ||
| 332 | if constexpr (isCoupled_) | ||
| 333 |
4/5✓ Branch 5 taken 217488 times.
✓ Branch 6 taken 192552 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 65908 times.
✓ Branch 7 taken 66968 times.
|
309682152 | return couplingManager_->effectiveViscosity(element, fvGeometry, scvf); |
| 334 | else | ||
| 335 | return asImp_().effectiveViscosityAtPos(scvf.ipGlobal()); | ||
| 336 | } | ||
| 337 | |||
| 338 | /*! | ||
| 339 | * \brief Returns the effective dynamic viscosity at a given sub control volume. | ||
| 340 | * \note Overload this if a fixed viscosity shall be prescribed. | ||
| 341 | */ | ||
| 342 | 90344 | Scalar effectiveViscosity(const Element& element, | |
| 343 | const FVElementGeometry& fvGeometry, | ||
| 344 | const SubControlVolume& scv) const | ||
| 345 | { | ||
| 346 | if constexpr (isCoupled_) | ||
| 347 | 90344 | return couplingManager_->effectiveViscosity(element, fvGeometry, scv); | |
| 348 | else | ||
| 349 | return asImp_().effectiveViscosityAtPos(scv.dofPosition()); | ||
| 350 | } | ||
| 351 | |||
| 352 | /*! | ||
| 353 | * \brief Returns the effective dynamic viscosity at a given position. | ||
| 354 | */ | ||
| 355 | Scalar effectiveViscosityAtPos(const GlobalPosition&) const | ||
| 356 | { | ||
| 357 | DUNE_THROW(Dune::NotImplemented, "effectiveViscosityAtPos not implemented"); | ||
| 358 | } | ||
| 359 | |||
| 360 | /*! | ||
| 361 | * \brief Applies the initial solution for all degrees of freedom of the grid. | ||
| 362 | * \param sol the initial solution vector | ||
| 363 | */ | ||
| 364 | template<class SolutionVector> | ||
| 365 | 33 | void applyInitialSolution(SolutionVector& sol) const | |
| 366 | { | ||
| 367 | 33 | sol.resize(this->gridGeometry().numDofs()); | |
| 368 |
1/2✓ Branch 3 taken 32 times.
✗ Branch 4 not taken.
|
33 | std::vector<bool> dofHandled(this->gridGeometry().numDofs(), false); |
| 369 |
1/2✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
|
33 | auto fvGeometry = localView(this->gridGeometry()); |
| 370 |
3/6✓ Branch 1 taken 33 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 199255 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1250 times.
✗ Branch 9 not taken.
|
597732 | for (const auto& element : elements(this->gridGeometry().gridView())) |
| 371 | { | ||
| 372 |
1/2✓ Branch 1 taken 1250 times.
✗ Branch 2 not taken.
|
199222 | fvGeometry.bindElement(element); |
| 373 |
4/4✓ Branch 0 taken 402729 times.
✓ Branch 1 taken 394159 times.
✓ Branch 2 taken 796888 times.
✓ Branch 3 taken 199222 times.
|
996110 | for (const auto& scv : scvs(fvGeometry)) |
| 374 | { | ||
| 375 | 796888 | const auto dofIdx = scv.dofIndex(); | |
| 376 |
2/2✓ Branch 0 taken 402729 times.
✓ Branch 1 taken 394159 times.
|
796888 | if (!dofHandled[dofIdx]) |
| 377 | { | ||
| 378 |
2/3✓ Branch 1 taken 81250 times.
✗ Branch 2 not taken.
✓ Branch 0 taken 5150 times.
|
402729 | dofHandled[dofIdx] = true; |
| 379 |
1/2✓ Branch 1 taken 10200 times.
✗ Branch 2 not taken.
|
478929 | sol[dofIdx] = asImp_().initial(scv)[scv.dofAxis()]; |
| 380 | } | ||
| 381 | } | ||
| 382 | } | ||
| 383 | 33 | } | |
| 384 | |||
| 385 | /*! | ||
| 386 | * \brief Evaluate the initial value at a sub control volume | ||
| 387 | */ | ||
| 388 |
2/2✓ Branch 0 taken 5150 times.
✓ Branch 1 taken 71050 times.
|
323205 | InitialValues initial(const SubControlVolume& scv) const |
| 389 | { | ||
| 390 |
2/2✓ Branch 0 taken 5150 times.
✓ Branch 1 taken 71050 times.
|
466869 | return asImp_().initialAtPos(scv.dofPosition()); |
| 391 | } | ||
| 392 | |||
| 393 | //! Convenience function for staggered grid implementation. | ||
| 394 | 7368 | Scalar pseudo3DWallFriction(const Element& element, | |
| 395 | const FVElementGeometry& fvGeometry, | ||
| 396 | const ElementVolumeVariables& elemVolVars, | ||
| 397 | const SubControlVolume& scv, | ||
| 398 | const Scalar height, | ||
| 399 | const Scalar factor = 8.0) const | ||
| 400 | { | ||
| 401 | 7368 | const Scalar velocity = elemVolVars[scv].velocity(); | |
| 402 | 7368 | const auto scvf = scvfs(fvGeometry, scv).begin(); | |
| 403 | 7368 | const Scalar viscosity = effectiveViscosity(element, fvGeometry, *scvf); | |
| 404 | 7368 | return pseudo3DWallFriction(velocity, viscosity, height, factor); | |
| 405 | } | ||
| 406 | |||
| 407 | /*! | ||
| 408 | * \brief An additional drag term can be included as source term for the momentum balance | ||
| 409 | * to mimic 3D flow behavior in 2D: | ||
| 410 | * \f[ | ||
| 411 | * f_{drag} = -(8 \mu / h^2)v | ||
| 412 | * \f] | ||
| 413 | * Here, \f$h\f$ corresponds to the extruded height that is | ||
| 414 | * bounded by the imaginary walls. See Flekkoy et al. (1995) \cite flekkoy1995a<BR> | ||
| 415 | * A value of 8.0 is used as a default factor, corresponding | ||
| 416 | * to the velocity profile at the center plane | ||
| 417 | * of the virtual height (maximum velocity). Setting this value to 12.0 corresponds | ||
| 418 | * to an depth-averaged velocity (Venturoli and Boek, 2006) \cite venturoli2006a. | ||
| 419 | */ | ||
| 420 | 7368 | Scalar pseudo3DWallFriction(const Scalar velocity, | |
| 421 | const Scalar viscosity, | ||
| 422 | const Scalar height, | ||
| 423 | const Scalar factor = 8.0) const | ||
| 424 | { | ||
| 425 | static_assert(dim == 2, "Pseudo 3D wall friction may only be used in 2D"); | ||
| 426 | 7368 | return -factor * velocity * viscosity / (height*height); | |
| 427 | } | ||
| 428 | |||
| 429 | /*! | ||
| 430 | * \brief Returns true if the scvf is located on a boundary with a slip condition. | ||
| 431 | * \note This member function must be overloaded in the problem implementation. | ||
| 432 | * Make sure to use scvf.center() for querying the spatial location. | ||
| 433 | */ | ||
| 434 | bool onSlipBoundary(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const | ||
| 435 | 338484 | { return asImp_().onSlipBoundaryAtPos(scvf.center()); } | |
| 436 | |||
| 437 | /*! | ||
| 438 | * \brief Returns true if the scvf is located on a boundary with a slip condition. | ||
| 439 | * \note This member function must be overloaded in the problem implementation. | ||
| 440 | */ | ||
| 441 | bool onSlipBoundaryAtPos(const GlobalPosition& pos) const | ||
| 442 | { return false; } | ||
| 443 | |||
| 444 | const CouplingManager& couplingManager() const | ||
| 445 | { | ||
| 446 | if constexpr (isCoupled_) | ||
| 447 | return *couplingManager_; | ||
| 448 | else | ||
| 449 | DUNE_THROW(Dune::InvalidStateException, | ||
| 450 | "Accessing coupling manager of an uncoupled problem is not possible." | ||
| 451 | ); | ||
| 452 | } | ||
| 453 | |||
| 454 | private: | ||
| 455 | //! Returns the implementation of the problem (i.e. static polymorphism) | ||
| 456 | Implementation& asImp_() | ||
| 457 | { return *static_cast<Implementation *>(this); } | ||
| 458 | |||
| 459 | //! \copydoc asImp_() | ||
| 460 | const Implementation& asImp_() const | ||
| 461 | { return *static_cast<const Implementation *>(this); } | ||
| 462 | |||
| 463 | GravityVector gravity_; | ||
| 464 | bool enableInertiaTerms_; | ||
| 465 | std::shared_ptr<CouplingManager> couplingManager_; | ||
| 466 | }; | ||
| 467 | |||
| 468 | /*! | ||
| 469 | * \brief Navier-Stokes default problem implementation for CVFE discretizations | ||
| 470 | */ | ||
| 471 | template<class TypeTag, class DM> | ||
| 472 | class NavierStokesMomentumProblemImpl<TypeTag, DiscretizationMethods::CVFE<DM>> | ||
| 473 | : public FVProblemWithSpatialParams<TypeTag> | ||
| 474 | { | ||
| 475 | using ParentType = FVProblemWithSpatialParams<TypeTag>; | ||
| 476 | using Implementation = GetPropType<TypeTag, Properties::Problem>; | ||
| 477 | |||
| 478 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 479 | using GridView = typename GridGeometry::GridView; | ||
| 480 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 481 | |||
| 482 | using GridVariables = GetPropType<TypeTag, Properties::GridVariables>; | ||
| 483 | using GridVolumeVariables = typename GridVariables::GridVolumeVariables; | ||
| 484 | using ElementVolumeVariables = typename GridVolumeVariables::LocalView; | ||
| 485 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 486 | |||
| 487 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 488 | using SubControlVolume = typename FVElementGeometry::SubControlVolume; | ||
| 489 | using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; | ||
| 490 | using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices; | ||
| 491 | |||
| 492 | static constexpr int dim = GridView::dimension; | ||
| 493 | static constexpr int dimWorld = GridView::dimensionworld; | ||
| 494 | |||
| 495 | using GlobalPosition = typename SubControlVolumeFace::GlobalPosition; | ||
| 496 | using VelocityVector = Dune::FieldVector<Scalar, dimWorld>; | ||
| 497 | using GravityVector = Dune::FieldVector<Scalar, dimWorld>; | ||
| 498 | using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>; | ||
| 499 | using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>; | ||
| 500 | |||
| 501 | public: | ||
| 502 | //! These types are used in place of the typical NumEqVector type. | ||
| 503 | //! In the numeqvector assembly type, only one equation per DOF (face) is considered | ||
| 504 | //! while the type here provides one entry for each world dimension. | ||
| 505 | using InitialValues = Dune::FieldVector<Scalar, dimWorld>; | ||
| 506 | using Sources = Dune::FieldVector<Scalar, dimWorld>; | ||
| 507 | using DirichletValues = Dune::FieldVector<Scalar, dimWorld>; | ||
| 508 | using BoundaryFluxes = Dune::FieldVector<Scalar, dimWorld>; | ||
| 509 | |||
| 510 | //! Export the boundary types. | ||
| 511 | using BoundaryTypes = NavierStokesMomentumBoundaryTypes<ModelTraits::dim()>; | ||
| 512 | |||
| 513 | //! This problem is used for the momentum balance model. | ||
| 514 | static constexpr bool isMomentumProblem() { return true; } | ||
| 515 | |||
| 516 | /*! | ||
| 517 | * \brief The constructor | ||
| 518 | * \param gridGeometry The finite volume grid geometry | ||
| 519 | * \param couplingManager A coupling manager (in case this problem is a coupled "mass-momentum"/full Navier-Stokes problem) | ||
| 520 | * \param paramGroup The parameter group in which to look for runtime parameters first (default is "") | ||
| 521 | */ | ||
| 522 | 30 | NavierStokesMomentumProblemImpl(std::shared_ptr<const GridGeometry> gridGeometry, | |
| 523 | std::shared_ptr<CouplingManager> couplingManager, | ||
| 524 | const std::string& paramGroup = "") | ||
| 525 | : ParentType(gridGeometry, paramGroup) | ||
| 526 | 60 | , gravity_(0.0) | |
| 527 |
1/2✓ Branch 2 taken 30 times.
✗ Branch 3 not taken.
|
90 | , couplingManager_(couplingManager) |
| 528 | { | ||
| 529 |
2/4✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 30 times.
|
30 | if (getParamFromGroup<bool>(paramGroup, "Problem.EnableGravity")) |
| 530 | ✗ | gravity_[dim-1] = -9.81; | |
| 531 | |||
| 532 |
1/2✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
|
30 | enableInertiaTerms_ = getParamFromGroup<bool>(paramGroup, "Problem.EnableInertiaTerms"); |
| 533 | 30 | } | |
| 534 | |||
| 535 | /*! | ||
| 536 | * \brief The constructor for usage without a coupling manager | ||
| 537 | * \param gridGeometry The finite volume grid geometry | ||
| 538 | * \param paramGroup The parameter group in which to look for runtime parameters first (default is "") | ||
| 539 | */ | ||
| 540 | 18 | NavierStokesMomentumProblemImpl(std::shared_ptr<const GridGeometry> gridGeometry, | |
| 541 | const std::string& paramGroup = "") | ||
| 542 |
1/2✓ Branch 2 taken 18 times.
✗ Branch 3 not taken.
|
36 | : NavierStokesMomentumProblemImpl(gridGeometry, {}, paramGroup) |
| 543 | 18 | {} | |
| 544 | |||
| 545 | /*! | ||
| 546 | * \brief Evaluate the source term for all phases within a given | ||
| 547 | * sub-control-volume. | ||
| 548 | * | ||
| 549 | * This is the method for the case where the source term is | ||
| 550 | * potentially solution dependent and requires some quantities that | ||
| 551 | * are specific to the fully-implicit method. | ||
| 552 | * | ||
| 553 | * \param element The finite element | ||
| 554 | * \param fvGeometry The finite-volume geometry | ||
| 555 | * \param elemVolVars All volume variables for the element | ||
| 556 | * \param scv The sub control volume | ||
| 557 | * | ||
| 558 | * For this method, the return parameter stores the conserved quantity rate | ||
| 559 | * generated or annihilate per volume unit. Positive values mean | ||
| 560 | * that the conserved quantity is created, negative ones mean that it vanishes. | ||
| 561 | * E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$. | ||
| 562 | */ | ||
| 563 | template<class ElementVolumeVariables> | ||
| 564 | 31591694 | Sources source(const Element &element, | |
| 565 | const FVElementGeometry& fvGeometry, | ||
| 566 | const ElementVolumeVariables& elemVolVars, | ||
| 567 | const SubControlVolume &scv) const | ||
| 568 | { | ||
| 569 | // forward to solution independent, fully-implicit specific interface | ||
| 570 | 31591694 | return asImp_().sourceAtPos(scv.center()); | |
| 571 | } | ||
| 572 | |||
| 573 | /*! | ||
| 574 | * \brief Evaluate the source term for all phases at a given position | ||
| 575 | * | ||
| 576 | * \param globalPos The position of the center of the finite volume | ||
| 577 | * for which the source term ought to be | ||
| 578 | * specified in global coordinates | ||
| 579 | * | ||
| 580 | * For this method, the values parameter stores the conserved quantity rate | ||
| 581 | * generated or annihilate per volume unit. Positive values mean | ||
| 582 | * that the conserved quantity is created, negative ones mean that it vanishes. | ||
| 583 | * E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$. | ||
| 584 | */ | ||
| 585 | 30436494 | Sources sourceAtPos(const GlobalPosition &globalPos) const | |
| 586 | { | ||
| 587 | //! As a default, i.e. if the user's problem does not overload any source method | ||
| 588 | //! return 0.0 (no source terms) | ||
| 589 |
2/2✓ Branch 0 taken 5976666 times.
✓ Branch 1 taken 1992222 times.
|
66849654 | return Sources(0.0); |
| 590 | } | ||
| 591 | |||
| 592 | /*! | ||
| 593 | * \brief Specifies which kind of boundary condition should be | ||
| 594 | * used for which equation on a given boundary segment. | ||
| 595 | * | ||
| 596 | * \param element The finite element | ||
| 597 | * \param scv The sub control volume | ||
| 598 | */ | ||
| 599 | 72813 | BoundaryTypes boundaryTypes(const Element& element, | |
| 600 | const SubControlVolume& scv) const | ||
| 601 | { | ||
| 602 | // forward it to the method which only takes the global coordinate | ||
| 603 |
2/2✓ Branch 0 taken 47412 times.
✓ Branch 1 taken 21130 times.
|
72813 | return asImp_().boundaryTypesAtPos(scv.dofPosition()); |
| 604 | } | ||
| 605 | |||
| 606 | /*! | ||
| 607 | * \brief Specifies which kind of boundary condition should be | ||
| 608 | * used for which equation on a given boundary segment. | ||
| 609 | * | ||
| 610 | * \param element The finite element | ||
| 611 | * \param scvf The sub control volume face | ||
| 612 | */ | ||
| 613 | BoundaryTypes boundaryTypes(const Element& element, | ||
| 614 | const SubControlVolumeFace& scvf) const | ||
| 615 | { | ||
| 616 | DUNE_THROW(Dune::InvalidStateException, "boundaryTypes(..., scvf) called for a CVFE method."); | ||
| 617 | } | ||
| 618 | |||
| 619 | /*! | ||
| 620 | * \brief Evaluate the boundary conditions for a Dirichlet | ||
| 621 | * control volume. | ||
| 622 | * | ||
| 623 | * \param element The finite element | ||
| 624 | * \param scv the sub control volume | ||
| 625 | */ | ||
| 626 | 86156 | DirichletValues dirichlet(const Element& element, const SubControlVolume& scv) const | |
| 627 | { | ||
| 628 | // forward it to the method which only takes the global coordinate | ||
| 629 |
4/4✓ Branch 0 taken 47274 times.
✓ Branch 1 taken 15758 times.
✓ Branch 2 taken 35556 times.
✓ Branch 3 taken 11852 times.
|
197676 | return asImp_().dirichletAtPos(scv.dofPosition()); |
| 630 | } | ||
| 631 | |||
| 632 | /*! | ||
| 633 | * \brief Evaluate the boundary conditions for a Dirichlet | ||
| 634 | * control volume face. | ||
| 635 | * | ||
| 636 | * \param element The finite element | ||
| 637 | * \param scvf the sub control volume face | ||
| 638 | */ | ||
| 639 | DirichletValues dirichlet(const Element& element, const SubControlVolumeFace& scvf) const | ||
| 640 | { | ||
| 641 | // forward it to the method which only takes the global coordinate | ||
| 642 | DUNE_THROW(Dune::InvalidStateException, "dirichlet(scvf) called for CVFE method."); | ||
| 643 | } | ||
| 644 | |||
| 645 | /*! | ||
| 646 | * \brief Returns the acceleration due to gravity. | ||
| 647 | * | ||
| 648 | * If the <tt>Problem.EnableGravity</tt> parameter is true, this means | ||
| 649 | * \f$\boldsymbol{g} = ( 0,\dots,\ -9.81)^T \f$, else \f$\boldsymbol{g} = ( 0,\dots, 0)^T \f$ | ||
| 650 | */ | ||
| 651 | const GravityVector& gravity() const | ||
| 652 | 32455694 | { return gravity_; } | |
| 653 | |||
| 654 | /*! | ||
| 655 | * \brief Returns whether inertia terms should be considered. | ||
| 656 | */ | ||
| 657 | 46901758 | bool enableInertiaTerms() const | |
| 658 |
6/6✓ Branch 0 taken 17582868 times.
✓ Branch 1 taken 29298040 times.
✓ Branch 2 taken 17370 times.
✓ Branch 3 taken 3474 times.
✓ Branch 4 taken 3 times.
✓ Branch 5 taken 3 times.
|
46901758 | { return enableInertiaTerms_; } |
| 659 | |||
| 660 | /*! | ||
| 661 | * \brief Returns a reference pressure | ||
| 662 | * This pressure is subtracted from the actual pressure for the momentum balance | ||
| 663 | * which potentially helps to improve numerical accuracy by avoiding issues related do floating point arithmetic. | ||
| 664 | * \note Overload this for reference pressures other than zero. | ||
| 665 | */ | ||
| 666 | Scalar referencePressure() const | ||
| 667 | { return 0.0; } | ||
| 668 | |||
| 669 | /*! | ||
| 670 | * \brief Returns the pressure at a given sub control volume face. | ||
| 671 | * \note Overload this if a fixed pressure shall be prescribed (e.g., given by an analytical solution). | ||
| 672 | */ | ||
| 673 | 46921262 | Scalar pressure(const Element& element, | |
| 674 | const FVElementGeometry& fvGeometry, | ||
| 675 | const SubControlVolumeFace& scvf) const | ||
| 676 | { | ||
| 677 | if constexpr (std::is_empty_v<CouplingManager>) | ||
| 678 | 6334400 | return asImp_().pressureAtPos(scvf.ipGlobal()); | |
| 679 | else | ||
| 680 |
1/2✓ Branch 2 taken 40354 times.
✗ Branch 3 not taken.
|
40586862 | return couplingManager_->pressure(element, fvGeometry, scvf); |
| 681 | } | ||
| 682 | |||
| 683 | /*! | ||
| 684 | * \brief Returns the pressure at a given sub control volume. | ||
| 685 | * \note Overload this if a fixed pressure shall be prescribed (e.g., given by an analytical solution). | ||
| 686 | */ | ||
| 687 | 3822000 | Scalar pressure(const Element& element, | |
| 688 | const FVElementGeometry& fvGeometry, | ||
| 689 | const SubControlVolume& scv, | ||
| 690 | const bool isPreviousTimeStep = false) const | ||
| 691 | { | ||
| 692 | if constexpr (std::is_empty_v<CouplingManager>) | ||
| 693 | 3822000 | return asImp_().pressureAtPos(scv.dofPosition()); | |
| 694 | else | ||
| 695 | return couplingManager_->pressure(element, fvGeometry, scv, isPreviousTimeStep); | ||
| 696 | } | ||
| 697 | |||
| 698 | /*! | ||
| 699 | * \brief Returns the pressure at a given position. | ||
| 700 | */ | ||
| 701 | Scalar pressureAtPos(const GlobalPosition&) const | ||
| 702 | { | ||
| 703 | DUNE_THROW(Dune::NotImplemented, "pressureAtPos not implemented"); | ||
| 704 | } | ||
| 705 | |||
| 706 | /*! | ||
| 707 | * \brief Returns the density at a given sub control volume face. | ||
| 708 | * \note Overload this if a fixed density shall be prescribed. | ||
| 709 | */ | ||
| 710 | ✗ | Scalar density(const Element& element, | |
| 711 | const FVElementGeometry& fvGeometry, | ||
| 712 | const SubControlVolumeFace& scvf) const | ||
| 713 | { | ||
| 714 | if constexpr (std::is_empty_v<CouplingManager>) | ||
| 715 | ✗ | return asImp_().densityAtPos(scvf.ipGlobal()); | |
| 716 | else | ||
| 717 | ✗ | return couplingManager_->density(element, fvGeometry, scvf); | |
| 718 | } | ||
| 719 | |||
| 720 | /*! | ||
| 721 | * \brief Returns the density at a given sub control volume. | ||
| 722 | * \note Overload this if a fixed density shall be prescribed. | ||
| 723 | */ | ||
| 724 | 5292054 | Scalar density(const Element& element, | |
| 725 | const FVElementGeometry& fvGeometry, | ||
| 726 | const SubControlVolume& scv, | ||
| 727 | const bool isPreviousTimeStep = false) const | ||
| 728 | { | ||
| 729 | if constexpr (std::is_empty_v<CouplingManager>) | ||
| 730 | 3822000 | return asImp_().densityAtPos(scv.dofPosition()); | |
| 731 | else | ||
| 732 | 1470054 | return couplingManager_->density(element, fvGeometry, scv, isPreviousTimeStep); | |
| 733 | } | ||
| 734 | |||
| 735 | |||
| 736 | /*! | ||
| 737 | * \brief Returns the density at a given position. | ||
| 738 | */ | ||
| 739 | Scalar densityAtPos(const GlobalPosition&) const | ||
| 740 | { | ||
| 741 | DUNE_THROW(Dune::NotImplemented, "densityAtPos not implemented"); | ||
| 742 | } | ||
| 743 | |||
| 744 | /*! | ||
| 745 | * \brief Returns the effective dynamic viscosity at a given sub control volume face. | ||
| 746 | * \note Overload this if a fixed viscosity shall be prescribed. | ||
| 747 | */ | ||
| 748 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 5207990 times.
|
7471188 | Scalar effectiveViscosity(const Element& element, |
| 749 | const FVElementGeometry& fvGeometry, | ||
| 750 | const SubControlVolumeFace& scvf) const | ||
| 751 | { | ||
| 752 | if constexpr (std::is_empty_v<CouplingManager>) | ||
| 753 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 5207990 times.
|
5208000 | return asImp_().effectiveViscosityAtPos(scvf.ipGlobal()); |
| 754 | else | ||
| 755 | 2263188 | return couplingManager_->effectiveViscosity(element, fvGeometry, scvf); | |
| 756 | } | ||
| 757 | |||
| 758 | /*! | ||
| 759 | * \brief Returns the effective dynamic viscosity at a given sub control volume. | ||
| 760 | * \note Overload this if a fixed viscosity shall be prescribed. | ||
| 761 | */ | ||
| 762 | 4066500 | Scalar effectiveViscosity(const Element& element, | |
| 763 | const FVElementGeometry& fvGeometry, | ||
| 764 | const SubControlVolume& scv, | ||
| 765 | const bool isPreviousTimeStep = false) const | ||
| 766 | { | ||
| 767 | if constexpr (std::is_empty_v<CouplingManager>) | ||
| 768 | 3822000 | return asImp_().effectiveViscosityAtPos(scv.dofPosition()); | |
| 769 | else | ||
| 770 | 244500 | return couplingManager_->effectiveViscosity(element, fvGeometry, scv, isPreviousTimeStep); | |
| 771 | } | ||
| 772 | |||
| 773 | /*! | ||
| 774 | * \brief Returns the effective dynamic viscosity at a given position. | ||
| 775 | */ | ||
| 776 | Scalar effectiveViscosityAtPos(const GlobalPosition&) const | ||
| 777 | { | ||
| 778 | DUNE_THROW(Dune::NotImplemented, "effectiveViscosityAtPos not implemented"); | ||
| 779 | } | ||
| 780 | |||
| 781 | /*! | ||
| 782 | * \brief Applies the initial solution for all degrees of freedom of the grid. | ||
| 783 | * \param sol the initial solution vector | ||
| 784 | */ | ||
| 785 | template<class SolutionVector> | ||
| 786 | void applyInitialSolution(SolutionVector& sol) const | ||
| 787 | { | ||
| 788 | static_assert(GridGeometry::discMethod == DiscretizationMethods::CVFE<DM>{}); | ||
| 789 | sol.resize(this->gridGeometry().numDofs()); | ||
| 790 | std::vector<bool> dofHandled(this->gridGeometry().numDofs(), false); | ||
| 791 | auto fvGeometry = localView(this->gridGeometry()); | ||
| 792 | for (const auto& element : elements(this->gridGeometry().gridView())) | ||
| 793 | { | ||
| 794 | fvGeometry.bindElement(element); | ||
| 795 | for (const auto& scv : scvs(fvGeometry)) | ||
| 796 | { | ||
| 797 | const auto dofIdx = scv.dofIndex(); | ||
| 798 | if (!dofHandled[dofIdx]) | ||
| 799 | { | ||
| 800 | dofHandled[dofIdx] = true; | ||
| 801 | sol[dofIdx] = asImp_().initial(scv); | ||
| 802 | } | ||
| 803 | } | ||
| 804 | } | ||
| 805 | } | ||
| 806 | |||
| 807 | /*! | ||
| 808 | * \brief Evaluate the initial value at an sub control volume | ||
| 809 | */ | ||
| 810 | InitialValues initial(const SubControlVolume& scv) const | ||
| 811 | { | ||
| 812 | static_assert(GridGeometry::discMethod == DiscretizationMethods::CVFE<DM>{}); | ||
| 813 | return asImp_().initialAtPos(scv.dofPosition()); | ||
| 814 | } | ||
| 815 | |||
| 816 | private: | ||
| 817 | //! Returns the implementation of the problem (i.e. static polymorphism) | ||
| 818 | Implementation &asImp_() | ||
| 819 | { return *static_cast<Implementation *>(this); } | ||
| 820 | |||
| 821 | //! \copydoc asImp_() | ||
| 822 | const Implementation &asImp_() const | ||
| 823 | { return *static_cast<const Implementation *>(this); } | ||
| 824 | |||
| 825 | GravityVector gravity_; | ||
| 826 | bool enableInertiaTerms_; | ||
| 827 | std::shared_ptr<CouplingManager> couplingManager_ = {}; | ||
| 828 | }; | ||
| 829 | |||
| 830 | /*! | ||
| 831 | * \ingroup NavierStokesModel | ||
| 832 | * \brief Navier-Stokes momentum problem class | ||
| 833 | * | ||
| 834 | * Inherit from this problem to implement Navier-Stokes momentum problems | ||
| 835 | */ | ||
| 836 | template<class TypeTag> | ||
| 837 | using NavierStokesMomentumProblem = NavierStokesMomentumProblemImpl< | ||
| 838 | TypeTag, typename GetPropType<TypeTag, Properties::GridGeometry>::DiscretizationMethod | ||
| 839 | >; | ||
| 840 | |||
| 841 | } // end namespace Dumux | ||
| 842 | |||
| 843 | #endif | ||
| 844 |