Albany: a Trilinos-based PDE code: Albany_HeatProblem.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_HEATPROBLEM_HPP
00008 #define ALBANY_HEATPROBLEM_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 #include "Albany_ProblemUtils.hpp"
00019 
00020 #include "QCAD_MaterialDatabase.hpp"
00021 
00022 namespace Albany {
00023 
00028   class HeatProblem : public AbstractProblem {
00029   public:
00030   
00032     HeatProblem(const Teuchos::RCP<Teuchos::ParameterList>& params,
00033     const Teuchos::RCP<ParamLib>& paramLib,
00034     const int numDim_,
00035     const Teuchos::RCP<const Epetra_Comm>& comm_);
00036 
00038     ~HeatProblem();
00039 
00041     virtual int spatialDimension() const { return numDim; }
00042 
00044     virtual void buildProblem(
00045       Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> >  meshSpecs,
00046       StateManager& stateMgr);
00047 
00048     // Build evaluators
00049     virtual Teuchos::Array< Teuchos::RCP<const PHX::FieldTag> >
00050     buildEvaluators(
00051       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00052       const Albany::MeshSpecsStruct& meshSpecs,
00053       Albany::StateManager& stateMgr,
00054       Albany::FieldManagerChoice fmchoice,
00055       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00056 
00058     Teuchos::RCP<const Teuchos::ParameterList> getValidProblemParameters() const;
00059 
00060   private:
00061 
00063     HeatProblem(const HeatProblem&);
00064     
00066     HeatProblem& operator=(const HeatProblem&);
00067 
00068   public:
00069 
00071     template <typename EvalT> 
00072     Teuchos::RCP<const PHX::FieldTag>
00073     constructEvaluators(
00074       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00075       const Albany::MeshSpecsStruct& meshSpecs,
00076       Albany::StateManager& stateMgr,
00077       Albany::FieldManagerChoice fmchoice,
00078       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00079 
00080     void constructDirichletEvaluators(const std::vector<std::string>& nodeSetIDs);
00081     void constructNeumannEvaluators(const Teuchos::RCP<Albany::MeshSpecsStruct>& meshSpecs);
00082 
00083   protected:
00084 
00086     bool periodic;
00087     bool haveSource;
00088     bool haveAbsorption;
00089     int numDim;
00090 
00091    Teuchos::RCP<QCAD::MaterialDatabase> materialDB;
00092    Teuchos::RCP<const Epetra_Comm> comm;
00093 
00094    Teuchos::RCP<Albany::Layouts> dl;
00095 
00096   };
00097 
00098 }
00099 
00100 #include "Intrepid_FieldContainer.hpp"
00101 #include "Intrepid_DefaultCubatureFactory.hpp"
00102 #include "Shards_CellTopology.hpp"
00103 #include "Albany_Utils.hpp"
00104 #include "Albany_ProblemUtils.hpp"
00105 #include "Albany_EvaluatorUtils.hpp"
00106 #include "Albany_ResponseUtilities.hpp"
00107 
00108 #include "PHAL_ThermalConductivity.hpp"
00109 #include "PHAL_Absorption.hpp"
00110 #include "PHAL_Source.hpp"
00111 //#include "PHAL_Neumann.hpp"
00112 #include "PHAL_HeatEqResid.hpp"
00113 
00114 
00115 template <typename EvalT>
00116 Teuchos::RCP<const PHX::FieldTag>
00117 Albany::HeatProblem::constructEvaluators(
00118   PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00119   const Albany::MeshSpecsStruct& meshSpecs,
00120   Albany::StateManager& stateMgr,
00121   Albany::FieldManagerChoice fieldManagerChoice,
00122   const Teuchos::RCP<Teuchos::ParameterList>& responseList)
00123 {
00124    using Teuchos::RCP;
00125    using Teuchos::rcp;
00126    using Teuchos::ParameterList;
00127    using PHX::DataLayout;
00128    using PHX::MDALayout;
00129    using std::vector;
00130    using std::string;
00131    using PHAL::AlbanyTraits;
00132 
00133    const CellTopologyData * const elem_top = &meshSpecs.ctd;
00134 
00135    RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > >
00136      intrepidBasis = Albany::getIntrepidBasis(*elem_top);
00137    RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology (elem_top));
00138 
00139 
00140    const int numNodes = intrepidBasis->getCardinality();
00141    const int worksetSize = meshSpecs.worksetSize;
00142 
00143    Intrepid::DefaultCubatureFactory<RealType> cubFactory;
00144    RCP <Intrepid::Cubature<RealType> > cellCubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree);
00145 
00146    const int numQPtsCell = cellCubature->getNumPoints();
00147    const int numVertices = cellType->getNodeCount();
00148 
00149 
00150    *out << "Field Dimensions: Workset=" << worksetSize 
00151         << ", Vertices= " << numVertices
00152         << ", Nodes= " << numNodes
00153         << ", QuadPts= " << numQPtsCell
00154         << ", Dim= " << numDim << std::endl;
00155 
00156    dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPtsCell,numDim));
00157    Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
00158 
00159   // Temporary variable used numerous times below
00160   Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
00161 
00162    Teuchos::ArrayRCP<string> dof_names(neq);
00163      dof_names[0] = "Temperature";
00164    Teuchos::ArrayRCP<string> dof_names_dot(neq);
00165      dof_names_dot[0] = "Temperature_dot";
00166    Teuchos::ArrayRCP<string> resid_names(neq);
00167      resid_names[0] = "Temperature Residual";
00168 
00169   fm0.template registerEvaluator<EvalT>
00170      (evalUtils.constructGatherSolutionEvaluator(false, dof_names, dof_names_dot));
00171 
00172   fm0.template registerEvaluator<EvalT>
00173      (evalUtils.constructScatterResidualEvaluator(false, resid_names));
00174 
00175   fm0.template registerEvaluator<EvalT>
00176     (evalUtils.constructGatherCoordinateVectorEvaluator());
00177 
00178   fm0.template registerEvaluator<EvalT>
00179     (evalUtils.constructMapToPhysicalFrameEvaluator( cellType, cellCubature));
00180 
00181   fm0.template registerEvaluator<EvalT>
00182     (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cellCubature));
00183 
00184   for (unsigned int i=0; i<neq; i++) {
00185     fm0.template registerEvaluator<EvalT>
00186       (evalUtils.constructDOFInterpolationEvaluator(dof_names[i]));
00187 
00188     fm0.template registerEvaluator<EvalT>
00189       (evalUtils.constructDOFInterpolationEvaluator(dof_names_dot[i]));
00190 
00191     fm0.template registerEvaluator<EvalT>
00192       (evalUtils.constructDOFGradInterpolationEvaluator(dof_names[i]));
00193   }
00194 
00195   { // Thermal conductivity
00196     RCP<ParameterList> p = rcp(new ParameterList);
00197 
00198     p->set<string>("QP Variable Name", "Thermal Conductivity");
00199     p->set<string>("QP Coordinate Vector Name", "Coord Vec");
00200     p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00201     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00202     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00203 
00204     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00205     Teuchos::ParameterList& paramList = params->sublist("Thermal Conductivity");
00206     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00207 
00208     // Here we assume that the instance of this problem applies on a single element block
00209     p->set<string>("Element Block Name", meshSpecs.ebName);
00210 
00211     if(materialDB != Teuchos::null)
00212       p->set< RCP<QCAD::MaterialDatabase> >("MaterialDB", materialDB);
00213 
00214     ev = rcp(new PHAL::ThermalConductivity<EvalT,AlbanyTraits>(*p));
00215     fm0.template registerEvaluator<EvalT>(ev);
00216   }
00217 
00218   if (haveAbsorption) { // Absorption
00219     RCP<ParameterList> p = rcp(new ParameterList);
00220 
00221     p->set<string>("QP Variable Name", "Absorption");
00222     p->set<string>("QP Coordinate Vector Name", "Coord Vec");
00223     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00224     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00225 
00226     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00227     Teuchos::ParameterList& paramList = params->sublist("Absorption");
00228     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00229 
00230     ev = rcp(new PHAL::Absorption<EvalT,AlbanyTraits>(*p));
00231     fm0.template registerEvaluator<EvalT>(ev);
00232   }
00233 
00234 // Check and see if a source term is specified for this problem in the main input file. 
00235   bool problemSpecifiesASource = params->isSublist("Source Functions");
00236 
00237   if(problemSpecifiesASource){
00238 
00239       // Sources the same everywhere if they are present at all
00240 
00241       haveSource = true;
00242       RCP<ParameterList> p = rcp(new ParameterList);
00243 
00244       p->set<string>("Source Name", "Source");
00245       p->set<string>("Variable Name", "Temperature");
00246       p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00247 
00248       p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00249       Teuchos::ParameterList& paramList = params->sublist("Source Functions");
00250       p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00251 
00252       ev = rcp(new PHAL::Source<EvalT,AlbanyTraits>(*p));
00253       fm0.template registerEvaluator<EvalT>(ev);
00254 
00255   }
00256   else if(materialDB != Teuchos::null){ // Sources can be specified in terms of materials or element blocks
00257 
00258       // Is the source function active for "this" element block?
00259 
00260       haveSource =  materialDB->isElementBlockSublist(meshSpecs.ebName, "Source Functions");
00261 
00262       if(haveSource){
00263 
00264         RCP<ParameterList> p = rcp(new ParameterList);
00265 
00266         p->set<string>("Source Name", "Source");
00267         p->set<string>("Variable Name", "Temperature");
00268         p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00269 
00270         p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00271         Teuchos::ParameterList& paramList = materialDB->getElementBlockSublist(meshSpecs.ebName, "Source Functions");
00272         p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00273 
00274         ev = rcp(new PHAL::Source<EvalT,AlbanyTraits>(*p));
00275         fm0.template registerEvaluator<EvalT>(ev);
00276     }
00277   }
00278 
00279   { // Temperature Resid
00280     RCP<ParameterList> p = rcp(new ParameterList("Temperature Resid"));
00281 
00282     //Input
00283     p->set<string>("Weighted BF Name", "wBF");
00284     p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
00285     p->set<string>("QP Variable Name", "Temperature");
00286 
00287     p->set<string>("QP Time Derivative Variable Name", "Temperature_dot");
00288 
00289     p->set<bool>("Have Source", haveSource);
00290     p->set<bool>("Have Absorption", haveAbsorption);
00291     p->set<string>("Source Name", "Source");
00292 
00293     p->set<string>("Thermal Conductivity Name", "Thermal Conductivity");
00294     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00295 
00296     p->set<string>("Absorption Name", "Thermal Conductivity");
00297     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00298 
00299     p->set<string>("Gradient QP Variable Name", "Temperature Gradient");
00300     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00301 
00302     p->set<string>("Weighted Gradient BF Name", "wGrad BF");
00303     p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00304     if (params->isType<string>("Convection Velocity"))
00305       p->set<string>("Convection Velocity",
00306                        params->get<string>("Convection Velocity"));
00307     if (params->isType<bool>("Have Rho Cp"))
00308       p->set<bool>("Have Rho Cp", params->get<bool>("Have Rho Cp"));
00309 
00310     //Output
00311     p->set<string>("Residual Name", "Temperature Residual");
00312     p->set< RCP<DataLayout> >("Node Scalar Data Layout", dl->node_scalar);
00313 
00314     ev = rcp(new PHAL::HeatEqResid<EvalT,AlbanyTraits>(*p));
00315     fm0.template registerEvaluator<EvalT>(ev);
00316   }
00317 
00318   if (fieldManagerChoice == Albany::BUILD_RESID_FM)  {
00319     PHX::Tag<typename EvalT::ScalarT> res_tag("Scatter", dl->dummy);
00320     fm0.requireField<EvalT>(res_tag);
00321     return res_tag.clone();
00322   }
00323 
00324   else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
00325     Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
00326     return respUtils.constructResponses(fm0, *responseList, Teuchos::null, stateMgr);
00327   }
00328 
00329   return Teuchos::null;
00330 }
00331 
00332 
00333 #endif // ALBANY_HEATNONLINEARSOURCEPROBLEM_HPP