Nektar++
Functions
ExtractSurface2DCFS.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <string>
#include <iostream>
#include <iomanip>
#include <MultiRegions/ExpList.h>
#include <MultiRegions/ExpList2DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous2D.h>
#include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
#include <MultiRegions/DisContField.h>
#include <LocalRegions/MatrixKey.h>
#include <LocalRegions/Expansion2D.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/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[10] Du_y: traceFieldsAdded[11] Dv_x: traceFieldsAdded[12] Dv_y: traceFieldsAdded[13]

Definition at line 68 of file ExtractSurface2DCFS.cpp.

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

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

Generated on Wed Nov 24 2021 10:33:29 for Nektar++ by doxygen 1.9.1