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
208  - psi_j * met_gamt.radial_vect().divergence(met_gamt)
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 }
Lorene::Metric::ricci_scal
virtual const Scalar & ricci_scal() const
Returns the Ricci scalar.
Definition: metric.C:350
Lorene::Time_slice_conf::hdirac
virtual const Vector & hdirac() const
Vector which vanishes in Dirac gauge.
Definition: time_slice_conf.C:815
Lorene::Tensor_sym::derive_cov
const Tensor_sym & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition: tensor_sym_calculus.C:192
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::radius
double radius
Radius of the horizon in LORENE's units.
Definition: isol_hor.h:266
Lorene::Scalar::inc_dzpuis
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...
Definition: scalar_r_manip.C:470
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::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::Tenseur::inc_dzpuis
void inc_dzpuis()
dzpuis += 1 ;
Definition: tenseur.C:1117
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::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::Mg3d::get_angu
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition: mg3d.C:473
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::Time_slice_conf::hh
virtual const Sym_tensor & hh(Param *=0x0, Param *=0x0) const
Deviation of the conformal metric from the flat metric : .
Definition: time_slice_conf.C:758
Lorene::Scalar::derive_cov
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
Lorene::Metric::cov
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:280
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::Isol_hor::mp
Map_af & mp
Affine mapping.
Definition: isol_hor.h:260
Lorene::Isol_hor::trK
Scalar trK
Trace of the extrinsic curvature.
Definition: isol_hor.h:332
Lorene::Map_af::integrale_surface
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
Definition: map_af_integ_surf.C:90
Lorene::Isol_hor::omega
double omega
Angular velocity in LORENE's units.
Definition: isol_hor.h:269
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Lorene::Vector::change_triad
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition: vector_change_triad.C:75
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::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::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::Metric::con
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:290
Lorene::Isol_hor::kappa_hor
double kappa_hor() const
Surface gravity
Definition: phys_param.C:218
Lorene::Tenseur::derive_con
const Tenseur & derive_con(const Metrique &) const
Returns the contravariant derivative of *this , with respect to met .
Definition: tenseur.C:1568
Lorene::Coord::set
void set(const Map *mp, Mtbl *(*construct)(const Map *))
Semi-constructor from a mapping and a method.
Definition: coord.C:134
Lorene::Isol_hor::trK_point
Scalar trK_point
Time derivative of the trace of the extrinsic curvature.
Definition: isol_hor.h:335
Lorene::Connection::get_delta
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
Lorene::log10
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition: cmp_math.C:322
Lorene::Time_slice::the_time
Evolution_std< double > the_time
Time label of each slice.
Definition: time_slice.h:193
Lorene::Tensor_sym::derive_con
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...
Definition: tensor_sym_calculus.C:204
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
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::Mg3d::std_base_vect_cart
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
Definition: mg3d_std_base.C:177
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Lorene::Map::x
Coord x
x coordinate centered on the grid
Definition: map.h:726
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::Map::xa
Coord xa
Absolute x coordinate.
Definition: map.h:730
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::Time_slice::jtime
int jtime
Time step index of the latest slice.
Definition: time_slice.h:190
Lorene::Metric::ricci
const Sym_tensor & ricci() const
Returns the Ricci tensor (given by the Connection p_connect )
Definition: metric.C:338
Lorene::Time_slice::beta
virtual const Vector & beta() const
shift vector at the current time step (jtime )
Definition: time_slice_access.C:87
Lorene::Metric_flat::get_triad
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the metric is defined.
Definition: metric.h:309
Lorene::Map::y
Coord y
y coordinate centered on the grid
Definition: map.h:727
Lorene::Sym_tensor::divergence
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Definition: sym_tensor.C:349
Lorene::Metric::connect
virtual const Connection & connect() const
Returns the connection.
Definition: metric.C:301
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition: tenseur_operateur.C:279
Lorene::Tenseur::dec_dzpuis
void dec_dzpuis()
dzpuis -= 1 ;
Definition: tenseur.C:1104
Lorene::Map::ya
Coord ya
Absolute y coordinate.
Definition: map.h:731