Albany: a Trilinos-based PDE code: Albany_ComprNSProblem.hpp Source File

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 ALBANY_COMPRNSPROBLEM_HPP
00008 #define ALBANY_COMPRNSPROBLEM_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 Albany {
00020 
00025   class ComprNSProblem : public AbstractProblem {
00026   public:
00027   
00029     ComprNSProblem(const Teuchos::RCP<Teuchos::ParameterList>& params,
00030      const Teuchos::RCP<ParamLib>& paramLib,
00031      const int numDim_);
00032 
00034     ~ComprNSProblem();
00035 
00037     virtual int spatialDimension() const { return numDim; }
00038 
00040     virtual void buildProblem(
00041       Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> >  meshSpecs,
00042       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     ComprNSProblem(const ComprNSProblem&);
00060     
00062     ComprNSProblem& operator=(const ComprNSProblem&);
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   protected:
00079     int numDim;
00080     Teuchos::RCP<Albany::Layouts> dl;
00081 
00082 
00083   };
00084 
00085 }
00086 
00087 #include "Intrepid_FieldContainer.hpp"
00088 #include "Intrepid_DefaultCubatureFactory.hpp"
00089 #include "Shards_CellTopology.hpp"
00090 
00091 #include "Albany_Utils.hpp"
00092 #include "Albany_ProblemUtils.hpp"
00093 #include "Albany_EvaluatorUtils.hpp"
00094 #include "Albany_ResponseUtilities.hpp"
00095 
00096 #include "PHAL_DOFVecGradInterpolation.hpp"
00097 
00098 #include "PHAL_ComprNSResid.hpp"
00099 
00100 #include "PHAL_Neumann.hpp"
00101 #include "PHAL_Source.hpp"
00102 #include "PHAL_ComprNSBodyForce.hpp"
00103 #include "PHAL_ComprNSViscosity.hpp"
00104 
00105 template <typename EvalT>
00106 Teuchos::RCP<const PHX::FieldTag>
00107 Albany::ComprNSProblem::constructEvaluators(
00108   PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00109   const Albany::MeshSpecsStruct& meshSpecs,
00110   Albany::StateManager& stateMgr,
00111   Albany::FieldManagerChoice fieldManagerChoice,
00112   const Teuchos::RCP<Teuchos::ParameterList>& responseList)
00113 {
00114   using Teuchos::RCP;
00115   using Teuchos::rcp;
00116   using Teuchos::ParameterList;
00117   using PHX::DataLayout;
00118   using PHX::MDALayout;
00119   using std::vector;
00120   using std::string;
00121   using std::map;
00122   using PHAL::AlbanyTraits;
00123   
00124   RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > >
00125     intrepidBasis = Albany::getIntrepidBasis(meshSpecs.ctd);
00126   RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology (&meshSpecs.ctd));
00127   
00128   const int numNodes = intrepidBasis->getCardinality();
00129   const int worksetSize = meshSpecs.worksetSize;
00130   
00131   Intrepid::DefaultCubatureFactory<RealType> cubFactory;
00132   RCP <Intrepid::Cubature<RealType> > cubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree);
00133   
00134   const int numQPts = cubature->getNumPoints();
00135   const int numVertices = cellType->getNodeCount();
00136   
00137   *out << "Field Dimensions: Workset=" << worksetSize 
00138        << ", Vertices= " << numVertices
00139        << ", Nodes= " << numNodes
00140        << ", QuadPts= " << numQPts
00141        << ", Dim= " << numDim << std::endl;
00142   
00143    int vecDim = neq;
00144 
00145    dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPts,numDim, vecDim));
00146    Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
00147    bool supportsTransient=true;
00148    int offset=0;
00149 
00150    // Temporary variable used numerous times below
00151    Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
00152 
00153    // Define Field Names
00154 
00155      Teuchos::ArrayRCP<string> dof_names(1);
00156      Teuchos::ArrayRCP<string> dof_names_dot(1);
00157      Teuchos::ArrayRCP<string> resid_names(1);
00158      dof_names[0] = "qFluct";
00159      dof_names_dot[0] = dof_names[0]+"_dot";
00160      resid_names[0] = "ComprNS Residual";
00161      fm0.template registerEvaluator<EvalT>
00162        (evalUtils.constructGatherSolutionEvaluator(true, dof_names, dof_names_dot, offset));
00163 
00164      fm0.template registerEvaluator<EvalT>
00165        (evalUtils.constructDOFVecInterpolationEvaluator(dof_names[0], offset));
00166 
00167      fm0.template registerEvaluator<EvalT>
00168        (evalUtils.constructDOFVecInterpolationEvaluator(dof_names_dot[0], offset));
00169 
00170      //     fm0.template registerEvaluator<EvalT>
00171      //  (evalUtils.constructDOFVecGradInterpolationEvaluator(dof_names[0], offset));
00172 
00173      fm0.template registerEvaluator<EvalT>
00174        (evalUtils.constructScatterResidualEvaluator(true, resid_names,offset, "Scatter ComprNS"));
00175 
00176    fm0.template registerEvaluator<EvalT>
00177      (evalUtils.constructGatherCoordinateVectorEvaluator());
00178 
00179    fm0.template registerEvaluator<EvalT>
00180      (evalUtils.constructMapToPhysicalFrameEvaluator(cellType, cubature));
00181 
00182    fm0.template registerEvaluator<EvalT>
00183      (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cubature));
00184 
00185    { // Specialized DofVecGrad Interpolation for this problem
00186     
00187      RCP<ParameterList> p = rcp(new ParameterList("DOFVecGrad Interpolation "+dof_names[0]));
00188      // Input
00189      p->set<string>("Variable Name", dof_names[0]);
00190      
00191      p->set<string>("Gradient BF Name", "Grad BF");
00192      p->set<int>("Offset of First DOF", offset);
00193      
00194      // Output (assumes same Name as input)
00195      p->set<string>("Gradient Variable Name", dof_names[0]+" Gradient");
00196      
00197      ev = rcp(new PHAL::DOFVecGradInterpolation<EvalT,AlbanyTraits>(*p,dl));
00198      fm0.template registerEvaluator<EvalT>(ev);
00199    }
00200 
00201   { // ComprNS Resid
00202     RCP<ParameterList> p = rcp(new ParameterList("ComprNS Resid"));
00203 
00204     //Input
00205     p->set<string>("Weighted BF Name", "wBF");
00206     p->set<string>("Weighted Gradient BF Name", "wGrad BF");
00207     p->set<string>("QP Variable Name", "qFluct");
00208     p->set<string>("QP Time Derivative Variable Name", "qFluct_dot");
00209     p->set<string>("Gradient QP Variable Name", "qFluct Gradient");
00210     p->set<string>("Body Force Name", "Body Force");   
00211  
00212     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00213     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_vecgradient);
00214     p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
00215     p->set< RCP<DataLayout> >("Node QP Gradient Data Layout", dl->node_qp_gradient);
00216 
00217     p->set<string>("Viscosity Mu QP Variable Name", "Viscosity Mu");
00218     p->set<string>("Viscosity Lambda QP Variable Name", "Viscosity Lambda");
00219     p->set<string>("Viscosity Kappa QP Variable Name", "Viscosity Kappa");
00220     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00221     
00222     p->set<string>("Stress Tau11 QP Variable Name", "Stress Tau11");
00223     p->set<string>("Stress Tau12 QP Variable Name", "Stress Tau12");
00224     p->set<string>("Stress Tau13 QP Variable Name", "Stress Tau13");
00225     p->set<string>("Stress Tau22 QP Variable Name", "Stress Tau22");
00226     p->set<string>("Stress Tau23 QP Variable Name", "Stress Tau23");
00227     p->set<string>("Stress Tau33 QP Variable Name", "Stress Tau33");
00228     
00229     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00230     Teuchos::ParameterList& paramList = params->sublist("Equation Set");
00231     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00232 
00233     //Output
00234     p->set<string>("Residual Name", "ComprNS Residual");
00235     p->set< RCP<DataLayout> >("Node Vector Data Layout", dl->node_vector);
00236 
00237     ev = rcp(new PHAL::ComprNSResid<EvalT,AlbanyTraits>(*p));
00238     fm0.template registerEvaluator<EvalT>(ev);
00239   }
00240   { // Viscosity
00241     RCP<ParameterList> p = rcp(new ParameterList("Viscosity"));
00242 
00243     //Input
00244     p->set<string>("Coordinate Vector Name", "Coord Vec");
00245     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00246     p->set< RCP<DataLayout> >("QP Gradient Data Layout", dl->qp_gradient); 
00247     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_vecgradient);
00248     p->set<string>("QP Variable Name", "qFluct");
00249     p->set<string>("Gradient QP Variable Name", "qFluct Gradient");
00250     
00251     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00252     Teuchos::ParameterList& paramList = params->sublist("Viscosity");
00253     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00254   
00255     //Output
00256     p->set<string>("Viscosity Mu QP Variable Name", "Viscosity Mu");
00257     p->set<string>("Viscosity Lambda QP Variable Name", "Viscosity Lambda");
00258     p->set<string>("Viscosity Kappa QP Variable Name", "Viscosity Kappa");
00259     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00260     p->set<string>("Stress Tau11 QP Variable Name", "Stress Tau11");
00261     p->set<string>("Stress Tau12 QP Variable Name", "Stress Tau12");
00262     p->set<string>("Stress Tau13 QP Variable Name", "Stress Tau13");
00263     p->set<string>("Stress Tau22 QP Variable Name", "Stress Tau22");
00264     p->set<string>("Stress Tau23 QP Variable Name", "Stress Tau23");
00265     p->set<string>("Stress Tau33 QP Variable Name", "Stress Tau33");
00266 
00267     ev = rcp(new PHAL::ComprNSViscosity<EvalT,AlbanyTraits>(*p));
00268     fm0.template registerEvaluator<EvalT>(ev);
00269     
00270   }
00271 
00272    { // ComprNS Body Force
00273     RCP<ParameterList> p = rcp(new ParameterList("Body Force"));
00274 
00275    //Input
00276     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00277     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00278     p->set< RCP<DataLayout> >("QP Gradient Data Layout", dl->qp_gradient);
00279     p->set<string>("Coordinate Vector Name", "Coord Vec");
00280 
00281     Teuchos::ParameterList& paramList = params->sublist("Body Force");
00282     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00283 
00284    //Output
00285      p->set<string>("Body Force Name", "Body Force");
00286 
00287     ev = rcp(new PHAL::ComprNSBodyForce<EvalT,AlbanyTraits>(*p));
00288     fm0.template registerEvaluator<EvalT>(ev);
00289   }
00290 
00291 
00292   if (fieldManagerChoice == Albany::BUILD_RESID_FM)  {
00293     PHX::Tag<typename EvalT::ScalarT> res_tag("Scatter ComprNS", dl->dummy);
00294     fm0.requireField<EvalT>(res_tag);
00295   }
00296   else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
00297     Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
00298     return respUtils.constructResponses(fm0, *responseList, Teuchos::null, stateMgr);
00299   }
00300 
00301   return Teuchos::null;
00302 }
00303 #endif // ALBANY_COMPRNSPROBLEM_HPP