flow_constraints.hpp

#ifndef __FLOWCONSTRAINTS_H__
#define __FLOWCONSTRAINTS_H__

#include <deal.II/optimization/rol/vector_adaptor.h>

#include "Ifpack.h"
#include "ROL_Constraint_SimOpt.hpp"
#include "design_parameterization/base_parameterization.hpp"
#include "dg/dg_base.hpp"
#include "linear_solver/linear_solver.h"
#include "parameters/all_parameters.h"

namespace PHiLiP {

using dealii_Vector = dealii::LinearAlgebra::distributed::Vector<double>;
using AdaptVector = dealii::Rol::VectorAdaptor<dealii_Vector>;


template<int dim>
class FlowConstraints : public ROL::Constraint_SimOpt<double> {
private:
    const int mpi_rank;
    const bool i_print;
    std::shared_ptr<DGBase<dim,double>> dg;

    std::shared_ptr<BaseParameterization<dim>> design_parameterization;


    Parameters::LinearSolverParam linear_solver_param;

    dealii::LinearAlgebra::distributed::Vector<double> design_var;


    Ifpack_Preconditioner *jacobian_prec;

    Ifpack_Preconditioner *adjoint_jacobian_prec;

protected:
    int i_out = 1000;
    int iupdate = 9000;

public:
    dealii::TrilinosWrappers::SparseMatrix dXvdXp;

    using ROL::Constraint_SimOpt<double>::value;

    double flow_CFL_;
    using ROL::Constraint_SimOpt<double>::applyAdjointJacobian_1;
        //(
        //ROL::Vector<double>& output_vector,
        //const ROL::Vector<double>& input_vector,
        //const ROL::Vector<double>& des_var_sim,
        //const ROL::Vector<double>& des_var_ctl,
        //const ROL::Vector<double>& dualv,
        //double& /*tol*/ 
        //);
    using ROL::Constraint_SimOpt<double>::applyAdjointJacobian_2;

    // using ROL::Constraint_SimOpt<double>::applyAdjointHessian_11;
    // using ROL::Constraint_SimOpt<double>::applyAdjointHessian_12;
    // using ROL::Constraint_SimOpt<double>::applyAdjointHessian_21;
    // using ROL::Constraint_SimOpt<double>::applyAdjointHessian_22;

    FlowConstraints(
        std::shared_ptr<DGBase<dim,double>> &_dg,
        std::shared_ptr<BaseParameterization<dim>> _design_parameterization,
        std::shared_ptr<dealii::TrilinosWrappers::SparseMatrix> precomputed_dXvdXp = nullptr);

    ~FlowConstraints();

    void update_1( const ROL::Vector<double>& des_var_sim, bool flag = true, int iter = -1 );


    void update_2( const ROL::Vector<double>& des_var_ctl, bool flag = true, int iter = -1 );


    void solve(
        ROL::Vector<double>& constraint_values,
        ROL::Vector<double>& des_var_sim,
        const ROL::Vector<double>& des_var_ctl,
        double& /*tol*/ 
        ) override;

    void value(
        ROL::Vector<double>& constraint_values,
        const ROL::Vector<double>& des_var_sim,
        const ROL::Vector<double>& des_var_ctl,
        double &/*tol*/ 
        ) override;
    
    void applyJacobian_1(
        ROL::Vector<double>& output_vector,
        const ROL::Vector<double>& input_vector,
        const ROL::Vector<double>& des_var_sim,
        const ROL::Vector<double>& des_var_ctl,
        double& /*tol*/ 
        ) override;

    void applyAdjointJacobian_1(
        ROL::Vector<double>& output_vector,
        const ROL::Vector<double>& input_vector,
        const ROL::Vector<double>& des_var_sim,
        const ROL::Vector<double>& des_var_ctl,
        double& /*tol*/ 
        ) override;

    void applyInverseJacobian_1(
        ROL::Vector<double>& output_vector,
        const ROL::Vector<double>& input_vector,
        const ROL::Vector<double>& des_var_sim,
        const ROL::Vector<double>& des_var_ctl,
        double& /*tol*/ 
        ) override;

    int construct_JacobianPreconditioner_1(
        const ROL::Vector<double>& des_var_sim,
        const ROL::Vector<double>& des_var_ctl);

    void destroy_JacobianPreconditioner_1();


    void applyInverseJacobianPreconditioner_1(
        ROL::Vector<double>& output_vector,
        const ROL::Vector<double>& input_vector,
        const ROL::Vector<double>& des_var_sim,
        const ROL::Vector<double>& des_var_ctl,
        double& /*tol*/ 
        );

    int construct_AdjointJacobianPreconditioner_1(
        const ROL::Vector<double>& des_var_sim,
        const ROL::Vector<double>& des_var_ctl);

    void destroy_AdjointJacobianPreconditioner_1();


    void applyInverseAdjointJacobianPreconditioner_1(
        ROL::Vector<double>& output_vector,
        const ROL::Vector<double>& input_vector,
        const ROL::Vector<double>& des_var_sim,
        const ROL::Vector<double>& des_var_ctl,
        double& /*tol*/ 
        );

    void applyInverseAdjointJacobian_1(
        ROL::Vector<double>& output_vector,
        const ROL::Vector<double>& input_vector,
        const ROL::Vector<double>& des_var_sim,
        const ROL::Vector<double>& des_var_ctl,
        double& /*tol*/ 
        ) override;

    void applyJacobian_2(
        ROL::Vector<double>& output_vector,
        const ROL::Vector<double>& input_vector,
        const ROL::Vector<double>& des_var_sim,
        const ROL::Vector<double>& des_var_ctl,
        double& /*tol*/ 
        ) override;


    void applyAdjointJacobian_2(
        ROL::Vector<double>& output_vector,
        const ROL::Vector<double>& input_vector,
        const ROL::Vector<double>& des_var_sim,
        const ROL::Vector<double>& des_var_ctl,
        double& /*tol*/
        ) override;


    void applyAdjointHessian_11 (
        ROL::Vector<double> &output_vector,
        const ROL::Vector<double> &dual,
        const ROL::Vector<double> &input_vector,
        const ROL::Vector<double> &des_var_sim,
        const ROL::Vector<double> &des_var_ctl,
        double &tol
        ) override;


    void applyAdjointHessian_12 (
        ROL::Vector<double> &output_vector,
        const ROL::Vector<double> &dual,
        const ROL::Vector<double> &input_vector,
        const ROL::Vector<double> &des_var_sim,
        const ROL::Vector<double> &des_var_ctl,
        double &tol
        ) override;


    void applyAdjointHessian_21 (
        ROL::Vector<double> &output_vector,
        const ROL::Vector<double> &dual,
        const ROL::Vector<double> &input_vector,
        const ROL::Vector<double> &des_var_sim,
        const ROL::Vector<double> &des_var_ctl,
        double &tol
        ) override;


    void applyAdjointHessian_22 (
        ROL::Vector<double> &output_vector,
        const ROL::Vector<double> &dual,
        const ROL::Vector<double> &input_vector,
        const ROL::Vector<double> &des_var_sim,
        const ROL::Vector<double> &des_var_ctl,
        double &tol
        ) override;

    // std::vector<double> solveAugmentedSystem(
    //     ROL::Vector<double> &v1,
    //     ROL::Vector<double> &v2,
    //     const ROL::Vector<double> &b1,
    //     const ROL::Vector<double> &b2,
    //     const ROL::Vector<double> &x,
    //     double & tol) override;
    //
    // void applyPreconditioner(ROL::Vector<double> &pv,
    //                          const ROL::Vector<double> &v,
    //                          const ROL::Vector<double> &x,
    //                          const ROL::Vector<double> &g,
    //                          double &tol) override;

};

} // PHiLiP namespace

#endif