| 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 | * \ingroup MultiDomain | ||
| 12 | * \brief A multidomain assembler for Jacobian and residual contribution per element (staggered method) | ||
| 13 | */ | ||
| 14 | #ifndef DUMUX_MULTIDOMAIN_STAGGERED_LOCAL_ASSEMBLER_HH | ||
| 15 | #define DUMUX_MULTIDOMAIN_STAGGERED_LOCAL_ASSEMBLER_HH | ||
| 16 | |||
| 17 | #include <dune/common/reservedvector.hh> | ||
| 18 | #include <dune/common/indices.hh> | ||
| 19 | #include <dune/common/hybridutilities.hh> | ||
| 20 | #include <dune/grid/common/gridenums.hh> // for GhostEntity | ||
| 21 | |||
| 22 | #include <dumux/common/reservedblockvector.hh> | ||
| 23 | #include <dumux/common/properties.hh> | ||
| 24 | #include <dumux/common/parameters.hh> | ||
| 25 | #include <dumux/common/numericdifferentiation.hh> | ||
| 26 | #include <dumux/common/typetraits/utility.hh> | ||
| 27 | #include <dumux/assembly/diffmethod.hh> | ||
| 28 | #include <dumux/assembly/fvlocalassemblerbase.hh> | ||
| 29 | #include <dumux/discretization/staggered/elementsolution.hh> | ||
| 30 | |||
| 31 | namespace Dumux { | ||
| 32 | |||
| 33 | /*! | ||
| 34 | * \ingroup Assembly | ||
| 35 | * \ingroup StaggeredDiscretization | ||
| 36 | * \ingroup MultiDomain | ||
| 37 | * \brief A base class for all multidomain local assemblers (staggered) | ||
| 38 | * \tparam id the id of the sub domain | ||
| 39 | * \tparam TypeTag the TypeTag | ||
| 40 | * \tparam Assembler the assembler type | ||
| 41 | * \tparam Implementation the actual assembler implementation | ||
| 42 | * \tparam implicit whether the assembly is explicit or implicit in time | ||
| 43 | */ | ||
| 44 | template<std::size_t id, class TypeTag, class Assembler, class Implementation, bool isImplicit = true> | ||
| 45 | 5441776 | class SubDomainStaggeredLocalAssemblerBase : public FVLocalAssemblerBase<TypeTag, Assembler,Implementation, isImplicit> | |
| 46 | { | ||
| 47 | using ParentType = FVLocalAssemblerBase<TypeTag, Assembler,Implementation, isImplicit>; | ||
| 48 | |||
| 49 | using Problem = GetPropType<TypeTag, Properties::Problem>; | ||
| 50 | using SolutionVector = typename Assembler::SolutionVector; | ||
| 51 | |||
| 52 | using GridVariables = GetPropType<TypeTag, Properties::GridVariables>; | ||
| 53 | using GridVolumeVariables = typename GridVariables::GridVolumeVariables; | ||
| 54 | using ElementVolumeVariables = typename GridVolumeVariables::LocalView; | ||
| 55 | using Scalar = typename GridVariables::Scalar; | ||
| 56 | |||
| 57 | using ElementFaceVariables = typename GetPropType<TypeTag, Properties::GridFaceVariables>::LocalView; | ||
| 58 | using CellCenterResidualValue = typename ParentType::LocalResidual::CellCenterResidualValue; | ||
| 59 | using FaceResidualValue = typename ParentType::LocalResidual::FaceResidualValue; | ||
| 60 | |||
| 61 | using GridGeometry = typename GridVariables::GridGeometry; | ||
| 62 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 63 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 64 | using GridView = typename GridGeometry::GridView; | ||
| 65 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 66 | |||
| 67 | using CouplingManager = typename Assembler::CouplingManager; | ||
| 68 | |||
| 69 | static constexpr auto numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq(); | ||
| 70 | |||
| 71 | public: | ||
| 72 | static constexpr auto domainId = typename Dune::index_constant<id>(); | ||
| 73 | static constexpr auto cellCenterId = GridGeometry::cellCenterIdx(); | ||
| 74 | static constexpr auto faceId = GridGeometry::faceIdx(); | ||
| 75 | |||
| 76 | static constexpr auto numEqCellCenter = CellCenterResidualValue::dimension; | ||
| 77 | static constexpr auto faceOffset = numEqCellCenter; | ||
| 78 | |||
| 79 | using ParentType::ParentType; | ||
| 80 | |||
| 81 | 5441776 | explicit SubDomainStaggeredLocalAssemblerBase(const Assembler& assembler, | |
| 82 | const Element& element, | ||
| 83 | const SolutionVector& curSol, | ||
| 84 | CouplingManager& couplingManager) | ||
| 85 | : ParentType(assembler, | ||
| 86 | element, | ||
| 87 | curSol, | ||
| 88 | localView(assembler.gridGeometry(domainId)), | ||
| 89 |
1/2✓ Branch 1 taken 2720888 times.
✗ Branch 2 not taken.
|
5441776 | localView(assembler.gridVariables(domainId).curGridVolVars()), |
| 90 |
1/2✓ Branch 1 taken 2720888 times.
✗ Branch 2 not taken.
|
5441776 | localView(assembler.gridVariables(domainId).prevGridVolVars()), |
| 91 |
3/6✓ Branch 1 taken 2720888 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2720888 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2720888 times.
✗ Branch 8 not taken.
|
16325328 | localView(assembler.gridVariables(domainId).gridFluxVarsCache()), |
| 92 | assembler.localResidual(domainId), | ||
| 93 | 10883552 | (element.partitionType() == Dune::GhostEntity)) | |
| 94 | 10883552 | , curElemFaceVars_(localView(assembler.gridVariables(domainId).curGridFaceVars())) | |
| 95 | 16325328 | , prevElemFaceVars_(localView(assembler.gridVariables(domainId).prevGridFaceVars())) | |
| 96 |
8/16✓ Branch 2 taken 2720888 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2720888 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 2720888 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 2720888 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 2720888 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 2720888 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 2720888 times.
✗ Branch 21 not taken.
✓ Branch 23 taken 2720888 times.
✗ Branch 24 not taken.
|
5441776 | , couplingManager_(couplingManager) |
| 97 | 5441776 | {} | |
| 98 | |||
| 99 | /*! | ||
| 100 | * \brief Computes the derivatives with respect to the given element and adds them | ||
| 101 | * to the global matrix. The element residual is written into the right hand side. | ||
| 102 | */ | ||
| 103 | template<class JacobianMatrixRow, class SubResidual, class GridVariablesTuple> | ||
| 104 | void assembleJacobianAndResidual(JacobianMatrixRow& jacRow, SubResidual& res, GridVariablesTuple& gridVariables) | ||
| 105 | { | ||
| 106 |
2/4✓ Branch 1 taken 1360444 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1360444 times.
✗ Branch 5 not taken.
|
2720888 | this->asImp_().bindLocalViews(); |
| 107 |
4/8✓ Branch 1 taken 1360444 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1360444 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1360444 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1360444 times.
✗ Branch 11 not taken.
|
5441776 | this->elemBcTypes().update(problem(), this->element(), this->fvGeometry()); |
| 108 | |||
| 109 |
2/4✓ Branch 1 taken 1360444 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1360444 times.
✗ Branch 5 not taken.
|
2720888 | assembleJacobianAndResidualImpl_(domainId, jacRow, res, gridVariables); |
| 110 | } | ||
| 111 | |||
| 112 | /*! | ||
| 113 | * \brief Assemble the residual only | ||
| 114 | */ | ||
| 115 | template<class SubResidual> | ||
| 116 | ✗ | void assembleResidual(SubResidual& res) | |
| 117 | { | ||
| 118 | ✗ | this->asImp_().bindLocalViews(); | |
| 119 | ✗ | this->elemBcTypes().update(problem(), this->element(), this->fvGeometry()); | |
| 120 | |||
| 121 | if constexpr (domainId == cellCenterId) | ||
| 122 | { | ||
| 123 | ✗ | const auto cellCenterGlobalI = problem().gridGeometry().elementMapper().index(this->element()); | |
| 124 | ✗ | res[cellCenterGlobalI] = this->asImp_().assembleCellCenterResidualImpl(); | |
| 125 | } | ||
| 126 | else | ||
| 127 | { | ||
| 128 | ✗ | for (auto&& scvf : scvfs(this->fvGeometry())) | |
| 129 | ✗ | res[scvf.dofIndex()] += this->asImp_().assembleFaceResidualImpl(scvf); | |
| 130 | } | ||
| 131 | ✗ | } | |
| 132 | |||
| 133 | /*! | ||
| 134 | * \brief Convenience function to evaluate the complete local residual for the current element. Automatically chooses the the appropriate | ||
| 135 | * element volume variables. | ||
| 136 | */ | ||
| 137 | 37900788 | CellCenterResidualValue evalLocalResidualForCellCenter() const | |
| 138 | { | ||
| 139 | if (!isImplicit) | ||
| 140 | if (this->assembler().isStationaryProblem()) | ||
| 141 | DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual"); | ||
| 142 | |||
| 143 |
3/6✓ Branch 0 taken 3620 times.
✓ Branch 1 taken 37900788 times.
✓ Branch 2 taken 5200 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
38983836 | if (this->elementIsGhost()) |
| 144 | { | ||
| 145 | ✗ | return CellCenterResidualValue(0.0); | |
| 146 | } | ||
| 147 | |||
| 148 |
3/5✓ Branch 1 taken 3620 times.
✓ Branch 2 taken 540936 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 1327520 times.
✗ Branch 8 not taken.
|
41653356 | return isImplicit ? evalLocalResidualForCellCenter(this->curElemVolVars(), this->curElemFaceVars()) |
| 149 | 37900788 | : evalLocalResidualForCellCenter(this->prevElemVolVars(), this->prevElemFaceVars()); | |
| 150 | } | ||
| 151 | |||
| 152 | /*! | ||
| 153 | * \brief Evaluates the complete local residual for the current cell center. | ||
| 154 | * \param elemVolVars The element volume variables | ||
| 155 | * \param elemFaceVars The element face variables | ||
| 156 | */ | ||
| 157 | 41653356 | CellCenterResidualValue evalLocalResidualForCellCenter(const ElementVolumeVariables& elemVolVars, | |
| 158 | const ElementFaceVariables& elemFaceVars) const | ||
| 159 | { | ||
| 160 | 83306712 | auto residual = evalLocalFluxAndSourceResidualForCellCenter(elemVolVars, elemFaceVars); | |
| 161 | |||
| 162 |
4/4✓ Branch 0 taken 38679336 times.
✓ Branch 1 taken 2974020 times.
✓ Branch 2 taken 38679336 times.
✓ Branch 3 taken 2974020 times.
|
83306712 | if (!this->assembler().isStationaryProblem()) |
| 163 | 43203972 | residual += evalLocalStorageResidualForCellCenter(); | |
| 164 | |||
| 165 | // handle cells with a fixed Dirichlet value | ||
| 166 | 124960068 | const auto cellCenterGlobalI = problem().gridGeometry().elementMapper().index(this->element()); | |
| 167 | 58063212 | const auto& scvI = this->fvGeometry().scv(cellCenterGlobalI); | |
| 168 |
2/2✓ Branch 0 taken 92426480 times.
✓ Branch 1 taken 25243500 times.
|
134079836 | for (int pvIdx = 0; pvIdx < numEqCellCenter; ++pvIdx) |
| 169 | { | ||
| 170 | static constexpr auto offset = numEq - numEqCellCenter; | ||
| 171 |
2/2✓ Branch 1 taken 6596088 times.
✓ Branch 2 taken 85369392 times.
|
92878520 | if (this->problem().isDirichletCell(this->element(), this->fvGeometry(), scvI, pvIdx + offset)) |
| 172 | 26411232 | residual[pvIdx] = this->curSol()[cellCenterId][cellCenterGlobalI][pvIdx] - this->problem().dirichlet(this->element(), scvI)[pvIdx + offset]; | |
| 173 | } | ||
| 174 | |||
| 175 | 41653356 | return residual; | |
| 176 | } | ||
| 177 | |||
| 178 | /*! | ||
| 179 | * \brief Convenience function to evaluate the flux and source terms (i.e, the terms without a time derivative) | ||
| 180 | * of the local residual for the current element. Automatically chooses the the appropriate | ||
| 181 | * element volume and face variables. | ||
| 182 | */ | ||
| 183 | CellCenterResidualValue evalLocalFluxAndSourceResidualForCellCenter() const | ||
| 184 | { | ||
| 185 | return isImplicit ? evalLocalFluxAndSourceResidualForCellCenter(this->curElemVolVars(), this->curElemFaceVars()) | ||
| 186 | : evalLocalFluxAndSourceResidualForCellCenter(this->prevElemVolVars(), this->prevElemFaceVars()); | ||
| 187 | } | ||
| 188 | |||
| 189 | /*! | ||
| 190 | * \brief Evaluates the flux and source terms (i.e, the terms without a time derivative) | ||
| 191 | * of the local residual for the current element. | ||
| 192 | * | ||
| 193 | * \param elemVolVars The element volume variables | ||
| 194 | * \param elemFaceVars The element face variables | ||
| 195 | */ | ||
| 196 | CellCenterResidualValue evalLocalFluxAndSourceResidualForCellCenter(const ElementVolumeVariables& elemVolVars, const ElementFaceVariables& elemFaceVars) const | ||
| 197 | { | ||
| 198 | 166613424 | return this->localResidual().evalFluxAndSourceForCellCenter(this->element(), this->fvGeometry(), elemVolVars, elemFaceVars, this->elemBcTypes(), this->elemFluxVarsCache()); | |
| 199 | } | ||
| 200 | |||
| 201 | /*! | ||
| 202 | * \brief Convenience function to evaluate storage term (i.e, the term with a time derivative) | ||
| 203 | * of the local residual for the current element. Automatically chooses the the appropriate | ||
| 204 | * element volume and face variables. | ||
| 205 | */ | ||
| 206 | CellCenterResidualValue evalLocalStorageResidualForCellCenter() const | ||
| 207 | { | ||
| 208 | 116038008 | return this->localResidual().evalStorageForCellCenter(this->element(), this->fvGeometry(), this->prevElemVolVars(), this->curElemVolVars()); | |
| 209 | } | ||
| 210 | |||
| 211 | /*! | ||
| 212 | * \brief Convenience function to evaluate the local residual for the current face. Automatically chooses the the appropriate | ||
| 213 | * element volume and face variables. | ||
| 214 | * \param scvf The sub control volume face | ||
| 215 | */ | ||
| 216 | FaceResidualValue evalLocalResidualForFace(const SubControlVolumeFace& scvf) const | ||
| 217 | { | ||
| 218 | if (!isImplicit) | ||
| 219 | if (this->assembler().isStationaryProblem()) | ||
| 220 | DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual"); | ||
| 221 | |||
| 222 |
2/6✓ Branch 0 taken 11332 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 22896 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
129791364 | if (this->elementIsGhost()) |
| 223 | { | ||
| 224 | return FaceResidualValue(0.0); | ||
| 225 | } | ||
| 226 | |||
| 227 |
4/8✓ Branch 1 taken 68798596 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3790128 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 6980320 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 1651648 times.
✗ Branch 10 not taken.
|
151435748 | return isImplicit ? evalLocalResidualForFace(scvf, this->curElemVolVars(), this->curElemFaceVars()) |
| 228 | : evalLocalResidualForFace(scvf, this->prevElemVolVars(), this->prevElemFaceVars()); | ||
| 229 | } | ||
| 230 | |||
| 231 | /*! | ||
| 232 | * \brief Evaluates the complete local residual for the current face. | ||
| 233 | * \param scvf The sub control volume face | ||
| 234 | * \param elemVolVars The element volume variables | ||
| 235 | * \param elemFaceVars The element face variables | ||
| 236 | */ | ||
| 237 | 151435748 | FaceResidualValue evalLocalResidualForFace(const SubControlVolumeFace& scvf, | |
| 238 | const ElementVolumeVariables& elemVolVars, | ||
| 239 | const ElementFaceVariables& elemFaceVars) const | ||
| 240 | { | ||
| 241 | 302871496 | auto residual = evalLocalFluxAndSourceResidualForFace(scvf, elemVolVars, elemFaceVars); | |
| 242 | |||
| 243 |
4/4✓ Branch 0 taken 136704428 times.
✓ Branch 1 taken 14731320 times.
✓ Branch 2 taken 136704428 times.
✓ Branch 3 taken 14731320 times.
|
302871496 | if (!this->assembler().isStationaryProblem()) |
| 244 | 273408856 | residual += evalLocalStorageResidualForFace(scvf); | |
| 245 | |||
| 246 | 172838528 | this->localResidual().evalDirichletBoundariesForFace(residual, this->problem(), this->element(), | |
| 247 | this->fvGeometry(), scvf, elemVolVars, elemFaceVars, | ||
| 248 | this->elemBcTypes(), this->elemFluxVarsCache()); | ||
| 249 | |||
| 250 | 151435748 | return residual; | |
| 251 | } | ||
| 252 | |||
| 253 | /*! | ||
| 254 | * \brief Convenience function to evaluate the flux and source terms (i.e, the terms without a time derivative) | ||
| 255 | * of the local residual for the current element. Automatically chooses the the appropriate | ||
| 256 | * element volume and face variables. | ||
| 257 | * \param scvf The sub control volume face | ||
| 258 | */ | ||
| 259 | FaceResidualValue evalLocalFluxAndSourceResidualForFace(const SubControlVolumeFace& scvf) const | ||
| 260 | { | ||
| 261 | return isImplicit ? evalLocalFluxAndSourceResidualForFace(scvf, this->curElemVolVars(), this->curElemFaceVars()) | ||
| 262 | : evalLocalFluxAndSourceResidualForFace(scvf, this->prevElemVolVars(), this->prevElemFaceVars()); | ||
| 263 | } | ||
| 264 | |||
| 265 | /*! | ||
| 266 | * \brief Evaluates the flux and source terms (i.e, the terms without a time derivative) | ||
| 267 | * of the local residual for the current face. | ||
| 268 | * \param scvf The sub control volume face | ||
| 269 | * \param elemVolVars The element volume variables | ||
| 270 | * \param elemFaceVars The element face variables | ||
| 271 | */ | ||
| 272 | FaceResidualValue evalLocalFluxAndSourceResidualForFace(const SubControlVolumeFace& scvf, | ||
| 273 | const ElementVolumeVariables& elemVolVars, | ||
| 274 | const ElementFaceVariables& elemFaceVars) const | ||
| 275 | { | ||
| 276 | 605742992 | return this->localResidual().evalFluxAndSourceForFace(this->element(), this->fvGeometry(), elemVolVars, elemFaceVars, this->elemBcTypes(), this->elemFluxVarsCache(), scvf); | |
| 277 | } | ||
| 278 | |||
| 279 | /*! | ||
| 280 | * \brief Convenience function to evaluate storage term (i.e, the term with a time derivative) | ||
| 281 | * of the local residual for the current face. Automatically chooses the the appropriate | ||
| 282 | * element volume and face variables. | ||
| 283 | * \param scvf The sub control volume face | ||
| 284 | */ | ||
| 285 | FaceResidualValue evalLocalStorageResidualForFace(const SubControlVolumeFace& scvf) const | ||
| 286 | { | ||
| 287 | 273408856 | return this->localResidual().evalStorageForFace(this->element(), this->fvGeometry(), this->prevElemVolVars(), this->curElemVolVars(), this->prevElemFaceVars(), this->curElemFaceVars(), scvf); | |
| 288 | } | ||
| 289 | |||
| 290 | const Problem& problem() const | ||
| 291 |
4/7✓ Branch 5 taken 738608 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 621836 times.
✓ Branch 8 taken 738608 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 621836 times.
✗ Branch 11 not taken.
|
294841520 | { return this->assembler().problem(domainId); } |
| 292 | |||
| 293 | //! The current element volume variables | ||
| 294 | ElementFaceVariables& curElemFaceVars() | ||
| 295 | 50982752 | { return curElemFaceVars_; } | |
| 296 | |||
| 297 | //! The element volume variables of the provious time step | ||
| 298 | ElementFaceVariables& prevElemFaceVars() | ||
| 299 | { return prevElemFaceVars_; } | ||
| 300 | |||
| 301 | //! The current element volume variables | ||
| 302 | const ElementFaceVariables& curElemFaceVars() const | ||
| 303 |
7/12✓ Branch 1 taken 68798596 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 3443776 times.
✓ Branch 6 taken 7521256 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 349972 times.
✓ Branch 10 taken 578048 times.
✗ Branch 11 not taken.
✓ Branch 16 taken 1327520 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1073600 times.
✗ Branch 20 not taken.
|
193089104 | { return curElemFaceVars_; } |
| 304 | |||
| 305 | //! The element volume variables of the provious time step | ||
| 306 | const ElementFaceVariables& prevElemFaceVars() const | ||
| 307 | 136704428 | { return prevElemFaceVars_; } | |
| 308 | |||
| 309 | ✗ | CouplingManager& couplingManager() | |
| 310 | ✗ | { return couplingManager_; } | |
| 311 | |||
| 312 | private: | ||
| 313 | |||
| 314 | //! Assembles the residuals and derivatives for the cell center dofs. | ||
| 315 | template<class JacobianMatrixRow, class SubResidual, class GridVariablesTuple> | ||
| 316 | 1360444 | auto assembleJacobianAndResidualImpl_(Dune::index_constant<cellCenterId>, JacobianMatrixRow& jacRow, SubResidual& res, GridVariablesTuple& gridVariables) | |
| 317 | { | ||
| 318 | 2720888 | auto& gridVariablesI = *std::get<domainId>(gridVariables); | |
| 319 | 5441776 | const auto cellCenterGlobalI = problem().gridGeometry().elementMapper().index(this->element()); | |
| 320 | 2720888 | const auto residual = this->asImp_().assembleCellCenterJacobianAndResidualImpl(jacRow[domainId], gridVariablesI); | |
| 321 | 1360444 | res[cellCenterGlobalI] = residual; | |
| 322 | |||
| 323 | |||
| 324 | // for the coupling blocks | ||
| 325 | using namespace Dune::Hybrid; | ||
| 326 | static constexpr auto otherDomainIds = makeIncompleteIntegerSequence<JacobianMatrixRow::size(), domainId>{}; | ||
| 327 | 1360444 | forEach(otherDomainIds, [&](auto&& domainJ) | |
| 328 | { | ||
| 329 |
0/13✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
1895064 | this->asImp_().assembleJacobianCellCenterCoupling(domainJ, jacRow[domainJ], residual, gridVariablesI); |
| 330 | }); | ||
| 331 | |||
| 332 | // handle cells with a fixed Dirichlet value | ||
| 333 | 2099052 | incorporateDirichletCells_(jacRow); | |
| 334 | 1360444 | } | |
| 335 | |||
| 336 | //! Assembles the residuals and derivatives for the face dofs. | ||
| 337 | template<class JacobianMatrixRow, class SubResidual, class GridVariablesTuple> | ||
| 338 | 1360444 | void assembleJacobianAndResidualImpl_(Dune::index_constant<faceId>, JacobianMatrixRow& jacRow, SubResidual& res, GridVariablesTuple& gridVariables) | |
| 339 | { | ||
| 340 | 2720888 | auto& gridVariablesI = *std::get<domainId>(gridVariables); | |
| 341 |
0/2✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
2720888 | const auto residual = this->asImp_().assembleFaceJacobianAndResidualImpl(jacRow[domainId], gridVariablesI); |
| 342 | |||
| 343 |
4/4✓ Branch 0 taken 5441776 times.
✓ Branch 1 taken 1360444 times.
✓ Branch 2 taken 5441776 times.
✓ Branch 3 taken 1360444 times.
|
8162664 | for(auto&& scvf : scvfs(this->fvGeometry())) |
| 344 | 32650656 | res[scvf.dofIndex()] += residual[scvf.localFaceIdx()]; | |
| 345 | |||
| 346 | // for the coupling blocks | ||
| 347 | using namespace Dune::Hybrid; | ||
| 348 | static constexpr auto otherDomainIds = makeIncompleteIntegerSequence<JacobianMatrixRow::size(), domainId>{}; | ||
| 349 |
2/6✓ Branch 1 taken 1360444 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1360444 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
1360444 | forEach(otherDomainIds, [&](auto&& domainJ) |
| 350 | { | ||
| 351 |
2/12✓ Branch 1 taken 947532 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 947532 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
|
1895064 | this->asImp_().assembleJacobianFaceCoupling(domainJ, jacRow[domainJ], residual, gridVariablesI); |
| 352 | }); | ||
| 353 | 1360444 | } | |
| 354 | |||
| 355 | //! If specified in the problem, a fixed Dirichlet value can be assigned to cell centered unknowns such as pressure | ||
| 356 | template<class JacobianMatrixRow> | ||
| 357 | 621836 | void incorporateDirichletCells_(JacobianMatrixRow& jacRow) | |
| 358 | { | ||
| 359 | 5441776 | const auto cellCenterGlobalI = problem().gridGeometry().elementMapper().index(this->element()); | |
| 360 | |||
| 361 | // overwrite the partial derivative with zero in case a fixed Dirichlet BC is used | ||
| 362 | static constexpr auto offset = numEq - numEqCellCenter; | ||
| 363 |
2/2✓ Branch 0 taken 2122216 times.
✓ Branch 1 taken 621836 times.
|
3482660 | for (int eqIdx = 0; eqIdx < numEqCellCenter; ++eqIdx) |
| 364 | { | ||
| 365 |
2/2✓ Branch 4 taken 155134 times.
✓ Branch 5 taken 1920582 times.
|
8534364 | if (problem().isDirichletCell(this->element(), this->fvGeometry(), this->fvGeometry().scv(cellCenterGlobalI), eqIdx + offset)) |
| 366 | { | ||
| 367 | using namespace Dune::Hybrid; | ||
| 368 | 4333440 | forEach(integralRange(Dune::Hybrid::size(jacRow)), [&, domainId = domainId](auto&& i) | |
| 369 | { | ||
| 370 | 624536 | auto& ccRowI = jacRow[i][cellCenterGlobalI]; | |
| 371 |
8/18✗ Branch 0 not taken.
✓ Branch 1 taken 839581 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 839581 times.
✓ Branch 4 taken 683447 times.
✓ Branch 5 taken 156134 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 780670 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 780670 times.
✓ Branch 10 taken 624536 times.
✓ Branch 11 taken 156134 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
|
2244787 | for (auto col = ccRowI.begin(); col != ccRowI.end(); ++col) |
| 372 | { | ||
| 373 |
4/6✓ Branch 2 taken 1000 times.
✓ Branch 3 taken 2960 times.
✓ Branch 4 taken 1000 times.
✓ Branch 5 taken 2960 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
2615966 | ccRowI[col.index()][eqIdx] = 0.0; |
| 374 | // set the diagonal entry to 1.0 | ||
| 375 |
4/4✓ Branch 0 taken 156134 times.
✓ Branch 1 taken 527313 times.
✓ Branch 2 taken 156134 times.
✓ Branch 3 taken 527313 times.
|
1991430 | if ((i == domainId) && (col.index() == cellCenterGlobalI)) |
| 376 | 156134 | ccRowI[col.index()][eqIdx][eqIdx] = 1.0; | |
| 377 | } | ||
| 378 | }); | ||
| 379 | } | ||
| 380 | } | ||
| 381 | 621836 | } | |
| 382 | |||
| 383 | ElementFaceVariables curElemFaceVars_; | ||
| 384 | ElementFaceVariables prevElemFaceVars_; | ||
| 385 | CouplingManager& couplingManager_; //!< the coupling manager | ||
| 386 | }; | ||
| 387 | |||
| 388 | /*! | ||
| 389 | * \ingroup Assembly | ||
| 390 | * \ingroup StaggeredDiscretization | ||
| 391 | * \ingroup MultiDomain | ||
| 392 | * \brief A base class for all implicit multidomain local assemblers (staggered) | ||
| 393 | * \tparam id the id of the sub domain | ||
| 394 | * \tparam TypeTag the TypeTag | ||
| 395 | * \tparam Assembler the assembler type | ||
| 396 | * \tparam Implementation the actual assembler implementation | ||
| 397 | */ | ||
| 398 | template<std::size_t id, class TypeTag, class Assembler, class Implementation> | ||
| 399 | 5441776 | class SubDomainStaggeredLocalAssemblerImplicitBase : public SubDomainStaggeredLocalAssemblerBase<id, TypeTag, Assembler, Implementation> | |
| 400 | { | ||
| 401 | using ParentType = SubDomainStaggeredLocalAssemblerBase<id, TypeTag, Assembler, Implementation>; | ||
| 402 | static constexpr auto domainId = Dune::index_constant<id>(); | ||
| 403 | public: | ||
| 404 | 2720888 | using ParentType::ParentType; | |
| 405 | |||
| 406 | 5441776 | void bindLocalViews() | |
| 407 | { | ||
| 408 | // get some references for convenience | ||
| 409 | 5441776 | auto& couplingManager = this->couplingManager(); | |
| 410 | 5441776 | const auto& element = this->element(); | |
| 411 | 5441776 | const auto& curSol = this->curSol(); | |
| 412 | 5441776 | auto&& fvGeometry = this->fvGeometry(); | |
| 413 | 5441776 | auto&& curElemVolVars = this->curElemVolVars(); | |
| 414 | 10883552 | auto&& curElemFaceVars = this->curElemFaceVars(); | |
| 415 | 5441776 | auto&& elemFluxVarsCache = this->elemFluxVarsCache(); | |
| 416 | |||
| 417 | // bind the caches | ||
| 418 | 5441776 | couplingManager.bindCouplingContext(domainId, element, this->assembler()); | |
| 419 | 5441776 | fvGeometry.bind(element); | |
| 420 | 5441776 | curElemVolVars.bind(element, fvGeometry, curSol); | |
| 421 | 5441776 | curElemFaceVars.bind(element, fvGeometry, curSol); | |
| 422 | 5441776 | elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars); | |
| 423 | 5441776 | if (!this->assembler().isStationaryProblem()) | |
| 424 | { | ||
| 425 | this->prevElemVolVars().bindElement(element, fvGeometry, this->assembler().prevSol()); | ||
| 426 | this->prevElemFaceVars().bindElement(element, fvGeometry, this->assembler().prevSol()); | ||
| 427 | } | ||
| 428 | 5441776 | } | |
| 429 | }; | ||
| 430 | |||
| 431 | /*! | ||
| 432 | * \ingroup Assembly | ||
| 433 | * \ingroup StaggeredDiscretization | ||
| 434 | * \ingroup MultiDomain | ||
| 435 | * \brief The staggered multidomain local assembler | ||
| 436 | * \tparam id the id of the sub domain | ||
| 437 | * \tparam TypeTag the TypeTag | ||
| 438 | * \tparam Assembler the assembler type | ||
| 439 | * \tparam DM the numeric differentiation method | ||
| 440 | * \tparam implicit whether the assembler is explicit or implicit in time | ||
| 441 | */ | ||
| 442 | template<std::size_t id, class TypeTag, class Assembler, DiffMethod DM = DiffMethod::numeric, bool implicit = true> | ||
| 443 | class SubDomainStaggeredLocalAssembler; | ||
| 444 | |||
| 445 | /*! | ||
| 446 | * \ingroup Assembly | ||
| 447 | * \ingroup StaggeredDiscretization | ||
| 448 | * \ingroup MultiDomain | ||
| 449 | * \brief Staggered scheme local assembler using numeric differentiation and implicit time discretization | ||
| 450 | * | ||
| 451 | * The assembly of the cellCenterResidual is done element-wise, the assembly of the face residual is done half-element-wise. | ||
| 452 | * | ||
| 453 | * A sketch of what this means can be found in the following image: | ||
| 454 | * | ||
| 455 | * \image html staggered_halfelementwise.png | ||
| 456 | * | ||
| 457 | * Half-element wise assembly means, that integrals are | ||
| 458 | split into contributions from the left and right part of the staggered control volume. For an example term \f$\int_{\Omega}\varrho u\text{d}\Omega\f$ this reads | ||
| 459 | \f$\int_{\Omega}\varrho u\text{d}\Omega = \frac{1}{2}\Omega_\text{left}\varrho_\text{left}u+\frac{1}{2}\Omega_\text{right}\varrho_\text{right}u\f$. | ||
| 460 | |||
| 461 | During assembly, \f$\frac{1}{2}\Omega_\text{left}\varrho_\text{left}u\f$ is added to the residual of the | ||
| 462 | staggered control volume \f$\Omega\f$, when the loops reach the scvf within element \f$\Omega_\text{left}\f$. \f$\frac{1}{2}\Omega_\text{right}\varrho_\text{right}u\f$ is added to the residual of \f$\Omega\f$, when the loops reach the scvf within in element \f$\Omega_\text{right}\f$. | ||
| 463 | Other terms are split analogously. | ||
| 464 | */ | ||
| 465 | template<std::size_t id, class TypeTag, class Assembler> | ||
| 466 | 5441776 | class SubDomainStaggeredLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/true> | |
| 467 | : public SubDomainStaggeredLocalAssemblerImplicitBase<id, TypeTag, Assembler, | ||
| 468 | SubDomainStaggeredLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, true> > | ||
| 469 | { | ||
| 470 | using ThisType = SubDomainStaggeredLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/true>; | ||
| 471 | using ParentType = SubDomainStaggeredLocalAssemblerImplicitBase<id, TypeTag, Assembler, ThisType>; | ||
| 472 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 473 | using CellCenterResidualValue = typename ParentType::LocalResidual::CellCenterResidualValue; | ||
| 474 | using FaceResidualValue = typename ParentType::LocalResidual::FaceResidualValue; | ||
| 475 | using Element = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView::template Codim<0>::Entity; | ||
| 476 | using GridFaceVariables = GetPropType<TypeTag, Properties::GridFaceVariables>; | ||
| 477 | using ElementFaceVariables = typename GetPropType<TypeTag, Properties::GridFaceVariables>::LocalView; | ||
| 478 | using FaceVariables = typename ElementFaceVariables::FaceVariables; | ||
| 479 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 480 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 481 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 482 | using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>; | ||
| 483 | using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>; | ||
| 484 | using FacePrimaryVariables = GetPropType<TypeTag, Properties::FacePrimaryVariables>; | ||
| 485 | using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>; | ||
| 486 | |||
| 487 | static constexpr bool enableGridFluxVarsCache = GetPropType<TypeTag, Properties::GridVariables>::GridFluxVariablesCache::cachingEnabled; | ||
| 488 | static constexpr int maxNeighbors = 4*(2*ModelTraits::dim()); | ||
| 489 | static constexpr auto domainI = Dune::index_constant<id>(); | ||
| 490 | static constexpr auto cellCenterId = GridGeometry::cellCenterIdx(); | ||
| 491 | static constexpr auto faceId = GridGeometry::faceIdx(); | ||
| 492 | |||
| 493 | static constexpr auto numEq = ModelTraits::numEq(); | ||
| 494 | static constexpr auto numEqCellCenter = CellCenterPrimaryVariables::dimension; | ||
| 495 | static constexpr auto numEqFace = FacePrimaryVariables::dimension; | ||
| 496 | |||
| 497 | public: | ||
| 498 |
2/8✓ Branch 2 taken 1360444 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 1360444 times.
✗ Branch 7 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
|
2720888 | using ParentType::ParentType; |
| 499 | |||
| 500 | CellCenterResidualValue assembleCellCenterResidualImpl() | ||
| 501 | { | ||
| 502 | ✗ | return this->evalLocalResidualForCellCenter(); | |
| 503 | } | ||
| 504 | |||
| 505 | FaceResidualValue assembleFaceResidualImpl(const SubControlVolumeFace& scvf) | ||
| 506 | { | ||
| 507 | ✗ | return this->evalLocalResidualForFace(scvf); | |
| 508 | } | ||
| 509 | |||
| 510 | /*! | ||
| 511 | * \brief Computes the derivatives with respect to the given element and adds them | ||
| 512 | * to the global matrix. | ||
| 513 | * | ||
| 514 | * \return The element residual at the current solution. | ||
| 515 | */ | ||
| 516 | template<class JacobianMatrixDiagBlock, class GridVariables> | ||
| 517 | 1360444 | CellCenterResidualValue assembleCellCenterJacobianAndResidualImpl(JacobianMatrixDiagBlock& A, GridVariables& gridVariables) | |
| 518 | { | ||
| 519 | assert(domainI == cellCenterId); | ||
| 520 | |||
| 521 | // get some aliases for convenience | ||
| 522 | 1360444 | const auto& element = this->element(); | |
| 523 | 1360444 | const auto& fvGeometry = this->fvGeometry(); | |
| 524 | 1360444 | auto&& curElemVolVars = this->curElemVolVars(); | |
| 525 | 2720888 | const auto& gridGeometry = this->problem().gridGeometry(); | |
| 526 | 1360444 | const auto& curSol = this->curSol()[domainI]; | |
| 527 | |||
| 528 | 2720888 | const auto cellCenterGlobalI = gridGeometry.elementMapper().index(element); | |
| 529 |
1/2✓ Branch 0 taken 354988 times.
✗ Branch 1 not taken.
|
1360444 | const auto origResidual = this->evalLocalResidualForCellCenter(); |
| 530 | |||
| 531 | ///////////////////////////////////////////////////////////////////////////////////////////////////////// | ||
| 532 | // Calculate derivatives of all cell center residuals in the element w.r.t. to other cell center dofs. // | ||
| 533 | ///////////////////////////////////////////////////////////////////////////////////////////////////////// | ||
| 534 | |||
| 535 | // lambda to evaluate the derivatives for cell center dofs with respect to neighbor cells | ||
| 536 | 6501890 | auto evaluateCellCenterDerivatives = [&](const std::size_t globalJ) | |
| 537 | { | ||
| 538 | // get the volVars of the element with respect to which we are going to build the derivative | ||
| 539 | 6501890 | auto&& scvJ = fvGeometry.scv(globalJ); | |
| 540 | 6501890 | const auto elementJ = fvGeometry.gridGeometry().element(globalJ); | |
| 541 | 13003780 | auto& curVolVars = this->getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scvJ); | |
| 542 | 6501890 | const auto origVolVars(curVolVars); | |
| 543 | |||
| 544 |
2/2✓ Branch 0 taken 17233802 times.
✓ Branch 1 taken 6501890 times.
|
23735692 | for (int pvIdx = 0; pvIdx < numEqCellCenter; ++pvIdx) |
| 545 | { | ||
| 546 | 17233802 | CellCenterPrimaryVariables cellCenterPriVars = curSol[globalJ]; | |
| 547 | using PrimaryVariables = typename VolumeVariables::PrimaryVariables; | ||
| 548 | 17233802 | PrimaryVariables priVars = makePriVarsFromCellCenterPriVars<PrimaryVariables>(cellCenterPriVars); | |
| 549 | |||
| 550 | 17233802 | constexpr auto offset = numEq - numEqCellCenter; | |
| 551 | |||
| 552 | 222104284 | auto evalResidual = [&](Scalar priVar) | |
| 553 | { | ||
| 554 | // update the volume variables | ||
| 555 | 31272780 | priVars[pvIdx + offset] = priVar; | |
| 556 |
1/4✓ Branch 1 taken 29404324 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
60677104 | auto elemSol = elementSolution<FVElementGeometry>(std::move(priVars)); |
| 557 |
2/4✓ Branch 1 taken 31272780 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 31272780 times.
✗ Branch 5 not taken.
|
62545560 | curVolVars.update(elemSol, this->problem(), elementJ, scvJ); |
| 558 | |||
| 559 | // update the coupling context | ||
| 560 |
2/2✓ Branch 0 taken 1564704 times.
✓ Branch 1 taken 27839620 times.
|
33141236 | cellCenterPriVars[pvIdx] = priVar; |
| 561 |
2/2✓ Branch 0 taken 2892224 times.
✓ Branch 1 taken 28380556 times.
|
31272780 | this->couplingManager().updateCouplingContext(domainI, *this, domainI, globalJ, cellCenterPriVars, pvIdx); |
| 562 | |||
| 563 | // compute element residual | ||
| 564 |
3/4✓ Branch 0 taken 1868456 times.
✓ Branch 1 taken 29404324 times.
✓ Branch 2 taken 1868456 times.
✗ Branch 3 not taken.
|
64414016 | return this->evalLocalResidualForCellCenter(); |
| 565 | }; | ||
| 566 | |||
| 567 | // create the vector storing the partial derivatives | ||
| 568 |
2/2✓ Branch 0 taken 14 times.
✓ Branch 1 taken 1749274 times.
|
17233802 | CellCenterResidualValue partialDeriv(0.0); |
| 569 | |||
| 570 | // derive the residuals numerically | ||
| 571 |
4/4✓ Branch 0 taken 51 times.
✓ Branch 1 taken 17233751 times.
✓ Branch 2 taken 51 times.
✓ Branch 3 taken 17233751 times.
|
34467604 | const auto& paramGroup = this->problem().paramGroup(); |
| 572 |
4/6✓ Branch 0 taken 51 times.
✓ Branch 1 taken 17233751 times.
✓ Branch 3 taken 51 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 51 times.
✗ Branch 7 not taken.
|
17233802 | static const int numDiffMethod = getParamFromGroup<int>(paramGroup, "Assembly.NumericDifferenceMethod"); |
| 573 |
4/6✓ Branch 0 taken 51 times.
✓ Branch 1 taken 17233751 times.
✓ Branch 3 taken 51 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 51 times.
✗ Branch 7 not taken.
|
17233802 | static const auto eps = this->couplingManager().numericEpsilon(domainI, paramGroup); |
| 574 |
2/4✓ Branch 1 taken 17233802 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 17233802 times.
✗ Branch 5 not taken.
|
34467604 | NumericDifferentiation::partialDerivative(evalResidual, priVars[pvIdx + offset], partialDeriv, origResidual, |
| 575 | 34467604 | eps(priVars[pvIdx + offset], pvIdx), numDiffMethod); | |
| 576 | |||
| 577 | // update the global jacobian matrix with the current partial derivatives | ||
| 578 |
1/2✓ Branch 1 taken 17233802 times.
✗ Branch 2 not taken.
|
17233802 | updateGlobalJacobian_(A, cellCenterGlobalI, globalJ, pvIdx, partialDeriv); |
| 579 | |||
| 580 | // restore the original volVars | ||
| 581 | 17233802 | curVolVars = origVolVars; | |
| 582 | |||
| 583 | // restore the undeflected state of the coupling context | ||
| 584 |
4/4✓ Branch 0 taken 2892224 times.
✓ Branch 1 taken 14341578 times.
✓ Branch 2 taken 2892224 times.
✓ Branch 3 taken 14341578 times.
|
37359828 | this->couplingManager().updateCouplingContext(domainI, *this, domainI, globalJ, curSol[globalJ], pvIdx); |
| 585 | } | ||
| 586 | }; | ||
| 587 | |||
| 588 | // get the list of cell center dofs that have an influence on the cell center resdiual of the current element | ||
| 589 | 1360444 | const auto& connectivityMap = gridGeometry.connectivityMap(); | |
| 590 | |||
| 591 | // evaluate derivatives w.r.t. own dof | ||
| 592 | 1360444 | evaluateCellCenterDerivatives(cellCenterGlobalI); | |
| 593 | |||
| 594 | // evaluate derivatives w.r.t. all other related cell center dofs | ||
| 595 |
4/4✓ Branch 0 taken 5141446 times.
✓ Branch 1 taken 1360444 times.
✓ Branch 2 taken 5141446 times.
✓ Branch 3 taken 1360444 times.
|
10583222 | for (const auto& globalJ : connectivityMap(cellCenterId, cellCenterId, cellCenterGlobalI)) |
| 596 | 5141446 | evaluateCellCenterDerivatives(globalJ); | |
| 597 | |||
| 598 | 1360444 | return origResidual; | |
| 599 | } | ||
| 600 | |||
| 601 | /*! | ||
| 602 | * \brief Computes the derivatives with respect to the given element and adds them | ||
| 603 | * to the global matrix. | ||
| 604 | * | ||
| 605 | * \return The element residual at the current solution. | ||
| 606 | */ | ||
| 607 | template<class JacobianMatrixDiagBlock, class GridVariables> | ||
| 608 | 1360444 | auto assembleFaceJacobianAndResidualImpl(JacobianMatrixDiagBlock& A, GridVariables& gridVariables) | |
| 609 | { | ||
| 610 | assert(domainI == faceId); | ||
| 611 | |||
| 612 | // get some aliases for convenience | ||
| 613 |
1/2✓ Branch 1 taken 1360444 times.
✗ Branch 2 not taken.
|
1360444 | const auto& problem = this->problem(); |
| 614 | 1360444 | const auto& element = this->element(); | |
| 615 |
1/2✓ Branch 1 taken 1360444 times.
✗ Branch 2 not taken.
|
1360444 | const auto& fvGeometry = this->fvGeometry(); |
| 616 |
2/4✓ Branch 1 taken 1360444 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1360444 times.
✗ Branch 5 not taken.
|
2720888 | const auto& gridGeometry = this->problem().gridGeometry(); |
| 617 |
1/2✓ Branch 1 taken 1360444 times.
✗ Branch 2 not taken.
|
1360444 | const auto& curSol = this->curSol()[domainI]; |
| 618 | |||
| 619 | using FaceSolutionVector = GetPropType<TypeTag, Properties::FaceSolutionVector>; // TODO: use reserved vector | ||
| 620 |
1/2✓ Branch 1 taken 1360444 times.
✗ Branch 2 not taken.
|
1360444 | FaceSolutionVector origResiduals; |
| 621 |
2/4✓ Branch 1 taken 1360444 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1360444 times.
✗ Branch 5 not taken.
|
2720888 | origResiduals.resize(fvGeometry.numScvf()); |
| 622 | 1360444 | origResiduals = 0.0; | |
| 623 | |||
| 624 | // treat the local residua of the face dofs: | ||
| 625 |
5/6✓ Branch 0 taken 5441776 times.
✓ Branch 1 taken 1360444 times.
✓ Branch 2 taken 5441776 times.
✓ Branch 3 taken 1360444 times.
✓ Branch 4 taken 5441776 times.
✗ Branch 5 not taken.
|
13604440 | for (auto&& scvf : scvfs(fvGeometry)) |
| 626 |
1/2✓ Branch 0 taken 5441776 times.
✗ Branch 1 not taken.
|
21767104 | origResiduals[scvf.localFaceIdx()] = this->evalLocalResidualForFace(scvf); |
| 627 | |||
| 628 | /////////////////////////////////////////////////////////////////////////////////////////////////// | ||
| 629 | // Calculate derivatives of all face residuals in the element w.r.t. to other face dofs. // | ||
| 630 | // Note that we do an element-wise assembly, therefore this is only the contribution of the // | ||
| 631 | // current element while the contribution of the element on the opposite side of the scvf will // | ||
| 632 | // be added separately. // | ||
| 633 | /////////////////////////////////////////////////////////////////////////////////////////////////// | ||
| 634 | |||
| 635 |
5/6✓ Branch 0 taken 5441776 times.
✓ Branch 1 taken 1360444 times.
✓ Branch 2 taken 5441776 times.
✓ Branch 3 taken 1360444 times.
✓ Branch 5 taken 5441776 times.
✗ Branch 6 not taken.
|
8162664 | for (auto&& scvf : scvfs(fvGeometry)) |
| 636 | { | ||
| 637 | // set the actual dof index | ||
| 638 |
1/2✓ Branch 1 taken 5441776 times.
✗ Branch 2 not taken.
|
5441776 | const auto faceGlobalI = scvf.dofIndex(); |
| 639 | |||
| 640 | using FaceSolution = GetPropType<TypeTag, Properties::StaggeredFaceSolution>; | ||
| 641 |
1/2✓ Branch 1 taken 5441776 times.
✗ Branch 2 not taken.
|
10883552 | const auto origFaceSolution = FaceSolution(scvf, curSol, gridGeometry); |
| 642 | |||
| 643 | // Lambda to evaluate the derivatives for faces | ||
| 644 | 42820088 | auto evaluateFaceDerivatives = [&](const std::size_t globalJ) | |
| 645 | { | ||
| 646 | // get the faceVars of the face with respect to which we are going to build the derivative | ||
| 647 | 42820088 | auto& faceVars = getFaceVarAccess_(gridVariables.curGridFaceVars(), this->curElemFaceVars(), scvf); | |
| 648 | 42820088 | const auto origFaceVars = faceVars; | |
| 649 | |||
| 650 |
2/2✓ Branch 0 taken 42820088 times.
✓ Branch 1 taken 42820088 times.
|
85640176 | for (int pvIdx = 0; pvIdx < numEqFace; ++pvIdx) |
| 651 | { | ||
| 652 |
1/2✓ Branch 1 taken 42820088 times.
✗ Branch 2 not taken.
|
42820088 | auto faceSolution = origFaceSolution; |
| 653 | |||
| 654 | 818457808 | auto evalResidual = [&](Scalar priVar) | |
| 655 | { | ||
| 656 | // update the volume variables | ||
| 657 | 70192160 | faceSolution[globalJ][pvIdx] = priVar; | |
| 658 | 70192160 | faceVars.update(faceSolution, problem, element, fvGeometry, scvf); | |
| 659 | |||
| 660 | // update the coupling context | ||
| 661 |
2/2✓ Branch 1 taken 13012416 times.
✓ Branch 2 taken 57179744 times.
|
70192160 | this->couplingManager().updateCouplingContext(domainI, *this, domainI, globalJ, faceSolution[globalJ], pvIdx); |
| 662 | |||
| 663 | // compute face residual | ||
| 664 |
1/2✓ Branch 0 taken 70192160 times.
✗ Branch 1 not taken.
|
70192160 | return this->evalLocalResidualForFace(scvf); |
| 665 | }; | ||
| 666 | |||
| 667 | // derive the residuals numerically | ||
| 668 |
2/2✓ Branch 0 taken 51 times.
✓ Branch 1 taken 42820037 times.
|
42820088 | FaceResidualValue partialDeriv(0.0); |
| 669 |
2/2✓ Branch 0 taken 51 times.
✓ Branch 1 taken 42820037 times.
|
42820088 | const auto& paramGroup = problem.paramGroup(); |
| 670 |
4/6✓ Branch 0 taken 51 times.
✓ Branch 1 taken 42820037 times.
✓ Branch 3 taken 51 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 51 times.
✗ Branch 7 not taken.
|
42820088 | static const int numDiffMethod = getParamFromGroup<int>(paramGroup, "Assembly.NumericDifferenceMethod"); |
| 671 |
4/6✓ Branch 0 taken 51 times.
✓ Branch 1 taken 42820037 times.
✓ Branch 3 taken 51 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 51 times.
✗ Branch 7 not taken.
|
42820088 | static const auto eps = this->couplingManager().numericEpsilon(domainI, paramGroup); |
| 672 |
1/2✓ Branch 4 taken 42820088 times.
✗ Branch 5 not taken.
|
128460264 | NumericDifferentiation::partialDerivative(evalResidual, faceSolution[globalJ][pvIdx], partialDeriv, origResiduals[scvf.localFaceIdx()], |
| 673 | 42820088 | eps(faceSolution[globalJ][pvIdx], pvIdx), numDiffMethod); | |
| 674 | |||
| 675 | // update the global jacobian matrix with the current partial derivatives | ||
| 676 |
1/2✓ Branch 1 taken 42820088 times.
✗ Branch 2 not taken.
|
42820088 | updateGlobalJacobian_(A, faceGlobalI, globalJ, pvIdx, partialDeriv); |
| 677 | |||
| 678 | // restore the original faceVars | ||
| 679 | 42820088 | faceVars = origFaceVars; | |
| 680 | |||
| 681 | // restore the undeflected state of the coupling context | ||
| 682 |
2/2✓ Branch 1 taken 13012416 times.
✓ Branch 2 taken 29807672 times.
|
55832504 | this->couplingManager().updateCouplingContext(domainI, *this, domainI, globalJ, origFaceSolution[globalJ], pvIdx); |
| 683 | } | ||
| 684 | }; | ||
| 685 | |||
| 686 | // evaluate derivatives w.r.t. own dof | ||
| 687 |
2/4✓ Branch 1 taken 5441776 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5441776 times.
✗ Branch 5 not taken.
|
10883552 | evaluateFaceDerivatives(scvf.dofIndex()); |
| 688 | |||
| 689 | // get the list of face dofs that have an influence on the resdiual of the current face | ||
| 690 | 5441776 | const auto& connectivityMap = gridGeometry.connectivityMap(); | |
| 691 | |||
| 692 | // evaluate derivatives w.r.t. all other related face dofs | ||
| 693 |
4/4✓ Branch 0 taken 37378312 times.
✓ Branch 1 taken 5441776 times.
✓ Branch 2 taken 37378312 times.
✓ Branch 3 taken 5441776 times.
|
59145416 | for (const auto& globalJ : connectivityMap(faceId, faceId, scvf.index())) |
| 694 |
1/2✓ Branch 1 taken 37378312 times.
✗ Branch 2 not taken.
|
37378312 | evaluateFaceDerivatives(globalJ); |
| 695 | } | ||
| 696 | |||
| 697 | 1360444 | return origResiduals; | |
| 698 | } | ||
| 699 | |||
| 700 | /*! | ||
| 701 | * \brief Computes the derivatives with respect to the given element and adds them | ||
| 702 | * to the global matrix. | ||
| 703 | */ | ||
| 704 | template<class JacobianBlock, class GridVariables> | ||
| 705 | 1360444 | void assembleJacobianCellCenterCoupling(Dune::index_constant<faceId> domainJ, JacobianBlock& A, | |
| 706 | const CellCenterResidualValue& origResidual, GridVariables& gridVariables) | ||
| 707 | { | ||
| 708 | ///////////////////////////////////////////////////////////////////////////////////////////////////////// | ||
| 709 | // Calculate derivatives of all cell center residuals in the element w.r.t. to all coupled faces dofs // | ||
| 710 | ///////////////////////////////////////////////////////////////////////////////////////////////////////// | ||
| 711 | |||
| 712 | // get some aliases for convenience | ||
| 713 | 1360444 | const auto& element = this->element(); | |
| 714 | 1360444 | const auto& fvGeometry = this->fvGeometry(); | |
| 715 | 2720888 | const auto& gridGeometry = this->problem().gridGeometry(); | |
| 716 | 1360444 | const auto& curSol = this->curSol()[domainJ]; | |
| 717 | // build derivatives with for cell center dofs w.r.t. cell center dofs | ||
| 718 | 2720888 | const auto cellCenterGlobalI = gridGeometry.elementMapper().index(element); | |
| 719 | |||
| 720 |
4/4✓ Branch 0 taken 5441776 times.
✓ Branch 1 taken 1360444 times.
✓ Branch 2 taken 5441776 times.
✓ Branch 3 taken 1360444 times.
|
8162664 | for (const auto& scvfJ : scvfs(fvGeometry)) |
| 721 | { | ||
| 722 | 5441776 | const auto globalJ = scvfJ.dofIndex(); | |
| 723 | |||
| 724 | // get the faceVars of the face with respect to which we are going to build the derivative | ||
| 725 | 5441776 | auto& faceVars = getFaceVarAccess_(gridVariables.curGridFaceVars(), this->curElemFaceVars(), scvfJ); | |
| 726 | 5441776 | const auto origFaceVars(faceVars); | |
| 727 | |||
| 728 |
2/2✓ Branch 0 taken 5441776 times.
✓ Branch 1 taken 5441776 times.
|
10883552 | for (int pvIdx = 0; pvIdx < numEqFace; ++pvIdx) |
| 729 | { | ||
| 730 |
2/2✓ Branch 0 taken 14 times.
✓ Branch 1 taken 1419938 times.
|
5441776 | auto facePriVars = curSol[globalJ]; |
| 731 | |||
| 732 | 32399488 | auto evalResidual = [&](Scalar priVar) | |
| 733 | { | ||
| 734 | // update the face variables | ||
| 735 | 17971808 | facePriVars[pvIdx] = priVar; | |
| 736 |
2/2✓ Branch 0 taken 1651648 times.
✓ Branch 1 taken 7334256 times.
|
8985904 | faceVars.updateOwnFaceOnly(facePriVars); |
| 737 | |||
| 738 | // update the coupling context | ||
| 739 |
2/2✓ Branch 0 taken 1651648 times.
✓ Branch 1 taken 7334256 times.
|
8985904 | this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, facePriVars, pvIdx); |
| 740 | |||
| 741 | // compute element residual | ||
| 742 |
1/2✓ Branch 0 taken 1520304 times.
✗ Branch 1 not taken.
|
8985904 | return this->evalLocalResidualForCellCenter(); |
| 743 | }; | ||
| 744 | |||
| 745 | // create the vector storing the partial derivatives | ||
| 746 |
2/2✓ Branch 0 taken 14 times.
✓ Branch 1 taken 1419938 times.
|
5441776 | CellCenterResidualValue partialDeriv(0.0); |
| 747 | |||
| 748 | // derive the residuals numerically | ||
| 749 |
4/4✓ Branch 0 taken 51 times.
✓ Branch 1 taken 5441725 times.
✓ Branch 2 taken 51 times.
✓ Branch 3 taken 5441725 times.
|
10883552 | const auto& paramGroup = this->assembler().problem(domainJ).paramGroup(); |
| 750 |
4/6✓ Branch 0 taken 51 times.
✓ Branch 1 taken 5441725 times.
✓ Branch 3 taken 51 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 51 times.
✗ Branch 7 not taken.
|
5441776 | static const int numDiffMethod = getParamFromGroup<int>(paramGroup, "Assembly.NumericDifferenceMethod"); |
| 751 |
4/6✓ Branch 0 taken 51 times.
✓ Branch 1 taken 5441725 times.
✓ Branch 3 taken 51 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 51 times.
✗ Branch 7 not taken.
|
5441776 | static const auto epsCoupl = this->couplingManager().numericEpsilon(domainJ, paramGroup); |
| 752 |
1/2✓ Branch 1 taken 5441776 times.
✗ Branch 2 not taken.
|
5441776 | NumericDifferentiation::partialDerivative(evalResidual, facePriVars[pvIdx], partialDeriv, origResidual, |
| 753 | 5441776 | epsCoupl(facePriVars[pvIdx], pvIdx), numDiffMethod); | |
| 754 | |||
| 755 | // update the global jacobian matrix with the current partial derivatives | ||
| 756 |
1/2✓ Branch 1 taken 5441776 times.
✗ Branch 2 not taken.
|
5441776 | updateGlobalJacobian_(A, cellCenterGlobalI, globalJ, pvIdx, partialDeriv); |
| 757 | |||
| 758 | // restore the original faceVars | ||
| 759 | 5441776 | faceVars = origFaceVars; | |
| 760 | |||
| 761 | // restore the undeflected state of the coupling context | ||
| 762 |
4/4✓ Branch 0 taken 1651648 times.
✓ Branch 1 taken 3790128 times.
✓ Branch 2 taken 1651648 times.
✓ Branch 3 taken 3790128 times.
|
12535200 | this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, curSol[globalJ], pvIdx); |
| 763 | } | ||
| 764 | } | ||
| 765 | 1360444 | } | |
| 766 | |||
| 767 | template<std::size_t otherId, class JacobianBlock, class GridVariables> | ||
| 768 | 412912 | void assembleJacobianCellCenterCoupling(Dune::index_constant<otherId> domainJ, JacobianBlock& A, | |
| 769 | const CellCenterResidualValue& res, GridVariables& gridVariables) | ||
| 770 | { | ||
| 771 | ///////////////////////////////////////////////////////////////////////////////////////////////////////////////// | ||
| 772 | // Calculate derivatives of all cell center residuals in the element w.r.t. all dofs in the coupling stencil. // | ||
| 773 | ///////////////////////////////////////////////////////////////////////////////////////////////////////////////// | ||
| 774 | |||
| 775 | // get some aliases for convenience | ||
| 776 | 412912 | const auto& element = this->element(); | |
| 777 | |||
| 778 | // get stencil information | ||
| 779 | 412912 | const auto& stencil = this->couplingManager().couplingStencil(domainI, element, domainJ); | |
| 780 | |||
| 781 |
4/4✓ Branch 0 taken 11332 times.
✓ Branch 1 taken 401580 times.
✓ Branch 2 taken 11332 times.
✓ Branch 3 taken 401580 times.
|
825824 | if (stencil.empty()) |
| 782 | return; | ||
| 783 | |||
| 784 |
4/4✓ Branch 0 taken 11332 times.
✓ Branch 1 taken 11332 times.
✓ Branch 2 taken 11332 times.
✓ Branch 3 taken 11332 times.
|
33996 | for (const auto globalJ : stencil) |
| 785 | { | ||
| 786 |
2/3✓ Branch 0 taken 3620 times.
✓ Branch 1 taken 7712 times.
✗ Branch 2 not taken.
|
11332 | const auto origResidual = this->couplingManager().evalCouplingResidual(domainI, *this, domainJ, globalJ); |
| 787 | 11332 | const auto& curSol = this->curSol()[domainJ]; | |
| 788 | 11332 | const auto origPriVarsJ = curSol[globalJ]; | |
| 789 | |||
| 790 |
2/2✓ Branch 0 taken 22896 times.
✓ Branch 1 taken 11332 times.
|
34228 | for (int pvIdx = 0; pvIdx < JacobianBlock::block_type::cols; ++pvIdx) |
| 791 | { | ||
| 792 | 112440 | auto evalCouplingResidual = [&](Scalar priVar) | |
| 793 | { | ||
| 794 | 20856 | auto deflectedPriVarsJ = origPriVarsJ; | |
| 795 | 22896 | deflectedPriVarsJ[pvIdx] = priVar; | |
| 796 | 22896 | this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, deflectedPriVarsJ, pvIdx); | |
| 797 |
1/2✓ Branch 0 taken 5200 times.
✗ Branch 1 not taken.
|
22896 | return this->couplingManager().evalCouplingResidual(domainI, *this, domainJ, globalJ); |
| 798 | }; | ||
| 799 | |||
| 800 | // create the vector storing the partial derivatives | ||
| 801 |
2/2✓ Branch 0 taken 9 times.
✓ Branch 1 taken 5191 times.
|
22896 | CellCenterResidualValue partialDeriv(0.0); |
| 802 | |||
| 803 | // derive the residuals numerically | ||
| 804 |
4/4✓ Branch 0 taken 17 times.
✓ Branch 1 taken 22879 times.
✓ Branch 2 taken 17 times.
✓ Branch 3 taken 22879 times.
|
45792 | const auto& paramGroup = this->assembler().problem(domainJ).paramGroup(); |
| 805 |
4/6✓ Branch 0 taken 17 times.
✓ Branch 1 taken 22879 times.
✓ Branch 3 taken 17 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 17 times.
✗ Branch 7 not taken.
|
22896 | static const int numDiffMethod = getParamFromGroup<int>(paramGroup, "Assembly.NumericDifferenceMethod"); |
| 806 |
4/6✓ Branch 0 taken 17 times.
✓ Branch 1 taken 22879 times.
✓ Branch 3 taken 17 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 17 times.
✗ Branch 7 not taken.
|
22896 | static const auto epsCoupl = this->couplingManager().numericEpsilon(domainJ, paramGroup); |
| 807 |
2/4✓ Branch 1 taken 22896 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20856 times.
✗ Branch 5 not taken.
|
43752 | NumericDifferentiation::partialDerivative(evalCouplingResidual, origPriVarsJ[pvIdx], partialDeriv, origResidual, |
| 808 | 43752 | epsCoupl(origPriVarsJ[pvIdx], pvIdx), numDiffMethod); | |
| 809 | |||
| 810 | // update the global stiffness matrix with the current partial derivatives | ||
| 811 | 91584 | const auto cellCenterGlobalI = this->problem().gridGeometry().elementMapper().index(element); | |
| 812 |
1/2✓ Branch 1 taken 22896 times.
✗ Branch 2 not taken.
|
22896 | updateGlobalJacobian_(A, cellCenterGlobalI, globalJ, pvIdx, partialDeriv); |
| 813 | |||
| 814 | // restore the undeflected state of the coupling context | ||
| 815 |
1/2✓ Branch 1 taken 22896 times.
✗ Branch 2 not taken.
|
22896 | this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, origPriVarsJ, pvIdx); |
| 816 | } | ||
| 817 | } | ||
| 818 | } | ||
| 819 | |||
| 820 | /*! | ||
| 821 | * \brief Computes the derivatives with respect to the given element and adds them | ||
| 822 | * to the global matrix. | ||
| 823 | */ | ||
| 824 | template<class JacobianBlock, class ElementResidualVector, class GridVariables> | ||
| 825 | 1360444 | void assembleJacobianFaceCoupling(Dune::index_constant<cellCenterId> domainJ, JacobianBlock& A, | |
| 826 | const ElementResidualVector& origResiduals, GridVariables& gridVariables) | ||
| 827 | { | ||
| 828 | ///////////////////////////////////////////////////////////////////////////////////////////////////// | ||
| 829 | // Calculate derivatives of all face residuals in the element w.r.t. all coupled cell center dofs // | ||
| 830 | ///////////////////////////////////////////////////////////////////////////////////////////////////// | ||
| 831 | |||
| 832 | // get some aliases for convenience | ||
| 833 | 1360444 | const auto& problem = this->problem(); | |
| 834 | 1360444 | const auto& fvGeometry = this->fvGeometry(); | |
| 835 | 2720888 | const auto& gridGeometry = this->problem().gridGeometry(); | |
| 836 | 1360444 | const auto& connectivityMap = gridGeometry.connectivityMap(); | |
| 837 | 1360444 | const auto& curSol = this->curSol()[domainJ]; | |
| 838 | |||
| 839 | // build derivatives with for cell center dofs w.r.t. cell center dofs | ||
| 840 |
4/4✓ Branch 0 taken 5441776 times.
✓ Branch 1 taken 1360444 times.
✓ Branch 2 taken 5441776 times.
✓ Branch 3 taken 1360444 times.
|
8162664 | for (auto&& scvf : scvfs(fvGeometry)) |
| 841 | { | ||
| 842 | // set the actual dof index | ||
| 843 | 5441776 | const auto faceGlobalI = scvf.dofIndex(); | |
| 844 | |||
| 845 | // build derivatives with for face dofs w.r.t. cell center dofs | ||
| 846 |
4/4✓ Branch 0 taken 15724668 times.
✓ Branch 1 taken 5441776 times.
✓ Branch 2 taken 15724668 times.
✓ Branch 3 taken 5441776 times.
|
37491772 | for (const auto& globalJ : connectivityMap(faceId, cellCenterId, scvf.index())) |
| 847 | { | ||
| 848 | // get the volVars of the element with respect to which we are going to build the derivative | ||
| 849 |
1/2✓ Branch 1 taken 15724668 times.
✗ Branch 2 not taken.
|
15724668 | auto&& scvJ = fvGeometry.scv(globalJ); |
| 850 |
1/2✓ Branch 1 taken 15724668 times.
✗ Branch 2 not taken.
|
15724668 | const auto elementJ = fvGeometry.gridGeometry().element(globalJ); |
| 851 | 47174004 | auto& curVolVars = this->getVolVarAccess(gridVariables.curGridVolVars(), this->curElemVolVars(), scvJ); | |
| 852 | 15724668 | const auto origVolVars(curVolVars); | |
| 853 | 15724668 | const auto origCellCenterPriVars = curSol[globalJ]; | |
| 854 | |||
| 855 |
2/2✓ Branch 0 taken 41738052 times.
✓ Branch 1 taken 15724668 times.
|
57462720 | for (int pvIdx = 0; pvIdx < numEqCellCenter; ++pvIdx) |
| 856 | { | ||
| 857 | using PrimaryVariables = typename VolumeVariables::PrimaryVariables; | ||
| 858 | 41738052 | PrimaryVariables priVars = makePriVarsFromCellCenterPriVars<PrimaryVariables>(origCellCenterPriVars); | |
| 859 | |||
| 860 | 41738052 | constexpr auto offset = PrimaryVariables::dimension - CellCenterPrimaryVariables::dimension; | |
| 861 | |||
| 862 | 344808388 | auto evalResidual = [&](Scalar priVar) | |
| 863 | { | ||
| 864 | // update the volume variables | ||
| 865 | 75767584 | priVars[pvIdx + offset] = priVar; | |
| 866 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
75767584 | auto elemSol = elementSolution<FVElementGeometry>(std::move(priVars)); |
| 867 |
1/2✓ Branch 1 taken 75767584 times.
✗ Branch 2 not taken.
|
75767584 | curVolVars.update(elemSol, problem, elementJ, scvJ); |
| 868 | |||
| 869 | // update the coupling context | ||
| 870 | 75767584 | auto deflectedCellCenterPriVars = origCellCenterPriVars; | |
| 871 |
2/2✓ Branch 0 taken 3788480 times.
✓ Branch 1 taken 67482040 times.
|
75767584 | deflectedCellCenterPriVars[pvIdx] = priVar; |
| 872 |
2/2✓ Branch 0 taken 6980320 times.
✓ Branch 1 taken 68787264 times.
|
75767584 | this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, deflectedCellCenterPriVars, pvIdx); |
| 873 | |||
| 874 | // compute face residual | ||
| 875 |
2/4✓ Branch 0 taken 75767584 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 75767584 times.
✗ Branch 3 not taken.
|
227302752 | return this->evalLocalResidualForFace(scvf); |
| 876 | }; | ||
| 877 | |||
| 878 | // derive the residuals numerically | ||
| 879 |
2/2✓ Branch 0 taken 51 times.
✓ Branch 1 taken 41738001 times.
|
41738052 | FaceResidualValue partialDeriv(0.0); |
| 880 |
4/4✓ Branch 0 taken 51 times.
✓ Branch 1 taken 41738001 times.
✓ Branch 2 taken 51 times.
✓ Branch 3 taken 41738001 times.
|
83476104 | const auto& paramGroup = this->assembler().problem(domainJ).paramGroup(); |
| 881 |
4/6✓ Branch 0 taken 51 times.
✓ Branch 1 taken 41738001 times.
✓ Branch 3 taken 51 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 51 times.
✗ Branch 7 not taken.
|
41738052 | static const int numDiffMethod = getParamFromGroup<int>(paramGroup, "Assembly.NumericDifferenceMethod"); |
| 882 |
4/6✓ Branch 0 taken 51 times.
✓ Branch 1 taken 41738001 times.
✓ Branch 3 taken 51 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 51 times.
✗ Branch 7 not taken.
|
41738052 | static const auto epsCoupl = this->couplingManager().numericEpsilon(domainJ, paramGroup); |
| 883 |
4/8✓ Branch 1 taken 41738052 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 41738052 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 41738052 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 41738052 times.
✗ Branch 11 not taken.
|
166952208 | NumericDifferentiation::partialDerivative(evalResidual, priVars[pvIdx + offset], partialDeriv, origResiduals[scvf.localFaceIdx()], |
| 884 | 83476104 | epsCoupl(priVars[pvIdx + offset], pvIdx), numDiffMethod); | |
| 885 | |||
| 886 | // update the global jacobian matrix with the current partial derivatives | ||
| 887 |
1/2✓ Branch 1 taken 41738052 times.
✗ Branch 2 not taken.
|
41738052 | updateGlobalJacobian_(A, faceGlobalI, globalJ, pvIdx, partialDeriv); |
| 888 | |||
| 889 | // restore the original volVars | ||
| 890 | 41738052 | curVolVars = origVolVars; | |
| 891 | |||
| 892 | // restore the undeflected state of the coupling context | ||
| 893 |
2/2✓ Branch 0 taken 6980320 times.
✓ Branch 1 taken 34757732 times.
|
48718372 | this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, origCellCenterPriVars, pvIdx); |
| 894 | } | ||
| 895 | } | ||
| 896 | } | ||
| 897 | 1360444 | } | |
| 898 | |||
| 899 | template<std::size_t otherId, class JacobianBlock, class ElementResidualVector, class GridVariables> | ||
| 900 | 412912 | void assembleJacobianFaceCoupling(Dune::index_constant<otherId> domainJ, JacobianBlock& A, | |
| 901 | const ElementResidualVector& res, GridVariables& gridVariables) | ||
| 902 | { | ||
| 903 | ////////////////////////////////////////////////////////////////////////////////////////////////////////// | ||
| 904 | // Calculate derivatives of all face residuals in the element w.r.t. all dofs in the coupling stencil. // | ||
| 905 | ////////////////////////////////////////////////////////////////////////////////////////////////////////// | ||
| 906 | |||
| 907 | // get some aliases for convenience | ||
| 908 | 412912 | const auto& fvGeometry = this->fvGeometry(); | |
| 909 | 412912 | const auto& curSol = this->curSol()[domainJ]; | |
| 910 | |||
| 911 | // build derivatives with for cell center dofs w.r.t. cell center dofs | ||
| 912 |
4/4✓ Branch 0 taken 1651648 times.
✓ Branch 1 taken 412912 times.
✓ Branch 2 taken 1651648 times.
✓ Branch 3 taken 412912 times.
|
2477472 | for (auto&& scvf : scvfs(fvGeometry)) |
| 913 | { | ||
| 914 | // set the actual dof index | ||
| 915 | 1651648 | const auto faceGlobalI = scvf.dofIndex(); | |
| 916 | |||
| 917 | // get stencil information | ||
| 918 | 1651648 | const auto& stencil = this->couplingManager().couplingStencil(domainI, scvf, domainJ); | |
| 919 | |||
| 920 |
4/4✓ Branch 0 taken 11332 times.
✓ Branch 1 taken 1640316 times.
✓ Branch 2 taken 11332 times.
✓ Branch 3 taken 1640316 times.
|
3303296 | if (stencil.empty()) |
| 921 | continue; | ||
| 922 | |||
| 923 | // build derivatives with for face dofs w.r.t. cell center dofs | ||
| 924 |
4/4✓ Branch 0 taken 11332 times.
✓ Branch 1 taken 11332 times.
✓ Branch 2 taken 11332 times.
✓ Branch 3 taken 11332 times.
|
33996 | for (const auto& globalJ : stencil) |
| 925 | { | ||
| 926 |
1/2✓ Branch 0 taken 11332 times.
✗ Branch 1 not taken.
|
11332 | const auto origPriVarsJ = curSol[globalJ]; |
| 927 |
1/2✓ Branch 0 taken 11332 times.
✗ Branch 1 not taken.
|
11332 | const auto origResidual = this->couplingManager().evalCouplingResidual(domainI, scvf, *this, domainJ, globalJ); |
| 928 | |||
| 929 |
2/2✓ Branch 0 taken 22896 times.
✓ Branch 1 taken 11332 times.
|
34228 | for (int pvIdx = 0; pvIdx < JacobianBlock::block_type::cols; ++pvIdx) |
| 930 | { | ||
| 931 | 112440 | auto evalCouplingResidual = [&](Scalar priVar) | |
| 932 | { | ||
| 933 | 20856 | auto deflectedPriVars = origPriVarsJ; | |
| 934 | 22896 | deflectedPriVars[pvIdx] = priVar; | |
| 935 | 22896 | this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, deflectedPriVars, pvIdx); | |
| 936 |
1/2✓ Branch 0 taken 22896 times.
✗ Branch 1 not taken.
|
22896 | return this->couplingManager().evalCouplingResidual(domainI, scvf, *this, domainJ, globalJ); |
| 937 | }; | ||
| 938 | |||
| 939 | // derive the residuals numerically | ||
| 940 |
2/2✓ Branch 0 taken 17 times.
✓ Branch 1 taken 22879 times.
|
22896 | FaceResidualValue partialDeriv(0.0); |
| 941 |
4/4✓ Branch 0 taken 17 times.
✓ Branch 1 taken 22879 times.
✓ Branch 2 taken 17 times.
✓ Branch 3 taken 22879 times.
|
45792 | const auto& paramGroup = this->assembler().problem(domainJ).paramGroup(); |
| 942 |
4/6✓ Branch 0 taken 17 times.
✓ Branch 1 taken 22879 times.
✓ Branch 3 taken 17 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 17 times.
✗ Branch 7 not taken.
|
22896 | static const int numDiffMethod = getParamFromGroup<int>(paramGroup, "Assembly.NumericDifferenceMethod"); |
| 943 |
4/6✓ Branch 0 taken 17 times.
✓ Branch 1 taken 22879 times.
✓ Branch 3 taken 17 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 17 times.
✗ Branch 7 not taken.
|
22896 | static const auto epsCoupl = this->couplingManager().numericEpsilon(domainJ, paramGroup); |
| 944 |
2/4✓ Branch 1 taken 22896 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20856 times.
✗ Branch 5 not taken.
|
43752 | NumericDifferentiation::partialDerivative(evalCouplingResidual, origPriVarsJ[pvIdx], partialDeriv, origResidual, |
| 945 | 43752 | epsCoupl(origPriVarsJ[pvIdx], pvIdx), numDiffMethod); | |
| 946 | |||
| 947 | // update the global stiffness matrix with the current partial derivatives | ||
| 948 |
1/2✓ Branch 1 taken 22896 times.
✗ Branch 2 not taken.
|
22896 | updateGlobalJacobian_(A, faceGlobalI, globalJ, pvIdx, partialDeriv); |
| 949 | |||
| 950 | // restore the undeflected state of the coupling context | ||
| 951 |
1/2✓ Branch 1 taken 22896 times.
✗ Branch 2 not taken.
|
22896 | this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, origPriVarsJ, pvIdx); |
| 952 | } | ||
| 953 | } | ||
| 954 | } | ||
| 955 | 412912 | } | |
| 956 | |||
| 957 | template<class JacobianMatrixDiagBlock, class GridVariables> | ||
| 958 | void evalAdditionalDerivatives(const std::vector<std::size_t>& additionalDofDependencies, | ||
| 959 | JacobianMatrixDiagBlock& A, GridVariables& gridVariables) | ||
| 960 | { } | ||
| 961 | |||
| 962 | /*! | ||
| 963 | * \brief Updates the current global Jacobian matrix with the | ||
| 964 | * partial derivatives of all equations in regard to the | ||
| 965 | * primary variable 'pvIdx' at dof 'col'. Specialization for cc methods. | ||
| 966 | */ | ||
| 967 | template<class SubMatrix, class CCOrFacePrimaryVariables> | ||
| 968 | 214559020 | static void updateGlobalJacobian_(SubMatrix& matrix, | |
| 969 | const int globalI, | ||
| 970 | const int globalJ, | ||
| 971 | const int pvIdx, | ||
| 972 | const CCOrFacePrimaryVariables& partialDeriv) | ||
| 973 | { | ||
| 974 |
2/2✓ Branch 0 taken 154845250 times.
✓ Branch 1 taken 107279510 times.
|
524249520 | for (int eqIdx = 0; eqIdx < partialDeriv.size(); eqIdx++) |
| 975 | { | ||
| 976 | // A[i][col][eqIdx][pvIdx] is the rate of change of | ||
| 977 | // the residual of equation 'eqIdx' at dof 'i' | ||
| 978 | // depending on the primary variable 'pvIdx' at dof | ||
| 979 | // 'col'. | ||
| 980 | |||
| 981 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 154845250 times.
|
309690500 | assert(pvIdx >= 0); |
| 982 | 619381000 | assert(eqIdx < matrix[globalI][globalJ].size()); | |
| 983 |
2/4✗ Branch 2 not taken.
✓ Branch 3 taken 154845250 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 154845250 times.
|
619381000 | assert(pvIdx < matrix[globalI][globalJ][eqIdx].size()); |
| 984 | 753560548 | matrix[globalI][globalJ][eqIdx][pvIdx] += partialDeriv[eqIdx]; | |
| 985 | } | ||
| 986 | 214559020 | } | |
| 987 | |||
| 988 | private: | ||
| 989 | |||
| 990 | ✗ | FaceVariables& getFaceVarAccess_(GridFaceVariables& gridFaceVariables, ElementFaceVariables& elemFaceVars, const SubControlVolumeFace& scvf) | |
| 991 | { | ||
| 992 | if constexpr (GetPropType<TypeTag, Properties::GridVariables>::GridFaceVariables::cachingEnabled) | ||
| 993 | 96523728 | return gridFaceVariables.faceVars(scvf.index()); | |
| 994 | else | ||
| 995 | return elemFaceVars[scvf]; | ||
| 996 | } | ||
| 997 | }; | ||
| 998 | |||
| 999 | } // end namespace Dumux | ||
| 1000 | |||
| 1001 | #endif | ||
| 1002 |