]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement the postprocessor. Will hook it in tomorrow.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 16 May 2008 03:41:39 +0000 (03:41 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 16 May 2008 03:41:39 +0000 (03:41 +0000)
git-svn-id: https://svn.dealii.org/trunk@16108 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-33/step-33.cc

index 4a2edfc2e8c9791e01ed641c1790bd852c4ac15f..99be279ecdb86d66a4c5869bbd358202f2492ebd 100644 (file)
@@ -158,6 +158,67 @@ struct EulerEquations
                                   // and $O_2$.
     static const double gas_gamma;
 
+    
+                                    // In the following, we will need to
+                                    // compute the kinetic energy and the
+                                    // pressure from a vector of conserved
+                                    // variables. This we can do 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$).
+                                    //
+                                    // There is one slight problem: We will
+                                    // need to call the following functions
+                                    // with input arguments of type
+                                    // <code>std::vector@<number@></code> and
+                                    // <code>Vector@<number@></code>. The
+                                    // problem is that the former has an
+                                    // access operator
+                                    // <code>operator[]</code> whereas the
+                                    // latter, for historical reasons, has
+                                    // <code>operator()</code>. We wouldn't
+                                    // be able to write the function in a
+                                    // generic way if we were to use one or
+                                    // the other of these. Fortunately, we
+                                    // can use the following trick: instead
+                                    // of writing <code>v[i]</code> or
+                                    // <code>v(i)</code>, we can use
+                                    // <code>*(v.begin() + i)</code>, i.e. we
+                                    // generate an iterator that points to
+                                    // the <code>i</code>th element, and then
+                                    // dereference it. This works for both
+                                    // kinds of vectors -- not the prettiest
+                                    // solution, but one that works.
+    template <typename number, typename InputVector>
+    static
+    number
+    compute_kinetic_energy (const InputVector &W)
+      {
+       number kinetic_energy = 0;
+       for (unsigned int d=0; d<dim; ++d)
+         kinetic_energy += *(W.begin()+first_momentum_component+d) *
+                           *(W.begin()+first_momentum_component+d);
+       kinetic_energy *= 1./(2 * *(W.begin() + density_component));
+
+       return kinetic_energy;
+      }
+
+
+    template <typename number, typename InputVector>
+    static
+    number
+    compute_pressure (const InputVector &W)
+      {
+       return ((gas_gamma-1.0) *
+               (*(W.begin() + energy_component) -
+                compute_kinetic_energy<number>(W)));
+      }        
+    
+
+    
                                     // We define the flux function
                                     // $F(W)$ as one large matrix.
                                     // Each row of this matrix
@@ -197,37 +258,21 @@ struct EulerEquations
     void flux_matrix (const std::vector<number> &W,
                      number (&flux)[n_components][dim])
       {
-                                        // 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]);
+                                        // First compute the pressure that
+                                        // appears in the flux matrix, and
+                                        // then compute the first
+                                        // <code>dim</code> columns of the
+                                        // matrix that correspond to the
+                                        // momentum terms:
+       const number pressure = compute_pressure<number> (W);
        
-       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)
          {
            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];
+             flux[first_momentum_component+d][e]
+               = W[first_momentum_component+d] *
+               W[first_momentum_component+e] /
+               W[density_component];
          
            flux[first_momentum_component+d][d] += pressure;
          }
@@ -280,12 +325,276 @@ struct EulerEquations
            normal_flux[di] += 0.5*alpha*(Wplus[di] - Wminus[di]);
          }
     }
+
+    
+                                    // Finally, we declare a class that
+                                    // implements a postprocessing of data
+                                    // components. The problem this class
+                                    // solves is that the variables in the
+                                    // formulation of the Euler equations we
+                                    // use are in conservative rather than
+                                    // physical form: they are momentum
+                                    // densities $\mathbf m=\rho\mathbf v$,
+                                    // density $\rho$, and energy density
+                                    // $E$. What we would like to also put
+                                    // into our output file are velocities
+                                    // $\mathbf v=\frac{\mathbf m}{\rho}$ and
+                                    // pressure $p=(\gamma-1)(E-\frac{1}{2}
+                                    // \rho |\mathbf v|^2)$.
+                                    //
+                                    // In addition, we would like to add the
+                                    // possibility to generate schlieren
+                                    // plots. Schlieren plots are a way to
+                                    // visualize shocks and other sharp
+                                    // interfaces. The word "schlieren" a
+                                    // German word that may be translated as
+                                    // "striae" -- it may be simpler to
+                                    // explain it by an example, however:
+                                    // schlieren is what you see when you,
+                                    // for example, pour highly concentrated
+                                    // alcohol, or a transparent saline
+                                    // solution into water; the two have the
+                                    // same color, but they have different
+                                    // refractive indices and so before they
+                                    // are fully mixed light goes through the
+                                    // mixture along bent rays that lead to
+                                    // brightness variations if you look at
+                                    // it. That's "schlieren". A similar
+                                    // effect happens in compressible flow
+                                    // due because the refractive index
+                                    // depends on the pressure (and therefore
+                                    // the density) of the gas.
+                                    //
+                                    // The origin of the word refers to
+                                    // two-dimensional projections of a
+                                    // three-dimensional volume (we see a 2d
+                                    // picture of the 3d fluid). In
+                                    // computational fluid dynamics, we can
+                                    // get an idea of this effect by
+                                    // considering what causes it: density
+                                    // variations. Schlieren plots are
+                                    // therefore produced by plotting
+                                    // $s=|\nabla \rho|^2$; obviously, $s$ is
+                                    // large in shocks and at other highly
+                                    // dynamic places. If so desired by the
+                                    // user (by specifying this in the input
+                                    // file), we would like to generate these
+                                    // schlieren plots in addition to the
+                                    // other derived quantities listed above.
+                                    //
+                                    // The implementation of the algorithms
+                                    // to compute derived quantities from the
+                                    // ones that solve our problem, and to
+                                    // output them into data file, rests on
+                                    // the DataPostprocessor class. It has
+                                    // extensive documentation, and other
+                                    // uses of the class can also be found in
+                                    // step-29. We therefore refrain from
+                                    // extensive comments.
+    class Postprocessor : public DataPostprocessor<dim>
+    {
+      public:
+       Postprocessor (const bool do_schlieren_plot);
+       
+       virtual
+       void
+       compute_derived_quantities_vector (const std::vector<Vector<double> >              &uh,
+                                          const std::vector<std::vector<Tensor<1,dim> > > &duh,
+                                          const std::vector<std::vector<Tensor<2,dim> > > &dduh,
+                                          const std::vector<Point<dim> >                  &normals,
+                                          std::vector<Vector<double> >                    &computed_quantities) const;
+
+       virtual std::vector<std::string> get_names () const;
+    
+       virtual
+       std::vector<DataComponentInterpretation::DataComponentInterpretation>
+       get_data_component_interpretation () const;
+    
+       virtual UpdateFlags get_needed_update_flags () const;
+
+       virtual unsigned int n_output_variables() const;
+
+      private:
+       const bool do_schlieren_plot;
+    };
 };
     
 
 template <int dim>
 const double EulerEquations<dim>::gas_gamma = 1.4;
 
+
+
+template <int dim>
+EulerEquations<dim>::Postprocessor::
+Postprocessor (const bool do_schlieren_plot)
+               :
+               do_schlieren_plot (do_schlieren_plot)
+{}
+
+
+                                // This is the only function worth commenting
+                                // on. When generating graphical output, the
+                                // DataOut and related classes will call this
+                                // function on each cell, with values,
+                                // gradients, hessians, and normal vectors
+                                // (in case we're working on faces) at each
+                                // quadrature point. Note that the data at
+                                // each quadrature point is itself
+                                // vector-valued, namely the conserved
+                                // variables. What we're going to do here is
+                                // to compute the quantities we're interested
+                                // in at each quadrature point. Note that for
+                                // this we can ignore the hessians ("dduh")
+                                // and normal vectors; to avoid compiler
+                                // warnings about unused variables, we
+                                // comment out their names.
+template <int dim>
+void
+EulerEquations<dim>::Postprocessor::
+compute_derived_quantities_vector (const std::vector<Vector<double> >              &uh,
+                                  const std::vector<std::vector<Tensor<1,dim> > > &duh,
+                                  const std::vector<std::vector<Tensor<2,dim> > > &/*dduh*/,
+                                  const std::vector<Point<dim> >                  &/*normals*/,
+                                  std::vector<Vector<double> >                    &computed_quantities) const
+{
+                                  // At the beginning of the function, let us
+                                  // make sure that all variables have the
+                                  // correct sizes, so that we can access
+                                  // individual vector elements without
+                                  // having to wonder whether we might read
+                                  // or write invalid elements; we also check
+                                  // that the <code>duh</code> vector only
+                                  // contains data if we really need it (the
+                                  // system knows about this because we say
+                                  // so in the
+                                  // <code>get_needed_update_flags()</code>
+                                  // function below). For the inner vectors,
+                                  // we check that at least the first element
+                                  // of the outer vector has the correct
+                                  // inner size:
+  const unsigned int n_quadrature_points = uh.size();
+
+  if (do_schlieren_plot == true)
+    Assert (duh.size() == n_quadrature_points,
+           ExcInternalError())
+  else
+    Assert (duh.size() == 0,
+           ExcInternalError());
+
+  Assert (computed_quantities.size() == n_quadrature_points,
+         ExcInternalError());
+
+  Assert (uh[0].size() == n_components,
+         ExcInternalError());
+
+  if (do_schlieren_plot == true)
+    Assert (computed_quantities[0].size() == dim+2,
+           ExcInternalError())
+  else
+    Assert (computed_quantities[0].size() == dim+1,
+           ExcInternalError());
+
+                                  // Then loop over all quadrature points and
+                                  // do our work there. The code should be
+                                  // pretty self-explanatory. The order of
+                                  // output variables is first
+                                  // <code>dim</code> velocities, then the
+                                  // pressure, and if so desired the
+                                  // schlieren plot. Note that we try to be
+                                  // generic about the order of variables in
+                                  // the input vector, using the
+                                  // <code>first_momentum_component</code>
+                                  // and <code>density_component</code>
+                                  // information:
+  for (unsigned int q=0; q<n_quadrature_points; ++q)
+    {
+      const double density = uh[q](density_component);
+      const double energy  = uh[q](energy_component);
+
+      for (unsigned int d=0; d<dim; ++d)
+       computed_quantities[q](d)
+         = uh[q](first_momentum_component+d) / density;
+
+      computed_quantities[q](dim) = compute_pressure<double> (uh[q]);
+
+      if (do_schlieren_plot == true)
+       computed_quantities[q](dim+1) = duh[q][density_component] *
+                                       duh[q][density_component];
+    }
+}
+
+
+template <int dim>
+std::vector<std::string>
+EulerEquations<dim>::Postprocessor::
+get_names () const
+{
+  std::vector<std::string> names;
+  for (unsigned int d=0; d<dim; ++d)
+    names.push_back ("velocity");
+  names.push_back ("pressure");
+
+  if (do_schlieren_plot == true)
+    names.push_back ("schlieren_plot");
+
+  return names;
+}
+
+
+template <int dim>
+std::vector<DataComponentInterpretation::DataComponentInterpretation>
+EulerEquations<dim>::Postprocessor::
+get_data_component_interpretation () const
+{
+  std::vector<DataComponentInterpretation::DataComponentInterpretation>
+    interpretation;
+
+  for (unsigned int d=0; d<dim; ++d)
+    interpretation.push_back (DataComponentInterpretation::
+                             component_is_part_of_vector);
+
+  interpretation.push_back (DataComponentInterpretation::
+                           component_is_scalar);
+
+  if (do_schlieren_plot == true)
+    interpretation.push_back (DataComponentInterpretation::
+                             component_is_scalar);
+
+  return interpretation;
+}
+
+
+
+template <int dim>
+UpdateFlags
+EulerEquations<dim>::Postprocessor::
+get_needed_update_flags () const
+{
+  if (do_schlieren_plot == true)  
+    return update_values | update_gradients;
+  else
+    return update_values;
+}
+
+
+
+template <int dim>
+unsigned int
+EulerEquations<dim>::Postprocessor::
+n_output_variables () const
+{
+  if (do_schlieren_plot == true)
+    return dim+2;
+  else
+    return dim+1;
+}
+
+
+template class EulerEquations<2>::Postprocessor;
+
+
+
                                 // @sect3{Run time parameter handling}
 
                                 // Our next job is to define a few
@@ -588,8 +897,10 @@ namespace Parameters
       const std::string stab = prm.get("stab");
       if (stab == "constant") 
        stabilization_kind = constant;
-      else if (stab == "mesh ") 
+      else if (stab == "mesh") 
        stabilization_kind = mesh_dependent;
+      else
+       AssertThrow (false, ExcNotImplemented());
   
       stabilization_value = prm.get_double("stab value");
     }
@@ -1633,10 +1944,10 @@ void ConsLaw<dim>::solve (Vector<double> &dsolution, int &niter, double &lin_res
   }
 }
 
-                                 // @sect3{Postprocessing and Output}
-                                 // Recover the physical variables from the conservative
-                                 // variables so that output will be (perhaps) more
-                                 // meaningfull.
+                                 // @sect3{Postprocessing and Output} Recover
+                                 // the physical variables from the
+                                 // conservative variables so that output will
+                                 // be (perhaps) more meaningfull.
 template <int dim>
 void ConsLaw<dim>::postprocess() {
   const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;

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.