Nektar++
Public Member Functions | Static Public Member Functions | Static Public Attributes | Protected Member Functions | Private Attributes | List of all members
Nektar::FieldUtils::ProcessWSS Class Reference

This processing module calculates the wall shear stress and adds it as an extra-field to the output file, and writes it to a surface output file. More...

#include <ProcessWSS.h>

Inheritance diagram for Nektar::FieldUtils::ProcessWSS:
[legend]

Public Member Functions

 ProcessWSS (FieldSharedPtr f)
 
virtual ~ProcessWSS ()
 
virtual void Process (po::variables_map &vm)
 Write mesh to output file. More...
 
virtual std::string GetModuleName ()
 
virtual std::string GetModuleDescription ()
 
- Public Member Functions inherited from Nektar::FieldUtils::ProcessBoundaryExtract
 ProcessBoundaryExtract (FieldSharedPtr f)
 
virtual ~ProcessBoundaryExtract ()
 
virtual ModulePriority GetModulePriority ()
 
- Public Member Functions inherited from Nektar::FieldUtils::ProcessModule
 ProcessModule ()
 
 ProcessModule (FieldSharedPtr p_f)
 
- Public Member Functions inherited from Nektar::FieldUtils::Module
FIELD_UTILS_EXPORT Module (FieldSharedPtr p_f)
 
virtual ~Module ()=default
 
const ConfigOptionGetConfigOption (const std::string &key) const
 
FIELD_UTILS_EXPORT void RegisterConfig (std::string key, std::string value="")
 Register a configuration option with a module. More...
 
FIELD_UTILS_EXPORT void PrintConfig ()
 Print out all configuration options for a module. More...
 
FIELD_UTILS_EXPORT void SetDefaults ()
 Sets default configuration options for those which have not been set. More...
 
FIELD_UTILS_EXPORT void AddFile (std::string fileType, std::string fileName)
 
FIELD_UTILS_EXPORT void EvaluateTriFieldAtEquiSpacedPts (LocalRegions::ExpansionSharedPtr &exp, const Array< OneD, const NekDouble > &infield, Array< OneD, NekDouble > &outfield)
 

Static Public Member Functions

static std::shared_ptr< Modulecreate (FieldSharedPtr f)
 Creates an instance of this class. More...
 
- Static Public Member Functions inherited from Nektar::FieldUtils::ProcessBoundaryExtract
static std::shared_ptr< Modulecreate (FieldSharedPtr f)
 Creates an instance of this class. More...
 

Static Public Attributes

static ModuleKey className
 
- Static Public Attributes inherited from Nektar::FieldUtils::ProcessBoundaryExtract
static ModuleKey className
 

Protected Member Functions

void GetViscosity (const Array< OneD, MultiRegions::ExpListSharedPtr > exp, Array< OneD, NekDouble > &mu, NekDouble &lambda)
 
void GetVelocity (const Array< OneD, MultiRegions::ExpListSharedPtr > exp, Array< OneD, Array< OneD, NekDouble >> &vel)
 
- Protected Member Functions inherited from Nektar::FieldUtils::Module
 Module ()
 

Private Attributes

int m_spacedim
 

Additional Inherited Members

- Public Attributes inherited from Nektar::FieldUtils::Module
FieldSharedPtr m_f
 Field object. More...
 
- Protected Attributes inherited from Nektar::FieldUtils::Module
std::map< std::string, ConfigOptionm_config
 List of configuration values. More...
 
std::set< std::string > m_allowedFiles
 List of allowed file formats. More...
 

Detailed Description

This processing module calculates the wall shear stress and adds it as an extra-field to the output file, and writes it to a surface output file.

Definition at line 49 of file ProcessWSS.h.

Constructor & Destructor Documentation

◆ ProcessWSS()

Nektar::FieldUtils::ProcessWSS::ProcessWSS ( FieldSharedPtr  f)

Definition at line 54 of file ProcessWSS.cpp.

◆ ~ProcessWSS()

Nektar::FieldUtils::ProcessWSS::~ProcessWSS ( )
virtual

Definition at line 58 of file ProcessWSS.cpp.

59 {
60 }

Member Function Documentation

◆ create()

static std::shared_ptr<Module> Nektar::FieldUtils::ProcessWSS::create ( FieldSharedPtr  f)
inlinestatic

Creates an instance of this class.

Definition at line 53 of file ProcessWSS.h.

54  {
56  }
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

◆ GetModuleDescription()

virtual std::string Nektar::FieldUtils::ProcessWSS::GetModuleDescription ( )
inlinevirtual

Reimplemented from Nektar::FieldUtils::ProcessBoundaryExtract.

Definition at line 70 of file ProcessWSS.h.

71  {
72  return "Calculating wall shear stress";
73  }

◆ GetModuleName()

virtual std::string Nektar::FieldUtils::ProcessWSS::GetModuleName ( )
inlinevirtual

Reimplemented from Nektar::FieldUtils::ProcessBoundaryExtract.

Definition at line 65 of file ProcessWSS.h.

66  {
67  return "ProcessWSS";
68  }

◆ GetVelocity()

void Nektar::FieldUtils::ProcessWSS::GetVelocity ( const Array< OneD, MultiRegions::ExpListSharedPtr exp,
Array< OneD, Array< OneD, NekDouble >> &  vel 
)
protected

Definition at line 400 of file ProcessWSS.cpp.

403 {
404  int npoints = exp[0]->GetNpoints();
405  if (boost::iequals(m_f->m_variables[0], "u"))
406  {
407  // IncNavierStokesSolver
408  for (int i = 0; i < m_spacedim; ++i)
409  {
410  vel[i] = Array<OneD, NekDouble>(npoints);
411  Vmath::Vcopy(npoints, exp[i]->GetPhys(), 1, vel[i], 1);
412  }
413  }
414  else if (boost::iequals(m_f->m_variables[0], "rho") &&
415  boost::iequals(m_f->m_variables[1], "rhou"))
416  {
417  // CompressibleFlowSolver
418  for (int i = 0; i < m_spacedim; ++i)
419  {
420  vel[i] = Array<OneD, NekDouble>(npoints);
421  Vmath::Vdiv(npoints, exp[i + 1]->GetPhys(), 1, exp[0]->GetPhys(), 1,
422  vel[i], 1);
423  }
424  }
425  else
426  {
427  // Unknown
428  ASSERTL0(false, "Could not identify velocity for WSS");
429  }
430 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
FieldSharedPtr m_f
Field object.
Definition: Module.h:225
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.
Definition: Vmath.cpp:284
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255

References ASSERTL0, Nektar::FieldUtils::Module::m_f, m_spacedim, Vmath::Vcopy(), and Vmath::Vdiv().

Referenced by Process().

◆ GetViscosity()

void Nektar::FieldUtils::ProcessWSS::GetViscosity ( const Array< OneD, MultiRegions::ExpListSharedPtr exp,
Array< OneD, NekDouble > &  mu,
NekDouble lambda 
)
protected

Definition at line 296 of file ProcessWSS.cpp.

299 {
300  NekDouble m_mu;
301  int npoints = exp[0]->GetNpoints();
302 
303  if (boost::iequals(m_f->m_variables[0], "u"))
304  {
305  // IncNavierStokesSolver
306  m_mu = m_f->m_session->GetParameter("Kinvis");
307  Vmath::Fill(npoints, m_mu, mu, 1);
308  lambda = 0;
309  }
310  else if (boost::iequals(m_f->m_variables[0], "rho") &&
311  boost::iequals(m_f->m_variables[1], "rhou"))
312  {
313  // CompressibleFlowSolver
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,
318  "Constant");
319 
320  if (m_ViscosityType == "Variable")
321  {
322  // Check equation of state
323  std::string eosType;
324  bool m_idealGas;
325  m_f->m_session->LoadSolverInfo("EquationOfState", eosType,
326  "IdealGas");
327  m_idealGas = boost::iequals(eosType, "IdealGas");
328  ASSERTL0(
329  m_idealGas,
330  "Only IdealGas EOS implemented for Variable ViscosityType");
331 
332  // Get relevant parameters
333  NekDouble m_gamma;
334  NekDouble m_pInf;
336  NekDouble m_gasConstant;
337  NekDouble cv_inv;
338  NekDouble m_Tref;
339  NekDouble m_TRatioSutherland;
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,
344  287.058);
345  m_f->m_session->LoadParameter("Tref", m_Tref, 288.15);
346  m_TRatioSutherland = 110.0 / m_Tref;
347 
348  // Get temperature from flowfield
349  cv_inv = (m_gamma - 1.0) / m_gasConstant;
350  // e = 1/rho ( E - 1/2 ( rhou^2/rho + ... ) )
351  Array<OneD, NekDouble> tmp(npoints, 0.0);
352  Array<OneD, NekDouble> energy(npoints, 0.0);
353  Array<OneD, NekDouble> temperature(npoints, 0.0);
354  Vmath::Vcopy(npoints, exp[m_spacedim + 1]->GetPhys(), 1, energy, 1);
355  for (int i = 0; i < m_spacedim; i++)
356  {
357  // rhou^2
358  Vmath::Vmul(npoints, exp[i + 1]->GetPhys(), 1,
359  exp[i + 1]->GetPhys(), 1, tmp, 1);
360  // rhou^2/rho
361  Vmath::Vdiv(npoints, tmp, 1, exp[0]->GetPhys(), 1, tmp, 1);
362  // 0.5 rhou^2/rho
363  Vmath::Smul(npoints, 0.5, tmp, 1, tmp, 1);
364  // E - 0.5 rhou^2/rho - ...
365  Vmath::Vsub(npoints, energy, 1, tmp, 1, energy, 1);
366  }
367  // rhoe/rho
368  Vmath::Vdiv(npoints, energy, 1, exp[0]->GetPhys(), 1, energy, 1);
369  // T = e/Cv
370  Vmath::Smul(npoints, cv_inv, energy, 1, temperature, 1);
371 
372  // Variable viscosity through the Sutherland's law
373  //
374  // WARNING, if this routine is modified the same must be done in the
375  // CompressibleFlowSolver function in VariableConverter.cpp
376  // (this class should be restructured).
377 
378  const NekDouble C = m_TRatioSutherland;
379  NekDouble mu_star = m_mu;
380  NekDouble T_star = m_pInf / (m_rhoInf * m_gasConstant);
381  NekDouble ratio;
382  for (int i = 0; i < npoints; ++i)
383  {
384  ratio = temperature[i] / T_star;
385  mu[i] = mu_star * ratio * sqrt(ratio) * (1 + C) / (ratio + C);
386  }
387  }
388  else
389  {
390  Vmath::Fill(npoints, m_mu, mu, 1);
391  }
392  }
393  else
394  {
395  // Unknown
396  ASSERTL0(false, "Invalid variables for WSS");
397  }
398 }
NekDouble m_mu
NekDouble m_rhoInf
double NekDouble
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.
Definition: Vmath.cpp:209
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
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.
Definition: Vmath.cpp:419
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291

References ASSERTL0, Vmath::Fill(), Nektar::FieldUtils::Module::m_f, m_mu, m_rhoInf, m_spacedim, Vmath::Smul(), tinysimd::sqrt(), Vmath::Vcopy(), Vmath::Vdiv(), Vmath::Vmul(), and Vmath::Vsub().

Referenced by Process().

◆ Process()

void Nektar::FieldUtils::ProcessWSS::Process ( po::variables_map &  vm)
virtual

Write mesh to output file.

Reimplemented from Nektar::FieldUtils::ProcessBoundaryExtract.

Definition at line 62 of file ProcessWSS.cpp.

63 {
65 
66  int i, j;
67  int nfields = m_f->m_variables.size();
68  int expdim = m_f->m_graph->GetSpaceDimension();
69  m_spacedim = expdim + m_f->m_numHomogeneousDir;
70 
71  if (m_f->m_exp[0]->GetNumElmts() == 0)
72  {
73  return;
74  }
75 
76  if (m_spacedim == 1)
77  {
78  ASSERTL0(false, "Error: wss for a 1D problem cannot "
79  "be computed");
80  }
81 
82  // Declare arrays
83  int nshear = m_spacedim + 1;
84  int nstress = m_spacedim * m_spacedim;
85  int ngrad = m_spacedim * m_spacedim;
86 
87  Array<OneD, Array<OneD, NekDouble>> velocity(nfields);
88  Array<OneD, Array<OneD, NekDouble>> grad(ngrad);
89  Array<OneD, Array<OneD, NekDouble>> stress(nstress), fstress(nstress);
90  Array<OneD, Array<OneD, NekDouble>> fshear(nshear);
91 
92  Array<OneD, MultiRegions::ExpListSharedPtr> BndExp(nshear);
93  Array<OneD, MultiRegions::ExpListSharedPtr> BndElmtExp(nfields);
94 
95  // will resuse nfields expansions to write shear components.
96  if (nshear > nfields)
97  {
98  m_f->m_exp.resize(nshear);
99  for (i = nfields; i < nshear; ++i)
100  {
101  m_f->m_exp[nfields + i] =
102  m_f->AppendExpList(m_f->m_numHomogeneousDir);
103  }
104  }
105 
106  // Create map of boundary ids for partitioned domains
107  SpatialDomains::BoundaryConditions bcs(m_f->m_session,
108  m_f->m_exp[0]->GetGraph());
110  bcs.GetBoundaryRegions();
111  map<int, int> BndRegionMap;
112  int cnt = 0;
113  for (auto &breg_it : bregions)
114  {
115  BndRegionMap[breg_it.first] = cnt++;
116  }
117 
118  // Loop over boundaries to Write
119  for (int b = 0; b < m_f->m_bndRegionsToWrite.size(); ++b)
120  {
121  if (BndRegionMap.count(m_f->m_bndRegionsToWrite[b]) == 1)
122  {
123  int bnd = BndRegionMap[m_f->m_bndRegionsToWrite[b]];
124  // Get expansion list for boundary and for elements containing this
125  // bnd
126  for (i = 0; i < nshear; i++)
127  {
128  BndExp[i] = m_f->m_exp[i]->UpdateBndCondExpansion(bnd);
129  }
130  for (i = 0; i < nfields; i++)
131  {
132  m_f->m_exp[i]->GetBndElmtExpansion(bnd, BndElmtExp[i]);
133  }
134 
135  // Get number of points in expansions
136  int nqb = BndExp[0]->GetTotPoints();
137  int nqe = BndElmtExp[0]->GetTotPoints();
138 
139  // Initialise local arrays for the velocity gradients, and
140  // stress components size of total number of quadrature
141  // points for elements in this bnd
142  for (i = 0; i < ngrad; ++i)
143  {
144  grad[i] = Array<OneD, NekDouble>(nqe);
145  }
146 
147  for (i = 0; i < nstress; ++i)
148  {
149  stress[i] = Array<OneD, NekDouble>(nqe);
150  }
151 
152  Array<OneD, NekDouble> div(nqe, 0.0);
153 
154  // initialise arrays in the boundary
155  for (i = 0; i < nstress; ++i)
156  {
157  fstress[i] = Array<OneD, NekDouble>(nqb);
158  }
159 
160  for (i = 0; i < nshear; ++i)
161  {
162  fshear[i] = Array<OneD, NekDouble>(nqb, 0.0);
163  }
164 
165  // Extract Velocities
166  GetVelocity(BndElmtExp, velocity);
167 
168  // Extract viscosity coefficients
169  NekDouble lambda;
170  Array<OneD, NekDouble> mu(nqe, 0.0);
171  GetViscosity(BndElmtExp, mu, lambda);
172 
173  // Compute gradients
174  for (i = 0; i < m_spacedim; ++i)
175  {
176  if (m_spacedim == 2)
177  {
178  BndElmtExp[i]->PhysDeriv(velocity[i],
179  grad[i * m_spacedim + 0],
180  grad[i * m_spacedim + 1]);
181  }
182  else
183  {
184  BndElmtExp[i]->PhysDeriv(
185  velocity[i], grad[i * m_spacedim + 0],
186  grad[i * m_spacedim + 1], grad[i * m_spacedim + 2]);
187  }
188  // Add contribution to div(u)
189  Vmath::Vadd(nqe, grad[i * m_spacedim + i], 1, div, 1, div, 1);
190  }
191 
192  // Velocity divergence scaled by lambda * mu
193  Vmath::Smul(nqe, lambda, div, 1, div, 1);
194  Vmath::Vmul(nqe, mu, 1, div, 1, div, 1);
195 
196  // Compute stress component terms
197  // tau_ij = mu*(u_i,j + u_j,i) + mu*lambda*delta_ij*div(u)
198  for (i = 0; i < m_spacedim; ++i)
199  {
200  for (j = i; j < m_spacedim; ++j)
201  {
202  Vmath::Vadd(nqe, grad[i * m_spacedim + j], 1,
203  grad[j * m_spacedim + i], 1,
204  stress[i * m_spacedim + j], 1);
205 
206  Vmath::Vmul(nqe, mu, 1, stress[i * m_spacedim + j], 1,
207  stress[i * m_spacedim + j], 1);
208 
209  if (i == j)
210  {
211  // Add divergence term to diagonal
212  Vmath::Vadd(nqe, stress[i * m_spacedim + j], 1, div, 1,
213  stress[i * m_spacedim + j], 1);
214  }
215  else
216  {
217  // Copy to make symmetric
218  Vmath::Vcopy(nqe, stress[i * m_spacedim + j], 1,
219  stress[j * m_spacedim + i], 1);
220  }
221  }
222  }
223 
224  // Get boundary stress values.
225  for (j = 0; j < nstress; ++j)
226  {
227  m_f->m_exp[0]->ExtractElmtToBndPhys(bnd, stress[j], fstress[j]);
228  }
229 
230  // Get normals
231  Array<OneD, Array<OneD, NekDouble>> normals;
232  m_f->m_exp[0]->GetBoundaryNormals(bnd, normals);
233  // Reverse normals, to get correct orientation for the body
234  for (i = 0; i < m_spacedim; ++i)
235  {
236  Vmath::Neg(nqb, normals[i], 1);
237  }
238 
239  // calculate wss, and update coeffs in the boundary expansion
240  // S = tau_ij * n_j
241  for (i = 0; i < m_spacedim; ++i)
242  {
243  for (j = 0; j < m_spacedim; ++j)
244  {
245  Vmath::Vvtvp(nqb, normals[j], 1,
246  fstress[i * m_spacedim + j], 1, fshear[i], 1,
247  fshear[i], 1);
248  }
249  }
250 
251  // T = S - (S.n)n
252  for (i = 0; i < m_spacedim; ++i)
253  {
254  Vmath::Vvtvp(nqb, normals[i], 1, fshear[i], 1,
255  fshear[nshear - 1], 1, fshear[nshear - 1], 1);
256  }
257  Vmath::Smul(nqb, -1.0, fshear[nshear - 1], 1, fshear[nshear - 1],
258  1);
259 
260  for (i = 0; i < m_spacedim; i++)
261  {
262  Vmath::Vvtvp(nqb, normals[i], 1, fshear[nshear - 1], 1,
263  fshear[i], 1, fshear[i], 1);
264  BndExp[i]->FwdTransLocalElmt(fshear[i],
265  BndExp[i]->UpdateCoeffs());
266  }
267 
268  // Tw
269  Vmath::Zero(nqb, fshear[nshear - 1], 1);
270  for (i = 0; i < m_spacedim; ++i)
271  {
272  Vmath::Vvtvp(nqb, fshear[i], 1, fshear[i], 1,
273  fshear[nshear - 1], 1, fshear[nshear - 1], 1);
274  }
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());
278  }
279  }
280 
281  if (m_spacedim == 2)
282  {
283  m_f->m_variables[0] = "Shear_s";
284  m_f->m_variables[1] = "Shear_n";
285  m_f->m_variables[2] = "Shear_mag";
286  }
287  else
288  {
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";
293  }
294 }
virtual void Process(po::variables_map &vm)
void GetViscosity(const Array< OneD, MultiRegions::ExpListSharedPtr > exp, Array< OneD, NekDouble > &mu, NekDouble &lambda)
Definition: ProcessWSS.cpp:296
void GetVelocity(const Array< OneD, MultiRegions::ExpListSharedPtr > exp, Array< OneD, Array< OneD, NekDouble >> &vel)
Definition: ProcessWSS.cpp:400
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:210
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:534
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
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
Definition: Vmath.cpp:574
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.
Definition: Vmath.cpp:359
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492

References ASSERTL0, Nektar::SpatialDomains::BoundaryConditions::GetBoundaryRegions(), GetVelocity(), GetViscosity(), Nektar::FieldUtils::Module::m_f, m_spacedim, Vmath::Neg(), Nektar::FieldUtils::ProcessBoundaryExtract::Process(), Vmath::Smul(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vvtvp(), and Vmath::Zero().

Member Data Documentation

◆ className

ModuleKey Nektar::FieldUtils::ProcessWSS::className
static
Initial value:
"Computes wall shear stress field.")
static std::shared_ptr< Module > create(FieldSharedPtr f)
Creates an instance of this class.
Definition: ProcessWSS.h:53
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
std::pair< ModuleType, std::string > ModuleKey
Definition: Module.h:285
ModuleFactory & GetModuleFactory()
Definition: Module.cpp:49

Definition at line 57 of file ProcessWSS.h.

◆ m_spacedim

int Nektar::FieldUtils::ProcessWSS::m_spacedim
private

Definition at line 83 of file ProcessWSS.h.

Referenced by GetVelocity(), GetViscosity(), and Process().