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 
28 
29 
30 /*
31  * $Id: strot_dirac_solvehij.C,v 1.11 2016/12/05 16:18:15 j_novak Exp $
32  * $Log: strot_dirac_solvehij.C,v $
33  * Revision 1.11 2016/12/05 16:18:15 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.10 2014/10/13 08:53:40 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.9 2010/10/11 10:21:31 j_novak
40  * Less output
41  *
42  * Revision 1.8 2005/11/28 14:45:16 j_novak
43  * Improved solution of the Poisson tensor equation in the case of a transverse
44  * tensor.
45  *
46  * Revision 1.7 2005/09/16 14:04:49 j_novak
47  * The equation for hij is now solved only for mer > mer_fix_omega. It uses the
48  * Poisson solver of the class Sym_tensor_trans.
49  *
50  * Revision 1.6 2005/04/20 14:26:29 j_novak
51  * Removed some outputs.
52  *
53  * Revision 1.5 2005/04/05 09:24:05 j_novak
54  * minor modifs
55  *
56  * Revision 1.4 2005/02/17 17:32:16 f_limousin
57  * Change the name of some quantities to be consistent with other classes
58  * (for instance nnn is changed to nn, shift to beta, beta to lnq...)
59  *
60  * Revision 1.3 2005/02/16 12:51:49 j_novak
61  * Corrected an error on the matter source terms.
62  *
63  * Revision 1.2 2005/02/09 13:34:56 lm_lin
64  *
65  * Remove the Laplacian of hij from the term source_hh and fix an overall
66  * minus error.
67  *
68  * Revision 1.1 2005/01/31 08:51:48 j_novak
69  * New files for rotating stars in Dirac gauge (still under developement).
70  *
71  *
72  * $Header: /cvsroot/Lorene/C++/Source/Star/strot_dirac_solvehij.C,v 1.11 2016/12/05 16:18:15 j_novak Exp $
73  *
74  */
75 
76 // Lorene headers
77 #include "star_rot_dirac.h"
78 #include "unites.h"
79 
80 namespace Lorene {
82 
83  using namespace Unites ;
84 
85  const Metric_flat& mets = mp.flat_met_spher() ;
86  const Base_vect_spher& bspher = mp.get_bvect_spher() ;
87 
88  //==================================
89  // Source for hij
90  //==================================
91 
92  const Vector& dln_psi = ln_psi.derive_cov(mets) ; // D_i ln(Psi)
93  const Vector& dqq = qqq.derive_cov(mets) ; // D_i Q
94  const Tensor_sym& dhh = hh.derive_cov(mets) ; // D_k h^{ij}
95  const Tensor_sym& dtgam = tgamma.cov().derive_cov(mets) ;
96  const Sym_tensor& tgam_dd = tgamma.cov() ; // {\tilde \gamma}_{ij}
97  const Sym_tensor& tgam_uu = tgamma.con() ; // {\tilde \gamma}^{ij}
98  const Vector& tdln_psi_u = ln_psi.derive_con(tgamma) ; // tD^i ln(Psi)
99  const Vector& tdnn_u = nn.derive_con(tgamma) ; // tD^i N
100 // const Scalar& div_beta = beta.divergence(mets) ; // D_k beta^k
101  Scalar div_beta(mp) ; //##
102  div_beta.set_etat_zero() ;
103 
104  Scalar tmp(mp) ;
105  Sym_tensor sym_tmp(mp, CON, bspher) ;
106 
107  // Quadratic part of the Ricci tensor of gam_tilde
108  // ------------------------------------------------
109 
110  Sym_tensor ricci_star(mp, CON, bspher) ;
111 
112  ricci_star = contract(hh, 0, 1, dhh.derive_cov(mets), 2, 3) ;
113 
114  ricci_star.inc_dzpuis() ; // dzpuis : 3 --> 4
115 
116  for (int i=1; i<=3; i++) {
117  for (int j=1; j<=i; j++) {
118  tmp = 0 ;
119  for (int k=1; k<=3; k++) {
120  for (int l=1; l<=3; l++) {
121  tmp += dhh(i,k,l) * dhh(j,l,k) ;
122  }
123  }
124  sym_tmp.set(i,j) = tmp ;
125  }
126  }
127  ricci_star -= sym_tmp ;
128 
129  for (int i=1; i<=3; i++) {
130  for (int j=1; j<=i; j++) {
131  tmp = 0 ;
132  for (int k=1; k<=3; k++) {
133  for (int l=1; l<=3; l++) {
134  for (int m=1; m<=3; m++) {
135  for (int n=1; n<=3; n++) {
136 
137  tmp += 0.5 * tgam_uu(i,k)* tgam_uu(j,l)
138  * dhh(m,n,k) * dtgam(m,n,l)
139  + tgam_dd(n,l) * dhh(m,n,k)
140  * (tgam_uu(i,k) * dhh(j,l,m) + tgam_uu(j,k) * dhh(i,l,m) )
141  - tgam_dd(k,l) *tgam_uu(m,n) * dhh(i,k,m) * dhh(j,l,n) ;
142  }
143  }
144  }
145  }
146  sym_tmp.set(i,j) = tmp ;
147  }
148  }
149  ricci_star += sym_tmp ;
150 
151  ricci_star = 0.5 * ricci_star ;
152 
153  // Curvature scalar of conformal metric :
154  // -------------------------------------
155 
156  Scalar tricci_scal =
157  0.25 * contract(tgam_uu, 0, 1,
158  contract(dhh, 0, 1, dtgam, 0, 1), 0, 1 )
159  - 0.5 * contract(tgam_uu, 0, 1,
160  contract(dhh, 0, 1, dtgam, 0, 2), 0, 1 ) ;
161 
162  // Full quadratic part of source for h : S^{ij}
163  // --------------------------------------------
164 
165  Sym_tensor ss(mp, CON, bspher) ;
166 
167  sym_tmp = nn * (ricci_star + 8.* tdln_psi_u * tdln_psi_u)
168  + 4.*( tdln_psi_u * tdnn_u + tdnn_u * tdln_psi_u )
169  - 0.3333333333333333 *
170  ( nn * (tricci_scal + 8.* contract(dln_psi, 0, tdln_psi_u, 0) )
171  + 8.* contract(dln_psi, 0, tdnn_u, 0) ) *tgam_uu ;
172 
173  ss = sym_tmp / psi4 ;
174 
175  sym_tmp = contract(tgam_uu, 1,
176  contract(tgam_uu, 1, dqq.derive_cov(mets), 0), 1) ;
177 
178  sym_tmp.inc_dzpuis() ; // dzpuis : 3 --> 4
179 
180  for (int i=1; i<=3; i++) {
181  for (int j=1; j<=i; j++) {
182  tmp = 0 ;
183  for (int k=1; k<=3; k++) {
184  for (int l=1; l<=3; l++) {
185  tmp += ( hh(i,k)*dhh(l,j,k) + hh(k,j)*dhh(i,l,k)
186  - hh(k,l)*dhh(i,j,k) ) * dqq(l) ;
187  }
188  }
189  sym_tmp.set(i,j) += 0.5 * tmp ;
190  }
191  }
192 
194  tmp.inc_dzpuis() ; // dzpuis : 3 --> 4
195 
196  sym_tmp -= 0.3333333333333333 * tmp *tgam_uu ;
197 
198  ss -= sym_tmp / (psi4*psi2) ;
199 
200  for (int i=1; i<=3; i++) {
201  for (int j=1; j<=i; j++) {
202  tmp = 0 ;
203  for (int k=1; k<=3; k++) {
204  for (int l=1; l<=3; l++) {
205  tmp += tgam_dd(k,l) * aa(i,k) * aa(j,l) ;
206  }
207  }
208  sym_tmp.set(i,j) = tmp ;
209  }
210  }
211 
212  ss += (2.*nn) * ( sym_tmp - qpig*( psi4* stress_euler
213  - 0.3333333333333333 * s_euler * tgam_uu )
214  ) ;
215 
216 // maxabs(ss, "ss tot") ;
217 
218  // Source for h^{ij}
219  // -----------------
220 
221  Sym_tensor lbh = hh.derive_lie(beta) ;
222 
223  Sym_tensor source_hh = - lbh.derive_lie(beta) ;
224  source_hh.inc_dzpuis() ;
225 
226  source_hh += 2.* nn * ss ;
227 
228  source_hh += - 1.3333333333333333 * div_beta* lbh
229  - 2. * nn.derive_lie(beta) * aa ;
230 
231 
232  for (int i=1; i<=3; i++) {
233  for (int j=1; j<=i; j++) {
234  tmp = 0 ;
235  for (int k=1; k<=3; k++) {
236  tmp += ( hh.derive_con(mets)(k,j,i)
237  + hh.derive_con(mets)(i,k,j)
238  - hh.derive_con(mets)(i,j,k) ) * dqq(k) ;
239  }
240  sym_tmp.set(i,j) = tmp ;
241  }
242  }
243 
244  source_hh -= nn / (psi4*psi2) * sym_tmp ;
245 
246  tmp = - div_beta.derive_lie(beta) ;
247  tmp.inc_dzpuis() ;
248  source_hh += 0.6666666666666666*
249  ( tmp - 0.6666666666666666* div_beta * div_beta ) * hh ;
250 
251 
252  // Term (d/dt - Lie_beta) (L beta)^{ij}--> sym_tmp
253  // ---------------------------------------
254  Sym_tensor l_beta = beta.ope_killing_conf(mets) ;
255 
256  sym_tmp = - l_beta.derive_lie(beta) ;
257 
258  sym_tmp.inc_dzpuis() ;
259 
260  // Final source:
261  // ------------
262  source_hh += 0.6666666666666666* div_beta * l_beta - sym_tmp ;
263 
264  source_hh = - ( psi4/nn/nn )*source_hh ;
265 
266  for (int i=1; i<=3; i++)
267  for (int j=i; j<=3; j++) {
268  source_hh.set(i,j).set_dzpuis(4) ;
269  }
270 
271  Sym_tensor_trans source_hht(mp, bspher, mets) ;
272  source_hht = source_hh ;
273 // cout << " Max( divergence of source_hh ) " << endl ;
274 // for (int i=1; i<=3; i++)
275 // cout << max(abs(source_htt.divergence(mets)(i))) ;
276 
277  Scalar h_prev = hh.trace(mets) ;
278  hij_new = source_hht.poisson(&h_prev) ;
279 
280  if (mp.get_mg()->get_np(0) == 1) { //Axial symmetry
281  hij_new.set(1,3).set_etat_zero() ;
282  hij_new.set(2,3).set_etat_zero() ;
283  }
284 
285 }
286 }
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:293
Map & mp
Mapping associated with the star.
Definition: star.h:180
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:330
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.
Sym_tensor stress_euler
Spatial part of the stress-energy tensor with respect to the Eulerian observer.
Definition: star.h:212
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
Sym_tensor_trans hh
is defined by .
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...
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Definition: star.h:201
Vector beta
Shift vector.
Definition: star.h:228
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
const Tensor_sym & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:814
Sym_tensor ope_killing_conf(const Metric &gam) const
Computes the conformal Killing operator associated with a given metric.
Definition: vector.C:468
Scalar psi4
Conformal factor .
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
Transverse symmetric tensors of rank 2.
Definition: sym_tensor.h:611
Scalar derive_lie(const Vector &v) const
Computes the derivative of this along a vector field v.
Definition: scalar_deriv.C:419
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:825
Spherical orthonormal vectorial bases (triads).
Definition: base_vect.h:308
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:283
void solve_hij(Sym_tensor_trans &hij_new) const
Solution of the tensor Poisson equation for rotating stars in Dirac gauge.
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:663
Scalar nn
Lapse function N .
Definition: star.h:225
Sym_tensor_trans poisson(const Scalar *h_guess=0x0) const
Computes the solution of a tensorial transverse Poisson equation with *this as a source: In particu...
Symmetric tensors (with respect to two of their arguments).
Definition: tensor.h:1050
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
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
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
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.