ryujin 2.5.0 revision c64172211cf21344cb9b9b0f3a1266aa6400d392
Loading...
Searching...
No Matches
parabolic_module.template.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2026 by the ryujin authors
4//
5
6#pragma once
7
8#include <ryujin/base/computing_timer.h>
9#include <ryujin/base/convenience_macros.h>
10#include <ryujin/base/loop.h>
11#include <ryujin/base/simd.h>
12#include <ryujin/euler_poisson/laplace_operator.h>
13#include <ryujin/euler_poisson/parabolic_module.h>
14
15#include <deal.II/dofs/dof_tools.h>
16#include <deal.II/lac/linear_operator.h>
17#include <deal.II/lac/precondition.h>
18#include <deal.II/lac/solver_cg.h>
19#include <deal.II/matrix_free/fe_evaluation.h>
20#include <deal.II/numerics/vector_tools.h>
21#include <deal.II/numerics/vector_tools.templates.h>
22
23
24namespace ryujin
25{
26 namespace EulerPoisson
27 {
28 using namespace dealii;
29
30 template <typename HyperbolicDescription, int dim, typename Number>
31 ParabolicModule<HyperbolicDescription, dim, Number>::ParabolicModule(
32 const MPIEnsemble &mpi_ensemble,
33 const OfflineData<dim, Number> &offline_data,
34 const HyperbolicSystem &hyperbolic_system,
35 const ParabolicSystem &parabolic_system,
36 const InitialValues<HyperbolicDescription, dim, Number> &initial_values,
37 const std::string &subsection)
38 : Interface::ParabolicModule<dim, Number>(subsection)
39 , mpi_ensemble_(mpi_ensemble)
40 , hyperbolic_system_(&hyperbolic_system)
41 , parabolic_system_(&parabolic_system)
42 , offline_data_(&offline_data)
43 , initial_values_(&initial_values)
44 , cycle_(0)
45 , n_iterations_gauss_(0)
46 , n_iterations_step_(0)
47 , potential_initialized_(false)
48 , t_background_density_(std::numeric_limits<Number>::lowest())
49 , t_magnetic_field_(std::numeric_limits<Number>::lowest())
50 {
51 gauss_law_restart_strategy_ = GaussLawRestartStrategy::no_restart;
52 this->add_parameter(
53 "gauss law restart strategy",
54 gauss_law_restart_strategy_,
55 "Strategy used when restarting the gauss law. Options are "
56 "\'no restart\', \'full restart\', \'correction\', "
57 "\'static no restart\', and \'static full restart\'.");
58
59 gmg_max_iter_ = 15;
60 this->add_parameter("multigrid - max iter",
61 gmg_max_iter_,
62 "Maximal number of CG iterations with GMG smoother");
63
64 gmg_smoother_range_ = 8.;
65 this->add_parameter("multigrid - chebyshev range",
66 gmg_smoother_range_,
67 "Chebyshev smoother: eigenvalue range parameter");
68
69 gmg_smoother_max_eig_ = 2.0;
70 this->add_parameter("multigrid - chebyshev max eig",
71 gmg_smoother_max_eig_,
72 "Chebyshev smoother: maximal eigenvalue");
73
74 gmg_smoother_degree_ = 3;
75 this->add_parameter("multigrid - chebyshev degree",
76 gmg_smoother_degree_,
77 "Chebyshev smoother: degree");
78
79 gmg_smoother_n_cg_iter_ = 10;
80 this->add_parameter(
81 "multigrid - chebyshev cg iter",
82 gmg_smoother_n_cg_iter_,
83 "Chebyshev smoother: number of CG iterations to approximate "
84 "eigenvalue");
85
86 gmg_min_level_ = 0;
87 this->add_parameter(
88 "multigrid - min level",
89 gmg_min_level_,
90 "Minimal mesh level to be visited in the geometric multigrid "
91 "cycle where the coarse grid solver (Chebyshev) is called");
92
93 tolerance_ = Number(1.0e-12);
94 this->add_parameter(
95 "tolerance", tolerance_, "Tolerance for linear solvers");
96
97 tolerance_linfty_norm_ = false;
98 this->add_parameter(
99 "tolerance linfty norm",
100 tolerance_linfty_norm_,
101 "Use the l_infty norm instead of the l_2 norm for the "
102 "stopping criterion");
103
104 electrostatic_configuration_list_ =
105 ElectrostaticConfigurationLibrary<dim, Number>::create(
106 parabolic_system_->subsection());
107
108 const auto populate = [this]() {
109 selected_electrostatic_configuration_ =
110 &ElectrostaticConfigurationLibrary<dim, Number>::select(
111 electrostatic_configuration_list_,
112 parabolic_system_->electrostatic_configuration());
113 };
114
115 this->parse_parameters_call_back.connect(populate);
116 populate();
117 }
118
119
120 template <typename HyperbolicDescription, int dim, typename Number>
121 void ParabolicModule<HyperbolicDescription, dim, Number>::prepare()
122 {
123#ifdef DEBUG_OUTPUT
124 std::cout << "ParabolicModule<dim, Number>::prepare()" << std::endl;
125#endif
126 /*
127 * The cycle_ variabe is only used for gmg reinitialization, simply
128 * reset it to zero on prepare().
129 */
130 cycle_ = 0;
131
132 const auto &discretization = offline_data_->discretization();
133 AssertThrow(discretization.ansatz() == Ansatz::dg_q1 ||
134 discretization.ansatz() == Ansatz::cg_q1,
135 dealii::ExcMessage("The Euler-Poisson module currently only "
136 "supports cG/dg Q1 finite elements."));
137
138 AssertThrow(!offline_data_->dof_handler().has_hp_capabilities(),
139 dealii::ExcMessage(
140 "The Euler-Poisson module currently does not support "
141 "DoFHandlers set up with hp capabilities."));
142
143 potential_initialized_ = false;
144
145 /*
146 * (Re)initialize matrix free object:
147 */
148
149 typename MatrixFree<dim, Number>::AdditionalData additional_data;
150 additional_data.tasks_parallel_scheme =
151 MatrixFree<dim, Number>::AdditionalData::none;
152
153 // First index CG, second index hyperbolic ansatz
154 std::vector<const dealii::DoFHandler<dim> *> dof_handlers = {
155 &offline_data_->dof_handler_cg(), &offline_data_->dof_handler()};
156
157 create_constraints();
158 std::vector<const dealii::AffineConstraints<Number> *>
159 affine_constraints = {&affine_constraints_potential_,
160 &offline_data_->affine_constraints()};
161
162 // First index full quadrature, second index lumped quadrature
163 std::vector<dealii::Quadrature<1>> quadratures = {
164 discretization.quadrature_1d()[0],
165 discretization.nodal_quadrature_1d()[0]};
166
167 matrix_free_.reinit(discretization.mapping(),
168 dof_handlers,
169 affine_constraints,
170 quadratures,
171 additional_data);
172
173 /*
174 * (Re)initialize operators and preconditioners:
175 */
176
177 laplace_operator_.initialize(matrix_free_);
178 laplace_operator_.compute_diagonal(diagonal_preconditioner_);
179 update_operator_.initialize(matrix_free_, density_, magnetic_field_);
180
181 typename decltype(multigrid_preconditioner_)::MultigridParameters
182 parameters{gmg_max_iter_,
183 gmg_smoother_range_,
184 gmg_smoother_max_eig_,
185 gmg_smoother_degree_,
186 gmg_smoother_n_cg_iter_,
187 gmg_min_level_,
188 tolerance_};
189
190 multigrid_preconditioner_.initialize(
191 *offline_data_,
192 selected_electrostatic_configuration_->dirichlet_boundaries(),
193 parameters);
194
195 /*
196 * (Re)initialize auxiliary vectors:
197 */
198
199 const auto &potential_partitioner =
200 matrix_free_.get_dof_info(0).vector_partitioner;
201 potential_rhs_.reinit(potential_partitioner);
202
203 const auto &scalar_partitioner =
204 matrix_free_.get_dof_info(1).vector_partitioner;
205 density_.reinit(scalar_partitioner);
206 background_density_.reinit(scalar_partitioner);
207
208 magnetic_field_.reinit(dim == 2 ? 1 : dim);
209 for (unsigned int i = 0; i < magnetic_field_.n_blocks(); ++i)
210 magnetic_field_.block(i).reinit(scalar_partitioner);
211
212 velocity_rhs_.reinit(dim);
213 for (unsigned int i = 0; i < dim; ++i)
214 velocity_rhs_.block(i).reinit(scalar_partitioner);
215
216 /*
217 * Populate background fields. The vectors have just been
218 * reinitialized, so invalidate the times they were populated for:
219 */
220
221 t_background_density_ = std::numeric_limits<Number>::lowest();
222 t_magnetic_field_ = std::numeric_limits<Number>::lowest();
223
224 if (!selected_electrostatic_configuration_->is_time_dependent()) {
225 update_background_density(Number(0.));
226 update_magnetic_field(Number(0.));
227 }
228 }
229
230
231 template <typename HyperbolicDescription, int dim, typename Number>
232 void
233 ParabolicModule<HyperbolicDescription, dim, Number>::reinit_state_vector(
234 StateVector<Number> &state_vector) const
235 {
236#ifdef DEBUG_OUTPUT
237 std::cout << "ParabolicModule<dim, Number>::reinit_state_vector()"
238 << std::endl;
239#endif
240
241 auto &[U, Q, U_dynamic, Q_dynamic] = state_vector;
242 U_dynamic.resize(1);
243
244 const auto &partitioner = matrix_free_.get_dof_info(0).vector_partitioner;
245 U_dynamic[0].reinit_with_scalar_partitioner(partitioner, 1);
246 U_dynamic[0].deal_ii_vector() = 0.;
247 }
248
249
250 template <typename HyperbolicDescription, int dim, typename Number>
251 void
252 ParabolicModule<HyperbolicDescription, dim, Number>::prepare_state_vector(
253 StateVector<Number> &state_vector, Number t) const
254 {
255#ifdef DEBUG_OUTPUT
256 std::cout << "ParabolicModule<dim, Number>::prepare_state_vector()"
257 << std::endl;
258#endif
259
260 AssertThrow(gauss_law_restart_strategy_ !=
261 GaussLawRestartStrategy::correction,
262 dealii::ExcNotImplemented());
263
264 /*
265 * Check whether the potential has been populated by a solution transfer:
266 */
267
268 if (!potential_initialized_) {
269 copy_to_memory_space<dealii::MemorySpace::Host>(
270 std::as_const(state_vector), StateVectorPart::dynamic_state);
271 const auto &phi = std::as_const(state_vector).dynamic_state[0];
272 potential_initialized_ = !phi.deal_ii_vector().all_zero();
273 }
274
275 /*
276 * We (re)compute the potential (and enforce the magnetic drift
277 * velocity) if it is not initialized, and recompute it on every
278 * step if the restart strategy is set to full_restart or
279 * static_full_restart.
280 */
281
282 if (!potential_initialized_ ||
283 (gauss_law_restart_strategy_ ==
284 GaussLawRestartStrategy::full_restart) ||
285 (gauss_law_restart_strategy_ ==
286 GaussLawRestartStrategy::static_full_restart)) {
287
288 compute_potential(t, state_vector);
289
290 if (!potential_initialized_ &&
291 parabolic_system_->magnetic_drift_limit())
292 enforce_magnetic_drift_velocity(state_vector);
293 potential_initialized_ = true;
294 }
295 }
296
297
298 template <typename HyperbolicDescription, int dim, typename Number>
299 Number
300 ParabolicModule<HyperbolicDescription, dim, Number>::backward_euler_step(
301 const StateVector<Number> &old_state_vector,
302 StageVectors stage_state_vectors,
303 StageWeights stage_weights,
304 StateVector<Number> &new_state_vector,
305 Number t,
306 Number tau,
307 Number /*tau_max*/) const
308 {
309 AssertThrow(tau != Number(0.),
310 dealii::ExcMessage("A parabolic module cannot choose the "
311 "time step size: tau must be nonzero."));
312
313 switch (stage_state_vectors.size()) {
314 case 0:
315 backward_euler_step<0>(old_state_vector,
316 stage_state_vectors,
317 stage_weights,
318 new_state_vector,
319 t,
320 tau);
321 break;
322 case 1:
323 backward_euler_step<1>(old_state_vector,
324 stage_state_vectors,
325 stage_weights,
326 new_state_vector,
327 t,
328 tau);
329 break;
330 case 2:
331 backward_euler_step<2>(old_state_vector,
332 stage_state_vectors,
333 stage_weights,
334 new_state_vector,
335 t,
336 tau);
337 break;
338 case 3:
339 backward_euler_step<3>(old_state_vector,
340 stage_state_vectors,
341 stage_weights,
342 new_state_vector,
343 t,
344 tau);
345 break;
346 case 4:
347 backward_euler_step<4>(old_state_vector,
348 stage_state_vectors,
349 stage_weights,
350 new_state_vector,
351 t,
352 tau);
353 break;
354 default:
355 AssertThrow(false, dealii::ExcMessage("Unsupported number of stages."));
356 }
357
358 return tau;
359 }
360
361
362 template <typename HyperbolicDescription, int dim, typename Number>
363 template <int stages>
364 void
365 ParabolicModule<HyperbolicDescription, dim, Number>::backward_euler_step(
366 const StateVector<Number> &old_state_vector,
367 StageVectors /*stage_state_vectors*/,
368 StageWeights /*stage_weights*/,
369 StateVector<Number> &new_state_vector,
370 Number t,
371 Number tau) const
372 {
373 step(old_state_vector,
374 new_state_vector,
375 t,
376 tau,
377 /*crank_nicolson_extrapolation = */ false);
378 }
379
380
381 template <typename HyperbolicDescription, int dim, typename Number>
382 Number
383 ParabolicModule<HyperbolicDescription, dim, Number>::crank_nicolson_step(
384 const StateVector<Number> &old_state_vector,
385 StateVector<Number> &new_state_vector,
386 Number t,
387 Number tau,
388 Number /*tau_max*/) const
389 {
390 AssertThrow(tau != Number(0.),
391 dealii::ExcMessage("A parabolic module cannot choose the "
392 "time step size: tau must be nonzero."));
393
394 try {
395 /* Backward Euler step to half time step, and extrapolate: */
396
397 step(old_state_vector,
398 new_state_vector,
399 t,
400 tau / Number(2.),
401 /*crank_nicolson_extrapolation = */ true);
402
403 } catch (Correction) {
404
405 /*
406 * Under very rare circumstances we might fail to perform a Crank
407 * Nicolson step because the extrapolation step produced
408 * inadmissible states. We could correct the update now by
409 * performing a limiting step (either convex limiting, or flux
410 * corrected transport)... but *meh*, just perform a backward Euler
411 * step:
412 */
413 step(old_state_vector,
414 new_state_vector,
415 t,
416 tau,
417 /*crank_nicolson_extrapolation = */ false);
418 }
419
420 return tau;
421 }
422
423
424 template <typename HyperbolicDescription, int dim, typename Number>
425 void ParabolicModule<HyperbolicDescription, dim, Number>::
426 print_solver_statistics(std::ostream &output) const
427 {
428 output << " [ " << std::setprecision(2) << std::fixed //
429 << n_iterations_gauss_ << " GMG gauss -- " //
430 << n_iterations_step_ << " GMG step ]" << std::endl;
431 }
432
433
434 template <typename HyperbolicDescription, int dim, typename Number>
435 void
436 ParabolicModule<HyperbolicDescription, dim, Number>::create_constraints()
437 {
438#ifdef DEBUG_OUTPUT
439 std::cout << "ParabolicModule<dim, Number>::create_constraints()"
440 << std::endl;
441#endif
442
443 const auto &discretization = offline_data_->discretization();
444 const auto &dof_handler = offline_data_->dof_handler_cg();
445
446 affine_constraints_potential_.clear();
447
448 const auto locally_relevant =
449 DoFTools::extract_locally_relevant_dofs(dof_handler);
450
451 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
452 affine_constraints_potential_.reinit(locally_owned, locally_relevant);
453
454 DoFTools::make_hanging_node_constraints(offline_data_->dof_handler_cg(),
455 affine_constraints_potential_);
456
457 /*
458 * Enforce periodic boundary conditions. We assume that the mesh is in
459 * "normal configuration."
460 */
461
462 const auto &periodic_faces =
463 discretization.triangulation().get_periodic_face_map();
464
465 for (const auto &[left, value] : periodic_faces) {
466 const auto &[right, orientation] = value;
467
468 typename DoFHandler<dim>::cell_iterator dof_cell_left(
469 &left.first->get_triangulation(),
470 left.first->level(),
471 left.first->index(),
472 &dof_handler);
473
474 typename DoFHandler<dim>::cell_iterator dof_cell_right(
475 &right.first->get_triangulation(),
476 right.first->level(),
477 right.first->index(),
478 &dof_handler);
479
480 if constexpr (std::is_same_v<Number, double>) {
481 DoFTools::make_periodicity_constraints(
482 dof_cell_left->face(left.second),
483 dof_cell_right->face(right.second),
484 affine_constraints_potential_,
485 ComponentMask(),
486 orientation);
487 } else {
488 AssertThrow(false, dealii::ExcNotImplemented());
489 __builtin_trap();
490 }
491 }
492
493 for (const auto &it :
494 selected_electrostatic_configuration_->dirichlet_boundaries())
495 DoFTools::make_zero_boundary_constraints(
496 offline_data_->dof_handler_cg(), it, affine_constraints_potential_);
497
498 affine_constraints_potential_.close();
499 }
500
501
502 template <typename HyperbolicDescription, int dim, typename Number>
503 void ParabolicModule<HyperbolicDescription, dim, Number>::
504 update_background_density(const Number t) const
505 {
506#ifdef DEBUG_OUTPUT
507 std::cout << "ParabolicModule<dim, Number>::update_background_density()"
508 << std::endl;
509#endif
510
511 /*
512 * Skip updating the background density if t > 0 and if the fields
513 * are time independent:
514 */
515 if (!selected_electrostatic_configuration_->is_time_dependent() &&
516 (t > Number(0.)))
517 return;
518
519 /*
520 * Skip updating if we have already populated the background density
521 * for the chosen time t.
522 */
523 if (std::abs(t_background_density_ - t) < 1.e-12)
524 return;
525
526#ifdef DEBUG_OUTPUT
527 std::cout << " updating to t = " << t << std::endl;
528#endif
529
530 ComputingTimer::Scope scope("time step [X] - interpolate data vectors");
531
532 const auto &discretization = offline_data_->discretization();
533 background_density_.zero_out_ghost_values();
534 dealii::VectorTools::interpolate(
535 discretization.mapping(),
536 offline_data_->dof_handler(),
537 dealii::ScalarFunctionFromFunctionObject<dim, Number>(
538 [&](const dealii::Point<dim> &p) {
539 return selected_electrostatic_configuration_
540 ->background_density(p, t);
541 }),
542 background_density_);
543 background_density_.update_ghost_values();
544
545 t_background_density_ = t;
546 }
547
548
549 template <typename HyperbolicDescription, int dim, typename Number>
550 void
551 ParabolicModule<HyperbolicDescription, dim, Number>::update_magnetic_field(
552 const Number t) const
553 {
554#ifdef DEBUG_OUTPUT
555 std::cout << "ParabolicModule<dim, Number>::update_magnetic_field()"
556 << std::endl;
557#endif
558
559 /*
560 * Skip updating the background density if t > 0 and if the fields
561 * are time independent:
562 */
563 if (!selected_electrostatic_configuration_->is_time_dependent() &&
564 (t > Number(0.)))
565 return;
566
567 /*
568 * Skip updating if we have already populated the background density
569 * for the chosen time t.
570 */
571 if (std::abs(t_magnetic_field_ - t) < 1.e-12)
572 return;
573
574#ifdef DEBUG_OUTPUT
575 std::cout << " updating to t = " << t << std::endl;
576#endif
577
578 ComputingTimer::Scope scope("time step [X] - interpolate data vectors");
579
580 const auto &discretization = offline_data_->discretization();
581 for (unsigned int k = 0; k < (dim == 2 ? 1 : dim); ++k) {
582 magnetic_field_.block(k).zero_out_ghost_values();
583 dealii::VectorTools::interpolate(
584 discretization.mapping(),
585 offline_data_->dof_handler(),
586 to_function<dim, Number>(
587 [&](const dealii::Point<dim> &p) {
588 return selected_electrostatic_configuration_->magnetic_field(
589 p, t);
590 },
591 k),
592 magnetic_field_.block(k));
593 }
594 magnetic_field_.update_ghost_values();
595
596 t_magnetic_field_ = t;
597 }
598
599
600 template <typename HyperbolicDescription, int dim, typename Number>
601 void ParabolicModule<HyperbolicDescription, dim, Number>::compute_potential(
602 const Number t, StateVector<Number> &state_vector) const
603 {
604#ifdef DEBUG_OUTPUT
605 std::cout << "ParabolicModule<dim, Number>::compute_potential()"
606 << std::endl;
607#endif
608 using HostSpace = dealii::MemorySpace::Host;
609
610 copy_to_memory_space<HostSpace>(std::as_const(state_vector),
611 StateVectorPart::state);
612 move_to_memory_space<HostSpace>(state_vector,
613 StateVectorPart::dynamic_state);
614
615 const auto U_view =
616 std::as_const(state_vector).state.template view<problem_dimension>();
617 auto &U_dynamic = state_vector.dynamic_state;
618 auto &potential = U_dynamic[0].deal_ii_vector();
619
620 const unsigned int n_owned = offline_data_->n_locally_owned();
621
622 constexpr unsigned int order_fe = 1;
623 constexpr unsigned int order_quad = 2;
624
625 /*
626 * -----------------------------------------------------------------------
627 * Step 1a: build right hand side for Gauss law
628 * -----------------------------------------------------------------------
629 */
630
631 ComputingTimer::Scope scope("time step [P] 1 - enforce Gauss law");
632
633 update_background_density(t);
634
635 const auto body_copy = [&](auto sentinel, unsigned int i) {
636 using T = decltype(sentinel);
637 const auto view = hyperbolic_system_->template view<dim, T>();
638 const auto U_i = U_view.template read_tensor<T>(i);
639 const auto rho_i = view.density(U_i);
640 write_entry<T>(density_, rho_i, i);
641 };
642
643 cpu_simd_loop<Number>(
644 "time_step_parabolic_1a", body_copy, 0, n_owned, n_owned);
645
646 density_.update_ghost_values();
647
648 const auto body_matrix_free = [this](const auto &data,
649 auto &dst,
650 const auto &src,
651 const auto range) {
652 FEEvaluation<dim, order_fe, order_quad, /*components*/ 1, Number>
653 fee_potential(data, /*CG*/ 0, /*lumped quadrature*/ 1);
654 FEEvaluation<dim, order_fe, order_quad, /*components*/ 1, Number>
655 fee_density(data, /*hyperbolic*/ 1, /*lumped quadrature*/ 1);
656 FEEvaluation<dim, order_fe, order_quad, /*components*/ 1, Number>
657 fee_background(data, /*hyperbolic*/ 1, /*lumped quadrature*/ 1);
658
659
660 const Number alpha = parabolic_system_->alpha();
661
662 for (unsigned int cell = range.first; cell < range.second; ++cell) {
663 fee_potential.reinit(cell);
664 fee_density.reinit(cell);
665 fee_background.reinit(cell);
666
667 fee_density.gather_evaluate(src, dealii::EvaluationFlags::values);
668 fee_background.gather_evaluate(background_density_,
669 dealii::EvaluationFlags::values);
670
671 for (unsigned int q = 0; q < fee_potential.n_q_points; ++q) {
672 const auto density_q = fee_density.get_value(q);
673 const auto background_q = fee_background.get_value(q);
674
675 const auto value = alpha * (density_q + background_q);
676 fee_potential.submit_value(value, q);
677 }
678 fee_potential.integrate_scatter(dealii::EvaluationFlags::values, dst);
679 }
680 };
681
682 matrix_free_.template cell_loop<ScalarHostVector<Number>,
683 ScalarHostVector<Number>>(
684 body_matrix_free,
685 potential_rhs_,
686 density_,
687 /*zero destination*/ true);
688
689 /*
690 * -----------------------------------------------------------------------
691 * Step 1b: solve Poisson problem
692 * -----------------------------------------------------------------------
693 */
694
695 matrix_free_.get_affine_constraints(0).distribute(potential);
696 matrix_free_.get_affine_constraints(0).set_zero(potential_rhs_);
697
698 const auto tolerance =
699 (tolerance_linfty_norm_ ? potential_rhs_.linfty_norm()
700 : potential_rhs_.l2_norm()) *
701 tolerance_;
702
703 typename dealii::SolverCG<ScalarHostVector<Number>>::AdditionalData
704 solver_data;
705
706 try {
707 SolverControl solver_control(gmg_max_iter_, tolerance);
708 dealii::SolverCG<ScalarHostVector<Number>> solver(solver_control,
709 solver_data);
710 solver.solve(laplace_operator_,
711 potential,
712 potential_rhs_,
713 multigrid_preconditioner_);
714
715
716 if (potential_initialized_) {
717 /* update exponential moving average */
718 n_iterations_gauss_ =
719 0.9 * n_iterations_gauss_ + 0.1 * solver_control.last_step();
720 } else {
721 n_iterations_gauss_ = solver_control.last_step();
722 }
723
724 } catch (SolverControl::NoConvergence &) {
725 SolverControl solver_control(1000, tolerance);
726 dealii::SolverCG<ScalarHostVector<Number>> solver(solver_control,
727 solver_data);
728
729 solver.solve(laplace_operator_,
730 potential,
731 potential_rhs_,
732 diagonal_preconditioner_);
733
734 if (potential_initialized_) {
735 /* update exponential moving average */
736 n_iterations_gauss_ *= 0.9;
737 n_iterations_gauss_ +=
738 0.1 * gmg_max_iter_ + 0.1 * solver_control.last_step();
739 } else {
740 n_iterations_gauss_ = gmg_max_iter_ + solver_control.last_step();
741 }
742
743 /* update exponential moving average, counting also GMG iterations */
744 }
745
746 matrix_free_.get_affine_constraints(0).distribute(potential);
747 }
748
749
750 template <typename HyperbolicDescription, int dim, typename Number>
751 void ParabolicModule<HyperbolicDescription, dim, Number>::
752 enforce_magnetic_drift_velocity(StateVector<Number> &state_vector) const
753 {
754#ifdef DEBUG_OUTPUT
755 std::cout
756 << "ParabolicModule<dim, Number>::enforce_magnetic_drift_velocity()"
757 << std::endl;
758#endif
759
760 using HostSpace = dealii::MemorySpace::Host;
761
762 move_to_memory_space<HostSpace>(state_vector, StateVectorPart::state);
763 copy_to_memory_space<HostSpace>(std::as_const(state_vector),
764 StateVectorPart::dynamic_state);
765
766 const auto U_view = state_vector.state.template view<problem_dimension>();
767 const auto &U_dynamic = std::as_const(state_vector).dynamic_state;
768 const auto &potential = U_dynamic[0].deal_ii_vector();
769
770 const unsigned int n_owned = offline_data_->n_locally_owned();
771
772 const auto lumped_mass_matrix_inverse_view =
773 offline_data_->lumped_mass_matrix_inverse().template view<1>();
774
775 constexpr unsigned int order_fe = 1;
776 constexpr unsigned int order_quad = 2;
777
778 /*
779 * -----------------------------------------------------------------------
780 * Step 1c: enforce magnetic drift velocity
781 * -----------------------------------------------------------------------
782 */
783
784 update_magnetic_field(Number(0.));
785
786 /* Project gradient of potential into velocity space: */
787
788 const auto body_velocity =
789 [](const auto &data, auto &dst, const auto &src, const auto range) {
790 FEEvaluation<dim, order_fe, order_quad, /*components*/ 1, Number>
791 fee_pot(data, /*CG*/ 0, /*lumped quadrature*/ 1);
792 FEEvaluation<dim, order_fe, order_quad, /*components*/ dim, Number>
793 fee_vel(data, /*hyperbolic*/ 1, /*lumped quadrature*/ 1);
794
795 for (unsigned int cell = range.first; cell < range.second; ++cell) {
796 fee_pot.reinit(cell);
797 fee_vel.reinit(cell);
798
799 fee_pot.gather_evaluate(src, dealii::EvaluationFlags::gradients);
800 for (unsigned int q = 0; q < fee_pot.n_q_points; ++q) {
801 fee_vel.submit_value(fee_pot.get_gradient(q), q);
802 }
803 fee_vel.integrate_scatter(dealii::EvaluationFlags::values, dst);
804 }
805 };
806
807 matrix_free_
808 .template cell_loop<BlockHostVector, ScalarHostVector<Number>>(
809 body_velocity,
810 velocity_rhs_,
811 potential,
812 /*zero destination*/ true);
813
814 const auto body = [&](auto sentinel, unsigned int i) {
815 using T = decltype(sentinel);
816 const auto view = hyperbolic_system_->template view<dim, T>();
817
818 const auto m_i_inv =
819 lumped_mass_matrix_inverse_view.template read_entry<T>(i);
820
821 auto U_i = U_view.template read_tensor<T>(i);
822 const auto rho_i = view.density(U_i);
823 const auto m_i = view.momentum(U_i);
824 const auto v_i = m_i / rho_i;
825
826 dealii::Tensor<1, (dim == 2 ? 1 : dim), T> magnetic_field;
827 for (unsigned int d = 0; d < (dim == 2 ? 1 : dim); ++d)
828 magnetic_field[d] = read_entry<T>(magnetic_field_.block(d), i);
829
830 dealii::Tensor<1, dim, T> grad_phi;
831 for (unsigned int d = 0; d < dim; ++d)
832 grad_phi[d] = m_i_inv * read_entry<T>(velocity_rhs_.block(d), i);
833
834 auto new_v_i = v_i;
835
836 if constexpr (dim == 2) {
837 new_v_i = -magnetic_field[0] * cross_product_2d(grad_phi) /
838 magnetic_field.norm_square();
839
840 } else if constexpr (dim == 3) {
841 new_v_i = -cross_product_3d(grad_phi, magnetic_field) /
842 magnetic_field.norm_square();
843 }
844
845 for (unsigned int d = 0; d < dim; ++d)
846 U_i[1 + d] = rho_i * new_v_i[d];
847
848 /* Update the total energy accordingly: */
849 if constexpr (view.have_energy_equation)
850 U_i[1 + dim] +=
851 Number(0.5) * rho_i * (new_v_i.norm_square() - v_i.norm_square());
852
853 U_view.template write_tensor<T>(U_i, i);
854 };
855
856 cpu_simd_loop<Number>(
857 "time_step_parabolic_1c", body, 0, n_owned, n_owned);
858 }
859
860
861 template <typename HyperbolicDescription, int dim, typename Number>
862 void ParabolicModule<HyperbolicDescription, dim, Number>::step(
863 const StateVector<Number> &old_state_vector,
864 StateVector<Number> &new_state_vector,
865 Number t,
866 Number tau [[maybe_unused]],
867 const bool crank_nicolson_extrapolation [[maybe_unused]]) const
868 {
869#ifdef DEBUG_OUTPUT
870 std::cout << "ParabolicModule<dim, Number>::step()" << std::endl;
871 std::cout << " perform time-step with tau = " << tau << std::endl;
872 if (crank_nicolson_extrapolation)
873 std::cout << " and extrapolate to t + 2 * tau" << std::endl;
874#endif
875
876 const Number alpha = parabolic_system_->alpha();
877
878 using HostSpace = dealii::MemorySpace::Host;
879
880 copy_to_memory_space<HostSpace>(old_state_vector,
881 StateVectorPart::state,
882 StateVectorPart::dynamic_state);
883 move_to_memory_space<HostSpace>(new_state_vector,
884 StateVectorPart::dynamic_state);
885
886 const auto &old_U = old_state_vector.state;
887 const auto old_U_view = old_U.template view<problem_dimension>();
888 const auto &old_U_dynamic = old_state_vector.dynamic_state;
889 const auto &old_potential = old_U_dynamic[0].deal_ii_vector();
890
891 auto &new_U = new_state_vector.state;
892 auto &new_U_dynamic = new_state_vector.dynamic_state;
893 auto &new_potential = new_U_dynamic[0].deal_ii_vector();
894
895 const unsigned int n_owned = offline_data_->n_locally_owned();
896
897 const auto lumped_mass_matrix_inverse_view =
898 offline_data_->lumped_mass_matrix_inverse().template view<1>();
899
900 constexpr unsigned int order_fe = 1;
901 constexpr unsigned int order_quad = 2;
902
903 /*
904 * Initialize the new potential with the old one:
905 */
906
907 new_potential = old_potential;
908
909 /*
910 * If the Gauss law restart strategy is "static full restart" or
911 * "static no restart", we skip updating the potential.
912 */
913 if ((gauss_law_restart_strategy_ !=
914 GaussLawRestartStrategy::static_no_restart) &&
915 (gauss_law_restart_strategy_ !=
916 GaussLawRestartStrategy::static_full_restart)) {
917
918 /*
919 * ---------------------------------------------------------------------
920 * Step 2a: build right hand side for potential update
921 *
922 * The right-hand side reads:
923 * (\nabla \varphi^n, \nabla \chi) +
924 * \tau \alpha \langle \rho^n B^{-1} v^n, \nabla \chi \rangle
925 *
926 * In case of a time-dependent background density, we add a term
927 * \theta \alpha \langle \rho_b^{n+1} - \rho_b^n, \chi \rangle to
928 * account for the time dependence. Here, t_{n+1} is the final time
929 * t_n + tau, or t_n + 2 * tau (in case of Crank Nicolson). This
930 * ensures that we are consistent with the Gauß law involution
931 * "-\Delta \varphi^{n+1} = \alpha \rho^{n+1}."
932 * ---------------------------------------------------------------------
933 */
934
935 ComputingTimer::Scope scope("time step [P] 2 - update potential");
936
937 /* Query the magnetic field at the time t + tau: */
938 update_magnetic_field(t + tau);
939
940 /*
941 * Write out density and assemble velocity part. We need density_
942 * to be set to the correct density for UpdateOperator::vmult()
943 */
944
945 const auto body_copy = [&](auto sentinel, unsigned int i) {
946 using T = decltype(sentinel);
947 const auto view = hyperbolic_system_->template view<dim, T>();
948
949 const auto U_i = old_U_view.template read_tensor<T>(i);
950 const auto rho_i = view.density(U_i);
951 const auto m_i = view.momentum(U_i);
952
953 dealii::Tensor<1, (dim == 2 ? 1 : dim), T> magnetic_field;
954 for (unsigned int d = 0; d < (dim == 2 ? 1 : dim); ++d)
955 magnetic_field[d] = read_entry<T>(magnetic_field_.block(d), i);
956
957 const auto velocity_rhs =
958 tau * alpha * rho_i *
959 apply_B_n_inverse(magnetic_field, tau, m_i / rho_i);
960
961 write_entry<T>(density_, rho_i, i);
962 for (unsigned int d = 0; d < dim; ++d)
963 write_entry<T>(velocity_rhs_.block(d), velocity_rhs[d], i);
964 };
965
966 cpu_simd_loop<Number>(
967 "time_step_parabolic_2a", body_copy, 0, n_owned, n_owned);
968
969 density_.update_ghost_values();
970
971 /* Apply Laplace operator to right hand side: */
972
973 const auto body_laplace = [](const auto &data,
974 auto &dst,
975 const auto &src,
976 const auto range) {
977 FEEvaluation<dim, order_fe, order_quad, /*components*/ 1, Number> fee(
978 data, /*CG*/ 0, /*full quadrature*/ 0);
979
980 for (unsigned int cell = range.first; cell < range.second; ++cell) {
981 fee.reinit(cell);
982 fee.gather_evaluate(src, dealii::EvaluationFlags::gradients);
983
984 for (unsigned int q = 0; q < fee.n_q_points; ++q) {
985 const auto grad_potential = fee.get_gradient(q);
986 fee.submit_gradient(grad_potential, q);
987 }
988 fee.integrate_scatter(dealii::EvaluationFlags::gradients, dst);
989 }
990 };
991
992 matrix_free_.template cell_loop<ScalarHostVector<Number>,
993 ScalarHostVector<Number>>(
994 body_laplace,
995 potential_rhs_,
996 old_potential,
997 /*zero destination*/ true);
998
999 /* Apply Velocity contribution to right hand side: */
1000
1001 const auto body_velocity = [](const auto &data,
1002 auto &dst,
1003 const auto &src,
1004 const auto range) {
1005 FEEvaluation<dim, order_fe, order_quad, /*components*/ 1, Number>
1006 fee_pot(data, /*CG*/ 0, /*lumped quadrature*/ 1);
1007 FEEvaluation<dim, order_fe, order_quad, /*components*/ dim, Number>
1008 fee_vel(data, /*hyperbolic*/ 1, /*lumped quadrature*/ 1);
1009
1010 for (unsigned int cell = range.first; cell < range.second; ++cell) {
1011 fee_pot.reinit(cell);
1012 fee_vel.reinit(cell);
1013
1014 fee_vel.gather_evaluate(src, dealii::EvaluationFlags::values);
1015
1016 for (unsigned int q = 0; q < fee_pot.n_q_points; ++q) {
1017 if constexpr (dim == 1) {
1018 decltype(fee_pot.get_gradient(q)) velocity_rhs;
1019 velocity_rhs[0] = fee_vel.get_value(q);
1020 fee_pot.submit_gradient(velocity_rhs, q);
1021 } else {
1022 fee_pot.submit_gradient(fee_vel.get_value(q), q);
1023 }
1024 }
1025 fee_pot.integrate_scatter(dealii::EvaluationFlags::gradients, dst);
1026 }
1027 };
1028
1029 matrix_free_
1030 .template cell_loop<ScalarHostVector<Number>, BlockHostVector>(
1031 body_velocity,
1032 potential_rhs_,
1033 velocity_rhs_,
1034 /*zero destination*/ false);
1035
1036 /* Time-dependent background density: */
1037
1038 if (selected_electrostatic_configuration_->is_time_dependent()) {
1039
1040 /*
1041 * Subtract background density at time t_n:
1042 */
1043
1044 update_background_density(t);
1045
1046 Number factor = (crank_nicolson_extrapolation ? -0.5 : -1.0) * alpha;
1047
1048 const auto body = [&factor](const auto &data,
1049 auto &dst,
1050 const auto &src,
1051 const auto range) {
1052 FEEvaluation<dim, order_fe, order_quad, /*components*/ 1, Number>
1053 fee_potential(data, /*CG*/ 0, /*lumped quadrature*/ 1);
1054 FEEvaluation<dim, order_fe, order_quad, /*components*/ 1, Number>
1055 fee_background(data, /*hyperbolic*/ 1, /*lumped quadrature*/ 1);
1056
1057 for (unsigned int cell = range.first; cell < range.second; ++cell) {
1058 fee_potential.reinit(cell);
1059 fee_background.reinit(cell);
1060 fee_background.gather_evaluate(src, EvaluationFlags::values);
1061
1062 for (unsigned int q = 0; q < fee_potential.n_q_points; ++q) {
1063 const auto background_q = fee_background.get_value(q);
1064 fee_potential.submit_value(factor * background_q, q);
1065 }
1066 fee_potential.integrate_scatter(EvaluationFlags::values, dst);
1067 }
1068 };
1069
1070 matrix_free_.template cell_loop<ScalarHostVector<Number>,
1071 ScalarHostVector<Number>>(
1072 body,
1073 potential_rhs_,
1074 background_density_,
1075 /*zero destination*/ false);
1076
1077 /*
1078 * Add background density at time t_{n+1}:
1079 */
1080
1081 update_background_density(
1082 t + (crank_nicolson_extrapolation ? 2. : 1.) * tau);
1083
1084 factor *= -1.;
1085
1086 matrix_free_.template cell_loop<ScalarHostVector<Number>,
1087 ScalarHostVector<Number>>(
1088 body,
1089 potential_rhs_,
1090 background_density_,
1091 /*zero destination*/ false);
1092 }
1093
1094 /*
1095 * ---------------------------------------------------------------------
1096 * Step 2b: solve modified poisson problem
1097 * ---------------------------------------------------------------------
1098 */
1099
1100 update_operator_.set_alpha(alpha);
1101 update_operator_.set_theta_tau(tau);
1102
1103 matrix_free_.get_affine_constraints(0).distribute(new_potential);
1104 matrix_free_.get_affine_constraints(0).set_zero(potential_rhs_);
1105
1106 const auto tolerance =
1107 (tolerance_linfty_norm_ ? potential_rhs_.linfty_norm()
1108 : potential_rhs_.l2_norm()) *
1109 tolerance_;
1110
1111 typename dealii::SolverCG<ScalarHostVector<Number>>::AdditionalData
1112 solver_data;
1113
1114 try {
1115 SolverControl solver_control(gmg_max_iter_, tolerance);
1116 dealii::SolverCG<ScalarHostVector<Number>> solver(solver_control,
1117 solver_data);
1118 solver.solve(update_operator_,
1119 new_potential,
1120 potential_rhs_,
1121 multigrid_preconditioner_);
1122
1123 /* update exponential moving average */
1124 n_iterations_step_ =
1125 0.9 * n_iterations_step_ + 0.1 * solver_control.last_step();
1126
1127 } catch (SolverControl::NoConvergence &) {
1128 SolverControl solver_control(1000, tolerance);
1129 dealii::SolverCG<ScalarHostVector<Number>> solver(solver_control,
1130 solver_data);
1131
1132 solver.solve(update_operator_,
1133 new_potential,
1134 potential_rhs_,
1135 diagonal_preconditioner_);
1136
1137 /* update exponential moving average, counting also GMG iterations */
1138 n_iterations_step_ *= 0.9;
1139 n_iterations_step_ +=
1140 0.1 * gmg_max_iter_ + 0.1 * solver_control.last_step();
1141 }
1142
1143 matrix_free_.get_affine_constraints(0).distribute(new_potential);
1144 }
1145
1146 /*
1147 * ---------------------------------------------------------------------
1148 * Step 2c: update velocity vector field; Crank-Nicolson extrapolation
1149 * ---------------------------------------------------------------------
1150 */
1151
1152 /* Project gradient of potential into velocity space: */
1153
1154 const auto body_velocity =
1155 [](const auto &data, auto &dst, const auto &src, const auto range) {
1156 FEEvaluation<dim, order_fe, order_quad, /*components*/ 1, Number>
1157 fee_pot(data, /*CG*/ 0, /*lumped quadrature*/ 1);
1158 FEEvaluation<dim, order_fe, order_quad, /*components*/ dim, Number>
1159 fee_vel(data, /*hyperbolic*/ 1, /*lumped quadrature*/ 1);
1160
1161 for (unsigned int cell = range.first; cell < range.second; ++cell) {
1162 fee_pot.reinit(cell);
1163 fee_vel.reinit(cell);
1164
1165 fee_pot.gather_evaluate(src, dealii::EvaluationFlags::gradients);
1166 for (unsigned int q = 0; q < fee_pot.n_q_points; ++q) {
1167 fee_vel.submit_value(fee_pot.get_gradient(q), q);
1168 }
1169 fee_vel.integrate_scatter(dealii::EvaluationFlags::values, dst);
1170 }
1171 };
1172
1173 matrix_free_
1174 .template cell_loop<BlockHostVector, ScalarHostVector<Number>>(
1175 body_velocity,
1176 velocity_rhs_,
1177 new_potential,
1178 /*zero destination*/ true);
1179
1180 /*
1181 * Now that we have written out the gradients, copy over the old
1182 * state vector and perform the Crank-Nicolson extrapolation step on
1183 * the potential:
1184 */
1185
1186 new_U = old_U;
1187 move_to_memory_space<HostSpace>(new_U); /* invalidates default */
1188 const auto new_U_view = new_U.template view<problem_dimension>();
1189
1190 if (crank_nicolson_extrapolation) {
1191 new_potential *= Number(2.);
1192 new_potential -= old_potential;
1193 }
1194
1195 /*
1196 * Update the momentum and total energy:
1197 */
1198
1199 const auto body = [&](auto sentinel, unsigned int i) {
1200 using T = decltype(sentinel);
1201 const auto view = hyperbolic_system_->template view<dim, T>();
1202
1203 const auto m_i_inv =
1204 lumped_mass_matrix_inverse_view.template read_entry<T>(i);
1205
1206 const auto old_U_i = old_U_view.template read_tensor<T>(i);
1207 const auto rho_i = view.density(old_U_i);
1208 const auto old_m_i = view.momentum(old_U_i);
1209 const auto old_v_i = old_m_i / rho_i;
1210
1211 dealii::Tensor<1, (dim == 2 ? 1 : dim), T> magnetic_field;
1212 for (unsigned int d = 0; d < (dim == 2 ? 1 : dim); ++d)
1213 magnetic_field[d] = read_entry<T>(magnetic_field_.block(d), i);
1214
1215 dealii::Tensor<1, dim, T> grad_phi;
1216 for (unsigned int d = 0; d < dim; ++d)
1217 grad_phi[d] = m_i_inv * read_entry<T>(velocity_rhs_.block(d), i);
1218
1219 auto new_v_i =
1220 apply_B_n_inverse(magnetic_field, tau, old_v_i - tau * grad_phi);
1221
1222 /* Perform an extrapolation step: */
1223 if (crank_nicolson_extrapolation)
1224 new_v_i = Number(2.) * new_v_i - old_v_i;
1225
1226 auto new_U_i = old_U_i;
1227 for (unsigned int d = 0; d < dim; ++d)
1228 new_U_i[1 + d] = rho_i * new_v_i[d];
1229
1230 /* Update the total energy accordingly: */
1231 if constexpr (view.have_energy_equation)
1232 new_U_i[1 + dim] += Number(0.5) * rho_i *
1233 (new_v_i.norm_square() - old_v_i.norm_square());
1234
1235 new_U_view.template write_tensor<T>(new_U_i, i);
1236 };
1237
1238 cpu_simd_loop<Number>(
1239 "time_step_parabolic_2c", body, 0, n_owned, n_owned);
1240 }
1241
1242 } // namespace EulerPoisson
1243} // namespace ryujin
ryujin::EulerPoisson::ParabolicModule::StageVectors
typename Interface::ParabolicModule< dim, Number >::StageVectors StageVectors
Definition parabolic_module.h:144
ryujin::EulerPoisson::ParabolicModule::ParabolicModule
ParabolicModule(const MPIEnsemble &mpi_ensemle, const OfflineData< dim, Number > &offline_data, const HyperbolicSystem &hyperbolic_system, const ParabolicSystem &parabolic_system, const InitialValues< HyperbolicDescription, dim, Number > &initial_values, const std::string &subsection="/ParabolicModule")
Definition parabolic_module.template.h:31
ryujin::EulerPoisson::ParabolicModule::StageWeights
typename Interface::ParabolicModule< dim, Number >::StageWeights StageWeights
Definition parabolic_module.h:147
ryujin::EulerPoisson::ParabolicModule::HyperbolicSystem
typename HyperbolicDescription::HyperbolicSystem HyperbolicSystem
Definition parabolic_module.h:149
ryujin::EulerPoisson::ParabolicModule::print_solver_statistics
void print_solver_statistics(std::ostream &output) const override
Definition parabolic_module.template.h:426
ryujin::EulerPoisson::ParabolicModule::crank_nicolson_step
Number crank_nicolson_step(const StateVector< Number > &old_state_vector, StateVector< Number > &new_state_vector, Number t, Number tau, Number tau_max=std::numeric_limits< Number >::max()) const override
Definition parabolic_module.template.h:383
ryujin::EulerPoisson::ParabolicModule::prepare_state_vector
void prepare_state_vector(StateVector< Number > &state_vector, Number t) const override
Definition parabolic_module.template.h:252
ryujin::EulerPoisson::ParabolicModule::backward_euler_step
Number backward_euler_step(const StateVector< Number > &old_state_vector, StageVectors stage_state_vectors, StageWeights stage_weights, StateVector< Number > &new_state_vector, Number t, Number tau, Number tau_max=std::numeric_limits< Number >::max()) const override
Definition parabolic_module.template.h:300
ryujin::EulerPoisson::ParabolicModule::reinit_state_vector
void reinit_state_vector(StateVector< Number > &state_vector) const override
Definition parabolic_module.template.h:233
ryujin::Library::select
static Base & select(const std::vector< std::unique_ptr< Base > > &objects, const std::string &name)
Definition library.h:189
ryujin::Library::create
static std::vector< std::unique_ptr< Base > > create(Args... args)
Definition library.h:148
computing_timer.h
convenience_macros.h
parabolic_module.h
ryujin::TaskPhase::step
@ step
ryujin::EulerPoisson::static_full_restart
@ static_full_restart
Definition parabolic_module.h:102
ryujin::EulerPoisson::full_restart
@ full_restart
Definition parabolic_module.h:69
ryujin::EulerPoisson::static_no_restart
@ static_no_restart
Definition parabolic_module.h:92
ryujin::EulerPoisson::correction
@ correction
Definition parabolic_module.h:82
ryujin::EulerPoisson::no_restart
@ no_restart
Definition parabolic_module.h:49
ryujin::StateVectorPart::dynamic_state
@ dynamic_state
ryujin::StateVectorPart::state
@ state
ryujin::Ansatz::cg_q1
@ cg_q1
ryujin::Ansatz::dg_q1
@ dg_q1
laplace_operator.h
ryujin::apply_B_n_inverse
DEAL_II_ALWAYS_INLINE dealii::Tensor< 1, dim, Number > apply_B_n_inverse(const dealii::Tensor< 1,(dim==2 ? 1 :dim), Number > &magnetic_field, const Number2 &theta_tau, const dealii::Tensor< 1, dim, Number > &velocity)
Definition laplace_operator.h:51
Generated by doxygen 1.9.8