Peano
CCZ4Kernels.cpph
Go to the documentation of this file.
1 #pragma once
2 
3 #include <algorithm>
4 
5 
6 double 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
25 void 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 
537 void 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 }
TP::min
static double min(double const x, double const y)
Definition: FuncAndJacobian.cpp:21
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.u
u
Definition: euler.py:113
euler.j
j
Definition: euler.py:95