Albany: a Trilinos-based PDE code: FELIX_StokesFOBodyForce

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 "Teuchos_VerboseObject.hpp"
00009 #include "Phalanx_DataLayout.hpp"
00010 #include "Sacado.hpp"
00011 
00012 #include "Intrepid_FunctionSpaceTools.hpp"
00013 
00014 //uncomment the following line if you want debug output to be printed to screen
00015 //#define OUTPUT_TO_SCREEN
00016 
00017 namespace FELIX {
00018 const double pi = 3.1415926535897932385;
00019 const double g = 9.8; //gravity for FELIX; hard-coded here for now
00020 const double rho = 910; //density for FELIX; hard-coded here for now
00021 const double rho_g = rho*g; //density for FELIX; hard-coded here for now
00022 
00023 //**********************************************************************
00024 
00025 template<typename EvalT, typename Traits>
00026 StokesFOBodyForce<EvalT, Traits>::
00027 StokesFOBodyForce(const Teuchos::ParameterList& p,
00028                   const Teuchos::RCP<Albany::Layouts>& dl) :
00029   force(p.get<std::string>("Body Force Name"), dl->qp_vector), 
00030   A(1.0), 
00031   n(3.0), 
00032   alpha(0.0)
00033 {
00034   Teuchos::RCP<Teuchos::FancyOStream> out(Teuchos::VerboseObjectBase::getDefaultOStream());
00035   Teuchos::ParameterList* bf_list = 
00036     p.get<Teuchos::ParameterList*>("Parameter List");
00037 
00038   std::string type = bf_list->get("Type", "None");
00039   A = bf_list->get("Glen's Law A", 1.0); 
00040   n = bf_list->get("Glen's Law n", 3.0); 
00041   alpha = bf_list->get("FELIX alpha", 0.0);
00042 #ifdef OUTPUT_TO_SCREEN
00043   *out << "alpha: " << alpha << std::endl;
00044 #endif 
00045   alpha *= pi/180.0; //convert alpha to radians 
00046   if (type == "None") {
00047     bf_type = NONE;
00048   }
00049   else if (type == "FO INTERP SURF GRAD") {
00050 #ifdef OUTPUT_TO_SCREEN
00051     *out << "INTERP SURFACE GRAD Source!" << std::endl;
00052 #endif
00053     surfaceGrad = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00054              p.get<std::string>("Surface Height Gradient Name"), dl->qp_gradient);
00055     this->addDependentField(surfaceGrad);
00056      bf_type = FO_INTERP_SURF_GRAD;
00057   }
00058   else if (type == "FOSinCos2D") {
00059     bf_type = FO_SINCOS2D;  
00060     muFELIX = PHX::MDField<ScalarT,Cell,QuadPoint>(
00061             p.get<std::string>("FELIX Viscosity QP Variable Name"), dl->qp_scalar);
00062     coordVec = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00063             p.get<std::string>("Coordinate Vector Name"), dl->qp_gradient);
00064     this->addDependentField(muFELIX); 
00065     this->addDependentField(coordVec);
00066   }
00067   else if (type == "FOSinExp2D") {
00068     bf_type = FO_SINEXP2D;  
00069     muFELIX = PHX::MDField<ScalarT,Cell,QuadPoint>(
00070             p.get<std::string>("FELIX Viscosity QP Variable Name"), dl->qp_scalar);
00071     coordVec = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00072             p.get<std::string>("Coordinate Vector Name"), dl->qp_gradient);
00073     this->addDependentField(muFELIX); 
00074     this->addDependentField(coordVec);
00075   }
00076   else if (type == "FOCosExp2D") {
00077     bf_type = FO_COSEXP2D;  
00078     muFELIX = PHX::MDField<ScalarT,Cell,QuadPoint>(
00079             p.get<std::string>("FELIX Viscosity QP Variable Name"), dl->qp_scalar);
00080     coordVec = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00081             p.get<std::string>("Coordinate Vector Name"), dl->qp_gradient);
00082     this->addDependentField(muFELIX); 
00083     this->addDependentField(coordVec);
00084   }
00085   else if (type == "FOCosExp2DFlip") {
00086     bf_type = FO_COSEXP2DFLIP;  
00087     muFELIX = PHX::MDField<ScalarT,Cell,QuadPoint>(
00088             p.get<std::string>("FELIX Viscosity QP Variable Name"), dl->qp_scalar);
00089     coordVec = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00090             p.get<std::string>("Coordinate Vector Name"), dl->qp_gradient);
00091     this->addDependentField(muFELIX); 
00092     this->addDependentField(coordVec);
00093   }
00094   else if (type == "FOCosExp2DAll") {
00095     bf_type = FO_COSEXP2DALL;  
00096     muFELIX = PHX::MDField<ScalarT,Cell,QuadPoint>(
00097             p.get<std::string>("FELIX Viscosity QP Variable Name"), dl->qp_scalar);
00098     coordVec = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00099             p.get<std::string>("Coordinate Vector Name"), dl->qp_gradient);
00100     this->addDependentField(muFELIX); 
00101     this->addDependentField(coordVec);
00102   }
00103   else if (type == "FOSinCosZ") {
00104     bf_type = FO_SINCOSZ;  
00105     muFELIX = PHX::MDField<ScalarT,Cell,QuadPoint>(
00106             p.get<std::string>("FELIX Viscosity QP Variable Name"), dl->qp_scalar);
00107     coordVec = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00108             p.get<std::string>("Coordinate Vector Name"), dl->qp_gradient);
00109     this->addDependentField(muFELIX); 
00110     this->addDependentField(coordVec);
00111   }
00112   else if (type == "Poisson") {
00113     bf_type = POISSON;  
00114     muFELIX = PHX::MDField<ScalarT,Cell,QuadPoint>(
00115             p.get<std::string>("FELIX Viscosity QP Variable Name"), dl->qp_scalar);
00116     coordVec = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00117             p.get<std::string>("Coordinate Vector Name"), dl->qp_gradient);
00118     this->addDependentField(muFELIX); 
00119     this->addDependentField(coordVec);
00120   }
00121   //kept for backward compatibility. Use type = "FO INTERP GRAD SURF" instead.
00122   else if ((type == "FO ISMIP-HOM Test A") || (type == "FO ISMIP-HOM Test B") || (type == "FO ISMIP-HOM Test C") || (type == "FO ISMIP-HOM Test D")) {
00123   *out << "ISMIP-HOM Tests A/B/C/D \n WARNING: computing INTERP SURFACE GRAD Source! \nPlease set  Force Type = FO INTERP GRAD SURF." << std::endl;
00124     surfaceGrad = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00125         p.get<std::string>("Surface Height Gradient Name"), dl->qp_gradient);
00126     this->addDependentField(surfaceGrad);
00127     bf_type = FO_INTERP_SURF_GRAD;
00128   }
00129   else if (type == "FO Dome") {
00130 #ifdef OUTPUT_TO_SCREEN
00131     *out << "Dome Source!" << std::endl;
00132 #endif 
00133     coordVec = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00134             p.get<std::string>("Coordinate Vector Name"), dl->qp_gradient);
00135     bf_type = FO_DOME; 
00136     this->addDependentField(coordVec);
00137   }
00138 
00139   this->addEvaluatedField(force);
00140 
00141   std::vector<PHX::DataLayout::size_type> dims;
00142   dl->qp_gradient->dimensions(dims);
00143   numQPs  = dims[1];
00144   numDims = dims[2];
00145   dl->qp_vector->dimensions(dims);
00146   vecDim  = dims[2];
00147 
00148 //*out << " in FELIX Stokes FO source! " << std::endl;
00149 //*out << " vecDim = " << vecDim << std::endl;
00150 //*out << " numDims = " << numDims << std::endl;
00151 //*out << " numQPs = " << numQPs << std::endl; 
00152 
00153   this->setName("StokesFOBodyForce"+PHX::TypeString<EvalT>::value);
00154 }
00155 
00156 //**********************************************************************
00157 template<typename EvalT, typename Traits>
00158 void StokesFOBodyForce<EvalT, Traits>::
00159 postRegistrationSetup(typename Traits::SetupData d,
00160                       PHX::FieldManager<Traits>& fm)
00161 {
00162   if (bf_type == FO_SINCOS2D || bf_type == FO_SINEXP2D || bf_type == FO_COSEXP2D || bf_type == FO_COSEXP2DFLIP || 
00163       bf_type == FO_COSEXP2DALL || bf_type == FO_SINCOSZ || bf_type == POISSON) {
00164     this->utils.setFieldData(muFELIX,fm);
00165     this->utils.setFieldData(coordVec,fm);
00166   }
00167   else if (bf_type == FO_DOME) {
00168     this->utils.setFieldData(coordVec,fm);
00169   }
00170   else if (bf_type == FO_INTERP_SURF_GRAD)
00171     this->utils.setFieldData(surfaceGrad,fm);
00172 
00173   this->utils.setFieldData(force,fm); 
00174 }
00175 
00176 //**********************************************************************
00177 template<typename EvalT, typename Traits>
00178 void StokesFOBodyForce<EvalT, Traits>::
00179 evaluateFields(typename Traits::EvalData workset)
00180 {
00181  if (bf_type == NONE) {
00182    for (std::size_t cell=0; cell < workset.numCells; ++cell) 
00183      for (std::size_t qp=0; qp < numQPs; ++qp)       
00184        for (std::size_t i=0; i < vecDim; ++i) 
00185      force(cell,qp,i) = 0.0;
00186  }
00187  //source using the gradient of the interpolated surface height
00188  else if (bf_type == FO_INTERP_SURF_GRAD) {
00189    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00190      for (std::size_t qp=0; qp < numQPs; ++qp) {
00191        force(cell,qp,0) = rho_g*surfaceGrad(cell,qp,0);
00192        force(cell,qp,1) = rho_g*surfaceGrad(cell,qp,1);
00193      }
00194    }
00195  }
00196  else if (bf_type == FO_SINCOS2D) {
00197    double xphase=0.0, yphase=0.0;
00198    double r = 3*pi;
00199    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00200      for (std::size_t qp=0; qp < numQPs; ++qp) {      
00201        ScalarT* f = &force(cell,qp,0);
00202        MeshScalarT x2pi = 2.0*pi*coordVec(cell,qp,0);
00203        MeshScalarT y2pi = 2.0*pi*coordVec(cell,qp,1);
00204        MeshScalarT muargt = 2.0*pi*cos(x2pi + xphase)*cos(y2pi + yphase) + r;  
00205        MeshScalarT muqp = 0.5*pow(A, -1.0/n)*pow(muargt, 1.0/n - 1.0);
00206        MeshScalarT dmuargtdx = -4.0*pi*pi*sin(x2pi + xphase)*cos(y2pi + yphase); 
00207        MeshScalarT dmuargtdy = -4.0*pi*pi*cos(x2pi + xphase)*sin(y2pi + yphase); 
00208        MeshScalarT exx = 2.0*pi*cos(x2pi + xphase)*cos(y2pi + yphase) + r; 
00209        MeshScalarT eyy = -2.0*pi*cos(x2pi + xphase)*cos(y2pi + yphase) - r; 
00210        MeshScalarT exy = 0.0;  
00211        f[0] = 2.0*muqp*(-4.0*pi*pi*sin(x2pi + xphase)*cos(y2pi + yphase))  
00212             + 2.0*0.5*pow(A, -1.0/n)*(1.0/n - 1.0)*pow(muargt, 1.0/n - 2.0)*(dmuargtdx*(2.0*exx + eyy) + dmuargtdy*exy); 
00213        f[1] = 2.0*muqp*(4.0*pi*pi*cos(x2pi + xphase)*sin(y2pi + yphase)) 
00214             + 2.0*0.5*pow(A, -1.0/n)*(1.0/n - 1.0)*pow(muargt, 1.0/n - 2.0)*(dmuargtdx*exy + dmuargtdy*(exx + 2.0*eyy));
00215      }
00216    }
00217  }
00218  else if (bf_type == FO_SINEXP2D) {
00219    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00220      for (std::size_t qp=0; qp < numQPs; ++qp) {      
00221        ScalarT* f = &force(cell,qp,0);
00222        MeshScalarT x2pi = 2.0*pi*coordVec(cell,qp,0);
00223        MeshScalarT y2pi = 2.0*pi*coordVec(cell,qp,1);
00224        MeshScalarT muqp = 1.0 ; //0.5*pow(A, -1.0/n)*pow(muargt, 1.0/n - 1.0);
00225        f[0] = -1.0*(-4.0*muqp*exp(coordVec(cell,qp,0)) - 12.0*muqp*pi*pi*cos(x2pi)*cos(y2pi) + 4.0*muqp*pi*pi*cos(y2pi));   
00226        f[1] =  -1.0*(20.0*muqp*pi*pi*sin(x2pi)*sin(y2pi)); 
00227      }
00228    }
00229  }
00230  else if (bf_type == POISSON) { //source term for debugging of Neumann BC
00231    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00232      for (std::size_t qp=0; qp < numQPs; ++qp) {      
00233        ScalarT* f = &force(cell,qp,0);
00234        MeshScalarT x = coordVec(cell,qp,0);
00235        f[0] = exp(x);  
00236      }
00237    }
00238  }
00239  else if (bf_type == FO_COSEXP2D) {
00240    const double a = 1.0; 
00241    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00242      for (std::size_t qp=0; qp < numQPs; ++qp) {      
00243        ScalarT* f = &force(cell,qp,0);
00244        MeshScalarT x2pi = 2.0*pi*coordVec(cell,qp,0);
00245        MeshScalarT x = coordVec(cell,qp,0);
00246        MeshScalarT y2pi = 2.0*pi*coordVec(cell,qp,1);
00247        MeshScalarT muqp = 1.0 ; //0.5*pow(A, -1.0/n)*pow(muargt, 1.0/n - 1.0);
00248        f[0] = 2.0*muqp*(2.0*a*a*exp(a*x)*cos(y2pi) + 6.0*pi*pi*cos(x2pi)*cos(y2pi) - 2.0*pi*pi*exp(a*x)*cos(y2pi)); 
00249        f[1] = 2.0*muqp*(-3.0*pi*a*exp(a*x)*sin(y2pi) - 10.0*pi*pi*sin(x2pi)*sin(y2pi)); 
00250      }
00251    }
00252  }
00253  else if (bf_type == FO_COSEXP2DFLIP) {
00254    const double a = 1.0; 
00255    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00256      for (std::size_t qp=0; qp < numQPs; ++qp) {      
00257        ScalarT* f = &force(cell,qp,0);
00258        MeshScalarT x2pi = 2.0*pi*coordVec(cell,qp,0);
00259        MeshScalarT x = coordVec(cell,qp,0);
00260        MeshScalarT y2pi = 2.0*pi*coordVec(cell,qp,1);
00261        MeshScalarT muqp = 1.0 ; //0.5*pow(A, -1.0/n)*pow(muargt, 1.0/n - 1.0);
00262        f[0] = 2.0*muqp*(-3.0*pi*a*exp(a*x)*sin(y2pi) - 10.0*pi*pi*sin(x2pi)*sin(y2pi)); 
00263        f[1] = 2.0*muqp*(1.0/2.0*a*a*exp(a*x)*cos(y2pi) + 6.0*pi*pi*cos(x2pi)*cos(y2pi) - 8.0*pi*pi*exp(a*x)*cos(y2pi)); 
00264      }
00265    }
00266  }
00267  else if (bf_type == FO_COSEXP2DALL) {
00268    const double a = 1.0; 
00269    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00270      for (std::size_t qp=0; qp < numQPs; ++qp) {
00271        ScalarT* f = &force(cell,qp,0);
00272        MeshScalarT x2pi = 2.0*pi*coordVec(cell,qp,0);
00273        MeshScalarT x = coordVec(cell,qp,0);
00274        MeshScalarT y2pi = 2.0*pi*coordVec(cell,qp,1);
00275        MeshScalarT muargt = (a*a + 4.0*pi*pi - 2.0*pi*a)*sin(y2pi)*sin(y2pi) + 1.0/4.0*(2.0*pi+a)*(2.0*pi+a)*cos(y2pi)*cos(y2pi);
00276        muargt = sqrt(muargt)*exp(a*x);
00277        MeshScalarT muqp = 1.0/2.0*pow(A, -1.0/n)*pow(muargt, 1.0/n - 1.0);
00278        MeshScalarT dmuargtdx = a*muargt;
00279        MeshScalarT dmuargtdy = 3.0/2.0*pi*(a*a+4.0*pi*pi-4.0*pi*a)*cos(y2pi)*sin(y2pi)*exp(a*x)/sqrt((a*a + 4.0*pi*pi - 2.0*pi*a)*sin(y2pi)*sin(y2pi) + 1.0/4.0*(2.0*pi+a)*(2.0*pi+a)*cos(y2pi)*cos(y2pi));
00280        MeshScalarT exx = a*exp(a*x)*sin(y2pi);
00281        MeshScalarT eyy = -2.0*pi*exp(a*x)*sin(y2pi);
00282        MeshScalarT exy = 1.0/2.0*(2.0*pi+a)*exp(a*x)*cos(y2pi);
00283        f[0] = 2.0*muqp*(2.0*a*a*exp(a*x)*sin(y2pi) - 3.0*pi*a*exp(a*x)*sin(y2pi) - 2.0*pi*pi*exp(a*x)*sin(y2pi))
00284             + 2.0*0.5*pow(A, -1.0/n)*(1.0/n-1.0)*pow(muargt, 1.0/n-2.0)*(dmuargtdx*(2.0*exx + eyy) + dmuargtdy*exy);
00285        f[1] = 2.0*muqp*(3.0*a*pi*exp(a*x)*cos(y2pi) + 1.0/2.0*a*a*exp(a*x)*cos(y2pi) - 8.0*pi*pi*exp(a*x)*cos(y2pi))
00286             + 2.0*0.5*pow(A, -1.0/n)*(1.0/n-1.0)*pow(muargt, 1.0/n-2.0)*(dmuargtdx*exy + dmuargtdy*(exx + 2.0*eyy));
00287      }
00288    }
00289  }
00290  // Doubly-periodic MMS with polynomial in Z for FO Stokes
00291  else if (bf_type == FO_SINCOSZ) {
00292    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00293      for (std::size_t qp=0; qp < numQPs; ++qp) {      
00294        ScalarT* f = &force(cell,qp,0);
00295        MeshScalarT x2pi = 2.0*pi*coordVec(cell,qp,0);
00296        MeshScalarT y2pi = 2.0*pi*coordVec(cell,qp,1);
00297        MeshScalarT& z = coordVec(cell,qp,2);
00298        MeshScalarT muqp = 1.0; //hard coded to constant for now 
00299        
00300        ScalarT t1 = z*(1.0-z)*(1.0-2.0*z); 
00301        ScalarT t2 = 2.0*z - 1.0; 
00302 
00303        f[0] = 2.0*muqp*(-16.0*pi*pi*t1 + 3.0*t2)*sin(x2pi)*sin(y2pi); 
00304        f[1] = 2.0*muqp*(16.0*pi*pi*t1 - 3.0*t2)*cos(x2pi)*cos(y2pi);  
00305      }
00306    }
00307  }
00308  //source for dome test case 
00309  else if (bf_type == FO_DOME) {
00310    for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00311      for (std::size_t qp=0; qp < numQPs; ++qp) {      
00312        ScalarT* f = &force(cell,qp,0);
00313        MeshScalarT x = coordVec(cell,qp,0);
00314        MeshScalarT y = coordVec(cell,qp,1);
00315        f[0] = -rho_g*x*0.7071/sqrt(450.0-x*x-y*y)/sqrt(450.0);  
00316        f[1] = -rho_g*y*0.7071/sqrt(450.0-x*x-y*y)/sqrt(450.0);  
00317      }
00318    }
00319  }
00320 }
00321 
00322 }
00323
FELIX_StokesFOBodyForce_Def.hpp
