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