Skip to content

Commit

Permalink
Merge pull request #190 from magomedovs/minor_fixes
Browse files Browse the repository at this point in the history
Minor code fixes
  • Loading branch information
bangerth authored Jul 11, 2024
2 parents 2cfb27c + 0122f2b commit e5ace12
Show file tree
Hide file tree
Showing 9 changed files with 63 additions and 64 deletions.
1 change: 0 additions & 1 deletion TravelingWaves/IntegrateSystem.h
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@

#ifndef INTEGRATE_SYSTEM
#define INTEGRATE_SYSTEM

Expand Down
24 changes: 12 additions & 12 deletions TravelingWaves/Parameters.cc
Original file line number Diff line number Diff line change
Expand Up @@ -7,13 +7,13 @@ namespace TravelingWave
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".

Expand All @@ -35,12 +35,12 @@ namespace TravelingWave
}

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()
Expand Down
18 changes: 9 additions & 9 deletions TravelingWaves/Parameters.h
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ namespace TravelingWave

struct FiniteElements : ParameterAcceptor
{
FiniteElements();
FiniteElements();

unsigned int poly_degree;
unsigned int quadrature_points_number;
Expand All @@ -35,25 +35,25 @@ namespace TravelingWave
{
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
Expand Down
12 changes: 6 additions & 6 deletions TravelingWaves/Solution.cc
Original file line number Diff line number Diff line change
Expand Up @@ -7,15 +7,15 @@ 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)
, lambda(ilambda)
, wave_speed(iwave_speed)
{}
SolutionStruct::SolutionStruct(const std::vector<double> &ix, const std::vector<double> &iu,
const std::vector<double> &iT, const std::vector<double> &ilambda)
const std::vector<double> &iT, const std::vector<double> &ilambda)
: SolutionStruct(ix, iu, iT, ilambda, 0.)
{}

Expand Down Expand Up @@ -59,7 +59,7 @@ namespace TravelingWave

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;
Expand All @@ -78,9 +78,9 @@ namespace TravelingWave
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;
}
Expand Down
18 changes: 9 additions & 9 deletions TravelingWaves/Solution.h
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
#ifndef SOLUTION_STRUCT
#define SOLUTION_STRUCT
#ifndef SOLUTION
#define SOLUTION

#include <deal.II/base/function.h>

Expand All @@ -21,20 +21,20 @@ namespace TravelingWave
{
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
Expand Down
2 changes: 1 addition & 1 deletion TravelingWaves/TravelingWaveSolver.cc
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ namespace TravelingWave
, 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)
Expand Down
36 changes: 18 additions & 18 deletions TravelingWaves/TravelingWaveSolver.h
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
#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>
Expand Down Expand Up @@ -46,9 +46,9 @@ namespace TravelingWave
{
public:
TravelingWaveSolver(const Parameters &parameters, const SolutionStruct &initial_guess_input);

void set_triangulation(const Triangulation<1> &itriangulation);

void run(const std::string filename="solution", const bool save_solution_to_file=true);
void get_solution(SolutionStruct &solution) const;
void get_triangulation(Triangulation<1> &otriangulation) const;
Expand All @@ -61,7 +61,7 @@ namespace TravelingWave
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();
Expand All @@ -77,37 +77,37 @@ namespace TravelingWave
std::map<std::string, unsigned int> boundary_and_centering_dof_numbers;

// Parameters of the problem, taken from a .prm file.
const Parameters &params;
const Problem &problem; // Reference variable, just for convenience.
const Parameters &params;
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
Expand Down
6 changes: 3 additions & 3 deletions TravelingWaves/calculate_profile.cc
Original file line number Diff line number Diff line change
Expand Up @@ -254,10 +254,10 @@ 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;
Expand Down
10 changes: 5 additions & 5 deletions TravelingWaves/calculate_profile.h
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
#ifndef INITIAL_GUESS
#define INITIAL_GUESS
#ifndef CALCULATE_PROFILE
#define CALCULATE_PROFILE

#include "Parameters.h"
#include "Solution.h"
Expand All @@ -21,9 +21,9 @@ namespace TravelingWave
void compute_initial_guess_deflagration(const Parameters &params, SolutionStruct &initial_guess);

void calculate_profile(Parameters& parameters,
const bool continuation_for_delta=false /* Compute with the continuation. */,
const double delta_start=0.01 /* The starting value of delta for the continuation method. */,
const unsigned int number_of_continuation_points=10);
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

Expand Down

0 comments on commit e5ace12

Please sign in to comment.