Peano 4
Loading...
Searching...
No Matches
InitialValues.cpp
Go to the documentation of this file.
1#include "InitialValues.h"
2#include "Properties.h"
7#include "Constants.h"
8
9#include "tarch/la/Vector.h"
10
12
13#include <stdio.h>
14#include <string.h>
15
16#ifdef IncludeTwoPunctures
18#endif
19
21 double * __restrict__ Q, // Q[64+0],
23 double t
24) {
25 constexpr int nVars = 59;
26 constexpr double pi = M_PI;
27 constexpr double peak_number = 2.0;
28 constexpr double ICA = 0.1;
29 double HH = 1.0 - ICA*sin( peak_number*pi*( x[0] - t));
30 double dxH = -peak_number*pi*ICA*cos( peak_number * pi*(x[0] - t));
31 double dxphi = - pow(HH,(-7.0/6.0))*dxH/6.0;
32 double phi = pow(( 1.0 / HH),(1.0/6.0));
33 double Kxx = - 0.5*peak_number*pi*ICA*cos( peak_number * pi*(x[0] - t))/sqrt( 1.0 - ICA*sin( peak_number*pi*( x[0] - t)) );
34 double traceK = Kxx/HH;
35 tarch::memset(Q, .0, nVars*sizeof(double));
36 Q[0] = phi*phi*HH ; //\tilde(\gamma)_xx
37 Q[3] = phi*phi ; //\tilde(\gamma)_yy
38 Q[5] = phi*phi ; //\tilde(\gamma)_zz
39 Q[6] = phi*phi*(Kxx - 1.0/3.0*traceK*HH ) ; //\tilde(A)_xx
40 Q[9] = phi*phi*(0.0 - 1.0/3.0*traceK*1.0) ; //\tilde(A)_yy
41 Q[11] = phi*phi*(0.0 - 1.0/3.0*traceK*1.0) ; //\tilde(A)_zz
42 Q[16] = sqrt(HH); //\alpha
43 Q[13] = 2.0/(3.0*pow(HH,(5.0/3.0)))*dxH ; //\hat(\Gamma)^x
44 Q[23] = 1.0/(2.0)*dxH*pow(HH,(-1.0/2.0)) ; //A_x
45 Q[35] = pow(HH,(-1.0/3.0))*dxH/3.0 ; //D_xxx
46 Q[38] = phi*dxphi ; //D_xyy
47 Q[40] = phi*dxphi ; //D_xzz
48 Q[53] = traceK; //K
49 Q[54] = phi; //\phi
50 Q[55] = dxphi; //P_x
51}
52
54 double * __restrict__ Q, // Q[64+0],
56 double t
57) {
58 constexpr int nVars = 59;
59 constexpr double pi = M_PI;
60 constexpr double peak_number = 2.0;
61 constexpr double ICA = 0.1;
62 double phase = peak_number*pi*( x[0] - x[1] - t);
63 double H=1+2.0*ICA*sin(phase);
64 tarch::memset(Q, .0, nVars*sizeof(double));
65 Q[0] = pow(H,(-1.0/3.0))*(1+ICA*sin(phase)); //\tilde(\gamma)_xx
66 Q[1] = -pow(H,(-1.0/3.0))*(ICA*sin(phase)); //\tilde(\gamma)_xy
67 Q[3] = pow(H,(-1.0/3.0))*(1+ICA*sin(phase)); //\tilde(\gamma)_yy
68 Q[5] = pow(H,(-1.0/3.0)); //\tilde(\gamma)_zz
69 Q[6] = pow(H,(-5.0/6.0))*(peak_number*pi*ICA/2.0)*cos(phase)*(1.0 - 2.0 * (1+ICA*sin(phase)) / (3.0*H) ); //\tilde(A)_xx
70 Q[7] = pow(H,(-5.0/6.0))*(peak_number*pi*ICA/2.0)*cos(phase)*(2.0 * ICA*sin(phase) / (3.0*H) - 1.0 ); //\tilde(A)_xy
71 Q[9] = pow(H,(-5.0/6.0))*(peak_number*pi*ICA/2.0)*cos(phase)*(1.0 - 2.0 * (1+ICA*sin(phase)) / (3.0*H) ); //\tilde(A)_yy
72 Q[11] = pow(H,(-11.0/6.0))*(peak_number*pi*ICA/3.0)*cos(phase); //\tilde(A)_zz
73 Q[13] = pow(H,(-5.0/3.0))*(4.0*peak_number*pi*ICA/3.0)*cos(phase); //\hat(\Gamma)^x
74 Q[14] = -pow(H,(-5.0/3.0))*(4.0*peak_number*pi*ICA/3.0)*cos(phase); //\hat(\Gamma)^y
75 Q[16] = 0.5*log(H); //ln(\alpha)
76 Q[23] = pow(H,-1.0)*(peak_number*pi*ICA)*cos(phase); //A_x
77 Q[24] = -pow(H,-1.0)*(peak_number*pi*ICA)*cos(phase); //A_y
78 Q[35] = pow(H,(-4.0/3.0))*(peak_number*pi*ICA/6.0)*cos(phase)*(1.0 + 4.0 * ICA*sin(phase) ); //D_xxx
79 Q[36] = -pow(H,(-4.0/3.0))*(peak_number*pi*ICA/6.0)*cos(phase)*(3.0 + 4.0 * ICA*sin(phase) ); //D_xxy
80 Q[38] = pow(H,(-4.0/3.0))*(peak_number*pi*ICA/6.0)*cos(phase)*(1.0 + 4.0 * ICA*sin(phase) ); //D_xyy
81 Q[40] = -pow(H,(-4.0/3.0))*(peak_number*pi*ICA/3.0)*cos(phase); //D_xzz
82 Q[41] = -pow(H,(-4.0/3.0))*(peak_number*pi*ICA/6.0)*cos(phase)*(1.0 + 4.0 * ICA*sin(phase) ); //D_yxx
83 Q[42] = pow(H,(-4.0/3.0))*(peak_number*pi*ICA/6.0)*cos(phase)*(3.0 + 4.0 * ICA*sin(phase) ); //D_yxy
84 Q[41] = -pow(H,(-4.0/3.0))*(peak_number*pi*ICA/6.0)*cos(phase)*(1.0 + 4.0 * ICA*sin(phase) ); //D_yyy
85 Q[46] = pow(H,(-4.0/3.0))*(peak_number*pi*ICA/3.0)*cos(phase); //D_yzz
86 Q[53] = pow(H,(-3.0/2.0))*(peak_number*pi*ICA)*cos(phase); //K
87 Q[54] = -(1.0/6.0)*log(H); //ln(\phi)
88 Q[55] = -pow(H,-1.0)*(peak_number*pi*ICA/3.0)*cos(phase); //P_x
89 Q[56] = pow(H,-1.0)*(peak_number*pi*ICA/3.0)*cos(phase); //P_x
90}
91
93 double * __restrict__ Q, // Q[64+0],
95 double t
96) {
97 constexpr int nVars = 59;
98 constexpr double pi = M_PI;
99 constexpr double peak_number = 2.0;
100 constexpr double ICA = 1e-4;
101 double HH = ICA*sin( peak_number*pi*( X[0] - t));
102 double dxHH = peak_number*pi*ICA*cos( peak_number * pi*(X[0] - t));
103 double dtHH = -peak_number*pi*ICA*cos( peak_number * pi*(X[0] - t));
104 tarch::memset(Q, .0, nVars*sizeof(double));
105 Q[0] = 1.0 ; //\tilde(\gamma)_xx
106 Q[3] = 1+HH ; //\tilde(\gamma)_yy
107 Q[5] = 1-HH ; //\tilde(\gamma)_zz
108 Q[6] = 0.0 ; //\tilde(A)_xx
109 Q[9] = -0.5*dtHH ; //\tilde(A)_yy
110 Q[11] = 0.5*dtHH ; //\tilde(A)_zz
111 Q[16] = log(1.0) ; //ln(\alpha)
112 Q[35] = 0.0 ; //D_xxx
113 Q[38] = 0.5*dxHH ; //D_xyy
114 Q[40] = -0.5*dxHH ; //D_xzz
115 Q[54] = log(1.0) ; //ln(\phi)
116}
117
119 double * __restrict__ Q, // Q[64+0],
121 double t
122) {
123 tarch::memset(Q, .0, 59*sizeof(double));
124 Q[0] = 1.0 ; //\tilde(\gamma)_xx
125 Q[3] = 1.0 ; //\tilde(\gamma)_yy
126 Q[5] = 1.0 ; //\tilde(\gamma)_zz
127 Q[16] = 1.0;
128 Q[54] = 1.0; //\phi
129}
130
131#ifdef IncludeTwoPunctures
132void applications::exahype2::ccz4::ApplyTwoPunctures(
133 double * __restrict__ Q, // Q[64+0],
135 double t,
137 bool low_res
138) {
139 constexpr int nVars = 59;
140 const double coor[3]={X[0],X[1],X[2]};
141 double LgradQ[3*nVars];
142 tarch::memset(Q, .0, nVars*sizeof(double));
143 tp->Interpolate(coor,Q,low_res); //do the interpolate
144 //std::cout << coor[0] <<coor[1] << coor[2] << "\n";
145 //std::cout << "real quantites without tilde" <<"\n";
146 //for (int i=0;i<nVars;i++){std::cout << i <<"\t"<< Q[i] << "\n";}
147 TP_bindding::SOCCZ4Cal(Q); //calculate corresponding soccz4 quantities
148 //std::cout << "after treatment" <<"\n";
149 //for (int i=0;i<nVars;i++){std::cout << i <<"\t" << Q[i] << "\n";}
150 TP_bindding::GradientCal(coor, Q, LgradQ, nVars, tp, low_res); //calculate gradient for auxiliary variables
151 //for (int d=0;d<3;d++)
152 //for (int i=0;i<nVars;i++) {std::cout << d <<"\t" << i <<"\t" << LgradQ[d*nVars+i] << "\n";}
153 TP_bindding::AuxiliaryCal(Q, LgradQ, nVars); //calculate the auxiliary variables
154 //std::cout << "after treatment" <<"\n";
155 //for (int i=0;i<nVars;i++){std::cout << i <<"\t" << Q[i] << "\n";}
156
157 //exit(0);
158}
159#endif
160
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
void linearWave(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &X, double t)
void diagonal_gaugeWave(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t)
void flat(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &X, double t)
void gaugeWave(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &x, double t)
double * memset(double *dest, double ch, size_t byteCount)
Alternative GPU-ready version of memset.
Simple vector class.
Definition Vector.h:134