LORENE
strot_dirac_solvenq.C
1 /*
2  * Solution of the two scalar Poisson equations for rotating stars
3  * in Dirac gauge.
4  *
5  * (see file star_rot_Dirac.h for documentation).
6  *
7  */
8 
9 /*
10  * Copyright (c) 2005 Lap-Ming Lin & Jerome Novak
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License version 2
16  * as published by the Free Software Foundation.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 
31 /*
32  * $Id: strot_dirac_solvenq.C,v 1.6 2016/12/05 16:18:15 j_novak Exp $
33  * $Log: strot_dirac_solvenq.C,v $
34  * Revision 1.6 2016/12/05 16:18:15 j_novak
35  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
36  *
37  * Revision 1.5 2014/10/13 08:53:40 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.4 2014/10/06 15:13:18 j_novak
41  * Modified #include directives to use c++ syntax.
42  *
43  * Revision 1.3 2005/02/17 17:31:29 f_limousin
44  * Change the name of some quantities to be consistent with other classes
45  * (for instance nnn is changed to nn, shift to beta, beta to lnq...)
46  *
47  * Revision 1.2 2005/02/02 09:23:20 lm_lin
48  *
49  * Correct a mistake in the calculation of log(N).
50  *
51  *
52  * $Header: /cvsroot/Lorene/C++/Source/Star/strot_dirac_solvenq.C,v 1.6 2016/12/05 16:18:15 j_novak Exp $
53  *
54  */
55 
56 // C headers
57 #include <cmath>
58 #include <cassert>
59 
60 // Lorene headers
61 #include "star_rot_dirac.h"
62 #include "unites.h"
63 
64 namespace Lorene {
65 void Star_rot_Dirac::solve_logn_f(Scalar& ln_f_new) const {
66 
67  using namespace Unites ;
68 
69  //================================================
70  // Source for log(n) containing matter field terms
71  //================================================
72 
73  Scalar source_logn = qpig* psi4* (ener_euler + s_euler) ;
74 
75  ln_f_new = source_logn.poisson() ;
76 
77 }
78 
79 void Star_rot_Dirac::solve_logn_q(Scalar& ln_q_new) const {
80 
81  const Metric_flat& mets = mp.flat_met_spher() ;
82  const Base_vect_spher& bspher = mp.get_bvect_spher() ;
83 
84  const Vector& dln_psi = ln_psi.derive_cov(mets) ; // D_i ln(Psi)
85  const Vector& dln = logn.derive_cov(mets) ; // D_i ln(N)
86 
87 
88 
89  //=============================================
90  // Source for log(n) containing quadratic terms
91  //=============================================
92 
93  Scalar source_logn = psi4* aa_quad
94  - contract(dln.up_down(mets), 0, dln, 0)
95  - 2.* contract(dln_psi, 0, logn.derive_con(tgamma), 0) ;
96 
97 
98  Tensor_sym tmp(mp, 2, COV, bspher, 0, 1) ;
99  tmp = dln.derive_cov(mets) ;
100  tmp.inc_dzpuis() ; //dzpuis 3 -> 4
101 
102  source_logn -= contract( hh, 0, 1, tmp+dln*dln, 0, 1 ) ;
103 
104  ln_q_new = source_logn.poisson() ;
105 
106 }
107 
108 void Star_rot_Dirac::solve_qqq(Scalar& q_new) const {
109 
110  using namespace Unites ;
111 
112  const Metric_flat& mets = mp.flat_met_spher() ;
113 
114  // Source for Q
115  // ------------
116 
117  const Vector& dln_psi = ln_psi.derive_cov(mets) ; // D_i ln(Psi)
118  const Vector& dqq = qqq.derive_cov(mets) ; // D_i Q
119  const Tensor_sym& dhh = hh.derive_cov(mets) ; // D_k h^{ij}
120  const Tensor_sym& dtgam = tgamma.cov().derive_cov(mets) ;
121  // D_k {\tilde \gamma}_{ij}
122  Scalar source_qq = psi4 * qqq * ( qpig* s_euler + 0.75* aa_quad ) ;
123 
124  Scalar tmp = contract( hh, 0, 1, dqq.derive_cov(mets), 0, 1 ) ;
125  tmp.inc_dzpuis() ;
126 
127  source_qq -= tmp ;
128 
129 
130  // tmp = \tilde{R}_{*}/4 + 2 D_i ln(Psi) \tilde{D}^i ln(Psi)
131  tmp = 0.0625 * contract( dhh, 0, 1, dtgam, 0, 1 ).trace(tgamma)
132  - 0.125 * contract( dhh, 0, 1, dtgam, 0, 2 ).trace(tgamma)
133  + 2.* contract( dln_psi, 0, ln_psi.derive_con(tgamma), 0) ;
134 
135  source_qq += psi2 * ( nn * tmp
136  + 2*contract(dln_psi, 0, nn.derive_con(tgamma), 0) ) ;
137 
138 
139  q_new = source_qq.poisson() + 1. ;
140 
141  if (q_new.get_etat() == ETATUN) q_new.std_spectral_base() ;
142 
143 }
144 }
void solve_logn_f(Scalar &ln_f_new) const
Solution of the two scalar Poisson equations for rotating stars in Dirac gauge.
Map & mp
Mapping associated with the star.
Definition: star.h:180
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.
Flat metric for tensor calculation.
Definition: metric.h:261
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition: scalar_pde.C:139
Sym_tensor_trans hh
is defined by .
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
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:560
Scalar s_euler
Trace of the stress scalar in the Eulerian frame.
Definition: star.h:201
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 inc_dzpuis()
dzpuis += 1 ;
Definition: tenseur.C:1133
Scalar psi4
Conformal factor .
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 logn
Logarithm of the lapse N .
Definition: star.h:222
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
Scalar nn
Lapse function N .
Definition: star.h:225
void solve_logn_q(Scalar &ln_q_new) const
Solution of the two scalar Poisson equations for rotating stars in Dirac gauge.
void solve_qqq(Scalar &q_new) const
Solution of the two scalar Poisson equations for rotating stars in Dirac gauge.
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
Scalar ener_euler
Total energy density in the Eulerian frame.
Definition: star.h:198
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