Albany: a Trilinos-based PDE code: Albany_NavierStokes.cpp 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 #include "Albany_NavierStokes.hpp"
00008 #include "AAdapt_InitialCondition.hpp"
00009 
00010 #include "Intrepid_FieldContainer.hpp"
00011 #include "Intrepid_DefaultCubatureFactory.hpp"
00012 #include "Shards_CellTopology.hpp"
00013 #include "PHAL_FactoryTraits.hpp"
00014 #include "Albany_Utils.hpp"
00015 #include "Albany_ProblemUtils.hpp"
00016 
00017 void
00018 Albany::NavierStokes::
00019 getVariableType(Teuchos::ParameterList& paramList,
00020     const std::string& defaultType,
00021     Albany::NavierStokes::NS_VAR_TYPE& variableType,
00022     bool& haveVariable,
00023     bool& haveEquation)
00024 {
00025   std::string type = paramList.get("Variable Type", defaultType);
00026   if (type == "None")
00027     variableType = NS_VAR_TYPE_NONE;
00028   else if (type == "Constant")
00029     variableType = NS_VAR_TYPE_CONSTANT;
00030   else if (type == "DOF")
00031     variableType = NS_VAR_TYPE_DOF;
00032   else
00033     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00034            "Unknown variable type " << type << std::endl);
00035   haveVariable = (variableType != NS_VAR_TYPE_NONE);
00036   haveEquation = (variableType == NS_VAR_TYPE_DOF);
00037 }
00038 
00039 std::string
00040 Albany::NavierStokes::
00041 variableTypeToString(Albany::NavierStokes::NS_VAR_TYPE variableType)
00042 {
00043   if (variableType == NS_VAR_TYPE_NONE)
00044     return "None";
00045   else if (variableType == NS_VAR_TYPE_CONSTANT)
00046     return "Constant";
00047   return "DOF";
00048 }
00049 
00050 Albany::NavierStokes::
00051 NavierStokes( const Teuchos::RCP<Teuchos::ParameterList>& params_,
00052              const Teuchos::RCP<ParamLib>& paramLib_,
00053              const int numDim_) :
00054   Albany::AbstractProblem(params_, paramLib_),
00055   haveFlow(false),
00056   haveHeat(false),
00057   haveNeut(false),
00058   haveFlowEq(false),
00059   haveHeatEq(false),
00060   haveNeutEq(false),
00061   haveSource(false),
00062   haveNeutSource(false),
00063   havePSPG(false),
00064   haveSUPG(false),
00065   porousMedia(false),
00066   numDim(numDim_)
00067 {
00068   if (numDim==1) periodic = params->get("Periodic BC", false);
00069   else           periodic = false;
00070   if (periodic) *out <<" Periodic Boundary Conditions being used." <<std::endl;
00071 
00072   getVariableType(params->sublist("Flow"), "DOF", flowType, 
00073       haveFlow, haveFlowEq);
00074   getVariableType(params->sublist("Heat"), "None", heatType, 
00075       haveHeat, haveHeatEq);
00076   getVariableType(params->sublist("Neutronics"), "None", neutType, 
00077       haveNeut, haveNeutEq);
00078 
00079   if (haveFlowEq) {
00080     havePSPG = params->get("Have Pressure Stabilization", true);
00081     porousMedia = params->get("Porous Media",false);
00082   }
00083 
00084   if (haveFlow && (haveFlowEq || haveHeatEq))
00085     haveSUPG = params->get("Have SUPG Stabilization", true);
00086 
00087   if (haveHeatEq)
00088     haveSource =  params->isSublist("Source Functions");
00089 
00090   if (haveNeutEq)
00091     haveNeutSource =  params->isSublist("Neutron Source");
00092 
00093   // Compute number of equations
00094   int num_eq = 0;
00095   if (haveFlowEq) num_eq += numDim+1;
00096   if (haveHeatEq) num_eq += 1;
00097   if (haveNeutEq) num_eq += 1;
00098   this->setNumEquations(num_eq);
00099 
00100   // Print out a summary of the problem
00101   *out << "Navier-Stokes problem:" << std::endl
00102        << "\tSpatial dimension:      " << numDim << std::endl
00103        << "\tFlow variables:         " << variableTypeToString(flowType) 
00104        << std::endl
00105        << "\tHeat variables:         " << variableTypeToString(heatType) 
00106        << std::endl
00107        << "\tNeutronics variables:   " << variableTypeToString(neutType) 
00108        << std::endl
00109        << "\tPressure stabilization: " << havePSPG << std::endl
00110        << "\tUpwind stabilization:   " << haveSUPG << std::endl
00111        << "\tPorous media:           " << porousMedia << std::endl;
00112 }
00113 
00114 Albany::NavierStokes::
00115 ~NavierStokes()
00116 {
00117 }
00118 
00119 void
00120 Albany::NavierStokes::
00121 buildProblem(
00122   Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> >  meshSpecs,
00123   Albany::StateManager& stateMgr)
00124 {
00125   using Teuchos::rcp;
00126 
00127  /* Construct All Phalanx Evaluators */
00128   TEUCHOS_TEST_FOR_EXCEPTION(meshSpecs.size()!=1,std::logic_error,"Problem supports one Material Block");
00129 
00130   fm.resize(1);
00131   fm[0]  = rcp(new PHX::FieldManager<PHAL::AlbanyTraits>);
00132   buildEvaluators(*fm[0], *meshSpecs[0], stateMgr, BUILD_RESID_FM, 
00133       Teuchos::null);
00134 
00135   if(meshSpecs[0]->nsNames.size() > 0) // Build a nodeset evaluator if nodesets are present
00136 
00137      constructDirichletEvaluators(meshSpecs[0]->nsNames);
00138 
00139   if(meshSpecs[0]->ssNames.size() > 0) // Build a sideset evaluator if sidesets are present
00140 
00141      constructNeumannEvaluators(meshSpecs[0]);
00142 
00143 }
00144 
00145 Teuchos::Array< Teuchos::RCP<const PHX::FieldTag> >
00146 Albany::NavierStokes::
00147 buildEvaluators(
00148   PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00149   const Albany::MeshSpecsStruct& meshSpecs,
00150   Albany::StateManager& stateMgr,
00151   Albany::FieldManagerChoice fmchoice,
00152   const Teuchos::RCP<Teuchos::ParameterList>& responseList)
00153 {
00154   // Call constructeEvaluators<EvalT>(*rfm[0], *meshSpecs[0], stateMgr);
00155   // for each EvalT in PHAL::AlbanyTraits::BEvalTypes
00156   ConstructEvaluatorsOp<NavierStokes> op(
00157     *this, fm0, meshSpecs, stateMgr, fmchoice, responseList);
00158   boost::mpl::for_each<PHAL::AlbanyTraits::BEvalTypes>(op);
00159   return *op.tags;
00160 }
00161 
00162 void
00163 Albany::NavierStokes::constructDirichletEvaluators(
00164         const std::vector<std::string>& nodeSetIDs)
00165 {
00166    // Construct Dirichlet evaluators for all nodesets and names
00167    std::vector<std::string> dirichletNames(neq);
00168    int index = 0;
00169    if (haveFlowEq) {
00170      dirichletNames[index++] = "ux";
00171      if (numDim>=2) dirichletNames[index++] = "uy";
00172      if (numDim==3) dirichletNames[index++] = "uz";
00173      dirichletNames[index++] = "p";
00174    }
00175    if (haveHeatEq) dirichletNames[index++] = "T";
00176    if (haveNeutEq) dirichletNames[index++] = "phi";
00177    Albany::BCUtils<Albany::DirichletTraits> dirUtils;
00178    dfm = dirUtils.constructBCEvaluators(nodeSetIDs, dirichletNames,
00179                                           this->params, this->paramLib);
00180 }
00181 
00182 // Neumann BCs
00183 void
00184 Albany::NavierStokes::constructNeumannEvaluators(const Teuchos::RCP<Albany::MeshSpecsStruct>& meshSpecs)
00185 {
00186 
00187    // Note: we only enter this function if sidesets are defined in the mesh file
00188    // i.e. meshSpecs.ssNames.size() > 0
00189 
00190    Albany::BCUtils<Albany::NeumannTraits> nbcUtils;
00191 
00192    // Check to make sure that Neumann BCs are given in the input file
00193 
00194    if(!nbcUtils.haveBCSpecified(this->params)) {
00195       return;
00196    }
00197 
00198 
00199    // Construct BC evaluators for all side sets and names
00200    // Note that the string index sets up the equation offset, so ordering is important
00201 
00202    // Currently we aren't exactly doing this right.  I think to do this
00203    // correctly we need different neumann evaluators for each DOF (velocity,
00204    // pressure, temperature, flux) since velocity is a vector and the 
00205    // others are scalars.  The dof_names stuff is only used
00206    // for robin conditions, so at this point, as long as we don't enable
00207    // robin conditions, this should work.
00208 
00209    std::vector<std::string> nbcNames;
00210    Teuchos::RCP< Teuchos::Array<std::string> > dof_names = 
00211      Teuchos::rcp(new Teuchos::Array<std::string>);
00212    Teuchos::Array<Teuchos::Array<int> > offsets;
00213    int idx = 0;
00214    if (haveFlowEq) {
00215      nbcNames.push_back("ux");
00216      offsets.push_back(Teuchos::Array<int>(1,idx++));
00217      if (numDim>=2) {
00218        nbcNames.push_back("uy");
00219        offsets.push_back(Teuchos::Array<int>(1,idx++));
00220      }
00221      if (numDim==3) {
00222        nbcNames.push_back("uz");
00223        offsets.push_back(Teuchos::Array<int>(1,idx++));
00224      }
00225      nbcNames.push_back("p");
00226      offsets.push_back(Teuchos::Array<int>(1,idx++));
00227      dof_names->push_back("Velocity");
00228      dof_names->push_back("Pressure");
00229    }
00230    if (haveHeatEq) {
00231      nbcNames.push_back("T");
00232      offsets.push_back(Teuchos::Array<int>(1,idx++));
00233      dof_names->push_back("Temperature");
00234    }
00235    if (haveNeutEq) {
00236      nbcNames.push_back("phi");
00237      offsets.push_back(Teuchos::Array<int>(1,idx++));
00238      dof_names->push_back("Neutron Flux");
00239    }
00240    
00241    // Construct BC evaluators for all possible names of conditions
00242    // Should only specify flux vector components (dudx, dudy, dudz), or dudn, not both
00243    std::vector<std::string> condNames(2); //dudx, dudy, dudz, dudn
00244 
00245    // Note that sidesets are only supported for two and 3D currently
00246    if(numDim == 2) 
00247     condNames[0] = "(dudx, dudy)";
00248    else if(numDim == 3)
00249     condNames[0] = "(dudx, dudy, dudz)";
00250    else
00251     TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00252        std::endl << "Error: Sidesets only supported in 2 and 3D." << std::endl);
00253 
00254    condNames[1] = "dudn";
00255 
00256    nfm.resize(1);
00257 
00258    nfm[0] = nbcUtils.constructBCEvaluators(meshSpecs, nbcNames, 
00259              Teuchos::arcp(dof_names), 
00260              false, 0, condNames, offsets, dl,
00261              this->params, this->paramLib);
00262 }
00263 
00264 
00265 Teuchos::RCP<const Teuchos::ParameterList>
00266 Albany::NavierStokes::getValidProblemParameters() const
00267 {
00268   Teuchos::RCP<Teuchos::ParameterList> validPL =
00269     this->getGenericProblemParams("ValidNavierStokesParams");
00270 
00271   if (numDim==1)
00272     validPL->set<bool>("Periodic BC", false, "Flag to indicate periodic BC for 1D problems");
00273   validPL->set<bool>("Have Pressure Stabilization", true);
00274   validPL->set<bool>("Have SUPG Stabilization", true);
00275   validPL->set<bool>("Porous Media", false, "Flag to use porous media equations");
00276   validPL->sublist("Flow", false, "");
00277   validPL->sublist("Heat", false, "");
00278   validPL->sublist("Neutronics", false, "");
00279   validPL->sublist("Thermal Conductivity", false, "");
00280   validPL->sublist("Density", false, "");
00281   validPL->sublist("Viscosity", false, "");
00282   validPL->sublist("Volumetric Expansion Coefficient", false, "");
00283   validPL->sublist("Specific Heat", false, "");
00284   validPL->sublist("Body Force", false, "");
00285   validPL->sublist("Porosity", false, "");
00286   validPL->sublist("Permeability", false, "");
00287   validPL->sublist("Forchheimer", false, "");
00288   
00289   validPL->sublist("Neutron Source", false, "");
00290   validPL->sublist("Neutron Diffusion Coefficient", false, "");
00291   validPL->sublist("Absorption Cross Section", false, "");
00292   validPL->sublist("Fission Cross Section", false, "");
00293   validPL->sublist("Neutrons per Fission", false, "");
00294   validPL->sublist("Scattering Cross Section", false, "");
00295   validPL->sublist("Average Scattering Angle", false, "");
00296   validPL->sublist("Energy Released per Fission", false, "");
00297 
00298   return validPL;
00299 }
00300