TravelingWaves.h

Go to the documentation of this file.

 
 *   
 *   #ifndef AUXILIARY_FUNCTIONS
 *   #define AUXILIARY_FUNCTIONS
 *   
 *   #include <vector>
 *   #include <cmath>
 *   #include <fstream>
 *   #include <string>
 *   #include <cstring>
 *   
 * @endcode
 * 
 * Comparison of numbers with a given tolerance.
 * 
 * @code
 *   template <typename T>
 *   bool isapprox(const T &a, const T &b, const double tol = 1e-10)
 *   {
 *     return (std::abs( a - b ) < tol);
 *   }
 *   
 * @endcode
 * 
 * Fill the std::vector with the values from the range [interval_begin, interval_end].
 * 
 * @code
 *   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;
 *     }
 *   }
 *   
 * @endcode
 * 
 * Check the file existence.
 * 
 * @code
 *   inline bool file_exists(const std::string &filename) 
 *   {
 *     std::ifstream f(filename.c_str());
 *     return f.good();
 *   }
 *   
 *   #endif
 * @endcode
 
 
<a name="ann-IntegrateSystem.h"></a>
<h1>Annotated version of IntegrateSystem.h</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Shamil Magomedov
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #ifndef INTEGRATE_SYSTEM
 *   #define INTEGRATE_SYSTEM
 *   
 *   #include <boost/numeric/odeint.hpp>
 *   
 *   #include <iostream>
 *   #include <fstream>
 *   #include <string>
 *   
 *   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";
 *     }
 *   }
 *   
 * @endcode
 * 
 * type of RK integrator 
 * 
 * @code
 *   enum class Integrator_Type 
 *   {
 *     dopri5,
 *     cash_karp54,
 *     fehlberg78
 *   };
 *   
 * @endcode
 * 
 * Observer
 * 
 * @code
 *   template <typename state_type>
 *   class Push_back_state_time
 *   {
 *   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)
 *       : m_states(states), m_times(times) 
 *     {}
 *   
 *     void operator() (const state_type& x, double t) 
 *     {
 *       m_states.push_back(x);
 *       m_times.push_back(t);
 *     }
 *   };
 *   
 *   
 * @endcode
 * 
 * Integrate system at specified points. 
 * 
 * @code
 *   template <typename ODE_obj_T, typename state_type, typename Iterable_type>
 *   void IntegrateSystemAtTimePoints(
 *     std::vector<state_type>& x_vec, std::vector<double>& t_vec, const Iterable_type& iterable_time_span,
 *     const ODE_obj_T& ode_system_obj, 
 *     state_type& x, const double dt,
 *     Integrator_Type integrator_type=Integrator_Type::dopri5, bool save_to_file_flag=false,
 *     const double abs_er_tol=1.0e-15, const double rel_er_tol=1.0e-15
 *     )
 *   {
 *     using namespace boost::numeric::odeint;
 *   
 *     if (integrator_type == Integrator_Type::dopri5) 
 *     {
 *       typedef runge_kutta_dopri5< state_type > error_stepper_type;
 *       integrate_times( make_controlled< error_stepper_type >(abs_er_tol, rel_er_tol),
 *                 ode_system_obj, x, iterable_time_span.begin(), iterable_time_span.end(), dt, Push_back_state_time< state_type >(x_vec, t_vec) );
 *     } 
 *     else if (integrator_type == Integrator_Type::cash_karp54) 
 *     {
 *       typedef runge_kutta_cash_karp54< state_type > error_stepper_type;
 *       integrate_times( make_controlled< error_stepper_type >(abs_er_tol, rel_er_tol),
 *                 ode_system_obj, x, iterable_time_span.begin(), iterable_time_span.end(), dt, Push_back_state_time< state_type >(x_vec, t_vec) );
 *     }
 *     else 
 *     { // Integrator_Type::fehlberg78
 *       typedef runge_kutta_fehlberg78< state_type > error_stepper_type;
 *       integrate_times( make_controlled< error_stepper_type >(abs_er_tol, rel_er_tol),
 *                 ode_system_obj, x, iterable_time_span.begin(), iterable_time_span.end(), dt, Push_back_state_time< state_type >(x_vec, t_vec) );
 *     }
 *   
 *     if (save_to_file_flag) 
 *     {
 *       SaveSolutionIntoFile(x_vec, t_vec);
 *     } 
 *   
 *   }
 *   
 *   #endif
 * @endcode
 
 
<a name="ann-LimitSolution.cc"></a>
<h1>Annotated version of LimitSolution.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Shamil Magomedov
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #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)
 *       , problem(params.problem)
 *       , wave_speed(problem.wave_speed_init)
 *       , lambda_0(ilambda_0)
 *       , u_0(iu_0)
 *       , T_0(iT_0)
 *       , root_sign(iroot_sign)
 *     {
 *       calculate_constants_A_B();
 *     }
 *   
 *     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;
 *     }
 *   
 *     void LimitSolution::calculate_constants_A_B()
 *     {
 *       B = T_0 - u_0;
 *       A = u_0 * (1 - wave_speed) + problem.epsilon * (u_0 * u_0  + T_0) / 2;
 *     }
 *   
 *     void LimitSolution::set_wave_speed(double iwave_speed)
 *     {
 *       wave_speed = iwave_speed;
 *       calculate_constants_A_B();
 *     }
 *   
 *     void LimitSolution::calculate_u_T_omega()
 *     {
 *       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
 * @endcode
 
 
<a name="ann-LimitSolution.h"></a>
<h1>Annotated version of LimitSolution.h</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Shamil Magomedov
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #ifndef LIMIT_SOLUTION
 *   #define LIMIT_SOLUTION
 *   
 *   #include "Parameters.h"
 *   #include <iostream>
 *   #include <vector>
 *   
 *   namespace TravelingWave
 *   {
 *     typedef std::vector< double > state_type;
 *   
 *     class LimitSolution
 *     {
 *     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 */);
 *       void calculate_u_T_omega();
 *       void set_wave_speed(double iwave_speed);
 *   
 *       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;
 *   
 *       void calculate_constants_A_B();
 *   
 *       const Parameters &params;
 *       const Problem    &problem;
 *       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
 * @endcode
 
 
<a name="ann-LinearInterpolator.h"></a>
<h1>Annotated version of LinearInterpolator.h</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Shamil Magomedov
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #ifndef LINEAR_INTERPOLATOR
 *   #define LINEAR_INTERPOLATOR
 *   
 *   #include <cmath>
 *   #include <algorithm>
 *   #include <vector>
 *   
 * @endcode
 * 
 * Linear interpolation class
 * 
 * @code
 *   template <typename Number_Type>
 *   class LinearInterpolator
 *   {
 *   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)
 *     : x_points(ix_points)
 *     , y_points(iy_points)
 *   {}
 *   
 *   template <typename Number_Type>
 *   Number_Type LinearInterpolator<Number_Type>::value(const Number_Type x) const
 *   {
 *     Number_Type res = 0.;
 *   
 *     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())
 *     {
 *       res = y_points[0];
 *     }
 *     else if (lower == x_points.end())
 *     {
 *       res = y_points[x_points.size()-1];
 *     }
 *     else
 *     {
 *       right_index = lower - x_points.begin();
 *       left_index = right_index - 1;
 *   
 *       Number_Type y_2 = y_points[right_index];
 *       Number_Type y_1 = y_points[left_index];
 *       Number_Type x_2 = x_points[right_index];
 *       Number_Type x_1 = x_points[left_index];
 *   
 *       res = (y_2 - y_1) / (x_2 - x_1) * (x - x_1) + y_1;
 *     }
 *   
 *     return res;
 *   }
 *   
 *   #endif
 * @endcode
 
 
<a name="ann-Parameters.cc"></a>
<h1>Annotated version of Parameters.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Shamil Magomedov
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "Parameters.h"
 *   
 *   namespace TravelingWave
 *   {
 *     using namespace dealii;
 *     
 *     Problem::Problem() 
 *       : ParameterAcceptor("Problem")
 *     {
 *       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.
 *     }
 *   
 *     FiniteElements::FiniteElements()
 *       : ParameterAcceptor("Finite elements")
 *     {
 *       add_parameter("Polynomial degree", poly_degree = 1);
 *       add_parameter("Number of quadrature points", quadrature_points_number = 3);
 *     }
 *   
 *     Mesh::Mesh()
 *       : ParameterAcceptor("Mesh")
 *     {
 *       add_parameter("Interval left boundary", interval_left = -50.0);  
 *       add_parameter("Interval right boundary", interval_right = 20.0);
 *       add_parameter<unsigned int>("Refinements number", refinements_number = 10);
 *       add_parameter("Adaptive mesh refinement", adaptive = 1);    // 1 for adaptive; 0 for global.
 *     }
 *   
 *     Solver::Solver()
 *       : ParameterAcceptor("Solver")
 *     {
 *       add_parameter("Tolerance", tol = 1e-10);
 *     }
 *   
 *   } // namespace TravelingWave
 * @endcode
 
 
<a name="ann-Parameters.h"></a>
<h1>Annotated version of Parameters.h</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Shamil Magomedov
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #ifndef PARAMETERS
 *   #define PARAMETERS
 *   
 *   #include <deal.II/base/parameter_acceptor.h>
 *   
 *   namespace TravelingWave
 *   {
 *     using namespace dealii;
 *   
 *     struct Problem : ParameterAcceptor
 *     {
 *       Problem();
 *   
 *       double delta, epsilon;
 *       double Pr, Le;
 *       double k, theta, q;
 *       double T_ign;
 *       int wave_type;
 *       int T_r_bc_type;
 *       double T_left, T_right;
 *       double u_left, u_right;
 *       
 *       double wave_speed_init;
 *     };
 *   
 *     struct FiniteElements : ParameterAcceptor
 *     {
 *       FiniteElements();
 *   
 *       unsigned int poly_degree;
 *       unsigned int quadrature_points_number;
 *     };
 *   
 *     struct Mesh : ParameterAcceptor
 *     {
 *       Mesh();
 *       
 *       double        interval_left;
 *       double        interval_right;
 *       unsigned int  refinements_number;
 *       int           adaptive;
 *     };
 *   
 *     struct Solver : ParameterAcceptor
 *     {
 *       Solver();
 *   
 *       double tol;
 *     };
 *   
 *     struct Parameters
 *     {
 *       Problem           problem;
 *       FiniteElements    fe;
 *       Mesh              mesh;
 *       Solver            solver;
 *     };
 *   
 *   } // namespace TravelingWave
 *   
 *   #endif
 * @endcode
 
 
<a name="ann-Solution.cc"></a>
<h1>Annotated version of Solution.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Shamil Magomedov
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "Solution.h"
 *   
 *   namespace TravelingWave
 *   {
 *   
 *     using namespace dealii;
 *     
 *     SolutionStruct::SolutionStruct() {}
 *     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)
 *       , wave_speed(iwave_speed)
 *     {}
 *     SolutionStruct::SolutionStruct(const std::vector<double> &ix, const std::vector<double> &iu, 
 *                                     const std::vector<double> &iT, const std::vector<double> &ilambda)
 *       : SolutionStruct(ix, iu, iT, ilambda, 0.)
 *     {}
 *   
 *     void SolutionStruct::reinit(const unsigned int number_of_elements)
 *     {
 *       wave_speed = 0.;
 *       x.clear();
 *       u.clear();
 *       T.clear();
 *       lambda.clear();
 *   
 *       x.resize(number_of_elements);
 *       u.resize(number_of_elements);
 *       T.resize(number_of_elements);
 *       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);
 *       file_for_wave_speed_output << wave_speed << std::endl;
 *       file_for_wave_speed_output.close();
 *     }
 *   
 *   
 *     Interpolant::Interpolant(const std::vector<double> &ix_points, const std::vector<double> &iy_points)
 *       : interpolant(ix_points, 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>
 *     SolutionVectorFunction<InterpolantType>::SolutionVectorFunction(InterpolantType iu_interpolant, InterpolantType iT_interpolant, InterpolantType ilambda_interpolant) 
 *       : Function<1>(3)
 *       , u_interpolant(iu_interpolant)
 *       , T_interpolant(iT_interpolant)
 *       , lambda_interpolant(ilambda_interpolant)
 *     {}
 *   
 *     template <typename InterpolantType>
 *     double SolutionVectorFunction<InterpolantType>::value(const Point<1> &p, const unsigned int component) const
 *     {
 *       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;
 *     }
 *   
 *     template class SolutionVectorFunction<Interpolant>;
 *   
 *   } // namespace TravelingWave
 * @endcode
 
 
<a name="ann-Solution.h"></a>
<h1>Annotated version of Solution.h</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Shamil Magomedov
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #ifndef SOLUTION
 *   #define SOLUTION
 *   
 *   #include <deal.II/base/function.h>
 *   
 *   #include "LinearInterpolator.h"
 *   
 *   #include <vector>
 *   #include <string>
 *   #include <fstream>
 *   #include <iostream>
 *   #include <iomanip>
 *   
 *   
 *   namespace TravelingWave
 *   {
 *     using namespace dealii;
 *   
 * @endcode
 * 
 * The structure for keeping the solution: arrays of coordinates @f$\xi@f$, solution @f$u@f$, @f$T@f$, @f$\lambda@f$, and the wave speed @f$c@f$. 
 * 
 * @code
 *     struct SolutionStruct
 *     {
 *       SolutionStruct();
 *       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
 *     };
 *   
 * @endcode
 * 
 * Interpolation class
 * 
 * @code
 *     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:
 *       LinearInterpolator<double> interpolant;
 *     };
 *   
 * @endcode
 * 
 * Vector function @f$(u(p), T(p), \lambda(p))@f$
 * 
 * @code
 *     template <typename InterpolantType>
 *     class SolutionVectorFunction : public Function<1>
 *     {
 *     public:
 *       SolutionVectorFunction(InterpolantType iu_interpolant, InterpolantType iT_interpolant, InterpolantType ilambda_interpolant);
 *       virtual double value(const Point<1> &p, const unsigned int component = 0) const override;
 *   
 *     private:
 *       InterpolantType u_interpolant;
 *       InterpolantType T_interpolant;
 *       InterpolantType lambda_interpolant;
 *     };
 *   
 *   } // namespace TravelingWave
 *   
 *   #endif
 * @endcode
 
 
<a name="ann-TravelingWaveSolver.cc"></a>
<h1>Annotated version of TravelingWaveSolver.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Shamil Magomedov
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "TravelingWaveSolver.h"
 *   
 *   namespace TravelingWave
 *   {
 *     using namespace dealii;
 *   
 * @endcode
 * 
 * Constructor of the class that takes parameters of the problem and an initial guess for Newton's iterations.
 * 
 * @code
 *     TravelingWaveSolver::TravelingWaveSolver(const Parameters &parameters, const SolutionStruct &initial_guess_input) 
 *       : params(parameters)
 *       , problem(params.problem)
 *       , number_of_quadrature_points((params.fe.quadrature_points_number > 0) ? params.fe.quadrature_points_number : (params.fe.poly_degree + 1))
 *       , triangulation_uploaded(false)
 *       , fe(FE_Q<1>(params.fe.poly_degree), 1, 
 *             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)
 *       , current_wave_speed(0.)
 *       , initial_guess(initial_guess_input)
 *       , computing_timer(std::cout, TimerOutput::never, TimerOutput::wall_times)
 *     {
 * @endcode
 * 
 * Table with values of some parameters to be written to the standard output before calculations.
 * 
 * @code
 *       TableHandler table;
 *       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";
 *       table.write_text(std::cout, TableHandler::TextOutputFormat::org_mode_table);
 *       std::cout << "\n";
 *     }
 *   
 * @endcode
 * 
 * A function that takes a triangulation and assigns it to the member variable <code>triangulation </code>.
 * 
 * @code
 *     void TravelingWaveSolver::set_triangulation(const Triangulation<1> &itriangulation)
 *     {
 *       triangulation.clear();
 *       triangulation.copy_triangulation(itriangulation);
 *       triangulation_uploaded = true;
 *     }
 *   
 * @endcode
 * 
 * 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 @f$\xi = 0@f$, and corresponding to temperature.
 * 
 * @code
 *     void TravelingWaveSolver::find_boundary_and_centering_dof_numbers()
 *     {
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *       {
 *         for (const auto &v_ind : cell->vertex_indices())
 *         {
 *           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);
 *           }
 *         }
 *       }
 *     }
 *   
 * @endcode
 * 
 * Set solution values, corresponding to Dirichlet boundary conditions and the centering condition @f$T(0) = T_{\mathrm{ign}}@f$.
 * 
 * @code
 *     void TravelingWaveSolver::set_boundary_and_centering_values()
 *     {
 *       current_solution[boundary_and_centering_dof_numbers["u_right"]] = problem.u_right;
 *   
 *       current_solution[boundary_and_centering_dof_numbers["T_left"]] = problem.T_left;
 *       if (problem.T_r_bc_type == 1)        // 1 for "Dirichlet"
 *       {
 *         current_solution[boundary_and_centering_dof_numbers["T_right"]] = problem.T_right;
 *       } // else is 0 for "Neumann"
 *       current_solution[boundary_and_centering_dof_numbers["T_zero"]] = problem.T_ign;
 *   
 *       current_solution[boundary_and_centering_dof_numbers["lambda_right"]] = 0.;
 *     }
 *   
 *   
 *     void TravelingWaveSolver::setup_system(const bool initial_step)
 *     {
 *       TimerOutput::Scope t(computing_timer, "set up");
 *   
 *       dof_handler.distribute_dofs(fe);
 *   
 *       std::cout << "Number of dofs : " << dof_handler.n_dofs() << std::endl;
 *   
 *       extended_solution_dim = dof_handler.n_dofs() + 1;
 *   
 *       find_boundary_and_centering_dof_numbers();
 *   
 * @endcode
 * 
 * Boundary condition constraints for @f$du@f$, @f$dT@f$ and @f$d\lambda@f$.
 * 
 * @code
 *       zero_boundary_constraints.clear();
 *   
 * @endcode
 * 
 * Dirichlet homogeneous boundary condition for @f$du@f$ at the right boundary.
 * 
 * @code
 *       zero_boundary_constraints.add_line(boundary_and_centering_dof_numbers["u_right"]);
 *   
 * @endcode
 * 
 * Dirichlet homogeneous boundary condition for @f$dT@f$ at the left boundary.
 * 
 * @code
 *       zero_boundary_constraints.add_line(boundary_and_centering_dof_numbers["T_left"]);
 * @endcode
 * 
 * For the temperature at the left boundary there are two possibilities:
 * 
 * @code
 *       if (problem.T_r_bc_type == 1)        // 1 for "Dirichlet"
 *       {
 *         std::cout << "Dirichlet condition for the temperature at the right boundary." << std::endl;
 *         zero_boundary_constraints.add_line(boundary_and_centering_dof_numbers["T_right"]);
 *       } // else is 0 for "Neumann"
 *       else
 *       {
 *         std::cout << "Neumann condition for the temperature at the right boundary." << std::endl;
 *       }
 *   
 * @endcode
 * 
 * Dirichlet homogeneous boundary condition for @f$d\lambda@f$ at the right boundary. (At the left boundary we consider the homogeneous Neumann boundary condition for @f$d\lambda@f$.)
 * 
 * @code
 *       zero_boundary_constraints.add_line(boundary_and_centering_dof_numbers["lambda_right"]);
 *       
 *       zero_boundary_constraints.close();
 *   
 * @endcode
 * 
 * We create extended dynamic sparsity pattern with an additional row and an additional column.
 * 
 * @code
 *       DynamicSparsityPattern dsp(extended_solution_dim);
 *       {
 *         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,
 *                                 dsp,
 *                                 /*keep_constrained_dofs*/ true);
 *         }
 *   
 * @endcode
 * 
 * Adding elements to the last column.
 * 
 * @code
 *         for (unsigned int i = 0; i < extended_solution_dim; ++i)
 *         {
 *           dsp.add(i, extended_solution_dim - 1);
 *         }
 * @endcode
 * 
 * Adding one element to the last row, corresponding to the T(0).
 * 
 * @code
 *         dsp.add(extended_solution_dim - 1, boundary_and_centering_dof_numbers["T_zero"]);
 *       }
 *       
 * @endcode
 * 
 * Initialization
 * 
 * @code
 *       sparsity_pattern_extended.copy_from(dsp);
 *       jacobian_matrix_extended.reinit(sparsity_pattern_extended); 
 *       jacobian_matrix_extended_factorization.reset();
 *   
 *       current_solution_extended.reinit(extended_solution_dim);
 *   
 *       if (initial_step)
 *       {
 *         current_solution.reinit(dof_handler.n_dofs());
 *       }
 *   
 *     }
 *   
 *   
 *     void TravelingWaveSolver::set_initial_guess()
 *     {
 *       current_wave_speed = initial_guess.wave_speed;
 *   
 * @endcode
 * 
 * The initial condition is a discrete set of coordinates @f$\xi@f$ and values of functions @f$u@f$, @f$T@f$ and @f$\lambda@f$. From the three sets we create three continuous functions using interpolation, which then form one continuous vector function of <code> SolutionVectorFunction </code> type.
 * 
 * @code
 *       Interpolant u_interpolant(initial_guess.x, initial_guess.u);
 *       Interpolant T_interpolant(initial_guess.x, initial_guess.T);
 *       Interpolant lambda_interpolant(initial_guess.x, initial_guess.lambda);
 *   
 *       SolutionVectorFunction init_guess_func(u_interpolant, T_interpolant, lambda_interpolant);
 *   
 *       VectorTools::interpolate(dof_handler, init_guess_func, current_solution); 
 *   
 *       set_boundary_and_centering_values();
 *   
 *       for (unsigned int i = 0; i < extended_solution_dim - 1; ++i)
 *       {
 *         current_solution_extended(i) = current_solution(i);
 *       }
 *       current_solution_extended(extended_solution_dim - 1) = current_wave_speed;
 *     }
 *   
 * @endcode
 * 
 * Heaviside function.
 * 
 * @code
 *     double TravelingWaveSolver::Heaviside_func(double x) const
 *     {
 *       if (x > 0) 
 *       {
 *         return 1.;
 *       }
 *       else
 *       {
 *         return 0.;
 *       }
 *     }
 *   
 *   
 *     void TravelingWaveSolver::compute_and_factorize_jacobian(const Vector<double> &evaluation_point_extended)
 *     {
 *       {
 *         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)
 *         {
 *           evaluation_point(i) = evaluation_point_extended(i);
 *         }
 *   
 *         const double wave_speed = evaluation_point_extended(extended_solution_dim - 1);
 *   
 *         std::cout << "Computing Jacobian matrix ... " << std::endl;
 *     
 *         const QGauss<1> quadrature_formula(number_of_quadrature_points);
 *     
 *         jacobian_matrix_extended = 0;
 *   
 *         FEValues<1> fe_values(fe,
 *                               quadrature_formula,
 *                               update_values | update_gradients | 
 *                               update_quadrature_points | update_JxW_values);
 *   
 *         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);
 *         Vector<double> row_last_element_vector(dofs_per_cell);
 *   
 *         std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *         
 *         const FEValuesExtractors::Scalar velocity(0);
 *         const FEValuesExtractors::Scalar temperature(1);
 *         const FEValuesExtractors::Scalar lambda(2);
 *   
 *         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())
 *         {
 *           cell_matrix = 0;
 *           row_last_element_vector = 0;
 *   
 *           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) += (
 *   
 *                   del_Pr_eps * (-grad_phi_u[i] * grad_phi_u[j])
 *                   + phi_u[i] * (
 *                     - (1 - wave_speed + problem.epsilon * current_velocity_values[q]) * grad_phi_u[j][0]
 *                     - problem.epsilon * current_velocity_gradients[q][0] * phi_u[j]
 *                     - problem.epsilon / 2. * grad_phi_T[j][0]
 *                   )
 *                   
 *                   + problem.delta * (-grad_phi_T[i] * grad_phi_T[j])
 *                   + phi_T[i] * (
 *                     - wave_speed * grad_phi_u[j][0]
 *                     + wave_speed * grad_phi_T[j][0]
 *                     + problem.q * kappa_1(current_temperature_values[q], current_lambda_values[q]) * phi_T[j]
 *                     + problem.q * kappa_2(current_temperature_values[q], current_lambda_values[q]) * phi_lambda[j]
 *                   )
 *   
 *                   + del_Le * (-grad_phi_lambda[i] * grad_phi_lambda[j])
 *                   + phi_lambda[i] * (
 *                     kappa_1(current_temperature_values[q], current_lambda_values[q]) * phi_T[j]
 *                     + wave_speed * grad_phi_lambda[j][0]
 *                     + kappa_2(current_temperature_values[q], current_lambda_values[q]) * phi_lambda[j]
 *                   )
 *   
 *                 ) * fe_values.JxW(q);
 *   
 *               }
 *   
 *               row_last_element_vector(i) += (
 *                 (phi_u[i] * current_velocity_gradients[q][0])
 *                 + (phi_T[i] * current_temperature_gradients[q][0])
 *                 - (phi_T[i] * current_velocity_gradients[q][0])
 *                 + (phi_lambda[i] * current_lambda_gradients[q][0])
 *               ) * 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],
 *                                             cell_matrix(i, j));
 *             }
 *             
 * @endcode
 * 
 * Adding elements to the last column.
 * 
 * @code
 *             jacobian_matrix_extended.add(local_dof_indices[i],
 *                                           extended_solution_dim - 1,
 *                                           row_last_element_vector(i));
 *           }
 *   
 *         }
 *   
 * @endcode
 * 
 * Global dof indices of dofs for @f$dT@f$ and @f$d\lambda@f$, associated with vertex @f$\xi = 0@f$.
 * 
 * @code
 *         types::global_dof_index T_zero_point_dof_ind(0), lambda_zero_point_dof_ind(0);        
 *         
 * @endcode
 * 
 * Approximating the derivative of @f$T@f$ at @f$\xi = 0@f$ as done in @ref step_14 "step-14".
 * 
 * @code
 *         double T_point_derivative(0.);
 *         double T_at_zero_point(0.);
 *         double lambda_at_zero_point(0.);
 *         {
 *           double derivative_evaluation_point = 0.;         // Point at which T = T_ign.
 *   
 *           const QTrapezoid<1> quadrature_formula;
 *           FEValues<1>        fe_values(fe,
 *                                 quadrature_formula,
 *                                 update_values | update_gradients | update_quadrature_points);
 *   
 *           const FEValuesExtractors::Scalar temperature(1);
 *           const FEValuesExtractors::Scalar lambda(2);
 *           
 *           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);
 *           
 *           unsigned int derivative_evaluation_point_hits = 0;
 *   
 *           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;
 *                   }
 *                 }
 *   
 *                 T_at_zero_point = current_temperature_values[q_point];
 *                 lambda_at_zero_point = current_lambda_values[q_point];
 *   
 *                 T_point_derivative += current_temperature_gradients[q_point][0];
 *                 ++derivative_evaluation_point_hits;
 *   
 *                 break;
 *               }
 *             }
 *           }
 *           T_point_derivative /= static_cast<double>(derivative_evaluation_point_hits);
 *         }
 *   
 * @endcode
 * 
 * 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).
 * 
 * @code
 *         double term_with_delta_func(0.);
 *         term_with_delta_func = problem.k * std::exp(-problem.theta / T_at_zero_point) * (1 - lambda_at_zero_point) / std::abs(T_point_derivative);
 *         jacobian_matrix_extended.add(T_zero_point_dof_ind, T_zero_point_dof_ind, problem.q * term_with_delta_func);
 *         jacobian_matrix_extended.add(lambda_zero_point_dof_ind, T_zero_point_dof_ind, term_with_delta_func);
 *   
 * @endcode
 * 
 * Add 1 to the position <code> T_zero_point_dof_ind </code> of the last row of the matrix.
 * 
 * @code
 *         jacobian_matrix_extended.add(extended_solution_dim - 1, T_zero_point_dof_ind, 1.);
 *         
 *         zero_boundary_constraints.condense(jacobian_matrix_extended);
 *       }
 *   
 *       {
 *         TimerOutput::Scope t(computing_timer, "factorizing the Jacobian");
 *   
 *         std::cout << "Factorizing Jacobian matrix" << std::endl;
 *   
 *         jacobian_matrix_extended_factorization = std::make_unique<SparseDirectUMFPACK>();
 *         jacobian_matrix_extended_factorization->factorize(jacobian_matrix_extended);
 *       }
 *   
 *     }
 *   
 *   
 *     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)
 *       {
 *         evaluation_point(i) = evaluation_point_extended(i);
 *       }
 *   
 *       const double wave_speed = evaluation_point_extended(extended_solution_dim - 1);
 *   
 *       const QGauss<1> quadrature_formula(number_of_quadrature_points);
 *       FEValues<1> fe_values(fe,
 *                   quadrature_formula,
 *                   update_values | update_gradients | 
 *                   update_quadrature_points | update_JxW_values);
 *    
 *       const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
 *       const unsigned int n_q_points    = quadrature_formula.size();
 *      
 *       Vector<double> cell_residual(dofs_per_cell); 
 *       std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *   
 *       const FEValuesExtractors::Scalar velocity(0);
 *       const FEValuesExtractors::Scalar temperature(1);
 *       const FEValuesExtractors::Scalar lambda(2);
 *   
 *       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())
 *       {
 *         cell_residual = 0;
 *   
 *         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)
 *           {
 *             cell_residual(i) += (
 *   
 *               del_Pr_eps * (-grad_phi_u[i] * current_velocity_gradients[q])
 *               + phi_u[i] * (
 *                 - current_velocity_gradients[q][0] * (1 - wave_speed + problem.epsilon * current_velocity_values[q]) 
 *                 - problem.epsilon / 2. * current_temperature_gradients[q][0]
 *               )
 *   
 *               + problem.delta * (-grad_phi_T[i] * current_temperature_gradients[q])
 *               + phi_T[i] * (
 *                 wave_speed * (current_temperature_gradients[q][0] - current_velocity_gradients[q][0])
 *                 + problem.q * omega(current_temperature_values[q], current_lambda_values[q])
 *               )
 *   
 *               + del_Le * (-grad_phi_lambda[i] * current_lambda_gradients[q])
 *               + phi_lambda[i] * (
 *                 wave_speed * current_lambda_gradients[q][0] + omega(current_temperature_values[q], current_lambda_values[q])
 *               )
 *   
 *             ) * 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;
 *   
 *       return residual_norm;
 *     }
 *   
 * @endcode
 * 
 * Split the solution vector into two parts: one part is the solution @f$u@f$, @f$T@f$ and @f$\lambda@f$, and another part is the wave speed.
 * 
 * @code
 *     void TravelingWaveSolver::split_extended_solution_vector()
 *     {
 *       for (unsigned int i = 0; i < extended_solution_dim - 1; ++i)
 *       {
 *         current_solution(i) = current_solution_extended(i);
 *       }
 *   
 *       current_wave_speed = current_solution_extended(extended_solution_dim - 1);
 *     }
 *   
 *   
 *     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;
 *   
 *       jacobian_matrix_extended_factorization->vmult(solution_extended, rhs);
 *    
 *       zero_boundary_constraints.distribute(solution_extended);
 *   
 *     }
 *   
 *   
 * @endcode
 * 
 * Function for adaptive mesh refinement based on <code> KellyErrorEstimator </code>.
 * 
 * @code
 *     void TravelingWaveSolver::refine_mesh()
 *     {
 *       Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
 *   
 *       const FEValuesExtractors::Scalar lambda(2);
 *   
 *       KellyErrorEstimator<1>::estimate(
 *         dof_handler,
 *         QGauss<0>( 0 /* number_of_quadrature_points */),
 *         {},
 *         current_solution,
 *         estimated_error_per_cell,
 *         fe.component_mask(lambda)
 *       );
 *   
 *       GridRefinement::refine_and_coarsen_fixed_number(triangulation,
 *                               estimated_error_per_cell,
 *                               0.1,
 *                               0.05);
 *   
 *       triangulation.prepare_coarsening_and_refinement();
 *   
 *       SolutionTransfer<1> solution_transfer(dof_handler);
 *       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());
 *       solution_transfer.interpolate(current_solution, tmp);
 *       current_solution = std::move(tmp);
 *   
 *       set_boundary_and_centering_values();
 *   
 *       for (unsigned int i = 0; i < extended_solution_dim - 1; ++i)
 *       {
 *         current_solution_extended(i) = current_solution(i);
 *       }
 *       current_solution_extended(extended_solution_dim - 1) = current_wave_speed;
 *   
 *     }
 *   
 *   
 *     double TravelingWaveSolver::run_newton_iterations(const double target_tolerance)
 *     {
 *       
 *       double residual_norm = 0.;
 *       {
 *         typename SUNDIALS::KINSOL< Vector<double> >::AdditionalData additional_data;
 *         additional_data.function_tolerance = target_tolerance;
 *     
 *         SUNDIALS::KINSOL<Vector<double>> nonlinear_solver(additional_data);
 *     
 *         nonlinear_solver.reinit_vector = [&](Vector<double> &x) {
 *           x.reinit(extended_solution_dim);
 *         };
 *   
 *         nonlinear_solver.residual = [&](const Vector<double> &evaluation_point, Vector<double> &residual) {
 *           residual_norm = compute_residual(evaluation_point, residual);
 *   
 *           return 0;
 *         };
 *     
 *         nonlinear_solver.setup_jacobian = [&](const Vector<double> &evaluation_point, const Vector<double> & /*current_f*/) {
 *           compute_and_factorize_jacobian(evaluation_point);
 *   
 *           return 0;
 *         };
 *   
 *         nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &rhs, Vector<double> &solution, const double tolerance) {
 *           this->solve(rhs, solution, tolerance);
 *   
 *           return 0;
 *         };
 *   
 *         nonlinear_solver.solve(current_solution_extended);
 *       }
 *   
 *       return residual_norm;
 *   
 *     }
 *   
 * @endcode
 * 
 * Output the solution (@f$u@f$, @f$T@f$ and @f$\lambda@f$) and the wave speed into two separate files with double precision. The files can be read by gnuplot.
 * 
 * @code
 *     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())
 *       {
 *         for (const auto &v_ind : cell->vertex_indices())
 *         {
 *           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);
 *       file_for_wave_speed_output << wave_speed << std::endl;
 *       file_for_wave_speed_output.close();
 *     }
 *   
 * @endcode
 * 
 * Copy the solution into the <code> SolutionStruct </code> object, that stores the solution in an ordered manner. 
 * 
 * @code
 *     void TravelingWaveSolver::get_solution(SolutionStruct &solution) const
 *     {
 * @endcode
 * 
 * To obtain an ordered solution array, we first create a set consisting of the elements <code> {x, u, T, lambda} </code> in which the sorting is done by coordinate, and then copy the contents of the set into the arrays of the <code> SolutionStruct </code> object.
 * 
 * @code
 *       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())
 *       {
 *         for (const auto &v_ind : cell->vertex_indices())
 *         {
 *           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.clear();
 *       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.
 *       {
 * @endcode
 * 
 * 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.
 * 
 * @code
 *         Triangulation<1> triangulation_left;
 *         GridGenerator::subdivided_hyper_cube(
 *           triangulation_left,
 *           static_cast<unsigned int>(std::abs( 0. - params.mesh.interval_left )),
 *           params.mesh.interval_left, 0.
 *         );
 *   
 *         Triangulation<1> triangulation_right;
 *         GridGenerator::subdivided_hyper_cube(
 *           triangulation_right,
 *           static_cast<unsigned int>(std::abs( params.mesh.interval_right - 0. )),
 *           0., params.mesh.interval_right
 *         );
 *   
 *         GridGenerator::merge_triangulations(triangulation_left, triangulation_right, triangulation);
 *   
 *       }
 *   
 *       if (triangulation_uploaded == false)
 *       {
 *         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.
 *         {
 *           triangulation.refine_global(n_refinements);
 *         }
 *       }
 *   
 *       setup_system(/*initial step*/ true);
 *       set_initial_guess();
 *   
 *       if (save_solution_to_file)
 *       {
 *         output_with_double_precision(current_solution, current_wave_speed, "solution_initial_data");
 *       }
 *   
 *       if (mesh_refinement_type == 1)    // Compute with ADAPTIVE mesh refinement.
 *       {
 *         double residual_norm = 0.;
 *         {
 *           Vector<double> tmp_residual(extended_solution_dim);
 *           residual_norm = compute_residual(current_solution_extended, tmp_residual);
 *         }
 *   
 *         unsigned int refinement_cycle = 0;
 *         while ((residual_norm > tol) && (refinement_cycle < n_refinements))
 *         {
 *           computing_timer.reset();
 *           std::cout << "Mesh refinement step " << refinement_cycle << std::endl;
 *   
 *           if (refinement_cycle != 0) { refine_mesh(); }
 *           
 *           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;
 *     
 *           residual_norm = run_newton_iterations(target_tolerance);
 *           split_extended_solution_vector();
 *   
 *           {
 *             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();
 *   
 *           ++refinement_cycle;
 *         }
 *         if (save_solution_to_file)
 *         {
 *           output_with_double_precision(current_solution, current_wave_speed, filename);
 *         }
 *   
 *       }
 *       else if (mesh_refinement_type == 0)   // Compute with GLOBAL mesh refinement.
 *       {
 *         run_newton_iterations(tol);
 *         split_extended_solution_vector();
 *         
 *         if (save_solution_to_file)
 *         {
 *           output_with_double_precision(current_solution, current_wave_speed, filename);
 *         }
 *         
 *         {
 *           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
 * @endcode
 
 
<a name="ann-TravelingWaveSolver.h"></a>
<h1>Annotated version of TravelingWaveSolver.h</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Shamil Magomedov
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #ifndef TRAVELING_WAVE_SOLVER
 *   #define TRAVELING_WAVE_SOLVER
 *   
 *   #include <deal.II/base/timer.h>
 *   #include <deal.II/base/function.h>
 *   #include <deal.II/base/table_handler.h>
 *   #include <deal.II/grid/tria.h>
 *   #include <deal.II/grid/grid_generator.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/lac/full_matrix.h>
 *   #include <deal.II/lac/sparse_matrix.h>
 *   #include <deal.II/lac/sparse_direct.h>
 *   #include <deal.II/lac/dynamic_sparsity_pattern.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"
 *   
 *   #include <cmath>
 *   #include <algorithm>
 *   #include <string>
 *   #include <fstream>
 *   #include <iostream>
 *   #include <iomanip>
 *   #include <set>
 *   
 * @endcode
 * 
 * Namespace of the program
 * 
 * @code
 *   namespace TravelingWave
 *   {
 *     using namespace dealii;
 *   
 * @endcode
 * 
 * The main class for construction of the traveling wave solutions.
 * 
 * @code
 *     class TravelingWaveSolver
 *     {
 *     public:
 *       TravelingWaveSolver(const Parameters &parameters, const SolutionStruct &initial_guess_input);
 *   
 *       void set_triangulation(const Triangulation<1> &itriangulation);
 *   
 *       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);
 *       void find_boundary_and_centering_dof_numbers();
 *       void set_boundary_and_centering_values();
 *   
 *       void set_initial_guess();
 *   
 *       double Heaviside_func(double x) const;
 *   
 *       void compute_and_factorize_jacobian(const Vector<double> &evaluation_point_extended);
 *       double compute_residual(const Vector<double> &evaluation_point_extended, Vector<double> &residual);
 *       void split_extended_solution_vector();
 *   
 *       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");
 *   
 * @endcode
 * 
 * The dimension of the finite element solution increased by one to account for the value corresponding to the wave speed.
 * 
 * @code
 *       unsigned int extended_solution_dim;
 *       std::map<std::string, unsigned int> boundary_and_centering_dof_numbers;
 *   
 * @endcode
 * 
 * Parameters of the problem, taken from a .prm file.
 * 
 * @code
 *       const Parameters  &params;
 *       const Problem     &problem;    // Reference variable, just for convenience.
 *   
 *       unsigned int number_of_quadrature_points;
 *   
 *       Triangulation<1> triangulation;
 * @endcode
 * 
 * The flag indicating whether the triangulation was uploaded externally or created within the <code> run </code> member function.
 * 
 * @code
 *       bool            triangulation_uploaded;
 *       FESystem<1>     fe;
 *       DoFHandler<1>   dof_handler;
 *   
 * @endcode
 * 
 * Constraints for Dirichlet boundary conditions.
 * 
 * @code
 *       AffineConstraints<double> zero_boundary_constraints;        
 *   
 *       SparsityPattern                       sparsity_pattern_extended;
 *       SparseMatrix<double>                  jacobian_matrix_extended;
 *       std::unique_ptr<SparseDirectUMFPACK>  jacobian_matrix_extended_factorization;
 *   
 * @endcode
 * 
 * Finite element solution of the problem.
 * 
 * @code
 *       Vector<double>  current_solution;
 *   
 * @endcode
 * 
 * Value of the wave speed @f$c@f$.
 * 
 * @code
 *       double          current_wave_speed;
 *   
 * @endcode
 * 
 * Solution with an additional term, corresponding to the variable wave_speed.
 * 
 * @code
 *       Vector<double>  current_solution_extended;
 *   
 * @endcode
 * 
 * Initial guess for Newton's iterations.
 * 
 * @code
 *       SolutionStruct  initial_guess;
 *   
 *       TimerOutput     computing_timer;
 *     };
 *   
 *   } // namespace TravelingWave
 *   
 *   #endif
 * @endcode
 
 
<a name="ann-calculate_profile.cc"></a>
<h1>Annotated version of calculate_profile.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Shamil Magomedov
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #include "TravelingWaveSolver.h"
 *   #include "calculate_profile.h"
 *   
 *   namespace TravelingWave
 *   {
 *     using namespace dealii;
 *   
 * @endcode
 * 
 * Computation of the limit case (ideal) solution, corresponding to @f$\delta = 0@f$, 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.
 * 
 * @code
 *     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, 
 *                                       SolutionStruct &LimitSol, 
 *                                       const double root_sign)
 *     {
 *       LimitSolution limit_sol(parameters, lambda_0, u_0, T_0, root_sign);
 *       limit_sol.set_wave_speed(wave_speed);
 *       
 *       {
 * @endcode
 * 
 * We take more integration points to better resolve the transition layer.
 * 
 * @code
 *         std::vector<double> t_span(static_cast<unsigned int>(std::abs( 0. - parameters.mesh.interval_left )));
 *         double finer_mesh_starting_value = -9.1;
 *         linspace(parameters.mesh.interval_left, finer_mesh_starting_value, t_span); 
 *         std::vector<double> fine_grid(10000);
 *         linspace(finer_mesh_starting_value + 1e-4, 0., fine_grid); 
 *         t_span.insert(t_span.end(), fine_grid.begin(), fine_grid.end());
 *   
 * @endcode
 * 
 * Reverse the order of the elements (because we need to perform back in time integration).
 * 
 * @code
 *         std::reverse(t_span.begin(), t_span.end());
 *   
 *         state_type lambda_val(1);
 *         lambda_val[0] = lambda_0;                 // initial value
 *         IntegrateSystemAtTimePoints(limit_sol.lambda_vec, limit_sol.t_vec, t_span,
 *           limit_sol, 
 *           lambda_val,
 *           -1e-2, Integrator_Type::dopri5);
 *       }
 *   
 *       limit_sol.calculate_u_T_omega();
 *   
 * @endcode
 * 
 * Reverse the order of elements
 * 
 * @code
 *       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];
 *       }
 *     }
 *   
 *   
 * @endcode
 * 
 * Construction of an initial guess for detonation wave solution. The ODE is solved for the ideal system with @f$\delta = 0@f$.
 * 
 * @code
 *     void compute_initial_guess_detonation(const Parameters &params, SolutionStruct &initial_guess, const double root_sign)
 *     {
 *       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.;
 *   
 *       compute_limit_sol_left_part(params, current_wave_speed, u_0, T_0, lambda_0, initial_guess, root_sign);
 *   
 *       initial_guess.wave_speed = current_wave_speed;
 *   
 *       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;
 *         }
 *       }
 *   
 * @endcode
 * 
 * Adding the points to the right part of the interval (w.r.t. @f$\xi = 0@f$).
 * 
 * @code
 *       unsigned int number_of_additional_points = 5;
 *       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.);
 *       }
 *   
 *     }
 *   
 *   
 * @endcode
 * 
 * Construction of a piecewise constant initial guess for deflagration wave solution.
 * 
 * @code
 *     void compute_initial_guess_deflagration(const Parameters &params, SolutionStruct &initial_guess)
 *     {
 *       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.)
 *         {
 *           return -std::exp(x * std::abs(1 - current_wave_speed) / del_Pr_eps) + 1;
 *         }
 *         else 
 *         {
 *           return 0.;
 *         }
 *       };
 *   
 *       unsigned int multiplier_for_number_of_points = 7;
 *       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)
 *       {
 *         u_init_arr[i] = u_init_guess_func(x_span[i]);
 *         T_init_arr[i] = T_init_guess_func(x_span[i]);
 *         lambda_init_arr[i] = lambda_init_guess_func(x_span[i]);
 *       }
 *   
 *       initial_guess.x = x_span;
 *       initial_guess.u = u_init_arr;
 *       initial_guess.T = T_init_arr;
 *       initial_guess.lambda = lambda_init_arr;
 *       initial_guess.wave_speed = current_wave_speed;
 *   
 *     }
 *   
 *   
 * @endcode
 * 
 * Compute the traveling-wave profile. The continuation method can be switched on by setting the argument <code> continuation_for_delta </code> as <code> true </code>.
 * 
 * @code
 *     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)
 *     {
 *       SolutionStruct sol;
 *   
 *       if (parameters.problem.wave_type == 1)                                 // detonation wave
 *       {
 *         compute_initial_guess_detonation(parameters, sol);
 *       }
 *       else if (parameters.problem.wave_type == 0)                // deflagration wave
 *       {
 *         compute_initial_guess_deflagration(parameters, sol);
 *       }
 *       
 *       if (continuation_for_delta == false)
 *       {
 *         TravelingWaveSolver wave(parameters, sol);
 *         std::string filename = "solution_delta-" + Utilities::to_string(parameters.problem.delta) + "_eps-" 
 *                                                           + Utilities::to_string(parameters.problem.epsilon);
 *         wave.run(filename);
 *         wave.get_solution(sol);
 *       }
 *       else        // Run with continuation_for_delta.
 *       {
 *         double delta_target = parameters.problem.delta;
 *         parameters.problem.delta = delta_start;
 *   
 *         std::vector<double> delta_span(number_of_continuation_points);
 *   
 * @endcode
 * 
 * Generate a sequence of delta values being uniformly distributed in log10 scale.
 * 
 * @code
 *         {
 *           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());
 *           linspace(delta_start_log10, delta_target_log10, delta_log_span);
 *   
 *           for (unsigned int i = 0; i < delta_span.size(); ++i)
 *           {
 *             delta_span[i] = std::pow(10, delta_log_span[i]);
 *           }
 *         }
 *   
 *         Triangulation<1> refined_triangulation;
 *         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);
 *   
 *           TravelingWaveSolver wave(parameters, sol);
 *   
 *           if (first_iter_flag)
 *           {
 *             first_iter_flag = false;
 *           }
 *           else
 *           {
 *             wave.set_triangulation(refined_triangulation);        
 *           }
 *   
 *           wave.run(filename);
 *           wave.get_solution(sol);
 *           wave.get_triangulation(refined_triangulation);
 *         }
 *   
 *       }
 *   
 * @endcode
 * 
 * Error estimation.
 * 
 * @code
 *       {
 *         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
 * @endcode
 
 
<a name="ann-calculate_profile.h"></a>
<h1>Annotated version of calculate_profile.h</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Shamil Magomedov
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #ifndef CALCULATE_PROFILE
 *   #define CALCULATE_PROFILE
 *   
 *   #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, 
 *                                       SolutionStruct &LimitSol, 
 *                                       const double root_sign = 1.);
 *   
 *           void compute_initial_guess_detonation(const Parameters &params, SolutionStruct &initial_guess, const double root_sign = 1.);
 *           void compute_initial_guess_deflagration(const Parameters &params, SolutionStruct &initial_guess);
 *   
 *     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
 * @endcode
 
 
<a name="ann-main.cc"></a>
<h1>Annotated version of main.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* -----------------------------------------------------------------------------
 *    *
 *    * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
 *    * Copyright (C) 2024 by Shamil Magomedov
 *    *
 *    * This file is part of the deal.II code gallery.
 *    *
 *    * -----------------------------------------------------------------------------
 *    */
 *   
 *   #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)
 *       {
 * @endcode
 * 
 * Check if file argv[1] exists.
 * 
 * @code
 *         if (file_exists(argv[1])) 
 *         {
 *           prm_filename = argv[1];
 *         }
 *         else
 *         {
 *           std::string errorMessage = "File \"" + std::string(argv[1]) + "\" is not found.";
 *           throw std::runtime_error(errorMessage);
 *         }
 *       }
 *       else
 *       {
 * @endcode
 * 
 * Check if the file "ParametersList.prm" exists in the current or in the parent directory.
 * 
 * @code
 *         if (!(file_exists(prm_filename) || file_exists("../" + prm_filename)))
 *         {
 *           std::string errorMessage = "File \"" + prm_filename + "\" is not found.";
 *           throw std::runtime_error(errorMessage);
 *         }
 *         else
 *         {
 *           if (!file_exists(prm_filename))
 *           {
 *             prm_filename = "../" + prm_filename;
 *           }
 *         }
 *       }
 *   
 *       std::cout << "Reading parameters... " << std::flush;
 *       ParameterAcceptor::initialize(prm_filename);
 *       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;
 *   }
 * @endcode
 
 
<a name="ann-plot.py"></a>
<h1>Annotated version of plot.py</h1>
@code{.py}
## -----------------------------------------------------------------------------
##
## SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
## Copyright (C) 2024 by Shamil Magomedov
##
## This file is part of the deal.II code gallery.
##
## -----------------------------------------------------------------------------
##
 
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)
@endcode
 
 
*/