Peano 4
Loading...
Searching...
No Matches
CCZ4Kernels.cpph
Go to the documentation of this file.
1#pragma once
2
3#include <algorithm>
4
5
7 const double* const Q,
8 int normal,
9 const double CCZ4e,
10 const double CCZ4ds,
11 const double CCZ4GLMc,
12 const double CCZ4GLMd
13) {
14 const double np = 1.0;
15 const double qmin = std::min({Q[0],Q[3],Q[5]});
16 const double alpha = std::max({1.0, Q[16]}) * std::max({1.0, pow(Q[54],1.0/np)})/std::sqrt(qmin);
17
18 constexpr double sqrtwo = 1.4142135623730951;
19 const double tempA = alpha * std::max({sqrtwo, CCZ4e, CCZ4ds, CCZ4GLMc/alpha, CCZ4GLMd/alpha});
20 const double tempB = Q[17+normal];//DOT_PRODUCT(Q(18:20),nv(:))
21 return std::max({1.0, std::abs(-tempA-tempB), std::abs(tempA-tempB)});
22}
23
24//this version use alpha and phi
25void applications::exahype2::ccz4::source(double* S, const double* const Q,
26 const int CCZ4LapseType,
27 const double CCZ4ds,
28 const double CCZ4c,
29 const double CCZ4e,
30 const double CCZ4f,
31 const double CCZ4bs,
32 const double CCZ4sk,
33 const double CCZ4xi,
34 const double CCZ4itau,
35 const double CCZ4eta,
36 const double CCZ4k1,
37 const double CCZ4k2,
38 const double CCZ4k3,
39 const double CCZ4SO
40 )
41{
42 // Baojiu-01-08-22: alpha now is larger than 0.0001, in agreement with 2112.10567
43 //const double alpha = std::exp(std::fmax(-20., std::fmin(20.,Q[16])));
44 const double alpha = std::fmax(1e-4,Q[16]);
45 //printf("alpha %f\n",alpha);
46 double fa = 1.0;
47 double faa = 0.0;
48 if (CCZ4LapseType==1)
49 {
50 fa = 2./alpha;
51 faa = -fa/alpha;
52 }
53
54 // Note g_cov is symmetric
55 const double g_cov[3][3] = { {Q[0], Q[1], Q[2]}, {Q[1], Q[3], Q[4]}, {Q[2], Q[4], Q[5]} };
56 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];
57 const double invdet = 1./det;
58
59 const double g_contr[3][3] = {
60 { ( 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},
61 {-( 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},
62 {-(-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}
63 };
64
65 // NOTE Aex is symmetric
66 double Aex[3][3] = { {Q[6], Q[7], Q[8]}, {Q[7], Q[9], Q[10]}, {Q[8], Q[10], Q[11]} };
67
68 double traceA = 0;
69 for (int i=0;i<3;i++)
70 for (int j=0;j<3;j++) traceA += g_contr[i][j]*Aex[i][j];
71
72 for (int i=0;i<3;i++)
73 for (int j=0;j<3;j++) Aex[i][j] -= 1./3. * traceA * g_cov[i][j];
74
75 // Matrix multiplications Amix = matmul(g_contr, Aex) Aup = matmul(g_contr, transpose(Amix))
76 double Amix[3][3]={0,0,0,0,0,0,0,0,0};
77 for (int i = 0; i < 3; i++)
78 for (int j = 0; j < 3; j++)
79 for (int u = 0; u < 3; u++) Amix[i][j] += g_contr[i][u] * Aex[u][j];
80 double Aup[3][3]={0,0,0,0,0,0,0,0,0};
81 for (int i = 0; i < 3; i++)
82 for (int j = 0; j < 3; j++)
83 for (int u = 0; u < 3; u++) Aup[i][j] += g_contr[i][u] * Amix[j][u]; // Note the transposition is in the indices
84
85 const double Theta = Q[12]; // needed later
86
87 const double DD[3][3][3] = {
88 {{Q[35], Q[36], Q[37]}, {Q[36], Q[38], Q[39]}, {Q[37], Q[39], Q[40]}},
89 {{Q[41], Q[42], Q[43]}, {Q[42], Q[44], Q[45]}, {Q[43], Q[45], Q[46]}},
90 {{Q[47], Q[48], Q[49]}, {Q[48], Q[50], Q[51]}, {Q[49], Q[51], Q[52]}}
91 };
92
93 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};
94 for (int k = 0; k < 3; k++)
95 for (int m = 0; m < 3; m++)
96 for (int l = 0; l < 3; l++)
97 for (int n = 0; n < 3; n++)
98 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];
99
100 // Baojiu-01-08-22: moved here because phi will be needed in the new Christoffel
101 // Baojiu-01-08-22: now phi^2 is bound to be larger than 0.0001, in agreement with 2112.10567
102 //const double phi = std::exp(std::fmax(-20., std::fmin(20.,Q[54])));
103 const double np = 1.0; // Baojiu-13-06-23: chi=phi^n, np is the index n here
104 const double phi = std::fmax(1e-2,pow(Q[54],(1.0/np))); // Baojiu-13-06-23: changed so that Q[54] is now chi=phi^n
105 const double phi2 = phi*phi;
106 const double chi = pow(phi,np); // Baojiu-13-06-23: new change of variable
107
108 const double PP[3] = {Q[55], Q[56], Q[57]};
109
110 // =====================================================
111 // Baojiu-13-06-23 COMMENTED OUT ALL LINES IN THIS BLOCK
112 // -----------------------------------------------------
113//
114// 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};
115// 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};
116//
117// for (int j = 0; j < 3; j++)
118// for (int i = 0; i < 3; i++)
119// for (int k = 0; k < 3; k++)
120// for (int l = 0; l < 3; l++)
121// {
122// Christoffel_tilde[i][j][k] += g_contr[k][l] * ( DD[i][j][l] + DD[j][i][l] - DD[l][i][j] );
123// // Baojiu-01-08-22: multiplied all PP by 1/phi
124// 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] );
125// }
126//
127//
128// 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};
129// 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};
130// for (int l = 0; l < 3; l++)
131// for (int m = 0; m < 3; m++)
132// for (int j = 0; j < 3; j++)
133// for (int i = 0; i < 3; i++)
134// for (int k = 0; k < 3; k++)
135// {
136// // Baojiu-01-08-22: added PP*PP terms, and multiplied all PP terms by 1/phi
137// dChristoffelSrc[k][i][j][m] += dgup[k][m][l]*(DD[i][j][l] +DD[j][i][l] -DD[l][i][j] )
138// - dgup[k][m][l]*(g_cov[j][l]*PP[i]+g_cov[i][l]*PP[j]-g_cov[i][j]*PP[l])/phi
139// -2.0*g_contr[m][l]*(DD[k][j][l]*PP[i]+DD[k][i][l]*PP[j]-DD[k][i][j]*PP[l])/phi
140// + 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;
141//
142// dChristoffel_tildeSrc[k][i][j][m] += dgup[k][m][l]*(DD[i][j][l]+DD[j][i][l]-DD[l][i][j]);
143// }
144//
145// 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};
146// for (int m = 0; m < 3; m++)
147// for (int j = 0; j < 3; j++)
148// for (int k = 0; k < 3; k++)
149// for (int i = 0; i < 3; i++)
150// {
151// RiemannSrc[i][k][j][m] = dChristoffelSrc[k][i][j][m] - dChristoffelSrc[j][i][k][m];
152// 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];
153// }
154//
155// double RicciSrc[3][3] = {0,0,0,0,0,0,0,0,0};
156// for (int l = 0; l < 3; l++)
157// for (int n = 0; n < 3; n++)
158// for (int m = 0; m < 3; m++) RicciSrc[m][n] += RiemannSrc[m][l][n][l];
159//
160// // Here we directly compute the derivative of Gtilde from its original definition as contracted Christoffel symbol,
161// // without assuming unit determinant of the conformal metric. Back to the roots, and as few assumptions as possible...
162// double dGtildeSrc[3][3] = {0,0,0,0,0,0,0,0,0};
163// for (int l = 0; l < 3; l++)
164// for (int j = 0; j < 3; j++)
165// for (int i = 0; i < 3; i++)
166// 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];
167//
168// double dZSrc[3][3] = {0,0,0,0,0,0,0,0,0};
169// for (int j = 0; j < 3; j++)
170// for (int i = 0; i < 3; i++)
171// for (int k = 0; k < 3; k++) dZSrc[k][i] += CCZ4ds*(DD[k][i][j]*(Ghat[j]-Gtilde[j]) - 0.5*g_cov[i][j]*dGtildeSrc[k][j]);
172//
173// double nablaZSrc[3][3] = {0,0,0,0,0,0,0,0,0};
174// for (int j = 0; j < 3; j++)
175// for (int i = 0; i < 3; i++)
176// {
177// nablaZSrc[i][j] = dZSrc[i][j];
178// for (int k = 0; k < 3; k++) nablaZSrc[i][j] -= Christoffel[i][j][k]*Z[k];
179// }
180//
181// double RicciPlusNablaZSrc[3][3];
182// for (int i = 0; i < 3; i++)
183// for (int j = 0; j < 3; j++) RicciPlusNablaZSrc[i][j] = RicciSrc[i][j] + nablaZSrc[i][j] + nablaZSrc[j][i];
184//
185 // -----------------------------------------------------
186 // Baojiu-13-06-23 COMMENTED OUT ALL LINES IN THIS BLOCK
187 // =====================================================
188
189 // ================================
190 // Baojiu-21-09-22 START OF CHANGES
191 // --------------------------------
192
193 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};
194 for (int j = 0; j < 3; j++)
195 for (int i = 0; i < 3; i++)
196 for (int k = 0; k < 3; k++)
197 for (int l = 0; l < 3; l++) Christoffel_tilde[i][j][k] += g_contr[k][l] * ( DD[i][j][l] + DD[j][i][l] - DD[l][i][j] );
198
199 double Christoffel_tilde_down[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};
200 for (int i = 0; i < 3; i++)
201 for (int j = 0; j < 3; j++)
202 for (int k = 0; k < 3; k++)
203 for (int l = 0; l < 3; l++) Christoffel_tilde_down[i][j][k] += g_cov[i][l]*Christoffel_tilde[j][k][l];
204
205 double Gtilde[3] = {0,0,0};
206 for (int l = 0; l < 3; l++)
207 for (int j = 0; j < 3; j++)
208 for (int i = 0; i < 3; i++) Gtilde[i] += g_contr[j][l] * Christoffel_tilde[j][l][i];
209
210 const double Ghat[3] = {Q[13], Q[14], Q[15]};
211
212 double Z[3] = {0,0,0}; // Matrix vector multiplications
213 for (int i=0;i<3;i++)
214 for (int j=0;j<3;j++) Z[i] += 0.5*CCZ4ds*( g_cov[i][j]* (Ghat[j] - Gtilde[j]));
215 double Zup[3] = {0,0,0};
216 for (int i=0;i<3;i++)
217 for (int j=0;j<3;j++) Zup[i] += phi2 * g_contr[i][j] * Z[j];
218
219 double RicciPlusNablaZSrc[3][3] = {0,0,0,0,0,0,0,0,0};
220 for (int i = 0; i < 3; i++)
221 for (int j = 0; j < 3; j++)
222 {
223 RicciPlusNablaZSrc[i][j] -= (np-1.0)*PP[i]*PP[j]/(np*np*chi*chi); // Baojiu-13-06-23: changed to chi=phi^n
224 for (int k = 0; k < 3; k++)
225 {
226 RicciPlusNablaZSrc[i][j] += Ghat[k]*DD[k][i][j];
227 RicciPlusNablaZSrc[i][j] -= PP[k]*Christoffel_tilde[i][j][k]/(np*chi); // Baojiu-13-06-23: changed to chi=phi^n
228 RicciPlusNablaZSrc[i][j] += 2.0*Zup[k]/(np*chi*phi2)*(g_cov[i][k]*PP[j]+g_cov[j][k]*PP[i]-g_cov[i][j]*PP[k]); // Baojiu-13-06-23: changed to chi=phi^n
229 for (int l = 0; l < 3; l++)
230 {
231 RicciPlusNablaZSrc[i][j] -= (np+1.0)*g_cov[i][j]*g_contr[k][l]*PP[k]*PP[l]/(np*np*chi*chi); // Baojiu-13-06-23: changed to chi=phi^n
232 for (int m = 0; m < 3; m++)
233 {
234 RicciPlusNablaZSrc[i][j] += g_contr[l][m]*( Christoffel_tilde[l][i][k] * Christoffel_tilde_down[j][k][m] +
235 Christoffel_tilde[l][j][k] * Christoffel_tilde_down[i][k][m] +
236 Christoffel_tilde[i][m][k] * Christoffel_tilde_down[k][l][j] );
237 RicciPlusNablaZSrc[i][j] -= g_cov[i][j]*g_contr[l][m]*Christoffel_tilde[l][m][k]*PP[k]/(np*chi); // Baojiu-13-06-23: changed to chi=phi^n
238 }
239 }
240 }
241 }
242 // ------------------------------
243 // Baojiu-21-09-22 END OF CHANGES
244 // ==============================
245
246 double RPlusTwoNablaZSrc = 0;
247 for (int i = 0; i < 3; i++)
248 for (int j = 0; j < 3; j++) RPlusTwoNablaZSrc += g_contr[i][j]*RicciPlusNablaZSrc[i][j];
249 RPlusTwoNablaZSrc*=phi2;
250
251 const double AA[3] = {Q[23], Q[24], Q[25]};
252
253 double nablaijalphaSrc[3][3] = {0,0,0,0,0,0,0,0,0};
254 for (int j = 0; j < 3; j++)
255 for (int i = 0; i < 3; i++)
256 {
257 // ================================
258 // Baojiu-21-09-22 START OF CHANGES
259 // --------------------------------
260// for (int k = 0; k < 3; k++) nablaijalphaSrc[i][j] -= Christoffel[i][j][k]*AA[k]; // This is the original line; all other lines are new
261 nablaijalphaSrc[i][j] += (PP[i]*AA[j] + PP[j]*AA[i])/(np*chi); // Baojiu-13-06-23: changed to chi=phi^n
262 for (int k = 0; k < 3; k++)
263 {
264 nablaijalphaSrc[i][j] -= Christoffel_tilde[i][j][k]*AA[k];
265 for (int l = 0; l < 3; l++)
266 {
267 nablaijalphaSrc[i][j] -= g_cov[i][j]*g_contr[k][l]*PP[k]*AA[l]/(np*chi); // Baojiu-13-06-23: changed to chi=phi^n
268 }
269 }
270 // --------------------------------
271 // Baojiu-21-09-22 END OF CHANGES
272 // ================================
273 }
274
275 double nablanablaalphaSrc=0;
276 for (int i = 0; i < 3; i++)
277 for (int j = 0; j < 3; j++) nablanablaalphaSrc += g_contr[i][j]*nablaijalphaSrc[i][j];
278 nablanablaalphaSrc *= phi2;
279
280 // ================================
281 // Baojiu-21-09-22 START OF CHANGES
282 // --------------------------------
283 // these lines won't make a difference
284// nablanablaalphaSrc = 0.0;
285// for (int i = 0; i < 3; i++)
286// {
287// nablanablaalphaSrc -= phi2*Gtilde[i]*AA[i];
288// for (int j = 0; j < 3; j++) nablanablaalphaSrc -= phi2*g_contr[i][j]*PP[i]*AA[j]/(np*chi); // Baojiu-13-06-23: changed to chi=phi^n
289// }
290 // --------------------------------
291 // Baojiu-21-09-22 END OF CHANGES
292 // ================================
293
294 double SecondOrderTermsSrc[3][3];
295 for (int i = 0; i < 3; i++)
296 for (int j = 0; j < 3; j++) SecondOrderTermsSrc[i][j] = -nablaijalphaSrc[i][j] + alpha*RicciPlusNablaZSrc[i][j];
297
298 double traceSrc = 0;
299 for (int i = 0; i < 3; i++)
300 for (int j = 0; j < 3; j++) traceSrc += g_contr[i][j]*SecondOrderTermsSrc[i][j];
301
302 for (int i = 0; i < 3; i++)
303 for (int j = 0; j < 3; j++) SecondOrderTermsSrc[i][j] -= 1./3 * traceSrc * g_cov[i][j];
304
305 double dtgamma[3][3];
306 for (int i = 0; i < 3; i++)
307 for (int j = 0; j < 3; j++) dtgamma[i][j] = -2.0 * alpha * Aex[i][j] - CCZ4itau*(det -1.0)*g_cov[i][j];
308
309 const double BB[3][3] = {
310 {Q[26], Q[27], Q[28]}, {Q[29], Q[30], Q[31]}, {Q[32], Q[33], Q[34]}
311 };
312
313 double traceB = BB[0][0] + BB[1][1] + BB[2][2];
314 const double beta[3] = {Q[17], Q[18], Q[19]};
315
316 for (int k = 0; k < 3; k++)
317 for (int j = 0; j < 3; j++)
318 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];
319
320 // MATMUL(Aex,Amix)
321 double Atemp[3][3]={0,0,0,0,0,0,0,0,0};
322 for (int i = 0; i < 3; i++)
323 for (int j = 0; j < 3; j++)
324 for (int u = 0; u < 3; u++) Atemp[i][j] += Aex[i][u] * Amix[u][j];
325
326 const double traceK = Q[53];
327
329 double dtK[3][3];
330 for (int i = 0; i < 3; i++)
331 for (int j = 0; j < 3; j++) dtK[i][j] = phi2*SecondOrderTermsSrc[i][j] + alpha*Aex[i][j]*(traceK-2.*Theta*CCZ4c) - 2.*alpha*Atemp[i][j] - CCZ4itau*g_cov[i][j]*traceA;
332
333 for (int j = 0; j < 3; j++)
334 for (int i = 0; i < 3; i++)
335 for (int k = 0; k < 3; k++) dtK[i][j] += Aex[k][i]*BB[j][k] + Aex[k][j]*BB[i][k] - 2./3.*Aex[i][j]*BB[k][k];
336
337 const double K0 = Q[58];
338
339 // Baojiu-01-08-22: alpha*k1 -> k1 below
340 double dtTraceK = -nablanablaalphaSrc + alpha*(RPlusTwoNablaZSrc + traceK*traceK - 2.0*CCZ4c*Theta*traceK) - 3.0*alpha*CCZ4k1*(1.+CCZ4k2)*Theta/alpha;
341
342 // Baojiu-01-08-22: multiplied the second term on the RHS by phi in the equation below
343 double dtchi = beta[0]*PP[0] + beta[1]*PP[1] + beta[2]*PP[2] + np/3.*(alpha*traceK-traceB)*chi; // Baojiu-13-06-23: changed to chi=phi^n
344
345 // Baojiu-01-08-22: multiplied the first term on the RHS by alpha in the equation below
346 double dtalpha = -alpha*alpha*fa*(traceK-K0-2.*CCZ4c*Theta)+beta[0]*AA[0]+beta[1]*AA[1]+beta[2]*AA[2];
347
348 double Aupdown = 0;
349 for (int i = 0; i < 3; i++)
350 for (int j = 0; j < 3; j++) Aupdown += Aex[i][j]*Aup[i][j];
351
352 double sumzupaa = 0.0;
353 for (int i = 0; i < 3; i++) sumzupaa += Zup[i]*AA[i];
354 // Baojiu-01-08-22: divided sumzupaa by alpha, and rescaled alpha*k1 -> k1, in the equation below
355 double dtTheta = 0.5*alpha*CCZ4e*CCZ4e*(RPlusTwoNablaZSrc - Aupdown + 2./3.*traceK*traceK) - alpha*(CCZ4c*Theta*traceK + sumzupaa/alpha + CCZ4k1*(2.+CCZ4k2)*Theta/alpha);
356
357 double dtGhat[3] = {0,0,0};
358 for (int i = 0; i < 3; i++)
359 {
360 double temp1=0, temp2=0, temp3=0, temp4=0, temp5=0, temp6=0;
361 for (int m = 0; m < 3; m++)
362 {
363 temp1 += Aup[i][m]*PP[m]/(np*chi); // Baojiu-13-06-23: changed to chi=phi^n
364 temp3 += Aup[i][m]*AA[m]/alpha; // Baojiu-01-08-22: divided AA by alpha
365 temp2 += g_contr[m][i]*(-Theta*AA[m]/alpha -2./3.*traceK*Z[m]); // Baojiu-01-08-22: divided AA by alpha
366 temp4 += g_contr[i][m]*Z[m];
367 temp5 += Gtilde[m]*BB[m][i];
368 for (int n = 0; n < 3; n++) temp6 += Christoffel_tilde[m][n][i]*Aup[m][n];
369 }
370 dtGhat[i] += 2.*alpha*(temp6 - 3.*temp1 + temp2 - temp3 - CCZ4k1*temp4/alpha) - temp5 + 2./3.*Gtilde[i]*traceB; // Baojiu-01-08-22: rescaled alpha*k1 -> k1
371
372 for (int l = 0; l < 3; l++)
373 for (int k = 0; k < 3; k++) dtGhat[i] += 2.*CCZ4k3*(2./3.*g_contr[i][l]*Z[l]*BB[k][k] - g_contr[l][k]*Z[l]*BB[k][i]);
374 }
375
376 double ov[3];
377 for (int k = 0; k < 3; k++)
378 {
379 double temp=0;
380 for (int i = 0; i < 3; i++)
381 for (int j = 0; j < 3; j++) temp += dgup[k][i][j]*Aex[i][j];
382 ov[k] = 2*alpha*temp;
383 }
384
385 // matrix vector multiplication in a loop and add result to existing vector
386 for (int i = 0; i < 3; i++)
387 for (int j = 0; j < 3; j++) dtGhat[i] += CCZ4sk*g_contr[i][j]*ov[j];
388
389 const double myb[3] = {Q[20], Q[21], Q[22]};
390
391 double dtbb[3];
392 for (int i = 0; i < 3; i++) dtbb[i] = CCZ4sk*(CCZ4xi*dtGhat[i] - CCZ4eta*myb[i]);
393
394 double dtbeta[3];
395 for (int i = 0; i < 3; i++) dtbeta[i] = CCZ4f*myb[i];
396 for (int i = 0; i < 3; i++) dtbeta[i] += CCZ4bs*(beta[0]*BB[0][i] + beta[1]*BB[1][i] + beta[2]*BB[2][i]);
397 for (int i = 0; i < 3; i++) dtbeta[i] *= CCZ4sk;
398
399 // Auxiliary variables
400 double dtA[3] ={0,0,0};
401 for (int i = 0; i < 3; i++)
402 {
403 // Baojiu-01-08-22: in the equation below: fa -> 2*fa
404 dtA[i] = -alpha*AA[i]*(2.*fa+alpha*faa)*(traceK - K0 - 2.*CCZ4c*Theta);
405 for (int j = 0; j < 3; j++) dtA[i] += BB[i][j]*AA[j];
406 }
407
408 for (int k = 0; k < 3; k++)
409 {
410 double temp = 0;
411 for (int i = 0; i < 3; i++)
412 for (int j = 0; j < 3; j++) temp+= dgup[k][i][j]*Aex[i][j];
413 dtA[k] += -CCZ4sk*alpha*alpha*fa*temp; // Baojiu-13-06-23: here there is only one alpha as in the Dumbser paper
414 }
415
416 double dtB[3][3] = {0,0,0,0,0,0,0,0,0};
417 for (int i = 0; i < 3; i++)
418 for (int j = 0; j < 3; j++)
419 for (int u = 0; u < 3; u++) dtB[i][j] += CCZ4sk*(BB[i][u] * BB[u][j]);
420
421 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};
422 for (int m = 0; m < 3; m++)
423 for (int j = 0; j < 3; j++)
424 for (int i = 0; i < 3; i++)
425 for (int k = 0; k < 3; k++)
426 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];
427
428 for (int j = 0; j < 3; j++)
429 for (int i = 0; i < 3; i++)
430 for (int k = 0; k < 3; k++)
431 {
432 dtD[k][i][j] -= AA[k]*Aex[i][j]; // Baojiu-01-08-22: alpha*AA -> alpha
433 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];
434 }
435
436 double dtP[3] = {0,0,0};
437 for (int i = 0; i < 3; i++)
438 for (int j = 0; j < 3; j++) dtP[i] += BB[i][j]*PP[j];
439
440 for (int k = 0; k < 3; k++)
441 {
442 double temp=0;
443 for (int i = 0; i < 3; i++)
444 for (int j = 0; j < 3; j++) temp += dgup[k][i][j]*Aex[i][j];
445 dtP[k] += np/3.*AA[k]*traceK*chi; // Baojiu-13-06-23: changed to chi=phi^n
446 dtP[k] += np/3.*(alpha*traceK-traceB)*PP[k]; // Baojiu-13-06-23: changed to chi=phi^n
447 dtP[k] += np/3.*alpha*CCZ4sk*temp*chi; // Baojiu-13-06-23: changed to chi=phi^n
448// dtP[k] += np/3.*alpha*CCZ4sk*temp; // Baojiu-13-06-23: changed to chi=phi^n // Baojiu-27-06-23
449 }
450
451 S[0] = dtgamma[0][0];
452 S[1] = dtgamma[0][1];
453 S[2] = dtgamma[0][2];
454 S[3] = dtgamma[1][1];
455 S[4] = dtgamma[1][2];
456 S[5] = dtgamma[2][2];
457 S[6] = dtK[0][0];
458 S[7] = dtK[0][1];
459 S[8] = dtK[0][2];
460 S[9] = dtK[1][1];
461 S[10] = dtK[1][2];
462 S[11] = dtK[2][2];
463 //S[12] = 0.0;
464 S[12] = dtTheta;
465 for (int i = 0; i < 3; i++) S[13+i] = dtGhat[i];
466 S[16] = dtalpha;
467 for (int i = 0; i < 3; i++) S[17+i] = dtbeta[i];
468 for (int i = 0; i < 3; i++) S[20+i] = dtbb[i];
469 for (int i = 0; i < 3; i++) S[23+i] = dtA[i];
470 S[26] = dtB[0][0];
471 S[27] = dtB[0][1];
472 S[28] = dtB[0][2];
473 S[29] = dtB[1][0];
474 S[30] = dtB[1][1];
475 S[31] = dtB[1][2];
476 S[32] = dtB[2][0];
477 S[33] = dtB[2][1];
478 S[34] = dtB[2][2];
479 S[35] = dtD[0][0][0];
480 S[36] = dtD[0][0][1];
481 S[37] = dtD[0][0][2];
482 S[38] = dtD[0][1][1];
483 S[39] = dtD[0][1][2];
484 S[40] = dtD[0][2][2];
485 S[41] = dtD[1][0][0];
486 S[42] = dtD[1][0][1];
487 S[43] = dtD[1][0][2];
488 S[44] = dtD[1][1][1];
489 S[45] = dtD[1][1][2];
490 S[46] = dtD[1][2][2];
491 S[47] = dtD[2][0][0];
492 S[48] = dtD[2][0][1];
493 S[49] = dtD[2][0][2];
494 S[50] = dtD[2][1][1];
495 S[51] = dtD[2][1][2];
496 S[52] = dtD[2][2][2];
497 S[53] = dtTraceK;
498 S[54] = dtchi;
499 for (int i = 0; i < 3; i++) S[55+i] = dtP[i];
500 S[58] = 0;
501 //used only in soccz4
502 if (CCZ4SO==1){
503 for (int i = 23; i < 53; i++) S[i] = 0.0;
504 for (int i = 0; i < 3; i++) S[55+i] = 0.0;
505 }
506}
507
508
509void applications::exahype2::ccz4::ncp(double* BgradQ, const double* const Q, const double* const gradQSerialised, const int normal,
510 const int CCZ4LapseType,
511 const double CCZ4ds,
512 const double CCZ4c,
513 const double CCZ4e,
514 const double CCZ4f,
515 const double CCZ4bs,
516 const double CCZ4sk,
517 const double CCZ4xi,
518 const double CCZ4mu,
519 const double CCZ4SO
520 )
521{
522 // Han 09-2022: add a parameter to switch off or on the advection here
523 const double advec=0.0;
524 const double advec_aux=0.0;
525
526// const double alpha = std::exp(std::fmax(-20., std::fmin(20.,Q[16])));
527 const double alpha = std::fmax(1e-4,Q[16]); // Baojiu-01-08-22: changed to make sure alpha>0.0001
528
529 //printf("alpha %f\n",alpha);
530 const double alpha2 = alpha*alpha;
531 double fa = 1.0;
532 double faa = 0.0;
533 if (CCZ4LapseType==1)
534 {
535 fa = 2./alpha;
536 faa = -fa/alpha;
537 }
538
539 constexpr int nVar(59);
540
541 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};
542
543 // De-serialise input data and fill static array
544 // FIXME the use of 2D arrays can be avoided: all terms not in the normal are 0
545 for (int i=0; i<nVar; i++) gradQin[i][normal] = gradQSerialised[i+normal*nVar];
546
547 // Note g_cov is symmetric
548 const double g_cov[3][3] = { {Q[0], Q[1], Q[2]}, {Q[1], Q[3], Q[4]}, {Q[2], Q[4], Q[5]} };
549 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]);
550
551 const double g_contr[3][3] = {
552 { ( 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},
553 {-( 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},
554 {-(-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}
555 };
556
557 // NOTE Aex is symmetric
558 double Aex[3][3] = { {Q[6], Q[7], Q[8]}, {Q[7], Q[9], Q[10]}, {Q[8], Q[10], Q[11]} };
559
560 double traceA = 0;
561 for (int i=0;i<3;i++)
562 for (int j=0;j<3;j++) traceA+=g_contr[i][j]*Aex[i][j];
563 //traceA *= 1./3;
564
565 for (int i=0;i<3;i++)
566 for (int j=0;j<3;j++) Aex[i][j] -= 1./3. * traceA * g_cov[i][j];
567
568 // Matrix multiplications Amix = matmul(g_contr, Aex) Aup = matmul(g_contr, mytranspose(Amix))
569 double Amix[3][3]={0,0,0,0,0,0,0,0,0};
570 for (int i = 0; i < 3; i++)
571 for (int j = 0; j < 3; j++)
572 for (int u = 0; u < 3; u++) Amix[i][j] += g_contr[i][u] * Aex[u][j];
573 double Aup[3][3]={0,0,0,0,0,0,0,0,0};
574 for (int i = 0; i < 3; i++)
575 for (int j = 0; j < 3; j++)
576 for (int u = 0; u < 3; u++) Aup[i][j] += g_contr[i][u] * Amix[j][u]; // Note the transposition is in the indices
577
578 const double DD[3][3][3] = {
579 {{Q[35], Q[36], Q[37]}, {Q[36], Q[38], Q[39]}, {Q[37], Q[39], Q[40]}},
580 {{Q[41], Q[42], Q[43]}, {Q[42], Q[44], Q[45]}, {Q[43], Q[45], Q[46]}},
581 {{Q[47], Q[48], Q[49]}, {Q[48], Q[50], Q[51]}, {Q[49], Q[51], Q[52]}}
582 };
583 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};
584 for (int k = 0; k < 3; k++)
585 for (int m = 0; m < 3; m++)
586 for (int l = 0; l < 3; l++)
587 for (int n = 0; n < 3; n++)
588 for (int j = 0; j < 3; j++) dgup[k][m][l] -= g_contr[m][n]*g_contr[j][l]*2*DD[k][n][j];
589
590 // Baojiu-01-08-22: changed so that phi^2>0.0001
591// const double phi = std::exp(std::fmax(-20., std::fmin(20.,Q[54])));
592 const double np = 1.0; // Baojiu-13-06-23: changed to chi=phi^n
593 const double phi = std::fmax(1e-2,pow(Q[54],1.0/np)); // Baojiu-13-06-23: changed to chi=phi^n, note Q[54] is chi
594 const double chi = pow(phi,np); // Baojiu-13-06-23: changed to chi=phi^n
595 const double phi2 = phi*phi;
596
597 const double PP[3] = {Q[55], Q[56], Q[57]};
598 double dPP[3][3] = {
599 {gradQin[55][0],gradQin[56][0],gradQin[57][0]},
600 {gradQin[55][1],gradQin[56][1],gradQin[57][1]},
601 {gradQin[55][2],gradQin[56][2],gradQin[57][2]}
602 };
603
604 const double AA[3] = {Q[23], Q[24], Q[25]};
605 double dAA[3][3] = {
606 {gradQin[23][0],gradQin[24][0],gradQin[25][0]},
607 {gradQin[23][1],gradQin[24][1],gradQin[25][1]},
608 {gradQin[23][2],gradQin[24][2],gradQin[25][2]}
609 };
610
611// 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};
612 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};
613// 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};
614
615 for (int j = 0; j < 3; j++)
616 for (int i = 0; i < 3; i++)
617 for (int k = 0; k < 3; k++)
618 {
619// Christoffel_kind1[i][j][k] = DD[k][i][j] + DD[j][i][k] - DD[i][j][k];
620 for (int l = 0; l < 3; l++)
621 {
622 Christoffel_tilde[i][j][k] += g_contr[k][l] * ( DD[i][j][l] + DD[j][i][l] - DD[l][i][j] );
623// 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] )/phi;
624 }
625 }
626
627 double Gtilde[3] = {0,0,0};
628 for (int l = 0; l < 3; l++)
629 for (int j = 0; j < 3; j++)
630 for (int i = 0; i < 3; i++) Gtilde[i] += g_contr[j][l] * Christoffel_tilde[j][l][i];
631
632 const double Ghat[3] = {Q[13], Q[14], Q[15]};
633
634 double Z[3] = {0,0,0};
635 for (int i=0;i<3;i++)
636 for (int j=0;j<3;j++) Z[i] += 0.5*( g_cov[i][j]* (Ghat[j] - Gtilde[j]));
637 double Zup[3] = {0,0,0};
638 for (int i=0;i<3;i++)
639 for (int j=0;j<3;j++) Zup[i] += phi2 * g_contr[i][j] * Z[j];
640
641 const double dDD[3][3][3][3] = {
642 {
643 {
644 {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]},
645 },
646 {
647 {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]},
648 },
649 {
650 {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]}
651 }
652 },
653 {
654 {
655 {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]},
656 },
657 {
658 {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]},
659 },
660 {
661 {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]}
662 }
663 },
664 {
665 {
666 {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]},
667 },
668 {
669 {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]},
670 },
671 {
672 {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]}
673 }
674 }
675 };
676
677 // =====================================================
678 // Baojiu-13-06-23 COMMENTED OUT ALL LINES IN THIS BLOCK
679 // -----------------------------------------------------
680//
681// 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};
682// 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};
683// for (int i = 0; i < 3; i++)
684// for (int j = 0; j < 3; j++)
685// for (int m = 0; m < 3; m++)
686// for (int k = 0; k < 3; k++)
687// {
688// dChristoffelNCP [k][i][j][m] = 0;
689// dChristoffel_tildeNCP[k][i][j][m] = 0;
690// for (int l = 0; l < 3; l++)
691// {
692// // Baojiu-01-08-22: multiplied the original dPP terms by 1/phi //a small issue
693// dChristoffelNCP[k][i][j][m] += 0.5*g_contr[m][l] * (
694// 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]
695// - 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 );
696//
697// 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]);
698//
699// }
700// }
701//
702// double RiemannNCP[3][3][3][3] = {0};
703// for (int i = 0; i < 3; i++)
704// for (int j = 0; j < 3; j++)
705// for (int m = 0; m < 3; m++)
706// for (int k = 0; k < 3; k++) RiemannNCP[i][k][j][m] = dChristoffelNCP[k][i][j][m] - dChristoffelNCP[j][i][k][m];
707//
708// double RicciNCP[3][3] = {0,0,0,0,0,0,0,0,0};
709// for (int m = 0; m < 3; m++)
710// for (int n = 0; n < 3; n++)
711// for (int l = 0; l < 3; l++) RicciNCP[m][n] += RiemannNCP[m][l][n][l];
712//
713// double dGtildeNCP[3][3] = {0,0,0,0,0,0,0,0,0};
714// for (int i = 0; i < 3; i++)
715// for (int k = 0; k < 3; k++)
716// for (int j = 0; j < 3; j++)
717// for (int l = 0; l < 3; l++) dGtildeNCP[k][i] += g_contr[j][l]*dChristoffel_tildeNCP[k][j][l][i];
718//
719// double dZNCP[3][3] = {0,0,0,0,0,0,0,0,0};
720// for (int j = 0; j < 3; j++)
721// for (int i = 0; i < 3; i++)
722// for (int k = 0; k < 3; k++) dZNCP[k][i] += CCZ4ds*0.5*g_cov[i][j]*(dGhat[k][j]-dGtildeNCP[k][j]);
723//
724// double RicciPlusNablaZNCP[3][3];
725// for (int i = 0; i < 3; i++)
726// for (int j = 0; j < 3; j++) RicciPlusNablaZNCP[i][j] = RicciNCP[i][j] + dZNCP[i][j] + dZNCP[j][i];
727//
728 // -----------------------------------------------------
729 // Baojiu-13-06-23 COMMENTED OUT ALL LINES IN THIS BLOCK
730 // =====================================================
731
732 const double dGhat[3][3] = {
733 {gradQin[13][0],gradQin[14][0],gradQin[15][0]},
734 {gradQin[13][1],gradQin[14][1],gradQin[15][1]},
735 {gradQin[13][2],gradQin[14][2],gradQin[15][2]}
736 };
737
738 // ================================
739 // Baojiu-21-09-22 START OF CHANGES
740 // --------------------------------
741 double RicciPlusNablaZNCP[3][3] = {0,0,0,0,0,0,0,0,0};
742 for (int i = 0; i < 3; i++)
743 for (int j = 0; j < 3; j++)
744 {
745 RicciPlusNablaZNCP[i][j] += 0.5*(dPP[i][j] + dPP[j][i])/(np*chi); // Baojiu-13-06-23: changed to chi=phi^n
746 for (int k = 0; k < 3; k++)
747 {
748 RicciPlusNablaZNCP[i][j] += 0.5*( g_cov[k][i]*dGhat[j][k] + g_cov[k][j]*dGhat[i][k] );
749 for (int l = 0; l < 3; l++)
750 {
751 RicciPlusNablaZNCP[i][j] -= 0.5*g_contr[k][l]*( dDD[k][l][i][j] + dDD[l][k][i][j] ); // explicit symmetrisation
752 RicciPlusNablaZNCP[i][j] += 0.5*g_cov[i][j]*g_contr[k][l]*( dPP[k][l] + dPP[l][k] )/(np*chi); // Baojiu-13-06-23: changed to chi=phi^n
753 }
754 }
755 }
756 // --------------------------------
757 // Baojiu-21-09-22 END OF CHANGES
758 // ================================
759
760 double RPlusTwoNablaZNCP = 0;
761 for (int i = 0; i < 3; i++)
762 for (int j = 0; j < 3; j++) RPlusTwoNablaZNCP += g_contr[i][j]*RicciPlusNablaZNCP[i][j];
763 RPlusTwoNablaZNCP*=phi2;
764
765 double nablaijalphaNCP[3][3];
766 for (int i = 0; i < 3; i++)
767 for (int j = 0; j < 3; j++) nablaijalphaNCP[i][j] = 0.5*( dAA[i][j] + dAA[j][i] ); // Baojiu-01-08-22: removed the alpha
768
769 double nablanablaalphaNCP = 0;
770 for (int i = 0; i < 3; i++)
771 for (int j = 0; j < 3; j++) nablanablaalphaNCP += g_contr[i][j]*nablaijalphaNCP[i][j];
772 nablanablaalphaNCP*=phi2;
773
774 double SecondOrderTermsNCP[3][3];
775 for (int i = 0; i < 3; i++)
776 for (int j = 0; j < 3; j++) SecondOrderTermsNCP[i][j] = -nablaijalphaNCP[i][j] + alpha*RicciPlusNablaZNCP[i][j];
777
778 double traceNCP = 0;
779 for (int i = 0; i < 3; i++)
780 for (int j = 0; j < 3; j++) traceNCP += g_contr[i][j]*SecondOrderTermsNCP[i][j];
781
782 for (int i = 0; i < 3; i++)
783 for (int j = 0; j < 3; j++) SecondOrderTermsNCP[i][j] -= 1./3 * traceNCP * g_cov[i][j];
784
785 const double beta[3] = {Q[17], Q[18], Q[19]};
786
787 const double dAex[3][3][3] = {
788 {{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]}},
789 {{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]}},
790 {{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]}}
791 };
792
793 double dtK[3][3];
794 for (int i = 0; i < 3; i++)
795 for (int j = 0; j < 3; j++) dtK[i][j] = phi2*SecondOrderTermsNCP[i][j] + advec*(beta[0]*dAex[0][i][j] + beta[1]*dAex[1][i][j] + beta[2]*dAex[2][i][j]); // extrinsic curvature
796
797 const double dtraceK[3] = {gradQin[53][0], gradQin[53][1], gradQin[53][2]};
798
799 double dtTraceK = -nablanablaalphaNCP + alpha*RPlusTwoNablaZNCP + advec*(beta[0]*dtraceK[0] + beta[1]*dtraceK[1] + beta[2]*dtraceK[2]);
800
801 const double BB[3][3] = {
802 {Q[26], Q[27], Q[28]}, {Q[29], Q[30], Q[31]}, {Q[32], Q[33], Q[34]}
803 };
804
805 double traceB = BB[0][0] + BB[1][1] + BB[2][2]; // TODO direct from Q!
806
807 double Aupdown = 0;
808 for (int i = 0; i < 3; i++)
809 for (int j = 0; j < 3; j++) Aupdown += Aex[i][j]*Aup[i][j];
810
811 const double dTheta[3] = {gradQin[12][0],gradQin[12][1],gradQin[12][2]};
812 double dtTheta = 0.5*alpha*CCZ4e*CCZ4e*RPlusTwoNablaZNCP + advec*(beta[0]*dTheta[0] + beta[1]*dTheta[1] + beta[2]*dTheta[2]);
813
814// double divAupNCP[3] = {0,0,0};
815// for (int i = 0; i < 3; i++)
816// for (int j = 0; j < 3; j++)
817// for (int l = 0; l < 3; l++)
818// for (int k = 0; k < 3; k++) divAupNCP[i] += g_contr[i][l]*g_contr[j][k]*dAex[j][l][k];
819
820 const double dBB[3][3][3] = {
821 {
822 {CCZ4sk*gradQin[26][0],CCZ4sk*gradQin[27][0],CCZ4sk*gradQin[28][0]}, {CCZ4sk*gradQin[29][0],CCZ4sk*gradQin[30][0],CCZ4sk*gradQin[31][0]}, {CCZ4sk*gradQin[32][0],CCZ4sk*gradQin[33][0],CCZ4sk*gradQin[34][0]}
823 },
824 {
825 {CCZ4sk*gradQin[26][1],CCZ4sk*gradQin[27][1],CCZ4sk*gradQin[28][1]}, {CCZ4sk*gradQin[29][1],CCZ4sk*gradQin[30][1],CCZ4sk*gradQin[31][1]}, {CCZ4sk*gradQin[32][1],CCZ4sk*gradQin[33][1],CCZ4sk*gradQin[34][1]}
826 },
827 {
828 {CCZ4sk*gradQin[26][2],CCZ4sk*gradQin[27][2],CCZ4sk*gradQin[28][2]}, {CCZ4sk*gradQin[29][2],CCZ4sk*gradQin[30][2],CCZ4sk*gradQin[31][2]}, {CCZ4sk*gradQin[32][2],CCZ4sk*gradQin[33][2],CCZ4sk*gradQin[34][2]}
829 }
830 };
831
832 double dtGhat[3] = {0,0,0};
833 for (int i = 0; i < 3; i++)
834 {
835 double temp=0, temp2=0;
836 for (int j = 0; j < 3; j++)
837 {
838 temp +=g_contr[i][j]*dtraceK[j];
839 temp2 +=g_contr[j][i]*dTheta[j];
840 }
841 dtGhat[i] = -4./3.*alpha*temp + 2*alpha*temp2 + advec*(beta[0]*dGhat[0][i] + beta[1]*dGhat[1][i] + beta[2]*dGhat[2][i]);
842 for (int l = 0; l < 3; l++)
843 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]);
844 }
845
846 double ov[3];
847 for (int k = 0; k < 3; k++)
848 {
849 double temp=0;
850 for (int i = 0; i < 3; i++)
851 for (int j = 0; j < 3; j++) temp += g_contr[i][j]*dAex[k][i][j];
852 ov[k] = 2*alpha*temp;
853 }
854 for (int i = 0; i < 3; i++)
855 for (int j = 0; j < 3; j++) dtGhat[i] += CCZ4sk*g_contr[i][j]*ov[j];
856
857 double dtbb[3];
858 for (int i = 0; i < 3; i++)
859 {
860 dtbb[i] = CCZ4xi*dtGhat[i] + CCZ4bs*advec*( 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]);
861 dtbb[i] *= CCZ4sk;
862 }
863
864 double dtA[3] = {0,0,0};
865 double dK0[3] = {0,0,0};
866 for (int i = 0; i < 3; i++)
867 {
868 dtA[i] = -alpha*alpha*fa*(dtraceK[i] - dK0[i] - CCZ4c*2*dTheta[i]) + advec_aux*(beta[0]*dAA[0][i] + beta[1]*dAA[1][i] + beta[2]*dAA[2][i]); // Baojiu-01-08-22: alpha*fa -> alpha^2*fa
869 }
870 for (int k = 0; k < 3; k++)
871 {
872 double temp=0;
873 for (int i = 0; i < 3; i++)
874 for (int j = 0; j < 3; j++) temp += g_contr[i][j]*dAex[k][i][j];
875 dtA[k] -= CCZ4sk*alpha*alpha*fa*temp; // Baojiu-13-06-23: here there is only one alpha as in the Dumbser paper
876 }
877
878 double dtB[3][3] = {
879 {CCZ4f*gradQin[20][0],CCZ4f*gradQin[21][0],CCZ4f*gradQin[22][0]},
880 {CCZ4f*gradQin[20][1],CCZ4f*gradQin[21][1],CCZ4f*gradQin[22][1]},
881 {CCZ4f*gradQin[20][2],CCZ4f*gradQin[21][2],CCZ4f*gradQin[22][2]}
882 };
883 for (int i = 0; i < 3; i++)
884 for (int k = 0; k < 3; k++)
885 for (int j = 0; j < 3; j++)
886 {
887// dtB[k][i] += CCZ4mu*alpha2 * g_contr[i][j]*( dPP[k][j] - dPP[j][k]); // Baojiu-13-06-23: Option 1 - similar form as in the Dumbser paper
888 dtB[k][i] += CCZ4mu*alpha2 * g_contr[i][j]*( dPP[k][j] - dPP[j][k])/(np*chi); // Baojiu-13-06-23: Option 2 - same values as in the Dumbser paper // Baojiu-27-06-23
889 for (int n = 0; n < 3; n++)
890 for (int l = 0; l < 3; l++) dtB[k][i] -= CCZ4mu*alpha2 * g_contr[i][j]*g_contr[n][l]*(dDD[k][l][j][n] - dDD[l][k][j][n]);
891 }
892 for (int i = 0; i < 3; i++)
893 for (int j = 0; j < 3; j++) dtB[i][j] += CCZ4bs*advec_aux*(beta[0]*dBB[0][i][j] + beta[1]*dBB[1][i][j] + beta[2]*dBB[2][i][j]);
894 for (int i = 0; i < 3; i++)
895 for (int j = 0; j < 3; j++) dtB[i][j] *= CCZ4sk;
896
897 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};
898 for (int i = 0; i < 3; i++)
899 for (int j = 0; j < 3; j++)
900 for (int k = 0; k < 3; k++) dtD[i][j][k] = -alpha*dAex[i][j][k];
901 for (int i = 0; i < 3; i++)
902 for (int j = 0; j < 3; j++)
903 for (int k = 0; k < 3; k++)
904 for (int m = 0; m < 3; m++)
905 {
906 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]);
907 for (int n = 0; n < 3; n++)
908 {
909 dtD[k][i][j] += 1./3*alpha*g_cov[i][j]*g_contr[n][m]*dAex[k][n][m];
910 }
911 }
912 for (int i = 0; i < 3; i++)
913 for (int j = 0; j < 3; j++)
914 for (int k = 0; k < 3; k++) dtD[i][j][k] += advec_aux*(beta[0]*dDD[0][i][j][k] + beta[1]*dDD[1][i][j][k] + beta[2]*dDD[2][i][j][k]);
915
916 double dtP[3] = {0,0,0};
917 for (int i = 0; i < 3; i++) dtP[i] = advec_aux*(beta[0]*dPP[0][i] + beta[1]*dPP[1][i] + beta[2]*dPP[2][i]);
918 for (int k = 0; k < 3; k++)
919 {
920 double temp=0;
921 for (int m = 0; m < 3; m++)
922 for (int n = 0; n < 3; n++) temp += g_contr[m][n]*dAex[k][m][n];
923 dtP[k] += np/3.*chi*alpha*dtraceK[k]; // Baojiu-13-06-23: changed to chi=phi^n
924 dtP[k] += np/3.*chi*alpha*CCZ4sk*temp; // Baojiu-13-06-23: changed to chi=phi^n
925// dtP[k] += np/3.* alpha*CCZ4sk*temp; // Baojiu-13-06-23: changed to chi=phi^n // Baojiu-27-06-23
926 for (int i = 0; i < 3; i++) dtP[k] -= np/6.*chi*(dBB[k][i][i] + dBB[i][k][i]); // Baojiu-13-06-23: changed to chi=phi^n
927 }
928
929 double dtgamma[3][3] = {0,0,0,0,0,0,0,0,0};
930 double dtalpha = 0;
931 double dtbeta[3] = {0,0,0};
932 double dtchi = 0;
933
934 BgradQ[0] = -dtgamma[0][0];
935 BgradQ[1] = -dtgamma[0][1];
936 BgradQ[2] = -dtgamma[0][2];
937 BgradQ[3] = -dtgamma[1][1];
938 BgradQ[4] = -dtgamma[1][2];
939 BgradQ[5] = -dtgamma[2][2];
940 BgradQ[6] = -dtK[0][0];
941 BgradQ[7] = -dtK[0][1];
942 BgradQ[8] = -dtK[0][2];
943 BgradQ[9] = -dtK[1][1];
944 BgradQ[10] = -dtK[1][2];
945 BgradQ[11] = -dtK[2][2]; // ok
946 //BgradQ[12] = 0.0;
947 BgradQ[12] = -dtTheta; // ok
948 for (int i = 0; i < 3; i++) BgradQ[13+i] = -dtGhat[i]; // buggy
949 BgradQ[16] = -dtalpha;
950 for (int i = 0; i < 3; i++) BgradQ[17+i] = -dtbeta[i];
951 for (int i = 0; i < 3; i++) BgradQ[20+i] = -dtbb[i];
952 for (int i = 0; i < 3; i++) BgradQ[23+i] = -dtA[i];
953 BgradQ[26] = -dtB[0][0]; // note: thes are all 0 for default CCZ4sk=0
954 BgradQ[27] = -dtB[0][1];
955 BgradQ[28] = -dtB[0][2];
956 BgradQ[29] = -dtB[1][0];
957 BgradQ[30] = -dtB[1][1];
958 BgradQ[31] = -dtB[1][2];
959 BgradQ[32] = -dtB[2][0];
960 BgradQ[33] = -dtB[2][1];
961 BgradQ[34] = -dtB[2][2];
962
963 BgradQ[35] = -dtD[0][0][0];
964 BgradQ[36] = -dtD[0][0][1];
965 BgradQ[37] = -dtD[0][0][2];
966 BgradQ[38] = -dtD[0][1][1];
967 BgradQ[39] = -dtD[0][1][2];
968 BgradQ[40] = -dtD[0][2][2];
969
970 BgradQ[41] = -dtD[1][0][0];
971 BgradQ[42] = -dtD[1][0][1];
972 BgradQ[43] = -dtD[1][0][2];
973 BgradQ[44] = -dtD[1][1][1];
974 BgradQ[45] = -dtD[1][1][2];
975 BgradQ[46] = -dtD[1][2][2];
976
977 BgradQ[47] = -dtD[2][0][0];
978 BgradQ[48] = -dtD[2][0][1];
979 BgradQ[49] = -dtD[2][0][2];
980 BgradQ[50] = -dtD[2][1][1];
981 BgradQ[51] = -dtD[2][1][2];
982 BgradQ[52] = -dtD[2][2][2];
983 BgradQ[53] = -dtTraceK;
984 BgradQ[54] = -dtchi;
985 for (int i = 0; i < 3; i++) BgradQ[55+i] = -dtP[i];
986 BgradQ[58] = 0;
987 //if (dtP[1]!=0) printf("dtP[1] = %g, dtP[2] = %g \n\n",dtP[1],dtP[2]);
988 //used only in soccz4
989 if (CCZ4SO==1) {
990 for (int i = 23; i < 53; i++) BgradQ[i] = 0.0;
991 for (int i = 0; i < 3; i++) BgradQ[55+i] = 0.0;
992 }
993}
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ j
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ k
KeywordToAvoidDuplicateSymbolsForInlinedFunctions void source(double *S, const double *const Q, const int CCZ4LapseType, const double CCZ4ds, const double CCZ4c, const double CCZ4e, const double CCZ4f, const double CCZ4bs, const double CCZ4sk, const double CCZ4xi, const double CCZ4itau, const double CCZ4eta, const double CCZ4k1, const double CCZ4k2, const double CCZ4k3, const double CCZ4SO) InlineMethod
The source term is one out of two terms that we use in our CCZ4 formulation.
KeywordToAvoidDuplicateSymbolsForInlinedFunctions double maxEigenvalue(const double *const Q, int normal, const double CCZ4e, const double CCZ4ds, const double CCZ4GLMc, const double CCZ4GLMd) InlineMethod
KeywordToAvoidDuplicateSymbolsForInlinedFunctions void ncp(double *BgradQ, const double *const Q, const double *const gradQSerialised, const int normal, const int CCZ4LapseType, const double CCZ4ds, const double CCZ4c, const double CCZ4e, const double CCZ4f, const double CCZ4bs, const double CCZ4sk, const double CCZ4xi, const double CCZ4mu, const double CCZ4SO) InlineMethod