LORENE
star_bin_upmetr.C
1 /*
2  * Methods of Star_bin::update_metric
3  *
4  * (see file star.h for documentation)
5  *
6  */
7 
8 /*
9  * Copyright (c) 2004 Francois Limousin
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License as published by
15  * the Free Software Foundation; either version 2 of the License, or
16  * (at your option) any later version.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 char star_binupmetr_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star/star_bin_upmetr.C,v 1.14 2014/10/13 08:53:38 j_novak Exp $" ;
30 
31 /*
32  * $Id: star_bin_upmetr.C,v 1.14 2014/10/13 08:53:38 j_novak Exp $
33  * $Log: star_bin_upmetr.C,v $
34  * Revision 1.14 2014/10/13 08:53:38 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.13 2005/09/13 19:38:31 f_limousin
38  * Reintroduction of the resolution of the equations in cartesian coordinates.
39  *
40  * Revision 1.12 2005/02/24 16:05:14 f_limousin
41  * Change the name of some variables (for instance dcov_logn --> dlogn).
42  *
43  * Revision 1.11 2005/02/18 13:14:18 j_novak
44  * Changing of malloc/free to new/delete + suppression of some unused variables
45  * (trying to avoid compilation warnings).
46  *
47  * Revision 1.10 2005/02/17 17:34:10 f_limousin
48  * Change the name of some quantities to be consistent with other classes
49  * (for instance nnn is changed to nn, shift to beta, beta to lnq...)
50  *
51  * Revision 1.9 2004/06/22 12:52:26 f_limousin
52  * Change qq, qq_auto and qq_comp to beta, beta_auto and beta_comp.
53  *
54  * Revision 1.8 2004/06/07 16:23:52 f_limousin
55  * New treatment for conformally flat metrics.
56  *
57  * Revision 1.7 2004/04/08 16:33:16 f_limousin
58  * The new variable is ln(Q) instead of Q=psi^2*N. It improves the
59  * convergence of the code.
60  *
61  * Revision 1.6 2004/03/23 09:58:55 f_limousin
62  * Add function Star::update_decouple()
63  *
64  * Revision 1.5 2004/02/27 10:54:27 f_limousin
65  * To avoid errors when merging versions of Lorene.
66  *
67  * Revision 1.4 2004/02/27 09:56:42 f_limousin
68  * Many minor changes.
69  *
70  * Revision 1.3 2004/02/21 17:05:13 e_gourgoulhon
71  * Method Scalar::point renamed Scalar::val_grid_point.
72  * Method Scalar::set_point renamed Scalar::set_grid_point.
73  *
74  * Revision 1.2 2004/01/20 15:20:08 f_limousin
75  * First version
76  *
77  *
78  * $Header: /cvsroot/Lorene/C++/Source/Star/star_bin_upmetr.C,v 1.14 2014/10/13 08:53:38 j_novak Exp $ *
79  */
80 
81 
82 // Headers Lorene
83 #include "cmp.h"
84 #include "star.h"
85 #include "graphique.h"
86 #include "utilitaires.h"
87 
88 //----------------------------------//
89 // Version without relaxation //
90 //----------------------------------//
91 
92 namespace Lorene {
93 void Star_bin::update_metric(const Star_bin& comp, double om) {
94 
95  // Computation of quantities coming from the companion
96  // ---------------------------------------------------
97 
98  int nz = mp.get_mg()->get_nzone() ;
99  int nr = mp.get_mg()->get_nr(0);
100  int nt = mp.get_mg()->get_nt(0);
101  int np = mp.get_mg()->get_np(0);
102 
103  const Map& mp_comp (comp.get_mp()) ;
104 
105  if ( (comp.logn_auto).get_etat() == ETATZERO ) {
107  }
108  else{
110  logn_comp.import( comp.logn_auto ) ;
112  }
113 
114 
117 
118  Vector comp_beta(comp.beta_auto) ;
119  comp_beta.change_triad(mp_comp.get_bvect_cart()) ;
120  comp_beta.change_triad(mp.get_bvect_cart()) ;
121 
122  assert ( *(beta_comp.get_triad()) == *(comp_beta.get_triad())) ;
123 
124  (beta_comp.set(1)).import( comp_beta(1) ) ;
125  (beta_comp.set(2)).import( comp_beta(2) ) ;
126  (beta_comp.set(3)).import( comp_beta(3) ) ;
127 
129 
130  if ( (comp.lnq_auto).get_etat() == ETATZERO ) {
132  }
133  else{
135  lnq_comp.import( comp.lnq_auto ) ;
137  }
138 
139 
141  Sym_tensor comp_hij(comp.hij_auto) ;
142  comp_hij.change_triad(mp_comp.get_bvect_cart()) ;
143  comp_hij.change_triad(mp.get_bvect_cart()) ;
144 
145  assert ( *(hij_comp.get_triad()) == *(comp_hij.get_triad())) ;
146 
147  for(int i=1; i<=3; i++)
148  for(int j=i; j<=3; j++) {
149 
150  hij_comp.set(i,j).set_etat_qcq() ;
151  hij_comp.set(i,j).import( (comp_hij)(i,j) ) ;
152  }
153 
154  hij_comp.std_spectral_base() ; // set the bases for spectral expansions
155 
156 // Lapse function N
157 // ----------------
158 
159  logn = logn_auto + logn_comp ;
160 
161  nn = exp( logn ) ;
162 
163  nn.std_spectral_base() ; // set the bases for spectral expansions
164 
165 // Quantity lnq = log(psi^2*N)
166 // ----------------------
167 
168  lnq = lnq_auto + lnq_comp ;
169 
170  psi4 = exp(2*lnq) / (nn * nn) ;
172 
173 // Shift vector
174 // -------------
175 
176  beta = beta_auto + beta_comp ;
177 
178 // Coefficients of the 3-metric tilde
179 // ----------------------------------
180 
181  Sym_tensor gtilde_con (mp, CON, mp.get_bvect_cart()) ;
182 
183  for(int i=1; i<=3; i++)
184  for(int j=i; j<=3; j++) {
185 
186  hij.set(i,j) = hij_auto(i,j) + hij_comp(i,j) ;
187  gtilde_con.set(i,j) = hij(i,j) + flat.con()(i,j) ;
188  }
189 
190  gtilde = gtilde_con ;
191 
192  Sym_tensor tens_gamma = gtilde_con / psi4 ;
193  gamma = tens_gamma ;
194 
195 
196  // For conformally flat metrics
197  // ----------------------------
198 
199  if (conf_flat) {
202  hij.set_etat_zero() ;
203  gtilde = flat ;
204  tens_gamma = flat.con() / psi4 ;
205  gamma = tens_gamma ;
206  }
207 
208 
209  // Determinant of gtilde
210 
211  Scalar det_gtilde (gtilde.determinant()) ;
212 
213  double* max_det = new double[nz] ;
214  double* min_det = new double[nz] ;
215  double* moy_det = new double[nz] ;
216 
217  for (int i=0; i<nz; i++){
218  min_det[i] = 2 ;
219  moy_det[i] = 0 ;
220  max_det[i] = 0 ;
221  }
222 
223  for (int l=0; l<nz; l++)
224  for (int k=0; k<np; k++)
225  for (int j=0; j<nt; j++)
226  for (int i=0; i<nr; i++){
227 
228  moy_det[l] = moy_det[l] + det_gtilde.val_grid_point(l,k,j,i) ;
229  if (det_gtilde.val_grid_point(l,k,j,i) > max_det[l]){
230  max_det[l] = det_gtilde.val_grid_point(l,k,j,i) ;
231  }
232  if (det_gtilde.val_grid_point(l,k,j,i) < min_det[l]){
233  min_det[l] = det_gtilde.val_grid_point(l,k,j,i) ;
234  }
235  }
236 
237  cout << "average determinant of gtilde in each zone : " << endl ;
238  for (int l=0; l<nz; l++){
239  cout << moy_det[l]/(nr*nt*np) << " " ;
240  }
241  cout << endl ;
242 
243 
244  cout << "maximum of the determinant of gtilde in each zone : " << endl ;
245  for (int l=0; l<nz; l++){
246  cout << max_det[l] << " " ;
247  }
248  cout << endl ;
249 
250  cout << "minimum of the determinant of gtilde in each zone : " << endl ;
251  for (int l=0; l<nz; l++){
252  cout << min_det[l] << " " ;
253  }
254  cout << endl ;
255 
256  // ... extrinsic curvature (tkij_auto and kcar_auto)
257  extrinsic_curvature(om) ;
258 
259 
260 // The derived quantities are obsolete
261 // -----------------------------------
262 
263  del_deriv() ;
264 
265  delete max_det ;
266  delete moy_det ;
267  delete min_det ;
268 }
269 
270 
271 
272 //----------------------------------//
273 // Version with relaxation //
274 //----------------------------------//
275 
277  const Star_bin& star_jm1,
278  double relax, double om) {
279 
280 
281  // Computation of quantities coming from the companion
282  // ---------------------------------------------------
283 
284  int nz = mp.get_mg()->get_nzone() ;
285  int nr = mp.get_mg()->get_nr(0);
286  int nt = mp.get_mg()->get_nt(0);
287  int np = mp.get_mg()->get_np(0);
288 
289  const Map& mp_comp (comp.get_mp()) ;
290 
291  if ( (comp.logn_auto).get_etat() == ETATZERO ) {
293  }
294  else{
296  logn_comp.import( comp.logn_auto ) ;
298  }
299 
300 
303 
304  Vector comp_beta(comp.beta_auto) ;
305  comp_beta.change_triad(mp_comp.get_bvect_cart()) ;
306  comp_beta.change_triad(mp.get_bvect_cart()) ;
307 
308  assert ( *(beta_comp.get_triad()) == *(comp_beta.get_triad())) ;
309 
310  (beta_comp.set(1)).import( comp_beta(1) ) ;
311  (beta_comp.set(2)).import( comp_beta(2) ) ;
312  (beta_comp.set(3)).import( comp_beta(3) ) ;
313 
315 
316 
317  if ( (comp.lnq_auto).get_etat() == ETATZERO ) {
319  }
320  else{
322  lnq_comp.import( comp.lnq_auto ) ;
324  }
325 
327 
328  Sym_tensor comp_hij(comp.hij_auto) ;
329  comp_hij.change_triad(mp_comp.get_bvect_cart()) ;
330  comp_hij.change_triad(mp.get_bvect_cart()) ;
331 
332  assert ( *(hij_comp.get_triad()) == *(comp_hij.get_triad())) ;
333 
334 
335  for(int i=1; i<=3; i++)
336  for(int j=i; j<=3; j++) {
337 
338  hij_comp.set(i,j).set_etat_qcq() ;
339  hij_comp.set(i,j).import( (comp_hij)(i,j) ) ;
340  }
341 
343 
344 // Relaxation on logn_comp, beta_comp, lnq_comp, hij_comp
345 // ---------------------------------------------------------------
346  double relaxjm1 = 1. - relax ;
347 
348  logn_comp = relax * logn_comp + relaxjm1 * (star_jm1.logn_comp) ;
349 
350  beta_comp = relax * beta_comp + relaxjm1
351  * (star_jm1.beta_comp) ;
352 
353  lnq_comp = relax * lnq_comp + relaxjm1 * (star_jm1.lnq_comp) ;
354 
355 
356  for(int i=1; i<=3; i++)
357  for(int j=i; j<=3; j++) {
358 
359  hij_comp.set(i,j) = relax * hij_comp(i,j)
360  + relaxjm1 * (star_jm1.hij_comp)(i,j) ;
361 
362  }
363 
364 // Lapse function N
365 // ----------------
366 
367  logn = logn_auto + logn_comp ;
368 
369  nn = exp( logn ) ;
370 
371  nn.std_spectral_base() ; // set the bases for spectral expansions
372 
373 
374 // Quantity lnq = log(psi^2 * N)
375 // --------------------------
376 
377  lnq = lnq_auto + lnq_comp ;
378 
379  psi4 = exp(2*lnq) / (nn * nn) ;
381 
382 // Shift vector
383 // ------------
384 
385  beta = beta_auto + beta_comp ;
386 
387 // Coefficients of the 3-metric tilde
388 // ----------------------------------
389 
390  Sym_tensor gtilde_con(mp, CON, mp.get_bvect_cart()) ;
391 
392  for(int i=1; i<=3; i++)
393  for(int j=i; j<=3; j++) {
394 
395  hij.set(i,j) = hij_auto(i,j) + hij_comp(i,j) ;
396  gtilde_con.set(i,j) = hij(i,j) + flat.con()(i,j) ;
397  }
398 
399 
400  gtilde = gtilde_con ;
401  Sym_tensor tens_gamma(gtilde_con / psi4) ;
402  gamma = tens_gamma ;
403 
404  // For conformally flat metrics
405  // ----------------------------
406 
407  if (conf_flat) {
410  hij.set_etat_zero() ;
411  gtilde = flat ;
412  tens_gamma = flat.con() / psi4 ;
413  gamma = tens_gamma ;
414  }
415 
416 // Computation of det(gtilde)
417 
418  Scalar det_gtilde (gtilde.determinant()) ;
419 
420  double* max_det = new double[nz] ;
421  double* min_det = new double[nz] ;
422  double* moy_det = new double[nz] ;
423 
424  for (int i=0; i<nz; i++){
425  min_det[i] = 2 ;
426  moy_det[i] = 0 ;
427  max_det[i] = 0 ;
428  }
429 
430  for (int l=0; l<nz; l++)
431  for (int k=0; k<np; k++)
432  for (int j=0; j<nt; j++)
433  for (int i=0; i<nr; i++){
434 
435  moy_det[l] = moy_det[l] + det_gtilde.val_grid_point(l,k,j,i) ;
436  if (det_gtilde.val_grid_point(l,k,j,i) > max_det[l]){
437  max_det[l] = det_gtilde.val_grid_point(l,k,j,i) ;
438  }
439  if (det_gtilde.val_grid_point(l,k,j,i) < min_det[l]){
440  min_det[l] = det_gtilde.val_grid_point(l,k,j,i) ;
441  }
442  }
443 
444  cout << "average determinant of gtilde in each zone : " << endl ;
445  for (int l=0; l<nz; l++){
446  cout << moy_det[l]/(nr*nt*np) << " " ;
447  }
448  cout << endl ;
449 
450 
451  cout << "maximum of the determinant of gtilde in each zone : " << endl ;
452  for (int l=0; l<nz; l++){
453  cout << max_det[l] << " " ;
454  }
455  cout << endl ;
456 
457  cout << "minimum of the determinant of gtilde in each zone : " << endl ;
458  for (int l=0; l<nz; l++){
459  cout << min_det[l] << " " ;
460  }
461  cout << endl ;
462 
463 
464  // ... extrinsic curvature (tkij_auto and kcar_auto)
465  extrinsic_curvature(om) ;
466 
467 
468 // The derived quantities are obsolete
469 // -----------------------------------
470 
471  del_deriv() ;
472 
473  delete max_det ;
474  delete moy_det ;
475  delete min_det ;
476 
477 }
478 
479 }
Base class for coordinate mappings.
Definition: map.h:670
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
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
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric_flat.C:153
virtual const Scalar & determinant() const
Returns the determinant.
Definition: metric.C:392
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
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: scalar.C:353
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
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition: scalar.h:637
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:324
void import(const Scalar &ci)
Assignment to another Scalar defined on a different mapping.
Definition: scalar_import.C:68
Class for stars in binary system.
Definition: star.h:483
Scalar lnq_auto
Scalar field generated principally by the star.
Definition: star.h:543
Sym_tensor hij_comp
Deviation of the inverse conformal metric from the inverse flat metric generated principally by the ...
Definition: star.h:594
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: star_bin.C:369
Vector beta_comp
Part of the shift vector generated principally by the star (Spherical components with respect to the ...
Definition: star.h:575
void update_metric(const Star_bin &comp, double omega)
Computes metric coefficients from known potentials, when the companion is another star.
Scalar logn_auto
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principally by...
Definition: star.h:527
Scalar logn_comp
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principally by...
Definition: star.h:532
void extrinsic_curvature(double omega)
Computes tkij_auto and akcar_auto from beta_auto, nn and Q.
Scalar lnq_comp
Scalar field generated principally by the companion star.
Definition: star.h:548
Vector beta_auto
Part of the shift vector generated principally by the star (Spherical components with respect to the ...
Definition: star.h:570
Sym_tensor hij
Total deviation of the inverse conformal metric from the inverse flat metric.
Definition: star.h:581
Sym_tensor hij_auto
Deviation of the inverse conformal metric from the inverse flat metric generated principally by the ...
Definition: star.h:588
Scalar psi4
Conformal factor .
Definition: star.h:552
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the star) .
Definition: star.h:562
Metric gtilde
Conformal metric .
Definition: star.h:565
bool conf_flat
true if the 3-metric is conformally flat, false for a more general metric.
Definition: star.h:681
Scalar logn
Logarithm of the lapse N .
Definition: star.h:222
Scalar nn
Lapse function N .
Definition: star.h:225
Metric gamma
3-metric
Definition: star.h:235
const Map & get_mp() const
Returns the mapping.
Definition: star.h:355
Map & mp
Mapping associated with the star.
Definition: star.h:180
Vector beta
Shift vector.
Definition: star.h:228
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:223
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
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition: tensor.C:519
virtual void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition: tensor.C:481
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition: tensor.C:497
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:654
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: tensor.C:926
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tensor.h:866
Lorene prototypes.
Definition: app_hor.h:64