| 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-FileCopyrightText: 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 Core | ||
| 10 | * \brief A point source class, | ||
| 11 | * i.e. sources located at a single point in space | ||
| 12 | */ | ||
| 13 | |||
| 14 | #ifndef DUMUX_POINTSOURCE_HH | ||
| 15 | #define DUMUX_POINTSOURCE_HH | ||
| 16 | |||
| 17 | #include <functional> | ||
| 18 | |||
| 19 | #include <dune/common/exceptions.hh> | ||
| 20 | #include <dune/common/reservedvector.hh> | ||
| 21 | #include <dumux/common/properties.hh> | ||
| 22 | #include <dumux/common/parameters.hh> | ||
| 23 | #include <dumux/common/numeqvector.hh> | ||
| 24 | #include <dumux/geometry/boundingboxtree.hh> | ||
| 25 | #include <dumux/geometry/intersectspointgeometry.hh> | ||
| 26 | #include <dumux/geometry/intersectingentities.hh> | ||
| 27 | |||
| 28 | #include <dumux/discretization/method.hh> | ||
| 29 | |||
| 30 | namespace Dumux { | ||
| 31 | |||
| 32 | /*! | ||
| 33 | * \ingroup Core | ||
| 34 | * \brief A point source base class | ||
| 35 | * \tparam PositionType the position type | ||
| 36 | * \tparam ValueType the a vector type storing the source for all equations | ||
| 37 | */ | ||
| 38 | template<class PositionType, class ValueType> | ||
| 39 | class PointSource | ||
| 40 | { | ||
| 41 | public: | ||
| 42 | //! Export the scalar type | ||
| 43 | using Scalar = std::decay_t<decltype(std::declval<ValueType>()[0])>; | ||
| 44 | //! Export the position type | ||
| 45 | using GlobalPosition = PositionType; | ||
| 46 | //! Export the value type | ||
| 47 | using Values = ValueType; | ||
| 48 | |||
| 49 | //! Constructor for constant point sources | ||
| 50 | 171177 | PointSource(GlobalPosition pos, Values values) | |
| 51 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
171177 | : values_(values), pos_(pos), embeddings_(1) {} |
| 52 | |||
| 53 | //! Constructor for sol dependent point sources, when there is no | ||
| 54 | // value known at the time of initialization | ||
| 55 | PointSource(GlobalPosition pos) | ||
| 56 | : values_(0.0), pos_(pos), embeddings_(1) {} | ||
| 57 | |||
| 58 | //! Convenience += operator overload modifying only the values | ||
| 59 | PointSource& operator+= (Scalar s) | ||
| 60 | { | ||
| 61 | values_ += s; | ||
| 62 | return *this; | ||
| 63 | } | ||
| 64 | |||
| 65 | //! Convenience -= operator overload modifying only the values | ||
| 66 | PointSource& operator-= (Scalar s) | ||
| 67 | { | ||
| 68 | values_ -= s; | ||
| 69 | return *this; | ||
| 70 | } | ||
| 71 | |||
| 72 | //! Convenience *= operator overload modifying only the values | ||
| 73 | 4080 | PointSource& operator*= (Scalar s) | |
| 74 | { | ||
| 75 | 4080 | values_ *= s; | |
| 76 | return *this; | ||
| 77 | } | ||
| 78 | |||
| 79 | //! Convenience /= operator overload modifying only the values | ||
| 80 | 38547497 | PointSource& operator/= (Scalar s) | |
| 81 | { | ||
| 82 |
1/2✓ Branch 0 taken 72 times.
✗ Branch 1 not taken.
|
42082051 | values_ /= s; |
| 83 | return *this; | ||
| 84 | } | ||
| 85 | |||
| 86 | //! Convenience = operator overload modifying only the values | ||
| 87 | 3524032 | PointSource& operator= (const Values& values) | |
| 88 | { | ||
| 89 | 3524032 | values_ = values; | |
| 90 | return *this; | ||
| 91 | } | ||
| 92 | |||
| 93 | //! Convenience = operator overload modifying only the values | ||
| 94 | 35008223 | PointSource& operator= (Scalar s) | |
| 95 | { | ||
| 96 | 35008223 | values_ = s; | |
| 97 | return *this; | ||
| 98 | } | ||
| 99 | |||
| 100 | //! return the source values | ||
| 101 | 3534554 | Values values() const | |
| 102 | 3534554 | { return values_; } | |
| 103 | |||
| 104 | //! return the source position | ||
| 105 | 239 | const GlobalPosition& position() const | |
| 106 | 11401 | { return pos_; } | |
| 107 | |||
| 108 | //! an update function called before adding the value | ||
| 109 | // to the local residual in the problem in scvPointSources | ||
| 110 | // to be overloaded by derived classes | ||
| 111 | template<class Problem, class FVElementGeometry, class ElementVolumeVariables> | ||
| 112 | void update(const Problem &problem, | ||
| 113 | const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity &element, | ||
| 114 | const FVElementGeometry &fvGeometry, | ||
| 115 | const ElementVolumeVariables &elemVolVars, | ||
| 116 | const typename FVElementGeometry::SubControlVolume &scv) | ||
| 117 | {} | ||
| 118 | |||
| 119 | //! set the number of embeddings for this point source | ||
| 120 | 173684 | void setEmbeddings(std::size_t embeddings) | |
| 121 | { | ||
| 122 |
3/5✓ Branch 1 taken 84736 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 33058 times.
✗ Branch 5 not taken.
✓ Branch 3 taken 48640 times.
|
89052 | embeddings_ = embeddings; |
| 123 | } | ||
| 124 | |||
| 125 | /*! | ||
| 126 | * \brief get the number of embeddings for this point source | ||
| 127 | * \note A point source might be located on the intersection between several scvs. | ||
| 128 | * If so, there are point sources for every neighboring scv with the same position. | ||
| 129 | * `embeddings` returns the number of neighboring scvs. | ||
| 130 | * Example: If I want to inject 1kg/s at a location that is on the inner face of an scv | ||
| 131 | * the point source exists in both scvs. Both have a value of 1kg/s. | ||
| 132 | * We then divide the value by the number of embeddings to not inject 2kg/s but 1kg/s. | ||
| 133 | * \note This division is done in the problem.scvPointSources() if this behaviour is not explicitly | ||
| 134 | * changed by e.g. overloading this function in the problem implementation. | ||
| 135 | */ | ||
| 136 | 37921398 | std::size_t embeddings() const | |
| 137 | { | ||
| 138 |
1/2✓ Branch 0 taken 72 times.
✗ Branch 1 not taken.
|
38554060 | return embeddings_; |
| 139 | } | ||
| 140 | |||
| 141 | protected: | ||
| 142 | Values values_; //!< value of the point source for each equation | ||
| 143 | private: | ||
| 144 | GlobalPosition pos_; //!< position of the point source | ||
| 145 | std::size_t embeddings_; //!< how many SCVs the point source is associated with | ||
| 146 | }; | ||
| 147 | |||
| 148 | /*! | ||
| 149 | * \ingroup Core | ||
| 150 | * \brief A point source class with an identifier to attach data | ||
| 151 | * \tparam GlobalPosition the position type | ||
| 152 | * \tparam SourceValues the a vector type storing the source for all equations | ||
| 153 | * \tparam I the ID type | ||
| 154 | */ | ||
| 155 | template<class GlobalPosition, class SourceValues, class I> | ||
| 156 | class IdPointSource : public PointSource<GlobalPosition, SourceValues> | ||
| 157 | { | ||
| 158 | using ParentType = PointSource<GlobalPosition, SourceValues>; | ||
| 159 | using Scalar = typename ParentType::Scalar; | ||
| 160 | |||
| 161 | public: | ||
| 162 | //! export the id type | ||
| 163 | using IdType = I; | ||
| 164 | |||
| 165 | //! Constructor for constant point sources | ||
| 166 | IdPointSource(GlobalPosition pos, SourceValues values, IdType id) | ||
| 167 | : ParentType(pos, values), id_(id) {} | ||
| 168 | |||
| 169 | //! Constructor for sol dependent point sources, when there is no | ||
| 170 | // value known at the time of initialization | ||
| 171 | 170938 | IdPointSource(GlobalPosition pos, IdType id) | |
| 172 | 175308 | : ParentType(pos, SourceValues(0.0)), id_(id) {} | |
| 173 | |||
| 174 | //! return the sources identifier | ||
| 175 | 38532255 | IdType id() const | |
| 176 |
4/4✓ Branch 0 taken 670560 times.
✓ Branch 1 taken 6541223 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 4377727 times.
|
25749805 | { return id_; } |
| 177 | |||
| 178 | //! Convenience = operator overload modifying only the values | ||
| 179 | 3524032 | IdPointSource& operator= (const SourceValues& values) | |
| 180 | { | ||
| 181 | 3524032 | ParentType::operator=(values); | |
| 182 | return *this; | ||
| 183 | } | ||
| 184 | |||
| 185 | //! Convenience = operator overload modifying only the values | ||
| 186 | 35008223 | IdPointSource& operator= (Scalar s) | |
| 187 | { | ||
| 188 | 35008223 | ParentType::operator=(s); | |
| 189 | return *this; | ||
| 190 | } | ||
| 191 | |||
| 192 | private: | ||
| 193 | IdType id_; | ||
| 194 | }; | ||
| 195 | |||
| 196 | /*! | ||
| 197 | * \ingroup Core | ||
| 198 | * \brief A point source class for time dependent point sources | ||
| 199 | */ | ||
| 200 | template<class TypeTag> | ||
| 201 |
15/43✓ Branch 1 taken 72 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 72 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 8 times.
✓ Branch 12 taken 4 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 17 taken 4 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✓ Branch 21 taken 19 times.
✗ Branch 22 not taken.
✗ Branch 24 not taken.
✓ Branch 25 taken 9368 times.
✗ Branch 26 not taken.
✓ Branch 27 taken 9368 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✓ Branch 30 taken 9 times.
✗ Branch 31 not taken.
✓ Branch 34 taken 3 times.
✗ Branch 35 not taken.
✓ Branch 36 taken 1 times.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✓ Branch 39 taken 1 times.
✓ Branch 40 taken 1 times.
✗ Branch 41 not taken.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 0 not taken.
✗ Branch 11 not taken.
✗ Branch 16 not taken.
✓ Branch 32 taken 2 times.
✗ Branch 33 not taken.
|
28407 | class SolDependentPointSource : public PointSource<Dune::FieldVector<typename GetPropType<TypeTag, Properties::GridGeometry>::GridView::ctype, |
| 202 | GetPropType<TypeTag, Properties::GridGeometry>::GridView::dimensionworld>, | ||
| 203 | Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>> | ||
| 204 | { | ||
| 205 | using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView; | ||
| 206 | using SourceValues = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>; | ||
| 207 | using Problem = GetPropType<TypeTag, Properties::Problem>; | ||
| 208 | using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView; | ||
| 209 | using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView; | ||
| 210 | using SubControlVolume = typename FVElementGeometry::SubControlVolume; | ||
| 211 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 212 | |||
| 213 | static const int dimworld = GridView::dimensionworld; | ||
| 214 | using GlobalPosition = typename Dune::FieldVector<typename GridView::ctype, dimworld>; | ||
| 215 | // returns the PointSource values as PrimaryVariables | ||
| 216 | using ValueFunction = typename std::function<SourceValues(const Problem &problem, | ||
| 217 | const Element &element, | ||
| 218 | const FVElementGeometry &fvGeometry, | ||
| 219 | const ElementVolumeVariables &elemVolVars, | ||
| 220 | const SubControlVolume &scv)>; | ||
| 221 | |||
| 222 | using ParentType = PointSource<GlobalPosition, SourceValues>; | ||
| 223 | using Scalar = typename ParentType::Scalar; | ||
| 224 | |||
| 225 | public: | ||
| 226 | //! Constructor for sol dependent point sources, when there is no | ||
| 227 | // value known at the time of initialization | ||
| 228 | 3 | SolDependentPointSource(GlobalPosition pos, | |
| 229 | ValueFunction valueFunction) | ||
| 230 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
5 | : ParentType(pos, SourceValues(0.0)), valueFunction_(valueFunction) {} |
| 231 | |||
| 232 | //! an update function called before adding the value | ||
| 233 | // to the local residual in the problem in scvPointSources | ||
| 234 | // to be overloaded by derived classes | ||
| 235 | 9440 | void update(const Problem &problem, | |
| 236 | const Element &element, | ||
| 237 | const FVElementGeometry &fvGeometry, | ||
| 238 | const ElementVolumeVariables &elemVolVars, | ||
| 239 | const SubControlVolume &scv) | ||
| 240 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 9440 times.
|
9440 | { this->values_ = valueFunction_(problem, element, fvGeometry, elemVolVars, scv); } |
| 241 | |||
| 242 | //! Convenience = operator overload modifying only the values | ||
| 243 | SolDependentPointSource& operator= (const SourceValues& values) | ||
| 244 | { | ||
| 245 | ParentType::operator=(values); | ||
| 246 | return *this; | ||
| 247 | } | ||
| 248 | |||
| 249 | //! Convenience = operator overload modifying only the values | ||
| 250 | SolDependentPointSource& operator= (Scalar s) | ||
| 251 | { | ||
| 252 | ParentType::operator=(s); | ||
| 253 | return *this; | ||
| 254 | } | ||
| 255 | |||
| 256 | private: | ||
| 257 | ValueFunction valueFunction_; | ||
| 258 | }; | ||
| 259 | |||
| 260 | /*! | ||
| 261 | * \ingroup Core | ||
| 262 | * \brief A helper class calculating a sub control volume to point source map | ||
| 263 | * This class uses the bounding box tree implementation to identify in which | ||
| 264 | * sub control volume(s) a point source falls. | ||
| 265 | */ | ||
| 266 | class BoundingBoxTreePointSourceHelper | ||
| 267 | { | ||
| 268 | public: | ||
| 269 | //! calculate a DOF index to point source map from given vector of point sources | ||
| 270 | template<class GridGeometry, class PointSource, class PointSourceMap> | ||
| 271 | 26 | static void computePointSourceMap(const GridGeometry& gridGeometry, | |
| 272 | const std::vector<PointSource>& sources, | ||
| 273 | PointSourceMap& pointSourceMap, | ||
| 274 | const std::string& paramGroup = "") | ||
| 275 | { | ||
| 276 | 26 | const auto& boundingBoxTree = gridGeometry.boundingBoxTree(); | |
| 277 | |||
| 278 |
4/4✓ Branch 0 taken 7 times.
✓ Branch 1 taken 236 times.
✓ Branch 2 taken 22 times.
✓ Branch 3 taken 3 times.
|
480 | for (const auto& s : sources) |
| 279 | { | ||
| 280 | // compute in which elements the point source falls | ||
| 281 | 239 | const auto entities = intersectingEntities(s.position(), boundingBoxTree); | |
| 282 | |||
| 283 | // continue with next point source if no intersection with the grid are found | ||
| 284 |
2/4✓ Branch 0 taken 20 times.
✓ Branch 1 taken 219 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
239 | if (entities.empty()) |
| 285 | continue; | ||
| 286 | |||
| 287 | // make local copy of point source for the map | ||
| 288 | 219 | auto source = s; | |
| 289 | |||
| 290 | // split the source values equally among all concerned entities | ||
| 291 | 219 | source.setEmbeddings(entities.size()*source.embeddings()); | |
| 292 | |||
| 293 | if constexpr (GridGeometry::discMethod == DiscretizationMethods::box | ||
| 294 | || GridGeometry::discMethod == DiscretizationMethods::fcdiamond) | ||
| 295 | { | ||
| 296 | // loop over all concerned elements | ||
| 297 | 23 | auto fvGeometry = localView(gridGeometry); | |
| 298 |
4/4✓ Branch 1 taken 52 times.
✓ Branch 2 taken 20 times.
✓ Branch 3 taken 20 times.
✓ Branch 4 taken 3 times.
|
75 | for (const auto eIdx : entities) |
| 299 | { | ||
| 300 | // check in which subcontrolvolume(s) we are | ||
| 301 |
1/2✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
|
52 | const auto element = boundingBoxTree.entitySet().entity(eIdx); |
| 302 |
1/2✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
|
52 | fvGeometry.bindElement(element); |
| 303 | |||
| 304 | 52 | const auto globalPos = source.position(); | |
| 305 | // loop over all sub control volumes and check if the point source is inside | ||
| 306 | 52 | constexpr int dim = GridGeometry::GridView::dimension; | |
| 307 |
2/2✓ Branch 0 taken 120 times.
✓ Branch 1 taken 40 times.
|
172 | Dune::ReservedVector<std::size_t, 1<<dim> scvIndices; |
| 308 |
2/2✓ Branch 0 taken 232 times.
✓ Branch 1 taken 52 times.
|
284 | for (const auto& scv : scvs(fvGeometry)) |
| 309 |
3/4✓ Branch 1 taken 232 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 52 times.
✓ Branch 5 taken 180 times.
|
232 | if (intersectsPointGeometry(globalPos, fvGeometry.geometry(scv))) |
| 310 | 52 | scvIndices.push_back(scv.indexInElement()); | |
| 311 | |||
| 312 | // for all scvs that tested positive add the point sources | ||
| 313 | // to the element/scv to point source map | ||
| 314 |
3/4✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 52 times.
✓ Branch 3 taken 52 times.
|
104 | for (const auto scvIdx : scvIndices) |
| 315 | { | ||
| 316 |
1/2✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
|
52 | const auto key = std::make_pair(eIdx, scvIdx); |
| 317 | 52 | if (pointSourceMap.count(key)) | |
| 318 | ✗ | pointSourceMap.at(key).push_back(source); | |
| 319 | else | ||
| 320 |
5/11✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 4 times.
✓ Branch 6 taken 4 times.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✓ Branch 4 taken 48 times.
|
104 | pointSourceMap.insert({key, {source}}); |
| 321 | // split equally on the number of matched scvs | ||
| 322 |
1/2✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
|
52 | auto& s = pointSourceMap.at(key).back(); |
| 323 | 52 | s.setEmbeddings(scvIndices.size()*s.embeddings()); | |
| 324 | } | ||
| 325 | } | ||
| 326 | 23 | } | |
| 327 | else if constexpr (GridGeometry::discMethod == DiscretizationMethods::cctpfa | ||
| 328 | || GridGeometry::discMethod == DiscretizationMethods::ccmpfa) | ||
| 329 | { | ||
| 330 |
2/2✓ Branch 0 taken 269 times.
✓ Branch 1 taken 196 times.
|
465 | for (const auto eIdx : entities) |
| 331 | { | ||
| 332 | // add the pointsource to the DOF map | ||
| 333 |
1/2✓ Branch 0 taken 269 times.
✗ Branch 1 not taken.
|
269 | const auto key = std::make_pair(eIdx, /*scvIdx=*/ 0); |
| 334 |
2/6✓ Branch 0 taken 269 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 269 times.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
269 | if (pointSourceMap.count(key)) |
| 335 | ✗ | pointSourceMap.at(key).push_back(source); | |
| 336 | else | ||
| 337 |
5/14✓ Branch 1 taken 269 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 8 times.
✓ Branch 6 taken 8 times.
✓ Branch 7 taken 8 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 4 taken 261 times.
✗ Branch 9 not taken.
|
538 | pointSourceMap.insert({key, {source}}); |
| 338 | } | ||
| 339 | } | ||
| 340 | else | ||
| 341 | DUNE_THROW(Dune::NotImplemented, "BoundingBoxTreePointSourceHelper for discretization method " << GridGeometry::discMethod); | ||
| 342 | } | ||
| 343 | 38 | } | |
| 344 | }; | ||
| 345 | |||
| 346 | } // end namespace Dumux | ||
| 347 | |||
| 348 | #endif | ||
| 349 |