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`, or `density_current`) | 26| `-stab` | Stabilization method | 27| `-implicit` | Use implicit time integartor formulation | 28| `-bc_wall` | Use wall boundary conditions on this list of faces | 29| `-bc_outflow` | Use outflow boundary conditions on this list of faces | 30| `-bc_slip_x` | Use slip boundary conditions, for the x component, on this list of faces | 31| `-bc_slip_y` | Use slip boundary conditions, for the y component, on this list of faces | 32| `-bc_slip_z` | Use slip boundary conditions, for the z component, on this list of faces | 33| `-viz_refine` | Use regular refinement for visualization | 34| `-degree` | Polynomial degree of tensor product basis (must be >= 1) | 35| `-units_meter` | 1 meter in scaled length units | 36| `-units_second` | 1 second in scaled time units | 37| `-units_kilogram` | 1 kilogram in scaled mass units | 38| `-units_Kelvin` | 1 Kelvin in scaled temperature units | 39| `-theta0` | Reference potential temperature | 40| `-thetaC` | Perturbation of potential temperature | 41| `-P0` | Atmospheric pressure | 42| `-N` | Brunt-Vaisala frequency | 43| `-cv` | Heat capacity at constant volume | 44| `-cp` | Heat capacity at constant pressure | 45| `-g` | Gravitational acceleration | 46| `-lambda` | Stokes hypothesis second viscosity coefficient | 47| `-mu` | Shear dynamic viscosity coefficient | 48| `-k` | Thermal conductivity | 49| `-CtauS` | Scale coefficient for stabilization tau (nondimensional) | 50| `-strong_form` | Strong (1) or weak/integrated by parts (0) advection residual | 51| `-lx` | Length scale in x direction | 52| `-ly` | Length scale in y direction | 53| `-lz` | Length scale in z direction | 54| `-rc` | Characteristic radius of thermal bubble | 55| `-resx` | Resolution in x | 56| `-resy` | Resolution in y | 57| `-resz` | Resolution in z | 58| `-center` | Location of bubble center | 59| `-dc_axis` | Axis of density current cylindrical anomaly, or {0,0,0} for spherically symmetric | 60| `-output_freq` | Frequency of output, in number of steps | 61| `-continue` | Continue from previous solution | 62| `-degree` | Polynomial degree of tensor product basis | 63| `-qextra` | Number of extra quadrature points | 64| `-of` | Output folder | 65 66For the case of a square/cubic mesh, the list of face indices to be used with `-bc_wall` and/or `-bc_slip_x`, 67`-bc_slip_y`, and `-bc_slip_z` or `bc_outflow` are: 68 69* 2D: 70 - faceMarkerBottom = 1; 71 - faceMarkerRight = 2; 72 - faceMarkerTop = 3; 73 - faceMarkerLeft = 4; 74* 3D: 75 - faceMarkerBottom = 1; 76 - faceMarkerTop = 2; 77 - faceMarkerFront = 3; 78 - faceMarkerBack = 4; 79 - faceMarkerRight = 5; 80 - faceMarkerLeft = 6; 81 82### Advection 83 84This problem solves the convection (advection) equation for the total (scalar) energy density, 85transported by the (vector) velocity field. 86 87This is 3D advection given in two formulations based upon the weak form. 88 89State Variables: 90 91 *q = ( rho, U<sub>1</sub>, U<sub>2</sub>, U<sub>3</sub>, E )* 92 93 *rho* - Mass Density 94 95 *U<sub>i</sub>* - Momentum Density , *U<sub>i</sub> = rho ui* 96 97 *E* - Total Energy Density, *E = rho Cv T + rho (u u) / 2 + rho g z* 98 99Advection Equation: 100 101 *dE/dt + div( E _u_ ) = 0* 102 103#### Initial Conditions 104 105Mass Density: 106 Constant mass density of 1.0 107 108Momentum Density: 109 Rotational field in x,y with no momentum in z 110 111Energy Density: 112 Maximum of 1. x0 decreasing linearly to 0. as radial distance increases 113 to 1/8, then 0. everywhere else 114 115#### Boundary Conditions 116 117Mass Density: 118 0.0 flux 119 120Momentum Density: 121 0.0 122 123Energy Density: 124 0.0 flux 125 126### Density Current 127 128This problem solves the full compressible Navier-Stokes equations, using 129operator composition and design of coupled solvers in the context of atmospheric 130modeling. This problem uses the formulation given in Semi-Implicit Formulations 131of the Navier-Stokes Equations: Application to Nonhydrostatic Atmospheric Modeling, 132Giraldo, Restelli, and Lauter (2010). 133 134The 3D compressible Navier-Stokes equations are formulated in conservation form with state 135variables of density, momentum density, and total energy density. 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 = rho c<sub>v</sub> T + rho (u u) / 2 + rho g z* 146 147Navier-Stokes Equations: 148 149 *drho/dt + div( U ) = 0* 150 151 *dU/dt + div( rho (u x u) + P I<sub>3</sub> ) + rho g khat = div( F<sub>u</sub> )* 152 153 *dE/dt + div( (E + P) u ) = div( F<sub>e</sub> )* 154 155Viscous Stress: 156 157 *F<sub>u</sub> = mu (grad( u ) + grad( u )^T + lambda div ( u ) I<sub>3</sub>)* 158 159Thermal Stress: 160 161 *F<sub>e</sub> = u F<sub>u</sub> + k grad( T )* 162 163Equation of State: 164 165 *P = (gamma - 1) (E - rho (u u) / 2 - rho g z)* 166 167Temperature: 168 169 *T = (E / rho - (u u) / 2 - g z) / c<sub>v</sub>* 170 171Constants: 172 173 *lambda = - 2 / 3*, From Stokes hypothesis 174 175 *mu* , Dynamic viscosity 176 177 *k* , Thermal conductivity 178 179 *c<sub>v</sub>* , Specific heat, constant volume 180 181 *c<sub>p</sub>* , Specific heat, constant pressure 182 183 *g* , Gravity 184 185 *gamma = c<sub>p</sub> / c<sub>v</sub>*, Specific heat ratio 186 187#### Initial Conditions 188 189Potential Temperature: 190 191 *theta = thetabar + deltatheta* 192 193 *thetabar = theta0 exp( N * * 2 z / g )* 194 195 *deltatheta = 196 r <= rc : theta0(1 + cos(pi r)) / 2 197 r > rc : 0* 198 199 *r = sqrt( (x - xc) * * 2 + (y - yc) * * 2 + (z - zc) * * 2 )* 200 with *(xc,yc,zc)* center of domain 201 202Exner Pressure: 203 204 *Pi = Pibar + deltaPi* 205 206 *Pibar = g * * 2 (exp( - N * * 2 z / g ) - 1) / (cp theta0 N * * 2)* 207 208 *deltaPi = 0* (hydrostatic balance) 209 210Velocity/Momentum Density: 211 212 *U<sub>i</sub> = u<sub>i</sub> = 0* 213 214Conversion to Conserved Variables: 215 216 *rho = P0 Pi**(c<sub>v</sub>/R<sub>d</sub>) / (R<sub>d</sub> theta)* 217 218 *E = rho (c<sub>v</sub> theta Pi + (u u)/2 + g z)* 219 220Constants: 221 222 *theta0* , Potential temperature constant 223 224 *thetaC* , Potential temperature perturbation 225 226 *P0* , Pressure at the surface 227 228 *N* , Brunt-Vaisala frequency 229 230 *c<sub>v</sub>* , Specific heat, constant volume 231 232 *c<sub>p</sub>* , Specific heat, constant pressure 233 234 *R<sub>d</sub>* = c<sub>p</sub> - c<sub>v</sub>, Specific heat difference 235 236 *g* , Gravity 237 238 *r<sub>c</sub>* , Characteristic radius of thermal bubble 239 240 *l<sub>x</sub>* , Characteristic length scale of domain in x 241 242 *l<sub>y</sub>* , Characteristic length scale of domain in y 243 244 *l<sub>z</sub>* , Characteristic length scale of domain in z 245 246 247#### Boundary Conditions 248 249Mass Density: 250 0.0 flux 251 252Momentum Density: 253 0.0 254 255Energy Density: 256 0.0 flux 257 258### Time Discretization 259 260For all different problems, the time integration is performed with an explicit 261or implicit formulation. 262 263### Space Discretization 264 265The geometric factors and coordinate transformations required for the integration of the weak form 266are described in the file [`common.h`](common.h) 267