Nektar++
Functions | Variables
CompressibleBL.cpp File Reference
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <string>
#include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
#include <MultiRegions/ContField.h>
#include <MultiRegions/DisContField.h>
#include <MultiRegions/ExpList2DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous2D.h>
#include <LocalRegions/Expansion2D.h>
#include <LocalRegions/Expansion3D.h>
#include <LocalRegions/MatrixKey.h>
#include <LibUtilities/BasicUtils/FieldIO.h>
#include <LibUtilities/BasicUtils/NekFactory.hpp>
#include <LibUtilities/BasicUtils/SessionReader.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/Memory/NekMemoryManager.hpp>
#include <SpatialDomains/MeshGraphIO.h>
#include <SolverUtils/SolverUtilsDeclspec.h>

Go to the source code of this file.

Functions

void COMPBL (Array< OneD, NekDouble > v, Array< OneD, NekDouble > dv)
 
void RK4 (Array< OneD, NekDouble > y, Array< OneD, NekDouble > dydx, int n, NekDouble x, NekDouble h, Array< OneD, NekDouble > yout)
 
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)
 
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)
 
int main (int argc, char *argv[])
 

Variables

NekDouble m_Re
 
NekDouble m_Mach
 
NekDouble L
 
NekDouble m_Tinf
 
NekDouble m_Suth
 
NekDouble m_Tw
 
NekDouble m_Twall
 
NekDouble m_Gamma
 
NekDouble m_Pr
 
NekDouble m_long
 
NekDouble m_uInf
 
NekDouble m_rhoInf
 
NekDouble m_R
 
NekDouble m_vInf
 
NekDouble m_mu
 
NekDouble m_To = 273.11
 
const int m_xpoints = 1000001
 
const NekDouble Nvisc = 1
 
const NekDouble Omega = 1
 
const NekDouble etamax = 10.0
 
const NekDouble errtol = 1e-5
 

Function Documentation

◆ COMPBL()

void COMPBL ( Array< OneD, NekDouble >  v,
Array< OneD, NekDouble >  dv 
)

Calculate the compressible boundary layer using the similarity solution

Definition at line 94 of file CompressibleBL.cpp.

95{
96 NekDouble c, dcdg, cp;
97
98 if (Nvisc == 1)
99 {
100 c = sqrt(v[3]) * (1.0 + m_Suth) / (v[3] + m_Suth);
101 dcdg = 1.0 / (2. * sqrt(v[3])) - sqrt(v[3]) / (v[3] + m_Suth);
102 dcdg = dcdg * (1.0 + m_Suth) / (v[3] + m_Suth);
103 cp = dcdg * v[4];
104 }
105 if (Nvisc == 2)
106 {
107 c = pow(v[3], (Omega - 1.0));
108 dcdg = (Omega - 1.0) * pow(v[3], (Omega - 2.0));
109 cp = dcdg * v[4];
110 }
111 if (Nvisc == 3)
112 {
113 c = sqrt(m_Twall) * (1.0 + m_Suth) / (m_Suth + m_Twall);
114 cp = 0.0;
115 }
116
117 dv[0] = v[1];
118 dv[1] = v[2];
119 dv[2] = -v[2] * (cp + v[0]) / c;
120 dv[3] = v[4];
121 dv[4] = -v[4] * (cp + m_Pr * v[0]) / c -
122 m_Pr * (m_Gamma - 1.0) * pow(m_Mach, 2.0) * pow(v[2], 2);
123}
m_Suth
NekDouble m_Suth
Definition: CompressibleBL.cpp:71
m_Gamma
NekDouble m_Gamma
Definition: CompressibleBL.cpp:74
m_Twall
NekDouble m_Twall
Definition: CompressibleBL.cpp:73
Nvisc
const NekDouble Nvisc
Definition: CompressibleBL.cpp:86
m_Mach
NekDouble m_Mach
Definition: CompressibleBL.cpp:68
Omega
const NekDouble Omega
Definition: CompressibleBL.cpp:87
m_Pr
NekDouble m_Pr
Definition: CompressibleBL.cpp:75
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:294

References m_Gamma, m_Mach, m_Pr, m_Suth, m_Twall, Nvisc, Omega, and tinysimd::sqrt().

Referenced by OUTPUT(), RK4(), and RKDUMB().

◆ main()

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

Calculate the compressible boundary layer solution for a given 3D mesh and dump the solution into a .rst file.

Definition at line 330 of file CompressibleBL.cpp.

331{
332 // Variable initialisation
333 int nmax = 5;
334 int maxit = 10;
335
336 string opt;
337
338 int i, j, k, numModes;
339
340 Array<OneD, NekDouble> xx(m_xpoints, 0.0);
341 Array<OneD, Array<OneD, NekDouble>> ff(nmax);
342 Array<OneD, NekDouble> parameter(9, 0.0);
343
344 for (i = 0; i < nmax; i++)
345 {
346 ff[i] = Array<OneD, NekDouble>(m_xpoints);
347 }
348
349 Array<OneD, NekDouble> vstart(nmax, 0.0);
350 Array<OneD, NekDouble> v(2);
351 Array<OneD, NekDouble> dv(2);
352 Array<OneD, NekDouble> f(2);
353 Array<OneD, NekDouble> f1(2);
354 Array<OneD, NekDouble> f2(2);
355
356 NekDouble al11, al21, al12, al22, det;
357
358 // Reading the session file
359 LibUtilities::SessionReaderSharedPtr vSession =
360 LibUtilities::SessionReader::CreateInstance(argc, argv);
361
362 // Read in mesh from input file and create an object of class MeshGraph
363 SpatialDomains::MeshGraphSharedPtr graphShPt =
364 SpatialDomains::MeshGraphIO::Read(vSession);
365
366 int expdim = graphShPt->GetMeshDimension();
367
368 int nElements, nQuadraturePts = 0;
369 Array<OneD, NekDouble> x_QuadraturePts;
370 Array<OneD, NekDouble> y_QuadraturePts;
371 Array<OneD, NekDouble> z_QuadraturePts;
372
373 MultiRegions::ContFieldSharedPtr Domain;
374 Domain = MemoryManager<MultiRegions::ContField>::AllocateSharedPtr(
375 vSession, graphShPt, vSession->GetVariable(0));
376
377 // Get the total number of elements
378 nElements = Domain->GetExpSize();
379 std::cout << "Number of elements = " << nElements
380 << std::endl;
381
382 // Get the total number of quadrature points (depends on n. modes)
383 nQuadraturePts = Domain->GetTotPoints();
384 std::cout << "Number of quadrature points = " << nQuadraturePts
385 << std::endl;
386
387 // Coordinates of the quadrature points
388 x_QuadraturePts = Array<OneD, NekDouble>(nQuadraturePts);
389 y_QuadraturePts = Array<OneD, NekDouble>(nQuadraturePts);
390 z_QuadraturePts = Array<OneD, NekDouble>(nQuadraturePts);
391 Domain->GetCoords(x_QuadraturePts, y_QuadraturePts, z_QuadraturePts);
392
393 if (expdim == 1)
394 {
395 ASSERTL0(false, "Routine available for 2D and 3D problem only.")
396 }
397
398 // Loading parameters from session file
399 vSession->LoadParameter("Re", m_Re, 1.0);
400 vSession->LoadParameter("Mach", m_Mach, 1.0);
401 vSession->LoadParameter("TInf", m_Tinf, 1.0);
402 vSession->LoadParameter("Twall", m_Twall, 1.0);
403 vSession->LoadParameter("Gamma", m_Gamma, 1.0);
404 vSession->LoadParameter("Pr", m_Pr, 1.0);
405 vSession->LoadParameter("L", m_long, 1.0);
406 vSession->LoadParameter("rhoInf", m_rhoInf, 1.0);
407 vSession->LoadParameter("uInf", m_uInf, 1.0);
408 vSession->LoadParameter("GasConstant", m_R, 1.0);
409 vSession->LoadParameter("vInf", m_vInf, 1.0);
410 vSession->LoadParameter("mu", m_mu, 1.0);
411
412 // Rescaling factors
413 m_Suth = 110.4 / m_Tinf;
414 m_Tw = m_Twall / m_Tinf;
415 m_Re = m_Re / m_long;
416
417 cout << "Number of points"
418 << " " << m_xpoints << endl;
419
420 // Defining the solution arrays
421 Array<OneD, NekDouble> u_QuadraturePts(nQuadraturePts, 0.0);
422 Array<OneD, NekDouble> v_QuadraturePts(nQuadraturePts, 0.0);
423 Array<OneD, NekDouble> rho_QuadraturePts(nQuadraturePts, 0.0);
424 Array<OneD, NekDouble> T_QuadraturePts(nQuadraturePts, 0.0);
425
426 // Calculation of the similarity variables
427 if (m_Tw > 0)
428 {
429 vstart[3] = m_Tw;
430 }
431 if (m_Tw < 0.0)
432 {
433 v[1] = 1.0 + 0.5 * 0.84 * (m_Gamma - 1) * (m_Mach * m_Mach);
434 v[0] = 0.47 * pow(v[1], 0.21);
435 }
436 else
437 {
438 v[1] = 0.062 * pow(m_Mach, 2) -
439 0.1 * (m_Tw - 1.0) * (10 + m_Mach) / (0.2 + m_Mach);
440 v[0] = 0.45 - 0.01 * m_Mach + (m_Tw - 1.0) * 0.06;
441 m_Twall = m_Tw;
442 }
443
444 dv[0] = v[0] * 0.01;
445
446 if (m_Tw < 0.0)
447 {
448 dv[1] = v[1] * 0.01;
449 }
450 else
451 {
452 dv[1] = 0.1;
453 }
454
455 vstart[2] = v[0];
456
457 if (m_Tw < 0)
458 {
459 vstart[3] = v[1];
460 m_Twall = vstart[3];
461 }
462 else
463 {
464 vstart[4] = v[1];
465 }
466
467 RKDUMB(vstart, 5, 0.0, etamax, m_xpoints, xx, ff);
468
469 for (k = 0; k < maxit; k++)
470 {
471 vstart[2] = v[0];
472
473 if (m_Tw < 0)
474 {
475 vstart[3] = v[1];
476 m_Twall = vstart[3];
477 }
478 else
479 {
480 vstart[4] = v[1];
481 }
482
483 RKDUMB(vstart, 5, 0.0, etamax, m_xpoints, xx, ff);
484
485 NekDouble err = fabs(ff[1][m_xpoints] - 1) + fabs(ff[3][m_xpoints] - 1);
486
487 cout << "err" << scientific << setprecision(9) << " " << err << endl;
488
489 if (expdim == 2)
490 {
491 if (err < errtol)
492 {
493 cout << "Calculating" << endl;
494 OUTPUT(m_xpoints, xx, ff, nQuadraturePts, x_QuadraturePts,
495 y_QuadraturePts, u_QuadraturePts, v_QuadraturePts,
496 rho_QuadraturePts, T_QuadraturePts);
497 break;
498 }
499 else
500 {
501 f[0] = ff[1][m_xpoints] - 1;
502 f[1] = ff[3][m_xpoints] - 1;
503 vstart[2] = v[0] + dv[0];
504
505 if (m_Tw < 0)
506 {
507 vstart[3] = v[1];
508 m_Twall = vstart[3];
509 }
510 else
511 {
512 vstart[4] = v[1];
513 }
514
515 RKDUMB(vstart, 5, 0.0, etamax, m_xpoints, xx, ff);
516
517 f1[0] = ff[1][m_xpoints] - 1;
518 f1[1] = ff[3][m_xpoints] - 1;
519
520 vstart[2] = v[0];
521
522 if (m_Tw < 0)
523 {
524 vstart[3] = v[1] + dv[1];
525 m_Twall = vstart[3];
526 }
527 else
528 {
529 vstart[4] = v[1] + dv[1];
530 }
531
532 RKDUMB(vstart, 5, 0.0, etamax, m_xpoints, xx, ff);
533
534 f2[0] = ff[1][m_xpoints] - 1;
535 f2[1] = ff[3][m_xpoints] - 1;
536
537 al11 = (f1[0] - f[0]) / dv[0];
538 al21 = (f1[1] - f[1]) / dv[0];
539 al12 = (f2[0] - f[0]) / dv[1];
540 al22 = (f2[1] - f[1]) / dv[1];
541 det = al11 * al22 - al21 * al12;
542
543 dv[0] = (-al22 * f[0] + al12 * f[1]) / det;
544 dv[1] = (al21 * f[0] - al11 * f[1]) / det;
545 v[0] = v[0] + dv[0];
546 v[1] = v[1] + dv[1];
547 }
548 }
549 else if (expdim == 3)
550 {
551 {
552 if (err < errtol)
553 {
554 cout << "Calculating" << endl;
555 OUTPUT(m_xpoints, xx, ff, nQuadraturePts, x_QuadraturePts,
556 z_QuadraturePts, u_QuadraturePts, v_QuadraturePts,
557 rho_QuadraturePts, T_QuadraturePts);
558 break;
559 }
560 else
561 {
562 f[0] = ff[1][m_xpoints] - 1;
563 f[1] = ff[3][m_xpoints] - 1;
564 vstart[2] = v[0] + dv[0];
565
566 if (m_Tw < 0)
567 {
568 vstart[3] = v[1];
569 m_Twall = vstart[3];
570 }
571 else
572 {
573 vstart[4] = v[1];
574 }
575
576 RKDUMB(vstart, 5, 0.0, etamax, m_xpoints, xx, ff);
577
578 f1[0] = ff[1][m_xpoints] - 1;
579 f1[1] = ff[3][m_xpoints] - 1;
580
581 vstart[2] = v[0];
582
583 if (m_Tw < 0)
584 {
585 vstart[3] = v[1] + dv[1];
586 m_Twall = vstart[3];
587 }
588 else
589 {
590 vstart[4] = v[1] + dv[1];
591 }
592
593 RKDUMB(vstart, 5, 0.0, etamax, m_xpoints, xx, ff);
594
595 f2[0] = ff[1][m_xpoints] - 1;
596 f2[1] = ff[3][m_xpoints] - 1;
597
598 al11 = (f1[0] - f[0]) / dv[0];
599 al21 = (f1[1] - f[1]) / dv[0];
600 al12 = (f2[0] - f[0]) / dv[1];
601 al22 = (f2[1] - f[1]) / dv[1];
602 det = al11 * al22 - al21 * al12;
603
604 dv[0] = (-al22 * f[0] + al12 * f[1]) / det;
605 dv[1] = (al21 * f[0] - al11 * f[1]) / det;
606 v[0] = v[0] + dv[0];
607 v[1] = v[1] + dv[1];
608 }
609 }
610 }
611 }
612
613 // Verification of the compressible similarity solution
614 ofstream verif;
615 verif.open("similarity_solution.dat");
616 for (i = 0; i < nQuadraturePts; i++)
617 {
618 verif << scientific << setprecision(9) << x_QuadraturePts[i] << " \t "
619 << y_QuadraturePts[i] << " \t ";
620 verif << scientific << setprecision(9) << u_QuadraturePts[i] << " \t "
621 << v_QuadraturePts[i] << " \t ";
622 verif << scientific << setprecision(9) << rho_QuadraturePts[i]
623 << " \t " << T_QuadraturePts[i] << endl;
624 }
625 verif.close();
626
627 // Calculation of the physical variables
628 for (i = 0; i < nQuadraturePts; i++)
629 {
630 rho_QuadraturePts[i] = rho_QuadraturePts[i] * m_rhoInf;
631 u_QuadraturePts[i] = u_QuadraturePts[i] * m_uInf;
632 v_QuadraturePts[i] = v_QuadraturePts[i] * m_uInf;
633 T_QuadraturePts[i] = T_QuadraturePts[i] * m_Tinf;
634
635 T_QuadraturePts[i] = T_QuadraturePts[i] * rho_QuadraturePts[i] * m_R;
636 T_QuadraturePts[i] = T_QuadraturePts[i] / (m_Gamma - 1);
637 T_QuadraturePts[i] =
638 T_QuadraturePts[i] +
639 0.5 * rho_QuadraturePts[i] *
640 (pow(u_QuadraturePts[i], 2.0) + pow(v_QuadraturePts[i], 2.0));
641
642 u_QuadraturePts[i] = u_QuadraturePts[i] * rho_QuadraturePts[i];
643 v_QuadraturePts[i] = v_QuadraturePts[i] * rho_QuadraturePts[i];
644 }
645 string file_name;
646 if (expdim == 2)
647 {
648 MultiRegions::ContFieldSharedPtr Domain;
649 Domain = MemoryManager<MultiRegions::ContField>::AllocateSharedPtr(
650 vSession, graphShPt, vSession->GetVariable(0));
651
652 Array<OneD, MultiRegions::ExpListSharedPtr> Exp(4);
653 MultiRegions::ExpListSharedPtr Exp2D_uk;
654 Exp2D_uk = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
655 vSession, graphShPt);
656
657 MultiRegions::ExpListSharedPtr Exp2D_vk;
658 Exp2D_vk = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
659 vSession, graphShPt);
660
661 MultiRegions::ExpListSharedPtr Exp2D_rhok;
662 Exp2D_rhok = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
663 vSession, graphShPt);
664
665 MultiRegions::ExpListSharedPtr Exp2D_Tk;
666 Exp2D_Tk = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
667 vSession, graphShPt);
668
669 // Filling the 2D expansion using a recursive algorithm based on the
670 // mesh ordering
671 LibUtilities::BasisSharedPtr Basis;
672 Basis = Domain->GetExp(0)->GetBasis(0);
673 numModes = Basis->GetNumModes();
674
675 std::cout << "Number of modes = " << numModes << std::endl;
676
677 // Copying the ukGlobal vector in m_phys (with the same pattern of
678 // m_phys)
679 Vmath::Vcopy(nQuadraturePts, u_QuadraturePts, 1, Exp2D_uk->UpdatePhys(),
680 1);
681 Vmath::Vcopy(nQuadraturePts, v_QuadraturePts, 1, Exp2D_vk->UpdatePhys(),
682 1);
683 Vmath::Vcopy(nQuadraturePts, rho_QuadraturePts, 1,
684 Exp2D_rhok->UpdatePhys(), 1);
685 Vmath::Vcopy(nQuadraturePts, T_QuadraturePts, 1, Exp2D_Tk->UpdatePhys(),
686 1);
687
688 // Initialisation of the ExpList Exp
689 Exp[0] = Exp2D_rhok;
690 Exp[1] = Exp2D_uk;
691 Exp[2] = Exp2D_vk;
692 Exp[3] = Exp2D_Tk;
693
694 // Expansion coefficient extraction (necessary to write the .fld file)
695 Exp[0]->FwdTrans(Exp2D_rhok->GetPhys(), Exp[0]->UpdateCoeffs());
696 Exp[1]->FwdTrans(Exp2D_uk->GetPhys(), Exp[1]->UpdateCoeffs());
697 Exp[2]->FwdTrans(Exp2D_vk->GetPhys(), Exp[2]->UpdateCoeffs());
698 Exp[3]->FwdTrans(Exp2D_Tk->GetPhys(), Exp[3]->UpdateCoeffs());
699
700 // Definition of the name of the .fld file
701 cout << argv[1] << endl;
702 string tmp = argv[1];
703 int len = tmp.size();
704 for (i = 0; i < len - 4; ++i)
705 {
706 file_name += argv[1][i];
707 }
708 file_name = file_name + ".rst";
709
710 // Definition of the Field
711 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
712 Exp[0]->GetFieldDefinitions();
713 std::vector<std::vector<NekDouble>> FieldData(FieldDef.size());
714
715 for (j = 0; j < 4; j++)
716 {
717 for (i = 0; i < FieldDef.size(); i++)
718 {
719 if (j == 0)
720 {
721 FieldDef[i]->m_fields.push_back("rho");
722 }
723 else if (j == 1)
724 {
725 FieldDef[i]->m_fields.push_back("rhou");
726 }
727 else if (j == 2)
728 {
729 FieldDef[i]->m_fields.push_back("rhov");
730 }
731 else if (j == 3)
732 {
733 FieldDef[i]->m_fields.push_back("E");
734 }
735 Exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
736 }
737 }
738
739 LibUtilities::Write(file_name, FieldDef, FieldData);
740 }
741 else if (expdim == 3)
742 {
743 MultiRegions::ContFieldSharedPtr Domain;
744 Domain = MemoryManager<MultiRegions::ContField>::AllocateSharedPtr(
745 vSession, graphShPt, vSession->GetVariable(0));
746
747 Array<OneD, NekDouble> w_QuadraturePts;
748 w_QuadraturePts = Array<OneD, NekDouble>(nQuadraturePts, 0.0);
749 Array<OneD, MultiRegions::ExpListSharedPtr> Exp(5);
750
751 MultiRegions::ExpListSharedPtr Exp3D_uk;
752 Exp3D_uk = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
753 vSession, graphShPt);
754
755 MultiRegions::ExpListSharedPtr Exp3D_vk;
756 Exp3D_vk = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
757 vSession, graphShPt);
758
759 MultiRegions::ExpListSharedPtr Exp3D_wk;
760 Exp3D_wk = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
761 vSession, graphShPt);
762
763 MultiRegions::ExpListSharedPtr Exp3D_rhok;
764 Exp3D_rhok = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
765 vSession, graphShPt);
766
767 MultiRegions::ExpListSharedPtr Exp3D_Tk;
768 Exp3D_Tk = MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(
769 vSession, graphShPt);
770
771 // Filling the 3D expansion using a recursive algorithm based
772 // on the mesh ordering
773 LibUtilities::BasisSharedPtr Basis;
774 Basis = Domain->GetExp(0)->GetBasis(0);
775 numModes = Basis->GetNumModes();
776
777 std::cout << "Number of modes = " << numModes << std::endl;
778
779 // Copying the ukGlobal vector in m_phys (with the same pattern
780 // of m_phys)
781 Vmath::Vcopy(nQuadraturePts, rho_QuadraturePts, 1,
782 Exp3D_rhok->UpdatePhys(), 1);
783 Vmath::Vcopy(nQuadraturePts, u_QuadraturePts, 1, Exp3D_uk->UpdatePhys(),
784 1);
785 Vmath::Vcopy(nQuadraturePts, w_QuadraturePts, 1, Exp3D_vk->UpdatePhys(),
786 1);
787 Vmath::Vcopy(nQuadraturePts, v_QuadraturePts, 1, Exp3D_wk->UpdatePhys(),
788 1);
789 Vmath::Vcopy(nQuadraturePts, T_QuadraturePts, 1, Exp3D_Tk->UpdatePhys(),
790 1);
791
792 // Initialisation of the ExpList Exp
793 Exp[0] = Exp3D_rhok;
794 Exp[1] = Exp3D_uk;
795 Exp[2] = Exp3D_vk;
796 Exp[3] = Exp3D_wk;
797 Exp[4] = Exp3D_Tk;
798
799 // Expansion coefficient extraction (necessary to write the .fld file)
800 Exp[0]->FwdTrans(Exp3D_rhok->GetPhys(), Exp[0]->UpdateCoeffs());
801 Exp[1]->FwdTrans(Exp3D_uk->GetPhys(), Exp[1]->UpdateCoeffs());
802 Exp[2]->FwdTrans(Exp3D_vk->GetPhys(), Exp[2]->UpdateCoeffs());
803 Exp[3]->FwdTrans(Exp3D_wk->GetPhys(), Exp[3]->UpdateCoeffs());
804 Exp[4]->FwdTrans(Exp3D_Tk->GetPhys(), Exp[4]->UpdateCoeffs());
805
806 // Definition of the name of the .fld file
807 cout << argv[1] << endl;
808 string tmp = argv[1];
809 int len = tmp.size();
810 for (i = 0; i < len - 4; ++i)
811 {
812 file_name += argv[1][i];
813 }
814 file_name = file_name + ".rst";
815
816 // Definition of the Field
817 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
818 Exp[0]->GetFieldDefinitions();
819 std::vector<std::vector<NekDouble>> FieldData(FieldDef.size());
820
821 for (j = 0; j < 5; j++)
822 {
823 for (i = 0; i < FieldDef.size(); i++)
824 {
825 if (j == 0)
826 {
827 FieldDef[i]->m_fields.push_back("rho");
828 }
829 else if (j == 1)
830 {
831 FieldDef[i]->m_fields.push_back("rhou");
832 }
833 else if (j == 2)
834 {
835 FieldDef[i]->m_fields.push_back("rhov");
836 }
837 else if (j == 3)
838 {
839 FieldDef[i]->m_fields.push_back("rhow");
840 }
841 else if (j == 4)
842 {
843 FieldDef[i]->m_fields.push_back("E");
844 }
845 Exp[j]->AppendFieldData(FieldDef[i], FieldData[i]);
846 }
847 }
848
849 LibUtilities::Write(file_name, FieldDef, FieldData);
850 }
851
852 std::cout << "----------------------------------------------------\n";
853 std::cout << "\n=================================================\n";
854 std::cout << "Similarity solution \n";
855 std::cout << "===================================================\n";
856 std::cout << "***************************************************\n";
857 std::cout << "DATA FROM THE SESSION FILE:\n";
858 std::cout << "Reynolds number = " << m_Re << "\t[-]"
859 << std::endl;
860 std::cout << "Mach number = " << m_Mach << "\t[-]"
861 << std::endl;
862 std::cout << "Characteristic length = " << m_long << "\t[m]"
863 << std::endl;
864 std::cout << "U_infinity = " << m_uInf << "\t[m/s]"
865 << std::endl;
866 std::cout << "***************************************************\n";
867 std::cout << "---------------------------------------------------\n";
868 std::cout << "MESH and EXPANSION DATA:\n";
869 std::cout << "Done." << std::endl;
870
871 return 0;
872}
m_mu
NekDouble m_mu
Definition: CompressibleBL.cpp:81
m_Tinf
NekDouble m_Tinf
Definition: CompressibleBL.cpp:70
m_R
NekDouble m_R
Definition: CompressibleBL.cpp:79
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:212
etamax
const NekDouble etamax
Definition: CompressibleBL.cpp:88
m_Tw
NekDouble m_Tw
Definition: CompressibleBL.cpp:72
errtol
const NekDouble errtol
Definition: CompressibleBL.cpp:89
m_rhoInf
NekDouble m_rhoInf
Definition: CompressibleBL.cpp:78
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:170
m_long
NekDouble m_long
Definition: CompressibleBL.cpp:76
m_xpoints
const int m_xpoints
Definition: CompressibleBL.cpp:84
m_vInf
NekDouble m_vInf
Definition: CompressibleBL.cpp:80
m_uInf
NekDouble m_uInf
Definition: CompressibleBL.cpp:77
m_Re
NekDouble m_Re
Definition: CompressibleBL.cpp:67
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
Nektar::LibUtilities::Basis
Represents a basis of a given type.
Definition: Basis.h:198
Nektar::LibUtilities::Basis::GetNumModes
int GetNumModes() const
Return order of basis from the basis specification.
Definition: Basis.h:207
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:81
Nektar::LibUtilities::BasisSharedPtr
std::shared_ptr< Basis > BasisSharedPtr
Definition: FoundationsFwd.hpp:69
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:113
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:245
Nektar::MultiRegions::ExpListSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:56
Nektar::MultiRegions::ContFieldSharedPtr
std::shared_ptr< ContField > ContFieldSharedPtr
Definition: ContField.h:268
Nektar::SpatialDomains::MeshGraphSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825

References ASSERTL0, errtol, etamax, Nektar::LibUtilities::Basis::GetNumModes(), m_Gamma, m_long, m_Mach, m_mu, m_Pr, m_R, m_Re, m_rhoInf, m_Suth, m_Tinf, m_Tw, m_Twall, m_uInf, m_vInf, m_xpoints, OUTPUT(), RKDUMB(), Vmath::Vcopy(), and Nektar::LibUtilities::Write().

◆ 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 
)

Create the output file

Definition at line 212 of file CompressibleBL.cpp.

220{
221 Array<OneD, NekDouble> z(m_xpoints, 0.0);
222 Array<OneD, NekDouble> v(m_xpoints, 0.0);
223 Array<OneD, NekDouble> dv(m_xpoints, 0.0);
224 Array<OneD, NekDouble> u(m_xpoints, 0.0);
225 Array<OneD, NekDouble> t(m_xpoints, 0.0);
226 Array<OneD, NekDouble> rho(m_xpoints, 0.0);
227 Array<OneD, NekDouble> mu(m_xpoints, 0.0);
228 Array<OneD, NekDouble> vv(m_xpoints, 0.0);
229 Array<OneD, NekDouble> velocity(m_xpoints, 0.0);
230 Array<OneD, NekDouble> test(m_xpoints, 0.0);
231
232 NekDouble dd, dm, scale;
233 NekDouble xcher, ycher;
234 int index = -1;
235
236 z[0] = 0.0;
237 NekDouble sumd = 0.0;
238
239 for (int i = 1; i < m_xpoints; i++)
240 {
241 z[i] = z[i - 1] + 0.5 * (xx[i] - xx[i - 1]) * (ff[3][i] + ff[3][i - 1]);
242 dm = ff[3][i - 1] - ff[1][i - 1];
243 dd = ff[3][i] - ff[1][i];
244 sumd = sumd + 0.5 * (xx[i] - xx[i - 1]) * (dd + dm);
245 }
246
247 scale = sumd;
248
249 ofstream file3;
250 file3.open("physical_data.dat");
251
252 NekDouble xin, rex, delsx, delta;
253
254 for (int i = 0; i < m_xpoints; i++)
255 {
256 for (int k = 0; k < 5; k++)
257 {
258 v[k] = ff[k][i];
259 }
260 COMPBL(v, dv);
261 u[i] = ff[1][i];
262 t[i] = ff[3][i];
263 rho[i] = (1.0 / ff[3][i]);
264 vv[i] = -ff[0][i] / sqrt(m_uInf);
265 mu[i] = pow(t[i], 1.5) * (1 + m_Suth) / (t[i] + m_Suth) / (m_Re);
266 velocity[i] = ff[0][i];
267 }
268
269 NekDouble scale2, coeff;
270
271 for (int i = 0; i < nQuadraturePts; i++)
272 {
273 if (i % 100000 == 0)
274 {
275 cout << "i"
276 << " " << i << "/" << nQuadraturePts << endl;
277 }
278
279 xcher = x_QuadraturePts[i];
280 ycher = y_QuadraturePts[i];
281
282 scale = sumd;
283 xin = xcher;
284 rex = 0.5 * pow(((m_Re) / scale), 2) + (m_Re)*xin;
285 delsx = sqrt(2.0 / rex) * scale * (xin)*m_Pr;
286 scale = scale / delsx;
287 delta = 4.91 * sqrt((xin * m_mu) / (m_rhoInf * m_uInf));
288 scale2 = ycher * (scale * delta) / sqrt(etamax);
289 coeff = 0.5 * sqrt(2 / (xcher * m_Re));
290
291 if (scale2 > z[m_xpoints - 3])
292 {
293 u_QuadraturePts[i] = 1;
294 rho_QuadraturePts[i] = 1;
295 T_QuadraturePts[i] = 1.0 / rho_QuadraturePts[i];
296 v_QuadraturePts[i] =
297 coeff * (z[m_xpoints - 3] - velocity[m_xpoints - 3]);
298
299 file3 << xcher << " " << ycher << " "
300 << velocity[m_xpoints - 3] << " " << z[m_xpoints - 3]
301 << " " << u[m_xpoints - 3] << endl;
302 }
303 else
304 {
305 for (int j = 0; j < m_xpoints - 1; j++)
306 {
307 if ((z[j] <= scale2) && (z[j + 1] > scale2))
308 {
309 index = j;
310 break;
311 }
312 }
313 if (index == -1)
314 {
315 ASSERTL0(false, "Could not determine index in CompressibleBL");
316 }
317
318 u_QuadraturePts[i] = u[index];
319 rho_QuadraturePts[i] = rho[index];
320 T_QuadraturePts[i] = 1.0 / rho_QuadraturePts[i];
321 v_QuadraturePts[i] = coeff * (u[index] * scale2 - velocity[index]);
322 }
323 }
324}
COMPBL
void COMPBL(Array< OneD, NekDouble > v, Array< OneD, NekDouble > dv)
Definition: CompressibleBL.cpp:94
Nektar::UnitTests::z
std::vector< double > z(NPUPPER)

References ASSERTL0, COMPBL(), etamax, m_mu, m_Pr, m_Re, m_rhoInf, m_Suth, m_uInf, m_xpoints, tinysimd::sqrt(), Nektar::UnitTests::test(), and Nektar::UnitTests::z().

Referenced by main().

◆ RK4()

void RK4 ( Array< OneD, NekDouble >  y,
Array< OneD, NekDouble >  dydx,
int  n,
NekDouble  x,
NekDouble  h,
Array< OneD, NekDouble >  yout 
)

Perform the RK4 integration

Definition at line 128 of file CompressibleBL.cpp.

130{
131 int nmax = 5;
132
133 Array<OneD, NekDouble> yt(nmax, 0.0);
134 Array<OneD, NekDouble> dyt(nmax, 0.0);
135 Array<OneD, NekDouble> dym(nmax, 0.0);
136 NekDouble hh = h * 0.5;
137 NekDouble h6 = h / 6;
138
139 for (int i = 0; i < n; i++)
140 {
141 yt[i] = y[i] + hh * dydx[i];
142 }
143
144 COMPBL(yt, dyt);
145
146 for (int i = 0; i < n; i++)
147 {
148 yt[i] = y[i] + hh * dyt[i];
149 }
150
151 COMPBL(yt, dym);
152
153 for (int i = 0; i < n; i++)
154 {
155 yt[i] = y[i] + h * dym[i];
156 dym[i] = dyt[i] + dym[i];
157 }
158
159 COMPBL(yt, dyt);
160
161 for (int i = 0; i < n; i++)
162 {
163 yout[i] = y[i] + h6 * (dydx[i] + dyt[i] + 2 * dym[i]);
164 }
165}

References COMPBL().

Referenced by RKDUMB().

◆ 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 
)

Calculate initial guess for RK4

Definition at line 170 of file CompressibleBL.cpp.

173{
174 int nmax = 5;
175 NekDouble x, h;
176 Array<OneD, NekDouble> v(nmax, 0.0);
177 Array<OneD, NekDouble> dv(nmax, 0.0);
178
179 for (int i = 0; i < nvar; i++)
180 {
181 v[i] = vstart[i];
182 y[i][0] = v[i];
183 }
184
185 xx[0] = x1;
186 x = x1;
187 h = (x2 - x1) / m_xpoints;
188
189 for (int k = 0; k < m_xpoints; k++)
190 {
191 COMPBL(v, dv);
192 RK4(v, dv, nvar, x, h, v);
193
194 if (x + h == x)
195 {
196 cout << "bug" << endl;
197 }
198
199 x = x + h;
200 xx[k + 1] = x;
201
202 for (int i = 0; i < nvar; i++)
203 {
204 y[i][k + 1] = v[i];
205 }
206 }
207}
RK4
void RK4(Array< OneD, NekDouble > y, Array< OneD, NekDouble > dydx, int n, NekDouble x, NekDouble h, Array< OneD, NekDouble > yout)
Definition: CompressibleBL.cpp:128

References COMPBL(), m_xpoints, and RK4().

Referenced by main().

Variable Documentation

◆ errtol

const NekDouble errtol = 1e-5

Definition at line 89 of file CompressibleBL.cpp.

Referenced by main().

◆ etamax

const NekDouble etamax = 10.0

Definition at line 88 of file CompressibleBL.cpp.

Referenced by main(), and OUTPUT().

◆ L

NekDouble L

Definition at line 69 of file CompressibleBL.cpp.

Referenced by Nektar::LibUtilities::FractionalInTimeIntegrationScheme::computeL(), Nektar::LibUtilities::FractionalInTimeIntegrationScheme::computeQML(), Nektar::LibUtilities::FractionalInTimeIntegrationScheme::computeTaus(), main(), and Nektar::LibUtilities::FractionalInTimeIntegrationScheme::v_TimeIntegrate().

◆ m_Gamma

NekDouble m_Gamma

Definition at line 74 of file CompressibleBL.cpp.

Referenced by COMPBL(), and main().

◆ m_long

NekDouble m_long

Definition at line 76 of file CompressibleBL.cpp.

Referenced by main().

◆ m_Mach

NekDouble m_Mach

Definition at line 68 of file CompressibleBL.cpp.

Referenced by COMPBL(), and main().

◆ m_mu

NekDouble m_mu

Definition at line 81 of file CompressibleBL.cpp.

Referenced by Nektar::FieldUtils::ProcessWSS::GetViscosity(), main(), Nektar::MMFDiffusion::Morphogenesis(), and OUTPUT().

◆ m_Pr

NekDouble m_Pr

Definition at line 75 of file CompressibleBL.cpp.

Referenced by COMPBL(), main(), and OUTPUT().

◆ m_R

NekDouble m_R

Definition at line 79 of file CompressibleBL.cpp.

Referenced by main().

◆ m_Re

NekDouble m_Re

Definition at line 67 of file CompressibleBL.cpp.

Referenced by main(), and OUTPUT().

◆ m_rhoInf

NekDouble m_rhoInf

Definition at line 78 of file CompressibleBL.cpp.

Referenced by Nektar::FieldUtils::ProcessWSS::GetViscosity(), main(), and OUTPUT().

◆ m_Suth

NekDouble m_Suth

Definition at line 71 of file CompressibleBL.cpp.

Referenced by COMPBL(), main(), and OUTPUT().

◆ m_Tinf

NekDouble m_Tinf

Definition at line 70 of file CompressibleBL.cpp.

Referenced by main().

◆ m_To

NekDouble m_To = 273.11

Definition at line 82 of file CompressibleBL.cpp.

◆ m_Tw

NekDouble m_Tw

Definition at line 72 of file CompressibleBL.cpp.

Referenced by main().

◆ m_Twall

NekDouble m_Twall

Definition at line 73 of file CompressibleBL.cpp.

Referenced by COMPBL(), and main().

◆ m_uInf

NekDouble m_uInf

Definition at line 77 of file CompressibleBL.cpp.

Referenced by main(), and OUTPUT().

◆ m_vInf

NekDouble m_vInf

Definition at line 80 of file CompressibleBL.cpp.

Referenced by main().

◆ m_xpoints

const int m_xpoints = 1000001

Definition at line 84 of file CompressibleBL.cpp.

Referenced by main(), OUTPUT(), and RKDUMB().

◆ Nvisc

const NekDouble Nvisc = 1

Definition at line 86 of file CompressibleBL.cpp.

Referenced by COMPBL().

◆ Omega

const NekDouble Omega = 1

Definition at line 87 of file CompressibleBL.cpp.

Referenced by COMPBL().

Generated on Wed Oct 30 2024 16:07:35 for Nektar++ by doxygen 1.9.4