Peano
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(WITH_GPU_OMP_TARGET)
14 #pragma omp declare target
15 #endif
16 void 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(WITH_GPU_OMP_TARGET)
106 #pragma omp end declare target
107 #endif
108 
109 #if defined(WITH_GPU_OMP_TARGET)
110 #pragma omp declare target
111 #endif
112 void 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(WITH_GPU_OMP_TARGET)
497 #pragma omp end declare target
498 #endif
499 
500 
501 #if defined(WITH_GPU_OMP_TARGET)
502 #pragma omp declare target
503 #endif
504 void 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(WITH_GPU_OMP_TARGET)
966 #pragma omp end declare target
967 #endif
968 
969 #if defined(WITH_GPU_OMP_TARGET)
970 #pragma omp declare target
971 #endif
972 void 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(WITH_GPU_OMP_TARGET)
1155 #pragma omp end declare target
1156 #endif
euler.u
u
Definition: euler.py:113
euler.j
j
Definition: euler.py:95
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