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
21char pseudo_misner_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bin_ns_bh/pseudo_misner.C,v 1.4 2014/10/13 08:52:43 j_novak Exp $" ;
22
23/*
24 * $Id: pseudo_misner.C,v 1.4 2014/10/13 08:52:43 j_novak Exp $
25 * $Log: pseudo_misner.C,v $
26 * Revision 1.4 2014/10/13 08:52:43 j_novak
27 * Lorene classes and functions now belong to the namespace Lorene.
28 *
29 * Revision 1.3 2014/10/06 15:13:02 j_novak
30 * Modified #include directives to use c++ syntax.
31 *
32 * Revision 1.2 2007/04/24 20:13:53 f_limousin
33 * Implementation of Dirichlet and Neumann BC for the lapse
34 *
35 * Revision 1.1 2005/09/09 09:00:02 p_grandclement
36 * add pseudo_misner
37 *
38
39 * $Header: /cvsroot/Lorene/C++/Source/Bin_ns_bh/pseudo_misner.C,v 1.4 2014/10/13 08:52:43 j_novak Exp $
40 *
41 */
42
43// Headers C
44#include <cmath>
45
46// Headers Lorene
47#include "bhole.h"
48#include "nbr_spx.h"
49#include "et_bin_nsbh.h"
50#include "etoile.h"
51#include "param.h"
52#include "bin_ns_bh.h"
53
54#include "graphique.h"
55#include "utilitaires.h"
56#include "unites.h"
57
58namespace Lorene {
59void Bin_ns_bh::pseudo_misner (int& ite, int itemax, double relax,
60 double precis, int bound_nn, double lim_nn) {
61
62 using namespace Unites ;
63
64 // Parameters for the function Map_et::poisson for n_auto
65 // ------------------------------------------------------
66 Cmp source_n_prev (star.get_mp()) ;
67 source_n_prev.set_etat_zero() ;
68
69 Param par_poisson1 ;
70 par_poisson1.add_int(itemax, 0) ; // maximum number of iterations
71 par_poisson1.add_double(relax, 0) ; // relaxation parameter
72 par_poisson1.add_double(precis, 1) ; // required precision
73 par_poisson1.add_int_mod(ite, 0) ; // number of iterations actually used
74 par_poisson1.add_cmp_mod(source_n_prev) ;
75
76 // Parameters for the function Map_et::poisson for confpsi_auto
77 // ------------------------------------------------------------
78 Cmp source_psi_prev (star.get_mp()) ;
79 source_psi_prev.set_etat_zero() ;
80
81 Param par_poisson2 ;
82 par_poisson2.add_int(itemax, 0) ; // maximum number of iterations
83 par_poisson2.add_double(relax, 0) ; // relaxation parameter
84 par_poisson2.add_double(precis, 1) ; // required precision
85 par_poisson2.add_int_mod(ite, 0) ; // number of iterations actually used
86 par_poisson2.add_cmp_mod(source_psi_prev) ;
87
88
89 bool loop = true ;
90
91 //=========================================================================
92 // Start of iteration
93 //=========================================================================
94 double erreur ;
95 int itere = 1 ;
96
97 while (loop) {
98
99 Tenseur n_auto_old (star.n_auto()) ;
100 Tenseur psi_auto_old (star.confpsi_auto()) ;
101 Tenseur n_auto_hole (hole.n_auto()) ;
102
103 Tenseur confpsi_q = pow(star.confpsi, 4.) ;
104 Tenseur confpsi_c = pow(star.confpsi, 5.) ;
105
106 // Lapse star
107 Tenseur source_n (qpig * star.nnn * confpsi_q * (star.ener_euler + star.s_euler)
109 source_n.set_std_base() ;
110 source_n().poisson(par_poisson1, star.n_auto.set()) ;
111 star.n_auto.set() = star.n_auto() + 0.5 ;
112 star.n_auto.set() = relax * star.n_auto() + (1-relax)*n_auto_old() ;
113 erreur = max(diffrelmax(star.n_auto(), n_auto_old())) ;
114
115 // Psi star
116 Tenseur source_psi (-0.5 * qpig * confpsi_c * star.ener_euler) ;
117 source_psi.set_std_base() ;
118 source_psi().poisson(par_poisson2, star.confpsi_auto.set()) ;
120 star.confpsi_auto.set() = relax*star.confpsi_auto() + (1-relax)*psi_auto_old() ;
121
122 // Trou noir :
123 hole.update_metric (star) ;
124 hole.solve_lapse_with_ns (relax, bound_nn, lim_nn) ;
125 //erreur = max(diffrelmax(hole.n_auto(), n_auto_hole())) ;
126
127 hole.solve_psi_with_ns (relax) ;
128
131
134 star.fait_d_psi() ;
135 star.hydro_euler() ;
136
137 cout << "Step " << itere << " " << erreur << endl ;
138
139 if ((itere==itemax) || (erreur<precis))
140 loop = false ;
141 itere ++ ;
142 }
143
144 ite = itere ;
145}
146}
Tenseur n_auto
Part of N generated by the hole.
Definition bhole.h:286
void solve_psi_with_ns(double relax)
Solves the equation for ~:
void solve_lapse_with_ns(double relax, int bound_nn, double lim_nn)
Solves the equation for N ~:
double get_omega() const
Returns the orbital velocity.
Definition bin_ns_bh.h:249
Et_bin_nsbh star
The neutron star.
Definition bin_ns_bh.h:128
Bhole hole
The black hole.
Definition bin_ns_bh.h:131
double get_x_axe() const
Returns a constant reference to the black hole.
Definition bin_ns_bh.h:253
Tenseur confpsi_auto
Part of the conformal factor $\Psi$ generated principaly by the star.
Tenseur d_confpsi_comp
Gradient of {\tt confpsi_comp} (Cartesian components with respect to {\tt ref_triad})
void update_metric(const Bhole &comp)
Computes metric coefficients from known potentials, when the companion is a black hole.
Tenseur d_n_auto
Gradient of {\tt n_auto} (Cartesian components with respect to {\tt ref_triad})
Definition et_bin_nsbh.h:93
virtual void kinematics(double omega, double x_axe)
Computes the quantities bsn and pot_centri .
Tenseur n_auto
Part of the lapse {\it N} generated principaly by the star.
Definition et_bin_nsbh.h:85
void update_metric_der_comp(const Bhole &comp)
Computes the derivative of metric functions related to the companion black hole.
Tenseur confpsi
Total conformal factor $\Psi$.
Tenseur d_confpsi_auto
Gradient of {\tt confpsi_auto} (Cartesian components with respect to {\tt ref_triad})
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
void fait_d_psi()
Computes the gradient of the total velocity potential .
Definition etoile_bin.C:764
Tenseur nnn
Total lapse function.
Definition etoile.h:509
const Map & get_mp() const
Returns the mapping.
Definition etoile.h:659
virtual void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition etoile.C:566
Tenseur ener_euler
Total energy density in the Eulerian frame.
Definition etoile.h:465
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
Definition etoile.h:468
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition cmp_math.C:539
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...
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.