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
31char 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
43namespace Lorene {
44void 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::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::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::omega
double omega
Angular velocity in LORENE's units.
Definition isol_hor.h:269
Lorene::Isol_hor::trK
Scalar trK
Trace of the extrinsic curvature.
Definition isol_hor.h:332
Lorene::Isol_hor::met_gamt
Metric met_gamt
3 metric tilde
Definition isol_hor.h:326
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::trK_point
Scalar trK_point
Time derivative of the trace of the extrinsic curvature.
Definition isol_hor.h:335
Lorene::Isol_hor::kappa_hor
double kappa_hor() const
Surface gravity
Definition phys_param.C:218
Lorene::Isol_hor::update_aa
void update_aa()
Conformal representation of the traceless part of the extrinsic curvature: .
Definition isol_hor.C:842
Lorene::Isol_hor::radius
double radius
Radius of the horizon in LORENE's units.
Definition isol_hor.h:266
Lorene::Isol_hor::mp
Map_af & mp
Affine mapping.
Definition isol_hor.h:260
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
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::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::Map::y
Coord y
y coordinate centered on the grid
Definition map.h:727
Lorene::Map::ya
Coord ya
Absolute y coordinate.
Definition map.h:731
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::x
Coord x
x coordinate centered on the grid
Definition map.h:726
Lorene::Map::xa
Coord xa
Absolute x coordinate.
Definition map.h:730
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
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::Metric::con
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition metric.C:290
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::Metric::cov
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition metric.C:280
Lorene::Metric::ricci
const Sym_tensor & ricci() const
Returns the Ricci tensor (given by the Connection p_connect )
Definition metric.C:338
Lorene::Metric::connect
virtual const Connection & connect() const
Returns the connection.
Definition metric.C:301
Lorene::Metric::ricci_scal
virtual const Scalar & ricci_scal() const
Returns the Ricci scalar.
Definition metric.C:350
Lorene::Mg3d::get_angu
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition mg3d.C:473
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::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::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::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::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::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::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::Tenseur::inc_dzpuis
void inc_dzpuis()
dzpuis += 1 ;
Definition tenseur.C:1117
Lorene::Tenseur::dec_dzpuis
void dec_dzpuis()
dzpuis -= 1 ;
Definition tenseur.C:1104
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::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::Time_slice_conf::nn
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )
Definition time_slice_conf.C:591
Lorene::Time_slice_conf::hdirac
virtual const Vector & hdirac() const
Vector which vanishes in Dirac gauge.
Definition time_slice_conf.C:815
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::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::Time_slice::jtime
int jtime
Time step index of the latest slice.
Definition time_slice.h:190
Lorene::Time_slice::beta
virtual const Vector & beta() const
shift vector at the current time step (jtime )
Definition time_slice_access.C:87
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::Time_slice::the_time
Evolution_std< double > the_time
Time label of each slice.
Definition time_slice.h:193
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::Vector::divergence
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition vector.C:381
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::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Lorene::log10
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition cmp_math.C:322
Lorene::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:504
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
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::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::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::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition tenseur_operateur.C:279
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
Lorene prototypes.
Definition app_hor.h:64
Unites
Standard units of space, time and mass.