Albany: a Trilinos-based PDE code: Albany_CahnHillProblem.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_CAHNHILLPROBLEM_HPP
00008 #define ALBANY_CAHNHILLPROBLEM_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 CahnHillProblem : public AbstractProblem {
00027   public:
00028   
00030     CahnHillProblem(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     ~CahnHillProblem();
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     CahnHillProblem(const CahnHillProblem&);
00062     
00064     CahnHillProblem& operator=(const CahnHillProblem&);
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 
00080   protected:
00081 
00082     int numDim;
00083 
00084     bool haveNoise; // Langevin noise present
00085 
00086     Teuchos::RCP<const Epetra_Comm> comm;
00087 
00088     Teuchos::RCP<Albany::Layouts> dl;
00089 
00090   };
00091 
00092 }
00093 
00094 #include "Intrepid_FieldContainer.hpp"
00095 #include "Intrepid_DefaultCubatureFactory.hpp"
00096 #include "Shards_CellTopology.hpp"
00097 #include "Albany_Utils.hpp"
00098 #include "Albany_ProblemUtils.hpp"
00099 #include "Albany_EvaluatorUtils.hpp"
00100 #include "Albany_ResponseUtilities.hpp"
00101 
00102 #include "PHAL_CahnHillChemTerm.hpp"
00103 #include "PHAL_LangevinNoiseTerm.hpp"
00104 #include "PHAL_CahnHillRhoResid.hpp"
00105 #include "PHAL_CahnHillWResid.hpp"
00106 
00107 
00108 template <typename EvalT>
00109 Teuchos::RCP<const PHX::FieldTag>
00110 Albany::CahnHillProblem::constructEvaluators(
00111   PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00112   const Albany::MeshSpecsStruct& meshSpecs,
00113   Albany::StateManager& stateMgr,
00114   Albany::FieldManagerChoice fieldManagerChoice,
00115   const Teuchos::RCP<Teuchos::ParameterList>& responseList)
00116 {
00117    using Teuchos::RCP;
00118    using Teuchos::rcp;
00119    using Teuchos::ParameterList;
00120    using PHX::DataLayout;
00121    using PHX::MDALayout;
00122    using std::vector;
00123    using std::string;
00124    using PHAL::AlbanyTraits;
00125 
00126    const CellTopologyData * const elem_top = &meshSpecs.ctd;
00127 
00128    RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > >
00129      intrepidBasis = Albany::getIntrepidBasis(*elem_top);
00130    RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology (elem_top));
00131 
00132 
00133    const int numNodes = intrepidBasis->getCardinality();
00134    const int worksetSize = meshSpecs.worksetSize;
00135 
00136    Intrepid::DefaultCubatureFactory<RealType> cubFactory;
00137    RCP <Intrepid::Cubature<RealType> > cellCubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree);
00138 
00139    const int numQPtsCell = cellCubature->getNumPoints();
00140    const int numVertices = cellType->getNodeCount();
00141 
00142 
00143    *out << "Field Dimensions: Workset=" << worksetSize 
00144         << ", Vertices= " << numVertices
00145         << ", Nodes= " << numNodes
00146         << ", QuadPts= " << numQPtsCell
00147         << ", Dim= " << numDim << std::endl;
00148 
00149    dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPtsCell,numDim));
00150    Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
00151 
00152   // Temporary variable used numerous times below
00153   Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
00154 
00155    Teuchos::ArrayRCP<string> dof_names(neq);
00156      dof_names[0] = "Rho"; // The concentration difference variable 0 \leq \rho \leq 1
00157      dof_names[1] = "W"; // The chemical potential difference variable
00158    Teuchos::ArrayRCP<string> dof_names_dot(neq);
00159      dof_names_dot[0] = "rhoDot";
00160      dof_names_dot[1] = "wDot"; // not currently used
00161    Teuchos::ArrayRCP<string> resid_names(neq);
00162      resid_names[0] = "Rho Residual";
00163      resid_names[1] = "W Residual";
00164 
00165   fm0.template registerEvaluator<EvalT>
00166      (evalUtils.constructGatherSolutionEvaluator(false, dof_names, dof_names_dot));
00167 
00168   fm0.template registerEvaluator<EvalT>
00169      (evalUtils.constructScatterResidualEvaluator(false, resid_names));
00170 
00171   fm0.template registerEvaluator<EvalT>
00172     (evalUtils.constructGatherCoordinateVectorEvaluator());
00173 
00174   fm0.template registerEvaluator<EvalT>
00175     (evalUtils.constructMapToPhysicalFrameEvaluator( cellType, cellCubature));
00176 
00177   fm0.template registerEvaluator<EvalT>
00178     (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cellCubature));
00179 
00180   for (unsigned int i=0; i<neq; i++) {
00181     fm0.template registerEvaluator<EvalT>
00182       (evalUtils.constructDOFInterpolationEvaluator(dof_names[i]));
00183 
00184     fm0.template registerEvaluator<EvalT>
00185       (evalUtils.constructDOFInterpolationEvaluator(dof_names_dot[i]));
00186 
00187     fm0.template registerEvaluator<EvalT>
00188       (evalUtils.constructDOFGradInterpolationEvaluator(dof_names[i]));
00189   }
00190 
00191 
00192   { // Form the Chemical Energy term in Eq. 2.2
00193 
00194     RCP<ParameterList> p = rcp(new ParameterList("Chem Energy Term"));
00195 
00196     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00197 
00198     // b value in Equation 1.1
00199     p->set<double>("b Value", params->get<double>("b"));
00200 
00201     //Input
00202     p->set<string>("Rho QP Variable Name", "Rho");
00203     p->set<string>("W QP Variable Name", "W");
00204 
00205     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00206     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00207 
00208     //Output
00209     p->set<string>("Chemical Energy Term", "Chemical Energy Term");
00210 
00211     ev = rcp(new PHAL::CahnHillChemTerm<EvalT,AlbanyTraits>(*p));
00212     fm0.template registerEvaluator<EvalT>(ev);
00213   }
00214 
00215   if(params->isParameter("Langevin Noise SD")){
00216 
00217    // Form the Langevin noise term
00218 
00219     haveNoise = true;
00220 
00221     RCP<ParameterList> p = rcp(new ParameterList("Langevin Noise Term"));
00222 
00223     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00224 
00225     // Standard deviation of the noise
00226     p->set<double>("SD Value", params->get<double>("Langevin Noise SD"));
00227     // Time period over which to apply the noise (-1 means over the whole time)
00228     p->set<Teuchos::Array<int> >("Langevin Noise Time Period", 
00229         params->get<Teuchos::Array<int> >("Langevin Noise Time Period", Teuchos::tuple<int>(-1, -1)));
00230 
00231     //Input
00232     p->set<string>("Rho QP Variable Name", "Rho");
00233 
00234     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00235     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00236 
00237     //Output
00238     p->set<string>("Langevin Noise Term", "Langevin Noise Term");
00239 
00240     ev = rcp(new PHAL::LangevinNoiseTerm<EvalT,AlbanyTraits>(*p));
00241     fm0.template registerEvaluator<EvalT>(ev);
00242   }
00243 
00244   { // Rho Resid
00245     RCP<ParameterList> p = rcp(new ParameterList("Rho Resid"));
00246 
00247     //Input
00248     p->set<string>("Weighted BF Name", "wBF");
00249     p->set<string>("Weighted Gradient BF Name", "wGrad BF");
00250     if(haveNoise)
00251       p->set<string>("Langevin Noise Term", "Langevin Noise Term");
00252     // Accumulate in the Langevin noise term?
00253     p->set<bool>("Have Noise", haveNoise);
00254 
00255     p->set<string>("Chemical Energy Term", "Chemical Energy Term");
00256     p->set<string>("Gradient QP Variable Name", "Rho Gradient");
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     p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
00261     p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00262 
00263     // gamma value in Equation 2.2
00264     p->set<double>("gamma Value", params->get<double>("gamma"));
00265 
00266     //Output
00267     p->set<string>("Residual Name", "Rho Residual");
00268     p->set< RCP<DataLayout> >("Node Scalar Data Layout", dl->node_scalar);
00269 
00270     ev = rcp(new PHAL::CahnHillRhoResid<EvalT,AlbanyTraits>(*p));
00271     fm0.template registerEvaluator<EvalT>(ev);
00272   }
00273 
00274   { // W Resid
00275     RCP<ParameterList> p = rcp(new ParameterList("W Resid"));
00276 
00277     //Input
00278     p->set<string>("Weighted BF Name", "wBF");
00279     p->set<string>("BF Name", "BF");
00280     p->set<string>("Rho QP Time Derivative Variable Name", "rhoDot");
00281     p->set<string>("Weighted Gradient BF Name", "wGrad BF");
00282     p->set<string>("Gradient QP Variable Name", "W Gradient");
00283 
00284     // Mass lump time term?
00285     p->set<bool>("Lump Mass", params->get<bool>("Lump Mass"));
00286 
00287     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00288     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00289     p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
00290     p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00291 
00292     //Output
00293     p->set<string>("Residual Name", "W Residual");
00294     p->set< RCP<DataLayout> >("Node Scalar Data Layout", dl->node_scalar);
00295 
00296     ev = rcp(new PHAL::CahnHillWResid<EvalT,AlbanyTraits>(*p));
00297     fm0.template registerEvaluator<EvalT>(ev);
00298   }
00299 
00300   if (fieldManagerChoice == Albany::BUILD_RESID_FM)  {
00301     PHX::Tag<typename EvalT::ScalarT> res_tag("Scatter", dl->dummy);
00302     fm0.requireField<EvalT>(res_tag);
00303     return res_tag.clone();
00304   }
00305 
00306   else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
00307     Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
00308     return respUtils.constructResponses(fm0, *responseList, Teuchos::null, stateMgr);
00309   }
00310 
00311   return Teuchos::null;
00312 }
00313 
00314 
00315 #endif // ALBANY_CAHNHILLPROBLEM_HPP