| 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 DualNetworkCoupling | ||
| 10 | * \ingroup PoreNetworkModels | ||
| 11 | * \brief Extended source stencil helper class for coupling managers | ||
| 12 | */ | ||
| 13 | |||
| 14 | #ifndef DUMUX_DUAL_NETWORK_EXTENDEDSOURCESTENCIL_HH | ||
| 15 | #define DUMUX_DUAL_NETWORK_EXTENDEDSOURCESTENCIL_HH | ||
| 16 | |||
| 17 | #include <vector> | ||
| 18 | |||
| 19 | #include <dune/common/indices.hh> | ||
| 20 | |||
| 21 | #include <dumux/common/properties.hh> | ||
| 22 | #include <dumux/common/parameters.hh> | ||
| 23 | #include <dumux/common/numericdifferentiation.hh> | ||
| 24 | #include <dumux/discretization/method.hh> | ||
| 25 | |||
| 26 | namespace Dumux::PoreNetwork { | ||
| 27 | |||
| 28 | /*! | ||
| 29 | * \ingroup DualNetworkCoupling | ||
| 30 | * \ingroup PoreNetworkModels | ||
| 31 | * \brief A class managing an extended source stencil | ||
| 32 | * \tparam CouplingManager the coupling manager type | ||
| 33 | */ | ||
| 34 | template<class CouplingManager> | ||
| 35 | class PNMHeatExtendedSourceStencil | ||
| 36 | { | ||
| 37 | using MDTraits = typename CouplingManager::MultiDomainTraits; | ||
| 38 | using Scalar = typename MDTraits::Scalar; | ||
| 39 | |||
| 40 | template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag; | ||
| 41 | template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>; | ||
| 42 | template<std::size_t id> using GridView = typename GridGeometry<id>::GridView; | ||
| 43 | template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity; | ||
| 44 | |||
| 45 | static constexpr auto solidDomainIdx = CouplingManager::solidDomainIdx; | ||
| 46 | static constexpr auto voidDomainIdx = CouplingManager::voidDomainIdx; | ||
| 47 | |||
| 48 | template<std::size_t id> | ||
| 49 | static constexpr bool isBox() | ||
| 50 | { return GridGeometry<id>::discMethod == DiscretizationMethods::box; } | ||
| 51 | public: | ||
| 52 | /*! | ||
| 53 | * \brief extend the jacobian pattern of the diagonal block of domain i | ||
| 54 | * by those entries that are not already in the uncoupled pattern | ||
| 55 | * \note Such additional dependencies can arise from the coupling, e.g. if a coupling source | ||
| 56 | * term depends on a non-local average of a quantity of the same domain | ||
| 57 | */ | ||
| 58 | template<std::size_t id, class JacobianPattern> | ||
| 59 | 4 | void extendJacobianPattern(const CouplingManager& couplingManager, Dune::index_constant<id> domainI, JacobianPattern& pattern) const | |
| 60 | { | ||
| 61 | // add additional dof dependencies | ||
| 62 |
4/4✓ Branch 2 taken 3191 times.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 3191 times.
✓ Branch 5 taken 2 times.
|
6390 | for (const auto& element : elements(couplingManager.gridView(domainI))) |
| 63 | { | ||
| 64 | 10750 | const auto& dofs = extendedSourceStencil_(couplingManager, domainI, element); | |
| 65 | if constexpr (isBox<domainI>()) | ||
| 66 | { | ||
| 67 |
4/4✓ Branch 0 taken 6382 times.
✓ Branch 1 taken 3191 times.
✓ Branch 2 taken 6382 times.
✓ Branch 3 taken 3191 times.
|
25528 | for (int i = 0; i < element.subEntities(GridView<domainI>::dimension); ++i) |
| 68 |
4/4✓ Branch 0 taken 17572 times.
✓ Branch 1 taken 6382 times.
✓ Branch 2 taken 17572 times.
✓ Branch 3 taken 6382 times.
|
108580 | for (const auto globalJ : dofs) |
| 69 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 17572 times.
|
35144 | pattern.add(couplingManager.problem(domainI).gridGeometry().vertexMapper().subIndex(element, i, GridView<domainI>::dimension), globalJ); |
| 70 | } | ||
| 71 | else | ||
| 72 | { | ||
| 73 | const auto globalI = couplingManager.problem(domainI).gridGeometry().elementMapper().index(element); | ||
| 74 | for (const auto globalJ : dofs) | ||
| 75 | pattern.add(globalI, globalJ); | ||
| 76 | } | ||
| 77 | } | ||
| 78 | 4 | } | |
| 79 | |||
| 80 | /*! | ||
| 81 | * \brief evaluate additional derivatives of the element residual of a domain with respect | ||
| 82 | * to dofs in the same domain that are not in the regular stencil (per default this is not the case) | ||
| 83 | * \note Such additional dependencies can arise from the coupling, e.g. if a coupling source | ||
| 84 | * term depends on a non-local average of a quantity of the same domain | ||
| 85 | * \note This is the same for box and cc | ||
| 86 | */ | ||
| 87 | template<std::size_t i, class LocalAssemblerI, class SolutionVector, class JacobianMatrixDiagBlock, class GridVariables> //TODO edit function s.t. curSol from couplingmanager is used like in dumux/dumux/multidomain/embedded/extendedsourcestencil.hh (commit 3257876) | ||
| 88 | 25528 | void evalAdditionalDomainDerivatives(CouplingManager& couplingManager, | |
| 89 | Dune::index_constant<i> domainI, | ||
| 90 | const LocalAssemblerI& localAssemblerI, | ||
| 91 | const SolutionVector& curSol, | ||
| 92 | JacobianMatrixDiagBlock& A, | ||
| 93 | GridVariables& gridVariables) const | ||
| 94 | { | ||
| 95 | // const auto& curSolI = couplingManager.curSol(domainI); //TODO:curSol declared protected in .../dumux/dumux/multidomain/couplingmanager.hh | ||
| 96 | 25528 | const auto& curSolI = curSol; | |
| 97 | 25528 | constexpr auto numEq = std::decay_t<decltype(curSolI[0])>::size(); | |
| 98 | 25528 | const auto& elementI = localAssemblerI.element(); | |
| 99 | |||
| 100 | // only do something if we have an extended stencil | ||
| 101 |
3/5✓ Branch 0 taken 8736 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 4028 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 4028 times.
|
25528 | if (extendedSourceStencil_(couplingManager, domainI, elementI).empty()) |
| 102 | 17472 | return; | |
| 103 | |||
| 104 | // compute the undeflected residual (source only!) | ||
| 105 | 8056 | const auto origResidual = localAssemblerI.evalLocalSourceResidual(elementI); | |
| 106 | |||
| 107 | // compute derivate for all additional dofs in the circle stencil | ||
| 108 |
4/5✗ Branch 0 not taken.
✓ Branch 1 taken 35144 times.
✓ Branch 2 taken 4028 times.
✓ Branch 3 taken 35144 times.
✓ Branch 4 taken 4028 times.
|
78344 | for (const auto dofIndex : extendedSourceStencil_(couplingManager, domainI, elementI)) |
| 109 | { | ||
| 110 | // std::cout << "deflecting " << dofIndex << std::endl; | ||
| 111 | 70288 | auto partialDerivs = origResidual; | |
| 112 | 70288 | const auto origPriVars = curSolI[dofIndex]; | |
| 113 | |||
| 114 | // calculate derivatives w.r.t to the privars at the dof at hand | ||
| 115 |
2/2✓ Branch 0 taken 70288 times.
✓ Branch 1 taken 35144 times.
|
210864 | for (int pvIdx = 0; pvIdx < numEq; pvIdx++) |
| 116 | { | ||
| 117 | // reset partial derivatives | ||
| 118 | 140576 | partialDerivs = 0.0; | |
| 119 | |||
| 120 | 281152 | auto evalResiduals = [&](Scalar priVar) | |
| 121 | { | ||
| 122 | // update the coupling context (solution vector and recompute element residual) | ||
| 123 | 70288 | auto priVars = origPriVars; | |
| 124 | 70288 | priVars[pvIdx] = priVar; | |
| 125 | 70288 | couplingManager.updateCouplingContext(domainI, localAssemblerI, domainI, dofIndex, priVars, pvIdx); | |
| 126 | 70288 | return localAssemblerI.evalLocalSourceResidual(elementI); | |
| 127 | }; | ||
| 128 | |||
| 129 | // derive the residuals numerically | ||
| 130 |
4/6✓ Branch 0 taken 1 times.
✓ Branch 1 taken 70287 times.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
|
140576 | static const int numDiffMethod = getParam<int>("Assembly.NumericDifferenceMethod"); |
| 131 |
3/6✓ Branch 1 taken 70288 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 70288 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 70288 times.
✗ Branch 8 not taken.
|
421728 | NumericDifferentiation::partialDerivative(evalResiduals, curSolI[dofIndex][pvIdx], |
| 132 | partialDerivs, origResidual, numDiffMethod); //TODO: maybe change like in dumux/dumux/md/embedded/extendedsourcestencil.hh line 142 | ||
| 133 | |||
| 134 | // update the global stiffness matrix with the current partial derivatives | ||
| 135 |
2/2✓ Branch 0 taken 140576 times.
✓ Branch 1 taken 70288 times.
|
702880 | for (const auto& scvJ : scvs(localAssemblerI.fvGeometry())) |
| 136 | { | ||
| 137 |
2/2✓ Branch 0 taken 281152 times.
✓ Branch 1 taken 140576 times.
|
843456 | for (int eqIdx = 0; eqIdx < numEq; eqIdx++) |
| 138 | { | ||
| 139 | // A[i][col][eqIdx][pvIdx] is the rate of change of | ||
| 140 | // the residual of equation 'eqIdx' at dof 'i' | ||
| 141 | // depending on the primary variable 'pvIdx' at dof | ||
| 142 | // 'col'. | ||
| 143 | // std::cout << "Adding " << partialDerivs[scvJ.indexInElement()][eqIdx] << ", pvIdx " << pvIdx << std::endl; | ||
| 144 |
4/8✓ Branch 1 taken 281152 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 281152 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 281152 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 281152 times.
✗ Branch 11 not taken.
|
2249216 | A[scvJ.dofIndex()][dofIndex][eqIdx][pvIdx] += partialDerivs[scvJ.indexInElement()][eqIdx]; |
| 145 | } | ||
| 146 | } | ||
| 147 | |||
| 148 | // restore the original coupling context | ||
| 149 |
1/2✓ Branch 1 taken 70288 times.
✗ Branch 2 not taken.
|
140576 | couplingManager.updateCouplingContext(domainI, localAssemblerI, domainI, dofIndex, origPriVars, pvIdx); |
| 150 | } | ||
| 151 | } | ||
| 152 | } | ||
| 153 | |||
| 154 | //! clear the internal data | ||
| 155 | 1 | void clear() { sourceStencils_.clear(); } | |
| 156 | |||
| 157 | //! return a reference to the stencil | ||
| 158 | auto& stencil() | ||
| 159 |
1/2✓ Branch 1 taken 9632 times.
✗ Branch 2 not taken.
|
10639 | { return sourceStencils_; } |
| 160 | |||
| 161 | private: | ||
| 162 | //! the extended source stencil due to the source average (always empty for lowdim, but may be filled for bulk) | ||
| 163 | template<std::size_t id> | ||
| 164 | ✗ | const auto& extendedSourceStencil_(const CouplingManager& couplingManager, Dune::index_constant<id> domainId, const Element<id>& element) const | |
| 165 | { | ||
| 166 | if constexpr (domainId == voidDomainIdx) | ||
| 167 | { | ||
| 168 | ✗ | const auto voidElementIdx = couplingManager.problem(voidDomainIdx).gridGeometry().elementMapper().index(element); | |
| 169 | ✗ | if (sourceStencils_.count(voidElementIdx)) | |
| 170 | ✗ | return sourceStencils_.at(voidElementIdx); | |
| 171 | } | ||
| 172 | |||
| 173 |
2/4✓ Branch 0 taken 8736 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 8736 times.
✗ Branch 3 not taken.
|
21840 | return couplingManager.emptyStencil(domainId); |
| 174 | } | ||
| 175 | |||
| 176 | //! the additional stencil for the kernel evaluations / circle averages | ||
| 177 | std::unordered_map<std::size_t, std::vector<std::size_t>> sourceStencils_; | ||
| 178 | }; | ||
| 179 | |||
| 180 | } // end namespace Dumux::PoreNetwork | ||
| 181 | |||
| 182 | #endif | ||
| 183 |