Albany: a Trilinos-based PDE code: PHAL_LinComprNSResid

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 LinComprNSResid<EvalT, Traits>::
00019 LinComprNSResid(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   Residual   (p.get<std::string>                   ("Residual Name"),
00033               p.get<Teuchos::RCP<PHX::DataLayout> >("Node Vector Data Layout") ) 
00034 {
00035 
00036 
00037 
00038   this->addDependentField(qFluct);
00039   this->addDependentField(qFluctGrad);
00040   this->addDependentField(qFluctDot);
00041   this->addDependentField(force);
00042   this->addDependentField(wBF);
00043   this->addDependentField(wGradBF);
00044 
00045   this->addEvaluatedField(Residual);
00046 
00047 
00048   this->setName("LinComprNSResid"+PHX::TypeString<EvalT>::value);
00049 
00050   std::vector<PHX::DataLayout::size_type> dims;
00051   wGradBF.fieldTag().dataLayout().dimensions(dims);
00052   numNodes = dims[1];
00053   numQPs   = dims[2];
00054   numDims  = dims[3];
00055 
00056 
00057   Teuchos::ParameterList* bf_list =
00058   p.get<Teuchos::ParameterList*>("Parameter List");
00059   std::string eqnType = bf_list->get("Type", "Euler");
00060   
00061   if (eqnType == "Euler") {
00062     std::cout << "setting euler equations!" << std::endl; 
00063     eqn_type = EULER; 
00064   }
00065   else if (eqnType == "Navier-Stokes") {
00066     std::cout << "setting n-s equations!" << std::endl; 
00067     eqn_type = NS; 
00068   }
00069 
00070 
00071   qFluct.fieldTag().dataLayout().dimensions(dims);
00072   vecDim  = dims[2];
00073 
00074   Teuchos::Array<double> defaultBaseFlowData(numDims+2);  
00075   baseFlowData = bf_list->get("Base Flow Data", defaultBaseFlowData); 
00076   //for EULER, baseFlowData = (ubar, vbar, wbar, zetabar, pbar)
00077   //for NS, baseFlowData = (ubar, vbar, wbar, Tbar, rhobar)
00078 
00079   gamma_gas = bf_list->get("Gamma", 1.4); 
00080   Rgas = bf_list->get("Gas constant R", 0.714285733);
00081   Pr = bf_list->get("Prandtl number Pr", 0.72); 
00082   Re = bf_list->get("Reynolds number Re", 1.0); 
00083   mu = bf_list->get("Viscocity mu", 0.0); 
00084   lambda = -2.0/3.0*mu; //Stokes' hypothesis
00085   kappa = bf_list->get("Diffusivity kappa", 0.0);  
00086   IBP_convect_terms = bf_list->get("IBP Convective Terms", false); 
00087 
00088   if (IBP_convect_terms == true)
00089     std::cout  << "Integrating convective terms by parts in weak form." << std::endl; 
00090 
00091 
00092 std::cout << " vecDim = " << vecDim << std::endl;
00093 std::cout << " numDims = " << numDims << std::endl;
00094 
00095 if (baseFlowData.size()!=numDims+2) {TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00096                                   std::endl << "Error in PHAL::LinComprNS constructor:  " <<
00097                                   "baseFlow data should have length numDims + 2 =  " << numDims+2 << "." << std::endl);} 
00098 
00099 
00100 if (eqn_type == EULER & vecDim != numDims+1) {TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00101                                   std::endl << "Error in PHAL::LinComprNS constructor:  " <<
00102                                   "Invalid Parameter vecDim.  vecDim should be numDims + 1 = " << numDims + 1 << " for Euler equations." << std::endl);}  
00103 
00104 if (eqn_type == NS & vecDim != numDims+2) {TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00105                                   std::endl << "Error in PHAL::LinComprNS constructor:  " <<
00106                                   "Invalid Parameter vecDim.  vecDim should be numDims + 2 = " << numDims + 2 << " for Navier-Stokes equations." << std::endl);}  
00107 
00108 }
00109 
00110 //**********************************************************************
00111 template<typename EvalT, typename Traits>
00112 void LinComprNSResid<EvalT, Traits>::
00113 postRegistrationSetup(typename Traits::SetupData d,
00114                       PHX::FieldManager<Traits>& fm)
00115 {
00116   this->utils.setFieldData(qFluct,fm);
00117   this->utils.setFieldData(qFluctGrad,fm);
00118   this->utils.setFieldData(qFluctDot,fm);
00119   this->utils.setFieldData(force,fm);
00120   this->utils.setFieldData(wBF,fm);
00121   this->utils.setFieldData(wGradBF,fm);
00122 
00123   this->utils.setFieldData(Residual,fm);
00124 }
00125 
00126 //**********************************************************************
00127 template<typename EvalT, typename Traits>
00128 void LinComprNSResid<EvalT, Traits>::
00129 evaluateFields(typename Traits::EvalData workset)
00130 {
00131   typedef Intrepid::FunctionSpaceTools FST;
00132 
00133   if (eqn_type == EULER) { //Euler equations
00134    if (numDims == 1) { //1D case
00135     double ubar = baseFlowData[0];
00136     double zetabar = baseFlowData[1]; 
00137     double pbar = baseFlowData[2];
00138     if (IBP_convect_terms == false) {//variational formulation in which the convective terms are not integrated by parts
00139       for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00140         for (std::size_t node=0; node < numNodes; ++node) {
00141           for (std::size_t i=0; i<vecDim; i++) 
00142              Residual(cell,node,i) = 0.0; 
00143           for (std::size_t qp=0; qp < numQPs; ++qp) {
00144              for (std::size_t i=0; i < vecDim; i++) {
00145                 Residual(cell,node,i) = qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
00146              }
00147           }
00148           for (std::size_t qp=0; qp < numQPs; ++qp) {
00149              Residual(cell, node, 0) += ubar*qFluctGrad(cell,qp,0,0)*wBF(cell,node,qp)  
00150                                      + zetabar*qFluctGrad(cell,qp,1,0)*wBF(cell,node,qp)  
00151                                      + force(cell,qp,0)*wBF(cell,node,qp); //ubar*du'/dx + zetabar*dp'/dx + f0
00152              Residual(cell, node, 1) += gamma_gas*pbar*qFluctGrad(cell,qp,0,0)*wBF(cell,node,qp) 
00153                                      + ubar*qFluctGrad(cell,qp,1,0)*wBF(cell,node,qp)  
00154                                      + force(cell,qp,1)*wBF(cell,node,qp); //gamma*pbar*du'/dx + ubar*dp'/dx + f2
00155              
00156             } 
00157           } 
00158         }
00159      }
00160      else { //variational formulation in which the convective terms are integrated by parts
00161        for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00162          for (std::size_t node=0; node < numNodes; ++node) {
00163            for (std::size_t i=0; i<vecDim; i++) 
00164              Residual(cell,node,i) = 0.0; 
00165            for (std::size_t qp=0; qp < numQPs; ++qp) {
00166               for (std::size_t i=0; i < vecDim; i++) {
00167                  Residual(cell,node,i) = qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
00168               }
00169            }
00170            for (std::size_t qp=0; qp < numQPs; ++qp) {
00171               Residual(cell, node, 0) += -1.0*ubar*qFluct(cell,qp,0)*wGradBF(cell,node,qp,0)  
00172                                       - zetabar*qFluct(cell,qp,1)*wGradBF(cell,node,qp,0) 
00173                                       + force(cell,qp,0)*wBF(cell,node,qp); //ubar*du'/dx + zetabar*dp'/dx + f0
00174               Residual(cell, node, 1) += -1.0*gamma_gas*pbar*qFluct(cell,qp,0)*wGradBF(cell,node,qp,0)  
00175                                       - ubar*qFluct(cell,qp,1)*wGradBF(cell,node,qp,0)  
00176                                       + force(cell,qp,1)*wBF(cell,node,qp); //gamma*pbar*du'/dx + ubar*dp'/dx  + f2
00177             } 
00178           } 
00179         }
00180      }
00181     }
00182    if (numDims == 2) { //2D case
00183     double ubar = baseFlowData[0]; 
00184     double vbar = baseFlowData[1]; 
00185     double zetabar = baseFlowData[2]; 
00186     double pbar = baseFlowData[3];
00187     if (IBP_convect_terms == false) {//variational formulation in which the convective terms are not integrated by parts
00188       for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00189         for (std::size_t node=0; node < numNodes; ++node) {
00190           for (std::size_t i=0; i<vecDim; i++) 
00191              Residual(cell,node,i) = 0.0; 
00192           for (std::size_t qp=0; qp < numQPs; ++qp) {
00193              for (std::size_t i=0; i < vecDim; i++) {
00194                 Residual(cell,node,i) = qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
00195              }
00196           }
00197           for (std::size_t qp=0; qp < numQPs; ++qp) {
00198              Residual(cell, node, 0) += ubar*qFluctGrad(cell,qp,0,0)*wBF(cell,node,qp) + vbar*qFluctGrad(cell,qp,0,1)*wBF(cell,node,qp) 
00199                                      + zetabar*qFluctGrad(cell,qp,2,0)*wBF(cell,node,qp) 
00200                                      + force(cell,qp,0)*wBF(cell,node,qp); //ubar*du'/dx + vbar*du'/dy + zetabar*dp'/dx + f0
00201              Residual(cell, node, 1) += ubar*qFluctGrad(cell,qp,1,0)*wBF(cell,node,qp) + vbar*qFluctGrad(cell,qp,1,1)*wBF(cell,node,qp) 
00202                                      + zetabar*qFluctGrad(cell,qp,2,1)*wBF(cell,node,qp) 
00203                                      + force(cell,qp,1)*wBF(cell,node,qp); //ubar*dv'/dx + vbar*dv'/dy + zetabar*dp'/dy + f1
00204              Residual(cell, node, 2) += gamma_gas*pbar*(qFluctGrad(cell,qp,0,0) + qFluctGrad(cell,qp,1,1))*wBF(cell,node,qp) 
00205                                      + ubar*qFluctGrad(cell,qp,2,0)*wBF(cell,node,qp) + vbar*qFluctGrad(cell,qp,2,1)*wBF(cell,node,qp) 
00206                                      + force(cell,qp,2)*wBF(cell,node,qp); //gamma*pbar*div(u') + ubar*dp'/dx + vbar*dp'/dy + f2
00207             } 
00208           } 
00209         }
00210      }
00211      else { //variational formulation in which the convective terms are integrated by parts
00212        for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00213          for (std::size_t node=0; node < numNodes; ++node) {
00214            for (std::size_t i=0; i<vecDim; i++) 
00215              Residual(cell,node,i) = 0.0; 
00216            for (std::size_t qp=0; qp < numQPs; ++qp) {
00217               for (std::size_t i=0; i < vecDim; i++) {
00218                  Residual(cell,node,i) = qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
00219               }
00220            }
00221            for (std::size_t qp=0; qp < numQPs; ++qp) {
00222               Residual(cell, node, 0) += -1.0*ubar*qFluct(cell,qp,0)*wGradBF(cell,node,qp,0) - vbar*qFluct(cell,qp,0)*wGradBF(cell,node,qp,1) 
00223                                       - zetabar*qFluct(cell,qp,2)*wGradBF(cell,node,qp,0) 
00224                                       + force(cell,qp,0)*wBF(cell,node,qp); //ubar*du'/dx + vbar*du'/dy + zetabar*dp'/dx + f0
00225               Residual(cell, node, 1) += -1.0*ubar*qFluct(cell,qp,1)*wGradBF(cell,node,qp,0) - vbar*qFluct(cell,qp,1)*wGradBF(cell,node,qp,1) 
00226                                       - zetabar*qFluct(cell,qp,2)*wGradBF(cell,node,qp,1) 
00227                                       + force(cell,qp,1)*wBF(cell,node,qp); //ubar*dv'/dx + vbar*dv'/dy + zetabar*dp'/dy + f1
00228               Residual(cell, node, 2) += -1.0*gamma_gas*pbar*(qFluct(cell,qp,0)*wGradBF(cell,node,qp,0) + qFluct(cell,qp,1)*wGradBF(cell,node,qp,1)) 
00229                                       - ubar*qFluct(cell,qp,2)*wGradBF(cell,node,qp,0) - vbar*qFluct(cell,qp,2)*wGradBF(cell,node,qp,1) 
00230                                       + force(cell,qp,2)*wBF(cell,node,qp); //gamma*pbar*div(u') + ubar*dp'/dx + vbar*dp'/dy + f2
00231             } 
00232           } 
00233         }
00234      }
00235     }
00236    else if (numDims == 3) { //3D case
00237     double ubar = baseFlowData[0]; 
00238     double vbar = baseFlowData[1]; 
00239     double wbar = baseFlowData[2]; 
00240     double zetabar = baseFlowData[3]; 
00241     double pbar = baseFlowData[4];
00242     if (IBP_convect_terms == false) {//variational formulation in which the convective terms are not integrated by parts
00243       for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00244         for (std::size_t node=0; node < numNodes; ++node) {
00245           for (std::size_t i=0; i<vecDim; i++) 
00246              Residual(cell,node,i) = 0.0; 
00247           for (std::size_t qp=0; qp < numQPs; ++qp) {
00248              for (std::size_t i=0; i < vecDim; i++) {
00249                 Residual(cell,node,i) = qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
00250              }
00251           }
00252           for (std::size_t qp=0; qp < numQPs; ++qp) {
00253              Residual(cell, node, 0) += ubar*qFluctGrad(cell,qp,0,0)*wBF(cell,node,qp) + vbar*qFluctGrad(cell,qp,0,1)*wBF(cell,node,qp) 
00254                                      +wbar*qFluctGrad(cell,qp,0,2)*wBF(cell,node,qp) + zetabar*qFluctGrad(cell,qp,3,0)*wBF(cell,node,qp) 
00255                                      + force(cell,qp,0)*wBF(cell,node,qp); //ubar*du'/dx + vbar*du'/dy + wbar*du'/dz + zetabar*dp'/dx + f0
00256              Residual(cell, node, 1) += ubar*qFluctGrad(cell,qp,1,0)*wBF(cell,node,qp) + vbar*qFluctGrad(cell,qp,1,1)*wBF(cell,node,qp) 
00257                                      + wbar*qFluctGrad(cell,qp,1,2)*wBF(cell,node,qp) + zetabar*qFluctGrad(cell,qp,3,1)*wBF(cell,node,qp) 
00258                                      + force(cell,qp,1)*wBF(cell,node,qp); //ubar*dv'/dx + vbar*dv'/dy + wbar*dv'/dz + zetabar*dp'/dy + f1
00259              Residual(cell, node, 2) += ubar*qFluctGrad(cell,qp,2,0)*wBF(cell,node,qp) + vbar*qFluctGrad(cell,qp,2,1)*wBF(cell,node,qp)
00260                                      + wbar*qFluctGrad(cell,qp,2,2)*wBF(cell,node,qp) + zetabar*qFluctGrad(cell,qp,3,2)*wBF(cell,node,qp)
00261                                      + force(cell,qp,2)*wBF(cell,node,qp); //ubar*dw'/dx + vbar*dw'/dy + wbar*dw'/dz + zetabar*dp'/dz + f2
00262              Residual(cell, node, 3) += gamma_gas*pbar*(qFluctGrad(cell,qp,0,0) + qFluctGrad(cell,qp,1,1) + qFluctGrad(cell,qp,2,2))*wBF(cell,node,qp) 
00263                                      + ubar*qFluctGrad(cell,qp,3,0)*wBF(cell,node,qp) + vbar*qFluctGrad(cell,qp,3,1)*wBF(cell,node,qp) + wbar*qFluctGrad(cell,qp,3,2)*wBF(cell,node,qp) 
00264                                      + force(cell,qp,3)*wBF(cell,node,qp); //gamma*pbar*div(u') + ubar*dp'/dx + vbar*dp'/dy + wbar*dp'/dz + f3
00265             } 
00266           } 
00267         }
00268      }
00269      else { //variational formulation in which the convective terms are integrated by parts
00270       for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00271         for (std::size_t node=0; node < numNodes; ++node) {
00272           for (std::size_t i=0; i<vecDim; i++) 
00273              Residual(cell,node,i) = 0.0; 
00274           for (std::size_t qp=0; qp < numQPs; ++qp) {
00275              for (std::size_t i=0; i < vecDim; i++) {
00276                 Residual(cell,node,i) = qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
00277              }
00278           }
00279           for (std::size_t qp=0; qp < numQPs; ++qp) {
00280              Residual(cell, node, 0) += -1.0*ubar*qFluct(cell,qp,0)*wGradBF(cell,node,qp,0) - vbar*qFluct(cell,qp,0)*wGradBF(cell,node,qp,1) 
00281                                      - wbar*qFluct(cell,qp,0)*wGradBF(cell,node,qp,2) - zetabar*qFluct(cell,qp,3)*wGradBF(cell,node,qp,0) 
00282                                      + force(cell,qp,0)*wBF(cell,node,qp); //ubar*du'/dx + vbar*du'/dy + wbar*du'/dz + zetabar*dp'/dx + f0
00283              Residual(cell, node, 1) += -1.0*ubar*qFluct(cell,qp,1)*wGradBF(cell,node,qp,0) - vbar*qFluct(cell,qp,1)*wGradBF(cell,node,qp,1) 
00284                                      - wbar*qFluct(cell,qp,1)*wGradBF(cell,node,qp,2) - zetabar*qFluct(cell,qp,3)*wGradBF(cell,node,qp,1) 
00285                                      + force(cell,qp,1)*wBF(cell,node,qp); //ubar*dv'/dx + vbar*dv'/dy + wbar*dv'/dz + zetabar*dp'/dy + f1
00286              Residual(cell, node, 2) += -1.0*ubar*qFluct(cell,qp,2)*wGradBF(cell,node,qp,0) - vbar*qFluct(cell,qp,2)*wGradBF(cell,node,qp,1)
00287                                      - wbar*qFluct(cell,qp,2)*wGradBF(cell,node,qp,2) - zetabar*qFluct(cell,qp,3)*wGradBF(cell,node,qp,2)
00288                                      + force(cell,qp,2)*wBF(cell,node,qp); //ubar*dw'/dx + vbar*dw'/dy + wbar*dw'/dz + zetabar*dp'/dz + f2
00289              Residual(cell, node, 3) += -1.0*gamma_gas*pbar*(qFluct(cell,qp,0)*wGradBF(cell,node,qp,0) + qFluct(cell,qp,1)*wGradBF(cell,node,qp,1) + qFluct(cell,qp,2)*wGradBF(cell,node,qp,2)) 
00290                                      - ubar*qFluct(cell,qp,3)*wGradBF(cell,node,qp,0) - vbar*qFluct(cell,qp,3)*wGradBF(cell,node,qp,1) - wbar*qFluct(cell,qp,3)*wGradBF(cell,node,qp,2) 
00291                                      + force(cell,qp,3)*wBF(cell,node,qp); //gamma*pbar*div(u') + ubar*dp'/dx + vbar*dp'/dy + wbar*dp'/dz + f3
00292             } 
00293           } 
00294         }
00295      }
00296    }
00297   }
00298   else if (eqn_type == NS) { //Navier-Stokes equations
00299     if (numDims == 2) { //2D case
00300       double ubar = baseFlowData[0]; 
00301       double vbar = baseFlowData[1]; 
00302       double Tbar = baseFlowData[2]; 
00303       double rhobar = baseFlowData[3];
00304       for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00305         for (std::size_t node=0; node < numNodes; ++node) {
00306           for (std::size_t i=0; i<vecDim; i++) 
00307              Residual(cell,node,i) = 0.0; 
00308           for (std::size_t qp=0; qp < numQPs; ++qp) {
00309              for (std::size_t i=0; i < vecDim; i++) {
00310                 Residual(cell,node,i) = rhobar*qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
00311                 if (i == vecDim-1)
00312                   Residual(cell,node,i) = qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
00313              }
00314           }
00315           for (std::size_t qp=0; qp < numQPs; ++qp) {
00316              Residual(cell, node, 0) += rhobar*ubar*qFluctGrad(cell,qp,0,0)*wBF(cell,node,qp) + rhobar*vbar*qFluctGrad(cell,qp,0,1)*wBF(cell,node,qp) 
00317                                      + rhobar*Rgas*qFluctGrad(cell,qp,2,0)*wBF(cell,node,qp) + Rgas*Tbar*qFluctGrad(cell,qp,3,0)*wBF(cell,node,qp) 
00318                                      + 1.0/Re*((2.0*mu + lambda)*qFluctGrad(cell,qp,0,0)*wGradBF(cell,node,qp,0) + lambda*qFluctGrad(cell,qp,1,1)*wGradBF(cell,node,qp,0) 
00319                                                + mu*qFluctGrad(cell,qp,1,0)*wGradBF(cell,node,qp,1) + mu*qFluctGrad(cell,qp,0,1)*wGradBF(cell,node,qp,1))
00320                                      + force(cell,qp,0)*wBF(cell,node,qp); //rhobar*ubar*du'/dx + rhobar*vbar*du'/dy + rhobar*R*dT'/dx + R*Tbar*drho'/dx + 
00321                                                                            //1/Re*(d/dx((2*mu+lambda)*du'/dx + lambda*dv'/dy) + d/dy*(mu*dv'/dx + mu*du'/dy)) + 
00322                                                                            //f0
00323              Residual(cell, node, 1) += rhobar*ubar*qFluctGrad(cell,qp,1,0)*wBF(cell,node,qp) + rhobar*vbar*qFluctGrad(cell,qp,1,1)*wBF(cell,node,qp) 
00324                                      + rhobar*Rgas*qFluctGrad(cell,qp,2,1)*wBF(cell,node,qp) + Rgas*Tbar*qFluctGrad(cell,qp,3,1)*wBF(cell,node,qp) 
00325                                      + 1.0/Re*(mu*qFluctGrad(cell,qp,1,0)*wGradBF(cell,node,qp,0) + mu*qFluctGrad(cell,qp,0,1)*wGradBF(cell,node,qp,0)
00326                                                + lambda*qFluctGrad(cell,qp,0,0)*wGradBF(cell,node,qp,1) + (2.0*mu + lambda)*qFluctGrad(cell,qp,1,1)*wGradBF(cell,node,qp,1))
00327                                      + force(cell,qp,1)*wBF(cell,node,qp); //rhobar*ubar*dv'/dx + rhobar*vbar*dv'/dy + rhobar*R*dT'/dy + R*Tbar*drho'/dy + 
00328                                                                            //1/Re*(d/dx(mu*dv'/dx + mu*du'/dy) + d/dy(lambda*du'/dx + (2*mu + lambda)*dv'/dy)) +
00329                                                                            //f1
00330              Residual(cell, node, 2) += rhobar*Tbar*(gamma_gas - 1.0)*(qFluctGrad(cell,qp,0,0) + qFluctGrad(cell,qp,1,1))*wBF(cell,node,qp) 
00331                                      + rhobar*ubar*qFluctGrad(cell,qp,2,0)*wBF(cell,node,qp) + rhobar*vbar*qFluctGrad(cell,qp,2,1)*wBF(cell,node,qp) 
00332                                      + 1.0/Re*(gamma_gas*kappa/Pr*(qFluctGrad(cell,qp,2,0)*wGradBF(cell,node,qp,0) + qFluctGrad(cell,qp,2,1)*wGradBF(cell,node,qp,1)))
00333                                      + force(cell,qp,2)*wBF(cell,node,qp); //rhobar*Tbar*(gamma-1)*div(u') + rhobar*ubar*dT'/dx + rhobar*vbar*dT'/dy +
00334                                                                            //1/Re*(d/dx(gamma*kappa/Pr*dT'/dx) + d/dy(gamma*kappa/Pr*dT'/dy) +
00335                                                                            //f2 
00336              Residual(cell, node, 3) += rhobar*(qFluctGrad(cell,qp,0,0) + qFluctGrad(cell,qp,1,1))*wBF(cell,node,qp) 
00337                                      + ubar*qFluctGrad(cell,qp,3,0)*wBF(cell,node,qp) + vbar*qFluctGrad(cell,qp,3,1)*wBF(cell,node,qp) 
00338                                      + force(cell,qp,3)*wBF(cell,node,qp); //rhobar*div(u') + ubar*drho'/dx + vbar*drho'/dy + f3 
00339                                      
00340             } 
00341           } 
00342         }
00343     }
00344     else if (numDims == 3) { //3D case
00345       double ubar = baseFlowData[0]; 
00346       double vbar = baseFlowData[1]; 
00347       double wbar = baseFlowData[2]; 
00348       double Tbar = baseFlowData[3]; 
00349       double rhobar = baseFlowData[4];
00350       for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00351         for (std::size_t node=0; node < numNodes; ++node) {
00352           for (std::size_t i=0; i<vecDim; i++) 
00353              Residual(cell,node,i) = 0.0; 
00354           for (std::size_t qp=0; qp < numQPs; ++qp) {
00355              for (std::size_t i=0; i < vecDim; i++) {
00356                 Residual(cell,node,i) = rhobar*qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
00357                 if (i == vecDim-1)
00358                   Residual(cell,node,i) = qFluctDot(cell,qp,i)*wBF(cell,node,qp); 
00359              }
00360           }
00361           for (std::size_t qp=0; qp < numQPs; ++qp) {
00362              Residual(cell, node, 0) += rhobar*ubar*qFluctGrad(cell,qp,0,0)*wBF(cell,node,qp) + rhobar*vbar*qFluctGrad(cell,qp,0,1)*wBF(cell,node,qp)
00363                                      + rhobar*wbar*qFluctGrad(cell,qp,0,2)*wBF(cell,node,qp) 
00364                                      + rhobar*Rgas*qFluctGrad(cell,qp,3,0)*wBF(cell,node,qp) + Rgas*Tbar*qFluctGrad(cell,qp,4,0)*wBF(cell,node,qp) 
00365                                      + 1.0/Re*((2.0*mu + lambda)*qFluctGrad(cell,qp,0,0)*wGradBF(cell,node,qp,0) + lambda*qFluctGrad(cell,qp,1,1)*wGradBF(cell,node,qp,0)
00366                                                + lambda*qFluctGrad(cell,qp,2,2)*wGradBF(cell,node,qp,0) 
00367                                                + mu*qFluctGrad(cell,qp,1,0)*wGradBF(cell,node,qp,1) + mu*qFluctGrad(cell,qp,0,1)*wGradBF(cell,node,qp,1)
00368                                                + mu*qFluctGrad(cell,qp,2,0)*wGradBF(cell,node,qp,2) + mu*qFluctGrad(cell,qp,0,2)*wGradBF(cell,node,qp,2))
00369                                      + force(cell,qp,0)*wBF(cell,node,qp); //rhobar*ubar*du'/dx + rhobar*vbar*du'/dy + rhobar*wbar*du'/dz + rhobar*R*dT'/dx + R*Tbar*drho'/dx + 
00370                                                                            //1/Re*(d/dx((2*mu+lambda)*du'/dx + lambda*dv'/dy + lambda*dw'/dz) + d/dy*(mu*dv'/dx + mu*du'/dy)) + 
00371                                                                            //d/dz(mu*dw'/dx + mu*du'/dz)) + f0
00372              Residual(cell, node, 1) += rhobar*ubar*qFluctGrad(cell,qp,1,0)*wBF(cell,node,qp) + rhobar*vbar*qFluctGrad(cell,qp,1,1)*wBF(cell,node,qp)
00373                                      + rhobar*wbar*qFluctGrad(cell,qp,1,2)*wBF(cell,node,qp)  
00374                                      + rhobar*Rgas*qFluctGrad(cell,qp,3,1)*wBF(cell,node,qp) + Rgas*Tbar*qFluctGrad(cell,qp,4,1)*wBF(cell,node,qp) 
00375                                      + 1.0/Re*(mu*qFluctGrad(cell,qp,1,0)*wGradBF(cell,node,qp,0) + mu*qFluctGrad(cell,qp,0,1)*wGradBF(cell,node,qp,0)
00376                                                + lambda*qFluctGrad(cell,qp,0,0)*wGradBF(cell,node,qp,1) + (2.0*mu + lambda)*qFluctGrad(cell,qp,1,1)*wGradBF(cell,node,qp,1)
00377                                                + lambda*qFluctGrad(cell,qp,2,2)*wGradBF(cell,node,qp,1) 
00378                                                + mu*qFluctGrad(cell,qp,2,1)*wGradBF(cell,node,qp,2) + mu*qFluctGrad(cell,qp,1,2)*wGradBF(cell,node,qp,2))
00379                                      + force(cell,qp,1)*wBF(cell,node,qp); //rhobar*ubar*dv'/dx + rhobar*vbar*dv'/dy + rhobar*wbar*dv'/dz + rhobar*R*dT'/dy + R*Tbar*drho'/dy + 
00380                                                                            //1/Re*(d/dx(mu*dv'/dx + mu*du'/dy) + d/dy(lambda*du'/dx + (2*mu + lambda)*dv'/dy + lambda*dw'/dz) +
00381                                                                            //d/dz(mu*dw'/dy + mu*dv'/dz)) + f1
00382              Residual(cell, node, 2) += rhobar*ubar*qFluctGrad(cell,qp,2,0)*wBF(cell,node,qp) + rhobar*vbar*qFluctGrad(cell,qp,2,1)*wBF(cell,node,qp) 
00383                                      + rhobar*wbar*qFluctGrad(cell,qp,2,2)*wBF(cell,node,qp) 
00384                                      + rhobar*Rgas*qFluctGrad(cell,qp,3,2)*wBF(cell,node,qp) + Rgas*Tbar*qFluctGrad(cell,qp,4,2)*wBF(cell,node,qp)
00385                                      + 1.0/Re*(mu*qFluctGrad(cell,qp,2,0)*wGradBF(cell,node,qp,0) + mu*qFluctGrad(cell,qp,0,2)*wGradBF(cell,node,qp,0) 
00386                                                + mu*qFluctGrad(cell,qp,2,1)*wGradBF(cell,node,qp,1) + mu*qFluctGrad(cell,qp,1,2)*wGradBF(cell,node,qp,1) 
00387                                                + lambda*qFluctGrad(cell,qp,0,0)*wGradBF(cell,node,qp,2) + lambda*qFluctGrad(cell,qp,1,1)*wGradBF(cell,node,qp,2) 
00388                                                + (2.0*mu + lambda)*qFluctGrad(cell,qp,2,2)*wGradBF(cell,node,qp,2)) 
00389                                      + force(cell,qp,2)*wBF(cell,node,qp); //rhobar*ubar*dw'/dx + rhobar*vbar*dw'/dy + rhobar*wbar*dw'/dz + rhobar*R*dT'/dz + R*Tbar*drho'/dz + 
00390                                                                            //1/Re*(d/dx(mu*dw'/dx + mu*du'/dz) + d/dy(mu*dw'/dy + mu*dv'/dz) + 
00391                                                                            //d/dz(lambda*du'/dx + lambda*dv'/dy + (2*mu + lambda)*dw'/dz)) + f2 
00392              Residual(cell, node, 3) += rhobar*Tbar*(gamma_gas - 1.0)*(qFluctGrad(cell,qp,0,0) + qFluctGrad(cell,qp,1,1) + qFluctGrad(cell,qp,2,2))*wBF(cell,node,qp) 
00393                                      + rhobar*ubar*qFluctGrad(cell,qp,3,0)*wBF(cell,node,qp) + rhobar*vbar*qFluctGrad(cell,qp,3,1)*wBF(cell,node,qp) 
00394                                      + 1.0/Re*(gamma_gas*kappa/Pr*(qFluctGrad(cell,qp,3,0)*wGradBF(cell,node,qp,0) + qFluctGrad(cell,qp,3,1)*wGradBF(cell,node,qp,1) + qFluctGrad(cell,qp,3,2)*wGradBF(cell,node,qp,2)))
00395                                      + force(cell,qp,3)*wBF(cell,node,qp); //rhobar*Tbar*(gamma-1)*div(u') + rhobar*ubar*dT'/dx + rhobar*vbar*dT'/dy + rhobar*wbar*dT'/dz +
00396                                                                            //1/Re*(d/dx(gamma*kappa/Pr*dT'/dx) + d/dy(gamma*kappa/Pr*dT'/dy) + d/dz(gamma*kappa/Pr*dT'/dz +
00397                                                                            //f3
00398              Residual(cell, node, 4) += rhobar*(qFluctGrad(cell,qp,0,0) + qFluctGrad(cell,qp,1,1) + qFluctGrad(cell,qp,2,2))*wBF(cell,node,qp) 
00399                                      + ubar*qFluctGrad(cell,qp,4,0)*wBF(cell,node,qp) + vbar*qFluctGrad(cell,qp,4,1)*wBF(cell,node,qp) + wbar*qFluctGrad(cell,qp,4,2)*wBF(cell,node,qp)
00400                                      + force(cell,qp,4)*wBF(cell,node,qp); //rhobar*div(u') + ubar*drho'/dx + vbar*drho'/dy + wbar*drho'/dz + f4 
00401                                      
00402             } 
00403           } 
00404         }
00405     }
00406   }
00407 }
00408 
00409 //**********************************************************************
00410 }
00411
PHAL_LinComprNSResid_Def.hpp
