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