From 727f836fbb397374399ade100084ed6980a159ce Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@mi.fu-berlin.de>
Date: Fri, 5 Jul 2019 17:52:01 +0200
Subject: [PATCH] .
---
.../levelpatchpreconditioner.hh | 7 ++
.../multilevelpatchpreconditioner.hh | 38 ++++++++--
.../preconditioners/nbodycontacttransfer.cc | 58 +++++++++------
src/spatial-solving/tnnmg/CMakeLists.txt | 1 +
src/spatial-solving/tnnmg/linearcorrection.hh | 74 +++++++++++++++++++
src/spatial-solving/tnnmg/linearization.hh | 8 +-
src/spatial-solving/tnnmg/tnnmgstep.hh | 31 ++++----
7 files changed, 168 insertions(+), 49 deletions(-)
create mode 100644 src/spatial-solving/tnnmg/linearcorrection.hh
diff --git a/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh b/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
index 52d5f1d1..a8216abd 100644
--- a/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
+++ b/src/spatial-solving/preconditioners/levelpatchpreconditioner.hh
@@ -176,6 +176,13 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
return patches_.size();
}
+ size_t level() const {
+ return level_;
+ }
+
+ size_t fineLevel() const {
+ return fineContactNetwork_.level();
+ }
};
#endif
diff --git a/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh b/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
index 3d1f9afc..9bab0af7 100644
--- a/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
+++ b/src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh
@@ -9,9 +9,8 @@
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
-#include "nbodytransfer.hh"
+#include "nbodycontacttransfer.hh"
#include "levelpatchpreconditioner.hh"
-#include <dune/faultnetworks/dgmgtransfer.hh>
template <class ContactNetwork, class PatchSolver, class MatrixType, class VectorType>
class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
@@ -24,7 +23,9 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
using LevelGlobalPreconditioner = LevelPatchPreconditioner<LevelContactNetwork, PatchSolver, MatrixType, VectorType>;
using BitVector = typename LinearIterationStep<MatrixType, VectorType>::BitVector;
- using MGTransfer = ;
+ using MGTransfer = NBodyContactTransfer<ContactNetwork, VectorType>;
+
+ static const int dim = ContactNetwork::dim;
const ContactNetwork& contactNetwork_;
const Dune::BitSetVector<1>& activeLevels_;
@@ -40,6 +41,7 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
PatchSolver coarseSolver_;
+ std::vector<BitVector> recompute_;
std::vector<std::shared_ptr<MGTransfer>> mgTransfer_; //TODO
public:
@@ -54,7 +56,6 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
// init level patch preconditioners and multigrid transfer
levelPatchPreconditioners_.resize(0);
- mgTransfer_.resize(0);
int maxLevel = -1;
for (size_t i=0; i<activeLevels_.size(); i++) {
@@ -89,13 +90,38 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
levelRhs_.resize(size());
// init multigrid transfer
- for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
- mgTransfer_.push_back(std::make_shared<MGTransfer>(...));
+ mgTransfer_.resize(levelPatchPreconditioners_.size());
+ for (size_t i=0; i<mgTransfer_.size(); i++) {
+ mgTransfer_[i] = std::make_shared<MGTransfer>();
+ mgTransfer_.setup(contactNetwork_, levelPatchPreconditioners_[i]->level(), levelPatchPreconditioners_[i]->fineLevel());
}
+
+ //TODO
mgTransfer_.push_back(std::make_shared<MGTransfer>(coarseGlobalPreconditioner_->basis(), *allFaultBasis_));
+
+ // Remove any recompute filed so that initially the full transferoperator is assembled
+ for (size_t i=0; i<mgTransfer_.size(); i++)
+ std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->setRecomputeBitField(nullptr);
+
+ // specify the subset of entries to be reassembled at each iteration
+ recompute_.resize(size());
+ recompute_[recompute_.size()-1] = contactNetwork_.nBodyAssembler().totalHasObstacle_;
+ for (size_t i=mgTransfer_.size(); i>=1; i--)
+ mgTransfer_[i-1]->restrict(recompute_[i], recompute_[i-1]);
+
+ for (size_t i=0; i<this->mgTransfer_.size(); i++)
+ std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->setRecomputeBitField(&recompute_[i]);
+
}
+ void setIgnore(const BitVector& i) override {
+ this->setIgnore(i);
+
+ auto mgTransfer = std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(this->mgTransfer_[level-1]);
+ mgTransfer->setCriticalBitField(&critical);
+ }
+
void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) override {
this->setProblem(mat, x, rhs);
diff --git a/src/spatial-solving/preconditioners/nbodycontacttransfer.cc b/src/spatial-solving/preconditioners/nbodycontacttransfer.cc
index 961ee3f5..0f0a58c5 100644
--- a/src/spatial-solving/preconditioners/nbodycontacttransfer.cc
+++ b/src/spatial-solving/preconditioners/nbodycontacttransfer.cc
@@ -11,20 +11,7 @@
#include <dune/solvers/transferoperators/densemultigridtransfer.hh>
template<class ContactNetwork, class VectorType>
-void NBodyContactTransfer<ContactNetwork, VectorType>::setup(const ContactNetwork& contactNetwork, const int coarseLevel, const int fineLevel)
-
- (const std::vector<const GridType*>& grids, int colevel,
- const std::vector<const MatrixType*>& mortarTransferOperator,
- const CoordSystemVector& fineLocalCoordSystems,
- const std::vector<const BitSetVector<1>*>& fineHasObstacle,
- const std::vector<std::array<int,2> >& gridIdx)
-
-transfer->setup(*grids[0], *grids[1], i,
- contactAssembler.contactCoupling_[0]->mortarLagrangeMatrix(),
- contactAssembler.localCoordSystems_,
- *contactAssembler.contactCoupling_[0]->nonmortarBoundary().getVertices());
-
-{
+void NBodyContactTransfer<ContactNetwork, VectorType>::setup(const ContactNetwork& contactNetwork, const int coarseLevel, const int fineLevel) {
const size_t nBodies = contactNetwork.nBodies();
const size_t nCouplings = contactNetwork.nCouplings();
@@ -69,23 +56,48 @@ transfer->setup(*grids[0], *grids[1], i,
// Combine the submatrices
// ///////////////////////////////////
- if (fineLevel == contactNetwork.maxLevel())
- combineSubMatrices(submat, mortarTransferOperator, fineLocalCoordSystems, fineHasObstacle, gridIdx);
- else
- Dune::MatrixVector::BlockMatrixView<MatrixType>::setupBlockMatrix(submat, this->matrix_);
+ /*
+ (const std::vector<const GridType*>& grids, int colevel,
+ const std::vector<const MatrixType*>& mortarTransferOperator,
+ const CoordSystemVector& fineLocalCoordSystems,
+ const std::vector<const BitSetVector<1>*>& fineHasObstacle,
+ const std::vector<std::array<int,2> >& gridIdx)
- for (size_t i=0; i<nBodies; i++) {
- if (fL[i]==0) {
- delete(submat[i]);
+ transfer->setup(*grids[0], *grids[1], i,
+ contactAssembler.contactCoupling_[0]->mortarLagrangeMatrix(),
+ contactAssembler.localCoordSystems_,
+ *);
+
+ */
+
+ if (fineLevel == contactNetwork.size()-1) {
+ const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+ const auto& contactCouplings = nBodyAssembler.getContactCouplings();
+
+ const std::vector<const MatrixType*> mortarTransferOperators(nBodyAssembler.nCouplings());
+ const std::vector<const BitSetVector<1>*> fineHasObstacle(nBodyAssembler.nCouplings());
+ const std::vector<std::array<int,2> > gridIdx(nBodyAssembler.nCouplings());
+
+ for (size_t i=0; i<nBodyAssembler.nCouplings(); i++) {
+ mortarTransferOperators[i] = &contactCouplings[i]->mortarLagrangeMatrix();
+ fineHasObstacle[i] = contactCouplings[i]->nonmortarBoundary().getVertices();
+ gridIdx[i] = nBodyAssembler.getCoupling(i).gridIdx_;
}
+
+ combineSubMatrices(submat, mortarTransferOperators, nBodyAssembler.getLocalCoordSystems(), fineHasObstacle, gridIdx);
+ } else {
+ Dune::MatrixVector::BlockMatrixView<MatrixType>::setupBlockMatrix(submat, this->matrix_);
}
- // Delete separate transfer objects
for (size_t i=0; i<nBodies; i++) {
+ if (fL[i] <= cL[i]) {
+ delete(submat[i]);
+ }
delete(gridTransfer[i]);
}
}
+/*
template<class VectorType>
template<class GridType>
void NBodyContactTransfer<VectorType>::setupHierarchy(std::vector<std::shared_ptr<NonSmoothNewtonContactTransfer<VectorType> > >& mgTransfers,
@@ -159,7 +171,7 @@ void NBodyContactTransfer<VectorType>::setupHierarchy(std::vector<std::shared_pt
for (size_t j=0; j<gridTransfer[i].size(); j++)
delete(gridTransfer[i][j]);
}
-
+*/
template<class VectorType>
void NBodyContactTransfer<VectorType>::combineSubMatrices(const std::vector<const MatrixType*>& submat,
diff --git a/src/spatial-solving/tnnmg/CMakeLists.txt b/src/spatial-solving/tnnmg/CMakeLists.txt
index 25682fc9..698f4e86 100644
--- a/src/spatial-solving/tnnmg/CMakeLists.txt
+++ b/src/spatial-solving/tnnmg/CMakeLists.txt
@@ -1,6 +1,7 @@
add_custom_target(dune-tectonic_spatial-solving_tnnmg_src SOURCES
functional.hh
linearization.hh
+ linearcorrection.hh
linesearchsolver.hh
localbisectionsolver.hh
localproblem.hh
diff --git a/src/spatial-solving/tnnmg/linearcorrection.hh b/src/spatial-solving/tnnmg/linearcorrection.hh
new file mode 100644
index 00000000..2f99d1ad
--- /dev/null
+++ b/src/spatial-solving/tnnmg/linearcorrection.hh
@@ -0,0 +1,74 @@
+#ifndef DUNE_TNNMG_ITERATIONSTEPS_LINEARCORRECTION_HH
+#define DUNE_TNNMG_ITERATIONSTEPS_LINEARCORRECTION_HH 1
+
+#include <functional>
+#include <memory>
+
+#include <dune/solvers/iterationsteps/lineariterationstep.hh>
+#include <dune/solvers/solvers/iterativesolver.hh>
+#include <dune/solvers/solvers/linearsolver.hh>
+#include <dune/solvers/common/canignore.hh>
+
+/**
+ * \brief linear correction step for use by \ref TNNMGStep
+ *
+ * The function object should solve the linear equation \f$ A x = b \f$.
+ * The ignore bitfield identifies the dofs to be truncated by the
+ * linear solver. The truncation dofs are passed by setting its ignore
+ * member.
+ * Be aware that full rows or columns of `A` might contain only zeroes.
+ */
+template<typename Matrix, typename Vector>
+using LinearCorrection = std::function<void(const Matrix& A, Vector& x, const Vector& b, const Dune::Solvers::DefaultBitVector_t<Vector>& ignore)>;
+
+template<typename Matrix, typename Vector>
+LinearCorrection<Matrix, Vector>
+makeLinearCorrection(std::shared_ptr< Dune::Solvers::LinearSolver<Matrix, Vector> > linearSolver)
+{
+ return [=](const Matrix& A, Vector& x, const Vector& b, const Dune::Solvers::DefaultBitVector_t<Vector>& ignore) {
+
+ auto canIgnoreCast = std::dynamic_pointer_cast<CanIgnore<Dune::Solvers::DefaultBitVector_t<Vector>>>( linearSolver );
+ if (canIgnoreCast)
+ canIgnoreCast->setIgnore(ignore);
+ else
+ DUNE_THROW(Dune::Exception, "LinearCorrection: linearSolver cannot set ignore member!");
+
+ linearSolver->setProblem(A, x, b);
+ linearSolver->preprocess();
+ linearSolver->solve();
+ };
+}
+
+template<typename Matrix, typename Vector>
+LinearCorrection<Matrix, Vector>
+makeLinearCorrection(std::shared_ptr< Dune::Solvers::IterativeSolver<Vector> > iterativeSolver)
+{
+ return [=](const Matrix& A, Vector& x, const Vector& b, const Dune::Solvers::DefaultBitVector_t<Vector>& ignore) {
+ using LinearIterationStep = Dune::Solvers::LinearIterationStep<Matrix, Vector>;
+
+ auto linearIterationStep = dynamic_cast<LinearIterationStep*>(&iterativeSolver->getIterationStep());
+ if (not linearIterationStep)
+ DUNE_THROW(Dune::Exception, "iterative solver must use a linear iteration step");
+
+ linearIterationStep->setIgnore(ignore);
+ linearIterationStep->setProblem(A, x, b);
+ iterativeSolver->preprocess();
+ iterativeSolver->solve();
+ };
+}
+
+template<typename Matrix, typename Vector>
+LinearCorrection<Matrix, Vector>
+makeLinearCorrection(std::shared_ptr< Dune::Solvers::LinearIterationStep<Matrix, Vector> > linearIterationStep, int nIterationSteps = 1)
+{
+ return [=](const Matrix& A, Vector& x, const Vector& b, const Dune::Solvers::DefaultBitVector_t<Vector>& ignore) {
+ linearIterationStep->setIgnore(ignore);
+ linearIterationStep->setProblem(A, x, b);
+ linearIterationStep->preprocess();
+
+ for (int i = 0; i < nIterationSteps; ++i)
+ linearIterationStep->iterate();
+ };
+}
+
+#endif
diff --git a/src/spatial-solving/tnnmg/linearization.hh b/src/spatial-solving/tnnmg/linearization.hh
index 73b90b42..8dcecd15 100644
--- a/src/spatial-solving/tnnmg/linearization.hh
+++ b/src/spatial-solving/tnnmg/linearization.hh
@@ -144,13 +144,11 @@ class Linearization {
// -grad is needed for Newton step
negativeGradient_ *= -1;
-
-
// truncate gradient and hessian
- truncateVector(negativeGradient_, truncationFlags_);
- truncateMatrix(hessian_, truncationFlags_, truncationFlags_);
+ // truncateVector(negativeGradient_, truncationFlags_);
+ // truncateMatrix(hessian_, truncationFlags_, truncationFlags_);
- //print(hessian_, "Hessian:");
+ // the actual truncation is done by the solver of linearized problem
}
diff --git a/src/spatial-solving/tnnmg/tnnmgstep.hh b/src/spatial-solving/tnnmg/tnnmgstep.hh
index 6846b234..dc2ca30d 100644
--- a/src/spatial-solving/tnnmg/tnnmgstep.hh
+++ b/src/spatial-solving/tnnmg/tnnmgstep.hh
@@ -15,10 +15,11 @@
#include <dune/solvers/solvers/linearsolver.hh>
#include <dune/solvers/common/canignore.hh>
-#include "localproblem.hh"
+#include <dune/solvers/solvers/umfpacksolver.hh>
+#include "localproblem.hh"
-#include <dune/tnnmg/iterationsteps/linearcorrection.hh>
+#include "linearcorrection.hh"
#include <dune/tectonic/../../src/utils/debugutils.hh>
@@ -65,7 +66,7 @@ class TNNMGStep :
: Base(x),
f_(&f),
nonlinearSmoother_(nonlinearSmoother),
- //linearCorrection_(makeLinearCorrection<ConstrainedMatrix>(iterativeSolver)),
+ linearCorrection_(makeLinearCorrection<ConstrainedMatrix>(iterativeSolver)),
projection_(projection),
lineSolver_(lineSolver)
{}
@@ -85,7 +86,7 @@ class TNNMGStep :
: Base(x),
f_(&f),
nonlinearSmoother_(nonlinearSmoother),
- linearSolver_(linearSolver),
+ linearCorrection_(makeLinearCorrection(linearSolver)),
projection_(projection),
lineSolver_(lineSolver)
{}
@@ -106,7 +107,7 @@ class TNNMGStep :
: Base(x),
f_(&f),
nonlinearSmoother_(nonlinearSmoother),
- //linearCorrection_(makeLinearCorrection(linearIterationStep, noOfLinearIterationSteps)),
+ linearCorrection_(makeLinearCorrection(linearIterationStep, noOfLinearIterationSteps)),
projection_(projection),
lineSolver_(lineSolver)
{}
@@ -165,21 +166,21 @@ class TNNMGStep :
// std::cout << "- computing linear correction..." << std::endl;
- /* auto canIgnoreCast = std::dynamic_pointer_cast<CanIgnore<Dune::Solvers::DefaultBitVector_t<Vector>>>( linearSolver_ );
- if (canIgnoreCast)
- canIgnoreCast->setIgnore(*this->ignoreNodes_);*/
+ linearCorrection_(A, constrainedCorrection_, r, linearization_->truncated());
- LocalProblem<ConstrainedMatrix, ConstrainedVector> localProblem(A, r, *this->ignoreNodes_);
- Vector newR;
- localProblem.getLocalRhs(constrainedCorrection_, newR);
+ LocalProblem<ConstrainedMatrix, ConstrainedVector> localProblem(A, r, linearization_->truncated());
+ Vector newR, constrainedCorrection_compare;
+ Solvers::resizeInitializeZero(constrainedCorrection_compare, r);
+ localProblem.getLocalRhs(constrainedCorrection_compare, newR);
/*print(*this->ignoreNodes_, "ignoreNodes:");
print(A, "A:");
print(localProblem.getMat(), "localMat:");*/
+ auto directSolver = std::make_shared<Dune::Solvers::UMFPackSolver<ConstrainedMatrix, ConstrainedVector>>();
- linearSolver_->setProblem(localProblem.getMat(), constrainedCorrection_, newR);
- linearSolver_->preprocess();
- linearSolver_->solve();
+ directSolver->setProblem(localProblem.getMat(), constrainedCorrection_compare, newR);
+ directSolver->preprocess();
+ directSolver->solve();
// linearCorrection_(A, constrainedCorrection_, r);
@@ -249,7 +250,7 @@ class TNNMGStep :
std::shared_ptr<Linearization> linearization_;
//! \brief linear correction
- std::shared_ptr<LinearSolver> linearSolver_;
+ LinearCorrection<ConstrainedMatrix, ConstrainedVector> linearCorrection_;
typename Linearization::ConstrainedVector constrainedCorrection_;
Vector correction_;
--
GitLab