Albany: a Trilinos-based PDE code: PHAL_SeparableScatterScalarResponse

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 
00010 // **********************************************************************
00011 // Base Class Generic Implemtation
00012 // **********************************************************************
00013 namespace PHAL {
00014 
00015 template<typename EvalT, typename Traits>
00016 SeparableScatterScalarResponseBase<EvalT, Traits>::
00017 SeparableScatterScalarResponseBase(const Teuchos::ParameterList& p,
00018            const Teuchos::RCP<Albany::Layouts>& dl)
00019 {
00020   setup(p, dl);
00021 }
00022 
00023 template<typename EvalT, typename Traits>
00024 void SeparableScatterScalarResponseBase<EvalT, Traits>::
00025 postRegistrationSetup(typename Traits::SetupData d,
00026                       PHX::FieldManager<Traits>& fm)
00027 {
00028   this->utils.setFieldData(local_response,fm);
00029 }
00030 
00031 template<typename EvalT, typename Traits>
00032 void
00033 SeparableScatterScalarResponseBase<EvalT, Traits>::
00034 setup(const Teuchos::ParameterList& p, const Teuchos::RCP<Albany::Layouts>& dl)
00035 {
00036   bool stand_alone = p.get<bool>("Stand-alone Evaluator");
00037 
00038   // Setup fields we require
00039   PHX::Tag<ScalarT> local_response_tag =
00040     p.get<PHX::Tag<ScalarT> >("Local Response Field Tag");
00041   local_response = PHX::MDField<ScalarT>(local_response_tag);
00042   if (stand_alone)
00043     this->addDependentField(local_response);
00044   else
00045     this->addEvaluatedField(local_response);
00046 }
00047 
00048 // **********************************************************************
00049 // Specialization: Jacobian
00050 // **********************************************************************
00051 
00052 template<typename Traits>
00053 SeparableScatterScalarResponse<PHAL::AlbanyTraits::Jacobian, Traits>::
00054 SeparableScatterScalarResponse(const Teuchos::ParameterList& p,
00055     const Teuchos::RCP<Albany::Layouts>& dl) 
00056 {
00057   this->setup(p,dl);
00058 }
00059 
00060 template<typename Traits>
00061 void SeparableScatterScalarResponse<PHAL::AlbanyTraits::Jacobian, Traits>::
00062 preEvaluate(typename Traits::PreEvalData workset)
00063 {
00064   // Initialize derivatives
00065   Teuchos::RCP<Epetra_MultiVector> dgdx = workset.dgdx;
00066   Teuchos::RCP<Epetra_MultiVector> overlapped_dgdx = workset.overlapped_dgdx;
00067   if (dgdx != Teuchos::null) {
00068     dgdx->PutScalar(0.0);
00069     overlapped_dgdx->PutScalar(0.0);
00070   }
00071 
00072   Teuchos::RCP<Epetra_MultiVector> dgdxdot = workset.dgdxdot;
00073   Teuchos::RCP<Epetra_MultiVector> overlapped_dgdxdot = 
00074     workset.overlapped_dgdxdot;
00075   if (dgdxdot != Teuchos::null) {
00076     dgdxdot->PutScalar(0.0);
00077     overlapped_dgdxdot->PutScalar(0.0);
00078   }
00079 }
00080 
00081 template<typename Traits>
00082 void SeparableScatterScalarResponse<PHAL::AlbanyTraits::Jacobian, Traits>::
00083 evaluateFields(typename Traits::EvalData workset)
00084 {
00085   // Here we scatter the *local* response derivative
00086   Teuchos::RCP<Epetra_MultiVector> dgdx = workset.overlapped_dgdx;
00087   Teuchos::RCP<Epetra_MultiVector> dgdxdot = workset.overlapped_dgdxdot;
00088   Teuchos::RCP<Epetra_MultiVector> dg;
00089   if (dgdx != Teuchos::null)
00090     dg = dgdx;
00091   else
00092     dg = dgdxdot;
00093    
00094   // Loop over cells in workset
00095   for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00096     const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID = 
00097       workset.wsElNodeEqID[cell];
00098 
00099     // Loop over responses
00100     for (std::size_t res = 0; res < this->global_response.size(); res++) {
00101       ScalarT& val = this->local_response(cell, res);
00102 
00103       // Loop over nodes in cell
00104       for (unsigned int node_dof=0; node_dof<numNodes; node_dof++) {
00105   int neq = nodeID[node_dof].size();
00106     
00107   // Loop over equations per node
00108   for (unsigned int eq_dof=0; eq_dof<neq; eq_dof++) {
00109       
00110     // local derivative component
00111     int deriv = neq * node_dof + eq_dof;
00112       
00113     // local DOF
00114     int dof = nodeID[node_dof][eq_dof];
00115             
00116     // Set dg/dx
00117     dg->SumIntoMyValue(dof, res, val.dx(deriv));
00118 
00119   } // column equations
00120       } // column nodes
00121     } // response
00122   } // cell 
00123 }
00124 
00125 template<typename Traits>
00126 void SeparableScatterScalarResponse<PHAL::AlbanyTraits::Jacobian, Traits>::
00127 postEvaluate(typename Traits::PostEvalData workset)
00128 {
00129   // Here we scatter the *global* response
00130   Teuchos::RCP<Epetra_Vector> g = workset.g;
00131   if (g != Teuchos::null)
00132     for (std::size_t res = 0; res < this->global_response.size(); res++) {
00133       (*g)[res] = this->global_response[res].val();
00134   }
00135   
00136   // Here we scatter the *global* response derivatives
00137   Teuchos::RCP<Epetra_MultiVector> dgdx = workset.dgdx;
00138   Teuchos::RCP<Epetra_MultiVector> overlapped_dgdx = workset.overlapped_dgdx;
00139   if (dgdx != Teuchos::null)
00140     dgdx->Export(*overlapped_dgdx, *workset.x_importer, Add);
00141 
00142   Teuchos::RCP<Epetra_MultiVector> dgdxdot = workset.dgdxdot;
00143   Teuchos::RCP<Epetra_MultiVector> overlapped_dgdxdot = 
00144     workset.overlapped_dgdxdot;
00145   if (dgdxdot != Teuchos::null) 
00146     dgdxdot->Export(*overlapped_dgdxdot, *workset.x_importer, Add);
00147 }
00148 
00149 // **********************************************************************
00150 // Specialization: Stochastic Galerkin Jacobian
00151 // **********************************************************************
00152 #ifdef ALBANY_SG_MP
00153 template<typename Traits>
00154 SeparableScatterScalarResponse<PHAL::AlbanyTraits::SGJacobian, Traits>::
00155 SeparableScatterScalarResponse(const Teuchos::ParameterList& p,
00156     const Teuchos::RCP<Albany::Layouts>& dl) 
00157 {
00158   this->setup(p,dl);
00159 }
00160 
00161 template<typename Traits>
00162 void SeparableScatterScalarResponse<PHAL::AlbanyTraits::SGJacobian, Traits>::
00163 preEvaluate(typename Traits::PreEvalData workset)
00164 {
00165   // Initialize derivatives
00166   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dgdx_sg = 
00167     workset.sg_dgdx;
00168   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> overlapped_dgdx_sg = 
00169     workset.overlapped_sg_dgdx;
00170   if (dgdx_sg != Teuchos::null) {
00171     dgdx_sg->init(0.0);
00172     overlapped_dgdx_sg->init(0.0);
00173   }
00174 
00175   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dgdxdot_sg = 
00176     workset.sg_dgdxdot;
00177   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> overlapped_dgdxdot_sg = 
00178     workset.overlapped_sg_dgdxdot;
00179   if (dgdxdot_sg != Teuchos::null) {
00180     dgdxdot_sg->init(0.0);
00181     overlapped_dgdxdot_sg->init(0.0);
00182   }
00183 }
00184 
00185 template<typename Traits>
00186 void SeparableScatterScalarResponse<PHAL::AlbanyTraits::SGJacobian, Traits>::
00187 evaluateFields(typename Traits::EvalData workset)
00188 {
00189   // Here we scatter the *local* SG response derivative
00190   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dgdx_sg = 
00191     workset.overlapped_sg_dgdx;
00192   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dgdxdot_sg = 
00193     workset.overlapped_sg_dgdxdot;
00194   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dg_sg;
00195   if (dgdx_sg != Teuchos::null)
00196     dg_sg = dgdx_sg;
00197   else
00198     dg_sg = dgdxdot_sg;
00199 
00200   // Loop over cells in workset
00201   for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00202     const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID = 
00203       workset.wsElNodeEqID[cell];
00204 
00205     // Loop over responses
00206     for (std::size_t res = 0; res < this->global_response.size(); res++) {
00207       ScalarT& val = this->local_response(cell, res);
00208 
00209       // Loop over nodes in cell
00210       for (unsigned int node_dof=0; node_dof<numNodes; node_dof++) {
00211   int neq = nodeID[node_dof].size();
00212     
00213   // Loop over equations per node
00214   for (unsigned int eq_dof=0; eq_dof<neq; eq_dof++) {
00215       
00216     // local derivative component
00217     int deriv = neq * node_dof + eq_dof;
00218       
00219     // local DOF
00220     int dof = nodeID[node_dof][eq_dof];
00221     
00222     // Set dg/dx
00223     for (int block=0; block<dg_sg->size(); block++)
00224       (*dg_sg)[block].SumIntoMyValue(dof, res, 
00225              val.dx(deriv).coeff(block));
00226     
00227   } // column equations
00228       } // response
00229     } // node
00230   } // cell 
00231 }
00232 
00233 template<typename Traits>
00234 void SeparableScatterScalarResponse<PHAL::AlbanyTraits::SGJacobian, Traits>::
00235 postEvaluate(typename Traits::PostEvalData workset)
00236 {
00237   // Here we scatter the *global* SG response
00238   Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > g_sg = workset.sg_g;
00239   if (g_sg != Teuchos::null) {
00240     for (std::size_t res = 0; res < this->global_response.size(); res++) {
00241       ScalarT& val = this->global_response[res];
00242       for (int block=0; block<g_sg->size(); block++)
00243   (*g_sg)[block][res] = val.val().coeff(block);
00244     }
00245   }
00246 
00247   // Here we scatter the *global* SG response derivatives
00248   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dgdx_sg = 
00249     workset.sg_dgdx;
00250   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> overlapped_dgdx_sg = 
00251     workset.overlapped_sg_dgdx;
00252   if (dgdx_sg != Teuchos::null)
00253     for (int block=0; block<dgdx_sg->size(); block++)
00254       (*dgdx_sg)[block].Export((*overlapped_dgdx_sg)[block], 
00255              *workset.x_importer, Add);
00256 
00257   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dgdxdot_sg = 
00258     workset.sg_dgdxdot;
00259   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> overlapped_dgdxdot_sg = 
00260     workset.overlapped_sg_dgdxdot;
00261   if (dgdxdot_sg != Teuchos::null) 
00262     for (int block=0; block<dgdxdot_sg->size(); block++)
00263       (*dgdxdot_sg)[block].Export((*overlapped_dgdxdot_sg)[block], 
00264           *workset.x_importer, Add);
00265 }
00266 
00267 // **********************************************************************
00268 // Specialization: Multi-point Jacobian
00269 // **********************************************************************
00270 
00271 template<typename Traits>
00272 SeparableScatterScalarResponse<PHAL::AlbanyTraits::MPJacobian, Traits>::
00273 SeparableScatterScalarResponse(const Teuchos::ParameterList& p,
00274        const Teuchos::RCP<Albany::Layouts>& dl) 
00275 {
00276   this->setup(p,dl);
00277 }
00278 
00279 template<typename Traits>
00280 void SeparableScatterScalarResponse<PHAL::AlbanyTraits::MPJacobian, Traits>::
00281 preEvaluate(typename Traits::PreEvalData workset)
00282 {
00283   // Initialize derivatives
00284   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dgdx_mp = 
00285     workset.mp_dgdx;
00286   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> overlapped_dgdx_mp = 
00287     workset.overlapped_mp_dgdx;
00288   if (dgdx_mp != Teuchos::null) {
00289     dgdx_mp->init(0.0);
00290     overlapped_dgdx_mp->init(0.0);
00291   }
00292 
00293   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dgdxdot_mp = 
00294     workset.mp_dgdxdot;
00295   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> overlapped_dgdxdot_mp = 
00296     workset.overlapped_mp_dgdxdot;
00297   if (dgdxdot_mp != Teuchos::null) {
00298     dgdxdot_mp->init(0.0);
00299     overlapped_dgdxdot_mp->init(0.0);
00300   }
00301 }
00302 
00303 template<typename Traits>
00304 void SeparableScatterScalarResponse<PHAL::AlbanyTraits::MPJacobian, Traits>::
00305 evaluateFields(typename Traits::EvalData workset)
00306 {
00307   // Here we scatter the *local* MP response derivative
00308   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dgdx_mp = 
00309     workset.overlapped_mp_dgdx;
00310   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dgdxdot_mp = 
00311     workset.overlapped_mp_dgdxdot;
00312   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dg_mp;
00313   if (dgdx_mp != Teuchos::null)
00314     dg_mp = dgdx_mp;
00315   else
00316     dg_mp = dgdxdot_mp;
00317 
00318   // Loop over cells in workset
00319   for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00320     const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID = 
00321       workset.wsElNodeEqID[cell];
00322 
00323     // Loop over responses
00324     for (std::size_t res = 0; res < this->global_response.size(); res++) {
00325       ScalarT& val = this->local_response(cell, res);
00326 
00327       // Loop over nodes in cell
00328       for (unsigned int node_dof=0; node_dof<numNodes; node_dof++) {
00329   int neq = nodeID[node_dof].size();
00330   
00331   // Loop over equations per node
00332   for (unsigned int eq_dof=0; eq_dof<neq; eq_dof++) {
00333     
00334     // local derivative component
00335     int deriv = neq * node_dof + eq_dof;
00336     
00337     // local DOF
00338     int dof = nodeID[node_dof][eq_dof];
00339           
00340     // Set dg/dx
00341     for (int block=0; block<dg_mp->size(); block++)
00342       (*dg_mp)[block].SumIntoMyValue(dof, res, 
00343              val.dx(deriv).coeff(block));
00344     
00345   } // column equations
00346       } // response
00347     } // node
00348   } // cell 
00349 }
00350 
00351 template<typename Traits>
00352 void SeparableScatterScalarResponse<PHAL::AlbanyTraits::MPJacobian, Traits>::
00353 postEvaluate(typename Traits::PostEvalData workset)
00354 {
00355   // Here we scatter the *global* MP response
00356   Teuchos::RCP<Stokhos::ProductEpetraVector> g_mp = workset.mp_g;
00357   if (g_mp != Teuchos::null) {
00358     for (std::size_t res = 0; res < this->global_response.size(); res++) {
00359       ScalarT& val = this->global_response[res];
00360       for (int block=0; block<g_mp->size(); block++)
00361   (*g_mp)[block][res] = val.val().coeff(block);
00362     }
00363   }
00364 
00365   // Here we scatter the *global* MP response derivatives
00366   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dgdx_mp = 
00367     workset.mp_dgdx;
00368   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> overlapped_dgdx_mp = 
00369     workset.overlapped_mp_dgdx;
00370   if (dgdx_mp != Teuchos::null)
00371     for (int block=0; block<dgdx_mp->size(); block++)
00372       (*dgdx_mp)[block].Export((*overlapped_dgdx_mp)[block], 
00373              *workset.x_importer, Add);
00374 
00375   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dgdxdot_mp = 
00376     workset.mp_dgdxdot;
00377   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> overlapped_dgdxdot_mp = 
00378     workset.overlapped_mp_dgdxdot;
00379   if (dgdxdot_mp != Teuchos::null) 
00380     for (int block=0; block<dgdxdot_mp->size(); block++)
00381       (*dgdxdot_mp)[block].Export((*overlapped_dgdxdot_mp)[block], 
00382           *workset.x_importer, Add);
00383 }
00384 #endif //ALBANY_SG_MP
00385 
00386 }
00387
PHAL_SeparableScatterScalarResponse_Def.hpp
