| 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 Flux | ||
| 10 | * \brief Forchheimer's law | ||
| 11 | * This file contains the calculation of the Forchheimer velocity | ||
| 12 | * for a given Darcy velocity | ||
| 13 | */ | ||
| 14 | #ifndef DUMUX_FLUX_FORCHHEIMER_VELOCITY_HH | ||
| 15 | #define DUMUX_FLUX_FORCHHEIMER_VELOCITY_HH | ||
| 16 | |||
| 17 | #include <dune/common/fvector.hh> | ||
| 18 | #include <dune/common/fmatrix.hh> | ||
| 19 | |||
| 20 | #include <dumux/common/math.hh> | ||
| 21 | #include <dumux/common/parameters.hh> | ||
| 22 | #include <dumux/common/properties.hh> | ||
| 23 | #include <dumux/common/typetraits/typetraits.hh> | ||
| 24 | |||
| 25 | #include <dumux/flux/facetensoraverage.hh> | ||
| 26 | |||
| 27 | namespace Dumux { | ||
| 28 | |||
| 29 | /*! | ||
| 30 | * \ingroup Flux | ||
| 31 | * \brief Forchheimer's law | ||
| 32 | * For a detailed description see dumux/flow/forchheimerslaw.hh | ||
| 33 | * | ||
| 34 | * This file contains the calculation of the Forchheimer velocity | ||
| 35 | * for a given Darcy velocity | ||
| 36 | */ | ||
| 37 | template<class Scalar, class GridGeometry, class FluxVariables> | ||
| 38 | class ForchheimerVelocity | ||
| 39 | { | ||
| 40 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 41 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 42 | using GridView = typename GridGeometry::GridView; | ||
| 43 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 44 | |||
| 45 | static constexpr int dim = GridView::dimension; | ||
| 46 | static constexpr int dimWorld = GridView::dimensionworld; | ||
| 47 | using GlobalPosition = typename Element::Geometry::GlobalCoordinate; | ||
| 48 | public: | ||
| 49 | using UpwindScheme = typename FluxVariables::UpwindScheme; | ||
| 50 | using DimWorldVector = Dune::FieldVector<Scalar, dimWorld>; | ||
| 51 | using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>; | ||
| 52 | |||
| 53 | |||
| 54 | template<class ElementVolumeVariables> | ||
| 55 | 3568 | static DimWorldVector velocity(const FVElementGeometry& fvGeometry, | |
| 56 | const ElementVolumeVariables& elemVolVars, | ||
| 57 | const SubControlVolumeFace& scvf, | ||
| 58 | const int phaseIdx, | ||
| 59 | const DimWorldMatrix sqrtK, | ||
| 60 | const DimWorldVector darcyVelocity) | ||
| 61 | { | ||
| 62 | 3568 | const auto& problem = elemVolVars.gridVolVars().problem(); | |
| 63 | |||
| 64 | // Obtain the Forchheimer coefficient from the spatial parameters | ||
| 65 | 3568 | const Scalar forchCoeff = problem.spatialParams().forchCoeff(scvf); | |
| 66 | |||
| 67 | // Initial guess of velocity: Darcy relation | ||
| 68 | 3568 | DimWorldVector velocity = darcyVelocity; | |
| 69 | |||
| 70 | 3568 | DimWorldVector deltaV(0.0); // the change in velocity between Newton iterations | |
| 71 | 3568 | DimWorldVector residual(10e10); // the residual (function value that is to be minimized) | |
| 72 | 3568 | DimWorldMatrix gradF(0.0); // slope of equation that is to be solved | |
| 73 | |||
| 74 | // Search by means of the Newton method for a root of Forchheimer equation | ||
| 75 |
5/8✓ Branch 0 taken 2 times.
✓ Branch 1 taken 3566 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
|
3568 | static const Scalar epsilon = getParamFromGroup<Scalar>(problem.paramGroup(), "Forchheimer.NewtonTolerance", 1e-12); |
| 76 |
5/8✓ Branch 0 taken 2 times.
✓ Branch 1 taken 3566 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
|
3568 | static const std::size_t maxNumIter = getParamFromGroup<std::size_t>(problem.paramGroup(), "Forchheimer.MaxIterations", 30); |
| 77 |
2/2✓ Branch 0 taken 25004 times.
✓ Branch 1 taken 3568 times.
|
57144 | for (int k = 0; residual.two_norm() > epsilon ; ++k) |
| 78 | { | ||
| 79 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 25004 times.
|
25004 | if (k >= maxNumIter) |
| 80 | ✗ | DUNE_THROW(NumericalProblem, "could not determine forchheimer velocity within " | |
| 81 | << k << " iterations"); | ||
| 82 | |||
| 83 | // calculate current value of Forchheimer equation | ||
| 84 | 25004 | forchheimerResidual_(residual, | |
| 85 | forchCoeff, | ||
| 86 | sqrtK, | ||
| 87 | darcyVelocity, | ||
| 88 | velocity, | ||
| 89 | elemVolVars, | ||
| 90 | scvf, | ||
| 91 | phaseIdx); | ||
| 92 | |||
| 93 | // newton's method: slope (gradF) and current value (residual) of function is known, | ||
| 94 | 25004 | forchheimerDerivative_(gradF, | |
| 95 | forchCoeff, | ||
| 96 | sqrtK, | ||
| 97 | velocity, | ||
| 98 | elemVolVars, | ||
| 99 | scvf, | ||
| 100 | phaseIdx); | ||
| 101 | |||
| 102 | // solve for change in velocity ("x-Axis intercept") | ||
| 103 | 25004 | gradF.solve(deltaV, residual); | |
| 104 | 25004 | velocity -= deltaV; | |
| 105 | } | ||
| 106 | |||
| 107 | // The fluxvariables expect a value on which an upwinding of the mobility is performed. | ||
| 108 | // We therefore divide by the mobility here. | ||
| 109 | 14272 | auto upwindTerm = [phaseIdx](const auto& volVars){ return volVars.mobility(phaseIdx); }; | |
| 110 | 12912 | const Scalar forchheimerUpwindMobility = UpwindScheme::multiplier(elemVolVars, scvf, upwindTerm, velocity * scvf.unitOuterNormal(), phaseIdx); | |
| 111 | |||
| 112 | // Do not divide by zero. If the mobility is zero, the resulting flux with upwinding will be zero anyway. | ||
| 113 |
1/2✓ Branch 0 taken 3568 times.
✗ Branch 1 not taken.
|
3568 | if (forchheimerUpwindMobility > 0.0) |
| 114 | velocity /= forchheimerUpwindMobility; | ||
| 115 | |||
| 116 | 3568 | return velocity; | |
| 117 | } | ||
| 118 | |||
| 119 | /*! \brief Returns the harmonic mean of \f$\sqrt{K_0}\f$ and \f$\sqrt{K_1}\f$. | ||
| 120 | * | ||
| 121 | * This is a specialization for scalar-valued permeabilities which returns a tensor with identical diagonal entries. | ||
| 122 | */ | ||
| 123 | template<class Problem, class ElementVolumeVariables, | ||
| 124 | bool scalarPerm = std::is_same<typename Problem::SpatialParams::PermeabilityType, Scalar>::value, | ||
| 125 | std::enable_if_t<scalarPerm, int> = 0> | ||
| 126 | 4080 | static DimWorldMatrix calculateHarmonicMeanSqrtPermeability(const Problem& problem, | |
| 127 | const ElementVolumeVariables& elemVolVars, | ||
| 128 | const SubControlVolumeFace& scvf) | ||
| 129 | { | ||
| 130 | using std::sqrt; | ||
| 131 | 4080 | Scalar harmonicMeanSqrtK(0.0); | |
| 132 | |||
| 133 |
1/2✓ Branch 0 taken 1360 times.
✗ Branch 1 not taken.
|
4080 | const auto insideScvIdx = scvf.insideScvIdx(); |
| 134 |
1/2✓ Branch 0 taken 1360 times.
✗ Branch 1 not taken.
|
4080 | const auto& insideVolVars = elemVolVars[insideScvIdx]; |
| 135 |
3/4✓ Branch 0 taken 3600 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 1360 times.
✗ Branch 3 not taken.
|
5440 | const Scalar Ki = getPermeability_(problem, insideVolVars, scvf.ipGlobal()); |
| 136 | 4080 | const Scalar sqrtKi = sqrt(Ki); | |
| 137 | |||
| 138 | 4080 | if (!scvf.boundary()) | |
| 139 | { | ||
| 140 |
1/2✓ Branch 0 taken 1360 times.
✗ Branch 1 not taken.
|
3600 | const auto outsideScvIdx = scvf.outsideScvIdx(); |
| 141 |
1/2✓ Branch 0 taken 1360 times.
✗ Branch 1 not taken.
|
3600 | const auto& outsideVolVars = elemVolVars[outsideScvIdx]; |
| 142 | 7200 | const Scalar Kj = getPermeability_(problem, outsideVolVars, scvf.ipGlobal()); | |
| 143 | 3600 | const Scalar sqrtKj = sqrt(Kj); | |
| 144 |
1/2✓ Branch 0 taken 3600 times.
✗ Branch 1 not taken.
|
3600 | harmonicMeanSqrtK = faceTensorAverage(sqrtKi, sqrtKj, scvf.unitOuterNormal()); |
| 145 | } | ||
| 146 | else | ||
| 147 | harmonicMeanSqrtK = sqrtKi; | ||
| 148 | |||
| 149 | // Convert to tensor | ||
| 150 | 4080 | DimWorldMatrix result(0.0); | |
| 151 |
2/2✓ Branch 0 taken 8160 times.
✓ Branch 1 taken 4080 times.
|
12240 | for (int i = 0; i < dimWorld; ++i) |
| 152 | 24480 | result[i][i] = harmonicMeanSqrtK; | |
| 153 | |||
| 154 | 4080 | return result; | |
| 155 | } | ||
| 156 | |||
| 157 | /*! \brief Returns the harmonic mean of \f$\sqrt{\mathbf{K_0}}\f$ and \f$\sqrt{\mathbf{K_1}}\f$. | ||
| 158 | * | ||
| 159 | * This is a specialization for tensor-valued permeabilities. | ||
| 160 | */ | ||
| 161 | template<class Problem, class ElementVolumeVariables, | ||
| 162 | bool scalarPerm = std::is_same<typename Problem::SpatialParams::PermeabilityType, Scalar>::value, | ||
| 163 | std::enable_if_t<!scalarPerm, int> = 0> | ||
| 164 | static DimWorldMatrix calculateHarmonicMeanSqrtPermeability(const Problem& problem, | ||
| 165 | const ElementVolumeVariables& elemVolVars, | ||
| 166 | const SubControlVolumeFace& scvf) | ||
| 167 | { | ||
| 168 | using std::sqrt; | ||
| 169 | DimWorldMatrix sqrtKi(0.0); | ||
| 170 | DimWorldMatrix sqrtKj(0.0); | ||
| 171 | DimWorldMatrix harmonicMeanSqrtK(0.0); | ||
| 172 | |||
| 173 | const auto insideScvIdx = scvf.insideScvIdx(); | ||
| 174 | const auto& insideVolVars = elemVolVars[insideScvIdx]; | ||
| 175 | const auto& Ki = getPermeability_(problem, insideVolVars, scvf.ipGlobal()); | ||
| 176 | // Make sure the permeability matrix does not have off-diagonal entries | ||
| 177 | // TODO implement method using eigenvalues and eigenvectors for general tensors | ||
| 178 | if (!isDiagonal_(Ki)) | ||
| 179 | DUNE_THROW(Dune::NotImplemented, "Only diagonal permeability tensors are supported."); | ||
| 180 | |||
| 181 | for (int i = 0; i < dim; ++i) | ||
| 182 | sqrtKi[i][i] = sqrt(Ki[i][i]); | ||
| 183 | |||
| 184 | if (!scvf.boundary()) | ||
| 185 | { | ||
| 186 | const auto outsideScvIdx = scvf.outsideScvIdx(); | ||
| 187 | const auto& outsideVolVars = elemVolVars[outsideScvIdx]; | ||
| 188 | const auto& Kj = getPermeability_(problem, outsideVolVars, scvf.ipGlobal()); | ||
| 189 | // Make sure the permeability matrix does not have off-diagonal entries | ||
| 190 | if (!isDiagonal_(Kj)) | ||
| 191 | DUNE_THROW(Dune::NotImplemented, "Only diagonal permeability tensors are supported."); | ||
| 192 | |||
| 193 | for (int i = 0; i < dim; ++i) | ||
| 194 | sqrtKj[i][i] = sqrt(Kj[i][i]); | ||
| 195 | |||
| 196 | harmonicMeanSqrtK = faceTensorAverage(sqrtKi, sqrtKj, scvf.unitOuterNormal()); | ||
| 197 | } | ||
| 198 | else | ||
| 199 | harmonicMeanSqrtK = sqrtKi; | ||
| 200 | |||
| 201 | return harmonicMeanSqrtK; | ||
| 202 | } | ||
| 203 | |||
| 204 | private: | ||
| 205 | |||
| 206 | //! calculate the residual of the Forchheimer equation | ||
| 207 | template<class ElementVolumeVariables> | ||
| 208 | 25004 | static void forchheimerResidual_(DimWorldVector& residual, | |
| 209 | const Scalar forchCoeff, | ||
| 210 | const DimWorldMatrix& sqrtK, | ||
| 211 | const DimWorldVector& darcyVelocity, | ||
| 212 | const DimWorldVector& oldVelocity, | ||
| 213 | const ElementVolumeVariables& elemVolVars, | ||
| 214 | const SubControlVolumeFace& scvf, | ||
| 215 | const int phaseIdx) | ||
| 216 | { | ||
| 217 | 25004 | residual = oldVelocity; | |
| 218 | // residual += k_r/mu K grad p | ||
| 219 | 25004 | residual -= darcyVelocity; | |
| 220 | // residual += c_F rho / mu abs(v) sqrt(K) v | ||
| 221 | 150024 | auto upwindTerm = [phaseIdx](const auto& volVars){ return volVars.density(phaseIdx) / volVars.viscosity(phaseIdx); }; | |
| 222 | 75012 | const Scalar rhoOverMu = UpwindScheme::multiplier(elemVolVars, scvf, upwindTerm, oldVelocity * scvf.unitOuterNormal(), phaseIdx); | |
| 223 | 50008 | sqrtK.usmv(forchCoeff * rhoOverMu * oldVelocity.two_norm(), oldVelocity, residual); | |
| 224 | 25004 | } | |
| 225 | |||
| 226 | //! The analytical derivative of the the Forcheimer equation's residual | ||
| 227 | template<class ElementVolumeVariables> | ||
| 228 | 25004 | static void forchheimerDerivative_(DimWorldMatrix& derivative, | |
| 229 | const Scalar forchCoeff, | ||
| 230 | const DimWorldMatrix& sqrtK, | ||
| 231 | const DimWorldVector& velocity, | ||
| 232 | const ElementVolumeVariables& elemVolVars, | ||
| 233 | const SubControlVolumeFace& scvf, | ||
| 234 | const int phaseIdx) | ||
| 235 | { | ||
| 236 | // Initialize because for low velocities, we add and do not set the entries. | ||
| 237 | 25004 | derivative = 0.0; | |
| 238 | |||
| 239 | 150024 | auto upwindTerm = [phaseIdx](const auto& volVars){ return volVars.density(phaseIdx) / volVars.viscosity(phaseIdx); }; | |
| 240 | |||
| 241 | 25004 | const Scalar absV = velocity.two_norm(); | |
| 242 | // This part of the derivative is only calculated if the velocity is sufficiently small: | ||
| 243 | // do not divide by zero! | ||
| 244 | // This should not be very bad because the derivative is only intended to give an | ||
| 245 | // approximation of the gradient of the | ||
| 246 | // function that goes into the Newton scheme. | ||
| 247 | // This is important in the case of a one-phase region in two-phase flow. The non-existing | ||
| 248 | // phase has a velocity of zero (kr=0). | ||
| 249 | // f = sqrtK * vTimesV (this is a matrix) | ||
| 250 | // f *= forchCoeff density / |velocity| / viscosity (multiply all entries with scalars) | ||
| 251 | 75012 | const Scalar rhoOverMu = UpwindScheme::multiplier(elemVolVars, scvf, upwindTerm, velocity * scvf.unitOuterNormal(), phaseIdx); | |
| 252 |
2/2✓ Branch 0 taken 24646 times.
✓ Branch 1 taken 358 times.
|
25004 | if (absV > 1e-20) |
| 253 | { | ||
| 254 |
2/2✓ Branch 0 taken 49292 times.
✓ Branch 1 taken 24646 times.
|
73938 | for (int i = 0; i < dim; i++) |
| 255 | { | ||
| 256 |
2/2✓ Branch 0 taken 73938 times.
✓ Branch 1 taken 49292 times.
|
123230 | for (int k = 0; k <= i; k++) |
| 257 | { | ||
| 258 |
10/10✓ Branch 0 taken 24646 times.
✓ Branch 1 taken 49292 times.
✓ Branch 2 taken 24646 times.
✓ Branch 3 taken 49292 times.
✓ Branch 4 taken 24646 times.
✓ Branch 5 taken 49292 times.
✓ Branch 6 taken 24646 times.
✓ Branch 7 taken 49292 times.
✓ Branch 8 taken 24646 times.
✓ Branch 9 taken 49292 times.
|
369690 | derivative[i][k] = sqrtK[i][i] * velocity[i] * velocity[k] * forchCoeff |
| 259 |
2/2✓ Branch 0 taken 24646 times.
✓ Branch 1 taken 49292 times.
|
73938 | / absV * rhoOverMu; |
| 260 | |||
| 261 |
2/2✓ Branch 0 taken 24646 times.
✓ Branch 1 taken 49292 times.
|
73938 | if (k < i) |
| 262 | 123230 | derivative[k][i] = derivative[i][k]; | |
| 263 | } | ||
| 264 | } | ||
| 265 | } | ||
| 266 | |||
| 267 | // add on the main diagonal: | ||
| 268 | // 1 + sqrtK_i forchCoeff density |v| / viscosity | ||
| 269 |
2/2✓ Branch 0 taken 50008 times.
✓ Branch 1 taken 25004 times.
|
75012 | for (int i = 0; i < dim; i++) |
| 270 | 250040 | derivative[i][i] += (1.0 + (sqrtK[i][i]*forchCoeff * absV * rhoOverMu)); | |
| 271 | 25004 | } | |
| 272 | |||
| 273 | /*! | ||
| 274 | * \brief Check whether all off-diagonal entries of a tensor are zero. | ||
| 275 | * | ||
| 276 | * \param K the tensor that is to be checked. | ||
| 277 | * | ||
| 278 | * \return True if all off-diagonals are zero. | ||
| 279 | * | ||
| 280 | */ | ||
| 281 | static bool isDiagonal_(const DimWorldMatrix & K) | ||
| 282 | { | ||
| 283 | using std::abs; | ||
| 284 | for (int i = 0; i < dim; i++) | ||
| 285 | for (int k = 0; k < dim; k++) | ||
| 286 | if ((i != k) && (abs(K[i][k]) >= 1e-25)) | ||
| 287 | return false; | ||
| 288 | return true; | ||
| 289 | } | ||
| 290 | |||
| 291 | template<class Problem, class VolumeVariables> | ||
| 292 | ✗ | static decltype(auto) getPermeability_(const Problem& problem, | |
| 293 | const VolumeVariables& volVars, | ||
| 294 | const GlobalPosition& scvfIpGlobal) | ||
| 295 | { | ||
| 296 | if constexpr (Problem::SpatialParams::evaluatePermeabilityAtScvfIP()) | ||
| 297 | return problem.spatialParams().permeabilityAtPos(scvfIpGlobal); | ||
| 298 | else | ||
| 299 |
3/4✓ Branch 0 taken 3600 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 3600 times.
✗ Branch 3 not taken.
|
6320 | return volVars.permeability(); |
| 300 | } | ||
| 301 | }; | ||
| 302 | |||
| 303 | } // end namespace Dumux | ||
| 304 | |||
| 305 | #endif | ||
| 306 |