LORENE
pseudo_misner.C
1 /*
2  * Copyright (c) 2003 Keisuke Taniguchi
3  *
4  * This file is part of LORENE.
5  *
6  * LORENE is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License version 2
8  * as published by the Free Software Foundation.
9  *
10  * LORENE is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with LORENE; if not, write to the Free Software
17  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18  *
19  */
20 
21 
22 
23 /*
24  * $Id: pseudo_misner.C,v 1.5 2016/12/05 16:17:46 j_novak Exp $
25  * $Log: pseudo_misner.C,v $
26  * Revision 1.5 2016/12/05 16:17:46 j_novak
27  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
28  *
29  * Revision 1.4 2014/10/13 08:52:43 j_novak
30  * Lorene classes and functions now belong to the namespace Lorene.
31  *
32  * Revision 1.3 2014/10/06 15:13:02 j_novak
33  * Modified #include directives to use c++ syntax.
34  *
35  * Revision 1.2 2007/04/24 20:13:53 f_limousin
36  * Implementation of Dirichlet and Neumann BC for the lapse
37  *
38  * Revision 1.1 2005/09/09 09:00:02 p_grandclement
39  * add pseudo_misner
40  *
41 
42  * $Header: /cvsroot/Lorene/C++/Source/Bin_ns_bh/pseudo_misner.C,v 1.5 2016/12/05 16:17:46 j_novak Exp $
43  *
44  */
45 
46 // Headers C
47 #include <cmath>
48 
49 // Headers Lorene
50 #include "bhole.h"
51 #include "nbr_spx.h"
52 #include "et_bin_nsbh.h"
53 #include "etoile.h"
54 #include "param.h"
55 #include "bin_ns_bh.h"
56 
57 #include "graphique.h"
58 #include "utilitaires.h"
59 #include "unites.h"
60 
61 namespace Lorene {
62 void Bin_ns_bh::pseudo_misner (int& ite, int itemax, double relax,
63  double precis, int bound_nn, double lim_nn) {
64 
65  using namespace Unites ;
66 
67  // Parameters for the function Map_et::poisson for n_auto
68  // ------------------------------------------------------
69  Cmp source_n_prev (star.get_mp()) ;
70  source_n_prev.set_etat_zero() ;
71 
72  Param par_poisson1 ;
73  par_poisson1.add_int(itemax, 0) ; // maximum number of iterations
74  par_poisson1.add_double(relax, 0) ; // relaxation parameter
75  par_poisson1.add_double(precis, 1) ; // required precision
76  par_poisson1.add_int_mod(ite, 0) ; // number of iterations actually used
77  par_poisson1.add_cmp_mod(source_n_prev) ;
78 
79  // Parameters for the function Map_et::poisson for confpsi_auto
80  // ------------------------------------------------------------
81  Cmp source_psi_prev (star.get_mp()) ;
82  source_psi_prev.set_etat_zero() ;
83 
84  Param par_poisson2 ;
85  par_poisson2.add_int(itemax, 0) ; // maximum number of iterations
86  par_poisson2.add_double(relax, 0) ; // relaxation parameter
87  par_poisson2.add_double(precis, 1) ; // required precision
88  par_poisson2.add_int_mod(ite, 0) ; // number of iterations actually used
89  par_poisson2.add_cmp_mod(source_psi_prev) ;
90 
91 
92  bool loop = true ;
93 
94  //=========================================================================
95  // Start of iteration
96  //=========================================================================
97  double erreur ;
98  int itere = 1 ;
99 
100  while (loop) {
101 
102  Tenseur n_auto_old (star.n_auto()) ;
103  Tenseur psi_auto_old (star.confpsi_auto()) ;
104  Tenseur n_auto_hole (hole.n_auto()) ;
105 
106  Tenseur confpsi_q = pow(star.confpsi, 4.) ;
107  Tenseur confpsi_c = pow(star.confpsi, 5.) ;
108 
109  // Lapse star
110  Tenseur source_n (qpig * star.nnn * confpsi_q * (star.ener_euler + star.s_euler)
112  source_n.set_std_base() ;
113  source_n().poisson(par_poisson1, star.n_auto.set()) ;
114  star.n_auto.set() = star.n_auto() + 0.5 ;
115  star.n_auto.set() = relax * star.n_auto() + (1-relax)*n_auto_old() ;
116  erreur = max(diffrelmax(star.n_auto(), n_auto_old())) ;
117 
118  // Psi star
119  Tenseur source_psi (-0.5 * qpig * confpsi_c * star.ener_euler) ;
120  source_psi.set_std_base() ;
121  source_psi().poisson(par_poisson2, star.confpsi_auto.set()) ;
122  star.confpsi_auto.set() = star.confpsi_auto() + 0.5 ;
123  star.confpsi_auto.set() = relax*star.confpsi_auto() + (1-relax)*psi_auto_old() ;
124 
125  // Trou noir :
126  hole.update_metric (star) ;
127  hole.solve_lapse_with_ns (relax, bound_nn, lim_nn) ;
128  //erreur = max(diffrelmax(hole.n_auto(), n_auto_hole())) ;
129 
130  hole.solve_psi_with_ns (relax) ;
131 
134 
137  star.fait_d_psi() ;
138  star.hydro_euler() ;
139 
140  cout << "Step " << itere << " " << erreur << endl ;
141 
142  if ((itere==itemax) || (erreur<precis))
143  loop = false ;
144  itere ++ ;
145  }
146 
147  ite = itere ;
148 }
149 }
const Map & get_mp() const
Returns the mapping.
Definition: etoile.h:662
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
virtual void kinematics(double omega, double x_axe)
Computes the quantities bsn and pot_centri .
double get_omega() const
Returns the orbital velocity.
Definition: bin_ns_bh.h:252
Tenseur nnn
Total lapse function.
Definition: etoile.h:512
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
Definition: etoile.h:471
Tenseur flat_scalar_prod(const Tenseur &t1, const Tenseur &t2)
Scalar product of two Tenseur when the metric is : performs the contraction of the last index of t1 w...
void update_metric_der_comp(const Bhole &comp)
Computes the derivative of metric functions related to the companion black hole.
Tenseur n_auto
Part of the lapse { N} generated principaly by the star.
Definition: et_bin_nsbh.h:85
Tenseur confpsi
Total conformal factor $$.
Definition: et_bin_nsbh.h:101
Tenseur d_confpsi_comp
Gradient of { confpsi_comp} (Cartesian components with respect to { ref_triad})
Definition: et_bin_nsbh.h:119
void solve_psi_with_ns(double relax)
Solves the equation for ~: with the condition that on the horizon, where f is the value of on th...
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:840
Et_bin_nsbh star
The neutron star.
Definition: bin_ns_bh.h:131
void fait_d_psi()
Computes the gradient of the total velocity potential .
Definition: etoile_bin.C:767
virtual void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition: etoile.C:569
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:438
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Tenseur d_confpsi_auto
Gradient of { confpsi_auto} (Cartesian components with respect to { ref_triad})
Definition: et_bin_nsbh.h:114
Tenseur d_n_auto
Gradient of { n_auto} (Cartesian components with respect to { ref_triad})
Definition: et_bin_nsbh.h:93
void update_metric(const Bhole &comp)
Computes metric coefficients from known potentials, when the companion is a black hole...
Tenseur confpsi_auto
Part of the conformal factor $$ generated principaly by the star.
Definition: et_bin_nsbh.h:104
double get_x_axe() const
Returns a constant reference to the black hole.
Definition: bin_ns_bh.h:256
void solve_lapse_with_ns(double relax, int bound_nn, double lim_nn)
Solves the equation for N ~: with the condition that N =0 on the horizon.
Definition: bhole_with_ns.C:92
Bhole hole
The black hole.
Definition: bin_ns_bh.h:134
Tenseur n_auto
Part of N generated by the hole.
Definition: bhole.h:286
Tenseur ener_euler
Total energy density in the Eulerian frame.
Definition: etoile.h:468
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition: cmp_math.C:542
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
Definition: et_bin_hydro.C:109