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 LinComprNSBodyForce<EvalT, Traits>::
00019 LinComprNSBodyForce(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 << "Lin 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 == "Steady Euler") {
00032 bf_type = STEADYEUL;
00033 coordVec = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00034 p.get<std::string>("Coordinate Vector Name"),
00035 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Gradient Data Layout") );
00036 this->addDependentField(coordVec);
00037 }
00038 else if (type == "Unsteady Euler MMS") {
00039 bf_type = UNSTEADYEULMMS;
00040 coordVec = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00041 p.get<std::string>("Coordinate Vector Name"),
00042 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Gradient Data Layout") );
00043 this->addDependentField(coordVec);
00044 }
00045 else if (type == "Driven Pulse") {
00046 bf_type = DRIVENPULSE;
00047 coordVec = PHX::MDField<MeshScalarT,Cell,QuadPoint,Dim>(
00048 p.get<std::string>("Coordinate Vector Name"),
00049 p.get<Teuchos::RCP<PHX::DataLayout> >("QP Gradient Data Layout") );
00050 this->addDependentField(coordVec);
00051 }
00052
00053 this->addEvaluatedField(force);
00054
00055 Teuchos::RCP<PHX::DataLayout> gradient_dl =
00056 p.get< Teuchos::RCP<PHX::DataLayout> >("QP Gradient Data Layout");
00057 std::vector<PHX::DataLayout::size_type> dims;
00058 gradient_dl->dimensions(dims);
00059 numQPs = dims[1];
00060 numDims = dims[2];
00061 Teuchos::RCP<PHX::DataLayout> vector_dl =
00062 p.get< Teuchos::RCP<PHX::DataLayout> >("QP Vector Data Layout");
00063 vector_dl->dimensions(dims);
00064 vecDim = dims[2];
00065
00066 std::cout << " in Lin Compr NS Stokes source! " << std::endl;
00067 std::cout << " vecDim = " << vecDim << std::endl;
00068 std::cout << " numDims = " << numDims << std::endl;
00069 std::cout << " numQPs = " << numQPs << std::endl;
00070
00071
00072 this->setName("LinComprNSBodyForce"+PHX::TypeString<EvalT>::value);
00073 }
00074
00075
00076 template<typename EvalT, typename Traits>
00077 void LinComprNSBodyForce<EvalT, Traits>::
00078 postRegistrationSetup(typename Traits::SetupData d,
00079 PHX::FieldManager<Traits>& fm)
00080 {
00081 if (bf_type == STEADYEUL || bf_type == UNSTEADYEULMMS || bf_type == DRIVENPULSE) {
00082 this->utils.setFieldData(coordVec,fm);
00083 }
00084
00085 this->utils.setFieldData(force,fm);
00086 }
00087
00088
00089 template<typename EvalT, typename Traits>
00090 void LinComprNSBodyForce<EvalT, Traits>::
00091 evaluateFields(typename Traits::EvalData workset)
00092 {
00093 if (bf_type == NONE) {
00094 for (std::size_t cell=0; cell < workset.numCells; ++cell)
00095 for (std::size_t qp=0; qp < numQPs; ++qp)
00096 for (std::size_t i=0; i < vecDim; ++i)
00097 force(cell,qp,i) = 0.0;
00098 }
00099 else if (bf_type == STEADYEUL) {
00100 const double ubar = 1.0;
00101 const double vbar = 1.0;
00102 const double zetabar = 1.0;
00103 const double pbar = 0.0;
00104 const double gamma_gas = 1.4;
00105 for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00106 for (std::size_t qp=0; qp < numQPs; ++qp) {
00107 ScalarT* f = &force(cell,qp,0);
00108 MeshScalarT x2pi = 2.0*pi*coordVec(cell,qp,0);
00109 MeshScalarT y2pi = 2.0*pi*coordVec(cell,qp,1);
00110 MeshScalarT x = coordVec(cell,qp,0);
00111 MeshScalarT y = coordVec(cell,qp,1);
00112 f[0] = -1.0*(ubar*(y - x*sin(x)) + vbar*x + zetabar*2.0*x*(0.5-y));
00113 f[1] = -1.0*(ubar*cos(x)*y + vbar*sin(x) - zetabar*x*x);
00114 f[2] = -1.0*(gamma_gas*pbar*(y - x*sin(x) + sin(x)) + ubar*2.0*x*(0.5-y) - vbar*x*x);
00115 }
00116 }
00117 }
00118 else if (bf_type == UNSTEADYEULMMS) {
00119 const double ubar = 0.0;
00120 const double vbar = 0.0;
00121 const double zetabar = 1.0;
00122 const double pbar = 0.7142857;
00123 const double a = 1.0;
00124 const double gamma_gas = 1.4;
00125 const RealType time = workset.current_time;
00126 for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00127 for (std::size_t qp=0; qp < numQPs; ++qp) {
00128 ScalarT* f = &force(cell,qp,0);
00129 MeshScalarT x2pi = 2.0*pi*coordVec(cell,qp,0);
00130 MeshScalarT y2pi = 2.0*pi*coordVec(cell,qp,1);
00131 MeshScalarT x = coordVec(cell,qp,0);
00132 MeshScalarT y = coordVec(cell,qp,1);
00133 f[0] = -1.0*exp(-a*time)*(-a*sin(x2pi)*cos(y2pi) + ubar*2.0*pi*cos(x2pi)*cos(y2pi)
00134 -vbar*2.0*pi*sin(x2pi)*sin(y2pi) + 2.0*pi*zetabar*cos(x2pi)*sin(y2pi));
00135 f[1] = -1.0*exp(-a*time)*(-a*cos(x2pi)*sin(y2pi) - 2.0*pi*ubar*sin(x2pi)*sin(y2pi)
00136 + vbar*2.0*pi*cos(x2pi)*cos(y2pi) + 2.0*pi*zetabar*sin(x2pi)*cos(y2pi));
00137 f[2] = -1.0*exp(-a*time)*(-a*sin(x2pi)*sin(y2pi) + gamma_gas*pbar*4.0*pi*cos(x2pi)*cos(y2pi) +
00138 ubar*2.0*pi*cos(x2pi)*sin(y2pi) + vbar*2.0*pi*sin(x2pi)*cos(y2pi));
00139 }
00140 }
00141 }
00142 else if (bf_type == DRIVENPULSE) {
00143 const RealType time = workset.current_time;
00144 const double tref = 1.0/347.9693;
00145 for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00146 for (std::size_t qp=0; qp < numQPs; ++qp) {
00147 ScalarT* f = &force(cell,qp,0);
00148 MeshScalarT x = coordVec(cell,qp,0);
00149 MeshScalarT y = coordVec(cell,qp,1);
00150 f[0] = 0.0;
00151 if (x >= 0.9 & x <= 1.0 & y >= 0.9 & y <= 1.0)
00152 f[1] = (1.0e-4)*cos(2.0*pi*1000*time*tref);
00153 else
00154 f[1] = 0.0;
00155 f[2] = 0.0;
00156 }
00157 }
00158
00159 }
00160 }
00161
00162 }
00163
