Albany: a Trilinos-based PDE code: PHAL_ComprNSResid

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 
00008 #include "Teuchos_TestForException.hpp"
00009 #include "Phalanx_DataLayout.hpp"
00010 
00011 #include "Intrepid_FunctionSpaceTools.hpp"
00012 
00013 namespace PHAL {
00014 
00015 
00016 //**********************************************************************
00017 template<typename EvalT, typename Traits>
00018 ComprNSResid<EvalT, Traits>::
00019 ComprNSResid(const Teuchos::ParameterList& p) :
00020   wBF     (p.get<std::string>                   ("Weighted BF Name"),
00021          p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Scalar Data Layout") ),
00022   wGradBF    (p.get<std::string>                   ("Weighted Gradient BF Name"),
00023          p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Gradient Data Layout") ),
00024   qFluct       (p.get<std::string>                   ("QP Variable Name"),
00025          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
00026   qFluctGrad      (p.get<std::string>                   ("Gradient QP Variable Name"),
00027          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Tensor Data Layout") ),
00028   qFluctDot       (p.get<std::string>                   ("QP Time Derivative Variable Name"),
00029          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
00030   force       (p.get<std::string>              ("Body Force Name"),
00031                p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
00032   mu          (p.get<std::string>                   ("Viscosity Mu QP Variable Name"),
00033          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),  
00034   lambda          (p.get<std::string>                   ("Viscosity Lambda QP Variable Name"),
00035          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),  
00036   kappa          (p.get<std::string>                   ("Viscosity Kappa QP Variable Name"),
00037          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),  
00038   tau11          (p.get<std::string>                   ("Stress Tau11 QP Variable Name"),
00039          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),  
00040   tau12       (p.get<std::string>                   ("Stress Tau12 QP Variable Name"),
00041          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
00042   tau13       (p.get<std::string>                   ("Stress Tau13 QP Variable Name"),
00043          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
00044   tau22       (p.get<std::string>                   ("Stress Tau22 QP Variable Name"),
00045          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
00046   tau23       (p.get<std::string>                   ("Stress Tau23 QP Variable Name"),
00047          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
00048   tau33       (p.get<std::string>                   ("Stress Tau33 QP Variable Name"),
00049          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
00050   Residual   (p.get<std::string>                   ("Residual Name"),
00051               p.get<Teuchos::RCP<PHX::DataLayout> >("Node Vector Data Layout") ) 
00052 {
00053 
00054   //TO DOs: 
00055   //3D 
00056 
00057   this->addDependentField(qFluct);
00058   this->addDependentField(qFluctGrad);
00059   this->addDependentField(qFluctDot);
00060   this->addDependentField(force);
00061   this->addDependentField(mu);
00062   this->addDependentField(kappa);
00063   this->addDependentField(lambda);
00064   this->addDependentField(wBF);
00065   this->addDependentField(wGradBF);
00066   this->addDependentField(tau11);
00067   this->addDependentField(tau12);
00068   this->addDependentField(tau13);
00069   this->addDependentField(tau22);
00070   this->addDependentField(tau23);
00071   this->addDependentField(tau33);
00072 
00073   this->addEvaluatedField(Residual);
00074 
00075 
00076   this->setName("ComprNSResid"+PHX::TypeString<EvalT>::value);
00077 
00078   std::vector<PHX::DataLayout::size_type> dims;
00079   wGradBF.fieldTag().dataLayout().dimensions(dims);
00080   numNodes = dims[1];
00081   numQPs   = dims[2];
00082   numDims  = dims[3];
00083 
00084 
00085   Teuchos::ParameterList* bf_list =
00086   p.get<Teuchos::ParameterList*>("Parameter List");
00087 
00088   qFluct.fieldTag().dataLayout().dimensions(dims);
00089   vecDim  = dims[2];
00090 
00091   gamma_gas = bf_list->get("Gamma", 1.4); 
00092   Rgas = bf_list->get("Gas constant R", 0.714285733);
00093   Pr = bf_list->get("Prandtl number Pr", 0.72); 
00094   Re = bf_list->get("Reynolds number Re", 1.0); 
00095 
00096 
00097 
00098 std::cout << " vecDim = " << vecDim << std::endl;
00099 std::cout << " numDims = " << numDims << std::endl;
00100 
00101 
00102 if (vecDim != numDims+2) {TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00103                                   std::endl << "Error in PHAL::ComprNS constructor:  " <<
00104                                   "Invalid Parameter vecDim.  vecDim should be numDims + 2 = " << numDims + 2 << "." << std::endl);}  
00105 
00106 
00107 }
00108 
00109 //**********************************************************************
00110 template<typename EvalT, typename Traits>
00111 void ComprNSResid<EvalT, Traits>::
00112 postRegistrationSetup(typename Traits::SetupData d,
00113                       PHX::FieldManager<Traits>& fm)
00114 {
00115   this->utils.setFieldData(qFluct,fm);
00116   this->utils.setFieldData(qFluctGrad,fm);
00117   this->utils.setFieldData(qFluctDot,fm);
00118   this->utils.setFieldData(force,fm);
00119   this->utils.setFieldData(mu,fm);
00120   this->utils.setFieldData(kappa,fm);
00121   this->utils.setFieldData(lambda,fm);
00122   this->utils.setFieldData(wBF,fm);
00123   this->utils.setFieldData(wGradBF,fm);
00124   this->utils.setFieldData(tau11,fm);
00125   this->utils.setFieldData(tau12,fm);
00126   this->utils.setFieldData(tau13,fm);
00127   this->utils.setFieldData(tau22,fm);
00128   this->utils.setFieldData(tau23,fm);
00129   this->utils.setFieldData(tau33,fm);
00130 
00131   this->utils.setFieldData(Residual,fm);
00132 }
00133 
00134 //**********************************************************************
00135 template<typename EvalT, typename Traits>
00136 void ComprNSResid<EvalT, Traits>::
00137 evaluateFields(typename Traits::EvalData workset)
00138 {
00139   typedef Intrepid::FunctionSpaceTools FST;
00140 
00141   if (numDims == 2) { //2D case; order of variables is rho, u, v, T
00142       for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00143         for (std::size_t node=0; node < numNodes; ++node) {
00144           for (std::size_t i=0; i<vecDim; i++) 
00145              Residual(cell,node,i) = 0.0; 
00146           for (std::size_t qp=0; qp < numQPs; ++qp) {
00147              Residual(cell,node,0) = qFluctDot(cell,qp,0)*wBF(cell,node,qp);  //d(rho)/dt
00148              for (std::size_t i=1; i < vecDim; i++) {
00149                 Residual(cell,node,i) = qFluct(cell,qp,0)*qFluctDot(cell,qp,i)*wBF(cell,node,qp); //rho*du_i/dt; rho*dT/dt 
00150              }
00151           }
00152           for (std::size_t qp=0; qp < numQPs; ++qp) {
00153              Residual(cell, node, 0) += qFluct(cell,qp,0)*(qFluctGrad(cell,qp,1,0)+qFluctGrad(cell,qp,2,1))*wBF(cell,node,qp) //rho*div(u)
00154                                       + (qFluct(cell,qp,1)*qFluctGrad(cell,qp,0,0) + qFluct(cell,qp,2)*qFluctGrad(cell,qp,0,1))*wBF(cell,node,qp) //u*d(rho)/dx + v*d(rho)/dy  
00155                                       + force(cell,qp,0)*wBF(cell,node,qp); //f0
00156              Residual(cell, node, 1) += qFluct(cell,qp,0)*(qFluct(cell,qp,1)*qFluctGrad(cell,qp,1,0) + qFluct(cell,qp,2)*qFluctGrad(cell,qp,1,1))*wBF(cell,node,qp) //rho*(u*du/dx + v*du/dy)
00157                                       + Rgas*(qFluct(cell,qp,0)*qFluctGrad(cell,qp,3,0) + qFluct(cell,qp,3)*qFluctGrad(cell,qp,0,0))*wBF(cell,node,qp) //R*(rho*dT/dx + T*d(rho)/dx) 
00158                                       + 1.0/Re*tau11(cell,qp)*wGradBF(cell,node,qp,0) //tau11
00159                                       + 1.0/Re*tau12(cell,qp)*wGradBF(cell,node,qp,1)//tau12
00160                                       + force(cell,qp,1)*wBF(cell,node,qp); //f1
00161              Residual(cell, node, 2) += qFluct(cell,qp,0)*(qFluct(cell,qp,1)*qFluctGrad(cell,qp,2,0) + qFluct(cell,qp,2)*qFluctGrad(cell,qp,2,1))*wBF(cell,node,qp) //rho*(u*dv/dx + v*dv/dy)
00162                                       + Rgas*(qFluct(cell,qp,0)*qFluctGrad(cell,qp,3,1) + qFluct(cell,qp,3)*qFluctGrad(cell,qp,0,1))*wBF(cell,node,qp) //R*(rho*dT/dy + T*d(rho)/dy)
00163                                       + 1.0/Re*tau12(cell,qp)*wGradBF(cell,node,qp,0) //tau21
00164                                       + 1.0/Re*tau22(cell,qp)*wGradBF(cell,node,qp,1) //tau22
00165                                       + force(cell,qp,2)*wBF(cell,node,qp); //f2
00166              Residual(cell, node, 3) += qFluct(cell,qp,0)*(qFluct(cell,qp,1)*qFluctGrad(cell,qp,3,0) + qFluct(cell,qp,2)*qFluctGrad(cell,qp,3,1) //rho*(u*dT/dx + v*dT/dy) 
00167                                       + (gamma_gas - 1.0)*qFluct(cell,qp,3)*(qFluctGrad(cell,qp,1,0) + qFluctGrad(cell,qp,2,1)))*wBF(cell,node,qp) //(gamma-1)*T*div(u)
00168                                       - (gamma_gas - 1.0)/Rgas*qFluctGrad(cell,qp,1,0)*1.0/Re*tau11(cell,qp)*wBF(cell,node,qp) //-(gamma-1)/R*du/dx*tau11 
00169                                       - (gamma_gas - 1.0)/Rgas*1.0/Re*qFluctGrad(cell,qp,2,0)*tau12(cell,qp)*wBF(cell,node,qp) //-(gamma-1)/R*(dv/dx*tau12)
00170                                       - (gamma_gas - 1.0)/Rgas*1.0/Re*qFluctGrad(cell,qp,1,1)*tau12(cell,qp)*wBF(cell,node,qp) //-(gamma-1)/R*(du/dy*tau12)
00171                                       - (gamma_gas - 1.0)/Rgas*qFluctGrad(cell,qp,2,1)*1.0/Re*tau22(cell,qp)*wBF(cell,node,qp) // -(gamma-1)/R*dv/dy*tau22
00172                                       + gamma_gas*kappa(cell,qp)/(Pr*Re)*(qFluctGrad(cell,qp,3,0)*wGradBF(cell,node,qp,0) + qFluctGrad(cell,qp,3,1)*wGradBF(cell,node,qp,1)) //gamma*kappa/(Pr*Re)*(Delta T)
00173                                       + force(cell,qp,3)*wBF(cell,node,qp);  //f3 
00174            } 
00175           } 
00176         }
00177      }
00178    else if (numDims == 3) { //3D case - TO IMPLEMENT
00179       for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00180         for (std::size_t node=0; node < numNodes; ++node) {
00181           for (std::size_t i=0; i<vecDim; i++) 
00182              Residual(cell,node,i) = 0.0; 
00183           for (std::size_t qp=0; qp < numQPs; ++qp) {
00184              for (std::size_t i=0; i < vecDim; i++) {
00185                 Residual(cell,node,i) = qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
00186              }
00187           }
00188           for (std::size_t qp=0; qp < numQPs; ++qp) {
00189              Residual(cell, node, 0) += 0.0;  
00190              Residual(cell, node, 1) += 0.0;  
00191              Residual(cell, node, 2) += 0.0;  
00192              Residual(cell, node, 3) += 0.0;  
00193              Residual(cell, node, 4) += 0.0;  
00194             } 
00195           } 
00196         }
00197      }
00198 }
00199 
00200 //**********************************************************************
00201 }
00202
PHAL_ComprNSResid_Def.hpp
