15#include <deal.II/lac/linear_operator.h>
16#include <deal.II/lac/precondition.h>
17#include <deal.II/lac/solver_cg.h>
18#include <deal.II/matrix_free/fe_evaluation.h>
19#include <deal.II/multigrid/mg_coarse.h>
20#include <deal.II/multigrid/mg_matrix.h>
21#include <deal.II/multigrid/mg_transfer.templates.h>
22#include <deal.II/multigrid/mg_transfer_matrix_free.h>
23#include <deal.II/multigrid/multigrid.h>
29 namespace NavierStokes
31 using namespace dealii;
33 template <
typename Description,
int dim,
typename Number>
35 const MPI_Comm &mpi_communicator,
36 std::map<std::string, dealii::Timer> &computing_timer,
41 const std::string &subsection )
42 : ParameterAcceptor(subsection)
43 , mpi_communicator_(mpi_communicator)
44 , computing_timer_(computing_timer)
45 , hyperbolic_system_(&hyperbolic_system)
46 , parabolic_system_(¶bolic_system)
47 , offline_data_(&offline_data)
48 , initial_values_(&initial_values)
51 , n_iterations_velocity_(0.)
52 , n_iterations_internal_energy_(0.)
54 use_gmg_velocity_ =
false;
55 add_parameter(
"multigrid velocity",
57 "Use geometric multigrid for velocity component");
59 gmg_max_iter_vel_ = 12;
60 add_parameter(
"multigrid velocity - max iter",
62 "Maximal number of CG iterations with GMG smoother");
64 gmg_smoother_range_vel_ = 8.;
65 add_parameter(
"multigrid velocity - chebyshev range",
66 gmg_smoother_range_vel_,
67 "Chebyshev smoother: eigenvalue range parameter");
69 gmg_smoother_max_eig_vel_ = 2.0;
70 add_parameter(
"multigrid velocity - chebyshev max eig",
71 gmg_smoother_max_eig_vel_,
72 "Chebyshev smoother: maximal eigenvalue");
74 use_gmg_internal_energy_ =
false;
75 add_parameter(
"multigrid energy",
76 use_gmg_internal_energy_,
77 "Use geometric multigrid for internal energy component");
79 gmg_max_iter_en_ = 15;
80 add_parameter(
"multigrid energy - max iter",
82 "Maximal number of CG iterations with GMG smoother");
84 gmg_smoother_range_en_ = 15.;
85 add_parameter(
"multigrid energy - chebyshev range",
86 gmg_smoother_range_en_,
87 "Chebyshev smoother: eigenvalue range parameter");
89 gmg_smoother_max_eig_en_ = 2.0;
90 add_parameter(
"multigrid energy - chebyshev max eig",
91 gmg_smoother_max_eig_en_,
92 "Chebyshev smoother: maximal eigenvalue");
94 gmg_smoother_degree_ = 3;
95 add_parameter(
"multigrid - chebyshev degree",
97 "Chebyshev smoother: degree");
99 gmg_smoother_n_cg_iter_ = 10;
101 "multigrid - chebyshev cg iter",
102 gmg_smoother_n_cg_iter_,
103 "Chebyshev smoother: number of CG iterations to approximate "
108 "multigrid - min level",
110 "Minimal mesh level to be visited in the geometric multigrid "
111 "cycle where the coarse grid solver (Chebyshev) is called");
113 tolerance_ = Number(1.0e-12);
114 add_parameter(
"tolerance", tolerance_,
"Tolerance for linear solvers");
116 tolerance_linfty_norm_ =
false;
117 add_parameter(
"tolerance linfty norm",
118 tolerance_linfty_norm_,
119 "Use the l_infty norm instead of the l_2 norm for the "
120 "stopping criterion");
124 template <
typename Description,
int dim,
typename Number>
128 std::cout <<
"ParabolicSolver<dim, Number>::prepare()" << std::endl;
131 const auto &discretization = offline_data_->discretization();
133 dealii::ExcMessage(
"The NavierStokes module currently only "
134 "supports cG Q1 finite elements."));
138 typename MatrixFree<dim, Number>::AdditionalData additional_data;
139 additional_data.tasks_parallel_scheme =
140 MatrixFree<dim, Number>::AdditionalData::none;
142 matrix_free_.reinit(discretization.mapping(),
143 offline_data_->dof_handler(),
144 offline_data_->affine_constraints(),
145 discretization.quadrature_1d(),
148 const auto &scalar_partitioner =
149 matrix_free_.get_dof_info(0).vector_partitioner;
151 velocity_.reinit(dim);
152 velocity_rhs_.reinit(dim);
153 for (
unsigned int i = 0; i < dim; ++i) {
154 velocity_.block(i).reinit(scalar_partitioner);
155 velocity_rhs_.block(i).reinit(scalar_partitioner);
158 internal_energy_.reinit(scalar_partitioner);
159 internal_energy_rhs_.reinit(scalar_partitioner);
161 density_.reinit(scalar_partitioner);
165 if (!use_gmg_velocity_ && !use_gmg_internal_energy_)
168 const unsigned int n_levels =
169 offline_data_->dof_handler().get_triangulation().n_global_levels();
170 const unsigned int min_level = std::min(gmg_min_level_, n_levels - 1);
171 MGLevelObject<IndexSet> relevant_sets(0, n_levels - 1);
172 for (
unsigned int level = 0; level < n_levels; ++level)
173 dealii::DoFTools::extract_locally_relevant_level_dofs(
174 offline_data_->dof_handler(), level, relevant_sets[level]);
175 mg_constrained_dofs_.initialize(offline_data_->dof_handler(),
177 std::set<types::boundary_id> boundary_ids;
180 mg_constrained_dofs_.make_zero_boundary_constraints(
181 offline_data_->dof_handler(), boundary_ids);
183 typename MatrixFree<dim, float>::AdditionalData additional_data_level;
184 additional_data_level.tasks_parallel_scheme =
185 MatrixFree<dim, float>::AdditionalData::none;
187 level_matrix_free_.resize(min_level, n_levels - 1);
188 level_density_.resize(min_level, n_levels - 1);
189 for (
unsigned int level = min_level; level < n_levels; ++level) {
190 additional_data_level.mg_level = level;
191 AffineConstraints<double> constraints(relevant_sets[level]);
195 level_matrix_free_[level].reinit(discretization.mapping(),
196 offline_data_->dof_handler(),
198 discretization.quadrature_1d(),
199 additional_data_level);
200 level_matrix_free_[level].initialize_dof_vector(level_density_[level]);
203 mg_transfer_velocity_.build(offline_data_->dof_handler(),
204 mg_constrained_dofs_,
206 mg_transfer_energy_.build(offline_data_->dof_handler(),
211 template <
typename Description,
int dim,
typename Number>
218 const bool reinitialize_gmg)
const
221 std::cout <<
"ParabolicSolver<dim, Number>::step()" << std::endl;
224 const auto &old_U = std::get<0>(old_state_vector);
225 auto &new_U = std::get<0>(new_state_vector);
229 using VA = VectorizedArray<Number>;
231 const auto &lumped_mass_matrix = offline_data_->lumped_mass_matrix();
232 const auto &affine_constraints = offline_data_->affine_constraints();
236 constexpr auto simd_length = VA::size();
237 const unsigned int n_owned = offline_data_->n_locally_owned();
238 const unsigned int n_regular = n_owned / simd_length * simd_length;
243 std::cout <<
" perform time-step with tau = " << tau << std::endl;
247 std::atomic<bool> restart_needed =
false;
257 Scope scope(computing_timer_,
"time step [P] 1 - update velocities");
261 auto loop = [&](
auto sentinel,
unsigned int left,
unsigned int right) {
262 using T =
decltype(sentinel);
263 unsigned int stride_size = get_stride_size<T>;
265 const auto view = hyperbolic_system_->template view<dim, T>();
268 for (
unsigned int i = left; i < right; i += stride_size) {
269 const auto U_i = old_U.template get_tensor<T>(i);
270 const auto rho_i = view.density(U_i);
271 const auto M_i = view.momentum(U_i);
272 const auto rho_e_i = view.internal_energy(U_i);
273 const auto m_i = get_entry<T>(lumped_mass_matrix, i);
275 write_entry<T>(density_, rho_i, i);
277 for (
unsigned int d = 0; d < dim; ++d) {
278 write_entry<T>(velocity_.block(d), M_i[d] / rho_i, i);
279 write_entry<T>(velocity_rhs_.block(d), m_i * (M_i[d]), i);
281 write_entry<T>(internal_energy_, rho_e_i / rho_i, i);
286 loop(Number(), n_regular, n_owned);
288 loop(VA(), 0, n_regular);
299 const auto &boundary_map = offline_data_->boundary_map();
301 for (
auto entry : boundary_map) {
302 const auto i = entry.first;
306 const auto normal = std::get<0>(entry.second);
307 const auto id = std::get<3>(entry.second);
308 const auto position = std::get<4>(entry.second);
312 Tensor<1, dim, Number> V_i;
313 Tensor<1, dim, Number> RHS_i;
314 for (
unsigned int d = 0; d < dim; ++d) {
315 V_i[d] = velocity_.block(d).local_element(i);
316 RHS_i[d] = velocity_rhs_.block(d).local_element(i);
318 V_i -= 1. * (V_i * normal) * normal;
319 RHS_i -= 1. * (RHS_i * normal) * normal;
320 for (
unsigned int d = 0; d < dim; ++d) {
321 velocity_.block(d).local_element(i) = V_i[d];
322 velocity_rhs_.block(d).local_element(i) = RHS_i[d];
328 for (
unsigned int d = 0; d < dim; ++d) {
329 velocity_.block(d).local_element(i) = Number(0.);
330 velocity_rhs_.block(d).local_element(i) = Number(0.);
336 const auto U_i = initial_values_->initial_state(position, t + tau);
337 const auto view = hyperbolic_system_->template view<dim, Number>();
338 const auto rho_i = view.density(U_i);
339 const auto V_i = view.momentum(U_i) / rho_i;
340 const auto e_i = view.internal_energy(U_i) / rho_i;
342 for (
unsigned int d = 0; d < dim; ++d) {
343 velocity_.block(d).local_element(i) = V_i[d];
344 velocity_rhs_.block(d).local_element(i) = V_i[d];
347 internal_energy_.local_element(i) = e_i;
358 affine_constraints.set_zero(density_);
359 affine_constraints.set_zero(internal_energy_);
360 for (
unsigned int d = 0; d < dim; ++d) {
361 affine_constraints.set_zero(velocity_.block(d));
362 affine_constraints.set_zero(velocity_rhs_.block(d));
368 lumped_mass_matrix, density_, affine_constraints);
375 if (use_gmg_velocity_ && reinitialize_gmg) {
376 MGLevelObject<
typename PreconditionChebyshev<
378 LinearAlgebra::distributed::BlockVector<float>,
380 smoother_data(level_matrix_free_.min_level(),
381 level_matrix_free_.max_level());
383 level_velocity_matrices_.resize(level_matrix_free_.min_level(),
384 level_matrix_free_.max_level());
385 mg_transfer_velocity_.interpolate_to_mg(
386 offline_data_->dof_handler(), level_density_, density_);
388 for (
unsigned int level = level_matrix_free_.min_level();
389 level <= level_matrix_free_.max_level();
391 level_velocity_matrices_[level].initialize(
394 level_matrix_free_[level],
395 level_density_[level],
398 level_velocity_matrices_[level].compute_diagonal(
399 smoother_data[level].preconditioner);
400 if (level == level_matrix_free_.min_level()) {
401 smoother_data[level].degree = numbers::invalid_unsigned_int;
402 smoother_data[level].eig_cg_n_iterations = 500;
403 smoother_data[level].smoothing_range = 1e-3;
405 smoother_data[level].degree = gmg_smoother_degree_;
406 smoother_data[level].eig_cg_n_iterations =
407 gmg_smoother_n_cg_iter_;
408 smoother_data[level].smoothing_range = gmg_smoother_range_vel_;
409 if (gmg_smoother_n_cg_iter_ == 0)
410 smoother_data[level].max_eigenvalue = gmg_smoother_max_eig_vel_;
413 mg_smoother_velocity_.initialize(level_velocity_matrices_,
424 *std::min_element(internal_energy_.begin(), internal_energy_.end());
425 e_min_old = Utilities::MPI::min(e_min_old, mpi_communicator_);
428 constexpr Number eps = std::numeric_limits<Number>::epsilon();
429 e_min_old *= (1. - 1000. * eps);
435 Scope scope(computing_timer_,
"time step [P] 1 - update velocities");
441 *parabolic_system_, *offline_data_, matrix_free_, density_, tau);
443 const auto tolerance_velocity =
444 (tolerance_linfty_norm_ ? velocity_rhs_.linfty_norm()
445 : velocity_rhs_.l2_norm()) *
454 if (!use_gmg_velocity_)
455 throw SolverControl::NoConvergence(0, 0.);
457 using bvt_float = LinearAlgebra::distributed::BlockVector<float>;
459 MGCoarseGridApplySmoother<bvt_float> mg_coarse;
460 mg_coarse.initialize(mg_smoother_velocity_);
462 mg::Matrix<bvt_float> mg_matrix(level_velocity_matrices_);
464 Multigrid<bvt_float> mg(mg_matrix,
466 mg_transfer_velocity_,
467 mg_smoother_velocity_,
468 mg_smoother_velocity_,
469 level_velocity_matrices_.min_level(),
470 level_velocity_matrices_.max_level());
472 const auto &dof_handler = offline_data_->dof_handler();
473 PreconditionMG<dim, bvt_float, MGTransferVelocity<dim, float>>
474 preconditioner(dof_handler, mg, mg_transfer_velocity_);
476 SolverControl solver_control(gmg_max_iter_vel_, tolerance_velocity);
477 SolverCG<BlockVector> solver(solver_control);
479 velocity_operator, velocity_, velocity_rhs_, preconditioner);
482 n_iterations_velocity_ =
483 0.9 * n_iterations_velocity_ + 0.1 * solver_control.last_step();
485 }
catch (SolverControl::NoConvergence &) {
487 SolverControl solver_control(1000, tolerance_velocity);
488 SolverCG<BlockVector> solver(solver_control);
490 velocity_operator, velocity_, velocity_rhs_, diagonal_matrix);
493 n_iterations_velocity_ *= 0.9;
494 n_iterations_velocity_ +=
495 0.1 * (use_gmg_velocity_ ? gmg_max_iter_vel_ : 0) +
496 0.1 * solver_control.last_step();
506 Scope scope(computing_timer_,
507 "time step [P] 2 - update internal energy");
512 matrix_free_.template cell_loop<ScalarVector, BlockVector>(
513 [
this](
const auto &data,
516 const auto cell_range) {
517 FEEvaluation<dim, order_fe, order_quad, dim, Number> velocity(
519 FEEvaluation<dim, order_fe, order_quad, 1, Number> energy(data);
521 const auto mu = parabolic_system_->mu();
522 const auto lambda = parabolic_system_->lambda();
524 for (
unsigned int cell = cell_range.first;
525 cell < cell_range.second;
527 velocity.reinit(cell);
529 velocity.gather_evaluate(src, EvaluationFlags::gradients);
531 for (
unsigned int q = 0; q < velocity.n_q_points; ++q) {
532 if constexpr (dim == 1) {
534 const auto gradient = velocity.get_gradient(q);
535 auto S = (4. / 3. * mu + lambda) * gradient;
536 energy.submit_value(gradient * S, q);
540 const auto symmetric_gradient =
541 velocity.get_symmetric_gradient(q);
542 const auto divergence = trace(symmetric_gradient);
543 auto S = 2. * mu * symmetric_gradient;
544 for (
unsigned int d = 0; d < dim; ++d)
545 S[d][d] += (lambda - 2. / 3. * mu) * divergence;
546 energy.submit_value(symmetric_gradient * S, q);
549 energy.integrate_scatter(EvaluationFlags::values, dst);
552 internal_energy_rhs_,
556 const auto &lumped_mass_matrix = offline_data_->lumped_mass_matrix();
560 auto loop = [&](
auto sentinel,
unsigned int left,
unsigned int right) {
561 using T =
decltype(sentinel);
562 unsigned int stride_size = get_stride_size<T>;
564 const auto view = hyperbolic_system_->template view<dim, T>();
567 for (
unsigned int i = left; i < right; i += stride_size) {
568 const auto rhs_i = get_entry<T>(internal_energy_rhs_, i);
569 const auto m_i = get_entry<T>(lumped_mass_matrix, i);
570 const auto rho_i = get_entry<T>(density_, i);
571 const auto e_i = get_entry<T>(internal_energy_, i);
573 const auto U_i = old_U.template get_tensor<T>(i);
574 const auto V_i = view.momentum(U_i) / rho_i;
576 dealii::Tensor<1, dim, T> V_i_new;
577 for (
unsigned int d = 0; d < dim; ++d) {
578 V_i_new[d] = get_entry<T>(velocity_.block(d), i);
585 const auto correction = Number(0.5) * (V_i - V_i_new).norm_square();
588 const auto result = m_i * rho_i * (e_i + correction) + tau * rhs_i;
589 write_entry<T>(internal_energy_rhs_, result, i);
594 loop(Number(), n_regular, n_owned);
596 loop(VA(), 0, n_regular);
607 const auto &boundary_map = offline_data_->boundary_map();
609 for (
auto entry : boundary_map) {
610 const auto i = entry.first;
614 const auto id = std::get<3>(entry.second);
615 const auto position = std::get<4>(entry.second);
619 const auto U_i = initial_values_->initial_state(position, t + tau);
620 const auto view = hyperbolic_system_->template view<dim, Number>();
621 const auto rho_i = view.density(U_i);
622 const auto e_i = view.internal_energy(U_i) / rho_i;
623 internal_energy_rhs_.local_element(i) = e_i;
633 affine_constraints.set_zero(internal_energy_rhs_);
640 if (use_gmg_internal_energy_ && reinitialize_gmg) {
641 MGLevelObject<
typename PreconditionChebyshev<
643 LinearAlgebra::distributed::Vector<float>>::AdditionalData>
644 smoother_data(level_matrix_free_.min_level(),
645 level_matrix_free_.max_level());
647 level_energy_matrices_.resize(level_matrix_free_.min_level(),
648 level_matrix_free_.max_level());
650 for (
unsigned int level = level_matrix_free_.min_level();
651 level <= level_matrix_free_.max_level();
653 level_energy_matrices_[level].initialize(
655 level_matrix_free_[level],
656 level_density_[level],
657 tau * parabolic_system_->cv_inverse_kappa(),
659 level_energy_matrices_[level].compute_diagonal(
660 smoother_data[level].preconditioner);
661 if (level == level_matrix_free_.min_level()) {
662 smoother_data[level].degree = numbers::invalid_unsigned_int;
663 smoother_data[level].eig_cg_n_iterations = 500;
664 smoother_data[level].smoothing_range = 1e-3;
666 smoother_data[level].degree = gmg_smoother_degree_;
667 smoother_data[level].eig_cg_n_iterations =
668 gmg_smoother_n_cg_iter_;
669 smoother_data[level].smoothing_range = gmg_smoother_range_en_;
670 if (gmg_smoother_n_cg_iter_ == 0)
671 smoother_data[level].max_eigenvalue = gmg_smoother_max_eig_en_;
674 mg_smoother_energy_.initialize(level_energy_matrices_, smoother_data);
684 Scope scope(computing_timer_,
685 "time step [P] 2 - update internal energy");
690 const auto &kappa = parabolic_system_->cv_inverse_kappa();
692 *offline_data_, matrix_free_, density_, tau * kappa);
694 const auto tolerance_internal_energy =
695 (tolerance_linfty_norm_ ? internal_energy_rhs_.linfty_norm()
696 : internal_energy_rhs_.l2_norm()) *
700 if (!use_gmg_internal_energy_)
701 throw SolverControl::NoConvergence(0, 0.);
703 using vt_float = LinearAlgebra::distributed::Vector<float>;
704 MGCoarseGridApplySmoother<vt_float> mg_coarse;
705 mg_coarse.initialize(mg_smoother_energy_);
706 mg::Matrix<vt_float> mg_matrix(level_energy_matrices_);
708 Multigrid<vt_float> mg(mg_matrix,
713 level_energy_matrices_.min_level(),
714 level_energy_matrices_.max_level());
716 const auto &dof_handler = offline_data_->dof_handler();
717 PreconditionMG<dim, vt_float, MGTransferEnergy<dim, float>>
718 preconditioner(dof_handler, mg, mg_transfer_energy_);
720 SolverControl solver_control(gmg_max_iter_en_,
721 tolerance_internal_energy);
722 SolverCG<ScalarVector> solver(solver_control);
723 solver.solve(energy_operator,
725 internal_energy_rhs_,
729 n_iterations_internal_energy_ = 0.9 * n_iterations_internal_energy_ +
730 0.1 * solver_control.last_step();
732 }
catch (SolverControl::NoConvergence &) {
734 SolverControl solver_control(1000, tolerance_internal_energy);
735 SolverCG<ScalarVector> solver(solver_control);
736 solver.solve(energy_operator,
738 internal_energy_rhs_,
742 n_iterations_internal_energy_ *= 0.9;
743 n_iterations_internal_energy_ +=
744 0.1 * (use_gmg_internal_energy_ ? gmg_max_iter_en_ : 0) +
745 0.1 * solver_control.last_step();
753 auto e_min_new = *std::min_element(internal_energy_.begin(),
754 internal_energy_.end());
755 e_min_new = Utilities::MPI::min(e_min_new, mpi_communicator_);
757 if (e_min_new < e_min_old) {
759 std::cout << std::fixed << std::setprecision(16);
760 std::cout <<
"Bounds violation: internal energy (critical)!\n"
761 <<
"\t\te_min_old: " << e_min_old <<
"\n"
762 <<
"\t\te_min_old (delta): "
764 <<
"\t\te_min_new: " << e_min_new <<
"\n"
767 restart_needed =
true;
780 Scope scope(computing_timer_,
"time step [P] 3 - write back vectors");
785 auto loop = [&](
auto sentinel,
unsigned int left,
unsigned int right) {
786 using T =
decltype(sentinel);
787 unsigned int stride_size = get_stride_size<T>;
789 const auto view = hyperbolic_system_->template view<dim, T>();
792 for (
unsigned int i = left; i < right; i += stride_size) {
793 auto U_i = old_U.template get_tensor<T>(i);
794 const auto rho_i = view.density(U_i);
796 Tensor<1, dim, T> m_i_new;
797 for (
unsigned int d = 0; d < dim; ++d) {
798 m_i_new[d] = rho_i * get_entry<T>(velocity_.block(d), i);
801 const auto rho_e_i_new = rho_i * get_entry<T>(internal_energy_, i);
803 const auto E_i_new = rho_e_i_new + 0.5 * m_i_new * m_i_new / rho_i;
805 for (
unsigned int d = 0; d < dim; ++d)
806 U_i[1 + d] = m_i_new[d];
807 U_i[1 + dim] = E_i_new;
809 new_U.template write_tensor<T>(U_i, i);
814 loop(Number(), n_regular, n_owned);
816 loop(VA(), 0, n_regular);
821 new_U.update_ghost_values();
826 if (restart_needed) {
827 switch (id_violation_strategy) {
839 template <
typename Description,
int dim,
typename Number>
841 std::ostream &output)
const
843 output <<
" [ " << std::setprecision(2) << std::fixed
844 << n_iterations_velocity_
845 << (use_gmg_velocity_ ?
" GMG vel -- " :
" CG vel -- ")
846 << n_iterations_internal_energy_
847 << (use_gmg_internal_energy_ ?
" GMG int ]" :
" CG int ]")
void reinit(const vector_type &lumped_mass_matrix, const vector_type &density, const dealii::AffineConstraints< Number > &affine_constraints)
void initialize(const OfflineData< dim, Number2 > &offline_data, const dealii::MatrixFree< dim, Number > &matrix_free, const dealii::LinearAlgebra::distributed::Vector< Number > &density, const Number time_factor, const unsigned int level=dealii::numbers::invalid_unsigned_int)
void backward_euler_step(const StateVector &old_state_vector, const Number old_t, StateVector &new_state_vector, Number tau, const IDViolationStrategy id_violation_strategy, const bool reinitialize_gmg) const
typename Description::ParabolicSystem ParabolicSystem
typename View::StateVector StateVector
void print_solver_statistics(std::ostream &output) const
typename Description::HyperbolicSystem HyperbolicSystem
ParabolicSolver(const MPI_Comm &mpi_communicator, std::map< std::string, dealii::Timer > &computing_timer, const HyperbolicSystem &hyperbolic_system, const ParabolicSystem ¶bolic_system, const OfflineData< dim, Number > &offline_data, const InitialValues< Description, dim, Number > &initial_values, const std::string &subsection="ParabolicSolver")
void initialize(const ParabolicSystem ¶bolic_system, const OfflineData< dim, Number2 > &offline_data, const dealii::MatrixFree< dim, Number > &matrix_free, const dealii::LinearAlgebra::distributed::Vector< Number > &density, const Number theta_x_tau, const unsigned int level=dealii::numbers::invalid_unsigned_int)
#define RYUJIN_PARALLEL_REGION_BEGIN
#define RYUJIN_PARALLEL_REGION_END
DEAL_II_ALWAYS_INLINE Number negative_part(const Number number)
#define LIKWID_MARKER_START(opt)
#define CALLGRIND_START_INSTRUMENTATION
#define LIKWID_MARKER_STOP(opt)
#define CALLGRIND_STOP_INSTRUMENTATION