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