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