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 char 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 
69 namespace Lorene {
70 void 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 ;
98  rr.std_spectral_base() ;
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) ;
158  lapconf_bh.raccord(1) ;
159 
160  lapconf = lapconf_bh ;
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.)) ;
221  lapse.annule_domain(0) ;
222  lapse.raccord(1) ;
223 
224  // confo = sqrt(r_are) ;
225  confo = sqrt(0.9*r_are) ;
227  confo.annule_domain(0) ;
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. ;
444  confo.annule_domain(0) ;
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 
551  shift.annule_domain(0) ;
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 }
Lorene::Scalar::raccord
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Definition: scalar_raccord.C:62
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:98
Lorene::Black_hole::extr_curv_bh
void extr_curv_bh()
Computes taij , taij_quad from shift , lapse , confo .
Definition: blackhole_extr_curv.C:62
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Lorene::Black_hole::mass_irr
virtual double mass_irr() const
Irreducible mass of the black hole.
Definition: blackhole_global.C:69
Lorene::Tenseur
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:301
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
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:470
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:70
Lorene::Map::cost
Coord cost
Definition: map.h:722
Lorene::Black_hole::lapconf_bh
Scalar lapconf_bh
Part of lapconf from the analytic background.
Definition: blackhole.h:103
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_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:211
Lorene::Black_hole::shift_bh
Vector shift_bh
Part of the shift vector from the analytic background.
Definition: blackhole.h:115
Lorene::Mg3d
Multi-domain grid.
Definition: grilles.h:273
Lorene::Tenseur::set
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
Lorene::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Mg3d::get_angu
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition: mg3d.C:473
Lorene::Black_hole::mass_kom
virtual double mass_kom() const
Komar mass.
Definition: blackhole_global.C:289
Lorene::Tensor::set_etat_zero
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition: tensor.C:497
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:281
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Black_hole::lapse
Scalar lapse
Lapse function generated by the black hole.
Definition: blackhole.h:106
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:151
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::mass_bh
double mass_bh
Gravitational mass of BH.
Definition: blackhole.h:88
Lorene::Map::sint
Coord sint
Definition: map.h:721
Lorene::Map::r
Coord r
r coordinate centered on the grid
Definition: map.h:718
Lorene::Black_hole::mp
Map & mp
Mapping associated with the black hole.
Definition: blackhole.h:80
Lorene::Black_hole::taij
Sym_tensor taij
Trace of the extrinsic curvature.
Definition: blackhole.h:127
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Lorene::Tensor::annule_domain
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:666
Lorene::Tensor::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition: tensor.C:481
Unites
Standard units of space, time and mass.
Lorene::Map::sinp
Coord sinp
Definition: map.h:723
Lorene::Scalar::poisson_dirichlet
Scalar poisson_dirichlet(const Valeur &limite, int num) const
Is identicall to Scalar::poisson() .
Definition: scalar_pde_frontiere.C:64
Lorene::Vector::set
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
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:350
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Lorene::Tbl
Basic array class.
Definition: tbl.h:161
Lorene::Black_hole::mass_adm
virtual double mass_adm() const
ADM mass.
Definition: blackhole_global.C:103
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Black_hole::taij_quad
Scalar taij_quad
Part of the scalar generated by .
Definition: blackhole.h:135
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:66
Lorene::Scalar::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:353
Lorene::Map::get_bvect_cart
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
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Lorene::Scalar::set_dzpuis
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:808
Lorene::Vector
Tensor field of valence 1.
Definition: vector.h:188
Lorene::Tenseur::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:636
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:784
Lorene::Black_hole::lapconf_rs
Scalar lapconf_rs
Part of lapconf from the numerical computation.
Definition: blackhole.h:100
Lorene::Black_hole::rad_ah
virtual double rad_ah() const
Radius of the apparent horizon.
Definition: blackhole_global.C:505
Lorene::Black_hole::del_deriv
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: blackhole.C:205
Lorene::Scalar::get_dzpuis
int get_dzpuis() const
Returns dzpuis.
Definition: scalar.h:557
Lorene::Black_hole::confo
Scalar confo
Conformal factor generated by the black hole.
Definition: blackhole.h:118
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:410
Lorene::Map::cosp
Coord cosp
Definition: map.h:724
Lorene::Black_hole::shift
Vector shift
Shift vector generated by the black hole.
Definition: blackhole.h:109
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::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:71
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition: tenseur_operateur.C:279
Lorene::Vector::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: vector.C:316