]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add the forgotten forcing term to the documentation in various places. Make it generi...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 20 May 2008 16:38:57 +0000 (16:38 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 20 May 2008 16:38:57 +0000 (16:38 +0000)
git-svn-id: https://svn.dealii.org/trunk@16137 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-33/doc/intro.dox
deal.II/examples/step-33/step-33.cc

index 79fdb7a1ee835ef5421932bd9cf729acfa0bc966..ef696f4c252bb644c1c079fbd7e73f16b039b908 100644 (file)
@@ -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}
index fb24fe3235852a98f16aea326fe2122543c3b0f7..44da0a71fdcb742686f0f18881eed41faae10694 100644 (file)
@@ -298,8 +298,8 @@ struct EulerEquations
                                     // well:
     template <typename InputVector, typename number>
     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<double> iflux[n_components][dim];
        Sacado::Fad::DFad<double> 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<n_components; ++di)
          {
@@ -369,6 +369,31 @@ struct EulerEquations
          }
       }
 
+
+    template <typename InputVector, typename number>
+    static
+    void compute_forcing_vector (const InputVector &W,
+                                number (&forcing)[n_components])
+      {
+       const double gravity = -1.0;
+
+       for (unsigned int c=0; c<n_components; ++c)
+         switch (c)
+           {
+             case first_momentum_component+dim-1:
+                   forcing[c] = gravity * W[density_component];
+                   break;
+             case energy_component:
+                   forcing[c] = gravity *
+                                W[density_component] *
+                                W[first_momentum_component+dim-1];
+                   break;
+             default:
+                   forcing[c] = 0;
+           }
+      }
+    
+    
     
                                     // Finally, we declare a class that
                                     // implements a postprocessing of data
@@ -1875,6 +1900,8 @@ void ConservationLaw<dim>::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
                                  // <code>W_theta</code>, $\mathbf{z}_i$ is
@@ -2062,18 +2089,26 @@ assemble_cell_term (const FEValues<dim>             &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<double> FluxMatrix[EulerEquations<dim>::n_components][dim];
   FluxMatrix *flux = new FluxMatrix[n_q_points];
+
+  typedef Sacado::Fad::DFad<double> ForcingVector[EulerEquations<dim>::n_components];
+  ForcingVector *forcing = new ForcingVector[n_q_points];
   
   for (unsigned int q=0; q<n_q_points; ++q)
-    EulerEquations<dim>::flux_matrix (W_theta[q], flux[q]);
+    {
+      EulerEquations<dim>::compute_flux_matrix (W_theta[q], flux[q]);
+      EulerEquations<dim>::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<dim>             &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<dim>             &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; d<dim; d++)
            F_i += 1.0*std::pow(fe_v.get_cell()->diameter(),
                                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<dim>::density_component] *
-                  fe_v.shape_value_component(i,point, component_i) *
-                  fe_v.JxW(point);
-         else if (component_i == EulerEquations<dim>::energy_component)
-           F_i += parameters.gravity *
-                  W_theta[point][EulerEquations<dim>::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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.