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 
30 
31 
32 /*
33  * $Header: /cvsroot/Lorene/C++/Source/App_hor/excision_surf.C,v 1.7 2016/12/05 16:17:44 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 
63 namespace Lorene {
64 Excision_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 
90 Excision_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 
132 {
133  del_deriv() ;
134 }
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 
162  p_get_BC_conf_fact_1 = 0x0 ;
163  p_get_BC_lapse_1 = 0x0 ;
164  p_get_BC_shift_1 = 0x0 ;
165  p_get_BC_Npsi_1 = 0x0 ;
166  p_get_BC_conf_fact_2 = 0x0 ;
167  p_get_BC_conf_fact_3 = 0x0 ;
168  p_get_BC_conf_fact_4 = 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 
195 void 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
280 const Scalar& Excision_surf::get_BC_conf_fact_1(bool isMOTS) const{
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
364 const 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
452 const 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
590 const 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
603 const 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.
647 const 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().
784 if (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
855 const 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
888 const 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().
1127 const 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..
1361 const 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)
1705 const 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.
1717 const 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.
1757 const 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.
1804 const 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.
1825 const 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 }
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
Definition: scalar.h:1328
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:675
const Scalar & theta_plus() const
Computes the outgoing null expansion .
Definition: spheroid.C:892
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_4
Source of Dirichlet condtion on , based on einstein equations (conservation of isotropic gauge) ...
Definition: excision_surf.h:93
Metric for tensor calculation.
Definition: metric.h:90
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...
double delta_t
The time step for evolution in parabolic drivers.
Definition: excision_surf.h:68
Sym_tensor Kij
The 3-d extrinsic curvature on the slice.
Definition: excision_surf.h:65
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:293
const Scalar & sqrt_q() const
Computes the normal vector field to the 2-surface.
Definition: spheroid.C:696
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:312
Scalar * p_get_BC_Npsi_2
Source of Dirichlet boundary condition on .
Definition: excision_surf.h:98
Scalar * p_get_BC_Npsi_4
Source of Dirichlet boundary condition on .
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. Does NOT modify the member dt_expa: do it by hand!
Scalar expa
The 2d expansion, directly evolved from the initial excision with Einstein Equations.
Definition: excision_surf.h:74
const Scalar & lapang() const
Returns the angular Laplacian of *this , where .
Definition: scalar_deriv.C:461
Vector * p_get_BC_shift_4
Source of Dirichlet BC for the shift vector , partly from projection of Einstein Equations.
Definition: excision_surf.h:97
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:330
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 ...
Scalar * p_get_BC_conf_fact_2
Source of Neumann boundary condition on ,.
Definition: excision_surf.h:88
Lorene prototypes.
Definition: app_hor.h:67
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...
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:141
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
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...
Flat metric for tensor calculation.
Definition: metric.h:261
Scalar & set_expa()
Sets the expansion function on the surface at time t (considering to protect this function) ...
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
void operator=(const Excision_surf &)
Assignment to another Excision_surf.
Spheroid sph
The associated Spheroid object.
Definition: excision_surf.h:50
void coef_i() const
Computes the physical value of *this.
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
virtual const Vector & radial_vect() const
Returns the radial vector normal to a spherical slicing and pointing toward spatial infinity...
Definition: metric.C:365
const Sym_tensor & shear() const
Computes the shear of the 2-surface .
Definition: spheroid.C:928
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:790
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition: scalar.C:371
Tensor field of valence 1.
Definition: vector.h:188
Scalar poisson_angu(double lambda=0) const
Solves the (generalized) angular Poisson equation with *this as source.
Definition: scalar_pde.C:203
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Scalar & set_dt_expa()
Sets the time derivative of the expansion function on the surface at time t (considering to protect t...
Scalar * p_get_BC_lapse_3
Source of Dirichlet condtion on , based on einstein equations.
Definition: excision_surf.h:92
Scalar * p_get_BC_Npsi_5
Source of Neumann boundary condition on .
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:814
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:643
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: vector.C:322
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:386
Scalar lapse
The lapse defined on the 3 slice.
Definition: excision_surf.h:56
int get_etat() const
Returns the logical state.
Definition: valeur.h:760
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...
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.
Definition: excision_surf.C:64
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric_flat.C:156
Scalar * p_get_BC_conf_fact_4
Source of Birichlet boundary condition on ,.
Definition: excision_surf.h:90
Scalar * p_get_BC_conf_fact_3
Source of Neumann boundary condition on ,.
Definition: excision_surf.h:89
virtual void del_deriv() const
Deletes all the derived quantities.
Tensor derive_cov2dflat(const Tensor &uu) const
Computes the round covariant derivative on the spheroid.
Definition: spheroid.C:957
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 ...
Scalar * p_get_BC_Npsi_3
Source of Dirichlet boundary condition on .
Definition: excision_surf.h:99
Tensor derive_cov2d(const Tensor &uu) const
Computes the total covariant derivative on the spheroid.
Definition: spheroid.C:1243
Scalar * p_derive_t_expa
Computation of an updated expansion scalar.
Definition: excision_surf.h:94
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:438
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:387
virtual void sauve(FILE *) const
Save in a file.
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_ricci() const
Returns the 2-ricci scalar .
Definition: spheroid.h:229
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
Vector * p_get_BC_shift_3
Source of Dirichlet BC for the shift vector , partly derived from kinematical relation.
Definition: excision_surf.h:96
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Tensor handling.
Definition: tensor.h:294
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Definition: tensor.C:1011
Active physical coordinates and mapping derivatives.
Definition: coord.h:90
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.
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:283
Scalar * p_get_BC_lapse_1
Source of Dirichlet boundary condition of .
Definition: excision_surf.h:85
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:690
double no_of_steps
The internal number of timesteps for one iteration.
Definition: excision_surf.h:71
Surface where boundary conditions for quantities in the bulk will be calculated It relies on geometri...
Definition: excision_surf.h:43
const Scalar & theta_minus() const
Computes the ingoing null expansion .
Definition: spheroid.C:912
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:663
virtual ~Excision_surf()
Destructor.
Scalar dt_expa
The time derivative of the expansion, derived from Einstein equations and arbitrary evolution...
Definition: excision_surf.h:77
Scalar conf_fact
The value of the conformal factor on the 3-slice.
Definition: excision_surf.h:53
const Vector & get_ll() const
Returns the vector .
Definition: spheroid.h:247
const Metric & get_qab() const
Returns the metric .
Definition: spheroid.h:226
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:196
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...
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
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 * p_get_BC_lapse_2
Source of Dirichlet boundary condition of .
Definition: excision_surf.h:91
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...
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:610
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:302
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...
Vector * p_get_BC_shift_1
Source of Dirichlet BC for the shift vector .
Definition: excision_surf.h:86
Metric gamij
The 3-d metric on the slice.
Definition: excision_surf.h:62
Vector * p_get_BC_shift_2
Source of Dirichlet BC for the shift vector .
Definition: excision_surf.h:95
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:874
Scalar * p_get_BC_Npsi_1
Source of Neumann boundary condition on .
Definition: excision_surf.h:87
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 void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: tensor.C:935
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:324
Vector shift
The Shift 3-vector on the slice.
Definition: excision_surf.h:59
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:226
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...
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:607
const Scalar & get_hsurf() const
Returns the field h_surf.
Definition: spheroid.h:223
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
const Tensor & delta() const
Computes the delta coefficients for covariant derivative.
Definition: spheroid.C:1206
Scalar * p_get_BC_conf_fact_1
Source of Neumann boundary condition on ,.
Definition: excision_surf.h:84
Tensor down(int ind, const Metric &gam) const
Computes a new tensor by lowering an index of *this.
Coord cost
Definition: map.h:740