Albany: a Trilinos-based PDE code: TimeDepBC

Go to the documentation of this file.
00001 //*****************************************************************//
00002 //    Albany 2.0:  Copyright 2012 Sandia Corporation               //
00003 //    This Software is released under the BSD license detailed     //
00004 //    in the file "license.txt" in the top-level Albany directory  //
00005 //*****************************************************************//
00006 
00007 #include "Teuchos_TestForException.hpp"
00008 #include "Phalanx_DataLayout.hpp"
00009 #include "Sacado_ParameterRegistration.hpp"
00010 
00011 // **********************************************************************
00012 // Genereric Template Code for Constructor and PostRegistrationSetup
00013 // **********************************************************************
00014 
00015 namespace LCM {
00016 
00017 template <typename EvalT, typename Traits> 
00018 TimeDepBC_Base<EvalT, Traits>::
00019 TimeDepBC_Base(Teuchos::ParameterList& p) :
00020   offset(p.get<int>("Equation Offset")),
00021   
00022   PHAL::DirichletBase<EvalT, Traits>(p)
00023 {
00024   timeValues = p.get<Teuchos::Array<RealType> >("Time Values").toVector();
00025   BCValues = p.get<Teuchos::Array<RealType> >("BC Values").toVector();
00026 
00027   TEUCHOS_TEST_FOR_EXCEPTION( !(timeValues.size() == BCValues.size()),
00028             Teuchos::Exceptions::InvalidParameter,
00029             "Dimension of \"Time Values\" and \"BC Values\" do not match" );
00030 }
00031 
00032 // **********************************************************************
00033 template<typename EvalT, typename Traits>
00034 typename TimeDepBC_Base<EvalT, Traits>::ScalarT
00035 TimeDepBC_Base<EvalT, Traits>::
00036 computeVal(RealType time)
00037 {
00038   TEUCHOS_TEST_FOR_EXCEPTION( time > timeValues.back(),
00039             Teuchos::Exceptions::InvalidParameter,
00040             "Time is growing unbounded!" );
00041   ScalarT Val;
00042   RealType slope;
00043   unsigned int Index(0);
00044 
00045   while( timeValues[Index] < time )
00046     Index++;
00047 
00048   if (Index == 0)
00049     Val = BCValues[Index];
00050   else
00051   {
00052     slope = ( BCValues[Index] - BCValues[Index - 1] ) / ( timeValues[Index] - timeValues[Index - 1] );
00053     Val = BCValues[Index-1] + slope * ( time - timeValues[Index - 1] );
00054   }
00055 
00056   return Val;
00057 }
00058 // **********************************************************************
00059 // Specialization: Residual
00060 // **********************************************************************
00061 template<typename Traits>
00062 TimeDepBC<PHAL::AlbanyTraits::Residual, Traits>::
00063 TimeDepBC(Teuchos::ParameterList& p) :
00064   TimeDepBC_Base<PHAL::AlbanyTraits::Residual, Traits>(p)
00065 {
00066 }
00067 
00068 // **********************************************************************
00069 template<typename Traits>
00070 void 
00071 TimeDepBC<PHAL::AlbanyTraits::Residual, Traits>::
00072 evaluateFields(typename Traits::EvalData dirichletWorkset)
00073 {
00074   Teuchos::RCP<Epetra_Vector> f = dirichletWorkset.f;
00075   Teuchos::RCP<const Epetra_Vector> x = dirichletWorkset.x;
00076   RealType time = dirichletWorkset.current_time;
00077   // Grab the vector off node GIDs for this Node Set ID from the std::map
00078   const std::vector<std::vector<int> >& nsNodes = 
00079     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00080 
00081   int lunk; // global and local indicies into unknown vector
00082   ScalarT Val;
00083   
00084   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00085     lunk = nsNodes[inode][this->offset];
00086 
00087     Val = this->computeVal(time);
00088 
00089     (*f)[lunk] = ((*x)[lunk] - Val);
00090   }
00091 }
00092 
00093 // **********************************************************************
00094 // Specialization: Jacobian
00095 // **********************************************************************
00096 template<typename Traits>
00097 TimeDepBC<PHAL::AlbanyTraits::Jacobian, Traits>::
00098 TimeDepBC(Teuchos::ParameterList& p) :
00099   TimeDepBC_Base<PHAL::AlbanyTraits::Jacobian, Traits>(p)
00100 {
00101 }
00102 // **********************************************************************
00103 template<typename Traits>
00104 void TimeDepBC<PHAL::AlbanyTraits::Jacobian, Traits>::
00105 evaluateFields(typename Traits::EvalData dirichletWorkset)
00106 {
00107 
00108   Teuchos::RCP<Epetra_Vector> f = dirichletWorkset.f;
00109   Teuchos::RCP<Epetra_CrsMatrix> jac = dirichletWorkset.Jac;
00110   Teuchos::RCP<const Epetra_Vector> x = dirichletWorkset.x;
00111   RealType time = dirichletWorkset.current_time;
00112   const RealType j_coeff = dirichletWorkset.j_coeff;
00113   const std::vector<std::vector<int> >& nsNodes = 
00114     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00115   const std::vector<double*>& nsNodeCoords = 
00116     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00117 
00118   RealType* matrixEntries;
00119   int*    matrixIndices;
00120   int     numEntries;
00121   RealType diag=j_coeff;
00122   bool fillResid = (f != Teuchos::null);
00123 
00124   int lunk; // local indicies into unknown vector
00125   ScalarT Val; 
00126   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) 
00127   {
00128     lunk = nsNodes[inode][this->offset];
00129     
00130     Val = this->computeVal(time);
00131     
00132     // replace jac values for the X dof 
00133     jac->ExtractMyRowView(lunk, numEntries, matrixEntries, matrixIndices);
00134     for (int i=0; i<numEntries; i++) matrixEntries[i]=0;
00135     jac->ReplaceMyValues(lunk, 1, &diag, &lunk);
00136     
00137     if (fillResid)  (*f)[lunk] = ((*x)[lunk] - Val.val());
00138   }
00139 }
00140 
00141 // **********************************************************************
00142 // Specialization: Tangent
00143 // **********************************************************************
00144 template<typename Traits>
00145 TimeDepBC<PHAL::AlbanyTraits::Tangent, Traits>::
00146 TimeDepBC(Teuchos::ParameterList& p) :
00147   TimeDepBC_Base<PHAL::AlbanyTraits::Tangent, Traits>(p)
00148 {
00149 }
00150 // **********************************************************************
00151 template<typename Traits>
00152 void TimeDepBC<PHAL::AlbanyTraits::Tangent, Traits>::
00153 evaluateFields(typename Traits::EvalData dirichletWorkset)
00154 {
00155   Teuchos::RCP<Epetra_Vector> f = dirichletWorkset.f;
00156   Teuchos::RCP<Epetra_MultiVector> fp = dirichletWorkset.fp;
00157   Teuchos::RCP<Epetra_MultiVector> JV = dirichletWorkset.JV;
00158   Teuchos::RCP<const Epetra_Vector> x = dirichletWorkset.x;
00159   Teuchos::RCP<const Epetra_MultiVector> Vx = dirichletWorkset.Vx;
00160   RealType time = dirichletWorkset.current_time;
00161   const RealType j_coeff = dirichletWorkset.j_coeff;
00162   const std::vector<std::vector<int> >& nsNodes = 
00163     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00164   const std::vector<double*>& nsNodeCoords = 
00165     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00166 
00167   int lunk; // global and local indicies into unknown vector
00168   ScalarT Val;
00169 
00170   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) 
00171   {
00172     lunk = nsNodes[inode][this->offset];
00173     
00174     Val = this->computeVal(time);
00175 
00176     if (f != Teuchos::null)
00177       (*f)[lunk] = ((*x)[lunk] - Val.val());
00178 
00179     if (JV != Teuchos::null) 
00180       for (int i=0; i<dirichletWorkset.num_cols_x; i++)
00181   (*JV)[i][lunk] = j_coeff*(*Vx)[i][lunk];
00182     
00183     if (fp != Teuchos::null) 
00184       for (int i=0; i<dirichletWorkset.num_cols_p; i++)
00185       (*fp)[i][lunk] = -Val.dx(dirichletWorkset.param_offset+i);
00186   }
00187 }
00188 
00189 // **********************************************************************
00190 // Specialization: Stochastic Galerkin Residual
00191 // **********************************************************************
00192 #ifdef ALBANY_SG_MP
00193 template<typename Traits>
00194 TimeDepBC<PHAL::AlbanyTraits::SGResidual, Traits>::
00195 TimeDepBC(Teuchos::ParameterList& p) :
00196   TimeDepBC_Base<PHAL::AlbanyTraits::SGResidual, Traits>(p)
00197 {
00198 }
00199 // **********************************************************************
00200 template<typename Traits>
00201 void TimeDepBC<PHAL::AlbanyTraits::SGResidual, Traits>::
00202 evaluateFields(typename Traits::EvalData dirichletWorkset)
00203 {
00204   Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> f = 
00205     dirichletWorkset.sg_f;
00206   Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly> x = 
00207     dirichletWorkset.sg_x;
00208   const std::vector<std::vector<int> >& nsNodes = 
00209     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00210   const std::vector<double*>& nsNodeCoords = 
00211     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00212 
00213   RealType time = dirichletWorkset.current_time;
00214   int lunk; // global and local indicies into unknown vector
00215   ScalarT Val;
00216   int nblock = x->size();
00217   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00218     lunk = nsNodes[inode][this->offset];
00219 
00220     Val = this->computeVal(time);
00221 
00222     for (int block=0; block<nblock; block++)
00223       (*f)[block][lunk] = ((*x)[block][lunk] - Val.coeff(block));
00224   }
00225 }
00226 
00227 // **********************************************************************
00228 // Specialization: Stochastic Galerkin Jacobian
00229 // **********************************************************************
00230 template<typename Traits>
00231 TimeDepBC<PHAL::AlbanyTraits::SGJacobian, Traits>::
00232 TimeDepBC(Teuchos::ParameterList& p) :
00233   TimeDepBC_Base<PHAL::AlbanyTraits::SGJacobian, Traits>(p)
00234 {
00235 }
00236 // **********************************************************************
00237 template<typename Traits>
00238 void TimeDepBC<PHAL::AlbanyTraits::SGJacobian, Traits>::
00239 evaluateFields(typename Traits::EvalData dirichletWorkset)
00240 {
00241   Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly> f = 
00242     dirichletWorkset.sg_f;
00243   Teuchos::RCP< Stokhos::VectorOrthogPoly<Epetra_CrsMatrix> > jac = 
00244     dirichletWorkset.sg_Jac;
00245   Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly> x = 
00246     dirichletWorkset.sg_x;
00247   const RealType j_coeff = dirichletWorkset.j_coeff;
00248   const std::vector<std::vector<int> >& nsNodes = 
00249     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00250   const std::vector<double*>& nsNodeCoords = 
00251     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00252 
00253   RealType* matrixEntries;
00254   int*    matrixIndices;
00255   int     numEntries;
00256   RealType diag=j_coeff;
00257   bool fillResid = (f != Teuchos::null);
00258 
00259   int nblock = 0;
00260   if (f != Teuchos::null)
00261     nblock = f->size();
00262   int nblock_jac = jac->size();
00263 
00264   RealType time = dirichletWorkset.current_time;
00265   int lunk; // local indicies into unknown vector
00266   ScalarT Val;
00267 
00268   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) 
00269   {
00270     lunk = nsNodes[inode][this->offset];
00271     
00272     Val = this->computeVal(time);
00273     
00274     // replace jac values for the X dof 
00275     for (int block=0; block<nblock_jac; block++) 
00276     {
00277       (*jac)[block].ExtractMyRowView(lunk, numEntries, matrixEntries, 
00278              matrixIndices);
00279       for (int i=0; i<numEntries; i++) matrixEntries[i]=0;
00280     }
00281     (*jac)[0].ReplaceMyValues(lunk, 1, &diag, &lunk);
00282     
00283     if (fillResid)
00284       for (int block=0; block<nblock; block++) 
00285   (*f)[block][lunk] = ((*x)[block][lunk] - Val.val().coeff(block));
00286   }
00287 }
00288 
00289 // **********************************************************************
00290 // Specialization: Stochastic Galerkin Tangent
00291 // **********************************************************************
00292 template<typename Traits>
00293 TimeDepBC<PHAL::AlbanyTraits::SGTangent, Traits>::
00294 TimeDepBC(Teuchos::ParameterList& p) :
00295   TimeDepBC_Base<PHAL::AlbanyTraits::SGTangent, Traits>(p)
00296 {
00297 }
00298 // **********************************************************************
00299 template<typename Traits>
00300 void TimeDepBC<PHAL::AlbanyTraits::SGTangent, Traits>::
00301 evaluateFields(typename Traits::EvalData dirichletWorkset)
00302 {
00303   Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> f = 
00304     dirichletWorkset.sg_f;
00305   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> fp = 
00306     dirichletWorkset.sg_fp;
00307   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> JV = 
00308     dirichletWorkset.sg_JV;
00309   Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly> x = 
00310     dirichletWorkset.sg_x;
00311   Teuchos::RCP<const Epetra_MultiVector> Vx = dirichletWorkset.Vx;
00312   const RealType j_coeff = dirichletWorkset.j_coeff;
00313   const std::vector<std::vector<int> >& nsNodes = 
00314     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00315   const std::vector<double*>& nsNodeCoords = 
00316     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00317 
00318   int nblock = x->size();
00319 
00320   if (JV != Teuchos::null)
00321     JV->init(0.0);
00322   if (fp != Teuchos::null)
00323     fp->init(0.0);
00324 
00325   RealType time = dirichletWorkset.current_time;
00326   int lunk; // global and local indicies into unknown vector
00327   ScalarT Val;
00328   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) 
00329   {
00330     lunk = nsNodes[inode][this->offset];
00331     
00332     Val = this->computeVal(time);
00333 
00334     if (f != Teuchos::null)
00335       for (int block=0; block<nblock; block++) 
00336   (*f)[block][lunk] = (*x)[block][lunk] - Val.val().coeff(block);
00337 
00338     if (JV != Teuchos::null) 
00339       for (int i=0; i<dirichletWorkset.num_cols_x; i++)
00340   (*JV)[0][i][lunk] = j_coeff*(*Vx)[i][lunk];
00341     
00342     if (fp != Teuchos::null) 
00343       for (int i=0; i<dirichletWorkset.num_cols_p; i++)
00344   for (int block=0; block<nblock; block++) 
00345     (*fp)[block][i][lunk] = 
00346       -Val.dx(dirichletWorkset.param_offset+i).coeff(block);
00347   }
00348 }
00349 
00350 
00351 // **********************************************************************
00352 // Specialization: Multi-point Residual
00353 // **********************************************************************
00354 template<typename Traits>
00355 TimeDepBC<PHAL::AlbanyTraits::MPResidual, Traits>::
00356 TimeDepBC(Teuchos::ParameterList& p) :
00357   TimeDepBC_Base<PHAL::AlbanyTraits::MPResidual, Traits>(p)
00358 {
00359 }
00360 // **********************************************************************
00361 template<typename Traits>
00362 void TimeDepBC<PHAL::AlbanyTraits::MPResidual, Traits>::
00363 evaluateFields(typename Traits::EvalData dirichletWorkset)
00364 {
00365   Teuchos::RCP<Stokhos::ProductEpetraVector> f = 
00366     dirichletWorkset.mp_f;
00367   Teuchos::RCP<const Stokhos::ProductEpetraVector> x = 
00368     dirichletWorkset.mp_x;
00369   const std::vector<std::vector<int> >& nsNodes = 
00370     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00371   const std::vector<double*>& nsNodeCoords = 
00372     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00373 
00374   RealType time = dirichletWorkset.current_time;
00375   int lunk; // global and local indicies into unknown vector
00376   ScalarT Val;
00377   
00378   int nblock = x->size();
00379   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) 
00380   {
00381     lunk = nsNodes[inode][0];
00382 
00383     Val = this->computeVal(time);
00384 
00385     for (int block=0; block<nblock; block++) 
00386       (*f)[block][lunk] = ((*x)[block][lunk] - Val.coeff(block));
00387   }
00388 }
00389 
00390 // **********************************************************************
00391 // Specialization: Multi-point Jacobian
00392 // **********************************************************************
00393 template<typename Traits>
00394 TimeDepBC<PHAL::AlbanyTraits::MPJacobian, Traits>::
00395 TimeDepBC(Teuchos::ParameterList& p) :
00396   TimeDepBC_Base<PHAL::AlbanyTraits::MPJacobian, Traits>(p)
00397 {
00398 }
00399 // **********************************************************************
00400 template<typename Traits>
00401 void TimeDepBC<PHAL::AlbanyTraits::MPJacobian, Traits>::
00402 evaluateFields(typename Traits::EvalData dirichletWorkset)
00403 {
00404   Teuchos::RCP<Stokhos::ProductEpetraVector> f = 
00405     dirichletWorkset.mp_f;
00406   Teuchos::RCP< Stokhos::ProductContainer<Epetra_CrsMatrix> > jac = 
00407     dirichletWorkset.mp_Jac;
00408   Teuchos::RCP<const Stokhos::ProductEpetraVector> x = 
00409     dirichletWorkset.mp_x;
00410   const RealType j_coeff = dirichletWorkset.j_coeff;
00411   const std::vector<std::vector<int> >& nsNodes = 
00412     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00413   const std::vector<double*>& nsNodeCoords = 
00414     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00415 
00416   RealType* matrixEntries;
00417   int*    matrixIndices;
00418   int     numEntries;
00419   RealType diag=j_coeff;
00420   bool fillResid = (f != Teuchos::null);
00421 
00422   int nblock = 0;
00423   if (f != Teuchos::null)
00424     nblock = f->size();
00425   int nblock_jac = jac->size();
00426 
00427   RealType time = dirichletWorkset.current_time;
00428   int lunk; // local indicies into unknown vector
00429   ScalarT Val;
00430   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) 
00431   {
00432     lunk = nsNodes[inode][this->offset];
00433     
00434     Val = this->computeVal(time);
00435     
00436     // replace jac values
00437     for (int block=0; block<nblock_jac; block++) 
00438     {
00439       (*jac)[block].ExtractMyRowView(lunk, numEntries, matrixEntries, 
00440              matrixIndices);
00441       for (int i=0; i<numEntries; i++) 
00442   matrixEntries[i]=0;
00443       (*jac)[block].ReplaceMyValues(lunk, 1, &diag, &lunk);
00444     }
00445     
00446     if (fillResid)
00447       for (int block=0; block<nblock; block++) 
00448   (*f)[block][lunk] = ((*x)[block][lunk] - Val.val().coeff(block));
00449   }
00450 }
00451 
00452 // **********************************************************************
00453 // Specialization: Multi-point Tangent
00454 // **********************************************************************
00455 template<typename Traits>
00456 TimeDepBC<PHAL::AlbanyTraits::MPTangent, Traits>::
00457 TimeDepBC(Teuchos::ParameterList& p) :
00458   TimeDepBC_Base<PHAL::AlbanyTraits::MPTangent, Traits>(p)
00459 {
00460 }
00461 // **********************************************************************
00462 template<typename Traits>
00463 void TimeDepBC<PHAL::AlbanyTraits::MPTangent, Traits>::
00464 evaluateFields(typename Traits::EvalData dirichletWorkset)
00465 {
00466   Teuchos::RCP<Stokhos::ProductEpetraVector> f = 
00467     dirichletWorkset.mp_f;
00468   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> fp = 
00469     dirichletWorkset.mp_fp;
00470   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> JV = 
00471     dirichletWorkset.mp_JV;
00472   Teuchos::RCP<const Stokhos::ProductEpetraVector> x = 
00473     dirichletWorkset.mp_x;
00474   Teuchos::RCP<const Epetra_MultiVector> Vx = dirichletWorkset.Vx;
00475   const RealType j_coeff = dirichletWorkset.j_coeff;
00476   const std::vector<std::vector<int> >& nsNodes = 
00477     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00478   const std::vector<double*>& nsNodeCoords = 
00479     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00480 
00481   int nblock = x->size();
00482 
00483   RealType time = dirichletWorkset.current_time;
00484   int lunk; // global and local indicies into unknown vector
00485   ScalarT Val;
00486   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) 
00487   {
00488     lunk = nsNodes[inode][this->offset];
00489     
00490     Val = this->computeVal(time);
00491 
00492     if (f != Teuchos::null)
00493       for (int block=0; block<nblock; block++) 
00494   (*f)[block][lunk] = (*x)[block][lunk] - Val.val().coeff(block);
00495 
00496 
00497     if (JV != Teuchos::null) 
00498       for (int i=0; i<dirichletWorkset.num_cols_x; i++)
00499   for (int block=0; block<nblock; block++) 
00500     (*JV)[block][i][lunk] = j_coeff*(*Vx)[i][lunk];
00501 
00502     
00503     if (fp != Teuchos::null) 
00504       for (int i=0; i<dirichletWorkset.num_cols_p; i++)
00505   for (int block=0; block<nblock; block++)
00506     (*fp)[block][i][lunk] = 
00507       -Val.dx(dirichletWorkset.param_offset+i).coeff(block);
00508   }
00509 }
00510 #endif //ALBANY_SG_MP
00511 
00512 } // namespace LCM
00513
TimeDepBC_Def.hpp
