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) {
303 const auto i = std::get<0>(entry);
307 const auto normal = std::get<1>(entry);
308 const auto id = std::get<4>(entry);
309 const auto position = std::get<5>(entry);
313 Tensor<1, dim, Number> V_i;
314 Tensor<1, dim, Number> RHS_i;
315 for (
unsigned int d = 0; d < dim; ++d) {
316 V_i[d] = velocity_.block(d).local_element(i);
317 RHS_i[d] = velocity_rhs_.block(d).local_element(i);
319 V_i -= 1. * (V_i * normal) * normal;
320 RHS_i -= 1. * (RHS_i * normal) * normal;
321 for (
unsigned int d = 0; d < dim; ++d) {
322 velocity_.block(d).local_element(i) = V_i[d];
323 velocity_rhs_.block(d).local_element(i) = RHS_i[d];
329 for (
unsigned int d = 0; d < dim; ++d) {
330 velocity_.block(d).local_element(i) = Number(0.);
331 velocity_rhs_.block(d).local_element(i) = Number(0.);
337 const auto U_i = initial_values_->initial_state(position, t + tau);
338 const auto view = hyperbolic_system_->template view<dim, Number>();
339 const auto rho_i = view.density(U_i);
340 const auto V_i = view.momentum(U_i) / rho_i;
341 const auto e_i = view.internal_energy(U_i) / rho_i;
343 for (
unsigned int d = 0; d < dim; ++d) {
344 velocity_.block(d).local_element(i) = V_i[d];
345 velocity_rhs_.block(d).local_element(i) = V_i[d];
348 internal_energy_.local_element(i) = e_i;
359 affine_constraints.set_zero(density_);
360 affine_constraints.set_zero(internal_energy_);
361 for (
unsigned int d = 0; d < dim; ++d) {
362 affine_constraints.set_zero(velocity_.block(d));
363 affine_constraints.set_zero(velocity_rhs_.block(d));
369 lumped_mass_matrix, density_, affine_constraints);
376 if (use_gmg_velocity_ && reinitialize_gmg) {
377 MGLevelObject<
typename PreconditionChebyshev<
379 LinearAlgebra::distributed::BlockVector<float>,
381 smoother_data(level_matrix_free_.min_level(),
382 level_matrix_free_.max_level());
384 level_velocity_matrices_.resize(level_matrix_free_.min_level(),
385 level_matrix_free_.max_level());
386 mg_transfer_velocity_.interpolate_to_mg(
387 offline_data_->dof_handler(), level_density_, density_);
389 for (
unsigned int level = level_matrix_free_.min_level();
390 level <= level_matrix_free_.max_level();
392 level_velocity_matrices_[level].initialize(
395 level_matrix_free_[level],
396 level_density_[level],
399 level_velocity_matrices_[level].compute_diagonal(
400 smoother_data[level].preconditioner);
401 if (level == level_matrix_free_.min_level()) {
402 smoother_data[level].degree = numbers::invalid_unsigned_int;
403 smoother_data[level].eig_cg_n_iterations = 500;
404 smoother_data[level].smoothing_range = 1e-3;
406 smoother_data[level].degree = gmg_smoother_degree_;
407 smoother_data[level].eig_cg_n_iterations =
408 gmg_smoother_n_cg_iter_;
409 smoother_data[level].smoothing_range = gmg_smoother_range_vel_;
410 if (gmg_smoother_n_cg_iter_ == 0)
411 smoother_data[level].max_eigenvalue = gmg_smoother_max_eig_vel_;
414 mg_smoother_velocity_.initialize(level_velocity_matrices_,
424 Scope scope(computing_timer_,
425 "time step [P] _ - synchronization barriers");
431 *std::min_element(internal_energy_.begin(), internal_energy_.end());
432 e_min_old = Utilities::MPI::min(e_min_old, mpi_communicator_);
435 constexpr Number eps = std::numeric_limits<Number>::epsilon();
436 e_min_old *= (1. - 1000. * eps);
443 Scope scope(computing_timer_,
"time step [P] 1 - update velocities");
449 *parabolic_system_, *offline_data_, matrix_free_, density_, tau);
451 const auto tolerance_velocity =
452 (tolerance_linfty_norm_ ? velocity_rhs_.linfty_norm()
453 : velocity_rhs_.l2_norm()) *
462 if (!use_gmg_velocity_)
463 throw SolverControl::NoConvergence(0, 0.);
465 using bvt_float = LinearAlgebra::distributed::BlockVector<float>;
467 MGCoarseGridApplySmoother<bvt_float> mg_coarse;
468 mg_coarse.initialize(mg_smoother_velocity_);
470 mg::Matrix<bvt_float> mg_matrix(level_velocity_matrices_);
472 Multigrid<bvt_float> mg(mg_matrix,
474 mg_transfer_velocity_,
475 mg_smoother_velocity_,
476 mg_smoother_velocity_,
477 level_velocity_matrices_.min_level(),
478 level_velocity_matrices_.max_level());
480 const auto &dof_handler = offline_data_->dof_handler();
481 PreconditionMG<dim, bvt_float, MGTransferVelocity<dim, float>>
482 preconditioner(dof_handler, mg, mg_transfer_velocity_);
484 SolverControl solver_control(gmg_max_iter_vel_, tolerance_velocity);
485 SolverCG<BlockVector> solver(solver_control);
487 velocity_operator, velocity_, velocity_rhs_, preconditioner);
490 n_iterations_velocity_ =
491 0.9 * n_iterations_velocity_ + 0.1 * solver_control.last_step();
493 }
catch (SolverControl::NoConvergence &) {
495 SolverControl solver_control(1000, tolerance_velocity);
496 SolverCG<BlockVector> solver(solver_control);
498 velocity_operator, velocity_, velocity_rhs_, diagonal_matrix);
501 n_iterations_velocity_ *= 0.9;
502 n_iterations_velocity_ +=
503 0.1 * (use_gmg_velocity_ ? gmg_max_iter_vel_ : 0) +
504 0.1 * solver_control.last_step();
514 Scope scope(computing_timer_,
515 "time step [P] 2 - update internal energy");
520 matrix_free_.template cell_loop<ScalarVector, BlockVector>(
521 [
this](
const auto &data,
524 const auto cell_range) {
525 FEEvaluation<dim, order_fe, order_quad, dim, Number> velocity(
527 FEEvaluation<dim, order_fe, order_quad, 1, Number> energy(data);
529 const auto mu = parabolic_system_->mu();
530 const auto lambda = parabolic_system_->lambda();
532 for (
unsigned int cell = cell_range.first;
533 cell < cell_range.second;
535 velocity.reinit(cell);
537 velocity.gather_evaluate(src, EvaluationFlags::gradients);
539 for (
unsigned int q = 0; q < velocity.n_q_points; ++q) {
540 if constexpr (dim == 1) {
542 const auto gradient = velocity.get_gradient(q);
543 auto S = (4. / 3. * mu + lambda) * gradient;
544 energy.submit_value(gradient * S, q);
548 const auto symmetric_gradient =
549 velocity.get_symmetric_gradient(q);
550 const auto divergence = trace(symmetric_gradient);
551 auto S = 2. * mu * symmetric_gradient;
552 for (
unsigned int d = 0; d < dim; ++d)
553 S[d][d] += (lambda - 2. / 3. * mu) * divergence;
554 energy.submit_value(symmetric_gradient * S, q);
557 energy.integrate_scatter(EvaluationFlags::values, dst);
560 internal_energy_rhs_,
564 const auto &lumped_mass_matrix = offline_data_->lumped_mass_matrix();
568 auto loop = [&](
auto sentinel,
unsigned int left,
unsigned int right) {
569 using T =
decltype(sentinel);
570 unsigned int stride_size = get_stride_size<T>;
572 const auto view = hyperbolic_system_->template view<dim, T>();
575 for (
unsigned int i = left; i < right; i += stride_size) {
576 const auto rhs_i = get_entry<T>(internal_energy_rhs_, i);
577 const auto m_i = get_entry<T>(lumped_mass_matrix, i);
578 const auto rho_i = get_entry<T>(density_, i);
579 const auto e_i = get_entry<T>(internal_energy_, i);
581 const auto U_i = old_U.template get_tensor<T>(i);
582 const auto V_i = view.momentum(U_i) / rho_i;
584 dealii::Tensor<1, dim, T> V_i_new;
585 for (
unsigned int d = 0; d < dim; ++d) {
586 V_i_new[d] = get_entry<T>(velocity_.block(d), i);
593 const auto correction = Number(0.5) * (V_i - V_i_new).norm_square();
596 const auto result = m_i * rho_i * (e_i + correction) + tau * rhs_i;
597 write_entry<T>(internal_energy_rhs_, result, i);
602 loop(Number(), n_regular, n_owned);
604 loop(VA(), 0, n_regular);
615 const auto &boundary_map = offline_data_->boundary_map();
617 for (
auto entry : boundary_map) {
619 const auto i = std::get<0>(entry);
623 const auto id = std::get<4>(entry);
624 const auto position = std::get<5>(entry);
628 const auto U_i = initial_values_->initial_state(position, t + tau);
629 const auto view = hyperbolic_system_->template view<dim, Number>();
630 const auto rho_i = view.density(U_i);
631 const auto e_i = view.internal_energy(U_i) / rho_i;
632 internal_energy_rhs_.local_element(i) = e_i;
642 affine_constraints.set_zero(internal_energy_rhs_);
649 if (use_gmg_internal_energy_ && reinitialize_gmg) {
650 MGLevelObject<
typename PreconditionChebyshev<
652 LinearAlgebra::distributed::Vector<float>>::AdditionalData>
653 smoother_data(level_matrix_free_.min_level(),
654 level_matrix_free_.max_level());
656 level_energy_matrices_.resize(level_matrix_free_.min_level(),
657 level_matrix_free_.max_level());
659 for (
unsigned int level = level_matrix_free_.min_level();
660 level <= level_matrix_free_.max_level();
662 level_energy_matrices_[level].initialize(
664 level_matrix_free_[level],
665 level_density_[level],
666 tau * parabolic_system_->cv_inverse_kappa(),
668 level_energy_matrices_[level].compute_diagonal(
669 smoother_data[level].preconditioner);
670 if (level == level_matrix_free_.min_level()) {
671 smoother_data[level].degree = numbers::invalid_unsigned_int;
672 smoother_data[level].eig_cg_n_iterations = 500;
673 smoother_data[level].smoothing_range = 1e-3;
675 smoother_data[level].degree = gmg_smoother_degree_;
676 smoother_data[level].eig_cg_n_iterations =
677 gmg_smoother_n_cg_iter_;
678 smoother_data[level].smoothing_range = gmg_smoother_range_en_;
679 if (gmg_smoother_n_cg_iter_ == 0)
680 smoother_data[level].max_eigenvalue = gmg_smoother_max_eig_en_;
683 mg_smoother_energy_.initialize(level_energy_matrices_, smoother_data);
693 Scope scope(computing_timer_,
694 "time step [P] 2 - update internal energy");
699 const auto &kappa = parabolic_system_->cv_inverse_kappa();
701 *offline_data_, matrix_free_, density_, tau * kappa);
703 const auto tolerance_internal_energy =
704 (tolerance_linfty_norm_ ? internal_energy_rhs_.linfty_norm()
705 : internal_energy_rhs_.l2_norm()) *
709 if (!use_gmg_internal_energy_)
710 throw SolverControl::NoConvergence(0, 0.);
712 using vt_float = LinearAlgebra::distributed::Vector<float>;
713 MGCoarseGridApplySmoother<vt_float> mg_coarse;
714 mg_coarse.initialize(mg_smoother_energy_);
715 mg::Matrix<vt_float> mg_matrix(level_energy_matrices_);
717 Multigrid<vt_float> mg(mg_matrix,
722 level_energy_matrices_.min_level(),
723 level_energy_matrices_.max_level());
725 const auto &dof_handler = offline_data_->dof_handler();
726 PreconditionMG<dim, vt_float, MGTransferEnergy<dim, float>>
727 preconditioner(dof_handler, mg, mg_transfer_energy_);
729 SolverControl solver_control(gmg_max_iter_en_,
730 tolerance_internal_energy);
731 SolverCG<ScalarVector> solver(solver_control);
732 solver.solve(energy_operator,
734 internal_energy_rhs_,
738 n_iterations_internal_energy_ = 0.9 * n_iterations_internal_energy_ +
739 0.1 * solver_control.last_step();
741 }
catch (SolverControl::NoConvergence &) {
743 SolverControl solver_control(1000, tolerance_internal_energy);
744 SolverCG<ScalarVector> solver(solver_control);
745 solver.solve(energy_operator,
747 internal_energy_rhs_,
751 n_iterations_internal_energy_ *= 0.9;
752 n_iterations_internal_energy_ +=
753 0.1 * (use_gmg_internal_energy_ ? gmg_max_iter_en_ : 0) +
754 0.1 * solver_control.last_step();
761 Scope scope(computing_timer_,
762 "time step [P] _ - synchronization barriers");
765 auto e_min_new = *std::min_element(internal_energy_.begin(),
766 internal_energy_.end());
767 e_min_new = Utilities::MPI::min(e_min_new, mpi_communicator_);
769 if (e_min_new < e_min_old) {
771 std::cout << std::fixed << std::setprecision(16);
772 std::cout <<
"Bounds violation: internal energy (critical)!\n"
773 <<
"\t\te_min_old: " << e_min_old <<
"\n"
774 <<
"\t\te_min_old (delta): "
776 <<
"\t\te_min_new: " << e_min_new <<
"\n"
779 restart_needed =
true;
792 Scope scope(computing_timer_,
"time step [P] 3 - write back vectors");
797 auto loop = [&](
auto sentinel,
unsigned int left,
unsigned int right) {
798 using T =
decltype(sentinel);
799 unsigned int stride_size = get_stride_size<T>;
801 const auto view = hyperbolic_system_->template view<dim, T>();
804 for (
unsigned int i = left; i < right; i += stride_size) {
805 auto U_i = old_U.template get_tensor<T>(i);
806 const auto rho_i = view.density(U_i);
808 Tensor<1, dim, T> m_i_new;
809 for (
unsigned int d = 0; d < dim; ++d) {
810 m_i_new[d] = rho_i * get_entry<T>(velocity_.block(d), i);
813 const auto rho_e_i_new = rho_i * get_entry<T>(internal_energy_, i);
815 const auto E_i_new = rho_e_i_new + 0.5 * m_i_new * m_i_new / rho_i;
817 for (
unsigned int d = 0; d < dim; ++d)
818 U_i[1 + d] = m_i_new[d];
819 U_i[1 + dim] = E_i_new;
821 new_U.template write_tensor<T>(U_i, i);
826 loop(Number(), n_regular, n_owned);
828 loop(VA(), 0, n_regular);
833 new_U.update_ghost_values();
838 if (restart_needed) {
839 switch (id_violation_strategy) {
851 template <
typename Description,
int dim,
typename Number>
853 std::ostream &output)
const
855 output <<
" [ " << std::setprecision(2) << std::fixed
856 << n_iterations_velocity_
857 << (use_gmg_velocity_ ?
" GMG vel -- " :
" CG vel -- ")
858 << n_iterations_internal_energy_
859 << (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