LORENE
strot_dirac_solvehij.C
1/*
2 * Solution of the tensor Poisson equation for rotating stars in Dirac gauge.
3 *
4 * (see file star_rot_dirac.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2005 Lap-Ming Lin & Jerome Novak
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License version 2
15 * as published by the Free Software Foundation.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28char strot_dirac_solvehij_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star/strot_dirac_solvehij.C,v 1.10 2014/10/13 08:53:40 j_novak Exp $" ;
29
30/*
31 * $Id: strot_dirac_solvehij.C,v 1.10 2014/10/13 08:53:40 j_novak Exp $
32 * $Log: strot_dirac_solvehij.C,v $
33 * Revision 1.10 2014/10/13 08:53:40 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.9 2010/10/11 10:21:31 j_novak
37 * Less output
38 *
39 * Revision 1.8 2005/11/28 14:45:16 j_novak
40 * Improved solution of the Poisson tensor equation in the case of a transverse
41 * tensor.
42 *
43 * Revision 1.7 2005/09/16 14:04:49 j_novak
44 * The equation for hij is now solved only for mer > mer_fix_omega. It uses the
45 * Poisson solver of the class Sym_tensor_trans.
46 *
47 * Revision 1.6 2005/04/20 14:26:29 j_novak
48 * Removed some outputs.
49 *
50 * Revision 1.5 2005/04/05 09:24:05 j_novak
51 * minor modifs
52 *
53 * Revision 1.4 2005/02/17 17:32:16 f_limousin
54 * Change the name of some quantities to be consistent with other classes
55 * (for instance nnn is changed to nn, shift to beta, beta to lnq...)
56 *
57 * Revision 1.3 2005/02/16 12:51:49 j_novak
58 * Corrected an error on the matter source terms.
59 *
60 * Revision 1.2 2005/02/09 13:34:56 lm_lin
61 *
62 * Remove the Laplacian of hij from the term source_hh and fix an overall
63 * minus error.
64 *
65 * Revision 1.1 2005/01/31 08:51:48 j_novak
66 * New files for rotating stars in Dirac gauge (still under developement).
67 *
68 *
69 * $Header: /cvsroot/Lorene/C++/Source/Star/strot_dirac_solvehij.C,v 1.10 2014/10/13 08:53:40 j_novak Exp $
70 *
71 */
72
73// Lorene headers
74#include "star_rot_dirac.h"
75#include "unites.h"
76
77namespace Lorene {
79
80 using namespace Unites ;
81
84
85 //==================================
86 // Source for hij
87 //==================================
88
89 const Vector& dln_psi = ln_psi.derive_cov(mets) ; // D_i ln(Psi)
90 const Vector& dqq = qqq.derive_cov(mets) ; // D_i Q
91 const Tensor_sym& dhh = hh.derive_cov(mets) ; // D_k h^{ij}
92 const Tensor_sym& dtgam = tgamma.cov().derive_cov(mets) ;
93 const Sym_tensor& tgam_dd = tgamma.cov() ; // {\tilde \gamma}_{ij}
94 const Sym_tensor& tgam_uu = tgamma.con() ; // {\tilde \gamma}^{ij}
95 const Vector& tdln_psi_u = ln_psi.derive_con(tgamma) ; // tD^i ln(Psi)
96 const Vector& tdnn_u = nn.derive_con(tgamma) ; // tD^i N
97// const Scalar& div_beta = beta.divergence(mets) ; // D_k beta^k
98 Scalar div_beta(mp) ; //##
99 div_beta.set_etat_zero() ;
100
101 Scalar tmp(mp) ;
102 Sym_tensor sym_tmp(mp, CON, bspher) ;
103
104 // Quadratic part of the Ricci tensor of gam_tilde
105 // ------------------------------------------------
106
108
109 ricci_star = contract(hh, 0, 1, dhh.derive_cov(mets), 2, 3) ;
110
111 ricci_star.inc_dzpuis() ; // dzpuis : 3 --> 4
112
113 for (int i=1; i<=3; i++) {
114 for (int j=1; j<=i; j++) {
115 tmp = 0 ;
116 for (int k=1; k<=3; k++) {
117 for (int l=1; l<=3; l++) {
118 tmp += dhh(i,k,l) * dhh(j,l,k) ;
119 }
120 }
121 sym_tmp.set(i,j) = tmp ;
122 }
123 }
125
126 for (int i=1; i<=3; i++) {
127 for (int j=1; j<=i; j++) {
128 tmp = 0 ;
129 for (int k=1; k<=3; k++) {
130 for (int l=1; l<=3; l++) {
131 for (int m=1; m<=3; m++) {
132 for (int n=1; n<=3; n++) {
133
134 tmp += 0.5 * tgam_uu(i,k)* tgam_uu(j,l)
135 * dhh(m,n,k) * dtgam(m,n,l)
136 + tgam_dd(n,l) * dhh(m,n,k)
137 * (tgam_uu(i,k) * dhh(j,l,m) + tgam_uu(j,k) * dhh(i,l,m) )
138 - tgam_dd(k,l) *tgam_uu(m,n) * dhh(i,k,m) * dhh(j,l,n) ;
139 }
140 }
141 }
142 }
143 sym_tmp.set(i,j) = tmp ;
144 }
145 }
147
148 ricci_star = 0.5 * ricci_star ;
149
150 // Curvature scalar of conformal metric :
151 // -------------------------------------
152
154 0.25 * contract(tgam_uu, 0, 1,
155 contract(dhh, 0, 1, dtgam, 0, 1), 0, 1 )
156 - 0.5 * contract(tgam_uu, 0, 1,
157 contract(dhh, 0, 1, dtgam, 0, 2), 0, 1 ) ;
158
159 // Full quadratic part of source for h : S^{ij}
160 // --------------------------------------------
161
162 Sym_tensor ss(mp, CON, bspher) ;
163
165 + 4.*( tdln_psi_u * tdnn_u + tdnn_u * tdln_psi_u )
166 - 0.3333333333333333 *
167 ( nn * (tricci_scal + 8.* contract(dln_psi, 0, tdln_psi_u, 0) )
168 + 8.* contract(dln_psi, 0, tdnn_u, 0) ) *tgam_uu ;
169
170 ss = sym_tmp / psi4 ;
171
173 contract(tgam_uu, 1, dqq.derive_cov(mets), 0), 1) ;
174
175 sym_tmp.inc_dzpuis() ; // dzpuis : 3 --> 4
176
177 for (int i=1; i<=3; i++) {
178 for (int j=1; j<=i; j++) {
179 tmp = 0 ;
180 for (int k=1; k<=3; k++) {
181 for (int l=1; l<=3; l++) {
182 tmp += ( hh(i,k)*dhh(l,j,k) + hh(k,j)*dhh(i,l,k)
183 - hh(k,l)*dhh(i,j,k) ) * dqq(l) ;
184 }
185 }
186 sym_tmp.set(i,j) += 0.5 * tmp ;
187 }
188 }
189
190 tmp = qqq.derive_con(tgamma).divergence(tgamma) ;
191 tmp.inc_dzpuis() ; // dzpuis : 3 --> 4
192
193 sym_tmp -= 0.3333333333333333 * tmp *tgam_uu ;
194
195 ss -= sym_tmp / (psi4*psi2) ;
196
197 for (int i=1; i<=3; i++) {
198 for (int j=1; j<=i; j++) {
199 tmp = 0 ;
200 for (int k=1; k<=3; k++) {
201 for (int l=1; l<=3; l++) {
202 tmp += tgam_dd(k,l) * aa(i,k) * aa(j,l) ;
203 }
204 }
205 sym_tmp.set(i,j) = tmp ;
206 }
207 }
208
209 ss += (2.*nn) * ( sym_tmp - qpig*( psi4* stress_euler
210 - 0.3333333333333333 * s_euler * tgam_uu )
211 ) ;
212
213// maxabs(ss, "ss tot") ;
214
215 // Source for h^{ij}
216 // -----------------
217
219
220 Sym_tensor source_hh = - lbh.derive_lie(beta) ;
221 source_hh.inc_dzpuis() ;
222
223 source_hh += 2.* nn * ss ;
224
225 source_hh += - 1.3333333333333333 * div_beta* lbh
226 - 2. * nn.derive_lie(beta) * aa ;
227
228
229 for (int i=1; i<=3; i++) {
230 for (int j=1; j<=i; j++) {
231 tmp = 0 ;
232 for (int k=1; k<=3; k++) {
233 tmp += ( hh.derive_con(mets)(k,j,i)
234 + hh.derive_con(mets)(i,k,j)
235 - hh.derive_con(mets)(i,j,k) ) * dqq(k) ;
236 }
237 sym_tmp.set(i,j) = tmp ;
238 }
239 }
240
241 source_hh -= nn / (psi4*psi2) * sym_tmp ;
242
243 tmp = - div_beta.derive_lie(beta) ;
244 tmp.inc_dzpuis() ;
245 source_hh += 0.6666666666666666*
246 ( tmp - 0.6666666666666666* div_beta * div_beta ) * hh ;
247
248
249 // Term (d/dt - Lie_beta) (L beta)^{ij}--> sym_tmp
250 // ---------------------------------------
252
253 sym_tmp = - l_beta.derive_lie(beta) ;
254
255 sym_tmp.inc_dzpuis() ;
256
257 // Final source:
258 // ------------
259 source_hh += 0.6666666666666666* div_beta * l_beta - sym_tmp ;
260
261 source_hh = - ( psi4/nn/nn )*source_hh ;
262
263 for (int i=1; i<=3; i++)
264 for (int j=i; j<=3; j++) {
265 source_hh.set(i,j).set_dzpuis(4) ;
266 }
267
270// cout << " Max( divergence of source_hh ) " << endl ;
271// for (int i=1; i<=3; i++)
272// cout << max(abs(source_htt.divergence(mets)(i))) ;
273
275 hij_new = source_hht.poisson(&h_prev) ;
276
277 if (mp.get_mg()->get_np(0) == 1) { //Axial symmetry
278 hij_new.set(1,3).set_etat_zero() ;
279 hij_new.set(2,3).set_etat_zero() ;
280 }
281
282}
283}
Spherical orthonormal vectorial bases (triads).
Definition base_vect.h:308
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
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
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.
Definition metric.C:290
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition metric.C:280
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
const 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.
Sym_tensor_trans hh
is defined by .
Scalar psi4
Conformal factor .
void solve_hij(Sym_tensor_trans &hij_new) const
Solution of the tensor Poisson equation for rotating stars in Dirac gauge.
Scalar nn
Lapse function N .
Definition star.h:225
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Definition star.h:201
Sym_tensor stress_euler
Spatial part of the stress-energy tensor with respect to the Eulerian observer.
Definition star.h:212
Map & mp
Mapping associated with the star.
Definition star.h:180
Vector beta
Shift vector.
Definition star.h:228
Transverse symmetric tensors of rank 2.
Definition sym_tensor.h:608
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
Symmetric tensors (with respect to two of their arguments).
Definition tensor.h:1037
Tensor field of valence 1.
Definition vector.h:188
Sym_tensor ope_killing_conf(const Metric &gam) const
Computes the conformal Killing operator associated with a given metric.
Definition vector.C:462
const Tensor_sym & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
const Tensor_sym & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
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.