Peano 4
Loading...
Searching...
No Matches
TwoPunctures.cpp
Go to the documentation of this file.
1/* TwoPunctures: File "TwoPunctures.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>
10
11namespace TP {
12using namespace Utilities;
13using namespace Z4VectorShortcuts;
14
15void
17{
18 int nvar = 1;
19
20 double *s_x, *s_y, *s_z;
21 double al, A, Am1, be, B, phi, R, r, X;
22 int ivar, i, j, k, i3D, indx;
23 derivs U;
24 FILE *debug_file;
25
27 nvar = 4;
28
29 s_x =new double[n1*n2*n3]();
30 s_y =new double[n1*n2*n3]();
31 s_z =new double[n1*n2*n3]();
33 for (ivar = 0; ivar < nvar; ivar++)
34 for (i = 0; i < n1; i++)
35 for (j = 0; j < n2; j++)
36 for (k = 0; k < n3; k++)
37 {
38 i3D = Index(ivar,i,j,k,1,n1,n2,n3);
39
40 al = Pih * (2 * i + 1) / n1;
41 A = -cos (al);
42 be = Pih * (2 * j + 1) / n2;
43 B = -cos (be);
44 phi = 2. * Pi * k / n3;
45
46 /* Calculation of (X,R)*/
47 AB_To_XR (nvar, A, B, &X, &R, U);
48 /* Calculation of (x,r)*/
49 C_To_c (nvar, X, R, &(s_x[i3D]), &r, U);
50 /* Calculation of (y,z)*/
51 rx3_To_xyz (nvar, s_x[i3D], r, phi, &(s_y[i3D]), &(s_z[i3D]), U);
52 }
53 // @TODO where is the function Set_Initial_Guess_for_u implemented!???
54 // When looking in EinsteinInitialData, I find it at
55 /*
56TOVSolver/interface.ccl
5735:CCTK_INT FUNCTION Set_Initial_Guess_for_u( \
5843:PROVIDES FUNCTION Set_Initial_Guess_for_u \
5944: WITH TOV_Set_Initial_Guess_for_u \
60 */
62 for (ivar = 0; ivar < nvar; ivar++)
63 for (i = 0; i < n1; i++)
64 for (j = 0; j < n2; j++)
65 for (k = 0; k < n3; k++)
66 {
67 indx = Index(ivar,i,j,k,1,n1,n2,n3);
68 v.d0[indx]/=(-cos(Pih * (2 * i + 1) / n1)-1.0);
69 }
72 {
73 debug_file=fopen("initial.dat", "w");
74 assert(debug_file);
75 for (ivar = 0; ivar < nvar; ivar++)
76 for (i = 0; i < n1; i++)
77 for (j = 0; j < n2; j++)
78 {
79 al = Pih * (2 * i + 1) / n1;
80 A = -cos (al);
81 Am1 = A -1.0;
82 be = Pih * (2 * j + 1) / n2;
83 B = -cos (be);
84 phi = 0.0;
85 indx = Index(ivar,i,j,0,1,n1,n2,n3);
86 U.d0[0] = Am1 * v.d0[indx]; /* U*/
87 U.d1[0] = v.d0[indx] + Am1 * v.d1[indx]; /* U_A*/
88 U.d2[0] = Am1 * v.d2[indx]; /* U_B*/
89 U.d3[0] = Am1 * v.d3[indx]; /* U_3*/
90 U.d11[0] = 2 * v.d1[indx] + Am1 * v.d11[indx]; /* U_AA*/
91 U.d12[0] = v.d2[indx] + Am1 * v.d12[indx]; /* U_AB*/
92 U.d13[0] = v.d3[indx] + Am1 * v.d13[indx]; /* U_AB*/
93 U.d22[0] = Am1 * v.d22[indx]; /* U_BB*/
94 U.d23[0] = Am1 * v.d23[indx]; /* U_B3*/
95 U.d33[0] = Am1 * v.d33[indx]; /* U_33*/
96 /* Calculation of (X,R)*/
97 AB_To_XR (nvar, A, B, &X, &R, U);
98 /* Calculation of (x,r)*/
99 C_To_c (nvar, X, R, &(s_x[indx]), &r, U);
100 /* Calculation of (y,z)*/
101 rx3_To_xyz (nvar, s_x[i3D], r, phi, &(s_y[indx]), &(s_z[indx]), U);
102 fprintf(debug_file,
103 "%.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g "
104 "%.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g\n",
105 (double)s_x[indx], (double)s_y[indx],
106 (double)A,(double)B,
107 (double)U.d0[0],
108 (double)(-cos(Pih * (2 * i + 1) / n1)-1.0),
109 (double)U.d1[0],
110 (double)U.d2[0],
111 (double)U.d3[0],
112 (double)U.d11[0],
113 (double)U.d22[0],
114 (double)U.d33[0],
115 (double)v.d0[indx],
116 (double)v.d1[indx],
117 (double)v.d2[indx],
118 (double)v.d3[indx],
119 (double)v.d11[indx],
120 (double)v.d22[indx],
121 (double)v.d33[indx]
122 );
123 }
124 fprintf(debug_file, "\n\n");
125 for (i=n2-10; i<n2; i++)
126 {
127 double d;
128 indx = Index(0,0,i,0,1,n1,n2,n3);
130 s_x[indx], 0.0, 0.0);
131 fprintf(debug_file, "%.16g %.16g\n",
132 (double)s_x[indx], (double)d);
133 }
134 fprintf(debug_file, "\n\n");
135 for (i=n2-10; i<n2-1; i++)
136 {
137 double d;
138 int ip= Index(0,0,i+1,0,1,n1,n2,n3);
139 indx = Index(0,0,i,0,1,n1,n2,n3);
140 for (j=-10; j<10; j++)
141 {
143 s_x[indx]+(s_x[ip]-s_x[indx])*j/10,
144 0.0, 0.0);
145 fprintf(debug_file, "%.16g %.16g\n",
146 (double)(s_x[indx]+(s_x[ip]-s_x[indx])*j/10), (double)d);
147 }
148 }
149 fprintf(debug_file, "\n\n");
150 for (i = 0; i < n1; i++)
151 for (j = 0; j < n2; j++)
152 {
153 X = 2*(2.0*i/n1-1.0);
154 R = 2*(1.0*j/n2);
155 if (X*X+R*R > 1.0)
156 {
157 C_To_c (nvar, X, R, &(s_x[indx]), &r, U);
158 rx3_To_xyz (nvar, s_x[i3D], r, phi, &(s_y[indx]), &(s_z[indx]), U);
159 *U.d0 = s_x[indx]*s_x[indx];
160 *U.d1 = 2*s_x[indx];
161 *U.d2 = 0.0;
162 *U.d3 = 0.0;
163 *U.d11 = 2.0;
164 *U.d22 = 0.0;
165 *U.d33 = *U.d12 = *U.d23 = *U.d13 = 0.0;
166 C_To_c (nvar, X, R, &(s_x[indx]), &r, U);
167 fprintf(debug_file,
168 "%.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g\n",
169 (double)s_x[indx], (double)r, (double)X, (double)R, (double)U.d0[0],
170 (double)U.d1[0],
171 (double)U.d2[0],
172 (double)U.d3[0],
173 (double)U.d11[0],
174 (double)U.d22[0],
175 (double)U.d33[0]);
176 }
177 }
178 fclose(debug_file);
179 }
180 free(s_z);
181 free(s_y);
182 free(s_x);
183 free_derivs (&U, nvar);
184 /*exit(0);*/
185}
186
187/* -------------------------------------------------------------------*/
188void
190{
191
194
195 mp = par_m_plus;
196 mm = par_m_minus;
197
198 int imin[3], imax[3];
199 int const ntotal = n1 * n2 * n3 * nvar;
200#if 0
201 int percent10 = 0;
202#endif
203 static double *F = NULL;
204 double admMass;
205
206 if (! F) {
207 double up, um;
208 /* Solve only when called for the first time */
209 F = dvector (0, ntotal - 1);
210 allocate_derivs (&u, ntotal);
211 allocate_derivs (&v, ntotal);
212 allocate_derivs (&cf_v, ntotal);
213
214 if (use_sources) {
215 TP_INFO ("Solving puncture equation for BH-NS/NS-NS system");
216 } else {
217 TP_INFO ("Solving puncture equation for BH-BH system");
218 }
219 TP_INFO ( "b = %g", par_b);
220
221 /* initialise to 0 */
222 for (int j = 0; j < ntotal; j++)
223 {
224 cf_v.d0[j] = 0.0;
225 cf_v.d1[j] = 0.0;
226 cf_v.d2[j] = 0.0;
227 cf_v.d3[j] = 0.0;
228 cf_v.d11[j] = 0.0;
229 cf_v.d12[j] = 0.0;
230 cf_v.d13[j] = 0.0;
231 cf_v.d22[j] = 0.0;
232 cf_v.d23[j] = 0.0;
233 cf_v.d33[j] = 0.0;
234 v.d0[j] = 0.0;
235 v.d1[j] = 0.0;
236 v.d2[j] = 0.0;
237 v.d3[j] = 0.0;
238 v.d11[j] = 0.0;
239 v.d12[j] = 0.0;
240 v.d13[j] = 0.0;
241 v.d22[j] = 0.0;
242 v.d23[j] = 0.0;
243 v.d33[j] = 0.0;
244 }
245 /* call for external initial guess */
247 {
249 }
250
251 /* If bare masses are not given, iteratively solve for them given the
252 target ADM masses target_M_plus and target_M_minus and with initial
253 guesses given by par_m_plus and par_m_minus. */
254 if(!(give_bare_mass)) {
255 double tmp, mp_adm_err, mm_adm_err;
256 char valbuf[100];
257
258 double M_p = target_M_plus;
259 double M_m = target_M_minus;
260
261 TP_INFO ( "Attempting to find bare masses.");
262 TP_INFO ( "Target ADM masses: M_p=%g and M_m=%g",
263 (double) M_p, (double) M_m);
264 TP_INFO ( "ADM mass tolerance: %g", (double) adm_tol);
265
266 /* Loop until both ADM masses are within adm_tol of their target */
267 do {
268 TP_INFO ( "Bare masses: mp=%.15g, mm=%.15g",
269 (double)mp, (double)mm);
270 Newton (nvar, n1, n2, n3, v, Newton_tol, 1);
271
272 F_of_v (nvar, n1, n2, n3, v, F, u);
273
274 up = PunctIntPolAtArbitPosition(0, nvar, n1, n2, n3, v, par_b, 0., 0.);
275 um = PunctIntPolAtArbitPosition(0, nvar, n1, n2, n3, v,-par_b, 0., 0.);
276
277 /* Calculate the ADM masses from the current bare mass guess */
278 mp_adm = (1 + up) * mp + mp * mm / (4. * par_b);
279 mm_adm = (1 + um) * mm + mp * mm / (4. * par_b);
280
281 /* Check how far the current ADM masses are from the target */
282 mp_adm_err = fabs(M_p-mp_adm);
283 mm_adm_err = fabs(M_m-mm_adm);
284 TP_INFO ( "ADM mass error: M_p_err=%.15g, M_m_err=%.15g",
285 (double)mp_adm_err, (double)mm_adm_err);
286
287 /* Invert the ADM mass equation and update the bare mass guess so that
288 it gives the correct target ADM masses */
289 tmp = -4*par_b*( 1 + um + up + um*up ) +
290 sqrt(16*par_b*M_m*(1 + um)*(1 + up) +
291 pow(-M_m + M_p + 4*par_b*(1 + um)*(1 + up),2));
292 mp = (tmp + M_p - M_m)/(2.*(1 + up));
293 mm = (tmp - M_p + M_m)/(2.*(1 + um));
294
295 /* Set the par_m_plus and par_m_minus parameters */
296 par_m_plus = mp;
297 par_m_minus = mm;
298
299 } while ( (mp_adm_err > adm_tol) ||
300 (mm_adm_err > adm_tol) );
301
302 TP_INFO ( "Found bare masses.");
303 } else {
304 TP_INFO ( "Use input bare masses." );
305 }
306
308
309 F_of_v (nvar, n1, n2, n3, v, F, u);
310
311 SpecCoef(n1, n2, n3, 0, v.d0, cf_v.d0);
312
313 // TODO destructor
314 // HS copy d0 for contiguous access --- note this is valied for nvar=1
315 _d0contig = Utilities::dvector(0,n1*n2*n3);
316 int ctr(0);
317 for (int i = 0; i < n1; i++)
318 for (int j = 0; j < n2; j++)
319 for (int k = 0; k < n3; k++) _d0contig[ctr++] = cf_v.d0[i + n1 * (j + n2 * k)];
320
321 // And the low res version --- note how we truncate coefficients here
322 // --- the loops use the truncated indices while the access uses the full info n1 and n2
323 _d0contig_low = Utilities::dvector(0,_n1_low*_n2_low*_n3_low);
324 ctr = 0;
325 for (int i = 0; i < _n1_low; i++)
326 for (int j = 0; j < _n2_low; j++)
327 for (int k = 0; k < _n3_low; k++) _d0contig_low[ctr++] = cf_v.d0[i + n1 * (j + n2 * k)];
328
329 TP_INFO (
330 "The two puncture masses are mp=%.17g and mm=%.17g",
331 (double) mp, (double) mm);
332
333 up = PunctIntPolAtArbitPosition(0, nvar, n1, n2, n3, v, par_b, 0., 0.);
334 um = PunctIntPolAtArbitPosition(0, nvar, n1, n2, n3, v,-par_b, 0., 0.);
335
336 /* Calculate the ADM masses from the current bare mass guess */
337 mp_adm = (1 + up) * mp + mp * mm / (4. * par_b);
338 mm_adm = (1 + um) * mm + mp * mm / (4. * par_b);
339
340 TP_INFO ( "Puncture 1 ADM mass is %g", (double) mp_adm);
341 TP_INFO ( "Puncture 2 ADM mass is %g", (double) mm_adm);
342
343 /* print out ADM mass, eq.: \Delta M_ADM=2*r*u=4*b*V for A=1,B=0,phi=0 */
344 admMass = (mp + mm
345 - 4*par_b*PunctEvalAtArbitPosition(v.d0, 0, 1, 0, 0, nvar, n1, n2, n3));
346 TP_INFO ( "The total ADM mass is %g", (double) admMass);
347
348 /*
349 Run this in Mathematica (version 8 or later) with
350 math -script <file>
351
352 Needs["SymbolicC`"];
353 co = Table["center_offset[" <> ToString[i] <> "]", {i, 0, 2}];
354 r1 = co + {"par_b", 0, 0};
355 r2 = co + {-"par_b", 0, 0};
356 {p1, p2} = Table["par_P_" <> bh <> "[" <> ToString[i] <> "]", {bh, {"plus", "minus"}}, {i, 0, 2}];
357 {s1, s2} = Table["par_S_" <> bh <> "[" <> ToString[i] <> "]", {bh, {"plus", "minus"}}, {i, 0, 2}];
358
359 J = Cross[r1, p1] + Cross[r2, p2] + s1 + s2;
360
361 JVar = Table["*J" <> ToString[i], {i, 1, 3}];
362 Print[OutputForm@StringReplace[
363 ToCCodeString@MapThread[CAssign[#1, CExpression[#2]] &, {JVar, J}],
364 "\"" -> ""]];
365 */
366
370 }
371
372 if (grid_setup_method == "Taylor expansion")
373 {
375 }
376 else if (grid_setup_method == "evaluation")
377 {
379 }
380 else
381 {
382 TP_ERROR ("Illegal value for grid_setup_method.");
383 std::abort();
384 }
385
386 if(initial_lapse == "twopunctures")
387 TP_ERROR("Please specify a lapse which we can use");
388 antisymmetric_lapse = (initial_lapse == "twopunctures-antisymmetric");
389 averaged_lapse = (initial_lapse == "twopunctures-averaged");
390 pmn_lapse = (initial_lapse == "psi^n");
391 if (pmn_lapse)
392 TP_INFO ("Setting initial lapse to psi^%f profile.",
394 brownsville_lapse = (initial_lapse == "brownsville");
396 TP_INFO ( "Setting initial lapse to a Brownsville-style profile "
397 "with exp %f.",
399
400 TP_INFO ("Preparing interpolation of result");
401 if (metric_type == "static conformal") {
402 if (conformal_storage == "factor") {
403 conformal_state = 1;
404 } else if (conformal_storage == "factor+derivs") {
405 conformal_state = 2;
406 } else if (conformal_storage == "factor+derivs+2nd derivs") {
407 conformal_state = 3;
408 }
409 } else {
410 conformal_state = 0;
411 }
412
413 runned = true;
414} /* End of TwoPointures_Setup() */
415
416
420void
421TwoPunctures::Interpolate (const double* const pos, double *Q, bool low_res) {
422 if(!runned) {
423 TP_ERROR ("You must call Run() before interpolating.");
424 std::abort();
425 }
426
427 double x=pos[0], y=pos[1], z=pos[2];
428 double xx, yy, zz;
429 xx = x - center_offset[0];
430 yy = y - center_offset[1];
431 zz = z - center_offset[2];
432
433 /* We implement swapping the x and z coordinates as follows.
434 The bulk of the code that performs the actual calculations
435 is unchanged. This code looks only at local variables.
436 Before the bulk --i.e., here-- we swap all x and z tensor
437 components, and after the code --i.e., at the end of this
438 main loop-- we swap everything back. */
439 if (swap_xz) {
440 /* Swap the x and z coordinates */
441 SWAP (xx, zz);
442 }
443
444 const double TP4 = TP_epsilon*TP_epsilon*TP_epsilon*TP_epsilon;
445
446 double r_plus
447 = sqrt((xx - par_b)*(xx - par_b) + yy*yy + zz*zz);
448 double r_minus
449 = sqrt((xx + par_b)*(xx + par_b) + yy*yy + zz*zz);
450
451 double U;
452 switch (gsm)
453 {
456 (0, nvar, n1, n2, n3, v, xx, yy, zz);
457 break;
458 case GSM_evaluation:
459 U = PunctIntPolAtArbitPositionFast(cf_v, xx, yy, zz, low_res);
460 break;
461 default:
462 assert (0);
463 }
464 r_plus = pow (pow (r_plus, 4) + TP4, 0.25);
465 r_minus = pow (pow (r_minus, 4) + TP4, 0.25);
466 if (r_plus < TP_Tiny)
467 r_plus = TP_Tiny;
468 if (r_minus < TP_Tiny)
469 r_minus = TP_Tiny;
470 double psi1 = 1
471 + 0.5 * mp / r_plus
472 + 0.5 * mm / r_minus + U;
473#define EXTEND(M,r) \
474 ( M * (3./8 * pow(r, 4) / pow(TP_Extend_Radius, 5) - \
475 5./4 * pow(r, 2) / pow(TP_Extend_Radius, 3) + \
476 15./8 / TP_Extend_Radius))
477 if (r_plus < TP_Extend_Radius) {
478 psi1 = 1
479 + 0.5 * EXTEND(mp,r_plus)
480 + 0.5 * mm / r_minus + U;
481 }
482 if (r_minus < TP_Extend_Radius) {
483 psi1 = 1
484 + 0.5 * EXTEND(mm,r_minus)
485 + 0.5 * mp / r_plus + U;
486 }
487 double static_psi = 1;
488
489 double Aij[3][3];
490 BY_Aijofxyz (xx, yy, zz, Aij);
491
492 double old_alp=1.0;
494 TP_WARN("@todo: Old Lapse is not be possible, Q as input vector is undefined.");
495
496
497 if ((conformal_state > 0) || (pmn_lapse) || (brownsville_lapse)) {
498
499 double xp, yp, zp, rp, ir;
500 double s1, s3, s5;
501 double p, px, py, pz, pxx, pxy, pxz, pyy, pyz, pzz;
502 p = 1.0;
503 px = py = pz = 0.0;
504 pxx = pxy = pxz = 0.0;
505 pyy = pyz = pzz = 0.0;
506
507 /* first puncture */
508 xp = xx - par_b;
509 yp = yy;
510 zp = zz;
511 rp = sqrt (xp*xp + yp*yp + zp*zp);
512 rp = pow (pow (rp, 4) + TP4, 0.25);
513 if (rp < TP_Tiny)
514 rp = TP_Tiny;
515 ir = 1.0/rp;
516
517 if (rp < TP_Extend_Radius) {
518 ir = EXTEND(1., rp);
519 }
520
521 s1 = 0.5* mp *ir;
522 s3 = -s1*ir*ir;
523 s5 = -3.0*s3*ir*ir;
524
525 p += s1;
526
527 px += xp*s3;
528 py += yp*s3;
529 pz += zp*s3;
530
531 pxx += xp*xp*s5 + s3;
532 pxy += xp*yp*s5;
533 pxz += xp*zp*s5;
534 pyy += yp*yp*s5 + s3;
535 pyz += yp*zp*s5;
536 pzz += zp*zp*s5 + s3;
537
538 /* second puncture */
539 xp = xx + par_b;
540 yp = yy;
541 zp = zz;
542 rp = sqrt (xp*xp + yp*yp + zp*zp);
543 rp = pow (pow (rp, 4) + TP4, 0.25);
544 if (rp < TP_Tiny)
545 rp = TP_Tiny;
546 ir = 1.0/rp;
547
548 if (rp < TP_Extend_Radius) {
549 ir = EXTEND(1., rp);
550 }
551
552 s1 = 0.5* mm *ir;
553 s3 = -s1*ir*ir;
554 s5 = -3.0*s3*ir*ir;
555
556 p += s1;
557
558 px += xp*s3;
559 py += yp*s3;
560 pz += zp*s3;
561
562 pxx += xp*xp*s5 + s3;
563 pxy += xp*yp*s5;
564 pxz += xp*zp*s5;
565 pyy += yp*yp*s5 + s3;
566 pyz += yp*zp*s5;
567 pzz += zp*zp*s5 + s3;
568
569 // @TODO: psix, psiy, psiz, psixx, psixy are not defined, currently.
570 // They were stored by Thorn einsteinbase/StaticConformal
571
572 if (conformal_state >= 1) {
573 static_psi = p;
575 }
576 if (conformal_state >= 2) {
580 }
581 if (conformal_state >= 3) {
588 }
589
590 if (pmn_lapse)
593 Q[lapse] = 2.0/(1.0+pow(p, initial_lapse_psi_exponent));
594
595 } /* if conformal-state > 0 */
596
598
599 // First, zero everything
600 // @TODO: it should be checked if this represents Vaccuum
601 for(int i=0; i<Qlen; i++) Q[i] = 0;
602
603 Q[g11] = pow (psi1 / static_psi, 4);
604 Q[g12] = 0;
605 Q[g13] = 0;
606 Q[g22] = pow (psi1 / static_psi, 4);
607 Q[g23]= 0;
608 Q[g33] = pow (psi1 / static_psi, 4);
609
610 const double oneoverpsisq = 1./(psi1*psi1);
611 Q[K11] = Aij[0][0] * oneoverpsisq;// / pow(psi1, 2);
612 Q[K12] = Aij[0][1] * oneoverpsisq;// / pow(psi1, 2);
613 Q[K13] = Aij[0][2] * oneoverpsisq;// / pow(psi1, 2);
614 Q[K22] = Aij[1][1] * oneoverpsisq;// / pow(psi1, 2);
615 Q[K23] = Aij[1][2] * oneoverpsisq;// / pow(psi1, 2);
616 Q[K33] = Aij[2][2] * oneoverpsisq;// / pow(psi1, 2);
617 //here provide the "physical" metric and ex curvature (i.e. without any bar), the real decomposition of CCZ4 is done in ADM class.
619 Q[lapse] =
620 ((1.0 -0.5* mp /r_plus -0.5* mm/r_minus)
621 /(1.0 +0.5* mp /r_plus +0.5* mm/r_minus));
622
623 if (r_plus < TP_Extend_Radius) {
624 Q[lapse] =
625 ((1.0 -0.5*EXTEND(mp, r_plus) -0.5* mm/r_minus)
626 /(1.0 +0.5*EXTEND(mp, r_plus) +0.5* mm/r_minus));
627 }
628 if (r_minus < TP_Extend_Radius) {
629 Q[lapse] =
630 ((1.0 -0.5*EXTEND(mm, r_minus) -0.5* mp/r_plus)
631 /(1.0 +0.5*EXTEND(mp, r_minus) +0.5* mp/r_plus));
632 }
633
634 if (averaged_lapse) {
635 Q[lapse] = 0.5 * (1.0 + Q[lapse]);
636 }
637 }
639 Q[lapse] *= old_alp;
640
641 if (swap_xz) {
642 /* Swap the x and z components of all tensors */
643 if (conformal_state >= 2) {
645 }
646 if (conformal_state >= 3) {
649 }
650 SWAP (Q[g11], Q[g33]);
651 SWAP (Q[g12], Q[g23]);
652 SWAP (Q[K11], Q[K33]);
653 SWAP (Q[K12], Q[K23]);
654 } /* if swap_xz */
655
657 {
658 TP_WARN("@TODO Rescale_Sources is in some other thorn and has not been copied");
659#if 0 == 1
660 Rescale_Sources(cctkGH,
661 cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2],
662 x, y, z,
663 (conformal_state > 0) ? psi : NULL,
664 gxx, gyy, gzz,
665 gxy, gxz, gyz);
666#endif
667 }
668
669#if 0 == 1
670 /* Cleanup not possible as we don't know when the grid queries are finished */
671 /* Keep the result around for the next time */
672 free_dvector (F, 0, ntotal - 1);
673 free_derivs (&u, ntotal);
674 free_derivs (&v, ntotal);
675 free_derivs (&cf_v, ntotal);
676#endif
677}
678
679} // namespace TP
680
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 Pi
#define SWAP(a, b)
#define Pih
#define EXTEND(M, r)
void free_derivs(derivs *v, int n)
GRID_SETUP_METHOD gsm
void Derivatives_AB3(int nvar, int n1, int n2, int n3, derivs v)
bool runned
an internal check
void SpecCoef(int n1, int n2, int n3, int ivar, double *v, double *cf)
void Interpolate(const double *const pos, double *Q, bool low_res=false)
Interpolation function for an external caller.
void F_of_v(int nvar, int n1, int n2, int n3, derivs v, double *F, derivs u)
void Newton(int nvar, int n1, int n2, int n3, derivs v, double tol, int itmax)
Definition Newton.cpp:452
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 set_initial_guess(derivs v)
void BY_Aijofxyz(double x, double y, double z, double Aij[3][3])
Definition Equations.cpp:83
double PunctTaylorExpandAtArbitPosition(int ivar, int nvar, int n1, int n2, int n3, derivs v, double x, double y, double z)
double PunctIntPolAtArbitPosition(int ivar, int nvar, int n1, int n2, int n3, derivs v, double x, double y, double z)
void allocate_derivs(derivs *v, int n)
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)
double * dvector(long nl, long nh)
void free_dvector(double *v, long nl, long nh)
This file contains aliases for making access to the long state vector Q as used eg.
@ GSM_Taylor_expansion
@ GSM_evaluation
int TP_MyProc()
double * d12
double * d23
double * d1
double * d11
double * d3
double * d33
double * d13
double * d0
double * d22
double * d2
void TP_INFO(const char *fmt,...)
void TP_ERROR(const char *fmt,...)
void TP_WARN(const char *fmt,...)
double par_S_plus[3]
std::string grid_setup_method
std::string metric_type
bool use_external_initial_guess
static constexpr int npoints_phi
double par_P_plus[3]
bool solve_momentum_constraint
std::string conformal_storage
double par_S_minus[3]
static constexpr int npoints_A
double TP_Extend_Radius
double center_offset[3]
double target_M_minus
static constexpr int npoints_B_low
double par_P_minus[3]
static constexpr int npoints_phi_low
static constexpr int npoints_B
static constexpr int npoints_A_low
double initial_lapse_psi_exponent
std::string initial_lapse
bool do_initial_debug_output