LORENE
hole_bhns_equilibrium.C
1 /*
2  * Method of class Hole_bhns to compute black-hole metric quantities
3  * in a black hole-neutron star binary
4  *
5  * (see file hole_bhns.h for documentation).
6  *
7  */
8 
9 /*
10  * Copyright (c) 2005-2007 Keisuke Taniguchi
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License version 2
16  * as published by the Free Software Foundation.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 
31 /*
32  * $Id: hole_bhns_equilibrium.C,v 1.5 2016/12/05 16:17:55 j_novak Exp $
33  * $Log: hole_bhns_equilibrium.C,v $
34  * Revision 1.5 2016/12/05 16:17:55 j_novak
35  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
36  *
37  * Revision 1.4 2014/10/13 08:53:00 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.3 2014/10/06 15:13:10 j_novak
41  * Modified #include directives to use c++ syntax.
42  *
43  * Revision 1.2 2008/05/15 19:05:12 k_taniguchi
44  * Change of some parameters.
45  *
46  * Revision 1.1 2007/06/22 01:24:36 k_taniguchi
47  * *** empty log message ***
48  *
49  *
50  * $Header: /cvsroot/Lorene/C++/Source/Hole_bhns/hole_bhns_equilibrium.C,v 1.5 2016/12/05 16:17:55 j_novak Exp $
51  *
52  */
53 
54 // C++ headers
55 //#include <>
56 
57 // C headers
58 #include <cmath>
59 
60 // Lorene headers
61 #include "hole_bhns.h"
62 #include "cmp.h"
63 #include "tenseur.h"
64 #include "param.h"
65 #include "eos.h"
66 #include "unites.h"
67 #include "proto.h"
68 #include "utilitaires.h"
69 //#include "graphique.h"
70 
71 namespace Lorene {
72 void Hole_bhns::equilibrium_bhns(int mer, int mermax_bh,
73  int filter_r, int filter_r_s, int filter_p_s,
74  double x_rot, double y_rot, double precis,
75  double omega_orb, double resize_bh,
76  const Tbl& fact_resize, Tbl& diff) {
77 
78  // Fundamental constants and units
79  // -------------------------------
80  using namespace Unites ;
81 
82  // Initializations
83  // ---------------
84 
85  const Mg3d* mg = mp.get_mg() ;
86  int nz = mg->get_nzone() ; // total number of domains
87 
88  // Re-adjustment of the boundary of domains
89  // ----------------------------------------
90 
91  double rr_in_1 = mp.val_r(1, -1., M_PI/2, 0.) ;
92 
93  /*
94  // Three shells outside the shell including NS
95  // -------------------------------------------
96 
97  // Resize of the outer boundary of the shell including the NS
98  double rr_out_nm5 = mp.val_r(nz-5, 1., M_PI/2., 0.) ;
99  mp.resize(nz-5, rr_in_1/rr_out_nm5 * fact_resize(1)) ;
100 
101  // Resize of the innner boundary of the shell including the NS
102  double rr_out_nm6 = mp.val_r(nz-6, 1., M_PI/2., 0.) ;
103  mp.resize(nz-6, rr_in_1/rr_out_nm6 * fact_resize(0)) ;
104 
105  if (mer % 2 == 0) {
106 
107  // Resize of the domain N-2
108  double rr_out_nm2 = mp.val_r(nz-2, 1., M_PI/2., 0.) ;
109  mp.resize(nz-2, 8. * rr_in_1 * fact_resize(1) / rr_out_nm2) ;
110 
111  // Resize of the domain N-3
112  double rr_out_nm3 = mp.val_r(nz-3, 1., M_PI/2., 0.) ;
113  mp.resize(nz-3, 4. * rr_in_1 * fact_resize(1) / rr_out_nm3) ;
114 
115  // Resize of the domain N-4
116  double rr_out_nm4 = mp.val_r(nz-4, 1., M_PI/2., 0.) ;
117  mp.resize(nz-4, 2. * rr_in_1 * fact_resize(1) / rr_out_nm4) ;
118 
119  if (nz > 7) {
120 
121  // Resize of the domain 1
122  double rr_out_1 = mp.val_r(1, 1., M_PI/2., 0.) ;
123  mp.resize(1, rr_in_1/rr_out_1 * resize_bh) ;
124 
125  if (nz > 8) {
126 
127  // Resize of the domain from 2 to N-7
128  double rr_out_1_new = mp.val_r(1, 1., M_PI/2., 0.) ;
129  double rr_out_nm6_new = mp.val_r(nz-6, 1., M_PI/2., 0.) ;
130  double dr = (rr_out_nm6_new - rr_out_1_new) / double(nz - 7) ;
131 
132  for (int i=1; i<nz-7; i++) {
133 
134  double rr = rr_out_1_new + i * dr ;
135  double rr_out_ip1 = mp.val_r(i+1, 1., M_PI/2., 0.) ;
136  mp.resize(i+1, rr/rr_out_ip1) ;
137 
138  }
139 
140  }
141 
142  }
143 
144  }
145  */
146 
147  /*
148  // Two shells outside the shell including NS
149  // -----------------------------------------
150 
151  // Resize of the outer boundary of the shell including the NS
152  double rr_out_nm4 = mp.val_r(nz-4, 1., M_PI/2., 0.) ;
153  mp.resize(nz-4, rr_in_1/rr_out_nm4 * fact_resize(1)) ;
154 
155  // Resize of the innner boundary of the shell including the NS
156  double rr_out_nm5 = mp.val_r(nz-5, 1., M_PI/2., 0.) ;
157  mp.resize(nz-5, rr_in_1/rr_out_nm5 * fact_resize(0)) ;
158 
159  // if (mer % 2 == 0) {
160 
161  // Resize of the domain N-2
162  double rr_out_nm2 = mp.val_r(nz-2, 1., M_PI/2., 0.) ;
163  mp.resize(nz-2, 3. * rr_in_1 * fact_resize(1) / rr_out_nm2) ;
164 
165  // Resize of the domain N-3
166  double rr_out_nm3 = mp.val_r(nz-3, 1., M_PI/2., 0.) ;
167  mp.resize(nz-3, 1.5 * rr_in_1 * fact_resize(1) / rr_out_nm3) ;
168 
169  if (nz > 6) {
170 
171  // Resize of the domain 1
172  double rr_out_1 = mp.val_r(1, 1., M_PI/2., 0.) ;
173  mp.resize(1, rr_in_1/rr_out_1 * resize_bh) ;
174 
175  if (nz > 7) {
176 
177  // Resize of the domain from 2 to N-6
178  double rr_out_nm5_new = mp.val_r(nz-5, 1., M_PI/2., 0.) ;
179 
180  for (int i=1; i<nz-6; i++) {
181 
182  double rr_out_i = mp.val_r(i, 1., M_PI/2., 0.) ;
183 
184  double rr_mid = rr_out_i
185  + (rr_out_nm5_new - rr_out_i) / double(nz - 5 - i) ;
186 
187  double rr_2timesi = 2. * rr_out_i ;
188 
189  if (rr_2timesi < rr_mid) {
190 
191  double rr_out_ip1 = mp.val_r(i+1, 1., M_PI/2., 0.) ;
192  mp.resize(i+1, rr_2timesi / rr_out_ip1) ;
193 
194  }
195  else {
196 
197  double rr_out_ip1 = mp.val_r(i+1, 1., M_PI/2., 0.) ;
198  mp.resize(i+1, rr_mid / rr_out_ip1) ;
199 
200  } // End of else
201 
202  } // End of i loop
203 
204  } // End of (nz > 7) loop
205 
206  } // End of (nz > 6) loop
207 
208  // } // End of (mer % 2) loop
209  */
210 
211  // One shell outside the shell including NS
212  // ----------------------------------------
213 
214  // Resize of the outer boundary of the shell including the NS
215  double rr_out_nm3 = mp.val_r(nz-3, 1., M_PI/2., 0.) ;
216  mp.resize(nz-3, rr_in_1/rr_out_nm3 * fact_resize(1)) ;
217 
218  // Resize of the innner boundary of the shell including the NS
219  double rr_out_nm4 = mp.val_r(nz-4, 1., M_PI/2., 0.) ;
220  mp.resize(nz-4, rr_in_1/rr_out_nm4 * fact_resize(0)) ;
221 
222  // if (mer % 2 == 0) {
223 
224  // Resize of the domain N-2
225  double rr_out_nm2 = mp.val_r(nz-2, 1., M_PI/2., 0.) ;
226  mp.resize(nz-2, 2. * rr_in_1 * fact_resize(1) / rr_out_nm2) ;
227 
228  if (nz > 5) {
229 
230  // Resize of the domain 1
231  double rr_out_1 = mp.val_r(1, 1., M_PI/2., 0.) ;
232  mp.resize(1, rr_in_1/rr_out_1 * resize_bh) ;
233 
234  if (nz > 6) {
235 
236  // Resize of the domain from 2 to N-5
237  double rr_out_nm4_new = mp.val_r(nz-4, 1., M_PI/2., 0.) ;
238 
239  for (int i=1; i<nz-5; i++) {
240 
241  double rr_out_i = mp.val_r(i, 1., M_PI/2., 0.) ;
242 
243  double rr_mid = rr_out_i
244  + (rr_out_nm4_new - rr_out_i) / double(nz - 4 - i) ;
245 
246  double rr_2timesi = 2. * rr_out_i ;
247 
248  if (rr_2timesi < rr_mid) {
249 
250  double rr_out_ip1 = mp.val_r(i+1, 1., M_PI/2., 0.) ;
251  mp.resize(i+1, rr_2timesi / rr_out_ip1) ;
252 
253  }
254  else {
255 
256  double rr_out_ip1 = mp.val_r(i+1, 1., M_PI/2., 0.) ;
257  mp.resize(i+1, rr_mid / rr_out_ip1) ;
258 
259  } // End of else
260 
261  } // End of i loop
262 
263  } // End of (nz > 6) loop
264 
265  } // End of (nz > 5) loop
266 
267  // } // End of (mer % 2) loop
268 
269 
270  // Inner boundary condition
271  // ------------------------
272 
273  Valeur bc_lapc(mg->get_angu()) ;
274  Valeur bc_conf(mg->get_angu()) ;
275 
276  Valeur bc_shif_x(mg->get_angu()) ;
277  Valeur bc_shif_y(mg->get_angu()) ;
278  Valeur bc_shif_z(mg->get_angu()) ;
279 
280  // Error indicators
281  // ----------------
282 
283  double& diff_lapconf = diff.set(0) ;
284  double& diff_confo = diff.set(1) ;
285  double& diff_shift_x = diff.set(2) ;
286  double& diff_shift_y = diff.set(3) ;
287  double& diff_shift_z = diff.set(4) ;
288 
289  Scalar lapconf_jm1 = lapconf_auto_rs ; // Lapconf function at previous step
290  Scalar confo_jm1 = confo_auto_rs ; // Conformal factor at preious step
291  Vector shift_jm1 = shift_auto_rs ; // Shift vector at previous step
292 
293  // Auxiliary quantities
294  // --------------------
295 
296  Scalar source_lapconf(mp) ;
297  Scalar source_confo(mp) ;
298  Vector source_shift(mp, CON, mp.get_bvect_cart()) ;
299 
300  Scalar lapconf_m1(mp) ; // = lapconf_auto_rs + 0.5
301  Scalar confo_m1(mp) ; // = confo_auto_rs + 0.5
302 
303  double mass = ggrav * mass_bh ;
304 
305  Scalar rr(mp) ;
306  rr = mp.r ;
307  rr.std_spectral_base() ;
308  Scalar st(mp) ;
309  st = mp.sint ;
310  st.std_spectral_base() ;
311  Scalar ct(mp) ;
312  ct = mp.cost ;
313  ct.std_spectral_base() ;
314  Scalar sp(mp) ;
315  sp = mp.sinp ;
316  sp.std_spectral_base() ;
317  Scalar cp(mp) ;
318  cp = mp.cosp ;
319  cp.std_spectral_base() ;
320 
321  Vector ll(mp, CON, mp.get_bvect_cart()) ;
322  ll.set_etat_qcq() ;
323  ll.set(1) = st % cp ;
324  ll.set(2) = st % sp ;
325  ll.set(3) = ct ;
326  ll.std_spectral_base() ;
327 
328  Vector dlappsi(mp, COV, mp.get_bvect_cart()) ;
329  for (int i=1; i<=3; i++) {
330  dlappsi.set(i) = lapconf_auto_rs.deriv(i) + d_lapconf_comp(i)
331  - 7. * lapconf_tot % (confo_auto_rs.deriv(i) + d_confo_comp(i))
332  / confo_tot ;
333  }
334 
335  dlappsi.std_spectral_base() ;
336 
337  //======================================//
338  // Start of iteration //
339  //======================================//
340 
341  for (int mer_bh=0; mer_bh<mermax_bh; mer_bh++) {
342 
343  cout << "--------------------------------------------------" << endl ;
344  cout << "step: " << mer_bh << endl ;
345  cout << "diff_lapconf = " << diff_lapconf << endl ;
346  cout << "diff_confo = " << diff_confo << endl ;
347  cout << "diff_shift : x = " << diff_shift_x
348  << " y = " << diff_shift_y << " z = " << diff_shift_z << endl ;
349 
350  if (kerrschild) {
351 
352  cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
353  abort() ;
354 
355  } // End of Kerr-Schild
356  else { // Isotropic coordinates with the maximal slicing
357 
358  // Sets C/M^2 for each case of the lapse boundary condition
359  // --------------------------------------------------------
360  double cc ;
361 
362  if (bc_lapconf_nd) { // Neumann boundary condition
363  if (bc_lapconf_fs) { // First condition
364  // d(\alpha \psi)/dr = 0
365  // ---------------------
366  cc = 2. * (sqrt(13.) - 1.) / 3. ;
367  }
368  else { // Second condition
369  // d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
370  // -----------------------------------------
371  cc = 4. / 3. ;
372  }
373  }
374  else { // Dirichlet boundary condition
375  if (bc_lapconf_fs) { // First condition
376  // (\alpha \psi) = 1/2
377  // -------------------
378  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
379  abort() ;
380  }
381  else { // Second condition
382  // (\alpha \psi) = 1/sqrt(2.) \psi_KS
383  // ----------------------------------
384  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
385  abort() ;
386  // cc = 2. * sqrt(2.) ;
387  }
388  }
389 
390  Scalar r_are(mp) ;
392  r_are.std_spectral_base() ;
393 
394  Scalar lapbh_iso(mp) ;
395  lapbh_iso = sqrt(1. - 2.*mass/r_are/rr
396  + cc*cc*pow(mass/r_are/rr,4.)) ;
397  lapbh_iso.std_spectral_base() ;
398  lapbh_iso.annule_domain(0) ;
399  lapbh_iso.raccord(1) ;
400 
401  Scalar psibh_iso(mp) ;
402  psibh_iso = sqrt(r_are) ;
403  psibh_iso.std_spectral_base() ;
404  psibh_iso.annule_domain(0) ;
405  psibh_iso.raccord(1) ;
406 
407  Scalar dlapbh_iso(mp) ;
408  dlapbh_iso = mass/r_are/rr - 2.*cc*cc*pow(mass/r_are/rr,4.) ;
409  dlapbh_iso.std_spectral_base() ;
410  dlapbh_iso.annule_domain(0) ;
411  dlapbh_iso.raccord(1) ;
412 
413  //---------------------------------------------------------------//
414  // Resolution of the Poisson equation for the lapconf function //
415  //---------------------------------------------------------------//
416 
417  // Source term
418  // -----------
419 
420  Scalar tmpl1(mp) ;
421  tmpl1 = 0.875 * lapconf_tot % taij_quad_auto
422  / pow(confo_tot, 8.) ;
423  tmpl1.std_spectral_base() ; // dzpuis = 4
424  tmpl1.annule_domain(0) ;
425  tmpl1.raccord(1) ;
426 
427  Scalar tmpl2(mp) ;
428  tmpl2 = 0.875 * (lapconf_comp+0.5) * taij_quad_comp
429  * (pow(confo_tot/(confo_comp+0.5),6.)*(lapconf_comp+0.5)
430  /lapconf_tot - 1.)
431  / pow(confo_comp+0.5,8.) ;
432  tmpl2.std_spectral_base() ;
433  tmpl2.annule_domain(0) ; // dzpuis = 4
434  tmpl2.raccord(1) ;
435 
436  Scalar tmpl3(mp) ;
437  tmpl3 = 5.25 * cc * cc * pow(mass,4.) * lapbh_iso
438  * (pow(confo_tot,6.)*lapbh_iso/lapconf_tot - pow(psibh_iso,5.))
439  / pow(r_are*rr,6.) ;
440  tmpl3.std_spectral_base() ;
441  tmpl3.annule_domain(0) ;
442  tmpl3.raccord(1) ;
443 
444  tmpl3.inc_dzpuis(4) ; // dzpuis : 0 -> 4
445 
446  source_lapconf = tmpl1 + tmpl2 + tmpl3 ;
447  source_lapconf.std_spectral_base() ;
448 
449  source_lapconf.annule_domain(0) ;
450  source_lapconf.raccord(1) ;
451  /*
452  if (source_lapconf.get_dzpuis() != 4) {
453  source_lapconf.set_dzpuis(4) ;
454  }
455  source_lapconf.std_spectral_base() ;
456  */
457  if (filter_r != 0) {
458  if (source_lapconf.get_etat() != ETATZERO) {
459  source_lapconf.filtre(filter_r) ;
460  }
461  }
462 
463  bc_lapc = bc_lapconf() ;
464 
465  lapconf_m1.set_etat_qcq() ;
466 
467  if (bc_lapconf_nd) {
468  lapconf_m1 = source_lapconf.poisson_neumann(bc_lapc, 0) ;
469  }
470  else {
471  lapconf_m1 = source_lapconf.poisson_dirichlet(bc_lapc, 0) ;
472  }
473 
474  // Re-construction of the lapconf function
475  // ---------------------------------------
476 
477  lapconf_auto_rs = lapconf_m1 - 0.5 ;
480 
482  lapconf_auto.annule_domain(0) ; // lapconf_auto,_comp->0.5 (r->inf)
483  lapconf_auto.raccord(1) ; // lapconf_tot -> 1 (r->inf)
484 
485 
486  //---------------------------------------------------------------//
487  // Resolution of the Poisson equation for the conformal factor //
488  //---------------------------------------------------------------//
489 
490  // Source term
491  // -----------
492 
493  Scalar tmpc1 = - 0.125 * taij_quad_auto / pow(confo_tot, 7.) ;
494  tmpc1.std_spectral_base() ; // dzpuis = 4
495  tmpc1.annule_domain(0) ;
496  tmpc1.raccord(1) ;
497 
498  Scalar tmpc2 = 0.75 * cc * cc * pow(mass,4.)
499  * (pow(psibh_iso,5.)
500  - pow(confo_tot,7.)*lapbh_iso*lapbh_iso
502  / pow(r_are*rr,6.) ;
503  tmpc2.std_spectral_base() ;
504  tmpc2.annule_domain(0) ;
505  tmpc2.raccord(1) ;
506 
507  tmpc2.inc_dzpuis(4) ; // dzpuis : 0 -> 4
508 
509  Scalar tmpc3 = 0.125 * taij_quad_comp
510  * (1. - pow(confo_tot/(confo_comp+0.5),7.)
511  *pow((lapconf_comp+0.5)/lapconf_tot,2.))
512  / pow(confo_comp+0.5, 7.) ;
513  tmpc3.std_spectral_base() ; // dzpuis = 4
514  tmpc3.annule_domain(0) ;
515  tmpc3.raccord(1) ;
516 
517  source_confo = tmpc1 + tmpc2 + tmpc3 ;
518  source_confo.std_spectral_base() ;
519 
520  source_confo.annule_domain(0) ;
521  source_confo.raccord(1) ;
522  /*
523  if (source_confo.get_dzpuis() != 4) {
524  source_confo.set_dzpuis(4) ;
525  }
526  source_confo.std_spectral_base() ;
527  */
528  if (filter_r != 0) {
529  if (source_confo.get_etat() != ETATZERO) {
530  source_confo.filtre(filter_r) ;
531  }
532  }
533 
534  bc_conf = bc_confo(omega_orb, x_rot, y_rot) ;
535 
536  confo_m1.set_etat_qcq() ;
537 
538  confo_m1 = source_confo.poisson_neumann(bc_conf, 0) ;
539 
540  // Re-construction of the conformal factor
541  // ---------------------------------------
542  confo_auto_rs = confo_m1 - 0.5 ;
545 
547  confo_auto.annule_domain(0) ; // confo_auto,_comp->0.5 (r->inf)
548  confo_auto.raccord(1) ; // confo_tot -> 1 (r->inf)
549 
550 
551  //-----------------------------------------------------------//
552  // Resolution of the Poisson equation for the shift vector //
553  //-----------------------------------------------------------//
554 
555  // Source term
556  // -----------
557 
558  Vector dlapconf(mp, COV, mp.get_bvect_cart()) ;
559  for (int i=1; i<=3; i++) {
560  dlapconf.set(i) = lapconf_auto_rs.deriv(i)
561  - 7. * lapconf_tot % confo_auto_rs.deriv(i)
562  / confo_tot ;
563  }
564 
565  dlapconf.std_spectral_base() ;
566 
567  Vector tmps1 = 2. * contract(taij_auto_rs, 1, dlappsi, 0)
568  / pow(confo_tot, 7.)
569  + 2. * contract(taij_comp, 1, dlapconf, 0)
570  * (lapconf_comp+0.5) / lapconf_tot / pow(confo_comp+0.5, 7.) ;
571  tmps1.std_spectral_base() ; // dzpuis = 4
572  tmps1.annule_domain(0) ;
573  for (int i=1; i<=3; i++) {
574  tmps1.set(i).raccord(1) ;
575  }
576 
577  Vector tmps2(mp, CON, mp.get_bvect_cart()) ;
578  tmps2.set_etat_qcq() ;
579  for (int i=1; i<=3; i++) {
580  tmps2.set(i) = 2. * psibh_iso
581  * (dlapbh_iso + 0.5*(lapbh_iso - 1.)
582  *(lapbh_iso - 7.*lapconf_tot/confo_tot))
583  * (taij_tot_rs(i,1)%ll(1) + taij_tot_rs(i,2)%ll(2)
584  + taij_tot_rs(i,3)%ll(3)) / pow(confo_tot,7.) / rr ;
585  }
586  tmps2.std_spectral_base() ;
587  tmps2.annule_domain(0) ;
588  for (int i=1; i<=3; i++) {
589  tmps2.set(i).raccord(1) ;
590  }
591  for (int i=1; i<=3; i++) {
592  (tmps2.set(i)).inc_dzpuis(2) ; // dzpuis : 2 -> 4
593  }
594 
595  Vector tmps3(mp, CON, mp.get_bvect_cart()) ;
596  tmps3.set_etat_qcq() ;
597  for (int i=1; i<=3; i++) {
598  tmps3.set(i) = 2. * cc * mass * mass * lapbh_iso
599  * (dlappsi(i) - 3.*ll(i)*(ll(1)%dlappsi(1)
600  + ll(2)%dlappsi(2)
601  + ll(3)%dlappsi(3)))
602  / lapconf_tot / pow(r_are*rr,3.) ;
603  }
604  tmps3.std_spectral_base() ;
605  tmps3.annule_domain(0) ;
606  for (int i=1; i<=3; i++) {
607  tmps3.set(i).raccord(1) ;
608  }
609  for (int i=1; i<=3; i++) {
610  (tmps3.set(i)).inc_dzpuis(2) ; // dzpuis : 2 -> 4
611  }
612 
613  Vector tmps4 = 4. * cc * mass * mass
614  * (dlapbh_iso * (1. - psibh_iso*lapbh_iso/lapconf_tot)
615  + 0.5 * lapbh_iso * (lapbh_iso - 1.)
616  * (6.*(psibh_iso/confo_tot - 1.)
617  + psibh_iso*(1./confo_tot - lapbh_iso/lapconf_tot)))
618  * ll / rr / pow(r_are*rr,3.) ;
619  tmps4.std_spectral_base() ;
620  tmps4.annule_domain(0) ;
621  for (int i=1; i<=3; i++) {
622  tmps4.set(i).raccord(1) ;
623  }
624  for (int i=1; i<=3; i++) {
625  (tmps4.set(i)).inc_dzpuis(4) ; // dzpuis : 0 -> 4
626  }
627 
628  Vector dlappsi_comp(mp, COV, mp.get_bvect_cart()) ;
629  for (int i=1; i<=3; i++) {
630  dlappsi_comp.set(i) = ((lapconf_comp+0.5)/lapconf_tot - 1.)
631  * d_lapconf_comp(i)
632  - 7. * (lapconf_comp+0.5) * ((confo_comp+0.5)/confo_tot - 1.)
633  * d_confo_comp(i) / (confo_comp+0.5) ;
634  }
635 
636  dlappsi_comp.std_spectral_base() ;
637 
638  Vector tmps5 = 2. * contract(taij_comp, 1, dlappsi_comp, 0)
639  / pow(confo_comp+0.5, 7.) ;
640  tmps5.std_spectral_base() ;
641  tmps5.annule_domain(0) ;
642  for (int i=1; i<=3; i++) {
643  tmps5.set(i).raccord(1) ;
644  }
645 
646  source_shift = tmps1 + tmps2 + tmps3 + tmps4 + tmps5 ;
647  source_shift.std_spectral_base() ;
648  source_shift.annule_domain(0) ;
649 
650  for (int i=1; i<=3; i++) {
651  source_shift.set(i).raccord(1) ;
652  }
653 
654  if (filter_r_s != 0) {
655  for (int i=1; i<=3; i++) {
656  if (source_shift(i).get_etat() != ETATZERO)
657  source_shift.set(i).filtre(filter_r_s) ;
658  }
659  }
660 
661  if (filter_p_s != 0) {
662  for (int i=1; i<=3; i++) {
663  if (source_shift(i).get_etat() != ETATZERO) {
664  source_shift.set(i).filtre_phi(filter_p_s, nz-1) ;
665  /*
666  for (int l=1; l<nz; l++) {
667  source_shift.set(i).filtre_phi(filter_p_s, l) ;
668  }
669  */
670  }
671  }
672  }
673 
674  /*
675  for (int i=1; i<=3; i++) {
676  if (source_shift(i).dz_nonzero()) {
677  assert( source_shift(i).get_dzpuis() == 4 ) ;
678  }
679  else {
680  (source_shift.set(i)).set_dzpuis(4) ;
681  }
682  }
683  */
684 
685  Tenseur source_p(mp, 1, CON, mp.get_bvect_cart()) ;
686  source_p.set_etat_qcq() ;
687  for (int i=0; i<3; i++) {
688  source_p.set(i) = Cmp(source_shift(i+1)) ;
689  }
690 
691  Tenseur resu_p(mp, 1, CON, mp.get_bvect_cart()) ;
692  resu_p.set_etat_qcq() ;
693 
694  for (int i=0; i<3; i++) {
695  resu_p.set(i) = Cmp(shift_auto_rs(i+1)) ;
696  }
697 
698  // Boundary condition
699  bc_shif_x = bc_shift_x(omega_orb, y_rot) ;
700  bc_shif_y = bc_shift_y(omega_orb, x_rot) ;
701  bc_shif_z = bc_shift_z() ;
702 
703  poisson_vect_frontiere(1./3., source_p, resu_p,
704  bc_shif_x, bc_shif_y, bc_shif_z,
705  0, precis, 7) ;
706 
707  for (int i=1; i<=3; i++) {
708  shift_auto_rs.set(i) = resu_p(i-1) ;
709  }
710 
713  for (int i=1; i<=3; i++) {
714  shift_auto_rs.set(i).raccord(1) ;
715  }
716 
720  for (int i=1; i<=3; i++) {
721  shift_auto.set(i).raccord(1) ;
722  }
723 
724  } // End of isotropic
725 
726  //------------------------------------------------//
727  // Relative difference in the metric quantities //
728  //------------------------------------------------//
729 
730  // Difference is calculated only outside the inner boundary.
731 
732  Tbl tdiff_lapconf = diffrel(lapconf_auto_rs, lapconf_jm1) ;
733  tdiff_lapconf.set(0) = 0. ;
734  cout << "Relative difference in the lapconf function : " << endl ;
735  for (int l=0; l<nz; l++) {
736  cout << tdiff_lapconf(l) << " " ;
737  }
738  cout << endl ;
739 
740  diff_lapconf = tdiff_lapconf(1) ;
741  for (int l=2; l<nz; l++) {
742  diff_lapconf += tdiff_lapconf(l) ;
743  }
744  diff_lapconf /= nz ;
745 
746  Tbl tdiff_confo = diffrel(confo_auto_rs, confo_jm1) ;
747  tdiff_confo.set(0) = 0. ;
748  cout << "Relative difference in the conformal factor : " << endl ;
749  for (int l=0; l<nz; l++) {
750  cout << tdiff_confo(l) << " " ;
751  }
752  cout << endl ;
753 
754  diff_confo = tdiff_confo(1) ;
755  for (int l=2; l<nz; l++) {
756  diff_confo += tdiff_confo(l) ;
757  }
758  diff_confo /= nz ;
759 
760  Tbl tdiff_shift_x = diffrel(shift_auto_rs(1), shift_jm1(1)) ;
761  tdiff_shift_x.set(0) = 0. ;
762  cout << "Relative difference in the shift vector (x) : " << endl ;
763  for (int l=0; l<nz; l++) {
764  cout << tdiff_shift_x(l) << " " ;
765  }
766  cout << endl ;
767 
768  diff_shift_x = tdiff_shift_x(1) ;
769  for (int l=2; l<nz; l++) {
770  diff_shift_x += tdiff_shift_x(l) ;
771  }
772  diff_shift_x /= nz ;
773 
774  Tbl tdiff_shift_y = diffrel(shift_auto_rs(2), shift_jm1(2)) ;
775  tdiff_shift_y.set(0) = 0. ;
776  cout << "Relative difference in the shift vector (y) : " << endl ;
777  for (int l=0; l<nz; l++) {
778  cout << tdiff_shift_y(l) << " " ;
779  }
780  cout << endl ;
781 
782  diff_shift_y = tdiff_shift_y(1) ;
783  for (int l=2; l<nz; l++) {
784  diff_shift_y += tdiff_shift_y(l) ;
785  }
786  diff_shift_y /= nz ;
787 
788  Tbl tdiff_shift_z = diffrel(shift_auto_rs(3), shift_jm1(3)) ;
789  tdiff_shift_z.set(0) = 0. ;
790  cout << "Relative difference in the shift vector (z) : " << endl ;
791  for (int l=0; l<nz; l++) {
792  cout << tdiff_shift_z(l) << " " ;
793  }
794  cout << endl ;
795 
796  diff_shift_z = tdiff_shift_z(1) ;
797  for (int l=2; l<nz; l++) {
798  diff_shift_z += tdiff_shift_z(l) ;
799  }
800  diff_shift_z /= nz ;
801 
802  /*
803  des_profile( lapconf_auto_rs, 0., 10.,
804  M_PI/2., 0., "Residual lapconf function of BH",
805  "Lapconf (theta=pi/2, phi=0)" ) ;
806 
807  des_profile( lapconf_auto_bh, 0., 10.,
808  M_PI/2., 0., "Analytic lapconf function of BH",
809  "Lapconf (theta=pi/2, phi=0)" ) ;
810 
811  des_profile( lapconf_auto, 0., 10.,
812  M_PI/2., 0., "Self lapconf function of BH",
813  "Lapconf (theta=pi/2, phi=0)" ) ;
814 
815  des_profile( lapconf_tot, 0., 10.,
816  M_PI/2., 0., "Total lapconf function of BH",
817  "Lapconf (theta=pi/2, phi=0)" ) ;
818 
819  des_profile( confo_auto_rs, 0., 10.,
820  M_PI/2., 0., "Residual conformal factor of BH",
821  "Confo (theta=pi/2, phi=0)" ) ;
822 
823  des_profile( confo_auto_bh, 0., 10.,
824  M_PI/2., 0., "Analytic conformal factor of BH",
825  "Confo (theta=pi/2, phi=0)" ) ;
826 
827  des_profile( confo_auto, 0., 10.,
828  M_PI/2., 0., "Self conformal factor of BH",
829  "Confo (theta=pi/2, phi=0)" ) ;
830 
831  des_profile( confo_tot, 0., 10.,
832  M_PI/2., 0., "Total conformal factor of BH",
833  "Confo (theta=pi/2, phi=0)" ) ;
834 
835  des_coupe_vect_z( shift_auto_rs, 0., -3., 0.5, 3,
836  "Residual shift vector of NS") ;
837 
838  des_coupe_vect_z( shift_auto_bh, 0., -3., 0.5, 3,
839  "Analytic shift vector of NS") ;
840 
841  des_coupe_vect_z( shift_auto, 0., -3., 0.5, 3,
842  "Self shift vector of NS") ;
843 
844  des_coupe_vect_z( shift_tot, 0., -3., 0.5, 3,
845  "Total Shift vector seen by NS") ;
846  */
847  } // End of main loop
848 
849  //====================================//
850  // End of iteration //
851  //====================================//
852 
853 }
854 }
Map & mp
Mapping associated with the black hole.
Definition: blackhole.h:80
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:675
Scalar lapconf_auto_bh
Part of the lapconf function from the analytic background.
Definition: hole_bhns.h:92
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_comp
Part of the scalar from the companion star.
Definition: hole_bhns.h:241
const Valeur bc_shift_z() const
Boundary condition on the apparent horizon of the black hole for the shift vector of the direction: ...
Definition: hole_bhns_bc.C:406
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
Scalar confo_auto
Conformal factor generated by the black hole.
Definition: hole_bhns.h:163
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
Vector shift_auto
Shift vector generated by the black hole.
Definition: hole_bhns.h:132
Sym_tensor taij_comp
Part of the extrinsic curvature tensor generated by the companion star.
Definition: hole_bhns.h:221
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
Values and coefficients of a (real-value) function.
Definition: valeur.h:297
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:560
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
Sym_tensor taij_auto_rs
Part of the extrinsic curvature tensor numericalty computed for the black hole.
Definition: hole_bhns.h:211
Scalar lapconf_comp
Lapconf function generated by the companion star.
Definition: hole_bhns.h:98
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: vector.C:322
const Valeur bc_shift_x(double ome_orb, double y_rot) const
Boundary condition on the apparent horizon of the black hole for the shift vector of the direction: ...
Definition: hole_bhns_bc.C:186
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Coord sint
Definition: map.h:739
const Valeur bc_shift_y(double ome_orb, double x_rot) const
Boundary condition on the apparent horizon of the black hole for the shift vector of the direction: ...
Definition: hole_bhns_bc.C:296
Scalar confo_auto_rs
Part of the conformal factor from the numerical computation.
Definition: hole_bhns.h:157
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
Vector d_lapconf_comp
Derivative of the lapconf function generated by the companion star.
Definition: hole_bhns.h:123
Vector shift_auto_bh
Part of the shift vector from the analytic background.
Definition: hole_bhns.h:129
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...
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition: mg3d.C:604
Scalar confo_comp
Conformal factor generated by the companion star.
Definition: hole_bhns.h:166
Scalar poisson_neumann(const Valeur &, int) const
Idem as Scalar::poisson_dirichlet , the boundary condition being on the radial derivative of the solu...
Scalar taij_quad_auto
Part of the scalar from the black hole.
Definition: hole_bhns.h:238
Vector shift_auto_rs
Part of the shift vector from the numerical computation.
Definition: hole_bhns.h:126
Coord sinp
Definition: map.h:741
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
const Scalar & deriv(int i) const
Returns of *this , where .
Definition: scalar_deriv.C:359
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
Definition: scalar_manip.C:157
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
virtual void resize(int l, double lambda)=0
Rescales the outer boundary of one domain.
Scalar lapconf_auto
Lapconf function generated by the black hole.
Definition: hole_bhns.h:95
Scalar confo_tot
Total conformal factor.
Definition: hole_bhns.h:169
Multi-domain grid.
Definition: grilles.h:279
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:724
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
Vector d_confo_comp
Derivative of the conformal factor generated by the companion star.
Definition: hole_bhns.h:185
Sym_tensor taij_tot_rs
Part of the extrinsic curvature tensor from the numerical computation.
Definition: hole_bhns.h:190
Coord cosp
Definition: map.h:742
Scalar lapconf_auto_rs
Part of the lapconf function from the numerical computation.
Definition: hole_bhns.h:89
bool bc_lapconf_fs
true for the first type BC for the lapconf function, false for the second type BH ...
Definition: hole_bhns.h:78
bool bc_lapconf_nd
true for the Neumann type BC for the lapconf function, false for the Dirichlet type BH ...
Definition: hole_bhns.h:73
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() .
const Valeur bc_lapconf() const
Boundary condition on the apparent horizon of the black hole for the lapconf function: 2-D Valeur...
Definition: hole_bhns_bc.C:71
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
Scalar lapconf_tot
Total lapconf function.
Definition: hole_bhns.h:101
Coord r
r coordinate centered on the grid
Definition: map.h:736
void equilibrium_bhns(int mer, int mermax_bh, int filter_r, int filter_r_s, int filter_p_s, double x_rot, double y_rot, double precis, double omega_orb, double resize_bh, const Tbl &fact_resize, Tbl &diff)
Computes a black-hole part in a black hole-neutron star binary by giving boundary conditions on the a...
Scalar confo_auto_bh
Part of the conformal factor from the analytic background.
Definition: hole_bhns.h:160
Coord cost
Definition: map.h:740