Albany: a Trilinos-based PDE code: HydMorphProblem.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 HYDRIDEMORPHPROBLEM_HPP
00008 #define HYDRIDEMORPHPROBLEM_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 HydMorphProblem : public AbstractProblem {
00029   public:
00030 
00032     HydMorphProblem(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     ~HydMorphProblem();
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     HydMorphProblem(const HydMorphProblem&);
00064 
00066     HydMorphProblem& operator=(const HydMorphProblem&);
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 
00085     bool haveHeatSource;
00086     int numDim;
00087 
00088     Teuchos::RCP<QCAD::MaterialDatabase> materialDB;
00089     Teuchos::RCP<const Epetra_Comm> comm;
00090 
00091     Teuchos::RCP<Albany::Layouts> dl;
00092 
00093   };
00094 
00095 }
00096 
00097 #include "Intrepid_FieldContainer.hpp"
00098 #include "Intrepid_DefaultCubatureFactory.hpp"
00099 #include "Shards_CellTopology.hpp"
00100 #include "Albany_Utils.hpp"
00101 #include "Albany_ProblemUtils.hpp"
00102 #include "Albany_EvaluatorUtils.hpp"
00103 #include "Albany_ResponseUtilities.hpp"
00104 
00105 #include "PHAL_ThermalConductivity.hpp"
00106 #include "JThermConductivity.hpp"
00107 #include "PHAL_Source.hpp"
00108 //#include "PHAL_Neumann.hpp"
00109 #include "PHAL_HeatEqResid.hpp"
00110 #include "HydFractionResid.hpp"
00111 
00112 
00113 template <typename EvalT>
00114 Teuchos::RCP<const PHX::FieldTag>
00115 Albany::HydMorphProblem::constructEvaluators(
00116   PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00117   const Albany::MeshSpecsStruct& meshSpecs,
00118   Albany::StateManager& stateMgr,
00119   Albany::FieldManagerChoice fieldManagerChoice,
00120   const Teuchos::RCP<Teuchos::ParameterList>& responseList)
00121 {
00122    using Teuchos::RCP;
00123    using Teuchos::rcp;
00124    using Teuchos::ParameterList;
00125    using PHX::DataLayout;
00126    using PHX::MDALayout;
00127    using std::vector;
00128    using std::string;
00129    using PHAL::AlbanyTraits;
00130 
00131    const CellTopologyData * const elem_top = &meshSpecs.ctd;
00132 
00133    RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > >
00134      intrepidBasis = Albany::getIntrepidBasis(*elem_top);
00135    RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology (elem_top));
00136 
00137 
00138    const int numNodes = intrepidBasis->getCardinality();
00139    const int worksetSize = meshSpecs.worksetSize;
00140 
00141    Intrepid::DefaultCubatureFactory<RealType> cubFactory;
00142    RCP <Intrepid::Cubature<RealType> > cellCubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree);
00143 
00144    const int numQPtsCell = cellCubature->getNumPoints();
00145    const int numVertices = cellType->getNodeCount();
00146 
00147 
00148    *out << "Field Dimensions: Workset=" << worksetSize
00149         << ", Vertices= " << numVertices
00150         << ", Nodes= " << numNodes
00151         << ", QuadPts= " << numQPtsCell
00152         << ", Dim= " << numDim << std::endl;
00153 
00154    dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPtsCell,numDim));
00155    Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
00156 
00157   // Temporary variable used numerous times below
00158   Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
00159 
00160 // The coupled heat and hydrogen diffusion equations
00161 
00162   Teuchos::ArrayRCP<std::string> dof_names(neq);
00163   Teuchos::ArrayRCP<std::string> dof_names_dot(neq);
00164   Teuchos::ArrayRCP<std::string> resid_names(neq);
00165 
00166   dof_names[0] = "Temperature";
00167   dof_names_dot[0] = "Temperature_dot";
00168   resid_names[0] = "Temperature Residual";
00169 
00170   dof_names[1] = "HydFraction";
00171   dof_names_dot[1] = "HydFraction_dot";
00172   resid_names[1] = "HydFraction Residual";
00173 
00174   fm0.template registerEvaluator<EvalT>
00175     (evalUtils.constructGatherSolutionEvaluator(false, dof_names, dof_names_dot));
00176 
00177   fm0.template registerEvaluator<EvalT>
00178     (evalUtils.constructScatterResidualEvaluator(false, resid_names));
00179 
00180   fm0.template registerEvaluator<EvalT>
00181     (evalUtils.constructGatherCoordinateVectorEvaluator());
00182 
00183   fm0.template registerEvaluator<EvalT>
00184     (evalUtils.constructMapToPhysicalFrameEvaluator( cellType, cellCubature));
00185 
00186   fm0.template registerEvaluator<EvalT>
00187     (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cellCubature));
00188 
00189   for (unsigned int i=0; i<neq; i++) {
00190     fm0.template registerEvaluator<EvalT>
00191       (evalUtils.constructDOFInterpolationEvaluator(dof_names[i]));
00192 
00193     fm0.template registerEvaluator<EvalT>
00194       (evalUtils.constructDOFInterpolationEvaluator(dof_names_dot[i]));
00195 
00196     fm0.template registerEvaluator<EvalT>
00197       (evalUtils.constructDOFGradInterpolationEvaluator(dof_names[i]));
00198   }
00199 
00200   { // Thermal conductivity
00201     RCP<ParameterList> p = rcp(new ParameterList);
00202 
00203     p->set<std::string>("QP Variable Name", "Thermal Conductivity");
00204     p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00205     p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00206     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00207     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00208 
00209     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00210     Teuchos::ParameterList& paramList = params->sublist("Thermal Conductivity");
00211     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00212 
00213     // Here we assume that the instance of this problem applies on a single element block
00214     p->set<std::string>("Element Block Name", meshSpecs.ebName);
00215 
00216     if(materialDB != Teuchos::null)
00217       p->set< RCP<QCAD::MaterialDatabase> >("MaterialDB", materialDB);
00218 
00219     ev = rcp(new PHAL::ThermalConductivity<EvalT,AlbanyTraits>(*p));
00220     fm0.template registerEvaluator<EvalT>(ev);
00221   }
00222 
00223 // Check and see if a source term is specified for this problem in the main input file.
00224   bool problemSpecifiesASource = params->isSublist("Source Functions");
00225 
00226   if(problemSpecifiesASource){
00227 
00228       // Sources the same everywhere if they are present at all
00229 
00230       haveHeatSource = true;
00231       RCP<ParameterList> p = rcp(new ParameterList);
00232 
00233       p->set<std::string>("Source Name", "Source");
00234       p->set<std::string>("Variable Name", "Temperature");
00235       p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00236 
00237       p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00238       Teuchos::ParameterList& paramList = params->sublist("Source Functions");
00239       p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00240 
00241       ev = rcp(new PHAL::Source<EvalT,AlbanyTraits>(*p));
00242       fm0.template registerEvaluator<EvalT>(ev);
00243 
00244   }
00245   else if(materialDB != Teuchos::null){ // Sources can be specified in terms of materials or element blocks
00246 
00247       // Is the source function active for "this" element block?
00248 
00249       haveHeatSource =  materialDB->isElementBlockSublist(meshSpecs.ebName, "Source Functions");
00250 
00251       if(haveHeatSource){
00252 
00253         RCP<ParameterList> p = rcp(new ParameterList);
00254 
00255         p->set<std::string>("Source Name", "Source");
00256         p->set<std::string>("Variable Name", "Temperature");
00257         p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00258 
00259         p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00260         Teuchos::ParameterList& paramList = materialDB->getElementBlockSublist(meshSpecs.ebName, "Source Functions");
00261         p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00262 
00263         ev = rcp(new PHAL::Source<EvalT,AlbanyTraits>(*p));
00264         fm0.template registerEvaluator<EvalT>(ev);
00265     }
00266   }
00267 
00268   { // Temperature Resid
00269     RCP<ParameterList> p = rcp(new ParameterList("Temperature Resid"));
00270 
00271     //Input
00272     p->set<std::string>("Weighted BF Name", "wBF");
00273     p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
00274     p->set<std::string>("QP Variable Name", "Temperature");
00275 
00276     p->set<std::string>("QP Time Derivative Variable Name", "Temperature_dot");
00277 
00278     p->set<bool>("Have Source", haveHeatSource);
00279     p->set<bool>("Have Absorption", false);
00280     p->set<std::string>("Source Name", "Source");
00281 
00282     p->set<std::string>("Thermal Conductivity Name", "Thermal Conductivity");
00283     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00284 
00285     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00286 
00287     p->set<std::string>("Gradient QP Variable Name", "Temperature Gradient");
00288     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00289 
00290     p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00291     p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00292     if (params->isType<std::string>("Convection Velocity"))
00293       p->set<std::string>("Convection Velocity",
00294                        params->get<std::string>("Convection Velocity"));
00295     if (params->isType<bool>("Have Rho Cp"))
00296       p->set<bool>("Have Rho Cp", params->get<bool>("Have Rho Cp"));
00297 
00298     //Output
00299     p->set<std::string>("Residual Name", "Temperature Residual");
00300     p->set< RCP<DataLayout> >("Node Scalar Data Layout", dl->node_scalar);
00301 
00302     ev = rcp(new PHAL::HeatEqResid<EvalT,AlbanyTraits>(*p));
00303     fm0.template registerEvaluator<EvalT>(ev);
00304   }
00305 
00306   { // The coefficient in front of Grad T
00307     RCP<ParameterList> p = rcp(new ParameterList);
00308 
00309     p->set<std::string>("Temperature Name", "Temperature");
00310     p->set<std::string>("QP Variable Name", "J Conductivity");
00311 
00312     Teuchos::ParameterList& paramList = params->sublist("Material Parameters");
00313     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00314 
00315     // Here we assume that the instance of this problem applies on a single element block
00316     p->set<std::string>("Element Block Name", meshSpecs.ebName);
00317 
00318     if(materialDB != Teuchos::null)
00319       p->set< RCP<QCAD::MaterialDatabase> >("MaterialDB", materialDB);
00320 
00321     ev = rcp(new PHAL::JThermConductivity<EvalT,AlbanyTraits>(*p, dl));
00322     fm0.template registerEvaluator<EvalT>(ev);
00323   }
00324 
00325   { // Hydrogen Concentration Resid
00326     RCP<ParameterList> p = rcp(new ParameterList("Hydrogen Concentration Resid"));
00327 
00328     //Input
00329     p->set<std::string>("Weighted BF Name", "wBF");
00330     p->set<std::string>("Temperature Name", "Temperature");
00331     p->set<std::string>("Temp Time Derivative Name", "Temperature_dot");
00332     p->set<std::string>("QP Time Derivative Variable Name", "HydFraction_dot");
00333     p->set<std::string>("J Conductivity Name", "J Conductivity");
00334     p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00335     p->set<std::string>("QP Variable Name", "HydFraction");
00336 
00337     Teuchos::ParameterList& paramList = params->sublist("Material Parameters");
00338     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00339 
00340     // Here we assume that the instance of this problem applies on a single element block
00341     p->set<std::string>("Element Block Name", meshSpecs.ebName);
00342 
00343     if(materialDB != Teuchos::null)
00344       p->set< RCP<QCAD::MaterialDatabase> >("MaterialDB", materialDB);
00345 
00346     //Output
00347     p->set<std::string>("Residual Name", "HydFraction Residual");
00348 
00349     ev = rcp(new PHAL::HydFractionResid<EvalT,AlbanyTraits>(*p, dl));
00350     fm0.template registerEvaluator<EvalT>(ev);
00351   }
00352 
00353   if (fieldManagerChoice == Albany::BUILD_RESID_FM)  {
00354     PHX::Tag<typename EvalT::ScalarT> res_tag("Scatter", dl->dummy);
00355     fm0.requireField<EvalT>(res_tag);
00356     return res_tag.clone();
00357   }
00358 
00359   else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
00360     Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
00361     return respUtils.constructResponses(fm0, *responseList, Teuchos::null, stateMgr);
00362   }
00363 
00364   return Teuchos::null;
00365 }
00366 
00367 
00368 #endif // ALBANY_HYDMORPHPROBLEM_HPP