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