From: bangerth Date: Tue, 20 May 2008 16:38:57 +0000 (+0000) Subject: Add the forgotten forcing term to the documentation in various places. Make it generi... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9884a5a11cff115eb3f57caa7910ff70b2e1647a;p=dealii-svn.git Add the forgotten forcing term to the documentation in various places. Make it generic in the program as well, rather than adding it by hand to the cell assembly part. git-svn-id: https://svn.dealii.org/trunk@16137 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-33/doc/intro.dox b/deal.II/examples/step-33/doc/intro.dox index 79fdb7a1ee..ef696f4c25 100644 --- a/deal.II/examples/step-33/doc/intro.dox +++ b/deal.II/examples/step-33/doc/intro.dox @@ -28,19 +28,21 @@ The equations that describe the movement of a compressible, inviscid gas (the so-called Euler equations of gas dynamics) are a basic system of conservation laws. In spatial dimension $d$ they read @f[ -\partial_t \mathbf{w} + \nabla \cdot \mathbf{F}(\mathbf{w}) = \mathbf{0}, +\partial_t \mathbf{w} + \nabla \cdot \mathbf{F}(\mathbf{w}) = +\mathbf{G}(\mathbf w), @f] with the solution $\mathbf{w}=(\rho v_1,\ldots,\rho v_d,\rho, E)^{\top}$ consisting of $\rho$ the fluid density, ${\mathbf v}=(v_1,\ldots v_d)^T$ the flow velocity (and thus $\rho\mathbf v$ being the linear momentum density), and $E$ the energy density of the gas. We interpret the equations above as -$\partial_t \mathbf{w}_i + \nabla \cdot \mathbf{F}_i(\mathbf{w}) = 0$, $i=1,\ldots,dim+2$. +$\partial_t \mathbf{w}_i + \nabla \cdot \mathbf{F}_i(\mathbf{w}) = \mathbf +G_i(\mathbf w)$, $i=1,\ldots,dim+2$. For the Euler equations, the flux matrix $\mathbf F$ (or system of flux functions) is defined as (shown here for the case $d=3$) @f{eqnarray*} - \mathbf F + \mathbf F(\mathbf w) = \left( \begin{array}{ccc} @@ -52,12 +54,29 @@ is defined as (shown here for the case $d=3$) \end{array} \right), @f} -such that the entire system of equations reads: +and we will choose as particular right hand side forcing only the effects of +gravity, described by +@f{eqnarray*} + \mathbf G(\mathbf w) + = + \left( + \begin{array}{c} + g_1\rho \\ + g_2\rho \\ + g_3\rho \\ + 0 \\ + \rho \mathbf g \cdot \mathbf v + \end{array} + \right), +@f} +where $\mathbf g=(g_1,g_2,g_3)^T$ denotes the gravity vector. +With this, the entire system of equations reads: @f{eqnarray*} - \partial_t \rho + \sum_{s=1}^d \frac{\partial(\rho v_s)}{\partial x_s} &=& 0, \\ \partial_t (\rho v_i) + \sum_{s=1}^d \frac{\partial(\rho v_i v_s + - \delta_{is} p)}{\partial x_s} &=& 0, \qquad i=1,\dots,d, \\ - \partial_t E + \sum_{s=1}^d \frac{\partial((E+p)v_s)}{\partial x_s} &=& 0. + \delta_{is} p)}{\partial x_s} &=& g_i \rho, \qquad i=1,\dots,d, \\ + \partial_t \rho + \sum_{s=1}^d \frac{\partial(\rho v_s)}{\partial x_s} &=& 0, \\ + \partial_t E + \sum_{s=1}^d \frac{\partial((E+p)v_s)}{\partial x_s} &=& + \rho \mathbf g \cdot \mathbf v. @f} These equations describe, respectively, the conservation of momentum, mass, and energy. @@ -113,9 +132,14 @@ R(\mathbf{W}_{n+1})(\mathbf z) &=& \mathbf{z}\right) - \int_{\Omega} \left(\mathbf{F}(\tilde{\mathbf{w}}), \nabla\mathbf{z}\right) + h^{\eta}(\nabla \mathbf{w} , \nabla \mathbf{z}) +\\ +&& \qquad + \int_{\partial \Omega} \left(\mathbf{H}(\tilde{\mathbf{w}}^+), \mathbf{w}^-(\tilde{\mathbf{w}}^+), \mathbf{n}), \mathbf{z}\right) +- +\int_{\partial \Omega} \left(\mathbf{G}(\mathbf{w}), +\mathbf{z}\right) \\ & = & 0 @f} diff --git a/deal.II/examples/step-33/step-33.cc b/deal.II/examples/step-33/step-33.cc index fb24fe3235..44da0a71fd 100644 --- a/deal.II/examples/step-33/step-33.cc +++ b/deal.II/examples/step-33/step-33.cc @@ -298,8 +298,8 @@ struct EulerEquations // well: template static - void flux_matrix (const InputVector &W, - number (&flux)[n_components][dim]) + void compute_flux_matrix (const InputVector &W, + number (&flux)[n_components][dim]) { // First compute the pressure that // appears in the flux matrix, and @@ -356,8 +356,8 @@ struct EulerEquations Sacado::Fad::DFad iflux[n_components][dim]; Sacado::Fad::DFad oflux[n_components][dim]; - flux_matrix (Wplus, iflux); - flux_matrix (Wminus, oflux); + compute_flux_matrix (Wplus, iflux); + compute_flux_matrix (Wminus, oflux); for (unsigned int di=0; di + static + void compute_forcing_vector (const InputVector &W, + number (&forcing)[n_components]) + { + const double gravity = -1.0; + + for (unsigned int c=0; c::assemble_system () // \left(\mathbf{F}(\tilde{\mathbf{w}}), // \nabla\mathbf{z}_i\right)_K + // h^{\eta}(\nabla \mathbf{w} , \nabla + // \mathbf{z}_i)_K - + // (\mathbf{G}(\tilde{\mathbf w}), // \mathbf{z}_i)_K$ where $\tilde{\mathbf w}$ // is represented by the variable // W_theta, $\mathbf{z}_i$ is @@ -2062,18 +2089,26 @@ assemble_cell_term (const FEValues &fe_v, // Next, in order to compute the cell // contributions, we need to evaluate - // $F(\tilde{\mathbf w})$ at all quadrature + // $F(\tilde{\mathbf w})$ and + // $G(\tilde{\mathbf w})$ at all quadrature // points. To store these, we also need to // allocate a bit of memory. Note that we - // compute the flux matrices in terms of - // autodifferentiation variables, so that - // the Jacobian contributions can later - // easily be computed from it: + // compute the flux matrices and right hand + // sides in terms of autodifferentiation + // variables, so that the Jacobian + // contributions can later easily be + // computed from it: typedef Sacado::Fad::DFad FluxMatrix[EulerEquations::n_components][dim]; FluxMatrix *flux = new FluxMatrix[n_q_points]; + + typedef Sacado::Fad::DFad ForcingVector[EulerEquations::n_components]; + ForcingVector *forcing = new ForcingVector[n_q_points]; for (unsigned int q=0; q::flux_matrix (W_theta[q], flux[q]); + { + EulerEquations::compute_flux_matrix (W_theta[q], flux[q]); + EulerEquations::compute_forcing_vector (W_theta[q], forcing[q]); + } // We now have all of the pieces in place, @@ -2108,9 +2143,13 @@ assemble_cell_term (const FEValues &fe_v, // \mathbf{w}_{\text{component\_i}}}{\partial // x_d} , \frac{\partial // (\mathbf{z}_i)_{\text{component\_i}}}{\partial - // x_d} \right)_K$, where integrals are - // understood to be evaluated through - // summation over quadrature points. + // x_d} \right)_K$ + // $-(\mathbf{G}(\tilde{\mathbf{w}} + // )_{\text{component\_i}}, + // (\mathbf{z}_i)_{\text{component\_i}})_K$, + // where integrals are understood to be + // evaluated through summation over + // quadrature points. // // We initialy sum all contributions of the // residual in the positive sense, so that @@ -2143,28 +2182,16 @@ assemble_cell_term (const FEValues &fe_v, fe_v.shape_grad_component(i, point, component_i)[d] * fe_v.JxW(point); - for (unsigned int d = 0; d < dim; d++) + for (unsigned int d=0; ddiameter(), parameters.diffusion_power) * grad_W[point][component_i][d] * fe_v.shape_grad_component(i, point, component_i)[d] * fe_v.JxW(point); - - // The gravity component only - // enters into the energy equation - // and into the vertical component - // of the velocity. - if (component_i == dim - 1) - F_i += parameters.gravity * - W_theta[point][EulerEquations::density_component] * - fe_v.shape_value_component(i,point, component_i) * - fe_v.JxW(point); - else if (component_i == EulerEquations::energy_component) - F_i += parameters.gravity * - W_theta[point][EulerEquations::density_component] * - W_theta[point][dim-1] * - fe_v.shape_value_component(i,point, component_i) * - fe_v.JxW(point); + + F_i -= forcing[point][component_i] * + fe_v.shape_value_component(i, point, component_i) * + fe_v.JxW(point); } // At the end of the loop, we have to