Nektar++
Functions
ExtractSurface3DCFS.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <string>
#include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
#include <MultiRegions/ExpList.h>
#include <LocalRegions/Expansion.h>
#include <LocalRegions/Expansion3D.h>
#include <LocalRegions/MatrixKey.h>
#include <MultiRegions/DisContField.h>
#include <LibUtilities/BasicUtils/FieldIO.h>
#include <LibUtilities/BasicUtils/NekFactory.hpp>
#include <LibUtilities/BasicUtils/SessionReader.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/Communication/Comm.h>
#include <LibUtilities/Memory/NekMemoryManager.hpp>
#include <MultiRegions/ContField.h>
#include <SpatialDomains/MeshGraph.h>
#include <SolverUtils/SolverUtilsDeclspec.h>

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Evaluation of the velocity gradient in the cartesian directions Du_x: traceFieldsAdded[17] Du_y: traceFieldsAdded[18] Du_z: traceFieldsAdded[19] Dv_x: traceFieldsAdded[20] Dv_y: traceFieldsAdded[21] Dv_z: traceFieldsAdded[22] Dw_x: traceFieldsAdded[23] Dw_y: traceFieldsAdded[24] Dw_z: traceFieldsAdded[25]

Definition at line 65 of file ExtractSurface3DCFS.cpp.

66 {
67  string fname = std::string(argv[2]);
68  int fdot = fname.find_last_of('.');
69  if (fdot != std::string::npos)
70  {
71  string ending = fname.substr(fdot);
72 
73  // If .chk or .fld we exchange the extension in the output file.
74  // For all other files (e.g. .bse) we append the extension to avoid
75  // conflicts.
76  if (ending == ".chk" || ending == ".fld")
77  {
78  fname = fname.substr(0, fdot);
79  }
80  }
81 
82  fname = fname + ".txt";
83 
84  int cnt;
85  int id1, id2;
86  int i, j, n, e, b;
87  Array<OneD, NekDouble> auxArray;
88 
89  int nBndEdgePts, nBndEdges, nBndRegions;
90 
91  if (argc < 3)
92  {
93  fprintf(stderr, "Usage: ExtractSurface3DCFS meshfile fieldFile\n");
94  fprintf(stderr,
95  "Extracts a surface from a 3D fld file"
96  "(only for CompressibleFlowSolver and purely 3D .fld files)\n");
97  exit(1);
98  }
99 
100  LibUtilities::SessionReaderSharedPtr vSession =
101  LibUtilities::SessionReader::CreateInstance(3, argv);
102  SpatialDomains::MeshGraphSharedPtr graphShPt =
103  SpatialDomains::MeshGraph::Read(vSession);
104 
105  std::string m_ViscosityType;
106 
107  NekDouble m_gamma;
108  NekDouble m_pInf;
109  NekDouble m_rhoInf;
110  NekDouble m_uInf;
111  NekDouble m_vInf;
112  NekDouble m_wInf;
113  NekDouble m_gasConstant;
114  NekDouble m_Twall;
115  NekDouble m_mu;
116 
117  int m_spacedim = 3;
118  int nDimensions = m_spacedim;
119  int phys_offset;
120 
121  // Get gamma parameter from session file.
122  ASSERTL0(vSession->DefinesParameter("Gamma"),
123  "Compressible flow sessions must define a Gamma parameter.");
124  vSession->LoadParameter("Gamma", m_gamma, 1.4);
125 
126  // Get E0 parameter from session file.
127  ASSERTL0(vSession->DefinesParameter("pInf"),
128  "Compressible flow sessions must define a pInf parameter.");
129  vSession->LoadParameter("pInf", m_pInf, 101325);
130 
131  // Get rhoInf parameter from session file.
132  ASSERTL0(vSession->DefinesParameter("rhoInf"),
133  "Compressible flow sessions must define a rhoInf parameter.");
134  vSession->LoadParameter("rhoInf", m_rhoInf, 1.225);
135 
136  // Get uInf parameter from session file.
137  ASSERTL0(vSession->DefinesParameter("uInf"),
138  "Compressible flow sessions must define a uInf parameter.");
139  vSession->LoadParameter("uInf", m_uInf, 0.1);
140 
141  // Get vInf parameter from session file.
142  if (m_spacedim == 2 || m_spacedim == 3)
143  {
144  ASSERTL0(vSession->DefinesParameter("vInf"),
145  "Compressible flow sessions must define a vInf parameter"
146  "for 2D/3D problems.");
147  vSession->LoadParameter("vInf", m_vInf, 0.0);
148  }
149 
150  // Get wInf parameter from session file.
151  if (m_spacedim == 3)
152  {
153  ASSERTL0(vSession->DefinesParameter("wInf"),
154  "Compressible flow sessions must define a wInf parameter"
155  "for 3D problems.");
156  vSession->LoadParameter("wInf", m_wInf, 0.0);
157  }
158 
159  vSession->LoadParameter("GasConstant", m_gasConstant, 287.058);
160  vSession->LoadParameter("Twall", m_Twall, 300.15);
161  vSession->LoadSolverInfo("ViscosityType", m_ViscosityType, "Constant");
162  vSession->LoadParameter("mu", m_mu, 1.78e-05);
163 
164  //--------------------------------------------------------------------------
165  // Import field file
166  string fieldFile(argv[2]);
167  vector<LibUtilities::FieldDefinitionsSharedPtr> fieldDef;
168  vector<vector<NekDouble>> fieldData;
169 
170  LibUtilities::Import(fieldFile, fieldDef, fieldData);
171  //--------------------------------------------------------------------------
172 
173  //--------------------------------------------------------------------------
174  // Set up Expansion information
175  vector<vector<LibUtilities::PointsType>> pointsType;
176  for (i = 0; i < fieldDef.size(); ++i)
177  {
178  vector<LibUtilities::PointsType> ptype;
179  for (j = 0; j < 3; ++j)
180  {
181  ptype.push_back(LibUtilities::ePolyEvenlySpaced);
182  }
183  pointsType.push_back(ptype);
184  }
185  graphShPt->SetExpansionInfo(fieldDef, pointsType);
186 
187  //--------------------------------------------------------------------------
188 
189  //--------------------------------------------------------------------------
190  // Define Expansion
191  int nfields = fieldDef[0]->m_fields.size();
192  Array<OneD, MultiRegions::ExpListSharedPtr> Exp(nfields);
193  Array<OneD, MultiRegions::ExpListSharedPtr> pFields(nfields);
194 
195  for (i = 0; i < pFields.size(); i++)
196  {
197  pFields[i] =
198  MemoryManager<MultiRegions ::DisContField>::AllocateSharedPtr(
199  vSession, graphShPt, vSession->GetVariable(i));
200  }
201 
202  MultiRegions::ExpListSharedPtr Exp3D;
203  Exp3D = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(vSession,
204  graphShPt);
205 
206  Exp[0] = Exp3D;
207 
208  for (i = 1; i < nfields; ++i)
209  {
210  Exp[i] =
211  MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(*Exp3D);
212  }
213 
214  // Count of the point on the surface
215  int nSurfacePts = 0;
216  if (pFields[0]->GetBndCondExpansions().size())
217  {
218  nSurfacePts = 0;
219  cnt = 0;
220  nBndRegions = pFields[0]->GetBndCondExpansions().size();
221  for (b = 0; b < nBndRegions; ++b)
222  {
223  nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
224  for (e = 0; e < nBndEdges; ++e)
225  {
226  nBndEdgePts = pFields[0]
227  ->GetBndCondExpansions()[b]
228  ->GetExp(e)
229  ->GetTotPoints();
230 
231  if (pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
232  "WallViscous" ||
233  pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
234  "WallAdiabatic" ||
235  pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
236  "Wall")
237  {
238  nSurfacePts += nBndEdgePts;
239  }
240  }
241  }
242  }
243 
244  int nSolutionPts = pFields[0]->GetNpoints();
245  int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
246  int nElements = pFields[0]->GetExpSize();
247 
248  Array<OneD, NekDouble> tmp(nSolutionPts, 0.0);
249 
250  Array<OneD, NekDouble> x(nSolutionPts);
251  Array<OneD, NekDouble> y(nSolutionPts);
252  Array<OneD, NekDouble> z(nSolutionPts);
253 
254  Array<OneD, NekDouble> traceX(nTracePts);
255  Array<OneD, NekDouble> traceY(nTracePts);
256  Array<OneD, NekDouble> traceZ(nTracePts);
257 
258  Array<OneD, NekDouble> surfaceX(nSurfacePts);
259  Array<OneD, NekDouble> surfaceY(nSurfacePts);
260  Array<OneD, NekDouble> surfaceZ(nSurfacePts);
261 
262  pFields[0]->GetCoords(x, y, z);
263 
264  pFields[0]->ExtractTracePhys(x, traceX);
265  pFields[0]->ExtractTracePhys(y, traceY);
266  pFields[0]->ExtractTracePhys(z, traceZ);
267  //--------------------------------------------------------------------------
268 
269  //--------------------------------------------------------------------------
270  // Copy data from field file
271  Array<OneD, Array<OneD, NekDouble>> uFields(nfields);
272  Array<OneD, Array<OneD, NekDouble>> traceFields(nfields);
273  Array<OneD, Array<OneD, NekDouble>> surfaceFields(nfields);
274 
275  // Extract the physical values of the solution at the boundaries
276  for (j = 0; j < nfields; ++j)
277  {
278  uFields[j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
279  traceFields[j] = Array<OneD, NekDouble>(nTracePts, 0.0);
280  surfaceFields[j] = Array<OneD, NekDouble>(nSurfacePts, 0.0);
281 
282  for (i = 0; i < fieldData.size(); ++i)
283  {
284  Exp[j]->ExtractDataToCoeffs(fieldDef[i], fieldData[i],
285  fieldDef[i]->m_fields[j],
286  Exp[j]->UpdateCoeffs());
287  }
288  Exp[j]->BwdTrans(Exp[j]->GetCoeffs(), Exp[j]->UpdatePhys());
289  Vmath::Vcopy(nSolutionPts, Exp[j]->GetPhys(), 1, uFields[j], 1);
290  pFields[0]->ExtractTracePhys(uFields[j], traceFields[j]);
291  }
292 
293  // Fields to add in the output file
294  int nfieldsAdded = 34;
295  Array<OneD, Array<OneD, NekDouble>> traceFieldsAdded(nfieldsAdded);
296  Array<OneD, Array<OneD, NekDouble>> surfaceFieldsAdded(nfieldsAdded);
297 
298  for (j = 0; j < nfieldsAdded; ++j)
299  {
300  traceFieldsAdded[j] = Array<OneD, NekDouble>(nTracePts, 0.0);
301  surfaceFieldsAdded[j] = Array<OneD, NekDouble>(nSurfacePts, 0.0);
302  }
303 
304  /******** Evaluation of normals and tangents on the trace *****************
305  * nx -> traceFieldsAdded[0];
306  * ny -> traceFieldsAdded[1];
307  * nz -> traceFieldsAdded[2];
308  * bx -> traceFieldsAdded[3];
309  * by -> traceFieldsAdded[4];
310  * bz -> traceFieldsAdded[5];
311  * tx -> traceFieldsAdded[6];
312  * ty -> traceFieldsAdded[7];
313  * tz -> traceFieldsAdded[8];
314 
315  ***************************************************************************/
316 
317  Array<OneD, Array<OneD, NekDouble>> m_traceNormals(nDimensions);
318  for (i = 0; i < nDimensions; ++i)
319  {
320  m_traceNormals[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
321  }
322  pFields[0]->GetTrace()->GetNormals(m_traceNormals);
323 
324  Array<OneD, Array<OneD, NekDouble>> m_traceTangents(nDimensions);
325  Array<OneD, Array<OneD, NekDouble>> m_traceBinormals(nDimensions);
326  Array<OneD, Array<OneD, NekDouble>> h(nDimensions);
327  Array<OneD, NekDouble> tmpNorm(nTracePts, 1.0);
328  Array<OneD, NekDouble> NormH(nTracePts, 0.0);
329  Array<OneD, NekDouble> tmpTrace(nTracePts, 0.0);
330 
331  for (i = 0; i < nDimensions; ++i)
332  {
333  m_traceTangents[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
334  m_traceBinormals[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
335  h[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
336  }
337 
338  // Normals
339 
340  // nx
341  Vmath::Vcopy(nTracePts, &m_traceNormals[0][0], 1, &traceFieldsAdded[0][0],
342  1);
343 
344  // ny
345  Vmath::Vcopy(nTracePts, &m_traceNormals[1][0], 1, &traceFieldsAdded[1][0],
346  1);
347 
348  // nz
349  Vmath::Vcopy(nTracePts, &m_traceNormals[2][0], 1, &traceFieldsAdded[2][0],
350  1);
351 
352  // Tangents and Binormals
353  // h1
354  Vmath::Vadd(nTracePts, &m_traceNormals[0][0], 1, &tmpNorm[0], 1, &h[0][0],
355  1);
356  // h2
357  Vmath::Vcopy(nTracePts, &m_traceNormals[1][0], 1, &h[1][0], 1);
358  // h3
359  Vmath::Vcopy(nTracePts, &m_traceNormals[2][0], 1, &h[2][0], 1);
360 
361  // Norm of h
362  for (i = 0; i < m_spacedim; i++)
363  {
364  Vmath::Vvtvp(nTracePts, &h[i][0], 1, &h[i][0], 1, &NormH[0], 1,
365  &NormH[0], 1);
366  }
367 
368  // b1
369  Vmath::Vmul(nTracePts, &h[0][0], 1, &h[1][0], 1, &tmpTrace[0], 1);
370 
371  Vmath::Vdiv(nTracePts, &tmpTrace[0], 1, &NormH[0], 1, &tmpTrace[0], 1);
372 
373  Vmath::Smul(nTracePts, -2.0, &tmpTrace[0], 1, &m_traceBinormals[0][0], 1);
374 
375  Vmath::Vcopy(nTracePts, &m_traceBinormals[0][0], 1, &traceFieldsAdded[3][0],
376  1);
377 
378  // b2
379  Vmath::Vmul(nTracePts, &h[1][0], 1, &h[1][0], 1, &tmpTrace[0], 1);
380 
381  Vmath::Vdiv(nTracePts, &tmpTrace[0], 1, &NormH[0], 1, &tmpTrace[0], 1);
382 
383  Vmath::Smul(nTracePts, -2.0, &tmpTrace[0], 1, &tmpTrace[0], 1);
384 
385  Vmath::Vadd(nTracePts, &tmpTrace[0], 1, &tmpNorm[0], 1,
386  &m_traceBinormals[1][0], 1);
387 
388  Vmath::Vcopy(nTracePts, &m_traceBinormals[1][0], 1, &traceFieldsAdded[4][0],
389  1);
390 
391  // b3
392  Vmath::Vmul(nTracePts, &h[1][0], 1, &h[2][0], 1, &tmpTrace[0], 1);
393 
394  Vmath::Vdiv(nTracePts, &tmpTrace[0], 1, &NormH[0], 1, &tmpTrace[0], 1);
395 
396  Vmath::Smul(nTracePts, -2.0, &tmpTrace[0], 1, &m_traceBinormals[2][0], 1);
397 
398  Vmath::Vcopy(nTracePts, &m_traceBinormals[2][0], 1, &traceFieldsAdded[5][0],
399  1);
400 
401  // t1
402  Vmath::Vmul(nTracePts, &h[0][0], 1, &h[2][0], 1, &tmpTrace[0], 1);
403 
404  Vmath::Vdiv(nTracePts, &tmpTrace[0], 1, &NormH[0], 1, &tmpTrace[0], 1);
405 
406  Vmath::Smul(nTracePts, -2.0, &tmpTrace[0], 1, &m_traceTangents[0][0], 1);
407 
408  Vmath::Vcopy(nTracePts, &m_traceTangents[0][0], 1, &traceFieldsAdded[6][0],
409  1);
410 
411  // t2
412  Vmath::Vcopy(nTracePts, &m_traceBinormals[2][0], 1, &m_traceTangents[1][0],
413  1);
414 
415  Vmath::Vcopy(nTracePts, &m_traceTangents[1][0], 1, &traceFieldsAdded[7][0],
416  1);
417 
418  // t3
419  Vmath::Vmul(nTracePts, &h[2][0], 1, &h[2][0], 1, &tmpTrace[0], 1);
420 
421  Vmath::Vdiv(nTracePts, &tmpTrace[0], 1, &NormH[0], 1, &tmpTrace[0], 1);
422 
423  Vmath::Smul(nTracePts, -2.0, &tmpTrace[0], 1, &tmpTrace[0], 1);
424 
425  Vmath::Vadd(nTracePts, &tmpTrace[0], 1, &tmpNorm[0], 1,
426  &m_traceTangents[2][0], 1);
427 
428  Vmath::Vcopy(nTracePts, &m_traceTangents[2][0], 1, &traceFieldsAdded[8][0],
429  1);
430 
431  /******** Evaluation of the pressure ***************************************
432  * P = (E-1/2.*rho.*((rhou./rho).^2+(rhov./rho).^2))*(gamma - 1);
433  * P -> traceFieldsAdded[9];
434  ***************************************************************************/
435 
436  Array<OneD, NekDouble> pressure(nSolutionPts, 0.0);
437  NekDouble gammaMinusOne = m_gamma - 1.0;
438 
439  for (i = 0; i < m_spacedim; i++)
440  {
441  Vmath::Vmul(nSolutionPts, &uFields[i + 1][0], 1, &uFields[i + 1][0], 1,
442  &tmp[0], 1);
443 
444  Vmath::Smul(nSolutionPts, 0.5, &tmp[0], 1, &tmp[0], 1);
445 
446  Vmath::Vadd(nSolutionPts, &pressure[0], 1, &tmp[0], 1, &pressure[0], 1);
447  }
448 
449  Vmath::Vdiv(nSolutionPts, &pressure[0], 1, &uFields[0][0], 1, &pressure[0],
450  1);
451 
452  Vmath::Vsub(nSolutionPts, &uFields[nfields - 1][0], 1, &pressure[0], 1,
453  &pressure[0], 1);
454 
455  Vmath::Smul(nSolutionPts, gammaMinusOne, &pressure[0], 1, &pressure[0], 1);
456 
457  // Extract trace
458  pFields[0]->ExtractTracePhys(pressure, traceFieldsAdded[9]);
459 
460  /******** Evaluation of the temperature ************************************
461  * T = P/(R*rho);
462  * T -> traceFieldsAdded[10];
463  ***************************************************************************/
464 
465  Array<OneD, NekDouble> temperature(nSolutionPts, 0.0);
466 
467  Vmath::Vdiv(nSolutionPts, &pressure[0], 1, &uFields[0][0], 1,
468  &temperature[0], 1);
469 
470  NekDouble GasConstantInv = 1.0 / m_gasConstant;
471  Vmath::Smul(nSolutionPts, GasConstantInv, &temperature[0], 1,
472  &temperature[0], 1);
473 
474  // Extract trace
475  pFields[0]->ExtractTracePhys(temperature, traceFieldsAdded[10]);
476 
477  /*** Evaluation of the temperature gradient in the normal direction ********
478  * DT_n -> traceFieldsAdded[11]
479  ***************************************************************************/
480 
481  Array<OneD, Array<OneD, NekDouble>> Dtemperature(nDimensions);
482  Array<OneD, Array<OneD, NekDouble>> traceDtemperature(nDimensions);
483 
484  for (i = 0; i < nDimensions; ++i)
485  {
486  Dtemperature[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
487  traceDtemperature[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
488  }
489 
490  for (i = 0; i < nDimensions; ++i)
491  {
492  for (n = 0; n < nElements; n++)
493  {
494  phys_offset = pFields[0]->GetPhys_Offset(n);
495 
496  pFields[i]->GetExp(n)->PhysDeriv(i, temperature + phys_offset,
497  auxArray =
498  Dtemperature[i] + phys_offset);
499  }
500  // Extract trace
501  pFields[0]->ExtractTracePhys(Dtemperature[i], traceDtemperature[i]);
502  }
503 
504  for (i = 0; i < nDimensions; ++i)
505  {
506  Vmath::Vmul(nTracePts, &m_traceNormals[i][0], 1,
507  &traceDtemperature[i][0], 1, &tmp[0], 1);
508 
509  Vmath::Vadd(nTracePts, &traceFieldsAdded[11][0], 1, &tmp[0], 1,
510  &traceFieldsAdded[11][0], 1);
511  }
512 
513  /*** Evaluation of the pressure gradient ***********************************
514  * DP_t -> traceFieldsAdded[12] tangent direction
515  * DP_b -> traceFieldsAdded[13] binormal direction
516  * DP_x -> traceFieldsAdded[14]
517  * DP_y -> traceFieldsAdded[15]
518  * DP_z -> traceFieldsAdded[16]
519  ***************************************************************************/
520 
521  Array<OneD, Array<OneD, NekDouble>> Dpressure(nDimensions);
522  Array<OneD, Array<OneD, NekDouble>> traceDpressure(nDimensions);
523 
524  for (i = 0; i < nDimensions; ++i)
525  {
526  Dpressure[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
527  traceDpressure[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
528  }
529 
530  for (i = 0; i < nDimensions; ++i)
531  {
532  for (n = 0; n < nElements; n++)
533  {
534  phys_offset = pFields[0]->GetPhys_Offset(n);
535 
536  pFields[i]->GetExp(n)->PhysDeriv(i, pressure + phys_offset,
537  auxArray =
538  Dpressure[i] + phys_offset);
539  }
540  // Extract trace
541  pFields[0]->ExtractTracePhys(Dpressure[i], traceDpressure[i]);
542  }
543 
544  // Dp_t
545  for (i = 0; i < nDimensions; ++i)
546  {
547  Vmath::Vmul(nTracePts, &m_traceTangents[i][0], 1, &traceDpressure[i][0],
548  1, &tmp[0], 1);
549 
550  Vmath::Vadd(nTracePts, &traceFieldsAdded[12][0], 1, &tmp[0], 1,
551  &traceFieldsAdded[12][0], 1);
552  }
553 
554  // Dp_b
555  for (i = 0; i < nDimensions; ++i)
556  {
557  Vmath::Vmul(nTracePts, &m_traceBinormals[i][0], 1,
558  &traceDpressure[i][0], 1, &tmp[0], 1);
559 
560  Vmath::Vadd(nTracePts, &traceFieldsAdded[13][0], 1, &tmp[0], 1,
561  &traceFieldsAdded[13][0], 1);
562  }
563 
564  // Dp_x
565  Vmath::Vcopy(nTracePts, &traceDpressure[0][0], 1, &traceFieldsAdded[14][0],
566  1);
567 
568  // Dp_y
569  Vmath::Vcopy(nTracePts, &traceDpressure[1][0], 1, &traceFieldsAdded[15][0],
570  1);
571 
572  // Dp_z
573  Vmath::Vcopy(nTracePts, &traceDpressure[2][0], 1, &traceFieldsAdded[16][0],
574  1);
575 
576  /** Evaluation of the velocity gradient in the cartesian directions
577  * Du_x: traceFieldsAdded[17]
578  * Du_y: traceFieldsAdded[18]
579  * Du_z: traceFieldsAdded[19]
580  * Dv_x: traceFieldsAdded[20]
581  * Dv_y: traceFieldsAdded[21]
582  * Dv_z: traceFieldsAdded[22]
583  * Dw_x: traceFieldsAdded[23]
584  * Dw_y: traceFieldsAdded[24]
585  * Dw_z: traceFieldsAdded[25]
586  **/
587 
588  Array<OneD, Array<OneD, Array<OneD, NekDouble>>> Dvelocity(nDimensions);
589  Array<OneD, Array<OneD, Array<OneD, NekDouble>>> traceDvelocity(
590  nDimensions);
591  Array<OneD, Array<OneD, NekDouble>> velocity(nDimensions);
592 
593  for (i = 0; i < nDimensions; ++i)
594  {
595  Dvelocity[i] = Array<OneD, Array<OneD, NekDouble>>(nDimensions);
596  traceDvelocity[i] = Array<OneD, Array<OneD, NekDouble>>(nDimensions);
597  velocity[i] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
598 
599  Vmath::Vdiv(nSolutionPts, uFields[i + 1], 1, uFields[0], 1, velocity[i],
600  1);
601 
602  for (j = 0; j < nDimensions; ++j)
603  {
604  Dvelocity[i][j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
605  traceDvelocity[i][j] = Array<OneD, NekDouble>(nTracePts, 0.0);
606  }
607  }
608 
609  for (i = 0; i < nDimensions; ++i)
610  {
611  for (j = 0; j < nDimensions; ++j)
612  {
613  for (n = 0; n < nElements; n++)
614  {
615  phys_offset = pFields[0]->GetPhys_Offset(n);
616 
617  pFields[i]->GetExp(n)->PhysDeriv(j, velocity[i] + phys_offset,
618  auxArray = Dvelocity[i][j] +
619  phys_offset);
620  }
621 
622  // Extract trace
623  pFields[0]->ExtractTracePhys(Dvelocity[i][j], traceDvelocity[i][j]);
624  }
625  }
626 
627  Vmath::Vcopy(nTracePts, &traceDvelocity[0][0][0], 1,
628  &traceFieldsAdded[17][0], 1);
629  Vmath::Vcopy(nTracePts, &traceDvelocity[0][1][0], 1,
630  &traceFieldsAdded[18][0], 1);
631  Vmath::Vcopy(nTracePts, &traceDvelocity[0][2][0], 1,
632  &traceFieldsAdded[19][0], 1);
633  Vmath::Vcopy(nTracePts, &traceDvelocity[1][0][0], 1,
634  &traceFieldsAdded[20][0], 1);
635  Vmath::Vcopy(nTracePts, &traceDvelocity[1][1][0], 1,
636  &traceFieldsAdded[21][0], 1);
637  Vmath::Vcopy(nTracePts, &traceDvelocity[1][2][0], 1,
638  &traceFieldsAdded[22][0], 1);
639  Vmath::Vcopy(nTracePts, &traceDvelocity[2][0][0], 1,
640  &traceFieldsAdded[23][0], 1);
641  Vmath::Vcopy(nTracePts, &traceDvelocity[2][1][0], 1,
642  &traceFieldsAdded[24][0], 1);
643  Vmath::Vcopy(nTracePts, &traceDvelocity[2][2][0], 1,
644  &traceFieldsAdded[25][0], 1);
645 
646  /*** Evaluation of shear stresses ******************************************
647  * tau_xx -> traceFieldsAdded[26]
648  * tau_yy -> traceFieldsAdded[27]
649  * tau_zz -> traceFieldsAdded[28]
650  * tau_xy -> traceFieldsAdded[29]
651  * tau_xz -> traceFieldsAdded[30]
652  * tau_yz -> traceFieldsAdded[31]
653  ***************************************************************************/
654 
655  // Stokes hypotesis
656  const NekDouble lambda = -2.0 / 3.0;
657 
658  // Auxiliary variables
659  Array<OneD, NekDouble> mu(nSolutionPts, 0.0);
660  Array<OneD, NekDouble> mu2(nSolutionPts, 0.0);
661  Array<OneD, NekDouble> divVel(nSolutionPts, 0.0);
662 
663  // Variable viscosity through the Sutherland's law
664  if (m_ViscosityType == "Variable")
665  {
666  NekDouble mu_star = m_mu;
667  NekDouble T_star = m_pInf / (m_rhoInf * m_gasConstant);
668  NekDouble ratio;
669 
670  for (int i = 0; i < nSolutionPts; ++i)
671  {
672  ratio = temperature[i] / T_star;
673  mu[i] = mu_star * ratio * sqrt(ratio) * (T_star + 110.0) /
674  (temperature[i] + 110.0);
675  }
676  }
677  else
678  {
679  Vmath::Fill(nSolutionPts, m_mu, &mu[0], 1);
680  }
681 
682  // Computing diagonal terms of viscous stress tensor
683  Array<OneD, Array<OneD, NekDouble>> temp(m_spacedim);
684  Array<OneD, Array<OneD, NekDouble>> Sgg(m_spacedim);
685 
686  // mu2 = 2 * mu
687  Vmath::Smul(nSolutionPts, 2.0, &mu[0], 1, &mu2[0], 1);
688 
689  // Velocity divergence
690  Vmath::Vadd(nSolutionPts, &divVel[0], 1, &Dvelocity[0][0][0], 1, &divVel[0],
691  1);
692  Vmath::Vadd(nSolutionPts, &divVel[0], 1, &Dvelocity[1][1][0], 1, &divVel[0],
693  1);
694 
695  // Velocity divergence scaled by lambda * mu
696  Vmath::Smul(nSolutionPts, lambda, &divVel[0], 1, &divVel[0], 1);
697  Vmath::Vmul(nSolutionPts, &mu[0], 1, &divVel[0], 1, &divVel[0], 1);
698 
699  // Diagonal terms of viscous stress tensor (Sxx, Syy)
700  // Sjj = 2 * mu * du_j/dx_j - (2 / 3) * mu * sum_j(du_j/dx_j)
701  for (j = 0; j < m_spacedim; ++j)
702  {
703  temp[j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
704  Sgg[j] = Array<OneD, NekDouble>(nSolutionPts, 0.0);
705 
706  Vmath::Vmul(nSolutionPts, &mu2[0], 1, &Dvelocity[j][j][0], 1,
707  &temp[j][0], 1);
708 
709  Vmath::Vadd(nSolutionPts, &temp[j][0], 1, &divVel[0], 1, &Sgg[j][0], 1);
710  }
711 
712  // Extra diagonal terms of viscous stress tensor
713  Array<OneD, NekDouble> Sxy(nSolutionPts, 0.0);
714  Array<OneD, NekDouble> Sxz(nSolutionPts, 0.0);
715  Array<OneD, NekDouble> Syz(nSolutionPts, 0.0);
716 
717  // Sxy = (du/dy + dv/dx)
718  Vmath::Vadd(nSolutionPts, &Dvelocity[0][1][0], 1, &Dvelocity[1][0][0], 1,
719  &Sxy[0], 1);
720 
721  // Sxz = (du/dz + dw/dx)
722  Vmath::Vadd(nSolutionPts, &Dvelocity[0][2][0], 1, &Dvelocity[2][0][0], 1,
723  &Sxz[0], 1);
724 
725  // Syz = (dv/dz + dw/dy)
726  Vmath::Vadd(nSolutionPts, &Dvelocity[1][2][0], 1, &Dvelocity[2][1][0], 1,
727  &Syz[0], 1);
728 
729  // Sxy = mu * (du/dy + dv/dx)
730  Vmath::Vmul(nSolutionPts, &mu[0], 1, &Sxy[0], 1, &Sxy[0], 1);
731 
732  // Sxz = mu * (du/dy + dv/dx)
733  Vmath::Vmul(nSolutionPts, &mu[0], 1, &Sxz[0], 1, &Sxz[0], 1);
734 
735  // Syz = mu * (du/dy + dv/dx)
736  Vmath::Vmul(nSolutionPts, &mu[0], 1, &Syz[0], 1, &Syz[0], 1);
737 
738  pFields[0]->ExtractTracePhys(Sgg[0], traceFieldsAdded[26]);
739  pFields[0]->ExtractTracePhys(Sgg[1], traceFieldsAdded[27]);
740  pFields[0]->ExtractTracePhys(Sgg[2], traceFieldsAdded[28]);
741  pFields[0]->ExtractTracePhys(Sxy, traceFieldsAdded[29]);
742  pFields[0]->ExtractTracePhys(Sxz, traceFieldsAdded[30]);
743  pFields[0]->ExtractTracePhys(Syz, traceFieldsAdded[31]);
744 
745  /*** Evaluation of dinamic viscosity ***************************************
746  * mu -> traceFieldsAdded[32]
747  ***************************************************************************/
748 
749  pFields[0]->ExtractTracePhys(mu, traceFieldsAdded[32]);
750 
751  /*** Evaluation of Mach number *********************************************
752  * M -> traceFieldsAdded[33]
753  ***************************************************************************/
754  NekDouble gamma = m_gamma;
755 
756  // Speed of sound
757  Array<OneD, NekDouble> soundspeed(nSolutionPts, 0.0);
758 
759  Vmath::Vdiv(nSolutionPts, pressure, 1, uFields[0], 1, soundspeed, 1);
760  Vmath::Smul(nSolutionPts, gamma, soundspeed, 1, soundspeed, 1);
761  Vmath::Vsqrt(nSolutionPts, soundspeed, 1, soundspeed, 1);
762 
763  // Mach
764  Array<OneD, NekDouble> mach(nSolutionPts, 0.0);
765 
766  for (int i = 0; i < m_spacedim; ++i)
767  {
768  Vmath::Vvtvp(nSolutionPts, uFields[i + 1], 1, uFields[i + 1], 1, mach,
769  1, mach, 1);
770  }
771 
772  Vmath::Vdiv(nSolutionPts, mach, 1, uFields[0], 1, mach, 1);
773  Vmath::Vdiv(nSolutionPts, mach, 1, uFields[0], 1, mach, 1);
774  Vmath::Vsqrt(nSolutionPts, mach, 1, mach, 1);
775  Vmath::Vdiv(nSolutionPts, mach, 1, soundspeed, 1, mach, 1);
776 
777  pFields[0]->ExtractTracePhys(mach, traceFieldsAdded[33]);
778 
779  /**************************************************************************/
780  // Extract coordinates
781 
782  if (pFields[0]->GetBndCondExpansions().size())
783  {
784  id1 = 0;
785  cnt = 0;
786  nBndRegions = pFields[0]->GetBndCondExpansions().size();
787  for (b = 0; b < nBndRegions; ++b)
788  {
789  nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
790  for (e = 0; e < nBndEdges; ++e)
791  {
792  nBndEdgePts = pFields[0]
793  ->GetBndCondExpansions()[b]
794  ->GetExp(e)
795  ->GetTotPoints();
796 
797  id2 = pFields[0]->GetTrace()->GetPhys_Offset(
798  pFields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
799  cnt++));
800 
801  if (pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
802  "WallViscous" ||
803  pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
804  "WallAdiabatic" ||
805  pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
806  "Wall")
807  {
808 
809  Vmath::Vcopy(nBndEdgePts, &traceX[id2], 1, &surfaceX[id1],
810  1);
811 
812  Vmath::Vcopy(nBndEdgePts, &traceY[id2], 1, &surfaceY[id1],
813  1);
814 
815  Vmath::Vcopy(nBndEdgePts, &traceZ[id2], 1, &surfaceZ[id1],
816  1);
817 
818  id1 += nBndEdgePts;
819  }
820  }
821  }
822  }
823 
824  // Extract fields
825  if (pFields[0]->GetBndCondExpansions().size())
826  {
827 
828  for (j = 0; j < nfields; ++j)
829  {
830  cout << "field " << j << endl;
831 
832  id1 = 0;
833  cnt = 0;
834  nBndRegions = pFields[j]->GetBndCondExpansions().size();
835  for (b = 0; b < nBndRegions; ++b)
836  {
837  nBndEdges = pFields[j]->GetBndCondExpansions()[b]->GetExpSize();
838  for (e = 0; e < nBndEdges; ++e)
839  {
840  nBndEdgePts = pFields[j]
841  ->GetBndCondExpansions()[b]
842  ->GetExp(e)
843  ->GetTotPoints();
844 
845  id2 = pFields[j]->GetTrace()->GetPhys_Offset(
846  pFields[j]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
847  cnt++));
848 
849  if (pFields[j]->GetBndConditions()[b]->GetUserDefined() ==
850  "WallViscous" ||
851  pFields[j]->GetBndConditions()[b]->GetUserDefined() ==
852  "WallAdiabatic" ||
853  pFields[j]->GetBndConditions()[b]->GetUserDefined() ==
854  "Wall")
855  {
856  Vmath::Vcopy(nBndEdgePts, &traceFields[j][id2], 1,
857  &surfaceFields[j][id1], 1);
858 
859  id1 += nBndEdgePts;
860  }
861  }
862  }
863  }
864  }
865 
866  // Extract fields added
867  if (pFields[0]->GetBndCondExpansions().size())
868  {
869  for (j = 0; j < nfieldsAdded; ++j)
870  {
871  cout << "field added " << j << endl;
872 
873  id1 = 0;
874  cnt = 0;
875  nBndRegions = pFields[0]->GetBndCondExpansions().size();
876  for (b = 0; b < nBndRegions; ++b)
877  {
878  nBndEdges = pFields[0]->GetBndCondExpansions()[b]->GetExpSize();
879  for (e = 0; e < nBndEdges; ++e)
880  {
881  nBndEdgePts = pFields[0]
882  ->GetBndCondExpansions()[b]
883  ->GetExp(e)
884  ->GetTotPoints();
885 
886  id2 = pFields[0]->GetTrace()->GetPhys_Offset(
887  pFields[0]->GetTraceMap()->GetBndCondIDToGlobalTraceID(
888  cnt++));
889 
890  if (pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
891  "WallViscous" ||
892  pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
893  "WallAdiabatic" ||
894  pFields[0]->GetBndConditions()[b]->GetUserDefined() ==
895  "Wall")
896  {
897  Vmath::Vcopy(nBndEdgePts, &traceFieldsAdded[j][id2], 1,
898  &surfaceFieldsAdded[j][id1], 1);
899 
900  id1 += nBndEdgePts;
901  }
902  }
903  }
904  }
905  }
906 
907  //==========================================================================
908  //==========================================================================
909  //==========================================================================
910 
911  // Print the surface coordinates and the surface solution in a .txt file
912  ofstream outfile;
913  outfile.open(fname.c_str());
914  outfile << "% x[m] "
915  << " \t"
916  << "y[m] "
917  << " \t"
918  << "z[m] "
919  << " \t"
920  << "nx[] "
921  << " \t"
922  << "ny[] "
923  << " \t"
924  << "nz[] "
925  << " \t"
926  << "bx[] "
927  << " \t"
928  << "by[] "
929  << " \t"
930  << "bz[] "
931  << " \t"
932  << "tx[] "
933  << " \t"
934  << "ty[] "
935  << " \t"
936  << "tz[] "
937  << " \t"
938  << "rho[kg/m^3] "
939  << " \t"
940  << "rhou[kg/(m^2 s)] "
941  << " \t"
942  << "rhov[kg/(m^2 s)] "
943  << " \t"
944  << "rhow[kg/(m^2 s)] "
945  << " \t"
946  << "E[Pa] "
947  << " \t"
948  << "p[Pa] "
949  << " \t"
950  << "T[k] "
951  << " \t"
952  << "dT/dn[k/m] "
953  << " \t"
954  << "dp/dT[Pa/m] "
955  << " \t"
956  << "dp/dB[Pa/m] "
957  << " \t"
958  << "dp/dx[Pa/m] "
959  << " \t"
960  << "dp/dy[Pa/m] "
961  << " \t"
962  << "dp/dz[Pa/m] "
963  << " \t"
964  << "du/dx[s^-1] "
965  << " \t"
966  << "du/dy[s^-1] "
967  << " \t"
968  << "du/dz[s^-1] "
969  << " \t"
970  << "dv/dx[s^-1] "
971  << " \t"
972  << "dv/dy[s^-1] "
973  << " \t"
974  << "dv/dz[s^-1] "
975  << " \t"
976  << "dw/dx[s^-1] "
977  << " \t"
978  << "dw/dy[s^-1] "
979  << " \t"
980  << "dw/dz[s^-1] "
981  << " \t"
982  << "tau_xx[Pa] "
983  << " \t"
984  << "tau_yy[Pa] "
985  << " \t"
986  << "tau_zz[Pa] "
987  << " \t"
988  << "tau_xy[Pa] "
989  << " \t"
990  << "tau_xz[Pa] "
991  << " \t"
992  << "tau_yz[Pa] "
993  << " \t"
994  << "mu[Pa s] "
995  << " \t"
996  << "M[] "
997  << " \t" << endl;
998  for (i = 0; i < nSurfacePts; ++i)
999  {
1000  outfile
1001  << scientific << setw(17) << setprecision(16) << surfaceX[i]
1002  << " \t " << surfaceY[i] << " \t " << surfaceZ[i] << " \t "
1003  << surfaceFieldsAdded[0][i] << " \t " << surfaceFieldsAdded[1][i]
1004  << " \t " << surfaceFieldsAdded[2][i] << " \t "
1005  << surfaceFieldsAdded[3][i] << " \t " << surfaceFieldsAdded[4][i]
1006  << " \t " << surfaceFieldsAdded[5][i] << " \t "
1007  << surfaceFieldsAdded[6][i] << " \t " << surfaceFieldsAdded[7][i]
1008  << " \t " << surfaceFieldsAdded[8][i] << " \t "
1009  << surfaceFields[0][i] << " \t " << surfaceFields[1][i] << " \t "
1010  << surfaceFields[2][i] << " \t " << surfaceFields[3][i] << " \t "
1011  << surfaceFields[4][i] << " \t " << surfaceFieldsAdded[9][i]
1012  << " \t " << surfaceFieldsAdded[10][i] << " \t "
1013  << surfaceFieldsAdded[11][i] << " \t " << surfaceFieldsAdded[12][i]
1014  << " \t " << surfaceFieldsAdded[13][i] << " \t "
1015  << surfaceFieldsAdded[14][i] << " \t " << surfaceFieldsAdded[15][i]
1016  << " \t " << surfaceFieldsAdded[16][i] << " \t "
1017  << surfaceFieldsAdded[17][i] << " \t " << surfaceFieldsAdded[18][i]
1018  << " \t " << surfaceFieldsAdded[19][i] << " \t "
1019  << surfaceFieldsAdded[20][i] << " \t " << surfaceFieldsAdded[21][i]
1020  << " \t " << surfaceFieldsAdded[22][i] << " \t "
1021  << surfaceFieldsAdded[23][i] << " \t " << surfaceFieldsAdded[24][i]
1022  << " \t " << surfaceFieldsAdded[25][i] << " \t "
1023  << surfaceFieldsAdded[26][i] << " \t " << surfaceFieldsAdded[27][i]
1024  << " \t " << surfaceFieldsAdded[28][i] << " \t "
1025  << surfaceFieldsAdded[29][i] << " \t " << surfaceFieldsAdded[30][i]
1026  << " \t " << surfaceFieldsAdded[31][i] << " \t "
1027  << surfaceFieldsAdded[32][i] << " \t " << surfaceFieldsAdded[33][i]
1028  << " \t " << endl;
1029  }
1030  outfile << endl << endl;
1031  outfile.close();
1032 
1033  return 0;
1034 }
m_mu
NekDouble m_mu
Definition: CompressibleBL.cpp:86
m_Twall
NekDouble m_Twall
Definition: CompressibleBL.cpp:78
m_rhoInf
NekDouble m_rhoInf
Definition: CompressibleBL.cpp:83
m_vInf
NekDouble m_vInf
Definition: CompressibleBL.cpp:85
m_uInf
NekDouble m_uInf
Definition: CompressibleBL.cpp:82
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:83
CG_Iterations.pressure
pressure
Definition: CG_Iterations.py:229
Nektar::LibUtilities::Import
void Import(const std::string &infilename, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble >> &fielddata, FieldMetaDataMap &fieldinfomap, const Array< OneD, int > &ElementIDs)
This function allows for data to be imported from an FLD file when a session and/or communicator is n...
Definition: FieldIO.cpp:290
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:115
Nektar::LibUtilities::ePolyEvenlySpaced
@ ePolyEvenlySpaced
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:75
Nektar::MultiRegions::ExpListSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:55
Nektar::SpatialDomains::MeshGraphSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:172
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vsqrt
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:534
Vmath::Vmul
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
Vmath::Vvtvp
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
Vmath::Vadd
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
Vmath::Smul
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
Vmath::Vdiv
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
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
Vmath::Vsub
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
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294

References ASSERTL0, Nektar::LibUtilities::ePolyEvenlySpaced, Vmath::Fill(), Nektar::LibUtilities::Import(), m_mu, m_rhoInf, m_Twall, m_uInf, m_vInf, CG_Iterations::pressure, Vmath::Smul(), tinysimd::sqrt(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vdiv(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vsub(), and Vmath::Vvtvp().

Generated on Sat Mar 25 2023 12:35:53 for Nektar++ by doxygen 1.9.1