LORENE
star_bhns_equilibrium.C
1/*
2 * Method of class Star_bhns to compute an equilibrium configuration
3 * of a neutron star in a black hole-neutron star binary
4 *
5 * (see file star_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 star_bhns_equilibrium_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star_bhns/star_bhns_equilibrium.C,v 1.4 2014/10/13 08:53:40 j_novak Exp $" ;
30
31/*
32 * $Id: star_bhns_equilibrium.C,v 1.4 2014/10/13 08:53:40 j_novak Exp $
33 * $Log: star_bhns_equilibrium.C,v $
34 * Revision 1.4 2014/10/13 08:53:40 j_novak
35 * Lorene classes and functions now belong to the namespace Lorene.
36 *
37 * Revision 1.3 2014/10/06 15:13:16 j_novak
38 * Modified #include directives to use c++ syntax.
39 *
40 * Revision 1.2 2008/05/15 19:13:45 k_taniguchi
41 * Change of some parameters.
42 *
43 * Revision 1.1 2007/06/22 01:30:45 k_taniguchi
44 * *** empty log message ***
45 *
46 *
47 * $Header: /cvsroot/Lorene/C++/Source/Star_bhns/star_bhns_equilibrium.C,v 1.4 2014/10/13 08:53:40 j_novak Exp $
48 *
49 */
50
51// C++ headers
52//#include <>
53
54// C headers
55#include <cmath>
56
57// Lorene headers
58#include "star_bhns.h"
59#include "param.h"
60#include "cmp.h"
61#include "tenseur.h"
62#include "utilitaires.h"
63#include "unites.h"
64//#include "graphique.h"
65
66namespace Lorene {
67void Star_bhns::equilibrium_bhns(double ent_c, const double& mass_bh,
68 const double& sepa, bool kerrschild,
69 int mer, int mermax_ns, int mermax_potvit,
70 int mermax_poisson, int filter_r,
71 int filter_r_s, int filter_p_s,
72 double relax_poisson, double relax_potvit,
73 double thres_adapt, double resize_ns,
74 const Tbl& fact_resize, Tbl& diff) {
75
76 // Fundamental constants and units
77 // -------------------------------
78 using namespace Unites ;
79
80 // Initializations
81 // ---------------
82
83 const Mg3d* mg = mp.get_mg() ;
84 int nz = mg->get_nzone() ; // Total number of domain
85 // int nt = mg->get_nt(0) ;
86
87 // The following is required to initialize mp_prev as a Map_et:
88 Map_et& mp_et = dynamic_cast<Map_et&>(mp) ;
89
90 // Domain and radial indices of points at the surface of the star:
91 int l_b = nzet - 1 ;
92 int i_b = mg->get_nr(l_b) - 1 ;
93 int k_b ;
94 int j_b ;
95
96 // Value of the enthalpy defining the surface of the star
97 double ent_b = 0 ;
98
99 // Error indicators
100 // ----------------
101
102 double& diff_ent = diff.set(0) ;
103 double& diff_vel_pot = diff.set(1) ;
104 double& diff_lapconf = diff.set(2) ;
105 double& diff_confo = diff.set(3) ;
106 double& diff_shift_x = diff.set(4) ;
107 double& diff_shift_y = diff.set(5) ;
108 double& diff_shift_z = diff.set(6) ;
109 double& diff_dHdr = diff.set(7) ;
110 double& diff_dHdr_min = diff.set(8) ;
111 double& diff_phi_min = diff.set(9) ;
112 double& diff_radius = diff.set(10) ;
113
114 // Parameters for the function Map_et::adapt
115 // -----------------------------------------
116
118 int nitermax = 100 ;
119 int niter ;
120 int adapt_flag = 1 ; // 1 = performs the full computation,
121 // 0 = performs only the rescaling by
122 // the factor alpha_r
123
124 // Number of domains for searching the enthalpy isosurfaces
125 int nz_search = nzet + 1 ;
126 // int nz_search = nzet ;
127
128 double precis_secant = 1.e-14 ;
129 double alpha_r ;
130 double reg_map = 1. ; // 1 = regular mapping, 0 = contracting mapping
131
132 Tbl ent_limit(nz) ;
133
134 par_adapt.add_int(nitermax, 0) ; // maximum number of iterations to
135 // locate zeros by the secant method
136 par_adapt.add_int(nzet, 1) ; // number of domains where the
137 // adjustment to the isosurfaces of
138 // ent is to be performed
139 par_adapt.add_int(nz_search, 2) ; // number of domains to search for
140 // the enthalpy isosurface
141 par_adapt.add_int(adapt_flag, 3) ; // 1 = performs the full computation,
142 // 0 = performs only the rescaling by
143 // the factor alpha_r
144 par_adapt.add_int(j_b, 4) ; // theta index of the collocation point
145 // (theta_*, phi_*)
146 par_adapt.add_int(k_b, 5) ; // theta index of the collocation point
147 // (theta_*, phi_*)
148 par_adapt.add_int_mod(niter, 0) ; // number of iterations actually used
149 // in the secant method
150 par_adapt.add_double(precis_secant, 0) ; // required absolute precision in
151 // the determination of zeros by
152 // the secant method
153 par_adapt.add_double(reg_map, 1) ; // 1. = regular mapping,
154 // 0 = contracting mapping
155 par_adapt.add_double(alpha_r, 2) ; // factor by which all the radial
156 // distances will be multiplied
157 par_adapt.add_tbl(ent_limit, 0) ; // array of values of the field ent
158 // to define the isosurfaces
159
162
163 ssjm1lapconf.set_etat_qcq() ;
164 ssjm1confo.set_etat_qcq() ;
165
166 double precis_poisson = 1.e-14 ;
167
168 // Parameters for the function Scalar::poisson for lapconf_auto
169 // ------------------------------------------------------------
170
172
173 par_lapconf.add_int(mermax_poisson, 0) ; // maximum number of iterations
174 par_lapconf.add_double(relax_poisson, 0) ; // relaxation parameter
175 par_lapconf.add_double(precis_poisson, 1) ; // required precision
176 par_lapconf.add_int_mod(niter, 0) ; // number of iterations actually used
177 par_lapconf.add_cmp_mod( ssjm1lapconf ) ;
178
179 // Parameters for the function Scalar::poisson for confo_auto
180 // ----------------------------------------------------------
181
183
184 par_confo.add_int(mermax_poisson, 0) ; // maximum number of iterations
185 par_confo.add_double(relax_poisson, 0) ; // relaxation parameter
186 par_confo.add_double(precis_poisson, 1) ; // required precision
187 par_confo.add_int_mod(niter, 0) ; // number of iterations actually used
188 par_confo.add_cmp_mod( ssjm1confo ) ;
189
190 // Parameters for the function Vector::poisson for shift method 2
191 // --------------------------------------------------------------
192
194
195 par_shift2.add_int(mermax_poisson, 0) ; // maximum number of iterations
196 par_shift2.add_double(relax_poisson, 0) ; // relaxation parameter
197 par_shift2.add_double(precis_poisson, 1) ; // required precision
198 par_shift2.add_int_mod(niter, 0) ; // number of iterations actually used
199
201 ssjm1khi.set_etat_qcq() ;
202
204 ssjm1wshift.set_etat_qcq() ;
205 for (int i=0; i<3; i++) {
206 ssjm1wshift.set(i) = Cmp(ssjm1_wshift(i+1)) ;
207 }
208
209 par_shift2.add_cmp_mod(ssjm1khi) ;
210 par_shift2.add_tenseur_mod(ssjm1wshift) ;
211
212 Scalar ent_jm1 = ent ; // Enthalpy at previous step
213
214 Scalar lapconf_m1(mp) ; // = lapconf_auto - 0.5
215 Scalar confo_m1(mp) ; // = confo_auto - 0.5
216
217 Scalar source_lapconf(mp) ; // Source term in the equation for lapconf_auto
218 source_lapconf.set_etat_qcq() ;
219 Scalar source_confo(mp) ; // Source term in the equation for confo_auto
220 source_confo.set_etat_qcq() ;
221 Vector source_shift(mp, CON, mp.get_bvect_cart()) ; // Source term
222 // in the equation for shift_auto
223 source_shift.set_etat_qcq() ;
224
225 //===========================================================
226 // Start of iteration
227 //===========================================================
228
229 for (int mer_ns=0; mer_ns<mermax_ns; mer_ns++) {
230
231 cout << "-----------------------------------------------" << endl ;
232 cout << "step: " << mer_ns << endl ;
233 cout << "diff_ent = " << diff_ent << endl ;
234
235 //------------------------------------------------------
236 // Resolution of the elliptic equation for the velocity
237 // scalar potential
238 //------------------------------------------------------
239
240 if (irrotational) {
241 diff_vel_pot = velo_pot_bhns(mass_bh, sepa, kerrschild,
243 relax_potvit) ;
244 }
245 else {
246 diff_vel_pot = 0. ; // to avoid the warning
247 }
248
249 //-------------------------------------
250 // Computation of the new radial scale
251 //-------------------------------------
252
253 // alpha_r (r = alpha_r r') is determined so that the enthalpy
254 // takes the requested value ent_b at the stellar surface
255
256 // Values at the center of the star:
257 // lapconf_auto = lapconf_auto=0(r->infty) + 0.5
258 // The term lapconf_auto=0(r->infty) should be rescaled.
259 double lapconf_auto_c = lapconf_auto.val_grid_point(0,0,0,0) - 0.5 ;
260 double lapconf_comp_c = lapconf_comp.val_grid_point(0,0,0,0) ;
261
262 double confo_c = confo_tot.val_grid_point(0,0,0,0) ;
263
264 double gam_c = gam.val_grid_point(0,0,0,0) ;
265 double gam0_c = gam0.val_grid_point(0,0,0,0) ;
266
267 double hhh_c = exp(ent_c) ;
268 double hhh_b = exp(ent_b) ;
269
270 // Search for the reference point (theta_*, phi_*) [notation of
271 // Bonazzola, Gourgoulhon & Marck PRD 58, 104020 (1998)]
272 // at the surface of the star
273 // ------------------------------------------------------------
274 double alpha_r2 = 0. ;
275
276 for (int k=0; k<mg->get_np(l_b); k++) {
277 for (int j=0; j<mg->get_nt(l_b); j++) {
278
279 double lapconf_auto_b =
281 double lapconf_comp_b =
283
285
286 double gam_b = gam.val_grid_point(l_b,k,j,i_b) ;
287 double gam0_b = gam0.val_grid_point(l_b,k,j,i_b) ;
288
289 double aaa = (gam0_c*gam_b*hhh_b*confo_c)
291
292 // See Eq (100) from Gourgoulhon et al. (2001)
294 + 0.5 * (aaa - 1.))
296
297 if (alpha_r2_jk > alpha_r2) {
299 k_b = k ;
300 j_b = j ;
301 }
302 }
303 }
304
306
307 cout << "k_b, j_b, alpha_r: " << k_b << " " << j_b << " "
308 << alpha_r << endl ;
309
310 // New value of lapconf_auto
311 // -------------------------
312
313 lapconf_auto = alpha_r2 * (lapconf_auto - 0.5) + 0.5 ;
315 lapconf_tot_tmp.std_spectral_base() ;
316
317 /*
318 confo_auto = alpha_r2 * (confo_auto - 0.5) + 0.5 ;
319 Scalar confo_tot_tmp = confo_auto + confo_comp ;
320 confo_tot_tmp.std_spectral_base() ;
321 */
322 //------------------------------------------------------------
323 // Change the values of the inner points of the second domain
324 // by those of the outer points of the first domain
325 //------------------------------------------------------------
326
328
329 //--------------------------------------------
330 // First integral --> enthalpy in all space
331 // See Eq (98) from Gourgoulhon et al. (2001)
332 //--------------------------------------------
333
334 double log_lapconf_c = log(lapconf_tot_tmp.val_grid_point(0,0,0,0)) ;
335 double log_confo_c = log(confo_tot.val_grid_point(0,0,0,0)) ;
336 double loggam_c = loggam.val_grid_point(0,0,0,0) ;
337 double pot_centri_c = pot_centri.val_grid_point(0,0,0,0) ;
338
342
343
344 //----------------------------------------------------------
345 // Change the enthalpy field to be set its maximum position
346 // at the coordinate center
347 //----------------------------------------------------------
348
349 double dentdx = ent.dsdx().val_grid_point(0,0,0,0) ;
350 double dentdy = ent.dsdy().val_grid_point(0,0,0,0) ;
351
352 cout << "dH/dx|_center = " << dentdx << endl ;
353 cout << "dH/dy|_center = " << dentdy << endl ;
354
355 double dec_fact_x = 1. ;
356 double dec_fact_y = 1. ;
357
358 Scalar func_in(mp) ;
359 func_in = 1. - dec_fact_x * (dentdx/ent_c) * mp.x
360 - dec_fact_y * (dentdy/ent_c) * mp.y ;
361
362 func_in.annule(nzet, nz-1) ;
363 func_in.std_spectral_base() ;
364
365 Scalar func_ex(mp) ;
366 func_ex = 1. ;
367 func_ex.annule(0, nzet-1) ;
368 func_ex.std_spectral_base() ;
369
370 // New enthalpy field
371 // ------------------
372 ent = ent * (func_in + func_ex) ;
373
375
376 double dentdx_new = ent.dsdx().val_grid_point(0,0,0,0) ;
377 double dentdy_new = ent.dsdy().val_grid_point(0,0,0,0) ;
378 cout << "dH/dx|_new = " << dentdx_new << endl ;
379 cout << "dH/dy|_new = " << dentdy_new << endl ;
380
381 //-----------------------------------------------------
382 // Adaptation of the mapping to the new enthalpy field
383 //-----------------------------------------------------
384
385 double dent_eq = ent.dsdr().val_point(ray_eq_pi(),M_PI/2.,M_PI) ;
386 double dent_pole = ent.dsdr().val_point(ray_pole(),0.,0.) ;
387 double rap_dent = fabs( dent_eq / dent_pole ) ;
388 cout << "| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
390
391 if ( rap_dent < thres_adapt ) {
392 adapt_flag = 0 ; // No adaptation of the mapping
393 cout << "******* FROZEN MAPPING *********" << endl ;
394 }
395 else {
396 adapt_flag = 1 ; // The adaptation of the mapping is to be
397 // performed
398 }
399
400 ent_limit.set_etat_qcq() ;
401 for (int l=0; l<nzet; l++) { // loop on domains inside the star
402 ent_limit.set(l) = ent.val_grid_point(l, k_b, j_b, i_b) ;
403 }
404 ent_limit.set(nzet-1) = ent_b ;
405
407
408 Cmp ent_cmp(ent) ;
410 ent = ent_cmp ;
411
412 // Readjustment of the external boundary of domain l=nzet
413 // to keep a fixed ratio with respect to star's surface
414
415 double rr_in_1 = mp.val_r(1, -1., M_PI/2., 0.) ;
416
417 // Resize of the outer boundary of the shell including the BH
418 double rr_out_nm2 = mp.val_r(nz-2, 1., M_PI/2., 0.) ;
419 mp.resize(nz-2, rr_in_1/rr_out_nm2 * fact_resize(1)) ;
420
421 // Resize of the inner boundary of the shell including the BH
422 double rr_out_nm3 = mp.val_r(nz-3, 1., M_PI/2., 0.) ;
423 mp.resize(nz-3, rr_in_1/rr_out_nm3 * fact_resize(0)) ;
424
425 // if (mer % 2 == 0) {
426
427 if (nz > 4) {
428
429 // Resize of the domain 1
430 double rr_out_1 = mp.val_r(1, 1., M_PI/2., 0.) ;
432
433 if (nz > 5) {
434
435 // Resize of the domain from 2 to N-4
436 double rr_out_nm3_new = mp.val_r(nz-3, 1., M_PI/2., 0.) ;
437
438 for (int i=1; i<nz-4; i++) {
439
440 double rr_out_i = mp.val_r(i, 1., M_PI/2., 0.) ;
441
442 double rr_mid = rr_out_i
443 + (rr_out_nm3_new - rr_out_i) / double(nz - 3 - i) ;
444
445 double rr_2timesi = 2. * rr_out_i ;
446
447 if (rr_2timesi < rr_mid) {
448
449 double rr_out_ip1 = mp.val_r(i+1, 1., M_PI/2., 0.) ;
451
452 }
453 else {
454
455 double rr_out_ip1 = mp.val_r(i+1, 1., M_PI/2., 0.) ;
456 mp.resize(i+1, rr_mid / rr_out_ip1) ;
457
458 } // End of else
459
460 } // End of i loop
461
462 } // End of (nz > 5) loop
463
464 } // End of (nz > 4) loop
465
466 // } // End of (mer % 2) loop
467
468 //----------------------------------------------------
469 // Computation of the enthalpy at the new grid points
470 //----------------------------------------------------
471
472 mp_prev.homothetie(alpha_r) ;
473
474 Cmp ent_cmp2 (ent) ;
476 ent = ent_cmp2 ;
477
478 double ent_s_max = -1 ;
479 int k_s_max = -1 ;
480 int j_s_max = -1 ;
481
482 for (int k=0; k<mg->get_np(l_b); k++) {
483 for (int j=0; j<mg->get_nt(l_b); j++) {
484 double xx = fabs( ent.val_grid_point(l_b, k, j, i_b) ) ;
485 if (xx > ent_s_max) {
486 ent_s_max = xx ;
487 k_s_max = k ;
488 j_s_max = j ;
489 }
490 }
491 }
492 cout << "Max. abs(enthalpy) at the boundary between domains nzet-1"
493 << " and nzet : " << endl ;
494 cout << " " << ent_s_max << " reached for k = " << k_s_max
495 << " and j = " << j_s_max << endl ;
496
497 //-------------------
498 // Equation of state
499 //-------------------
500
501 equation_of_state() ; // computes new values for nbar (n), ener (e)
502 // and press (p) from the new ent (H)
503
504 //----------------------------------------------------------
505 // Matter source terms in the gravitational field equations
506 //----------------------------------------------------------
507
508 hydro_euler_bhns(kerrschild, mass_bh, sepa) ;
509 // computes new values for ener_euler (E),
510 // s_euler (S) and u_euler (U^i)
511
512
513 //-------------------------------------------------
514 // Computation of the minimum of the indicator chi
515 //-------------------------------------------------
516 double azimu_min = phi_min() ;
517 double rad_chi_min = radius_p(azimu_min) ;
519
520 cout << "| dH/dr_eq / dH/dr_pole | (minimum) = " << chi_min << endl ;
521 cout << " phi = " << azimu_min / M_PI << " [M_PI]" << endl ;
522 cout << " radius = " << rad_chi_min / km << " [km]" << endl ;
523
527
528 //-----------------------------------
529 // Poisson equation for lapconf_auto
530 //-----------------------------------
531
532 // Source
533 //--------
534
535 Scalar sou_lap1(mp) ; // dzpuis = 0
537 * (0.5*ener_euler + s_euler) ;
538
539 sou_lap1.std_spectral_base() ;
540 sou_lap1.annule(nzet,nz-1) ;
541 sou_lap1.inc_dzpuis(4) ; // dzpuis : 0 -> 4
542
543 Scalar sou_lap2(mp) ; // dzpuis = 4
544 sou_lap2 = 0.875 * (lapconf_auto+0.5) * taij_quad_auto
545 / pow(confo_auto+0.5,8.) ;
546 sou_lap2.std_spectral_base() ;
547
549
550 source_lapconf.std_spectral_base() ;
551 // source_lapse.annule(nzet,nz-1) ;
552
553 if (filter_r != 0) {
554 if (source_lapconf.get_etat() != ETATZERO) {
555 source_lapconf.filtre(filter_r) ;
556 // source_lapse.filtre_r(filter_r,0) ;
557 }
558 }
559
560 assert(source_lapconf.get_dzpuis() == 4) ;
561
562 // Resolution of the Poisson equation (Outer BC : lapconf_m1 -> 0)
563 // ----------------------------------
564
565 lapconf_m1.set_etat_qcq() ;
566 lapconf_m1 = lapconf_auto - 0.5 ;
569
570 // Check: has the Poisson equation been correctly solved ?
571 // -------------------------------------------------------
572
574 cout <<
575 "Relative error in the resolution of the equation for lapconf_auto : "
576 << endl ;
577 for (int l=0; l<nz; l++) {
578 cout << tdiff_lapconf(l) << " " ;
579 }
580 cout << endl ;
582
583 // Re-construction of the lapconf function
584 // ---------------------------------------
585 lapconf_auto = lapconf_m1 + 0.5 ;
586 // lapconf_tot = lapconf_auto + lapconf_comp
587 // lapconf_auto, _comp -> 0.5 (r -> inf)
588 // lapconf_tot -> 1 (r -> inf)
589
590 //---------------------------------
591 // Poisson equation for confo_auto
592 //---------------------------------
593
594 // Source
595 //--------
596
597 Scalar sou_con1(mp) ; // dzpuis = 0
598 sou_con1 = - 0.5 * qpig * pow(confo_tot,5.) * ener_euler ;
599 sou_con1.std_spectral_base() ;
600 sou_con1.annule(nzet,nz-1) ;
601 sou_con1.inc_dzpuis(4) ; // dzpuis : 0 -> 4
602
603 Scalar sou_con2(mp) ; // dzpuis = 4
604 sou_con2 = - 0.125 * taij_quad_auto / pow(confo_auto+0.5,7.) ;
605 sou_con2.std_spectral_base() ;
606
608
609 source_confo.std_spectral_base() ;
610 // source_confo.annule(nzet,nz-1) ;
611
612 if (filter_r != 0) {
613 if (source_confo.get_etat() != ETATZERO) {
614 source_confo.filtre(filter_r) ;
615 // source_confo.filtre_r(filter_r,0) ;
616 }
617 }
618
619 assert(source_confo.get_dzpuis() == 4) ;
620
621 // Resolution of the Poisson equation (Outer BC : confo_m1 -> 0)
622 // ----------------------------------
623
624 confo_m1.set_etat_qcq() ;
625 confo_m1 = confo_auto - 0.5 ;
628
629 // Check: has the Poisson equation been correctly solved ?
630 // -------------------------------------------------------
631
633 cout <<
634 "Relative error in the resolution of the equation for confo_auto : "
635 << endl ;
636 for (int l=0; l<nz; l++) {
637 cout << tdiff_confo(l) << " " ;
638 }
639 cout << endl ;
641
642 // Re-construction of the conformal factor
643 // ---------------------------------------
644 confo_auto = confo_m1 + 0.5 ; // confo_tot = confo_auto + confo_comp
645 // confo_auto, _comp -> 0.5 (r -> inf)
646 // confo_tot -> 1 (r -> inf)
647
648 //----------------------------------------
649 // Vector Poisson equation for shift_auto
650 //----------------------------------------
651
652 // Source
653 // ------
654
655 Vector sou_shif1(mp, CON, mp.get_bvect_cart()) ; // dzpuis = 0
656 sou_shif1.set_etat_qcq() ;
657
658 for (int i=1; i<=3; i++) {
659 sou_shif1.set(i) = 4.*qpig * lapconf_tot_tmp
660 * pow(confo_tot, 3.)
661 * (ener_euler + press) * u_euler(i) ;
662 }
663
664 sou_shif1.std_spectral_base() ;
665 sou_shif1.annule(nzet, nz-1) ;
666
667 for (int i=1; i<=3; i++) {
668 (sou_shif1.set(i)).inc_dzpuis(4) ; // dzpuis: 0 -> 4
669 }
670
671 Vector sou_shif2(mp, CON, mp.get_bvect_cart()) ; // dzpuis = 4
672 sou_shif2.set_etat_qcq() ;
673 for (int i=1; i<=3; i++) {
674 sou_shif2.set(i) = 2. *
676 -7.*(lapconf_auto+0.5)*d_confo_auto(1)
677 /(confo_auto+0.5))
679 -7.*(lapconf_auto+0.5)*d_confo_auto(2)
680 /(confo_auto+0.5))
682 -7.*(lapconf_auto+0.5)*d_confo_auto(3)
683 /(confo_auto+0.5))
684 ) / pow(confo_auto+0.5,7.) ;
685 }
686 sou_shif2.std_spectral_base() ;
687
689
690 source_shift.std_spectral_base() ;
691 // source_shift.annule(nzet, nz-1) ;
692
693 // Resolution of the Poisson equation
694 // ----------------------------------
695
696 // Filter for the source of shift vector
697
698 if (filter_r_s != 0) {
699 for (int i=1; i<=3; i++) {
700 if (source_shift(i).get_etat() != ETATZERO) {
701 source_shift.set(i).filtre(filter_r_s) ;
702 // source_shift.set(i).filtre_r(filter_r_s, 0) ;
703 }
704 }
705 }
706
707 if (filter_p_s != 0) {
708 for (int i=1; i<=3; i++) {
709 if (source_shift(i).get_etat() != ETATZERO) {
710 (source_shift.set(i)).filtre_phi(filter_p_s, nz-1) ;
711 // (source_shift.set(i)).filtre_phi(filter_p_s, 0) ;
712 }
713 }
714 }
715
716 for (int i=1; i<=3; i++) {
717 if(source_shift(i).dz_nonzero()) {
718 assert( source_shift(i).get_dzpuis() == 4 ) ;
719 }
720 else {
721 (source_shift.set(i)).set_dzpuis(4) ;
722 }
723 }
724
725 double lambda = double(1) / double(3) ;
726
727 Tenseur source_p(mp, 1, CON, mp.get_bvect_cart() ) ;
728 source_p.set_etat_qcq() ;
729 for (int i=0; i<3; i++) {
730 source_p.set(i) = Cmp(source_shift(i+1)) ;
731 }
732
733 Tenseur vect_auxi(mp, 1, CON, mp.get_bvect_cart()) ;
734 vect_auxi.set_etat_qcq() ;
735 for (int i=0; i<3 ;i++) {
736 vect_auxi.set(i) = 0. ;
737 }
739 scal_auxi.set_etat_qcq() ;
740 scal_auxi.set().annule_hard() ;
741 scal_auxi.set_std_base() ;
742
743 Tenseur resu_p(mp, 1, CON, mp.get_bvect_cart() ) ;
744 resu_p.set_etat_qcq() ;
745
746 source_p.poisson_vect(lambda, par_shift2, resu_p, vect_auxi,
747 scal_auxi) ;
748
749 for (int i=1; i<=3; i++)
750 shift_auto.set(i) = resu_p(i-1) ;
751
753
754 for (int i=0; i<3; i++) {
756 }
757
758 // Check: has the equation for shift_auto been correctly solved ?
759 // --------------------------------------------------------------
760
762 + lambda * shift_auto.divergence(flat).derive_con(flat) ;
763
764 source_shift.dec_dzpuis() ;
765
769
770 cout <<
771 "Relative error in the resolution of the equation for shift_auto : "
772 << endl ;
773 cout << "x component : " ;
774 for (int l=0; l<nz; l++) {
775 cout << tdiff_shift_x(l) << " " ;
776 }
777 cout << endl ;
778 cout << "y component : " ;
779 for (int l=0; l<nz; l++) {
780 cout << tdiff_shift_y(l) << " " ;
781 }
782 cout << endl ;
783 cout << "z component : " ;
784 for (int l=0; l<nz; l++) {
785 cout << tdiff_shift_z(l) << " " ;
786 }
787 cout << endl ;
788
792
793
794 //-----------------------------
795 // Relative change in enthalpy
796 //-----------------------------
797
800 for (int l=0; l<nzet; l++) {
802 }
803 diff_ent /= nzet ;
804
805 ent_jm1 = ent ;
806
807 /*
808 des_profile( lapconf_auto, 0., 10.,
809 M_PI/2., M_PI, "Self lapconf function of NS",
810 "Lapconf (theta=pi/2, phi=0)" ) ;
811
812 des_profile( lapconf_tot, 0., 10.,
813 M_PI/2., M_PI, "Total lapconf function seen by NS",
814 "Lapconf (theta=pi/2, phi=0)" ) ;
815
816 des_profile( confo_auto, 0., 10.,
817 M_PI/2., M_PI, "Self conformal factor of NS",
818 "Confo (theta=pi/2, phi=0)" ) ;
819
820 des_profile( confo_tot, 0., 10.,
821 M_PI/2., M_PI, "Total conformal factor seen by NS",
822 "Confo (theta=pi/2, phi=0)" ) ;
823
824 des_coupe_vect_z( shift_auto, 0., -3., 0.5, 3,
825 "Self shift vector of NS") ;
826
827 des_coupe_vect_z( shift_tot, 0., -3., 0.5, 3,
828 "Total shift vector seen by NS") ;
829 */
830 } // End of main loop
831
832 //=========================================================
833 // End of iteration
834 //=========================================================
835
836
837}
838}
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Radial mapping of rather general form.
Definition map.h:2752
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
Coord y
y coordinate centered on the grid
Definition map.h:727
virtual void reevaluate(const Map *mp_prev, int nzet, Cmp &uu) const =0
Recomputes the values of a Cmp at the collocation points after a change in the mapping.
virtual void resize(int l, double lambda)=0
Rescales the outer boundary of one domain.
virtual void adapt(const Cmp &ent, const Param &par, int nbr=0)=0
Adaptation of the mapping to a given scalar field.
Coord x
x coordinate centered on the grid
Definition map.h:726
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.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Multi-domain grid.
Definition grilles.h:273
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Parameter storage.
Definition param.h:125
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
const Scalar & dsdy() const
Returns of *this , where .
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
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition scalar.h:637
const Scalar & dsdx() const
Returns of *this , where .
const Scalar & dsdr() const
Returns of *this .
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the star ).
Definition star_bhns.h:192
Sym_tensor taij_auto
Part of the extrinsic curvature tensor generated by shift_auto , lapse_auto , and confo_auto .
Definition star_bhns.h:182
Scalar lapconf_auto
Lapconf function generated by the star.
Definition star_bhns.h:113
double velo_pot_bhns(const double &mass_bh, const double &sepa, bool kerrschild, int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
void equilibrium_bhns(double ent_c, const double &mass_bh, const double &sepa, bool kerrschild, int mer, int mermax_ns, int mermax_potvit, int mermax_poisson, int filter_r, int filter_r_s, int filter_p_s, double relax_poisson, double relax_potvit, double thres_adapt, double resize_ns, const Tbl &fact_resize, Tbl &diff)
Computes an equilibrium configuration.
Vector ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for by means...
Definition star_bhns.h:220
Scalar confo_tot
Total conformal factor.
Definition star_bhns.h:163
Vector d_lapconf_auto
Derivative of the lapconf function generated by the star .
Definition star_bhns.h:130
Vector d_confo_auto
Derivative of the conformal factor generated by the star .
Definition star_bhns.h:168
double radius_p(double phi)
Radius of the star to the direction of and .
Scalar gam0
Lorentz factor between the co-orbiting observer and the Eulerian one.
Definition star_bhns.h:107
void hydro_euler_bhns(bool kerrschild, const double &mass_bh, const double &sepa)
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Scalar gam
Lorentz factor between the fluid and the co-orbiting observer.
Definition star_bhns.h:102
Scalar ssjm1_lapconf
Effective source at the previous step for the resolution of the Poisson equation for lapconf_auto .
Definition star_bhns.h:197
bool irrotational
true for an irrotational star, false for a corotating one
Definition star_bhns.h:72
Scalar pot_centri
Centrifugal potential.
Definition star_bhns.h:110
double chi_rp(double radius, double phi)
Sensitive indicator of the mass-shedding to the direction of , , .
Scalar lapconf_comp
Lapconf function generated by the companion black hole.
Definition star_bhns.h:116
Scalar taij_quad_auto
Part of the scalar generated by .
Definition star_bhns.h:187
double phi_min()
Azimuthal angle when the indicator of the mass-shedding takes its minimum chi_min.
Vector shift_auto
Shift vector generated by the star.
Definition star_bhns.h:138
Scalar ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Definition star_bhns.h:210
Scalar confo_auto
Conformal factor generated by the star.
Definition star_bhns.h:157
Scalar ssjm1_confo
Effective source at the previous step for the resolution of the Poisson equation for confo_auto .
Definition star_bhns.h:202
Scalar loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition star_bhns.h:93
Scalar ener_euler
Total energy density in the Eulerian frame.
Definition star.h:198
void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition star.C:462
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Definition star.h:201
Scalar press
Fluid pressure.
Definition star.h:194
double ray_eq_pi() const
Coordinate radius at , [r_unit].
Scalar ent
Log-enthalpy.
Definition star.h:190
Vector u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition star.h:207
Map & mp
Mapping associated with the star.
Definition star.h:180
int nzet
Number of domains of *mp occupied by the star.
Definition star.h:183
double ray_pole() const
Coordinate radius at [r_unit].
Basic array class.
Definition tbl.h:161
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
void smooth(int nzet, Valeur &uuva) const
Changes the function *this as a smooth one when there exists a discontinuity between the nucleus and ...
Tensor field of valence 1.
Definition vector.h:188
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition vector.C:381
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:270
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:504
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Cmp abs(const Cmp &)
Absolute value.
Definition cmp_math.C:410
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition tensor.C:1014
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.