Nektar++
Functions
ExtractSurface2DCFS.cpp File Reference
#include <cstdio>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <string>
#include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
#include <MultiRegions/ExpList.h>
#include <MultiRegions/ExpList2DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous2D.h>
#include <LocalRegions/Expansion.h>
#include <LocalRegions/Expansion2D.h>
#include <LocalRegions/MatrixKey.h>
#include <MultiRegions/DisContField.h>
#include <LibUtilities/BasicUtils/ErrorUtil.hpp>
#include <LibUtilities/BasicUtils/FieldIO.h>
#include <LibUtilities/BasicUtils/NekFactory.hpp>
#include <LibUtilities/BasicUtils/SessionReader.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/Communication/Comm.h>
#include <LibUtilities/Memory/NekMemoryManager.hpp>
#include <MultiRegions/ContField.h>
#include <SpatialDomains/MeshGraph.h>
#include <SolverUtils/SolverUtilsDeclspec.h>
#include <boost/program_options.hpp>

Go to the source code of this file.

Functions

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

Function Documentation

◆ main()

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

Evaluation of the velocity gradient in the cartesian directions Du_x: traceFieldsAdded[10] Du_y: traceFieldsAdded[11] Dv_x: traceFieldsAdded[12] Dv_y: traceFieldsAdded[13]

Definition at line 73 of file ExtractSurface2DCFS.cpp.

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

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(), CellMLToNektar.cellml_metadata::p, CG_Iterations::pressure, CellMLToNektar.translators::run(), Vmath::Smul(), tinysimd::sqrt(), Vmath::Vadd(), Vmath::Vcopy(), Vmath::Vdiv(), Vmath::Vmul(), Vmath::Vsqrt(), Vmath::Vsub(), and Vmath::Vvtvp().

Generated on Sat Mar 25 2023 12:35:53 for Nektar++ by doxygen 1.9.1