| 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 TwoPModel | ||
| 10 | * \brief Element-wise calculation of the residual and its derivatives | ||
| 11 | * for a two-phase, incompressible test problem. | ||
| 12 | */ | ||
| 13 | |||
| 14 | #ifndef DUMUX_2P_INCOMPRESSIBLE_TEST_LOCAL_RESIDUAL_HH | ||
| 15 | #define DUMUX_2P_INCOMPRESSIBLE_TEST_LOCAL_RESIDUAL_HH | ||
| 16 | |||
| 17 | #include <cmath> | ||
| 18 | |||
| 19 | #include <dumux/common/properties.hh> | ||
| 20 | #include <dumux/common/parameters.hh> | ||
| 21 | |||
| 22 | #include <dumux/discretization/method.hh> | ||
| 23 | #include <dumux/discretization/elementsolution.hh> | ||
| 24 | #include <dumux/discretization/extrusion.hh> | ||
| 25 | |||
| 26 | #include <dumux/porousmediumflow/immiscible/localresidual.hh> | ||
| 27 | #include <dumux/porousmediumflow/2p/formulation.hh> | ||
| 28 | |||
| 29 | namespace Dumux { | ||
| 30 | |||
| 31 | /*! | ||
| 32 | * \ingroup TwoPModel | ||
| 33 | * \brief Element-wise calculation of the residual and its derivatives | ||
| 34 | * for a two-phase, incompressible test problem. | ||
| 35 | */ | ||
| 36 | template<class TypeTag> | ||
| 37 | class TwoPIncompressibleLocalResidual : public ImmiscibleLocalResidual<TypeTag> | ||
| 38 | { | ||
| 39 | using ParentType = ImmiscibleLocalResidual<TypeTag>; | ||
| 40 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 41 | using Problem = GetPropType<TypeTag, Properties::Problem>; | ||
| 42 | using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>; | ||
| 43 | using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView; | ||
| 44 | using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; | ||
| 45 | using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>; | ||
| 46 | using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView; | ||
| 47 | using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>; | ||
| 48 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 49 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 50 | using SubControlVolume = typename GridGeometry::SubControlVolume; | ||
| 51 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 52 | using Extrusion = Extrusion_t<GridGeometry>; | ||
| 53 | using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView; | ||
| 54 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 55 | using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>; | ||
| 56 | using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>; | ||
| 57 | |||
| 58 | static constexpr int numPhases = ModelTraits::numFluidPhases(); | ||
| 59 | static constexpr int pressureIdx = ModelTraits::Indices::pressureIdx; | ||
| 60 | static constexpr int saturationIdx = ModelTraits::Indices::saturationIdx; | ||
| 61 | static constexpr int conti0EqIdx = ModelTraits::Indices::conti0EqIdx; | ||
| 62 | |||
| 63 | public: | ||
| 64 |
2/4✓ Branch 1 taken 3715381 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3715381 times.
✗ Branch 5 not taken.
|
7430762 | using ParentType::ParentType; |
| 65 | |||
| 66 | /*! | ||
| 67 | * \brief Adds storage derivatives for wetting and nonwetting phase | ||
| 68 | * | ||
| 69 | * Compute storage derivatives for the wetting and the nonwetting phase with respect to \f$p_w\f$ | ||
| 70 | * and \f$S_n\f$. | ||
| 71 | * | ||
| 72 | * \param partialDerivatives The partial derivatives | ||
| 73 | * \param problem The problem | ||
| 74 | * \param element The element | ||
| 75 | * \param fvGeometry The finite volume element geometry | ||
| 76 | * \param curVolVars The current volume variables | ||
| 77 | * \param scv The sub control volume | ||
| 78 | */ | ||
| 79 | template<class PartialDerivativeMatrix> | ||
| 80 | ✗ | void addStorageDerivatives(PartialDerivativeMatrix& partialDerivatives, | |
| 81 | const Problem& problem, | ||
| 82 | const Element& element, | ||
| 83 | const FVElementGeometry& fvGeometry, | ||
| 84 | const VolumeVariables& curVolVars, | ||
| 85 | const SubControlVolume& scv) const | ||
| 86 | { | ||
| 87 | static_assert(!FluidSystem::isCompressible(0), | ||
| 88 | "2p/incompressiblelocalresidual.hh: Only incompressible fluids are allowed!"); | ||
| 89 | static_assert(!FluidSystem::isCompressible(1), | ||
| 90 | "2p/incompressiblelocalresidual.hh: Only incompressible fluids are allowed!"); | ||
| 91 | static_assert(ModelTraits::numFluidPhases() == 2, | ||
| 92 | "2p/incompressiblelocalresidual.hh: Only two-phase models are allowed!"); | ||
| 93 | static_assert(ModelTraits::priVarFormulation() == TwoPFormulation::p0s1, | ||
| 94 | "2p/incompressiblelocalresidual.hh: Analytic differentiation has to be checked for p1-s0 formulation!"); | ||
| 95 | |||
| 96 | ✗ | const auto poreVolume = Extrusion::volume(fvGeometry, scv)*curVolVars.porosity(); | |
| 97 | ✗ | const auto dt = this->timeLoop().timeStepSize(); | |
| 98 | |||
| 99 | // partial derivative of the phase storage terms w.r.t. S_n | ||
| 100 | ✗ | partialDerivatives[conti0EqIdx+0][saturationIdx] -= poreVolume*curVolVars.density(0)/dt; | |
| 101 | ✗ | partialDerivatives[conti0EqIdx+1][saturationIdx] += poreVolume*curVolVars.density(1)/dt; | |
| 102 | ✗ | } | |
| 103 | |||
| 104 | /*! | ||
| 105 | * \brief Adds source derivatives for wetting and nonwetting phase. | ||
| 106 | * | ||
| 107 | * \param partialDerivatives The partial derivatives | ||
| 108 | * \param problem The problem | ||
| 109 | * \param element The element | ||
| 110 | * \param fvGeometry The finite volume element geometry | ||
| 111 | * \param curVolVars The current volume variables | ||
| 112 | * \param scv The sub control volume | ||
| 113 | */ | ||
| 114 | template<class PartialDerivativeMatrix> | ||
| 115 | ✗ | void addSourceDerivatives(PartialDerivativeMatrix& partialDerivatives, | |
| 116 | const Problem& problem, | ||
| 117 | const Element& element, | ||
| 118 | const FVElementGeometry& fvGeometry, | ||
| 119 | const VolumeVariables& curVolVars, | ||
| 120 | const SubControlVolume& scv) const | ||
| 121 | ✗ | { /* TODO maybe forward to problem for the user to implement the source derivatives?*/ } | |
| 122 | |||
| 123 | /*! | ||
| 124 | * \brief Adds flux derivatives for wetting and nonwetting phase for cell-centered FVM using TPFA | ||
| 125 | * | ||
| 126 | * Compute derivatives for the wetting and the nonwetting phase flux with respect to \f$p_w\f$ | ||
| 127 | * and \f$S_n\f$. | ||
| 128 | * | ||
| 129 | * \param derivativeMatrices The partial derivatives | ||
| 130 | * \param problem The problem | ||
| 131 | * \param element The element | ||
| 132 | * \param fvGeometry The finite volume element geometry | ||
| 133 | * \param curElemVolVars The current element volume variables | ||
| 134 | * \param elemFluxVarsCache The element flux variables cache | ||
| 135 | * \param scvf The sub control volume face | ||
| 136 | */ | ||
| 137 | template<class PartialDerivativeMatrices, class T = TypeTag> | ||
| 138 | std::enable_if_t<GetPropType<T, Properties::GridGeometry>::discMethod == DiscretizationMethods::cctpfa, void> | ||
| 139 | 317152 | addFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices, | |
| 140 | const Problem& problem, | ||
| 141 | const Element& element, | ||
| 142 | const FVElementGeometry& fvGeometry, | ||
| 143 | const ElementVolumeVariables& curElemVolVars, | ||
| 144 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 145 | const SubControlVolumeFace& scvf) const | ||
| 146 | { | ||
| 147 | static_assert(!FluidSystem::isCompressible(0), | ||
| 148 | "2p/incompressiblelocalresidual.hh: Only incompressible fluids are allowed!"); | ||
| 149 | static_assert(!FluidSystem::isCompressible(1), | ||
| 150 | "2p/incompressiblelocalresidual.hh: Only incompressible fluids are allowed!"); | ||
| 151 | static_assert(FluidSystem::viscosityIsConstant(0), | ||
| 152 | "2p/incompressiblelocalresidual.hh: Only fluids with constant viscosities are allowed!"); | ||
| 153 | static_assert(FluidSystem::viscosityIsConstant(1), | ||
| 154 | "2p/incompressiblelocalresidual.hh: Only fluids with constant viscosities are allowed!"); | ||
| 155 | static_assert(ModelTraits::numFluidPhases() == 2, | ||
| 156 | "2p/incompressiblelocalresidual.hh: Only two-phase models are allowed!"); | ||
| 157 | static_assert(ModelTraits::priVarFormulation() == TwoPFormulation::p0s1, | ||
| 158 | "2p/incompressiblelocalresidual.hh: Analytic differentiation has to be checked for p1-s0 formulation!"); | ||
| 159 | |||
| 160 | using AdvectionType = GetPropType<TypeTag, Properties::AdvectionType>; | ||
| 161 | |||
| 162 | // evaluate the current wetting phase Darcy flux and resulting upwind weights | ||
| 163 |
5/8✓ Branch 0 taken 1 times.
✓ Branch 1 taken 317151 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
|
317152 | static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight"); |
| 164 | 317152 | const auto flux_w = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache); | |
| 165 | 317152 | const auto flux_n = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 1, elemFluxVarsCache); | |
| 166 |
4/4✓ Branch 0 taken 156880 times.
✓ Branch 1 taken 160272 times.
✓ Branch 2 taken 156880 times.
✓ Branch 3 taken 160272 times.
|
634304 | const auto insideWeight_w = std::signbit(flux_w) ? (1.0 - upwindWeight) : upwindWeight; |
| 167 | 317152 | const auto outsideWeight_w = 1.0 - insideWeight_w; | |
| 168 |
4/4✓ Branch 0 taken 157072 times.
✓ Branch 1 taken 160080 times.
✓ Branch 2 taken 157072 times.
✓ Branch 3 taken 160080 times.
|
634304 | const auto insideWeight_n = std::signbit(flux_n) ? (1.0 - upwindWeight) : upwindWeight; |
| 169 | 317152 | const auto outsideWeight_n = 1.0 - insideWeight_n; | |
| 170 | |||
| 171 | // get references to the two participating vol vars & parameters | ||
| 172 | 317152 | const auto insideScvIdx = scvf.insideScvIdx(); | |
| 173 | 317152 | const auto outsideScvIdx = scvf.outsideScvIdx(); | |
| 174 | 317152 | const auto outsideElement = fvGeometry.gridGeometry().element(outsideScvIdx); | |
| 175 |
1/2✓ Branch 0 taken 317152 times.
✗ Branch 1 not taken.
|
317152 | const auto& insideScv = fvGeometry.scv(insideScvIdx); |
| 176 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 317152 times.
|
317152 | const auto& outsideScv = fvGeometry.scv(outsideScvIdx); |
| 177 | 317152 | const auto& insideVolVars = curElemVolVars[insideScvIdx]; | |
| 178 | 317152 | const auto& outsideVolVars = curElemVolVars[outsideScvIdx]; | |
| 179 | |||
| 180 | 317152 | const auto insidefluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction( | |
| 181 | 951456 | element, insideScv, elementSolution<FVElementGeometry>(insideVolVars.priVars()) | |
| 182 | ); | ||
| 183 | 317152 | const auto outsidefluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction( | |
| 184 | 951456 | outsideElement, outsideScv, elementSolution<FVElementGeometry>(outsideVolVars.priVars()) | |
| 185 | ); | ||
| 186 | |||
| 187 | // get references to the two participating derivative matrices | ||
| 188 | 317152 | auto& dI_dI = derivativeMatrices[insideScvIdx]; | |
| 189 | 317152 | auto& dI_dJ = derivativeMatrices[outsideScvIdx]; | |
| 190 | |||
| 191 | // some quantities to be reused (rho & mu are constant and thus equal for all cells) | ||
| 192 |
3/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 317151 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
317152 | static const auto rho_w = insideVolVars.density(0); |
| 193 |
3/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 317151 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
317152 | static const auto rho_n = insideVolVars.density(1); |
| 194 |
3/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 317151 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
317152 | static const auto rhow_muw = rho_w/insideVolVars.viscosity(0); |
| 195 |
3/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 317151 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
317152 | static const auto rhon_mun = rho_n/insideVolVars.viscosity(1); |
| 196 | 317152 | const auto rhowKrw_muw_inside = rho_w*insideVolVars.mobility(0); | |
| 197 | 317152 | const auto rhonKrn_mun_inside = rho_n*insideVolVars.mobility(1); | |
| 198 | 317152 | const auto rhowKrw_muw_outside = rho_w*outsideVolVars.mobility(0); | |
| 199 | 317152 | const auto rhonKrn_mun_outside = rho_n*outsideVolVars.mobility(1); | |
| 200 | |||
| 201 | // derivative w.r.t. to Sn is the negative of the one w.r.t. Sw | ||
| 202 | 317152 | const auto insideSw = insideVolVars.saturation(0); | |
| 203 | 317152 | const auto outsideSw = outsideVolVars.saturation(0); | |
| 204 | 317152 | const auto dKrw_dSn_inside = -1.0*insidefluidMatrixInteraction.dkrw_dsw(insideSw); | |
| 205 | 317152 | const auto dKrw_dSn_outside = -1.0*outsidefluidMatrixInteraction.dkrw_dsw(outsideSw); | |
| 206 | 317152 | const auto dKrn_dSn_inside = -1.0*insidefluidMatrixInteraction.dkrn_dsw(insideSw); | |
| 207 | 317152 | const auto dKrn_dSn_outside = -1.0*outsidefluidMatrixInteraction.dkrn_dsw(outsideSw); | |
| 208 | 317152 | const auto dpc_dSn_inside = -1.0*insidefluidMatrixInteraction.dpc_dsw(insideSw); | |
| 209 | 317152 | const auto dpc_dSn_outside = -1.0*outsidefluidMatrixInteraction.dpc_dsw(outsideSw); | |
| 210 | |||
| 211 | 317152 | const auto tij = elemFluxVarsCache[scvf].advectionTij(); | |
| 212 | |||
| 213 | // precalculate values | ||
| 214 | 317152 | const auto up_w = rhowKrw_muw_inside*insideWeight_w + rhowKrw_muw_outside*outsideWeight_w; | |
| 215 | 317152 | const auto up_n = rhonKrn_mun_inside*insideWeight_n + rhonKrn_mun_outside*outsideWeight_n; | |
| 216 | 317152 | const auto rho_mu_flux_w = rhow_muw*flux_w; | |
| 217 | 317152 | const auto rho_mu_flux_n = rhon_mun*flux_n; | |
| 218 | 317152 | const auto tij_up_w = tij*up_w; | |
| 219 | 317152 | const auto tij_up_n = tij*up_n; | |
| 220 | |||
| 221 | // partial derivative of the wetting phase flux w.r.t. p_w | ||
| 222 | 634304 | dI_dI[conti0EqIdx+0][pressureIdx] += tij_up_w; | |
| 223 | 634304 | dI_dJ[conti0EqIdx+0][pressureIdx] -= tij_up_w; | |
| 224 | |||
| 225 | // partial derivative of the wetting phase flux w.r.t. S_n | ||
| 226 | 634304 | dI_dI[conti0EqIdx+0][saturationIdx] += rho_mu_flux_w*dKrw_dSn_inside*insideWeight_w; | |
| 227 | 634304 | dI_dJ[conti0EqIdx+0][saturationIdx] += rho_mu_flux_w*dKrw_dSn_outside*outsideWeight_w; | |
| 228 | |||
| 229 | // partial derivative of the nonwetting phase flux w.r.t. p_w | ||
| 230 | 634304 | dI_dI[conti0EqIdx+1][pressureIdx] += tij_up_n; | |
| 231 | 634304 | dI_dJ[conti0EqIdx+1][pressureIdx] -= tij_up_n; | |
| 232 | |||
| 233 | // partial derivative of the nonwetting phase flux w.r.t. S_n (relative permeability derivative contribution) | ||
| 234 | 634304 | dI_dI[conti0EqIdx+1][saturationIdx] += rho_mu_flux_n*dKrn_dSn_inside*insideWeight_n; | |
| 235 | 634304 | dI_dJ[conti0EqIdx+1][saturationIdx] += rho_mu_flux_n*dKrn_dSn_outside*outsideWeight_n; | |
| 236 | |||
| 237 | // partial derivative of the nonwetting phase flux w.r.t. S_n (capillary pressure derivative contribution) | ||
| 238 | 634304 | dI_dI[conti0EqIdx+1][saturationIdx] += tij_up_n*dpc_dSn_inside; | |
| 239 | 634304 | dI_dJ[conti0EqIdx+1][saturationIdx] -= tij_up_n*dpc_dSn_outside; | |
| 240 | 317152 | } | |
| 241 | |||
| 242 | /*! | ||
| 243 | * \brief Adds flux derivatives for wetting and nonwetting phase for box method | ||
| 244 | * | ||
| 245 | * Compute derivatives for the wetting and the nonwetting phase flux with respect to \f$p_w\f$ | ||
| 246 | * and \f$S_n\f$. | ||
| 247 | * | ||
| 248 | * \param A The Jacobian Matrix | ||
| 249 | * \param problem The problem | ||
| 250 | * \param element The element | ||
| 251 | * \param fvGeometry The finite volume element geometry | ||
| 252 | * \param curElemVolVars The current element volume variables | ||
| 253 | * \param elemFluxVarsCache The element flux variables cache | ||
| 254 | * \param scvf The sub control volume face | ||
| 255 | */ | ||
| 256 | template<class JacobianMatrix, class T = TypeTag> | ||
| 257 | std::enable_if_t<GetPropType<T, Properties::GridGeometry>::discMethod == DiscretizationMethods::box, void> | ||
| 258 | 276480 | addFluxDerivatives(JacobianMatrix& A, | |
| 259 | const Problem& problem, | ||
| 260 | const Element& element, | ||
| 261 | const FVElementGeometry& fvGeometry, | ||
| 262 | const ElementVolumeVariables& curElemVolVars, | ||
| 263 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 264 | const SubControlVolumeFace& scvf) const | ||
| 265 | { | ||
| 266 | static_assert(!FluidSystem::isCompressible(0), | ||
| 267 | "2p/incompressiblelocalresidual.hh: Only incompressible fluids are allowed!"); | ||
| 268 | static_assert(!FluidSystem::isCompressible(1), | ||
| 269 | "2p/incompressiblelocalresidual.hh: Only incompressible fluids are allowed!"); | ||
| 270 | static_assert(FluidSystem::viscosityIsConstant(0), | ||
| 271 | "2p/incompressiblelocalresidual.hh: Only fluids with constant viscosities are allowed!"); | ||
| 272 | static_assert(FluidSystem::viscosityIsConstant(1), | ||
| 273 | "2p/incompressiblelocalresidual.hh: Only fluids with constant viscosities are allowed!"); | ||
| 274 | static_assert(ModelTraits::numFluidPhases() == 2, | ||
| 275 | "2p/incompressiblelocalresidual.hh: Only two-phase models are allowed!"); | ||
| 276 | static_assert(ModelTraits::priVarFormulation() == TwoPFormulation::p0s1, | ||
| 277 | "2p/incompressiblelocalresidual.hh: Analytic differentiation has to be checked for p0-s1 formulation!"); | ||
| 278 | |||
| 279 | using AdvectionType = GetPropType<TypeTag, Properties::AdvectionType>; | ||
| 280 | |||
| 281 | // evaluate the current wetting phase Darcy flux and resulting upwind weights | ||
| 282 |
5/8✓ Branch 0 taken 1 times.
✓ Branch 1 taken 276479 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
|
276480 | static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight"); |
| 283 | 276480 | const auto flux_w = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache); | |
| 284 | 276480 | const auto flux_n = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 1, elemFluxVarsCache); | |
| 285 |
4/4✓ Branch 0 taken 14247 times.
✓ Branch 1 taken 262233 times.
✓ Branch 2 taken 14247 times.
✓ Branch 3 taken 262233 times.
|
552960 | const auto insideWeight_w = std::signbit(flux_w) ? (1.0 - upwindWeight) : upwindWeight; |
| 286 | 276480 | const auto outsideWeight_w = 1.0 - insideWeight_w; | |
| 287 |
4/4✓ Branch 0 taken 150008 times.
✓ Branch 1 taken 126472 times.
✓ Branch 2 taken 150008 times.
✓ Branch 3 taken 126472 times.
|
552960 | const auto insideWeight_n = std::signbit(flux_n) ? (1.0 - upwindWeight) : upwindWeight; |
| 288 | 276480 | const auto outsideWeight_n = 1.0 - insideWeight_n; | |
| 289 | |||
| 290 | // get references to the two participating vol vars & parameters | ||
| 291 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 276480 times.
|
276480 | const auto insideScvIdx = scvf.insideScvIdx(); |
| 292 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 276480 times.
|
276480 | const auto outsideScvIdx = scvf.outsideScvIdx(); |
| 293 | 276480 | const auto& insideScv = fvGeometry.scv(insideScvIdx); | |
| 294 | 276480 | const auto& outsideScv = fvGeometry.scv(outsideScvIdx); | |
| 295 | 276480 | const auto& insideVolVars = curElemVolVars[insideScv]; | |
| 296 | 276480 | const auto& outsideVolVars = curElemVolVars[outsideScv]; | |
| 297 | |||
| 298 | 276480 | const auto elemSol = elementSolution(element, curElemVolVars, fvGeometry); | |
| 299 | |||
| 300 | 552960 | const auto insidefluidMatrixInteraction | |
| 301 | 552960 | = problem.spatialParams().fluidMatrixInteraction(element, insideScv, elemSol); | |
| 302 | 552960 | const auto outsidefluidMatrixInteraction | |
| 303 | 552960 | = problem.spatialParams().fluidMatrixInteraction(element, outsideScv, elemSol); | |
| 304 | |||
| 305 | // some quantities to be reused (rho & mu are constant and thus equal for all cells) | ||
| 306 |
3/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 276479 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
276480 | static const auto rho_w = insideVolVars.density(0); |
| 307 |
3/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 276479 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
276480 | static const auto rho_n = insideVolVars.density(1); |
| 308 |
3/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 276479 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
276480 | static const auto rhow_muw = rho_w/insideVolVars.viscosity(0); |
| 309 |
3/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 276479 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
276480 | static const auto rhon_mun = rho_n/insideVolVars.viscosity(1); |
| 310 | 276480 | const auto rhowKrw_muw_inside = rho_w*insideVolVars.mobility(0); | |
| 311 | 276480 | const auto rhonKrn_mun_inside = rho_n*insideVolVars.mobility(1); | |
| 312 | 276480 | const auto rhowKrw_muw_outside = rho_w*outsideVolVars.mobility(0); | |
| 313 | 276480 | const auto rhonKrn_mun_outside = rho_n*outsideVolVars.mobility(1); | |
| 314 | |||
| 315 | // let the Law for the advective fluxes calculate the transmissibilities | ||
| 316 |
1/4✓ Branch 0 taken 276480 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
276480 | const auto ti = AdvectionType::calculateTransmissibilities(problem, |
| 317 | element, | ||
| 318 | fvGeometry, | ||
| 319 | curElemVolVars, | ||
| 320 | scvf, | ||
| 321 | 552960 | elemFluxVarsCache[scvf]); | |
| 322 | |||
| 323 | // get the rows of the jacobian matrix for the inside/outside scv | ||
| 324 | 276480 | auto& dI_dJ_inside = A[insideScv.dofIndex()]; | |
| 325 | 276480 | auto& dI_dJ_outside = A[outsideScv.dofIndex()]; | |
| 326 | |||
| 327 | // precalculate values | ||
| 328 | 276480 | const auto up_w = rhowKrw_muw_inside*insideWeight_w + rhowKrw_muw_outside*outsideWeight_w; | |
| 329 | 276480 | const auto up_n = rhonKrn_mun_inside*insideWeight_n + rhonKrn_mun_outside*outsideWeight_n; | |
| 330 | 276480 | const auto rho_mu_flux_w = rhow_muw*flux_w; | |
| 331 | 276480 | const auto rho_mu_flux_n = rhon_mun*flux_n; | |
| 332 | |||
| 333 | // add the partial derivatives w.r.t all scvs in the element | ||
| 334 |
4/4✓ Branch 0 taken 1105920 times.
✓ Branch 1 taken 276480 times.
✓ Branch 2 taken 1105920 times.
✓ Branch 3 taken 276480 times.
|
1658880 | for (const auto& scvJ : scvs(fvGeometry)) |
| 335 | { | ||
| 336 | 1105920 | const auto globalJ = scvJ.dofIndex(); | |
| 337 | 1105920 | const auto localJ = scvJ.indexInElement(); | |
| 338 | |||
| 339 | // the transmissibily associated with the scvJ | ||
| 340 |
1/2✓ Branch 1 taken 1105920 times.
✗ Branch 2 not taken.
|
1105920 | const auto tj = ti[localJ]; |
| 341 | |||
| 342 | // partial derivative of the wetting phase flux w.r.t. p_w | ||
| 343 | 1105920 | const auto tj_up_w = tj*up_w; | |
| 344 |
3/6✓ Branch 1 taken 1105920 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1105920 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1105920 times.
✗ Branch 8 not taken.
|
1105920 | dI_dJ_inside[globalJ][conti0EqIdx+0][pressureIdx] += tj_up_w; |
| 345 |
3/6✓ Branch 1 taken 1105920 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1105920 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1105920 times.
✗ Branch 8 not taken.
|
1105920 | dI_dJ_outside[globalJ][conti0EqIdx+0][pressureIdx] -= tj_up_w; |
| 346 | |||
| 347 | // partial derivative of the nonwetting phase flux w.r.t. p_w | ||
| 348 | 1105920 | const auto tj_up_n = tj*up_n; | |
| 349 |
3/6✓ Branch 1 taken 1105920 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1105920 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1105920 times.
✗ Branch 8 not taken.
|
1105920 | dI_dJ_inside[globalJ][conti0EqIdx+1][pressureIdx] += tj_up_n; |
| 350 |
5/6✓ Branch 1 taken 1105920 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 276480 times.
✓ Branch 4 taken 829440 times.
✓ Branch 5 taken 276480 times.
✓ Branch 6 taken 829440 times.
|
1105920 | dI_dJ_outside[globalJ][conti0EqIdx+1][pressureIdx] -= tj_up_n; |
| 351 | |||
| 352 | // partial derivatives w.r.t. S_n (are the negative of those w.r.t sw) | ||
| 353 | // relative permeability contributions only for inside/outside | ||
| 354 |
2/2✓ Branch 0 taken 276480 times.
✓ Branch 1 taken 829440 times.
|
1105920 | if (localJ == insideScvIdx) |
| 355 | { | ||
| 356 | // partial derivative of the wetting phase flux w.r.t. S_n | ||
| 357 | 276480 | const auto insideSw = insideVolVars.saturation(0); | |
| 358 | 276480 | const auto dKrw_dSn_inside = -1.0*insidefluidMatrixInteraction.dkrw_dsw(insideSw); | |
| 359 | 276480 | const auto dFluxW_dSnJ = rho_mu_flux_w*dKrw_dSn_inside*insideWeight_w; | |
| 360 |
3/6✓ Branch 1 taken 276480 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 276480 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 276480 times.
✗ Branch 8 not taken.
|
276480 | dI_dJ_inside[globalJ][conti0EqIdx+0][saturationIdx] += dFluxW_dSnJ; |
| 361 |
1/2✓ Branch 1 taken 276480 times.
✗ Branch 2 not taken.
|
276480 | dI_dJ_outside[globalJ][conti0EqIdx+0][saturationIdx] -= dFluxW_dSnJ; |
| 362 | |||
| 363 | // partial derivative of the nonwetting phase flux w.r.t. S_n (k_rn contribution) | ||
| 364 | 276480 | const auto dKrn_dSn_inside = -1.0*insidefluidMatrixInteraction.dkrn_dsw(insideSw); | |
| 365 | 276480 | const auto dFluxN_dSnJ_krn = rho_mu_flux_n*dKrn_dSn_inside*insideWeight_n; | |
| 366 |
3/6✓ Branch 1 taken 276480 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 276480 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 276480 times.
✗ Branch 8 not taken.
|
276480 | dI_dJ_inside[globalJ][conti0EqIdx+1][saturationIdx] += dFluxN_dSnJ_krn; |
| 367 |
1/2✓ Branch 1 taken 276480 times.
✗ Branch 2 not taken.
|
276480 | dI_dJ_outside[globalJ][conti0EqIdx+1][saturationIdx] -= dFluxN_dSnJ_krn; |
| 368 | |||
| 369 | // partial derivative of the nonwetting phase flux w.r.t. S_n (p_c contribution) | ||
| 370 | 276480 | const auto dFluxN_dSnJ_pc = -1.0*tj_up_n*insidefluidMatrixInteraction.dpc_dsw(insideSw); | |
| 371 |
3/6✓ Branch 1 taken 276480 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 276480 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 276480 times.
✗ Branch 8 not taken.
|
276480 | dI_dJ_inside[globalJ][conti0EqIdx+1][saturationIdx] += dFluxN_dSnJ_pc; |
| 372 |
1/2✓ Branch 1 taken 276480 times.
✗ Branch 2 not taken.
|
276480 | dI_dJ_outside[globalJ][conti0EqIdx+1][saturationIdx] -= dFluxN_dSnJ_pc; |
| 373 | } | ||
| 374 |
2/2✓ Branch 0 taken 276480 times.
✓ Branch 1 taken 552960 times.
|
829440 | else if (localJ == outsideScvIdx) |
| 375 | { | ||
| 376 | // see comments for (globalJ == insideScvIdx) | ||
| 377 | 276480 | const auto outsideSw = outsideVolVars.saturation(0); | |
| 378 | 276480 | const auto dKrw_dSn_outside = -1.0*outsidefluidMatrixInteraction.dkrw_dsw(outsideSw); | |
| 379 | 276480 | const auto dFluxW_dSnJ = rho_mu_flux_w*dKrw_dSn_outside*outsideWeight_w; | |
| 380 |
3/6✓ Branch 1 taken 276480 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 276480 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 276480 times.
✗ Branch 8 not taken.
|
276480 | dI_dJ_inside[globalJ][conti0EqIdx+0][saturationIdx] += dFluxW_dSnJ; |
| 381 |
1/2✓ Branch 1 taken 276480 times.
✗ Branch 2 not taken.
|
276480 | dI_dJ_outside[globalJ][conti0EqIdx+0][saturationIdx] -= dFluxW_dSnJ; |
| 382 | |||
| 383 | 276480 | const auto dKrn_dSn_outside = -1.0*outsidefluidMatrixInteraction.dkrn_dsw(outsideSw); | |
| 384 | 276480 | const auto dFluxN_dSnJ_krn = rho_mu_flux_n*dKrn_dSn_outside*outsideWeight_n; | |
| 385 |
3/6✓ Branch 1 taken 276480 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 276480 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 276480 times.
✗ Branch 8 not taken.
|
276480 | dI_dJ_inside[globalJ][conti0EqIdx+1][saturationIdx] += dFluxN_dSnJ_krn; |
| 386 |
1/2✓ Branch 1 taken 276480 times.
✗ Branch 2 not taken.
|
276480 | dI_dJ_outside[globalJ][conti0EqIdx+1][saturationIdx] -= dFluxN_dSnJ_krn; |
| 387 | |||
| 388 | 276480 | const auto dFluxN_dSnJ_pc = -1.0*tj_up_n*outsidefluidMatrixInteraction.dpc_dsw(outsideSw); | |
| 389 |
3/6✓ Branch 1 taken 276480 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 276480 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 276480 times.
✗ Branch 8 not taken.
|
276480 | dI_dJ_inside[globalJ][conti0EqIdx+1][saturationIdx] += dFluxN_dSnJ_pc; |
| 390 |
1/2✓ Branch 1 taken 276480 times.
✗ Branch 2 not taken.
|
276480 | dI_dJ_outside[globalJ][conti0EqIdx+1][saturationIdx] -= dFluxN_dSnJ_pc; |
| 391 | } | ||
| 392 | else | ||
| 393 | { | ||
| 394 | 1105920 | const auto& fluidMatrixInteractionJ | |
| 395 | 1105920 | = problem.spatialParams().fluidMatrixInteraction(element, scvJ, elemSol); | |
| 396 | |||
| 397 | 1105920 | const auto swJ = curElemVolVars[scvJ].saturation(0); | |
| 398 | 552960 | const auto dFluxN_dSnJ_pc = tj_up_n*fluidMatrixInteractionJ.dpc_dsw(swJ); | |
| 399 |
3/6✓ Branch 1 taken 552960 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 552960 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 552960 times.
✗ Branch 8 not taken.
|
552960 | dI_dJ_inside[globalJ][conti0EqIdx+1][saturationIdx] -= dFluxN_dSnJ_pc; |
| 400 |
1/2✓ Branch 1 taken 552960 times.
✗ Branch 2 not taken.
|
552960 | dI_dJ_outside[globalJ][conti0EqIdx+1][saturationIdx] += dFluxN_dSnJ_pc; |
| 401 | } | ||
| 402 | } | ||
| 403 | 276480 | } | |
| 404 | |||
| 405 | /*! | ||
| 406 | * \brief Adds cell-centered Dirichlet flux derivatives for wetting and nonwetting phase | ||
| 407 | * | ||
| 408 | * Compute derivatives for the wetting and the nonwetting phase flux with respect to \f$p_w\f$ | ||
| 409 | * and \f$S_n\f$. | ||
| 410 | * | ||
| 411 | * \param derivativeMatrices The matrices containing the derivatives | ||
| 412 | * \param problem The problem | ||
| 413 | * \param element The element | ||
| 414 | * \param fvGeometry The finite volume element geometry | ||
| 415 | * \param curElemVolVars The current element volume variables | ||
| 416 | * \param elemFluxVarsCache The element flux variables cache | ||
| 417 | * \param scvf The sub control volume face | ||
| 418 | */ | ||
| 419 | template<class PartialDerivativeMatrices> | ||
| 420 | 3392 | void addCCDirichletFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices, | |
| 421 | const Problem& problem, | ||
| 422 | const Element& element, | ||
| 423 | const FVElementGeometry& fvGeometry, | ||
| 424 | const ElementVolumeVariables& curElemVolVars, | ||
| 425 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 426 | const SubControlVolumeFace& scvf) const | ||
| 427 | { | ||
| 428 | using AdvectionType = GetPropType<TypeTag, Properties::AdvectionType>; | ||
| 429 | |||
| 430 | // evaluate the current wetting phase Darcy flux and resulting upwind weights | ||
| 431 |
5/8✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3391 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
|
3392 | static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight"); |
| 432 | 3392 | const auto flux_w = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache); | |
| 433 | 3392 | const auto flux_n = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 1, elemFluxVarsCache); | |
| 434 |
4/4✓ Branch 0 taken 1124 times.
✓ Branch 1 taken 2268 times.
✓ Branch 2 taken 1124 times.
✓ Branch 3 taken 2268 times.
|
6784 | const auto insideWeight_w = std::signbit(flux_w) ? (1.0 - upwindWeight) : upwindWeight; |
| 435 | 3392 | const auto outsideWeight_w = 1.0 - insideWeight_w; | |
| 436 |
4/4✓ Branch 0 taken 1124 times.
✓ Branch 1 taken 2268 times.
✓ Branch 2 taken 1124 times.
✓ Branch 3 taken 2268 times.
|
6784 | const auto insideWeight_n = std::signbit(flux_n) ? (1.0 - upwindWeight) : upwindWeight; |
| 437 | 3392 | const auto outsideWeight_n = 1.0 - insideWeight_n; | |
| 438 | |||
| 439 | // get references to the two participating vol vars & parameters | ||
| 440 |
1/2✓ Branch 0 taken 3392 times.
✗ Branch 1 not taken.
|
3392 | const auto insideScvIdx = scvf.insideScvIdx(); |
| 441 |
1/2✓ Branch 0 taken 3392 times.
✗ Branch 1 not taken.
|
3392 | const auto& insideScv = fvGeometry.scv(insideScvIdx); |
| 442 | 3392 | const auto& insideVolVars = curElemVolVars[insideScvIdx]; | |
| 443 | 6784 | const auto& outsideVolVars = curElemVolVars[scvf.outsideScvIdx()]; | |
| 444 | |||
| 445 | 3392 | const auto insidefluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction( | |
| 446 | 10176 | element, insideScv, elementSolution<FVElementGeometry>(insideVolVars.priVars()) | |
| 447 | ); | ||
| 448 | |||
| 449 | // some quantities to be reused (rho & mu are constant and thus equal for all cells) | ||
| 450 |
3/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3391 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
3392 | static const auto rho_w = insideVolVars.density(0); |
| 451 |
3/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3391 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
3392 | static const auto rho_n = insideVolVars.density(1); |
| 452 |
3/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3391 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
3392 | static const auto rhow_muw = rho_w/insideVolVars.viscosity(0); |
| 453 |
3/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3391 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
3392 | static const auto rhon_mun = rho_n/insideVolVars.viscosity(1); |
| 454 | 3392 | const auto rhowKrw_muw_inside = rho_w*insideVolVars.mobility(0); | |
| 455 | 3392 | const auto rhonKrn_mun_inside = rho_n*insideVolVars.mobility(1); | |
| 456 | 3392 | const auto rhowKrw_muw_outside = rho_w*outsideVolVars.mobility(0); | |
| 457 | 3392 | const auto rhonKrn_mun_outside = rho_n*outsideVolVars.mobility(1); | |
| 458 | |||
| 459 | // get reference to the inside derivative matrix | ||
| 460 | 3392 | auto& dI_dI = derivativeMatrices[insideScvIdx]; | |
| 461 | |||
| 462 | // derivative w.r.t. to Sn is the negative of the one w.r.t. Sw | ||
| 463 | 3392 | const auto insideSw = insideVolVars.saturation(0); | |
| 464 | 3392 | const auto dKrw_dSn_inside = -1.0*insidefluidMatrixInteraction.dkrw_dsw(insideSw); | |
| 465 | 3392 | const auto dKrn_dSn_inside = -1.0*insidefluidMatrixInteraction.dkrn_dsw(insideSw); | |
| 466 | 3392 | const auto dpc_dSn_inside = -1.0*insidefluidMatrixInteraction.dpc_dsw(insideSw); | |
| 467 | |||
| 468 | 3392 | const auto tij = elemFluxVarsCache[scvf].advectionTij(); | |
| 469 | // partial derivative of the wetting phase flux w.r.t. p_w | ||
| 470 | 3392 | const auto up_w = rhowKrw_muw_inside*insideWeight_w + rhowKrw_muw_outside*outsideWeight_w; | |
| 471 | 6784 | dI_dI[conti0EqIdx+0][pressureIdx] += tij*up_w; | |
| 472 | |||
| 473 | // partial derivative of the wetting phase flux w.r.t. S_n | ||
| 474 | 6784 | dI_dI[conti0EqIdx+0][saturationIdx] += rhow_muw*flux_w*dKrw_dSn_inside*insideWeight_w; | |
| 475 | |||
| 476 | // partial derivative of the nonwetting phase flux w.r.t. p_w | ||
| 477 | 3392 | const auto up_n = rhonKrn_mun_inside*insideWeight_n + rhonKrn_mun_outside*outsideWeight_n; | |
| 478 | 6784 | dI_dI[conti0EqIdx+1][pressureIdx] += tij*up_n; | |
| 479 | |||
| 480 | // partial derivative of the nonwetting phase flux w.r.t. S_n (relative permeability derivative contribution) | ||
| 481 | 6784 | dI_dI[conti0EqIdx+1][saturationIdx] += rhon_mun*flux_n*dKrn_dSn_inside*insideWeight_n; | |
| 482 | |||
| 483 | // partial derivative of the nonwetting phase flux w.r.t. S_n (capillary pressure derivative contribution) | ||
| 484 | 6784 | dI_dI[conti0EqIdx+1][saturationIdx] += tij*dpc_dSn_inside*up_n; | |
| 485 | 3392 | } | |
| 486 | |||
| 487 | /*! | ||
| 488 | * \brief Adds Robin flux derivatives for wetting and nonwetting phase | ||
| 489 | * | ||
| 490 | * \param derivativeMatrices The matrices containing the derivatives | ||
| 491 | * \param problem The problem | ||
| 492 | * \param element The element | ||
| 493 | * \param fvGeometry The finite volume element geometry | ||
| 494 | * \param curElemVolVars The current element volume variables | ||
| 495 | * \param elemFluxVarsCache The element flux variables cache | ||
| 496 | * \param scvf The sub control volume face | ||
| 497 | */ | ||
| 498 | template<class PartialDerivativeMatrices> | ||
| 499 | ✗ | void addRobinFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices, | |
| 500 | const Problem& problem, | ||
| 501 | const Element& element, | ||
| 502 | const FVElementGeometry& fvGeometry, | ||
| 503 | const ElementVolumeVariables& curElemVolVars, | ||
| 504 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 505 | const SubControlVolumeFace& scvf) const | ||
| 506 | ✗ | { /* TODO maybe forward to problem for the user to implement the Robin derivatives?*/ } | |
| 507 | }; | ||
| 508 | |||
| 509 | } // end namespace Dumux | ||
| 510 | |||
| 511 | #endif | ||
| 512 |