Albany: a Trilinos-based PDE code: FELIX

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 #ifndef FELIX_STOKES_HPP
00008 #define FELIX_STOKES_HPP
00009 
00010 #include "Teuchos_RCP.hpp"
00011 #include "Teuchos_ParameterList.hpp"
00012 
00013 #include "Albany_AbstractProblem.hpp"
00014 
00015 #include "Phalanx.hpp"
00016 #include "PHAL_Workset.hpp"
00017 #include "PHAL_Dimension.hpp"
00018 
00019 namespace FELIX {
00020 
00025   class Stokes : public Albany::AbstractProblem {
00026   public:
00027   
00029     Stokes(const Teuchos::RCP<Teuchos::ParameterList>& params,
00030      const Teuchos::RCP<ParamLib>& paramLib,
00031      const int numDim_);
00032 
00034     ~Stokes();
00035 
00037     virtual int spatialDimension() const { return numDim; }
00038 
00040     virtual void buildProblem(
00041       Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> >  meshSpecs,
00042       Albany::StateManager& stateMgr);
00043 
00044     // Build evaluators
00045     virtual Teuchos::Array< Teuchos::RCP<const PHX::FieldTag> >
00046     buildEvaluators(
00047       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00048       const Albany::MeshSpecsStruct& meshSpecs,
00049       Albany::StateManager& stateMgr,
00050       Albany::FieldManagerChoice fmchoice,
00051       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00052 
00054     Teuchos::RCP<const Teuchos::ParameterList> getValidProblemParameters() const;
00055 
00056   private:
00057 
00059     Stokes(const Stokes&);
00060     
00062     Stokes& operator=(const Stokes&);
00063 
00064   public:
00065 
00067     template <typename EvalT> Teuchos::RCP<const PHX::FieldTag>
00068     constructEvaluators(
00069       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00070       const Albany::MeshSpecsStruct& meshSpecs,
00071       Albany::StateManager& stateMgr,
00072       Albany::FieldManagerChoice fmchoice,
00073       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00074 
00075     void constructDirichletEvaluators(const Albany::MeshSpecsStruct& meshSpecs);
00076     void constructNeumannEvaluators(const Teuchos::RCP<Albany::MeshSpecsStruct>& meshSpecs);
00077 
00078 
00079   protected:
00080 
00082     enum NS_VAR_TYPE {
00083       NS_VAR_TYPE_NONE,      
00084       NS_VAR_TYPE_CONSTANT,  
00085       NS_VAR_TYPE_DOF        
00086     };
00087 
00088     void getVariableType(Teuchos::ParameterList& paramList,
00089        const std::string& defaultType,
00090        NS_VAR_TYPE& variableType,
00091        bool& haveVariable,
00092        bool& haveEquation);
00093     std::string variableTypeToString(const NS_VAR_TYPE variableType);
00094 
00095   protected:
00096     
00097     int numDim;        
00098 
00099     NS_VAR_TYPE flowType; 
00100 
00101     bool haveFlow;     
00102 
00103     bool haveFlowEq;     
00104 
00105     bool haveSource;   
00106     bool havePSPG;     
00107 
00108    Teuchos::RCP<Albany::Layouts> dl; 
00109     
00110   };
00111 
00112 }
00113 
00114 #include "Intrepid_FieldContainer.hpp"
00115 #include "Intrepid_DefaultCubatureFactory.hpp"
00116 #include "Shards_CellTopology.hpp"
00117 
00118 #include "Albany_Utils.hpp"
00119 #include "Albany_ProblemUtils.hpp"
00120 #include "Albany_EvaluatorUtils.hpp"
00121 #include "Albany_ResponseUtilities.hpp"
00122 
00123 #include "FELIX_StokesContravarientMetricTensor.hpp"
00124 #include "PHAL_Neumann.hpp"
00125 #include "FELIX_StokesBodyForce.hpp"
00126 #include "FELIX_StokesRm.hpp"
00127 #include "FELIX_StokesTauM.hpp"
00128 #include "FELIX_StokesMomentumResid.hpp"
00129 #include "FELIX_StokesContinuityResid.hpp"
00130 #include "FELIX_Viscosity.hpp"
00131 
00132 template <typename EvalT>
00133 Teuchos::RCP<const PHX::FieldTag>
00134 FELIX::Stokes::constructEvaluators(
00135   PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00136   const Albany::MeshSpecsStruct& meshSpecs,
00137   Albany::StateManager& stateMgr,
00138   Albany::FieldManagerChoice fieldManagerChoice,
00139   const Teuchos::RCP<Teuchos::ParameterList>& responseList)
00140 {
00141   using Teuchos::RCP;
00142   using Teuchos::rcp;
00143   using Teuchos::ParameterList;
00144   using PHX::DataLayout;
00145   using PHX::MDALayout;
00146   using std::vector;
00147   using std::string;
00148   using std::map;
00149   using PHAL::AlbanyTraits;
00150   
00151   RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > >
00152     intrepidBasis = Albany::getIntrepidBasis(meshSpecs.ctd);
00153   RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology (&meshSpecs.ctd));
00154   
00155   const int numNodes = intrepidBasis->getCardinality();
00156   const int worksetSize = meshSpecs.worksetSize;
00157   
00158   Intrepid::DefaultCubatureFactory<RealType> cubFactory;
00159   RCP <Intrepid::Cubature<RealType> > cubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree);
00160   
00161   const int numQPts = cubature->getNumPoints();
00162   const int numVertices = cellType->getNodeCount();
00163   
00164   *out << "Field Dimensions: Workset=" << worksetSize 
00165        << ", Vertices= " << numVertices
00166        << ", Nodes= " << numNodes
00167        << ", QuadPts= " << numQPts
00168        << ", Dim= " << numDim << std::endl;
00169   
00170 
00171    dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPts,numDim, numDim));
00172    TEUCHOS_TEST_FOR_EXCEPTION(dl->vectorAndGradientLayoutsAreEquivalent==false, std::logic_error,
00173                               "Data Layout Usage in Stokes problem assumes vecDim = numDim");
00174    Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
00175    bool supportsTransient=true;
00176    int offset=0;
00177 
00178    // Temporary variable used numerous times below
00179    Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
00180 
00181    // Define Field Names
00182 
00183    if (haveFlowEq) {
00184      Teuchos::ArrayRCP<std::string> dof_names(1);
00185      Teuchos::ArrayRCP<std::string> dof_names_dot(1);
00186      Teuchos::ArrayRCP<std::string> resid_names(1);
00187      dof_names[0] = "Velocity";
00188      dof_names_dot[0] = dof_names[0]+"_dot";
00189      resid_names[0] = "Momentum Residual";
00190      fm0.template registerEvaluator<EvalT>
00191        (evalUtils.constructGatherSolutionEvaluator(true, dof_names, dof_names_dot, offset));
00192 
00193      fm0.template registerEvaluator<EvalT>
00194        (evalUtils.constructDOFVecInterpolationEvaluator(dof_names[0], offset));
00195 
00196      fm0.template registerEvaluator<EvalT>
00197        (evalUtils.constructDOFVecInterpolationEvaluator(dof_names_dot[0], offset));
00198 
00199      fm0.template registerEvaluator<EvalT>
00200        (evalUtils.constructDOFVecGradInterpolationEvaluator(dof_names[0], offset));
00201 
00202      fm0.template registerEvaluator<EvalT>
00203        (evalUtils.constructScatterResidualEvaluator(true, resid_names,offset, "Scatter Momentum"));
00204      offset += numDim;
00205    }
00206 
00207    if (haveFlowEq) {
00208      Teuchos::ArrayRCP<std::string> dof_names(1);
00209      Teuchos::ArrayRCP<std::string> dof_names_dot(1);
00210      Teuchos::ArrayRCP<std::string> resid_names(1);
00211      dof_names[0] = "Pressure";
00212      dof_names_dot[0] = dof_names[0]+"_dot";
00213      resid_names[0] = "Continuity Residual";
00214      fm0.template registerEvaluator<EvalT>
00215        (evalUtils.constructGatherSolutionEvaluator(false, dof_names, dof_names_dot, offset));
00216 
00217      fm0.template registerEvaluator<EvalT>
00218        (evalUtils.constructDOFInterpolationEvaluator(dof_names[0], offset));
00219 
00220      fm0.template registerEvaluator<EvalT>
00221        (evalUtils.constructDOFInterpolationEvaluator(dof_names_dot[0], offset));
00222 
00223      fm0.template registerEvaluator<EvalT>
00224        (evalUtils.constructDOFGradInterpolationEvaluator(dof_names[0], offset));
00225 
00226      fm0.template registerEvaluator<EvalT>
00227        (evalUtils.constructScatterResidualEvaluator(false, resid_names,offset, "Scatter Continuity"));
00228      offset ++;
00229    }
00230 
00231 
00232    fm0.template registerEvaluator<EvalT>
00233      (evalUtils.constructGatherCoordinateVectorEvaluator());
00234 
00235    fm0.template registerEvaluator<EvalT>
00236      (evalUtils.constructMapToPhysicalFrameEvaluator(cellType, cubature));
00237 
00238    fm0.template registerEvaluator<EvalT>
00239      (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cubature));
00240 
00241   if (havePSPG) { // Compute Contravarient Metric Tensor
00242     RCP<ParameterList> p = 
00243       rcp(new ParameterList("Contravarient Metric Tensor"));
00244 
00245     // Inputs: X, Y at nodes, Cubature, and Basis
00246     p->set<std::string>("Coordinate Vector Name","Coord Vec");
00247     p->set< RCP<Intrepid::Cubature<RealType> > >("Cubature", cubature);
00248 
00249     p->set<RCP<shards::CellTopology> >("Cell Type", cellType);
00250 
00251     // Outputs: BF, weightBF, Grad BF, weighted-Grad BF, all in physical space
00252     p->set<std::string>("Contravarient Metric Tensor Name", "Gc");
00253 
00254     ev = rcp(new FELIX::StokesContravarientMetricTensor<EvalT,AlbanyTraits>(*p,dl));
00255     fm0.template registerEvaluator<EvalT>(ev);
00256   }
00257 
00258   
00259 
00260   if (haveFlowEq) { // Body Force
00261     RCP<ParameterList> p = rcp(new ParameterList("Body Force"));
00262 
00263     //Input
00264     p->set<std::string>("FELIX Viscosity QP Variable Name", "FELIX Viscosity");
00265     p->set<std::string>("Coordinate Vector Name", "Coord Vec");
00266 
00267     Teuchos::ParameterList& paramList = params->sublist("Body Force");
00268     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00269       
00270 
00271     //Output
00272     p->set<std::string>("Body Force Name", "Body Force");
00273 
00274     ev = rcp(new FELIX::StokesBodyForce<EvalT,AlbanyTraits>(*p,dl));
00275     fm0.template registerEvaluator<EvalT>(ev);
00276   }
00277 
00278 
00279 
00280   if (haveFlowEq) { // Rm
00281     RCP<ParameterList> p = rcp(new ParameterList("Rm"));
00282 
00283     //Input
00284     p->set<std::string>("Velocity QP Variable Name", "Velocity");
00285     p->set<std::string>("Velocity Dot QP Variable Name", "Velocity_dot");
00286     p->set<std::string>("Velocity Gradient QP Variable Name", "Velocity Gradient");
00287     p->set<std::string>("Pressure Gradient QP Variable Name", "Pressure Gradient");
00288     p->set<std::string>("Body Force QP Variable Name", "Body Force");
00289     p->set<std::string>("Coordinate Vector Name", "Coord Vec");
00290 
00291     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00292   
00293     //Output
00294     p->set<std::string>("Rm Name", "Rm");
00295 
00296     ev = rcp(new FELIX::StokesRm<EvalT,AlbanyTraits>(*p,dl));
00297     fm0.template registerEvaluator<EvalT>(ev);
00298   }
00299   
00300   //IK, 7/24/2012
00301   if (haveFlowEq) { // FELIX viscosity
00302     RCP<ParameterList> p = rcp(new ParameterList("FELIX Viscosity"));
00303 
00304     //Input
00305     p->set<std::string>("Velocity Gradient QP Variable Name", "Velocity Gradient");
00306     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00307     Teuchos::ParameterList& paramList = params->sublist("FELIX Viscosity");
00308     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00309     p->set<std::string>("Coordinate Vector Name", "Coord Vec");
00310   
00311     //Output
00312     p->set<std::string>("FELIX Viscosity QP Variable Name", "FELIX Viscosity");
00313 
00314     ev = rcp(new FELIX::Viscosity<EvalT,AlbanyTraits>(*p,dl));
00315     fm0.template registerEvaluator<EvalT>(ev);
00316   }
00317 
00318 
00319   if (haveFlowEq && havePSPG) { // Tau M
00320     RCP<ParameterList> p = rcp(new ParameterList("Tau M"));
00321 
00322     //Input
00323     p->set<std::string>("Velocity QP Variable Name", "Velocity");
00324     p->set<std::string>("Contravarient Metric Tensor Name", "Gc"); 
00325     p->set<std::string>("Jacobian Det Name", "Jacobian Det"); 
00326     p->set<std::string>("FELIX Viscosity QP Variable Name", "FELIX Viscosity");
00327 
00328     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00329     Teuchos::ParameterList& paramList = params->sublist("Tau M");
00330     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00331 
00332     //Output
00333     p->set<std::string>("Tau M Name", "Tau M");
00334 
00335     ev = rcp(new FELIX::StokesTauM<EvalT,AlbanyTraits>(*p,dl));
00336     fm0.template registerEvaluator<EvalT>(ev);
00337   }
00338 
00339   if (haveFlowEq) { // Momentum Resid
00340     RCP<ParameterList> p = rcp(new ParameterList("Momentum Resid"));
00341 
00342     //Input
00343     p->set<std::string>("Weighted BF Name", "wBF");
00344     p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00345     p->set<std::string>("Velocity QP Variable Name", "Velocity");
00346     p->set<std::string>("Velocity Gradient QP Variable Name", "Velocity Gradient");
00347     p->set<std::string>("Pressure QP Variable Name", "Pressure");
00348     p->set<std::string>("Pressure Gradient QP Variable Name", "Pressure Gradient");
00349     p->set<std::string>("FELIX Viscosity QP Variable Name", "FELIX Viscosity");
00350     p->set<std::string>("Body Force Name", "Body Force");
00351 
00352     p->set<std::string>("Velocity QP Variable Name", "Velocity");
00353     p->set<std::string>("Density QP Variable Name", "Density");
00354     p->set<std::string> ("Tau M Name", "Tau M");
00355  
00356     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00357   
00358     //Output
00359     p->set<std::string>("Residual Name", "Momentum Residual");
00360 
00361     ev = rcp(new FELIX::StokesMomentumResid<EvalT,AlbanyTraits>(*p,dl));
00362     fm0.template registerEvaluator<EvalT>(ev);
00363   }
00364   
00365 
00366   if (haveFlowEq) { // Continuity Resid
00367     RCP<ParameterList> p = rcp(new ParameterList("Continuity Resid"));
00368 
00369     //Input
00370     p->set<std::string>("Weighted BF Name", "wBF");
00371     p->set<std::string>("Velocity QP Variable Name", "Velocity");
00372     p->set<std::string>("Gradient QP Variable Name", "Velocity Gradient");
00373     p->set<std::string>("Density QP Variable Name", "Density");
00374 
00375     p->set<bool>("Have PSPG", havePSPG);
00376     p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00377     p->set<std::string> ("Tau M Name", "Tau M");
00378     p->set<std::string> ("Rm Name", "Rm");
00379 
00380     //Output
00381     p->set<std::string>("Residual Name", "Continuity Residual");
00382 
00383     ev = rcp(new FELIX::StokesContinuityResid<EvalT,AlbanyTraits>(*p,dl));
00384     fm0.template registerEvaluator<EvalT>(ev);
00385   }
00386 
00387 
00388   if (fieldManagerChoice == Albany::BUILD_RESID_FM)  {
00389     Teuchos::RCP<const PHX::FieldTag> ret_tag;
00390     if (haveFlowEq) {
00391       PHX::Tag<typename EvalT::ScalarT> mom_tag("Scatter Momentum", dl->dummy);
00392       fm0.requireField<EvalT>(mom_tag);
00393       PHX::Tag<typename EvalT::ScalarT> con_tag("Scatter Continuity", dl->dummy);
00394       fm0.requireField<EvalT>(con_tag);
00395       ret_tag = mom_tag.clone();
00396     }
00397     return ret_tag;
00398   }
00399   else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
00400     Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
00401     return respUtils.constructResponses(fm0, *responseList, stateMgr);
00402   }
00403 
00404   return Teuchos::null;
00405 }
00406 #endif // FELIX_STOKES_HPP
FELIX_Stokes.hpp
