| 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 Data associated with a point source | ||
| 11 | */ | ||
| 12 | |||
| 13 | #ifndef DUMUX_MULTIDOMAIN_EMBEDDED_POINTSOURCEDATA_HH | ||
| 14 | #define DUMUX_MULTIDOMAIN_EMBEDDED_POINTSOURCEDATA_HH | ||
| 15 | |||
| 16 | #include <vector> | ||
| 17 | #include <dune/common/fvector.hh> | ||
| 18 | #include <dumux/common/properties.hh> | ||
| 19 | #include <dumux/common/indextraits.hh> | ||
| 20 | #include <dumux/discretization/method.hh> | ||
| 21 | |||
| 22 | namespace Dumux { | ||
| 23 | |||
| 24 | /*! | ||
| 25 | * \ingroup EmbeddedCoupling | ||
| 26 | * \brief A point source data class used for integration in multidimensional models | ||
| 27 | * \note The point source and related data are connected via an identifier (id) | ||
| 28 | */ | ||
| 29 | template<class MDTraits> | ||
| 30 |
7/12✓ Branch 0 taken 48640 times.
✓ Branch 1 taken 37019 times.
✓ Branch 2 taken 48640 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 85659 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 37019 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 26709 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 26709 times.
✗ Branch 14 not taken.
|
310395 | class PointSourceData |
| 31 | { | ||
| 32 | using Scalar = typename MDTraits::Scalar; | ||
| 33 | using ShapeValues = typename std::vector<Dune::FieldVector<Scalar, 1> >; | ||
| 34 | |||
| 35 | template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag; | ||
| 36 | template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>; | ||
| 37 | template<std::size_t id> using GridView = typename GridGeometry<id>::GridView; | ||
| 38 | template<std::size_t id> using GridIndex = typename IndexTraits<GridView<id>>::GridIndex; | ||
| 39 | template<std::size_t id> using SolutionVector = GetPropType<SubDomainTypeTag<id>, Properties::SolutionVector>; | ||
| 40 | template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>; | ||
| 41 | |||
| 42 | static constexpr auto bulkIdx = typename MDTraits::template SubDomain<0>::Index(); | ||
| 43 | static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index(); | ||
| 44 | |||
| 45 | template<std::size_t id> | ||
| 46 | static constexpr bool isBox() | ||
| 47 | { return GridGeometry<id>::discMethod == DiscretizationMethods::box; } | ||
| 48 | |||
| 49 | public: | ||
| 50 | void addBulkInterpolation(const ShapeValues& shapeValues, | ||
| 51 | const std::vector<GridIndex<bulkIdx>>& cornerIndices, | ||
| 52 | GridIndex<bulkIdx> eIdx) | ||
| 53 | { | ||
| 54 | static_assert(isBox<bulkIdx>(), "This interface is only available for the box method."); | ||
| 55 | 4494 | bulkElementIdx_ = eIdx; | |
| 56 |
1/2✓ Branch 1 taken 4494 times.
✗ Branch 2 not taken.
|
4494 | bulkCornerIndices_ = cornerIndices; |
| 57 |
1/2✓ Branch 1 taken 4494 times.
✗ Branch 2 not taken.
|
4494 | bulkShapeValues_ = shapeValues; |
| 58 | } | ||
| 59 | |||
| 60 | void addLowDimInterpolation(const ShapeValues& shapeValues, | ||
| 61 | const std::vector<GridIndex<lowDimIdx>>& cornerIndices, | ||
| 62 | GridIndex<lowDimIdx> eIdx) | ||
| 63 | { | ||
| 64 | static_assert(isBox<lowDimIdx>(), "This interface is only available for the box method."); | ||
| 65 | 236 | lowDimElementIdx_ = eIdx; | |
| 66 |
1/2✓ Branch 1 taken 236 times.
✗ Branch 2 not taken.
|
236 | lowDimCornerIndices_ = cornerIndices; |
| 67 |
1/2✓ Branch 1 taken 236 times.
✗ Branch 2 not taken.
|
236 | lowDimShapeValues_ = shapeValues; |
| 68 | } | ||
| 69 | |||
| 70 | ✗ | void addBulkInterpolation(GridIndex<bulkIdx> eIdx) | |
| 71 | { | ||
| 72 | static_assert(!isBox<bulkIdx>(), "This interface is not available for the box method."); | ||
| 73 |
1/2✓ Branch 1 taken 81165 times.
✗ Branch 2 not taken.
|
81165 | bulkElementIdx_ = eIdx; |
| 74 | ✗ | } | |
| 75 | |||
| 76 | ✗ | void addLowDimInterpolation(GridIndex<lowDimIdx> eIdx) | |
| 77 | { | ||
| 78 | static_assert(!isBox<lowDimIdx>(), "This interface is not available for the box method."); | ||
| 79 |
2/3✓ Branch 0 taken 48640 times.
✓ Branch 1 taken 36783 times.
✗ Branch 2 not taken.
|
85423 | lowDimElementIdx_ = eIdx; |
| 80 | ✗ | } | |
| 81 | |||
| 82 | 13917199 | PrimaryVariables<bulkIdx> interpolateBulk(const SolutionVector<bulkIdx>& sol) const | |
| 83 | { | ||
| 84 |
4/4✓ Branch 0 taken 24 times.
✓ Branch 1 taken 1018714 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 464198 times.
|
170153426 | PrimaryVariables<bulkIdx> bulkPriVars(0.0); |
| 85 | if (isBox<bulkIdx>()) | ||
| 86 | { | ||
| 87 |
4/4✓ Branch 0 taken 55668796 times.
✓ Branch 1 taken 13917199 times.
✓ Branch 2 taken 55668796 times.
✓ Branch 3 taken 13917199 times.
|
69585995 | for (int i = 0; i < bulkCornerIndices_.size(); ++i) |
| 88 |
2/2✓ Branch 0 taken 55668796 times.
✓ Branch 1 taken 55668796 times.
|
111337592 | for (int priVarIdx = 0; priVarIdx < bulkPriVars.size(); ++priVarIdx) |
| 89 | 334012776 | bulkPriVars[priVarIdx] += sol[bulkCornerIndices_[i]][priVarIdx]*bulkShapeValues_[i]; | |
| 90 | } | ||
| 91 | else | ||
| 92 | { | ||
| 93 |
8/8✓ Branch 0 taken 24 times.
✓ Branch 1 taken 1018714 times.
✓ Branch 2 taken 24 times.
✓ Branch 3 taken 1018714 times.
✓ Branch 4 taken 10 times.
✓ Branch 5 taken 464198 times.
✓ Branch 6 taken 10 times.
✓ Branch 7 taken 464198 times.
|
312472454 | bulkPriVars = sol[bulkElementIdx()]; |
| 94 | } | ||
| 95 | 13917199 | return bulkPriVars; | |
| 96 | } | ||
| 97 | |||
| 98 | 67680 | PrimaryVariables<lowDimIdx> interpolateLowDim(const SolutionVector<lowDimIdx>& sol) const | |
| 99 | { | ||
| 100 |
4/4✓ Branch 0 taken 24 times.
✓ Branch 1 taken 7375757 times.
✓ Branch 2 taken 1227794 times.
✓ Branch 3 taken 464198 times.
|
170107710 | PrimaryVariables<lowDimIdx> lowDimPriVars(0.0); |
| 101 | if (isBox<lowDimIdx>()) | ||
| 102 | { | ||
| 103 |
4/4✓ Branch 0 taken 135360 times.
✓ Branch 1 taken 67680 times.
✓ Branch 2 taken 135360 times.
✓ Branch 3 taken 67680 times.
|
203040 | for (int i = 0; i < lowDimCornerIndices_.size(); ++i) |
| 104 |
2/2✓ Branch 0 taken 135360 times.
✓ Branch 1 taken 135360 times.
|
270720 | for (int priVarIdx = 0; priVarIdx < lowDimPriVars.size(); ++priVarIdx) |
| 105 | 676800 | lowDimPriVars[priVarIdx] += sol[lowDimCornerIndices_[i]][priVarIdx]*lowDimShapeValues_[i]; | |
| 106 | } | ||
| 107 | else | ||
| 108 | { | ||
| 109 |
8/8✓ Branch 0 taken 670572 times.
✓ Branch 1 taken 7375757 times.
✓ Branch 2 taken 24 times.
✓ Branch 3 taken 11457025 times.
✓ Branch 4 taken 10 times.
✓ Branch 5 taken 1691982 times.
✓ Branch 6 taken 10 times.
✓ Branch 7 taken 464198 times.
|
336556028 | lowDimPriVars = sol[lowDimElementIdx()]; |
| 110 | } | ||
| 111 |
2/4✓ Branch 0 taken 670548 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2853484 times.
|
3591712 | return lowDimPriVars; |
| 112 | } | ||
| 113 | |||
| 114 | ✗ | GridIndex<lowDimIdx> lowDimElementIdx() const | |
| 115 | ✗ | { return lowDimElementIdx_; } | |
| 116 | |||
| 117 | ✗ | GridIndex<bulkIdx> bulkElementIdx() const | |
| 118 | ✗ | { return bulkElementIdx_; } | |
| 119 | |||
| 120 | private: | ||
| 121 | ShapeValues bulkShapeValues_, lowDimShapeValues_; | ||
| 122 | std::vector<GridIndex<bulkIdx>> bulkCornerIndices_; | ||
| 123 | std::vector<GridIndex<lowDimIdx>> lowDimCornerIndices_; | ||
| 124 | GridIndex<bulkIdx> bulkElementIdx_; | ||
| 125 | GridIndex<lowDimIdx> lowDimElementIdx_; | ||
| 126 | }; | ||
| 127 | |||
| 128 | /*! | ||
| 129 | * \ingroup EmbeddedCoupling | ||
| 130 | * \brief A point source data class used for integration in multidimensional models | ||
| 131 | * \note The point source and related data are connected via an identifier (id) | ||
| 132 | * When explicitly computing the circle average, i.e. the pressure for | ||
| 133 | * the source term is computed as an integral over the circle describing | ||
| 134 | * the surface of the one-dimensional tube. This exact determination of the bulk pressure | ||
| 135 | * is necessary for convergence studies. | ||
| 136 | */ | ||
| 137 | template<class MDTraits> | ||
| 138 | class PointSourceDataCircleAverage : public PointSourceData<MDTraits> | ||
| 139 | { | ||
| 140 | using ParentType = PointSourceData<MDTraits>; | ||
| 141 | using Scalar = typename MDTraits::Scalar; | ||
| 142 | using ShapeValues = typename std::vector<Dune::FieldVector<Scalar, 1> >; | ||
| 143 | |||
| 144 | template<std::size_t id> using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag; | ||
| 145 | template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>; | ||
| 146 | template<std::size_t id> using GridView = typename GridGeometry<id>::GridView; | ||
| 147 | template<std::size_t id> using GridIndex = typename IndexTraits<GridView<id>>::GridIndex; | ||
| 148 | template<std::size_t id> using SolutionVector = GetPropType<SubDomainTypeTag<id>, Properties::SolutionVector>; | ||
| 149 | template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>; | ||
| 150 | |||
| 151 | static constexpr auto bulkIdx = typename MDTraits::template SubDomain<0>::Index(); | ||
| 152 | static constexpr auto lowDimIdx = typename MDTraits::template SubDomain<1>::Index(); | ||
| 153 | |||
| 154 | template<std::size_t id> | ||
| 155 | static constexpr bool isBox() | ||
| 156 | { return GridGeometry<id>::discMethod == DiscretizationMethods::box; } | ||
| 157 | |||
| 158 | public: | ||
| 159 |
10/16✓ Branch 0 taken 48640 times.
✓ Branch 1 taken 10310 times.
✓ Branch 2 taken 48640 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 58950 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 48640 times.
✓ Branch 7 taken 10310 times.
✓ Branch 8 taken 48640 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 58950 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 10310 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 10310 times.
✗ Branch 17 not taken.
|
353700 | PointSourceDataCircleAverage() : enableBulkCircleInterpolation_(false) {} |
| 160 | |||
| 161 | 14121796 | PrimaryVariables<bulkIdx> interpolateBulk(const SolutionVector<bulkIdx>& sol) const | |
| 162 | { | ||
| 163 | // bulk interpolation on the circle is only enabled for source in the | ||
| 164 | // lower dimensional domain if we use a circle distributed bulk sources | ||
| 165 | // (see coupling manager circle sources) | ||
| 166 | 14121796 | PrimaryVariables<bulkIdx> bulkPriVars(sol[0]); | |
| 167 |
1/2✓ Branch 0 taken 10597764 times.
✗ Branch 1 not taken.
|
14121796 | bulkPriVars = 0.0; |
| 168 |
1/2✓ Branch 0 taken 14121796 times.
✗ Branch 1 not taken.
|
14121796 | if (enableBulkCircleInterpolation_) |
| 169 | { | ||
| 170 | // compute the average of the bulk privars over the circle around the integration point | ||
| 171 | // this computes $\bar{p} = \frac{1}{2\pi R} int_0^2*\pi p R \text{d}\theta. | ||
| 172 | if (isBox<bulkIdx>()) | ||
| 173 | { | ||
| 174 |
3/6✓ Branch 0 taken 4379284 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 4379284 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 4379284 times.
✗ Branch 5 not taken.
|
13137852 | assert(circleCornerIndices_.size() == circleShapeValues_.size()); |
| 175 | |||
| 176 | Scalar weightSum = 0.0; | ||
| 177 |
4/4✓ Branch 0 taken 71220086 times.
✓ Branch 1 taken 4379284 times.
✓ Branch 2 taken 71220086 times.
✓ Branch 3 taken 4379284 times.
|
151198740 | for (std::size_t j = 0; j < circleStencil_.size(); ++j) |
| 178 | { | ||
| 179 | 71220086 | PrimaryVariables<bulkIdx> priVars(0.0); | |
| 180 | 71220086 | const auto& cornerIndices = *(circleCornerIndices_[j]); | |
| 181 | 71220086 | const auto& shapeValues = circleShapeValues_[j]; | |
| 182 |
4/4✓ Branch 0 taken 569760688 times.
✓ Branch 1 taken 71220086 times.
✓ Branch 2 taken 569760688 times.
✓ Branch 3 taken 71220086 times.
|
640980774 | for (int i = 0; i < cornerIndices.size(); ++i) |
| 183 | { | ||
| 184 | 1139521376 | const auto& localSol = sol[cornerIndices[i]]; | |
| 185 | 569760688 | const auto& shapeValue = shapeValues[i]; | |
| 186 |
2/2✓ Branch 0 taken 569760688 times.
✓ Branch 1 taken 569760688 times.
|
1139521376 | for (int priVarIdx = 0; priVarIdx < PrimaryVariables<bulkIdx>::size(); ++priVarIdx) |
| 187 | 1139521376 | priVars[priVarIdx] += localSol[priVarIdx]*shapeValue; | |
| 188 | } | ||
| 189 | // multiply with weight and add | ||
| 190 | 142440172 | priVars *= circleIpWeight_[j]; | |
| 191 | 71220086 | weightSum += circleIpWeight_[j]; | |
| 192 | 142440172 | bulkPriVars += priVars; | |
| 193 | } | ||
| 194 | 4379284 | bulkPriVars /= weightSum; | |
| 195 | } | ||
| 196 | else | ||
| 197 | { | ||
| 198 | Scalar weightSum = 0.0; | ||
| 199 |
4/4✓ Branch 0 taken 902061139 times.
✓ Branch 1 taken 9742512 times.
✓ Branch 2 taken 902061139 times.
✓ Branch 3 taken 9742512 times.
|
1823607302 | for (int j = 0; j < circleStencil_.size(); ++j) |
| 200 | { | ||
| 201 |
2/2✓ Branch 0 taken 937277183 times.
✓ Branch 1 taken 902061139 times.
|
1839338322 | for (int priVarIdx = 0; priVarIdx < bulkPriVars.size(); ++priVarIdx) |
| 202 | 3889972908 | bulkPriVars[priVarIdx] += sol[circleStencil_[j]][priVarIdx]*circleIpWeight_[j]; | |
| 203 | |||
| 204 | 1804122278 | weightSum += circleIpWeight_[j]; | |
| 205 | } | ||
| 206 | 9742512 | bulkPriVars /= weightSum; | |
| 207 | } | ||
| 208 | 14121796 | return bulkPriVars; | |
| 209 | } | ||
| 210 | else | ||
| 211 | { | ||
| 212 | ✗ | return ParentType::interpolateBulk(sol); | |
| 213 | } | ||
| 214 | } | ||
| 215 | |||
| 216 | 2242 | void addCircleInterpolation(const std::vector<const std::vector<GridIndex<bulkIdx>>*>& circleCornerIndices, | |
| 217 | const std::vector<ShapeValues>& circleShapeValues, | ||
| 218 | const std::vector<Scalar>& circleIpWeight, | ||
| 219 | const std::vector<GridIndex<bulkIdx>>& circleStencil) | ||
| 220 | { | ||
| 221 | 2242 | circleCornerIndices_ = circleCornerIndices; | |
| 222 | 2242 | circleShapeValues_ = circleShapeValues; | |
| 223 | 2242 | circleIpWeight_ = circleIpWeight; | |
| 224 | 2242 | circleStencil_ = circleStencil; | |
| 225 | 2242 | enableBulkCircleInterpolation_ = true; | |
| 226 | |||
| 227 | 2242 | } | |
| 228 | |||
| 229 | void addCircleInterpolation(const std::vector<Scalar>& circleIpWeight, | ||
| 230 | const std::vector<GridIndex<bulkIdx>>& circleStencil) | ||
| 231 | { | ||
| 232 |
1/2✓ Branch 1 taken 56708 times.
✗ Branch 2 not taken.
|
56708 | circleIpWeight_ = circleIpWeight; |
| 233 |
1/2✓ Branch 1 taken 56708 times.
✗ Branch 2 not taken.
|
56708 | circleStencil_ = circleStencil; |
| 234 |
1/2✓ Branch 1 taken 8068 times.
✗ Branch 2 not taken.
|
56708 | enableBulkCircleInterpolation_ = true; |
| 235 | } | ||
| 236 | |||
| 237 | const std::vector<GridIndex<bulkIdx>>& circleStencil() const | ||
| 238 | { return circleStencil_; } | ||
| 239 | |||
| 240 | private: | ||
| 241 | std::vector<const std::vector<GridIndex<bulkIdx>>*> circleCornerIndices_; | ||
| 242 | std::vector<ShapeValues> circleShapeValues_; | ||
| 243 | std::vector<Scalar> circleIpWeight_; | ||
| 244 | std::vector<GridIndex<bulkIdx>> circleStencil_; | ||
| 245 | bool enableBulkCircleInterpolation_; | ||
| 246 | }; | ||
| 247 | |||
| 248 | } // end namespace Dumux | ||
| 249 | |||
| 250 | #endif | ||
| 251 |