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