LORENE
Isol_hole/isol_hole_compute_metric.C
1/*
2 * Method of class Isol_Hole to compute metric data associated to a quasistationary single black hole spacetime slice.
3 *
4 * (see file isol_hole.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2009 Nicolas Vasset
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License as published by
15 * the Free Software Foundation; either version 2 of the License, or
16 * (at your option) any later version.
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// Headers C
30#include "math.h"
31
32// Headers Lorene
33#include "tenseur.h"
34#include "star.h"
35#include "param.h"
36#include "graphique.h"
37#include "unites.h"
38#include "spheroid.h"
39#include "excision_surf.h"
40#include "excision_hor.h"
41#include "param_elliptic.h"
42#include "vector.h"
43#include "scalar.h"
44#include "diff.h"
45#include "proto.h"
46#include "isol_hole.h"
47
48
49namespace Lorene {
50void Isol_hole::compute_stat_metric(double precis, double relax, int mer_max, int mer_max2, bool isvoid){
51
52 // Fundamental constants and units
53 // -------------------------------
54 using namespace Unites ;
55
56 // Noticing the user that the iteration has started
57
58 cout << "================================================================" << endl;
59 cout << "STARTING THE MAIN ITERATION FOR COMPUTING METRIC FIELDS" << endl;
60 cout << " iteration parameters are the following: " << endl;
61 cout << " convergence precision required:" << precis << endl;
62 cout << " max number of global steps :" << mer_max << endl;
63 cout << " relaxation parameter :" << relax << endl;
64 cout << "================================================================" << endl;
65
66
67 // Construction of a multi-grid (Mg3d) and an affine mapping from the class mapping
68 // --------------------------------------------------------------------------------
69 const Map_af* map = dynamic_cast<const Map_af*>(&mp) ;
70 const Mg3d* mgrid = (*map).get_mg();
71
72 // Construct angular grid for h(theta,phi)
73 const Mg3d* g_angu = (*mgrid).get_angu_1dom() ;
74
75 const int nz = (*mgrid).get_nzone(); // Number of domains
76 int nt = (*mgrid).get_nt(1); // Number of collocation points in theta in each domain
77 const int np = (*mgrid).get_np(1);
78 const Coord& rr = (*map).r;
79 Scalar rrr (*map) ;
80 rrr = rr ;
81 rrr.std_spectral_base();
82 assert((rrr.val_grid_point(1,0,0,0) - 1.) <= 1.e-9); // For now the code handles only horizons at r=1, corresponding to the first shell inner boundary. This test assures this is the case with our mapping.
83
84 // Angular mapping defined as well
85 //--------------------------------
86 double r_limits2[] = {rrr.val_grid_point(1,0,0,0), rrr.val_grid_point(2,0,0,0)} ;
87 const Map_af map_2(*g_angu, r_limits2); //2D mapping; check if this is useful.
88 const Metric_flat& mets = (*map).flat_met_spher() ;
89
90 //----------------
91 // Initializations
92 // ---------------
93
94
95
96 Scalar logn (*map) ; logn = log(lapse) ;
97 logn.std_spectral_base();
98
99
100 Scalar logpsi(*map) ; logpsi = log(conf_fact) ;
101 logpsi.std_spectral_base();
102 Scalar psi4 (*map) ; psi4 = conf_fact*conf_fact*conf_fact*conf_fact ;
103 Scalar npsi (*map) ; npsi =conf_fact*lapse ;
104
105
106 Scalar conf_fact_new(*map) ; conf_fact_new.annule_hard(); conf_fact_new.std_spectral_base();
107 Scalar npsi_new(*map); npsi_new.annule_hard(); npsi_new.std_spectral_base();
108
109 Vector shift_new (*map, CON, (*map).get_bvect_spher());
110 for(int i=1; i<=3; i++){
111 shift_new.set(i)=0;
112 }
113 shift_new.std_spectral_base();
114
115 // Non conformally flat variables
116
117
118 Sym_tensor gamtuu = mets.con() + hij;
119 Metric gamt(gamtuu);
120 Metric gam(gamt.cov()*psi4) ;
121 Sym_tensor gamma = gam.cov();
122
123
124 // Extrinsic curvature variables
125
126 Sym_tensor aa(*map, CON, (*map).get_bvect_spher());
127 for (int iii= 1; iii<=3; iii++){
128 for(int j=1; j<=3; j++){
129 aa.set(iii,j)= 0;
130 }
131 }
132 aa.std_spectral_base();
133 Scalar aa_quad_scal(*map) ; aa_quad_scal = 0. ;
134
135 Sym_tensor aa_hat(*map, CON, (*map).get_bvect_spher());
136 for (int iii= 1; iii<=3; iii++){
137 for(int j=1; j<=3; j++){
138 aa_hat.set(iii,j)= 0;
139 }
140 }
141
142 Sym_tensor kuu = aa/psi4 ;
143 Sym_tensor kuu2 = aa_hat/(psi4*psi4*sqrt(psi4));
144 Sym_tensor kdd = contract (gamma, 0, contract(gamma, 1, kuu, 0),1);
145
146
147 // (2,1)-rank delta tensor: difference between ricci rotation coefficients.
148
149 Tensor delta = -0.5*contract( hij, 1, gamt.cov().derive_cov(mets), 2);
150 Scalar tmp(*map);
151
152 for (int i=1; i<=3; i++) {
153 for (int j=1; j<=3; j++) {
154 for (int k=1; k<=3; k++) {
155 tmp = 0.;
156 tmp = -0.5 *(gamt.cov().derive_con(mets))(i,j,k);
157 for (int l=1; l<=3; l++) {
158 tmp += -0.5*( gamt.cov()(i,l)*(hij.derive_cov(mets))(k,l,j) + gamt.cov()(l,j)*(hij.derive_cov(mets))(k,l,i));
159 }
160 delta.set(k,i,j) += tmp ;
161 }
162 }
163 }
164
165
166 // Conformal Rstar scalar(eq 61, Bonazzola et al. 2003)
167 Scalar Rstar =
168 0.25 * contract(gamt.con(), 0, 1,
169 contract(hij.derive_cov(mets), 0, 1, gamt.cov().derive_cov(mets), 0, 1), 0, 1 )
170 - 0.5 * contract(gamt.con(), 0, 1,
171 contract(hij.derive_cov(mets), 0, 1, gamt.cov().derive_cov(mets), 0, 2), 0, 1 ) ;
172
173 Scalar norm(*map);
174 Scalar norm3(*map);
175
176
177 if (isvoid == false){
178
179 cout <<"FAIL: case of non-void spacetime not treated yet" << endl;
180 }
181
182 else
183 {
184
185 // Parameters for the iteration
186
187
188 double diff_ent = 1 ; // initialization of difference marker between two iterations;
189
190 int util = 0; // Tool used to stop tensorial iteration at any wished step "util"
191
195
196
197 for(int mer=0 ;(diff_ent > precis) && (mer<mer_max) ; mer++) {
198
199 //Global relaxation coefficient
200
201 // Scalar variables linked to the norm of normal vector to horizon.
202 norm = sqrt(1. + hij(1,1)); norm.std_spectral_base();
203 norm3 = sqrt(1. + hij(3,3)); norm3.std_spectral_base();
204
205
206
210
211
212
214 // Solving for (Psi-1) //
216
217
218 // Setting of the boundary
219
220 double diric= 0.;
221 double neum = 1.;
222
223 Vector ssalt = rrr.derive_cov(gam);
224 Vector ssaltcon = ssalt.up_down(gam);
225 Scalar ssnormalt = sqrt(contract (ssalt,0, ssaltcon, 0));
226 ssnormalt.std_spectral_base();
227
228 ssalt.annule_domain(nz-1);
229 ssalt.annule_domain(0);
230 ssaltcon.annule_domain(nz-1);
231 ssaltcon.annule_domain(0);
232
233 ssalt = ssalt/ssnormalt;
234 ssaltcon = ssaltcon/ssnormalt;
235 Vector ssconalt = ssaltcon*conf_fact*conf_fact; // \tilde{s} in the notations of Gourgoulhon and Jaramillo, 2006
236 ssconalt.std_spectral_base();
237 ssconalt.annule_domain(nz-1);
238 Scalar bound3bis = -((1./conf_fact)*contract((contract(kdd,1,ssconalt,0)),0, ssconalt,0));
239
240 bound3bis.annule_domain(nz-1);
241 bound3bis += -conf_fact*ssconalt.divergence(gamt);
242 bound3bis.annule_domain(nz-1);
243 bound3bis = 0.25*bound3bis;
244 bound3bis += -contract(conf_fact.derive_cov(gamt),0,ssconalt,0) + conf_fact.dsdr();
245 bound3bis.annule_domain(nz-1);
246 bound3bis.std_spectral_base();
247 bound3bis.set_spectral_va().ylm();
248
249 Mtbl_cf *boundd3bis = bound3bis.set_spectral_va().c_cf;
250
251
253
254
255
256 // Computing the source
257
258 Scalar source_conf_fact(*map) ; source_conf_fact=3. ; // Pour le fun...
259 source_conf_fact.std_spectral_base();
260
261 Scalar d2logpsi = contract(conf_fact.derive_cov(mets).derive_cov(mets), 0, 1, hij, 0,1);
262 d2logpsi.inc_dzpuis(1);
263
264 source_conf_fact = -(0.125* aa_quad_scal )/(psi4*conf_fact*conf_fact*conf_fact) + conf_fact* 0.125* Rstar - d2logpsi;
265
266 source_conf_fact.std_spectral_base();
267
268 if (source_conf_fact.get_etat() == ETATZERO) {
269 source_conf_fact.annule_hard() ;
270 source_conf_fact.set_dzpuis(4) ;
271 source_conf_fact.std_spectral_base() ;
272 }
273 source_conf_fact.set_spectral_va().ylm();
274
275
276
278
279 // System inversion
280
281 Param_elliptic source11(source_conf_fact);
282 conf_fact_new = source_conf_fact.sol_elliptic_boundary(source11, *boundd3bis, diric , neum) + 1 ; // Resolution has been done for quantity Q-1, because our solver gives a vanishing solution at infinity!
283
284
286
287 // tests for resolution
288
289 Scalar baba2 = (conf_fact_new-1).laplacian();
290// cout << "psi+1-résolution" << endl;
291// maxabs (baba2 - source_conf_fact);
292
293 Scalar psinewbis = conf_fact_new -1. ; psinewbis.annule_domain(nz -1);
294 psinewbis.std_spectral_base();
295 psinewbis = psinewbis.dsdr();
296 Scalar psinewfin2 (map_2) ;
297 psinewfin2.allocate_all();
298 psinewfin2.set_etat_qcq();
299 psinewfin2.std_spectral_base();
300
301 for (int k=0; k<np; k++)
302 for (int j=0; j<nt; j++) {
303
304 psinewfin2.set_grid_point(0, k, j, 0) = psinewbis.val_grid_point(1, k,j,0) - bound3bis.val_grid_point(1, k, j, 0);
305
306 }
307
308 // maxabs (psinewfin2);
309
310
312
313 // Update during the loop
314
315 conf_fact = conf_fact_new* (1-relax) + conf_fact* relax ;
316 psi4 = conf_fact*conf_fact*conf_fact*conf_fact;
317 logpsi = log(conf_fact) ;
318 logpsi.std_spectral_base();
319
320
321
322
326
327
328
330 // Solving for (N*Psi -1)/
332
333
334 // Setting of the boundary
335 assert (NorKappa == false) ;
336 Scalar bound(*map);
337 bound = (boundNoK)*conf_fact -1;
338 bound.annule_domain(nz -1);
339 bound.std_spectral_base();
340 bound.set_spectral_va().ylm();
341 Mtbl_cf *boundd = bound.get_spectral_va().c_cf;
342
343 diric =1;
344 neum = 0 ;
345
346
347
349
350 // Computing the source ...
351 Scalar d2lognpsi = contract(npsi.derive_cov(mets).derive_cov(mets), 0, 1, hij, 0,1);
352 d2lognpsi.inc_dzpuis(1); // dzpuis correction.
353
354 Scalar source_npsi = npsi*(aa_quad_scal*(7./8.)/(psi4*psi4) + Rstar/8.) - d2lognpsi; // ?
355 source_npsi.std_spectral_base();
356 if (source_npsi.get_etat() == ETATZERO) {
357 source_npsi.annule_hard() ;
358 source_npsi.set_dzpuis(4) ;
359 source_npsi.std_spectral_base() ;
360 }
361
362
364
365 // Inversion of the operator
366 Param_elliptic source1 (source_npsi);
367 npsi_new = source_npsi.sol_elliptic_boundary(source1, *boundd, diric, neum) ;
368
369 npsi_new = npsi_new +1;
370
371
373
374 // Resolution tests in npsi
375 Scalar baba = npsi_new.laplacian();
376// cout << "resolution_npsi" << endl;
377// maxabs (baba - source_npsi);
378
379
380 // cout << "bound_npsi" << endl;
381 Scalar npsibound2 (map_2) ;
382 npsibound2.allocate_all();
383 npsibound2.set_etat_qcq();
384 npsibound2.std_spectral_base();
385 for (int k=0; k<np; k++)
386 for (int j=0; j<nt; j++) {
387
388 npsibound2.set_grid_point(0, k, j, 0) = npsi_new.val_grid_point(1, k,j,0) - bound.val_grid_point(1, k, j, 0) -1.;
389
390 }
391
392 // maxabs (npsibound2);
393
394
396
397 // Update during the loop
398
399
400 npsi = npsi_new*(1-relax) + npsi* relax;
401 lapse = npsi/conf_fact;
402 logn = log(lapse);
403 logn.std_spectral_base();
404
405
409
411 //Resolution in Beta //
413
414
415 // Setting of the boundary
416
417 bound = (boundNoK)/(conf_fact*conf_fact) ;
418 bound.annule_domain(nz -1);
419
420
421 Scalar hor_rot(*map); hor_rot.annule_hard(); hor_rot = Omega; // Rotation parameter for spacetime
422 hor_rot.std_spectral_base() ; hor_rot.mult_rsint();
423 hor_rot.annule_domain(nz -1);
424
425 Vector limit = shift_new;
426 Vector ephi(*map, CON, (*map).get_bvect_spher());
427 ephi.set(1).annule_hard();
428 ephi.set(2).annule_hard();
429 ephi.set(3) = 1;
430 ephi.std_spectral_base();
431 ephi.annule_domain(nz -1);
432
433 limit = bound*ssconalt + hor_rot*ephi;
434 limit.std_spectral_base(); // Boundary is fixed by value of 3 components of a vector (rather than value of potentials)
435
436 Scalar Vrb = limit(1); Vrb.set_spectral_va().ylm();
437 Scalar mmub = limit.mu(); mmub.set_spectral_va().ylm();
438 Scalar etab = limit.eta(); etab.set_spectral_va().ylm();
439
440
441
443
444 // Computing the source
445 Vector deltaA = - 2*lapse*contract(delta, 1,2, aa, 0,1);
446 Vector hijddb = - contract (shift.derive_cov(mets).derive_cov(mets), 1,2, hij, 0,1) ;
447 Vector hijddivb = - 0.3333333333333* contract (shift.divergence(mets).derive_cov(mets),0, hij,1);
448 hijddb.inc_dzpuis(); // dzpuis fixing patch...
449 hijddivb.inc_dzpuis();
450
451 Vector sourcevect2(*map,CON, (*map).get_bvect_spher());
452 sourcevect2= (2.* contract(aa, 1, lapse.derive_cov(mets),0)-12*lapse*contract(aa, 1, logpsi.derive_cov(mets), 0)) + deltaA + hijddb + hijddivb ;
453
454 sourcevect2.set(1).set_dzpuis(4);
455 sourcevect2.set(2).set_dzpuis(4);
456 sourcevect2.set(3).set_dzpuis(4);
457 sourcevect2.std_spectral_base();
458 if(sourcevect2.eta().get_etat() == ETATZERO)
459 { sourcevect2.set(2).annule_hard();}
460
461
462 double lam = (1./3.);
463
464
465
467
468 // System inversion
469 sourcevect2.poisson_boundary2(lam, shift_new, Vrb, etab, mmub, 1., 0., 1. ,0. ,1. ,0.) ;
470
471
472
474
475 // resolution tests
476 Vector source2 = contract(shift_new.derive_con(mets).derive_cov(mets), 1,2) + lam* contract(shift_new.derive_cov(mets), 0,1).derive_con(mets);
477 source2.inc_dzpuis(1);
478 // maxabs (source2 - sourcevect2);
479
480 Scalar mufin = shift_new.mu();
481 mufin.set_spectral_va().coef();
482
483 Scalar mufin2 (map_2) ;
484 mufin2.allocate_all();
485 mufin2.set_etat_qcq();
486 mufin2.std_spectral_base();
487
488 for (int k=0; k<np; k++)
489 for (int j=0; j<nt; j++) {
490
491 mufin2.set_grid_point(0, k, j, 0) = mufin.val_grid_point(1, k,j,0) - mmub.val_grid_point(1, k, j, 0);
492
493 }
494
495 // maxabs (mufin2);
496
497
498 Scalar brfin = shift_new(1);
499 brfin.set_spectral_va().coef();
500
501 Scalar brfin2 (map_2) ;
502 brfin2.allocate_all();
503 brfin2.set_etat_qcq();
504 brfin2.std_spectral_base();
505
506 for (int k=0; k<np; k++)
507 for (int j=0; j<nt; j++) {
508
509 brfin2.set_grid_point(0, k, j, 0) = brfin.val_grid_point(1, k,j,0) - Vrb.val_grid_point(1, k, j, 0);
510
511 }
512 // maxabs (brfin2);
513
514
515
517
518 // Update during the loop
519 for (int ii=1; ii <=3; ii++){
520 shift.set(ii) = shift_new(ii)*(1-relax) + shift(ii)* relax;
521 }
522
523
524 diff_ent = max(maxabs(npsi_new - npsi )); // Convergence parameter (discutable relevance...)
525
526
527
531
532
533
535 // Tensor hij resolution ///
537
538 if (isCF == false){
539
540 if (diff_ent <=5.e-3) { // No resolution until we are close to the result.
541
542 //WARNING; parameter maybe to be changed according to convergence parameters in Poisson-Hole/Kerr1.9/pplncp.C
543 util = util+1; // Loop marker for NCF equation.
544
545
549
550 // Local convergence can be asked for hij equation. Allows to satisfy integrability conditions.
551 double precis2 = 1.e5*precis ; // Here we ask for a local convergence; this can be improved later. WARNING; parameter maybe to be changed according to convergence parameters in Poisson-Hole/Kerr1.9/pplncp.C
552
553 double diff_ent2 = 1 ; // Local convergence marker
554 double relax2 = 0.5; // Local relaxation parameter. If not 1, the determinant condition won't be satisfied on a particular iteration.
555
556
557 Sym_tensor sourcehij = hij; // Random initialization...
558
559 for(int mer2=0 ;(diff_ent2 > precis2) && (mer2<mer_max2) ; mer2++) {
560
562
563 // Calculation of the source
564
565 // The double Lie derivative term is taken care of in the subroutine.
566
567 secmembre_kerr(sourcehij);
568
569
571 //System inversion (note that no boundary condition is imposed)
572
573
574 Sym_tensor hij_new = hij;
575
576 hij_new = boundfree_tensBC (sourcehij, shift , conf_fact, lapse, hij, precis2);
577
578 cout << "maximum of convergence marker for the subiteration" << endl;
579
580 diff_ent2 = max(maxabs(hij - hij_new));
581 hij = relax2*hij_new + (1 - relax2)*hij;
582
583 cout << "mer2, diffent2" << endl;
584
585 cout << mer2 << endl;
586 cout << diff_ent2 << endl;
587
588
589 }
590
592 // Resolution tests
593
594
595 Sym_tensor gammatilde = mets.con() + hij;
596 Metric gammatilde2(gammatilde); Scalar detgam = gammatilde2.determinant();
597 // cout << "determinant of result" << endl;
598// maxabs (detgam-1.);
599
600// cout << "comment l'equation en hij est elle vérifiée?" << endl;
601
602 Sym_tensor test =contract (hij.derive_cov(mets).derive_con(mets), 2,3);
603 test.annule(nz-1, nz-1);
604 test = test - hij.derive_lie(shift).derive_lie(shift)/ ((lapse/(conf_fact*conf_fact))*(lapse/(conf_fact*conf_fact)));
605 test.annule(nz-1, nz-1);
606 Sym_tensor youps = test - sourcehij/((lapse/(conf_fact*conf_fact))*(lapse/(conf_fact*conf_fact)));
607// maxabs (youps);
608// maxabs((youps).trace(mets));
609// cout << " AAABBB" << endl;
610// maxabs((youps).compute_A());
611// maxabs((youps).compute_tilde_B());
612
613
614 }
615 }
616
617
618
619
623
624
625
627 // Global variable update after an entire loop ////
629
630
631 gamtuu = mets.con() + hij;
632 gamt = gamtuu; // Métrique.
633 gam = gamt.cov()*psi4;
634 gamma = gam.cov();
635
636
637 for (int i=1; i<=3; i++) {
638 for (int j=1; j<=3; j++) {
639
640 tmp = 0;
641 tmp = ((shift.derive_con(mets))(i,j) + (shift.derive_con(mets))(j,i) - (2./3.)*shift.divergence(mets) * mets.con()(i,j))*(1./(2.*lapse));
642 aa.set(i,j) = tmp ;
643
644 }
645 }
646
647 aa = aa - (hij.derive_lie(shift) + (2./3.)*shift.divergence(mets)*hij)*(1./(2.*lapse)); //Non conformally flat correction; we assume here dhij/dt = 0.
648
649 aa_hat = aa*psi4*sqrt(psi4); // Rescaling of traceless exrinsic curvature.
650 aa_hat.std_spectral_base();
651
652
653 hatA = aa_hat; // Probably obsolete very soon... replace aa_hat by hatA.
654
655 Sym_tensor aaud = aa.up_down(gamt);
656 Sym_tensor aaud_hat = aa_hat.up_down(gamt);
657 aa_quad_scal = contract(contract (aa_hat, 0, aaud_hat, 0), 0,1);
658
659 delta = -0.5*contract( hij, 1, gamt.cov().derive_cov(mets), 2);
660
661
662 for (int i=1; i<=3; i++) {
663 for (int j=1; j<=3; j++) {
664 for (int k=1; k<=3; k++) {
665 tmp = 0.;
666 tmp = -0.5 *(gamt.cov().derive_con(mets))(i,j,k);
667 for (int l=1; l<=3; l++) {
668 tmp += -0.5*( gamt.cov()(i,l)*(hij.derive_cov(mets))(k,l,j) + gamt.cov()(l,j)*(hij.derive_cov(mets))(k,l,i));
669 }
670 delta.set(k,i,j) += tmp ;
671 }
672 }
673 }
674
675 Rstar =
676 0.25 * contract(gamt.con(), 0, 1,
677 contract(gamt.con().derive_cov(mets), 0, 1, gamt.cov().derive_cov(mets), 0, 1), 0, 1 )
678 - 0.5 * contract(gamt.con(), 0, 1,
679 contract(gamt.con().derive_cov(mets), 0, 1, gamt.cov().derive_cov(mets), 0, 2), 0, 1 ) ;
680
681 kuu = aa/(psi4);
682 kdd = contract (gamma, 0, contract(gamma, 1, kuu, 0),1);
683
684
685
687
688 // Convergence markers
689 cout << "diffent" << endl;
690 cout<< diff_ent << endl;
691 cout <<"mer" << mer << endl;
692
693
694
695
699
700 //------------------------------------------------
701 // Check of Einstein equations (3+1 form)
702 //------------------------------------------------
703
704
705 // Lapse
706 //------
707 Scalar lapse2(*map) ;
708 lapse2 = lapse ;
709 lapse2.std_spectral_base();
710
711 // 3-metric
712 //---------
713 // const Metric gam(mets.cov()*psi4) ;
714
715 Sym_tensor gamt2(*map, COV, (*map).get_bvect_spher());
716 for (int i=1; i<=3; i++)
717 for (int j=1; j<=3; j++)
718 { gamt2.set(i,j)=gam.cov()(i,j);
719 }
720 //Shift
721 //-----
722 Vector beta = shift ;
723
724 // Extrinsic curvature
725 //--------------------
726
727 Scalar TrK3(*map);
728 Sym_tensor k_uu = aa/(psi4) ;
729 k_uu.dec_dzpuis(k_uu(1,1).get_dzpuis());
730 Sym_tensor k_dd = k_uu.up_down(gam); // Another way of computing the same thing, just to be sure...
731
732 TrK3 = k_uu.trace(gam);
733 // TrK3.spectral_display("TraceKvraie", 1.e-10);
734
735 // Hamiltonian constraint
736 //-----------------------
737 Scalar ham_constr = gam.ricci_scal() ;
738 ham_constr.dec_dzpuis(3) ;
739 ham_constr += TrK3*TrK3 - contract(k_uu, 0, 1, k_dd, 0, 1) ;
740 // maxabs(ham_constr, "Hamiltonian constraint: ") ;
741
742 ham_constr.set_spectral_va().ylm();
743 // ham_constr.spectral_display("ham_constr", 1.e-9);
744
745
746
747 // Momentum constraint
748 //-------------------
749 Vector mom_constr = k_uu.divergence(gam) - TrK3.derive_con(gam) ;
750 mom_constr.dec_dzpuis(2) ;
751// maxabs(mom_constr, "Momentum constraint: ") ;
752// mom_constr(1).spectral_display("mom1", 1.e-9) ;
753// mom_constr(2).spectral_display("mom2", 1.e-9) ;
754// mom_constr(3).spectral_display("mom3", 1.e-9) ;
755
756
757 // Evolution equations
758 //--------------------
759 Sym_tensor evol_eq = lapse2*gam.ricci()
760 - lapse2.derive_cov(gam).derive_cov(gam);
761 evol_eq.dec_dzpuis() ;
762 evol_eq += k_dd.derive_lie(beta) ;
763 evol_eq.dec_dzpuis(2) ;
764 evol_eq += lapse2*(TrK3*k_dd - 2*contract(k_dd, 1, k_dd.up(0, gam), 0) ) ;
765// maxabs(evol_eq, "Evolution equations: ") ;
766
767// evol_eq.trace(gam).spectral_display("evoltrace", 1.e-10);
768// maxabs (evol_eq.trace(gam));
769
770
771 // evol_eq.spectral_display("evol", 1.e-10);
772
773
774 }
775
776
777
778 cout << "================================================================" << endl;
779 cout << " THE ITERATION HAS NOW CONVERGED" << endl;
780 cout << "================================================================" << endl;
781 }
782 return;
783}
784
785
786
787
788
789
790
791
792
793
794
795}
Active physical coordinates and mapping derivatives.
Definition coord.h:90
Scalar boundNoK
Indicates if the boundary value for the lapse or the surface gravity is used.
Definition isol_hole.h:65
Scalar lapse
Lapse function.
Definition isol_hole.h:76
bool NorKappa
Rotation rate of the horizon in the azimuthal direction.
Definition isol_hole.h:61
Sym_tensor hij
Deviation tensor( non-flat part of the conformal 3-metric on the slice; see Bonazzola et al.
Definition isol_hole.h:87
void secmembre_kerr(Sym_tensor &source_hh)
Computes the rhs of hyperbolic equation for conformal metric assuming statioarity; WARNING; up to now...
Vector shift
Shift vector.
Definition isol_hole.h:82
Sym_tensor hatA
Rescaled tracefree extrinsic curvature tensor: see Cordero et al.
Definition isol_hole.h:92
const Map & mp
Mapping associated with the star.
Definition isol_hole.h:55
bool isCF
Stores the boundary value of the lapse or surface gravity.
Definition isol_hole.h:69
void compute_stat_metric(double precis, double relax, int mer_max, int mer_max2, bool isvoid=true)
Computes a quasi-stationary 3-slice from the chosen parameters.
Affine radial mapping.
Definition map.h:2027
Flat metric for tensor calculation.
Definition metric.h:261
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Metric for tensor calculation.
Definition metric.h:90
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition metric.C:290
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition metric.C:280
const Sym_tensor & ricci() const
Returns the Ricci tensor (given by the Connection p_connect )
Definition metric.C:338
virtual const Scalar & determinant() const
Returns the determinant.
Definition metric.C:392
virtual const Scalar & ricci_scal() const
Returns the Ricci scalar.
Definition metric.C:350
Multi-domain grid.
Definition grilles.h:273
const Mg3d * get_angu_1dom() const
Returns the pointer on the associated mono-domain angular grid.
Definition mg3d.C:494
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
This class contains the parameters needed to call the general elliptic solver.
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:353
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition scalar.C:784
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
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
const Scalar & dsdr() const
Returns of *this .
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition scalar.h:601
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition scalar.C:365
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
Scalar sol_elliptic_boundary(Param_elliptic &params, const Mtbl_cf &bound, double fact_dir, double fact_neu) const
Resolution of a general elliptic equation, putting zero at infinity and with inner boundary condition...
Definition scalar_pde.C:256
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
Definition scalar.h:684
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 void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
Sym_tensor derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
Definition sym_tensor.C:360
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
Definition sym_tensor.C:349
const Tenseur & derive_con(const Metrique &) const
Returns the contravariant derivative of *this , with respect to met .
Definition tenseur.C:1568
Tensor handling.
Definition tensor.h:288
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
void coef() const
Computes the coeffcients of *this.
Tensor field of valence 1.
Definition vector.h:188
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:316
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition vector.C:381
virtual const Scalar & mu() const
Gives the field such that the angular components of the vector are written:
virtual const Scalar & eta() const
Gives the field such that the angular components of the vector are written:
void poisson_boundary2(double lam, Vector &resu, Scalar boundvr, Scalar boundeta, Scalar boundmu, double dir_vr, double neum_vr, double dir_eta, double neum_eta, double dir_mu, double neum_mu) const
Alternative to previous poisson_boundary method for vectors ; this uses method 6 for vectorial solvin...
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
const Tensor_sym & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition tensor.C:808
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition tensor.C:1014
virtual void annule(int l_min, int l_max)
Sets the Tensor to zero in several domains.
Definition tensor.C:671
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_sym & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition tensor.C:1002
Tensor up(int ind, const Metric &gam) const
Computes a new tensor by raising an index of *this.
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:654
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition tensor.C:926
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.