Reference documentation for deal.II version 9.6.0
\(\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
TravelingWaves.h
Go to the documentation of this file.
1
374 *  
375 *   #ifndef AUXILIARY_FUNCTIONS
376 *   #define AUXILIARY_FUNCTIONS
377 *  
378 *   #include <vector>
379 *   #include <cmath>
380 *   #include <fstream>
381 *   #include <string>
382 *   #include <cstring>
383 *  
384 * @endcode
385 *
386 * Comparison of numbers with a given tolerance.
387 *
388 * @code
389 *   template <typename T>
390 *   bool isapprox(const T &a, const T &b, const double tol = 1e-10)
391 *   {
392 *   return (std::abs( a - b ) < tol);
393 *   }
394 *  
395 * @endcode
396 *
397 * Fill the std::vector with the values from the range [interval_begin, interval_end].
398 *
399 * @code
400 *   template <typename T>
401 *   void linspace(T interval_begin, T interval_end, std::vector<T> &arr)
402 *   {
403 *   const size_t SIZE = arr.size();
404 *   const T step = (interval_end - interval_begin) / static_cast<T>(SIZE - 1);
405 *   for (size_t i = 0; i < SIZE; ++i)
406 *   {
407 *   arr[i] = interval_begin + i * step;
408 *   }
409 *   }
410 *  
411 * @endcode
412 *
413 * Check the file existence.
414 *
415 * @code
416 *   inline bool file_exists(const std::string &filename)
417 *   {
418 *   std::ifstream f(filename.c_str());
419 *   return f.good();
420 *   }
421 *  
422 *   #endif
423 * @endcode
424
425
426<a name="ann-IntegrateSystem.h"></a>
427<h1>Annotated version of IntegrateSystem.h</h1>
428 *
429 *
430 *
431 *
432 * @code
433 *   /* -----------------------------------------------------------------------------
434 *   *
435 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
436 *   * Copyright (C) 2024 by Shamil Magomedov
437 *   *
438 *   * This file is part of the deal.II code gallery.
439 *   *
440 *   * -----------------------------------------------------------------------------
441 *   */
442 *  
443 *   #ifndef INTEGRATE_SYSTEM
444 *   #define INTEGRATE_SYSTEM
445 *  
446 *   #include <boost/numeric/odeint.hpp>
447 *  
448 *   #include <iostream>
449 *   #include <fstream>
450 *   #include <string>
451 *  
452 *   template <typename state_T, typename time_T>
453 *   void SaveSolutionIntoFile(const std::vector<state_T>& x_vec, const std::vector<time_T>& t_vec, std::string filename="output_ode_sol.txt")
454 *   {
455 *   if (!x_vec.empty() && !t_vec.empty())
456 *   {
457 *   std::ofstream output(filename);
458 *   output << std::setprecision(16);
459 *  
460 *   size_t dim = x_vec[0].size();
461 *   for (size_t i = 0; i < t_vec.size(); ++i)
462 *   {
463 *   output << std::fixed << t_vec[i];
464 *   for (size_t j = 0; j < dim; ++j)
465 *   {
466 *   output << std::scientific << " " << x_vec[i][j];
467 *   }
468 *   output << "\n";
469 *   }
470 *   output.close();
471 *   }
472 *   else
473 *   {
474 *   std::cout << "Solution is not saved into file.\n";
475 *   }
476 *   }
477 *  
478 * @endcode
479 *
480 * type of RK integrator
481 *
482 * @code
483 *   enum class Integrator_Type
484 *   {
485 *   dopri5,
486 *   cash_karp54,
487 *   fehlberg78
488 *   };
489 *  
490 * @endcode
491 *
492 * Observer
493 *
494 * @code
495 *   template <typename state_type>
496 *   class Push_back_state_time
497 *   {
498 *   public:
499 *   std::vector<state_type>& m_states;
500 *   std::vector<double>& m_times;
501 *  
502 *   Push_back_state_time(std::vector<state_type>& states, std::vector<double>& times)
503 *   : m_states(states), m_times(times)
504 *   {}
505 *  
506 *   void operator() (const state_type& x, double t)
507 *   {
508 *   m_states.push_back(x);
509 *   m_times.push_back(t);
510 *   }
511 *   };
512 *  
513 *  
514 * @endcode
515 *
516 * Integrate system at specified points.
517 *
518 * @code
519 *   template <typename ODE_obj_T, typename state_type, typename Iterable_type>
520 *   void IntegrateSystemAtTimePoints(
521 *   std::vector<state_type>& x_vec, std::vector<double>& t_vec, const Iterable_type& iterable_time_span,
522 *   const ODE_obj_T& ode_system_obj,
523 *   state_type& x, const double dt,
524 *   Integrator_Type integrator_type=Integrator_Type::dopri5, bool save_to_file_flag=false,
525 *   const double abs_er_tol=1.0e-15, const double rel_er_tol=1.0e-15
526 *   )
527 *   {
528 *   using namespace boost::numeric::odeint;
529 *  
530 *   if (integrator_type == Integrator_Type::dopri5)
531 *   {
532 *   typedef runge_kutta_dopri5< state_type > error_stepper_type;
533 *   integrate_times( make_controlled< error_stepper_type >(abs_er_tol, rel_er_tol),
534 *   ode_system_obj, x, iterable_time_span.begin(), iterable_time_span.end(), dt, Push_back_state_time< state_type >(x_vec, t_vec) );
535 *   }
536 *   else if (integrator_type == Integrator_Type::cash_karp54)
537 *   {
538 *   typedef runge_kutta_cash_karp54< state_type > error_stepper_type;
539 *   integrate_times( make_controlled< error_stepper_type >(abs_er_tol, rel_er_tol),
540 *   ode_system_obj, x, iterable_time_span.begin(), iterable_time_span.end(), dt, Push_back_state_time< state_type >(x_vec, t_vec) );
541 *   }
542 *   else
543 *   { // Integrator_Type::fehlberg78
544 *   typedef runge_kutta_fehlberg78< state_type > error_stepper_type;
545 *   integrate_times( make_controlled< error_stepper_type >(abs_er_tol, rel_er_tol),
546 *   ode_system_obj, x, iterable_time_span.begin(), iterable_time_span.end(), dt, Push_back_state_time< state_type >(x_vec, t_vec) );
547 *   }
548 *  
549 *   if (save_to_file_flag)
550 *   {
551 *   SaveSolutionIntoFile(x_vec, t_vec);
552 *   }
553 *  
554 *   }
555 *  
556 *   #endif
557 * @endcode
558
559
560<a name="ann-LimitSolution.cc"></a>
561<h1>Annotated version of LimitSolution.cc</h1>
562 *
563 *
564 *
565 *
566 * @code
567 *   /* -----------------------------------------------------------------------------
568 *   *
569 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
570 *   * Copyright (C) 2024 by Shamil Magomedov
571 *   *
572 *   * This file is part of the deal.II code gallery.
573 *   *
574 *   * -----------------------------------------------------------------------------
575 *   */
576 *  
577 *   #include "LimitSolution.h"
578 *  
579 *   namespace TravelingWave
580 *   {
581 *  
582 *   LimitSolution::LimitSolution(const Parameters &parameters, const double ilambda_0, const double iu_0, const double iT_0, const double iroot_sign)
583 *   : params(parameters)
584 *   , problem(params.problem)
585 *   , wave_speed(problem.wave_speed_init)
586 *   , lambda_0(ilambda_0)
587 *   , u_0(iu_0)
588 *   , T_0(iT_0)
589 *   , root_sign(iroot_sign)
590 *   {
591 *   calculate_constants_A_B();
592 *   }
593 *  
594 *   double LimitSolution::omega_func(const double lambda, const double T) const
595 *   {
596 *   return problem.k * (1. - lambda) * std::exp(-problem.theta / T);
597 *   }
598 *  
599 *   void LimitSolution::operator() (const state_type &x , state_type &dxdt , const double /* t */)
600 *   {
601 *   dxdt[0] = -1. / wave_speed * omega_func(x[0], T_func(x[0]));
602 *   }
603 *  
604 *   double LimitSolution::u_func(const double lambda) const
605 *   {
606 *   double coef = 2 * (wave_speed - 1) / problem.epsilon - 1;
607 *   return (coef + root_sign * std::sqrt(coef * coef - 4 * (problem.q * lambda + B - 2 * A / problem.epsilon))) / 2;
608 *   }
609 *  
610 *   double LimitSolution::T_func(const double lambda) const
611 *   {
612 *   return u_func(lambda) + problem.q * lambda + B;
613 *   }
614 *  
615 *   void LimitSolution::calculate_constants_A_B()
616 *   {
617 *   B = T_0 - u_0;
618 *   A = u_0 * (1 - wave_speed) + problem.epsilon * (u_0 * u_0 + T_0) / 2;
619 *   }
620 *  
621 *   void LimitSolution::set_wave_speed(double iwave_speed)
622 *   {
623 *   wave_speed = iwave_speed;
624 *   calculate_constants_A_B();
625 *   }
626 *  
627 *   void LimitSolution::calculate_u_T_omega()
628 *   {
629 *   if (!t_vec.empty() && !lambda_vec.empty())
630 *   {
631 *   u_vec.resize(lambda_vec.size());
632 *   T_vec.resize(lambda_vec.size());
633 *   omega_vec.resize(lambda_vec.size());
634 *   for (unsigned int i = 0; i < lambda_vec.size(); ++i)
635 *   {
636 *   u_vec[i].resize(1);
637 *   T_vec[i].resize(1);
638 *   omega_vec[i].resize(1);
639 *  
640 *   u_vec[i][0] = u_func(lambda_vec[i][0]);
641 *   T_vec[i][0] = T_func(lambda_vec[i][0]);
642 *   omega_vec[i][0] = omega_func(lambda_vec[i][0], T_vec[i][0]);
643 *   }
644 *   }
645 *   else
646 *   {
647 *   std::cout << "t_vec or lambda_vec vector is empty!" << std::endl;
648 *   }
649 *   }
650 *  
651 *   } // namespace TravelingWave
652 * @endcode
653
654
655<a name="ann-LimitSolution.h"></a>
656<h1>Annotated version of LimitSolution.h</h1>
657 *
658 *
659 *
660 *
661 * @code
662 *   /* -----------------------------------------------------------------------------
663 *   *
664 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
665 *   * Copyright (C) 2024 by Shamil Magomedov
666 *   *
667 *   * This file is part of the deal.II code gallery.
668 *   *
669 *   * -----------------------------------------------------------------------------
670 *   */
671 *  
672 *   #ifndef LIMIT_SOLUTION
673 *   #define LIMIT_SOLUTION
674 *  
675 *   #include "Parameters.h"
676 *   #include <iostream>
677 *   #include <vector>
678 *  
679 *   namespace TravelingWave
680 *   {
681 *   typedef std::vector< double > state_type;
682 *  
683 *   class LimitSolution
684 *   {
685 *   public:
686 *   LimitSolution(const Parameters &parameters, const double ilambda_0, const double iu_0, const double iT_0, const double root_sign = 1.);
687 *  
688 *   void operator() (const state_type &x , state_type &dxdt , const double /* t */);
689 *   void calculate_u_T_omega();
690 *   void set_wave_speed(double iwave_speed);
691 *  
692 *   std::vector<double> t_vec;
693 *   std::vector<state_type> omega_vec;
694 *   std::vector<state_type> lambda_vec;
695 *   std::vector<state_type> u_vec;
696 *   std::vector<state_type> T_vec;
697 *  
698 *   private:
699 *   double omega_func(const double lambda, const double T) const;
700 *   double u_func(const double lambda) const;
701 *   double T_func(const double lambda) const;
702 *  
703 *   void calculate_constants_A_B();
704 *  
705 *   const Parameters &params;
706 *   const Problem &problem;
707 *   double wave_speed;
708 *  
709 *   const double lambda_0, u_0, T_0; // Initial values.
710 *   double A, B; // Integration constants.
711 *  
712 *   const double root_sign; // Plus or minus one.
713 *   };
714 *  
715 *  
716 *   } // namespace TravelingWave
717 *  
718 *   #endif
719 * @endcode
720
721
722<a name="ann-LinearInterpolator.h"></a>
723<h1>Annotated version of LinearInterpolator.h</h1>
724 *
725 *
726 *
727 *
728 * @code
729 *   /* -----------------------------------------------------------------------------
730 *   *
731 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
732 *   * Copyright (C) 2024 by Shamil Magomedov
733 *   *
734 *   * This file is part of the deal.II code gallery.
735 *   *
736 *   * -----------------------------------------------------------------------------
737 *   */
738 *  
739 *   #ifndef LINEAR_INTERPOLATOR
740 *   #define LINEAR_INTERPOLATOR
741 *  
742 *   #include <cmath>
743 *   #include <algorithm>
744 *   #include <vector>
745 *  
746 * @endcode
747 *
748 * Linear interpolation class
749 *
750 * @code
751 *   template <typename Number_Type>
752 *   class LinearInterpolator
753 *   {
754 *   public:
755 *   LinearInterpolator(const std::vector<Number_Type> &ix_points, const std::vector<Number_Type> &iy_points);
756 *   Number_Type value(const Number_Type x) const;
757 *  
758 *   private:
759 *   const std::vector<Number_Type> x_points; // Must be an increasing sequence, i.e. x[i] < x[i+1]
760 *   const std::vector<Number_Type> y_points;
761 *   };
762 *  
763 *   template <typename Number_Type>
764 *   LinearInterpolator<Number_Type>::LinearInterpolator(const std::vector<Number_Type> &ix_points, const std::vector<Number_Type> &iy_points)
765 *   : x_points(ix_points)
766 *   , y_points(iy_points)
767 *   {}
768 *  
769 *   template <typename Number_Type>
770 *   Number_Type LinearInterpolator<Number_Type>::value(const Number_Type x) const
771 *   {
772 *   Number_Type res = 0.;
773 *  
774 *   auto lower = std::lower_bound(x_points.begin(), x_points.end(), x);
775 *   unsigned int right_index = 0;
776 *   unsigned int left_index = 0;
777 *   if (lower == x_points.begin())
778 *   {
779 *   res = y_points[0];
780 *   }
781 *   else if (lower == x_points.end())
782 *   {
783 *   res = y_points[x_points.size()-1];
784 *   }
785 *   else
786 *   {
787 *   right_index = lower - x_points.begin();
788 *   left_index = right_index - 1;
789 *  
790 *   Number_Type y_2 = y_points[right_index];
791 *   Number_Type y_1 = y_points[left_index];
792 *   Number_Type x_2 = x_points[right_index];
793 *   Number_Type x_1 = x_points[left_index];
794 *  
795 *   res = (y_2 - y_1) / (x_2 - x_1) * (x - x_1) + y_1;
796 *   }
797 *  
798 *   return res;
799 *   }
800 *  
801 *   #endif
802 * @endcode
803
804
805<a name="ann-Parameters.cc"></a>
806<h1>Annotated version of Parameters.cc</h1>
807 *
808 *
809 *
810 *
811 * @code
812 *   /* -----------------------------------------------------------------------------
813 *   *
814 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
815 *   * Copyright (C) 2024 by Shamil Magomedov
816 *   *
817 *   * This file is part of the deal.II code gallery.
818 *   *
819 *   * -----------------------------------------------------------------------------
820 *   */
821 *  
822 *   #include "Parameters.h"
823 *  
824 *   namespace TravelingWave
825 *   {
826 *   using namespace dealii;
827 *  
828 *   Problem::Problem()
829 *   : ParameterAcceptor("Problem")
830 *   {
831 *   add_parameter("delta", delta = 0.01);
832 *   add_parameter("epsilon", epsilon = 0.1);
833 *   add_parameter("Prandtl number", Pr = 0.75);
834 *   add_parameter("Lewis number", Le = 1.0);
835 *   add_parameter("Constant of reaction rate", k = 1.0);
836 *   add_parameter("Activation energy", theta = 1.65);
837 *   add_parameter("Heat release", q = 1.7);
838 *   add_parameter("Ignition Temperature", T_ign = 1.0);
839 *   add_parameter("Type of the wave (deflagration / detonation)", wave_type = 1); // 0 for "deflagration"; 1 for "detonation".
840 *  
841 *   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).
842 *  
843 *   add_parameter("T_left", T_left = 5.3); // Dirichlet boundary condition.
844 *   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.
845 *   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.
846 *   add_parameter("u_right", u_right = 0.); // Dirichlet boundary condition.
847 *  
848 *   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.
849 *   }
850 *  
851 *   FiniteElements::FiniteElements()
852 *   : ParameterAcceptor("Finite elements")
853 *   {
854 *   add_parameter("Polynomial degree", poly_degree = 1);
855 *   add_parameter("Number of quadrature points", quadrature_points_number = 3);
856 *   }
857 *  
858 *   Mesh::Mesh()
859 *   : ParameterAcceptor("Mesh")
860 *   {
861 *   add_parameter("Interval left boundary", interval_left = -50.0);
862 *   add_parameter("Interval right boundary", interval_right = 20.0);
863 *   add_parameter<unsigned int>("Refinements number", refinements_number = 10);
864 *   add_parameter("Adaptive mesh refinement", adaptive = 1); // 1 for adaptive; 0 for global.
865 *   }
866 *  
867 *   Solver::Solver()
868 *   : ParameterAcceptor("Solver")
869 *   {
870 *   add_parameter("Tolerance", tol = 1e-10);
871 *   }
872 *  
873 *   } // namespace TravelingWave
874 * @endcode
875
876
877<a name="ann-Parameters.h"></a>
878<h1>Annotated version of Parameters.h</h1>
879 *
880 *
881 *
882 *
883 * @code
884 *   /* -----------------------------------------------------------------------------
885 *   *
886 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
887 *   * Copyright (C) 2024 by Shamil Magomedov
888 *   *
889 *   * This file is part of the deal.II code gallery.
890 *   *
891 *   * -----------------------------------------------------------------------------
892 *   */
893 *  
894 *   #ifndef PARAMETERS
895 *   #define PARAMETERS
896 *  
897 *   #include <deal.II/base/parameter_acceptor.h>
898 *  
899 *   namespace TravelingWave
900 *   {
901 *   using namespace dealii;
902 *  
903 *   struct Problem : ParameterAcceptor
904 *   {
905 *   Problem();
906 *  
907 *   double delta, epsilon;
908 *   double Pr, Le;
909 *   double k, theta, q;
910 *   double T_ign;
911 *   int wave_type;
912 *   int T_r_bc_type;
913 *   double T_left, T_right;
914 *   double u_left, u_right;
915 *  
916 *   double wave_speed_init;
917 *   };
918 *  
919 *   struct FiniteElements : ParameterAcceptor
920 *   {
921 *   FiniteElements();
922 *  
923 *   unsigned int poly_degree;
924 *   unsigned int quadrature_points_number;
925 *   };
926 *  
927 *   struct Mesh : ParameterAcceptor
928 *   {
929 *   Mesh();
930 *  
931 *   double interval_left;
932 *   double interval_right;
933 *   unsigned int refinements_number;
934 *   int adaptive;
935 *   };
936 *  
937 *   struct Solver : ParameterAcceptor
938 *   {
939 *   Solver();
940 *  
941 *   double tol;
942 *   };
943 *  
944 *   struct Parameters
945 *   {
946 *   Problem problem;
947 *   FiniteElements fe;
948 *   Mesh mesh;
949 *   Solver solver;
950 *   };
951 *  
952 *   } // namespace TravelingWave
953 *  
954 *   #endif
955 * @endcode
956
957
958<a name="ann-Solution.cc"></a>
959<h1>Annotated version of Solution.cc</h1>
960 *
961 *
962 *
963 *
964 * @code
965 *   /* -----------------------------------------------------------------------------
966 *   *
967 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
968 *   * Copyright (C) 2024 by Shamil Magomedov
969 *   *
970 *   * This file is part of the deal.II code gallery.
971 *   *
972 *   * -----------------------------------------------------------------------------
973 *   */
974 *  
975 *   #include "Solution.h"
976 *  
977 *   namespace TravelingWave
978 *   {
979 *  
980 *   using namespace dealii;
981 *  
982 *   SolutionStruct::SolutionStruct() {}
983 *   SolutionStruct::SolutionStruct(const std::vector<double> &ix, const std::vector<double> &iu,
984 *   const std::vector<double> &iT, const std::vector<double> &ilambda, double iwave_speed)
985 *   : x(ix)
986 *   , u(iu)
987 *   , T(iT)
988 *   , lambda(ilambda)
989 *   , wave_speed(iwave_speed)
990 *   {}
991 *   SolutionStruct::SolutionStruct(const std::vector<double> &ix, const std::vector<double> &iu,
992 *   const std::vector<double> &iT, const std::vector<double> &ilambda)
993 *   : SolutionStruct(ix, iu, iT, ilambda, 0.)
994 *   {}
995 *  
996 *   void SolutionStruct::reinit(const unsigned int number_of_elements)
997 *   {
998 *   wave_speed = 0.;
999 *   x.clear();
1000 *   u.clear();
1001 *   T.clear();
1002 *   lambda.clear();
1003 *  
1004 *   x.resize(number_of_elements);
1005 *   u.resize(number_of_elements);
1006 *   T.resize(number_of_elements);
1007 *   lambda.resize(number_of_elements);
1008 *   }
1009 *  
1010 *   void SolutionStruct::save_to_file(std::string filename = "sol") const
1011 *   {
1012 *   const std::string file_for_solution = filename + ".txt";
1013 *   std::ofstream output(file_for_solution);
1014 *  
1015 *   output << std::scientific << std::setprecision(16);
1016 *   for (unsigned int i = 0; i < x.size(); ++i)
1017 *   {
1018 *   output << std::fixed << x[i];
1019 *   output << std::scientific << " " << u[i] << " " << T[i] << " " << lambda[i] << "\n";
1020 *   }
1021 *   output.close();
1022 *  
1023 *   std::ofstream file_for_wave_speed_output("wave_speed-" + file_for_solution);
1024 *   file_for_wave_speed_output << std::scientific << std::setprecision(16);
1025 *   file_for_wave_speed_output << wave_speed << std::endl;
1026 *   file_for_wave_speed_output.close();
1027 *   }
1028 *  
1029 *  
1030 *   Interpolant::Interpolant(const std::vector<double> &ix_points, const std::vector<double> &iy_points)
1031 *   : interpolant(ix_points, iy_points)
1032 *   {}
1033 *  
1034 *   double Interpolant::value(const Point<1> &p, const unsigned int component) const
1035 *   {
1036 *   double x = p[0];
1037 *   double res = interpolant.value(x);
1038 *  
1039 *   return res;
1040 *   }
1041 *  
1042 *  
1043 *   template <typename InterpolantType>
1044 *   SolutionVectorFunction<InterpolantType>::SolutionVectorFunction(InterpolantType iu_interpolant, InterpolantType iT_interpolant, InterpolantType ilambda_interpolant)
1045 *   : Function<1>(3)
1046 *   , u_interpolant(iu_interpolant)
1047 *   , T_interpolant(iT_interpolant)
1048 *   , lambda_interpolant(ilambda_interpolant)
1049 *   {}
1050 *  
1051 *   template <typename InterpolantType>
1052 *   double SolutionVectorFunction<InterpolantType>::value(const Point<1> &p, const unsigned int component) const
1053 *   {
1054 *   double res = 0.;
1055 *   if (component == 0) { res = u_interpolant.value(p); }
1056 *   else if (component == 1) { res = T_interpolant.value(p); }
1057 *   else if (component == 2) { res = lambda_interpolant.value(p); }
1058 *  
1059 *   return res;
1060 *   }
1061 *  
1062 *   template class SolutionVectorFunction<Interpolant>;
1063 *  
1064 *   } // namespace TravelingWave
1065 * @endcode
1066
1067
1068<a name="ann-Solution.h"></a>
1069<h1>Annotated version of Solution.h</h1>
1070 *
1071 *
1072 *
1073 *
1074 * @code
1075 *   /* -----------------------------------------------------------------------------
1076 *   *
1077 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1078 *   * Copyright (C) 2024 by Shamil Magomedov
1079 *   *
1080 *   * This file is part of the deal.II code gallery.
1081 *   *
1082 *   * -----------------------------------------------------------------------------
1083 *   */
1084 *  
1085 *   #ifndef SOLUTION
1086 *   #define SOLUTION
1087 *  
1088 *   #include <deal.II/base/function.h>
1089 *  
1090 *   #include "LinearInterpolator.h"
1091 *  
1092 *   #include <vector>
1093 *   #include <string>
1094 *   #include <fstream>
1095 *   #include <iostream>
1096 *   #include <iomanip>
1097 *  
1098 *  
1099 *   namespace TravelingWave
1100 *   {
1101 *   using namespace dealii;
1102 *  
1103 * @endcode
1104 *
1105 * 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$.
1106 *
1107 * @code
1108 *   struct SolutionStruct
1109 *   {
1110 *   SolutionStruct();
1111 *   SolutionStruct(const std::vector<double> &ix, const std::vector<double> &iu,
1112 *   const std::vector<double> &iT, const std::vector<double> &ilambda, const double iwave_speed);
1113 *   SolutionStruct(const std::vector<double> &ix, const std::vector<double> &iu,
1114 *   const std::vector<double> &iT, const std::vector<double> &ilambda);
1115 *  
1116 *   void reinit(const unsigned int number_of_elements);
1117 *  
1118 *   void save_to_file(std::string filename) const;
1119 *  
1120 *   std::vector<double> x; // mesh coordinates (must be an increasing sequence)
1121 *   std::vector<double> u; // array of u components
1122 *   std::vector<double> T; // array of T components
1123 *   std::vector<double> lambda; // array of lambda components
1124 *  
1125 *   double wave_speed; // speed of the wave
1126 *   };
1127 *  
1128 * @endcode
1129 *
1130 * Interpolation class
1131 *
1132 * @code
1133 *   class Interpolant : public Function<1>
1134 *   {
1135 *   public:
1136 *   Interpolant(const std::vector<double> &ix_points, const std::vector<double> &iy_points);
1137 *   virtual double value(const Point<1> &p, const unsigned int component = 0) const override;
1138 *  
1139 *   private:
1140 *   LinearInterpolator<double> interpolant;
1141 *   };
1142 *  
1143 * @endcode
1144 *
1145 * Vector function @f$(u(p), T(p), \lambda(p))@f$
1146 *
1147 * @code
1148 *   template <typename InterpolantType>
1149 *   class SolutionVectorFunction : public Function<1>
1150 *   {
1151 *   public:
1152 *   SolutionVectorFunction(InterpolantType iu_interpolant, InterpolantType iT_interpolant, InterpolantType ilambda_interpolant);
1153 *   virtual double value(const Point<1> &p, const unsigned int component = 0) const override;
1154 *  
1155 *   private:
1156 *   InterpolantType u_interpolant;
1157 *   InterpolantType T_interpolant;
1158 *   InterpolantType lambda_interpolant;
1159 *   };
1160 *  
1161 *   } // namespace TravelingWave
1162 *  
1163 *   #endif
1164 * @endcode
1165
1166
1167<a name="ann-TravelingWaveSolver.cc"></a>
1168<h1>Annotated version of TravelingWaveSolver.cc</h1>
1169 *
1170 *
1171 *
1172 *
1173 * @code
1174 *   /* -----------------------------------------------------------------------------
1175 *   *
1176 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
1177 *   * Copyright (C) 2024 by Shamil Magomedov
1178 *   *
1179 *   * This file is part of the deal.II code gallery.
1180 *   *
1181 *   * -----------------------------------------------------------------------------
1182 *   */
1183 *  
1184 *   #include "TravelingWaveSolver.h"
1185 *  
1186 *   namespace TravelingWave
1187 *   {
1188 *   using namespace dealii;
1189 *  
1190 * @endcode
1191 *
1192 * Constructor of the class that takes parameters of the problem and an initial guess for Newton's iterations.
1193 *
1194 * @code
1195 *   TravelingWaveSolver::TravelingWaveSolver(const Parameters &parameters, const SolutionStruct &initial_guess_input)
1196 *   : params(parameters)
1197 *   , problem(params.problem)
1198 *   , number_of_quadrature_points((params.fe.quadrature_points_number > 0) ? params.fe.quadrature_points_number : (params.fe.poly_degree + 1))
1199 *   , triangulation_uploaded(false)
1200 *   , fe(FE_Q<1>(params.fe.poly_degree), 1,
1201 *   FE_Q<1>(params.fe.poly_degree), 1,
1202 *   FE_Q<1>(params.fe.poly_degree), 1) // 3 fe basis sets, corresponding to du, dT, dlambda
1203 *   , dof_handler(triangulation)
1204 *   , current_wave_speed(0.)
1205 *   , initial_guess(initial_guess_input)
1206 *   , computing_timer(std::cout, TimerOutput::never, TimerOutput::wall_times)
1207 *   {
1208 * @endcode
1209 *
1210 * Table with values of some parameters to be written to the standard output before calculations.
1211 *
1212 * @code
1213 *   TableHandler table;
1214 *   table.add_value("Parameter name", "number of quadrature points");
1215 *   table.add_value("value", number_of_quadrature_points);
1216 *  
1217 *   table.add_value("Parameter name", "delta");
1218 *   table.add_value("value", params.problem.delta);
1219 *  
1220 *   table.add_value("Parameter name", "epsilon");
1221 *   table.add_value("value", params.problem.epsilon);
1222 *  
1223 *   table.add_value("Parameter name", "params.problem.wave_speed_init");
1224 *   table.add_value("value", params.problem.wave_speed_init);
1225 *  
1226 *   table.add_value("Parameter name", "initial_guess.wave_speed");
1227 *   table.add_value("value", initial_guess.wave_speed);
1228 *  
1229 *   table.add_value("Parameter name", "T_left");
1230 *   table.add_value("value", params.problem.T_left);
1231 *  
1232 *   table.set_precision("value", 2);
1233 *   table.set_scientific("value", true);
1234 *  
1235 *   std::cout << "\n";
1236 *   table.write_text(std::cout, TableHandler::TextOutputFormat::org_mode_table);
1237 *   std::cout << "\n";
1238 *   }
1239 *  
1240 * @endcode
1241 *
1242 * A function that takes a triangulation and assigns it to the member variable <code>triangulation </code>.
1243 *
1244 * @code
1245 *   void TravelingWaveSolver::set_triangulation(const Triangulation<1> &itriangulation)
1246 *   {
1247 *   triangulation.clear();
1248 *   triangulation.copy_triangulation(itriangulation);
1249 *   triangulation_uploaded = true;
1250 *   }
1251 *  
1252 * @endcode
1253 *
1254 * 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.
1255 *
1256 * @code
1257 *   void TravelingWaveSolver::find_boundary_and_centering_dof_numbers()
1258 *   {
1259 *   for (const auto &cell : dof_handler.active_cell_iterators())
1260 *   {
1261 *   for (const auto &v_ind : cell->vertex_indices())
1262 *   {
1263 *   if (isapprox(cell->vertex(v_ind)[0], params.mesh.interval_left))
1264 *   {
1265 *   boundary_and_centering_dof_numbers["u_left"] = cell->vertex_dof_index(v_ind, 0);
1266 *   boundary_and_centering_dof_numbers["T_left"] = cell->vertex_dof_index(v_ind, 1);
1267 *   boundary_and_centering_dof_numbers["lambda_left"] = cell->vertex_dof_index(v_ind, 2);
1268 *   }
1269 *   else if (isapprox(cell->vertex(v_ind)[0], params.mesh.interval_right))
1270 *   {
1271 *   boundary_and_centering_dof_numbers["u_right"] = cell->vertex_dof_index(v_ind, 0);
1272 *   boundary_and_centering_dof_numbers["T_right"] = cell->vertex_dof_index(v_ind, 1);
1273 *   boundary_and_centering_dof_numbers["lambda_right"] = cell->vertex_dof_index(v_ind, 2);
1274 *   }
1275 *   else if (isapprox(cell->vertex(v_ind)[0], 0.))
1276 *   {
1277 *   boundary_and_centering_dof_numbers["T_zero"] = cell->vertex_dof_index(v_ind, 1);
1278 *   }
1279 *   }
1280 *   }
1281 *   }
1282 *  
1283 * @endcode
1284 *
1285 * Set solution values, corresponding to Dirichlet boundary conditions and the centering condition @f$T(0) = T_{\mathrm{ign}}@f$.
1286 *
1287 * @code
1288 *   void TravelingWaveSolver::set_boundary_and_centering_values()
1289 *   {
1290 *   current_solution[boundary_and_centering_dof_numbers["u_right"]] = problem.u_right;
1291 *  
1292 *   current_solution[boundary_and_centering_dof_numbers["T_left"]] = problem.T_left;
1293 *   if (problem.T_r_bc_type == 1) // 1 for "Dirichlet"
1294 *   {
1295 *   current_solution[boundary_and_centering_dof_numbers["T_right"]] = problem.T_right;
1296 *   } // else is 0 for "Neumann"
1297 *   current_solution[boundary_and_centering_dof_numbers["T_zero"]] = problem.T_ign;
1298 *  
1299 *   current_solution[boundary_and_centering_dof_numbers["lambda_right"]] = 0.;
1300 *   }
1301 *  
1302 *  
1303 *   void TravelingWaveSolver::setup_system(const bool initial_step)
1304 *   {
1305 *   TimerOutput::Scope t(computing_timer, "set up");
1306 *  
1307 *   dof_handler.distribute_dofs(fe);
1308 *  
1309 *   std::cout << "Number of dofs : " << dof_handler.n_dofs() << std::endl;
1310 *  
1311 *   extended_solution_dim = dof_handler.n_dofs() + 1;
1312 *  
1313 *   find_boundary_and_centering_dof_numbers();
1314 *  
1315 * @endcode
1316 *
1317 * Boundary condition constraints for @f$du@f$, @f$dT@f$ and @f$d\lambda@f$.
1318 *
1319 * @code
1320 *   zero_boundary_constraints.clear();
1321 *  
1322 * @endcode
1323 *
1324 * Dirichlet homogeneous boundary condition for @f$du@f$ at the right boundary.
1325 *
1326 * @code
1327 *   zero_boundary_constraints.add_line(boundary_and_centering_dof_numbers["u_right"]);
1328 *  
1329 * @endcode
1330 *
1331 * Dirichlet homogeneous boundary condition for @f$dT@f$ at the left boundary.
1332 *
1333 * @code
1334 *   zero_boundary_constraints.add_line(boundary_and_centering_dof_numbers["T_left"]);
1335 * @endcode
1336 *
1337 * For the temperature at the left boundary there are two possibilities:
1338 *
1339 * @code
1340 *   if (problem.T_r_bc_type == 1) // 1 for "Dirichlet"
1341 *   {
1342 *   std::cout << "Dirichlet condition for the temperature at the right boundary." << std::endl;
1343 *   zero_boundary_constraints.add_line(boundary_and_centering_dof_numbers["T_right"]);
1344 *   } // else is 0 for "Neumann"
1345 *   else
1346 *   {
1347 *   std::cout << "Neumann condition for the temperature at the right boundary." << std::endl;
1348 *   }
1349 *  
1350 * @endcode
1351 *
1352 * 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$.)
1353 *
1354 * @code
1355 *   zero_boundary_constraints.add_line(boundary_and_centering_dof_numbers["lambda_right"]);
1356 *  
1357 *   zero_boundary_constraints.close();
1358 *  
1359 * @endcode
1360 *
1361 * We create extended dynamic sparsity pattern with an additional row and an additional column.
1362 *
1363 * @code
1364 *   DynamicSparsityPattern dsp(extended_solution_dim);
1365 *   {
1366 *   std::vector<types::global_dof_index> dofs_on_this_cell;
1367 *   dofs_on_this_cell.reserve(dof_handler.get_fe_collection().max_dofs_per_cell());
1368 *  
1369 *   for (const auto &cell : dof_handler.active_cell_iterators())
1370 *   {
1371 *   const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1372 *   dofs_on_this_cell.resize(dofs_per_cell);
1373 *   cell->get_dof_indices(dofs_on_this_cell);
1374 *  
1375 *   zero_boundary_constraints.add_entries_local_to_global(dofs_on_this_cell,
1376 *   dsp,
1377 *   /*keep_constrained_dofs*/ true);
1378 *   }
1379 *  
1380 * @endcode
1381 *
1382 * Adding elements to the last column.
1383 *
1384 * @code
1385 *   for (unsigned int i = 0; i < extended_solution_dim; ++i)
1386 *   {
1387 *   dsp.add(i, extended_solution_dim - 1);
1388 *   }
1389 * @endcode
1390 *
1391 * Adding one element to the last row, corresponding to the T(0).
1392 *
1393 * @code
1394 *   dsp.add(extended_solution_dim - 1, boundary_and_centering_dof_numbers["T_zero"]);
1395 *   }
1396 *  
1397 * @endcode
1398 *
1399 * Initialization
1400 *
1401 * @code
1402 *   sparsity_pattern_extended.copy_from(dsp);
1403 *   jacobian_matrix_extended.reinit(sparsity_pattern_extended);
1404 *   jacobian_matrix_extended_factorization.reset();
1405 *  
1406 *   current_solution_extended.reinit(extended_solution_dim);
1407 *  
1408 *   if (initial_step)
1409 *   {
1410 *   current_solution.reinit(dof_handler.n_dofs());
1411 *   }
1412 *  
1413 *   }
1414 *  
1415 *  
1416 *   void TravelingWaveSolver::set_initial_guess()
1417 *   {
1418 *   current_wave_speed = initial_guess.wave_speed;
1419 *  
1420 * @endcode
1421 *
1422 * 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.
1423 *
1424 * @code
1425 *   Interpolant u_interpolant(initial_guess.x, initial_guess.u);
1426 *   Interpolant T_interpolant(initial_guess.x, initial_guess.T);
1427 *   Interpolant lambda_interpolant(initial_guess.x, initial_guess.lambda);
1428 *  
1429 *   SolutionVectorFunction init_guess_func(u_interpolant, T_interpolant, lambda_interpolant);
1430 *  
1431 *   VectorTools::interpolate(dof_handler, init_guess_func, current_solution);
1432 *  
1433 *   set_boundary_and_centering_values();
1434 *  
1435 *   for (unsigned int i = 0; i < extended_solution_dim - 1; ++i)
1436 *   {
1437 *   current_solution_extended(i) = current_solution(i);
1438 *   }
1439 *   current_solution_extended(extended_solution_dim - 1) = current_wave_speed;
1440 *   }
1441 *  
1442 * @endcode
1443 *
1444 * Heaviside function.
1445 *
1446 * @code
1447 *   double TravelingWaveSolver::Heaviside_func(double x) const
1448 *   {
1449 *   if (x > 0)
1450 *   {
1451 *   return 1.;
1452 *   }
1453 *   else
1454 *   {
1455 *   return 0.;
1456 *   }
1457 *   }
1458 *  
1459 *  
1460 *   void TravelingWaveSolver::compute_and_factorize_jacobian(const Vector<double> &evaluation_point_extended)
1461 *   {
1462 *   {
1463 *   TimerOutput::Scope t(computing_timer, "assembling the Jacobian");
1464 *  
1465 *   Vector<double> evaluation_point(dof_handler.n_dofs());
1466 *   for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
1467 *   {
1468 *   evaluation_point(i) = evaluation_point_extended(i);
1469 *   }
1470 *  
1471 *   const double wave_speed = evaluation_point_extended(extended_solution_dim - 1);
1472 *  
1473 *   std::cout << "Computing Jacobian matrix ... " << std::endl;
1474 *  
1475 *   const QGauss<1> quadrature_formula(number_of_quadrature_points);
1476 *  
1477 *   jacobian_matrix_extended = 0;
1478 *  
1479 *   FEValues<1> fe_values(fe,
1480 *   quadrature_formula,
1483 *  
1484 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1485 *   const unsigned int n_q_points = quadrature_formula.size();
1486 *  
1487 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1488 *   Vector<double> row_last_element_vector(dofs_per_cell);
1489 *  
1490 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1491 *  
1492 *   const FEValuesExtractors::Scalar velocity(0);
1493 *   const FEValuesExtractors::Scalar temperature(1);
1494 *   const FEValuesExtractors::Scalar lambda(2);
1495 *  
1496 *   std::vector<double> current_velocity_values(n_q_points);
1497 *   std::vector<double> current_temperature_values(n_q_points);
1498 *   std::vector<double> current_lambda_values(n_q_points);
1499 *  
1500 *   std::vector<Tensor<1, 1>> current_velocity_gradients(n_q_points);
1501 *   std::vector<Tensor<1, 1>> current_temperature_gradients(n_q_points);
1502 *   std::vector<Tensor<1, 1>> current_lambda_gradients(n_q_points);
1503 *  
1504 *   std::vector<double> phi_u(dofs_per_cell);
1505 *   std::vector<Tensor<1, 1>> grad_phi_u(dofs_per_cell);
1506 *   std::vector<double> phi_T(dofs_per_cell);
1507 *   std::vector<Tensor<1, 1>> grad_phi_T(dofs_per_cell);
1508 *   std::vector<double> phi_lambda(dofs_per_cell);
1509 *   std::vector<Tensor<1, 1>> grad_phi_lambda(dofs_per_cell);
1510 *  
1511 *   for (const auto &cell : dof_handler.active_cell_iterators())
1512 *   {
1513 *   cell_matrix = 0;
1514 *   row_last_element_vector = 0;
1515 *  
1516 *   fe_values.reinit(cell);
1517 *  
1518 *   fe_values[velocity].get_function_values(evaluation_point, current_velocity_values);
1519 *   fe_values[temperature].get_function_values(evaluation_point, current_temperature_values);
1520 *   fe_values[lambda].get_function_values(evaluation_point, current_lambda_values);
1521 *  
1522 *   fe_values[velocity].get_function_gradients(evaluation_point, current_velocity_gradients);
1523 *   fe_values[temperature].get_function_gradients(evaluation_point, current_temperature_gradients);
1524 *   fe_values[lambda].get_function_gradients(evaluation_point, current_lambda_gradients);
1525 *  
1526 *   auto kappa_1 = [=](double T, double lambda){
1527 *   return problem.k * (1 - lambda) * std::exp(-problem.theta / T) * (
1528 *   problem.theta / (T * T) * Heaviside_func(T - problem.T_ign) /* + Delta_function(T - problem.T_ign) */
1529 *   );
1530 *   };
1531 *  
1532 *   auto kappa_2 = [=](double T, double lambda){
1533 *   return -problem.k * std::exp(-problem.theta / T) * Heaviside_func(T - problem.T_ign);
1534 *   };
1535 *  
1536 *   for (unsigned int q = 0; q < n_q_points; ++q)
1537 *   {
1538 *   for (unsigned int k = 0; k < dofs_per_cell; ++k)
1539 *   {
1540 *   phi_u[k] = fe_values[velocity].value(k, q);
1541 *   grad_phi_u[k] = fe_values[velocity].gradient(k, q);
1542 *   phi_T[k] = fe_values[temperature].value(k, q);
1543 *   grad_phi_T[k] = fe_values[temperature].gradient(k, q);
1544 *   phi_lambda[k] = fe_values[lambda].value(k, q);
1545 *   grad_phi_lambda[k] = fe_values[lambda].gradient(k, q);
1546 *   }
1547 *  
1548 *   const double del_Pr_eps = (problem.Pr * 4 * problem.delta / (3 * problem.epsilon));
1549 *   const double del_Le = (problem.delta / problem.Le);
1550 *  
1551 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1552 *   {
1553 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1554 *   {
1555 *   cell_matrix(i, j) += (
1556 *  
1557 *   del_Pr_eps * (-grad_phi_u[i] * grad_phi_u[j])
1558 *   + phi_u[i] * (
1559 *   - (1 - wave_speed + problem.epsilon * current_velocity_values[q]) * grad_phi_u[j][0]
1560 *   - problem.epsilon * current_velocity_gradients[q][0] * phi_u[j]
1561 *   - problem.epsilon / 2. * grad_phi_T[j][0]
1562 *   )
1563 *  
1564 *   + problem.delta * (-grad_phi_T[i] * grad_phi_T[j])
1565 *   + phi_T[i] * (
1566 *   - wave_speed * grad_phi_u[j][0]
1567 *   + wave_speed * grad_phi_T[j][0]
1568 *   + problem.q * kappa_1(current_temperature_values[q], current_lambda_values[q]) * phi_T[j]
1569 *   + problem.q * kappa_2(current_temperature_values[q], current_lambda_values[q]) * phi_lambda[j]
1570 *   )
1571 *  
1572 *   + del_Le * (-grad_phi_lambda[i] * grad_phi_lambda[j])
1573 *   + phi_lambda[i] * (
1574 *   kappa_1(current_temperature_values[q], current_lambda_values[q]) * phi_T[j]
1575 *   + wave_speed * grad_phi_lambda[j][0]
1576 *   + kappa_2(current_temperature_values[q], current_lambda_values[q]) * phi_lambda[j]
1577 *   )
1578 *  
1579 *   ) * fe_values.JxW(q);
1580 *  
1581 *   }
1582 *  
1583 *   row_last_element_vector(i) += (
1584 *   (phi_u[i] * current_velocity_gradients[q][0])
1585 *   + (phi_T[i] * current_temperature_gradients[q][0])
1586 *   - (phi_T[i] * current_velocity_gradients[q][0])
1587 *   + (phi_lambda[i] * current_lambda_gradients[q][0])
1588 *   ) * fe_values.JxW(q);
1589 *   }
1590 *  
1591 *   }
1592 *  
1593 *   cell->get_dof_indices(local_dof_indices);
1594 *  
1595 *   for (const unsigned int i : fe_values.dof_indices())
1596 *   {
1597 *   for (const unsigned int j : fe_values.dof_indices())
1598 *   {
1599 *   jacobian_matrix_extended.add(local_dof_indices[i],
1600 *   local_dof_indices[j],
1601 *   cell_matrix(i, j));
1602 *   }
1603 *  
1604 * @endcode
1605 *
1606 * Adding elements to the last column.
1607 *
1608 * @code
1609 *   jacobian_matrix_extended.add(local_dof_indices[i],
1610 *   extended_solution_dim - 1,
1611 *   row_last_element_vector(i));
1612 *   }
1613 *  
1614 *   }
1615 *  
1616 * @endcode
1617 *
1618 * Global dof indices of dofs for @f$dT@f$ and @f$d\lambda@f$, associated with vertex @f$\xi = 0@f$.
1619 *
1620 * @code
1621 *   types::global_dof_index T_zero_point_dof_ind(0), lambda_zero_point_dof_ind(0);
1622 *  
1623 * @endcode
1624 *
1625 * Approximating the derivative of @f$T@f$ at @f$\xi = 0@f$ as done in @ref step_14 "step-14".
1626 *
1627 * @code
1628 *   double T_point_derivative(0.);
1629 *   double T_at_zero_point(0.);
1630 *   double lambda_at_zero_point(0.);
1631 *   {
1632 *   double derivative_evaluation_point = 0.; // Point at which T = T_ign.
1633 *  
1634 *   const QTrapezoid<1> quadrature_formula;
1635 *   FEValues<1> fe_values(fe,
1636 *   quadrature_formula,
1638 *  
1639 *   const FEValuesExtractors::Scalar temperature(1);
1640 *   const FEValuesExtractors::Scalar lambda(2);
1641 *  
1642 *   const unsigned int n_q_points = quadrature_formula.size();
1643 *   std::vector<double> current_temperature_values(n_q_points);
1644 *   std::vector<Tensor<1, 1>> current_temperature_gradients(n_q_points);
1645 *   std::vector<double> current_lambda_values(n_q_points);
1646 *  
1647 *   unsigned int derivative_evaluation_point_hits = 0;
1648 *  
1649 *   for (const auto &cell : dof_handler.active_cell_iterators())
1650 *   {
1651 *   for (const auto &vertex : cell->vertex_indices())
1652 *   {
1653 *   if (isapprox(cell->vertex(vertex)[0], derivative_evaluation_point))
1654 *   {
1655 *   T_zero_point_dof_ind = cell->vertex_dof_index(vertex, 1);
1656 *   lambda_zero_point_dof_ind = cell->vertex_dof_index(vertex, 2);
1657 *  
1658 *   fe_values.reinit(cell);
1659 *   fe_values[temperature].get_function_values(current_solution, current_temperature_values);
1660 *   fe_values[temperature].get_function_gradients(current_solution, current_temperature_gradients);
1661 *   fe_values[lambda].get_function_values(current_solution, current_lambda_values);
1662 *  
1663 *   unsigned int q_point = 0;
1664 *   for (; q_point < n_q_points; ++q_point)
1665 *   {
1666 *   if (isapprox(fe_values.quadrature_point(q_point)[0], derivative_evaluation_point))
1667 *   {
1668 *   break;
1669 *   }
1670 *   }
1671 *  
1672 *   T_at_zero_point = current_temperature_values[q_point];
1673 *   lambda_at_zero_point = current_lambda_values[q_point];
1674 *  
1675 *   T_point_derivative += current_temperature_gradients[q_point][0];
1676 *   ++derivative_evaluation_point_hits;
1677 *  
1678 *   break;
1679 *   }
1680 *   }
1681 *   }
1682 *   T_point_derivative /= static_cast<double>(derivative_evaluation_point_hits);
1683 *   }
1684 *  
1685 * @endcode
1686 *
1687 * 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).
1688 *
1689 * @code
1690 *   double term_with_delta_func(0.);
1691 *   term_with_delta_func = problem.k * std::exp(-problem.theta / T_at_zero_point) * (1 - lambda_at_zero_point) / std::abs(T_point_derivative);
1692 *   jacobian_matrix_extended.add(T_zero_point_dof_ind, T_zero_point_dof_ind, problem.q * term_with_delta_func);
1693 *   jacobian_matrix_extended.add(lambda_zero_point_dof_ind, T_zero_point_dof_ind, term_with_delta_func);
1694 *  
1695 * @endcode
1696 *
1697 * Add 1 to the position <code> T_zero_point_dof_ind </code> of the last row of the matrix.
1698 *
1699 * @code
1700 *   jacobian_matrix_extended.add(extended_solution_dim - 1, T_zero_point_dof_ind, 1.);
1701 *  
1702 *   zero_boundary_constraints.condense(jacobian_matrix_extended);
1703 *   }
1704 *  
1705 *   {
1706 *   TimerOutput::Scope t(computing_timer, "factorizing the Jacobian");
1707 *  
1708 *   std::cout << "Factorizing Jacobian matrix" << std::endl;
1709 *  
1710 *   jacobian_matrix_extended_factorization = std::make_unique<SparseDirectUMFPACK>();
1711 *   jacobian_matrix_extended_factorization->factorize(jacobian_matrix_extended);
1712 *   }
1713 *  
1714 *   }
1715 *  
1716 *  
1717 *   double TravelingWaveSolver::compute_residual(const Vector<double> &evaluation_point_extended, Vector<double> &residual)
1718 *   {
1719 *   TimerOutput::Scope t(computing_timer, "assembling the residual");
1720 *  
1721 *   std::cout << "Computing residual vector ... " << std::endl;
1722 *  
1723 *   residual = 0;
1724 *  
1725 *   Vector<double> evaluation_point(dof_handler.n_dofs());
1726 *   for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
1727 *   {
1728 *   evaluation_point(i) = evaluation_point_extended(i);
1729 *   }
1730 *  
1731 *   const double wave_speed = evaluation_point_extended(extended_solution_dim - 1);
1732 *  
1733 *   const QGauss<1> quadrature_formula(number_of_quadrature_points);
1734 *   FEValues<1> fe_values(fe,
1735 *   quadrature_formula,
1738 *  
1739 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1740 *   const unsigned int n_q_points = quadrature_formula.size();
1741 *  
1742 *   Vector<double> cell_residual(dofs_per_cell);
1743 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1744 *  
1745 *   const FEValuesExtractors::Scalar velocity(0);
1746 *   const FEValuesExtractors::Scalar temperature(1);
1747 *   const FEValuesExtractors::Scalar lambda(2);
1748 *  
1749 *   std::vector<double> current_velocity_values(n_q_points);
1750 *   std::vector<Tensor<1, 1>> current_velocity_gradients(n_q_points);
1751 *   std::vector<double> current_temperature_values(n_q_points);
1752 *   std::vector<Tensor<1, 1>> current_temperature_gradients(n_q_points);
1753 *   std::vector<double> current_lambda_values(n_q_points);
1754 *   std::vector<Tensor<1, 1>> current_lambda_gradients(n_q_points);
1755 *  
1756 *   std::vector<double> phi_u(dofs_per_cell);
1757 *   std::vector<Tensor<1, 1>> grad_phi_u(dofs_per_cell);
1758 *   std::vector<double> phi_T(dofs_per_cell);
1759 *   std::vector<Tensor<1, 1>> grad_phi_T(dofs_per_cell);
1760 *   std::vector<double> phi_lambda(dofs_per_cell);
1761 *   std::vector<Tensor<1, 1>> grad_phi_lambda(dofs_per_cell);
1762 *  
1763 *   for (const auto &cell : dof_handler.active_cell_iterators())
1764 *   {
1765 *   cell_residual = 0;
1766 *  
1767 *   fe_values.reinit(cell);
1768 *  
1769 *   fe_values[velocity].get_function_values(evaluation_point, current_velocity_values);
1770 *   fe_values[velocity].get_function_gradients(evaluation_point, current_velocity_gradients);
1771 *   fe_values[temperature].get_function_values(evaluation_point, current_temperature_values);
1772 *   fe_values[temperature].get_function_gradients(evaluation_point, current_temperature_gradients);
1773 *   fe_values[lambda].get_function_values(evaluation_point, current_lambda_values);
1774 *   fe_values[lambda].get_function_gradients(evaluation_point, current_lambda_gradients);
1775 *  
1776 *   auto omega = [=](double T, double lambda){
1777 *   return problem.k * (1 - lambda) * std::exp(-problem.theta / T) * Heaviside_func(T - problem.T_ign);
1778 *   };
1779 *  
1780 *   for (unsigned int q = 0; q < n_q_points; ++q)
1781 *   {
1782 *   for (unsigned int k = 0; k < dofs_per_cell; ++k)
1783 *   {
1784 *   phi_u[k] = fe_values[velocity].value(k, q);
1785 *   grad_phi_u[k] = fe_values[velocity].gradient(k, q);
1786 *   phi_T[k] = fe_values[temperature].value(k, q);
1787 *   grad_phi_T[k] = fe_values[temperature].gradient(k, q);
1788 *   phi_lambda[k] = fe_values[lambda].value(k, q);
1789 *   grad_phi_lambda[k] = fe_values[lambda].gradient(k, q);
1790 *   }
1791 *  
1792 *   double del_Pr_eps = (problem.Pr * 4 * problem.delta / (3 * problem.epsilon));
1793 *   double del_Le = (problem.delta / problem.Le);
1794 *  
1795 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1796 *   {
1797 *   cell_residual(i) += (
1798 *  
1799 *   del_Pr_eps * (-grad_phi_u[i] * current_velocity_gradients[q])
1800 *   + phi_u[i] * (
1801 *   - current_velocity_gradients[q][0] * (1 - wave_speed + problem.epsilon * current_velocity_values[q])
1802 *   - problem.epsilon / 2. * current_temperature_gradients[q][0]
1803 *   )
1804 *  
1805 *   + problem.delta * (-grad_phi_T[i] * current_temperature_gradients[q])
1806 *   + phi_T[i] * (
1807 *   wave_speed * (current_temperature_gradients[q][0] - current_velocity_gradients[q][0])
1808 *   + problem.q * omega(current_temperature_values[q], current_lambda_values[q])
1809 *   )
1810 *  
1811 *   + del_Le * (-grad_phi_lambda[i] * current_lambda_gradients[q])
1812 *   + phi_lambda[i] * (
1813 *   wave_speed * current_lambda_gradients[q][0] + omega(current_temperature_values[q], current_lambda_values[q])
1814 *   )
1815 *  
1816 *   ) * fe_values.JxW(q);
1817 *   }
1818 *  
1819 *   }
1820 *  
1821 *   cell->get_dof_indices(local_dof_indices);
1822 *  
1823 *   for (const unsigned int i : fe_values.dof_indices())
1824 *   {
1825 *   residual(local_dof_indices[i]) += cell_residual(i);
1826 *   }
1827 *   }
1828 *  
1829 *   residual(extended_solution_dim - 1) = 0.;
1830 *  
1831 *   zero_boundary_constraints.condense(residual);
1832 *  
1833 *   double residual_norm = residual.l2_norm();
1834 *  
1835 *   std::cout << std::defaultfloat;
1836 *   std::cout << "norm of residual = " << residual_norm << std::endl;
1837 *  
1838 *   return residual_norm;
1839 *   }
1840 *  
1841 * @endcode
1842 *
1843 * 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.
1844 *
1845 * @code
1846 *   void TravelingWaveSolver::split_extended_solution_vector()
1847 *   {
1848 *   for (unsigned int i = 0; i < extended_solution_dim - 1; ++i)
1849 *   {
1850 *   current_solution(i) = current_solution_extended(i);
1851 *   }
1852 *  
1853 *   current_wave_speed = current_solution_extended(extended_solution_dim - 1);
1854 *   }
1855 *  
1856 *  
1857 *   void TravelingWaveSolver::solve(const Vector<double> &rhs, Vector<double> &solution_extended, const double /*tolerance*/)
1858 *   {
1859 *   TimerOutput::Scope t(computing_timer, "linear system solve");
1860 *  
1861 *   std::cout << "Solving linear system ... " << std::endl;
1862 *  
1863 *   jacobian_matrix_extended_factorization->vmult(solution_extended, rhs);
1864 *  
1865 *   zero_boundary_constraints.distribute(solution_extended);
1866 *  
1867 *   }
1868 *  
1869 *  
1870 * @endcode
1871 *
1872 * Function for adaptive mesh refinement based on <code> KellyErrorEstimator </code>.
1873 *
1874 * @code
1875 *   void TravelingWaveSolver::refine_mesh()
1876 *   {
1877 *   Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
1878 *  
1879 *   const FEValuesExtractors::Scalar lambda(2);
1880 *  
1882 *   dof_handler,
1883 *   QGauss<0>( 0 /* number_of_quadrature_points */),
1884 *   {},
1885 *   current_solution,
1886 *   estimated_error_per_cell,
1887 *   fe.component_mask(lambda)
1888 *   );
1889 *  
1891 *   estimated_error_per_cell,
1892 *   0.1,
1893 *   0.05);
1894 *  
1896 *  
1897 *   SolutionTransfer<1> solution_transfer(dof_handler);
1898 *   solution_transfer.prepare_for_coarsening_and_refinement(current_solution);
1899 *  
1901 *  
1902 *   setup_system(/*initial_step=*/ false);
1903 *  
1904 *   Vector<double> tmp(dof_handler.n_dofs());
1905 *   solution_transfer.interpolate(current_solution, tmp);
1906 *   current_solution = std::move(tmp);
1907 *  
1908 *   set_boundary_and_centering_values();
1909 *  
1910 *   for (unsigned int i = 0; i < extended_solution_dim - 1; ++i)
1911 *   {
1912 *   current_solution_extended(i) = current_solution(i);
1913 *   }
1914 *   current_solution_extended(extended_solution_dim - 1) = current_wave_speed;
1915 *  
1916 *   }
1917 *  
1918 *  
1919 *   double TravelingWaveSolver::run_newton_iterations(const double target_tolerance)
1920 *   {
1921 *  
1922 *   double residual_norm = 0.;
1923 *   {
1924 *   typename SUNDIALS::KINSOL< Vector<double> >::AdditionalData additional_data;
1925 *   additional_data.function_tolerance = target_tolerance;
1926 *  
1927 *   SUNDIALS::KINSOL<Vector<double>> nonlinear_solver(additional_data);
1928 *  
1929 *   nonlinear_solver.reinit_vector = [&](Vector<double> &x) {
1930 *   x.reinit(extended_solution_dim);
1931 *   };
1932 *  
1933 *   nonlinear_solver.residual = [&](const Vector<double> &evaluation_point, Vector<double> &residual) {
1934 *   residual_norm = compute_residual(evaluation_point, residual);
1935 *  
1936 *   return 0;
1937 *   };
1938 *  
1939 *   nonlinear_solver.setup_jacobian = [&](const Vector<double> &evaluation_point, const Vector<double> & /*current_f*/) {
1940 *   compute_and_factorize_jacobian(evaluation_point);
1941 *  
1942 *   return 0;
1943 *   };
1944 *  
1945 *   nonlinear_solver.solve_with_jacobian = [&](const Vector<double> &rhs, Vector<double> &solution, const double tolerance) {
1946 *   this->solve(rhs, solution, tolerance);
1947 *  
1948 *   return 0;
1949 *   };
1950 *  
1951 *   nonlinear_solver.solve(current_solution_extended);
1952 *   }
1953 *  
1954 *   return residual_norm;
1955 *  
1956 *   }
1957 *  
1958 * @endcode
1959 *
1960 * 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.
1961 *
1962 * @code
1963 *   void TravelingWaveSolver::output_with_double_precision(const Vector<double> &solution, const double wave_speed, const std::string filename)
1964 *   {
1965 *   TimerOutput::Scope t(computing_timer, "graphical output txt");
1966 *  
1967 *   const std::string file_for_solution = filename + ".txt";
1968 *   std::ofstream output(file_for_solution);
1969 *  
1970 *   for (const auto &cell : dof_handler.active_cell_iterators())
1971 *   {
1972 *   for (const auto &v_ind : cell->vertex_indices())
1973 *   {
1974 *   double u = solution(cell->vertex_dof_index(v_ind, 0));
1975 *   double T = solution(cell->vertex_dof_index(v_ind, 1));
1976 *   double lambda = solution(cell->vertex_dof_index(v_ind, 2));
1977 *  
1978 *   output << std::scientific << std::setprecision(16);
1979 *   output << cell->vertex(v_ind)[0];
1980 *  
1981 *   output << std::scientific << std::setprecision(16);
1982 *   output << std::scientific << " " << u << " " << T << " " << lambda << "\n";
1983 *   }
1984 *   output << "\n";
1985 *   }
1986 *  
1987 *   output.close();
1988 *  
1989 *   std::ofstream file_for_wave_speed_output("wave_speed-" + file_for_solution);
1990 *   file_for_wave_speed_output << std::scientific << std::setprecision(16);
1991 *   file_for_wave_speed_output << wave_speed << std::endl;
1992 *   file_for_wave_speed_output.close();
1993 *   }
1994 *  
1995 * @endcode
1996 *
1997 * Copy the solution into the <code> SolutionStruct </code> object, that stores the solution in an ordered manner.
1998 *
1999 * @code
2000 *   void TravelingWaveSolver::get_solution(SolutionStruct &solution) const
2001 *   {
2002 * @endcode
2003 *
2004 * 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.
2005 *
2006 * @code
2007 *   auto comp = [](const std::vector<double> &a, const std::vector<double> &b) {
2008 *   return a[0] < b[0];
2009 *   };
2010 *   std::set<std::vector<double>, decltype(comp)> solution_set(comp);
2011 *   for (const auto &cell : dof_handler.active_cell_iterators())
2012 *   {
2013 *   for (const auto &v_ind : cell->vertex_indices())
2014 *   {
2015 *   double x = cell->vertex(v_ind)[0];
2016 *   double u = current_solution(cell->vertex_dof_index(v_ind, 0));
2017 *   double T = current_solution(cell->vertex_dof_index(v_ind, 1));
2018 *   double lambda = current_solution(cell->vertex_dof_index(v_ind, 2));
2019 *   solution_set.insert({x, u, T, lambda});
2020 *   }
2021 *   }
2022 *  
2023 *   solution.x.clear();
2024 *   solution.u.clear();
2025 *   solution.T.clear();
2026 *   solution.lambda.clear();
2027 *  
2028 *   solution.x.reserve(solution_set.size());
2029 *   solution.u.reserve(solution_set.size());
2030 *   solution.T.reserve(solution_set.size());
2031 *   solution.lambda.reserve(solution_set.size());
2032 *  
2033 *   for (auto it = solution_set.begin(); it != solution_set.end(); ++it)
2034 *   {
2035 *   solution.x.push_back((*it)[0]);
2036 *   solution.u.push_back((*it)[1]);
2037 *   solution.T.push_back((*it)[2]);
2038 *   solution.lambda.push_back((*it)[3]);
2039 *   }
2040 *  
2041 *   solution.wave_speed = current_wave_speed;
2042 *  
2043 *   }
2044 *  
2045 *  
2046 *   void TravelingWaveSolver::get_triangulation(Triangulation<1> &otriangulation) const
2047 *   {
2048 *   otriangulation.clear();
2049 *   otriangulation.copy_triangulation(triangulation);
2050 *   }
2051 *  
2052 *  
2053 *   void TravelingWaveSolver::run(const std::string filename, const bool save_solution_to_file)
2054 *   {
2055 *   const int mesh_refinement_type = params.mesh.adaptive;
2056 *   const unsigned int n_refinements = params.mesh.refinements_number;
2057 *   const double tol = params.solver.tol;
2058 *  
2059 *   if (triangulation_uploaded == false) // If the triangulation is not loaded from outside, we will create one.
2060 *   {
2061 * @endcode
2062 *
2063 * 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.
2064 *
2065 * @code
2066 *   Triangulation<1> triangulation_left;
2068 *   triangulation_left,
2069 *   static_cast<unsigned int>(std::abs( 0. - params.mesh.interval_left )),
2070 *   params.mesh.interval_left, 0.
2071 *   );
2072 *  
2073 *   Triangulation<1> triangulation_right;
2075 *   triangulation_right,
2076 *   static_cast<unsigned int>(std::abs( params.mesh.interval_right - 0. )),
2077 *   0., params.mesh.interval_right
2078 *   );
2079 *  
2080 *   GridGenerator::merge_triangulations(triangulation_left, triangulation_right, triangulation);
2081 *  
2082 *   }
2083 *  
2084 *   if (triangulation_uploaded == false)
2085 *   {
2086 *   if (mesh_refinement_type == 1) // For ADAPTIVE mesh refinement.
2087 *   {
2088 *   triangulation.refine_global(1); // refine initial mesh globally, before adaptive refinement cycles.
2089 *   }
2090 *   else if (mesh_refinement_type == 0) // For GLOBAL mesh refinement.
2091 *   {
2092 *   triangulation.refine_global(n_refinements);
2093 *   }
2094 *   }
2095 *  
2096 *   setup_system(/*initial step*/ true);
2097 *   set_initial_guess();
2098 *  
2099 *   if (save_solution_to_file)
2100 *   {
2101 *   output_with_double_precision(current_solution, current_wave_speed, "solution_initial_data");
2102 *   }
2103 *  
2104 *   if (mesh_refinement_type == 1) // Compute with ADAPTIVE mesh refinement.
2105 *   {
2106 *   double residual_norm = 0.;
2107 *   {
2108 *   Vector<double> tmp_residual(extended_solution_dim);
2109 *   residual_norm = compute_residual(current_solution_extended, tmp_residual);
2110 *   }
2111 *  
2112 *   unsigned int refinement_cycle = 0;
2113 *   while ((residual_norm > tol) && (refinement_cycle < n_refinements))
2114 *   {
2115 *   computing_timer.reset();
2116 *   std::cout << "Mesh refinement step " << refinement_cycle << std::endl;
2117 *  
2118 *   if (refinement_cycle != 0) { refine_mesh(); }
2119 *  
2120 *   const double target_tolerance = 0.1 * std::pow(0.1, refinement_cycle); // Decrease tolerance for Newton solver at each refinement step.
2121 *   std::cout << " Target_tolerance: " << target_tolerance << std::endl;
2122 *  
2123 *   residual_norm = run_newton_iterations(target_tolerance);
2124 *   split_extended_solution_vector();
2125 *  
2126 *   {
2127 *   std::cout << std::scientific << std::setprecision(16);
2128 *   std::cout << "current_wave_speed = " << current_wave_speed << std::endl;
2129 *   std::cout << std::defaultfloat;
2130 *   }
2131 *  
2132 *   computing_timer.print_summary();
2133 *  
2134 *   ++refinement_cycle;
2135 *   }
2136 *   if (save_solution_to_file)
2137 *   {
2138 *   output_with_double_precision(current_solution, current_wave_speed, filename);
2139 *   }
2140 *  
2141 *   }
2142 *   else if (mesh_refinement_type == 0) // Compute with GLOBAL mesh refinement.
2143 *   {
2144 *   run_newton_iterations(tol);
2145 *   split_extended_solution_vector();
2146 *  
2147 *   if (save_solution_to_file)
2148 *   {
2149 *   output_with_double_precision(current_solution, current_wave_speed, filename);
2150 *   }
2151 *  
2152 *   {
2153 *   std::cout << std::scientific << std::setprecision(16);
2154 *   std::cout << "current_wave_speed = " << current_wave_speed << std::endl;
2155 *   std::cout << std::defaultfloat;
2156 *   }
2157 *  
2158 *   computing_timer.print_summary();
2159 *  
2160 *   }
2161 *  
2162 *   }
2163 *  
2164 *  
2165 *   } // namespace TravelingWave
2166 * @endcode
2167
2168
2169<a name="ann-TravelingWaveSolver.h"></a>
2170<h1>Annotated version of TravelingWaveSolver.h</h1>
2171 *
2172 *
2173 *
2174 *
2175 * @code
2176 *   /* -----------------------------------------------------------------------------
2177 *   *
2178 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
2179 *   * Copyright (C) 2024 by Shamil Magomedov
2180 *   *
2181 *   * This file is part of the deal.II code gallery.
2182 *   *
2183 *   * -----------------------------------------------------------------------------
2184 *   */
2185 *  
2186 *   #ifndef TRAVELING_WAVE_SOLVER
2187 *   #define TRAVELING_WAVE_SOLVER
2188 *  
2189 *   #include <deal.II/base/timer.h>
2190 *   #include <deal.II/base/function.h>
2191 *   #include <deal.II/base/table_handler.h>
2192 *   #include <deal.II/grid/tria.h>
2193 *   #include <deal.II/grid/grid_generator.h>
2194 *   #include <deal.II/grid/grid_refinement.h>
2195 *   #include <deal.II/grid/grid_out.h>
2196 *   #include <deal.II/dofs/dof_handler.h>
2197 *   #include <deal.II/dofs/dof_tools.h>
2198 *   #include <deal.II/fe/fe_q.h>
2199 *   #include <deal.II/fe/fe_system.h>
2200 *   #include <deal.II/fe/fe_values.h>
2201 *   #include <deal.II/lac/vector.h>
2202 *   #include <deal.II/lac/full_matrix.h>
2203 *   #include <deal.II/lac/sparse_matrix.h>
2204 *   #include <deal.II/lac/sparse_direct.h>
2205 *   #include <deal.II/lac/dynamic_sparsity_pattern.h>
2206 *   #include <deal.II/numerics/vector_tools.h>
2207 *   #include <deal.II/numerics/error_estimator.h>
2208 *   #include <deal.II/numerics/solution_transfer.h>
2209 *   #include <deal.II/numerics/data_out.h>
2210 *   #include <deal.II/sundials/kinsol.h>
2211 *  
2212 *   #include "Parameters.h"
2213 *   #include "Solution.h"
2214 *   #include "AuxiliaryFunctions.h"
2215 *  
2216 *   #include <cmath>
2217 *   #include <algorithm>
2218 *   #include <string>
2219 *   #include <fstream>
2220 *   #include <iostream>
2221 *   #include <iomanip>
2222 *   #include <set>
2223 *  
2224 * @endcode
2225 *
2226 * Namespace of the program
2227 *
2228 * @code
2229 *   namespace TravelingWave
2230 *   {
2231 *   using namespace dealii;
2232 *  
2233 * @endcode
2234 *
2235 * The main class for construction of the traveling wave solutions.
2236 *
2237 * @code
2238 *   class TravelingWaveSolver
2239 *   {
2240 *   public:
2241 *   TravelingWaveSolver(const Parameters &parameters, const SolutionStruct &initial_guess_input);
2242 *  
2243 *   void set_triangulation(const Triangulation<1> &itriangulation);
2244 *  
2245 *   void run(const std::string filename="solution", const bool save_solution_to_file=true);
2246 *   void get_solution(SolutionStruct &solution) const;
2247 *   void get_triangulation(Triangulation<1> &otriangulation) const;
2248 *  
2249 *   private:
2250 *   void setup_system(const bool initial_step);
2251 *   void find_boundary_and_centering_dof_numbers();
2252 *   void set_boundary_and_centering_values();
2253 *  
2254 *   void set_initial_guess();
2255 *  
2256 *   double Heaviside_func(double x) const;
2257 *  
2258 *   void compute_and_factorize_jacobian(const Vector<double> &evaluation_point_extended);
2259 *   double compute_residual(const Vector<double> &evaluation_point_extended, Vector<double> &residual);
2260 *   void split_extended_solution_vector();
2261 *  
2262 *   void solve(const Vector<double> &rhs, Vector<double> &solution, const double /*tolerance*/);
2263 *   void refine_mesh();
2264 *   double run_newton_iterations(const double target_tolerance=1e-5);
2265 *  
2266 *   void output_with_double_precision(const Vector<double> &solution, const double wave_speed, const std::string filename="solution");
2267 *  
2268 * @endcode
2269 *
2270 * The dimension of the finite element solution increased by one to account for the value corresponding to the wave speed.
2271 *
2272 * @code
2273 *   unsigned int extended_solution_dim;
2274 *   std::map<std::string, unsigned int> boundary_and_centering_dof_numbers;
2275 *  
2276 * @endcode
2277 *
2278 * Parameters of the problem, taken from a .prm file.
2279 *
2280 * @code
2281 *   const Parameters &params;
2282 *   const Problem &problem; // Reference variable, just for convenience.
2283 *  
2284 *   unsigned int number_of_quadrature_points;
2285 *  
2287 * @endcode
2288 *
2289 * The flag indicating whether the triangulation was uploaded externally or created within the <code> run </code> member function.
2290 *
2291 * @code
2292 *   bool triangulation_uploaded;
2293 *   FESystem<1> fe;
2294 *   DoFHandler<1> dof_handler;
2295 *  
2296 * @endcode
2297 *
2298 * Constraints for Dirichlet boundary conditions.
2299 *
2300 * @code
2301 *   AffineConstraints<double> zero_boundary_constraints;
2302 *  
2303 *   SparsityPattern sparsity_pattern_extended;
2304 *   SparseMatrix<double> jacobian_matrix_extended;
2305 *   std::unique_ptr<SparseDirectUMFPACK> jacobian_matrix_extended_factorization;
2306 *  
2307 * @endcode
2308 *
2309 * Finite element solution of the problem.
2310 *
2311 * @code
2312 *   Vector<double> current_solution;
2313 *  
2314 * @endcode
2315 *
2316 * Value of the wave speed @f$c@f$.
2317 *
2318 * @code
2319 *   double current_wave_speed;
2320 *  
2321 * @endcode
2322 *
2323 * Solution with an additional term, corresponding to the variable wave_speed.
2324 *
2325 * @code
2326 *   Vector<double> current_solution_extended;
2327 *  
2328 * @endcode
2329 *
2330 * Initial guess for Newton's iterations.
2331 *
2332 * @code
2333 *   SolutionStruct initial_guess;
2334 *  
2335 *   TimerOutput computing_timer;
2336 *   };
2337 *  
2338 *   } // namespace TravelingWave
2339 *  
2340 *   #endif
2341 * @endcode
2342
2343
2344<a name="ann-calculate_profile.cc"></a>
2345<h1>Annotated version of calculate_profile.cc</h1>
2346 *
2347 *
2348 *
2349 *
2350 * @code
2351 *   /* -----------------------------------------------------------------------------
2352 *   *
2353 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
2354 *   * Copyright (C) 2024 by Shamil Magomedov
2355 *   *
2356 *   * This file is part of the deal.II code gallery.
2357 *   *
2358 *   * -----------------------------------------------------------------------------
2359 *   */
2360 *  
2361 *   #include "TravelingWaveSolver.h"
2362 *   #include "calculate_profile.h"
2363 *  
2364 *   namespace TravelingWave
2365 *   {
2366 *   using namespace dealii;
2367 *  
2368 * @endcode
2369 *
2370 * 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.
2371 *
2372 * @code
2373 *   void compute_limit_sol_left_part(const Parameters &parameters,
2374 *   const double wave_speed,
2375 *   const double u_0,
2376 *   const double T_0,
2377 *   const double lambda_0,
2378 *   SolutionStruct &LimitSol,
2379 *   const double root_sign)
2380 *   {
2381 *   LimitSolution limit_sol(parameters, lambda_0, u_0, T_0, root_sign);
2382 *   limit_sol.set_wave_speed(wave_speed);
2383 *  
2384 *   {
2385 * @endcode
2386 *
2387 * We take more integration points to better resolve the transition layer.
2388 *
2389 * @code
2390 *   std::vector<double> t_span(static_cast<unsigned int>(std::abs( 0. - parameters.mesh.interval_left )));
2391 *   double finer_mesh_starting_value = -9.1;
2392 *   linspace(parameters.mesh.interval_left, finer_mesh_starting_value, t_span);
2393 *   std::vector<double> fine_grid(10000);
2394 *   linspace(finer_mesh_starting_value + 1e-4, 0., fine_grid);
2395 *   t_span.insert(t_span.end(), fine_grid.begin(), fine_grid.end());
2396 *  
2397 * @endcode
2398 *
2399 * Reverse the order of the elements (because we need to perform back in time integration).
2400 *
2401 * @code
2402 *   std::reverse(t_span.begin(), t_span.end());
2403 *  
2404 *   state_type lambda_val(1);
2405 *   lambda_val[0] = lambda_0; // initial value
2406 *   IntegrateSystemAtTimePoints(limit_sol.lambda_vec, limit_sol.t_vec, t_span,
2407 *   limit_sol,
2408 *   lambda_val,
2409 *   -1e-2, Integrator_Type::dopri5);
2410 *   }
2411 *  
2412 *   limit_sol.calculate_u_T_omega();
2413 *  
2414 * @endcode
2415 *
2416 * Reverse the order of elements
2417 *
2418 * @code
2419 *   std::reverse(limit_sol.t_vec.begin(), limit_sol.t_vec.end());
2420 *   std::reverse(limit_sol.lambda_vec.begin(), limit_sol.lambda_vec.end());
2421 *   std::reverse(limit_sol.u_vec.begin(), limit_sol.u_vec.end());
2422 *   std::reverse(limit_sol.T_vec.begin(), limit_sol.T_vec.end());
2423 *   std::reverse(limit_sol.omega_vec.begin(), limit_sol.omega_vec.end());
2424 *  
2425 *   SaveSolutionIntoFile(limit_sol.lambda_vec, limit_sol.t_vec, "solution_lambda_limit.txt");
2426 *   SaveSolutionIntoFile(limit_sol.u_vec, limit_sol.t_vec, "solution_u_limit.txt");
2427 *   SaveSolutionIntoFile(limit_sol.T_vec, limit_sol.t_vec, "solution_T_limit.txt");
2428 *   SaveSolutionIntoFile(limit_sol.omega_vec, limit_sol.t_vec, "solution_omega_limit.txt");
2429 *  
2430 *   LimitSol.reinit(limit_sol.t_vec.size());
2431 *   LimitSol.wave_speed = wave_speed;
2432 *   for (unsigned int i=0; i < limit_sol.t_vec.size(); ++i)
2433 *   {
2434 *   LimitSol.x[i] = limit_sol.t_vec[i];
2435 *   LimitSol.u[i] = limit_sol.u_vec[i][0];
2436 *   LimitSol.T[i] = limit_sol.T_vec[i][0];
2437 *   LimitSol.lambda[i] = limit_sol.lambda_vec[i][0];
2438 *   }
2439 *   }
2440 *  
2441 *  
2442 * @endcode
2443 *
2444 * Construction of an initial guess for detonation wave solution. The ODE is solved for the ideal system with @f$\delta = 0@f$.
2445 *
2446 * @code
2447 *   void compute_initial_guess_detonation(const Parameters &params, SolutionStruct &initial_guess, const double root_sign)
2448 *   {
2449 *   const Problem &problem = params.problem;
2450 *   double current_wave_speed(problem.wave_speed_init);
2451 *  
2452 *   { // 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.
2453 *   double DeltaT = problem.T_left - problem.T_right;
2454 *   double qDT = problem.q - DeltaT;
2455 *   current_wave_speed = 1. + problem.epsilon * (problem.u_right - (qDT * qDT + DeltaT) / (2 * qDT));
2456 *   }
2457 *  
2458 *   double u_0 = problem.u_right;
2459 *   double T_0 = problem.T_right;
2460 *   double lambda_0 = 0.;
2461 *  
2462 *   compute_limit_sol_left_part(params, current_wave_speed, u_0, T_0, lambda_0, initial_guess, root_sign);
2463 *  
2464 *   initial_guess.wave_speed = current_wave_speed;
2465 *  
2466 *   for (int i = initial_guess.x.size() - 1; i > - 1; --i)
2467 *   {
2468 *   if (isapprox(initial_guess.x[i], 0.))
2469 *   {
2470 *   initial_guess.u[i] = problem.u_right;
2471 *   initial_guess.T[i] = problem.T_ign;
2472 *   initial_guess.lambda[i] = 0.;
2473 *   break;
2474 *   }
2475 *   }
2476 *  
2477 * @endcode
2478 *
2479 * Adding the points to the right part of the interval (w.r.t. @f$\xi = 0@f$).
2480 *
2481 * @code
2482 *   unsigned int number_of_additional_points = 5;
2483 *   for (unsigned int i = 0; i < number_of_additional_points; ++i)
2484 *   {
2485 *   initial_guess.x.push_back(params.mesh.interval_right / (std::pow(2., number_of_additional_points - 1 - i)));
2486 *   initial_guess.u.push_back(problem.u_right);
2487 *   initial_guess.T.push_back(problem.T_right);
2488 *   initial_guess.lambda.push_back(0.);
2489 *   }
2490 *  
2491 *   }
2492 *  
2493 *  
2494 * @endcode
2495 *
2496 * Construction of a piecewise constant initial guess for deflagration wave solution.
2497 *
2498 * @code
2499 *   void compute_initial_guess_deflagration(const Parameters &params, SolutionStruct &initial_guess)
2500 *   {
2501 *   const Problem &problem = params.problem;
2502 *   double current_wave_speed(problem.wave_speed_init);
2503 *  
2504 *   double del_Pr_eps = (problem.Pr * 4 * problem.delta / (3 * problem.epsilon));
2505 *   double del_Le = (problem.delta / problem.Le);
2506 *  
2507 *   auto u_init_guess_func = [&](double x) {
2508 *   if (x < 0.)
2509 *   {
2510 *   return problem.u_left;
2511 *   }
2512 *   else
2513 *   {
2514 *   return problem.u_right;
2515 *   }
2516 *   };
2517 *  
2518 *   auto T_init_guess_func = [&](double x) {
2519 *   if (x < 0.)
2520 *   {
2521 *   return problem.T_left;
2522 *   }
2523 *   else if (isapprox(x, 0.))
2524 *   {
2525 *   return problem.T_ign;
2526 *   }
2527 *   else
2528 *   {
2529 *   return problem.T_right;
2530 *   }
2531 *   };
2532 *  
2533 *   auto lambda_init_guess_func = [=](double x) {
2534 *   if (x < 0.)
2535 *   {
2536 *   return -std::exp(x * std::abs(1 - current_wave_speed) / del_Pr_eps) + 1;
2537 *   }
2538 *   else
2539 *   {
2540 *   return 0.;
2541 *   }
2542 *   };
2543 *  
2544 *   unsigned int multiplier_for_number_of_points = 7;
2545 *   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 )));
2546 *   std::vector<double> x_span(number_of_points);
2547 *   linspace(params.mesh.interval_left, params.mesh.interval_right, x_span);
2548 *  
2549 *   std::vector<double> u_init_arr(number_of_points);
2550 *   std::vector<double> T_init_arr(number_of_points);
2551 *   std::vector<double> lambda_init_arr(number_of_points);
2552 *  
2553 *   for (unsigned int i = 0; i < number_of_points; ++i)
2554 *   {
2555 *   u_init_arr[i] = u_init_guess_func(x_span[i]);
2556 *   T_init_arr[i] = T_init_guess_func(x_span[i]);
2557 *   lambda_init_arr[i] = lambda_init_guess_func(x_span[i]);
2558 *   }
2559 *  
2560 *   initial_guess.x = x_span;
2561 *   initial_guess.u = u_init_arr;
2562 *   initial_guess.T = T_init_arr;
2563 *   initial_guess.lambda = lambda_init_arr;
2564 *   initial_guess.wave_speed = current_wave_speed;
2565 *  
2566 *   }
2567 *  
2568 *  
2569 * @endcode
2570 *
2571 * 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>.
2572 *
2573 * @code
2574 *   void calculate_profile(Parameters& parameters,
2575 *   const bool continuation_for_delta /* Compute with the continuation. */,
2576 *   const double delta_start /* The starting value of delta for the continuation method. */,
2577 *   const unsigned int number_of_continuation_points)
2578 *   {
2579 *   SolutionStruct sol;
2580 *  
2581 *   if (parameters.problem.wave_type == 1) // detonation wave
2582 *   {
2583 *   compute_initial_guess_detonation(parameters, sol);
2584 *   }
2585 *   else if (parameters.problem.wave_type == 0) // deflagration wave
2586 *   {
2587 *   compute_initial_guess_deflagration(parameters, sol);
2588 *   }
2589 *  
2590 *   if (continuation_for_delta == false)
2591 *   {
2592 *   TravelingWaveSolver wave(parameters, sol);
2593 *   std::string filename = "solution_delta-" + Utilities::to_string(parameters.problem.delta) + "_eps-"
2594 *   + Utilities::to_string(parameters.problem.epsilon);
2595 *   wave.run(filename);
2596 *   wave.get_solution(sol);
2597 *   }
2598 *   else // Run with continuation_for_delta.
2599 *   {
2600 *   double delta_target = parameters.problem.delta;
2601 *   parameters.problem.delta = delta_start;
2602 *  
2603 *   std::vector<double> delta_span(number_of_continuation_points);
2604 *  
2605 * @endcode
2606 *
2607 * Generate a sequence of delta values being uniformly distributed in log10 scale.
2608 *
2609 * @code
2610 *   {
2611 *   double delta_start_log10 = std::log10(delta_start);
2612 *   double delta_target_log10 = std::log10(delta_target);
2613 *  
2614 *   std::vector<double> delta_log_span(delta_span.size());
2615 *   linspace(delta_start_log10, delta_target_log10, delta_log_span);
2616 *  
2617 *   for (unsigned int i = 0; i < delta_span.size(); ++i)
2618 *   {
2619 *   delta_span[i] = std::pow(10, delta_log_span[i]);
2620 *   }
2621 *   }
2622 *  
2623 *   Triangulation<1> refined_triangulation;
2624 *   bool first_iter_flag = true;
2625 *  
2626 *   for (double delta : delta_span)
2627 *   {
2628 *   parameters.problem.delta = delta;
2629 *   std::string filename = "solution_delta-" + Utilities::to_string(parameters.problem.delta) + "_eps-"
2630 *   + Utilities::to_string(parameters.problem.epsilon);
2631 *  
2632 *   TravelingWaveSolver wave(parameters, sol);
2633 *  
2634 *   if (first_iter_flag)
2635 *   {
2636 *   first_iter_flag = false;
2637 *   }
2638 *   else
2639 *   {
2640 *   wave.set_triangulation(refined_triangulation);
2641 *   }
2642 *  
2643 *   wave.run(filename);
2644 *   wave.get_solution(sol);
2645 *   wave.get_triangulation(refined_triangulation);
2646 *   }
2647 *  
2648 *   }
2649 *  
2650 * @endcode
2651 *
2652 * Error estimation.
2653 *
2654 * @code
2655 *   {
2656 *   unsigned int sol_length = sol.x.size();
2657 *   double u_r = sol.u[sol_length-1]; // Dirichlet boundary condition
2658 *   double T_r = sol.T[sol_length-1]; // Dirichlet condition only for detonation case
2659 *   double u_l = sol.u[0];
2660 *   double T_l = sol.T[0]; // Dirichlet boundary condition
2661 *   double wave_speed = sol.wave_speed;
2662 *  
2663 *   std::cout << "Error estimates:" << std::endl;
2664 *   double DeltaT = T_l - T_r;
2665 *   double qDT = parameters.problem.q - DeltaT;
2666 *  
2667 *   double wave_speed_formula = 1. + parameters.problem.epsilon * (u_r - (qDT * qDT + DeltaT) / (2 * qDT));
2668 *   std::cout << std::setw(18) << std::left << "For wave speed" << " : " << std::setw(5) << wave_speed - wave_speed_formula << std::endl;
2669 *  
2670 *   double u_l_formula = DeltaT + u_r - parameters.problem.q;
2671 *   std::cout << std::setw(18) << std::left << "For u_l" << " : " << std::setw(5) << u_l - u_l_formula << std::endl;
2672 *   }
2673 *  
2674 *   }
2675 *  
2676 *   } // namespace TravelingWave
2677 * @endcode
2678
2679
2680<a name="ann-calculate_profile.h"></a>
2681<h1>Annotated version of calculate_profile.h</h1>
2682 *
2683 *
2684 *
2685 *
2686 * @code
2687 *   /* -----------------------------------------------------------------------------
2688 *   *
2689 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
2690 *   * Copyright (C) 2024 by Shamil Magomedov
2691 *   *
2692 *   * This file is part of the deal.II code gallery.
2693 *   *
2694 *   * -----------------------------------------------------------------------------
2695 *   */
2696 *  
2697 *   #ifndef CALCULATE_PROFILE
2698 *   #define CALCULATE_PROFILE
2699 *  
2700 *   #include "Parameters.h"
2701 *   #include "Solution.h"
2702 *   #include "LimitSolution.h"
2703 *   #include "IntegrateSystem.h"
2704 *   #include "AuxiliaryFunctions.h"
2705 *  
2706 *   namespace TravelingWave
2707 *   {
2708 *   void compute_limit_sol_left_part(const Parameters &parameters,
2709 *   const double wave_speed,
2710 *   const double u_0,
2711 *   const double T_0,
2712 *   const double lambda_0,
2713 *   SolutionStruct &LimitSol,
2714 *   const double root_sign = 1.);
2715 *  
2716 *   void compute_initial_guess_detonation(const Parameters &params, SolutionStruct &initial_guess, const double root_sign = 1.);
2717 *   void compute_initial_guess_deflagration(const Parameters &params, SolutionStruct &initial_guess);
2718 *  
2719 *   void calculate_profile(Parameters& parameters,
2720 *   const bool continuation_for_delta=false /* Compute with the continuation. */,
2721 *   const double delta_start=0.01 /* The starting value of delta for the continuation method. */,
2722 *   const unsigned int number_of_continuation_points=10);
2723 *  
2724 *   } // namespace TravelingWave
2725 *  
2726 *   #endif
2727 * @endcode
2728
2729
2730<a name="ann-main.cc"></a>
2731<h1>Annotated version of main.cc</h1>
2732 *
2733 *
2734 *
2735 *
2736 * @code
2737 *   /* -----------------------------------------------------------------------------
2738 *   *
2739 *   * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
2740 *   * Copyright (C) 2024 by Shamil Magomedov
2741 *   *
2742 *   * This file is part of the deal.II code gallery.
2743 *   *
2744 *   * -----------------------------------------------------------------------------
2745 *   */
2746 *  
2747 *   #include "calculate_profile.h"
2748 *  
2749 *   int main(int argc, char *argv[])
2750 *   {
2751 *  
2752 *   try
2753 *   {
2754 *   using namespace TravelingWave;
2755 *  
2756 *   Parameters parameters;
2757 *  
2758 *   std::string prm_filename = "ParametersList.prm";
2759 *   if (argc > 1)
2760 *   {
2761 * @endcode
2762 *
2763 * Check if file argv[1] exists.
2764 *
2765 * @code
2766 *   if (file_exists(argv[1]))
2767 *   {
2768 *   prm_filename = argv[1];
2769 *   }
2770 *   else
2771 *   {
2772 *   std::string errorMessage = "File \"" + std::string(argv[1]) + "\" is not found.";
2773 *   throw std::runtime_error(errorMessage);
2774 *   }
2775 *   }
2776 *   else
2777 *   {
2778 * @endcode
2779 *
2780 * Check if the file "ParametersList.prm" exists in the current or in the parent directory.
2781 *
2782 * @code
2783 *   if (!(file_exists(prm_filename) || file_exists("../" + prm_filename)))
2784 *   {
2785 *   std::string errorMessage = "File \"" + prm_filename + "\" is not found.";
2786 *   throw std::runtime_error(errorMessage);
2787 *   }
2788 *   else
2789 *   {
2790 *   if (!file_exists(prm_filename))
2791 *   {
2792 *   prm_filename = "../" + prm_filename;
2793 *   }
2794 *   }
2795 *   }
2796 *  
2797 *   std::cout << "Reading parameters... " << std::flush;
2798 *   ParameterAcceptor::initialize(prm_filename);
2799 *   std::cout << "done" << std::endl;
2800 *  
2801 *   calculate_profile(parameters, /* With continuation_for_delta */ false, 0.1, 3);
2802 *  
2803 *   }
2804 *   catch (std::exception &exc)
2805 *   {
2806 *   std::cerr << std::endl
2807 *   << std::endl
2808 *   << "----------------------------------------------------"
2809 *   << std::endl;
2810 *   std::cerr << "Exception on processing: " << std::endl
2811 *   << exc.what() << std::endl
2812 *   << "Aborting!" << std::endl
2813 *   << "----------------------------------------------------"
2814 *   << std::endl;
2815 *   return 1;
2816 *   }
2817 *   catch (...)
2818 *   {
2819 *   std::cerr << std::endl
2820 *   << std::endl
2821 *   << "----------------------------------------------------"
2822 *   << std::endl;
2823 *   std::cerr << "Unknown exception!" << std::endl
2824 *   << "Aborting!" << std::endl
2825 *   << "----------------------------------------------------"
2826 *   << std::endl;
2827 *   return 1;
2828 *   }
2829 *  
2830 *   return 0;
2831 *   }
2832 * @endcode
2833
2834
2835<a name="ann-plot.py"></a>
2836<h1>Annotated version of plot.py</h1>
2837@code{.py}
2838## -----------------------------------------------------------------------------
2839##
2840## SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
2841## Copyright (C) 2024 by Shamil Magomedov
2842##
2843## This file is part of the deal.II code gallery.
2844##
2845## -----------------------------------------------------------------------------
2846##
2847
2848import numpy as np
2849import matplotlib.pyplot as plt
2850import os
2851import sys
2852
2853plot_params = {
2854 #'backend': 'pdf',
2855 # 'lines.marker' : 'x',
2856 'scatter.marker' : 'x',
2857 'lines.markersize' : 4,
2858 'lines.linewidth' : 1,
2859 'axes.labelsize': 16,
2860 # 'textfontsize': 12,
2861 'font.size' : 16,
2862 'legend.fontsize': 16,
2863 'xtick.labelsize': 14,
2864 'ytick.labelsize': 14,
2865 'text.usetex': True,
2866 'figure.figsize': [9,6],
2867 'axes.grid': True
2868}
2869
2870plt.rcParams.update(plot_params)
2871
2872
2873if len(sys.argv) > 1:
2874
2875 filename = sys.argv[1]
2876
2877 if os.path.exists(filename):
2878 data = np.loadtxt(filename, np.float64)
2879 data_unique = np.unique(data, axis=0)
2880 data_unique = np.array(sorted(data_unique, key=lambda x : x[0]))
2881 x = data_unique[:, 0]
2882 u_sol = data_unique[:, 1]
2883 T_sol = data_unique[:, 2]
2884 lambda_sol = data_unique[:, 3]
2885
2886 fig, ax = plt.subplots(nrows=1, ncols=1)
2887
2888 ax.scatter(x, u_sol, label=r"@f$u@f$", color='blue')
2889 ax.scatter(x, T_sol, label=r"@f$T@f$", color='red')
2890 ax.scatter(x, lambda_sol, label=r"@f$\lambda@f$", color='green')
2891
2892
2893 # Plot of limit solutions for the detonation case. Uncomment, if needed.
2894 #===============================================================#
2895 '''
2896
2897 path_to_solution_files = os.path.split(filename)[0]
2898 u_limit_path = os.path.join(path_to_solution_files, 'solution_u_limit.txt')
2899 T_limit_path = os.path.join(path_to_solution_files, 'solution_T_limit.txt')
2900 lambda_limit_path = os.path.join(path_to_solution_files, 'solution_lambda_limit.txt')
2901
2902 if os.path.exists(u_limit_path):
2903 u_limit = np.loadtxt(u_limit_path, np.float64)
2904 ax.plot(u_limit[:, 0], u_limit[:, 1], label=r"@f$u_{\mathrm{lim}}@f$", color='blue')
2905 ax.plot([0, x[-1]], [u_sol[-1], u_sol[-1]], color='blue')
2906 else:
2907 print("No such file:", u_limit_path)
2908
2909 if os.path.exists(T_limit_path):
2910 T_limit = np.loadtxt(T_limit_path, np.float64)
2911 ax.plot(T_limit[:, 0], T_limit[:, 1], label=r"@f$T_{\mathrm{lim}}@f$", color='red')
2912 ax.plot([0, x[-1]], [T_sol[-1], T_sol[-1]], color='red')
2913 else:
2914 print("No such file:", T_limit_path)
2915
2916 if os.path.exists(lambda_limit_path):
2917 lambda_limit = np.loadtxt(lambda_limit_path, np.float64)
2918 ax.plot(lambda_limit[:, 0], lambda_limit[:, 1], label=r"@f$\lambda_{\mathrm{lim}}@f$", color='green')
2919 ax.plot([0, x[-1]], [lambda_sol[-1], lambda_sol[-1]], color='green')
2920 else:
2921 print("No such file:", lambda_limit_path)
2922
2923
2924 '''
2925 #===============================================================#
2926
2927
2928 ax.set_xlabel(r"@f$\xi@f$")
2929 ax.set_ylabel(r"@f$u, T, \lambda@f$")
2930 ax.legend()
2931
2932 # plt.savefig("fast_deflagration_delta_0.01.png", bbox_inches='tight', dpi=500)
2933 # plt.savefig('slow_deflagration_delta_0.01.png', bbox_inches='tight', dpi=500)
2934 # plt.savefig('detonation_delta_0.01.png', bbox_inches='tight', dpi=500)
2935
2936 plt.show()
2937 else:
2938 print("No such file:", filename)
2939@endcode
2940
2941
2942*/
Definition fe_q.h:554
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
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)
Definition point.h:111
void add_value(const std::string &key, const T value)
unsigned int n_active_cells() const
void refine_global(const unsigned int times=1)
virtual void execute_coarsening_and_refinement() override
Definition tria.cc:3320
virtual void copy_triangulation(const ::Triangulation< dim, spacedim > &other_tria) override
Definition tria.cc:4004
virtual bool prepare_coarsening_and_refinement() override
Definition tria.cc:2805
virtual void clear() override
Definition tria.cc:1864
Point< 3 > vertices[4]
Point< 2 > first
Definition grid_out.cc:4623
unsigned int vertex_indices[2]
__global__ void set(Number *val, const Number s, const size_type N)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
const Event initial
Definition event.cc:64
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)
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())
@ matrix
Contents is actually a matrix.
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
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
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
T scatter(const MPI_Comm comm, const std::vector< T > &objects_to_send, const unsigned int root_process=0)
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={})
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
DEAL_II_HOST constexpr TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
void copy(const T *begin, const T *end, U *dest)
int(& functions)(const void *v1, const void *v2)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
STL namespace.
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation