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 
72 namespace Lorene {
73 void 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) ;
161  lapconf_bh.raccord(1) ;
162 
163  lapconf = lapconf_bh ;
165  lapconf_rs = 0. ;
167 
168  lapse = lapconf ;
170 
171  confo = 1. ;
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  }
184 
185  shift = shift_bh ;
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) ;
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.)) ;
224  lapse.annule_domain(0) ;
225  lapse.raccord(1) ;
226 
227  // confo = sqrt(r_are) ;
228  confo = sqrt(0.9*r_are) ;
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 
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. ;
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 
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 }
Map & mp
Mapping associated with the black hole.
Definition: blackhole.h:80
virtual double mass_irr() const
Irreducible mass of the black hole.
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:675
Scalar lapconf_rs
Part of lapconf from the numerical computation.
Definition: blackhole.h:100
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition: tensor.C:490
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Scalar taij_quad
Part of the scalar generated by .
Definition: blackhole.h:135
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
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: blackhole.C:208
void extr_curv_bh()
Computes taij , taij_quad from shift , lapse , confo .
Lorene prototypes.
Definition: app_hor.h:67
double mass_bh
Gravitational mass of BH.
Definition: blackhole.h:88
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
bool kerrschild
true for a Kerr-Schild background, false for a conformally flat background
Definition: blackhole.h:85
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
virtual double mass_adm() const
ADM mass.
Values and coefficients of a (real-value) function.
Definition: valeur.h:297
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:359
Tensor field of valence 1.
Definition: vector.h:188
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:507
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:814
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: vector.C:322
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Scalar confo
Conformal factor generated by the black hole.
Definition: blackhole.h:118
Coord sint
Definition: map.h:739
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 direction: ...
Definition: blackhole_bc.C:214
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...
Scalar lapconf_bh
Part of lapconf from the analytic background.
Definition: blackhole.h:103
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition: mg3d.C:604
Scalar poisson_neumann(const Valeur &, int) const
Idem as Scalar::poisson_dirichlet , the boundary condition being on the radial derivative of the solu...
Vector shift_bh
Part of the shift vector from the analytic background.
Definition: blackhole.h:115
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:563
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...
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 direction: ...
Definition: blackhole_bc.C:284
Coord sinp
Definition: map.h:741
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
virtual double rad_ah() const
Radius of the apparent horizon.
const Scalar & deriv(int i) const
Returns of *this , where .
Definition: scalar_deriv.C:359
Vector shift_rs
Part of the shift vector from the numerical computation.
Definition: blackhole.h:112
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
const Valeur bc_shift_z() const
Boundary condition on the apparent horizon of the black hole for the shift vector of the direction: ...
Definition: blackhole_bc.C:353
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
virtual double mass_kom() const
Komar mass.
Multi-domain grid.
Definition: grilles.h:279
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:809
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:663
Scalar lapse
Lapse function generated by the black hole.
Definition: blackhole.h:106
Vector shift
Shift vector generated by the black hole.
Definition: blackhole.h:109
Coord cosp
Definition: map.h:742
Sym_tensor taij
Trace of the extrinsic curvature.
Definition: blackhole.h:127
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition: tensor.C:506
Basic array class.
Definition: tbl.h:164
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:302
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:652
const Scalar r_coord(bool neumann, bool first) const
Expresses the areal radial coordinate by that in spatially isotropic coordinates. ...
Scalar poisson_dirichlet(const Valeur &limite, int num) const
Is identicall to Scalar::poisson() .
Scalar lapconf
A function (lapse function * conformal factor) lapconf generated by the black hole.
Definition: blackhole.h:97
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
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:304
Coord r
r coordinate centered on the grid
Definition: map.h:736
Coord cost
Definition: map.h:740