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