Albany: a Trilinos-based PDE code: PHAL_ComprNSViscosity

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.hpp"
00010 
00011 #include "Intrepid_FunctionSpaceTools.hpp"
00012 
00013 namespace PHAL {
00014 //**********************************************************************
00015 
00016 template<typename EvalT, typename Traits>
00017 ComprNSViscosity<EvalT, Traits>::
00018 ComprNSViscosity(const Teuchos::ParameterList& p) :
00019   qFluct       (p.get<std::string>                   ("QP Variable Name"),
00020          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
00021   qFluctGrad      (p.get<std::string>                   ("Gradient QP Variable Name"),
00022          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Tensor Data Layout") ),
00023   mu          (p.get<std::string>                   ("Viscosity Mu QP Variable Name"),
00024          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),  
00025   kappa       (p.get<std::string>                   ("Viscosity Kappa QP Variable Name"),
00026          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
00027   lambda       (p.get<std::string>                   ("Viscosity Lambda QP Variable Name"),
00028          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
00029   tau11          (p.get<std::string>                   ("Stress Tau11 QP Variable Name"),
00030          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),  
00031   tau12       (p.get<std::string>                   ("Stress Tau12 QP Variable Name"),
00032          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
00033   tau13       (p.get<std::string>                   ("Stress Tau13 QP Variable Name"),
00034          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
00035   tau22       (p.get<std::string>                   ("Stress Tau22 QP Variable Name"),
00036          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
00037   tau23       (p.get<std::string>                   ("Stress Tau23 QP Variable Name"),
00038          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ), 
00039   tau33       (p.get<std::string>                   ("Stress Tau33 QP Variable Name"),
00040          p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ) 
00041 {
00042   Teuchos::ParameterList* visc_list = 
00043    p.get<Teuchos::ParameterList*>("Parameter List");
00044 
00045   std::string viscType = visc_list->get("Type", "Constant");
00046 
00047   if (viscType == "Constant"){ 
00048     std::cout << "Constant viscosity!" << std::endl;
00049     visc_type = CONSTANT;
00050   }
00051   else if (viscType == "Sutherland") {
00052    std::cout << "Sutherland viscosity!" << std::endl; 
00053     visc_type = SUTHERLAND; 
00054   }
00055   
00056   muref = visc_list->get("Mu_ref", 1.0); 
00057   kapparef = visc_list->get("Kappa_ref", 1.0);  
00058   Tref = visc_list->get("T_ref", 1.0);  
00059   Pr = visc_list->get("Prandtl number Pr", 0.72); 
00060   Cp = visc_list->get("Specific heat Cp", 1.0); 
00061   
00062   coordVec = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00063             p.get<std::string>("Coordinate Vector Name"),
00064       p.get<Teuchos::RCP<PHX::DataLayout> >("QP Gradient Data Layout") );
00065 
00066   this->addDependentField(qFluct);
00067   this->addDependentField(qFluctGrad);
00068   this->addDependentField(coordVec);
00069   this->addEvaluatedField(mu);
00070   this->addEvaluatedField(kappa);
00071   this->addEvaluatedField(lambda);
00072   this->addEvaluatedField(tau11);
00073   this->addEvaluatedField(tau12);
00074   this->addEvaluatedField(tau13);
00075   this->addEvaluatedField(tau22);
00076   this->addEvaluatedField(tau23);
00077   this->addEvaluatedField(tau33);
00078 
00079   std::vector<PHX::DataLayout::size_type> dims;
00080   qFluctGrad.fieldTag().dataLayout().dimensions(dims);
00081   numQPs   = dims[2];
00082   numDims  = dims[3];
00083 
00084   qFluct.fieldTag().dataLayout().dimensions(dims);
00085   vecDim  = dims[2];
00086 
00087   std::cout << "vecdim in viscosity evaluator: " << vecDim << std::endl; 
00088   std::cout << "numDims in viscosity evaluator: " << numDims << std::endl; 
00089   std::cout << "numQPs in viscosity evaluator: " << numQPs << std::endl; 
00090   std::cout << "Mu_ref: " << muref << std::endl; 
00091   std::cout << "Kappa_ref: " << kapparef << std::endl; 
00092 
00093   this->setName("ComprNSViscosity"+PHX::TypeString<EvalT>::value);
00094 }
00095 
00096 //**********************************************************************
00097 template<typename EvalT, typename Traits>
00098 void ComprNSViscosity<EvalT, Traits>::
00099 postRegistrationSetup(typename Traits::SetupData d,
00100                       PHX::FieldManager<Traits>& fm)
00101 {
00102   this->utils.setFieldData(qFluct,fm);
00103   this->utils.setFieldData(qFluctGrad,fm);
00104   this->utils.setFieldData(mu,fm); 
00105   this->utils.setFieldData(kappa,fm); 
00106   this->utils.setFieldData(lambda,fm); 
00107   this->utils.setFieldData(coordVec,fm); 
00108   this->utils.setFieldData(tau11,fm); 
00109   this->utils.setFieldData(tau12,fm); 
00110   this->utils.setFieldData(tau13,fm); 
00111   this->utils.setFieldData(tau22,fm); 
00112   this->utils.setFieldData(tau23,fm); 
00113   this->utils.setFieldData(tau33,fm); 
00114 }
00115 
00116 //**********************************************************************
00117 template<typename EvalT, typename Traits>
00118 void ComprNSViscosity<EvalT, Traits>::
00119 evaluateFields(typename Traits::EvalData workset)
00120 {
00121   //Visocisity coefficients
00122   if (visc_type == CONSTANT){
00123     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00124       for (std::size_t qp=0; qp < numQPs; ++qp) {
00125         mu(cell,qp) = 1.0;
00126         kappa(cell,qp) = mu(cell,qp)*Cp/Pr/kapparef; 
00127         mu(cell,qp) = 1.0/muref; //non-dimensionalize mu 
00128         lambda(cell,qp) = -2.0/3.0*mu(cell,qp); //Stokes' hypothesis 
00129       }
00130     }
00131   }
00132   else if (visc_type == SUTHERLAND){
00133     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00134       for (std::size_t qp=0; qp < numQPs; ++qp) {
00135         ScalarT T = qFluct(cell,qp,vecDim-1)*Tref; //temperature (dimensional)
00136         mu(cell,qp) = (1.458e-6)*sqrt(T*T*T)/(T + 110.4); //mu = (1.458e-6)*T^(1/5)/(T + 110.4) 
00137         kappa(cell,qp) = mu(cell,qp)*Cp/Pr/kapparef; 
00138         mu(cell,qp) = mu(cell,qp)/muref; //non-dimensionalize mu 
00139         lambda(cell,qp) = -2.0/3.0*mu(cell,qp); //Stokes' hypothesis
00140       }
00141     }
00142   }
00143   //Viscous stresses
00144   for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00145     for (std::size_t qp=0; qp < numQPs; ++qp) {
00146       tau11(cell,qp) = mu(cell,qp)*2.0*qFluctGrad(cell,qp,1,0) + lambda(cell,qp)*(qFluctGrad(cell,qp,1,0) + qFluctGrad(cell,qp,2,1)); //mu*2*du/dx + lambda*div(u) 
00147       tau12(cell,qp) = mu(cell,qp)*(qFluctGrad(cell,qp,1,1) + qFluctGrad(cell,qp,2,0)); //mu*(du/dy + dv/dx)
00148       tau13(cell,qp) = 0.0; 
00149       tau22(cell,qp) = mu(cell,qp)*2.0*qFluctGrad(cell,qp,2,1) + lambda(cell,qp)*(qFluctGrad(cell,qp,1,0) + qFluctGrad(cell,qp,2,1)); //mu*2*dv/dy + lambda*div(u) 
00150       tau23(cell,qp) = 0.0; 
00151       tau33(cell,qp) = 0.0; 
00152     }
00153   }
00154   if (numDims == 3) {//3D case 
00155     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00156       for (std::size_t qp=0; qp < numQPs; ++qp) {
00157         tau11(cell,qp) += lambda(cell,qp)*qFluctGrad(cell,qp,3,2); //+lambda*dw/dz 
00158         tau13(cell,qp) += mu(cell,qp)*(qFluctGrad(cell,qp,1,2) + qFluctGrad(cell,qp,3,0)); //mu*(du/dz + dw/dx)
00159         tau22(cell,qp) += lambda(cell,qp)*qFluctGrad(cell,qp,3,2); //+lambda*dw/dz 
00160         tau23(cell,qp) += mu(cell,qp)*(qFluctGrad(cell,qp,2,3) + qFluctGrad(cell,qp,3,1)); //mu*(dv/dz + dw/dy)
00161         tau33(cell,qp) += mu(cell,qp)*2.0*qFluctGrad(cell,qp,3,2) + lambda(cell,qp)*(qFluctGrad(cell,qp,1,0) + qFluctGrad(cell,qp,2,1) + qFluct(cell,qp,3,2)); //mu*2*dw/dz + lambda*div(u) 
00162       }
00163     }
00164   }
00165 }
00166 
00167 }
00168
PHAL_ComprNSViscosity_Def.hpp
