Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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/MeshGraph2D.h>
#include <SolverUtils/SolverUtilsDeclspec.h>
Include dependency graph for ExtractSurface2DCFS.cpp:

Go to the source code of this file.

Functions

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

Function Documentation

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

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