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