| 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 Assembly | ||
| 10 | * \ingroup StaggeredDiscretization | ||
| 11 | * \brief Calculates the element-wise residual for the staggered FV scheme | ||
| 12 | */ | ||
| 13 | #ifndef DUMUX_STAGGERED_LOCAL_RESIDUAL_HH | ||
| 14 | #define DUMUX_STAGGERED_LOCAL_RESIDUAL_HH | ||
| 15 | |||
| 16 | #include <dumux/common/timeloop.hh> | ||
| 17 | #include <dumux/common/properties.hh> | ||
| 18 | #include <dumux/discretization/extrusion.hh> | ||
| 19 | |||
| 20 | namespace Dumux { | ||
| 21 | |||
| 22 | /*! | ||
| 23 | * \ingroup Assembly | ||
| 24 | * \ingroup StaggeredDiscretization | ||
| 25 | * \brief Calculates the element-wise residual for the staggered FV scheme | ||
| 26 | */ | ||
| 27 | template<class TypeTag> | ||
| 28 | class StaggeredLocalResidual | ||
| 29 | { | ||
| 30 | using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView; | ||
| 31 | |||
| 32 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 33 | using Implementation = GetPropType<TypeTag, Properties::LocalResidual>; | ||
| 34 | using Problem = GetPropType<TypeTag, Properties::Problem>; | ||
| 35 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 36 | using ElementBoundaryTypes = GetPropType<TypeTag, Properties::ElementBoundaryTypes>; | ||
| 37 | using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView; | ||
| 38 | using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView; | ||
| 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 FaceSubControlVolume = typename GridGeometry::Traits::FaceSubControlVolume; | ||
| 44 | using Extrusion = Extrusion_t<GridGeometry>; | ||
| 45 | using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>; | ||
| 46 | |||
| 47 | using CellCenterResidual = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>; | ||
| 48 | using FaceResidual = GetPropType<TypeTag, Properties::FacePrimaryVariables>; | ||
| 49 | using ElementFaceVariables = typename GetPropType<TypeTag, Properties::GridFaceVariables>::LocalView; | ||
| 50 | |||
| 51 | using TimeLoop = TimeLoopBase<Scalar>; | ||
| 52 | |||
| 53 | public: | ||
| 54 | using CellCenterResidualValue = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>; | ||
| 55 | using FaceResidualValue = GetPropType<TypeTag, Properties::FacePrimaryVariables>; | ||
| 56 | using ElementResidualVector = CellCenterResidualValue; | ||
| 57 | |||
| 58 | //! the constructor | ||
| 59 | 2628 | StaggeredLocalResidual(const Problem* problem, | |
| 60 | const TimeLoop* timeLoop = nullptr) | ||
| 61 | : problem_(problem) | ||
| 62 |
2/2✓ Branch 0 taken 2578 times.
✓ Branch 1 taken 50 times.
|
2628 | , timeLoop_(timeLoop) |
| 63 | {} | ||
| 64 | |||
| 65 | //! Convenience function to evaluate the flux and source terms for the cell center residual | ||
| 66 | 41373356 | CellCenterResidualValue evalFluxAndSourceForCellCenter(const Element& element, | |
| 67 | const FVElementGeometry& fvGeometry, | ||
| 68 | const ElementVolumeVariables& elemVolVars, | ||
| 69 | const ElementFaceVariables& elemFaceVars, | ||
| 70 | const ElementBoundaryTypes& bcTypes, | ||
| 71 | const ElementFluxVariablesCache& elemFluxVarsCache) const | ||
| 72 | { | ||
| 73 | 41373356 | CellCenterResidualValue residual(0.0); | |
| 74 | |||
| 75 | // evaluate the source term | ||
| 76 |
4/4✓ Branch 0 taken 41373356 times.
✓ Branch 1 taken 41373356 times.
✓ Branch 2 taken 41373356 times.
✓ Branch 3 taken 41373356 times.
|
124120068 | for (auto&& scv : scvs(fvGeometry)) |
| 77 | 41373356 | asImp().evalSourceForCellCenter(residual, this->problem(), element, fvGeometry, elemVolVars, elemFaceVars, scv); | |
| 78 | |||
| 79 | // evaluate the flux term | ||
| 80 |
4/4✓ Branch 0 taken 165493424 times.
✓ Branch 1 taken 41373356 times.
✓ Branch 2 taken 165493424 times.
✓ Branch 3 taken 41373356 times.
|
248240136 | for (auto&& scvf : scvfs(fvGeometry)) |
| 81 | 165493424 | asImp().evalFluxForCellCenter(residual, this->problem(), element, fvGeometry, elemVolVars, elemFaceVars, bcTypes, elemFluxVarsCache, scvf); | |
| 82 | |||
| 83 | 41373356 | return residual; | |
| 84 | } | ||
| 85 | |||
| 86 | //! Evaluate the flux terms for a cell center residual | ||
| 87 | 165493424 | void evalFluxForCellCenter(CellCenterResidualValue& residual, | |
| 88 | const Problem& problem, | ||
| 89 | const Element& element, | ||
| 90 | const FVElementGeometry& fvGeometry, | ||
| 91 | const ElementVolumeVariables& elemVolVars, | ||
| 92 | const ElementFaceVariables& elemFaceVars, | ||
| 93 | const ElementBoundaryTypes& elemBcTypes, | ||
| 94 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 95 | const SubControlVolumeFace& scvf) const | ||
| 96 | { | ||
| 97 |
2/2✓ Branch 0 taken 156137274 times.
✓ Branch 1 taken 9356150 times.
|
165493424 | if (!scvf.boundary()) |
| 98 | 297551780 | residual += asImp_().computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache); | |
| 99 | else | ||
| 100 | 18424796 | residual += asImp_().computeBoundaryFluxForCellCenter(problem, element, fvGeometry, scvf, elemVolVars, elemFaceVars, elemBcTypes, elemFluxVarsCache); | |
| 101 | 165493424 | } | |
| 102 | |||
| 103 | //! Evaluate the source terms for a cell center residual | ||
| 104 | 37620788 | void evalSourceForCellCenter(CellCenterResidualValue& residual, | |
| 105 | const Problem& problem, | ||
| 106 | const Element& element, | ||
| 107 | const FVElementGeometry& fvGeometry, | ||
| 108 | const ElementVolumeVariables& curElemVolVars, | ||
| 109 | const ElementFaceVariables& curElemFaceVars, | ||
| 110 | const SubControlVolume& scv) const | ||
| 111 | { | ||
| 112 | 37620788 | const auto curExtrusionFactor = curElemVolVars[scv].extrusionFactor(); | |
| 113 | |||
| 114 | // subtract the source term from the local rate | ||
| 115 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
37620788 | auto source = asImp_().computeSourceForCellCenter(problem, element, fvGeometry, curElemVolVars, curElemFaceVars, scv); |
| 116 | 75241576 | source *= Extrusion::volume(fvGeometry, scv)*curExtrusionFactor; | |
| 117 | 37620788 | residual -= source; | |
| 118 | 37620788 | } | |
| 119 | |||
| 120 | //! Evaluate the storage terms for a cell center residual | ||
| 121 | 38399336 | CellCenterResidualValue evalStorageForCellCenter(const Element &element, | |
| 122 | const FVElementGeometry& fvGeometry, | ||
| 123 | const ElementVolumeVariables& prevElemVolVars, | ||
| 124 | const ElementVolumeVariables& curElemVolVars) const | ||
| 125 | { | ||
| 126 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 38399336 times.
|
38399336 | assert(timeLoop_ && "no time loop set for storage term evaluation"); |
| 127 | 38399336 | CellCenterResidualValue storage(0.0); | |
| 128 | |||
| 129 |
4/4✓ Branch 0 taken 38399336 times.
✓ Branch 1 taken 38399336 times.
✓ Branch 2 taken 38399336 times.
✓ Branch 3 taken 38399336 times.
|
115198008 | for (auto&& scv : scvs(fvGeometry)) |
| 130 | 38399336 | asImp().evalStorageForCellCenter(storage, problem(), element, fvGeometry, prevElemVolVars, curElemVolVars, scv); | |
| 131 | |||
| 132 | 38399336 | return storage; | |
| 133 | } | ||
| 134 | |||
| 135 | //! Evaluate the storage terms for a cell center residual | ||
| 136 | 38399336 | void evalStorageForCellCenter(CellCenterResidualValue& residual, | |
| 137 | const Problem& problem, | ||
| 138 | const Element &element, | ||
| 139 | const FVElementGeometry& fvGeometry, | ||
| 140 | const ElementVolumeVariables& prevElemVolVars, | ||
| 141 | const ElementVolumeVariables& curElemVolVars, | ||
| 142 | const SubControlVolume& scv) const | ||
| 143 | { | ||
| 144 | 38399336 | CellCenterResidualValue storage(0.0); | |
| 145 | 38399336 | const auto& curVolVars = curElemVolVars[scv]; | |
| 146 | 38399336 | const auto& prevVolVars = prevElemVolVars[scv]; | |
| 147 | |||
| 148 | // mass balance within the element. this is the | ||
| 149 | // \f$\frac{m}{\partial t}\f$ term if using implicit | ||
| 150 | // euler as time discretization. | ||
| 151 | |||
| 152 | // We might need a more explicit way for | ||
| 153 | // doing the time discretization... | ||
| 154 | 38399336 | auto prevCCStorage = asImp_().computeStorageForCellCenter(problem, scv, prevVolVars); | |
| 155 | 38399336 | auto curCCStorage = asImp_().computeStorageForCellCenter(problem, scv, curVolVars); | |
| 156 | |||
| 157 | 38399336 | prevCCStorage *= prevVolVars.extrusionFactor(); | |
| 158 | 38399336 | curCCStorage *= curVolVars.extrusionFactor(); | |
| 159 | |||
| 160 | 38399336 | storage = std::move(curCCStorage); | |
| 161 | 38399336 | storage -= std::move(prevCCStorage); | |
| 162 | 76798672 | storage *= Extrusion::volume(fvGeometry, scv); | |
| 163 | 38399336 | storage /= timeLoop_->timeStepSize(); | |
| 164 | |||
| 165 | 38399336 | residual += storage; | |
| 166 | 38399336 | } | |
| 167 | |||
| 168 | //! for compatibility with FVLocalAssemblerBase | ||
| 169 | template<class... Args> | ||
| 170 | CellCenterResidualValue evalFluxAndSource(Args&&... args) const | ||
| 171 | { | ||
| 172 | return CellCenterResidualValue(0.0); | ||
| 173 | } | ||
| 174 | |||
| 175 | //! for compatibility with FVLocalAssemblerBase | ||
| 176 | template<class... Args> | ||
| 177 | CellCenterResidualValue evalStorage(Args&&... args) const | ||
| 178 | { | ||
| 179 | return CellCenterResidualValue(0.0); | ||
| 180 | } | ||
| 181 | |||
| 182 | /*! | ||
| 183 | * \name User interface | ||
| 184 | * \note The following methods are usually expensive to evaluate | ||
| 185 | * They are useful for outputting residual information. | ||
| 186 | */ | ||
| 187 | // \{ | ||
| 188 | |||
| 189 | //! Convenience function to evaluate the flux and source terms for the face residual | ||
| 190 | 150537748 | FaceResidualValue evalFluxAndSourceForFace(const Element& element, | |
| 191 | const FVElementGeometry& fvGeometry, | ||
| 192 | const ElementVolumeVariables& elemVolVars, | ||
| 193 | const ElementFaceVariables& elemFaceVars, | ||
| 194 | const ElementBoundaryTypes& bcTypes, | ||
| 195 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 196 | const SubControlVolumeFace& scvf) const | ||
| 197 | { | ||
| 198 | 150537748 | FaceResidualValue residual(0.0); | |
| 199 | 150537748 | asImp().evalSourceForFace(residual, this->problem(), element, fvGeometry, elemVolVars, elemFaceVars, scvf); | |
| 200 | 150537748 | asImp().evalFluxForFace(residual, this->problem(), element, fvGeometry, elemVolVars, elemFaceVars, bcTypes, elemFluxVarsCache, scvf); | |
| 201 | |||
| 202 | 150537748 | return residual; | |
| 203 | } | ||
| 204 | |||
| 205 | //! Evaluate the flux terms for a face residual | ||
| 206 | 150537748 | void evalFluxForFace(FaceResidualValue& residual, | |
| 207 | const Problem& problem, | ||
| 208 | const Element& element, | ||
| 209 | const FVElementGeometry& fvGeometry, | ||
| 210 | const ElementVolumeVariables& elemVolVars, | ||
| 211 | const ElementFaceVariables& elemFaceVars, | ||
| 212 | const ElementBoundaryTypes& elemBcTypes, | ||
| 213 | const ElementFluxVariablesCache& elemFluxVarsCache, | ||
| 214 | const SubControlVolumeFace& scvf) const | ||
| 215 | { | ||
| 216 |
2/2✓ Branch 0 taken 141926782 times.
✓ Branch 1 taken 8610966 times.
|
150537748 | if (!scvf.boundary()) |
| 217 | 141926782 | residual += asImp_().computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache); | |
| 218 | else | ||
| 219 | 8610966 | residual += asImp_().computeBoundaryFluxForFace(problem, element, fvGeometry, scvf, elemVolVars, elemFaceVars, elemBcTypes, elemFluxVarsCache); | |
| 220 | 150537748 | } | |
| 221 | |||
| 222 | //! Evaluate the source terms for a face residual | ||
| 223 | 150537748 | void evalSourceForFace(FaceResidualValue& residual, | |
| 224 | const Problem& problem, | ||
| 225 | const Element& element, | ||
| 226 | const FVElementGeometry& fvGeometry, | ||
| 227 | const ElementVolumeVariables& elemVolVars, | ||
| 228 | const ElementFaceVariables& elemFaceVars, | ||
| 229 | const SubControlVolumeFace& scvf) const | ||
| 230 | { | ||
| 231 | // the source term: | ||
| 232 | 150537748 | auto source = asImp_().computeSourceForFace(problem, element, fvGeometry, scvf, elemVolVars, elemFaceVars); | |
| 233 | 301075496 | const auto& scv = fvGeometry.scv(scvf.insideScvIdx()); | |
| 234 | 150537748 | const auto extrusionFactor = elemVolVars[scv].extrusionFactor(); | |
| 235 | |||
| 236 | // construct staggered scv (half of the element) | ||
| 237 | 602150992 | auto faceScvCenter = scvf.center() + scv.center(); | |
| 238 | 301075496 | faceScvCenter *= 0.5; | |
| 239 | 301075496 | FaceSubControlVolume faceScv(faceScvCenter, 0.5*scv.volume()); | |
| 240 | |||
| 241 | 301075496 | source *= Extrusion::volume(fvGeometry, faceScv)*extrusionFactor; | |
| 242 | 150537748 | residual -= source; | |
| 243 | 150537748 | } | |
| 244 | |||
| 245 | //! Evaluate the storage terms for a face residual | ||
| 246 | 135806428 | FaceResidualValue evalStorageForFace(const Element& element, | |
| 247 | const FVElementGeometry& fvGeometry, | ||
| 248 | const ElementVolumeVariables& prevElemVolVars, | ||
| 249 | const ElementVolumeVariables& curElemVolVars, | ||
| 250 | const ElementFaceVariables& prevElemFaceVars, | ||
| 251 | const ElementFaceVariables& curElemFaceVars, | ||
| 252 | const SubControlVolumeFace& scvf) const | ||
| 253 | { | ||
| 254 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 135806428 times.
|
135806428 | assert(timeLoop_ && "no time loop set for storage term evaluation"); |
| 255 | 135806428 | FaceResidualValue storage(0.0); | |
| 256 | 135806428 | asImp().evalStorageForFace(storage, problem(), element, fvGeometry, prevElemVolVars, curElemVolVars, prevElemFaceVars, curElemFaceVars, scvf); | |
| 257 | 135806428 | return storage; | |
| 258 | } | ||
| 259 | |||
| 260 | //! Evaluate the storage terms for a face residual | ||
| 261 | 135806428 | void evalStorageForFace(FaceResidualValue& residual, | |
| 262 | const Problem& problem, | ||
| 263 | const Element& element, | ||
| 264 | const FVElementGeometry& fvGeometry, | ||
| 265 | const ElementVolumeVariables& prevElemVolVars, | ||
| 266 | const ElementVolumeVariables& curElemVolVars, | ||
| 267 | const ElementFaceVariables& prevElemFaceVars, | ||
| 268 | const ElementFaceVariables& curElemFaceVars, | ||
| 269 | const SubControlVolumeFace& scvf) const | ||
| 270 | { | ||
| 271 | 271612856 | const auto& scv = fvGeometry.scv(scvf.insideScvIdx()); | |
| 272 | |||
| 273 | 135806428 | auto storage = asImp_().computeStorageForFace(problem, scvf, curElemVolVars[scv], curElemFaceVars); | |
| 274 | 135806428 | storage -= asImp_().computeStorageForFace(problem, scvf, prevElemVolVars[scv], prevElemFaceVars); | |
| 275 | |||
| 276 | 135806428 | const auto extrusionFactor = curElemVolVars[scv].extrusionFactor(); | |
| 277 | |||
| 278 | // construct staggered scv (half of the element) | ||
| 279 | 407419284 | auto faceScvCenter = scvf.center() + scv.center(); | |
| 280 | 135806428 | faceScvCenter *= 0.5; | |
| 281 | 135806428 | FaceSubControlVolume faceScv(faceScvCenter, 0.5*scv.volume()); | |
| 282 | |||
| 283 | 271612856 | storage *= Extrusion::volume(fvGeometry, faceScv)*extrusionFactor; | |
| 284 | 135806428 | storage /= timeLoop_->timeStepSize(); | |
| 285 | |||
| 286 | 135806428 | residual += storage; | |
| 287 | 135806428 | } | |
| 288 | |||
| 289 | //! If no solution has been set, we treat the problem as stationary. | ||
| 290 | bool isStationary() const | ||
| 291 | { return !timeLoop_; } | ||
| 292 | |||
| 293 | //! the problem | ||
| 294 | ✗ | const Problem& problem() const | |
| 295 | ✗ | { return *problem_; } | |
| 296 | |||
| 297 | protected: | ||
| 298 | |||
| 299 | Implementation &asImp_() | ||
| 300 | { return *static_cast<Implementation*>(this); } | ||
| 301 | |||
| 302 | const Implementation &asImp_() const | ||
| 303 | { return *static_cast<const Implementation*>(this); } | ||
| 304 | |||
| 305 | |||
| 306 | TimeLoop& timeLoop() | ||
| 307 | { return *timeLoop_; } | ||
| 308 | |||
| 309 | const TimeLoop& timeLoop() const | ||
| 310 | { return *timeLoop_; } | ||
| 311 | |||
| 312 | Implementation &asImp() | ||
| 313 | { return *static_cast<Implementation*>(this); } | ||
| 314 | |||
| 315 | const Implementation &asImp() const | ||
| 316 | { return *static_cast<const Implementation*>(this); } | ||
| 317 | |||
| 318 | private: | ||
| 319 | const Problem* problem_; //!< the problem we are assembling this residual for | ||
| 320 | const TimeLoop* timeLoop_; | ||
| 321 | |||
| 322 | }; | ||
| 323 | |||
| 324 | } // end namespace Dumux | ||
| 325 | |||
| 326 | #endif | ||
| 327 |