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