LORENE
star_bhns_vel_pot.C
1/*
2 * Method of class Star_bhns to compute the velocity scalar potential $\psi$
3 * by solving the continuity equation.
4 *
5 * (see file star_bhns.h for documentation).
6 *
7 */
8
9/*
10 * Copyright (c) 2005-2007 Keisuke Taniguchi
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
29char star_bhns_vel_pot_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star_bhns/star_bhns_vel_pot.C,v 1.3 2014/10/13 08:53:41 j_novak Exp $" ;
30
31/*
32 * $Id: star_bhns_vel_pot.C,v 1.3 2014/10/13 08:53:41 j_novak Exp $
33 * $Log: star_bhns_vel_pot.C,v $
34 * Revision 1.3 2014/10/13 08:53:41 j_novak
35 * Lorene classes and functions now belong to the namespace Lorene.
36 *
37 * Revision 1.2 2008/05/15 19:20:29 k_taniguchi
38 * Change of some parameters.
39 *
40 * Revision 1.1 2007/06/22 01:33:14 k_taniguchi
41 * *** empty log message ***
42 *
43 *
44 * $Header: /cvsroot/Lorene/C++/Source/Star_bhns/star_bhns_vel_pot.C,v 1.3 2014/10/13 08:53:41 j_novak Exp $
45 *
46 */
47
48// C++ headers
49//#include <>
50
51// C headers
52//#include <>
53
54// Lorene headers
55#include "star_bhns.h"
56#include "eos.h"
57#include "param.h"
58#include "cmp.h"
59#include "tenseur.h"
60#include "utilitaires.h"
61#include "unites.h"
62
63// Local prototype
64namespace Lorene {
65Cmp raccord_c1(const Cmp& uu, int l1) ;
66
67double Star_bhns::velo_pot_bhns(const double& mass_bh, const double& sepa,
68 bool kerrschild, int mermax, double precis,
69 double relax) {
70
71 // Fundamental constants and units
72 // -------------------------------
73 using namespace Unites ;
74
75 int nzm1 = mp.get_mg()->get_nzone() - 1 ;
76
77 //--------------------------------
78 // Specific relativistic enthalpy ---> hhh
79 //--------------------------------
80
81 Scalar hhh = exp(ent) ;
82 hhh.std_spectral_base() ;
83
84 //---------------------------------------------------
85 // Computation of W^i = psi4 h Gamma_n B^i/lapse_tot
86 // See Eq (62) from Gourgoulhon et al. (2001)
87 //---------------------------------------------------
88
89 Vector www = hhh * gam_euler * bsn * psi4 ;
90
91 www.change_triad( mp.get_bvect_cart() ) ; // components on the mapping
92 // Cartesian basis
93
94 //-------------------------------------------------
95 // Constant value of W^i at the center of the star
96 //-------------------------------------------------
97
99 v_orb.set_etat_qcq() ;
100
101 for (int i=1; i<=3; i++) {
102 v_orb.set(i) = www(i).val_grid_point(0, 0, 0, 0) ;
103 }
104
105 v_orb.set_triad( *(www.get_triad()) ) ;
106 v_orb.std_spectral_base() ;
107
108 //-------------------------------------------------
109 // Source and coefficients a,b for poisson_compact (independent from psi0)
110 //-------------------------------------------------
111
112 Scalar dndh_log = eos.der_nbar_ent(ent, nzet) ;
113
114 // In order to avoid any division by zero in the computation of zeta_h
115 // the value of dndh_log is set to 1 in the external domains:
116 for (int l=nzet; l<=nzm1; l++) {
117 dndh_log.set_domain(l) = 1. ;
118 }
119
121 zeta_h.std_spectral_base() ;
122
123 double erreur ;
124
125 Scalar source(mp) ;
126 Vector bb(mp, CON, mp.get_bvect_spher()) ;
128
129 if (kerrschild) {
130
131 cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
132 abort() ;
133
134 } // End of Kerr-Schild
135 else { // Isotropic coordinates with the maximal slicing
136
139 lnlappsi.std_spectral_base() ;
140
141 bb = (1. - zeta_h) * ent.derive_con(flat_spher)
142 + zeta_h * lnlappsi.derive_con(flat_spher) ;
143
144 Vector dentmb(mp, COV, mp.get_bvect_cart()) ;
145 dentmb.set_etat_qcq() ;
146 for (int i=1; i<=3; i++) {
147 dentmb.set(i) = ent.deriv(i)
150 }
151 dentmb.std_spectral_base() ;
152
153 source = contract(www - v_orb, 0, ent.derive_cov(flat), 0)
154 +zeta_h*(contract(v_orb, 0, dentmb, 0)
156 ) ;
157
158 source.annule(nzet, nzm1) ;
159
160 } // End of Isotropic coordinates
161
162
163 //----------------------------------------------------
164 // Resolution by means of Map_radial::poisson_compact
165 //----------------------------------------------------
166
167 Param par ;
168 int niter ;
169 par.add_int(mermax) ;
170 par.add_double(precis, 0) ;
171 par.add_double(relax, 1) ;
172 par.add_int_mod(niter) ;
173
174 if (psi0.get_etat() == ETATZERO) {
175 psi0 = 0. ;
176 }
177
180 Cmp psi0_cmp(psi0) ;
181
182 Tenseur bb_cmp(mp, 1, CON, mp.get_bvect_spher()) ;
183 bb_cmp.set_etat_qcq() ;
184 Cmp bb_cmp1(bb(1)) ;
185 Cmp bb_cmp2(bb(2)) ;
186 Cmp bb_cmp3(bb(3)) ;
187 bb_cmp.set(0) = bb_cmp1 ;
188 bb_cmp.set(1) = bb_cmp2 ;
189 bb_cmp.set(2) = bb_cmp3 ;
190
191 source_cmp.va.ylm() ;
192
194
195 psi0 = psi0_cmp ;
196
197 //-----------------------
198 // Check of the solution
199 //-----------------------
200
202
204
205 source.set_spectral_va().ylm_i() ;
206
207 erreur = diffrel(oper, source)(0) ;
208
209 cout << "Check of the resolution of the continuity equation : "
210 << endl ;
211 cout << " norme(source) : " << norme(source)(0)
212 << " diff oper/source : " << erreur << endl ;
213
214 //--------------------------
215 // Computation of grad(psi)
216 //--------------------------
217
218 for (int i=1; i<=3; i++)
219 d_psi.set(i) = (psi0.derive_cov(flat))(i) + v_orb(i) ;
220
221
222 // C^1 continuation of d_psi outside the star
223 // (to ensure a smooth enthalpy field accross the stellar surface)
224 // ----------------------------------------------------------------
225
227
228 for (int i=1; i<=3; i++) {
229 Cmp d_psi_i( d_psi(i) ) ;
230 d_psi_i = raccord_c1(d_psi_i, nzet) ;
231 d_psi.set(i) = d_psi_i ;
232 }
233
234 return erreur ;
235
236}
237}
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 & deriv(int i) const
Returns of *this , where .
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...
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the star ).
Definition star_bhns.h:192
double velo_pot_bhns(const double &mass_bh, const double &sepa, bool kerrschild, int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
Vector bsn
3-vector shift, divided by N , of the rotating coordinates, .
Definition star_bhns.h:99
Vector d_confo_comp
Derivative of the conformal factor generated by the companion black hole.
Definition star_bhns.h:173
Scalar confo_tot
Total conformal factor.
Definition star_bhns.h:163
Scalar psi4
Fourth power of the total conformal factor.
Definition star_bhns.h:176
Vector d_lapconf_auto
Derivative of the lapconf function generated by the star .
Definition star_bhns.h:130
Vector d_psi
Gradient of (in the irrotational case) (Spherical components with respect to the mapping of the star...
Definition star_bhns.h:82
Vector d_confo_auto
Derivative of the conformal factor generated by the star .
Definition star_bhns.h:168
Scalar lapconf_tot
Total lapconf function.
Definition star_bhns.h:119
Vector d_lapconf_comp
Derivative of the lapconf function generated by the companion black hole.
Definition star_bhns.h:135
Scalar psi0
Scalar potential of the non-translational part of fluid 4-velocity (in the irrotational case)
Definition star_bhns.h:77
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 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
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.