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 char 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 
78 namespace 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
725  dn.change_triad(ref_triad) ;
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
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.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
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 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
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
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
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.