LORENE
single_regul.C
1/*
2 * Copyright (c) 2005 Francois Limousin
3 * Jose Luis Jaramillo
4 *
5 * This file is part of LORENE.
6 *
7 * LORENE is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * LORENE is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with LORENE; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 *
21 */
22
23
24char single_regul_C[] = "$Header: /cvsroot/Lorene/C++/Source/Isol_hor/single_regul.C,v 1.4 2014/10/13 08:53:01 j_novak Exp $" ;
25
26/*
27 * $Id: single_regul.C,v 1.4 2014/10/13 08:53:01 j_novak Exp $
28 * $Log: single_regul.C,v $
29 * Revision 1.4 2014/10/13 08:53:01 j_novak
30 * Lorene classes and functions now belong to the namespace Lorene.
31 *
32 * Revision 1.3 2014/10/06 15:13:11 j_novak
33 * Modified #include directives to use c++ syntax.
34 *
35 * Revision 1.2 2008/08/19 06:42:00 j_novak
36 * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
37 * cast-type operations, and constant strings that must be defined as const char*
38 *
39 * Revision 1.1 2007/04/13 15:28:35 f_limousin
40 * Lots of improvements, generalisation to an arbitrary state of
41 * rotation, implementation of the spatial metric given by Samaya.
42 *
43 *
44 * $Header: /cvsroot/Lorene/C++/Source/Isol_hor/single_regul.C,v 1.4 2014/10/13 08:53:01 j_novak Exp $
45 *
46 */
47
48
49//Standard
50#include <cstdlib>
51#include <cmath>
52
53//Lorene
54#include "isol_hor.h"
55#include "nbr_spx.h"
56#include "tensor.h"
57
58namespace Lorene {
59double Single_hor::regularisation (const Vector& shift_auto_temp,
60 const Vector& shift_comp_temp, double om) {
61
62 Vector shift_auto(shift_auto_temp) ;
63 shift_auto.change_triad(shift_auto.get_mp().get_bvect_cart()) ;
64 Vector shift_comp(shift_comp_temp) ;
65 shift_comp.change_triad(shift_comp.get_mp().get_bvect_cart()) ;
66 Vector shift_old (shift_auto) ;
67
68 double orientation = shift_auto.get_mp().get_rot_phi() ;
69 assert ((orientation==0) || (orientation == M_PI)) ;
70 double orientation_autre = shift_comp.get_mp().get_rot_phi() ;
71 assert ((orientation_autre==0) || (orientation_autre == M_PI)) ;
72
73 int alignes = (orientation == orientation_autre) ? 1 : -1 ;
74
75 int np = shift_auto.get_mp().get_mg()->get_np(1) ;
76 int nt = shift_auto.get_mp().get_mg()->get_nt(1) ;
77 int nr = shift_auto.get_mp().get_mg()->get_nr(1) ;
78
79 // Minimisation of the derivative of the shift on r
80 Vector shift_tot (shift_auto.get_mp(), CON, *shift_auto.get_triad()) ;
81 shift_tot.set(1).import(alignes*shift_comp(1)) ;
82 shift_tot.set(2).import(alignes*shift_comp(2)) ;
83 shift_tot.set(3).import(shift_comp(3)) ;
84
85 shift_tot = shift_tot + shift_auto ;
86
87 double indic = (orientation == 0) ? 1 : -1 ;
88
89 Vector tbi (shift_tot) ;
90 if (om != 0) {
91 for (int i=1 ; i<=3 ; i++) {
92 tbi.set(i).set_spectral_va().coef_i() ;
94 }
95
96 tbi.set(1) = *shift_tot(1).get_spectral_va().c - indic *om * shift_tot.get_mp().ya ;
97 tbi.set(2) = *shift_tot(2).get_spectral_va().c + indic *om * shift_tot.get_mp().xa ;
98 tbi.std_spectral_base() ;
99 tbi.set(1).annule_domain(nz-1) ;
100 tbi.set(2).annule_domain(nz-1) ;
101 }
102
103 Vector derive_r (shift_auto.get_mp(), CON, *shift_auto.get_triad()) ;
104 for (int i=1 ; i<=3 ; i++)
105 derive_r.set(i) = tbi(i).dsdr() ;
106
107
108 // We substract a function in order that Kij is regular
109
110 Valeur val_hor (shift_auto.get_mp().get_mg()) ;
111 Valeur fonction_radiale (shift_auto.get_mp().get_mg()) ;
112 Scalar enleve (shift_auto.get_mp()) ;
113
114 double erreur = 0 ;
115 for (int comp=1 ; comp<=3 ; comp++) {
116 val_hor.annule_hard() ;
117 for (int k=0 ; k<np ; k++)
118 for (int j=0 ; j<nt ; j++)
119 for (int i=0 ; i<nr ; i++)
120 val_hor.set(1, k, j, i) = derive_r(comp).
121 val_grid_point(1, k, j, 0) ;
122
123 double r_0 = shift_auto.get_mp().val_r (1, -1, 0, 0) ;
124 double r_1 = shift_auto.get_mp().val_r (1, 1, 0, 0) ;
125
126 fonction_radiale = pow(r_1-shift_auto.get_mp().r, 3.)*
127 (shift_auto.get_mp().r-r_0)/pow(r_1-r_0, 3.) ;
128 fonction_radiale.annule(0) ;
129 fonction_radiale.annule(2, nz-1) ;
130
131 enleve = fonction_radiale * val_hor ;
132 enleve.set_spectral_va().set_base (shift_auto(comp).
133 get_spectral_va().get_base()) ;
134
135 if (norme(enleve)(1) != 0)
136 shift_auto.set(comp) = shift_auto(comp) - enleve ;
137 if (norme(shift_auto(comp))(1) > 1e-5) {
138 Tbl diff (diffrelmax (shift_auto(comp), shift_old(comp))) ;
139 if (erreur < diff(1))
140 erreur = diff(1) ;
141 }
142 }
143
144 shift_auto.change_triad(shift_auto.get_mp().get_bvect_spher()) ;
145
146 beta_auto = shift_auto ;
147
148 return erreur ;
149}
150
151
152// Regularisation if only one black hole :
154
155 Vector shift (beta) ;
156
158 // Vector B (without boost and rotation)
159 Vector tbi (shift) ;
160
161 for (int i=1 ; i<=3 ; i++) {
162 tbi.set(i).set_spectral_va().coef_i() ;
164 }
165
166 for (int i=1 ; i<=3 ; i++)
167 shift(i).get_spectral_va().coef_i() ;
168
169 tbi.set(1) = *shift(1).get_spectral_va().c - omega*mp.y ;
170 tbi.set(2) = *shift(2).get_spectral_va().c + omega*mp.x ;
171 if (shift(3).get_etat() != ETATZERO)
172 tbi.set(3) = *shift(3).get_spectral_va().c ;
173 else
174 tbi.set(3) = 0. ;
175 tbi.std_spectral_base() ;
176
177 // We only need values at the horizon
178 tbi.set(1).annule_domain(mp.get_mg()->get_nzone()-1) ;
179 tbi.set(2).annule_domain(mp.get_mg()->get_nzone()-1) ;
180
181 Vector derive_r (mp, CON, mp.get_bvect_cart()) ;
182 derive_r.set_etat_qcq() ;
183 for (int i=1 ; i<=3 ; i++)
184 derive_r.set(i) = tbi(i).dsdr() ;
185
186 Valeur val_hor (mp.get_mg()) ;
187 Valeur fonction_radiale (mp.get_mg()) ;
188 Scalar enleve (mp) ;
189
190 double erreur = 0 ;
191 int np = mp.get_mg()->get_np(1) ;
192 int nt = mp.get_mg()->get_nt(1) ;
193 int nr = mp.get_mg()->get_nr(1) ;
194
195 double r_0 = mp.val_r(1, -1, 0, 0) ;
196 double r_1 = mp.val_r(1, 1, 0, 0) ;
197
198 for (int comp=1 ; comp<=3 ; comp++) {
199 val_hor.annule_hard() ;
200 for (int k=0 ; k<np ; k++)
201 for (int j=0 ; j<nt ; j++)
202 for (int i=0 ; i<nr ; i++)
203 val_hor.set(1, k, j, i) = derive_r(comp).val_grid_point(1, k, j, 0) ;
204
205 fonction_radiale = pow(r_1-mp.r, 3.)* (mp.r-r_0)/pow(r_1-r_0, 3.) ;
206 fonction_radiale.annule(0) ;
207 fonction_radiale.annule(2, nz-1) ;
208
209 enleve = fonction_radiale*val_hor ;
210 enleve.set_spectral_va().base = shift(comp).get_spectral_va().base ;
211
212 Scalar copie (shift(comp)) ;
213 shift.set(comp) = shift(comp)-enleve ;
214 shift.std_spectral_base() ;
215
216 assert (shift(comp).check_dzpuis(0)) ;
217
218 // Intensity of the correction (if nonzero)
219 Tbl norm (norme(shift(comp))) ;
220 if (norm(1) > 1e-5) {
221 Tbl diff (diffrelmax (copie, shift(comp))) ;
222 if (erreur<diff(1))
223 erreur = diff(1) ;
224 }
225 }
226
228 beta = shift ;
229
230 return erreur ;
231}
232}
virtual double val_r(int l, double xi, double theta, double pphi) const
Returns the value of the radial coordinate r for a given in a given domain.
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
Coord y
y coordinate centered on the grid
Definition map.h:727
Coord ya
Absolute y coordinate.
Definition map.h:731
Coord r
r coordinate centered on the grid
Definition map.h:718
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
Coord x
x coordinate centered on the grid
Definition map.h:726
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition map.h:775
Coord xa
Absolute x coordinate.
Definition map.h:730
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
void import(const Scalar &ci)
Assignment to another Scalar defined on a different mapping.
Vector beta_auto
Shift function .
Definition isol_hor.h:944
Map_af & mp
Affine mapping.
Definition isol_hor.h:900
int nz
Number of zones.
Definition isol_hor.h:903
double regularisation(const Vector &shift_auto, const Vector &shift_comp, double ang_vel)
Corrects shift_auto in such a way that the total is equal to zero in the horizon,...
double omega
Angular velocity in LORENE's units.
Definition isol_hor.h:909
double regularise_one()
Corrects the shift in the innermost shell, so that it remains and that equals zero on the horizon.
Vector beta
Shift function .
Definition isol_hor.h:950
Basic array class.
Definition tbl.h:161
Values and coefficients of a (real-value) function.
Definition valeur.h:287
void set_etat_c_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl c ).
Definition valeur.C:701
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition valeur.C:810
void annule(int l)
Sets the Valeur to zero in a given domain.
Definition valeur.C:744
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition valeur.h:363
void coef_i() const
Computes the physical value of *this.
Base_val base
Bases on which the spectral expansion is performed.
Definition valeur.h:305
void annule_hard()
Sets the Valeur to zero in a hard way.
Definition valeur.C:723
Tensor field of valence 1.
Definition vector.h:188
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:316
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
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 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
const Map & get_mp() const
Returns the mapping.
Definition tensor.h:861
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tensor.h:866
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition tensor.C:481
Lorene prototypes.
Definition app_hor.h:64