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
29char 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
175namespace Lorene {
176Isol_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
193Isol_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
222Isol_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
253Isol_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
343void Isol_hor::operator=(const Isol_hor& isolhor_in) {
344
345Time_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
378ostream& 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
401void 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
455const Scalar& Isol_hor::n_auto() const {
456
457 assert( n_auto_evol.is_known(jtime) ) ;
458 return n_auto_evol[jtime] ;
459}
460
461const Scalar& Isol_hor::n_comp() const {
462
463 assert( n_comp_evol.is_known(jtime) ) ;
464 return n_comp_evol[jtime] ;
465}
466
468
469 assert( psi_auto_evol.is_known(jtime) ) ;
470 return psi_auto_evol[jtime] ;
471}
472
474
475 assert( psi_comp_evol.is_known(jtime) ) ;
476 return psi_comp_evol[jtime] ;
477}
478
479const Vector& Isol_hor::dnn() const {
480
481 assert( dn_evol.is_known(jtime) ) ;
482 return dn_evol[jtime] ;
483}
484
485const Vector& Isol_hor::dpsi() const {
486
487 assert( dpsi_evol.is_known(jtime) ) ;
488 return dpsi_evol[jtime] ;
489}
490
492
493 assert( beta_auto_evol.is_known(jtime) ) ;
494 return beta_auto_evol[jtime] ;
495}
496
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
515void 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
540void 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
550void 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
577void 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
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
669void 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
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
712void 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
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 {
863 cout << "regul = " << regul << endl ;
864 Scalar nn_sxpun (division_xpun (Cmp(nn()), 0)) ;
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
884const 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
918void 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) ;
941 ww_pert.set_spectral_va().set_base_t(T_COSSIN_CI) ;
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
1087double 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
1116double 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}
1125void 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}
Bases of the spectral expansions.
Definition base_val.h:322
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
double * phi
Array of values of at the np collocation points.
Definition grilles.h:213
double * tet
Array of values of at the nt collocation points.
Definition grilles.h:211
Spacelike time-slice of an Isolated Horizon in a 3+1 spacetime with conformal decomposition.
Definition isol_hor.h:254
virtual const Vector & dnn() const
Covariant derivative of the lapse function at the current time step jtime
Definition isol_hor.C:479
Scalar expansion() const
Expansion of the outgoing null normal ( )
Definition phys_param.C:263
void met_kerr_perturb()
Initialisation of the metric tilde from equation (15) of Dain (2002).
Definition isol_hor.C:895
virtual const Scalar & aa_quad() const
Conformal representation .
Definition isol_hor.C:884
Evolution_std< Sym_tensor > aa_nn
Values at successive time steps of the components .
Definition isol_hor.h:320
void init_bhole()
Sets the values of the fields to :
Definition isol_hor.C:733
Sym_tensor gamt_point
Time derivative of the 3-metric tilde.
Definition isol_hor.h:329
double regularise_one()
Corrects the shift in the innermost shell, so that it remains and that equals zero on the horizon.
Evolution_std< Scalar > aa_quad_evol
Values at successive time steps of the components .
Definition isol_hor.h:323
double omega
Angular velocity in LORENE's units.
Definition isol_hor.h:269
double regul
Intensity of the correction on the shift vector.
Definition isol_hor.h:278
virtual const Scalar & n_comp() const
Lapse function at the current time step jtime.
Definition isol_hor.C:461
Scalar trK
Trace of the extrinsic curvature.
Definition isol_hor.h:332
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator<<).
Definition isol_hor.C:378
const Scalar darea_hor() const
Element of area of the horizon.
Definition phys_param.C:146
Metric met_gamt
3 metric tilde
Definition isol_hor.h:326
Evolution_std< Vector > dpsi_evol
Values at successive time steps of the covariant derivative of the conformal factor .
Definition isol_hor.h:298
Evolution_std< Scalar > psi_auto_evol
Values at successive time steps of the conformal factor .
Definition isol_hor.h:287
int nz
Number of zones.
Definition isol_hor.h:263
double ang_mom_hor() const
Angular momentum (modulo)
Definition phys_param.C:178
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
void set_psi(const Scalar &psi_in)
Sets the conformal factor relating the physical metric to the conformal one: .
Definition isol_hor.C:515
virtual const Vector & beta_comp() const
Shift function at the current time step jtime.
Definition isol_hor.C:497
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
Scalar trK_point
Time derivative of the trace of the extrinsic curvature.
Definition isol_hor.h:335
Evolution_std< Vector > beta_auto_evol
Values at successive time steps of the shift function .
Definition isol_hor.h:301
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
void set_nn(const Scalar &nn_in)
Sets the lapse.
Definition isol_hor.C:540
double boost_z
Boost velocity in z-direction.
Definition isol_hor.h:275
virtual const Scalar & psi_auto() const
Conformal factor at the current time step jtime.
Definition isol_hor.C:467
double ang_mom_adm() const
ADM angular Momentum
Definition phys_param.C:248
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
Scalar decouple
Function used to construct from the total .
Definition isol_hor.h:345
void update_aa()
Conformal representation of the traceless part of the extrinsic curvature: .
Definition isol_hor.C:842
virtual const Scalar & n_auto() const
Lapse function at the current time step jtime.
Definition isol_hor.C:455
virtual ~Isol_hor()
Destructor.
Definition isol_hor.C:336
double radius
Radius of the horizon in LORENE's units.
Definition isol_hor.h:266
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
double mass_hor() const
Mass computed at the horizon
Definition phys_param.C:206
void init_bhole_seul()
Initiates for a single black hole.
Definition isol_hor.C:797
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
void set_gamt(const Metric &gam_tilde)
Sets the conformal metric to gam_tilde.
Definition isol_hor.C:550
double radius_hor() const
Radius of the horizon.
Definition phys_param.C:167
Map_af & mp
Affine mapping.
Definition isol_hor.h:260
virtual const Scalar & psi_comp() const
Conformal factor at the current time step jtime.
Definition isol_hor.C:473
double axi_break() const
Breaking of the axial symmetry on the horizon.
Definition isol_hor.C:1087
Evolution_std< Scalar > n_auto_evol
Values at successive time steps of the lapse function .
Definition isol_hor.h:281
Evolution_std< Scalar > n_comp_evol
Values at successive time steps of the lapse function .
Definition isol_hor.h:284
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
Isol_hor(Map_af &mpi, int depth_in=3)
Standard constructor.
Definition isol_hor.C:176
double boost_x
Boost velocity in x-direction.
Definition isol_hor.h:272
void operator=(const Isol_hor &)
Assignment to another Isol_hor.
Definition isol_hor.C:343
Evolution_std< Vector > beta_comp_evol
Values at successive time steps of the shift function .
Definition isol_hor.h:304
double area_hor() const
Area of the horizon.
Definition phys_param.C:157
virtual const Vector & beta_auto() const
Shift function at the current time step jtime.
Definition isol_hor.C:491
virtual void sauve(FILE *fich, bool partial_save) const
Total or partial saves in a binary file.
Definition isol_hor.C:401
double omega_hor() const
Orbital velocity
Definition phys_param.C:233
Evolution_std< Scalar > psi_comp_evol
Values at successive time steps of the lapse function .
Definition isol_hor.h:290
Affine radial mapping.
Definition map.h:2027
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Definition map.h:791
Coord sint
Definition map.h:721
Coord r
r coordinate centered on the grid
Definition map.h:718
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
Coord cost
Definition map.h:722
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
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
Flat metric for tensor calculation.
Definition metric.h:261
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the metric is defined.
Definition metric.h:309
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Metric for tensor calculation.
Definition metric.h:90
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition metric.C:290
const Map & get_mp() const
Returns the mapping.
Definition metric.h:202
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition metric.C:280
virtual const Scalar & determinant() const
Returns the determinant.
Definition metric.C:392
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for domain no. l.
Definition grilles.h:500
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
Parameter storage.
Definition param.h:125
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition param.C:361
const Scalar & get_scalar(int position=0) const
Returns the reference of a Scalar stored in the list.
Definition param.C:1393
void add_scalar(const Scalar &ti, int position=0)
Adds the address of a new Scalar to the list.
Definition param.C:1348
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
virtual void sauve(FILE *) const
Save in a file.
Definition scalar.C:686
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
void div_rsint_dzpuis(int ced_mult_r)
Division by but with the output flag dzpuis set to ced_mult_r .
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
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Definition scalar.h:615
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition scalar.C:391
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition scalar.h:601
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
void div_rsint()
Division by everywhere; dzpuis is not changed.
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
void set_outer_boundary(int l, double x)
Sets the value of the Scalar at the outer boundary of a given domain.
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...
void import(const Scalar &ci)
Assignment to another Scalar defined on a different mapping.
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
Sym_tensor derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
Definition sym_tensor.C:360
Spacelike time slice of a 3+1 spacetime with conformal decomposition.
Definition time_slice.h:498
Evolution_std< Sym_tensor > hh_evol
Values at successive time steps of the components .
Definition time_slice.h:530
virtual const Sym_tensor & gam_dd() const
Induced metric (covariant components ) at the current time step (jtime )
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )
virtual const Sym_tensor & k_uu() const
Extrinsic curvature tensor (contravariant components ) at the current time step (jtime )
Evolution_std< Sym_tensor > hata_evol
Values at successive time steps of the components .
Definition time_slice.h:542
void operator=(const Time_slice_conf &)
Assignment to another Time_slice_conf.
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
Metric * p_tgamma
Pointer on the conformal metric at the current time step (jtime)
Definition time_slice.h:560
virtual Sym_tensor aa() const
Conformal representation of the traceless part of the extrinsic curvature: .
virtual void set_hata(const Sym_tensor &hata_in)
Sets the conformal representation of the traceless part of the extrinsic curvature: .
Vector * p_hdirac
Pointer on the vector (which vanishes in Dirac gauge), at the current time step (jtime).
Definition time_slice.h:571
const Scalar & psi4() const
Factor at the current time step (jtime ).
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
virtual const Sym_tensor & gam_uu() const
Induced metric (contravariant components ) at the current time step (jtime )
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
Definition time_slice.h:507
Scalar * p_ln_psi
Pointer on the logarithm of at the current time step (jtime)
Definition time_slice.h:566
Scalar * p_psi4
Pointer on the factor at the current time step (jtime)
Definition time_slice.h:563
int jtime
Time step index of the latest slice.
Definition time_slice.h:190
Evolution_std< Scalar > trk_evol
Values at successive time steps of the trace K of the extrinsic curvature.
Definition time_slice.h:224
Metric * p_gamma
Pointer on the induced metric at the current time step (jtime)
Definition time_slice.h:239
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
virtual void sauve(FILE *fich, bool partial_save) const
Total or partial saves in a binary file.
Definition time_slice.C:507
Evolution_full< Tbl > adm_mass_evol
ADM mass at each time step, since the creation of the slice.
Definition time_slice.h:233
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
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
virtual const Vector & beta() const
shift vector at the current time step (jtime )
const Metric & gam() const
Induced metric at the current time step (jtime )
int depth
Number of stored time slices.
Definition time_slice.h:179
Evolution_std< Scalar > n_evol
Values at successive time steps of the lapse function N.
Definition time_slice.h:216
virtual double adm_mass() const
Returns the ADM mass (geometrical units) at the current step.
Evolution_std< double > the_time
Time label of each slice.
Definition time_slice.h:193
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator<<).
Definition time_slice.C:411
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
Evolution_std< Vector > beta_evol
Values at successive time steps of the shift vector .
Definition time_slice.h:219
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
Definition valeur.h:480
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition valeur.C:810
Tensor field of valence 1.
Definition vector.h:188
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:316
Sym_tensor ope_killing_conf(const Metric &gam) const
Computes the conformal Killing operator associated with a given metric.
Definition vector.C:462
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition cmp_math.C:481
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#define T_COSSIN_SI
sin impair-cos pair alternes, sin pour m=0
#define T_COSSIN_CI
cos impair-sin pair alternes, cos pour m=0
#define P_COSSIN
dev. standart
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
#define R_CHEB
base de Chebychev ordinaire (fin)
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
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.
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
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
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
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
virtual void sauve(FILE *) const
Save in a binary file.
Definition tensor_sym.C:372
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
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tensor.h:866
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition tensor.C:481
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition tensor.C:497
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:654
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition app_hor.h:64