GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/multidomain/embedded/extendedsourcestencil.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 45 45 100.0%
Functions: 44 44 100.0%
Branches: 71 98 72.4%
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 EmbeddedCoupling
10 * \brief Extended source stencil helper class for coupling managers
11 */
12
13 #ifndef DUMUX_MULTIDOMAIN_EMBEDDED_EXTENDEDSOURCESTENCIL_HH
14 #define DUMUX_MULTIDOMAIN_EMBEDDED_EXTENDEDSOURCESTENCIL_HH
15
16 #include <vector>
17
18 #include <dune/common/indices.hh>
19
20 #include <dumux/common/properties.hh>
21 #include <dumux/common/parameters.hh>
22 #include <dumux/common/numericdifferentiation.hh>
23 #include <dumux/assembly/numericepsilon.hh>
24
25 #include <dumux/discretization/method.hh>
26
27 namespace Dumux::EmbeddedCoupling {
28
29 /*!
30 * \ingroup EmbeddedCoupling
31 * \brief A class managing an extended source stencil
32 * \tparam CouplingManager the coupling manager type
33 */
34 template<class CouplingManager>
35
2/4
✓ Branch 0 taken 17 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 17 times.
✗ Branch 3 not taken.
 34 class ExtendedSourceStencil
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 bulkIdx = typename MDTraits::template SubDomain<0>::Index();
46 static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index();
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 92 void extendJacobianPattern(const CouplingManager& couplingManager, Dune::index_constant<id> domainI, JacobianPattern& pattern) const
60 {
61 // add additional dof dependencies
62
6/8
✗ Branch 0 not taken.
✓ Branch 1 taken 46 times.
✓ Branch 3 taken 954038 times.
✓ Branch 4 taken 24 times.
✓ Branch 5 taken 11662 times.
✓ Branch 6 taken 24 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 8424 times.
 3792832 for (const auto& element : elements(couplingManager.gridView(domainI)))
63 {
64 1924800 const auto& dofs = extendedSourceStencil_(couplingManager, domainI, element);
65 if constexpr (isBox<domainI>())
66 {
67
5/5
✓ Branch 0 taken 80 times.
✓ Branch 1 taken 123608 times.
✓ Branch 2 taken 15526 times.
✓ Branch 3 taken 123608 times.
✓ Branch 4 taken 15446 times.
 278348 for (int i = 0; i < element.subEntities(GridView<domainI>::dimension); ++i)
68
4/4
✓ Branch 0 taken 45776 times.
✓ Branch 1 taken 123648 times.
✓ Branch 2 taken 45776 times.
✓ Branch 3 taken 123648 times.
 924992 for (const auto globalJ : dofs)
69
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 45776 times.
 91552 pattern.add(couplingManager.problem(domainI).gridGeometry().vertexMapper().subIndex(element, i, GridView<domainI>::dimension), globalJ);
70 }
71 else
72 {
73
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 930106 times.
 1910756 const auto globalI = couplingManager.problem(domainI).gridGeometry().elementMapper().index(element);
74
4/4
✓ Branch 0 taken 67718 times.
✓ Branch 1 taken 938530 times.
✓ Branch 2 taken 67718 times.
✓ Branch 3 taken 938530 times.
 5902052 for (const auto globalJ : dofs)
75 135436 pattern.add(globalI, globalJ);
76 }
77 }
78 92 }
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 JacobianMatrixDiagBlock, class GridVariables>
88 13157336 void evalAdditionalDomainDerivatives(CouplingManager& couplingManager,
89 Dune::index_constant<i> domainI,
90 const LocalAssemblerI& localAssemblerI,
91 JacobianMatrixDiagBlock& A,
92 GridVariables& gridVariables) const
93 {
94
1/2
✓ Branch 0 taken 856990 times.
✗ Branch 1 not taken.
 13157336 const auto& curSolI = couplingManager.curSol(domainI);
95 13157336 constexpr auto numEq = std::decay_t<decltype(curSolI[0])>::size();
96 13157336 const auto& elementI = localAssemblerI.element();
97
98 // only do something if we have an extended stencil
99
5/5
✓ Branch 0 taken 856990 times.
✓ Branch 1 taken 5516791 times.
✓ Branch 2 taken 204887 times.
✓ Branch 3 taken 5516791 times.
✓ Branch 4 taken 204887 times.
 13157336 if (extendedSourceStencil_(couplingManager, domainI, elementI).empty())
100 12747562 return;
101
102 // compute the undeflected residual (source only!)
103 409774 const auto origResidual = localAssemblerI.evalLocalSourceResidual(elementI);
104
105 // compute derivate for all additional dofs in the circle stencil
106
4/5
✗ Branch 0 not taken.
✓ Branch 1 taken 654492 times.
✓ Branch 2 taken 204887 times.
✓ Branch 3 taken 654492 times.
✓ Branch 4 taken 204887 times.
 1718758 for (const auto dofIndex : extendedSourceStencil_(couplingManager, domainI, elementI))
107 {
108 1308984 auto partialDerivs = origResidual;
109 1308984 const auto origPriVars = curSolI[dofIndex];
110 1308984 auto priVars = origPriVars;
111
112 // calculate derivatives w.r.t to the privars at the dof at hand
113
2/2
✓ Branch 0 taken 684582 times.
✓ Branch 1 taken 654492 times.
 2678148 for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
114 {
115 // reset partial derivatives
116 1369164 partialDerivs = 0.0;
117
118 2738328 const auto evalResiduals = [&](const Scalar priVar)
119 {
120 // update the coupling context (solution vector and recompute element residual)
121 1308984 priVars[pvIdx] = priVar;
122 684582 couplingManager.updateCouplingContext(domainI, localAssemblerI, domainI, dofIndex, priVars, pvIdx);
123 684582 return localAssemblerI.evalLocalSourceResidual(elementI);
124 };
125
126 // derive the residuals numerically
127
7/10
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 684557 times.
✓ Branch 3 taken 23 times.
✓ Branch 4 taken 2 times.
✓ Branch 6 taken 23 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 23 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 23 times.
✗ Branch 13 not taken.
 1369164 static const NumericEpsilon<Scalar, numEq> eps_{localAssemblerI.problem().paramGroup()};
128
7/10
✓ Branch 0 taken 25 times.
✓ Branch 1 taken 684557 times.
✓ Branch 3 taken 23 times.
✓ Branch 4 taken 2 times.
✓ Branch 6 taken 23 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 23 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 23 times.
✗ Branch 13 not taken.
 1369164 static const int numDiffMethod = getParamFromGroup<int>(localAssemblerI.problem().paramGroup(), "Assembly.NumericDifferenceMethod");
129 4107492 NumericDifferentiation::partialDerivative(
130
3/5
✓ Branch 2 taken 624402 times.
✓ Branch 3 taken 60180 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 60180 times.
✗ Branch 7 not taken.
 1489524 evalResiduals, priVars[pvIdx], partialDerivs, origResidual, eps_(priVars[pvIdx], pvIdx), numDiffMethod
131 );
132
133 // update the global stiffness matrix with the current partial derivatives
134
4/4
✓ Branch 0 taken 1116622 times.
✓ Branch 1 taken 684582 times.
✓ Branch 2 taken 1116622 times.
✓ Branch 3 taken 684582 times.
 6340736 for (const auto& scvJ : scvs(localAssemblerI.fvGeometry()))
135 {
136
2/2
✓ Branch 0 taken 1176802 times.
✓ Branch 1 taken 1116622 times.
 4586848 for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
137 {
138 // A[i][col][eqIdx][pvIdx] is the rate of change of
139 // the residual of equation 'eqIdx' at dof 'i'
140 // depending on the primary variable 'pvIdx' at dof
141 // 'col'.
142
5/10
✓ Branch 1 taken 1176802 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1176802 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1176802 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 683042 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 120360 times.
✗ Branch 14 not taken.
 8667616 A[scvJ.dofIndex()][dofIndex][eqIdx][pvIdx] += partialDerivs[scvJ.indexInElement()][eqIdx];
143 }
144 }
145
146 // restore the current element solution
147 1489524 priVars[pvIdx] = origPriVars[pvIdx];
148
149 // restore the original coupling context
150 2738328 couplingManager.updateCouplingContext(domainI, localAssemblerI, domainI, dofIndex, origPriVars, pvIdx);
151 }
152 }
153 }
154
155 //! clear the internal data
156 12 void clear() { sourceStencils_.clear(); }
157
158 //! return a reference to the stencil
159 typename CouplingManager::template CouplingStencils<bulkIdx>& stencil()
160
4/6
✓ Branch 1 taken 67865 times.
✓ Branch 2 taken 48640 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 67865 times.
✓ Branch 5 taken 48640 times.
✗ Branch 6 not taken.
 443517 { return sourceStencils_; }
161
162 //! return a const reference to the stencil
163 const typename CouplingManager::template CouplingStencils<bulkIdx>& stencil() const
164 1875744 { return sourceStencils_; }
165
166 private:
167 //! the extended source stencil due to the source average (always empty for lowdim, but may be filled for bulk)
168 template<std::size_t id>
169 6872117 const auto& extendedSourceStencil_(const CouplingManager& couplingManager, Dune::index_constant<id> domainId, const Element<id>& element) const
170 {
171 if constexpr (domainId == bulkIdx)
172 {
173 6872117 const auto bulkElementIdx = couplingManager.problem(bulkIdx).gridGeometry().elementMapper().index(element);
174 20616351 if (auto stencil = sourceStencils_.find(bulkElementIdx); stencil != sourceStencils_.end())
175 10863430 return stencil->second;
176 }
177
178
4/8
✓ Branch 0 taken 856990 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 856990 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 8424 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 8424 times.
 4611632 return couplingManager.emptyStencil(domainId);
179 }
180
181 //! the additional stencil for the kernel evaluations / circle averages
182 typename CouplingManager::template CouplingStencils<bulkIdx> sourceStencils_;
183 };
184
185 } // end namespace Dumux::EmbeddedCoupling
186
187 #endif
188