LORENE
blackhole_eq_spher.C
1/*
2 * Method of class Black_hole to compute a single black hole
3 *
4 * (see file blackhole.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2005-2007 Keisuke Taniguchi
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License version 2
15 * as published by the Free Software Foundation.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28
29
30/*
31 * $Id: blackhole_eq_spher.C,v 1.6 2016/12/05 16:17:48 j_novak Exp $
32 * $Log: blackhole_eq_spher.C,v $
33 * Revision 1.6 2016/12/05 16:17:48 j_novak
34 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35 *
36 * Revision 1.5 2014/10/13 08:52:46 j_novak
37 * Lorene classes and functions now belong to the namespace Lorene.
38 *
39 * Revision 1.4 2014/10/06 15:13:02 j_novak
40 * Modified #include directives to use c++ syntax.
41 *
42 * Revision 1.3 2008/07/02 20:42:53 k_taniguchi
43 * Modification of the argument and so on.
44 *
45 * Revision 1.2 2008/05/15 19:26:30 k_taniguchi
46 * Change of some parameters.
47 *
48 * Revision 1.1 2007/06/22 01:19:11 k_taniguchi
49 * *** empty log message ***
50 *
51 *
52 * $Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_eq_spher.C,v 1.6 2016/12/05 16:17:48 j_novak Exp $
53 *
54 */
55
56// C++ headers
57//#include <>
58
59// C headers
60#include <cmath>
61
62// Lorene headers
63#include "blackhole.h"
64#include "cmp.h"
65#include "tenseur.h"
66#include "param.h"
67#include "unites.h"
68#include "proto.h"
69#include "utilitaires.h"
70#include "graphique.h"
71
72namespace Lorene {
73void Black_hole::equilibrium_spher(bool neumann, bool first,
74 double spin_omega, double precis,
75 double precis_shift) {
76
77 // Fundamental constants and units
78 // -------------------------------
79 using namespace Unites ;
80
81 // Initializations
82 // ---------------
83
84 const Mg3d* mg = mp.get_mg() ;
85 int nz = mg->get_nzone() ; // total number of domains
86
87 double mass = ggrav * mass_bh ;
88
89 // Inner boundary condition
90 // ------------------------
91
92 Valeur bc_lpcnf(mg->get_angu()) ;
93 Valeur bc_conf(mg->get_angu()) ;
94
95 Valeur bc_shif_x(mg->get_angu()) ;
96 Valeur bc_shif_y(mg->get_angu()) ;
97 Valeur bc_shif_z(mg->get_angu()) ;
98
99 Scalar rr(mp) ;
100 rr = mp.r ;
101 rr.std_spectral_base() ;
102 Scalar st(mp) ;
103 st = mp.sint ;
104 st.std_spectral_base() ;
105 Scalar ct(mp) ;
106 ct = mp.cost ;
107 ct.std_spectral_base() ;
108 Scalar sp(mp) ;
109 sp = mp.sinp ;
110 sp.std_spectral_base() ;
111 Scalar cp(mp) ;
112 cp = mp.cosp ;
113 cp.std_spectral_base() ;
114
115 Vector ll(mp, CON, mp.get_bvect_cart()) ;
116 ll.set_etat_qcq() ;
117 ll.set(1) = st * cp ;
118 ll.set(2) = st * sp ;
119 ll.set(3) = ct ;
120 ll.std_spectral_base() ;
121
122 // Sets C/M^2 for each case of the lapse boundary condition
123 // --------------------------------------------------------
124 double cc ;
125
126 if (neumann) { // Neumann boundary condition
127 if (first) { // First condition
128 // d(\alpha \psi)/dr = 0
129 // ---------------------
130 cc = 2. * (sqrt(13.) - 1.) / 3. ;
131 }
132 else { // Second condition
133 // d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
134 // -----------------------------------------
135 cc = 4. / 3. ;
136 }
137 }
138 else { // Dirichlet boundary condition
139 if (first) { // First condition
140 // (\alpha \psi) = 1/2
141 // -------------------
142 cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
143 abort() ;
144 }
145 else { // Second condition
146 // (\alpha \psi) = 1/sqrt(2.) \psi_KS
147 // ----------------------------------
148 cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
149 abort() ;
150 // cc = 2. * sqrt(2.) ;
151 }
152 }
153
154
155 // Ininital metric
156 if (kerrschild) {
157
158 lapconf_bh = 1. / sqrt(1. + 2. * mass / rr) ;
159 lapconf_bh.std_spectral_base() ;
160 lapconf_bh.annule_domain(0) ;
161 lapconf_bh.raccord(1) ;
162
163 lapconf = lapconf_bh ;
164 lapconf.std_spectral_base() ;
165 lapconf_rs = 0. ;
166 lapconf_rs.std_spectral_base() ;
167
168 lapse = lapconf ;
169 lapse.std_spectral_base() ;
170
171 confo = 1. ;
172 confo.std_spectral_base() ;
173
174 Scalar lapse_bh(mp) ;
175 lapse_bh = 1. / sqrt(1. + 2. * mass / rr) ;
176 lapse_bh.std_spectral_base() ;
177 lapse_bh.annule_domain(0) ;
178 lapse_bh.raccord(1) ;
179
180 for (int i=1; i<=3; i++) {
181 shift_bh.set(i) = 2. * lapse_bh * lapse_bh * mass * ll(i) / rr ;
182 }
183 shift_bh.std_spectral_base() ;
184
185 shift = shift_bh ;
186 shift.std_spectral_base() ;
187 shift_rs.set_etat_zero() ;
188
189 }
190 else { // Isotropic coordinates
191
192 Scalar r_are(mp) ;
193 r_are = r_coord(neumann, first) ;
194 r_are.std_spectral_base() ;
195 r_are.annule_domain(0) ;
196 r_are.raccord(1) ;
197 /*
198 cout << "r_are:" << endl ;
199 for (int l=0; l<nz; l++) {
200 cout << r_are.val_grid_point(l,0,0,0) << endl ;
201 }
202 */
203
204 // Exact, non-spinning case
205 /*
206 lapconf = sqrt(1. - 2.*mass/r_are/rr
207 + cc*cc*pow(mass/r_are/rr,4.))
208 * sqrt(r_are) ;
209 */
210 lapconf = sqrt(1. - 1.9*mass/r_are/rr
211 + cc*cc*pow(mass/r_are/rr,4.))
212 * sqrt(r_are) ;
213 lapconf.std_spectral_base() ;
214 lapconf.annule_domain(0) ;
215 lapconf.raccord(1) ;
216
217 /*
218 lapse = sqrt(1. - 2.*mass/r_are/rr
219 + cc*cc*pow(mass/r_are/rr,4.)) ;
220 */
221 lapse = sqrt(1. - 1.9*mass/r_are/rr
222 + cc*cc*pow(mass/r_are/rr,4.)) ;
223 lapse.std_spectral_base() ;
224 lapse.annule_domain(0) ;
225 lapse.raccord(1) ;
226
227 // confo = sqrt(r_are) ;
228 confo = sqrt(0.9*r_are) ;
229 confo.std_spectral_base() ;
230 confo.annule_domain(0) ;
231 confo.raccord(1) ;
232
233 for (int i=1; i<=3; i++) {
234 shift.set(i) = mass * mass * cc * ll(i) / rr / rr
235 / pow(r_are,3.) ;
236 }
237
238 shift.std_spectral_base() ;
239
240 for (int i=1; i<=3; i++) {
241 shift.set(i).annule_domain(0) ;
242 shift.set(i).raccord(1) ;
243 }
244
245 /*
246 des_profile( r_are, 0., 20, M_PI/2., 0.,
247 "Areal coordinate",
248 "Areal (theta=pi/2, phi=0)" ) ;
249
250 des_profile( lapse, 0., 20, M_PI/2., 0.,
251 "Initial lapse function of BH",
252 "Lapse (theta=pi/2, phi=0)" ) ;
253
254 des_profile( confo, 0., 20, M_PI/2., 0.,
255 "Initial conformal factor of BH",
256 "Confo (theta=pi/2, phi=0)" ) ;
257
258 des_profile( shift(1), 0., 20, M_PI/2., 0.,
259 "Initial shift vector (X) of BH",
260 "Shift (theta=pi/2, phi=0)" ) ;
261
262 des_coupe_vect_z( shift, 0., -3., 0.5, 4,
263 "Shift vector of BH") ;
264 */
265 }
266
267 // Auxiliary quantities
268 // --------------------
269
270 Scalar source_lapconf(mp) ;
271 Scalar source_confo(mp) ;
272 Vector source_shift(mp, CON, mp.get_bvect_cart()) ;
273
274 Scalar lapconf_m1(mp) ; // = lapconf - 1 (only for the isotropic case)
275 Scalar confo_m1(mp) ; // = confo - 1
276
277 Scalar lapconf_jm1(mp) ;
278 Scalar confo_jm1(mp) ;
279 Vector shift_jm1(mp, CON, mp.get_bvect_cart()) ;
280
281 double diff_lp = 1. ;
282 double diff_cf = 1. ;
283 double diff_sx = 1. ;
284 double diff_sy = 1. ;
285 double diff_sz = 1. ;
286
287 int mermax = 200 ; // max number of iterations
288
289 //======================================//
290 // Start of iteration //
291 //======================================//
292 /*
293 for (int mer=0;
294 (diff_lp > precis) || (diff_cf > precis) && (mer < mermax); mer++) {
295
296 for (int mer=0;
297 (diff_sx > precis) || (diff_sy > precis) || (diff_sz > precis)
298 && (mer < mermax); mer++) {
299 */
300 for (int mer=0; (diff_lp > precis) && (diff_cf > precis) && (diff_sx > precis) && (diff_sy > precis) && (diff_sz > precis) && (mer < mermax); mer++) {
301
302 cout << "--------------------------------------------------" << endl ;
303 cout << "step: " << mer << endl ;
304 cout << "diff_lapconf = " << diff_lp << endl ;
305 cout << "diff_confo = " << diff_cf << endl ;
306 cout << "diff_shift : x = " << diff_sx
307 << " y = " << diff_sy << " z = " << diff_sz << endl ;
308
309 if (kerrschild) {
310 lapconf_jm1 = lapconf_rs ;
311 confo_jm1 = confo ;
312 shift_jm1 = shift_rs ;
313 }
314 else {
315 lapconf_jm1 = lapconf ;
316 confo_jm1 = confo ;
317 shift_jm1 = shift ;
318 }
319
320 //------------------------------------------//
321 // Computation of the extrinsic curvature //
322 //------------------------------------------//
323
324 extr_curv_bh() ;
325
326 //---------------------------------------------------------------//
327 // Resolution of the Poisson equation for the lapconf function //
328 //---------------------------------------------------------------//
329
330 // Source term
331 // -----------
332
333 if (kerrschild) {
334
335 cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
336 abort() ;
337
338 }
339 else { // Isotropic coordinates with the maximal slicing
340
341 source_lapconf = 7. * lapconf_jm1 % taij_quad
342 / pow(confo_jm1, 8.) / 8. ;
343
344 }
345
346 source_lapconf.annule_domain(0) ;
347 source_lapconf.set_dzpuis(4) ;
348 source_lapconf.std_spectral_base() ;
349
350 /*
351 Scalar tmp_source = source_lapse ;
352 tmp_source.dec_dzpuis(4) ;
353 tmp_source.std_spectral_base() ;
354
355 des_profile( tmp_source, 0., 20, M_PI/2., 0.,
356 "Source term of lapse",
357 "source_lapse (theta=pi/2, phi=0)" ) ;
358 */
359
360 bc_lpcnf = bc_lapconf(neumann, first) ;
361
362
363 if (kerrschild) {
364
365 cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
366 abort() ;
367 /*
368 lapconf_rs.set_etat_qcq() ;
369
370 if (neumann) {
371 lapconf_rs = source_lapconf.poisson_neumann(bc_lpcnf, 0) ;
372 }
373 else {
374 lapconf_rs = source_lapconf.poisson_dirichlet(bc_lpcnf, 0) ;
375 }
376
377 // Re-construction of the lapse function
378 // -------------------------------------
379 lapconf_rs.annule_domain(0) ;
380 lapconf_rs.raccord(1) ;
381
382 lapconf = lapconf_rs + lapconf_bh ;
383 lapconf.annule_domain(0) ;
384 lapconf.raccord(1) ;
385 */
386 }
387 else { // Isotropic coordinates with the maximal slicing
388
389 lapconf_m1.set_etat_qcq() ;
390
391 if (neumann) {
392 lapconf_m1 = source_lapconf.poisson_neumann(bc_lpcnf, 0) ;
393 }
394 else {
395 lapconf_m1 = source_lapconf.poisson_dirichlet(bc_lpcnf, 0) ;
396 }
397
398 // Re-construction of the lapse function
399 // -------------------------------------
400 lapconf = lapconf_m1 + 1. ;
401 lapconf.annule_domain(0) ;
402 lapconf.raccord(1) ;
403 /*
404 des_profile( lapse, 0., 20, M_PI/2., 0.,
405 "Lapse function of BH",
406 "Lapse (theta=pi/2, phi=0)" ) ;
407 */
408 }
409
410 //---------------------------------------------------------------//
411 // Resolution of the Poisson equation for the conformal factor //
412 //---------------------------------------------------------------//
413
414 // Source term
415 // -----------
416
417 if (kerrschild) {
418
419 cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
420 abort() ;
421
422 }
423 else { // Isotropic coordinates with the maximal slicing
424
425 Scalar tmp1 = - 0.125 * taij_quad / pow(confo_jm1, 7.) ;
426 tmp1.std_spectral_base() ;
427 tmp1.inc_dzpuis(4-tmp1.get_dzpuis()) ;
428
429 source_confo = tmp1 ;
430
431 }
432
433 source_confo.annule_domain(0) ;
434 source_confo.set_dzpuis(4) ;
435 source_confo.std_spectral_base() ;
436
437 bc_conf = bc_confo() ;
438
439 confo_m1.set_etat_qcq() ;
440
441 confo_m1 = source_confo.poisson_neumann(bc_conf, 0) ;
442
443 // Re-construction of the conformal factor
444 // ---------------------------------------
445
446 confo = confo_m1 + 1. ;
447 confo.annule_domain(0) ;
448 confo.raccord(1) ;
449
450 //-----------------------------------------------------------//
451 // Resolution of the Poisson equation for the shift vector //
452 //-----------------------------------------------------------//
453
454 // Source term
455 // -----------
456
457 Scalar confo7(mp) ;
458 confo7 = pow(confo_jm1, 7.) ;
459 confo7.std_spectral_base() ;
460
461 Vector dlappsi(mp, COV, mp.get_bvect_cart()) ;
462 for (int i=1; i<=3; i++)
463 dlappsi.set(i) = (lapconf_jm1.deriv(i)
464 - 7.*lapconf*confo_jm1.deriv(i)/confo_jm1)
465 / confo7 ;
466
467 dlappsi.std_spectral_base() ;
468 dlappsi.annule_domain(0) ;
469
470
471 if (kerrschild) {
472
473 cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
474 abort() ;
475
476 }
477 else { // Isotropic coordinates with the maximal slicing
478
479 source_shift = 2. * contract(taij, 1, dlappsi, 0) ;
480
481 }
482
483 source_shift.annule_domain(0) ;
484 /*
485 for (int i=1; i<=3; i++)
486 (source_shift.set(i)).raccord(1) ;
487 */
488 /*
489 for (int i=1; i<=3; i++) {
490 (source_shift.set(i)).set_dzpuis(4) ;
491 }
492 */
493 source_shift.std_spectral_base() ;
494
495 for (int i=1; i<=3; i++) {
496 if (source_shift(i).get_etat() != ETATZERO)
497 source_shift.set(i).filtre(4) ;
498 }
499
500 /*
501 for (int i=1; i<=3; i++) {
502 if (source_shift(i).dz_nonzero()) {
503 assert( source_shift(i).get_dzpuis() == 4 ) ;
504 }
505 else {
506 (source_shift.set(i)).set_dzpuis(4) ;
507 }
508 }
509 */
510
511 Tenseur source_p(mp, 1, CON, mp.get_bvect_cart()) ;
512 source_p.set_etat_qcq() ;
513 for (int i=0; i<3; i++) {
514 source_p.set(i) = Cmp(source_shift(i+1)) ;
515 }
516
517 Tenseur resu_p(mp, 1, CON, mp.get_bvect_cart()) ;
518 resu_p.set_etat_qcq() ;
519
520 for (int i=0; i<3; i++) {
521 resu_p.set(i) = Cmp(shift_jm1(i+1)) ;
522 }
523
524 bc_shif_x = bc_shift_x(spin_omega) ; // Rotating BH
525 bc_shif_y = bc_shift_y(spin_omega) ; // Rotating BH
526 bc_shif_z = bc_shift_z() ;
527 /*
528 cout << bc_shif_x << endl ;
529 arrete() ;
530 cout << bc_shif_y << endl ;
531 arrete() ;
532 cout << bc_shif_z << endl ;
533 arrete() ;
534 */
535 poisson_vect_frontiere(1./3., source_p, resu_p,
536 bc_shif_x, bc_shif_y, bc_shif_z,
537 0, precis_shift, 14) ;
538
539
540 if (kerrschild) {
541 for (int i=1; i<=3; i++)
542 shift_rs.set(i) = resu_p(i-1) ;
543
544 for (int i=1; i<=3; i++)
545 shift.set(i) = shift_rs(i) + shift_bh(i) ;
546
547 shift_rs.annule_domain(0) ;
548 }
549 else { // Isotropic coordinates with the maximal slicing
550 for (int i=1; i<=3; i++)
551 shift.set(i) = resu_p(i-1) ;
552 }
553
554 shift.annule_domain(0) ;
555
556 for (int i=1; i<=3; i++)
557 shift.set(i).raccord(1) ;
558
559
560 /*
561 Tbl diff_shftx = diffrel(shift(1), shift_ex(1)) ;
562 double diff_shfx = diff_shftx(1) ;
563 for (int l=2; l<nz; l++) {
564 diff_shfx += diff_shftx(l) ;
565 }
566 diff_shfx /= nz ;
567
568 cout << "diff_shfx : " << diff_shfx << endl ;
569 */
570
571 //------------------------------------------------//
572 // Relative difference in the metric quantities //
573 //------------------------------------------------//
574
575 /*
576 des_profile( lapse, 0., 20, M_PI/2., 0.,
577 "Lapse function of BH",
578 "Lapse (theta=pi/2, phi=0)" ) ;
579
580 des_profile( confo, 0., 20, M_PI/2., 0.,
581 "Conformal factor of BH",
582 "Confo (theta=pi/2, phi=0)" ) ;
583
584 des_coupe_vect_z( shift, 0., -3., 0.5, 4,
585 "Shift vector of BH") ;
586 */
587 // Difference is calculated only outside the inner boundary.
588
589 // Lapconf function
590 // ----------------
591 Tbl diff_lapconf(nz) ;
592
593 if (kerrschild) {
594
595 diff_lapconf = diffrel(lapconf_rs, lapconf_jm1) ;
596
597 }
598 else { // Isotropic coordinates with the maximal slicing
599
600 diff_lapconf = diffrel(lapconf, lapconf_jm1) ;
601
602 }
603 /*
604 cout << "lapse: " << endl ;
605 for (int l=0; l<nz; l++) {
606 cout << lapse.val_grid_point(l,0,0,0) << endl ;
607 }
608
609 cout << "lapse_jm1: " << endl ;
610 for (int l=0; l<nz; l++) {
611 cout << lapse_jm1.val_grid_point(l,0,0,0) << endl ;
612 }
613 */
614
615 cout << "Relative difference in the lapconf function : " << endl ;
616 for (int l=0; l<nz; l++) {
617 cout << diff_lapconf(l) << " " ;
618 }
619 cout << endl ;
620
621 diff_lp = diff_lapconf(1) ;
622 for (int l=2; l<nz; l++) {
623 diff_lp += diff_lapconf(l) ;
624 }
625 diff_lp /= nz ;
626
627 // Conformal factor
628 // ----------------
629 Tbl diff_confo = diffrel(confo, confo_jm1) ;
630
631 cout << "Relative difference in the conformal factor : " << endl ;
632 for (int l=0; l<nz; l++) {
633 cout << diff_confo(l) << " " ;
634 }
635 cout << endl ;
636
637 diff_cf = diff_confo(1) ;
638 for (int l=2; l<nz; l++) {
639 diff_cf += diff_confo(l) ;
640 }
641 diff_cf /= nz ;
642
643 // Shift vector
644 // ------------
645 Tbl diff_shift_x(nz) ;
646 Tbl diff_shift_y(nz) ;
647 Tbl diff_shift_z(nz) ;
648
649 if (kerrschild) {
650
651 diff_shift_x = diffrel(shift_rs(1), shift_jm1(1)) ;
652
653 }
654 else { // Isotropic coordinates with the maximal slicing
655
656 diff_shift_x = diffrel(shift(1), shift_jm1(1)) ;
657
658 }
659
660 cout << "Relative difference in the shift vector (x) : " << endl ;
661 for (int l=0; l<nz; l++) {
662 cout << diff_shift_x(l) << " " ;
663 }
664 cout << endl ;
665
666 diff_sx = diff_shift_x(1) ;
667 for (int l=2; l<nz; l++) {
668 diff_sx += diff_shift_x(l) ;
669 }
670 diff_sx /= nz ;
671
672 if (kerrschild) {
673
674 diff_shift_y = diffrel(shift_rs(2), shift_jm1(2)) ;
675
676 }
677 else { // Isotropic coordinates with the maximal slicing
678
679 diff_shift_y = diffrel(shift(2), shift_jm1(2)) ;
680
681 }
682
683 cout << "Relative difference in the shift vector (y) : " << endl ;
684 for (int l=0; l<nz; l++) {
685 cout << diff_shift_y(l) << " " ;
686 }
687 cout << endl ;
688
689 diff_sy = diff_shift_y(1) ;
690 for (int l=2; l<nz; l++) {
691 diff_sy += diff_shift_y(l) ;
692 }
693 diff_sy /= nz ;
694
695 if (kerrschild) {
696
697 diff_shift_z = diffrel(shift_rs(3), shift_jm1(3)) ;
698
699 }
700 else { // Isotropic coordinates with the maximal slicing
701
702 diff_shift_z = diffrel(shift(3), shift_jm1(3)) ;
703
704 }
705
706 cout << "Relative difference in the shift vector (z) : " << endl ;
707 for (int l=0; l<nz; l++) {
708 cout << diff_shift_z(l) << " " ;
709 }
710 cout << endl ;
711
712 diff_sz = diff_shift_z(1) ;
713 for (int l=2; l<nz; l++) {
714 diff_sz += diff_shift_z(l) ;
715 }
716 diff_sz /= nz ;
717
718 // Mass
719 if (kerrschild) {
720
721 cout << "Mass_bh : " << mass_bh / msol << " [M_sol]" << endl ;
722 double rad_apphor = rad_ah() ;
723 cout << " : " << 0.5 * rad_apphor / ggrav / msol
724 << " [M_sol]" << endl ;
725
726 }
727 else { // Isotropic coordinates with the maximal slicing
728
729 cout << "Mass_bh : " << mass_bh / msol
730 << " [M_sol]" << endl ;
731
732 }
733
734 // ADM mass, Komar mass
735 // --------------------
736
737 double irr_gm, adm_gm, kom_gm ;
738 irr_gm = mass_irr() / mass_bh - 1. ;
739 adm_gm = mass_adm() / mass_bh - 1. ;
740 kom_gm = mass_kom() / mass_bh - 1. ;
741 cout << "Irreducible mass : " << mass_irr() / msol << endl ;
742 cout << "Gravitaitonal mass : " << mass_bh / msol << endl ;
743 cout << "ADM mass : " << mass_adm() / msol << endl ;
744 cout << "Komar mass : " << mass_kom() / msol << endl ;
745 cout << "Diff. (Madm-Mirr)/Mirr : " << mass_adm()/mass_irr() - 1.
746 << endl ;
747 cout << "Diff. (Mkom-Mirr)/Mirr : " << mass_kom()/mass_irr() - 1.
748 << endl ;
749 cout << "Diff. (Madm-Mkom)/Madm : " << 1. - mass_kom()/mass_adm()
750 << endl ;
751 cout << "Diff. (Mirr-Mg)/Mg : " << irr_gm << endl ;
752 cout << "Diff. (Madm-Mg)/Mg : " << adm_gm << endl ;
753 cout << "Diff. (Mkom-Mg)/Mg : " << kom_gm << endl ;
754
755 cout << endl ;
756
757 del_deriv() ;
758
759 // Relaxation
760 /*
761 lapse = 0.75 * lapse + 0.25 * lapse_jm1 ;
762 confo = 0.75 * confo + 0.25 * confo_jm1 ;
763 shift = 0.75 * shift + 0.25 * shift_jm1 ;
764 */
765 /*
766 des_profile( lapse, 0., 20, M_PI/2., 0.,
767 "Lapse function of BH",
768 "Lapse (theta=pi/2, phi=0)" ) ;
769
770 des_profile( confo, 0., 20, M_PI/2., 0.,
771 "Conformal factor of BH",
772 "Confo (theta=pi/2, phi=0)" ) ;
773
774 des_profile( shift(1), 0., 20, M_PI/2., 0.,
775 "Shift vector (X) of BH",
776 "Shift (theta=pi/2, phi=0)" ) ;
777 */
778
779 } // End of iteration loop
780
781 //====================================//
782 // End of iteration //
783 //====================================//
784
785 // Exact solution
786 // --------------
787 Scalar lapconf_exact(mp) ;
788 Scalar confo_exact(mp) ;
789 Vector shift_exact(mp, CON, mp.get_bvect_cart()) ;
790
791 if (kerrschild) {
792
793 lapconf_exact = 1./sqrt(1.+2.*mass/rr) ;
794
795 confo_exact = 1. ;
796
797 for (int i=1; i<=3; i++)
798 shift_exact.set(i) =
799 2.*mass*lapconf_exact%lapconf_exact%ll(i)/rr ;
800
801 }
802 else {
803
804 Scalar rare(mp) ;
805 rare = r_coord(neumann, first) ;
806 rare.std_spectral_base() ;
807
808 lapconf_exact = sqrt(1. - 2.*mass/rare/rr
809 + cc*cc*pow(mass/rare/rr,4.)) * sqrt(rare) ;
810
811 confo_exact = sqrt(rare) ;
812
813 for (int i=1; i<=3; i++) {
814 shift_exact.set(i) = mass * mass * cc * ll(i) / rr / rr
815 / pow(rare,3.) ;
816 }
817
818 }
819
820 lapconf_exact.annule_domain(0) ;
821 lapconf_exact.std_spectral_base() ;
822 lapconf_exact.raccord(1) ;
823
824 confo_exact.annule_domain(0) ;
825 confo_exact.std_spectral_base() ;
826 confo_exact.raccord(1) ;
827
828 shift_exact.annule_domain(0) ;
829 shift_exact.std_spectral_base() ;
830 for (int i=1; i<=3; i++)
831 shift_exact.set(i).raccord(1) ;
832
833 Scalar lapconf_resi = lapconf - lapconf_exact ;
834 Scalar confo_resi = confo - confo_exact ;
835 Vector shift_resi = shift - shift_exact ;
836 /*
837 des_profile( lapse, 0., 20, M_PI/2., 0.,
838 "Lapse function",
839 "Lapse (theta=pi/2, phi=0)" ) ;
840
841 des_profile( lapse_exact, 0., 20, M_PI/2., 0.,
842 "Exact lapse function",
843 "Exact lapse (theta=pi/2, phi=0)" ) ;
844
845 des_profile( lapse_resi, 0., 20, M_PI/2., 0.,
846 "Residual of the lapse function",
847 "Delta Lapse (theta=pi/2, phi=0)" ) ;
848
849 des_profile( confo, 0., 20, M_PI/2., 0.,
850 "Conformal factor",
851 "Confo (theta=pi/2, phi=0)" ) ;
852
853 des_profile( confo_exact, 0., 20, M_PI/2., 0.,
854 "Exact conformal factor",
855 "Exact confo (theta=pi/2, phi=0)" ) ;
856
857 des_profile( confo_resi, 0., 20, M_PI/2., 0.,
858 "Residual of the conformal factor",
859 "Delta Confo (theta=pi/2, phi=0)" ) ;
860
861 des_profile( shift(1), 0., 20, M_PI/2., 0.,
862 "Shift vector (X)",
863 "Shift (X) (theta=pi/2, phi=0)" ) ;
864
865 des_profile( shift_exact(1), 0., 20, M_PI/2., 0.,
866 "Exact shift vector (X)",
867 "Exact shift (X) (theta=pi/2, phi=0)" ) ;
868
869 des_profile( shift_resi(1), 0., 20, M_PI/2., 0.,
870 "Residual of the shift vector X",
871 "Delta shift (X) (theta=pi/2, phi=0)" ) ;
872 */
873 /*
874 des_coupe_vect_z( shift_resi, 0., -3., 0.5, 4,
875 "Delta Shift vector of BH") ;
876 */
877
878 // Relative difference in the lapconf function
879 Tbl diff_lapconf_exact = diffrel(lapconf, lapconf_exact) ;
880 diff_lapconf_exact.set(0) = 0. ;
881 cout << "Relative difference in the lapconf function : " << endl ;
882 for (int l=0; l<nz; l++) {
883 cout << diff_lapconf_exact(l) << " " ;
884 }
885 cout << endl ;
886
887 // Relative difference in the conformal factor
888 Tbl diff_confo_exact = diffrel(confo, confo_exact) ;
889 diff_confo_exact.set(0) = 0. ;
890 cout << "Relative difference in the conformal factor : " << endl ;
891 for (int l=0; l<nz; l++) {
892 cout << diff_confo_exact(l) << " " ;
893 }
894 cout << endl ;
895
896 // Relative difference in the shift vector
897 Tbl diff_shift_exact_x = diffrel(shift(1), shift_exact(1)) ;
898 Tbl diff_shift_exact_y = diffrel(shift(2), shift_exact(2)) ;
899 Tbl diff_shift_exact_z = diffrel(shift(3), shift_exact(3)) ;
900
901 cout << "Relative difference in the shift vector (x) : " << endl ;
902 for (int l=0; l<nz; l++) {
903 cout << diff_shift_exact_x(l) << " " ;
904 }
905 cout << endl ;
906 cout << "Relative difference in the shift vector (y) : " << endl ;
907 for (int l=0; l<nz; l++) {
908 cout << diff_shift_exact_y(l) << " " ;
909 }
910 cout << endl ;
911 cout << "Relative difference in the shift vector (z) : " << endl ;
912 for (int l=0; l<nz; l++) {
913 cout << diff_shift_exact_z(l) << " " ;
914 }
915 cout << endl ;
916
917 //---------------------------------//
918 // Info printing //
919 //---------------------------------//
920
921
922}
923}
Lorene::Black_hole::shift_bh
Vector shift_bh
Part of the shift vector from the analytic background.
Definition blackhole.h:115
Lorene::Black_hole::mass_kom
virtual double mass_kom() const
Komar mass.
Definition blackhole_global.C:292
Lorene::Black_hole::shift_rs
Vector shift_rs
Part of the shift vector from the numerical computation.
Definition blackhole.h:112
Lorene::Black_hole::bc_confo
const Valeur bc_confo() const
Boundary condition on the apparent horizon of the black hole for the conformal factor: 2-D Valeur.
Definition blackhole_bc.C:413
Lorene::Black_hole::taij_quad
Scalar taij_quad
Part of the scalar generated by .
Definition blackhole.h:135
Lorene::Black_hole::bc_shift_x
const Valeur bc_shift_x(double omega_r) const
Boundary condition on the apparent horizon of the black hole for the shift vector of the \fx\f direct...
Definition blackhole_bc.C:214
Lorene::Black_hole::taij
Sym_tensor taij
Trace of the extrinsic curvature.
Definition blackhole.h:127
Lorene::Black_hole::r_coord
const Scalar r_coord(bool neumann, bool first) const
Expresses the areal radial coordinate by that in spatially isotropic coordinates.
Definition blackhole_r_coord.C:69
Lorene::Black_hole::mass_irr
virtual double mass_irr() const
Irreducible mass of the black hole.
Definition blackhole_global.C:72
Lorene::Black_hole::mp
Map & mp
Mapping associated with the black hole.
Definition blackhole.h:80
Lorene::Black_hole::lapconf
Scalar lapconf
A function (lapse function * conformal factor) lapconf generated by the black hole.
Definition blackhole.h:97
Lorene::Black_hole::bc_lapconf
const Valeur bc_lapconf(bool neumann, bool first) const
Boundary condition on the apparent horizon of the black hole for the lapse function: 2-D Valeur.
Definition blackhole_bc.C:74
Lorene::Black_hole::shift
Vector shift
Shift vector generated by the black hole.
Definition blackhole.h:109
Lorene::Black_hole::rad_ah
virtual double rad_ah() const
Radius of the apparent horizon.
Definition blackhole_global.C:508
Lorene::Black_hole::lapconf_rs
Scalar lapconf_rs
Part of lapconf from the numerical computation.
Definition blackhole.h:100
Lorene::Black_hole::equilibrium_spher
void equilibrium_spher(bool neumann, bool first, double spin_omega, double precis=1.e-14, double precis_shift=1.e-8)
Computes a spherical, static black-hole by giving boundary conditions on the apparent horizon.
Definition blackhole_eq_spher.C:73
Lorene::Black_hole::mass_adm
virtual double mass_adm() const
ADM mass.
Definition blackhole_global.C:106
Lorene::Black_hole::bc_shift_z
const Valeur bc_shift_z() const
Boundary condition on the apparent horizon of the black hole for the shift vector of the \fz\f direct...
Definition blackhole_bc.C:353
Lorene::Black_hole::lapconf_bh
Scalar lapconf_bh
Part of lapconf from the analytic background.
Definition blackhole.h:103
Lorene::Black_hole::bc_shift_y
const Valeur bc_shift_y(double omega_r) const
Boundary condition on the apparent horizon of the black hole for the shift vector of the \fy\f direct...
Definition blackhole_bc.C:284
Lorene::Black_hole::lapse
Scalar lapse
Lapse function generated by the black hole.
Definition blackhole.h:106
Lorene::Black_hole::kerrschild
bool kerrschild
true for a Kerr-Schild background, false for a conformally flat background
Definition blackhole.h:85
Lorene::Black_hole::confo
Scalar confo
Conformal factor generated by the black hole.
Definition blackhole.h:118
Lorene::Black_hole::mass_bh
double mass_bh
Gravitational mass of BH.
Definition blackhole.h:88
Lorene::Black_hole::extr_curv_bh
void extr_curv_bh()
Computes taij , taij_quad from shift , lapse , confo .
Definition blackhole_extr_curv.C:65
Lorene::Black_hole::del_deriv
virtual void del_deriv() const
Deletes all the derived quantities.
Definition blackhole.C:208
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:279
Lorene::Mg3d::get_angu
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition mg3d.C:604
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:465
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:393
Lorene::Scalar::poisson_neumann
Scalar poisson_neumann(const Valeur &, int) const
Idem as Scalar::poisson_dirichlet , the boundary condition being on the radial derivative of the solu...
Definition scalar_pde_frontiere.C:101
Lorene::Scalar::raccord
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Definition scalar_raccord.C:65
Lorene::Scalar::deriv
const Scalar & deriv(int i) const
Returns of *this , where .
Definition scalar_deriv.C:359
Lorene::Scalar::filtre
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
Definition scalar_manip.C:157
Lorene::Scalar::get_dzpuis
int get_dzpuis() const
Returns dzpuis.
Definition scalar.h:563
Lorene::Scalar::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:359
Lorene::Scalar::std_spectral_base
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition scalar.C:790
Lorene::Scalar::inc_dzpuis
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Definition scalar_r_manip.C:473
Lorene::Scalar::set_dzpuis
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:814
Lorene::Scalar::poisson_dirichlet
Scalar poisson_dirichlet(const Valeur &limite, int num) const
Is identicall to Scalar::poisson() .
Definition scalar_pde_frontiere.C:67
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case).
Definition tbl.h:281
Lorene::Tenseur
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:304
Lorene::Tenseur::set
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:830
Lorene::Tenseur::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tenseur.C:642
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition valeur.h:297
Lorene::Vector
Tensor field of valence 1.
Definition vector.h:188
Lorene::Vector::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:319
Lorene::Vector::set
Scalar & set(int)
Read/write access to a component.
Definition vector.C:299
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:223
Lorene::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:507
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:351
Lorene::Tensor::annule_domain
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:675
Lorene::Tensor::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition tensor.C:490
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition tenseur_operateur.C:282
Lorene
Lorene prototypes.
Definition app_hor.h:67
Unites
Standard units of space, time and mass.