Albany: a Trilinos-based PDE code: Albany_NavierStokes.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_NAVIERSTOKES_HPP
00008 #define ALBANY_NAVIERSTOKES_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 NavierStokes : public AbstractProblem {
00026   public:
00027   
00029     NavierStokes(const Teuchos::RCP<Teuchos::ParameterList>& params,
00030      const Teuchos::RCP<ParamLib>& paramLib,
00031      const int numDim_);
00032 
00034     ~NavierStokes();
00035 
00037     virtual int spatialDimension() const { return numDim; }
00038 
00040     virtual void buildProblem(
00041       Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> >  meshSpecs,
00042       StateManager& stateMgr);
00043 
00044     // Build evaluators
00045     virtual Teuchos::Array< Teuchos::RCP<const PHX::FieldTag> >
00046     buildEvaluators(
00047       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00048       const Albany::MeshSpecsStruct& meshSpecs,
00049       Albany::StateManager& stateMgr,
00050       Albany::FieldManagerChoice fmchoice,
00051       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00052 
00054     Teuchos::RCP<const Teuchos::ParameterList> getValidProblemParameters() const;
00055 
00056   private:
00057 
00059     NavierStokes(const NavierStokes&);
00060     
00062     NavierStokes& operator=(const NavierStokes&);
00063 
00064   public:
00065 
00067     template <typename EvalT> Teuchos::RCP<const PHX::FieldTag>
00068     constructEvaluators(
00069       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00070       const Albany::MeshSpecsStruct& meshSpecs,
00071       Albany::StateManager& stateMgr,
00072       Albany::FieldManagerChoice fmchoice,
00073       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00074     
00075     void constructDirichletEvaluators(const std::vector<std::string>& nodeSetIDs);
00076     void constructNeumannEvaluators(const Teuchos::RCP<Albany::MeshSpecsStruct>& meshSpecs); 
00077 
00078   protected:
00079 
00081     enum NS_VAR_TYPE {
00082       NS_VAR_TYPE_NONE,      
00083       NS_VAR_TYPE_CONSTANT,  
00084       NS_VAR_TYPE_DOF        
00085     };
00086 
00087     void getVariableType(Teuchos::ParameterList& paramList,
00088        const std::string& defaultType,
00089        NS_VAR_TYPE& variableType,
00090        bool& haveVariable,
00091        bool& haveEquation);
00092     std::string variableTypeToString(const NS_VAR_TYPE variableType);
00093 
00094   protected:
00095     
00096     bool periodic;     
00097     int numDim;        
00098 
00099     NS_VAR_TYPE flowType; 
00100     NS_VAR_TYPE heatType; 
00101     NS_VAR_TYPE neutType; 
00102 
00103     bool haveFlow;     
00104     bool haveHeat;     
00105     bool haveNeut;     
00106 
00107     bool haveFlowEq;     
00108     bool haveHeatEq;     
00109     bool haveNeutEq;     
00110 
00111     bool haveSource;   
00112     bool haveNeutSource;   
00113     bool havePSPG;     
00114     bool haveSUPG;     
00115     bool porousMedia;  
00116     
00117     Teuchos::RCP<Albany::Layouts> dl;
00118   };
00119 
00120 }
00121 
00122 #include <boost/type_traits/is_same.hpp>
00123 
00124 #include "Intrepid_FieldContainer.hpp"
00125 #include "Intrepid_DefaultCubatureFactory.hpp"
00126 #include "Shards_CellTopology.hpp"
00127 
00128 #include "Albany_Utils.hpp"
00129 #include "Albany_ProblemUtils.hpp"
00130 #include "Albany_EvaluatorUtils.hpp"
00131 #include "Albany_ResponseUtilities.hpp"
00132 
00133 #include "PHAL_NSContravarientMetricTensor.hpp"
00134 #include "PHAL_NSMaterialProperty.hpp"
00135 #include "PHAL_Source.hpp"
00136 #include "PHAL_NSBodyForce.hpp"
00137 #include "PHAL_NSPermeabilityTerm.hpp"
00138 #include "PHAL_NSForchheimerTerm.hpp"
00139 #include "PHAL_NSRm.hpp"
00140 #include "PHAL_NSTauM.hpp"
00141 #include "PHAL_NSTauT.hpp"
00142 #include "PHAL_NSMomentumResid.hpp"
00143 #include "PHAL_NSContinuityResid.hpp"
00144 #include "PHAL_NSThermalEqResid.hpp"
00145 #include "PHAL_NSNeutronEqResid.hpp"
00146 
00147 template <typename EvalT>
00148 Teuchos::RCP<const PHX::FieldTag>
00149 Albany::NavierStokes::constructEvaluators(
00150   PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00151   const Albany::MeshSpecsStruct& meshSpecs,
00152   Albany::StateManager& stateMgr,
00153   Albany::FieldManagerChoice fieldManagerChoice,
00154   const Teuchos::RCP<Teuchos::ParameterList>& responseList)
00155 {
00156   using Teuchos::RCP;
00157   using Teuchos::rcp;
00158   using Teuchos::ParameterList;
00159   using PHX::DataLayout;
00160   using PHX::MDALayout;
00161   using std::vector;
00162   using std::string;
00163   using std::map;
00164   using PHAL::AlbanyTraits;
00165  
00166   const CellTopologyData * const elem_top = &meshSpecs.ctd;
00167  
00168   RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > >
00169     intrepidBasis = Albany::getIntrepidBasis(meshSpecs.ctd);
00170   RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology (&meshSpecs.ctd));
00171   
00172   const int numNodes = intrepidBasis->getCardinality();
00173   const int worksetSize = meshSpecs.worksetSize;
00174   
00175   Intrepid::DefaultCubatureFactory<RealType> cubFactory;
00176   RCP <Intrepid::Cubature<RealType> > cubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree);
00177   
00178   const int numQPts = cubature->getNumPoints();
00179   const int numVertices = cellType->getNodeCount();
00180 
00181   // Print only for the residual specialization
00182   
00183   if(boost::is_same<EvalT, PHAL::AlbanyTraits::Residual>::value)
00184 
00185     *out << "Field Dimensions: Workset=" << worksetSize 
00186        << ", Vertices= " << numVertices
00187        << ", Nodes= " << numNodes
00188        << ", QuadPts= " << numQPts
00189        << ", Dim= " << numDim << std::endl;
00190   
00191    dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPts,numDim));
00192    TEUCHOS_TEST_FOR_EXCEPTION(dl->vectorAndGradientLayoutsAreEquivalent==false, std::logic_error,
00193                               "Data Layout Usage in NavierStokes problem assumes vecDim = numDim");
00194 
00195    Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
00196    bool supportsTransient=true;
00197    int offset=0;
00198 
00199    // Temporary variable used numerous times below
00200    Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
00201 
00202    // Define Field Names
00203 
00204    if (haveFlowEq) {
00205      Teuchos::ArrayRCP<string> dof_names(1);
00206      Teuchos::ArrayRCP<string> dof_names_dot(1);
00207      Teuchos::ArrayRCP<string> resid_names(1);
00208      dof_names[0] = "Velocity";
00209      dof_names_dot[0] = dof_names[0]+"_dot";
00210      resid_names[0] = "Momentum Residual";
00211      fm0.template registerEvaluator<EvalT>
00212        (evalUtils.constructGatherSolutionEvaluator(true, dof_names, dof_names_dot, offset));
00213 
00214      fm0.template registerEvaluator<EvalT>
00215        (evalUtils.constructDOFVecInterpolationEvaluator(dof_names[0], offset));
00216 
00217      fm0.template registerEvaluator<EvalT>
00218        (evalUtils.constructDOFVecInterpolationEvaluator(dof_names_dot[0], offset));
00219 
00220      fm0.template registerEvaluator<EvalT>
00221        (evalUtils.constructDOFVecGradInterpolationEvaluator(dof_names[0], offset));
00222 
00223      fm0.template registerEvaluator<EvalT>
00224        (evalUtils.constructScatterResidualEvaluator(true, resid_names,offset, "Scatter Momentum"));
00225      offset += numDim;
00226    }
00227    else if (haveFlow) { // Constant velocity
00228     RCP<ParameterList> p = rcp(new ParameterList);
00229 
00230     p->set<string>("Material Property Name", "Velocity");
00231     p->set< RCP<DataLayout> >("Data Layout", dl->qp_vector);
00232     p->set<string>("Coordinate Vector Name", "Coord Vec");
00233     p->set< RCP<DataLayout> >("Coordinate Vector Data Layout", dl->qp_vector);
00234 
00235     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00236     Teuchos::ParameterList& paramList = params->sublist("Flow");
00237     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00238 
00239     ev = rcp(new PHAL::NSMaterialProperty<EvalT,AlbanyTraits>(*p));
00240     fm0.template registerEvaluator<EvalT>(ev);
00241   }
00242 
00243    if (haveFlowEq) {
00244      Teuchos::ArrayRCP<string> dof_names(1);
00245      Teuchos::ArrayRCP<string> dof_names_dot(1);
00246      Teuchos::ArrayRCP<string> resid_names(1);
00247      dof_names[0] = "Pressure";
00248      dof_names_dot[0] = dof_names[0]+"_dot";
00249      resid_names[0] = "Continuity Residual";
00250      fm0.template registerEvaluator<EvalT>
00251        (evalUtils.constructGatherSolutionEvaluator(false, dof_names, dof_names_dot, offset));
00252 
00253      fm0.template registerEvaluator<EvalT>
00254        (evalUtils.constructDOFInterpolationEvaluator(dof_names[0], offset));
00255 
00256      fm0.template registerEvaluator<EvalT>
00257        (evalUtils.constructDOFInterpolationEvaluator(dof_names_dot[0], offset));
00258 
00259      fm0.template registerEvaluator<EvalT>
00260        (evalUtils.constructDOFGradInterpolationEvaluator(dof_names[0], offset));
00261 
00262      fm0.template registerEvaluator<EvalT>
00263        (evalUtils.constructScatterResidualEvaluator(false, resid_names,offset, "Scatter Continuity"));
00264      offset ++;
00265    }
00266 
00267    if (haveHeatEq) { // Gather Solution Temperature
00268      Teuchos::ArrayRCP<string> dof_names(1);
00269      Teuchos::ArrayRCP<string> dof_names_dot(1);
00270      Teuchos::ArrayRCP<string> resid_names(1);
00271      dof_names[0] = "Temperature";
00272      dof_names_dot[0] = dof_names[0]+"_dot";
00273      resid_names[0] = dof_names[0]+" Residual";
00274      fm0.template registerEvaluator<EvalT>
00275        (evalUtils.constructGatherSolutionEvaluator(false, dof_names, dof_names_dot, offset));
00276 
00277      fm0.template registerEvaluator<EvalT>
00278        (evalUtils.constructDOFInterpolationEvaluator(dof_names[0], offset));
00279 
00280      fm0.template registerEvaluator<EvalT>
00281        (evalUtils.constructDOFInterpolationEvaluator(dof_names_dot[0], offset));
00282 
00283      fm0.template registerEvaluator<EvalT>
00284        (evalUtils.constructDOFGradInterpolationEvaluator(dof_names[0], offset));
00285 
00286      fm0.template registerEvaluator<EvalT>
00287        (evalUtils.constructScatterResidualEvaluator(false, resid_names, offset, "Scatter Temperature"));
00288      offset ++;
00289    }
00290    else if (haveHeat) { // Constant temperature
00291     RCP<ParameterList> p = rcp(new ParameterList);
00292 
00293     p->set<string>("Material Property Name", "Temperature");
00294     p->set< RCP<DataLayout> >("Data Layout", dl->qp_scalar);
00295     p->set<string>("Coordinate Vector Name", "Coord Vec");
00296     p->set< RCP<DataLayout> >("Coordinate Vector Data Layout", dl->qp_vector);
00297 
00298     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00299     Teuchos::ParameterList& paramList = params->sublist("Heat");
00300     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00301 
00302     ev = rcp(new PHAL::NSMaterialProperty<EvalT,AlbanyTraits>(*p));
00303     fm0.template registerEvaluator<EvalT>(ev);
00304   }
00305 
00306    if (haveNeutEq) { // Gather Solution Neutron
00307      Teuchos::ArrayRCP<string> dof_names(1);
00308      Teuchos::ArrayRCP<string> dof_names_dot(1);
00309      Teuchos::ArrayRCP<string> resid_names(1);
00310      dof_names[0] = "Neutron Flux";
00311      dof_names_dot[0] = dof_names[0]+"_dot";
00312      resid_names[0] = dof_names[0]+" Residual";
00313      fm0.template registerEvaluator<EvalT>
00314        (evalUtils.constructGatherSolutionEvaluator(false, dof_names, dof_names_dot, offset));
00315 
00316      fm0.template registerEvaluator<EvalT>
00317        (evalUtils.constructDOFInterpolationEvaluator(dof_names[0], offset));
00318 
00319      fm0.template registerEvaluator<EvalT>
00320        (evalUtils.constructDOFInterpolationEvaluator(dof_names_dot[0], offset));
00321 
00322      fm0.template registerEvaluator<EvalT>
00323        (evalUtils.constructDOFGradInterpolationEvaluator(dof_names[0], offset));
00324 
00325      fm0.template registerEvaluator<EvalT>
00326        (evalUtils.constructScatterResidualEvaluator(false, resid_names, offset, "Scatter Neutron"));
00327      offset ++;
00328    }
00329    else if (haveNeut) { // Constant neutron flux
00330      RCP<ParameterList> p = rcp(new ParameterList);
00331 
00332      p->set<string>("Material Property Name", "Neutron Flux");
00333      p->set< RCP<DataLayout> >("Data Layout", dl->qp_scalar);
00334      p->set<string>("Coordinate Vector Name", "Coord Vec");
00335      p->set< RCP<DataLayout> >("Coordinate Vector Data Layout", dl->qp_vector);
00336      
00337      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00338      Teuchos::ParameterList& paramList = params->sublist("Neutronics");
00339      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00340      
00341      ev = rcp(new PHAL::NSMaterialProperty<EvalT,AlbanyTraits>(*p));
00342      fm0.template registerEvaluator<EvalT>(ev);
00343    }
00344 
00345    fm0.template registerEvaluator<EvalT>
00346      (evalUtils.constructGatherCoordinateVectorEvaluator());
00347 
00348    fm0.template registerEvaluator<EvalT>
00349      (evalUtils.constructMapToPhysicalFrameEvaluator(cellType, cubature));
00350 
00351    fm0.template registerEvaluator<EvalT>
00352      (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cubature));
00353 
00354   if (havePSPG || haveSUPG) { // Compute Contravarient Metric Tensor
00355     RCP<ParameterList> p = 
00356       rcp(new ParameterList("Contravarient Metric Tensor"));
00357 
00358     // Inputs: X, Y at nodes, Cubature, and Basis
00359     p->set<string>("Coordinate Vector Name","Coord Vec");
00360     p->set< RCP<DataLayout> >("Coordinate Data Layout", dl->vertices_vector);
00361     p->set< RCP<Intrepid::Cubature<RealType> > >("Cubature", cubature);
00362 
00363     p->set<RCP<shards::CellTopology> >("Cell Type", cellType);
00364 
00365     // Outputs: BF, weightBF, Grad BF, weighted-Grad BF, all in physical space
00366     p->set<string>("Contravarient Metric Tensor Name", "Gc");
00367     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00368 
00369     ev = rcp(new PHAL::NSContravarientMetricTensor<EvalT,AlbanyTraits>(*p));
00370     fm0.template registerEvaluator<EvalT>(ev);
00371   }
00372 
00373   if (haveFlowEq || haveHeatEq) { // Density
00374     RCP<ParameterList> p = rcp(new ParameterList);
00375 
00376     p->set<string>("Material Property Name", "Density");
00377     p->set< RCP<DataLayout> >("Data Layout", dl->qp_scalar);
00378     p->set<string>("Coordinate Vector Name", "Coord Vec");
00379     p->set< RCP<DataLayout> >("Coordinate Vector Data Layout", dl->qp_vector);
00380 
00381     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00382     Teuchos::ParameterList& paramList = params->sublist("Density");
00383     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00384 
00385     ev = rcp(new PHAL::NSMaterialProperty<EvalT,AlbanyTraits>(*p));
00386     fm0.template registerEvaluator<EvalT>(ev);
00387   }
00388 
00389   if (haveFlowEq) { // Viscosity
00390     RCP<ParameterList> p = rcp(new ParameterList);
00391 
00392     p->set<string>("Material Property Name", "Viscosity");
00393     p->set< RCP<DataLayout> >("Data Layout", dl->qp_scalar);
00394     p->set<string>("Coordinate Vector Name", "Coord Vec");
00395     p->set< RCP<DataLayout> >("Coordinate Vector Data Layout", dl->qp_vector);
00396 
00397     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00398     Teuchos::ParameterList& paramList = params->sublist("Viscosity");
00399     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00400 
00401     ev = rcp(new PHAL::NSMaterialProperty<EvalT,AlbanyTraits>(*p));
00402     fm0.template registerEvaluator<EvalT>(ev);
00403   }
00404 
00405   if (haveHeatEq) { // Specific Heat
00406     RCP<ParameterList> p = rcp(new ParameterList);
00407 
00408     p->set<string>("Material Property Name", "Specific Heat");
00409      p->set< RCP<DataLayout> >("Data Layout", dl->qp_scalar);
00410     p->set<string>("Coordinate Vector Name", "Coord Vec");
00411     p->set< RCP<DataLayout> >("Coordinate Vector Data Layout", dl->qp_vector);
00412 
00413     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00414     Teuchos::ParameterList& paramList = params->sublist("Specific Heat");
00415     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00416 
00417     ev = rcp(new PHAL::NSMaterialProperty<EvalT,AlbanyTraits>(*p));
00418     fm0.template registerEvaluator<EvalT>(ev);
00419   }
00420 
00421   if (haveHeatEq) { // Thermal conductivity
00422     RCP<ParameterList> p = rcp(new ParameterList);
00423 
00424     p->set<string>("Material Property Name", "Thermal Conductivity");
00425     p->set< RCP<DataLayout> >("Data Layout", dl->qp_scalar);
00426     p->set<string>("Coordinate Vector Name", "Coord Vec");
00427     p->set< RCP<DataLayout> >("Coordinate Vector Data Layout", dl->qp_vector);
00428 
00429     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00430     Teuchos::ParameterList& paramList = params->sublist("Thermal Conductivity");
00431     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00432 
00433     ev = rcp(new PHAL::NSMaterialProperty<EvalT,AlbanyTraits>(*p));
00434     fm0.template registerEvaluator<EvalT>(ev);
00435   }
00436 
00437   if (haveNeutEq) { // Neutron diffusion
00438     RCP<ParameterList> p = rcp(new ParameterList);
00439 
00440     p->set<string>("Material Property Name", "Neutron Diffusion Coefficient");
00441     p->set< RCP<DataLayout> >("Data Layout", dl->qp_scalar);
00442     p->set<string>("Coordinate Vector Name", "Coord Vec");
00443     p->set< RCP<DataLayout> >("Coordinate Vector Data Layout", dl->qp_vector);
00444     p->set<string>("Temperature Variable Name", "Temperature");
00445     p->set<string>("Absorption Cross Section Name", 
00446        "Absorption Cross Section");
00447     p->set<string>("Scattering Cross Section Name", 
00448        "Scattering Cross Section");
00449     p->set<string>("Average Scattering Angle Name", 
00450        "Average Scattering Angle");
00451 
00452     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00453     Teuchos::ParameterList& paramList = 
00454       params->sublist("Neutron Diffusion Coefficient");
00455     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00456 
00457     ev = rcp(new PHAL::NSMaterialProperty<EvalT,AlbanyTraits>(*p));
00458     fm0.template registerEvaluator<EvalT>(ev);
00459   }
00460 
00461   if (haveNeutEq) { // Reference temperature
00462     RCP<ParameterList> p = rcp(new ParameterList);
00463 
00464     p->set<string>("Material Property Name", "Reference Temperature");
00465     p->set< RCP<DataLayout> >("Data Layout", dl->qp_scalar);
00466     p->set<string>("Coordinate Vector Name", "Coord Vec");
00467     p->set< RCP<DataLayout> >("Coordinate Vector Data Layout", dl->qp_vector);
00468 
00469     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00470     Teuchos::ParameterList& paramList = 
00471       params->sublist("Reference Temperature");
00472     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00473 
00474     ev = rcp(new PHAL::NSMaterialProperty<EvalT,AlbanyTraits>(*p));
00475     fm0.template registerEvaluator<EvalT>(ev);
00476   }
00477 
00478   if (haveNeutEq) { // Neutron absorption cross section
00479     RCP<ParameterList> p = rcp(new ParameterList);
00480 
00481     p->set<string>("Material Property Name", "Absorption Cross Section");
00482     p->set< RCP<DataLayout> >("Data Layout", dl->qp_scalar);
00483     p->set<string>("Coordinate Vector Name", "Coord Vec");
00484     p->set< RCP<DataLayout> >("Coordinate Vector Data Layout", dl->qp_vector);
00485     p->set<string>("Temperature Variable Name", "Temperature");
00486 
00487     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00488     Teuchos::ParameterList& paramList = 
00489       params->sublist("Absorption Cross Section");
00490     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00491 
00492     ev = rcp(new PHAL::NSMaterialProperty<EvalT,AlbanyTraits>(*p));
00493     fm0.template registerEvaluator<EvalT>(ev);
00494   }
00495 
00496   if (haveNeutEq) { // Neutron scattering cross section
00497     RCP<ParameterList> p = rcp(new ParameterList);
00498 
00499     p->set<string>("Material Property Name", 
00500        "Scattering Cross Section");
00501     p->set< RCP<DataLayout> >("Data Layout", dl->qp_scalar);
00502     p->set<string>("Coordinate Vector Name", "Coord Vec");
00503     p->set< RCP<DataLayout> >("Coordinate Vector Data Layout", dl->qp_vector);
00504     p->set<string>("Temperature Variable Name", "Temperature");
00505 
00506     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00507     Teuchos::ParameterList& paramList = 
00508       params->sublist("Scattering Cross Section");
00509     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00510 
00511     ev = rcp(new PHAL::NSMaterialProperty<EvalT,AlbanyTraits>(*p));
00512     fm0.template registerEvaluator<EvalT>(ev);
00513   }
00514 
00515   if (haveNeutEq || (haveNeut && haveHeatEq)) { // Neutron fission cross section
00516     RCP<ParameterList> p = rcp(new ParameterList);
00517 
00518     p->set<string>("Material Property Name", "Fission Cross Section");
00519     p->set< RCP<DataLayout> >("Data Layout", dl->qp_scalar);
00520     p->set<string>("Coordinate Vector Name", "Coord Vec");
00521     p->set< RCP<DataLayout> >("Coordinate Vector Data Layout", dl->qp_vector);
00522     p->set<string>("Temperature Variable Name", "Temperature");
00523 
00524     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00525     Teuchos::ParameterList& paramList = 
00526       params->sublist("Fission Cross Section");
00527     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00528 
00529     ev = rcp(new PHAL::NSMaterialProperty<EvalT,AlbanyTraits>(*p));
00530     fm0.template registerEvaluator<EvalT>(ev);
00531   }
00532 
00533   if (haveNeutEq) { // Neutron released per fission
00534     RCP<ParameterList> p = rcp(new ParameterList);
00535 
00536     p->set<string>("Material Property Name", "Neutrons per Fission");
00537     p->set< RCP<DataLayout> >("Data Layout", dl->qp_scalar);
00538     p->set<string>("Coordinate Vector Name", "Coord Vec");
00539     p->set< RCP<DataLayout> >("Coordinate Vector Data Layout", dl->qp_vector);
00540     p->set<string>("Temperature Variable Name", "Temperature");
00541 
00542     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00543     Teuchos::ParameterList& paramList = 
00544       params->sublist("Neutrons per Fission");
00545     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00546 
00547     ev = rcp(new PHAL::NSMaterialProperty<EvalT,AlbanyTraits>(*p));
00548     fm0.template registerEvaluator<EvalT>(ev);
00549   }
00550 
00551   if (haveNeutEq) { // Average scattering angle
00552     RCP<ParameterList> p = rcp(new ParameterList);
00553 
00554     p->set<string>("Material Property Name", "Average Scattering Angle");
00555     p->set< RCP<DataLayout> >("Data Layout", dl->qp_scalar);
00556     p->set<string>("Coordinate Vector Name", "Coord Vec");
00557     p->set< RCP<DataLayout> >("Coordinate Vector Data Layout", dl->qp_vector);
00558     p->set<string>("Temperature Variable Name", "Temperature");
00559     p->set("Default Value", 0.0);
00560 
00561     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00562     Teuchos::ParameterList& paramList = 
00563       params->sublist("Average Scattering Angle");
00564     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00565 
00566     ev = rcp(new PHAL::NSMaterialProperty<EvalT,AlbanyTraits>(*p));
00567     fm0.template registerEvaluator<EvalT>(ev);
00568   }
00569 
00570   if (haveNeut && haveHeatEq) { // Energy released per fission
00571     RCP<ParameterList> p = rcp(new ParameterList);
00572 
00573     p->set<string>("Material Property Name", "Energy Released per Fission");
00574     p->set< RCP<DataLayout> >("Data Layout", dl->qp_scalar);
00575     p->set<string>("Coordinate Vector Name", "Coord Vec");
00576     p->set< RCP<DataLayout> >("Coordinate Vector Data Layout", dl->qp_vector);
00577 
00578     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00579     Teuchos::ParameterList& paramList = 
00580       params->sublist("Energy Released per Fission");
00581     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00582 
00583     ev = rcp(new PHAL::NSMaterialProperty<EvalT,AlbanyTraits>(*p));
00584     fm0.template registerEvaluator<EvalT>(ev);
00585   }
00586 
00587   if (haveFlowEq && haveHeat) { // Volumetric Expansion Coefficient
00588     RCP<ParameterList> p = rcp(new ParameterList);
00589 
00590     p->set<string>("Material Property Name", 
00591        "Volumetric Expansion Coefficient");
00592     p->set< RCP<DataLayout> >("Data Layout", dl->qp_scalar);
00593     p->set<string>("Coordinate Vector Name", "Coord Vec");
00594     p->set< RCP<DataLayout> >("Coordinate Vector Data Layout", dl->qp_vector);
00595 
00596     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00597     Teuchos::ParameterList& paramList = 
00598       params->sublist("Volumetric Expansion Coefficient");
00599     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00600 
00601     ev = rcp(new PHAL::NSMaterialProperty<EvalT,AlbanyTraits>(*p));
00602     fm0.template registerEvaluator<EvalT>(ev);
00603   }
00604 
00605   if (porousMedia) { // Porosity, \phi
00606     RCP<ParameterList> p = rcp(new ParameterList);
00607 
00608     p->set<string>("Material Property Name", 
00609        "Porosity");
00610     p->set< RCP<DataLayout> >("Data Layout", dl->qp_scalar);
00611     p->set<string>("Coordinate Vector Name", "Coord Vec");
00612     p->set< RCP<DataLayout> >("Coordinate Vector Data Layout", dl->qp_vector);
00613 
00614     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00615     Teuchos::ParameterList& paramList = 
00616       params->sublist("Porosity");
00617     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00618 
00619     ev = rcp(new PHAL::NSMaterialProperty<EvalT,AlbanyTraits>(*p));
00620     fm0.template registerEvaluator<EvalT>(ev);
00621   }
00622 
00623   if (porousMedia) { // Permeability Tensor, K
00624     RCP<ParameterList> p = rcp(new ParameterList);
00625 
00626     p->set<string>("Material Property Name", 
00627        "Permeability");
00628     p->set< RCP<DataLayout> >("Data Layout", dl->qp_scalar);
00629     p->set<string>("Coordinate Vector Name", "Coord Vec");
00630     p->set< RCP<DataLayout> >("Coordinate Vector Data Layout", dl->qp_vector);
00631 
00632     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00633     Teuchos::ParameterList& paramList = 
00634       params->sublist("Permeability");
00635     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00636 
00637     ev = rcp(new PHAL::NSMaterialProperty<EvalT,AlbanyTraits>(*p));
00638     fm0.template registerEvaluator<EvalT>(ev);
00639   }
00640 
00641   if (porousMedia) { // Forchheimer Tensor, F
00642     RCP<ParameterList> p = rcp(new ParameterList);
00643 
00644     p->set<string>("Material Property Name", 
00645        "Forchheimer");
00646     p->set< RCP<DataLayout> >("Data Layout", dl->qp_scalar);
00647     p->set<string>("Coordinate Vector Name", "Coord Vec");
00648     p->set< RCP<DataLayout> >("Coordinate Vector Data Layout", dl->qp_vector);
00649 
00650     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00651     Teuchos::ParameterList& paramList = 
00652       params->sublist("Forchheimer");
00653     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00654 
00655     ev = rcp(new PHAL::NSMaterialProperty<EvalT,AlbanyTraits>(*p));
00656     fm0.template registerEvaluator<EvalT>(ev);
00657   } 
00658 
00659   if (haveHeatEq && haveSource) { // Source
00660     RCP<ParameterList> p = rcp(new ParameterList);
00661 
00662     p->set<string>("Source Name", "Heat Source");
00663     p->set<string>("Variable Name", "Temperature");
00664     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00665     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00666     p->set<string>("QP Coordinate Vector Name", "Coord Vec");
00667     p->set<string>("Neutron Flux Name", "Neutron Flux");
00668     p->set<string>("Fission Cross Section Name", "Fission Cross Section");
00669     p->set<string>("Energy Released per Fission Name", 
00670        "Energy Released per Fission");
00671 
00672     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00673     Teuchos::ParameterList& paramList = params->sublist("Source Functions");
00674     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00675 
00676     ev = rcp(new PHAL::Source<EvalT,AlbanyTraits>(*p));
00677     fm0.template registerEvaluator<EvalT>(ev);
00678   }
00679 
00680   if (haveNeutEq && haveNeutSource) { // Source
00681     RCP<ParameterList> p = rcp(new ParameterList);
00682 
00683     p->set<string>("Source Name", "Neutron Source");
00684     p->set<string>("Variable Name", "Neutron Flux");
00685    p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00686     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00687     p->set<string>("QP Coordinate Vector Name", "Coord Vec");
00688 
00689     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00690     Teuchos::ParameterList& paramList = params->sublist("Neutron Source");
00691     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00692 
00693     ev = rcp(new PHAL::Source<EvalT,AlbanyTraits>(*p));
00694     fm0.template registerEvaluator<EvalT>(ev);
00695   }
00696 
00697   if (haveFlowEq) { // Body Force
00698     RCP<ParameterList> p = rcp(new ParameterList("Body Force"));
00699 
00700     //Input
00701     p->set<bool>("Have Heat", haveHeat);
00702     p->set<string>("Temperature QP Variable Name", "Temperature");
00703     p->set<string>("Density QP Variable Name", "Density");
00704     p->set<string>("Volumetric Expansion Coefficient QP Variable Name", 
00705        "Volumetric Expansion Coefficient");
00706 
00707     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00708     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector); 
00709 
00710     Teuchos::ParameterList& paramList = params->sublist("Body Force");
00711     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00712   
00713     //Output
00714     p->set<string>("Body Force Name", "Body Force");
00715 
00716     ev = rcp(new PHAL::NSBodyForce<EvalT,AlbanyTraits>(*p));
00717     fm0.template registerEvaluator<EvalT>(ev);
00718   }
00719 
00720 /*  if (haveFlowEq && haveNeumannx) { // Neumann BC's listed in input file and sidesets are present
00721 
00722     RCP<ParameterList> p = rcp(new ParameterList);
00723     // cell side
00724       
00725     const CellTopologyData * const side_top = elem_top->side[0].topology;
00726 
00727     p->set<string>("Side Set ID", meshSpecs.ssNames[0]);
00728     
00729     RCP<shards::CellTopology> sideType = rcp(new shards::CellTopology(side_top));
00730     RCP <Intrepid::Cubature<RealType> > sideCubature = cubFactory.create(*sideType, meshSpecs.cubatureDegree);
00731   
00732     // Inputs: X, Y at nodes, Cubature, and Basis
00733     p->set<string>("Node Variable Name", "Neumannx");
00734     p->set<string>("Coordinate Vector Name", "Coord Vec");
00735     p->set< RCP<DataLayout> >("Coordinate Data Layout", dl->vertices_vector);
00736     p->set< RCP<Intrepid::Cubature<RealType> > >("Cubature", cubature);
00737     
00738     p->set< RCP<Intrepid::Cubature<RealType> > >("Side Cubature", sideCubature);
00739     
00740     p->set< RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > > >
00741         ("Intrepid Basis", intrepidBasis);
00742 
00743     p->set<RCP<shards::CellTopology> >("Cell Type", cellType);
00744     p->set<RCP<shards::CellTopology> >("Side Type", sideType);
00745 
00746     // Output
00747     p->set< RCP<DataLayout> >("Node Scalar Data Layout", dl->node_scalar);
00748   
00749 //    p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00750     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00751 
00752     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00753     Teuchos::ParameterList& paramList = params->sublist("Neumann BCs");
00754     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00755     
00756     ev = rcp(new PHAL::Neumann<EvalT,AlbanyTraits>(*p));
00757     fm0.template registerEvaluator<EvalT>(ev);
00758   }
00759 
00760  if (haveFlowEq && haveNeumanny) { // Neumann BC's listed in input file and sidesets are present
00761 
00762     RCP<ParameterList> p = rcp(new ParameterList);
00763     // cell side
00764 
00765     const CellTopologyData * const side_top = elem_top->side[0].topology;
00766 
00767     p->set<string>("Side Set ID", meshSpecs.ssNames[0]);
00768 
00769     RCP<shards::CellTopology> sideType = rcp(new shards::CellTopology(side_top));
00770     RCP <Intrepid::Cubature<RealType> > sideCubature = cubFactory.create(*sideType, meshSpecs.cubatureDegree);
00771 
00772     // Inputs: X, Y at nodes, Cubature, and Basis
00773     p->set<string>("Node Variable Name", "Neumanny");
00774     p->set<string>("Coordinate Vector Name", "Coord Vec");
00775     p->set< RCP<DataLayout> >("Coordinate Data Layout", dl->vertices_vector);
00776     p->set< RCP<Intrepid::Cubature<RealType> > >("Cubature", cubature);
00777 
00778     p->set< RCP<Intrepid::Cubature<RealType> > >("Side Cubature", sideCubature);
00779 
00780     p->set< RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > > >
00781         ("Intrepid Basis", intrepidBasis);
00782 
00783     p->set<RCP<shards::CellTopology> >("Cell Type", cellType);
00784     p->set<RCP<shards::CellTopology> >("Side Type", sideType);
00785 
00786     // Output
00787     p->set< RCP<DataLayout> >("Node Scalar Data Layout", dl->node_scalar);
00788 
00789 //    p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00790     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00791 
00792     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00793     Teuchos::ParameterList& paramList = params->sublist("Neumann BCs");
00794     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00795 
00796     ev = rcp(new PHAL::Neumann<EvalT,AlbanyTraits>(*p));
00797     fm0.template registerEvaluator<EvalT>(ev);
00798   }
00799 
00800   if (haveFlowEq && haveNeumannz) { // Neumann BC's listed in input file and sidesets are present
00801 
00802     RCP<ParameterList> p = rcp(new ParameterList);
00803     // cell side
00804 
00805     const CellTopologyData * const side_top = elem_top->side[0].topology;
00806 
00807     p->set<string>("Side Set ID", meshSpecs.ssNames[0]);
00808 
00809     RCP<shards::CellTopology> sideType = rcp(new shards::CellTopology(side_top));
00810     RCP <Intrepid::Cubature<RealType> > sideCubature = cubFactory.create(*sideType, meshSpecs.cubatureDegree);
00811 
00812     // Inputs: X, Y at nodes, Cubature, and Basis
00813     p->set<string>("Node Variable Name", "Neumannz");
00814     p->set<string>("Coordinate Vector Name", "Coord Vec");
00815     p->set< RCP<DataLayout> >("Coordinate Data Layout", dl->vertices_vector);
00816     p->set< RCP<Intrepid::Cubature<RealType> > >("Cubature", cubature);
00817 
00818     p->set< RCP<Intrepid::Cubature<RealType> > >("Side Cubature", sideCubature);
00819 
00820     p->set< RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > > >
00821         ("Intrepid Basis", intrepidBasis);
00822 
00823     p->set<RCP<shards::CellTopology> >("Cell Type", cellType);
00824     p->set<RCP<shards::CellTopology> >("Side Type", sideType);
00825 
00826     // Output
00827     p->set< RCP<DataLayout> >("Node Scalar Data Layout", dl->node_scalar);
00828 
00829 //    p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00830     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00831 
00832     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00833     Teuchos::ParameterList& paramList = params->sublist("Neumann BCs");
00834     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00835 
00836     ev = rcp(new PHAL::Neumann<EvalT,AlbanyTraits>(*p));
00837     fm0.template registerEvaluator<EvalT>(ev);
00838    }
00839 */
00840 
00841    if (porousMedia) { // Permeability term in momentum equation
00842     std::cout <<"This is flow through porous media problem" << std::endl;
00843     RCP<ParameterList> p = rcp(new ParameterList("Permeability Term"));
00844 
00845     //Input
00846     p->set<string>("Velocity QP Variable Name", "Velocity");
00847     p->set<string>("Density QP Variable Name", "Density");
00848     p->set<string>("Viscosity QP Variable Name", "Viscosity");
00849     p->set<string>("Porosity QP Variable Name", "Porosity");
00850     p->set<string>("Permeability QP Variable Name", "Permeability");
00851 
00852     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00853     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00854     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00855     p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00856 
00857     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00858   
00859     //Output
00860     p->set<string>("Permeability Term", "Permeability Term");
00861 
00862     ev = rcp(new PHAL::NSPermeabilityTerm<EvalT,AlbanyTraits>(*p));
00863     fm0.template registerEvaluator<EvalT>(ev);
00864   } 
00865 
00866   if (porousMedia) { // Forchheimer term in momentum equation
00867     RCP<ParameterList> p = rcp(new ParameterList("Forchheimer Term"));
00868 
00869     //Input
00870     p->set<string>("Velocity QP Variable Name", "Velocity");
00871     p->set<string>("Density QP Variable Name", "Density");
00872     p->set<string>("Viscosity QP Variable Name", "Viscosity");
00873     p->set<string>("Porosity QP Variable Name", "Porosity");
00874     p->set<string>("Permeability QP Variable Name", "Permeability");
00875     p->set<string>("Forchheimer QP Variable Name", "Forchheimer");
00876 
00877     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00878     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00879     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00880     p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00881 
00882     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00883   
00884     //Output
00885     p->set<string>("Forchheimer Term", "Forchheimer Term");
00886 
00887     ev = rcp(new PHAL::NSForchheimerTerm<EvalT,AlbanyTraits>(*p));
00888     fm0.template registerEvaluator<EvalT>(ev);
00889   } 
00890 
00891   if (haveFlowEq) { // Rm
00892     RCP<ParameterList> p = rcp(new ParameterList("Rm"));
00893 
00894     //Input
00895     p->set<string>("Velocity QP Variable Name", "Velocity");
00896     p->set<string>("Velocity Dot QP Variable Name", "Velocity_dot");
00897     p->set<string>("Velocity Gradient QP Variable Name", "Velocity Gradient");
00898     p->set<string>("Pressure Gradient QP Variable Name", "Pressure Gradient");
00899     p->set<string>("Density QP Variable Name", "Density");
00900     p->set<string>("Body Force QP Variable Name", "Body Force");
00901     p->set<string>("Porosity QP Variable Name", "Porosity");
00902     p->set<string>("Permeability Term", "Permeability Term");
00903     p->set<string>("Forchheimer Term", "Forchheimer Term");
00904  
00905     p->set<bool>("Porous Media", porousMedia);
00906 
00907     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00908     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00909     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00910     p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00911 
00912     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00913   
00914     //Output
00915     p->set<string>("Rm Name", "Rm");
00916 
00917     ev = rcp(new PHAL::NSRm<EvalT,AlbanyTraits>(*p));
00918     fm0.template registerEvaluator<EvalT>(ev);
00919   }
00920 
00921   if (haveFlowEq && (haveSUPG || havePSPG)) { // Tau M
00922     RCP<ParameterList> p = rcp(new ParameterList("Tau M"));
00923 
00924     //Input
00925     p->set<string>("Velocity QP Variable Name", "Velocity");
00926     p->set<std::string>("Contravarient Metric Tensor Name", "Gc"); 
00927     p->set<string>("Density QP Variable Name", "Density");
00928     p->set<string>("Viscosity QP Variable Name", "Viscosity");
00929 
00930     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00931     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00932     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00933 
00934     //Output
00935     p->set<string>("Tau M Name", "Tau M");
00936 
00937     ev = rcp(new PHAL::NSTauM<EvalT,AlbanyTraits>(*p));
00938     fm0.template registerEvaluator<EvalT>(ev);
00939   }
00940 
00941   if (haveHeatEq && haveFlow && haveSUPG) { // Tau T
00942     RCP<ParameterList> p = rcp(new ParameterList("Tau T"));
00943 
00944     //Input
00945     p->set<string>("Velocity QP Variable Name", "Velocity");
00946     p->set<std::string>("Contravarient Metric Tensor Name", "Gc"); 
00947     p->set<string>("Thermal Conductivity Name", "Thermal Conductivity");
00948     p->set<string>("Density QP Variable Name", "Density");
00949     p->set<string>("Specific Heat QP Variable Name", "Specific Heat");
00950 
00951     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00952     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00953     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00954 
00955     //Output
00956     p->set<string>("Tau T Name", "Tau T");
00957 
00958     ev = rcp(new PHAL::NSTauT<EvalT,AlbanyTraits>(*p));
00959     fm0.template registerEvaluator<EvalT>(ev);
00960   }
00961 
00962   if (haveFlowEq) { // Momentum Resid
00963     RCP<ParameterList> p = rcp(new ParameterList("Momentum Resid"));
00964 
00965     //Input
00966     p->set<string>("Weighted BF Name", "wBF");
00967     p->set<string>("Weighted Gradient BF Name", "wGrad BF");
00968     p->set<string>("Velocity Gradient QP Variable Name", "Velocity Gradient");
00969     p->set<string>("Pressure QP Variable Name", "Pressure");
00970     p->set<string>("Pressure Gradient QP Variable Name", "Pressure Gradient");
00971     p->set<string>("Viscosity QP Variable Name", "Viscosity");
00972     p->set<string>("Rm Name", "Rm");
00973 
00974     p->set<bool>("Have SUPG", haveSUPG);
00975     p->set<string>("Velocity QP Variable Name", "Velocity");
00976     p->set<string>("Density QP Variable Name", "Density");
00977     p->set<string> ("Tau M Name", "Tau M");
00978 
00979     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00980     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);    
00981     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00982     p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
00983     p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00984     
00985     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00986   
00987     //Output
00988     p->set<string>("Residual Name", "Momentum Residual");
00989     p->set< RCP<DataLayout> >("Node Vector Data Layout", dl->node_vector);
00990 
00991     ev = rcp(new PHAL::NSMomentumResid<EvalT,AlbanyTraits>(*p));
00992     fm0.template registerEvaluator<EvalT>(ev);
00993   }
00994 
00995   if (haveFlowEq) { // Continuity Resid
00996     RCP<ParameterList> p = rcp(new ParameterList("Continuity Resid"));
00997 
00998     //Input
00999     p->set<string>("Weighted BF Name", "wBF");
01000     p->set<string>("Gradient QP Variable Name", "Velocity Gradient");
01001     p->set<string>("Density QP Variable Name", "Density");
01002 
01003     p->set<bool>("Have PSPG", havePSPG);
01004     p->set<string>("Weighted Gradient BF Name", "wGrad BF");
01005     p->set<std::string> ("Tau M Name", "Tau M");
01006     p->set<std::string> ("Rm Name", "Rm");
01007 
01008     p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
01009     p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
01010     p->set< RCP<PHX::DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
01011     p->set< RCP<PHX::DataLayout> >("QP Vector Data Layout", dl->qp_vector);
01012     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
01013 
01014     //Output
01015     p->set<string>("Residual Name", "Continuity Residual");
01016     p->set< RCP<DataLayout> >("Node Scalar Data Layout", dl->node_scalar);
01017 
01018     ev = rcp(new PHAL::NSContinuityResid<EvalT,AlbanyTraits>(*p));
01019     fm0.template registerEvaluator<EvalT>(ev);
01020   }
01021 
01022    if (haveHeatEq) { // Temperature Resid
01023     RCP<ParameterList> p = rcp(new ParameterList("Temperature Resid"));
01024 
01025     //Input
01026     p->set<string>("Weighted BF Name", "wBF");
01027     p->set<string>("Weighted Gradient BF Name", "wGrad BF");
01028     p->set<string>("QP Variable Name", "Temperature");
01029     p->set<string>("QP Time Derivative Variable Name", "Temperature_dot");
01030     p->set<string>("Gradient QP Variable Name", "Temperature Gradient");
01031     p->set<string>("Density QP Variable Name", "Density");
01032     p->set<string>("Specific Heat QP Variable Name", "Specific Heat");
01033     p->set<string>("Thermal Conductivity Name", "Thermal Conductivity");
01034     
01035     p->set<bool>("Have Source", haveSource);
01036     p->set<string>("Source Name", "Heat Source");
01037 
01038     p->set<bool>("Have Flow", haveFlow);
01039     p->set<string>("Velocity QP Variable Name", "Velocity");
01040     
01041     p->set<bool>("Have SUPG", haveSUPG);
01042     p->set<string> ("Tau T Name", "Tau T");
01043  
01044     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
01045     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
01046     p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
01047     p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
01048 
01049     //Output
01050     p->set<string>("Residual Name", "Temperature Residual");
01051     p->set< RCP<DataLayout> >("Node Scalar Data Layout", dl->node_scalar);
01052 
01053     ev = rcp(new PHAL::NSThermalEqResid<EvalT,AlbanyTraits>(*p));
01054     fm0.template registerEvaluator<EvalT>(ev);
01055   }
01056 
01057    if (haveNeutEq) { // Neutron Resid
01058     RCP<ParameterList> p = rcp(new ParameterList("Neutron Resid"));
01059 
01060     //Input
01061     p->set<string>("Weighted BF Name", "wBF");
01062     p->set<string>("Weighted Gradient BF Name", "wGrad BF");
01063     p->set<string>("QP Variable Name", "Neutron Flux");
01064     p->set<string>("Gradient QP Variable Name", "Neutron Flux Gradient");
01065     p->set<string>("Neutron Diffusion Name", "Neutron Diffusion Coefficient");
01066     p->set<string>("Neutron Absorption Name", "Absorption Cross Section");
01067     p->set<string>("Neutron Fission Name", "Fission Cross Section");
01068     p->set<string>("Neutrons per Fission Name", "Neutrons per Fission");
01069     
01070     p->set<bool>("Have Neutron Source", haveNeutSource);
01071     p->set<string>("Source Name", "Neutron Source");
01072  
01073     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
01074     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
01075     p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
01076     p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
01077 
01078     //Output
01079     p->set<string>("Residual Name", "Neutron Flux Residual");
01080     p->set< RCP<DataLayout> >("Node Scalar Data Layout", dl->node_scalar);
01081 
01082     ev = rcp(new PHAL::NSNeutronEqResid<EvalT,AlbanyTraits>(*p));
01083     fm0.template registerEvaluator<EvalT>(ev);
01084   }
01085 
01086   if (fieldManagerChoice == Albany::BUILD_RESID_FM)  {
01087     Teuchos::RCP<const PHX::FieldTag> ret_tag;
01088     if (haveFlowEq) {
01089       PHX::Tag<typename EvalT::ScalarT> mom_tag("Scatter Momentum", dl->dummy);
01090       fm0.requireField<EvalT>(mom_tag);
01091       PHX::Tag<typename EvalT::ScalarT> con_tag("Scatter Continuity", dl->dummy);
01092       fm0.requireField<EvalT>(con_tag);
01093       ret_tag = mom_tag.clone();
01094     }
01095     if (haveHeatEq) {
01096       PHX::Tag<typename EvalT::ScalarT> heat_tag("Scatter Temperature", dl->dummy);
01097       fm0.requireField<EvalT>(heat_tag);
01098       ret_tag = heat_tag.clone();
01099     }
01100     if (haveNeutEq) {
01101       PHX::Tag<typename EvalT::ScalarT> neut_tag("Scatter Neutron", dl->dummy);
01102       fm0.requireField<EvalT>(neut_tag);
01103       ret_tag = neut_tag.clone();
01104     }
01105     return ret_tag;
01106   }
01107   else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
01108     Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
01109     return respUtils.constructResponses(fm0, *responseList, Teuchos::null, stateMgr);
01110   }
01111 
01112   return Teuchos::null;
01113 }
01114 #endif // ALBANY_NAVIERSTOKES_HPP