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