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