388 *
template <
typename T>
389 *
bool isapprox(
const T &a,
const T &b,
const double tol = 1e-10)
399 *
template <
typename T>
402 *
const size_t SIZE =
arr.size();
404 *
for (
size_t i = 0; i <
SIZE; ++i)
417 *
std::ifstream f(
filename.c_str());
425<a name=
"ann-IntegrateSystem.h"></a>
451 *
template <
typename state_T,
typename time_T>
457 *
output << std::setprecision(16);
459 *
size_t dim =
x_vec[0].size();
460 *
for (
size_t i = 0; i <
t_vec.size(); ++i)
462 *
output << std::fixed <<
t_vec[i];
463 *
for (
size_t j = 0;
j < dim; ++
j)
465 *
output << std::scientific <<
" " <<
x_vec[i][
j];
473 *
std::cout <<
"Solution is not saved into file.\n";
494 *
template <
typename state_type>
498 *
std::vector<state_type>&
m_states;
499 *
std::vector<double>&
m_times;
518 *
template <
typename ODE_obj_T,
typename state_type,
typename Iterable_type>
527 *
using namespace boost::numeric::odeint;
559<a name=
"ann-LimitSolution.cc"></a>
633 *
for (
unsigned int i = 0; i <
lambda_vec.size(); ++i)
635 *
u_vec[i].resize(1);
636 *
T_vec[i].resize(1);
646 *
std::cout <<
"t_vec or lambda_vec vector is empty!" << std::endl;
654<a name=
"ann-LimitSolution.h"></a>
691 *
std::vector<double>
t_vec;
694 *
std::vector<state_type>
u_vec;
695 *
std::vector<state_type>
T_vec;
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;
721<a name=
"ann-LinearInterpolator.h"></a>
750 *
template <
typename Number_Type>
758 *
const std::vector<Number_Type>
x_points;
759 *
const std::vector<Number_Type>
y_points;
762 *
template <
typename Number_Type>
768 *
template <
typename Number_Type>
804<a name=
"ann-Parameters.cc"></a>
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);
840 *
add_parameter(
"Type of boundary condition for the temperature at the right boundary",
T_r_bc_type = 1);
842 *
add_parameter(
"T_left",
T_left = 5.3);
843 *
add_parameter(
"T_right",
T_right = 0.9);
844 *
add_parameter(
"u_left",
u_left = -0.2);
845 *
add_parameter(
"u_right",
u_right = 0.);
847 *
add_parameter(
"Initial guess for the wave speed",
wave_speed_init = 1.2);
853 *
add_parameter(
"Polynomial degree",
poly_degree = 1);
860 *
add_parameter(
"Interval left boundary",
interval_left = -50.0);
861 *
add_parameter(
"Interval right boundary",
interval_right = 20.0);
863 *
add_parameter(
"Adaptive mesh refinement", adaptive = 1);
869 *
add_parameter(
"Tolerance", tol = 1e-10);
876<a name=
"ann-Parameters.h"></a>
908 *
double k, theta,
q;
957<a name=
"ann-Solution.cc"></a>
991 *
const std::vector<double> &
iT,
const std::vector<double> &
ilambda)
1014 *
output << std::scientific << std::setprecision(16);
1015 *
for (
unsigned int i = 0; i <
x.size(); ++i)
1017 *
output << std::fixed <<
x[i];
1018 *
output << std::scientific <<
" " <<
u[i] <<
" " <<
T[i] <<
" " <<
lambda[i] <<
"\n";
1033 *
double Interpolant::value(
const Point<1> &p,
const unsigned int component)
const
1042 *
template <
typename InterpolantType>
1050 *
template <
typename InterpolantType>
1067<a name=
"ann-Solution.h"></a>
1089 *
#
include "LinearInterpolator.h"
1100 *
using namespace dealii;
1111 *
const std::vector<double> &
iT,
const std::vector<double> &
ilambda,
const double iwave_speed);
1113 *
const std::vector<double> &
iT,
const std::vector<double> &
ilambda);
1119 *
std::vector<double>
x;
1120 *
std::vector<double>
u;
1121 *
std::vector<double>
T;
1122 *
std::vector<double>
lambda;
1136 *
virtual double value(
const Point<1> &p,
const unsigned int component = 0)
const override;
1147 *
template <
typename InterpolantType>
1152 *
virtual double value(
const Point<1> &p,
const unsigned int component = 0)
const override;
1166<a name=
"ann-TravelingWaveSolver.cc"></a>
1183 *
#
include "TravelingWaveSolver.h"
1187 *
using namespace dealii;
1213 *
table.add_value(
"Parameter name",
"number of quadrature points");
1216 *
table.add_value(
"Parameter name",
"delta");
1217 *
table.add_value(
"value",
params.problem.delta);
1219 *
table.add_value(
"Parameter name",
"epsilon");
1220 *
table.add_value(
"value",
params.problem.epsilon);
1222 *
table.add_value(
"Parameter name",
"params.problem.wave_speed_init");
1223 *
table.add_value(
"value",
params.problem.wave_speed_init);
1225 *
table.add_value(
"Parameter name",
"initial_guess.wave_speed");
1228 *
table.add_value(
"Parameter name",
"T_left");
1229 *
table.add_value(
"value",
params.problem.T_left);
1231 *
table.set_precision(
"value", 2);
1232 *
table.set_scientific(
"value",
true);
1234 *
std::cout <<
"\n";
1236 *
std::cout <<
"\n";
1258 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1292 *
if (
problem.T_r_bc_type == 1)
1306 *
dof_handler.distribute_dofs(fe);
1308 *
std::cout <<
"Number of dofs : " << dof_handler.n_dofs() << std::endl;
1339 *
if (
problem.T_r_bc_type == 1)
1341 *
std::cout <<
"Dirichlet condition for the temperature at the right boundary." << std::endl;
1346 *
std::cout <<
"Neumann condition for the temperature at the right boundary." << std::endl;
1366 *
dofs_on_this_cell.reserve(dof_handler.get_fe_collection().max_dofs_per_cell());
1368 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1370 *
const unsigned
int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1465 *
for (
unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
1472 *
std::cout <<
"Computing Jacobian matrix ... " << std::endl;
1483 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1489 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
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);
1510 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1515 *
fe_values.reinit(cell);
1531 *
auto kappa_2 = [=](
double T,
double ){
1535 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
1537 *
for (
unsigned int k = 0;
k < dofs_per_cell; ++
k)
1550 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1552 *
for (
unsigned int j = 0;
j < dofs_per_cell; ++
j)
1578 *
) * fe_values.JxW(
q);
1587 *
) * fe_values.JxW(
q);
1592 *
cell->get_dof_indices(local_dof_indices);
1594 *
for (
const unsigned int i : fe_values.dof_indices())
1599 *
local_dof_indices[
j],
1648 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1657 *
fe_values.reinit(cell);
1707 *
std::cout <<
"Factorizing Jacobian matrix" << std::endl;
1720 *
std::cout <<
"Computing residual vector ... " << std::endl;
1725 *
for (
unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
1738 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1742 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
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);
1762 *
for (
const auto &cell : dof_handler.active_cell_iterators())
1766 *
fe_values.reinit(cell);
1775 *
auto omega = [=](
double T,
double lambda){
1779 *
for (
unsigned int q = 0;
q < n_q_points; ++
q)
1781 *
for (
unsigned int k = 0;
k < dofs_per_cell; ++
k)
1794 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
1815 *
) * fe_values.JxW(
q);
1820 *
cell->get_dof_indices(local_dof_indices);
1822 *
for (
const unsigned int i : fe_values.dof_indices())
1834 *
std::cout << std::defaultfloat;
1835 *
std::cout <<
"norm of residual = " <<
residual_norm << std::endl;
1860 *
std::cout <<
"Solving linear system ... " << std::endl;
1886 *
fe.component_mask(lambda)
1950 *
this->solve(
rhs, solution, tolerance);
1974 *
for (
const auto &cell : dof_handler.active_cell_iterators())
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));
1982 *
output << std::scientific << std::setprecision(16);
1983 *
output << cell->vertex(
v_ind)[0];
1985 *
output << std::scientific << std::setprecision(16);
1986 *
output << std::scientific <<
" " <<
u <<
" " <<
T <<
" " <<
lambda <<
"\n";
2011 *
auto comp = [](
const std::vector<double> &a,
const std::vector<double> &
b) {
2012 *
return a[0] <
b[0];
2015 *
for (
const auto &cell : dof_handler.active_cell_iterators())
2019 *
double
x = cell->vertex(
v_ind)[0];
2027 *
solution.x.
clear();
2028 *
solution.u.
clear();
2029 *
solution.T.
clear();
2030 *
solution.lambda.
clear();
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]);
2061 *
const double tol =
params.solver.tol;
2073 *
static_cast<unsigned int>(
std::abs( 0. -
params.mesh.interval_left )),
2074 *
params.mesh.interval_left, 0.
2080 *
static_cast<unsigned int>(
std::abs(
params.mesh.interval_right - 0. )),
2081 *
0.,
params.mesh.interval_right
2131 *
std::cout << std::scientific << std::setprecision(16);
2133 *
std::cout << std::defaultfloat;
2157 *
std::cout << std::scientific << std::setprecision(16);
2159 *
std::cout << std::defaultfloat;
2173<a name=
"ann-TravelingWaveSolver.h"></a>
2218 *
#
include "AuxiliaryFunctions.h"
2235 *
using namespace dealii;
2337 * SolutionStruct initial_guess;
2339 * TimerOutput computing_timer;
2342 * } // namespace TravelingWave
2348<a name="ann-calculate_profile.cc"></a>
2349<h1>Annotated version of calculate_profile.cc</h1>
2355 * /* -----------------------------------------------------------------------------
2357 * * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
2358 * * Copyright (C) 2024 by Shamil Magomedov
2360 * * This file is part of the deal.II code gallery.
2362 * * -----------------------------------------------------------------------------
2365 * #include "TravelingWaveSolver.h"
2366 * #include "calculate_profile.h"
2368 * namespace TravelingWave
2370 * using namespace dealii;
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.
2377 * void compute_limit_sol_left_part(const Parameters ¶meters,
2378 * const double wave_speed,
2381 * const double lambda_0,
2382 * SolutionStruct &LimitSol,
2383 * const double root_sign)
2385 * LimitSolution limit_sol(parameters, lambda_0, u_0, T_0, root_sign);
2386 * limit_sol.set_wave_speed(wave_speed);
2391 * We take more integration points to better resolve the transition layer.
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());
2403 * Reverse the order of the elements (because we need to perform back in time integration).
2406 * std::reverse(t_span.begin(), t_span.end());
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,
2413 * -1e-2, Integrator_Type::dopri5);
2416 * limit_sol.calculate_u_T_omega();
2420 * Reverse the order of elements
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());
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");
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)
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];
2448 * Construction of an initial guess for detonation wave solution. The ODE is solved for the ideal system with @f$\delta = 0@f$.
2451 * void compute_initial_guess_detonation(const Parameters ¶ms, SolutionStruct &initial_guess, const double root_sign)
2453 * const Problem &problem = params.problem;
2454 * double current_wave_speed(problem.wave_speed_init);
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));
2462 * double u_0 = problem.u_right;
2463 * double T_0 = problem.T_right;
2464 * double lambda_0 = 0.;
2466 * compute_limit_sol_left_part(params, current_wave_speed, u_0, T_0, lambda_0, initial_guess, root_sign);
2468 * initial_guess.wave_speed = current_wave_speed;
2470 * for (int i = initial_guess.x.size() - 1; i > - 1; --i)
2472 * if (isapprox(initial_guess.x[i], 0.))
2474 * initial_guess.u[i] = problem.u_right;
2475 * initial_guess.T[i] = problem.T_ign;
2476 * initial_guess.lambda[i] = 0.;
2483 * Adding the points to the right part of the interval (w.r.t. @f$\xi = 0@f$).
2486 * unsigned int number_of_additional_points = 5;
2487 * for (unsigned int i = 0; i < number_of_additional_points; ++i)
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.);
2500 * Construction of a piecewise constant initial guess for deflagration wave solution.
2503 * void compute_initial_guess_deflagration(const Parameters ¶ms, SolutionStruct &initial_guess)
2505 * const Problem &problem = params.problem;
2506 * double current_wave_speed(problem.wave_speed_init);
2508 * double del_Pr_eps = (problem.Pr * 4 * problem.delta / (3 * problem.epsilon));
2509 * double del_Le = (problem.delta / problem.Le);
2511 * auto u_init_guess_func = [&](double x) {
2514 * return problem.u_left;
2518 * return problem.u_right;
2522 * auto T_init_guess_func = [&](double x) {
2525 * return problem.T_left;
2527 * else if (isapprox(x, 0.))
2529 * return problem.T_ign;
2533 * return problem.T_right;
2537 * auto lambda_init_guess_func = [=](double x) {
2540 * return -std::exp(x * std::abs(1 - current_wave_speed) / del_Pr_eps) + 1;
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);
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);
2557 * for (unsigned int i = 0; i < number_of_points; ++i)
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]);
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;
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>.
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)
2583 * SolutionStruct sol;
2585 * if (parameters.problem.wave_type == 1) // detonation wave
2587 * compute_initial_guess_detonation(parameters, sol);
2589 * else if (parameters.problem.wave_type == 0) // deflagration wave
2591 * compute_initial_guess_deflagration(parameters, sol);
2594 * if (continuation_for_delta == false)
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);
2602 * else // Run with continuation_for_delta.
2604 * double delta_target = parameters.problem.delta;
2605 * parameters.problem.delta = delta_start;
2607 * std::vector<double> delta_span(number_of_continuation_points);
2611 * Generate a sequence of delta values being uniformly distributed in log10 scale.
2615 * double delta_start_log10 = std::log10(delta_start);
2616 * double delta_target_log10 = std::log10(delta_target);
2618 * std::vector<double> delta_log_span(delta_span.size());
2619 * linspace(delta_start_log10, delta_target_log10, delta_log_span);
2621 * for (unsigned int i = 0; i < delta_span.size(); ++i)
2623 * delta_span[i] = std::pow(10, delta_log_span[i]);
2627 * Triangulation<1> refined_triangulation;
2628 * bool first_iter_flag = true;
2630 * for (double delta : delta_span)
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);
2636 * TravelingWaveSolver wave(parameters, sol);
2638 * if (first_iter_flag)
2640 * first_iter_flag = false;
2644 * wave.set_triangulation(refined_triangulation);
2647 * wave.run(filename);
2648 * wave.get_solution(sol);
2649 * wave.get_triangulation(refined_triangulation);
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;
2667 * std::cout << "Error estimates:" << std::endl;
2668 * double DeltaT = T_l - T_r;
2669 * double qDT = parameters.problem.q - DeltaT;
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;
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;
2680 * } // namespace TravelingWave
2684<a name="ann-calculate_profile.h"></a>
2685<h1>Annotated version of calculate_profile.h</h1>
2691 * /* -----------------------------------------------------------------------------
2693 * * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
2694 * * Copyright (C) 2024 by Shamil Magomedov
2696 * * This file is part of the deal.II code gallery.
2698 * * -----------------------------------------------------------------------------
2701 * #ifndef CALCULATE_PROFILE
2702 * #define CALCULATE_PROFILE
2704 * #include "Parameters.h"
2705 * #include "Solution.h"
2706 * #include "LimitSolution.h"
2707 * #include "IntegrateSystem.h"
2708 * #include "AuxiliaryFunctions.h"
2710 * namespace TravelingWave
2712 * void compute_limit_sol_left_part(const Parameters ¶meters,
2713 * const double wave_speed,
2716 * const double lambda_0,
2717 * SolutionStruct &LimitSol,
2718 * const double root_sign = 1.);
2720 * void compute_initial_guess_detonation(const Parameters ¶ms, SolutionStruct &initial_guess, const double root_sign = 1.);
2721 * void compute_initial_guess_deflagration(const Parameters ¶ms, SolutionStruct &initial_guess);
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);
2728 * } // namespace TravelingWave
2734<a name="ann-main.cc"></a>
2735<h1>Annotated version of main.cc</h1>
2741 * /* -----------------------------------------------------------------------------
2743 * * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
2744 * * Copyright (C) 2024 by Shamil Magomedov
2746 * * This file is part of the deal.II code gallery.
2748 * * -----------------------------------------------------------------------------
2751 * #include "calculate_profile.h"
2753 * int main(int argc, char *argv[])
2758 * using namespace TravelingWave;
2760 * Parameters parameters;
2762 * std::string prm_filename = "ParametersList.prm";
2767 * Check if file argv[1] exists.
2770 * if (file_exists(argv[1]))
2772 * prm_filename = argv[1];
2776 * std::string errorMessage = "File \"" + std::string(argv[1]) + "\" is not found.";
2777 * throw std::runtime_error(errorMessage);
2784 * Check if the file "ParametersList.prm" exists in the current or in the parent directory.
2787 * if (!(file_exists(prm_filename) || file_exists("../" + prm_filename)))
2789 * std::string errorMessage = "File \"" + prm_filename + "\" is not found.";
2790 * throw std::runtime_error(errorMessage);
2794 * if (!file_exists(prm_filename))
2796 * prm_filename = "../" + prm_filename;
2801 * std::cout << "Reading parameters... " << std::flush;
2802 * ParameterAcceptor::initialize(prm_filename);
2803 * std::cout << "done" << std::endl;
2805 * calculate_profile(parameters, /* With continuation_for_delta */ false, 0.1, 3);
2808 * catch (std::exception &exc)
2810 * std::cerr << std::endl
2812 * << "----------------------------------------------------"
2814 * std::cerr << "Exception on processing: " << std::endl
2815 * << exc.what() << std::endl
2816 * << "Aborting!" << std::endl
2817 * << "----------------------------------------------------"
2823 * std::cerr << std::endl
2825 * << "----------------------------------------------------"
2827 * std::cerr << "Unknown exception!" << std::endl
2828 * << "Aborting!" << std::endl
2829 * << "----------------------------------------------------"
2839<a name="ann-plot.py"></a>
2840<h1>Annotated version of plot.py</h1>
2842## -----------------------------------------------------------------------------
2844## SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
2845## Copyright (C) 2024 by Shamil Magomedov
2847## This file is part of the deal.II code gallery.
2849## -----------------------------------------------------------------------------
2853import matplotlib.pyplot as plt
2859 # 'lines.marker
' : 'x',
2861 'lines.markersize
' : 4,
2862 'lines.linewidth
' : 1,
2863 'axes.labelsize
': 16,
2867 'xtick.labelsize
': 14,
2868 'ytick.labelsize
': 14,
2869 'text.usetex
': True,
2874plt.rcParams.update(plot_params)
2877if len(sys.argv) > 1:
2879 filename = sys.argv[1]
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]
2890 fig, ax = plt.subplots(nrows=1, ncols=1)
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
')
2897 # Plot of limit solutions for the detonation case. Uncomment, if needed.
2898 #===============================================================#
2929 #===============================================================#
2932 ax.set_xlabel(r"@f$\xi@f$")
2933 ax.set_ylabel(r"@f$u, T, \lambda@f$")
2936 # plt.savefig("fast_deflagration_delta_0.01.png", bbox_inches='tight', dpi=500)
2942 print("No such file:", filename)
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)
static constexpr unsigned int dimension
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
#define DEAL_II_VERSION_GTE(major, minor, subminor)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
void subdivided_hyper_cube(Triangulation< dim, spacedim > &tria, const unsigned int repetitions, const double left=0., const double right=1., const bool colorize=false)
void merge_triangulations(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result, const double duplicated_vertex_tolerance=1.0e-12, const bool copy_manifold_ids=false, const bool copy_boundary_ids=false)
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
@ matrix
Contents is actually a matrix.
constexpr types::blas_int zero
constexpr types::blas_int one
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.)
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.)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
T scatter(const MPI_Comm comm, const std::vector< T > &objects_to_send, const unsigned int root_process=0)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
DEAL_II_HOST constexpr TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation