Peano 4
Loading...
Searching...
No Matches
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>
9
10namespace TP {
11using 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
24double
25TwoPunctures::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
82void
83TwoPunctures::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/*-----------------------------------------------------------*/
152void
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/*-----------------------------------------------------------*/
180void
181TwoPunctures::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
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$ our matrix elements will nabla phi_i dx f By this will be a sparse as these basis functions are chosen to not overlap with each other almost everywhere In other they have only local support We can read off the right hand side values
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
#define Pi
double BY_KKofxyz(double x, double y, double z)
Definition Equations.cpp:25
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)
void BY_Aijofxyz(double x, double y, double z, double Aij[3][3])
Definition Equations.cpp:83
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)
This file contains aliases for making access to the long state vector Q as used eg.
double * d11
double * d33
double * d0
double * d22
double par_S_plus[3]
double par_P_plus[3]
double par_S_minus[3]
double par_P_minus[3]