LORENE
isol_hor.C
1 /*
2  * Methods of class Isol_hor
3  *
4  * (see file isol_hor.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2004 Jose Luis Jaramillo
10  * Francois Limousin
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License version 2
16  * as published by the Free Software Foundation.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 
31 /*
32  * $Id: isol_hor.C,v 1.36 2016/12/05 16:17:56 j_novak Exp $
33  * $Log: isol_hor.C,v $
34  * Revision 1.36 2016/12/05 16:17:56 j_novak
35  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
36  *
37  * Revision 1.35 2014/10/13 08:53:01 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.34 2014/10/06 15:13:11 j_novak
41  * Modified #include directives to use c++ syntax.
42  *
43  * Revision 1.33 2009/05/18 22:04:27 j_novak
44  * Changed pow(psi_in, 6) to psi*...*psi in the call to Time_slice_conf constructor. This is to get a well-defined basis.
45  *
46  * Revision 1.32 2008/12/02 15:02:21 j_novak
47  * Implementation of the new constrained formalism, following Cordero et al. 2009
48  * paper. The evolution eqs. are solved as a first-order system. Not tested yet!
49  *
50  * Revision 1.31 2006/05/24 16:55:31 f_limousin
51  * Improvement of dn_comp() and dpsi_comp()
52  *
53  * Revision 1.30 2005/10/21 16:20:55 jl_jaramillo
54  * Version for the paper JaramL05
55  *
56  * Revision 1.29 2005/09/13 18:33:17 f_limousin
57  * New function vv_bound_cart_bin(double) for computing binaries with
58  * berlin condition for the shift vector.
59  * Suppress all the symy and asymy in the importations.
60  *
61  * Revision 1.28 2005/09/12 12:33:54 f_limousin
62  * Compilation Warning - Change of convention for the angular velocity
63  * Add Berlin boundary condition in the case of binary horizons.
64  *
65  * Revision 1.27 2005/07/08 13:15:23 f_limousin
66  * Improvements of boundary_vv_cart(), boundary_nn_lapl().
67  * Add a fonction to compute the departure of axisymmetry.
68  *
69  * Revision 1.26 2005/05/12 14:48:07 f_limousin
70  * New boundary condition for the lapse : boundary_nn_lapl().
71  *
72  * Revision 1.25 2005/04/15 09:34:16 jl_jaramillo
73  * Function "adapt_hor" for adapting the the excised surface to
74  * a given surface. The zero expansion surface is not properly implemented
75  *
76  * Revision 1.24 2005/04/08 12:16:52 f_limousin
77  * Function set_psi(). And dependance in phi.
78  *
79  * Revision 1.23 2005/04/03 19:48:22 f_limousin
80  * Implementation of set_psi(psi_in). And minor changes to avoid warnings.
81  *
82  * Revision 1.22 2005/04/02 15:49:21 f_limousin
83  * New choice (Lichnerowicz) for aaquad. New member data nz.
84  *
85  * Revision 1.21 2005/03/31 09:45:31 f_limousin
86  * New functions compute_ww(...) and aa_kerr_ww().
87  *
88  * Revision 1.20 2005/03/30 12:08:20 f_limousin
89  * Implementation of K^{ij} (Eq.(13) Of Sergio (2002)).
90  *
91  * Revision 1.19 2005/03/28 19:42:39 f_limousin
92  * Implement the metric and A^{ij}A_{ij} of Sergio for pertubations
93  * of Kerr black holes.
94  *
95  * Revision 1.18 2005/03/24 17:05:34 f_limousin
96  * Small change
97  *
98  * Revision 1.17 2005/03/24 16:50:28 f_limousin
99  * Add parameters solve_shift and solve_psi in par_isol.d and in function
100  * init_dat(...). Implement Isolhor::kerr_perturb().
101  *
102  * Revision 1.16 2005/03/10 10:19:42 f_limousin
103  * Add the regularisation of the shift in the case of a single black hole
104  * and lapse zero on the horizon.
105  *
106  * Revision 1.15 2005/03/09 10:29:53 f_limousin
107  * New function update_aa().
108  *
109  * Revision 1.14 2005/03/06 16:59:14 f_limousin
110  * New function Isol_hor::aa() (the one belonging to the class
111  * Time_slice_conf need to compute the time derivative of hh and thus
112  * cannot work in the class Isol_hor).
113  *
114  * Revision 1.13 2005/03/03 15:12:17 f_limousin
115  * Implement function operator>>
116  *
117  * Revision 1.12 2005/03/03 10:05:36 f_limousin
118  * Introduction of members boost_x and boost_z.
119  *
120  * Revision 1.11 2005/02/07 10:35:05 f_limousin
121  * Add the regularisation of the shift for the case N=0 on the horizon.
122  *
123  * Revision 1.10 2004/12/31 15:36:43 f_limousin
124  * Add the constructor from a file and change the standard constructor.
125  *
126  * Revision 1.9 2004/12/29 16:14:22 f_limousin
127  * Add new function beta_comp(const Isol_hor& comp).
128  *
129  * Revision 1.7 2004/11/05 10:57:03 f_limousin
130  * Delete argument partial_save in the function sauve.
131  *
132  * Revision 1.6 2004/11/05 10:10:21 f_limousin
133  * Construction of an isolhor with the Metric met_gamt instead
134  * of a Sym_tensor.
135  *
136  * Revision 1.5 2004/11/03 17:16:06 f_limousin
137  * Change the standart constructor. Add 4 memebers : trK, trK_point,
138  * gamt and gamt_point.
139  * Add also a constructor from a file.
140  *
141  * Revision 1.3 2004/10/29 15:44:45 jl_jaramillo
142  * Remove two members
143  *
144  * Revision 1.2 2004/09/28 16:07:16 f_limousin
145  * Remove all unused functions.
146  *
147  * Revision 1.1 2004/09/09 14:07:26 jl_jaramillo
148  * First version
149  *
150  * Revision 1.1 2004/03/30 14:00:31 jl_jaramillo
151  * New class Isol_hor (first version).
152  *
153  *
154  * $Header: /cvsroot/Lorene/C++/Source/Isol_hor/isol_hor.C,v 1.36 2016/12/05 16:17:56 j_novak Exp $
155  *
156  */
157 
158 // C headers
159 #include <cstdlib>
160 #include <cassert>
161 
162 // Lorene headers
163 #include "param.h"
164 #include "utilitaires.h"
165 #include "time_slice.h"
166 #include "isol_hor.h"
167 #include "tensor.h"
168 #include "metric.h"
169 #include "evolution.h"
170 //#include "graphique.h"
171 
172 //--------------//
173 // Constructors //
174 //--------------//
175 // Standard constructor
176 // --------------------
177 
178 namespace Lorene {
179 Isol_hor::Isol_hor(Map_af& mpi, int depth_in) :
180  Time_slice_conf(mpi, mpi.get_bvect_spher(), mpi.flat_met_spher()),
181  mp(mpi), nz(mpi.get_mg()->get_nzone()), radius ((mpi.get_alpha())[0]),
182  omega(0), boost_x(0), boost_z(0), regul(0),
183  n_auto_evol(depth_in), n_comp_evol(depth_in),
184  psi_auto_evol(depth_in), psi_comp_evol(depth_in),
185  dn_evol(depth_in), dpsi_evol(depth_in),
186  beta_auto_evol(depth_in), beta_comp_evol(depth_in),
187  aa_auto_evol(depth_in), aa_comp_evol(depth_in),
188  aa_nn(depth_in), aa_quad_evol(depth_in),
189  met_gamt(mpi.flat_met_spher()), gamt_point(mpi, CON, mpi.get_bvect_spher()),
190  trK(mpi), trK_point(mpi), decouple(mpi){
191 }
192 
193 // Constructor from conformal decomposition
194 // ----------------------------------------
195 
196 Isol_hor::Isol_hor(Map_af& mpi, const Scalar& lapse_in,
197  const Scalar& psi_in, const Vector& shift_in,
198  const Sym_tensor& aa_in,
199  const Metric& metgamt, const Sym_tensor& gamt_point_in,
200  const Scalar& trK_in, const Scalar& trK_point_in,
201  const Metric_flat& ff_in, int depth_in)
202  : Time_slice_conf(lapse_in, shift_in, ff_in, psi_in, metgamt.con() -
203  ff_in.con(), psi_in*psi_in*psi_in*psi_in*psi_in*psi_in*aa_in,
204  trK_in, depth_in),
205  mp(mpi), nz(mpi.get_mg()->get_nzone()), radius ((mpi.get_alpha())[0]),
206  omega(0), boost_x(0), boost_z(0), regul(0),
207  n_auto_evol(depth_in), n_comp_evol(depth_in),
208  psi_auto_evol(depth_in), psi_comp_evol(depth_in),
209  dn_evol(depth_in), dpsi_evol(depth_in),
210  beta_auto_evol(depth_in), beta_comp_evol(depth_in),
211  aa_auto_evol(depth_in), aa_comp_evol(depth_in),
212  aa_nn(depth_in), aa_quad_evol(depth_in),
213  met_gamt(metgamt), gamt_point(gamt_point_in),
214  trK(trK_in), trK_point(trK_point_in), decouple(lapse_in.get_mp()){
215 
216  // hh_evol, trk_evol
217  hh_evol.update(met_gamt.con() - ff.con(), jtime, the_time[jtime]) ;
218  trk_evol.update(trK, jtime, the_time[jtime]) ;
219 
220 }
221 
222 // Copy constructor
223 // ----------------
224 
225 Isol_hor::Isol_hor(const Isol_hor& isolhor_in)
226  : Time_slice_conf(isolhor_in),
227  mp(isolhor_in.mp),
228  nz(isolhor_in.nz),
229  radius(isolhor_in.radius),
230  omega(isolhor_in.omega),
231  boost_x(isolhor_in.boost_x),
232  boost_z(isolhor_in.boost_z),
233  regul(isolhor_in.regul),
234  n_auto_evol(isolhor_in.n_auto_evol),
235  n_comp_evol(isolhor_in.n_comp_evol),
236  psi_auto_evol(isolhor_in.psi_auto_evol),
237  psi_comp_evol(isolhor_in.psi_comp_evol),
238  dn_evol(isolhor_in.dn_evol),
239  dpsi_evol(isolhor_in.dpsi_evol),
240  beta_auto_evol(isolhor_in.beta_auto_evol),
241  beta_comp_evol(isolhor_in.beta_comp_evol),
242  aa_auto_evol(isolhor_in.aa_auto_evol),
243  aa_comp_evol(isolhor_in.aa_comp_evol),
244  aa_nn(isolhor_in.aa_nn),
245  aa_quad_evol(isolhor_in.aa_quad_evol),
246  met_gamt(isolhor_in.met_gamt),
247  gamt_point(isolhor_in.gamt_point),
248  trK(isolhor_in.trK),
249  trK_point(isolhor_in.trK_point),
250  decouple(isolhor_in.decouple){
251 }
252 
253 // Constructor from a file
254 // -----------------------
255 
256 Isol_hor::Isol_hor(Map_af& mpi, FILE* fich,
257  bool partial_read, int depth_in)
258  : Time_slice_conf(mpi, mpi.get_bvect_spher(), mpi.flat_met_spher(),
259  fich, partial_read, depth_in),
260  mp(mpi), nz(mpi.get_mg()->get_nzone()), radius ((mpi.get_alpha())[0]),
261  omega(0), boost_x(0), boost_z(0), regul(0),
262  n_auto_evol(depth_in), n_comp_evol(depth_in),
263  psi_auto_evol(depth_in), psi_comp_evol(depth_in),
264  dn_evol(depth_in), dpsi_evol(depth_in),
265  beta_auto_evol(depth_in), beta_comp_evol(depth_in),
266  aa_auto_evol(depth_in), aa_comp_evol(depth_in),
267  aa_nn(depth_in), aa_quad_evol(depth_in),
268  met_gamt(mpi.flat_met_spher()),
269  gamt_point(mpi, CON, mpi.get_bvect_spher()),
270  trK(mpi), trK_point(mpi), decouple(mpi){
271 
272  fread_be(&omega, sizeof(double), 1, fich) ;
273  fread_be(&boost_x, sizeof(double), 1, fich) ;
274  fread_be(&boost_z, sizeof(double), 1, fich) ;
275 
276  int jmin = jtime - depth + 1 ;
277  int indicator ;
278 
279  // psi_evol
280  for (int j=jmin; j<=jtime; j++) {
281  fread_be(&indicator, sizeof(int), 1, fich) ;
282  if (indicator == 1) {
283  Scalar psi_file(mp, *(mp.get_mg()), fich) ;
284  psi_evol.update(psi_file, j, the_time[j]) ;
285  }
286  }
287 
288  // n_auto_evol
289  for (int j=jmin; j<=jtime; j++) {
290  fread_be(&indicator, sizeof(int), 1, fich) ;
291  if (indicator == 1) {
292  Scalar nn_auto_file(mp, *(mp.get_mg()), fich) ;
293  n_auto_evol.update(nn_auto_file, j, the_time[j]) ;
294  }
295  }
296 
297  // psi_auto_evol
298  for (int j=jmin; j<=jtime; j++) {
299  fread_be(&indicator, sizeof(int), 1, fich) ;
300  if (indicator == 1) {
301  Scalar psi_auto_file(mp, *(mp.get_mg()), fich) ;
302  psi_auto_evol.update(psi_auto_file, j, the_time[j]) ;
303  }
304  }
305 
306  // beta_auto_evol
307  for (int j=jmin; j<=jtime; j++) {
308  fread_be(&indicator, sizeof(int), 1, fich) ;
309  if (indicator == 1) {
310  Vector beta_auto_file(mp, mpi.get_bvect_spher(), fich) ;
311  beta_auto_evol.update(beta_auto_file, j, the_time[j]) ;
312  }
313  }
314 
315  // met_gamt, gamt_point, trK, trK_point
316 
317  Sym_tensor met_file (mp, mp.get_bvect_spher(), fich) ;
318  met_gamt = met_file ;
319 
320  Sym_tensor gamt_point_file (mp, mp.get_bvect_spher(), fich) ;
321  gamt_point = gamt_point_file ;
322 
323  Scalar trK_file (mp, *(mp.get_mg()), fich) ;
324  trK = trK_file ;
325 
326  Scalar trK_point_file (mp, *(mp.get_mg()), fich) ;
327  trK_point = trK_point_file ;
328 
329  // hh_evol, trk_evol
330  hh_evol.update(met_gamt.con() - ff.con(), jtime, the_time[jtime]) ;
331  trk_evol.update(trK, jtime, the_time[jtime]) ;
332 
333 }
334 
335  //--------------//
336  // Destructor //
337  //--------------//
338 
340 
341 
342  //-----------------------//
343  // Mutators / assignment //
344  //-----------------------//
345 
346 void Isol_hor::operator=(const Isol_hor& isolhor_in) {
347 
348 Time_slice_conf::operator=(isolhor_in) ;
349  mp = isolhor_in.mp ;
350  nz = isolhor_in.nz ;
351  radius = isolhor_in.radius ;
352  omega = isolhor_in.omega ;
353  boost_x = isolhor_in.boost_x ;
354  boost_z = isolhor_in.boost_z ;
355  regul = isolhor_in.regul ;
356  n_auto_evol = isolhor_in.n_auto_evol ;
357  n_comp_evol = isolhor_in.n_comp_evol ;
358  psi_auto_evol = isolhor_in.psi_auto_evol ;
359  psi_comp_evol = isolhor_in.psi_comp_evol ;
360  dn_evol = isolhor_in.dn_evol ;
361  dpsi_evol = isolhor_in.dpsi_evol ;
362  beta_auto_evol = isolhor_in.beta_auto_evol ;
363  beta_comp_evol = isolhor_in.beta_comp_evol ;
364  aa_auto_evol = isolhor_in.aa_auto_evol ;
365  aa_comp_evol = isolhor_in.aa_comp_evol ;
366  aa_nn = isolhor_in.aa_nn ;
367  aa_quad_evol = isolhor_in.aa_quad_evol ;
368  met_gamt = isolhor_in.met_gamt ;
369  gamt_point = isolhor_in.gamt_point ;
370  trK = isolhor_in.trK ;
371  trK_point = isolhor_in.trK_point ;
372  decouple = isolhor_in.decouple ;
373 }
374 
375 
376  //------------------//
377  // output //
378  //------------------//
379 
380 
381 ostream& Isol_hor::operator>>(ostream& flux) const {
382 
384 
385  flux << '\n' << "radius of the horizon : " << radius << '\n' ;
386  flux << "boost in x-direction : " << boost_x << '\n' ;
387  flux << "boost in z-direction : " << boost_z << '\n' ;
388  flux << "angular velocity omega : " << omega_hor() << '\n' ;
389  flux << "area of the horizon : " << area_hor() << '\n' ;
390  flux << "ang. mom. of horizon : " << ang_mom_hor() << '\n' ;
391  flux << "ADM ang. mom. : " << ang_mom_adm() << '\n' ;
392  flux << "Mass of the horizon : " << mass_hor() << '\n' ;
393  flux << "ADM Mass : " << adm_mass() << '\n' ;
394 
395  return flux ;
396 }
397 
398 
399  //--------------------------//
400  // Save in a file //
401  //--------------------------//
402 
403 
404 void Isol_hor::sauve(FILE* fich, bool partial_save) const {
405 
406 
407  // Writing of quantities common to all derived classes of Time_slice
408  // -----------------------------------------------------------------
409 
410  Time_slice_conf::sauve(fich, partial_save) ;
411 
412  fwrite_be (&omega, sizeof(double), 1, fich) ;
413  fwrite_be (&boost_x, sizeof(double), 1, fich) ;
414  fwrite_be (&boost_z, sizeof(double), 1, fich) ;
415 
416  // Writing of quantities common to Isol_hor
417  // -----------------------------------------
418 
419  int jmin = jtime - depth + 1 ;
420 
421  // psi_evol
422  for (int j=jmin; j<=jtime; j++) {
423  int indicator = (psi_evol.is_known(j)) ? 1 : 0 ;
424  fwrite_be(&indicator, sizeof(int), 1, fich) ;
425  if (indicator == 1) psi_evol[j].sauve(fich) ;
426  }
427 
428  // n_auto_evol
429  for (int j=jmin; j<=jtime; j++) {
430  int indicator = (n_auto_evol.is_known(j)) ? 1 : 0 ;
431  fwrite_be(&indicator, sizeof(int), 1, fich) ;
432  if (indicator == 1) n_auto_evol[j].sauve(fich) ;
433  }
434 
435  // psi_auto_evol
436  for (int j=jmin; j<=jtime; j++) {
437  int indicator = (psi_auto_evol.is_known(j)) ? 1 : 0 ;
438  fwrite_be(&indicator, sizeof(int), 1, fich) ;
439  if (indicator == 1) psi_auto_evol[j].sauve(fich) ;
440  }
441 
442  // beta_auto_evol
443  for (int j=jmin; j<=jtime; j++) {
444  int indicator = (beta_auto_evol.is_known(j)) ? 1 : 0 ;
445  fwrite_be(&indicator, sizeof(int), 1, fich) ;
446  if (indicator == 1) beta_auto_evol[j].sauve(fich) ;
447  }
448 
449  met_gamt.con().sauve(fich) ;
450  gamt_point.sauve(fich) ;
451  trK.sauve(fich) ;
452  trK_point.sauve(fich) ;
453 }
454 
455 // Accessors
456 // ---------
457 
458 const Scalar& Isol_hor::n_auto() const {
459 
460  assert( n_auto_evol.is_known(jtime) ) ;
461  return n_auto_evol[jtime] ;
462 }
463 
464 const Scalar& Isol_hor::n_comp() const {
465 
466  assert( n_comp_evol.is_known(jtime) ) ;
467  return n_comp_evol[jtime] ;
468 }
469 
470 const Scalar& Isol_hor::psi_auto() const {
471 
472  assert( psi_auto_evol.is_known(jtime) ) ;
473  return psi_auto_evol[jtime] ;
474 }
475 
476 const Scalar& Isol_hor::psi_comp() const {
477 
478  assert( psi_comp_evol.is_known(jtime) ) ;
479  return psi_comp_evol[jtime] ;
480 }
481 
482 const Vector& Isol_hor::dnn() const {
483 
484  assert( dn_evol.is_known(jtime) ) ;
485  return dn_evol[jtime] ;
486 }
487 
488 const Vector& Isol_hor::dpsi() const {
489 
490  assert( dpsi_evol.is_known(jtime) ) ;
491  return dpsi_evol[jtime] ;
492 }
493 
494 const Vector& Isol_hor::beta_auto() const {
495 
496  assert( beta_auto_evol.is_known(jtime) ) ;
497  return beta_auto_evol[jtime] ;
498 }
499 
500 const Vector& Isol_hor::beta_comp() const {
501 
502  assert( beta_comp_evol.is_known(jtime) ) ;
503  return beta_comp_evol[jtime] ;
504 }
505 
507 
508  assert( aa_auto_evol.is_known(jtime) ) ;
509  return aa_auto_evol[jtime] ;
510 }
511 
513 
514  assert( aa_comp_evol.is_known(jtime) ) ;
515  return aa_comp_evol[jtime] ;
516 }
517 
518 void Isol_hor::set_psi(const Scalar& psi_in) {
519 
520  psi_evol.update(psi_in, jtime, the_time[jtime]) ;
521 
522  // Reset of quantities depending on Psi:
523  if (p_psi4 != 0x0) {
524  delete p_psi4 ;
525  p_psi4 = 0x0 ;
526  }
527  if (p_ln_psi != 0x0) {
528  delete p_ln_psi ;
529  p_ln_psi = 0x0 ;
530  }
531  if (p_gamma != 0x0) {
532  delete p_gamma ;
533  p_gamma = 0x0 ;
534  }
535  gam_dd_evol.downdate(jtime) ;
536  gam_uu_evol.downdate(jtime) ;
537  k_dd_evol.downdate(jtime) ;
538  k_uu_evol.downdate(jtime) ;
539  adm_mass_evol.downdate(jtime) ;
540 
541 }
542 
543 void Isol_hor::set_nn(const Scalar& nn_in) {
544 
545  n_evol.update(nn_in, jtime, the_time[jtime]) ;
546 
547  hata_evol.downdate(jtime) ;
548  aa_quad_evol.downdate(jtime) ;
549  k_dd_evol.downdate(jtime) ;
550  k_uu_evol.downdate(jtime) ;
551 }
552 
553 void Isol_hor::set_gamt(const Metric& gam_tilde) {
554 
555  if (p_tgamma != 0x0) {
556  delete p_tgamma ;
557  p_tgamma = 0x0 ;
558  }
559  if (p_gamma != 0x0) {
560  delete p_gamma ;
561  p_gamma = 0x0 ;
562  }
563 
564  met_gamt = gam_tilde ;
565 
566  gam_dd_evol.downdate(jtime) ;
567  gam_uu_evol.downdate(jtime) ;
568  k_dd_evol.downdate(jtime) ;
569  k_uu_evol.downdate(jtime) ;
570  hh_evol.downdate(jtime) ;
571 
572  hh_evol.update(gam_tilde.con() - ff.con(), jtime, the_time[jtime]) ;
573 
574 }
575 
576 
577 
578 // Import the lapse from the companion (Bhole case)
579 
580 void Isol_hor::n_comp(const Isol_hor& comp) {
581 
582  double ttime = the_time[jtime] ;
583 
584  Scalar temp (mp) ;
585  temp.import(comp.n_auto()) ;
586  temp.std_spectral_base() ;
587  n_comp_evol.update(temp, jtime, ttime) ;
588  n_evol.update(temp + n_auto(), jtime, ttime) ;
589 
590  Vector dn_comp (mp, COV, mp.get_bvect_cart()) ;
591  dn_comp.set_etat_qcq() ;
592  Vector auxi (comp.n_auto().derive_cov(comp.ff)) ;
593  auxi.dec_dzpuis(2) ;
594  auxi.change_triad(auxi.get_mp().get_bvect_cart()) ;
595  for (int i=1 ; i<=3 ; i++){
596  if (auxi(i).get_etat() != ETATZERO)
597  auxi.set(i).raccord(3) ;
598  }
599 
600  auxi.change_triad(mp.get_bvect_cart()) ;
601  assert ( *(auxi.get_triad()) == *(dn_comp.get_triad())) ;
602 
603  for (int i=1 ; i<=3 ; i++){
604  dn_comp.set(i).import(auxi(i)) ;
605  dn_comp.set(i).set_spectral_va().set_base(auxi(i).get_spectral_va().
606  get_base()) ;
607  }
608  dn_comp.inc_dzpuis(2) ;
609  dn_comp.change_triad(mp.get_bvect_spher()) ;
610  /*
611  Vector dn_comp_zec (n_comp().derive_cov(ff)) ;
612  for (int i=1 ; i<=3 ; i++)
613  for (int l=nz-1 ; l<=nz-1 ; l++) {
614  if (dn_comp.set(i).get_etat() == ETATQCQ)
615  dn_comp.set(i).set_domain(l) = dn_comp_zec(i).domain(l) ;
616  }
617  */
618  dn_evol.update(n_auto().derive_cov(ff) + dn_comp, jtime, ttime) ;
619 
620 
621  /*
622  Scalar tr_K (mp) ;
623 
624  Mtbl mxabs (mp.xa) ;
625  Mtbl myabs (mp.ya) ;
626  Mtbl mzabs (mp.za) ;
627  Scalar comp_r (mp) ;
628  comp_r.annule_hard() ;
629  for (int l=1 ; l<mp.get_mg()->get_nzone() ; l++) {
630  int nr = mp.get_mg()->get_nr (l) ;
631  if (l==mp.get_mg()->get_nzone()-1)
632  nr -- ;
633  int np = mp.get_mg()->get_np (l) ;
634  int nt = mp.get_mg()->get_nt (l) ;
635  double xabs, yabs, zabs, air, theta, phi ;
636 
637  for (int k=0 ; k<np ; k++)
638  for (int j=0 ; j<nt ; j++)
639  for (int i=0 ; i<nr ; i++) {
640 
641  xabs = mxabs (l, k, j, i) ;
642  yabs = myabs (l, k, j, i) ;
643  zabs = mzabs (l, k, j, i) ;
644 
645  // coordinates of the point in 2 :
646  comp.mp.convert_absolute
647  (xabs, yabs, zabs, air, theta, phi) ;
648  comp_r.set_grid_point(l,k,j,i) = air ;
649  }
650  }
651 
652  Scalar fact(mp) ;
653  fact = 0.0000000000000001 ;
654 
655 // Scalar fact(mp) ;
656 // fact = 0.7*gam().radial_vect().divergence(gam()) ;
657 // fact.dec_dzpuis(2) ;
658 
659  tr_K = 1/mp.r/mp.r ;
660  tr_K = tr_K * fact ;
661  tr_K += fact/comp_r/comp_r ;
662  tr_K.std_spectral_base() ;
663  tr_K.annule(0, 0) ;
664  tr_K.raccord(1) ;
665  tr_K.inc_dzpuis(2) ;
666  trk_evol.update(tr_K, jtime, the_time[jtime]) ;
667 */
668 }
669 
670 // Import the conformal factor from the companion (Bhole case)
671 
672 void Isol_hor::psi_comp (const Isol_hor& comp) {
673 
674  double ttime = the_time[jtime] ;
675 
676  Scalar temp (mp) ;
677  temp.import(comp.psi_auto()) ;
678  temp.std_spectral_base() ;
679  psi_comp_evol.update(temp, jtime, ttime) ;
680  psi_evol.update(temp + psi_auto(), jtime, ttime) ;
681 
682  Vector dpsi_comp (mp, COV, mp.get_bvect_cart()) ;
683  dpsi_comp.set_etat_qcq() ;
684  Vector auxi (comp.psi_auto().derive_cov(comp.ff)) ;
685  auxi.dec_dzpuis(2) ;
686  auxi.change_triad(auxi.get_mp().get_bvect_cart()) ;
687  for (int i=1 ; i<=3 ; i++){
688  if (auxi(i).get_etat() != ETATZERO)
689  auxi.set(i).raccord(3) ;
690  }
691 
692  auxi.change_triad(mp.get_bvect_cart()) ;
693  assert ( *(auxi.get_triad()) == *(dpsi_comp.get_triad())) ;
694 
695  for (int i=1 ; i<=3 ; i++){
696  dpsi_comp.set(i).import(auxi(i)) ;
697  dpsi_comp.set(i).set_spectral_va().set_base(auxi(i).get_spectral_va().
698  get_base()) ;
699  }
700  dpsi_comp.inc_dzpuis(2) ;
701  dpsi_comp.change_triad(mp.get_bvect_spher()) ;
702  /*
703  Vector dpsi_comp_zec (psi_comp().derive_cov(ff)) ;
704  for (int i=1 ; i<=3 ; i++)
705  for (int l=nz-1 ; l<=nz-1 ; l++) {
706  if (dpsi_comp.set(i).get_etat() == ETATQCQ)
707  dpsi_comp.set(i).set_domain(l) = dpsi_comp_zec(i).domain(l) ;
708  }
709  */
710 
711  dpsi_evol.update(psi_auto().derive_cov(ff) + dpsi_comp, jtime, ttime) ;
712 
713 }
714 
715 void Isol_hor::beta_comp (const Isol_hor& comp) {
716 
717  double ttime = the_time[jtime] ;
718 
719  Vector tmp_vect (mp, CON, mp.get_bvect_cart()) ;
720  Vector shift_comp (comp.beta_auto()) ;
721  shift_comp.change_triad(comp.mp.get_bvect_cart()) ;
722  shift_comp.change_triad(mp.get_bvect_cart()) ;
723  assert (*(shift_comp.get_triad()) == *(tmp_vect.get_triad())) ;
724 
725  tmp_vect.set(1).import(shift_comp(1)) ;
726  tmp_vect.set(2).import(shift_comp(2)) ;
727  tmp_vect.set(3).import(shift_comp(3)) ;
728  tmp_vect.std_spectral_base() ;
729  tmp_vect.change_triad(mp.get_bvect_spher()) ;
730 
731  beta_comp_evol.update(tmp_vect, jtime,ttime) ;
732  beta_evol.update(beta_auto() + beta_comp(), jtime, ttime) ;
733 }
734 
735 //Initialisation to Schwartzchild
737 
738  double ttime = the_time[jtime] ;
739  Scalar auxi(mp) ;
740 
741  // Initialisation of the lapse different of zero on the horizon
742  // at the first step
743  auxi = 0.5 - 0.5/mp.r ;
744  auxi.annule(0, 0);
745  auxi.set_dzpuis(0) ;
746 
747  Scalar temp(mp) ;
748  temp = auxi;
749  temp.std_spectral_base() ;
750  temp.raccord(1) ;
751  n_auto_evol.update(temp, jtime, ttime) ;
752 
753  temp = 0.5 ;
754  temp.std_spectral_base() ;
755  n_comp_evol.update(temp, jtime, ttime) ;
756  n_evol.update(n_auto() + n_comp(), jtime, ttime) ;
757 
758  auxi = 0.5 + radius/mp.r ;
759  auxi.annule(0, 0);
760  auxi.set_dzpuis(0) ;
761  temp = auxi;
762  temp.std_spectral_base() ;
763  temp.raccord(1) ;
764  psi_auto_evol.update(temp, jtime, ttime) ;
765 
766  temp = 0.5 ;
767  temp.std_spectral_base() ;
768  psi_comp_evol.update(temp, jtime, ttime) ;
769  psi_evol.update(psi_auto() + psi_comp(), jtime, ttime) ;
770 
771  dn_evol.update(nn().derive_cov(ff), jtime, ttime) ;
772  dpsi_evol.update(psi().derive_cov(ff), jtime, ttime) ;
773 
774  Vector temp_vect1(mp, CON, mp.get_bvect_spher()) ;
775  temp_vect1.set(1) = 0.0/mp.r/mp.r ;
776  temp_vect1.set(2) = 0. ;
777  temp_vect1.set(3) = 0. ;
778  temp_vect1.std_spectral_base() ;
779 
780  Vector temp_vect2(mp, CON, mp.get_bvect_spher()) ;
781  temp_vect2.set_etat_zero() ;
782 
783  beta_auto_evol.update(temp_vect1, jtime, ttime) ;
784  beta_comp_evol.update(temp_vect2, jtime, ttime) ;
785  beta_evol.update(temp_vect1, jtime, ttime) ;
786 }
787 
789 
790  Metric flat (mp.flat_met_spher()) ;
791  met_gamt = flat ;
792 
794  trK.set_etat_zero() ;
796 
797 }
798 
799 
801 
802  double ttime = the_time[jtime] ;
803  Scalar auxi(mp) ;
804 
805  auxi = (1-radius/mp.r)/(1+radius/mp.r) ;
806  auxi.annule(0, 0);
807  auxi.set_outer_boundary((*mp.get_mg()).get_nzone(), 1.) ;
808  auxi.set_dzpuis(0) ;
809 
810  Scalar temp(mp) ;
811  temp = auxi;
812  temp.std_spectral_base() ;
813  temp.raccord(1) ;
814  n_auto_evol.update(temp, jtime, ttime) ;
815 
816  temp.set_etat_zero() ;
817  n_comp_evol.update(temp, jtime, ttime) ;
818  n_evol.update(temp, jtime, ttime) ;
819 
820 
821  auxi = 1 + radius/mp.r ;
822  auxi.annule(0, 0);
823  auxi.set_outer_boundary((*mp.get_mg()).get_nzone(), 1.) ;
824  auxi.set_dzpuis(0) ;
825 
826  temp = auxi;
827  temp.std_spectral_base() ;
828  temp.raccord(1) ;
829  psi_auto_evol.update(temp, jtime, ttime) ;
830  temp.set_etat_zero() ;
831  psi_comp_evol.update(temp, jtime, ttime) ;
832  psi_evol.update(temp, jtime, ttime) ;
833 
834  dn_evol.update(nn().derive_cov(ff), jtime, ttime) ;
835  dpsi_evol.update(psi().derive_cov(ff), jtime, ttime) ;
836 
837  Vector temp_vect(mp, CON, mp.get_bvect_spher()) ;
838  temp_vect.set_etat_zero() ;
839  beta_auto_evol.update(temp_vect, jtime, ttime) ;
840  beta_comp_evol.update(temp_vect, jtime, ttime) ;
841  beta_evol.update(temp_vect, jtime, ttime) ;
842 
843 }
844 
846 
847  Sym_tensor aa_new (mp, CON, mp.get_bvect_spher()) ;
848  int nnt = mp.get_mg()->get_nt(1) ;
849  int nnp = mp.get_mg()->get_np(1) ;
850 
851  int check ;
852  check = 0 ;
853  for (int k=0; k<nnp; k++)
854  for (int j=0; j<nnt; j++){
855  if (nn().val_grid_point(1, k, j , 0) < 1e-12){
856  check = 1 ;
857  break ;
858  }
859  }
860 
861  if (check == 0)
862  aa_new = ( beta().ope_killing_conf(met_gamt) + gamt_point )
863  / (2.* nn()) ;
864  else {
865  regul = regularise_one() ;
866  cout << "regul = " << regul << endl ;
867  Scalar nn_sxpun (division_xpun (Cmp(nn()), 0)) ;
868  aa_new = beta().ope_killing_conf(met_gamt) + gamt_point ;
869 
870  Scalar auxi (mp) ;
871  for (int i=1 ; i<=3 ; i++)
872  for (int j=i ; j<=3 ; j++) {
873  auxi = aa_new(i, j) ;
874  auxi = division_xpun (auxi, 0) ;
875  aa_new.set(i,j) = auxi / nn_sxpun / 2. ;
876  }
877  }
878 
879  Sym_tensor hata_new = aa_new*psi4()*psi()*psi() ;
880  set_hata(hata_new) ; // set aa to aa_new and delete values of
881  // k_uu and k_dd.
882  Sym_tensor aa_dd (aa_new.up_down(met_gamt)) ;
883  Scalar aquad (contract(aa_dd, 0, 1, aa_new, 0, 1)*psi4()*psi4()*psi4()) ;
884  aa_quad_evol.update(aquad, jtime, the_time[jtime]) ;
885 }
886 
887 const Scalar& Isol_hor::aa_quad() const {
888 
889  if (!aa_quad_evol.is_known(jtime) ) {
890  Sym_tensor aa_dd (aa().up_down(met_gamt)) ;
891  Scalar aquad (contract(aa_dd, 0, 1, aa(), 0, 1)*psi4()*psi4()*psi4()) ;
892  aa_quad_evol.update(aquad, jtime, the_time[jtime]) ;
893  }
894 
895  return aa_quad_evol[jtime] ;
896 }
897 
899 
900  Sym_tensor gamm (gam().cov()) ;
901  Sym_tensor gamt (gamm / gamm(3,3)) ;
902 
903  Metric metgamt (gamt) ;
904  met_gamt = metgamt ;
905 
906  Scalar psi_perturb (pow(gamm(3,3), 0.25)) ;
907  psi_perturb.std_spectral_base() ;
908  set_psi(psi_perturb) ;
909 
910  cout << "met_gamt" << endl << norme(met_gamt.cov()(1,1)) << endl
911  << norme(met_gamt.cov()(2,1)) << endl << norme(met_gamt.cov()(3,1))
912  << endl << norme(met_gamt.cov()(2,2)) << endl
913  << norme(met_gamt.cov()(3,2)) << endl << norme(met_gamt.cov()(3,3))
914  << endl ;
915  cout << "determinant" << norme(met_gamt.determinant()) << endl ;
916 
917  hh_evol.update(met_gamt.con() - ff.con(), jtime, the_time[jtime]) ;
918 }
919 
920 
921 void Isol_hor::aa_kerr_ww(double mm, double aaa) {
922 
923  Scalar rr(mp) ;
924  rr = mp.r ;
925  Scalar cost (mp) ;
926  cost = mp.cost ;
927  Scalar sint (mp) ;
928  sint = mp.sint ;
929 
930  // rbl
931  Scalar rbl = rr + mm + (mm*mm - aaa*aaa) / (4*rr) ;
932 
933  // sigma inverse
934  Scalar sigma_inv = 1. / (rbl*rbl + aaa*aaa*cost*cost) ;
935  sigma_inv.set_domain(0) = 1. ;
936 
937  // ww perturbation
938  Scalar ww_pert (mp) ;
939  ww_pert = - 1*(mm*aaa*aaa*aaa*pow(sint, 4.)*cost) * sigma_inv ;
940  ww_pert.set_spectral_va().set_base_r(0,R_CHEBPIM_P) ;
941  for (int l=1; l<nz-1; l++)
942  ww_pert.set_spectral_va().set_base_r(l,R_CHEB) ;
943  ww_pert.set_spectral_va().set_base_r(nz-1,R_CHEBU) ;
944  ww_pert.set_spectral_va().set_base_t(T_COSSIN_CI) ;
945  ww_pert.set_spectral_va().set_base_p(P_COSSIN) ;
946 
947  // Quadratic part A^{ij]A_{ij}
948  // Lichnerowicz choice
949  //----------------------------
950 
951  // BY - BY
952  Vector dw_by (mp, COV, mp.get_bvect_spher()) ;
953  dw_by.set(1) = 0. ;
954  dw_by.set(2) = 3 * aaa*mm*sint*sint*sint / rr ;
955  dw_by.set(3) = 0. ;
956  dw_by.set(2).set_spectral_va().set_base_r(0,R_CHEBPIM_P) ;
957  for (int l=1; l<nz-1; l++)
958  dw_by.set(2).set_spectral_va().set_base_r(l,R_CHEB) ;
959  dw_by.set(2).set_spectral_va().set_base_r(nz-1,R_CHEBU) ;
960  dw_by.set(2).set_spectral_va().set_base_t(T_COSSIN_SI) ;
961  dw_by.set(2).set_spectral_va().set_base_p(P_COSSIN) ;
962 
963  Scalar aquad_1 = 2*contract(dw_by, 0, dw_by.up_down(ff), 0) *
964  gam_dd()(3,3) / gam_dd()(1,1) ;
965  aquad_1.div_rsint_dzpuis(1) ;
966  aquad_1.div_rsint_dzpuis(2) ;
967  aquad_1.div_rsint_dzpuis(3) ;
968  aquad_1.div_rsint_dzpuis(4) ;
969 
970  // BY - dw_pert
971  Vector dw_pert(ww_pert.derive_con(ff)) ;
972  Scalar aquad_2 = 4*contract(dw_by, 0, dw_pert, 0) *
973  gam_dd()(3,3) / gam_dd()(1,1) ;
974 
975  aquad_2.div_rsint_dzpuis(3) ;
976  aquad_2.div_rsint_dzpuis(4) ;
977  aquad_2.div_rsint() ;
978  aquad_2.div_rsint() ;
979 
980  // dw_pert - dw_pert
981  Scalar aquad_3 = 2*contract(dw_pert, 0, dw_pert.up_down(ff), 0) *
982  gam_dd()(3,3) / gam_dd()(1,1) ;
983 
984  aquad_3.div_rsint() ;
985  aquad_3.div_rsint() ;
986  aquad_3.div_rsint() ;
987  aquad_3.div_rsint() ;
988 
989  // Total
990  Scalar aquad = aquad_1 + aquad_2 + aquad_3 ;
991 
992  aquad.set_domain(0) = 0. ;
993  Base_val sauve_base (aquad.get_spectral_va().get_base()) ;
994 
995  aquad = aquad * pow(gam_dd()(1,1), 2.) * pow(gam_dd()(3,3), -2.) ;
996  aquad.set_spectral_va().set_base(sauve_base) ;
997 
998 /*
999  cout << "norme de aquad" << endl << norme(aquad) << endl ;
1000  cout << "norme de aa_quad" << endl << norme(aa_quad()) << endl ;
1001 
1002  des_meridian (aquad, 0, 4, "aquad", 1) ;
1003  des_meridian (aa_quad(), 0, 4, "aa_quad()", 2) ;
1004  des_meridian (aa_quad()-aquad, 0, 4, "diff aa_quad", 3) ;
1005  arrete() ;
1006 */
1007 
1008  aa_quad_evol.update(aquad, jtime, the_time[jtime]) ;
1009 
1010 
1011  // Extrinsic curvature A^{ij} and A_{ij}
1012  // Dynamical choice
1013  // -------------------------------------
1014 
1015  Scalar s_r (ww_pert.derive_cov(ff)(2)) ;
1016  s_r = - s_r * gam().cov()(3,3) / gam().cov()(1,1) ;
1017  s_r.div_rsint() ;
1018 
1019  Scalar temp = dw_by(2) ;
1020  temp = - temp * gam().cov()(3,3) / gam().cov()(1,1) ;
1021  temp.div_rsint_dzpuis(2) ;
1022 
1023  s_r = s_r + temp ;
1024  s_r.annule_domain(0) ;
1025 
1026  Scalar s_t (ww_pert.derive_cov(ff)(1)) ;
1027  s_t = s_t * gam().cov()(3,3) / gam().cov()(1,1) ;
1028  s_t.div_rsint() ;
1029 
1030  temp = dw_by(1) ;
1031  temp = temp * gam().cov()(3,3) / gam().cov()(1,1) ;
1032  temp.div_rsint_dzpuis(2) ;
1033 
1034  s_t = s_t + temp ;
1035  s_t.annule_domain(0) ;
1036 
1037 
1038  Vector ss (mp, CON, mp.get_bvect_spher()) ;
1039  ss.set(1) = s_r ;
1040  ss.set(2) = s_t ;
1041  ss.set(3) = 0. ;
1042 
1043  Sym_tensor aij (mp, CON, mp.get_bvect_spher()) ;
1044  aij.set(1,1) = 0. ;
1045  aij.set(2,1) = 0. ;
1046  aij.set(2,2) = 0. ;
1047  aij.set(3,3) = 0. ;
1048  aij.set(3,1) = s_r ;
1049  aij.set(3,1).div_rsint() ;
1050  aij.set(3,2) = s_t ;
1051  aij.set(3,2).div_rsint() ;
1052 
1053  Base_val base_31 (aij(3,1).get_spectral_va().get_base()) ;
1054  Base_val base_32 (aij(3,2).get_spectral_va().get_base()) ;
1055 
1056  aij.set(3,1) = aij(3,1) * pow(gam_dd()(1,1), 5./3.)
1057  * pow(gam_dd()(3,3), -5./3.) ;
1058  aij.set(3,1) = aij(3,1) * pow(psi(), -6.) ;
1059  aij.set(3,1).set_spectral_va().set_base(base_31) ;
1060  aij.set(3,2) = aij(3,2) * pow(gam_dd()(1,1), 5./3.)
1061  * pow(gam_dd()(3,3), -5./3.) ;
1062  aij.set(3,2) = aij(3,2) * pow(psi(), -6.) ;
1063  aij.set(3,2).set_spectral_va().set_base(base_32) ;
1064 
1065  /*
1066  cout << "norme de A(3,1)" << endl << norme(aij(3,1)) << endl ;
1067  cout << "norme de A(3,2)" << endl << norme(aij(3,2)) << endl ;
1068 
1069  cout << "norme de A_init(3,1)" << endl << norme(aa()(3,1)) << endl ;
1070  cout << "norme de A_init(3,2)" << endl << norme(aa()(3,2)) << endl ;
1071 
1072  des_meridian(aij(3,1), 0., 4., "aij(3,1)", 0) ;
1073  des_meridian(aa()(3,1), 0., 4., "aa_init(3,1)", 1) ;
1074  des_meridian(aa()(3,1)-aij(3,1), 0., 4., "diff_aa(3,1)", 2) ;
1075  des_meridian(aij(3,2), 0., 4., "aij(3,2)", 3) ;
1076  des_meridian(aa()(3,2), 0., 4., "aa_init(3,2)", 4) ;
1077  des_meridian(aa()(3,2)-aij(3,2), 0., 4., "diff_aa(3,2)", 5) ;
1078  arrete() ;
1079  */
1080  Sym_tensor hataij = aij*psi4()*psi()*psi() ;
1081  hata_evol.update(hataij, jtime, the_time[jtime]) ;
1082  Sym_tensor kij (aij) ;
1083  kij = kij * pow(gam().determinant(), -1./3.) ;
1084  kij.std_spectral_base() ;
1085  k_uu_evol.update(kij, jtime, the_time[jtime]) ;
1086  k_dd_evol.update(kij.up_down(gam()), jtime, the_time[jtime]) ;
1087 
1088 }
1089 
1090 double Isol_hor::axi_break() const {
1091 
1092  Vector phi (ff.get_mp(), CON, *(ff.get_triad()) ) ;
1093 
1094  Scalar tmp (ff.get_mp() ) ;
1095  tmp = 1 ;
1096  tmp.std_spectral_base() ;
1097  tmp.mult_rsint() ;
1098 
1099  phi.set(1) = 0. ;
1100  phi.set(2) = 0. ;
1101  phi.set(3) = tmp ;
1102 
1103  Sym_tensor q_uu ( gam_uu() - gam().radial_vect() * gam().radial_vect() ) ;
1104  Sym_tensor q_dd ( q_uu.up_down(gam()) ) ;
1105 
1106  Sym_tensor L_phi_q_dd ( q_dd.derive_lie( phi) ) ;
1107  Sym_tensor L_phi_q_uu ( contract(contract(L_phi_q_dd, 0, q_uu, 0), 0, q_uu,0) ) ;
1108 
1109 
1110  Scalar integrand ( contract( L_phi_q_dd, 0, 1, L_phi_q_uu, 0, 1 ) * darea_hor() ) ;
1111 
1112  double axibreak = mp.integrale_surface(integrand, 1.0000000001)/
1113  (4 * M_PI * radius_hor()* radius_hor() ) ;
1114 
1115  return axibreak ;
1116 
1117 }
1118 
1119 double fonc_expansion(double rr, const Param& par_expansion) {
1120 
1121  Scalar expa = par_expansion.get_scalar(0) ;
1122  double theta = par_expansion.get_double(0) ;
1123  double phi = par_expansion.get_double(1) ;
1124 
1125  return expa.val_point(rr, theta, phi) ;
1126 
1127 }
1128 void Isol_hor::adapt_hor(double c_min, double c_max) {
1129 
1130  Scalar expa (expansion()) ;
1131  Scalar app_hor(mp) ;
1132  app_hor.annule_hard() ;
1133  int nitmax = 200 ;
1134  int nit ;
1135 
1136  double precis = 1.e-13 ;
1137 
1138  // Calculation of the radius of the apparent horizon for each (theta, phi)
1139  // -----------------------------------------------------------------------
1140 
1141  for (int nt=0; nt<mp.get_mg()->get_nt(1); nt++)
1142  for (int np=0; np<mp.get_mg()->get_np(1); np++) {
1143 
1144  double theta = mp.get_mg()->get_grille3d(1)->tet[nt] ;
1145  double phi = mp.get_mg()->get_grille3d(1)->phi[np] ;
1146 
1147  Param par_expansion ;
1148  par_expansion.add_scalar(expa, 0) ;
1149  par_expansion.add_double(theta, 0) ;
1150  par_expansion.add_double(phi, 1) ;
1151  double r_app_hor = zerosec_b(fonc_expansion, par_expansion, c_min, c_max,
1152  precis, nitmax, nit) ;
1153 
1154  app_hor.set_grid_point(1, np, nt, 0) = r_app_hor ;
1155  }
1156 
1157  // Transformation of the 3-metric and extrinsic curvature
1158  // ------------------------------------------------------
1159 
1160  Scalar rr (mp) ;
1161  rr = mp.r ;
1162 
1163  Scalar trans_11 (mp) ;
1164  Scalar r_new (mp) ;
1165  r_new.annule_hard() ;
1166  // trans_11.annule_hard() ;
1167  for (int l=1; l<nz; l++)
1168  for (int nr=0; nr<mp.get_mg()->get_nr(1); nr++)
1169  for (int nt=0; nt<mp.get_mg()->get_nt(1); nt++)
1170  for (int np=0; np<mp.get_mg()->get_np(1); np++) {
1171  r_new.set_grid_point(l, np, nt, nr) = rr.val_grid_point(l, np, nt, nr) -
1172  app_hor.val_grid_point(1, np, nt, 0) + 1 ;
1173  // trans_11.set_grid_point(l, np, nt, nr) = 1. /
1174  // app_hor.val_grid_point(1, np, nt, 0) ; // !
1175  }
1176  r_new.std_spectral_base() ;
1177 
1178  Itbl comp(2) ;
1179  comp.set(0) = CON ;
1180  comp.set(1) = COV ;
1181 
1182  Scalar trans_12 (r_new.dsdt()) ;
1183  trans_12.div_r() ;
1184  Scalar trans_13 (r_new.stdsdp()) ;
1185  trans_13.div_r() ;
1186  for (int nr=0; nr<mp.get_mg()->get_nr(1); nr++)
1187  for (int nt=0; nt<mp.get_mg()->get_nt(1); nt++)
1188  for (int np=0; np<mp.get_mg()->get_np(1); np++) {
1189  trans_12.set_grid_point(nz-1, np, nt, nr) = trans_12.val_grid_point(1, np, nt, nr) ;
1190  trans_13.set_grid_point(nz-1, np, nt, nr) = trans_13.val_grid_point(1, np, nt, nr) ;
1191  }
1192 
1193  // Transformation matrix
1194  Tensor trans (mp, 2, comp, mp.get_bvect_spher()) ;
1195  trans.set(1,1) = 1 ;
1196  trans.set(1,2) = 0;//trans_12 ;
1197  trans.set(1,3) = 0;//trans_13 ;
1198  trans.set(2,2) = 1. ;
1199  trans.set(3,3) = 1. ;
1200  trans.set(2,1) = 0. ;
1201  trans.set(3,1) = 0. ;
1202  trans.set(3,2) = 0. ;
1203  trans.set(2,3) = 0. ;
1204  trans.std_spectral_base() ;
1205 
1206  cout << "trans(1,3)" << endl << norme(trans(1,3)) << endl ;
1207 
1208  Sym_tensor gamma_uu (gam().con()) ;
1209  Sym_tensor kk_uu (k_uu()) ;
1210 
1211  gamma_uu = contract(gamma_uu, 0, 1, trans * trans, 1, 3) ;
1212  kk_uu = contract(kk_uu, 0, 1, trans * trans, 1, 3) ;
1213 
1214  Sym_tensor copie_gamma (gamma_uu) ;
1215  Sym_tensor copie_kk (kk_uu) ;
1216 
1217  // dz_puis set to zero
1218  kk_uu.dec_dzpuis(2) ;
1219  for(int i=1; i<=3; i++)
1220  for(int j=i; j<=3; j++){
1221  kk_uu.set(i,j).annule_hard() ;
1222  gamma_uu.set(i,j).annule_hard() ;
1223  }
1224 
1225  copie_kk.dec_dzpuis(2) ;
1226 
1227  Scalar expa_trans(mp) ;
1228  expa_trans.annule_hard() ;
1229  expa.dec_dzpuis(2) ;
1230 
1231  /*
1232  copie_gamma.set(2,2).div_r() ;
1233  copie_gamma.set(2,2).div_r() ;
1234  copie_gamma.set(3,3).div_rsint() ;
1235  copie_gamma.set(3,3).div_rsint() ;
1236  copie_gamma.set(1,2).div_r() ;
1237  copie_gamma.set(1,3).div_rsint() ;
1238  // gamma_uu.set(2,3).div_r() ;
1239  // gamma_uu.set(2,3).div_rsint() ;
1240  */
1241 
1242  //Importation
1243  for(int i=1; i<=3; i++)
1244  for(int j=i; j<=3; j++)
1245  for (int l=1; l<nz; l++)
1246  for (int nr=0; nr<mp.get_mg()->get_nr(1); nr++)
1247  for (int nt=0; nt<mp.get_mg()->get_nt(1); nt++)
1248  for (int np=0; np<mp.get_mg()->get_np(1); np++) {
1249 
1250  double theta = mp.get_mg()->get_grille3d(1)->tet[nt] ;
1251  double phi = mp.get_mg()->get_grille3d(1)->phi[np] ;
1252  double r_inv = rr.val_grid_point(l, np, nt, nr) +
1253  app_hor.val_grid_point(1, np, nt, 0) - 1. ;
1254 
1255  gamma_uu.set(i,j).set_grid_point(l, np, nt, nr) =
1256  copie_gamma(i,j).val_point(r_inv, theta, phi) ;
1257  kk_uu.set(i,j).set_grid_point(l, np, nt, nr) =
1258  copie_kk(i,j).val_point(r_inv, theta, phi) ;
1259 
1260  expa_trans.set_grid_point(l, np, nt, nr) = expa.val_point(r_inv, theta, phi) ;
1261  }
1262  kk_uu.std_spectral_base() ; // Save the base?
1263  gamma_uu.std_spectral_base() ;
1264  expa_trans.std_spectral_base() ;
1265 
1266  for (int l=1; l<nz; l++)
1267  for (int nr=0; nr<mp.get_mg()->get_nr(1); nr++)
1268  for (int nt=0; nt<mp.get_mg()->get_nt(1); nt++)
1269  for (int np=0; np<mp.get_mg()->get_np(1); np++) {
1270  gamma_uu.set(1,2).set_grid_point(l,np,nt,nr) = gamma_uu.set(1,2).val_grid_point(l,np,nt,nr)
1271  / (1 +app_hor.val_grid_point(1, np, nt, 0)/ rr.val_grid_point(l, np, nt, nr)
1272  - 1 / rr.val_grid_point(l, np, nt, nr) ) ;
1273  gamma_uu.set(1,3).set_grid_point(l,np,nt,nr) = gamma_uu.set(1,3).val_grid_point(l,np,nt,nr)
1274  / (1 +app_hor.val_grid_point(1, np, nt, 0)/ rr.val_grid_point(l, np, nt, nr)
1275  - 1 / rr.val_grid_point(l, np, nt, nr) ) ;
1276  gamma_uu.set(2,2).set_grid_point(l,np,nt,nr) = gamma_uu.set(2,2).val_grid_point(l,np,nt,nr)
1277  / (1 +app_hor.val_grid_point(1, np, nt, 0)/ rr.val_grid_point(l, np, nt, nr)
1278  - 1 / rr.val_grid_point(l, np, nt, nr) )
1279  / (1 +app_hor.val_grid_point(1, np, nt, 0)/ rr.val_grid_point(l, np, nt, nr)
1280  - 1 / rr.val_grid_point(l, np, nt, nr) ) ;
1281  gamma_uu.set(2,3).set_grid_point(l,np,nt,nr) = gamma_uu.set(2,3).val_grid_point(l,np,nt,nr)
1282  / (1 +app_hor.val_grid_point(1, np, nt, 0)/ rr.val_grid_point(l, np, nt, nr)
1283  - 1 / rr.val_grid_point(l, np, nt, nr) )
1284  / (1 +app_hor.val_grid_point(1, np, nt, 0)/ rr.val_grid_point(l, np, nt, nr)
1285  - 1 / rr.val_grid_point(l, np, nt, nr) ) ;
1286  gamma_uu.set(3,3).set_grid_point(l,np,nt,nr) = gamma_uu.set(3,3).val_grid_point(l,np,nt,nr)
1287  / (1 +app_hor.val_grid_point(1, np, nt, 0)/ rr.val_grid_point(l, np, nt, nr)
1288  - 1 / rr.val_grid_point(l, np, nt, nr) )
1289  / (1 +app_hor.val_grid_point(1, np, nt, 0)/ rr.val_grid_point(l, np, nt, nr)
1290  - 1 / rr.val_grid_point(l, np, nt, nr) ) ;
1291  }
1292 
1293  /*
1294  gamma_uu.set(2,2).mult_r() ;
1295  gamma_uu.set(2,2).mult_r() ;
1296  gamma_uu.set(3,3).mult_rsint() ;
1297  gamma_uu.set(3,3).mult_rsint() ;
1298  gamma_uu.set(1,2).mult_r() ;
1299  gamma_uu.set(1,3).mult_rsint() ;
1300  // gamma_uu.set(2,3).mult_r() ;
1301  // gamma_uu.set(2,3).mult_rsint() ;
1302  */
1303 
1304 
1305 
1306  cout << "gamma_uu(1,1)" << endl << norme(gamma_uu(1,1)) << endl ;
1307  cout << "gamma_uu(2,1)" << endl << norme(gamma_uu(2,1)) << endl ;
1308  cout << "gamma_uu(3,1)" << endl << norme(gamma_uu(3,1)) << endl ;
1309  cout << "gamma_uu(2,2)" << endl << norme(gamma_uu(2,2)) << endl ;
1310  cout << "gamma_uu(3,2)" << endl << norme(gamma_uu(3,2)) << endl ;
1311  cout << "gamma_uu(3,3)" << endl << norme(gamma_uu(3,3)) << endl ;
1312 
1313  kk_uu.inc_dzpuis(2) ;
1314  expa_trans.inc_dzpuis(2) ;
1315 
1316  Metric gamm (gamma_uu) ;
1317 
1318  // Updates
1319  gam_uu_evol.update(gamma_uu, jtime, the_time[jtime]) ;
1320  gam_dd_evol.update(gamm.cov(), jtime, the_time[jtime]) ;
1321  k_uu_evol.update(kk_uu, jtime, the_time[jtime]) ;
1322 
1323  if (p_psi4 != 0x0) {
1324  delete p_psi4 ;
1325  p_psi4 = 0x0 ;
1326  }
1327  if (p_ln_psi != 0x0) {
1328  delete p_ln_psi ;
1329  p_ln_psi = 0x0 ;
1330  }
1331  if (p_gamma != 0x0) {
1332  delete p_gamma ;
1333  p_gamma = 0x0 ;
1334  }
1335  if (p_tgamma != 0x0) {
1336  delete p_tgamma ;
1337  p_tgamma = 0x0 ;
1338  }
1339  if (p_hdirac != 0x0) {
1340  delete p_hdirac ;
1341  p_hdirac = 0x0 ;
1342  }
1343 
1344  k_dd_evol.downdate(jtime) ;
1345  psi_evol.downdate(jtime) ;
1346  hata_evol.downdate(jtime) ;
1347  aa_quad_evol.downdate(jtime) ;
1348  beta_evol.downdate(jtime) ;
1349  n_evol.downdate(jtime) ;
1350  hh_evol.downdate(jtime) ;
1351 
1352 
1353  Scalar new_expa (expansion()) ;
1354  //des_meridian(expa_trans, 1., 6., "Expansion trans", 1) ;
1355  //des_meridian(new_expa, 1.000000001, 6., "Expansion new", 2) ;
1356  //des_meridian(expa_trans- new_expa, 1.000000001, 4., "diff Expansion trans", 3) ;
1357 
1358 
1359 
1360 }
1361 
1362 }
virtual const Vector & beta_comp() const
Shift function at the current time step jtime.
Definition: isol_hor.C:500
double area_hor() const
Area of the horizon.
Definition: phys_param.C:160
virtual const Vector & beta() const
shift vector at the current time step (jtime )
Evolution_std< Scalar > n_comp_evol
Values at successive time steps of the lapse function .
Definition: isol_hor.h:284
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:675
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
const Scalar & psi4() const
Factor at the current time step (jtime ).
Metric for tensor calculation.
Definition: metric.h:90
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for domain no. l.
Definition: grilles.h:517
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition: tensor.C:490
double regul
Intensity of the correction on the shift vector.
Definition: isol_hor.h:278
virtual void sauve(FILE *) const
Save in a binary file.
Definition: tensor_sym.C:375
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:293
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
double omega
Angular velocity in LORENE&#39;s units.
Definition: isol_hor.h:269
virtual const Sym_tensor & gam_uu() const
Induced metric (contravariant components ) at the current time step (jtime )
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
void operator=(const Time_slice_conf &)
Assignment to another Time_slice_conf.
virtual ~Isol_hor()
Destructor.
Definition: isol_hor.C:339
Evolution_std< Scalar > n_auto_evol
Values at successive time steps of the lapse function .
Definition: isol_hor.h:281
Evolution_std< Scalar > psi_comp_evol
Values at successive time steps of the lapse function .
Definition: isol_hor.h:290
Evolution_std< Vector > beta_auto_evol
Values at successive time steps of the shift function .
Definition: isol_hor.h:301
#define P_COSSIN
dev. standart
Definition: type_parite.h:245
Evolution_std< Scalar > psi_evol
Values at successive time steps of the conformal factor relating the physical metric to the conform...
Definition: time_slice.h:520
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
void init_bhole()
Sets the values of the fields to :
Definition: isol_hor.C:736
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:330
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition: scalar.C:397
Sym_tensor gamt_point
Time derivative of the 3-metric tilde.
Definition: isol_hor.h:329
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:801
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator<<).
Scalar * p_ln_psi
Pointer on the logarithm of at the current time step (jtime)
Definition: time_slice.h:569
Lorene prototypes.
Definition: app_hor.h:67
double ang_mom_hor() const
Angular momentum (modulo)
Definition: phys_param.C:181
double radius
Radius of the horizon in LORENE&#39;s units.
Definition: isol_hor.h:266
Evolution_full< Tbl > adm_mass_evol
ADM mass at each time step, since the creation of the slice.
Definition: time_slice.h:236
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
void update_aa()
Conformal representation of the traceless part of the extrinsic curvature: .
Definition: isol_hor.C:845
double * phi
Array of values of at the np collocation points.
Definition: grilles.h:219
Flat metric for tensor calculation.
Definition: metric.h:261
void operator=(const Isol_hor &)
Assignment to another Isol_hor.
Definition: isol_hor.C:346
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the metric is defined.
Definition: metric.h:309
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
double zerosec_b(double(*f)(double, const Param &), const Param &par, double a, double b, double precis, int nitermax, int &niter)
Finding the zero a function on a bounded domain.
Definition: zerosec_borne.C:71
int jtime
Time step index of the latest slice.
Definition: time_slice.h:193
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:790
Evolution_std< Sym_tensor > gam_dd_evol
Values at successive time steps of the covariant components of the induced metric ...
Definition: time_slice.h:201
void init_bhole_seul()
Initiates for a single black hole.
Definition: isol_hor.C:800
Scalar trK
Trace of the extrinsic curvature.
Definition: isol_hor.h:332
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
virtual const Scalar & psi_comp() const
Conformal factor at the current time step jtime.
Definition: isol_hor.C:476
void set(const Map *mp, Mtbl *(*construct)(const Map *))
Semi-constructor from a mapping and a method.
Definition: coord.C:137
Tensor field of valence 1.
Definition: vector.h:188
int nz
Number of zones.
Definition: isol_hor.h:263
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Definition: scalar.h:621
Evolution_std< Scalar > psi_auto_evol
Values at successive time steps of the conformal factor .
Definition: isol_hor.h:287
Vector * p_hdirac
Pointer on the vector (which vanishes in Dirac gauge), at the current time step (jtime).
Definition: time_slice.h:574
virtual const Vector & beta_auto() const
Shift function at the current time step jtime.
Definition: isol_hor.C:494
const Map & get_mp() const
Returns the mapping.
Definition: metric.h:202
double omega_hor() const
Orbital velocity.
Definition: phys_param.C:236
virtual const Vector & dnn() const
Covariant derivative of the lapse function at the current time step jtime.
Definition: isol_hor.C:482
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:814
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition: tensor.C:817
virtual const Sym_tensor & gam_dd() const
Induced metric (covariant components ) at the current time step (jtime )
Evolution_std< Sym_tensor > k_uu_evol
Values at successive time steps of the contravariant components of the extrinsic curvature tensor ...
Definition: time_slice.h:216
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
void met_kerr_perturb()
Initialisation of the metric tilde from equation (15) of Dain (2002).
Definition: isol_hor.C:898
Sym_tensor ope_killing_conf(const Metric &gam) const
Computes the conformal Killing operator associated with a given metric.
Definition: vector.C:468
Coord sint
Definition: map.h:739
const Scalar darea_hor() const
Element of area of the horizon.
Definition: phys_param.C:149
virtual void set_hata(const Sym_tensor &hata_in)
Sets the conformal representation of the traceless part of the extrinsic curvature: ...
virtual const Sym_tensor & aa_comp() const
Conformal representation of the traceless part of the extrinsic curvature: Returns the value at the ...
Definition: isol_hor.C:512
const Scalar & get_scalar(int position=0) const
Returns the reference of a Scalar stored in the list.
Definition: param.C:1396
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric_flat.C:156
Map_af & mp
Affine mapping.
Definition: isol_hor.h:260
virtual void sauve(FILE *) const
Save in a file.
Definition: scalar.C:692
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:484
Scalar expansion() const
Expansion of the outgoing null normal ( )
Definition: phys_param.C:266
double boost_x
Boost velocity in x-direction.
Definition: isol_hor.h:272
void set_psi(const Scalar &psi_in)
Sets the conformal factor relating the physical metric to the conformal one: .
Definition: isol_hor.C:518
double axi_break() const
Breaking of the axial symmetry on the horizon.
Definition: isol_hor.C:1090
void init_met_trK()
Sets the 3-metric tilde to the flat metric and gamt_point, trK and trK_point to zero.
Definition: isol_hor.C:788
void set_outer_boundary(int l, double x)
Sets the value of the Scalar at the outer boundary of a given domain.
Definition: scalar_manip.C:321
Evolution_std< Sym_tensor > gam_uu_evol
Values at successive time steps of the contravariant components of the induced metric ...
Definition: time_slice.h:206
double * tet
Array of values of at the nt collocation points.
Definition: grilles.h:217
Evolution_std< Sym_tensor > aa_comp_evol
Values at successive time steps of the components of the conformal representation of the traceless p...
Definition: isol_hor.h:316
Metric met_gamt
3 metric tilde
Definition: isol_hor.h:326
Scalar decouple
Function used to construct from the total .
Definition: isol_hor.h:345
Parameter storage.
Definition: param.h:125
double ang_mom_adm() const
ADM angular Momentum.
Definition: phys_param.C:251
Evolution_std< Sym_tensor > k_dd_evol
Values at successive time steps of the covariant components of the extrinsic curvature tensor ...
Definition: time_slice.h:211
const Metric & gam() const
Induced metric at the current time step (jtime )
Spacelike time slice of a 3+1 spacetime with conformal decomposition.
Definition: time_slice.h:501
Scalar * p_psi4
Pointer on the factor at the current time step (jtime)
Definition: time_slice.h:566
virtual const Scalar & aa_quad() const
Conformal representation .
Definition: isol_hor.C:887
void div_rsint_dzpuis(int ced_mult_r)
Division by but with the output flag dzpuis set to ced_mult_r .
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Definition: fwrite_be.C:73
Metric * p_gamma
Pointer on the induced metric at the current time step (jtime)
Definition: time_slice.h:242
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Evolution_std< Scalar > aa_quad_evol
Values at successive time steps of the components .
Definition: isol_hor.h:323
virtual Sym_tensor aa() const
Conformal representation of the traceless part of the extrinsic curvature: .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
#define T_COSSIN_SI
sin impair-cos pair alternes, sin pour m=0
Definition: type_parite.h:214
Evolution_std< Vector > beta_comp_evol
Values at successive time steps of the shift function .
Definition: isol_hor.h:304
Metric * p_tgamma
Pointer on the conformal metric at the current time step (jtime)
Definition: time_slice.h:563
double radius_hor() const
Radius of the horizon.
Definition: phys_param.C:170
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
Definition: type_parite.h:176
virtual const Scalar & determinant() const
Returns the determinant.
Definition: metric.C:395
double val_point(double r, double theta, double phi) const
Computes the value of the field at an arbitrary point , by means of the spectral expansion.
Definition: scalar.C:896
Evolution_std< double > the_time
Time label of each slice.
Definition: time_slice.h:196
virtual const Scalar & n_auto() const
Lapse function at the current time step jtime.
Definition: isol_hor.C:458
virtual const Sym_tensor & k_uu() const
Extrinsic curvature tensor (contravariant components ) at the current time step (jtime ) ...
Scalar trK_point
Time derivative of the trace of the extrinsic curvature.
Definition: isol_hor.h:335
Evolution_std< Vector > dn_evol
Values at successive time steps of the covariant derivative of the lapse with respect to the flat met...
Definition: isol_hor.h:294
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
Definition: fread_be.C:72
void import(const Scalar &ci)
Assignment to another Scalar defined on a different mapping.
Definition: scalar_import.C:71
void set_gamt(const Metric &gam_tilde)
Sets the conformal metric to gam_tilde.
Definition: isol_hor.C:553
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:283
Bases of the spectral expansions.
Definition: base_val.h:325
void add_scalar(const Scalar &ti, int position=0)
Adds the address of a new Scalar to the list.
Definition: param.C:1351
void set_base_r(int l, int base_r)
Sets the expansion basis for r ( ) functions in a given domain.
Definition: valeur.C:839
Affine radial mapping.
Definition: map.h:2048
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
Definition: time_slice.h:510
void set_nn(const Scalar &nn_in)
Sets the lapse.
Definition: isol_hor.C:543
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:809
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:663
#define T_COSSIN_CI
cos impair-sin pair alternes, cos pour m=0
Definition: type_parite.h:212
Evolution_std< Scalar > trk_evol
Values at successive time steps of the trace K of the extrinsic curvature.
Definition: time_slice.h:227
virtual void sauve(FILE *fich, bool partial_save) const
Total or partial saves in a binary file.
Definition: isol_hor.C:404
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator<<).
Definition: isol_hor.C:381
Evolution_std< Sym_tensor > aa_auto_evol
Values at successive time steps of the components of the conformal representation of the traceless p...
Definition: isol_hor.h:310
Isol_hor(Map_af &mpi, int depth_in=3)
Standard constructor.
Definition: isol_hor.C:179
void div_rsint()
Division by everywhere; dzpuis is not changed.
Evolution_std< Scalar > n_evol
Values at successive time steps of the lapse function N.
Definition: time_slice.h:219
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition: tensor.C:506
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition: param.C:364
virtual double adm_mass() const
Returns the ADM mass (geometrical units) at the current step.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:610
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:302
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
Evolution_std< Sym_tensor > hata_evol
Values at successive time steps of the components .
Definition: time_slice.h:545
Evolution_std< Sym_tensor > aa_nn
Values at successive time steps of the components .
Definition: isol_hor.h:320
virtual const Scalar & psi_auto() const
Conformal factor at the current time step jtime.
Definition: isol_hor.C:470
Evolution_std< Vector > beta_evol
Values at successive time steps of the shift vector .
Definition: time_slice.h:222
double regularise_one()
Corrects the shift in the innermost shell, so that it remains and that equals zero on the horizon...
virtual void sauve(FILE *fich, bool partial_save) const
Total or partial saves in a binary file.
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
Spacelike time-slice of an Isolated Horizon in a 3+1 spacetime with conformal decomposition.
Definition: isol_hor.h:254
int depth
Number of stored time slices.
Definition: time_slice.h:182
double mass_hor() const
Mass computed at the horizon.
Definition: phys_param.C:209
Evolution_std< Vector > dpsi_evol
Values at successive time steps of the covariant derivative of the conformal factor ...
Definition: isol_hor.h:298
virtual const Vector & dpsi() const
Covariant derivative with respect to the flat metric of the conformal factor at the current time ste...
Definition: isol_hor.C:488
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Evolution_std< Sym_tensor > hh_evol
Values at successive time steps of the components .
Definition: time_slice.h:533
const Metric_flat & flat_met_spher() const
Returns the flat metric associated with the spherical coordinates and with components expressed in th...
Definition: map.C:324
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:226
double boost_z
Boost velocity in z-direction.
Definition: isol_hor.h:275
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:607
Coord r
r coordinate centered on the grid
Definition: map.h:736
virtual const Sym_tensor & aa_auto() const
Conformal representation of the traceless part of the extrinsic curvature: Returns the value at the ...
Definition: isol_hor.C:506
virtual const Scalar & n_comp() const
Lapse function at the current time step jtime.
Definition: isol_hor.C:464
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
Definition: valeur.h:490
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166
Coord cost
Definition: map.h:740