LORENE
star_bin_vel_pot.C
1 /*
2  * Method of class Star_bin to compute the velocity scalar potential $\psi$
3  * by solving the continuity equation.
4  *
5  * (see file star.h for documentation).
6  *
7  */
8 
9 /*
10  * Copyright (c) 2004 Francois Limousin
11  * 2000-2001 Eric Gourgoulhon (for precedent version)
12  *
13  * This file is part of LORENE.
14  *
15  * LORENE is free software; you can redistribute it and/or modify
16  * it under the terms of the GNU General Public License as published by
17  * the Free Software Foundation; either version 2 of the License, or
18  * (at your option) any later version.
19  *
20  * LORENE is distributed in the hope that it will be useful,
21  * but WITHOUT ANY WARRANTY; without even the implied warranty of
22  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23  * GNU General Public License for more details.
24  *
25  * You should have received a copy of the GNU General Public License
26  * along with LORENE; if not, write to the Free Software
27  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
28  *
29  */
30 
31 
32 
33 /*
34  * $Id: star_bin_vel_pot.C,v 1.8 2016/12/05 16:18:15 j_novak Exp $
35  * $Log: star_bin_vel_pot.C,v $
36  * Revision 1.8 2016/12/05 16:18:15 j_novak
37  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
38  *
39  * Revision 1.7 2014/10/13 08:53:39 j_novak
40  * Lorene classes and functions now belong to the namespace Lorene.
41  *
42  * Revision 1.6 2005/09/13 19:38:31 f_limousin
43  * Reintroduction of the resolution of the equations in cartesian coordinates.
44  *
45  * Revision 1.5 2005/02/24 16:07:57 f_limousin
46  * Change the name of some variables (for instance dcov_logn --> dlogn).
47  *
48  * Revision 1.4 2005/02/17 17:33:11 f_limousin
49  * Change the name of some quantities to be consistent with other classes
50  * (for instance nnn is changed to nn, shift to beta, beta to lnq...)
51  *
52  * Revision 1.3 2004/06/22 12:53:09 f_limousin
53  * Change qq, qq_auto and qq_comp to beta, beta_auto and beta_comp.
54  *
55  * Revision 1.2 2004/06/07 16:26:01 f_limousin
56  * Many modif...
57  *
58  * Revision 1.1 2004/04/08 16:34:08 f_limousin
59  * First version
60  *
61  *
62  * $Header: /cvsroot/Lorene/C++/Source/Star/star_bin_vel_pot.C,v 1.8 2016/12/05 16:18:15 j_novak Exp $
63  *
64  */
65 
66 // Headers Lorene
67 #include "star.h"
68 #include "eos.h"
69 #include "param.h"
70 #include "cmp.h"
71 #include "tenseur.h"
72 #include "graphique.h"
73 #include "utilitaires.h"
74 
75 // Local prototype
76 namespace Lorene {
77 Cmp raccord_c1(const Cmp& uu, int l1) ;
78 
79 double Star_bin::velocity_potential(int mermax, double precis, double relax) {
80 
81  int nzm1 = mp.get_mg()->get_nzone() - 1 ;
82 
83  //----------------------------------
84  // Specific relativistic enthalpy ---> hhh
85  //----------------------------------
86 
87  Scalar hhh = exp(ent) ; // = 1 at the Newtonian limit
88  hhh.std_spectral_base() ;
89 
90  //----------------------------------------------
91  // Computation of W^i = - A^2 h Gamma_n B^i/N
92  // See Eq (62) from Gourgoulhon et al. (2001)
93  //----------------------------------------------
94 
95  Vector www = hhh * gam_euler * bsn * psi4 ;
96 
97  www.change_triad( mp.get_bvect_cart() ) ; // components on the mapping
98  // Cartesian basis
99 
100  //-------------------------------------------------
101  // Constant value of W^i at the center of the star
102  //-------------------------------------------------
103 
104  Vector v_orb(mp, CON, mp.get_bvect_cart()) ;
105 
106  v_orb.set_etat_qcq() ;
107  for (int i=1; i<=3; i++) {
108  v_orb.set(i) = www(i).val_grid_point(0, 0, 0, 0) ;
109  }
110 
111  v_orb.set_triad( *(www.get_triad()) ) ;
112  v_orb.std_spectral_base() ;
113 
114  //-------------------------------------------------
115  // Source and coefficients a,b for poisson_compact (independent from psi0)
116  //-------------------------------------------------
117 
118  Scalar dndh_log = eos.der_nbar_ent(ent, nzet) ;
119 
120  // In order to avoid any division by zero in the computation of zeta_h
121  // the value of dndh_log is set to 1 in the external domains:
122  for (int l=nzet; l <= nzm1; l++) {
123  dndh_log.set_domain(l) = 1 ;
124  }
125 
126  Scalar zeta_h( ent / dndh_log ) ;
127  zeta_h.std_spectral_base() ;
128 
129  Metric_flat flat_spher (mp.flat_met_spher()) ;
130  Vector bb = (1 - zeta_h) * ent.derive_con(flat_spher) +
131  zeta_h * lnq.derive_con(flat_spher) ;
132 
133  Scalar entmb = ent - lnq ;
134 
136  v_orb.change_triad(mp.get_bvect_cart()) ;
137 
138  Tensor dcovdcov_psi0 = psi0.derive_cov(flat).derive_cov(flat) ;
139  dcovdcov_psi0.inc_dzpuis() ;
140 
141  // See Eq (63) from Gourgoulhon et al. (2001)
142  Scalar source = contract(www - v_orb, 0, ent.derive_cov(flat), 0)
143  + zeta_h * ( contract(v_orb, 0, entmb.derive_cov(flat), 0) +
145  + contract(hij, 0, 1, (psi0.derive_cov(flat)
146  + v_orb.down(0, flat))*entmb.derive_cov(flat),
147  0, 1)
148  - contract(hij, 0, 1, dcovdcov_psi0, 0, 1))
150  + v_orb.down(0, flat)),
151  0, 1) ;
152 
153 /*
154  des_meridian(zeta_h,0., 4., "zeta_h", 10) ;
155  arrete() ;
156  des_meridian(bb(1),0., 4., "bb(1)", 10) ;
157  arrete() ;
158  des_meridian(bb(2),0., 4., "bb(2)", 10) ;
159  arrete() ;
160  des_meridian(bb(3),0., 4., "bb(3)", 10) ;
161  arrete() ;
162  des_meridian(psi0,0., 4., "psi0", 10) ;
163  arrete() ;
164  des_meridian(source,0., 4., "source", 10) ;
165  arrete() ;
166 */
167 
168 
170  v_orb.change_triad(mp.get_bvect_cart()) ;
171 
172  source.annule(nzet, nzm1) ;
173 
174  //---------------------------------------------------
175  // Resolution by means of Map_radial::poisson_compact
176  //---------------------------------------------------
177 
178  Param par ;
179  int niter ;
180  par.add_int(mermax) ;
181  par.add_double(precis, 0) ;
182  par.add_double(relax, 1) ;
183  par.add_int_mod(niter) ;
184 
185 
186  if (psi0.get_etat() == ETATZERO) {
187  psi0 = 0 ;
188  }
189 
190  Cmp source_cmp (source) ;
191  Cmp zeta_h_cmp (zeta_h) ;
192  Cmp psi0_cmp (psi0) ;
193 
194  Tenseur bb_cmp(mp, 1, CON, mp.get_bvect_spher()) ;
195  bb_cmp.set_etat_qcq() ;
196  Cmp bb_cmp1 (bb(1)) ;
197  Cmp bb_cmp2 (bb(2)) ;
198  Cmp bb_cmp3 (bb(3)) ;
199  bb_cmp.set(0) = bb_cmp1 ;
200  bb_cmp.set(1) = bb_cmp2 ;
201  bb_cmp.set(2) = bb_cmp3 ;
202 
203  source_cmp.va.ylm() ;
204 
205  cout << "source" << endl << norme(source_cmp)<< endl ;
206  cout << "zeta_h " << endl << norme(zeta_h_cmp) << endl ;
207  cout << "bb(1)" << endl << norme(bb_cmp(0)) << endl ;
208  cout << "bb(2)" << endl << norme(bb_cmp(1)) << endl ;
209  cout << "bb(3)" << endl << norme(bb_cmp(2)) << endl ;
210  cout << "psiO" << endl << norme(psi0_cmp) << endl ;
211 
212  mp.poisson_compact(source_cmp, zeta_h_cmp, bb_cmp, par, psi0_cmp ) ;
213 
214  psi0 = psi0_cmp ;
215 
216  cout << "psiO apres" << endl << norme(psi0) << endl ;
217 
218  //---------------------------------------------------
219  // Check of the solution
220  //---------------------------------------------------
221 
222  Scalar bb_dpsi0 = contract(bb, 0, psi0.derive_cov(flat_spher), 0) ;
223 
224  Scalar oper = zeta_h * psi0.laplacian() + bb_dpsi0 ;
225 
226  source.set_spectral_va().ylm_i() ;
227 
228  double erreur = diffrel(oper, source)(0) ;
229 
230  cout << "Check of the resolution of the continuity equation : "
231  << endl ;
232  cout << " norme(source) : " << norme(source)(0)
233  << " diff oper/source : " << erreur << endl ;
234 
235 
236  //--------------------------------
237  // Computation of grad(psi)
238  //--------------------------------
239 
240  v_orb.change_triad(mp.get_bvect_cart()) ;
242 
243  for (int i=1; i<=3; i++)
244  d_psi.set(i) = (psi0.derive_cov(flat))(i) + v_orb(i) ;
245 
246  v_orb.change_triad(mp.get_bvect_cart()) ;
248 
249  // C^1 continuation of d_psi outside the star
250  // (to ensure a smooth enthalpy field accross the stellar surface)
251  // ----------------------------------------------------------------
252 
253  d_psi.annule(nzet, nzm1) ;
254 
255  for (int i=1; i<=3; i++) {
256  Cmp d_psi_i (d_psi(i)) ;
257  d_psi_i = raccord_c1(d_psi_i, nzet) ;
258  d_psi.set(i) = d_psi_i ;
259  }
260 
261  return erreur ;
262 
263 }
264 }
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition: tensor.C:490
Scalar psi0
Scalar potential of the non-translational part of fluid 4-velocity (in the irrotational case) ...
Definition: star.h:496
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:273
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition: param.C:249
Map & mp
Mapping associated with the star.
Definition: star.h:180
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:134
Cmp der_nbar_ent(const Cmp &ent, int nzet, int l_min=0, Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy and extra parameters.
Definition: eos.C:431
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition: scalar.C:397
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:801
Vector d_psi
Gradient of (in the irrotational case) (Spherical components with respect to the mapping of the star...
Definition: star.h:501
Lorene prototypes.
Definition: app_hor.h:67
const Eos & eos
Equation of state of the stellar matter.
Definition: star.h:185
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:141
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
Definition: scalar_deriv.C:436
Flat metric for tensor calculation.
Definition: metric.h:261
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
virtual void poisson_compact(const Cmp &source, const Cmp &aa, const Tenseur &bb, const Param &par, Cmp &psi) const =0
Resolution of the elliptic equation in the case where the stellar interior is covered by a single do...
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
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:560
Scalar ent
Log-enthalpy.
Definition: star.h:190
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
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Definition: scalar.h:621
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:507
int nzet
Number of domains of *mp occupied by the star.
Definition: star.h:183
Scalar gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition: star.h:204
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:484
Sym_tensor hij
Total deviation of the inverse conformal metric from the inverse flat metric.
Definition: star.h:581
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tensor.h:879
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the star) ...
Definition: star.h:562
Parameter storage.
Definition: param.h:125
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Tensor handling.
Definition: tensor.h:294
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
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
double velocity_potential(int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:724
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:809
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:318
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition: tensor.C:528
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:610
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:302
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:652
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
virtual void annule(int l_min, int l_max)
Sets the Tensor to zero in several domains.
Definition: tensor.C:680
Vector bsn
3-vector shift, divided by N, of the rotating coordinates, .
Definition: star.h:518
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:464
Scalar psi4
Conformal factor .
Definition: star.h:552
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
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:304
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition: param.C:388
Tensor down(int ind, const Metric &gam) const
Computes a new tensor by lowering an index of *this.