Peano
Loading...
Searching...
No Matches
FWave.py
Go to the documentation of this file.
1# This file is part of the ExaHyPE2 project. For conditions of distribution and
2# use, please see the copyright notice at www.peano-framework.org
3
4fWave = r"""
5 /*******************************
6 * Flux and source preparation *
7 ******************************/
8
9 double fluxRight[NumberOfUnknowns] {0.0};
10 double fluxLeft[NumberOfUnknowns] {0.0};
11 double sourceRight[NumberOfUnknowns] {0.0};
12 double sourceLeft[NumberOfUnknowns] {0.0};
13
14 flux(QL, x, h, t, dt, normal, fluxLeft);
15 flux(QR, x, h, t, dt, normal, fluxRight);
16
17 if (normal == 1) {
18 sourceLeft[Shortcuts::rhoU + 1] = -GRAVITY * QL[Shortcuts::rho];
19 sourceLeft[Shortcuts::rhoE] = -GRAVITY * QL[Shortcuts::rhoU + 1];
20
21 sourceRight[Shortcuts::rhoU + 1] = -GRAVITY * QR[Shortcuts::rho];
22 sourceRight[Shortcuts::rhoE] = -GRAVITY * QR[Shortcuts::rhoU + 1];
23 }
24
25 std::fill_n(FL, NumberOfUnknowns, 0.0);
26 std::fill_n(FR, NumberOfUnknowns, 0.0);
27
28 /********************
29 * State extraction *
30 *******************/
31
32 // left state
33 const auto rhoL = QL[Shortcuts::rho];
34 const auto sqrtRhoL = sqrt(rhoL);
35 const auto uL = QL[Shortcuts::rhoU + 0] / rhoL;
36 const auto vL = QL[Shortcuts::rhoU + 1] / rhoL;
37 #if DIMENSIONS == 3
38 const auto wL = QL[Shortcuts::rhoU + 2] / rhoL;
39 const auto uSqL = uL * uL + vL * vL + wL * wL;
40 #else
41 const auto uSqL = uL * uL + vL * vL;
42 #endif
43 const auto EL = QL[Shortcuts::rhoE];
44 const auto pL = (GAMMA - 1.0) * (EL - 0.5 * rhoL * uSqL);
45
46 #if defined(GPUOffloadingOff)
47 assertion1(rhoL > 0, rhoL);
48 assertion1(EL > 0, EL);
49 assertion1(pL > 0, pL);
50 #endif
51
52 // right state
53 const auto rhoR = QR[Shortcuts::rho];
54 const auto sqrtRhoR = sqrt(rhoR);
55 const auto uR = QR[Shortcuts::rhoU + 0] / rhoR;
56 const auto vR = QR[Shortcuts::rhoU + 1] / rhoR;
57 #if DIMENSIONS == 3
58 const auto wR = QR[Shortcuts::rhoU + 2] / rhoR;
59 const auto uSqR = uR * uR + vR * vR + wR * wR;
60 #else
61 const auto uSqR = uR * uR + vR * vR;
62 #endif
63 const auto ER = QR[Shortcuts::rhoE];
64 const auto pR = (GAMMA - 1.0) * (ER - 0.5 * rhoR * uSqR);
65
66 #if defined(GPUOffloadingOff)
67 assertion1(rhoR > 0, rhoR);
68 assertion1(ER > 0, ER);
69 assertion1(pR > 0, pR);
70 #endif
71
72 /********************************************
73 * Roe averages from LeVeque Chapter 15.3.4 *
74 *******************************************/
75
76 const auto iSumSqrtRho = 1.0 / (sqrtRhoL + sqrtRhoR);
77
78 const auto gammaAvg = GAMMA;
79 const auto HRoe = (
80 (pR + ER) / sqrtRhoR
81 +
82 (pL + EL) / sqrtRhoL
83 ) * iSumSqrtRho;
84
85 tarch::la::Vector<DIMENSIONS, double> roeAveragedVelocities(0.0);
86 roeAveragedVelocities(0) = (sqrtRhoL * uL + sqrtRhoR * uR) * iSumSqrtRho;
87 roeAveragedVelocities(1) = (sqrtRhoL * vL + sqrtRhoR * vR) * iSumSqrtRho;
88 #if DIMENSIONS == 3
89 roeAveragedVelocities(2) = (sqrtRhoL * wL + sqrtRhoR * wR) * iSumSqrtRho;
90 #endif
91 const auto qSquared = 0.5 * (roeAveragedVelocities * roeAveragedVelocities);
92 const auto speedOfSoundRoe = sqrt((gammaAvg - 1) * (HRoe - qSquared));
93
94 /*************************************************************
95 * Roe Eigenvalues from LeVeque Chapter 18.8, Equation 18.47 *
96 ************************************************************/
97
98 // in x-direction, these are uRoe - speedOfSoundRoe, DIMENSIONS x uRoe, and uRoe + speedOfSoundRoe
99 tarch::la::Vector<NumberOfUnknowns, double> roeEigenvalues(0.0);
100
101 // first eigenvalue is left uRoe - speedOfSoundRoe
102 roeEigenvalues(0) = roeAveragedVelocities(normal) - speedOfSoundRoe;
103
104 // DIMENSIONS x uRoe
105 for (int eigenvalue = 1; eigenvalue <= DIMENSIONS; eigenvalue++) {
106 roeEigenvalues(eigenvalue) = roeAveragedVelocities(normal);
107 }
108
109 // last eigenvalue is right uRoe + speedOfSoundRoe
110 roeEigenvalues(DIMENSIONS + 1) = roeAveragedVelocities(normal) + speedOfSoundRoe;
111
112 /***********************************************
113 * Delta values from Bale et al. Equation 1.17 *
114 **********************************************/
115
116 tarch::la::Vector<NumberOfUnknowns, double> delta(0.0);
117
118 for (int unknown = 0; unknown < NumberOfUnknowns; ++unknown) {
119 delta(unknown) = fluxRight[unknown] - fluxLeft[unknown];
120 delta(unknown) -= 0.5 * h(normal) * (sourceLeft[unknown] + sourceRight[unknown]);
121 }
122
123 /****************************************************************
124 * Right Eigenvectors from LeVeque Chapter 18.8, Equation 18.48 *
125 ***************************************************************/
126
127 // @todo figure out for 3D
128 // => github.com/clawpack/riemann/blob/master/src/rpn3_euler.f90
129
130 const int transverse = normal == 0 ? 1 : 0;
131
132 // first right Eigenvector
133 tarch::la::Vector<NumberOfUnknowns, double> r1(0.0);
134 r1(0) = 1.0;
135 r1(1 + normal) = roeAveragedVelocities(normal) - speedOfSoundRoe;
136 r1(1 + transverse) = roeAveragedVelocities(transverse);
137 r1(3) = HRoe - roeAveragedVelocities(normal) * speedOfSoundRoe;
138
139 // second right Eigenvector
140 tarch::la::Vector<NumberOfUnknowns, double> r2(0.0);
141 r2(0) = 0.0;
142 r2(1 + normal) = 0.0;
143 r2(1 + transverse) = 1.0;
144 r2(3) = roeAveragedVelocities(transverse);
145
146 // third right Eigenvector
147 tarch::la::Vector<NumberOfUnknowns, double> r3(0.0);
148 r3(0) = 1.0;
149 r3(1) = roeAveragedVelocities(0);
150 r3(2) = roeAveragedVelocities(1);
151 r3(3) = qSquared;
152
153 // fourth right Eigenvector
154 tarch::la::Vector<NumberOfUnknowns, double> r4(0.0);
155 r4(0) = 1.0;
156 r4(1 + normal) = roeAveragedVelocities(normal) + speedOfSoundRoe;
157 r4(1 + transverse) = roeAveragedVelocities(transverse);
158 r4(3) = HRoe + roeAveragedVelocities(normal) * speedOfSoundRoe;
159
160 /***************************************************
161 * Beta values following Bale et al. Equation 1.14 *
162 **************************************************/
163 // inverted Eigenvectors taken from https://github.com/clawpack/riemann/blob/master/src/rpn2_euler_4wave.f90
164
165 // @todo figure out for 3D
166 tarch::la::Vector<NumberOfUnknowns, double> beta(0.0);
167
168 beta(2) = (
169 (HRoe - 2 * qSquared) * delta(0)
170 + roeAveragedVelocities(normal) * delta(1 + normal)
171 + roeAveragedVelocities(transverse) * delta(1 + transverse)
172 - delta(3)
173 ) / (HRoe - qSquared);
174
175 beta(1) = delta(1 + transverse) - roeAveragedVelocities(transverse) * delta(0);
176
177 beta(3) = (
178 delta(1 + normal)
179 + (speedOfSoundRoe - roeAveragedVelocities(normal)) * delta(0)
180 - speedOfSoundRoe * beta(2)
181 ) / (2 * speedOfSoundRoe);
182
183 beta(0) = delta(0) - beta(2) - beta(3);
184
185 /****************
186 * fWave fluxes *
187 ***************/
188
189 const tarch::la::Vector<NumberOfUnknowns, double> fWaves[NumberOfUnknowns] = {
190 beta(0) * r1,
191 beta(1) * r2,
192 beta(2) * r3,
193 beta(3) * r4
194 };
195
196 for (int wave = 0; wave < NumberOfUnknowns; ++wave) {
197 const auto eigenvalue = roeEigenvalues(wave);
198 if (eigenvalue < 0.0) {
199 // left going wave
200 for (int unknown = 0; unknown < NumberOfUnknowns; ++unknown) {
201 FL[unknown] += fWaves[wave](unknown);
202 }
203 } else {
204 // right going wave
205 for (int unknown = 0; unknown < NumberOfUnknowns; ++unknown) {
206 FR[unknown] -= fWaves[wave](unknown);
207 }
208 }
209 }
210
211 return tarch::la::maxAbs(roeEigenvalues);
212"""