Peano
Loading...
Searching...
No Matches
riemannsolver_diffuse.h
Go to the documentation of this file.
1#include <algorithm>
2
3#include "kernels/GaussLegendreBasis.h"
4#include "kernels/KernelUtils.h"
5
6namespace Numerics {
7 template <class Shortcuts>
8 void rotate_dofs_inverse(double* Q, const int direction_1) {
9 Shortcuts s;
10
11 int direction_2 = direction_1 == 0 ? 1 : 0;
12 int direction_3 = direction_1 == 2 ? 1 : 2;
13
14 // velocities
15 double u = Q[s.v + 0];
16 double v = Q[s.v + 1];
17 double w = Q[s.v + 2];
18
19 // stresses
20 double sigma_xx = Q[s.sigma + 0];
21 double sigma_yy = Q[s.sigma + 1];
22 double sigma_zz = Q[s.sigma + 2];
23
24 double sigma_xy = Q[s.sigma + 3 + 0];
25 double sigma_xz = Q[s.sigma + 3 + 1];
26 double sigma_yz = Q[s.sigma + 3 + 2];
27
28 Q[s.v + direction_1] = u;
29 Q[s.v + direction_2] = v;
30 Q[s.v + direction_3] = w;
31 Q[s.sigma + direction_1] = sigma_xx;
32 Q[s.sigma + direction_2] = sigma_yy;
33 Q[s.sigma + direction_3] = sigma_zz;
34 Q[s.sigma + 2 + direction_1 + direction_2] = sigma_xy;
35 Q[s.sigma + 2 + direction_1 + direction_3] = sigma_xz;
36 Q[s.sigma + 2 + direction_2 + direction_3] = sigma_yz;
37 };
38
39 template <class Shortcuts>
40 void rotate_dofs(double* Q, const int direction_1) {
41 Shortcuts s;
42
43 int direction_2 = direction_1 == 0 ? 1 : 0;
44 int direction_3 = direction_1 == 2 ? 1 : 2;
45
46 // velocities
47 double u = Q[s.v + direction_1];
48 double v = Q[s.v + direction_2];
49 double w = Q[s.v + direction_3];
50
51 // stresses
52 double sigma_xx = Q[s.sigma + direction_1];
53 double sigma_yy = Q[s.sigma + direction_2];
54 double sigma_zz = Q[s.sigma + direction_3];
55
56 double sigma_xy = Q[s.sigma + 2 + direction_1 + direction_2];
57 double sigma_xz = Q[s.sigma + 2 + direction_1 + direction_3];
58 double sigma_yz = Q[s.sigma + 2 + direction_2 + direction_3];
59
60 Q[s.v + 0] = u;
61 Q[s.v + 1] = v;
62 Q[s.v + 2] = w;
63 Q[s.sigma + 0] = sigma_xx;
64 Q[s.sigma + 1] = sigma_yy;
65 Q[s.sigma + 2] = sigma_zz;
66 Q[s.sigma + 3 + 0] = sigma_xy;
67 Q[s.sigma + 3 + 1] = sigma_xz;
68 Q[s.sigma + 3 + 2] = sigma_yz;
69 };
70
71 template <class Shortcuts>
72 void right_eigenvectors(const double Q[], const int direction, double eigenvectors[]) {
73 Shortcuts s;
74 std::fill_n(eigenvectors, (s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables) * (s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables), 0.0);
75
76 kernels::idx2 idx2(s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables, s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables);
77 double rho = Q[s.rho];
78 double cp = Q[s.cp];
79 double cs = Q[s.cs];
80 double i_alpha = Q[s.alpha] <= 0 ? 0 : 1.0 / Q[s.alpha];
81
82 eigenvectors[idx2(0, s.sigma + 0)] = rho * (cp * cp);
83 eigenvectors[idx2(0, s.sigma + 1)] = rho * ((cp * cp) - 2 * (cs * cs));
84 eigenvectors[idx2(0, s.sigma + 2)] = rho * ((cp * cp) - 2 * (cs * cs));
85 eigenvectors[idx2(0, s.v + 0)] = cp;
86
87 eigenvectors[idx2(1, s.sigma + 3)] = rho * (cs * cs);
88 eigenvectors[idx2(1, s.v + 1)] = cs;
89
90 eigenvectors[idx2(2, s.sigma + 4)] = rho * (cs * cs);
91 eigenvectors[idx2(2, s.v + 2)] = cs;
92
93 eigenvectors[idx2(3, s.sigma + 1)] = 1;
94
95 eigenvectors[idx2(4, s.sigma + 2)] = 1;
96
97 eigenvectors[idx2(5, s.sigma + 5)] = 1;
98
99 eigenvectors[idx2(6, s.cs)] = 1;
100
101 eigenvectors[idx2(7, s.cp)] = 1;
102
103 eigenvectors[idx2(8, s.rho)] = 1;
104
105 eigenvectors[idx2(9, s.sigma + 0)] = -2.0 * Q[s.sigma + 0];
106 eigenvectors[idx2(9, s.sigma + 3)] = -2.0 * Q[s.sigma + 3];
107 eigenvectors[idx2(9, s.sigma + 4)] = -2.0 * Q[s.sigma + 4];
108 eigenvectors[idx2(9, s.v + 0)] = Q[s.v + 0] * i_alpha;
109 eigenvectors[idx2(9, s.v + 1)] = Q[s.v + 1] * i_alpha;
110 eigenvectors[idx2(9, s.v + 2)] = Q[s.v + 2] * i_alpha;
111 eigenvectors[idx2(9, s.alpha)] = 1;
112
113 eigenvectors[idx2(10, s.sigma + 4)] = rho * (cs * cs);
114 eigenvectors[idx2(10, s.v + 2)] = -cs;
115
116 eigenvectors[idx2(11, s.sigma + 3)] = rho * (cs * cs);
117 eigenvectors[idx2(11, s.v + 1)] = -cs;
118
119 eigenvectors[idx2(12, s.sigma + 0)] = rho * (cp * cp);
120 eigenvectors[idx2(12, s.sigma + 1)] = rho * ((cp * cp) - 2 * (cs * cs));
121 eigenvectors[idx2(12, s.sigma + 2)] = rho * ((cp * cp) - 2 * (cs * cs));
122 eigenvectors[idx2(12, s.v + 0)] = -cp;
123
124#ifdef Asserts
125 for (int i = 0; i < s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables; i++) {
126 for (int j = 0; j < s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables; j++) {
127 if (!std::isfinite(eigenvectors[idx2(i, j)])) {
128 for (int k = 0; k < s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables; k++) {
129 std::cout << Q[k] << std::endl;
130 }
131 }
132 assertion2(std::isfinite(eigenvectors[idx2(i, j)]), i, j);
133 }
134 }
135#endif
136 };
137
138 template <class Shortcuts>
139 void right_eigenvectors_inverse(const double Q[], const int direction, double eigenvectors[]) {
140 Shortcuts s;
141 std::fill_n(eigenvectors, (s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables) * (s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables), 0.0);
142 kernels::idx2 idx2(s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables, s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables);
143
144 double rho = Q[s.rho];
145 double cp = Q[s.cp];
146 double cs = Q[s.cs];
147 double i_alpha = Q[s.alpha] <= 0 ? 0 : 1.0 / Q[s.alpha];
148
149 eigenvectors[idx2(0, s.sigma + 0)] = 1.0 / (rho * (cp * cp) * 2.0);
150 eigenvectors[idx2(0, s.v + 0)] = 1.0 / (cp * 2.0);
151 eigenvectors[idx2(0, s.alpha)] = -(cp * rho * Q[s.v + 0] * i_alpha - 2 * Q[s.sigma + 0]) / (rho * (cp * cp) * 2.0);
152
153 eigenvectors[idx2(1, s.sigma + 3)] = 1.0 / (rho * (cs * cs) * 2.0);
154 eigenvectors[idx2(1, s.v + 1)] = 1.0 / (cs * 2.0);
155 eigenvectors[idx2(1, s.alpha)] = -(cs * rho * Q[s.v + 1] * i_alpha - 2 * Q[s.sigma + 3]) / (rho * (cs * cs) * 2.0);
156
157 eigenvectors[idx2(2, s.sigma + 4)] = 1.0 / (2.0 * rho * (cs * cs));
158 eigenvectors[idx2(2, s.v + 2)] = 1.0 / (cs * 2.0);
159 eigenvectors[idx2(2, s.alpha)] = -(cs * rho * Q[s.v + 2] * i_alpha - 2 * Q[s.sigma + 4]) / (rho * (cs * cs) * 2.0);
160
161 eigenvectors[idx2(3, s.sigma + 0)] = (-(cp * cp) + 2 * (cs * cs)) / (cp * cp);
162 eigenvectors[idx2(3, s.sigma + 1)] = 1.0;
163 eigenvectors[idx2(3, s.alpha)] = 2 * Q[s.sigma + 0] * ((-cp * cp) + 2 * (cs * cs)) / (cp * cp);
164
165 eigenvectors[idx2(4, s.sigma + 0)] = (-(cp * cp) + 2 * (cs * cs)) / (cp * cp);
166 eigenvectors[idx2(4, s.sigma + 2)] = 1.0;
167 eigenvectors[idx2(4, s.alpha)] = 2 * Q[s.sigma + 0] * ((-cp * cp) + 2 * (cs * cs)) / (cp * cp);
168
169 eigenvectors[idx2(5, s.sigma + 5)] = 1.0;
170
171 eigenvectors[idx2(6, s.cs)] = 1.0;
172
173 eigenvectors[idx2(7, s.cp)] = 1.0;
174
175 eigenvectors[idx2(8, s.rho)] = 1.0;
176
177 eigenvectors[idx2(9, s.alpha)] = 1.0;
178
179 eigenvectors[idx2(10, s.sigma + 4)] = 1.0 / ((cs * cs) * rho * 2.0);
180 eigenvectors[idx2(10, s.v + 2)] = -1.0 / (cs * 2.0);
181 eigenvectors[idx2(10, s.alpha)] = (cs * rho * Q[s.v + 2] * i_alpha + 2 * Q[s.sigma + 4]) / (rho * (cs * cs) * 2.0);
182
183 eigenvectors[idx2(11, s.sigma + 3)] = 1.0 / ((cs * cs) * rho * 2.0);
184 eigenvectors[idx2(11, s.v + 1)] = -1.0 / (cs * 2.0);
185 eigenvectors[idx2(11, s.alpha)] = (cs * rho * Q[s.v + 1] * i_alpha + 2 * Q[s.sigma + 3]) / (rho * (cs * cs) * 2.0);
186
187 eigenvectors[idx2(12, s.sigma + 0)] = 1.0 / ((cp * cp) * rho * 2.0);
188 eigenvectors[idx2(12, s.v + 0)] = -1.0 / (cp * 2.0);
189 eigenvectors[idx2(12, s.alpha)] = (cp * rho * Q[s.v + 0] * i_alpha + 2 * Q[s.sigma + 0]) / (rho * (cp * cp) * 2.0);
190
191#ifdef Asserts
192 for (int i = 0; i < s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables; i++) {
193 for (int j = 0; j < s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables; j++) {
194 if (!std::isfinite(eigenvectors[idx2(i, j)])) {
195 for (int k = 0; k < s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables; k++) {
196 std::cout << Q[k] << std::endl;
197 }
198 }
199 assertion2(std::isfinite(eigenvectors[idx2(i, j)]), i, j);
200 }
201 }
202#endif
203 };
204
205 template <class Shortcuts>
206 void eigenvalues(const double Q[], double lambda[]) {
207 Shortcuts s;
208 double cp = Q[s.cp];
209 double cs = Q[s.cs];
210 lambda[0] = -cp;
211 lambda[1] = -cs;
212 lambda[2] = -cs;
213 lambda[3] = 0;
214 lambda[4] = 0;
215 lambda[5] = 0;
216 lambda[6] = 0;
217 lambda[7] = 0;
218 lambda[8] = 0;
219 lambda[9] = 0;
220 lambda[10] = cs;
221 lambda[11] = cs;
222 lambda[12] = cp;
223 }
224
225 template <class Shortcuts>
226 double computeSmax(const double* const Q_L, const double* const Q_R) {
227 Shortcuts s;
228
229 // eigenvalues left and right
230 double Ev_R[(s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables)];
231 double Ev_L[(s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables)];
232 std::fill_n(Ev_R, (s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables), 0.0);
233 std::fill_n(Ev_L, (s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables), 0.0);
234
235 eigenvalues<Shortcuts>(Q_R, Ev_R);
236 eigenvalues<Shortcuts>(Q_L, Ev_L);
237
238 double ev_max = std::numeric_limits<double>::min();
239
240 for (int i = 0; i < s.NumberOfUnknowns; i++) {
241 ev_max = std::max(std::abs(Ev_R[i]), ev_max);
242 ev_max = std::max(std::abs(Ev_L[i]), ev_max);
243 }
244 return ev_max;
245 }
246
250 template <class Shortcuts>
251 void computeAbsA(const double* const Q_L, const double* const Q_R, double absA[], int direction) {
252 Shortcuts s;
253 kernels::idx2 idx2(Dimensions, s.NumberOfUnknowns);
254 kernels::idx2 idxA(s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables, s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables);
255 std::fill_n(absA, (s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables) * (s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables), 0.0);
256
257 // transposed eigenvectors and inverse eigenvectors
258 double R[(s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables) * (s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables)];
259 std::fill_n(R, (s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables) * (s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables), 0.0);
260 double R_inv[(s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables) * (s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables)];
261 std::fill_n(R_inv, (s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables) * (s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables), 0.0);
262
263 // Average Q
264 double Qavg[s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables];
265 std::fill_n(Qavg, (s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables), 0.0);
266
267 // Eigenvalues for the average
268 double Ev[(s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables)];
269 std::fill_n(Ev, (s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables), 0.0);
270
271 double ev_max = computeSmax<Shortcuts>(Q_L, Q_R);
272
273 // Compute average of both sides
274 for (int k = 0; k < s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables; k++) {
275 Qavg[k] = (Q_L[k] + Q_R[k]) * 0.5;
276 }
277
278 rotate_dofs<Shortcuts>(Qavg, direction);
279
280 assertion(!tarch::la::equals(Qavg[s.alpha], 0.0));
281 assertion(!tarch::la::equals(Qavg[s.rho], 0.0));
282
283 // compute eigenvectors for transformed third state
284 right_eigenvectors<Shortcuts>(Qavg, direction, R);
285 right_eigenvectors_inverse<Shortcuts>(Qavg, direction, R_inv);
286
287 // We store T R = R^T T^-1 in R (nice it's already transposed there ;))
288 for (int i = 0; i < s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables; i++) {
289 rotate_dofs_inverse<Shortcuts>(R + idxA(i, 0), direction);
290 }
291
292 // we store R^-1*T^-1 in R^-1
293 for (int i = 0; i < s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables; i++) {
294 rotate_dofs_inverse<Shortcuts>(R_inv + idxA(i, 0), direction);
295 }
296
297 eigenvalues<Shortcuts>(Qavg, Ev);
298
299 // Compute HLLEM anti diffusion
300 for (int i = 0; i < s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables; i++) {
301 Ev[i] = 1.0 - std::abs(Ev[i]) / ev_max; // we store deltaLambda directly in EV
302 }
303
304 // Add Rusanov diffusion
305 for (int i = 0; i < s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables; i++) {
306 absA[idxA(i, i)] = -0.5 * ev_max;
307 }
308 // std::cout << ev_max << std::endl;
309 // std::cout <<"absA" << std::endl;
310 /* for(int i=0; i< s.NumberOfAuxiliaryVariables+s.NumberOfUnknowns; i++){ */
311 /* for(int j=0; j< s.NumberOfAuxiliaryVariables+s.NumberOfUnknowns; j++){ */
312 /* //std::cout << absA[idxA(i,j)] << ";"; */
313 /* } */
314 /* //std::cout << std::endl; */
315 /* } */
316
317 for (int i = 0; i < s.NumberOfAuxiliaryVariables + s.NumberOfUnknowns; i++) {
318 for (int j = 0; j < s.NumberOfAuxiliaryVariables + s.NumberOfUnknowns; j++) {
319 R[idxA(j, i)] = R[idxA(j, i)] * Ev[j];
320 }
321 }
322
323 for (int i = 0; i < s.NumberOfAuxiliaryVariables + s.NumberOfUnknowns; i++) {
324 for (int j = 0; j < s.NumberOfAuxiliaryVariables + s.NumberOfUnknowns; j++) {
325 // for(int k=0; k< s.NumberOfAuxiliaryVariables+s.NumberOfUnknowns; k++){
326 for (int k = 6; k < 10; k++) { // This needs clarification: Maurizio only uses EV 7-10 and not all of them
327 absA[idxA(i, j)] += 0.5 * ev_max * R[idxA(k, i)] * R_inv[idxA(k, j)];
328 }
329 }
330 }
331
332 /* //std::cout <<"absA" << std::endl; */
333 /* for(int i=0; i< s.NumberOfAuxiliaryVariables+s.NumberOfUnknowns; i++){ */
334 /* for(int j=0; j< s.NumberOfAuxiliaryVariables+s.NumberOfUnknowns; j++){ */
335 /* //std::cout << absA[idxA(i,j)] << ";"; */
336 /* } */
337 /* //std::cout << std::endl; */
338 /* } */
339 }
340
341 template <class Solver, class Shortcuts>
342 void hllem(
343 Solver* solver, double* FL_, double* FR_, const double* const QL_, const double* const QR_, const int direction
344 ) {
345 Shortcuts s;
346
347 // approximate graidient
348 double gradQ[(s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables) * (Dimensions)];
349 std::fill_n(gradQ, ((s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables) * Dimensions), 0.0);
350
351 // average state
352 double ncp[s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables];
353 std::fill_n(ncp, (s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables), 0.0);
354
355 // average state
356 double Qavg[s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables];
357 std::fill_n(Qavg, (s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables), 0.0);
358
359 kernels::idx2 idx2(Dimensions, s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables);
360 kernels::idx2 idxA(s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables, s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables);
361
362 // HLLEM diffusion matrices;
363 double absA[(s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables) * (s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables)];
364 std::fill_n(absA, (s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables) * (s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables), 0.0);
365
366 computeAbsA<Shortcuts>(QL_, QR_, absA, direction);
367
368 // approximate gradient on node
369 for (int k = 0; k < s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables; k++) {
370 gradQ[idx2(direction, k)] = (QR_[k] - QL_[k]);
371 }
372
373 // compute average on node
374 for (int k = 0; k < s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables; k++) {
375 Qavg[k] = (QL_[k] + QR_[k]) * 0.5;
376 }
377 assertion2(!tarch::la::equals(Qavg[s.rho], 0.0), QL_[s.rho], QR_[s.rho]);
378
379 // compute for middle state ncp
380 solver->nonConservativeProduct(Qavg, gradQ, ncp);
381
382 /* std::cout <<"ncp" << std::endl; */
383 /* for(int j=0; j< s.NumberOfAuxiliaryVariables+s.NumberOfUnknowns; j++){ */
384 /* std::cout << ncp[j] << ";"; */
385 /* } */
386 /* std::cout << std::endl; */
387
388 /* std::cout <<"Qavg" << std::endl; */
389 /* for(int j=0; j< s.NumberOfAuxiliaryVariables+s.NumberOfUnknowns; j++){ */
390 /* std::cout << Qavg[j] << ";"; */
391 /* } */
392 /* std::cout << std::endl; */
393
394 /* std::cout <<"gradQ" << std::endl; */
395 /* for(int j=0; j< s.NumberOfAuxiliaryVariables+s.NumberOfUnknowns; j++){ */
396 /* std::cout << gradQ[idx2(direction,j)] << ";"; */
397 /* } */
398 /* std::cout << std::endl; */
399
400 for (int k = 0; k < s.NumberOfUnknowns; k++) {
401 for (int l = 0; l < s.NumberOfUnknowns + s.NumberOfAuxiliaryVariables; l++) {
402 FR_[k] += absA[idxA(k, l)] * gradQ[idx2(direction, l)];
403 }
404 }
405
406 // std::cout <<"FR" << std::endl;
407 /* for(int j=0; j< s.NumberOfUnknowns; j++){ */
408 /* std::cout << FR_[j] << ";"; */
409 /* } */
410 /* std::cout << std::endl; */
411
412 std::copy_n(FR_, s.NumberOfUnknowns, FL_);
413
414 for (int k = 0; k < s.NumberOfUnknowns; k++) {
415 FR_[k] -= 0.5 * ncp[k];
416 FL_[k] += 0.5 * ncp[k];
417 }
418 }
419}; // namespace Numerics
assertion2
#define assertion2(expr, param0, param1)
Definition Assertions.h:1019
assertion
#define assertion(expr)
Definition Assertions.h:1017
Numerics::rotate_dofs_inverse
void rotate_dofs_inverse(double *Q, const int direction_1)
Definition riemannsolver_diffuse.h:8
Numerics::right_eigenvectors_inverse
void right_eigenvectors_inverse(const double Q[], const int direction, double eigenvectors[])
Definition riemannsolver_diffuse.h:139
Numerics::right_eigenvectors
void right_eigenvectors(const double Q[], const int direction, double eigenvectors[])
Definition riemannsolver_diffuse.h:72
Numerics::computeAbsA
void computeAbsA(const double *const Q_L, const double *const Q_R, double absA[], int direction)
Computes AT = 0.5 * smax * T * (I + R*deltaLambda + R^-1)
Definition riemannsolver_diffuse.h:251
Numerics::rotate_dofs
void rotate_dofs(double *Q, const int direction_1)
Definition riemannsolver_diffuse.h:40
Numerics::eigenvalues
void eigenvalues(const double Q[], double lambda[])
Definition riemannsolver_diffuse.h:206
Numerics::computeSmax
double computeSmax(const double *const Q_L, const double *const Q_R)
Definition riemannsolver_diffuse.h:226
Numerics::hllem
void hllem(Solver *solver, double *FL_, double *FR_, const double *const QL_, const double *const QR_, const int direction)
Definition riemannsolver_diffuse.h:342
tarch::la::equals
bool equals(const Matrix< Rows, Cols, Scalar > &lhs, const Matrix< Rows, Cols, Scalar > &rhs, const Scalar &tolerance=NUMERICAL_ZERO_DIFFERENCE)
Compares to matrices on equality by means of a numerical accuracy.