1## libCEED: Navier-Stokes Example 2 3This page provides a description of the Navier-Stokes example for the libCEED library, based on PETSc. 4 5The Navier-Stokes problem solves the compressible Navier-Stokes equations in three dimensions using an explicit time integration. 6The state variables are mass density, momentum density, and energy density. 7 8The main Navier-Stokes solver for libCEED is defined in [`navierstokes.c`](navierstokes.c) with different problem definitions according to the application of interest. 9 10Build by using 11 12`make` 13 14and run with 15 16`./navierstokes` 17 18Available runtime options for all problem cases are: 19 20| Option | Meaning | Default | 21| :-------------------------------| :-----------------------------------------------------------------------------------| :-----------------------| 22| `-ceed` | CEED resource specifier | `/cpu/self/opt/blocked` | 23| `-test` | Run in test mode | `false` | 24| `-compare_final_state_atol` | Test absolute tolerance | `1E-11` | 25| `-compare_final_state_filename` | Test filename | | 26| `-problem` | Problem to solve (`advection`, `advection2d`, `density_current`, or `euler_vortex`) | `density_current` | 27| `-implicit` | Use implicit time integartor formulation | `false` (`explicit`) | 28| `-viz_refine` | Use regular refinement for visualization | `0` | 29| `-degree` | Polynomial degree of tensor product basis (must be >= 1) | `1` | 30| `-q_extra` | Number of extra quadrature points | `2` | 31| `-output_freq` | Frequency of output, in number of steps | `10` | 32| `-continue` | Continue from previous solution | `0` | 33| `-output_dir` | Output directory | `.` | 34 35### Advection 36 37This problem solves the convection (advection) equation for the total (scalar) energy density, transported by the (vector) velocity field. 38 39This is 3D advection given in two formulations based upon the weak form. 40 41State Variables: 42 43 *q = ( rho, U<sub>1</sub>, U<sub>2</sub>, U<sub>3</sub>, E )* 44 45 *rho* - Mass Density 46 47 *U<sub>i</sub>* - Momentum Density , *U<sub>i</sub> = rho ui* 48 49 *E* - Total Energy Density, *E = rho Cv T + rho (u u) / 2 + rho g z* 50 51Advection Equation: 52 53 *dE/dt + div( E _u_ ) = 0* 54 55#### Initial Conditions 56 57Mass Density: 58 Constant mass density of 1.0 59 60Momentum Density: 61 Rotational field in x,y with no momentum in z 62 63Energy Density: 64 Maximum of 1. x0 decreasing linearly to 0. as radial distance increases 65 to 1/8, then 0. everywhere else 66 67#### Boundary Conditions 68 69This problem is solved for two test cases with different BCs. 70 71##### Rotation 72 73Mass Density: 74 0.0 flux 75 76Momentum Density: 77 0.0 78 79Energy Density: 80 0.0 flux 81 82##### Translation 83 84Mass Density: 85 0.0 flux 86 87Momentum Density: 88 0.0 89 90Energy Density: 91 92Inflow BCs: 93 *E = E(wind)* 94 95Outflow BCs: 96 *E = E(boundary)* 97 98Both In/Outflow BCs for E are applied weakly. 99 100#### Runtime options for Advection 2D (-problem advection2d) 101 102| Option | Meaning | Default | Unit | 103| :-------------------| :--------------------------------------------------------| :----------| :----| 104| `-lx` | Length scale in x direction | `8000` | `m` | 105| `-ly` | Length scale in y direction | `8000` | `m` | 106| `-rc` | Characteristic radius of thermal bubble | `1000` | `m` | 107| `-units_meter` | 1 meter in scaled length units | `1E-2` | | 108| `-units_second` | 1 second in scaled time units | `1E-2` | | 109| `-units_kilogram` | 1 kilogram in scaled mass units | `1E-6` | | 110| `-strong_form` | Strong (1) or weak/integrated by parts (0) residual | `0` | | 111| `-stab` | Stabilization method (`none`, `su`, or `supg`) | `none` | | 112| `-CtauS` | Scale coefficient for stabilization tau (nondimensional) | `0` | | 113| `-wind_type` | Wind type in Advection (`rotation` or `translation`) | `rotation` | | 114| `-wind_translation` | Constant wind vector when `-wind_type translation` | `1,0,0` | | 115| `-E_wind` | Total energy of inflow wind when `-wind_type translation`| `1E6` | `J` | 116 117 118#### Runtime options for Advection 3D (-problem advection) 119 120| Option | Meaning | Default | Unit | 121| :--------------------| :---------------------------------------------------------| :----------| :----| 122| `-lx` | Length scale in x direction | `8000` | `m` | 123| `-ly` | Length scale in y direction | `8000` | `m` | 124| `-lz` | Length scale in z direction | `4000` | `m` | 125| `-rc` | Characteristic radius of thermal bubble | `1000` | `m` | 126| `-units_meter` | 1 meter in scaled length units | `1E-2` | | 127| `-units_second` | 1 second in scaled time units | `1E-2` | | 128| `-units_kilogram` | 1 kilogram in scaled mass units | `1E-6` | | 129| `-strong_form` | Strong (1) or weak/integrated by parts (0) residual | `0` | | 130| `-stab` | Stabilization method (`none`, `su`, or `supg`) | `none` | | 131| `-CtauS` | Scale coefficient for stabilization tau (nondimensional) | `0` | | 132| `-wind_type` | Wind type in Advection (`rotation` or `translation`) | `rotation` | | 133| `-wind_translation` | Constant wind vector when `-wind_type translation` | `1,0,0` | | 134| `-E_wind` | Total energy of inflow wind when `-wind_type translation` | `1E6` | `J` | 135| `-bubble_type` | `sphere` (3D) or `cylinder` (2D) | `shpere` | | 136| `-bubble_continuity` | `smooth`, `back_sharp`, or `thick` | `smooth` | | 137 138 139### Euler Traveling Vortex 140 141This problem solves the 3D Euler equations for vortex evolution provided in On the Order of Accuracy and Numerical Performance of Two Classes of Finite Volume WENO Schemes, Zhang, Zhang, and Shu (2011). 142 143State Variables: 144 145 *q = ( rho, U<sub>1</sub>, U<sub>2</sub>, U<sub>3</sub>, E )* 146 147 *rho* - Mass Density 148 149 *U<sub>i</sub>* - Momentum Density , *U<sub>i</sub> = rho u<sub>i</sub>* 150 151 *E* - Total Energy Density, *E = P / (gamma - 1) + rho (u u) / 2* 152 153Euler Equations: 154 155 *drho/dt + div( U ) = 0* 156 157 *dU/dt + div( rho (u x u) + P I<sub>3</sub> ) = 0* 158 159 *dE/dt + div( (E + P) u ) = 0* 160 161Constants: 162 163 *c<sub>v</sub>* , Specific heat, constant volume 164 165 *c<sub>p</sub>* , Specific heat, constant pressure 166 167 *gamma = c<sub>p</sub> / c<sub>v</sub>*, Specific heat ratio 168 169 *epsilon* , Vortex Strength 170 171#### Initial Conditions 172 173Temperature: 174 175 *T = 1 - (gamma - 1) epsilon^2 exp(1 - r^2) / (8 gamma pi^2)* 176 177Entropy: 178 179 *S = 1* , Constant entropy 180 181Density: 182 183 *rho = (T/S)^(1 / (gamma - 1))* 184 185Pressure: 186 187 *P = rho T* 188 189Velocity: 190 191 *u<sub>i</sub> = 1 + epsilon exp((1 - r^2)/2) [yc - y, x - xc, 0] / (2 pi)* 192 193 *r = sqrt( (x - xc)^2 + (y - yc)^2 )* 194 with *(xc,yc)* center of the xy-plane in the domain 195 196#### Boundary Conditions 197 198For this problem, in/outflow BCs are implemented where the validity of the weak 199form of the governing equations is extended to the outflow. 200For the inflow fluxes, prescribed T_inlet and P_inlet are converted to 201conservative variables and applied weakly. 202 203#### Runtime options for Euler Traveling Vortex (-problem euler_vortex) 204 205| Option | Meaning | Default | Unit | 206| :------------------| :----------------------------------------| :--------------| :---------| 207| `-lx` | Length scale in x direction | `1000` | `m` | 208| `-ly` | Length scale in y direction | `1000` | `m` | 209| `-lz` | Length scale in z direction | `1` | `m` | 210| `-center` | Location of vortex center | `(lx,ly,lz)/2` | `(m,m,m)` | 211| `-units_meter` | 1 meter in scaled length units | `1E-2` | | 212| `-units_second` | 1 second in scaled time units | `1E-2` | | 213| `-mean_velocity` | Background velocity vector | `(1,1,0)` | | 214| `-vortex_strength` | Strength of vortex | `5` | | 215| `-euler_test` | Euler test option (`t1`-`t4` and `none`) | `none` | | 216 217### Density Current 218 219This problem solves the full compressible Navier-Stokes equations, using operator composition and design of coupled solvers in the context of atmospheric modeling. 220This problem uses the formulation given in Semi-Implicit Formulations of the Navier-Stokes Equations: Application to Nonhydrostatic Atmospheric Modeling, Giraldo, Restelli, and Lauter (2010). 221 222The 3D compressible Navier-Stokes equations are formulated in conservation form with state variables of density, momentum density, and total energy density. 223 224State Variables: 225 226 *q = ( rho, U<sub>1</sub>, U<sub>2</sub>, U<sub>3</sub>, E )* 227 228 *rho* - Mass Density 229 230 *U<sub>i</sub>* - Momentum Density , *U<sub>i</sub> = rho u<sub>i</sub>* 231 232 *E* - Total Energy Density, *E = rho c<sub>v</sub> T + rho (u u) / 2 + rho g z* 233 234Navier-Stokes Equations: 235 236 *drho/dt + div( U ) = 0* 237 238 *dU/dt + div( rho (u x u) + P I<sub>3</sub> ) + rho g khat = div( F<sub>u</sub> )* 239 240 *dE/dt + div( (E + P) u ) = div( F<sub>e</sub> )* 241 242Viscous Stress: 243 244 *F<sub>u</sub> = mu (grad( u ) + grad( u )^T + lambda div ( u ) I<sub>3</sub>)* 245 246Thermal Stress: 247 248 *F<sub>e</sub> = u F<sub>u</sub> + k grad( T )* 249 250Equation of State: 251 252 *P = (gamma - 1) (E - rho (u u) / 2 - rho g z)* 253 254Temperature: 255 256 *T = (E / rho - (u u) / 2 - g z) / c<sub>v</sub>* 257 258Constants: 259 260 *lambda = - 2 / 3*, From Stokes hypothesis 261 262 *mu* , Dynamic viscosity 263 264 *k* , Thermal conductivity 265 266 *c<sub>v</sub>* , Specific heat, constant volume 267 268 *c<sub>p</sub>* , Specific heat, constant pressure 269 270 *g* , Gravity 271 272 *gamma = c<sub>p</sub> / c<sub>v</sub>*, Specific heat ratio 273 274#### Initial Conditions 275 276Potential Temperature: 277 278 *theta = thetabar + deltatheta* 279 280 *thetabar = theta0 exp( N * * 2 z / g )* 281 282 *deltatheta = 283 r <= rc : theta0(1 + cos(pi r)) / 2 284 r > rc : 0* 285 286 *r = sqrt( (x - xc) * * 2 + (y - yc) * * 2 + (z - zc) * * 2 )* 287 with *(xc,yc,zc)* center of domain 288 289Exner Pressure: 290 291 *Pi = Pibar + deltaPi* 292 293 *Pibar = g * * 2 (exp( - N * * 2 z / g ) - 1) / (cp theta0 N * * 2)* 294 295 *deltaPi = 0* (hydrostatic balance) 296 297Velocity/Momentum Density: 298 299 *U<sub>i</sub> = u<sub>i</sub> = 0* 300 301Conversion to Conserved Variables: 302 303 *rho = P0 Pi**(c<sub>v</sub>/R<sub>d</sub>) / (R<sub>d</sub> theta)* 304 305 *E = rho (c<sub>v</sub> theta Pi + (u u)/2 + g z)* 306 307Constants: 308 309 *theta0* , Potential temperature constant 310 311 *thetaC* , Potential temperature perturbation 312 313 *P0* , Pressure at the surface 314 315 *N* , Brunt-Vaisala frequency 316 317 *c<sub>v</sub>* , Specific heat, constant volume 318 319 *c<sub>p</sub>* , Specific heat, constant pressure 320 321 *R<sub>d</sub>* = c<sub>p</sub> - c<sub>v</sub>, Specific heat difference 322 323 *g* , Gravity 324 325 *r<sub>c</sub>* , Characteristic radius of thermal bubble 326 327 *l<sub>x</sub>* , Characteristic length scale of domain in x 328 329 *l<sub>y</sub>* , Characteristic length scale of domain in y 330 331 *l<sub>z</sub>* , Characteristic length scale of domain in z 332 333 334#### Boundary Conditions 335 336Mass Density: 337 0.0 flux 338 339Momentum Density: 340 0.0 341 342Energy Density: 343 0.0 flux 344 345#### Runtime options for Density Current (-problem density_current) 346 347| Option | Meaning | Default | Unit | 348| :-----------------| :-----------------------------------------------------------------------------------| :--------------| :----------| 349| `-lx` | Length scale in x direction | `8000` | `m` | 350| `-ly` | Length scale in y direction | `8000` | `m` | 351| `-lz` | Length scale in z direction | `4000` | `m` | 352| `-center` | Location of bubble center | `(lx,ly,lz)/2` | `(m,m,m)` | 353| `-dc_axis` | Axis of density current cylindrical anomaly, or `(0,0,0)` for spherically symmetric | `(0,0,0)` | | 354| `-rc` | Characteristic radius of thermal bubble | `1000` | `m` | 355| `-bc_wall` | Use wall boundary conditions on this list of faces | `-` | | 356| `-bc_slip_x` | Use slip boundary conditions, for the x component, on this list of faces | `5,6` | | 357| `-bc_slip_y` | Use slip boundary conditions, for the y component, on this list of faces | `3,4` | | 358| `-bc_slip_z` | Use slip boundary conditions, for the z component, on this list of faces | `1,2` | | 359| `-units_meter` | 1 meter in scaled length units | `1E-2` | | 360| `-units_second` | 1 second in scaled time units | `1E-2` | | 361| `-units_kilogram` | 1 kilogram in scaled mass units | `1E-6` | | 362| `-units_Kelvin` | 1 Kelvin in scaled temperature units | `1` | | 363| `-stab` | Stabilization method (`none`, `su`, or `supg`) | `none` | | 364| `-theta0` | Reference potential temperature | `300` | `K` | 365| `-thetaC` | Perturbation of potential temperature | `-15` | `K` | 366| `-P0` | Atmospheric pressure | `1E5` | `Pa` | 367| `-N` | Brunt-Vaisala frequency | `0.01` | `1/s` | 368| `-cv` | Heat capacity at constant volume | `717` | `J/(kg K)` | 369| `-cp` | Heat capacity at constant pressure | `1004 ` | `J/(kg K)` | 370| `-g` | Gravitational acceleration | `9.81` | `m/s^2` | 371| `-lambda` | Stokes hypothesis second viscosity coefficient | `-2/3` | | 372| `-mu` | Shear dynamic viscosity coefficient | `75` | `Pa s` | 373| `-k` | Thermal conductivity | `0.02638` | `W/(m K)` | 374 375For the case of a square/cubic mesh, the list of face indices to be used with `-bc_wall` and/or `-bc_slip_x`, `-bc_slip_y`, and `-bc_slip_z` are: 376 377* 2D: 378 - faceMarkerBottom = 1; 379 - faceMarkerRight = 2; 380 - faceMarkerTop = 3; 381 - faceMarkerLeft = 4; 382* 3D: 383 - faceMarkerBottom = 1; 384 - faceMarkerTop = 2; 385 - faceMarkerFront = 3; 386 - faceMarkerBack = 4; 387 - faceMarkerRight = 5; 388 - faceMarkerLeft = 6; 389 390### Time Discretization 391 392For all different problems, the time integration is performed with an explicit or implicit formulation. 393 394### Space Discretization 395 396The geometric factors and coordinate transformations required for the integration of the weak form are described in the files [`qfunctions/setupgeo.h`](qfunctions/setupgeo.h) and [`qfunctions/setupgeo2d.h`](qfunctions/setupgeo2d.h). 397