LORENE
Excised_slice/isol_hole_compute_metric.C
1/*
2 * Method of class Isol_Hole to compute metric data associated to a quasistationary single
3 * black hole spacetime slice.
4 *
5 * (see file isol_hole.h for documentation).
6 *
7 */
8
9/*
10 * Copyright (c) 2009 Nicolas Vasset
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 as published by
16 * the Free Software Foundation; either version 2 of the License, or
17 * (at your option) any later version.
18 *
19 * LORENE is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License
25 * along with LORENE; if not, write to the Free Software
26 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27 *
28 */
29
30// Headers Lorene
31#include "param_elliptic.h"
32#include "proto.h"
33#include "excised_slice.h"
34#include "unites.h"
35
36
37namespace Lorene {
38void Excised_slice::compute_stat_metric(double precis, double Omega, bool NorKappa,
39 Scalar boundNoK, bool isCF,double relax, int mer_max,
40 int mer_max2, bool isvoid) {
41
42 // Fundamental constants and units
43 // -------------------------------
44 using namespace Unites ;
45
46 // Verifying that we are in the right case
47 assert(type_hor==1);
48
49 // Noticing the user that the iteration has started
50
51 cout << "================================================================" << endl;
52 cout << "STARTING THE MAIN ITERATION FOR COMPUTING METRIC FIELDS" << endl;
53 cout << " iteration parameters are the following: " << endl;
54 cout << " convergence precision required:" << precis << endl;
55 cout << " max number of global steps :" << mer_max << endl;
56 cout << " relaxation parameter :" << relax << endl;
57 cout << "================================================================" << endl;
58
59
60 // Construction of a multi-grid (Mg3d) and an affine mapping from the class mapping
61 // --------------------------------------------------------------------------------
62 const Map_af* map = dynamic_cast<const Map_af*>(&mp) ;
63 const Mg3d* mgrid = (*map).get_mg();
64
65 // Construct angular grid for h(theta,phi)
66 const Mg3d* g_angu = (*mgrid).get_angu_1dom() ;
67
68 const int nz = (*mgrid).get_nzone(); // Number of domains
69 int nt = (*mgrid).get_nt(1); // Number of collocation points in theta in each domain
70 const int np = (*mgrid).get_np(1);
71 const Coord& rr = (*map).r;
72 Scalar rrr (*map) ;
73 rrr = rr ;
74 rrr.std_spectral_base();
75
76 // For now the code handles only horizons at r=1, corresponding to the first shell
77 // inner boundary. This test assures this is the case with our mapping.
78 assert((rrr.val_grid_point(1,0,0,0) - 1.) <= 1.e-9);
79
80 // Angular mapping defined as well
81 //--------------------------------
82 double r_limits2[] = {rrr.val_grid_point(1,0,0,0), rrr.val_grid_point(2,0,0,0)} ;
83 const Map_af map_2(*g_angu, r_limits2); //2D mapping; check if this is useful.
84 const Metric_flat& mets = (*map).flat_met_spher() ;
85
86
87 //----------------
88 // Initializations
89 // ---------------
90 Scalar logn (*map) ; logn = log(lapse) ;
91 logn.std_spectral_base();
92
93
94 Scalar logpsi(*map) ; logpsi = log(conf_fact) ;
95 logpsi.std_spectral_base();
96 Scalar psi4 (*map) ; psi4 = conf_fact*conf_fact*conf_fact*conf_fact ;
97 Scalar npsi (*map) ; npsi =conf_fact*lapse ;
98
99
100 Scalar conf_fact_new(*map) ; conf_fact_new.annule_hard(); conf_fact_new.std_spectral_base();
101 Scalar npsi_new(*map); npsi_new.annule_hard(); npsi_new.std_spectral_base();
102
103 Vector shift_new (*map, CON, (*map).get_bvect_spher());
104 for(int i=1; i<=3; i++){
105 shift_new.set(i)=0;
106 }
107 shift_new.std_spectral_base();
108
109 // Non conformally flat variables
110 //-------------------------------
111 Sym_tensor gamtuu = mets.con() + hij;
112 Metric gamt(gamtuu);
113 Metric gam(gamt.cov()*psi4) ;
114 Sym_tensor gamma = gam.cov();
115
116
117 // Extrinsic curvature variables
118 //------------------------------
119 Sym_tensor aa(*map, CON, (*map).get_bvect_spher());
120 for (int iii= 1; iii<=3; iii++){
121 for(int j=1; j<=3; j++){
122 aa.set(iii,j)= 0;
123 }
124 }
125 aa.std_spectral_base();
126 Scalar aa_quad_scal(*map) ; aa_quad_scal = 0. ;
127
128 Sym_tensor aa_hat(*map, CON, (*map).get_bvect_spher());
129 for (int iii= 1; iii<=3; iii++){
130 for(int j=1; j<=3; j++){
131 aa_hat.set(iii,j)= 0;
132 }
133 }
134
135 Sym_tensor kuu = aa/psi4 ;
136 Sym_tensor kuu2 = aa_hat/(psi4*psi4*sqrt(psi4));
137 Sym_tensor kdd = contract (gamma, 0, contract(gamma, 1, kuu, 0),1);
138
139
140 // (2,1)-rank delta tensor: difference between ricci rotation coefficients.
141 //-------------------------------------------------------------------------
142 Tensor delta = -0.5*contract( hij, 1, gamt.cov().derive_cov(mets), 2);
143 Scalar tmp(*map);
144
145 for (int i=1; i<=3; i++) {
146 for (int j=1; j<=3; j++) {
147 for (int k=1; k<=3; k++) {
148 tmp = 0.;
149 tmp = -0.5 *(gamt.cov().derive_con(mets))(i,j,k);
150 for (int l=1; l<=3; l++) {
151 tmp += -0.5*( gamt.cov()(i,l)*(hij.derive_cov(mets))(k,l,j)
152 + gamt.cov()(l,j)*(hij.derive_cov(mets))(k,l,i));
153 }
154 delta.set(k,i,j) += tmp ;
155 }
156 }
157 }
158
159
160 // Conformal Rstar scalar(eq 61, Bonazzola et al. 2003)
161 Scalar Rstar =
162 0.25 * contract(gamt.con(), 0, 1,
163 contract(hij.derive_cov(mets), 0, 1, gamt.cov().derive_cov(mets), 0, 1), 0, 1 )
164 - 0.5 * contract(gamt.con(), 0, 1,
165 contract(hij.derive_cov(mets), 0, 1, gamt.cov().derive_cov(mets), 0, 2), 0, 1 ) ;
166
167 Scalar norm(*map);
168 Scalar norm3(*map);
169
170
171 if (isvoid == false){
172 cout <<"FAIL: case of non-void spacetime not treated yet" << endl;
173 }
174 else {
175
176 // Parameters for the iteration
177 //-----------------------------
178
179 double diff_ent = 1 ; // initialization of difference marker between two iterations;
180
181 int util = 0; // Tool used to stop tensorial iteration at any wished step "util"
182
186
187
188 for(int mer=0 ;(diff_ent > precis) && (mer<mer_max) ; mer++) {
189
190 //Global relaxation coefficient
191
192 // Scalar variables linked to the norm of normal vector to horizon.
193 norm = sqrt(1. + hij(1,1)); norm.std_spectral_base();
194 norm3 = sqrt(1. + hij(3,3)); norm3.std_spectral_base();
195
197 // Solving for (Psi-1) //
199
200
201 // Setting of the boundary
202 //------------------------
203 double diric= 0.;
204 double neum = 1.;
205
206 Vector ssalt = rrr.derive_cov(gam);
207 Vector ssaltcon = ssalt.up_down(gam);
208 Scalar ssnormalt = sqrt(contract (ssalt,0, ssaltcon, 0));
209 ssnormalt.std_spectral_base();
210
211 ssalt.annule_domain(nz-1);
212 ssalt.annule_domain(0);
213 ssaltcon.annule_domain(nz-1);
214 ssaltcon.annule_domain(0);
215
216 ssalt = ssalt/ssnormalt;
217 ssaltcon = ssaltcon/ssnormalt;
218 // \tilde{s} in the notations of Gourgoulhon and Jaramillo, 2006
219 Vector ssconalt = ssaltcon*conf_fact*conf_fact;
220 ssconalt.std_spectral_base();
221 ssconalt.annule_domain(nz-1);
222 Scalar bound3bis = -((1./conf_fact)*contract((contract(kdd,1,ssconalt,0)),0, ssconalt,0));
223
224 bound3bis.annule_domain(nz-1);
225 bound3bis += -conf_fact*ssconalt.divergence(gamt);
226 bound3bis.annule_domain(nz-1);
227 bound3bis = 0.25*bound3bis;
228 bound3bis += -contract(conf_fact.derive_cov(gamt),0,ssconalt,0) + conf_fact.dsdr();
229 bound3bis.annule_domain(nz-1);
230 bound3bis.std_spectral_base();
231 bound3bis.set_spectral_va().ylm();
232
233 Mtbl_cf *boundd3bis = bound3bis.set_spectral_va().c_cf;
234
235 // Computing the source
236 //---------------------
237 Scalar source_conf_fact(*map) ; source_conf_fact=3. ; // Pour le fun...
238 source_conf_fact.std_spectral_base();
239
240 Scalar d2logpsi = contract(conf_fact.derive_cov(mets).derive_cov(mets), 0, 1, hij, 0,1);
241 d2logpsi.inc_dzpuis(1);
242
243 source_conf_fact = -(0.125* aa_quad_scal )/(psi4*conf_fact*conf_fact*conf_fact)
244 + conf_fact* 0.125* Rstar - d2logpsi;
245
246 source_conf_fact.std_spectral_base();
247
248 if (source_conf_fact.get_etat() == ETATZERO) {
249 source_conf_fact.annule_hard() ;
250 source_conf_fact.set_dzpuis(4) ;
251 source_conf_fact.std_spectral_base() ;
252 }
253 source_conf_fact.set_spectral_va().ylm();
254
255 // System inversion
256 //-----------------
257 Param_elliptic source11(source_conf_fact);
258 // Resolution has been done for quantity Q-1,
259 // because our solver gives a vanishing solution at infinity!
260 conf_fact_new =
261 source_conf_fact.sol_elliptic_boundary(source11, *boundd3bis, diric , neum) + 1 ;
262
263 // tests for resolution
264 //---------------------
265 Scalar baba2 = (conf_fact_new-1).laplacian();
266// cout << "psi+1-resolution" << endl;
267// maxabs (baba2 - source_conf_fact);
268
269 Scalar psinewbis = conf_fact_new -1. ; psinewbis.annule_domain(nz -1);
270 psinewbis.std_spectral_base();
271 psinewbis = psinewbis.dsdr();
272 Scalar psinewfin2 (map_2) ;
273 psinewfin2.allocate_all();
274 psinewfin2.set_etat_qcq();
275 psinewfin2.std_spectral_base();
276
277 for (int k=0; k<np; k++)
278 for (int j=0; j<nt; j++) {
279 psinewfin2.set_grid_point(0, k, j, 0) =
280 psinewbis.val_grid_point(1, k,j,0) - bound3bis.val_grid_point(1, k, j, 0);
281 }
282 // maxabs (psinewfin2);
283
284
285 // Update during the loop
286 //-----------------------
287 conf_fact = conf_fact_new* (1-relax) + conf_fact* relax ;
288 psi4 = conf_fact*conf_fact*conf_fact*conf_fact;
289 logpsi = log(conf_fact) ;
290 logpsi.std_spectral_base();
291
292
293
295 // Solving for (N*Psi -1)/
297
298
299 // Setting of the boundary
300 //------------------------
301 assert (NorKappa == false) ;
302 Scalar bound(*map);
303 bound = (boundNoK)*conf_fact -1;
304 bound.annule_domain(nz -1);
305 bound.std_spectral_base();
306 bound.set_spectral_va().ylm();
307 Mtbl_cf *boundd = bound.get_spectral_va().c_cf;
308
309 diric =1;
310 neum = 0 ;
311
312 // Computing the source ...
313 //-------------------------
314 Scalar d2lognpsi = contract(npsi.derive_cov(mets).derive_cov(mets), 0, 1, hij, 0,1);
315 d2lognpsi.inc_dzpuis(1); // dzpuis correction.
316
317 Scalar source_npsi = npsi*(aa_quad_scal*(7./8.)/(psi4*psi4) + Rstar/8.) - d2lognpsi;
318 source_npsi.std_spectral_base();
319 if (source_npsi.get_etat() == ETATZERO) {
320 source_npsi.annule_hard() ;
321 source_npsi.set_dzpuis(4) ;
322 source_npsi.std_spectral_base() ;
323 }
324
325
326 // Inversion of the operator
327 //--------------------------
328 Param_elliptic source1 (source_npsi);
329 npsi_new = source_npsi.sol_elliptic_boundary(source1, *boundd, diric, neum) ;
330
331 npsi_new = npsi_new +1;
332
333
334 // Resolution tests in npsi
335 //-------------------------
336 Scalar baba = npsi_new.laplacian();
337 // cout << "resolution_npsi" << endl;
338 // maxabs (baba - source_npsi);
339 // cout << "bound_npsi" << endl;
340 Scalar npsibound2 (map_2) ;
341 npsibound2.allocate_all();
342 npsibound2.set_etat_qcq();
343 npsibound2.std_spectral_base();
344 for (int k=0; k<np; k++)
345 for (int j=0; j<nt; j++) {
346 npsibound2.set_grid_point(0, k, j, 0) =
347 npsi_new.val_grid_point(1, k,j,0) - bound.val_grid_point(1, k, j, 0) -1.;
348 }
349 // maxabs (npsibound2);
350
351
352 // Update during the loop
353 //-----------------------
354 npsi = npsi_new*(1-relax) + npsi* relax;
355 lapse = npsi/conf_fact;
356 logn = log(lapse);
357 logn.std_spectral_base();
358
359
360
362 //Resolution in Beta //
364
365
366 // Setting of the boundary
367 //------------------------
368 bound = (boundNoK)/(conf_fact*conf_fact) ;
369 bound.annule_domain(nz -1);
370
371 // Rotation parameter for spacetime
372 Scalar hor_rot(*map); hor_rot.annule_hard(); hor_rot = Omega;
373 hor_rot.std_spectral_base() ; hor_rot.mult_rsint();
374 hor_rot.annule_domain(nz -1);
375
376 Vector limit = shift_new;
377 Vector ephi(*map, CON, (*map).get_bvect_spher());
378 ephi.set(1).annule_hard();
379 ephi.set(2).annule_hard();
380 ephi.set(3) = 1;
381 ephi.std_spectral_base();
382 ephi.annule_domain(nz -1);
383
384 limit = bound*ssconalt + hor_rot*ephi;
385 // Boundary is fixed by value of 3 components of a vector (rather than value of potentials)
386 limit.std_spectral_base();
387
388 Scalar Vrb = limit(1); Vrb.set_spectral_va().ylm();
389 Scalar mmub = limit.mu(); mmub.set_spectral_va().ylm();
390 Scalar etab = limit.eta(); etab.set_spectral_va().ylm();
391
392 // Computing the source
393 //---------------------
394 Vector deltaA = - 2*lapse*contract(delta, 1,2, aa, 0,1);
395 Vector hijddb = - contract (shift.derive_cov(mets).derive_cov(mets), 1,2, hij, 0,1) ;
396 Vector hijddivb =
397 - 0.3333333333333* contract (shift.divergence(mets).derive_cov(mets),0, hij,1);
398 hijddb.inc_dzpuis(); // dzpuis fixing patch...
399 hijddivb.inc_dzpuis();
400
401 Vector sourcevect2(*map,CON, (*map).get_bvect_spher());
402 sourcevect2 = 2.* contract(aa, 1, lapse.derive_cov(mets),0)
403 - 12*lapse*contract(aa, 1, logpsi.derive_cov(mets), 0)
404 + deltaA + hijddb + hijddivb ;
405
406 sourcevect2.set(1).set_dzpuis(4);
407 sourcevect2.set(2).set_dzpuis(4);
408 sourcevect2.set(3).set_dzpuis(4);
409 sourcevect2.std_spectral_base();
410 if(sourcevect2.eta().get_etat() == ETATZERO)
411 { sourcevect2.set(2).annule_hard();}
412
413 double lam = (1./3.);
414
415 // System inversion
416 //-----------------
417 sourcevect2.poisson_boundary2(lam, shift_new, Vrb, etab, mmub, 1., 0., 1. ,0. ,1. ,0.) ;
418
419
420 // resolution tests
421 //-----------------
422 Vector source2 = contract(shift_new.derive_con(mets).derive_cov(mets), 1,2)
423 + lam* contract(shift_new.derive_cov(mets), 0,1).derive_con(mets);
424 source2.inc_dzpuis(1);
425 // maxabs (source2 - sourcevect2);
426
427 Scalar mufin = shift_new.mu();
428 mufin.set_spectral_va().coef();
429
430 Scalar mufin2 (map_2) ;
431 mufin2.allocate_all();
432 mufin2.set_etat_qcq();
433 mufin2.std_spectral_base();
434
435 for (int k=0; k<np; k++)
436 for (int j=0; j<nt; j++) {
437 mufin2.set_grid_point(0, k, j, 0) =
438 mufin.val_grid_point(1, k,j,0) - mmub.val_grid_point(1, k, j, 0);
439 }
440 // maxabs (mufin2);
441
442 Scalar brfin = shift_new(1);
443 brfin.set_spectral_va().coef();
444
445 Scalar brfin2 (map_2) ;
446 brfin2.allocate_all();
447 brfin2.set_etat_qcq();
448 brfin2.std_spectral_base();
449
450 for (int k=0; k<np; k++)
451 for (int j=0; j<nt; j++) {
452 brfin2.set_grid_point(0, k, j, 0) =
453 brfin.val_grid_point(1, k,j,0) - Vrb.val_grid_point(1, k, j, 0);
454 }
455 // maxabs (brfin2);
456
457
458 // Update during the loop
459 //-----------------------
460 for (int ii=1; ii <=3; ii++){
461 shift.set(ii) = shift_new(ii)*(1-relax) + shift(ii)* relax;
462 }
463
464 diff_ent = max(maxabs(npsi_new - npsi )); // Convergence parameter (discutable relevance...)
465
466
468 // Tensor hij resolution ///
470
471 if (isCF == false){
472
473 if (diff_ent <=5.e-3) { // No resolution until we are close to the result.
474
475 //WARNING; parameter maybe to be changed according to convergence parameters
476 //in Poisson-Hole/Kerr1.9/pplncp.C
477 util = util+1; // Loop marker for NCF equation.
478
479
483
484 // Local convergence can be asked for hij equation.
485 // Allows to satisfy integrability conditions.
486
487 // Here we ask for a local convergence; this can be improved later.
488 // WARNING; parameter maybe to be changed according to convergence parameters
489 // in Poisson-Hole/Kerr1.9/pplncp.C
490 double precis2 = 1.e5*precis ;
491
492 double diff_ent2 = 1 ; // Local convergence marker
493 // Local relaxation parameter. If not 1, the determinant condition won't be satisfied
494 // on a particular iteration.
495 double relax2 = 0.5;
496
497 Sym_tensor sourcehij = hij; // Random initialization...
498 // cuts off high spherical harmonics with threshold being the last argument.
499 // coupe_l_tous (hij, aa, nn, ppsi, bb, nt, 6);
500
501
502 for(int mer2=0 ;(diff_ent2 > precis2) && (mer2<mer_max2) ; mer2++) {
503
504 // Calculation of the source
505 //--------------------------
506
507 // The double Lie derivative term is taken care of in the subroutine.
508 secmembre_kerr(sourcehij);
509
510 //System inversion (note that no boundary condition is imposed)
511 //-------------------------------------------------------------
512 Sym_tensor hij_new = hij;
513
514 hij_new = boundfree_tensBC (sourcehij, shift , conf_fact, lapse, hij, precis2);
515
516 cout << "maximum of convergence marker for the subiteration" << endl;
517
518 diff_ent2 = max(maxabs(hij - hij_new));
519 hij = relax2*hij_new + (1 - relax2)*hij;
520
521 cout << "mer2, diffent2" << endl;
522 cout << mer2 << endl;
523 cout << diff_ent2 << endl;
524 }
525
526 // Resolution tests
527 //-----------------
528
529 Sym_tensor gammatilde = mets.con() + hij;
530 Metric gammatilde2(gammatilde); Scalar detgam = gammatilde2.determinant();
531 // cout << "determinant of result" << endl;
532 // maxabs (detgam-1.);
533 // cout << "comment l'equation en hij est elle verifiee?" << endl;
534
535 Sym_tensor test =contract (hij.derive_cov(mets).derive_con(mets), 2,3);
536 test.annule(nz-1, nz-1);
537 test = test
539 / ((lapse/(conf_fact*conf_fact))*(lapse/(conf_fact*conf_fact)));
540 test.annule(nz-1, nz-1);
541 Sym_tensor youps = test
542 - sourcehij
543 / ((lapse/(conf_fact*conf_fact))*(lapse/(conf_fact*conf_fact)));
544 // maxabs (youps);
545 // maxabs((youps).trace(mets));
546 // cout << " AAABBB" << endl;
547 // maxabs((youps).compute_A());
548 // maxabs((youps).compute_tilde_B());
549 }
550 }
551
552
554 // Global variable update after an entire loop ////
556
557 gamtuu = mets.con() + hij;
558 gamt = gamtuu; // Metric
559 gam = gamt.cov()*psi4;
560 gamma = gam.cov();
561
562 for (int i=1; i<=3; i++) {
563 for (int j=1; j<=3; j++) {
564 tmp = 0;
565 tmp = ((shift.derive_con(mets))(i,j) + (shift.derive_con(mets))(j,i)
566 - (2./3.)*shift.divergence(mets) * mets.con()(i,j))*(1./(2.*lapse));
567 aa.set(i,j) = tmp ;
568 }
569 }
570
571 //Non conformally flat correction; we assume here dhij/dt = 0.
572 aa = aa - (hij.derive_lie(shift) + (2./3.)*shift.divergence(mets)*hij)*(1./(2.*lapse));
573
574 aa_hat = aa*psi4*sqrt(psi4); // Rescaling of traceless exrinsic curvature.
575 aa_hat.std_spectral_base();
576
577 hatA = aa_hat; // Probably obsolete very soon... replace aa_hat by hatA.
578
579 Sym_tensor aaud = aa.up_down(gamt);
580 Sym_tensor aaud_hat = aa_hat.up_down(gamt);
581 aa_quad_scal = contract(contract (aa_hat, 0, aaud_hat, 0), 0,1);
582
583 delta = -0.5*contract( hij, 1, gamt.cov().derive_cov(mets), 2);
584
585 for (int i=1; i<=3; i++) {
586 for (int j=1; j<=3; j++) {
587 for (int k=1; k<=3; k++) {
588 tmp = 0.;
589 tmp = -0.5 *(gamt.cov().derive_con(mets))(i,j,k);
590 for (int l=1; l<=3; l++) {
591 tmp += -0.5*( gamt.cov()(i,l)*(hij.derive_cov(mets))(k,l,j)
592 + gamt.cov()(l,j)*(hij.derive_cov(mets))(k,l,i));
593 }
594 delta.set(k,i,j) += tmp ;
595 }
596 }
597 }
598
599 Rstar =
600 0.25 * contract(gamt.con(), 0, 1,
601 contract(gamt.con().derive_cov(mets), 0, 1,
602 gamt.cov().derive_cov(mets), 0, 1), 0, 1 )
603 - 0.5 * contract(gamt.con(), 0, 1,
604 contract(gamt.con().derive_cov(mets), 0, 1,
605 gamt.cov().derive_cov(mets), 0, 2), 0, 1 ) ;
606 kuu = aa/(psi4);
607 kdd = contract (gamma, 0, contract(gamma, 1, kuu, 0),1);
608
609 // Convergence markers
610 //--------------------
611 cout << "diffent" << endl;
612 cout<< diff_ent << endl;
613 cout <<"mer" << mer << endl;
614
615
616 //------------------------------------------------
617 // Check of Einstein equations (3+1 form)
618 //------------------------------------------------
619
620
621 // Lapse
622 //------
623 Scalar lapse2(*map) ;
624 lapse2 = lapse ;
625 lapse2.std_spectral_base();
626
627 // 3-metric
628 //---------
629 // const Metric gam(mets.cov()*psi4) ;
630
631 Sym_tensor gamt2(*map, COV, (*map).get_bvect_spher());
632 for (int i=1; i<=3; i++)
633 for (int j=1; j<=3; j++)
634 { gamt2.set(i,j)=gam.cov()(i,j);
635 }
636
637 //Shift
638 //-----
639 Vector beta = shift ;
640
641 // Extrinsic curvature
642 //--------------------
643 Scalar TrK3(*map);
644 Sym_tensor k_uu = aa/(psi4) ;
645 k_uu.dec_dzpuis(k_uu(1,1).get_dzpuis());
646 // Another way of computing the same thing, just to be sure...
647 Sym_tensor k_dd = k_uu.up_down(gam);
648 TrK3 = k_uu.trace(gam);
649 // TrK3.spectral_display("TraceKvraie", 1.e-10);
650
651 // Hamiltonian constraint
652 //-----------------------
653 Scalar ham_constr = gam.ricci_scal() ;
654 ham_constr.dec_dzpuis(3) ;
655 ham_constr += TrK3*TrK3 - contract(k_uu, 0, 1, k_dd, 0, 1) ;
656 // maxabs(ham_constr, "Hamiltonian constraint: ") ;
657
658 ham_constr.set_spectral_va().ylm();
659 // ham_constr.spectral_display("ham_constr", 1.e-9);
660
661 // Momentum constraint
662 //-------------------
663 Vector mom_constr = k_uu.divergence(gam) - TrK3.derive_con(gam) ;
664 mom_constr.dec_dzpuis(2) ;
665 // maxabs(mom_constr, "Momentum constraint: ") ;
666 // mom_constr(1).spectral_display("mom1", 1.e-9) ;
667 // mom_constr(2).spectral_display("mom2", 1.e-9) ;
668 // mom_constr(3).spectral_display("mom3", 1.e-9) ;
669
670 // Evolution equations
671 //--------------------
672 Sym_tensor evol_eq = lapse2*gam.ricci() - lapse2.derive_cov(gam).derive_cov(gam);
673 evol_eq.dec_dzpuis() ;
674 evol_eq += k_dd.derive_lie(beta) ;
675 evol_eq.dec_dzpuis(2) ;
676 evol_eq += lapse2*(TrK3*k_dd - 2*contract(k_dd, 1, k_dd.up(0, gam), 0) ) ;
677 // maxabs(evol_eq, "Evolution equations: ") ;
678 // evol_eq.trace(gam).spectral_display("evoltrace", 1.e-10);
679 // maxabs (evol_eq.trace(gam));
680 // evol_eq.spectral_display("evol", 1.e-10);
681
682 }
683
684 cout << "================================================================" << endl;
685 cout << " THE ITERATION HAS NOW CONVERGED" << endl;
686 cout << "================================================================" << endl;
687 }
688 return;
689}
690}
Active physical coordinates and mapping derivatives.
Definition coord.h:90
Sym_tensor hatA
Rescaled tracefree extrinsic curvature tensor: rescaled same way as Cordero et al.
void compute_stat_metric(double precis, double Omega_i, bool NorKappa_i, Scalar NoK_i, bool isCF_i=true, double relax=0.2, int mer_max=2000, int mer_max2=200, bool isvoid=true)
If hor_type=1, computes a quasi-stationary 3-slice from the chosen parameters.
int type_hor
Chosen horizon type.
void secmembre_kerr(Sym_tensor &source_hh)
If hor_type=1, computes the rhs of hyperbolic equation for conformal metric assuming stationarity; WA...
Sym_tensor hij
Deviation tensor( non-flat part of the conformal 3-metric on the slice; see Bonazzola et al.
Scalar lapse
Lapse function.
const Map & mp
Mapping associated with the slice.
Vector shift
Shift vector.
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.