Add a semi-smooth primal dual contact solver

Based on the work of Alexander Popp et al.

dune/contact/solvers/primaldualsemismoothnewtonsolver.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set ts=4 sw=2 et sts=2:
+
+#include <limits>
+#include <dune/common/timer.hh>
+#include <dune/istl/io.hh>
+
+#include <dune/solvers/solvers/quadraticipopt.hh>
+#include <dune/solvers/common/boxconstraint.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/umfpack.hh>
+#include <dune/istl/solvers.hh>
+
+namespace Dune {
+namespace Contact {
+
+template<class ProblemType>
+void PrimalDualSemiSmoothNewtonSolver<ProblemType>::solve()
+{
+  auto& iterates = *iterates_;
+
+  size_t totalSize(0);
+  for (size_t i=0; i<iterates.size(); i++)
+    totalSize += iterates[i].size();
+
+  // If no Dirichlet conditions are prescribed we have to setup the initial correction
+  if (correction_.size() != dim*totalSize) {
+    correction_.resize(dim*totalSize);
+    correction_ = 0;
+  }
+
+  problem_->contactAssembler()->getContactCouplings()[0]->path_ = debugPath_;
+  problem_->contactAssembler()->getContactCouplings()[0]->it_ = -1;
+
+  // compute energy infeasibility of initial iterate and build the contact mapping
+  energy_ = problem_->energy(iterates);
+  infeasibility_ = problem_->infeasibility(iterates);
+
+  // if there are non-homogeneous Dirichlet conditions,
+  // then we have to adjust the initial infeasibility
+  infeasibility_ = std::max(infeasibility_, correction_.infinity_norm());
+
+  if (this->verbosity_ == Solver::REDUCED or this->verbosity_== Solver::FULL)
+    std::cout<< " Initial energy: " << energy_
+             << " and infeasibility: " << infeasibility_ << std::endl;
+
+  // initialise the iteration counter and error
+  int it(0);
+  field_type error = std::numeric_limits<field_type>::max();
+
+  // assemble linearisations of initial deformation
+  problem_->assembleLinearisations(iterates);
+  // initialise Lagrange multiplier with zero
+  lagrangeMultiplier_.resize(problem_->getConstraints().size());
+ // lagrangeMultiplier_ = -1000;
+  lagrangeMultiplier_ = problem_->leastSquaresLagrangeMultiplier();
+
+  problem_->initialiseEmptyActiveSet();
+
+  // create IpOpt solver
+  QuadraticIPOptSolver<ScalarMatrix, ScalarVector> ipoptSolver(localTol_, localMaxIterations_,
+                                                            Solver::QUIET, linearSolverType_);
+
+  bool activeSetChanged(true);
+  // Loop until desired tolerance or maximum number of iterations is reached
+  for (; activeSetChanged or (it<this->maxIterations_ and (error>this->tolerance_ or std::isnan(error))); it++)
+  {
+
+    if (this->verbosity_ == Solver::REDUCED or this->verbosity_== Solver::FULL)
+      std::cout << "************************\n"
+                << "*  PriDual iteration " << it << "\n"
+                << "************************\n";
+
+    problem_->contactAssembler()->getContactCouplings()[0]->path_ = debugPath_;
+    problem_->contactAssembler()->getContactCouplings()[0]->it_ = it;
+
+    // Measure norm of the old iterate
+    field_type oldNorm(0);
+    for (size_t i=0; i<iterates_->size(); i++)
+      oldNorm += (*errorNorms_[i])(iterates[i]);
+
+    // assemble the local quadratic problem
+    problem_->assemblePrimalDualProblem(lagrangeMultiplier_);
+
+    ScalarVector priDualIterate = combineVector(correction_, ScalarVector(lagrangeMultiplier_.size(),0));
+
+    field_type newEnergy;
+    field_type newInfeasibility;
+    std::vector<BlockVector> newIterates;
+
+    // setup the linearised problem
+    problem_->addHardDirichlet(correction_, dirichletNodes_);
+    ipoptSolver.setProblem(problem_->priDualMatrix(), priDualIterate, problem_->priDualRhs());
+
+    Timer timer;
+
+#if HAVE_SUITESPARSE_UMFPACK
+    UMFPack<ScalarMatrix> umfPack(problem_->priDualMatrix(), 0);
+    std::cout<<"took "<<timer.elapsed()<<" to decompose matrix\n";
+
+    // solve
+    InverseOperatorResult res;
+    auto copyRhs = problem_->priDualRhs();
+    umfPack.apply(priDualIterate, copyRhs, res);
+    if (!res.converged)
+      std::cout<<"Warning: Not converged! Reduction "<<res.reduction<<" iterations "<<res.iterations<<std::endl;
+#endif
+
+    auto residual = problem_->priDualRhs();
+    problem_->priDualMatrix().mmv(priDualIterate, residual);
+    std::cout<<"Residual " << residual.infinity_norm()<<std::endl;
+
+
+
+    if (this->verbosity_ == Solver::FULL)
+      std::cout << "Solved local problem in " << timer.stop() << " secs\n";
+
+    std::tie(correction_, lagrangeMultiplier_) = decomposeVector(priDualIterate);
+
+    // new iterate
+    int counter(0);
+    newIterates = iterates;
+    for (size_t i=0; i<newIterates.size(); i++)
+      for (size_t j=0; j<newIterates[i].size(); j++)
+        for (int k = 0; k < dim; k++)
+          newIterates[i][j][k] += correction_[counter++][0];
+
+    std::cout<<"iterates inf "<<newIterates[0].infinity_norm() << " and " << newIterates[1].infinity_norm()<<std::endl;
+    // compute new energy and infeasibility
+    newEnergy = problem_->energy(newIterates);
+    newInfeasibility = problem_->infeasibility(newIterates);
+
+    // update the active set, needs to be called after the deformation is updated
+    // which is done in "infeasibility"
+    problem_->updateActiveSetAndLagrangeMultiplier(lagrangeMultiplier_);
+    activeSetChanged = (problem_->changedActiveNodes() > 0);
+
+    if (this->verbosity_ == Solver::FULL) {
+      std::cout << "Active nodes changed " << problem_->changedActiveNodes()
+                << " active now " << problem_->getActiveSet().count() << " \n";
+      std::cout << "New energy of potential iterate " << newEnergy
+                << " and infeasibility " << newInfeasibility << std::endl;
+    }
+
+    // update the Newton problem
+    problem_->assembleLinearisations(newIterates);
+
+    // if the correction is very small stop the loop
+    field_type infNorm = correction_.infinity_norm();
+    if(this->verbosity_==Solver::FULL)
+      std::cout << "Correction infinite norm: " << infNorm << std::endl;
+
+    // compute error estimate
+    std::vector<BlockVector> corrections(iterates.size());
+    counter = 0;
+    for (size_t i=0; i<iterates.size();i++) {
+      corrections[i].resize(iterates[i].size());
+      for (size_t j=0; j<corrections[i].size(); j++)
+        for (int k = 0; k < dim; k++)
+          corrections[i][j][k] = correction_[counter++][0];
+    }
+
+    auto residuals = problem_->residuals(lagrangeMultiplier_);
+    for (int k=0; k<2; k++)
+      for (size_t j=0; j<residuals[k].size();j++)
+        if (dN_[k][j].any())
+          residuals[k][j] = 0;
+
+    field_type absError = 0;
+    for (size_t i=0; i < corrections.size(); ++i)
+      absError += (*errorNorms_[i])(residuals[i]);
+    absError += newInfeasibility;
+
+    error = absError;
+    if (this->verbosity_ == Solver::REDUCED or this->verbosity_== Solver::FULL)
+      std::cout << "In filter iteration "<< it
+                << ", absolut norm of correction : " << absError
+                << " relative error " << error
+                << "\nEnergy " << newEnergy
+                << " and infeasibility " << newInfeasibility << std::endl;
+
+    // update iterates, energy and infeasibility
+    iterates = newIterates;
+    energy_ = newEnergy;
+    infeasibility_ = newInfeasibility;
+
+    // set correction to zero to remove any Dirichlet values
+    correction_  = 0;
+
+  } // end of filter loop
+
+  if (this->verbosity_ == Solver::REDUCED or this->verbosity_== Solver::FULL)
+    std::cout << "Solved problem in " << it << " iterations with final energy " << energy_
+              << " and infeasibility " << infeasibility_ << std::endl;
+}
+
+} /* namespace Contact */
+} /* namespace Dune */

dune/contact/solvers/primaldualsemismoothnewtonsolver.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set ts=4 sw=2 et sts=2:
+#ifndef DUNE_CONTACT_SOLVERS_PRIMAL_DUAL_SEMI_SMOOTH_NEWTON_SOLVER_HH
+#define DUNE_CONTACT_SOLVERS_PRIMAL_DUAL_SEMI_SMOOTH_NEWTON_SOLVER_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include <dune/solvers/norms/norm.hh>
+#include <dune/solvers/solvers/iterativesolver.hh>
+
+namespace Dune {
+namespace Contact {
+
+/** \brief Primal-Dual Semi-Smooth Newton method
+ *
+ *  See Popp, Wall et al. 2011
+*/
+template <class Problem>
+class PrimalDualSemiSmoothNewtonSolver
+    : public Dune::Solvers::IterativeSolver<typename Problem::BlockVector> {
+
+  using BlockVector = typename Problem::BlockVector;
+  using field_type = typename Dune::FieldTraits<BlockVector>::field_type;
+  using ScalarVector = typename Problem::ScalarVector;
+  using ScalarMatrix = typename Problem::ScalarMatrix;
+  using Matrix = typename Problem::MatrixType;
+  using JacobianType = typename Problem::JacobianType;
+
+  static constexpr int dim = BlockVector::block_type::dimension;
+  using Base = Dune::Solvers::IterativeSolver<BlockVector>;
+
+public:
+  PrimalDualSemiSmoothNewtonSolver(
+      const Dune::ParameterTree& config, std::vector<BlockVector>& x,
+      Problem& problem,
+      const std::vector<std::shared_ptr<Norm<BlockVector> > >& errorNorms)
+      : Base(config.get<field_type>("tolerance"),
+             config.get<int>("maxIterations"),
+             config.get("verbosity", NumProc::FULL)) {
+    setupParameter(config);
+    setProblem(x, problem);
+    setErrorNorm(errorNorms);
+
+    size_t totalSize(0);
+    for (size_t i = 0; i < x.size(); i++)
+      totalSize += x[i].size();
+
+    dirichletNodes_.resize(totalSize*dim, false);
+  }
+
+  //! Setup the trust region parameter from a config file
+  void setupParameter(const Dune::ParameterTree& config) {
+    localTol_ = config.get<field_type>("localTolerance");
+    localMaxIterations_ = config.get<int>("localMaxIterations");
+    linearSolverType_ = config.get("linearSolverType", "");
+  }
+
+  /** \brief Set the Dirichlet values.
+     *
+     * The Dirichlet conditions are enforced during the solution
+     * of the first linearised problem rather than modifying the
+     * initial iterate of the filter scheme. This avoids possible problems
+     *in the energy evaluation due to inverted elements.
+     * Hence Dirichlet values for the correction must be handed over!
+     */
+  void setDirichletValues(
+      const std::vector<BlockVector>& dirichletValues,
+      const std::vector<Dune::BitSetVector<dim> >& dirichletNodes) {
+    // only adjust correction Dirichlet conditions are prescribed
+    correction_.resize(dirichletNodes_.size());
+    correction_ = 0;
+    int counter(0);
+    for (size_t i = 0; i < dirichletNodes.size(); i++)
+      for (size_t j = 0; j < dirichletNodes[i].size(); counter++, j++)
+        if (dirichletNodes[i][j].any())
+          for (int k = 0; k < dim; k++) {
+            dirichletNodes_[dim*counter + k] = dirichletNodes[i][j][k];
+            correction_[dim*counter + k][0] = dirichletValues[i][j][k];
+          }
+
+    dN_ = dirichletNodes;
+  }
+
+  //! Set initial iterate and the problem class
+  void setProblem(std::vector<BlockVector>& x, Problem& problem) {
+    iterates_ = Dune::stackobject_to_shared_ptr<std::vector<BlockVector> >(x);
+    problem_ = Dune::stackobject_to_shared_ptr<Problem>(problem);
+  }
+
+  //! Set the error norms to be used within the filter method
+  void
+  setErrorNorm(const std::vector<std::shared_ptr<Norm<BlockVector> > >& errorNorms) {
+    errorNorms_ = errorNorms;
+  }
+
+  //! Solve the problem
+  void solve();
+
+  //! Return solution object
+  const std::vector<BlockVector>& getSol() const { return *iterates_; }
+
+  std::string debugPath_;
+  std::vector<Dune::BitSetVector<dim> > dN_;
+
+private:
+  ScalarVector combineVector(const ScalarVector& primal, const ScalarVector& dual) const {
+    size_t index(0);
+    ScalarVector combinedVector(primal.size() + dual.size());
+    for (size_t i = 0; i < primal.size(); i++)
+      combinedVector[index++] = primal[i];
+
+    for (size_t i = 0; i < dual.size(); i++)
+      combinedVector[index++] = dual[i];
+    return combinedVector;
+  }
+
+  std::tuple<ScalarVector, ScalarVector> decomposeVector(const ScalarVector& primalDual) const {
+    size_t index(0);
+    ScalarVector primal(correction_.size());
+    for (size_t i = 0; i < primal.size(); i++)
+        primal[i] = primalDual[index++];
+
+    ScalarVector dual(lagrangeMultiplier_.size());
+    for (size_t i = 0; i < dual.size(); i++)
+      dual[i] = primalDual[index++];
+    return std::tie(primal, dual);
+  }
+
+  //! The iterates of the global scheme
+  std::shared_ptr<std::vector<BlockVector> > iterates_;
+
+  ScalarVector lagrangeMultiplier_;
+
+  //! Dirichlet degrees of freedom
+  Dune::BitSetVector<1> dirichletNodes_;
+
+  //! The iterate for the local problems, containing possible Dirichlet values
+  ScalarVector correction_;
+
+  //! The problem type, providing methods to assemble the nonlinear and
+  //! linearised
+  //! energy and constraints
+  std::shared_ptr<Problem> problem_;
+
+  //! The norms to estimate the errors
+  std::vector<std::shared_ptr<Norm<BlockVector> > > errorNorms_;
+
+  //! Tolerance for the solution of the local problem
+  field_type localTol_;
+  //! Max iterations for the solution of the local problem
+  int localMaxIterations_;
+
+  //! Store the old energy
+  field_type energy_;
+  //! Store the infeasibility of the old iterate
+  field_type infeasibility_;
+
+  //! The linear solver to be used within IpOpt
+  std::string linearSolverType_;
+};
+
+} /* namespace Contact */
+} /* namespace Dune */
+
+#include "primaldualsemismoothnewtonsolver.cc"
+
+#endif