Peano 4
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>
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
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
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
268 double psi, psi2, psi4, psi7, r_plus, r_minus;
269 FILE *debugfile = NULL;
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
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
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
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
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
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 be
#define FAC
#define Pi
#define Piq
#define Pih
void free_derivs(derivs *v, int n)
void Derivatives_AB3(int nvar, int n1, int n2, int n3, derivs v)
void JFD_times_dv(int i, int j, int k, int nvar, int n1, int n2, int n3, derivs dv, derivs u, double *values)
void SpecCoef(int n1, int n2, int n3, int ivar, double *v, double *cf)
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)
double PunctEvalAtArbitPositionFasterLowRes(double A, double B, double phi)
void F_of_v(int nvar, int n1, int n2, int n3, derivs v, double *F, derivs u)
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.
void AB_To_XR(int nvar, double A, double B, double *X, double *R, derivs U)
void rx3_To_xyz(int nvar, double x, double r, double phi, double *y, double *z, derivs U)
void calculate_derivs(int i, int j, int k, int ivar, int nvar, int n1, int n2, int n3, derivs v, derivs vv)
double PunctTaylorExpandAtArbitPosition(int ivar, int nvar, int n1, int n2, int n3, derivs v, double x, double y, double z)
double PunctEvalAtArbitPositionFaster(double A, double B, double phi)
void J_times_dv(int nvar, int n1, int n2, int n3, derivs dv, double *Jdv, derivs u)
double PunctIntPolAtArbitPosition(int ivar, int nvar, int n1, int n2, int n3, derivs v, double x, double y, double z)
void SetMatrix_JFD(int nvar, int n1, int n2, int n3, derivs u, int *ncols, int **cols, double **Matrix)
double interpol(double a, double b, double c, derivs v)
void allocate_derivs(derivs *v, int n)
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)
double PunctEvalAtArbitPosition(double *v, int ivar, double A, double B, double phi, int nvar, int n1, int n2, int n3)
int Index(int ivar, int i, int j, int k, int nvar, int n1, int n2, int n3)
double PunctIntPolAtArbitPositionFast(derivs v, double x, double y, double z, bool low_res=false)
void C_To_c(int nvar, double X, double R, double *x, double *r, derivs U)
void fourft(double *u, int N, int inv)
double chebev(double a, double b, double c[], int m, double x)
int maximum3(int i, int j, int k)
int maximum2(int i, int j)
double * dvector(long nl, long nh)
void free_dvector(double *v, long nl, long nh)
void free_ivector(int *v, long nl, long nh)
void chebft_Zeros(double u[], int n, int inv)
void fourder2(double u[], double d2u[], int N)
int minimum2(int i, int j)
void chder(double *c, double *cder, int n)
double *** d3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
int * ivector(long nl, long nh)
void fourder(double u[], double du[], int N)
double fourev(double *u, int N, double x)
double ** dmatrix(long nrl, long nrh, long ncl, long nch)
void free_d3tensor(double ***t, long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
This file contains aliases for making access to the long state vector Q as used eg.
double chebev_wrec(double recvec[], double coeff[], int m)
static double min(double const x, double const y)
void recurrence(double x, int m, double *recvec)
int TP_MyProc()
double * d12
double * d23
double * d1
double * d11
double * d3
double * d33
double * d13
double * d0
double * d22
double * d2
static constexpr int npoints_phi
static constexpr int npoints_A
static constexpr int npoints_B_low
static constexpr int npoints_phi_low
static constexpr int npoints_B
static constexpr int npoints_A_low
bool do_residuum_debug_output