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