Reference documentation for deal.II version GIT relicensing-1182-g1782d71672 2024-07-22 00:00:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
The 'Traveling-wave solutions of a qualitative model for combustion waves.' code gallery program

This program was contributed by Shamil Magomedov <magomedov.shamil.m@gmail.com>.
It comes without any warranty or support by its authors or the authors of deal.II.

This program is part of the deal.II code gallery and consists of the following files (click to inspect):

Pictures from this code gallery program

Annotated version of README.md

Traveling-wave solutions of a qualitative model for combustion waves

This program demonstrates the use of adaptive finite elements to compute traveling-wave solutions of partial differential equations. One of the challenges in solving this type of problem is the presence of an unknown wave speed. Below we show how to overcome this.

Building, compiling and running

To run the program, enter the following commands from the current directory:

mkdir build && cd build
cmake ..
make
./main ../ParametersList.prm

Problem

To illustrate the algorithm for the computation of traveling-wave profiles, we consider a combustion model described in [1]. In a moving reference frame of a steadily propagating wave we have the following nondimensionalized system:

\begin{align*} - c u_{\xi} + (1 + \epsilon u) u_{\xi} &= -\dfrac{\epsilon}{2} T_{\xi} + \dfrac{4 \delta}{3 \epsilon} \Pr u_{\xi \xi}, \\ - c(T_{\xi} - u_{\xi}) &= q \omega + \delta T_{\xi \xi}, \\ - c \lambda_{\xi} &= \omega + \dfrac{\delta}{\text{Le}} \lambda_{\xi \xi}. \end{align*}

Here, \(u\) is the pressure, \(T\) is the temperature, \(\lambda\) is the reaction-progress variable, varying from \(0\) in the fresh mixture to \(1\) in the burnt products, \(c > 0\) is the unknown wave speed, \(\Pr\) and \(\mathrm{Le}\) are the Prandtl and Lewis numbers, \(q\) is the energy of heat release. The model parameters \(\epsilon\) and \(\delta\) determine the strength of nonlinearity and dissipative effects. The reaction rate \(\omega\) is taken as

\begin{align*} \omega = k (1 - \lambda) \exp(-\theta / T) \, \mathrm{H}(T - T_{\mathrm{ign}}), \end{align*}

with activation energy \(\theta\), ignition temperature \(T_{\mathrm{ign}}\), constant of the reaction rate \(k\), and the Heaviside step function \(\mathrm{H}\).

The boundary conditions at \(\xi = -\infty\) are

\begin{align*} u_{\xi} = 0, \ T = T_l, \ \lambda = 1, \end{align*}

and at \(\xi = +\infty\) are

\begin{align*} u = u_r, \ \lambda = 0. \end{align*}

The right boundary condition for temperature is \(T = T_r\) for detonation waves (supersonic regime, i.e. \(c > 1\)) and \(T_{\xi} = 0\) for deflagration waves (subsonic regime, i.e. \(c < 1\)).

Because of translational invariance, we need to impose another constraint on the system to fix a particular solution. So we choose the following centering condition: \(T(0) = T_{\mathrm{ign}}\).

Numerical algorithm

Newton–Raphson iteration scheme

The nonlinear boundary value problem is solved numerically on a finite interval \(I = [l, r]\) \(\left(|l|, |r| \gg 1 \right)\), using a Newton–Raphson iteration scheme, similar to one, described in deal.II tutorial step-15. The main difference from step-15 is that we have an additional scalar unknown, the front velocity \(c\). So the algorithm has to be modified to take this feature into account.

Rewriting the system in a vector form

\begin{align*} \mathbf{F}(u, T, \lambda, c) = \left(\begin{array}{c} \dfrac{4 \delta}{3 \epsilon} \Pr u_{\xi \xi} - u_{\xi} (1 - c + \epsilon u) -\dfrac{\epsilon}{2} T_{\xi} \\[7pt] \delta T_{\xi \xi} + c(T_{\xi} - u_{\xi}) + q \omega \\[5pt] \dfrac{\delta}{\text{Le}} \lambda_{\xi \xi} + c \lambda_{\xi} + \omega \end{array} \right) = 0, \end{align*}

we define a Newton–Raphson iteration as

\begin{align*} \mathbf{F'}(\mathbf{x}^k, \mathbf{dx}^k) = - \mathbf{F}(\mathbf{x}^k), \end{align*}

\begin{align*} \mathbf{x}^{n+1} = \mathbf{x}^{n} + \alpha^k \mathbf{dx}^k, \end{align*}

where \(k\) is the step number, \(\mathbf{x}^k = (u^k, T^k, \lambda^k, c^k)^{\top}\) is a vector argument, \(\mathbf{dx}^k = (du^k, dT^k, d\lambda^k, dc^k)^{\top}\) is an increment, \(\alpha^k\) is some damping parameter for managing the global convergence behavior and \(\mathbf{F'}(\mathbf{x}^k, \mathbf{dx}^k)\) is the directional derivative, defined as

\begin{align*} \mathbf{F'}(\mathbf{x}, \mathbf{dx}) = \dfrac{\mathrm{d}}{\mathrm{d} \varepsilon} \Big|_{\varepsilon=0} \mathbf{F}(\mathbf{x} &+ \varepsilon \mathbf{dx}). \end{align*}

The system to be solved at every iteration step to obtain the increment \(\mathbf{dx}^k = (du^k, dT^k, d\lambda^k, dc^k)^{\top}\) can be represented in matrix vector notation as follows

\begin{align*} \begin{pmatrix} \dfrac{4 \delta}{3 \epsilon} \Pr \partial_{\xi \xi} - (1 - c + \epsilon u)\partial_{\xi} - \epsilon u_{\xi} & -\dfrac{\epsilon}{2} \partial_{\xi} & 0 & u_{\xi} \\[9pt] -c \partial_{\xi} & \delta \partial_{\xi \xi} + c \partial_{\xi} + q \kappa_1 & q \kappa_2 & T_{\xi} - u_{\xi} \\[9pt] 0 & \kappa_1 & \dfrac{\delta}{\text{Le}} \partial_{\xi \xi} + c \partial_{\xi} + \kappa_2 & \lambda_{\xi} \end{pmatrix} \begin{pmatrix} du \\[9pt] dT \\[9pt] d\lambda \\[9pt] dc \end{pmatrix} = -\begin{pmatrix} f_1 \\[9pt] f_2 \\[9pt] f_3 \end{pmatrix}, \end{align*}

where

\begin{align*} \kappa_1 &= k (1 - \lambda) \exp(-\theta / T) \left[ \dfrac{\theta }{T^2} \, \text{H}(T - T_{\text{ign}}) + \delta(T - T_{\text{ign}}) \right], \\ \kappa_2 &= - k \exp(-\theta / T) \, \text{H}(T - T_{\text{ign}}), \end{align*}

in which \(\delta(\cdot)\) is a Dirac delta function, and \(f_i \, (i=1,2,3)\) are the components of the vector function \(\mathbf{F}(u, T, \lambda, c)\). The term \(\delta(T - T_{\text{ign}})\) can be rewritten as

\begin{align*} \delta(T - T_{\text{ign}}) = \frac{\delta(\xi)}{|T'(0)|}. \end{align*}

We choose the initial guess \(\mathbf{x}^0\) to include the appropriate boundary values, therefore the update \(\mathbf{dx}^k\) uses homogeneous Dirichlet or Neumann boundary conditions.

Weak formulation

We multiply both sides of the equation for \(\mathbf{dx}^k\) with vector valued test function \(\mathbf{v} = (v_1, v_2, v_3)^{\top}\) and integrate over the domain \(\Omega\) to obtain a scalar equation

\begin{align*} J(\mathbf{dx}, \mathbf{v}) = -b(\mathbf{v}). \end{align*}

\begin{align*} J(\mathbf{dx}, &\mathbf{v}) = \int \limits_{\Omega} \mathbf{v} \cdot \mathbf{F'}(\mathbf{x}, \mathbf{dx}) \, d\xi = \\ = &\dfrac{4 \delta}{3 \epsilon} \Pr (-\partial_{\xi} v_1, \partial_{\xi} du) + (v_1, - (1 - c + \epsilon u)\partial_{\xi} du - \epsilon u_{\xi} du -\dfrac{\epsilon}{2} \partial_{\xi} dT + u_{\xi} dc) + \\ &+ \delta (-\partial_{\xi} v_2, \partial_{\xi} dT) + (v_2, -c \, \partial_{\xi} du + c \, \partial_{\xi} dT + q \kappa_1 dT + q \kappa_2 d\lambda + T_{\xi} dc - u_{\xi} dc) + \\ &+ \dfrac{\delta}{\text{Le}} (-\partial_{\xi} v_3, \partial_{\xi} d\lambda) + (v_3, \kappa_1 dT + c \partial_{\xi} d\lambda + \kappa_2 d\lambda + \lambda_{\xi} dc). \end{align*}

\begin{align*} b(\mathbf{v}) = \int \limits_{\Omega} &\mathbf{v} \cdot \mathbf{F}(\mathbf{x}) \, d\xi = \\ = &\dfrac{4 \delta}{3 \epsilon} \Pr (-\partial_{\xi} v_1, u_{\xi}) + (v_1, - u_{\xi} (1 - c + \epsilon u) -\dfrac{\epsilon}{2} T_{\xi}) + \\ &+ \delta (-\partial_{\xi} v_2, T_{\xi}) + (v_2, c(T_{\xi} - u_{\xi}) + q \omega) + \\ &+ \dfrac{\delta}{\text{Le}} (-\partial_{\xi} v_3, \lambda_{\xi}) + (v_3, c \lambda_{\xi} + \omega). \end{align*}

In the above expressions second derivatives disappear due to integration by parts with homogeneous Dirichlet and Neumann boundary conditions. The solution is sought as an expansion

\begin{align*} \begin{pmatrix} du \\[9pt] dT \\[9pt] d\lambda \\[9pt] dc \end{pmatrix} = \sum \limits_{i = 1}^{3N} U_{i}\begin{pmatrix} \phi_i^1 \\[9pt] \phi_i^2 \\[9pt] \phi_i^3 \\[9pt] 0 \end{pmatrix} + U_{3N + 1}\begin{pmatrix} 0 \\[9pt] 0 \\[9pt] 0 \\[9pt] 1 \end{pmatrix} \quad \in \quad V_p \times V_p \times V_p \times \mathbb{R}. \end{align*}

where \(V_p\) is a finite element space of continuous, piecewise polynomials of degree \(p\). The set of vector functions \((\phi_i^1, \phi_i^2, \phi_i^3)^{\top} \in V_p^3\) form the basis of the corresponding space. We then choose test functions \(\mathbf{v}\) to be the same as the basis functions, and obtain the linear system \(J U = b\). Elements of the matrix and the right-hand side are computed according to following formulas:

\begin{align*} J_{ij} = &\dfrac{4 \delta}{3 \epsilon} \Pr (-\partial_{\xi} \phi_i^1, \partial_{\xi} \phi_j^1) + (\phi_i^1, - (1 - c + \epsilon u)\partial_{\xi} \phi_j^1 - \epsilon u_{\xi} \phi_j^1 -\dfrac{\epsilon}{2} \partial_{\xi} \phi_j^2) + \\ &+ \delta (-\partial_{\xi} \phi_i^2, \partial_{\xi} \phi_j^2) + (\phi_i^2, -c \, \partial_{\xi} \phi_j^1 + c \, \partial_{\xi} \phi_j^2 + q \kappa_1 \phi_j^2 + q \kappa_2 \phi_j^3) + \\ &+ \dfrac{\delta}{\text{Le}} (-\partial_{\xi} \phi_i^3, \partial_{\xi} \phi_j^3) + (\phi_i^3, \kappa_1 \phi_j^2 + c \partial_{\xi} \phi_j^3 + \kappa_2 \phi_j^3), \end{align*}

\begin{align*} J_{i, 3N + 1} = (\phi_i^1, u_{\xi}) + (\phi_i^2, T_{\xi} - u_{\xi}) + (\phi_i^3, \lambda_{\xi}), \end{align*}

\begin{align*} b_{i} = &\dfrac{4 \delta}{3 \epsilon} \Pr (-\partial_{\xi} \phi_i^1, u_{\xi}) + (\phi_i^1, - u_{\xi} (1 - c + \epsilon u) -\dfrac{\epsilon}{2} T_{\xi}) + \\ &+ \delta (-\partial_{\xi} \phi_i^2, T_{\xi}) + (\phi_i^2, c(T_{\xi} - u_{\xi}) + q \omega) + \\ &+ \dfrac{\delta}{\text{Le}} (-\partial_{\xi} \phi_i^3, \lambda_{\xi}) + (\phi_i^3, c \lambda_{\xi} + \omega), \end{align*}

for \(i, j < 3N + 1\).

In order for the system to have a unique solution, we need to supplement it with one more equation, corresponding to the constraint \(T(0) = T_{\text{ign}}\). The initial approximation to the solution is set so as to satisfy this condition, so we just need the computed increment function \(dT\) to be zero at the specified point. Thus, we add a row of zeros with a value of 1 in the position corresponding to \(dT(0)\) to the matrix \(J\) and set \(b_{3N + 1} = 0\).

The resulting sparsity pattern structure has the form shown in the figure below. Sparsity pattern corresponding to basis polynomial degree @f$p=1@f$.

The integration of the terms with \(\kappa_1\) need special attention because of the Dirac delta function. If \(U_n\) and \(U_m\) are the degrees of freedom, associated with the vertex \(\xi = 0\) (i.e., \(\phi_n^2(0) = 1\) and \(\phi_m^3(0) = 1\)), we get

\begin{align*} (\phi_n^2, q (k (1 - \lambda) \exp(-\theta / T) \delta(T - T_{\text{ign}})) \phi_n^2) = \dfrac{q k (1 - \lambda(0)) \exp(-\theta / T(0))}{|T'(0)|} \end{align*}

and

\begin{align*} (\phi_m^3, (k (1 - \lambda) \exp(-\theta / T) \delta(T - T_{\text{ign}})) \phi_n^2) = \dfrac{k (1 - \lambda(0)) \exp(-\theta / T(0))}{|T'(0)|}. \end{align*}

Initial guess

The initial guess for detonation wave is obtained from the following problem

\begin{align*} u_{\xi} (- c + 1 + \epsilon u) &= -\dfrac{\epsilon}{2} T_{\xi} , \\ - c(T_{\xi} - u_{\xi}) &= q \omega, \\ - c \lambda_{\xi} &= \omega, \end{align*}

which is the limiting case of the system at \(\delta = 0\). The problem reduces to the nonlinear initial value problem

\begin{align*} \lambda_{\xi} = -\dfrac{k}{c} (1 - \lambda) \exp \left( \dfrac{-\theta}{T(\lambda)} \right), \end{align*}

with initial condition \(\lambda(0) = 0\); see [1] for details.

For the deflagration case, the initial guess is taken piecewise constant for \(u\) and \(T\), and

\begin{align*} \lambda(\xi) = \begin{cases} -\exp \left(\xi (1 - c) \Big/ \left(\dfrac{4 \delta}{3 \epsilon} \Pr \right) \right) + 1 \quad &\mathrm{for}\ \xi \in [l, 0], \\ 0 \quad &\mathrm{for}\ \xi \in (0, r] \end{cases} \end{align*}

for the reaction-progress variable. The value in the interval \((0, 1)\) is chosen as the initial guess for the front velocity \(c\).

Boundary conditions

In the numerical solution, the boundary conditions described in the beginning are imposed at the ends of the interval \(I\). In addition, a homogeneous Neumann condition is applied to the function \(d\lambda\) at the left boundary.

Program

Parameters

The calculation parameters are set in the ParametersList.prm file. To reproduce the results obtained below, you can run the program with the parameter files ParametersListDeflagrationSlow.prm, ParametersListDeflagrationFast.prm and ParametersListDetonation.prm.

Class <tt>TravelingWaveSolver</tt>

TravelingWaveSolver is the main class for computation of the traveling-wave profiles.

The implementation of Newton's method is based on that described in step-77 and relies on SUNDIALS' KINSOL package. Because of the additional unknown, the front velocity, we expand the Jacobi matrix by one column and one row (jacobian_matrix_extended), and add one more element to the solution vector (current_solution_extended). After completing the Newton iterations, we split the resulting extended solution vector current_solution_extended into two parts: the solution vector current_solution, corresponding to \((u, T, \lambda)\), and the front velocity current_wave_speed. After that the adaptive mesh refinement is performed using the current_solution vector, which is very important for resolving a narrow transition layer with a large solution gradient in the vicinity of zero. The KellyErrorEstimator is used as a refinement indicator.

Function <tt>calculate_profile</tt>

The full calculation cycle is done in the calculate_profile function. First, we construct an initial guess to the solution depending on the selected wave type and store the result as an object of type SolutionStruct. This object, along with the problem parameters, is then passed to the constructor of the TravelingWaveSolver class to calculate the traveling wave.

Decreasing the dissipation parameter \(\delta\) leads to the appearance of large gradients in solutions in the neighborhood of zero. As a consequence, Newton's method becomes more sensitive to the initial data and ceases to converge. To solve this problem, the calculate_profile function implements the method of continuation by the \(\delta\) parameter (for an example, see step-57). The solution and the refined triangulation are saved after each step of the method using the get_solution and get_triangulation functions and then passed to the next step.

Error estimation

Integration of the governing equations over the real line gives the following relations:

\begin{align*} u_l (1 - c) + \frac{\epsilon}{2} u_l^2 + \frac{\epsilon}{2} T_l &= u_r (1 - c) + \frac{\epsilon}{2} u_r^2 + \frac{\epsilon}{2} T_r , \\ T_l - u_l &= T_r - u_r + q. \end{align*}

These relations let us express any two parameters of \(c, T_l, T_r, u_l, u_r\) in terms of the remaining three. Thus, we can write

\begin{align*} u_l &= (T_l - T_r) + u_r - q, \\ c &= 1 + \epsilon \left( u_r - \dfrac{(q - (T_l - T_r))^2 + (T_l - T_r)}{2 (q - (T_l - T_r))} \right). \end{align*}

This means that since we choose the three parameters \(T_l, T_r, u_r\) for the detonation case ourselves, the above formulas give us the exact values of \(c\) and \(u_l\). These can be used to obtain the value of the error in the calculated \(c\) and \(u_l\).

For the deflagration case, however, we can only choose two parameters, \(T_l\) and \(u_r\). The remaining three are determined during the solution, so the formulas can only give us an error estimate.

Initial guess

To get the initial condition for detonation, we have to solve the nonlinear initial value problem for \(\lambda\) we mentioned earlier. This is done in the LimitSolution class. Numerical integration is performed using the odeint library of the Boost with its interface in IntegrateSystem.h.

Results

To visualize the computed profiles, one can use gnuplot typing

plot for [i=2:4] "solution_filename" using 1:i w p title word("u T lambda", i-1)

or execute the python script plot.py

python plot.py "solution_filename"

Slow deflagration for @f$\delta = 0.01@f$

The calculated wave speed is \(c = 0.0909\). Slow deflagration, @f$\delta = 0.01@f$, @f$c = 0.0909@f$

Fast deflagration for @f$\delta = 0.01@f$

The calculated wave speed is \(c = 0.8252\). Fast deflagration, @f$\delta = 0.01@f$, @f$0.8252@f$

Detonation for @f$\delta = 0.01@f$ and @f$\delta = 0.001@f$

The calculated wave speed in both cases is the same \(c = 1.216481\), as expected. Solid lines represent the detonation profile for the ideal case, when \(\delta=0\).

Detonation, @f$\delta = 0.01@f$, @f$c = 1.216481@f$.

Detonation, @f$\delta = 0.001@f$, @f$c = 1.216481@f$

Acknowledgments

I would like to thank my friend Oleg Rogozin for introducing me to the deal.II library and the world of finite elements.

References

  1. Goldin A.Y., Magomedov S.M., Faria L.M., Kasimov A.R. Study of a qualitative model for combustion waves: Flames, detonations, and deflagration-to-detonation transition. Computers & Fluids 2024; 273:106213.

Annotated version of AuxiliaryFunctions.h

Comparison of numbers with a given tolerance.

  template <typename T>
  bool isapprox(const T &a, const T &b, const double tol = 1e-10)
  {
  return (std::abs( a - b ) < tol);
  }
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)

Fill the std::vector with the values from the range [interval_begin, interval_end].

  template <typename T>
  void linspace(T interval_begin, T interval_end, std::vector<T> &arr)
  {
  const size_t SIZE = arr.size();
  const T step = (interval_end - interval_begin) / static_cast<T>(SIZE - 1);
  for (size_t i = 0; i < SIZE; ++i)
  {
  arr[i] = interval_begin + i * step;
  }
  }

Check the file existence.

  inline bool file_exists(const std::string &filename)
  {
  std::ifstream f(filename.c_str());
  return f.good();
  }

Annotated version of IntegrateSystem.h

  template <typename state_T, typename time_T>
  void SaveSolutionIntoFile(const std::vector<state_T>& x_vec, const std::vector<time_T>& t_vec, std::string filename="output_ode_sol.txt")
  {
  if (!x_vec.empty() && !t_vec.empty())
  {
  std::ofstream output(filename);
  output << std::setprecision(16);
  size_t dim = x_vec[0].size();
  for (size_t i = 0; i < t_vec.size(); ++i)
  {
  output << std::fixed << t_vec[i];
  for (size_t j = 0; j < dim; ++j)
  {
  output << std::scientific << " " << x_vec[i][j];
  }
  output << "\n";
  }
  output.close();
  }
  else
  {
  std::cout << "Solution is not saved into file.\n";
  }
  }

type of RK integrator

Observer

  template <typename state_type>
  {
  public:
  std::vector<state_type>& m_states;
  std::vector<double>& m_times;
  Push_back_state_time(std::vector<state_type>& states, std::vector<double>& times)
  {}
  void operator() (const state_type& x, double t)
  {
  m_states.push_back(x);
  m_times.push_back(t);
  }
  };
friend class Tensor
Definition tensor.h:882

Integrate system at specified points.

  template <typename ODE_obj_T, typename state_type, typename Iterable_type>
  std::vector<state_type>& x_vec, std::vector<double>& t_vec, const Iterable_type& iterable_time_span,
  state_type& x, const double dt,
  const double abs_er_tol=1.0e-15, const double rel_er_tol=1.0e-15
  )
  {
  using namespace boost::numeric::odeint;
  {
  }
  {
  }
  else
  { // Integrator_Type::fehlberg78
  }
  {
  }
  }
  #endif *

Annotated version of LimitSolution.cc

  #include "LimitSolution.h"
  namespace TravelingWave
  {
  LimitSolution::LimitSolution(const Parameters &parameters, const double ilambda_0, const double iu_0, const double iT_0, const double iroot_sign)
  : params(parameters)
  {
  }
  double LimitSolution::omega_func(const double lambda, const double T) const
  {
  return problem.k * (1. - lambda) * std::exp(-problem.theta / T);
  }
  void LimitSolution::operator() (const state_type &x , state_type &dxdt , const double /* t */)
  {
  dxdt[0] = -1. / wave_speed * omega_func(x[0], T_func(x[0]));
  }
  double LimitSolution::u_func(const double lambda) const
  {
  double coef = 2 * (wave_speed - 1) / problem.epsilon - 1;
  return (coef + root_sign * std::sqrt(coef * coef - 4 * (problem.q * lambda + B - 2 * A / problem.epsilon))) / 2;
  }
  double LimitSolution::T_func(const double lambda) const
  {
  return u_func(lambda) + problem.q * lambda + B;
  }
  {
  B = T_0 - u_0;
  A = u_0 * (1 - wave_speed) + problem.epsilon * (u_0 * u_0 + T_0) / 2;
  }
  {
  }
  {
  if (!t_vec.empty() && !lambda_vec.empty())
  {
  u_vec.resize(lambda_vec.size());
  T_vec.resize(lambda_vec.size());
  omega_vec.resize(lambda_vec.size());
  for (unsigned int i = 0; i < lambda_vec.size(); ++i)
  {
  u_vec[i].resize(1);
  T_vec[i].resize(1);
  omega_vec[i].resize(1);
  u_vec[i][0] = u_func(lambda_vec[i][0]);
  T_vec[i][0] = T_func(lambda_vec[i][0]);
  omega_vec[i][0] = omega_func(lambda_vec[i][0], T_vec[i][0]);
  }
  }
  else
  {
  std::cout << "t_vec or lambda_vec vector is empty!" << std::endl;
  }
  }
  } // namespace TravelingWave *
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)

Annotated version of LimitSolution.h

  #include "Parameters.h"
  namespace TravelingWave
  {
  typedef std::vector< double > state_type;
  {
  public:
  LimitSolution(const Parameters &parameters, const double ilambda_0, const double iu_0, const double iT_0, const double root_sign = 1.);
  void operator() (const state_type &x , state_type &dxdt , const double /* t */);
  std::vector<double> t_vec;
  std::vector<state_type> omega_vec;
  std::vector<state_type> lambda_vec;
  std::vector<state_type> u_vec;
  std::vector<state_type> T_vec;
  private:
  double omega_func(const double lambda, const double T) const;
  double u_func(const double lambda) const;
  double T_func(const double lambda) const;
  double wave_speed;
  const double lambda_0, u_0, T_0; // Initial values.
  double A, B; // Integration constants.
  const double root_sign; // Plus or minus one.
  };
  } // namespace TravelingWave
  #endif *

Annotated version of LinearInterpolator.h

Linear interpolation class

  template <typename Number_Type>
  {
  public:
  LinearInterpolator(const std::vector<Number_Type> &ix_points, const std::vector<Number_Type> &iy_points);
  Number_Type value(const Number_Type x) const;
  private:
  const std::vector<Number_Type> x_points; // Must be an increasing sequence, i.e. x[i] < x[i+1]
  const std::vector<Number_Type> y_points;
  };
  template <typename Number_Type>
  LinearInterpolator<Number_Type>::LinearInterpolator(const std::vector<Number_Type> &ix_points, const std::vector<Number_Type> &iy_points)
  {}
  {
  auto lower = std::lower_bound(x_points.begin(), x_points.end(), x);
  unsigned int right_index = 0;
  unsigned int left_index = 0;
  if (lower == x_points.begin())
  {
  }
  else if (lower == x_points.end())
  {
  res = y_points[x_points.size()-1];
  }
  else
  {
  res = (y_2 - y_1) / (x_2 - x_1) * (x - x_1) + y_1;
  }
  return res;
  }
  #endif *

Annotated version of Parameters.cc

  #include "Parameters.h"
  namespace TravelingWave
  {
  using namespace dealii;
  {
  add_parameter("delta", delta = 0.01);
  add_parameter("epsilon", epsilon = 0.1);
  add_parameter("Prandtl number", Pr = 0.75);
  add_parameter("Lewis number", Le = 1.0);
  add_parameter("Constant of reaction rate", k = 1.0);
  add_parameter("Activation energy", theta = 1.65);
  add_parameter("Heat release", q = 1.7);
  add_parameter("Ignition Temperature", T_ign = 1.0);
  add_parameter("Type of the wave (deflagration / detonation)", wave_type = 1); // 0 for "deflagration"; 1 for "detonation".
  add_parameter("Type of boundary condition for the temperature at the right boundary", T_r_bc_type = 1); // 0 for "Neumann" (deflagration); 1 for "Dirichlet" (detonation).
  add_parameter("T_left", T_left = 5.3); // Dirichlet boundary condition.
  add_parameter("T_right", T_right = 0.9); // For detonation waves the value serves as a Dirichlet boundary condition. For deflagration waves it serves for construction of the piecewise constant initial guess.
  add_parameter("u_left", u_left = -0.2); // For detonation waves the value is ignored. For deflagration waves it serves for construction of the piecewise constant initial guess.
  add_parameter("u_right", u_right = 0.); // Dirichlet boundary condition.
  add_parameter("Initial guess for the wave speed", wave_speed_init = 1.2); // For detonation waves the value is ignored. For deflagration waves it serves as an initial guess for the wave speed.
  }
  {
  add_parameter("Polynomial degree", poly_degree = 1);
  add_parameter("Number of quadrature points", quadrature_points_number = 3);
  }
  {
  add_parameter("Interval left boundary", interval_left = -50.0);
  add_parameter("Interval right boundary", interval_right = 20.0);
  add_parameter("Adaptive mesh refinement", adaptive = 1); // 1 for adaptive; 0 for global.
  }
  {
  add_parameter("Tolerance", tol = 1e-10);
  }
  } // namespace TravelingWave

Annotated version of Parameters.h

  namespace TravelingWave
  {
  using namespace dealii;
  {
  double delta, epsilon;
  double Pr, Le;
  double k, theta, q;
  double T_ign;
  double T_left, T_right;
  double u_left, u_right;
  };
  {
  unsigned int poly_degree;
  };
  {
  Mesh();
  unsigned int refinements_number;
  int adaptive;
  };
  {
  double tol;
  };
  struct Parameters
  {
  Solver solver;
  };
  } // namespace TravelingWave

Annotated version of Solution.cc

  #include "Solution.h"
  namespace TravelingWave
  {
  using namespace dealii;
  SolutionStruct::SolutionStruct(const std::vector<double> &ix, const std::vector<double> &iu,
  const std::vector<double> &iT, const std::vector<double> &ilambda, double iwave_speed)
  : x(ix)
  , u(iu)
  , T(iT)
  , lambda(ilambda)
  {}
  SolutionStruct::SolutionStruct(const std::vector<double> &ix, const std::vector<double> &iu,
  const std::vector<double> &iT, const std::vector<double> &ilambda)
  {}
  {
  T.clear();
  lambda.clear();
  lambda.resize(number_of_elements);
  }
  void SolutionStruct::save_to_file(std::string filename = "sol") const
  {
  const std::string file_for_solution = filename + ".txt";
  std::ofstream output(file_for_solution);
  output << std::scientific << std::setprecision(16);
  for (unsigned int i = 0; i < x.size(); ++i)
  {
  output << std::fixed << x[i];
  output << std::scientific << " " << u[i] << " " << T[i] << " " << lambda[i] << "\n";
  }
  output.close();
  std::ofstream file_for_wave_speed_output("wave_speed-" + file_for_solution);
  file_for_wave_speed_output << std::scientific << std::setprecision(16);
  }
  Interpolant::Interpolant(const std::vector<double> &ix_points, const std::vector<double> &iy_points)
  {}
  double Interpolant::value(const Point<1> &p, const unsigned int component) const
  {
  double x = p[0];
  double res = interpolant.value(x);
  return res;
  }
  template <typename InterpolantType>
  : Function<1>(3)
  {}
  {
  double res = 0.;
  if (component == 0) { res = u_interpolant.value(p); }
  else if (component == 1) { res = T_interpolant.value(p); }
  else if (component == 2) { res = lambda_interpolant.value(p); }
  return res;
  }
  } // namespace TravelingWave *
Definition point.h:111
constexpr void clear()
STL namespace.

Annotated version of Solution.h

The structure for keeping the solution: arrays of coordinates \(\xi\), solution \(u\), \(T\), \(\lambda\), and the wave speed \(c\).

  {
  SolutionStruct(const std::vector<double> &ix, const std::vector<double> &iu,
  const std::vector<double> &iT, const std::vector<double> &ilambda, const double iwave_speed);
  SolutionStruct(const std::vector<double> &ix, const std::vector<double> &iu,
  const std::vector<double> &iT, const std::vector<double> &ilambda);
  void reinit(const unsigned int number_of_elements);
  void save_to_file(std::string filename) const;
  std::vector<double> x; // mesh coordinates (must be an increasing sequence)
  std::vector<double> u; // array of u components
  std::vector<double> T; // array of T components
  std::vector<double> lambda; // array of lambda components
  double wave_speed; // speed of the wave
  };

Interpolation class

  class Interpolant : public Function<1>
  {
  public:
  Interpolant(const std::vector<double> &ix_points, const std::vector<double> &iy_points);
  virtual double value(const Point<1> &p, const unsigned int component = 0) const override;
  private:
  };

Vector function \((u(p), T(p), \lambda(p))\)

  template <typename InterpolantType>
  {
  public:
  virtual double value(const Point<1> &p, const unsigned int component = 0) const override;
  private:
  };
  } // namespace TravelingWave
  #endif *

Annotated version of TravelingWaveSolver.cc

  #include "TravelingWaveSolver.h"
  namespace TravelingWave
  {
  using namespace dealii;

Constructor of the class that takes parameters of the problem and an initial guess for Newton's iterations.

  : params(parameters)
  , fe(FE_Q<1>(params.fe.poly_degree), 1,
  FE_Q<1>(params.fe.poly_degree), 1) // 3 fe basis sets, corresponding to du, dT, dlambda
  , dof_handler(triangulation)
  {
Definition fe_q.h:554
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation

Table with values of some parameters to be written to the standard output before calculations.

  table.add_value("Parameter name", "number of quadrature points");
  table.add_value("value", number_of_quadrature_points);
  table.add_value("Parameter name", "delta");
  table.add_value("value", params.problem.delta);
  table.add_value("Parameter name", "epsilon");
  table.add_value("value", params.problem.epsilon);
  table.add_value("Parameter name", "params.problem.wave_speed_init");
  table.add_value("value", params.problem.wave_speed_init);
  table.add_value("Parameter name", "initial_guess.wave_speed");
  table.add_value("value", initial_guess.wave_speed);
  table.add_value("Parameter name", "T_left");
  table.add_value("value", params.problem.T_left);
  table.set_precision("value", 2);
  table.set_scientific("value", true);
  std::cout << "\n";
  std::cout << "\n";
  }

A function that takes a triangulation and assigns it to the member variable triangulation .

Here we find the indices of the degrees of freedom, associated with the boundary vertices, and the degree of freedom, associated with the vertex with coordinate \(\xi = 0\), and corresponding to temperature.

  {
  for (const auto &cell : dof_handler.active_cell_iterators())
  {
  {
  if (isapprox(cell->vertex(v_ind)[0], params.mesh.interval_left))
  {
  boundary_and_centering_dof_numbers["u_left"] = cell->vertex_dof_index(v_ind, 0);
  boundary_and_centering_dof_numbers["T_left"] = cell->vertex_dof_index(v_ind, 1);
  boundary_and_centering_dof_numbers["lambda_left"] = cell->vertex_dof_index(v_ind, 2);
  }
  else if (isapprox(cell->vertex(v_ind)[0], params.mesh.interval_right))
  {
  boundary_and_centering_dof_numbers["u_right"] = cell->vertex_dof_index(v_ind, 0);
  boundary_and_centering_dof_numbers["T_right"] = cell->vertex_dof_index(v_ind, 1);
  boundary_and_centering_dof_numbers["lambda_right"] = cell->vertex_dof_index(v_ind, 2);
  }
  else if (isapprox(cell->vertex(v_ind)[0], 0.))
  {
  boundary_and_centering_dof_numbers["T_zero"] = cell->vertex_dof_index(v_ind, 1);
  }
  }
  }
  }
unsigned int vertex_indices[2]

Set solution values, corresponding to Dirichlet boundary conditions and the centering condition \(T(0) = T_{\mathrm{ign}}\).

  {
  if (problem.T_r_bc_type == 1) // 1 for "Dirichlet"
  {
  } // else is 0 for "Neumann"
  }
  {
  dof_handler.distribute_dofs(fe);
  std::cout << "Number of dofs : " << dof_handler.n_dofs() << std::endl;
  extended_solution_dim = dof_handler.n_dofs() + 1;

Boundary condition constraints for \(du\), \(dT\) and \(d\lambda\).

Dirichlet homogeneous boundary condition for \(du\) at the right boundary.

Dirichlet homogeneous boundary condition for \(dT\) at the left boundary.

For the temperature at the left boundary there are two possibilities:

  if (problem.T_r_bc_type == 1) // 1 for "Dirichlet"
  {
  std::cout << "Dirichlet condition for the temperature at the right boundary." << std::endl;
  } // else is 0 for "Neumann"
  else
  {
  std::cout << "Neumann condition for the temperature at the right boundary." << std::endl;
  }

Dirichlet homogeneous boundary condition for \(d\lambda\) at the right boundary. (At the left boundary we consider the homogeneous Neumann boundary condition for \(d\lambda\).)

We create extended dynamic sparsity pattern with an additional row and an additional column.

  {
  std::vector<types::global_dof_index> dofs_on_this_cell;
  dofs_on_this_cell.reserve(dof_handler.get_fe_collection().max_dofs_per_cell());
  for (const auto &cell : dof_handler.active_cell_iterators())
  {
  const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
  dofs_on_this_cell.resize(dofs_per_cell);
  cell->get_dof_indices(dofs_on_this_cell);
  zero_boundary_constraints.add_entries_local_to_global(dofs_on_this_cell,
  /*keep_constrained_dofs*/ true);
  }

Adding elements to the last column.

  for (unsigned int i = 0; i < extended_solution_dim; ++i)
  {
  }

Adding one element to the last row, corresponding to the T(0).

Initialization

The initial condition is a discrete set of coordinates \(\xi\) and values of functions \(u\), \(T\) and \(\lambda\). From the three sets we create three continuous functions using interpolation, which then form one continuous vector function of SolutionVectorFunction type.

  for (unsigned int i = 0; i < extended_solution_dim - 1; ++i)
  {
  }
  }
void interpolate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const ComponentMask &component_mask={})

Heaviside function.

  {
  if (x > 0)
  {
  return 1.;
  }
  else
  {
  return 0.;
  }
  }
  {
  {
  TimerOutput::Scope t(computing_timer, "assembling the Jacobian");
  Vector<double> evaluation_point(dof_handler.n_dofs());
  for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
  {
  }
  std::cout << "Computing Jacobian matrix ... " << std::endl;
  FEValues<1> fe_values(fe,
  const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
  const unsigned int n_q_points = quadrature_formula.size();
  FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
  std::vector<double> current_velocity_values(n_q_points);
  std::vector<double> current_temperature_values(n_q_points);
  std::vector<double> current_lambda_values(n_q_points);
  std::vector<Tensor<1, 1>> current_velocity_gradients(n_q_points);
  std::vector<Tensor<1, 1>> current_temperature_gradients(n_q_points);
  std::vector<Tensor<1, 1>> current_lambda_gradients(n_q_points);
  std::vector<double> phi_u(dofs_per_cell);
  std::vector<Tensor<1, 1>> grad_phi_u(dofs_per_cell);
  std::vector<double> phi_T(dofs_per_cell);
  std::vector<Tensor<1, 1>> grad_phi_T(dofs_per_cell);
  std::vector<double> phi_lambda(dofs_per_cell);
  std::vector<Tensor<1, 1>> grad_phi_lambda(dofs_per_cell);
  for (const auto &cell : dof_handler.active_cell_iterators())
  {
  fe_values.reinit(cell);
  fe_values[velocity].get_function_values(evaluation_point, current_velocity_values);
  fe_values[temperature].get_function_values(evaluation_point, current_temperature_values);
  fe_values[lambda].get_function_values(evaluation_point, current_lambda_values);
  fe_values[velocity].get_function_gradients(evaluation_point, current_velocity_gradients);
  fe_values[temperature].get_function_gradients(evaluation_point, current_temperature_gradients);
  fe_values[lambda].get_function_gradients(evaluation_point, current_lambda_gradients);
  auto kappa_1 = [=](double T, double lambda){
  return problem.k * (1 - lambda) * std::exp(-problem.theta / T) * (
  problem.theta / (T * T) * Heaviside_func(T - problem.T_ign) /* + Delta_function(T - problem.T_ign) */
  );
  };
  auto kappa_2 = [=](double T, double lambda){
  return -problem.k * std::exp(-problem.theta / T) * Heaviside_func(T - problem.T_ign);
  };
  for (unsigned int q = 0; q < n_q_points; ++q)
  {
  for (unsigned int k = 0; k < dofs_per_cell; ++k)
  {
  phi_u[k] = fe_values[velocity].value(k, q);
  grad_phi_u[k] = fe_values[velocity].gradient(k, q);
  phi_T[k] = fe_values[temperature].value(k, q);
  grad_phi_T[k] = fe_values[temperature].gradient(k, q);
  phi_lambda[k] = fe_values[lambda].value(k, q);
  grad_phi_lambda[k] = fe_values[lambda].gradient(k, q);
  }
  const double del_Pr_eps = (problem.Pr * 4 * problem.delta / (3 * problem.epsilon));
  const double del_Le = (problem.delta / problem.Le);
  for (unsigned int i = 0; i < dofs_per_cell; ++i)
  {
  for (unsigned int j = 0; j < dofs_per_cell; ++j)
  {
  cell_matrix(i, j) += (
  + phi_u[i] * (
  - problem.epsilon / 2. * grad_phi_T[j][0]
  )
  + problem.delta * (-grad_phi_T[i] * grad_phi_T[j])
  + phi_T[i] * (
  )
  + phi_lambda[i] * (
  )
  ) * fe_values.JxW(q);
  }
  ) * fe_values.JxW(q);
  }
  }
  cell->get_dof_indices(local_dof_indices);
  for (const unsigned int i : fe_values.dof_indices())
  {
  for (const unsigned int j : fe_values.dof_indices())
  {
  jacobian_matrix_extended.add(local_dof_indices[i],
  local_dof_indices[j],
  }
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:74

Adding elements to the last column.

Global dof indices of dofs for \(dT\) and \(d\lambda\), associated with vertex \(\xi = 0\).

Approximating the derivative of \(T\) at \(\xi = 0\) as done in step-14.

  double T_at_zero_point(0.);
  {
  double derivative_evaluation_point = 0.; // Point at which T = T_ign.
  FEValues<1> fe_values(fe,
  const unsigned int n_q_points = quadrature_formula.size();
  std::vector<double> current_temperature_values(n_q_points);
  std::vector<Tensor<1, 1>> current_temperature_gradients(n_q_points);
  std::vector<double> current_lambda_values(n_q_points);
  for (const auto &cell : dof_handler.active_cell_iterators())
  {
  for (const auto &vertex : cell->vertex_indices())
  {
  if (isapprox(cell->vertex(vertex)[0], derivative_evaluation_point))
  {
  T_zero_point_dof_ind = cell->vertex_dof_index(vertex, 1);
  lambda_zero_point_dof_ind = cell->vertex_dof_index(vertex, 2);
  fe_values.reinit(cell);
  fe_values[temperature].get_function_values(current_solution, current_temperature_values);
  fe_values[temperature].get_function_gradients(current_solution, current_temperature_gradients);
  fe_values[lambda].get_function_values(current_solution, current_lambda_values);
  unsigned int q_point = 0;
  for (; q_point < n_q_points; ++q_point)
  {
  if (isapprox(fe_values.quadrature_point(q_point)[0], derivative_evaluation_point))
  {
  break;
  }
  }
  break;
  }
  }
  }
  }

Here we add to the matrix the terms that appear after integrating the terms with the Dirac delta function (which we skipped inside the loop).

Add 1 to the position T_zero_point_dof_ind of the last row of the matrix.

  }
  {
  TimerOutput::Scope t(computing_timer, "factorizing the Jacobian");
  std::cout << "Factorizing Jacobian matrix" << std::endl;
  jacobian_matrix_extended_factorization = std::make_unique<SparseDirectUMFPACK>();
  }
  }
  double TravelingWaveSolver::compute_residual(const Vector<double> &evaluation_point_extended, Vector<double> &residual)
  {
  TimerOutput::Scope t(computing_timer, "assembling the residual");
  std::cout << "Computing residual vector ... " << std::endl;
  residual = 0;
  Vector<double> evaluation_point(dof_handler.n_dofs());
  for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
  {
  }
  FEValues<1> fe_values(fe,
  const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
  const unsigned int n_q_points = quadrature_formula.size();
  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
  std::vector<double> current_velocity_values(n_q_points);
  std::vector<Tensor<1, 1>> current_velocity_gradients(n_q_points);
  std::vector<double> current_temperature_values(n_q_points);
  std::vector<Tensor<1, 1>> current_temperature_gradients(n_q_points);
  std::vector<double> current_lambda_values(n_q_points);
  std::vector<Tensor<1, 1>> current_lambda_gradients(n_q_points);
  std::vector<double> phi_u(dofs_per_cell);
  std::vector<Tensor<1, 1>> grad_phi_u(dofs_per_cell);
  std::vector<double> phi_T(dofs_per_cell);
  std::vector<Tensor<1, 1>> grad_phi_T(dofs_per_cell);
  std::vector<double> phi_lambda(dofs_per_cell);
  std::vector<Tensor<1, 1>> grad_phi_lambda(dofs_per_cell);
  for (const auto &cell : dof_handler.active_cell_iterators())
  {
  fe_values.reinit(cell);
  fe_values[velocity].get_function_values(evaluation_point, current_velocity_values);
  fe_values[velocity].get_function_gradients(evaluation_point, current_velocity_gradients);
  fe_values[temperature].get_function_values(evaluation_point, current_temperature_values);
  fe_values[temperature].get_function_gradients(evaluation_point, current_temperature_gradients);
  fe_values[lambda].get_function_values(evaluation_point, current_lambda_values);
  fe_values[lambda].get_function_gradients(evaluation_point, current_lambda_gradients);
  auto omega = [=](double T, double lambda){
  return problem.k * (1 - lambda) * std::exp(-problem.theta / T) * Heaviside_func(T - problem.T_ign);
  };
  for (unsigned int q = 0; q < n_q_points; ++q)
  {
  for (unsigned int k = 0; k < dofs_per_cell; ++k)
  {
  phi_u[k] = fe_values[velocity].value(k, q);
  grad_phi_u[k] = fe_values[velocity].gradient(k, q);
  phi_T[k] = fe_values[temperature].value(k, q);
  grad_phi_T[k] = fe_values[temperature].gradient(k, q);
  phi_lambda[k] = fe_values[lambda].value(k, q);
  grad_phi_lambda[k] = fe_values[lambda].gradient(k, q);
  }
  double del_Pr_eps = (problem.Pr * 4 * problem.delta / (3 * problem.epsilon));
  double del_Le = (problem.delta / problem.Le);
  for (unsigned int i = 0; i < dofs_per_cell; ++i)
  {
  + phi_u[i] * (
  - problem.epsilon / 2. * current_temperature_gradients[q][0]
  )
  + phi_T[i] * (
  )
  + phi_lambda[i] * (
  )
  ) * fe_values.JxW(q);
  }
  }
  cell->get_dof_indices(local_dof_indices);
  for (const unsigned int i : fe_values.dof_indices())
  {
  residual(local_dof_indices[i]) += cell_residual(i);
  }
  }
  residual(extended_solution_dim - 1) = 0.;
  zero_boundary_constraints.condense(residual);
  double residual_norm = residual.l2_norm();
  std::cout << std::defaultfloat;
  std::cout << "norm of residual = " << residual_norm << std::endl;
  }
void cell_residual(Vector< double > &result, const FEValuesBase< dim > &fe, const std::vector< Tensor< 1, dim > > &input, const ArrayView< const std::vector< double > > &velocity, double factor=1.)
Definition advection.h:130

Split the solution vector into two parts: one part is the solution \(u\), \(T\) and \(\lambda\), and another part is the wave speed.

  {
  for (unsigned int i = 0; i < extended_solution_dim - 1; ++i)
  {
  }
  }
  void TravelingWaveSolver::solve(const Vector<double> &rhs, Vector<double> &solution_extended, const double /*tolerance*/)
  {
  TimerOutput::Scope t(computing_timer, "linear system solve");
  std::cout << "Solving linear system ... " << std::endl;
  }

Function for adaptive mesh refinement based on KellyErrorEstimator .

  {
  dof_handler,
  QGauss<0>( 0 /* number_of_quadrature_points */),
  {},
  fe.component_mask(lambda)
  );
  0.1,
  0.05);
  triangulation.prepare_coarsening_and_refinement();
  solution_transfer.prepare_for_coarsening_and_refinement(current_solution);
  triangulation.execute_coarsening_and_refinement();
  setup_system(/*initial_step=*/ false);
  Vector<double> tmp(dof_handler.n_dofs());
  current_solution = std::move(tmp);
  for (unsigned int i = 0; i < extended_solution_dim - 1; ++i)
  {
  }
  }
  {
  double residual_norm = 0.;
  {
  typename SUNDIALS::KINSOL< Vector<double> >::AdditionalData additional_data;
  additional_data.function_tolerance = target_tolerance;
  nonlinear_solver.reinit_vector = [&](Vector<double> &x) {
  };
  residual_norm = compute_residual(evaluation_point, residual);
  return 0;
  };
  nonlinear_solver.setup_jacobian = [&](const Vector<double> &evaluation_point, const Vector<double> & /*current_f*/) {
  return 0;
  };
  nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &rhs, Vector<double> &solution, const double tolerance) {
  this->solve(rhs, solution, tolerance);
  return 0;
  };
  }
  }
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())

Output the solution ( \(u\), \(T\) and \(\lambda\)) and the wave speed into two separate files with double precision. The files can be read by gnuplot.

  void TravelingWaveSolver::output_with_double_precision(const Vector<double> &solution, const double wave_speed, const std::string filename)
  {
  TimerOutput::Scope t(computing_timer, "graphical output txt");
  const std::string file_for_solution = filename + ".txt";
  std::ofstream output(file_for_solution);
  for (const auto &cell : dof_handler.active_cell_iterators())
  {
  {
  double u = solution(cell->vertex_dof_index(v_ind, 0));
  double T = solution(cell->vertex_dof_index(v_ind, 1));
  double lambda = solution(cell->vertex_dof_index(v_ind, 2));
  output << std::scientific << std::setprecision(16);
  output << cell->vertex(v_ind)[0];
  output << std::scientific << std::setprecision(16);
  output << std::scientific << " " << u << " " << T << " " << lambda << "\n";
  }
  output << "\n";
  }
  output.close();
  std::ofstream file_for_wave_speed_output("wave_speed-" + file_for_solution);
  file_for_wave_speed_output << std::scientific << std::setprecision(16);
  }

Copy the solution into the SolutionStruct object, that stores the solution in an ordered manner.

To obtain an ordered solution array, we first create a set consisting of the elements {x, u, T, lambda} in which the sorting is done by coordinate, and then copy the contents of the set into the arrays of the SolutionStruct object.

  auto comp = [](const std::vector<double> &a, const std::vector<double> &b) {
  return a[0] < b[0];
  };
  std::set<std::vector<double>, decltype(comp)> solution_set(comp);
  for (const auto &cell : dof_handler.active_cell_iterators())
  {
  {
  double x = cell->vertex(v_ind)[0];
  double u = current_solution(cell->vertex_dof_index(v_ind, 0));
  double T = current_solution(cell->vertex_dof_index(v_ind, 1));
  double lambda = current_solution(cell->vertex_dof_index(v_ind, 2));
  solution_set.insert({x, u, T, lambda});
  }
  }
  solution.x.clear();
  solution.u.clear();
  solution.T.clear();
  solution.lambda.clear();
  solution.x.reserve(solution_set.size());
  solution.u.reserve(solution_set.size());
  solution.T.reserve(solution_set.size());
  solution.lambda.reserve(solution_set.size());
  for (auto it = solution_set.begin(); it != solution_set.end(); ++it)
  {
  solution.x.push_back((*it)[0]);
  solution.u.push_back((*it)[1]);
  solution.T.push_back((*it)[2]);
  solution.lambda.push_back((*it)[3]);
  }
  solution.wave_speed = current_wave_speed;
  }
  void TravelingWaveSolver::get_triangulation(Triangulation<1> &otriangulation) const
  {
  otriangulation.copy_triangulation(triangulation);
  }
  void TravelingWaveSolver::run(const std::string filename, const bool save_solution_to_file)
  {
  const int mesh_refinement_type = params.mesh.adaptive;
  const unsigned int n_refinements = params.mesh.refinements_number;
  const double tol = params.solver.tol;
  if (triangulation_uploaded == false) // If the triangulation is not loaded from outside, we will create one.
  {

We create two triangulations: one to the left and one to the right of zero coordinate. After that we merge them to obtain one triangulation, which contains zero point.

  static_cast<unsigned int>(std::abs( 0. - params.mesh.interval_left )),
  params.mesh.interval_left, 0.
  );
  static_cast<unsigned int>(std::abs( params.mesh.interval_right - 0. )),
  0., params.mesh.interval_right
  );
  }
  {
  if (mesh_refinement_type == 1) // For ADAPTIVE mesh refinement.
  {
  triangulation.refine_global(1); // refine initial mesh globally, before adaptive refinement cycles.
  }
  else if (mesh_refinement_type == 0) // For GLOBAL mesh refinement.
  {
  }
  }
  setup_system(/*initial step*/ true);
  {
  }
  if (mesh_refinement_type == 1) // Compute with ADAPTIVE mesh refinement.
  {
  double residual_norm = 0.;
  {
  }
  unsigned int refinement_cycle = 0;
  {
  std::cout << "Mesh refinement step " << refinement_cycle << std::endl;
  const double target_tolerance = 0.1 * std::pow(0.1, refinement_cycle); // Decrease tolerance for Newton solver at each refinement step.
  std::cout << " Target_tolerance: " << target_tolerance << std::endl;
  {
  std::cout << std::scientific << std::setprecision(16);
  std::cout << "current_wave_speed = " << current_wave_speed << std::endl;
  std::cout << std::defaultfloat;
  }
  computing_timer.print_summary();
  }
  {
  }
  }
  else if (mesh_refinement_type == 0) // Compute with GLOBAL mesh refinement.
  {
  {
  }
  {
  std::cout << std::scientific << std::setprecision(16);
  std::cout << "current_wave_speed = " << current_wave_speed << std::endl;
  std::cout << std::defaultfloat;
  }
  computing_timer.print_summary();
  }
  }
  } // namespace TravelingWave
void subdivided_hyper_cube(Triangulation< dim, spacedim > &tria, const unsigned int repetitions, const double left=0., const double right=1., const bool colorize=false)
void merge_triangulations(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result, const double duplicated_vertex_tolerance=1.0e-12, const bool copy_manifold_ids=false, const bool copy_boundary_ids=false)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)

Annotated version of TravelingWaveSolver.h

  #include <deal.II/base/timer.h>
  #include <deal.II/base/function.h>
  #include <deal.II/grid/tria.h>
  #include <deal.II/grid/grid_refinement.h>
  #include <deal.II/grid/grid_out.h>
  #include <deal.II/dofs/dof_handler.h>
  #include <deal.II/dofs/dof_tools.h>
  #include <deal.II/fe/fe_q.h>
  #include <deal.II/fe/fe_system.h>
  #include <deal.II/fe/fe_values.h>
  #include <deal.II/lac/vector.h>
  #include <deal.II/numerics/vector_tools.h>
  #include <deal.II/numerics/error_estimator.h>
  #include <deal.II/numerics/solution_transfer.h>
  #include <deal.II/numerics/data_out.h>
  #include <deal.II/sundials/kinsol.h>
  #include "Parameters.h"
  #include "Solution.h"
  #include "AuxiliaryFunctions.h"

Namespace of the program

  namespace TravelingWave
  {
  using namespace dealii;

The main class for construction of the traveling wave solutions.

  {
  public:
  void run(const std::string filename="solution", const bool save_solution_to_file=true);
  void get_solution(SolutionStruct &solution) const;
  void get_triangulation(Triangulation<1> &otriangulation) const;
  private:
  void setup_system(const bool initial_step);
  double Heaviside_func(double x) const;
  double compute_residual(const Vector<double> &evaluation_point_extended, Vector<double> &residual);
  void solve(const Vector<double> &rhs, Vector<double> &solution, const double /*tolerance*/);
  void refine_mesh();
  double run_newton_iterations(const double target_tolerance=1e-5);
  void output_with_double_precision(const Vector<double> &solution, const double wave_speed, const std::string filename="solution");

The dimension of the finite element solution increased by one to account for the value corresponding to the wave speed.

  std::map<std::string, unsigned int> boundary_and_centering_dof_numbers;

Parameters of the problem, taken from a .prm file.

  const Problem &problem; // Reference variable, just for convenience.

The flag indicating whether the triangulation was uploaded externally or created within the run member function.

Constraints for Dirichlet boundary conditions.

Finite element solution of the problem.

Value of the wave speed \(c\).

Solution with an additional term, corresponding to the variable wave_speed.

Initial guess for Newton's iterations.

  };
  } // namespace TravelingWave
  #endif *

Annotated version of calculate_profile.cc

  #include "TravelingWaveSolver.h"
  #include "calculate_profile.h"
  namespace TravelingWave
  {
  using namespace dealii;

Computation of the limit case (ideal) solution, corresponding to \(\delta = 0\), by solving the ODE. The output is the part of the solution to the left of zero. Here u_0, T_0, lambda_0 are the values of the medium state to the right of zero.

  void compute_limit_sol_left_part(const Parameters &parameters,
  const double wave_speed,
  const double u_0,
  const double T_0,
  const double lambda_0,
  const double root_sign)
  {
  limit_sol.set_wave_speed(wave_speed);
  {

We take more integration points to better resolve the transition layer.

  std::vector<double> t_span(static_cast<unsigned int>(std::abs( 0. - parameters.mesh.interval_left )));
  linspace(parameters.mesh.interval_left, finer_mesh_starting_value, t_span);
  std::vector<double> fine_grid(10000);
  t_span.insert(t_span.end(), fine_grid.begin(), fine_grid.end());

Reverse the order of the elements (because we need to perform back in time integration).

  std::reverse(t_span.begin(), t_span.end());
  lambda_val[0] = lambda_0; // initial value
  }
  limit_sol.calculate_u_T_omega();

Reverse the order of elements

  std::reverse(limit_sol.t_vec.begin(), limit_sol.t_vec.end());
  std::reverse(limit_sol.lambda_vec.begin(), limit_sol.lambda_vec.end());
  std::reverse(limit_sol.u_vec.begin(), limit_sol.u_vec.end());
  std::reverse(limit_sol.T_vec.begin(), limit_sol.T_vec.end());
  std::reverse(limit_sol.omega_vec.begin(), limit_sol.omega_vec.end());
  SaveSolutionIntoFile(limit_sol.lambda_vec, limit_sol.t_vec, "solution_lambda_limit.txt");
  SaveSolutionIntoFile(limit_sol.u_vec, limit_sol.t_vec, "solution_u_limit.txt");
  SaveSolutionIntoFile(limit_sol.T_vec, limit_sol.t_vec, "solution_T_limit.txt");
  SaveSolutionIntoFile(limit_sol.omega_vec, limit_sol.t_vec, "solution_omega_limit.txt");
  LimitSol.reinit(limit_sol.t_vec.size());
  LimitSol.wave_speed = wave_speed;
  for (unsigned int i=0; i < limit_sol.t_vec.size(); ++i)
  {
  LimitSol.x[i] = limit_sol.t_vec[i];
  LimitSol.u[i] = limit_sol.u_vec[i][0];
  LimitSol.T[i] = limit_sol.T_vec[i][0];
  LimitSol.lambda[i] = limit_sol.lambda_vec[i][0];
  }
  }

Construction of an initial guess for detonation wave solution. The ODE is solved for the ideal system with \(\delta = 0\).

  {
  const Problem &problem = params.problem;
  double current_wave_speed(problem.wave_speed_init);
  { // Here we compute the exact value of the wave speed c for the detonation case. We can do this because we have the Dirichlet boundary conditions T_l, T_r and u_r. Exact values of u_l and c are obtained using the integral relations.
  double DeltaT = problem.T_left - problem.T_right;
  double qDT = problem.q - DeltaT;
  current_wave_speed = 1. + problem.epsilon * (problem.u_right - (qDT * qDT + DeltaT) / (2 * qDT));
  }
  double u_0 = problem.u_right;
  double T_0 = problem.T_right;
  double lambda_0 = 0.;
  for (int i = initial_guess.x.size() - 1; i > - 1; --i)
  {
  if (isapprox(initial_guess.x[i], 0.))
  {
  initial_guess.u[i] = problem.u_right;
  initial_guess.T[i] = problem.T_ign;
  initial_guess.lambda[i] = 0.;
  break;
  }
  }

Adding the points to the right part of the interval (w.r.t. \(\xi = 0\)).

  for (unsigned int i = 0; i < number_of_additional_points; ++i)
  {
  initial_guess.x.push_back(params.mesh.interval_right / (std::pow(2., number_of_additional_points - 1 - i)));
  initial_guess.u.push_back(problem.u_right);
  initial_guess.T.push_back(problem.T_right);
  initial_guess.lambda.push_back(0.);
  }
  }

Construction of a piecewise constant initial guess for deflagration wave solution.

  {
  const Problem &problem = params.problem;
  double current_wave_speed(problem.wave_speed_init);
  double del_Pr_eps = (problem.Pr * 4 * problem.delta / (3 * problem.epsilon));
  double del_Le = (problem.delta / problem.Le);
  auto u_init_guess_func = [&](double x) {
  if (x < 0.)
  {
  return problem.u_left;
  }
  else
  {
  return problem.u_right;
  }
  };
  auto T_init_guess_func = [&](double x) {
  if (x < 0.)
  {
  return problem.T_left;
  }
  else if (isapprox(x, 0.))
  {
  return problem.T_ign;
  }
  else
  {
  return problem.T_right;
  }
  };
  auto lambda_init_guess_func = [=](double x) {
  if (x < 0.)
  {
  }
  else
  {
  return 0.;
  }
  };
  unsigned int number_of_points = multiplier_for_number_of_points * static_cast<unsigned int>(std::trunc(std::abs( params.mesh.interval_right - params.mesh.interval_left )));
  std::vector<double> x_span(number_of_points);
  linspace(params.mesh.interval_left, params.mesh.interval_right, x_span);
  std::vector<double> u_init_arr(number_of_points);
  std::vector<double> T_init_arr(number_of_points);
  std::vector<double> lambda_init_arr(number_of_points);
  for (unsigned int i = 0; i < number_of_points; ++i)
  {
  }
  }

Compute the traveling-wave profile. The continuation method can be switched on by setting the argument continuation_for_delta as true .

  void calculate_profile(Parameters& parameters,
  const bool continuation_for_delta /* Compute with the continuation. */,
  const double delta_start /* The starting value of delta for the continuation method. */,
  const unsigned int number_of_continuation_points)
  {
  if (parameters.problem.wave_type == 1) // detonation wave
  {
  }
  else if (parameters.problem.wave_type == 0) // deflagration wave
  {
  }
  {
  std::string filename = "solution_delta-" + Utilities::to_string(parameters.problem.delta) + "_eps-"
  + Utilities::to_string(parameters.problem.epsilon);
  wave.get_solution(sol);
  }
  else // Run with continuation_for_delta.
  {
  double delta_target = parameters.problem.delta;
  parameters.problem.delta = delta_start;
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:479

Generate a sequence of delta values being uniformly distributed in log10 scale.

  {
  double delta_start_log10 = std::log10(delta_start);
  double delta_target_log10 = std::log10(delta_target);
  std::vector<double> delta_log_span(delta_span.size());
  for (unsigned int i = 0; i < delta_span.size(); ++i)
  {
  }
  }
  bool first_iter_flag = true;
  for (double delta : delta_span)
  {
  parameters.problem.delta = delta;
  std::string filename = "solution_delta-" + Utilities::to_string(parameters.problem.delta) + "_eps-"
  + Utilities::to_string(parameters.problem.epsilon);
  {
  }
  else
  {
  wave.set_triangulation(refined_triangulation);
  }
  wave.get_solution(sol);
  wave.get_triangulation(refined_triangulation);
  }
  }

Error estimation.

  {
  unsigned int sol_length = sol.x.size();
  double u_r = sol.u[sol_length-1]; // Dirichlet boundary condition
  double T_r = sol.T[sol_length-1]; // Dirichlet condition only for detonation case
  double u_l = sol.u[0];
  double T_l = sol.T[0]; // Dirichlet boundary condition
  double wave_speed = sol.wave_speed;
  std::cout << "Error estimates:" << std::endl;
  double DeltaT = T_l - T_r;
  double qDT = parameters.problem.q - DeltaT;
  double wave_speed_formula = 1. + parameters.problem.epsilon * (u_r - (qDT * qDT + DeltaT) / (2 * qDT));
  std::cout << std::setw(18) << std::left << "For wave speed" << " : " << std::setw(5) << wave_speed - wave_speed_formula << std::endl;
  double u_l_formula = DeltaT + u_r - parameters.problem.q;
  std::cout << std::setw(18) << std::left << "For u_l" << " : " << std::setw(5) << u_l - u_l_formula << std::endl;
  }
  }
  } // namespace TravelingWave *

Annotated version of calculate_profile.h

  #include "Parameters.h"
  #include "Solution.h"
  #include "LimitSolution.h"
  #include "IntegrateSystem.h"
  #include "AuxiliaryFunctions.h"
  namespace TravelingWave
  {
  void compute_limit_sol_left_part(const Parameters &parameters,
  const double wave_speed,
  const double u_0,
  const double T_0,
  const double lambda_0,
  const double root_sign = 1.);
  void calculate_profile(Parameters& parameters,
  const bool continuation_for_delta=false /* Compute with the continuation. */,
  const double delta_start=0.01 /* The starting value of delta for the continuation method. */,
  const unsigned int number_of_continuation_points=10);
  } // namespace TravelingWave
  #endif *

Annotated version of main.cc

  #include "calculate_profile.h"
  int main(int argc, char *argv[])
  {
  try
  {
  using namespace TravelingWave;
  Parameters parameters;
  std::string prm_filename = "ParametersList.prm";
  if (argc > 1)
  {

Check if file argv[1] exists.

  if (file_exists(argv[1]))
  {
  }
  else
  {
  std::string errorMessage = "File \"" + std::string(argv[1]) + "\" is not found.";
  throw std::runtime_error(errorMessage);
  }
  }
  else
  {

Check if the file "ParametersList.prm" exists in the current or in the parent directory.

  {
  std::string errorMessage = "File \"" + prm_filename + "\" is not found.";
  throw std::runtime_error(errorMessage);
  }
  else
  {
  {
  }
  }
  }
  std::cout << "Reading parameters... " << std::flush;
  std::cout << "done" << std::endl;
  calculate_profile(parameters, /* With continuation_for_delta */ false, 0.1, 3);
  }
  catch (std::exception &exc)
  {
  std::cerr << std::endl
  << std::endl
  << "----------------------------------------------------"
  << std::endl;
  std::cerr << "Exception on processing: " << std::endl
  << exc.what() << std::endl
  << "Aborting!" << std::endl
  << "----------------------------------------------------"
  << std::endl;
  return 1;
  }
  catch (...)
  {
  std::cerr << std::endl
  << std::endl
  << "----------------------------------------------------"
  << std::endl;
  std::cerr << "Unknown exception!" << std::endl
  << "Aborting!" << std::endl
  << "----------------------------------------------------"
  << std::endl;
  return 1;
  }
  return 0;
  }
static void initialize(const std::string &filename="", const std::string &output_filename="", const ParameterHandler::OutputStyle output_style_for_output_filename=ParameterHandler::Short, ParameterHandler &prm=ParameterAcceptor::prm, const ParameterHandler::OutputStyle output_style_for_filename=ParameterHandler::DefaultStyle)

Annotated version of plot.py

import numpy as np
import matplotlib.pyplot as plt
import os
import sys
plot_params = {
#'backend': 'pdf',
# 'lines.marker' : 'x',
'scatter.marker' : 'x',
'lines.markersize' : 4,
'lines.linewidth' : 1,
'axes.labelsize': 16,
# 'textfontsize': 12,
'font.size' : 16,
'legend.fontsize': 16,
'xtick.labelsize': 14,
'ytick.labelsize': 14,
'text.usetex': True,
'figure.figsize': [9,6],
'axes.grid': True
}
plt.rcParams.update(plot_params)
if len(sys.argv) > 1:
filename = sys.argv[1]
if os.path.exists(filename):
data = np.loadtxt(filename, np.float64)
data_unique = np.unique(data, axis=0)
data_unique = np.array(sorted(data_unique, key=lambda x : x[0]))
x = data_unique[:, 0]
u_sol = data_unique[:, 1]
T_sol = data_unique[:, 2]
lambda_sol = data_unique[:, 3]
fig, ax = plt.subplots(nrows=1, ncols=1)
ax.scatter(x, u_sol, label=r"@f$u@f$", color='blue')
ax.scatter(x, T_sol, label=r"@f$T@f$", color='red')
ax.scatter(x, lambda_sol, label=r"@f$\lambda@f$", color='green')
# Plot of limit solutions for the detonation case. Uncomment, if needed.
#===============================================================#
'''
path_to_solution_files = os.path.split(filename)[0]
u_limit_path = os.path.join(path_to_solution_files, 'solution_u_limit.txt')
T_limit_path = os.path.join(path_to_solution_files, 'solution_T_limit.txt')
lambda_limit_path = os.path.join(path_to_solution_files, 'solution_lambda_limit.txt')
if os.path.exists(u_limit_path):
u_limit = np.loadtxt(u_limit_path, np.float64)
ax.plot(u_limit[:, 0], u_limit[:, 1], label=r"@f$u_{\mathrm{lim}}@f$", color='blue')
ax.plot([0, x[-1]], [u_sol[-1], u_sol[-1]], color='blue')
else:
print("No such file:", u_limit_path)
if os.path.exists(T_limit_path):
T_limit = np.loadtxt(T_limit_path, np.float64)
ax.plot(T_limit[:, 0], T_limit[:, 1], label=r"@f$T_{\mathrm{lim}}@f$", color='red')
ax.plot([0, x[-1]], [T_sol[-1], T_sol[-1]], color='red')
else:
print("No such file:", T_limit_path)
if os.path.exists(lambda_limit_path):
lambda_limit = np.loadtxt(lambda_limit_path, np.float64)
ax.plot(lambda_limit[:, 0], lambda_limit[:, 1], label=r"@f$\lambda_{\mathrm{lim}}@f$", color='green')
ax.plot([0, x[-1]], [lambda_sol[-1], lambda_sol[-1]], color='green')
else:
print("No such file:", lambda_limit_path)
'''
#===============================================================#
ax.set_xlabel(r"@f$\xi@f$")
ax.set_ylabel(r"@f$u, T, \lambda@f$")
ax.legend()
# plt.savefig("fast_deflagration_delta_0.01.png", bbox_inches='tight', dpi=500)
# plt.savefig('slow_deflagration_delta_0.01.png', bbox_inches='tight', dpi=500)
# plt.savefig('detonation_delta_0.01.png', bbox_inches='tight', dpi=500)
plt.show()
else:
print("No such file:", filename)