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