922 * <a name=
"step_67-Equationdata"></a>
943 *
const unsigned int component = 0)
const override;
975 *
const unsigned int component)
const
977 *
const double t = this->
get_time();
983 *
Assert(dim == 2, ExcNotImplemented());
984 *
const double beta = 5;
990 *
const double factor =
993 *
std::abs(1. - (gamma - 1.) / gamma * 0.25 * factor * factor));
995 *
const double u = 1. - factor * (
x[1] -
x0[1]);
996 *
const double v = factor * (
x[0] - t -
x0[0]);
998 *
if (component == 0)
1000 *
else if (component == 1)
1002 *
else if (component == 2)
1007 *
std::exp2(
density_log * (gamma / (gamma - 1.)));
1015 *
if (component == 0)
1017 *
else if (component == 1)
1019 *
else if (component == dim + 1)
1020 *
return 3.097857142857143;
1036 * <a name=
"step_67-LowstorageexplicitRungeKuttatimeintegrators"></a>
1156 *
unsigned int n_stages()
const
1185 *
template <
typename VectorType,
typename Operator>
1187 *
const double current_time,
1189 *
VectorType &solution,
1191 *
VectorType &
vec_ki)
const
1219 *
std::vector<double>
bi;
1220 *
std::vector<double>
ai;
1221 *
std::vector<double>
ci;
1229 * <a name=
"step_67-ImplementationofpointwiseoperationsoftheEulerequations"></a>
1248 *
href=
"https://en.wikipedia.org/wiki/Inline_function">
inline</a>
to where
1257 * context in code that comes after the place where the function was called.
1258 * (We note that compilers are generally quite good at figuring out which
1259 * functions to inline by themselves. Here is a place where compilers may or
1260 * may not have figured it out by themselves but where we know for sure that
1261 * inlining is a win.)
1265 * Another trick we apply is a separate variable for the inverse density
1266 * @f$\frac{1}{\rho}@f$. This enables the compiler to only perform a single
1267 * division for the flux, despite the division being used at several
1268 * places. As divisions are around ten to twenty times as expensive as
1269 * multiplications or additions, avoiding redundant divisions is crucial for
1270 * performance. We note that taking the inverse first and later multiplying
1271 * with it is not equivalent to a division in floating point arithmetic due
1272 * to roundoff effects, so the compiler is not allowed to exchange one way by
1273 * the other with standard optimization flags. However, it is also not
1274 * particularly difficult to write the code in the right way.
1278 * To summarize, the chosen strategy of always inlining and careful
1279 * definition of expensive arithmetic operations allows us to write compact
1280 * code without passing all intermediate results around, despite making sure
1281 * that the code maps to excellent machine code.
1284 * template <int dim, typename Number>
1285 * inline DEAL_II_ALWAYS_INLINE
1286 * Tensor<1, dim, Number>
1287 * euler_velocity(const Tensor<1, dim + 2, Number> &conserved_variables)
1289 * const Number inverse_density = Number(1.) / conserved_variables[0];
1291 * Tensor<1, dim, Number> velocity;
1292 * for (unsigned int d = 0; d < dim; ++d)
1293 * velocity[d] = conserved_variables[1 + d] * inverse_density;
1300 * The next function computes the pressure from the vector of conserved
1301 * variables, using the formula @f$p = (\gamma - 1) \left(E - \frac 12 \rho
1302 * \mathbf{u}\cdot \mathbf{u}\right)@f$. As explained above, we use the
1303 * velocity from the `euler_velocity()` function. Note that we need to
1304 * specify the first template argument `dim` here because the compiler is
1305 * not able to deduce it from the arguments of the tensor, whereas the
1306 * second argument (number type) can be automatically deduced.
1309 * template <int dim, typename Number>
1310 * inline DEAL_II_ALWAYS_INLINE
1312 * euler_pressure(const Tensor<1, dim + 2, Number> &conserved_variables)
1314 * const Tensor<1, dim, Number> velocity =
1315 * euler_velocity<dim>(conserved_variables);
1317 * Number rho_u_dot_u = conserved_variables[1] * velocity[0];
1318 * for (unsigned int d = 1; d < dim; ++d)
1319 * rho_u_dot_u += conserved_variables[1 + d] * velocity[d];
1321 * return (gamma - 1.) * (conserved_variables[dim + 1] - 0.5 * rho_u_dot_u);
1326 * Here is the definition of the Euler flux function, i.e., the definition
1327 * of the actual equation. Given the velocity and pressure (that the
1328 * compiler optimization will make sure are done only once), this is
1329 * straight-forward given the equation stated in the introduction.
1332 * template <int dim, typename Number>
1333 * inline DEAL_II_ALWAYS_INLINE
1334 * Tensor<1, dim + 2, Tensor<1, dim, Number>>
1335 * euler_flux(const Tensor<1, dim + 2, Number> &conserved_variables)
1337 * const Tensor<1, dim, Number> velocity =
1338 * euler_velocity<dim>(conserved_variables);
1339 * const Number pressure = euler_pressure<dim>(conserved_variables);
1341 * Tensor<1, dim + 2, Tensor<1, dim, Number>> flux;
1342 * for (unsigned int d = 0; d < dim; ++d)
1344 * flux[0][d] = conserved_variables[1 + d];
1345 * for (unsigned int e = 0; e < dim; ++e)
1346 * flux[e + 1][d] = conserved_variables[e + 1] * velocity[d];
1347 * flux[d + 1][d] += pressure;
1348 * flux[dim + 1][d] =
1349 * velocity[d] * (conserved_variables[dim + 1] + pressure);
1357 * This next function is a helper to simplify the implementation of the
1358 * numerical flux, implementing the action of a tensor of tensors (with
1359 * non-standard outer dimension of size `dim + 2`, so the standard overloads
1371 *
for (
unsigned int d = 0;
d < n_components; ++
d)
1372 *
result[d] = matrix[d] * vector;
1432 *
DG schemes in
that case.
1456 * at a neighbor cell.
1459 *
template <
int dim,
typename Number>
1486 *
0.5 * lambda * (
u_m -
u_p);
1496 *
const Number
s_pos =
1498 *
const Number
s_neg =
1528 *
template <
int dim,
typename Number>
1532 *
const unsigned int component)
1538 *
for (
unsigned int d = 0;
d < dim; ++
d)
1540 *
result[v] = function.value(p, component);
1546 *
template <
int dim,
typename Number,
int n_components = dim + 2>
1556 *
for (
unsigned int d = 0;
d < dim; ++
d)
1558 *
for (
unsigned int d = 0;
d < n_components; ++
d)
1559 *
result[d][v] = function.value(p, d);
1569 * <a name=
"step_67-TheEulerOperationclass"></a>
1599 *
template <
int dim,
int degree,
int n_po
ints_1d>
1621 *
void apply(
const double current_time,
1652 *
std::map<types::boundary_id, std::unique_ptr<Function<dim>>>
1654 *
std::map<types::boundary_id, std::unique_ptr<Function<dim>>>
1663 *
const std::pair<unsigned int, unsigned int> &
cell_range)
const;
1665 *
void local_apply_cell(
1669 *
const std::pair<unsigned int, unsigned int> &
cell_range)
const;
1675 *
const std::pair<unsigned int, unsigned int> &
face_range)
const;
1681 *
const std::pair<unsigned int, unsigned int> &
face_range)
const;
1686 *
template <
int dim,
int degree,
int n_po
ints_1d>
1720 *
const std::vector<const DoFHandler<dim> *> dof_handlers = {&dof_handler};
1722 *
const std::vector<const AffineConstraints<double> *> constraints = {&dummy};
1727 *
additional_data.mapping_update_flags =
1730 *
additional_data.mapping_update_flags_inner_faces =
1733 *
additional_data.mapping_update_flags_boundary_faces =
1736 *
additional_data.tasks_parallel_scheme =
1740 *
mapping, dof_handlers, constraints,
quadratures, additional_data);
1745 *
template <
int dim,
int degree,
int n_po
ints_1d>
1749 *
data.initialize_dof_vector(vector);
1781 *
template <
int dim,
int degree,
int n_po
ints_1d>
1789 *
ExcMessage(
"You already set the boundary with id " +
1790 *
std::to_string(
static_cast<int>(boundary_id)) +
1791 *
" to another type of boundary before now setting " +
1794 *
ExcMessage(
"Expected function with dim+2 components"));
1800 *
template <
int dim,
int degree,
int n_po
ints_1d>
1808 *
ExcMessage(
"You already set the boundary with id " +
1809 *
std::to_string(
static_cast<int>(boundary_id)) +
1810 *
" to another type of boundary before now setting " +
1811 *
"it as subsonic outflow"));
1813 *
ExcMessage(
"Expected function with dim+2 components"));
1819 *
template <
int dim,
int degree,
int n_po
ints_1d>
1827 *
ExcMessage(
"You already set the boundary with id " +
1828 *
std::to_string(
static_cast<int>(boundary_id)) +
1829 *
" to another type of boundary before now setting " +
1830 *
"it as wall boundary"));
1836 *
template <
int dim,
int degree,
int n_po
ints_1d>
1850 * <a name=
"step_67-Localevaluators"></a>
1861 * affine element shapes),
we now set
the number quadrature points
as a
1909 * quadrature point
data.
1976 *
for (
unsigned int d = 0;
d < dim; ++
d)
1978 *
for (
unsigned int d = 0;
d < dim; ++
d)
2046 *
const std::pair<unsigned int, unsigned int> &
face_range)
const
2055 *
phi_p.reinit(face);
2058 *
phi_m.reinit(face);
2061 *
for (
const unsigned int q :
phi_m.quadrature_point_indices())
2066 *
phi_m.normal_vector(
q));
2099 * direction
of the normal vector.
2141 *
template <
int dim,
int degree,
int n_po
ints_1d>
2146 *
const std::pair<unsigned int, unsigned int> &
face_range)
const
2155 *
for (
const unsigned int q :
phi.quadrature_point_indices())
2158 *
const auto normal =
phi.normal_vector(
q);
2161 *
for (
unsigned int d = 1;
d < dim; ++
d)
2171 *
for (
unsigned int d = 0;
d < dim; ++
d)
2173 *
w_p[dim + 1] =
w_m[dim + 1];
2179 *
phi.quadrature_point(
q));
2186 *
phi.quadrature_point(
q),
2193 *
"you set a boundary condition for "
2194 *
"this part of the domain boundary?"));
2202 *
for (
unsigned int d = 0;
d < dim; ++
d)
2203 *
flux[d + 1][v] = 0.;
2249 *
template <
int dim,
int degree,
int n_po
ints_1d>
2254 *
const std::pair<unsigned int, unsigned int> &
cell_range)
const
2263 *
phi.read_dof_values(src);
2265 *
inverse.apply(
phi.begin_dof_values(),
phi.begin_dof_values());
2267 *
phi.set_dof_values(dst);
2276 * <a name=
"step_67-Theapplyandrelatedfunctions"></a>
2321 *
const double current_time,
2329 *
i.
second->set_time(current_time);
2331 *
i.
second->set_time(current_time);
2333 *
data.loop(&EulerOperator::local_apply_cell,
2405 *
const Number current_time,
2417 *
i.
second->set_time(current_time);
2419 *
i.
second->set_time(current_time);
2421 *
data.loop(&EulerOperator::local_apply_cell,
2440 *
std::function<
void(
const unsigned int,
const unsigned int)>(),
2444 *
if (
ai == Number())
2449 *
const Number
k_i =
next_ri.local_element(i);
2450 *
const Number
sol_i = solution.local_element(i);
2451 *
solution.local_element(i) =
sol_i +
bi *
k_i;
2459 *
const Number
k_i =
next_ri.local_element(i);
2460 *
const Number
sol_i = solution.local_element(i);
2461 *
solution.local_element(i) =
sol_i +
bi *
k_i;
2504 *
\left(J^K\right)^{-1} S^{-1} S J^
K
2529 *
solution.zero_out_ghost_values();
2537 *
inverse.transform_from_q_points_to_basis(dim + 2,
2538 *
phi.begin_dof_values(),
2539 *
phi.begin_dof_values());
2540 *
phi.set_dof_values(solution);
2582 *
for (
unsigned int cell = 0; cell <
data.n_cell_batches(); ++cell)
2587 *
for (
const unsigned int q :
phi.quadrature_point_indices())
2592 *
const auto JxW =
phi.JxW(
q);
2595 *
for (
unsigned int d = 0;
d < dim; ++
d)
2599 *
for (
unsigned int v = 0; v <
data.n_active_entries_per_cell_batch(cell);
2601 *
for (
unsigned int d = 0;
d < 3; ++
d)
2607 *
std::array<double, 3>
errors;
2608 *
for (
unsigned int d = 0;
d < 3; ++
d)
2664 *
template <
int dim,
int degree,
int n_po
ints_1d>
2672 *
for (
unsigned int cell = 0; cell <
data.n_cell_batches(); ++cell)
2677 *
for (
const unsigned int q :
phi.quadrature_point_indices())
2683 *
const auto inverse_jacobian =
phi.inverse_jacobian(
q);
2686 *
for (
unsigned int d = 0;
d < dim; ++
d)
2694 *
for (
unsigned int d = 0;
d < dim; ++
d)
2696 *
for (
unsigned int i = 0; i < 5; ++i)
2701 *
for (
unsigned int d = 0;
d < dim; ++
d)
2707 *
const auto max_eigenvalue =
std::sqrt(
2720 *
for (
unsigned int v = 0; v <
data.n_active_entries_per_cell_batch(cell);
2735 * <a name=
"step_67-TheEulerProblemclass"></a>
2805 *
virtual void evaluate_vector_field(
2809 *
virtual std::vector<std::string> get_names()
const override;
2811 *
virtual std::vector<
2813 *
get_data_component_interpretation()
const override;
2815 *
virtual UpdateFlags get_needed_update_flags()
const override;
2824 *
template <
int dim>
2845 *
template <
int dim>
2854 *
ExcInternalError());
2857 *
ExcInternalError());
2858 *
Assert(
inputs.solution_values[0].size() == dim + 2, ExcInternalError());
2861 *
ExcInternalError());
2866 *
for (
unsigned int d = 0;
d < dim + 2; ++
d)
2867 *
solution[d] =
inputs.solution_values[p](d);
2869 *
const double density = solution[0];
2873 *
for (
unsigned int d = 0;
d < dim; ++
d)
2880 *
inputs.solution_gradients[p][0] *
inputs.solution_gradients[p][0];
2886 *
template <
int dim>
2889 *
std::vector<std::string> names;
2890 *
for (
unsigned int d = 0;
d < dim; ++
d)
2891 *
names.emplace_back(
"velocity");
2892 *
names.emplace_back(
"pressure");
2893 *
names.emplace_back(
"speed_of_sound");
2896 *
names.emplace_back(
"schlieren_plot");
2910 *
template <
int dim>
2911 *
std::vector<DataComponentInterpretation::DataComponentInterpretation>
2914 *
std::vector<DataComponentInterpretation::DataComponentInterpretation>
2916 *
for (
unsigned int d = 0;
d < dim; ++
d)
2938 *
template <
int dim>
2960 *
:
pcout(std::cout, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
2964 *
, fe(FE_DGQ<dim>(fe_degree) ^ (dim + 2))
2965 *
, mapping(fe_degree)
2966 *
, dof_handler(triangulation)
2967 *
, timer(pcout, TimerOutput::never, TimerOutput::wall_times)
2997 *
template <
int dim>
3005 *
for (
unsigned int d = 1;
d < dim; ++
d)
3010 *
for (
unsigned int d = 1;
d < dim; ++
d)
3040 *
std::vector<double>({0., 0., -0.2})));
3051 *
dof_handler.distribute_dofs(fe);
3066 *
std::locale s =
pcout.get_stream().getloc();
3067 *
pcout.get_stream().imbue(std::locale(
""));
3068 *
pcout <<
"Number of degrees of freedom: " << dof_handler.n_dofs()
3069 *
<<
" ( = " << (dim + 2) <<
" [vars] x "
3073 *
pcout.get_stream().imbue(s);
3101 * tutorial programs. A particularly nice feature of the
3102 * `write_vtu_in_parallel()` function is the fact that it can combine output
3103 * from all MPI ranks into a single file, making it unnecessary to have a
3104 * central record of all such files (namely, the "pvtu" file).
3108 * For parallel programs, it is often instructive to look at the partitioning
3109 * of cells among processors. To this end, one can pass a vector of numbers
3110 * to DataOut::add_data_vector() that contains as many entries as the
3111 * current processor has active cells; these numbers should then be the
3112 * rank of the processor that owns each of these cells. Such a vector
3113 * could, for example, be obtained from
3114 * GridTools::get_subdomain_association(). On the other hand, on each MPI
3115 * process, DataOut will only read those entries that correspond to locally
3116 * owned cells, and these of course all have the same value: namely, the rank
3117 * of the current process. What is in the remaining entries of the vector
3136 *
template <
int dim>
3139 *
const std::array<double, 3>
errors =
3143 *
pcout <<
"Time:" << std::setw(8) << std::setprecision(3) << time
3144 *
<<
", dt: " << std::setw(8) << std::setprecision(2) <<
time_step
3145 *
<<
", " <<
quantity_name <<
" rho: " << std::setprecision(4)
3146 *
<< std::setw(10) <<
errors[0] <<
", rho * u: " << std::setprecision(4)
3147 *
<< std::setw(10) <<
errors[1] <<
", energy:" << std::setprecision(4)
3148 *
<< std::setw(10) <<
errors[2] << std::endl;
3157 *
flags.write_higher_order_cells =
true;
3158 *
data_out.set_flags(flags);
3160 *
data_out.attach_dof_handler(dof_handler);
3162 *
std::vector<std::string> names;
3163 *
names.emplace_back(
"density");
3164 *
for (
unsigned int d = 0;
d < dim; ++
d)
3165 *
names.emplace_back(
"momentum");
3166 *
names.emplace_back(
"energy");
3168 *
std::vector<DataComponentInterpretation::DataComponentInterpretation>
3172 *
for (
unsigned int d = 0;
d < dim; ++
d)
3178 *
data_out.add_data_vector(dof_handler, solution, names,
interpretation);
3180 *
data_out.add_data_vector(solution, postprocessor);
3185 *
reference.reinit(solution);
3187 *
reference.sadd(-1., 1, solution);
3188 *
std::vector<std::string> names;
3189 *
names.emplace_back(
"error_density");
3190 *
for (
unsigned int d = 0;
d < dim; ++
d)
3191 *
names.emplace_back(
"error_momentum");
3192 *
names.emplace_back(
"error_energy");
3194 *
std::vector<DataComponentInterpretation::DataComponentInterpretation>
3198 *
for (
unsigned int d = 0;
d < dim; ++
d)
3204 *
data_out.add_data_vector(dof_handler,
3212 *
data_out.add_data_vector(
mpi_owner,
"owner");
3214 *
data_out.build_patches(mapping,
3243 *
template <
int dim>
3250 *
pcout <<
"Running with "
3252 *
<<
" MPI processes" << std::endl;
3254 *
<< (std::is_same_v<Number, double> ?
"doubles" :
"floats") <<
" = "
3272 *
for (
const auto &cell :
triangulation.active_cell_iterators())
3273 *
if (cell->is_locally_owned())
3283 *
<<
", initial transport scaling: "
3311 *
while (time < final_time - 1e-12)
3334 *
time >= final_time - 1e-12)
3336 *
static_cast<unsigned int>(std::round(time /
output_tick)));
3340 *
pcout << std::endl;
3359 *
using namespace dealii;
3368 *
catch (std::exception &exc)
3370 *
std::cerr << std::endl
3372 *
<<
"----------------------------------------------------"
3374 *
std::cerr <<
"Exception on processing: " << std::endl
3375 *
<< exc.what() << std::endl
3376 *
<<
"Aborting!" << std::endl
3377 *
<<
"----------------------------------------------------"
3384 *
std::cerr << std::endl
3386 *
<<
"----------------------------------------------------"
3388 *
std::cerr <<
"Unknown exception!" << std::endl
3389 *
<<
"Aborting!" << std::endl
3390 *
<<
"----------------------------------------------------"
3401<a name=
"step_67-Programoutput"></a><
h3>
Program output</
h3>
3409Number
of degrees
of freedom: 147,456 ( = 4 [vars]
x 1,024 [cells]
x 36 [dofs/cell/
var] )
3424+-------------------------------------------+------------------+------------+------------------+
3428+-------------------------------------------+------------------+------------+------------------+
3429| compute
errors | 11 | 0.008066s 13 | 0.00952s | 0.01041s 20 |
3430| compute
transport speed | 258 | 0.01012s 13 | 0.05392s | 0.08574s 25 |
3431| output | 11 | 0.9597s 13 | 0.9613s | 0.9623s 6 |
3435+-------------------------------------------+------------------+------------+------------------+
3473+-------------------------------------------+------------------+------------+------------------+
3477+-------------------------------------------+------------------+------------+------------------+
3478| compute
errors | 11 | 0.4239s 12 | 0.4318s | 0.4408s 9 |
3479| compute
transport speed | 2053 | 3.962s 12 | 6.727s | 10.12s 7 |
3480| output | 11 | 30.35s 12 | 30.36s | 30.37s 9 |
3484+-------------------------------------------+------------------+------------+------------------+
3658<table
align=
"center" class=
"doxtable">
3806 constexpr unsigned int testcase = 1;
3884Number
of degrees
of freedom: 3,686,400 ( = 4 [vars]
x 25,600 [cells]
x 36 [dofs/cell/
var] )
3887Time: 0,
dt: 7.4e-05, norm
rho: 4.17e-16,
rho *
u: 1.629e-16, energy: 1.381e-15
3888Time: 0.05,
dt: 6.3e-05, norm
rho: 0.02075,
rho *
u: 0.03801, energy: 0.08772
3889Time: 0.1,
dt: 5.9e-05, norm
rho: 0.02211,
rho *
u: 0.04515, energy: 0.08953
3890Time: 0.15,
dt: 5.7e-05, norm
rho: 0.02261,
rho *
u: 0.04592, energy: 0.08967
3891Time: 0.2,
dt: 5.8e-05, norm
rho: 0.02058,
rho *
u: 0.04361, energy: 0.08222
3892Time: 0.25,
dt: 5.9e-05, norm
rho: 0.01695,
rho *
u: 0.04203, energy: 0.06873
3893Time: 0.3,
dt: 5.9e-05, norm
rho: 0.01653,
rho *
u: 0.0401, energy: 0.06604
3894Time: 0.35,
dt: 5.7e-05, norm
rho: 0.01774,
rho *
u: 0.04264, energy: 0.0706
3898Time: 1.95,
dt: 5.8e-05, norm
rho: 0.01488,
rho *
u: 0.03923, energy: 0.05185
3899Time: 2,
dt: 5.7e-05, norm
rho: 0.01432,
rho *
u: 0.03969, energy: 0.04889
3901+-------------------------------------------+------------------+------------+------------------+
3904| Section |
no.
calls | min time rank | avg time | max time rank |
3905+-------------------------------------------+------------------+------------+------------------+
3906| compute
errors | 41 | 0.01112s 35 | 0.0672s | 0.1337s 0 |
3907| compute
transport speed | 6914 | 5.422s 35 | 15.96s | 29.99s 1 |
3908| output | 41 | 37.24s 35 | 37.3s | 37.37s 0 |
3912+-------------------------------------------+------------------+------------+------------------+
3928Number
of degrees
of freedom: 14,745,600 ( = 4 [vars]
x 102,400 [cells]
x 36 [dofs/cell/
var] )
3932+-------------------------------------------+------------------+------------+------------------+
3935| Section |
no.
calls | min time rank | avg time | max time rank |
3936+-------------------------------------------+------------------+------------+------------------+
3937| compute
errors | 41 | 0.04537s 32 | 0.173s | 0.3489s 0 |
3938| compute
transport speed | 13858 | 40.75s 32 | 85.99s | 149.8s 0 |
3939| output | 41 | 153.8s 32 | 153.9s | 154.1s 0 |
3943+-------------------------------------------+------------------+------------+------------------+
3980Number
of degrees
of freedom: 58,982,400 ( = 4 [vars]
x 409,600 [cells]
x 36 [dofs/cell/
var] )
3984Time: 1.95,
dt: 1.4e-05, norm
rho: 0.01488,
rho *
u: 0.03923, energy: 0.05183
3985Time: 2,
dt: 1.4e-05, norm
rho: 0.01431,
rho *
u: 0.03969, energy: 0.04887
3987+-------------------------------------------+------------------+------------+------------------+
3990| Section |
no.
calls | min time rank | avg time | max time rank |
3991+-------------------------------------------+------------------+------------+------------------+
3992| compute
errors | 41 | 0.1758s 30 | 0.672s | 1.376s 1 |
3993| compute
transport speed | 27748 | 321.3s 34 | 678.8s | 1202s 1 |
3994| output | 41 | 616.3s 32 | 616.4s | 616.4s 34 |
3995|
rk time
stepping total | 138733 | 1.983e+04s 1 | 2.036e+04s | 2.072e+04s 34 |
3998+-------------------------------------------+------------------+------------+------------------+
4015Number
of degrees
of freedom: 221,184,000 ( = 5 [vars]
x 204,800 [cells]
x 216 [dofs/cell/
var] )
4019Time: 1.95,
dt: 0.00011, norm
rho: 0.01131,
rho *
u: 0.03056, energy: 0.04091
4020Time: 2,
dt: 0.00011, norm
rho: 0.0119,
rho *
u: 0.03142, energy: 0.04425
4022+-------------------------------------------+------------------+------------+------------------+
4025| Section |
no.
calls | min time rank | avg time | max time rank |
4026+-------------------------------------------+------------------+------------+------------------+
4027| compute
errors | 41 | 0.6551s 34 | 3.216s | 7.281s 0 |
4029| output | 41 | 1350s 34 | 1353s | 1357s 0 |
4030|
rk time
stepping total | 17723 | 1.519e+04s 0 | 1.558e+04s | 1.582e+04s 34 |
4033+-------------------------------------------+------------------+------------+------------------+
4062Number
of degrees
of freedom: 1,769,472,000 ( = 5 [vars]
x 1,638,400 [cells]
x 216 [dofs/cell/
var] )
4069+-------------------------------------------+------------------+------------+------------------+
4073+-------------------------------------------+------------------+------------+------------------+
4074| compute
errors | 41 | 2.632s 178 | 7.221s | 16.56s 0 |
4076| output | 41 | 8065s 176 | 8070s | 8079s 0 |
4077|
rk time
stepping total | 35350 | 4.25e+04s 0 | 4.43e+04s | 4.515e+04s 193 |
4080+-------------------------------------------+------------------+------------+------------------+
4087<a name=
"step_67-Possibilitiesforextensions"></a><
h3>Possibilities
for extensions</
h3>
4213href=
"https://en.wikipedia.org/wiki/Perfectly_matched_layer">
perfectly
4246<a name=
"step_67-PlainProg"></a>
void write_vtu_in_parallel(const std::string &filename, const MPI_Comm comm) const
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
void submit_dof_value(const value_type val_in, const unsigned int dof)
void set_dof_values(VectorType &dst, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip()) const
value_type get_value(const unsigned int q_point) const
const ShapeInfoType * data
Tensor< 2, dim, VectorizedArrayType > inverse_jacobian(const unsigned int q_point) const
std_cxx20::ranges::iota_view< unsigned int, unsigned int > quadrature_point_indices() const
Point< dim, VectorizedArrayType > quadrature_point(const unsigned int q) const
const unsigned int n_q_points
void gather_evaluate(const VectorType &input_vector, const EvaluationFlags::EvaluationFlags evaluation_flag)
void evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
Abstract base class for mapping classes.
void transform_from_q_points_to_basis(const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &inner_face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_face_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
Tensor< rank, dim, Number > sum(const Tensor< rank, dim, Number > &local, const MPI_Comm mpi_communicator)
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
void print_wall_time_statistics(const MPI_Comm mpi_comm, const double print_quantile=0.) const
static constexpr std::size_t size()
#define DEAL_II_ALWAYS_INLINE
#define DEAL_II_OPENMP_SIMD_PRAGMA
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
#define DEAL_II_NOT_IMPLEMENTED()
std::vector< index_type > data
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
DataComponentInterpretation
@ component_is_part_of_vector
void approximate(const SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
The namespace for the EvaluationFlags enum.
void hyper_rectangle(Triangulation< dim, spacedim > &tria, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void cylinder(Triangulation< dim > &tria, const double radius=1., const double half_length=1.)
void channel_with_cylinder(Triangulation< dim > &tria, const double shell_region_width=0.03, const unsigned int n_shells=2, const double skewness=2.0, const bool colorize=false)
@ valid
Iterator points to a valid object.
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
@ diagonal
Matrix is diagonal.
@ general
No special properties.
constexpr types::blas_int zero
constexpr types::blas_int one
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
void apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_data, const ViewTypeIn in, ViewTypeOut out)
@ LOW_STORAGE_RK_STAGE9_ORDER5
@ LOW_STORAGE_RK_STAGE3_ORDER3
@ LOW_STORAGE_RK_STAGE7_ORDER4
@ LOW_STORAGE_RK_STAGE5_ORDER4
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
T max(const T &t, const MPI_Comm mpi_communicator)
T min(const T &t, const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
std::string get_current_vectorization_level()
Number truncate_to_n_digits(const Number number, const unsigned int n_digits)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
constexpr T pow(const T base, const int iexp)
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)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
long double gamma(const unsigned int n)
DEAL_II_HOST constexpr TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
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 > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, 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
****code * * MPI_Finalize()
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)