Peano 4
Loading...
Searching...
No Matches
TP_bindding.h
Go to the documentation of this file.
1#pragma once
2
4//#include "AbstractFiniteVolumeCCZ4.h"
5//The functions in this namespace works as binddings of two puncture library
6
7namespace TP_bindding {
8
9 //calculate the soccz4 quantites use the real physics quantities given in TP lib
10 inline void SOCCZ4Cal(double* __restrict__ Q){
11 //take care: the metric below are all non-conformal, i.e. without tilde
12 double g_cov[3][3] = { {Q[0], Q[1], Q[2]}, {Q[1], Q[3], Q[4]}, {Q[2], Q[4], Q[5]} };
13 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];
14 double invdet = 1./det;
15
16 double g_contr[3][3] = {
17 { ( 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},
18 {-( 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},
19 {-(-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}
20 };
21
22 double phi=std::pow(det,-1./6.);
23 double phisq=phi*phi;
24
25 double traceK=0;
26 double Kex[3][3] = { {Q[6], Q[7], Q[8]}, {Q[7], Q[9], Q[10]}, {Q[8], Q[10], Q[11]} };
27 for (int i=0;i<3;i++)
28 for (int j=0;j<3;j++) traceK += g_contr[i][j]*Kex[i][j];
29
30 //now adjust the solution
31 for (int i=0;i<6;i++) {
32 Q[i]=phisq*Q[i]; //\tilde{\gamma}_ij
33 Q[i+6]= phisq*Q[i+6]-1./3. * traceK * Q[i]; //\tilde{A}_ij
34 }
35
36 //Q[53]=traceK; Q[54]=std::log(phi); Q[16]=std::log(Q[16]);
37 const double np = 1.0;
38 Q[53]=traceK; Q[54]=std::pow(phi,np); //Q[16]=phi;
39 }
40
41 //use this to calculate the gradient
42 inline void GradientCal(const double* X, double* __restrict__ Q, double* LgradQ, int nVars, TP::TwoPunctures* tp, bool low_res=false)
43 {
44 constexpr double epsilon = 1e-4;
45 double Qp1[nVars],Qm1[nVars],Qp2[nVars],Qm2[nVars];
46
47 for (int d=0;d<3;d++)
48 {
49 double xp1[3]={X[0],X[1],X[2]}; double xp2[3]={X[0],X[1],X[2]};
50 double xm1[3]={X[0],X[1],X[2]}; double xm2[3]={X[0],X[1],X[2]};
51 xp1[d]+=epsilon; xp2[d]+=2*epsilon; xm1[d]-=epsilon; xm2[d]-=2*epsilon;
52
53 tp->Interpolate(xp1,Qp1, low_res); SOCCZ4Cal(Qp1);
54 tp->Interpolate(xp2,Qp2, low_res); SOCCZ4Cal(Qp2);
55 tp->Interpolate(xm1,Qm1, low_res); SOCCZ4Cal(Qm1);
56 tp->Interpolate(xm2,Qm2, low_res); SOCCZ4Cal(Qm2);
57
58 for (int i=0; i<nVars; i++) {
59 LgradQ[d*nVars+i] = ( 8.0*Qp1[i] - 8.0*Qm1[i] + Qm2[i] - Qp2[i] )/(12.0*epsilon);
60 }
61 }
62 }
63
64 //calculate the auxiliary variables using Q and its gradient
65 inline void AuxiliaryCal(double* __restrict__ Q, double* LgradQ, int nVars){
66 double gradQ[3][59]={0};
67 for (int d=0;d<3;d++)
68 for (int i=0;i<nVars;i++) {gradQ[d][i]=LgradQ[d*nVars+i];}
69
70 //here all quantites are coformal, i.e. with the tilde
71 double g_cov[3][3] = { {Q[0], Q[1], Q[2]}, {Q[1], Q[3], Q[4]}, {Q[2], Q[4], Q[5]} };
72 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];
73 double invdet = 1./det;
74
75 double g_contr[3][3] = {
76 { ( 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},
77 {-( 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},
78 {-(-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}
79 };
80
81 //A[i]=A_i=\partial_i ln\alpha
82 for (int i=0;i<3;i++) Q[23+i]=gradQ[i][16];
83
84 //B[i][j]=B_i^j=\partial_i \beta^j //the component order of this variable is different between fortran and C++!
85 for (int i=0;i<3;i++)
86 for (int j=0;j<3;j++) Q[26+i*3+j]=gradQ[i][17+j];
87
88 //D[k][i][j]=D_ijk= 0.5 * \partial_k \tilde{\gamma}_ij
89 for (int i=0;i<3;i++)
90 for (int j=0;j<6;j++) Q[35+i*6+j]=0.5*gradQ[i][0+j];
91
92 //P[i]=P_i=\partial_i ln\phi
93 for (int i=0;i<3;i++) Q[55+i]=gradQ[i][54];
94
95 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};
96 double DD[3][3][3] = {
97 {{Q[35], Q[36], Q[37]}, {Q[36], Q[38], Q[39]}, {Q[37], Q[39], Q[40]}},
98 {{Q[41], Q[42], Q[43]}, {Q[42], Q[44], Q[45]}, {Q[43], Q[45], Q[46]}},
99 {{Q[47], Q[48], Q[49]}, {Q[48], Q[50], Q[51]}, {Q[49], Q[51], Q[52]}}};
100 for (int j = 0; j < 3; j++)
101 for (int i = 0; i < 3; i++)
102 for (int k = 0; k < 3; k++)
103 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] );}
104
105 for (int i = 0; i < 3; i++)
106 for (int l = 0; l < 3; l++)
107 for (int j = 0; j < 3; j++) {Q[13+i] += g_contr[j][l] * Christoffel_tilde[j][l][i];}
108 }
109}
110
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ j
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ k
void Interpolate(const double *const pos, double *Q, bool low_res=false)
Interpolation function for an external caller.
void AuxiliaryCal(double *__restrict__ Q, double *LgradQ, int nVars)
Definition TP_bindding.h:65
void SOCCZ4Cal(double *__restrict__ Q)
Definition TP_bindding.h:10
void GradientCal(const double *X, double *__restrict__ Q, double *LgradQ, int nVars, TP::TwoPunctures *tp, bool low_res=false)
Definition TP_bindding.h:42