Albany: a Trilinos-based PDE code: QCAD_SchrodingerDirichlet

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 QCAD {
00016 
00017 template<typename EvalT,typename Traits>
00018 SchrodingerDirichletBase<EvalT, Traits>::
00019 SchrodingerDirichletBase(Teuchos::ParameterList& p) :
00020   offset(p.get<int>("Equation Offset")),
00021   nodeSetID(p.get<std::string>("Node Set ID"))
00022 {
00023   value = p.get<RealType>("Dirichlet Value");
00024 
00025   std::string name = p.get< std::string >("Dirichlet Name");
00026   PHX::Tag<ScalarT> fieldTag(name, p.get< Teuchos::RCP<PHX::DataLayout> >("Data Layout"));
00027 
00028   this->addEvaluatedField(fieldTag);
00029 
00030   this->setName(name+PHX::TypeString<EvalT>::value);
00031 
00032   // Set up values as parameters for parameter library
00033   Teuchos::RCP<ParamLib> paramLib = p.get< Teuchos::RCP<ParamLib> >
00034                ("Parameter Library", Teuchos::null);
00035 
00036   new Sacado::ParameterRegistration<EvalT, SPL_Traits> (name, this, paramLib);
00037 }
00038 
00039 template<typename EvalT, typename Traits>
00040 void SchrodingerDirichletBase<EvalT, Traits>::
00041 postRegistrationSetup(typename Traits::SetupData d,
00042                       PHX::FieldManager<Traits>& fm)
00043 {
00044 }
00045 
00046 // **********************************************************************
00047 // Specialization: Residual
00048 // **********************************************************************
00049 template<typename Traits>
00050 SchrodingerDirichlet<PHAL::AlbanyTraits::Residual, Traits>::
00051 SchrodingerDirichlet(Teuchos::ParameterList& p) :
00052   SchrodingerDirichletBase<PHAL::AlbanyTraits::Residual, Traits>(p)
00053 {
00054 }
00055 
00056 // **********************************************************************
00057 template<typename Traits>
00058 void SchrodingerDirichlet<PHAL::AlbanyTraits::Residual, Traits>::
00059 evaluateFields(typename Traits::EvalData dirichletWorkset)
00060 {
00061   Teuchos::RCP<Epetra_Vector> f = dirichletWorkset.f;
00062   Teuchos::RCP<const Epetra_Vector> x = dirichletWorkset.x;
00063   // Grab the vector off node GIDs for this Node Set ID from the std::map
00064   const std::vector<std::vector<int> >& nsNodes = dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00065 
00066   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00067       int lunk = nsNodes[inode][this->offset];
00068       (*f)[lunk] = ((*x)[lunk] - this->value);
00069   }
00070 }
00071 
00072 // **********************************************************************
00073 // Specialization: Jacobian
00074 // **********************************************************************
00075 template<typename Traits>
00076 SchrodingerDirichlet<PHAL::AlbanyTraits::Jacobian, Traits>::
00077 SchrodingerDirichlet(Teuchos::ParameterList& p) :
00078   SchrodingerDirichletBase<PHAL::AlbanyTraits::Jacobian, Traits>(p)
00079 {
00080 }
00081 
00082 // **********************************************************************
00083 template<typename Traits>
00084 void SchrodingerDirichlet<PHAL::AlbanyTraits::Jacobian, Traits>::
00085 evaluateFields(typename Traits::EvalData dirichletWorkset)
00086 {
00087 
00088   Teuchos::RCP<Epetra_Vector> f = dirichletWorkset.f;
00089   Teuchos::RCP<Epetra_CrsMatrix> jac = dirichletWorkset.Jac;
00090   Teuchos::RCP<const Epetra_Vector> x = dirichletWorkset.x;
00091   const RealType j_coeff = dirichletWorkset.j_coeff;
00092   const std::vector<std::vector<int> >& nsNodes = dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00093 
00094   RealType* matrixEntries;
00095   int*    matrixIndices;
00096   int     numEntries;
00097   RealType diag=j_coeff;
00098   bool fillResid = (f != Teuchos::null);
00099 
00100   int nMyRows = jac->NumMyRows();
00101   RealType zero = 0.0;
00102   //std::vector<int> globalRows(nMyRows);
00103   //jac->RowMap().MyGlobalElements(&globalRows[0]);
00104 
00105   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00106       int lunk = nsNodes[inode][this->offset];
00107       jac->ExtractMyRowView(lunk, numEntries, matrixEntries, matrixIndices);
00108       for (int i=0; i<numEntries; i++) matrixEntries[i]=0; // zero out row
00109       for (int i=0; i<nMyRows; i++) jac->ReplaceMyValues(i, 1, &zero, &lunk); //zero out col
00110       //int gunk = globalRows[lunk]; // convert local row index -> global index
00111 
00112       jac->ReplaceMyValues(lunk, 1, &diag, &lunk); //set diagonal element = j_coeff
00113 
00114       if (fillResid) (*f)[lunk] = ((*x)[lunk] - this->value.val());
00115   }
00116 }
00117 
00118 // **********************************************************************
00119 // Specialization: Tangent
00120 // **********************************************************************
00121 template<typename Traits>
00122 SchrodingerDirichlet<PHAL::AlbanyTraits::Tangent, Traits>::
00123 SchrodingerDirichlet(Teuchos::ParameterList& p) :
00124   SchrodingerDirichletBase<PHAL::AlbanyTraits::Tangent, Traits>(p)
00125 {
00126 }
00127 
00128 // **********************************************************************
00129 template<typename Traits>
00130 void SchrodingerDirichlet<PHAL::AlbanyTraits::Tangent, Traits>::
00131 evaluateFields(typename Traits::EvalData dirichletWorkset)
00132 {
00133 
00134   Teuchos::RCP<Epetra_Vector> f = dirichletWorkset.f;
00135   Teuchos::RCP<Epetra_MultiVector> fp = dirichletWorkset.fp;
00136   Teuchos::RCP<Epetra_MultiVector> JV = dirichletWorkset.JV;
00137   Teuchos::RCP<const Epetra_Vector> x = dirichletWorkset.x;
00138   Teuchos::RCP<const Epetra_MultiVector> Vx = dirichletWorkset.Vx;
00139   const RealType j_coeff = dirichletWorkset.j_coeff;
00140   const std::vector<std::vector<int> >& nsNodes = 
00141     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00142 
00143   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00144     int lunk = nsNodes[inode][this->offset];
00145 
00146     if (f != Teuchos::null) 
00147       (*f)[lunk] = ((*x)[lunk] - this->value.val());
00148     
00149     if (JV != Teuchos::null)
00150       for (int i=0; i<dirichletWorkset.num_cols_x; i++)
00151   (*JV)[i][lunk] = j_coeff*(*Vx)[i][lunk];
00152     
00153     if (fp != Teuchos::null)
00154       for (int i=0; i<dirichletWorkset.num_cols_p; i++)
00155   (*fp)[i][lunk] = -this->value.dx(dirichletWorkset.param_offset+i);
00156   }
00157 }
00158 
00159 // **********************************************************************
00160 // Specialization: Stochastic Galerkin Residual
00161 // **********************************************************************
00162 
00163 #ifdef ALBANY_SG_MP
00164 template<typename Traits>
00165 SchrodingerDirichlet<PHAL::AlbanyTraits::SGResidual, Traits>::
00166 SchrodingerDirichlet(Teuchos::ParameterList& p) :
00167   SchrodingerDirichletBase<PHAL::AlbanyTraits::SGResidual, Traits>(p)
00168 {
00169 }
00170 // **********************************************************************
00171 template<typename Traits>
00172 void SchrodingerDirichlet<PHAL::AlbanyTraits::SGResidual, Traits>::
00173 evaluateFields(typename Traits::EvalData dirichletWorkset)
00174 {
00175   Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> f = 
00176     dirichletWorkset.sg_f;
00177   Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly> x = 
00178     dirichletWorkset.sg_x;
00179   const std::vector<std::vector<int> >& nsNodes = 
00180     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00181 
00182   int nblock = x->size();
00183   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00184       int lunk = nsNodes[inode][this->offset];
00185       for (int block=0; block<nblock; block++)
00186   (*f)[block][lunk] = ((*x)[block][lunk] - this->value.coeff(block));
00187   }
00188 }
00189 
00190 // **********************************************************************
00191 // Specialization: Stochastic Galerkin Jacobian
00192 // **********************************************************************
00193 template<typename Traits>
00194 SchrodingerDirichlet<PHAL::AlbanyTraits::SGJacobian, Traits>::
00195 SchrodingerDirichlet(Teuchos::ParameterList& p) :
00196   SchrodingerDirichletBase<PHAL::AlbanyTraits::SGJacobian, Traits>(p)
00197 {
00198 }
00199 
00200 // **********************************************************************
00201 template<typename Traits>
00202 void SchrodingerDirichlet<PHAL::AlbanyTraits::SGJacobian, Traits>::
00203 evaluateFields(typename Traits::EvalData dirichletWorkset)
00204 {
00205   Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly> f = 
00206     dirichletWorkset.sg_f;
00207   Teuchos::RCP< Stokhos::VectorOrthogPoly<Epetra_CrsMatrix> > jac = 
00208     dirichletWorkset.sg_Jac;
00209   Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly> x = 
00210     dirichletWorkset.sg_x;
00211   const RealType j_coeff = dirichletWorkset.j_coeff;
00212   const std::vector<std::vector<int> >& nsNodes = 
00213     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00214 
00215   RealType* matrixEntries;
00216   int*    matrixIndices;
00217   int     numEntries;
00218   int nblock = 0;
00219   if (f != Teuchos::null)
00220     nblock = f->size();
00221   int nblock_jac = jac->size();
00222   RealType diag=j_coeff;
00223   bool fillResid = (f != Teuchos::null);
00224 
00225   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00226       int lunk = nsNodes[inode][this->offset];
00227       for (int block=0; block<nblock_jac; block++) {
00228   (*jac)[block].ExtractMyRowView(lunk, numEntries, matrixEntries, 
00229                matrixIndices);
00230   for (int i=0; i<numEntries; i++) 
00231     matrixEntries[i]=0;
00232       }
00233       (*jac)[0].ReplaceMyValues(lunk, 1, &diag, &lunk);
00234       if (fillResid) {
00235   for (int block=0; block<nblock; block++)
00236     (*f)[block][lunk] = 
00237       (*x)[block][lunk] - this->value.val().coeff(block);
00238       }
00239   }
00240 }
00241 
00242 // **********************************************************************
00243 // Specialization: Stochastic Galerkin Tangent
00244 // **********************************************************************
00245 template<typename Traits>
00246 SchrodingerDirichlet<PHAL::AlbanyTraits::SGTangent, Traits>::
00247 SchrodingerDirichlet(Teuchos::ParameterList& p) :
00248   SchrodingerDirichletBase<PHAL::AlbanyTraits::SGTangent, Traits>(p)
00249 {
00250 }
00251 
00252 // **********************************************************************
00253 template<typename Traits>
00254 void SchrodingerDirichlet<PHAL::AlbanyTraits::SGTangent, Traits>::
00255 evaluateFields(typename Traits::EvalData dirichletWorkset)
00256 {
00257 
00258   Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> f = 
00259     dirichletWorkset.sg_f;
00260   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> fp = 
00261     dirichletWorkset.sg_fp;
00262   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> JV = 
00263     dirichletWorkset.sg_JV;
00264   Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly> x = 
00265     dirichletWorkset.sg_x;
00266   Teuchos::RCP<const Epetra_MultiVector> Vx = dirichletWorkset.Vx;
00267   const RealType j_coeff = dirichletWorkset.j_coeff;
00268   const std::vector<std::vector<int> >& nsNodes = 
00269     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00270 
00271   int nblock = x->size();
00272 
00273   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00274     int lunk = nsNodes[inode][this->offset];
00275 
00276     if (f != Teuchos::null) 
00277       for (int block=0; block<nblock; block++)
00278   (*f)[block][lunk] = 
00279     (*x)[block][lunk] - this->value.val().coeff(block);
00280 
00281     if (JV != Teuchos::null)
00282       for (int i=0; i<dirichletWorkset.num_cols_x; i++)
00283   (*JV)[0][i][lunk] = j_coeff*(*Vx)[i][lunk];
00284 
00285     if (fp != Teuchos::null)
00286       for (int i=0; i<dirichletWorkset.num_cols_p; i++)
00287   for (int block=0; block<nblock; block++)
00288     (*fp)[block][i][lunk] = 
00289       -this->value.dx(dirichletWorkset.param_offset+i).coeff(block);
00290   }
00291 }
00292 
00293 // **********************************************************************
00294 // Specialization: Multi-point Residual
00295 // **********************************************************************
00296 
00297 template<typename Traits>
00298 SchrodingerDirichlet<PHAL::AlbanyTraits::MPResidual, Traits>::
00299 SchrodingerDirichlet(Teuchos::ParameterList& p) :
00300   SchrodingerDirichletBase<PHAL::AlbanyTraits::MPResidual, Traits>(p)
00301 {
00302 }
00303 // **********************************************************************
00304 template<typename Traits>
00305 void SchrodingerDirichlet<PHAL::AlbanyTraits::MPResidual, Traits>::
00306 evaluateFields(typename Traits::EvalData dirichletWorkset)
00307 {
00308   Teuchos::RCP<Stokhos::ProductEpetraVector> f = 
00309     dirichletWorkset.mp_f;
00310   Teuchos::RCP<const Stokhos::ProductEpetraVector> x = 
00311     dirichletWorkset.mp_x;
00312   const std::vector<std::vector<int> >& nsNodes = 
00313     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00314 
00315   int nblock = x->size();
00316   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00317       int lunk = nsNodes[inode][this->offset];
00318       for (int block=0; block<nblock; block++)
00319   (*f)[block][lunk] = ((*x)[block][lunk] - this->value.coeff(block));
00320   }
00321 }
00322 
00323 // **********************************************************************
00324 // Specialization: Multi-point Jacobian
00325 // **********************************************************************
00326 template<typename Traits>
00327 SchrodingerDirichlet<PHAL::AlbanyTraits::MPJacobian, Traits>::
00328 SchrodingerDirichlet(Teuchos::ParameterList& p) :
00329   SchrodingerDirichletBase<PHAL::AlbanyTraits::MPJacobian, Traits>(p)
00330 {
00331 }
00332 
00333 // **********************************************************************
00334 template<typename Traits>
00335 void SchrodingerDirichlet<PHAL::AlbanyTraits::MPJacobian, Traits>::
00336 evaluateFields(typename Traits::EvalData dirichletWorkset)
00337 {
00338 
00339   Teuchos::RCP<Stokhos::ProductEpetraVector> f = 
00340     dirichletWorkset.mp_f;
00341   Teuchos::RCP< Stokhos::ProductContainer<Epetra_CrsMatrix> > jac = 
00342     dirichletWorkset.mp_Jac;
00343   Teuchos::RCP<const Stokhos::ProductEpetraVector> x = 
00344     dirichletWorkset.mp_x;
00345   const RealType j_coeff = dirichletWorkset.j_coeff;
00346   const std::vector<std::vector<int> >& nsNodes = 
00347     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00348 
00349   RealType* matrixEntries;
00350   int*    matrixIndices;
00351   int     numEntries;
00352   int nblock = 0;
00353   if (f != Teuchos::null)
00354     nblock = f->size();
00355   int nblock_jac = jac->size();
00356   RealType diag=j_coeff;
00357   bool fillResid = (f != Teuchos::null);
00358 
00359   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00360       int lunk = nsNodes[inode][this->offset];
00361       for (int block=0; block<nblock_jac; block++) {
00362   (*jac)[block].ExtractMyRowView(lunk, numEntries, matrixEntries, 
00363                matrixIndices);
00364   for (int i=0; i<numEntries; i++) 
00365     matrixEntries[i]=0;
00366   (*jac)[block].ReplaceMyValues(lunk, 1, &diag, &lunk);
00367       }
00368       if (fillResid) {
00369   for (int block=0; block<nblock; block++)
00370     (*f)[block][lunk] = 
00371       (*x)[block][lunk] - this->value.val().coeff(block);
00372       }
00373   }
00374 }
00375 
00376 // **********************************************************************
00377 // Specialization: Multi-point Tangent
00378 // **********************************************************************
00379 template<typename Traits>
00380 SchrodingerDirichlet<PHAL::AlbanyTraits::MPTangent, Traits>::
00381 SchrodingerDirichlet(Teuchos::ParameterList& p) :
00382   SchrodingerDirichletBase<PHAL::AlbanyTraits::MPTangent, Traits>(p)
00383 {
00384 }
00385 
00386 // **********************************************************************
00387 template<typename Traits>
00388 void SchrodingerDirichlet<PHAL::AlbanyTraits::MPTangent, Traits>::
00389 evaluateFields(typename Traits::EvalData dirichletWorkset)
00390 {
00391 
00392   Teuchos::RCP<Stokhos::ProductEpetraVector> f = 
00393     dirichletWorkset.mp_f;
00394   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> fp = 
00395     dirichletWorkset.mp_fp;
00396   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> JV = 
00397     dirichletWorkset.mp_JV;
00398   Teuchos::RCP<const Stokhos::ProductEpetraVector> x = 
00399     dirichletWorkset.mp_x;
00400   Teuchos::RCP<const Epetra_MultiVector> Vx = dirichletWorkset.Vx;
00401   const RealType j_coeff = dirichletWorkset.j_coeff;
00402   const std::vector<std::vector<int> >& nsNodes = 
00403     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00404 
00405   int nblock = x->size();
00406 
00407   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00408     int lunk = nsNodes[inode][this->offset];
00409 
00410     if (f != Teuchos::null) 
00411       for (int block=0; block<nblock; block++)
00412   (*f)[block][lunk] = 
00413     (*x)[block][lunk] - this->value.val().coeff(block);
00414 
00415     if (JV != Teuchos::null)
00416       for (int i=0; i<dirichletWorkset.num_cols_x; i++)
00417   for (int block=0; block<nblock; block++)
00418     (*JV)[block][i][lunk] = j_coeff*(*Vx)[i][lunk];
00419 
00420     if (fp != Teuchos::null)
00421       for (int i=0; i<dirichletWorkset.num_cols_p; i++)
00422   for (int block=0; block<nblock; block++)
00423     (*fp)[block][i][lunk] = 
00424       -this->value.dx(dirichletWorkset.param_offset+i).coeff(block);
00425   }
00426 }
00427 #endif //ALBANY_SG_MP
00428 
00429 // **********************************************************************
00430 // Simple evaluator to aggregate all SchrodingerDirichlet BCs into one "field"
00431 // **********************************************************************
00432 
00433 template<typename EvalT, typename Traits>
00434 SchrodingerDirichletAggregator<EvalT, Traits>::
00435 SchrodingerDirichletAggregator(Teuchos::ParameterList& p) 
00436 {
00437   Teuchos::RCP<PHX::DataLayout> dl =  p.get< Teuchos::RCP<PHX::DataLayout> >("Data Layout");
00438 
00439   std::vector<std::string>& dbcs = *(p.get<std::vector<std::string>* >("DBC Names"));
00440 
00441   for (unsigned int i=0; i<dbcs.size(); i++) {
00442     PHX::Tag<ScalarT> fieldTag(dbcs[i], dl);
00443     this->addDependentField(fieldTag);
00444   }
00445 
00446   PHX::Tag<ScalarT> fieldTag(p.get<std::string>("DBC Aggregator Name"), dl);
00447   this->addEvaluatedField(fieldTag);
00448 
00449   this->setName("SchrodingerDirichlet Aggregator"+PHX::TypeString<EvalT>::value);
00450 }
00451 
00452 // **********************************************************************
00453 }
00454
QCAD_SchrodingerDirichlet_Def.hpp
