full_space_step.hpp

#ifndef __FULLSPACE_STEP_H__
#define __FULLSPACE_STEP_H__

#include "ROL_Types.hpp"
#include "ROL_Step.hpp"
#include "ROL_LineSearch.hpp"

// Unconstrained Methods
#include "ROL_GradientStep.hpp"
#include "ROL_NonlinearCGStep.hpp"
#include "ROL_SecantStep.hpp"
#include "ROL_NewtonStep.hpp"
#include "ROL_NewtonKrylovStep.hpp"

#include <sstream>
#include <iomanip>

#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/solver_control.h>

#include "optimization/rol_to_dealii_vector.hpp"
#include "optimization/dealii_solver_rol_vector.hpp"
#include "optimization/flow_constraints.hpp"
#include "optimization/rol_objective.hpp"

#include "optimization/kkt_operator.hpp"
#include "optimization/kkt_birosghattas_preconditioners.hpp"

namespace ROL {


template <class Real>
class FullSpace_BirosGhattas : public Step<Real> {
private:

    ROL::Ptr<Vector<Real> > design_variable_cloner_;
    ROL::Ptr<Vector<Real> > lagrange_variable_cloner_;


    ROL::Ptr<Objective<Real>> merit_function_;
    Real penalty_value_;

    ROL::Ptr<Vector<Real>> lagrange_mult_search_direction_;

    ROL::Ptr<Vector<Real>> previous_reduced_gradient_;

    ROL::Ptr<Secant<Real> > secant_;

    ROL::Ptr<LineSearch<Real> >  lineSearch_;
  
    ESecant esec_;
    ELineSearch els_;
    ECurvatureCondition econd_;
  

    int verbosity_;

    bool acceptLastAlpha_;
  
    ROL::ParameterList parlist_;
  
    std::string lineSearchName_;  
    std::string secantName_;

    bool use_approximate_full_space_preconditioner_;

    std::string preconditioner_name_;
  
    double search_ctl_norm;
    double search_sim_norm;
    double search_adj_norm;

    int n_linesearches;
public:
  
    using Step<Real>::initialize; 
    using Step<Real>::compute; 
    using Step<Real>::update; 
  
    FullSpace_BirosGhattas(
        ROL::ParameterList &parlist,
        const ROL::Ptr<LineSearch<Real> > &lineSearch = ROL::nullPtr,
        const ROL::Ptr<Secant<Real> > &secant = ROL::nullPtr);

    void computeLagrangianGradient(
        Vector<Real> &lagrangian_gradient,
        const Vector<Real> &design_variables,
        const Vector<Real> &lagrange_mult,
        const Vector<Real> &objective_gradient,
        Constraint<Real> &equal_constraints) const;

    void computeInitialLagrangeMultiplier(
        Vector<Real> &lagrange_mult,
        const Vector<Real> &design_variables,
        const Vector<Real> &objective_gradient,
        Constraint<Real> &equal_constraints) const;

  

    virtual void initialize(
        Vector<Real> &design_variables,
        const Vector<Real> &gradient,
        Vector<Real> &lagrange_mult,
        const Vector<Real> &equal_constraints_values,
        Objective<Real> &objective,
        Constraint<Real> &equal_constraints,
        AlgorithmState<Real> &algo_state ) override;


    void initialize(
        Vector<Real> &design_variables,
        const Vector<Real> &gradient,
        Vector<Real> &lagrange_mult,
        const Vector<Real> &equal_constraints_values,
        Objective<Real> &objective,
        Constraint<Real> &equal_constraints,
        BoundConstraint<Real> &bound_constraints,
        AlgorithmState<Real> &algo_state ) override;

    Real computeAugmentedLagrangianPenalty(
        const Vector<Real> &search_direction,
        const Vector<Real> &lagrange_mult_search_direction,
        const Vector<Real> &design_variables,
        const Vector<Real> &objective_gradient,
        const Vector<Real> &equal_constraints_values,
        const Vector<Real> &adjoint_jacobian_lagrange,
        Constraint<Real> &equal_constraints,
        const Real offset);

    template<typename MatrixType, typename VectorType, typename PreconditionerType>
    std::vector<double>
    solve_linear (
        MatrixType &matrix_A,
        VectorType &right_hand_side,
        VectorType &solution,
        PreconditionerType &preconditioner);
        //const PHiLiP::Parameters::LinearSolverParam & param = );

    std::vector<Real> solve_KKT_system(
        Vector<Real> &search_direction,
        Vector<Real> &lag_search_direction,
        const Vector<Real> &design_variables,
        const Vector<Real> &lagrange_mult,
        Objective<Real> &objective,
        Constraint<Real> &equal_constraints);


    void compute(
        Vector<Real> &search_direction,
        const Vector<Real> &design_variables,
        const Vector<Real> &lagrange_mult,
        Objective<Real> &objective,
        Constraint<Real> &equal_constraints,
        AlgorithmState<Real> &algo_state ) override;

    void compute(
        Vector<Real> &search_direction,
        const Vector<Real> &design_variables,
        const Vector<Real> &lagrange_mult,
        Objective<Real> &objective,
        Constraint<Real> &equal_constraints,
        BoundConstraint<Real> &bound_constraints,
        AlgorithmState<Real> &algo_state ) override;
  
    void update(
        Vector<Real> &design_variables,
        Vector<Real> &lagrange_mult,
        const Vector<Real> &search_direction,
        Objective<Real> &objective,
        Constraint<Real> &equal_constraints,
        AlgorithmState<Real> &algo_state ) override;

    void update (
        Vector<Real> &design_variables,
        Vector<Real> &lagrange_mult,
        const Vector<Real> &search_direction,
        Objective<Real> &objective,
        Constraint<Real> &equal_constraints,
        BoundConstraint< Real > &bound_constraints,
        AlgorithmState<Real> &algo_state ) override;
  
    std::string printHeader( void ) const override;
    
    std::string printName( void ) const override;
  
    std::string print( AlgorithmState<Real> & algo_state, bool print_header = false ) const override;

private:
    dealii::ConditionalOStream pcout; 

}; // class FullSpace_BirosGhattas

} // namespace ROL
#endif