A bit further into the program.

git-svn-id: https://svn.dealii.org/trunk@16080 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 6ede0eab970d7bac0f6b2afc3abe79e3082dec14..30cce446e6c432af7b9f02a9d999c4b407ca5472 100644 (file)
--- 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 3a241557365edba8ec68ab559a2cbb05caa5704b..a49ae86710c2f46ccdcdaeda8de3c2e783decc6c 100644 (file)
--- 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 <typename number>
     static
-    void flux_matrix(number (&flux)[n_components][dim],
-                    const std::vector<number> &W)
+    void flux_matrix (const std::vector<number> &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; d<dim; ++d)
+         kinetic_energy += W[first_momentum_component+d] *
+                           W[first_momentum_component+d];
+       kinetic_energy *= 1./(2 * W[density_component]);
+       
+       const number pressure = (gas_gamma-1.0)*(W[energy_component] - kinetic_energy);
+
+                                        // Then compute the first
+                                        // <code>dim</code> columns
+                                        // of the matrix that
+                                        // correspond to the momentum
+                                        // terms:
        for (unsigned int d=0; d<dim; ++d)
-         rho_normVsqr += W[d]*W[d];
-                                        // Since W are $\rho v$, we
-                                        // get a $\rho^2$ in the
-                                        // numerator, so dividing a
-                                        // $\rho$ out gives the
-                                        // desired $ \rho |v|^2$.
-       rho_normVsqr /= W[density_component];
-
-       number pressure = (gas_gamma-1.0)*(W[energy_component] - number(0.5)*(rho_normVsqr));
-
-                                        // We compute the momentum terms.  We divide by the
-                                        // density here to get $v_i \rho v_j$
-       for (unsigned int d = 0; d < dim; d++)
          {
-           for (unsigned int d1 = 0; d1 < dim; d1++)
-             flux[d][d1] = W[d]*W[d1]/W[density_component];
+           for (unsigned int e=0; e<dim; ++e)
+             flux[first_momentum_component+d][e] = W[first_momentum_component+d] *
+                                                   W[first_momentum_component+e] /
+                                                   W[density_component];
          
-                                            // The pressure contribution, along the diagonal:
-           flux[d][d] += pressure;
-                                            // Advection/conservation of density:
-           flux[density_component][d] = W[d]; 
-                                            // And, lastly, conservation of energy.
-           flux[energy_component][d] = W[d]/W[density_component]*
-                                       (W[energy_component] + pressure); // energy
+           flux[first_momentum_component+d][d] += pressure;
          }
+       
+                                        // Then the terms for the
+                                        // density (i.e. mass
+                                        // conservation):
+       for (unsigned int d=0; d<dim; ++d)
+         flux[density_component][d] = W[first_momentum_component+d]; 
+
+                                        // And, lastly, conservation
+                                        // of energy:
+       for (unsigned int d=0; d<dim; ++d)
+         flux[energy_component][d] = W[first_momentum_component+d] /
+                                     W[density_component] *
+                                     (W[energy_component] + pressure);
       }
 
 
@@ -214,32 +240,25 @@ struct EulerEquations
                                     // $\alpha$.
     template <typename number>
     static
-    void LFNumFlux(std::vector<std::vector<Sacado::Fad::DFad<double> > > &nflux,
-                  const std::vector<Point<dim> > &points, 
-                  const std::vector<Point<dim> > &normals,
-                  const std::vector<std::vector<number> > &Wplus,
-                  const std::vector<std::vector<number> > &Wminus,
-                  double alpha)
+    void numerical_normal_flux(const Point<dim> &normal,
+                              const std::vector<number> &Wplus,
+                              const std::vector<number> &Wminus,
+                              const double alpha,
+                              Sacado::Fad::DFad<double> (&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<double> iflux[n_components][dim];
-           Sacado::Fad::DFad<double> oflux[n_components][dim];
+       Sacado::Fad::DFad<double> iflux[n_components][dim];
+       Sacado::Fad::DFad<double> 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<n_components; ++di)
-             {
-               nflux[q][di] = 0;
-               for (unsigned int d=0; d<dim; ++d) 
-                 nflux[q][di] += 0.5*(iflux[di][d] + oflux[di][d])*normals[q](d);
+       for (unsigned int di=0; di<n_components; ++di)
+         {
+           normal_flux[di] = 0;
+           for (unsigned int d=0; d<dim; ++d) 
+             normal_flux[di] += 0.5*(iflux[di][d] + oflux[di][d]) * normal(d);
              
-               nflux[q][di] += 0.5*alpha*(Wplus[q][di] - Wminus[q][di]);
-             }
+           normal_flux[di] += 0.5*alpha*(Wplus[di] - Wminus[di]);
          }
     }
 };
@@ -662,7 +681,7 @@ void ConsLaw<dim>::assemble_cell_term(
   FluxMatrix *flux = new FluxMatrix[n_q_points];
   
   for (unsigned int q=0; q < n_q_points; ++q)
-    EulerEquations<dim>::flux_matrix(flux[q], Wcn[q]);
+    EulerEquations<dim>::flux_matrix(Wcn[q], flux[q]);
   
 
                                   // We now have all of the function values/grads/fluxes,
@@ -900,7 +919,9 @@ void ConsLaw<dim>::assemble_face_term(
    
                                   // Determine the Lax-Friedrich's stability parameter,
                                   // and evaluate the numerical flux function at the quadrature points
-  std::vector<std::vector<Sacado::Fad::DFad<double> > > nflux(n_q_points, std::vector<Sacado::Fad::DFad<double> >(EulerEquations<dim>::n_components, 0));
+  typedef Sacado::Fad::DFad<double> NormalFlux[EulerEquations<dim>::n_components];
+  NormalFlux *normal_fluxes = new NormalFlux[n_q_points];
+
   double alpha = 1;
 
   switch(flux_params.LF_stab) {
@@ -912,43 +933,46 @@ void ConsLaw<dim>::assemble_face_term(
          break;
   }
 
-  EulerEquations<dim>::LFNumFlux(nflux, fe_v.get_quadrature_points(), normals, Wplus, Wminus,
-                            alpha);
+  for (unsigned int q=0; q<n_q_points; ++q)
+    EulerEquations<dim>::numerical_normal_flux(normals[q], Wplus[q], Wminus[q], alpha,
+                                              normal_fluxes[q]);
 
                                   // Now assemble the face term
-  for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i) {
-    if (!fe_v.get_fe().has_support_on_face(i, face_no)) continue;
-    F_i = 0;
-    for (unsigned int point=0; point<n_q_points; ++point)
-      {
-       const unsigned int
-         component_i = fe_v.get_fe().system_to_component_index(i).first;
-
-       F_i += nflux[point][component_i]*fe_v.shape_value_component(i, point, component_i)*fe_v.JxW(point);
-
-      } 
-
-                                    // Retrieve a pointer to the jacobian.
-    double *values = &(F_i.fastAccessDx(0));
+  for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
+    {
+      if (!fe_v.get_fe().has_support_on_face(i, face_no))
+       continue;
+      
+      F_i = 0;
+      for (unsigned int point=0; point<n_q_points; ++point)
+       {
+         const unsigned int
+           component_i = fe_v.get_fe().system_to_component_index(i).first;
+         
+         F_i += normal_fluxes[point][component_i] *
+                fe_v.shape_value_component(i, point, component_i) *
+                fe_v.JxW(point);
+       } 
 
-                                    // Honestly, I forget why this can happen, but 
-                                    // for some reason it can!!
-    if (!values) continue;
+                                      // Retrieve a pointer to the jacobian.
+      double *values = &(F_i.fastAccessDx(0));
+      Assert (values != 0, ExcInternalError());
 
-                                    // 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[0], reinterpret_cast<int*>(&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<int*>(&dofs_neighbor[0]));
-    }
+                                 dofs_per_cell, &values[0], reinterpret_cast<int*>(&dofs[0]));
+      if (boundary < 0) {
+       Matrix->SumIntoGlobalValues(dofs[i],
+                                   dofs_per_cell, &values[dofs_per_cell], reinterpret_cast<int*>(&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