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