LORENE
et_bin_bhns_extr_equil.C
1/*
2 * Method of class Et_bin_bhns_extr to compute an equilibrium configuration
3 * of a BH-NS binary system with an extreme mass ratio
4 *
5 * (see file et_bin_bhns_extr.h for documentation).
6 *
7 */
8
9/*
10 * Copyright (c) 2004-2005 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 et_bin_bhns_extr_equil_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_bhns_extr_equil.C,v 1.11 2014/10/13 08:52:54 j_novak Exp $" ;
30
31/*
32 * $Id: et_bin_bhns_extr_equil.C,v 1.11 2014/10/13 08:52:54 j_novak Exp $
33 * $Log: et_bin_bhns_extr_equil.C,v $
34 * Revision 1.11 2014/10/13 08:52:54 j_novak
35 * Lorene classes and functions now belong to the namespace Lorene.
36 *
37 * Revision 1.10 2014/10/06 15:13:07 j_novak
38 * Modified #include directives to use c++ syntax.
39 *
40 * Revision 1.9 2005/02/28 23:11:05 k_taniguchi
41 * Modification to include the case of the conformally flat background metric
42 *
43 * Revision 1.8 2005/01/25 17:33:19 k_taniguchi
44 * Suppression of the filter for the source term of the shift vector.
45 *
46 * Revision 1.7 2005/01/03 18:01:12 k_taniguchi
47 * Addition of the method to fix the position of the neutron star
48 * in the coordinate system.
49 *
50 * Revision 1.6 2004/12/29 16:29:55 k_taniguchi
51 * Suppression of "dzpius" for the shift vector.
52 *
53 * Revision 1.5 2004/12/22 18:26:53 k_taniguchi
54 * Change an argument of poisson_vect_falloff.
55 *
56 * Revision 1.4 2004/12/06 17:59:50 k_taniguchi
57 * Change the position of resize.
58 *
59 * Revision 1.3 2004/12/02 21:31:56 k_taniguchi
60 * Set a filter for the shift vector.
61 *
62 * Revision 1.2 2004/12/02 15:05:36 k_taniguchi
63 * Modification of the procedure for resize.
64 *
65 * Revision 1.1 2004/11/30 20:48:45 k_taniguchi
66 * *** empty log message ***
67 *
68 *
69 * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_bhns_extr_equil.C,v 1.11 2014/10/13 08:52:54 j_novak Exp $
70 *
71 */
72
73// C headers
74#include <cmath>
75
76// Lorene headers
77#include "et_bin_bhns_extr.h"
78#include "etoile.h"
79#include "map.h"
80#include "coord.h"
81#include "param.h"
82#include "eos.h"
83#include "graphique.h"
84#include "utilitaires.h"
85#include "unites.h"
86
87namespace Lorene {
88void Et_bin_bhns_extr::equil_bhns_extr_ks(double ent_c, const double& mass,
89 const double& sepa, int mermax,
90 int mermax_poisson,
91 double relax_poisson,
92 int mermax_potvit,
93 double relax_potvit, int np_filter,
94 double thres_adapt,
95 Tbl& diff) {
96
97 // Fundamental constants and units
98 // -------------------------------
99 using namespace Unites ;
100
101 assert( kerrschild == true ) ;
102
103 // Initializations
104 // ---------------
105
106 const Mg3d* mg = mp.get_mg() ;
107 int nz = mg->get_nzone() ; // total number of domains
108
109 // The following is required to initialize mp_prev as a Map_et:
110 Map_et& mp_et = dynamic_cast<Map_et&>(mp) ;
111
112 // Domain and radial indices of points at the surface of the star:
113 int l_b = nzet - 1 ;
114 int i_b = mg->get_nr(l_b) - 1 ;
115 int k_b ;
116 int j_b ;
117
118 // Value of the enthalpy defining the surface of the star
119 double ent_b = 0 ;
120
121 // Error indicators
122 // ----------------
123
124 double& diff_ent = diff.set(0) ;
125 double& diff_vel_pot = diff.set(1) ;
126 double& diff_logn = diff.set(2) ;
127 double& diff_beta = diff.set(3) ;
128 double& diff_shift_x = diff.set(4) ;
129 double& diff_shift_y = diff.set(5) ;
130 double& diff_shift_z = diff.set(6) ;
131
132 // Parameters for the function Map_et::adapt
133 // -----------------------------------------
134
135 Param par_adapt ;
136 int nitermax = 100 ;
137 int niter ;
138 int adapt_flag = 1 ; // 1 = performs the full computation,
139 // 0 = performs only the rescaling by
140 // the factor alpha_r
141 int nz_search = nzet ; // Number of domains for searching
142 // the enthalpy isosurfaces
143 double precis_secant = 1.e-14 ;
144 double alpha_r ;
145 double reg_map = 1. ; // 1 = regular mapping, 0 = contracting mapping
146
147 Tbl ent_limit(nz) ;
148
149 par_adapt.add_int(nitermax, 0) ; // maximum number of iterations to
150 // locate zeros by the secant method
151 par_adapt.add_int(nzet, 1) ; // number of domains where the adjustment
152 // to the isosurfaces of ent is to be
153 // performed
154 par_adapt.add_int(nz_search, 2) ; // number of domains to search for
155 // the enthalpy isosurface
156 par_adapt.add_int(adapt_flag, 3) ; // 1 = performs the full computation,
157 // 0 = performs only the rescaling by
158 // the factor alpha_r
159 par_adapt.add_int(j_b, 4) ; // theta index of the collocation point
160 // (theta_*, phi_*)
161 par_adapt.add_int(k_b, 5) ; // theta index of the collocation point
162 // (theta_*, phi_*)
163 par_adapt.add_int_mod(niter, 0) ; // number of iterations actually
164 // used in the secant method
165 par_adapt.add_double(precis_secant, 0) ; // required absolute precision in
166 // the determination of zeros by
167 // the secant method
168 par_adapt.add_double(reg_map, 1) ; // 1 = regular mapping,
169 // 0 = contracting mapping
170 par_adapt.add_double(alpha_r, 2) ; // factor by which all the radial
171 // distances will be multiplied
172 par_adapt.add_tbl(ent_limit, 0) ; // array of values of the field ent
173 // to define the isosurfaces
174
175 // Parameters for the function Map_et::poisson for logn_auto
176 // ---------------------------------------------------------
177
178 double precis_poisson = 1.e-16 ;
179
180 Param par_poisson1 ;
181
182 par_poisson1.add_int(mermax_poisson, 0) ; // maximum number of iterations
183 par_poisson1.add_double(relax_poisson, 0) ; // relaxation parameter
184 par_poisson1.add_double(precis_poisson, 1) ; // required precision
185 par_poisson1.add_int_mod(niter, 0) ; // number of iterations actually
186 // used
187 par_poisson1.add_cmp_mod( ssjm1_logn ) ;
188
189 // Parameters for the function Map_et::poisson for beta_auto
190 // ---------------------------------------------------------
191
192 Param par_poisson2 ;
193
194 par_poisson2.add_int(mermax_poisson, 0) ; // maximum number of iterations
195 par_poisson2.add_double(relax_poisson, 0) ; // relaxation parameter
196 par_poisson2.add_double(precis_poisson, 1) ; // required precision
197 par_poisson2.add_int_mod(niter, 0) ; // number of iterations actually
198 // used
199 par_poisson2.add_cmp_mod( ssjm1_beta ) ;
200
201 // Parameters for the function Tenseur::poisson_vect
202 // -------------------------------------------------
203
204 Param par_poisson_vect ;
205
206 par_poisson_vect.add_int(mermax_poisson, 0) ; // maximum number of
207 // iterations
208 par_poisson_vect.add_double(relax_poisson, 0) ; // relaxation parameter
209 par_poisson_vect.add_double(precis_poisson, 1) ; // required precision
210 par_poisson_vect.add_cmp_mod( ssjm1_khi ) ;
211 par_poisson_vect.add_tenseur_mod( ssjm1_wshift ) ;
212 par_poisson_vect.add_int_mod(niter, 0) ;
213
214 // External potential
215 // See Eq (99) from Gourgoulhon et al. (2001)
216 // -----------------------------------------
217
218 Tenseur pot_ext = logn_comp + pot_centri + loggam ;
219
220 Tenseur ent_jm1 = ent ; // Enthalpy at previous step
221
222 Tenseur source(mp) ; // source term in the equation for logn_auto
223 // and beta_auto
224
225 Tenseur source_shift(mp, 1, CON, ref_triad) ; // source term in the
226 // equation for shift_auto
227
228 //==========================================================//
229 // Start of iteration //
230 //==========================================================//
231
232 for(int mer=0 ; mer<mermax ; mer++ ) {
233
234 cout << "-----------------------------------------------" << endl ;
235 cout << "step: " << mer << endl ;
236 cout << "diff_ent = " << diff_ent << endl ;
237
238 //------------------------------------------------------
239 // Resolution of the elliptic equation for the velocity
240 // scalar potential
241 //------------------------------------------------------
242
243 if (irrotational) {
244 diff_vel_pot = velocity_pot_extr(mass, sepa, mermax_potvit,
245 precis_poisson, relax_potvit) ;
246 }
247
248 //-------------------------------------
249 // Computation of the new radial scale
250 //-------------------------------------
251
252 // alpha_r (r = alpha_r r') is determined so that the enthalpy
253 // takes the requested value ent_b at the stellar surface
254
255 // Values at the center of the star:
256 double logn_auto_c = logn_auto()(0, 0, 0, 0) ;
257 double pot_ext_c = pot_ext()(0, 0, 0, 0) ;
258
259 // Search for the reference point (theta_*, phi_*)
260 // [notation of Bonazzola, Gourgoulhon & Marck PRD 58, 104020 (1998)]
261 // at the surface of the star
262 // ------------------------------------------------------------------
263 double alpha_r2 = 0 ;
264 for (int k=0; k<mg->get_np(l_b); k++) {
265 for (int j=0; j<mg->get_nt(l_b); j++) {
266
267 double pot_ext_b = pot_ext()(l_b, k, j, i_b) ;
268 double logn_auto_b = logn_auto()(l_b, k, j, i_b) ;
269
270 // See Eq (100) from Gourgoulhon et al. (2001)
271 double alpha_r2_jk = ( ent_c - ent_b + pot_ext_c - pot_ext_b)
272 / ( logn_auto_b - logn_auto_c ) ;
273
274 if (alpha_r2_jk > alpha_r2) {
275 alpha_r2 = alpha_r2_jk ;
276 k_b = k ;
277 j_b = j ;
278 }
279
280 }
281 }
282
283 alpha_r = sqrt(alpha_r2) ;
284
285 cout << "k_b, j_b, alpha_r: " << k_b << " " << j_b << " "
286 << alpha_r << endl ;
287
288 // New value of logn_auto
289 // ----------------------
290
291 logn_auto = alpha_r2 * logn_auto ;
292 logn_auto_regu = alpha_r2 * logn_auto_regu ;
293 logn_auto_c = logn_auto()(0, 0, 0, 0) ;
294
295 //------------------------------------------------------------
296 // Change the values of the inner points of the second domain
297 // by those of the outer points of the first domain
298 //------------------------------------------------------------
299
300 (logn_auto().va).smooth(nzet, (logn_auto.set()).va) ;
301
302 //--------------------------------------------
303 // First integral --> enthalpy in all space
304 // See Eq (98) from Gourgoulhon et al. (2001)
305 //--------------------------------------------
306
307 ent = (ent_c + logn_auto_c + pot_ext_c) - logn_auto - pot_ext ;
308
309 //----------------------------------------------------------
310 // Change the enthalpy field to be set its maximum position
311 // at the coordinate center
312 //----------------------------------------------------------
313
314 double dentdx = ent().dsdx()(0, 0, 0, 0) ;
315 double dentdy = ent().dsdy()(0, 0, 0, 0) ;
316
317 cout << "dH/dx|_center = " << dentdx << endl ;
318 cout << "dH/dy|_center = " << dentdy << endl ;
319
320 double dec_fact = 1. ;
321
322 Tenseur func_in(mp) ;
323 func_in.set_etat_qcq() ;
324 func_in.set() = 1. - dec_fact * (dentdx/ent_c) * mp.x
325 - dec_fact * (dentdy/ent_c) * mp.y ;
326 func_in.set().annule(nzet, nz-1) ;
327 func_in.set_std_base() ;
328
329 Tenseur func_ex(mp) ;
330 func_ex.set_etat_qcq() ;
331 func_ex.set() = 1. ;
332 func_ex.set().annule(0, nzet-1) ;
333 func_ex.set_std_base() ;
334
335 // New enthalpy field
336 // ------------------
337 ent.set() = ent() * (func_in() + func_ex()) ;
338
339 (ent().va).smooth(nzet, (ent.set()).va) ;
340
341 double dentdx_new = ent().dsdx()(0, 0, 0, 0) ;
342 double dentdy_new = ent().dsdy()(0, 0, 0, 0) ;
343 cout << "dH/dx|_new = " << dentdx_new << endl ;
344 cout << "dH/dy|_new = " << dentdy_new << endl ;
345
346 //-----------------------------------------------------
347 // Adaptation of the mapping to the new enthalpy field
348 //----------------------------------------------------
349
350 // Shall the adaptation be performed (cusp) ?
351 // ------------------------------------------
352
353 double dent_eq = ent().dsdr().val_point(ray_eq_pi(),M_PI/2.,M_PI) ;
354 double dent_pole = ent().dsdr().val_point(ray_pole(),0.,0.) ;
355 double rap_dent = fabs( dent_eq / dent_pole ) ;
356 cout << "| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
357
358 if ( rap_dent < thres_adapt ) {
359 adapt_flag = 0 ; // No adaptation of the mapping
360 cout << "******* FROZEN MAPPING *********" << endl ;
361 }
362 else{
363 adapt_flag = 1 ; // The adaptation of the mapping is to be
364 // performed
365 }
366
367 ent_limit.set_etat_qcq() ;
368 for (int l=0; l<nzet; l++) { // loop on domains inside the star
369 ent_limit.set(l) = ent()(l, k_b, j_b, i_b) ;
370 }
371
372 ent_limit.set(nzet-1) = ent_b ;
373 Map_et mp_prev = mp_et ;
374
375 mp.adapt(ent(), par_adapt) ;
376
377 //----------------------------------------------------
378 // Computation of the enthalpy at the new grid points
379 //----------------------------------------------------
380
381 mp_prev.homothetie(alpha_r) ;
382
383 mp.reevaluate(&mp_prev, nzet+1, ent.set()) ;
384
385 double fact_resize = 1. / alpha_r ;
386 for (int l=nzet; l<nz-1; l++) {
387 mp_et.resize(l, fact_resize) ;
388 }
389 mp_et.resize_extr(fact_resize) ;
390
391 double ent_s_max = -1 ;
392 int k_s_max = -1 ;
393 int j_s_max = -1 ;
394 for (int k=0; k<mg->get_np(l_b); k++) {
395 for (int j=0; j<mg->get_nt(l_b); j++) {
396 double xx = fabs( ent()(l_b, k, j, i_b) ) ;
397 if (xx > ent_s_max) {
398 ent_s_max = xx ;
399 k_s_max = k ;
400 j_s_max = j ;
401 }
402 }
403 }
404 cout << "Max. abs(enthalpy) at the boundary between domains nzet-1"
405 << " and nzet : " << endl ;
406 cout << " " << ent_s_max << " reached for k = " << k_s_max
407 << " and j = " << j_s_max << endl ;
408
409 //-------------------
410 // Equation of state
411 //-------------------
412
413 equation_of_state() ; // computes new values for nbar (n), ener (e)
414 // and press (p) from the new ent (H)
415
416 //----------------------------------------------------------
417 // Matter source terms in the gravitational field equations
418 //---------------------------------------------------------
419
420 hydro_euler_extr(mass, sepa) ; // computes new values for
421 // ener_euler (E), s_euler (S)
422 // and u_euler (U^i)
423
424
425 //----------------------------------------------------
426 // Auxiliary terms for the source of Poisson equation
427 //----------------------------------------------------
428
429 const Coord& xx = mp.x ;
430 const Coord& yy = mp.y ;
431 const Coord& zz = mp.z ;
432
433 Tenseur r_bh(mp) ;
434 r_bh.set_etat_qcq() ;
435 r_bh.set() = pow( (xx+sepa)*(xx+sepa) + yy*yy + zz*zz, 0.5) ;
436 r_bh.set_std_base() ;
437
438 Tenseur xx_cov(mp, 1, COV, ref_triad) ;
439 xx_cov.set_etat_qcq() ;
440 xx_cov.set(0) = xx + sepa ;
441 xx_cov.set(1) = yy ;
442 xx_cov.set(2) = zz ;
443 xx_cov.set_std_base() ;
444
445 Tenseur xsr_cov(mp, 1, COV, ref_triad) ;
446 xsr_cov = xx_cov / r_bh ;
447 xsr_cov.set_std_base() ;
448
449 Tenseur msr(mp) ;
450 msr = ggrav * mass / r_bh ;
451 msr.set_std_base() ;
452
453 Tenseur lapse_bh(mp) ;
454 lapse_bh = 1. / sqrt( 1.+2.*msr ) ;
455 lapse_bh.set_std_base() ;
456
457 Tenseur lapse_bh2(mp) ; // lapse_bh * lapse_bh
458 lapse_bh2 = 1. / (1.+2.*msr) ;
459 lapse_bh2.set_std_base() ;
460
461 Tenseur ldnu(mp) ;
462 ldnu = flat_scalar_prod_desal(xsr_cov, d_logn_auto) ;
463 ldnu.set_std_base() ;
464
465 Tenseur ldbeta(mp) ;
466 ldbeta = flat_scalar_prod_desal(xsr_cov, d_beta_auto) ;
467 ldbeta.set_std_base() ;
468
469 Tenseur lltkij(mp) ;
470 lltkij.set_etat_qcq() ;
471 lltkij.set() = 0 ;
472 lltkij.set_std_base() ;
473
474 for (int i=0; i<3; i++)
475 for (int j=0; j<3; j++)
476 lltkij.set() += xsr_cov(i) * xsr_cov(j) * tkij_auto(i, j) ;
477
478 Tenseur lshift(mp) ;
479 lshift = flat_scalar_prod_desal(xsr_cov, shift_auto) ;
480 lshift.set_std_base() ;
481
482 Tenseur d_ldnu(mp, 1, COV, ref_triad) ;
483 d_ldnu = ldnu.gradient() ; // (d/dx, d/dy, d/dz)
484 d_ldnu.change_triad(ref_triad) ; // --> (d/dX, d/dY, d/dZ)
485
486 Tenseur ldldnu(mp) ;
487 ldldnu = flat_scalar_prod_desal(xsr_cov, d_ldnu) ;
488 ldldnu.set_std_base() ;
489
490 Tenseur d_ldbeta(mp, 1, COV, ref_triad) ;
491 d_ldbeta = ldbeta.gradient() ; // (d/dx, d/dy, d/dz)
492 d_ldbeta.change_triad(ref_triad) ; // --> (d/dX, d/dY, d/dZ)
493
494 Tenseur ldldbeta(mp) ;
495 ldldbeta = flat_scalar_prod_desal(xsr_cov, d_ldbeta) ;
496 ldldbeta.set_std_base() ;
497
498 //------------------------------------------
499 // Poisson equation for logn_auto (nu_auto)
500 //------------------------------------------
501
502 // Source
503 // ------
504
505 if (relativistic) {
506
507 source = qpig * a_car % (ener_euler + s_euler) + akcar_auto
509 + 2.*lapse_bh2 % msr % (ldnu % ldbeta + ldldnu)
510 + lapse_bh2 % lapse_bh2 % msr % (2.*(ldnu + 4.*msr % ldnu)
511 - ldbeta) / r_bh
512 - (4.*a_car % lapse_bh2 % lapse_bh2 % msr / 3. / nnn / r_bh)
513 * (2.+3.*msr) * (3.+4.*msr) * lltkij
514 + (2.*a_car % lapse_bh2 % lapse_bh2 % lapse_bh % msr
515 / nnn / r_bh / r_bh) * (2.+10.*msr+9.*msr%msr) * lshift
516 + (4.*pow(lapse_bh2, 3.) % msr % msr / 3. / r_bh / r_bh)
517 % (2.*(a_car%lapse_bh2/nnn/nnn - 1.) * pow(2.+3.*msr, 2.)
518 + (a_car - 1.) % pow(1.+3.*msr, 2.)
519 - 3.*(a_car%lapse_bh/nnn - 1.)*(2.+10.*msr+9.*msr%msr)) ;
520
521 }
522 else {
523 cout << "The computation of BH-NS binary systems"
524 << " should be relativistic !!!" << endl ;
525 abort() ;
526 }
527
528 source.set_std_base() ;
529
530 // Resolution of the Poisson equation
531 // ----------------------------------
532
533 int k_falloff = 1 ;
534
535 source().poisson_falloff(par_poisson1, logn_auto.set(), k_falloff) ;
536
537 // Construct logn_auto_regu for et_bin_upmetr_extr.C
538 // -------------------------------------------------
539
541
542 // Check: has the Poisson equation been correctly solved ?
543 // -----------------------------------------------------
544
545 Tbl tdiff_logn = diffrel(logn_auto().laplacien(), source()) ;
546
547 cout <<
548 "Relative error in the resolution of the equation for logn_auto : "
549 << endl ;
550
551 for (int l=0; l<nz; l++) {
552 cout << tdiff_logn(l) << " " ;
553 }
554 cout << endl ;
555 diff_logn = max(abs(tdiff_logn)) ;
556
557 if (relativistic) {
558
559 //--------------------------------
560 // Poisson equation for beta_auto
561 //--------------------------------
562
563 // Source
564 // ------
565
566 source = qpig * a_car % s_euler + 0.75 * akcar_auto
569 + lapse_bh2 % msr % (ldnu%ldnu + ldbeta%ldbeta + 2.*ldldbeta)
570 + lapse_bh2 % lapse_bh2 % msr % (2.*(1.+4.*msr) * ldbeta
571 - ldnu) / r_bh
572 - (a_car % lapse_bh2 % lapse_bh2 % msr / nnn / r_bh)
573 * (2.+3.*msr) * (3.+4.*msr) * lltkij
574 + (2.*a_car % lapse_bh2 % lapse_bh2 % lapse_bh % msr
575 / nnn / r_bh / r_bh) * (2.+10.*msr+9.*msr%msr) * lshift
576 + (2.*pow(lapse_bh2, 3.) % msr % msr / r_bh / r_bh)
577 % ((a_car%lapse_bh2/nnn/nnn - 1.) * pow(2.+3.*msr, 2.)
578 + (a_car - 1.) * pow(1.+3.*msr, 2.)
579 - 2.*(a_car%lapse_bh/nnn - 1.)*(2.+10.*msr+9.*msr%msr)) ;
580
581 source.set_std_base() ;
582
583 // Resolution of the Poisson equation
584 // ----------------------------------
585
586 source().poisson_falloff(par_poisson2, beta_auto.set(),
587 k_falloff) ;
588
589 // Check: has the Poisson equation been correctly solved ?
590 // -----------------------------------------------------
591
592 Tbl tdiff_beta = diffrel(beta_auto().laplacien(), source()) ;
593
594 cout << "Relative error in the resolution of the equation for "
595 << "beta_auto : " << endl ;
596 for (int l=0; l<nz; l++) {
597 cout << tdiff_beta(l) << " " ;
598 }
599 cout << endl ;
600 diff_beta = max(abs(tdiff_beta)) ;
601
602 //----------------------------------------
603 // Vector Poisson equation for shift_auto
604 //----------------------------------------
605
606 // Some auxiliary terms for the source
607 // -----------------------------------
608
609 Tenseur xx_con(mp, 1, CON, ref_triad) ;
610 xx_con.set_etat_qcq() ;
611 xx_con.set(0) = xx + sepa ;
612 xx_con.set(1) = yy ;
613 xx_con.set(2) = zz ;
614 xx_con.set_std_base() ;
615
616 Tenseur xsr_con(mp, 1, CON, ref_triad) ;
617 xsr_con = xx_con / r_bh ;
618 xsr_con.set_std_base() ;
619
620 // Components of shift_auto with respect to the Cartesian triad
621 // (d/dx, d/dy, d/dz) of the mapping :
622
623 Tenseur shift_auto_local = shift_auto ;
624 shift_auto_local.change_triad( mp.get_bvect_cart() ) ;
625
626 // Gradient (partial derivatives with respect to the Cartesian
627 // coordinates of the mapping)
628 // dn(i, j) = D_i N^j
629
630 Tenseur dn = shift_auto_local.gradient() ;
631
632 // Return to the absolute reference frame
634
635 // Trace of D_i N^j = divergence of N^j :
636 Tenseur divn = contract(dn, 0, 1) ;
637
638 // l^j D_j N^i
639 Tenseur ldn_con = contract(xsr_con, 0, dn, 0) ;
640
641 // D_j (l^k D_k N^i): dldn(j, i)
642 Tenseur ldn_local = ldn_con ;
643 ldn_local.change_triad( mp.get_bvect_cart() ) ;
644 Tenseur dldn = ldn_local.gradient() ;
645 dldn.change_triad(ref_triad) ;
646
647 // l^j D_j (l^k D_k N^i)
648 Tenseur ldldn = contract(xsr_con, 0, dldn, 0) ;
649
650 // l_k D_j N^k
651 Tenseur ldn_cov = contract(xsr_cov, 0, dn, 1) ;
652
653 // l^j l_k D_j N^k
654 Tenseur lldn_cov = contract(xsr_con, 0, ldn_cov, 0) ;
655
656 // eta^{ij} l_k D_j N^k
657 Tenseur eldn(mp, 1, CON, ref_triad) ;
658 eldn.set_etat_qcq() ;
659 eldn.set(0) = ldn_cov(0) ;
660 eldn.set(1) = ldn_cov(1) ;
661 eldn.set(2) = ldn_cov(2) ;
662 eldn.set_std_base() ;
663
664 // l^i D_j N^j
665 Tenseur ldivn = xsr_con % divn ;
666
667 // D_j (l^i D_k N^k): dldivn(j, i)
668 Tenseur ldivn_local = ldivn ;
669 ldivn_local.change_triad( mp.get_bvect_cart() ) ;
670 Tenseur dldivn = ldivn_local.gradient() ;
671 dldivn.change_triad(ref_triad) ;
672
673 // l^j D_j (l^i D_k N^k)
674 Tenseur ldldivn = contract(xsr_con, 0, dldivn, 0) ;
675
676 // l_j N^j
677 Tenseur ln = contract(xsr_cov, 0, shift_auto, 0) ;
678
679 Tenseur vtmp = 6. * d_beta_auto - 8. * d_logn_auto ;
680
681 Tenseur lvtmp = contract(xsr_con, 0, vtmp, 0) ;
682
683 // eta^{ij} vtmp_j
684 Tenseur evtmp(mp, 1, CON, ref_triad) ;
685 evtmp.set_etat_qcq() ;
686 evtmp.set(0) = vtmp(0) ;
687 evtmp.set(1) = vtmp(1) ;
688 evtmp.set(2) = vtmp(2) ;
689 evtmp.set_std_base() ;
690
691 // lapse_ns
692 Tenseur lapse_ns(mp) ;
693 lapse_ns = exp(logn_auto) ;
694 lapse_ns.set_std_base() ;
695
696 // Source
697 // ------
698
699 source_shift = (-4.*qpig) * nnn % a_car % (ener_euler + press)
700 % u_euler
702 - 2.*nnn % lapse_bh2 * msr / r_bh
704 + 2.*lapse_bh2 * msr * (3.*ldldn + ldldivn) / 3.
705 - lapse_bh2 * msr / r_bh
706 * (4.*ldivn - lapse_bh2 % (3.*ldn_con + 8.*msr * ldn_con)
707 - (eldn + 2.*lapse_bh2*(9.+11.*msr)*lldn_cov%xsr_con) / 3.)
708 - 2.*lapse_bh2 % lapse_bh2 * msr / r_bh / r_bh
709 * ( (4.+11.*msr) * shift_auto
710 - lapse_bh2 * (12.+51.*msr+46.*msr*msr) * ln % xsr_con )
711 / 3.
712 + 8.*pow(lapse_bh2, 4.) * msr / r_bh / r_bh
713 % (lapse_ns - 1.) * (2.+10.*msr+9.*msr*msr) * xsr_con / 3.
714 + 2.*pow(lapse_bh2, 3.) * msr / r_bh * (2.+3.*msr)
715 * ( (1.+2.*msr) * evtmp - (3.+2.*msr) * lvtmp * xsr_con) / 3. ;
716
717 source_shift.set_std_base() ;
718
719 // Resolution of the Poisson equation
720 // ----------------------------------
721
722 // Filter for the source of shift vector :
723
724 for (int i=0; i<3; i++) {
725 for (int l=0; l<nz; l++) {
726 if (source_shift(i).get_etat() != ETATZERO)
727 source_shift.set(i).filtre_phi(np_filter, l) ;
728 }
729 }
730
731 // For Tenseur::poisson_vect, the triad must be the mapping
732 // triad, not the reference one:
733
734 source_shift.change_triad( mp.get_bvect_cart() ) ;
735 /*
736 for (int i=0; i<3; i++) {
737 if(source_shift(i).dz_nonzero()) {
738 assert( source_shift(i).get_dzpuis() == 4 ) ;
739 }
740 else {
741 (source_shift.set(i)).set_dzpuis(4) ;
742 }
743 }
744
745 source_shift.dec2_dzpuis() ; // dzpuis 4 -> 2
746 */
747 double lambda_shift = double(1) / double(3) ;
748
749 int* shift_falloff ;
750 shift_falloff = new int[4] ;
751 shift_falloff[0] = 1 ;
752 shift_falloff[1] = 1 ;
753 shift_falloff[2] = 2 ;
754 shift_falloff[3] = 1 ;
755
756 source_shift.poisson_vect_falloff(lambda_shift, par_poisson_vect,
758 khi_shift, shift_falloff) ;
759
760 delete[] shift_falloff ;
761
762 // Check: has the equation for shift_auto been correctly solved ?
763 // --------------------------------------------------------------
764
765 // Divergence of shift_auto :
766 Tenseur divna = contract(shift_auto.gradient(), 0, 1) ;
767 // divna.dec2_dzpuis() ; // dzpuis 2 -> 0
768
769 // Grad(div) :
770 Tenseur graddivn = divna.gradient() ;
771 // graddivn.inc2_dzpuis() ; // dzpuis 2 -> 4
772
773 // Full operator :
774 Tenseur lap_shift(mp, 1, CON, mp.get_bvect_cart() ) ;
775 lap_shift.set_etat_qcq() ;
776 for (int i=0; i<3; i++) {
777 lap_shift.set(i) = shift_auto(i).laplacien()
778 + lambda_shift * graddivn(i) ;
779 }
780
781 Tbl tdiff_shift_x = diffrel(lap_shift(0), source_shift(0)) ;
782 Tbl tdiff_shift_y = diffrel(lap_shift(1), source_shift(1)) ;
783 Tbl tdiff_shift_z = diffrel(lap_shift(2), source_shift(2)) ;
784
785 cout << "Relative error in the resolution of the equation for "
786 << "shift_auto : " << endl ;
787 cout << "x component : " ;
788 for (int l=0; l<nz; l++) {
789 cout << tdiff_shift_x(l) << " " ;
790 }
791 cout << endl ;
792 cout << "y component : " ;
793 for (int l=0; l<nz; l++) {
794 cout << tdiff_shift_y(l) << " " ;
795 }
796 cout << endl ;
797 cout << "z component : " ;
798 for (int l=0; l<nz; l++) {
799 cout << tdiff_shift_z(l) << " " ;
800 }
801 cout << endl ;
802
803 diff_shift_x = max(abs(tdiff_shift_x)) ;
804 diff_shift_y = max(abs(tdiff_shift_y)) ;
805 diff_shift_z = max(abs(tdiff_shift_z)) ;
806
807 // Final result
808 // ------------
809 // The output of Tenseur::poisson_vect_falloff is on the mapping
810 // triad, it should therefore be transformed to components on the
811 // reference triad :
812
814
815 } // End of relativistic equations
816
817 //------------------------------
818 // Relative change in enthalpy
819 //------------------------------
820
821 Tbl diff_ent_tbl = diffrel( ent(), ent_jm1() ) ;
822 diff_ent = diff_ent_tbl(0) ;
823 for (int l=1; l<nzet; l++) {
824 diff_ent += diff_ent_tbl(l) ;
825 }
826 diff_ent /= nzet ;
827
828 ent_jm1 = ent ;
829
830 } // End of main loop
831
832 //========================================================//
833 // End of iteration //
834 //========================================================//
835
836}
837
838
839void Et_bin_bhns_extr::equil_bhns_extr_cf(double ent_c, const double& mass,
840 const double& sepa, int mermax,
841 int mermax_poisson,
842 double relax_poisson,
843 int mermax_potvit,
844 double relax_potvit, int np_filter,
845 double thres_adapt,
846 Tbl& diff) {
847
848 // Fundamental constants and units
849 // -------------------------------
850 using namespace Unites ;
851
852 assert( kerrschild == false ) ;
853
854 // Initializations
855 // ---------------
856
857 const Mg3d* mg = mp.get_mg() ;
858 int nz = mg->get_nzone() ; // total number of domains
859
860 // The following is required to initialize mp_prev as a Map_et:
861 Map_et& mp_et = dynamic_cast<Map_et&>(mp) ;
862
863 // Domain and radial indices of points at the surface of the star:
864 int l_b = nzet - 1 ;
865 int i_b = mg->get_nr(l_b) - 1 ;
866 int k_b ;
867 int j_b ;
868
869 // Value of the enthalpy defining the surface of the star
870 double ent_b = 0 ;
871
872 // Error indicators
873 // ----------------
874
875 double& diff_ent = diff.set(0) ;
876 double& diff_vel_pot = diff.set(1) ;
877 double& diff_logn = diff.set(2) ;
878 double& diff_beta = diff.set(3) ;
879 double& diff_shift_x = diff.set(4) ;
880 double& diff_shift_y = diff.set(5) ;
881 double& diff_shift_z = diff.set(6) ;
882
883 // Parameters for the function Map_et::adapt
884 // -----------------------------------------
885
886 Param par_adapt ;
887 int nitermax = 100 ;
888 int niter ;
889 int adapt_flag = 1 ; // 1 = performs the full computation,
890 // 0 = performs only the rescaling by
891 // the factor alpha_r
892 int nz_search = nzet ; // Number of domains for searching
893 // the enthalpy isosurfaces
894 double precis_secant = 1.e-14 ;
895 double alpha_r ;
896 double reg_map = 1. ; // 1 = regular mapping, 0 = contracting mapping
897
898 Tbl ent_limit(nz) ;
899
900 par_adapt.add_int(nitermax, 0) ; // maximum number of iterations to
901 // locate zeros by the secant method
902 par_adapt.add_int(nzet, 1) ; // number of domains where the adjustment
903 // to the isosurfaces of ent is to be
904 // performed
905 par_adapt.add_int(nz_search, 2) ; // number of domains to search for
906 // the enthalpy isosurface
907 par_adapt.add_int(adapt_flag, 3) ; // 1 = performs the full computation,
908 // 0 = performs only the rescaling by
909 // the factor alpha_r
910 par_adapt.add_int(j_b, 4) ; // theta index of the collocation point
911 // (theta_*, phi_*)
912 par_adapt.add_int(k_b, 5) ; // theta index of the collocation point
913 // (theta_*, phi_*)
914 par_adapt.add_int_mod(niter, 0) ; // number of iterations actually
915 // used in the secant method
916 par_adapt.add_double(precis_secant, 0) ; // required absolute precision in
917 // the determination of zeros by
918 // the secant method
919 par_adapt.add_double(reg_map, 1) ; // 1 = regular mapping,
920 // 0 = contracting mapping
921 par_adapt.add_double(alpha_r, 2) ; // factor by which all the radial
922 // distances will be multiplied
923 par_adapt.add_tbl(ent_limit, 0) ; // array of values of the field ent
924 // to define the isosurfaces
925
926 // Parameters for the function Map_et::poisson for logn_auto
927 // ---------------------------------------------------------
928
929 double precis_poisson = 1.e-16 ;
930
931 Param par_poisson1 ;
932
933 par_poisson1.add_int(mermax_poisson, 0) ; // maximum number of iterations
934 par_poisson1.add_double(relax_poisson, 0) ; // relaxation parameter
935 par_poisson1.add_double(precis_poisson, 1) ; // required precision
936 par_poisson1.add_int_mod(niter, 0) ; // number of iterations actually
937 // used
938 par_poisson1.add_cmp_mod( ssjm1_logn ) ;
939
940 // Parameters for the function Map_et::poisson for beta_auto
941 // ---------------------------------------------------------
942
943 Param par_poisson2 ;
944
945 par_poisson2.add_int(mermax_poisson, 0) ; // maximum number of iterations
946 par_poisson2.add_double(relax_poisson, 0) ; // relaxation parameter
947 par_poisson2.add_double(precis_poisson, 1) ; // required precision
948 par_poisson2.add_int_mod(niter, 0) ; // number of iterations actually
949 // used
950 par_poisson2.add_cmp_mod( ssjm1_beta ) ;
951
952 // Parameters for the function Tenseur::poisson_vect
953 // -------------------------------------------------
954
955 Param par_poisson_vect ;
956
957 par_poisson_vect.add_int(mermax_poisson, 0) ; // maximum number of
958 // iterations
959 par_poisson_vect.add_double(relax_poisson, 0) ; // relaxation parameter
960 par_poisson_vect.add_double(precis_poisson, 1) ; // required precision
961 par_poisson_vect.add_cmp_mod( ssjm1_khi ) ;
962 par_poisson_vect.add_tenseur_mod( ssjm1_wshift ) ;
963 par_poisson_vect.add_int_mod(niter, 0) ;
964
965 // External potential
966 // See Eq (99) from Gourgoulhon et al. (2001)
967 // -----------------------------------------
968
969 Tenseur pot_ext = logn_comp + pot_centri + loggam ;
970
971 Tenseur ent_jm1 = ent ; // Enthalpy at previous step
972
973 Tenseur source(mp) ; // source term in the equation for logn_auto
974 // and beta_auto
975
976 Tenseur source_shift(mp, 1, CON, ref_triad) ; // source term in the
977 // equation for shift_auto
978
979 //==========================================================//
980 // Start of iteration //
981 //==========================================================//
982
983 for(int mer=0 ; mer<mermax ; mer++ ) {
984
985 cout << "-----------------------------------------------" << endl ;
986 cout << "step: " << mer << endl ;
987 cout << "diff_ent = " << diff_ent << endl ;
988
989 //------------------------------------------------------
990 // Resolution of the elliptic equation for the velocity
991 // scalar potential
992 //------------------------------------------------------
993
994 if (irrotational) {
995 diff_vel_pot = velocity_pot_extr(mass, sepa, mermax_potvit,
996 precis_poisson, relax_potvit) ;
997 }
998
999 //-------------------------------------
1000 // Computation of the new radial scale
1001 //-------------------------------------
1002
1003 // alpha_r (r = alpha_r r') is determined so that the enthalpy
1004 // takes the requested value ent_b at the stellar surface
1005
1006 // Values at the center of the star:
1007 double logn_auto_c = logn_auto()(0, 0, 0, 0) ;
1008 double pot_ext_c = pot_ext()(0, 0, 0, 0) ;
1009
1010 // Search for the reference point (theta_*, phi_*)
1011 // [notation of Bonazzola, Gourgoulhon & Marck PRD 58, 104020 (1998)]
1012 // at the surface of the star
1013 // ------------------------------------------------------------------
1014 double alpha_r2 = 0 ;
1015 for (int k=0; k<mg->get_np(l_b); k++) {
1016 for (int j=0; j<mg->get_nt(l_b); j++) {
1017
1018 double pot_ext_b = pot_ext()(l_b, k, j, i_b) ;
1019 double logn_auto_b = logn_auto()(l_b, k, j, i_b) ;
1020
1021 // See Eq (100) from Gourgoulhon et al. (2001)
1022 double alpha_r2_jk = ( ent_c - ent_b + pot_ext_c - pot_ext_b)
1023 / ( logn_auto_b - logn_auto_c ) ;
1024
1025 if (alpha_r2_jk > alpha_r2) {
1026 alpha_r2 = alpha_r2_jk ;
1027 k_b = k ;
1028 j_b = j ;
1029 }
1030
1031 }
1032 }
1033
1034 alpha_r = sqrt(alpha_r2) ;
1035
1036 cout << "k_b, j_b, alpha_r: " << k_b << " " << j_b << " "
1037 << alpha_r << endl ;
1038
1039 // New value of logn_auto
1040 // ----------------------
1041
1042 logn_auto = alpha_r2 * logn_auto ;
1043 logn_auto_regu = alpha_r2 * logn_auto_regu ;
1044 logn_auto_c = logn_auto()(0, 0, 0, 0) ;
1045
1046 //------------------------------------------------------------
1047 // Change the values of the inner points of the second domain
1048 // by those of the outer points of the first domain
1049 //------------------------------------------------------------
1050
1051 (logn_auto().va).smooth(nzet, (logn_auto.set()).va) ;
1052
1053 //--------------------------------------------
1054 // First integral --> enthalpy in all space
1055 // See Eq (98) from Gourgoulhon et al. (2001)
1056 //--------------------------------------------
1057
1058 ent = (ent_c + logn_auto_c + pot_ext_c) - logn_auto - pot_ext ;
1059
1060 //---------------------------------------------------------
1061 // Change the enthalpy field to accelerate the convergence
1062 //---------------------------------------------------------
1063
1064 double dentdx = ent().dsdx()(0, 0, 0, 0) ;
1065 double dentdy = ent().dsdy()(0, 0, 0, 0) ;
1066
1067 cout << "dH/dx|_center = " << dentdx << endl ;
1068 cout << "dH/dy|_center = " << dentdy << endl ;
1069
1070 double dec_fact = 1. ;
1071
1072 Tenseur func_in(mp) ;
1073 func_in.set_etat_qcq() ;
1074 func_in.set() = 1. - dec_fact * (dentdx/ent_c) * mp.x ;
1075 func_in.set().annule(nzet, nz-1) ;
1076 func_in.set_std_base() ;
1077
1078 Tenseur func_ex(mp) ;
1079 func_ex.set_etat_qcq() ;
1080 func_ex.set() = 1. ;
1081 func_ex.set().annule(0, nzet-1) ;
1082 func_ex.set_std_base() ;
1083
1084 // New enthalpy field
1085 // ------------------
1086 ent.set() = ent() * (func_in() + func_ex()) ;
1087
1088 (ent().va).smooth(nzet, (ent.set()).va) ;
1089
1090 double dentdx_new = ent().dsdx()(0, 0, 0, 0) ;
1091
1092 cout << "dH/dx|_new = " << dentdx_new << endl ;
1093
1094 //-----------------------------------------------------
1095 // Adaptation of the mapping to the new enthalpy field
1096 //----------------------------------------------------
1097
1098 // Shall the adaptation be performed (cusp) ?
1099 // ------------------------------------------
1100
1101 double dent_eq = ent().dsdr().val_point(ray_eq_pi(),M_PI/2.,M_PI) ;
1102 double dent_pole = ent().dsdr().val_point(ray_pole(),0.,0.) ;
1103 double rap_dent = fabs( dent_eq / dent_pole ) ;
1104 cout << "| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
1105
1106 if ( rap_dent < thres_adapt ) {
1107 adapt_flag = 0 ; // No adaptation of the mapping
1108 cout << "******* FROZEN MAPPING *********" << endl ;
1109 }
1110 else{
1111 adapt_flag = 1 ; // The adaptation of the mapping is to be
1112 // performed
1113 }
1114
1115 ent_limit.set_etat_qcq() ;
1116 for (int l=0; l<nzet; l++) { // loop on domains inside the star
1117 ent_limit.set(l) = ent()(l, k_b, j_b, i_b) ;
1118 }
1119
1120 ent_limit.set(nzet-1) = ent_b ;
1121 Map_et mp_prev = mp_et ;
1122
1123 mp.adapt(ent(), par_adapt) ;
1124
1125 //----------------------------------------------------
1126 // Computation of the enthalpy at the new grid points
1127 //----------------------------------------------------
1128
1129 mp_prev.homothetie(alpha_r) ;
1130
1131 mp.reevaluate_symy(&mp_prev, nzet+1, ent.set()) ;
1132
1133 double fact_resize = 1. / alpha_r ;
1134 for (int l=nzet; l<nz-1; l++) {
1135 mp_et.resize(l, fact_resize) ;
1136 }
1137 mp_et.resize_extr(fact_resize) ;
1138
1139 double ent_s_max = -1 ;
1140 int k_s_max = -1 ;
1141 int j_s_max = -1 ;
1142 for (int k=0; k<mg->get_np(l_b); k++) {
1143 for (int j=0; j<mg->get_nt(l_b); j++) {
1144 double xx = fabs( ent()(l_b, k, j, i_b) ) ;
1145 if (xx > ent_s_max) {
1146 ent_s_max = xx ;
1147 k_s_max = k ;
1148 j_s_max = j ;
1149 }
1150 }
1151 }
1152 cout << "Max. abs(enthalpy) at the boundary between domains nzet-1"
1153 << " and nzet : " << endl ;
1154 cout << " " << ent_s_max << " reached for k = " << k_s_max
1155 << " and j = " << j_s_max << endl ;
1156
1157 //-------------------
1158 // Equation of state
1159 //-------------------
1160
1161 equation_of_state() ; // computes new values for nbar (n), ener (e)
1162 // and press (p) from the new ent (H)
1163
1164 //----------------------------------------------------------
1165 // Matter source terms in the gravitational field equations
1166 //---------------------------------------------------------
1167
1168 hydro_euler_extr(mass, sepa) ; // computes new values for
1169 // ener_euler (E), s_euler (S)
1170 // and u_euler (U^i)
1171
1172
1173 //----------------------------------------------------
1174 // Auxiliary terms for the source of Poisson equation
1175 //----------------------------------------------------
1176
1177 const Coord& xx = mp.x ;
1178 const Coord& yy = mp.y ;
1179 const Coord& zz = mp.z ;
1180
1181 Tenseur r_bh(mp) ;
1182 r_bh.set_etat_qcq() ;
1183 r_bh.set() = pow( (xx+sepa)*(xx+sepa) + yy*yy + zz*zz, 0.5) ;
1184 r_bh.set_std_base() ;
1185
1186 Tenseur xx_cov(mp, 1, COV, ref_triad) ;
1187 xx_cov.set_etat_qcq() ;
1188 xx_cov.set(0) = xx + sepa ;
1189 xx_cov.set(1) = yy ;
1190 xx_cov.set(2) = zz ;
1191 xx_cov.set_std_base() ;
1192
1193 Tenseur msr(mp) ;
1194 msr = ggrav * mass / r_bh ;
1195 msr.set_std_base() ;
1196
1197 Tenseur tmp(mp) ;
1198 tmp = 1. / ( 1. - 0.25*msr*msr ) ;
1199 tmp.set_std_base() ;
1200
1201 Tenseur xdnu(mp) ;
1202 xdnu = flat_scalar_prod_desal(xx_cov, d_logn_auto) ;
1203 xdnu.set_std_base() ;
1204
1205 Tenseur xdbeta(mp) ;
1206 xdbeta = flat_scalar_prod_desal(xx_cov, d_beta_auto) ;
1207 xdbeta.set_std_base() ;
1208
1209 //------------------------------------------
1210 // Poisson equation for logn_auto (nu_auto)
1211 //------------------------------------------
1212
1213 // Source
1214 // ------
1215
1216 if (relativistic) {
1217
1218 source = qpig * a_car % (ener_euler + s_euler) + akcar_auto
1220 - 0.5 * tmp % msr % msr % xdnu / r_bh / r_bh
1221 - tmp % msr % xdbeta / r_bh / r_bh ;
1222
1223 }
1224 else {
1225 cout << "The computation of BH-NS binary systems"
1226 << " should be relativistic !!!" << endl ;
1227 abort() ;
1228 }
1229
1230 source.set_std_base() ;
1231
1232 // Resolution of the Poisson equation
1233 // ----------------------------------
1234
1235 int k_falloff = 1 ;
1236
1237 source().poisson_falloff(par_poisson1, logn_auto.set(), k_falloff) ;
1238
1239 // Construct logn_auto_regu for et_bin_upmetr_extr.C
1240 // -------------------------------------------------
1241
1243
1244 // Check: has the Poisson equation been correctly solved ?
1245 // -----------------------------------------------------
1246
1247 Tbl tdiff_logn = diffrel(logn_auto().laplacien(), source()) ;
1248
1249 cout <<
1250 "Relative error in the resolution of the equation for logn_auto : "
1251 << endl ;
1252
1253 for (int l=0; l<nz; l++) {
1254 cout << tdiff_logn(l) << " " ;
1255 }
1256 cout << endl ;
1257 diff_logn = max(abs(tdiff_logn)) ;
1258
1259 if (relativistic) {
1260
1261 //--------------------------------
1262 // Poisson equation for beta_auto
1263 //--------------------------------
1264
1265 // Source
1266 // ------
1267
1268 source = qpig * a_car % s_euler + 0.75 * akcar_auto
1271 - tmp % msr % xdnu / r_bh / r_bh
1272 - 0.5 * tmp % msr %msr % xdbeta / r_bh / r_bh ;
1273
1274 source.set_std_base() ;
1275
1276 // Resolution of the Poisson equation
1277 // ----------------------------------
1278
1279 source().poisson_falloff(par_poisson2, beta_auto.set(),
1280 k_falloff) ;
1281
1282 // Check: has the Poisson equation been correctly solved ?
1283 // -----------------------------------------------------
1284
1285 Tbl tdiff_beta = diffrel(beta_auto().laplacien(), source()) ;
1286
1287 cout << "Relative error in the resolution of the equation for "
1288 << "beta_auto : " << endl ;
1289 for (int l=0; l<nz; l++) {
1290 cout << tdiff_beta(l) << " " ;
1291 }
1292 cout << endl ;
1293 diff_beta = max(abs(tdiff_beta)) ;
1294
1295 //----------------------------------------
1296 // Vector Poisson equation for shift_auto
1297 //----------------------------------------
1298
1299 // Some auxiliary terms for the source
1300 // -----------------------------------
1301
1302 Tenseur bhtmp(mp, 1, COV, ref_triad) ;
1303 bhtmp.set_etat_qcq() ;
1304 bhtmp = tmp % msr % (3.*msr-8.) % xx_cov / r_bh / r_bh ;
1305 bhtmp.set_std_base() ;
1306
1307 Tenseur vtmp = 6. * d_beta_auto - 8. * d_logn_auto ;
1308
1309 // Source
1310 // ------
1311
1312 source_shift = (-4.*qpig) * nnn % a_car % (ener_euler + press)
1313 % u_euler
1314 + nnn % flat_scalar_prod_desal(tkij_auto, vtmp+bhtmp) ;
1315
1316 source_shift.set_std_base() ;
1317
1318 // Resolution of the Poisson equation
1319 // ----------------------------------
1320
1321 // Filter for the source of shift vector :
1322
1323 for (int i=0; i<3; i++) {
1324 for (int l=0; l<nz; l++) {
1325 if (source_shift(i).get_etat() != ETATZERO)
1326 source_shift.set(i).filtre_phi(np_filter, l) ;
1327 }
1328 }
1329
1330 // For Tenseur::poisson_vect, the triad must be the mapping
1331 // triad, not the reference one:
1332
1333 source_shift.change_triad( mp.get_bvect_cart() ) ;
1334 /*
1335 for (int i=0; i<3; i++) {
1336 if(source_shift(i).dz_nonzero()) {
1337 assert( source_shift(i).get_dzpuis() == 4 ) ;
1338 }
1339 else {
1340 (source_shift.set(i)).set_dzpuis(4) ;
1341 }
1342 }
1343
1344 source_shift.dec2_dzpuis() ; // dzpuis 4 -> 2
1345 */
1346 double lambda_shift = double(1) / double(3) ;
1347
1348 int* shift_falloff ;
1349 shift_falloff = new int[4] ;
1350 shift_falloff[0] = 1 ;
1351 shift_falloff[1] = 1 ;
1352 shift_falloff[2] = 2 ;
1353 shift_falloff[3] = 1 ;
1354
1355 source_shift.poisson_vect_falloff(lambda_shift, par_poisson_vect,
1357 khi_shift, shift_falloff) ;
1358
1359 delete[] shift_falloff ;
1360
1361 // Check: has the equation for shift_auto been correctly solved ?
1362 // --------------------------------------------------------------
1363
1364 // Divergence of shift_auto :
1365 Tenseur divna = contract(shift_auto.gradient(), 0, 1) ;
1366 // divna.dec2_dzpuis() ; // dzpuis 2 -> 0
1367
1368 // Grad(div) :
1369 Tenseur graddivn = divna.gradient() ;
1370 // graddivn.inc2_dzpuis() ; // dzpuis 2 -> 4
1371
1372 // Full operator :
1373 Tenseur lap_shift(mp, 1, CON, mp.get_bvect_cart() ) ;
1374 lap_shift.set_etat_qcq() ;
1375 for (int i=0; i<3; i++) {
1376 lap_shift.set(i) = shift_auto(i).laplacien()
1377 + lambda_shift * graddivn(i) ;
1378 }
1379
1380 Tbl tdiff_shift_x = diffrel(lap_shift(0), source_shift(0)) ;
1381 Tbl tdiff_shift_y = diffrel(lap_shift(1), source_shift(1)) ;
1382 Tbl tdiff_shift_z = diffrel(lap_shift(2), source_shift(2)) ;
1383
1384 cout << "Relative error in the resolution of the equation for "
1385 << "shift_auto : " << endl ;
1386 cout << "x component : " ;
1387 for (int l=0; l<nz; l++) {
1388 cout << tdiff_shift_x(l) << " " ;
1389 }
1390 cout << endl ;
1391 cout << "y component : " ;
1392 for (int l=0; l<nz; l++) {
1393 cout << tdiff_shift_y(l) << " " ;
1394 }
1395 cout << endl ;
1396 cout << "z component : " ;
1397 for (int l=0; l<nz; l++) {
1398 cout << tdiff_shift_z(l) << " " ;
1399 }
1400 cout << endl ;
1401
1402 diff_shift_x = max(abs(tdiff_shift_x)) ;
1403 diff_shift_y = max(abs(tdiff_shift_y)) ;
1404 diff_shift_z = max(abs(tdiff_shift_z)) ;
1405
1406 // Final result
1407 // ------------
1408 // The output of Tenseur::poisson_vect_falloff is on the mapping
1409 // triad, it should therefore be transformed to components on the
1410 // reference triad :
1411
1413
1414 } // End of relativistic equations
1415
1416 //------------------------------
1417 // Relative change in enthalpy
1418 //------------------------------
1419
1420 Tbl diff_ent_tbl = diffrel( ent(), ent_jm1() ) ;
1421 diff_ent = diff_ent_tbl(0) ;
1422 for (int l=1; l<nzet; l++) {
1423 diff_ent += diff_ent_tbl(l) ;
1424 }
1425 diff_ent /= nzet ;
1426
1427 ent_jm1 = ent ;
1428
1429 } // End of main loop
1430
1431 //========================================================//
1432 // End of iteration //
1433 //========================================================//
1434
1435}
1436
1437}
void filtre_phi(int n, int zone)
Sets the n lasts coefficients in to 0 in the domain zone .
Definition cmp_manip.C:101
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:348
Active physical coordinates and mapping derivatives.
Definition coord.h:90
void equil_bhns_extr_ks(double ent_c, const double &mass, const double &sepa, int mermax, int mermax_poisson, double relax_poisson, int mermax_potvit, double relax_potvit, int np_filter, double thres_adapt, Tbl &diff)
Computes an equilibrium configuration of a BH-NS binary system in the Kerr-Schild background metric u...
bool kerrschild
Indicator of the background metric: true for the Kerr-Shild metric, false for the conformally flat on...
double velocity_pot_extr(const double &mass, const double &sepa, int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
void equil_bhns_extr_cf(double ent_c, const double &mass, const double &sepa, int mermax, int mermax_poisson, double relax_poisson, int mermax_potvit, double relax_potvit, int np_filter, double thres_adapt, Tbl &diff)
Computes an equilibrium configuration of a BH-NS binary system in the conformally flat background met...
void hydro_euler_extr(const double &mass, const double &sepa)
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Tenseur d_beta_auto
Gradient of beta_auto (Cartesian components with respect to ref_triad )
Definition etoile.h:879
Tenseur d_logn_auto
Gradient of logn_auto (Cartesian components with respect to ref_triad )
Definition etoile.h:859
Cmp ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Definition etoile.h:973
const Base_vect & ref_triad
Reference triad ("absolute frame"), with respect to which the components of all the member Tenseur 's...
Definition etoile.h:828
bool irrotational
true for an irrotational star, false for a corotating one
Definition etoile.h:822
Tenseur w_shift
Vector used in the decomposition of shift_auto , following Shibata's prescription [Prog.
Definition etoile.h:908
Tenseur logn_comp
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principaly by ...
Definition etoile.h:854
Tenseur ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for by means...
Definition etoile.h:983
Cmp ssjm1_logn
Effective source at the previous step for the resolution of the Poisson equation for logn_auto by mea...
Definition etoile.h:959
Tenseur pot_centri
Centrifugal potential.
Definition etoile.h:953
Tenseur loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition etoile.h:849
Tenseur shift_auto
Part of the shift vector generated principaly by the star.
Definition etoile.h:889
Tenseur_sym tkij_auto
Part of the extrinsic curvature tensor generated by shift_auto .
Definition etoile.h:925
Tenseur akcar_auto
Part of the scalar generated by shift_auto , i.e.
Definition etoile.h:938
Cmp ssjm1_beta
Effective source at the previous step for the resolution of the Poisson equation for beta_auto by mea...
Definition etoile.h:965
Tenseur khi_shift
Scalar used in the decomposition of shift_auto , following Shibata's prescription [Prog.
Definition etoile.h:918
int nzet
Number of domains of *mp occupied by the star.
Definition etoile.h:432
Tenseur logn_auto_regu
Regular part of the logarithm of the part of the lapse N generated principaly by the star.
Definition etoile.h:491
double ray_eq_pi() const
Coordinate radius at , [r_unit].
Tenseur nnn
Total lapse function.
Definition etoile.h:509
Tenseur logn_auto
Total of the logarithm of the part of the lapse N generated principaly by the star.
Definition etoile.h:484
virtual void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition etoile.C:566
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition etoile.h:474
Map & mp
Mapping associated with the star.
Definition etoile.h:429
Tenseur press
Fluid pressure.
Definition etoile.h:461
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one.
Definition etoile.h:437
Tenseur ener_euler
Total energy density in the Eulerian frame.
Definition etoile.h:465
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
Definition etoile.h:468
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case)
Definition etoile.h:457
Tenseur beta_auto
Logarithm of the part of the product AN generated principaly by by the star.
Definition etoile.h:506
Tenseur a_car
Total conformal factor .
Definition etoile.h:515
double ray_pole() const
Coordinate radius at [r_unit].
Radial mapping of rather general form.
Definition map.h:2752
void resize_extr(double lambda)
Rescales the outer boundary of the outermost domain in the case of non-compactified external domain.
virtual void homothetie(double lambda)
Sets a new radial scale.
Definition map_et.C:905
virtual void resize(int l, double lambda)
Rescales the outer boundary of one domain.
Definition map_et.C:928
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 reevaluate_symy(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 adapt(const Cmp &ent, const Param &par, int nbr=0)=0
Adaptation of the mapping to a given scalar field.
Coord z
z coordinate centered on the grid
Definition map.h:728
Coord x
x coordinate centered on the grid
Definition map.h:726
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
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition param.C:315
void add_cmp_mod(Cmp &ti, int position=0)
Adds the address of a new modifiable Cmp to the list.
Definition param.C:1004
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition param.C:385
void add_tenseur_mod(Tenseur &ti, int position=0)
Adds the address of a new modifiable Tenseur to the list.
Definition param.C:1142
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition param.C:246
void add_tbl(const Tbl &ti, int position=0)
Adds the address of a new Tbl to the list.
Definition param.C:522
Basic array class.
Definition tbl.h:161
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tenseur.C:636
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition tenseur.C:1542
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition tenseur.C:1170
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition tenseur.C:668
int get_etat() const
Returns the logical state.
Definition tenseur.h:707
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
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Tenseur flat_scalar_prod_desal(const Tenseur &t1, const Tenseur &t2)
Same as flat_scalar_prod but with desaliasing.
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.