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