| 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 OnePModel | ||
| 10 | * \brief Element-wise calculation of the residual and its derivatives | ||
| 11 | * for a single-phase, incompressible, test problem. | ||
| 12 | */ | ||
| 13 | |||
| 14 | #ifndef DUMUX_1P_INCOMPRESSIBLE_LOCAL_RESIDUAL_HH | ||
| 15 | #define DUMUX_1P_INCOMPRESSIBLE_LOCAL_RESIDUAL_HH | ||
| 16 | |||
| 17 | #include <dumux/discretization/method.hh> | ||
| 18 | #include <dumux/porousmediumflow/immiscible/localresidual.hh> | ||
| 19 | |||
| 20 | namespace Dumux { | ||
| 21 | |||
| 22 | /*! | ||
| 23 | * \ingroup OnePModel | ||
| 24 | * \brief Element-wise calculation of the residual and its derivatives | ||
| 25 | * for a single-phase, incompressible, test problem. | ||
| 26 | */ | ||
| 27 | template<class TypeTag> | ||
| 28 | class OnePIncompressibleLocalResidual : public ImmiscibleLocalResidual<TypeTag> | ||
| 29 | { | ||
| 30 | using ParentType = ImmiscibleLocalResidual<TypeTag>; | ||
| 31 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 32 | using Problem = GetPropType<TypeTag, Properties::Problem>; | ||
| 33 | using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>; | ||
| 34 | using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView; | ||
| 35 | using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>; | ||
| 36 | using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>; | ||
| 37 | using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView; | ||
| 38 | using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView; | ||
| 39 | using SubControlVolume = typename FVElementGeometry::SubControlVolume; | ||
| 40 | using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; | ||
| 41 | using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView; | ||
| 42 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 43 | using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices; | ||
| 44 | // first index for the mass balance | ||
| 45 | enum { conti0EqIdx = Indices::conti0EqIdx }; | ||
| 46 | enum { pressureIdx = Indices::pressureIdx }; | ||
| 47 | |||
| 48 | static constexpr int dim = GridView::dimension; | ||
| 49 | static constexpr int dimWorld = GridView::dimensionworld; | ||
| 50 | |||
| 51 | public: | ||
| 52 |
4/8✓ Branch 1 taken 199046 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 199046 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 921521 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 921521 times.
✗ Branch 11 not taken.
|
4455204 | using ParentType::ParentType; |
| 53 | |||
| 54 | template<class PartialDerivativeMatrix> | ||
| 55 | ✗ | void addStorageDerivatives(PartialDerivativeMatrix& partialDerivatives, | |
| 56 | const Problem& problem, | ||
| 57 | const Element& element, | ||
| 58 | const FVElementGeometry& fvGeometry, | ||
| 59 | const VolumeVariables& curVolVars, | ||
| 60 | ✗ | const SubControlVolume& scv) const {} | |
| 61 | |||
| 62 | template<class PartialDerivativeMatrix> | ||
| 63 | ✗ | void addSourceDerivatives(PartialDerivativeMatrix& partialDerivatives, | |
| 64 | const Problem& problem, | ||
| 65 | const Element& element, | ||
| 66 | const FVElementGeometry& fvGeometry, | ||
| 67 | const VolumeVariables& curVolVars, | ||
| 68 | const SubControlVolume& scv) const | ||
| 69 | { | ||
| 70 | 170900 | problem.addSourceDerivatives(partialDerivatives, element, fvGeometry, curVolVars, scv); | |
| 71 | ✗ | } | |
| 72 | |||
| 73 | //! Flux derivatives for the cell-centered tpfa scheme | ||
| 74 | template<class PartialDerivativeMatrices, class T = TypeTag> | ||
| 75 | std::enable_if_t<GetPropType<T, Properties::GridGeometry>::discMethod == DiscretizationMethods::cctpfa, void> | ||
| 76 | 205320 | addFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices, | |
| 77 | const Problem& problem, | ||
| 78 | const Element& element, | ||
| 79 | const FVElementGeometry& fvGeometry, | ||
| 80 | const ElementVolumeVariables& curElemVolVars, | ||
| 81 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 82 | const SubControlVolumeFace& scvf) const | ||
| 83 | { | ||
| 84 | static_assert(!FluidSystem::isCompressible(0), | ||
| 85 | "1p/incompressiblelocalresidual.hh: Only incompressible fluids are allowed!"); | ||
| 86 | static_assert(FluidSystem::viscosityIsConstant(0), | ||
| 87 | "1p/incompressiblelocalresidual.hh: Only fluids with constant viscosities are allowed!"); | ||
| 88 | |||
| 89 | // we know the "upwind factor" is constant, get inner one here and compute derivatives | ||
| 90 |
3/4✓ Branch 0 taken 12 times.
✓ Branch 1 taken 205308 times.
✓ Branch 3 taken 12 times.
✗ Branch 4 not taken.
|
205321 | static const Scalar up = curElemVolVars[scvf.insideScvIdx()].density() |
| 91 | 35 | / curElemVolVars[scvf.insideScvIdx()].viscosity(); | |
| 92 | |||
| 93 | // for network grids we have to do something else | ||
| 94 |
4/4✓ Branch 0 taken 588 times.
✓ Branch 1 taken 6100 times.
✓ Branch 2 taken 588 times.
✓ Branch 3 taken 6100 times.
|
13376 | if (dim < dimWorld && scvf.numOutsideScvs() > 1) |
| 95 | { | ||
| 96 | 588 | const auto insideTi = elemFluxVarsCache[scvf].advectionTij(); | |
| 97 | 588 | const auto sumTi = [&] | |
| 98 | { | ||
| 99 | 588 | Scalar sumTi(insideTi); | |
| 100 |
4/4✓ Branch 0 taken 1248 times.
✓ Branch 1 taken 588 times.
✓ Branch 2 taken 1248 times.
✓ Branch 3 taken 588 times.
|
3084 | for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i) |
| 101 | { | ||
| 102 | 1248 | const auto& flippedScvf = fvGeometry.flipScvf(scvf.index(), i); | |
| 103 | 1248 | const auto& outsideFluxVarsCache = elemFluxVarsCache[flippedScvf]; | |
| 104 | 1248 | sumTi += outsideFluxVarsCache.advectionTij(); | |
| 105 | } | ||
| 106 | 588 | return sumTi; | |
| 107 | 588 | }(); | |
| 108 | |||
| 109 | // add partial derivatives to the respective given matrices | ||
| 110 | 1176 | derivativeMatrices[scvf.insideScvIdx()][conti0EqIdx][pressureIdx] | |
| 111 | 588 | += up*insideTi*(1.0 - insideTi/sumTi); | |
| 112 |
4/4✓ Branch 0 taken 1248 times.
✓ Branch 1 taken 588 times.
✓ Branch 2 taken 1248 times.
✓ Branch 3 taken 588 times.
|
3084 | for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i) |
| 113 | 1248 | derivativeMatrices[scvf.outsideScvIdx(i)][conti0EqIdx][pressureIdx] | |
| 114 | 3744 | += -up*insideTi*elemFluxVarsCache[fvGeometry.flipScvf(scvf.index(), i)].advectionTij()/sumTi; | |
| 115 | |||
| 116 | } | ||
| 117 | else | ||
| 118 | { | ||
| 119 | 204732 | const auto deriv = elemFluxVarsCache[scvf].advectionTij()*up; | |
| 120 | // add partial derivatives to the respective given matrices | ||
| 121 | 409464 | derivativeMatrices[scvf.insideScvIdx()][conti0EqIdx][pressureIdx] += deriv; | |
| 122 | 409464 | derivativeMatrices[scvf.outsideScvIdx()][conti0EqIdx][pressureIdx] -= deriv; | |
| 123 | } | ||
| 124 | 205320 | } | |
| 125 | |||
| 126 | //! Flux derivatives for the cell-centered mpfa scheme | ||
| 127 | template<class PartialDerivativeMatrices, class T = TypeTag> | ||
| 128 | std::enable_if_t<GetPropType<T, Properties::GridGeometry>::discMethod == DiscretizationMethods::ccmpfa, void> | ||
| 129 | 2000 | addFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices, | |
| 130 | const Problem& problem, | ||
| 131 | const Element& element, | ||
| 132 | const FVElementGeometry& fvGeometry, | ||
| 133 | const ElementVolumeVariables& curElemVolVars, | ||
| 134 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 135 | const SubControlVolumeFace& scvf) const | ||
| 136 | { | ||
| 137 | static_assert(!FluidSystem::isCompressible(0), | ||
| 138 | "1p/incompressiblelocalresidual.hh: Only incompressible fluids are allowed!"); | ||
| 139 | static_assert(FluidSystem::viscosityIsConstant(0), | ||
| 140 | "1p/incompressiblelocalresidual.hh: Only fluids with constant viscosities are allowed!"); | ||
| 141 | |||
| 142 | // we know the "upwind factor" is constant, get inner one here and compute derivatives | ||
| 143 |
3/4✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1997 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
|
2000 | static const Scalar up = curElemVolVars[scvf.insideScvIdx()].density() |
| 144 | 6 | / curElemVolVars[scvf.insideScvIdx()].viscosity(); | |
| 145 | |||
| 146 | 2000 | const auto& fluxVarsCache = elemFluxVarsCache[scvf]; | |
| 147 | 2000 | const auto localFaceIdx = fluxVarsCache.ivLocalFaceIndex(); | |
| 148 | 2000 | const auto usesSecondary = fluxVarsCache.usesSecondaryIv(); | |
| 149 | 2000 | const auto switchSign = fluxVarsCache.advectionSwitchFluxSign(); | |
| 150 | |||
| 151 | 2000 | const auto& stencil = fluxVarsCache.advectionStencil(); | |
| 152 | 2000 | const auto& tij = usesSecondary ? fluxVarsCache.advectionSecondaryDataHandle().T()[localFaceIdx] | |
| 153 |
2/4✓ Branch 0 taken 2000 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2000 times.
✗ Branch 3 not taken.
|
4000 | : fluxVarsCache.advectionPrimaryDataHandle().T()[localFaceIdx]; |
| 154 | |||
| 155 | // We assume same the tij are order as the stencil up to stencil.size() | ||
| 156 | // any contribution of Dirichlet BCs is assumed to be placed afterwards | ||
| 157 |
3/6✓ Branch 0 taken 2000 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2000 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2000 times.
✗ Branch 5 not taken.
|
6000 | assert(stencil.size() <= tij.size()); |
| 158 |
4/4✓ Branch 0 taken 7364 times.
✓ Branch 1 taken 2000 times.
✓ Branch 2 taken 7364 times.
✓ Branch 3 taken 2000 times.
|
18728 | for (unsigned int i = 0; i < stencil.size();++i) |
| 159 |
2/2✓ Branch 0 taken 3576 times.
✓ Branch 1 taken 3788 times.
|
14940 | derivativeMatrices[stencil[i]][conti0EqIdx][pressureIdx] += switchSign ? -tij[i]*up |
| 160 | 7576 | : tij[i]*up; | |
| 161 | 2000 | } | |
| 162 | |||
| 163 | //! Flux derivatives for the box scheme | ||
| 164 | template<class JacobianMatrix, class T = TypeTag> | ||
| 165 | std::enable_if_t<GetPropType<T, Properties::GridGeometry>::discMethod == DiscretizationMethods::box, void> | ||
| 166 | 113161 | addFluxDerivatives(JacobianMatrix& A, | |
| 167 | const Problem& problem, | ||
| 168 | const Element& element, | ||
| 169 | const FVElementGeometry& fvGeometry, | ||
| 170 | const ElementVolumeVariables& curElemVolVars, | ||
| 171 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 172 | const SubControlVolumeFace& scvf) const | ||
| 173 | { | ||
| 174 | static_assert(!FluidSystem::isCompressible(0), | ||
| 175 | "1p/incompressiblelocalresidual.hh: Only incompressible fluids are allowed!"); | ||
| 176 | static_assert(FluidSystem::viscosityIsConstant(0), | ||
| 177 | "1p/incompressiblelocalresidual.hh: Only fluids with constant viscosities are allowed!"); | ||
| 178 | |||
| 179 | using AdvectionType = GetPropType<T, Properties::AdvectionType>; | ||
| 180 |
1/4✓ Branch 0 taken 112556 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
113766 | const auto ti = AdvectionType::calculateTransmissibilities(problem, |
| 181 | element, | ||
| 182 | fvGeometry, | ||
| 183 | curElemVolVars, | ||
| 184 | scvf, | ||
| 185 | 226322 | elemFluxVarsCache[scvf]); | |
| 186 | |||
| 187 |
4/4✓ Branch 0 taken 2 times.
✓ Branch 1 taken 113159 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 113159 times.
|
226322 | const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); |
| 188 |
4/4✓ Branch 0 taken 2 times.
✓ Branch 1 taken 113159 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 113151 times.
|
113766 | const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx()); |
| 189 | |||
| 190 | // we know the "upwind factor" is constant, get inner one here and compute derivatives | ||
| 191 |
3/4✓ Branch 0 taken 10 times.
✓ Branch 1 taken 113151 times.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
|
113161 | static const Scalar up = curElemVolVars[scvf.insideScvIdx()].density() |
| 192 | 50 | / curElemVolVars[scvf.insideScvIdx()].viscosity(); | |
| 193 |
4/4✓ Branch 0 taken 445234 times.
✓ Branch 1 taken 113161 times.
✓ Branch 2 taken 444024 times.
✓ Branch 3 taken 112556 times.
|
1115580 | for (const auto& scv : scvs(fvGeometry)) |
| 194 | { | ||
| 195 |
1/2✓ Branch 1 taken 444024 times.
✗ Branch 2 not taken.
|
445234 | auto d = up*ti[scv.indexInElement()]; |
| 196 |
3/6✓ Branch 1 taken 444024 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 444024 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 444024 times.
✗ Branch 8 not taken.
|
890468 | A[insideScv.dofIndex()][scv.dofIndex()][conti0EqIdx][pressureIdx] += d; |
| 197 |
2/4✓ Branch 1 taken 444024 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 444024 times.
✗ Branch 5 not taken.
|
890468 | A[outsideScv.dofIndex()][scv.dofIndex()][conti0EqIdx][pressureIdx] -= d; |
| 198 | } | ||
| 199 | 113161 | } | |
| 200 | |||
| 201 | //! Dirichlet flux derivatives for the cell-centered tpfa scheme | ||
| 202 | template<class PartialDerivativeMatrices, class T = TypeTag> | ||
| 203 | std::enable_if_t<GetPropType<T, Properties::GridGeometry>::discMethod == DiscretizationMethods::cctpfa, void> | ||
| 204 | 2220 | addCCDirichletFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices, | |
| 205 | const Problem& problem, | ||
| 206 | const Element& element, | ||
| 207 | const FVElementGeometry& fvGeometry, | ||
| 208 | const ElementVolumeVariables& curElemVolVars, | ||
| 209 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 210 | const SubControlVolumeFace& scvf) const | ||
| 211 | { | ||
| 212 | // we know the "upwind factor" is constant, get inner one here | ||
| 213 |
3/4✓ Branch 0 taken 11 times.
✓ Branch 1 taken 2209 times.
✓ Branch 3 taken 11 times.
✗ Branch 4 not taken.
|
2221 | static const Scalar up = curElemVolVars[scvf.insideScvIdx()].density() |
| 214 | 32 | / curElemVolVars[scvf.insideScvIdx()].viscosity(); | |
| 215 | 2220 | const auto deriv = elemFluxVarsCache[scvf].advectionTij()*up; | |
| 216 | |||
| 217 | // compute and add partial derivative to the respective given matrices | ||
| 218 | 4440 | derivativeMatrices[scvf.insideScvIdx()][conti0EqIdx][pressureIdx] += deriv; | |
| 219 | 2220 | } | |
| 220 | |||
| 221 | //! Dirichlet flux derivatives for the cell-centered mpfa scheme | ||
| 222 | template<class PartialDerivativeMatrices, class T = TypeTag> | ||
| 223 | std::enable_if_t<GetPropType<T, Properties::GridGeometry>::discMethod == DiscretizationMethods::ccmpfa, void> | ||
| 224 | addCCDirichletFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices, | ||
| 225 | const Problem& problem, | ||
| 226 | const Element& element, | ||
| 227 | const FVElementGeometry& fvGeometry, | ||
| 228 | const ElementVolumeVariables& curElemVolVars, | ||
| 229 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 230 | const SubControlVolumeFace& scvf) const | ||
| 231 | { | ||
| 232 | 112 | addFluxDerivatives(derivativeMatrices, problem, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf); | |
| 233 | } | ||
| 234 | |||
| 235 | //! Robin-type flux derivatives | ||
| 236 | template<class PartialDerivativeMatrices> | ||
| 237 | ✗ | void addRobinFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices, | |
| 238 | const Problem& problem, | ||
| 239 | const Element& element, | ||
| 240 | const FVElementGeometry& fvGeometry, | ||
| 241 | const ElementVolumeVariables& curElemVolVars, | ||
| 242 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 243 | const SubControlVolumeFace& scvf) const | ||
| 244 | { | ||
| 245 | //! Robin-type boundary conditions are problem-specific. | ||
| 246 | //! We can't put a general implementation here - users defining Robin-type BCs | ||
| 247 | //! while using analytical Jacobian assembly must overload this function! | ||
| 248 | ✗ | } | |
| 249 | }; | ||
| 250 | |||
| 251 | } // end namespace Dumux | ||
| 252 | |||
| 253 | #endif | ||
| 254 |