Albany: a Trilinos-based PDE code: FELIX_StokesBodyForce

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 FELIX {
00014 const double pi = 3.1415926535897932385;
00015 const double g = 9.8; //gravity for FELIX; hard-coded here for now
00016 const double rho = 910; //density for FELIX; hard-coded here for now
00017 const double L = 5; 
00018 const double alpha = 0.5*pi/180; 
00019 const double Z = 1; 
00020 const double U = L*pow(2*rho*g*Z,3); 
00021 const double cx = 1.0; //1e-9*U; 
00022 const double cy = 1.0; //1e-9*U; 
00023 
00024 //**********************************************************************
00025 //function returning u solution to Poly source test case
00026 template<typename ScalarT>
00027 ScalarT u2d(const ScalarT& x, const ScalarT& y){
00028   ScalarT u = 20.0*(x-1.0)*(x-1.0)*x*x*y*(2.0*y*y-3.0*y+1.0);
00029   return u; 
00030 }
00031 
00032 //function returning v solution to Poly source test case
00033 template<typename ScalarT>
00034 ScalarT v2d(const ScalarT& x, const ScalarT& y){
00035   ScalarT v = -20.0*x*(2.0*x*x-3.0*x+1.0)*y*y*(1.0-y)*(1.0-y); 
00036   return v; 
00037 }
00038 
00039 //function returning p solution to Poly source test case
00040 template<typename ScalarT>
00041 ScalarT p2d(const ScalarT& x, const ScalarT& y){
00042   ScalarT p = (5.0/2.0*y*y-10.0*x) + 4.1613;
00043   return p; 
00044 }
00045 
00046 //function returning tau11 stress for Poly source term test case
00047 template<typename ScalarT>
00048 ScalarT tau112d(const ScalarT& x, const ScalarT& y){
00049   int num_deriv = 1;  
00050   Sacado::Fad::DFad<ScalarT> xfad(num_deriv, 0, x);
00051   Sacado::Fad::DFad<ScalarT> yfad(y); 
00052   Sacado::Fad::DFad<ScalarT> ufad = u2d(xfad, yfad); 
00053   ScalarT dudx_ad = ufad.dx(0); 
00054   ScalarT p = p2d(x,y); 
00055   ScalarT tau = 2.0*dudx_ad - p;   
00056   return tau; 
00057 }
00058 
00059 //function returning tau12 stress for Poly source term test case
00060 template<typename ScalarT>
00061 ScalarT tau122d(const ScalarT& x, const ScalarT& y){
00062   int num_deriv = 2;  
00063   Sacado::Fad::DFad<ScalarT> xfad(num_deriv, 0, x);
00064   Sacado::Fad::DFad<ScalarT> yfad(num_deriv, 1, y); 
00065   Sacado::Fad::DFad<ScalarT> ufad = u2d(xfad, yfad); 
00066   Sacado::Fad::DFad<ScalarT> vfad = v2d(xfad, yfad); 
00067   ScalarT dudy_ad = ufad.dx(1); 
00068   ScalarT dvdx_ad = vfad.dx(0); 
00069   ScalarT tau = (dudy_ad + dvdx_ad);   
00070   return tau; 
00071 }
00072 
00073 //function returning tau22 stress for Poly source term test case
00074 template<typename ScalarT>
00075 ScalarT tau222d(const ScalarT& x, const ScalarT& y){
00076   int num_deriv = 1;  
00077   Sacado::Fad::DFad<ScalarT> xfad(x);
00078   Sacado::Fad::DFad<ScalarT> yfad(num_deriv, 0, y); 
00079   Sacado::Fad::DFad<ScalarT> vfad = v2d(xfad, yfad); 
00080   ScalarT dvdy_ad = vfad.dx(0); 
00081   ScalarT p = p2d(x,y);  
00082   ScalarT tau = 2.0*dvdy_ad - p;   
00083   return tau; 
00084 }
00085 
00086 //returns function giving top surface boundary for MMF Test A ISMIP-HOM test case
00087 template<typename ScalarT>
00088 ScalarT top_surf(const ScalarT& x, const ScalarT& y){
00089   ScalarT s = -x*tan(alpha); 
00090   return s; 
00091 }
00092 
00093 //returns function basal surface boundary for MMF TestA ISMIP-HOM test case
00094 template<typename ScalarT>
00095 ScalarT basal_surf(const ScalarT& x, const ScalarT& y){
00096   ScalarT s = top_surf(x, y);
00097   ScalarT b = s + Z/2.0*sin(2.0*pi*x/L)*sin(2.0*pi*y/L) - Z; 
00098   return b; 
00099 }
00100 
00101 //returns expression for u-velocity for MMF Test A ISMIP-HOM test case 
00102 template<typename ScalarT>
00103 ScalarT u_vel(const ScalarT& x, const ScalarT& y, const ScalarT& z){
00104   ScalarT b = basal_surf(x, y); 
00105   ScalarT s = top_surf(x, y); 
00106   ScalarT a = (s-z)/(s-b); 
00107   ScalarT u = cx*(1.0-a*a*a*a); 
00108   return u; 
00109 }
00110 
00111 
00112 //returns expression for v-velocity for MMF Test A ISMIP-HOM test case 
00113 template<typename ScalarT>
00114 ScalarT v_vel(const ScalarT& x, const ScalarT& y, const ScalarT& z){
00115   ScalarT b = basal_surf(x, y); 
00116   ScalarT s = top_surf(x, y); 
00117   ScalarT a = (s-z)/(s-b); 
00118   ScalarT a4 = a*a*a*a; 
00119   ScalarT v = cy/(s-b)*(1.0-a4) - 0.5*cx/(s-b)*(1.0-a4)*Z*cos(2.0*pi*x/L)*cos(2*pi*y/L); //exp(ct*t) term??
00120   return v; 
00121 }
00122 
00123 //returns expression for w-velocity for MMF Test A ISMIP-HOM test case 
00124 template<typename ScalarT>
00125 ScalarT w_vel(const ScalarT& x, const ScalarT& y, const ScalarT& z){
00126   int num_deriv = 2;  
00127   Sacado::Fad::DFad<ScalarT> xfad(num_deriv, 0, x);
00128   Sacado::Fad::DFad<ScalarT> yfad(num_deriv, 1, y); 
00129   Sacado::Fad::DFad<ScalarT> bfad = basal_surf(xfad, yfad); 
00130   Sacado::Fad::DFad<ScalarT> sfad = top_surf(xfad, yfad); 
00131   ScalarT dbdx_ad = bfad.dx(0); 
00132   ScalarT dbdy_ad = bfad.dx(1); 
00133   ScalarT dsdx_ad = sfad.dx(0); 
00134   ScalarT dsdy_ad = sfad.dx(1); 
00135   ScalarT u = u_vel(x,y,z); 
00136   ScalarT v = v_vel(x,y,z); 
00137   ScalarT s = top_surf(x,y); 
00138   ScalarT b = basal_surf(x,y);
00139   ScalarT w = u*(dbdx_ad*(s-z)/(s-b) + dsdx_ad*(z-b)/(s-b)) + 
00140               v*(dbdy_ad*(s-z)/(s-b) + dsdy_ad*(z-b)/(s-b)); 
00141   return w; 
00142 }
00143 
00144 //returns expression for mu for MMF Test A ISMIP-HOM test case 
00145 template<typename ScalarT>
00146 ScalarT mu_eval(const ScalarT& x, const ScalarT& y, const ScalarT& z){
00147   /*int num_deriv = 3;  
00148   Sacado::Fad::DFad<ScalarT> xfad(num_deriv, 0, x);
00149   Sacado::Fad::DFad<ScalarT> yfad(num_deriv, 1, y); 
00150   Sacado::Fad::DFad<ScalarT> zfad(num_deriv, 2, y);
00151   Sacado::Fad::DFad<ScalarT> ufad = u_vel(xfad, yfad, zfad); 
00152   Sacado::Fad::DFad<ScalarT> vfad = v_vel(xfad, yfad, zfad); 
00153   Sacado::Fad::DFad<ScalarT> wfad = w_vel(xfad, yfad, zfad); 
00154   ScalarT dudx_ad = ufad.dx(0); 
00155   ScalarT dudy_ad = ufad.dx(1); 
00156   ScalarT dudz_ad = ufad.dx(2); 
00157   ScalarT dvdx_ad = vfad.dx(0); 
00158   ScalarT dvdy_ad = vfad.dx(1); 
00159   ScalarT dvdz_ad = vfad.dx(2); 
00160   ScalarT dwdx_ad = wfad.dx(0); 
00161   ScalarT dwdy_ad = wfad.dx(1); 
00162   ScalarT dwdz_ad = wfad.dx(2); 
00163   ScalarT mu = 0.25*(dudy_ad + dvdx_ad)*(dudy_ad + dvdx_ad) + 
00164                0.25*(dudz_ad + dwdx_ad)*(dudz_ad + dwdx_ad) + 
00165                0.25*(dvdz_ad + dwdy_ad)*(dvdz_ad + dwdy_ad) - 
00166                dudx_ad*dvdy_ad - dudx_ad*dwdz_ad - dvdy_ad*dwdz_ad; 
00167   mu = pow(mu, (1.0-n)/(2.0*n)); 
00168   mu *= 0.25*pow(A, -1.0/n);
00169   */ 
00170   ScalarT mu = 1.0;  
00171   return mu; 
00172 }
00173 
00174 //returns expression for pressure for MMF Test A ISMIP-HOM test case 
00175 template<typename ScalarT>
00176 ScalarT press(const ScalarT& x, const ScalarT& y, const ScalarT& z){
00177   int num_deriv = 2;  
00178   Sacado::Fad::DFad<ScalarT> xfad(num_deriv, 0, x);
00179   Sacado::Fad::DFad<ScalarT> yfad(num_deriv, 1, y); 
00180   Sacado::Fad::DFad<ScalarT> zfad(z); 
00181   Sacado::Fad::DFad<ScalarT> ufad = u_vel(xfad, yfad, zfad); 
00182   Sacado::Fad::DFad<ScalarT> vfad = v_vel(xfad, yfad, zfad); 
00183   ScalarT dudx_ad = ufad.dx(0); 
00184   ScalarT dvdy_ad = vfad.dx(1);
00185   ScalarT mu = mu_eval(x, y, z); 
00186   ScalarT s = top_surf(x,y); 
00187   ScalarT p = 2.0*mu*dudx_ad + 2.0*mu*dvdy_ad - rho*g*(s-z);  
00188   return p; 
00189 }
00190 
00191 //returns expression for tau11 for MMF Test A ISMIP-HOM test case 
00192 template<typename ScalarT>
00193 ScalarT tau11(const ScalarT& x, const ScalarT& y, const ScalarT& z){
00194   int num_deriv = 1;  
00195   Sacado::Fad::DFad<ScalarT> xfad(num_deriv, 0, x);
00196   Sacado::Fad::DFad<ScalarT> yfad(y); 
00197   Sacado::Fad::DFad<ScalarT> zfad(z); 
00198   Sacado::Fad::DFad<ScalarT> ufad = u_vel(xfad, yfad, zfad); 
00199   ScalarT dudx_ad = ufad.dx(0); 
00200   ScalarT mu = mu_eval(x, y, z); 
00201   ScalarT p = press(x,y,z); 
00202   ScalarT tau = 2.0*mu*dudx_ad - p;   
00203   return tau; 
00204 }
00205 
00206 //returns expression for tau12 for MMF Test A ISMIP-HOM test case 
00207 template<typename ScalarT>
00208 ScalarT tau12(const ScalarT& x, const ScalarT& y, const ScalarT& z){
00209   int num_deriv = 2;  
00210   Sacado::Fad::DFad<ScalarT> xfad(num_deriv, 0, x);
00211   Sacado::Fad::DFad<ScalarT> yfad(num_deriv, 1, y); 
00212   Sacado::Fad::DFad<ScalarT> zfad(z); 
00213   Sacado::Fad::DFad<ScalarT> ufad = u_vel(xfad, yfad, zfad); 
00214   Sacado::Fad::DFad<ScalarT> vfad = v_vel(xfad, yfad, zfad); 
00215   ScalarT dudy_ad = ufad.dx(1); 
00216   ScalarT dvdx_ad = vfad.dx(0); 
00217   ScalarT mu = mu_eval(x, y, z); 
00218   ScalarT tau = mu*(dudy_ad + dvdx_ad);   
00219   return tau; 
00220 }
00221 
00222 //returns expression for tau13 for MMF Test A ISMIP-HOM test case 
00223 template<typename ScalarT>
00224 ScalarT tau13(const ScalarT& x, const ScalarT& y, const ScalarT& z){
00225   int num_deriv = 2;  
00226   Sacado::Fad::DFad<ScalarT> xfad(num_deriv, 0, x);
00227   Sacado::Fad::DFad<ScalarT> zfad(num_deriv, 1, z); 
00228   Sacado::Fad::DFad<ScalarT> yfad(y); 
00229   Sacado::Fad::DFad<ScalarT> ufad = u_vel(xfad, yfad, zfad); 
00230   Sacado::Fad::DFad<ScalarT> wfad = w_vel(xfad, yfad, zfad); 
00231   ScalarT dudz_ad = ufad.dx(1); 
00232   ScalarT dwdx_ad = wfad.dx(0); 
00233   ScalarT mu = mu_eval(x, y, z); 
00234   ScalarT tau = mu*(dudz_ad + dwdx_ad);   
00235   return tau; 
00236 }
00237 
00238 //returns expression for tau22 for MMF Test A ISMIP-HOM test case 
00239 template<typename ScalarT>
00240 ScalarT tau22(const ScalarT& x, const ScalarT& y, const ScalarT& z){
00241   int num_deriv = 1;  
00242   Sacado::Fad::DFad<ScalarT> xfad(x);
00243   Sacado::Fad::DFad<ScalarT> yfad(num_deriv, 0, y); 
00244   Sacado::Fad::DFad<ScalarT> zfad(z); 
00245   Sacado::Fad::DFad<ScalarT> vfad = v_vel(xfad, yfad, zfad); 
00246   ScalarT dvdy_ad = vfad.dx(0); 
00247   ScalarT mu = mu_eval(x, y, z);
00248   ScalarT p = press(x,y,z);  
00249   ScalarT tau = mu*(2.0*mu*dvdy_ad - p);   
00250   return tau; 
00251 }
00252 
00253 //returns expression for tau23 for MMF Test A ISMIP-HOM test case 
00254 template<typename ScalarT>
00255 ScalarT tau23(const ScalarT& x, const ScalarT& y, const ScalarT& z){
00256   int num_deriv = 2;  
00257   Sacado::Fad::DFad<ScalarT> xfad(x);
00258   Sacado::Fad::DFad<ScalarT> yfad(num_deriv, 0, y); 
00259   Sacado::Fad::DFad<ScalarT> zfad(num_deriv, 1, z); 
00260   Sacado::Fad::DFad<ScalarT> vfad = v_vel(xfad, yfad, zfad); 
00261   Sacado::Fad::DFad<ScalarT> wfad = w_vel(xfad, yfad, zfad); 
00262   ScalarT dvdz_ad = vfad.dx(1); 
00263   ScalarT dwdy_ad = wfad.dx(0); 
00264   ScalarT mu = mu_eval(x, y, z);
00265   ScalarT tau = mu*(dvdz_ad + dwdy_ad);   
00266   return tau; 
00267 }
00268 
00269 //returns expression for tau33 for MMF Test A ISMIP-HOM test case 
00270 template<typename ScalarT>
00271 ScalarT tau33(const ScalarT& x, const ScalarT& y, const ScalarT& z){
00272   int num_deriv = 1;  
00273   Sacado::Fad::DFad<ScalarT> xfad(x);
00274   Sacado::Fad::DFad<ScalarT> yfad(y); 
00275   Sacado::Fad::DFad<ScalarT> zfad(num_deriv, 0, z); 
00276   Sacado::Fad::DFad<ScalarT> wfad = w_vel(xfad, yfad, zfad); 
00277   ScalarT dwdz_ad = wfad.dx(0); 
00278   ScalarT mu = mu_eval(x, y, z);
00279   ScalarT p = press(x,y,z);  
00280   ScalarT tau = 2.0*mu*dwdz_ad - p;   
00281   return tau; 
00282 }
00283 
00284 
00285 template<typename EvalT, typename Traits>
00286 StokesBodyForce<EvalT, Traits>::
00287 StokesBodyForce(const Teuchos::ParameterList& p,
00288                 const Teuchos::RCP<Albany::Layouts>& dl) :
00289   force(p.get<std::string>("Body Force Name"),dl->qp_vector),
00290   A(1.0), 
00291   n(3.0)
00292 {
00293 
00294   Teuchos::ParameterList* bf_list = 
00295     p.get<Teuchos::ParameterList*>("Parameter List");
00296 
00297   std::string type = bf_list->get("Type", "None");
00298   A = bf_list->get("Glen's Law A", 1.0); 
00299   n = bf_list->get("Glen's Law n", 3.0); 
00300   if (type == "None") {
00301     bf_type = NONE;
00302   }
00303   else if (type == "Gravity") {
00304     bf_type = GRAVITY;
00305   }
00306   else if (type == "Poly") {
00307     bf_type = POLY;  
00308     muFELIX = PHX::MDField<ScalarT,Cell,QuadPoint>(
00309             p.get<std::string>("FELIX Viscosity QP Variable Name"),dl->qp_scalar);
00310     coordVec = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00311             p.get<std::string>("Coordinate Vector Name"),dl->qp_gradient);
00312     this->addDependentField(muFELIX); 
00313     this->addDependentField(coordVec);
00314   }
00315   else if (type == "PolySacado") {
00316     bf_type = POLYSACADO;  
00317     muFELIX = PHX::MDField<ScalarT,Cell,QuadPoint>(
00318             p.get<std::string>("FELIX Viscosity QP Variable Name"), dl->qp_scalar);
00319     coordVec = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00320             p.get<std::string>("Coordinate Vector Name"), dl->qp_gradient);
00321     this->addDependentField(muFELIX); 
00322     this->addDependentField(coordVec);
00323   }
00324   else if (type == "SinSin") {
00325     bf_type = SINSIN;  
00326     muFELIX = PHX::MDField<ScalarT,Cell,QuadPoint>(
00327             p.get<std::string>("FELIX Viscosity QP Variable Name"),dl->qp_scalar);
00328     coordVec = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00329             p.get<std::string>("Coordinate Vector Name"), dl->qp_gradient);
00330     this->addDependentField(muFELIX); 
00331     this->addDependentField(coordVec);
00332   }
00333   else if (type == "FullStokesBasal") {
00334     bf_type = FULLSTOKESBASAL;  
00335     muFELIX = PHX::MDField<ScalarT,Cell,QuadPoint>(
00336             p.get<std::string>("FELIX Viscosity QP Variable Name"),dl->qp_scalar);
00337     coordVec = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00338             p.get<std::string>("Coordinate Vector Name"), dl->qp_gradient);
00339     this->addDependentField(muFELIX); 
00340     this->addDependentField(coordVec);
00341   }
00342   else if (type == "SinSinGlen") {
00343     bf_type = SINSINGLEN;  
00344     muFELIX = PHX::MDField<ScalarT,Cell,QuadPoint>(
00345             p.get<std::string>("FELIX Viscosity QP Variable Name"),dl->qp_scalar);
00346     coordVec = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00347             p.get<std::string>("Coordinate Vector Name"), dl->qp_gradient);
00348     this->addDependentField(muFELIX); 
00349     this->addDependentField(coordVec);
00350   }
00351   else if (type == "SinCosZ") {
00352     bf_type = SINCOSZ;  
00353     muFELIX = PHX::MDField<ScalarT,Cell,QuadPoint>(
00354             p.get<std::string>("FELIX Viscosity QP Variable Name"),dl->qp_scalar);
00355     coordVec = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00356             p.get<std::string>("Coordinate Vector Name"), dl->qp_gradient);
00357     this->addDependentField(muFELIX); 
00358     this->addDependentField(coordVec);
00359   }
00360   else if (type == "TestAMMF") {
00361     bf_type = TESTAMMF;  
00362     muFELIX = PHX::MDField<ScalarT,Cell,QuadPoint>(
00363             p.get<std::string>("FELIX Viscosity QP Variable Name"),dl->qp_scalar);
00364     coordVec = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00365             p.get<std::string>("Coordinate Vector Name"), dl->qp_gradient);
00366     this->addDependentField(muFELIX); 
00367     this->addDependentField(coordVec);
00368   }
00369 
00370   this->addEvaluatedField(force);
00371 
00372   std::vector<PHX::DataLayout::size_type> dims;
00373   dl->qp_gradient->dimensions(dims);
00374   numQPs  = dims[1];
00375   numDims = dims[2];
00376 
00377   if (bf_type == GRAVITY) {
00378     if (bf_list->isType<Teuchos::Array<double> >("Gravity Vector"))
00379       gravity = bf_list->get<Teuchos::Array<double> >("Gravity Vector");
00380     else {
00381       gravity.resize(numDims);
00382       gravity[numDims-1] = -1.0;
00383     }
00384   }
00385   this->setName("StokesBodyForce"+PHX::TypeString<EvalT>::value);
00386 }
00387 
00388 //**********************************************************************
00389 template<typename EvalT, typename Traits>
00390 void StokesBodyForce<EvalT, Traits>::
00391 postRegistrationSetup(typename Traits::SetupData d,
00392                       PHX::FieldManager<Traits>& fm)
00393 {
00394   if (bf_type == GRAVITY) {
00395   }
00396   else if (bf_type == POLY || bf_type == SINSIN || bf_type == SINSINGLEN || bf_type == FULLSTOKESBASAL || bf_type == SINCOSZ || bf_type == POLYSACADO || bf_type == TESTAMMF) {
00397     this->utils.setFieldData(muFELIX,fm);
00398     this->utils.setFieldData(coordVec,fm);
00399   }
00400 
00401   this->utils.setFieldData(force,fm); 
00402 }
00403 
00404 //**********************************************************************
00405 template<typename EvalT, typename Traits>
00406 void StokesBodyForce<EvalT, Traits>::
00407 evaluateFields(typename Traits::EvalData workset)
00408 {
00409  if (bf_type == NONE) {
00410    for (std::size_t cell=0; cell < workset.numCells; ++cell) 
00411      for (std::size_t qp=0; qp < numQPs; ++qp)       
00412        for (std::size_t i=0; i < numDims; ++i) 
00413      force(cell,qp,i) = 0.0;
00414  }
00415  else if (bf_type == GRAVITY) {
00416  for (std::size_t cell=0; cell < workset.numCells; ++cell) 
00417    for (std::size_t qp=0; qp < numQPs; ++qp) {
00418      for (std::size_t i=0; i < numDims; ++i) { 
00419    force(cell,qp,i) = -1.0*g*rho*gravity[i];
00420       }
00421    }
00422  }
00423 
00424  //The following is hard-coded for a 2D Stokes problem with manufactured solution
00425  else if (bf_type == POLY) {
00426    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00427      for (std::size_t qp=0; qp < numQPs; ++qp) {      
00428        ScalarT* f = &force(cell,qp,0);
00429        MeshScalarT* X = &coordVec(cell,qp,0);
00430        ScalarT& muqp = muFELIX(cell,qp);
00431        f[0] =  40.0*muqp*(2.0*X[1]*X[1] - 3.0*X[1]+1.0)*X[1]*(6.0*X[0]*X[0] -6.0*X[0] + 1.0)
00432               + 120*muqp*(X[0]-1.0)*(X[0]-1.0)*X[0]*X[0]*(2.0*X[1]-1.0) 
00433               + 10.0*muqp;      
00434        f[1] = - 120.0*muqp*(1.0-X[1])*(1.0-X[1])*X[1]*X[1]*(2.0*X[0]-1.0)
00435               - 40.0*muqp*(2.0*X[0]*X[0] - 3.0*X[0] + 1.0)*X[0]*(6.0*X[1]*X[1] - 6.0*X[1] + 1.0)
00436               - 5*muqp*X[1];
00437      }
00438    }
00439  }
00440  // Doubly-periodic MMS derived by Irina. 
00441  else if (bf_type == SINSIN) {
00442    double xphase=0.0, yphase=0.0; // Expose as parameters for verification
00443    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00444      for (std::size_t qp=0; qp < numQPs; ++qp) {      
00445        ScalarT* f = &force(cell,qp,0);
00446        MeshScalarT x2pi = 2.0*pi*coordVec(cell,qp,0);
00447        MeshScalarT y2pi = 2.0*pi*coordVec(cell,qp,1);
00448        ScalarT& muqp = muFELIX(cell,qp);
00449 
00450        f[0] = -4.0*muqp*pi*(2*pi-1)*sin(x2pi + xphase)*sin(y2pi + yphase);
00451        f[1] = -4.0*muqp*pi*(2*pi+1)*cos(x2pi + xphase)*cos(y2pi + yphase);
00452      }
00453    }
00454  }
00455  else if (bf_type == FULLSTOKESBASAL) {
00456    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00457      for (std::size_t qp=0; qp < numQPs; ++qp) {      
00458        ScalarT* f = &force(cell,qp,0);
00459        MeshScalarT x = coordVec(cell,qp,0);
00460        MeshScalarT x2pi = 2.0*pi*coordVec(cell,qp,0);
00461        MeshScalarT y2pi = 2.0*pi*coordVec(cell,qp,1);
00462        ScalarT& muqp = muFELIX(cell,qp);
00463        f[0] = -2.0*pi*sin(y2pi)*(-1.0*muqp*exp(x) + cos(x2pi) + 4.0*muqp*pi*pi*exp(x));
00464        f[1] = -1.0*cos(y2pi)*(4.0*muqp*pi*pi*exp(x) - muqp*exp(x)  + 2.0*pi*sin(x2pi)); 
00465      }
00466    }
00467  }
00468  //MMS test case for 2D nonlinear Stokes with Glen's law on unit square 
00469  else if (bf_type == SINSINGLEN) {
00470    double xphase=0.0, yphase=0.0; 
00471    double r = 3*pi;
00472    //cout << "muqp force: " << endl; 
00473    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00474      for (std::size_t qp=0; qp < numQPs; ++qp) {      
00475        ScalarT* f = &force(cell,qp,0);
00476        MeshScalarT x2pi = 2.0*pi*coordVec(cell,qp,0);
00477        MeshScalarT y2pi = 2.0*pi*coordVec(cell,qp,1);
00478        //MeshScalarT muqp = x2pi*y2pi + r;
00479        MeshScalarT muargt = 2.0*pi*cos(x2pi + xphase)*cos(y2pi + yphase) + r;  
00480        MeshScalarT muqp = 0.5*pow(A, -1.0/n)*pow(muargt, 1.0/n - 1.0);
00481        MeshScalarT dudx = 2.0*pi*cos(x2pi + xphase)*cos(y2pi + yphase) + r; 
00482        MeshScalarT dudy = -2.0*pi*sin(x2pi + xphase)*cos(y2pi + yphase); 
00483        MeshScalarT dvdx = 2.0*pi*sin(x2pi + xphase)*sin(y2pi+yphase); 
00484        MeshScalarT dvdy = -2.0*pi*cos(x2pi + xphase)*cos(y2pi + yphase) - r;
00485        MeshScalarT dmuargtdx = -4.0*pi*pi*sin(x2pi + xphase)*cos(y2pi + yphase); 
00486        MeshScalarT dmuargtdy = -4.0*pi*pi*cos(x2pi + xphase)*sin(y2pi + yphase);  
00487        //cout.precision(10);
00488        //if (cell == 10)  
00489        //   cout << muqp << endl; 
00490        f[0] = -8.0*pi*pi*sin(x2pi + xphase)*cos(y2pi + yphase)*(muqp - 1.0) 
00491             + 0.5*pow(A, -1.0/n)*(1.0/n - 1.0)*pow(muargt, 1.0/n - 2.0)*(dmuargtdx*dudx + dmuargtdy*dudy);
00492        f[1] = 8.0*pi*pi*cos(x2pi + xphase)*sin(y2pi + yphase)*(muqp + 1.0) 
00493             + 0.5*pow(A, -1.0/n)*(1.0/n - 1.0)*pow(muargt, 1.0/n - 2.0)*(dmuargtdx*dvdx + dmuargtdy*dvdy);
00494      }
00495    }
00496  }
00497  else if (bf_type == POLYSACADO) {
00498    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00499      for (std::size_t qp=0; qp < numQPs; ++qp) {      
00500        ScalarT* f = &force(cell,qp,0);
00501        MeshScalarT& x = coordVec(cell,qp,0);
00502        MeshScalarT& y = coordVec(cell,qp,1);
00503        ScalarT X = x; ScalarT Y = y; 
00504        int num_deriv = 2; 
00505        Sacado::Fad::DFad<ScalarT> xfad(num_deriv, 0, X); 
00506        Sacado::Fad::DFad<ScalarT> yfad(num_deriv, 1, Y); 
00507        Sacado::Fad::DFad<ScalarT> tau11fad = tau112d(xfad, yfad); 
00508        Sacado::Fad::DFad<ScalarT> tau12fad = tau122d(xfad, yfad); 
00509        Sacado::Fad::DFad<ScalarT> tau22fad = tau222d(xfad, yfad); 
00510        f[0] = 1.0*(tau11fad.dx(0) + tau12fad.dx(1)); 
00511        f[1] = 1.0*(tau12fad.dx(0) + tau22fad.dx(1));
00512      }
00513    }
00514  }
00515  else if (bf_type == TESTAMMF) {
00516    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00517      for (std::size_t qp=0; qp < numQPs; ++qp) {      
00518        ScalarT* f = &force(cell,qp,0);
00519        MeshScalarT& x = coordVec(cell,qp,0);
00520        MeshScalarT& y = coordVec(cell,qp,1);
00521        MeshScalarT& z = coordVec(cell,qp,2);
00522        ScalarT X = x; ScalarT Y = y; ScalarT Z = z; 
00523        int num_deriv = 3; 
00524        Sacado::Fad::DFad<ScalarT> xfad(num_deriv, 0, X); 
00525        Sacado::Fad::DFad<ScalarT> yfad(num_deriv, 1, Y); 
00526        Sacado::Fad::DFad<ScalarT> zfad(num_deriv, 2, Z);
00527        Sacado::Fad::DFad<ScalarT> tau11fad = tau11(xfad, yfad, zfad); 
00528        Sacado::Fad::DFad<ScalarT> tau12fad = tau12(xfad, yfad, zfad); 
00529        Sacado::Fad::DFad<ScalarT> tau13fad = tau13(xfad, yfad, zfad); 
00530        Sacado::Fad::DFad<ScalarT> tau22fad = tau22(xfad, yfad, zfad); 
00531        Sacado::Fad::DFad<ScalarT> tau23fad = tau23(xfad, yfad, zfad); 
00532        Sacado::Fad::DFad<ScalarT> tau33fad = tau33(xfad, yfad, zfad); 
00533        f[0] = tau11fad.dx(0) + tau12fad.dx(1) + tau13fad.dx(2);  
00534        f[1] = tau12fad.dx(0) + tau22fad.dx(1) + tau23fad.dx(2); 
00535        f[2] = tau13fad.dx(0) + tau23fad.dx(1) + tau33fad.dx(2) + rho*g; 
00536      }
00537    }
00538  }
00539  // Doubly-periodic MMS with polynomial in Z
00540  else if (bf_type == SINCOSZ) {
00541    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00542      for (std::size_t qp=0; qp < numQPs; ++qp) {      
00543        ScalarT* f = &force(cell,qp,0);
00544        MeshScalarT x2pi = 2.0*pi*coordVec(cell,qp,0);
00545        MeshScalarT y2pi = 2.0*pi*coordVec(cell,qp,1);
00546        MeshScalarT& z = coordVec(cell,qp,2);
00547        ScalarT& muqp = muFELIX(cell,qp);
00548 
00549        /*ScalarT t1 = muqp*z*(1-z)*(1-2*z);
00550        ScalarT t2 = muqp*(1-6*z+6*z*z);
00551        ScalarT t3 = muqp*z*z*(1-z)*(1-z);
00552 
00553        f[0] =  (-8*pi*pi*t1 + t2 + 4*pi*pi)*sin(x2pi)*sin(y2pi);
00554        f[1] = -(-8*pi*pi*t1 + t2 + 4*pi*pi)*cos(x2pi)*cos(y2pi);
00555        f[2] = -2*pi*(-8*pi*pi*t3 + 2*t1 + 1)*cos(x2pi)*sin(y2pi);
00556        */
00557        ScalarT t1 = muqp*z*(1.0-z)*(1.0-2.0*z); 
00558        ScalarT t2 = muqp*(12.0*z - 6.0); 
00559        ScalarT t3 = muqp*z*z*(1.0-z)*(1.0-z); 
00560        ScalarT t4 = muqp*(1.0-6.0*z+6.0*z*z); 
00561      
00562        f[0] = (-8.0*pi*pi*t1 + t2 + 4.0*pi*pi*z)*sin(x2pi)*sin(y2pi); 
00563        f[1] = (8.0*pi*pi*t1 - t2 - 4.0*pi*pi*z)*cos(x2pi)*cos(y2pi); 
00564        f[2] = (16.0*pi*pi*pi*t3 - 4.0*pi*t4 - 2.0*pi)*cos(x2pi)*sin(y2pi); 
00565      }
00566    }
00567  }
00568 }
00569 
00570 }
00571
FELIX_StokesBodyForce_Def.hpp
