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