diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh index 53e647e594d54014599d6d30fa6a6e1c62b7ed31..316047d549294f8c0c896ad3a1a81dcfa851b90d 100644 --- a/dumux/multidomain/embedded/couplingmanager1d3d.hh +++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh @@ -26,12 +26,16 @@ #ifndef DUMUX_MULTIDOMAIN_EMBEDDED_COUPLINGMANAGER_1D3D_HH #define DUMUX_MULTIDOMAIN_EMBEDDED_COUPLINGMANAGER_1D3D_HH +#include #include +#include #include +#include #include #include +#include #include #include #include @@ -39,6 +43,7 @@ #include #include #include +#include namespace Dumux { @@ -47,7 +52,26 @@ namespace Dumux { * \brief The coupling mode */ enum class EmbeddedCouplingMode -{ line, average, cylindersources, kernel }; +{ line, average, cylindersources, kernel, kernelAnisotropic, wellIndex, projection }; + +/*! + * \ingroup EmbeddedCoupling + * \brief Convert the coupling mode to a human readable name + */ +std::string toString(EmbeddedCouplingMode mode) +{ + switch (mode) + { + case EmbeddedCouplingMode::line: return "line"; + case EmbeddedCouplingMode::average: return "average"; + case EmbeddedCouplingMode::cylindersources: return "cylinder"; + case EmbeddedCouplingMode::kernel: return "kernel"; + case EmbeddedCouplingMode::kernelAnisotropic: return "kernelAnisotropic"; + case EmbeddedCouplingMode::wellIndex: return "wellIndex"; + case EmbeddedCouplingMode::projection: return "projection"; + default: return "unknown"; + } +} //! point source traits for the circle average coupling mode template @@ -296,15 +320,20 @@ public: this->precomputeVertexIndices(bulkIdx); this->precomputeVertexIndices(lowDimIdx); - // iterate over all lowdim elements - const auto& lowDimProblem = this->problem(lowDimIdx); - for (const auto& lowDimElement : elements(this->gridView(lowDimIdx))) + // intersect the bounding box trees + this->glueGrids(); + + // iterate over all intersection and add point sources + for (const auto& is : intersections(this->glue())) { - // get the Gaussian quadrature rule for the low dim element - const auto lowDimGeometry = lowDimElement.geometry(); - const auto& quad = Dune::QuadratureRules::rule(lowDimGeometry.type(), order); + // all inside elements are identical... + const auto& lowDimElement = is.targetEntity(0); + const auto lowDimElementIdx = lowDimFvGridGeometry.elementMapper().index(lowDimElement); - const auto lowDimElementIdx = lowDimGridGeometry.elementMapper().index(lowDimElement); + // get the intersection geometry + const auto intersectionGeometry = is.geometry(); + // get the Gaussian quadrature rule for the local intersection + const auto& quad = Dune::QuadratureRules::rule(intersectionGeometry.type(), order); // apply the Gaussian quadrature rule and define point sources at each quadrature point // note that the approximation is not optimal if @@ -316,7 +345,7 @@ public: for (auto&& qp : quad) { // global position of the quadrature point - const auto globalPos = lowDimGeometry.global(qp.position()); + const auto globalPos = intersectionGeometry.global(qp.position()); const auto bulkElementIndices = intersectingEntities(globalPos, bulkTree); @@ -330,14 +359,15 @@ public: ////////////////////////////////////////////////////////// static const auto numIp = getParam("MixedDimension.NumCircleSegments"); - const auto radius = lowDimProblem.spatialParams().radius(lowDimElementIdx); - const auto normal = lowDimGeometry.corner(1)-lowDimGeometry.corner(0); + const auto radius = this->problem(lowDimIdx).spatialParams().radius(lowDimElementIdx); + const auto normal = intersectionGeometry.corner(1)-intersectionGeometry.corner(0); const auto circleAvgWeight = 2*M_PI*radius/numIp; + const auto integrationElement = intersectionGeometry.integrationElement(qp.position()); + const auto qpweight = qp.weight(); const auto circlePoints = EmbeddedCoupling::circlePoints(globalPos, normal, radius, numIp); std::vector circleIpWeight; circleIpWeight.reserve(circlePoints.size()); std::vector> circleStencil; circleStencil.reserve(circlePoints.size()); - // for box std::vector>*> circleCornerIndices; using ShapeValues = std::vector >; @@ -358,16 +388,16 @@ public: circleStencil.push_back(bulkElementIdx); circleIpWeight.push_back(localCircleAvgWeight); - // precompute interpolation data for box scheme for each cut element - if (isBox()) + // precompute interpolation data for box scheme for each cut bulk element + if constexpr (isBox()) { - const auto bulkElement = bulkGridGeometry.element(bulkElementIdx); + const auto bulkElement = bulkFvGridGeometry.element(bulkElementIdx); circleCornerIndices.push_back(&(this->vertexIndices(bulkIdx, bulkElementIdx))); // evaluate shape functions at the integration point const auto bulkGeometry = bulkElement.geometry(); ShapeValues shapeValues; - this->getShapeValues(bulkIdx, bulkGridGeometry, bulkGeometry, circlePoints[k], shapeValues); + this->getShapeValues(bulkIdx, bulkFvGridGeometry, bulkGeometry, circlePoints[k], shapeValues); circleShapeValues.emplace_back(std::move(shapeValues)); } } @@ -390,16 +420,16 @@ public: circleStencil.begin(), circleStencil.end()); } + // surface fraction that is inside the domain (only different than 1.0 on the boundary) + const auto surfaceFraction = Scalar(insideCirclePoints)/Scalar(circlePoints.size()); // loop over the bulk elements at the integration points (usually one except when it is on a face or edge or vertex) for (auto bulkElementIdx : bulkElementIndices) { const auto id = this->idCounter_++; - const auto ie = lowDimGeometry.integrationElement(qp.position()); - const auto qpweight = qp.weight(); - this->pointSources(bulkIdx).emplace_back(globalPos, id, qpweight, ie, std::vector({bulkElementIdx})); + this->pointSources(bulkIdx).emplace_back(globalPos, id, qpweight, integrationElement*surfaceFraction, bulkElementIdx); this->pointSources(bulkIdx).back().setEmbeddings(bulkElementIndices.size()); - this->pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, ie, std::vector({lowDimElementIdx})); + this->pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, integrationElement*surfaceFraction, lowDimElementIdx); this->pointSources(lowDimIdx).back().setEmbeddings(bulkElementIndices.size()); // pre compute additional data used for the evaluation of @@ -409,7 +439,7 @@ public: if constexpr (isBox()) { ShapeValues shapeValues; - this->getShapeValues(lowDimIdx, lowDimGridGeometry, lowDimGeometry, globalPos, shapeValues); + this->getShapeValues(lowDimIdx, lowDimFvGridGeometry, intersectionGeometry, globalPos, shapeValues); psData.addLowDimInterpolation(shapeValues, this->vertexIndices(lowDimIdx, lowDimElementIdx), lowDimElementIdx); } else @@ -422,9 +452,9 @@ public: { psData.addCircleInterpolation(circleCornerIndices, circleShapeValues, circleIpWeight, circleStencil); - const auto bulkGeometry = bulkGridGeometry.element(bulkElementIdx).geometry(); + const auto bulkGeometry = bulkFvGridGeometry.element(bulkElementIdx).geometry(); ShapeValues shapeValues; - this->getShapeValues(bulkIdx, bulkGridGeometry, bulkGeometry, globalPos, shapeValues); + this->getShapeValues(bulkIdx, bulkFvGridGeometry, bulkGeometry, globalPos, shapeValues); psData.addBulkInterpolation(shapeValues, this->vertexIndices(bulkIdx, bulkElementIdx), bulkElementIdx); } else @@ -436,6 +466,21 @@ public: // publish point source data in the global vector this->pointSourceData().emplace_back(std::move(psData)); + // compute the average flux scaling factor + const auto outsideGeometry = bulkFvGridGeometry.element(bulkElementIdx).geometry(); + const auto& quadBulk = Dune::QuadratureRules::rule(outsideGeometry.type(), 5); + // iterate over all quadrature points + Scalar avgMinDist = 0.0; + for (auto&& qpBulk : quadBulk) + { + const auto globalPosBulk = outsideGeometry.global(qpBulk.position()); + const auto d = computeMinDistance_(intersectionGeometry.corner(0), intersectionGeometry.corner(1), globalPosBulk); + avgMinDist += d*qpBulk.weight(); + } + + // add it to the internal data vector + this->averageDistanceToBulkCell().push_back(avgMinDist); + // export the bulk coupling stencil if (isBox()) { @@ -565,6 +610,16 @@ public: // \} private: + inline Scalar computeMinDistance_(const GlobalPosition& a, const GlobalPosition& b, const GlobalPosition& x3d) const + { + // r is the minimum distance to the line through a and b + const auto ab = b - a; + const auto t = (x3d - a)*ab/ab.two_norm2(); + auto proj = a; proj.axpy(t, ab); + const auto r = (proj - x3d).two_norm(); + return r; + } + //! the extended source stencil object EmbeddedCoupling::ExtendedSourceStencil extendedSourceStencil_; @@ -582,10 +637,11 @@ private: */ template class EmbeddedCouplingManager1d3d -: public EmbeddedCouplingManagerBase> +: public EmbeddedCouplingManagerBase, + CircleAveragePointSourceTraits> { using ThisType = EmbeddedCouplingManager1d3d; - using ParentType = EmbeddedCouplingManagerBase; + using ParentType = EmbeddedCouplingManagerBase>; using Scalar = typename MDTraits::Scalar; using SolutionVector = typename MDTraits::SolutionVector; using PointSourceData = typename ParentType::PointSourceTraits::PointSourceData; @@ -625,6 +681,35 @@ public: computeLowDimVolumeFractions(); } + /*! + * \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 + void extendJacobianPattern(Dune::index_constant domainI, JacobianPattern& pattern) const + { + extendedSourceStencil_.extendJacobianPattern(*this, domainI, pattern); + } + + /*! + * \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 + void evalAdditionalDomainDerivatives(Dune::index_constant domainI, + const LocalAssemblerI& localAssemblerI, + const typename LocalAssemblerI::LocalResidual::ElementResidualVector&, + JacobianMatrixDiagBlock& A, + GridVariables& gridVariables) + { + extendedSourceStencil_.evalAdditionalDomainDerivatives(*this, domainI, localAssemblerI, this->curSol(), A, gridVariables); + } + /* \brief Compute integration point point sources and associated data * * This method uses grid glue to intersect the given grids. Over each intersection @@ -636,6 +721,9 @@ public: */ void computePointSourceData(std::size_t order = 1, bool verbose = false) { + // if we use the circle average as the 3D values or a point evaluation + static const bool useCircleAverage = getParam("MixedDimension.UseCircleAverage", true); + // Initialize the bulk bounding box tree const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry(); const auto& lowDimGridGeometry = this->problem(lowDimIdx).gridGeometry(); @@ -649,20 +737,28 @@ public: // clear all internal members like pointsource vectors and stencils // initializes the point source id counter this->clear(); + extendedSourceStencil_.stencil().clear(); // precompute the vertex indices for efficiency this->precomputeVertexIndices(bulkIdx); this->precomputeVertexIndices(lowDimIdx); + // intersect the bounding box trees + this->glueGrids(); + // iterate over all lowdim elements const auto& lowDimProblem = this->problem(lowDimIdx); - for (const auto& lowDimElement : elements(this->gridView(lowDimIdx))) + // iterate over all intersection and add point sources + for (const auto& is : intersections(this->glue())) { - // get the Gaussian quadrature rule for the low dim element - const auto lowDimGeometry = lowDimElement.geometry(); - const auto& quad = Dune::QuadratureRules::rule(lowDimGeometry.type(), order); + // all inside elements are identical... + const auto& lowDimElement = is.inside(0); + const auto lowDimElementIdx = lowDimFvGridGeometry.elementMapper().index(lowDimElement); - const auto lowDimElementIdx = lowDimGridGeometry.elementMapper().index(lowDimElement); + // get the intersection geometry + const auto intersectionGeometry = is.geometry(); + // get the Gaussian quadrature rule for the local intersection + const auto& quad = Dune::QuadratureRules::rule(intersectionGeometry.type(), order); // apply the Gaussian quadrature rule and define point sources at each quadrature point // note that the approximation is not optimal if @@ -674,7 +770,7 @@ public: for (auto&& qp : quad) { // global position of the quadrature point - const auto globalPos = lowDimGeometry.global(qp.position()); + const auto globalPos = intersectionGeometry.global(qp.position()); const auto bulkElementIndices = intersectingEntities(globalPos, bulkTree); @@ -687,29 +783,85 @@ public: // get points on the cylinder surface at the integration point //////////////////////////////////////////////////////////////// - static const auto numIp = getParam("MixedDimension.NumCircleSegments", 25); + static const auto numIp = getParam("MixedDimension.NumCircleSegments"); const auto radius = lowDimProblem.spatialParams().radius(lowDimElementIdx); - const auto normal = lowDimGeometry.corner(1)-lowDimGeometry.corner(0); - const auto integrationElement = lowDimGeometry.integrationElement(qp.position())*2*M_PI*radius/Scalar(numIp); - const auto weight = qp.weight()/(2*M_PI*radius); + const auto normal = intersectionGeometry.corner(1)-intersectionGeometry.corner(0); + const auto integrationElement = intersectionGeometry.integrationElement(qp.position())*2*M_PI*radius/Scalar(numIp); + const auto qpweight = qp.weight()/(2*M_PI*radius); + const auto circleAvgWeight = 2*M_PI*radius/numIp; const auto circlePoints = EmbeddedCoupling::circlePoints(globalPos, normal, radius, numIp); + std::vector> circleBulkElementIndices(circlePoints.size()); + std::vector circleIpWeight; circleIpWeight.reserve(circlePoints.size()); + std::vector> circleStencil; circleStencil.reserve(circlePoints.size()); + // for box + std::vector>*> circleCornerIndices; + using ShapeValues = std::vector >; + std::vector circleShapeValues; + + // go over all points of the average operator + for (int k = 0; k < circlePoints.size(); ++k) + { + circleBulkElementIndices[k] = intersectingEntities(circlePoints[k], bulkTree); + if (circleBulkElementIndices[k].empty()) + continue; + + const auto localCircleAvgWeight = circleAvgWeight / circleBulkElementIndices[k].size(); + for (const auto bulkElementIdx : circleBulkElementIndices[k]) + { + circleStencil.push_back(bulkElementIdx); + circleIpWeight.push_back(localCircleAvgWeight); + + // precompute interpolation data for box scheme for each cut bulk element + if (isBox()) + { + const auto bulkElement = bulkFvGridGeometry.element(bulkElementIdx); + circleCornerIndices.push_back(&(this->vertexIndices(bulkIdx, bulkElementIdx))); + + // evaluate shape functions at the integration point + const auto bulkGeometry = bulkElement.geometry(); + ShapeValues shapeValues; + this->getShapeValues(bulkIdx, bulkFvGridGeometry, bulkGeometry, circlePoints[k], shapeValues); + circleShapeValues.emplace_back(std::move(shapeValues)); + } + } + } + + // export low dim circle stencil + if (isBox()) + { + // we insert all vertices and make it unique later + for (const auto& vertices : circleCornerIndices) + { + this->couplingStencils(lowDimIdx)[lowDimElementIdx].insert(this->couplingStencils(lowDimIdx)[lowDimElementIdx].end(), + vertices->begin(), vertices->end()); + + } + } + else + { + this->couplingStencils(lowDimIdx)[lowDimElementIdx].insert(this->couplingStencils(lowDimIdx)[lowDimElementIdx].end(), + circleStencil.begin(), circleStencil.end()); + } + for (int k = 0; k < circlePoints.size(); ++k) { const auto& circlePos = circlePoints[k]; - const auto circleBulkElementIndices = intersectingEntities(circlePos, bulkTree); - if (circleBulkElementIndices.empty()) + + // find bulk elements intersection with the circle elements + if (circleBulkElementIndices[k].empty()) continue; // loop over the bulk elements at the integration points (usually one except when it is on a face or edge or vertex) // and add a point source at every point on the circle - for (const auto bulkElementIdx : circleBulkElementIndices) + for (const auto bulkElementIdx : circleBulkElementIndices[k]) { const auto id = this->idCounter_++; - this->pointSources(bulkIdx).emplace_back(circlePos, id, weight, integrationElement, std::vector({bulkElementIdx})); - this->pointSources(bulkIdx).back().setEmbeddings(circleBulkElementIndices.size()); - this->pointSources(lowDimIdx).emplace_back(globalPos, id, weight, integrationElement, std::vector({lowDimElementIdx})); - this->pointSources(lowDimIdx).back().setEmbeddings(circleBulkElementIndices.size()); + + this->pointSources(bulkIdx).emplace_back(circlePos, id, qpweight, integrationElement, bulkElementIdx); + this->pointSources(bulkIdx).back().setEmbeddings(circleBulkElementIndices[k].size()); + this->pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, integrationElement, lowDimElementIdx); + this->pointSources(lowDimIdx).back().setEmbeddings(circleBulkElementIndices[k].size()); // pre compute additional data used for the evaluation of // the actual solution dependent source term @@ -719,7 +871,7 @@ public: { using ShapeValues = std::vector >; ShapeValues shapeValues; - this->getShapeValues(lowDimIdx, lowDimGridGeometry, lowDimGeometry, globalPos, shapeValues); + this->getShapeValues(lowDimIdx, lowDimGridGeometry, intersectionGeometry, globalPos, shapeValues); psData.addLowDimInterpolation(shapeValues, this->vertexIndices(lowDimIdx, lowDimElementIdx), lowDimElementIdx); } else @@ -730,6 +882,9 @@ public: // add data needed to compute integral over the circle if constexpr (isBox()) { + if (useCircleAverage) + psData.addCircleInterpolation(circleCornerIndices, circleShapeValues, circleIpWeight, circleStencil); + using ShapeValues = std::vector >; const auto bulkGeometry = bulkGridGeometry.element(bulkElementIdx).geometry(); ShapeValues shapeValues; @@ -738,25 +893,14 @@ public: } else { + if (useCircleAverage) + psData.addCircleInterpolation(circleIpWeight, circleStencil); psData.addBulkInterpolation(bulkElementIdx); } // publish point source data in the global vector this->pointSourceData().emplace_back(std::move(psData)); - // export the lowdim coupling stencil - // we insert all vertices / elements and make it unique later - if (isBox()) - { - const auto& vertices = this->vertexIndices(bulkIdx, bulkElementIdx); - this->couplingStencils(lowDimIdx)[lowDimElementIdx].insert(this->couplingStencils(lowDimIdx)[lowDimElementIdx].end(), - vertices.begin(), vertices.end()); - } - else - { - this->couplingStencils(lowDimIdx)[lowDimElementIdx].push_back(bulkElementIdx); - } - // export the bulk coupling stencil // we insert all vertices / elements and make it unique later if (isBox()) @@ -771,11 +915,53 @@ public: this->couplingStencils(bulkIdx)[bulkElementIdx].push_back(lowDimElementIdx); } + // export bulk circle stencil (only needed for circle average) + if (useCircleAverage) + { + if (isBox()) + { + // we insert all vertices and make it unique later + for (const auto& vertices : circleCornerIndices) + { + extendedSourceStencil_.stencil()[bulkElementIdx].insert(extendedSourceStencil_.stencil()[bulkElementIdx].end(), + vertices->begin(), vertices->end()); + + } + } + else + { + extendedSourceStencil_.stencil()[bulkElementIdx].insert(extendedSourceStencil_.stencil()[bulkElementIdx].end(), + circleStencil.begin(), circleStencil.end()); + } + } } } } } + // make the circle stencil unique (for source derivatives) + for (auto&& stencil : extendedSourceStencil_.stencil()) + { + std::sort(stencil.second.begin(), stencil.second.end()); + stencil.second.erase(std::unique(stencil.second.begin(), stencil.second.end()), stencil.second.end()); + + // remove the vertices element (box) + if (isBox()) + { + const auto& indices = this->vertexIndices(bulkIdx, stencil.first); + stencil.second.erase(std::remove_if(stencil.second.begin(), stencil.second.end(), + [&](auto i){ return std::find(indices.begin(), indices.end(), i) != indices.end(); }), + stencil.second.end()); + } + // remove the own element (cell-centered) + else + { + stencil.second.erase(std::remove_if(stencil.second.begin(), stencil.second.end(), + [&](auto i){ return i == stencil.first; }), + stencil.second.end()); + } + } + // make stencils unique using namespace Dune::Hybrid; forEach(integralRange(Dune::index_constant<2>{}), [&](const auto domainIdx) @@ -849,6 +1035,8 @@ public: // \} private: + //! the extended source stencil object + EmbeddedCoupling::ExtendedSourceStencil extendedSourceStencil_; //! vector for the volume fraction of the lowdim domain in the bulk domain cells std::vector lowDimVolumeInBulkElement_; @@ -889,11 +1077,23 @@ class EmbeddedCouplingManager1d3d static constexpr bool isBox() { return GridGeometry::discMethod == DiscretizationMethod::box; } + static_assert(!isBox() && !isBox(), "The kernel coupling method is only implemented for the tpfa method"); + static_assert(Dune::Capabilities::isCartesian::Grid>::v, "The kernel coupling method is only implemented for structured grids"); + enum { bulkDim = GridView::dimension, lowDimDim = GridView::dimension, dimWorld = GridView::dimensionworld }; + + // detect if a class (the spatial params) has a kernelWidthFactor() function + template + using VariableKernelWidthDetector = decltype(std::declval().kernelWidthFactor(std::declval()...)); + + template + static constexpr bool hasKernelWidthFactor() + { return Dune::Std::is_detected::value; } + public: static constexpr EmbeddedCouplingMode couplingMode = EmbeddedCouplingMode::kernel; @@ -905,6 +1105,10 @@ public: { ParentType::init(bulkProblem, lowDimProblem, curSol); computeLowDimVolumeFractions(); + + const auto refinement = getParamFromGroup(bulkProblem->paramGroup(), "Grid.Refinement", 0); + if (refinement > 0) + DUNE_THROW(Dune::NotImplemented, "The current intersection detection may likely fail for refined grids."); } /*! @@ -952,146 +1156,133 @@ public: Dune::Timer watch; std::cout << "Initializing the point sources..." << std::endl; - // clear all internal members like pointsource vectors and stencils - // initializes the point source id counter - this->clear(); - bulkSourceIds_.clear(); - bulkSourceWeights_.clear(); - extendedSourceStencil_.clear(); - - // precompute the vertex indices for efficiency for the box method - this->precomputeVertexIndices(bulkIdx); - this->precomputeVertexIndices(lowDimIdx); + // prepare the internal data structures + prepareDataStructures_(); const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry(); const auto& lowDimGridGeometry = this->problem(lowDimIdx).gridGeometry(); - bulkSourceIds_.resize(this->gridView(bulkIdx).size(0)); - bulkSourceWeights_.resize(this->gridView(bulkIdx).size(0)); + // generate a bunch of random vectors and values for + // Monte-carlo integration on the cylinder defined by line and radius + static const int samples = getParam("MixedDimension.NumMCSamples"); + CylinderIntegration cylIntegration(samples, 1); - // intersect the bounding box trees - this->glueGrids(); + static const auto writeIntegrationPointsToFile = getParam("MixedDimension.WriteIntegrationPointsToFile", false); + if (writeIntegrationPointsToFile) + { + std::ofstream ipPointFile("kernel_points.log", std::ios::trunc); + ipPointFile << "x,y,z\n"; + } - this->pointSourceData().reserve(this->glue().size()); - this->averageDistanceToBulkCell().reserve(this->glue().size()); - const Scalar kernelWidth = getParam("MixedDimension.KernelWidth"); for (const auto& is : intersections(this->glue())) { // all inside elements are identical... const auto& inside = is.targetEntity(0); // get the intersection geometry for integrating over it const auto intersectionGeometry = is.geometry(); - // get the Gaussian quadrature rule for the local intersection - const auto& quad = Dune::QuadratureRules::rule(intersectionGeometry.type(), order); - const auto lowDimElementIdx = lowDimGridGeometry.elementMapper().index(inside); + const auto lowDimElementIdx = lowDimFvGridGeometry.elementMapper().index(inside); - // iterate over all quadrature points and place a source - // for 1d: make a new point source - // for 3d: make a new kernel volume source - for (auto&& qp : quad) - { - // compute the coupling stencils - for (int outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx) - { - const auto& outside = is.domainEntity(outsideIdx); - const auto bulkElementIdx = bulkGridGeometry.elementMapper().index(outside); + // for each intersection integrate kernel and add: + // * 1d: a new point source + // * 3d: a new kernel volume source + const auto radius = this->problem(lowDimIdx).spatialParams().radius(lowDimElementIdx); + const auto kernelWidthFactor = kernelWidthFactor_(this->problem(lowDimIdx).spatialParams(), lowDimElementIdx); + const auto kernelWidth = kernelWidthFactor*radius; + const auto a = intersectionGeometry.corner(0); + const auto b = intersectionGeometry.corner(1); + cylIntegration.setGeometry(a, b, kernelWidth); - // each quadrature point will be a point source for the sub problem - const auto globalPos = intersectionGeometry.global(qp.position()); - const auto id = this->idCounter_++; - const auto qpweight = qp.weight(); - const auto ie = intersectionGeometry.integrationElement(qp.position()); - this->pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, ie, std::vector({lowDimElementIdx})); - this->pointSources(lowDimIdx).back().setEmbeddings(is.numDomainNeighbors()); - computeBulkSource(globalPos, kernelWidth, id, lowDimElementIdx, bulkElementIdx, qpweight*ie/is.numDomainNeighbors()); + // we can have multiple 3d elements per 1d intersection, for evaluation taking one of them should be fine + for (int outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx) + { + // compute source id + // for each id we have + // (1) a point source for the 1d domain + // (2) a bulk source weight for the each element in the 3d domain (the fraction of the total source/sink for the 1d element with the point source) + // TODO: i.e. one integration of the kernel should be enough (for each of the weights it's weight/embeddings) + // (3) the flux scaling factor for each outside element, i.e. each id + const auto id = this->idCounter_++; + + // compute the weights for the bulk volume sources + const auto& outside = is.domainEntity(outsideIdx); + const auto bulkElementIdx = bulkFvGridGeometry.elementMapper().index(outside); + const auto surfaceFactor = computeBulkSource(intersectionGeometry, radius, kernelWidth, id, lowDimElementIdx, bulkElementIdx, cylIntegration, is.numDomainNeighbors()); - // pre compute additional data used for the evaluation of - // the actual solution dependent source term - PointSourceData psData; + // place a point source at the intersection center + const auto center = intersectionGeometry.center(); + this->pointSources(lowDimIdx).emplace_back(center, id, surfaceFactor, intersectionGeometry.volume(), lowDimElementIdx); + this->pointSources(lowDimIdx).back().setEmbeddings(is.numDomainNeighbors()); - if constexpr (isBox()) - { - using ShapeValues = std::vector >; - const auto lowDimGeometry = lowDimGridGeometry.element(lowDimElementIdx).geometry(); - ShapeValues shapeValues; - this->getShapeValues(lowDimIdx, lowDimGridGeometry, lowDimGeometry, globalPos, shapeValues); - psData.addLowDimInterpolation(shapeValues, this->vertexIndices(lowDimIdx, lowDimElementIdx), lowDimElementIdx); - } - else - { - psData.addLowDimInterpolation(lowDimElementIdx); - } + // pre compute additional data used for the evaluation of + // the actual solution dependent source term + PointSourceData psData; - // add data needed to compute integral over the circle - if constexpr (isBox()) - { - using ShapeValues = std::vector >; - const auto bulkGeometry = bulkGridGeometry.element(bulkElementIdx).geometry(); - ShapeValues shapeValues; - this->getShapeValues(bulkIdx, bulkGridGeometry, bulkGeometry, globalPos, shapeValues); - psData.addBulkInterpolation(shapeValues, this->vertexIndices(bulkIdx, bulkElementIdx), bulkElementIdx); - } - else - { - psData.addBulkInterpolation(bulkElementIdx); - } + // lowdim interpolation (evaluate at center) + if (isBox()) + { + using ShapeValues = std::vector >; + const auto lowDimGeometry = this->problem(lowDimIdx).fvGridGeometry().element(lowDimElementIdx).geometry(); + ShapeValues shapeValues; + this->getShapeValues(lowDimIdx, this->problem(lowDimIdx).fvGridGeometry(), lowDimGeometry, center, shapeValues); + psData.addLowDimInterpolation(shapeValues, this->vertexIndices(lowDimIdx, lowDimElementIdx), lowDimElementIdx); + } + else + { + psData.addLowDimInterpolation(lowDimElementIdx); + } - // publish point source data in the global vector - this->pointSourceData().emplace_back(std::move(psData)); + // bulk interpolation (evaluate at center) + if (isBox()) + { + using ShapeValues = std::vector >; + const auto bulkGeometry = this->problem(bulkIdx).fvGridGeometry().element(bulkElementIdx).geometry(); + ShapeValues shapeValues; + this->getShapeValues(bulkIdx, this->problem(bulkIdx).fvGridGeometry(), bulkGeometry, center, shapeValues); + psData.addBulkInterpolation(shapeValues, this->vertexIndices(bulkIdx, bulkElementIdx), bulkElementIdx); + } + else + { + psData.addBulkInterpolation(bulkElementIdx); + } - // compute average distance to bulk cell - this->averageDistanceToBulkCell().push_back(averageDistancePointGeometry(globalPos, outside.geometry())); + // publish point source data in the global vector + this->pointSourceData().emplace_back(std::move(psData)); - // export the lowdim coupling stencil - // we insert all vertices / elements and make it unique later - if (isBox()) - { - const auto& vertices = this->vertexIndices(bulkIdx, bulkElementIdx); - this->couplingStencils(lowDimIdx)[lowDimElementIdx].insert(this->couplingStencils(lowDimIdx)[lowDimElementIdx].end(), - vertices.begin(), vertices.end()); - } - else - { - this->couplingStencils(lowDimIdx)[lowDimElementIdx].push_back(bulkElementIdx); - } + // compute the average flux scaling factor + const auto outsideGeometry = outside.geometry(); + static const int fluxScalingFactorIntOrder = getParam("MixedDimension.FluxScalingIntegrationOrder", 5); + const auto& quadBulk = Dune::QuadratureRules::rule(outsideGeometry.type(), fluxScalingFactorIntOrder); + // iterate over all quadrature points + Scalar avgMinDist = 0.0; + for (auto&& qpBulk : quadBulk) + { + const auto globalPosBulk = outsideGeometry.global(qpBulk.position()); + const auto d = computeMinDistance_(a, b, globalPosBulk); + avgMinDist += d*qpBulk.weight(); } - } - } - // make extra stencils unique - for (auto&& stencil : extendedSourceStencil_.stencil()) - { - std::sort(stencil.second.begin(), stencil.second.end()); - stencil.second.erase(std::unique(stencil.second.begin(), stencil.second.end()), stencil.second.end()); + // add it to the internal data vector + this->averageDistanceToBulkCell().push_back(avgMinDist); + fluxScalingFactor_.push_back(this->problem(bulkIdx).fluxScalingFactor(avgMinDist, radius, kernelWidth)); - // remove the vertices element (box) - if (isBox()) - { - const auto& indices = this->vertexIndices(bulkIdx, stencil.first); - stencil.second.erase(std::remove_if(stencil.second.begin(), stencil.second.end(), - [&](auto i){ return std::find(indices.begin(), indices.end(), i) != indices.end(); }), - stencil.second.end()); - } - // remove the own element (cell-centered) - else - { - stencil.second.erase(std::remove_if(stencil.second.begin(), stencil.second.end(), - [&](auto i){ return i == stencil.first; }), - stencil.second.end()); + // export the lowdim coupling stencil + // we insert all vertices / elements and make it unique later + if (isBox()) + { + const auto& vertices = this->vertexIndices(bulkIdx, bulkElementIdx); + this->couplingStencils(lowDimIdx)[lowDimElementIdx].insert(this->couplingStencils(lowDimIdx)[lowDimElementIdx].end(), + vertices.begin(), vertices.end()); + } + else + { + this->couplingStencils(lowDimIdx)[lowDimElementIdx].push_back(bulkElementIdx); + } } } - // make stencils unique - using namespace Dune::Hybrid; - forEach(integralRange(Dune::index_constant<2>{}), [&](const auto domainIdx) - { - for (auto&& stencil : this->couplingStencils(domainIdx)) - { - std::sort(stencil.second.begin(), stencil.second.end()); - stencil.second.erase(std::unique(stencil.second.begin(), stencil.second.end()), stencil.second.end()); - } - }); + // make the stencils unique + makeUniqueStencil_(); if (!this->pointSources(bulkIdx).empty()) DUNE_THROW(Dune::InvalidStateException, "Kernel method shouldn't have point sources in the bulk domain but only volume sources!"); @@ -1156,95 +1347,263 @@ public: } //! return all source ids for a bulk elements - const std::vector bulkSourceIds(GridIndex eIdx) const - { return bulkSourceIds_[eIdx]; } + const std::vector& bulkSourceIds(GridIndex eIdx, int scvIdx = 0) const + { return bulkSourceIds_[eIdx][scvIdx]; } //! return all source ids for a bulk elements - const std::vector bulkSourceWeights(GridIndex eIdx) const - { return bulkSourceWeights_[eIdx]; } + const std::vector& bulkSourceWeights(GridIndex eIdx, int scvIdx = 0) const + { return bulkSourceWeights_[eIdx][scvIdx]; } + + //! The flux scaling factor for a source with id + Scalar fluxScalingFactor(std::size_t id) const + { return fluxScalingFactor_[id]; } // \} private: - //! TODO: How to optimize this? - void computeBulkSource(const GlobalPosition& globalPos, const Scalar kernelWidth, - std::size_t id, GridIndex lowDimElementIdx, GridIndex coupledBulkElementIdx, - Scalar pointSourceWeight) + //! compute the kernel distributed sources and add stencils + template + Scalar computeBulkSource(const Line& line, const Scalar radius, const Scalar kernelWidth, + std::size_t id, GridIndex lowDimElementIdx, GridIndex coupledBulkElementIdx, + const CylIntegration& cylIntegration, int embeddings) { - // make sure it is mass conservative - // i.e. the point source in the 1d domain needs to have the exact same integral as the distributed - // kernel source integral in the 3d domain. Correct the integration formula by balancing the error - // by scaling the kernel with the checksum inverse - Scalar checkSum = 0.0; - std::vector mask(this->gridView(bulkIdx).size(0), false); - for (const auto& element : elements(this->gridView(bulkIdx))) + // Monte-carlo integration on the cylinder defined by line and radius + static const auto bulkParamGroup = this->problem(bulkIdx).paramGroup(); + static const auto min = getParamFromGroup(bulkParamGroup, "Grid.LowerLeft"); + static const auto max = getParamFromGroup(bulkParamGroup, "Grid.UpperRight"); + static const auto cells = getParamFromGroup(bulkParamGroup, "Grid.Cells"); + const auto cylSamples = cylIntegration.size(); + const auto& a = line.corner(0); + const auto& b = line.corner(1); + + // optionally write debugging / visual output of the integration points + static const auto writeIntegrationPointsToFile = getParam("MixedDimension.WriteIntegrationPointsToFile", false); + if (writeIntegrationPointsToFile) { - Scalar weight = 0.0; - const auto geometry = element.geometry(); - const auto& quad = Dune::QuadratureRules::rule(geometry.type(), 3); - for (auto&& qp : quad) + std::ofstream ipPointFile("kernel_points.log", std::ios::app); + for (int i = 0; i < cylSamples; ++i) { - const auto qpweight = qp.weight(); - const auto ie = geometry.integrationElement(qp.position()); - weight += evalKernel(globalPos, geometry.global(qp.position()), kernelWidth)*qpweight*ie; + const auto& point = cylIntegration.getIntegrationPoint(i); + const auto [hasIntersection, bulkElementIdx] = intersectingEntityCartesian(point, min, max, cells); + if (hasIntersection) + ipPointFile << point[0] << "," << point[1] << "," << point[2] << '\n'; } + } - if (weight > 1e-13) + Scalar integral = 0.0; + for (int i = 0; i < cylSamples; ++i) + { + const auto& point = cylIntegration.getIntegrationPoint(i); + // TODO: below only works for Cartesian grids with ijk numbering (e.g. level 0 YaspGrid (already fails when refined)) + // more general is the bounding box tree solution which always works, however it's much slower + //const auto bulkIndices = intersectingEntities(point, this->problem(bulkIdx).fvGridGeometry().boundingBoxTree(), true); + if (const auto [hasIntersection, bulkElementIdx] = intersectingEntityCartesian(point, min, max, cells); hasIntersection) { - const auto& bulkGridGeometry = this->problem(bulkIdx).gridGeometry(); - const auto bulkElementIdx = bulkGridGeometry.elementMapper().index(element); - bulkSourceIds_[bulkElementIdx].push_back(id); - bulkSourceWeights_[bulkElementIdx].push_back(weight*pointSourceWeight); - mask[bulkElementIdx] = true; - - // add lowDim dofs that the source is related to to the bulk stencil - if (isBox()) + assert(bulkElementIdx < this->problem(bulkIdx).fvGridGeometry().gridView().size(0)); + const auto localWeight = evalConstKernel_(a, b, point, radius, kernelWidth)*cylIntegration.integrationElement(i)/Scalar(embeddings); + integral += localWeight; + if (!bulkSourceIds_[bulkElementIdx][0].empty() && id == bulkSourceIds_[bulkElementIdx][0].back()) { - const auto& vertices = this->vertexIndices(lowDimIdx, lowDimElementIdx); - this->couplingStencils(bulkIdx)[bulkElementIdx].insert(this->couplingStencils(bulkIdx)[bulkElementIdx].end(), - vertices.begin(), vertices.end()); - + bulkSourceWeights_[bulkElementIdx][0].back() += localWeight; } else { - this->couplingStencils(bulkIdx)[bulkElementIdx].push_back(lowDimElementIdx); + bulkSourceIds_[bulkElementIdx][0].emplace_back(id); + bulkSourceWeights_[bulkElementIdx][0].emplace_back(localWeight); + addBulkSourceStencils_(bulkElementIdx, lowDimElementIdx, coupledBulkElementIdx); } + } + } + + // the surface factor (which fraction of the source is inside the domain and needs to be considered) + const auto length = (a-b).two_norm()/Scalar(embeddings); + return integral/length; + } + + void prepareDataStructures_() + { + // clear all internal members like pointsource vectors and stencils + // initializes the point source id counter + this->clear(); + bulkSourceIds_.clear(); + bulkSourceWeights_.clear(); + extendedSourceStencil_.stencil().clear(); + + // precompute the vertex indices for efficiency for the box method + this->preComputeVertexIndices(bulkIdx); + this->preComputeVertexIndices(lowDimIdx); + + bulkSourceIds_.resize(this->gridView(bulkIdx).size(0)); + bulkSourceWeights_.resize(this->gridView(bulkIdx).size(0)); + + // intersect the bounding box trees + this->glueGrids(); - // tpfa - extendedSourceStencil_.stencil()[bulkElementIdx].push_back(coupledBulkElementIdx); + // reserve memory for data + this->pointSourceData().reserve(this->glue().size()); + this->averageDistanceToBulkCell().reserve(this->glue().size()); + fluxScalingFactor_.reserve(this->glue().size()); - // compute check sum -> should sum up to 1.0 to be mass conservative - checkSum += weight; + // reserve memory for stencils + const auto numBulkElements = this->gridView(bulkIdx).size(0); + for (GridIndex bulkElementIdx = 0; bulkElementIdx < numBulkElements; ++bulkElementIdx) + { + this->couplingStencils(bulkIdx)[bulkElementIdx].reserve(50); + extendedSourceStencil_.stencil()[bulkElementIdx].reserve(50); + bulkSourceIds_[bulkElementIdx][0].reserve(20); + bulkSourceWeights_[bulkElementIdx][0].reserve(20); + } + } + + //! Make the stencils unique + void makeUniqueStencil_() + { + // make extra stencils unique + for (auto&& stencil : extendedSourceStencil_.stencil()) + { + std::sort(stencil.second.begin(), stencil.second.end()); + stencil.second.erase(std::unique(stencil.second.begin(), stencil.second.end()), stencil.second.end()); + + // remove the vertices element (box) + if (isBox()) + { + const auto& indices = this->vertexIndices(bulkIdx, stencil.first); + stencil.second.erase(std::remove_if(stencil.second.begin(), stencil.second.end(), + [&](auto i){ return std::find(indices.begin(), indices.end(), i) != indices.end(); }), + stencil.second.end()); } + // remove the own element (cell-centered) + else + { + stencil.second.erase(std::remove_if(stencil.second.begin(), stencil.second.end(), + [&](auto i){ return i == stencil.first; }), + stencil.second.end()); + } + } + + // make stencils unique + using namespace Dune::Hybrid; + forEach(integralRange(Dune::index_constant<2>{}), [&](const auto domainIdx) + { + for (auto&& stencil : this->couplingStencils(domainIdx)) + { + std::sort(stencil.second.begin(), stencil.second.end()); + stencil.second.erase(std::unique(stencil.second.begin(), stencil.second.end()), stencil.second.end()); + } + }); + } + + //! add additional stencil entries for the bulk element + void addBulkSourceStencils_(GridIndex bulkElementIdx, GridIndex coupledLowDimElementIdx, GridIndex coupledBulkElementIdx) + { + // add the lowdim element to the coupling stencil of this bulk element + if (isBox()) + { + const auto& vertices = this->vertexIndices(lowDimIdx, coupledLowDimElementIdx); + this->couplingStencils(bulkIdx)[bulkElementIdx].insert(this->couplingStencils(bulkIdx)[bulkElementIdx].end(), + vertices.begin(), vertices.end()); + + } + else + { + auto& s = this->couplingStencils(bulkIdx)[bulkElementIdx]; + s.push_back(coupledLowDimElementIdx); + } + + // the extended source stencil, every 3d element with a source is coupled to + // the element/dofs where the 3d quantities are measured + if (isBox()) + { + const auto& vertices = this->vertexIndices(bulkIdx, coupledBulkElementIdx); + extendedSourceStencil_.stencil()[bulkElementIdx].insert(extendedSourceStencil_.stencil()[bulkElementIdx].end(), + vertices.begin(), vertices.end()); + } + else + { + auto& s = extendedSourceStencil_.stencil()[bulkElementIdx]; + s.push_back(coupledBulkElementIdx); } + } + + //! a cylindrical kernel around the segment a->b + Scalar evalConstKernel_(const GlobalPosition& a, + const GlobalPosition& b, + const GlobalPosition& point, + const Scalar R, + const Scalar rho) const noexcept + { + // projection of point onto line a + t*(b-a) + const auto ab = b - a; + const auto t = (point - a)*ab/ab.two_norm2(); - for (GridIndex eIdx = 0; eIdx < bulkSourceWeights_.size(); ++eIdx) - if (mask[eIdx]) bulkSourceWeights_[eIdx].back() /= checkSum; + // return 0 if we are outside cylinder + if (t < 0.0 || t > 1.0) + return 0.0; + + // compute distance + auto proj = a; proj.axpy(t, ab); + const auto r = (proj - point).two_norm(); + + if (r > rho) + return 0.0; - // balance error of the quadrature rule -> TODO: what to do at boundaries - // const auto diff = 1.0 - checkSum; - // std::cout << "Integrated kernel with integration error of " << diff << std::endl; + return 1.0/(M_PI*rho*rho); } - //! an isotropic cubic kernel with derivatives 0 at r=origin and r=width and domain integral 1 - inline Scalar evalKernel(const GlobalPosition& origin, - const GlobalPosition& pos, - const Scalar width) const noexcept + //! a cylindrical kernel around the segment a->b + Scalar evalCubicKernel_(const GlobalPosition& a, + const GlobalPosition& b, + const GlobalPosition& point, + const Scalar R, + const Scalar rho) const noexcept { - const auto r = (pos-origin).two_norm(); - const auto r2 = r*r; - const auto r3 = r2*r; + // projection of point onto line a + t*(b-a) + const auto ab = b - a; + const auto t = (point - a)*ab/ab.two_norm2(); - if (r > width) + // return 0 if we are outside cylinder + if (t < 0.0 || t > 1.0) return 0.0; - const Scalar w2 = width*width; - const Scalar w3 = w2*width; - const Scalar k = 15.0/(4*M_PI*w3); - const Scalar a = 2.0/w3; - const Scalar b = 3.0/w2; + // compute distance + auto proj = a; proj.axpy(t, ab); + const auto r = (proj - point).two_norm(); + + if (r > rho) + return 0.0; - return k*(a*r3 - b*r2 + 1.0); + const auto r2 = r*r; const auto r3 = r*r*r; + const auto rho2 = rho*rho; const auto rho3 = rho2*rho; + return 10.0*(2.0*r3/rho3 - 3.0*r2/rho2 + 1.0)/(3.0*M_PI*rho2); + } + + Scalar computeMinDistance_(const GlobalPosition& a, const GlobalPosition& b, const GlobalPosition& x3d) const noexcept + { + // r is the minimum distance to the line through a and b + const auto ab = b - a; + const auto t = (x3d - a)*ab/ab.two_norm2(); + auto proj = a; proj.axpy(t, ab); + const auto r = (proj - x3d).two_norm(); + return r; + } + + /*! + * \brief Get the kernel width factor from the spatial params (if possible) + */ + template + auto kernelWidthFactor_(const SpatialParams& spatialParams, unsigned int eIdx) + -> std::enable_if_t(), Scalar> + { return spatialParams.kernelWidthFactor(eIdx); } + + /*! + * \brief Get the kernel width factor (constant) from the input file (if possible) + */ + template + auto kernelWidthFactor_(const SpatialParams& spatialParams, unsigned int eIdx) + -> std::enable_if_t(), Scalar> + { + static const Scalar kernelWidthFactor = getParam("MixedDimension.KernelWidthFactor"); + return kernelWidthFactor; } //! the extended source stencil object for kernel coupling @@ -1252,9 +1611,11 @@ private: //! vector for the volume fraction of the lowdim domain in the bulk domain cells std::vector lowDimVolumeInBulkElement_; //! kernel sources to integrate for each bulk element - std::vector> bulkSourceIds_; + std::vector, isBox() ? 1<> bulkSourceIds_; //! the integral of the kernel for each point source / integration point, i.e. weight for the source - std::vector> bulkSourceWeights_; + std::vector, isBox() ? 1<> bulkSourceWeights_; + //! the flux scaling factor for the respective source id + std::vector fluxScalingFactor_; }; } // end namespace Dumux