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 char 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 
77 namespace Lorene {
79 
80  using namespace Unites ;
81 
82  const Metric_flat& mets = mp.flat_met_spher() ;
83  const Base_vect_spher& bspher = mp.get_bvect_spher() ;
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 
107  Sym_tensor ricci_star(mp, CON, bspher) ;
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  }
124  ricci_star -= sym_tmp ;
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  }
146  ricci_star += sym_tmp ;
147 
148  ricci_star = 0.5 * ricci_star ;
149 
150  // Curvature scalar of conformal metric :
151  // -------------------------------------
152 
153  Scalar tricci_scal =
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 
164  sym_tmp = nn * (ricci_star + 8.* tdln_psi_u * tdln_psi_u)
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 
172  sym_tmp = contract(tgam_uu, 1,
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 
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 
218  Sym_tensor lbh = hh.derive_lie(beta) ;
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  // ---------------------------------------
251  Sym_tensor l_beta = beta.ope_killing_conf(mets) ;
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 
268  Sym_tensor_trans source_hht(mp, bspher, mets) ;
269  source_hht = source_hh ;
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 
274  Scalar h_prev = hh.trace(mets) ;
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 }
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:290
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:462
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:324
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:783
Lorene prototypes.
Definition: app_hor.h:64
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:765
Flat metric for tensor calculation.
Definition: metric.h:261
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
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:808
Sym_tensor ope_killing_conf(const Metric &gam) const
Computes the conformal Killing operator associated with a given metric.
Definition: vector.C:462
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:381
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:1002
Transverse symmetric tensors of rank 2.
Definition: sym_tensor.h:608
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:816
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:280
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:654
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:1037
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:360
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
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:223
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.