LORENE
excision_lapse_4.C
1 
2 // Header Lorene:
3 #include "nbr_spx.h"
4 #include "utilitaires.h"
5 #include "graphique.h"
6 #include "math.h"
7 #include "metric.h"
8 #include "param.h"
9 #include "param_elliptic.h"
10 #include "tensor.h"
11 #include "unites.h"
12 #include "excision_surf.h"
13 namespace Lorene {
14 //Sym_tensor secmembre_kerr ( const Sym_tensor& hij, const Sym_tensor& aa,const Scalar& nn,const Scalar& ppsi,const Vector& bb) {
15 const Scalar& Excision_surf::get_BC_lapse_4 (Scalar& old_nn, Vector& beta_point, Sym_tensor& strain_tens) const{
16 
17  using namespace Unites;
18  if (p_get_BC_lapse_4 == 0x0){
19 
20  Scalar nn = lapse;
21  Scalar ppsi = conf_fact;
22  Vector bb = shift;
23  Sym_tensor aa = Kij.up_down(gamij)*ppsi*ppsi*ppsi*ppsi;
24 
25  const int nz = (*(aa.get_mp().get_mg())).get_nzone();
26 
27  Sym_tensor hij = gamij.con();
28  for (int ii=1; ii<=3; ii++)
29  for(int jj=1; jj<=3; jj++)
30  { hij.set(ii,jj).annule_hard();
31  }
32  hij.annule_domain(nz-1);
33  hij.std_spectral_base(); // Set non conformally flat part to zero.
34 
35 
36  const Vector& beta = bb;
37 
38  const Scalar& psi4 = ppsi*ppsi*ppsi*ppsi;
39  Scalar ln_psi = log(ppsi); ln_psi.std_spectral_base();
40  ln_psi.annule_domain(nz-1);
41 
42  const Scalar qq = nn*ppsi*ppsi;
43 
44 
45  const Metric_flat& ff = (hij.get_mp()).flat_met_spher() ;
46 
47  const Sym_tensor& tgam_uu = ff.con();
48 
49  const Base_vect_spher& otriad = hij.get_mp().get_bvect_spher();
50 
51  Scalar nn_point(hij.get_mp());
52  nn_point.annule_hard();
53  nn_point.annule_domain(nz -1);
54 
55  nn_point.std_spectral_base();
56 
57 
58 
59  const Sym_tensor& tgam_dd = ff.cov() ; // {\tilde \gamma}_{ij}
60  const Vector& dln_psi = ln_psi.derive_cov(ff) ; // D_i ln(Psi)
61  const Vector& tdln_psi_u = ln_psi.derive_con(ff) ; // tD^i ln(Psi)
62  const Vector& tdnn_u = nn.derive_con(ff) ; // tD^i N
63  const Vector& dqq = qq.derive_cov(ff) ; // D_i Q
64  const Scalar& div_beta = beta.divergence(ff) ; // D_k beta^k
65  Sym_tensor l_beta = beta.ope_killing_conf(ff) ; // Attention aux headers a inclure
66 
67  //==================================
68  // Source for hij
69  //==================================
70 
71 
72  Scalar tmp(hij.get_mp()) ;
73  Sym_tensor sym_tmp(hij.get_mp(), CON, otriad) ;
74 
75  // Full quadratic part of source for h : S^{ij}
76  // --------------------------------------------
77 
78  Sym_tensor ss(hij.get_mp(), CON, otriad) ; // Source secondaire
79 
80  sym_tmp = nn * (8.* tdln_psi_u * tdln_psi_u)
81  + 4.*( tdln_psi_u * tdnn_u + tdnn_u * tdln_psi_u )
82  - 0.3333333333333333 *
83  ( nn * ( 8.* contract(dln_psi, 0, tdln_psi_u, 0) )
84  + 8.* contract(dln_psi, 0, tdnn_u, 0) ) *tgam_uu ;
85 
86  ss = sym_tmp / psi4 ;
87 
88  sym_tmp = contract(tgam_uu, 1,
89  contract(tgam_uu, 1, dqq.derive_cov(ff), 0), 1) ;
90 
91  sym_tmp.inc_dzpuis() ; // dzpuis : 3 --> 4
92 
93  tmp = qq.derive_con(ff).divergence(ff) ;
94  tmp.inc_dzpuis() ; // dzpuis : 3 --> 4
95 
96  sym_tmp -= 0.3333333333333333 * tmp *tgam_uu ;
97 
98  ss -= sym_tmp / (psi4*ppsi*ppsi) ; // Voir dans quel sens sont construits psi et psi4 (eviter les multiplications d'erreurs)
99 
100  for (int i=1; i<=3; i++) {
101  for (int j=1; j<=i; j++) {
102  tmp = 0 ;
103  for (int k=1; k<=3; k++) {
104  for (int l=1; l<=3; l++) {
105  tmp += tgam_dd(k,l) * aa(i,k) * aa(j,l) ;
106  }
107  }
108  sym_tmp.set(i,j) = tmp ;
109  }
110  }
111 
112  tmp = psi4 * strain_tens.trace(ff) ; // S = S_i^i
113 
114  ss += (2.*nn) * ( sym_tmp );
115  Sym_tensor ss2 =2.*nn*( qpig*(psi4*strain_tens - 0.33333333333333 * tmp * tgam_uu));
116 
117  ss2.annule_domain(nz-1); //Dont care, i want only data on the horizon...
118 
119  ss += -ss2;
120 
121  // Source for h^{ij}
122  // -----------------
123 
124 
125  Sym_tensor source_hh = 2.* nn * ss ;
126  source_hh += 2. * nn.derive_lie(beta) * aa ; //HERE
127 
128  // Term (d/dt - Lie_beta) (L beta)^{ij}--> sym_tmp
129  // ---------------------------------------
130 
131  sym_tmp = beta_point.ope_killing_conf(ff);
132  sym_tmp.annule_domain(nz-1);
133  sym_tmp = sym_tmp - l_beta.derive_lie(beta) ;
134 
135  sym_tmp.annule_domain(nz-1);
136 
137 
138  // Final source:
139  // ------------
140  source_hh += 0.6666666666666666* div_beta * l_beta - sym_tmp ;
141 
142 Scalar dNdt(hij.get_mp());
143 
144  dNdt = -source_hh(1,1)/2.*aa(1,1); // Take any component as long as the aa part does not vanish...
145 
146  dNdt.annule_domain(nz-1);
147 
148 Scalar bound_N = old_nn + dNdt*delta_t; bound_N.std_spectral_base();
149  bound_N.set_spectral_va().ylm();
150 
151  p_get_BC_lapse_4 = new Scalar(bound_N);
152 
153 }
154 return *p_get_BC_lapse_4 ;
155 }
156 
157 
158 }
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:675
Scalar * p_get_BC_lapse_4
Source of Dirichlet condtion on , based on einstein equations (conservation of isotropic gauge) ...
Definition: excision_surf.h:93
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
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:293
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:801
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
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
Flat metric for tensor calculation.
Definition: metric.h:261
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
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
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Definition: scalar_deriv.C:402
Tensor field of valence 1.
Definition: vector.h:188
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
Sym_tensor ope_killing_conf(const Metric &gam) const
Computes the conformal Killing operator associated with a given metric.
Definition: vector.C:468
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:387
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
Scalar derive_lie(const Vector &v) const
Computes the derivative of this along a vector field v.
Definition: scalar_deriv.C:419
Spherical orthonormal vectorial bases (triads).
Definition: base_vect.h:308
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:663
Scalar conf_fact
The value of the conformal factor on the 3-slice.
Definition: excision_surf.h:53
const Scalar & get_BC_lapse_4(Scalar &old_nn, Vector &beta_point, Sym_tensor &strain_tens) const
Source for Dirichlet BC on the lapse, based on einstein equations (conservation of isotropic gauge) ...
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:610
Metric gamij
The 3-d metric on the slice.
Definition: excision_surf.h:62
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
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: tensor.C:935
Sym_tensor derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
Definition: sym_tensor.C:363
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
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.