Albany: a Trilinos-based PDE code: HydrideProblem.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 HYDRIDEPROBLEM_HPP
00008 #define HYDRIDEPROBLEM_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 namespace Albany {
00021 
00026   class HydrideProblem : public AbstractProblem {
00027   public:
00028   
00030     HydrideProblem(const Teuchos::RCP<Teuchos::ParameterList>& params,
00031     const Teuchos::RCP<ParamLib>& paramLib,
00032     const int numDim_,
00033     const Teuchos::RCP<const Epetra_Comm>& comm_);
00034 
00036     ~HydrideProblem();
00037 
00039     virtual int spatialDimension() const { return numDim; }
00040 
00042     virtual void buildProblem(
00043       Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> >  meshSpecs,
00044       StateManager& stateMgr);
00045 
00046     // Build evaluators
00047     virtual Teuchos::Array< Teuchos::RCP<const PHX::FieldTag> >
00048     buildEvaluators(
00049       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00050       const Albany::MeshSpecsStruct& meshSpecs,
00051       Albany::StateManager& stateMgr,
00052       Albany::FieldManagerChoice fmchoice,
00053       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00054 
00056     Teuchos::RCP<const Teuchos::ParameterList> getValidProblemParameters() const;
00057 
00058   private:
00059 
00061     HydrideProblem(const HydrideProblem&);
00062     
00064     HydrideProblem& operator=(const HydrideProblem&);
00065 
00066   public:
00067 
00069     template <typename EvalT> 
00070     Teuchos::RCP<const PHX::FieldTag>
00071     constructEvaluators(
00072       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00073       const Albany::MeshSpecsStruct& meshSpecs,
00074       Albany::StateManager& stateMgr,
00075       Albany::FieldManagerChoice fmchoice,
00076       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00077 
00078     void constructDirichletEvaluators(const std::vector<std::string>& nodeSetIDs);
00079     void constructNeumannEvaluators(const Teuchos::RCP<Albany::MeshSpecsStruct>& meshSpecs);
00080 
00081   protected:
00082 
00083     int numDim;
00084 
00085     bool haveNoise; // Langevin noise present
00086 
00087     Teuchos::RCP<const Epetra_Comm> comm;
00088 
00089     Teuchos::RCP<Albany::Layouts> dl;
00090 
00091   };
00092 
00093 }
00094 
00095 #include "Intrepid_FieldContainer.hpp"
00096 #include "Intrepid_DefaultCubatureFactory.hpp"
00097 #include "Shards_CellTopology.hpp"
00098 #include "Albany_Utils.hpp"
00099 #include "Albany_ProblemUtils.hpp"
00100 #include "Albany_EvaluatorUtils.hpp"
00101 #include "Albany_ResponseUtilities.hpp"
00102 
00103 #include "HydrideChemTerm.hpp"
00104 #include "HydrideStressTerm.hpp"
00105 #include "HydrideStress.hpp"
00106 #include "PHAL_LangevinNoiseTerm.hpp"
00107 #include "HydrideCResid.hpp"
00108 #include "HydrideWResid.hpp"
00109 
00110 #include "PHAL_SaveStateField.hpp"
00111 #include "Strain.hpp"
00112 #include "ElasticityResid.hpp"
00113 
00114 #include "ElasticModulus.hpp"
00115 #include "PoissonsRatio.hpp"
00116 #include "Stress.hpp"
00117 
00118 
00119 template <typename EvalT>
00120 Teuchos::RCP<const PHX::FieldTag>
00121 Albany::HydrideProblem::constructEvaluators(
00122   PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00123   const Albany::MeshSpecsStruct& meshSpecs,
00124   Albany::StateManager& stateMgr,
00125   Albany::FieldManagerChoice fieldManagerChoice,
00126   const Teuchos::RCP<Teuchos::ParameterList>& responseList)
00127 {
00128    using Teuchos::RCP;
00129    using Teuchos::rcp;
00130    using Teuchos::ParameterList;
00131    using PHX::DataLayout;
00132    using PHX::MDALayout;
00133    using std::vector;
00134    using std::string;
00135    using PHAL::AlbanyTraits;
00136 
00137    // get the name of the current element block
00138    std::string elementBlockName = meshSpecs.ebName;
00139 
00140    const CellTopologyData * const elem_top = &meshSpecs.ctd;
00141 
00142    RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > >
00143      intrepidBasis = Albany::getIntrepidBasis(*elem_top);
00144    RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology (elem_top));
00145 
00146 
00147    const int numNodes = intrepidBasis->getCardinality();
00148    const int worksetSize = meshSpecs.worksetSize;
00149 
00150    Intrepid::DefaultCubatureFactory<RealType> cubFactory;
00151    RCP <Intrepid::Cubature<RealType> > cellCubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree);
00152 
00153    const int numQPtsCell = cellCubature->getNumPoints();
00154    const int numVertices = cellType->getNodeCount();
00155 
00156 
00157    *out << "Field Dimensions: Workset=" << worksetSize 
00158         << ", Vertices= " << numVertices
00159         << ", Nodes= " << numNodes
00160         << ", QuadPts= " << numQPtsCell
00161         << ", Dim= " << numDim << std::endl;
00162 
00163    dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPtsCell,numDim));
00164    Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
00165 
00166   // Temporary variable used numerous times below
00167   Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
00168 
00169   int offset=0;
00170 
00171 // The displacement equations
00172 
00173   {
00174      Teuchos::ArrayRCP<std::string> dof_names(1);
00175      Teuchos::ArrayRCP<std::string> dof_names_dot(1);
00176      Teuchos::ArrayRCP<std::string> resid_names(1);
00177      dof_names[0] = "Displacement";
00178      dof_names_dot[0] = dof_names[0]+"DotDot";
00179      resid_names[0] = "Displacement Residual";
00180 
00181      fm0.template registerEvaluator<EvalT>
00182        (evalUtils.constructGatherSolutionEvaluator(true, dof_names, dof_names_dot, offset));
00183 
00184      fm0.template registerEvaluator<EvalT>
00185        (evalUtils.constructDOFVecInterpolationEvaluator(dof_names[0], offset));
00186 
00187      fm0.template registerEvaluator<EvalT>
00188        (evalUtils.constructDOFVecInterpolationEvaluator(dof_names_dot[0], offset));
00189 
00190      fm0.template registerEvaluator<EvalT>
00191        (evalUtils.constructDOFVecGradInterpolationEvaluator(dof_names[0], offset));
00192 
00193      fm0.template registerEvaluator<EvalT>
00194        (evalUtils.constructScatterResidualEvaluator(true, resid_names, offset, "Scatter Displacement"));
00195 
00196      offset += numDim;
00197 
00198    }
00199 
00200 // Equations for c and w
00201 
00202    {
00203      int nscalars = 2;
00204      Teuchos::ArrayRCP<std::string> dof_names(nscalars);
00205        dof_names[0] = "C"; // The concentration difference variable 0 \leq C \leq 1
00206        dof_names[1] = "W"; // The chemical potential difference variable
00207      Teuchos::ArrayRCP<std::string> dof_names_dot(nscalars);
00208        dof_names_dot[0] = dof_names[0]+"Dot";
00209        dof_names_dot[1] = dof_names[1]+"Dot"; // not currently used
00210      Teuchos::ArrayRCP<std::string> resid_names(nscalars);
00211        resid_names[0] = "C Residual";
00212        resid_names[1] = "W Residual";
00213   
00214     fm0.template registerEvaluator<EvalT>
00215        (evalUtils.constructGatherSolutionEvaluator(false, dof_names, dof_names_dot, offset));
00216   
00217     fm0.template registerEvaluator<EvalT>
00218        (evalUtils.constructScatterResidualEvaluator(false, resid_names, offset, "Scatter c and w"));
00219   
00220     for (unsigned int i=0; i<nscalars; i++) {
00221       fm0.template registerEvaluator<EvalT>
00222         (evalUtils.constructDOFInterpolationEvaluator(dof_names[i], offset));
00223   
00224       fm0.template registerEvaluator<EvalT>
00225         (evalUtils.constructDOFInterpolationEvaluator(dof_names_dot[i], offset));
00226   
00227       fm0.template registerEvaluator<EvalT>
00228         (evalUtils.constructDOFGradInterpolationEvaluator(dof_names[i], offset));
00229     }
00230 
00231     offset += nscalars;
00232 
00233   }
00234 
00235   fm0.template registerEvaluator<EvalT>
00236     (evalUtils.constructGatherCoordinateVectorEvaluator());
00237 
00238   fm0.template registerEvaluator<EvalT>
00239     (evalUtils.constructMapToPhysicalFrameEvaluator( cellType, cellCubature));
00240 
00241   fm0.template registerEvaluator<EvalT>
00242     (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cellCubature));
00243 
00244 
00245   { // Form the Chemical Energy term in Eq. 2.2
00246 
00247     RCP<ParameterList> p = rcp(new ParameterList("Chem Energy Term"));
00248 
00249     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00250 
00251     // b value in Equation 1.1
00252     p->set<double>("b Value", params->get<double>("b"));
00253 
00254     //Input
00255     p->set<std::string>("C QP Variable Name", "C");
00256     p->set<std::string>("W QP Variable Name", "W");
00257 
00258     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00259     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00260 
00261     //Output
00262     p->set<std::string>("Chemical Energy Term", "Chemical Energy Term");
00263 
00264     ev = rcp(new HYD::HydrideChemTerm<EvalT,AlbanyTraits>(*p));
00265     fm0.template registerEvaluator<EvalT>(ev);
00266   }
00267 
00268   { // Form the Stress term in Eq. 2.2
00269 
00270     RCP<ParameterList> p = rcp(new ParameterList("Stress Term"));
00271 
00272     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00273 
00274     // b value in Equation 1.1
00275     p->set<double>("e Value", params->get<double>("e"));
00276 
00277     //Input
00278     p->set<std::string>("Stress Name", "Stress"); 
00279     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00280 
00281     //Output
00282     p->set<std::string>("Stress Term", "Stress Term");
00283     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00284 
00285     ev = rcp(new HYD::HydrideStressTerm<EvalT,AlbanyTraits>(*p));
00286     fm0.template registerEvaluator<EvalT>(ev);
00287   }
00288 
00289   if(params->isParameter("Langevin Noise SD")){
00290 
00291    // Form the Langevin noise term
00292 
00293     haveNoise = true;
00294 
00295     RCP<ParameterList> p = rcp(new ParameterList("Langevin Noise Term"));
00296 
00297     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00298 
00299     // Standard deviation of the noise
00300     p->set<double>("SD Value", params->get<double>("Langevin Noise SD"));
00301     // Time period over which to apply the noise (-1 means over the whole time)
00302     p->set<Teuchos::Array<int> >("Langevin Noise Time Period", 
00303         params->get<Teuchos::Array<int> >("Langevin Noise Time Period", Teuchos::tuple<int>(-1, -1)));
00304 
00305     //Input
00306     p->set<std::string>("C QP Variable Name", "C");
00307 
00308     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00309     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00310 
00311     //Output
00312     p->set<std::string>("Langevin Noise Term", "Langevin Noise Term");
00313 
00314     ev = rcp(new PHAL::LangevinNoiseTerm<EvalT,AlbanyTraits>(*p));
00315     fm0.template registerEvaluator<EvalT>(ev);
00316   }
00317 
00318   { // Strain
00319     RCP<ParameterList> p = rcp(new ParameterList("Strain"));
00320 
00321     //Input
00322     p->set<std::string>("Gradient QP Variable Name", "Displacement Gradient");
00323     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00324 
00325     //Output
00326     p->set<std::string>("Strain Name", "Strain"); //dl->qp_tensor also
00327 
00328     ev = rcp(new LCM::Strain<EvalT,AlbanyTraits>(*p, dl));
00329     fm0.template registerEvaluator<EvalT>(ev);
00330   }
00331 
00332 
00333 #if 0
00334   { // Hydride stress
00335       RCP<ParameterList> p = rcp(new ParameterList("HydrideStress"));
00336 
00337       //Input
00338       p->set<std::string>("Strain Name", "Strain");
00339       p->set<std::string>("C QP Variable Name", "C");
00340       p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00341       p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00342 
00343       // e value in Equation between 1.2 and 1.3
00344       p->set<double>("e Value", params->get<double>("e"));
00345 
00346       //Output
00347       p->set<std::string>("Stress Name", "Stress"); //dl->qp_tensor also
00348 
00349       ev = rcp(new HYD::HydrideStress<EvalT,AlbanyTraits>(*p));
00350       fm0.template registerEvaluator<EvalT>(ev);
00351       p = stateMgr.registerStateVariable("Stress",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0);
00352       ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00353       fm0.template registerEvaluator<EvalT>(ev);
00354 
00355   }
00356 #else
00357   { // Elastic Modulus
00358     RCP<ParameterList> p = rcp(new ParameterList);
00359 
00360     p->set<std::string>("QP Variable Name", "Elastic Modulus");
00361     p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00362     p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00363     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00364     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00365 
00366     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00367     Teuchos::ParameterList& paramList = params->sublist("Elastic Modulus");
00368     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00369 
00370     ev = rcp(new LCM::ElasticModulus<EvalT,AlbanyTraits>(*p));
00371     fm0.template registerEvaluator<EvalT>(ev);
00372   }
00373 
00374   { // Poissons Ratio 
00375     RCP<ParameterList> p = rcp(new ParameterList);
00376 
00377     p->set<std::string>("QP Variable Name", "Poissons Ratio");
00378     p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00379     p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00380     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00381     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00382 
00383     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00384     Teuchos::ParameterList& paramList = params->sublist("Poissons Ratio");
00385     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00386 
00387     ev = rcp(new LCM::PoissonsRatio<EvalT,AlbanyTraits>(*p));
00388     fm0.template registerEvaluator<EvalT>(ev);
00389   }
00390   { // Linear elasticity stress
00391       RCP<ParameterList> p = rcp(new ParameterList("Stress"));
00392 
00393       //Input
00394       p->set<std::string>("Strain Name", "Strain");
00395       p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00396 
00397       p->set<std::string>("Elastic Modulus Name", "Elastic Modulus");
00398       p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00399 
00400       p->set<std::string>("Poissons Ratio Name", "Poissons Ratio");  // dl->qp_scalar also
00401 
00402       //Output
00403       p->set<std::string>("Stress Name", "Stress"); //dl->qp_tensor also
00404 
00405       ev = rcp(new LCM::Stress<EvalT,AlbanyTraits>(*p));
00406       fm0.template registerEvaluator<EvalT>(ev);
00407       p = stateMgr.registerStateVariable("Stress",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0);
00408       ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00409       fm0.template registerEvaluator<EvalT>(ev);
00410   }
00411 #endif
00412 
00413 
00414   { // C Resid
00415     RCP<ParameterList> p = rcp(new ParameterList("C Resid"));
00416 
00417     //Input
00418     p->set<std::string>("Weighted BF Name", "wBF");
00419     p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00420     if(haveNoise)
00421       p->set<std::string>("Langevin Noise Term", "Langevin Noise Term");
00422     // Accumulate in the Langevin noise term?
00423     p->set<bool>("Have Noise", haveNoise);
00424 
00425     p->set<std::string>("Chemical Energy Term", "Chemical Energy Term");
00426     p->set<std::string>("Gradient QP Variable Name", "C Gradient");
00427     p->set<std::string>("Stress Term", "Stress Term");
00428 
00429     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00430     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00431     p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
00432     p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00433 
00434     // gamma value in Equation 2.2
00435     p->set<double>("gamma Value", params->get<double>("gamma"));
00436 
00437     //Output
00438     p->set<std::string>("Residual Name", "C Residual");
00439     p->set< RCP<DataLayout> >("Node Scalar Data Layout", dl->node_scalar);
00440 
00441     ev = rcp(new HYD::HydrideCResid<EvalT,AlbanyTraits>(*p));
00442     fm0.template registerEvaluator<EvalT>(ev);
00443   }
00444 
00445   { // W Resid
00446     RCP<ParameterList> p = rcp(new ParameterList("W Resid"));
00447 
00448     //Input
00449     p->set<std::string>("Weighted BF Name", "wBF");
00450     p->set<std::string>("BF Name", "BF");
00451     p->set<std::string>("C QP Time Derivative Variable Name", "CDot");
00452     p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00453     p->set<std::string>("Gradient QP Variable Name", "W Gradient");
00454 
00455     // Mass lump time term?
00456     p->set<bool>("Lump Mass", params->get<bool>("Lump Mass"));
00457 
00458     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00459     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00460     p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
00461     p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00462 
00463     //Output
00464     p->set<std::string>("Residual Name", "W Residual");
00465     p->set< RCP<DataLayout> >("Node Scalar Data Layout", dl->node_scalar);
00466 
00467     ev = rcp(new HYD::HydrideWResid<EvalT,AlbanyTraits>(*p));
00468     fm0.template registerEvaluator<EvalT>(ev);
00469   }
00470 
00471   { // Displacement Resid
00472     RCP<ParameterList> p = rcp(new ParameterList("Displacement Resid"));
00473 
00474     //Input
00475     p->set<std::string>("Stress Name", "Stress");
00476     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00477 
00478     p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00479     p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00480     p->set<bool>("Disable Transient", true);
00481 
00482     //Output
00483     p->set<std::string>("Residual Name", "Displacement Residual");
00484     p->set< RCP<DataLayout> >("Node Vector Data Layout", dl->node_vector);
00485 
00486     ev = rcp(new LCM::ElasticityResid<EvalT,AlbanyTraits>(*p));
00487     fm0.template registerEvaluator<EvalT>(ev);
00488   }
00489 
00490   if (fieldManagerChoice == Albany::BUILD_RESID_FM)  {
00491     PHX::Tag<typename EvalT::ScalarT> disp_tag("Scatter Displacement", dl->dummy);
00492     fm0.requireField<EvalT>(disp_tag);
00493     PHX::Tag<typename EvalT::ScalarT> c_w_tag("Scatter c and w", dl->dummy);
00494     fm0.requireField<EvalT>(c_w_tag);
00495     return disp_tag.clone();
00496   }
00497 
00498   else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
00499     Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
00500     return respUtils.constructResponses(fm0, *responseList, Teuchos::null, stateMgr);
00501   }
00502 
00503   return Teuchos::null;
00504 }
00505 
00506 
00507 #endif // ALBANY_CAHNHILLPROBLEM_HPP