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 
24 
25 
26 /*
27  * $Id: single_regul.C,v 1.5 2016/12/05 16:17:56 j_novak Exp $
28  * $Log: single_regul.C,v $
29  * Revision 1.5 2016/12/05 16:17:56 j_novak
30  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
31  *
32  * Revision 1.4 2014/10/13 08:53:01 j_novak
33  * Lorene classes and functions now belong to the namespace Lorene.
34  *
35  * Revision 1.3 2014/10/06 15:13:11 j_novak
36  * Modified #include directives to use c++ syntax.
37  *
38  * Revision 1.2 2008/08/19 06:42:00 j_novak
39  * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
40  * cast-type operations, and constant strings that must be defined as const char*
41  *
42  * Revision 1.1 2007/04/13 15:28:35 f_limousin
43  * Lots of improvements, generalisation to an arbitrary state of
44  * rotation, implementation of the spatial metric given by Samaya.
45  *
46  *
47  * $Header: /cvsroot/Lorene/C++/Source/Isol_hor/single_regul.C,v 1.5 2016/12/05 16:17:56 j_novak Exp $
48  *
49  */
50 
51 
52 //Standard
53 #include <cstdlib>
54 #include <cmath>
55 
56 //Lorene
57 #include "isol_hor.h"
58 #include "nbr_spx.h"
59 #include "tensor.h"
60 
61 namespace Lorene {
62 double Single_hor::regularisation (const Vector& shift_auto_temp,
63  const Vector& shift_comp_temp, double om) {
64 
65  Vector shift_auto(shift_auto_temp) ;
66  shift_auto.change_triad(shift_auto.get_mp().get_bvect_cart()) ;
67  Vector shift_comp(shift_comp_temp) ;
68  shift_comp.change_triad(shift_comp.get_mp().get_bvect_cart()) ;
69  Vector shift_old (shift_auto) ;
70 
71  double orientation = shift_auto.get_mp().get_rot_phi() ;
72  assert ((orientation==0) || (orientation == M_PI)) ;
73  double orientation_autre = shift_comp.get_mp().get_rot_phi() ;
74  assert ((orientation_autre==0) || (orientation_autre == M_PI)) ;
75 
76  int alignes = (orientation == orientation_autre) ? 1 : -1 ;
77 
78  int np = shift_auto.get_mp().get_mg()->get_np(1) ;
79  int nt = shift_auto.get_mp().get_mg()->get_nt(1) ;
80  int nr = shift_auto.get_mp().get_mg()->get_nr(1) ;
81 
82  // Minimisation of the derivative of the shift on r
83  Vector shift_tot (shift_auto.get_mp(), CON, *shift_auto.get_triad()) ;
84  shift_tot.set(1).import(alignes*shift_comp(1)) ;
85  shift_tot.set(2).import(alignes*shift_comp(2)) ;
86  shift_tot.set(3).import(shift_comp(3)) ;
87 
88  shift_tot = shift_tot + shift_auto ;
89 
90  double indic = (orientation == 0) ? 1 : -1 ;
91 
92  Vector tbi (shift_tot) ;
93  if (om != 0) {
94  for (int i=1 ; i<=3 ; i++) {
95  tbi.set(i).set_spectral_va().coef_i() ;
97  }
98 
99  tbi.set(1) = *shift_tot(1).get_spectral_va().c - indic *om * shift_tot.get_mp().ya ;
100  tbi.set(2) = *shift_tot(2).get_spectral_va().c + indic *om * shift_tot.get_mp().xa ;
101  tbi.std_spectral_base() ;
102  tbi.set(1).annule_domain(nz-1) ;
103  tbi.set(2).annule_domain(nz-1) ;
104  }
105 
106  Vector derive_r (shift_auto.get_mp(), CON, *shift_auto.get_triad()) ;
107  for (int i=1 ; i<=3 ; i++)
108  derive_r.set(i) = tbi(i).dsdr() ;
109 
110 
111  // We substract a function in order that Kij is regular
112 
113  Valeur val_hor (shift_auto.get_mp().get_mg()) ;
114  Valeur fonction_radiale (shift_auto.get_mp().get_mg()) ;
115  Scalar enleve (shift_auto.get_mp()) ;
116 
117  double erreur = 0 ;
118  for (int comp=1 ; comp<=3 ; comp++) {
119  val_hor.annule_hard() ;
120  for (int k=0 ; k<np ; k++)
121  for (int j=0 ; j<nt ; j++)
122  for (int i=0 ; i<nr ; i++)
123  val_hor.set(1, k, j, i) = derive_r(comp).
124  val_grid_point(1, k, j, 0) ;
125 
126  double r_0 = shift_auto.get_mp().val_r (1, -1, 0, 0) ;
127  double r_1 = shift_auto.get_mp().val_r (1, 1, 0, 0) ;
128 
129  fonction_radiale = pow(r_1-shift_auto.get_mp().r, 3.)*
130  (shift_auto.get_mp().r-r_0)/pow(r_1-r_0, 3.) ;
131  fonction_radiale.annule(0) ;
132  fonction_radiale.annule(2, nz-1) ;
133 
134  enleve = fonction_radiale * val_hor ;
135  enleve.set_spectral_va().set_base (shift_auto(comp).
136  get_spectral_va().get_base()) ;
137 
138  if (norme(enleve)(1) != 0)
139  shift_auto.set(comp) = shift_auto(comp) - enleve ;
140  if (norme(shift_auto(comp))(1) > 1e-5) {
141  Tbl diff (diffrelmax (shift_auto(comp), shift_old(comp))) ;
142  if (erreur < diff(1))
143  erreur = diff(1) ;
144  }
145  }
146 
147  shift_auto.change_triad(shift_auto.get_mp().get_bvect_spher()) ;
148 
149  beta_auto = shift_auto ;
150 
151  return erreur ;
152 }
153 
154 
155 // Regularisation if only one black hole :
157 
158  Vector shift (beta) ;
159 
160  shift.change_triad(mp.get_bvect_cart()) ;
161  // Vector B (without boost and rotation)
162  Vector tbi (shift) ;
163 
164  for (int i=1 ; i<=3 ; i++) {
165  tbi.set(i).set_spectral_va().coef_i() ;
166  tbi.set(i).set_spectral_va().set_etat_c_qcq() ;
167  }
168 
169  for (int i=1 ; i<=3 ; i++)
170  shift(i).get_spectral_va().coef_i() ;
171 
172  tbi.set(1) = *shift(1).get_spectral_va().c - omega*mp.y ;
173  tbi.set(2) = *shift(2).get_spectral_va().c + omega*mp.x ;
174  if (shift(3).get_etat() != ETATZERO)
175  tbi.set(3) = *shift(3).get_spectral_va().c ;
176  else
177  tbi.set(3) = 0. ;
178  tbi.std_spectral_base() ;
179 
180  // We only need values at the horizon
181  tbi.set(1).annule_domain(mp.get_mg()->get_nzone()-1) ;
182  tbi.set(2).annule_domain(mp.get_mg()->get_nzone()-1) ;
183 
184  Vector derive_r (mp, CON, mp.get_bvect_cart()) ;
185  derive_r.set_etat_qcq() ;
186  for (int i=1 ; i<=3 ; i++)
187  derive_r.set(i) = tbi(i).dsdr() ;
188 
189  Valeur val_hor (mp.get_mg()) ;
190  Valeur fonction_radiale (mp.get_mg()) ;
191  Scalar enleve (mp) ;
192 
193  double erreur = 0 ;
194  int np = mp.get_mg()->get_np(1) ;
195  int nt = mp.get_mg()->get_nt(1) ;
196  int nr = mp.get_mg()->get_nr(1) ;
197 
198  double r_0 = mp.val_r(1, -1, 0, 0) ;
199  double r_1 = mp.val_r(1, 1, 0, 0) ;
200 
201  for (int comp=1 ; comp<=3 ; comp++) {
202  val_hor.annule_hard() ;
203  for (int k=0 ; k<np ; k++)
204  for (int j=0 ; j<nt ; j++)
205  for (int i=0 ; i<nr ; i++)
206  val_hor.set(1, k, j, i) = derive_r(comp).val_grid_point(1, k, j, 0) ;
207 
208  fonction_radiale = pow(r_1-mp.r, 3.)* (mp.r-r_0)/pow(r_1-r_0, 3.) ;
209  fonction_radiale.annule(0) ;
210  fonction_radiale.annule(2, nz-1) ;
211 
212  enleve = fonction_radiale*val_hor ;
213  enleve.set_spectral_va().base = shift(comp).get_spectral_va().base ;
214 
215  Scalar copie (shift(comp)) ;
216  shift.set(comp) = shift(comp)-enleve ;
217  shift.std_spectral_base() ;
218 
219  assert (shift(comp).check_dzpuis(0)) ;
220 
221  // Intensity of the correction (if nonzero)
222  Tbl norm (norme(shift(comp))) ;
223  if (norm(1) > 1e-5) {
224  Tbl diff (diffrelmax (copie, shift(comp))) ;
225  if (erreur<diff(1))
226  erreur = diff(1) ;
227  }
228  }
229 
230  shift.change_triad(mp.get_bvect_spher()) ;
231  beta = shift ;
232 
233  return erreur ;
234 }
235 }
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:675
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition: tensor.C:490
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:351
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:801
Lorene prototypes.
Definition: app_hor.h:67
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
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, which should ensure the regularity of .
Definition: single_regul.C:62
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
void coef_i() const
Computes the physical value of *this.
double omega
Angular velocity in LORENE&#39;s units.
Definition: isol_hor.h:909
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.
Definition: map_af_radius.C:99
Values and coefficients of a (real-value) function.
Definition: valeur.h:297
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition: map.h:793
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:322
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:386
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:484
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tensor.h:879
Mtbl * c
Values of the function at the points of the multi-grid.
Definition: valeur.h:309
Vector beta_auto
Shift function .
Definition: isol_hor.h:944
Vector beta
Shift function .
Definition: isol_hor.h:950
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
double regularise_one()
Corrects the shift in the innermost shell, so that it remains and that equals zero on the horizon...
Definition: single_regul.C:156
void import(const Scalar &ci)
Assignment to another Scalar defined on a different mapping.
Definition: scalar_import.C:71
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:809
Map_af & mp
Affine mapping.
Definition: isol_hor.h:900
Coord y
y coordinate centered on the grid
Definition: map.h:745
Coord x
x coordinate centered on the grid
Definition: map.h:744
void set_etat_c_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl c )...
Definition: valeur.C:704
Basic array class.
Definition: tbl.h:164
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:610
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:302
int nz
Number of zones.
Definition: isol_hor.h:903
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:874
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition: cmp_math.C:542
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:607
Coord r
r coordinate centered on the grid
Definition: map.h:736