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