pathfollowingfunctions.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers mueller@ibb.uni-stuttgart.de
// SPDX-License-Identifier: LGPL-3.0-or-later
 
#pragma once
 
#include <cmath>
#include <optional>
#include <string>
#include <utility>
#include <vector>
 
#include <dune/common/exceptions.hh>
 
#include <Eigen/Core>
 
#include <ikarus/solver/linearsolver/linearsolver.hh>
#include <ikarus/utils/concepts.hh>
#include <ikarus/utils/defaultfunctions.hh>
#include <ikarus/utils/traits.hh>
 
namespace Ikarus {
struct SubsidiaryArgs
{
  double stepSize;              
  Eigen::VectorX<double> DD;    
  double Dlambda;               
  double f;                     
  Eigen::VectorX<double> dfdDD; 
  double dfdDlambda;            
  int currentStep;              
 
  void setZero(const Concepts::EigenType auto& firstParameter) {
    stepSize = 0.0;
    DD.resizeLike(firstParameter);
    DD.setZero();
    Dlambda = 0.0;
    f       = 0.0;
    dfdDD.resizeLike(firstParameter);
    dfdDD.setZero();
    dfdDlambda  = 0.0;
    currentStep = 0;
  }
};
 
struct ArcLength
{
  void operator()(SubsidiaryArgs& args) const {
    if (psi) {
      const auto root = sqrt(args.DD.squaredNorm() + psi.value() * psi.value() * args.Dlambda * args.Dlambda);
      args.f          = root - args.stepSize;
      args.dfdDD      = args.DD / root;
      args.dfdDlambda = (psi.value() * psi.value() * args.Dlambda) / root;
    } else {
      DUNE_THROW(Dune::InvalidStateException, "You have to call initialPrediction first. Otherwise psi is not defined");
    }
  }
 
  template <typename F>
  requires(Ikarus::Concepts::LinearSolverCheck<Ikarus::LinearSolver, typename F::Traits::template Range<1>,
                                               typename F::Domain::SolutionVectorType>)
  void initialPrediction(typename F::Domain& req, F& residual, SubsidiaryArgs& args) {
    SolverTypeTag solverTag;
    using JacobianType = std::remove_cvref_t<typename F::Traits::template Range<1>>;
    static_assert((traits::isSpecializationTypeAndNonTypes<Eigen::Matrix, JacobianType>::value) or
                      (traits::isSpecializationTypeNonTypeAndType<Eigen::SparseMatrix, JacobianType>::value),
                  "Linear solver not implemented for the chosen derivative type of the non-linear operator");
 
    if constexpr (traits::isSpecializationTypeAndNonTypes<Eigen::Matrix, JacobianType>::value)
      solverTag = SolverTypeTag::d_LDLT;
    else
      solverTag = SolverTypeTag::sd_SimplicialLDLT;
 
    req.parameter()  = 1.0;
    decltype(auto) R = residual(req);
    decltype(auto) K = derivative(residual)(req);
 
    auto linearSolver = LinearSolver(solverTag); // for the linear predictor step
    linearSolver.analyzePattern(K);
    linearSolver.factorize(K);
    linearSolver.solve(args.DD, -R);
 
    const auto DD2 = args.DD.squaredNorm();
 
    psi    = sqrt(DD2);
    auto s = sqrt(psi.value() * psi.value() + DD2);
 
    args.DD      = args.DD * args.stepSize / s;
    args.Dlambda = args.stepSize / s;
 
    req.globalSolution() = args.DD;
    req.parameter()      = args.Dlambda;
  }
 
  template <typename F>
  void intermediatePrediction(typename F::Domain& req, F& residual, SubsidiaryArgs& args) {
    req.globalSolution() += args.DD;
    req.parameter() += args.Dlambda;
  }
 
  constexpr auto name() const { return std::string("Arc length"); }
 
private:
  std::optional<double> psi;
};
 
struct LoadControlSubsidiaryFunction
{
  void operator()(SubsidiaryArgs& args) const {
    args.f = args.Dlambda - args.stepSize;
    args.dfdDD.setZero();
    args.dfdDlambda = 1.0;
  }
 
  template <typename F>
  void initialPrediction(typename F::Domain& req, F& residual, SubsidiaryArgs& args) {
    args.Dlambda    = args.stepSize;
    req.parameter() = args.Dlambda;
  }
 
  template <typename F>
  void intermediatePrediction(typename F::Domain& req, F& residual, SubsidiaryArgs& args) {
    req.parameter() += args.Dlambda;
  }
 
  constexpr auto name() const { return std::string("Load Control"); }
};
 
struct DisplacementControl
{
  explicit DisplacementControl(std::vector<int> p_controlledIndices)
      : controlledIndices{std::move(p_controlledIndices)} {}
 
  void operator()(SubsidiaryArgs& args) const {
    const auto controlledDOFsNorm = args.DD(controlledIndices).norm();
    args.f                        = controlledDOFsNorm - args.stepSize;
    args.dfdDlambda               = 0.0;
    args.dfdDD.setZero();
    args.dfdDD(controlledIndices) = args.DD(controlledIndices) / controlledDOFsNorm;
  }
 
  template <typename F>
  void initialPrediction(typename F::Domain& req, F& residual, SubsidiaryArgs& args) {
    args.DD(controlledIndices).array() = args.stepSize;
    req.globalSolution()               = args.DD;
  }
 
  template <typename F>
  void intermediatePrediction(typename F::Domain& req, F& residual, SubsidiaryArgs& args) {
    req.globalSolution() += args.DD;
  }
 
  constexpr auto name() const { return std::string("Displacement Control"); }
 
private:
  std::vector<int> controlledIndices; 
};
} // namespace Ikarus