Peano
Loading...
Searching...
No Matches
FuncAndJacobian.cpp
Go to the documentation of this file.
1/* TwoPunctures: File "FuncAndJacobian.c"*/
2
3#include <assert.h>
4#include <stdio.h>
5#include <stdlib.h>
6#include <string.h>
7#include <math.h>
8#include <ctype.h>
9#include <time.h>
10#include "libtwopunctures/TwoPunctures.h"
11#include "libtwopunctures/TP_Parameters.h"
12
13
14namespace TP {
15using namespace Utilities;
16
17#define FAC sin(al)*sin(be)*sin(al)*sin(be)*sin(al)*sin(be)
18/*#define FAC sin(al)*sin(be)*sin(al)*sin(be)*/
19/*#define FAC 1*/
20
21static inline double min (double const x, double const y)
22{
23 return x<y ? x : y;
24}
25
26/* --------------------------------------------------------------------------*/
27int
28TwoPunctures::Index (int ivar, int i, int j, int k, int nvar, int n1, int n2, int n3)
29{
30 int i1 = i, j1 = j, k1 = k;
31
32 if (i1 < 0)
33 i1 = -(i1 + 1);
34 if (i1 >= n1)
35 i1 = 2 * n1 - (i1 + 1);
36
37 if (j1 < 0)
38 j1 = -(j1 + 1);
39 if (j1 >= n2)
40 j1 = 2 * n2 - (j1 + 1);
41
42 if (k1 < 0)
43 k1 = k1 + n3;
44 if (k1 >= n3)
45 k1 = k1 - n3;
46
47 return ivar + nvar * (i1 + n1 * (j1 + n2 * k1));
48}
49
50/* --------------------------------------------------------------------------*/
51void
52TwoPunctures::allocate_derivs (derivs * v, int n)
53{
54 int m = n - 1;
55 (*v).d0 = dvector (0, m);
56 (*v).d1 = dvector (0, m);
57 (*v).d2 = dvector (0, m);
58 (*v).d3 = dvector (0, m);
59 (*v).d11 = dvector (0, m);
60 (*v).d12 = dvector (0, m);
61 (*v).d13 = dvector (0, m);
62 (*v).d22 = dvector (0, m);
63 (*v).d23 = dvector (0, m);
64 (*v).d33 = dvector (0, m);
65}
66
67/* --------------------------------------------------------------------------*/
68void
69TwoPunctures::free_derivs (derivs * v, int n)
70{
71 int m = n - 1;
72 free_dvector ((*v).d0, 0, m);
73 free_dvector ((*v).d1, 0, m);
74 free_dvector ((*v).d2, 0, m);
75 free_dvector ((*v).d3, 0, m);
76 free_dvector ((*v).d11, 0, m);
77 free_dvector ((*v).d12, 0, m);
78 free_dvector ((*v).d13, 0, m);
79 free_dvector ((*v).d22, 0, m);
80 free_dvector ((*v).d23, 0, m);
81 free_dvector ((*v).d33, 0, m);
82}
83
84/* --------------------------------------------------------------------------*/
85void
86TwoPunctures::Derivatives_AB3 (int nvar, int n1, int n2, int n3, derivs v)
87{
88 int i, j, k, ivar, N, *indx;
89 double *p, *dp, *d2p, *q, *dq, *r, *dr;
90
91 N = maximum3 (n1, n2, n3);
92 p = dvector (0, N);
93 dp = dvector (0, N);
94 d2p = dvector (0, N);
95 q = dvector (0, N);
96 dq = dvector (0, N);
97 r = dvector (0, N);
98 dr = dvector (0, N);
99 indx = ivector (0, N);
100
101 for (ivar = 0; ivar < nvar; ivar++)
102 {
103 for (k = 0; k < n3; k++)
104 { /* Calculation of Derivatives w.r.t. A-Dir. */
105 for (j = 0; j < n2; j++)
106 { /* (Chebyshev_Zeros)*/
107 for (i = 0; i < n1; i++)
108 {
109 indx[i] = Index (ivar, i, j, k, nvar, n1, n2, n3);
110 p[i] = v.d0[indx[i]];
111 }
112 chebft_Zeros (p, n1, 0);
113 chder (p, dp, n1);
114 chder (dp, d2p, n1);
115 chebft_Zeros (dp, n1, 1);
116 chebft_Zeros (d2p, n1, 1);
117 for (i = 0; i < n1; i++)
118 {
119 v.d1[indx[i]] = dp[i];
120 v.d11[indx[i]] = d2p[i];
121 }
122 }
123 }
124 for (k = 0; k < n3; k++)
125 { /* Calculation of Derivatives w.r.t. B-Dir. */
126 for (i = 0; i < n1; i++)
127 { /* (Chebyshev_Zeros)*/
128 for (j = 0; j < n2; j++)
129 {
130 indx[j] = Index (ivar, i, j, k, nvar, n1, n2, n3);
131 p[j] = v.d0[indx[j]];
132 q[j] = v.d1[indx[j]];
133 }
134 chebft_Zeros (p, n2, 0);
135 chebft_Zeros (q, n2, 0);
136 chder (p, dp, n2);
137 chder (dp, d2p, n2);
138 chder (q, dq, n2);
139 chebft_Zeros (dp, n2, 1);
140 chebft_Zeros (d2p, n2, 1);
141 chebft_Zeros (dq, n2, 1);
142 for (j = 0; j < n2; j++)
143 {
144 v.d2[indx[j]] = dp[j];
145 v.d22[indx[j]] = d2p[j];
146 v.d12[indx[j]] = dq[j];
147 }
148 }
149 }
150 for (i = 0; i < n1; i++)
151 { /* Calculation of Derivatives w.r.t. phi-Dir. (Fourier)*/
152 for (j = 0; j < n2; j++)
153 {
154 for (k = 0; k < n3; k++)
155 {
156 indx[k] = Index (ivar, i, j, k, nvar, n1, n2, n3);
157 p[k] = v.d0[indx[k]];
158 q[k] = v.d1[indx[k]];
159 r[k] = v.d2[indx[k]];
160 }
161 fourft (p, n3, 0);
162 fourder (p, dp, n3);
163 fourder2 (p, d2p, n3);
164 fourft (dp, n3, 1);
165 fourft (d2p, n3, 1);
166 fourft (q, n3, 0);
167 fourder (q, dq, n3);
168 fourft (dq, n3, 1);
169 fourft (r, n3, 0);
170 fourder (r, dr, n3);
171 fourft (dr, n3, 1);
172 for (k = 0; k < n3; k++)
173 {
174 v.d3[indx[k]] = dp[k];
175 v.d33[indx[k]] = d2p[k];
176 v.d13[indx[k]] = dq[k];
177 v.d23[indx[k]] = dr[k];
178 }
179 }
180 }
181 }
182 free_dvector (p, 0, N);
183 free_dvector (dp, 0, N);
184 free_dvector (d2p, 0, N);
185 free_dvector (q, 0, N);
186 free_dvector (dq, 0, N);
187 free_dvector (r, 0, N);
188 free_dvector (dr, 0, N);
189 free_ivector (indx, 0, N);
190}
191
192/* --------------------------------------------------------------------------*/
193void
194TwoPunctures::F_of_v (int nvar, int n1, int n2, int n3, derivs v, double *F,
195 derivs u)
196{
197 /* Calculates the left hand sides of the non-linear equations F_m(v_n)=0*/
198 /* and the function u (u.d0[]) as well as its derivatives*/
199 /* (u.d1[], u.d2[], u.d3[], u.d11[], u.d12[], u.d13[], u.d22[], u.d23[], u.d33[])*/
200 /* at interior points and at the boundaries "+/-"*/
201
202 int i, j, k, ivar, indx;
203 double al, be, A, B, X, R, x, r, phi, y, z, Am1, *values;
204 derivs U;
205 double *sources;
206
207 values = dvector (0, nvar - 1);
208 allocate_derivs (&U, nvar);
209
210 sources=new double[n1*n2*n3]();
211 if (use_sources)
212 {
213 double *s_x, *s_y, *s_z;
214 int i3D;
215 s_x =new double[n1*n2*n3]();
216 s_y =new double[n1*n2*n3]();
217 s_z =new double[n1*n2*n3]();
218 for (i = 0; i < n1; i++)
219 for (j = 0; j < n2; j++)
220 for (k = 0; k < n3; k++)
221 {
222 i3D = Index(0,i,j,k,1,n1,n2,n3);
223
224 al = Pih * (2 * i + 1) / n1;
225 A = -cos (al);
226 be = Pih * (2 * j + 1) / n2;
227 B = -cos (be);
228 phi = 2. * Pi * k / n3;
229
230 Am1 = A - 1;
231 for (ivar = 0; ivar < nvar; ivar++)
232 {
233 indx = Index (ivar, i, j, k, nvar, n1, n2, n3);
234 U.d0[ivar] = Am1 * v.d0[indx]; /* U*/
235 U.d1[ivar] = v.d0[indx] + Am1 * v.d1[indx]; /* U_A*/
236 U.d2[ivar] = Am1 * v.d2[indx]; /* U_B*/
237 U.d3[ivar] = Am1 * v.d3[indx]; /* U_3*/
238 U.d11[ivar] = 2 * v.d1[indx] + Am1 * v.d11[indx]; /* U_AA*/
239 U.d12[ivar] = v.d2[indx] + Am1 * v.d12[indx]; /* U_AB*/
240 U.d13[ivar] = v.d3[indx] + Am1 * v.d13[indx]; /* U_AB*/
241 U.d22[ivar] = Am1 * v.d22[indx]; /* U_BB*/
242 U.d23[ivar] = Am1 * v.d23[indx]; /* U_B3*/
243 U.d33[ivar] = Am1 * v.d33[indx]; /* U_33*/
244 }
245 /* Calculation of (X,R) and*/
246 /* (U_X, U_R, U_3, U_XX, U_XR, U_X3, U_RR, U_R3, U_33)*/
247 AB_To_XR (nvar, A, B, &X, &R, U);
248 /* Calculation of (x,r) and*/
249 /* (U, U_x, U_r, U_3, U_xx, U_xr, U_x3, U_rr, U_r3, U_33)*/
250 C_To_c (nvar, X, R, &(s_x[i3D]), &r, U);
251 /* Calculation of (y,z) and*/
252 /* (U, U_x, U_y, U_z, U_xx, U_xy, U_xz, U_yy, U_yz, U_zz)*/
253 rx3_To_xyz (nvar, s_x[i3D], r, phi, &(s_y[i3D]), &(s_z[i3D]), U);
254 }
255 // @TODO: Set_Rho_ADM is an external call to another Cactuss thorn or so.
256 /* Set_Rho_ADM(cctkGH, n1*n2*n3, sources, s_x, s_y, s_z); */ // COMMENTED AS WE DO NOT HAVE CCTK HERE
257 free(s_z);
258 free(s_y);
259 free(s_x);
260 }
261 else
262 for (i = 0; i < n1; i++)
263 for (j = 0; j < n2; j++)
264 for (k = 0; k < n3; k++)
265 sources[Index(0,i,j,k,1,n1,n2,n3)]=0.0;
266
267 Derivatives_AB3 (nvar, n1, n2, n3, v);
268 double psi, psi2, psi4, psi7, r_plus, r_minus;
269 FILE *debugfile = NULL;
270 if (do_residuum_debug_output && TP_MyProc() == 0)
271 {
272 debugfile = fopen("res.dat", "w");
273 assert(debugfile);
274 }
275 for (i = 0; i < n1; i++)
276 {
277 for (j = 0; j < n2; j++)
278 {
279 for (k = 0; k < n3; k++)
280 {
281
282 al = Pih * (2 * i + 1) / n1;
283 A = -cos (al);
284 be = Pih * (2 * j + 1) / n2;
285 B = -cos (be);
286 phi = 2. * Pi * k / n3;
287
288 Am1 = A - 1;
289 for (ivar = 0; ivar < nvar; ivar++)
290 {
291 indx = Index (ivar, i, j, k, nvar, n1, n2, n3);
292 U.d0[ivar] = Am1 * v.d0[indx]; /* U*/
293 U.d1[ivar] = v.d0[indx] + Am1 * v.d1[indx]; /* U_A*/
294 U.d2[ivar] = Am1 * v.d2[indx]; /* U_B*/
295 U.d3[ivar] = Am1 * v.d3[indx]; /* U_3*/
296 U.d11[ivar] = 2 * v.d1[indx] + Am1 * v.d11[indx]; /* U_AA*/
297 U.d12[ivar] = v.d2[indx] + Am1 * v.d12[indx]; /* U_AB*/
298 U.d13[ivar] = v.d3[indx] + Am1 * v.d13[indx]; /* U_AB*/
299 U.d22[ivar] = Am1 * v.d22[indx]; /* U_BB*/
300 U.d23[ivar] = Am1 * v.d23[indx]; /* U_B3*/
301 U.d33[ivar] = Am1 * v.d33[indx]; /* U_33*/
302 }
303 /* Calculation of (X,R) and*/
304 /* (U_X, U_R, U_3, U_XX, U_XR, U_X3, U_RR, U_R3, U_33)*/
305 AB_To_XR (nvar, A, B, &X, &R, U);
306 /* Calculation of (x,r) and*/
307 /* (U, U_x, U_r, U_3, U_xx, U_xr, U_x3, U_rr, U_r3, U_33)*/
308 C_To_c (nvar, X, R, &x, &r, U);
309 /* Calculation of (y,z) and*/
310 /* (U, U_x, U_y, U_z, U_xx, U_xy, U_xz, U_yy, U_yz, U_zz)*/
311 rx3_To_xyz (nvar, x, r, phi, &y, &z, U);
312 NonLinEquations (sources[Index(0,i,j,k,1,n1,n2,n3)],
313 A, B, X, R, x, r, phi, y, z, U, values);
314 for (ivar = 0; ivar < nvar; ivar++)
315 {
316 indx = Index (ivar, i, j, k, nvar, n1, n2, n3);
317 F[indx] = values[ivar] * FAC;
318 /* if ((i<5) && ((j<5) || (j>n2-5)))*/
319 /* F[indx] = 0.0;*/
320 u.d0[indx] = U.d0[ivar]; /* U*/
321 u.d1[indx] = U.d1[ivar]; /* U_x*/
322 u.d2[indx] = U.d2[ivar]; /* U_y*/
323 u.d3[indx] = U.d3[ivar]; /* U_z*/
324 u.d11[indx] = U.d11[ivar]; /* U_xx*/
325 u.d12[indx] = U.d12[ivar]; /* U_xy*/
326 u.d13[indx] = U.d13[ivar]; /* U_xz*/
327 u.d22[indx] = U.d22[ivar]; /* U_yy*/
328 u.d23[indx] = U.d23[ivar]; /* U_yz*/
329 u.d33[indx] = U.d33[ivar]; /* U_zz*/
330 }
331 if (debugfile && (k==0))
332 {
333 r_plus = sqrt ((x - par_b) * (x - par_b) + y * y + z * z);
334 r_minus = sqrt ((x + par_b) * (x + par_b) + y * y + z * z);
335 psi = 1.+
336 0.5 * par_m_plus / r_plus +
337 0.5 * par_m_minus / r_minus +
338 U.d0[0];
339 psi2 = psi * psi;
340 psi4 = psi2 * psi2;
341 psi7 = psi * psi2 * psi4;
342 fprintf(debugfile,
343 "%.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g\n",
344 (double)x, (double)y, (double)A, (double)B,
345 (double)(U.d11[0] +
346 U.d22[0] +
347 U.d33[0] +
348/* 0.125 * BY_KKofxyz (x, y, z) / psi7 +*/
349 (2.0 * Pi / psi2/psi * sources[indx]) * FAC),
350 (double)((U.d11[0] +
351 U.d22[0] +
352 U.d33[0])*FAC),
353 (double)(-(2.0 * Pi / psi2/psi * sources[indx]) * FAC),
354 (double)sources[indx]
355 /*(double)F[indx]*/
356 );
357 }
358 }
359 }
360 }
361 if (debugfile)
362 {
363 fclose(debugfile);
364 }
365 free(sources);
366 free_dvector (values, 0, nvar - 1);
367 free_derivs (&U, nvar);
368}
369
370/* --------------------------------------------------------------------------*/
371void
372TwoPunctures::J_times_dv (int nvar, int n1, int n2, int n3, derivs dv,
373 double *Jdv, derivs u)
374{ /* Calculates the left hand sides of the non-linear equations F_m(v_n)=0*/
375 /* and the function u (u.d0[]) as well as its derivatives*/
376 /* (u.d1[], u.d2[], u.d3[], u.d11[], u.d12[], u.d13[], u.d22[], u.d23[], u.d33[])*/
377 /* at interior points and at the boundaries "+/-"*/
378
379 int i, j, k, ivar, indx;
380 double al, be, A, B, X, R, x, r, phi, y, z, Am1, *values;
381 derivs dU, U;
382
383 Derivatives_AB3 (nvar, n1, n2, n3, dv);
384
385
386 for (i = 0; i < n1; i++)
387 {
388 values = dvector (0, nvar - 1);
389 allocate_derivs (&dU, nvar);
390 allocate_derivs (&U, nvar);
391 for (j = 0; j < n2; j++)
392 {
393 for (k = 0; k < n3; k++)
394 {
395
396 al = Pih * (2 * i + 1) / n1;
397 A = -cos (al);
398 be = Pih * (2 * j + 1) / n2;
399 B = -cos (be);
400 phi = 2. * Pi * k / n3;
401
402 Am1 = A - 1;
403 for (ivar = 0; ivar < nvar; ivar++)
404 {
405 indx = Index (ivar, i, j, k, nvar, n1, n2, n3);
406 dU.d0[ivar] = Am1 * dv.d0[indx]; /* dU*/
407 dU.d1[ivar] = dv.d0[indx] + Am1 * dv.d1[indx]; /* dU_A*/
408 dU.d2[ivar] = Am1 * dv.d2[indx]; /* dU_B*/
409 dU.d3[ivar] = Am1 * dv.d3[indx]; /* dU_3*/
410 dU.d11[ivar] = 2 * dv.d1[indx] + Am1 * dv.d11[indx]; /* dU_AA*/
411 dU.d12[ivar] = dv.d2[indx] + Am1 * dv.d12[indx]; /* dU_AB*/
412 dU.d13[ivar] = dv.d3[indx] + Am1 * dv.d13[indx]; /* dU_AB*/
413 dU.d22[ivar] = Am1 * dv.d22[indx]; /* dU_BB*/
414 dU.d23[ivar] = Am1 * dv.d23[indx]; /* dU_B3*/
415 dU.d33[ivar] = Am1 * dv.d33[indx]; /* dU_33*/
416 U.d0[ivar] = u.d0[indx]; /* U */
417 U.d1[ivar] = u.d1[indx]; /* U_x*/
418 U.d2[ivar] = u.d2[indx]; /* U_y*/
419 U.d3[ivar] = u.d3[indx]; /* U_z*/
420 U.d11[ivar] = u.d11[indx]; /* U_xx*/
421 U.d12[ivar] = u.d12[indx]; /* U_xy*/
422 U.d13[ivar] = u.d13[indx]; /* U_xz*/
423 U.d22[ivar] = u.d22[indx]; /* U_yy*/
424 U.d23[ivar] = u.d23[indx]; /* U_yz*/
425 U.d33[ivar] = u.d33[indx]; /* U_zz*/
426 }
427 /* Calculation of (X,R) and*/
428 /* (dU_X, dU_R, dU_3, dU_XX, dU_XR, dU_X3, dU_RR, dU_R3, dU_33)*/
429 AB_To_XR (nvar, A, B, &X, &R, dU);
430 /* Calculation of (x,r) and*/
431 /* (dU, dU_x, dU_r, dU_3, dU_xx, dU_xr, dU_x3, dU_rr, dU_r3, dU_33)*/
432 C_To_c (nvar, X, R, &x, &r, dU);
433 /* Calculation of (y,z) and*/
434 /* (dU, dU_x, dU_y, dU_z, dU_xx, dU_xy, dU_xz, dU_yy, dU_yz, dU_zz)*/
435 rx3_To_xyz (nvar, x, r, phi, &y, &z, dU);
436 LinEquations (A, B, X, R, x, r, phi, y, z, dU, U, values);
437 for (ivar = 0; ivar < nvar; ivar++)
438 {
439 indx = Index (ivar, i, j, k, nvar, n1, n2, n3);
440 Jdv[indx] = values[ivar] * FAC;
441 }
442 }
443 }
444 free_dvector (values, 0, nvar - 1);
445 free_derivs (&dU, nvar);
446 free_derivs (&U, nvar);
447 }
448}
449
450/* --------------------------------------------------------------------------*/
451void
452TwoPunctures::JFD_times_dv (int i, int j, int k, int nvar, int n1, int n2,
453 int n3, derivs dv, derivs u, double *values)
454{ /* Calculates rows of the vector 'J(FD)*dv'.*/
455 /* First row to be calculated: row = Index(0, i, j, k; nvar, n1, n2, n3)*/
456 /* Last row to be calculated: row = Index(nvar-1, i, j, k; nvar, n1, n2, n3)*/
457 /* These rows are stored in the vector JFDdv[0] ... JFDdv[nvar-1].*/
458
459 int ivar, indx;
460 double al, be, A, B, X, R, x, r, phi, y, z, Am1;
461 double sin_al, sin_al_i1, sin_al_i2, sin_al_i3, cos_al;
462 double sin_be, sin_be_i1, sin_be_i2, sin_be_i3, cos_be;
463 double dV0, dV1, dV2, dV3, dV11, dV12, dV13, dV22, dV23, dV33,
464 ha, ga, ga2, hb, gb, gb2, hp, gp, gp2, gagb, gagp, gbgp;
465 derivs dU, U;
466
467 allocate_derivs (&dU, nvar);
468 allocate_derivs (&U, nvar);
469
470 if (k < 0)
471 k = k + n3;
472 if (k >= n3)
473 k = k - n3;
474
475 ha = Pi / n1; /* ha: Stepsize with respect to (al)*/
476 al = ha * (i + 0.5);
477 A = -cos (al);
478 ga = 1 / ha;
479 ga2 = ga * ga;
480
481 hb = Pi / n2; /* hb: Stepsize with respect to (be)*/
482 be = hb * (j + 0.5);
483 B = -cos (be);
484 gb = 1 / hb;
485 gb2 = gb * gb;
486 gagb = ga * gb;
487
488 hp = 2 * Pi / n3; /* hp: Stepsize with respect to (phi)*/
489 phi = hp * k;
490 gp = 1 / hp;
491 gp2 = gp * gp;
492 gagp = ga * gp;
493 gbgp = gb * gp;
494
495
496 sin_al = sin (al);
497 sin_be = sin (be);
498 sin_al_i1 = 1 / sin_al;
499 sin_be_i1 = 1 / sin_be;
500 sin_al_i2 = sin_al_i1 * sin_al_i1;
501 sin_be_i2 = sin_be_i1 * sin_be_i1;
502 sin_al_i3 = sin_al_i1 * sin_al_i2;
503 sin_be_i3 = sin_be_i1 * sin_be_i2;
504 cos_al = -A;
505 cos_be = -B;
506
507 Am1 = A - 1;
508 for (ivar = 0; ivar < nvar; ivar++)
509 {
510 int iccc = Index (ivar, i, j, k, nvar, n1, n2, n3),
511 ipcc = Index (ivar, i + 1, j, k, nvar, n1, n2, n3),
512 imcc = Index (ivar, i - 1, j, k, nvar, n1, n2, n3),
513 icpc = Index (ivar, i, j + 1, k, nvar, n1, n2, n3),
514 icmc = Index (ivar, i, j - 1, k, nvar, n1, n2, n3),
515 iccp = Index (ivar, i, j, k + 1, nvar, n1, n2, n3),
516 iccm = Index (ivar, i, j, k - 1, nvar, n1, n2, n3),
517 icpp = Index (ivar, i, j + 1, k + 1, nvar, n1, n2, n3),
518 icmp = Index (ivar, i, j - 1, k + 1, nvar, n1, n2, n3),
519 icpm = Index (ivar, i, j + 1, k - 1, nvar, n1, n2, n3),
520 icmm = Index (ivar, i, j - 1, k - 1, nvar, n1, n2, n3),
521 ipcp = Index (ivar, i + 1, j, k + 1, nvar, n1, n2, n3),
522 imcp = Index (ivar, i - 1, j, k + 1, nvar, n1, n2, n3),
523 ipcm = Index (ivar, i + 1, j, k - 1, nvar, n1, n2, n3),
524 imcm = Index (ivar, i - 1, j, k - 1, nvar, n1, n2, n3),
525 ippc = Index (ivar, i + 1, j + 1, k, nvar, n1, n2, n3),
526 impc = Index (ivar, i - 1, j + 1, k, nvar, n1, n2, n3),
527 ipmc = Index (ivar, i + 1, j - 1, k, nvar, n1, n2, n3),
528 immc = Index (ivar, i - 1, j - 1, k, nvar, n1, n2, n3);
529 /* Derivatives of (dv) w.r.t. (al,be,phi):*/
530 dV0 = dv.d0[iccc];
531 dV1 = 0.5 * ga * (dv.d0[ipcc] - dv.d0[imcc]);
532 dV2 = 0.5 * gb * (dv.d0[icpc] - dv.d0[icmc]);
533 dV3 = 0.5 * gp * (dv.d0[iccp] - dv.d0[iccm]);
534 dV11 = ga2 * (dv.d0[ipcc] + dv.d0[imcc] - 2 * dv.d0[iccc]);
535 dV22 = gb2 * (dv.d0[icpc] + dv.d0[icmc] - 2 * dv.d0[iccc]);
536 dV33 = gp2 * (dv.d0[iccp] + dv.d0[iccm] - 2 * dv.d0[iccc]);
537 dV12 =
538 0.25 * gagb * (dv.d0[ippc] - dv.d0[ipmc] + dv.d0[immc] - dv.d0[impc]);
539 dV13 =
540 0.25 * gagp * (dv.d0[ipcp] - dv.d0[imcp] + dv.d0[imcm] - dv.d0[ipcm]);
541 dV23 =
542 0.25 * gbgp * (dv.d0[icpp] - dv.d0[icpm] + dv.d0[icmm] - dv.d0[icmp]);
543 /* Derivatives of (dv) w.r.t. (A,B,phi):*/
544 dV11 = sin_al_i3 * (sin_al * dV11 - cos_al * dV1);
545 dV12 = sin_al_i1 * sin_be_i1 * dV12;
546 dV13 = sin_al_i1 * dV13;
547 dV22 = sin_be_i3 * (sin_be * dV22 - cos_be * dV2);
548 dV23 = sin_be_i1 * dV23;
549 dV1 = sin_al_i1 * dV1;
550 dV2 = sin_be_i1 * dV2;
551 /* Derivatives of (dU) w.r.t. (A,B,phi):*/
552 dU.d0[ivar] = Am1 * dV0;
553 dU.d1[ivar] = dV0 + Am1 * dV1;
554 dU.d2[ivar] = Am1 * dV2;
555 dU.d3[ivar] = Am1 * dV3;
556 dU.d11[ivar] = 2 * dV1 + Am1 * dV11;
557 dU.d12[ivar] = dV2 + Am1 * dV12;
558 dU.d13[ivar] = dV3 + Am1 * dV13;
559 dU.d22[ivar] = Am1 * dV22;
560 dU.d23[ivar] = Am1 * dV23;
561 dU.d33[ivar] = Am1 * dV33;
562
563 indx = Index (ivar, i, j, k, nvar, n1, n2, n3);
564 U.d0[ivar] = u.d0[indx]; /* U */
565 U.d1[ivar] = u.d1[indx]; /* U_x*/
566 U.d2[ivar] = u.d2[indx]; /* U_y*/
567 U.d3[ivar] = u.d3[indx]; /* U_z*/
568 U.d11[ivar] = u.d11[indx]; /* U_xx*/
569 U.d12[ivar] = u.d12[indx]; /* U_xy*/
570 U.d13[ivar] = u.d13[indx]; /* U_xz*/
571 U.d22[ivar] = u.d22[indx]; /* U_yy*/
572 U.d23[ivar] = u.d23[indx]; /* U_yz*/
573 U.d33[ivar] = u.d33[indx]; /* U_zz*/
574 }
575 /* Calculation of (X,R) and*/
576 /* (dU_X, dU_R, dU_3, dU_XX, dU_XR, dU_X3, dU_RR, dU_R3, dU_33)*/
577 AB_To_XR (nvar, A, B, &X, &R, dU);
578 /* Calculation of (x,r) and*/
579 /* (dU, dU_x, dU_r, dU_3, dU_xx, dU_xr, dU_x3, dU_rr, dU_r3, dU_33)*/
580 C_To_c (nvar, X, R, &x, &r, dU);
581 /* Calculation of (y,z) and*/
582 /* (dU, dU_x, dU_y, dU_z, dU_xx, dU_xy, dU_xz, dU_yy, dU_yz, dU_zz)*/
583 rx3_To_xyz (nvar, x, r, phi, &y, &z, dU);
584 LinEquations (A, B, X, R, x, r, phi, y, z, dU, U, values);
585 for (ivar = 0; ivar < nvar; ivar++)
586 values[ivar] *= FAC;
587
588 free_derivs (&dU, nvar);
589 free_derivs (&U, nvar);
590}
591
592/* --------------------------------------------------------------------------*/
593void
594TwoPunctures::SetMatrix_JFD (int nvar, int n1, int n2, int n3, derivs u,
595 int *ncols, int **cols, double **Matrix)
596{
597
598 int column, row, mcol;
599 int i, i1, i_0, i_1, j, j1, j_0, j_1, k, k1, k_0, k_1, N1, N2, N3,
600 ivar, ivar1, ntotal = nvar * n1 * n2 * n3;
601 double *values;
602 derivs dv;
603
604 values = dvector (0, nvar - 1);
605 allocate_derivs (&dv, ntotal);
606
607 N1 = n1 - 1;
608 N2 = n2 - 1;
609 N3 = n3 - 1;
610
611 for (i = 0; i < n1; i++)
612 {
613 for (j = 0; j < n2; j++)
614 {
615 for (k = 0; k < n3; k++)
616 {
617 for (ivar = 0; ivar < nvar; ivar++)
618 {
619 row = Index (ivar, i, j, k, nvar, n1, n2, n3);
620 ncols[row] = 0;
621 dv.d0[row] = 0;
622 }
623 }
624 }
625 }
626 for (i = 0; i < n1; i++)
627 {
628 for (j = 0; j < n2; j++)
629 {
630 for (k = 0; k < n3; k++)
631 {
632 for (ivar = 0; ivar < nvar; ivar++)
633 {
634 column = Index (ivar, i, j, k, nvar, n1, n2, n3);
635 dv.d0[column] = 1;
636
637 i_0 = maximum2 (0, i - 1);
638 i_1 = minimum2 (N1, i + 1);
639 j_0 = maximum2 (0, j - 1);
640 j_1 = minimum2 (N2, j + 1);
641 k_0 = k - 1;
642 k_1 = k + 1;
643/* i_0 = 0;
644 i_1 = N1;
645 j_0 = 0;
646 j_1 = N2;
647 k_0 = 0;
648 k_1 = N3;*/
649
650 for (i1 = i_0; i1 <= i_1; i1++)
651 {
652 for (j1 = j_0; j1 <= j_1; j1++)
653 {
654 for (k1 = k_0; k1 <= k_1; k1++)
655 {
656 JFD_times_dv (i1, j1, k1, nvar, n1, n2, n3,
657 dv, u, values);
658 for (ivar1 = 0; ivar1 < nvar; ivar1++)
659 {
660 if (values[ivar1] != 0)
661 {
662 row = Index (ivar1, i1, j1, k1, nvar, n1, n2, n3);
663 mcol = ncols[row];
664 cols[row][mcol] = column;
665 Matrix[row][mcol] = values[ivar1];
666 ncols[row] += 1;
667 }
668 }
669 }
670 }
671 }
672
673 dv.d0[column] = 0;
674 }
675 }
676 }
677 }
678 free_derivs (&dv, ntotal);
679 free_dvector (values, 0, nvar - 1);
680}
681
682/* --------------------------------------------------------------------------*/
683/* Calculates the value of v at an arbitrary position (A,B,phi)*/
684double
685TwoPunctures::PunctEvalAtArbitPosition (double *v, int ivar, double A, double B, double phi,
686 int nvar, int n1, int n2, int n3)
687{
688 int i, j, k, N;
689 double *p, *values1, **values2, result;
690
691 N = maximum3 (n1, n2, n3);
692 p = dvector (0, N);
693 values1 = dvector (0, N);
694 values2 = dmatrix (0, N, 0, N);
695
696 for (k = 0; k < n3; k++)
697 {
698 for (j = 0; j < n2; j++)
699 {
700 for (i = 0; i < n1; i++)
701 p[i] = v[ivar + nvar * (i + n1 * (j + n2 * k))];
702 chebft_Zeros (p, n1, 0);
703 values2[j][k] = chebev (-1, 1, p, n1, A);
704 }
705 }
706
707 for (k = 0; k < n3; k++)
708 {
709 for (j = 0; j < n2; j++)
710 p[j] = values2[j][k];
711 chebft_Zeros (p, n2, 0);
712 values1[k] = chebev (-1, 1, p, n2, B);
713 }
714
715 fourft (values1, n3, 0);
716 result = fourev (values1, n3, phi);
717
718 free_dvector (p, 0, N);
719 free_dvector (values1, 0, N);
720 free_dmatrix (values2, 0, N, 0, N);
721
722 return result;
723}
724
725/* --------------------------------------------------------------------------*/
726void
727TwoPunctures::calculate_derivs (int i, int j, int k, int ivar, int nvar, int n1, int n2,
728 int n3, derivs v, derivs vv)
729{
730 double al = Pih * (2 * i + 1) / n1, be = Pih * (2 * j + 1) / n2,
731 sin_al = sin (al), sin2_al = sin_al * sin_al, cos_al = cos (al),
732 sin_be = sin (be), sin2_be = sin_be * sin_be, cos_be = cos (be);
733
734 vv.d0[0] = v.d0[Index (ivar, i, j, k, nvar, n1, n2, n3)];
735 vv.d1[0] = v.d1[Index (ivar, i, j, k, nvar, n1, n2, n3)] * sin_al;
736 vv.d2[0] = v.d2[Index (ivar, i, j, k, nvar, n1, n2, n3)] * sin_be;
737 vv.d3[0] = v.d3[Index (ivar, i, j, k, nvar, n1, n2, n3)];
738 vv.d11[0] = v.d11[Index (ivar, i, j, k, nvar, n1, n2, n3)] * sin2_al
739 + v.d1[Index (ivar, i, j, k, nvar, n1, n2, n3)] * cos_al;
740 vv.d12[0] =
741 v.d12[Index (ivar, i, j, k, nvar, n1, n2, n3)] * sin_al * sin_be;
742 vv.d13[0] = v.d13[Index (ivar, i, j, k, nvar, n1, n2, n3)] * sin_al;
743 vv.d22[0] = v.d22[Index (ivar, i, j, k, nvar, n1, n2, n3)] * sin2_be
744 + v.d2[Index (ivar, i, j, k, nvar, n1, n2, n3)] * cos_be;
745 vv.d23[0] = v.d23[Index (ivar, i, j, k, nvar, n1, n2, n3)] * sin_be;
746 vv.d33[0] = v.d33[Index (ivar, i, j, k, nvar, n1, n2, n3)];
747}
748
749/* --------------------------------------------------------------------------*/
750double
751TwoPunctures::interpol (double a, double b, double c, derivs v)
752{
753 return v.d0[0]
754 + a * v.d1[0] + b * v.d2[0] + c * v.d3[0]
755 + 0.5 * a * a * v.d11[0] + a * b * v.d12[0] + a * c * v.d13[0]
756 + 0.5 * b * b * v.d22[0] + b * c * v.d23[0] + 0.5 * c * c * v.d33[0];
757}
758
759/* --------------------------------------------------------------------------*/
760/* Calculates the value of v at an arbitrary position (x,y,z)*/
761double
762TwoPunctures::PunctTaylorExpandAtArbitPosition (int ivar, int nvar, int n1,
763 int n2, int n3, derivs v, double x, double y,
764 double z)
765{
766
767 double xs, ys, zs, rs2, phi, X, R, A, B, al, be, aux1, aux2, a, b, c,
768 result, Ui;
769 int i, j, k;
770 derivs vv;
771 allocate_derivs (&vv, 1);
772
773 xs = x / par_b;
774 ys = y / par_b;
775 zs = z / par_b;
776 rs2 = ys * ys + zs * zs;
777 phi = atan2 (z, y);
778 if (phi < 0)
779 phi += 2 * Pi;
780
781 aux1 = 0.5 * (xs * xs + rs2 - 1);
782 aux2 = sqrt (aux1 * aux1 + rs2);
783 X = asinh (sqrt (aux1 + aux2));
784 R = asin (min(1.0, sqrt (-aux1 + aux2)));
785 if (x < 0)
786 R = Pi - R;
787
788 A = 2 * tanh (0.5 * X) - 1;
789 B = tan (0.5 * R - Piq);
790 al = Pi - acos (A);
791 be = Pi - acos (B);
792
793 i = rint (al * n1 / Pi - 0.5);
794 j = rint (be * n2 / Pi - 0.5);
795 k = rint (0.5 * phi * n3 / Pi);
796
797 a = al - Pi * (i + 0.5) / n1;
798 b = be - Pi * (j + 0.5) / n2;
799 c = phi - 2 * Pi * k / n3;
800
801 calculate_derivs (i, j, k, ivar, nvar, n1, n2, n3, v, vv);
802 result = interpol (a, b, c, vv);
803 free_derivs (&vv, 1);
804
805 Ui = (A - 1) * result;
806
807 return Ui;
808}
809
810/* --------------------------------------------------------------------------*/
811/* Calculates the value of v at an arbitrary position (x,y,z)*/
812double
813TwoPunctures::PunctIntPolAtArbitPosition (const int ivar, const int nvar, const int n1,
814 const int n2, const int n3, derivs v, double x, double y,
815 double z)
816{
817
818 double xs, ys, zs, rs2, phi, X, R, A, B, aux1, aux2, result, Ui;
819
820 xs = x / par_b;
821 ys = y / par_b;
822 zs = z / par_b;
823 rs2 = ys * ys + zs * zs;
824 phi = atan2 (z, y);
825 if (phi < 0)
826 phi += 2 * Pi;
827
828 aux1 = 0.5 * (xs * xs + rs2 - 1);
829 aux2 = sqrt (aux1 * aux1 + rs2);
830 X = asinh (sqrt (aux1 + aux2));
831 R = asin (min(1.0, sqrt (-aux1 + aux2)));
832 if (x < 0)
833 R = Pi - R;
834
835 A = 2 * tanh (0.5 * X) - 1;
836 B = tan (0.5 * R - Piq);
837
838 result = PunctEvalAtArbitPosition (v.d0, ivar, A, B, phi, nvar, n1, n2, n3);
839
840 Ui = (A - 1) * result;
841
842 return Ui;
843}
844
845
849
850
851/* Calculates the value of v at an arbitrary position (A,B,phi)* using the fast routine */
852double
853TwoPunctures::PunctEvalAtArbitPositionFast (double *v, const int ivar, double A, double B, double phi, const int nvar, const int n1, const int n2, const int n3)
854{
855 int i, j, k, N;
856 double *p, *values1, **values2, result;
857 // VASILIS: Nothing should be changed in this routine. This is used by PunctIntPolAtArbitPositionFast
858
859 N = maximum3 (n1, n2, n3);
860
861 p = dvector (0, N);
862 values1 = dvector (0, N);
863 values2 = dmatrix (0, N, 0, N);
864
865 for (k = 0; k < n3; k++)
866 {
867 for (j = 0; j < n2; j++)
868 {
869 for (i = 0; i < n1; i++) p[i] = v[ivar + nvar * (i + n1 * (j + n2 * k))];
870 // chebft_Zeros (p, n1, 0);
871 values2[j][k] = chebev (-1, 1, p, n1, A);
872 }
873 }
874
875 for (k = 0; k < n3; k++)
876 {
877 for (j = 0; j < n2; j++) p[j] = values2[j][k];
878 // chebft_Zeros (p, n2, 0);
879 values1[k] = chebev (-1, 1, p, n2, B);
880 }
881
882 // fourft (values1, n3, 0);
883 result = fourev (values1, n3, phi);
884
885 free_dvector (p, 0, N);
886 free_dvector (values1, 0, N);
887 free_dmatrix (values2, 0, N, 0, N);
888
889 return result;
890 // */
891 // return 0.;
892}
893
894
895// Chebychev recurrence, x is expected to be scaled to be within [-1,1]
896void inline recurrence(double x, int m, double * recvec)
897{
898 recvec[0] = 1;
899 recvec[1] = x;
900 for (int i=2;i<m;i++)
901 {
902 recvec[i] = 2*x*recvec[i-1] - recvec[i-2];
903 }
904}
905
906inline double chebev_wrec(double recvec[], double coeff[], int m)
907{
908 double result = -0.5*coeff[0];
909
910#pragma omp simd reduction(+:result)
911 for (int i=0;i<m;i++)
912 {
913 result += recvec[i]*coeff[i];
914 }
915 return result;
916}
917
918
919// Tweaked version NOTE A and B are expected to be in [-1,1] --- seems to be fulfilled
920double
921TwoPunctures::PunctEvalAtArbitPositionFaster (double A, double B, double phi)
922{
923 double p[npoints_B];
924 double values1[npoints_phi];
925 double values2[npoints_B][npoints_phi];
926 double _recvec[npoints_A];
927
928 // Eval and cache recurrence for point A
929 recurrence(A, npoints_A, _recvec);
930
931 // scalar coeff
932 int myglob(0);
933 for (int j = 0; j < npoints_B; j++)
934 for (int k = 0; k < npoints_phi; k++) values2[j][k] = -0.5* _d0contig[myglob++];
935
936 // rest of the chebychev evaluation
937 myglob=0;
938 for (int i = 0; i < npoints_A; i++)
939 for (int j = 0; j < npoints_B; j++)
940#pragma omp simd safelen(4)
941 for (int k = 0; k < npoints_phi; k++)
942 {
943 const int index = k + j*npoints_phi + i*npoints_phi*npoints_B;
944 values2[j][k] += _recvec[i]*_d0contig[index];
945 }
946
947 // Eval and cache recurrence for point B
948 recurrence(B, npoints_A, _recvec);
949
950 // Another chebychev evaluation
951 for (int k = 0; k < npoints_phi; k++)
952 {
953 for (int j = 0; j < npoints_B; j++) p[j] = values2[j][k];
954 values1[k]=chebev_wrec(_recvec, p, npoints_B);
955 }
956
957 double result = fourev (values1, npoints_phi, phi);
958
959
960 return result;
961 // */
962 // return 0.;
963}
964
965 double
966TwoPunctures::PunctEvalAtArbitPositionFasterLowRes (double A, double B, double phi)
967{
968
969 // TODO is it threadsafe to make them members?
970 //double * p = dvector (0, npoints_B_low);
971 //double * values1 = dvector (0, npoints_phi_low);
972 //double ** values2 = dmatrix (0, npoints_B_low, 0, npoints_phi_low);
973 //double * _recvec = dvector (0, npoints_A_low);
974 double p[npoints_B_low];
975 double values1[npoints_phi_low];
976 double values2[npoints_B_low][npoints_phi_low];
977 double _recvec[npoints_A_low];
978
979 // Eval and cache recurrence for point A
980 recurrence(A, npoints_A_low, _recvec);
981
982 // scalar coeff
983 int myglob(0);
984 for (int j = 0; j < npoints_B_low; j++)
985 for (int k = 0; k < npoints_phi_low; k++) values2[j][k] = -0.5* _d0contig_low[myglob++];
986
987 // rest of the chebychev evaluation
988 myglob=0;
989 for (int i = 0; i < npoints_A_low; i++)
990 for (int j = 0; j < npoints_B_low; j++)
991#pragma omp simd safelen(4)
992 for (int k = 0; k < npoints_phi_low; k++)
993 {
994
995 const int index = k + j*n3 + i*n3*n2;
996 //assert(index== myglob++);
997 values2[j][k] += _recvec[i]*_d0contig[index];
998 //values2[j][k] += _recvec[i]*_d0contig_low[myglob++];
999 }
1000
1001 // Eval and cache recurrence for point B
1002 recurrence(B, npoints_A_low, _recvec);
1003
1004 // Another chebychev evaluation
1005 for (int k = 0; k < npoints_phi_low; k++)
1006 {
1007 for (int j = 0; j < npoints_B_low; j++) p[j] = values2[j][k];
1008 values1[k]=chebev_wrec(_recvec, p, npoints_B_low);
1009 }
1010
1011 double result = fourev (values1, npoints_phi_low, phi);
1012
1013 //free_dvector (_recvec, 0, _n1_low);
1014 //free_dvector (p, 0, _n2_low);
1015 //free_dvector (values1, 0, _n3_low);
1016 //free_dmatrix (values2, 0, _n2_low, 0, _n3_low);
1017
1018 return result;
1019 // */
1020 // return 0.;
1021}
1022
1023// --------------------------------------------------------------------------*/
1024// Calculates the value of v at an arbitrary position (x,y,z) if the spectral coefficients are known //
1025/* --------------------------------------------------------------------------*/
1026double
1027TwoPunctures::PunctIntPolAtArbitPositionFast (derivs v, double x, double y, double z, bool low_res)
1028{
1029
1030 double xs, ys, zs, rs2, phi, X, R, A, B, aux1, aux2, result, Ui;
1031 // VASILIS: Here the struct derivs v refers to the spectral coeffiecients of variable v not the variable v itself
1032
1033 xs = x / par_b;
1034 ys = y / par_b;
1035 zs = z / par_b;
1036 rs2 = ys * ys + zs * zs;
1037 phi = atan2 (z, y);
1038 if (phi < 0)
1039 phi += 2 * Pi;
1040
1041 aux1 = 0.5 * (xs * xs + rs2 - 1);
1042 aux2 = sqrt (aux1 * aux1 + rs2);
1043 X = asinh (sqrt (aux1 + aux2));
1044 R = asin (min(1.0, sqrt (-aux1 + aux2)));
1045 if (x < 0)
1046 R = Pi - R;
1047
1048 A = 2 * tanh (0.5 * X) - 1;
1049 B = tan (0.5 * R - Piq);
1050
1051 if (low_res) result = PunctEvalAtArbitPositionFasterLowRes(A, B, phi);
1052 else result = PunctEvalAtArbitPositionFaster (A, B, phi);
1053
1054 Ui = (A - 1) * result;
1055
1056 return Ui;
1057}
1058
1059// Evaluates the spectral expansion coefficients of v
1060void
1061TwoPunctures::SpecCoef (int n1, int n2, int n3, int ivar, double *v, double *cf)
1062{
1063
1064 // VASILIS: Here v is a pointer to the values of the variable v at the collocation points and cf_v a pointer to the spectral coefficients that this routine calculates
1065
1066 int i, j, k, N, n, l, m;
1067 double *p, ***values3, ***values4;
1068
1069 N=maximum3(n1,n2,n3);
1070 p=dvector(0,N);
1071 values3=d3tensor(0,n1,0,n2,0,n3);
1072 values4=d3tensor(0,n1,0,n2,0,n3);
1073
1074
1075
1076 // Caclulate values3[n,j,k] = a_n^{j,k} = (sum_i^(n1-1) f(A_i,B_j,phi_k) Tn(-A_i))/k_n , k_n = N/2 or N
1077 for(k=0;k<n3;k++) {
1078 for(j=0;j<n2;j++) {
1079
1080 for(i=0;i<n1;i++) p[i]=v[ivar + (i + n1 * (j + n2 * k))];
1081
1082 chebft_Zeros(p,n1,0);
1083 for (n=0;n<n1;n++) {
1084 values3[n][j][k] = p[n];
1085 }
1086 }
1087 }
1088
1089 // Caclulate values4[n,l,k] = a_{n,l}^{k} = (sum_j^(n2-1) a_n^{j,k} Tn(B_j))/k_l , k_l = N/2 or N
1090
1091 for (n = 0; n < n1; n++){
1092 for(k=0;k<n3;k++) {
1093 for(j=0;j<n2;j++) p[j]=values3[n][j][k];
1094 chebft_Zeros(p,n2,0);
1095 for (l = 0; l < n2; l++){
1096 values4[n][l][k] = p[l];
1097 }
1098 }
1099 }
1100
1101 // Caclulate coefficients a_{n,l,m} = (sum_k^(n3-1) a_{n,m}^{k} fourier(phi_k))/k_m , k_m = N/2 or N
1102 for (i = 0; i < n1; i++){
1103 for (j = 0; j < n2; j++){
1104 for(k=0;k<n3;k++) p[k]=values4[i][j][k];
1105 fourft(p,n3,0);
1106 for (k = 0; k<n3; k++){
1107 cf[ivar + (i + n1 * (j + n2 * k))] = p[k];
1108 }
1109 }
1110 }
1111
1112 free_dvector(p,0,N);
1113 free_d3tensor(values3,0,n1,0,n2,0,n3);
1114 free_d3tensor(values4,0,n1,0,n2,0,n3);
1115
1116}
1117
1118} // namespace TP
TP_Parameters.h
TP::TwoPunctures::free_derivs
void free_derivs(derivs *v, int n)
Definition FuncAndJacobian.cpp:69
TP::TwoPunctures::Derivatives_AB3
void Derivatives_AB3(int nvar, int n1, int n2, int n3, derivs v)
Definition FuncAndJacobian.cpp:86
TP::TwoPunctures::JFD_times_dv
void JFD_times_dv(int i, int j, int k, int nvar, int n1, int n2, int n3, derivs dv, derivs u, double *values)
Definition FuncAndJacobian.cpp:452
TP::TwoPunctures::SpecCoef
void SpecCoef(int n1, int n2, int n3, int ivar, double *v, double *cf)
Definition FuncAndJacobian.cpp:1061
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::PunctEvalAtArbitPositionFasterLowRes
double PunctEvalAtArbitPositionFasterLowRes(double A, double B, double phi)
Definition FuncAndJacobian.cpp:966
TP::TwoPunctures::F_of_v
void F_of_v(int nvar, int n1, int n2, int n3, derivs v, double *F, derivs u)
Definition FuncAndJacobian.cpp:194
TP::TwoPunctures::PunctEvalAtArbitPositionFast
double PunctEvalAtArbitPositionFast(double *v, int ivar, double A, double B, double phi, int nvar, int n1, int n2, int n3)
Fast Spectral Interpolation Routine Stuff.
Definition FuncAndJacobian.cpp:853
TP::TwoPunctures::AB_To_XR
void AB_To_XR(int nvar, double A, double B, double *X, double *R, derivs U)
Definition CoordTransf.cpp:18
TP::TwoPunctures::rx3_To_xyz
void rx3_To_xyz(int nvar, double x, double r, double phi, double *y, double *z, derivs U)
Definition CoordTransf.cpp:116
TP::TwoPunctures::calculate_derivs
void calculate_derivs(int i, int j, int k, int ivar, int nvar, int n1, int n2, int n3, derivs v, derivs vv)
Definition FuncAndJacobian.cpp:727
TP::TwoPunctures::PunctTaylorExpandAtArbitPosition
double PunctTaylorExpandAtArbitPosition(int ivar, int nvar, int n1, int n2, int n3, derivs v, double x, double y, double z)
Definition FuncAndJacobian.cpp:762
TP::TwoPunctures::PunctEvalAtArbitPositionFaster
double PunctEvalAtArbitPositionFaster(double A, double B, double phi)
Definition FuncAndJacobian.cpp:921
TP::TwoPunctures::J_times_dv
void J_times_dv(int nvar, int n1, int n2, int n3, derivs dv, double *Jdv, derivs u)
Definition FuncAndJacobian.cpp:372
TP::TwoPunctures::PunctIntPolAtArbitPosition
double PunctIntPolAtArbitPosition(int ivar, int nvar, int n1, int n2, int n3, derivs v, double x, double y, double z)
Definition FuncAndJacobian.cpp:813
TP::TwoPunctures::SetMatrix_JFD
void SetMatrix_JFD(int nvar, int n1, int n2, int n3, derivs u, int *ncols, int **cols, double **Matrix)
Definition FuncAndJacobian.cpp:594
TP::TwoPunctures::interpol
double interpol(double a, double b, double c, derivs v)
Definition FuncAndJacobian.cpp:751
TP::TwoPunctures::allocate_derivs
void allocate_derivs(derivs *v, int n)
Definition FuncAndJacobian.cpp:52
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::TwoPunctures::PunctEvalAtArbitPosition
double PunctEvalAtArbitPosition(double *v, int ivar, double A, double B, double phi, int nvar, int n1, int n2, int n3)
Definition FuncAndJacobian.cpp:685
TP::TwoPunctures::Index
int Index(int ivar, int i, int j, int k, int nvar, int n1, int n2, int n3)
Definition FuncAndJacobian.cpp:28
TP::TwoPunctures::PunctIntPolAtArbitPositionFast
double PunctIntPolAtArbitPositionFast(derivs v, double x, double y, double z, bool low_res=false)
Definition FuncAndJacobian.cpp:1027
TP::TwoPunctures::C_To_c
void C_To_c(int nvar, double X, double R, double *x, double *r, derivs U)
Definition CoordTransf.cpp:53
TP::Utilities::fourft
void fourft(double *u, int N, int inv)
Definition TP_Utilities.cpp:378
TP::Utilities::chebev
double chebev(double a, double b, double c[], int m, double x)
Definition TP_Utilities.cpp:354
TP::Utilities::maximum3
int maximum3(int i, int j, int k)
Definition TP_Utilities.cpp:226
TP::Utilities::maximum2
int maximum2(int i, int j)
Definition TP_Utilities.cpp:216
TP::Utilities::dvector
double * dvector(long nl, long nh)
Definition TP_Utilities.cpp:30
TP::Utilities::free_dvector
void free_dvector(double *v, long nl, long nh)
Definition TP_Utilities.cpp:157
TP::Utilities::free_ivector
void free_ivector(int *v, long nl, long nh)
Definition TP_Utilities.cpp:149
TP::Utilities::chebft_Zeros
void chebft_Zeros(double u[], int n, int inv)
Definition TP_Utilities.cpp:250
TP::Utilities::fourder2
void fourder2(double u[], double d2u[], int N)
Definition TP_Utilities.cpp:458
TP::Utilities::minimum2
int minimum2(int i, int j)
Definition TP_Utilities.cpp:194
TP::Utilities::chder
void chder(double *c, double *cder, int n)
Definition TP_Utilities.cpp:342
TP::Utilities::d3tensor
double *** d3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
Definition TP_Utilities.cpp:102
TP::Utilities::ivector
int * ivector(long nl, long nh)
Definition TP_Utilities.cpp:16
TP::Utilities::fourder
void fourder(double u[], double du[], int N)
Definition TP_Utilities.cpp:441
TP::Utilities::fourev
double fourev(double *u, int N, double x)
Definition TP_Utilities.cpp:476
TP::Utilities::dmatrix
double ** dmatrix(long nrl, long nrh, long ncl, long nch)
Definition TP_Utilities.cpp:73
TP::Utilities::free_d3tensor
void free_d3tensor(double ***t, long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
Definition TP_Utilities.cpp:183
TP::Utilities::free_dmatrix
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
Definition TP_Utilities.cpp:174
TP
This file contains aliases for making access to the long state vector Q as used eg.
Definition CoordTransf.cpp:13
TP::chebev_wrec
double chebev_wrec(double recvec[], double coeff[], int m)
Definition FuncAndJacobian.cpp:906
TP::derivs
struct TP::DERIVS derivs
TP::min
static double min(double const x, double const y)
Definition FuncAndJacobian.cpp:21
TP::recurrence
void recurrence(double x, int m, double *recvec)
Definition FuncAndJacobian.cpp:896
TP::TP_MyProc
int TP_MyProc()
Definition TwoPunctures.h:24
TP::DERIVS::d12
double * d12
Definition TwoPunctures.h:31
TP::DERIVS::d23
double * d23
Definition TwoPunctures.h:31
TP::DERIVS::d1
double * d1
Definition TwoPunctures.h:31
TP::DERIVS::d11
double * d11
Definition TwoPunctures.h:31
TP::DERIVS::d3
double * d3
Definition TwoPunctures.h:31
TP::DERIVS::d33
double * d33
Definition TwoPunctures.h:31
TP::DERIVS::d13
double * d13
Definition TwoPunctures.h:31
TP::DERIVS::d0
double * d0
Definition TwoPunctures.h:31
TP::DERIVS::d22
double * d22
Definition TwoPunctures.h:31
TP::DERIVS::d2
double * d2
Definition TwoPunctures.h:31
TP::Parameters::npoints_phi
static constexpr int npoints_phi
Definition TP_Parameters.h:28
TP::Parameters::npoints_A
static constexpr int npoints_A
Definition TP_Parameters.h:26
TP::Parameters::npoints_B_low
static constexpr int npoints_B_low
Definition TP_Parameters.h:30
TP::Parameters::npoints_phi_low
static constexpr int npoints_phi_low
Definition TP_Parameters.h:31
TP::Parameters::npoints_B
static constexpr int npoints_B
Definition TP_Parameters.h:27
TP::Parameters::npoints_A_low
static constexpr int npoints_A_low
Definition TP_Parameters.h:29
TP::Parameters::do_residuum_debug_output
bool do_residuum_debug_output
Definition TP_Parameters.h:54