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 
27 char 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
67 namespace Lorene {
68 Cmp 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 
96  Vector v_orb(mp, CON, mp.get_bvect_cart()) ;
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 
110  Scalar dndh_log(mp) ;
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 
128  Scalar zeta_h( ent / dndh_log ) ;
129  zeta_h.std_spectral_base() ;
130 
131  Scalar Psichi = Psi % chi ;
132  Psichi.std_spectral_base() ;
133 
134  Scalar lnPsi2N = log(Psichi) ;
135  lnPsi2N.std_spectral_base() ;
136 
137  Metric_flat flat_spher (mp.flat_met_spher()) ;
138 
139  Vector bb = (1 - zeta_h) * ent.derive_con(flat_spher)
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 
182  Cmp source_cmp (source) ;
183  Cmp zeta_h_cmp (zeta_h) ;
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 
199  mp.poisson_compact(nzet, source_cmp, zeta_h_cmp, bb_cmp, par, psi0_cmp ) ;
200 
201  psi0 = psi0_cmp ;
202 
203  cout << "psiO apres" << endl << norme(psi0) << endl ;
204 
205  //---------------------------------------------------
206  // Check of the solution
207  //---------------------------------------------------
208 
209  Scalar bb_dpsi0 = contract(bb, 0, psi0.derive_cov(flat_spher), 0) ;
210 
211  Scalar oper = zeta_h * psi0.laplacian() + bb_dpsi0 ;
212 
213  source.set_spectral_va().ylm_i() ;
214 
215  cout << "Check of the resolution of the continuity equation : " << endl ;
216  Tbl terr = diffrel(oper, source) ;
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 
229  d_psi.set_etat_qcq() ;
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 
240  d_psi.annule(nzet, nzm1) ;
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 }
Lorene::Tenseur
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:301
Lorene::Star::gam_euler
Scalar gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition: star.h:204
Lorene::Tenseur::set
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
Lorene::log
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:296
Lorene::Param::add_double
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:315
Lorene::exp
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
Lorene::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Star_bin_xcts::chi
Scalar chi
Total function .
Definition: star.h:1155
Lorene::Star::nzet
int nzet
Number of domains of *mp occupied by the star.
Definition: star.h:183
Lorene::Cmp::va
Valeur va
The numerical value of the Cmp
Definition: cmp.h:464
Lorene::Map::get_bvect_spher
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:783
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Map::flat_met_spher
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
Lorene::Tensor::annule
virtual void annule(int l_min, int l_max)
Sets the Tensor to zero in several domains.
Definition: tensor.C:671
Lorene::Scalar::derive_cov
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
Lorene::norme
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:481
Lorene::Map::poisson_compact
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...
Lorene::Metric_flat
Flat metric for tensor calculation.
Definition: metric.h:261
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::Star_bin_xcts::psi0
Scalar psi0
Scalar potential of the non-translational part of fluid 4-velocity (in the irrotational case)
Definition: star.h:1104
Lorene::Vector::change_triad
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition: vector_change_triad.C:75
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
Lorene::Tensor::set_etat_qcq
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition: tensor.C:481
Lorene::Star_bin_xcts::psi4
Scalar psi4
Conformal factor .
Definition: star.h:1158
Lorene::Star::mp
Map & mp
Mapping associated with the star.
Definition: star.h:180
Lorene::Eos::der_nbar_ent
Cmp der_nbar_ent(const Cmp &ent, int nzet, int l_min=0, const Param *par=0x0) const
Computes the logarithmic derivative from the log-enthalpy and extra parameters.
Definition: eos.C:407
Lorene::Vector::set
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
Lorene::Star::ent
Scalar ent
Log-enthalpy.
Definition: star.h:190
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
Lorene::Tbl
Basic array class.
Definition: tbl.h:161
Lorene::Scalar::laplacian
const Scalar & laplacian(int ced_mult_r=4) const
Returns the Laplacian of *this.
Definition: scalar_deriv.C:436
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Star_bin_xcts::bsn
Vector bsn
3-vector shift, divided by N, of the rotating coordinates, .
Definition: star.h:1126
Lorene::Param::add_int_mod
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition: param.C:385
Lorene::Map::get_bvect_cart
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
Lorene::Scalar::annule
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition: scalar.C:391
Lorene::Vector
Tensor field of valence 1.
Definition: vector.h:188
Lorene::Star_bin_xcts::flat
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the star) .
Definition: star.h:1177
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
Lorene::Tenseur::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:636
Lorene::Scalar::derive_con
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
Lorene::Tensor::get_triad
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tensor.h:866
Lorene::Scalar::set_domain
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Definition: scalar.h:615
Lorene::Star_bin_xcts::Psi
Scalar Psi
Total conformal factor .
Definition: star.h:1152
Lorene::Scalar::std_spectral_base
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
Lorene::Param
Parameter storage.
Definition: param.h:125
Lorene::Star::eos
const Eos & eos
Equation of state of the stellar matter.
Definition: star.h:185
Lorene::Valeur::ylm
void ylm()
Computes the coefficients of *this.
Definition: valeur_ylm.C:138
Lorene::Tensor::set_triad
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition: tensor.C:519
Lorene::Star_bin_xcts::d_psi
Vector d_psi
Gradient of (in the irrotational case) (Spherical components with respect to the mapping of the star...
Definition: star.h:1109
Lorene::Param::add_int
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition: param.C:246
Lorene::Valeur::ylm_i
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:131
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition: tenseur_operateur.C:279
Lorene::Star_bin_xcts::velocity_potential
double velocity_potential(int mermax, double precis, double relax)
Computes the non-translational part of the velocity scalar potential by solving the continuity equat...
Definition: star_bin_vel_pot_xcts.C:70
Lorene::Vector::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: vector.C:316