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