DiFfRG
Loading...
Searching...
No Matches
ldg.hh
Go to the documentation of this file.
1#pragma once
2
3// standard library
4#include <memory>
5#include <sstream>
6
7// external libraries
8#include <deal.II/base/multithread_info.h>
9#include <deal.II/base/quadrature_lib.h>
10#include <deal.II/base/timer.h>
11#include <deal.II/dofs/dof_handler.h>
12#include <deal.II/dofs/dof_tools.h>
13#include <deal.II/fe/fe_interface_values.h>
14#include <deal.II/fe/fe_values.h>
15#include <deal.II/lac/block_sparse_matrix.h>
16#include <deal.II/lac/block_vector.h>
17#include <deal.II/lac/full_matrix.h>
18#include <deal.II/lac/vector.h>
19#include <deal.II/lac/vector_memory.h>
20#include <deal.II/meshworker/mesh_loop.h>
21#include <deal.II/numerics/matrix_tools.h>
22#include <deal.II/numerics/vector_tools.h>
23#include <tbb/tbb.h>
24
25// DiFfRG
26#include <DiFfRG/common/utils.hh>
27#include <DiFfRG/discretization/common/EoM.hh>
28#include <DiFfRG/discretization/common/abstract_assembler.hh>
29#include <DiFfRG/discretization/data/data_output.hh>
30
31namespace DiFfRG
32{
33 namespace LDG
34 {
35 using namespace dealii;
36 using std::array, std::vector, std::unique_ptr;
37
38 template <typename Discretization_, typename Model_>
39 class LDGAssemblerBase : public AbstractAssembler<typename Discretization_::VectorType,
40 typename Discretization_::SparseMatrixType, Discretization_::dim>
41 {
42 public:
43 using Discretization = Discretization_;
44 using Model = Model_;
45 using NumberType = typename Discretization::NumberType;
46 using VectorType = typename Discretization::VectorType;
47
48 using Components = typename Discretization::Components;
49 static constexpr uint dim = Discretization::dim;
50
51 LDGAssemblerBase(Discretization &discretization, Model &model, const JSONValue &json)
52 : discretization(discretization), model(model), fe(discretization.get_fe()),
53 dof_handler(discretization.get_dof_handler()), mapping(discretization.get_mapping()),
54 threads(json.get_uint("/discretization/threads")), batch_size(json.get_uint("/discretization/batch_size")),
55 EoM_cell(*(dof_handler.active_cell_iterators().end())),
56 old_EoM_cell(*(dof_handler.active_cell_iterators().end())),
57 EoM_abs_tol(json.get_double("/discretization/EoM_abs_tol")),
58 EoM_max_iter(json.get_uint("/discretization/EoM_max_iter"))
59 {
60 if (this->threads == 0) this->threads = dealii::MultithreadInfo::n_threads() / 2;
61 spdlog::get("log")->info("FEM: Using {} threads for assembly.", threads);
62 }
63
64 virtual IndexSet get_differential_indices() const override
65 {
66 ComponentMask component_mask(model.template differential_components<dim>());
67 return DoFTools::extract_dofs(dof_handler, component_mask);
68 }
69
70 const auto &get_discretization() const { return discretization; }
71 auto &get_discretization() { return discretization; }
72
73 virtual void reinit() override
74 {
75 // Boundary dofs
76 std::vector<IndexSet> component_boundary_dofs(Components::count_fe_functions());
77 for (uint c = 0; c < Components::count_fe_functions(); ++c) {
78 ComponentMask component_mask(Components::count_fe_functions(), false);
79 component_mask.set(c, true);
80 component_boundary_dofs[c] = DoFTools::extract_boundary_dofs(dof_handler, component_mask);
81 }
82 std::vector<std::vector<Point<dim>>> component_boundary_points(Components::count_fe_functions());
83 for (uint c = 0; c < Components::count_fe_functions(); ++c) {
84 component_boundary_points[c].resize(component_boundary_dofs[c].n_elements());
85 for (uint i = 0; i < component_boundary_dofs[c].n_elements(); ++i)
86 component_boundary_points[c][i] =
87 discretization.get_support_point(component_boundary_dofs[c].nth_index_in_set(i));
88 }
89
90 auto &constraints = discretization.get_constraints();
91 constraints.clear();
92 DoFTools::make_hanging_node_constraints(dof_handler, constraints);
93 model.affine_constraints(constraints, component_boundary_dofs, component_boundary_points);
94 constraints.close();
95 }
96
97 virtual void rebuild_jacobian_sparsity() = 0;
98
99 virtual void set_time(double t) override { model.set_time(t); }
100
101 virtual void refinement_indicator(Vector<double> & /*indicator*/, const VectorType & /*solution*/) = 0;
102
103 double average_time_variable_residual_assembly()
104 {
105 double t = 0.;
106 double n = timings_variable_residual.size();
107 for (const auto &t_ : timings_variable_residual)
108 t += t_ / n;
109 return t;
110 }
111 uint num_variable_residuals() const { return timings_variable_residual.size(); }
112
113 double average_time_variable_jacobian_assembly()
114 {
115 double t = 0.;
116 double n = timings_variable_jacobian.size();
117 for (const auto &t_ : timings_variable_jacobian)
118 t += t_ / n;
119 return t;
120 }
121 uint num_variable_jacobians() const { return timings_variable_jacobian.size(); }
122
123 protected:
124 Discretization &discretization;
125 Model &model;
126 const FiniteElement<dim> &fe;
127 const DoFHandler<dim> &dof_handler;
128 const Mapping<dim> &mapping;
129
130 uint threads;
131 uint batch_size;
132
133 mutable typename DoFHandler<dim>::cell_iterator EoM_cell;
134 typename DoFHandler<dim>::cell_iterator old_EoM_cell;
135 const double EoM_abs_tol;
136 const uint EoM_max_iter;
137 mutable Point<dim> EoM;
138 FullMatrix<NumberType> extractor_jacobian;
139 FullMatrix<NumberType> extractor_jacobian_u;
140 FullMatrix<NumberType> extractor_jacobian_du;
141 FullMatrix<NumberType> extractor_jacobian_ddu;
142 std::vector<types::global_dof_index> extractor_dof_indices;
143
144 std::vector<double> timings_variable_residual;
145 std::vector<double> timings_variable_jacobian;
146 };
147
148 namespace internal
149 {
154 template <typename Discretization> struct ScratchData {
155 static constexpr uint dim = Discretization::dim;
156 using NumberType = typename Discretization::NumberType;
157 static constexpr uint n_fe_subsystems = Discretization::Components::count_fe_subsystems();
158 using Iterator = typename DoFHandler<dim>::active_cell_iterator;
159 using t_Iterator = typename Triangulation<dim>::active_cell_iterator;
160
161 ScratchData(const Mapping<dim> &mapping, const vector<const DoFHandler<dim> *> &dofh,
162 const Quadrature<dim> &quadrature, const Quadrature<dim - 1> &quadrature_face,
163 const UpdateFlags update_flags = update_values | update_gradients | update_quadrature_points |
164 update_JxW_values,
165 const UpdateFlags interface_update_flags = update_values | update_gradients |
166 update_quadrature_points | update_JxW_values |
167 update_normal_vectors)
168 {
169 AssertThrow(dofh.size() >= n_fe_subsystems,
170 StandardExceptions::ExcDimensionMismatch(dofh.size(), n_fe_subsystems));
171
172 for (uint i = 0; i < n_fe_subsystems; ++i) {
173 const auto &fe = dofh[i]->get_fe();
174 fe_values[i] = unique_ptr<FEValues<dim>>(new FEValues<dim>(mapping, fe, quadrature, update_flags));
175 fe_interface_values[i] = unique_ptr<FEInterfaceValues<dim>>(
176 new FEInterfaceValues<dim>(mapping, fe, quadrature_face, interface_update_flags));
177 fe_boundary_values[i] = unique_ptr<FEFaceValues<dim>>(
178 new FEFaceValues<dim>(mapping, fe, quadrature_face, interface_update_flags));
179
180 n_components[i] = fe.n_components();
181 solution[i].resize(quadrature.size(), Vector<NumberType>(n_components[i]));
182 solution_interface[0][i].resize(quadrature_face.size(), Vector<NumberType>(n_components[i]));
183 solution_interface[1][i].resize(quadrature_face.size(), Vector<NumberType>(n_components[i]));
184
185 cell[i] = dofh[i]->begin_active();
186 ncell[i] = dofh[i]->begin_active();
187 }
188 solution_dot.resize(quadrature.size(), Vector<NumberType>(n_components[0]));
189 }
190
191 ScratchData(const ScratchData<Discretization> &scratch_data)
192 {
193 for (uint i = 0; i < n_fe_subsystems; ++i) {
194 const auto &old_fe = scratch_data.fe_values[i];
195 const auto &old_fe_i = scratch_data.fe_interface_values[i];
196 const auto &old_fe_b = scratch_data.fe_boundary_values[i];
197
198 fe_values[i] = unique_ptr<FEValues<dim>>(new FEValues<dim>(
199 old_fe->get_mapping(), old_fe->get_fe(), old_fe->get_quadrature(), old_fe->get_update_flags()));
200 fe_interface_values[i] = unique_ptr<FEInterfaceValues<dim>>(new FEInterfaceValues<dim>(
201 old_fe_i->get_mapping(), old_fe_i->get_fe(), old_fe_i->get_quadrature(), old_fe_i->get_update_flags()));
202 fe_boundary_values[i] = unique_ptr<FEFaceValues<dim>>(new FEFaceValues<dim>(
203 old_fe_b->get_mapping(), old_fe_b->get_fe(), old_fe_b->get_quadrature(), old_fe_b->get_update_flags()));
204
205 n_components[i] = scratch_data.n_components[i];
206 solution[i].resize(scratch_data.solution[i].size(), Vector<NumberType>(n_components[i]));
207 solution_interface[0][i].resize(scratch_data.solution_interface[0][i].size(),
208 Vector<NumberType>(n_components[i]));
209 solution_interface[1][i].resize(scratch_data.solution_interface[1][i].size(),
210 Vector<NumberType>(n_components[i]));
211
212 cell[i] = scratch_data.cell[i];
213 ncell[i] = scratch_data.ncell[i];
214 }
215 solution_dot.resize(scratch_data.solution_dot.size(), Vector<NumberType>(n_components[0]));
216 }
217
218 const auto &new_fe_values(const t_Iterator &t_cell)
219 {
220 for (uint i = 0; i < n_fe_subsystems; ++i) {
221 cell[i]->copy_from(*t_cell);
222 fe_values[i]->reinit(cell[i]);
223 }
224 return fe_values;
225 }
226 const auto &new_fe_interface_values(const t_Iterator &t_cell, uint f, uint sf, const t_Iterator &t_ncell,
227 uint nf, unsigned int nsf)
228 {
229 for (uint i = 0; i < n_fe_subsystems; ++i) {
230 cell[i]->copy_from(*t_cell);
231 ncell[i]->copy_from(*t_ncell);
232 fe_interface_values[i]->reinit(cell[i], f, sf, ncell[i], nf, nsf);
233 }
234 return fe_interface_values;
235 }
236 const auto &new_fe_boundary_values(const t_Iterator &t_cell, uint face_no)
237 {
238 for (uint i = 0; i < n_fe_subsystems; ++i) {
239 cell[i]->copy_from(*t_cell);
240 fe_boundary_values[i]->reinit(cell[i], face_no);
241 }
242 return fe_boundary_values;
243 }
244
245 array<uint, n_fe_subsystems> n_components;
246 array<Iterator, n_fe_subsystems> cell;
247 array<Iterator, n_fe_subsystems> ncell;
248
249 array<unique_ptr<FEValues<dim>>, n_fe_subsystems> fe_values;
250 array<unique_ptr<FEInterfaceValues<dim>>, n_fe_subsystems> fe_interface_values;
251 array<unique_ptr<FEFaceValues<dim>>, n_fe_subsystems> fe_boundary_values;
252
253 array<vector<Vector<NumberType>>, n_fe_subsystems> solution;
254 vector<Vector<NumberType>> solution_dot;
255 array<array<vector<Vector<NumberType>>, n_fe_subsystems>, 2> solution_interface;
256 };
257
258 template <typename NumberType> struct CopyData_R {
259 struct CopyDataFace_R {
260 Vector<NumberType> cell_residual;
261 std::vector<types::global_dof_index> joint_dof_indices;
262 };
263
264 Vector<NumberType> cell_residual;
265 Vector<NumberType> cell_mass;
266 std::vector<types::global_dof_index> local_dof_indices;
267 std::vector<CopyDataFace_R> face_data;
268
269 template <class Iterator> void reinit(const Iterator &cell, uint dofs_per_cell)
270 {
271 cell_residual.reinit(dofs_per_cell);
272 cell_mass.reinit(dofs_per_cell);
273 local_dof_indices.resize(dofs_per_cell);
274 cell->get_dof_indices(local_dof_indices);
275 }
276
277 template <int dim> CopyDataFace_R &new_face_data(const FEInterfaceValues<dim> &fe_iv)
278 {
279 face_data.emplace_back();
280 auto &copy_data_face = face_data.back();
281 copy_data_face.cell_residual.reinit(fe_iv.n_current_interface_dofs());
282 copy_data_face.joint_dof_indices = fe_iv.get_interface_dof_indices();
283 return copy_data_face;
284 }
285 };
286
287 template <typename NumberType> struct CopyData_J {
288 struct CopyDataFace_J {
289 FullMatrix<NumberType> cell_jacobian;
290 std::vector<types::global_dof_index> joint_dof_indices_from;
291 std::vector<types::global_dof_index> joint_dof_indices_to;
292 };
293
294 FullMatrix<NumberType> cell_jacobian;
295 std::vector<types::global_dof_index> local_dof_indices_from;
296 std::vector<types::global_dof_index> local_dof_indices_to;
297 std::vector<CopyDataFace_J> face_data;
298
299 template <class Iterator>
300 void reinit(const Iterator &cell_from, const Iterator &cell_to, uint dofs_per_cell_from, uint dofs_per_cell_to)
301 {
302 cell_jacobian.reinit(dofs_per_cell_to, dofs_per_cell_from);
303 local_dof_indices_from.resize(dofs_per_cell_from);
304 local_dof_indices_to.resize(dofs_per_cell_to);
305 cell_from->get_dof_indices(local_dof_indices_from);
306 cell_to->get_dof_indices(local_dof_indices_to);
307 }
308
309 template <int dim>
310 CopyDataFace_J &new_face_data(const FEInterfaceValues<dim> &fe_iv_from, const FEInterfaceValues<dim> &fe_iv_to)
311 {
312 auto &copy_data_face = face_data.emplace_back();
313 copy_data_face.cell_jacobian.reinit(fe_iv_to.n_current_interface_dofs(),
314 fe_iv_from.n_current_interface_dofs());
315 copy_data_face.joint_dof_indices_from = fe_iv_from.get_interface_dof_indices();
316 copy_data_face.joint_dof_indices_to = fe_iv_to.get_interface_dof_indices();
317 return copy_data_face;
318 }
319 };
320
321 template <typename NumberType, uint n_fe_subsystems> struct CopyData_J_full {
322 struct CopyDataFace_J {
323 array<FullMatrix<NumberType>, n_fe_subsystems> cell_jacobian;
324 FullMatrix<NumberType> extractor_cell_jacobian;
325 array<vector<types::global_dof_index>, n_fe_subsystems> joint_dof_indices;
326 };
327
328 array<FullMatrix<NumberType>, n_fe_subsystems> cell_jacobian;
329 FullMatrix<NumberType> cell_mass_jacobian;
330 FullMatrix<NumberType> extractor_cell_jacobian;
331 array<vector<types::global_dof_index>, n_fe_subsystems> local_dof_indices;
332 vector<CopyDataFace_J> face_data;
333
334 uint dofs_per_cell;
335
336 template <class Iterator> void reinit(const array<Iterator, n_fe_subsystems> &cell, const uint n_extractors)
337 {
338 const uint n_dofs = cell[0]->get_fe().n_dofs_per_cell();
339 dofs_per_cell = n_dofs;
340 for (uint i = 0; i < n_fe_subsystems; ++i) {
341 const uint from_n_dofs = cell[i]->get_fe().n_dofs_per_cell();
342 if (i == 0) cell_mass_jacobian.reinit(n_dofs, from_n_dofs);
343 cell_jacobian[i].reinit(n_dofs, from_n_dofs);
344 local_dof_indices[i].resize(from_n_dofs);
345 cell[i]->get_dof_indices(local_dof_indices[i]);
346 }
347 if (n_extractors > 0) extractor_cell_jacobian.reinit(dofs_per_cell, n_extractors);
348 }
349
350 template <int dim>
351 CopyDataFace_J &new_face_data(const array<unique_ptr<FEInterfaceValues<dim>>, n_fe_subsystems> &fe_iv,
352 const uint n_extractors)
353 {
354 auto &copy_data_face = face_data.emplace_back();
355 for (uint i = 0; i < n_fe_subsystems; ++i) {
356 copy_data_face.cell_jacobian[i].reinit(fe_iv[0]->n_current_interface_dofs(),
357 fe_iv[i]->n_current_interface_dofs());
358 copy_data_face.joint_dof_indices[i] = fe_iv[i]->get_interface_dof_indices();
359 }
360 if (n_extractors > 0)
361 copy_data_face.extractor_cell_jacobian.reinit(fe_iv[0]->n_current_interface_dofs(), n_extractors);
362 return copy_data_face;
363 }
364 };
365
366 template <typename NumberType> struct CopyData_I {
367 struct CopyFaceData_I {
368 std::array<uint, 2> cell_indices;
369 std::array<double, 2> values;
370 };
371 std::vector<CopyFaceData_I> face_data;
372 double value;
373 uint cell_index;
374 };
375 } // namespace internal
376
383 template <typename Discretization_, typename Model_>
384 class Assembler : public LDGAssemblerBase<Discretization_, Model_>
385 {
386 using Base = LDGAssemblerBase<Discretization_, Model_>;
387
388 public:
389 using Discretization = Discretization_;
390 using Model = Model_;
391 using NumberType = typename Discretization::NumberType;
392 using VectorType = typename Discretization::VectorType;
393
394 using Components = typename Discretization::Components;
395 static constexpr uint dim = Discretization::dim;
396 static constexpr uint stencil = Components::count_fe_subsystems();
397
398 private:
399 template <typename... T> auto fe_conv(std::tuple<T &...> &t) const
400 {
401 if constexpr (stencil == 2)
402 return named_tuple<std::tuple<T &...>, "fe_functions", "LDG1", "extractors", "variables">(t);
403 else if constexpr (stencil == 3)
404 return named_tuple<std::tuple<T &...>, "fe_functions", "LDG1", "LDG2", "extractors", "variables">(t);
405 else if constexpr (stencil == 4)
406 return named_tuple<std::tuple<T &...>, "fe_functions", "LDG1", "LDG2", "LDG3", "extractors", "variables">(t);
407 else
408 throw std::runtime_error("Only <= 3 LDG subsystems are supported.");
409 }
410
411 template <typename... T> auto fe_more_conv(std::tuple<T &...> &t) const
412 {
413 if constexpr (stencil == 2)
414 return named_tuple<std::tuple<T &...>, "fe_functions", "LDG1", "fe_derivatives", "fe_hessians", "extractors",
415 "variables">(t);
416 else if constexpr (stencil == 3)
417 return named_tuple<std::tuple<T &...>, "fe_functions", "LDG1", "LDG2", "fe_derivatives", "fe_hessians",
418 "extractors", "variables">(t);
419 else if constexpr (stencil == 4)
420 return named_tuple<std::tuple<T &...>, "fe_functions", "LDG1", "LDG2", "LDG3", "fe_derivatives",
421 "fe_hessians", "extractors", "variables">(t);
422 else
423 throw std::runtime_error("Only <= 3 LDG subsystems are supported.");
424 }
425
426 template <typename... T> auto ref_conv(std::tuple<T &...> &t) const
427 {
428 if constexpr (stencil == 2)
429 return named_tuple<std::tuple<T &...>, "fe_functions", "LDG1">(t);
430 else if constexpr (stencil == 3)
431 return named_tuple<std::tuple<T &...>, "fe_functions", "LDG1", "LDG2">(t);
432 else if constexpr (stencil == 4)
433 return named_tuple<std::tuple<T &...>, "fe_functions", "LDG1", "LDG2", "LDG3">(t);
434 else
435 throw std::runtime_error("Only <= 3 LDG subsystems are supported.");
436 }
437
438 public:
439 Assembler(Discretization &discretization, Model &model, const JSONValue &json)
440 : Base(discretization, model, json),
441 quadrature(fe.degree + 1 + json.get_uint("/discretization/overintegration")),
442 quadrature_face(fe.degree + 1 + json.get_uint("/discretization/overintegration")),
443 dof_handler_list(discretization.get_dof_handler_list())
444 {
445 static_assert(Components::count_fe_subsystems() > 1, "LDG must have a submodel with index 1.");
446 reinit();
447 }
448
449 virtual void reinit_vector(VectorType &vec) const override { vec.reinit(dof_handler.n_dofs()); }
450
457 virtual void attach_data_output(DataOutput<dim, VectorType> &data_out, const VectorType &solution,
458 const VectorType &variables, const VectorType &dt_solution = VectorType(),
459 const VectorType &residual = VectorType()) override
460 {
461 rebuild_ldg_vectors(solution);
462 readouts(data_out, solution, variables);
463
464 const auto fe_function_names = Components::FEFunction_Descriptor::get_names_vector();
465 std::vector<std::string> fe_function_names_residual;
466 for (const auto &name : fe_function_names)
467 fe_function_names_residual.push_back(name + "_residual");
468 std::vector<std::string> fe_function_names_dot;
469 for (const auto &name : fe_function_names)
470 fe_function_names_dot.push_back(name + "_dot");
471
472 auto &fe_out = data_out.fe_output();
473 fe_out.attach(*dof_handler_list[0], solution, fe_function_names);
474 if (dt_solution.size() > 0) fe_out.attach(dof_handler, dt_solution, fe_function_names_dot);
475 if (residual.size() > 0) fe_out.attach(dof_handler, residual, fe_function_names_residual);
476 for (uint k = 1; k < Components::count_fe_subsystems(); ++k) {
477 sol_vector_vec_tmp[k] = sol_vector[k];
478 fe_out.attach(*dof_handler_list[k], sol_vector_vec_tmp[k], "LDG" + std::to_string(k));
479 }
480 }
481
482 virtual void reinit() override
483 {
484 const auto init_mass = [&](uint i) {
485 // build the sparsity of the mass matrix and mass matrix of all ldg levels
486 auto dofs_per_component = DoFTools::count_dofs_per_fe_component(*(dof_handler_list[i]));
487
488 if (i == 0) {
489 auto n_fe = dofs_per_component.size();
490 for (uint j = 1; j < n_fe; ++j)
491 if (dofs_per_component[j] != dofs_per_component[0])
492 throw std::runtime_error("For LDG the FE basis of all systems must be equal!");
493
494 BlockDynamicSparsityPattern dsp(n_fe, n_fe);
495 for (uint i = 0; i < n_fe; ++i)
496 for (uint j = 0; j < n_fe; ++j)
497 dsp.block(i, j).reinit(dofs_per_component[0], dofs_per_component[0]);
498 dsp.collect_sizes();
499 DoFTools::make_sparsity_pattern(*(dof_handler_list[0]), dsp, discretization.get_constraints(0), true);
500 sparsity_pattern_mass.copy_from(dsp);
501
502 // i do not understand why this is needed.
503 component_mass_matrix_inverse.reinit(sparsity_pattern_mass.block(0, 0));
504 mass_matrix.reinit(sparsity_pattern_mass);
505
506 MatrixCreator::create_mass_matrix(*(dof_handler_list[0]), quadrature, mass_matrix,
507 (const Function<dim, NumberType> *const)nullptr,
508 discretization.get_constraints(0));
509 build_inverse(mass_matrix.block(0, 0), component_mass_matrix_inverse);
510 } else {
511 sol_vector[i].reinit(dofs_per_component);
512 sol_vector_tmp[i].reinit(dofs_per_component);
513 ldg_matrix_built[i] = false;
514 jacobian_tmp_built[i] = false;
515 }
516 };
517
518 const auto init_jacobian = [&](uint i) {
519 // build the jacobian and subjacobians
520 if (i == 0) {
521 build_ldg_sparsity(sparsity_pattern_jacobian, *(dof_handler_list[0]), *(dof_handler_list[0]), stencil,
522 true);
523 for (uint k = 1; k < Components::count_fe_subsystems(); ++k)
524 jacobian_tmp[k].reinit(sparsity_pattern_jacobian);
525 } else {
526 build_ldg_sparsity(sparsity_pattern_ug[i], *(dof_handler_list[0]), *(dof_handler_list[i]), 1);
527 j_ug[i].reinit(sparsity_pattern_ug[i]);
528 }
529 };
530
531 auto init_ldg = [&](uint i) {
532 // build the subjacobian sparsity patterns of all matrices that contribute to the jacobian = uu + ug*gu
533 build_ldg_sparsity(sparsity_pattern_gu[i], *(dof_handler_list[i]), *(dof_handler_list[0]), i);
534 j_gu[i].reinit(sparsity_pattern_gu[i]);
535
536 // these are the "in-between" dependencies of the ldg levels
537 build_ldg_sparsity(sparsity_pattern_wg[i], *(dof_handler_list[i]), *(dof_handler_list[i - 1]), 1);
538 j_wg[i].reinit(sparsity_pattern_wg[i]);
539 j_wg_tmp[i].reinit(sparsity_pattern_wg[i]);
540 };
541
542 Timer timer;
543
544 Base::reinit();
545
546 vector<std::thread> threads;
547 for (uint i = 0; i < Components::count_fe_subsystems(); ++i)
548 threads.emplace_back(std::thread(init_mass, i));
549 for (uint i = 0; i < Components::count_fe_subsystems(); ++i)
550 threads.emplace_back(std::thread(init_jacobian, i));
551 for (uint i = 1; i < Components::count_fe_subsystems(); ++i)
552 threads.emplace_back(std::thread(init_ldg, i));
553 for (auto &t : threads)
554 t.join();
555
556 timings_reinit.push_back(timer.wall_time());
557 }
558
559 virtual void rebuild_jacobian_sparsity() override
560 {
561 build_ldg_sparsity(sparsity_pattern_jacobian, *(dof_handler_list[0]), *(dof_handler_list[0]), stencil, true);
562 for (uint k = 1; k < Components::count_fe_subsystems(); ++k)
563 jacobian_tmp[k].reinit(sparsity_pattern_jacobian);
564 }
565
566 virtual void refinement_indicator(Vector<double> &indicator, const VectorType &solution_global) override
567 {
568 using Iterator = typename Triangulation<dim>::active_cell_iterator;
569 using Scratch = internal::ScratchData<Discretization>;
570 using CopyData = internal::CopyData_I<NumberType>;
571
572 const auto cell_worker = [&](const Iterator &t_cell, Scratch &scratch_data, CopyData &copy_data) {
573 const auto &fe_v = scratch_data.new_fe_values(t_cell);
574 copy_data.cell_index = t_cell->active_cell_index();
575 copy_data.value = 0;
576
577 const auto &JxW = fe_v[0]->get_JxW_values();
578 const auto &q_points = fe_v[0]->get_quadrature_points();
579 const auto &q_indices = fe_v[0]->quadrature_point_indices();
580
581 auto &solution = scratch_data.solution;
582 fe_v[0]->get_function_values(solution_global, solution[0]);
583 for (uint i = 1; i < Components::count_fe_subsystems(); ++i)
584 fe_v[i]->get_function_values(sol_vector[i], solution[i]);
585
586 double local_indicator = 0.;
587 for (const auto &q_index : q_indices) {
588 const auto &x_q = q_points[q_index];
589 auto sol_q = local_sol_q(solution, q_index);
590 model.cell_indicator(local_indicator, x_q, ref_conv(sol_q));
591
592 copy_data.value += JxW[q_index] * local_indicator;
593 }
594 };
595 const auto face_worker = [&](const Iterator &t_cell, const uint &f, const uint &sf, const Iterator &t_ncell,
596 const uint &nf, const unsigned int &nsf, Scratch &scratch_data,
597 CopyData &copy_data) {
598 const auto &fe_iv = scratch_data.new_fe_interface_values(t_cell, f, sf, t_ncell, nf, nsf);
599
600 auto &copy_data_face = copy_data.face_data.emplace_back();
601 copy_data_face.cell_indices[0] = t_cell->active_cell_index();
602 copy_data_face.cell_indices[1] = t_ncell->active_cell_index();
603 copy_data_face.values[0] = 0;
604 copy_data_face.values[1] = 0;
605
606 const auto &JxW = fe_iv[0]->get_JxW_values();
607 const auto &q_points = fe_iv[0]->get_quadrature_points();
608 const auto &q_indices = fe_iv[0]->quadrature_point_indices();
609 const std::vector<Tensor<1, dim>> &normals = fe_iv[0]->get_normal_vectors();
610 array<double, 2> local_indicator{};
611
612 auto &solution = scratch_data.solution_interface;
613 fe_iv[0]->get_fe_face_values(0).get_function_values(solution_global, solution[0][0]);
614 fe_iv[0]->get_fe_face_values(1).get_function_values(solution_global, solution[1][0]);
615 for (uint i = 1; i < Components::count_fe_subsystems(); ++i) {
616 fe_iv[i]->get_fe_face_values(0).get_function_values(sol_vector[i], solution[0][i]);
617 fe_iv[i]->get_fe_face_values(1).get_function_values(sol_vector[i], solution[1][i]);
618 }
619
620 for (const auto &q_index : q_indices) {
621 const auto &x_q = q_points[q_index];
622 auto sol_q_s = local_sol_q(solution[0], q_index);
623 auto sol_q_n = local_sol_q(solution[1], q_index);
624 model.face_indicator(local_indicator, normals[q_index], x_q, ref_conv(sol_q_s), ref_conv(sol_q_n));
625
626 copy_data_face.values[0] += JxW[q_index] * local_indicator[0] * (1. + t_cell->at_boundary());
627 copy_data_face.values[1] += JxW[q_index] * local_indicator[1] * (1. + t_ncell->at_boundary());
628 }
629 };
630 const auto copier = [&](const CopyData &c) {
631 for (auto &cdf : c.face_data)
632 for (uint j = 0; j < 2; ++j)
633 indicator[cdf.cell_indices[j]] += cdf.values[j];
634 indicator[c.cell_index] += c.value;
635 };
636
637 const UpdateFlags update_flags = update_values | update_quadrature_points | update_JxW_values;
638 Scratch scratch_data(mapping, dof_handler_list, quadrature, quadrature_face, update_flags);
639 CopyData copy_data;
640 MeshWorker::AssembleFlags assemble_flags =
641 MeshWorker::assemble_own_cells | MeshWorker::assemble_own_interior_faces_once;
642
643 rebuild_ldg_vectors(solution_global);
644 MeshWorker::mesh_loop(dof_handler.begin_active(), dof_handler.end(), cell_worker, copier, scratch_data,
645 copy_data, assemble_flags, nullptr, face_worker, threads, batch_size);
646 }
647
648 virtual const BlockSparsityPattern &get_sparsity_pattern_jacobian() const override
649 {
650 return sparsity_pattern_jacobian;
651 }
652 virtual const BlockSparseMatrix<NumberType> &get_mass_matrix() const override { return mass_matrix; }
653
661 virtual void mass(VectorType &residual, const VectorType &solution_global, const VectorType &solution_global_dot,
662 NumberType weight) override
663 {
664 using Iterator = typename Triangulation<dim>::active_cell_iterator;
665 using Scratch = internal::ScratchData<Discretization>;
666 using CopyData = internal::CopyData_R<NumberType>;
667 const auto &constraints = discretization.get_constraints();
668
669 const auto cell_worker = [&](const Iterator &t_cell, Scratch &scratch_data, CopyData &copy_data) {
670 const auto &fe_v = scratch_data.new_fe_values(t_cell);
671 const uint n_dofs = fe_v[0]->get_fe().n_dofs_per_cell();
672 copy_data.reinit(scratch_data.cell[0], n_dofs);
673
674 const auto &JxW = fe_v[0]->get_JxW_values();
675 const auto &q_points = fe_v[0]->get_quadrature_points();
676 const auto &q_indices = fe_v[0]->quadrature_point_indices();
677
678 auto &solution = scratch_data.solution;
679 auto &solution_dot = scratch_data.solution_dot;
680 fe_v[0]->get_function_values(solution_global, solution[0]);
681 fe_v[0]->get_function_values(solution_global_dot, solution_dot);
682
683 array<NumberType, Components::count_fe_functions(0)> mass{};
684 for (const auto &q_index : q_indices) {
685 const auto &x_q = q_points[q_index];
686 model.mass(mass, x_q, solution[0][q_index], solution_dot[q_index]);
687
688 for (uint i = 0; i < n_dofs; ++i) {
689 const auto component_i = fe_v[0]->get_fe().system_to_component_index(i).first;
690 copy_data.cell_residual(i) += weight * JxW[q_index] * // dx
691 fe_v[0]->shape_value_component(i, q_index, component_i) *
692 mass[component_i]; // phi_i(x_q) * mass(x_q, u_q)
693 }
694 }
695 };
696 const auto copier = [&](const CopyData &c) {
697 constraints.distribute_local_to_global(c.cell_residual, c.local_dof_indices, residual);
698 };
699
700 const UpdateFlags update_flags = update_values | update_quadrature_points | update_JxW_values;
701 const MeshWorker::AssembleFlags assemble_flags = MeshWorker::assemble_own_cells;
702 Scratch scratch_data(mapping, dof_handler_list, quadrature, quadrature_face, update_flags);
703 CopyData copy_data;
704
705 MeshWorker::mesh_loop(dof_handler.begin_active(), dof_handler.end(), cell_worker, copier, scratch_data,
706 copy_data, assemble_flags, nullptr, nullptr, threads, batch_size);
707 }
708
716 virtual void residual(VectorType &residual, const VectorType &solution_global, NumberType weight,
717 const VectorType &solution_global_dot, NumberType weight_mass,
718 const VectorType &variables = VectorType()) override
719 {
720 using Iterator = typename Triangulation<dim>::active_cell_iterator;
721 using Scratch = internal::ScratchData<Discretization>;
722 using CopyData = internal::CopyData_R<NumberType>;
723 const auto &constraints = discretization.get_constraints();
724
725 // Find the EoM and extract whatever data is needed for the model.
726 std::array<NumberType, Components::count_extractors()> __extracted_data{{}};
727 if constexpr (Components::count_extractors() > 0)
728 this->extract(__extracted_data, solution_global, variables, true, false, true);
729 const auto &extracted_data = __extracted_data;
730
731 const auto cell_worker = [&](const Iterator &t_cell, Scratch &scratch_data, CopyData &copy_data) {
732 const auto &fe_v = scratch_data.new_fe_values(t_cell);
733 const uint n_dofs = fe_v[0]->get_fe().n_dofs_per_cell();
734 copy_data.reinit(scratch_data.cell[0], n_dofs);
735
736 const auto &JxW = fe_v[0]->get_JxW_values();
737 const auto &q_points = fe_v[0]->get_quadrature_points();
738 const auto &q_indices = fe_v[0]->quadrature_point_indices();
739
740 auto &solution = scratch_data.solution;
741 auto &solution_dot = scratch_data.solution_dot;
742 fe_v[0]->get_function_values(solution_global, solution[0]);
743 fe_v[0]->get_function_values(solution_global_dot, solution_dot);
744 for (uint i = 1; i < Components::count_fe_subsystems(); ++i)
745 fe_v[i]->get_function_values(sol_vector[i], solution[i]);
746
747 array<NumberType, Components::count_fe_functions(0)> mass{};
748 array<Tensor<1, dim, NumberType>, Components::count_fe_functions(0)> flux{};
749 array<NumberType, Components::count_fe_functions(0)> source{};
750 for (const auto &q_index : q_indices) {
751 const auto &x_q = q_points[q_index];
752 auto sol_q = std::tuple_cat(local_sol_q(solution, q_index), std::tie(extracted_data, variables));
753 model.mass(mass, x_q, solution[0][q_index], solution_dot[q_index]);
754 model.flux(flux, x_q, fe_conv(sol_q));
755 model.source(source, x_q, fe_conv(sol_q));
756
757 for (uint i = 0; i < n_dofs; ++i) {
758 const auto component_i = fe_v[0]->get_fe().system_to_component_index(i).first;
759 copy_data.cell_residual(i) += weight * JxW[q_index] * // dx
760 (-scalar_product(fe_v[0]->shape_grad_component(i, q_index, component_i),
761 flux[component_i]) // -dphi_i(x_q) * flux(x_q, u_q)
762 + fe_v[0]->shape_value_component(i, q_index, component_i) *
763 (source[component_i])); // -phi_i(x_q) * source(x_q, u_q)
764 copy_data.cell_mass(i) += weight_mass * JxW[q_index] * // dx
765 fe_v[0]->shape_value_component(i, q_index, component_i) *
766 mass[component_i]; // phi_i(x_q) * mass(x_q, u_q)
767 }
768 }
769 };
770 const auto boundary_worker = [&](const Iterator &t_cell, const uint &face_no, Scratch &scratch_data,
771 CopyData &copy_data) {
772 const auto &fe_fv = scratch_data.new_fe_boundary_values(t_cell, face_no);
773 const uint n_dofs = fe_fv[0]->get_fe().n_dofs_per_cell();
774
775 const auto &JxW = fe_fv[0]->get_JxW_values();
776 const auto &q_points = fe_fv[0]->get_quadrature_points();
777 const auto &q_indices = fe_fv[0]->quadrature_point_indices();
778 const auto &normals = fe_fv[0]->get_normal_vectors();
779
780 auto &solution = scratch_data.solution;
781 fe_fv[0]->get_function_values(solution_global, solution[0]);
782 for (uint i = 1; i < Components::count_fe_subsystems(); ++i)
783 fe_fv[i]->get_function_values(sol_vector[i], solution[i]);
784
785 array<Tensor<1, dim, NumberType>, Components::count_fe_functions(0)> numflux{};
786 for (const auto &q_index : q_indices) {
787 const auto &x_q = q_points[q_index];
788 auto sol_q = std::tuple_cat(local_sol_q(solution, q_index), std::tie(extracted_data, variables));
789 model.boundary_numflux(numflux, normals[q_index], x_q, fe_conv(sol_q));
790
791 for (uint i = 0; i < n_dofs; ++i) {
792 const auto component_i = fe_fv[0]->get_fe().system_to_component_index(i).first;
793 copy_data.cell_residual(i) +=
794 weight * JxW[q_index] * // weight * dx
795 (fe_fv[0]->shape_value_component(i, q_index, component_i) *
796 scalar_product(numflux[component_i], normals[q_index])); // phi_i(x_q) * numflux(x_q, u_q) * n(x_q)
797 }
798 }
799 };
800 const auto face_worker = [&](const Iterator &t_cell, const uint &f, const uint &sf, const Iterator &t_ncell,
801 const uint &nf, const unsigned int &nsf, Scratch &scratch_data,
802 CopyData &copy_data) {
803 const auto &fe_iv = scratch_data.new_fe_interface_values(t_cell, f, sf, t_ncell, nf, nsf);
804 const uint n_dofs = fe_iv[0]->n_current_interface_dofs();
805 auto &copy_data_face = copy_data.new_face_data(*(fe_iv[0]));
806
807 const auto &JxW = fe_iv[0]->get_JxW_values();
808 const auto &q_points = fe_iv[0]->get_quadrature_points();
809 const auto &q_indices = fe_iv[0]->quadrature_point_indices();
810 const auto &normals = fe_iv[0]->get_normal_vectors();
811
812 auto &solution = scratch_data.solution_interface;
813 fe_iv[0]->get_fe_face_values(0).get_function_values(solution_global, solution[0][0]);
814 fe_iv[0]->get_fe_face_values(1).get_function_values(solution_global, solution[1][0]);
815 for (uint i = 1; i < Components::count_fe_subsystems(); ++i) {
816 fe_iv[i]->get_fe_face_values(0).get_function_values(sol_vector[i], solution[0][i]);
817 fe_iv[i]->get_fe_face_values(1).get_function_values(sol_vector[i], solution[1][i]);
818 }
819
820 array<Tensor<1, dim, NumberType>, Components::count_fe_functions(0)> numflux{};
821 for (const auto &q_index : q_indices) {
822 const auto &x_q = q_points[q_index];
823 auto sol_q_s = std::tuple_cat(local_sol_q(solution[0], q_index), std::tie(extracted_data, variables));
824 auto sol_q_n = std::tuple_cat(local_sol_q(solution[1], q_index), std::tie(extracted_data, variables));
825 model.numflux(numflux, normals[q_index], x_q, fe_conv(sol_q_s), fe_conv(sol_q_n));
826
827 for (uint i = 0; i < n_dofs; ++i) {
828 const auto &cd_i = fe_iv[0]->interface_dof_to_dof_indices(i);
829 const auto component_i = cd_i[0] == numbers::invalid_unsigned_int
830 ? fe_iv[0]->get_fe().system_to_component_index(cd_i[1]).first
831 : fe_iv[0]->get_fe().system_to_component_index(cd_i[0]).first;
832 copy_data_face.cell_residual(i) +=
833 weight * JxW[q_index] * // weight * dx
834 (fe_iv[0]->jump_in_shape_values(i, q_index, component_i) *
835 scalar_product(numflux[component_i],
836 normals[q_index])); // [[phi_i(x_q)]] * numflux(x_q, u_q) * n(x_q)
837 }
838 }
839 };
840 const auto copier = [&](const CopyData &c) {
841 constraints.distribute_local_to_global(c.cell_residual, c.local_dof_indices, residual);
842 constraints.distribute_local_to_global(c.cell_mass, c.local_dof_indices, residual);
843 for (auto &cdf : c.face_data)
844 constraints.distribute_local_to_global(cdf.cell_residual, cdf.joint_dof_indices, residual);
845 };
846
847 const UpdateFlags update_flags =
848 update_values | update_gradients | update_quadrature_points | update_JxW_values;
849 const MeshWorker::AssembleFlags assemble_flags = MeshWorker::assemble_own_cells |
850 MeshWorker::assemble_boundary_faces |
851 MeshWorker::assemble_own_interior_faces_once;
852 Scratch scratch_data(mapping, dof_handler_list, quadrature, quadrature_face, update_flags);
853 CopyData copy_data;
854
855 Timer timer;
856
857 rebuild_ldg_vectors(solution_global);
858 MeshWorker::mesh_loop(dof_handler.begin_active(), dof_handler.end(), cell_worker, copier, scratch_data,
859 copy_data, assemble_flags, boundary_worker, face_worker, threads, batch_size);
860
861 timings_residual.push_back(timer.wall_time());
862 }
863
864 virtual void jacobian_mass(BlockSparseMatrix<NumberType> &jacobian, const VectorType &solution_global,
865 const VectorType &solution_global_dot, NumberType alpha = 1.,
866 NumberType beta = 1.) override
867 {
868 using Iterator = typename Triangulation<dim>::active_cell_iterator;
869 using Scratch = internal::ScratchData<Discretization>;
870 using CopyData = internal::CopyData_J_full<NumberType, Components::count_fe_subsystems()>;
871 const auto &constraints = discretization.get_constraints();
872
873 const auto cell_worker = [&](const Iterator &t_cell, Scratch &scratch_data, CopyData &copy_data) {
874 const auto &fe_v = scratch_data.new_fe_values(t_cell);
875 const uint to_n_dofs = fe_v[0]->get_fe().n_dofs_per_cell();
876 copy_data.reinit(scratch_data.cell, Components::count_extractors());
877
878 const auto &JxW = fe_v[0]->get_JxW_values();
879 const auto &q_points = fe_v[0]->get_quadrature_points();
880 const auto &q_indices = fe_v[0]->quadrature_point_indices();
881
882 auto &solution = scratch_data.solution;
883 auto &solution_dot = scratch_data.solution_dot;
884
885 fe_v[0]->get_function_values(solution_global, solution[0]);
886 fe_v[0]->get_function_values(solution_global_dot, solution_dot);
887
888 SimpleMatrix<NumberType, Components::count_fe_functions(0)> j_mass;
889 SimpleMatrix<NumberType, Components::count_fe_functions(0)> j_mass_dot;
890
891 const uint from_n_dofs = fe_v[0]->get_fe().n_dofs_per_cell();
892 for (const auto &q_index : q_indices) {
893 const auto &x_q = q_points[q_index];
894
895 model.template jacobian_mass<0>(j_mass, x_q, solution[0][q_index], solution_dot[q_index]);
896 model.template jacobian_mass<1>(j_mass_dot, x_q, solution[0][q_index], solution_dot[q_index]);
897
898 for (uint i = 0; i < to_n_dofs; ++i) {
899 const auto component_i = fe_v[0]->get_fe().system_to_component_index(i).first;
900 for (uint j = 0; j < from_n_dofs; ++j) {
901 const auto component_j = fe_v[0]->get_fe().system_to_component_index(j).first;
902 copy_data.cell_jacobian[0](i, j) +=
903 JxW[q_index] * fe_v[0]->shape_value_component(j, q_index, component_j) * // weight * dx * phi_j(x_q)
904 fe_v[0]->shape_value_component(i, q_index, component_i) *
905 (alpha * j_mass_dot(component_i, component_j) +
906 beta * j_mass(component_i, component_j)); // -phi_i(x_q) * jsource(x_q, u_q)
907 }
908 }
909 }
910 };
911 const auto copier = [&](const CopyData &c) {
912 constraints.distribute_local_to_global(c.cell_jacobian[0], c.local_dof_indices[0], jacobian);
913 };
914
915 const UpdateFlags update_flags = update_values | update_quadrature_points | update_JxW_values;
916 const MeshWorker::AssembleFlags assemble_flags = MeshWorker::assemble_own_cells;
917 Scratch scratch_data(mapping, dof_handler_list, quadrature, quadrature_face, update_flags);
918 CopyData copy_data;
919
920 MeshWorker::mesh_loop(dof_handler.begin_active(), dof_handler.end(), cell_worker, copier, scratch_data,
921 copy_data, assemble_flags, nullptr, nullptr, threads, batch_size);
922 }
923
931 virtual void jacobian(BlockSparseMatrix<NumberType> &jacobian, const VectorType &solution_global,
932 NumberType weight, const VectorType &solution_global_dot, NumberType alpha, NumberType beta,
933 const VectorType &variables = VectorType()) override
934 {
935 if (is_close(weight, 0.))
936 throw std::runtime_error("Please call jacobian_mass instead of jacobian for weight == 0).");
937 using Iterator = typename Triangulation<dim>::active_cell_iterator;
938 using Scratch = internal::ScratchData<Discretization>;
939 using CopyData = internal::CopyData_J_full<NumberType, Components::count_fe_subsystems()>;
940 const auto &constraints = discretization.get_constraints();
941
942 // Find the EoM and extract whatever data is needed for the model.
943 std::array<NumberType, Components::count_extractors()> __extracted_data{{}};
944 if constexpr (Components::count_extractors() > 0) {
945 this->extract(__extracted_data, solution_global, variables, true, true, true);
946 if (this->jacobian_extractors(this->extractor_jacobian, solution_global, variables))
947 jacobian.reinit(sparsity_pattern_jacobian);
948 }
949 const auto &extracted_data = __extracted_data;
950
951 bool exception = false;
952
953 const auto cell_worker = [&](const Iterator &t_cell, Scratch &scratch_data, CopyData &copy_data) {
954 const auto &fe_v = scratch_data.new_fe_values(t_cell);
955 const uint to_n_dofs = fe_v[0]->get_fe().n_dofs_per_cell();
956 copy_data.reinit(scratch_data.cell, Components::count_extractors());
957
958 const auto &JxW = fe_v[0]->get_JxW_values();
959 const auto &q_points = fe_v[0]->get_quadrature_points();
960 const auto &q_indices = fe_v[0]->quadrature_point_indices();
961
962 auto &solution = scratch_data.solution;
963 auto &solution_dot = scratch_data.solution_dot;
964
965 fe_v[0]->get_function_values(solution_global, solution[0]);
966 fe_v[0]->get_function_values(solution_global_dot, solution_dot);
967 for (uint i = 1; i < Components::count_fe_subsystems(); ++i)
968 fe_v[i]->get_function_values(sol_vector[i], solution[i]);
969
970 SimpleMatrix<NumberType, Components::count_fe_functions(0)> j_mass;
971 SimpleMatrix<NumberType, Components::count_fe_functions(0)> j_mass_dot;
972 auto j_flux = jacobian_tuple<Tensor<1, dim, NumberType>, Model>();
973 auto j_source = jacobian_tuple<NumberType, Model>();
974 SimpleMatrix<Tensor<1, dim, NumberType>, Components::count_fe_functions(), Components::count_extractors()>
975 j_extr_flux;
976 SimpleMatrix<NumberType, Components::count_fe_functions(), Components::count_extractors()> j_extr_source;
977 constexpr_for<0, Components::count_fe_subsystems(), 1>([&](auto k) {
978 if (jacobian_tmp_built[k] && model.get_components().jacobians_constant(0, k)) return;
979
980 const uint from_n_dofs = fe_v[k]->get_fe().n_dofs_per_cell();
981 for (const auto &q_index : q_indices) {
982 const auto &x_q = q_points[q_index];
983 auto sol_q = std::tuple_cat(local_sol_q(solution, q_index), std::tie(extracted_data, variables));
984
985 if constexpr (k == 0) {
986 this->model.template jacobian_mass<0>(j_mass, x_q, solution[0][q_index], solution_dot[q_index]);
987 this->model.template jacobian_mass<1>(j_mass_dot, x_q, solution[0][q_index], solution_dot[q_index]);
988 if constexpr (Components::count_extractors() > 0) {
989 this->model.template jacobian_flux_extr<stencil>(j_extr_flux, x_q, fe_conv(sol_q));
990 this->model.template jacobian_source_extr<stencil>(j_extr_source, x_q, fe_conv(sol_q));
991 }
992 }
993 model.template jacobian_flux<k, 0>(std::get<k>(j_flux), x_q, fe_conv(sol_q));
994 model.template jacobian_source<k, 0>(std::get<k>(j_source), x_q, fe_conv(sol_q));
995
996 if (!std::get<k>(j_flux).is_finite() || !std::get<k>(j_source).is_finite()) exception = true;
997
998 for (uint i = 0; i < to_n_dofs; ++i) {
999 const auto component_i = fe_v[0]->get_fe().system_to_component_index(i).first;
1000 for (uint j = 0; j < from_n_dofs; ++j) {
1001 const auto component_j = fe_v[k]->get_fe().system_to_component_index(j).first;
1002
1003 copy_data.cell_jacobian[k](i, j) +=
1004 JxW[q_index] *
1005 fe_v[k]->shape_value_component(j, q_index, component_j) * // weight * dx * phi_j(x_q)
1006 (-scalar_product(fe_v[0]->shape_grad_component(i, q_index, component_i),
1007 std::get<k>(j_flux)(component_i, component_j)) // -dphi_i(x_q) * jflux(x_q, u_q)
1008 + fe_v[0]->shape_value_component(i, q_index, component_i) *
1009 std::get<k>(j_source)(component_i, component_j)); // -phi_i(x_q) * jsource(x_q, u_q)
1010 if constexpr (k == 0) {
1011 copy_data.cell_mass_jacobian(i, j) +=
1012 JxW[q_index] *
1013 fe_v[0]->shape_value_component(j, q_index, component_j) * // weight * dx * phi_j(x_q)
1014 fe_v[0]->shape_value_component(i, q_index, component_i) *
1015 (alpha / weight * j_mass_dot(component_i, component_j) +
1016 beta / weight * j_mass(component_i, component_j)); // -phi_i(x_q) * jsource(x_q, u_q)
1017 }
1018 }
1019
1020 // extractor contribution
1021 if constexpr (k == 0)
1022 if constexpr (Components::count_extractors() > 0)
1023 for (uint e = 0; e < Components::count_extractors(); ++e)
1024 copy_data.extractor_cell_jacobian(i, e) +=
1025 weight * JxW[q_index] * // dx * phi_j * (
1026 (-scalar_product(fe_v[0]->shape_grad_component(i, q_index, component_i),
1027 j_extr_flux(component_i, e)) // -dphi_i * jflux
1028 + fe_v[0]->shape_value_component(i, q_index, component_i) *
1029 j_extr_source(component_i, e)); // -phi_i * jsource)
1030 }
1031 }
1032 });
1033 };
1034 const auto boundary_worker = [&](const Iterator &t_cell, const uint &face_no, Scratch &scratch_data,
1035 CopyData &copy_data) {
1036 const auto &fe_fv = scratch_data.new_fe_boundary_values(t_cell, face_no);
1037 const uint to_n_dofs = fe_fv[0]->get_fe().n_dofs_per_cell();
1038
1039 const auto &JxW = fe_fv[0]->get_JxW_values();
1040 const auto &q_points = fe_fv[0]->get_quadrature_points();
1041 const auto &q_indices = fe_fv[0]->quadrature_point_indices();
1042 const auto &normals = fe_fv[0]->get_normal_vectors();
1043 auto &solution = scratch_data.solution;
1044
1045 fe_fv[0]->get_function_values(solution_global, solution[0]);
1046 for (uint i = 1; i < Components::count_fe_subsystems(); ++i)
1047 fe_fv[i]->get_function_values(sol_vector[i], solution[i]);
1048
1049 auto j_boundary_numflux = jacobian_tuple<Tensor<1, dim>, Model>();
1050 constexpr_for<0, Components::count_fe_subsystems(), 1>([&](auto k) {
1051 if (jacobian_tmp_built[k] && model.get_components().jacobians_constant(0, k)) return;
1052 for (const auto &q_index : q_indices) {
1053 const auto &x_q = q_points[q_index];
1054 auto sol_q = std::tuple_cat(local_sol_q(solution, q_index), std::tie(extracted_data, variables));
1055
1056 const uint from_n_dofs = fe_fv[k]->get_fe().n_dofs_per_cell();
1057
1058 model.template jacobian_boundary_numflux<k, 0>(std::get<k>(j_boundary_numflux), normals[q_index], x_q,
1059 fe_conv(sol_q));
1060
1061 if (!std::get<k>(j_boundary_numflux).is_finite()) exception = true;
1062
1063 for (uint i = 0; i < to_n_dofs; ++i) {
1064 const auto component_i = fe_fv[0]->get_fe().system_to_component_index(i).first;
1065 for (uint j = 0; j < from_n_dofs; ++j) {
1066 const auto component_j = fe_fv[k]->get_fe().system_to_component_index(j).first;
1067
1068 copy_data.cell_jacobian[k](i, j) +=
1069 JxW[q_index] *
1070 fe_fv[k]->shape_value_component(j, q_index, component_j) * // weight * dx * phi_j(x_q)
1071 (fe_fv[0]->shape_value_component(i, q_index, component_i) *
1072 scalar_product(std::get<k>(j_boundary_numflux)(component_i, component_j),
1073 normals[q_index])); // phi_i(x_q) * j_numflux(x_q, u_q) * n(x_q)
1074 }
1075 }
1076 }
1077 });
1078 };
1079 const auto face_worker = [&](const Iterator &t_cell, const uint &f, const uint &sf, const Iterator &t_ncell,
1080 const uint &nf, const unsigned int &nsf, Scratch &scratch_data,
1081 CopyData &copy_data) {
1082 const auto &fe_iv = scratch_data.new_fe_interface_values(t_cell, f, sf, t_ncell, nf, nsf);
1083 const uint to_n_dofs = fe_iv[0]->n_current_interface_dofs();
1084 auto &copy_data_face = copy_data.new_face_data(fe_iv, Components::count_extractors());
1085
1086 const auto &JxW = fe_iv[0]->get_JxW_values();
1087 const auto &q_points = fe_iv[0]->get_quadrature_points();
1088 const auto &q_indices = fe_iv[0]->quadrature_point_indices();
1089 const auto &normals = fe_iv[0]->get_normal_vectors();
1090
1091 auto &solution = scratch_data.solution_interface;
1092 fe_iv[0]->get_fe_face_values(0).get_function_values(solution_global, solution[0][0]);
1093 fe_iv[0]->get_fe_face_values(1).get_function_values(solution_global, solution[1][0]);
1094 for (uint i = 1; i < Components::count_fe_subsystems(); ++i) {
1095 fe_iv[i]->get_fe_face_values(0).get_function_values(sol_vector[i], solution[0][i]);
1096 fe_iv[i]->get_fe_face_values(1).get_function_values(sol_vector[i], solution[1][i]);
1097 }
1098
1099 auto j_numflux = jacobian_2_tuple<Tensor<1, dim>, Model>();
1100 constexpr_for<0, Components::count_fe_subsystems(), 1>([&](auto k) {
1101 if (jacobian_tmp_built[k] && model.get_components().jacobians_constant(0, k)) return;
1102 for (const auto &q_index : q_indices) {
1103 const auto &x_q = q_points[q_index];
1104 auto sol_q_s = std::tuple_cat(local_sol_q(solution[0], q_index), std::tie(extracted_data, variables));
1105 auto sol_q_n = std::tuple_cat(local_sol_q(solution[1], q_index), std::tie(extracted_data, variables));
1106
1107 const uint from_n_dofs = fe_iv[k]->n_current_interface_dofs();
1108
1109 model.template jacobian_numflux<k, 0>(std::get<k>(j_numflux), normals[q_index], x_q, fe_conv(sol_q_s),
1110 fe_conv(sol_q_n));
1111
1112 if (!std::get<k>(j_numflux)[0].is_finite() || !std::get<k>(j_numflux)[1].is_finite()) exception = true;
1113
1114 for (uint i = 0; i < to_n_dofs; ++i) {
1115 const auto &cd_i = fe_iv[0]->interface_dof_to_dof_indices(i);
1116 const uint face_no_i = cd_i[0] == numbers::invalid_unsigned_int ? 1 : 0;
1117 const auto &component_i = fe_iv[0]->get_fe().system_to_component_index(cd_i[face_no_i]).first;
1118 for (uint j = 0; j < from_n_dofs; ++j) {
1119 const auto &cd_j = fe_iv[k]->interface_dof_to_dof_indices(j);
1120 const uint face_no_j = cd_j[0] == numbers::invalid_unsigned_int ? 1 : 0;
1121 const auto &component_j = fe_iv[k]->get_fe().system_to_component_index(cd_j[face_no_j]).first;
1122
1123 copy_data_face.cell_jacobian[k](i, j) +=
1124 JxW[q_index] *
1125 fe_iv[k]->get_fe_face_values(face_no_j).shape_value_component(
1126 cd_j[face_no_j], q_index, component_j) * // weight * dx * phi_j(x_q)
1127 (fe_iv[0]->jump_in_shape_values(i, q_index, component_i) *
1128 scalar_product(std::get<k>(j_numflux)[face_no_j](component_i, component_j),
1129 normals[q_index])); // [[phi_i(x_q)]] * j_numflux(x_q, u_q)
1130 }
1131 }
1132 }
1133 });
1134 };
1135 const auto copier = [&](const CopyData &c) {
1136 try {
1137 constraints.distribute_local_to_global(c.cell_jacobian[0], c.local_dof_indices[0], jacobian);
1138 constraints.distribute_local_to_global(c.cell_mass_jacobian, c.local_dof_indices[0], jacobian);
1139 for (auto &cdf : c.face_data) {
1140 constraints.distribute_local_to_global(cdf.cell_jacobian[0], cdf.joint_dof_indices[0], jacobian);
1141 if constexpr (Components::count_extractors() > 0) {
1142 FullMatrix<NumberType> extractor_dependence(cdf.joint_dof_indices[0].size(),
1143 extractor_dof_indices.size());
1144 cdf.extractor_cell_jacobian.mmult(extractor_dependence, this->extractor_jacobian);
1145 constraints.distribute_local_to_global(extractor_dependence, cdf.joint_dof_indices[0],
1146 extractor_dof_indices, jacobian);
1147 }
1148 }
1149 if constexpr (Components::count_extractors() > 0) {
1150 FullMatrix<NumberType> extractor_dependence(c.local_dof_indices[0].size(), extractor_dof_indices.size());
1151 c.extractor_cell_jacobian.mmult(extractor_dependence, this->extractor_jacobian);
1152 constraints.distribute_local_to_global(extractor_dependence, c.local_dof_indices[0],
1153 extractor_dof_indices, jacobian);
1154 }
1155
1156 // LDG things
1157
1158 for (uint i = 1; i < Components::count_fe_subsystems(); ++i) {
1159 if (jacobian_tmp_built[i] && model.get_components().jacobians_constant(0, i)) continue;
1160 j_ug[i].add(c.local_dof_indices[0], c.local_dof_indices[i], c.cell_jacobian[i]);
1161 }
1162 for (auto &cdf : c.face_data) {
1163 for (uint i = 1; i < Components::count_fe_subsystems(); ++i) {
1164 if (jacobian_tmp_built[i] && model.get_components().jacobians_constant(0, i)) continue;
1165 j_ug[i].add(cdf.joint_dof_indices[0], cdf.joint_dof_indices[i], cdf.cell_jacobian[i]);
1166 }
1167 }
1168 } catch (...) {
1169 exception = true;
1170 }
1171 };
1172
1173 const UpdateFlags update_flags =
1174 update_values | update_gradients | update_quadrature_points | update_JxW_values;
1175 const MeshWorker::AssembleFlags assemble_flags = MeshWorker::assemble_own_cells |
1176 MeshWorker::assemble_boundary_faces |
1177 MeshWorker::assemble_own_interior_faces_once;
1178 Scratch scratch_data(mapping, dof_handler_list, quadrature, quadrature_face, update_flags);
1179 CopyData copy_data;
1180
1181 Timer timer;
1182
1183 constexpr_for<1, Components::count_fe_subsystems(), 1>([&](auto k) {
1184 if (!ldg_matrix_built[k] || !model.get_components().jacobians_constant(k, k - 1))
1185 rebuild_ldg_jacobian<k>(solution_global);
1186 if (!jacobian_tmp_built[k] || !model.get_components().jacobians_constant(0, k)) j_ug[k] = 0;
1187 });
1188 rebuild_ldg_vectors(solution_global);
1189
1190 MeshWorker::mesh_loop(dof_handler.begin_active(), dof_handler.end(), cell_worker, copier, scratch_data,
1191 copy_data, assemble_flags, boundary_worker, face_worker, threads, batch_size);
1192
1193 if (exception) throw std::runtime_error("Infinity encountered in jacobian construction");
1194
1195 tbb::parallel_for(
1196 tbb::blocked_range<uint>(1, Components::count_fe_subsystems()), [&](tbb::blocked_range<uint> rk) {
1197 for (uint k = rk.begin(); k < rk.end(); ++k) {
1198 if (!jacobian_tmp_built[k] || !model.get_components().jacobians_constant(0, k)) {
1199 jacobian_tmp[k] = 0;
1200 tbb::parallel_for(
1201 tbb::blocked_range<uint>(0, Components::count_fe_functions(0)), [&](tbb::blocked_range<uint> r) {
1202 for (uint q = r.begin(); q < r.end(); ++q)
1203 for (const auto &c : model.get_components().ldg_couplings(k, 0))
1204 j_ug[k].block(q, c[0]).mmult(jacobian_tmp[k].block(q, c[1]), j_gu[k].block(c[0], c[1]),
1205 Vector<NumberType>(), false);
1206 });
1207 jacobian_tmp_built[k] = true;
1208 }
1209 }
1210 });
1211
1212 tbb::parallel_for(
1213 tbb::blocked_range<uint>(0, Components::count_fe_functions(0)), [&](tbb::blocked_range<uint> rc1) {
1214 tbb::parallel_for(tbb::blocked_range<uint>(0, Components::count_fe_functions(0)),
1215 [&](tbb::blocked_range<uint> rc2) {
1216 for (uint c1 = rc1.begin(); c1 < rc1.end(); ++c1)
1217 for (uint c2 = rc2.begin(); c2 < rc2.end(); ++c2) {
1218 for (uint k = 1; k < Components::count_fe_subsystems(); ++k)
1219 jacobian.block(c1, c2).add(NumberType(1.), jacobian_tmp[k].block(c1, c2));
1220 jacobian.block(c1, c2) *= weight;
1221 }
1222 });
1223 });
1224
1225 timings_jacobian.push_back(timer.wall_time());
1226 }
1227
1228 void log(const std::string logger)
1229 {
1230 std::stringstream ss;
1231 ss << "LDG Assembler: " << std::endl;
1232 ss << " Reinit: " << average_time_reinit() * 1000 << "ms (" << num_reinits() << ")" << std::endl;
1233 ss << " Residual: " << average_time_residual_assembly() * 1000 << "ms (" << num_residuals() << ")"
1234 << std::endl;
1235 ss << " Jacobian: " << average_time_jacobian_assembly() * 1000 << "ms (" << num_jacobians() << ")"
1236 << std::endl;
1237 spdlog::get(logger)->info(ss.str());
1238 }
1239
1240 double average_time_reinit() const
1241 {
1242 double t = 0.;
1243 double n = timings_reinit.size();
1244 for (const auto &t_ : timings_reinit)
1245 t += t_ / n;
1246 return t;
1247 }
1248 uint num_reinits() const { return timings_reinit.size(); }
1249
1250 double average_time_residual_assembly() const
1251 {
1252 double t = 0.;
1253 double n = timings_residual.size();
1254 for (const auto &t_ : timings_residual)
1255 t += t_ / n;
1256 return t;
1257 }
1258 uint num_residuals() const { return timings_residual.size(); }
1259
1260 double average_time_jacobian_assembly() const
1261 {
1262 double t = 0.;
1263 double n = timings_jacobian.size();
1264 for (const auto &t_ : timings_jacobian)
1265 t += t_ / n;
1266 return t;
1267 }
1268 uint num_jacobians() const { return timings_jacobian.size(); }
1269
1270 protected:
1271 using Base::discretization;
1272 using Base::dof_handler;
1273 using Base::fe;
1274 using Base::mapping;
1275 using Base::model;
1276
1277 QGauss<dim> quadrature;
1278 QGauss<dim - 1> quadrature_face;
1279 using Base::batch_size;
1280 using Base::threads;
1281
1282 std::vector<const DoFHandler<dim> *> dof_handler_list;
1283
1284 mutable array<BlockVector<NumberType>, Components::count_fe_subsystems()> sol_vector;
1285 mutable array<BlockVector<NumberType>, Components::count_fe_subsystems()> sol_vector_tmp;
1286 mutable array<Vector<NumberType>, Components::count_fe_subsystems()> sol_vector_vec_tmp;
1287
1288 BlockSparsityPattern sparsity_pattern_jacobian;
1289 BlockSparsityPattern sparsity_pattern_mass;
1290 array<BlockSparsityPattern, Components::count_fe_subsystems()> sparsity_pattern_ug;
1291 array<BlockSparsityPattern, Components::count_fe_subsystems()> sparsity_pattern_gu;
1292 array<BlockSparsityPattern, Components::count_fe_subsystems()> sparsity_pattern_wg;
1293
1294 array<BlockSparseMatrix<NumberType>, Components::count_fe_subsystems()> jacobian_tmp;
1295
1296 BlockSparseMatrix<NumberType> mass_matrix;
1297 SparseMatrix<NumberType> component_mass_matrix_inverse;
1298 array<BlockSparseMatrix<NumberType>, Components::count_fe_subsystems()> j_ug;
1299 mutable array<BlockSparseMatrix<NumberType>, Components::count_fe_subsystems()> j_gu;
1300 mutable array<BlockSparseMatrix<NumberType>, Components::count_fe_subsystems()> j_wg;
1301 mutable array<BlockSparseMatrix<NumberType>, Components::count_fe_subsystems()> j_wg_tmp;
1302
1303 std::vector<double> timings_reinit;
1304 std::vector<double> timings_residual;
1305 std::vector<double> timings_jacobian;
1306
1307 mutable array<bool, Components::count_fe_subsystems()> ldg_matrix_built;
1308 mutable array<bool, Components::count_fe_subsystems()> jacobian_tmp_built;
1309
1310 using Base::EoM_abs_tol;
1311 using Base::EoM_max_iter;
1312 using Base::extractor_dof_indices;
1313
1314 void rebuild_ldg_vectors(const VectorType &sol) const
1315 {
1316 constexpr_for<1, Components::count_fe_subsystems(), 1>([&](auto k) {
1317 if (!model.get_components().jacobians_constant(k, k - 1)) {
1318 if (k == 1)
1319 build_ldg_vector<k - 1, k>(sol, sol_vector[k], sol_vector_tmp[k]);
1320 else
1321 build_ldg_vector<k - 1, k>(sol_vector[k - 1], sol_vector[k], sol_vector_tmp[k]);
1322 } else {
1323 if (!ldg_matrix_built[k]) rebuild_ldg_jacobian<k>(sol);
1324
1325 j_gu[k].vmult(sol_vector[k], sol);
1326 }
1327 });
1328 }
1329
1330 template <int k> void rebuild_ldg_jacobian(const VectorType &sol) const
1331 {
1332 static_assert(k > 0);
1333 ldg_matrix_built[k] = true;
1334 if constexpr (k == 1)
1335 build_ldg_jacobian<0, 1>(sol, j_gu[1], j_wg_tmp[1]);
1336 else {
1337 if (!ldg_matrix_built[k - 1]) rebuild_ldg_jacobian<k - 1>(sol);
1338 build_ldg_jacobian<k - 1, k>(sol_vector[k - 1], j_wg[k], j_wg_tmp[k]);
1339 for (const auto &c : model.get_components().ldg_couplings(k, 0))
1340 for (const auto &b : model.get_components().ldg_couplings(k, k - 1))
1341 if (b[0] == c[0])
1342 j_wg[k]
1343 .block(b[0], b[1])
1344 .mmult(j_gu[k].block(c[0], c[1]), j_gu[k - 1].block(b[1], c[1]), Vector<NumberType>(), false);
1345 }
1346 }
1347
1357 template <int from, int to, typename VectorType, typename VectorTypeldg>
1358 void build_ldg_vector(const VectorType &solution_global, VectorTypeldg &ldg_vector,
1359 VectorTypeldg &ldg_vector_tmp) const
1360 {
1361 static_assert(to - from == 1, "can only build LDG from last level!");
1362 using Iterator = typename Triangulation<dim>::active_cell_iterator;
1363 using Scratch = internal::ScratchData<Discretization>;
1364 using CopyData = internal::CopyData_R<NumberType>;
1365 const auto &constraints = discretization.get_constraints(to);
1366
1367 const auto cell_worker = [&](const Iterator &t_cell, Scratch &scratch_data, CopyData &copy_data) {
1368 const auto &fe_v = scratch_data.new_fe_values(t_cell);
1369 const uint to_n_dofs = fe_v[to]->get_fe().n_dofs_per_cell();
1370 copy_data.reinit(scratch_data.cell[to], to_n_dofs);
1371
1372 const auto &JxW = fe_v[to]->get_JxW_values();
1373 const auto &q_points = fe_v[to]->get_quadrature_points();
1374 const auto &q_indices = fe_v[to]->quadrature_point_indices();
1375 auto &solution = scratch_data.solution;
1376
1377 fe_v[from]->get_function_values(solution_global, solution[from]);
1378
1379 array<Tensor<1, dim, NumberType>, Components::count_fe_functions(to)> flux{};
1380 array<NumberType, Components::count_fe_functions(to)> source{};
1381 for (const auto &q_index : q_indices) {
1382 const auto &x_q = q_points[q_index];
1383 model.template ldg_flux<to>(flux, x_q, solution[from][q_index]);
1384 model.template ldg_source<to>(source, x_q, solution[from][q_index]);
1385
1386 for (uint i = 0; i < to_n_dofs; ++i) {
1387 const auto component_i = fe_v[to]->get_fe().system_to_component_index(i).first;
1388 copy_data.cell_residual(i) += JxW[q_index] * // dx
1389 (-scalar_product(fe_v[to]->shape_grad_component(i, q_index, component_i),
1390 flux[component_i]) // -dphi_i(x_q) * flux(x_q, u_q)
1391 + fe_v[to]->shape_value_component(i, q_index, component_i) *
1392 source[component_i]); // -phi_i(x_q) * source(x_q, u_q)
1393 }
1394 }
1395 };
1396 const auto boundary_worker = [&](const Iterator &t_cell, const uint &face_no, Scratch &scratch_data,
1397 CopyData &copy_data) {
1398 const auto &fe_fv = scratch_data.new_fe_boundary_values(t_cell, face_no);
1399 const uint to_n_dofs = fe_fv[to]->get_fe().n_dofs_per_cell();
1400
1401 const auto &JxW = fe_fv[to]->get_JxW_values();
1402 const auto &q_points = fe_fv[to]->get_quadrature_points();
1403 const auto &q_indices = fe_fv[to]->quadrature_point_indices();
1404 const std::vector<Tensor<1, dim>> &normals = fe_fv[from]->get_normal_vectors();
1405 auto &solution = scratch_data.solution;
1406
1407 fe_fv[from]->get_function_values(solution_global, solution[from]);
1408
1409 array<Tensor<1, dim, NumberType>, Components::count_fe_functions(to)> numflux{};
1410 for (const auto &q_index : q_indices) {
1411 const auto &x_q = q_points[q_index];
1412 model.template ldg_boundary_numflux<to>(numflux, normals[q_index], x_q, solution[from][q_index]);
1413
1414 for (uint i = 0; i < to_n_dofs; ++i) {
1415 const auto component_i = fe_fv[to]->get_fe().system_to_component_index(i).first;
1416 copy_data.cell_residual(i) +=
1417 JxW[q_index] * // weight * dx
1418 (fe_fv[to]->shape_value_component(i, q_index, component_i) *
1419 scalar_product(numflux[component_i], normals[q_index])); // phi_i(x_q) * numflux(x_q, u_q) * n(x_q)
1420 }
1421 }
1422 };
1423 const auto face_worker = [&](const Iterator &t_cell, const uint &f, const uint &sf, const Iterator &t_ncell,
1424 const uint &nf, const unsigned int &nsf, Scratch &scratch_data,
1425 CopyData &copy_data) {
1426 const auto &fe_iv = scratch_data.new_fe_interface_values(t_cell, f, sf, t_ncell, nf, nsf);
1427 const uint to_n_dofs = fe_iv[to]->n_current_interface_dofs();
1428 auto &copy_data_face = copy_data.new_face_data(*(fe_iv[to]));
1429
1430 const auto &JxW = fe_iv[to]->get_JxW_values();
1431 const auto &q_points = fe_iv[to]->get_quadrature_points();
1432 const auto &q_indices = fe_iv[to]->quadrature_point_indices();
1433 // normals are facing outwards!
1434 const std::vector<Tensor<1, dim>> &normals = fe_iv[to]->get_normal_vectors();
1435 auto &solution = scratch_data.solution_interface;
1436
1437 fe_iv[from]->get_fe_face_values(0).get_function_values(solution_global, solution[0][from]);
1438 fe_iv[from]->get_fe_face_values(1).get_function_values(solution_global, solution[1][from]);
1439
1440 array<Tensor<1, dim, NumberType>, Components::count_fe_functions(to)> numflux{};
1441 for (const auto &q_index : q_indices) {
1442 const auto &x_q = q_points[q_index];
1443 model.template ldg_numflux<to>(numflux, normals[q_index], x_q, solution[0][from][q_index],
1444 solution[1][from][q_index]);
1445
1446 for (uint i = 0; i < to_n_dofs; ++i) {
1447 const auto &cd_i = fe_iv[to]->interface_dof_to_dof_indices(i);
1448 const auto component_i = cd_i[0] == numbers::invalid_unsigned_int
1449 ? fe_iv[to]->get_fe().system_to_component_index(cd_i[1]).first
1450 : fe_iv[to]->get_fe().system_to_component_index(cd_i[0]).first;
1451 copy_data_face.cell_residual(i) +=
1452 JxW[q_index] * // weight * dx
1453 (fe_iv[to]->jump_in_shape_values(i, q_index, component_i) *
1454 scalar_product(numflux[component_i],
1455 normals[q_index])); // [[phi_i(x_q)]] * numflux(x_q, u_q) * n(x_q)
1456 }
1457 }
1458 };
1459 const auto copier = [&](const CopyData &c) {
1460 constraints.distribute_local_to_global(c.cell_residual, c.local_dof_indices, ldg_vector_tmp);
1461 for (auto &cdf : c.face_data)
1462 constraints.distribute_local_to_global(cdf.cell_residual, cdf.joint_dof_indices, ldg_vector_tmp);
1463 };
1464
1465 const MeshWorker::AssembleFlags assemble_flags = MeshWorker::assemble_own_cells |
1466 MeshWorker::assemble_boundary_faces |
1467 MeshWorker::assemble_own_interior_faces_once;
1468 Scratch scratch_data(mapping, dof_handler_list, quadrature, quadrature_face);
1469 CopyData copy_data;
1470
1471 ldg_vector_tmp = 0;
1472 MeshWorker::mesh_loop(dof_handler.begin_active(), dof_handler.end(), cell_worker, copier, scratch_data,
1473 copy_data, assemble_flags, boundary_worker, face_worker, threads, batch_size);
1474
1475 for (uint i = 0; i < Components::count_fe_functions(to); ++i)
1476 component_mass_matrix_inverse.vmult(ldg_vector.block(i), ldg_vector_tmp.block(i));
1477 }
1478
1488 template <int from, int to, typename VectorType>
1489 void build_ldg_jacobian(const VectorType &solution_global, BlockSparseMatrix<NumberType> &ldg_jacobian,
1490 BlockSparseMatrix<NumberType> &ldg_jacobian_tmp) const
1491 {
1492 static_assert(to - from == 1, "can only build LDG from last level!");
1493 using Iterator = typename Triangulation<dim>::active_cell_iterator;
1494 using Scratch = internal::ScratchData<Discretization>;
1495 using CopyData = internal::CopyData_J<NumberType>;
1496
1497 const auto cell_worker = [&](const Iterator &t_cell, Scratch &scratch_data, CopyData &copy_data) {
1498 const auto &fe_v = scratch_data.new_fe_values(t_cell);
1499 const uint to_n_dofs = fe_v[to]->get_fe().n_dofs_per_cell();
1500 const uint from_n_dofs = fe_v[from]->get_fe().n_dofs_per_cell();
1501 copy_data.reinit(scratch_data.cell[from], scratch_data.cell[to], from_n_dofs, to_n_dofs);
1502
1503 const auto &JxW = fe_v[to]->get_JxW_values();
1504 const auto &q_points = fe_v[to]->get_quadrature_points();
1505 const auto &q_indices = fe_v[to]->quadrature_point_indices();
1506 auto &solution = scratch_data.solution;
1507
1508 fe_v[from]->get_function_values(solution_global, solution[from]);
1509
1510 SimpleMatrix<Tensor<1, dim>, Components::count_fe_functions(to), Components::count_fe_functions(from)> j_flux;
1511 SimpleMatrix<NumberType, Components::count_fe_functions(to), Components::count_fe_functions(from)> j_source;
1512 for (const auto &q_index : q_indices) {
1513 const auto &x_q = q_points[q_index];
1514 model.template jacobian_flux<from, to>(j_flux, x_q, solution[from][q_index]);
1515 model.template jacobian_source<from, to>(j_source, x_q, solution[from][q_index]);
1516
1517 for (uint i = 0; i < to_n_dofs; ++i) {
1518 const auto component_i = fe_v[to]->get_fe().system_to_component_index(i).first;
1519 for (uint j = 0; j < from_n_dofs; ++j) {
1520 const auto component_j = fe_v[from]->get_fe().system_to_component_index(j).first;
1521 copy_data.cell_jacobian(i, j) +=
1522 JxW[q_index] *
1523 fe_v[from]->shape_value_component(j, q_index, component_j) * // weight * dx * phi_j(x_q)
1524 (-scalar_product(fe_v[to]->shape_grad_component(i, q_index, component_i),
1525 j_flux(component_i, component_j)) // -dphi_i(x_q) * jflux(x_q, u_q)
1526 + fe_v[to]->shape_value_component(i, q_index, component_i) *
1527 j_source(component_i, component_j)); // -phi_i(x_q) * jsource(x_q, u_q)
1528 }
1529 }
1530 }
1531 };
1532 const auto boundary_worker = [&](const Iterator &t_cell, const uint &face_no, Scratch &scratch_data,
1533 CopyData &copy_data) {
1534 const auto &fe_fv = scratch_data.new_fe_boundary_values(t_cell, face_no);
1535 const uint to_n_dofs = fe_fv[to]->get_fe().n_dofs_per_cell();
1536 const uint from_n_dofs = fe_fv[from]->get_fe().n_dofs_per_cell();
1537
1538 const auto &JxW = fe_fv[to]->get_JxW_values();
1539 const auto &q_points = fe_fv[to]->get_quadrature_points();
1540 const auto &q_indices = fe_fv[to]->quadrature_point_indices();
1541 const std::vector<Tensor<1, dim>> &normals = fe_fv[to]->get_normal_vectors();
1542 auto &solution = scratch_data.solution;
1543
1544 fe_fv[from]->get_function_values(solution_global, solution[from]);
1545
1546 SimpleMatrix<Tensor<1, dim>, Components::count_fe_functions(to), Components::count_fe_functions(from)>
1547 j_boundary_numflux;
1548 for (const auto &q_index : q_indices) {
1549 const auto &x_q = q_points[q_index];
1550 model.template jacobian_boundary_numflux<from, to>(j_boundary_numflux, normals[q_index], x_q,
1551 solution[from][q_index]);
1552
1553 for (uint i = 0; i < to_n_dofs; ++i) {
1554 const auto component_i = fe_fv[to]->get_fe().system_to_component_index(i).first;
1555 for (uint j = 0; j < from_n_dofs; ++j) {
1556 const auto component_j = fe_fv[from]->get_fe().system_to_component_index(j).first;
1557 copy_data.cell_jacobian(i, j) +=
1558 JxW[q_index] *
1559 fe_fv[from]->shape_value_component(j, q_index, component_j) * // weight * dx * phi_j(x_q)
1560 (fe_fv[to]->shape_value_component(i, q_index, component_i) *
1561 scalar_product(j_boundary_numflux(component_i, component_j),
1562 normals[q_index])); // phi_i(x_q) * j_numflux(x_q, u_q) * n(x_q)
1563 }
1564 }
1565 }
1566 };
1567 const auto face_worker = [&](const Iterator &t_cell, const uint &f, const uint &sf, const Iterator &t_ncell,
1568 const uint &nf, const unsigned int &nsf, Scratch &scratch_data,
1569 CopyData &copy_data) {
1570 const auto &fe_iv = scratch_data.new_fe_interface_values(t_cell, f, sf, t_ncell, nf, nsf);
1571 const uint to_n_dofs = fe_iv[to]->n_current_interface_dofs();
1572 const uint from_n_dofs = fe_iv[from]->n_current_interface_dofs();
1573 auto &copy_data_face = copy_data.new_face_data(*(fe_iv[from]), *(fe_iv[to]));
1574
1575 const auto &JxW = fe_iv[to]->get_JxW_values();
1576 const auto &q_points = fe_iv[to]->get_quadrature_points();
1577 const auto &q_indices = fe_iv[to]->quadrature_point_indices();
1578 const std::vector<Tensor<1, dim>> &normals = fe_iv[to]->get_normal_vectors();
1579 auto &solution = scratch_data.solution_interface;
1580
1581 fe_iv[from]->get_fe_face_values(0).get_function_values(solution_global, solution[0][from]);
1582 fe_iv[from]->get_fe_face_values(1).get_function_values(solution_global, solution[1][from]);
1583
1584 array<SimpleMatrix<Tensor<1, dim>, Components::count_fe_functions(to), Components::count_fe_functions(from)>,
1585 2>
1586 j_numflux;
1587 for (const auto &q_index : q_indices) {
1588 const auto &x_q = q_points[q_index];
1589 model.template jacobian_numflux<from, to>(j_numflux, normals[q_index], x_q, solution[0][from][q_index],
1590 solution[1][from][q_index]);
1591
1592 for (uint i = 0; i < to_n_dofs; ++i) {
1593 const auto &cd_i = fe_iv[to]->interface_dof_to_dof_indices(i);
1594 const uint face_no_i = cd_i[0] == numbers::invalid_unsigned_int ? 1 : 0;
1595 const auto &component_i = fe_iv[to]->get_fe().system_to_component_index(cd_i[face_no_i]).first;
1596 for (uint j = 0; j < from_n_dofs; ++j) {
1597 const auto &cd_j = fe_iv[from]->interface_dof_to_dof_indices(j);
1598 const uint face_no_j = cd_j[0] == numbers::invalid_unsigned_int ? 1 : 0;
1599 const auto &component_j = fe_iv[from]->get_fe().system_to_component_index(cd_j[face_no_j]).first;
1600
1601 copy_data_face.cell_jacobian(i, j) +=
1602 JxW[q_index] *
1603 fe_iv[from]->get_fe_face_values(face_no_j).shape_value_component(
1604 cd_j[face_no_j], q_index, component_j) * // weight * dx * phi_j(x_q)
1605 (fe_iv[to]->jump_in_shape_values(i, q_index, component_i) *
1606 scalar_product(j_numflux[face_no_j](component_i, component_j),
1607 normals[q_index])); // [[phi_i(x_q)]] * j_numflux(x_q, u_q)
1608 }
1609 }
1610 }
1611 };
1612 const auto copier = [&](const CopyData &c) {
1613 ldg_jacobian_tmp.add(c.local_dof_indices_to, c.local_dof_indices_from, c.cell_jacobian);
1614 for (auto &cdf : c.face_data)
1615 ldg_jacobian_tmp.add(cdf.joint_dof_indices_to, cdf.joint_dof_indices_from, cdf.cell_jacobian);
1616 };
1617
1618 const MeshWorker::AssembleFlags assemble_flags = MeshWorker::assemble_own_cells |
1619 MeshWorker::assemble_boundary_faces |
1620 MeshWorker::assemble_own_interior_faces_once;
1621 Scratch scratch_data(mapping, dof_handler_list, quadrature, quadrature_face);
1622 CopyData copy_data;
1623
1624 ldg_jacobian_tmp = 0;
1625 ldg_jacobian = 0;
1626 MeshWorker::mesh_loop(dof_handler.begin_active(), dof_handler.end(), cell_worker, copier, scratch_data,
1627 copy_data, assemble_flags, boundary_worker, face_worker, threads, batch_size);
1628 for (const auto &c : model.get_components().ldg_couplings(to, from))
1629 component_mass_matrix_inverse.mmult(ldg_jacobian.block(c[0], c[1]), ldg_jacobian_tmp.block(c[0], c[1]),
1630 Vector<NumberType>(), false);
1631 }
1632
1641 void build_ldg_sparsity(BlockSparsityPattern &sparsity_pattern, const DoFHandler<dim> &to_dofh,
1642 const DoFHandler<dim> &from_dofh, const int stencil = 1,
1643 bool add_extractor_dofs = false) const
1644 {
1645 const auto &triangulation = discretization.get_triangulation();
1646 auto to_dofs_per_component = DoFTools::count_dofs_per_fe_component(to_dofh);
1647 auto from_dofs_per_component = DoFTools::count_dofs_per_fe_component(from_dofh);
1648 auto to_n_fe = to_dofs_per_component.size();
1649 auto from_n_fe = from_dofs_per_component.size();
1650 for (uint j = 1; j < from_dofs_per_component.size(); ++j)
1651 if (from_dofs_per_component[j] != from_dofs_per_component[0])
1652 throw std::runtime_error("For LDG the FE basis of all systems must be equal!");
1653 for (uint j = 1; j < to_dofs_per_component.size(); ++j)
1654 if (to_dofs_per_component[j] != to_dofs_per_component[0])
1655 throw std::runtime_error("For LDG the FE basis of all systems must be equal!");
1656
1657 BlockDynamicSparsityPattern dsp(to_n_fe, from_n_fe);
1658 for (uint i = 0; i < to_n_fe; ++i)
1659 for (uint j = 0; j < from_n_fe; ++j)
1660 dsp.block(i, j).reinit(to_dofs_per_component[i], from_dofs_per_component[j]);
1661 dsp.collect_sizes();
1662
1663 const auto to_dofs_per_cell = to_dofh.get_fe().dofs_per_cell;
1664 const auto from_dofs_per_cell = from_dofh.get_fe().dofs_per_cell;
1665
1666 for (const auto &t_cell : triangulation.active_cell_iterators()) {
1667 std::vector<types::global_dof_index> to_dofs(to_dofs_per_cell);
1668 std::vector<types::global_dof_index> from_dofs(from_dofs_per_cell);
1669 const auto to_cell = t_cell->as_dof_handler_iterator(to_dofh);
1670 const auto from_cell = t_cell->as_dof_handler_iterator(from_dofh);
1671 to_cell->get_dof_indices(to_dofs);
1672 from_cell->get_dof_indices(from_dofs);
1673
1674 std::function<void(decltype(from_cell) &, const int)> add_all_neighbor_dofs = [&](const auto &from_cell,
1675 const int stencil_level) {
1676 for (const auto face_no : from_cell->face_indices()) {
1677 const auto face = from_cell->face(face_no);
1678 if (!face->at_boundary()) {
1679 auto neighbor_cell = from_cell->neighbor(face_no);
1680
1681 if (dim == 1)
1682 while (neighbor_cell->has_children())
1683 neighbor_cell = neighbor_cell->child(face_no == 0 ? 1 : 0);
1684
1685 // add all children
1686 else if (neighbor_cell->has_children()) {
1687 throw std::runtime_error("not yet implemented lol");
1688 }
1689
1690 if (!neighbor_cell->is_active()) continue;
1691
1692 std::vector<types::global_dof_index> tmp(from_dofs_per_cell);
1693 neighbor_cell->get_dof_indices(tmp);
1694
1695 from_dofs.insert(std::end(from_dofs), std::begin(tmp), std::end(tmp));
1696
1697 if (stencil_level < stencil) add_all_neighbor_dofs(neighbor_cell, stencil_level + 1);
1698 }
1699 }
1700 };
1701
1702 add_all_neighbor_dofs(from_cell, 1);
1703
1704 for (const auto i : to_dofs)
1705 for (const auto j : from_dofs)
1706 dsp.add(i, j);
1707 }
1708
1709 if (add_extractor_dofs)
1710 for (uint row = 0; row < dsp.n_rows(); ++row)
1711 dsp.add_row_entries(row, extractor_dof_indices);
1712
1713 sparsity_pattern.copy_from(dsp);
1714 }
1715
1722 void build_inverse(const SparseMatrix<NumberType> &in, SparseMatrix<NumberType> &out) const
1723 {
1724 GrowingVectorMemory<Vector<NumberType>> mem;
1725 SparseDirectUMFPACK inverse;
1726 inverse.initialize(in);
1727 // out is m x n
1728 // we go row-wise, i.e. we keep one n fixed and insert a row
1729 tbb::parallel_for(tbb::blocked_range<int>(0, out.n()), [&](tbb::blocked_range<int> r) {
1730 typename VectorMemory<Vector<NumberType>>::Pointer tmp(mem);
1731 tmp->reinit(out.m());
1732 for (int n = r.begin(); n < r.end(); ++n) {
1733 *tmp = 0;
1734 (*tmp)[n] = 1.;
1735 inverse.solve(*tmp);
1736 for (auto it = out.begin(n); it != out.end(n); ++it)
1737 it->value() = (*tmp)[it->column()];
1738 }
1739 });
1740 }
1741
1742 protected:
1743 constexpr static int nothing = 0;
1744 using Base::EoM;
1745 using Base::EoM_cell;
1746 using Base::extractor_jacobian_u;
1747 using Base::old_EoM_cell;
1748
1749 void readouts(DataOutput<dim, VectorType> &data_out, const VectorType &solution_global,
1750 const VectorType &variables) const
1751 {
1752 auto helper = [&](auto EoMfun, auto outputter) {
1753 auto EoM_cell = this->EoM_cell;
1754 auto EoM = get_EoM_point(
1755 EoM_cell, solution_global, dof_handler, mapping, EoMfun, [&](const auto &p, const auto &) { return p; },
1756 EoM_abs_tol, EoM_max_iter);
1757 auto EoM_unit = mapping.transform_real_to_unit_cell(EoM_cell, EoM);
1758
1759 using t_Iterator = typename Triangulation<dim>::active_cell_iterator;
1760
1761 std::vector<std::shared_ptr<FEValues<dim>>> fe_v;
1762 for (uint k = 0; k < Components::count_fe_subsystems(); ++k) {
1763 fe_v.emplace_back(std::make_shared<FEValues<dim>>(
1764 mapping, discretization.get_fe(k), EoM_unit,
1765 update_values | update_gradients | update_quadrature_points | update_JxW_values | update_hessians));
1766
1767 auto cell = dof_handler_list[k]->begin_active();
1768 cell->copy_from(*t_Iterator(EoM_cell));
1769 fe_v[k]->reinit(cell);
1770 }
1771
1772 std::vector<std::vector<Vector<NumberType>>> solutions;
1773 for (uint k = 0; k < Components::count_fe_subsystems(); ++k) {
1774 solutions.push_back({Vector<NumberType>(Components::count_fe_functions(k))});
1775 if (k == 0)
1776 fe_v[0]->get_function_values(solution_global, solutions[k]);
1777 else
1778 fe_v[k]->get_function_values(sol_vector[k], solutions[k]);
1779 }
1780 std::vector<Vector<NumberType>> solutions_vector;
1781 for (uint k = 0; k < Components::count_fe_subsystems(); ++k)
1782 solutions_vector.push_back(solutions[k][0]);
1783
1784 std::array<NumberType, Components::count_extractors()> __extracted_data{{}};
1785 if constexpr (Components::count_extractors() > 0)
1786 extract(__extracted_data, solution_global, variables, true, false, false);
1787 const auto &extracted_data = __extracted_data;
1788
1789 std::vector<std::vector<Tensor<1, dim, NumberType>>> solution_grad{
1790 std::vector<Tensor<1, dim, NumberType>>(Components::count_fe_functions())};
1791 std::vector<std::vector<Tensor<2, dim, NumberType>>> solution_hess{
1792 std::vector<Tensor<2, dim, NumberType>>(Components::count_fe_functions())};
1793 fe_v[0]->get_function_gradients(solution_global, solution_grad);
1794 fe_v[0]->get_function_hessians(solution_global, solution_hess);
1795
1796 auto solution_tuple = std::tuple_cat(vector_to_tuple<Components::count_fe_subsystems()>(solutions_vector),
1797 std::tie(solution_grad[0], solution_hess[0], extracted_data, variables));
1798
1799 outputter(data_out, EoM, fe_more_conv(solution_tuple));
1800 };
1801 model.template readouts_multiple(helper, data_out);
1802 }
1803
1804 void extract(std::array<NumberType, Components::count_extractors()> &data, const VectorType &solution_global,
1805 const VectorType &variables, bool search_EoM, bool set_EoM, bool postprocess) const
1806 {
1807 auto EoM = this->EoM;
1808 auto EoM_cell = this->EoM_cell;
1809 if (search_EoM || EoM_cell == *(dof_handler.active_cell_iterators().end()))
1810 EoM = get_EoM_point(
1811 EoM_cell, solution_global, dof_handler, mapping,
1812 [&](const auto &p, const auto &values) { return model.EoM(p, values); },
1813 [&](const auto &p, const auto &values) { return postprocess ? model.EoM_postprocess(p, values) : p; },
1814 EoM_abs_tol, EoM_max_iter);
1815 if (set_EoM) {
1816 this->EoM = EoM;
1817 this->EoM_cell = EoM_cell;
1818 }
1819 auto EoM_unit = mapping.transform_real_to_unit_cell(EoM_cell, EoM);
1820
1821 rebuild_ldg_vectors(solution_global);
1822
1823 using t_Iterator = typename Triangulation<dim>::active_cell_iterator;
1824
1825 std::vector<std::shared_ptr<FEValues<dim>>> fe_v;
1826 for (uint k = 0; k < Components::count_fe_subsystems(); ++k) {
1827 fe_v.emplace_back(std::make_shared<FEValues<dim>>(
1828 mapping, discretization.get_fe(k), EoM_unit,
1829 update_values | update_gradients | update_quadrature_points | update_JxW_values | update_hessians));
1830
1831 auto cell = dof_handler_list[k]->begin_active();
1832 cell->copy_from(*t_Iterator(EoM_cell));
1833 fe_v[k]->reinit(cell);
1834 }
1835
1836 std::vector<std::vector<Vector<NumberType>>> solutions;
1837 for (uint k = 0; k < Components::count_fe_subsystems(); ++k) {
1838 solutions.push_back({Vector<NumberType>(Components::count_fe_functions(k))});
1839 if (k == 0)
1840 fe_v[0]->get_function_values(solution_global, solutions[k]);
1841 else
1842 fe_v[k]->get_function_values(sol_vector[k], solutions[k]);
1843 }
1844 std::vector<Vector<NumberType>> solutions_vector;
1845 for (uint k = 0; k < Components::count_fe_subsystems(); ++k)
1846 solutions_vector.push_back(solutions[k][0]);
1847
1848 std::vector<std::vector<Tensor<1, dim, NumberType>>> solution_grad{
1849 std::vector<Tensor<1, dim, NumberType>>(Components::count_fe_functions())};
1850 std::vector<std::vector<Tensor<2, dim, NumberType>>> solution_hess{
1851 std::vector<Tensor<2, dim, NumberType>>(Components::count_fe_functions())};
1852 fe_v[0]->get_function_gradients(solution_global, solution_grad);
1853 fe_v[0]->get_function_hessians(solution_global, solution_hess);
1854
1855 auto solution_tuple = std::tuple_cat(vector_to_tuple<Components::count_fe_subsystems()>(solutions_vector),
1856 std::tie(solution_grad[0], solution_hess[0], this->nothing, variables));
1857
1858 model.template extract(data, EoM, fe_more_conv(solution_tuple));
1859 }
1860
1861 bool jacobian_extractors(FullMatrix<NumberType> &extractor_jacobian, const VectorType &solution_global,
1862 const VectorType &variables)
1863 {
1864 if (extractor_jacobian_u.m() != Components::count_extractors() ||
1865 extractor_jacobian_u.n() != Components::count_fe_functions())
1866 extractor_jacobian_u =
1867 FullMatrix<NumberType>(Components::count_extractors(), Components::count_fe_functions());
1868
1869 EoM = get_EoM_point(
1870 EoM_cell, solution_global, dof_handler, mapping,
1871 [&](const auto &p, const auto &values) { return model.EoM(p, values); },
1872 [&](const auto &p, const auto &values) { return model.EoM_postprocess(p, values); }, EoM_abs_tol,
1873 EoM_max_iter);
1874 auto EoM_unit = mapping.transform_real_to_unit_cell(EoM_cell, EoM);
1875 bool new_cell = (old_EoM_cell != EoM_cell);
1876 old_EoM_cell = EoM_cell;
1877
1878 using t_Iterator = typename Triangulation<dim>::active_cell_iterator;
1879
1880 std::vector<std::shared_ptr<FEValues<dim>>> fe_v;
1881 for (uint k = 0; k < Components::count_fe_subsystems(); ++k) {
1882 fe_v.emplace_back(std::make_shared<FEValues<dim>>(
1883 mapping, discretization.get_fe(k), EoM_unit,
1884 update_values | update_gradients | update_quadrature_points | update_JxW_values | update_hessians));
1885
1886 auto cell = dof_handler_list[k]->begin_active();
1887 cell->copy_from(*t_Iterator(EoM_cell));
1888 fe_v[k]->reinit(cell);
1889 }
1890
1891 const uint n_dofs = fe_v[0]->get_fe().n_dofs_per_cell();
1892 if (new_cell) {
1893 // spdlog::get("log")->info("FEM: Rebuilding the jacobian sparsity pattern");
1894 extractor_dof_indices.resize(n_dofs);
1895 EoM_cell->get_dof_indices(extractor_dof_indices);
1896 rebuild_jacobian_sparsity();
1897 }
1898
1899 std::vector<std::vector<Vector<NumberType>>> solutions;
1900 for (uint k = 0; k < Components::count_fe_subsystems(); ++k) {
1901 solutions.push_back({Vector<NumberType>(Components::count_fe_functions(k))});
1902 if (k == 0)
1903 fe_v[0]->get_function_values(solution_global, solutions[k]);
1904 else
1905 fe_v[k]->get_function_values(sol_vector[k], solutions[k]);
1906 }
1907 std::vector<Vector<NumberType>> solutions_vector;
1908 for (uint k = 0; k < Components::count_fe_subsystems(); ++k)
1909 solutions_vector.push_back(solutions[k][0]);
1910
1911 std::vector<std::vector<Tensor<1, dim, NumberType>>> solution_grad{
1912 std::vector<Tensor<1, dim, NumberType>>(Components::count_fe_functions())};
1913 std::vector<std::vector<Tensor<2, dim, NumberType>>> solution_hess{
1914 std::vector<Tensor<2, dim, NumberType>>(Components::count_fe_functions())};
1915 fe_v[0]->get_function_gradients(solution_global, solution_grad);
1916 fe_v[0]->get_function_hessians(solution_global, solution_hess);
1917
1918 auto solution_tuple = std::tuple_cat(vector_to_tuple<Components::count_fe_subsystems()>(solutions_vector),
1919 std::tie(solution_grad[0], solution_hess[0], this->nothing, variables));
1920
1921 extractor_jacobian_u = 0;
1922 model.template jacobian_extractors<0>(extractor_jacobian_u, EoM, fe_more_conv(solution_tuple));
1923
1924 if (extractor_jacobian.m() != Components::count_extractors() || extractor_jacobian.n() != n_dofs)
1925 extractor_jacobian = FullMatrix<NumberType>(Components::count_extractors(), n_dofs);
1926
1927 for (uint e = 0; e < Components::count_extractors(); ++e)
1928 for (uint i = 0; i < n_dofs; ++i) {
1929 const auto component_i = fe_v[0]->get_fe().system_to_component_index(i).first;
1930 extractor_jacobian(e, i) =
1931 extractor_jacobian_u(e, component_i) * fe_v[0]->shape_value_component(i, 0, component_i);
1932 }
1933
1934 return new_cell;
1935 }
1936
1937 using Base::timings_variable_jacobian;
1938 using Base::timings_variable_residual;
1939 template <typename... T> static constexpr auto v_tie(T &&...t)
1940 {
1941 return named_tuple<std::tuple<T &...>, "variables", "extractors">(std::tie(t...));
1942 }
1943
1944 template <typename... T> static constexpr auto e_tie(T &&...t)
1945 {
1946 return named_tuple<std::tuple<T &...>, "fe_functions", "fe_derivatives", "fe_hessians", "extractors",
1947 "variables">(std::tie(t...));
1948 }
1949
1950 virtual void residual_variables(VectorType &residual, const VectorType &variables,
1951 const VectorType &spatial_solution) override
1952 {
1953 Timer timer;
1954 std::array<NumberType, Components::count_extractors()> __extracted_data{{}};
1955 if constexpr (Components::count_extractors() > 0)
1956 extract(__extracted_data, spatial_solution, variables, true, false, false);
1957 const auto &extracted_data = __extracted_data;
1958 model.dt_variables(residual, v_tie(variables, extracted_data));
1959 timings_variable_residual.push_back(timer.wall_time());
1960 }
1961
1962 virtual void jacobian_variables(FullMatrix<NumberType> &jacobian, const VectorType &variables,
1963 const VectorType &spatial_solution) override
1964 {
1965 Timer timer;
1966 std::array<NumberType, Components::count_extractors()> __extracted_data{{}};
1967 if constexpr (Components::count_extractors() > 0)
1968 extract(__extracted_data, spatial_solution, variables, true, false, false);
1969 const auto &extracted_data = __extracted_data;
1970 model.template jacobian_variables<0>(jacobian, v_tie(variables, extracted_data));
1971 timings_variable_jacobian.push_back(timer.wall_time());
1972 }
1973 };
1974 } // namespace LDG
1975} // namespace DiFfRG
abstract_assembler.hh
DiFfRG::AbstractAssembler
This is the general assembler interface for any kind of discretization. An assembler is responsible f...
Definition abstract_assembler.hh:39
DiFfRG::DataOutput
Class to manage writing to files. FEM functions are written to vtk files and other data is written to...
Definition data_output.hh:20
DiFfRG::DataOutput::fe_output
FEOutput< dim, VectorType > & fe_output()
Returns a reference to the FEOutput object used to write FEM functions to .vtu files and ....
DiFfRG::JSONValue
A wrapper around the boost json value class.
Definition json.hh:19
DiFfRG::LDG::Assembler
The LDG assembler that can be used for any LDG scheme, with as many levels as one wants.
Definition ldg.hh:385
DiFfRG::LDG::Assembler::extract
void extract(std::array< NumberType, Components::count_extractors()> &data, const VectorType &solution_global, const VectorType &variables, bool search_EoM, bool set_EoM, bool postprocess) const
Definition ldg.hh:1804
DiFfRG::LDG::Assembler::average_time_residual_assembly
double average_time_residual_assembly() const
Definition ldg.hh:1250
DiFfRG::LDG::Assembler::sol_vector
array< BlockVector< NumberType >, Components::count_fe_subsystems()> sol_vector
Definition ldg.hh:1284
DiFfRG::LDG::Assembler::build_ldg_jacobian
void build_ldg_jacobian(const VectorType &solution_global, BlockSparseMatrix< NumberType > &ldg_jacobian, BlockSparseMatrix< NumberType > &ldg_jacobian_tmp) const
Build the LDG jacobian at level 'to', which takes information from level 'from'.
Definition ldg.hh:1489
DiFfRG::LDG::Assembler::fe
const FiniteElement< dim > & fe
Definition ldg.hh:126
DiFfRG::LDG::Assembler::component_mass_matrix_inverse
SparseMatrix< NumberType > component_mass_matrix_inverse
Definition ldg.hh:1297
DiFfRG::LDG::Assembler::sparsity_pattern_jacobian
BlockSparsityPattern sparsity_pattern_jacobian
Definition ldg.hh:1288
DiFfRG::LDG::Assembler::residual
virtual void residual(VectorType &residual, const VectorType &solution_global, NumberType weight, const VectorType &solution_global_dot, NumberType weight_mass, const VectorType &variables=VectorType()) override
Construct the system residual, i.e. Res = grad(flux) - source.
Definition ldg.hh:716
DiFfRG::LDG::Assembler::jacobian_mass
virtual void jacobian_mass(BlockSparseMatrix< NumberType > &jacobian, const VectorType &solution_global, const VectorType &solution_global_dot, NumberType alpha=1., NumberType beta=1.) override
Definition ldg.hh:864
DiFfRG::LDG::Assembler::num_residuals
uint num_residuals() const
Definition ldg.hh:1258
DiFfRG::LDG::Assembler::j_wg_tmp
array< BlockSparseMatrix< NumberType >, Components::count_fe_subsystems()> j_wg_tmp
Definition ldg.hh:1301
DiFfRG::LDG::Assembler::get_mass_matrix
virtual const BlockSparseMatrix< NumberType > & get_mass_matrix() const override
Obtain the mass matrix.
Definition ldg.hh:652
DiFfRG::LDG::Assembler::build_ldg_sparsity
void build_ldg_sparsity(BlockSparsityPattern &sparsity_pattern, const DoFHandler< dim > &to_dofh, const DoFHandler< dim > &from_dofh, const int stencil=1, bool add_extractor_dofs=false) const
Create a sparsity pattern for matrices between the DoFs of two DoFHandlers, with given stencil size.
Definition ldg.hh:1641
DiFfRG::LDG::Assembler::average_time_jacobian_assembly
double average_time_jacobian_assembly() const
Definition ldg.hh:1260
DiFfRG::LDG::Assembler::readouts
void readouts(DataOutput< dim, VectorType > &data_out, const VectorType &solution_global, const VectorType &variables) const
Definition ldg.hh:1749
DiFfRG::LDG::Assembler::NumberType
typename Discretization::NumberType NumberType
Definition ldg.hh:391
DiFfRG::LDG::Assembler::v_tie
static constexpr auto v_tie(T &&...t)
Definition ldg.hh:1939
DiFfRG::LDG::Assembler::mapping
const Mapping< dim > & mapping
Definition ldg.hh:128
DiFfRG::LDG::Assembler::timings_residual
std::vector< double > timings_residual
Definition ldg.hh:1304
DiFfRG::LDG::Assembler::sparsity_pattern_ug
array< BlockSparsityPattern, Components::count_fe_subsystems()> sparsity_pattern_ug
Definition ldg.hh:1290
DiFfRG::LDG::Assembler::sol_vector_tmp
array< BlockVector< NumberType >, Components::count_fe_subsystems()> sol_vector_tmp
Definition ldg.hh:1285
DiFfRG::LDG::Assembler::stencil
static constexpr uint stencil
Definition ldg.hh:396
DiFfRG::LDG::Assembler::fe_more_conv
auto fe_more_conv(std::tuple< T &... > &t) const
Definition ldg.hh:411
DiFfRG::LDG::Assembler::residual_variables
virtual void residual_variables(VectorType &residual, const VectorType &variables, const VectorType &spatial_solution) override
Definition ldg.hh:1950
DiFfRG::LDG::Assembler::VectorType
typename Discretization::VectorType VectorType
Definition ldg.hh:392
DiFfRG::LDG::Assembler::sparsity_pattern_mass
BlockSparsityPattern sparsity_pattern_mass
Definition ldg.hh:1289
DiFfRG::LDG::Assembler::reinit_vector
virtual void reinit_vector(VectorType &vec) const override
Definition ldg.hh:449
DiFfRG::LDG::Assembler::jacobian_variables
virtual void jacobian_variables(FullMatrix< NumberType > &jacobian, const VectorType &variables, const VectorType &spatial_solution) override
Definition ldg.hh:1962
DiFfRG::LDG::Assembler::Components
typename Discretization::Components Components
Definition ldg.hh:394
DiFfRG::LDG::Assembler::e_tie
static constexpr auto e_tie(T &&...t)
Definition ldg.hh:1944
DiFfRG::LDG::Assembler::num_jacobians
uint num_jacobians() const
Definition ldg.hh:1268
DiFfRG::LDG::Assembler::attach_data_output
virtual void attach_data_output(DataOutput< dim, VectorType > &data_out, const VectorType &solution, const VectorType &variables, const VectorType &dt_solution=VectorType(), const VectorType &residual=VectorType()) override
Attach all intermediate (ldg) vectors to the data output.
Definition ldg.hh:457
DiFfRG::LDG::Assembler::mass
virtual void mass(VectorType &residual, const VectorType &solution_global, const VectorType &solution_global_dot, NumberType weight) override
Construct the mass.
Definition ldg.hh:661
DiFfRG::LDG::Assembler::dof_handler_list
std::vector< const DoFHandler< dim > * > dof_handler_list
Definition ldg.hh:1282
DiFfRG::LDG::Assembler::fe_conv
auto fe_conv(std::tuple< T &... > &t) const
Definition ldg.hh:399
DiFfRG::LDG::Assembler::model
Model & model
Definition ldg.hh:125
DiFfRG::LDG::Assembler::quadrature_face
QGauss< dim - 1 > quadrature_face
Definition ldg.hh:1278
DiFfRG::LDG::Assembler::sparsity_pattern_wg
array< BlockSparsityPattern, Components::count_fe_subsystems()> sparsity_pattern_wg
Definition ldg.hh:1292
DiFfRG::LDG::Assembler::rebuild_ldg_jacobian
void rebuild_ldg_jacobian(const VectorType &sol) const
Definition ldg.hh:1330
DiFfRG::LDG::Assembler::dim
static constexpr uint dim
Definition ldg.hh:395
DiFfRG::LDG::Assembler::log
void log(const std::string logger)
Definition ldg.hh:1228
DiFfRG::LDG::Assembler::ldg_matrix_built
array< bool, Components::count_fe_subsystems()> ldg_matrix_built
Definition ldg.hh:1307
DiFfRG::LDG::Assembler::j_ug
array< BlockSparseMatrix< NumberType >, Components::count_fe_subsystems()> j_ug
Definition ldg.hh:1298
DiFfRG::LDG::Assembler::reinit
virtual void reinit() override
Reinitialize the assembler. This is necessary if the mesh has changed, e.g. after a mesh refinement.
Definition ldg.hh:482
DiFfRG::LDG::Assembler::rebuild_jacobian_sparsity
virtual void rebuild_jacobian_sparsity() override
Definition ldg.hh:559
DiFfRG::LDG::Assembler::ref_conv
auto ref_conv(std::tuple< T &... > &t) const
Definition ldg.hh:426
DiFfRG::LDG::Assembler::jacobian_extractors
bool jacobian_extractors(FullMatrix< NumberType > &extractor_jacobian, const VectorType &solution_global, const VectorType &variables)
Definition ldg.hh:1861
DiFfRG::LDG::Assembler::build_inverse
void build_inverse(const SparseMatrix< NumberType > &in, SparseMatrix< NumberType > &out) const
Build the inverse of matrix in and save the result to out.
Definition ldg.hh:1722
DiFfRG::LDG::Assembler::Discretization
Discretization_ Discretization
Definition ldg.hh:389
DiFfRG::LDG::Assembler::sparsity_pattern_gu
array< BlockSparsityPattern, Components::count_fe_subsystems()> sparsity_pattern_gu
Definition ldg.hh:1291
DiFfRG::LDG::Assembler::discretization
Discretization & discretization
Definition ldg.hh:124
DiFfRG::LDG::Assembler::batch_size
uint batch_size
Definition ldg.hh:131
DiFfRG::LDG::Assembler::threads
uint threads
Definition ldg.hh:130
DiFfRG::LDG::Assembler::jacobian_tmp
array< BlockSparseMatrix< NumberType >, Components::count_fe_subsystems()> jacobian_tmp
Definition ldg.hh:1294
DiFfRG::LDG::Assembler::get_sparsity_pattern_jacobian
virtual const BlockSparsityPattern & get_sparsity_pattern_jacobian() const override
Obtain the sparsity pattern of the jacobian matrix.
Definition ldg.hh:648
DiFfRG::LDG::Assembler::Assembler
Assembler(Discretization &discretization, Model &model, const JSONValue &json)
Definition ldg.hh:439
DiFfRG::LDG::Assembler::mass_matrix
BlockSparseMatrix< NumberType > mass_matrix
Definition ldg.hh:1296
DiFfRG::LDG::Assembler::jacobian
virtual void jacobian(BlockSparseMatrix< NumberType > &jacobian, const VectorType &solution_global, NumberType weight, const VectorType &solution_global_dot, NumberType alpha, NumberType beta, const VectorType &variables=VectorType()) override
Construct the system jacobian, i.e. dRes/du.
Definition ldg.hh:931
DiFfRG::LDG::Assembler::j_wg
array< BlockSparseMatrix< NumberType >, Components::count_fe_subsystems()> j_wg
Definition ldg.hh:1300
DiFfRG::LDG::Assembler::extractor_dof_indices
std::vector< types::global_dof_index > extractor_dof_indices
Definition ldg.hh:142
DiFfRG::LDG::Assembler::num_reinits
uint num_reinits() const
Definition ldg.hh:1248
DiFfRG::LDG::Assembler::rebuild_ldg_vectors
void rebuild_ldg_vectors(const VectorType &sol) const
Definition ldg.hh:1314
DiFfRG::LDG::Assembler::jacobian_tmp_built
array< bool, Components::count_fe_subsystems()> jacobian_tmp_built
Definition ldg.hh:1308
DiFfRG::LDG::Assembler::quadrature
QGauss< dim > quadrature
Definition ldg.hh:1277
DiFfRG::LDG::Assembler::timings_reinit
std::vector< double > timings_reinit
Definition ldg.hh:1303
DiFfRG::LDG::Assembler::sol_vector_vec_tmp
array< Vector< NumberType >, Components::count_fe_subsystems()> sol_vector_vec_tmp
Definition ldg.hh:1286
DiFfRG::LDG::Assembler::timings_jacobian
std::vector< double > timings_jacobian
Definition ldg.hh:1305
DiFfRG::LDG::Assembler::Model
Model_ Model
Definition ldg.hh:390
DiFfRG::LDG::Assembler::average_time_reinit
double average_time_reinit() const
Definition ldg.hh:1240
DiFfRG::LDG::Assembler::j_gu
array< BlockSparseMatrix< NumberType >, Components::count_fe_subsystems()> j_gu
Definition ldg.hh:1299
DiFfRG::LDG::Assembler::build_ldg_vector
void build_ldg_vector(const VectorType &solution_global, VectorTypeldg &ldg_vector, VectorTypeldg &ldg_vector_tmp) const
Build the LDG vector at level 'to', which takes information from level 'from'.
Definition ldg.hh:1358
DiFfRG::LDG::Assembler::refinement_indicator
virtual void refinement_indicator(Vector< double > &indicator, const VectorType &solution_global) override
Definition ldg.hh:566
DiFfRG::LDG::Assembler::dof_handler
const DoFHandler< dim > & dof_handler
Definition ldg.hh:127
DiFfRG::LDG::Discretization::NumberType
NumberType_ NumberType
Definition ldg.hh:33
DiFfRG::LDG::Discretization::Components
Components_ Components
Definition ldg.hh:32
DiFfRG::LDG::Discretization::dim
static constexpr uint dim
Definition ldg.hh:37
DiFfRG::LDG::Discretization::VectorType
Vector< NumberType > VectorType
Definition ldg.hh:34
DiFfRG::LDG::LDGAssemblerBase
Definition ldg.hh:41
DiFfRG::LDG::LDGAssemblerBase::extractor_jacobian_du
FullMatrix< NumberType > extractor_jacobian_du
Definition ldg.hh:140
DiFfRG::LDG::LDGAssemblerBase::EoM_abs_tol
const double EoM_abs_tol
Definition ldg.hh:135
DiFfRG::LDG::LDGAssemblerBase::fe
const FiniteElement< dim > & fe
Definition ldg.hh:126
DiFfRG::LDG::LDGAssemblerBase::num_variable_jacobians
uint num_variable_jacobians() const
Definition ldg.hh:121
DiFfRG::LDG::LDGAssemblerBase::extractor_jacobian_u
FullMatrix< NumberType > extractor_jacobian_u
Definition ldg.hh:139
DiFfRG::LDG::LDGAssemblerBase::extractor_jacobian_ddu
FullMatrix< NumberType > extractor_jacobian_ddu
Definition ldg.hh:141
DiFfRG::LDG::LDGAssemblerBase::NumberType
typename Discretization::NumberType NumberType
Definition ldg.hh:45
DiFfRG::LDG::LDGAssemblerBase::get_discretization
auto & get_discretization()
Definition ldg.hh:71
DiFfRG::LDG::LDGAssemblerBase::average_time_variable_residual_assembly
double average_time_variable_residual_assembly()
Definition ldg.hh:103
DiFfRG::LDG::LDGAssemblerBase::mapping
const Mapping< dim > & mapping
Definition ldg.hh:128
DiFfRG::LDG::LDGAssemblerBase::dim
static constexpr uint dim
Definition ldg.hh:49
DiFfRG::LDG::LDGAssemblerBase::extractor_jacobian
FullMatrix< NumberType > extractor_jacobian
Definition ldg.hh:138
DiFfRG::LDG::LDGAssemblerBase::old_EoM_cell
DoFHandler< dim >::cell_iterator old_EoM_cell
Definition ldg.hh:134
DiFfRG::LDG::LDGAssemblerBase::Discretization
Discretization_ Discretization
Definition ldg.hh:43
DiFfRG::LDG::LDGAssemblerBase::VectorType
typename Discretization::VectorType VectorType
Definition ldg.hh:46
DiFfRG::LDG::LDGAssemblerBase::set_time
virtual void set_time(double t) override
Set the current time. The assembler should usually just forward this to the numerical model.
Definition ldg.hh:99
DiFfRG::LDG::LDGAssemblerBase::get_discretization
const auto & get_discretization() const
Definition ldg.hh:70
DiFfRG::LDG::LDGAssemblerBase::model
Model & model
Definition ldg.hh:125
DiFfRG::LDG::LDGAssemblerBase::refinement_indicator
virtual void refinement_indicator(Vector< double > &, const VectorType &)=0
DiFfRG::LDG::LDGAssemblerBase::reinit
virtual void reinit() override
Reinitialize the assembler. This is necessary if the mesh has changed, e.g. after a mesh refinement.
Definition ldg.hh:73
DiFfRG::LDG::LDGAssemblerBase::average_time_variable_jacobian_assembly
double average_time_variable_jacobian_assembly()
Definition ldg.hh:113
DiFfRG::LDG::LDGAssemblerBase::EoM_max_iter
const uint EoM_max_iter
Definition ldg.hh:136
DiFfRG::LDG::LDGAssemblerBase::get_differential_indices
virtual IndexSet get_differential_indices() const override
Obtain the dofs which contain time derivatives.
Definition ldg.hh:64
DiFfRG::LDG::LDGAssemblerBase::timings_variable_residual
std::vector< double > timings_variable_residual
Definition ldg.hh:144
DiFfRG::LDG::LDGAssemblerBase::timings_variable_jacobian
std::vector< double > timings_variable_jacobian
Definition ldg.hh:145
DiFfRG::LDG::LDGAssemblerBase::Components
typename Discretization::Components Components
Definition ldg.hh:48
DiFfRG::LDG::LDGAssemblerBase::discretization
Discretization & discretization
Definition ldg.hh:124
DiFfRG::LDG::LDGAssemblerBase::batch_size
uint batch_size
Definition ldg.hh:131
DiFfRG::LDG::LDGAssemblerBase::threads
uint threads
Definition ldg.hh:130
DiFfRG::LDG::LDGAssemblerBase::LDGAssemblerBase
LDGAssemblerBase(Discretization &discretization, Model &model, const JSONValue &json)
Definition ldg.hh:51
DiFfRG::LDG::LDGAssemblerBase::EoM
Point< dim > EoM
Definition ldg.hh:137
DiFfRG::LDG::LDGAssemblerBase::EoM_cell
DoFHandler< dim >::cell_iterator EoM_cell
Definition ldg.hh:133
DiFfRG::LDG::LDGAssemblerBase::extractor_dof_indices
std::vector< types::global_dof_index > extractor_dof_indices
Definition ldg.hh:142
DiFfRG::LDG::LDGAssemblerBase::num_variable_residuals
uint num_variable_residuals() const
Definition ldg.hh:111
DiFfRG::LDG::LDGAssemblerBase::Model
Model_ Model
Definition ldg.hh:44
DiFfRG::LDG::LDGAssemblerBase::rebuild_jacobian_sparsity
virtual void rebuild_jacobian_sparsity()=0
DiFfRG::LDG::LDGAssemblerBase::dof_handler
const DoFHandler< dim > & dof_handler
Definition ldg.hh:127
DiFfRG::Quadrature
Definition quadrature.hh:50
DiFfRG::Quadrature::size
size_t size() const
DiFfRG::SimpleMatrix
A simple NxM-matrix class, which is used for cell-wise Jacobians.
Definition tuples.hh:153
data_output.hh
DiFfRG
Definition complex_math.hh:14
DiFfRG::get_EoM_point
dealii::Point< dim > get_EoM_point(typename dealii::DoFHandler< dim >::cell_iterator &EoM_cell, const VectorType &sol, const dealii::DoFHandler< dim > &dof_handler, const dealii::Mapping< dim > &mapping, const EoMFUN &get_EoM, const EoMPFUN &EoM_postprocess=[](const auto &p, const auto &values) { return p;}, const double EoM_abs_tol=1e-5, const uint max_iter=100)
Get the EoM point for a given solution and model.
Definition EoM.hh:662
DiFfRG::local_sol_q
auto local_sol_q(const std::array< T, N > &a, uint q_index)
Definition tuples.hh:278
DiFfRG::is_close
bool __forceinline__ __host__ __device__ is_close(T1 a, T2 b, T3 eps_)
Function to evaluate whether two floats are equal to numerical precision. Tests for both relative and...
Definition math.hh:160
DiFfRG::jacobian_tuple
auto jacobian_tuple()
Definition tuples.hh:288
DiFfRG::vector_to_tuple
auto vector_to_tuple(const std::vector< T > &v)
Definition tuples.hh:216
DiFfRG::constexpr_for
constexpr void constexpr_for(F &&f)
A compile-time for loop, which calls the lambda f of signature void(integer) for each index.
Definition utils.hh:27
DiFfRG::jacobian_2_tuple
auto jacobian_2_tuple()
Definition tuples.hh:299
DiFfRG::uint
unsigned int uint
Definition utils.hh:22
std::get
constexpr auto & get(DiFfRG::named_tuple< tuple_type, strs... > &ob)
Definition tuples.hh:119
DiFfRG::LDG::internal::CopyData_I::CopyFaceData_I
Definition ldg.hh:367
DiFfRG::LDG::internal::CopyData_I::CopyFaceData_I::cell_indices
std::array< uint, 2 > cell_indices
Definition ldg.hh:368
DiFfRG::LDG::internal::CopyData_I::CopyFaceData_I::values
std::array< double, 2 > values
Definition ldg.hh:369
DiFfRG::LDG::internal::CopyData_I
Definition ldg.hh:366
DiFfRG::LDG::internal::CopyData_I::cell_index
uint cell_index
Definition ldg.hh:373
DiFfRG::LDG::internal::CopyData_I::value
double value
Definition ldg.hh:372
DiFfRG::LDG::internal::CopyData_I::face_data
std::vector< CopyFaceData_I > face_data
Definition ldg.hh:371
DiFfRG::LDG::internal::CopyData_J::CopyDataFace_J
Definition ldg.hh:288
DiFfRG::LDG::internal::CopyData_J::CopyDataFace_J::cell_jacobian
FullMatrix< NumberType > cell_jacobian
Definition ldg.hh:289
DiFfRG::LDG::internal::CopyData_J::CopyDataFace_J::joint_dof_indices_from
std::vector< types::global_dof_index > joint_dof_indices_from
Definition ldg.hh:290
DiFfRG::LDG::internal::CopyData_J::CopyDataFace_J::joint_dof_indices_to
std::vector< types::global_dof_index > joint_dof_indices_to
Definition ldg.hh:291
DiFfRG::LDG::internal::CopyData_J_full::CopyDataFace_J
Definition ldg.hh:322
DiFfRG::LDG::internal::CopyData_J_full::CopyDataFace_J::extractor_cell_jacobian
FullMatrix< NumberType > extractor_cell_jacobian
Definition ldg.hh:324
DiFfRG::LDG::internal::CopyData_J_full::CopyDataFace_J::cell_jacobian
array< FullMatrix< NumberType >, n_fe_subsystems > cell_jacobian
Definition ldg.hh:323
DiFfRG::LDG::internal::CopyData_J_full::CopyDataFace_J::joint_dof_indices
array< vector< types::global_dof_index >, n_fe_subsystems > joint_dof_indices
Definition ldg.hh:325
DiFfRG::LDG::internal::CopyData_J_full
Definition ldg.hh:321
DiFfRG::LDG::internal::CopyData_J_full::reinit
void reinit(const array< Iterator, n_fe_subsystems > &cell, const uint n_extractors)
Definition ldg.hh:336
DiFfRG::LDG::internal::CopyData_J_full::new_face_data
CopyDataFace_J & new_face_data(const array< unique_ptr< FEInterfaceValues< dim > >, n_fe_subsystems > &fe_iv, const uint n_extractors)
Definition ldg.hh:351
DiFfRG::LDG::internal::CopyData_J_full::cell_mass_jacobian
FullMatrix< NumberType > cell_mass_jacobian
Definition ldg.hh:329
DiFfRG::LDG::internal::CopyData_J_full::cell_jacobian
array< FullMatrix< NumberType >, n_fe_subsystems > cell_jacobian
Definition ldg.hh:328
DiFfRG::LDG::internal::CopyData_J_full::face_data
vector< CopyDataFace_J > face_data
Definition ldg.hh:332
DiFfRG::LDG::internal::CopyData_J_full::extractor_cell_jacobian
FullMatrix< NumberType > extractor_cell_jacobian
Definition ldg.hh:330
DiFfRG::LDG::internal::CopyData_J_full::local_dof_indices
array< vector< types::global_dof_index >, n_fe_subsystems > local_dof_indices
Definition ldg.hh:331
DiFfRG::LDG::internal::CopyData_J_full::dofs_per_cell
uint dofs_per_cell
Definition ldg.hh:334
DiFfRG::LDG::internal::CopyData_J
Definition ldg.hh:287
DiFfRG::LDG::internal::CopyData_J::new_face_data
CopyDataFace_J & new_face_data(const FEInterfaceValues< dim > &fe_iv_from, const FEInterfaceValues< dim > &fe_iv_to)
Definition ldg.hh:310
DiFfRG::LDG::internal::CopyData_J::local_dof_indices_to
std::vector< types::global_dof_index > local_dof_indices_to
Definition ldg.hh:296
DiFfRG::LDG::internal::CopyData_J::cell_jacobian
FullMatrix< NumberType > cell_jacobian
Definition ldg.hh:294
DiFfRG::LDG::internal::CopyData_J::face_data
std::vector< CopyDataFace_J > face_data
Definition ldg.hh:297
DiFfRG::LDG::internal::CopyData_J::local_dof_indices_from
std::vector< types::global_dof_index > local_dof_indices_from
Definition ldg.hh:295
DiFfRG::LDG::internal::CopyData_J::reinit
void reinit(const Iterator &cell_from, const Iterator &cell_to, uint dofs_per_cell_from, uint dofs_per_cell_to)
Definition ldg.hh:300
DiFfRG::LDG::internal::CopyData_R::CopyDataFace_R
Definition ldg.hh:259
DiFfRG::LDG::internal::CopyData_R::CopyDataFace_R::cell_residual
Vector< NumberType > cell_residual
Definition ldg.hh:260
DiFfRG::LDG::internal::CopyData_R::CopyDataFace_R::joint_dof_indices
std::vector< types::global_dof_index > joint_dof_indices
Definition ldg.hh:261
DiFfRG::LDG::internal::CopyData_R
Definition ldg.hh:258
DiFfRG::LDG::internal::CopyData_R::face_data
std::vector< CopyDataFace_R > face_data
Definition ldg.hh:267
DiFfRG::LDG::internal::CopyData_R::reinit
void reinit(const Iterator &cell, uint dofs_per_cell)
Definition ldg.hh:269
DiFfRG::LDG::internal::CopyData_R::cell_mass
Vector< NumberType > cell_mass
Definition ldg.hh:265
DiFfRG::LDG::internal::CopyData_R::local_dof_indices
std::vector< types::global_dof_index > local_dof_indices
Definition ldg.hh:266
DiFfRG::LDG::internal::CopyData_R::cell_residual
Vector< NumberType > cell_residual
Definition ldg.hh:264
DiFfRG::LDG::internal::CopyData_R::new_face_data
CopyDataFace_R & new_face_data(const FEInterfaceValues< dim > &fe_iv)
Definition ldg.hh:277
DiFfRG::LDG::internal::ScratchData
Class to hold data for each assembly thread, i.e. FEValues for cells, interfaces, as well as pre-allo...
Definition ldg.hh:154
DiFfRG::LDG::internal::ScratchData::n_components
array< uint, n_fe_subsystems > n_components
Definition ldg.hh:245
DiFfRG::LDG::internal::ScratchData::NumberType
typename Discretization::NumberType NumberType
Definition ldg.hh:156
DiFfRG::LDG::internal::ScratchData::solution_dot
vector< Vector< NumberType > > solution_dot
Definition ldg.hh:254
DiFfRG::LDG::internal::ScratchData::n_fe_subsystems
static constexpr uint n_fe_subsystems
Definition ldg.hh:157
DiFfRG::LDG::internal::ScratchData::new_fe_values
const auto & new_fe_values(const t_Iterator &t_cell)
Definition ldg.hh:218
DiFfRG::LDG::internal::ScratchData::dim
static constexpr uint dim
Definition ldg.hh:155
DiFfRG::LDG::internal::ScratchData::cell
array< Iterator, n_fe_subsystems > cell
Definition ldg.hh:246
DiFfRG::LDG::internal::ScratchData::fe_boundary_values
array< unique_ptr< FEFaceValues< dim > >, n_fe_subsystems > fe_boundary_values
Definition ldg.hh:251
DiFfRG::LDG::internal::ScratchData::t_Iterator
typename Triangulation< dim >::active_cell_iterator t_Iterator
Definition ldg.hh:159
DiFfRG::LDG::internal::ScratchData::ScratchData
ScratchData(const ScratchData< Discretization > &scratch_data)
Definition ldg.hh:191
DiFfRG::LDG::internal::ScratchData::Iterator
typename DoFHandler< dim >::active_cell_iterator Iterator
Definition ldg.hh:158
DiFfRG::LDG::internal::ScratchData::fe_interface_values
array< unique_ptr< FEInterfaceValues< dim > >, n_fe_subsystems > fe_interface_values
Definition ldg.hh:250
DiFfRG::LDG::internal::ScratchData::solution_interface
array< array< vector< Vector< NumberType > >, n_fe_subsystems >, 2 > solution_interface
Definition ldg.hh:255
DiFfRG::LDG::internal::ScratchData::ScratchData
ScratchData(const Mapping< dim > &mapping, const vector< const DoFHandler< dim > * > &dofh, const Quadrature< dim > &quadrature, const Quadrature< dim - 1 > &quadrature_face, const UpdateFlags update_flags=update_values|update_gradients|update_quadrature_points|update_JxW_values, const UpdateFlags interface_update_flags=update_values|update_gradients|update_quadrature_points|update_JxW_values|update_normal_vectors)
Definition ldg.hh:161
DiFfRG::LDG::internal::ScratchData::new_fe_interface_values
const auto & new_fe_interface_values(const t_Iterator &t_cell, uint f, uint sf, const t_Iterator &t_ncell, uint nf, unsigned int nsf)
Definition ldg.hh:226
DiFfRG::LDG::internal::ScratchData::ncell
array< Iterator, n_fe_subsystems > ncell
Definition ldg.hh:247
DiFfRG::LDG::internal::ScratchData::fe_values
array< unique_ptr< FEValues< dim > >, n_fe_subsystems > fe_values
Definition ldg.hh:249
DiFfRG::LDG::internal::ScratchData::solution
array< vector< Vector< NumberType > >, n_fe_subsystems > solution
Definition ldg.hh:253
DiFfRG::LDG::internal::ScratchData::new_fe_boundary_values
const auto & new_fe_boundary_values(const t_Iterator &t_cell, uint face_no)
Definition ldg.hh:236
DiFfRG::named_tuple
A class to store a tuple with elements that can be accessed by name. The names are stored as FixedStr...
Definition tuples.hh:27