| 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-FileCopyrightText: 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 TracerModel | ||
| 10 | * \brief Element-wise calculation of the local residual for problems | ||
| 11 | * using fully implicit tracer model. | ||
| 12 | */ | ||
| 13 | |||
| 14 | #ifndef DUMUX_TRACER_LOCAL_RESIDUAL_HH | ||
| 15 | #define DUMUX_TRACER_LOCAL_RESIDUAL_HH | ||
| 16 | |||
| 17 | #include <dune/common/exceptions.hh> | ||
| 18 | |||
| 19 | #include <dumux/common/properties.hh> | ||
| 20 | #include <dumux/common/parameters.hh> | ||
| 21 | #include <dumux/common/numeqvector.hh> | ||
| 22 | #include <dumux/discretization/method.hh> | ||
| 23 | #include <dumux/discretization/extrusion.hh> | ||
| 24 | #include <dumux/discretization/defaultlocaloperator.hh> | ||
| 25 | #include <dumux/flux/referencesystemformulation.hh> | ||
| 26 | |||
| 27 | namespace Dumux { | ||
| 28 | |||
| 29 | /*! | ||
| 30 | * \ingroup TracerModel | ||
| 31 | * \brief Element-wise calculation of the local residual for problems | ||
| 32 | * using fully implicit tracer model. | ||
| 33 | */ | ||
| 34 | template<class TypeTag> | ||
| 35 | class TracerLocalResidual | ||
| 36 | : public DiscretizationDefaultLocalOperator<TypeTag> | ||
| 37 | { | ||
| 38 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 39 | using ParentType = DiscretizationDefaultLocalOperator<TypeTag>; | ||
| 40 | using Problem = GetPropType<TypeTag, Properties::Problem>; | ||
| 41 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 42 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 43 | using SubControlVolume = typename GridGeometry::SubControlVolume; | ||
| 44 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 45 | using Extrusion = Extrusion_t<GridGeometry>; | ||
| 46 | using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>; | ||
| 47 | using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>; | ||
| 48 | using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView; | ||
| 49 | using GridView = typename GridGeometry::GridView; | ||
| 50 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 51 | using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView; | ||
| 52 | using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>; | ||
| 53 | using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>; | ||
| 54 | using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>; | ||
| 55 | |||
| 56 | static constexpr int numComponents = ModelTraits::numFluidComponents(); | ||
| 57 | static constexpr bool useMoles = getPropValue<TypeTag, Properties::UseMoles>(); | ||
| 58 | static constexpr int phaseIdx = 0; | ||
| 59 | |||
| 60 | public: | ||
| 61 |
3/6✓ Branch 1 taken 12660300 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 950000 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2892000 times.
✗ Branch 8 not taken.
|
17512300 | using ParentType::ParentType; |
| 62 | |||
| 63 | /*! | ||
| 64 | * \brief Evaluates the amount of all conservation quantities | ||
| 65 | * (e.g. phase mass) within a sub-control volume. | ||
| 66 | * | ||
| 67 | * The result should be averaged over the volume (e.g. phase mass | ||
| 68 | * inside a sub control volume divided by the volume) | ||
| 69 | * | ||
| 70 | * \param problem The problem | ||
| 71 | * \param scv The sub control volume | ||
| 72 | * \param volVars The primary and secondary variables on the scv | ||
| 73 | */ | ||
| 74 | 76608070 | NumEqVector computeStorage(const Problem& problem, | |
| 75 | const SubControlVolume& scv, | ||
| 76 | const VolumeVariables& volVars) const | ||
| 77 | { | ||
| 78 |
5/8✓ Branch 0 taken 20302035 times.
✓ Branch 1 taken 12000000 times.
✓ Branch 2 taken 2202035 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 21052000 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 21052000 times.
✗ Branch 7 not taken.
|
76608070 | NumEqVector storage(0.0); |
| 79 | |||
| 80 | // regularize saturation so we don't get singular matrices when the saturation is zero | ||
| 81 | // note that the fluxes will still be zero (zero effective diffusion coefficient), | ||
| 82 | // and we still solve the equation storage = 0 yielding the correct result | ||
| 83 | using std::max; | ||
| 84 |
3/4✓ Branch 0 taken 20302035 times.
✓ Branch 1 taken 12000000 times.
✓ Branch 2 taken 21052000 times.
✗ Branch 3 not taken.
|
53354035 | const Scalar saturation = max(1e-8, volVars.saturation(phaseIdx)); |
| 85 | |||
| 86 | // formulation with mole balances | ||
| 87 | if (useMoles) | ||
| 88 | { | ||
| 89 |
3/4✓ Branch 0 taken 5225235 times.
✓ Branch 1 taken 2525000 times.
✓ Branch 2 taken 808000 times.
✗ Branch 3 not taken.
|
9541470 | for (int compIdx = 0; compIdx < numComponents; ++compIdx) |
| 90 |
2/4✓ Branch 0 taken 175235 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 808000 times.
✗ Branch 3 not taken.
|
7016470 | storage[compIdx] += volVars.porosity() |
| 91 |
2/4✓ Branch 0 taken 175235 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 808000 times.
✗ Branch 3 not taken.
|
7016470 | * volVars.molarDensity(phaseIdx) |
| 92 | 7016470 | * volVars.moleFraction(phaseIdx, compIdx) | |
| 93 |
2/4✓ Branch 0 taken 175235 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 808000 times.
✗ Branch 3 not taken.
|
7016470 | * saturation; |
| 94 | } | ||
| 95 | // formulation with mass balances | ||
| 96 | else | ||
| 97 | { | ||
| 98 |
3/4✓ Branch 0 taken 57176800 times.
✓ Branch 1 taken 27575000 times.
✓ Branch 2 taken 20244000 times.
✗ Branch 3 not taken.
|
127266600 | for (int compIdx = 0; compIdx < numComponents; ++compIdx) |
| 99 |
2/4✓ Branch 0 taken 2026800 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 20244000 times.
✗ Branch 3 not taken.
|
99691600 | storage[compIdx] += volVars.porosity() |
| 100 |
2/4✓ Branch 0 taken 2026800 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 20244000 times.
✗ Branch 3 not taken.
|
99691600 | * volVars.density(phaseIdx) |
| 101 | 99691600 | * volVars.massFraction(phaseIdx, compIdx) | |
| 102 |
2/4✓ Branch 0 taken 2026800 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 20244000 times.
✗ Branch 3 not taken.
|
99691600 | * saturation; |
| 103 | } | ||
| 104 | |||
| 105 | 30100000 | return storage; | |
| 106 | } | ||
| 107 | |||
| 108 | /*! | ||
| 109 | * \brief Evaluates the total flux of all conservation quantities | ||
| 110 | * over a face of a sub-control volume. | ||
| 111 | * | ||
| 112 | * \param problem The problem | ||
| 113 | * \param element The element | ||
| 114 | * \param fvGeometry The finite volume geometry context | ||
| 115 | * \param elemVolVars The volume variables for all flux stencil elements | ||
| 116 | * \param scvf The sub control volume face | ||
| 117 | * \param elemFluxVarsCache The cache related to flux computation | ||
| 118 | */ | ||
| 119 | 89571492 | NumEqVector computeFlux(const Problem& problem, | |
| 120 | const Element& element, | ||
| 121 | const FVElementGeometry& fvGeometry, | ||
| 122 | const ElementVolumeVariables& elemVolVars, | ||
| 123 | const SubControlVolumeFace& scvf, | ||
| 124 | const ElementFluxVariablesCache& elemFluxVarsCache) const | ||
| 125 | { | ||
| 126 | 89571492 | FluxVariables fluxVars; | |
| 127 | 89571492 | fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache); | |
| 128 | static constexpr auto referenceSystemFormulationDiffusion = FluxVariables::MolecularDiffusionType::referenceSystemFormulation(); | ||
| 129 | |||
| 130 | 89571492 | NumEqVector flux(0.0); | |
| 131 | |||
| 132 | 156502911 | const auto massOrMoleDensity = [](const auto& volVars, const int phaseIdx) | |
| 133 | 156502911 | { return useMoles ? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); }; | |
| 134 | |||
| 135 | 156502911 | const auto massOrMoleFraction = [](const auto& volVars, const int phaseIdx, const int compIdx) | |
| 136 | 156502911 | { return useMoles ? volVars.moleFraction(phaseIdx, compIdx) : volVars.massFraction(phaseIdx, compIdx); }; | |
| 137 | |||
| 138 | // Check for the reference system and adapt units of the flux accordingly. | ||
| 139 | 156490692 | const auto adaptFluxUnits = [](const Scalar referenceFlux, const Scalar molarMass, | |
| 140 | const ReferenceSystemFormulation referenceSystemFormulation) | ||
| 141 | { | ||
| 142 | if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged) | ||
| 143 | return useMoles ? referenceFlux/molarMass | ||
| 144 | 9041492 | : referenceFlux; | |
| 145 | else if (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged) | ||
| 146 | return useMoles ? referenceFlux | ||
| 147 | : referenceFlux*molarMass; | ||
| 148 | else | ||
| 149 | DUNE_THROW(Dune::NotImplemented, "other reference systems than mass and molar averaged are not implemented"); | ||
| 150 | }; | ||
| 151 | |||
| 152 | 77163492 | const auto diffusiveFluxes = fluxVars.molecularDiffusionFlux(phaseIdx); | |
| 153 | 1999800 | auto dispersiveFluxes = decltype(diffusiveFluxes)(0.0); | |
| 154 | if constexpr (ModelTraits::enableCompositionalDispersion()) | ||
| 155 | 1999800 | dispersiveFluxes = fluxVars.compositionalDispersionFlux(phaseIdx); | |
| 156 | |||
| 157 |
3/3✓ Branch 1 taken 137041000 times.
✓ Branch 2 taken 70121800 times.
✓ Branch 0 taken 3999600 times.
|
246062184 | for (int compIdx = 0; compIdx < numComponents; ++compIdx) |
| 158 | { | ||
| 159 | // the physical quantities for which we perform upwinding | ||
| 160 | 312993603 | auto upwindTerm = [&massOrMoleDensity, &massOrMoleFraction, compIdx](const auto& volVars) | |
| 161 | 156502911 | { return massOrMoleDensity(volVars, phaseIdx)*massOrMoleFraction(volVars, phaseIdx, compIdx); }; | |
| 162 | |||
| 163 | // advective fluxes | ||
| 164 | 156490692 | flux[compIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm); | |
| 165 | // diffusive and dispersive fluxes | ||
| 166 | 156490692 | flux[compIdx] += adaptFluxUnits(diffusiveFluxes[compIdx], | |
| 167 | FluidSystem::molarMass(compIdx), | ||
| 168 | referenceSystemFormulationDiffusion); | ||
| 169 | if constexpr (ModelTraits::enableCompositionalDispersion()) | ||
| 170 | { | ||
| 171 | static constexpr auto referenceSystemFormulationDispersion = | ||
| 172 | FluxVariables::DispersionFluxType::referenceSystemFormulation(); | ||
| 173 | 3999600 | flux[compIdx] += adaptFluxUnits(dispersiveFluxes[compIdx], | |
| 174 | FluidSystem::molarMass(compIdx), | ||
| 175 | referenceSystemFormulationDispersion); | ||
| 176 | } | ||
| 177 | } | ||
| 178 | |||
| 179 | 89571492 | return flux; | |
| 180 | } | ||
| 181 | |||
| 182 | /*! | ||
| 183 | * \brief TODO docme! | ||
| 184 | * | ||
| 185 | * \param partialDerivatives TODO docme! | ||
| 186 | * \param problem The problem | ||
| 187 | * \param element The element | ||
| 188 | * \param fvGeometry The finite volume geometry context | ||
| 189 | * \param curVolVars The current volume variables | ||
| 190 | * \param scv The sub control volume | ||
| 191 | */ | ||
| 192 | template<class PartialDerivativeMatrix> | ||
| 193 |
1/2✓ Branch 0 taken 1376800 times.
✗ Branch 1 not taken.
|
1376800 | void addStorageDerivatives(PartialDerivativeMatrix& partialDerivatives, |
| 194 | const Problem& problem, | ||
| 195 | const Element& element, | ||
| 196 | const FVElementGeometry& fvGeometry, | ||
| 197 | const VolumeVariables& curVolVars, | ||
| 198 | const SubControlVolume& scv) const | ||
| 199 | { | ||
| 200 | // regularize saturation so we don't get singular matrices when the saturation is zero | ||
| 201 | // note that the fluxes will still be zero (zero effective diffusion coefficient), | ||
| 202 | // and we still solve the equation storage = 0 yielding the correct result | ||
| 203 | using std::max; | ||
| 204 |
1/2✓ Branch 0 taken 1376800 times.
✗ Branch 1 not taken.
|
2753600 | const auto saturation = max(1e-8, curVolVars.saturation(phaseIdx)); |
| 205 | |||
| 206 | 1376800 | const auto porosity = curVolVars.porosity(); | |
| 207 | 1376800 | const auto rho = useMoles ? curVolVars.molarDensity() : curVolVars.density(); | |
| 208 | 1376800 | const auto d_storage = Extrusion::volume(fvGeometry, scv)*porosity*rho*saturation/this->timeLoop().timeStepSize(); | |
| 209 | |||
| 210 |
2/2✓ Branch 0 taken 1426800 times.
✓ Branch 1 taken 1376800 times.
|
2803600 | for (int compIdx = 0; compIdx < numComponents; ++compIdx) |
| 211 | 1426800 | partialDerivatives[compIdx][compIdx] += d_storage; | |
| 212 | 1376800 | } | |
| 213 | |||
| 214 | /*! | ||
| 215 | * \brief TODO docme! | ||
| 216 | * | ||
| 217 | * \param partialDerivatives TODO docme! | ||
| 218 | * \param problem The problem | ||
| 219 | * \param element The element | ||
| 220 | * \param fvGeometry The finite volume geometry context | ||
| 221 | * \param curVolVars The current volume variables | ||
| 222 | * \param scv The sub control volume | ||
| 223 | */ | ||
| 224 | template<class PartialDerivativeMatrix> | ||
| 225 | void addSourceDerivatives(PartialDerivativeMatrix& partialDerivatives, | ||
| 226 | const Problem& problem, | ||
| 227 | const Element& element, | ||
| 228 | const FVElementGeometry& fvGeometry, | ||
| 229 | const VolumeVariables& curVolVars, | ||
| 230 | const SubControlVolume& scv) const | ||
| 231 | { | ||
| 232 | // TODO maybe forward to the problem? -> necessary for reaction terms | ||
| 233 | } | ||
| 234 | |||
| 235 | template<class PartialDerivativeMatrices, class T = TypeTag> | ||
| 236 | std::enable_if_t<GetPropType<T, Properties::GridGeometry>::discMethod != DiscretizationMethods::box, void> | ||
| 237 | 19600 | addFluxDerivatives(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 | if constexpr (FVElementGeometry::GridGeometry::discMethod != DiscretizationMethods::cctpfa) | ||
| 246 | DUNE_THROW(Dune::NotImplemented, "Analytic flux differentiation only implemented for tpfa"); | ||
| 247 | |||
| 248 | // advective term: we do the same for all tracer components | ||
| 249 | 19600 | auto rho = [](const VolumeVariables& volVars) | |
| 250 | 19600 | { return useMoles ? volVars.molarDensity() : volVars.density(); }; | |
| 251 | |||
| 252 | // the volume flux | ||
| 253 | 19600 | const auto volFlux = problem.spatialParams().volumeFlux(element, fvGeometry, curElemVolVars, scvf); | |
| 254 | |||
| 255 | // the upwind weight | ||
| 256 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 19598 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
|
19600 | static const Scalar upwindWeight = getParam<Scalar>("Flux.UpwindWeight"); |
| 257 | |||
| 258 | // get the inside and outside volvars | ||
| 259 | 19600 | const auto& insideVolVars = curElemVolVars[scvf.insideScvIdx()]; | |
| 260 |
2/2✓ Branch 1 taken 9800 times.
✓ Branch 2 taken 9800 times.
|
19600 | const auto& outsideVolVars = curElemVolVars[scvf.outsideScvIdx()]; |
| 261 | |||
| 262 |
2/2✓ Branch 0 taken 9800 times.
✓ Branch 1 taken 9800 times.
|
19600 | const Scalar insideWeight = std::signbit(volFlux) ? (1.0 - upwindWeight) : upwindWeight; |
| 263 | 19600 | const Scalar outsideWeight = 1.0 - insideWeight; | |
| 264 | 19600 | const auto advDerivII = volFlux*rho(insideVolVars)*insideWeight; | |
| 265 | 19600 | const auto advDerivIJ = volFlux*rho(outsideVolVars)*outsideWeight; | |
| 266 | |||
| 267 | // diffusive term | ||
| 268 | static constexpr auto referenceSystemFormulation = FluxVariables::MolecularDiffusionType::referenceSystemFormulation(); | ||
| 269 | 19600 | const auto& fluxCache = elemFluxVarsCache[scvf]; | |
| 270 | 19600 | const Scalar rhoInside = massOrMolarDensity(insideVolVars, referenceSystemFormulation, phaseIdx); | |
| 271 | 19600 | const Scalar rhoOutside = massOrMolarDensity(outsideVolVars, referenceSystemFormulation, phaseIdx); | |
| 272 | 19600 | const Scalar massOrMolarDensity = 0.5*(rhoInside + rhoOutside); | |
| 273 | |||
| 274 |
2/2✓ Branch 0 taken 39200 times.
✓ Branch 1 taken 19600 times.
|
58800 | for (int compIdx = 0; compIdx < numComponents; ++compIdx) |
| 275 | { | ||
| 276 | // diffusive term | ||
| 277 | 39200 | Scalar diffDeriv = 0.0; | |
| 278 | if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged) | ||
| 279 | { | ||
| 280 | 39200 | diffDeriv = useMoles ? massOrMolarDensity*fluxCache.diffusionTij(phaseIdx, compIdx)/FluidSystem::molarMass(compIdx) | |
| 281 | 39200 | : massOrMolarDensity*fluxCache.diffusionTij(phaseIdx, compIdx); | |
| 282 | } | ||
| 283 | else if (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged) | ||
| 284 | { | ||
| 285 | diffDeriv = useMoles ? massOrMolarDensity*fluxCache.diffusionTij(phaseIdx, compIdx) | ||
| 286 | : massOrMolarDensity*fluxCache.diffusionTij(phaseIdx, compIdx)*FluidSystem::molarMass(compIdx); | ||
| 287 | } | ||
| 288 | else | ||
| 289 | 39200 | DUNE_THROW(Dune::NotImplemented, "other reference systems than mass and molar averaged are not implemented"); | |
| 290 | |||
| 291 |
1/2✓ Branch 1 taken 39200 times.
✗ Branch 2 not taken.
|
39200 | derivativeMatrices[scvf.insideScvIdx()][compIdx][compIdx] += (advDerivII + diffDeriv); |
| 292 |
1/2✓ Branch 0 taken 39200 times.
✗ Branch 1 not taken.
|
39200 | if (!scvf.boundary()) |
| 293 | 39200 | derivativeMatrices[scvf.outsideScvIdx()][compIdx][compIdx] += (advDerivIJ - diffDeriv); | |
| 294 | } | ||
| 295 | 19600 | } | |
| 296 | |||
| 297 | template<class JacobianMatrix, class T = TypeTag> | ||
| 298 | std::enable_if_t<GetPropType<T, Properties::GridGeometry>::discMethod == DiscretizationMethods::box, void> | ||
| 299 | 20000 | addFluxDerivatives(JacobianMatrix& A, | |
| 300 | const Problem& problem, | ||
| 301 | const Element& element, | ||
| 302 | const FVElementGeometry& fvGeometry, | ||
| 303 | const ElementVolumeVariables& curElemVolVars, | ||
| 304 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 305 | const SubControlVolumeFace& scvf) const | ||
| 306 | { | ||
| 307 | |||
| 308 | // advective term: we do the same for all tracer components | ||
| 309 | 20000 | auto rho = [](const VolumeVariables& volVars) | |
| 310 | 20000 | { return useMoles ? volVars.molarDensity() : volVars.density(); }; | |
| 311 | |||
| 312 | // the volume flux | ||
| 313 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 19996 times.
|
20000 | const auto volFlux = problem.spatialParams().volumeFlux(element, fvGeometry, curElemVolVars, scvf); |
| 314 | |||
| 315 | // the upwind weight | ||
| 316 |
5/6✓ Branch 0 taken 4 times.
✓ Branch 1 taken 19996 times.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
|
20000 | static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight"); |
| 317 | |||
| 318 | // get the inside and outside volvars | ||
| 319 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 20000 times.
|
20000 | const auto& insideVolVars = curElemVolVars[scvf.insideScvIdx()]; |
| 320 |
2/2✓ Branch 0 taken 10000 times.
✓ Branch 1 taken 10000 times.
|
20000 | const auto& outsideVolVars = curElemVolVars[scvf.outsideScvIdx()]; |
| 321 | |||
| 322 |
2/2✓ Branch 0 taken 10000 times.
✓ Branch 1 taken 10000 times.
|
20000 | const auto insideWeight = std::signbit(volFlux) ? (1.0 - upwindWeight) : upwindWeight; |
| 323 | 20000 | const auto outsideWeight = 1.0 - insideWeight; | |
| 324 | 20000 | const auto advDerivII = volFlux*rho(insideVolVars)*insideWeight; | |
| 325 | 20000 | const auto advDerivIJ = volFlux*rho(outsideVolVars)*outsideWeight; | |
| 326 | |||
| 327 | // diffusive term | ||
| 328 | static constexpr auto referenceSystemFormulation = FluxVariables::MolecularDiffusionType::referenceSystemFormulation(); | ||
| 329 | using DiffusionType = GetPropType<T, Properties::MolecularDiffusionType>; | ||
| 330 | 20000 | const auto ti = DiffusionType::calculateTransmissibilities(problem, | |
| 331 | element, | ||
| 332 | fvGeometry, | ||
| 333 | curElemVolVars, | ||
| 334 | scvf, | ||
| 335 | 20000 | elemFluxVarsCache[scvf], | |
| 336 | phaseIdx); | ||
| 337 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 20000 times.
|
20000 | const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); |
| 338 | 20000 | const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx()); | |
| 339 | |||
| 340 |
2/2✓ Branch 0 taken 40000 times.
✓ Branch 1 taken 20000 times.
|
60000 | for (int compIdx = 0; compIdx < numComponents; ++compIdx) |
| 341 | { | ||
| 342 |
2/2✓ Branch 0 taken 160000 times.
✓ Branch 1 taken 40000 times.
|
200000 | for (const auto& scv : scvs(fvGeometry)) |
| 343 | { | ||
| 344 | // diffusive term | ||
| 345 | 160000 | auto diffDeriv = 0.0; | |
| 346 | if (referenceSystemFormulation == ReferenceSystemFormulation::massAveraged) | ||
| 347 | 160000 | diffDeriv += useMoles ? ti[compIdx][scv.indexInElement()]/FluidSystem::molarMass(compIdx) | |
| 348 |
1/2✓ Branch 1 taken 160000 times.
✗ Branch 2 not taken.
|
160000 | : ti[compIdx][scv.indexInElement()]; |
| 349 | else if (referenceSystemFormulation == ReferenceSystemFormulation::molarAveraged) | ||
| 350 | diffDeriv += useMoles ? ti[compIdx][scv.indexInElement()] | ||
| 351 | : ti[compIdx][scv.indexInElement()]*FluidSystem::molarMass(compIdx); | ||
| 352 | else | ||
| 353 |
1/2✓ Branch 1 taken 160000 times.
✗ Branch 2 not taken.
|
160000 | DUNE_THROW(Dune::NotImplemented, "other reference systems than mass and molar averaged are not implemented"); |
| 354 |
2/4✓ Branch 1 taken 160000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 160000 times.
✗ Branch 5 not taken.
|
160000 | A[insideScv.dofIndex()][scv.dofIndex()][compIdx][compIdx] += diffDeriv; |
| 355 |
1/2✓ Branch 1 taken 160000 times.
✗ Branch 2 not taken.
|
160000 | A[outsideScv.dofIndex()][scv.dofIndex()][compIdx][compIdx] -= diffDeriv; |
| 356 | } | ||
| 357 | |||
| 358 |
2/4✓ Branch 1 taken 40000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 40000 times.
✗ Branch 5 not taken.
|
40000 | A[insideScv.dofIndex()][insideScv.dofIndex()][compIdx][compIdx] += advDerivII; |
| 359 |
2/4✓ Branch 1 taken 40000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 40000 times.
✗ Branch 5 not taken.
|
40000 | A[insideScv.dofIndex()][outsideScv.dofIndex()][compIdx][compIdx] += advDerivIJ; |
| 360 |
2/4✓ Branch 1 taken 40000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 40000 times.
✗ Branch 5 not taken.
|
40000 | A[outsideScv.dofIndex()][outsideScv.dofIndex()][compIdx][compIdx] -= advDerivII; |
| 361 |
1/2✓ Branch 1 taken 40000 times.
✗ Branch 2 not taken.
|
40000 | A[outsideScv.dofIndex()][insideScv.dofIndex()][compIdx][compIdx] -= advDerivIJ; |
| 362 | } | ||
| 363 | 20000 | } | |
| 364 | |||
| 365 | template<class PartialDerivativeMatrices> | ||
| 366 | ✗ | void addCCDirichletFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices, | |
| 367 | const Problem& problem, | ||
| 368 | const Element& element, | ||
| 369 | const FVElementGeometry& fvGeometry, | ||
| 370 | const ElementVolumeVariables& curElemVolVars, | ||
| 371 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 372 | const SubControlVolumeFace& scvf) const | ||
| 373 | { | ||
| 374 | // do the same as for inner facets | ||
| 375 | ✗ | addFluxDerivatives(derivativeMatrices, problem, element, fvGeometry, | |
| 376 | curElemVolVars, elemFluxVarsCache, scvf); | ||
| 377 | ✗ | } | |
| 378 | |||
| 379 | template<class PartialDerivativeMatrices> | ||
| 380 | void addRobinFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices, | ||
| 381 | const Problem& problem, | ||
| 382 | const Element& element, | ||
| 383 | const FVElementGeometry& fvGeometry, | ||
| 384 | const ElementVolumeVariables& curElemVolVars, | ||
| 385 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 386 | const SubControlVolumeFace& scvf) const | ||
| 387 | { | ||
| 388 | // TODO maybe forward to the problem? | ||
| 389 | } | ||
| 390 | }; | ||
| 391 | |||
| 392 | } // end namespace Dumux | ||
| 393 | |||
| 394 | #endif | ||
| 395 |