diff --git a/dumux/common/fvproblem.hh b/dumux/common/fvproblem.hh index a9857b41d0d03bea6798584faa01cd715af67d9e..dd8b8956a86594b09e94b5d4655f05172ac4836b 100644 --- a/dumux/common/fvproblem.hh +++ b/dumux/common/fvproblem.hh @@ -25,6 +25,7 @@ #define DUMUX_FV_PROBLEM_HH #include +#include #include #include @@ -75,6 +76,9 @@ class FVProblem static constexpr bool isBox = GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::Box; + using PointSourceMap = std::map, + std::vector >; + public: /*! * \brief Constructor @@ -418,7 +422,7 @@ public: // call convienience problem interface function asImp_().pointSource(pointSource, element, fvGeometry, elemVolVars, scv); // at last take care about multiplying with the correct volume - pointSource /= volume; + pointSource /= volume*pointSource.embeddings(); // add the point source values to the local residual source += pointSource.values(); } @@ -431,7 +435,7 @@ public: * \brief Compute the point source map, i.e. which scvs have point source contributions * \note Call this on the problem before assembly if you want to enable point sources set * via the addPointSources member function. - */ + */ void computePointSourceMap() { // clear the given point source maps in case it's not empty @@ -451,10 +455,16 @@ public: } } + /*! + * \brief Get the point source map. It stores the point sources per scv + */ + const PointSourceMap& pointSourceMap() const + { return pointSourceMap_; } + /*! * \brief Applies the initial solution for all degrees of freedom of the grid. * \param sol the initial solution vector - */ + */ void applyInitialSolution(SolutionVector& sol) const { // set the initial values by forwarding to a specialized method @@ -517,8 +527,8 @@ public: */ Scalar extrusionFactorAtPos(const GlobalPosition &globalPos) const { - //! As a default, i.e. if the user's problem does not overload any extrusion factor method - //! return 1.0 + // As a default, i.e. if the user's problem does not overload + // any extrusion factor method, return 1.0 return 1.0; } @@ -569,8 +579,7 @@ private: std::string problemName_; //! A map from an scv to a vector of point sources - std::map, - std::vector > pointSourceMap_; + PointSourceMap pointSourceMap_; }; } // end namespace Dumux diff --git a/dumux/common/pointsource.hh b/dumux/common/pointsource.hh index ac4efb84f9b836800c49ba0fd715fffa4d894384..271228fe8e1c6a34e92ca007e659d641025e9216 100644 --- a/dumux/common/pointsource.hh +++ b/dumux/common/pointsource.hh @@ -115,14 +115,8 @@ public: } //! return the source values - // don't forget to call this when it's overloaded - // in the derived class PrimaryVariables values() const - { - auto values = PrimaryVariables(values_); - values /= embeddings_; - return values; - } + { return values_; } //! return the source position const GlobalPosition& position() const @@ -144,7 +138,17 @@ public: embeddings_ = embeddings; } - //! get the number of embeddings for this point source + /*! + * \brief get the number of embeddings for this point source + * \note A point source might be located on the intersection between several scvs. + * If so, there are point sources for every neighboring scv with the same position. + * `embeddings` returns the number of neighboring scvs. + * Example: If I want to inject 1kg/s at a location that is on the inner face of an scv + * the point source exists in both scvs. Both have a value of 1kg/s. + * We then divide the value by the number of embeddings to not inject 2kg/s but 1kg/s. + * \note This division is done in the problem.scvPointSources() if this behaviour is not explicitly + * changed by e.g. overloading this function in the problem implementation. + */ std::size_t embeddings() const { return embeddings_; @@ -286,9 +290,10 @@ class BoundingBoxTreePointSourceHelper public: //! calculate a DOF index to point source map from given vector of point sources + template static void computePointSourceMap(const FVGridGeometry& fvGridGeometry, std::vector& sources, - std::map, std::vector >& pointSourceMap) + PointSourceMap& pointSourceMap) { const auto& boundingBoxTree = fvGridGeometry.boundingBoxTree();