diff --git a/dune/tectonic/time-stepping/step.hh b/dune/tectonic/time-stepping/step.hh index a6433c33e6321f631ea389736dd1bbdc89c45b75..86d4b5e24c86646b3402f5b1f6c384772ef3eeaf 100644 --- a/dune/tectonic/time-stepping/step.hh +++ b/dune/tectonic/time-stepping/step.hh @@ -12,7 +12,7 @@ #include <dune/solvers/iterationsteps/cgstep.hh> #include <dune/solvers/solvers/loopsolver.hh> -#include "../spatial-solving/preconditioners/multilevelpatchpreconditioner.hh" +#include "../spatial-solving/makelinearsolver.hh" #include <dune/tectonic/utils/reductionfactors.hh> @@ -130,71 +130,7 @@ class Step : protected StepBase<Factory, ContactNetwork, Updaters, ErrorNorms> { using Vector = typename Factory::Vector; using Matrix = typename Factory::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); - */ - - // set multigrid solver - auto smoother = TruncatedBlockGSStep<Matrix, Vector>(); - - // 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(this->contactNetwork_.nLevels()-1); - for (size_t i=0; i<transfer.size(); i++) { - transfer[i] = std::make_shared<TransferOperator>(); - transfer[i]->setup(this->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(smoother); - linearMultigridStep->setTransferOperators(transfer); - - Norm norm(*linearMultigridStep); - auto linearSolver = std::make_shared<LinearSolver>(linearMultigridStep, this->parset_.template get<int>("solver.tnnmg.main.multi"), this->parset_.template get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"), norm, Solver::QUIET); + auto linearSolver = makeLinearSolver<ContactNetwork, Vector>(this->parset_, this->contactNetwork_); Vector x; x.resize(nBodyAssembler_.totalHasObstacle_.size());