Peano
Loading...
Searching...
No Matches
MGCCZ4Kernels.cpp
Go to the documentation of this file.
1
2#include "tarch/multicore/multicore.h"
3#include "MGCCZ4Kernels.h"
4
5#include "Constants.h"
6
7#include <limits>
8#include <cmath>
9#include <stdio.h>
10#include <string.h>
11
12
13#if defined(GPUOffloadingOMP)
14#pragma omp declare target
15#endif
16void examples::exahype2::mgccz4::enforceMGCCZ4constraints(double * luh)
17{
18 double g_cov[3][3] = { {luh[0], luh[1], luh[2]}, {luh[1], luh[3], luh[4]}, {luh[2], luh[4], luh[5]} };
19 const double det = luh[0]*luh[3]*luh[5] - luh[0]*luh[4]*luh[4] - luh[1]*luh[1]*luh[5] + 2*luh[1]*luh[2]*luh[4] -luh[2]*luh[2]*luh[3];
20
21 const double phisq = 1./std::cbrt(det);
22 for (int i = 0; i < 3; i++)
23 for (int j = 0; j < 3; j++) g_cov[i][j] *= phisq;
24
25 const double det2 = g_cov[0][0]*g_cov[1][1]*g_cov[2][2] -
26 g_cov[0][0]*g_cov[1][2]*g_cov[2][1] - g_cov[1][0]*g_cov[0][1]*g_cov[2][2] +
27 g_cov[1][0]*g_cov[0][2]*g_cov[2][1] + g_cov[2][0]*g_cov[0][1]*g_cov[1][2] -
28 g_cov[2][0]*g_cov[0][2]*g_cov[1][1];
29
30
31 const double invdet = 1./det2;
32 const double g_contr[3][3] = {
33 { ( g_cov[1][1]*g_cov[2][2]-g_cov[1][2]*g_cov[2][1])*invdet, -(g_cov[0][1]*g_cov[2][2]-g_cov[0][2]*g_cov[2][1])*invdet, -(-g_cov[0][1]*g_cov[1][2]+g_cov[0][2]*g_cov[1][1])*invdet},
34 {-( g_cov[1][0]*g_cov[2][2]-g_cov[1][2]*g_cov[2][0])*invdet, (g_cov[0][0]*g_cov[2][2]-g_cov[0][2]*g_cov[2][0])*invdet, -( g_cov[0][0]*g_cov[1][2]-g_cov[0][2]*g_cov[1][0])*invdet},
35 {-(-g_cov[1][0]*g_cov[2][1]+g_cov[1][1]*g_cov[2][0])*invdet, -(g_cov[0][0]*g_cov[2][1]-g_cov[0][1]*g_cov[2][0])*invdet, ( g_cov[0][0]*g_cov[1][1]-g_cov[0][1]*g_cov[1][0])*invdet}
36 };
37 double Aex[3][3] = { {luh[6], luh[7], luh[8]}, {luh[7], luh[9], luh[10]}, {luh[8], luh[10], luh[11]} };
38
39 double traceA = 0;
40 for (int i=0;i<3;i++)
41 for (int j=0;j<3;j++) traceA += g_contr[i][j]*Aex[i][j];
42
43 for (int i=0;i<3;i++)
44 for (int j=0;j<3;j++) Aex[i][j] -= 1./3. * traceA * g_cov[i][j];
45
46 luh[ 0] = g_cov[0][0];
47 luh[ 1] = g_cov[0][1];
48 luh[ 2] = g_cov[0][2];
49 luh[ 3] = g_cov[1][1];
50 luh[ 4] = g_cov[1][2];
51 luh[ 5] = g_cov[2][2];
52
53 luh[ 6] = Aex[0][0];
54 luh[ 7] = Aex[0][1];
55 luh[ 8] = Aex[0][2];
56 luh[ 9] = Aex[1][1];
57 luh[10] = Aex[1][2];
58 luh[11] = Aex[2][2];
59 //
60 // As suggested by our PRD referee, we also enforce the algebraic constraint that results from the first spatial derivative of the constraint
61 // det \tilde{\gamma}_ij = 0, which leads to
62 //
63 // \tilde{\gamma}^{ij} D_kij = 0
64 //
65 // and is thus a condition of trace-freeness on all submatrices D_kij for k=1,2,3.
66 //
67
68 double DD[3][3][3] = {
69 {{luh[35], luh[36], luh[37]}, {luh[36], luh[38], luh[39]}, {luh[37], luh[39], luh[40]}},
70 {{luh[41], luh[42], luh[43]}, {luh[42], luh[44], luh[45]}, {luh[43], luh[45], luh[46]}},
71 {{luh[47], luh[48], luh[49]}, {luh[48], luh[50], luh[51]}, {luh[49], luh[51], luh[52]}}
72 };
73
74 for (int l=0;l<3;l++)
75 {
76 double traceDk = 0;
77 for (int i=0;i<3;i++)
78 for (int j=0;j<3;j++) traceDk += g_contr[i][j]*DD[l][i][j];
79
80 for (int i=0;i<3;i++)
81 for (int j=0;j<3;j++) DD[l][i][j] -= 1./3 * g_cov[i][j]*traceDk;
82 }
83
84 luh[35] = DD[0][0][0];
85 luh[36] = DD[0][0][1];
86 luh[37] = DD[0][0][2];
87 luh[38] = DD[0][1][1];
88 luh[39] = DD[0][1][2];
89 luh[40] = DD[0][2][2];
90
91 luh[41] = DD[1][0][0];
92 luh[42] = DD[1][0][1];
93 luh[43] = DD[1][0][2];
94 luh[44] = DD[1][1][1];
95 luh[45] = DD[1][1][2];
96 luh[46] = DD[1][2][2];
97
98 luh[47] = DD[2][0][0];
99 luh[48] = DD[2][0][1];
100 luh[49] = DD[2][0][2];
101 luh[50] = DD[2][1][1];
102 luh[51] = DD[2][1][2];
103 luh[52] = DD[2][2][2];
104}
105#if defined(GPUOffloadingOMP)
106#pragma omp end declare target
107#endif
108
109#if defined(GPUOffloadingOMP)
110#pragma omp declare target
111#endif
112void examples::exahype2::mgccz4::source(double* S, const double* const Q,
113 const int MGCCZ4LapseType,
114 const double MGCCZ4ds,
115 const double MGCCZ4c,
116 const double MGCCZ4e,
117 const double MGCCZ4f,
118 const double MGCCZ4bs,
119 const double MGCCZ4sk,
120 const double MGCCZ4xi,
121 const double MGCCZ4itau,
122 const double MGCCZ4eta,
123 const double MGCCZ4k1,
124 const double MGCCZ4k2,
125 const double MGCCZ4k3
126 )
127{
128 const double alpha = std::exp(std::fmax(-20., std::fmin(20.,Q[16])));
129 //printf("alpha %f\n",alpha);
130 double fa = 1.0;
131 double faa = 0.0;
132 if (MGCCZ4LapseType==1)
133 {
134 fa = 2./alpha;
135 faa = -fa/alpha;
136 }
137
138 // Note g_cov is symmetric
139 const double g_cov[3][3] = { {Q[0], Q[1], Q[2]}, {Q[1], Q[3], Q[4]}, {Q[2], Q[4], Q[5]} };
140 const double det = Q[0]*Q[3]*Q[5] - Q[0]*Q[4]*Q[4] - Q[1]*Q[1]*Q[5] + 2*Q[1]*Q[2]*Q[4] -Q[2]*Q[2]*Q[3];
141 const double invdet = 1./det;
142
143 const double g_contr[3][3] = {
144 { ( Q[3]*Q[5]-Q[4]*Q[4])*invdet, -( Q[1]*Q[5]-Q[2]*Q[4])*invdet, -(-Q[1]*Q[4]+Q[2]*Q[3])*invdet},
145 {-( Q[1]*Q[5]-Q[4]*Q[2])*invdet, ( Q[0]*Q[5]-Q[2]*Q[2])*invdet, -( Q[0]*Q[4]-Q[2]*Q[1])*invdet},
146 {-(-Q[1]*Q[4]+Q[3]*Q[2])*invdet, -( Q[0]*Q[4]-Q[1]*Q[2])*invdet, ( Q[0]*Q[3]-Q[1]*Q[1])*invdet}
147 };
148
149 // NOTE Aex is symmetric
150 double Aex[3][3] = { {Q[6], Q[7], Q[8]}, {Q[7], Q[9], Q[10]}, {Q[8], Q[10], Q[11]} };
151
152 double traceA = 0;
153 for (int i=0;i<3;i++)
154 for (int j=0;j<3;j++) traceA += g_contr[i][j]*Aex[i][j];
155
156 for (int i=0;i<3;i++)
157 for (int j=0;j<3;j++) Aex[i][j] -= 1./3. * traceA * g_cov[i][j];
158
159 // Matrix multiplications Amix = matmul(g_contr, Aex) Aup = matmul(g_contr, transpose(Amix))
160 double Amix[3][3]={0,0,0,0,0,0,0,0,0};
161 for (int i = 0; i < 3; i++)
162 for (int j = 0; j < 3; j++)
163 for (int u = 0; u < 3; u++) Amix[i][j] += g_contr[i][u] * Aex[u][j];
164 double Aup[3][3]={0,0,0,0,0,0,0,0,0};
165 for (int i = 0; i < 3; i++)
166 for (int j = 0; j < 3; j++)
167 for (int u = 0; u < 3; u++) Aup[i][j] += g_contr[i][u] * Amix[j][u]; // Note the transposition is in the indices
168
169 const double Theta = Q[12]; // needed later
170
171 const double DD[3][3][3] = {
172 {{Q[35], Q[36], Q[37]}, {Q[36], Q[38], Q[39]}, {Q[37], Q[39], Q[40]}},
173 {{Q[41], Q[42], Q[43]}, {Q[42], Q[44], Q[45]}, {Q[43], Q[45], Q[46]}},
174 {{Q[47], Q[48], Q[49]}, {Q[48], Q[50], Q[51]}, {Q[49], Q[51], Q[52]}}
175 };
176
177 double dgup[3][3][3] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
178 for (int k = 0; k < 3; k++)
179 for (int m = 0; m < 3; m++)
180 for (int l = 0; l < 3; l++)
181 for (int n = 0; n < 3; n++)
182 for (int j = 0; j < 3; j++) dgup[k][m][l] -= 2.0*g_contr[m][n]*g_contr[j][l]*DD[k][n][j];
183
184 const double PP[3] = {Q[55], Q[56], Q[57]};
185
186 double Christoffel[3][3][3] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
187 double Christoffel_tilde[3][3][3] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
188
189 for (int j = 0; j < 3; j++)
190 for (int i = 0; i < 3; i++)
191 for (int k = 0; k < 3; k++)
192 for (int l = 0; l < 3; l++)
193 {
194 Christoffel_tilde[i][j][k] += g_contr[k][l] * ( DD[i][j][l] + DD[j][i][l] - DD[l][i][j] );
195 Christoffel[i][j][k] += g_contr[k][l] * ( DD[i][j][l] + DD[j][i][l] - DD[l][i][j] ) - g_contr[k][l] * ( g_cov[j][l] * PP[i] + g_cov[i][l] * PP[j] - g_cov[i][j] * PP[l] );
196 }
197
198 double Gtilde[3] = {0,0,0};
199 for (int l = 0; l < 3; l++)
200 for (int j = 0; j < 3; j++)
201 for (int i = 0; i < 3; i++) Gtilde[i] += g_contr[j][l] * Christoffel_tilde[j][l][i];
202
203 const double Ghat[3] = {Q[13], Q[14], Q[15]};
204
205 const double phi = std::exp(std::fmax(-20., std::fmin(20.,Q[54])));
206 const double phi2 = phi*phi;
207
208 double Z[3] = {0,0,0}; // Matrix vector multiplications
209 for (int i=0;i<3;i++)
210 for (int j=0;j<3;j++) Z[i] += 0.5*MGCCZ4ds*( g_cov[i][j]* (Ghat[j] - Gtilde[j]));
211 double Zup[3] = {0,0,0};
212 for (int i=0;i<3;i++)
213 for (int j=0;j<3;j++) Zup[i] += phi2 * g_contr[i][j] * Z[j];
214
215 double dChristoffelSrc[3][3][3][3] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
216 double dChristoffel_tildeSrc[3][3][3][3] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
217 for (int l = 0; l < 3; l++)
218 for (int m = 0; m < 3; m++)
219 for (int j = 0; j < 3; j++)
220 for (int i = 0; i < 3; i++)
221 for (int k = 0; k < 3; k++)
222 {
223 dChristoffelSrc[k][i][j][m] += dgup[k][m][l]*( DD[i][j][l]+ DD[j][i][l]-DD[l][i][j])
224 - dgup[k][m][l]*(g_cov[j][l]*PP[i]+g_cov[i][l]*PP[j]-g_cov[i][j]*PP[l])
225 -2.0*g_contr[m][l]*(DD[k][j][l]*PP[i]+DD[k][i][l]*PP[j]-DD[k][i][j]*PP[l]);
226
227 dChristoffel_tildeSrc[k][i][j][m] += dgup[k][m][l]*(DD[i][j][l]+DD[j][i][l]-DD[l][i][j]);
228 }
229
230
231 double RiemannSrc[3][3][3][3] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
232 for (int m = 0; m < 3; m++)
233 for (int j = 0; j < 3; j++)
234 for (int k = 0; k < 3; k++)
235 for (int i = 0; i < 3; i++)
236 {
237 RiemannSrc[i][k][j][m] = dChristoffelSrc[k][i][j][m] - dChristoffelSrc[j][i][k][m];
238 for (int l = 0; l < 3; l++) RiemannSrc[i][k][j][m] += Christoffel[i][j][l]*Christoffel[l][k][m] - Christoffel[i][k][l]*Christoffel[l][j][m];
239 }
240
241 double RicciSrc[3][3] = {0,0,0,0,0,0,0,0,0};
242 for (int l = 0; l < 3; l++)
243 for (int n = 0; n < 3; n++)
244 for (int m = 0; m < 3; m++) RicciSrc[m][n] += RiemannSrc[m][l][n][l];
245
246 // Here we directly compute the derivative of Gtilde from its original definition as contracted Christoffel symbol,
247 // without assuming unit determinant of the conformal metric. Back to the roots, and as few assumptions as possible...
248 double dGtildeSrc[3][3] = {0,0,0,0,0,0,0,0,0};
249 for (int l = 0; l < 3; l++)
250 for (int j = 0; j < 3; j++)
251 for (int i = 0; i < 3; i++)
252 for (int k = 0; k < 3; k++) dGtildeSrc[k][i] += dgup[k][j][l]*Christoffel_tilde[j][l][i] + g_contr[j][l]*dChristoffel_tildeSrc[k][j][l][i];
253
254 double dZSrc[3][3] = {0,0,0,0,0,0,0,0,0};
255 for (int j = 0; j < 3; j++)
256 for (int i = 0; i < 3; i++)
257 for (int k = 0; k < 3; k++) dZSrc[k][i] += MGCCZ4ds*(DD[k][i][j]*(Ghat[j]-Gtilde[j]) - 0.5*g_cov[i][j]*dGtildeSrc[k][j]);
258
259 double nablaZSrc[3][3] = {0,0,0,0,0,0,0,0,0};
260 for (int j = 0; j < 3; j++)
261 for (int i = 0; i < 3; i++)
262 {
263 nablaZSrc[i][j] = dZSrc[i][j];
264 for (int k = 0; k < 3; k++) nablaZSrc[i][j] -= Christoffel[i][j][k]*Z[k];
265 }
266
267 double RicciPlusNablaZSrc[3][3];
268 for (int i = 0; i < 3; i++)
269 for (int j = 0; j < 3; j++) RicciPlusNablaZSrc[i][j] = RicciSrc[i][j] + nablaZSrc[i][j] + nablaZSrc[j][i];
270
271 double RPlusTwoNablaZSrc = 0;
272 for (int i = 0; i < 3; i++)
273 for (int j = 0; j < 3; j++) RPlusTwoNablaZSrc += g_contr[i][j]*RicciPlusNablaZSrc[i][j]; // TODO fuse these steps
274 RPlusTwoNablaZSrc*=phi2;
275
276
277 const double AA[3] = {Q[23], Q[24], Q[25]};
278
279 double nablaijalphaSrc[3][3];
280 for (int j = 0; j < 3; j++)
281 for (int i = 0; i < 3; i++)
282 {
283 nablaijalphaSrc[i][j] = alpha*AA[j]*AA[i];
284 for (int k = 0; k < 3; k++) nablaijalphaSrc[i][j] -= alpha*Christoffel[i][j][k]*AA[k];
285 }
286
287 double nablanablaalphaSrc=0;
288 for (int i = 0; i < 3; i++)
289 for (int j = 0; j < 3; j++) nablanablaalphaSrc += g_contr[i][j]*nablaijalphaSrc[i][j];
290 nablanablaalphaSrc *= phi2;
291
292
293 double SecondOrderTermsSrc[3][3];
294 for (int i = 0; i < 3; i++)
295 for (int j = 0; j < 3; j++) SecondOrderTermsSrc[i][j] = -nablaijalphaSrc[i][j] + alpha*RicciPlusNablaZSrc[i][j];
296
297 double traceSrc = 0;
298 for (int i = 0; i < 3; i++)
299 for (int j = 0; j < 3; j++) traceSrc += g_contr[i][j]*SecondOrderTermsSrc[i][j];
300
301 for (int i = 0; i < 3; i++)
302 for (int j = 0; j < 3; j++) SecondOrderTermsSrc[i][j] -= 1./3 * traceSrc * g_cov[i][j];
303
304 double dtgamma[3][3];
305 for (int i = 0; i < 3; i++)
306 for (int j = 0; j < 3; j++) dtgamma[i][j] = -2.0 * alpha * Aex[i][j] - MGCCZ4itau*(det -1.0)*g_cov[i][j];
307
308 const double BB[3][3] = {
309 {Q[26], Q[27], Q[28]}, {Q[29], Q[30], Q[31]}, {Q[32], Q[33], Q[34]}
310 };
311
312 double traceB = BB[0][0] + BB[1][1] + BB[2][2];
313 const double beta[3] = {Q[17], Q[18], Q[19]};
314
315 for (int k = 0; k < 3; k++)
316 for (int j = 0; j < 3; j++)
317 for (int i = 0; i < 3; i++) dtgamma[i][j] += g_cov[k][i]*BB[j][k] + g_cov[k][j]*BB[i][k] - 2./3. *g_cov[i][j]*BB[k][k] + 2.*beta[k]*DD[k][i][j];
318
319
320 // MATMUL(Aex,Amix)
321 double Atemp[3][3]={0,0,0,0,0,0,0,0,0};
322 for (int i = 0; i < 3; i++)
323 for (int j = 0; j < 3; j++)
324 for (int u = 0; u < 3; u++) Atemp[i][j] += Aex[i][u] * Amix[u][j];
325
326 const double traceK = Q[53];
327
329 double dtK[3][3];
330 for (int i = 0; i < 3; i++)
331 for (int j = 0; j < 3; j++) dtK[i][j] = phi2*SecondOrderTermsSrc[i][j] + alpha*Aex[i][j]*(traceK-2.*Theta) - 2.*alpha*Atemp[i][j] - MGCCZ4itau*g_cov[i][j]*traceA;
332
333 for (int j = 0; j < 3; j++)
334 for (int i = 0; i < 3; i++)
335 for (int k = 0; k < 3; k++) dtK[i][j] += Aex[k][i]*BB[j][k] + Aex[k][j]*BB[i][k] - 2./3.*Aex[i][j]*BB[k][k];
336
337 const double K0 = Q[58];
338
339 double dtTraceK = -nablanablaalphaSrc + alpha*(RPlusTwoNablaZSrc + traceK*traceK - 2.0*MGCCZ4c*Theta*traceK) -3.0*alpha*MGCCZ4k1*(1.+MGCCZ4k2)*Theta;
340 double dtphi = beta[0]*PP[0] + beta[1]*PP[1] + beta[2]*PP[2] + 1./3.*(alpha*traceK-traceB);
341 double dtalpha = -alpha*fa*(traceK-K0-2.*MGCCZ4c*Theta)+beta[0]*AA[0]+beta[1]*AA[1]+beta[2]*AA[2];
342
343
344 double Aupdown = 0;
345 for (int i = 0; i < 3; i++)
346 for (int j = 0; j < 3; j++) Aupdown += Aex[i][j]*Aup[i][j];
347
348
349 double sumzupaa = 0.0;
350 for (int i = 0; i < 3; i++) sumzupaa += Zup[i]*AA[i];
351 const double dtTheta = 0.5*alpha*MGCCZ4e*MGCCZ4e*(RPlusTwoNablaZSrc - Aupdown + 2./3.*traceK*traceK) - alpha*(MGCCZ4c*Theta*traceK + sumzupaa + MGCCZ4k1*(2.+MGCCZ4k2)*Theta); // Baojiu
352
353
354 double dtGhat[3] = {0,0,0};
355 for (int i = 0; i < 3; i++)
356 {
357 double temp1=0, temp2=0, temp3=0, temp4=0, temp5=0, temp6=0;
358 for (int m = 0; m < 3; m++)
359 {
360 temp1 += Aup[i][m]*PP[m];
361 temp3 += Aup[i][m]*AA[m];
362 temp2 += g_contr[m][i]*(-Theta*AA[m] -2./3.*traceK*Z[m]);
363 temp4 += g_contr[i][m]*Z[m];
364 temp5 += Gtilde[m]*BB[m][i];
365 for (int n = 0; n < 3; n++) temp6 += Christoffel_tilde[m][n][i]*Aup[m][n];
366 }
367 dtGhat[i] += 2.*alpha*(temp6 - 3.*temp1 + temp2 - temp3 - MGCCZ4k1*temp4) - temp5 + 2./3.*Gtilde[i]*traceB;
368
369 for (int l = 0; l < 3; l++)
370 for (int k = 0; k < 3; k++) dtGhat[i] += 2.*MGCCZ4k3*(2./3.*g_contr[i][l]*Z[l]*BB[k][k] - g_contr[l][k]*Z[l]*BB[k][i]);
371 }
372
373 double ov[3];
374 for (int k = 0; k < 3; k++)
375 {
376 double temp=0;
377 for (int i = 0; i < 3; i++)
378 for (int j = 0; j < 3; j++) temp += dgup[k][i][j]*Aex[i][j];
379 ov[k] = 2*alpha*temp;
380 }
381
382 // matrix vector multiplication in a loop and add result to existing vector
383 for (int i = 0; i < 3; i++)
384 for (int j = 0; j < 3; j++) dtGhat[i] += MGCCZ4sk*g_contr[i][j]*ov[j];
385
386 const double myb[3] = {Q[20], Q[21], Q[22]};
387
388 double dtbb[3];
389 for (int i = 0; i < 3; i++) dtbb[i] = MGCCZ4sk*(MGCCZ4xi*dtGhat[i] - MGCCZ4eta*myb[i]);
390
391 double dtbeta[3];
392 for (int i = 0; i < 3; i++) dtbeta[i] = MGCCZ4f*myb[i];
393 for (int i = 0; i < 3; i++) dtbeta[i] += MGCCZ4bs*(beta[0]*BB[0][i] + beta[1]*BB[1][i] + beta[2]*BB[2][i]);
394 for (int i = 0; i < 3; i++) dtbeta[i] *= MGCCZ4sk;
395
396
397 // Auxiliary variables
398 double dtA[3] ={0,0,0};
399 for (int i = 0; i < 3; i++)
400 {
401 dtA[i] = -alpha*AA[i]*(fa+alpha*faa)*(traceK - K0 - 2.*MGCCZ4c*Theta);
402 for (int j = 0; j < 3; j++) dtA[i] += BB[i][j]*AA[j];
403 }
404
405 for (int k = 0; k < 3; k++)
406 {
407 double temp = 0;
408 for (int i = 0; i < 3; i++)
409 for (int j = 0; j < 3; j++) temp+= dgup[k][i][j]*Aex[i][j];
410 dtA[k] += -MGCCZ4sk*alpha*fa*temp;
411 }
412
413 double dtB[3][3] ={0,0,0,0,0,0,0,0,0};
414 for (int i = 0; i < 3; i++)
415 for (int j = 0; j < 3; j++)
416 for (int u = 0; u < 3; u++) dtB[i][j] += MGCCZ4sk*(BB[i][u] * BB[u][j]);
417
418 double dtD[3][3][3] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
419 for (int m = 0; m < 3; m++)
420 for (int j = 0; j < 3; j++)
421 for (int i = 0; i < 3; i++)
422 for (int k = 0; k < 3; k++)
423 for (int n = 0; n < 3; n++) dtD[k][i][j] += 1./3*alpha*g_cov[i][j]*dgup[k][n][m]*Aex[n][m]; // explicitly remove the trace of tilde A again
424
425 for (int j = 0; j < 3; j++)
426 for (int i = 0; i < 3; i++)
427 for (int k = 0; k < 3; k++)
428 {
429 dtD[k][i][j] -= alpha*AA[k]*Aex[i][j];
430 for (int l = 0; l < 3; l++) dtD[k][i][j] += BB[k][l]*DD[l][i][j] + BB[j][l]*DD[k][l][i] + BB[i][l]*DD[k][l][j] - 2./3.*BB[l][l]*DD[k][i][j];
431 }
432
433 double dtP[3] = {0,0,0};
434 for (int i = 0; i < 3; i++)
435 for (int j = 0; j < 3; j++) dtP[i] += BB[i][j]*PP[j];
436
437 for (int k = 0; k < 3; k++)
438 {
439 double temp=0;
440 for (int i = 0; i < 3; i++)
441 for (int j = 0; j < 3; j++) temp += dgup[k][i][j]*Aex[i][j];
442 dtP[k] += 1./3.*alpha*(AA[k]*traceK + MGCCZ4sk*temp);
443 }
444
445
446 S[0] = dtgamma[0][0];
447 S[1] = dtgamma[0][1];
448 S[2] = dtgamma[0][2];
449 S[3] = dtgamma[1][1];
450 S[4] = dtgamma[1][2];
451 S[5] = dtgamma[2][2];
452 S[6] = dtK[0][0];
453 S[7] = dtK[0][1];
454 S[8] = dtK[0][2];
455 S[9] = dtK[1][1];
456 S[10] = dtK[1][2];
457 S[11] = dtK[2][2];
458 S[12] = dtTheta;
459 for (int i = 0; i < 3; i++) S[13+i] = dtGhat[i];
460 S[16] = dtalpha;
461 for (int i = 0; i < 3; i++) S[17+i] = dtbeta[i];
462 for (int i = 0; i < 3; i++) S[20+i] = dtbb[i];
463 for (int i = 0; i < 3; i++) S[23+i] = dtA[i];
464 S[26] = dtB[0][0];
465 S[27] = dtB[1][0];
466 S[28] = dtB[2][0];
467 S[29] = dtB[0][1];
468 S[30] = dtB[1][1];
469 S[31] = dtB[2][1];
470 S[32] = dtB[0][2];
471 S[33] = dtB[1][2];
472 S[34] = dtB[2][2];
473 S[35] = dtD[0][0][0];
474 S[36] = dtD[0][0][1];
475 S[37] = dtD[0][0][2];
476 S[38] = dtD[0][1][1];
477 S[39] = dtD[0][1][2];
478 S[40] = dtD[0][2][2];
479 S[41] = dtD[1][0][0];
480 S[42] = dtD[1][0][1];
481 S[43] = dtD[1][0][2];
482 S[44] = dtD[1][1][1];
483 S[45] = dtD[1][1][2];
484 S[46] = dtD[1][2][2];
485 S[47] = dtD[2][0][0];
486 S[48] = dtD[2][0][1];
487 S[49] = dtD[2][0][2];
488 S[50] = dtD[2][1][1];
489 S[51] = dtD[2][1][2];
490 S[52] = dtD[2][2][2];
491 S[53] = dtTraceK;
492 S[54] = dtphi;
493 for (int i = 0; i < 3; i++) S[55+i] = dtP[i];
494 S[58] = 0;
495}
496#if defined(GPUOffloadingOMP)
497#pragma omp end declare target
498#endif
499
500
501#if defined(GPUOffloadingOMP)
502#pragma omp declare target
503#endif
504void examples::exahype2::mgccz4::ncp(double* BgradQ, const double* const Q, const double* const gradQSerialised, const int normal,
505 const int MGCCZ4LapseType,
506 const double MGCCZ4ds,
507 const double MGCCZ4c,
508 const double MGCCZ4e,
509 const double MGCCZ4f,
510 const double MGCCZ4bs,
511 const double MGCCZ4sk,
512 const double MGCCZ4xi,
513 const double MGCCZ4mu
514 )
515{
516 const double alpha = std::exp(std::fmax(-20., std::fmin(20.,Q[16])));
517 //printf("alpha %f\n",alpha);
518 const double alpha2 = alpha*alpha;
519 double fa = 1.0;
520 double faa = 0.0;
521 if (MGCCZ4LapseType==1)
522 {
523 fa = 2./alpha;
524 faa = -fa/alpha;
525 }
526
527 constexpr int nVar(64);
528
529 double gradQin[64][3]={0}; //={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
530
531 // De-serialise input data and fill static array
532 // FIXME the use of 2D arrays can be avoided: all terms not in the normal are 0
533 for (int i=0; i<nVar; i++) gradQin[i][normal] = gradQSerialised[i+normal*nVar];
534
535
536 // Note g_cov is symmetric
537 const double g_cov[3][3] = { {Q[0], Q[1], Q[2]}, {Q[1], Q[3], Q[4]}, {Q[2], Q[4], Q[5]} };
538 const double invdet = 1./( Q[0]*Q[3]*Q[5] - Q[0]*Q[4]*Q[4] - Q[1]*Q[1]*Q[5] + 2*Q[1]*Q[2]*Q[4] -Q[2]*Q[2]*Q[3]);
539
540 const double g_contr[3][3] = {
541 { ( Q[3]*Q[5]-Q[4]*Q[4])*invdet, -( Q[1]*Q[5]-Q[2]*Q[4])*invdet, -(-Q[1]*Q[4]+Q[2]*Q[3])*invdet},
542 {-( Q[1]*Q[5]-Q[4]*Q[2])*invdet, ( Q[0]*Q[5]-Q[2]*Q[2])*invdet, -( Q[0]*Q[4]-Q[2]*Q[1])*invdet},
543 {-(-Q[1]*Q[4]+Q[3]*Q[2])*invdet, -( Q[0]*Q[4]-Q[1]*Q[2])*invdet, ( Q[0]*Q[3]-Q[1]*Q[1])*invdet}
544 };
545
546
547 // NOTE Aex is symmetric
548 double Aex[3][3] = { {Q[6], Q[7], Q[8]}, {Q[7], Q[9], Q[10]}, {Q[8], Q[10], Q[11]} };
549
550 double traceA = 0;
551 for (int i=0;i<3;i++)
552 for (int j=0;j<3;j++) traceA+=g_contr[i][j]*Aex[i][j];
553 traceA *= 1./3;
554
555 for (int i=0;i<3;i++)
556 for (int j=0;j<3;j++) Aex[i][j] -= traceA * g_cov[i][j];
557
558 // Matrix multiplications Amix = matmul(g_contr, Aex) Aup = matmul(g_contr, mytranspose(Amix))
559 double Amix[3][3]={0,0,0,0,0,0,0,0,0};
560 for (int i = 0; i < 3; i++)
561 for (int j = 0; j < 3; j++)
562 for (int u = 0; u < 3; u++) Amix[i][j] += g_contr[i][u] * Aex[u][j];
563 double Aup[3][3]={0,0,0,0,0,0,0,0,0};
564 for (int i = 0; i < 3; i++)
565 for (int j = 0; j < 3; j++)
566 for (int u = 0; u < 3; u++) Aup[i][j] += g_contr[i][u] * Amix[j][u]; // Note the transposition is in the indices
567
568 const double DD[3][3][3] = {
569 {{Q[35], Q[36], Q[37]}, {Q[36], Q[38], Q[39]}, {Q[37], Q[39], Q[40]}},
570 {{Q[41], Q[42], Q[43]}, {Q[42], Q[44], Q[45]}, {Q[43], Q[45], Q[46]}},
571 {{Q[47], Q[48], Q[49]}, {Q[48], Q[50], Q[51]}, {Q[49], Q[51], Q[52]}}
572 };
573 double dgup[3][3][3] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
574 for (int k = 0; k < 3; k++)
575 for (int m = 0; m < 3; m++)
576 for (int l = 0; l < 3; l++)
577 for (int n = 0; n < 3; n++)
578 for (int j = 0; j < 3; j++) dgup[k][m][l] -= g_contr[m][n]*g_contr[j][l]*2*DD[k][n][j];
579
580 const double PP[3] = {Q[55], Q[56], Q[57]};
581
582 double Christoffel[3][3][3] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
583 double Christoffel_tilde[3][3][3] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
584 //double Christoffel_kind1[3][3][3] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
585
586 for (int j = 0; j < 3; j++)
587 for (int i = 0; i < 3; i++)
588 for (int k = 0; k < 3; k++)
589 {
590 //Christoffel_kind1[i][j][k] = DD[k][i][j] + DD[j][i][k] - DD[i][j][k];
591 for (int l = 0; l < 3; l++)
592 {
593 Christoffel_tilde[i][j][k] += g_contr[k][l] * ( DD[i][j][l] + DD[j][i][l] - DD[l][i][j] );
594 Christoffel[i][j][k] += g_contr[k][l] * ( DD[i][j][l] + DD[j][i][l] - DD[l][i][j] ) - g_contr[k][l] * ( g_cov[j][l] * PP[i] + g_cov[i][l] * PP[j] - g_cov[i][j] * PP[l] );
595 }
596 }
597
598 double Gtilde[3] = {0,0,0};
599 for (int l = 0; l < 3; l++)
600 for (int j = 0; j < 3; j++)
601 for (int i = 0; i < 3; i++) Gtilde[i] += g_contr[j][l] * Christoffel_tilde[j][l][i];
602
603 const double Ghat[3] = {Q[13], Q[14], Q[15]};
604
605 const double phi = std::exp(std::fmax(-20., std::fmin(20.,Q[54])));
606 const double phi2 = phi*phi;
607
608 double Z[3] = {0,0,0};
609 for (int i=0;i<3;i++)
610 for (int j=0;j<3;j++) Z[i] += 0.5*( g_cov[i][j]* (Ghat[j] - Gtilde[j]));
611 double Zup[3] = {0,0,0};
612 for (int i=0;i<3;i++)
613 for (int j=0;j<3;j++) Zup[i] += phi2 * g_contr[i][j] * Z[j];
614
615
616 const double dDD[3][3][3][3] = {
617 {
618 {
619 {gradQin[35][0],gradQin[36][0],gradQin[37][0]}, {gradQin[36][0],gradQin[38][0],gradQin[39][0]}, {gradQin[37][0],gradQin[39][0],gradQin[40][0]},
620 },
621 {
622 {gradQin[41][0],gradQin[42][0],gradQin[43][0]}, {gradQin[42][0],gradQin[44][0],gradQin[45][0]}, {gradQin[43][0],gradQin[45][0],gradQin[46][0]},
623 },
624 {
625 {gradQin[47][0],gradQin[48][0],gradQin[49][0]}, {gradQin[48][0],gradQin[50][0],gradQin[51][0]}, {gradQin[49][0],gradQin[51][0],gradQin[52][0]}
626 }
627 },
628 {
629 {
630 {gradQin[35][1],gradQin[36][1],gradQin[37][1]}, {gradQin[36][1],gradQin[38][1],gradQin[39][1]}, {gradQin[37][1],gradQin[39][1],gradQin[40][1]},
631 },
632 {
633 {gradQin[41][1],gradQin[42][1],gradQin[43][1]}, {gradQin[42][1],gradQin[44][1],gradQin[45][1]}, {gradQin[43][1],gradQin[45][1],gradQin[46][1]},
634 },
635 {
636 {gradQin[47][1],gradQin[48][1],gradQin[49][1]}, {gradQin[48][1],gradQin[50][1],gradQin[51][1]}, {gradQin[49][1],gradQin[51][1],gradQin[52][1]}
637 }
638 },
639 {
640 {
641 {gradQin[35][2],gradQin[36][2],gradQin[37][2]}, {gradQin[36][2],gradQin[38][2],gradQin[39][2]}, {gradQin[37][2],gradQin[39][2],gradQin[40][2]},
642 },
643 {
644 {gradQin[41][2],gradQin[42][2],gradQin[43][2]}, {gradQin[42][2],gradQin[44][2],gradQin[45][2]}, {gradQin[43][2],gradQin[45][2],gradQin[46][2]},
645 },
646 {
647 {gradQin[47][2],gradQin[48][2],gradQin[49][2]}, {gradQin[48][2],gradQin[50][2],gradQin[51][2]}, {gradQin[49][2],gradQin[51][2],gradQin[52][2]}
648 }
649 }
650 };
651
652
653 const double dPP[3][3] = {
654 {gradQin[55][0],gradQin[56][0],gradQin[57][0]},
655 {gradQin[55][1],gradQin[56][1],gradQin[57][1]},
656 {gradQin[55][2],gradQin[56][2],gradQin[57][2]}
657 };
658
659
660 double dChristoffelNCP[3][3][3][3] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
661 double dChristoffel_tildeNCP[3][3][3][3] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
662 for (int i = 0; i < 3; i++)
663 for (int j = 0; j < 3; j++)
664 for (int m = 0; m < 3; m++)
665 for (int k = 0; k < 3; k++)
666 {
667 dChristoffelNCP [k][i][j][m] = 0;
668 dChristoffel_tildeNCP[k][i][j][m] = 0;
669 for (int l = 0; l < 3; l++)
670 {
671 dChristoffelNCP[k][i][j][m] += 0.5*g_contr[m][l] * (
672 dDD[k][i][j][l] + dDD[i][k][j][l] + dDD[k][j][i][l] + dDD[j][k][i][l] - dDD[k][l][i][j] - dDD[l][k][i][j]
673 - g_cov[j][l]*(dPP[k][i] + dPP[i][k]) - g_cov[i][l]*(dPP[k][j]+dPP[j][k]) + g_cov[i][j]*(dPP[k][l]+dPP[l][k]) );
674
675 dChristoffel_tildeNCP[k][i][j][m] += 0.5*g_contr[m][l]*(dDD[k][i][j][l] + dDD[i][k][j][l] + dDD[k][j][i][l] + dDD[j][k][i][l] - dDD[k][l][i][j] - dDD[l][k][i][j]);
676
677 }
678 }
679
680 double RiemannNCP[3][3][3][3] = {0};
681 for (int i = 0; i < 3; i++)
682 for (int j = 0; j < 3; j++)
683 for (int m = 0; m < 3; m++)
684 for (int k = 0; k < 3; k++) RiemannNCP[i][k][j][m] = dChristoffelNCP[k][i][j][m] - dChristoffelNCP[j][i][k][m];
685
686 double RicciNCP[3][3] = {0,0,0,0,0,0,0,0,0};
687 for (int m = 0; m < 3; m++)
688 for (int n = 0; n < 3; n++)
689 for (int l = 0; l < 3; l++) RicciNCP[m][n] += RiemannNCP[m][l][n][l];
690
691 double dGtildeNCP[3][3] = {0,0,0,0,0,0,0,0,0};
692 for (int i = 0; i < 3; i++)
693 for (int k = 0; k < 3; k++)
694 for (int j = 0; j < 3; j++)
695 for (int l = 0; l < 3; l++) dGtildeNCP[k][i] += g_contr[j][l]*dChristoffel_tildeNCP[k][j][l][i];
696
697
698 const double dGhat[3][3] = {
699 {gradQin[13][0],gradQin[14][0],gradQin[15][0]},
700 {gradQin[13][1],gradQin[14][1],gradQin[15][1]},
701 {gradQin[13][2],gradQin[14][2],gradQin[15][2]}
702 };
703
704
705
706 double dZNCP[3][3] = {0,0,0,0,0,0,0,0,0};
707 for (int j = 0; j < 3; j++)
708 for (int i = 0; i < 3; i++)
709 for (int k = 0; k < 3; k++) dZNCP[k][i] += MGCCZ4ds*0.5*g_cov[i][j]*(dGhat[k][j]-dGtildeNCP[k][j]);
710
711
712 double RicciPlusNablaZNCP[3][3];
713 for (int i = 0; i < 3; i++)
714 for (int j = 0; j < 3; j++) RicciPlusNablaZNCP[i][j] = RicciNCP[i][j] + dZNCP[i][j] + dZNCP[j][i];
715
716 double RPlusTwoNablaZNCP = 0;
717 for (int i = 0; i < 3; i++)
718 for (int j = 0; j < 3; j++) RPlusTwoNablaZNCP += g_contr[i][j]*RicciPlusNablaZNCP[i][j]; // TODO fuse these steps
719 RPlusTwoNablaZNCP*=phi2;
720
721
722 const double AA[3] = {Q[23], Q[24], Q[25]};
723 const double dAA[3][3] = {
724 {gradQin[23][0],gradQin[24][0],gradQin[25][0]},
725 {gradQin[23][1],gradQin[24][1],gradQin[25][1]},
726 {gradQin[23][2],gradQin[24][2],gradQin[25][2]}
727 };
728
729 double nablaijalphaNCP[3][3];
730 for (int i = 0; i < 3; i++)
731 for (int j = 0; j < 3; j++) nablaijalphaNCP[i][j] = alpha*0.5*( dAA[i][j] + dAA[j][i] );
732
733 double nablanablaalphaNCP = 0;
734 for (int i = 0; i < 3; i++)
735 for (int j = 0; j < 3; j++) nablanablaalphaNCP += g_contr[i][j]*nablaijalphaNCP[i][j];
736 nablanablaalphaNCP*=phi2;
737
738
739 double SecondOrderTermsNCP[3][3];
740 for (int i = 0; i < 3; i++)
741 for (int j = 0; j < 3; j++) SecondOrderTermsNCP[i][j] = -nablaijalphaNCP[i][j] + alpha*RicciPlusNablaZNCP[i][j];
742
743 double traceNCP = 0;
744 for (int i = 0; i < 3; i++)
745 for (int j = 0; j < 3; j++) traceNCP += g_contr[i][j]*SecondOrderTermsNCP[i][j];
746
747 for (int i = 0; i < 3; i++)
748 for (int j = 0; j < 3; j++) SecondOrderTermsNCP[i][j] -= 1./3 * traceNCP * g_cov[i][j];
749
750 const double beta[3] = {Q[17], Q[18], Q[19]};
751
752 const double dAex[3][3][3] = {
753 {{gradQin[6][0],gradQin[7][0],gradQin[8][0]}, {gradQin[7][0], gradQin[9][0], gradQin[10][0]}, {gradQin[8][0], gradQin[10][0], gradQin[11][0]}},
754 {{gradQin[6][1],gradQin[7][1],gradQin[8][1]}, {gradQin[7][1], gradQin[9][1], gradQin[10][1]}, {gradQin[8][1], gradQin[10][1], gradQin[11][1]}},
755 {{gradQin[6][2],gradQin[7][2],gradQin[8][2]}, {gradQin[7][2], gradQin[9][2], gradQin[10][2]}, {gradQin[8][2], gradQin[10][2], gradQin[11][2]}}
756 };
760 double dtK[3][3];
761 for (int i = 0; i < 3; i++)
762 for (int j = 0; j < 3; j++) dtK[i][j] = phi2*SecondOrderTermsNCP[i][j] + beta[0] * dAex[0][i][j] + beta[1] * dAex[1][i][j] + beta[2] * dAex[2][i][j]; // extrinsic curvature
763
764 const double dtraceK[3] = {gradQin[53][0], gradQin[53][1], gradQin[53][2]};
765
766 double dtTraceK = -nablanablaalphaNCP + alpha*RPlusTwoNablaZNCP + beta[0]*dtraceK[0] + beta[1]*dtraceK[1] + beta[2]*dtraceK[2];
767
768 const double BB[3][3] = {
769 {Q[26], Q[27], Q[28]}, {Q[29], Q[30], Q[31]}, {Q[32], Q[33], Q[34]}
770 };
771
772 double traceB = BB[0][0] + BB[1][1] + BB[2][2]; // TODO direct from Q!
773
774 double Aupdown = 0;
775 for (int i = 0; i < 3; i++)
776 for (int j = 0; j < 3; j++) Aupdown += Aex[i][j]*Aup[i][j];
777
778 const double dTheta[3] = {gradQin[12][0],gradQin[12][1],gradQin[12][2]};
779 const double dtTheta = 0.5*alpha*MGCCZ4e*MGCCZ4e*( RPlusTwoNablaZNCP ) + beta[0]*dTheta[0] + beta[1]*dTheta[1] + beta[2]*dTheta[2]; // *** original cleaning ***
780
781 double divAupNCP[3] = {0,0,0};
782
783 for (int i = 0; i < 3; i++)
784 for (int j = 0; j < 3; j++)
785 for (int l = 0; l < 3; l++)
786 for (int k = 0; k < 3; k++) divAupNCP[i] += g_contr[i][l]*g_contr[j][k]*dAex[j][l][k];
787
788
789 const double dBB[3][3][3] = {
790 {
791 {MGCCZ4sk*gradQin[26][0],MGCCZ4sk*gradQin[27][0],MGCCZ4sk*gradQin[28][0]}, {MGCCZ4sk*gradQin[29][0],MGCCZ4sk*gradQin[30][0],MGCCZ4sk*gradQin[31][0]}, {MGCCZ4sk*gradQin[32][0],MGCCZ4sk*gradQin[33][0],MGCCZ4sk*gradQin[34][0]}
792 },
793 {
794 {MGCCZ4sk*gradQin[26][1],MGCCZ4sk*gradQin[27][1],MGCCZ4sk*gradQin[28][1]}, {MGCCZ4sk*gradQin[29][1],MGCCZ4sk*gradQin[30][1],MGCCZ4sk*gradQin[31][1]}, {MGCCZ4sk*gradQin[32][1],MGCCZ4sk*gradQin[33][1],MGCCZ4sk*gradQin[34][1]}
795 },
796 {
797 {MGCCZ4sk*gradQin[26][2],MGCCZ4sk*gradQin[27][2],MGCCZ4sk*gradQin[28][2]}, {MGCCZ4sk*gradQin[29][2],MGCCZ4sk*gradQin[30][2],MGCCZ4sk*gradQin[31][2]}, {MGCCZ4sk*gradQin[32][2],MGCCZ4sk*gradQin[33][2],MGCCZ4sk*gradQin[34][2]}
798 }
799 };
800
801 double dtGhat[3] = {0,0,0};
802 for (int i = 0; i < 3; i++)
803 {
804 double temp=0, temp2=0;
805 for (int j = 0; j < 3; j++)
806 {
807 temp +=g_contr[i][j]*dtraceK[j];
808 temp2 +=g_contr[j][i]*dTheta[j];
809 }
810 dtGhat[i] = -4./3.*alpha*temp + 2*alpha*temp2 + beta[0]*dGhat[0][i] + beta[1]*dGhat[1][i] + beta[2]*dGhat[2][i];
811 for (int l = 0; l < 3; l++)
812 for (int k = 0; k < 3; k++) dtGhat[i] += g_contr[k][l]*0.5*(dBB[k][l][i] + dBB[l][k][i]) + 1./3.*g_contr[i][k]*0.5*(dBB[k][l][l] + dBB[l][k][l]);
813 }
814
815 double ov[3];
816 for (int k = 0; k < 3; k++)
817 {
818 double temp=0;
819 for (int i = 0; i < 3; i++)
820 for (int j = 0; j < 3; j++) temp += g_contr[i][j]*dAex[k][i][j];
821 ov[k] = 2*alpha*temp;
822 }
823 for (int i = 0; i < 3; i++)
824 for (int j = 0; j < 3; j++) dtGhat[i] += MGCCZ4sk*g_contr[i][j]*ov[j];
825
826 double dtbb[3];
827 for (int i = 0; i < 3; i++)
828 {
829 dtbb[i] = MGCCZ4xi*dtGhat[i] + MGCCZ4bs * ( beta[0]*gradQin[20+i][0] + beta[1]*gradQin[20+i][1] + beta[2]*gradQin[20+i][2] - beta[0]*gradQin[13+i][0] - beta[1]*gradQin[13+i][1] - beta[2]*gradQin[13+i][2]);
830 dtbb[i] *= MGCCZ4sk;
831 }
832
833 // Auxiliary variables
834 double dtA[3];
835 double dK0[3] = {0,0,0};
836 for (int i = 0; i < 3; i++)
837 {
838 dtA[i] = -alpha*fa*(dtraceK[i] - dK0[i] - MGCCZ4c*2*dTheta[i]) + beta[0]*dAA[0][i] + beta[1]*dAA[1][i] + beta[2]*dAA[2][i];
839 }
840 /*if ((dtA[1]+dtA[2])!=0 && dtA[2]>1e-9) {printf("cpp dtA[1]=%g, dtA[2]=%g\n",dtA[1],dtA[2]);
841 printf("cpp dtheta[1]=%g, dtheta[2]=%g \n",dTheta[1],dTheta[2]);
842 printf("cpp dtracek[1]=%g, dtracek[2]=%g \n",dtraceK[1],dtraceK[2]);
843 }*/
844 //if ((dTheta[0]+dTheta[1]+dTheta[2])!=0) printf("cpp dtheta[0] = %g, dtheta[1] = %g, dtheta[2] = %g \n\n",dTheta[0],dTheta[1],dTheta[2]);
845 for (int k = 0; k < 3; k++)
846 {
847 double temp=0;
848 for (int i = 0; i < 3; i++)
849 for (int j = 0; j < 3; j++) temp += g_contr[i][j]*dAex[k][i][j]; // TODO we computed this quantity few lines earlier alrady
850 dtA[k] -= MGCCZ4sk*alpha*fa*temp;
851 }
852
853 double dtB[3][3] = {
854 {MGCCZ4f*gradQin[20][0],MGCCZ4f*gradQin[21][0],MGCCZ4f*gradQin[22][0]},
855 {MGCCZ4f*gradQin[20][1],MGCCZ4f*gradQin[21][1],MGCCZ4f*gradQin[22][1]},
856 {MGCCZ4f*gradQin[20][2],MGCCZ4f*gradQin[21][2],MGCCZ4f*gradQin[22][2]}
857 };
858
859
860 for (int i = 0; i < 3; i++)
861 for (int k = 0; k < 3; k++)
862 for (int j = 0; j < 3; j++)
863 {
864 dtB[k][i] += MGCCZ4mu*alpha2 * g_contr[i][j]*( dPP[k][j] - dPP[j][k]);
865 for (int n = 0; n < 3; n++)
866 for (int l = 0; l < 3; l++) dtB[k][i] -= MGCCZ4mu*alpha2 * g_contr[i][j]*g_contr[n][l]*( dDD[k][l][j][n] - dDD[l][k][j][n]);
867 }
868
869 for (int i = 0; i < 3; i++)
870 for (int j = 0; j < 3; j++) dtB[i][j] += MGCCZ4bs*(beta[0]*dBB[0][i][j] + beta[1]*dBB[1][i][j] + beta[2]*dBB[2][i][j]);
871
872 // NOTE 0 value param
873 for (int i = 0; i < 3; i++)
874 for (int j = 0; j < 3; j++) dtB[i][j] *= MGCCZ4sk;
875
876 double dtD[3][3][3];
877 for (int i = 0; i < 3; i++)
878 for (int j = 0; j < 3; j++)
879 for (int k = 0; k < 3; k++) dtD[i][j][k] = -alpha*dAex[i][j][k];
880
881 for (int i = 0; i < 3; i++)
882 for (int j = 0; j < 3; j++)
883 for (int k = 0; k < 3; k++)
884 for (int m = 0; m < 3; m++)
885 {
886 dtD[k][i][j] += ( 0.25*(g_cov[m][i]*(dBB[k][j][m] + dBB[j][k][m]) + g_cov[m][j]*(dBB[k][i][m]+dBB[i][k][m])) - 1./6.*g_cov[i][j]*(dBB[k][m][m]+dBB[m][k][m]) );
887 for (int n = 0; n < 3; n++)
888 dtD[k][i][j] += 1./3*alpha*g_cov[i][j]*g_contr[n][m]*dAex[k][n][m]; // explicitly remove the trace of tilde A again
889 }
890
891 for (int i = 0; i < 3; i++)
892 for (int j = 0; j < 3; j++)
893 for (int k = 0; k < 3; k++) dtD[i][j][k] += beta[0]*dDD[0][i][j][k] + beta[1]*dDD[1][i][j][k] + beta[2]*dDD[2][i][j][k];
894
895 double dtP[3];
896 for (int i = 0; i < 3; i++) dtP[i] = beta[0]*dPP[0][i] + beta[1]*dPP[1][i] + beta[2]*dPP[2][i];
897 for (int k = 0; k < 3; k++)
898 {
899 double temp=0;
900 for (int m = 0; m < 3; m++)
901 for (int n = 0; n < 3; n++) temp += g_contr[m][n]*dAex[k][m][n];
902 dtP[k] += 1./3*alpha*(dtraceK[k] + MGCCZ4sk*temp);
903 for (int i = 0; i < 3; i++) dtP[k] -= 1./6*(dBB[k][i][i] + dBB[i][k][i]);
904 }
905
906 double dtgamma[3][3] = {0,0,0,0,0,0,0,0,0};
907 double dtalpha = 0;
908 double dtbeta[3] = {0,0,0};
909 double dtphi = 0;
910
911 BgradQ[0] = -dtgamma[0][0];
912 BgradQ[1] = -dtgamma[0][1];
913 BgradQ[2] = -dtgamma[0][2];
914 BgradQ[3] = -dtgamma[1][1];
915 BgradQ[4] = -dtgamma[1][2];
916 BgradQ[5] = -dtgamma[2][2];
917 BgradQ[6] = -dtK[0][0];
918 BgradQ[7] = -dtK[0][1];
919 BgradQ[8] = -dtK[0][2];
920 BgradQ[9] = -dtK[1][1];
921 BgradQ[10] = -dtK[1][2];
922 BgradQ[11] = -dtK[2][2]; // ok
923 BgradQ[12] = -dtTheta; // ok
924 for (int i = 0; i < 3; i++) BgradQ[13+i] = -dtGhat[i]; // buggy
925 BgradQ[16] = -dtalpha;
926 for (int i = 0; i < 3; i++) BgradQ[17+i] = -dtbeta[i];
927 for (int i = 0; i < 3; i++) BgradQ[20+i] = -dtbb[i];
928 for (int i = 0; i < 3; i++) BgradQ[23+i] = -dtA[i];
929 BgradQ[26] = -dtB[0][0]; // note: thes are all 0 for default MGCCZ4sk=0
930 BgradQ[27] = -dtB[1][0];
931 BgradQ[28] = -dtB[2][0];
932 BgradQ[29] = -dtB[0][1];
933 BgradQ[30] = -dtB[1][1];
934 BgradQ[31] = -dtB[2][1];
935 BgradQ[32] = -dtB[0][2];
936 BgradQ[33] = -dtB[1][2];
937 BgradQ[34] = -dtB[2][2];
938
939 BgradQ[35] = -dtD[0][0][0];
940 BgradQ[36] = -dtD[0][0][1];
941 BgradQ[37] = -dtD[0][0][2];
942 BgradQ[38] = -dtD[0][1][1];
943 BgradQ[39] = -dtD[0][1][2];
944 BgradQ[40] = -dtD[0][2][2];
945
946 BgradQ[41] = -dtD[1][0][0];
947 BgradQ[42] = -dtD[1][0][1];
948 BgradQ[43] = -dtD[1][0][2];
949 BgradQ[44] = -dtD[1][1][1];
950 BgradQ[45] = -dtD[1][1][2];
951 BgradQ[46] = -dtD[1][2][2];
952
953 BgradQ[47] = -dtD[2][0][0];
954 BgradQ[48] = -dtD[2][0][1];
955 BgradQ[49] = -dtD[2][0][2];
956 BgradQ[50] = -dtD[2][1][1];
957 BgradQ[51] = -dtD[2][1][2];
958 BgradQ[52] = -dtD[2][2][2];
959 BgradQ[53] = -dtTraceK;
960 BgradQ[54] = -dtphi;
961 for (int i = 0; i < 3; i++) BgradQ[55+i] = -dtP[i];
962 BgradQ[58] = 0;
963 //if (dtP[1]!=0) printf("dtP[1] = %g, dtP[2] = %g \n\n",dtP[1],dtP[2]);
964}
965#if defined(GPUOffloadingOMP)
966#pragma omp end declare target
967#endif
968
969#if defined(GPUOffloadingOMP)
970#pragma omp declare target
971#endif
972void examples::exahype2::mgccz4::admconstraints(double* constraints, const double* const Q, const double* const gradQSerialised)
973{
974 constexpr int nVar(64);
975 double gradQin[64][3]={0};
976
977 // De-serialise input data and fill static array
978 for (int normal=0; normal<3; normal++)
979 for (int i=0; i<nVar; i++) gradQin[i][normal] = gradQSerialised[i+normal*nVar];
980
981 // Note g_cov is symmetric
982 const double g_cov[3][3] = { {Q[0], Q[1], Q[2]}, {Q[1], Q[3], Q[4]}, {Q[2], Q[4], Q[5]} };
983 const double invdet = 1./( Q[0]*Q[3]*Q[5] - Q[0]*Q[4]*Q[4] - Q[1]*Q[1]*Q[5] + 2*Q[1]*Q[2]*Q[4] -Q[2]*Q[2]*Q[3]);
984
985 const double g_contr[3][3] = {
986 { ( Q[3]*Q[5]-Q[4]*Q[4])*invdet, -( Q[1]*Q[5]-Q[2]*Q[4])*invdet, -(-Q[1]*Q[4]+Q[2]*Q[3])*invdet},
987 {-( Q[1]*Q[5]-Q[4]*Q[2])*invdet, ( Q[0]*Q[5]-Q[2]*Q[2])*invdet, -( Q[0]*Q[4]-Q[2]*Q[1])*invdet},
988 {-(-Q[1]*Q[4]+Q[3]*Q[2])*invdet, -( Q[0]*Q[4]-Q[1]*Q[2])*invdet, ( Q[0]*Q[3]-Q[1]*Q[1])*invdet}
989 };
990
991 // NOTE Aex is symmetric
992 double Aex[3][3] = { {Q[6], Q[7], Q[8]}, {Q[7], Q[9], Q[10]}, {Q[8], Q[10], Q[11]} };
993
994 double traceA = 0;
995 for (int i=0;i<3;i++)
996 for (int j=0;j<3;j++) traceA+=g_contr[i][j]*Aex[i][j];
997 traceA *= 1./3;
998
999 for (int i=0;i<3;i++)
1000 for (int j=0;j<3;j++) Aex[i][j] -= traceA * g_cov[i][j];
1001
1002 const double dAex[3][3][3] = {
1003 {{gradQin[6][0],gradQin[7][0],gradQin[8][0]}, {gradQin[7][0], gradQin[9][0], gradQin[10][0]}, {gradQin[8][0], gradQin[10][0], gradQin[11][0]}},
1004 {{gradQin[6][1],gradQin[7][1],gradQin[8][1]}, {gradQin[7][1], gradQin[9][1], gradQin[10][1]}, {gradQin[8][1], gradQin[10][1], gradQin[11][1]}},
1005 {{gradQin[6][2],gradQin[7][2],gradQin[8][2]}, {gradQin[7][2], gradQin[9][2], gradQin[10][2]}, {gradQin[8][2], gradQin[10][2], gradQin[11][2]}}
1006 };
1007
1008 const double traceK = Q[53];
1009 const double dtraceK[3] = {gradQin[53][0], gradQin[53][1], gradQin[53][2]};
1010
1011 const double phi = std::exp(std::fmax(-20., std::fmin(20.,Q[54])));
1012 const double phi2 = phi*phi;
1013 const double PP[3] = {Q[55], Q[56], Q[57]};
1014
1015 const double dPP[3][3] = {
1016 {gradQin[55][0],gradQin[56][0],gradQin[57][0]},
1017 {gradQin[55][1],gradQin[56][1],gradQin[57][1]},
1018 {gradQin[55][2],gradQin[56][2],gradQin[57][2]}
1019 };
1020
1021 const double DD[3][3][3] = {
1022 {{Q[35], Q[36], Q[37]}, {Q[36], Q[38], Q[39]}, {Q[37], Q[39], Q[40]}},
1023 {{Q[41], Q[42], Q[43]}, {Q[42], Q[44], Q[45]}, {Q[43], Q[45], Q[46]}},
1024 {{Q[47], Q[48], Q[49]}, {Q[48], Q[50], Q[51]}, {Q[49], Q[51], Q[52]}}
1025 };
1026 const double dDD[3][3][3][3] = {
1027 {
1028 {
1029 {gradQin[35][0],gradQin[36][0],gradQin[37][0]}, {gradQin[36][0],gradQin[38][0],gradQin[39][0]}, {gradQin[37][0],gradQin[39][0],gradQin[40][0]},
1030 },
1031 {
1032 {gradQin[41][0],gradQin[42][0],gradQin[43][0]}, {gradQin[42][0],gradQin[44][0],gradQin[45][0]}, {gradQin[43][0],gradQin[45][0],gradQin[46][0]},
1033 },
1034 {
1035 {gradQin[47][0],gradQin[48][0],gradQin[49][0]}, {gradQin[48][0],gradQin[50][0],gradQin[51][0]}, {gradQin[49][0],gradQin[51][0],gradQin[52][0]}
1036 }
1037 },
1038 {
1039 {
1040 {gradQin[35][1],gradQin[36][1],gradQin[37][1]}, {gradQin[36][1],gradQin[38][1],gradQin[39][1]}, {gradQin[37][1],gradQin[39][1],gradQin[40][1]},
1041 },
1042 {
1043 {gradQin[41][1],gradQin[42][1],gradQin[43][1]}, {gradQin[42][1],gradQin[44][1],gradQin[45][1]}, {gradQin[43][1],gradQin[45][1],gradQin[46][1]},
1044 },
1045 {
1046 {gradQin[47][1],gradQin[48][1],gradQin[49][1]}, {gradQin[48][1],gradQin[50][1],gradQin[51][1]}, {gradQin[49][1],gradQin[51][1],gradQin[52][1]}
1047 }
1048 },
1049 {
1050 {
1051 {gradQin[35][2],gradQin[36][2],gradQin[37][2]}, {gradQin[36][2],gradQin[38][2],gradQin[39][2]}, {gradQin[37][2],gradQin[39][2],gradQin[40][2]},
1052 },
1053 {
1054 {gradQin[41][2],gradQin[42][2],gradQin[43][2]}, {gradQin[42][2],gradQin[44][2],gradQin[45][2]}, {gradQin[43][2],gradQin[45][2],gradQin[46][2]},
1055 },
1056 {
1057 {gradQin[47][2],gradQin[48][2],gradQin[49][2]}, {gradQin[48][2],gradQin[50][2],gradQin[51][2]}, {gradQin[49][2],gradQin[51][2],gradQin[52][2]}
1058 }
1059 }
1060 };
1061
1062 double dgup[3][3][3] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1063 for (int k = 0; k < 3; k++)
1064 for (int m = 0; m < 3; m++)
1065 for (int l = 0; l < 3; l++)
1066 for (int n = 0; n < 3; n++)
1067 for (int j = 0; j < 3; j++) dgup[k][m][l] -= g_contr[m][n]*g_contr[j][l]*2*DD[k][n][j];
1068
1069 double Kex[3][3]={ 0 };
1070 for (int i=0;i<3;i++)
1071 for (int j=0;j<3;j++) Kex[i][j]=Aex[i][j]/phi2 + (1./3)*traceK*g_cov[i][j]/phi2;
1072 double Kmix[3][3]={0,0,0,0,0,0,0,0,0};
1073 for (int i = 0; i < 3; i++)
1074 for (int j = 0; j < 3; j++)
1075 for (int u = 0; u < 3; u++) Kmix[i][j] += phi2*g_contr[i][u] * Kex[u][j];
1076 double Kup[3][3]={0,0,0,0,0,0,0,0,0};
1077 for (int i = 0; i < 3; i++)
1078 for (int j = 0; j < 3; j++)
1079 for (int u = 0; u < 3; u++) Kup[i][j] += phi2*g_contr[i][u] * Kmix[j][u]; // Note the transposition is in the indices
1080
1081 double Christoffel[3][3][3] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1082 for (int j = 0; j < 3; j++)
1083 for (int i = 0; i < 3; i++)
1084 for (int k = 0; k < 3; k++)
1085 for (int l = 0; l < 3; l++) Christoffel[i][j][k] += g_contr[k][l] * ( DD[i][j][l] + DD[j][i][l] - DD[l][i][j] ) - g_contr[k][l] * ( g_cov[j][l] * PP[i] + g_cov[i][l] * PP[j] - g_cov[i][j] * PP[l] );
1086
1087 double dChristoffel[3][3][3][3] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
1088
1089 for (int i = 0; i < 3; i++)
1090 for (int j = 0; j < 3; j++)
1091 for (int m = 0; m < 3; m++)
1092 for (int k = 0; k < 3; k++)
1093 for (int l = 0; l < 3; l++)
1094 {
1095 dChristoffel[k][i][j][m] += 0.5*g_contr[m][l] * (
1096 dDD[k][i][j][l] + dDD[i][k][j][l] + dDD[k][j][i][l] + dDD[j][k][i][l] - dDD[k][l][i][j] - dDD[l][k][i][j]
1097 - g_cov[j][l]*(dPP[k][i] + dPP[i][k]) - g_cov[i][l]*(dPP[k][j]+dPP[j][k]) + g_cov[i][j]*(dPP[k][l]+dPP[l][k]) )
1098 +dgup[k][m][l]*(DD[i][j][l]+DD[j][i][l]-DD[l][i][j])
1099 -dgup[k][m][l]*(g_cov[j][l]*PP[i]+g_cov[i][l]*PP[j]-g_cov[i][j]*PP[l])
1100 -2*g_contr[m][l]*(DD[k][j][l]*PP[i]+DD[k][i][l]*PP[j]-DD[k][i][j]*PP[l]);
1101 }
1102
1103 double Riemann[3][3][3][3] = {0};
1104 for (int i = 0; i < 3; i++)
1105 for (int j = 0; j < 3; j++)
1106 for (int m = 0; m < 3; m++)
1107 for (int k = 0; k < 3; k++){
1108 Riemann[i][k][j][m] = dChristoffel[k][i][j][m] - dChristoffel[j][i][k][m];
1109 for (int l = 0; l < 3; l++){
1110 Riemann[i][k][j][m] += Christoffel[i][j][l]*Christoffel[l][k][m]-Christoffel[i][k][l]*Christoffel[l][j][m];
1111 }
1112 }
1113
1114 double Ricci[3][3] = {0,0,0,0,0,0,0,0,0};
1115 for (int m = 0; m < 3; m++)
1116 for (int n = 0; n < 3; n++)
1117 for (int l = 0; l < 3; l++) Ricci[m][n] += Riemann[m][l][n][l];
1118
1119 double R=0;
1120 for (int i = 0; i < 3; i++)
1121 for (int j = 0; j < 3; j++) R += phi2*g_contr[i][j]*Ricci[i][j];
1122 double Kupdown=0;
1123 for (int i = 0; i < 3; i++)
1124 for (int j = 0; j < 3; j++) Kupdown += Kex[i][j]*Kup[i][j];
1125
1126 double Ham=0;
1127 Ham = R-Kupdown+traceK*traceK;
1128
1129 double dKex[3][3][3]={0};
1130 for (int i = 0; i < 3; i++)
1131 for (int j = 0; j < 3; j++)
1132 for (int k = 0; k < 3; k++) dKex[k][i][j]= (1.0/phi2)*(
1133 dAex[k][i][j]-2*Aex[i][j]*PP[k]+(1.0/3)*dtraceK[k]*g_cov[i][j]
1134 +(2.0/3)*traceK*DD[k][i][j]-(2.0/3)*traceK*g_cov[i][j]*PP[k]);
1135
1136 double Mom[3]={0,0,0};
1137 for (int i = 0; i < 3; i++)
1138 for (int j = 0; j < 3; j++)
1139 for (int l = 0; l < 3; l++){
1140 Mom[i] += phi2*g_contr[j][l]*(dKex[l][i][j]-dKex[i][j][l]);
1141 for (int m = 0; m < 3; m++){
1142 Mom[i] += phi2*g_contr[j][l]*(-Christoffel[j][l][m]*Kex[m][i]+Christoffel[j][i][m]*Kex[m][l]);
1143 }
1144 }
1145
1146 memset(constraints, 0, 6*sizeof(double));
1147 constraints[0] = Ham;
1148 constraints[1] = Mom[0];
1149 constraints[2] = Mom[1];
1150 constraints[3] = Mom[2];
1151 constraints[4] = 1.0/invdet-1.0;
1152 constraints[5] = traceA;
1153}
1154#if defined(GPUOffloadingOMP)
1155#pragma omp end declare target
1156#endif
u
float u
Definition equation_of_state.h:107
temp
const int temp
Definition kernel_hydro.h:264
m
float m
Definition hydro.h:257
examples::exahype2::mgccz4::ncp
void ncp(double *BgradQ, const double *const Q, const double *const gradQSerialised, const int normal, const int MGCCZ4LapseType, const double MGCCZ4ds, const double MGCCZ4c, const double MGCCZ4e, const double MGCCZ4f, const double MGCCZ4bs, const double MGCCZ4sk, const double MGCCZ4xi, const double MGCCZ4mu)
Definition MGCCZ4Kernels.cpp:504
examples::exahype2::mgccz4::source
void source(double *S, const double *const Q, const int MGCCZ4LapseType, const double MGCCZ4ds, const double MGCCZ4c, const double MGCCZ4e, const double MGCCZ4f, const double MGCCZ4bs, const double MGCCZ4sk, const double MGCCZ4xi, const double MGCCZ4itau, const double MGCCZ4eta, const double MGCCZ4k1, const double MGCCZ4k2, const double MGCCZ4k3)
Definition MGCCZ4Kernels.cpp:112
examples::exahype2::mgccz4::admconstraints
void admconstraints(double *constraints, const double *const Q, const double *const gradQSerialised)
Definition MGCCZ4Kernels.cpp:972
examples::exahype2::mgccz4::enforceMGCCZ4constraints
void enforceMGCCZ4constraints(double *luh)
Definition MGCCZ4Kernels.cpp:16
beta
const struct part *restrict const struct part *restrict const float const float beta
Definition signal_velocity.h:45