LORENE
excision_surf.C
1/*
2 * Definition of methods for the class Excision_surf and its subclass App_hor
3 *
4 */
6// LOCAL VERSION
7// Several additional functions that may be temporary (on the testing pad...)
9
10/*
11 * Copyright (c) 2008 Jose-Luis Jaramillo & Jerome Novak & Nicolas Vasset
12 *
13 * This file is part of LORENE.
14 *
15 * LORENE is free software; you can redistribute it and/or modify
16 * it under the terms of the GNU General Public License version
17 * as published by the Free Software Foundation.
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
30char excision_surf_C[] = "$Header: /cvsroot/Lorene/C++/Source/App_hor/excision_surf.C,v 1.6 2014/10/13 08:52:38 j_novak Exp $" ;
31
32/*
33 * $Header: /cvsroot/Lorene/C++/Source/App_hor/excision_surf.C,v 1.6 2014/10/13 08:52:38 j_novak Exp $
34 *
35 */
36
37// C headers
38#include <cmath>
39
40// Lorene headers
41#include "metric.h"
42#include "spheroid.h"
43#include "utilitaires.h"
44#include "param.h"
45#include "itbl.h"
46#include "map.h"
47#include <cassert>
48#include "nbr_spx.h"
49#include "math.h"
50#include "param_elliptic.h"
51#include "tensor.h"
52#include "sym_tensor.h"
53#include "diff.h"
54#include "proto.h"
55#include "param.h"
56#include "excision_surf.h"
57
58//---------------//
59// Constructors //
60//--------------//
61
62
63namespace Lorene {
64Excision_surf::Excision_surf(const Scalar& h_in, const Metric& gij, const Sym_tensor& Kij2, const Scalar& ppsi, const Scalar& nn, const Vector& beta, double timestep, int int_nos = 1):
65 sph(h_in, gij, Kij2),
66 conf_fact(ppsi),
67 lapse(nn),
68 shift(beta),
69 gamij (gij),
70 Kij(Kij2),
71 delta_t(timestep),
72 no_of_steps(int_nos),
73 expa(sph.theta_plus()),
74 dt_expa(sph.theta_plus())
75
76{
77
79
80 set_der_0x0() ;
81
82}
83
84
85
86
87
88//Copy constructor//
89
90Excision_surf::Excision_surf (const Excision_surf &exc_in) :sph(exc_in.sph),
91 conf_fact(exc_in.conf_fact),
92 lapse(exc_in.lapse),
93 shift(exc_in.shift),
94 gamij (exc_in.gamij),
95 Kij (exc_in.Kij),
96 delta_t(exc_in.delta_t),
97 no_of_steps(exc_in.no_of_steps),
98 expa(exc_in.expa),
99 dt_expa(exc_in.dt_expa)
100
101{
102 set_der_0x0() ;
103
104}
105
106
107
108
109// Assignment to another Excision_surf
111{
112
113 sph = surf_in.sph ;
114 conf_fact = surf_in.conf_fact ;
115 lapse = surf_in.lapse ;
116 shift = surf_in.shift ;
117 gamij = surf_in.gamij ;
118 Kij = surf_in.Kij ;
119 delta_t = surf_in.delta_t ;
120 no_of_steps = surf_in.no_of_steps ;
121 expa = surf_in.expa;
122 dt_expa = surf_in.dt_expa;
123 del_deriv() ; // Deletes all derived quantities
124
125}
126
127//------------//
128//Destructor //
129//-----------//
130
135
136// -----------------//
137// Memory management//
138//------------------//
140 if (p_get_BC_conf_fact_1 != 0x0) delete p_get_BC_conf_fact_1 ;
141 if (p_get_BC_lapse_1 != 0x0) delete p_get_BC_lapse_1 ;
142 if (p_get_BC_shift_1 != 0x0) delete p_get_BC_shift_1 ;
143 if (p_get_BC_Npsi_1 != 0x0) delete p_get_BC_Npsi_1 ;
144 if (p_get_BC_conf_fact_2 != 0x0) delete p_get_BC_conf_fact_2 ;
145 if (p_get_BC_conf_fact_3 != 0x0) delete p_get_BC_conf_fact_3 ;
146 if (p_get_BC_conf_fact_4 != 0x0) delete p_get_BC_conf_fact_4 ;
147 if (p_get_BC_lapse_2 != 0x0) delete p_get_BC_lapse_2 ;
148 if (p_get_BC_lapse_3 != 0x0) delete p_get_BC_lapse_3 ;
149 if (p_get_BC_lapse_4 != 0x0) delete p_get_BC_lapse_4 ;
150 if (p_derive_t_expa != 0x0) delete p_derive_t_expa ;
151 if (p_get_BC_shift_2 != 0x0) delete p_get_BC_shift_2 ;
152 if (p_get_BC_shift_3 != 0x0) delete p_get_BC_shift_3 ;
153 if (p_get_BC_shift_4 != 0x0) delete p_get_BC_shift_4 ;
154 if (p_get_BC_Npsi_2 != 0x0) delete p_get_BC_Npsi_2 ;
155 if (p_get_BC_Npsi_3 != 0x0) delete p_get_BC_Npsi_3 ;
156 if (p_get_BC_Npsi_4 != 0x0) delete p_get_BC_Npsi_4 ;
157 if (p_get_BC_Npsi_5 != 0x0) delete p_get_BC_Npsi_5 ;
158 set_der_0x0() ;
159}
160
163 p_get_BC_lapse_1 = 0x0 ;
164 p_get_BC_shift_1 = 0x0 ;
165 p_get_BC_Npsi_1 = 0x0 ;
169 p_get_BC_lapse_2 = 0x0 ;
170 p_get_BC_lapse_3 = 0x0 ;
171 p_get_BC_lapse_4 = 0x0 ;
172 p_derive_t_expa = 0x0 ;
173 p_get_BC_shift_2 = 0x0 ;
174 p_get_BC_shift_3 = 0x0 ;
175 p_get_BC_shift_4 = 0x0 ;
176 p_get_BC_Npsi_2 = 0x0 ;
177 p_get_BC_Npsi_3 = 0x0 ;
178 p_get_BC_Npsi_4 = 0x0 ;
179 p_get_BC_Npsi_5 = 0x0 ;
180
181}
182
183
184
185 //--------------//
186 // Outputs //
187 //--------------//
188
189
190
191// Provides the three parameters to use for hyperbolic evolution of the expansion,
192// Using expansion and its time derivative on initial data.
193// WARNING: to be called once and for all on initial data. Re-calling it does not make sense.
194
195void Excision_surf::get_evol_params_from_ID(double alpha, double beta, double gamma, Scalar& Ee, Vector& Jj, Sym_tensor& Ss){
196 cout << "===========================================================================================" << endl;
197 cout << "Starting the routine that computes parameters for the hyperbolic evolution of the expansion" << endl;
198 cout << "Those parameters should be set once and for all: do not call that routine twice" <<endl;
199 cout << "===========================================================================================" << endl;
200
201 Scalar expa_0= sph.theta_plus();
202 expa_0.set_spectral_va().ylm();
203
204 set_expa() = expa_0;
205
206 if ((max(expa_0))(0)<=1.e-5)
207
208 {
209 cout << "=============================================================================" << endl;
210 cout << " WARNING: Routine failure..." << endl;
211 cout << "the horizon expansion is already pretty close to zero..."<< endl;
212 cout << "We advise to switch from hyperbolic evolution to parabolic evolution" << endl;
213 cout << "This routine will not do anything relevant on parameters alpha, beta, gamma" << endl;
214 cout << "=============================================================================" << endl;
215 double is_expansion_adapted = 1.;
216 assert (is_expansion_adapted ==2.);
217
218 return;
219 }
220 else{
221
222 Scalar fbN = derive_t_expa(Ee, Jj, Ss); // Basic "right" time derivative of the expansion from ID
223 set_dt_expa()=fbN;
224
225 Scalar tf_zero = 2.*fbN/expa_0; tf_zero.std_spectral_base();
226 tf_zero.set_spectral_va().ylm(); // Relevant quantity for parameter characterization (see notes).
227
228 Mtbl_cf* tfz = tf_zero.set_spectral_va().c_cf;// Mtbl_cf storing the spectral coefficients of tf_zero.
229
230 // Choice for \beta (see notes for all choices)
231 beta = 2.*(*tfz).val_in_bound_jk(0,0,0);
232
233 cout << "beta: " << beta << endl;
234
235 // Choice for \gamma
236 gamma = -beta*beta* (1./8.); // Any choice between zero and -beta*beta*(1./4.) is acceptable.
237
238 cout << "gamma: " << gamma << endl;
239
240 // Choice for \alpha, with an additional assertion.
241 alpha = 0.;
242 double alpha_aux=0;
243 double nb_l = (*tfz).set(0).get_dim(1);
244 double nb_m = (*tfz).set(0).get_dim(2);
245
246 cout << "nb_l: " << nb_l << endl;
247 cout << "nb_m: " << nb_m << endl;
248
249 if (nb_m==1){
250 alpha = 2.*fabs(beta - (*tfz).val_in_bound_jk(0,1,0));
251 }
252 else{
253 for (int ii=0; ii <nb_m; ii++){
254 alpha_aux = 2.*fabs(beta - (*tfz).val_in_bound_jk(0,1,ii));
255 if (alpha_aux >=alpha){
256 alpha = alpha_aux;
257 }
258 }
259 }
260 cout << "alpha: " << alpha << endl;
261 // Test for higher harmonics, supposedly of lower amplitude, but...
262 for (int ii=2; ii <nb_l; ii++)
263 for (int jj=0; jj <nb_m; jj++){
264 alpha_aux = 2.*(beta - (*tfz).val_in_bound_jk(0,ii,jj))/(double(ii*(ii+1)));
265 if (alpha_aux >= alpha){
266 cout << "higher Ylm harmonics are predominant in the expansion variation on the horizon!" << endl;
267 cout << "changing the values of the parameters accordingly..." << endl;
268 alpha = alpha_aux;}
269 }
270
271 return;
272}
273}
274
275//---------//
276//Accessors//
277//---------//
278
279// Source for the Neumann BC on the conformal factor, with arbitrary expansion
281
282
283 if (p_get_BC_conf_fact_1 == 0x0){
284
285 //Initialization of expa in the trivial case
286 Scalar exppa(sph.theta_plus().get_mp());
287 if ( isMOTS == true)
288 {
289 exppa.set_etat_zero();
290
291 }
292 else {
293 assert (expa.get_spectral_va().get_etat() !=0);
294 exppa = expa;
295
296
297// expa.spectral_display();
298// int tgh; cin >> tgh;
299 }
300
301 // 3d interpolation for expa
302 Scalar expa_3 (lapse.get_mp());
303
304 expa_3.allocate_all();
305 expa_3.std_spectral_base();
306
307 int nz = (*lapse.get_mp().get_mg()).get_nzone();
308 int nr = (*lapse.get_mp().get_mg()).get_nr(1);
309 int nt = (*lapse.get_mp().get_mg()).get_nt(1);
310 int np = (*lapse.get_mp().get_mg()).get_np(1);
311
312
313 for (int f= 0; f<nz; f++)
314 for (int k=0; k<np; k++)
315 for (int j=0; j<nt; j++) {
316 for (int l=0; l<nr; l++) {
317
318 expa_3.set_grid_point(f,k,j,l) = exppa.val_grid_point(0,k,j,0);
319
320 }
321 }
322
323 if (nz >2){
324
325 expa_3.annule_domain(0);
326 expa_3.annule_domain(nz - 1);
327 }
328 expa_3.std_spectral_base();
329 expa_3.set_spectral_va().ylm();
330 // End Mapping interpolation
331
332
333 Sym_tensor gamconfcov = gamij.cov()/pow(conf_fact, 4);
334 gamconfcov.std_spectral_base();
335 Metric gamconf(gamconfcov);
336 Vector tilde_s = gamconf.radial_vect();
337 Scalar bound_psi = -((1./conf_fact)*contract((contract(Kij,1,tilde_s,0)),0, tilde_s,0));
338 bound_psi.annule_domain(nz-1);
339 bound_psi += -conf_fact*tilde_s.divergence(gamconf);
340 bound_psi.annule_domain(nz-1);
341 bound_psi += (conf_fact*conf_fact*conf_fact)*Kij.trace(gamij);
342 bound_psi = 0.25*bound_psi;
343 bound_psi.annule_domain(nz-1);
344 bound_psi = (bound_psi -contract(conf_fact.derive_cov(gamconf),0,tilde_s,0)) + conf_fact.dsdr();
345 Scalar bound_psi2 = expa_3*((conf_fact*conf_fact*conf_fact))/4.;
346 // Remark: the used expansion term is actually \frac{\theta^{+}}{N}, as it is not normalized by the lapse in the Spheroid class.
347 bound_psi2.annule_domain(nz -1);
348 bound_psi = bound_psi + bound_psi2;
349 bound_psi.std_spectral_base();
350 bound_psi.set_spectral_va().ylm();
351
352 p_get_BC_conf_fact_1 = new Scalar(bound_psi);
353
354
355}
356 return *p_get_BC_conf_fact_1 ;
357
358}
359
360
361
362
363// Source for the Dirichlet BC on the conformal factor, based on a parabolic driver
364const Scalar& Excision_surf::get_BC_conf_fact_2(double c_psi_lap, double c_psi_fin, Scalar& expa_fin) const{
365 if (p_get_BC_conf_fact_2 == 0x0){
366
367 // Definition of ff
368 // ================
369
370 // Start Mapping interpolation
371 Scalar thetaplus = sph.theta_plus();
372 Scalar theta_plus3 (lapse.get_mp());
373
374 theta_plus3.allocate_all();
375 theta_plus3.std_spectral_base();
376
377 Scalar expa_fin3 (lapse.get_mp());
378
379 expa_fin3.allocate_all();
380 expa_fin3.std_spectral_base();
381
382 int nz = (*lapse.get_mp().get_mg()).get_nzone();
383 int nr = (*lapse.get_mp().get_mg()).get_nr(1);
384 int nt = (*lapse.get_mp().get_mg()).get_nt(1);
385 int np = (*lapse.get_mp().get_mg()).get_np(1);
386
387
388 for (int f= 0; f<nz; f++)
389 for (int k=0; k<np; k++)
390 for (int j=0; j<nt; j++) {
391 for (int l=0; l<nr; l++) {
392
393 theta_plus3.set_grid_point(f,k,j,l) = thetaplus.val_grid_point(0,k,j,0);
394 expa_fin3.set_grid_point(f,k,j,l) = expa_fin.val_grid_point(0,k,j,0);
395
396 }
397 }
398 if (nz >2){
399 theta_plus3.annule_domain(0);
400 theta_plus3.annule_domain(nz - 1);
401 expa_fin3.annule_domain(0);
402 expa_fin3.annule_domain(nz - 1);
403 }
404 // End Mapping interpolation
405
406
407 Scalar ff = lapse*(c_psi_lap*theta_plus3.lapang() + c_psi_fin*(theta_plus3 - expa_fin3));
408 ff.std_spectral_base();
409
410 // Definition of k_1
411 Scalar k_1 =delta_t*ff;
412
413 // Intermediate value of the expansion, for Runge-Kutta 2nd order scheme
414 Scalar psi_int = conf_fact + k_1; psi_int.std_spectral_base();
415 Sym_tensor gamconfcov = gamij.cov()/pow(psi_int, 4); // think about the consistency of redifining the conformal metric
416 // since in this manner the unimodular conditions is lost
417 gamconfcov.std_spectral_base();
418 Metric gamconf(gamconfcov);
419 Vector tilde_s = gamconf.radial_vect();
420 Scalar theta_int = ((1./psi_int)*contract((contract(Kij,1,tilde_s,0)),0, tilde_s,0));
421 theta_int += psi_int*tilde_s.divergence(gamconf);
422 theta_int += 4.*contract(psi_int.derive_cov(gamconf),0,tilde_s,0);
423 theta_int = theta_int/pow(psi_int,3);
424 theta_int += -Kij.trace(gamij);
425 theta_int.std_spectral_base();
426
427 // Recalculation of ff with intermediate values.
428 Scalar ff_int = lapse*(c_psi_lap*theta_int.lapang() + c_psi_fin*(theta_int - expa_fin3));
429
430 // Definition of k_2
431 Scalar k_2 = delta_t*ff_int;
432 k_2.std_spectral_base();
433
434 // Result of RK2 evolution
435 if (nz >2){
436 k_2.annule_domain(nz-1);}
437 Scalar bound_psi = conf_fact + k_2;
438 bound_psi.std_spectral_base();
439 bound_psi.set_spectral_va().ylm();
440
441 // Assignment of output
442 p_get_BC_conf_fact_2 = new Scalar(bound_psi);
443
444}
445 return *p_get_BC_conf_fact_2 ;
446}
447
448
449
450
451// Source for the Neumann BC on the conformal factor, based on a parabolic driver for the expansion
452const Scalar& Excision_surf::get_BC_conf_fact_3(double c_theta_lap, double c_theta_fin, Scalar& expa_fin) const{
453 if (p_get_BC_conf_fact_3 == 0x0){
454
455 // Definition of ff
456 // ================
457
458 // Start Mapping interpolation
459 Scalar thetaplus = sph.theta_plus();
460 Scalar theta_plus3 (lapse.get_mp());
461
462 theta_plus3.allocate_all();
463 theta_plus3.std_spectral_base();
464
465 Scalar expa_fin3 (lapse.get_mp());
466
467 expa_fin3.allocate_all();
468 expa_fin3.std_spectral_base();
469
470 int nz = (*lapse.get_mp().get_mg()).get_nzone();
471 int nr = (*lapse.get_mp().get_mg()).get_nr(1);
472 int nt = (*lapse.get_mp().get_mg()).get_nt(1);
473 int np = (*lapse.get_mp().get_mg()).get_np(1);
474
475
476 for (int f= 0; f<nz; f++)
477 for (int k=0; k<np; k++)
478 for (int j=0; j<nt; j++) {
479 for (int l=0; l<nr; l++) {
480
481 theta_plus3.set_grid_point(f,k,j,l) = thetaplus.val_grid_point(0,k,j,0);
482 expa_fin3.set_grid_point(f,k,j,l) = expa_fin.val_grid_point(0,k,j,0);
483
484 }
485 }
486 if (nz >2){
487 theta_plus3.annule_domain(0);
488 theta_plus3.annule_domain(nz - 1);
489 expa_fin3.annule_domain(0);
490 expa_fin3.annule_domain(nz - 1);
491 }
492
493 // End Mapping interpolation
494// cout << "convergence?" << endl;
495// cout << expa_fin3.val_grid_point(1,0,0,0) << endl;
496// cout << theta_plus3.val_grid_point(1,0,0,0) << endl;
497
498
499 Scalar ff = lapse*(c_theta_lap*theta_plus3.lapang() + c_theta_fin*(theta_plus3 - expa_fin3));
500 ff.std_spectral_base();
501
502 // Definition of k_1
503 Scalar k_1 =delta_t*ff;
504
505 // Intermediate value of the expansion, for Runge-Kutta 2nd order scheme
506 Scalar theta_int = theta_plus3 + k_1; theta_int.std_spectral_base();
507
508 // Recalculation of ff with intermediate values.
509 Scalar ff_int = lapse*(c_theta_lap*theta_int.lapang() + c_theta_fin*(theta_int - expa_fin3));
510
511 // Definition of k_2
512 Scalar k_2 = delta_t*ff_int;
513 k_2.std_spectral_base();
514
515 // Result of RK2 evolution
516 Scalar bound_theta = theta_plus3 + k_2;
517 bound_theta.std_spectral_base();
518
519 // Calculation of Neumann BC for Psi
520 Scalar bound_psi = get_BC_conf_fact_1(true) + bound_theta*pow(conf_fact,3)/4.;
521 bound_psi.std_spectral_base();
522 bound_psi.set_spectral_va().ylm();
523
524 // Assignment of output
525 p_get_BC_conf_fact_3 = new Scalar(bound_psi);
526
527}
528 return *p_get_BC_conf_fact_3 ;
529}
530
531// Source for the Dirchlet BC on the conformal factor, based on the consistency condition derived from the trace
532// of the kinematical relation
534 if (p_get_BC_conf_fact_4 == 0x0){
535
536 // Definition of ff
537 // ================
538 const Metric_flat& flat = lapse.get_mp().flat_met_spher() ;
539
540 int nz = (*lapse.get_mp().get_mg()).get_nzone();
541 Scalar ff = contract(shift, 0, conf_fact.derive_cov(flat), 0)
542 + 1./6. * (conf_fact * (shift.divergence(flat) - lapse*Kij.trace(gamij))) ; // Add he N K term
543 // Divergence with respect to the conformal metric coincides with divergence withh respect to the
544 // flat metric (from the unimodular condition on the conformal metric)
545 // In this way, we do not need to recalculate a conformal metric in the intermediate RK step that would violated
546 // the unimodular condition
547
548 ff.std_spectral_base() ;
549
550 // Definition of k_1
551 Scalar k_1 =delta_t*ff;
552 k_1.annule_domain(nz-1);
553
554
555 // Intermediate value of the expansion, for Runge-Kutta 2nd order scheme
556 Scalar psi_int = conf_fact + k_1; psi_int.std_spectral_base();
557 psi_int.annule_domain(nz-1);
558
559 // Recalculation of ff with intermediate values.
560 Scalar ff_int = contract(shift, 0, psi_int.derive_cov(flat), 0)
561 + 1./6. * psi_int * (shift.divergence(flat) - lapse*Kij.trace(gamij)) ; // Add he N K term
562
563
564 // Definition of k_2
565 Scalar k_2 = delta_t*ff_int;
566
567 // k_2 = k_2*1000; //TO REMOVE
568
569 k_2.std_spectral_base();
570
571 k_2.annule_domain(nz-1);
572 // Result of RK2 evolution
573 Scalar bound_psi = conf_fact + k_2;
574
575
576 bound_psi.std_spectral_base();
577 bound_psi.set_spectral_va().ylm();
578
579 // Assignment of output
580 p_get_BC_conf_fact_4 = new Scalar(bound_psi);
581
582}
583 return *p_get_BC_conf_fact_4 ;
584}
585
586
587
588
589// Source for the Dirichlet BC on the lapse
590const Scalar& Excision_surf::get_BC_lapse_1(double value) const{
591 if (p_get_BC_lapse_1 == 0x0){
592
593 Scalar bound_lapse(lapse.get_mp());
594 bound_lapse = value; bound_lapse.std_spectral_base();
595 bound_lapse.set_spectral_va().ylm();
596 p_get_BC_lapse_1 = new Scalar(bound_lapse);
597
598}
599 return *p_get_BC_lapse_1 ;
600}
601
602// Source for Dirichlet BC on the lapse, based on a parabolic driver
603const Scalar& Excision_surf::get_BC_lapse_2(double lapse_fin, double c_lapse_lap, double c_lapse_fin) const{
604 if (p_get_BC_lapse_2 == 0x0){
605
606
607
608 Scalar ff = lapse*(c_lapse_lap*lapse.lapang() + c_lapse_fin*lapse);
609 ff.std_spectral_base();
610 ff = ff -lapse*c_lapse_fin*lapse_fin;
611 ff.std_spectral_base();
612
613
614 // Definition of k_1
615 Scalar k_1 =delta_t*ff;
616
617 // Intermediate value of lapse, for Runge-Kutta 2nd order scheme
618 Scalar lapse_int = lapse + k_1; lapse_int.std_spectral_base();
619
620 // Recalculation of ff with intermediate values.
621 Scalar ff_int = lapse*(c_lapse_lap*lapse_int.lapang() + c_lapse_fin*lapse_int);
622 ff_int.std_spectral_base();
623 ff_int = ff_int -lapse*c_lapse_fin*lapse_fin;
624 ff_int.std_spectral_base();
625
626
627 // Definition of k_2
628 Scalar k_2 = delta_t*ff_int;
629 k_2.std_spectral_base();
630
631 // Result of RK2 evolution
632 Scalar bound_lapse = lapse + k_2;
633 bound_lapse.std_spectral_base();
634 bound_lapse.set_spectral_va().ylm();
635
636
637 p_get_BC_lapse_2 = new Scalar(bound_lapse);
638
639}
640 return *p_get_BC_lapse_2 ;
641}
642
643
644
645
646// Source for Dirichlet BC on the lapse, based on einstein equations!! dttheta is 2d scalar, and matter quantities are 3d.
647const Scalar& Excision_surf::get_BC_lapse_3(Scalar& dttheta, Scalar& Ee, Vector& Jj, Sym_tensor& Sij, bool sph_sym) const{
648 if (p_get_BC_lapse_3 == 0x0){
649 int nz2 = (*lapse.get_mp().get_mg()).get_nzone();
650 // Radial vector for the full 3-metric.
651 Vector sss = gamij.radial_vect();
652 Vector sss_down = sss.up_down(gamij);
653 Scalar bb = contract (shift,0, sss_down,0);
654
655 Scalar bound_N3(lapse.get_mp()); bound_N3.allocate_all(); bound_N3.std_spectral_base();// Result to be stored there.
656
657 Scalar Kappa = bb*contract(contract(Kij,0,sss,0),0,sss,0) - contract(sss,0, lapse.derive_cov(gamij),0);
658 Scalar Matter = 8.*M_PI*(bb*Ee);
659 Matter.annule_domain(nz2-1);
660 Matter += - 8.*M_PI*(bb + lapse)*contract(Jj,0,sss_down,0);
661 Matter.annule_domain(nz2-1);
662 Matter += 8.*M_PI*lapse*contract(contract(Sij,0,sss_down,0),0,sss_down,0);
663
664
665 // 2d interpolation of the Kappa constant and the matter terms.
666
667 Scalar Kappa2 = sph.get_ricci();
668 Kappa2.annule_hard();
669 Scalar bb2 = Kappa2;
670 Scalar Matter2 = Kappa2;
671 Scalar nn2 = Kappa2;
672
673 const Metric_flat& flat2 = Kappa2.get_mp().flat_met_spher();
674
675 int nt = (*Kappa2.get_mp().get_mg()).get_nt(0);
676 int np = (*Kappa2.get_mp().get_mg()).get_np(0);
677 for (int k=0; k<np; k++)
678 for (int j=0; j<nt; j++) {
679 Kappa2.set_grid_point(0,k,j,0) = Kappa.val_grid_point(1,k,j,0);
680 bb2.set_grid_point(0,k,j,0) = bb.val_grid_point(1,k,j,0);
681 Matter2.set_grid_point(0,k,j,0) = Matter.val_grid_point(1,k,j,0);
682 nn2.set_grid_point(0,k,j,0) = lapse.val_grid_point(1,k,j,0);
683 }
684
685 Kappa2.std_spectral_base();
686 bb2.std_spectral_base();
687 Matter2.std_spectral_base();
688
689 Scalar Tplus = sph.theta_plus();
690 Scalar Tminus = sph.theta_minus();
691 Scalar Ricci = sph.get_ricci();
692 Vector Ll = sph.get_ll();
693 Sym_tensor Shear = sph.shear();
694 Metric qab = sph.get_qab();
695
696 Scalar source = 2.*contract(Ll.up_down(qab),0,sph.derive_cov2d(bb2),0) + bb2*(contract(sph.derive_cov2d(Ll.up_down(qab)),0,1)+ contract(Ll,0,Ll.up_down(qab),0)- 0.5*(Ricci + Tplus*Tminus) + 0.25*Tplus*Tplus + 0.5*contract(Shear.up_down(qab),0,1,Shear,0,1)) + Matter2 + Tplus*Kappa2 + dttheta;
697
698 Scalar lapb = contract(sph.derive_cov2d(sph.derive_cov2d(bb2)).up(0,qab),0,1);
699
700 source = source + lapb;
701
702
703 Scalar source_add = 2.*contract(Ll.up_down(qab),0,sph.derive_cov2d(nn2),0);
704
705 source = source - source_add;
706
707 Scalar sqrt_q_h2 = sph.sqrt_q() * sph.get_hsurf() * sph.get_hsurf() ;
708 Tensor Delta2 = sph.delta();
709 Sym_tensor qab_con = sph.get_qab().con() / (sph.get_hsurf() * sph.get_hsurf()) ;
710 qab_con.set(1,1) = 1. ;
711 qab_con.set(1,2) = 0. ;
712 qab_con.set(1,3) = 0. ;
713 qab_con.std_spectral_base() ;
714
715 Sym_tensor hab =(qab_con*sqrt_q_h2 - flat2.con()) / (sph.get_hsurf()*sph.get_hsurf()) ;
716 hab.set(1,1) = 1. ;
717 hab.set(1,2) = 0. ;
718 hab.set(1,3) = 0. ;
719 hab.std_spectral_base() ;
720
721 Vector dN = sph.derive_cov2dflat(nn2);
722 Tensor ddN = sph.derive_cov2dflat(dN);
723
724 Scalar lap_rem = contract(qab_con, 0,1, contract(Delta2,0,dN,0),0,1)*sqrt_q_h2 - sph.get_hsurf()*sph.get_hsurf()*contract(hab,0,1,ddN,0,1);// What remains of the laplacian
725
726 source = sqrt_q_h2*source + lap_rem;
727
728 Scalar bound_N=nn2;
729
730 if (sph_sym == true){
731 Scalar lapang_s_par = (contract(sph.derive_cov2d(Ll.up_down(qab)),0,1)+ contract(Ll,0,Ll.up_down(qab),0)- 0.5*(Ricci + Tplus*Tminus) - 0.25*Tplus*Tplus - 0.5*contract(Shear.up_down(qab),0,1,Shear,0,1))*sqrt_q_h2;
732 double lapang_par = lapang_s_par.val_grid_point(0,0,0,0);
733
734 bound_N = source.poisson_angu(lapang_par);
735 bound_N.set_spectral_va().coef_i();
736 bound_N.set_spectral_va().ylm();
737 }
738
739 else{
740 cout << "non_spherical case has not been treated yet!" << endl;}
741
742 // Interpolation to get a 3d value (as poisson solvers take this for boundary...)
743
744
745
746
747 int nr2 = (*lapse.get_mp().get_mg()).get_nr(1);
748 int nt2 = (*lapse.get_mp().get_mg()).get_nt(1);
749 int np2 = (*lapse.get_mp().get_mg()).get_np(1);
750
751
752 for (int f= 0; f<nz2; f++)
753 for (int k=0; k<np2; k++)
754 for (int j=0; j<nt2; j++) {
755 for (int l=0; l<nr2; l++) {
756
757 bound_N3.set_grid_point(f,k,j,l) = bound_N.val_grid_point(0,k,j,0);
758
759
760 }
761 }
762 if (nz2 >2){
763 bound_N3.annule_domain(nz2 - 1);
764
765 }
766
767 bound_N3.std_spectral_base();
768 bound_N3.set_spectral_va().ylm();
769
770 // End interpolation
771
772 p_get_BC_lapse_3 = new Scalar(bound_N3);
773
774 }
775 return *p_get_BC_lapse_3 ;
776}
777
778
779
780
781
782// Used to retrieve d (theta)/dt over only initial data. Probably obsolete in a practical use on CoCoNuT. Uses the same formalism as get_BC_lapse_3().
784if (p_derive_t_expa == 0x0){
785
786 int nz2 = (*lapse.get_mp().get_mg()).get_nzone();
787 // Radial vector for the full 3-metric.
788 Vector sss = gamij.radial_vect();
789 Vector sss_down = sss.up_down(gamij);
790 Scalar bb = contract (shift,0, sss_down,0);
791
792 Scalar Kappa = bb*contract(contract(Kij,0,sss,0),0,sss,0) - contract(sss,0, lapse.derive_cov(gamij),0);
793 Scalar Matter = 8.*M_PI*(bb*Ee);
794 Matter.annule_domain(nz2-1);
795 Matter += - 8.*M_PI*(bb + lapse)*contract(Jj,0,sss_down,0);
796 Matter.annule_domain(nz2-1);
797 Matter += 8.*M_PI*lapse*contract(contract(Sij,0,sss_down,0),0,sss_down,0);
798
799
800 // 2d interpolation of the Kappa constant and the matter terms.
801
802 Scalar Kappa2 = sph.get_ricci();
803 Kappa2.annule_hard();
804 Scalar bb2 = Kappa2;
805 Scalar Matter2 = Kappa2;
806 Scalar nn2 = Kappa2;
807
808 int nt = (*Kappa2.get_mp().get_mg()).get_nt(0);
809 int np = (*Kappa2.get_mp().get_mg()).get_np(0);
810 for (int k=0; k<np; k++)
811 for (int j=0; j<nt; j++) {
812 Kappa2.set_grid_point(0,k,j,0) = Kappa.val_grid_point(1,k,j,0);
813 bb2.set_grid_point(0,k,j,0) = bb.val_grid_point(1,k,j,0);
814 Matter2.set_grid_point(0,k,j,0) = Matter.val_grid_point(1,k,j,0);
815 nn2.set_grid_point(0,k,j,0) = lapse.val_grid_point(1,k,j,0);
816 }
817
818 Kappa2.std_spectral_base();
819 bb2.std_spectral_base();
820 Matter2.std_spectral_base();
821
822 Scalar Tplus = sph.theta_plus();
823 Scalar Tminus = sph.theta_minus();
824 Scalar Ricci2 = sph.get_ricci();
825 Vector Ll = sph.get_ll();
826 Sym_tensor Shear = sph.shear();
827 Metric qab = sph.get_qab();
828
829 Scalar source = 2.*contract(Ll.up_down(qab),0,sph.derive_cov2d(nn2 -bb2),0) + (nn2- bb2)*(contract(sph.derive_cov2d(Ll.up_down(qab)),0,1)+ contract(Ll,0,Ll.up_down(qab),0)- 0.5*(Ricci2 + Tplus*Tminus)) -(bb2 +nn2)*(0.25*Tplus*Tplus + 0.5*contract(Shear.up_down(qab),0,1,Shear,0,1)) - Matter2 - Tplus*Kappa2 ;
830
831 Scalar lapNmb = contract(sph.derive_cov2d(sph.derive_cov2d(nn2-bb2)).up(0,qab),0,1);
832
833 source = source + lapNmb;
834 source.set_spectral_va().ylm();
835
836 Scalar difftheta = source*delta_t;
837 cout << "mean difference between old and new expansion" << endl;
838 cout << difftheta.val_grid_point(0,0,0,0) << endl;
839
840
841 p_derive_t_expa = new Scalar (source);
842 }
843
844 return *p_derive_t_expa ;
845}
846
847
848
849
850
851
852
853
854// Source for the Dirichlet BC on the shift
855const Vector& Excision_surf::get_BC_shift_1(double Omega) const{
856 if (p_get_BC_shift_1 == 0x0){
857
858 // Radial vector for the full 3-metric.
859 Vector sss = gamij.radial_vect();
860
861 // Boundary value for the radial part of the shift
862 Scalar bound = lapse ;
863 // bound = bound + 0.05; // REMOVE real quick
864
865 // Tangent part of the shift
866 // (For now, only axisymmetric configurations are envisaged)
867
868 Vector V_par = shift;
869 Scalar V_phi = lapse; V_phi.annule_hard(); V_phi = 1.; // Rotation parameter for spacetime
870 V_phi.std_spectral_base() ; V_phi.mult_rsint();
871 V_par.set(1).annule_hard();
872 V_par.set(2).annule_hard();
873 V_par.set(3) = V_phi;
874
875 V_par = V_par*Omega;
876
877
878 // Construction of the total shift boundary condition
879 Vector bound_shift = bound*sss + V_par;
880 bound_shift.std_spectral_base(); // Boundary is fixed by value of 3 components of a vector (rather than value of potentials)
881 p_get_BC_shift_1 = new Vector(bound_shift);
882 }
883 return *p_get_BC_shift_1 ;
884}
885
886
887// Source for the Dirichlet BC on the shift, based on a Parabolic driver
888const Vector& Excision_surf::get_BC_shift_2(double c_bb_lap, double c_bb_fin, double c_V_lap, double epsilon) const{
889 if (p_get_BC_shift_2 == 0x0){
890
891 // Radial vector for the full 3-metric.
892 Vector sss = gamij.radial_vect();
893 Vector sss_down = sss.up_down(gamij);
894
895// // Boundary value for the radial part of the shift: parabolic driver for (b-N)
896 // Scalar bound = lapse ;
897 Scalar bb = contract (shift,0, sss_down,0);
898
899 Scalar b_min_N = bb - lapse;
900 Scalar ff = lapse*(c_bb_lap*b_min_N.lapang() + c_bb_fin*b_min_N);
901
902 ff.std_spectral_base();
903
904
905 // Definition of k_1
906 Scalar k_1 =delta_t*ff;
907
908 // Intermediate value of b-N, for Runge-Kutta 2nd order scheme
909 Scalar b_min_N_int = b_min_N + k_1; b_min_N_int.std_spectral_base();
910
911 // Recalculation of ff with intermediate values.
912 Scalar ff_int = lapse*(c_bb_lap*b_min_N_int.lapang() + c_bb_fin*b_min_N_int);
913
914 // Definition of k_2
915 Scalar k_2 = delta_t*ff_int;
916 k_2.std_spectral_base();
917
918 // Result of RK2 evolution
919 Scalar bound_b_min_N = b_min_N + k_2;
920 bound_b_min_N.std_spectral_base();
921 bound_b_min_N.set_spectral_va().ylm();
922
923 Scalar bb2 = bound_b_min_N + lapse; // Look out for additional term for time variation of lapse and sss
924
925
926 // Tangent part of the shift, with parabolic driver
927
928
929 Vector V_par = shift - bb*sss;
930 Sym_tensor q_upup = gamij.con() - sss*sss;
931 Sym_tensor q_downdown = gamij.cov() - sss_down*sss_down;
932 Tensor q_updown = q_upup.down(1, gamij);
933
934
935
936
937
938 // Calculation of the conformal 2d laplacian of V
939 Tensor d_V = contract(q_updown, 1, contract(q_updown,0 , V_par.derive_cov(gamij),1 ),1);
940 Tensor d_V_con = contract(d_V,1,q_upup,1);
941 Tensor dd_V = d_V_con.derive_cov(gamij);
942 // Vector lap_V = contract(dd_V.derive_cov(gamij),1,2);
943 dd_V = contract(q_updown, 1, contract(q_updown,1 ,contract(q_updown,0,dd_V, 2), 2), 2);
944 Tensor dd_Vdown = contract (dd_V,1,q_downdown,1);
945 Vector Ricci1 = contract(dd_Vdown,0,1) - contract(dd_Vdown,0,2);
946 Vector Ricci2 = contract (q_upup,1,Ricci1,0);
947 Vector lap_V = contract(dd_V, 1,2);
948
949
950
951// Tensor dd_V = V_par.derive_con(gamij);
952// dd_V = contract(q_updown, 1, contract(q_updown,1 ,dd_V, 0), 1);
953// Vector lap_V = contract(q_updown, 1, contract(dd_V.derive_cov(gamij),1,2), 0);
954
955
956
957 // 3d interpolation of the Ricci scalar on the surface.
958
959 Scalar ricci2 = sph.get_ricci();
960
961 // Start Mapping interpolation
962
963 Scalar ricci3 (lapse.get_mp());
964
965 ricci3.allocate_all();
966 ricci3.std_spectral_base();
967
968 int nz = (*lapse.get_mp().get_mg()).get_nzone();
969 int nr = (*lapse.get_mp().get_mg()).get_nr(1);
970 int nt = (*lapse.get_mp().get_mg()).get_nt(1);
971 int np = (*lapse.get_mp().get_mg()).get_np(1);
972
973
974 for (int f= 0; f<nz; f++)
975 for (int k=0; k<np; k++)
976 for (int j=0; j<nt; j++) {
977 for (int l=0; l<nr; l++) {
978
979 ricci3.set_grid_point(f,k,j,l) = ricci2.val_grid_point(0,k,j,0);
980
981
982 }
983 }
984 if (nz >2){
985 // ricci3.annule_domain(0);
986 ricci3.annule_domain(nz - 1);
987
988 }
989 // End Mapping interpolation
990
991 // Construction of the Ricci COV tensor on the sphere
992
993 Sym_tensor ricci_t = gamij.cov() - sss_down*sss_down;
994 ricci_t = 0.5*ricci3*ricci_t;
995 ricci_t.std_spectral_base();
996
997 Tensor ricci_t_updown = contract(q_upup,0, ricci_t,0);
998
999 // Calculation of ff
1000
1001 // Vector ffV = c_V_lap*lapse*(lap_V + contract(ricci_t_updown,1, V_par,0));
1002 Vector ffV = c_V_lap*lapse*(lap_V + Ricci2);
1003
1004 ffV.annule_domain(nz-1);
1005
1006
1007// cout << "verification of vanishing shear" << endl;
1008// cout << "points on the surface" << endl;
1009
1010// Tbl val_ffV(3, nt, np);
1011// val_ffV.set_etat_qcq();
1012// for(int mm=0; mm <np; mm++)
1013// for (int pp=0; pp<nt; pp++){
1014// val_ffV.set(0,pp,mm) = ffV(1).val_grid_point(1,mm,pp,0);
1015// val_ffV.set(1,pp,mm) = ffV(2).val_grid_point(1,mm,pp,0);
1016// val_ffV.set(2,pp,mm) = ffV(3).val_grid_point(1,mm,pp,0);
1017// }
1018
1019// // cout << val_ffV<< endl;
1020// cout << max(val_ffV) << endl;
1021// val_ffV.set_etat_nondef();
1022
1023
1024
1025 ffV = ffV + epsilon*lapse*contract(V_par,0,V_par.derive_cov(gamij),1); // Add of dragging term for time transport.
1026 ffV.annule_domain(nz-1);
1027
1028 ffV.std_spectral_base();
1029
1030// cout << "verification of vanishing shear with dragging" << endl;
1031// cout << "points on the surface" << endl;
1032// val_ffV.set_etat_qcq();
1033// for(int mm=0; mm <np; mm++)
1034// for (int pp=0; pp<nt; pp++){
1035// val_ffV.set(0,pp,mm) = ffV(1).val_grid_point(1,mm,pp,0);
1036// val_ffV.set(1,pp,mm) = ffV(2).val_grid_point(1,mm,pp,0);
1037// val_ffV.set(2,pp,mm) = ffV(3).val_grid_point(1,mm,pp,0);
1038// }
1039
1040// // cout << val_ffV<< endl;
1041// cout << max(val_ffV) << endl;
1042// val_ffV.set_etat_nondef();
1043
1044
1045
1046
1047 // Definition of k_1
1048 Vector k_1V =delta_t*ffV;
1049
1050 // Intermediate value of Npsi, for Runge-Kutta 2nd order scheme
1051 if (nz >2){
1052 k_1V.annule_domain(nz-1);
1053 } // Patch to avoid dzpuis problems if existent.
1054 Vector V_par_int = V_par + k_1V;// V_par_int.std_spectral_base();
1055
1056 // Calculation of the conformal 2d laplacian of V
1057 Tensor d_V_int = contract(q_updown, 1, contract(q_updown,0 , V_par_int.derive_cov(gamij),1 ),1);
1058 Tensor d_V_con_int = contract(d_V_int,1,q_upup,1);
1059 Tensor dd_V_int = d_V_con_int.derive_cov(gamij);
1060 // Vector lap_V = contract(dd_V.derive_cov(gamij),1,2);
1061 dd_V_int = contract(q_updown, 1, contract(q_updown,1 ,contract(q_updown,0,dd_V_int, 2), 2), 2);
1062 Tensor dd_Vdown_int = contract (dd_V_int,1,q_downdown,1);
1063 Vector Ricci1_int = contract(dd_Vdown_int,0,1) - contract(dd_Vdown_int,0,2);
1064 Vector Ricci2_int = contract (q_upup,1,Ricci1_int,0);
1065 Vector lap_V_int = contract(dd_V_int, 1,2);
1066
1067
1068 // Recalculation of ff with intermediate values.
1069// Tensor dd_V_int = V_par_int.derive_con(gamij).derive_cov(gamij);
1070// // Vector lap_V = contract(dd_V.derive_cov(gamij),1,2);
1071// dd_V_int = contract(q_updown, 1, contract(q_updown,1 ,contract(q_updown,0,dd_V_int, 2), 2), 2);
1072// Vector lap_V_int = contract(dd_V_int, 1,2);
1073
1074 Vector ffV_int = c_V_lap*lapse*(lap_V_int + contract(ricci_t_updown,1, V_par_int,0));
1075 ffV_int.annule_domain(nz-1);
1076 ffV_int = ffV_int + epsilon*lapse*contract(V_par_int,0,V_par_int.derive_cov(gamij),1); // Add of dragging term for time transport.
1077
1078 // ffV_int = -ffV_int; // Only a test..
1079
1080
1081 ffV_int.annule_domain(nz-1);
1082
1083
1084// cout << "verification of vanishing shear for intermediate RK value" << endl;
1085// cout << "points on the surface" << endl;
1086
1087
1088// val_ffV.set_etat_qcq();
1089// for(int mm=0; mm <np; mm++)
1090// for (int pp=0; pp<nt; pp++){
1091// val_ffV.set(0,pp,mm) = ffV(1).val_grid_point(1,mm,pp,0);
1092// val_ffV.set(1,pp,mm) = ffV(2).val_grid_point(1,mm,pp,0);
1093// val_ffV.set(2,pp,mm) = ffV(3).val_grid_point(1,mm,pp,0);
1094// }
1095
1096// // cout << val_ffV<< endl;
1097// cout << max(val_ffV) << endl;
1098// val_ffV.set_etat_nondef();
1099
1100
1101
1102 // Definition of k_2
1103 Vector k_2V = delta_t*ffV_int;
1104 // k_2.std_spectral_base();
1105
1106 // Result of RK2 evolution
1107 if (nz >2){
1108 k_2V.annule_domain(nz-1);
1109 }
1110 Vector bound_V = V_par + k_2V;
1111 // bound_V.std_spectral_base();
1112
1113 // Construction of the total shift boundary condition
1114 Vector bound_shift = bb2*sss + bound_V;
1115 bound_shift.std_spectral_base();
1116 p_get_BC_shift_2 = new Vector(bound_shift);
1117}
1118 return *p_get_BC_shift_2 ;
1119}
1120
1121
1122
1123
1124
1125// Source for the Dirichlet BC on the shift, based on kinematical relation for the radial part, and a parabolic evolution for the tangent part.
1126// It takes as a fixed argument the time derivative of the psi factor, that can be deduced from the function.get_BC_conf_fact_4().
1127const Vector& Excision_surf::get_BC_shift_3(Scalar& dtpsi, double c_V_lap, double epsilon) const{
1128 if (p_get_BC_shift_3 == 0x0){
1129
1130 // Radial vector for the full 3-metric.
1131 Vector sss = gamij.radial_vect();
1132 Vector sss_down = sss.up_down(gamij);
1133
1134 const Metric_flat& flat = lapse.get_mp().flat_met_spher() ;
1135
1136 int nz = (*lapse.get_mp().get_mg()).get_nzone();
1137
1138// // Boundary value for the radial part of the shift: parabolic driver for (b-N)
1139 // Scalar bound = lapse ;
1140 Scalar bb = contract (shift,0, sss_down,0);
1141
1142 Scalar pure_source = contract(sss,0,bb.derive_cov(flat),0)*conf_fact*1./6.;
1143 pure_source.annule_domain(nz-1); // CAREFUL
1144 pure_source = dtpsi - pure_source;
1145 pure_source.annule_domain(nz-1);
1146
1147 Scalar factor = conf_fact*sss.divergence(flat)*1./6.;
1148 factor.annule_domain(nz-1);
1149 factor = factor + contract (sss, 0, conf_fact.derive_cov(flat),0);
1150
1151 Scalar source = pure_source/factor;
1152
1153
1154
1155 Scalar bb2 = source; // Look out for additional term for time variation of lapse and sss
1156
1157
1158 // Tangent part of the shift, with parabolic driver
1159
1160
1161 Vector V_par = shift - bb*sss;
1162 Sym_tensor q_upup = gamij.con() - sss*sss;
1163 Sym_tensor q_downdown = gamij.cov() - sss_down*sss_down;
1164 Tensor q_updown = q_upup.down(1, gamij);
1165
1166
1167 // Calculation of the conformal 2d laplacian of V
1168 Tensor d_V = contract(q_updown, 1, contract(q_updown,0 , V_par.derive_cov(gamij),1 ),1);
1169 Tensor d_V_con = contract(d_V,1,q_upup,1);
1170 Tensor dd_V = d_V_con.derive_cov(gamij);
1171 // Vector lap_V = contract(dd_V.derive_cov(gamij),1,2);
1172 dd_V = contract(q_updown, 1, contract(q_updown,1 ,contract(q_updown,0,dd_V, 2), 2), 2);
1173 Tensor dd_Vdown = contract (dd_V,1,q_downdown,1);
1174 Vector Ricci1 = contract(dd_Vdown,0,1) - contract(dd_Vdown,0,2);
1175 Vector Ricci2 = contract (q_upup,1,Ricci1,0);
1176 Vector lap_V = contract(dd_V, 1,2);
1177
1178
1179
1180// Tensor dd_V = V_par.derive_con(gamij);
1181// dd_V = contract(q_updown, 1, contract(q_updown,1 ,dd_V, 0), 1);
1182// Vector lap_V = contract(q_updown, 1, contract(dd_V.derive_cov(gamij),1,2), 0);
1183
1184
1185
1186 // 3d interpolation of the Ricci scalar on the surface.
1187
1188 Scalar ricci2 = sph.get_ricci();
1189
1190 // Start Mapping interpolation
1191
1192 Scalar ricci3 (lapse.get_mp());
1193
1194 ricci3.allocate_all();
1195 ricci3.std_spectral_base();
1196
1197
1198 int nr = (*lapse.get_mp().get_mg()).get_nr(1);
1199 int nt = (*lapse.get_mp().get_mg()).get_nt(1);
1200 int np = (*lapse.get_mp().get_mg()).get_np(1);
1201
1202
1203 for (int f= 0; f<nz; f++)
1204 for (int k=0; k<np; k++)
1205 for (int j=0; j<nt; j++) {
1206 for (int l=0; l<nr; l++) {
1207
1208 ricci3.set_grid_point(f,k,j,l) = ricci2.val_grid_point(0,k,j,0);
1209
1210
1211 }
1212 }
1213 if (nz >2){
1214 // ricci3.annule_domain(0);
1215 ricci3.annule_domain(nz - 1);
1216
1217 }
1218 // End Mapping interpolation
1219
1220 // Construction of the Ricci COV tensor on the sphere
1221
1222 Sym_tensor ricci_t = gamij.cov() - sss_down*sss_down;
1223 ricci_t = 0.5*ricci3*ricci_t;
1224 ricci_t.std_spectral_base();
1225
1226 Tensor ricci_t_updown = contract(q_upup,0, ricci_t,0);
1227
1228 // Calculation of ff
1229
1230 // Vector ffV = c_V_lap*lapse*(lap_V + contract(ricci_t_updown,1, V_par,0));
1231 Vector ffV = c_V_lap*lapse*(lap_V + Ricci2);
1232
1233 ffV.annule_domain(nz-1);
1234
1235
1236// cout << "verification of vanishing shear" << endl;
1237// cout << "points on the surface" << endl;
1238
1239// Tbl val_ffV(3, nt, np);
1240// val_ffV.set_etat_qcq();
1241// for(int mm=0; mm <np; mm++)
1242// for (int pp=0; pp<nt; pp++){
1243// val_ffV.set(0,pp,mm) = ffV(1).val_grid_point(1,mm,pp,0);
1244// val_ffV.set(1,pp,mm) = ffV(2).val_grid_point(1,mm,pp,0);
1245// val_ffV.set(2,pp,mm) = ffV(3).val_grid_point(1,mm,pp,0);
1246// }
1247
1248// // cout << val_ffV<< endl;
1249// cout << max(val_ffV) << endl;
1250// val_ffV.set_etat_nondef();
1251
1252
1253
1254 ffV = ffV + epsilon*lapse*contract(V_par,0,V_par.derive_cov(gamij),1); // Add of dragging term for time transport.
1255 ffV.annule_domain(nz-1);
1256
1257 ffV.std_spectral_base();
1258
1259// cout << "verification of vanishing shear with dragging" << endl;
1260// cout << "points on the surface" << endl;
1261// val_ffV.set_etat_qcq();
1262// for(int mm=0; mm <np; mm++)
1263// for (int pp=0; pp<nt; pp++){
1264// val_ffV.set(0,pp,mm) = ffV(1).val_grid_point(1,mm,pp,0);
1265// val_ffV.set(1,pp,mm) = ffV(2).val_grid_point(1,mm,pp,0);
1266// val_ffV.set(2,pp,mm) = ffV(3).val_grid_point(1,mm,pp,0);
1267// }
1268
1269// // cout << val_ffV<< endl;
1270// cout << max(val_ffV) << endl;
1271// val_ffV.set_etat_nondef();
1272
1273
1274
1275
1276 // Definition of k_1
1277 Vector k_1V =delta_t*ffV;
1278
1279 // Intermediate value of Npsi, for Runge-Kutta 2nd order scheme
1280 if (nz >2){
1281 k_1V.annule_domain(nz-1);
1282 } // Patch to avoid dzpuis problems if existent.
1283 Vector V_par_int = V_par + k_1V;// V_par_int.std_spectral_base();
1284
1285 // Calculation of the conformal 2d laplacian of V
1286 Tensor d_V_int = contract(q_updown, 1, contract(q_updown,0 , V_par_int.derive_cov(gamij),1 ),1);
1287 Tensor d_V_con_int = contract(d_V_int,1,q_upup,1);
1288 Tensor dd_V_int = d_V_con_int.derive_cov(gamij);
1289 // Vector lap_V = contract(dd_V.derive_cov(gamij),1,2);
1290 dd_V_int = contract(q_updown, 1, contract(q_updown,1 ,contract(q_updown,0,dd_V_int, 2), 2), 2);
1291 Tensor dd_Vdown_int = contract (dd_V_int,1,q_downdown,1);
1292 Vector Ricci1_int = contract(dd_Vdown_int,0,1) - contract(dd_Vdown_int,0,2);
1293 Vector Ricci2_int = contract (q_upup,1,Ricci1_int,0);
1294 Vector lap_V_int = contract(dd_V_int, 1,2);
1295
1296
1297 // Recalculation of ff with intermediate values.
1298// Tensor dd_V_int = V_par_int.derive_con(gamij).derive_cov(gamij);
1299// // Vector lap_V = contract(dd_V.derive_cov(gamij),1,2);
1300// dd_V_int = contract(q_updown, 1, contract(q_updown,1 ,contract(q_updown,0,dd_V_int, 2), 2), 2);
1301// Vector lap_V_int = contract(dd_V_int, 1,2);
1302
1303 Vector ffV_int = c_V_lap*lapse*(lap_V_int + contract(ricci_t_updown,1, V_par_int,0));
1304 ffV_int.annule_domain(nz-1);
1305 ffV_int = ffV_int + epsilon*lapse*contract(V_par_int,0,V_par_int.derive_cov(gamij),1); // Add of dragging term for time transport.
1306
1307 // ffV_int = -ffV_int; // Only a test..
1308
1309
1310 ffV_int.annule_domain(nz-1);
1311
1312
1313// cout << "verification of vanishing shear for intermediate RK value" << endl;
1314// cout << "points on the surface" << endl;
1315
1316
1317// val_ffV.set_etat_qcq();
1318// for(int mm=0; mm <np; mm++)
1319// for (int pp=0; pp<nt; pp++){
1320// val_ffV.set(0,pp,mm) = ffV(1).val_grid_point(1,mm,pp,0);
1321// val_ffV.set(1,pp,mm) = ffV(2).val_grid_point(1,mm,pp,0);
1322// val_ffV.set(2,pp,mm) = ffV(3).val_grid_point(1,mm,pp,0);
1323// }
1324
1325// // cout << val_ffV<< endl;
1326// cout << max(val_ffV) << endl;
1327// val_ffV.set_etat_nondef();
1328
1329
1330
1331 // Definition of k_2
1332 Vector k_2V = delta_t*ffV_int;
1333 // k_2.std_spectral_base();
1334
1335 // Result of RK2 evolution
1336 if (nz >2){
1337 k_2V.annule_domain(nz-1);
1338 }
1339 Vector bound_V = V_par + k_2V;
1340 // bound_V.std_spectral_base();
1341
1342 // Construction of the total shift boundary condition
1343 Vector bound_shift = bb2*sss + bound_V;
1344 bound_shift.std_spectral_base();
1345 p_get_BC_shift_3 = new Vector(bound_shift);
1346}
1347 return *p_get_BC_shift_3 ;
1348}
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359// Source for the Dirichlet BC on the shift, based on projection of Einstein Equations for the radial part, and a parabolic evolution for the tangent part.
1360// It takes as fixed arguments the time derivative of the expansion, the matter terms and a boolean specifying whether or not we are in spherical symmetry..
1361const Vector& Excision_surf::get_BC_shift_4(Scalar& dttheta, Scalar& Ee, Vector& Jj, Sym_tensor& Sij, double c_V_lap, double epsilon, bool sph_sym) const{
1362 if (p_get_BC_shift_4 == 0x0){
1363
1364 // Radial vector for the full 3-metric.
1365 Vector sss = gamij.radial_vect();
1366 Vector sss_down = sss.up_down(gamij);
1367
1368 int nz = (*lapse.get_mp().get_mg()).get_nzone();
1369
1370// // Boundary value for the radial part of the shift: parabolic driver for (b-N)
1371 // Scalar bound = lapse ;
1372 Scalar bb = contract (shift,0, sss_down,0);
1373
1374
1375
1376 Scalar bound_b3(lapse.get_mp()); bound_b3.allocate_all(); bound_b3.set_spectral_base(bb.get_spectral_base());// Result to be stored there.
1377
1378 Scalar Kappa = bb*contract(contract(Kij,0,sss,0),0,sss,0) - contract(sss,0, lapse.derive_cov(gamij),0);
1379 Scalar Matter = 8.*M_PI*(bb*Ee);
1380 Matter.annule_domain(nz-1);
1381 Matter += - 8.*M_PI*(bb + lapse)*contract(Jj,0,sss_down,0);
1382 Matter.annule_domain(nz-1);
1383 Matter += 8.*M_PI*lapse*contract(contract(Sij,0,sss_down,0),0,sss_down,0);
1384
1385
1386 // 2d interpolation of the Kappa constant and the matter terms.
1387
1388 Scalar Kappa2 = sph.get_ricci();
1389 Kappa2.annule_hard();
1390 Scalar bb2 = Kappa2;
1391 Scalar Matter2 = Kappa2;
1392 Scalar nn2 = Kappa2;
1393
1394 const Metric_flat& flat2 = Kappa2.get_mp().flat_met_spher();
1395
1396 int nt = (*Kappa2.get_mp().get_mg()).get_nt(0);
1397 int np = (*Kappa2.get_mp().get_mg()).get_np(0);
1398 for (int k=0; k<np; k++)
1399 for (int j=0; j<nt; j++) {
1400 Kappa2.set_grid_point(0,k,j,0) = Kappa.val_grid_point(1,k,j,0);
1401 bb2.set_grid_point(0,k,j,0) = bb.val_grid_point(1,k,j,0);
1402 Matter2.set_grid_point(0,k,j,0) = Matter.val_grid_point(1,k,j,0);
1403 nn2.set_grid_point(0,k,j,0) = lapse.val_grid_point(1,k,j,0);
1404 }
1405
1406 Kappa2.std_spectral_base();
1407 bb2.std_spectral_base();
1408 Matter2.std_spectral_base();
1409 nn2.std_spectral_base();
1410
1411
1412 Scalar Tplus = sph.theta_plus();
1413 Scalar Tminus = sph.theta_minus();
1414 Scalar Ricci = sph.get_ricci();
1415 Vector Ll = sph.get_ll();
1416 Sym_tensor Shear = sph.shear();
1417 Metric qab = sph.get_qab();
1418
1419 Scalar source = 2.*contract(Ll.up_down(qab),0,sph.derive_cov2d(nn2),0) + nn2*(contract(sph.derive_cov2d(Ll.up_down(qab)),0,1)+ contract(Ll,0,Ll.up_down(qab),0)- 0.5*(Ricci + Tplus*Tminus) - 0.25*Tplus*Tplus - 0.5*contract(Shear.up_down(qab),0,1,Shear,0,1)) - Matter2 - Tplus*Kappa2 - dttheta;
1420
1421 Scalar lapN = contract(sph.derive_cov2d(sph.derive_cov2d(nn2)).up(0,qab),0,1);
1422
1423 source = source + lapN;
1424
1425
1426 Scalar source_add = 2.*contract(Ll.up_down(qab),0,sph.derive_cov2d(bb2),0);
1427
1428 source = source - source_add;
1429
1430 Scalar sqrt_q_h2 = sph.sqrt_q() * sph.get_hsurf() * sph.get_hsurf() ;
1431 Tensor Delta2 = sph.delta();
1432 Sym_tensor qab_con = sph.get_qab().con() / (sph.get_hsurf() * sph.get_hsurf()) ;
1433 qab_con.set(1,1) = 1. ;
1434 qab_con.set(1,2) = 0. ;
1435 qab_con.set(1,3) = 0. ;
1436 qab_con.std_spectral_base() ;
1437
1438 Sym_tensor hab =(qab_con*sqrt_q_h2 - flat2.con()) / (sph.get_hsurf()*sph.get_hsurf()) ;
1439 hab.set(1,1) = 1. ;
1440 hab.set(1,2) = 0. ;
1441 hab.set(1,3) = 0. ;
1442 hab.std_spectral_base() ;
1443
1444 Vector db = sph.derive_cov2dflat(bb2);
1445 Tensor ddb = sph.derive_cov2dflat(db);
1446
1447 Scalar lap_rem = contract(qab_con, 0,1, contract(Delta2,0,db,0),0,1)*sqrt_q_h2 - sph.get_hsurf()*sph.get_hsurf()*contract(hab,0,1,ddb,0,1);// What remains of the laplacian
1448
1449 source = sqrt_q_h2*source + lap_rem;
1450
1451 Scalar bound_b=bb2;
1452 Scalar lapang_s_par = (contract(sph.derive_cov2d(Ll.up_down(qab)),0,1)+ contract(Ll,0,Ll.up_down(qab),0)- 0.5*(Ricci + Tplus*Tminus) + 0.25*Tplus*Tplus + 0.5*contract(Shear.up_down(qab),0,1,Shear,0,1))*sqrt_q_h2;
1453 double lapang_par = lapang_s_par.val_grid_point(0,0,0,0);
1454
1455 if (sph_sym == true){
1456 bound_b = source/lapang_par;
1457 bound_b.set_spectral_va().coef_i();
1458 bound_b.set_spectral_va().ylm();
1459 }
1460
1461 else{
1462 cout << "non spherical case has not been treated yet!" << endl;
1463 }
1464
1465 // Interpolation to get a 3d value (as poisson solvers take this for boundary...)
1466
1467
1468
1469
1470 int nr2 = (*lapse.get_mp().get_mg()).get_nr(1);
1471 int nt2 = (*lapse.get_mp().get_mg()).get_nt(1);
1472 int np2 = (*lapse.get_mp().get_mg()).get_np(1);
1473
1474
1475 for (int f= 0; f<nz; f++)
1476 for (int k=0; k<np2; k++)
1477 for (int j=0; j<nt2; j++) {
1478 for (int l=0; l<nr2; l++) {
1479
1480 bound_b3.set_grid_point(f,k,j,l) = bound_b.val_grid_point(0,k,j,0);
1481
1482
1483 }
1484 }
1485 if (nz >2){
1486 bound_b3.annule_domain(nz - 1);
1487
1488 }
1489
1490 bound_b3.std_spectral_base();
1491 bound_b3.set_spectral_va().ylm();
1492
1493
1494
1495
1496 Scalar bbb2 = bound_b3; // Look out for additional term for time variation of lapse and sss
1497
1498
1499 // Tangent part of the shift, with parabolic driver
1500
1501
1502 Vector V_par = shift - bb*sss;
1503 Sym_tensor q_upup = gamij.con() - sss*sss;
1504 Sym_tensor q_downdown = gamij.cov() - sss_down*sss_down;
1505 Tensor q_updown = q_upup.down(1, gamij);
1506
1507
1508 // Calculation of the conformal 2d laplacian of V
1509 Tensor d_V = contract(q_updown, 1, contract(q_updown,0 , V_par.derive_cov(gamij),1 ),1);
1510 Tensor d_V_con = contract(d_V,1,q_upup,1);
1511 Tensor dd_V = d_V_con.derive_cov(gamij);
1512 // Vector lap_V = contract(dd_V.derive_cov(gamij),1,2);
1513 dd_V = contract(q_updown, 1, contract(q_updown,1 ,contract(q_updown,0,dd_V, 2), 2), 2);
1514 Tensor dd_Vdown = contract (dd_V,1,q_downdown,1);
1515 Vector Ricci1 = contract(dd_Vdown,0,1) - contract(dd_Vdown,0,2);
1516 Vector Ricci2 = contract (q_upup,1,Ricci1,0);
1517 Vector lap_V = contract(dd_V, 1,2);
1518
1519
1520
1521// Tensor dd_V = V_par.derive_con(gamij);
1522// dd_V = contract(q_updown, 1, contract(q_updown,1 ,dd_V, 0), 1);
1523// Vector lap_V = contract(q_updown, 1, contract(dd_V.derive_cov(gamij),1,2), 0);
1524
1525
1526
1527 // 3d interpolation of the Ricci scalar on the surface.
1528
1529 Scalar ricci2 = sph.get_ricci();
1530
1531 // Start Mapping interpolation
1532
1533 Scalar ricci3 (lapse.get_mp());
1534
1535 ricci3.allocate_all();
1536 ricci3.std_spectral_base();
1537
1538 for (int f= 0; f<nz; f++)
1539 for (int k=0; k<np2; k++)
1540 for (int j=0; j<nt2; j++) {
1541 for (int l=0; l<nr2; l++) {
1542
1543 ricci3.set_grid_point(f,k,j,l) = ricci2.val_grid_point(0,k,j,0);
1544
1545
1546 }
1547 }
1548 if (nz >2){
1549 // ricci3.annule_domain(0);
1550 ricci3.annule_domain(nz - 1);
1551
1552 }
1553 // End Mapping interpolation
1554
1555 // Construction of the Ricci COV tensor on the sphere
1556
1557 Sym_tensor ricci_t = gamij.cov() - sss_down*sss_down;
1558 ricci_t = 0.5*ricci3*ricci_t;
1559 ricci_t.std_spectral_base();
1560
1561 Tensor ricci_t_updown = contract(q_upup,0, ricci_t,0);
1562
1563 // Calculation of ff
1564
1565 // Vector ffV = c_V_lap*lapse*(lap_V + contract(ricci_t_updown,1, V_par,0));
1566 Vector ffV = c_V_lap*lapse*(lap_V + Ricci2);
1567
1568 ffV.annule_domain(nz-1);
1569
1570
1571// cout << "verification of vanishing shear" << endl;
1572// cout << "points on the surface" << endl;
1573
1574// Tbl val_ffV(3, nt, np);
1575// val_ffV.set_etat_qcq();
1576// for(int mm=0; mm <np; mm++)
1577// for (int pp=0; pp<nt; pp++){
1578// val_ffV.set(0,pp,mm) = ffV(1).val_grid_point(1,mm,pp,0);
1579// val_ffV.set(1,pp,mm) = ffV(2).val_grid_point(1,mm,pp,0);
1580// val_ffV.set(2,pp,mm) = ffV(3).val_grid_point(1,mm,pp,0);
1581// }
1582
1583// // cout << val_ffV<< endl;
1584// cout << max(val_ffV) << endl;
1585// val_ffV.set_etat_nondef();
1586
1587
1588
1589 ffV = ffV + epsilon*lapse*contract(V_par,0,V_par.derive_cov(gamij),1); // Add of dragging term for time transport.
1590 ffV.annule_domain(nz-1);
1591
1592 ffV.std_spectral_base();
1593
1594// cout << "verification of vanishing shear with dragging" << endl;
1595// cout << "points on the surface" << endl;
1596// val_ffV.set_etat_qcq();
1597// for(int mm=0; mm <np; mm++)
1598// for (int pp=0; pp<nt; pp++){
1599// val_ffV.set(0,pp,mm) = ffV(1).val_grid_point(1,mm,pp,0);
1600// val_ffV.set(1,pp,mm) = ffV(2).val_grid_point(1,mm,pp,0);
1601// val_ffV.set(2,pp,mm) = ffV(3).val_grid_point(1,mm,pp,0);
1602// }
1603
1604// // cout << val_ffV<< endl;
1605// cout << max(val_ffV) << endl;
1606// val_ffV.set_etat_nondef();
1607
1608
1609
1610
1611 // Definition of k_1
1612 Vector k_1V =delta_t*ffV;
1613
1614 // Intermediate value of Npsi, for Runge-Kutta 2nd order scheme
1615 if (nz >2){
1616 k_1V.annule_domain(nz-1);
1617 } // Patch to avoid dzpuis problems if existent.
1618 Vector V_par_int = V_par + k_1V;// V_par_int.std_spectral_base();
1619
1620 // Calculation of the conformal 2d laplacian of V
1621 Tensor d_V_int = contract(q_updown, 1, contract(q_updown,0 , V_par_int.derive_cov(gamij),1 ),1);
1622 Tensor d_V_con_int = contract(d_V_int,1,q_upup,1);
1623 Tensor dd_V_int = d_V_con_int.derive_cov(gamij);
1624 // Vector lap_V = contract(dd_V.derive_cov(gamij),1,2);
1625 dd_V_int = contract(q_updown, 1, contract(q_updown,1 ,contract(q_updown,0,dd_V_int, 2), 2), 2);
1626 Tensor dd_Vdown_int = contract (dd_V_int,1,q_downdown,1);
1627 Vector Ricci1_int = contract(dd_Vdown_int,0,1) - contract(dd_Vdown_int,0,2);
1628 Vector Ricci2_int = contract (q_upup,1,Ricci1_int,0);
1629 Vector lap_V_int = contract(dd_V_int, 1,2);
1630
1631
1632 // Recalculation of ff with intermediate values.
1633// Tensor dd_V_int = V_par_int.derive_con(gamij).derive_cov(gamij);
1634// // Vector lap_V = contract(dd_V.derive_cov(gamij),1,2);
1635// dd_V_int = contract(q_updown, 1, contract(q_updown,1 ,contract(q_updown,0,dd_V_int, 2), 2), 2);
1636// Vector lap_V_int = contract(dd_V_int, 1,2);
1637
1638 Vector ffV_int = c_V_lap*lapse*(lap_V_int + contract(ricci_t_updown,1, V_par_int,0));
1639 ffV_int.annule_domain(nz-1);
1640 ffV_int = ffV_int + epsilon*lapse*contract(V_par_int,0,V_par_int.derive_cov(gamij),1); // Add of dragging term for time transport.
1641
1642 // ffV_int = -ffV_int; // Only a test..
1643
1644
1645 ffV_int.annule_domain(nz-1);
1646
1647
1648// cout << "verification of vanishing shear for intermediate RK value" << endl;
1649// cout << "points on the surface" << endl;
1650
1651
1652// val_ffV.set_etat_qcq();
1653// for(int mm=0; mm <np; mm++)
1654// for (int pp=0; pp<nt; pp++){
1655// val_ffV.set(0,pp,mm) = ffV(1).val_grid_point(1,mm,pp,0);
1656// val_ffV.set(1,pp,mm) = ffV(2).val_grid_point(1,mm,pp,0);
1657// val_ffV.set(2,pp,mm) = ffV(3).val_grid_point(1,mm,pp,0);
1658// }
1659
1660// // cout << val_ffV<< endl;
1661// cout << max(val_ffV) << endl;
1662// val_ffV.set_etat_nondef();
1663
1664
1665
1666 // Definition of k_2
1667 Vector k_2V = delta_t*ffV_int;
1668 // k_2.std_spectral_base();
1669
1670 // Result of RK2 evolution
1671 if (nz >2){
1672 k_2V.annule_domain(nz-1);
1673 }
1674 Vector bound_V = V_par + k_2V;
1675 // bound_V.std_spectral_base();
1676
1677 // Construction of the total shift boundary condition
1678 Vector bound_shift = bbb2*sss + bound_V;
1679 bound_shift.std_spectral_base();
1680 p_get_BC_shift_4 = new Vector(bound_shift);
1681}
1682 return *p_get_BC_shift_4 ;
1683}
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704// Source for the Dirichlet BC on (N*Psi1)
1705const Scalar& Excision_surf::get_BC_Npsi_1(double value) const{
1706 if (p_get_BC_Npsi_1 == 0x0){
1707
1708 Scalar bound_Npsi = value*conf_fact;
1709 bound_Npsi.set_spectral_va().ylm();
1710 p_get_BC_Npsi_1 = new Scalar(bound_Npsi);
1711
1712}
1713 return *p_get_BC_Npsi_1 ;
1714}
1715
1716// Source for the Dirichlet BC on (N*Psi1), based on a parabolic driver.
1717const Scalar& Excision_surf::get_BC_Npsi_2(double npsi_fin, double c_npsi_lap, double c_npsi_fin) const{
1718 if (p_get_BC_Npsi_2 == 0x0){
1719
1720
1721 Scalar npsi = lapse*conf_fact; npsi.std_spectral_base();
1722 Scalar ff = lapse*(c_npsi_lap*npsi.lapang() + c_npsi_fin*npsi);
1723 ff.std_spectral_base();
1724 ff = ff -lapse*c_npsi_fin*npsi_fin;
1725 ff.std_spectral_base();
1726
1727
1728 // Definition of k_1
1729 Scalar k_1 =delta_t*ff;
1730
1731 // Intermediate value of Npsi, for Runge-Kutta 2nd order scheme
1732 Scalar npsi_int = npsi + k_1; npsi_int.std_spectral_base();
1733
1734 // Recalculation of ff with intermediate values.
1735 Scalar ff_int = lapse*(c_npsi_lap*npsi_int.lapang() + c_npsi_fin*(npsi_int - npsi_fin));
1736
1737 // Definition of k_2
1738 Scalar k_2 = delta_t*ff_int;
1739 k_2.std_spectral_base();
1740
1741 // Result of RK2 evolution
1742 Scalar bound_npsi = npsi + k_2;
1743 bound_npsi.std_spectral_base();
1744 bound_npsi.set_spectral_va().ylm();
1745
1746
1747
1748// Scalar bound_Npsi = value*conf_fact;
1749// bound_Npsi.set_spectral_va().ylm();
1750 p_get_BC_Npsi_2 = new Scalar(bound_npsi);
1751
1752}
1753 return *p_get_BC_Npsi_2 ;
1754}
1755
1756// Source for the Dirichlet BC on (N*Psi1), with a Kerr_Schild-like form for the lapse.
1757const Scalar& Excision_surf::get_BC_Npsi_3(double n_0, double beta) const{
1758 if (p_get_BC_Npsi_3 == 0x0){
1759
1760
1761 const Coord& ct = lapse.get_mp().cost;
1762 Scalar boundN(lapse.get_mp()); boundN = sqrt(n_0 + beta*ct*ct); // Kerr_Schild form for the lapse.
1763 boundN.std_spectral_base();
1764 Scalar bound_npsi = boundN*conf_fact;
1765 bound_npsi.set_spectral_va().ylm();
1766
1767// Scalar npsi = lapse*conf_fact; npsi.std_spectral_base();
1768// Scalar ff = lapse*(c_npsi_lap*npsi.lapang() + c_npsi_fin*npsi);
1769// ff.std_spectral_base();
1770// ff += -lapse*c_npsi_fin*npsi_fin;
1771// ff.std_spectral_base();
1772
1773
1774// // Definition of k_1
1775// Scalar k_1 =delta_t*ff;
1776
1777// // Intermediate value of Npsi, for Runge-Kutta 2nd order scheme
1778// Scalar npsi_int = npsi + k_1; npsi_int.std_spectral_base();
1779
1780// // Recalculation of ff with intermediate values.
1781// Scalar ff_int = lapse*(c_npsi_lap*npsi_int.lapang() + c_npsi_fin*(npsi_int - npsi_fin));
1782
1783// // Definition of k_2
1784// Scalar k_2 = delta_t*ff_int;
1785// k_2.std_spectral_base();
1786
1787// // Result of RK2 evolution
1788// Scalar bound_npsi = npsi + k_2;
1789// bound_npsi.std_spectral_base();
1790// bound_npsi.set_spectral_va().ylm();
1791
1792
1793
1794// Scalar bound_Npsi = value*conf_fact;
1795// bound_Npsi.set_spectral_va().ylm();
1796 p_get_BC_Npsi_3 = new Scalar(bound_npsi);
1797
1798}
1799 return *p_get_BC_Npsi_3 ;
1800}
1801
1802
1803// Source for a Dirichlet BC on (N*Psi1), fixing the surface gravity as a constant in space and time.
1804const Scalar& Excision_surf::get_BC_Npsi_4(double Kappa) const{
1805 if (p_get_BC_Npsi_4 == 0x0){
1806
1807
1808 const Vector s_i = gamij.radial_vect();
1809 Scalar bound_npsi = contract(s_i,0, lapse.derive_cov(gamij),0); bound_npsi.set_dzpuis(0);
1810 bound_npsi = bound_npsi - Kappa;
1811 bound_npsi.std_spectral_base();
1812 bound_npsi = bound_npsi*(conf_fact/contract(contract(Kij,0,s_i,0),0,s_i,0));
1813 bound_npsi.set_dzpuis(0); bound_npsi.std_spectral_base();
1814 bound_npsi.set_spectral_va().ylm();
1815
1816 p_get_BC_Npsi_4 = new Scalar(bound_npsi);
1817
1818}
1819 return *p_get_BC_Npsi_4 ;
1820}
1821
1822
1823
1824// Source for a Neumann BC on (N*Psi1), fixing the surface gravity as a constant in space and time. // Check redundancy with the previous one.
1825const Scalar& Excision_surf::get_BC_Npsi_5(double Kappa) const{
1826 if (p_get_BC_Npsi_5 == 0x0){
1827
1828
1829 const Vector s_i = gamij.radial_vect();
1830 int nz = (*lapse.get_mp().get_mg()).get_nzone();
1831 Scalar bound_npsi = lapse*contract(s_i,0, conf_fact.derive_cov(gamij),0); bound_npsi.annule_domain(nz -1);
1832 Scalar bound_npsi2 = pow(conf_fact,3)*lapse*contract(contract(Kij,0,s_i,0),0,s_i,0);
1833 bound_npsi2.annule_domain(nz-1);
1834 bound_npsi += bound_npsi2;
1835 bound_npsi = bound_npsi + pow(conf_fact,3)*(Kappa); bound_npsi.annule_domain(nz -1);
1836 bound_npsi = bound_npsi -(conf_fact*conf_fact*contract(s_i,0, (conf_fact*lapse).derive_cov(gamij),0) - (conf_fact*lapse).dsdr());
1837 // bound_npsi = bound_npsi*(conf_fact/contract(contract(Kij,0,s_i,0),0,s_i,0));
1838 bound_npsi.annule_domain(nz-1); bound_npsi.std_spectral_base();
1839 bound_npsi.set_spectral_va().ylm();
1840
1841 p_get_BC_Npsi_5 = new Scalar(bound_npsi);
1842
1843}
1844 return *p_get_BC_Npsi_5 ;
1845}
1846
1847
1848
1849
1850
1851 void Excision_surf::sauve(FILE* ) const {
1852
1853 cout << "c'est pas fait!" << endl ;
1854 return ;
1855
1856 }
1857}
Active physical coordinates and mapping derivatives.
Definition coord.h:90
Surface where boundary conditions for quantities in the bulk will be calculated It relies on geometri...
Scalar * p_get_BC_Npsi_5
Source of Neumann boundary condition on .
Scalar * p_get_BC_lapse_3
Source of Dirichlet condtion on , based on einstein equations.
Scalar * p_get_BC_Npsi_3
Source of Dirichlet boundary condition on .
Sym_tensor Kij
The 3-d extrinsic curvature on the slice.
Scalar * p_get_BC_Npsi_4
Source of Dirichlet boundary condition on .
void operator=(const Excision_surf &)
Assignment to another Excision_surf.
Scalar & set_expa()
Sets the expansion function on the surface at time t (considering to protect this function)
const Vector & get_BC_shift_3(Scalar &dtpsi, double c_V_lap, double epsilon) const
Source for a Dirichlet BC on the shift, based on a Parabolic driver; Radial part is dealt with using ...
virtual void sauve(FILE *) const
Save in a file.
const Scalar & get_BC_conf_fact_2(double c_psi_lap, double c_psi_fin, Scalar &expa_fin) const
Source for the Dirichlet BC on the conformal factor, based on a parabolic driver for the conformal fa...
Scalar * p_get_BC_lapse_1
Source of Dirichlet boundary condition of .
Scalar * p_get_BC_conf_fact_1
Source of Neumann boundary condition on ,.
const Vector & get_BC_shift_2(double c_bb_lap, double c_bb_fin, double c_V_lap, double epsilon) const
Source for a Dirichlet BC on the shift, based on a Parabolic driver; no assumptions are made except a...
Scalar lapse
The lapse defined on the 3 slice.
Scalar * p_get_BC_lapse_4
Source of Dirichlet condtion on , based on einstein equations (conservation of isotropic gauge)
Vector shift
The Shift 3-vector on the slice.
const Scalar & get_BC_lapse_3(Scalar &dttheta, Scalar &Ee, Vector &Jj, Sym_tensor &Sij, bool sph_sym=true) const
Source for Dirichlet BC on the lapse, based on einstein equations.
void get_evol_params_from_ID(double alpha, double beta, double gamma, Scalar &Ee, Vector &Jj, Sym_tensor &Ss)
Computes the parameters for the hyperbolic evolution in set_expa_hyperb(), so that the expansion has ...
const Vector & get_BC_shift_4(Scalar &dttheta, Scalar &Ee, Vector &Jj, Sym_tensor &Sij, double c_V_lap, double epsilon, bool sph_sym=true) const
Source for a Dirichlet BC on the shift, based on a Parabolic driver; Radial part is dealt with using ...
virtual ~Excision_surf()
Destructor.
const Scalar & get_BC_conf_fact_1(bool isMOTS=false) const
Source for a Neumann BC on the conformal factor. If boolean isMOTS is false, it is based on expansion...
Scalar * p_get_BC_conf_fact_3
Source of Neumann boundary condition on ,.
Scalar * p_get_BC_Npsi_2
Source of Dirichlet boundary condition on .
const Scalar & get_BC_conf_fact_4() const
Source for the Dirchlet BC on the conformal factor, based on the consistency condition derived from t...
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Scalar * p_get_BC_conf_fact_2
Source of Neumann boundary condition on ,
Scalar & set_dt_expa()
Sets the time derivative of the expansion function on the surface at time t (considering to protect t...
Vector * p_get_BC_shift_4
Source of Dirichlet BC for the shift vector , partly from projection of Einstein Equations.
const Scalar & get_BC_conf_fact_3(double c_theta_lap, double c_theta_fin, Scalar &expa_fin) const
Source for the Neumann BC on the conformal factor, based on a parabolic driver for the expansio.
const Scalar & get_BC_Npsi_5(double Kappa) const
Source for a Neumann BC on (N*Psi1), fixing a constant surface gravity in space and time.
Scalar * p_get_BC_lapse_2
Source of Dirichlet boundary condition of .
const Scalar & derive_t_expa(Scalar &Ee, Vector &Jj, Sym_tensor &Sij) const
Forms the prospective time derivative for the expansion using projected Einstein equations....
Spheroid sph
The associated Spheroid object.
Vector * p_get_BC_shift_2
Source of Dirichlet BC for the shift vector .
const Scalar & get_BC_Npsi_3(double n_0, double beta) const
Source for the Dirichlet BC on (N*Psi1), with Kerr_Schild-like form for the lapse boundary.
Vector * p_get_BC_shift_3
Source of Dirichlet BC for the shift vector , partly derived from kinematical relation.
virtual void del_deriv() const
Deletes all the derived quantities.
const Scalar & get_BC_Npsi_4(double Kappa) const
Source for a Dirichlet BC on (N*Psi1), fixing a constant surface gravity in space and time.
Scalar * p_get_BC_Npsi_1
Source of Neumann boundary condition on .
Scalar conf_fact
The value of the conformal factor on the 3-slice.
double no_of_steps
The internal number of timesteps for one iteration.
Metric gamij
The 3-d metric on the slice.
Scalar * p_derive_t_expa
Computation of an updated expansion scalar.
Excision_surf(const Scalar &h_in, const Metric &gij, const Sym_tensor &Kij2, const Scalar &ppsi, const Scalar &nn, const Vector &beta, double timestep, int int_nos)
Constructor of an excision surface embedded in a 3-slice (Time_slice ) of 3+1 formalism.
const Scalar & get_BC_Npsi_2(double value, double c_npsi_lap, double c_npsi_fin) const
Source for the Dirichlet BC on (N*Psi1), based on a parabolic driver.
const Scalar & get_BC_lapse_2(double lapse_fin, double c_lapse_lap, double c_lapse_fi) const
Source for Dirichlet BC on the lapse, based on a parabolic driver towards arbitrary constant value.
double delta_t
The time step for evolution in parabolic drivers.
Scalar expa
The 2d expansion, directly evolved from the initial excision with Einstein Equations.
Vector * p_get_BC_shift_1
Source of Dirichlet BC for the shift vector .
Scalar * p_get_BC_conf_fact_4
Source of Birichlet boundary condition on ,.
Scalar dt_expa
The time derivative of the expansion, derived from Einstein equations and arbitrary evolution.
Coord cost
Definition map.h:722
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
const Metric_flat & flat_met_spher() const
Returns the flat metric associated with the spherical coordinates and with components expressed in th...
Definition map.C:321
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 Vector & radial_vect() const
Returns the radial vector normal to a spherical slicing and pointing toward spatial infinity.
Definition metric.C:362
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition metric.C:280
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
double val_in_bound_jk(int l, int j, int k) const
Computes the angular coefficient of index j,k of the field represented by *this at by means of the s...
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
const Scalar & lapang() const
Returns the angular Laplacian of *this , where .
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Scalar poisson_angu(double lambda=0) const
Solves the (generalized) angular Poisson equation with *this as source.
Definition scalar_pde.C:200
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 set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
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
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
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va
Definition scalar.h:1294
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
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
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
Definition scalar.C:797
Tensor derive_cov2dflat(const Tensor &uu) const
Computes the round covariant derivative on the spheroid.
Definition spheroid.C:954
const Scalar & sqrt_q() const
Computes the normal vector field to the 2-surface.
Definition spheroid.C:693
const Scalar & theta_minus() const
Computes the ingoing null expansion .
Definition spheroid.C:909
const Metric & get_qab() const
Returns the metric .
Definition spheroid.h:226
const Tensor & delta() const
Computes the delta coefficients for covariant derivative.
Definition spheroid.C:1203
const Vector & get_ll() const
Returns the vector .
Definition spheroid.h:247
const Scalar & get_hsurf() const
Returns the field h_surf.
Definition spheroid.h:223
const Scalar & get_ricci() const
Returns the 2-ricci scalar .
Definition spheroid.h:229
const Sym_tensor & shear() const
Computes the shear of the 2-surface .
Definition spheroid.C:925
Tensor derive_cov2d(const Tensor &uu) const
Computes the total covariant derivative on the spheroid.
Definition spheroid.C:1240
const Scalar & theta_plus() const
Computes the outgoing null expansion .
Definition spheroid.C:889
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
Tensor handling.
Definition tensor.h:288
int get_etat() const
Returns the logical state.
Definition valeur.h:726
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
void coef_i() const
Computes the physical value of *this.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
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
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 pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
Tensor down(int ind, const Metric &gam) const
Computes a new tensor by lowering an index of *this.
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 .
Lorene prototypes.
Definition app_hor.h:64