Nektar++
Functions
ExtractSurface2DCFS.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/ExpList2DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous2D.h>
#include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
#include <MultiRegions/DisContField2D.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/ContField2D.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 70 of file ExtractSurface2DCFS.cpp.

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

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