extendedsourcestencil.hh 8.99 KB
Newer Older
1
2
3
4
5
6
7
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
8
 *   the Free Software Foundation, either version 3 of the License, or       *
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
/*!
 * \file
 * \ingroup EmbeddedCoupling
 * \brief Extended source stencil helper class for coupling managers
 */

#ifndef DUMUX_MULTIDOMAIN_EMBEDDED_EXTENDEDSOURCESTENCIL_HH
#define DUMUX_MULTIDOMAIN_EMBEDDED_EXTENDEDSOURCESTENCIL_HH

#include <vector>
Timo Koch's avatar
Timo Koch committed
29
30
31

#include <dune/common/indices.hh>

32
#include <dumux/common/properties.hh>
Timo Koch's avatar
Timo Koch committed
33
34
35
#include <dumux/common/parameters.hh>
#include <dumux/common/numericdifferentiation.hh>
#include <dumux/discretization/method.hh>
36

37
namespace Dumux::EmbeddedCoupling {
38
39
40
41
42
43
44
45
46
47
48
49

/*!
 * \ingroup EmbeddedCoupling
 * \brief A class managing an extended source stencil
 * \tparam CouplingManager the coupling manager type
 */
template<class CouplingManager>
class ExtendedSourceStencil
{
    using MDTraits = typename CouplingManager::MultiDomainTraits;
    using Scalar = typename MDTraits::Scalar;

50
    template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
51
52
    template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
    template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
53
54
    template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;

55
    static constexpr auto bulkIdx = typename MDTraits::template SubDomain<0>::Index();
Timo Koch's avatar
Timo Koch committed
56
    static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index();
57
58
59

    template<std::size_t id>
    static constexpr bool isBox()
60
    { return GridGeometry<id>::discMethod == DiscretizationMethod::box; }
61
62
63
64
65
66
67
68
69
70
71
72
73
74
public:
    /*!
     * \brief extend the jacobian pattern of the diagonal block of domain i
     *        by those entries that are not already in the uncoupled pattern
     * \note Such additional dependencies can arise from the coupling, e.g. if a coupling source
     *       term depends on a non-local average of a quantity of the same domain
     */
    template<std::size_t id, class JacobianPattern>
    void extendJacobianPattern(const CouplingManager& couplingManager, Dune::index_constant<id> domainI, JacobianPattern& pattern) const
    {
        // add additional dof dependencies
        for (const auto& element : elements(couplingManager.gridView(domainI)))
        {
            const auto& dofs = extendedSourceStencil_(couplingManager, domainI, element);
75
            if constexpr (isBox<domainI>())
76
77
78
            {
                for (int i = 0; i < element.subEntities(GridView<domainI>::dimension); ++i)
                    for (const auto globalJ : dofs)
79
                        pattern.add(couplingManager.problem(domainI).gridGeometry().vertexMapper().subIndex(element, i, GridView<domainI>::dimension), globalJ);
80
81
82
            }
            else
            {
83
                const auto globalI = couplingManager.problem(domainI).gridGeometry().elementMapper().index(element);
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
                for (const auto globalJ : dofs)
                    pattern.add(globalI, globalJ);
            }
        }
    }

    /*!
     * \brief evaluate additional derivatives of the element residual of a domain with respect
     *        to dofs in the same domain that are not in the regular stencil (per default this is not the case)
     * \note Such additional dependencies can arise from the coupling, e.g. if a coupling source
     *       term depends on a non-local average of a quantity of the same domain
     * \note This is the same for box and cc
     */
    template<std::size_t i, class LocalAssemblerI, class SolutionVector, class JacobianMatrixDiagBlock, class GridVariables>
    void evalAdditionalDomainDerivatives(CouplingManager& couplingManager,
                                         Dune::index_constant<i> domainI,
                                         const LocalAssemblerI& localAssemblerI,
                                         const SolutionVector& curSol,
                                         JacobianMatrixDiagBlock& A,
                                         GridVariables& gridVariables) const
    {
105
        constexpr auto numEq = std::decay_t<decltype(curSol[domainI][0])>::size();
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
        const auto& elementI = localAssemblerI.element();

        // only do something if we have an extended stencil
        if (extendedSourceStencil_(couplingManager, domainI, elementI).empty())
            return;

        // compute the undeflected residual (source only!)
        const auto origResidual = localAssemblerI.evalLocalSourceResidual(elementI);

        // compute derivate for all additional dofs in the circle stencil
        for (const auto dofIndex : extendedSourceStencil_(couplingManager, domainI, elementI))
        {
            auto partialDerivs = origResidual;
            const auto origPriVars = curSol[domainI][dofIndex];

            // calculate derivatives w.r.t to the privars at the dof at hand
            for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
            {
                // reset partial derivatives
                partialDerivs = 0.0;

                auto evalResiduals = [&](Scalar priVar)
                {
                    // update the coupling context (solution vector and recompute element residual)
                    auto priVars = origPriVars;
                    priVars[pvIdx] = priVar;
                    couplingManager.updateCouplingContext(domainI, localAssemblerI, domainI, dofIndex, priVars, pvIdx);
                    return localAssemblerI.evalLocalSourceResidual(elementI);
                };

                // derive the residuals numerically
                static const int numDiffMethod = getParam<int>("Assembly.NumericDifferenceMethod");
                NumericDifferentiation::partialDerivative(evalResiduals, curSol[domainI][dofIndex][pvIdx],
                                                          partialDerivs, origResidual, numDiffMethod);

                // update the global stiffness matrix with the current partial derivatives
142
                for (const auto& scvJ : scvs(localAssemblerI.fvGeometry()))
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
                {
                    for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
                    {
                        // A[i][col][eqIdx][pvIdx] is the rate of change of
                        // the residual of equation 'eqIdx' at dof 'i'
                        // depending on the primary variable 'pvIdx' at dof
                        // 'col'.
                        A[scvJ.dofIndex()][dofIndex][eqIdx][pvIdx] += partialDerivs[scvJ.indexInElement()][eqIdx];
                    }
                }

                // restore the original coupling context
                couplingManager.updateCouplingContext(domainI, localAssemblerI, domainI, dofIndex, origPriVars, pvIdx);
            }
        }
    }

160
161
162
    //! clear the internal data
    void clear() { sourceStencils_.clear(); }

163
    //! return a reference to the stencil
164
    typename CouplingManager::template CouplingStencils<bulkIdx>& stencil()
165
166
167
    { return sourceStencils_; }

private:
168
169
    //! the extended source stencil due to the source average (always empty for lowdim, but may be filled for bulk)
    template<std::size_t id>
170
    const auto& extendedSourceStencil_(const CouplingManager& couplingManager, Dune::index_constant<id> domainId, const Element<id>& element) const
171
    {
Timo Koch's avatar
Timo Koch committed
172
        if constexpr (domainId == bulkIdx)
173
        {
174
            const auto bulkElementIdx = couplingManager.problem(bulkIdx).gridGeometry().elementMapper().index(element);
175
176
            if (auto stencil = sourceStencils_.find(bulkElementIdx); stencil != sourceStencils_.end())
                return stencil->second;
177
        }
178
179
        
        return couplingManager.emptyStencil(domainId);
180
    }
181
182

    //! the additional stencil for the kernel evaluations / circle averages
183
    typename CouplingManager::template CouplingStencils<bulkIdx> sourceStencils_;
184
185
};

186
} // end namespace Dumux::EmbeddedCoupling
187
188

#endif