From c263d7e815484190fc01120b1a666861ef525363 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Mon, 27 Apr 2009 16:16:32 +0000
Subject: [PATCH] new method to setup a solver using a truncated cg as the
inner solver. This compiles, but it doesn't work yet
[[Imported from SVN: r4105]]
---
src/riemanniantrsolver.cc | 102 ++++++++++++++++++++++++++++++++------
src/riemanniantrsolver.hh | 13 +++++
2 files changed, 100 insertions(+), 15 deletions(-)
diff --git a/src/riemanniantrsolver.cc b/src/riemanniantrsolver.cc
index 461f3f15..bc833fe8 100644
--- a/src/riemanniantrsolver.cc
+++ b/src/riemanniantrsolver.cc
@@ -6,16 +6,20 @@
#include <dune/disc/miscoperators/laplace.hh>
#include <dune/disc/operators/p1operator.hh>
+// For using a monotone multigrid as the inner solver
#include <dune-solvers/iterationsteps/trustregiongsstep.hh>
#include <dune-solvers/iterationsteps/mmgstep.hh>
#include <dune-solvers/transferoperators/truncatedcompressedmgtransfer.hh>
#include <dune-solvers/transferoperators/mandelobsrestrictor.hh>
#include <dune-solvers/solvers/iterativesolver.hh>
+#include "maxnormtrustregion.hh"
+
+// For using a truncated cg as the inner solver
+#include <dune-solvers/solvers/tcgsolver.hh>
#include <dune-solvers/norms/energynorm.hh>
#include <dune-solvers/norms/h1seminorm.hh>
-#include "maxnormtrustregion.hh"
// for debugging
#include "../test/fdcheck.hh"
@@ -23,20 +27,20 @@
template <class GridType, class TargetSpace>
void RiemannianTrustRegionSolver<GridType,TargetSpace>::
setup(const GridType& grid,
- const GeodesicFEAssembler<typename GridType::LeafGridView,TargetSpace>* assembler,
- const SolutionType& x,
- const Dune::BitSetVector<blocksize>& dirichletNodes,
- double tolerance,
- int maxTrustRegionSteps,
- double initialTrustRegionRadius,
- int multigridIterations,
- double mgTolerance,
- int mu,
- int nu1,
- int nu2,
- int baseIterations,
- double baseTolerance,
- bool instrumented)
+ const GeodesicFEAssembler<typename GridType::LeafGridView,TargetSpace>* assembler,
+ const SolutionType& x,
+ const Dune::BitSetVector<blocksize>& dirichletNodes,
+ double tolerance,
+ int maxTrustRegionSteps,
+ double initialTrustRegionRadius,
+ int multigridIterations,
+ double mgTolerance,
+ int mu,
+ int nu1,
+ int nu2,
+ int baseIterations,
+ double baseTolerance,
+ bool instrumented)
{
grid_ = &grid;
assembler_ = assembler;
@@ -139,6 +143,74 @@ setup(const GridType& grid,
}
+
+template <class GridType, class TargetSpace>
+void RiemannianTrustRegionSolver<GridType,TargetSpace>::
+setupTCG(const GridType& grid,
+ const GeodesicFEAssembler<typename GridType::LeafGridView,TargetSpace>* assembler,
+ const SolutionType& x,
+ const Dune::BitSetVector<blocksize>& dirichletNodes,
+ double tolerance,
+ int maxTrustRegionSteps,
+ double initialTrustRegionRadius,
+ int innerIterations,
+ double innerTolerance,
+ bool instrumented)
+{
+ grid_ = &grid;
+ assembler_ = assembler;
+ x_ = x;
+ this->tolerance_ = tolerance;
+ maxTrustRegionSteps_ = maxTrustRegionSteps;
+ initialTrustRegionRadius_ = initialTrustRegionRadius;
+ innerIterations_ = innerIterations;
+ innerTolerance_ = innerTolerance;
+ instrumented_ = instrumented;
+
+ // ////////////////////////////////////////////////////////////
+ // Create Hessian matrix and its occupation structure
+ // ////////////////////////////////////////////////////////////
+
+ hessianMatrix_ = std::auto_ptr<MatrixType>(new MatrixType);
+ Dune::MatrixIndexSet indices(grid_->size(1), grid_->size(1));
+ assembler_->getNeighborsPerVertex(indices);
+ indices.exportIdx(*hessianMatrix_);
+
+ // //////////////////////////////////////////////////////////////////////////////////////
+ // Assemble a Laplace matrix to create a norm that's equivalent to the H1-norm
+ // //////////////////////////////////////////////////////////////////////////////////////
+ Dune::LeafP1Function<GridType,double> u(grid),f(grid);
+ Dune::LaplaceLocalStiffness<typename GridType::LeafGridView,double> laplaceStiffness;
+ Dune::LeafP1OperatorAssembler<GridType,double,1>* A = new Dune::LeafP1OperatorAssembler<GridType,double,1>(grid);
+ A->assemble(laplaceStiffness,u,f);
+
+ typedef typename Dune::LeafP1OperatorAssembler<GridType,double,1>::RepresentationType LaplaceMatrixType;
+
+ if (h1SemiNorm_)
+ delete h1SemiNorm_;
+
+ h1SemiNorm_ = new H1SemiNorm<CorrectionType>(**A);
+
+ // ////////////////////////////////////////////////////
+ // Create a truncated conjugate gradient solver
+ // ////////////////////////////////////////////////////
+
+ innerSolver_ = new TruncatedCGSolver<MatrixType,CorrectionType>(*hessianMatrix_,
+ x,
+ this->rhs_,
+ innerIterations_,
+ innerTolerance_,
+ h1SemiNorm_,
+ initialTrustRegionRadius,
+ Solver::QUIET);
+
+ // Write all intermediate solutions, if requested
+ if (instrumented_
+ && dynamic_cast<IterativeSolver<CorrectionType>*>(innerSolver_))
+ dynamic_cast<IterativeSolver<CorrectionType>*>(innerSolver_)->historyBuffer_ = "tmp/mgHistory";
+
+}
+
template <class GridType, class TargetSpace>
void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
{
diff --git a/src/riemanniantrsolver.hh b/src/riemanniantrsolver.hh
index 67decc48..2f6e7528 100644
--- a/src/riemanniantrsolver.hh
+++ b/src/riemanniantrsolver.hh
@@ -41,6 +41,7 @@ public:
hessianMatrix_(std::auto_ptr<MatrixType>(NULL)), h1SemiNorm_(NULL)
{}
+ /** \brief Set up the solver using a monotone multigrid method as the inner solver */
void setup(const GridType& grid,
const GeodesicFEAssembler<typename GridType::LeafGridView, TargetSpace>* rodAssembler,
const SolutionType& x,
@@ -57,6 +58,18 @@ public:
double baseTolerance,
bool instrumented);
+ /** \brief Set up the solver using a truncated cg method as the inner solver */
+ void setupTCG(const GridType& grid,
+ const GeodesicFEAssembler<typename GridType::LeafGridView, TargetSpace>* rodAssembler,
+ const SolutionType& x,
+ const Dune::BitSetVector<blocksize>& dirichletNodes,
+ double tolerance,
+ int maxTrustRegionSteps,
+ double initialTrustRegionRadius,
+ int innerIterations,
+ double innerTolerance,
+ bool instrumented);
+
void solve();
void setInitialSolution(const SolutionType& x) {
--
GitLab