LORENE
bound_hor.C
1 /*
2  * Method of class Isol_hor to compute boundary conditions
3  *
4  * (see file isol_hor.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2004 Jose Luis Jaramillo
10  * Francois Limousin
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License version 2
16  * as published by the Free Software Foundation.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 
31 /*
32  * $Id: bound_hor.C,v 1.37 2016/12/05 16:17:56 j_novak Exp $
33  * $Log: bound_hor.C,v $
34  * Revision 1.37 2016/12/05 16:17:56 j_novak
35  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
36  *
37  * Revision 1.36 2014/10/13 08:53:00 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.35 2014/10/06 15:13:10 j_novak
41  * Modified #include directives to use c++ syntax.
42  *
43  * Revision 1.34 2008/08/27 11:22:25 j_novak
44  * Minor modifications
45  *
46  * Revision 1.33 2008/08/19 06:42:00 j_novak
47  * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
48  * cast-type operations, and constant strings that must be defined as const char*
49  *
50  * Revision 1.32 2007/04/13 15:28:35 f_limousin
51  * Lots of improvements, generalisation to an arbitrary state of
52  * rotation, implementation of the spatial metric given by Samaya.
53  *
54  * Revision 1.31 2006/08/01 14:35:48 f_limousin
55  * Many small modifs
56  *
57  * Revision 1.30 2006/02/22 17:02:04 f_limousin
58  * Removal of warnings
59  *
60  * Revision 1.29 2006/02/22 16:29:55 jl_jaramillo
61  * corrections on the relaxation and boundary conditions
62  *
63  * Revision 1.27 2005/10/24 16:44:40 jl_jaramillo
64  * Cook boundary condition ans minot bound of kss
65  *
66  * Revision 1.26 2005/10/21 16:20:55 jl_jaramillo
67  * Version for the paper JaramL05
68  *
69  * Revision 1.25 2005/09/13 18:33:17 f_limousin
70  * New function vv_bound_cart_bin(double) for computing binaries with
71  * berlin condition for the shift vector.
72  * Suppress all the symy and asymy in the importations.
73  *
74  * Revision 1.24 2005/09/12 12:33:54 f_limousin
75  * Compilation Warning - Change of convention for the angular velocity
76  * Add Berlin boundary condition in the case of binary horizons.
77  *
78  * Revision 1.23 2005/07/08 13:15:23 f_limousin
79  * Improvements of boundary_vv_cart(), boundary_nn_lapl().
80  * Add a fonction to compute the departure of axisymmetry.
81  *
82  * Revision 1.22 2005/06/09 08:05:32 f_limousin
83  * Implement a new function sol_elliptic_boundary() and
84  * Vector::poisson_boundary(...) which solve the vectorial poisson
85  * equation (method 6) with an inner boundary condition.
86  *
87  * Revision 1.21 2005/05/12 14:48:07 f_limousin
88  * New boundary condition for the lapse : boundary_nn_lapl().
89  *
90  * Revision 1.20 2005/04/29 14:04:17 f_limousin
91  * Implementation of boundary_vv_x (y,z) for binary black holes.
92  *
93  * Revision 1.19 2005/04/19 16:40:51 jl_jaramillo
94  * change of sign of ang_vel in vv_bound_cart. Convention of Phys. Rep.
95  *
96  * Revision 1.18 2005/04/08 12:16:52 f_limousin
97  * Function set_psi(). And dependance in phi.
98  *
99  * Revision 1.17 2005/04/02 15:49:21 f_limousin
100  * New choice (Lichnerowicz) for aaquad. New member data nz.
101  *
102  * Revision 1.16 2005/03/22 13:25:36 f_limousin
103  * Small changes. The angular velocity and A^{ij} are computed
104  * with a differnet sign.
105  *
106  * Revision 1.15 2005/03/10 10:19:42 f_limousin
107  * Add the regularisation of the shift in the case of a single black hole
108  * and lapse zero on the horizon.
109  *
110  * Revision 1.14 2005/03/03 10:00:55 f_limousin
111  * The funtions beta_boost_x() and beta_boost_z() have been added.
112  *
113  * Revision 1.13 2005/02/24 17:21:04 f_limousin
114  * Suppression of the function beta_bound_cart() and direct computation
115  * of boundary_beta_x, y and z.
116  *
117  * Revision 1.12 2004/12/31 15:34:37 f_limousin
118  * Modifications to avoid warnings
119  *
120  * Revision 1.11 2004/12/22 18:15:16 f_limousin
121  * Many different changes.
122  *
123  * Revision 1.10 2004/11/24 19:30:58 jl_jaramillo
124  * Berlin boundary conditions vv_bound_cart
125  *
126  * Revision 1.9 2004/11/18 09:49:44 jl_jaramillo
127  * Some new conditions for the shift (Neumann + Dirichlet)
128  *
129  * Revision 1.8 2004/11/05 10:52:26 f_limousin
130  * Replace double aa by double cc in argument of boundary_beta_x
131  * boundary_beta_y and boundary_beta_z to avoid warnings.
132  *
133  * Revision 1.7 2004/10/29 15:42:14 jl_jaramillo
134  * Static shift boundary conbdition
135  *
136  * Revision 1.6 2004/10/01 16:46:51 f_limousin
137  * Added a pure Dirichlet boundary condition
138  *
139  * Revision 1.5 2004/09/28 16:06:41 f_limousin
140  * Correction of an error when taking the bases of the boundary
141  * condition for the shift.
142  *
143  * Revision 1.4 2004/09/17 13:36:23 f_limousin
144  * Add some new boundary conditions
145  *
146  * Revision 1.2 2004/09/09 16:53:49 f_limousin
147  * Add the two lines $Id: bound_hor.C,v 1.37 2016/12/05 16:17:56 j_novak Exp $Log: for CVS.
148  *
149  *
150  * $Header: /cvsroot/Lorene/C++/Source/Isol_hor/bound_hor.C,v 1.37 2016/12/05 16:17:56 j_novak Exp $
151  *
152  */
153 
154 // C++ headers
155 #include "headcpp.h"
156 
157 // C headers
158 #include <cstdlib>
159 #include <cassert>
160 
161 // Lorene headers
162 #include "time_slice.h"
163 #include "isol_hor.h"
164 #include "metric.h"
165 #include "evolution.h"
166 #include "unites.h"
167 #include "graphique.h"
168 #include "utilitaires.h"
169 #include "param.h"
170 
171 
172 // Dirichlet boundary condition for Psi
173 //-------------------------------------
174 // ONE HAS TO GUARANTEE THAT BETA IS NOT ZERO, BUT IT IS PROPORTIONAL TO THE RADIAL VECTOR
175 
176 namespace Lorene {
178 
179  Scalar tmp = - 6 * contract(beta(), 0, psi().derive_cov(ff), 0) ;
180  tmp = tmp / (contract(beta().derive_cov(ff), 0, 1) - nn() * trk() ) - 1 ;
181 
182  // We have substracted 1, since we solve for zero condition at infinity
183  //and then we add 1 to the solution
184 
185  Valeur psi_bound (mp.get_mg()->get_angu() ) ;
186 
187  int nnp = mp.get_mg()->get_np(1) ;
188  int nnt = mp.get_mg()->get_nt(1) ;
189 
190  psi_bound = 1 ;
191 
192  for (int k=0 ; k<nnp ; k++)
193  for (int j=0 ; j<nnt ; j++)
194  psi_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
195 
196  psi_bound.std_base_scal() ;
197 
198  return psi_bound ;
199 
200 }
201 
202 
203 // Neumann boundary condition for Psi
204 //-------------------------------------
205 
207 
208  // Introduce 2-trace gamma tilde dot
209  Scalar tmp = - 1./ 6. * psi() * (beta().divergence(ff) - nn() * trk() )
210  - beta()(2)* psi().derive_cov(ff)(2) - beta()(3)* psi().derive_cov(ff)(3) ;
211 
212  tmp = tmp / beta()(1) ;
213 
214  // in this case you don't have to substract any value
215 
216  Valeur psi_bound (mp.get_mg()->get_angu() ) ;
217 
218  int nnp = mp.get_mg()->get_np(1) ;
219  int nnt = mp.get_mg()->get_nt(1) ;
220 
221  psi_bound = 1 ;
222 
223  for (int k=0 ; k<nnp ; k++)
224  for (int j=0 ; j<nnt ; j++)
225  psi_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
226 
227  psi_bound.std_base_scal() ;
228 
229  return psi_bound ;
230 
231 }
232 
233 
235 
236  Scalar tmp = psi() * psi() * psi() * trk()
237  - contract(k_dd(), 0, 1, tradial_vect_hor() * tradial_vect_hor(), 0, 1)
238  / psi()
239  - 4.* contract(tradial_vect_hor(), 0, psi().derive_cov(ff), 0) ;
240 
241  // rho = 1 is the standart case.
242  double rho = 1. ;
243  tmp += (tradial_vect_hor().divergence(ff)) * (rho - 1.)*psi() ;
244 
245  tmp = tmp / (rho * tradial_vect_hor().divergence(ff)) - 1. ;
246 
247  tmp.std_spectral_base() ;
248  tmp.inc_dzpuis(2) ;
249 
250  // We have substracted 1, since we solve for zero condition at infinity
251  //and then we add 1 to the solution
252 
253  Valeur psi_bound (mp.get_mg()->get_angu() ) ;
254 
255  int nnp = mp.get_mg()->get_np(1) ;
256  int nnt = mp.get_mg()->get_nt(1) ;
257 
258  psi_bound = 1 ;
259 
260  for (int k=0 ; k<nnp ; k++)
261  for (int j=0 ; j<nnt ; j++)
262  psi_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
263 
264  psi_bound.std_base_scal() ;
265 
266  return psi_bound ;
267 
268 }
269 
271 
272  Scalar tmp = psi() * psi() * psi() * trk()
273  - contract(k_dd(), 0, 1, tradial_vect_hor() * tradial_vect_hor(), 0, 1)
274  / psi()
276  - 4 * ( tradial_vect_hor()(2) * psi().derive_cov(ff)(2)
277  + tradial_vect_hor()(3) * psi().derive_cov(ff)(3) ) ;
278 
279  tmp = tmp / (4 * tradial_vect_hor()(1)) ;
280 
281  // in this case you don't have to substract any value
282 
283  Valeur psi_bound (mp.get_mg()->get_angu() ) ;
284 
285  int nnp = mp.get_mg()->get_np(1) ;
286  int nnt = mp.get_mg()->get_nt(1) ;
287 
288  psi_bound = 1 ;
289 
290  for (int k=0 ; k<nnp ; k++)
291  for (int j=0 ; j<nnt ; j++)
292  psi_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
293 
294  psi_bound.std_base_scal() ;
295 
296  return psi_bound ;
297 
298 }
299 
301 
302  int nnp = mp.get_mg()->get_np(1) ;
303  int nnt = mp.get_mg()->get_nt(1) ;
304 
305  Valeur psi_bound (mp.get_mg()->get_angu()) ;
306 
307  psi_bound = 1 ; // Juste pour affecter dans espace des configs ;
308 
309 // if (psi_comp_evol.is_known(jtime)) {
310 // for (int k=0 ; k<nnp ; k++)
311 // for (int j=0 ; j<nnt ; j++)
312 // psi_bound.set(0, k, j, 0) = - 0.5/radius*(psi_auto().val_grid_point(1, k, j, 0) + psi_comp().val_grid_point(1, k, j, 0)) ;
313 // }
314 // else {
315  for (int k=0 ; k<nnp ; k++)
316  for (int j=0 ; j<nnt ; j++)
317  psi_bound.set(0, k, j, 0) = - 0.5/radius*psi().val_grid_point(1, k, j, 0) ;
318 // }
319 
320 
321  psi_bound.std_base_scal() ;
322 
323  return psi_bound ;
324 
325 }
326 
328 
329  Scalar rho (mp) ;
330  rho = 5. ;
331  rho.std_spectral_base() ;
332 
333 
334  Scalar tmp(mp) ;
335  // tmp = nn()+1. - 1 ;
336 
337 
338  //Scalar b_tilde_tmp = b_tilde() ;
339  //b_tilde_tmp.set_domain(0) = 1. ;
340  //tmp = pow(nn()/b_tilde_tmp, 0.5) ;
341 
342 
343  tmp = 1.5 ;
344  tmp.std_spectral_base() ;
345 
346  //tmp = 1/ (2* nn()) ;
347 
348  tmp = (tmp + rho * psi())/(1 + rho) ;
349 
350  tmp = tmp - 1. ;
351 
352 
353 
354  Valeur psi_bound (mp.get_mg()->get_angu()) ;
355 
356  psi_bound = 1 ; // Juste pour affecter dans espace des configs ;
357 
358  int nnp = mp.get_mg()->get_np(1) ;
359  int nnt = mp.get_mg()->get_nt(1) ;
360 
361  for (int k=0 ; k<nnp ; k++)
362  for (int j=0 ; j<nnt ; j++)
363  psi_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
364 
365  psi_bound.std_base_scal() ;
366 
367  return psi_bound ;
368 
369 }
370 
371 // Dirichlet boundary condition on nn using the extrinsic curvature
372 // (No time evolution taken into account! Make this)
373 //--------------------------------------------------------------------------
375 
376  Scalar tmp(mp) ;
377  double rho = 0. ;
378 
379  // Scalar kk_rr = 0.8*psi().derive_cov(mp.flat_met_spher())(1) / psi() ;
380  Scalar kk_rr = contract( gam().radial_vect() * gam().radial_vect(), 0, 1
381  , k_dd(), 0, 1 ) ;
382 
383  Scalar k_kerr (mp) ;
384  k_kerr = 0.1 ;//1.*kappa_hor() ;
385  k_kerr.std_spectral_base() ;
386  k_kerr.inc_dzpuis(2) ;
387 
388  Scalar temp (rho*nn()) ;
389  temp.inc_dzpuis(2) ;
390 
391  tmp = contract(gam().radial_vect(), 0, nn().derive_cov(ff), 0) + temp
392  - k_kerr ;
393 
394  tmp = tmp / (kk_rr + rho) - 1;
395 
396  Scalar diN (contract(gam().radial_vect(), 0, nn().derive_cov(ff), 0)) ;
397  cout << "k_kerr = " << k_kerr.val_grid_point(1, 0, 0, 0) << endl ;
398  cout << "D^i N = " << diN.val_grid_point(1, 0, 0, 0) << endl ;
399  cout << "kss = " << kk_rr.val_grid_point(1, 0, 0, 0) << endl ;
400 
401  // We have substracted 1, since we solve for zero condition at infinity
402  //and then we add 1 to the solution
403 
404  int nnp = mp.get_mg()->get_np(1) ;
405  int nnt = mp.get_mg()->get_nt(1) ;
406 
407  Valeur nn_bound (mp.get_mg()->get_angu()) ;
408 
409  nn_bound = 1 ; // Juste pour affecter dans espace des configs ;
410 
411 
412  for (int k=0 ; k<nnp ; k++)
413  for (int j=0 ; j<nnt ; j++)
414  nn_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
415 
416  nn_bound.std_base_scal() ;
417 
418  return nn_bound ;
419 
420 }
421 
422 
423 const Valeur Isol_hor::boundary_nn_Neu_kk(int step) const {
424 
425  const Vector& dnnn = nn().derive_cov(ff) ;
426  double rho = 5. ;
427 
428  step = 100 ; // Truc bidon pour eviter warning
429 
430  Scalar kk_rr = contract( gam().radial_vect() * gam().radial_vect(), 0, 1
431  , k_dd(), 0, 1 ) ;
432 
433  Scalar k_kerr (mp) ;
434  k_kerr = kappa_hor() ;
435  //k_kerr = 0.06 ;
436  k_kerr.std_spectral_base() ;
437  k_kerr.inc_dzpuis(2) ;
438 
439  Scalar k_hor (mp) ;
440  k_hor = kappa_hor() ;
441  k_hor.std_spectral_base() ;
442 
443  // Vector beta_tilde_d (beta().down(0, met_gamt)) ;
444  //Scalar naass = 1./2. * contract( met_gamt.radial_vect() * met_gamt.radial_vect(), 0, 1, beta_tilde_d.ope_killing_conf(met_gamt) , 0, 1 ) ;
445 
447  0, 1, aa_auto().up_down(met_gamt), 0, 1) ;
448 
449  Scalar traceK = 1./3. * nn() * trk() *
451  , met_gamt.cov() , 0, 1 ) ;
452 
453  Scalar sdN (contract(gam().radial_vect(), 0, nn().derive_cov(ff), 0)) ;
454 
455 
456 
457  Scalar tmp = psi() * psi() * ( k_kerr + naass + 1.* traceK)
458  - met_gamt.radial_vect()(2) * dnnn(2)
459  - met_gamt.radial_vect()(3) * dnnn(3)
460  + ( rho - 1 ) * met_gamt.radial_vect()(1) * dnnn(1) ;
461 
462 
463  tmp = tmp / (rho * met_gamt.radial_vect()(1)) ;
464 
465 
466 
467 
468  Scalar rhs ( sdN - nn() * kk_rr) ;
469  cout << "kappa_pres = " << k_kerr.val_grid_point(1, 0, 0, 0) << endl ;
470  cout << "kappa_hor = " << k_hor.val_grid_point(1, 0, 0, 0) << endl ;
471  cout << "sDN = " << sdN.val_grid_point(1, 0, 0, 0) << endl ;
472 
473  // in this case you don't have to substract any value
474 
475  int nnp = mp.get_mg()->get_np(1) ;
476  int nnt = mp.get_mg()->get_nt(1) ;
477 
478  Valeur nn_bound (mp.get_mg()->get_angu()) ;
479 
480  nn_bound = 1 ; // Juste pour affecter dans espace des configs ;
481 
482  for (int k=0 ; k<nnp ; k++)
483  for (int j=0 ; j<nnt ; j++)
484  nn_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
485 
486  nn_bound.std_base_scal() ;
487 
488  return nn_bound ;
489 
490 }
491 
493 
494  const Vector& dnnn = nn().derive_cov(ff) ;
495  double rho = 1. ;
496 
497 
498  Scalar kk_rr = contract( gam().radial_vect() * gam().radial_vect(), 0, 1
499  , k_dd(), 0, 1 ) ;
500 
501  Sym_tensor qq_uu = gam_uu() - gam().radial_vect() * gam().radial_vect() ;
502 
503  // Free function L_theta = h
504  //--------------------------
505  Scalar L_theta (mp) ;
506  L_theta = .0;
507  L_theta.std_spectral_base() ;
508  L_theta.inc_dzpuis(4) ;
509 
510  //Divergence of Omega
511  //-------------------
512  Vector hom1 = nn().derive_cov(met_gamt)/nn() ;
513  hom1 = contract(qq_uu.down(1, gam()), 0, hom1, 0) ;
514 
515  Vector hom2 = -contract( k_dd(), 1, gam().radial_vect() , 0) ;
516  hom2 = contract(qq_uu.down(1, gam()), 0, hom2, 0) ;
517 
518  Vector hom = hom1 + hom2 ;
519 
520  Scalar div_omega = 1.*contract( qq_uu, 0, 1, (1.*hom1 + 1.*hom2).derive_cov(gam()), 0, 1) ;
521  div_omega.inc_dzpuis() ;
522 
523  //---------------------
524 
525 
526  // Two-dimensional Ricci scalar
527  //-----------------------------
528 
529  Scalar rr (mp) ;
530  rr = mp.r ;
531 
532  Scalar Ricci_conf(mp) ;
533  Ricci_conf = 2 / rr / rr ;
534  Ricci_conf.std_spectral_base() ;
535 
536  Scalar Ricci(mp) ;
537  Scalar log_psi (log(psi())) ;
538  log_psi.std_spectral_base() ;
539  Ricci = pow(psi(), -4) * (Ricci_conf - 4*log_psi.lapang())/rr /rr ;
540  Ricci.std_spectral_base() ;
541  Ricci.inc_dzpuis(4) ;
542  //-------------------------------
543 
544  Scalar theta_k = -1/(2*nn()) * (gam().radial_vect().divergence(gam()) -
545  kk_rr + trk() ) ;
546 
547  int nnp = mp.get_mg()->get_np(1) ;
548  int nnt = mp.get_mg()->get_nt(1) ;
549  /*
550  for (int k=0 ; k<nnp ; k++)
551  for (int j=0 ; j<nnt ; j++){
552  cout << theta_k.val_grid_point(1, k, j, 0) << endl ;
553  cout << nn().val_grid_point(1, k, j, 0) << endl ;
554  }
555  */
556 
557 
558 
559  //Source
560  //------
561  Scalar source = div_omega + contract( qq_uu, 0, 1, hom * hom , 0, 1) - 0.5 * Ricci - L_theta;
562  source = source / theta_k ;
563 
564  Scalar tmp = ( source + nn() * kk_rr + rho * contract(gam().radial_vect(), 0,
565  nn().derive_cov(gam()), 0))/(1+rho)
566  - gam().radial_vect()(2) * dnnn(2) - gam().radial_vect()(3) * dnnn(3) ;
567 
568  tmp = tmp / gam().radial_vect()(1) ;
569 
570  // in this case you don't have to substract any value
571 
572  Valeur nn_bound (mp.get_mg()->get_angu()) ;
573 
574  nn_bound = 1 ; // Juste pour affecter dans espace des configs ;
575 
576 
577  for (int k=0 ; k<nnp ; k++)
578  for (int j=0 ; j<nnt ; j++)
579  nn_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
580 
581  nn_bound.std_base_scal() ;
582 
583  return nn_bound ;
584 
585 }
586 
587 
588 
589 const Valeur Isol_hor::boundary_nn_Dir_eff(double cc)const {
590 
591  Scalar tmp(mp) ;
592 
593  tmp = - nn().derive_cov(ff)(1) ;
594 
595  // rho = 1 is the standart case
596  double rho = 1. ;
597  tmp.dec_dzpuis(2) ;
598  tmp += cc * (rho -1)*nn() ;
599  tmp = tmp / (rho*cc) ;
600 
601  tmp = tmp - 1. ;
602 
603  // We have substracted 1, since we solve for zero condition at infinity
604  //and then we add 1 to the solution
605 
606  int nnp = mp.get_mg()->get_np(1) ;
607  int nnt = mp.get_mg()->get_nt(1) ;
608 
609  Valeur nn_bound (mp.get_mg()->get_angu()) ;
610 
611  nn_bound = 1 ; // Juste pour affecter dans espace des configs ;
612 
613  for (int k=0 ; k<nnp ; k++)
614  for (int j=0 ; j<nnt ; j++)
615  nn_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
616 
617  nn_bound.std_base_scal() ;
618 
619  return nn_bound ;
620 
621 }
622 
623 
624 
625 const Valeur Isol_hor::boundary_nn_Neu_eff(double cc)const {
626 
627  Scalar tmp = - cc * nn() ;
628  // Scalar tmp = - nn()/psi()*psi().dsdr() ;
629 
630  // in this case you don't have to substract any value
631 
632  int nnp = mp.get_mg()->get_np(1) ;
633  int nnt = mp.get_mg()->get_nt(1) ;
634 
635  Valeur nn_bound (mp.get_mg()->get_angu()) ;
636 
637  nn_bound = 1 ; // Juste pour affecter dans espace des configs ;
638 
639  for (int k=0 ; k<nnp ; k++)
640  for (int j=0 ; j<nnt ; j++)
641  nn_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
642 
643  nn_bound.std_base_scal() ;
644 
645  return nn_bound ;
646 
647 }
648 
649 
650 const Valeur Isol_hor::boundary_nn_Dir(double cc)const {
651 
652  Scalar rho(mp);
653  rho = 0. ; // 0 is the standard case
654 
655  Scalar tmp(mp) ;
656  tmp = cc ;
657 
658  /*
659  if (b_tilde().val_grid_point(1, 0, 0, 0) < 0.08)
660  tmp = 0.25 ;
661  else {
662  cout << "OK" << endl ;
663  // des_profile(nn(), 0, 20, M_PI/2, M_PI) ;
664  rho = 5. ;
665  tmp = b_tilde()*psi()*psi() ;
666  }
667  */
668 
669  //tmp = 1./(2*psi()) ;
670  // tmp = - psi() * nn().dsdr() / (psi().dsdr()) ;
671 
672  // We have substracted 1, since we solve for zero condition at infinity
673  //and then we add 1 to the solution
674 
675  tmp = (tmp + rho * nn())/(1 + rho) ;
676 
677  tmp = tmp - 1 ;
678 
679  int nnp = mp.get_mg()->get_np(1) ;
680  int nnt = mp.get_mg()->get_nt(1) ;
681 
682  Valeur nn_bound (mp.get_mg()->get_angu()) ;
683 
684  nn_bound = 1 ; // Juste pour affecter dans espace des configs ;
685 
686  for (int k=0 ; k<nnp ; k++)
687  for (int j=0 ; j<nnt ; j++)
688  nn_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
689 
690  nn_bound.std_base_scal() ;
691 
692  return nn_bound ;
693 
694 }
695 
697 
698  int nnp = mp.get_mg()->get_np(1) ;
699  int nnt = mp.get_mg()->get_nt(1) ;
700 
701 
702  //Preliminaries
703  //-------------
704 
705  //Metrics on S^2
706  //--------------
707 
708  Sym_tensor q_uu = gam_uu() - gam().radial_vect() * gam().radial_vect() ;
709  Sym_tensor q_dd = q_uu.up_down(gam()) ;
710 
711  Scalar det_q = q_dd(2,2) * q_dd(3,3) - q_dd(2,3) * q_dd(3,2) ;
712  Scalar square_q (pow(det_q, 0.5)) ;
713  square_q.std_spectral_base() ;
714 
715  Sym_tensor qhat_uu = square_q * q_uu ;
716 
717  Sym_tensor ffr_uu = ff.con() - ff.radial_vect() * ff.radial_vect() ;
718  Sym_tensor ffr_dd = ff.cov() - ff.radial_vect().down(0, ff) * ff.radial_vect().down(0,ff) ;
719 
720  Sym_tensor h_uu = qhat_uu - ffr_uu ;
721 
722  //------------------
723 
724 
725  // Source of the "round" laplacian:
726  //---------------------------------
727 
728  //Source 1: "Divergence" term
729  //---------------------------
730  Vector kqs = contract(k_dd(), 1, gam().radial_vect(), 0 ) ;
731  kqs = contract( q_uu.down(0, gam()), 1, kqs, 0) ;
732  Scalar div_kqs = contract( q_uu, 0, 1, kqs.derive_cov(gam()), 0, 1) ;
733 
734  Scalar source1 = div_kqs ;
735  source1 *= square_q ;
736 
737 
738  // Source 2: correction term
739  //--------------------------
740 
741  Vector corr = - contract(h_uu, 1, nn().derive_cov(ff), 0) / nn() ;
742  Scalar source2 = contract( ffr_dd, 0, 1, corr.derive_con(ff), 0, 1 ) ;
743 
744 
745  // Source 3: (physical) divergence of Omega
746  //-----------------------------------------
747 
748  Scalar div_omega(mp) ;
749 
750  //Source from (L_l\theta_k=0)
751 
752  Scalar L_theta (mp) ;
753  L_theta = 0. ;
754  L_theta.std_spectral_base() ;
755 
756  Scalar kk_rr = contract( gam().radial_vect() * gam().radial_vect(), 0, 1
757  , k_dd(), 0, 1 ) ;
758 
759  // Scalar kappa (mp) ;
760  // kappa = kappa_hor() ;
761  // kappa.std_spectral_base() ;
762  // kappa.set_dzpuis(2) ;
763 
764 
765  Scalar kappa = contract(gam().radial_vect(), 0, nn().derive_cov(gam()), 0) - nn() * kk_rr ;
766  Scalar theta_k = -1/(2*nn()) * (gam().radial_vect().divergence(gam()) -
767  kk_rr + trk() ) ;
768 
769  Sym_tensor qqq = gam_uu() - gam().radial_vect() * gam().radial_vect() ;
770 
771  Vector hom = nn().derive_cov(met_gamt) - contract( k_dd(), 1, gam().radial_vect() , 0) ;
772  hom = contract(qqq.down(1, gam()), 0, hom, 0) ;
773 
774  Scalar rr(mp) ;
775  rr = mp.r ;
776 
777  Scalar Ricci_conf = 2 / rr /rr ;
778  Ricci_conf.std_spectral_base() ;
779 
780  Scalar log_psi (log(psi())) ;
781  log_psi.std_spectral_base() ;
782  Scalar Ricci = pow(psi(), -4) * (Ricci_conf - 4*log_psi.lapang())/rr /rr ;
783  Ricci.std_spectral_base() ;
784  Ricci.inc_dzpuis(4) ;
785 
786 
787  div_omega = L_theta - contract(qqq, 0, 1, hom * hom, 0, 1) + 0.5 * Ricci
788  + theta_k * kappa ;
789  div_omega.dec_dzpuis() ;
790 
791  //End of Source from (L_l\theta_k=0)
792  //----------------------------------
793 
794  div_omega = 0. ;
795  // div_omega = 0.01*log(1/(2*psi())) ;
796  div_omega.std_spectral_base() ;
797  div_omega.inc_dzpuis(3) ;
798 
799 
800  //Construction of source3 (correction of div_omega by the square root of the determinant)
801  Scalar source3 = div_omega ;
802  source3 *= square_q ;
803 
804 
805  //"Correction" to the div_omega term (no Y_00 component must be present)
806 
807  double corr_const = mp.integrale_surface(source3, radius + 1e-15)/(4. * M_PI) ;
808  cout << "Constant part of div_omega = " << corr_const <<endl ;
809 
810  Scalar corr_div_omega (mp) ;
811  corr_div_omega = corr_const ;
812  corr_div_omega.std_spectral_base() ;
813  corr_div_omega.set_dzpuis(3) ;
814 
815  source3 -= corr_div_omega ;
816 
817 
818 
819  //Source
820  //------
821 
822  Scalar source = source1 + source2 + source3 ;
823 
824 
825 
826  //Resolution of round angular laplacian
827  //-------------------------------------
828 
829  Scalar source_tmp = source ;
830 
831  Scalar logn (mp) ;
832  logn = source.poisson_angu() ;
833 
834  double cc = 0.2 ; // Integration constant
835 
836  Scalar lapse (mp) ;
837  lapse = (exp(logn)*cc) ;
838  lapse.std_spectral_base() ;
839 
840 
841  //Test of Divergence of Omega
842  //---------------------------
843 
844  ofstream err ;
845  err.open ("err_laplBC.d", ofstream::app ) ;
846 
847  hom = nn().derive_cov(gam())/nn()
848  - contract( k_dd(), 1, gam().radial_vect() , 0) ;
849  hom = contract(q_uu.down(1, gam()), 0, hom, 0) ;
850 
851  Scalar div_omega_test = contract( q_uu, 0, 1, hom.derive_cov(gam()), 0, 1) ;
852 
853 
854  Scalar err_lapl = div_omega_test - div_omega ;
855 
856  double max_err = err_lapl.val_grid_point(1, 0, 0, 0) ;
857  double min_err = err_lapl.val_grid_point(1, 0, 0, 0) ;
858 
859  for (int k=0 ; k<nnp ; k++)
860  for (int j=0 ; j<nnt ; j++){
861  if (err_lapl.val_grid_point(1, k, j, 0) > max_err)
862  max_err = err_lapl.val_grid_point(1, k, j, 0) ;
863  if (err_lapl.val_grid_point(1, k, j, 0) < min_err)
864  min_err = err_lapl.val_grid_point(1, k, j, 0) ;
865  }
866 
867  err << mer << " " << max_err << " " << min_err << endl ;
868 
869  err.close() ;
870 
871 
872 
873  //Construction of the Valeur
874  //--------------------------
875 
876  lapse = lapse - 1. ;
877 
878  // int nnp = mp.get_mg()->get_np(1) ;
879  // int nnt = mp.get_mg()->get_nt(1) ;
880 
881  Valeur nn_bound (mp.get_mg()->get_angu()) ;
882 
883  nn_bound = 1 ; // Juste pour affecter dans espace des configs ;
884 
885  for (int k=0 ; k<nnp ; k++)
886  for (int j=0 ; j<nnt ; j++)
887  nn_bound.set(0, k, j, 0) = lapse.val_grid_point(1, k, j, 0) ;
888 
889  nn_bound.std_base_scal() ;
890 
891  return nn_bound ;
892 
893 }
894 
895 
896 // Component r of boundary value of beta (using expression in terms
897 // of radial vector)
898 // --------------------------------------
899 
901 
902  Scalar tmp (mp) ;
903 
904  tmp = nn() * radial_vect_hor()(1) ;
905 
906  Base_val base_tmp (radial_vect_hor()(1).get_spectral_va().base) ;
907 
908  int nnp = mp.get_mg()->get_np(1) ;
909  int nnt = mp.get_mg()->get_nt(1) ;
910 
911  Valeur bnd_beta_r (mp.get_mg()->get_angu()) ;
912 
913  bnd_beta_r = 1. ; // Juste pour affecter dans espace des configs ;
914 
915  for (int k=0 ; k<nnp ; k++)
916  for (int j=0 ; j<nnt ; j++)
917  bnd_beta_r.set(1, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
918 
919  bnd_beta_r.set_base_r(0, base_tmp.get_base_r(0)) ;
920  for (int l=0 ; l<(*mp.get_mg()).get_nzone()-1 ; l++)
921  bnd_beta_r.set_base_r(l, base_tmp.get_base_r(l)) ;
922  bnd_beta_r.set_base_r((*mp.get_mg()).get_nzone()-1, base_tmp.get_base_r((*mp.get_mg()).get_nzone()-1)) ;
923  bnd_beta_r.set_base_t(tmp.get_spectral_va().get_base().get_base_t(1)) ;
924  bnd_beta_r.set_base_p(tmp.get_spectral_va().get_base().get_base_p(1)) ;
925 
926 // bnd_beta_r.set_base(*(mp.get_mg()->std_base_vect_spher()[0])) ;
927 
928  cout << "norme de lim_vr" << endl << norme(bnd_beta_r) << endl ;
929  cout << "bases" << endl << bnd_beta_r.base << endl ;
930 
931  return bnd_beta_r ;
932 
933 
934 }
935 
936 
937 // Component theta of boundary value of beta (using expression in terms
938 // of radial vector)
939 // ------------------------------------------
940 
942 
943  Scalar tmp(mp) ;
944 
945  tmp = nn() * radial_vect_hor()(2) ;
946  Base_val base_tmp (radial_vect_hor()(2).get_spectral_va().base) ;
947 
948 
949  int nnp = mp.get_mg()->get_np(1) ;
950  int nnt = mp.get_mg()->get_nt(1) ;
951 
952  Valeur bnd_beta_theta (mp.get_mg()->get_angu()) ;
953 
954  bnd_beta_theta = 1. ; // Juste pour affecter dans espace des configs ;
955 
956  for (int k=0 ; k<nnp ; k++)
957  for (int j=0 ; j<nnt ; j++)
958  bnd_beta_theta.set(1, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
959 
960 // bnd_beta_theta.set_base(*(mp.get_mg()->std_base_vect_spher()[1])) ;
961 
962  bnd_beta_theta.set_base_r(0, base_tmp.get_base_r(0)) ;
963  for (int l=0 ; l<(*mp.get_mg()).get_nzone()-1 ; l++)
964  bnd_beta_theta.set_base_r(l, base_tmp.get_base_r(l)) ;
965  bnd_beta_theta.set_base_r((*mp.get_mg()).get_nzone()-1, base_tmp.get_base_r((*mp.get_mg()).get_nzone()-1)) ;
966  bnd_beta_theta.set_base_t(base_tmp.get_base_t(1)) ;
967  bnd_beta_theta.set_base_p(base_tmp.get_base_p(1)) ;
968 
969  cout << "bases" << endl << bnd_beta_theta.base << endl ;
970 
971 
972  return bnd_beta_theta ;
973 
974 
975 }
976 
977 // Component phi of boundary value of beta (using expression in terms
978 // of radial vector)
979 // -----------------------------------------------------------
980 
981 const Valeur Isol_hor::boundary_beta_phi(double om)const {
982 
983  Scalar tmp (mp) ;
984 
985  Scalar ang_vel(mp) ;
986  ang_vel = om ;
987  ang_vel.std_spectral_base() ;
988  ang_vel.mult_rsint() ;
989 
990  tmp = nn() * radial_vect_hor()(3) - ang_vel ;
991  Base_val base_tmp (ang_vel.get_spectral_va().base) ;
992 
993  int nnp = mp.get_mg()->get_np(1) ;
994  int nnt = mp.get_mg()->get_nt(1) ;
995 
996  Valeur bnd_beta_phi (mp.get_mg()->get_angu()) ;
997 
998  bnd_beta_phi = 1. ; // Juste pour affecter dans espace des configs ;
999 
1000  for (int k=0 ; k<nnp ; k++)
1001  for (int j=0 ; j<nnt ; j++)
1002  bnd_beta_phi.set(1, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
1003 
1004  for (int l=0 ; l<(*mp.get_mg()).get_nzone()-1 ; l++)
1005 
1006 // bnd_beta_phi.set_base(*(mp.get_mg()->std_base_vect_spher()[2])) ;
1007 
1008  bnd_beta_phi.set_base_r(0, base_tmp.get_base_r(0)) ;
1009  for (int l=0 ; l<(*mp.get_mg()).get_nzone()-1 ; l++)
1010  bnd_beta_phi.set_base_r(l, base_tmp.get_base_r(l)) ;
1011  bnd_beta_phi.set_base_r((*mp.get_mg()).get_nzone()-1, base_tmp.get_base_r((*mp.get_mg()).get_nzone()-1)) ;
1012  bnd_beta_phi.set_base_t(base_tmp.get_base_t(1)) ;
1013  bnd_beta_phi.set_base_p(base_tmp.get_base_p(1)) ;
1014 
1015 // cout << "bound beta_phi" << endl << bnd_beta_phi << endl ;
1016 
1017  return bnd_beta_phi ;
1018 
1019 }
1020 
1021 // Component x of boundary value of beta
1022 //--------------------------------------
1023 
1024 const Valeur Isol_hor:: boundary_beta_x(double om)const {
1025 
1026  // Les alignemenents pour le signe des CL.
1027  double orientation = mp.get_rot_phi() ;
1028  assert ((orientation == 0) || (orientation == M_PI)) ;
1029  int aligne = (orientation == 0) ? 1 : -1 ;
1030 
1031  int nnp = mp.get_mg()->get_np(1) ;
1032  int nnt = mp.get_mg()->get_nt(1) ;
1033 
1034  Vector tmp_vect = nn() * gam().radial_vect() ;
1035  tmp_vect.change_triad(mp.get_bvect_cart() ) ;
1036 
1037  //Isol_hor boundary conditions
1038 
1039  Valeur lim_x (mp.get_mg()->get_angu()) ;
1040 
1041  lim_x = 1 ; // Juste pour affecter dans espace des configs ;
1042 
1043  Mtbl ya_mtbl (mp.get_mg()) ;
1044  ya_mtbl.set_etat_qcq() ;
1045  ya_mtbl = mp.ya ;
1046 
1047  Scalar cosp (mp) ;
1048  cosp = mp.cosp ;
1049  Scalar cost (mp) ;
1050  cost = mp.cost ;
1051  Scalar sinp (mp) ;
1052  sinp = mp.sinp ;
1053  Scalar sint (mp) ;
1054  sint = mp.sint ;
1055 
1056  Scalar dep_phi (mp) ;
1057  dep_phi = 0.0*sint*cosp ;
1058 
1059  for (int k=0 ; k<nnp ; k++)
1060  for (int j=0 ; j<nnt ; j++)
1061  lim_x.set(0, k, j, 0) = aligne * om * ya_mtbl(1, k, j, 0) * (1 +
1062  dep_phi.val_grid_point(1, k, j, 0))
1063  + tmp_vect(1).val_grid_point(1, k, j, 0) ;
1064 
1065  lim_x.set_base(*(mp.get_mg()->std_base_vect_cart()[0])) ;
1066 
1067  return lim_x ;
1068 }
1069 
1070 
1071 // Component y of boundary value of beta
1072 //--------------------------------------
1073 
1074 const Valeur Isol_hor:: boundary_beta_y(double om)const {
1075 
1076  // Les alignemenents pour le signe des CL.
1077  double orientation = mp.get_rot_phi() ;
1078  assert ((orientation == 0) || (orientation == M_PI)) ;
1079  int aligne = (orientation == 0) ? 1 : -1 ;
1080 
1081 
1082  int nnp = mp.get_mg()->get_np(1) ;
1083  int nnt = mp.get_mg()->get_nt(1) ;
1084 
1085  Vector tmp_vect = nn() * gam().radial_vect() ;
1086  tmp_vect.change_triad(mp.get_bvect_cart() ) ;
1087 
1088  //Isol_hor boundary conditions
1089 
1090  Valeur lim_y (mp.get_mg()->get_angu()) ;
1091 
1092  lim_y = 1 ; // Juste pour affecter dans espace des configs ;
1093 
1094  Mtbl xa_mtbl (mp.get_mg()) ;
1095  xa_mtbl.set_etat_qcq() ;
1096  xa_mtbl = mp.xa ;
1097 
1098  Scalar cosp (mp) ;
1099  cosp = mp.cosp ;
1100  Scalar cost (mp) ;
1101  cost = mp.cost ;
1102  Scalar sinp (mp) ;
1103  sinp = mp.sinp ;
1104  Scalar sint (mp) ;
1105  sint = mp.sint ;
1106 
1107  Scalar dep_phi (mp) ;
1108  dep_phi = 0.0*sint*cosp ;
1109 
1110  for (int k=0 ; k<nnp ; k++)
1111  for (int j=0 ; j<nnt ; j++)
1112  lim_y.set(0, k, j, 0) = - aligne * om * xa_mtbl(1, k, j, 0) *(1 +
1113  dep_phi.val_grid_point(1, k, j, 0))
1114  + tmp_vect(2).val_grid_point(1, k, j, 0) ;
1115 
1116  lim_y.set_base(*(mp.get_mg()->std_base_vect_cart()[1])) ;
1117 
1118  return lim_y ;
1119 }
1120 
1121 // Component z of boundary value of beta
1122 //--------------------------------------
1123 
1125 
1126  int nnp = mp.get_mg()->get_np(1) ;
1127  int nnt = mp.get_mg()->get_nt(1) ;
1128 
1129  Vector tmp_vect = nn() * gam().radial_vect() ;
1130  tmp_vect.change_triad(mp.get_bvect_cart() ) ;
1131 
1132  //Isol_hor boundary conditions
1133 
1134  Valeur lim_z (mp.get_mg()->get_angu()) ;
1135 
1136  lim_z = 1 ; // Juste pour affecter dans espace des configs ;
1137 
1138  for (int k=0 ; k<nnp ; k++)
1139  for (int j=0 ; j<nnt ; j++)
1140  lim_z.set(0, k, j, 0) = tmp_vect(3).val_grid_point(1, k, j, 0) ;
1141 
1142  lim_z.set_base(*(mp.get_mg()->std_base_vect_cart()[2])) ;
1143 
1144  return lim_z ;
1145 }
1146 
1148 
1149  Valeur lim_x (mp.get_mg()->get_angu()) ;
1150  lim_x = boost_x ; // Juste pour affecter dans espace des configs ;
1151 
1152  lim_x.set_base(*(mp.get_mg()->std_base_vect_cart()[0])) ;
1153 
1154  return lim_x ;
1155 }
1156 
1157 
1159 
1160  Valeur lim_z (mp.get_mg()->get_angu()) ;
1161  lim_z = boost_z ; // Juste pour affecter dans espace des configs ;
1162 
1163  lim_z.set_base(*(mp.get_mg()->std_base_vect_cart()[2])) ;
1164 
1165  return lim_z ;
1166 }
1167 
1168 
1169 // Neumann boundary condition for b_tilde
1170 //---------------------------------------
1171 
1173 
1174  // Introduce 2-trace gamma tilde dot
1175 
1176  Vector s_tilde = met_gamt.radial_vect() ;
1177 
1178  Scalar hh_tilde = contract(s_tilde.derive_cov(met_gamt), 0, 1) ;
1179 
1180  //des_profile(hh_tilde, 1.00001, 10, M_PI/2., 0., "H_tilde") ;
1181 
1182  Scalar tmp (mp) ;
1183 
1184  tmp = + b_tilde() * hh_tilde - 2 * ( s_tilde(2) * b_tilde().derive_cov(ff)(2)
1185  + s_tilde(3) * b_tilde().derive_cov(ff)(3) ) ;
1186 
1187  Scalar constant (mp) ;
1188  constant = 0. ;
1189  constant.std_spectral_base() ;
1190  constant.inc_dzpuis(2) ;
1191 
1192  tmp += constant ;
1193  tmp = tmp / (2 * s_tilde(1) ) ;
1194 
1195  // in this case you don't have to substract any value
1196 
1197  Valeur b_tilde_bound (mp.get_mg()->get_angu() ) ;
1198 
1199  int nnp = mp.get_mg()->get_np(1) ;
1200  int nnt = mp.get_mg()->get_nt(1) ;
1201 
1202  b_tilde_bound = 1 ;
1203 
1204  for (int k=0 ; k<nnp ; k++)
1205  for (int j=0 ; j<nnt ; j++)
1206  b_tilde_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
1207 
1208  b_tilde_bound.std_base_scal() ;
1209 
1210  return b_tilde_bound ;
1211 
1212 }
1213 
1214 
1216 
1217  Vector s_tilde = met_gamt.radial_vect() ;
1218 
1219  Scalar hh_tilde = contract(s_tilde.derive_cov(met_gamt), 0, 1) ;
1220 
1221  Scalar tmp = 2 * contract (s_tilde, 0, b_tilde().derive_cov(ff) , 0) ;
1222 
1223  Scalar constant (mp) ;
1224  constant = -1. ;
1225  constant.std_spectral_base() ;
1226  constant.inc_dzpuis(2) ;
1227 
1228  tmp -= constant ;
1229 
1230  tmp = tmp / hh_tilde ;
1231 
1232  // des_profile(tmp, 1.00001, 10, M_PI/2., 0., "tmp") ;
1233 
1234  // We have substracted 1, since we solve for zero condition at infinity
1235  //and then we add 1 to the solution
1236 
1237  Valeur b_tilde_bound (mp.get_mg()->get_angu() ) ;
1238 
1239  int nnp = mp.get_mg()->get_np(1) ;
1240  int nnt = mp.get_mg()->get_nt(1) ;
1241 
1242  b_tilde_bound = 1 ;
1243 
1244  for (int k=0 ; k<nnp ; k++)
1245  for (int j=0 ; j<nnt ; j++)
1246  b_tilde_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
1247 
1248  b_tilde_bound.std_base_scal() ;
1249 
1250  return b_tilde_bound ;
1251 
1252 }
1253 
1254 const Vector Isol_hor::vv_bound_cart(double om) const{
1255 
1256  // Preliminaries
1257  //--------------
1258 
1259  // HH_tilde
1260  Vector s_tilde = met_gamt.radial_vect() ;
1261 
1262  Scalar hh_tilde = contract(s_tilde.derive_cov(met_gamt), 0, 1) ;
1263  hh_tilde.dec_dzpuis(2) ;
1264 
1265  // Tangential part of the shift
1266  Vector tmp_vect = b_tilde() * s_tilde ;
1267 
1268  Vector VV = b_tilde() * s_tilde - beta() ;
1269 
1270  // "Acceleration" term V^a \tilde{D}_a ln M
1271  Scalar accel = 2 * contract( VV, 0, contract( s_tilde, 0, s_tilde.down(0,
1272  met_gamt).derive_cov(met_gamt), 1), 0) ;
1273 
1274 
1275  // Divergence of V^a
1276  Sym_tensor qq_spher = met_gamt.con() - s_tilde * s_tilde ;
1277  Scalar div_VV = contract( qq_spher.down(0, met_gamt), 0, 1, VV.derive_cov(met_gamt), 0, 1) ;
1278 
1279 
1280  Scalar cosp (mp) ;
1281  cosp = mp.cosp ;
1282  Scalar cost (mp) ;
1283  cost = mp.cost ;
1284  Scalar sinp (mp) ;
1285  sinp = mp.sinp ;
1286  Scalar sint (mp) ;
1287  sint = mp.sint ;
1288 
1289  Scalar dep_phi (mp) ;
1290  dep_phi = 0.0*sint*cosp ;
1291 
1292  // Les alignemenents pour le signe des CL.
1293  double orientation = mp.get_rot_phi() ;
1294  assert ((orientation == 0) || (orientation == M_PI)) ;
1295  int aligne = (orientation == 0) ? 1 : -1 ;
1296 
1297  Vector angular (mp, CON, mp.get_bvect_cart()) ;
1298  Scalar yya (mp) ;
1299  yya = mp.ya ;
1300  Scalar xxa (mp) ;
1301  xxa = mp.xa ;
1302 
1303  angular.set(1) = - aligne * om * yya * (1 + dep_phi) ;
1304  angular.set(2) = aligne * om * xxa * (1 + dep_phi) ;
1305  angular.set(3).annule_hard() ;
1306 
1307 
1308  angular.set(1).set_spectral_va()
1309  .set_base(*(mp.get_mg()->std_base_vect_cart()[0])) ;
1310  angular.set(2).set_spectral_va()
1311  .set_base(*(mp.get_mg()->std_base_vect_cart()[1])) ;
1312  angular.set(3).set_spectral_va()
1313  .set_base(*(mp.get_mg()->std_base_vect_cart()[2])) ;
1314 
1315  angular.change_triad(mp.get_bvect_spher()) ;
1316 
1317 /*
1318  Scalar ang_vel (mp) ;
1319  ang_vel = om * (1 + dep_phi) ;
1320  ang_vel.std_spectral_base() ;
1321  ang_vel.mult_rsint() ;
1322 */
1323 
1324  Scalar kss (mp) ;
1325  kss = - 0.2 ;
1326  // kss.std_spectral_base() ;
1327  // kss.inc_dzpuis(2) ;
1328 
1329 
1330  //kss from from L_lN=0 condition
1331  //------------------------------
1332  kss = - 0.15 / nn() ;
1333  kss.inc_dzpuis(2) ;
1334  kss += contract(gam().radial_vect(), 0, nn().derive_cov(ff), 0) / nn() ;
1335 
1336 
1337  cout << "kappa_hor = " << kappa_hor() << endl ;
1338 
1339  /*
1340  // Apparent horizon condition
1341  //---------------------------
1342  kss = trk() ;
1343  kss -= (4* contract(psi().derive_cov(met_gamt), 0, met_gamt.radial_vect(),
1344  0)/psi() +
1345  contract(met_gamt.radial_vect().derive_cov(met_gamt), 0, 1)) /
1346  (psi() * psi()) ;
1347  */
1348 
1349 
1350  // beta^r component
1351  //-----------------
1352  double rho = 5. ; // rho>1 ; rho=1 "pure Dirichlet" version
1353  Scalar beta_r_corr = (rho - 1) * b_tilde() * hh_tilde;
1354  beta_r_corr.inc_dzpuis(2) ;
1355  Scalar nn_KK = nn() * trk() ;
1356  nn_KK.inc_dzpuis(2) ;
1357 
1358  beta_r_corr.set_dzpuis(2) ;
1359  nn_KK.set_dzpuis(2) ;
1360  accel.set_dzpuis(2) ;
1361  div_VV.set_dzpuis(2) ;
1362 
1363  Scalar b_tilde_new (mp) ;
1364  b_tilde_new = 2 * contract(s_tilde, 0, b_tilde().derive_cov(ff), 0)
1365  + beta_r_corr
1366  - 3 * nn() * kss + nn_KK + accel + div_VV ;
1367 
1368  b_tilde_new = b_tilde_new / (hh_tilde * rho) ;
1369 
1370  b_tilde_new.dec_dzpuis(2) ;
1371 
1372  tmp_vect.set(1) = b_tilde_new * s_tilde(1) + angular(1) ;
1373  tmp_vect.set(2) = b_tilde_new * s_tilde(2) + angular(2) ;
1374  tmp_vect.set(3) = b_tilde_new * s_tilde(3) + angular(3) ;
1375 
1376  tmp_vect.change_triad(mp.get_bvect_cart() ) ;
1377 
1378  return tmp_vect ;
1379 }
1380 
1381 
1382 const Vector Isol_hor::vv_bound_cart_bin(double, int ) const{
1383  /*
1384  // Les alignemenents pour le signe des CL.
1385  double orientation = mp.get_rot_phi() ;
1386  assert ((orientation == 0) || (orientation == M_PI)) ;
1387  int aligne = (orientation == 0) ? 1 : -1 ;
1388 
1389  Vector angular (mp, CON, mp.get_bvect_cart()) ;
1390  Scalar yya (mp) ;
1391  yya = mp.ya ;
1392  Scalar xxa (mp) ;
1393  xxa = mp.xa ;
1394 
1395  angular.set(1) = aligne * om * yya ;
1396  angular.set(2) = - aligne * om * xxa ;
1397  angular.set(3).annule_hard() ;
1398 
1399  angular.set(1).set_spectral_va()
1400  .set_base(*(mp.get_mg()->std_base_vect_cart()[0])) ;
1401  angular.set(2).set_spectral_va()
1402  .set_base(*(mp.get_mg()->std_base_vect_cart()[1])) ;
1403  angular.set(3).set_spectral_va()
1404  .set_base(*(mp.get_mg()->std_base_vect_cart()[2])) ;
1405 
1406  angular.change_triad(mp.get_bvect_spher()) ;
1407 
1408  // HH_tilde
1409  Vector s_tilde = met_gamt.radial_vect() ;
1410 
1411  Scalar hh_tilde = contract(s_tilde.derive_cov(met_gamt), 0, 1) ;
1412  hh_tilde.dec_dzpuis(2) ;
1413 
1414  Scalar btilde = b_tilde() - contract(angular, 0, s_tilde.up_down(met_gamt), 0) ;
1415 
1416  // Tangential part of the shift
1417  Vector tmp_vect = btilde * s_tilde ;
1418 
1419  // Value of kss
1420  // --------------
1421 
1422  Scalar kss (mp) ;
1423  kss = - 0.26 ;
1424  kss.std_spectral_base() ;
1425  kss.inc_dzpuis(2) ;
1426 
1427  // Apparent horizon boundary condition
1428  // -----------------------------------
1429 
1430 
1431  // kss frome a fich
1432 
1433  // Construction of the binary
1434  // --------------------------
1435 
1436  char* name_fich = "/home/francois/resu/bin_hor/Test/b11_9x9x8/bin.dat" ;
1437 
1438  FILE* fich_s = fopen(name_fich, "r") ;
1439  Mg3d grid_s (fich_s) ;
1440  Map_af map_un_s (grid_s, fich_s) ;
1441  Map_af map_deux_s (grid_s, fich_s) ;
1442  Bin_hor bin (map_un_s, map_deux_s, fich_s) ;
1443  fclose(fich_s) ;
1444 
1445  // Inititialisation of fields :
1446  // ----------------------------
1447 
1448  bin.set(1).n_comp_import (bin(2)) ;
1449  bin.set(1).psi_comp_import (bin(2)) ;
1450  bin.set(2).n_comp_import (bin(1)) ;
1451  bin.set(2).psi_comp_import (bin(1)) ;
1452  bin.decouple() ;
1453  bin.extrinsic_curvature() ;
1454 
1455  kss = contract(bin(jj).get_k_dd(), 0, 1, bin(jj).get_gam().radial_vect()*
1456  bin(jj).get_gam().radial_vect(), 0, 1) ;
1457 
1458 
1459  cout << "kappa_hor = " << kappa_hor() << endl ;
1460 
1461  // beta^r component
1462  //-----------------
1463  double rho = 5. ; // rho=0 "pure Dirichlet" version
1464  Scalar beta_r_corr = rho * btilde * hh_tilde;
1465  beta_r_corr.inc_dzpuis(2) ;
1466  Scalar nn_KK = nn() * trk() ;
1467  nn_KK.inc_dzpuis(2) ;
1468 
1469  beta_r_corr.set_dzpuis(2) ;
1470  nn_KK.set_dzpuis(2) ;
1471 
1472  Scalar b_tilde_new (mp) ;
1473  b_tilde_new = 2 * contract(s_tilde, 0, btilde.derive_cov(ff), 0)
1474  + beta_r_corr - 3 * nn() * kss + nn_KK ;
1475 
1476  b_tilde_new = b_tilde_new / (hh_tilde * (rho+1)) ;
1477 
1478  b_tilde_new.dec_dzpuis(2) ;
1479 
1480  tmp_vect.set(1) = b_tilde_new * s_tilde(1) + angular(1) ;
1481  tmp_vect.set(2) = b_tilde_new * s_tilde(2) + angular(2) ;
1482  tmp_vect.set(3) = b_tilde_new * s_tilde(3) + angular(3) ;
1483 
1484  tmp_vect.change_triad(mp.get_bvect_cart() ) ;
1485 
1486  return tmp_vect ;
1487  */
1488  Vector pipo(mp, CON, mp.get_bvect_cart()) ;
1489  return pipo ;
1490 }
1491 
1492 
1493 // Component x of boundary value of V^i
1494 //-------------------------------------
1495 
1496 const Valeur Isol_hor:: boundary_vv_x(double om)const {
1497 
1498  int nnp = mp.get_mg()->get_np(1) ;
1499  int nnt = mp.get_mg()->get_nt(1) ;
1500 
1501  //Isol_hor boundary conditions
1502 
1503  Valeur lim_x (mp.get_mg()->get_angu()) ;
1504 
1505  lim_x = 1 ; // Juste pour affecter dans espace des configs ;
1506 
1507  Scalar vv_x = vv_bound_cart(om)(1) ;
1508 
1509  for (int k=0 ; k<nnp ; k++)
1510  for (int j=0 ; j<nnt ; j++)
1511  lim_x.set(0, k, j, 0) = vv_x.val_grid_point(1, k, j, 0) ;
1512 
1513  lim_x.set_base(*(mp.get_mg()->std_base_vect_cart()[0])) ;
1514 
1515  return lim_x ;
1516 
1517 
1518 }
1519 
1520 
1521 // Component y of boundary value of V^i
1522 //--------------------------------------
1523 
1524 const Valeur Isol_hor:: boundary_vv_y(double om)const {
1525 
1526  int nnp = mp.get_mg()->get_np(1) ;
1527  int nnt = mp.get_mg()->get_nt(1) ;
1528 
1529  // Isol_hor boundary conditions
1530 
1531  Valeur lim_y (mp.get_mg()->get_angu()) ;
1532 
1533  lim_y = 1 ; // Juste pour affecter dans espace des configs ;
1534 
1535  Scalar vv_y = vv_bound_cart(om)(2) ;
1536 
1537  for (int k=0 ; k<nnp ; k++)
1538  for (int j=0 ; j<nnt ; j++)
1539  lim_y.set(0, k, j, 0) = vv_y.val_grid_point(1, k, j, 0) ;
1540 
1541  lim_y.set_base(*(mp.get_mg()->std_base_vect_cart()[1])) ;
1542 
1543  return lim_y ;
1544 }
1545 
1546 
1547 // Component z of boundary value of V^i
1548 //-------------------------------------
1549 
1550 const Valeur Isol_hor:: boundary_vv_z(double om)const {
1551 
1552  int nnp = mp.get_mg()->get_np(1) ;
1553  int nnt = mp.get_mg()->get_nt(1) ;
1554 
1555  // Isol_hor boundary conditions
1556 
1557  Valeur lim_z (mp.get_mg()->get_angu()) ;
1558 
1559  lim_z = 1 ; // Juste pour affecter dans espace des configs ;
1560 
1561  Scalar vv_z = vv_bound_cart(om)(3) ;
1562 
1563  for (int k=0 ; k<nnp ; k++)
1564  for (int j=0 ; j<nnt ; j++)
1565  lim_z.set(0, k, j, 0) = vv_z.val_grid_point(1, k, j, 0) ;
1566 
1567  lim_z.set_base(*(mp.get_mg()->std_base_vect_cart()[2])) ;
1568 
1569  return lim_z ;
1570 }
1571 
1572 // Component x of boundary value of V^i
1573 //-------------------------------------
1574 
1575 const Valeur Isol_hor:: boundary_vv_x_bin(double om, int jj)const {
1576 
1577  int nnp = mp.get_mg()->get_np(1) ;
1578  int nnt = mp.get_mg()->get_nt(1) ;
1579 
1580  //Isol_hor boundary conditions
1581 
1582  Valeur lim_x (mp.get_mg()->get_angu()) ;
1583 
1584  lim_x = 1 ; // Juste pour affecter dans espace des configs ;
1585 
1586  Scalar vv_x = vv_bound_cart_bin(om, jj)(1) ;
1587 
1588  for (int k=0 ; k<nnp ; k++)
1589  for (int j=0 ; j<nnt ; j++)
1590  lim_x.set(0, k, j, 0) = vv_x.val_grid_point(1, k, j, 0) ;
1591 
1592  lim_x.set_base(*(mp.get_mg()->std_base_vect_cart()[0])) ;
1593 
1594  return lim_x ;
1595 
1596 
1597 }
1598 
1599 // Component y of boundary value of V^i
1600 //--------------------------------------
1601 
1602 const Valeur Isol_hor:: boundary_vv_y_bin(double om, int jj)const {
1603 
1604  int nnp = mp.get_mg()->get_np(1) ;
1605  int nnt = mp.get_mg()->get_nt(1) ;
1606 
1607  // Isol_hor boundary conditions
1608 
1609  Valeur lim_y (mp.get_mg()->get_angu()) ;
1610 
1611  lim_y = 1 ; // Juste pour affecter dans espace des configs ;
1612 
1613  Scalar vv_y = vv_bound_cart_bin(om, jj)(2) ;
1614 
1615  for (int k=0 ; k<nnp ; k++)
1616  for (int j=0 ; j<nnt ; j++)
1617  lim_y.set(0, k, j, 0) = vv_y.val_grid_point(1, k, j, 0) ;
1618 
1619  lim_y.set_base(*(mp.get_mg()->std_base_vect_cart()[1])) ;
1620 
1621  return lim_y ;
1622 }
1623 
1624 
1625 // Component z of boundary value of V^i
1626 //-------------------------------------
1627 
1628 const Valeur Isol_hor:: boundary_vv_z_bin(double om, int jj)const {
1629 
1630  int nnp = mp.get_mg()->get_np(1) ;
1631  int nnt = mp.get_mg()->get_nt(1) ;
1632 
1633  // Isol_hor boundary conditions
1634 
1635  Valeur lim_z (mp.get_mg()->get_angu()) ;
1636 
1637  lim_z = 1 ; // Juste pour affecter dans espace des configs ;
1638 
1639  Scalar vv_z = vv_bound_cart_bin(om, jj)(3) ;
1640 
1641  for (int k=0 ; k<nnp ; k++)
1642  for (int j=0 ; j<nnt ; j++)
1643  lim_z.set(0, k, j, 0) = vv_z.val_grid_point(1, k, j, 0) ;
1644 
1645  lim_z.set_base(*(mp.get_mg()->std_base_vect_cart()[2])) ;
1646 
1647  return lim_z ;
1648 }
1649 }
Coord xa
Absolute x coordinate.
Definition: map.h:748
virtual const Vector & beta() const
shift vector at the current time step (jtime )
const Vector radial_vect_hor() const
Vector radial normal.
Definition: phys_param.C:98
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
const Valeur boundary_psi_Dir_evol() const
Dirichlet boundary condition for (evolution)
Definition: bound_hor.C:177
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
const Valeur boundary_psi_Dir_spat() const
Dirichlet boundary condition for (spatial)
Definition: bound_hor.C:234
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:293
const Valeur boundary_vv_z(double om) const
Component z of boundary value of .
Definition: bound_hor.C:1550
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:273
virtual const Sym_tensor & gam_uu() const
Induced metric (contravariant components ) at the current time step (jtime )
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
const Valeur boundary_nn_Dir(double aa) const
Dirichlet boundary condition .
Definition: bound_hor.C:650
const Scalar & lapang() const
Returns the angular Laplacian of *this , where .
Definition: scalar_deriv.C:461
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
const Valeur boundary_vv_z_bin(double om, int hole=0) const
Component z of boundary value of .
Definition: bound_hor.C:1628
const Valeur boundary_nn_Dir_lapl(int mer=1) const
Dirichlet boundary condition for N fixing the divergence of the connection form . ...
Definition: bound_hor.C:696
const Valeur boundary_nn_Dir_eff(double aa) const
Dirichlet boundary condition for N (effectif) .
Definition: bound_hor.C:589
Multi-domain array.
Definition: mtbl.h:118
const Valeur boundary_beta_x(double om) const
Component x of boundary value of .
Definition: bound_hor.C:1024
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
double radius
Radius of the horizon in LORENE&#39;s units.
Definition: isol_hor.h:266
int get_base_t(int l) const
Returns the expansion basis for functions in the domain of index l (e.g.
Definition: base_val.h:414
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
const Valeur boundary_nn_Dir_kk() const
Dirichlet boundary condition for N using the extrinsic curvature.
Definition: bound_hor.C:374
const Vector vv_bound_cart(double om) const
Vector for boundary conditions in cartesian.
Definition: bound_hor.C:1254
const Valeur boundary_vv_y(double om) const
Component y of boundary value of .
Definition: bound_hor.C:1524
double kappa_hor() const
Surface gravity.
Definition: phys_param.C:221
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
const Valeur boundary_psi_app_hor() const
Neumann boundary condition for (spatial)
Definition: bound_hor.C:300
virtual const Vector & radial_vect() const
Returns the radial vector normal to a spherical slicing and pointing toward spatial infinity...
Definition: metric.C:365
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
void set(const Map *mp, Mtbl *(*construct)(const Map *))
Semi-constructor from a mapping and a method.
Definition: coord.C:137
const Valeur boundary_psi_Dir() const
Dirichlet boundary condition for (spatial)
Definition: bound_hor.C:327
Tensor field of valence 1.
Definition: vector.h:188
Scalar poisson_angu(double lambda=0) const
Solves the (generalized) angular Poisson equation with *this as source.
Definition: scalar_pde.C:203
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:814
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition: scalar.h:643
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:386
int get_base_r(int l) const
Returns the expansion basis for r ( ) functions in the domain of index l (e.g.
Definition: base_val.h:403
Coord sint
Definition: map.h:739
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric_flat.C:156
Map_af & mp
Affine mapping.
Definition: isol_hor.h:260
const Valeur boundary_nn_Neu_eff(double aa) const
Neumann boundary condition on nn (effectif) .
Definition: bound_hor.C:625
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:484
virtual const Scalar & trk() const
Trace K of the extrinsic curvature at the current time step (jtime )
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
double boost_x
Boost velocity in x-direction.
Definition: isol_hor.h:272
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition: mg3d.C:604
Base_val base
Bases on which the spectral expansion is performed.
Definition: valeur.h:315
const Scalar b_tilde() const
Radial component of the shift with respect to the conformal metric.
Definition: phys_param.C:139
Metric met_gamt
3 metric tilde
Definition: isol_hor.h:326
const Valeur boundary_beta_phi(double om) const
Component phi of boundary value of .
Definition: bound_hor.C:981
Coord sinp
Definition: map.h:741
const Metric & gam() const
Induced metric at the current time step (jtime )
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:387
int get_base_p(int l) const
Returns the expansion basis for functions in the domain of index l (e.g.
Definition: base_val.h:425
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
const Valeur boundary_vv_x(double om) const
Component x of boundary value of .
Definition: bound_hor.C:1496
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
const Valeur boundary_beta_z() const
Component z of boundary value of .
Definition: bound_hor.C:1124
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
const Valeur boundary_beta_r() const
Component r of boundary value of .
Definition: bound_hor.C:900
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl.C:302
virtual const Sym_tensor & k_dd() const
Extrinsic curvature tensor (covariant components ) at the current time step (jtime ) ...
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:283
Coord ya
Absolute y coordinate.
Definition: map.h:749
const Valeur boundary_nn_Neu_kk(int nn=1) const
Neumann boundary condition for N using the extrinsic curvature.
Definition: bound_hor.C:423
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
Bases of the spectral expansions.
Definition: base_val.h:325
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric_flat.C:137
const Valeur boundary_vv_x_bin(double om, int hole=0) const
Component x of boundary value of .
Definition: bound_hor.C:1575
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
Definition: time_slice.h:510
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
const Valeur boundary_b_tilde_Dir() const
Dirichlet boundary condition for b_tilde.
Definition: bound_hor.C:1215
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:663
Coord cosp
Definition: map.h:742
const Valeur boundary_psi_Neu_spat() const
Neumann boundary condition for (spatial)
Definition: bound_hor.C:270
const Valeur boundary_beta_y(double om) const
Component y of boundary value of .
Definition: bound_hor.C:1074
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
const Valeur boundary_beta_theta() const
Component theta of boundary value of .
Definition: bound_hor.C:941
const Vector tradial_vect_hor() const
Vector radial normal tilde.
Definition: phys_param.C:119
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
const Valeur beta_boost_z() const
Boundary value for a boost in z-direction.
Definition: bound_hor.C:1158
const Valeur boundary_vv_y_bin(double om, int hole=0) const
Component y of boundary value of .
Definition: bound_hor.C:1602
const Vector vv_bound_cart_bin(double om, int hole=0) const
Vector for boundary conditions in cartesian for binary systems.
Definition: bound_hor.C:1382
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...
const Valeur boundary_b_tilde_Neu() const
Neumann boundary condition for b_tilde.
Definition: bound_hor.C:1172
const Valeur beta_boost_x() const
Boundary value for a boost in x-direction.
Definition: bound_hor.C:1147
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:226
double boost_z
Boost velocity in z-direction.
Definition: isol_hor.h:275
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:607
const Valeur boundary_nn_Neu_Cook() const
Neumann boundary condition for N using Cook&#39;s boundary condition.
Definition: bound_hor.C:492
Coord r
r coordinate centered on the grid
Definition: map.h:736
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:506
const Valeur boundary_psi_Neu_evol() const
Neumann boundary condition for (evolution)
Definition: bound_hor.C:206
const Base_val & get_base() const
Return the bases for spectral expansions (member base )
Definition: valeur.h:490
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition: valeur.h:373
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )
Tensor down(int ind, const Metric &gam) const
Computes a new tensor by lowering an index of *this.
Coord cost
Definition: map.h:740