| 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 BoundaryTests | ||
| 10 | * \brief The Darcy sub-problem of coupled Stokes-Darcy convergence test | ||
| 11 | */ | ||
| 12 | |||
| 13 | #ifndef DUMUX_DARCY_SUBPROBLEM_HH | ||
| 14 | #define DUMUX_DARCY_SUBPROBLEM_HH | ||
| 15 | |||
| 16 | #include <dune/common/fvector.hh> | ||
| 17 | #include <dune/grid/yaspgrid.hh> | ||
| 18 | |||
| 19 | #include <dumux/discretization/cctpfa.hh> | ||
| 20 | |||
| 21 | #include <dumux/common/boundarytypes.hh> | ||
| 22 | #include <dumux/common/numeqvector.hh> | ||
| 23 | |||
| 24 | #include <dumux/material/components/constant.hh> | ||
| 25 | #include <dumux/material/fluidsystems/1pliquid.hh> | ||
| 26 | |||
| 27 | #include <dumux/porousmediumflow/1p/model.hh> | ||
| 28 | #include <dumux/porousmediumflow/problem.hh> | ||
| 29 | |||
| 30 | #include "spatialparams.hh" | ||
| 31 | #include "testcase.hh" | ||
| 32 | |||
| 33 | namespace Dumux { | ||
| 34 | template <class TypeTag> | ||
| 35 | class DarcySubProblem; | ||
| 36 | |||
| 37 | /*! | ||
| 38 | * \ingroup BoundaryTests | ||
| 39 | * \brief The Darcy sub-problem of coupled Stokes-Darcy convergence test | ||
| 40 | */ | ||
| 41 | template <class TypeTag> | ||
| 42 | class DarcySubProblem : public PorousMediumFlowProblem<TypeTag> | ||
| 43 | { | ||
| 44 | using ParentType = PorousMediumFlowProblem<TypeTag>; | ||
| 45 | using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView; | ||
| 46 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 47 | using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; | ||
| 48 | using NumEqVector = Dumux::NumEqVector<PrimaryVariables>; | ||
| 49 | using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>; | ||
| 50 | using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>; | ||
| 51 | using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView; | ||
| 52 | using SubControlVolume = typename FVElementGeometry::SubControlVolume; | ||
| 53 | using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; | ||
| 54 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 55 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 56 | using GlobalPosition = typename Element::Geometry::GlobalCoordinate; | ||
| 57 | |||
| 58 | using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>; | ||
| 59 | |||
| 60 | static constexpr auto velocityXIdx = 0; | ||
| 61 | static constexpr auto velocityYIdx = 1; | ||
| 62 | static constexpr auto pressureIdx = 2; | ||
| 63 | |||
| 64 | public: | ||
| 65 | //! export the Indices | ||
| 66 | using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices; | ||
| 67 | |||
| 68 | 6 | DarcySubProblem(std::shared_ptr<const GridGeometry> gridGeometry, | |
| 69 | std::shared_ptr<CouplingManager> couplingManager, | ||
| 70 | std::shared_ptr<typename ParentType::SpatialParams> spatialParams, | ||
| 71 | const TestCase testCase, | ||
| 72 | const std::string& name) | ||
| 73 | : ParentType(gridGeometry, spatialParams, "Darcy") | ||
| 74 | , couplingManager_(couplingManager) | ||
| 75 |
7/22✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 6 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 6 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 6 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 6 times.
✓ Branch 20 taken 6 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
|
24 | , testCase_(testCase) |
| 76 | { | ||
| 77 |
7/20✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 6 times.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✓ Branch 14 taken 6 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 6 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 6 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
|
12 | problemName_ = name + "_" + getParamFromGroup<std::string>(this->paramGroup(), "Problem.Name"); |
| 78 | 6 | } | |
| 79 | |||
| 80 | /*! | ||
| 81 | * \brief The problem name. | ||
| 82 | */ | ||
| 83 | const std::string& name() const | ||
| 84 | { | ||
| 85 |
2/4✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
|
12 | return problemName_; |
| 86 | } | ||
| 87 | |||
| 88 | /*! | ||
| 89 | * \name Boundary conditions | ||
| 90 | */ | ||
| 91 | // \{ | ||
| 92 | |||
| 93 | /*! | ||
| 94 | * \brief Specifies which kind of boundary condition should be | ||
| 95 | * used for which equation on a given boundary control volume. | ||
| 96 | * | ||
| 97 | * \param element The element | ||
| 98 | * \param scvf The boundary sub control volume face | ||
| 99 | */ | ||
| 100 | BoundaryTypes boundaryTypes(const Element &element, const SubControlVolumeFace &scvf) const | ||
| 101 | { | ||
| 102 |
8/8✓ Branch 0 taken 1120 times.
✓ Branch 1 taken 3360 times.
✓ Branch 2 taken 1120 times.
✓ Branch 3 taken 3360 times.
✓ Branch 4 taken 11760 times.
✓ Branch 5 taken 11928 times.
✓ Branch 6 taken 3080 times.
✓ Branch 7 taken 9240 times.
|
44968 | BoundaryTypes values; |
| 103 | |||
| 104 |
8/8✓ Branch 0 taken 1120 times.
✓ Branch 1 taken 3360 times.
✓ Branch 2 taken 1120 times.
✓ Branch 3 taken 3360 times.
✓ Branch 4 taken 11760 times.
✓ Branch 5 taken 11928 times.
✓ Branch 6 taken 3080 times.
✓ Branch 7 taken 9240 times.
|
44968 | if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf)) |
| 105 | values.setAllCouplingNeumann(); | ||
| 106 | else | ||
| 107 | values.setAllDirichlet(); | ||
| 108 | |||
| 109 | return values; | ||
| 110 | } | ||
| 111 | |||
| 112 | /*! | ||
| 113 | * \brief Evaluates the boundary conditions for a Dirichlet control volume. | ||
| 114 | * | ||
| 115 | * \param element The element for which the Dirichlet boundary condition is set | ||
| 116 | * \param scvf The boundary subcontrolvolumeface | ||
| 117 | * | ||
| 118 | * For this method, the \a values parameter stores primary variables. | ||
| 119 | */ | ||
| 120 | ✗ | PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const | |
| 121 | { | ||
| 122 |
0/4✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
18480 | const auto p = analyticalSolution(scvf.center())[pressureIdx]; |
| 123 | 9240 | return PrimaryVariables(p); | |
| 124 | } | ||
| 125 | |||
| 126 | /*! | ||
| 127 | * \brief Evaluates the boundary conditions for a Neumann control volume. | ||
| 128 | * | ||
| 129 | * \param element The element for which the Neumann boundary condition is set | ||
| 130 | * \param fvGeometry The fvGeometry | ||
| 131 | * \param elemVolVars The element volume variables | ||
| 132 | * \param elemFluxVarsCache Flux variables caches for all faces in stencil | ||
| 133 | * \param scvf The boundary sub control volume face | ||
| 134 | * | ||
| 135 | * For this method, the \a values variable stores primary variables. | ||
| 136 | */ | ||
| 137 | template<class ElementVolumeVariables, class ElementFluxVarsCache> | ||
| 138 | ✗ | NumEqVector neumann(const Element& element, | |
| 139 | const FVElementGeometry& fvGeometry, | ||
| 140 | const ElementVolumeVariables& elemVolVars, | ||
| 141 | const ElementFluxVarsCache& elemFluxVarsCache, | ||
| 142 | const SubControlVolumeFace& scvf) const | ||
| 143 | { | ||
| 144 | ✗ | NumEqVector values(0.0); | |
| 145 | |||
| 146 | ✗ | if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf)) | |
| 147 | ✗ | values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf); | |
| 148 | |||
| 149 | ✗ | return values; | |
| 150 | } | ||
| 151 | |||
| 152 | // \} | ||
| 153 | |||
| 154 | /*! | ||
| 155 | * \name Volume terms | ||
| 156 | */ | ||
| 157 | // \{ | ||
| 158 | /*! | ||
| 159 | * \brief Evaluates the source term for all phases within a given | ||
| 160 | * sub control volume. | ||
| 161 | * | ||
| 162 | * \param globalPos The global position | ||
| 163 | */ | ||
| 164 | 478240 | NumEqVector sourceAtPos(const GlobalPosition& globalPos) const | |
| 165 | { | ||
| 166 |
2/5✗ Branch 0 not taken.
✓ Branch 1 taken 136640 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 341600 times.
✗ Branch 4 not taken.
|
478240 | switch (testCase_) |
| 167 | { | ||
| 168 | ✗ | case TestCase::ShiueExampleOne: | |
| 169 | ✗ | return rhsShiueEtAlExampleOne_(globalPos); | |
| 170 | case TestCase::ShiueExampleTwo: | ||
| 171 | 136640 | return rhsShiueEtAlExampleTwo_(globalPos); | |
| 172 | ✗ | case TestCase::Rybak: | |
| 173 | ✗ | return rhsRybak_(globalPos); | |
| 174 | 341600 | case TestCase::Schneider: | |
| 175 | 341600 | return rhsSchneiderEtAl_(globalPos); | |
| 176 | ✗ | default: | |
| 177 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid test case"); | |
| 178 | } | ||
| 179 | } | ||
| 180 | |||
| 181 | // \} | ||
| 182 | |||
| 183 | /*! | ||
| 184 | * \brief Evaluates the initial value for a control volume. | ||
| 185 | * | ||
| 186 | * \param element The element | ||
| 187 | * | ||
| 188 | * For this method, the \a priVars parameter stores primary | ||
| 189 | * variables. | ||
| 190 | */ | ||
| 191 | PrimaryVariables initial(const Element &element) const | ||
| 192 | { | ||
| 193 | return PrimaryVariables(0.0); | ||
| 194 | } | ||
| 195 | |||
| 196 | /*! | ||
| 197 | * \brief Returns the analytical solution of the problem at a given position. | ||
| 198 | * | ||
| 199 | * \param globalPos The global position | ||
| 200 | */ | ||
| 201 | 143640 | auto analyticalSolution(const GlobalPosition& globalPos) const | |
| 202 | { | ||
| 203 |
2/5✗ Branch 0 not taken.
✓ Branch 1 taken 70560 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 73080 times.
✗ Branch 4 not taken.
|
143640 | switch (testCase_) |
| 204 | { | ||
| 205 | ✗ | case TestCase::ShiueExampleOne: | |
| 206 | ✗ | return analyticalSolutionShiueEtAlExampleOne_(globalPos); | |
| 207 | 70560 | case TestCase::ShiueExampleTwo: | |
| 208 | 70560 | return analyticalSolutionShiueEtAlExampleTwo_(globalPos); | |
| 209 | ✗ | case TestCase::Rybak: | |
| 210 | ✗ | return analyticalSolutionRybak_(globalPos); | |
| 211 | 73080 | case TestCase::Schneider: | |
| 212 | 73080 | return analyticalSolutionSchneiderEtAl_(globalPos); | |
| 213 | ✗ | default: | |
| 214 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid test case"); | |
| 215 | } | ||
| 216 | } | ||
| 217 | |||
| 218 | // \} | ||
| 219 | |||
| 220 | //! Get the coupling manager | ||
| 221 | const CouplingManager& couplingManager() const | ||
| 222 |
16/24✓ Branch 0 taken 1120 times.
✓ Branch 1 taken 3360 times.
✓ Branch 2 taken 1120 times.
✓ Branch 3 taken 3360 times.
✓ Branch 4 taken 1120 times.
✓ Branch 5 taken 3360 times.
✓ Branch 6 taken 1120 times.
✓ Branch 7 taken 3360 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✓ Branch 18 taken 11760 times.
✓ Branch 19 taken 11928 times.
✓ Branch 20 taken 11760 times.
✓ Branch 21 taken 11928 times.
✓ Branch 22 taken 3080 times.
✓ Branch 23 taken 9240 times.
✓ Branch 24 taken 3080 times.
✓ Branch 25 taken 9240 times.
|
89936 | { return *couplingManager_; } |
| 223 | |||
| 224 | private: | ||
| 225 | |||
| 226 | // see Rybak et al., 2015: "Multirate time integration for coupled saturated/unsaturated porous medium and free flow systems" | ||
| 227 | ✗ | Dune::FieldVector<Scalar, 3> analyticalSolutionRybak_(const GlobalPosition& globalPos) const | |
| 228 | { | ||
| 229 | ✗ | Dune::FieldVector<Scalar, 3> sol(0.0); | |
| 230 | ✗ | const Scalar x = globalPos[0]; | |
| 231 | ✗ | const Scalar y = globalPos[1]; | |
| 232 | |||
| 233 | using std::exp; using std::sin; using std::cos; | ||
| 234 | ✗ | sol[velocityXIdx] = -0.5*M_PI*y*y*cos(M_PI*x); | |
| 235 | ✗ | sol[velocityYIdx] = -1.0*y*sin(M_PI*x); | |
| 236 | ✗ | sol[pressureIdx] = 0.5*y*y*sin(M_PI*x); | |
| 237 | ✗ | return sol; | |
| 238 | } | ||
| 239 | |||
| 240 | // see Rybak et al., 2015: "Multirate time integration for coupled saturated/unsaturated porous medium and free flow systems" | ||
| 241 | ✗ | NumEqVector rhsRybak_(const GlobalPosition& globalPos) const | |
| 242 | { | ||
| 243 | ✗ | const Scalar x = globalPos[0]; | |
| 244 | ✗ | const Scalar y = globalPos[1]; | |
| 245 | using std::sin; | ||
| 246 | ✗ | return NumEqVector((0.5*M_PI*y*M_PI*y - 1)*sin(M_PI*x)); | |
| 247 | } | ||
| 248 | |||
| 249 | // see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem" | ||
| 250 | ✗ | Dune::FieldVector<Scalar, 3> analyticalSolutionShiueEtAlExampleOne_(const GlobalPosition& globalPos) const | |
| 251 | { | ||
| 252 | ✗ | Dune::FieldVector<Scalar, 3> sol(0.0); | |
| 253 | ✗ | const Scalar x = globalPos[0]; | |
| 254 | ✗ | const Scalar y = globalPos[1]; | |
| 255 | |||
| 256 | using std::exp; using std::sin; using std::cos; | ||
| 257 | ✗ | sol[pressureIdx] = (exp(y) - y*exp(1)) * cos(M_PI*x); | |
| 258 | ✗ | sol[velocityXIdx] = M_PI*(-exp(1)*y + exp(y))*sin(M_PI*x); | |
| 259 | ✗ | sol[velocityYIdx] = (exp(1) - exp(y))*cos(M_PI*x); | |
| 260 | |||
| 261 | ✗ | return sol; | |
| 262 | } | ||
| 263 | |||
| 264 | // see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem" | ||
| 265 | ✗ | NumEqVector rhsShiueEtAlExampleOne_(const GlobalPosition& globalPos) const | |
| 266 | { | ||
| 267 | ✗ | const Scalar x = globalPos[0]; | |
| 268 | ✗ | const Scalar y = globalPos[1]; | |
| 269 | using std::exp; using std::sin; using std::cos; | ||
| 270 | ✗ | return NumEqVector(cos(M_PI*x) * (M_PI*M_PI*(exp(y) - y*exp(1)) - exp(y))); | |
| 271 | } | ||
| 272 | |||
| 273 | // see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem" | ||
| 274 | 70560 | Dune::FieldVector<Scalar, 3> analyticalSolutionShiueEtAlExampleTwo_(const GlobalPosition& globalPos) const | |
| 275 | { | ||
| 276 | 70560 | Dune::FieldVector<Scalar, 3> sol(0.0); | |
| 277 | 141120 | const Scalar x = globalPos[0]; | |
| 278 | 141120 | const Scalar y = globalPos[1]; | |
| 279 | |||
| 280 | 141120 | sol[pressureIdx] = x*(1.0-x)*(y-1.0) + (y-1.0)*(y-1.0)*(y-1.0)/3.0 + 2.0*x + 2.0*y + 4.0; | |
| 281 | 141120 | sol[velocityXIdx] = x*(y - 1.0) + (x - 1.0)*(y - 1.0) - 2.0; | |
| 282 | 141120 | sol[velocityYIdx] = x*(x - 1.0) - 1.0*(y - 1.0)*(y - 1.0) - 2.0; | |
| 283 | |||
| 284 | 70560 | return sol; | |
| 285 | } | ||
| 286 | |||
| 287 | // see Shiue et al., 2018: "Convergence of the MAC Scheme for the Stokes/Darcy Coupling Problem" | ||
| 288 | ✗ | NumEqVector rhsShiueEtAlExampleTwo_(const GlobalPosition& globalPos) const | |
| 289 | 136640 | { return NumEqVector(0.0); } | |
| 290 | |||
| 291 | // see Schneider et al., 2019: "Coupling staggered-grid and MPFA finite volume methods for | ||
| 292 | // free flow/porous-medium flow problems" | ||
| 293 | 73080 | Dune::FieldVector<Scalar, 3> analyticalSolutionSchneiderEtAl_(const GlobalPosition& globalPos) const | |
| 294 | { | ||
| 295 | 73080 | Dune::FieldVector<Scalar, 3> sol(0.0); | |
| 296 | 146160 | const Scalar x = globalPos[0]; | |
| 297 | 146160 | const Scalar y = globalPos[1]; | |
| 298 | static constexpr Scalar omega = M_PI; | ||
| 299 | static constexpr Scalar c = 0.0; | ||
| 300 | using std::exp; using std::sin; using std::cos; | ||
| 301 | 73080 | const Scalar sinOmegaX = sin(omega*x); | |
| 302 | 73080 | const Scalar cosOmegaX = cos(omega*x); | |
| 303 | static const Scalar expTwo = exp(2); | ||
| 304 | 73080 | const Scalar expYPlusOne = exp(y+1); | |
| 305 | |||
| 306 | 146160 | sol[pressureIdx] = (expYPlusOne + 2 - expTwo)*sinOmegaX; | |
| 307 | 219240 | sol[velocityXIdx] = c/(2*omega)*expYPlusOne*sinOmegaX*sinOmegaX | |
| 308 | 73080 | -omega*(expYPlusOne + 2 - expTwo)*cosOmegaX; | |
| 309 | 219240 | sol[velocityYIdx] = (0.5*c*(expYPlusOne + 2 - expTwo)*cosOmegaX | |
| 310 | 73080 | -(c*cosOmegaX + 1)*exp(y-1))*sinOmegaX; | |
| 311 | |||
| 312 | 73080 | return sol; | |
| 313 | } | ||
| 314 | |||
| 315 | // see Schneider et al., 2019: "Coupling staggered-grid and MPFA finite volume methods for | ||
| 316 | // free flow/porous-medium flow problems (with c = 0)" | ||
| 317 | 341600 | NumEqVector rhsSchneiderEtAl_(const GlobalPosition& globalPos) const | |
| 318 | { | ||
| 319 | 683200 | const Scalar x = globalPos[0]; | |
| 320 | 683200 | const Scalar y = globalPos[1]; | |
| 321 | using std::exp; using std::sin; using std::cos; | ||
| 322 | static constexpr Scalar omega = M_PI; | ||
| 323 | static constexpr Scalar c = 0.0; | ||
| 324 | 341600 | const Scalar cosOmegaX = cos(omega*x); | |
| 325 | static const Scalar expTwo = exp(2); | ||
| 326 | 341600 | const Scalar expYPlusOne = exp(y+1); | |
| 327 | |||
| 328 | 683200 | const Scalar result = (-(c*cosOmegaX + 1)*exp(y - 1) | |
| 329 | 341600 | + 1.5*c*expYPlusOne*cosOmegaX | |
| 330 | 341600 | + omega*omega*(expYPlusOne - expTwo + 2)) | |
| 331 | 341600 | *sin(omega*x); | |
| 332 | 341600 | return NumEqVector(result); | |
| 333 | } | ||
| 334 | |||
| 335 | static constexpr Scalar eps_ = 1e-7; | ||
| 336 | std::shared_ptr<CouplingManager> couplingManager_; | ||
| 337 | std::string problemName_; | ||
| 338 | TestCase testCase_; | ||
| 339 | }; | ||
| 340 | } // end namespace Dumux | ||
| 341 | |||
| 342 | #endif | ||
| 343 |