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