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 
24 char 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 
123 namespace Lorene {
124 void 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 
163  -2.*contract(hole1.dn, 0, hole1.psi_auto
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 
171  tmp_un = hole1.get_psi4()* contract(hole1.beta_auto, 0, hole1.trK.
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 
200  -2.*contract(hole2.dn, 0, hole2.psi_auto
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 
301 void 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 
442  hole1.set_der_0x0() ;
443  hole2.set_der_0x0() ;
444 
445  //set_hh_Samaya() ;
446 
447 }
448 
449 
450 // Resolution for shift with omega fixed.
451 // --------------------------------------
452 
453 void 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()
695  .set_base(*(hole1.mp.get_mg()->std_base_vect_cart()[0])) ;
696  omdsdp1.set(2).set_spectral_va()
697  .set_base(*(hole1.mp.get_mg()->std_base_vect_cart()[1])) ;
698  omdsdp1.set(3).set_spectral_va()
699  .set_base(*(hole1.mp.get_mg()->std_base_vect_cart()[2])) ;
700 
701  omdsdp1.annule_domain(nz-1) ;
702  omdsdp1.change_triad(hole1.mp.get_bvect_spher()) ;
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()
723  .set_base(*(hole2.mp.get_mg()->std_base_vect_cart()[0])) ;
724  omdsdp2.set(2).set_spectral_va()
725  .set_base(*(hole2.mp.get_mg()->std_base_vect_cart()[1])) ;
726  omdsdp2.set(3).set_spectral_va()
727  .set_base(*(hole2.mp.get_mg()->std_base_vect_cart()[2])) ;
728 
729  omdsdp2.annule_domain(nz-1) ;
730  omdsdp2.change_triad(hole2.mp.get_bvect_spher()) ;
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,
760  hole2.beta_auto, omega) ;
761  double diff_deux = hole2.regularisation (hole2.beta_auto,
762  hole1.beta_auto, omega) ;
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 }
Coord xa
Absolute x coordinate.
Definition: map.h:730
void beta_comp_import(const Single_hor &comp)
Imports the part of due to the companion hole comp.
Definition: single_hor.C:465
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
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...
double omega
Angular velocity.
Definition: isol_hor.h:1348
Vector dn
Covariant derivative of the lapse with respect to the flat metric .
Definition: isol_hor.h:937
Metric tgam
3 metric tilde
Definition: isol_hor.h:977
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition: single_hor.C:231
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
Scalar decouple
Function used to construct from the total .
Definition: isol_hor.h:1002
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:783
Lorene prototypes.
Definition: app_hor.h:64
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Scalar n_comp
Lapse function .
Definition: isol_hor.h:918
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, which should ensure the regularity of .
Definition: single_regul.C:59
Tensor derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Single_hor hole1
Black hole one.
Definition: isol_hor.h:1342
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
const Scalar & get_psi4() const
Conformal factor .
Definition: single_hor.C:288
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
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition: map.h:775
Sym_tensor aa_auto
Components of the conformal representation of the traceless part of the extrinsic curvature: ...
Definition: isol_hor.h:959
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
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
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
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
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:380
Sym_tensor aa
Components of the conformal representation of the traceless part of the extrinsic curvature: ...
Definition: isol_hor.h:971
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Definition: sym_tensor.C:349
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:481
const Valeur boundary_psi_app_hor() const
Neumann boundary condition for.
Definition: single_bound.C:69
Scalar psi
Conformal factor .
Definition: isol_hor.h:930
Scalar trK
Trace of the extrinsic curvature.
Definition: isol_hor.h:989
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 regul
Intensity of the correction on the shift vector.
Definition: isol_hor.h:912
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...
Sym_tensor hh
Deviation metric.
Definition: isol_hor.h:983
virtual const Connection & connect() const
Returns the connection.
Definition: metric.C:301
Vector beta_auto
Shift function .
Definition: isol_hor.h:944
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:381
Tensor handling.
Definition: tensor.h:288
Scalar n_auto
Lapse function .
Definition: isol_hor.h:915
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition: tensor.C:1002
Single_hor hole2
Black hole two.
Definition: isol_hor.h:1343
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
Vector dpsi
Covariant derivative of the conformal factor .
Definition: isol_hor.h:941
Coord ya
Absolute y coordinate.
Definition: map.h:731
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
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 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 ...
Map_af & mp
Affine mapping.
Definition: isol_hor.h:900
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:654
void n_comp_import(const Single_hor &comp)
Imports the part of N due to the companion hole comp .
Definition: single_hor.C:390
Sym_tensor gamt_point
Time derivative of the 3-metric tilde.
Definition: isol_hor.h:986
void psi_comp_import(const Single_hor &comp)
Imports the part of due to the companion hole comp .
Definition: single_hor.C:424
Scalar trK_point
Time derivative of the trace of the extrinsic curvature.
Definition: isol_hor.h:992
Scalar psi_auto
Conformal factor .
Definition: isol_hor.h:924
Basic array class.
Definition: tbl.h:161
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
Metric_flat ff
3 metric flat
Definition: isol_hor.h:980
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this.
Definition: scalar_deriv.C:390
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: tensor.C:926
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:223
Scalar nn
Lapse function .
Definition: isol_hor.h:921