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