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