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
29char 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
92namespace Lorene {
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
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
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
177
178// Coefficients of the 3-metric tilde
179// ----------------------------------
180
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
191
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
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)
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
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
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
386
387// Coefficients of the 3-metric tilde
388// ----------------------------------
389
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
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
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)
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}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Base class for coordinate mappings.
Definition map.h:670
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
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
virtual const Scalar & determinant() const
Returns the determinant.
Definition metric.C:392
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
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.
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
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
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
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tensor.h:866
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
Lorene prototypes.
Definition app_hor.h:64