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