1 
2 #ifndef MISFIT_ROBJ_H
3 #define MISFIT_ROBJ_H
4 
5 #include "ROL_Objective.hpp"
6 #include "ROL_SimController.hpp"
7 #include "femdata.hpp"
8 
9 template <class Real>
10 class Misfit_Objective : public ROL::Objective<Real> {
11 private:
12   // Advection Diffusion PDE
13   const ROL::Ptr<FEMdata<Real>> fem_;
14 
15   // Vector Storage
16   ROL::Ptr<ROL::SimController<Real>> stateStore_, adjointStore_;
17   ROL::Ptr<ROL::Vector<Real>> state_, adjoint_, state_sens_, adjoint_sens_;
18   ROL::Ptr<ROL::Vector<Real>> dualadjoint_;
19 
20   bool storage_;
21 
22   unsigned nstat_, nadjo_, nsens_, nsadj_;
23   unsigned nupda_, nfval_, ngrad_, nhess_, nprec_;
24 
25 public:
Misfit_Objective(const ROL::Ptr<FEMdata<Real>> & fem,ROL::ParameterList & list)26   Misfit_Objective(const ROL::Ptr<FEMdata<Real>> &fem,
27                    ROL::ParameterList            &list)
28     : fem_(fem), nstat_(0), nadjo_(0), nsens_(0), nsadj_(0),
29       nupda_(0), nfval_(0), ngrad_(0), nhess_(0), nprec_(0) {
30     stateStore_   = ROL::makePtr<ROL::SimController<Real>>();
31     adjointStore_ = ROL::makePtr<ROL::SimController<Real>>();
32 
33     // Vector Storage
34     state_        = fem_->createStateVector(list);
35     adjoint_      = fem_->createStateVector(list);
36     state_sens_   = fem_->createStateVector(list);
37     adjoint_sens_ = fem_->createStateVector(list);
38     dualadjoint_  = fem_->createResidualVector(list);
39 
40     storage_  = list.sublist("Problem").get("Use state storage", true);
41   }
42 
getAssembler(void) const43   const ROL::Ptr<Assembler<Real>> getAssembler(void) const {
44     return fem_->getAssembler();
45   }
46 
update(const ROL::Vector<Real> & z,bool flag=true,int iter=-1)47   void update( const ROL::Vector<Real> &z, bool flag = true, int iter = -1 ) {
48     nupda_++;
49     stateStore_->objectiveUpdate(true);
50     adjointStore_->objectiveUpdate(flag);
51   }
52 
solvePDE(ROL::Vector<Real> & u,const ROL::Vector<Real> & z)53   void solvePDE(ROL::Vector<Real> &u, const ROL::Vector<Real> &z) {
54     update(z,false);
55     solve_state_equation(u,z);
56   }
57 
value(const ROL::Vector<Real> & z,Real & tol)58   Real value( const ROL::Vector<Real> &z, Real &tol ) {
59     nfval_++;
60     solve_state_equation(*state_,z);
61     fem_->applyObjectiveHessian(*dualadjoint_,*state_); // Hu
62     dualadjoint_->scale(static_cast<Real>(0.5));        // 0.5 Hu
63     fem_->addObjectiveGradient(*dualadjoint_);          // 0.5 Hu + g
64     // 0.5 uHu + gu + c0
65     return state_->dot(dualadjoint_->dual()) + fem_->getObjectiveConstant();
66   }
67 
gradient(ROL::Vector<Real> & g,const ROL::Vector<Real> & z,Real & tol)68   void gradient( ROL::Vector<Real> &g, const ROL::Vector<Real> &z, Real &tol ) {
69     ngrad_++;
70     solve_state_equation(*state_,z);
71     solve_adjoint_equation(*adjoint_,*state_);
72     fem_->applyControlJacobian(g,*adjoint_,true);
73   }
74 
hessVec(ROL::Vector<Real> & hv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & z,Real & tol)75   void hessVec( ROL::Vector<Real> &hv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &z, Real &tol ) {
76     nhess_++;
77     solve_state_sensitivity(*state_sens_,v);
78     solve_adjoint_sensitivity(*adjoint_sens_,*state_sens_);
79     fem_->applyControlJacobian(hv,*adjoint_sens_,true);
80   }
81 
82   /** \brief Apply a reduced Hessian preconditioner.
83   */
precond(ROL::Vector<Real> & Pv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & z,Real & tol)84   virtual void precond( ROL::Vector<Real> &Pv, const ROL::Vector<Real> &v, const ROL::Vector<Real> &z, Real &tol ) {
85     nprec_++;
86     Pv.set(v.dual());
87   }
88 
summarize(std::ostream & stream) const89   void summarize(std::ostream &stream) const {
90     stream << std::endl;
91     stream << std::string(114,'=') << std::endl;
92     stream << "  Misfit_Objective::summarize" << std::endl;
93     stream << "    Number of calls to update:            " << nupda_ << std::endl;
94     stream << "    Number of calls to value:             " << nfval_ << std::endl;
95     stream << "    Number of calls to gradient:          " << ngrad_ << std::endl;
96     stream << "    Number of calls to hessVec:           " << nhess_ << std::endl;
97     stream << "    Number of calls to precond:           " << nprec_ << std::endl;
98     stream << "    Number of state solves:               " << nstat_ << std::endl;
99     stream << "    Number of adjoint solves:             " << nadjo_ << std::endl;
100     stream << "    Number of state sensitivity solves:   " << nsens_ << std::endl;
101     stream << "    Number of adjoint sensitivity solves: " << nsadj_ << std::endl;
102     stream << std::string(114,'=') << std::endl;
103     stream << std::endl;
104     fem_->summarize(stream);
105   }
106 
107 private:
108 
solve_state_equation(ROL::Vector<Real> & state,const ROL::Vector<Real> & control)109   void solve_state_equation(ROL::Vector<Real> &state, const ROL::Vector<Real> &control) {
110     bool isComputed = false;
111     if (storage_) {
112       isComputed = stateStore_->get(state,ROL::Objective<Real>::getParameter());
113     }
114     if (!isComputed || !storage_) {
115       fem_->applyControlJacobian(*dualadjoint_,control,false);  // Bz
116       fem_->addPDErhs(*dualadjoint_);                           // Bz + f
117       fem_->applyInversePDEJacobian(state,*dualadjoint_,false); // inv(A)(Bz + f)
118       state.scale(static_cast<Real>(-1));                       // -inv(A)(Bz + f)
119       nstat_++;
120       if (storage_) {
121         stateStore_->set(state,ROL::Objective<Real>::getParameter());
122       }
123     }
124   }
125 
solve_adjoint_equation(ROL::Vector<Real> & adjoint,const ROL::Vector<Real> & state)126   void solve_adjoint_equation(ROL::Vector<Real> &adjoint, const ROL::Vector<Real> &state) {
127     bool isComputed = false;
128     if (storage_) {
129       isComputed = adjointStore_->get(adjoint,ROL::Objective<Real>::getParameter());
130     }
131     if (!isComputed || !storage_) {
132       fem_->applyObjectiveHessian(*dualadjoint_,state);          // Hu
133       fem_->addObjectiveGradient(*dualadjoint_);                 // Hu + g
134       fem_->applyInversePDEJacobian(adjoint,*dualadjoint_,true); // inv(A')(Hu + g)
135       adjoint.scale(static_cast<Real>(-1));                      // -inv(A')(Hu + g)
136       nadjo_++;
137       if (storage_) {
138         adjointStore_->set(adjoint,ROL::Objective<Real>::getParameter());
139       }
140     }
141   }
142 
solve_state_sensitivity(ROL::Vector<Real> & state_sens,const ROL::Vector<Real> & v)143   void solve_state_sensitivity(ROL::Vector<Real> &state_sens, const ROL::Vector<Real> &v) {
144     fem_->applyControlJacobian(*dualadjoint_,v,false);               // Bv
145     fem_->applyInversePDEJacobian(state_sens,*dualadjoint_,false); // inv(A)Bv
146     nsens_++;
147   }
148 
solve_adjoint_sensitivity(ROL::Vector<Real> & adjoint_sens,const ROL::Vector<Real> & state_sens)149   void solve_adjoint_sensitivity(ROL::Vector<Real> &adjoint_sens, const ROL::Vector<Real> &state_sens) {
150     fem_->applyObjectiveHessian(*dualadjoint_,state_sens);          // Hv
151     fem_->applyInversePDEJacobian(adjoint_sens,*dualadjoint_,true); // inv(A')Hv
152     nsadj_++;
153   }
154 }; // class Misfit_Objective
155 
156 #endif
157