From 4bbc4feb727dcb59fe2ed6d9dccb392ad8ff3fc4 Mon Sep 17 00:00:00 2001
From: "Dennis.Glaeser" <dennis.glaeser@iws.uni-stuttgart.de>
Date: Wed, 26 Jun 2019 14:20:22 +0200
Subject: [PATCH] [projector] reduce system size for solve where possible
---
dumux/discretization/projection/projector.hh | 215 +++++++++++++++++--
1 file changed, 198 insertions(+), 17 deletions(-)
diff --git a/dumux/discretization/projection/projector.hh b/dumux/discretization/projection/projector.hh
index 9b74c78099..0222a16b5d 100644
--- a/dumux/discretization/projection/projector.hh
+++ b/dumux/discretization/projection/projector.hh
@@ -87,7 +87,40 @@ public:
Projector(Matrix&& massMatrix, Matrix&& projectionMatrix)
: massMat_(std::move(massMatrix))
, projMat_(std::move(projectionMatrix))
- {}
+ , numDofsTarget_(massMat_.N())
+ {
+ if (massMat_.N() != projMat_.N())
+ DUNE_THROW(Dune::InvalidStateException, "Matrix row size mismatch: " << massMat_.N() << " vs " << projMat_.N());
+ }
+
+ /*!
+ * \brief Constructor for projection into a target space that occupies
+ * a larger geometric region than the domain space. In this case,
+ * the mass matrix can be chosen such that is is solved for only
+ * those dofs which will be populated with values from the projection.
+ * This requires an additional index map that maps the entries of the
+ * projected solution into the solution vector for the target space.
+ * Furthermore, the number of degrees of freedom must be specified to
+ * set up the coefficient vector with correct size for the target space.
+ */
+ Projector(Matrix&& massMatrix,
+ Matrix&& projectionMatrix,
+ std::vector<std::size_t>&& indexMap,
+ std::size_t numDofsTarget)
+ : massMat_(std::move(massMatrix))
+ , projMat_(std::move(projectionMatrix))
+ , indexMapTarget_(std::move(indexMap))
+ , numDofsTarget_(numDofsTarget)
+ {
+ if (indexMapTarget_.size() != massMat_.N())
+ DUNE_THROW(Dune::InvalidStateException, "Target index map size mismatch: " << indexMapTarget_.size() << " vs " << massMat_.N());
+
+ if (massMat_.N() != projMat_.N())
+ DUNE_THROW(Dune::InvalidStateException, "Matrix row size mismatch: " << massMat_.N() << " vs " << projMat_.N());
+
+ if (*std::max_element(indexMapTarget_.begin(), indexMapTarget_.end()) > numDofsTarget_)
+ DUNE_THROW(Dune::InvalidStateException, "Index map exceeds provided number of dofs in target domain!");
+ }
/*!
* \brief Project a solution u into up
@@ -99,7 +132,7 @@ public:
{
// be picky about size of u
if ( u.size() != projMat_.M())
- DUNE_THROW(Dune::InvalidStateException, "Vector size mismatch" );
+ DUNE_THROW(Dune::InvalidStateException, "Vector size mismatch");
Dune::BlockVector<BlockType> up(massMat_.N());
@@ -111,6 +144,18 @@ public:
solver.setMaxIter(params.maxIterations);
solver.solve(massMat_, up, rhs);
+ // if target space occupies a larger region, fill missing entries with zero
+ if (!indexMapTarget_.empty())
+ {
+ Dune::BlockVector<BlockType> result(numDofsTarget_);
+
+ result = 0.0;
+ for (std::size_t i = 0; i < indexMapTarget_.size(); ++i)
+ result[indexMapTarget_[i]] = up[i];
+
+ return result;
+ }
+
return up;
}
@@ -122,7 +167,7 @@ public:
template< class BlockType, std::enable_if_t<!std::is_convertible<BlockType, ScalarType>::value, int> = 0 >
Dune::BlockVector<BlockType> project(const Dune::BlockVector<BlockType>& u, const Params& params = Params{}) const
{
- Dune::BlockVector<BlockType> result(massMat_.N());
+ Dune::BlockVector<BlockType> result(numDofsTarget_);
for (int pvIdx = 0; pvIdx < BlockType::size(); ++pvIdx)
{
@@ -146,6 +191,9 @@ public:
private:
Matrix massMat_;
Matrix projMat_;
+
+ std::vector<std::size_t> indexMapTarget_;
+ std::size_t numDofsTarget_;
};
/*!
@@ -167,6 +215,75 @@ public:
// Projector construction details
namespace Impl {
+/*!
+ * \brief Reduces a mass matrix and projection matrix such that they are composed
+ * of only those dofs that actually take part in the projection. Simultaneously,
+ * a container with the index map into the complete target space is filled so that
+ * the entries after projection can be assigned to the corresponding dof in the
+ * overall target space.
+ */
+template<class Matrix>
+void setupReducedMatrices(const Matrix& massMatrix, const Matrix& projMatrix, const std::vector<bool>& dofIsVoid,
+ Matrix& reducedM, Matrix& reducedP, std::vector<std::size_t>& expansionMap)
+{
+ const std::size_t numNonVoidDofs = std::count_if(dofIsVoid.begin(), dofIsVoid.end(), [] (bool v) { return !v; });
+
+ // reduce matrices to only dofs that take part and create index map
+ std::vector<std::size_t> reductionMap(massMatrix.N());
+ expansionMap.resize(numNonVoidDofs);
+
+ std::size_t idxInReducedSpace = 0;
+ for (std::size_t dofIdx = 0; dofIdx < dofIsVoid.size(); ++dofIdx)
+ if (!dofIsVoid[dofIdx])
+ {
+ reductionMap[dofIdx] = idxInReducedSpace;
+ expansionMap[idxInReducedSpace] = dofIdx;
+ idxInReducedSpace++;
+ }
+
+ // create reduced M/P matrix
+ Dune::MatrixIndexSet patternMReduced, patternPReduced;
+ patternMReduced.resize(numNonVoidDofs, numNonVoidDofs);
+ patternPReduced.resize(numNonVoidDofs, projMatrix.M());
+ for (auto rowIt = massMatrix.begin(); rowIt != massMatrix.end(); ++rowIt)
+ if (!dofIsVoid[rowIt.index()])
+ {
+ const auto reducedRowIdx = reductionMap[rowIt.index()];
+ for (auto colIt = (*rowIt).begin(); colIt != (*rowIt).end(); ++colIt)
+ if (!dofIsVoid[colIt.index()])
+ patternMReduced.add(reducedRowIdx, reductionMap[colIt.index()]);
+ }
+
+ for (auto rowIt = projMatrix.begin(); rowIt != projMatrix.end(); ++rowIt)
+ if (!dofIsVoid[rowIt.index()])
+ {
+ const auto reducedRowIdx = reductionMap[rowIt.index()];
+ for (auto colIt = (*rowIt).begin(); colIt != (*rowIt).end(); ++colIt)
+ patternPReduced.add(reducedRowIdx, colIt.index());
+ }
+
+ patternMReduced.exportIdx(reducedM);
+ patternPReduced.exportIdx(reducedP);
+
+ // fill matrix entries
+ for (auto rowIt = massMatrix.begin(); rowIt != massMatrix.end(); ++rowIt)
+ if (!dofIsVoid[rowIt.index()])
+ {
+ const auto reducedRowIdx = reductionMap[rowIt.index()];
+ for (auto colIt = (*rowIt).begin(); colIt != (*rowIt).end(); ++colIt)
+ if (!dofIsVoid[colIt.index()])
+ reducedM[reducedRowIdx][reductionMap[colIt.index()]] = *colIt;
+ }
+
+ for (auto rowIt = projMatrix.begin(); rowIt != projMatrix.end(); ++rowIt)
+ if (!dofIsVoid[rowIt.index()])
+ {
+ const auto reducedRowIdx = reductionMap[rowIt.index()];
+ for (auto colIt = (*rowIt).begin(); colIt != (*rowIt).end(); ++colIt)
+ reducedP[reducedRowIdx][colIt.index()] = *colIt;
+ }
+}
+
/*!
* \brief Creates a projector class between two function space bases
* \tparam doBidirectional If false, the backward projection matrix is not assembled
@@ -325,26 +442,90 @@ makeProjectorPair(const FEBasisDomain& feBasisDomain,
}
}
- // On those dofs that to not take part in any intersection we will have zeroes
- // in the mass matrices. The right hand side will be zero because the projection
- // matrix has no entries for those dofs as there was no intersection. Thus, we set
- // 1.0 on the diagonal for those dofs, such that after projection, the projected
- // solution is 0 on them.
- for (std::size_t dofIdx = 0; dofIdx < forwardM.N(); ++dofIdx)
- if (forwardM[dofIdx][dofIdx] == 0.0)
- forwardM[dofIdx][dofIdx] = 1.0;
+ // determine the dofs that do not take part in intersections
+ std::vector<bool> isVoidTarget(forwardM.N(), false);
+ for (std::size_t dofIdxTarget = 0; dofIdxTarget < forwardM.N(); ++dofIdxTarget)
+ if (forwardM[dofIdxTarget][dofIdxTarget] == 0.0)
+ isVoidTarget[dofIdxTarget] = true;
+ std::vector<bool> isVoidDomain;
if (doBidirectional)
{
- for (std::size_t dofIdx = 0; dofIdx < backwardM.N(); ++dofIdx)
- if (backwardM[dofIdx][dofIdx] == 0.0)
- backwardM[dofIdx][dofIdx] = 1.0;
+ isVoidDomain.resize(backwardM.N(), false);
+ for (std::size_t dofIdxDomain = 0; dofIdxDomain < backwardM.N(); ++dofIdxDomain)
+ if (backwardM[dofIdxDomain][dofIdxDomain] == 0.0)
+ isVoidDomain[dofIdxDomain] = true;
}
- ForwardProjector fw(std::move(forwardM), std::move(forwardP));
- BackwardProjector bw(std::move(backwardM), std::move(backwardP));
+ const bool hasVoidTarget = std::any_of(isVoidTarget.begin(), isVoidTarget.end(), [] (bool v) { return v; });
+ const bool hasVoidDomain = std::any_of(isVoidDomain.begin(), isVoidDomain.end(), [] (bool v) { return v; });
+ if (!hasVoidDomain && !hasVoidTarget)
+ {
+ return std::make_pair(ForwardProjector(std::move(forwardM), std::move(forwardP)),
+ BackwardProjector(std::move(backwardM), std::move(backwardP)));
+ }
+ else if (!hasVoidDomain && hasVoidTarget)
+ {
+ std::vector<std::size_t> expansionMapTarget;
+ ForwardProjectionMatrix forwardMReduced, forwardPReduced;
+ setupReducedMatrices(forwardM, forwardP, isVoidTarget,
+ forwardMReduced, forwardPReduced, expansionMapTarget);
+
+ return std::make_pair( ForwardProjector(std::move(forwardMReduced),
+ std::move(forwardPReduced),
+ std::move(expansionMapTarget),
+ forwardM.N()),
+ BackwardProjector(std::move(backwardM), std::move(backwardP)) );
+ }
+ else if (hasVoidDomain && !hasVoidTarget)
+ {
+ if (doBidirectional)
+ {
+ std::vector<std::size_t> expansionMapDomain;
+ BackwardProjectionMatrix backwardMReduced, backwardPReduced;
+ setupReducedMatrices(backwardM, backwardP, isVoidDomain,
+ backwardMReduced, backwardPReduced, expansionMapDomain);
+
+ return std::make_pair( ForwardProjector(std::move(forwardM), std::move(forwardP)),
+ BackwardProjector(std::move(backwardMReduced),
+ std::move(backwardPReduced),
+ std::move(expansionMapDomain),
+ backwardM.N()) );
+ }
+ else
+ return std::make_pair( ForwardProjector(std::move(forwardM), std::move(forwardP)),
+ BackwardProjector(std::move(backwardM), std::move(backwardP)) );
+ }
+ else
+ {
+ std::vector<std::size_t> expansionMapTarget;
+ ForwardProjectionMatrix forwardMReduced, forwardPReduced;
+ setupReducedMatrices(forwardM, forwardP, isVoidTarget,
+ forwardMReduced, forwardPReduced, expansionMapTarget);
- return std::make_pair(std::move(fw), std::move(bw));
+ if (doBidirectional)
+ {
+ std::vector<std::size_t> expansionMapDomain;
+ BackwardProjectionMatrix backwardMReduced, backwardPReduced;
+ setupReducedMatrices(backwardM, backwardP, isVoidDomain,
+ backwardMReduced, backwardPReduced, expansionMapDomain);
+
+ return std::make_pair( ForwardProjector(std::move(forwardMReduced),
+ std::move(forwardPReduced),
+ std::move(expansionMapTarget),
+ forwardM.N()),
+ BackwardProjector(std::move(backwardMReduced),
+ std::move(backwardPReduced),
+ std::move(expansionMapDomain),
+ backwardM.N()) );
+ }
+ else
+ return std::make_pair( ForwardProjector(std::move(forwardMReduced),
+ std::move(forwardPReduced),
+ std::move(expansionMapTarget),
+ forwardM.N()),
+ BackwardProjector(std::move(backwardM), std::move(backwardP)) );
+ }
}
} // end namespace Implementation
--
GitLab