LORENE
binhor_equations.C
1/*
2 * Copyright- (c) 2004-2005 Francois Limousin
3 * Jose-Luis Jaramillo
4 *
5 * This file is part of LORENE.
6 *
7 * LORENE is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * LORENE is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with LORENE; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 *
21 */
22
23
24char binhor_equations_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bin_hor/binhor_equations.C,v 1.21 2014/10/13 08:52:42 j_novak Exp $" ;
25
26/*
27 * $Id: binhor_equations.C,v 1.21 2014/10/13 08:52:42 j_novak Exp $
28 * $Log: binhor_equations.C,v $
29 * Revision 1.21 2014/10/13 08:52:42 j_novak
30 * Lorene classes and functions now belong to the namespace Lorene.
31 *
32 * Revision 1.20 2014/10/06 15:13:01 j_novak
33 * Modified #include directives to use c++ syntax.
34 *
35 * Revision 1.19 2008/02/06 18:20:02 f_limousin
36 * Fixed an error in the triad
37 *
38 * Revision 1.18 2007/08/22 16:12:33 f_limousin
39 * Changed the name of the function computing \tilde{\gamma}_{ij}
40 *
41 * Revision 1.17 2007/04/13 15:28:55 f_limousin
42 * Lots of improvements, generalisation to an arbitrary state of
43 * rotation, implementation of the spatial metric given by Samaya.
44 *
45 * Revision 1.16 2006/08/01 14:37:19 f_limousin
46 * New version
47 *
48 * Revision 1.15 2006/06/29 08:51:00 f_limousin
49 * *** empty log message ***
50 *
51 * Revision 1.14 2006/05/24 16:56:37 f_limousin
52 * Many small modifs.
53 *
54 * Revision 1.13 2005/11/15 14:04:00 f_limousin
55 * Minor change to control the resolution of the equation for psi.
56 *
57 * Revision 1.12 2005/10/23 16:39:43 f_limousin
58 * Simplification of the equation in the case of a conformally
59 * flat metric and maximal slicing
60 *
61 * Revision 1.11 2005/09/13 18:33:15 f_limousin
62 * New function vv_bound_cart_bin(double) for computing binaries with
63 * berlin condition for the shift vector.
64 * Suppress all the symy and asymy in the importations.
65 *
66 * Revision 1.10 2005/07/11 08:21:57 f_limousin
67 * Implementation of a new boundary condition for the lapse in the binary
68 * case : boundary_nn_Dir_lapl().
69 *
70 * Revision 1.9 2005/06/09 16:12:04 f_limousin
71 * Implementation of amg_mom_adm().
72 *
73 * Revision 1.8 2005/04/29 14:02:44 f_limousin
74 * Important changes : manage the dependances between quantities (for
75 * instance psi and psi4). New function write_global(ost).
76 *
77 * Revision 1.7 2005/03/10 17:21:52 f_limousin
78 * Add the Berlin boundary condition for the shift.
79 * Some changes to avoid warnings.
80 *
81 * Revision 1.6 2005/03/03 10:29:02 f_limousin
82 * Delete omega in the parameters of the function boundary_beta_z().
83 *
84 * Revision 1.5 2005/02/24 17:25:23 f_limousin
85 * The boundary conditions for psi, N and beta are now parameters in
86 * par_init.d and par_coal.d.
87 *
88 * Revision 1.4 2005/02/11 18:21:38 f_limousin
89 * Dirichlet_binaire and neumann_binaire are taking Scalars as arguments
90 * instead of Cmps.
91 *
92 * Revision 1.3 2005/02/07 10:46:28 f_limousin
93 * Many changes !! The sources are written differently to minimize the
94 * numerical error, the boundary conditions are changed...
95 *
96 * Revision 1.2 2004/12/31 15:41:26 f_limousin
97 * Correction of an error
98 *
99 * Revision 1.1 2004/12/29 16:11:34 f_limousin
100 * First version
101 *
102 *
103 * $Header: /cvsroot/Lorene/C++/Source/Bin_hor/binhor_equations.C,v 1.21 2014/10/13 08:52:42 j_novak Exp $
104 *
105 */
106
107//standard
108#include <cstdlib>
109#include <cmath>
110
111// Lorene
112#include "nbr_spx.h"
113#include "tensor.h"
114#include "tenseur.h"
115#include "isol_hor.h"
116#include "proto.h"
117#include "utilitaires.h"
118//#include "graphique.h"
119
120// Resolution for the lapse
121// ------------------------
122
123namespace Lorene {
124void Bin_hor::solve_lapse (double precision, double relax, int bound_nn,
125 double lim_nn) {
126
127 assert ((relax >0) && (relax<=1)) ;
128
129 cout << "-----------------------------------------------" << endl ;
130 cout << "Resolution LAPSE" << endl ;
131
132 Scalar lapse_un_old (hole1.n_auto) ;
133 Scalar lapse_deux_old (hole2.n_auto) ;
134
135 Sym_tensor taa_un = hole1.aa.up_down(hole1.tgam) ;
136 Scalar aa_quad_un = contract(taa_un, 0, 1, hole1.aa_auto, 0, 1) ;
137
138 Sym_tensor taa_deux = hole2.aa.up_down(hole2.tgam) ;
139 Scalar aa_quad_deux = contract(taa_deux, 0, 1, hole2.aa_auto, 0, 1) ;
140
141 Tensor hdirac1 (contract((hole1.hh).derive_cov(hole1.ff),0,2)) ;
142 Tensor hdirac2 (contract((hole2.hh).derive_cov(hole2.ff),0,2)) ;
143
144 // Source 1
145 // --------
146
147 Scalar source_un (hole1.mp) ;
148
149 // Conformally flat
150 /*
151 source_un = hole1.get_psi4()*hole1.nn*aa_quad_un
152 -2.*contract(hole1.dn, 0, hole1.psi_auto
153 .derive_con(hole1.ff), 0)/hole1.psi ;
154 */
155
156 source_un = hole1.get_psi4()*( hole1.nn*( aa_quad_un + 0.3333333333333333 *
159
161 .derive_cov(hole1.ff), 0), 0, hole1.dpsi, 0)/hole1.psi
162
164 .derive_con(hole1.ff), 0)/hole1.psi ;
165
166 - contract(hdirac1, 0, hole1.n_auto.derive_cov(hole1.ff), 0) ;
167
168
169 Scalar tmp_un (hole1.mp) ;
170
172 derive_cov(hole1.ff), 0)
174 .derive_cov(hole1.ff), 0, 1 ) ;
175
176 tmp_un.inc_dzpuis() ; // dzpuis: 3 -> 4
177
178 source_un += tmp_un ;
179
180
181 // Source 2
182 // ---------
183
184 Scalar source_deux (hole2.mp) ;
185
186 // Conformally flat
187 /*
188 source_deux = hole2.get_psi4()*hole2.nn*aa_quad_deux
189 -2.*contract(hole2.dn, 0, hole2.psi_auto
190 .derive_con(hole2.ff), 0)/hole2.psi ;
191 */
192
193 source_deux = hole2.get_psi4()*( hole2.nn*( aa_quad_deux + 0.3333333333333333
196
198 .derive_cov(hole2.ff), 0), 0, hole2.dpsi, 0)/hole2.psi
199
201 .derive_con(hole2.ff), 0)/hole2.psi ;
202
203 - contract(hdirac2, 0, hole2.n_auto.derive_cov(hole2.ff), 0) ;
204
205 Scalar tmp_deux (hole2.mp) ;
206
207 tmp_deux = hole2.get_psi4()* contract(hole2.beta_auto, 0, hole2.trK.
208 derive_cov(hole2.ff), 0)
210 .derive_cov(hole2.ff), 0, 1 ) ;
211
212 tmp_deux.inc_dzpuis() ; // dzpuis: 3 -> 4
213
214 source_deux += tmp_deux ;
215
216 cout << "source lapse" << endl << norme(source_un) << endl ;
217
218 // Boundary conditions and resolution
219 // -----------------------------------
220
221 Valeur lim_un (hole1.mp.get_mg()-> get_angu()) ;
222 Valeur lim_deux (hole1.mp.get_mg()-> get_angu()) ;
223
224 Scalar n_un_temp (hole1.n_auto) ;
225 Scalar n_deux_temp (hole2.n_auto) ;
226
227 switch (bound_nn) {
228
229 case 0 : {
230
231 lim_un = hole1.boundary_nn_Dir(lim_nn) ;
232 lim_deux = hole2.boundary_nn_Dir(lim_nn) ;
233
234 n_un_temp = n_un_temp - 1./2. ;
235 n_deux_temp = n_deux_temp - 1./2. ;
236
237 dirichlet_binaire (source_un, source_deux, lim_un, lim_deux,
238 n_un_temp, n_deux_temp, 0, precision) ;
239 break ;
240 }
241
242 case 1 : {
243
244 lim_un = hole1.boundary_nn_Neu(lim_nn) ;
245 lim_deux = hole2.boundary_nn_Neu(lim_nn) ;
246
247 neumann_binaire (source_un, source_deux, lim_un, lim_deux,
248 n_un_temp, n_deux_temp, 0, precision) ;
249 break ;
250 }
251
252 default : {
253 cout << "Unexpected type of boundary conditions for the lapse!"
254 << endl
255 << " bound_nn = " << bound_nn << endl ;
256 abort() ;
257 break ;
258 }
259
260 } // End of switch
261
262 n_un_temp = n_un_temp + 1./2. ;
263 n_deux_temp = n_deux_temp + 1./2. ;
264
265 n_un_temp.raccord(3) ;
266 n_deux_temp.raccord(3) ;
267
268 // Check: has the Poisson equation been correctly solved ?
269 // -----------------------------------------------------
270
271 int nz = hole1.mp.get_mg()->get_nzone() ;
272 cout << "lapse auto" << endl << norme (n_un_temp) << endl ;
273 Tbl tdiff_nn = diffrel(n_un_temp.laplacian(), source_un) ;
274
275 cout <<
276 "Relative error in the resolution of the equation for the lapse : "
277 << endl ;
278 for (int l=0; l<nz; l++) {
279 cout << tdiff_nn(l) << " " ;
280 }
281 cout << endl ;
282
283 // Relaxation :
284 // -------------
285
286 n_un_temp = relax*n_un_temp + (1-relax)*lapse_un_old ;
287 n_deux_temp = relax*n_deux_temp + (1-relax)*lapse_deux_old ;
288
289 hole1.n_auto = n_un_temp ;
290 hole2.n_auto = n_deux_temp ;
291
294
295}
296
297
298// Resolution for Psi
299// -------------------
300
301void Bin_hor::solve_psi (double precision, double relax, int bound_psi) {
302
303 assert ((relax>0) && (relax<=1)) ;
304
305 cout << "-----------------------------------------------" << endl ;
306 cout << "Resolution PSI" << endl ;
307
308 Scalar psi_un_old (hole1.psi_auto) ;
309 Scalar psi_deux_old (hole2.psi_auto) ;
310
311 Sym_tensor taa_un = hole1.aa.up_down(hole1.tgam) ;
312 Scalar aa_quad_un = contract(taa_un, 0, 1, hole1.aa_auto, 0, 1) ;
313
314 Sym_tensor taa_deux = hole2.aa.up_down(hole2.tgam) ;
315 Scalar aa_quad_deux = contract(taa_deux, 0, 1, hole2.aa_auto, 0, 1) ;
316
317 Tensor hdirac1 (contract((hole1.hh).derive_cov(hole1.ff),0,2)) ;
318 Tensor hdirac2 (contract((hole2.hh).derive_cov(hole2.ff),0,2)) ;
319
320 // Source 1
321 // ---------
322
323 Scalar source_un (hole1.mp) ;
324 /*
325 // Conformally flat source
326 source_un.annule_hard() ;
327 source_un.set_dzpuis(4) ;
328 source_un += - hole1.psi*hole1.get_psi4()* 0.125* aa_quad_un ;
329 source_un.std_spectral_base() ;
330 */
331
332 Scalar tmp_un (hole1.mp) ;
333
334 tmp_un = 0.125* hole1.psi_auto * (hole1.tgam).ricci_scal()
336 .derive_cov(hole1.ff), 0, 1 ) ;
337 tmp_un.inc_dzpuis() ; // dzpuis : 3 -> 4
338
339 tmp_un -= contract(hdirac1, 0, hole1.psi_auto
340 .derive_cov(hole1.ff), 0) ;
341
342 source_un = tmp_un - hole1.psi*hole1.get_psi4()* ( 0.125* aa_quad_un
343 - 8.33333333333333e-2* hole1.trK*hole1.trK*hole1.decouple ) ;
344 source_un.std_spectral_base() ;
345
346
347
348 // Source 2
349 // ---------
350
351 Scalar source_deux (hole2.mp) ;
352 /*
353 // Conformally flat source
354 source_deux.annule_hard() ;
355 source_deux.set_dzpuis(4) ;
356 source_deux += - hole2.psi*hole2.get_psi4()* 0.125* aa_quad_deux ;
357 source_deux.std_spectral_base() ;
358 */
359
360
361 Scalar tmp_deux (hole2.mp) ;
362
363 tmp_deux = 0.125* hole2.psi_auto * (hole2.tgam).ricci_scal()
365 .derive_cov(hole2.ff), 0, 1 ) ;
366 tmp_deux.inc_dzpuis() ; // dzpuis : 3 -> 4
367
368 tmp_deux -= contract(hdirac2, 0, hole2.psi_auto
369 .derive_cov(hole2.ff), 0) ;
370
371 source_deux = tmp_deux - hole2.psi*hole2.get_psi4()* ( 0.125* aa_quad_deux
372 - 8.33333333333333e-2* hole2.trK*hole2.trK*hole2.decouple ) ;
373 source_deux.std_spectral_base() ;
374
375
376 cout << "source psi" << endl << norme(source_un) << endl ;
377
378 // Boundary conditions and resolution :
379 // ------------------------------------
380
381 Valeur lim_un (hole1.mp.get_mg()-> get_angu()) ;
382 Valeur lim_deux (hole1.mp.get_mg()-> get_angu()) ;
383
384 Scalar psi_un_temp (hole1.psi_auto) ;
385 Scalar psi_deux_temp (hole2.psi_auto) ;
386
387 switch (bound_psi) {
388
389 case 0 : {
390
391 lim_un = hole1.boundary_psi_app_hor() ;
392 lim_deux = hole2.boundary_psi_app_hor() ;
393
394 neumann_binaire (source_un, source_deux, lim_un, lim_deux,
395 psi_un_temp, psi_deux_temp, 0, precision) ;
396 break ;
397 }
398
399 default : {
400 cout << "Unexpected type of boundary conditions for psi!"
401 << endl
402 << " bound_psi = " << bound_psi << endl ;
403 abort() ;
404 break ;
405 }
406
407 } // End of switch
408
409 psi_un_temp = psi_un_temp + 1./2. ;
410 psi_deux_temp = psi_deux_temp + 1./2. ;
411
412 psi_un_temp.raccord(3) ;
413 psi_deux_temp.raccord(3) ;
414
415 // Check: has the Poisson equation been correctly solved ?
416 // -----------------------------------------------------
417
418 int nz = hole1.mp.get_mg()->get_nzone() ;
419 cout << "psi auto" << endl << norme (psi_un_temp) << endl ;
420 Tbl tdiff_psi = diffrel(psi_un_temp.laplacian(), source_un) ;
421
422 cout <<
423 "Relative error in the resolution of the equation for psi : "
424 << endl ;
425 for (int l=0; l<nz; l++) {
426 cout << tdiff_psi(l) << " " ;
427 }
428 cout << endl ;
429
430 // Relaxation :
431 // -------------
432
433 psi_un_temp = relax*psi_un_temp + (1-relax)*psi_un_old ;
434 psi_deux_temp = relax*psi_deux_temp + (1-relax)*psi_deux_old ;
435
436 hole1.psi_auto = psi_un_temp ;
437 hole2.psi_auto = psi_deux_temp ;
438
441
444
445 //set_hh_Samaya() ;
446
447}
448
449
450// Resolution for shift with omega fixed.
451// --------------------------------------
452
453void Bin_hor::solve_shift (double precision, double relax, int bound_beta,
454 double omega_eff) {
455
456 cout << "------------------------------------------------" << endl ;
457 cout << "Resolution shift : Omega = " << omega << endl ;
458
459 Sym_tensor taa_un = hole1.aa.up_down(hole1.tgam) ;
460 Scalar aa_quad_un = contract(taa_un, 0, 1, hole1.aa_auto, 0, 1) ;
461
462 Sym_tensor taa_deux = hole2.aa.up_down(hole2.tgam) ;
463 Scalar aa_quad_deux = contract(taa_deux, 0, 1, hole2.aa_auto, 0, 1) ;
464
465 Tensor hdirac1 (contract((hole1.hh).derive_cov(hole1.ff),0,2)) ;
466 Tensor hdirac2 (contract((hole2.hh).derive_cov(hole2.ff),0,2)) ;
467
468 // Source 1
469 // ---------
470
471 Vector source_un (hole1.mp, CON, hole1.mp.get_bvect_spher()) ;
472 /*
473 // Conformally flat source
474 source_un = 2.* contract(hole1.aa_auto, 1, hole1.dn, 0)
475 - 12.*hole1.nn*contract(hole1.aa_auto, 1, hole1.dpsi, 0)
476 /hole1.psi;
477 */
478
479
480 Vector tmp_vect_un (hole1.mp, CON, hole1.mp.get_bvect_spher()) ;
481
482 source_un = 2.* contract(hole1.aa_auto, 1, hole1.dn, 0)
483 - 12.*hole1.nn*contract(hole1.aa_auto, 1, hole1.dpsi, 0)
484 /hole1.psi;
485
486
487 tmp_vect_un = 2./3.* hole1.trK.derive_con(hole1.tgam)
488 * hole1.decouple ;
489 tmp_vect_un.inc_dzpuis() ;
490
491 source_un += 2.* hole1.nn * ( tmp_vect_un
492 - contract(hole1.tgam.connect().get_delta(), 1, 2,
493 hole1.aa_auto, 0, 1) ) ;
494
495 Vector vtmp_un = contract(hole1.hh, 0, 1,
497 .derive_cov(hole1.ff), 1, 2)
498 + 1./3.*contract(hole1.hh, 1, hole1.beta_auto
500 - hdirac1.derive_lie(hole1.beta_auto)
502 vtmp_un.inc_dzpuis() ; // dzpuis: 3 -> 4
503
504 source_un -= vtmp_un ;
505
506 source_un += 2./3.* hole1.beta_auto.divergence(hole1.ff)
507 * hdirac1 ;
508
509 source_un.std_spectral_base() ;
510
511
512 // Source 2
513 // ---------
514
515 Vector source_deux (hole2.mp, CON, hole2.mp.get_bvect_spher()) ;
516 /*
517 // Conformally flat source
518 source_deux = 2.* contract(hole2.aa_auto, 1, hole2.dn, 0)
519 - 12.*hole2.nn*contract(hole2.aa_auto, 1, hole2.dpsi, 0)
520 /hole2.psi;
521 */
522
523
524 Vector tmp_vect_deux (hole2.mp, CON, hole2.mp.get_bvect_spher()) ;
525
526 source_deux = 2.* contract(hole2.aa_auto, 1, hole2.dn, 0)
527 - 12.*hole2.nn*contract(hole2.aa_auto, 1, hole2.dpsi, 0)
528 /hole2.psi;
529
530 tmp_vect_deux = 2./3.* hole2.trK.derive_con(hole2.tgam)
531 * hole2.decouple ;
532 tmp_vect_deux.inc_dzpuis() ;
533
534 source_deux += 2.* hole2.nn * ( tmp_vect_deux
535 - contract(hole2.tgam.connect().get_delta(), 1, 2,
536 hole2.aa*hole2.decouple, 0, 1) ) ;
537
538 Vector vtmp_deux = contract(hole2.hh, 0, 1,
540 .derive_cov(hole2.ff), 1, 2)
541 + 1./3.*contract(hole2.hh, 1, hole2.beta_auto
543 - hdirac2.derive_lie(hole2.beta_auto)
545 vtmp_deux.inc_dzpuis() ; // dzpuis: 3 -> 4
546
547 source_deux -= vtmp_deux ;
548
549 source_deux += 2./3.* hole2.beta_auto.divergence(hole2.ff)
550 * hdirac2 ;
551
552 source_deux.std_spectral_base() ;
553
554
555 Vector source_1 (source_un) ;
556 Vector source_2 (source_deux) ;
557 source_1.change_triad(hole1.mp.get_bvect_cart()) ;
558 source_2.change_triad(hole2.mp.get_bvect_cart()) ;
559
560 cout << "source shift_x" << endl << norme(source_1(1)) << endl ;
561 cout << "source shift_y" << endl << norme(source_1(2)) << endl ;
562 cout << "source shift_z" << endl << norme(source_1(3)) << endl ;
563
564 // Filter for high frequencies.
565 for (int i=1 ; i<=3 ; i++) {
566 source_un.set(i).filtre(4) ;
567 source_deux.set(i).filtre(4) ;
568 }
569
570 // Boundary conditions
571 // --------------------
572
573 Valeur lim_x_un (hole1.mp.get_mg()-> get_angu()) ;
574 Valeur lim_y_un (hole1.mp.get_mg()-> get_angu()) ;
575 Valeur lim_z_un (hole1.mp.get_mg()-> get_angu()) ;
576
577 Valeur lim_x_deux (hole2.mp.get_mg()-> get_angu()) ;
578 Valeur lim_y_deux (hole2.mp.get_mg()-> get_angu()) ;
579 Valeur lim_z_deux (hole2.mp.get_mg()-> get_angu()) ;
580
581 switch (bound_beta) {
582
583 case 0 : {
584
585 lim_x_un = hole1.boundary_beta_x(omega, omega_eff) ;
586 lim_y_un = hole1.boundary_beta_y(omega, omega_eff) ;
587 lim_z_un = hole1.boundary_beta_z() ;
588
589 lim_x_deux = hole2.boundary_beta_x(omega, omega_eff) ;
590 lim_y_deux = hole2.boundary_beta_y(omega, omega_eff) ;
591 lim_z_deux = hole2.boundary_beta_z() ;
592 break ;
593 }
594
595 default : {
596 cout << "Unexpected type of boundary conditions for beta!"
597 << endl
598 << " bound_beta = " << bound_beta << endl ;
599 abort() ;
600 break ;
601 }
602
603 } // End of switch
604
605
606 // We solve :
607 // -----------
608
609 Vector beta_un_old (hole1.beta_auto) ;
610 Vector beta_deux_old (hole2.beta_auto) ;
611 Vector beta1 (hole1.beta_auto) ;
612 Vector beta2 (hole2.beta_auto) ;
613
614 poisson_vect_binaire (1./3., source_un, source_deux,
615 lim_x_un, lim_y_un, lim_z_un,
616 lim_x_deux, lim_y_deux, lim_z_deux,
617 beta1, beta2, 0, precision) ;
618
619
622
623 for (int i=1 ; i<=3 ; i++) {
624 beta1.set(i).raccord(3) ;
625 beta2.set(i).raccord(3) ;
626 }
627
628 cout << "shift_auto x" << endl << norme(beta1(1)) << endl ;
629 cout << "shift_auto y" << endl << norme(beta1(2)) << endl ;
630 cout << "shift_auto z" << endl << norme(beta1(3)) << endl ;
631
634
635 // Check: has the Poisson equation been correctly solved ?
636 // -----------------------------------------------------
637
638 int nz = hole1.mp.get_mg()->get_nzone() ;
639 Vector lap_beta = (beta1.derive_con(hole1.ff)).divergence(hole1.ff)
640 + 1./3.* beta1.divergence(hole1.ff).derive_con(hole1.ff) ;
641 source_un.dec_dzpuis() ;
642
643 Tbl tdiff_beta_r = diffrel(lap_beta(1), source_un(1)) ;
644 Tbl tdiff_beta_t = diffrel(lap_beta(2), source_un(2)) ;
645 Tbl tdiff_beta_p = diffrel(lap_beta(3), source_un(3)) ;
646
647 cout <<
648 "Relative error in the resolution of the equation for beta : "
649 << endl ;
650 cout << "r component : " ;
651 for (int l=0; l<nz; l++) {
652 cout << tdiff_beta_r(l) << " " ;
653 }
654 cout << endl ;
655 cout << "t component : " ;
656 for (int l=0; l<nz; l++) {
657 cout << tdiff_beta_t(l) << " " ;
658 }
659 cout << endl ;
660 cout << "p component : " ;
661 for (int l=0; l<nz; l++) {
662 cout << tdiff_beta_p(l) << " " ;
663 }
664 cout << endl ;
665
666
667 // Relaxation
668 // -----------
669
670 Vector beta1_new (hole1.mp, CON, hole1.mp.get_bvect_spher()) ;
671 Vector beta2_new (hole2.mp, CON, hole2.mp.get_bvect_spher()) ;
672
673
674 // Construction of Omega d/dphi
675 // ----------------------------
676
677 Vector omdsdp1 (hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
678 Scalar yya1 (hole1.mp) ;
679 yya1 = hole1.mp.ya ;
680 Scalar xxa1 (hole1.mp) ;
681 xxa1 = hole1.mp.xa ;
682
683 if (fabs(hole1.mp.get_rot_phi()) < 1e-10){
684 omdsdp1.set(1) = - omega * yya1 ;
685 omdsdp1.set(2) = omega * xxa1 ;
686 omdsdp1.set(3).annule_hard() ;
687 }
688 else{
689 omdsdp1.set(1) = omega * yya1 ;
690 omdsdp1.set(2) = - omega * xxa1 ;
691 omdsdp1.set(3).annule_hard() ;
692 }
693
694 omdsdp1.set(1).set_spectral_va()
696 omdsdp1.set(2).set_spectral_va()
698 omdsdp1.set(3).set_spectral_va()
700
701 omdsdp1.annule_domain(nz-1) ;
703
704
705 Vector omdsdp2 (hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
706 Scalar yya2 (hole2.mp) ;
707 yya2 = hole2.mp.ya ;
708 Scalar xxa2 (hole2.mp) ;
709 xxa2 = hole2.mp.xa ;
710
711 if (fabs(hole2.mp.get_rot_phi()) < 1e-10){
712 omdsdp2.set(1) = - omega * yya2 ;
713 omdsdp2.set(2) = omega * xxa2 ;
714 omdsdp2.set(3).annule_hard() ;
715 }
716 else{
717 omdsdp2.set(1) = omega * yya2 ;
718 omdsdp2.set(2) = - omega * xxa2 ;
719 omdsdp2.set(3).annule_hard() ;
720 }
721
722 omdsdp2.set(1).set_spectral_va()
724 omdsdp2.set(2).set_spectral_va()
726 omdsdp2.set(3).set_spectral_va()
728
729 omdsdp2.annule_domain(nz-1) ;
731
732 // New shift
733 beta1_new = relax*(beta1+hole1.decouple*omdsdp1) + (1-relax)*beta_un_old ;
734 beta2_new = relax*(beta2+hole2.decouple*omdsdp2) + (1-relax)*beta_deux_old ;
735
736 hole1.beta_auto = beta1_new ;
737 hole2.beta_auto = beta2_new ;
738
741
742 // Regularisation of the shifts if necessary
743 // -----------------------------------------
744
745 int nnt = hole1.mp.get_mg()->get_nt(1) ;
746 int nnp = hole1.mp.get_mg()->get_np(1) ;
747
748 int check ;
749 check = 0 ;
750 for (int k=0; k<nnp; k++)
751 for (int j=0; j<nnt; j++){
752 if (fabs((hole1.n_auto+hole1.n_comp).val_grid_point(1, k, j , 0)) < 1e-8){
753 check = 1 ;
754 break ;
755 }
756 }
757
758 if (check == 1){
759 double diff_un = hole1.regularisation (hole1.beta_auto,
761 double diff_deux = hole2.regularisation (hole2.beta_auto,
763 hole1.regul = diff_un ;
764 hole2.regul = diff_deux ;
765 }
766
767 else {
768 hole1.regul = 0. ;
769 hole2.regul = 0. ;
770 }
771
772}
773}
void solve_psi(double precis, double relax, int bound_psi)
Solves the equation for the conformal factor : The fields are the total values excpet those with subs...
double omega
Angular velocity.
Definition isol_hor.h:1348
Single_hor hole1
Black hole one.
Definition isol_hor.h:1342
Single_hor hole2
Black hole two.
Definition isol_hor.h:1343
void solve_shift(double precis, double relax, int bound_beta, double omega_eff)
Solves the equation for the shift, using the Oohara-Nakarmure scheme : The fields are the total value...
void solve_lapse(double precis, double relax, int bound_nn, double lim_nn)
Solves the equation for the lapse : The fields are the total values except those with subscript ,...
const Tensor_sym & get_delta() const
Returns the tensor which defines the connection with respect to the flat one: is the difference bet...
Definition connection.h:271
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 ya
Absolute y coordinate.
Definition map.h:731
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition map.h:775
Coord xa
Absolute x coordinate.
Definition map.h:730
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
virtual const Connection & connect() const
Returns the connection.
Definition metric.C:301
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
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 filtre(int n)
Sets the n lasts coefficients in r to 0 in the external domain.
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
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 inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Vector dpsi
Covariant derivative of the conformal factor .
Definition isol_hor.h:941
Scalar psi_auto
Conformal factor .
Definition isol_hor.h:924
Scalar trK_point
Time derivative of the trace of the extrinsic curvature.
Definition isol_hor.h:992
const Valeur boundary_psi_app_hor() const
Neumann boundary condition for.
Scalar decouple
Function used to construct from the total .
Definition isol_hor.h:1002
Vector beta_auto
Shift function .
Definition isol_hor.h:944
void beta_comp_import(const Single_hor &comp)
Imports the part of due to the companion hole comp.
Definition single_hor.C:465
Scalar n_comp
Lapse function .
Definition isol_hor.h:918
Sym_tensor aa_auto
Components of the conformal representation of the traceless part of the extrinsic curvature:
Definition isol_hor.h:959
const Scalar & get_psi4() const
Conformal factor .
Definition single_hor.C:288
Metric_flat ff
3 metric flat
Definition isol_hor.h:980
Scalar psi
Conformal factor .
Definition isol_hor.h:930
Map_af & mp
Affine mapping.
Definition isol_hor.h:900
Vector dn
Covariant derivative of the lapse with respect to the flat metric .
Definition isol_hor.h:937
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition single_hor.C:231
double regularisation(const Vector &shift_auto, const Vector &shift_comp, double ang_vel)
Corrects shift_auto in such a way that the total is equal to zero in the horizon,...
double regul
Intensity of the correction on the shift vector.
Definition isol_hor.h:912
Sym_tensor gamt_point
Time derivative of the 3-metric tilde.
Definition isol_hor.h:986
void n_comp_import(const Single_hor &comp)
Imports the part of N due to the companion hole comp .
Definition single_hor.C:390
Scalar n_auto
Lapse function .
Definition isol_hor.h:915
Scalar nn
Lapse function .
Definition isol_hor.h:921
Sym_tensor hh
Deviation metric.
Definition isol_hor.h:983
Sym_tensor aa
Components of the conformal representation of the traceless part of the extrinsic curvature:
Definition isol_hor.h:971
Scalar trK
Trace of the extrinsic curvature.
Definition isol_hor.h:989
void psi_comp_import(const Single_hor &comp)
Imports the part of due to the companion hole comp .
Definition single_hor.C:424
Metric tgam
3 metric tilde
Definition isol_hor.h:977
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Definition sym_tensor.C:349
Basic array class.
Definition tbl.h:161
Tensor handling.
Definition tensor.h:288
Values and coefficients of a (real-value) function.
Definition valeur.h:287
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
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition vector.C:381
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 diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:504
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition cmp_math.C:481
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 .
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition tensor.C:1014
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
Tensor derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition tensor.C:1002
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition app_hor.h:64