Albany: a Trilinos-based PDE code: Albany_ThermoElectrostaticsProblem.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_THERMOELECTROSTATICSPROBLEM_HPP
00008 #define ALBANY_THERMOELECTROSTATICSPROBLEM_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 ThermoElectrostaticsProblem : public AbstractProblem {
00026   public:
00027   
00029     ThermoElectrostaticsProblem(
00030       const Teuchos::RCP<Teuchos::ParameterList>& params,
00031       const Teuchos::RCP<ParamLib>& paramLib,
00032       const int numDim_);
00033 
00035     virtual int spatialDimension() const { return numDim; }
00036 
00038     ~ThermoElectrostaticsProblem();
00039 
00041     virtual void buildProblem(
00042       Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> >  meshSpecs,
00043       StateManager& stateMgr);
00044 
00045     // Build evaluators
00046     virtual Teuchos::Array< Teuchos::RCP<const PHX::FieldTag> >
00047     buildEvaluators(
00048       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00049       const Albany::MeshSpecsStruct& meshSpecs,
00050       Albany::StateManager& stateMgr,
00051       Albany::FieldManagerChoice fmchoice,
00052       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00053 
00055     Teuchos::RCP<const Teuchos::ParameterList> getValidProblemParameters() const;
00056 
00057   private:
00058 
00060     ThermoElectrostaticsProblem(const ThermoElectrostaticsProblem&);
00061     
00063     ThermoElectrostaticsProblem& operator=(const ThermoElectrostaticsProblem&);
00064 
00065   public:
00066 
00068     template <typename EvalT> Teuchos::RCP<const PHX::FieldTag>
00069     constructEvaluators(
00070       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00071       const Albany::MeshSpecsStruct& meshSpecs,
00072       Albany::StateManager& stateMgr,
00073       Albany::FieldManagerChoice fmchoice,
00074       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00075 
00076     void constructDirichletEvaluators(const Albany::MeshSpecsStruct& meshSpecs);
00077 
00078   protected:
00079 
00081     int numDim;
00082   };
00083 
00084 }
00085 
00086 #include "Intrepid_FieldContainer.hpp"
00087 #include "Intrepid_DefaultCubatureFactory.hpp"
00088 #include "Shards_CellTopology.hpp"
00089 #include "Albany_Utils.hpp"
00090 #include "Albany_ProblemUtils.hpp"
00091 #include "Albany_EvaluatorUtils.hpp"
00092 #include "Albany_ResponseUtilities.hpp"
00093 
00094 #ifdef ALBANY_QCAD
00095 #include "PHAL_TEProp.hpp"
00096 #include "PHAL_JouleHeating.hpp"
00097 #include "QCAD_PoissonResid.hpp"
00098 #endif
00099 #include "PHAL_HeatEqResid.hpp"
00100 
00101 
00102 
00103 template <typename EvalT>
00104 Teuchos::RCP<const PHX::FieldTag>
00105 Albany::ThermoElectrostaticsProblem::constructEvaluators(
00106   PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00107   const Albany::MeshSpecsStruct& meshSpecs,
00108   Albany::StateManager& stateMgr,
00109   Albany::FieldManagerChoice fieldManagerChoice,
00110   const Teuchos::RCP<Teuchos::ParameterList>& responseList)
00111 {
00112    using Teuchos::RCP;
00113    using Teuchos::rcp;
00114    using Teuchos::ParameterList;
00115    using PHX::DataLayout;
00116    using PHX::MDALayout;
00117    using std::vector;
00118    using std::string;
00119    using PHAL::AlbanyTraits;
00120 
00121    RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology (&meshSpecs.ctd));
00122    RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > >
00123      intrepidBasis = Albany::getIntrepidBasis(meshSpecs.ctd);
00124 
00125    const int numNodes = intrepidBasis->getCardinality();
00126    const int worksetSize = meshSpecs.worksetSize;
00127 
00128    Intrepid::DefaultCubatureFactory<RealType> cubFactory;
00129    RCP <Intrepid::Cubature<RealType> > cubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree);
00130 
00131    const int numQPts = cubature->getNumPoints();
00132    const int numVertices = cellType->getNodeCount();
00133 
00134    *out << "Field Dimensions: Workset=" << worksetSize 
00135         << ", Vertices= " << numVertices
00136         << ", Nodes= " << numNodes
00137         << ", QuadPts= " << numQPts
00138         << ", Dim= " << numDim << std::endl;
00139 
00140 
00141    RCP<Albany::Layouts> dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPts,numDim));
00142    Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
00143    bool supportsTransient=false;
00144 
00145    // Temporary variable used numerous times below
00146    Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
00147 
00148    // Define Field Names
00149 
00150    Teuchos::ArrayRCP<string> dof_names(neq);
00151      dof_names[0] = "Potential";
00152      dof_names[1] = "Temperature";
00153 
00154    Teuchos::ArrayRCP<string> dof_names_dot(neq);
00155    if (supportsTransient) {
00156      for (int i=0; i<neq; i++) dof_names_dot[i] = dof_names[i]+"_dot";
00157    }
00158 
00159    Teuchos::ArrayRCP<string> resid_names(neq);
00160      for (int i=0; i<neq; i++) resid_names[i] = dof_names[i]+" Residual";
00161 
00162    if (supportsTransient) fm0.template registerEvaluator<EvalT>
00163        (evalUtils.constructGatherSolutionEvaluator(false, dof_names, dof_names_dot));
00164    else fm0.template registerEvaluator<EvalT>
00165        (evalUtils.constructGatherSolutionEvaluator_noTransient(false, dof_names));
00166 
00167    fm0.template registerEvaluator<EvalT>
00168      (evalUtils.constructScatterResidualEvaluator(false, resid_names));
00169 
00170    fm0.template registerEvaluator<EvalT>
00171      (evalUtils.constructGatherCoordinateVectorEvaluator());
00172 
00173    fm0.template registerEvaluator<EvalT>
00174      (evalUtils.constructMapToPhysicalFrameEvaluator(cellType, cubature));
00175 
00176    fm0.template registerEvaluator<EvalT>
00177      (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cubature));
00178 
00179    for (int i=0; i<neq; i++) {
00180      fm0.template registerEvaluator<EvalT>
00181        (evalUtils.constructDOFInterpolationEvaluator(dof_names[i], i));
00182 
00183      if (supportsTransient)
00184      fm0.template registerEvaluator<EvalT>
00185          (evalUtils.constructDOFInterpolationEvaluator(dof_names_dot[i], i));
00186 
00187      fm0.template registerEvaluator<EvalT>
00188          (evalUtils.constructDOFGradInterpolationEvaluator(dof_names[i], i));
00189   }
00190 
00191 
00192   { // Thermal conductivity
00193     RCP<ParameterList> p = rcp(new ParameterList);
00194 
00195     p->set<string>("QP Variable Name", "Thermal Conductivity");
00196     p->set<string>("QP Variable Name 2", "Permittivity");  // really electrical conductivity
00197     p->set<string>("QP Variable Name 3", "Rho Cp"); 
00198     p->set<string>("Coordinate Vector Name", "Coord Vec");
00199     p->set<string>("Temperature Variable Name", "Temperature");
00200     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00201     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00202 
00203     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00204     Teuchos::ParameterList& paramList = params->sublist("TE Properties");
00205     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00206 
00207     ev = rcp(new PHAL::TEProp<EvalT,AlbanyTraits>(*p));
00208     fm0.template registerEvaluator<EvalT>(ev);
00209   }
00210 
00211   {
00212     RCP<ParameterList> p = rcp(new ParameterList);
00213 
00214     //Input
00215     p->set<string>("Gradient Variable Name", "Potential Gradient");
00216     p->set<string>("Flux Variable Name", "Potential Flux");
00217     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00218 
00219     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00220 
00221     //Output
00222     p->set<string>("Source Name", "Joule");
00223     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00224 
00225     ev = rcp(new PHAL::JouleHeating<EvalT,AlbanyTraits>(*p));
00226     fm0.template registerEvaluator<EvalT>(ev);
00227   }
00228 
00229   { // Potential Resid
00230     RCP<ParameterList> p = rcp(new ParameterList("Potential Resid"));
00231 
00232     //Input
00233     p->set<string>("Weighted BF Name", "wBF");
00234     p->set<string>("QP Variable Name", "Potential");
00235 
00236     p->set<string>("Permittivity Name", "Permittivity");
00237 
00238     p->set<string>("Gradient QP Variable Name", "Potential Gradient");
00239     p->set<string>("Flux QP Variable Name", "Potential Flux");
00240 
00241     p->set<string>("Weighted Gradient BF Name", "wGrad BF");
00242 
00243     p->set<bool>("Have Source", false);
00244     p->set<string>("Source Name", "None");
00245 
00246     //Output
00247     p->set<string>("Residual Name", "Potential Residual");
00248 
00249     ev = rcp(new QCAD::PoissonResid<EvalT,AlbanyTraits>(*p,dl));
00250     fm0.template registerEvaluator<EvalT>(ev);
00251   }
00252 
00253   { // Temperature Resid
00254     RCP<ParameterList> p = rcp(new ParameterList("Temperature Resid"));
00255 
00256     //Input
00257     p->set<string>("Weighted BF Name", "wBF");
00258     p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
00259     p->set<string>("QP Variable Name", "Temperature");
00260 
00261     p->set<bool>("Have Source", true);
00262     p->set<string>("Source Name", "Joule");
00263 
00264     p->set<bool>("Have Absorption", false);
00265 
00266     p->set<string>("Thermal Conductivity Name", "Thermal Conductivity");
00267     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00268 
00269     p->set<string>("Gradient QP Variable Name", "Temperature Gradient");
00270     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00271 
00272     p->set<string>("Weighted Gradient BF Name", "wGrad BF");
00273     p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00274  
00275     // Poisson solve does not have transient terms
00276     p->set<bool>("Disable Transient", true);
00277     p->set<string>("QP Time Derivative Variable Name", "Temperature_dot");
00278 
00279     if (params->isType<string>("Convection Velocity")) {
00280       p->set<string>("Convection Velocity",params->get<string>("Convection Velocity"));
00281       p->set<string>("Rho Cp Name", "Rho Cp"); 
00282     }
00283 
00284     //Output
00285     p->set<string>("Residual Name", "Temperature Residual");
00286     p->set< RCP<DataLayout> >("Node Scalar Data Layout", dl->node_scalar);
00287 
00288     ev = rcp(new PHAL::HeatEqResid<EvalT,AlbanyTraits>(*p));
00289     fm0.template registerEvaluator<EvalT>(ev);
00290   }
00291 
00292   if (fieldManagerChoice == Albany::BUILD_RESID_FM)  {
00293     PHX::Tag<typename EvalT::ScalarT> res_tag("Scatter", dl->dummy);
00294     fm0.requireField<EvalT>(res_tag);
00295     return res_tag.clone();
00296   }
00297 
00298   else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
00299     Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
00300     return respUtils.constructResponses(fm0, *responseList, Teuchos::null, stateMgr);
00301   }
00302 
00303   return Teuchos::null;
00304 }
00305 #endif // ALBANY_HEATNONLINEARSOURCEPROBLEM_HPP