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