Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007 #include "Teuchos_TestForException.hpp"
00008 #include "Phalanx_DataLayout.hpp"
00009 #include "Sacado.hpp"
00010
00011 #include "Intrepid_FunctionSpaceTools.hpp"
00012
00013 namespace PHAL {
00014 const double pi = 3.1415926535897932385;
00015
00016
00017 template<typename EvalT, typename Traits>
00018 ComprNSBodyForce<EvalT, Traits>::
00019 ComprNSBodyForce(const Teuchos::ParameterList& p) :
00020 force(p.get<std::string>("Body Force Name"),
00021 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout") )
00022 {
00023 std::cout << "Compr NS body force constructor!" << std::endl;
00024 Teuchos::ParameterList* bf_list =
00025 p.get<Teuchos::ParameterList*>("Parameter List");
00026
00027 std::string type = bf_list->get("Type", "None");
00028 if (type == "None") {
00029 bf_type = NONE;
00030 }
00031 else if (type == "Taylor-Green Vortex") {
00032 std::cout << "Taylor-Green Vortex source" << std::endl;
00033 bf_type = TAYLOR_GREEN_VORTEX;
00034 coordVec = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00035 p.get<std::string>("Coordinate Vector Name"),
00036 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Gradient Data Layout") );
00037 this->addDependentField(coordVec);
00038 }
00039
00040 this->addEvaluatedField(force);
00041
00042 Teuchos::RCP<PHX::DataLayout> gradient_dl =
00043 p.get< Teuchos::RCP<PHX::DataLayout> >("QP Gradient Data Layout");
00044 std::vector<PHX::DataLayout::size_type> dims;
00045 gradient_dl->dimensions(dims);
00046 numQPs = dims[1];
00047 numDims = dims[2];
00048 Teuchos::RCP<PHX::DataLayout> vector_dl =
00049 p.get< Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout");
00050 vector_dl->dimensions(dims);
00051 vecDim = dims[2];
00052
00053 std::cout << " in Compr NS Stokes source! " << std::endl;
00054 std::cout << " vecDim = " << vecDim << std::endl;
00055 std::cout << " numDims = " << numDims << std::endl;
00056 std::cout << " numQPs = " << numQPs << std::endl;
00057
00058
00059 this->setName("ComprNSBodyForce"+PHX::TypeString<EvalT>::value);
00060 }
00061
00062
00063 template<typename EvalT, typename Traits>
00064 void ComprNSBodyForce<EvalT, Traits>::
00065 postRegistrationSetup(typename Traits::SetupData d,
00066 PHX::FieldManager<Traits>& fm)
00067 {
00068 if (bf_type == TAYLOR_GREEN_VORTEX) {
00069 this->utils.setFieldData(coordVec,fm);
00070 }
00071 this->utils.setFieldData(force,fm);
00072 }
00073
00074
00075 template<typename EvalT, typename Traits>
00076 void ComprNSBodyForce<EvalT, Traits>::
00077 evaluateFields(typename Traits::EvalData workset)
00078 {
00079 if (bf_type == NONE) {
00080 for (std::size_t cell=0; cell < workset.numCells; ++cell)
00081 for (std::size_t qp=0; qp < numQPs; ++qp)
00082 for (std::size_t i=0; i < vecDim; ++i)
00083 force(cell,qp,i) = 0.0;
00084 }
00085 else if (bf_type == TAYLOR_GREEN_VORTEX) {
00086 const RealType time = workset.current_time;
00087 const double Re = 1.0;
00088 const double Pr = 0.72;
00089 const double gamma_gas = 1.4;
00090 const double kappa = 1.0;
00091 const double mu = 1.0/Re;
00092 const double Rgas = 0.714285733;
00093 for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00094 for (std::size_t qp=0; qp < numQPs; ++qp) {
00095 ScalarT* f = &force(cell,qp,0);
00096 MeshScalarT x2pi = 2.0*pi*coordVec(cell,qp,0);
00097 MeshScalarT y2pi = 2.0*pi*coordVec(cell,qp,1);
00098 f[0] = 0.0;
00099 f[1] = 2.0*exp(-2.0*time)*cos(x2pi)*(-sin(y2pi) + exp(-2.0*time)*sin(x2pi)*pi + 4.0*mu*pi*pi*sin(y2pi)) + 2.0*Rgas*pi*exp(-4.0*time)*sin(x2pi);
00100 f[2] = 2.0*exp(-2.0*time)*cos(y2pi)*(sin(x2pi) + exp(-2.0*time)*sin(y2pi)*pi - 4.0*mu*pi*pi*sin(x2pi)) + 2.0*Rgas*pi*exp(-4.0*time)*sin(y2pi);
00101 f[3] = -2.0*exp(-4.0*time)*(-2.0*cos(x2pi) - 2.0*cos(y2pi) + exp(-2.0*time)*cos(x2pi)*sin(y2pi)*sin(x2pi)*pi
00102 - exp(-2.0*time)*sin(x2pi)*cos(y2pi)*sin(y2pi)*pi)
00103 + (gamma_gas-1.0)/Rgas*4.0*mu*exp(-2.0*time)*sin(x2pi)*sin(y2pi)*pi*2.0*exp(-2.0*time)*sin(x2pi)*pi*sin(y2pi)
00104 + (gamma_gas-1.0)/Rgas*4.0*mu*exp(-2.0*time)*sin(x2pi)*pi*sin(y2pi)*2.0*exp(-2.0*time)*sin(x2pi)*pi*sin(y2pi)
00105 - gamma_gas*kappa/(Pr*Re)*4.0*exp(-4.0*time)*pi*pi*(cos(x2pi) + cos(y2pi));
00106 }
00107 }
00108 }
00109 }
00110
00111 }
00112
