LORENE
et_bin_bhns_extr_eq_ylm.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 * 2004 Joshua A. Faber
12 *
13 * This file is part of LORENE.
14 *
15 * LORENE is free software; you can redistribute it and/or modify
16 * it under the terms of the GNU General Public License version 2
17 * as published by the Free Software Foundation.
18 *
19 * LORENE is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License
25 * along with LORENE; if not, write to the Free Software
26 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27 *
28 */
29
30
31
32/*
33 * $Id: et_bin_bhns_extr_eq_ylm.C,v 1.8 2016/12/05 16:17:52 j_novak Exp $
34 * $Log: et_bin_bhns_extr_eq_ylm.C,v $
35 * Revision 1.8 2016/12/05 16:17:52 j_novak
36 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
37 *
38 * Revision 1.7 2014/10/13 08:52:54 j_novak
39 * Lorene classes and functions now belong to the namespace Lorene.
40 *
41 * Revision 1.6 2014/10/06 15:13:07 j_novak
42 * Modified #include directives to use c++ syntax.
43 *
44 * Revision 1.5 2005/02/28 23:12:28 k_taniguchi
45 * Modification to include the case of the conformally flat background metric
46 *
47 * Revision 1.4 2005/01/31 20:28:14 k_taniguchi
48 * Addition of a number of coefficients which are deleted in the filtre_phi
49 * in the argument.
50 *
51 * Revision 1.3 2005/01/25 17:42:02 k_taniguchi
52 * Suppression of the filter for the source term of the shift vector.
53 *
54 * Revision 1.2 2005/01/03 18:03:39 k_taniguchi
55 * Addition of the method to fix the position of the neutron star
56 * in the coordinate system.
57 *
58 * Revision 1.1 2004/12/29 16:31:24 k_taniguchi
59 * *** empty log message ***
60 *
61 *
62 *
63 * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_bhns_extr_eq_ylm.C,v 1.8 2016/12/05 16:17:52 j_novak Exp $
64 *
65 */
66
67// C headers
68#include <cmath>
69
70// Lorene headers
71#include "et_bin_bhns_extr.h"
72#include "etoile.h"
73#include "map.h"
74#include "coord.h"
75#include "param.h"
76#include "eos.h"
77#include "graphique.h"
78#include "utilitaires.h"
79#include "unites.h"
80
81namespace Lorene {
82void Et_bin_bhns_extr::equil_bhns_extr_ylm_ks(double ent_c,
83 const double& mass,
84 const double& sepa,
85 double* nu_int,
86 double* beta_int,
87 double* shift_int,
88 int mermax,
89 int mermax_poisson,
90 double relax_poisson,
91 double relax_ylm,
92 int mermax_potvit,
93 double relax_potvit,
94 int np_filter,
95 double thres_adapt,
96 Tbl& diff) {
97
98 // Fundamental constants and units
99 // -------------------------------
100 using namespace Unites ;
101
102 assert( kerrschild == true ) ;
103
104 // Initializations
105 // ---------------
106
107 const Mg3d* mg = mp.get_mg() ;
108 int nz = mg->get_nzone() ; // total number of domains
109
110 // The following is required to initialize mp_prev as a Map_et:
111 Map_et& mp_et = dynamic_cast<Map_et&>(mp) ;
112
113 // Domain and radial indices of points at the surface of the star:
114 int l_b = nzet - 1 ;
115 int i_b = mg->get_nr(l_b) - 1 ;
116 int k_b ;
117 int j_b ;
118
119 // Value of the enthalpy defining the surface of the star
120 double ent_b = 0 ;
121
122 // Error indicators
123 // ----------------
124
125 double& diff_ent = diff.set(0) ;
126 double& diff_vel_pot = diff.set(1) ;
127 double& diff_logn = diff.set(2) ;
128 double& diff_beta = diff.set(3) ;
129 double& diff_shift_x = diff.set(4) ;
130 double& diff_shift_y = diff.set(5) ;
131 double& diff_shift_z = diff.set(6) ;
132
133 // Parameters for the function Map_et::adapt
134 // -----------------------------------------
135
136 Param par_adapt ;
137 int nitermax = 100 ;
138 int niter ;
139 int adapt_flag = 1 ; // 1 = performs the full computation,
140 // 0 = performs only the rescaling by
141 // the factor alpha_r
142 int nz_search = nzet ; // Number of domains for searching
143 // the enthalpy isosurfaces
144 double precis_secant = 1.e-14 ;
145 double alpha_r ;
146 double reg_map = 1. ; // 1 = regular mapping, 0 = contracting mapping
147
148 Tbl ent_limit(nz) ;
149
150 par_adapt.add_int(nitermax, 0) ; // maximum number of iterations to
151 // locate zeros by the secant method
152 par_adapt.add_int(nzet, 1) ; // number of domains where the adjustment
153 // to the isosurfaces of ent is to be
154 // performed
155 par_adapt.add_int(nz_search, 2) ; // number of domains to search for
156 // the enthalpy isosurface
157 par_adapt.add_int(adapt_flag, 3) ; // 1 = performs the full computation,
158 // 0 = performs only the rescaling by
159 // the factor alpha_r
160 par_adapt.add_int(j_b, 4) ; // theta index of the collocation point
161 // (theta_*, phi_*)
162 par_adapt.add_int(k_b, 5) ; // theta index of the collocation point
163 // (theta_*, phi_*)
164 par_adapt.add_int_mod(niter, 0) ; // number of iterations actually
165 // used in the secant method
166 par_adapt.add_double(precis_secant, 0) ; // required absolute precision in
167 // the determination of zeros by
168 // the secant method
169 par_adapt.add_double(reg_map, 1) ; // 1 = regular mapping,
170 // 0 = contracting mapping
171 par_adapt.add_double(alpha_r, 2) ; // factor by which all the radial
172 // distances will be multiplied
173 par_adapt.add_tbl(ent_limit, 0) ; // array of values of the field ent
174 // to define the isosurfaces
175
176 // Parameters for the function Map_et::poisson for logn_auto
177 // ---------------------------------------------------------
178
179 double precis_poisson = 1.e-16 ;
180
181 Param par_poisson1 ;
182
183 par_poisson1.add_int(mermax_poisson, 0) ; // maximum number of iterations
184 par_poisson1.add_double(relax_poisson, 0) ; // relaxation parameter
185 par_poisson1.add_double(precis_poisson, 1) ; // required precision
186 par_poisson1.add_int_mod(niter, 0) ; // number of iterations actually
187 // used
188 par_poisson1.add_cmp_mod( ssjm1_logn ) ;
189
190 // Parameters for the function Map_et::poisson for beta_auto
191 // ---------------------------------------------------------
192
193 Param par_poisson2 ;
194
195 par_poisson2.add_int(mermax_poisson, 0) ; // maximum number of iterations
196 par_poisson2.add_double(relax_poisson, 0) ; // relaxation parameter
197 par_poisson2.add_double(precis_poisson, 1) ; // required precision
198 par_poisson2.add_int_mod(niter, 0) ; // number of iterations actually
199 // used
200 par_poisson2.add_cmp_mod( ssjm1_beta ) ;
201
202 // Parameters for the function Tenseur::poisson_vect
203 // -------------------------------------------------
204
205 Param par_poisson_vect ;
206
207 par_poisson_vect.add_int(mermax_poisson, 0) ; // maximum number of
208 // iterations
209 par_poisson_vect.add_double(relax_poisson, 0) ; // relaxation parameter
210 par_poisson_vect.add_double(precis_poisson, 1) ; // required precision
211 par_poisson_vect.add_cmp_mod( ssjm1_khi ) ;
212 par_poisson_vect.add_tenseur_mod( ssjm1_wshift ) ;
213 par_poisson_vect.add_int_mod(niter, 0) ;
214
215 // External potential
216 // See Eq (99) from Gourgoulhon et al. (2001)
217 // -----------------------------------------
218
219 Tenseur pot_ext = logn_comp + pot_centri + loggam ;
220
221 Tenseur ent_jm1 = ent ; // Enthalpy at previous step
222
223 Tenseur source(mp) ; // source term in the equation for logn_auto
224 // and beta_auto
225
226 Tenseur source_shift(mp, 1, CON, ref_triad) ; // source term in the
227 // equation for shift_auto
228
229 // maximum l to evaluate multipole moments y_lm
230 int l_max = 4 ;
231
232 if (l_max > 6) {
233 cout << "We don't have those ylm programmed yet!!!!" << endl ;
234 abort() ;
235 }
236
237 int nylm = (l_max+1) * (l_max+1) ;
238
239 Cmp** ylmvec = new Cmp* [nylm] ;
240
241 //==========================================================//
242 // Start of iteration //
243 //==========================================================//
244
245 for(int mer=0 ; mer<mermax ; mer++ ) {
246
247 cout << "-----------------------------------------------" << endl ;
248 cout << "step: " << mer << endl ;
249 cout << "diff_ent = " << diff_ent << endl ;
250
251 //------------------------------------------------------
252 // Resolution of the elliptic equation for the velocity
253 // scalar potential
254 //------------------------------------------------------
255
256 if (irrotational) {
257 diff_vel_pot = velocity_pot_extr(mass, sepa, mermax_potvit,
258 precis_poisson, relax_potvit) ;
259 }
260
261 //-------------------------------------
262 // Computation of the new radial scale
263 //-------------------------------------
264
265 // alpha_r (r = alpha_r r') is determined so that the enthalpy
266 // takes the requested value ent_b at the stellar surface
267
268 // Values at the center of the star:
269 double logn_auto_c = logn_auto()(0, 0, 0, 0) ;
270 double pot_ext_c = pot_ext()(0, 0, 0, 0) ;
271
272 // Search for the reference point (theta_*, phi_*)
273 // [notation of Bonazzola, Gourgoulhon & Marck PRD 58, 104020 (1998)]
274 // at the surface of the star
275 // ------------------------------------------------------------------
276 double alpha_r2 = 0 ;
277 for (int k=0; k<mg->get_np(l_b); k++) {
278 for (int j=0; j<mg->get_nt(l_b); j++) {
279
280 double pot_ext_b = pot_ext()(l_b, k, j, i_b) ;
281 double logn_auto_b = logn_auto()(l_b, k, j, i_b) ;
282
283 // See Eq (100) from Gourgoulhon et al. (2001)
284 double alpha_r2_jk = ( ent_c - ent_b + pot_ext_c - pot_ext_b)
285 / ( logn_auto_b - logn_auto_c ) ;
286
287 if (alpha_r2_jk > alpha_r2) {
288 alpha_r2 = alpha_r2_jk ;
289 k_b = k ;
290 j_b = j ;
291 }
292
293 }
294 }
295
296 alpha_r = sqrt(alpha_r2) ;
297
298 cout << "k_b, j_b, alpha_r: " << k_b << " " << j_b << " "
299 << alpha_r << endl ;
300
301 // New value of logn_auto
302 // ----------------------
303
304 logn_auto = alpha_r2 * logn_auto ;
305 logn_auto_regu = alpha_r2 * logn_auto_regu ;
306 logn_auto_c = logn_auto()(0, 0, 0, 0) ;
307
308 //------------------------------------------------------------
309 // Change the values of the inner points of the second domain
310 // by those of the outer points of the first domain
311 //------------------------------------------------------------
312
313 (logn_auto().va).smooth(nzet, (logn_auto.set()).va) ;
314
315 //--------------------------------------------
316 // First integral --> enthalpy in all space
317 // See Eq (98) from Gourgoulhon et al. (2001)
318 //--------------------------------------------
319
320 ent = (ent_c + logn_auto_c + pot_ext_c) - logn_auto - pot_ext ;
321
322 //----------------------------------------------------------
323 // Change the enthalpy field to be set its maximum position
324 // at the coordinate center
325 //----------------------------------------------------------
326
327 double dentdx = ent().dsdx()(0, 0, 0, 0) ;
328 double dentdy = ent().dsdy()(0, 0, 0, 0) ;
329
330 cout << "dH/dx|_center = " << dentdx << endl ;
331 cout << "dH/dy|_center = " << dentdy << endl ;
332
333 double dec_fact = 1. ;
334
335 Tenseur func_in(mp) ;
336 func_in.set_etat_qcq() ;
337 func_in.set() = 1. - dec_fact * (dentdx/ent_c) * mp.x
338 - dec_fact * (dentdy/ent_c) * mp.y ;
339 func_in.set().annule(nzet, nz-1) ;
340 func_in.set_std_base() ;
341
342 Tenseur func_ex(mp) ;
343 func_ex.set_etat_qcq() ;
344 func_ex.set() = 1. ;
345 func_ex.set().annule(0, nzet-1) ;
346 func_ex.set_std_base() ;
347
348 // New enthalpy field
349 // ------------------
350 ent.set() = ent() * (func_in() + func_ex()) ;
351
352 (ent().va).smooth(nzet, (ent.set()).va) ;
353
354 double dentdx_new = ent().dsdx()(0, 0, 0, 0) ;
355 double dentdy_new = ent().dsdy()(0, 0, 0, 0) ;
356 cout << "dH/dx|_new = " << dentdx_new << endl ;
357 cout << "dH/dy|_new = " << dentdy_new << endl ;
358
359 //-----------------------------------------------------
360 // Adaptation of the mapping to the new enthalpy field
361 //----------------------------------------------------
362
363 // Shall the adaptation be performed (cusp) ?
364 // ------------------------------------------
365
366 double dent_eq = ent().dsdr().val_point(ray_eq_pi(),M_PI/2.,M_PI) ;
367 double dent_pole = ent().dsdr().val_point(ray_pole(),0.,0.) ;
368 double rap_dent = fabs( dent_eq / dent_pole ) ;
369 cout << "| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
370
371 if ( rap_dent < thres_adapt ) {
372 adapt_flag = 0 ; // No adaptation of the mapping
373 cout << "******* FROZEN MAPPING *********" << endl ;
374 }
375 else{
376 adapt_flag = 1 ; // The adaptation of the mapping is to be
377 // performed
378 }
379
380 ent_limit.set_etat_qcq() ;
381 for (int l=0; l<nzet; l++) { // loop on domains inside the star
382 ent_limit.set(l) = ent()(l, k_b, j_b, i_b) ;
383 }
384
385 ent_limit.set(nzet-1) = ent_b ;
386 Map_et mp_prev = mp_et ;
387
388 mp.adapt(ent(), par_adapt) ;
389
390 //----------------------------------------------------
391 // Computation of the enthalpy at the new grid points
392 //----------------------------------------------------
393
394 mp_prev.homothetie(alpha_r) ;
395
396 mp.reevaluate(&mp_prev, nzet+1, ent.set()) ;
397
398 double fact_resize = 1. / alpha_r ;
399 for (int l=nzet; l<nz-1; l++) {
400 mp_et.resize(l, fact_resize) ;
401 }
402 mp_et.resize_extr(fact_resize) ;
403
404 double ent_s_max = -1 ;
405 int k_s_max = -1 ;
406 int j_s_max = -1 ;
407 for (int k=0; k<mg->get_np(l_b); k++) {
408 for (int j=0; j<mg->get_nt(l_b); j++) {
409 double xx = fabs( ent()(l_b, k, j, i_b) ) ;
410 if (xx > ent_s_max) {
411 ent_s_max = xx ;
412 k_s_max = k ;
413 j_s_max = j ;
414 }
415 }
416 }
417 cout << "Max. abs(enthalpy) at the boundary between domains nzet-1"
418 << " and nzet : " << endl ;
419 cout << " " << ent_s_max << " reached for k = " << k_s_max
420 << " and j = " << j_s_max << endl ;
421
422
423 //******************************************
424 // Call this each time the mapping changes,
425 // but no more often than that
426 //******************************************
427
428 int oldindex = -1 ;
429 for (int l=0; l<=l_max; l++) {
430 for (int m=0; m<= l; m++) {
431
432 int index = l*l + 2*m ;
433 if(m > 0) index -= 1 ;
434 if(index != oldindex+1) {
435 cout << "error!! " << l << " " << m << " " << index
436 << " " << oldindex << endl ;
437 abort() ;
438 }
439 // set up Cmp's in ylmvec with proper parity
440
441 if ((l+m) % 2 == 0) {
442 ylmvec[index] = new Cmp (logn_auto()) ;
443 } else {
444 ylmvec[index] = new Cmp (shift_auto(2)) ;
445 }
446 oldindex = index ;
447 if (m > 0) {
448 index += 1 ;
449 if((l+m) % 2 == 0) {
450 ylmvec[index] = new Cmp (logn_auto()) ;
451 } else {
452 ylmvec[index] = new Cmp (shift_auto(2)) ;
453 }
454 oldindex = index ;
455 }
456 }
457 }
458 if(oldindex+1 != nylm) {
459 cout << "ERROR!!! " << oldindex << " "<< nylm << endl ;
460 abort() ;
461 }
462
463 // This may need to be in some class definition,
464 // it references the relevant map_et
465
466 get_ylm(nylm, ylmvec);
467
468 //-------------------
469 // Equation of state
470 //-------------------
471
472 equation_of_state() ; // computes new values for nbar (n), ener (e)
473 // and press (p) from the new ent (H)
474
475 //----------------------------------------------------------
476 // Matter source terms in the gravitational field equations
477 //---------------------------------------------------------
478
479 hydro_euler_extr(mass, sepa) ; // computes new values for
480 // ener_euler (E), s_euler (S)
481 // and u_euler (U^i)
482
483 //----------------------------------------------------
484 // Auxiliary terms for the source of Poisson equation
485 //----------------------------------------------------
486
487 const Coord& xx = mp.x ;
488 const Coord& yy = mp.y ;
489 const Coord& zz = mp.z ;
490
491 Tenseur r_bh(mp) ;
492 r_bh.set_etat_qcq() ;
493 r_bh.set() = pow( (xx+sepa)*(xx+sepa) + yy*yy + zz*zz, 0.5) ;
494 r_bh.set_std_base() ;
495
496 Tenseur xx_cov(mp, 1, COV, ref_triad) ;
497 xx_cov.set_etat_qcq() ;
498 xx_cov.set(0) = xx + sepa ;
499 xx_cov.set(1) = yy ;
500 xx_cov.set(2) = zz ;
501 xx_cov.set_std_base() ;
502
503 Tenseur xsr_cov(mp, 1, COV, ref_triad) ;
504 xsr_cov = xx_cov / r_bh ;
505 xsr_cov.set_std_base() ;
506
507 Tenseur msr(mp) ;
508 msr = ggrav * mass / r_bh ;
509 msr.set_std_base() ;
510
511 Tenseur lapse_bh(mp) ;
512 lapse_bh = 1. / sqrt( 1.+2.*msr ) ;
513 lapse_bh.set_std_base() ;
514
515 Tenseur lapse_bh2(mp) ; // lapse_bh * lapse_bh
516 lapse_bh2 = 1. / (1.+2.*msr) ;
517 lapse_bh2.set_std_base() ;
518
519 Tenseur ldnu(mp) ;
520 ldnu = flat_scalar_prod_desal(xsr_cov, d_logn_auto) ;
521 ldnu.set_std_base() ;
522
523 Tenseur ldbeta(mp) ;
524 ldbeta = flat_scalar_prod_desal(xsr_cov, d_beta_auto) ;
525 ldbeta.set_std_base() ;
526
527 Tenseur lltkij(mp) ;
528 lltkij.set_etat_qcq() ;
529 lltkij.set() = 0 ;
530 lltkij.set_std_base() ;
531
532 for (int i=0; i<3; i++)
533 for (int j=0; j<3; j++)
534 lltkij.set() += xsr_cov(i) * xsr_cov(j) * tkij_auto(i, j) ;
535
536 Tenseur lshift(mp) ;
537 lshift = flat_scalar_prod_desal(xsr_cov, shift_auto) ;
538 lshift.set_std_base() ;
539
540 Tenseur d_ldnu(mp, 1, COV, ref_triad) ;
541 d_ldnu = ldnu.gradient() ; // (d/dx, d/dy, d/dz)
542 d_ldnu.change_triad(ref_triad) ; // --> (d/dX, d/dY, d/dZ)
543
544 Tenseur ldldnu(mp) ;
545 ldldnu = flat_scalar_prod_desal(xsr_cov, d_ldnu) ;
546 ldldnu.set_std_base() ;
547
548 Tenseur d_ldbeta(mp, 1, COV, ref_triad) ;
549 d_ldbeta = ldbeta.gradient() ; // (d/dx, d/dy, d/dz)
550 d_ldbeta.change_triad(ref_triad) ; // --> (d/dX, d/dY, d/dZ)
551
552 Tenseur ldldbeta(mp) ;
553 ldldbeta = flat_scalar_prod_desal(xsr_cov, d_ldbeta) ;
554 ldldbeta.set_std_base() ;
555
556
557 //------------------------------------------
558 // Poisson equation for logn_auto (nu_auto)
559 //------------------------------------------
560
561 // Source
562 // ------
563
564 if (relativistic) {
565
566 source = qpig * a_car % (ener_euler + s_euler) + akcar_auto
567 - flat_scalar_prod_desal(d_logn_auto, d_beta_auto)
568 + 2.*lapse_bh2 % msr % (ldnu % ldbeta + ldldnu)
569 + lapse_bh2 % lapse_bh2 % msr % (2.*(ldnu + 4.*msr % ldnu)
570 - ldbeta) / r_bh
571 - (4.*a_car % lapse_bh2 % lapse_bh2 % msr / 3. / nnn / r_bh)
572 * (2.+3.*msr) * (3.+4.*msr) * lltkij
573 + (2.*a_car % lapse_bh2 % lapse_bh2 % lapse_bh % msr
574 / nnn / r_bh / r_bh) * (2.+10.*msr+9.*msr%msr) * lshift
575 + (4.*pow(lapse_bh2, 3.) % msr % msr / 3. / r_bh / r_bh)
576 % (2.*(a_car%lapse_bh2/nnn/nnn - 1.) * pow(2.+3.*msr, 2.)
577 + (a_car - 1.) % pow(1.+3.*msr, 2.)
578 - 3.*(a_car%lapse_bh/nnn - 1.)*(2.+10.*msr+9.*msr%msr)) ;
579
580 }
581 else {
582 cout << "The computation of BH-NS binary systems"
583 << " should be relativistic !!!" << endl ;
584 abort() ;
585 }
586
587 source.set_std_base() ;
588
589 // Resolution of the Poisson equation
590 // ----------------------------------
591
592 double* intvec_logn = new double[nylm] ;
593 // get_integrals needs access to map_et
594
595 get_integrals(nylm, intvec_logn, ylmvec, source()) ;
596
597 // relax if you have a previously calculated set of moments
598 if(nu_int[0] != -1.0) {
599 for (int i=0; i<nylm; i++) {
600 intvec_logn[i] *= relax_ylm ;
601 intvec_logn[i] += (1.0 - relax_ylm) * nu_int[i] ;
602 }
603 }
604
605 source().poisson_ylm(par_poisson1, logn_auto.set(), nylm,
606 intvec_logn) ;
607
608 for (int i=0; i<nylm; i++) {
609 nu_int[i] = intvec_logn[i] ;
610 }
611
612 delete [] intvec_logn ;
613
614 // Construct logn_auto_regu for et_bin_upmetr_extr.C
615 // -------------------------------------------------
616
617 logn_auto_regu = logn_auto ;
618
619 // Check: has the Poisson equation been correctly solved ?
620 // -----------------------------------------------------
621
622 Tbl tdiff_logn = diffrel(logn_auto().laplacien(), source()) ;
623
624 cout <<
625 "Relative error in the resolution of the equation for logn_auto : "
626 << endl ;
627
628 for (int l=0; l<nz; l++) {
629 cout << tdiff_logn(l) << " " ;
630 }
631 cout << endl ;
632 diff_logn = max(abs(tdiff_logn)) ;
633
634 if (relativistic) {
635
636 //--------------------------------
637 // Poisson equation for beta_auto
638 //--------------------------------
639
640 // Source
641 // ------
642
643 source = qpig * a_car % s_euler + 0.75 * akcar_auto
644 - 0.5 * flat_scalar_prod_desal(d_logn_auto, d_logn_auto)
645 - 0.5 * flat_scalar_prod_desal(d_beta_auto, d_beta_auto)
646 + lapse_bh2 % msr % (ldnu%ldnu + ldbeta%ldbeta + 2.*ldldbeta)
647 + lapse_bh2 % lapse_bh2 % msr % (2.*(1.+4.*msr) * ldbeta
648 - ldnu) / r_bh
649 - (a_car % lapse_bh2 % lapse_bh2 % msr / nnn / r_bh)
650 * (2.+3.*msr) * (3.+4.*msr) * lltkij
651 + (2.*a_car % lapse_bh2 % lapse_bh2 % lapse_bh % msr
652 / nnn / r_bh / r_bh) * (2.+10.*msr+9.*msr%msr) * lshift
653 + (2.*pow(lapse_bh2, 3.) % msr % msr / r_bh / r_bh)
654 % ((a_car%lapse_bh2/nnn/nnn - 1.) * pow(2.+3.*msr, 2.)
655 + (a_car - 1.) * pow(1.+3.*msr, 2.)
656 - 2.*(a_car%lapse_bh/nnn - 1.)*(2.+10.*msr+9.*msr%msr)) ;
657
658 source.set_std_base() ;
659
660 // Resolution of the Poisson equation
661 // ----------------------------------
662
663 double* intvec_beta = new double[nylm] ;
664 // get_integrals needs access to map_et
665
666 get_integrals(nylm, intvec_beta, ylmvec, source()) ;
667 // relax if you have a previously calculated set of moments
668 if(beta_int[0] != -1.0) {
669 for (int i=0; i<nylm; i++) {
670 intvec_beta[i] *= relax_ylm ;
671 intvec_beta[i] += (1.0 - relax_ylm) * beta_int[i] ;
672 }
673 }
674
675 source().poisson_ylm(par_poisson2, beta_auto.set(),
676 nylm, intvec_beta) ;
677
678 for (int i=0; i<nylm; i++) {
679 beta_int[i] = intvec_beta[i] ;
680 }
681
682 delete [] intvec_beta ;
683
684 // Check: has the Poisson equation been correctly solved ?
685 // -----------------------------------------------------
686
687 Tbl tdiff_beta = diffrel(beta_auto().laplacien(), source()) ;
688
689 cout << "Relative error in the resolution of the equation for "
690 << "beta_auto : " << endl ;
691 for (int l=0; l<nz; l++) {
692 cout << tdiff_beta(l) << " " ;
693 }
694 cout << endl ;
695 diff_beta = max(abs(tdiff_beta)) ;
696
697 //----------------------------------------
698 // Vector Poisson equation for shift_auto
699 //----------------------------------------
700
701 // Some auxiliary terms for the source
702 // -----------------------------------
703
704 Tenseur xx_con(mp, 1, CON, ref_triad) ;
705 xx_con.set_etat_qcq() ;
706 xx_con.set(0) = xx + sepa ;
707 xx_con.set(1) = yy ;
708 xx_con.set(2) = zz ;
709 xx_con.set_std_base() ;
710
711 Tenseur xsr_con(mp, 1, CON, ref_triad) ;
712 xsr_con = xx_con / r_bh ;
713 xsr_con.set_std_base() ;
714
715 // Components of shift_auto with respect to the Cartesian triad
716 // (d/dx, d/dy, d/dz) of the mapping :
717
718 Tenseur shift_auto_local = shift_auto ;
719 shift_auto_local.change_triad( mp.get_bvect_cart() ) ;
720
721 // Gradient (partial derivatives with respect to the Cartesian
722 // coordinates of the mapping)
723 // dn(i, j) = D_i N^j
724
725 Tenseur dn = shift_auto_local.gradient() ;
726
727 // Return to the absolute reference frame
728 dn.change_triad(ref_triad) ;
729
730 // Trace of D_i N^j = divergence of N^j :
731 Tenseur divn = contract(dn, 0, 1) ;
732
733 // l^j D_j N^i
734 Tenseur ldn_con = contract(xsr_con, 0, dn, 0) ;
735
736 // D_j (l^k D_k N^i): dldn(j, i)
737 Tenseur ldn_local = ldn_con ;
738 ldn_local.change_triad( mp.get_bvect_cart() ) ;
739 Tenseur dldn = ldn_local.gradient() ;
740 dldn.change_triad(ref_triad) ;
741
742 // l^j D_j (l^k D_k N^i)
743 Tenseur ldldn = contract(xsr_con, 0, dldn, 0) ;
744
745 // l_k D_j N^k
746 Tenseur ldn_cov = contract(xsr_cov, 0, dn, 1) ;
747
748 // l^j l_k D_j N^k
749 Tenseur lldn_cov = contract(xsr_con, 0, ldn_cov, 0) ;
750
751 // eta^{ij} l_k D_j N^k
752 Tenseur eldn(mp, 1, CON, ref_triad) ;
753 eldn.set_etat_qcq() ;
754 eldn.set(0) = ldn_cov(0) ;
755 eldn.set(1) = ldn_cov(1) ;
756 eldn.set(2) = ldn_cov(2) ;
757 eldn.set_std_base() ;
758
759 // l^i D_j N^j
760 Tenseur ldivn = xsr_con % divn ;
761
762 // D_j (l^i D_k N^k): dldivn(j, i)
763 Tenseur ldivn_local = ldivn ;
764 ldivn_local.change_triad( mp.get_bvect_cart() ) ;
765 Tenseur dldivn = ldivn_local.gradient() ;
766 dldivn.change_triad(ref_triad) ;
767
768 // l^j D_j (l^i D_k N^k)
769 Tenseur ldldivn = contract(xsr_con, 0, dldivn, 0) ;
770
771 // l_j N^j
772 Tenseur ln = contract(xsr_cov, 0, shift_auto, 0) ;
773
774 Tenseur vtmp = 6. * d_beta_auto - 8. * d_logn_auto ;
775
776 Tenseur lvtmp = contract(xsr_con, 0, vtmp, 0) ;
777
778 // eta^{ij} vtmp_j
779 Tenseur evtmp(mp, 1, CON, ref_triad) ;
780 evtmp.set_etat_qcq() ;
781 evtmp.set(0) = vtmp(0) ;
782 evtmp.set(1) = vtmp(1) ;
783 evtmp.set(2) = vtmp(2) ;
784 evtmp.set_std_base() ;
785
786 // lapse_ns
787 Tenseur lapse_ns(mp) ;
788 lapse_ns = exp(logn_auto) ;
789 lapse_ns.set_std_base() ;
790
791 // Source
792 // ------
793
794 source_shift = (-4.*qpig) * nnn % a_car % (ener_euler + press)
795 % u_euler
796 + nnn % flat_scalar_prod_desal(tkij_auto, vtmp)
797 - 2.*nnn % lapse_bh2 * msr / r_bh
798 % flat_scalar_prod_desal(tkij_auto, xsr_cov)
799 + 2.*lapse_bh2 * msr * (3.*ldldn + ldldivn) / 3.
800 - lapse_bh2 * msr / r_bh
801 * (4.*ldivn - lapse_bh2 % (3.*ldn_con + 8.*msr * ldn_con)
802 - (eldn + 2.*lapse_bh2*(9.+11.*msr)*lldn_cov%xsr_con) / 3.)
803 - 2.*lapse_bh2 % lapse_bh2 * msr / r_bh / r_bh
804 * ( (4.+11.*msr) * shift_auto
805 - lapse_bh2 * (12.+51.*msr+46.*msr*msr) * ln % xsr_con )
806 / 3.
807 + 8.*pow(lapse_bh2, 4.) * msr / r_bh / r_bh
808 % (lapse_ns - 1.) * (2.+10.*msr+9.*msr*msr) * xsr_con / 3.
809 + 2.*pow(lapse_bh2, 3.) * msr / r_bh * (2.+3.*msr)
810 * ( (1.+2.*msr) * evtmp - (3.+2.*msr) * lvtmp * xsr_con) / 3. ;
811
812 source_shift.set_std_base() ;
813
814 // Resolution of the Poisson equation
815 // ----------------------------------
816
817 // Filter for the source of shift vector :
818
819 for (int i=0; i<3; i++) {
820 for (int l=0; l<nz; l++) {
821 if (source_shift(i).get_etat() != ETATZERO)
822 source_shift.set(i).filtre_phi(np_filter, l) ;
823 }
824 }
825
826 // For Tenseur::poisson_vect, the triad must be the mapping
827 // triad, not the reference one:
828
829 source_shift.change_triad( mp.get_bvect_cart() ) ;
830 /*
831 for (int i=0; i<3; i++) {
832 if(source_shift(i).dz_nonzero()) {
833 assert( source_shift(i).get_dzpuis() == 4 ) ;
834 }
835 else {
836 (source_shift.set(i)).set_dzpuis(4) ;
837 }
838 }
839
840 source_shift.dec2_dzpuis() ; // dzpuis 4 -> 2
841 */
842 double lambda_shift = double(1) / double(3) ;
843
844 double* intvec_shift = new double[4*nylm] ;
845 for (int i=0; i<3; i++) {
846 double* intvec = new double[nylm];
847 get_integrals(nylm, intvec, ylmvec, source_shift(i));
848
849 for (int j=0; j<nylm; j++) {
850 intvec_shift[i*nylm+j] = intvec[j] ;
851 }
852 if(shift_int[i*nylm] != -1.0) {
853 for (int j=0; j<nylm; j++) {
854 intvec_shift[i*nylm+j] *= relax_ylm ;
855 intvec_shift[i*nylm+j] += (1.0 - relax_ylm)
856 * shift_int[i*nylm+j] ;
857 }
858 }
859 delete [] intvec ;
860 }
861 Tenseur source_scal (-skxk(source_shift)) ;
862 Cmp soscal (source_scal()) ;
863 double* intvec2 = new double[nylm] ;
864 get_integrals(nylm, intvec2, ylmvec, soscal) ;
865
866 for (int j=0; j<nylm; j++) {
867 intvec_shift[3*nylm+j] = intvec2[j] ;
868 }
869 if(shift_int[3*nylm] != -1.0) {
870 for (int i=0; i<nylm; i++) {
871 intvec_shift[3*nylm+i] *= relax_ylm ;
872 intvec_shift[3*nylm+i] += (1.0 - relax_ylm)
873 *shift_int[3*nylm+i] ;
874 }
875 }
876
877 delete [] intvec2 ;
878
879 source_shift.poisson_vect_ylm(lambda_shift, par_poisson_vect,
880 shift_auto, w_shift,
881 khi_shift, nylm, intvec_shift) ;
882
883 for (int i=0; i<4*nylm; i++) {
884 shift_int[i] = intvec_shift[i] ;
885 }
886
887 delete[] intvec_shift ;
888
889 // Check: has the equation for shift_auto been correctly solved ?
890 // --------------------------------------------------------------
891
892 // Divergence of shift_auto :
893 Tenseur divna = contract(shift_auto.gradient(), 0, 1) ;
894 // divna.dec2_dzpuis() ; // dzpuis 2 -> 0
895
896 // Grad(div) :
897 Tenseur graddivn = divna.gradient() ;
898 // graddivn.inc2_dzpuis() ; // dzpuis 2 -> 4
899
900 // Full operator :
901 Tenseur lap_shift(mp, 1, CON, mp.get_bvect_cart() ) ;
902 lap_shift.set_etat_qcq() ;
903 for (int i=0; i<3; i++) {
904 lap_shift.set(i) = shift_auto(i).laplacien()
905 + lambda_shift * graddivn(i) ;
906 }
907
908 Tbl tdiff_shift_x = diffrel(lap_shift(0), source_shift(0)) ;
909 Tbl tdiff_shift_y = diffrel(lap_shift(1), source_shift(1)) ;
910 Tbl tdiff_shift_z = diffrel(lap_shift(2), source_shift(2)) ;
911
912 cout << "Relative error in the resolution of the equation for "
913 << "shift_auto : " << endl ;
914 cout << "x component : " ;
915 for (int l=0; l<nz; l++) {
916 cout << tdiff_shift_x(l) << " " ;
917 }
918 cout << endl ;
919 cout << "y component : " ;
920 for (int l=0; l<nz; l++) {
921 cout << tdiff_shift_y(l) << " " ;
922 }
923 cout << endl ;
924 cout << "z component : " ;
925 for (int l=0; l<nz; l++) {
926 cout << tdiff_shift_z(l) << " " ;
927 }
928 cout << endl ;
929
930 diff_shift_x = max(abs(tdiff_shift_x)) ;
931 diff_shift_y = max(abs(tdiff_shift_y)) ;
932 diff_shift_z = max(abs(tdiff_shift_z)) ;
933
934 // Final result
935 // ------------
936 // The output of Tenseur::poisson_vect_falloff is on the mapping
937 // triad, it should therefore be transformed to components on the
938 // reference triad :
939
940 shift_auto.change_triad( ref_triad ) ;
941
942 } // End of relativistic equations
943
944 //------------------------------
945 // Relative change in enthalpy
946 //------------------------------
947
948 Tbl diff_ent_tbl = diffrel( ent(), ent_jm1() ) ;
949 diff_ent = diff_ent_tbl(0) ;
950 for (int l=1; l<nzet; l++) {
951 diff_ent += diff_ent_tbl(l) ;
952 }
953 diff_ent /= nzet ;
954
955 ent_jm1 = ent ;
956
957 for (int i=0; i<nylm; i++) {
958 delete ylmvec[i];
959 }
960
961 } // End of main loop
962
963 //========================================================//
964 // End of iteration //
965 //========================================================//
966
967 delete[] ylmvec ;
968
969}
970
971
972void Et_bin_bhns_extr::equil_bhns_extr_ylm_cf(double ent_c,
973 const double& mass,
974 const double& sepa,
975 double* nu_int,
976 double* beta_int,
977 double* shift_int,
978 int mermax,
979 int mermax_poisson,
980 double relax_poisson,
981 double relax_ylm,
982 int mermax_potvit,
983 double relax_potvit,
984 int np_filter,
985 double thres_adapt,
986 Tbl& diff) {
987
988 // Fundamental constants and units
989 // -------------------------------
990 using namespace Unites ;
991
992 assert( kerrschild == false ) ;
993
994 // Initializations
995 // ---------------
996
997 const Mg3d* mg = mp.get_mg() ;
998 int nz = mg->get_nzone() ; // total number of domains
999
1000 // The following is required to initialize mp_prev as a Map_et:
1001 Map_et& mp_et = dynamic_cast<Map_et&>(mp) ;
1002
1003 // Domain and radial indices of points at the surface of the star:
1004 int l_b = nzet - 1 ;
1005 int i_b = mg->get_nr(l_b) - 1 ;
1006 int k_b ;
1007 int j_b ;
1008
1009 // Value of the enthalpy defining the surface of the star
1010 double ent_b = 0 ;
1011
1012 // Error indicators
1013 // ----------------
1014
1015 double& diff_ent = diff.set(0) ;
1016 double& diff_vel_pot = diff.set(1) ;
1017 double& diff_logn = diff.set(2) ;
1018 double& diff_beta = diff.set(3) ;
1019 double& diff_shift_x = diff.set(4) ;
1020 double& diff_shift_y = diff.set(5) ;
1021 double& diff_shift_z = diff.set(6) ;
1022
1023 // Parameters for the function Map_et::adapt
1024 // -----------------------------------------
1025
1026 Param par_adapt ;
1027 int nitermax = 100 ;
1028 int niter ;
1029 int adapt_flag = 1 ; // 1 = performs the full computation,
1030 // 0 = performs only the rescaling by
1031 // the factor alpha_r
1032 int nz_search = nzet ; // Number of domains for searching
1033 // the enthalpy isosurfaces
1034 double precis_secant = 1.e-14 ;
1035 double alpha_r ;
1036 double reg_map = 1. ; // 1 = regular mapping, 0 = contracting mapping
1037
1038 Tbl ent_limit(nz) ;
1039
1040 par_adapt.add_int(nitermax, 0) ; // maximum number of iterations to
1041 // locate zeros by the secant method
1042 par_adapt.add_int(nzet, 1) ; // number of domains where the adjustment
1043 // to the isosurfaces of ent is to be
1044 // performed
1045 par_adapt.add_int(nz_search, 2) ; // number of domains to search for
1046 // the enthalpy isosurface
1047 par_adapt.add_int(adapt_flag, 3) ; // 1 = performs the full computation,
1048 // 0 = performs only the rescaling by
1049 // the factor alpha_r
1050 par_adapt.add_int(j_b, 4) ; // theta index of the collocation point
1051 // (theta_*, phi_*)
1052 par_adapt.add_int(k_b, 5) ; // theta index of the collocation point
1053 // (theta_*, phi_*)
1054 par_adapt.add_int_mod(niter, 0) ; // number of iterations actually
1055 // used in the secant method
1056 par_adapt.add_double(precis_secant, 0) ; // required absolute precision in
1057 // the determination of zeros by
1058 // the secant method
1059 par_adapt.add_double(reg_map, 1) ; // 1 = regular mapping,
1060 // 0 = contracting mapping
1061 par_adapt.add_double(alpha_r, 2) ; // factor by which all the radial
1062 // distances will be multiplied
1063 par_adapt.add_tbl(ent_limit, 0) ; // array of values of the field ent
1064 // to define the isosurfaces
1065
1066 // Parameters for the function Map_et::poisson for logn_auto
1067 // ---------------------------------------------------------
1068
1069 double precis_poisson = 1.e-16 ;
1070
1071 Param par_poisson1 ;
1072
1073 par_poisson1.add_int(mermax_poisson, 0) ; // maximum number of iterations
1074 par_poisson1.add_double(relax_poisson, 0) ; // relaxation parameter
1075 par_poisson1.add_double(precis_poisson, 1) ; // required precision
1076 par_poisson1.add_int_mod(niter, 0) ; // number of iterations actually
1077 // used
1078 par_poisson1.add_cmp_mod( ssjm1_logn ) ;
1079
1080 // Parameters for the function Map_et::poisson for beta_auto
1081 // ---------------------------------------------------------
1082
1083 Param par_poisson2 ;
1084
1085 par_poisson2.add_int(mermax_poisson, 0) ; // maximum number of iterations
1086 par_poisson2.add_double(relax_poisson, 0) ; // relaxation parameter
1087 par_poisson2.add_double(precis_poisson, 1) ; // required precision
1088 par_poisson2.add_int_mod(niter, 0) ; // number of iterations actually
1089 // used
1090 par_poisson2.add_cmp_mod( ssjm1_beta ) ;
1091
1092 // Parameters for the function Tenseur::poisson_vect
1093 // -------------------------------------------------
1094
1095 Param par_poisson_vect ;
1096
1097 par_poisson_vect.add_int(mermax_poisson, 0) ; // maximum number of
1098 // iterations
1099 par_poisson_vect.add_double(relax_poisson, 0) ; // relaxation parameter
1100 par_poisson_vect.add_double(precis_poisson, 1) ; // required precision
1101 par_poisson_vect.add_cmp_mod( ssjm1_khi ) ;
1102 par_poisson_vect.add_tenseur_mod( ssjm1_wshift ) ;
1103 par_poisson_vect.add_int_mod(niter, 0) ;
1104
1105 // External potential
1106 // See Eq (99) from Gourgoulhon et al. (2001)
1107 // -----------------------------------------
1108
1109 Tenseur pot_ext = logn_comp + pot_centri + loggam ;
1110
1111 Tenseur ent_jm1 = ent ; // Enthalpy at previous step
1112
1113 Tenseur source(mp) ; // source term in the equation for logn_auto
1114 // and beta_auto
1115
1116 Tenseur source_shift(mp, 1, CON, ref_triad) ; // source term in the
1117 // equation for shift_auto
1118
1119 // maximum l to evaluate multipole moments y_lm
1120 int l_max = 4 ;
1121
1122 if (l_max > 6) {
1123 cout << "We don't have those ylm programmed yet!!!!" << endl ;
1124 abort() ;
1125 }
1126
1127 int nylm = (l_max+1) * (l_max+1) ;
1128
1129 Cmp** ylmvec = new Cmp* [nylm] ;
1130
1131 //==========================================================//
1132 // Start of iteration //
1133 //==========================================================//
1134
1135 for(int mer=0 ; mer<mermax ; mer++ ) {
1136
1137 cout << "-----------------------------------------------" << endl ;
1138 cout << "step: " << mer << endl ;
1139 cout << "diff_ent = " << diff_ent << endl ;
1140
1141 //------------------------------------------------------
1142 // Resolution of the elliptic equation for the velocity
1143 // scalar potential
1144 //------------------------------------------------------
1145
1146 if (irrotational) {
1147 diff_vel_pot = velocity_pot_extr(mass, sepa, mermax_potvit,
1148 precis_poisson, relax_potvit) ;
1149 }
1150
1151 //-------------------------------------
1152 // Computation of the new radial scale
1153 //-------------------------------------
1154
1155 // alpha_r (r = alpha_r r') is determined so that the enthalpy
1156 // takes the requested value ent_b at the stellar surface
1157
1158 // Values at the center of the star:
1159 double logn_auto_c = logn_auto()(0, 0, 0, 0) ;
1160 double pot_ext_c = pot_ext()(0, 0, 0, 0) ;
1161
1162 // Search for the reference point (theta_*, phi_*)
1163 // [notation of Bonazzola, Gourgoulhon & Marck PRD 58, 104020 (1998)]
1164 // at the surface of the star
1165 // ------------------------------------------------------------------
1166 double alpha_r2 = 0 ;
1167 for (int k=0; k<mg->get_np(l_b); k++) {
1168 for (int j=0; j<mg->get_nt(l_b); j++) {
1169
1170 double pot_ext_b = pot_ext()(l_b, k, j, i_b) ;
1171 double logn_auto_b = logn_auto()(l_b, k, j, i_b) ;
1172
1173 // See Eq (100) from Gourgoulhon et al. (2001)
1174 double alpha_r2_jk = ( ent_c - ent_b + pot_ext_c - pot_ext_b)
1175 / ( logn_auto_b - logn_auto_c ) ;
1176
1177 if (alpha_r2_jk > alpha_r2) {
1178 alpha_r2 = alpha_r2_jk ;
1179 k_b = k ;
1180 j_b = j ;
1181 }
1182
1183 }
1184 }
1185
1186 alpha_r = sqrt(alpha_r2) ;
1187
1188 cout << "k_b, j_b, alpha_r: " << k_b << " " << j_b << " "
1189 << alpha_r << endl ;
1190
1191 // New value of logn_auto
1192 // ----------------------
1193
1194 logn_auto = alpha_r2 * logn_auto ;
1195 logn_auto_regu = alpha_r2 * logn_auto_regu ;
1196 logn_auto_c = logn_auto()(0, 0, 0, 0) ;
1197
1198 //------------------------------------------------------------
1199 // Change the values of the inner points of the second domain
1200 // by those of the outer points of the first domain
1201 //------------------------------------------------------------
1202
1203 (logn_auto().va).smooth(nzet, (logn_auto.set()).va) ;
1204
1205 //--------------------------------------------
1206 // First integral --> enthalpy in all space
1207 // See Eq (98) from Gourgoulhon et al. (2001)
1208 //--------------------------------------------
1209
1210 ent = (ent_c + logn_auto_c + pot_ext_c) - logn_auto - pot_ext ;
1211
1212 //---------------------------------------------------------
1213 // Change the enthalpy field to accelerate the convergence
1214 //---------------------------------------------------------
1215
1216 double dentdx = ent().dsdx()(0, 0, 0, 0) ;
1217 double dentdy = ent().dsdy()(0, 0, 0, 0) ;
1218
1219 cout << "dH/dx|_center = " << dentdx << endl ;
1220 cout << "dH/dy|_center = " << dentdy << endl ;
1221
1222 double dec_fact = 1. ;
1223
1224 Tenseur func_in(mp) ;
1225 func_in.set_etat_qcq() ;
1226 func_in.set() = 1. - dec_fact * (dentdx/ent_c) * mp.x ;
1227 func_in.set().annule(nzet, nz-1) ;
1228 func_in.set_std_base() ;
1229
1230 Tenseur func_ex(mp) ;
1231 func_ex.set_etat_qcq() ;
1232 func_ex.set() = 1. ;
1233 func_ex.set().annule(0, nzet-1) ;
1234 func_ex.set_std_base() ;
1235
1236 // New enthalpy field
1237 // ------------------
1238 ent.set() = ent() * (func_in() + func_ex()) ;
1239
1240 (ent().va).smooth(nzet, (ent.set()).va) ;
1241
1242 double dentdx_new = ent().dsdx()(0, 0, 0, 0) ;
1243
1244 cout << "dH/dx|_new = " << dentdx_new << endl ;
1245
1246 //-----------------------------------------------------
1247 // Adaptation of the mapping to the new enthalpy field
1248 //----------------------------------------------------
1249
1250 // Shall the adaptation be performed (cusp) ?
1251 // ------------------------------------------
1252
1253 double dent_eq = ent().dsdr().val_point(ray_eq_pi(),M_PI/2.,M_PI) ;
1254 double dent_pole = ent().dsdr().val_point(ray_pole(),0.,0.) ;
1255 double rap_dent = fabs( dent_eq / dent_pole ) ;
1256 cout << "| dH/dr_eq / dH/dr_pole | = " << rap_dent << endl ;
1257
1258 if ( rap_dent < thres_adapt ) {
1259 adapt_flag = 0 ; // No adaptation of the mapping
1260 cout << "******* FROZEN MAPPING *********" << endl ;
1261 }
1262 else{
1263 adapt_flag = 1 ; // The adaptation of the mapping is to be
1264 // performed
1265 }
1266
1267 ent_limit.set_etat_qcq() ;
1268 for (int l=0; l<nzet; l++) { // loop on domains inside the star
1269 ent_limit.set(l) = ent()(l, k_b, j_b, i_b) ;
1270 }
1271
1272 ent_limit.set(nzet-1) = ent_b ;
1273 Map_et mp_prev = mp_et ;
1274
1275 mp.adapt(ent(), par_adapt) ;
1276
1277 //----------------------------------------------------
1278 // Computation of the enthalpy at the new grid points
1279 //----------------------------------------------------
1280
1281 mp_prev.homothetie(alpha_r) ;
1282
1283 mp.reevaluate_symy(&mp_prev, nzet+1, ent.set()) ;
1284
1285 double fact_resize = 1. / alpha_r ;
1286 for (int l=nzet; l<nz-1; l++) {
1287 mp_et.resize(l, fact_resize) ;
1288 }
1289 mp_et.resize_extr(fact_resize) ;
1290
1291 double ent_s_max = -1 ;
1292 int k_s_max = -1 ;
1293 int j_s_max = -1 ;
1294 for (int k=0; k<mg->get_np(l_b); k++) {
1295 for (int j=0; j<mg->get_nt(l_b); j++) {
1296 double xx = fabs( ent()(l_b, k, j, i_b) ) ;
1297 if (xx > ent_s_max) {
1298 ent_s_max = xx ;
1299 k_s_max = k ;
1300 j_s_max = j ;
1301 }
1302 }
1303 }
1304 cout << "Max. abs(enthalpy) at the boundary between domains nzet-1"
1305 << " and nzet : " << endl ;
1306 cout << " " << ent_s_max << " reached for k = " << k_s_max
1307 << " and j = " << j_s_max << endl ;
1308
1309
1310 //******************************************
1311 // Call this each time the mapping changes,
1312 // but no more often than that
1313 //******************************************
1314
1315 int oldindex = -1 ;
1316 for (int l=0; l<=l_max; l++) {
1317 for (int m=0; m<= l; m++) {
1318
1319 int index = l*l + 2*m ;
1320 if(m > 0) index -= 1 ;
1321 if(index != oldindex+1) {
1322 cout << "error!! " << l << " " << m << " " << index
1323 << " " << oldindex << endl ;
1324 abort() ;
1325 }
1326 // set up Cmp's in ylmvec with proper parity
1327
1328 if ((l+m) % 2 == 0) {
1329 ylmvec[index] = new Cmp (logn_auto()) ;
1330 } else {
1331 ylmvec[index] = new Cmp (shift_auto(2)) ;
1332 }
1333 oldindex = index ;
1334 if (m > 0) {
1335 index += 1 ;
1336 if((l+m) % 2 == 0) {
1337 ylmvec[index] = new Cmp (logn_auto()) ;
1338 } else {
1339 ylmvec[index] = new Cmp (shift_auto(2)) ;
1340 }
1341 oldindex = index ;
1342 }
1343 }
1344 }
1345 if(oldindex+1 != nylm) {
1346 cout << "ERROR!!! " << oldindex << " "<< nylm << endl ;
1347 abort() ;
1348 }
1349
1350 // This may need to be in some class definition,
1351 // it references the relevant map_et
1352
1353 get_ylm(nylm, ylmvec);
1354
1355 //-------------------
1356 // Equation of state
1357 //-------------------
1358
1359 equation_of_state() ; // computes new values for nbar (n), ener (e)
1360 // and press (p) from the new ent (H)
1361
1362 //----------------------------------------------------------
1363 // Matter source terms in the gravitational field equations
1364 //---------------------------------------------------------
1365
1366 hydro_euler_extr(mass, sepa) ; // computes new values for
1367 // ener_euler (E), s_euler (S)
1368 // and u_euler (U^i)
1369
1370 //----------------------------------------------------
1371 // Auxiliary terms for the source of Poisson equation
1372 //----------------------------------------------------
1373
1374 const Coord& xx = mp.x ;
1375 const Coord& yy = mp.y ;
1376 const Coord& zz = mp.z ;
1377
1378 Tenseur r_bh(mp) ;
1379 r_bh.set_etat_qcq() ;
1380 r_bh.set() = pow( (xx+sepa)*(xx+sepa) + yy*yy + zz*zz, 0.5) ;
1381 r_bh.set_std_base() ;
1382
1383 Tenseur xx_cov(mp, 1, COV, ref_triad) ;
1384 xx_cov.set_etat_qcq() ;
1385 xx_cov.set(0) = xx + sepa ;
1386 xx_cov.set(1) = yy ;
1387 xx_cov.set(2) = zz ;
1388 xx_cov.set_std_base() ;
1389
1390 Tenseur msr(mp) ;
1391 msr = ggrav * mass / r_bh ;
1392 msr.set_std_base() ;
1393
1394 Tenseur tmp(mp) ;
1395 tmp = 1. / ( 1. - 0.25*msr*msr ) ;
1396 tmp.set_std_base() ;
1397
1398 Tenseur xdnu(mp) ;
1399 xdnu = flat_scalar_prod_desal(xx_cov, d_logn_auto) ;
1400 xdnu.set_std_base() ;
1401
1402 Tenseur xdbeta(mp) ;
1403 xdbeta = flat_scalar_prod_desal(xx_cov, d_beta_auto) ;
1404 xdbeta.set_std_base() ;
1405
1406 //------------------------------------------
1407 // Poisson equation for logn_auto (nu_auto)
1408 //------------------------------------------
1409
1410 // Source
1411 // ------
1412
1413 if (relativistic) {
1414
1415 source = qpig * a_car % (ener_euler + s_euler) + akcar_auto
1416 - flat_scalar_prod_desal(d_logn_auto, d_beta_auto)
1417 - 0.5 * tmp % msr % msr % xdnu / r_bh / r_bh
1418 - tmp % msr % xdbeta / r_bh / r_bh ;
1419
1420 }
1421 else {
1422 cout << "The computation of BH-NS binary systems"
1423 << " should be relativistic !!!" << endl ;
1424 abort() ;
1425 }
1426
1427 source.set_std_base() ;
1428
1429 // Resolution of the Poisson equation
1430 // ----------------------------------
1431
1432 double* intvec_logn = new double[nylm] ;
1433 // get_integrals needs access to map_et
1434
1435 get_integrals(nylm, intvec_logn, ylmvec, source()) ;
1436
1437 // relax if you have a previously calculated set of moments
1438 if(nu_int[0] != -1.0) {
1439 for (int i=0; i<nylm; i++) {
1440 intvec_logn[i] *= relax_ylm ;
1441 intvec_logn[i] += (1.0 - relax_ylm) * nu_int[i] ;
1442 }
1443 }
1444
1445 source().poisson_ylm(par_poisson1, logn_auto.set(), nylm,
1446 intvec_logn) ;
1447
1448 for (int i=0; i<nylm; i++) {
1449 nu_int[i] = intvec_logn[i] ;
1450 }
1451
1452 delete [] intvec_logn ;
1453
1454 // Construct logn_auto_regu for et_bin_upmetr_extr.C
1455 // -------------------------------------------------
1456
1457 logn_auto_regu = logn_auto ;
1458
1459 // Check: has the Poisson equation been correctly solved ?
1460 // -----------------------------------------------------
1461
1462 Tbl tdiff_logn = diffrel(logn_auto().laplacien(), source()) ;
1463
1464 cout <<
1465 "Relative error in the resolution of the equation for logn_auto : "
1466 << endl ;
1467
1468 for (int l=0; l<nz; l++) {
1469 cout << tdiff_logn(l) << " " ;
1470 }
1471 cout << endl ;
1472 diff_logn = max(abs(tdiff_logn)) ;
1473
1474 if (relativistic) {
1475
1476 //--------------------------------
1477 // Poisson equation for beta_auto
1478 //--------------------------------
1479
1480 // Source
1481 // ------
1482
1483 source = qpig * a_car % s_euler + 0.75 * akcar_auto
1484 - 0.5 * flat_scalar_prod_desal(d_logn_auto, d_logn_auto)
1485 - 0.5 * flat_scalar_prod_desal(d_beta_auto, d_beta_auto)
1486 - tmp % msr % xdnu / r_bh / r_bh
1487 - 0.5 * tmp % msr %msr % xdbeta / r_bh / r_bh ;
1488
1489 source.set_std_base() ;
1490
1491 // Resolution of the Poisson equation
1492 // ----------------------------------
1493
1494 double* intvec_beta = new double[nylm] ;
1495 // get_integrals needs access to map_et
1496
1497 get_integrals(nylm, intvec_beta, ylmvec, source()) ;
1498 // relax if you have a previously calculated set of moments
1499 if(beta_int[0] != -1.0) {
1500 for (int i=0; i<nylm; i++) {
1501 intvec_beta[i] *= relax_ylm ;
1502 intvec_beta[i] += (1.0 - relax_ylm) * beta_int[i] ;
1503 }
1504 }
1505
1506 source().poisson_ylm(par_poisson2, beta_auto.set(),
1507 nylm, intvec_beta) ;
1508
1509 for (int i=0; i<nylm; i++) {
1510 beta_int[i] = intvec_beta[i] ;
1511 }
1512
1513 delete [] intvec_beta ;
1514
1515 // Check: has the Poisson equation been correctly solved ?
1516 // -----------------------------------------------------
1517
1518 Tbl tdiff_beta = diffrel(beta_auto().laplacien(), source()) ;
1519
1520 cout << "Relative error in the resolution of the equation for "
1521 << "beta_auto : " << endl ;
1522 for (int l=0; l<nz; l++) {
1523 cout << tdiff_beta(l) << " " ;
1524 }
1525 cout << endl ;
1526 diff_beta = max(abs(tdiff_beta)) ;
1527
1528 //----------------------------------------
1529 // Vector Poisson equation for shift_auto
1530 //----------------------------------------
1531
1532 // Some auxiliary terms for the source
1533 // -----------------------------------
1534
1535 Tenseur bhtmp(mp, 1, COV, ref_triad) ;
1536 bhtmp.set_etat_qcq() ;
1537 bhtmp = tmp % msr % (3.*msr-8.) % xx_cov / r_bh / r_bh ;
1538 bhtmp.set_std_base() ;
1539
1540 Tenseur vtmp = 6. * d_beta_auto - 8. * d_logn_auto ;
1541
1542 // Source
1543 // ------
1544
1545 source_shift = (-4.*qpig) * nnn % a_car % (ener_euler + press)
1546 % u_euler
1547 + nnn % flat_scalar_prod_desal(tkij_auto, vtmp+bhtmp) ;
1548
1549 source_shift.set_std_base() ;
1550
1551 // Resolution of the Poisson equation
1552 // ----------------------------------
1553
1554 // Filter for the source of shift vector :
1555
1556 for (int i=0; i<3; i++) {
1557 for (int l=0; l<nz; l++) {
1558 if (source_shift(i).get_etat() != ETATZERO)
1559 source_shift.set(i).filtre_phi(np_filter, l) ;
1560 }
1561 }
1562
1563 // For Tenseur::poisson_vect, the triad must be the mapping
1564 // triad, not the reference one:
1565
1566 source_shift.change_triad( mp.get_bvect_cart() ) ;
1567 /*
1568 for (int i=0; i<3; i++) {
1569 if(source_shift(i).dz_nonzero()) {
1570 assert( source_shift(i).get_dzpuis() == 4 ) ;
1571 }
1572 else {
1573 (source_shift.set(i)).set_dzpuis(4) ;
1574 }
1575 }
1576
1577 source_shift.dec2_dzpuis() ; // dzpuis 4 -> 2
1578 */
1579 double lambda_shift = double(1) / double(3) ;
1580
1581 double* intvec_shift = new double[4*nylm] ;
1582 for (int i=0; i<3; i++) {
1583 double* intvec = new double[nylm];
1584 get_integrals(nylm, intvec, ylmvec, source_shift(i));
1585
1586 for (int j=0; j<nylm; j++) {
1587 intvec_shift[i*nylm+j] = intvec[j] ;
1588 }
1589 if(shift_int[i*nylm] != -1.0) {
1590 for (int j=0; j<nylm; j++) {
1591 intvec_shift[i*nylm+j] *= relax_ylm ;
1592 intvec_shift[i*nylm+j] += (1.0 - relax_ylm)
1593 * shift_int[i*nylm+j] ;
1594 }
1595 }
1596 delete [] intvec ;
1597 }
1598 Tenseur source_scal (-skxk(source_shift)) ;
1599 Cmp soscal (source_scal()) ;
1600 double* intvec2 = new double[nylm] ;
1601 get_integrals(nylm, intvec2, ylmvec, soscal) ;
1602
1603 for (int j=0; j<nylm; j++) {
1604 intvec_shift[3*nylm+j] = intvec2[j] ;
1605 }
1606 if(shift_int[3*nylm] != -1.0) {
1607 for (int i=0; i<nylm; i++) {
1608 intvec_shift[3*nylm+i] *= relax_ylm ;
1609 intvec_shift[3*nylm+i] += (1.0 - relax_ylm)
1610 *shift_int[3*nylm+i] ;
1611 }
1612 }
1613
1614 delete [] intvec2 ;
1615
1616 source_shift.poisson_vect_ylm(lambda_shift, par_poisson_vect,
1617 shift_auto, w_shift,
1618 khi_shift, nylm, intvec_shift) ;
1619
1620 for (int i=0; i<4*nylm; i++) {
1621 shift_int[i] = intvec_shift[i] ;
1622 }
1623
1624 delete[] intvec_shift ;
1625
1626 // Check: has the equation for shift_auto been correctly solved ?
1627 // --------------------------------------------------------------
1628
1629 // Divergence of shift_auto :
1630 Tenseur divna = contract(shift_auto.gradient(), 0, 1) ;
1631 // divna.dec2_dzpuis() ; // dzpuis 2 -> 0
1632
1633 // Grad(div) :
1634 Tenseur graddivn = divna.gradient() ;
1635 // graddivn.inc2_dzpuis() ; // dzpuis 2 -> 4
1636
1637 // Full operator :
1638 Tenseur lap_shift(mp, 1, CON, mp.get_bvect_cart() ) ;
1639 lap_shift.set_etat_qcq() ;
1640 for (int i=0; i<3; i++) {
1641 lap_shift.set(i) = shift_auto(i).laplacien()
1642 + lambda_shift * graddivn(i) ;
1643 }
1644
1645 Tbl tdiff_shift_x = diffrel(lap_shift(0), source_shift(0)) ;
1646 Tbl tdiff_shift_y = diffrel(lap_shift(1), source_shift(1)) ;
1647 Tbl tdiff_shift_z = diffrel(lap_shift(2), source_shift(2)) ;
1648
1649 cout << "Relative error in the resolution of the equation for "
1650 << "shift_auto : " << endl ;
1651 cout << "x component : " ;
1652 for (int l=0; l<nz; l++) {
1653 cout << tdiff_shift_x(l) << " " ;
1654 }
1655 cout << endl ;
1656 cout << "y component : " ;
1657 for (int l=0; l<nz; l++) {
1658 cout << tdiff_shift_y(l) << " " ;
1659 }
1660 cout << endl ;
1661 cout << "z component : " ;
1662 for (int l=0; l<nz; l++) {
1663 cout << tdiff_shift_z(l) << " " ;
1664 }
1665 cout << endl ;
1666
1667 diff_shift_x = max(abs(tdiff_shift_x)) ;
1668 diff_shift_y = max(abs(tdiff_shift_y)) ;
1669 diff_shift_z = max(abs(tdiff_shift_z)) ;
1670
1671 // Final result
1672 // ------------
1673 // The output of Tenseur::poisson_vect_falloff is on the mapping
1674 // triad, it should therefore be transformed to components on the
1675 // reference triad :
1676
1677 shift_auto.change_triad( ref_triad ) ;
1678
1679 } // End of relativistic equations
1680
1681 //------------------------------
1682 // Relative change in enthalpy
1683 //------------------------------
1684
1685 Tbl diff_ent_tbl = diffrel( ent(), ent_jm1() ) ;
1686 diff_ent = diff_ent_tbl(0) ;
1687 for (int l=1; l<nzet; l++) {
1688 diff_ent += diff_ent_tbl(l) ;
1689 }
1690 diff_ent /= nzet ;
1691
1692 ent_jm1 = ent ;
1693
1694 for (int i=0; i<nylm; i++) {
1695 delete ylmvec[i];
1696 }
1697
1698 } // End of main loop
1699
1700 //========================================================//
1701 // End of iteration //
1702 //========================================================//
1703
1704 delete[] ylmvec ;
1705
1706}
1707
1708}
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
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_ylm_ks
void equil_bhns_extr_ylm_ks(double ent_c, const double &mass, const double &sepa, double *nu_int, double *beta_int, double *shift_int, int mermax, int mermax_poisson, double relax_poisson, double relax_ylm, 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_eq_ylm.C:82
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::get_ylm
void get_ylm(int nylm, Cmp **ylmvec) const
Constructs spherical harmonics.
Definition et_bin_bhns_extr_ylm.C:66
Lorene::Et_bin_bhns_extr::get_integrals
void get_integrals(int nylm, double *intvec, Cmp **ylmvec, Cmp source) const
Computes multipole moments.
Definition et_bin_bhns_extr_ylm.C:258
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::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::Et_bin_bhns_extr::equil_bhns_extr_ylm_cf
void equil_bhns_extr_ylm_cf(double ent_c, const double &mass, const double &sepa, double *nu_int, double *beta_int, double *shift_int, int mermax, int mermax_poisson, double relax_poisson, double relax_ylm, 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_eq_ylm.C:972
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::skxk
Tenseur skxk(const Tenseur &)
Contraction of the last index of (*this) with or , depending on the type of S .
Definition tenseur_operateur.C:583
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.