-
#ifndef INTEGRATE_SYSTEM
#define INTEGRATE_SYSTEM
Problem::Problem()
: ParameterAcceptor("Problem")
{
- add_parameter("delta", delta = 0.01);
- add_parameter("epsilon", epsilon = 0.1);
- add_parameter("Prandtl number", Pr = 0.75);
- add_parameter("Lewis number", Le = 1.0);
- add_parameter("Constant of reaction rate", k = 1.0);
- add_parameter("Activation energy", theta = 1.65);
- add_parameter("Heat release", q = 1.7);
+ add_parameter("delta", delta = 0.01);
+ add_parameter("epsilon", epsilon = 0.1);
+ add_parameter("Prandtl number", Pr = 0.75);
+ add_parameter("Lewis number", Le = 1.0);
+ add_parameter("Constant of reaction rate", k = 1.0);
+ add_parameter("Activation energy", theta = 1.65);
+ add_parameter("Heat release", q = 1.7);
add_parameter("Ignition Temperature", T_ign = 1.0);
add_parameter("Type of the wave (deflagration / detonation)", wave_type = 1); // 0 for "deflagration"; 1 for "detonation".
}
Mesh::Mesh()
- : ParameterAcceptor("Mesh")
+ : ParameterAcceptor("Mesh")
{
- add_parameter("Interval left boundary", interval_left = -50.0);
- add_parameter("Interval right boundary", interval_right = 20.0);
- add_parameter<unsigned int>("Refinements number", refinements_number = 10);
- add_parameter("Adaptive mesh refinement", adaptive = 1); // 1 for adaptive; 0 for global.
+ add_parameter("Interval left boundary", interval_left = -50.0);
+ add_parameter("Interval right boundary", interval_right = 20.0);
+ add_parameter<unsigned int>("Refinements number", refinements_number = 10);
+ add_parameter("Adaptive mesh refinement", adaptive = 1); // 1 for adaptive; 0 for global.
}
Solver::Solver()
struct FiniteElements : ParameterAcceptor
{
- FiniteElements();
+ FiniteElements();
unsigned int poly_degree;
unsigned int quadrature_points_number;
{
Mesh();
- double interval_left;
- double interval_right;
- unsigned int refinements_number;
- int adaptive;
+ double interval_left;
+ double interval_right;
+ unsigned int refinements_number;
+ int adaptive;
};
struct Solver : ParameterAcceptor
{
Solver();
- double tol;
+ double tol;
};
struct Parameters
{
Problem problem;
- FiniteElements fe;
- Mesh mesh;
- Solver solver;
+ FiniteElements fe;
+ Mesh mesh;
+ Solver solver;
};
} // namespace TravelingWave
SolutionStruct::SolutionStruct() {}
SolutionStruct::SolutionStruct(const std::vector<double> &ix, const std::vector<double> &iu,
- const std::vector<double> &iT, const std::vector<double> &ilambda, double iwave_speed)
+ const std::vector<double> &iT, const std::vector<double> &ilambda, double iwave_speed)
: x(ix)
, u(iu)
, T(iT)
, wave_speed(iwave_speed)
{}
SolutionStruct::SolutionStruct(const std::vector<double> &ix, const std::vector<double> &iu,
- const std::vector<double> &iT, const std::vector<double> &ilambda)
+ const std::vector<double> &iT, const std::vector<double> &ilambda)
: SolutionStruct(ix, iu, iT, ilambda, 0.)
{}
double Interpolant::value(const Point<1> &p, const unsigned int component) const
{
- double x = p[0];
+ double x = p[0];
double res = interpolant.value(x);
return res;
double SolutionVectorFunction<InterpolantType>::value(const Point<1> &p, const unsigned int component) const
{
double res = 0.;
- if (component == 0) { res = u_interpolant.value(p); }
- else if (component == 1) { res = T_interpolant.value(p); }
- else if (component == 2) { res = lambda_interpolant.value(p); }
+ if (component == 0) { res = u_interpolant.value(p); }
+ else if (component == 1) { res = T_interpolant.value(p); }
+ else if (component == 2) { res = lambda_interpolant.value(p); }
return res;
}
-#ifndef SOLUTION_STRUCT
-#define SOLUTION_STRUCT
+#ifndef SOLUTION
+#define SOLUTION
#include <deal.II/base/function.h>
{
SolutionStruct();
SolutionStruct(const std::vector<double> &ix, const std::vector<double> &iu,
- const std::vector<double> &iT, const std::vector<double> &ilambda, const double iwave_speed);
+ const std::vector<double> &iT, const std::vector<double> &ilambda, const double iwave_speed);
SolutionStruct(const std::vector<double> &ix, const std::vector<double> &iu,
- const std::vector<double> &iT, const std::vector<double> &ilambda);
+ const std::vector<double> &iT, const std::vector<double> &ilambda);
void reinit(const unsigned int number_of_elements);
void save_to_file(std::string filename) const;
- std::vector<double> x; // mesh coordinates (must be an increasing sequence)
- std::vector<double> u; // array of u components
- std::vector<double> T; // array of T components
- std::vector<double> lambda; // array of lambda components
+ std::vector<double> x; // mesh coordinates (must be an increasing sequence)
+ std::vector<double> u; // array of u components
+ std::vector<double> T; // array of T components
+ std::vector<double> lambda; // array of lambda components
- double wave_speed; // speed of the wave
+ double wave_speed; // speed of the wave
};
// Interpolation class
, triangulation_uploaded(false)
, fe(FE_Q<1>(params.fe.poly_degree), 1,
FE_Q<1>(params.fe.poly_degree), 1,
- FE_Q<1>(params.fe.poly_degree), 1) // 3 fe basis sets, corresponding to du, dT, dlambda
+ FE_Q<1>(params.fe.poly_degree), 1) // 3 fe basis sets, corresponding to du, dT, dlambda
, dof_handler(triangulation)
, current_wave_speed(0.)
, initial_guess(initial_guess_input)
-#ifndef WAVE_CONSTRUCTOR
-#define WAVE_CONSTRUCTOR
+#ifndef TRAVELING_WAVE_SOLVER
+#define TRAVELING_WAVE_SOLVER
#include <deal.II/base/timer.h>
#include <deal.II/base/function.h>
{
public:
TravelingWaveSolver(const Parameters ¶meters, const SolutionStruct &initial_guess_input);
-
+
void set_triangulation(const Triangulation<1> &itriangulation);
-
+
void run(const std::string filename="solution", const bool save_solution_to_file=true);
void get_solution(SolutionStruct &solution) const;
void get_triangulation(Triangulation<1> &otriangulation) const;
void set_initial_guess();
double Heaviside_func(double x) const;
-
+
void compute_and_factorize_jacobian(const Vector<double> &evaluation_point_extended);
double compute_residual(const Vector<double> &evaluation_point_extended, Vector<double> &residual);
void split_extended_solution_vector();
std::map<std::string, unsigned int> boundary_and_centering_dof_numbers;
// Parameters of the problem, taken from a .prm file.
- const Parameters ¶ms;
- const Problem &problem; // Reference variable, just for convenience.
+ const Parameters ¶ms;
+ const Problem &problem; // Reference variable, just for convenience.
unsigned int number_of_quadrature_points;
Triangulation<1> triangulation;
// The flag indicating whether the triangulation was uploaded externally or created within the <code> run </code> member function.
- bool triangulation_uploaded;
- FESystem<1> fe;
- DoFHandler<1> dof_handler;
+ bool triangulation_uploaded;
+ FESystem<1> fe;
+ DoFHandler<1> dof_handler;
// Constraints for Dirichlet boundary conditions.
AffineConstraints<double> zero_boundary_constraints;
- SparsityPattern sparsity_pattern_extended;
- SparseMatrix<double> jacobian_matrix_extended;
- std::unique_ptr<SparseDirectUMFPACK> jacobian_matrix_extended_factorization;
+ SparsityPattern sparsity_pattern_extended;
+ SparseMatrix<double> jacobian_matrix_extended;
+ std::unique_ptr<SparseDirectUMFPACK> jacobian_matrix_extended_factorization;
// Finite element solution of the problem.
- Vector<double> current_solution;
+ Vector<double> current_solution;
// Value of the wave speed $c$.
- double current_wave_speed;
+ double current_wave_speed;
// Solution with an additional term, corresponding to the variable wave_speed.
- Vector<double> current_solution_extended;
+ Vector<double> current_solution_extended;
// Initial guess for Newton's iterations.
- SolutionStruct initial_guess;
+ SolutionStruct initial_guess;
- TimerOutput computing_timer;
+ TimerOutput computing_timer;
};
} // namespace TravelingWave
// Error estimation.
{
unsigned int sol_length = sol.x.size();
- double u_r = sol.u[sol_length-1]; // Dirichlet boundary condition
- double T_r = sol.T[sol_length-1]; // Dirichlet condition only for detonation case
+ double u_r = sol.u[sol_length-1]; // Dirichlet boundary condition
+ double T_r = sol.T[sol_length-1]; // Dirichlet condition only for detonation case
double u_l = sol.u[0];
- double T_l = sol.T[0]; // Dirichlet boundary condition
+ double T_l = sol.T[0]; // Dirichlet boundary condition
double wave_speed = sol.wave_speed;
std::cout << "Error estimates:" << std::endl;
-#ifndef INITIAL_GUESS
-#define INITIAL_GUESS
+#ifndef CALCULATE_PROFILE
+#define CALCULATE_PROFILE
#include "Parameters.h"
#include "Solution.h"
void compute_initial_guess_deflagration(const Parameters ¶ms, SolutionStruct &initial_guess);
void calculate_profile(Parameters& parameters,
- const bool continuation_for_delta=false /* Compute with the continuation. */,
- const double delta_start=0.01 /* The starting value of delta for the continuation method. */,
- const unsigned int number_of_continuation_points=10);
+ const bool continuation_for_delta=false /* Compute with the continuation. */,
+ const double delta_start=0.01 /* The starting value of delta for the continuation method. */,
+ const unsigned int number_of_continuation_points=10);
} // namespace TravelingWave