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}};
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}};
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}}};
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}}};
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}}};
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}};
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}};
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] + (1-CCZ4SO)*2.*beta[k]*DD[k][i][j];
330 //Han-24-03-2025: if SO is enabled this term is advection. thus need to be treated differently
331
332 // MATMUL(Aex,Amix)
333 double Atemp[3][3]= {{0.0}};
334#pragma clang loop unroll(disable)
335 for (int i = 0; i < 3; i++)
336 for (int j = 0; j < 3; j++)
337 for (int u = 0; u < 3; u++) Atemp[i][j] += Aex[i][u] * Amix[u][j];
338
339 const double traceK = Q[53];
340
342 double dtK[3][3];
343 #pragma clang loop unroll(disable)
344 for (int i = 0; i < 3; i++)
345 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;
346
347#pragma clang loop unroll(disable)
348 for (int j = 0; j < 3; j++)
349 for (int i = 0; i < 3; i++)
350 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];
351
352 const double K0 = Q[58];
353
354 // Baojiu-01-08-22: alpha*k1 -> k1 below
355 double dtTraceK = -nablanablaalphaSrc + alpha*(RPlusTwoNablaZSrc + traceK*traceK - 2.0*CCZ4c*Theta*traceK) - 3.0*alpha*CCZ4k1*(1.+CCZ4k2)*Theta/alpha;
356
357 // Baojiu-01-08-22: multiplied the second term on the RHS by phi in the equation below
358 double dtchi = (1-CCZ4SO)*(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
359 //Han-24-03-2025: if SO is enabled this term is advection. thus need to be treated differently
360
361 // Baojiu-01-08-22: multiplied the first term on the RHS by alpha in the equation below
362 double dtalpha = -alpha*alpha*fa*(traceK-K0-2.*CCZ4c*Theta)+ (1-CCZ4SO)*(beta[0]*AA[0]+beta[1]*AA[1]+beta[2]*AA[2]); //Han-24-03-2025: if SO is enabled this term is advection. thus need to be treated differently
363
364 double Aupdown = 0;
365 #pragma clang loop unroll(disable)
366 for (int i = 0; i < 3; i++)
367 for (int j = 0; j < 3; j++) Aupdown += Aex[i][j]*Aup[i][j];
368
369 double sumzupaa = 0.0;
370 for (int i = 0; i < 3; i++) sumzupaa += Zup[i]*AA[i];
371 // Baojiu-01-08-22: divided sumzupaa by alpha, and rescaled alpha*k1 -> k1, in the equation below
372 double dtTheta = 0.5*alpha*CCZ4e*CCZ4e*(RPlusTwoNablaZSrc - Aupdown + 2./3.*traceK*traceK) - alpha*(CCZ4c*Theta*traceK + sumzupaa/alpha + CCZ4k1*(2.+CCZ4k2)*Theta/alpha);
373
374 double dtGhat[3] = {0,0,0};
375 #pragma clang loop unroll(disable)
376 for (int i = 0; i < 3; i++)
377 {
378 double temp1=0, temp2=0, temp3=0, temp4=0, temp5=0, temp6=0;
379 #pragma clang loop unroll(disable)
380 for (int m = 0; m < 3; m++)
381 {
382 temp1 += Aup[i][m]*PP[m]/(np*chi); // Baojiu-13-06-23: changed to chi=phi^n
383 temp3 += Aup[i][m]*AA[m]/alpha; // Baojiu-01-08-22: divided AA by alpha
384 temp2 += g_contr[m][i]*(-Theta*AA[m]/alpha -2./3.*traceK*Z[m]); // Baojiu-01-08-22: divided AA by alpha
385 temp4 += g_contr[i][m]*Z[m];
386 temp5 += Gtilde[m]*BB[m][i];
387 for (int n = 0; n < 3; n++) temp6 += Christoffel_tilde[m][n][i]*Aup[m][n];
388 }
389 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
390
391 for (int l = 0; l < 3; l++)
392 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]);
393 }
394
395 double ov[3];
396 #pragma clang loop unroll(disable)
397 for (int k = 0; k < 3; k++)
398 {
399 double temp=0;
400 for (int i = 0; i < 3; i++)
401 for (int j = 0; j < 3; j++) temp += dgup[k][i][j]*Aex[i][j];
402 ov[k] = 2*alpha*temp;
403 }
404
405 // matrix vector multiplication in a loop and add result to existing vector
406 #pragma clang loop unroll(disable)
407 for (int i = 0; i < 3; i++)
408 for (int j = 0; j < 3; j++) dtGhat[i] += CCZ4helper*CCZ4sk*g_contr[i][j]*ov[j]; //19-02-24: add control parameter CCZ4helper
409
410 const double myb[3] = {Q[20], Q[21], Q[22]};
411
412 double dtbb[3];
413 for (int i = 0; i < 3; i++) dtbb[i] = CCZ4sk*(CCZ4xi*dtGhat[i] - CCZ4eta*myb[i]);
414
415 double dtbeta[3];
416#pragma clang loop unroll(disable)
417 for (int i = 0; i < 3; i++) dtbeta[i] = CCZ4f*myb[i];
418 for (int i = 0; i < 3; i++) dtbeta[i] += CCZ4bs*(1-CCZ4SO)*(beta[0]*BB[0][i] + beta[1]*BB[1][i] + beta[2]*BB[2][i]); //Han-24-03-2025: if SO is enabled this term is advection. thus need to be treated differently
419 for (int i = 0; i < 3; i++) dtbeta[i] *= CCZ4sk;
420
421 // Auxiliary variables
422 double dtA[3] ={0,0,0};
423#pragma clang loop unroll(disable)
424 for (int i = 0; i < 3; i++)
425 {
426 // Baojiu-01-08-22: in the equation below: fa -> 2*fa
427 dtA[i] = -alpha*AA[i]*(2.*fa+alpha*faa)*(traceK - K0 - 2.*CCZ4c*Theta);
428 for (int j = 0; j < 3; j++) dtA[i] += BB[i][j]*AA[j];
429 }
430
431 for (int k = 0; k < 3; k++)
432 {
433 double temp = 0;
434 for (int i = 0; i < 3; i++)
435 for (int j = 0; j < 3; j++) temp+= dgup[k][i][j]*Aex[i][j];
436 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
437 }
438
439 double dtB[3][3] = {{0.0}};
440 #pragma clang loop unroll(disable)
441 for (int i = 0; i < 3; i++)
442 for (int j = 0; j < 3; j++)
443 for (int u = 0; u < 3; u++) dtB[i][j] += CCZ4sk*(BB[i][u] * BB[u][j]);
444
445 double dtD[3][3][3] = {{{0.0}}};
446#pragma clang loop unroll(disable)
447 for (int m = 0; m < 3; m++)
448 for (int j = 0; j < 3; j++)
449 for (int i = 0; i < 3; i++)
450 for (int k = 0; k < 3; k++)
451 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
452
453#pragma clang loop unroll(disable)
454 for (int j = 0; j < 3; j++)
455 for (int i = 0; i < 3; i++)
456 for (int k = 0; k < 3; k++)
457 {
458 dtD[k][i][j] -= AA[k]*Aex[i][j]; // Baojiu-01-08-22: alpha*AA -> alpha
459 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];
460 }
461
462 double dtP[3] = {0,0,0};
463 #pragma clang loop unroll(disable)
464 for (int i = 0; i < 3; i++)
465 for (int j = 0; j < 3; j++) dtP[i] += BB[i][j]*PP[j];
466
467#pragma clang loop unroll(disable)
468 for (int k = 0; k < 3; k++)
469 {
470 double temp=0;
471 for (int i = 0; i < 3; i++)
472 for (int j = 0; j < 3; j++) temp += dgup[k][i][j]*Aex[i][j];
473 dtP[k] += np/3.*AA[k]*traceK*chi; // Baojiu-13-06-23: changed to chi=phi^n
474 dtP[k] += np/3.*(alpha*traceK-traceB)*PP[k]; // Baojiu-13-06-23: changed to chi=phi^n
475 dtP[k] += np/3.*CCZ4helper*alpha*CCZ4sk*temp*chi; // Baojiu-13-06-23: changed to chi=phi^n //30-07-24: add CCZ4helper
476// dtP[k] += np/3.*alpha*CCZ4sk*temp; // Baojiu-13-06-23: changed to chi=phi^n // Baojiu-27-06-23
477 }
478
479 S[0] = dtgamma[0][0];
480 S[1] = dtgamma[0][1];
481 S[2] = dtgamma[0][2];
482 S[3] = dtgamma[1][1];
483 S[4] = dtgamma[1][2];
484 S[5] = dtgamma[2][2];
485 S[6] = dtK[0][0];
486 S[7] = dtK[0][1];
487 S[8] = dtK[0][2];
488 S[9] = dtK[1][1];
489 S[10] = dtK[1][2];
490 S[11] = dtK[2][2];
491 //S[12] = 0.0;
492 S[12] = dtTheta;
493 for (int i = 0; i < 3; i++) S[13+i] = dtGhat[i];
494 S[16] = dtalpha;
495 for (int i = 0; i < 3; i++) S[17+i] = dtbeta[i];
496 for (int i = 0; i < 3; i++) S[20+i] = dtbb[i];
497 for (int i = 0; i < 3; i++) S[23+i] = dtA[i];
498 S[26] = dtB[0][0];
499 S[27] = dtB[0][1];
500 S[28] = dtB[0][2];
501 S[29] = dtB[1][0];
502 S[30] = dtB[1][1];
503 S[31] = dtB[1][2];
504 S[32] = dtB[2][0];
505 S[33] = dtB[2][1];
506 S[34] = dtB[2][2];
507 S[35] = dtD[0][0][0];
508 S[36] = dtD[0][0][1];
509 S[37] = dtD[0][0][2];
510 S[38] = dtD[0][1][1];
511 S[39] = dtD[0][1][2];
512 S[40] = dtD[0][2][2];
513 S[41] = dtD[1][0][0];
514 S[42] = dtD[1][0][1];
515 S[43] = dtD[1][0][2];
516 S[44] = dtD[1][1][1];
517 S[45] = dtD[1][1][2];
518 S[46] = dtD[1][2][2];
519 S[47] = dtD[2][0][0];
520 S[48] = dtD[2][0][1];
521 S[49] = dtD[2][0][2];
522 S[50] = dtD[2][1][1];
523 S[51] = dtD[2][1][2];
524 S[52] = dtD[2][2][2];
525 S[53] = dtTraceK;
526 S[54] = dtchi;
527 for (int i = 0; i < 3; i++) S[55+i] = dtP[i];
528 S[58] = 0;
529 //used only in soccz4
530 if (CCZ4SO==1){
531 for (int i = 23; i < 53; i++) S[i] = 0.0;
532 for (int i = 0; i < 3; i++) S[55+i] = 0.0;
533 }
534}
535
536
537void applications::exahype2::ccz4::ncp(double* BgradQ, const double* const Q, const double* const gradQSerialised, const int normal,
538 const int CCZ4LapseType,
539 const double CCZ4ds,
540 const double CCZ4c,
541 const double CCZ4e,
542 const double CCZ4f,
543 const double CCZ4bs,
544 const double CCZ4sk,
545 const double CCZ4xi,
546 const double CCZ4mu,
547 const double CCZ4SO
548 )
549{
550 const double CCZ4helper = (CCZ4SO==1)? 0.0 : 1.0; //19-02-24: add control parameter CCZ4helperSet
551 // Han 09-2022: add a parameter to switch off or on the advection here
552 const double advec=0.0;
553 const double advec_aux=0.0;
554
555// const double alpha = std::exp(std::fmax(-20., std::fmin(20.,Q[16])));
556 const double alpha = std::fmax(1e-4,Q[16]); // Baojiu-01-08-22: changed to make sure alpha>0.0001
557
558 //printf("alpha %f\n",alpha);
559 const double alpha2 = alpha*alpha;
560 double fa = 1.0;
561 double faa = 0.0;
562 if (CCZ4LapseType==1)
563 {
564 fa = 2./alpha;
565 faa = -fa/alpha;
566 }
567
568 constexpr int nVar(59);
569
570 double gradQin[59][3] = {{0}};
571
572 // De-serialise input data and fill static array
573 // FIXME the use of 2D arrays can be avoided: all terms not in the normal are 0
574 for (int i=0; i<nVar; i++) gradQin[i][normal] = gradQSerialised[i+normal*nVar];
575
576 // Note g_cov is symmetric
577 const double g_cov[3][3] = { {Q[0], Q[1], Q[2]}, {Q[1], Q[3], Q[4]}, {Q[2], Q[4], Q[5]} };
578 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]);
579
580 const double g_contr[3][3] = {
581 { ( 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},
582 {-( 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},
583 {-(-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}
584 };
585
586 // NOTE Aex is symmetric
587 double Aex[3][3] = { {Q[6], Q[7], Q[8]}, {Q[7], Q[9], Q[10]}, {Q[8], Q[10], Q[11]} };
588
589 double traceA = 0;
590 for (int i=0;i<3;i++)
591 for (int j=0;j<3;j++) traceA+=g_contr[i][j]*Aex[i][j];
592 //traceA *= 1./3;
593
594 for (int i=0;i<3;i++)
595 for (int j=0;j<3;j++) Aex[i][j] -= 1./3. * traceA * g_cov[i][j];
596
597 // Matrix multiplications Amix = matmul(g_contr, Aex) Aup = matmul(g_contr, mytranspose(Amix))
598 double Amix[3][3]={{0.0}};
599 for (int i = 0; i < 3; i++)
600 for (int j = 0; j < 3; j++)
601 for (int u = 0; u < 3; u++) Amix[i][j] += g_contr[i][u] * Aex[u][j];
602 double Aup[3][3]={{0.0}};
603 for (int i = 0; i < 3; i++)
604 for (int j = 0; j < 3; j++)
605 for (int u = 0; u < 3; u++) Aup[i][j] += g_contr[i][u] * Amix[j][u]; // Note the transposition is in the indices
606
607 const double DD[3][3][3] = {
608 {{Q[35], Q[36], Q[37]}, {Q[36], Q[38], Q[39]}, {Q[37], Q[39], Q[40]}},
609 {{Q[41], Q[42], Q[43]}, {Q[42], Q[44], Q[45]}, {Q[43], Q[45], Q[46]}},
610 {{Q[47], Q[48], Q[49]}, {Q[48], Q[50], Q[51]}, {Q[49], Q[51], Q[52]}}
611 };
612 double dgup[3][3][3] = {{{0.0}}};
613 for (int k = 0; k < 3; k++)
614 for (int m = 0; m < 3; m++)
615 for (int l = 0; l < 3; l++)
616 for (int n = 0; n < 3; n++)
617 for (int j = 0; j < 3; j++) dgup[k][m][l] -= g_contr[m][n]*g_contr[j][l]*2*DD[k][n][j];
618
619 // Baojiu-01-08-22: changed so that phi^2>0.0001
620// const double phi = std::exp(std::fmax(-20., std::fmin(20.,Q[54])));
621 const double np = 1.0; // Baojiu-13-06-23: changed to chi=phi^n
622 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
623 const double chi = pow(phi,np); // Baojiu-13-06-23: changed to chi=phi^n
624 const double phi2 = phi*phi;
625
626 //const double PP[3] = {Q[55], Q[56], Q[57]};
627 double dPP[3][3] = {
628 {gradQin[55][0],gradQin[56][0],gradQin[57][0]},
629 {gradQin[55][1],gradQin[56][1],gradQin[57][1]},
630 {gradQin[55][2],gradQin[56][2],gradQin[57][2]}
631 };
632
633 //const double AA[3] = {Q[23], Q[24], Q[25]};
634 double dAA[3][3] = {
635 {gradQin[23][0],gradQin[24][0],gradQin[25][0]},
636 {gradQin[23][1],gradQin[24][1],gradQin[25][1]},
637 {gradQin[23][2],gradQin[24][2],gradQin[25][2]}
638 };
639
640// double Christoffel[3][3][3] = {{{0.0}}};
641 double Christoffel_tilde[3][3][3] = {{{0.0}}};
642// double Christoffel_kind1[3][3][3] = {{{0.0}}};
643
644 for (int j = 0; j < 3; j++)
645 for (int i = 0; i < 3; i++)
646 for (int k = 0; k < 3; k++)
647 {
648// Christoffel_kind1[i][j][k] = DD[k][i][j] + DD[j][i][k] - DD[i][j][k];
649 for (int l = 0; l < 3; l++)
650 {
651 Christoffel_tilde[i][j][k] += g_contr[k][l] * ( DD[i][j][l] + DD[j][i][l] - DD[l][i][j] );
652// 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;
653 }
654 }
655
656 double Gtilde[3] = {0,0,0};
657 for (int l = 0; l < 3; l++)
658 for (int j = 0; j < 3; j++)
659 for (int i = 0; i < 3; i++) Gtilde[i] += g_contr[j][l] * Christoffel_tilde[j][l][i];
660
661 const double Ghat[3] = {Q[13], Q[14], Q[15]};
662
663 double Z[3] = {0,0,0};
664 for (int i=0;i<3;i++)
665 for (int j=0;j<3;j++) Z[i] += 0.5*( g_cov[i][j]* (Ghat[j] - Gtilde[j]));
666 double Zup[3] = {0,0,0};
667 for (int i=0;i<3;i++)
668 for (int j=0;j<3;j++) Zup[i] += phi2 * g_contr[i][j] * Z[j];
669
670 const double dDD[3][3][3][3] = {
671 {
672 {
673 {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]},
674 },
675 {
676 {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]},
677 },
678 {
679 {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]}
680 }
681 },
682 {
683 {
684 {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]},
685 },
686 {
687 {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]},
688 },
689 {
690 {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]}
691 }
692 },
693 {
694 {
695 {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]},
696 },
697 {
698 {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]},
699 },
700 {
701 {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]}
702 }
703 }
704 };
705
706 // =====================================================
707 // Baojiu-13-06-23 COMMENTED OUT ALL LINES IN THIS BLOCK
708 // -----------------------------------------------------
709//
710// 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};
711// 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};
712// for (int i = 0; i < 3; i++)
713// for (int j = 0; j < 3; j++)
714// for (int m = 0; m < 3; m++)
715// for (int k = 0; k < 3; k++)
716// {
717// dChristoffelNCP [k][i][j][m] = 0;
718// dChristoffel_tildeNCP[k][i][j][m] = 0;
719// for (int l = 0; l < 3; l++)
720// {
721// // Baojiu-01-08-22: multiplied the original dPP terms by 1/phi //a small issue
722// dChristoffelNCP[k][i][j][m] += 0.5*g_contr[m][l] * (
723// 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]
724// - 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 );
725//
726// 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]);
727//
728// }
729// }
730//
731// double RiemannNCP[3][3][3][3] = {0};
732// for (int i = 0; i < 3; i++)
733// for (int j = 0; j < 3; j++)
734// for (int m = 0; m < 3; m++)
735// for (int k = 0; k < 3; k++) RiemannNCP[i][k][j][m] = dChristoffelNCP[k][i][j][m] - dChristoffelNCP[j][i][k][m];
736//
737// double RicciNCP[3][3] = {0,0,0,0,0,0,0,0,0};
738// for (int m = 0; m < 3; m++)
739// for (int n = 0; n < 3; n++)
740// for (int l = 0; l < 3; l++) RicciNCP[m][n] += RiemannNCP[m][l][n][l];
741//
742// double dGtildeNCP[3][3] = {0,0,0,0,0,0,0,0,0};
743// for (int i = 0; i < 3; i++)
744// for (int k = 0; k < 3; k++)
745// for (int j = 0; j < 3; j++)
746// for (int l = 0; l < 3; l++) dGtildeNCP[k][i] += g_contr[j][l]*dChristoffel_tildeNCP[k][j][l][i];
747//
748// double dZNCP[3][3] = {0,0,0,0,0,0,0,0,0};
749// for (int j = 0; j < 3; j++)
750// for (int i = 0; i < 3; i++)
751// for (int k = 0; k < 3; k++) dZNCP[k][i] += CCZ4ds*0.5*g_cov[i][j]*(dGhat[k][j]-dGtildeNCP[k][j]);
752//
753// double RicciPlusNablaZNCP[3][3];
754// for (int i = 0; i < 3; i++)
755// for (int j = 0; j < 3; j++) RicciPlusNablaZNCP[i][j] = RicciNCP[i][j] + dZNCP[i][j] + dZNCP[j][i];
756//
757 // -----------------------------------------------------
758 // Baojiu-13-06-23 COMMENTED OUT ALL LINES IN THIS BLOCK
759 // =====================================================
760
761 const double dGhat[3][3] = {
762 {gradQin[13][0],gradQin[14][0],gradQin[15][0]},
763 {gradQin[13][1],gradQin[14][1],gradQin[15][1]},
764 {gradQin[13][2],gradQin[14][2],gradQin[15][2]}
765 };
766
767 // ================================
768 // Baojiu-21-09-22 START OF CHANGES
769 // --------------------------------
770 double RicciPlusNablaZNCP[3][3] = {{0.0}};
771 for (int i = 0; i < 3; i++)
772 for (int j = 0; j < 3; j++)
773 {
774 RicciPlusNablaZNCP[i][j] += 0.5*(dPP[i][j] + dPP[j][i])/(np*chi); // Baojiu-13-06-23: changed to chi=phi^n
775 for (int k = 0; k < 3; k++)
776 {
777 RicciPlusNablaZNCP[i][j] += 0.5*( g_cov[k][i]*dGhat[j][k] + g_cov[k][j]*dGhat[i][k] );
778 for (int l = 0; l < 3; l++)
779 {
780 RicciPlusNablaZNCP[i][j] -= 0.5*g_contr[k][l]*( dDD[k][l][i][j] + dDD[l][k][i][j] ); // explicit symmetrisation
781 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
782 }
783 }
784 }
785 // --------------------------------
786 // Baojiu-21-09-22 END OF CHANGES
787 // ================================
788
789 double RPlusTwoNablaZNCP = 0;
790 for (int i = 0; i < 3; i++)
791 for (int j = 0; j < 3; j++) RPlusTwoNablaZNCP += g_contr[i][j]*RicciPlusNablaZNCP[i][j];
792 RPlusTwoNablaZNCP*=phi2;
793
794 double nablaijalphaNCP[3][3];
795 for (int i = 0; i < 3; i++)
796 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
797
798 double nablanablaalphaNCP = 0;
799 for (int i = 0; i < 3; i++)
800 for (int j = 0; j < 3; j++) nablanablaalphaNCP += g_contr[i][j]*nablaijalphaNCP[i][j];
801 nablanablaalphaNCP*=phi2;
802
803 double SecondOrderTermsNCP[3][3];
804 for (int i = 0; i < 3; i++)
805 for (int j = 0; j < 3; j++) SecondOrderTermsNCP[i][j] = -nablaijalphaNCP[i][j] + alpha*RicciPlusNablaZNCP[i][j];
806
807 double traceNCP = 0;
808 for (int i = 0; i < 3; i++)
809 for (int j = 0; j < 3; j++) traceNCP += g_contr[i][j]*SecondOrderTermsNCP[i][j];
810
811 for (int i = 0; i < 3; i++)
812 for (int j = 0; j < 3; j++) SecondOrderTermsNCP[i][j] -= 1./3 * traceNCP * g_cov[i][j];
813
814 const double beta[3] = {Q[17], Q[18], Q[19]};
815
816 const double dAex[3][3][3] = {
817 {{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]}},
818 {{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]}},
819 {{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]}}
820 };
821
822 double dtK[3][3];
823 for (int i = 0; i < 3; i++)
824 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
825
826 const double dtraceK[3] = {gradQin[53][0], gradQin[53][1], gradQin[53][2]};
827
828 double dtTraceK = -nablanablaalphaNCP + alpha*RPlusTwoNablaZNCP + advec*(beta[0]*dtraceK[0] + beta[1]*dtraceK[1] + beta[2]*dtraceK[2]);
829
830 const double BB[3][3] = {
831 {Q[26], Q[27], Q[28]}, {Q[29], Q[30], Q[31]}, {Q[32], Q[33], Q[34]}
832 };
833
834 double traceB = BB[0][0] + BB[1][1] + BB[2][2]; // TODO direct from Q!
835
836 double Aupdown = 0;
837 for (int i = 0; i < 3; i++)
838 for (int j = 0; j < 3; j++) Aupdown += Aex[i][j]*Aup[i][j];
839
840 const double dTheta[3] = {gradQin[12][0],gradQin[12][1],gradQin[12][2]};
841 double dtTheta = 0.5*alpha*CCZ4e*CCZ4e*RPlusTwoNablaZNCP + advec*(beta[0]*dTheta[0] + beta[1]*dTheta[1] + beta[2]*dTheta[2]);
842
843// double divAupNCP[3] = {0,0,0};
844// for (int i = 0; i < 3; i++)
845// for (int j = 0; j < 3; j++)
846// for (int l = 0; l < 3; l++)
847// for (int k = 0; k < 3; k++) divAupNCP[i] += g_contr[i][l]*g_contr[j][k]*dAex[j][l][k];
848
849 const double dBB[3][3][3] = {
850 {
851 {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]}
852 },
853 {
854 {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]}
855 },
856 {
857 {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]}
858 }
859 };
860
861 double dtGhat[3] = {0,0,0};
862 for (int i = 0; i < 3; i++)
863 {
864 double temp=0, temp2=0;
865 for (int j = 0; j < 3; j++)
866 {
867 temp +=g_contr[i][j]*dtraceK[j];
868 temp2 +=g_contr[j][i]*dTheta[j];
869 }
870 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]);
871 for (int l = 0; l < 3; l++)
872 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]);
873 }
874
875 double ov[3];
876 for (int k = 0; k < 3; k++)
877 {
878 double temp=0;
879 for (int i = 0; i < 3; i++)
880 for (int j = 0; j < 3; j++) temp += g_contr[i][j]*dAex[k][i][j];
881 ov[k] = 2*alpha*temp;
882 }
883 for (int i = 0; i < 3; i++)
884 for (int j = 0; j < 3; j++) dtGhat[i] += CCZ4helper*CCZ4sk*g_contr[i][j]*ov[j]; //19-02-24: add control parameter CCZ4helper
885
886 double dtbb[3];
887 for (int i = 0; i < 3; i++)
888 {
889 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]);
890 dtbb[i] *= CCZ4sk;
891 }
892
893 double dtA[3] = {0,0,0};
894 double dK0[3] = {0,0,0};
895 for (int i = 0; i < 3; i++)
896 {
897 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
898 }
899 for (int k = 0; k < 3; k++)
900 {
901 double temp=0;
902 for (int i = 0; i < 3; i++)
903 for (int j = 0; j < 3; j++) temp += g_contr[i][j]*dAex[k][i][j];
904 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
905 }
906
907 double dtB[3][3] = {
908 {CCZ4f*gradQin[20][0],CCZ4f*gradQin[21][0],CCZ4f*gradQin[22][0]},
909 {CCZ4f*gradQin[20][1],CCZ4f*gradQin[21][1],CCZ4f*gradQin[22][1]},
910 {CCZ4f*gradQin[20][2],CCZ4f*gradQin[21][2],CCZ4f*gradQin[22][2]}
911 };
912 for (int i = 0; i < 3; i++)
913 for (int k = 0; k < 3; k++)
914 for (int j = 0; j < 3; j++)
915 {
916// 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
917 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
918 for (int n = 0; n < 3; n++)
919 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
920 }
921 for (int i = 0; i < 3; i++)
922 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]);
923 for (int i = 0; i < 3; i++)
924 for (int j = 0; j < 3; j++) dtB[i][j] *= CCZ4sk;
925
926 double dtD[3][3][3] = {{{0.0}}};
927 for (int i = 0; i < 3; i++)
928 for (int j = 0; j < 3; j++)
929 for (int k = 0; k < 3; k++) dtD[i][j][k] = -alpha*dAex[i][j][k];
930 for (int i = 0; i < 3; i++)
931 for (int j = 0; j < 3; j++)
932 for (int k = 0; k < 3; k++)
933 for (int m = 0; m < 3; m++)
934 {
935 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]);
936 for (int n = 0; n < 3; n++)
937 {
938 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
939 }
940 }
941 for (int i = 0; i < 3; i++)
942 for (int j = 0; j < 3; j++)
943 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]);
944
945 double dtP[3] = {0,0,0};
946 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]);
947 for (int k = 0; k < 3; k++)
948 {
949 double temp=0;
950 for (int m = 0; m < 3; m++)
951 for (int n = 0; n < 3; n++) temp += g_contr[m][n]*dAex[k][m][n];
952 dtP[k] += np/3.*chi*alpha*dtraceK[k]; // Baojiu-13-06-23: changed to chi=phi^n
953 dtP[k] += np/3.*CCZ4helper*chi*alpha*CCZ4sk*temp; // Baojiu-13-06-23: changed to chi=phi^n //30-07-24: add CCZ4helper
954// dtP[k] += np/3.* alpha*CCZ4sk*temp; // Baojiu-13-06-23: changed to chi=phi^n // Baojiu-27-06-23
955 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
956 }
957
958 double dtgamma[3][3] = {{0.0}};
959 double dtalpha = 0;
960 double dtbeta[3] = {0.0};
961 double dtchi = 0;
962
963 BgradQ[0] = -dtgamma[0][0];
964 BgradQ[1] = -dtgamma[0][1];
965 BgradQ[2] = -dtgamma[0][2];
966 BgradQ[3] = -dtgamma[1][1];
967 BgradQ[4] = -dtgamma[1][2];
968 BgradQ[5] = -dtgamma[2][2];
969 BgradQ[6] = -dtK[0][0];
970 BgradQ[7] = -dtK[0][1];
971 BgradQ[8] = -dtK[0][2];
972 BgradQ[9] = -dtK[1][1];
973 BgradQ[10] = -dtK[1][2];
974 BgradQ[11] = -dtK[2][2]; // ok
975 //BgradQ[12] = 0.0;
976 BgradQ[12] = -dtTheta; // ok
977 for (int i = 0; i < 3; i++) BgradQ[13+i] = -dtGhat[i]; // buggy
978 BgradQ[16] = -dtalpha;
979 for (int i = 0; i < 3; i++) BgradQ[17+i] = -dtbeta[i];
980 for (int i = 0; i < 3; i++) BgradQ[20+i] = -dtbb[i];
981 for (int i = 0; i < 3; i++) BgradQ[23+i] = -dtA[i];
982 BgradQ[26] = -dtB[0][0]; // note: thes are all 0 for default CCZ4sk=0
983 BgradQ[27] = -dtB[0][1];
984 BgradQ[28] = -dtB[0][2];
985 BgradQ[29] = -dtB[1][0];
986 BgradQ[30] = -dtB[1][1];
987 BgradQ[31] = -dtB[1][2];
988 BgradQ[32] = -dtB[2][0];
989 BgradQ[33] = -dtB[2][1];
990 BgradQ[34] = -dtB[2][2];
991
992 BgradQ[35] = -dtD[0][0][0];
993 BgradQ[36] = -dtD[0][0][1];
994 BgradQ[37] = -dtD[0][0][2];
995 BgradQ[38] = -dtD[0][1][1];
996 BgradQ[39] = -dtD[0][1][2];
997 BgradQ[40] = -dtD[0][2][2];
998
999 BgradQ[41] = -dtD[1][0][0];
1000 BgradQ[42] = -dtD[1][0][1];
1001 BgradQ[43] = -dtD[1][0][2];
1002 BgradQ[44] = -dtD[1][1][1];
1003 BgradQ[45] = -dtD[1][1][2];
1004 BgradQ[46] = -dtD[1][2][2];
1005
1006 BgradQ[47] = -dtD[2][0][0];
1007 BgradQ[48] = -dtD[2][0][1];
1008 BgradQ[49] = -dtD[2][0][2];
1009 BgradQ[50] = -dtD[2][1][1];
1010 BgradQ[51] = -dtD[2][1][2];
1011 BgradQ[52] = -dtD[2][2][2];
1012 BgradQ[53] = -dtTraceK;
1013 BgradQ[54] = -dtchi;
1014 for (int i = 0; i < 3; i++) BgradQ[55+i] = -dtP[i];
1015 BgradQ[58] = 0;
1016 //if (dtP[1]!=0) printf("dtP[1] = %g, dtP[2] = %g \n\n",dtP[1],dtP[2]);
1017 //used only in soccz4
1018 if (CCZ4SO==1) {
1019 for (int i = 23; i < 53; i++) BgradQ[i] = 0.0;
1020 for (int i = 0; i < 3; i++) BgradQ[55+i] = 0.0;
1021 }
1022}
applications::exahype2::ccz4::source
static 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)
The source term is one out of two terms that we use in our CCZ4 formulation.
applications::exahype2::ccz4::ncp
static 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)
applications::exahype2::ccz4::maxEigenvalue
static double maxEigenvalue(const double *const Q, int normal, const double CCZ4e, const double CCZ4ds, const double CCZ4GLMc, const double CCZ4GLMd)
euler.j
j
Definition euler.py:95