00001
00002
00003
00004
00005
00006
00007 #include "Teuchos_TestForException.hpp"
00008 #include "Phalanx_DataLayout.hpp"
00009
00010 #include "Intrepid_FunctionSpaceTools.hpp"
00011
00012 namespace PHAL {
00013
00014
00015 template<typename EvalT, typename Traits>
00016 NSThermalEqResid<EvalT, Traits>::
00017 NSThermalEqResid(const Teuchos::ParameterList& p) :
00018 wBF (p.get<std::string> ("Weighted BF Name"),
00019 p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Scalar Data Layout") ),
00020 Temperature (p.get<std::string> ("QP Variable Name"),
00021 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00022 Tdot (p.get<std::string> ("QP Time Derivative Variable Name"),
00023 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00024 ThermalCond (p.get<std::string> ("Thermal Conductivity Name"),
00025 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00026 wGradBF (p.get<std::string> ("Weighted Gradient BF Name"),
00027 p.get<Teuchos::RCP<PHX::DataLayout> >("Node QP Vector Data Layout") ),
00028 TGrad (p.get<std::string> ("Gradient QP Variable Name"),
00029 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") ),
00030 rho (p.get<std::string> ("Density QP Variable Name"),
00031 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00032 Cp (p.get<std::string> ("Specific Heat QP Variable Name"),
00033 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") ),
00034
00035 TResidual (p.get<std::string> ("Residual Name"),
00036 p.get<Teuchos::RCP<PHX::DataLayout> >("Node Scalar Data Layout") ),
00037 haveSource (p.get<bool>("Have Source")),
00038 haveFlow (p.get<bool>("Have Flow")),
00039 haveSUPG (p.get<bool>("Have SUPG"))
00040 {
00041 if (p.isType<bool>("Disable Transient"))
00042 enableTransient = !p.get<bool>("Disable Transient");
00043 else enableTransient = true;
00044
00045 this->addDependentField(wBF);
00046 this->addDependentField(wGradBF);
00047 this->addDependentField(Temperature);
00048 this->addDependentField(TGrad);
00049 if (enableTransient) this->addDependentField(Tdot);
00050 this->addDependentField(ThermalCond);
00051 this->addDependentField(rho);
00052 this->addDependentField(Cp);
00053
00054 if (haveSource) {
00055 Source = PHX::MDField<ScalarT,Cell,QuadPoint>(
00056 p.get<std::string>("Source Name"),
00057 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") );
00058 this->addDependentField(Source);
00059 }
00060
00061 if (haveFlow) {
00062 V = PHX::MDField<ScalarT,Cell,QuadPoint,Dim>(
00063 p.get<std::string>("Velocity QP Variable Name"),
00064 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") );
00065 this->addDependentField(V);
00066 }
00067
00068 if (haveSUPG) {
00069 TauT = PHX::MDField<ScalarT,Cell,QuadPoint>(
00070 p.get<std::string>("Tau T Name"),
00071 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Scalar Data Layout") );
00072 this->addDependentField(TauT);
00073 }
00074
00075 this->addEvaluatedField(TResidual);
00076
00077 Teuchos::RCP<PHX::DataLayout> vector_dl =
00078 p.get< Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout");
00079 std::vector<PHX::DataLayout::size_type> dims;
00080 vector_dl->dimensions(dims);
00081 numQPs = dims[1];
00082 numDims = dims[2];
00083
00084 Teuchos::RCP<PHX::DataLayout> node_dl =
00085 p.get< Teuchos::RCP<PHX::DataLayout> >("Node Scalar Data Layout");
00086 std::vector<PHX::DataLayout::size_type> ndims;
00087 node_dl->dimensions(ndims);
00088 numNodes = ndims[1];
00089
00090
00091 flux.resize(dims[0], numQPs, numDims);
00092 convection.resize(dims[0], numQPs);
00093
00094 this->setName("NSThermalEqResid"+PHX::TypeString<EvalT>::value);
00095 }
00096
00097
00098 template<typename EvalT, typename Traits>
00099 void NSThermalEqResid<EvalT, Traits>::
00100 postRegistrationSetup(typename Traits::SetupData d,
00101 PHX::FieldManager<Traits>& fm)
00102 {
00103 this->utils.setFieldData(wBF,fm);
00104 this->utils.setFieldData(wGradBF,fm);
00105 this->utils.setFieldData(Temperature,fm);
00106 this->utils.setFieldData(TGrad,fm);
00107 if (enableTransient) this->utils.setFieldData(Tdot,fm);
00108 this->utils.setFieldData(ThermalCond,fm);
00109 this->utils.setFieldData(rho,fm);
00110 this->utils.setFieldData(Cp,fm);
00111 if (haveSource) this->utils.setFieldData(Source,fm);
00112 if (haveFlow) this->utils.setFieldData(V,fm);
00113 if (haveSUPG) this->utils.setFieldData(TauT,fm);
00114
00115 this->utils.setFieldData(TResidual,fm);
00116 }
00117
00118
00119 template<typename EvalT, typename Traits>
00120 void NSThermalEqResid<EvalT, Traits>::
00121 evaluateFields(typename Traits::EvalData workset)
00122 {
00123 typedef Intrepid::FunctionSpaceTools FST;
00124
00125 FST::scalarMultiplyDataData<ScalarT> (flux, ThermalCond, TGrad);
00126
00127 FST::integrate<ScalarT>(TResidual, flux, wGradBF, Intrepid::COMP_CPP, false);
00128
00129 for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00130 for (std::size_t qp=0; qp < numQPs; ++qp) {
00131 convection(cell,qp) = 0.0;
00132 if (haveSource) convection(cell,qp) -= Source(cell,qp);
00133 if (workset.transientTerms && enableTransient)
00134 convection(cell,qp) += rho(cell,qp) * Cp(cell,qp) * Tdot(cell,qp);
00135 if (haveFlow) {
00136 for (std::size_t i=0; i < numDims; ++i) {
00137 convection(cell,qp) +=
00138 rho(cell,qp) * Cp(cell,qp) * V(cell,qp,i) * TGrad(cell,qp,i);
00139 }
00140 }
00141 }
00142 }
00143
00144 FST::integrate<ScalarT>(TResidual, convection, wBF, Intrepid::COMP_CPP, true);
00145
00146 if (haveSUPG) {
00147 for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00148 for (std::size_t node=0; node < numNodes; ++node) {
00149 for (std::size_t qp=0; qp < numQPs; ++qp) {
00150 for (std::size_t j=0; j < numDims; ++j) {
00151 TResidual(cell,node) +=
00152 rho(cell,qp) * Cp(cell,qp) * TauT(cell,qp) * convection(cell,qp) *
00153 V(cell,qp,j) * wGradBF(cell,node,qp,j);
00154 }
00155 }
00156 }
00157 }
00158 }
00159
00160 }
00161
00162
00163 }
00164
