Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
#ifndef DUNE_TECTONIC_SPATIAL_SOLVING_MAKELINEARSOLVER_HH
#define DUNE_TECTONIC_SPATIAL_SOLVING_MAKELINEARSOLVER_HH
#include <dune/common/parametertree.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/iterationsteps/multigridstep.hh>
#include <dune/solvers/iterationsteps/cgstep.hh>
#include <dune/solvers/solvers/loopsolver.hh>
#include "preconditioners/multilevelpatchpreconditioner.hh"
template <class ContactNetwork, class VectorType>
auto makeLinearSolver(const Dune::ParameterTree& parset, const ContactNetwork& contactNetwork) {
// make linear solver for linear correction in TNNMGStep
using Vector = VectorType;
using Matrix = typename ContactNetwork::Matrix;
/* old, pre multi threading, was unused
using Norm = EnergyNorm<Matrix, Vector>;
using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, Matrix, Vector>;
using LinearSolver = typename Dune::Solvers::LoopSolver<Vector>;
const auto& preconditionerParset = parset_.sub("solver.tnnmg.linear.preconditioner");
Dune::BitSetVector<1> activeLevels(contactNetwork_.nLevels(), true);
Preconditioner preconditioner(preconditionerParset, contactNetwork_, activeLevels);
preconditioner.setPatchDepth(preconditionerParset.template get<size_t>("patchDepth"));
preconditioner.build();
auto cgStep = std::make_shared<Dune::Solvers::CGStep<Matrix, Vector>>();
cgStep->setPreconditioner(preconditioner);
Norm norm(*cgStep);
return std::make_shared<LinearSolver>(cgStep, parset_.template get<int>("solver.tnnmg.main.multi"), parset_.template get<double>("solver.tnnmg.linear.tolerance"), norm, Solver::QUIET);
*/
// patch preconditioner only needs to be computed once per advance()
// make linear solver for linear correction in TNNMGStep
using Norm = EnergyNorm<Matrix, Vector>;
using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, Matrix, Vector>;
using LinearSolver = typename Dune::Solvers::LoopSolver<Vector>;
/*const auto& preconditionerParset = parset_.sub("solver.tnnmg.preconditioner");
Dune::BitSetVector<1> activeLevels(contactNetwork_.nLevels(), true);
Preconditioner preconditioner(preconditionerParset, contactNetwork_, activeLevels);
preconditioner.setPatchDepth(preconditionerParset.template get<size_t>("patchDepth"));
preconditioner.build();
auto cgStep = std::make_shared<Dune::Solvers::CGStep<Matrix, Vector>>();
cgStep->setPreconditioner(preconditioner);
Norm norm(*cgStep);
auto linearSolver = std::make_shared<LinearSolver>(cgStep, parset_.template get<int>("solver.tnnmg.main.multi"), parset_.template get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"), norm, Solver::QUIET);
*/
// transfer operators need to be recomputed on change due to setDeformation()
using TransferOperator = NBodyContactTransfer<ContactNetwork, Vector>;
using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
TransferOperators transfer(contactNetwork.nLevels()-1);
for (size_t i=0; i<transfer.size(); i++) {
transfer[i] = std::make_shared<TransferOperator>();
transfer[i]->setup(contactNetwork, i, i+1);
}
// Remove any recompute filed so that initially the full transferoperator is assembled
for (size_t i=0; i<transfer.size(); i++)
std::dynamic_pointer_cast<TruncatedMGTransfer<Vector> >(transfer[i])->setRecomputeBitField(nullptr);
auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >();
linearMultigridStep->setMGType(1, 3, 3);
linearMultigridStep->setSmoother(TruncatedBlockGSStep<Matrix, Vector>());
linearMultigridStep->setTransferOperators(transfer);
return std::make_shared<LinearSolver>(linearMultigridStep, parset.template get<int>("solver.tnnmg.main.multi"), parset.template get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"), Norm(*linearMultigridStep), Solver::QUIET);