| 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 RichardsModel | ||
| 10 | * \brief Element-wise calculation of the Jacobian matrix for problems | ||
| 11 | * using the Richards fully implicit models. | ||
| 12 | */ | ||
| 13 | |||
| 14 | #ifndef DUMUX_RICHARDS_LOCAL_RESIDUAL_HH | ||
| 15 | #define DUMUX_RICHARDS_LOCAL_RESIDUAL_HH | ||
| 16 | |||
| 17 | #include <dumux/common/properties.hh> | ||
| 18 | #include <dumux/common/parameters.hh> | ||
| 19 | #include <dumux/common/typetraits/typetraits.hh> | ||
| 20 | #include <dumux/common/numeqvector.hh> | ||
| 21 | #include <dumux/discretization/method.hh> | ||
| 22 | #include <dumux/discretization/extrusion.hh> | ||
| 23 | #include <dumux/flux/referencesystemformulation.hh> | ||
| 24 | |||
| 25 | namespace Dumux::Detail { | ||
| 26 | // helper structs and functions detecting if the user-defined problem class | ||
| 27 | // implements addRobinFluxDerivatives | ||
| 28 | template <typename T, typename ...Ts> | ||
| 29 | using RobinDerivDetector = decltype(std::declval<T>().addRobinFluxDerivatives(std::declval<Ts>()...)); | ||
| 30 | |||
| 31 | template<class T, typename ...Args> | ||
| 32 | static constexpr bool hasAddRobinFluxDerivatives() | ||
| 33 | { return Dune::Std::is_detected<RobinDerivDetector, T, Args...>::value; } | ||
| 34 | |||
| 35 | } // end namespace Dumux::Detail | ||
| 36 | |||
| 37 | namespace Dumux { | ||
| 38 | |||
| 39 | /*! | ||
| 40 | * \ingroup RichardsModel | ||
| 41 | * \brief Element-wise calculation of the Jacobian matrix for problems | ||
| 42 | * using the Richards fully implicit models. | ||
| 43 | */ | ||
| 44 | template<class TypeTag> | ||
| 45 | class RichardsLocalResidual : public GetPropType<TypeTag, Properties::BaseLocalResidual> | ||
| 46 | { | ||
| 47 | using Implementation = GetPropType<TypeTag, Properties::LocalResidual>; | ||
| 48 | |||
| 49 | using ParentType = GetPropType<TypeTag, Properties::BaseLocalResidual>; | ||
| 50 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 51 | using Problem = GetPropType<TypeTag, Properties::Problem>; | ||
| 52 | using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>; | ||
| 53 | using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>; | ||
| 54 | using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView; | ||
| 55 | using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>; | ||
| 56 | using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView; | ||
| 57 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 58 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 59 | using Extrusion = Extrusion_t<GridGeometry>; | ||
| 60 | using SubControlVolume = typename GridGeometry::SubControlVolume; | ||
| 61 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 62 | using GridView = typename GridGeometry::GridView; | ||
| 63 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 64 | using EnergyLocalResidual = GetPropType<TypeTag, Properties::EnergyLocalResidual>; | ||
| 65 | using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>; | ||
| 66 | using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices; | ||
| 67 | // first index for the mass balance | ||
| 68 | enum { conti0EqIdx = Indices::conti0EqIdx }; | ||
| 69 | |||
| 70 | // checks if the fluid system uses the Richards model index convention | ||
| 71 | static constexpr auto fsCheck = GetPropType<TypeTag, Properties::ModelTraits>::checkFluidSystem(FluidSystem{}); | ||
| 72 | |||
| 73 | // phase & component indices | ||
| 74 | static constexpr auto liquidPhaseIdx = FluidSystem::phase0Idx; | ||
| 75 | static constexpr auto gasPhaseIdx = FluidSystem::phase1Idx; | ||
| 76 | static constexpr auto liquidCompIdx = FluidSystem::comp0Idx; | ||
| 77 | |||
| 78 | //! An element solution that does not compile if the [] operator is used | ||
| 79 | struct InvalidElemSol | ||
| 80 | { | ||
| 81 | template<class Index> | ||
| 82 | double operator[] (const Index i) const | ||
| 83 | { static_assert(AlwaysFalse<Index>::value, "Solution-dependent material parameters not supported with analytical differentiation"); return 0.0; } | ||
| 84 | }; | ||
| 85 | |||
| 86 | public: | ||
| 87 |
2/4✓ Branch 1 taken 28471615 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28471615 times.
✗ Branch 5 not taken.
|
58355838 | using ParentType::ParentType; |
| 88 | |||
| 89 | /*! | ||
| 90 | * \brief Evaluates the rate of change of all conservation | ||
| 91 | * quantites (e.g. phase mass) within a sub-control | ||
| 92 | * volume of a finite volume element for the immiscible models. | ||
| 93 | * \param problem The problem | ||
| 94 | * \param scv The sub control volume | ||
| 95 | * \param volVars The current or previous volVars | ||
| 96 | * \note This function should not include the source and sink terms. | ||
| 97 | * \note The volVars can be different to allow computing | ||
| 98 | * the implicit euler time derivative here | ||
| 99 | */ | ||
| 100 | 3015660 | NumEqVector computeStorage(const Problem& problem, | |
| 101 | const SubControlVolume& scv, | ||
| 102 | const VolumeVariables& volVars) const | ||
| 103 | { | ||
| 104 | // partial time derivative of the phase mass | ||
| 105 | 107182328 | NumEqVector storage(0.0); | |
| 106 | 110197988 | storage[conti0EqIdx] = volVars.porosity() | |
| 107 | 217380316 | * volVars.density(liquidPhaseIdx) | |
| 108 | 110197988 | * volVars.saturation(liquidPhaseIdx); | |
| 109 | |||
| 110 | //! The energy storage in the water, air and solid phase | ||
| 111 | 110197988 | EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, liquidPhaseIdx); | |
| 112 | 110197988 | EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, gasPhaseIdx); | |
| 113 | 110197988 | EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars); | |
| 114 | |||
| 115 | 3015660 | return storage; | |
| 116 | } | ||
| 117 | |||
| 118 | |||
| 119 | /*! | ||
| 120 | * \brief Evaluates the mass flux over a face of a sub control volume. | ||
| 121 | * | ||
| 122 | * \param problem The problem | ||
| 123 | * \param element The current element. | ||
| 124 | * \param fvGeometry The finite-volume geometry | ||
| 125 | * \param elemVolVars The volume variables of the current element | ||
| 126 | * \param scvf The sub control volume face to compute the flux on | ||
| 127 | * \param elemFluxVarsCache The cache related to flux computation | ||
| 128 | */ | ||
| 129 | 2028255 | NumEqVector computeFlux(const Problem& problem, | |
| 130 | const Element& element, | ||
| 131 | const FVElementGeometry& fvGeometry, | ||
| 132 | const ElementVolumeVariables& elemVolVars, | ||
| 133 | const SubControlVolumeFace& scvf, | ||
| 134 | const ElementFluxVariablesCache& elemFluxVarsCache) const | ||
| 135 | { | ||
| 136 | 2028255 | FluxVariables fluxVars; | |
| 137 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
2028255 | fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache); |
| 138 | |||
| 139 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
2028255 | NumEqVector flux(0.0); |
| 140 | // the physical quantities for which we perform upwinding | ||
| 141 | ✗ | auto upwindTerm = [](const auto& volVars) | |
| 142 | 1846419276 | { return volVars.density(liquidPhaseIdx)*volVars.mobility(liquidPhaseIdx); }; | |
| 143 | |||
| 144 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
2028255 | flux[conti0EqIdx] = fluxVars.advectiveFlux(liquidPhaseIdx, upwindTerm); |
| 145 | |||
| 146 | //! Add advective phase energy fluxes for the water phase only. For isothermal model the contribution is zero. | ||
| 147 | 2028255 | EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, liquidPhaseIdx); | |
| 148 | |||
| 149 | //! Add diffusive energy fluxes. For isothermal model the contribution is zero. | ||
| 150 | //! The effective lambda is averaged over both fluid phases and the solid phase | ||
| 151 | 2028255 | EnergyLocalResidual::heatConductionFlux(flux, fluxVars); | |
| 152 | |||
| 153 | 2028255 | return flux; | |
| 154 | } | ||
| 155 | |||
| 156 | /*! | ||
| 157 | * \brief Adds the storage derivative | ||
| 158 | * | ||
| 159 | * \param partialDerivatives The partial derivatives | ||
| 160 | * \param problem The problem | ||
| 161 | * \param element The element | ||
| 162 | * \param fvGeometry The finite volume element geometry | ||
| 163 | * \param curVolVars The current volume variables | ||
| 164 | * \param scv The sub control volume | ||
| 165 | */ | ||
| 166 | template<class PartialDerivativeMatrix> | ||
| 167 | 8490028 | void addStorageDerivatives(PartialDerivativeMatrix& partialDerivatives, | |
| 168 | const Problem& problem, | ||
| 169 | const Element& element, | ||
| 170 | const FVElementGeometry& fvGeometry, | ||
| 171 | const VolumeVariables& curVolVars, | ||
| 172 | const SubControlVolume& scv) const | ||
| 173 | { | ||
| 174 | static_assert(!FluidSystem::isCompressible(0), | ||
| 175 | "richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!"); | ||
| 176 | |||
| 177 |
4/4✓ Branch 0 taken 11 times.
✓ Branch 1 taken 8490017 times.
✓ Branch 2 taken 11 times.
✓ Branch 3 taken 8490017 times.
|
16980056 | const auto poreVolume = Extrusion::volume(fvGeometry, scv)*curVolVars.porosity()*curVolVars.extrusionFactor(); |
| 178 |
3/4✓ Branch 0 taken 11 times.
✓ Branch 1 taken 8490017 times.
✓ Branch 3 taken 11 times.
✗ Branch 4 not taken.
|
8490028 | static const auto rho = curVolVars.density(0); |
| 179 | |||
| 180 | // partial derivative of storage term w.r.t. p_w | ||
| 181 | // d(Sw*rho*phi*V/dt)/dpw = rho*phi*V/dt*dsw/dpw = rho*phi*V/dt*dsw/dpc*dpc/dpw = -rho*phi*V/dt*dsw/dpc | ||
| 182 | 16980056 | const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, InvalidElemSol{}); | |
| 183 |
2/2✓ Branch 1 taken 7728281 times.
✓ Branch 2 taken 761747 times.
|
16218309 | partialDerivatives[conti0EqIdx][0] += -rho*poreVolume/this->timeLoop().timeStepSize()*fluidMatrixInteraction.dsw_dpc(curVolVars.capillaryPressure()); |
| 184 | 8490028 | } | |
| 185 | |||
| 186 | /*! | ||
| 187 | * \brief Adds source derivatives for wetting and nonwetting phase. | ||
| 188 | * | ||
| 189 | * \param partialDerivatives The partial derivatives | ||
| 190 | * \param problem The problem | ||
| 191 | * \param element The element | ||
| 192 | * \param fvGeometry The finite volume element geometry | ||
| 193 | * \param curVolVars The current volume variables | ||
| 194 | * \param scv The sub control volume | ||
| 195 | * | ||
| 196 | * \todo Maybe forward to problem for the user to implement the source derivatives? | ||
| 197 | */ | ||
| 198 | template<class PartialDerivativeMatrix> | ||
| 199 | ✗ | void addSourceDerivatives(PartialDerivativeMatrix& partialDerivatives, | |
| 200 | const Problem& problem, | ||
| 201 | const Element& element, | ||
| 202 | const FVElementGeometry& fvGeometry, | ||
| 203 | const VolumeVariables& curVolVars, | ||
| 204 | const SubControlVolume& scv) const | ||
| 205 | ✗ | { /* TODO maybe forward to problem for the user to implement the source derivatives?*/ } | |
| 206 | |||
| 207 | /*! | ||
| 208 | * \brief Adds flux derivatives for wetting and nonwetting phase for cell-centered FVM using TPFA | ||
| 209 | * | ||
| 210 | * Compute derivatives for the wetting and the nonwetting phase flux with respect to \f$p_w\f$ | ||
| 211 | * and \f$S_n\f$. | ||
| 212 | * | ||
| 213 | * \param derivativeMatrices The partial derivatives | ||
| 214 | * \param problem The problem | ||
| 215 | * \param element The element | ||
| 216 | * \param fvGeometry The finite volume element geometry | ||
| 217 | * \param curElemVolVars The current element volume variables | ||
| 218 | * \param elemFluxVarsCache The element flux variables cache | ||
| 219 | * \param scvf The sub control volume face | ||
| 220 | */ | ||
| 221 | template<class PartialDerivativeMatrices, class T = TypeTag> | ||
| 222 | std::enable_if_t<GetPropType<T, Properties::GridGeometry>::discMethod == DiscretizationMethods::cctpfa, void> | ||
| 223 | 16338378 | addFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices, | |
| 224 | const Problem& problem, | ||
| 225 | const Element& element, | ||
| 226 | const FVElementGeometry& fvGeometry, | ||
| 227 | const ElementVolumeVariables& curElemVolVars, | ||
| 228 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 229 | const SubControlVolumeFace& scvf) const | ||
| 230 | { | ||
| 231 | static_assert(!FluidSystem::isCompressible(0), | ||
| 232 | "richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!"); | ||
| 233 | static_assert(FluidSystem::viscosityIsConstant(0), | ||
| 234 | "richards/localresidual.hh: Analytic Jacobian only supports fluids with constant viscosity!"); | ||
| 235 | |||
| 236 | // get references to the two participating vol vars & parameters | ||
| 237 | 16338378 | const auto insideScvIdx = scvf.insideScvIdx(); | |
| 238 | 16338378 | const auto outsideScvIdx = scvf.outsideScvIdx(); | |
| 239 | 16338378 | const auto outsideElement = fvGeometry.gridGeometry().element(outsideScvIdx); | |
| 240 |
2/2✓ Branch 0 taken 53872 times.
✓ Branch 1 taken 16284506 times.
|
16338378 | const auto& insideScv = fvGeometry.scv(insideScvIdx); |
| 241 |
2/2✓ Branch 0 taken 16284506 times.
✓ Branch 1 taken 53872 times.
|
16338378 | const auto& outsideScv = fvGeometry.scv(outsideScvIdx); |
| 242 | 16338378 | const auto& insideVolVars = curElemVolVars[insideScvIdx]; | |
| 243 | 16338378 | const auto& outsideVolVars = curElemVolVars[outsideScvIdx]; | |
| 244 | |||
| 245 | // some quantities to be reused (rho & mu are constant and thus equal for all cells) | ||
| 246 |
3/4✓ Branch 0 taken 8 times.
✓ Branch 1 taken 16338370 times.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
|
16338378 | static const auto rho = insideVolVars.density(0); |
| 247 |
3/4✓ Branch 0 taken 8 times.
✓ Branch 1 taken 16338370 times.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
|
16338378 | static const auto mu = insideVolVars.viscosity(0); |
| 248 |
3/4✓ Branch 0 taken 8 times.
✓ Branch 1 taken 16338370 times.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
|
16338378 | static const auto rho_mu = rho/mu; |
| 249 | |||
| 250 | // upwind term | ||
| 251 | // evaluate the current wetting phase Darcy flux and resulting upwind weights | ||
| 252 | using AdvectionType = GetPropType<TypeTag, Properties::AdvectionType>; | ||
| 253 |
6/8✓ Branch 0 taken 9 times.
✓ Branch 1 taken 16338369 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 1 times.
✓ Branch 6 taken 8 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 8 times.
✗ Branch 10 not taken.
|
16338378 | static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight"); |
| 254 | 16338378 | const auto flux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache); | |
| 255 |
4/4✓ Branch 0 taken 7825063 times.
✓ Branch 1 taken 8513315 times.
✓ Branch 2 taken 7825063 times.
✓ Branch 3 taken 8513315 times.
|
32676756 | const auto insideWeight = std::signbit(flux) ? (1.0 - upwindWeight) : upwindWeight; |
| 256 | 16338378 | const auto outsideWeight = 1.0 - insideWeight; | |
| 257 |
4/4✓ Branch 0 taken 14800917 times.
✓ Branch 1 taken 1483589 times.
✓ Branch 2 taken 14800917 times.
✓ Branch 3 taken 1483589 times.
|
32676756 | const auto upwindTerm = rho*insideVolVars.mobility(0)*insideWeight + rho*outsideVolVars.mobility(0)*outsideWeight; |
| 258 | |||
| 259 |
4/4✓ Branch 0 taken 14800917 times.
✓ Branch 1 taken 1483589 times.
✓ Branch 2 taken 14800917 times.
✓ Branch 3 taken 1483589 times.
|
32676756 | const auto insideFluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, insideScv, InvalidElemSol{}); |
| 260 |
4/4✓ Branch 0 taken 14800917 times.
✓ Branch 1 taken 1483589 times.
✓ Branch 2 taken 14800917 times.
✓ Branch 3 taken 1483589 times.
|
32676756 | const auto outsideFluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(outsideElement, outsideScv, InvalidElemSol{}); |
| 261 | |||
| 262 | // material law derivatives | ||
| 263 |
2/2✓ Branch 0 taken 14854789 times.
✓ Branch 1 taken 1483589 times.
|
16338378 | const auto insideSw = insideVolVars.saturation(0); |
| 264 |
2/2✓ Branch 0 taken 14854789 times.
✓ Branch 1 taken 1483589 times.
|
16338378 | const auto outsideSw = outsideVolVars.saturation(0); |
| 265 |
2/2✓ Branch 0 taken 14854789 times.
✓ Branch 1 taken 1483589 times.
|
16338378 | const auto insidePc = insideVolVars.capillaryPressure(); |
| 266 |
2/2✓ Branch 0 taken 14854789 times.
✓ Branch 1 taken 1483589 times.
|
16338378 | const auto outsidePc = outsideVolVars.capillaryPressure(); |
| 267 | 16338378 | const auto dkrw_dsw_inside = insideFluidMatrixInteraction.dkrw_dsw(insideSw); | |
| 268 | 16338378 | const auto dkrw_dsw_outside = outsideFluidMatrixInteraction.dkrw_dsw(outsideSw); | |
| 269 | 16338378 | const auto dsw_dpw_inside = -insideFluidMatrixInteraction.dsw_dpc(insidePc); | |
| 270 | 16338378 | const auto dsw_dpw_outside = -outsideFluidMatrixInteraction.dsw_dpc(outsidePc); | |
| 271 | |||
| 272 | // the transmissibility | ||
| 273 | 16338378 | const auto tij = elemFluxVarsCache[scvf].advectionTij(); | |
| 274 | |||
| 275 | // get references to the two participating derivative matrices | ||
| 276 | 16338378 | auto& dI_dI = derivativeMatrices[insideScvIdx]; | |
| 277 | 16338378 | auto& dI_dJ = derivativeMatrices[outsideScvIdx]; | |
| 278 | |||
| 279 | // partial derivative of the wetting phase flux w.r.t. pw | ||
| 280 | 16338378 | dI_dI[conti0EqIdx][0] += tij*upwindTerm + rho_mu*flux*insideWeight*dkrw_dsw_inside*dsw_dpw_inside; | |
| 281 | 16338378 | dI_dJ[conti0EqIdx][0] += -tij*upwindTerm + rho_mu*flux*outsideWeight*dkrw_dsw_outside*dsw_dpw_outside; | |
| 282 | 16338378 | } | |
| 283 | |||
| 284 | /*! | ||
| 285 | * \brief Adds flux derivatives for box method | ||
| 286 | * | ||
| 287 | * \param A The Jacobian Matrix | ||
| 288 | * \param problem The problem | ||
| 289 | * \param element The element | ||
| 290 | * \param fvGeometry The finite volume element geometry | ||
| 291 | * \param curElemVolVars The current element volume variables | ||
| 292 | * \param elemFluxVarsCache The element flux variables cache | ||
| 293 | * \param scvf The sub control volume face | ||
| 294 | */ | ||
| 295 | template<class JacobianMatrix, class T = TypeTag> | ||
| 296 | std::enable_if_t<GetPropType<T, Properties::GridGeometry>::discMethod == DiscretizationMethods::box, void> | ||
| 297 | 200520 | addFluxDerivatives(JacobianMatrix& A, | |
| 298 | const Problem& problem, | ||
| 299 | const Element& element, | ||
| 300 | const FVElementGeometry& fvGeometry, | ||
| 301 | const ElementVolumeVariables& curElemVolVars, | ||
| 302 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 303 | const SubControlVolumeFace& scvf) const | ||
| 304 | { | ||
| 305 | static_assert(!FluidSystem::isCompressible(0), | ||
| 306 | "richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!"); | ||
| 307 | static_assert(FluidSystem::viscosityIsConstant(0), | ||
| 308 | "richards/localresidual.hh: Analytic Jacobian only supports fluids with constant viscosity!"); | ||
| 309 | |||
| 310 | // get references to the two participating vol vars & parameters | ||
| 311 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 200520 times.
|
200520 | const auto insideScvIdx = scvf.insideScvIdx(); |
| 312 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 200520 times.
|
200520 | const auto outsideScvIdx = scvf.outsideScvIdx(); |
| 313 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 200515 times.
|
200520 | const auto& insideScv = fvGeometry.scv(insideScvIdx); |
| 314 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 200515 times.
|
200520 | const auto& outsideScv = fvGeometry.scv(outsideScvIdx); |
| 315 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 200515 times.
|
200520 | const auto& insideVolVars = curElemVolVars[insideScvIdx]; |
| 316 |
2/2✓ Branch 0 taken 5 times.
✓ Branch 1 taken 200515 times.
|
200520 | const auto& outsideVolVars = curElemVolVars[outsideScvIdx]; |
| 317 | |||
| 318 | // some quantities to be reused (rho & mu are constant and thus equal for all cells) | ||
| 319 |
3/4✓ Branch 0 taken 5 times.
✓ Branch 1 taken 200515 times.
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
|
200520 | static const auto rho = insideVolVars.density(0); |
| 320 |
3/4✓ Branch 0 taken 5 times.
✓ Branch 1 taken 200515 times.
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
|
200520 | static const auto mu = insideVolVars.viscosity(0); |
| 321 |
3/4✓ Branch 0 taken 5 times.
✓ Branch 1 taken 200515 times.
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
|
200520 | static const auto rho_mu = rho/mu; |
| 322 | |||
| 323 | // upwind term | ||
| 324 | // evaluate the current wetting phase Darcy flux and resulting upwind weights | ||
| 325 | using AdvectionType = GetPropType<TypeTag, Properties::AdvectionType>; | ||
| 326 |
5/8✓ Branch 0 taken 5 times.
✓ Branch 1 taken 200515 times.
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 5 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 5 times.
✗ Branch 10 not taken.
|
200520 | static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight"); |
| 327 | 200520 | const auto flux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache); | |
| 328 |
4/4✓ Branch 0 taken 192011 times.
✓ Branch 1 taken 8509 times.
✓ Branch 2 taken 192011 times.
✓ Branch 3 taken 8509 times.
|
401040 | const auto insideWeight = std::signbit(flux) ? (1.0 - upwindWeight) : upwindWeight; |
| 329 | 200520 | const auto outsideWeight = 1.0 - insideWeight; | |
| 330 |
4/4✓ Branch 0 taken 123121 times.
✓ Branch 1 taken 8279 times.
✓ Branch 2 taken 123121 times.
✓ Branch 3 taken 8279 times.
|
401040 | const auto upwindTerm = rho*insideVolVars.mobility(0)*insideWeight + rho*outsideVolVars.mobility(0)*outsideWeight; |
| 331 | |||
| 332 |
4/4✓ Branch 0 taken 123121 times.
✓ Branch 1 taken 8279 times.
✓ Branch 2 taken 123121 times.
✓ Branch 3 taken 8279 times.
|
401040 | const auto insideFluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, insideScv, InvalidElemSol{}); |
| 333 |
4/4✓ Branch 0 taken 123121 times.
✓ Branch 1 taken 8279 times.
✓ Branch 2 taken 123121 times.
✓ Branch 3 taken 8279 times.
|
401040 | const auto outsideFluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, outsideScv, InvalidElemSol{}); |
| 334 | |||
| 335 | // material law derivatives | ||
| 336 |
2/2✓ Branch 0 taken 192241 times.
✓ Branch 1 taken 8279 times.
|
200520 | const auto insideSw = insideVolVars.saturation(0); |
| 337 |
2/2✓ Branch 0 taken 192241 times.
✓ Branch 1 taken 8279 times.
|
200520 | const auto outsideSw = outsideVolVars.saturation(0); |
| 338 |
2/2✓ Branch 0 taken 192241 times.
✓ Branch 1 taken 8279 times.
|
200520 | const auto insidePc = insideVolVars.capillaryPressure(); |
| 339 |
2/2✓ Branch 0 taken 192232 times.
✓ Branch 1 taken 8288 times.
|
200520 | const auto outsidePc = outsideVolVars.capillaryPressure(); |
| 340 | 200520 | const auto dkrw_dsw_inside = insideFluidMatrixInteraction.dkrw_dsw(insideSw); | |
| 341 | 200520 | const auto dkrw_dsw_outside = outsideFluidMatrixInteraction.dkrw_dsw(outsideSw); | |
| 342 | 200520 | const auto dsw_dpw_inside = -insideFluidMatrixInteraction.dsw_dpc(insidePc); | |
| 343 | 200520 | const auto dsw_dpw_outside = -outsideFluidMatrixInteraction.dsw_dpc(outsidePc); | |
| 344 | |||
| 345 | // so far it was the same as for tpfa | ||
| 346 | // the transmissibilities (flux derivatives with respect to all pw-dofs on the element) | ||
| 347 | ✗ | const auto ti = AdvectionType::calculateTransmissibilities( | |
| 348 | 401040 | problem, element, fvGeometry, curElemVolVars, scvf, elemFluxVarsCache[scvf] | |
| 349 | ); | ||
| 350 | |||
| 351 | // get the rows of the jacobian matrix for the inside/outside scv | ||
| 352 | 200520 | auto& dI_dJ_inside = A[insideScv.dofIndex()]; | |
| 353 | 200520 | auto& dI_dJ_outside = A[outsideScv.dofIndex()]; | |
| 354 | |||
| 355 | // add the partial derivatives w.r.t all scvs in the element | ||
| 356 |
4/4✓ Branch 0 taken 539280 times.
✓ Branch 1 taken 200520 times.
✓ Branch 2 taken 539280 times.
✓ Branch 3 taken 200520 times.
|
1479600 | for (const auto& scvJ : scvs(fvGeometry)) |
| 357 | { | ||
| 358 | // the transmissibily associated with the scvJ | ||
| 359 |
1/2✓ Branch 1 taken 539280 times.
✗ Branch 2 not taken.
|
539280 | const auto& tj = ti[scvJ.indexInElement()]; |
| 360 | 539280 | const auto globalJ = scvJ.dofIndex(); | |
| 361 | |||
| 362 | // partial derivative of the wetting phase flux w.r.t. p_w | ||
| 363 |
2/4✓ Branch 1 taken 539280 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 539280 times.
✗ Branch 5 not taken.
|
539280 | dI_dJ_inside[globalJ][conti0EqIdx][0] += tj*upwindTerm; |
| 364 |
1/2✓ Branch 1 taken 539280 times.
✗ Branch 2 not taken.
|
539280 | dI_dJ_outside[globalJ][conti0EqIdx][0] += -tj*upwindTerm; |
| 365 | } | ||
| 366 | |||
| 367 | 200520 | const auto upwindContributionInside = rho_mu*flux*insideWeight*dkrw_dsw_inside*dsw_dpw_inside; | |
| 368 | 200520 | const auto upwindContributionOutside = rho_mu*flux*outsideWeight*dkrw_dsw_outside*dsw_dpw_outside; | |
| 369 | |||
| 370 | // additional contribution of the upwind term only for inside and outside dof | ||
| 371 |
3/6✓ Branch 1 taken 200520 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 200520 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 200520 times.
✗ Branch 8 not taken.
|
401040 | A[insideScv.dofIndex()][insideScv.dofIndex()][conti0EqIdx][0] += upwindContributionInside; |
| 372 |
3/6✓ Branch 1 taken 200520 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 200520 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 200520 times.
✗ Branch 8 not taken.
|
401040 | A[insideScv.dofIndex()][outsideScv.dofIndex()][conti0EqIdx][0] += upwindContributionOutside; |
| 373 | |||
| 374 |
3/6✓ Branch 1 taken 200520 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 200520 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 200520 times.
✗ Branch 8 not taken.
|
401040 | A[outsideScv.dofIndex()][insideScv.dofIndex()][conti0EqIdx][0] -= upwindContributionInside; |
| 375 |
4/8✓ Branch 1 taken 200520 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 200520 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 200520 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 200520 times.
✗ Branch 9 not taken.
|
401040 | A[outsideScv.dofIndex()][outsideScv.dofIndex()][conti0EqIdx][0] -= upwindContributionOutside; |
| 376 | 200520 | } | |
| 377 | |||
| 378 | /*! | ||
| 379 | * \brief Adds cell-centered Dirichlet flux derivatives for wetting and nonwetting phase | ||
| 380 | * | ||
| 381 | * Compute derivatives for the wetting and the nonwetting phase flux with respect to \f$p_w\f$ | ||
| 382 | * and \f$S_n\f$. | ||
| 383 | * | ||
| 384 | * \param derivativeMatrices The matrices containing the derivatives | ||
| 385 | * \param problem The problem | ||
| 386 | * \param element The element | ||
| 387 | * \param fvGeometry The finite volume element geometry | ||
| 388 | * \param curElemVolVars The current element volume variables | ||
| 389 | * \param elemFluxVarsCache The element flux variables cache | ||
| 390 | * \param scvf The sub control volume face | ||
| 391 | */ | ||
| 392 | template<class PartialDerivativeMatrices> | ||
| 393 | 1184 | void addCCDirichletFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices, | |
| 394 | const Problem& problem, | ||
| 395 | const Element& element, | ||
| 396 | const FVElementGeometry& fvGeometry, | ||
| 397 | const ElementVolumeVariables& curElemVolVars, | ||
| 398 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 399 | const SubControlVolumeFace& scvf) const | ||
| 400 | { | ||
| 401 | static_assert(!FluidSystem::isCompressible(0), | ||
| 402 | "richards/localresidual.hh: Analytic Jacobian only supports incompressible fluids!"); | ||
| 403 | static_assert(FluidSystem::viscosityIsConstant(0), | ||
| 404 | "richards/localresidual.hh: Analytic Jacobian only supports fluids with constant viscosity!"); | ||
| 405 | |||
| 406 | |||
| 407 | // get references to the two participating vol vars & parameters | ||
| 408 |
1/2✓ Branch 0 taken 1184 times.
✗ Branch 1 not taken.
|
1184 | const auto insideScvIdx = scvf.insideScvIdx(); |
| 409 |
1/2✓ Branch 0 taken 1184 times.
✗ Branch 1 not taken.
|
1184 | const auto& insideScv = fvGeometry.scv(insideScvIdx); |
| 410 | 1184 | const auto& insideVolVars = curElemVolVars[insideScvIdx]; | |
| 411 | 2368 | const auto& outsideVolVars = curElemVolVars[scvf.outsideScvIdx()]; | |
| 412 |
0/4✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
2368 | const auto insideFluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, insideScv, InvalidElemSol{}); |
| 413 | |||
| 414 | // some quantities to be reused (rho & mu are constant and thus equal for all cells) | ||
| 415 |
3/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1183 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
1184 | static const auto rho = insideVolVars.density(0); |
| 416 |
3/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1183 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
1184 | static const auto mu = insideVolVars.viscosity(0); |
| 417 |
3/4✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1183 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
1184 | static const auto rho_mu = rho/mu; |
| 418 | |||
| 419 | // upwind term | ||
| 420 | // evaluate the current wetting phase Darcy flux and resulting upwind weights | ||
| 421 | using AdvectionType = GetPropType<TypeTag, Properties::AdvectionType>; | ||
| 422 |
5/8✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1183 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.
|
1184 | static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight"); |
| 423 | 1184 | const auto flux = AdvectionType::flux(problem, element, fvGeometry, curElemVolVars, scvf, 0, elemFluxVarsCache); | |
| 424 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 1184 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1184 times.
|
2368 | const auto insideWeight = std::signbit(flux) ? (1.0 - upwindWeight) : upwindWeight; |
| 425 | 1184 | const auto outsideWeight = 1.0 - insideWeight; | |
| 426 |
2/4✓ Branch 0 taken 1184 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1184 times.
✗ Branch 3 not taken.
|
2368 | const auto upwindTerm = rho*insideVolVars.mobility(0)*insideWeight + rho*outsideVolVars.mobility(0)*outsideWeight; |
| 427 | |||
| 428 | // material law derivatives | ||
| 429 |
1/2✓ Branch 0 taken 1184 times.
✗ Branch 1 not taken.
|
1184 | const auto insideSw = insideVolVars.saturation(0); |
| 430 |
1/2✓ Branch 0 taken 1184 times.
✗ Branch 1 not taken.
|
1184 | const auto insidePc = insideVolVars.capillaryPressure(); |
| 431 | 1184 | const auto dkrw_dsw_inside = insideFluidMatrixInteraction.dkrw_dsw(insideSw); | |
| 432 | 1184 | const auto dsw_dpw_inside = -insideFluidMatrixInteraction.dsw_dpc(insidePc); | |
| 433 | |||
| 434 | // the transmissibility | ||
| 435 | 1184 | const auto tij = elemFluxVarsCache[scvf].advectionTij(); | |
| 436 | |||
| 437 | // partial derivative of the wetting phase flux w.r.t. pw | ||
| 438 | 1184 | derivativeMatrices[insideScvIdx][conti0EqIdx][0] += tij*upwindTerm + rho_mu*flux*insideWeight*dkrw_dsw_inside*dsw_dpw_inside; | |
| 439 | 1184 | } | |
| 440 | |||
| 441 | /*! | ||
| 442 | * \brief Adds Robin flux derivatives for wetting and nonwetting phase | ||
| 443 | * | ||
| 444 | * \param derivativeMatrices The matrices containing the derivatives | ||
| 445 | * \param problem The problem | ||
| 446 | * \param element The element | ||
| 447 | * \param fvGeometry The finite volume element geometry | ||
| 448 | * \param curElemVolVars The current element volume variables | ||
| 449 | * \param elemFluxVarsCache The element flux variables cache | ||
| 450 | * \param scvf The sub control volume face | ||
| 451 | */ | ||
| 452 | template<class PartialDerivativeMatrices> | ||
| 453 | ✗ | void addRobinFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices, | |
| 454 | const Problem& problem, | ||
| 455 | const Element& element, | ||
| 456 | const FVElementGeometry& fvGeometry, | ||
| 457 | const ElementVolumeVariables& curElemVolVars, | ||
| 458 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 459 | const SubControlVolumeFace& scvf) const | ||
| 460 | { | ||
| 461 | if constexpr(Detail::hasAddRobinFluxDerivatives<Problem, | ||
| 462 | PartialDerivativeMatrices&, Element, FVElementGeometry, | ||
| 463 | ElementVolumeVariables, ElementFluxVariablesCache, SubControlVolumeFace>() | ||
| 464 | ) | ||
| 465 | 38094 | problem.addRobinFluxDerivatives(derivativeMatrices, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf); | |
| 466 | ✗ | } | |
| 467 | |||
| 468 | private: | ||
| 469 | Implementation *asImp_() | ||
| 470 | { return static_cast<Implementation *> (this); } | ||
| 471 | |||
| 472 | const Implementation *asImp_() const | ||
| 473 | { return static_cast<const Implementation *> (this); } | ||
| 474 | }; | ||
| 475 | |||
| 476 | } // end namespace Dumux | ||
| 477 | |||
| 478 | #endif | ||
| 479 |