qpenalty.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
namespace FEAT
{
  namespace Solver
  {
    template<typename Functional_>
    class QPenalty : public IterativeSolver<typename Functional_::VectorTypeR>
    {
      public:
        typedef Functional_ FunctionalType;
        typedef typename Functional_::VectorTypeR VectorType;
        typedef typename VectorType::DataType DataType;
        typedef IterativeSolver<VectorType> BaseClass;
 
      private:
        FunctionalType& _functional;
        std::shared_ptr<IterativeSolver<VectorType>> _inner_solver;
        DataType _initial_penalty_param;
        DataType _tol_penalty;
 
      public:
        explicit QPenalty(FunctionalType& functional, std::shared_ptr<IterativeSolver<VectorType>> inner_solver,
          DataType initial_penalty_param = DataType(1)) :
          BaseClass("QPenalty"),
          _functional(functional),
          _inner_solver(inner_solver),
          _initial_penalty_param(initial_penalty_param),
          _tol_penalty(Math::pow(Math::huge<DataType>(), DataType(0.25)))
        {
          // set communicator by functional (same interface as matrix)
          this->_set_comm_by_matrix(functional);
        }
 
        explicit QPenalty(const String& section_name, const PropertyMap* section,
          FunctionalType& functional, std::shared_ptr<IterativeSolver<VectorType>> inner_solver) :
          BaseClass("QPenalty", section_name, section),
          _functional(functional),
          _inner_solver(inner_solver),
          _initial_penalty_param(1),
          _tol_penalty(Math::pow(Math::huge<DataType>(), DataType(0.25)))
        {
          // set communicator by functional (same interface as matrix)
          this->_set_comm_by_matrix(functional);
 
          auto initial_penalty_param_p = section->query("initial_penalty_param");
          if(initial_penalty_param_p.second)
          {
            set_initial_penalty_param(DataType(std::stod(initial_penalty_param_p.first)));
          }
 
          auto tol_penalty_p = section->query("tol_penalty");
          if(tol_penalty_p.second)
          {
            set_tol_penalty(DataType(std::stod(tol_penalty_p.first)));
          }
        }
 
        QPenalty(const QPenalty&) = delete;
 
        virtual ~QPenalty()
        {
          _inner_solver = nullptr;
        }
 
        virtual String name() const override
        {
          return "QPenalty";
        }
 
        virtual void init_symbolic() override
        {
          BaseClass::init_symbolic();
          _inner_solver->init_symbolic();
        }
 
        virtual void done_symbolic() override
        {
          _inner_solver->done_symbolic();
          BaseClass::done_symbolic();
        }
 
        void set_initial_penalty_param(DataType initial_penalty_param)
        {
          XASSERT(initial_penalty_param > DataType(0));
 
          _initial_penalty_param = initial_penalty_param;
        }
 
        void set_tol_penalty(DataType tol_penalty)
        {
          XASSERT(tol_penalty > DataType(0));
 
          _tol_penalty = tol_penalty;
        }
 
        virtual Status apply(VectorType& vec_cor, const VectorType& vec_def) override
        {
          _functional.set_penalty_param(_initial_penalty_param);
 
          // clear solution vector
          vec_cor.format();
 
          // apply
          this->_status = _apply_intern(vec_cor, vec_def);
          this->plot_summary();
          return this->_status;
        }
 
        virtual Status correct(VectorType& vec_sol, const VectorType& vec_rhs) override
        {
          _functional.set_penalty_param(_initial_penalty_param);
 
          // apply
          this->_status = _apply_intern(vec_sol, vec_rhs);
          this->plot_summary();
          return this->_status;
        }
 
      protected:
        virtual Status _apply_intern(VectorType& vec_sol, const VectorType& vec_rhs)
        {
          Statistics::add_solver_expression(std::make_shared<ExpressionStartSolve>(this->name()));
 
          const Index inner_iter_digits(Math::ilog10(_inner_solver->get_max_iter()));
 
          Status st(_set_initial_defect(_functional.get_constraint()));
 
          IterationStats stat(*this);
 
          // The Penalty Iteration(TM)
          while(st == Status::progress)
          {
            ++(this->_num_iter);
 
            DataType def_old(this->_def_cur);
            //DataType penalty_param_old(penalty_param);
 
            Status inner_st(_inner_solver->correct(vec_sol, vec_rhs));
            if(inner_st == Status::aborted)
            {
              Statistics::add_solver_expression(std::make_shared<ExpressionEndSolve>(this->name(), Status::aborted, this->get_num_iter()));
              return Status::aborted;
            }
 
            // After the inner solver was called, the constraint is already there so we can just get() it
            this->_def_cur = _functional.get_constraint();
            DataType penalty_param(_functional.get_penalty_param());
 
            if(this->_plot_iter())
            {
              String msg = this->_plot_name
                + ": " + stringify(this->_num_iter).pad_front(this->_iter_digits)
                + " (" + stringify(this->_inner_solver->get_num_iter()).pad_front(inner_iter_digits)
                + ":" + stringify(inner_st) + ")"
                + " : " + stringify_fp_sci(this->_def_cur)
                + " / " + stringify_fp_sci(this->_def_cur / this->_def_init)
                + " / " + stringify_fp_sci(DataType(0.5)*Math::sqr(this->_def_cur))
                + " : " + stringify_fp_sci(penalty_param);
 
              this->_print_line(msg);
            }
 
            // ensure that the defect is neither NaN nor infinity
            if(!Math::isfinite(this->_def_cur))
            {
              Statistics::add_solver_expression(std::make_shared<ExpressionEndSolve>(this->name(), Status::aborted, this->get_num_iter()));
              return Status::aborted;
            }
 
            // is diverged?
            if(this->is_diverged())
            {
              Statistics::add_solver_expression(std::make_shared<ExpressionEndSolve>(this->name(), Status::diverged, this->get_num_iter()));
              return Status::diverged;
            }
 
            // minimum number of iterations performed?
            if(this->_num_iter < this->_min_iter)
            {
              Statistics::add_solver_expression(std::make_shared<ExpressionEndSolve>(this->name(), Status::progress, this->get_num_iter()));
              return Status::progress;
            }
 
            // maximum number of iterations performed?
            if(this->_num_iter >= this->_max_iter)
            {
              Statistics::add_solver_expression(std::make_shared<ExpressionEndSolve>(this->name(), Status::max_iter, this->get_num_iter()));
              return Status::max_iter;
            }
 
            // Check for convergence
            if(this->is_converged())
            {
              Statistics::add_solver_expression(std::make_shared<ExpressionEndSolve>(this->name(), Status::success, this->get_num_iter()));
              return Status::success;
            }
 
            // Increment for the penalty factor by at least factor 5, very arbitrary
            penalty_param = DataType(5)*Math::sqr(penalty_param)
              *Math::max(Math::sqr( def_old/this->_def_cur),DataType(1));
 
            if(penalty_param >= _tol_penalty)
            {
              Statistics::add_solver_expression(std::make_shared<ExpressionEndSolve>(this->name(), Status::stagnated, this->get_num_iter()));
              return Status::stagnated;
            }
 
            _functional.set_penalty_param(penalty_param);
          }
 
          // We should never come to this point
          Statistics::add_solver_expression(std::make_shared<ExpressionEndSolve>(this->name(), Status::undefined, this->get_num_iter()));
          return Status::undefined;
        }
 
      Status _set_initial_defect(const DataType def_init_)
      {
        // store new defect
        this->_def_init = this->_def_cur = def_init_;
        // insert special toe to signal new start of solver
        //Statistics::add_solver_toe(this->_branch, double(-1));
        //Statistics::add_solver_mpi_execute(this->_branch, double(-1));
        //Statistics::add_solver_mpi_wait(this->_branch, double(-1));
        //insert -1 as first defect, to signalize a new starting solver iteration run
        //Statistics::add_solver_defect(this->_branch, double(-1));
        //Statistics::add_solver_defect(this->_branch, double(this->_def_init));
        this->_num_iter = Index(0);
        this->_num_stag_iter = Index(0);
        Statistics::add_solver_expression(std::make_shared<ExpressionDefect>(this->name(), this->_def_init, this->get_num_iter()));
 
        // Plot?
        if(this->_plot_iter())
        {
          String msg = this->_plot_name
            +  ": " + stringify(0).pad_front(this->_iter_digits)
            + " : " + stringify_fp_sci(this->_def_init);
 
          this->_print_line(msg);
        }
 
        // Ensure that the defect is neither NaN nor infinity
        if(!Math::isfinite(this->_def_init))
        {
          return Status::aborted;
        }
 
        // Check if the initial defect lower than the absolute tolerance
        if(this->_def_init <= this->_tol_abs)
        {
          return Status::success;
        }
 
        // continue iterating
        return Status::progress;
      }
 
      using BaseClass::_set_initial_defect;
 
    };
 
    template<typename Functional_>
    inline std::shared_ptr<QPenalty<Functional_>> new_qpenalty( Functional_& functional,
    std::shared_ptr<IterativeSolver<typename Functional_::VectorTypeR>> inner_solver,
    typename Functional_::VectorTypeR::DataType initial_penalty_param =
    typename Functional_::VectorTypeR::DataType(1))
    {
      return std::make_shared<QPenalty<Functional_>>(functional, inner_solver, initial_penalty_param);
    }
 
    template<typename Functional_>
    inline std::shared_ptr<QPenalty<Functional_>> new_qpenalty(
      const String& section_name, const PropertyMap* section,
      Functional_& functional, std::shared_ptr<IterativeSolver<typename Functional_::VectorTypeR>> inner_solver)
    {
      return std::make_shared<QPenalty<Functional_>>(section_name, section, functional, inner_solver);
    }
 
  }
} // namespace FEAT