Loading [MathJax]/extensions/tex2jax.js
OGS
FluxCorrectedTransport.cpp
Go to the documentation of this file.
1
10#include "FluxCorrectedTransport.h"
11
12#include <Eigen/Core>
13#include <Eigen/Sparse>
14#include <ranges>
15
16#ifdef USE_PETSC
17#include <petscerror.h>
18#include <petscmat.h>
19#include <petscvec.h>
20
21#include "BaseLib/MPI.h"
22#endif
23
24namespace NumLib
25{
26namespace detail
27{
28
29template <typename MatrixVectorType>
30std::unique_ptr<MatrixVectorType> newZeroedInstance(
31 MathLib::MatrixSpecifications const& matrix_specification)
32{
33 auto result = MathLib::MatrixVectorTraits<MatrixVectorType>::newInstance(
34 matrix_specification);
35 result->setZero();
36 return result;
37}
38
39#ifndef USE_PETSC
40void calculateFluxCorrectedTransportSerial(
41 [[maybe_unused]] const double t, const double dt,
42 std::vector<GlobalVector*> const& x,
43 std::vector<GlobalVector*> const& x_prev, int const process_id,
44 const MathLib::MatrixSpecifications& matrix_specification, GlobalMatrix& M,
45 GlobalMatrix& K, GlobalVector& b)
46{
47 auto D = newZeroedInstance<GlobalMatrix>(matrix_specification);
48 auto F = newZeroedInstance<GlobalMatrix>(matrix_specification);
49
50 // compute artificial diffusion operator D
51 using RawMatrixType = Eigen::SparseMatrix<double, Eigen::RowMajor>;
52 auto& K_raw = K.getRawMatrix();
53 K_raw *= -1;
54 for (int k = 0; k < K_raw.outerSize(); ++k)
55 {
56 for (RawMatrixType::InnerIterator it(K_raw, k); it; ++it)
57 {
58 if (it.row() != it.col())
59 {
60 double const kij = it.value();
61 double const kji = K.get(it.col(), it.row());
62 double const dij = std::max({-kij, 0., -kji});
63 D->setValue(it.row(), it.col(), dij);
64 }
65 }
66 }
67 auto& D_raw = D->getRawMatrix();
68 D_raw -= (D_raw * Eigen::VectorXd::Ones(D_raw.cols())).eval().asDiagonal();
69
70 // compute F
71 for (int k = 0; k < D_raw.outerSize(); ++k)
72 {
73 for (RawMatrixType::InnerIterator it(D_raw, k); it; ++it)
74 {
75 double const dij = it.value();
76 double const xj = x[process_id]->get(it.col());
77 double const xi = x[process_id]->get(it.row());
78 double const fij = -dij * (xj - xi);
79 F->setValue(it.row(), it.col(), fij);
80 }
81 }
82
83 auto& M_raw = M.getRawMatrix();
84 for (int k = 0; k < M_raw.outerSize(); ++k)
85 {
86 for (RawMatrixType::InnerIterator it(M_raw, k); it; ++it)
87 {
88 double const mij = it.value();
89 double const xdotj = (x[process_id]->get(it.col()) -
90 x_prev[process_id]->get(it.col())) /
91 dt;
92 double const xdoti = (x[process_id]->get(it.row()) -
93 x_prev[process_id]->get(it.row())) /
94 dt;
95 double const fij = -mij * (xdotj - xdoti);
96 F->add(it.row(), it.col(), fij);
97 }
98 }
99
100 auto P_plus = newZeroedInstance<GlobalVector>(matrix_specification);
101 auto P_minus = newZeroedInstance<GlobalVector>(matrix_specification);
102 auto Q_plus = newZeroedInstance<GlobalVector>(matrix_specification);
103 auto Q_minus = newZeroedInstance<GlobalVector>(matrix_specification);
104 auto R_plus = newZeroedInstance<GlobalVector>(matrix_specification);
105 auto R_minus = newZeroedInstance<GlobalVector>(matrix_specification);
106
107 auto& F_raw = F->getRawMatrix();
108 for (int k = 0; k < F_raw.outerSize(); ++k)
109 {
110 for (RawMatrixType::InnerIterator it(F_raw, k); it; ++it)
111 {
112 if (it.row() != it.col())
113 {
114 double const fij = it.value();
115 P_plus->add(it.row(), std::max(0., fij));
116 P_minus->add(it.row(), std::min(0., fij));
117
118 double const x_prev_i = x_prev[process_id]->get(it.row());
119 double const x_prev_j = x_prev[process_id]->get(it.col());
120
121 double const Q_plus_i = Q_plus->get(it.row());
122 double Q_plus_i_tmp =
123 std::max({Q_plus_i, 0., x_prev_j - x_prev_i});
124 Q_plus->set(it.row(), Q_plus_i_tmp);
125
126 double const Q_minus_i = Q_minus->get(it.row());
127 double Q_minus_i_tmp =
128 std::min({Q_minus_i, 0., x_prev_j - x_prev_i});
129 Q_minus->set(it.row(), Q_minus_i_tmp);
130 }
131 }
132 }
133
134 Eigen::VectorXd const M_L =
135 (M_raw * Eigen::VectorXd::Ones(M_raw.cols())).eval();
136 for (auto k = R_plus->getRangeBegin(); k < R_plus->getRangeEnd(); ++k)
137 {
138 double const P_plus_i = P_plus->get(k);
139
140 if (P_plus_i == 0.0)
141 {
142 continue;
143 }
144
145 double const mi = M_L(k);
146 double const Q_plus_i = Q_plus->get(k);
147 R_plus->set(k, std::min(1.0, mi * Q_plus_i / dt / P_plus_i));
148 }
149
150 for (auto k = R_minus->getRangeBegin(); k < R_minus->getRangeEnd(); ++k)
151 {
152 double const P_minus_i = P_minus->get(k);
153 if (P_minus_i == 0.0)
154 {
155 continue;
156 }
157 double const mi = M_L(k);
158 double const Q_minus_i = Q_minus->get(k);
159 R_minus->set(k, std::min(1.0, mi * Q_minus_i / dt / P_minus_i));
160 }
161
162 for (int k = 0; k < F_raw.outerSize(); ++k)
163 {
164 for (RawMatrixType::InnerIterator it(F_raw, k); it; ++it)
165 {
166 double const fij = it.value();
167 if (fij > 0.)
168 {
169 double const R_plus_i = R_plus->get(it.row());
170 double const R_minus_j = R_minus->get(it.col());
171 double const alpha_ij = std::min(R_plus_i, R_minus_j);
172
173 b.add(it.row(), alpha_ij * fij);
174 }
175 else
176 {
177 double const R_minus_i = R_minus->get(it.row());
178 double const R_plus_j = R_plus->get(it.col());
179 double const alpha_ij = std::min(R_minus_i, R_plus_j);
180
181 b.add(it.row(), alpha_ij * fij);
182 }
183 }
184 }
185
186 // compute low-order operator
187 K_raw += D_raw;
188 K_raw *= -1;
189
190 // overwrite with the lumped mass matrix
191 M.setZero();
192 for (int k = 0; k < M.getNumberOfRows(); ++k)
193 {
194 M.setValue(k, k, M_L(k));
195 }
196}
197#endif // end of ifndef USE_PETSC
198
199#ifdef USE_PETSC
200void finalize(Mat& M)
201{
202 PetscCallAbort(PETSC_COMM_WORLD,
203 MatAssemblyBegin(M, MatAssemblyType::MAT_FINAL_ASSEMBLY));
204 PetscCallAbort(PETSC_COMM_WORLD,
205 MatAssemblyEnd(M, MatAssemblyType::MAT_FINAL_ASSEMBLY));
206}
207
208void calculateFluxCorrectedTransportPETSc(
209 [[maybe_unused]] const double t, const double dt,
210 std::vector<GlobalVector*> const& x,
211 std::vector<GlobalVector*> const& x_prev, int const process_id,
212 const MathLib::MatrixSpecifications& matrix_specification, GlobalMatrix& M,
213 GlobalMatrix& K, GlobalVector& b)
214{
215 BaseLib::MPI::Mpi mpi(PETSC_COMM_WORLD);
216 if (mpi.size > 1)
217 {
218 OGS_FATAL(
219 "The current parallel implementation of flux corrected transport "
220 "supports only one MPI rank, but the simulation uses {} ranks",
221 mpi.size);
222 }
223 Mat K_raw = K.getRawMatrix();
224 Mat D;
225 PetscCallAbort(PETSC_COMM_WORLD,
226 MatDuplicate(K_raw, MAT_DO_NOT_COPY_VALUES, &D));
227 PetscCallAbort(PETSC_COMM_WORLD, MatZeroEntries(D));
228 finalize(D);
229 PetscInt start, end;
230 PetscCallAbort(PETSC_COMM_WORLD,
231 MatGetOwnershipRange(K.getRawMatrix(), &start, &end));
232
233 PetscCallAbort(PETSC_COMM_WORLD, MatScale(K_raw, -1.0));
234 K.finalizeAssembly();
235 // compute artificial diffusion operator D
236 for (PetscInt i = start; i < end; ++i)
237 {
238 PetscInt number_of_columns;
239 PetscInt const* columns;
240 PetscScalar const* values;
241
242 PetscCallAbort(PETSC_COMM_WORLD, MatGetRow(K_raw, i, &number_of_columns,
243 &columns, &values));
244
245 auto const columns_span =
246 std::span<PetscInt const>(columns, number_of_columns);
247 auto const values_span =
248 std::span<PetscScalar const>(values, number_of_columns);
249 for (auto const [j, kij] : std::views::zip(columns_span, values_span))
250 {
251 if (i == j) // skip diagonal entries
252 {
253 continue;
254 }
255 PetscScalar kji;
256 // MatGetValue is not supported for all PETSc matrix types
257 PetscCallAbort(PETSC_COMM_WORLD, MatGetValue(K_raw, j, i, &kji));
258 PetscScalar dij = std::max({-kij, 0.0, -kji});
259 PetscCallAbort(PETSC_COMM_WORLD,
260 MatSetValue(D, i, j, dij, INSERT_VALUES));
261 }
262 // MatRestoreRow() is required after MatGetRow()
263 PetscCallAbort(
264 PETSC_COMM_WORLD,
265 MatRestoreRow(K_raw, i, &number_of_columns, &columns, &values));
266 }
267 finalize(D);
268 auto row_sums = newZeroedInstance<GlobalVector>(matrix_specification);
269 Vec row_sums_raw = row_sums->getRawVector();
270 PetscCallAbort(PETSC_COMM_WORLD, MatGetRowSum(D, row_sums_raw));
271 PetscCallAbort(PETSC_COMM_WORLD, VecScale(row_sums_raw, -1.0));
272 PetscCallAbort(PETSC_COMM_WORLD, VecAssemblyBegin(row_sums_raw));
273 PetscCallAbort(PETSC_COMM_WORLD, VecAssemblyEnd(row_sums_raw));
274 PetscCallAbort(PETSC_COMM_WORLD,
275 MatDiagonalSet(D, row_sums_raw, ADD_VALUES));
276 finalize(D);
277 // need to access x and x_prev in the following loops
278 x[process_id]->setLocalAccessibleVector();
279 x_prev[process_id]->setLocalAccessibleVector();
280
281 // compute F
282 Mat F;
283 PetscCallAbort(PETSC_COMM_WORLD,
284 MatDuplicate(K.getRawMatrix(), MAT_DO_NOT_COPY_VALUES, &F));
285 PetscCallAbort(PETSC_COMM_WORLD, MatZeroEntries(F));
286 for (PetscInt i = start; i < end; ++i)
287 {
288 PetscInt cols;
289 const PetscInt* col_indices;
290 const PetscScalar* values;
291
292 PetscCallAbort(PETSC_COMM_WORLD,
293 MatGetRow(D, i, &cols, &col_indices, &values));
294
295 auto const columns = std::span<PetscInt const>(col_indices, cols);
296 auto const values_D = std::span<PetscScalar const>(values, cols);
297
298 for (auto const [j, dij] : std::views::zip(columns, values_D))
299 {
300 PetscScalar const xi = x[process_id]->get(i);
301 PetscScalar const xj = x[process_id]->get(j);
302 PetscScalar const fij = -dij * (xj - xi);
303 PetscCallAbort(PETSC_COMM_WORLD,
304 MatSetValue(F, i, j, fij, INSERT_VALUES));
305 }
306 PetscCallAbort(PETSC_COMM_WORLD,
307 MatRestoreRow(D, i, &cols, &col_indices, &values));
308 }
309 finalize(F);
310 auto& M_raw = M.getRawMatrix();
311
312 auto xdot = newZeroedInstance<GlobalVector>(matrix_specification);
313 xdot->setLocalAccessibleVector();
314 for (PetscInt i = start; i < end; ++i)
315 {
316 xdot->add(i, (x[process_id]->get(i) - x_prev[process_id]->get(i)) / dt);
317 }
318 xdot->finalizeAssembly();
319 xdot->setLocalAccessibleVector();
320
321 for (PetscInt i = start; i < end; ++i)
322 {
323 PetscInt cols;
324 PetscInt const* col_indices;
325 PetscScalar const* values;
326
327 PetscCallAbort(PETSC_COMM_WORLD,
328 MatGetRow(M_raw, i, &cols, &col_indices, &values));
329
330 auto const columns = std::span<PetscInt const>(col_indices, cols);
331 auto const values_M = std::span<PetscScalar const>(values, cols);
332
333 for (auto const [j, mij] : std::views::zip(columns, values_M))
334 {
335 PetscScalar const fij = -mij * (xdot->get(j) - xdot->get(i));
336 PetscCallAbort(PETSC_COMM_WORLD,
337 MatSetValue(F, i, j, fij, ADD_VALUES));
338 }
339 PetscCallAbort(PETSC_COMM_WORLD,
340 MatRestoreRow(M_raw, i, &cols, &col_indices, &values));
341 }
342 finalize(F);
343 auto P_plus = newZeroedInstance<GlobalVector>(matrix_specification);
344 auto P_minus = newZeroedInstance<GlobalVector>(matrix_specification);
345
346 P_plus->setLocalAccessibleVector();
347 P_minus->setLocalAccessibleVector();
348 // Fill P_plus, P_minus, Q_plus, Q_minus
349 for (PetscInt i = start; i < end; ++i)
350 {
351 PetscInt cols;
352 PetscInt const* col_indices;
353 PetscScalar const* values;
354
355 PetscCallAbort(PETSC_COMM_WORLD,
356 MatGetRow(F, i, &cols, &col_indices, &values));
357
358 auto const column_indices =
359 std::span<PetscInt const>(col_indices, cols);
360 auto const values_span = std::span<PetscScalar const>(values, cols);
361 for (auto const [j, fij] : std::views::zip(column_indices, values_span))
362 {
363 // skip diagonal
364 if (i == j)
365 {
366 continue;
367 }
368 P_plus->add(i, std::max(0., fij));
369 P_minus->add(i, std::min(0., fij));
370 }
371 PetscCallAbort(PETSC_COMM_WORLD,
372 MatRestoreRow(F, i, &cols, &col_indices, &values));
373 }
374 P_plus->finalizeAssembly();
375 P_plus->setLocalAccessibleVector();
376 P_minus->finalizeAssembly();
377 P_minus->setLocalAccessibleVector();
378 auto Q_plus = newZeroedInstance<GlobalVector>(matrix_specification);
379 auto Q_minus = newZeroedInstance<GlobalVector>(matrix_specification);
380 Q_plus->setLocalAccessibleVector();
381 Q_minus->setLocalAccessibleVector();
382
383 for (PetscInt i = start; i < end; ++i)
384 {
385 PetscInt cols;
386 PetscInt const* col_indices;
387
388 PetscScalar Q_plus_i = 0.0;
389 PetscScalar Q_minus_i = 0.0;
390
391 PetscCallAbort(PETSC_COMM_WORLD,
392 MatGetRow(F, i, &cols, &col_indices, nullptr));
393 auto const column_indices =
394 std::span<PetscInt const>(col_indices, cols);
395 for (auto const j : column_indices)
396 {
397 // skip diagonal
398 if (i == j)
399 {
400 continue;
401 }
402 PetscScalar const x_prev_i = x_prev[process_id]->get(i);
403 PetscScalar const x_prev_j = x_prev[process_id]->get(j);
404
405 Q_plus_i = std::max({Q_plus_i, 0., x_prev_j - x_prev_i});
406 Q_minus_i = std::min({Q_minus_i, 0., x_prev_j - x_prev_i});
407 }
408 Q_plus->set(i, Q_plus_i);
409 Q_minus->set(i, Q_minus_i);
410 PetscCallAbort(PETSC_COMM_WORLD,
411 MatRestoreRow(F, i, &cols, &col_indices, nullptr));
412 }
413 Q_plus->finalizeAssembly();
414 Q_plus->setLocalAccessibleVector();
415 Q_minus->finalizeAssembly();
416 Q_minus->setLocalAccessibleVector();
417
418 auto R_plus = newZeroedInstance<GlobalVector>(matrix_specification);
419 R_plus->setLocalAccessibleVector();
420
421 auto row_sums_M = newZeroedInstance<GlobalVector>(matrix_specification);
422 Vec M_L = row_sums_M->getRawVector();
423 PetscCallAbort(PETSC_COMM_WORLD, MatGetRowSum(M_raw, M_L));
424 row_sums_M->finalizeAssembly();
425 row_sums_M->setLocalAccessibleVector();
426
427 for (auto k = R_plus->getRangeBegin(); k < R_plus->getRangeEnd(); ++k)
428 {
429 PetscScalar const P_plus_i = P_plus->get(k);
430 if (P_plus_i == 0.0)
431 {
432 continue;
433 }
434 PetscScalar const mi = row_sums_M->get(k);
435 PetscScalar const Q_plus_i = Q_plus->get(k);
436 R_plus->set(k, std::min(1.0, mi * Q_plus_i / dt / P_plus_i));
437 }
438 R_plus->finalizeAssembly();
439 R_plus->setLocalAccessibleVector();
440 auto R_minus = newZeroedInstance<GlobalVector>(matrix_specification);
441 R_minus->setLocalAccessibleVector();
442
443 for (auto k = R_minus->getRangeBegin(); k < R_minus->getRangeEnd(); ++k)
444 {
445 PetscScalar const P_minus_i = P_minus->get(k);
446 if (P_minus_i == 0.0)
447 {
448 continue;
449 }
450 PetscScalar const mi = row_sums_M->get(k);
451 PetscScalar const Q_minus_i = Q_minus->get(k);
452 R_minus->set(k, std::min(1.0, mi * Q_minus_i / dt / P_minus_i));
453 }
454 R_minus->finalizeAssembly();
455 R_minus->setLocalAccessibleVector();
456 // walk over F, R_plus, and R_minus and compute alpha values; set entries in
457 // the rhs b that limit the antidiffusive flux
458 for (PetscInt i = start; i < end; ++i)
459 {
460 PetscInt cols;
461 PetscInt const* col_indices;
462 PetscScalar const* values;
463
464 PetscCallAbort(PETSC_COMM_WORLD,
465 MatGetRow(F, i, &cols, &col_indices, &values));
466
467 auto col_indices_span = std::span<PetscInt const>(col_indices, cols);
468 auto values_span = std::span<PetscScalar const>(values, cols);
469
470 double alpha_ij;
471
472 for (auto [j, fij] : std::views::zip(col_indices_span, values_span))
473 {
474 if (i == j)
475 {
476 continue;
477 }
478 if (fij > 0.)
479 {
480 alpha_ij = std::min(R_plus->get(i), R_minus->get(j));
481 }
482 else
483 {
484 alpha_ij = std::min(R_minus->get(i), R_plus->get(j));
485 }
486 b.add(i, alpha_ij * fij);
487 }
488 PetscCallAbort(PETSC_COMM_WORLD,
489 MatRestoreRow(F, i, &cols, &col_indices, &values));
490 }
491
492 // compute low-order operator
493 // update K_raw with D: K_raw += D;
494 for (PetscInt i = start; i < end; ++i)
495 {
496 PetscInt cols_D;
497 PetscInt const* col_indices_D;
498 PetscScalar const* values_D;
499 PetscCallAbort(PETSC_COMM_WORLD,
500 MatGetRow(D, i, &cols_D, &col_indices_D, &values_D));
501
502 auto const column_indices =
503 std::span<PetscInt const>(col_indices_D, cols_D);
504 auto const values = std::span<PetscScalar const>(values_D, cols_D);
505 for (auto [j, Dij] : std::views::zip(column_indices, values))
506 {
507 PetscCallAbort(PETSC_COMM_WORLD,
508 MatSetValues(K_raw, 1, &i, 1, &j, &Dij, ADD_VALUES));
509 }
510 PetscCallAbort(PETSC_COMM_WORLD,
511 MatRestoreRow(D, i, &cols_D, &col_indices_D, &values_D));
512 }
513 finalize(K_raw);
514 PetscCallAbort(PETSC_COMM_WORLD, MatScale(K_raw, -1.0));
515 K.finalizeAssembly();
516 M.setZero();
517 for (PetscInt i = start; i < end; ++i)
518 {
519 PetscScalar const row_sum_M_i = row_sums_M->get(i);
520 MatSetValues(M.getRawMatrix(), 1, &i, 1, &i, &row_sum_M_i,
521 INSERT_VALUES);
522 }
523 M.finalizeAssembly();
524 MatDestroy(&D);
525 MatDestroy(&F);
526}
527#endif // end of ifdef USE_PETSC
528
529void calculateFluxCorrectedTransport(
530 const double t, const double dt, std::vector<GlobalVector*> const& x,
531 std::vector<GlobalVector*> const& x_prev, int const process_id,
532 const MathLib::MatrixSpecifications& matrix_specification, GlobalMatrix& M,
533 GlobalMatrix& K, GlobalVector& b)
534{
535#ifndef USE_PETSC
536 calculateFluxCorrectedTransportSerial(t, dt, x, x_prev, process_id,
537 matrix_specification, M, K, b);
538#endif // end of ifndef USE_PETSC
539
540#ifdef USE_PETSC
541 calculateFluxCorrectedTransportPETSc(t, dt, x, x_prev, process_id,
542 matrix_specification, M, K, b);
543#endif // end of ifdef USE_PETSC
544}
545
546} // namespace detail
547
548void computeFluxCorrectedTransport(
549 NumericalStabilization const& stabilizer, const double t, const double dt,
550 std::vector<GlobalVector*> const& x,
551 std::vector<GlobalVector*> const& x_prev, int const process_id,
552 const MathLib::MatrixSpecifications& matrix_specification, GlobalMatrix& M,
553 GlobalMatrix& K, GlobalVector& b)
554{
555 // if needed, calling the Flux-Corrected-Transport function
556 return std::visit(
557 [&](auto&& stabilizer)
558 {
559 using Stabilizer = std::decay_t<decltype(stabilizer)>;
560 if constexpr (std::is_same_v<Stabilizer,
561 NumLib::FluxCorrectedTransport>)
562 {
563 return detail::calculateFluxCorrectedTransport(
564 t, dt, x, x_prev, process_id, matrix_specification, M, K,
565 b);
566 }
567 },
568 stabilizer);
569}
570} // namespace NumLib
OGS_FATAL
#define OGS_FATAL(...)
Definition Error.h:26
MathLib::EigenMatrix::setZero
void setZero()
reset data entries to zero.
Definition EigenMatrix.h:65
MathLib::EigenMatrix::getNumberOfRows
IndexType getNumberOfRows() const
return the number of rows
Definition EigenMatrix.h:53
MathLib::EigenMatrix::setValue
int setValue(IndexType row, IndexType col, double val)
Definition EigenMatrix.h:77
MathLib::EigenMatrix::get
double get(IndexType row, IndexType col) const
get value. This function returns zero if the element doesn't exist.
Definition EigenMatrix.h:131
MathLib::EigenMatrix::getRawMatrix
RawMatrixType & getRawMatrix()
Definition EigenMatrix.h:145
MathLib::EigenVector
Global vector based on Eigen vector.
Definition EigenVector.h:25
MathLib::EigenVector::add
void add(IndexType rowId, double v)
add entry
Definition EigenVector.h:76
NumLib::detail::calculateFluxCorrectedTransport
void calculateFluxCorrectedTransport(const double t, const double dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, const MathLib::MatrixSpecifications &matrix_specification, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b)
Definition FluxCorrectedTransport.cpp:529
NumLib::detail::newZeroedInstance
std::unique_ptr< MatrixVectorType > newZeroedInstance(MathLib::MatrixSpecifications const &matrix_specification)
Definition FluxCorrectedTransport.cpp:30
NumLib::detail::calculateFluxCorrectedTransportPETSc
void calculateFluxCorrectedTransportPETSc(const double t, const double dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, const MathLib::MatrixSpecifications &matrix_specification, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b)
Definition FluxCorrectedTransport.cpp:208
NumLib::computeFluxCorrectedTransport
void computeFluxCorrectedTransport(NumericalStabilization const &stabilizer, const double t, const double dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, const MathLib::MatrixSpecifications &matrix_specification, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b)
Definition FluxCorrectedTransport.cpp:548
NumLib::NumericalStabilization
std::variant< NoStabilization, IsotropicDiffusionStabilization, FullUpwind, FluxCorrectedTransport > NumericalStabilization
Definition CreateNumericalStabilization.h:31
ProcessLib::Graph::get
auto & get(Tuples &... ts)
Definition Get.h:59
ProcessLib::Graph::eval
auto eval(Function &f, Tuples &... ts) -> typename detail::GetFunctionReturnType< decltype(&Function::eval)>::type
Definition Apply.h:274
finalize
void finalize()
Definition ogs_python_module.cpp:192