95class InverseMassOperator
98 typedef dealii::LinearAlgebra::distributed::Vector<Number> VectorType;
100 typedef InverseMassOperator<dim, n_components, Number> This;
102 typedef CellIntegrator<dim, n_components, Number> Integrator;
105 typedef dealii::MatrixFreeOperators::CellwiseInverseMassMatrix<dim, -1, n_components, Number>
106 InverseMassAsMatrixFreeOperator;
108 typedef std::pair<unsigned int, unsigned int> Range;
111 InverseMassOperator() : matrix_free(
nullptr), dof_index(0), quad_index(0)
116 initialize(dealii::MatrixFree<dim, Number>
const & matrix_free_in,
118 dealii::AffineConstraints<Number>
const * constraints =
nullptr)
120 this->matrix_free = &matrix_free_in;
121 dof_index = inverse_mass_operator_data.dof_index;
122 quad_index = inverse_mass_operator_data.quad_index;
124 data = inverse_mass_operator_data.parameters;
126 coefficient_is_variable = inverse_mass_operator_data.coefficient_is_variable;
127 consider_inverse_coefficient = inverse_mass_operator_data.consider_inverse_coefficient;
128 variable_coefficients = inverse_mass_operator_data.variable_coefficients;
131 AssertThrow(not coefficient_is_variable or variable_coefficients !=
nullptr,
132 dealii::ExcMessage(
"Pointer to variable coefficients not set properly."));
134 dealii::FiniteElement<dim>
const & fe = matrix_free->get_dof_handler(dof_index).get_fe();
139 if(data.implementation_type == InverseMassType::MatrixfreeOperator)
142 fe.base_element(0).dofs_per_cell == dealii::Utilities::pow(fe.degree + 1, dim),
144 "The matrix-free cell-wise inverse mass operator is currently only available for isotropic tensor-product elements."));
147 this->matrix_free->get_shape_info(0, quad_index).data[0].n_q_points_1d == fe.degree + 1,
149 "The matrix-free cell-wise inverse mass operator is currently only available if n_q_points_1d = n_nodes_1d."));
152 fe.conforms(dealii::FiniteElementData<dim>::L2),
154 "The matrix-free cell-wise inverse mass operator is only available for L2-conforming elements."));
161 mass_operator_data.dof_index = dof_index;
162 mass_operator_data.quad_index = quad_index;
163 mass_operator_data.coefficient_is_variable = coefficient_is_variable;
164 mass_operator_data.variable_coefficients = variable_coefficients;
165 mass_operator_data.consider_inverse_coefficient = consider_inverse_coefficient;
167 if(data.implementation_type == InverseMassType::ElementwiseKrylovSolver)
169 mass_operator_data.implement_block_diagonal_preconditioner_matrix_free =
true;
170 mass_operator_data.solver_block_diagonal = Elementwise::Solver::CG;
171 if(inverse_mass_operator_data.parameters.preconditioner == PreconditionerMass::None)
173 mass_operator_data.preconditioner_block_diagonal = Elementwise::Preconditioner::None;
175 else if(inverse_mass_operator_data.parameters.preconditioner ==
176 PreconditionerMass::PointJacobi)
178 mass_operator_data.preconditioner_block_diagonal =
179 Elementwise::Preconditioner::PointJacobi;
181 mass_operator_data.solver_data_block_diagonal =
182 inverse_mass_operator_data.parameters.solver_data;
186 if(constraints ==
nullptr)
188 dealii::AffineConstraints<Number> dummy_constraints;
189 dummy_constraints.close();
190 mass_operator.initialize(*matrix_free, dummy_constraints, mass_operator_data);
194 mass_operator.initialize(*matrix_free, *constraints, mass_operator_data);
197 if(data.implementation_type == InverseMassType::ElementwiseKrylovSolver or
198 data.implementation_type == InverseMassType::BlockMatrices)
203 fe.conforms(dealii::FiniteElementData<dim>::L2),
205 "The cell-wise inverse mass operator is only available for L2-conforming elements."));
207 block_jacobi_preconditioner =
208 std::make_shared<BlockJacobiPreconditioner<MassOperator<dim, n_components, Number>>>(
209 mass_operator,
true );
211 else if(data.implementation_type == InverseMassType::GlobalKrylovSolver)
214 solver_data.max_iter = inverse_mass_operator_data.parameters.solver_data.max_iter;
215 solver_data.solver_tolerance_abs =
216 inverse_mass_operator_data.parameters.solver_data.abs_tol;
217 solver_data.solver_tolerance_rel =
218 inverse_mass_operator_data.parameters.solver_data.rel_tol;
220 solver_data.use_preconditioner =
221 inverse_mass_operator_data.parameters.preconditioner != PreconditionerMass::None;
222 if(inverse_mass_operator_data.parameters.preconditioner == PreconditionerMass::None)
226 else if(inverse_mass_operator_data.parameters.preconditioner ==
227 PreconditionerMass::PointJacobi)
229 global_preconditioner =
230 std::make_shared<JacobiPreconditioner<MassOperator<dim, n_components, Number>>>(
231 mass_operator,
true );
233 else if(inverse_mass_operator_data.parameters.preconditioner ==
234 PreconditionerMass::BlockJacobi)
236 global_preconditioner =
237 std::make_shared<BlockJacobiPreconditioner<MassOperator<dim, n_components, Number>>>(
238 mass_operator,
true );
242 AssertThrow(
false, dealii::ExcMessage(
"This `PreconditionerMass` is not implemented."));
245 global_solver = std::make_shared<Krylov::SolverCG<MassOperator<dim, n_components, Number>,
247 VectorType>>(mass_operator,
248 *global_preconditioner,
254 dealii::ExcMessage(
"The specified `InverseMassType` is not implemented."));
266 if(data.implementation_type == InverseMassType::MatrixfreeOperator)
271 else if(data.implementation_type == InverseMassType::ElementwiseKrylovSolver or
272 data.implementation_type == InverseMassType::BlockMatrices)
278 block_jacobi_preconditioner->update();
280 else if(data.implementation_type == InverseMassType::GlobalKrylovSolver)
282 global_preconditioner->update();
286 AssertThrow(
false, dealii::ExcMessage(
"The specified InverseMassType is not implemented."));
294 if(data.implementation_type == InverseMassType::GlobalKrylovSolver)
296 AssertThrow(global_solver.get() != 0,
297 dealii::ExcMessage(
"Global mass solver has not been initialized."));
298 global_solver->solve(dst, src);
302 dst.zero_out_ghost_values();
304 if(data.implementation_type == InverseMassType::MatrixfreeOperator)
306 matrix_free->cell_loop(&This::cell_loop_matrix_free_operator,
this, dst, src);
310 AssertThrow(block_jacobi_preconditioner.get() != 0,
312 "Cell-wise iterative/direct block-Jacobi solver has not been initialized."));
313 block_jacobi_preconditioner->vmult(dst, src);
320 apply_scale(VectorType & dst,
double const scaling_factor, VectorType
const & src)
const
322 if(data.implementation_type == InverseMassType::MatrixfreeOperator)
328 dst.zero_out_ghost_values();
330 matrix_free->cell_loop(
331 &This::cell_loop_matrix_free_operator,
336 [&](
const unsigned int start_range,
const unsigned int end_range) {
337 for(
unsigned int i = start_range; i < end_range; ++i)
338 dst.local_element(i) *= scaling_factor;
345 dst *= scaling_factor;
352 cell_loop_matrix_free_operator(dealii::MatrixFree<dim, Number>
const &,
354 VectorType
const & src,
355 Range
const & cell_range)
const
357 Integrator integrator(*matrix_free, dof_index, quad_index);
358 InverseMassAsMatrixFreeOperator inverse_mass(integrator);
360 if(coefficient_is_variable)
362 for(
unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
364 integrator.reinit(cell);
365 integrator.read_dof_values(src, 0);
367 dealii::AlignedVector<dealii::VectorizedArray<Number>> inverse_JxW_times_coefficient(
368 integrator.n_q_points);
369 inverse_mass.fill_inverse_JxW_values(inverse_JxW_times_coefficient);
371 if(consider_inverse_coefficient)
376 for(
unsigned int q = 0; q < integrator.n_q_points; ++q)
378 inverse_JxW_times_coefficient[q] *=
379 this->variable_coefficients->get_coefficient_cell(cell, q);
387 for(
unsigned int q = 0; q < integrator.n_q_points; ++q)
389 inverse_JxW_times_coefficient[q] /=
390 this->variable_coefficients->get_coefficient_cell(cell, q);
394 inverse_mass.apply(inverse_JxW_times_coefficient,
396 integrator.begin_dof_values(),
397 integrator.begin_dof_values());
399 integrator.set_dof_values(dst, 0);
404 for(
unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
406 integrator.reinit(cell);
407 integrator.read_dof_values(src, 0);
409 inverse_mass.apply(integrator.begin_dof_values(), integrator.begin_dof_values());
411 integrator.set_dof_values(dst, 0);
416 dealii::MatrixFree<dim, Number>
const * matrix_free;
418 unsigned int dof_index, quad_index;
420 InverseMassParameters data;
423 bool coefficient_is_variable;
424 bool consider_inverse_coefficient;
426 VariableCoefficients<dealii::VectorizedArray<Number>>
const * variable_coefficients;
430 std::shared_ptr<PreconditionerBase<Number>> global_preconditioner;
431 std::shared_ptr<Krylov::SolverBase<VectorType>> global_solver;
437 std::shared_ptr<BlockJacobiPreconditioner<MassOperator<dim, n_components, Number>>>
438 block_jacobi_preconditioner;
441 MassOperator<dim, n_components, Number> mass_operator;