Nektar++
CompressibleBL.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: CompressibleBL.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: Generate the compressible boundary layer similarity solution
32 // on a 2D/3D mesh.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <cmath>
37 #include <cstdio>
38 #include <cstdlib>
39 #include <iomanip>
40 #include <iostream>
41 #include <string>
42 
43 #include <boost/core/ignore_unused.hpp>
44 
45 #include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
46 #include <MultiRegions/ContField.h>
47 #include <MultiRegions/DisContField.h>
48 #include <MultiRegions/ExpList.h>
49 #include <MultiRegions/ExpList2DHomogeneous1D.h>
50 #include <MultiRegions/ExpList3DHomogeneous1D.h>
51 #include <MultiRegions/ExpList3DHomogeneous2D.h>
52 
53 #include <LocalRegions/Expansion.h>
54 #include <LocalRegions/Expansion2D.h>
55 #include <LocalRegions/Expansion3D.h>
56 #include <LocalRegions/MatrixKey.h>
57 
58 #include <LibUtilities/BasicUtils/FieldIO.h>
59 #include <LibUtilities/BasicUtils/NekFactory.hpp>
60 #include <LibUtilities/BasicUtils/SessionReader.h>
61 #include <LibUtilities/BasicUtils/SharedArray.hpp>
62 #include <LibUtilities/Communication/Comm.h>
63 #include <LibUtilities/Memory/NekMemoryManager.hpp>
64 
65 #include <SpatialDomains/MeshGraph.h>
66 
67 #include <SolverUtils/SolverUtilsDeclspec.h>
68 
69 using namespace std;
70 using namespace Nektar;
71 
72 NekDouble m_Re;
73 NekDouble m_Mach;
74 NekDouble L;
75 NekDouble m_Tinf;
76 NekDouble m_Suth;
77 NekDouble m_Tw;
78 NekDouble m_Twall;
79 NekDouble m_Gamma;
80 NekDouble m_Pr;
81 NekDouble m_long;
82 NekDouble m_uInf;
83 NekDouble m_rhoInf;
84 NekDouble m_R;
85 NekDouble m_vInf;
86 NekDouble m_mu;
87 NekDouble m_To = 273.11;
88 
89 const int m_xpoints = 1000001;
90 
91 const NekDouble Nvisc = 1;
92 const NekDouble Omega = 1;
93 const NekDouble etamax = 10.0;
94 const NekDouble errtol = 1e-5;
95 
96 /**
97  * Calculate the compressible boundary layer using the similarity solution
98  */
99 void COMPBL(Array<OneD, NekDouble> v, Array<OneD, NekDouble> dv)
100 {
101  NekDouble c, dcdg, cp;
102 
103  if (Nvisc == 1)
104  {
105  c = sqrt(v[3]) * (1.0 + m_Suth) / (v[3] + m_Suth);
106  dcdg = 1.0 / (2. * sqrt(v[3])) - sqrt(v[3]) / (v[3] + m_Suth);
107  dcdg = dcdg * (1.0 + m_Suth) / (v[3] + m_Suth);
108  cp = dcdg * v[4];
109  }
110  if (Nvisc == 2)
111  {
112  c = pow(v[3], (Omega - 1.0));
113  dcdg = (Omega - 1.0) * pow(v[3], (Omega - 2.0));
114  cp = dcdg * v[4];
115  }
116  if (Nvisc == 3)
117  {
118  c = sqrt(m_Twall) * (1.0 + m_Suth) / (m_Suth + m_Twall);
119  cp = 0.0;
120  }
121 
122  dv[0] = v[1];
123  dv[1] = v[2];
124  dv[2] = -v[2] * (cp + v[0]) / c;
125  dv[3] = v[4];
126  dv[4] = -v[4] * (cp + m_Pr * v[0]) / c -
127  m_Pr * (m_Gamma - 1.0) * pow(m_Mach, 2.0) * pow(v[2], 2);
128 }
129 
130 /**
131  * Perform the RK4 integration
132  */
133 void RK4(Array<OneD, NekDouble> y, Array<OneD, NekDouble> dydx, int n,
134  NekDouble x, NekDouble h, Array<OneD, NekDouble> yout)
135 {
136  boost::ignore_unused(x);
137 
138  int nmax = 5;
139 
140  Array<OneD, NekDouble> yt(nmax, 0.0);
141  Array<OneD, NekDouble> dyt(nmax, 0.0);
142  Array<OneD, NekDouble> dym(nmax, 0.0);
143  NekDouble hh = h * 0.5;
144  NekDouble h6 = h / 6;
145 
146  for (int i = 0; i < n; i++)
147  {
148  yt[i] = y[i] + hh * dydx[i];
149  }
150 
151  COMPBL(yt, dyt);
152 
153  for (int i = 0; i < n; i++)
154  {
155  yt[i] = y[i] + hh * dyt[i];
156  }
157 
158  COMPBL(yt, dym);
159 
160  for (int i = 0; i < n; i++)
161  {
162  yt[i] = y[i] + h * dym[i];
163  dym[i] = dyt[i] + dym[i];
164  }
165 
166  COMPBL(yt, dyt);
167 
168  for (int i = 0; i < n; i++)
169  {
170  yout[i] = y[i] + h6 * (dydx[i] + dyt[i] + 2 * dym[i]);
171  }
172 }
173 
174 /**
175  * Calculate initial guess for RK4
176  */
177 void RKDUMB(Array<OneD, NekDouble> vstart, int nvar, NekDouble x1, NekDouble x2,
178  int m_xpoints, Array<OneD, NekDouble> xx,
179  Array<OneD, Array<OneD, NekDouble>> y)
180 {
181  int nmax = 5;
182  NekDouble x, h;
183  Array<OneD, NekDouble> v(nmax, 0.0);
184  Array<OneD, NekDouble> dv(nmax, 0.0);
185 
186  for (int i = 0; i < nvar; i++)
187  {
188  v[i] = vstart[i];
189  y[i][0] = v[i];
190  }
191 
192  xx[0] = x1;
193  x = x1;
194  h = (x2 - x1) / m_xpoints;
195 
196  for (int k = 0; k < m_xpoints; k++)
197  {
198  COMPBL(v, dv);
199  RK4(v, dv, nvar, x, h, v);
200 
201  if (x + h == x)
202  {
203  cout << "bug" << endl;
204  }
205 
206  x = x + h;
207  xx[k + 1] = x;
208 
209  for (int i = 0; i < nvar; i++)
210  {
211  y[i][k + 1] = v[i];
212  }
213  }
214 }
215 
216 /**
217  * Create the output file
218  */
219 void OUTPUT(int m_xpoints, Array<OneD, NekDouble> xx,
220  Array<OneD, Array<OneD, NekDouble>> ff, int nQuadraturePts,
221  Array<OneD, NekDouble> x_QuadraturePts,
222  Array<OneD, NekDouble> y_QuadraturePts,
223  Array<OneD, NekDouble> u_QuadraturePts,
224  Array<OneD, NekDouble> v_QuadraturePts,
225  Array<OneD, NekDouble> rho_QuadraturePts,
226  Array<OneD, NekDouble> T_QuadraturePts)
227 {
228  Array<OneD, NekDouble> z(m_xpoints, 0.0);
229  Array<OneD, NekDouble> v(m_xpoints, 0.0);
230  Array<OneD, NekDouble> dv(m_xpoints, 0.0);
231  Array<OneD, NekDouble> u(m_xpoints, 0.0);
232  Array<OneD, NekDouble> t(m_xpoints, 0.0);
233  Array<OneD, NekDouble> rho(m_xpoints, 0.0);
234  Array<OneD, NekDouble> mu(m_xpoints, 0.0);
235  Array<OneD, NekDouble> vv(m_xpoints, 0.0);
236  Array<OneD, NekDouble> velocity(m_xpoints, 0.0);
237  Array<OneD, NekDouble> test(m_xpoints, 0.0);
238 
239  NekDouble dd, dm, scale;
240  NekDouble xcher, ycher;
241  int index = -1;
242 
243  z[0] = 0.0;
244  NekDouble sumd = 0.0;
245 
246  for (int i = 1; i < m_xpoints; i++)
247  {
248  z[i] = z[i - 1] + 0.5 * (xx[i] - xx[i - 1]) * (ff[3][i] + ff[3][i - 1]);
249  dm = ff[3][i - 1] - ff[1][i - 1];
250  dd = ff[3][i] - ff[1][i];
251  sumd = sumd + 0.5 * (xx[i] - xx[i - 1]) * (dd + dm);
252  }
253 
254  scale = sumd;
255 
256  ofstream file3;
257  file3.open("physical_data.dat");
258 
259  NekDouble xin, rex, delsx, delta;
260 
261  for (int i = 0; i < m_xpoints; i++)
262  {
263  for (int k = 0; k < 5; k++)
264  {
265  v[k] = ff[k][i];
266  }
267  COMPBL(v, dv);
268  u[i] = ff[1][i];
269  t[i] = ff[3][i];
270  rho[i] = (1.0 / ff[3][i]);
271  vv[i] = -ff[0][i] / sqrt(m_uInf);
272  mu[i] = pow(t[i], 1.5) * (1 + m_Suth) / (t[i] + m_Suth) / (m_Re);
273  velocity[i] = ff[0][i];
274  }
275 
276  NekDouble scale2, coeff;
277 
278  for (int i = 0; i < nQuadraturePts; i++)
279  {
280  if (i % 100000 == 0)
281  {
282  cout << "i"
283  << " " << i << "/" << nQuadraturePts << endl;
284  }
285 
286  xcher = x_QuadraturePts[i];
287  ycher = y_QuadraturePts[i];
288 
289  scale = sumd;
290  xin = xcher;
291  rex = 0.5 * pow(((m_Re) / scale), 2) + (m_Re)*xin;
292  delsx = sqrt(2.0 / rex) * scale * (xin)*m_Pr;
293  scale = scale / delsx;
294  delta = 4.91 * sqrt((xin * m_mu) / (m_rhoInf * m_uInf));
295  scale2 = ycher * (scale * delta) / sqrt(etamax);
296  coeff = 0.5 * sqrt(2 / (xcher * m_Re));
297 
298  if (scale2 > z[m_xpoints - 3])
299  {
300  u_QuadraturePts[i] = 1;
301  rho_QuadraturePts[i] = 1;
302  T_QuadraturePts[i] = 1.0 / rho_QuadraturePts[i];
303  v_QuadraturePts[i] =
304  coeff * (z[m_xpoints - 3] - velocity[m_xpoints - 3]);
305 
306  file3 << xcher << " " << ycher << " "
307  << velocity[m_xpoints - 3] << " " << z[m_xpoints - 3]
308  << " " << u[m_xpoints - 3] << endl;
309  }
310  else
311  {
312  for (int j = 0; j < m_xpoints - 1; j++)
313  {
314  if ((z[j] <= scale2) && (z[j + 1] > scale2))
315  {
316  index = j;
317  break;
318  }
319  }
320  if (index == -1)
321  {
322  ASSERTL0(false, "Could not determine index in CompressibleBL");
323  }
324 
325  u_QuadraturePts[i] = u[index];
326  rho_QuadraturePts[i] = rho[index];
327  T_QuadraturePts[i] = 1.0 / rho_QuadraturePts[i];
328  v_QuadraturePts[i] = coeff * (u[index] * scale2 - velocity[index]);
329  }
330  }
331 }
332 
333 /**
334  * Calculate the compressible boundary layer solution for a given 3D mesh
335  * and dump the solution into a .rst file.
336  */
337 int main(int argc, char *argv[])
338 {
339  // Variable initialisation
340  int nmax = 5;
341  int maxit = 10;
342 
343  string opt;
344 
345  int i, j, k, numModes;
346 
347  Array<OneD, NekDouble> xx(m_xpoints, 0.0);
348  Array<OneD, Array<OneD, NekDouble>> ff(nmax);
349  Array<OneD, NekDouble> parameter(9, 0.0);
350 
351  for (i = 0; i < nmax; i++)
352  {
353  ff[i] = Array<OneD, NekDouble>(m_xpoints);
354  }
355 
356  Array<OneD, NekDouble> vstart(nmax, 0.0);
357  Array<OneD, NekDouble> v(2);
358  Array<OneD, NekDouble> dv(2);
359  Array<OneD, NekDouble> f(2);
360  Array<OneD, NekDouble> f1(2);
361  Array<OneD, NekDouble> f2(2);
362 
363  NekDouble al11, al21, al12, al22, det;
364 
365  // Reading the session file
366  LibUtilities::SessionReaderSharedPtr vSession =
367  LibUtilities::SessionReader::CreateInstance(argc, argv);
368 
369  // Read in mesh from input file and create an object of class MeshGraph
370  SpatialDomains::MeshGraphSharedPtr graphShPt =
371  SpatialDomains::MeshGraph::Read(vSession);
372 
373  int expdim = graphShPt->GetMeshDimension();
374 
375  int nElements, nQuadraturePts = 0;
376  Array<OneD, NekDouble> x_QuadraturePts;
377  Array<OneD, NekDouble> y_QuadraturePts;
378  Array<OneD, NekDouble> z_QuadraturePts;
379 
380  MultiRegions::ContFieldSharedPtr Domain;
381  Domain = MemoryManager<MultiRegions::ContField>::AllocateSharedPtr(
382  vSession, graphShPt, vSession->GetVariable(0));
383 
384  // Get the total number of elements
385  nElements = Domain->GetExpSize();
386  std::cout << "Number of elements = " << nElements
387  << std::endl;
388 
389  // Get the total number of quadrature points (depends on n. modes)
390  nQuadraturePts = Domain->GetTotPoints();
391  std::cout << "Number of quadrature points = " << nQuadraturePts
392  << std::endl;
393 
394  // Coordinates of the quadrature points
395  x_QuadraturePts = Array<OneD, NekDouble>(nQuadraturePts);
396  y_QuadraturePts = Array<OneD, NekDouble>(nQuadraturePts);
397  z_QuadraturePts = Array<OneD, NekDouble>(nQuadraturePts);
398  Domain->GetCoords(x_QuadraturePts, y_QuadraturePts, z_QuadraturePts);
399 
400  if (expdim == 1)
401  {
402  ASSERTL0(false, "Routine available for 2D and 3D problem only.")
403  }
404 
405  // Loading parameters from session file
406  vSession->LoadParameter("Re", m_Re, 1.0);
407  vSession->LoadParameter("Mach", m_Mach, 1.0);
408  vSession->LoadParameter("TInf", m_Tinf, 1.0);
409  vSession->LoadParameter("Twall", m_Twall, 1.0);
410  vSession->LoadParameter("Gamma", m_Gamma, 1.0);
411  vSession->LoadParameter("Pr", m_Pr, 1.0);
412  vSession->LoadParameter("L", m_long, 1.0);
413  vSession->LoadParameter("rhoInf", m_rhoInf, 1.0);
414  vSession->LoadParameter("uInf", m_uInf, 1.0);
415  vSession->LoadParameter("GasConstant", m_R, 1.0);
416  vSession->LoadParameter("vInf", m_vInf, 1.0);
417  vSession->LoadParameter("mu", m_mu, 1.0);
418 
419  // Rescaling factors
420  m_Suth = 110.4 / m_Tinf;
421  m_Tw = m_Twall / m_Tinf;
422  m_Re = m_Re / m_long;
423 
424  cout << "Number of points"
425  << " " << m_xpoints << endl;
426 
427  // Defining the solution arrays
428  Array<OneD, NekDouble> u_QuadraturePts(nQuadraturePts, 0.0);
429  Array<OneD, NekDouble> v_QuadraturePts(nQuadraturePts, 0.0);
430  Array<OneD, NekDouble> rho_QuadraturePts(nQuadraturePts, 0.0);
431  Array<OneD, NekDouble> T_QuadraturePts(nQuadraturePts, 0.0);
432 
433  // Calculation of the similarity variables
434  if (m_Tw > 0)
435  {
436  vstart[3] = m_Tw;
437  }
438  if (m_Tw < 0.0)
439  {
440  v[1] = 1.0 + 0.5 * 0.84 * (m_Gamma - 1) * (m_Mach * m_Mach);
441  v[0] = 0.47 * pow(v[1], 0.21);
442  }
443  else
444  {
445  v[1] = 0.062 * pow(m_Mach, 2) -
446  0.1 * (m_Tw - 1.0) * (10 + m_Mach) / (0.2 + m_Mach);
447  v[0] = 0.45 - 0.01 * m_Mach + (m_Tw - 1.0) * 0.06;
448  m_Twall = m_Tw;
449  }
450 
451  dv[0] = v[0] * 0.01;
452 
453  if (m_Tw < 0.0)
454  {
455  dv[1] = v[1] * 0.01;
456  }
457  else
458  {
459  dv[1] = 0.1;
460  }
461 
462  vstart[2] = v[0];
463 
464  if (m_Tw < 0)
465  {
466  vstart[3] = v[1];
467  m_Twall = vstart[3];
468  }
469  else
470  {
471  vstart[4] = v[1];
472  }
473 
474  RKDUMB(vstart, 5, 0.0, etamax, m_xpoints, xx, ff);
475 
476  for (k = 0; k < maxit; k++)
477  {
478  vstart[2] = v[0];
479 
480  if (m_Tw < 0)
481  {
482  vstart[3] = v[1];
483  m_Twall = vstart[3];
484  }
485  else
486  {
487  vstart[4] = v[1];
488  }
489 
490  RKDUMB(vstart, 5, 0.0, etamax, m_xpoints, xx, ff);
491 
492  NekDouble err = fabs(ff[1][m_xpoints] - 1) + fabs(ff[3][m_xpoints] - 1);
493 
494  cout << "err" << scientific << setprecision(9) << " " << err << endl;
495 
496  if (expdim == 2)
497  {
498  if (err < errtol)
499  {
500  cout << "Calculating" << endl;
501  OUTPUT(m_xpoints, xx, ff, nQuadraturePts, x_QuadraturePts,
502  y_QuadraturePts, u_QuadraturePts, v_QuadraturePts,
503  rho_QuadraturePts, T_QuadraturePts);
504  break;
505  }
506  else
507  {
508  f[0] = ff[1][m_xpoints] - 1;
509  f[1] = ff[3][m_xpoints] - 1;
510  vstart[2] = v[0] + dv[0];
511 
512  if (m_Tw < 0)
513  {
514  vstart[3] = v[1];
515  m_Twall = vstart[3];
516  }
517  else
518  {
519  vstart[4] = v[1];
520  }
521 
522  RKDUMB(vstart, 5, 0.0, etamax, m_xpoints, xx, ff);
523 
524  f1[0] = ff[1][m_xpoints] - 1;
525  f1[1] = ff[3][m_xpoints] - 1;
526 
527  vstart[2] = v[0];
528 
529  if (m_Tw < 0)
530  {
531  vstart[3] = v[1] + dv[1];
532  m_Twall = vstart[3];
533  }
534  else
535  {
536  vstart[4] = v[1] + dv[1];
537  }
538 
539  RKDUMB(vstart, 5, 0.0, etamax, m_xpoints, xx, ff);
540 
541  f2[0] = ff[1][m_xpoints] - 1;
542  f2[1] = ff[3][m_xpoints] - 1;
543 
544  al11 = (f1[0] - f[0]) / dv[0];
545  al21 = (f1[1] - f[1]) / dv[0];
546  al12 = (f2[0] - f[0]) / dv[1];
547  al22 = (f2[1] - f[1]) / dv[1];
548  det = al11 * al22 - al21 * al12;
549 
550  dv[0] = (-al22 * f[0] + al12 * f[1]) / det;
551  dv[1] = (al21 * f[0] - al11 * f[1]) / det;
552  v[0] = v[0] + dv[0];
553  v[1] = v[1] + dv[1];
554  }
555  }
556  else if (expdim == 3)
557  {
558  {
559  if (err < errtol)
560  {
561  cout << "Calculating" << endl;
562  OUTPUT(m_xpoints, xx, ff, nQuadraturePts, x_QuadraturePts,
563  z_QuadraturePts, u_QuadraturePts, v_QuadraturePts,
564  rho_QuadraturePts, T_QuadraturePts);
565  break;
566  }
567  else
568  {
569  f[0] = ff[1][m_xpoints] - 1;
570  f[1] = ff[3][m_xpoints] - 1;
571  vstart[2] = v[0] + dv[0];
572 
573  if (m_Tw < 0)
574  {
575  vstart[3] = v[1];
576  m_Twall = vstart[3];
577  }
578  else
579  {
580  vstart[4] = v[1];
581  }
582 
583  RKDUMB(vstart, 5, 0.0, etamax, m_xpoints, xx, ff);
584 
585  f1[0] = ff[1][m_xpoints] - 1;
586  f1[1] = ff[3][m_xpoints] - 1;
587 
588  vstart[2] = v[0];
589 
590  if (m_Tw < 0)
591  {
592  vstart[3] = v[1] + dv[1];
593  m_Twall = vstart[3];
594  }
595  else
596  {
597  vstart[4] = v[1] + dv[1];
598  }
599 
600  RKDUMB(vstart, 5, 0.0, etamax, m_xpoints, xx, ff);
601 
602  f2[0] = ff[1][m_xpoints] - 1;
603  f2[1] = ff[3][m_xpoints] - 1;
604 
605  al11 = (f1[0] - f[0]) / dv[0];
606  al21 = (f1[1] - f[1]) / dv[0];
607  al12 = (f2[0] - f[0]) / dv[1];
608  al22 = (f2[1] - f[1]) / dv[1];
609  det = al11 * al22 - al21 * al12;
610 
611  dv[0] = (-al22 * f[0] + al12 * f[1]) / det;
612  dv[1] = (al21 * f[0] - al11 * f[1]) / det;
613  v[0] = v[0] + dv[0];
614  v[1] = v[1] + dv[1];
615  }
616  }
617  }
618  }
619 
620  // Verification of the compressible similarity solution
621  ofstream verif;
622  verif.open("similarity_solution.dat");
623  for (i = 0; i < nQuadraturePts; i++)
624  {
625  verif << scientific << setprecision(9) << x_QuadraturePts[i] << " \t "
626  << y_QuadraturePts[i] << " \t ";
627  verif << scientific << setprecision(9) << u_QuadraturePts[i] << " \t "
628  << v_QuadraturePts[i] << " \t ";
629  verif << scientific << setprecision(9) << rho_QuadraturePts[i]
630  << " \t " << T_QuadraturePts[i] << endl;
631  }
632  verif.close();
633 
634  // Calculation of the physical variables
635  for (i = 0; i < nQuadraturePts; i++)
636  {
637  rho_QuadraturePts[i] = rho_QuadraturePts[i] * m_rhoInf;
638  u_QuadraturePts[i] = u_QuadraturePts[i] * m_uInf;
639  v_QuadraturePts[i] = v_QuadraturePts[i] * m_uInf;
640  T_QuadraturePts[i] = T_QuadraturePts[i] * m_Tinf;
641 
642  T_QuadraturePts[i] = T_QuadraturePts[i] * rho_QuadraturePts[i] * m_R;
643  T_QuadraturePts[i] = T_QuadraturePts[i] / (m_Gamma - 1);
644  T_QuadraturePts[i] =
645  T_QuadraturePts[i] +
646  0.5 * rho_QuadraturePts[i] *
647  (pow(u_QuadraturePts[i], 2.0) + pow(v_QuadraturePts[i], 2.0));
648 
649  u_QuadraturePts[i] = u_QuadraturePts[i] * rho_QuadraturePts[i];
650  v_QuadraturePts[i] = v_QuadraturePts[i] * rho_QuadraturePts[i];
651  }
652  string file_name;
653  if (expdim == 2)
654  {
655  MultiRegions::ContFieldSharedPtr Domain;
656  Domain = MemoryManager<MultiRegions::ContField>::AllocateSharedPtr(
657  vSession, graphShPt, vSession->GetVariable(0));
658 
659  Array<OneD, MultiRegions::ExpListSharedPtr> Exp(4);
660  MultiRegions::ExpListSharedPtr Exp2D_uk;
661  Exp2D_uk = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
662  vSession, graphShPt);
663 
664  MultiRegions::ExpListSharedPtr Exp2D_vk;
665  Exp2D_vk = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
666  vSession, graphShPt);
667 
668  MultiRegions::ExpListSharedPtr Exp2D_rhok;
669  Exp2D_rhok = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
670  vSession, graphShPt);
671 
672  MultiRegions::ExpListSharedPtr Exp2D_Tk;
673  Exp2D_Tk = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
674  vSession, graphShPt);
675 
676  // Filling the 2D expansion using a recursive algorithm based on the
677  // mesh ordering
678  LibUtilities::BasisSharedPtr Basis;
679  Basis = Domain->GetExp(0)->GetBasis(0);
680  numModes = Basis->GetNumModes();
681 
682  std::cout << "Number of modes = " << numModes << std::endl;
683 
684  // Copying the ukGlobal vector in m_phys (with the same pattern of
685  // m_phys)
686  Vmath::Vcopy(nQuadraturePts, u_QuadraturePts, 1, Exp2D_uk->UpdatePhys(),
687  1);
688  Vmath::Vcopy(nQuadraturePts, v_QuadraturePts, 1, Exp2D_vk->UpdatePhys(),
689  1);
690  Vmath::Vcopy(nQuadraturePts, rho_QuadraturePts, 1,
691  Exp2D_rhok->UpdatePhys(), 1);
692  Vmath::Vcopy(nQuadraturePts, T_QuadraturePts, 1, Exp2D_Tk->UpdatePhys(),
693  1);
694 
695  // Initialisation of the ExpList Exp
696  Exp[0] = Exp2D_rhok;
697  Exp[1] = Exp2D_uk;
698  Exp[2] = Exp2D_vk;
699  Exp[3] = Exp2D_Tk;
700 
701  // Expansion coefficient extraction (necessary to write the .fld file)
702  Exp[0]->FwdTrans(Exp2D_rhok->GetPhys(), Exp[0]->UpdateCoeffs());
703  Exp[1]->FwdTrans(Exp2D_uk->GetPhys(), Exp[1]->UpdateCoeffs());
704  Exp[2]->FwdTrans(Exp2D_vk->GetPhys(), Exp[2]->UpdateCoeffs());
705  Exp[3]->FwdTrans(Exp2D_Tk->GetPhys(), Exp[3]->UpdateCoeffs());
706 
707  // Definition of the name of the .fld file
708  cout << argv[1] << endl;
709  string tmp = argv[1];
710  int len = tmp.size();
711  for (i = 0; i < len - 4; ++i)
712  {
713  file_name += argv[1][i];
714  }
715  file_name = file_name + ".rst";
716 
717  // Definition of the Field
718  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
719  Exp[0]->GetFieldDefinitions();
720  std::vector<std::vector<NekDouble>> FieldData(FieldDef.size());
721 
722  for (j = 0; j < 4; j++)
723  {
724  for (i = 0; i < FieldDef.size(); i++)
725  {
726  if (j == 0)
727  {
728  FieldDef[i]->m_fields.push_back("rho");
729  }
730  else if (j == 1)
731  {
732  FieldDef[i]->m_fields.push_back("rhou");
733  }
734  else if (j == 2)
735  {
736  FieldDef[i]->m_fields.push_back("rhov");
737  }
738  else if (j == 3)
739  {
740  FieldDef[i]->m_fields.push_back("E");
741  }
742  Exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
743  }
744  }
745 
746  LibUtilities::Write(file_name, FieldDef, FieldData);
747  }
748  else if (expdim == 3)
749  {
750  MultiRegions::ContFieldSharedPtr Domain;
751  Domain = MemoryManager<MultiRegions::ContField>::AllocateSharedPtr(
752  vSession, graphShPt, vSession->GetVariable(0));
753 
754  Array<OneD, NekDouble> w_QuadraturePts;
755  w_QuadraturePts = Array<OneD, NekDouble>(nQuadraturePts, 0.0);
756  Array<OneD, MultiRegions::ExpListSharedPtr> Exp(5);
757 
758  MultiRegions::ExpListSharedPtr Exp3D_uk;
759  Exp3D_uk = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
760  vSession, graphShPt);
761 
762  MultiRegions::ExpListSharedPtr Exp3D_vk;
763  Exp3D_vk = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
764  vSession, graphShPt);
765 
766  MultiRegions::ExpListSharedPtr Exp3D_wk;
767  Exp3D_wk = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
768  vSession, graphShPt);
769 
770  MultiRegions::ExpListSharedPtr Exp3D_rhok;
771  Exp3D_rhok = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
772  vSession, graphShPt);
773 
774  MultiRegions::ExpListSharedPtr Exp3D_Tk;
775  Exp3D_Tk = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
776  vSession, graphShPt);
777 
778  // Filling the 3D expansion using a recursive algorithm based
779  // on the mesh ordering
780  LibUtilities::BasisSharedPtr Basis;
781  Basis = Domain->GetExp(0)->GetBasis(0);
782  numModes = Basis->GetNumModes();
783 
784  std::cout << "Number of modes = " << numModes << std::endl;
785 
786  // Copying the ukGlobal vector in m_phys (with the same pattern
787  // of m_phys)
788  Vmath::Vcopy(nQuadraturePts, rho_QuadraturePts, 1,
789  Exp3D_rhok->UpdatePhys(), 1);
790  Vmath::Vcopy(nQuadraturePts, u_QuadraturePts, 1, Exp3D_uk->UpdatePhys(),
791  1);
792  Vmath::Vcopy(nQuadraturePts, w_QuadraturePts, 1, Exp3D_vk->UpdatePhys(),
793  1);
794  Vmath::Vcopy(nQuadraturePts, v_QuadraturePts, 1, Exp3D_wk->UpdatePhys(),
795  1);
796  Vmath::Vcopy(nQuadraturePts, T_QuadraturePts, 1, Exp3D_Tk->UpdatePhys(),
797  1);
798 
799  // Initialisation of the ExpList Exp
800  Exp[0] = Exp3D_rhok;
801  Exp[1] = Exp3D_uk;
802  Exp[2] = Exp3D_vk;
803  Exp[3] = Exp3D_wk;
804  Exp[4] = Exp3D_Tk;
805 
806  // Expansion coefficient extraction (necessary to write the .fld file)
807  Exp[0]->FwdTrans(Exp3D_rhok->GetPhys(), Exp[0]->UpdateCoeffs());
808  Exp[1]->FwdTrans(Exp3D_uk->GetPhys(), Exp[1]->UpdateCoeffs());
809  Exp[2]->FwdTrans(Exp3D_vk->GetPhys(), Exp[2]->UpdateCoeffs());
810  Exp[3]->FwdTrans(Exp3D_wk->GetPhys(), Exp[3]->UpdateCoeffs());
811  Exp[4]->FwdTrans(Exp3D_Tk->GetPhys(), Exp[4]->UpdateCoeffs());
812 
813  // Definition of the name of the .fld file
814  cout << argv[1] << endl;
815  string tmp = argv[1];
816  int len = tmp.size();
817  for (i = 0; i < len - 4; ++i)
818  {
819  file_name += argv[1][i];
820  }
821  file_name = file_name + ".rst";
822 
823  // Definition of the Field
824  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
825  Exp[0]->GetFieldDefinitions();
826  std::vector<std::vector<NekDouble>> FieldData(FieldDef.size());
827 
828  for (j = 0; j < 5; j++)
829  {
830  for (i = 0; i < FieldDef.size(); i++)
831  {
832  if (j == 0)
833  {
834  FieldDef[i]->m_fields.push_back("rho");
835  }
836  else if (j == 1)
837  {
838  FieldDef[i]->m_fields.push_back("rhou");
839  }
840  else if (j == 2)
841  {
842  FieldDef[i]->m_fields.push_back("rhov");
843  }
844  else if (j == 3)
845  {
846  FieldDef[i]->m_fields.push_back("rhow");
847  }
848  else if (j == 4)
849  {
850  FieldDef[i]->m_fields.push_back("E");
851  }
852  Exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
853  }
854  }
855 
856  LibUtilities::Write(file_name, FieldDef, FieldData);
857  }
858 
859  std::cout << "----------------------------------------------------\n";
860  std::cout << "\n=================================================\n";
861  std::cout << "Similarity solution \n";
862  std::cout << "===================================================\n";
863  std::cout << "***************************************************\n";
864  std::cout << "DATA FROM THE SESSION FILE:\n";
865  std::cout << "Reynolds number = " << m_Re << "\t[-]"
866  << std::endl;
867  std::cout << "Mach number = " << m_Mach << "\t[-]"
868  << std::endl;
869  std::cout << "Characteristic length = " << m_long << "\t[m]"
870  << std::endl;
871  std::cout << "U_infinity = " << m_uInf << "\t[m/s]"
872  << std::endl;
873  std::cout << "***************************************************\n";
874  std::cout << "---------------------------------------------------\n";
875  std::cout << "MESH and EXPANSION DATA:\n";
876  std::cout << "Done." << std::endl;
877 
878  return 0;
879 }
m_Suth
NekDouble m_Suth
Definition: CompressibleBL.cpp:76
m_mu
NekDouble m_mu
Definition: CompressibleBL.cpp:86
main
int main(int argc, char *argv[])
Definition: CompressibleBL.cpp:337
m_Gamma
NekDouble m_Gamma
Definition: CompressibleBL.cpp:79
m_Tinf
NekDouble m_Tinf
Definition: CompressibleBL.cpp:75
m_R
NekDouble m_R
Definition: CompressibleBL.cpp:84
m_Twall
NekDouble m_Twall
Definition: CompressibleBL.cpp:78
etamax
const NekDouble etamax
Definition: CompressibleBL.cpp:93
Nvisc
const NekDouble Nvisc
Definition: CompressibleBL.cpp:91
m_Tw
NekDouble m_Tw
Definition: CompressibleBL.cpp:77
errtol
const NekDouble errtol
Definition: CompressibleBL.cpp:94
RKDUMB
void RKDUMB(Array< OneD, NekDouble > vstart, int nvar, NekDouble x1, NekDouble x2, int m_xpoints, Array< OneD, NekDouble > xx, Array< OneD, Array< OneD, NekDouble >> y)
Definition: CompressibleBL.cpp:177
m_Mach
NekDouble m_Mach
Definition: CompressibleBL.cpp:73
m_rhoInf
NekDouble m_rhoInf
Definition: CompressibleBL.cpp:83
m_long
NekDouble m_long
Definition: CompressibleBL.cpp:81
COMPBL
void COMPBL(Array< OneD, NekDouble > v, Array< OneD, NekDouble > dv)
Definition: CompressibleBL.cpp:99
m_xpoints
const int m_xpoints
Definition: CompressibleBL.cpp:89
m_vInf
NekDouble m_vInf
Definition: CompressibleBL.cpp:85
OUTPUT
void OUTPUT(int m_xpoints, Array< OneD, NekDouble > xx, Array< OneD, Array< OneD, NekDouble >> ff, int nQuadraturePts, Array< OneD, NekDouble > x_QuadraturePts, Array< OneD, NekDouble > y_QuadraturePts, Array< OneD, NekDouble > u_QuadraturePts, Array< OneD, NekDouble > v_QuadraturePts, Array< OneD, NekDouble > rho_QuadraturePts, Array< OneD, NekDouble > T_QuadraturePts)
Definition: CompressibleBL.cpp:219
m_uInf
NekDouble m_uInf
Definition: CompressibleBL.cpp:82
m_Re
NekDouble m_Re
Definition: CompressibleBL.cpp:72
m_To
NekDouble m_To
Definition: CompressibleBL.cpp:87
Omega
const NekDouble Omega
Definition: CompressibleBL.cpp:92
L
NekDouble L
Definition: CompressibleBL.cpp:74
m_Pr
NekDouble m_Pr
Definition: CompressibleBL.cpp:80
RK4
void RK4(Array< OneD, NekDouble > y, Array< OneD, NekDouble > dydx, int n, NekDouble x, NekDouble h, Array< OneD, NekDouble > yout)
Definition: CompressibleBL.cpp:133
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
Nektar::LibUtilities::Basis
Represents a basis of a given type.
Definition: Basis.h:213
Nektar::LibUtilities::Basis::GetNumModes
int GetNumModes() const
Return order of basis from the basis specification.
Definition: Basis.h:222
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:83
Nektar::LibUtilities::Write
void Write(const std::string &outFile, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble >> &fielddata, const FieldMetaDataMap &fieldinfomap, const bool backup)
This function allows for data to be written to an FLD file when a session and/or communicator is not ...
Definition: FieldIO.cpp:247
Nektar::LibUtilities::BasisSharedPtr
std::shared_ptr< Basis > BasisSharedPtr
Definition: FoundationsFwd.hpp:71
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:115
Nektar::MultiRegions::ExpListSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:55
Nektar::MultiRegions::ContFieldSharedPtr
std::shared_ptr< ContField > ContFieldSharedPtr
Definition: ContField.h:289
Nektar::SpatialDomains::MeshGraphSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:172
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:2
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294
Generated on Sat Mar 25 2023 12:35:49 for Nektar++ by doxygen 1.9.1