From 469ab6b7bc4dfa5940261a0ff24a8819a32bfad2 Mon Sep 17 00:00:00 2001 From: bangerth Date: Mon, 12 May 2008 20:26:00 +0000 Subject: [PATCH] A bit further into the program. git-svn-id: https://svn.dealii.org/trunk@16080 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-33/doc/intro.dox | 20 ++- deal.II/examples/step-33/step-33.cc | 212 ++++++++++++++----------- 2 files changed, 134 insertions(+), 98 deletions(-) diff --git a/deal.II/examples/step-33/doc/intro.dox b/deal.II/examples/step-33/doc/intro.dox index 6ede0eab97..30cce446e6 100644 --- a/deal.II/examples/step-33/doc/intro.dox +++ b/deal.II/examples/step-33/doc/intro.dox @@ -30,20 +30,32 @@ 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}, @f] -with the solution $\mathbf{w}=(\rho,\rho v_1,\ldots,\rho v_d, +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. The flux matrix $\mathbf F$ (or system of flux functions) -is defined such that the entire system of equations are +is defined as (shown here for the case $d=3$) +@f{eqnarray*} + \mathbf F + = + \left( + \begin{array}{ccc} + \rho v_1^2+p & \rho v_2v_1 & \rho v_3v_1 & \rho v_1 & (E+p)v_1 \\ + \rho v_1v_2 & \rho v_2^2+p & \rho v_3v_2 & \rho v_2 & (E+p)v_2 \\ + \rho v_1v_3 & \rho v_2v_3 & \rho v_3^2+p & \rho v_3 & (E+p)v_3 + \end{array} + \right), +@f} +such that 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. @f} -These equations describe, respectively, the conservation of mass, -momentum, and energy. +These equations describe, respectively, the conservation of momentum, +mass, and energy. The system is closed by a relation that defines the pressure: $p = (\gamma -1)(E-\frac{1}{2} \rho |\mathbf v|^2)$. For the constituents of air (mainly nitrogen and oxygen) and other diatomic gases, the ratio of diff --git a/deal.II/examples/step-33/step-33.cc b/deal.II/examples/step-33/step-33.cc index 3a24155736..a49ae86710 100644 --- a/deal.II/examples/step-33/step-33.cc +++ b/deal.II/examples/step-33/step-33.cc @@ -157,54 +157,80 @@ struct EulerEquations // and $O_2$. static const double gas_gamma; - // We define the flux function $F(W)$ as one large - // matrix. Each row of this matrix - // represents a scalar conservation law for - // the component in that row. We templatize - // the numerical type of the flux function so + // We define the flux function + // $F(W)$ as one large matrix. + // Each row of this matrix + // represents a scalar + // conservation law for the + // component in that row. The + // exact form of this matrix is + // given in the introduction. + // + // We templatize the numerical + // type of the flux function so // that we may use the automatic - // differentiation type here. The flux - // functions are defined in terms of the - // conserved variables $\rho w_0, \dots, \rho - // w_{d-1}, \rho, E$, so they do not look - // exactly like the Euler equations one is - // used to seeing. We evaluate the flux at a - // single quadrature point. + // differentiation type here. + // The flux functions are defined + // in terms of the conserved + // variables $\rho w_0, \dots, + // \rho w_{d-1}, \rho, E$, so + // they do not look exactly like + // the Euler equations one is + // used to seeing. We evaluate + // the flux at a single + // quadrature point. template static - void flux_matrix(number (&flux)[n_components][dim], - const std::vector &W) + void flux_matrix (const std::vector &W, + number (&flux)[n_components][dim]) { - - // Pressure is a dependent variable: $p = - // (\gas_gamma - 1)(E-\frac{1}{2} \rho |v|^2)$. - number rho_normVsqr; + // First compute the pressure + // that appears in the flux + // matrix, based on the + // energy density and the + // kinetic energy $\frac 12 + // \rho |\mathbf v|^2 = + // \frac{|\rho \mathbf + // v|^2}{2\rho}$ (note that + // the independent variables + // contain the momentum + // components $\rho v_i$, not + // the velocities $v_i$): + number kinetic_energy = 0; + for (unsigned int d=0; ddim columns + // of the matrix that + // correspond to the momentum + // terms: for (unsigned int d=0; d static - void LFNumFlux(std::vector > > &nflux, - const std::vector > &points, - const std::vector > &normals, - const std::vector > &Wplus, - const std::vector > &Wminus, - double alpha) + void numerical_normal_flux(const Point &normal, + const std::vector &Wplus, + const std::vector &Wminus, + const double alpha, + Sacado::Fad::DFad (&normal_flux)[n_components]) { - const unsigned int n_q_points = points.size(); - - // We evaluate the flux at each of the quadrature points. - for (unsigned int q = 0; q < n_q_points; q++) - { - Sacado::Fad::DFad iflux[n_components][dim]; - Sacado::Fad::DFad oflux[n_components][dim]; + Sacado::Fad::DFad iflux[n_components][dim]; + Sacado::Fad::DFad oflux[n_components][dim]; - flux_matrix(iflux, Wplus[q]); - flux_matrix(oflux, Wminus[q]); + flux_matrix(Wplus, iflux); + flux_matrix(Wminus, oflux); - for (unsigned int di=0; di::assemble_cell_term( FluxMatrix *flux = new FluxMatrix[n_q_points]; for (unsigned int q=0; q < n_q_points; ++q) - EulerEquations::flux_matrix(flux[q], Wcn[q]); + EulerEquations::flux_matrix(Wcn[q], flux[q]); // We now have all of the function values/grads/fluxes, @@ -900,7 +919,9 @@ void ConsLaw::assemble_face_term( // Determine the Lax-Friedrich's stability parameter, // and evaluate the numerical flux function at the quadrature points - std::vector > > nflux(n_q_points, std::vector >(EulerEquations::n_components, 0)); + typedef Sacado::Fad::DFad NormalFlux[EulerEquations::n_components]; + NormalFlux *normal_fluxes = new NormalFlux[n_q_points]; + double alpha = 1; switch(flux_params.LF_stab) { @@ -912,43 +933,46 @@ void ConsLaw::assemble_face_term( break; } - EulerEquations::LFNumFlux(nflux, fe_v.get_quadrature_points(), normals, Wplus, Wminus, - alpha); + for (unsigned int q=0; q::numerical_normal_flux(normals[q], Wplus[q], Wminus[q], alpha, + normal_fluxes[q]); // Now assemble the face term - for (unsigned int i=0; iSumIntoGlobalValues(dofs[i], - dofs_per_cell, &values[0], reinterpret_cast(&dofs[0])); - if (boundary < 0) { + // Update the matrix. Depending on whether there + // is/isn't a neighboring cell, we add more/less + // entries. Matrix->SumIntoGlobalValues(dofs[i], - dofs_per_cell, &values[dofs_per_cell], reinterpret_cast(&dofs_neighbor[0])); - } + dofs_per_cell, &values[0], reinterpret_cast(&dofs[0])); + if (boundary < 0) { + Matrix->SumIntoGlobalValues(dofs[i], + dofs_per_cell, &values[dofs_per_cell], reinterpret_cast(&dofs_neighbor[0])); + } - // And add into the residual - right_hand_side(dofs[i]) -= F_i.val(); - } + // And add into the residual + right_hand_side(dofs[i]) -= F_i.val(); + } + delete[] normal_fluxes; } // @sect4{Assembling the whole system} // Now we put all of the assembly pieces together -- 2.39.5