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
29char 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
68namespace Lorene {
69void 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) ;
388 r_are = r_coord(bc_lapconf_nd, bc_lapconf_fs) ;
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 ;
475 lapconf_auto_rs.annule_domain(0) ;
476 lapconf_auto_rs.raccord(1) ;
477
478 lapconf_auto = lapconf_auto_rs + lapconf_auto_bh ;
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
498 /lapconf_tot/lapconf_tot)
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 ;
540 confo_auto_rs.annule_domain(0) ;
541 confo_auto_rs.raccord(1) ;
542
543 confo_auto = confo_auto_rs + confo_auto_bh ;
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
708 shift_auto_rs.std_spectral_base() ;
709 shift_auto_rs.annule_domain(0) ;
710 for (int i=1; i<=3; i++) {
711 shift_auto_rs.set(i).raccord(1) ;
712 }
713
714 shift_auto = shift_auto_rs + shift_auto_bh ;
715 shift_auto.std_spectral_base() ;
716 shift_auto.annule_domain(0) ;
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::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::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::Black_hole::mp
Map & mp
Mapping associated with the black hole.
Definition blackhole.h:80
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::mass_bh
double mass_bh
Gravitational mass of BH.
Definition blackhole.h:88
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Lorene::Hole_bhns::confo_auto
Scalar confo_auto
Conformal factor generated by the black hole.
Definition hole_bhns.h:163
Lorene::Hole_bhns::lapconf_auto
Scalar lapconf_auto
Lapconf function generated by the black hole.
Definition hole_bhns.h:95
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::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::confo_auto_bh
Scalar confo_auto_bh
Part of the conformal factor from the analytic background.
Definition hole_bhns.h:160
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_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::taij_quad_auto
Scalar taij_quad_auto
Part of the scalar from the black hole.
Definition hole_bhns.h:238
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::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::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::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::Hole_bhns::confo_auto_rs
Scalar confo_auto_rs
Part of the conformal factor from the numerical computation.
Definition hole_bhns.h:157
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::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::Hole_bhns::shift_auto_bh
Vector shift_auto_bh
Part of the shift vector from the analytic background.
Definition hole_bhns.h:129
Lorene::Hole_bhns::confo_comp
Scalar confo_comp
Conformal factor generated by the companion star.
Definition hole_bhns.h:166
Lorene::Hole_bhns::lapconf_comp
Scalar lapconf_comp
Lapconf function generated by the companion star.
Definition hole_bhns.h:98
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::Hole_bhns::shift_auto
Vector shift_auto
Shift vector generated by the black hole.
Definition hole_bhns.h:132
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::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::Hole_bhns::taij_quad_comp
Scalar taij_quad_comp
Part of the scalar from the companion star.
Definition hole_bhns.h:241
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_tot
Scalar lapconf_tot
Total lapconf function.
Definition hole_bhns.h:101
Lorene::Hole_bhns::confo_tot
Scalar confo_tot
Total conformal factor.
Definition hole_bhns.h:169
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::Mg3d
Multi-domain grid.
Definition grilles.h:273
Lorene::Mg3d::get_angu
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition mg3d.C:473
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
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::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::filtre
void filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
Definition scalar_manip.C:151
Lorene::Scalar::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:353
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::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::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
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::Scalar::poisson_dirichlet
Scalar poisson_dirichlet(const Valeur &limite, int num) const
Is identicall to Scalar::poisson() .
Definition scalar_pde_frontiere.C:64
Lorene::Tbl
Basic array class.
Definition tbl.h:161
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::Tenseur::set
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
Lorene::Tenseur::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tenseur.C:636
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition valeur.h:287
Lorene::Vector
Tensor field of valence 1.
Definition vector.h:188
Lorene::Vector::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:316
Lorene::Vector::set
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Lorene::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:504
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
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition tenseur_operateur.C:279
Lorene
Lorene prototypes.
Definition app_hor.h:64
Unites
Standard units of space, time and mass.