Peano
Equations.cpp
Go to the documentation of this file.
1 /* TwoPunctures: File "Equations.c"*/
2 
3 #include <stdio.h>
4 #include <stdlib.h>
5 #include <string.h>
6 #include <math.h>
7 #include <ctype.h>
8 #include "libtwopunctures/TwoPunctures.h"
9 
10 namespace TP {
11 using namespace Utilities;
12 
13 /* U.d0[ivar] = U[ivar]; (ivar = 0..nvar-1) */
14 /* U.d1[ivar] = U[ivar]_x; */
15 /* U.d2[ivar] = U[ivar]_y; */
16 /* U.d3[ivar] = U[ivar]_z; */
17 /* U.d11[ivar] = U[ivar]_xx; */
18 /* U.d12[ivar] = U[ivar]_xy; */
19 /* U.d13[ivar] = U[ivar]_xz;*/
20 /* U.d22[ivar] = U[ivar]_yy;*/
21 /* U.d23[ivar] = U[ivar]_yz;*/
22 /* U.d33[ivar] = U[ivar]_zz;*/
23 
24 double
25 TwoPunctures::BY_KKofxyz (double x, double y, double z)
26 {
27 
28  int i, j;
29  double r_plus, r2_plus, r3_plus, r_minus, r2_minus, r3_minus, np_Pp, nm_Pm,
30  Aij, AijAij, n_plus[3], n_minus[3], np_Sp[3], nm_Sm[3];
31 
32  r2_plus = (x - par_b) * (x - par_b) + y * y + z * z;
33  r2_minus = (x + par_b) * (x + par_b) + y * y + z * z;
34  r_plus = sqrt (r2_plus);
35  r_minus = sqrt (r2_minus);
36  r3_plus = r_plus * r2_plus;
37  r3_minus = r_minus * r2_minus;
38 
39  n_plus[0] = (x - par_b) / r_plus;
40  n_minus[0] = (x + par_b) / r_minus;
41  n_plus[1] = y / r_plus;
42  n_minus[1] = y / r_minus;
43  n_plus[2] = z / r_plus;
44  n_minus[2] = z / r_minus;
45 
46  /* dot product: np_Pp = (n_+).(P_+); nm_Pm = (n_-).(P_-) */
47  np_Pp = 0;
48  nm_Pm = 0;
49  for (i = 0; i < 3; i++)
50  {
51  np_Pp += n_plus[i] * par_P_plus[i];
52  nm_Pm += n_minus[i] * par_P_minus[i];
53  }
54  /* cross product: np_Sp[i] = [(n_+) x (S_+)]_i; nm_Sm[i] = [(n_-) x (S_-)]_i*/
55  np_Sp[0] = n_plus[1] * par_S_plus[2] - n_plus[2] * par_S_plus[1];
56  np_Sp[1] = n_plus[2] * par_S_plus[0] - n_plus[0] * par_S_plus[2];
57  np_Sp[2] = n_plus[0] * par_S_plus[1] - n_plus[1] * par_S_plus[0];
58  nm_Sm[0] = n_minus[1] * par_S_minus[2] - n_minus[2] * par_S_minus[1];
59  nm_Sm[1] = n_minus[2] * par_S_minus[0] - n_minus[0] * par_S_minus[2];
60  nm_Sm[2] = n_minus[0] * par_S_minus[1] - n_minus[1] * par_S_minus[0];
61  AijAij = 0;
62  for (i = 0; i < 3; i++)
63  {
64  for (j = 0; j < 3; j++)
65  { /* Bowen-York-Curvature :*/
66  Aij =
67  + 1.5 * (par_P_plus[i] * n_plus[j] + par_P_plus[j] * n_plus[i]
68  + np_Pp * n_plus[i] * n_plus[j]) / r2_plus
69  + 1.5 * (par_P_minus[i] * n_minus[j] + par_P_minus[j] * n_minus[i]
70  + nm_Pm * n_minus[i] * n_minus[j]) / r2_minus
71  - 3.0 * (np_Sp[i] * n_plus[j] + np_Sp[j] * n_plus[i]) / r3_plus
72  - 3.0 * (nm_Sm[i] * n_minus[j] + nm_Sm[j] * n_minus[i]) / r3_minus;
73  if (i == j)
74  Aij -= +1.5 * (np_Pp / r2_plus + nm_Pm / r2_minus);
75  AijAij += Aij * Aij;
76  }
77  }
78 
79  return AijAij;
80 }
81 
82 void
83 TwoPunctures::BY_Aijofxyz (double x, double y, double z, double Aij[3][3])
84 {
85 
86  int i, j;
87  double r_plus, r2_plus, r3_plus, r_minus, r2_minus, r3_minus, np_Pp, nm_Pm,
88  n_plus[3], n_minus[3], np_Sp[3], nm_Sm[3];
89 
90  const double TP4 = TP_epsilon*TP_epsilon*TP_epsilon*TP_epsilon;
91  r2_plus = (x - par_b) * (x - par_b) + y * y + z * z;
92  r2_minus = (x + par_b) * (x + par_b) + y * y + z * z;
93  r2_plus = sqrt (r2_plus*r2_plus + TP4);
94  r2_minus = sqrt (r2_minus*r2_minus + TP4);
95 
96  const double TPT2 = TP_Tiny*TP_Tiny;
97  if (r2_plus < TPT2)
98  r2_plus = TPT2;
99  if (r2_minus < TPT2)
100  r2_minus = TPT2;
101  r_plus = sqrt (r2_plus);
102  r_minus = sqrt (r2_minus);
103  r3_plus = r_plus * r2_plus;
104  r3_minus = r_minus * r2_minus;
105 
106  n_plus[0] = (x - par_b) / r_plus;
107  n_minus[0] = (x + par_b) / r_minus;
108  n_plus[1] = y / r_plus;
109  n_minus[1] = y / r_minus;
110  n_plus[2] = z / r_plus;
111  n_minus[2] = z / r_minus;
112 
113  /* dot product: np_Pp = (n_+).(P_+); nm_Pm = (n_-).(P_-) */
114  np_Pp = 0;
115  nm_Pm = 0;
116  for (i = 0; i < 3; i++)
117  {
118  np_Pp += n_plus[i] * par_P_plus[i];
119  nm_Pm += n_minus[i] * par_P_minus[i];
120  }
121  /* cross product: np_Sp[i] = [(n_+) x (S_+)]_i; nm_Sm[i] = [(n_-) x (S_-)]_i*/
122  np_Sp[0] = n_plus[1] * par_S_plus[2] - n_plus[2] * par_S_plus[1];
123  np_Sp[1] = n_plus[2] * par_S_plus[0] - n_plus[0] * par_S_plus[2];
124  np_Sp[2] = n_plus[0] * par_S_plus[1] - n_plus[1] * par_S_plus[0];
125  nm_Sm[0] = n_minus[1] * par_S_minus[2] - n_minus[2] * par_S_minus[1];
126  nm_Sm[1] = n_minus[2] * par_S_minus[0] - n_minus[0] * par_S_minus[2];
127  nm_Sm[2] = n_minus[0] * par_S_minus[1] - n_minus[1] * par_S_minus[0];
128  for (i = 0; i < 3; i++)
129  {
130  for (j = 0; j < 3; j++)
131  { /* Bowen-York-Curvature :*/
132  Aij[i][j] =
133  + 1.5 * (par_P_plus[i] * n_plus[j] + par_P_plus[j] * n_plus[i]
134  + np_Pp * n_plus[i] * n_plus[j]) / r2_plus
135  + 1.5 * (par_P_minus[i] * n_minus[j] + par_P_minus[j] * n_minus[i]
136  + nm_Pm * n_minus[i] * n_minus[j]) / r2_minus
137  - 3.0 * (np_Sp[i] * n_plus[j] + np_Sp[j] * n_plus[i]) / r3_plus
138  - 3.0 * (nm_Sm[i] * n_minus[j] + nm_Sm[j] * n_minus[i]) / r3_minus;
139  }
140  }
141 
142  Aij[0][0] -= +1.5 * (np_Pp / r2_plus + nm_Pm / r2_minus);
143  Aij[1][1] -= +1.5 * (np_Pp / r2_plus + nm_Pm / r2_minus);
144  Aij[2][2] -= +1.5 * (np_Pp / r2_plus + nm_Pm / r2_minus);
145 
146 
147 }
148 
149 /*-----------------------------------------------------------*/
150 /******** Nonlinear Equations ***********/
151 /*-----------------------------------------------------------*/
152 void
153 TwoPunctures::NonLinEquations (double rho_adm,
154  double A, double B, double X, double R,
155  double x, double r, double phi,
156  double y, double z, derivs U, double *values)
157 {
158 
159  double r_plus, r_minus, psi, psi2, psi4, psi7;
160  double mu;
161 
162  r_plus = sqrt ((x - par_b) * (x - par_b) + y * y + z * z);
163  r_minus = sqrt ((x + par_b) * (x + par_b) + y * y + z * z);
164 
165  psi =
166  1. + 0.5 * par_m_plus / r_plus + 0.5 * par_m_minus / r_minus + U.d0[0];
167  psi2 = psi * psi;
168  psi4 = psi2 * psi2;
169  psi7 = psi * psi2 * psi4;
170 
171  values[0] =
172  U.d11[0] + U.d22[0] + U.d33[0] + 0.125 * BY_KKofxyz (x, y, z) / psi7 +
173  2.0 * Pi / psi2/psi * rho_adm;
174 
175 }
176 
177 /*-----------------------------------------------------------*/
178 /******** Linear Equations ***********/
179 /*-----------------------------------------------------------*/
180 void
181 TwoPunctures::LinEquations (double A, double B, double X, double R,
182  double x, double r, double phi,
183  double y, double z, derivs dU, derivs U, double *values)
184 {
185 
186  double r_plus, r_minus, psi, psi2, psi4, psi8;
187 
188  r_plus = sqrt ((x - par_b) * (x - par_b) + y * y + z * z);
189  r_minus = sqrt ((x + par_b) * (x + par_b) + y * y + z * z);
190 
191  psi =
192  1. + 0.5 * par_m_plus / r_plus + 0.5 * par_m_minus / r_minus + U.d0[0];
193  psi2 = psi * psi;
194  psi4 = psi2 * psi2;
195  psi8 = psi4 * psi4;
196 
197  values[0] = dU.d11[0] + dU.d22[0] + dU.d33[0]
198  - 0.875 * BY_KKofxyz (x, y, z) / psi8 * dU.d0[0];
199 }
200 
201 /*-----------------------------------------------------------*/
202 
203 } // namespace TP
TP::TwoPunctures::BY_KKofxyz
double BY_KKofxyz(double x, double y, double z)
Definition: Equations.cpp:25
TP::TwoPunctures::NonLinEquations
void NonLinEquations(double rho_adm, double A, double B, double X, double R, double x, double r, double phi, double y, double z, derivs U, double *values)
Definition: Equations.cpp:153
TP::TwoPunctures::BY_Aijofxyz
void BY_Aijofxyz(double x, double y, double z, double Aij[3][3])
Definition: Equations.cpp:83
TP::TwoPunctures::LinEquations
void LinEquations(double A, double B, double X, double R, double x, double r, double phi, double y, double z, derivs dU, derivs U, double *values)
Definition: Equations.cpp:181
TP
This file contains aliases for making access to the long state vector Q as used eg.
Definition: CoordTransf.cpp:13
csv_to_plot.y
list y
Definition: csv_to_plot.py:98
euler.j
j
Definition: euler.py:95
TP::DERIVS::d11
double * d11
Definition: TwoPunctures.h:31
TP::DERIVS::d33
double * d33
Definition: TwoPunctures.h:31
TP::DERIVS::d0
double * d0
Definition: TwoPunctures.h:31
TP::DERIVS::d22
double * d22
Definition: TwoPunctures.h:31