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"
13namespace Lorene {
14//Sym_tensor secmembre_kerr ( const Sym_tensor& hij, const Sym_tensor& aa,const Scalar& nn,const Scalar& ppsi,const Vector& bb) {
15const 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
142Scalar 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
148Scalar 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}
154return *p_get_BC_lapse_4 ;
155}
156
157
158}
Spherical orthonormal vectorial bases (triads).
Definition base_vect.h:308
Sym_tensor Kij
The 3-d extrinsic curvature on the slice.
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)
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.
Scalar conf_fact
The value of the conformal factor on the 3-slice.
Metric gamij
The 3-d metric on the slice.
double delta_t
The time step for evolution in parabolic drivers.
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Flat metric for tensor calculation.
Definition metric.h:261
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition metric.C:290
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
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
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...
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition scalar.C:380
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Scalar derive_lie(const Vector &v) const
Computes the derivative of this along a vector field v.
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
Sym_tensor derive_lie(const Vector &v) const
Computes the Lie derivative of this with respect to some vector field v.
Definition sym_tensor.C:360
void ylm()
Computes the coefficients of *this.
Definition valeur_ylm.C:138
Tensor field of valence 1.
Definition vector.h:188
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition vector.C:381
Sym_tensor ope_killing_conf(const Metric &gam) const
Computes the conformal Killing operator associated with a given metric.
Definition vector.C:462
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
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 .
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition tensor.C:816
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 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
Standard units of space, time and mass.