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