1## libCEED: Navier-Stokes Example 2 3This page provides a description of the Navier-Stokes example for the libCEED library, based on PETSc. 4PETSc v3.17 or a development version of PETSc at commit 0e95d842 or later is required. 5 6The Navier-Stokes problem solves the compressible Navier-Stokes equations in three dimensions using an explicit time integration. 7The state variables are mass density, momentum density, and energy density. 8 9The main Navier-Stokes solver for libCEED is defined in [`navierstokes.c`](navierstokes.c) with different problem definitions according to the application of interest. 10 11Build by using: 12 13`make` 14 15and run with: 16 17``` 18./navierstokes -ceed [ceed] -problem [problem type] -degree [degree] 19``` 20 21## Runtime options 22 23% inclusion-fluids-marker 24 25The Navier-Stokes mini-app is controlled via command-line options. 26The following options are common among all problem types: 27 28:::{list-table} Common Runtime Options 29:header-rows: 1 30 31* - Option 32 - Description 33 - Default value 34 35* - `-ceed` 36 - CEED resource specifier 37 - `/cpu/self/opt/blocked` 38 39* - `-test` 40 - Run in test mode 41 - `false` 42 43* - `-compare_final_state_atol` 44 - Test absolute tolerance 45 - `1E-11` 46 47* - `-compare_final_state_filename` 48 - Test filename 49 - 50 51* - `-problem` 52 - Problem to solve (`advection`, `advection2d`, `density_current`, or `euler_vortex`) 53 - `density_current` 54 55* - `-implicit` 56 - Use implicit time integartor formulation 57 - 58 59* - `-degree` 60 - Polynomial degree of tensor product basis (must be >= 1) 61 - `1` 62 63* - `-q_extra` 64 - Number of extra quadrature points 65 - `0` 66 67* - `-ts_monitor_solution` 68 - PETSc output format, such as `cgns:output-%d.cgns` (requires PETSc `--download-cgns`) 69 - 70 71* - `-ts_monitor_solution_interval` 72 - Number of time steps between visualization output frames. 73 - `1` 74 75* - `-viewer_cgns_batch_size` 76 - Number of frames written per CGNS file if the CGNS file name includes a format specifier (`%d`). 77 - `20` 78 79* - `-checkpoint_interval` 80 - Number of steps between writing binary checkpoints. `0` has no output, `-1` outputs final state only 81 - `10` 82 83* - `-checkpoint_vtk` 84 - Checkpoints include VTK (`*.vtu`) files for visualization. Consider `-ts_monitor_solution`instead. 85 - `false` 86 87* - `-viz_refine` 88 - Use regular refinement for VTK visualization 89 - `0` 90 91* - `-output_dir` 92 - Output directory for binary checkpoints and VTK files (if enabled). 93 - `.` 94 95* - `-output_add_stepnum2bin` 96 - Whether to add step numbers to output binary files 97 - `false` 98 99* - `-continue` 100 - Continue from previous solution (input is step number of previous solution) 101 - `0` 102 103* - `-continue_filename` 104 - Path to solution binary file from which to continue from 105 - `[output_dir]/ns-solution.bin` 106 107* - `-continue_time_filename` 108 - Path to time stamp binary file (only for legacy checkpoints) 109 - `[output_dir]/ns-time.bin` 110 111* - `-bc_wall` 112 - Use wall boundary conditions on this list of faces 113 - 114 115* - `-wall_comps` 116 - An array of constrained component numbers for wall BCs 117 - 118 119* - `-bc_slip_x` 120 - Use slip boundary conditions, for the x component, on this list of faces 121 - 122 123* - `-bc_slip_y` 124 - Use slip boundary conditions, for the y component, on this list of faces 125 - 126 127* - `-bc_slip_z` 128 - Use slip boundary conditions, for the z component, on this list of faces 129 - 130 131* - `-bc_inflow` 132 - Use inflow boundary conditions on this list of faces 133 - 134 135* - `-bc_outflow` 136 - Use outflow boundary conditions on this list of faces 137 - 138 139* - `-bc_freestream` 140 - Use freestream boundary conditions on this list of faces 141 - 142 143* - `-snes_view` 144 - View PETSc `SNES` nonlinear solver configuration 145 - 146 147* - `-log_view` 148 - View PETSc performance log 149 - 150 151* - `-help` 152 - View comprehensive information about run-time options 153 - 154::: 155 156For the case of a square/cubic mesh, the list of face indices to be used with `-bc_wall`, `bc_inflow`, `bc_outflow`, `bc_freestream` and/or `-bc_slip_x`, `-bc_slip_y`, and `-bc_slip_z` are: 157 158:::{list-table} 2D Face ID Labels 159:header-rows: 1 160* - PETSc Face Name 161 - Cartesian direction 162 - Face ID 163 164* - faceMarkerBottom 165 - -z 166 - 1 167 168* - faceMarkerRight 169 - +x 170 - 2 171 172* - faceMarkerTop 173 - +z 174 - 3 175 176* - faceMarkerLeft 177 - -x 178 - 4 179::: 180 181:::{list-table} 2D Face ID Labels 182:header-rows: 1 183* - PETSc Face Name 184 - Cartesian direction 185 - Face ID 186 187* - faceMarkerBottom 188 - -z 189 - 1 190 191* - faceMarkerTop 192 - +z 193 - 2 194 195* - faceMarkerFront 196 - -y 197 - 3 198 199* - faceMarkerBack 200 - +y 201 - 4 202 203* - faceMarkerRight 204 - +x 205 - 5 206 207* - faceMarkerLeft 208 - -x 209 - 6 210::: 211 212### Advection 213 214For testing purposes, there is a reduced mode for pure advection, which holds density $\rho$ and momentum density $\rho \bm u$ constant while advecting "total energy density" $E$. 215These are available in 2D and 3D. 216 217#### 2D advection 218 219For the 2D advection problem, the following additional command-line options are available: 220 221:::{list-table} Advection2D Runtime Options 222:header-rows: 1 223 224* - Option 225 - Description 226 - Default value 227 - Unit 228 229* - `-rc` 230 - Characteristic radius of thermal bubble 231 - `1000` 232 - `m` 233 234* - `-units_meter` 235 - 1 meter in scaled length units 236 - `1E-2` 237 - 238 239* - `-units_second` 240 - 1 second in scaled time units 241 - `1E-2` 242 - 243 244* - `-units_kilogram` 245 - 1 kilogram in scaled mass units 246 - `1E-6` 247 - 248 249* - `-strong_form` 250 - Strong (1) or weak/integrated by parts (0) residual 251 - `0` 252 - 253 254* - `-stab` 255 - Stabilization method (`none`, `su`, or `supg`) 256 - `none` 257 - 258 259* - `-CtauS` 260 - Scale coefficient for stabilization tau (nondimensional) 261 - `0` 262 - 263 264* - `-wind_type` 265 - Wind type in Advection (`rotation` or `translation`) 266 - `rotation` 267 - 268 269* - `-wind_translation` 270 - Constant wind vector when `-wind_type translation` 271 - `1,0,0` 272 - 273 274* - `-E_wind` 275 - Total energy of inflow wind when `-wind_type translation` 276 - `1E6` 277 - `J` 278::: 279 280An example of the `rotation` mode can be run with: 281 282``` 283./navierstokes -problem advection2d -dm_plex_box_faces 20,20 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1000,1000 -bc_wall 1,2,3,4 -wall_comps 4 -wind_type rotation -implicit -stab supg 284``` 285 286and the `translation` mode with: 287 288``` 289./navierstokes -problem advection2d -dm_plex_box_faces 20,20 -dm_plex_box_lower 0,0 -dm_plex_box_upper 1000,1000 -units_meter 1e-4 -wind_type translation -wind_translation 1,-.5 -bc_inflow 1,2,3,4 290``` 291Note the lengths in `-dm_plex_box_upper` are given in meters, and will be nondimensionalized according to `-units_meter`. 292 293#### 3D advection 294 295For the 3D advection problem, the following additional command-line options are available: 296 297:::{list-table} Advection3D Runtime Options 298:header-rows: 1 299 300* - Option 301 - Description 302 - Default value 303 - Unit 304 305* - `-rc` 306 - Characteristic radius of thermal bubble 307 - `1000` 308 - `m` 309 310* - `-units_meter` 311 - 1 meter in scaled length units 312 - `1E-2` 313 - 314 315* - `-units_second` 316 - 1 second in scaled time units 317 - `1E-2` 318 - 319 320* - `-units_kilogram` 321 - 1 kilogram in scaled mass units 322 - `1E-6` 323 - 324 325* - `-strong_form` 326 - Strong (1) or weak/integrated by parts (0) residual 327 - `0` 328 - 329 330* - `-stab` 331 - Stabilization method (`none`, `su`, or `supg`) 332 - `none` 333 - 334 335* - `-CtauS` 336 - Scale coefficient for stabilization tau (nondimensional) 337 - `0` 338 - 339 340* - `-wind_type` 341 - Wind type in Advection (`rotation` or `translation`) 342 - `rotation` 343 - 344 345* - `-wind_translation` 346 - Constant wind vector when `-wind_type translation` 347 - `1,0,0` 348 - 349 350* - `-E_wind` 351 - Total energy of inflow wind when `-wind_type translation` 352 - `1E6` 353 - `J` 354 355* - `-bubble_type` 356 - `sphere` (3D) or `cylinder` (2D) 357 - `shpere` 358 - 359 360* - `-bubble_continuity` 361 - `smooth`, `back_sharp`, or `thick` 362 - `smooth` 363 - 364::: 365 366An example of the `rotation` mode can be run with: 367 368``` 369./navierstokes -problem advection -dm_plex_box_faces 10,10,10 -dm_plex_dim 3 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 8000,8000,8000 -bc_wall 1,2,3,4,5,6 -wall_comps 4 -wind_type rotation -implicit -stab su 370``` 371 372and the `translation` mode with: 373 374``` 375./navierstokes -problem advection -dm_plex_box_faces 10,10,10 -dm_plex_dim 3 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 8000,8000,8000 -wind_type translation -wind_translation .5,-1,0 -bc_inflow 1,2,3,4,5,6 376``` 377 378### Inviscid Ideal Gas 379 380#### Isentropic Euler vortex 381 382For the Isentropic Vortex problem, the following additional command-line options are available: 383 384:::{list-table} Isentropic Vortex Runtime Options 385:header-rows: 1 386 387* - Option 388 - Description 389 - Default value 390 - Unit 391 392* - `-center` 393 - Location of vortex center 394 - `(lx,ly,lz)/2` 395 - `(m,m,m)` 396 397* - `-units_meter` 398 - 1 meter in scaled length units 399 - `1E-2` 400 - 401 402* - `-units_second` 403 - 1 second in scaled time units 404 - `1E-2` 405 - 406 407* - `-mean_velocity` 408 - Background velocity vector 409 - `(1,1,0)` 410 - 411 412* - `-vortex_strength` 413 - Strength of vortex < 10 414 - `5` 415 - 416 417* - `-c_tau` 418 - Stabilization constant 419 - `0.5` 420 - 421::: 422 423This problem can be run with: 424 425``` 426./navierstokes -problem euler_vortex -dm_plex_box_faces 20,20,1 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1000,1000,50 -dm_plex_dim 3 -bc_inflow 4,6 -bc_outflow 3,5 -bc_slip_z 1,2 -mean_velocity .5,-.8,0. 427``` 428 429#### Sod shock tube 430 431For the Shock Tube problem, the following additional command-line options are available: 432 433:::{list-table} Shock Tube Runtime Options 434:header-rows: 1 435 436* - Option 437 - Description 438 - Default value 439 - Unit 440 441* - `-units_meter` 442 - 1 meter in scaled length units 443 - `1E-2` 444 - 445 446* - `-units_second` 447 - 1 second in scaled time units 448 - `1E-2` 449 - 450 451* - `-yzb` 452 - Use YZB discontinuity capturing 453 - `none` 454 - 455 456* - `-stab` 457 - Stabilization method (`none`, `su`, or `supg`) 458 - `none` 459 - 460::: 461 462This problem can be run with: 463 464``` 465./navierstokes -problem shocktube -yzb -stab su -bc_slip_z 3,4 -bc_slip_y 1,2 -bc_wall 5,6 -dm_plex_dim 3 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 1000,100,100 -dm_plex_box_faces 200,1,1 -units_second 0.1 466``` 467 468### Newtonian viscosity, Ideal Gas 469 470For the Density Current, Channel, and Blasius problems, the following common command-line options are available: 471 472:::{list-table} Newtonian Ideal Gas problems Runtime Options 473:header-rows: 1 474 475* - Option 476 - Description 477 - Default value 478 - Unit 479 480* - `-units_meter` 481 - 1 meter in scaled length units 482 - `1` 483 - 484 485* - `-units_second` 486 - 1 second in scaled time units 487 - `1` 488 - 489 490* - `-units_kilogram` 491 - 1 kilogram in scaled mass units 492 - `1` 493 - 494 495* - `-units_Kelvin` 496 - 1 Kelvin in scaled temperature units 497 - `1` 498 - 499 500* - `-stab` 501 - Stabilization method (`none`, `su`, or `supg`) 502 - `none` 503 - 504 505* - `-c_tau` 506 - Stabilization constant, $c_\tau$ 507 - `0.5` 508 - 509 510* - `-Ctau_t` 511 - Stabilization time constant, $C_t$ 512 - `1.0` 513 - 514 515* - `-Ctau_v` 516 - Stabilization viscous constant, $C_v$ 517 - `36, 60, 128 for degree = 1, 2, 3` 518 - 519 520* - `-Ctau_C` 521 - Stabilization continuity constant, $C_c$ 522 - `1.0` 523 - 524 525* - `-Ctau_M` 526 - Stabilization momentum constant, $C_m$ 527 - `1.0` 528 - 529 530* - `-Ctau_E` 531 - Stabilization energy constant, $C_E$ 532 - `1.0` 533 - 534 535* - `-cv` 536 - Heat capacity at constant volume 537 - `717` 538 - `J/(kg K)` 539 540* - `-cp` 541 - Heat capacity at constant pressure 542 - `1004` 543 - `J/(kg K)` 544 545* - `-g` 546 - Gravitational acceleration 547 - `9.81` 548 - `m/s^2` 549 550* - `-lambda` 551 - Stokes hypothesis second viscosity coefficient 552 - `-2/3` 553 - 554 555* - `-mu` 556 - Shear dynamic viscosity coefficient 557 - `75` 558 - `Pa s` 559 560* - `-k` 561 - Thermal conductivity 562 - `0.02638` 563 - `W/(m K)` 564 565* - `-newtonian_unit_tests` 566 - Developer option to test properties 567 - `false` 568 - boolean 569 570* - `-state_var` 571 - State variables to solve solution with. `conservative` ($\rho, \rho \bm{u}, \rho e$) or `primitive` ($P, \bm{u}, T$) 572 - `conservative` 573 - string 574::: 575 576#### Newtonian Wave 577 578The newtonian wave problem has the following command-line options in addition to the Newtonian Ideal Gas options: 579 580:::{list-table} Newtonian Wave Runtime Options 581:header-rows: 1 582 583* - Option 584 - Description 585 - Default value 586 - Unit 587 588* - `-freestream_riemann` 589 - Riemann solver for boundaries (HLL or HLLC) 590 - `hllc` 591 - 592 593* - `-freestream_velocity` 594 - Freestream velocity vector 595 - `0,0,0` 596 - `m/s` 597 598* - `-freestream_temperature` 599 - Freestream temperature 600 - `288` 601 - `K` 602 603* - `-freestream_pressure` 604 - Freestream pressure 605 - `1.01e5` 606 - `Pa` 607 608* - `-epicenter` 609 - Coordinates of center of perturbation 610 - `0,0,0` 611 - `m` 612 613* - `-amplitude` 614 - Amplitude of the perturbation 615 - `0.1` 616 - 617 618* - `-width` 619 - Width parameter of the perturbation 620 - `0.002` 621 - `m` 622 623::: 624 625This problem can be run with the `newtonianwave.yaml` file via: 626 627``` 628./navierstokes -options_file newtonianwave.yaml 629``` 630 631```{literalinclude} ../../../../../examples/fluids/newtonianwave.yaml 632:language: yaml 633``` 634 635#### Vortex Shedding - Flow past Cylinder 636 637The vortex shedding, flow past cylinder problem has the following command-line options in addition to the Newtonian Ideal Gas options: 638 639:::{list-table} Vortex Shedding Runtime Options 640:header-rows: 1 641 642* - Option 643 - Description 644 - Default value 645 - Unit 646 647* - `-freestream_velocity` 648 - Freestream velocity vector 649 - `0,0,0` 650 - `m/s` 651 652* - `-freestream_temperature` 653 - Freestream temperature 654 - `288` 655 - `K` 656 657* - `-freestream_pressure` 658 - Freestream pressure 659 - `1.01e5` 660 - `Pa` 661 662::: 663 664The initial condition is taken from `-reference_temperature` and `-reference_pressure`. 665To run this problem, first generate a mesh: 666 667```console 668$ make -C examples/fluids/meshes 669``` 670 671Then run by building the executable and running: 672 673```console 674$ make build/fluids-navierstokes 675$ mpiexec -n 6 build/fluids-navierstokes -options_file vortexshedding.yaml 676``` 677 678The vortex shedding period is roughly 6 and this problem runs until time 100 (2000 time steps). 679 680```{literalinclude} ../../../../../examples/fluids/vortexshedding.yaml 681:language: yaml 682``` 683 684#### Density current 685 686The Density Current problem has the following command-line options in addition to the Newtonian Ideal Gas options: 687 688:::{list-table} Density Current Runtime Options 689:header-rows: 1 690 691* - Option 692 - Description 693 - Default value 694 - Unit 695 696* - `-center` 697 - Location of bubble center 698 - `(lx,ly,lz)/2` 699 - `(m,m,m)` 700 701* - `-dc_axis` 702 - Axis of density current cylindrical anomaly, or `(0,0,0)` for spherically symmetric 703 - `(0,0,0)` 704 - 705 706* - `-rc` 707 - Characteristic radius of thermal bubble 708 - `1000` 709 - `m` 710 711* - `-theta0` 712 - Reference potential temperature 713 - `300` 714 - `K` 715 716* - `-thetaC` 717 - Perturbation of potential temperature 718 - `-15` 719 - `K` 720 721* - `-P0` 722 - Atmospheric pressure 723 - `1E5` 724 - `Pa` 725 726* - `-N` 727 - Brunt-Vaisala frequency 728 - `0.01` 729 - `1/s` 730::: 731 732This problem can be run with: 733 734``` 735./navierstokes -problem density_current -dm_plex_box_faces 16,1,8 -degree 1 -dm_plex_box_lower 0,0,0 -dm_plex_box_upper 2000,125,1000 -dm_plex_dim 3 -rc 400. -bc_wall 1,2,5,6 -wall_comps 1,2,3 -bc_slip_y 3,4 -mu 75 736``` 737 738#### Channel flow 739 740The Channel problem has the following command-line options in addition to the Newtonian Ideal Gas options: 741 742:::{list-table} Channel Runtime Options 743:header-rows: 1 744 745* - Option 746 - Description 747 - Default value 748 - Unit 749 750* - `-umax` 751 - Maximum/centerline velocity of the flow 752 - `10` 753 - `m/s` 754 755* - `-theta0` 756 - Reference potential temperature 757 - `300` 758 - `K` 759 760* - `-P0` 761 - Atmospheric pressure 762 - `1E5` 763 - `Pa` 764 765* - `-body_force_scale` 766 - Multiplier for body force (`-1` for flow reversal) 767 - 1 768 - 769::: 770 771This problem can be run with the `channel.yaml` file via: 772 773``` 774./navierstokes -options_file channel.yaml 775``` 776```{literalinclude} ../../../../../examples/fluids/channel.yaml 777:language: yaml 778``` 779 780#### Blasius boundary layer 781 782The Blasius problem has the following command-line options in addition to the Newtonian Ideal Gas options: 783 784:::{list-table} Blasius Runtime Options 785:header-rows: 1 786 787* - Option 788 - Description 789 - Default value 790 - Unit 791 792* - `-velocity_infinity` 793 - Freestream velocity 794 - `40` 795 - `m/s` 796 797* - `-temperature_infinity` 798 - Freestream temperature 799 - `288` 800 - `K` 801 802* - `-temperature_wall` 803 - Wall temperature 804 - `288` 805 - `K` 806 807* - `-delta0` 808 - Boundary layer height at the inflow 809 - `4.2e-3` 810 - `m` 811 812* - `-P0` 813 - Atmospheric pressure 814 - `1.01E5` 815 - `Pa` 816 817* - `-platemesh_refine_height` 818 - Height at which `-platemesh_Ndelta` number of elements should refined into 819 - `5.9E-4` 820 - `m` 821 822* - `-platemesh_Ndelta` 823 - Number of elements to keep below `-platemesh_refine_height` 824 - `45` 825 - 826 827* - `-platemesh_growth` 828 - Growth rate of the elements in the refinement region 829 - `1.08` 830 - 831 832* - `-platemesh_top_angle` 833 - Downward angle of the top face of the domain. This face serves as an outlet. 834 - `5` 835 - `degrees` 836 837* - `-stg_use` 838 - Whether to use stg for the inflow conditions 839 - `false` 840 - 841 842* - `-platemesh_y_node_locs_path` 843 - Path to file with y node locations. If empty, will use mesh warping instead. 844 - `""` 845 - 846 847* - `-n_chebyshev` 848 - Number of Chebyshev terms 849 - `20` 850 - 851 852* - `-chebyshev_` 853 - Prefix for Chebyshev snes solve 854 - 855 - 856 857::: 858 859This problem can be run with the `blasius.yaml` file via: 860 861``` 862./navierstokes -options_file blasius.yaml 863``` 864 865```{literalinclude} ../../../../../examples/fluids/blasius.yaml 866:language: yaml 867``` 868 869#### STG Inflow for Flat Plate 870 871Using the STG Inflow for the blasius problem adds the following command-line options: 872 873:::{list-table} Blasius Runtime Options 874:header-rows: 1 875 876* - Option 877 - Description 878 - Default value 879 - Unit 880 881* - `-stg_inflow_path` 882 - Path to the STGInflow file 883 - `./STGInflow.dat` 884 - 885 886* - `-stg_rand_path` 887 - Path to the STGRand file 888 - `./STGRand.dat` 889 - 890 891* - `-stg_alpha` 892 - Growth rate of the wavemodes 893 - `1.01` 894 - 895 896* - `-stg_u0` 897 - Convective velocity, $U_0$ 898 - `0.0` 899 - `m/s` 900 901* - `-stg_mean_only` 902 - Only impose the mean velocity (no fluctutations) 903 - `false` 904 - 905 906* - `-stg_strong` 907 - Strongly enforce the STG inflow boundary condition 908 - `false` 909 - 910 911* - `-stg_fluctuating_IC` 912 - "Extrude" the fluctuations through the domain as an initial condition 913 - `false` 914 - 915 916::: 917 918This problem can be run with the `blasius.yaml` file via: 919 920``` 921./navierstokes -options_file blasius.yaml -stg_use true 922``` 923 924Note the added `-stg_use true` flag 925This overrides the `stg: use: false` setting in the `blasius.yaml` file, enabling the use of the STG inflow. 926