1 // @HEADER 2 // ************************************************************************ 3 // 4 // Rapid Optimization Library (ROL) Package 5 // Copyright (2014) Sandia Corporation 6 // 7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 8 // license for use of this work by or on behalf of the U.S. Government. 9 // 10 // Redistribution and use in source and binary forms, with or without 11 // modification, are permitted provided that the following conditions are 12 // met: 13 // 14 // 1. Redistributions of source code must retain the above copyright 15 // notice, this list of conditions and the following disclaimer. 16 // 17 // 2. Redistributions in binary form must reproduce the above copyright 18 // notice, this list of conditions and the following disclaimer in the 19 // documentation and/or other materials provided with the distribution. 20 // 21 // 3. Neither the name of the Corporation nor the names of the 22 // contributors may be used to endorse or promote products derived from 23 // this software without specific prior written permission. 24 // 25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY 26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE 29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 36 // 37 // Questions? Contact lead developers: 38 // Drew Kouri (dpkouri@sandia.gov) and 39 // Denis Ridzal (dridzal@sandia.gov) 40 // 41 // ************************************************************************ 42 // @HEADER 43 44 /*! \file constraint.hpp 45 \brief Defines the SimOpt constraint for the 'poisson' example. 46 */ 47 48 #ifndef ROL_PDEOPT_DYNCONSTRAINT_H 49 #define ROL_PDEOPT_DYNCONSTRAINT_H 50 51 #include "ROL_DynamicConstraint.hpp" 52 #include "dynpde.hpp" 53 #include "assembler.hpp" 54 #include "solver.hpp" 55 #include "pdevector.hpp" 56 57 // Do not instantiate the template in this translation unit. 58 extern template class Assembler<double>; 59 60 //// Global Timers. 61 #ifdef ROL_TIMERS 62 namespace ROL { 63 namespace PDEOPT { 64 ROL::Ptr<Teuchos::Time> DynConstraintSolverConstruct_Jacobian_un = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Solver Construction Time for Jacobian un"); 65 ROL::Ptr<Teuchos::Time> DynConstraintSolverSolve_Jacobian_un = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Solver Solution Time for Jacobian un"); 66 ROL::Ptr<Teuchos::Time> DynConstraintSolverSolve_AdjointJacobian_un = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Solver Solution Time for Adjoint Jacobian un"); 67 ROL::Ptr<Teuchos::Time> DynConstraintApplyJacobian_uo = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Jacobian uo"); 68 ROL::Ptr<Teuchos::Time> DynConstraintApplyJacobian_un = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Jacobian un"); 69 ROL::Ptr<Teuchos::Time> DynConstraintApplyJacobian_zf = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Jacobian zf"); 70 ROL::Ptr<Teuchos::Time> DynConstraintApplyJacobian_zp = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Jacobian zp"); 71 ROL::Ptr<Teuchos::Time> DynConstraintApplyAdjointJacobian_uo = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Adjoint Jacobian uo"); 72 ROL::Ptr<Teuchos::Time> DynConstraintApplyAdjointJacobian_un = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Adjoint Jacobian un"); 73 ROL::Ptr<Teuchos::Time> DynConstraintApplyAdjointJacobian_zf = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Adjoint Jacobian zf"); 74 ROL::Ptr<Teuchos::Time> DynConstraintApplyAdjointJacobian_zp = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Adjoint Jacobian zp"); 75 ROL::Ptr<Teuchos::Time> DynConstraintApplyHessian_uo_uo = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Hessian_uo_uo"); 76 ROL::Ptr<Teuchos::Time> DynConstraintApplyHessian_uo_un = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Hessian_uo_un"); 77 ROL::Ptr<Teuchos::Time> DynConstraintApplyHessian_un_uo = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Hessian_un_uo"); 78 ROL::Ptr<Teuchos::Time> DynConstraintApplyHessian_un_un = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Hessian_un_un"); 79 ROL::Ptr<Teuchos::Time> DynConstraintApplyHessian_uo_zf = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Hessian_uo_zf"); 80 ROL::Ptr<Teuchos::Time> DynConstraintApplyHessian_uo_zp = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Hessian_uo_zp"); 81 ROL::Ptr<Teuchos::Time> DynConstraintApplyHessian_zf_uo = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Hessian_zf_uo"); 82 ROL::Ptr<Teuchos::Time> DynConstraintApplyHessian_zp_uo = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Hessian_zp_uo"); 83 ROL::Ptr<Teuchos::Time> DynConstraintApplyHessian_un_zf = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Hessian_un_zf"); 84 ROL::Ptr<Teuchos::Time> DynConstraintApplyHessian_un_zp = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Hessian_un_zp"); 85 ROL::Ptr<Teuchos::Time> DynConstraintApplyHessian_zf_un = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Hessian_zf_un"); 86 ROL::Ptr<Teuchos::Time> DynConstraintApplyHessian_zp_un = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Hessian_zp_un"); 87 ROL::Ptr<Teuchos::Time> DynConstraintApplyHessian_zp_zp = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Hessian_zp_zp"); 88 ROL::Ptr<Teuchos::Time> DynConstraintApplyHessian_zp_zf = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Hessian_zp_zf"); 89 ROL::Ptr<Teuchos::Time> DynConstraintApplyHessian_zf_zp = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Hessian_zf_zp"); 90 ROL::Ptr<Teuchos::Time> DynConstraintApplyHessian_zf_zf = Teuchos::TimeMonitor::getNewCounter("ROL::PDEOPT: Dyn Constraint Apply Hessian_zf_zf"); 91 } 92 } 93 #endif 94 95 96 template<class Real> 97 class DynConstraint : public ROL::DynamicConstraint<Real> { 98 private: 99 const ROL::Ptr<DynamicPDE<Real>> pde_; 100 ROL::Ptr<Solver<Real>> solver_; 101 ROL::Ptr<Assembler<Real>> assembler_; 102 103 mutable ROL::Ptr<Tpetra::MultiVector<>> cvec_; 104 105 // Storage for spatial PDE residual 106 mutable ROL::Ptr<Tpetra::MultiVector<>> vecR_; 107 mutable bool isResAssembled_; 108 109 // Storage for spatial PDE jacobian at old time 110 mutable ROL::Ptr<Tpetra::CrsMatrix<>> matJuo_; 111 mutable bool isJuoAssembled_, isJuoZero_, isJuoNotImplemented_; 112 113 // Storage for spatial PDE jacobian at new time 114 mutable ROL::Ptr<Tpetra::CrsMatrix<>> matJun_; 115 mutable bool isJunAssembled_, isJunZero_, isJunNotImplemented_; 116 117 // Storage for field control jacobian 118 mutable ROL::Ptr<Tpetra::CrsMatrix<>> matJzf_; 119 mutable bool isJzfAssembled_, isJzfZero_, isJzfNotImplemented_; 120 121 // Storage for parameter control jacobian 122 mutable ROL::Ptr<Tpetra::MultiVector<>> vecJzp_; 123 mutable bool isJzpAssembled_, isJzpZero_, isJzpNotImplemented_; 124 125 // Storage for spatial PDE old/old hessian 126 mutable ROL::Ptr<Tpetra::CrsMatrix<>> matHuo_uo_; 127 mutable bool isHuo_uoAssembled_, isHuo_uoZero_, isHuo_uoNotImplemented_; 128 129 // Storage for spatial PDE old/field hessian 130 mutable ROL::Ptr<Tpetra::CrsMatrix<>> matHuo_zf_; 131 mutable ROL::Ptr<Tpetra::CrsMatrix<>> matHzf_uo_; 132 mutable bool isHuo_zfAssembled_, isHuo_zfZero_, isHuo_zfNotImplemented_; 133 mutable bool isHzf_uoAssembled_, isHzf_uoZero_, isHzf_uoNotImplemented_; 134 135 // Storage for spatial PDE old/parameter hessian 136 mutable ROL::Ptr<Tpetra::MultiVector<>> vecHuo_zp_; 137 mutable ROL::Ptr<Tpetra::MultiVector<>> vecHzp_uo_; 138 mutable bool isHuo_zpAssembled_, isHuo_zpZero_, isHuo_zpNotImplemented_; 139 mutable bool isHzp_uoAssembled_, isHzp_uoZero_, isHzp_uoNotImplemented_; 140 141 // Storage for spatial PDE new/old hessian 142 mutable ROL::Ptr<Tpetra::CrsMatrix<>> matHuo_un_; 143 mutable ROL::Ptr<Tpetra::CrsMatrix<>> matHun_uo_; 144 mutable bool isHuo_unAssembled_, isHuo_unZero_, isHuo_unNotImplemented_; 145 mutable bool isHun_uoAssembled_, isHun_uoZero_, isHun_uoNotImplemented_; 146 147 // Storage for spatial PDE new/new hessian 148 mutable ROL::Ptr<Tpetra::CrsMatrix<>> matHun_un_; 149 mutable bool isHun_unAssembled_, isHun_unZero_, isHun_unNotImplemented_; 150 151 // Storage for spatial PDE new/field hessian 152 mutable ROL::Ptr<Tpetra::CrsMatrix<>> matHun_zf_; 153 mutable ROL::Ptr<Tpetra::CrsMatrix<>> matHzf_un_; 154 mutable bool isHun_zfAssembled_, isHun_zfZero_, isHun_zfNotImplemented_; 155 mutable bool isHzf_unAssembled_, isHzf_unZero_, isHzf_unNotImplemented_; 156 157 // Storage for spatial PDE new/parameter hessian 158 mutable ROL::Ptr<Tpetra::MultiVector<>> vecHun_zp_; 159 mutable ROL::Ptr<Tpetra::MultiVector<>> vecHzp_un_; 160 mutable bool isHun_zpAssembled_, isHun_zpZero_, isHun_zpNotImplemented_; 161 mutable bool isHzp_unAssembled_, isHzp_unZero_, isHzp_unNotImplemented_; 162 163 // Storage for field/field hessian 164 mutable ROL::Ptr<Tpetra::CrsMatrix<>> matHzf_zf_; 165 mutable bool isHzf_zfAssembled_, isHzf_zfZero_, isHzf_zfNotImplemented_; 166 167 // Storage for parameter/parameter hessian 168 mutable ROL::Ptr<std::vector<std::vector<Real>>> matHzp_zp_; 169 mutable bool isHzp_zpAssembled_, isHzp_zpZero_, isHzp_zpNotImplemented_; 170 171 // Storage for field/parameter hessian 172 mutable ROL::Ptr<Tpetra::MultiVector<>> vecHzf_zp_; 173 mutable ROL::Ptr<Tpetra::MultiVector<>> vecHzp_zf_; 174 mutable bool isHzf_zpAssembled_, isHzf_zpZero_, isHzf_zpNotImplemented_; 175 mutable bool isHzp_zfAssembled_, isHzp_zfZero_, isHzp_zfNotImplemented_; 176 assembleResidual(const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const177 void assembleResidual(const ROL::Vector<Real> &uo, 178 const ROL::Vector<Real> &un, 179 const ROL::Vector<Real> &z, 180 const ROL::TimeStamp<Real> &ts) const { 181 if (!isResAssembled_) { 182 ROL::Ptr<const Tpetra::MultiVector<>> uof = getConstField(uo); 183 ROL::Ptr<const Tpetra::MultiVector<>> unf = getConstField(un); 184 ROL::Ptr<const Tpetra::MultiVector<>> zf = getConstField(z); 185 ROL::Ptr<const std::vector<Real>> zp = getConstParameter(z); 186 187 assembler_->assembleDynPDEResidual(vecR_,pde_,ts,uof,unf,zf,zp); 188 isResAssembled_ = true; 189 } 190 } 191 setSolver(void) const192 void setSolver(void) const { 193 #ifdef ROL_TIMERS 194 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintSolverConstruct_Jacobian_un); 195 #endif 196 solver_->setA(matJun_); 197 } 198 199 // Assemble the PDE Jacobian. assembleJuo(const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const200 void assembleJuo(const ROL::Vector<Real> &uo, 201 const ROL::Vector<Real> &un, 202 const ROL::Vector<Real> &z, 203 const ROL::TimeStamp<Real> &ts) const { 204 if ( !isJuoZero_ && !isJuoNotImplemented_ ) { 205 try { 206 if (!isJuoAssembled_) { 207 ROL::Ptr<const Tpetra::MultiVector<>> uof = getConstField(uo); 208 ROL::Ptr<const Tpetra::MultiVector<>> unf = getConstField(un); 209 ROL::Ptr<const Tpetra::MultiVector<>> zf = getConstField(z); 210 ROL::Ptr<const std::vector<Real>> zp = getConstParameter(z); 211 212 assembler_->assembleDynPDEJacobian_uo(matJuo_,pde_,ts,uof,unf,zf,zp); 213 isJuoAssembled_ = true; 214 } 215 } 216 catch ( Exception::Zero & ez ) { 217 isJuoZero_ = true; 218 } 219 catch ( Exception::NotImplemented & eni ) { 220 isJuoNotImplemented_ = true; 221 } 222 } 223 } 224 assembleJun(const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const225 void assembleJun(const ROL::Vector<Real> &uo, 226 const ROL::Vector<Real> &un, 227 const ROL::Vector<Real> &z, 228 const ROL::TimeStamp<Real> &ts) const { 229 if ( !isJunZero_ && !isJunNotImplemented_ ) { 230 try { 231 if (!isJunAssembled_) { 232 ROL::Ptr<const Tpetra::MultiVector<>> uof = getConstField(uo); 233 ROL::Ptr<const Tpetra::MultiVector<>> unf = getConstField(un); 234 ROL::Ptr<const Tpetra::MultiVector<>> zf = getConstField(z); 235 ROL::Ptr<const std::vector<Real>> zp = getConstParameter(z); 236 237 assembler_->assembleDynPDEJacobian_un(matJun_,pde_,ts,uof,unf,zf,zp); 238 isJunAssembled_ = true; 239 setSolver(); 240 } 241 } 242 catch ( Exception::Zero & ez ) { 243 isJunZero_ = true; 244 } 245 catch ( Exception::NotImplemented & eni ) { 246 isJunNotImplemented_ = true; 247 } 248 } 249 } 250 251 // Assemble the PDE Jacobian. assembleJzf(const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const252 void assembleJzf(const ROL::Vector<Real> &uo, 253 const ROL::Vector<Real> &un, 254 const ROL::Vector<Real> &z, 255 const ROL::TimeStamp<Real> &ts) const { 256 if ( !isJzfZero_ && !isJzfNotImplemented_ ) { 257 try { 258 if (!isJzfAssembled_) { 259 ROL::Ptr<const Tpetra::MultiVector<>> uof = getConstField(uo); 260 ROL::Ptr<const Tpetra::MultiVector<>> unf = getConstField(un); 261 ROL::Ptr<const Tpetra::MultiVector<>> zf = getConstField(z); 262 ROL::Ptr<const std::vector<Real>> zp = getConstParameter(z); 263 264 assembler_->assembleDynPDEJacobian_zf(matJzf_,pde_,ts,uof,unf,zf,zp); 265 isJzfAssembled_ = true; 266 } 267 } 268 catch ( Exception::Zero & ez ) { 269 isJzfZero_ = true; 270 } 271 catch ( Exception::NotImplemented & eni ) { 272 isJzfNotImplemented_ = true; 273 } 274 } 275 } 276 277 // Assemble the PDE Jacobian. assembleJzp(const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const278 void assembleJzp(const ROL::Vector<Real> &uo, 279 const ROL::Vector<Real> &un, 280 const ROL::Vector<Real> &z, 281 const ROL::TimeStamp<Real> &ts) const { 282 if ( !isJzpZero_ && !isJzpNotImplemented_ ) { 283 try { 284 if (!isJzpAssembled_) { 285 ROL::Ptr<const Tpetra::MultiVector<>> uof = getConstField(uo); 286 ROL::Ptr<const Tpetra::MultiVector<>> unf = getConstField(un); 287 ROL::Ptr<const Tpetra::MultiVector<>> zf = getConstField(z); 288 ROL::Ptr<const std::vector<Real>> zp = getConstParameter(z); 289 290 assembler_->assembleDynPDEJacobian_zp(vecJzp_,pde_,ts,uof,unf,zf,zp); 291 isJzpAssembled_ = true; 292 } 293 } 294 catch ( Exception::Zero & ez ) { 295 isJzpZero_ = true; 296 } 297 catch ( Exception::NotImplemented & eni ) { 298 isJzpNotImplemented_ = true; 299 } 300 } 301 } 302 303 // Assemble the PDE adjoint Hessian. assembleHuo_uo(const ROL::Vector<Real> & w,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const304 void assembleHuo_uo(const ROL::Vector<Real> &w, 305 const ROL::Vector<Real> &uo, 306 const ROL::Vector<Real> &un, 307 const ROL::Vector<Real> &z, 308 const ROL::TimeStamp<Real> &ts) const { 309 if ( !isHuo_uoZero_ && !isHuo_uoNotImplemented_ ) { 310 try { 311 if (!isHuo_uoAssembled_) { 312 ROL::Ptr<const Tpetra::MultiVector<>> wf = getConstField(w); 313 ROL::Ptr<const Tpetra::MultiVector<>> uof = getConstField(uo); 314 ROL::Ptr<const Tpetra::MultiVector<>> unf = getConstField(un); 315 ROL::Ptr<const Tpetra::MultiVector<>> zf = getConstField(z); 316 ROL::Ptr<const std::vector<Real>> zp = getConstParameter(z); 317 318 assembler_->assembleDynPDEHessian_uo_uo(matHuo_uo_,pde_,ts,wf,uof,unf,zf,zp); 319 isHuo_uoAssembled_ = true; 320 } 321 } 322 catch ( Exception::Zero & ez ) { 323 isHuo_uoZero_ = true; 324 } 325 catch ( Exception::NotImplemented & eni ) { 326 isHuo_uoNotImplemented_ = true; 327 } 328 } 329 } 330 331 // Assemble the PDE adjoint Hessian. assembleHuo_un(const ROL::Vector<Real> & w,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const332 void assembleHuo_un(const ROL::Vector<Real> &w, 333 const ROL::Vector<Real> &uo, 334 const ROL::Vector<Real> &un, 335 const ROL::Vector<Real> &z, 336 const ROL::TimeStamp<Real> &ts) const { 337 if ( !isHuo_unZero_ && !isHuo_unNotImplemented_ ) { 338 try { 339 if (!isHuo_uoAssembled_) { 340 ROL::Ptr<const Tpetra::MultiVector<>> wf = getConstField(w); 341 ROL::Ptr<const Tpetra::MultiVector<>> uof = getConstField(uo); 342 ROL::Ptr<const Tpetra::MultiVector<>> unf = getConstField(un); 343 ROL::Ptr<const Tpetra::MultiVector<>> zf = getConstField(z); 344 ROL::Ptr<const std::vector<Real>> zp = getConstParameter(z); 345 346 assembler_->assembleDynPDEHessian_uo_un(matHuo_uo_,pde_,ts,wf,uof,unf,zf,zp); 347 isHuo_unAssembled_ = true; 348 } 349 } 350 catch ( Exception::Zero & ez ) { 351 isHuo_unZero_ = true; 352 } 353 catch ( Exception::NotImplemented & eni ) { 354 isHuo_unNotImplemented_ = true; 355 } 356 } 357 } 358 359 // Assemble the PDE adjoint Hessian. assembleHuo_zf(const ROL::Vector<Real> & w,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const360 void assembleHuo_zf(const ROL::Vector<Real> &w, 361 const ROL::Vector<Real> &uo, 362 const ROL::Vector<Real> &un, 363 const ROL::Vector<Real> &z, 364 const ROL::TimeStamp<Real> &ts) const { 365 if ( !isHuo_zfZero_ && !isHuo_zfNotImplemented_ ) { 366 try { 367 if (!isHuo_zfAssembled_) { 368 ROL::Ptr<const Tpetra::MultiVector<>> wf = getConstField(w); 369 ROL::Ptr<const Tpetra::MultiVector<>> uof = getConstField(uo); 370 ROL::Ptr<const Tpetra::MultiVector<>> unf = getConstField(un); 371 ROL::Ptr<const Tpetra::MultiVector<>> zf = getConstField(z); 372 ROL::Ptr<const std::vector<Real>> zp = getConstParameter(z); 373 374 assembler_->assembleDynPDEHessian_uo_zf(matHuo_zf_,pde_,ts,wf,uof,unf,zf,zp); 375 isHuo_zfAssembled_ = true; 376 } 377 } 378 catch ( Exception::Zero & ez ) { 379 isHuo_zfZero_ = true; 380 } 381 catch ( Exception::NotImplemented & eni ) { 382 isHuo_zfNotImplemented_ = true; 383 } 384 } 385 } 386 387 // Assemble the PDE adjoint Hessian. assembleHzf_uo(const ROL::Vector<Real> & w,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const388 void assembleHzf_uo(const ROL::Vector<Real> &w, 389 const ROL::Vector<Real> &uo, 390 const ROL::Vector<Real> &un, 391 const ROL::Vector<Real> &z, 392 const ROL::TimeStamp<Real> &ts) const { 393 if ( !isHzf_uoZero_ && !isHzf_uoNotImplemented_ ) { 394 try { 395 if (!isHzf_uoAssembled_) { 396 ROL::Ptr<const Tpetra::MultiVector<>> wf = getConstField(w); 397 ROL::Ptr<const Tpetra::MultiVector<>> uof = getConstField(uo); 398 ROL::Ptr<const Tpetra::MultiVector<>> unf = getConstField(un); 399 ROL::Ptr<const Tpetra::MultiVector<>> zf = getConstField(z); 400 ROL::Ptr<const std::vector<Real>> zp = getConstParameter(z); 401 402 assembler_->assembleDynPDEHessian_zf_uo(matHzf_uo_,pde_,ts,wf,uof,unf,zf,zp); 403 isHzf_uoAssembled_ = true; 404 } 405 } 406 catch ( Exception::Zero & ez ) { 407 isHzf_uoZero_ = true; 408 } 409 catch ( Exception::NotImplemented & eni ) { 410 isHzf_uoNotImplemented_ = true; 411 } 412 } 413 } 414 415 // Assemble the PDE adjoint Hessian. assembleHuo_zp(const ROL::Vector<Real> & w,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const416 void assembleHuo_zp(const ROL::Vector<Real> &w, 417 const ROL::Vector<Real> &uo, 418 const ROL::Vector<Real> &un, 419 const ROL::Vector<Real> &z, 420 const ROL::TimeStamp<Real> &ts) const { 421 if ( !isHuo_zpZero_ && !isHuo_zpNotImplemented_ ) { 422 try { 423 if (!isHuo_zpAssembled_) { 424 ROL::Ptr<const Tpetra::MultiVector<>> wf = getConstField(w); 425 ROL::Ptr<const Tpetra::MultiVector<>> uof = getConstField(uo); 426 ROL::Ptr<const Tpetra::MultiVector<>> unf = getConstField(un); 427 ROL::Ptr<const Tpetra::MultiVector<>> zf = getConstField(z); 428 ROL::Ptr<const std::vector<Real>> zp = getConstParameter(z); 429 430 assembler_->assembleDynPDEHessian_uo_zp(vecHuo_zp_,pde_,ts,wf,uof,unf,zf,zp); 431 isHuo_zpAssembled_ = true; 432 } 433 } 434 catch ( Exception::Zero & ez ) { 435 isHuo_zpZero_ = true; 436 } 437 catch ( Exception::NotImplemented & eni ) { 438 isHuo_zpNotImplemented_ = true; 439 } 440 } 441 } 442 443 // Assemble the PDE adjoint Hessian. assembleHzp_uo(const ROL::Vector<Real> & w,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const444 void assembleHzp_uo(const ROL::Vector<Real> &w, 445 const ROL::Vector<Real> &uo, 446 const ROL::Vector<Real> &un, 447 const ROL::Vector<Real> &z, 448 const ROL::TimeStamp<Real> &ts) const { 449 if ( !isHzp_uoZero_ && !isHzp_uoNotImplemented_ ) { 450 try { 451 if (!isHzp_uoAssembled_) { 452 ROL::Ptr<const Tpetra::MultiVector<>> wf = getConstField(w); 453 ROL::Ptr<const Tpetra::MultiVector<>> uof = getConstField(uo); 454 ROL::Ptr<const Tpetra::MultiVector<>> unf = getConstField(un); 455 ROL::Ptr<const Tpetra::MultiVector<>> zf = getConstField(z); 456 ROL::Ptr<const std::vector<Real>> zp = getConstParameter(z); 457 458 assembler_->assembleDynPDEHessian_zp_uo(vecHzp_uo_,pde_,ts,wf,uof,unf,zf,zp); 459 isHzp_uoAssembled_ = true; 460 } 461 } 462 catch ( Exception::Zero & ez ) { 463 isHzp_uoZero_ = true; 464 } 465 catch ( Exception::NotImplemented & eni ) { 466 isHzp_uoNotImplemented_ = true; 467 } 468 } 469 } 470 471 // Assemble the PDE adjoint Hessian. assembleHun_un(const ROL::Vector<Real> & w,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const472 void assembleHun_un(const ROL::Vector<Real> &w, 473 const ROL::Vector<Real> &uo, 474 const ROL::Vector<Real> &un, 475 const ROL::Vector<Real> &z, 476 const ROL::TimeStamp<Real> &ts) const { 477 if ( !isHun_unZero_ && !isHun_unNotImplemented_ ) { 478 try { 479 if (!isHun_unAssembled_) { 480 ROL::Ptr<const Tpetra::MultiVector<>> wf = getConstField(w); 481 ROL::Ptr<const Tpetra::MultiVector<>> uof = getConstField(uo); 482 ROL::Ptr<const Tpetra::MultiVector<>> unf = getConstField(un); 483 ROL::Ptr<const Tpetra::MultiVector<>> zf = getConstField(z); 484 ROL::Ptr<const std::vector<Real>> zp = getConstParameter(z); 485 486 assembler_->assembleDynPDEHessian_un_un(matHun_un_,pde_,ts,wf,uof,unf,zf,zp); 487 isHun_unAssembled_ = true; 488 } 489 } 490 catch ( Exception::Zero & ez ) { 491 isHun_unZero_ = true; 492 } 493 catch ( Exception::NotImplemented & eni ) { 494 isHun_unNotImplemented_ = true; 495 } 496 } 497 } 498 499 // Assemble the PDE adjoint Hessian. assembleHun_uo(const ROL::Vector<Real> & w,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const500 void assembleHun_uo(const ROL::Vector<Real> &w, 501 const ROL::Vector<Real> &uo, 502 const ROL::Vector<Real> &un, 503 const ROL::Vector<Real> &z, 504 const ROL::TimeStamp<Real> &ts) const { 505 if ( !isHun_uoZero_ && !isHun_uoNotImplemented_ ) { 506 try { 507 if (!isHun_uoAssembled_) { 508 ROL::Ptr<const Tpetra::MultiVector<>> wf = getConstField(w); 509 ROL::Ptr<const Tpetra::MultiVector<>> uof = getConstField(uo); 510 ROL::Ptr<const Tpetra::MultiVector<>> unf = getConstField(un); 511 ROL::Ptr<const Tpetra::MultiVector<>> zf = getConstField(z); 512 ROL::Ptr<const std::vector<Real>> zp = getConstParameter(z); 513 514 assembler_->assembleDynPDEHessian_un_uo(matHun_un_,pde_,ts,wf,uof,unf,zf,zp); 515 isHun_uoAssembled_ = true; 516 } 517 } 518 catch ( Exception::Zero & ez ) { 519 isHun_uoZero_ = true; 520 } 521 catch ( Exception::NotImplemented & eni ) { 522 isHun_uoNotImplemented_ = true; 523 } 524 } 525 } 526 527 // Assemble the PDE adjoint Hessian. assembleHun_zf(const ROL::Vector<Real> & w,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const528 void assembleHun_zf(const ROL::Vector<Real> &w, 529 const ROL::Vector<Real> &uo, 530 const ROL::Vector<Real> &un, 531 const ROL::Vector<Real> &z, 532 const ROL::TimeStamp<Real> &ts) const { 533 if ( !isHun_zfZero_ && !isHun_zfNotImplemented_ ) { 534 try { 535 if (!isHun_zfAssembled_) { 536 ROL::Ptr<const Tpetra::MultiVector<>> wf = getConstField(w); 537 ROL::Ptr<const Tpetra::MultiVector<>> uof = getConstField(uo); 538 ROL::Ptr<const Tpetra::MultiVector<>> unf = getConstField(un); 539 ROL::Ptr<const Tpetra::MultiVector<>> zf = getConstField(z); 540 ROL::Ptr<const std::vector<Real>> zp = getConstParameter(z); 541 542 assembler_->assembleDynPDEHessian_un_zf(matHun_zf_,pde_,ts,wf,uof,unf,zf,zp); 543 isHun_zfAssembled_ = true; 544 } 545 } 546 catch ( Exception::Zero & ez ) { 547 isHun_zfZero_ = true; 548 } 549 catch ( Exception::NotImplemented & eni ) { 550 isHun_zfNotImplemented_ = true; 551 } 552 } 553 } 554 555 // Assemble the PDE adjoint Hessian. assembleHzf_un(const ROL::Vector<Real> & w,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const556 void assembleHzf_un(const ROL::Vector<Real> &w, 557 const ROL::Vector<Real> &uo, 558 const ROL::Vector<Real> &un, 559 const ROL::Vector<Real> &z, 560 const ROL::TimeStamp<Real> &ts) const { 561 if ( !isHzf_unZero_ && !isHzf_unNotImplemented_ ) { 562 try { 563 if (!isHzf_unAssembled_) { 564 ROL::Ptr<const Tpetra::MultiVector<>> wf = getConstField(w); 565 ROL::Ptr<const Tpetra::MultiVector<>> uof = getConstField(uo); 566 ROL::Ptr<const Tpetra::MultiVector<>> unf = getConstField(un); 567 ROL::Ptr<const Tpetra::MultiVector<>> zf = getConstField(z); 568 ROL::Ptr<const std::vector<Real>> zp = getConstParameter(z); 569 570 assembler_->assembleDynPDEHessian_zf_un(matHzf_un_,pde_,ts,wf,uof,unf,zf,zp); 571 isHzf_unAssembled_ = true; 572 } 573 } 574 catch ( Exception::Zero & ez ) { 575 isHzf_unZero_ = true; 576 } 577 catch ( Exception::NotImplemented & eni ) { 578 isHzf_unNotImplemented_ = true; 579 } 580 } 581 } 582 583 // Assemble the PDE adjoint Hessian. assembleHun_zp(const ROL::Vector<Real> & w,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const584 void assembleHun_zp(const ROL::Vector<Real> &w, 585 const ROL::Vector<Real> &uo, 586 const ROL::Vector<Real> &un, 587 const ROL::Vector<Real> &z, 588 const ROL::TimeStamp<Real> &ts) const { 589 if ( !isHun_zpZero_ && !isHun_zpNotImplemented_ ) { 590 try { 591 if (!isHun_zpAssembled_) { 592 ROL::Ptr<const Tpetra::MultiVector<>> wf = getConstField(w); 593 ROL::Ptr<const Tpetra::MultiVector<>> uof = getConstField(uo); 594 ROL::Ptr<const Tpetra::MultiVector<>> unf = getConstField(un); 595 ROL::Ptr<const Tpetra::MultiVector<>> zf = getConstField(z); 596 ROL::Ptr<const std::vector<Real>> zp = getConstParameter(z); 597 598 assembler_->assembleDynPDEHessian_un_zp(vecHun_zp_,pde_,ts,wf,uof,unf,zf,zp); 599 isHun_zpAssembled_ = true; 600 } 601 } 602 catch ( Exception::Zero & ez ) { 603 isHun_zpZero_ = true; 604 } 605 catch ( Exception::NotImplemented & eni ) { 606 isHun_zpNotImplemented_ = true; 607 } 608 } 609 } 610 611 // Assemble the PDE adjoint Hessian. assembleHzp_un(const ROL::Vector<Real> & w,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const612 void assembleHzp_un(const ROL::Vector<Real> &w, 613 const ROL::Vector<Real> &uo, 614 const ROL::Vector<Real> &un, 615 const ROL::Vector<Real> &z, 616 const ROL::TimeStamp<Real> &ts) const { 617 if ( !isHzp_unZero_ && !isHzp_unNotImplemented_ ) { 618 try { 619 if (!isHzp_unAssembled_) { 620 ROL::Ptr<const Tpetra::MultiVector<>> wf = getConstField(w); 621 ROL::Ptr<const Tpetra::MultiVector<>> uof = getConstField(uo); 622 ROL::Ptr<const Tpetra::MultiVector<>> unf = getConstField(un); 623 ROL::Ptr<const Tpetra::MultiVector<>> zf = getConstField(z); 624 ROL::Ptr<const std::vector<Real>> zp = getConstParameter(z); 625 626 assembler_->assembleDynPDEHessian_zp_un(vecHzp_un_,pde_,ts,wf,uof,unf,zf,zp); 627 isHzp_unAssembled_ = true; 628 } 629 } 630 catch ( Exception::Zero & ez ) { 631 isHzp_unZero_ = true; 632 } 633 catch ( Exception::NotImplemented & eni ) { 634 isHzp_unNotImplemented_ = true; 635 } 636 } 637 } 638 639 // Assemble the PDE adjoint Hessian. assembleHzf_zf(const ROL::Vector<Real> & w,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const640 void assembleHzf_zf(const ROL::Vector<Real> &w, 641 const ROL::Vector<Real> &uo, 642 const ROL::Vector<Real> &un, 643 const ROL::Vector<Real> &z, 644 const ROL::TimeStamp<Real> &ts) const { 645 if ( !isHzf_zfZero_ && !isHzf_zfNotImplemented_ ) { 646 try { 647 if (!isHzf_zfAssembled_) { 648 ROL::Ptr<const Tpetra::MultiVector<>> wf = getConstField(w); 649 ROL::Ptr<const Tpetra::MultiVector<>> uof = getConstField(uo); 650 ROL::Ptr<const Tpetra::MultiVector<>> unf = getConstField(un); 651 ROL::Ptr<const Tpetra::MultiVector<>> zf = getConstField(z); 652 ROL::Ptr<const std::vector<Real>> zp = getConstParameter(z); 653 654 assembler_->assembleDynPDEHessian_zf_zf(matHzf_zf_,pde_,ts,wf,uof,unf,zf,zp); 655 isHzf_zfAssembled_ = true; 656 } 657 } 658 catch ( Exception::Zero & ez ) { 659 isHzf_zfZero_ = true; 660 } 661 catch ( Exception::NotImplemented & eni ) { 662 isHzf_zfNotImplemented_ = true; 663 } 664 } 665 } 666 667 // Assemble the PDE adjoint Hessian. assembleHzf_zp(const ROL::Vector<Real> & w,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const668 void assembleHzf_zp(const ROL::Vector<Real> &w, 669 const ROL::Vector<Real> &uo, 670 const ROL::Vector<Real> &un, 671 const ROL::Vector<Real> &z, 672 const ROL::TimeStamp<Real> &ts) const { 673 if ( !isHzf_zpZero_ && !isHzf_zpNotImplemented_ ) { 674 try { 675 if (!isHzf_zpAssembled_) { 676 ROL::Ptr<const Tpetra::MultiVector<>> wf = getConstField(w); 677 ROL::Ptr<const Tpetra::MultiVector<>> uof = getConstField(uo); 678 ROL::Ptr<const Tpetra::MultiVector<>> unf = getConstField(un); 679 ROL::Ptr<const Tpetra::MultiVector<>> zf = getConstField(z); 680 ROL::Ptr<const std::vector<Real>> zp = getConstParameter(z); 681 682 assembler_->assembleDynPDEHessian_zf_zp(vecHzf_zp_,pde_,ts,wf,uof,unf,zf,zp); 683 isHzf_zpAssembled_ = true; 684 } 685 } 686 catch ( Exception::Zero & ez ) { 687 isHzf_zpZero_ = true; 688 } 689 catch ( Exception::NotImplemented & eni ) { 690 isHzf_zpNotImplemented_ = true; 691 } 692 } 693 } 694 695 // Assemble the PDE adjoint Hessian. assembleHzp_zf(const ROL::Vector<Real> & w,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const696 void assembleHzp_zf(const ROL::Vector<Real> &w, 697 const ROL::Vector<Real> &uo, 698 const ROL::Vector<Real> &un, 699 const ROL::Vector<Real> &z, 700 const ROL::TimeStamp<Real> &ts) const { 701 if ( !isHzp_zfZero_ && !isHzp_zfNotImplemented_ ) { 702 try { 703 if (!isHzp_zfAssembled_) { 704 ROL::Ptr<const Tpetra::MultiVector<>> wf = getConstField(w); 705 ROL::Ptr<const Tpetra::MultiVector<>> uof = getConstField(uo); 706 ROL::Ptr<const Tpetra::MultiVector<>> unf = getConstField(un); 707 ROL::Ptr<const Tpetra::MultiVector<>> zf = getConstField(z); 708 ROL::Ptr<const std::vector<Real>> zp = getConstParameter(z); 709 710 assembler_->assembleDynPDEHessian_zp_zf(vecHzp_zf_,pde_,ts,wf,uof,unf,zf,zp); 711 isHzp_zfAssembled_ = true; 712 } 713 } 714 catch ( Exception::Zero & ez ) { 715 isHzp_zfZero_ = true; 716 } 717 catch ( Exception::NotImplemented & eni ) { 718 isHzp_zfNotImplemented_ = true; 719 } 720 } 721 } 722 723 // Assemble the PDE adjoint Hessian. assembleHzp_zp(const ROL::Vector<Real> & w,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const724 void assembleHzp_zp(const ROL::Vector<Real> &w, 725 const ROL::Vector<Real> &uo, 726 const ROL::Vector<Real> &un, 727 const ROL::Vector<Real> &z, 728 const ROL::TimeStamp<Real> &ts) const { 729 if ( !isHzp_zpZero_ && !isHzp_zpNotImplemented_ ) { 730 try { 731 if (!isHzp_zpAssembled_) { 732 ROL::Ptr<const Tpetra::MultiVector<>> wf = getConstField(w); 733 ROL::Ptr<const Tpetra::MultiVector<>> uof = getConstField(uo); 734 ROL::Ptr<const Tpetra::MultiVector<>> unf = getConstField(un); 735 ROL::Ptr<const Tpetra::MultiVector<>> zf = getConstField(z); 736 ROL::Ptr<const std::vector<Real>> zp = getConstParameter(z); 737 738 assembler_->assembleDynPDEHessian_zp_zp(matHzp_zp_,pde_,ts,wf,uof,unf,zf,zp); 739 isHzp_zpAssembled_ = true; 740 } 741 } 742 catch ( Exception::Zero & ez ) { 743 isHzp_zpZero_ = true; 744 } 745 catch ( Exception::NotImplemented & eni ) { 746 isHzp_zpNotImplemented_ = true; 747 } 748 } 749 } 750 solveForward(ROL::Ptr<Tpetra::MultiVector<>> & x,const ROL::Ptr<const Tpetra::MultiVector<>> & b) const751 void solveForward(ROL::Ptr<Tpetra::MultiVector<>> &x, 752 const ROL::Ptr<const Tpetra::MultiVector<>> &b) const { 753 #ifdef ROL_TIMERS 754 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintSolverSolve_Jacobian_un); 755 #endif 756 solver_->solve(x,b,false); 757 } 758 solveAdjoint(ROL::Ptr<Tpetra::MultiVector<>> & x,const ROL::Ptr<const Tpetra::MultiVector<>> & b) const759 void solveAdjoint(ROL::Ptr<Tpetra::MultiVector<>> &x, 760 const ROL::Ptr<const Tpetra::MultiVector<>> &b) const { 761 #ifdef ROL_TIMERS 762 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintSolverSolve_AdjointJacobian_un); 763 #endif 764 solver_->solve(x,b,true); 765 } 766 applyJacobian_uo(const ROL::Ptr<Tpetra::MultiVector<>> & Jv,const ROL::Ptr<const Tpetra::MultiVector<>> & v,const bool trans=false) const767 void applyJacobian_uo(const ROL::Ptr<Tpetra::MultiVector<>> &Jv, 768 const ROL::Ptr<const Tpetra::MultiVector<>> &v, 769 const bool trans = false) const { 770 if (!trans) { 771 #ifdef ROL_TIMERS 772 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyJacobian_uo); 773 #endif 774 matJuo_->apply(*v,*Jv,Teuchos::NO_TRANS); 775 } 776 else { 777 #ifdef ROL_TIMERS 778 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyAdjointJacobian_uo); 779 #endif 780 matJuo_->apply(*v,*Jv,Teuchos::TRANS); 781 } 782 } 783 applyJacobian_un(const ROL::Ptr<Tpetra::MultiVector<>> & Jv,const ROL::Ptr<const Tpetra::MultiVector<>> & v,const bool trans=false) const784 void applyJacobian_un(const ROL::Ptr<Tpetra::MultiVector<>> &Jv, 785 const ROL::Ptr<const Tpetra::MultiVector<>> &v, 786 const bool trans = false) const { 787 if (!trans) { 788 #ifdef ROL_TIMERS 789 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyJacobian_un); 790 #endif 791 matJun_->apply(*v,*Jv,Teuchos::NO_TRANS); 792 } 793 else { 794 #ifdef ROL_TIMERS 795 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyAdjointJacobian_un); 796 #endif 797 matJun_->apply(*v,*Jv,Teuchos::TRANS); 798 } 799 } 800 applyJacobian_zf(const ROL::Ptr<Tpetra::MultiVector<>> & Jv,const ROL::Ptr<const Tpetra::MultiVector<>> & v,const bool trans=false) const801 void applyJacobian_zf(const ROL::Ptr<Tpetra::MultiVector<>> &Jv, 802 const ROL::Ptr<const Tpetra::MultiVector<>> &v, 803 const bool trans = false) const { 804 if (!trans) { 805 #ifdef ROL_TIMERS 806 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyJacobian_zf); 807 #endif 808 matJzf_->apply(*v,*Jv,Teuchos::NO_TRANS); 809 } 810 else { 811 #ifdef ROL_TIMERS 812 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyAdjointJacobian_zf); 813 #endif 814 matJzf_->apply(*v,*Jv,Teuchos::TRANS); 815 } 816 } 817 applyJacobian_zp(const ROL::Ptr<Tpetra::MultiVector<>> & Jv,const ROL::Ptr<const std::vector<Real>> & v,const bool zeroOut=true) const818 void applyJacobian_zp(const ROL::Ptr<Tpetra::MultiVector<>> &Jv, 819 const ROL::Ptr<const std::vector<Real>> &v, 820 const bool zeroOut = true) const { 821 #ifdef ROL_TIMERS 822 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyJacobian_zp); 823 #endif 824 if ( zeroOut ) { 825 Jv->putScalar(static_cast<Real>(0)); 826 } 827 const size_t size = static_cast<size_t>(v->size()); 828 for (size_t i = 0; i < size; ++i) { 829 Teuchos::ArrayView<const size_t> col(&i,1); 830 Jv->update((*v)[i],*(vecJzp_->subView(col)),static_cast<Real>(1)); 831 } 832 } 833 applyAdjointJacobian_zp(const ROL::Ptr<std::vector<Real>> & Jv,const ROL::Ptr<const Tpetra::MultiVector<>> & v) const834 void applyAdjointJacobian_zp(const ROL::Ptr<std::vector<Real>> &Jv, 835 const ROL::Ptr<const Tpetra::MultiVector<>> &v) const { 836 #ifdef ROL_TIMERS 837 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyAdjointJacobian_zp); 838 #endif 839 Teuchos::Array<Real> val(1,0); 840 const size_t size = static_cast<size_t>(Jv->size()); 841 for (size_t i = 0; i < size; ++i) { 842 Teuchos::ArrayView<const size_t> col(&i,1); 843 vecJzp_->subView(col)->dot(*v, val.view(0,1)); 844 (*Jv)[i] = val[0]; 845 } 846 } 847 applyHessian_uo_uo(const ROL::Ptr<Tpetra::MultiVector<>> & Hv,const ROL::Ptr<const Tpetra::MultiVector<>> & v) const848 void applyHessian_uo_uo(const ROL::Ptr<Tpetra::MultiVector<>> &Hv, 849 const ROL::Ptr<const Tpetra::MultiVector<>> &v) const { 850 #ifdef ROL_TIMERS 851 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyHessian_uo_uo); 852 #endif 853 if (!isHuo_uoNotImplemented_) { 854 if ( isHuo_uoZero_ ) { 855 Hv->putScalar(static_cast<Real>(0)); 856 } 857 else { 858 matHuo_uo_->apply(*v,*Hv); 859 } 860 } 861 else { 862 throw ROL::Exception::NotImplemented(">>> DynConstraint<Real>::applyHessian_uo_uo: Not Implemented!"); 863 } 864 } 865 applyHessian_uo_un(const ROL::Ptr<Tpetra::MultiVector<>> & Hv,const ROL::Ptr<const Tpetra::MultiVector<>> & v) const866 void applyHessian_uo_un(const ROL::Ptr<Tpetra::MultiVector<>> &Hv, 867 const ROL::Ptr<const Tpetra::MultiVector<>> &v) const { 868 #ifdef ROL_TIMERS 869 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyHessian_uo_un); 870 #endif 871 if (!isHuo_unNotImplemented_) { 872 if ( isHuo_unZero_ ) { 873 Hv->putScalar(static_cast<Real>(0)); 874 } 875 else { 876 matHuo_un_->apply(*v,*Hv); 877 } 878 } 879 else { 880 throw ROL::Exception::NotImplemented(">>> DynConstraint<Real>::applyHessian_uo_un: Not Implemented!"); 881 } 882 } 883 applyHessian_un_uo(const ROL::Ptr<Tpetra::MultiVector<>> & Hv,const ROL::Ptr<const Tpetra::MultiVector<>> & v) const884 void applyHessian_un_uo(const ROL::Ptr<Tpetra::MultiVector<>> &Hv, 885 const ROL::Ptr<const Tpetra::MultiVector<>> &v) const { 886 #ifdef ROL_TIMERS 887 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyHessian_un_uo); 888 #endif 889 if (!isHun_uoNotImplemented_) { 890 if ( isHun_uoZero_ ) { 891 Hv->putScalar(static_cast<Real>(0)); 892 } 893 else { 894 matHun_uo_->apply(*v,*Hv); 895 } 896 } 897 else { 898 throw ROL::Exception::NotImplemented(">>> DynConstraint<Real>::applyHessian_un_uo: Not Implemented!"); 899 } 900 } 901 applyHessian_un_un(const ROL::Ptr<Tpetra::MultiVector<>> & Hv,const ROL::Ptr<const Tpetra::MultiVector<>> & v) const902 void applyHessian_un_un(const ROL::Ptr<Tpetra::MultiVector<>> &Hv, 903 const ROL::Ptr<const Tpetra::MultiVector<>> &v) const { 904 #ifdef ROL_TIMERS 905 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyHessian_un_un); 906 #endif 907 if (!isHun_unNotImplemented_) { 908 if ( isHun_unZero_ ) { 909 Hv->putScalar(static_cast<Real>(0)); 910 } 911 else { 912 matHun_un_->apply(*v,*Hv); 913 } 914 } 915 else { 916 throw ROL::Exception::NotImplemented(">>> DynConstraint<Real>::applyHessian_un_un: Not Implemented!"); 917 } 918 } 919 applyHessian_uo_zf(const ROL::Ptr<Tpetra::MultiVector<>> & Hv,const ROL::Ptr<const Tpetra::MultiVector<>> & v) const920 void applyHessian_uo_zf(const ROL::Ptr<Tpetra::MultiVector<>> &Hv, 921 const ROL::Ptr<const Tpetra::MultiVector<>> &v) const { 922 #ifdef ROL_TIMERS 923 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyHessian_uo_zf); 924 #endif 925 if ( Hv != ROL::nullPtr ) { 926 if ( !isHuo_zfNotImplemented_ ) { 927 if ( isHuo_zfZero_ ) { 928 Hv->putScalar(static_cast<Real>(0)); 929 } 930 else { 931 matHuo_zf_->apply(*v,*Hv); 932 } 933 } 934 else { 935 throw ROL::Exception::NotImplemented(">>> DynConstraint<Real>::applyHessian_uo_zf: Not Implemented!"); 936 } 937 } 938 } 939 applyHessian_un_zf(const ROL::Ptr<Tpetra::MultiVector<>> & Hv,const ROL::Ptr<const Tpetra::MultiVector<>> & v) const940 void applyHessian_un_zf(const ROL::Ptr<Tpetra::MultiVector<>> &Hv, 941 const ROL::Ptr<const Tpetra::MultiVector<>> &v) const { 942 #ifdef ROL_TIMERS 943 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyHessian_un_zf); 944 #endif 945 if ( Hv != ROL::nullPtr ) { 946 if ( !isHun_zfNotImplemented_ ) { 947 if ( isHun_zfZero_ ) { 948 Hv->putScalar(static_cast<Real>(0)); 949 } 950 else { 951 matHun_zf_->apply(*v,*Hv); 952 } 953 } 954 else { 955 throw ROL::Exception::NotImplemented(">>> DynConstraint<Real>::applyHessian_un_zf: Not Implemented!"); 956 } 957 } 958 } 959 applyHessian_uo_zp(const ROL::Ptr<std::vector<Real>> & Hv,const ROL::Ptr<const Tpetra::MultiVector<>> & v) const960 void applyHessian_uo_zp(const ROL::Ptr<std::vector<Real>> &Hv, 961 const ROL::Ptr<const Tpetra::MultiVector<>> &v) const { 962 #ifdef ROL_TIMERS 963 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyHessian_uo_zp); 964 #endif 965 if ( Hv != ROL::nullPtr ) { 966 if ( !isHuo_zpNotImplemented_ ) { 967 const size_t size = static_cast<size_t>(Hv->size()); 968 if ( isHuo_zpZero_ ) { 969 Hv->assign(size,static_cast<Real>(0)); 970 } 971 else { 972 Teuchos::Array<Real> val(1,0); 973 for (size_t i = 0; i < size; ++i) { 974 Teuchos::ArrayView<const size_t> col(&i,1); 975 vecHuo_zp_->subView(col)->dot(*v, val.view(0,1)); 976 (*Hv)[i] += val[0]; 977 } 978 } 979 } 980 else { 981 throw ROL::Exception::NotImplemented(">>> DynConstraint<Real>::applyHessian_uo_zp: Not Implemented!"); 982 } 983 } 984 } 985 applyHessian_un_zp(const ROL::Ptr<std::vector<Real>> & Hv,const ROL::Ptr<const Tpetra::MultiVector<>> & v) const986 void applyHessian_un_zp(const ROL::Ptr<std::vector<Real>> &Hv, 987 const ROL::Ptr<const Tpetra::MultiVector<>> &v) const { 988 #ifdef ROL_TIMERS 989 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyHessian_un_zp); 990 #endif 991 if ( Hv != ROL::nullPtr ) { 992 if ( !isHun_zpNotImplemented_ ) { 993 const size_t size = static_cast<size_t>(Hv->size()); 994 if ( isHun_zpZero_ ) { 995 Hv->assign(size,static_cast<Real>(0)); 996 } 997 else { 998 Teuchos::Array<Real> val(1,0); 999 for (size_t i = 0; i < size; ++i) { 1000 Teuchos::ArrayView<const size_t> col(&i,1); 1001 vecHun_zp_->subView(col)->dot(*v, val.view(0,1)); 1002 (*Hv)[i] += val[0]; 1003 } 1004 } 1005 } 1006 else { 1007 throw ROL::Exception::NotImplemented(">>> DynConstraint<Real>::applyHessian_un_zp: Not Implemented!"); 1008 } 1009 } 1010 } 1011 applyHessian_zf_uo(const ROL::Ptr<Tpetra::MultiVector<>> & Hv,const ROL::Ptr<const Tpetra::MultiVector<>> & v) const1012 void applyHessian_zf_uo(const ROL::Ptr<Tpetra::MultiVector<>> &Hv, 1013 const ROL::Ptr<const Tpetra::MultiVector<>> &v) const { 1014 #ifdef ROL_TIMERS 1015 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyHessian_zf_uo); 1016 #endif 1017 if ( v != ROL::nullPtr ) { 1018 if ( !isHzf_uoNotImplemented_ ) { 1019 if ( isHzf_uoZero_ ) { 1020 Hv->putScalar(static_cast<Real>(0)); 1021 } 1022 else { 1023 matHzf_uo_->apply(*v,*Hv); 1024 } 1025 } 1026 else { 1027 throw ROL::Exception::NotImplemented(">>> DynConstraint<Real>::applyHessian_zf_uo: Not Implemented!"); 1028 } 1029 } 1030 } 1031 applyHessian_zf_un(const ROL::Ptr<Tpetra::MultiVector<>> & Hv,const ROL::Ptr<const Tpetra::MultiVector<>> & v) const1032 void applyHessian_zf_un(const ROL::Ptr<Tpetra::MultiVector<>> &Hv, 1033 const ROL::Ptr<const Tpetra::MultiVector<>> &v) const { 1034 #ifdef ROL_TIMERS 1035 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyHessian_zf_un); 1036 #endif 1037 if ( v != ROL::nullPtr ) { 1038 if ( !isHzf_unNotImplemented_ ) { 1039 if ( isHzf_unZero_ ) { 1040 Hv->putScalar(static_cast<Real>(0)); 1041 } 1042 else { 1043 matHzf_un_->apply(*v,*Hv); 1044 } 1045 } 1046 else { 1047 throw ROL::Exception::NotImplemented(">>> DynConstraint<Real>::applyHessian_zf_un: Not Implemented!"); 1048 } 1049 } 1050 } 1051 applyHessian_zp_uo(const ROL::Ptr<Tpetra::MultiVector<>> & Hv,const ROL::Ptr<const std::vector<Real>> & v,const bool zeroOut=true) const1052 void applyHessian_zp_uo(const ROL::Ptr<Tpetra::MultiVector<>> &Hv, 1053 const ROL::Ptr<const std::vector<Real>> &v, 1054 const bool zeroOut = true) const { 1055 #ifdef ROL_TIMERS 1056 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyHessian_zp_uo); 1057 #endif 1058 if ( v != ROL::nullPtr ) { 1059 if ( !isHzp_uoNotImplemented_ ) { 1060 const size_t size = static_cast<size_t>(v->size()); 1061 if ( zeroOut || (isHzp_uoZero_ && !zeroOut) ) { 1062 Hv->putScalar(static_cast<Real>(0)); 1063 } 1064 if ( !isHzp_uoZero_ ) { 1065 for (size_t i = 0; i < size; ++i) { 1066 Teuchos::ArrayView<const size_t> col(&i,1); 1067 Hv->update((*v)[i],*(vecHzp_uo_->subView(col)),static_cast<Real>(1)); 1068 } 1069 } 1070 } 1071 else { 1072 throw ROL::Exception::NotImplemented(">>> DynConstraint<Real>::applyHessian_zp_uo: Not Implemented!"); 1073 } 1074 } 1075 } 1076 applyHessian_zp_un(const ROL::Ptr<Tpetra::MultiVector<>> & Hv,const ROL::Ptr<const std::vector<Real>> & v,const bool zeroOut=true) const1077 void applyHessian_zp_un(const ROL::Ptr<Tpetra::MultiVector<>> &Hv, 1078 const ROL::Ptr<const std::vector<Real>> &v, 1079 const bool zeroOut = true) const { 1080 #ifdef ROL_TIMERS 1081 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyHessian_zp_un); 1082 #endif 1083 if ( v != ROL::nullPtr ) { 1084 if ( !isHzp_unNotImplemented_ ) { 1085 const size_t size = static_cast<size_t>(v->size()); 1086 if ( zeroOut || (isHzp_unZero_ && !zeroOut) ) { 1087 Hv->putScalar(static_cast<Real>(0)); 1088 } 1089 if ( !isHzp_unZero_ ) { 1090 for (size_t i = 0; i < size; ++i) { 1091 Teuchos::ArrayView<const size_t> col(&i,1); 1092 Hv->update((*v)[i],*(vecHzp_un_->subView(col)),static_cast<Real>(1)); 1093 } 1094 } 1095 } 1096 else { 1097 throw ROL::Exception::NotImplemented(">>> DynConstraint<Real>::applyHessian_zp_un: Not Implemented!"); 1098 } 1099 } 1100 } 1101 applyHessian_zf_zf(const ROL::Ptr<Tpetra::MultiVector<>> & Hv,const ROL::Ptr<const Tpetra::MultiVector<>> & v) const1102 void applyHessian_zf_zf(const ROL::Ptr<Tpetra::MultiVector<>> &Hv, 1103 const ROL::Ptr<const Tpetra::MultiVector<>> &v) const { 1104 #ifdef ROL_TIMERS 1105 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyHessian_zf_zf); 1106 #endif 1107 if ( v != ROL::nullPtr ) { 1108 if ( !isHzf_zfNotImplemented_ ) { 1109 if ( isHzf_zfZero_ ) { 1110 Hv->putScalar(static_cast<Real>(0)); 1111 } 1112 else { 1113 matHzf_zf_->apply(*v,*Hv); 1114 } 1115 } 1116 else { 1117 throw ROL::Exception::NotImplemented(">>> DynConstraint<Real>::applyHessian_zf_zf: Not Implemented!"); 1118 } 1119 } 1120 } 1121 applyHessian_zf_zp(const ROL::Ptr<std::vector<Real>> & Hv,const ROL::Ptr<const Tpetra::MultiVector<>> & v,const bool zeroOut=true) const1122 void applyHessian_zf_zp(const ROL::Ptr<std::vector<Real>> &Hv, 1123 const ROL::Ptr<const Tpetra::MultiVector<>> &v, 1124 const bool zeroOut = true) const { 1125 #ifdef ROL_TIMERS 1126 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyHessian_zf_zp); 1127 #endif 1128 if ( Hv != ROL::nullPtr && v != ROL::nullPtr ) { 1129 if ( !isHzf_zpNotImplemented_ ) { 1130 const size_t size = static_cast<size_t>(Hv->size()); 1131 if ( zeroOut || (isHzf_zpZero_ && !zeroOut) ) { 1132 Hv->assign(size,static_cast<Real>(0)); 1133 } 1134 if ( !isHzf_zpZero_ ) { 1135 Teuchos::Array<Real> val(1,0); 1136 for (size_t i = 0; i < size; ++i) { 1137 Teuchos::ArrayView<const size_t> col(&i,1); 1138 vecHzf_zp_->subView(col)->dot(*v, val.view(0,1)); 1139 (*Hv)[i] += val[0]; 1140 } 1141 } 1142 } 1143 else { 1144 throw ROL::Exception::NotImplemented(">>> DynConstraint<Real>::applyHessian_zf_zp: Not Implemented!"); 1145 } 1146 } 1147 } 1148 applyHessian_zp_zf(const ROL::Ptr<Tpetra::MultiVector<>> & Hv,const ROL::Ptr<const std::vector<Real>> & v,const bool zeroOut=true) const1149 void applyHessian_zp_zf(const ROL::Ptr<Tpetra::MultiVector<>> &Hv, 1150 const ROL::Ptr<const std::vector<Real>> &v, 1151 const bool zeroOut = true) const { 1152 #ifdef ROL_TIMERS 1153 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyHessian_zp_zf); 1154 #endif 1155 if ( Hv != ROL::nullPtr && v != ROL::nullPtr ) { 1156 if ( !isHzp_zfNotImplemented_ ) { 1157 const size_t size = static_cast<size_t>(v->size()); 1158 if ( zeroOut || (isHzp_zfZero_ && !zeroOut) ) { 1159 Hv->putScalar(static_cast<Real>(0)); 1160 } 1161 if ( !isHzp_zfZero_ ) { 1162 for (size_t i = 0; i < size; ++i) { 1163 Teuchos::ArrayView<const size_t> col(&i,1); 1164 Hv->update((*v)[i],*(vecHzp_zf_->subView(col)),static_cast<Real>(1)); 1165 } 1166 } 1167 } 1168 else { 1169 throw ROL::Exception::NotImplemented(">>> DynConstraint<Real>::applyHessian_zp_zf: Not Implemented!"); 1170 } 1171 } 1172 } 1173 applyHessian_zp_zp(const ROL::Ptr<std::vector<Real>> & Hv,const ROL::Ptr<const std::vector<Real>> & v,const bool zeroOut=true) const1174 void applyHessian_zp_zp(const ROL::Ptr<std::vector<Real>> &Hv, 1175 const ROL::Ptr<const std::vector<Real>> &v, 1176 const bool zeroOut = true ) const { 1177 #ifdef ROL_TIMERS 1178 Teuchos::TimeMonitor LocalTimer(*ROL::PDEOPT::DynConstraintApplyHessian_zp_zp); 1179 #endif 1180 if ( Hv != ROL::nullPtr ) { 1181 if ( !isHzp_zpNotImplemented_ ) { 1182 const size_t size = static_cast<size_t>(Hv->size()); 1183 if ( zeroOut || (isHzp_zpZero_ && !zeroOut) ) { 1184 Hv->assign(size,static_cast<Real>(0)); 1185 } 1186 if ( !isHzp_zpZero_ ) { 1187 for (size_t i = 0; i < size; ++i) { 1188 for (size_t j = 0; j < size; ++j) { 1189 (*Hv)[i] += (*matHzp_zp_)[i][j]*(*v)[j]; 1190 } 1191 } 1192 } 1193 } 1194 else { 1195 throw ROL::Exception::NotImplemented(">>> DynConstraint<Real>::applyHessian_zp_zp: Not Implemented!"); 1196 } 1197 } 1198 } 1199 1200 public: 1201 DynConstraint(const ROL::Ptr<DynamicPDE<Real>> & pde,const ROL::Ptr<MeshManager<Real>> & meshMgr,const ROL::Ptr<const Teuchos::Comm<int>> & comm,Teuchos::ParameterList & parlist,std::ostream & outStream=std::cout)1202 DynConstraint(const ROL::Ptr<DynamicPDE<Real>> &pde, 1203 const ROL::Ptr<MeshManager<Real>> &meshMgr, 1204 const ROL::Ptr<const Teuchos::Comm<int>> &comm, 1205 Teuchos::ParameterList &parlist, 1206 std::ostream &outStream = std::cout) 1207 : pde_ ( pde ), 1208 isResAssembled_ ( false ), 1209 isJuoAssembled_ ( false ), 1210 isJuoZero_ ( false ), 1211 isJuoNotImplemented_ ( false ), 1212 isJunAssembled_ ( false ), 1213 isJunZero_ ( false ), 1214 isJunNotImplemented_ ( false ), 1215 isJzfAssembled_ ( false ), 1216 isJzfZero_ ( false ), 1217 isJzfNotImplemented_ ( false ), 1218 isJzpAssembled_ ( false ), 1219 isJzpZero_ ( false ), 1220 isJzpNotImplemented_ ( false ), 1221 isHuo_uoAssembled_ ( false ), 1222 isHuo_uoZero_ ( false ), 1223 isHuo_uoNotImplemented_ ( false ), 1224 isHuo_zfAssembled_ ( false ), 1225 isHuo_zfZero_ ( false ), 1226 isHuo_zfNotImplemented_ ( false ), 1227 isHzf_uoAssembled_ ( false ), 1228 isHzf_uoZero_ ( false ), 1229 isHzf_uoNotImplemented_ ( false ), 1230 isHuo_zpAssembled_ ( false ), 1231 isHuo_zpZero_ ( false ), 1232 isHuo_zpNotImplemented_ ( false ), 1233 isHzp_uoAssembled_ ( false ), 1234 isHzp_uoZero_ ( false ), 1235 isHzp_uoNotImplemented_ ( false ), 1236 isHuo_unAssembled_ ( false ), 1237 isHuo_unZero_ ( false ), 1238 isHuo_unNotImplemented_ ( false ), 1239 isHun_uoAssembled_ ( false ), 1240 isHun_uoZero_ ( false ), 1241 isHun_uoNotImplemented_ ( false ), 1242 isHun_unAssembled_ ( false ), 1243 isHun_unZero_ ( false ), 1244 isHun_unNotImplemented_ ( false ), 1245 isHun_zfAssembled_ ( false ), 1246 isHun_zfZero_ ( false ), 1247 isHun_zfNotImplemented_ ( false ), 1248 isHzf_unAssembled_ ( false ), 1249 isHzf_unZero_ ( false ), 1250 isHzf_unNotImplemented_ ( false ), 1251 isHun_zpAssembled_ ( false ), 1252 isHun_zpZero_ ( false ), 1253 isHun_zpNotImplemented_ ( false ), 1254 isHzp_unAssembled_ ( false ), 1255 isHzp_unZero_ ( false ), 1256 isHzp_unNotImplemented_ ( false ), 1257 isHzf_zfAssembled_ ( false ), 1258 isHzf_zfZero_ ( false ), 1259 isHzf_zfNotImplemented_ ( false ), 1260 isHzp_zpAssembled_ ( false ), 1261 isHzp_zpZero_ ( false ), 1262 isHzp_zpNotImplemented_ ( false ), 1263 isHzf_zpAssembled_ ( false ), 1264 isHzf_zpZero_ ( false ), 1265 isHzf_zpNotImplemented_ ( false ), 1266 isHzp_zfAssembled_ ( false ), 1267 isHzp_zfZero_ ( false ), 1268 isHzp_zfNotImplemented_ ( false ) { 1269 // Construct assembler. 1270 assembler_ = ROL::makePtr<Assembler<Real>>(pde_->getFields(),meshMgr,comm,parlist,outStream); 1271 assembler_->setCellNodes(*pde_); 1272 // Construct solver. 1273 solver_ = ROL::makePtr<Solver<Real>>(parlist.sublist("Solver")); 1274 // Initialize state and constraint vectors. 1275 cvec_ = assembler_->createResidualVector(); 1276 } 1277 getAssembler(void) const1278 const ROL::Ptr<Assembler<Real>> getAssembler(void) const { 1279 return assembler_; 1280 } 1281 getPDE(void) const1282 const ROL::Ptr<DynamicPDE<Real>> getPDE(void) const { 1283 return pde_; 1284 } 1285 update_uo(const ROL::Vector<Real> & uo,const ROL::TimeStamp<Real> & ts)1286 void update_uo(const ROL::Vector<Real> &uo, 1287 const ROL::TimeStamp<Real> &ts) { 1288 isResAssembled_ = false; 1289 isJuoAssembled_ = (isJuoZero_ ? isJuoAssembled_ : false); 1290 isJunAssembled_ = (isJunZero_ ? isJunAssembled_ : false); 1291 isJzfAssembled_ = (isJzfZero_ ? isJzfAssembled_ : false); 1292 isJzpAssembled_ = (isJzpZero_ ? isJzpAssembled_ : false); 1293 isHuo_uoAssembled_ = (isHuo_uoZero_ ? isHuo_uoAssembled_ : false); 1294 isHuo_zfAssembled_ = (isHuo_zfZero_ ? isHuo_zfAssembled_ : false); 1295 isHuo_zpAssembled_ = (isHuo_zpZero_ ? isHuo_zpAssembled_ : false); 1296 isHzf_uoAssembled_ = (isHzf_uoZero_ ? isHzf_uoAssembled_ : false); 1297 isHzp_uoAssembled_ = (isHzp_uoZero_ ? isHzp_uoAssembled_ : false); 1298 isHun_unAssembled_ = (isHun_unZero_ ? isHun_unAssembled_ : false); 1299 isHuo_unAssembled_ = (isHuo_unZero_ ? isHuo_unAssembled_ : false); 1300 isHun_uoAssembled_ = (isHun_uoZero_ ? isHun_uoAssembled_ : false); 1301 isHun_zfAssembled_ = (isHun_zfZero_ ? isHun_zfAssembled_ : false); 1302 isHun_zpAssembled_ = (isHun_zpZero_ ? isHun_zpAssembled_ : false); 1303 isHzf_unAssembled_ = (isHzf_unZero_ ? isHzf_unAssembled_ : false); 1304 isHzp_unAssembled_ = (isHzp_unZero_ ? isHzp_unAssembled_ : false); 1305 isHzf_zfAssembled_ = (isHzf_zfZero_ ? isHzf_zfAssembled_ : false); 1306 isHzf_zpAssembled_ = (isHzf_zpZero_ ? isHzf_zpAssembled_ : false); 1307 isHzp_zfAssembled_ = (isHzp_zfZero_ ? isHzp_zfAssembled_ : false); 1308 isHzp_zpAssembled_ = (isHzp_zpZero_ ? isHzp_zpAssembled_ : false); 1309 } 1310 update_un(const ROL::Vector<Real> & un,const ROL::TimeStamp<Real> & ts)1311 void update_un(const ROL::Vector<Real> &un, 1312 const ROL::TimeStamp<Real> &ts) { 1313 isResAssembled_ = false; 1314 isJuoAssembled_ = (isJuoZero_ ? isJuoAssembled_ : false); 1315 isJunAssembled_ = (isJunZero_ ? isJunAssembled_ : false); 1316 isJzfAssembled_ = (isJzfZero_ ? isJzfAssembled_ : false); 1317 isJzpAssembled_ = (isJzpZero_ ? isJzpAssembled_ : false); 1318 isHuo_uoAssembled_ = (isHuo_uoZero_ ? isHuo_uoAssembled_ : false); 1319 isHuo_zfAssembled_ = (isHuo_zfZero_ ? isHuo_zfAssembled_ : false); 1320 isHuo_zpAssembled_ = (isHuo_zpZero_ ? isHuo_zpAssembled_ : false); 1321 isHzf_uoAssembled_ = (isHzf_uoZero_ ? isHzf_uoAssembled_ : false); 1322 isHzp_uoAssembled_ = (isHzp_uoZero_ ? isHzp_uoAssembled_ : false); 1323 isHun_unAssembled_ = (isHun_unZero_ ? isHun_unAssembled_ : false); 1324 isHuo_unAssembled_ = (isHuo_unZero_ ? isHuo_unAssembled_ : false); 1325 isHun_uoAssembled_ = (isHun_uoZero_ ? isHun_uoAssembled_ : false); 1326 isHun_zfAssembled_ = (isHun_zfZero_ ? isHun_zfAssembled_ : false); 1327 isHun_zpAssembled_ = (isHun_zpZero_ ? isHun_zpAssembled_ : false); 1328 isHzf_unAssembled_ = (isHzf_unZero_ ? isHzf_unAssembled_ : false); 1329 isHzp_unAssembled_ = (isHzp_unZero_ ? isHzp_unAssembled_ : false); 1330 isHzf_zfAssembled_ = (isHzf_zfZero_ ? isHzf_zfAssembled_ : false); 1331 isHzf_zpAssembled_ = (isHzf_zpZero_ ? isHzf_zpAssembled_ : false); 1332 isHzp_zfAssembled_ = (isHzp_zfZero_ ? isHzp_zfAssembled_ : false); 1333 isHzp_zpAssembled_ = (isHzp_zpZero_ ? isHzp_zpAssembled_ : false); 1334 } 1335 update_z(const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts)1336 void update_z(const ROL::Vector<Real> &z, 1337 const ROL::TimeStamp<Real> &ts) { 1338 isResAssembled_ = false; 1339 isJuoAssembled_ = (isJuoZero_ ? isJuoAssembled_ : false); 1340 isJunAssembled_ = (isJunZero_ ? isJunAssembled_ : false); 1341 isJzfAssembled_ = (isJzfZero_ ? isJzfAssembled_ : false); 1342 isJzpAssembled_ = (isJzpZero_ ? isJzpAssembled_ : false); 1343 isHuo_uoAssembled_ = (isHuo_uoZero_ ? isHuo_uoAssembled_ : false); 1344 isHuo_zfAssembled_ = (isHuo_zfZero_ ? isHuo_zfAssembled_ : false); 1345 isHuo_zpAssembled_ = (isHuo_zpZero_ ? isHuo_zpAssembled_ : false); 1346 isHzf_uoAssembled_ = (isHzf_uoZero_ ? isHzf_uoAssembled_ : false); 1347 isHzp_uoAssembled_ = (isHzp_uoZero_ ? isHzp_uoAssembled_ : false); 1348 isHun_unAssembled_ = (isHun_unZero_ ? isHun_unAssembled_ : false); 1349 isHuo_unAssembled_ = (isHuo_unZero_ ? isHuo_unAssembled_ : false); 1350 isHun_uoAssembled_ = (isHun_uoZero_ ? isHun_uoAssembled_ : false); 1351 isHun_zfAssembled_ = (isHun_zfZero_ ? isHun_zfAssembled_ : false); 1352 isHun_zpAssembled_ = (isHun_zpZero_ ? isHun_zpAssembled_ : false); 1353 isHzf_unAssembled_ = (isHzf_unZero_ ? isHzf_unAssembled_ : false); 1354 isHzp_unAssembled_ = (isHzp_unZero_ ? isHzp_unAssembled_ : false); 1355 isHzf_zfAssembled_ = (isHzf_zfZero_ ? isHzf_zfAssembled_ : false); 1356 isHzf_zpAssembled_ = (isHzf_zpZero_ ? isHzf_zpAssembled_ : false); 1357 isHzp_zfAssembled_ = (isHzp_zfZero_ ? isHzp_zfAssembled_ : false); 1358 isHzp_zpAssembled_ = (isHzp_zpZero_ ? isHzp_zpAssembled_ : false); 1359 } 1360 update(const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts)1361 void update(const ROL::Vector<Real> &uo, 1362 const ROL::Vector<Real> &un, 1363 const ROL::Vector<Real> &z, 1364 const ROL::TimeStamp<Real> &ts) { 1365 update_uo(uo,ts); 1366 update_un(un,ts); 1367 update_z(z,ts); 1368 } 1369 value(ROL::Vector<Real> & c,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const1370 void value(ROL::Vector<Real> &c, 1371 const ROL::Vector<Real> &uo, 1372 const ROL::Vector<Real> &un, 1373 const ROL::Vector<Real> &z, 1374 const ROL::TimeStamp<Real> &ts) const { 1375 const Real zero(0), one(1); 1376 ROL::Ptr<Tpetra::MultiVector<>> cf = getField(c); 1377 assembleResidual(uo,un,z,ts); 1378 cf->update(one, *vecR_, zero); 1379 } 1380 solve(ROL::Vector<Real> & c,const ROL::Vector<Real> & uo,ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts)1381 void solve(ROL::Vector<Real> &c, 1382 const ROL::Vector<Real> &uo, 1383 ROL::Vector<Real> &un, 1384 const ROL::Vector<Real> &z, 1385 const ROL::TimeStamp<Real> &ts) { 1386 ROL::DynamicConstraint<Real>::solve(c,uo,un,z,ts); 1387 } 1388 applyJacobian_uo(ROL::Vector<Real> & jv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const1389 void applyJacobian_uo(ROL::Vector<Real> &jv, 1390 const ROL::Vector<Real> &v, 1391 const ROL::Vector<Real> &uo, 1392 const ROL::Vector<Real> &un, 1393 const ROL::Vector<Real> &z, 1394 const ROL::TimeStamp<Real> &ts) const { 1395 ROL::Ptr<Tpetra::MultiVector<>> jvf = getField(jv); 1396 ROL::Ptr<const Tpetra::MultiVector<>> vf = getConstField(v); 1397 assembleJuo(uo,un,z,ts); 1398 applyJacobian_uo(jvf,vf,false); 1399 } 1400 applyJacobian_un(ROL::Vector<Real> & jv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const1401 void applyJacobian_un(ROL::Vector<Real> &jv, 1402 const ROL::Vector<Real> &v, 1403 const ROL::Vector<Real> &uo, 1404 const ROL::Vector<Real> &un, 1405 const ROL::Vector<Real> &z, 1406 const ROL::TimeStamp<Real> &ts) const { 1407 ROL::Ptr<Tpetra::MultiVector<>> jvf = getField(jv); 1408 ROL::Ptr<const Tpetra::MultiVector<>> vf = getConstField(v); 1409 assembleJun(uo,un,z,ts); 1410 applyJacobian_un(jvf,vf,false); 1411 } 1412 1413 applyJacobian_z(ROL::Vector<Real> & jv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const1414 void applyJacobian_z(ROL::Vector<Real> &jv, 1415 const ROL::Vector<Real> &v, 1416 const ROL::Vector<Real> &uo, 1417 const ROL::Vector<Real> &un, 1418 const ROL::Vector<Real> &z, 1419 const ROL::TimeStamp<Real> &ts) const { 1420 jv.zero(); 1421 ROL::Ptr<Tpetra::MultiVector<>> jvf = getField(jv); 1422 ROL::Ptr<const Tpetra::MultiVector<>> vf = getConstField(v); 1423 ROL::Ptr<const std::vector<Real>> vp = getConstParameter(v); 1424 if (vf != ROL::nullPtr) { 1425 assembleJzf(uo,un,z,ts); 1426 applyJacobian_zf(jvf,vf,false); 1427 } 1428 if (vp != ROL::nullPtr) { 1429 assembleJzp(uo,un,z,ts); 1430 applyJacobian_zp(jvf,vp); 1431 } 1432 } 1433 1434 applyAdjointJacobian_uo(ROL::Vector<Real> & ajv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const1435 void applyAdjointJacobian_uo(ROL::Vector<Real> &ajv, 1436 const ROL::Vector<Real> &v, 1437 const ROL::Vector<Real> &uo, 1438 const ROL::Vector<Real> &un, 1439 const ROL::Vector<Real> &z, 1440 const ROL::TimeStamp<Real> &ts) const { 1441 ROL::Ptr<Tpetra::MultiVector<>> ajvf = getField(ajv); 1442 ROL::Ptr<const Tpetra::MultiVector<>> vf = getConstField(v); 1443 assembleJuo(uo,un,z,ts); 1444 applyJacobian_uo(ajvf,vf,true); 1445 } 1446 1447 applyAdjointJacobian_un(ROL::Vector<Real> & ajv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const1448 void applyAdjointJacobian_un(ROL::Vector<Real> &ajv, 1449 const ROL::Vector<Real> &v, 1450 const ROL::Vector<Real> &uo, 1451 const ROL::Vector<Real> &un, 1452 const ROL::Vector<Real> &z, 1453 const ROL::TimeStamp<Real> &ts) const { 1454 ROL::Ptr<Tpetra::MultiVector<>> ajvf = getField(ajv); 1455 ROL::Ptr<const Tpetra::MultiVector<>> vf = getConstField(v); 1456 assembleJun(uo,un,z,ts); 1457 applyJacobian_un(ajvf,vf,true); 1458 } 1459 1460 applyAdjointJacobian_z(ROL::Vector<Real> & ajv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const1461 void applyAdjointJacobian_z(ROL::Vector<Real> &ajv, 1462 const ROL::Vector<Real> &v, 1463 const ROL::Vector<Real> &uo, 1464 const ROL::Vector<Real> &un, 1465 const ROL::Vector<Real> &z, 1466 const ROL::TimeStamp<Real> &ts) const { 1467 ROL::Ptr<Tpetra::MultiVector<>> ajvf = getField(ajv); 1468 ROL::Ptr<std::vector<Real>> ajvp = getParameter(ajv); 1469 ROL::Ptr<const Tpetra::MultiVector<>> vf = getConstField(v); 1470 if (ajvf != ROL::nullPtr) { 1471 assembleJzf(uo,un,z,ts); 1472 applyJacobian_zf(ajvf,vf,true); 1473 } 1474 if (ajvp != ROL::nullPtr) { 1475 assembleJzp(uo,un,z,ts); 1476 applyAdjointJacobian_zp(ajvp,vf); 1477 } 1478 } 1479 1480 applyInverseJacobian_un(ROL::Vector<Real> & ijv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const1481 void applyInverseJacobian_un(ROL::Vector<Real> &ijv, 1482 const ROL::Vector<Real> &v, 1483 const ROL::Vector<Real> &uo, 1484 const ROL::Vector<Real> &un, 1485 const ROL::Vector<Real> &z, 1486 const ROL::TimeStamp<Real> &ts) const { 1487 ROL::Ptr<Tpetra::MultiVector<>> ijvf = getField(ijv); 1488 ROL::Ptr<const Tpetra::MultiVector<>> vf = getConstField(v); 1489 assembleJun(uo,un,z,ts); 1490 solveForward(ijvf,vf); 1491 } 1492 1493 applyInverseAdjointJacobian_un(ROL::Vector<Real> & iajv,const ROL::Vector<Real> & v,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const1494 void applyInverseAdjointJacobian_un(ROL::Vector<Real> &iajv, 1495 const ROL::Vector<Real> &v, 1496 const ROL::Vector<Real> &uo, 1497 const ROL::Vector<Real> &un, 1498 const ROL::Vector<Real> &z, 1499 const ROL::TimeStamp<Real> &ts) const { 1500 ROL::Ptr<Tpetra::MultiVector<>> iajvf = getField(iajv); 1501 ROL::Ptr<const Tpetra::MultiVector<>> vf = getConstField(v); 1502 assembleJun(uo,un,z,ts); 1503 solveAdjoint(iajvf,vf); 1504 } 1505 1506 applyAdjointHessian_uo_uo(ROL::Vector<Real> & ahwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const1507 void applyAdjointHessian_uo_uo(ROL::Vector<Real> &ahwv, 1508 const ROL::Vector<Real> &w, 1509 const ROL::Vector<Real> &v, 1510 const ROL::Vector<Real> &uo, 1511 const ROL::Vector<Real> &un, 1512 const ROL::Vector<Real> &z, 1513 const ROL::TimeStamp<Real> &ts) const { 1514 ROL::Ptr<Tpetra::MultiVector<>> ahwvf = getField(ahwv); 1515 ROL::Ptr<const Tpetra::MultiVector<>> vf = getConstField(v); 1516 assembleHuo_uo(w,uo,un,z,ts); 1517 applyHessian_uo_uo(ahwvf,vf); 1518 } 1519 1520 applyAdjointHessian_uo_un(ROL::Vector<Real> & ahwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const1521 void applyAdjointHessian_uo_un(ROL::Vector<Real> &ahwv, 1522 const ROL::Vector<Real> &w, 1523 const ROL::Vector<Real> &v, 1524 const ROL::Vector<Real> &uo, 1525 const ROL::Vector<Real> &un, 1526 const ROL::Vector<Real> &z, 1527 const ROL::TimeStamp<Real> &ts) const { 1528 ROL::Ptr<Tpetra::MultiVector<>> ahwvf = getField(ahwv); 1529 ROL::Ptr<const Tpetra::MultiVector<>> vf = getConstField(v); 1530 assembleHuo_un(w,uo,un,z,ts); 1531 applyHessian_uo_un(ahwvf,vf); 1532 } 1533 1534 applyAdjointHessian_uo_z(ROL::Vector<Real> & ahwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const1535 void applyAdjointHessian_uo_z(ROL::Vector<Real> &ahwv, 1536 const ROL::Vector<Real> &w, 1537 const ROL::Vector<Real> &v, 1538 const ROL::Vector<Real> &uo, 1539 const ROL::Vector<Real> &un, 1540 const ROL::Vector<Real> &z, 1541 const ROL::TimeStamp<Real> &ts) const { 1542 ROL::Ptr<std::vector<Real>> ahwvp = getParameter(ahwv); 1543 ROL::Ptr<Tpetra::MultiVector<>> ahwvf = getField(ahwv); 1544 ROL::Ptr<const Tpetra::MultiVector<>> vf = getConstField(v); 1545 if (ahwvf != ROL::nullPtr) { 1546 assembleHuo_zf(w,uo,un,z,ts); 1547 applyHessian_uo_zf(ahwvf,vf); 1548 } 1549 if (ahwvp != ROL::nullPtr) { 1550 assembleHuo_zp(w,uo,un,z,ts); 1551 applyHessian_uo_zp(ahwvp,vf); 1552 } 1553 } 1554 1555 applyAdjointHessian_un_uo(ROL::Vector<Real> & ahwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const1556 void applyAdjointHessian_un_uo(ROL::Vector<Real> &ahwv, 1557 const ROL::Vector<Real> &w, 1558 const ROL::Vector<Real> &v, 1559 const ROL::Vector<Real> &uo, 1560 const ROL::Vector<Real> &un, 1561 const ROL::Vector<Real> &z, 1562 const ROL::TimeStamp<Real> &ts) const { 1563 assembleHun_uo(w,uo,un,z,ts); 1564 ROL::Ptr<Tpetra::MultiVector<>> ahwvf = getField(ahwv); 1565 ROL::Ptr<const Tpetra::MultiVector<>> vf = getConstField(v); 1566 applyHessian_un_uo(ahwvf,vf); 1567 } 1568 1569 applyAdjointHessian_un_un(ROL::Vector<Real> & ahwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const1570 void applyAdjointHessian_un_un(ROL::Vector<Real> &ahwv, 1571 const ROL::Vector<Real> &w, 1572 const ROL::Vector<Real> &v, 1573 const ROL::Vector<Real> &uo, 1574 const ROL::Vector<Real> &un, 1575 const ROL::Vector<Real> &z, 1576 const ROL::TimeStamp<Real> &ts) const { 1577 assembleHun_un(w,uo,un,z,ts); 1578 ROL::Ptr<Tpetra::MultiVector<>> ahwvf = getField(ahwv); 1579 ROL::Ptr<const Tpetra::MultiVector<>> vf = getConstField(v); 1580 applyHessian_un_un(ahwvf,vf); 1581 } 1582 1583 applyAdjointHessian_un_z(ROL::Vector<Real> & ahwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const1584 void applyAdjointHessian_un_z(ROL::Vector<Real> &ahwv, 1585 const ROL::Vector<Real> &w, 1586 const ROL::Vector<Real> &v, 1587 const ROL::Vector<Real> &uo, 1588 const ROL::Vector<Real> &un, 1589 const ROL::Vector<Real> &z, 1590 const ROL::TimeStamp<Real> &ts) const { 1591 ROL::Ptr<std::vector<Real>> ahwvp = getParameter(ahwv); 1592 ROL::Ptr<Tpetra::MultiVector<>> ahwvf = getField(ahwv); 1593 ROL::Ptr<const Tpetra::MultiVector<>> vf = getConstField(v); 1594 if (ahwvf != ROL::nullPtr) { 1595 assembleHun_zf(w,uo,un,z,ts); 1596 applyHessian_un_zf(ahwvf,vf); 1597 } 1598 if (ahwvp != ROL::nullPtr) { 1599 assembleHun_zp(w,uo,un,z,ts); 1600 applyHessian_un_zp(ahwvp,vf); 1601 } 1602 } 1603 1604 applyAdjointHessian_z_uo(ROL::Vector<Real> & ahwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const1605 void applyAdjointHessian_z_uo(ROL::Vector<Real> &ahwv, 1606 const ROL::Vector<Real> &w, 1607 const ROL::Vector<Real> &v, 1608 const ROL::Vector<Real> &uo, 1609 const ROL::Vector<Real> &un, 1610 const ROL::Vector<Real> &z, 1611 const ROL::TimeStamp<Real> &ts) const { 1612 ROL::Ptr<Tpetra::MultiVector<>> ahwvf = getField(ahwv); 1613 ROL::Ptr<const Tpetra::MultiVector<>> vf = getConstField(v); 1614 ROL::Ptr<const std::vector<Real>> vp = getConstParameter(v); 1615 if (vf != ROL::nullPtr) { 1616 assembleHzf_uo(w,uo,un,z,ts); 1617 applyHessian_zf_uo(ahwvf,vf); 1618 } 1619 if (vp != ROL::nullPtr) { 1620 bool zeroOut = (vf == ROL::nullPtr); 1621 assembleHzp_uo(w,uo,un,z,ts); 1622 applyHessian_zp_uo(ahwvf,vp,zeroOut); 1623 } 1624 } 1625 1626 applyAdjointHessian_z_un(ROL::Vector<Real> & ahwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const1627 void applyAdjointHessian_z_un(ROL::Vector<Real> &ahwv, 1628 const ROL::Vector<Real> &w, 1629 const ROL::Vector<Real> &v, 1630 const ROL::Vector<Real> &uo, 1631 const ROL::Vector<Real> &un, 1632 const ROL::Vector<Real> &z, 1633 const ROL::TimeStamp<Real> &ts) const { 1634 ROL::Ptr<Tpetra::MultiVector<>> ahwvf = getField(ahwv); 1635 ROL::Ptr<const Tpetra::MultiVector<>> vf = getConstField(v); 1636 ROL::Ptr<const std::vector<Real>> vp = getConstParameter(v); 1637 if (vf != ROL::nullPtr) { 1638 assembleHzf_un(w,uo,un,z,ts); 1639 applyHessian_zf_un(ahwvf,vf); 1640 } 1641 if (vp != ROL::nullPtr) { 1642 bool zeroOut = (vf == ROL::nullPtr); 1643 assembleHzp_un(w,uo,un,z,ts); 1644 applyHessian_zp_un(ahwvf,vp,zeroOut); 1645 } 1646 } 1647 1648 applyAdjointHessian_z_z(ROL::Vector<Real> & ahwv,const ROL::Vector<Real> & w,const ROL::Vector<Real> & v,const ROL::Vector<Real> & uo,const ROL::Vector<Real> & un,const ROL::Vector<Real> & z,const ROL::TimeStamp<Real> & ts) const1649 void applyAdjointHessian_z_z(ROL::Vector<Real> &ahwv, 1650 const ROL::Vector<Real> &w, 1651 const ROL::Vector<Real> &v, 1652 const ROL::Vector<Real> &uo, 1653 const ROL::Vector<Real> &un, 1654 const ROL::Vector<Real> &z, 1655 const ROL::TimeStamp<Real> &ts) const { 1656 ROL::Ptr<Tpetra::MultiVector<>> ahwvf = getField(ahwv); 1657 ROL::Ptr<std::vector<Real>> ahwvp = getParameter(ahwv); 1658 ROL::Ptr<const Tpetra::MultiVector<>> vf = getConstField(v); 1659 ROL::Ptr<const std::vector<Real>> vp = getConstParameter(v); 1660 1661 bool zeroOut = (vf == ROL::nullPtr); 1662 if (vf != ROL::nullPtr) { 1663 assembleHzf_zf(w,uo,un,z,ts); 1664 applyHessian_zf_zf(ahwvf,vf); 1665 } 1666 if (vp != ROL::nullPtr) { 1667 assembleHzp_zp(w,uo,un,z,ts); 1668 applyHessian_zp_zp(ahwvp,vp,zeroOut); 1669 } 1670 if (vf != ROL::nullPtr && vp != ROL::nullPtr) { 1671 assembleHzf_zp(w,uo,un,z,ts); 1672 applyHessian_zf_zp(ahwvp,vf,zeroOut); 1673 assembleHzp_zf(w,uo,un,z,ts); 1674 applyHessian_zp_zf(ahwvf,vp); 1675 } 1676 } 1677 1678 /***************************************************************************/ 1679 /* Output routines. */ 1680 /***************************************************************************/ printMeshData(std::ostream & outStream) const1681 void printMeshData(std::ostream &outStream) const { 1682 assembler_->printMeshData(outStream); 1683 } 1684 outputTpetraData() const1685 void outputTpetraData() const { 1686 Tpetra::MatrixMarket::Writer< Tpetra::CrsMatrix<>> matWriter; 1687 if (matJuo_ != ROL::nullPtr) { 1688 matWriter.writeSparseFile("jacobian_uo.txt", matJuo_); 1689 } 1690 else { 1691 std::ofstream emptyfile; 1692 emptyfile.open("jacobian_uo.txt"); 1693 emptyfile.close(); 1694 } 1695 if (matJun_ != ROL::nullPtr) { 1696 matWriter.writeSparseFile("jacobian_un.txt", matJun_); 1697 } 1698 else { 1699 std::ofstream emptyfile; 1700 emptyfile.open("jacobian_un.txt"); 1701 emptyfile.close(); 1702 } 1703 if (matJzf_ != ROL::nullPtr) { 1704 matWriter.writeSparseFile("jacobian_zf.txt", matJzf_); 1705 } 1706 else { 1707 std::ofstream emptyfile; 1708 emptyfile.open("jacobian_zf.txt"); 1709 emptyfile.close(); 1710 } 1711 if (vecR_ != ROL::nullPtr) { 1712 matWriter.writeDenseFile("residual.txt", vecR_); 1713 } 1714 else { 1715 std::ofstream emptyfile; 1716 emptyfile.open("residual.txt"); 1717 emptyfile.close(); 1718 } 1719 } 1720 outputTpetraVector(const ROL::Ptr<const Tpetra::MultiVector<>> & vec,const std::string & filename) const1721 void outputTpetraVector(const ROL::Ptr<const Tpetra::MultiVector<>> &vec, 1722 const std::string &filename) const { 1723 assembler_->outputTpetraVector(vec, filename); 1724 } 1725 inputTpetraVector(ROL::Ptr<Tpetra::MultiVector<>> & vec,const std::string & filename) const1726 void inputTpetraVector(ROL::Ptr<Tpetra::MultiVector<>> &vec, 1727 const std::string &filename) const { 1728 assembler_->inputTpetraVector(vec, filename); 1729 } 1730 /***************************************************************************/ 1731 /* End of output routines. */ 1732 /***************************************************************************/ 1733 1734 private: // Vector accessor functions 1735 getConstField(const ROL::Vector<Real> & x) const1736 ROL::Ptr<const Tpetra::MultiVector<> > getConstField(const ROL::Vector<Real> &x) const { 1737 ROL::Ptr<const Tpetra::MultiVector<> > xp; 1738 try { 1739 xp = dynamic_cast<const ROL::TpetraMultiVector<Real>&>(x).getVector(); 1740 } 1741 catch (std::exception &e) { 1742 try { 1743 ROL::Ptr<const ROL::TpetraMultiVector<Real> > xvec 1744 = dynamic_cast<const PDE_OptVector<Real>&>(x).getField(); 1745 if (xvec == ROL::nullPtr) { 1746 xp = ROL::nullPtr; 1747 } 1748 else { 1749 xp = xvec->getVector(); 1750 } 1751 } 1752 catch (std::exception &ee) { 1753 xp = ROL::nullPtr; 1754 } 1755 } 1756 return xp; 1757 } 1758 getField(ROL::Vector<Real> & x) const1759 ROL::Ptr<Tpetra::MultiVector<> > getField(ROL::Vector<Real> &x) const { 1760 ROL::Ptr<Tpetra::MultiVector<> > xp; 1761 try { 1762 xp = dynamic_cast<ROL::TpetraMultiVector<Real>&>(x).getVector(); 1763 } 1764 catch (std::exception &e) { 1765 try { 1766 ROL::Ptr<ROL::TpetraMultiVector<Real> > xvec 1767 = dynamic_cast<PDE_OptVector<Real>&>(x).getField(); 1768 if ( xvec == ROL::nullPtr ) { 1769 xp = ROL::nullPtr; 1770 } 1771 else { 1772 xp = xvec->getVector(); 1773 } 1774 } 1775 catch (std::exception &ee) { 1776 xp = ROL::nullPtr; 1777 } 1778 } 1779 return xp; 1780 } 1781 getConstParameter(const ROL::Vector<Real> & x) const1782 ROL::Ptr<const std::vector<Real> > getConstParameter(const ROL::Vector<Real> &x) const { 1783 ROL::Ptr<const std::vector<Real> > xp; 1784 try { 1785 xp = dynamic_cast<const ROL::StdVector<Real>&>(x).getVector(); 1786 } 1787 catch (std::exception &e) { 1788 try { 1789 ROL::Ptr<const ROL::StdVector<Real> > xvec 1790 = dynamic_cast<const PDE_OptVector<Real>&>(x).getParameter(); 1791 if ( xvec == ROL::nullPtr ) { 1792 xp = ROL::nullPtr; 1793 } 1794 else { 1795 xp = xvec->getVector(); 1796 } 1797 } 1798 catch (std::exception &ee) { 1799 xp = ROL::nullPtr; 1800 } 1801 } 1802 return xp; 1803 } 1804 getParameter(ROL::Vector<Real> & x) const1805 ROL::Ptr<std::vector<Real> > getParameter(ROL::Vector<Real> &x) const { 1806 ROL::Ptr<std::vector<Real> > xp; 1807 try { 1808 xp = dynamic_cast<ROL::StdVector<Real>&>(x).getVector(); 1809 } 1810 catch (std::exception &e) { 1811 try { 1812 ROL::Ptr<ROL::StdVector<Real> > xvec 1813 = dynamic_cast<PDE_OptVector<Real>&>(x).getParameter(); 1814 if ( xvec == ROL::nullPtr ) { 1815 xp = ROL::nullPtr; 1816 } 1817 else { 1818 xp = xvec->getVector(); 1819 } 1820 } 1821 catch (std::exception &ee) { 1822 xp = ROL::nullPtr; 1823 } 1824 } 1825 return xp; 1826 } 1827 }; 1828 1829 #endif 1830