From 16065f32f39dbc400b8e39a38aaa940c7fe7c739 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Mon, 12 Mar 2012 15:07:03 +0100
Subject: [PATCH] Add and Externalise MySolver class (mostly)
---
src/mysolver.cc | 65 +++++++++++++++++++++++++++++++++++++++
src/mysolver.hh | 44 +++++++++++++++++++++++++++
src/mysolver_tmpl.cc | 38 +++++++++++++++++++++++
src/one-body-sample.cc | 69 ++++++++----------------------------------
4 files changed, 159 insertions(+), 57 deletions(-)
create mode 100644 src/mysolver.cc
create mode 100644 src/mysolver.hh
create mode 100644 src/mysolver_tmpl.cc
diff --git a/src/mysolver.cc b/src/mysolver.cc
new file mode 100644
index 00000000..7606e56d
--- /dev/null
+++ b/src/mysolver.cc
@@ -0,0 +1,65 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#ifdef HAVE_IPOPT
+#undef HAVE_IPOPT
+#endif
+
+#include <dune/solvers/common/numproc.hh> // Solver::FULL
+
+#include "mysolver.hh"
+
+template <int dim, class VectorType, class OperatorType, class GridType,
+ class BlockProblemType>
+MySolver<dim, VectorType, OperatorType, GridType, BlockProblemType>::MySolver(
+ Dune::ParameterTree parset, int refinements, double solver_tolerance,
+ GridType const &grid, Dune::BitSetVector<dim> const &ignoreNodes)
+ : baseEnergyNorm(linearBaseSolverStep),
+ linearBaseSolver(&linearBaseSolverStep,
+ parset.get<int>("linear.maxiterations"),
+ solver_tolerance, &baseEnergyNorm, Solver::QUIET),
+ transferOperators(refinements),
+ multigridStep(new SolverType(linearIterationStep, nonlinearSmoother)) {
+ // linear iteration step
+ linearIterationStep.setNumberOfLevels(refinements + 1);
+ linearIterationStep.setMGType(parset.get<int>("linear.mu"),
+ parset.get<int>("linear.nu1"),
+ parset.get<int>("linear.nu2"));
+ linearIterationStep.basesolver_ = &linearBaseSolver;
+ linearIterationStep.setSmoother(&linearPresmoother, &linearPostsmoother);
+
+ // transfer operators
+ for (auto &x : transferOperators)
+ x = new CompressedMultigridTransfer<VectorType>;
+ TransferOperatorAssembler<GridType>(grid)
+ .assembleOperatorPointerHierarchy(transferOperators);
+ linearIterationStep.setTransferOperators(transferOperators);
+
+ // tnnmg iteration step
+ multigridStep->setSmoothingSteps(parset.get<int>("main.nu1"),
+ parset.get<int>("main.mu"),
+ parset.get<int>("main.nu2"));
+ multigridStep->ignoreNodes_ = &ignoreNodes;
+}
+
+template <int dim, class VectorType, class OperatorType, class GridType,
+ class BlockProblemType>
+MySolver<dim, VectorType, OperatorType, GridType,
+ BlockProblemType>::~MySolver() {
+ for (auto &x : transferOperators)
+ delete x;
+
+ delete multigridStep;
+}
+
+template <int dim, class VectorType, class OperatorType, class GridType,
+ class BlockProblemType>
+typename MySolver<dim, VectorType, OperatorType, GridType,
+ BlockProblemType>::SolverType *
+MySolver<dim, VectorType, OperatorType, GridType,
+ BlockProblemType>::getSolver() {
+ return multigridStep;
+}
+
+//#include "mysolver_tmpl.cc"
diff --git a/src/mysolver.hh b/src/mysolver.hh
new file mode 100644
index 00000000..d80b755b
--- /dev/null
+++ b/src/mysolver.hh
@@ -0,0 +1,44 @@
+#ifndef MYSOLVER_HH
+#define MYSOLVER_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
+#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
+#include <dune/fufem/assemblers/transferoperatorassembler.hh>
+#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
+
+template <int dim, class VectorType, class OperatorType, class GridType,
+ class BlockProblemType>
+class MySolver {
+private:
+ typedef GenericNonlinearGS<BlockProblemType> NonlinearSmootherType;
+ typedef TruncatedNonsmoothNewtonMultigrid<BlockProblemType,
+ NonlinearSmootherType> SolverType;
+
+public:
+ MySolver(Dune::ParameterTree parset, int refinements, double solver_tolerance,
+ GridType const &grid, Dune::BitSetVector<dim> const &ignoreNodes);
+
+ ~MySolver();
+
+ SolverType *getSolver();
+
+private:
+ TruncatedBlockGSStep<OperatorType, VectorType> linearBaseSolverStep;
+ EnergyNorm<OperatorType, VectorType> baseEnergyNorm;
+ LoopSolver<VectorType> linearBaseSolver;
+ TruncatedBlockGSStep<OperatorType, VectorType> linearPresmoother;
+ TruncatedBlockGSStep<OperatorType, VectorType> linearPostsmoother;
+ MultigridStep<OperatorType, VectorType> linearIterationStep;
+ std::vector<CompressedMultigridTransfer<VectorType> *> transferOperators;
+ NonlinearSmootherType nonlinearSmoother;
+ SolverType *multigridStep;
+};
+
+#endif
diff --git a/src/mysolver_tmpl.cc b/src/mysolver_tmpl.cc
new file mode 100644
index 00000000..188501d4
--- /dev/null
+++ b/src/mysolver_tmpl.cc
@@ -0,0 +1,38 @@
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/tectonic/myblockproblem.hh>
+#include <dune/tectonic/myconvexproblem.hh>
+
+// {{{ 2D
+typedef Dune::FieldVector<double, 2> SmallVector2;
+typedef Dune::FieldMatrix<double, 2, 2> SmallMatrix2;
+typedef Dune::BlockVector<SmallVector2> VectorType2;
+typedef Dune::BCRSMatrix<SmallMatrix2> OperatorType2;
+
+typedef MyConvexProblem<OperatorType2, VectorType2> MyConvexProblemType2;
+typedef MyBlockProblem<MyConvexProblemType2> MyBlockProblemType2;
+
+typedef Dune::YaspGrid<2> GridType2;
+
+template class MySolver<2, VectorType2, OperatorType2, GridType2,
+ MyBlockProblemType2>;
+// }}}
+
+// {{{ 3D
+typedef Dune::FieldVector<double, 3> SmallVector3;
+typedef Dune::FieldMatrix<double, 3, 3> SmallMatrix3;
+typedef Dune::BlockVector<SmallVector3> VectorType3;
+typedef Dune::BCRSMatrix<SmallMatrix3> OperatorType3;
+
+typedef MyConvexProblem<OperatorType3, VectorType3> MyConvexProblemType3;
+typedef MyBlockProblem<MyConvexProblemType3> MyBlockProblemType3;
+
+typedef Dune::YaspGrid<3> GridType3;
+
+template class MySolver<3, VectorType3, OperatorType3, GridType3,
+ MyBlockProblemType3>;
+// }}}
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index ca674939..952e247a 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -6,10 +6,6 @@
#error srcdir unset
#endif
-#ifdef HAVE_IPOPT
-#undef HAVE_IPOPT
-#endif
-
#ifndef HAVE_PYTHON
#error Python is required
#endif
@@ -42,20 +38,14 @@
#include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
#include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
#include <dune/fufem/assemblers/operatorassembler.hh>
-#include <dune/fufem/assemblers/transferoperatorassembler.hh>
#include <dune/fufem/dunepython.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/fufem/functionspacebases/p0basis.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
#include <dune/fufem/sharedpointermap.hh>
#include <dune/solvers/common/numproc.hh> // Solver::FULL
-#include <dune/solvers/iterationsteps/multigridstep.hh>
-#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
-#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
-#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
#include <dune/tectonic/myblockproblem.hh>
#include <dune/tectonic/myconvexproblem.hh>
@@ -65,6 +55,8 @@
#include "print.hh"
#include "vtk.hh"
+#include "mysolver.cc" // TODO
+
int const dim = 2;
template <class GridView>
@@ -139,7 +131,6 @@ int main(int argc, char *argv[]) {
auto const solver_tolerance = parset.get<double>("solver.tolerance");
auto const refinements = parset.get<size_t>("grid.refinements");
- size_t const levels = refinements + 1;
auto const verbose = parset.get<bool>("verbose");
Solver::VerbosityMode const verbosity =
@@ -218,45 +209,11 @@ int main(int argc, char *argv[]) {
auto const nodalIntegrals =
assemble_frictional<GridType, GridView, SmallVector, P1Basis>(
leafView, p1Basis, frictionalNodes);
+
// {{{ Set up TNNMG solver
- // linear iteration step components
- TruncatedBlockGSStep<OperatorType, VectorType> linearBaseSolverStep;
- EnergyNorm<OperatorType, VectorType> baseEnergyNorm(linearBaseSolverStep);
- LoopSolver<VectorType> linearBaseSolver(
- &linearBaseSolverStep,
- parset.get<int>("solver.tnnmg.linear.maxiterations"), solver_tolerance,
- &baseEnergyNorm, Solver::QUIET);
- TruncatedBlockGSStep<OperatorType, VectorType> linearPresmoother;
- TruncatedBlockGSStep<OperatorType, VectorType> linearPostsmoother;
-
- // linear iteration step
- MultigridStep<OperatorType, VectorType> linearIterationStep;
- linearIterationStep.setNumberOfLevels(levels);
- linearIterationStep.setMGType(parset.get<int>("solver.tnnmg.linear.mu"),
- parset.get<int>("solver.tnnmg.linear.nu1"),
- parset.get<int>("solver.tnnmg.linear.nu2"));
- linearIterationStep.basesolver_ = &linearBaseSolver;
- linearIterationStep.setSmoother(&linearPresmoother, &linearPostsmoother);
-
- // transfer operators
- std::vector<CompressedMultigridTransfer<VectorType> *> transferOperators(
- refinements);
- for (auto &x : transferOperators)
- x = new CompressedMultigridTransfer<VectorType>;
- TransferOperatorAssembler<GridType>(*grid)
- .assembleOperatorPointerHierarchy(transferOperators);
- linearIterationStep.setTransferOperators(transferOperators);
-
- // tnnmg iteration step
- typedef GenericNonlinearGS<MyBlockProblemType> NonlinearSmootherType;
- NonlinearSmootherType nonlinearSmoother;
- TruncatedNonsmoothNewtonMultigrid<MyBlockProblemType, NonlinearSmootherType>
- multigridStep(linearIterationStep, nonlinearSmoother);
- multigridStep.setSmoothingSteps(parset.get<int>("solver.tnnmg.main.nu1"),
- parset.get<int>("solver.tnnmg.main.mu"),
- parset.get<int>("solver.tnnmg.main.nu2"));
- multigridStep.ignoreNodes_ = &ignoreNodes;
- // }}}
+ MySolver<dim, VectorType, OperatorType, GridType, MyBlockProblemType>
+ mySolver(parset.sub("solver.tnnmg"), refinements, solver_tolerance, *grid,
+ ignoreNodes);
std::fstream octave_writer("data", std::fstream::out);
timer.reset();
@@ -286,10 +243,11 @@ int main(int argc, char *argv[]) {
*myGlobalNonlinearity, b4);
MyBlockProblemType myBlockProblem(parset, myConvexProblem);
- multigridStep.setProblem(u4_diff, myBlockProblem);
+ auto multigridStep = mySolver.getSolver();
+ multigridStep->setProblem(u4_diff, myBlockProblem);
LoopSolver<VectorType> overallSolver(
- &multigridStep, parset.get<size_t>("solver.tnnmg.maxiterations"),
+ multigridStep, parset.get<size_t>("solver.tnnmg.maxiterations"),
solver_tolerance, &energyNorm, verbosity);
overallSolver.solve();
@@ -330,10 +288,11 @@ int main(int argc, char *argv[]) {
*myGlobalNonlinearity, b5);
MyBlockProblemType myBlockProblem(parset, myConvexProblem);
- multigridStep.setProblem(u5_diff, myBlockProblem);
+ auto multigridStep = mySolver.getSolver();
+ multigridStep->setProblem(u5_diff, myBlockProblem);
LoopSolver<VectorType> overallSolver(
- &multigridStep, parset.get<size_t>("solver.tnnmg.maxiterations"),
+ multigridStep, parset.get<size_t>("solver.tnnmg.maxiterations"),
solver_tolerance, &energyNorm, verbosity);
overallSolver.solve();
@@ -376,10 +335,6 @@ int main(int argc, char *argv[]) {
std::cout << "Making " << timesteps << " time steps took "
<< timer.elapsed() << "s" << std::endl;
- // Clean up after the TNNMG solver
- for (auto &x : transferOperators)
- delete x;
-
octave_writer.close();
if (parset.get<bool>("printFrictionalBoundary")) {
--
GitLab