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