52 "Computes wall shear stress field.");
67 int nfields =
m_f->m_variables.size();
68 int expdim =
m_f->m_graph->GetSpaceDimension();
71 if (
m_f->m_exp[0]->GetNumElmts() == 0)
78 ASSERTL0(
false,
"Error: wss for a 1D problem cannot "
98 m_f->m_exp.resize(nshear);
99 for (i = nfields; i < nshear; ++i)
101 m_f->m_exp[nfields + i] =
102 m_f->AppendExpList(
m_f->m_numHomogeneousDir);
108 m_f->m_exp[0]->GetGraph());
111 map<int, int> BndRegionMap;
113 for (
auto &breg_it : bregions)
115 BndRegionMap[breg_it.first] = cnt++;
119 for (
int b = 0; b <
m_f->m_bndRegionsToWrite.size(); ++b)
121 if (BndRegionMap.count(
m_f->m_bndRegionsToWrite[b]) == 1)
123 int bnd = BndRegionMap[
m_f->m_bndRegionsToWrite[b]];
126 for (i = 0; i < nshear; i++)
128 BndExp[i] =
m_f->m_exp[i]->UpdateBndCondExpansion(bnd);
130 for (i = 0; i < nfields; i++)
132 m_f->m_exp[i]->GetBndElmtExpansion(bnd, BndElmtExp[i]);
136 int nqb = BndExp[0]->GetTotPoints();
137 int nqe = BndElmtExp[0]->GetTotPoints();
142 for (i = 0; i < ngrad; ++i)
147 for (i = 0; i < nstress; ++i)
155 for (i = 0; i < nstress; ++i)
160 for (i = 0; i < nshear; ++i)
178 BndElmtExp[i]->PhysDeriv(velocity[i],
184 BndElmtExp[i]->PhysDeriv(
225 for (j = 0; j < nstress; ++j)
227 m_f->m_exp[0]->ExtractElmtToBndPhys(bnd, stress[j], fstress[j]);
232 m_f->m_exp[0]->GetBoundaryNormals(bnd, normals);
255 fshear[nshear - 1], 1, fshear[nshear - 1], 1);
257 Vmath::Smul(nqb, -1.0, fshear[nshear - 1], 1, fshear[nshear - 1],
263 fshear[i], 1, fshear[i], 1);
264 BndExp[i]->FwdTransLocalElmt(fshear[i],
265 BndExp[i]->UpdateCoeffs());
273 fshear[nshear - 1], 1, fshear[nshear - 1], 1);
275 Vmath::Vsqrt(nqb, fshear[nshear - 1], 1, fshear[nshear - 1], 1);
276 BndExp[nshear - 1]->FwdTransLocalElmt(
277 fshear[nshear - 1], BndExp[nshear - 1]->UpdateCoeffs());
283 m_f->m_variables[0] =
"Shear_s";
284 m_f->m_variables[1] =
"Shear_n";
285 m_f->m_variables[2] =
"Shear_mag";
289 m_f->m_variables[0] =
"Shear_s";
290 m_f->m_variables[1] =
"Shear_n";
291 m_f->m_variables[2] =
"Shear_z";
292 m_f->m_variables[3] =
"Shear_mag";
301 int npoints = exp[0]->GetNpoints();
303 if (boost::iequals(
m_f->m_variables[0],
"u"))
306 m_mu =
m_f->m_session->GetParameter(
"Kinvis");
310 else if (boost::iequals(
m_f->m_variables[0],
"rho") &&
311 boost::iequals(
m_f->m_variables[1],
"rhou"))
314 std::string m_ViscosityType;
315 m_f->m_session->LoadParameter(
"mu",
m_mu, 1.78e-05);
316 m_f->m_session->LoadParameter(
"lambda", lambda, -2.0 / 3.0);
317 m_f->m_session->LoadSolverInfo(
"ViscosityType", m_ViscosityType,
320 if (m_ViscosityType ==
"Variable")
325 m_f->m_session->LoadSolverInfo(
"EquationOfState", eosType,
327 m_idealGas = boost::iequals(eosType,
"IdealGas");
330 "Only IdealGas EOS implemented for Variable ViscosityType");
340 m_f->m_session->LoadParameter(
"Gamma", m_gamma, 1.4);
341 m_f->m_session->LoadParameter(
"pInf", m_pInf, 101325);
342 m_f->m_session->LoadParameter(
"rhoInf",
m_rhoInf, 1.225);
343 m_f->m_session->LoadParameter(
"GasConstant", m_gasConstant,
345 m_f->m_session->LoadParameter(
"Tref", m_Tref, 288.15);
346 m_TRatioSutherland = 110.0 / m_Tref;
349 cv_inv = (m_gamma - 1.0) / m_gasConstant;
359 exp[i + 1]->GetPhys(), 1, tmp, 1);
361 Vmath::Vdiv(npoints, tmp, 1, exp[0]->GetPhys(), 1, tmp, 1);
365 Vmath::Vsub(npoints, energy, 1, tmp, 1, energy, 1);
368 Vmath::Vdiv(npoints, energy, 1, exp[0]->GetPhys(), 1, energy, 1);
370 Vmath::Smul(npoints, cv_inv, energy, 1, temperature, 1);
382 for (
int i = 0; i < npoints; ++i)
384 ratio = temperature[i] / T_star;
385 mu[i] = mu_star * ratio *
sqrt(ratio) * (1 + C) / (ratio + C);
396 ASSERTL0(
false,
"Invalid variables for WSS");
404 int npoints = exp[0]->GetNpoints();
405 if (boost::iequals(
m_f->m_variables[0],
"u"))
414 else if (boost::iequals(
m_f->m_variables[0],
"rho") &&
415 boost::iequals(
m_f->m_variables[1],
"rhou"))
421 Vmath::Vdiv(npoints, exp[i + 1]->GetPhys(), 1, exp[0]->GetPhys(), 1,
428 ASSERTL0(
false,
"Could not identify velocity for WSS");
#define ASSERTL0(condition, msg)
FieldSharedPtr m_f
Field object.
void GetViscosity(const Array< OneD, MultiRegions::ExpListSharedPtr > exp, Array< OneD, NekDouble > &mu, NekDouble &lambda)
virtual void Process(po::variables_map &vm)
Write mesh to output file.
void GetVelocity(const Array< OneD, MultiRegions::ExpListSharedPtr > exp, Array< OneD, Array< OneD, NekDouble >> &vel)
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
const BoundaryRegionCollection & GetBoundaryRegions(void) const
std::shared_ptr< Field > FieldSharedPtr
std::pair< ModuleType, std::string > ModuleKey
ModuleFactory & GetModuleFactory()
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
The above copyright notice and this permission notice shall be included.
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
void Neg(int n, T *x, const int incx)
Negate x = -x.
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Vdiv(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x/y.
void Zero(int n, T *x, const int incx)
Zero vector.
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
scalarT< T > sqrt(scalarT< T > in)