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
31char star_bin_vel_pot_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star/star_bin_vel_pot.C,v 1.7 2014/10/13 08:53:39 j_novak Exp $" ;
32
33/*
34 * $Id: star_bin_vel_pot.C,v 1.7 2014/10/13 08:53:39 j_novak Exp $
35 * $Log: star_bin_vel_pot.C,v $
36 * Revision 1.7 2014/10/13 08:53:39 j_novak
37 * Lorene classes and functions now belong to the namespace Lorene.
38 *
39 * Revision 1.6 2005/09/13 19:38:31 f_limousin
40 * Reintroduction of the resolution of the equations in cartesian coordinates.
41 *
42 * Revision 1.5 2005/02/24 16:07:57 f_limousin
43 * Change the name of some variables (for instance dcov_logn --> dlogn).
44 *
45 * Revision 1.4 2005/02/17 17:33:11 f_limousin
46 * Change the name of some quantities to be consistent with other classes
47 * (for instance nnn is changed to nn, shift to beta, beta to lnq...)
48 *
49 * Revision 1.3 2004/06/22 12:53:09 f_limousin
50 * Change qq, qq_auto and qq_comp to beta, beta_auto and beta_comp.
51 *
52 * Revision 1.2 2004/06/07 16:26:01 f_limousin
53 * Many modif...
54 *
55 * Revision 1.1 2004/04/08 16:34:08 f_limousin
56 * First version
57 *
58 *
59 * $Header: /cvsroot/Lorene/C++/Source/Star/star_bin_vel_pot.C,v 1.7 2014/10/13 08:53:39 j_novak Exp $
60 *
61 */
62
63// Headers Lorene
64#include "star.h"
65#include "eos.h"
66#include "param.h"
67#include "cmp.h"
68#include "tenseur.h"
69#include "graphique.h"
70#include "utilitaires.h"
71
72// Local prototype
73namespace Lorene {
74Cmp raccord_c1(const Cmp& uu, int l1) ;
75
76double Star_bin::velocity_potential(int mermax, double precis, double relax) {
77
78 int nzm1 = mp.get_mg()->get_nzone() - 1 ;
79
80 //----------------------------------
81 // Specific relativistic enthalpy ---> hhh
82 //----------------------------------
83
84 Scalar hhh = exp(ent) ; // = 1 at the Newtonian limit
85 hhh.std_spectral_base() ;
86
87 //----------------------------------------------
88 // Computation of W^i = - A^2 h Gamma_n B^i/N
89 // See Eq (62) from Gourgoulhon et al. (2001)
90 //----------------------------------------------
91
92 Vector www = hhh * gam_euler * bsn * psi4 ;
93
94 www.change_triad( mp.get_bvect_cart() ) ; // components on the mapping
95 // Cartesian basis
96
97 //-------------------------------------------------
98 // Constant value of W^i at the center of the star
99 //-------------------------------------------------
100
101 Vector v_orb(mp, CON, mp.get_bvect_cart()) ;
102
103 v_orb.set_etat_qcq() ;
104 for (int i=1; i<=3; i++) {
105 v_orb.set(i) = www(i).val_grid_point(0, 0, 0, 0) ;
106 }
107
108 v_orb.set_triad( *(www.get_triad()) ) ;
109 v_orb.std_spectral_base() ;
110
111 //-------------------------------------------------
112 // Source and coefficients a,b for poisson_compact (independent from psi0)
113 //-------------------------------------------------
114
115 Scalar dndh_log = eos.der_nbar_ent(ent, nzet) ;
116
117 // In order to avoid any division by zero in the computation of zeta_h
118 // the value of dndh_log is set to 1 in the external domains:
119 for (int l=nzet; l <= nzm1; l++) {
120 dndh_log.set_domain(l) = 1 ;
121 }
122
124 zeta_h.std_spectral_base() ;
125
127 Vector bb = (1 - zeta_h) * ent.derive_con(flat_spher) +
129
130 Scalar entmb = ent - lnq ;
131
133 v_orb.change_triad(mp.get_bvect_cart()) ;
134
136 dcovdcov_psi0.inc_dzpuis() ;
137
138 // See Eq (63) from Gourgoulhon et al. (2001)
139 Scalar source = contract(www - v_orb, 0, ent.derive_cov(flat), 0)
140 + zeta_h * ( contract(v_orb, 0, entmb.derive_cov(flat), 0) +
142 + contract(hij, 0, 1, (psi0.derive_cov(flat)
143 + v_orb.down(0, flat))*entmb.derive_cov(flat),
144 0, 1)
145 - contract(hij, 0, 1, dcovdcov_psi0, 0, 1))
147 + v_orb.down(0, flat)),
148 0, 1) ;
149
150/*
151 des_meridian(zeta_h,0., 4., "zeta_h", 10) ;
152 arrete() ;
153 des_meridian(bb(1),0., 4., "bb(1)", 10) ;
154 arrete() ;
155 des_meridian(bb(2),0., 4., "bb(2)", 10) ;
156 arrete() ;
157 des_meridian(bb(3),0., 4., "bb(3)", 10) ;
158 arrete() ;
159 des_meridian(psi0,0., 4., "psi0", 10) ;
160 arrete() ;
161 des_meridian(source,0., 4., "source", 10) ;
162 arrete() ;
163*/
164
165
167 v_orb.change_triad(mp.get_bvect_cart()) ;
168
169 source.annule(nzet, nzm1) ;
170
171 //---------------------------------------------------
172 // Resolution by means of Map_radial::poisson_compact
173 //---------------------------------------------------
174
175 Param par ;
176 int niter ;
177 par.add_int(mermax) ;
178 par.add_double(precis, 0) ;
179 par.add_double(relax, 1) ;
180 par.add_int_mod(niter) ;
181
182
183 if (psi0.get_etat() == ETATZERO) {
184 psi0 = 0 ;
185 }
186
189 Cmp psi0_cmp (psi0) ;
190
191 Tenseur bb_cmp(mp, 1, CON, mp.get_bvect_spher()) ;
192 bb_cmp.set_etat_qcq() ;
193 Cmp bb_cmp1 (bb(1)) ;
194 Cmp bb_cmp2 (bb(2)) ;
195 Cmp bb_cmp3 (bb(3)) ;
196 bb_cmp.set(0) = bb_cmp1 ;
197 bb_cmp.set(1) = bb_cmp2 ;
198 bb_cmp.set(2) = bb_cmp3 ;
199
200 source_cmp.va.ylm() ;
201
202 cout << "source" << endl << norme(source_cmp)<< endl ;
203 cout << "zeta_h " << endl << norme(zeta_h_cmp) << endl ;
204 cout << "bb(1)" << endl << norme(bb_cmp(0)) << endl ;
205 cout << "bb(2)" << endl << norme(bb_cmp(1)) << endl ;
206 cout << "bb(3)" << endl << norme(bb_cmp(2)) << endl ;
207 cout << "psiO" << endl << norme(psi0_cmp) << endl ;
208
210
211 psi0 = psi0_cmp ;
212
213 cout << "psiO apres" << endl << norme(psi0) << endl ;
214
215 //---------------------------------------------------
216 // Check of the solution
217 //---------------------------------------------------
218
220
222
223 source.set_spectral_va().ylm_i() ;
224
225 double erreur = diffrel(oper, source)(0) ;
226
227 cout << "Check of the resolution of the continuity equation : "
228 << endl ;
229 cout << " norme(source) : " << norme(source)(0)
230 << " diff oper/source : " << erreur << endl ;
231
232
233 //--------------------------------
234 // Computation of grad(psi)
235 //--------------------------------
236
237 v_orb.change_triad(mp.get_bvect_cart()) ;
239
240 for (int i=1; i<=3; i++)
241 d_psi.set(i) = (psi0.derive_cov(flat))(i) + v_orb(i) ;
242
243 v_orb.change_triad(mp.get_bvect_cart()) ;
245
246 // C^1 continuation of d_psi outside the star
247 // (to ensure a smooth enthalpy field accross the stellar surface)
248 // ----------------------------------------------------------------
249
250 d_psi.annule(nzet, nzm1) ;
251
252 for (int i=1; i<=3; i++) {
253 Cmp d_psi_i (d_psi(i)) ;
254 d_psi_i = raccord_c1(d_psi_i, nzet) ;
255 d_psi.set(i) = d_psi_i ;
256 }
257
258 return erreur ;
259
260}
261}
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Definition map.h:791
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...
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
Parameter storage.
Definition param.h:125
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 Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian 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:784
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Vector bsn
3-vector shift, divided by N, of the rotating coordinates, .
Definition star.h:518
Vector d_psi
Gradient of (in the irrotational case) (Spherical components with respect to the mapping of the star...
Definition star.h:501
double velocity_potential(int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
Scalar psi0
Scalar potential of the non-translational part of fluid 4-velocity (in the irrotational case)
Definition star.h:496
Sym_tensor hij
Total deviation of the inverse conformal metric from the inverse flat metric.
Definition star.h:581
Scalar psi4
Conformal factor .
Definition star.h:552
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the star) .
Definition star.h:562
const Eos & eos
Equation of state of the stellar matter.
Definition star.h:185
Scalar gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition star.h:204
Scalar ent
Log-enthalpy.
Definition star.h:190
Map & mp
Mapping associated with the star.
Definition star.h:180
int nzet
Number of domains of *mp occupied by the star.
Definition star.h:183
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
Tensor handling.
Definition tensor.h:288
Tensor field of valence 1.
Definition vector.h:188
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:270
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:504
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition cmp_math.C:481
virtual void annule(int l_min, int l_max)
Sets the Tensor to zero in several domains.
Definition tensor.C:671
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tensor.h:866
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition app_hor.h:64