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 pde.hpp 45 \brief Implements the local PDE interface for the Poisson control problem. 46 */ 47 48 #ifndef PDEOPT_DYNAMIC_ADV_DIFF_HPP 49 #define PDEOPT_DYNAMIC_ADV_DIFF_HPP 50 51 #include "../../TOOLS/dynpde.hpp" 52 #include "../../TOOLS/fe.hpp" 53 54 #include "Intrepid_HGRAD_QUAD_C1_FEM.hpp" 55 #include "Intrepid_HGRAD_QUAD_C2_FEM.hpp" 56 #include "Intrepid_DefaultCubatureFactory.hpp" 57 #include "Intrepid_FunctionSpaceTools.hpp" 58 #include "Intrepid_CellTools.hpp" 59 60 #include "ROL_Ptr.hpp" 61 62 #include "pde_adv_diff.hpp" 63 64 template <class Real> 65 class DynamicPDE_adv_diff : public DynamicPDE<Real> { 66 private: 67 // Cell node information 68 std::vector<std::vector<std::vector<int>>> bdryCellLocIds_; 69 // Finite element definition 70 ROL::Ptr<FE<Real>> fe_vol_; 71 // Local degrees of freedom on boundary, for each side of the reference cell (first index). 72 std::vector<std::vector<int>> fidx_; 73 // Coordinates of degrees freedom on boundary cells. 74 // Indexing: [sideset number][local side id](cell number, value at dof) 75 std::vector<std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>>> bdryCellDofValues_; 76 77 // Steady PDE without Dirichlet BC 78 ROL::Ptr<PDE_adv_diff<Real>> pde_; 79 Real theta_; 80 Real useDBC_; 81 82 public: DynamicPDE_adv_diff(Teuchos::ParameterList & parlist)83 DynamicPDE_adv_diff(Teuchos::ParameterList &parlist) { 84 pde_ = ROL::makePtr<PDE_adv_diff<Real>>(parlist); 85 // Time-dependent coefficients 86 theta_ = parlist.sublist("Time Discretization").get("Theta",1.0); 87 useDBC_ = parlist.sublist("Problem").get("Use Dirichlet",false); 88 } 89 residual(ROL::Ptr<Intrepid::FieldContainer<Real>> & res,const ROL::TimeStamp<Real> & ts,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)90 void residual(ROL::Ptr<Intrepid::FieldContainer<Real>> & res, 91 const ROL::TimeStamp<Real> & ts, 92 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff, 93 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff, 94 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff = ROL::nullPtr, 95 const ROL::Ptr<const std::vector<Real>> & z_param = ROL::nullPtr) { 96 const Real one(1); 97 // GET DIMENSIONS 98 int c = fe_vol_->gradN()->dimension(0); 99 int f = fe_vol_->gradN()->dimension(1); 100 int p = fe_vol_->gradN()->dimension(2); 101 // GET TIME STEP INFORMATION 102 Real told = ts.t[0], tnew = ts.t[1], dt = tnew-told; 103 // INITIALIZE STORAGE 104 ROL::Ptr<Intrepid::FieldContainer<Real>> pde 105 = ROL::makePtr<Intrepid::FieldContainer<Real>>(c, f); 106 ROL::Ptr<Intrepid::FieldContainer<Real>> valU_eval 107 = ROL::makePtr<Intrepid::FieldContainer<Real>>(c, p); 108 // COMPUTE OLD RESIDUAL 109 pde_->setTime(told); 110 pde_->residual(res,uo_coeff,z_coeff,z_param); 111 // Integrate Uold * N 112 fe_vol_->evaluateValue(valU_eval, uo_coeff); 113 Intrepid::FunctionSpaceTools::integrate<Real>(*pde, 114 *valU_eval, 115 *(fe_vol_->NdetJ()), 116 Intrepid::COMP_CPP, false); 117 Intrepid::RealSpaceTools<Real>::scale(*res, (one-theta_)*dt); 118 Intrepid::RealSpaceTools<Real>::subtract(*res, *pde); 119 // COMPUTE NEW RESIDUAL 120 pde_->setTime(tnew); 121 pde_->residual(pde,un_coeff,z_coeff,z_param); 122 // Integrate Uold * N 123 valU_eval->initialize(); 124 fe_vol_->evaluateValue(valU_eval, un_coeff); 125 Intrepid::FunctionSpaceTools::integrate<Real>(*res, 126 *valU_eval, 127 *(fe_vol_->NdetJ()), 128 Intrepid::COMP_CPP, true); 129 Intrepid::RealSpaceTools<Real>::scale(*pde, theta_*dt); 130 Intrepid::RealSpaceTools<Real>::add(*res, *pde); 131 // APPLY DIRICHLET CONDITIONS 132 if (useDBC_) { 133 int numLocalSideIds = bdryCellLocIds_[0].size(); 134 for (int j = 0; j < numLocalSideIds; ++j) { 135 int numCellsSide = bdryCellLocIds_[0][j].size(); 136 int numBdryDofs = fidx_[j].size(); 137 for (int k = 0; k < numCellsSide; ++k) { 138 int cidx = bdryCellLocIds_[0][j][k]; 139 for (int l = 0; l < numBdryDofs; ++l) { 140 (*res)(cidx,fidx_[j][l]) 141 = (*un_coeff)(cidx,fidx_[j][l]) - (*bdryCellDofValues_[0][j])(k,fidx_[j][l]); 142 } 143 } 144 } 145 } 146 } 147 Jacobian_uo(ROL::Ptr<Intrepid::FieldContainer<Real>> & jac,const ROL::TimeStamp<Real> & ts,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)148 void Jacobian_uo(ROL::Ptr<Intrepid::FieldContainer<Real>> & jac, 149 const ROL::TimeStamp<Real> & ts, 150 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff, 151 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff, 152 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff = ROL::nullPtr, 153 const ROL::Ptr<const std::vector<Real>> & z_param = ROL::nullPtr) { 154 const Real one(1); 155 // GET DIMENSIONS 156 int c = fe_vol_->gradN()->dimension(0); 157 int f = fe_vol_->gradN()->dimension(1); 158 // GET TIME STEP INFORMATION 159 Real told = ts.t[0], tnew = ts.t[1], dt = tnew-told; 160 // INITILAIZE JACOBIAN 161 ROL::Ptr<Intrepid::FieldContainer<Real>> pde 162 = ROL::makePtr<Intrepid::FieldContainer<Real>>(c, f, f); 163 // COMPUTE OLD RESIDUAL 164 pde_->setTime(told); 165 pde_->Jacobian_1(jac,uo_coeff,z_coeff,z_param); 166 // Integrate N * N 167 Intrepid::FunctionSpaceTools::integrate<Real>(*pde, 168 *(fe_vol_->N()), 169 *(fe_vol_->NdetJ()), 170 Intrepid::COMP_CPP, false); 171 Intrepid::RealSpaceTools<Real>::scale(*jac, (one-theta_)*dt); 172 Intrepid::RealSpaceTools<Real>::subtract(*jac, *pde); 173 // APPLY DIRICHLET CONDITIONS 174 if (useDBC_) { 175 int numLocalSideIds = bdryCellLocIds_[0].size(); 176 for (int j = 0; j < numLocalSideIds; ++j) { 177 int numCellsSide = bdryCellLocIds_[0][j].size(); 178 int numBdryDofs = fidx_[j].size(); 179 for (int k = 0; k < numCellsSide; ++k) { 180 int cidx = bdryCellLocIds_[0][j][k]; 181 for (int l = 0; l < numBdryDofs; ++l) { 182 //std::cout << "\n j=" << j << " l=" << l << " " << fidx[j][l]; 183 for (int m = 0; m < f; ++m) { 184 (*jac)(cidx,fidx_[j][l],m) = static_cast<Real>(0); 185 } 186 } 187 } 188 } 189 } 190 } 191 Jacobian_un(ROL::Ptr<Intrepid::FieldContainer<Real>> & jac,const ROL::TimeStamp<Real> & ts,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)192 void Jacobian_un(ROL::Ptr<Intrepid::FieldContainer<Real>> & jac, 193 const ROL::TimeStamp<Real> & ts, 194 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff, 195 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff, 196 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff = ROL::nullPtr, 197 const ROL::Ptr<const std::vector<Real>> & z_param = ROL::nullPtr) { 198 // GET DIMENSIONS 199 int c = fe_vol_->gradN()->dimension(0); 200 int f = fe_vol_->gradN()->dimension(1); 201 // GET TIME STEP INFORMATION 202 Real told = ts.t[0], tnew = ts.t[1], dt = tnew-told; 203 // INITILAIZE JACOBIAN 204 jac = ROL::makePtr<Intrepid::FieldContainer<Real>>(c, f, f); 205 ROL::Ptr<Intrepid::FieldContainer<Real>> pde; 206 // COMPUTE NEW RESIDUAL 207 pde_->setTime(tnew); 208 pde_->Jacobian_1(pde,un_coeff,z_coeff,z_param); 209 // Integrate N * N 210 Intrepid::FunctionSpaceTools::integrate<Real>(*jac, 211 *(fe_vol_->N()), 212 *(fe_vol_->NdetJ()), 213 Intrepid::COMP_CPP, false); 214 Intrepid::RealSpaceTools<Real>::scale(*pde, theta_*dt); 215 Intrepid::RealSpaceTools<Real>::add(*jac, *pde); 216 // APPLY DIRICHLET CONDITIONS 217 if (useDBC_) { 218 int numLocalSideIds = bdryCellLocIds_[0].size(); 219 for (int j = 0; j < numLocalSideIds; ++j) { 220 int numCellsSide = bdryCellLocIds_[0][j].size(); 221 int numBdryDofs = fidx_[j].size(); 222 for (int k = 0; k < numCellsSide; ++k) { 223 int cidx = bdryCellLocIds_[0][j][k]; 224 for (int l = 0; l < numBdryDofs; ++l) { 225 //std::cout << "\n j=" << j << " l=" << l << " " << fidx[j][l]; 226 for (int m = 0; m < f; ++m) { 227 (*jac)(cidx,fidx_[j][l],m) = static_cast<Real>(0); 228 } 229 (*jac)(cidx,fidx_[j][l],fidx_[j][l]) = static_cast<Real>(1); 230 } 231 } 232 } 233 } 234 } 235 Jacobian_zf(ROL::Ptr<Intrepid::FieldContainer<Real>> & jac,const ROL::TimeStamp<Real> & ts,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)236 void Jacobian_zf(ROL::Ptr<Intrepid::FieldContainer<Real>> & jac, 237 const ROL::TimeStamp<Real> & ts, 238 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff, 239 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff, 240 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff = ROL::nullPtr, 241 const ROL::Ptr<const std::vector<Real>> & z_param = ROL::nullPtr) { 242 throw Exception::Zero(">>> (PDE_adv_diff::Jacobian_zf): Jacobian is zero."); 243 } 244 Jacobian_zp(std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>> & jac,const ROL::TimeStamp<Real> & ts,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)245 void Jacobian_zp(std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>> & jac, 246 const ROL::TimeStamp<Real> & ts, 247 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff, 248 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff, 249 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff = ROL::nullPtr, 250 const ROL::Ptr<const std::vector<Real>> & z_param = ROL::nullPtr) { 251 const Real one(1); 252 // GET TIME STEP INFORMATION 253 Real told = ts.t[0], tnew = ts.t[1], dt = tnew-told; 254 // ADD CONTROL TERM TO RESIDUAL 255 int size = z_param->size(); 256 std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>> pde(size,ROL::nullPtr); 257 // EVALUATE OLD COMPONENT 258 pde_->setTime(told); 259 pde_->Jacobian_3(jac,uo_coeff,z_coeff,z_param); 260 // EVALUATE NEW COMPONENT 261 pde_->setTime(tnew); 262 pde_->Jacobian_3(pde,un_coeff,z_coeff,z_param); 263 for (int i = 0; i < size; ++i) { 264 Intrepid::RealSpaceTools<Real>::scale(*(jac[i]), (one-theta_)*dt); 265 Intrepid::RealSpaceTools<Real>::scale(*(pde[i]), theta_*dt); 266 Intrepid::RealSpaceTools<Real>::add(*(jac[i]), *(pde[i])); 267 // APPLY DIRICHLET CONDITIONS 268 if (useDBC_) { 269 int numLocalSideIds = bdryCellLocIds_[0].size(); 270 for (int j = 0; j < numLocalSideIds; ++j) { 271 int numCellsSide = bdryCellLocIds_[0][j].size(); 272 int numBdryDofs = fidx_[j].size(); 273 for (int k = 0; k < numCellsSide; ++k) { 274 int cidx = bdryCellLocIds_[0][j][k]; 275 for (int l = 0; l < numBdryDofs; ++l) { 276 (*(jac[i]))(cidx,fidx_[j][l]) = static_cast<Real>(0); 277 } 278 } 279 } 280 } 281 } 282 } 283 Hessian_uo_uo(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess,const ROL::TimeStamp<Real> & ts,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)284 void Hessian_uo_uo(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess, 285 const ROL::TimeStamp<Real> & ts, 286 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff, 287 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff, 288 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff, 289 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff = ROL::nullPtr, 290 const ROL::Ptr<const std::vector<Real>> & z_param = ROL::nullPtr) { 291 throw Exception::Zero(">>> (PDE_adv_diff::Hessian_uo_uo): Hessian is zero."); 292 } 293 Hessian_uo_un(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess,const ROL::TimeStamp<Real> & ts,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)294 void Hessian_uo_un(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess, 295 const ROL::TimeStamp<Real> & ts, 296 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff, 297 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff, 298 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff, 299 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff = ROL::nullPtr, 300 const ROL::Ptr<const std::vector<Real>> & z_param = ROL::nullPtr) { 301 throw Exception::Zero(">>> (PDE_adv_diff::Hessian_uo_un): Hessian is zero."); 302 } 303 Hessian_uo_zf(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess,const ROL::TimeStamp<Real> & ts,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)304 void Hessian_uo_zf(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess, 305 const ROL::TimeStamp<Real> & ts, 306 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff, 307 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff, 308 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff, 309 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff = ROL::nullPtr, 310 const ROL::Ptr<const std::vector<Real>> & z_param = ROL::nullPtr) { 311 throw Exception::Zero(">>> (PDE_adv_diff::Hessian_uo_zf): Hessian is zero."); 312 } 313 Hessian_uo_zp(std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>> & hess,const ROL::TimeStamp<Real> & ts,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)314 void Hessian_uo_zp(std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>> & hess, 315 const ROL::TimeStamp<Real> & ts, 316 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff, 317 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff, 318 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff, 319 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff = ROL::nullPtr, 320 const ROL::Ptr<const std::vector<Real>> & z_param = ROL::nullPtr) { 321 throw Exception::Zero(">>> (PDE_adv_diff::Hessian_uo_zp): Hessian is zero."); 322 } 323 Hessian_un_uo(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess,const ROL::TimeStamp<Real> & ts,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)324 void Hessian_un_uo(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess, 325 const ROL::TimeStamp<Real> & ts, 326 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff, 327 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff, 328 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff, 329 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff = ROL::nullPtr, 330 const ROL::Ptr<const std::vector<Real>> & z_param = ROL::nullPtr) { 331 throw Exception::Zero(">>> (PDE_adv_diff::Hessian_un_uo): Hessian is zero."); 332 } 333 Hessian_un_un(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess,const ROL::TimeStamp<Real> & ts,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)334 void Hessian_un_un(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess, 335 const ROL::TimeStamp<Real> & ts, 336 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff, 337 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff, 338 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff, 339 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff = ROL::nullPtr, 340 const ROL::Ptr<const std::vector<Real>> & z_param = ROL::nullPtr) { 341 throw Exception::Zero(">>> (PDE_adv_diff::Hessian_un_un): Hessian is zero."); 342 } 343 Hessian_un_zf(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess,const ROL::TimeStamp<Real> & ts,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)344 void Hessian_un_zf(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess, 345 const ROL::TimeStamp<Real> & ts, 346 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff, 347 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff, 348 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff, 349 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff = ROL::nullPtr, 350 const ROL::Ptr<const std::vector<Real>> & z_param = ROL::nullPtr) { 351 throw Exception::Zero(">>> (PDE_adv_diff::Hessian_un_zf): Hessian is zero."); 352 } 353 Hessian_un_zp(std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>> & hess,const ROL::TimeStamp<Real> & ts,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)354 void Hessian_un_zp(std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>> & hess, 355 const ROL::TimeStamp<Real> & ts, 356 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff, 357 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff, 358 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff, 359 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff = ROL::nullPtr, 360 const ROL::Ptr<const std::vector<Real>> & z_param = ROL::nullPtr) { 361 throw Exception::Zero(">>> (PDE_adv_diff::Hessian_un_zp): Hessian is zero."); 362 } 363 Hessian_zf_uo(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess,const ROL::TimeStamp<Real> & ts,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)364 void Hessian_zf_uo(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess, 365 const ROL::TimeStamp<Real> & ts, 366 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff, 367 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff, 368 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff, 369 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff = ROL::nullPtr, 370 const ROL::Ptr<const std::vector<Real>> & z_param = ROL::nullPtr) { 371 throw Exception::Zero(">>> (PDE_adv_diff::Hessian_zf_uo): Hessian is zero."); 372 } 373 Hessian_zf_un(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess,const ROL::TimeStamp<Real> & ts,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)374 void Hessian_zf_un(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess, 375 const ROL::TimeStamp<Real> & ts, 376 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff, 377 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff, 378 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff, 379 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff = ROL::nullPtr, 380 const ROL::Ptr<const std::vector<Real>> & z_param = ROL::nullPtr) { 381 throw Exception::Zero(">>> (PDE_adv_diff::Hessian_zf_un): Hessian is zero."); 382 } 383 Hessian_zf_zf(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess,const ROL::TimeStamp<Real> & ts,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)384 void Hessian_zf_zf(ROL::Ptr<Intrepid::FieldContainer<Real>> & hess, 385 const ROL::TimeStamp<Real> & ts, 386 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff, 387 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff, 388 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff, 389 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff = ROL::nullPtr, 390 const ROL::Ptr<const std::vector<Real>> & z_param = ROL::nullPtr) { 391 throw Exception::Zero(">>> (PDE_adv_diff::Hessian_zf_zf): Hessian is zero."); 392 } 393 Hessian_zf_zp(std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>> & hess,const ROL::TimeStamp<Real> & ts,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)394 void Hessian_zf_zp(std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>> & hess, 395 const ROL::TimeStamp<Real> & ts, 396 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff, 397 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff, 398 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff, 399 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff = ROL::nullPtr, 400 const ROL::Ptr<const std::vector<Real>> & z_param = ROL::nullPtr) { 401 throw Exception::Zero(">>> (PDE_adv_diff::Hessian_zf_zp): Hessian is zero."); 402 } 403 Hessian_zp_uo(std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>> & hess,const ROL::TimeStamp<Real> & ts,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)404 void Hessian_zp_uo(std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>> & hess, 405 const ROL::TimeStamp<Real> & ts, 406 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff, 407 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff, 408 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff, 409 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff = ROL::nullPtr, 410 const ROL::Ptr<const std::vector<Real>> & z_param = ROL::nullPtr) { 411 throw Exception::Zero(">>> (PDE_adv_diff::Hessian_zp_uo): Hessian is zero."); 412 } 413 Hessian_zp_un(std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>> & hess,const ROL::TimeStamp<Real> & ts,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)414 void Hessian_zp_un(std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>> & hess, 415 const ROL::TimeStamp<Real> & ts, 416 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff, 417 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff, 418 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff, 419 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff = ROL::nullPtr, 420 const ROL::Ptr<const std::vector<Real>> & z_param = ROL::nullPtr) { 421 throw Exception::Zero(">>> (PDE_adv_diff::Hessian_zp_un): Hessian is zero."); 422 } 423 Hessian_zp_zf(std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>> & hess,const ROL::TimeStamp<Real> & ts,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)424 void Hessian_zp_zf(std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>> & hess, 425 const ROL::TimeStamp<Real> & ts, 426 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff, 427 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff, 428 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff, 429 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff = ROL::nullPtr, 430 const ROL::Ptr<const std::vector<Real>> & z_param = ROL::nullPtr) { 431 throw Exception::Zero(">>> (PDE_adv_diff::Hessian_zp_zf): Hessian is zero."); 432 } 433 Hessian_zp_zp(std::vector<std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>>> & hess,const ROL::TimeStamp<Real> & ts,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff,const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff=ROL::nullPtr,const ROL::Ptr<const std::vector<Real>> & z_param=ROL::nullPtr)434 void Hessian_zp_zp(std::vector<std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>>> & hess, 435 const ROL::TimeStamp<Real> & ts, 436 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & l_coeff, 437 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & uo_coeff, 438 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & un_coeff, 439 const ROL::Ptr<const Intrepid::FieldContainer<Real>> & z_coeff = ROL::nullPtr, 440 const ROL::Ptr<const std::vector<Real>> & z_param = ROL::nullPtr) { 441 throw Exception::Zero(">>> (PDE_adv_diff::Hessian_zp_zp): Hessian is zero."); 442 } 443 RieszMap_1(ROL::Ptr<Intrepid::FieldContainer<Real>> & riesz)444 void RieszMap_1(ROL::Ptr<Intrepid::FieldContainer<Real>> & riesz) { 445 pde_->RieszMap_1(riesz); 446 } 447 RieszMap_2(ROL::Ptr<Intrepid::FieldContainer<Real>> & riesz)448 void RieszMap_2(ROL::Ptr<Intrepid::FieldContainer<Real>> & riesz) { 449 pde_->RieszMap_2(riesz); 450 } 451 getFields()452 std::vector<ROL::Ptr<Intrepid::Basis<Real, Intrepid::FieldContainer<Real>>>> getFields() { 453 return pde_->getFields(); 454 } 455 setCellNodes(const ROL::Ptr<Intrepid::FieldContainer<Real>> & volCellNodes,const std::vector<std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>>> & bdryCellNodes,const std::vector<std::vector<std::vector<int>>> & bdryCellLocIds)456 void setCellNodes(const ROL::Ptr<Intrepid::FieldContainer<Real>> &volCellNodes, 457 const std::vector<std::vector<ROL::Ptr<Intrepid::FieldContainer<Real>>>> &bdryCellNodes, 458 const std::vector<std::vector<std::vector<int>>> &bdryCellLocIds) { 459 pde_->setCellNodes(volCellNodes,bdryCellNodes,bdryCellLocIds); 460 bdryCellLocIds_ = bdryCellLocIds; 461 // Finite element definition. 462 fe_vol_ = pde_->getFE(); 463 if (useDBC_) { 464 // Set local boundary DOFs. 465 fidx_ = fe_vol_->getBoundaryDofs(); 466 // Compute Dirichlet values at DOFs. 467 int d = pde_->getFields()[0]->getBaseCellTopology().getDimension(); 468 int numSidesets = bdryCellLocIds_.size(); 469 bdryCellDofValues_.resize(numSidesets); 470 for (int i=0; i<numSidesets; ++i) { 471 int numLocSides = bdryCellLocIds_[i].size(); 472 bdryCellDofValues_[i].resize(numLocSides); 473 for (int j=0; j<numLocSides; ++j) { 474 int c = bdryCellLocIds_[i][j].size(); 475 int f = pde_->getFields()[0]->getCardinality(); 476 bdryCellDofValues_[i][j] = ROL::makePtr<Intrepid::FieldContainer<Real>>(c, f); 477 ROL::Ptr<Intrepid::FieldContainer<Real>> coords = 478 ROL::makePtr<Intrepid::FieldContainer<Real>>(c, f, d); 479 if (c > 0) { 480 fe_vol_->computeDofCoords(coords, bdryCellNodes[i][j]); 481 } 482 for (int k=0; k<c; ++k) { 483 for (int l=0; l<f; ++l) { 484 std::vector<Real> dofpoint(d); 485 for (int m=0; m<d; ++m) { 486 dofpoint[m] = (*coords)(k, l, m); 487 } 488 (*bdryCellDofValues_[i][j])(k, l) = evaluateDirichlet(dofpoint, i, j); 489 } 490 } 491 } 492 } 493 } 494 } 495 getFE(void) const496 const ROL::Ptr<FE<Real>> getFE(void) const { 497 return fe_vol_; 498 } 499 500 private: evaluateDirichlet(const std::vector<Real> & coords,int sideset,int locSideId) const501 Real evaluateDirichlet(const std::vector<Real> & coords, int sideset, int locSideId) const { 502 return static_cast<Real>(0); 503 } 504 }; // DynamicPDE_adv_diff 505 506 #endif 507