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