Peano 4
Loading...
Searching...
No Matches
riemannSolverLinear.cpph
Go to the documentation of this file.
1
14#include <algorithm>
15
18
19namespace kernels {
20namespace aderdg {
21namespace generic {
22namespace c {
23
24template <bool useFlux, bool useNCP, bool useMM, typename SolverType>
25void riemannSolverLinear(SolverType& solver, double* FL, double* FR,
26 const double* const QL, const double* const QR,
27 const double t,
28 const double dt,
31 const int direction) {
32 /*
33 * For the linear kernels, we need the material parameters in the
34 * extrapolated predictor.
35 * We compute the averages of the material parameters but
36 * do not use them in the max eigenvalue calculation.
37 *
38 * Consider numberOfParameters in:
39 * QL,QR,QavL,QavR,Qdiff
40 *
41 * Not necessary to consider numberOfParameters in
42 * FL,FR,LL,LR
43 */
44
45 constexpr int numberOfVariables = SolverType::NumberOfVariables;
46 constexpr int numberOfParameters = SolverType::NumberOfParameters;
47 constexpr int numberOfData = numberOfVariables+numberOfParameters;
48 constexpr int basisSize = SolverType::Order+1;
49
50 // Compute the average variables and parameters from the left and the right
51 double QavL[numberOfData] = {0.0}; // size: numberOfVariables+numberOfParameters
52 double QavR[numberOfData] = {0.0}; // size: numberOfVariables+numberOfParameters
53 kernels::riemannsolvers::util::averageRiemannInputs<basisSize,numberOfData>(
55 kernels::riemannsolvers::util::averageRiemannInputs<basisSize,numberOfData>(
57
58 //get abs max eigenvalue
59 const double maxHyperbolicEigenvalueL = solver.maxEigenvalue(QavL, faceCentre, dx, t, dt, direction);
60 const double maxHyperbolicEigenvalueR = solver.maxEigenvalue(QavR, faceCentre, dx, t, dt, direction);
61 const double smax = std::max(maxHyperbolicEigenvalueL, maxHyperbolicEigenvalueR);
62
63 double Qdiff[numberOfData] = {0.0};
64 double Qavg[numberOfData] = {0.0};
65
66 idx2 idx_2d(Dimensions, numberOfVariables);
67
68 double flux[numberOfVariables];
69
70 double ncp[numberOfVariables] = {0.0};
71
72 for(int l = 0; l < numberOfData; l++) {
73 Qavg[l] = 0.5 * (QavR[l] + QavL[l]);
74 }
75
76 //We copy the averaged material parameters to Qdiff as they are used in the flux term
77 //These must not be overritten !
78 for(int l = numberOfVariables; l < numberOfData; l++) {
79 Qdiff[l] = Qavg[l];
80 }
81
82 {
83 idx2 idx_FLR(basisSize, numberOfVariables);
84 idx2 idx_QLR(basisSize, numberOfData);
85 std::fill_n (FL, basisSize * numberOfVariables, 0.0);
86 std::fill_n (FR, basisSize * numberOfVariables, 0.0);
87
88 for (int j = 0; j < basisSize; j++) {
89 for(int l = 0 ; l< numberOfVariables ; l++){
90 Qdiff[l] = 0.5 * (QR[idx_QLR(j, l)] - QL[idx_QLR(j, l)]);
91 }
92
93 if(useNCP) {
94
95 solver.nonconservativeProduct(
96 Qavg,
97 Qdiff,
98 faceCentre,
99 dx,
100 t,
101 dt,
102 direction,
103 ncp
104 );
105
106/* removed when the nonconservative product format was changed since it no longer fits with new form of ncp (MML)
107 if(useMM){
108 solver.multiplyMaterialParameterMatrix(Qavg, ncp);
109 }
110*/
111
112 for(int l=0; l < numberOfVariables; l++) {
113 FL[idx_FLR(j, l)] += ncp[l];
114 FR[idx_FLR(j, l)] += ncp[l];
115 }
116 }
117
118 if(useFlux){
119
120 solver.flux(
121 Qdiff,
122 faceCentre,
123 dx,
124 t,
125 dt,
126 direction,
127 flux
128 );
129
130// if(useMM){
131// solver.multiplyMaterialParameterMatrix(Qavg, flux);
132// }
133
134 for(int l=0; l < numberOfVariables; l++) {
135 FL[idx_FLR(j, l)] += flux[l];
136 FR[idx_FLR(j, l)] += flux[l];
137 }
138 }
139
140 for(int l=0; l < numberOfVariables; l++) {
141 FL[idx_FLR(j, l)] -= smax*Qdiff[l];
142 FR[idx_FLR(j, l)] += smax*Qdiff[l];
143 }
144 }
145 }
146
147}
148
149
150} // namespace c
151} // namespace generic
152} // namespace aderdg
153} // namespace kernels
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ j
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ so we are integrating with f$ phi_k phi_k dx
-with-ExaHyPE2-benchmark
void riemannSolverLinear(SolverType &solver, double *FL, double *FR, const double *const QL, const double *const QR, const double t, const double dt, const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &dx, const int direction)
This file is part of the ExaHyPE project.
Simple vector class.
Definition Vector.h:134