This program was written for fun by David Neckels (NCAR) while working at Sandia (on the Wyoming Express bus to and from Corrales each day). The main purpose was to better understand Euler flow. The code solves the basic Euler equations of gas dynamics, by using a fully implicit Newton iteration (inspired by Sandia's Aria code). The code may be configured by an input deck to run different simulations on different meshes, with differing boundary conditions.
The original code and documentation was later slightly modified by Wolfgang Bangerth to make it more modular and allow replacing the parts that are specific to the Euler equations by other hyperbolic conservation laws without too much trouble.
Note:The program uses the Trilinos linear solvers (which are part of the Aztec/Amesos package of Trilinos) and an automatic differentiation package, Sacado, also part of Trilinos. deal.II must be configured to use this package. Refer to the ReadMe file for instructions how to do this.
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
they read
with the solution
consisting of
the fluid density,
the flow velocity (and thus
being the linear momentum density), and
the energy density of the gas. We interpret the equations above as
,
.
For the Euler equations, the flux matrix
(or system of flux functions) is defined as (shown here for the case
)
and we will choose as particular right hand side forcing only the effects of gravity, described by
where
denotes the gravity vector. With this, the entire system of equations reads:
These equations describe, respectively, the conservation of momentum, mass, and energy. The system is closed by a relation that defines the pressure:
. For the constituents of air (mainly nitrogen and oxygen) and other diatomic gases, the ratio of specific heats is
.
This problem obviously falls into the class of vector-valued problem. A general overview of how to deal with these problems in deal.II can be found in the Handling vector valued problems module.
Discretization happens in the usual way, taking into account that this is a hyperbolic problem in the same style as the simple one discussed in step-12: We choose a finite element space
, and integrate our conservation law against our (vector-valued) test function
. We then integrate by parts and approximate the boundary flux with a numerical flux
,
where a superscript
denotes the interior trace of a function, and
represents the outer trace. The diffusion term
is introduced strictly for stability, where
is the mesh size and
is a parameter prescribing how much diffusion to add.
On the boundary, we have to say what the outer trace
is. Depending on the boundary condition, we prescribe either of the following:
is prescribed to be the desired value.
except that the energy variable is modified to support a prescribed pressure
, i.e.
so that
and
. More information on these issues can be found, for example, in Ralf Hartmann's PhD thesis ("Adaptive Finite Element Methods for the Compressible Euler Equations", PhD thesis, University of Heidelberg, 2002).
We use a time stepping scheme to substitute the time derivative in the above equations. At each time step, our full discretization is thus that the residual applied to any test function
equals zero:
where
for
and
. Choosing
results in the explicit (forward) Euler scheme,
in the stable implicit (backward) Euler scheme, and
in the Crank-Nicolson scheme.
In the implementation below, we choose the Lax-Friedrichs flux for the function
, i.e.
, where
is either a fixed number specified in the input file, or where
is a mesh dependent value. In the latter case, it is chosen as
with
the diameter of the face to which the flux is applied, and
the current time step.
With these choices, equating the residual to zero results in a nonlinear system of equations which we solve the nonlinear system by a Newton iteration, i.e. by iterating
until
(the residual) is sufficiently small. By testing with the nodal basis of a finite element space instead of all
, we arrive at a linear system for
:
This linear system is, in general, neither symmetric nor has any particular definiteness properties. We will either use a direct solver or Trilinos' GMRES implementation to solve it. As will become apparent from the results shown below, this fully implicit iteration converges very rapidly (typically in 3 steps) and with the quadratic convergence order expected from a Newton method.
Since computing the Jacobian matrix
is a terrible beast, we use an automatic differentiation package, Sacado, to do this. Sacado is a package within the Trilinos framework and offers a C++ template class Sacado::Fad::DFad (Fad standing for "floating point automatic differentiation") that supports basic arithmetic operators and functions such as sqrt, sin, cos, pow, etc. In order to use this feature, one declares a collection of variables of this type and then denotes some of this collection as degrees of freedom, the rest of the variables being functions of the independent variables. These variables are used in an algorithm, and as the variables are used, their sensitivities with respect to the degrees of freedom are continuously updated.
One can imagine that for the full Jacobian, this could be prohibitively expensive: the number of independent variables are the
, the dependent variables the elements of the vector
. Both of these vectors can easily have tens of thousands of elements or more. However, it is important to note that not all elements of
depend on all elements of
: in fact, an entry in
only depends on an element of
if the two corresponding shape functions overlap and couple in the weak form. This means that it is enough if we do not use the Sacado type for the entire matrix computation, but only element by element.
The author has used this approach side by side with a hand coded Jacobian for the incompressible Navier-Stokes problem and found the Sacado approach to be just as fast as using a hand coded Jacobian, but infinitely simpler and less error prone: Since using the auto-differentiation requires only that one code the residual
, ensuring code correctness and maintaining code becomes tremendously more simple -- the Jacobian matrix
is computed by essentially the same code that also computes the residual
.
All this said, here's a very simple example showing how Sacado can be used:
#include <Sacado.hpp> #include <iostream> typedef Sacado::Fad::DFad<double> fad_double; main() { fad_double a,b,c; a = 1; b = 2; a.diff(0,2); // Set a to be dof 0, in a 2-dof system. b.diff(1,2); // Set b to be dof 1, in a 2-dof system. c = 2*a+cos(a*b); double *derivs = &c.fastAccessDx(0); // Access derivatives std::cout << "dc/da = " << derivs[0] << ", dc/db=" << derivs[1] << std::endl; }
The output are the derivatives
of
at
.
It should be noted that Sacado provides more auto-differentation capabilities than the small subset used in this program. However, understanding the example above is enough to understand the use of Sacado in this Euler flow program.
The program uses either the Aztec iterative solvers, or the Amesos sparse direct solver, both provided by the Trilinos package. This package is inherently designed to be used in a parallel program, however, it may be used in serial just as easily, as is done here. The Epetra package is the basic vector/matrix library upon which the solvers are built. This very powerful package can be used to describe the parallel distribution of a vector, and to define sparse matrices that operate on these vectors. Please view the commented code for more details on how these solvers are used within the example.
The example uses an ad-hoc refinement indicator that shows some usefulness in shock-type problems, and in the downhill flow example included. We refine according to the squared gradient of the density. Hanging nodes are handled by computing the numerical flux across cells that are of differing refinement levels, rather than using the ConstraintMatrix class as in all other tutorial programs so far. In this way, the example combines the continuous and DG methodologies. It also simplifies the generation of the Jacobian because we do not have to track constrained degrees of freedom through the automatic differentiation used to compute it.
Further, we enforce a maximum number of refinement levels to keep refinement under check. It is the author's experience that for adaptivity for a time dependent problem, refinement can easily lead the simulation to a screeching halt, because of time step restrictions if the mesh becomes too fine in any part of the domain, if care is not taken. The amount of refinement is limited in the example by letting the user specify the maximum level of refinement that will be present anywhere in the mesh. In this way, refinement tends not to slow the simulation to a halt. This, of course, is purely a heuristic strategy, and if the author's advisor heard about it, the author would likely be exiled forever from the finite element error estimation community.
We use an input file deck to drive the simulation. In this way, we can alter the boundary conditions and other important properties of the simulation without having to recompile. For more information on the format, look at the results section, where we describe an example input file in more detail.
In previous example programs, we have usually hard-coded the initial and boundary conditions. In this program, we instead use the expression parser class FunctionParser so that we can specify a generic expression in the input file and have it parsed at run time — this way, we can change initial conditions without the need to recompile the program. Consequently, no classes named InitialConditions or BoundaryConditions will be declared in the program below.
The implementation of this program is split into three essential parts:
EulerEquations class that encapsulates everything that completely describes the specifics of the Euler equations. This includes the flux matrix
, the numerical flux
, the right hand side
, boundary conditions, refinement indicators, postprocessing the output, and similar things that require knowledge of the meaning of the individual components of the solution vectors and the equations.
ConservationLaw class that deals with time stepping, outer nonlinear and inner linear solves, assembling the linear systems, and the top-level logic that drives all this.
The reason for this approach is that it separates the various concerns in a program: the ConservationLaw is written in such a way that it would be relatively straightforward to adapt it to a different set of equations: One would simply re-implement the members of the EulerEquations class for some other hyperbolic equation, or augment the existing equations by additional ones (for example by advecting additional variables, or by adding chemistry, etc). Such modifications, however, would not affect the time stepping, or the nonlinear solvers if correctly done, and consequently nothing in the ConservationLaw would have to be modified.
Similarly, if we wanted to improve on the linear or nonlinear solvers, or on the time stepping scheme (as hinted at at the end of the results section), then this would not require changes in the EulerEquations at all.
First a standard set of deal.II includes. Nothing special to comment on here:
#include <base/quadrature_lib.h> #include <base/function.h> #include <base/parameter_handler.h> #include <base/function_parser.h> #include <base/utilities.h> #include <base/conditional_ostream.h> #include <lac/vector.h> #include <lac/compressed_sparsity_pattern.h> #include <grid/tria.h> #include <grid/grid_generator.h> #include <grid/grid_out.h> #include <grid/grid_refinement.h> #include <grid/tria_accessor.h> #include <grid/tria_iterator.h> #include <grid/grid_in.h> #include <dofs/dof_handler.h> #include <dofs/dof_accessor.h> #include <dofs/dof_tools.h> #include <fe/fe_values.h> #include <fe/fe_system.h> #include <fe/mapping_q1.h> #include <fe/fe_q.h> #include <numerics/data_out.h> #include <numerics/vectors.h> #include <numerics/solution_transfer.h>
Then, as mentioned in the introduction, we use various Trilinos packages as linear solvers as well as for automatic differentiation. These are in the following include files.
In particular, Epetra is the basic trilinos vector/matrix library and comes with several header files pertaining to individual aspects of it that will become clear later on:
#include <Epetra_SerialComm.h> #include <Epetra_Map.h> #include <Epetra_CrsGraph.h> #include <Epetra_CrsMatrix.h> #include <Epetra_Vector.h>
Next, Teuchos is a Trilinos utility library that is used to set parameters within the Aztec solver library:
#include <Teuchos_ParameterList.hpp>
Aztec itself is the iterative solver library:
#include <AztecOO.h> #include <AztecOO_Operator.h>
Amesos is a direct solver package within Trilinos:
#include <Amesos.h>
Finally, Sacado is the automatic differentiation package, which is used to find the Jacobian for a fully implicit Newton iteration:
#include <Sacado.hpp>
And this again is C++, as well as a header file from BOOST that declares a class representing an array of fixed size:
#include <iostream> #include <fstream> #include <vector> #include <memory> #include <boost/array.hpp>
To end this section, introduce everythin in the dealii library into the current namespace:
using namespace dealii;
Here we define the flux function for this particular system of conservation laws, as well as pretty much everything else that's specific to the Euler equations for gas dynamics, for reasons discussed in the introduction. We group all this into a structure that defines everything that has to do with the flux. All members of this structure are static, i.e. the structure has no actual state specified by instance member variables. The better way to do this, rather than a structure with all static members would be to use a namespace -- but namespaces can't be templatized and we want some of the member variables of the structure to depend on the space dimension, which we in our usual way introduce using a template parameter.
template <int dim> struct EulerEquations {
First a few variables that describe the various components of our solution vector in a generic way. This includes the number of components in the system (Euler's equations have one entry for momenta in each spatial direction, plus the energy and density components, for a total of dim+2 components), as well as functions that describe the index within the solution vector of the first momentum component, the density component, and the energy density component. Note that all these numbers depend on the space dimension; defining them in a generic way (rather than by implicit convention) makes our code more flexible and makes it easier to later extend it, for example by adding more components to the equations.
static const unsigned int n_components = dim + 2; static const unsigned int first_momentum_component = 0; static const unsigned int density_component = dim; static const unsigned int energy_component = dim+1;
When generating graphical output way down in this program, we need to specify the names of the solution variables as well as how the various components group into vector and scalar fields. We could describe this there, but in order to keep things that have to do with the Euler equation localized here and the rest of the program as generic as possible, we provide this sort of information in the following two functions:
static std::vector<std::string> component_names () { std::vector<std::string> names (dim, "momentum"); names.push_back ("density"); names.push_back ("energy_density"); return names; } static std::vector<DataComponentInterpretation::DataComponentInterpretation> component_interpretation () { std::vector<DataComponentInterpretation::DataComponentInterpretation> data_component_interpretation (dim, DataComponentInterpretation::component_is_part_of_vector); data_component_interpretation .push_back (DataComponentInterpretation::component_is_scalar); data_component_interpretation .push_back (DataComponentInterpretation::component_is_scalar); return data_component_interpretation; }
Next, we define the gas constant. We will set it to 1.4 in its definition immediately following the declaration of this class (unlike integer variables, like the ones above, static const floating point member variables cannot be initialized within the class declaration in C++). This value of 1.4 is representative of a gas that consists of molecules composed of two atoms, such as air which consists up to small traces almost entirely of
and
.
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
(note that the independent variables contain the momentum components
, not the velocities
).
There is one slight problem: We will need to call the following functions with input arguments of type std::vector<number> and Vector<number>. The problem is that the former has an access operator operator[] whereas the latter, for historical reasons, has operator(). 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 v[i] or v(i), we can use *(v.begin() + i), i.e. we generate an iterator that points to the ith 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
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. Note that we know the size of the matrix: it has as many rows as the system has components, and dim columns; rather than using a FullMatrix object for such a matrix (which has a variable number of rows and columns and must therefore allocate memory on the heap each time such a matrix is created), we use a rectangular array of numbers right away.
We templatize the numerical type of the flux function so that we may use the automatic differentiation type here. Similarly, we will call the function with different input vector data types, so we templatize on it as well:
template <typename InputVector, typename number> static void compute_flux_matrix (const InputVector &W, number (&flux)[n_components][dim]) {
First compute the pressure that appears in the flux matrix, and then compute the first dim columns of the matrix that correspond to the momentum terms:
const number pressure = compute_pressure<number> (W); 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][d] += pressure; }
Then the terms for the density (i.e. mass conservation), and, lastly, conservation of energy:
for (unsigned int d=0; d<dim; ++d) flux[density_component][d] = W[first_momentum_component+d]; for (unsigned int d=0; d<dim; ++d) flux[energy_component][d] = W[first_momentum_component+d] / W[density_component] * (W[energy_component] + pressure); }
On the boundaries of the domain and across hanging nodes we use a numerical flux function to enforce boundary conditions. This routine is the basic Lax-Friedrich's flux with a stabilization parameter
. It's form has also been given already in the introduction:
template <typename InputVector> static void numerical_normal_flux (const Point<dim> &normal, const InputVector &Wplus, const InputVector &Wminus, const double alpha, Sacado::Fad::DFad<double> (&normal_flux)[n_components]) { Sacado::Fad::DFad<double> iflux[n_components][dim]; Sacado::Fad::DFad<double> oflux[n_components][dim]; compute_flux_matrix (Wplus, iflux); compute_flux_matrix (Wminus, oflux); 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]; normal_flux[di] += 0.5*alpha*(Wplus[di] - Wminus[di]); } }
In the same way as describing the flux function
, we also need to have a way to describe the right hand side forcing term. As mentioned in the introduction, we consider only gravity here, which leads to the specific form
, shown here for the 3d case. More specifically, we will consider only
in 3d, or
in 2d. This naturally leads to the following function:
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; } }
Another thing we have to deal with is boundary conditions. To this end, let us first define the kinds of boundary conditions we currently know how to deal with:
enum BoundaryKind
{
inflow_boundary,
outflow_boundary,
no_penetration_boundary,
pressure_boundary
};
The next part is to actually decide what to do at each kind of boundary. To this end, remember from the introduction that boundary conditions are specified by choosing a value
on the outside of a boundary given an inhomogeneity
and possibly the solution's value
on the inside. Both are then passed to the numerical flux
to define boundary contributions to the bilinear form.
Boundary conditions can in some cases be specified for each component of the solution vector independently. For example, if component
is marked for inflow, then
. If it is an outflow, then
. These two simple cases are handled first in the function below.
There is a little snag that makes this function unpleasant from a C++ language viewpoint: The output vector Wminus will of course be modified, so it shouldn't be a const argument. Yet it is in the implementation below, and needs to be in order to allow the code to compile. The reason is that we call this function at a place where Wminus is of type Table<2,Sacado::Fad::DFad<double> >, this being 2d table with indices representing the quadrature point and the vector component, respectively. We call this function with Wminus[q] as last argument; subscripting a 2d table yields a temporary accessor object representing a 1d vector, just what we want here. The problem is that a temporary accessor object can't be bound to a non-const reference argument of a function, as we would like here, according to the C++ 1998 and 2003 standards (something that will be fixed with the next standard in the form of rvalue references). We get away with making the output argument here a constant because it is the accessor object that's constant, not the table it points to: that one can still be written to. The hack is unpleasant nevertheless because it restricts the kind of data types that may be used as template argument to this function: a regular vector isn't going to do because that one can not be written to when marked const. With no good solution around at the moment, we'll go with the pragmatic, even if not pretty, solution shown here:
template <typename DataVector> static void compute_Wminus (const BoundaryKind (&boundary_kind)[n_components], const Point<dim> &normal_vector, const DataVector &Wplus, const Vector<double> &boundary_values, const DataVector &Wminus) { for (unsigned int c = 0; c < n_components; c++) switch (boundary_kind[c]) { case inflow_boundary: { Wminus[c] = boundary_values(c); break; } case outflow_boundary: { Wminus[c] = Wplus[c]; break; }
Prescribed pressure boundary conditions are a bit more complicated by the fact that even though the pressure is prescribed, we really are setting the energy component here, which will depend on velocity and pressure. So even though this seems like a Dirichlet type boundary condition, we get sensitivities of energy to velocity and density (unless these are also prescribed):
case pressure_boundary: { const typename DataVector::value_type density = (boundary_kind[density_component] == inflow_boundary ? boundary_values(density_component) : Wplus[density_component]); typename DataVector::value_type kinetic_energy = 0; for (unsigned int d=0; d<dim; ++d) if (boundary_kind[d] == inflow_boundary) kinetic_energy += boundary_values(d)*boundary_values(d); else kinetic_energy += Wplus[d]*Wplus[d]; kinetic_energy *= 1./2./density; Wminus[c] = boundary_values(c) / (gas_gamma-1.0) + kinetic_energy; break; } case no_penetration_boundary: {
We prescribe the velocity (we are dealing with a particular component here so that the average of the velocities is orthogonal to the surface normal. This creates sensitivies of across the velocity components.
Sacado::Fad::DFad<double> vdotn = 0;
for (unsigned int d = 0; d < dim; d++) {
vdotn += Wplus[d]*normal_vector[d];
}
Wminus[c] = Wplus[c] - 2.0*vdotn*normal_vector[c];
break;
}
default:
Assert (false, ExcNotImplemented());
}
}
In this class, we also want to specify how to refine the mesh. The class ConservationLaw that will use all the information we provide here in the EulerEquation class is pretty agnostic about the particular conservation law it solves: as doesn't even really care how many components a solution vector has. Consequently, it can't know what a reasonable refinement indicator would be. On the other hand, here we do, or at least we can come up with a reasonable choice: we simply look at the gradient of the density, and compute
, where
is the center of cell
.
There are certainly a number of equally reasonable refinement indicators, but this one does, and it is easy to compute:
static void compute_refinement_indicators (const DoFHandler<dim> &dof_handler, const Mapping<dim> &mapping, const Vector<double> &solution, Vector<double> &refinement_indicators) { const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; std::vector<unsigned int> dofs (dofs_per_cell); const QMidpoint<dim> quadrature_formula; const UpdateFlags update_flags = update_gradients; FEValues<dim> fe_v (mapping, dof_handler.get_fe(), quadrature_formula, update_flags); std::vector<std::vector<Tensor<1,dim> > > dU (1, std::vector<Tensor<1,dim> >(n_components)); typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) { fe_v.reinit(cell); fe_v.get_function_grads (solution, dU); refinement_indicators(cell_no) = std::log(1+ std::sqrt(dU[0][density_component] * dU[0][density_component])); } }
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
, density
, and energy density
. What we would like to also put into our output file are velocities
and pressure
.
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" is 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 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
; obviously,
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 duh vector only contains data if we really need it (the system knows about this because we say so in the get_needed_update_flags() 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 dim 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 first_momentum_component and density_component information:
for (unsigned int q=0; q<n_quadrature_points; ++q) { const double density = uh[q](density_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 (dim, 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; }
Our next job is to define a few classes that will contain run-time parameters (for example solver tolerances, number of iterations, stabilization parameter, and the like). One could do this in the main class, but we separate it from that one to make the program more modular and easier to read: Everything that has to do with run-time parameters will be in the following namespace, whereas the program logic is in the main class.
We will split the run-time parameters into a few separate structures, which we will all put into a namespace Parameters. Of these classes, there are a few that group the parameters for individual groups, such as for solvers, mesh refinement, or output. Each of these classes have functions declare_parameters() and parse_parameters() that declare parameter subsections and entries in a ParameterHandler object, and retrieve actual parameter values from such an object, respectively. These classes declare all their parameters in subsections of the ParameterHandler.
The final class of the following namespace combines all the previous classes by deriving from them and taking care of a few more entries at the top level of the input file, as well as a few odd other entries in subsections that are too short to warrent a structure by themselves.
It is worth pointing out one thing here: None of the classes below have a constructor that would initialize the various member variables. This isn't a problem, however, since we will read all variables declared in these classes from the input file (or indirectly: a ParameterHandler object will read it from there, and we will get the values from this object), and they will be initialized this way. In case a certain variable is not specified at all in the input file, this isn't a problem either: The ParameterHandler class will in this case simply take the default value that was specified when declaring an entry in the declare_parameters() functions of the classes below.
namespace Parameters
{
The first of these classes deals with parameters for the linear inner solver. It offers parameters that indicate which solver to use (GMRES as a solver for general non-symmetric indefinite systems, or a sparse direct solver), the amount of output to be produced, as well as various parameters that tweak the thresholded incomplete LU decomposition (ILUT) that we use as a preconditioner for GMRES.
In particular, the ILUT takes the following parameters:
.
The meaning of each parameter is also briefly described in the third argument of the ParameterHandler::declare_entry call in declare_parameters().
struct Solver { enum SolverType { gmres, direct }; SolverType solver; enum OutputType { quiet, verbose }; OutputType output; double linear_residual; int max_iterations; double ilut_fill; double ilut_atol; double ilut_rtol; double ilut_drop; static void declare_parameters (ParameterHandler &prm); void parse_parameters (ParameterHandler &prm); }; void Solver::declare_parameters (ParameterHandler &prm) { prm.enter_subsection("linear solver"); { prm.declare_entry("output", "quiet", Patterns::Selection("quiet|verbose"), "State whether output from solver runs should be printed. " "Choices are <quiet|verbose>."); prm.declare_entry("method", "gmres", Patterns::Selection("gmres|direct"), "The kind of solver for the linear system. " "Choices are <gmres|direct>."); prm.declare_entry("residual", "1e-10", Patterns::Double(), "Linear solver residual"); prm.declare_entry("max iters", "300", Patterns::Integer(), "Maximum solver iterations"); prm.declare_entry("ilut fill", "2", Patterns::Double(), "Ilut preconditioner fill"); prm.declare_entry("ilut absolute tolerance", "1e-9", Patterns::Double(), "Ilut preconditioner tolerance"); prm.declare_entry("ilut relative tolerance", "1.1", Patterns::Double(), "Ilut relative tolerance"); prm.declare_entry("ilut drop tolerance", "1e-10", Patterns::Double(), "Ilut drop tolerance"); } prm.leave_subsection(); } void Solver::parse_parameters (ParameterHandler &prm) { prm.enter_subsection("linear solver"); { const std::string op = prm.get("output"); if (op == "verbose") output = verbose; if (op == "quiet") output = quiet; const std::string sv = prm.get("method"); if (sv == "direct") solver = direct; else if (sv == "gmres") solver = gmres; linear_residual = prm.get_double("residual"); max_iterations = prm.get_integer("max iters"); ilut_fill = prm.get_double("ilut fill"); ilut_atol = prm.get_double("ilut absolute tolerance"); ilut_rtol = prm.get_double("ilut relative tolerance"); ilut_drop = prm.get_double("ilut drop tolerance"); } prm.leave_subsection(); }
Similarly, here are a few parameters that determine how the mesh is to be refined (and if it is to be refined at all). For what exactly the shock parameters do, see the mesh refinement functions further down.
struct Refinement { bool do_refine; double shock_val; double shock_levels; static void declare_parameters (ParameterHandler &prm); void parse_parameters (ParameterHandler &prm); }; void Refinement::declare_parameters (ParameterHandler &prm) { prm.enter_subsection("refinement"); { prm.declare_entry("refinement", "true", Patterns::Bool(), "Whether to perform mesh refinement or not"); prm.declare_entry("refinement fraction", "0.1", Patterns::Double(), "Fraction of high refinement"); prm.declare_entry("unrefinement fraction", "0.1", Patterns::Double(), "Fraction of low unrefinement"); prm.declare_entry("max elements", "1000000", Patterns::Double(), "maximum number of elements"); prm.declare_entry("shock value", "4.0", Patterns::Double(), "value for shock indicator"); prm.declare_entry("shock levels", "3.0", Patterns::Double(), "number of shock refinement levels"); } prm.leave_subsection(); } void Refinement::parse_parameters (ParameterHandler &prm) { prm.enter_subsection("refinement"); { do_refine = prm.get_bool ("refinement"); shock_val = prm.get_double("shock value"); shock_levels = prm.get_double("shock levels"); } prm.leave_subsection(); }
Next a section on flux modifications to make it more stable. In particular, two options are offered to stabilize the Lax-Friedrichs flux: either choose
where
is either a fixed number specified in the input file, or where
is a mesh dependent value. In the latter case, it is chosen as
with
the diameter of the face to which the flux is applied, and
the current time step.
struct Flux { enum StabilizationKind { constant, mesh_dependent }; StabilizationKind stabilization_kind; double stabilization_value; static void declare_parameters (ParameterHandler &prm); void parse_parameters (ParameterHandler &prm); }; void Flux::declare_parameters (ParameterHandler &prm) { prm.enter_subsection("flux"); { prm.declare_entry("stab", "mesh", Patterns::Selection("constant|mesh"), "Whether to use a constant stabilization parameter or " "a mesh-dependent one"); prm.declare_entry("stab value", "1", Patterns::Double(), "alpha stabilization"); } prm.leave_subsection(); } void Flux::parse_parameters (ParameterHandler &prm) { prm.enter_subsection("flux"); { const std::string stab = prm.get("stab"); if (stab == "constant") stabilization_kind = constant; else if (stab == "mesh") stabilization_kind = mesh_dependent; else AssertThrow (false, ExcNotImplemented()); stabilization_value = prm.get_double("stab value"); } prm.leave_subsection(); }
Then a section on output parameters. We offer to produce Schlieren plots (the squared gradient of the density, a tool to visualize shock fronts), and a time interval between graphical output in case we don't want an output file every time step.
struct Output { bool schlieren_plot; double output_step; static void declare_parameters (ParameterHandler &prm); void parse_parameters (ParameterHandler &prm); }; void Output::declare_parameters (ParameterHandler &prm) { prm.enter_subsection("output"); { prm.declare_entry("schlieren plot", "true", Patterns::Bool (), "Whether or not to produce schlieren plots"); prm.declare_entry("step", "-1", Patterns::Double(), "Output once per this period"); } prm.leave_subsection(); } void Output::parse_parameters (ParameterHandler &prm) { prm.enter_subsection("output"); { schlieren_plot = prm.get_bool("schlieren plot"); output_step = prm.get_double("step"); } prm.leave_subsection(); }
Finally the class that brings it all together. It declares a number of parameters itself, mostly ones at the top level of the parameter file as well as several in section too small to warrant their own classes. It also contains everything that is actually space dimension dependent, like initial or boundary conditions.
Since this class is derived from all the ones above, the declare_parameters() and parse_parameters() functions call the respective functions of the base classes as well.
Note that this class also handles the declaration of initial and boundary conditions specified in the input file. To this end, in both cases, there are entries like "w_0 value" which represent an expression in terms of
that describe the initial or boundary condition as a formula that will later be parsed by the FunctionParser class. Similar expressions exist for "w_1", "w_2", etc, denoting the dim+2 conserved variables of the Euler system. Similarly, we allow up to max_n_boundaries boundary indicators to be used in the input file, and each of these boundary indicators can be associated with an inflow, outflow, or pressure boundary condition, with inhomogenous boundary conditions being specified for each component and each boundary indicator separately.
The data structure used to store the boundary indicators is a bit complicated. It is an array of max_n_boundaries elements indicating the range of boundary indicators that will be accepted. For each entry in this array, we store a pair of data in the BoundaryCondition structure: first, an array of size n_components that for each component of the solution vector indicates whether it is an inflow, outflow, or other kind of boundary, and second a FunctionParser object that describes all components of the solution vector for this boundary id at once.
The BoundaryCondition structure requires a constructor since we need to tell the function parser object at construction time how many vector components it is to describe. This initialization can therefore not wait till we actually set the formulas the FunctionParser object represents later in AllParameters::parse_parameters()
For the same reason of having to tell Function objects their vector size at construction time, we have to have a constructor of the AllParameters class that at least initializes the other FunctionParser object, i.e. the one describing initial conditions.
template <int dim> struct AllParameters : public Solver, public Refinement, public Flux, public Output { static const unsigned int max_n_boundaries = 10; struct BoundaryConditions { typename EulerEquations<dim>::BoundaryKind kind[EulerEquations<dim>::n_components]; FunctionParser<dim> values; BoundaryConditions (); }; AllParameters (); double diffusion_power; double time_step, final_time; double theta; bool is_stationary; std::string mesh_filename; FunctionParser<dim> initial_conditions; BoundaryConditions boundary_conditions[max_n_boundaries]; static void declare_parameters (ParameterHandler &prm); void parse_parameters (ParameterHandler &prm); }; template <int dim> AllParameters<dim>::BoundaryConditions::BoundaryConditions () : values (EulerEquations<dim>::n_components) {} template <int dim> AllParameters<dim>::AllParameters () : initial_conditions (EulerEquations<dim>::n_components) {} template <int dim> void AllParameters<dim>::declare_parameters (ParameterHandler &prm) { prm.declare_entry("mesh", "grid.inp", Patterns::Anything(), "intput file name"); prm.declare_entry("diffusion power", "2.0", Patterns::Double(), "power of mesh size for diffusion"); prm.enter_subsection("time stepping"); { prm.declare_entry("time step", "0.1", Patterns::Double(0), "simulation time step"); prm.declare_entry("final time", "10.0", Patterns::Double(0), "simulation end time"); prm.declare_entry("theta scheme value", "0.5", Patterns::Double(0,1), "value for theta that interpolated between explicit " "Euler (theta=0), Crank-Nicolson (theta=0.5), and " "implicit Euler (theta=1)."); } prm.leave_subsection(); for (unsigned int b=0; b<max_n_boundaries; ++b) { prm.enter_subsection("boundary_" + Utilities::int_to_string(b)); { prm.declare_entry("no penetration", "false", Patterns::Bool(), "whether the named boundary allows gas to " "penetrate or is a rigid wall"); for (unsigned int di=0; di<EulerEquations<dim>::n_components; ++di) { prm.declare_entry("w_" + Utilities::int_to_string(di), "outflow", Patterns::Selection("inflow|outflow|pressure"), "<inflow|outflow|pressure>"); prm.declare_entry("w_" + Utilities::int_to_string(di) + " value", "0.0", Patterns::Anything(), "expression in x,y,z"); } } prm.leave_subsection(); } prm.enter_subsection("initial condition"); { for (unsigned int di=0; di<EulerEquations<dim>::n_components; ++di) prm.declare_entry("w_" + Utilities::int_to_string(di) + " value", "0.0", Patterns::Anything(), "expression in x,y,z"); } prm.leave_subsection(); Parameters::Solver::declare_parameters (prm); Parameters::Refinement::declare_parameters (prm); Parameters::Flux::declare_parameters (prm); Parameters::Output::declare_parameters (prm); } template <int dim> void AllParameters<dim>::parse_parameters (ParameterHandler &prm) { mesh_filename = prm.get("mesh"); diffusion_power = prm.get_double("diffusion power"); prm.enter_subsection("time stepping"); { time_step = prm.get_double("time step"); if (time_step == 0) { is_stationary = true; time_step = 1.0; final_time = 1.0; } else is_stationary = false; final_time = prm.get_double("final time"); theta = prm.get_double("theta scheme value"); } prm.leave_subsection(); for (unsigned int boundary_id=0; boundary_id<max_n_boundaries; ++boundary_id) { prm.enter_subsection("boundary_" + Utilities::int_to_string(boundary_id)); { std::vector<std::string> expressions(EulerEquations<dim>::n_components, "0.0"); const bool no_penetration = prm.get_bool("no penetration"); for (unsigned int di=0; di<EulerEquations<dim>::n_components; ++di) { const std::string boundary_type = prm.get("w_" + Utilities::int_to_string(di)); if ((di < dim) && (no_penetration == true)) boundary_conditions[boundary_id].kind[di] = EulerEquations<dim>::no_penetration_boundary; else if (boundary_type == "inflow") boundary_conditions[boundary_id].kind[di] = EulerEquations<dim>::inflow_boundary; else if (boundary_type == "pressure") boundary_conditions[boundary_id].kind[di] = EulerEquations<dim>::pressure_boundary; else if (boundary_type == "outflow") boundary_conditions[boundary_id].kind[di] = EulerEquations<dim>::outflow_boundary; else AssertThrow (false, ExcNotImplemented()); expressions[di] = prm.get("w_" + Utilities::int_to_string(di) + " value"); } boundary_conditions[boundary_id].values .initialize (FunctionParser<dim>::default_variable_names(), expressions, std::map<std::string, double>()); } prm.leave_subsection(); } prm.enter_subsection("initial condition"); { std::vector<std::string> expressions (EulerEquations<dim>::n_components, "0.0"); for (unsigned int di = 0; di < EulerEquations<dim>::n_components; di++) expressions[di] = prm.get("w_" + Utilities::int_to_string(di) + " value"); initial_conditions.initialize (FunctionParser<dim>::default_variable_names(), expressions, std::map<std::string, double>()); } prm.leave_subsection(); Parameters::Solver::parse_parameters (prm); Parameters::Refinement::parse_parameters (prm); Parameters::Flux::parse_parameters (prm); Parameters::Output::parse_parameters (prm); } }
Here finally comes the class that actually does something with all the Euler equation and parameter specifics we've defined above. The public interface is pretty much the same as always (the constructor now takes the name of a file from which to read parameters, which is passed on the command line). The private function interface is also pretty similar to the usual arrangement, with the assemble_system function split into three parts: one that contains the main loop over all cells and that then calls the other two for integrals over cells and faces, respectively.
template <int dim> class ConservationLaw { public: ConservationLaw (const char *input_filename); void run (); private: void setup_system (); void assemble_system (); void assemble_cell_term (const FEValues<dim> &fe_v, const std::vector<unsigned int> &dofs); void assemble_face_term (const unsigned int face_no, const FEFaceValuesBase<dim> &fe_v, const FEFaceValuesBase<dim> &fe_v_neighbor, const std::vector<unsigned int> &dofs, const std::vector<unsigned int> &dofs_neighbor, const bool external_face, const unsigned int boundary_id, const double face_diameter); std::pair<unsigned int, double> solve (Vector<double> &solution); void compute_refinement_indicators (Vector<double> &indicator) const; void refine_grid (const Vector<double> &indicator); void output_results () const;
The first few member variables are also rather standard. Note that we define a mapping object to be used throughout the program when assembling terms (we will hand it to every FEValues and FEFaceValues object); the mapping we use is just the standard
mapping -- nothing fancy, in other words -- but declaring one here and using it throughout the program will make it simpler later on to change it if that should become necessary. This is, in fact, rather pertinent: it is known that for transsonic simulations with the Euler equations, computations do not converge even as
if the boundary approximation is not of sufficiently high order.
Triangulation<dim> triangulation; const MappingQ1<dim> mapping; const FESystem<dim> fe; DoFHandler<dim> dof_handler; const QGauss<dim> quadrature; const QGauss<dim-1> face_quadrature;
Next come a number of data vectors that correspond to the solution of the previous time step (old_solution), the best guess of the current solution (current_solution; we say guess because the Newton iteration to compute it may not have converged yet, whereas old_solution refers to the fully converged final result of the previous time step), and a predictor for the solution at the next time step, computed by extrapolating the current and previous solution one time step into the future:
Vector<double> old_solution; Vector<double> current_solution; Vector<double> predictor; Vector<double> right_hand_side;
This final set of member variables (except for the object holding all run-time parameters at the very bottom and a screen output stream that only prints something if verbose output has been requested) deals with the interface we have in this program to the Trilinos library that provides us with linear solvers.
Trilinos is designed to be a library that also runs in parallel on distributed memory systems, so matrices and vectors need two things: (i) a communicator object that facilitates sending messages to remote machines, and (ii) a description which elements of a vector or matrix reside locally on a machine and which are stored remotely.
We do not actually run the current program in parallel, and so the objects we use here are pretty much dummy objects for this purpose: the communicator below represents a system that includes only a single machine, and the index map encodes that all elements are stored locally. Nevertheless, we need them.
Furthermore, we need a matrix object for the system matrix to be used in each Newton step. Note that map and matrix need to be updated for their sizes whenever we refine the mesh. In Trilinos, this is easiest done by simply deleting the previous object and creating a new one. To minimize hassle and avoid memory leaks, we use a std::auto_ptr instead of a plain pointer for this.
Epetra_SerialComm communicator;
std::auto_ptr<Epetra_Map> Map;
std::auto_ptr<Epetra_CrsMatrix> Matrix;
Parameters::AllParameters<dim> parameters;
ConditionalOStream verbose_cout;
};
There is nothing much to say about the constructor. Essentially, it reads the input file and fills the parameter object with the parsed values:
template <int dim> ConservationLaw<dim>::ConservationLaw (const char *input_filename) : mapping (), fe (FE_Q<dim>(1), EulerEquations<dim>::n_components), dof_handler (triangulation), quadrature (2), face_quadrature (2), verbose_cout (std::cout, false) { ParameterHandler prm; Parameters::AllParameters<dim>::declare_parameters (prm); prm.read_input (input_filename); parameters.parse_parameters (prm); verbose_cout.set_condition (parameters.output == Parameters::Solver::verbose); }
The following function is called each time the mesh is changed. Essentially what it does is to resize the Trilinos matrix. In addition to just resizing it, it also builds a sparsity pattern, initializes the row lengths of the matrix with the ones from this sparsity pattern, and finally puts zero entries into the places where nonzero entries will later be found. This will make subsequent operations on the matrix faster, because no new memory will need to be allocated:
template <int dim> void ConservationLaw<dim>::setup_system () { Map.reset (new Epetra_Map(dof_handler.n_dofs(), 0, communicator));
Now create a sparsity pattern, condense it, and count the number of nonzero entries per row:
CompressedSparsityPattern sparsity_pattern (dof_handler.n_dofs(), dof_handler.n_dofs()); DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); sparsity_pattern.compress(); std::vector<int> row_lengths (dof_handler.n_dofs()); for (unsigned int i=0; i<dof_handler.n_dofs(); ++i) row_lengths[i] = sparsity_pattern.row_length (i);
Next we build the matrix, using the constructor that optimizes with the existing lengths per row variable. After this, loop over the individual rows of the deal.II sparsity pattern and create entries in the Trilinos matrix in the corresponding places. At the end, call the FillComplete() function that indicates that no other matrix entries will be needed:
Matrix.reset (new Epetra_CrsMatrix(Copy, *Map, &row_lengths[0], true)); const unsigned int max_nonzero_entries = *std::max_element (row_lengths.begin(), row_lengths.end()); std::vector<double> values(max_nonzero_entries, 0); std::vector<int> row_indices(max_nonzero_entries); for (unsigned int row=0; row<do