LORENE
star_bin_vel_pot_xcts.C
1/*
2 * Method of class Star_bin_xcts to compute the velocity scalar
3 * potential $\psi$ by solving the continuity equation
4 * (see file star.h for documentation).
5 */
6
7/*
8 * Copyright (c) 2010 Michal Bejger
9 *
10 * This file is part of LORENE.
11 *
12 * LORENE is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU General Public License version 2
14 * as published by the Free Software Foundation.
15 *
16 * LORENE is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with LORENE; if not, write to the Free Software
23 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24 *
25 */
26
27char star_bin_vel_pot_xcts_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star/star_bin_vel_pot_xcts.C,v 1.7 2014/10/13 08:53:39 j_novak Exp $" ;
28
29/*
30 * $Id: star_bin_vel_pot_xcts.C,v 1.7 2014/10/13 08:53:39 j_novak Exp $
31 * $Log: star_bin_vel_pot_xcts.C,v $
32 * Revision 1.7 2014/10/13 08:53:39 j_novak
33 * Lorene classes and functions now belong to the namespace Lorene.
34 *
35 * Revision 1.6 2010/12/09 10:47:31 m_bejger
36 * Small changes, definition of lnPsi2N
37 *
38 * Revision 1.5 2010/10/26 20:01:06 m_bejger
39 * Cleanup
40 *
41 * Revision 1.4 2010/10/18 20:56:15 m_bejger
42 * Generalization for many domains in the star: buggy
43 *
44 * Revision 1.3 2010/06/17 15:07:10 m_bejger
45 * Psi4, lnPsi2N corrected
46 *
47 * Revision 1.2 2010/06/15 08:05:55 m_bejger
48 * Various fields were lacking bases
49 *
50 * Revision 1.1 2010/05/04 07:51:05 m_bejger
51 * Initial version
52 *
53 * $Header: /cvsroot/Lorene/C++/Source/Star/star_bin_vel_pot_xcts.C,v 1.7 2014/10/13 08:53:39 j_novak Exp $
54 *
55 */
56
57// Headers Lorene
58#include "star.h"
59#include "eos.h"
60#include "param.h"
61#include "cmp.h"
62#include "tenseur.h"
63#include "graphique.h"
64#include "utilitaires.h"
65
66// Local prototype
67namespace Lorene {
68Cmp raccord_c1(const Cmp& uu, int l1) ;
69
71 double precis,
72 double relax) {
73
74 int nzm1 = mp.get_mg()->get_nzone() - 1 ;
75
76 //----------------------------------
77 // Specific relativistic enthalpy ---> hhh
78 //----------------------------------
79
80 Scalar hhh = exp(ent) ; // = 1 at the Newtonian limit
81 hhh.std_spectral_base() ;
82
83 //------------------------------------------------------------------
84 // Computation of W^i = Psi^4 h Gamma_n B^i/N
85 // See Eq (62) from Gourgoulhon et al. (2001)
86 //------------------------------------------------------------------
87
88 Vector www = hhh * gam_euler * psi4 * bsn ;
89 www.change_triad( mp.get_bvect_cart() ) ; // components on the mapping
90 // Cartesian basis
91
92 //-------------------------------------------------
93 // Constant value of W^i at the center of the star
94 //-------------------------------------------------
95
97 v_orb.set_etat_qcq() ;
98
99 for (int i=1; i<=3; i++)
100 v_orb.set(i) = www(i).val_grid_point(0, 0, 0, 0) ;
101
102 v_orb.set_triad( *(www.get_triad()) ) ;
103 v_orb.std_spectral_base() ;
104
105 //-------------------------------------------------
106 // Source and coefficients a, b for poisson_compact
107 // (independent from psi0)
108 //-------------------------------------------------
109
111 dndh_log = 0 ;
112
113 for (int l=0; l<nzet; l++) {
114
115 Param par ; // Paramater for multi-domain equation of state
116 par.add_int(l) ;
117
118 dndh_log = dndh_log + eos.der_nbar_ent(ent, 1, l, &par) ;
119
120 }
121
122 // In order to avoid any division by zero in the computation of zeta_h
123 // the value of dndh_log is set to 1 in the external domains:
124
125 for (int l=nzet; l <= nzm1; l++)
126 dndh_log.set_domain(l) = 1 ;
127
129 zeta_h.std_spectral_base() ;
130
131 Scalar Psichi = Psi % chi ;
132 Psichi.std_spectral_base() ;
133
135 lnPsi2N.std_spectral_base() ;
136
138
140 + zeta_h * lnPsi2N.derive_con(flat_spher) ;
141
142 Scalar entmb = ent - lnPsi2N ;
143
145 v_orb.change_triad(mp.get_bvect_cart()) ;
146
147 // Eq. 63 of Gourgoulhon et al. (2001)
148 Scalar source = contract(www - v_orb, 0, ent.derive_cov(flat), 0)
149 + zeta_h * ( contract(v_orb, 0, entmb.derive_cov(flat), 0)
150 + contract(www/gam_euler, 0, gam_euler.derive_cov(flat), 0) ) ;
151
152 /*
153 des_meridian(zeta_h,0., 4., "zeta_h", 10) ;
154 arrete() ;
155 des_meridian(bb(1),0., 4., "bb(1)", 10) ;
156 arrete() ;
157 des_meridian(bb(2),0., 4., "bb(2)", 10) ;
158 arrete() ;
159 des_meridian(bb(3),0., 4., "bb(3)", 10) ;
160 arrete() ;
161 des_meridian(psi0,0., 4., "psi0", 10) ;
162 arrete() ;
163 des_meridian(source,0., 4., "source", 10) ;
164 arrete() ;
165 */
166
167 source.annule(nzet, nzm1) ;
168
169 //---------------------------------------------------
170 // Resolution by means of Map_radial::poisson_compact
171 //---------------------------------------------------
172
173 Param par ;
174 int niter ;
175 par.add_int(mermax) ;
176 par.add_double(precis, 0) ;
177 par.add_double(relax, 1) ;
178 par.add_int_mod(niter) ;
179
180 if (psi0.get_etat() == ETATZERO) psi0 = 0 ;
181
184 Cmp psi0_cmp (psi0) ;
185
186 Tenseur bb_cmp(mp, 1, CON, mp.get_bvect_spher()) ;
187 bb_cmp.set_etat_qcq() ;
188 Cmp bb_cmp1 (bb(1)) ;
189 Cmp bb_cmp2 (bb(2)) ;
190 Cmp bb_cmp3 (bb(3)) ;
191 bb_cmp.set(0) = bb_cmp1 ;
192 bb_cmp.set(1) = bb_cmp2 ;
193 bb_cmp.set(2) = bb_cmp3 ;
194
195 source_cmp.va.ylm() ;
196
197 //cout << "psiO prevu" << endl << norme(psi0) << endl ;
198
200
201 psi0 = psi0_cmp ;
202
203 cout << "psiO apres" << endl << norme(psi0) << endl ;
204
205 //---------------------------------------------------
206 // Check of the solution
207 //---------------------------------------------------
208
210
212
213 source.set_spectral_va().ylm_i() ;
214
215 cout << "Check of the resolution of the continuity equation : " << endl ;
217 double erreur = 0 ;
218 for (int l=0; l<nzet; l++) {
219 double err = terr(l) ;
220 cout << " domain " << l << " : norme(source) : " << norme(source)(l)
221 << " relative error : " << err << endl ;
222 if (err > erreur) erreur = err ;
223 }
224
225 //--------------------------------
226 // Computation of grad(psi)
227 //--------------------------------
228
230
231 for (int i=1; i<=3; i++)
232 d_psi.set(i) = (psi0.derive_cov(flat))(i) + v_orb(i) ;
233
235
236 // C^1 continuation of d_psi outside the star
237 // (to ensure a smooth enthalpy field accross the stellar surface)
238 // ----------------------------------------------------------------
239
241
242 for (int i=1; i<=3; i++)
243 d_psi.set(i) = raccord_c1(d_psi(i), nzet) ;
244
245 return erreur ;
246}
247}
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...
Scalar Psi
Total conformal factor .
Definition star.h:1152
Scalar psi0
Scalar potential of the non-translational part of fluid 4-velocity (in the irrotational case)
Definition star.h:1104
Scalar psi4
Conformal factor .
Definition star.h:1158
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the star) .
Definition star.h:1177
Scalar chi
Total function .
Definition star.h:1155
Vector bsn
3-vector shift, divided by N, of the rotating coordinates, .
Definition star.h:1126
Vector d_psi
Gradient of (in the irrotational case) (Spherical components with respect to the mapping of the star...
Definition star.h:1109
double velocity_potential(int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
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
Basic array class.
Definition tbl.h:161
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
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
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
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
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition tensor.C:481
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition app_hor.h:64