LORENE
CTS_gen.C
1 /*
2  * Method of class Isol_hor to compute intital data using generalized
3  * Conformal Thni Sandwich equations with N = Ntilde \Psi ^a
4  * and K_ij = \Psi ^\zeta A_ij + 1/3 K \gamma_ij
5  *
6  * (see file isol_hor.h for documentation).
7  *
8  */
9
10 /*
11  * Copyright (c) 2004 Jose Luis Jaramillo
12  * Francois Limousin
13  *
14  * This file is part of LORENE.
15  *
16  * LORENE is free software; you can redistribute it and/or modify
17  * it under the terms of the GNU General Public License version 2
18  * as published by the Free Software Foundation.
19  *
20  * LORENE is distributed in the hope that it will be useful,
21  * but WITHOUT ANY WARRANTY; without even the implied warranty of
22  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23  * GNU General Public License for more details.
24  *
25  * You should have received a copy of the GNU General Public License
26  * along with LORENE; if not, write to the Free Software
27  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
28  *
29  */
30
31 char CTS_gen[] = "$Header: /cvsroot/Lorene/C++/Source/Isol_hor/CTS_gen.C,v 1.7 2014/10/13 08:53:00 j_novak Exp$" ;
32
33 // Lorene headers
34 #include "time_slice.h"
35 #include "isol_hor.h"
36 #include "metric.h"
37 #include "evolution.h"
38 #include "unites.h"
39 #include "graphique.h"
40 #include "utilitaires.h"
41 #include "param.h"
42
43 namespace Lorene {
44 void Isol_hor::init_data_CTS_gen(int, double, int, int,
45  int solve_lapse, int solve_psi, int solve_shift,
46  double precis, double relax_nn ,
47  double relax_psi, double relax_beta ,
48  int niter, double a, double ) {
49
50  using namespace Unites ;
51
52  // Initialisations
53  // ===============
54
55  double lambda = 0. ; // to be used in the beta source
56
57  double ttime = the_time[jtime] ;
58  const Base_vect& triad = *(ff.get_triad()) ;
59
60  // Output file
61  ofstream conv("resconv.d") ;
62  ofstream kss_out("kss.d") ;
63  conv << " # diff_nn diff_psi diff_beta " << endl ;
64
65  // Initial values of the unknowns
66  Scalar nntilde_j = nn() * pow(psi(), a) ;
67  nntilde_j.std_spectral_base() ;
68  Scalar psi_j = psi() ;
69  Vector beta_j = beta() ;
70
71
72  // Iteration
73  //==========
74
75  for (int mer=0; mer<niter; mer++) {
76
77
78  // Useful functions
79  //-----------------
80  Scalar temp_scal_0 (mp) ;
81  Scalar temp_scal_2 (mp) ;
82  Scalar temp_scal_3 (mp) ;
83  Scalar temp_scal_4 (mp) ;
84
85  Vector temp_vect_2 (mp, CON, triad) ;
86  Vector temp_vect_3 (mp, CON, triad) ;
87  Vector temp_vect_4 (mp, CON, triad) ;
88
89  Scalar psi_j_a = pow(psi_j, a) ;
90  psi_j_a.std_spectral_base() ;
91  Scalar psi_j_ma = pow(psi_j, -a) ;
92  psi_j_ma.std_spectral_base() ;
93  Vector beta_j_d = beta_j.down(0, met_gamt) ;
94  Scalar nnpsia_j = nntilde_j * psi_j_a /( psi_j*psi_j*psi_j*psi_j*psi_j*psi_j) ; // to be used in the beta source
95
96
97  //Dynamical relaxation
98  //--------------------
99  // double relax_init = 0.05 ;
100  // double relax_speed = 0.005 ;
101
102  // relax_nn = relax_nn_fin - ( relax_nn_fin - relax_init ) * exp (- relax_speed * mer) ;
103  // relax_psi = relax_psi_fin - ( relax_psi_fin - relax_init ) * exp (- relax_speed * mer) ;
104  // relax_beta = relax_beta_fin - ( relax_beta_fin - relax_init ) * exp (- relax_speed * mer) ;
105  cout << " relax_nn = " << relax_nn << endl ;
106  cout << " relax_psi = " << relax_psi << endl ;
107  cout << " relax_beta = " << relax_beta << endl ;
108
109  //========
110  // Sources
111  //========
112
113  // Source for psi
114  //---------------
115
116  Scalar sour_psi (mp) ;
117
118  // Source physique
119
120  temp_scal_3 = 1./8. * met_gamt.ricci_scal() * psi_j ;
121  temp_scal_3.inc_dzpuis() ;
122  temp_scal_4 = 1./12. * trK*trK * psi_j*psi_j*psi_j*psi_j*psi_j
123  - 1./32. * psi_j*psi_j*psi_j*psi_j*psi_j * psi_j_ma*psi_j_ma / (nntilde_j*nntilde_j) *
124  contract(beta_j_d.ope_killing_conf(met_gamt), 0, 1, beta_j.ope_killing_conf(met_gamt), 0, 1) ;
125
126  sour_psi = temp_scal_3 + temp_scal_4 ;
127
128
129  // Correction Source
130
131  temp_scal_3 = - contract (hh(), 0, 1, psi_j.derive_cov(ff).derive_cov(ff), 0, 1) ;
132  temp_scal_3.inc_dzpuis() ;
133  temp_scal_4 = - contract (hdirac(), 0, psi_j.derive_cov(ff), 0) ;
134
135  sour_psi += temp_scal_3 + temp_scal_4 ;
136
137  sour_psi.annule_domain(0) ;
138
139  // Source for nntilde
140  //--------------------
141
142  Scalar sour_nntilde (mp) ;
143
144  // Source physique
145
146  temp_scal_2 = - psi_j*psi_j*psi_j*psi_j * psi_j_ma * trK_point ;
147  temp_scal_2.inc_dzpuis(2) ;
148
149  temp_scal_3 = psi_j*psi_j*psi_j*psi_j * psi_j_ma * contract (beta_j, 0, trK.derive_cov(ff), 0)
150  - a/8. * nntilde_j * met_gamt.ricci_scal() ;
151  temp_scal_3.inc_dzpuis() ;
152
153  temp_scal_4 = nntilde_j * (4-a)/12. * psi_j*psi_j*psi_j*psi_j * trK*trK
154  - 2 * (a+1) * contract(psi_j.derive_cov(ff), 0, nntilde_j.derive_con(met_gamt), 0) / psi_j
155  - a*(a+1) * nntilde_j * contract(psi_j.derive_cov(met_gamt), 0, psi_j.derive_con(met_gamt), 0) / (psi_j * psi_j)
156  + (a+8)/32. * psi_j*psi_j*psi_j*psi_j * psi_j_ma*psi_j_ma/nntilde_j
157  * contract(beta_j_d.ope_killing_conf(met_gamt), 0, 1, beta_j.ope_killing_conf(met_gamt), 0, 1) ;
158
159  sour_nntilde = temp_scal_2 + temp_scal_3 + temp_scal_4 ;
160
161
162  // Correction Source
163  temp_scal_3 = - contract (hh(), 0, 1, nntilde_j.derive_cov(ff).derive_cov(ff), 0, 1) ;
164  temp_scal_3.inc_dzpuis() ;
165
166  temp_scal_4 = - contract (hdirac(), 0, nntilde_j.derive_cov(ff), 0) ;
167
168  sour_nntilde += temp_scal_3 + temp_scal_4 ;
169  sour_nntilde.annule_domain(0) ;
170
171  // Source for beta
172  //----------------
173
174  Vector sour_beta(mp, CON, triad) ;
175
176  // Source physique
177
178  temp_vect_3 = 4./3. * psi_j_a * nntilde_j * trK.derive_con(met_gamt)
179  - contract(met_gamt.ricci(), 1, beta_j, 0).up(0, met_gamt) ;
180  temp_vect_3.inc_dzpuis() ;
181
182  temp_vect_4 = contract(beta_j.ope_killing_conf(met_gamt), 1, nnpsia_j.derive_cov(ff), 0) / nnpsia_j ;
183
184  sour_beta = temp_vect_3 + temp_vect_4 ;
185
186  // Correction Source
187
188
189  temp_vect_3 = (lambda - 1./3.) * contract (beta_j.derive_cov(ff).derive_con(ff), 0, 1)
190  - contract(contract(met_gamt.connect().get_delta(), 1, beta_j, 0).derive_con(ff), 1, 2)
191  - contract(hh(), 0, 1, beta_j.derive_cov(met_gamt).derive_cov(ff), 1, 2)
192  - 1./3. * contract (hh(), 1, beta_j.divergence(ff).derive_cov(ff), 0) ;
193  temp_vect_3.inc_dzpuis() ;
194
195  temp_vect_4 = - contract(hdirac(), 0, beta_j.derive_cov(met_gamt), 1)
196  - contract(met_gamt.connect().get_delta(), 1, 2, beta_j.derive_con(met_gamt), 0, 1) ;
197
198  sour_beta += temp_vect_3 + temp_vect_4 ;
199
200  //====================
201  // Boundary conditions
202  //====================
203
204
205  // BC Psi
206  //-------
207  Scalar tmp = psi_j * psi_j * psi_j * trK
209  - psi_j*psi_j*psi_j*psi_j_ma/(2. * nntilde_j)
210  * contract( beta_j_d.ope_killing_conf(met_gamt), 0, 1, met_gamt.radial_vect()
211  * met_gamt.radial_vect(), 0, 1)
212  - 1./3. * psi_j*psi_j*psi_j * trK * contract( met_gamt.cov(), 0, 1, met_gamt.radial_vect()
213  * met_gamt.radial_vect(), 0, 1)
214  - 4 * ( met_gamt.radial_vect()(2) * psi_j.derive_cov(ff)(2)
215  + met_gamt.radial_vect()(3) * psi_j.derive_cov(ff)(3) ) ;
216
217  tmp = tmp / (4 * met_gamt.radial_vect()(1)) ;
218
219  // in this case you don't have to substract any value
220  Valeur psi_bound (mp.get_mg()->get_angu() ) ;
221
222  int nnp = mp.get_mg()->get_np(1) ;
223  int nnt = mp.get_mg()->get_nt(1) ;
224
225  psi_bound = 1 ;
226
227  for (int k=0 ; k<nnp ; k++)
228  for (int j=0 ; j<nnt ; j++)
229  psi_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
230
231  psi_bound.std_base_scal() ;
232
233
234
235  // BC beta
236  //-------
237  const Coord& y = mp.y ;
238  const Coord& x = mp.x ;
239
240  Scalar xx(mp) ;
241  xx = x ;
242  xx.std_spectral_base() ;
243
244  Scalar yy(mp) ;
245  yy = y ;
246  yy.std_spectral_base() ;
247
248  Vector tmp_vect = nntilde_j * psi_j_a/(psi_j * psi_j) * met_gamt.radial_vect() ;
249  tmp_vect.change_triad(mp.get_bvect_cart() ) ;
250
251  // Beta-x
252
253  Valeur lim_x (mp.get_mg()->get_angu()) ;
254
255  lim_x = 1 ; // Juste pour affecter dans espace des configs ;
256
257  for (int k=0 ; k<nnp ; k++)
258  for (int j=0 ; j<nnt ; j++)
259  lim_x.set(0, k, j, 0) = omega * yy.val_grid_point(1, k, j, 0)
260  + tmp_vect(1).val_grid_point(1, k, j, 0) ;
261
262  lim_x.set_base(*(mp.get_mg()->std_base_vect_cart()[0])) ;
263
264
265  // Beta-y
266
267  Valeur lim_y (mp.get_mg()->get_angu()) ;
268
269  lim_y = 1 ; // Juste pour affecter dans espace des configs ;
270
271  for (int k=0 ; k<nnp ; k++)
272  for (int j=0 ; j<nnt ; j++)
273  lim_y.set(0, k, j, 0) = - omega * xx.val_grid_point(1, k, j, 0)
274  + tmp_vect(2).val_grid_point(1, k, j, 0) ;
275
276  lim_y.set_base(*(mp.get_mg()->std_base_vect_cart()[1])) ;
277
278  // Beta-z
279
280  Valeur lim_z (mp.get_mg()->get_angu()) ;
281
282  lim_z = 1 ; // Juste pour affecter dans espace des configs ;
283
284  for (int k=0 ; k<nnp ; k++)
285  for (int j=0 ; j<nnt ; j++)
286  lim_z.set(0, k, j, 0) = tmp_vect(3).val_grid_point(1, k, j, 0) ;
287
288  lim_z.set_base(*(mp.get_mg()->std_base_vect_cart()[2])) ;
289
290
291  // START of AEI
292  //============================
293
294
295  // BC AEI
296
297  Vector stilde_j = met_gamt.radial_vect() ;
298  Scalar btilde_j = contract (beta_j_d, 0, stilde_j, 0) ;
299
300
301  // HH_tilde
302  Scalar hh_tilde = contract(stilde_j.derive_cov(met_gamt), 0, 1) ;
303  hh_tilde.dec_dzpuis(2) ;
304
305  // Tangential part of the shift
306  tmp_vect = btilde_j * stilde_j ;
307  Vector VV_j = btilde_j * stilde_j - beta_j ;
308
309  // "Acceleration" term V^a \tilde{D}_a ln M
310  Scalar accel = 2 * contract( VV_j, 0, contract( stilde_j, 0, stilde_j.down(0,
311  met_gamt).derive_cov(met_gamt), 1), 0) ;
312
313
314  // Divergence of V^a
315  Sym_tensor qq_spher = met_gamt.con() - stilde_j * stilde_j ;
316  Scalar div_VV = contract( qq_spher.down(0, met_gamt), 0, 1, VV_j.derive_cov(met_gamt), 0, 1) ;
317
318  Vector angular (mp, CON, mp.get_bvect_cart()) ;
319  Scalar yya (mp) ;
320  yya = mp.ya ;
321  Scalar xxa (mp) ;
322  xxa = mp.xa ;
323
324  angular.set(1) = - omega * yya ;
325  angular.set(2) = omega * xxa ;
326  angular.set(3).annule_hard() ;
327
328
329  angular.set(1).set_spectral_va()
330  .set_base(*(mp.get_mg()->std_base_vect_cart()[0])) ;
331  angular.set(2).set_spectral_va()
332  .set_base(*(mp.get_mg()->std_base_vect_cart()[1])) ;
333  angular.set(3).set_spectral_va()
334  .set_base(*(mp.get_mg()->std_base_vect_cart()[2])) ;
335
336  angular.change_triad(mp.get_bvect_spher()) ;
337
338  //kss from from L_lN=0 condition
339  //------------------------------
340  Scalar corr_nn_kappa (mp) ;
341  corr_nn_kappa = nntilde_j * psi_j_a - psi_j*psi_j*contract (beta_j_d, 0, met_gamt.radial_vect(), 0) ;
342  // corr_nn_kappa.inc_dzpuis(2) ;
343  corr_nn_kappa = sqrt(sqrt (corr_nn_kappa*corr_nn_kappa)) ;
344  corr_nn_kappa.std_spectral_base() ;
345  // corr_nn_kappa.inc_dzpuis(2) ;
346
347  Scalar kss = - kappa_hor() * psi_j_ma / nntilde_j ;
348  kss.inc_dzpuis(2) ;
349  kss += contract(stilde_j, 0, nntilde_j.derive_cov(ff), 0)/ (psi_j*psi_j*nntilde_j)
350  + a * contract(stilde_j, 0, psi_j.derive_cov(ff), 0) / (psi_j*psi_j*psi_j)
351  + 0.1 * contract (nntilde_j.derive_cov(ff), 0, met_gamt.radial_vect(), 0) * corr_nn_kappa ;
352
353
354
355  // beta^r component
356  //-----------------
357  double rho = 5. ; // rho>1 ; rho=1 "pure Dirichlet" version
358  Scalar beta_r_corr = (rho - 1) * btilde_j * hh_tilde;
359  beta_r_corr.inc_dzpuis(2) ;
360  Scalar nn_KK = nntilde_j * psi_j_a * trK ;
361  nn_KK.inc_dzpuis(2) ;
362
363  beta_r_corr.set_dzpuis(2) ;
364  nn_KK.set_dzpuis(2) ;
365  accel.set_dzpuis(2) ;
366  div_VV.set_dzpuis(2) ;
367
368  Scalar b_tilde_new (mp) ;
369  b_tilde_new = 2 * contract(stilde_j, 0, btilde_j.derive_cov(ff), 0)
370  + beta_r_corr
371  - 3 * nntilde_j * psi_j_a * kss + nn_KK + accel + div_VV ;
372
373  b_tilde_new = b_tilde_new / (hh_tilde * rho) ;
374
375  b_tilde_new.dec_dzpuis(2) ;
376
377  tmp_vect.set(1) = b_tilde_new * stilde_j(1) + angular(1) ;
378  tmp_vect.set(2) = b_tilde_new * stilde_j(2) + angular(2) ;
379  tmp_vect.set(3) = b_tilde_new * stilde_j(3) + angular(3) ;
380
381  tmp_vect.change_triad(mp.get_bvect_cart() ) ;
382
383  // Beta-x
384  Valeur lim_x_AEI (mp.get_mg()->get_angu()) ;
385
386  lim_x_AEI = 1 ; // Juste pour affecter dans espace des configs ;
387
388  for (int k=0 ; k<nnp ; k++)
389  for (int j=0 ; j<nnt ; j++)
390  lim_x_AEI.set(0, k, j, 0) = tmp_vect(1).val_grid_point(1, k, j, 0) ;
391
392  lim_x_AEI.set_base(*(mp.get_mg()->std_base_vect_cart()[0])) ;
393
394
395  // Beta-y
396  Valeur lim_y_AEI (mp.get_mg()->get_angu()) ;
397
398  lim_y_AEI = 1 ; // Juste pour affecter dans espace des configs ;
399
400  for (int k=0 ; k<nnp ; k++)
401  for (int j=0 ; j<nnt ; j++)
402  lim_y_AEI.set(0, k, j, 0) = tmp_vect(2).val_grid_point(1, k, j, 0) ;
403
404  lim_y_AEI.set_base(*(mp.get_mg()->std_base_vect_cart()[1])) ;
405
406  // Beta-z
407  Valeur lim_z_AEI (mp.get_mg()->get_angu()) ;
408
409  lim_z_AEI = 1 ; // Juste pour affecter dans espace des configs ;
410
411  for (int k=0 ; k<nnp ; k++)
412  for (int j=0 ; j<nnt ; j++)
413  lim_z_AEI.set(0, k, j, 0) = tmp_vect(3).val_grid_point(1, k, j, 0) ;
414
415  lim_z_AEI.set_base(*(mp.get_mg()->std_base_vect_cart()[2])) ;
416
417
418  // End of AEI
419  //============================
420
421
422
423
424  // BC nn_tilde
425  //------------
426
427  // BC kappa=const
428
429  //Area
430  Scalar darea = psi_j* psi_j* psi_j* psi_j
431  * sqrt( met_gamt.cov()(2,2) * met_gamt.cov()(3,3) -
432  met_gamt.cov()(2,3) * met_gamt.cov()(2,3) ) ;
433
434  darea.std_spectral_base() ;
435
436  Scalar area_int (darea) ;
437  area_int.raccord(1) ;
438  double area = mp.integrale_surface(area_int, radius + 1e-15) ;
439
440  //Radius
441  double radius = area / (4. * M_PI);
442  radius = pow(radius, 1./2.) ;
443
444  // Angular momentum
445  Vector phi (mp, CON, *(ff.get_triad()) ) ;
446  tmp = 1 ;
447  tmp.std_spectral_base() ;
448  tmp.mult_rsint() ;
449  phi.set(1) = 0. ;
450  phi.set(2) = 0. ;
451  phi.set(3) = tmp ;
452
453  Scalar kksphi = psi_j*psi_j*psi_j_ma/(2.*nntilde_j) *
454  contract( beta_j_d.ope_killing_conf(met_gamt), 0, 1, met_gamt.radial_vect()* phi, 0, 1) +
455  1./3. * trK * psi_j*psi_j *
456  contract( met_gamt.cov(), 0, 1, met_gamt.radial_vect()* phi, 0, 1) ;
457
458  kksphi = kksphi / (8. * M_PI) ;
459  Scalar kkspsi_int = kksphi * darea ; // we correct with the curved
460  // element of area
461  double ang_mom = mp.integrale_surface(kkspsi_int, radius + 1e-15) ;
462
463  // Surface gravity
464  double kappa_kerr = (pow( radius, 4) - 4 * pow( ang_mom, 2)) / ( 2 * pow( radius, 3)
465  * sqrt( pow( radius, 4) + 4 * pow( ang_mom, 2) ) ) ;
466
467
468  // BC condition itself
469  rho = 5. ;
470
471  Scalar kappa (mp) ;
472  if (mer < 0)
473  kappa = kappa_kerr ;
474  else
475  kappa = 0.15 ;
476  kappa.std_spectral_base() ;
477  kappa.inc_dzpuis(2) ;
478
479
480
481  tmp = - a / psi_j * nntilde_j
482  * contract(met_gamt.radial_vect(), 0, psi_j.derive_cov(met_gamt), 0)
483  + 1./2. * psi_j * psi_j * psi_j_ma
484  * contract( beta_j_d.ope_killing_conf(met_gamt), 0, 1, met_gamt.radial_vect()
485  * met_gamt.radial_vect(), 0, 1)
486  + 1./3. * nntilde_j * psi_j * psi_j * trK
487  * contract( met_gamt.cov(), 0, 1, met_gamt.radial_vect()* met_gamt.radial_vect(), 0, 1)
488  + psi_j * psi_j * psi_j_ma * kappa
489  - met_gamt.radial_vect()(2) * nntilde_j.derive_cov(ff)(2)
490  - met_gamt.radial_vect()(3) * nntilde_j.derive_cov(ff)(3)
491  + ( rho - 1 ) * met_gamt.radial_vect()(1) * nntilde_j.derive_cov(ff)(1)
492  + 0. * contract (nntilde_j.derive_cov(ff), 0, met_gamt.radial_vect(), 0) * corr_nn_kappa ;
493
494  tmp = tmp / (rho * met_gamt.radial_vect()(1)) ;
495
496  // in this case you don't have to substract any value
497  Valeur nn_bound_kappa (mp.get_mg()->get_angu() ) ;
498
499  nn_bound_kappa = 1 ;
500
501  for (int k=0 ; k<nnp ; k++)
502  for (int j=0 ; j<nnt ; j++)
503  nn_bound_kappa.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
504
505  nn_bound_kappa.std_base_scal() ;
506
507
508  /*
509  // BC kappa const Dir
510
511
512  rho = 10. ;
513  tmp = 1./(a * contract(met_gamt.radial_vect(), 0, psi_j.derive_cov(met_gamt), 0)) *
514  (psi_j* psi_j*psi_j * psi_j_ma* kappa
515  - psi_j * contract(met_gamt.radial_vect(), 0, nntilde_j.derive_cov(met_gamt), 0)
516  + psi_j* psi_j*psi_j * psi_j_ma/2. *
517  contract(beta_j_d.ope_killing_conf(met_gamt), 0, 1, met_gamt.radial_vect()
518  * met_gamt.radial_vect(), 0, 1)
519  + 1./3. * nntilde_j * psi_j * psi_j *psi_j * trK
520  * contract( met_gamt.cov(), 0, 1, met_gamt.radial_vect()* met_gamt.radial_vect(), 0, 1)) ;
521
522  tmp = (tmp + rho * nn())/(1 + rho) ;
523  tmp = tmp - 1 ;
524
525  // in this case you don't have to substract any value
526  Valeur nn_bound_kappa_Dir (mp.get_mg()->get_angu() ) ;
527
528  nn_bound_kappa_Dir = 1 ;
529
530  for (int k=0 ; k<nnp ; k++)
531  for (int j=0 ; j<nnt ; j++)
532  nn_bound_kappa_Dir.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
533
534  nn_bound_kappa_Dir.std_base_scal() ;
535  */
536
537  // BC nn = const
538
539  rho = 0. ;
540  tmp = 0.2 * psi_j_ma ;
541  tmp = (tmp + rho * nntilde_j)/(1 + rho) ;
542  tmp = tmp - 1 ;
543
544  // in this case you don't have to substract any value
545  Valeur nn_bound_const (mp.get_mg()->get_angu() ) ;
546
547  nn_bound_const = 1 ;
548
549  for (int k=0 ; k<nnp ; k++)
550  for (int j=0 ; j<nnt ; j++)
551  nn_bound_const.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
552
553  nn_bound_const.std_base_scal() ;
554
555  boundary_nn_Neu_kk(0) ;
556
557  // BC nn AEI
558
559  rho = 5. ;
560  tmp = btilde_j * psi_j*psi_j*psi_j_ma ;
561  tmp = (tmp + rho * nntilde_j)/(1 + rho) ;
562  tmp = tmp - 1 ;
563
564  // in this case you don't have to substract any value
565  Valeur nn_bound_AEI (mp.get_mg()->get_angu() ) ;
566
567  nn_bound_AEI = 1 ;
568
569  for (int k=0 ; k<nnp ; k++)
570  for (int j=0 ; j<nnt ; j++)
571  nn_bound_AEI.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
572
573  nn_bound_AEI.std_base_scal() ;
574
575
576
577  //============================
578  // Resolution of the equations
579  //============================
580
581  Scalar psi_jp1(mp) ;
582  Scalar nntilde_jp1(mp) ;
583  Vector beta_jp1(beta_j) ;
584
585
586  if (solve_psi == 1) {
587  psi_jp1 = sour_psi.poisson_neumann(psi_bound, 0) + 1. ;
588
589  // Test:
590  maxabs(psi_jp1.laplacian() - sour_psi,
591  "Absolute error in the resolution of the equation for Psi") ;
592  // Relaxation (relax=1 -> new ; relax=0 -> old )
593  psi_jp1 = relax_psi * psi_jp1 + (1 - relax_psi) * psi_j ;
594  }
595
596
597  if (solve_lapse == 1) {
598  // nntilde_jp1 = sour_nntilde.poisson_neumann(nn_bound_kappa, 0) + 1. ;
599  // nntilde_jp1 = sour_nntilde.poisson_dirichlet(nn_bound_const, 0) + 1. ;
600  // nntilde_jp1 = sour_nntilde.poisson_dirichlet(nn_bound_kappa_Dir, 0) + 1. ;
601  nntilde_jp1 = sour_nntilde.poisson_dirichlet(nn_bound_AEI, 0) + 1. ;
602
603  // Test:
604  maxabs(nntilde_jp1.laplacian() - sour_nntilde,
605  "Absolute error in the resolution of the equation for nntilde") ;
606  // Relaxation (relax=1 -> new ; relax=0 -> old )
607  nntilde_jp1 = relax_nn * nntilde_jp1 + (1 - relax_nn) * nntilde_j ;
608  }
609
610  if (solve_shift == 1) {
611  double precision = 1e-8 ;
612  // poisson_vect_boundary(lambda, sour_beta, beta_jp1, lim_x,
613  // lim_y, lim_z, 0, precision, 20) ;
614  poisson_vect_boundary(lambda, sour_beta, beta_jp1, lim_x_AEI,
615  lim_y_AEI, lim_z_AEI, 0, precision, 20) ;
616
617
618  // Test
619  sour_beta.dec_dzpuis() ;
620  maxabs(beta_jp1.derive_con(ff).divergence(ff)
621  + lambda * beta_jp1.divergence(ff)
622  .derive_con(ff) - sour_beta,
623  "Absolute error in the resolution of the equation for beta") ;
624
625  cout << endl ;
626
627  // Relaxation (relax=1 -> new ; relax=0 -> old )
628  beta_jp1 = relax_beta * beta_jp1 + (1 - relax_beta) * beta_j ;
629  }
630
631
632
633  //====================
634  // Convergence control
635  //====================
636
637  double diff_nn, diff_psi, diff_beta ;
638  diff_nn = 1.e-16 ;
639  diff_psi = 1.e-16 ;
640  diff_beta = 1.e-16 ;
641  if (solve_lapse == 1)
642  diff_nn = max( diffrel(nntilde_j, nntilde_jp1) ) ;
643  if (solve_psi == 1)
644  diff_psi = max( diffrel(psi_j, psi_jp1) ) ;
645  if (solve_shift == 1)
646  diff_beta = max( maxabs(beta_jp1 - beta_j) ) ;
647
648  cout << "step = " << mer << " : diff_psi = " << diff_psi
649  << " diff_nn = " << diff_nn
650  << " diff_beta = " << diff_beta << endl ;
651  cout << "----------------------------------------------" << endl ;
652  if ((diff_psi<precis) && (diff_nn<precis) && (diff_beta<precis))
653  break ;
654
655  if (mer>0) {conv << mer << " " << log10(diff_nn) << " " << log10(diff_psi)
656  << " " << log10(diff_beta) << endl ; } ;
657
658
659
660  //======================
661  // Updates for next step
662  //======================
663
664  psi_j = psi_jp1 ;
665  nntilde_j = nntilde_jp1 ;
666  beta_j = beta_jp1 ;
667
668
669  //========================
670  // Savings of output files
671  //========================
672
673
674  Scalar kkss ( psi_j_ma/(2. * nntilde_j )
675  * contract( beta_j_d.ope_killing_conf(met_gamt), 0, 1, met_gamt.radial_vect()
676  * met_gamt.radial_vect(), 0, 1)
677  + 1./3. * trK * contract( met_gamt.cov(), 0, 1, met_gamt.radial_vect()
678  * met_gamt.radial_vect(), 0, 1) ) ;
679
680  double max_kss = kkss.val_grid_point(1, 0, 0, 0) ;
681  double min_kss = kkss.val_grid_point(1, 0, 0, 0) ;
682
683  hh_tilde = contract(met_gamt.radial_vect().derive_cov(met_gamt), 0, 1) ;
684  double max_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
685  double min_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
686
687
688  for (int k=0 ; k<nnp ; k++)
689  for (int j=0 ; j<nnt ; j++){
690  if (kkss.val_grid_point(1, k, j, 0) > max_kss)
691  max_kss = kkss.val_grid_point(1, k, j, 0) ;
692  if (kkss.val_grid_point(1, k, j, 0) < min_kss)
693  min_kss = kkss.val_grid_point(1, k, j, 0) ;
694
695  if (hh_tilde.val_grid_point(1, k, j, 0) > max_hh_tilde)
696  max_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
697  if (hh_tilde.val_grid_point(1, k, j, 0) < min_hh_tilde)
698  min_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
699  }
700  kss_out << mer << " " << max_kss << " " << min_kss
701  << " " << -1 * max_hh_tilde << " " << -1 * min_hh_tilde << endl ;
702
703  //Intermediate updates for tests
704  if (solve_psi == 1)
705  set_psi(psi_j) ;
706  if (solve_lapse == 1)
707  n_evol.update(nntilde_j * psi_j_a, jtime, ttime) ;
708  if (solve_shift == 1)
709  beta_evol.update(beta_j, jtime, ttime) ;
710
711  if (solve_shift == 1)
712  update_aa() ;
713
714  }
715
716  // Global updates
717  //--------------
718  Scalar psi_j_a = pow(psi_j, a) ;
719  psi_j_a.std_spectral_base() ;
720
721
722  if (solve_psi == 1)
723  set_psi(psi_j) ;
724  if (solve_lapse == 1)
725  n_evol.update(nntilde_j * psi_j_a, jtime, ttime) ;
726  if (solve_shift == 1)
727  beta_evol.update(beta_j, jtime, ttime) ;
728
729  if (solve_shift == 1)
730  update_aa() ;
731
732 /*
733  Vector beta_j_d = beta().down(0, met_gamt) ;
734  Scalar check ( -psi() * psi() * psi() * trK
735  + psi() * met_gamt.radial_vect().divergence(met_gamt)
736  + psi() * psi() * psi() * pow(psi(), -a)/(2. * nntilde_j)
737  * contract( beta_j_d.ope_killing_conf(met_gamt), 0, 1, met_gamt.radial_vect()
738  * met_gamt.radial_vect(), 0, 1)
739  + 1./3. * psi() * psi() * psi() * trK
740  * contract( met_gamt.cov(), 0, 1, met_gamt.radial_vect()
741  * met_gamt.radial_vect(), 0, 1)
742  + 4 * contract (met_gamt.radial_vect(), 0, psi().derive_cov(met_gamt),0) ) ;
743
744
745  des_meridian(check, 1, 4., "CHECK", 1) ;
746  arrete() ;
747
748  Scalar exp = contract(gam().radial_vect().derive_cov(gam()), 0,1)
749  + contract(contract(k_dd(), 0, gam().radial_vect(), 0),
750  0, gam().radial_vect(), 0)
751  - trk() ;
752
753  des_meridian(exp, 1, 4., "Expansion", 1) ;
754  arrete() ;
755 */
756  conv.close() ;
757  kss_out.close() ;
758
759
760
761
762 }
763
764
765 }
Coord xa
Absolute x coordinate.
Definition: map.h:748
virtual const Vector & beta() const
shift vector at the current time step (jtime )
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
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.
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:293
double omega
Angular velocity in LORENE&#39;s units.
Definition: isol_hor.h:269
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
virtual const Vector & hdirac() const
Vector which vanishes in Dirac gauge.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:801
Lorene prototypes.
Definition: app_hor.h:67
double radius
Radius of the horizon in LORENE&#39;s units.
Definition: isol_hor.h:266
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
void update_aa()
Conformal representation of the traceless part of the extrinsic curvature: .
Definition: isol_hor.C:845
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the metric is defined.
Definition: metric.h:309
double kappa_hor() const
Surface gravity.
Definition: phys_param.C:221
const Tensor_sym & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
int jtime
Time step index of the latest slice.
Definition: time_slice.h:193
virtual const Vector & radial_vect() const
Returns the radial vector normal to a spherical slicing and pointing toward spatial infinity...
Definition: metric.C:365
const Sym_tensor & ricci() const
Returns the Ricci tensor (given by the Connection p_connect )
Definition: metric.C:341
Scalar trK
Trace of the extrinsic curvature.
Definition: isol_hor.h:332
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
void set(const Map *mp, Mtbl *(*construct)(const Map *))
Semi-constructor from a mapping and a method.
Definition: coord.C:137
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Definition: scalar_deriv.C:402
const Tensor_sym & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
void inc_dzpuis()
dzpuis += 1 ;
Definition: tenseur.C:1133
const Tensor_sym & get_delta() const
Returns the tensor which defines the connection with respect to the flat one: is the difference bet...
Definition: connection.h:271
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:507
virtual const Sym_tensor & hh(Param *=0x0, Param *=0x0) const
Deviation of the conformal metric from the flat metric : .
Map_af & mp
Affine mapping.
Definition: isol_hor.h:260
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:840
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Definition: sym_tensor.C:352
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
virtual const Scalar & ricci_scal() const
Returns the Ricci scalar.
Definition: metric.C:353
void set_psi(const Scalar &psi_in)
Sets the conformal factor relating the physical metric to the conformal one: .
Definition: isol_hor.C:518
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition: mg3d.C:604
virtual const Connection & connect() const
Returns the connection.
Definition: metric.C:304
Metric met_gamt
3 metric tilde
Definition: isol_hor.h:326
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:438
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:387
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition: tensor.C:1011
Evolution_std< double > the_time
Time label of each slice.
Definition: time_slice.h:196
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition: cmp.C:647
Scalar trK_point
Time derivative of the trace of the extrinsic curvature.
Definition: isol_hor.h:335
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:283
Coord ya
Absolute y coordinate.
Definition: map.h:749
const Valeur boundary_nn_Neu_kk(int nn=1) const
Neumann boundary condition for N using the extrinsic curvature.
Definition: bound_hor.C:423
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
Definition: time_slice.h:510
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:809
void dec_dzpuis()
dzpuis -= 1 ;
Definition: tenseur.C:1120
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition: cmp_math.C:325
Coord y
y coordinate centered on the grid
Definition: map.h:745
Coord x
x coordinate centered on the grid
Definition: map.h:744
Evolution_std< Scalar > n_evol
Values at successive time steps of the lapse function N.
Definition: time_slice.h:219
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
Evolution_std< Vector > beta_evol
Values at successive time steps of the shift vector .
Definition: time_slice.h:222
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )