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 122Mass Density: 123 0.0 flux 124 125Momentum Density: 126 0.0 127 128Energy Density: 129 0.0 flux 130 131### Euler Traveling Vortex 132 133This problem solves the 3D Euler equations for vortex evolution provided 134in On the Order of Accuracy and Numerical Performance of Two Classes of 135Finite Volume WENO Schemes, Zhang, Zhang, and Shu (2011). 136 137State Variables: 138 139 *q = ( rho, U<sub>1</sub>, U<sub>2</sub>, U<sub>3</sub>, E )* 140 141 *rho* - Mass Density 142 143 *U<sub>i</sub>* - Momentum Density , *U<sub>i</sub> = rho u<sub>i</sub>* 144 145 *E* - Total Energy Density, *E = P / (gamma - 1) + rho (u u) / 2* 146 147Euler Equations: 148 149 *drho/dt + div( U ) = 0* 150 151 *dU/dt + div( rho (u x u) + P I<sub>3</sub> ) = 0* 152 153 *dE/dt + div( (E + P) u ) = 0* 154 155Constants: 156 157 *c<sub>v</sub>* , Specific heat, constant volume 158 159 *c<sub>p</sub>* , Specific heat, constant pressure 160 161 *gamma = c<sub>p</sub> / c<sub>v</sub>*, Specific heat ratio 162 163 *epsilon* , Vortex Strength 164 165#### Initial Conditions 166 167Temperature: 168 169 *T = 1 - (gamma - 1) epsilon^2 exp(1 - r^2) / (8 gamma pi^2)* 170 171Entropy: 172 173 *S = 1* , Constant entropy 174 175Density: 176 177 *rho = (T/S)^(1 / (gamma - 1))* 178 179Pressure: 180 181 *P = rho T* 182 183Velocity: 184 185 *u<sub>i</sub> = 1 + epsilon exp((1 - r^2)/2) [yc - y, x - xc, 0] / (2 pi)* 186 187 *r = sqrt( (x - xc)^2 + (y - yc)^2 )* 188 with *(xc,yc)* center of the xy-plane in the domain 189 190#### Boundary Conditions 191 192For this problem, in/outflow BCs are implemented where the validity of the weak 193form of the governing equations is extended to the outflow. 194For the inflow fluxes, prescribed T_inlet and P_inlet are converted to 195conservative variables and applied weakly. 196 197### Density Current 198 199This problem solves the full compressible Navier-Stokes equations, using 200operator composition and design of coupled solvers in the context of atmospheric 201modeling. This problem uses the formulation given in Semi-Implicit Formulations 202of the Navier-Stokes Equations: Application to Nonhydrostatic Atmospheric Modeling, 203Giraldo, Restelli, and Lauter (2010). 204 205The 3D compressible Navier-Stokes equations are formulated in conservation form with state 206variables of density, momentum density, and total energy density. 207 208State Variables: 209 210 *q = ( rho, U<sub>1</sub>, U<sub>2</sub>, U<sub>3</sub>, E )* 211 212 *rho* - Mass Density 213 214 *U<sub>i</sub>* - Momentum Density , *U<sub>i</sub> = rho u<sub>i</sub>* 215 216 *E* - Total Energy Density, *E = rho c<sub>v</sub> T + rho (u u) / 2 + rho g z* 217 218Navier-Stokes Equations: 219 220 *drho/dt + div( U ) = 0* 221 222 *dU/dt + div( rho (u x u) + P I<sub>3</sub> ) + rho g khat = div( F<sub>u</sub> )* 223 224 *dE/dt + div( (E + P) u ) = div( F<sub>e</sub> )* 225 226Viscous Stress: 227 228 *F<sub>u</sub> = mu (grad( u ) + grad( u )^T + lambda div ( u ) I<sub>3</sub>)* 229 230Thermal Stress: 231 232 *F<sub>e</sub> = u F<sub>u</sub> + k grad( T )* 233 234Equation of State: 235 236 *P = (gamma - 1) (E - rho (u u) / 2 - rho g z)* 237 238Temperature: 239 240 *T = (E / rho - (u u) / 2 - g z) / c<sub>v</sub>* 241 242Constants: 243 244 *lambda = - 2 / 3*, From Stokes hypothesis 245 246 *mu* , Dynamic viscosity 247 248 *k* , Thermal conductivity 249 250 *c<sub>v</sub>* , Specific heat, constant volume 251 252 *c<sub>p</sub>* , Specific heat, constant pressure 253 254 *g* , Gravity 255 256 *gamma = c<sub>p</sub> / c<sub>v</sub>*, Specific heat ratio 257 258#### Initial Conditions 259 260Potential Temperature: 261 262 *theta = thetabar + deltatheta* 263 264 *thetabar = theta0 exp( N * * 2 z / g )* 265 266 *deltatheta = 267 r <= rc : theta0(1 + cos(pi r)) / 2 268 r > rc : 0* 269 270 *r = sqrt( (x - xc) * * 2 + (y - yc) * * 2 + (z - zc) * * 2 )* 271 with *(xc,yc,zc)* center of domain 272 273Exner Pressure: 274 275 *Pi = Pibar + deltaPi* 276 277 *Pibar = g * * 2 (exp( - N * * 2 z / g ) - 1) / (cp theta0 N * * 2)* 278 279 *deltaPi = 0* (hydrostatic balance) 280 281Velocity/Momentum Density: 282 283 *U<sub>i</sub> = u<sub>i</sub> = 0* 284 285Conversion to Conserved Variables: 286 287 *rho = P0 Pi**(c<sub>v</sub>/R<sub>d</sub>) / (R<sub>d</sub> theta)* 288 289 *E = rho (c<sub>v</sub> theta Pi + (u u)/2 + g z)* 290 291Constants: 292 293 *theta0* , Potential temperature constant 294 295 *thetaC* , Potential temperature perturbation 296 297 *P0* , Pressure at the surface 298 299 *N* , Brunt-Vaisala frequency 300 301 *c<sub>v</sub>* , Specific heat, constant volume 302 303 *c<sub>p</sub>* , Specific heat, constant pressure 304 305 *R<sub>d</sub>* = c<sub>p</sub> - c<sub>v</sub>, Specific heat difference 306 307 *g* , Gravity 308 309 *r<sub>c</sub>* , Characteristic radius of thermal bubble 310 311 *l<sub>x</sub>* , Characteristic length scale of domain in x 312 313 *l<sub>y</sub>* , Characteristic length scale of domain in y 314 315 *l<sub>z</sub>* , Characteristic length scale of domain in z 316 317 318#### Boundary Conditions 319 320Mass Density: 321 0.0 flux 322 323Momentum Density: 324 0.0 325 326Energy Density: 327 0.0 flux 328 329### Time Discretization 330 331For all different problems, the time integration is performed with an explicit 332or implicit formulation. 333 334### Space Discretization 335 336The geometric factors and coordinate transformations required for the integration of the weak form 337are described in the file [`common.h`](common.h) 338