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