Nektar++
Functions
ExtractSurface3DCFS.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <string>
#include <iostream>
#include <iomanip>
#include <MultiRegions/ExpList.h>
#include <MultiRegions/ExpList1D.h>
#include <MultiRegions/ExpList2D.h>
#include <MultiRegions/ExpList3D.h>
#include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
#include <MultiRegions/DisContField3D.h>
#include <LocalRegions/MatrixKey.h>
#include <LocalRegions/Expansion3D.h>
#include <LocalRegions/Expansion.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/ContField3D.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 68 of file ExtractSurface3DCFS.cpp.

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(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vdiv(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vsub(), and Vmath::Vvtvp().

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