LORENE
init_data.C
1 /*
2  * Method of class Hor_isol to compute valid initial data for standard boundary
3  * conditions
4  *
5  * (see file isol_hor.h for documentation).
6  *
7  */
8 
9 /*
10  * Copyright (c) 2004 Jose Luis Jaramillo
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 char init_data_C[] = "$Header: /cvsroot/Lorene/C++/Source/Isol_hor/init_data.C,v 1.30 2014/10/13 08:53:01 j_novak Exp $" ;
30 
31 /*
32  * $Id: init_data.C,v 1.30 2014/10/13 08:53:01 j_novak Exp $
33  * $Log: init_data.C,v $
34  * Revision 1.30 2014/10/13 08:53:01 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.29 2014/10/06 15:13:10 j_novak
38  * Modified #include directives to use c++ syntax.
39  *
40  * Revision 1.28 2008/08/19 06:42:00 j_novak
41  * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
42  * cast-type operations, and constant strings that must be defined as const char*
43  *
44  * Revision 1.27 2006/02/22 17:02:04 f_limousin
45  * Removal of warnings
46  *
47  * Revision 1.26 2006/02/22 16:29:55 jl_jaramillo
48  * corrections on the relaxation and boundary conditions
49  *
50  * Revision 1.24 2006/01/18 09:04:27 f_limousin
51  * Minor modifs (warnings and errors at the compilation with gcc-3.4)
52  *
53  * Revision 1.23 2006/01/16 17:13:40 jl_jaramillo
54  * function for solving the spherical case
55  *
56  * Revision 1.22 2005/11/02 16:09:44 jl_jaramillo
57  * changes in boundary_nn_Dir_lapl
58  *
59  * Revision 1.21 2005/10/24 16:44:40 jl_jaramillo
60  * Cook boundary condition ans minot bound of kss
61  *
62  * Revision 1.20 2005/10/21 16:20:55 jl_jaramillo
63  * Version for the paper JaramL05
64  *
65  * Revision 1.19 2005/07/08 13:15:23 f_limousin
66  * Improvements of boundary_vv_cart(), boundary_nn_lapl().
67  * Add a fonction to compute the departure of axisymmetry.
68  *
69  * Revision 1.18 2005/06/09 08:05:32 f_limousin
70  * Implement a new function sol_elliptic_boundary() and
71  * Vector::poisson_boundary(...) which solve the vectorial poisson
72  * equation (method 6) with an inner boundary condition.
73  *
74  * Revision 1.17 2005/05/12 14:48:07 f_limousin
75  * New boundary condition for the lapse : boundary_nn_lapl().
76  *
77  * Revision 1.16 2005/04/08 12:16:52 f_limousin
78  * Function set_psi(). And dependance in phi.
79  *
80  * Revision 1.15 2005/04/03 19:48:22 f_limousin
81  * Implementation of set_psi(psi_in). And minor changes to avoid warnings.
82  *
83  * Revision 1.14 2005/04/02 15:49:21 f_limousin
84  * New choice (Lichnerowicz) for aaquad. New member data nz.
85  *
86  * Revision 1.13 2005/03/31 09:45:31 f_limousin
87  * New functions compute_ww(...) and aa_kerr_ww().
88  *
89  * Revision 1.12 2005/03/24 16:50:28 f_limousin
90  * Add parameters solve_shift and solve_psi in par_isol.d and in function
91  * init_dat(...). Implement Isolhor::kerr_perturb().
92  *
93  * Revision 1.11 2005/03/22 13:25:36 f_limousin
94  * Small changes. The angular velocity and A^{ij} are computed
95  * with a differnet sign.
96  *
97  * Revision 1.10 2005/03/09 10:18:08 f_limousin
98  * Save K_{ij}s^is^j in a file. Add solve_lapse in a file
99  *
100  * Revision 1.9 2005/03/06 16:56:13 f_limousin
101  * The computation of A^{ij} is no more necessary here thanks to the new
102  * function Isol_hor::aa().
103  *
104  * Revision 1.8 2005/03/04 17:04:57 jl_jaramillo
105  * Addition of boost to the shift after solving the shift equation
106  *
107  * Revision 1.7 2005/03/03 10:03:55 f_limousin
108  * The boundary conditions for the lapse, psi and shift are now
109  * parameters (in file par_hor.d).
110  *
111  * Revision 1.6 2004/12/22 18:15:30 f_limousin
112  * Many different changes.
113  *
114  * Revision 1.5 2004/11/08 14:51:21 f_limousin
115  * A regularisation for the computation of A^{ij } is done in the
116  * case lapse equal to zero on the horizon.
117  *
118  * Revision 1.1 2004/10/29 12:54:53 jl_jaramillo
119  * First version
120  *
121  * Revision 1.4 2004/10/01 16:47:51 f_limousin
122  * Case \alpha=0 included
123  *
124  * Revision 1.3 2004/09/28 16:10:05 f_limousin
125  * Many improvements. Now the resolution for the shift is working !
126  *
127  * Revision 1.1 2004/09/09 16:41:50 f_limousin
128  * First version
129  *
130  *
131  * $Header: /cvsroot/Lorene/C++/Source/Isol_hor/init_data.C,v 1.30 2014/10/13 08:53:01 j_novak Exp $
132  *
133  */
134 
135 // C++ headers
136 #include "headcpp.h"
137 
138 // C headers
139 #include <cstdlib>
140 #include <cassert>
141 
142 // Lorene headers
143 #include "isol_hor.h"
144 #include "metric.h"
145 #include "unites.h"
146 #include "graphique.h"
147 #include "cmp.h"
148 #include "tenseur.h"
149 #include "utilitaires.h"
150 #include "param.h"
151 
152 namespace Lorene {
153 void Isol_hor::init_data(int bound_nn, double lim_nn, int bound_psi,
154  int bound_beta, int solve_lapse, int solve_psi,
155  int solve_shift, double precis,
156  double relax_nn, double relax_psi, double relax_beta, int niter) {
157 
158  using namespace Unites ;
159 
160  // Initialisations
161  // ---------------
162  double ttime = the_time[jtime] ;
163 
164  ofstream conv("resconv.d") ;
165  ofstream kss("kss.d") ;
166  conv << " # diff_nn diff_psi diff_beta " << endl ;
167 
168  // Iteration
169  // ---------
170 // double relax_nn_fin = relax_nn ;
171 // double relax_psi_fin = relax_psi ;
172 // double relax_beta_fin = relax_beta ;
173 
174 
175  for (int mer=0; mer<niter; mer++) {
176 
177  //=========================================================
178  // Boundary conditions and resolution of elliptic equations
179  //=========================================================
180 
181  // Resolution of the Poisson equation for the lapse
182  // ------------------------------------------------
183 
184  double relax_init = 0.05 ;
185  double relax_speed = 0.005 ;
186 
187  double corr = 1 - (1 - relax_init) * exp (- relax_speed * mer) ;
188 
189  // relax_nn = relax_nn_fin - ( relax_nn_fin - relax_init ) * exp (- relax_speed * mer) ;
190  // relax_psi = relax_psi_fin - ( relax_psi_fin - relax_init ) * exp (- relax_speed * mer) ;
191  // relax_beta = relax_beta_fin - ( relax_beta_fin - relax_init ) * exp (- relax_speed * mer) ;
192 
193  cout << "nn = " << mer << " corr = " << corr << endl ;
194 
195 
196 
197  cout << " relax_nn = " << relax_nn << endl ;
198  cout << " relax_psi = " << relax_psi << endl ;
199  cout << " relax_beta = " << relax_beta << endl ;
200 
201 
202  Scalar sou_nn (source_nn()) ;
203  Scalar nn_jp1 (mp) ;
204  if (solve_lapse == 1) {
205  Valeur nn_bound (mp.get_mg()-> get_angu()) ;
206 
207  switch (bound_nn) {
208 
209  case 0 : {
210  nn_bound = boundary_nn_Dir(lim_nn) ;
211  nn_jp1 = sou_nn.poisson_dirichlet(nn_bound, 0) + 1. ;
212  break ;
213  }
214  case 1 : {
215  nn_bound = boundary_nn_Neu_eff(lim_nn) ;
216  nn_jp1 = sou_nn.poisson_neumann(nn_bound, 0) + 1. ;
217  break ;
218  }
219  case 2 : {
220  nn_bound = boundary_nn_Dir_eff(lim_nn) ;
221  nn_jp1 = sou_nn.poisson_dirichlet(nn_bound, 0) + 1. ;
222  break ;
223  }
224  case 3 : {
225  nn_bound = boundary_nn_Neu_kk(mer) ;
226  nn_jp1 = sou_nn.poisson_neumann(nn_bound, 0) + 1. ;
227  break ;
228  }
229  case 4 : {
230  nn_bound = boundary_nn_Dir_kk() ;
231  nn_jp1 = sou_nn.poisson_dirichlet(nn_bound, 0) + 1. ;
232  break ;
233  }
234  case 5 : {
235  nn_bound = boundary_nn_Dir_lapl(mer) ;
236  nn_jp1 = sou_nn.poisson_dirichlet(nn_bound, 0) + 1. ;
237  break ;
238  }
239  case 6 : {
240  nn_bound = boundary_nn_Neu_Cook() ;
241  nn_jp1 = sou_nn.poisson_neumann(nn_bound, 0) + 1. ;
242  break ;
243  }
244 
245 
246 
247 
248  default : {
249  cout <<"Unexpected type of boundary conditions for the lapse!"
250  << endl
251  << " bound_nn = " << bound_nn << endl ;
252  abort() ;
253  break ;
254  }
255 
256  } // End of switch
257 
258  // Test:
259  maxabs(nn_jp1.laplacian() - sou_nn,
260  "Absolute error in the resolution of the equation for N") ;
261 
262  // Relaxation (relax=1 -> new ; relax=0 -> old )
263  if (mer==0)
264  n_evol.update(nn_jp1, jtime, ttime) ;
265  else
266  nn_jp1 = relax_nn * nn_jp1 + (1 - relax_nn) * nn() ;
267  }
268 
269 
270  // Resolution of the Poisson equation for Psi
271  // ------------------------------------------
272 
273  Scalar sou_psi (source_psi()) ;
274  Scalar psi_jp1 (mp) ;
275  if (solve_psi == 1) {
276  Valeur psi_bound (mp.get_mg()-> get_angu()) ;
277 
278  switch (bound_psi) {
279 
280  case 0 : {
281  psi_bound = boundary_psi_app_hor() ;
282  psi_jp1 = sou_psi.poisson_neumann(psi_bound, 0) + 1. ;
283  break ;
284  }
285  case 1 : {
286  psi_bound = boundary_psi_Neu_spat() ;
287  psi_jp1 = sou_psi.poisson_neumann(psi_bound, 0) + 1. ;
288  break ;
289  }
290  case 2 : {
291  psi_bound = boundary_psi_Dir_spat() ;
292  psi_jp1 = sou_psi.poisson_dirichlet(psi_bound, 0) + 1. ;
293  break ;
294  }
295  case 3 : {
296  psi_bound = boundary_psi_Neu_evol() ;
297  psi_jp1 = sou_psi.poisson_neumann(psi_bound, 0) + 1. ;
298  break ;
299  }
300  case 4 : {
301  psi_bound = boundary_psi_Dir_evol() ;
302  psi_jp1 = sou_psi.poisson_dirichlet(psi_bound, 0) + 1. ;
303  break ;
304  }
305  case 5 : {
306  psi_bound = boundary_psi_Dir() ;
307  psi_jp1 = sou_psi.poisson_dirichlet(psi_bound, 0) + 1. ;
308  break ;
309  }
310  default : {
311  cout <<"Unexpected type of boundary conditions for psi!"
312  << endl
313  << " bound_psi = " << bound_psi << endl ;
314  abort() ;
315  break ;
316  }
317 
318  } // End of switch
319 
320  // Test:
321  maxabs(psi_jp1.laplacian() - sou_psi,
322  "Absolute error in the resolution of the equation for Psi") ;
323  // Relaxation (relax=1 -> new ; relax=0 -> old )
324  psi_jp1 = relax_psi * psi_jp1 + (1 - relax_psi) * psi() ;
325  }
326 
327  // Resolution of the vector Poisson equation for the shift
328  //---------------------------------------------------------
329 
330  // Source
331 
332  Vector beta_jp1(beta()) ;
333 
334  if (solve_shift == 1) {
335  Vector source_vector ( source_beta() ) ;
336  double lambda = 0; //1./3.;
337  Vector source_reg = - (1./3. - lambda) * beta().divergence(ff)
338  .derive_con(ff) ;
339  source_reg.inc_dzpuis() ;
340  source_vector = source_vector + source_reg ;
341 
342 
343  // CARTESIAN CASE
344  // #################################
345 
346  // Boundary values
347 
348  Valeur boundary_x (mp.get_mg()-> get_angu()) ;
349  Valeur boundary_y (mp.get_mg()-> get_angu()) ;
350  Valeur boundary_z (mp.get_mg()-> get_angu()) ;
351 
352  switch (bound_beta) {
353 
354  case 0 : {
355  boundary_x = boundary_beta_x(omega) ;
356  boundary_y = boundary_beta_y(omega) ;
357  boundary_z = boundary_beta_z() ;
358  break ;
359  }
360  case 1 : {
361  boundary_x = boundary_vv_x(omega) ;
362  boundary_y = boundary_vv_y(omega) ;
363  boundary_z = boundary_vv_z(omega) ;
364  break ;
365  }
366  default : {
367  cout <<"Unexpected type of boundary conditions for psi!"
368  << endl
369  << " bound_psi = " << bound_psi << endl ;
370  abort() ;
371  break ;
372  }
373  } // End of switch
374 
375  if (boost_x != 0.)
376  boundary_x -= beta_boost_x() ;
377  if (boost_z != 0.)
378  boundary_z -= beta_boost_z() ;
379 
380  // Resolution
381  //-----------
382 
383  double precision = 1e-8 ;
384  poisson_vect_boundary(lambda, source_vector, beta_jp1, boundary_x,
385  boundary_y, boundary_z, 0, precision, 20) ;
386 
387 
388 /*
389  // SPHERICAL CASE
390  // #################################
391 
392  // Boundary values
393 
394  Valeur boundary_r (mp.get_mg()-> get_angu()) ;
395  Valeur boundary_t (mp.get_mg()-> get_angu()) ;
396  Valeur boundary_p (mp.get_mg()-> get_angu()) ;
397 
398  switch (bound_beta) {
399 
400  case 0 : {
401  boundary_r = boundary_beta_r() ;
402  boundary_t = boundary_beta_theta() ;
403  boundary_p = boundary_beta_phi(omega) ;
404  break ;
405  }
406  case 1 : {
407  boundary_r = boundary_vv_x(omega) ;
408  boundary_t = boundary_vv_y(omega) ;
409  boundary_p = boundary_vv_z(omega) ;
410  break ;
411  }
412  default : {
413  cout <<"Unexpected type of boundary conditions for psi!"
414  << endl
415  << " bound_psi = " << bound_psi << endl ;
416  abort() ;
417  break ;
418  }
419  } // End of switch
420 
421  // Resolution
422  //-----------
423 
424  beta_jp1 = source_vector.poisson_dirichlet(lambda, boundary_r,
425  boundary_t, boundary_p, 0) ;
426 
427 
428  des_meridian(beta_jp1(1), 1.0000001, 10., "beta_r", 0) ;
429  des_meridian(beta_jp1(2), 1.0000001, 10., "beta_t", 1) ;
430  des_meridian(beta_jp1(3), 1.0000001, 10., "beta_p", 2) ;
431  arrete() ;
432  // #########################
433  // End of spherical case
434  // #########################
435 
436 
437 */
438  // Test
439  source_vector.dec_dzpuis() ;
440  maxabs(beta_jp1.derive_con(ff).divergence(ff)
441  + lambda * beta_jp1.divergence(ff)
442  .derive_con(ff) - source_vector,
443  "Absolute error in the resolution of the equation for beta") ;
444 
445  cout << endl ;
446 
447  // Boost
448  // -----
449 
450  Vector boost_vect(mp, CON, mp.get_bvect_cart()) ;
451  if (boost_x != 0.) {
452  boost_vect.set(1) = boost_x ;
453  boost_vect.set(2) = 0. ;
454  boost_vect.set(3) = 0. ;
455  boost_vect.std_spectral_base() ;
456  boost_vect.change_triad(mp.get_bvect_spher()) ;
457  beta_jp1 = beta_jp1 + boost_vect ;
458  }
459 
460  if (boost_z != 0.) {
461  boost_vect.set(1) = boost_z ;
462  boost_vect.set(2) = 0. ;
463  boost_vect.set(3) = 0. ;
464  boost_vect.std_spectral_base() ;
465  boost_vect.change_triad(mp.get_bvect_spher()) ;
466  beta_jp1 = beta_jp1 + boost_vect ;
467  }
468 
469  // Relaxation (relax=1 -> new ; relax=0 -> old )
470  beta_jp1 = relax_beta * beta_jp1 + (1 - relax_beta) * beta() ;
471  }
472 
473  //===========================================
474  // Convergence control
475  //===========================================
476 
477  double diff_nn, diff_psi, diff_beta ;
478  diff_nn = 1.e-16 ;
479  diff_psi = 1.e-16 ;
480  diff_beta = 1.e-16 ;
481  if (solve_lapse == 1)
482  diff_nn = max( diffrel(nn(), nn_jp1) ) ;
483  if (solve_psi == 1)
484  diff_psi = max( diffrel(psi(), psi_jp1) ) ;
485  if (solve_shift == 1)
486  diff_beta = max( maxabs(beta_jp1 - beta()) ) ;
487 
488  cout << "step = " << mer << " : diff_psi = " << diff_psi
489  << " diff_nn = " << diff_nn
490  << " diff_beta = " << diff_beta << endl ;
491  cout << "----------------------------------------------" << endl ;
492  if ((diff_psi<precis) && (diff_nn<precis) && (diff_beta<precis))
493  break ;
494 
495  if (mer>0) {conv << mer << " " << log10(diff_nn) << " " << log10(diff_psi)
496  << " " << log10(diff_beta) << endl ; } ;
497 
498  //=============================================
499  // Updates for next step
500  //=============================================
501 
502 
503  if (solve_psi == 1)
504  set_psi(psi_jp1) ;
505  if (solve_lapse == 1)
506  n_evol.update(nn_jp1, jtime, ttime) ;
507  if (solve_shift == 1)
508  beta_evol.update(beta_jp1, jtime, ttime) ;
509 
510  if (solve_shift == 1)
511  update_aa() ;
512 
513  // Saving ok K_{ij}s^is^j
514  // -----------------------
515 
516  Scalar kkss (contract(k_dd(), 0, 1, gam().radial_vect()*
517  gam().radial_vect(), 0, 1)) ;
518  double max_kss = kkss.val_grid_point(1, 0, 0, 0) ;
519  double min_kss = kkss.val_grid_point(1, 0, 0, 0) ;
520 
521  Scalar aaLss (pow(psi(), 6) * kkss) ;
522  double max_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
523  double min_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
524 
525  Scalar hh_tilde (contract(met_gamt.radial_vect().derive_cov(met_gamt), 0, 1)) ;
526  double max_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
527  double min_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
528 
529 
530  int nnp = mp.get_mg()->get_np(1) ;
531  int nnt = mp.get_mg()->get_nt(1) ;
532  for (int k=0 ; k<nnp ; k++)
533  for (int j=0 ; j<nnt ; j++){
534  if (kkss.val_grid_point(1, k, j, 0) > max_kss)
535  max_kss = kkss.val_grid_point(1, k, j, 0) ;
536  if (kkss.val_grid_point(1, k, j, 0) < min_kss)
537  min_kss = kkss.val_grid_point(1, k, j, 0) ;
538 
539  if (aaLss.val_grid_point(1, k, j, 0) > max_aaLss)
540  max_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
541  if (aaLss.val_grid_point(1, k, j, 0) < min_aaLss)
542  min_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
543 
544  if (hh_tilde.val_grid_point(1, k, j, 0) > max_hh_tilde)
545  max_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
546  if (hh_tilde.val_grid_point(1, k, j, 0) < min_hh_tilde)
547  min_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
548 
549  }
550 
551 
552  kss << mer << " " << max_kss << " " << min_kss << " " << max_aaLss << " " << min_aaLss
553  << " " << -1 * max_hh_tilde << " " << -1 * min_hh_tilde << endl ;
554  }
555 
556  conv.close() ;
557  kss.close() ;
558 
559 }
560 
561 
562 /*
563 
564 void Isol_hor::init_data_loop(int bound_nn, double lim_nn, int bound_psi,
565  int bound_beta, int solve_lapse, int solve_psi,
566  int solve_shift, double precis, double precis_loop,
567  double relax_nn, double relax_psi, double relax_beta, double relax_loop, int niter) {
568 
569  using namespace Unites ;
570 
571  // Initialisations
572  // ---------------
573  double ttime = the_time[jtime] ;
574 
575  ofstream conv("resconv.d") ;
576  ofstream kss("kss.d") ;
577  conv << " # diff_nn diff_psi diff_beta " << endl ;
578 
579  // Iteration
580  // ---------
581  for (int mer=0; mer<niter; mer++) {
582 
583  //=========================================================
584  // Boundary conditions and resolution of elliptic equations
585  //=========================================================
586 
587  // Resolution of the Poisson equation for the lapse
588  // ------------------------------------------------
589 
590 
591 
592  Scalar sou_nn (source_nn()) ;
593  Scalar nn_jp1 (mp) ;
594  if (solve_lapse == 1) {
595  Valeur nn_bound (mp.get_mg()-> get_angu()) ;
596 
597  switch (bound_nn) {
598 
599  case 0 : {
600  nn_bound = boundary_nn_Dir(lim_nn) ;
601  nn_jp1 = sou_nn.poisson_dirichlet(nn_bound, 0) + 1. ;
602  break ;
603  }
604  case 1 : {
605  nn_bound = boundary_nn_Neu_eff(lim_nn) ;
606  nn_jp1 = sou_nn.poisson_neumann(nn_bound, 0) + 1. ;
607  break ;
608  }
609  case 2 : {
610  nn_bound = boundary_nn_Dir_eff(lim_nn) ;
611  nn_jp1 = sou_nn.poisson_dirichlet(nn_bound, 0) + 1. ;
612  break ;
613  }
614  case 3 : {
615  nn_bound = boundary_nn_Neu_kk() ;
616  nn_jp1 = sou_nn.poisson_neumann(nn_bound, 0) + 1. ;
617  break ;
618  }
619  case 4 : {
620  nn_bound = boundary_nn_Dir_kk() ;
621  nn_jp1 = sou_nn.poisson_dirichlet(nn_bound, 0) + 1. ;
622  break ;
623  }
624  case 5 : {
625  nn_bound = boundary_nn_Dir_lapl(mer) ;
626  nn_jp1 = sou_nn.poisson_dirichlet(nn_bound, 0) + 1. ;
627  break ;
628  }
629  case 6 : {
630  nn_bound = boundary_nn_Neu_Cook() ;
631  nn_jp1 = sou_nn.poisson_neumann(nn_bound, 0) + 1. ;
632  break ;
633  }
634 
635 
636 
637 
638  default : {
639  cout <<"Unexpected type of boundary conditions for the lapse!"
640  << endl
641  << " bound_nn = " << bound_nn << endl ;
642  abort() ;
643  break ;
644  }
645 
646  } // End of switch
647 
648  // Test:
649  maxabs(nn_jp1.laplacian() - sou_nn,
650  "Absolute error in the resolution of the equation for N") ;
651 
652  // Relaxation (relax=1 -> new ; relax=0 -> old )
653  if (mer==0)
654  n_evol.update(nn_jp1, jtime, ttime) ;
655  else
656  nn_jp1 = relax_nn * nn_jp1 + (1 - relax_nn) * nn() ;
657  }
658 
659 
660  // Resolution of the Poisson equation for Psi
661  // ------------------------------------------
662 
663  Scalar sou_psi (source_psi()) ;
664  Scalar psi_jp1 (mp) ;
665  if (solve_psi == 1) {
666  Valeur psi_bound (mp.get_mg()-> get_angu()) ;
667 
668  switch (bound_psi) {
669 
670  case 0 : {
671  psi_bound = boundary_psi_app_hor() ;
672  psi_jp1 = sou_psi.poisson_neumann(psi_bound, 0) + 1. ;
673  break ;
674  }
675  case 1 : {
676  psi_bound = boundary_psi_Neu_spat() ;
677  psi_jp1 = sou_psi.poisson_neumann(psi_bound, 0) + 1. ;
678  break ;
679  }
680  case 2 : {
681  psi_bound = boundary_psi_Dir_spat() ;
682  psi_jp1 = sou_psi.poisson_dirichlet(psi_bound, 0) + 1. ;
683  break ;
684  }
685  case 3 : {
686  psi_bound = boundary_psi_Neu_evol() ;
687  psi_jp1 = sou_psi.poisson_neumann(psi_bound, 0) + 1. ;
688  break ;
689  }
690  case 4 : {
691  psi_bound = boundary_psi_Dir_evol() ;
692  psi_jp1 = sou_psi.poisson_dirichlet(psi_bound, 0) + 1. ;
693  break ;
694  }
695  case 5 : {
696  psi_bound = boundary_psi_Dir() ;
697  psi_jp1 = sou_psi.poisson_dirichlet(psi_bound, 0) + 1. ;
698  break ;
699  }
700  default : {
701  cout <<"Unexpected type of boundary conditions for psi!"
702  << endl
703  << " bound_psi = " << bound_psi << endl ;
704  abort() ;
705  break ;
706  }
707 
708  } // End of switch
709 
710  // Test:
711  maxabs(psi_jp1.laplacian() - sou_psi,
712  "Absolute error in the resolution of the equation for Psi") ;
713  // Relaxation (relax=1 -> new ; relax=0 -> old )
714  psi_jp1 = relax_psi * psi_jp1 + (1 - relax_psi) * psi() ;
715  }
716 
717  // Resolution of the vector Poisson equation for the shift
718  //---------------------------------------------------------
719 
720  // Source
721 
722  Vector beta_j(beta()) ;
723 
724  if (solve_shift == 1) {
725 
726  double lambda = 1./3.;
727  Vector beta_jp1 (beta()) ;
728  double thresh_loop = 1;
729  int n_loop = 0 ;
730 
731  while( thresh_loop > precis_loop ){
732 
733  Vector source_vector ( source_beta() ) ;
734  Vector source_reg = - (1./3. - lambda) * beta().divergence(ff)
735  .derive_con(ff) ;
736  source_reg.inc_dzpuis() ;
737  source_vector = source_vector + source_reg ;
738 
739 
740 
741  // Boundary values
742  // ===============
743 
744  Valeur boundary_x (mp.get_mg()-> get_angu()) ;
745  Valeur boundary_y (mp.get_mg()-> get_angu()) ;
746  Valeur boundary_z (mp.get_mg()-> get_angu()) ;
747 
748  switch (bound_beta) {
749 
750  case 0 : {
751  boundary_x = boundary_beta_x(omega) ;
752  boundary_y = boundary_beta_y(omega) ;
753  boundary_z = boundary_beta_z() ;
754  break ;
755  }
756  case 1 : {
757  boundary_x = boundary_vv_x(omega) ;
758  boundary_y = boundary_vv_y(omega) ;
759  boundary_z = boundary_vv_z(omega) ;
760  break ;
761  }
762  default : {
763  cout <<"Unexpected type of boundary conditions for beta!"
764  << endl
765  << " bound_beta = " << bound_beta << endl ;
766  abort() ;
767  break ;
768  }
769  } // End of switch
770 
771  if (boost_x != 0.)
772  boundary_x -= beta_boost_x() ;
773  if (boost_z != 0.)
774  boundary_z -= beta_boost_z() ;
775 
776  // Resolution
777  //-----------
778 
779  double precision = 1e-8 ;
780  poisson_vect_boundary(lambda, source_vector, beta_jp1, boundary_x,
781  boundary_y, boundary_z, 0, precision, 20) ;
782 
783  // Test
784  source_vector.dec_dzpuis() ;
785  maxabs(beta_jp1.derive_con(ff).divergence(ff)
786  + lambda * beta_jp1.divergence(ff)
787  .derive_con(ff) - source_vector,
788  "Absolute error in the resolution of the equation for beta") ;
789 
790  cout << endl ;
791 
792 
793 
794  // Relaxation_loop (relax=1 -> new ; relax=0 -> old )
795  beta_jp1 = relax_loop * beta_jp1 + (1 - relax_loop) * beta() ;
796 
797 
798  // Convergence loop
799  //=================
800 
801  double diff_beta_loop ;
802  diff_beta_loop = 1.e-16 ;
803  if (solve_shift == 1)
804  diff_beta_loop = max( maxabs(beta_jp1 - beta()) ) ;
805  cout << "step_loop = " << n_loop <<
806  << " diff_beta_loop = " << diff_beta_loop << endl ;
807  cout << "----------------------------------------------" << endl ;
808  thresh_loop = diff_beta_loop ;
809 
810  //Update loop
811  //===========
812  beta_evol.update(beta_jp1, jtime, ttime) ;
813  update_aa() ;
814  n_loop += 1 ;
815 
816  // End internal loop
817  }
818 
819  // Test for resolution of beta at this setp mer is already done in the internal loop
820 
821  // Relaxation beta (relax=1 -> new ; relax=0 -> old )
822  beta_jp1 = relax_beta * beta_jp1 + (1 - relax_loop) * beta_j ;
823 
824  }
825 
826 
827  //===========================================
828  // Convergence control
829  //===========================================
830 
831  double diff_nn, diff_psi, diff_beta ;
832  diff_nn = 1.e-16 ;
833  diff_psi = 1.e-16 ;
834  diff_beta = 1.e-16 ;
835  if (solve_lapse == 1)
836  diff_nn = max( diffrel(nn(), nn_jp1) ) ;
837  if (solve_psi == 1)
838  diff_psi = max( diffrel(psi(), psi_jp1) ) ;
839  if (solve_shift == 1)
840  diff_beta = max( maxabs(beta_jp1 - beta_j) ) ;
841 
842  cout << "step = " << mer << " : diff_psi = " << diff_psi
843  << " diff_nn = " << diff_nn
844  << " diff_beta = " << diff_beta << endl ;
845  cout << "----------------------------------------------" << endl ;
846  if ((diff_psi<precis) && (diff_nn<precis) && (diff_beta<precis))
847  break ;
848 
849  if (mer>0) {conv << mer << " " << log10(diff_nn) << " " << log10(diff_psi)
850  << " " << log10(diff_beta) << endl ; } ;
851 
852  //=============================================
853  // Updates for next step
854  //=============================================
855 
856 
857  if (solve_psi == 1)
858  set_psi(psi_jp1) ;
859  if (solve_lapse == 1)
860  n_evol.update(nn_jp1, jtime, ttime) ;
861  if (solve_shift == 1)
862  beta_evol.update(beta_jp1, jtime, ttime) ;
863 
864  if (solve_shift == 1)
865  update_aa() ;
866 
867  // Saving ok K_{ij}s^is^j
868  // -----------------------
869 
870  Scalar kkss (contract(k_dd(), 0, 1, gam().radial_vect()*
871  gam().radial_vect(), 0, 1)) ;
872  double max_kss = kkss.val_grid_point(1, 0, 0, 0) ;
873  double min_kss = kkss.val_grid_point(1, 0, 0, 0) ;
874 
875  Scalar aaLss (pow(psi(), 6) * kkss) ;
876  double max_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
877  double min_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
878 
879  Scalar hh_tilde (contract(met_gamt.radial_vect().derive_cov(met_gamt), 0, 1)) ;
880  double max_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
881  double min_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
882 
883 
884  int nnp = mp.get_mg()->get_np(1) ;
885  int nnt = mp.get_mg()->get_nt(1) ;
886  for (int k=0 ; k<nnp ; k++)
887  for (int j=0 ; j<nnt ; j++){
888  if (kkss.val_grid_point(1, k, j, 0) > max_kss)
889  max_kss = kkss.val_grid_point(1, k, j, 0) ;
890  if (kkss.val_grid_point(1, k, j, 0) < min_kss)
891  min_kss = kkss.val_grid_point(1, k, j, 0) ;
892 
893  if (aaLss.val_grid_point(1, k, j, 0) > max_aaLss)
894  max_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
895  if (aaLss.val_grid_point(1, k, j, 0) < min_aaLss)
896  min_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
897 
898  if (hh_tilde.val_grid_point(1, k, j, 0) > max_hh_tilde)
899  max_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
900  if (hh_tilde.val_grid_point(1, k, j, 0) < min_hh_tilde)
901  min_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
902 
903  }
904 
905 
906  kss << mer << " " << max_kss << " " << min_kss << " " << max_aaLss << " " << min_aaLss
907  << " " << -1 * max_hh_tilde << " " << -1 * min_hh_tilde << endl ;
908  }
909 
910  conv.close() ;
911  kss.close() ;
912 
913 }
914 
915 
916 */
917 
918 
919 
920 
921 
922 void Isol_hor::init_data_spher(int bound_nn, double lim_nn, int bound_psi,
923  int bound_beta, int solve_lapse, int solve_psi,
924  int solve_shift, double precis,
925  double relax, int niter) {
926 
927  using namespace Unites ;
928 
929  // Initialisations
930  // ---------------
931  double ttime = the_time[jtime] ;
932 
933  ofstream conv("resconv.d") ;
934  ofstream kss("kss.d") ;
935  conv << " # diff_nn diff_psi diff_beta " << endl ;
936 
937  // Iteration
938  // ---------
939  for (int mer=0; mer<niter; mer++) {
940 
941 
942  // des_meridian(psi(), 1, 10., "psi", 0) ;
943  // des_meridian(b_tilde(), 1, 10., "b_tilde", 1) ;
944  // des_meridian(nn(), 1, 10., "nn", 2) ;
945  // arrete() ;
946 
947 
948  //========
949  // Sources
950  //========
951 
952  // Useful functions
953  // ----------------
954  Vector tem_vect (beta() ) ;
955  Scalar dif_b = b_tilde() - tem_vect.set(1) ;
956  // cout << "dif_b = " << dif_b << endl ;
957  // arrete() ;
958 
959  Scalar dbdr ( b_tilde().dsdr() ) ;
960 
961  Scalar bsr (b_tilde()) ;
962  bsr.div_r() ;
963  bsr.inc_dzpuis(2) ;
964 
965  Scalar bsr2 ( bsr) ;
966  bsr2.div_r() ;
967  bsr2.inc_dzpuis(2) ;
968 
969  Scalar psisr (psi()) ;
970  psisr.div_r() ;
971  psisr.inc_dzpuis(2) ;
972 
973 
974 
975  // Source Psi
976  // ----------
977  Scalar source_psi_spher(mp) ;
978  source_psi_spher = -1./12. * psi4()*psi()/(nn() * nn()) * (dbdr - bsr) * (dbdr - bsr) ;
979 
980  // Source N
981  //---------
982  Scalar source_nn_spher(mp) ;
983  source_nn_spher = 2./3. * psi4() /nn() * (dbdr - bsr) * (dbdr - bsr)
984  - 2 * ln_psi().dsdr() * nn().dsdr() ;
985 
986  // Source b_tilde
987  //---------------
988  Scalar source_btilde_spher(mp) ;
989 
990  Scalar tmp ( -1./3. * (dbdr + 2 * bsr).dsdr() ) ;
991  tmp.std_spectral_base() ;
992  tmp.inc_dzpuis() ;
993 
994  source_btilde_spher = tmp + 2 * bsr2
995  + 4./3. * (dbdr - bsr) * ( nn().dsdr()/nn() - 6 * psi().dsdr()/psi() ) ;
996 
997  Scalar source_btilde_trun(mp) ;
998 
999  source_btilde_trun = tmp +
1000  4./3. * (dbdr - bsr) * ( nn().dsdr()/nn() - 6 * psi().dsdr()/psi() ) ;
1001 
1002 
1003  // Scalar diff_dbeta ( (dbdr + 2 * bsr).dsdr() - beta().divergence(ff).derive_con(ff)(1) ) ;
1004 
1005 
1006 
1007  // Parallel calculation
1008  //---------------------
1009 
1010  Scalar sourcepsi (source_psi()) ;
1011  Scalar sourcenn (source_nn()) ;
1012 
1013  Vector sourcebeta (source_beta()) ;
1014  Vector source_reg = 1./3. * beta().divergence(ff).derive_con(ff) ;
1015  source_reg.inc_dzpuis() ;
1016  sourcebeta -= source_reg ;
1017  Scalar source_btilde (sourcebeta(1) ) ;
1018 
1019  // Scalar diff_div = source_reg(1) + tmp ; ;
1020 
1021  Scalar mag_sou_psi ( source_psi_spher ) ;
1022  mag_sou_psi.dec_dzpuis(4) ;
1023  Scalar mag_sou_nn ( source_nn_spher ) ;
1024  mag_sou_nn.dec_dzpuis(4) ;
1025  Scalar mag_sou_btilde ( source_btilde_trun ) ;
1026  mag_sou_btilde.dec_dzpuis(4) ;
1027 
1028  Scalar diff_sou_psi ( source_psi_spher - sourcepsi) ;
1029  diff_sou_psi.dec_dzpuis(4) ;
1030  Scalar diff_sou_nn ( source_nn_spher - sourcenn) ;
1031  diff_sou_nn.dec_dzpuis(4) ;
1032  Scalar diff_sou_btilde ( source_btilde_trun - source_btilde) ;
1033  diff_sou_btilde.dec_dzpuis(4) ;
1034 
1035  /*
1036  cout << "dzpuis mag_btilde =" << mag_sou_btilde.get_dzpuis()<<endl ;
1037  des_meridian(diff_sou_psi, 1, 10., "diff_psi", 0) ;
1038  des_meridian(diff_sou_nn, 1, 10., "diff_nn", 1) ;
1039  des_meridian(diff_sou_btilde, 1, 10., "diff_btilde", 2) ;
1040  des_meridian(mag_sou_psi, 1, 10., "mag_psi", 3) ;
1041  des_meridian(mag_sou_nn, 1, 10., "mag_nn", 4) ;
1042  des_meridian(mag_sou_btilde, 1, 10., "mag_btilde", 5) ;
1043  // des_meridian(diff_dbeta, 1, 10., "diff_dbeta", 6) ;
1044 
1045  arrete() ;
1046  */
1047 
1048 
1049  //====================
1050  // Boundary conditions
1051  //====================
1052 
1053  // To avoid warnings;
1054  bound_nn = 1 ; lim_nn = 1. ; bound_psi = 1 ; bound_beta = 1 ;
1055 
1056  double kappa_0 = 0.2 - 1. ;
1057 
1058  Scalar kappa (mp) ;
1059  kappa = kappa_0 ;
1060  kappa.std_spectral_base() ;
1061  kappa.inc_dzpuis(2) ;
1062 
1063 
1064  int nnp = mp.get_mg()->get_np(1) ;
1065  int nnt = mp.get_mg()->get_nt(1) ;
1066 
1067 
1068  Valeur psi_bound (mp.get_mg()-> get_angu()) ;
1069  Valeur nn_bound (mp.get_mg()-> get_angu()) ;
1070  Valeur btilde_bound (mp.get_mg()-> get_angu()) ;
1071  psi_bound = 1. ; // Juste pour affecter dans espace des configs ;
1072  nn_bound = 1. ; // Juste pour affecter dans espace des configs ;
1073  btilde_bound = 1. ; // Juste pour affecter dans espace des configs ;
1074 
1075  Scalar tmp_psi = -1./4. * (2 * psisr +
1076  2./3. * psi4()/(psi() * nn()) * (dbdr - bsr) ) ;
1077 
1078  Scalar tmp_nn = kappa ; //+ 2./3. * psi() * psi() * (dbdr - bsr) ;
1079 
1080  Scalar tmp_btilde = nn() / (psi() * psi()) ;
1081 
1082 
1083  for (int k=0 ; k<nnp ; k++)
1084  for (int j=0 ; j<nnt ; j++){
1085  psi_bound.set(0, k, j, 0) = tmp_psi.val_grid_point(1, k, j, 0) ; // BC Psi
1086  nn_bound.set(0, k, j, 0) = tmp_nn.val_grid_point(1, k, j, 0) ; // BC N
1087  btilde_bound.set(0, k, j, 0) = tmp_btilde.val_grid_point(1, k, j, 0) ; // BC b_tilde
1088  }
1089 
1090  psi_bound.std_base_scal() ;
1091  nn_bound.std_base_scal() ;
1092  btilde_bound.std_base_scal() ;
1093 
1094 
1095  //=================================
1096  // Resolution of elliptic equations
1097  //=================================
1098 
1099  // Resolution of the Poisson equation for Psi
1100  // ------------------------------------------
1101  Scalar psi_jp1 (mp) ;
1102  if (solve_psi == 1) {
1103 
1104  psi_jp1 = source_psi_spher.poisson_neumann(psi_bound, 0) + 1. ;
1105 
1106  // Test:
1107  maxabs(psi_jp1.laplacian() - source_psi_spher,
1108  "Absolute error in the resolution of the equation for Psi") ;
1109  // Relaxation (relax=1 -> new ; relax=0 -> old )
1110  psi_jp1 = relax * psi_jp1 + (1 - relax) * psi() ;
1111  }
1112 
1113  // Resolution of the Poisson equation for the lapse
1114  // ------------------------------------------------
1115  Scalar nn_jp1 (mp) ;
1116  if (solve_lapse == 1) {
1117 
1118  nn_jp1 = source_nn_spher.poisson_dirichlet(nn_bound, 0) + 1. ;
1119 
1120  // Test:
1121  maxabs(nn_jp1.laplacian() - source_nn_spher,
1122  "Absolute error in the resolution of the equation for N") ;
1123 
1124  // Relaxation (relax=1 -> new ; relax=0 -> old )
1125  if (mer==0)
1126  n_evol.update(nn_jp1, jtime, ttime) ;
1127  else
1128  nn_jp1 = relax * nn_jp1 + (1 - relax) * nn() ;
1129 
1130  }
1131 
1132  // Resolution of the Poisson equation for b_tilde
1133  // ----------------------------------------------
1134  Scalar btilde_jp1 (mp) ;
1135  if (solve_shift == 1) {
1136 
1137  btilde_jp1 = source_btilde_spher.poisson_dirichlet(btilde_bound, 0) ;
1138 
1139  // Test:
1140  maxabs(btilde_jp1.laplacian() - source_btilde_spher,
1141  "Absolute error in the resolution of the equation for btilde") ;
1142  // Relaxation (relax=1 -> new ; relax=0 -> old )
1143  btilde_jp1 = relax * btilde_jp1 + (1 - relax) * b_tilde() ;
1144  }
1145 
1146 
1147  //===========================================
1148  // Convergence control
1149  //===========================================
1150 
1151  double diff_nn, diff_psi, diff_btilde ;
1152  diff_nn = 1.e-16 ;
1153  diff_psi = 1.e-16 ;
1154  diff_btilde = 1.e-16 ;
1155  if (solve_lapse == 1)
1156  diff_nn = max( diffrel(nn(), nn_jp1) ) ;
1157  if (solve_psi == 1)
1158  diff_psi = max( diffrel(psi(), psi_jp1) ) ;
1159  if (solve_shift == 1)
1160  diff_btilde = max( diffrel(btilde_jp1, b_tilde()) ) ;
1161 
1162  cout << "step = " << mer << " : diff_psi = " << diff_psi
1163  << " diff_nn = " << diff_nn
1164  << " diff_btilde = " << diff_btilde << endl ;
1165  cout << "----------------------------------------------" << endl ;
1166  if ((diff_psi<precis) && (diff_nn<precis) && (diff_btilde<precis))
1167  break ;
1168 
1169  if (mer>0) {conv << mer << " " << log10(diff_nn) << " " << log10(diff_psi)
1170  << " " << log10(diff_btilde) << endl ; } ;
1171 
1172  //=============================================
1173  // Updates for next step
1174  //=============================================
1175 
1176 
1177  if (solve_psi == 1)
1178  set_psi(psi_jp1) ;
1179  if (solve_lapse == 1)
1180  n_evol.update(nn_jp1, jtime, ttime) ;
1181  if (solve_shift == 1)
1182  {
1183  Vector beta_jp1 (btilde_jp1 * tgam().radial_vect()) ;
1184  cout << tgam().radial_vect() << endl ;
1185  beta_evol.update(beta_jp1, jtime, ttime) ;
1186  }
1187  if (solve_shift == 1 || solve_lapse == 1)
1188  {
1189  update_aa() ;
1190  }
1191 
1192  // Saving ok K_{ij}s^is^j
1193  // -----------------------
1194 
1195  Scalar kkss (contract(k_dd(), 0, 1, gam().radial_vect()*
1196  gam().radial_vect(), 0, 1)) ;
1197  double max_kss = kkss.val_grid_point(1, 0, 0, 0) ;
1198  double min_kss = kkss.val_grid_point(1, 0, 0, 0) ;
1199 
1200  Scalar aaLss (pow(psi(), 6) * kkss) ;
1201  double max_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
1202  double min_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
1203 
1204  Scalar hh_tilde (contract(met_gamt.radial_vect().derive_cov(met_gamt), 0, 1)) ;
1205  double max_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
1206  double min_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
1207 
1208 
1209 
1210  for (int k=0 ; k<nnp ; k++)
1211  for (int j=0 ; j<nnt ; j++){
1212  if (kkss.val_grid_point(1, k, j, 0) > max_kss)
1213  max_kss = kkss.val_grid_point(1, k, j, 0) ;
1214  if (kkss.val_grid_point(1, k, j, 0) < min_kss)
1215  min_kss = kkss.val_grid_point(1, k, j, 0) ;
1216 
1217  if (aaLss.val_grid_point(1, k, j, 0) > max_aaLss)
1218  max_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
1219  if (aaLss.val_grid_point(1, k, j, 0) < min_aaLss)
1220  min_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
1221 
1222  if (hh_tilde.val_grid_point(1, k, j, 0) > max_hh_tilde)
1223  max_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
1224  if (hh_tilde.val_grid_point(1, k, j, 0) < min_hh_tilde)
1225  min_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
1226 
1227  }
1228 
1229 
1230  kss << mer << " " << max_kss << " " << min_kss << " " << max_aaLss << " " << min_aaLss
1231  << " " << -1 * max_hh_tilde << " " << -1 * min_hh_tilde << endl ;
1232  }
1233 
1234  conv.close() ;
1235  kss.close() ;
1236 
1237 }
1238 
1239 
1240 
1241 void Isol_hor::init_data_alt(int, double, int,
1242  int, int solve_lapse, int solve_psi,
1243  int solve_shift, double precis,
1244  double relax, int niter) {
1245 
1246  using namespace Unites ;
1247 
1248  // Initialisations
1249  // ---------------
1250  double ttime = the_time[jtime] ;
1251 
1252  ofstream conv("resconv.d") ;
1253  ofstream kss("kss.d") ;
1254  conv << " # diff_nn diff_psi diff_beta " << endl ;
1255 
1256  Scalar psi_j (psi()) ;
1257  Scalar nn_j (nn()) ;
1258  Scalar btilde_j (b_tilde()) ;
1259  Scalar diffb ( btilde_j - b_tilde() ) ;
1260  Scalar theta_j (beta().divergence(ff)) ;
1261  theta_j.dec_dzpuis(2) ;
1262  Scalar chi_j (b_tilde()) ;
1263  chi_j.mult_r() ;
1264 
1265 
1266  // Iteration
1267  // ---------
1268 
1269  for (int mer=0; mer<niter; mer++) {
1270 
1271 
1272  des_meridian(psi_j, 1, 10., "psi", 0) ;
1273  des_meridian(nn_j, 1, 10., "nn", 1) ;
1274  des_meridian(theta_j, 1, 10., "Theta", 2) ;
1275  des_meridian(chi_j, 1, 10., "chi", 3) ;
1276  arrete() ;
1277 
1278 
1279  //========
1280  // Sources
1281  //========
1282 
1283  // Useful functions
1284  // ----------------
1285 
1286  Scalar psisr (psi_j) ;
1287  psisr.div_r() ;
1288  psisr.inc_dzpuis(2) ;
1289 
1290  Scalar dchidr ( chi_j.dsdr() ) ;
1291 
1292  Scalar chisr (chi_j) ;
1293  chisr.div_r() ;
1294  chisr.inc_dzpuis(2) ;
1295 
1296  Scalar rdthetadr (theta_j.dsdr() ) ;
1297  rdthetadr.mult_r() ;
1298  rdthetadr.inc_dzpuis(2) ;
1299 
1300  Scalar theta_dz4 (theta_j) ;
1301  theta_dz4.inc_dzpuis(4) ;
1302 
1303  Scalar dbmb (dchidr - 2 * chisr) ;
1304  dbmb.div_r() ;
1305 
1306 
1307  // Source Psi
1308  // ----------
1309  Scalar source_psi_spher(mp) ;
1310 
1311  source_psi_spher = -1./12. * psi_j*psi_j*psi_j*psi_j*psi_j/(nn_j * nn_j)
1312  * dbmb *dbmb ;
1313 
1314 
1315  // Source N
1316  //---------
1317  Scalar source_nn_spher(mp) ;
1318  source_nn_spher = 2./3. * psi_j*psi_j*psi_j*psi_j/nn_j * dbmb *dbmb
1319  - 2 * psi_j.dsdr()/psi_j * nn_j.dsdr() ;
1320 
1321 
1322  //====================
1323  // Boundary conditions
1324  //====================
1325  double kappa_0 = 0.2 - 1. ;
1326 
1327  Scalar kappa (mp) ;
1328  kappa = kappa_0 ;
1329  kappa.std_spectral_base() ;
1330  kappa.inc_dzpuis(2) ;
1331 
1332  int nnp = mp.get_mg()->get_np(1) ;
1333  int nnt = mp.get_mg()->get_nt(1) ;
1334 
1335  Valeur psi_bound (mp.get_mg()-> get_angu()) ;
1336  Valeur nn_bound (mp.get_mg()-> get_angu()) ;
1337 
1338  psi_bound = 1. ; // Juste pour affecter dans espace des configs ;
1339  nn_bound = 1. ; // Juste pour affecter dans espace des configs ;
1340 
1341  //psi
1342 
1343  Scalar tmp_psi = -1./4. * (2 * psisr +
1344  2./3. * psi_j*psi_j*psi_j/ nn_j * dbmb ) ;
1345 
1346  // tmp_psi = 2./3. * psi_j*psi_j*psi_j/ nn_j * (dchidr - 2 * chisr) ;
1347  // tmp_psi.div_r() ;
1348  // tmp_psi = -1./4. * (2 * psisr + tmp_psi) ;
1349 
1350  //nn
1351  Scalar tmp_nn = kappa ; //+ 2./3. * psi_j*psi_j * dbmb ;
1352 
1353 
1354 
1355  for (int k=0 ; k<nnp ; k++)
1356  for (int j=0 ; j<nnt ; j++){
1357  psi_bound.set(0, k, j, 0) = tmp_psi.val_grid_point(1, k, j, 0) ; // BC Psi
1358  nn_bound.set(0, k, j, 0) = tmp_nn.val_grid_point(1, k, j, 0) ; // BC N
1359  }
1360 
1361  psi_bound.std_base_scal() ;
1362  nn_bound.std_base_scal() ;
1363 
1364 
1365 
1366  //=================================
1367  // Resolution of elliptic equations
1368  //=================================
1369 
1370  // Resolution of the Poisson equation for Psi
1371  // ------------------------------------------
1372  Scalar psi_jp1 (mp) ;
1373  if (solve_psi == 1) {
1374 
1375  psi_jp1 = source_psi_spher.poisson_neumann(psi_bound, 0) + 1. ;
1376 
1377  // Test:
1378  maxabs(psi_jp1.laplacian() - source_psi_spher,
1379  "Absolute error in the resolution of the equation for Psi") ;
1380  // Relaxation (relax=1 -> new ; relax=0 -> old )
1381  psi_jp1 = relax * psi_jp1 + (1 - relax) * psi_j ;
1382  }
1383 
1384  // Resolution of the Poisson equation for the lapse
1385  // ------------------------------------------------
1386  Scalar nn_jp1 (mp) ;
1387  if (solve_lapse == 1) {
1388 
1389  nn_jp1 = source_nn_spher.poisson_dirichlet(nn_bound, 0) + 1. ;
1390 
1391  // Test:
1392  maxabs(nn_jp1.laplacian() - source_nn_spher,
1393  "Absolute error in the resolution of the equation for N") ;
1394 
1395  // Relaxation (relax=1 -> new ; relax=0 -> old )
1396  if (mer==0)
1397  n_evol.update(nn_jp1, jtime, ttime) ;
1398  else
1399  nn_jp1 = relax * nn_jp1 + (1 - relax) * nn_j ;
1400 
1401  }
1402 
1403 
1404  // Resolution for chi and Theta
1405  // ----------------------------
1406  Scalar theta_jp1 (mp) ;
1407  Scalar chi_jp1 (mp) ;
1408 
1409  if (solve_shift == 1) {
1410 
1411  // Initialisations loop on theta/chi
1412  Scalar theta_i(theta_j) ;
1413  Scalar chi_i(chi_j) ;
1414 
1415 
1416  // Iteration in theta/chi
1417  for (int i=0 ; i<niter ; i++) {
1418 
1419  des_meridian(theta_i, 1, 10., "Theta", 2) ;
1420  des_meridian(chi_i, 1, 10., "chi", 3) ;
1421  arrete() ;
1422 
1423 
1424 
1425  //Sources
1426 
1427  // Source_theta
1428  //-------------
1429  Scalar source_theta_spher(mp) ;
1430  source_theta_spher = (dbmb * (nn_j.dsdr()/nn_j - 6 * psi_j.dsdr()/psi_j)).dsdr() ;
1431  source_theta_spher.dec_dzpuis() ;
1432 
1433  // Source chi
1434  //-----------
1435  Scalar source_chi_spher(mp) ;
1436  source_chi_spher = 4./3. * (dchidr - 2 * chisr) * ( nn_j.dsdr()/nn_j - 6 * psi_j.dsdr()/psi_j )
1437  - 1./3. * rdthetadr + 2 * theta_dz4 ;
1438 
1439  //Boundaries
1440  Valeur theta_bound (mp.get_mg()-> get_angu()) ;
1441  Valeur chi_bound (mp.get_mg()-> get_angu()) ;
1442 
1443  theta_bound = 1. ; // Juste pour affecter dans espace des configs ;
1444  chi_bound = 1. ; // Juste pour affecter dans espace des configs ;
1445 
1446  //theta
1447  Scalar tmp_theta = dchidr ;
1448  tmp_theta.dec_dzpuis(2) ;
1449  tmp_theta += nn_j/(psi_j*psi_j) ;
1450  tmp_theta.div_r() ;
1451 
1452  //chi
1453  Scalar tmp_chi = nn_j/(psi_j*psi_j) ;
1454  tmp_chi.mult_r() ;
1455 
1456  for (int k=0 ; k<nnp ; k++)
1457  for (int j=0 ; j<nnt ; j++){
1458  theta_bound.set(0, k, j, 0) = tmp_theta.val_grid_point(1, k, j, 0) ; // BC Theta
1459  chi_bound.set(0, k, j, 0) = tmp_chi.val_grid_point(1, k, j, 0) ; // BC chi
1460  }
1461  theta_bound.std_base_scal() ;
1462  chi_bound.std_base_scal() ;
1463 
1464  //Resolution equations
1465  Scalar theta_ip1(mp) ;
1466  Scalar chi_ip1(mp) ;
1467 
1468  theta_ip1 = source_theta_spher.poisson_dirichlet(theta_bound, 0) ;
1469  chi_ip1 = source_chi_spher.poisson_dirichlet(chi_bound, 0) ;
1470 
1471  // Test:
1472  maxabs(theta_ip1.laplacian() - source_theta_spher,
1473  "Absolute error in the resolution of the equation for Theta") ;
1474  maxabs(chi_ip1.laplacian() - source_chi_spher,
1475  "Absolute error in the resolution of the equation for chi") ;
1476 
1477  // Relaxation (relax=1 -> new ; relax=0 -> old )
1478  theta_ip1 = relax * theta_ip1 + (1 - relax) * theta_i ;
1479  chi_ip1 = relax * chi_ip1 + (1 - relax) * chi_i ;
1480 
1481  // Convergence control of loop in theta/chi
1482  double diff_theta_int, diff_chi_int, int_precis ;
1483  diff_theta_int = 1.e-16 ;
1484  diff_chi_int = 1.e-16 ;
1485  int_precis = 1.e-3 ;
1486 
1487  diff_theta_int = max( diffrel(theta_ip1, theta_i) ) ;
1488  diff_chi_int = max( diffrel(chi_ip1, chi_i) ) ;
1489 
1490 
1491  cout << "internal step = " << i
1492  << " diff_theta_int = " << diff_theta_int
1493  << " diff_chi_int = " << diff_chi_int << endl ;
1494  cout << "----------------------------------------------" << endl ;
1495  if ((diff_theta_int<int_precis) && (diff_chi_int<int_precis))
1496  {
1497  theta_jp1 = theta_ip1 ;
1498  chi_jp1 = chi_ip1 ;
1499  break ;
1500  }
1501  // Updates of internal loop in theta/chi
1502  theta_i = theta_ip1 ;
1503  chi_i = chi_ip1 ;
1504  }
1505  }
1506 
1507 
1508  //===========================================
1509  // Convergence control
1510  //===========================================
1511 
1512  double diff_nn, diff_psi, diff_theta, diff_chi ;
1513  diff_nn = 1.e-16 ;
1514  diff_psi = 1.e-16 ;
1515  diff_theta = 1.e-16 ;
1516  diff_chi = 1.e-16 ;
1517 
1518  if (solve_lapse == 1)
1519  diff_nn = max( diffrel(nn_j, nn_jp1) ) ;
1520  if (solve_psi == 1)
1521  diff_psi = max( diffrel(psi_j, psi_jp1) ) ;
1522  if (solve_shift == 1)
1523  diff_theta = max( diffrel(theta_jp1, theta_j) ) ;
1524  if (solve_shift == 1)
1525  diff_chi = max( diffrel(chi_jp1, chi_j) ) ;
1526 
1527 
1528  cout << "step = " << mer << " : diff_psi = " << diff_psi
1529  << " diff_nn = " << diff_nn
1530  << " diff_theta = " << diff_theta
1531  << " diff_chi = " << diff_chi << endl ;
1532  cout << "----------------------------------------------" << endl ;
1533  if ((diff_psi<precis) && (diff_nn<precis) && (diff_theta<precis) && (diff_chi<precis))
1534  break ;
1535 
1536  if (mer>0) {conv << mer << " " << log10(diff_nn) << " " << log10(diff_psi)
1537  << " " << log10(diff_theta)
1538  << " " << log10(diff_chi) << endl ; } ;
1539 
1540  //=============================================
1541  // Updates for next step
1542  //=============================================
1543 
1544 
1545  if (solve_psi == 1)
1546  set_psi(psi_jp1) ;
1547  psi_j = psi_jp1 ;
1548  if (solve_lapse == 1)
1549  n_evol.update(nn_jp1, jtime, ttime) ;
1550  nn_j = nn_jp1 ;
1551  if (solve_shift == 1)
1552  {
1553  theta_j = theta_jp1 ;
1554  chi_j = chi_jp1 ;
1555  chi_jp1.mult_r() ;
1556  Vector beta_jp1 (chi_jp1 * tgam().radial_vect()) ;
1557  beta_evol.update(beta_jp1, jtime, ttime) ;
1558  }
1559 
1560  if (solve_shift == 1 || solve_lapse == 1)
1561  {
1562  update_aa() ;
1563  }
1564 
1565  // Saving ok K_{ij}s^is^j
1566  // -----------------------
1567 
1568  Scalar kkss (contract(k_dd(), 0, 1, gam().radial_vect()*
1569  gam().radial_vect(), 0, 1)) ;
1570  double max_kss = kkss.val_grid_point(1, 0, 0, 0) ;
1571  double min_kss = kkss.val_grid_point(1, 0, 0, 0) ;
1572 
1573  Scalar aaLss (pow(psi(), 6) * kkss) ;
1574  double max_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
1575  double min_aaLss = aaLss.val_grid_point(1, 0, 0, 0) ;
1576 
1577  Scalar hh_tilde (contract(met_gamt.radial_vect().derive_cov(met_gamt), 0, 1)) ;
1578  double max_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
1579  double min_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
1580 
1581 
1582 
1583  for (int k=0 ; k<nnp ; k++)
1584  for (int j=0 ; j<nnt ; j++){
1585  if (kkss.val_grid_point(1, k, j, 0) > max_kss)
1586  max_kss = kkss.val_grid_point(1, k, j, 0) ;
1587  if (kkss.val_grid_point(1, k, j, 0) < min_kss)
1588  min_kss = kkss.val_grid_point(1, k, j, 0) ;
1589 
1590  if (aaLss.val_grid_point(1, k, j, 0) > max_aaLss)
1591  max_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
1592  if (aaLss.val_grid_point(1, k, j, 0) < min_aaLss)
1593  min_aaLss = aaLss.val_grid_point(1, k, j, 0) ;
1594 
1595  if (hh_tilde.val_grid_point(1, k, j, 0) > max_hh_tilde)
1596  max_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
1597  if (hh_tilde.val_grid_point(1, k, j, 0) < min_hh_tilde)
1598  min_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
1599 
1600  }
1601 
1602 
1603  kss << mer << " " << max_kss << " " << min_kss << " " << max_aaLss << " " << min_aaLss
1604  << " " << -1 * max_hh_tilde << " " << -1 * min_hh_tilde << endl ;
1605  }
1606 
1607  conv.close() ;
1608  kss.close() ;
1609 
1610 }
1611 
1612 
1613 
1614 }
Lorene::Time_slice_conf::k_dd
virtual const Sym_tensor & k_dd() const
Extrinsic curvature tensor (covariant components ) at the current time step (jtime )
Definition: time_slice_conf.C:630
Lorene::Isol_hor::source_beta
const Vector source_beta() const
Source for .
Definition: sources_hor.C:190
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
Lorene::Isol_hor::boundary_nn_Neu_Cook
const Valeur boundary_nn_Neu_Cook() const
Neumann boundary condition for N using Cook's boundary condition.
Definition: bound_hor.C:489
Lorene::Tensor::derive_cov
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition: tensor.C:1002
Lorene::Isol_hor::boundary_beta_z
const Valeur boundary_beta_z() const
Component z of boundary value of .
Definition: bound_hor.C:1121
Lorene::Time_slice_conf::nn
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )
Definition: time_slice_conf.C:591
Lorene::Time_slice_conf::ln_psi
const Scalar & ln_psi() const
Logarithm of at the current time step (jtime ).
Definition: time_slice_conf.C:719
Lorene::Time_slice::beta_evol
Evolution_std< Vector > beta_evol
Values at successive time steps of the shift vector .
Definition: time_slice.h:219
Lorene::Scalar::dsdr
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
Lorene::Isol_hor::source_psi
const Scalar source_psi() const
Source for .
Definition: sources_hor.C:108
Lorene::Isol_hor::boundary_vv_z
const Valeur boundary_vv_z(double om) const
Component z of boundary value of .
Definition: bound_hor.C:1547
Lorene::exp
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
Lorene::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Isol_hor::boundary_nn_Dir
const Valeur boundary_nn_Dir(double aa) const
Dirichlet boundary condition .
Definition: bound_hor.C:647
Lorene::Time_slice::n_evol
Evolution_std< Scalar > n_evol
Values at successive time steps of the lapse function N.
Definition: time_slice.h:216
Lorene::Map::get_bvect_spher
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::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Isol_hor::boundary_nn_Dir_kk
const Valeur boundary_nn_Dir_kk() const
Dirichlet boundary condition for N using the extrinsic curvature.
Definition: bound_hor.C:371
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:435
Lorene::Vector::divergence
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:381
Lorene::arrete
void arrete(int a=0)
Setting a stop point in a code.
Definition: arrete.C:61
Lorene::Isol_hor::boundary_psi_Dir_spat
const Valeur boundary_psi_Dir_spat() const
Dirichlet boundary condition for (spatial)
Definition: bound_hor.C:231
Lorene::Isol_hor::boundary_vv_x
const Valeur boundary_vv_x(double om) const
Component x of boundary value of .
Definition: bound_hor.C:1493
Lorene::Isol_hor::mp
Map_af & mp
Affine mapping.
Definition: isol_hor.h:260
Lorene::Isol_hor::omega
double omega
Angular velocity in LORENE's units.
Definition: isol_hor.h:269
Lorene::Isol_hor::boundary_nn_Neu_eff
const Valeur boundary_nn_Neu_eff(double aa) const
Neumann boundary condition on nn (effectif)
Definition: bound_hor.C:622
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Lorene::Isol_hor::boundary_nn_Dir_eff
const Valeur boundary_nn_Dir_eff(double aa) const
Dirichlet boundary condition for N (effectif)
Definition: bound_hor.C:586
Lorene::Isol_hor::set_psi
void set_psi(const Scalar &psi_in)
Sets the conformal factor relating the physical metric to the conformal one: .
Definition: isol_hor.C:515
Lorene::Isol_hor::b_tilde
const Scalar b_tilde() const
Radial component of the shift with respect to the conformal metric.
Definition: phys_param.C:136
Lorene::Isol_hor::boundary_nn_Neu_kk
const Valeur boundary_nn_Neu_kk(int nn=1) const
Neumann boundary condition for N using the extrinsic curvature.
Definition: bound_hor.C:420
Unites
Standard units of space, time and mass.
Lorene::Isol_hor::met_gamt
Metric met_gamt
3 metric tilde
Definition: isol_hor.h:326
Lorene::Isol_hor::update_aa
void update_aa()
Conformal representation of the traceless part of the extrinsic curvature: .
Definition: isol_hor.C:842
Lorene::Tensor::set
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:654
Lorene::Time_slice_conf::psi
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
Definition: time_slice_conf.C:693
Lorene::Isol_hor::boundary_vv_y
const Valeur boundary_vv_y(double om) const
Component y of boundary value of .
Definition: bound_hor.C:1521
Lorene::Isol_hor::boundary_psi_Dir_evol
const Valeur boundary_psi_Dir_evol() const
Dirichlet boundary condition for (evolution)
Definition: bound_hor.C:174
Lorene::Isol_hor::boost_x
double boost_x
Boost velocity in x-direction.
Definition: isol_hor.h:272
Lorene::Isol_hor::boundary_psi_Neu_spat
const Valeur boundary_psi_Neu_spat() const
Neumann boundary condition for (spatial)
Definition: bound_hor.C:267
Lorene::Isol_hor::boundary_psi_Dir
const Valeur boundary_psi_Dir() const
Dirichlet boundary condition for (spatial)
Definition: bound_hor.C:324
Lorene::Isol_hor::boundary_nn_Dir_lapl
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:693
Lorene::Isol_hor::tgam
virtual const Metric & tgam() const
Conformal metric Returns the value at the current time step (jtime ).
Definition: isol_hor.h:514
Lorene::log10
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition: cmp_math.C:322
Lorene::Isol_hor::source_nn
const Scalar source_nn() const
Source for N.
Definition: sources_hor.C:145
Lorene::Time_slice::the_time
Evolution_std< double > the_time
Time label of each slice.
Definition: time_slice.h:193
Lorene::Map::get_bvect_cart
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Definition: map.h:791
des_meridian
void des_meridian(const Scalar &uu, double r_min, double r_max, const char *nomy, int ngraph, const char *device=0x0, bool closeit=false, bool draw_bound=true)
Draws 5 profiles of a scalar field along various radial axes in two meridional planes and .
Lorene::maxabs
Tbl maxabs(const Tensor &aa, const char *comment=0x0, ostream &ost=cout, bool verb=true)
Maxima in each domain of the absolute values of the tensor components.
Definition: tensor_calculus_ext.C:606
Lorene::Time_slice_conf::ff
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
Definition: time_slice.h:507
Lorene::Isol_hor::boost_z
double boost_z
Boost velocity in z-direction.
Definition: isol_hor.h:275
Lorene::Time_slice_conf::psi4
const Scalar & psi4() const
Factor at the current time step (jtime ).
Definition: time_slice_conf.C:707
Lorene::Metric::radial_vect
virtual const Vector & radial_vect() const
Returns the radial vector normal to a spherical slicing and pointing toward spatial infinity.
Definition: metric.C:362
Lorene::Scalar::derive_con
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
Lorene::Isol_hor::boundary_psi_app_hor
const Valeur boundary_psi_app_hor() const
Neumann boundary condition for (spatial)
Definition: bound_hor.C:297
Lorene::Tensor::inc_dzpuis
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
Lorene::Time_slice::jtime
int jtime
Time step index of the latest slice.
Definition: time_slice.h:190
Lorene::Isol_hor::beta_boost_z
const Valeur beta_boost_z() const
Boundary value for a boost in z-direction.
Definition: bound_hor.C:1155
Lorene::Isol_hor::boundary_beta_x
const Valeur boundary_beta_x(double om) const
Component x of boundary value of .
Definition: bound_hor.C:1021
Lorene::Time_slice::beta
virtual const Vector & beta() const
shift vector at the current time step (jtime )
Definition: time_slice_access.C:87
Lorene::Isol_hor::beta_boost_x
const Valeur beta_boost_x() const
Boundary value for a boost in x-direction.
Definition: bound_hor.C:1144
Lorene::Scalar::val_grid_point
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:637
Lorene::Isol_hor::boundary_psi_Neu_evol
const Valeur boundary_psi_Neu_evol() const
Neumann boundary condition for (evolution)
Definition: bound_hor.C:203
Lorene::Isol_hor::boundary_beta_y
const Valeur boundary_beta_y(double om) const
Component y of boundary value of .
Definition: bound_hor.C:1071
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition: tenseur_operateur.C:279
Lorene::Time_slice::gam
const Metric & gam() const
Induced metric at the current time step (jtime )
Definition: time_slice_access.C:95