Traveling waves

TravelingWaves/AuxiliaryFunctions.h

@@ -0,0 +1,36 @@
+#ifndef AUXILIARY_FUNCTIONS
+#define AUXILIARY_FUNCTIONS
+
+#include <vector>
+#include <cmath>
+#include <fstream>
+#include <string>
+#include <cstring>
+
+// Comparison of numbers with a given tolerance.
+template <typename T>
+bool isapprox(const T &a, const T &b, const double tol = 1e-10)
+{
+  return (std::abs( a - b ) < tol);
+}
+
+// Fill the std::vector with the values from the range [interval_begin, interval_end].
+template <typename T>
+void linspace(T interval_begin, T interval_end, std::vector<T> &arr) 
+{
+  const size_t SIZE = arr.size();
+  const T step = (interval_end - interval_begin) / static_cast<T>(SIZE - 1);
+  for (size_t i = 0; i < SIZE; ++i) 
+  {
+    arr[i] = interval_begin + i * step;
+  }
+}
+
+// Check the file existence.
+inline bool file_exists(const std::string &filename) 
+{
+  std::ifstream f(filename.c_str());
+  return f.good();
+}
+
+#endif

TravelingWaves/CMakeLists.txt

@@ -0,0 +1,43 @@
+##
+#  CMake script for the TravelingWaves program:
+##
+
+# The name of the project and target:
+SET(TARGET "main")
+
+SET(TARGET_SRC
+  ${TARGET}.cc calculate_profile.cc Solution.cc TravelingWaveSolver.cc Parameters.cc LimitSolution.cc 
+  )
+
+CMAKE_MINIMUM_REQUIRED(VERSION 3.13.4)
+
+FIND_PACKAGE(deal.II 9.5.0
+  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+  )
+IF(NOT ${deal.II_FOUND})
+  MESSAGE(FATAL_ERROR "\n"
+    "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
+    "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
+    "or set an environment variable \"DEAL_II_DIR\" that contains this path."
+    )
+ENDIF()
+
+#
+# Are all dependencies fulfilled?
+#
+IF(NOT (DEAL_II_WITH_UMFPACK AND DEAL_II_WITH_SUNDIALS)) # keep in one line
+  MESSAGE(FATAL_ERROR "
+    Error! This program requires a deal.II library that was configured with the following options:
+        DEAL_II_WITH_UMFPACK  = ON
+        DEAL_II_WITH_SUNDIALS = ON
+    However, the deal.II library found at ${DEAL_II_PATH} was configured with these options:
+        DEAL_II_WITH_UMFPACK  = ${DEAL_II_WITH_UMFPACK}
+        DEAL_II_WITH_SUNDIALS = ${DEAL_II_WITH_SUNDIALS}
+    This conflicts with the requirements."
+    )
+ENDIF()
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+PROJECT(${TARGET})
+
+DEAL_II_INVOKE_AUTOPILOT()

TravelingWaves/IntegrateSystem.h

@@ -0,0 +1,103 @@
+
+#ifndef INTEGRATE_SYSTEM
+#define INTEGRATE_SYSTEM
+
+#include <boost/numeric/odeint.hpp>
+
+#include <iostream>
+#include <fstream>
+#include <string>
+
+template <typename state_T, typename time_T>
+void SaveSolutionIntoFile(const std::vector<state_T>& x_vec, const std::vector<time_T>& t_vec, std::string filename="output_ode_sol.txt") 
+{
+  if (!x_vec.empty() && !t_vec.empty()) 
+  {
+    std::ofstream output(filename);
+    output << std::setprecision(16);
+
+    size_t dim = x_vec[0].size();
+    for (size_t i = 0; i < t_vec.size(); ++i) 
+    {
+      output << std::fixed << t_vec[i];
+      for (size_t j = 0; j < dim; ++j) 
+      {
+        output << std::scientific << " " << x_vec[i][j];
+      }
+      output << "\n";
+    }
+    output.close();
+  } 
+  else 
+  {
+    std::cout << "Solution is not saved into file.\n";
+  }
+}
+
+// type of RK integrator 
+enum class Integrator_Type 
+{
+  dopri5,
+  cash_karp54,
+  fehlberg78
+};
+
+// Observer
+template <typename state_type>
+class Push_back_state_time
+{
+public:
+  std::vector<state_type>& m_states;
+  std::vector<double>& m_times;
+
+  Push_back_state_time(std::vector<state_type>& states, std::vector<double>& times)
+    : m_states(states), m_times(times) 
+  {}
+
+  void operator() (const state_type& x, double t) 
+  {
+    m_states.push_back(x);
+    m_times.push_back(t);
+  }
+};
+
+
+// Integrate system at specified points. 
+template <typename ODE_obj_T, typename state_type, typename Iterable_type>
+void IntegrateSystemAtTimePoints(
+  std::vector<state_type>& x_vec, std::vector<double>& t_vec, const Iterable_type& iterable_time_span,
+  const ODE_obj_T& ode_system_obj, 
+  state_type& x, const double dt,
+  Integrator_Type integrator_type=Integrator_Type::dopri5, bool save_to_file_flag=false,
+  const double abs_er_tol=1.0e-15, const double rel_er_tol=1.0e-15
+  )
+{
+  using namespace boost::numeric::odeint;
+
+  if (integrator_type == Integrator_Type::dopri5) 
+  {
+    typedef runge_kutta_dopri5< state_type > error_stepper_type;
+    integrate_times( make_controlled< error_stepper_type >(abs_er_tol, rel_er_tol),
+              ode_system_obj, x, iterable_time_span.begin(), iterable_time_span.end(), dt, Push_back_state_time< state_type >(x_vec, t_vec) );
+  } 
+  else if (integrator_type == Integrator_Type::cash_karp54) 
+  {
+    typedef runge_kutta_cash_karp54< state_type > error_stepper_type;
+    integrate_times( make_controlled< error_stepper_type >(abs_er_tol, rel_er_tol),
+              ode_system_obj, x, iterable_time_span.begin(), iterable_time_span.end(), dt, Push_back_state_time< state_type >(x_vec, t_vec) );
+  }
+  else 
+  { // Integrator_Type::fehlberg78
+    typedef runge_kutta_fehlberg78< state_type > error_stepper_type;
+    integrate_times( make_controlled< error_stepper_type >(abs_er_tol, rel_er_tol),
+              ode_system_obj, x, iterable_time_span.begin(), iterable_time_span.end(), dt, Push_back_state_time< state_type >(x_vec, t_vec) );
+  }
+
+  if (save_to_file_flag) 
+  {
+    SaveSolutionIntoFile(x_vec, t_vec);
+  } 
+
+}
+
+#endif

TravelingWaves/LimitSolution.cc

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

TravelingWaves/LimitSolution.h

@@ -0,0 +1,47 @@
+#ifndef LIMIT_SOLUTION
+#define LIMIT_SOLUTION
+
+#include "Parameters.h"
+#include <iostream>
+#include <vector>
+
+namespace TravelingWave
+{
+  typedef std::vector< double > state_type;
+
+  class LimitSolution
+  {
+  public:
+    LimitSolution(const Parameters &parameters, const double ilambda_0, const double iu_0, const double iT_0, const double root_sign = 1.);
+
+    void operator() (const state_type &x , state_type &dxdt , const double /* t */);
+    void calculate_u_T_omega();
+    void set_wave_speed(double iwave_speed);
+
+    std::vector<double> t_vec;
+    std::vector<state_type> omega_vec;
+    std::vector<state_type> lambda_vec;
+    std::vector<state_type> u_vec;
+    std::vector<state_type> T_vec;
+
+  private:
+    double omega_func(const double lambda, const double T) const;
+    double u_func(const double lambda) const;
+    double T_func(const double lambda) const;
+
+    void calculate_constants_A_B();
+
+    const Parameters &params;
+    const Problem    &problem;
+    double wave_speed;
+
+    const double lambda_0, u_0, T_0;  // Initial values.
+    double A, B;                      // Integration constants.
+
+    const double root_sign;           // Plus or minus one.
+  };
+
+
+} // namespace TravelingWave
+
+#endif

TravelingWaves/LinearInterpolator.h

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