LORENE
binaire_constr.C
1 /*
2  * Methods of class Binaire for estimating the error in the Hamiltionian
3  * and momentum constraints
4  *
5  * (see file binaire.h for documentation).
6  */
7 
8 /*
9  * Copyright (c) 2000-2001 Eric Gourgoulhon
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 char binaire_constr_C[] = "$Header: /cvsroot/Lorene/C++/Source/Binaire/binaire_constr.C,v 1.3 2014/10/13 08:52:44 j_novak Exp $" ;
31 
32 /*
33  * $Id: binaire_constr.C,v 1.3 2014/10/13 08:52:44 j_novak Exp $
34  * $Log: binaire_constr.C,v $
35  * Revision 1.3 2014/10/13 08:52:44 j_novak
36  * Lorene classes and functions now belong to the namespace Lorene.
37  *
38  * Revision 1.2 2004/03/25 10:28:59 j_novak
39  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
40  *
41  * Revision 1.1.1.1 2001/11/20 15:19:30 e_gourgoulhon
42  * LORENE
43  *
44  * Revision 2.1 2000/03/13 17:05:34 eric
45  * *** empty log message ***
46  *
47  * Revision 2.0 2000/03/13 14:26:08 eric
48  * *** empty log message ***
49  *
50  *
51  * $Header: /cvsroot/Lorene/C++/Source/Binaire/binaire_constr.C,v 1.3 2014/10/13 08:52:44 j_novak Exp $
52  *
53  */
54 
55 // Headers C
56 #include "math.h"
57 
58 // Headers Lorene
59 #include "binaire.h"
60 #include "unites.h"
61 
62 
63 
64  //----------------------------------------------//
65  // Hamiltonian constraint //
66  //----------------------------------------------//
67 
68 namespace Lorene {
69 double Binaire::ham_constr() const {
70 
71  using namespace Unites ;
72 
73  if (p_ham_constr == 0x0) { // A new computation is required
74 
75 
76  Tenseur lap_alpha1( star1.get_mp() ) ;
77  Tenseur lap_alpha2( star2.get_mp() ) ;
78 
79  Tenseur source1( star1.get_mp() ) ;
80  Tenseur source2( star2.get_mp() ) ;
81 
82  Tenseur* p_lap_alpha[2] ;
83  Tenseur* p_source[2] ;
84  p_lap_alpha[0] = &lap_alpha1 ;
85  p_lap_alpha[1] = &lap_alpha2 ;
86  p_source[0] = &source1 ;
87  p_source[1] = &source2 ;
88 
89 
90  // Computation of the l.h.s. and r.h.s. of the Hamiltonian
91  // constraint in each star.
92  // -------------------------------------------------------
93 
94  double som = 0 ;
95 
96  for (int i=0; i<2; i++) {
97 
98  // Laplacian of alpha = ln(A)
99  // --------------------------
100 
101  Tenseur alpha_auto = et[i]->get_beta_auto()
102  - et[i]->get_logn_auto() ;
103 
104  *(p_lap_alpha[i]) = alpha_auto().laplacien() ;
105 
106  // Right-hand-side of the Hamiltonian constraint
107  // ---------------------------------------------
108 
109  const Tenseur& a_car = et[i]->get_a_car() ;
110  const Tenseur& ener_euler = et[i]->get_ener_euler() ;
111 
112  Tenseur d_alpha_auto = et[i]->get_d_beta_auto()
113  - et[i]->get_d_logn_auto() ;
114 
115  Tenseur d_alpha_comp = et[i]->get_d_beta_comp()
116  - et[i]->get_d_logn_comp() ;
117 
118  const Tenseur& akcar_auto = et[i]->get_akcar_auto() ;
119  const Tenseur& akcar_comp = et[i]->get_akcar_comp() ;
120 
121  *(p_source[i]) = - qpig * a_car * ener_euler
122  - 0.25 * ( akcar_auto + akcar_comp )
123  - 0.5 *
124  ( flat_scalar_prod(d_alpha_auto, d_alpha_auto)
125  + flat_scalar_prod(d_alpha_auto, d_alpha_comp)
126  ) ;
127 
128  // Relative difference
129  // -------------------
130  Tbl diff = diffrel( (*(p_lap_alpha[i]))(), (*(p_source[i]))() ) ;
131 
132  cout <<
133  "Binaire::ham_constr : relative difference Lap(alpha) <-> source : "
134  << endl << diff << endl ;
135 
136  som += max( abs(diff) ) ;
137 
138  }
139 
140 
141  // Total error
142  // -----------
143  p_ham_constr = new double ;
144 
145  *p_ham_constr = 0.5 * som ;
146 
147  }
148 
149  return *p_ham_constr ;
150 
151 }
152 
153 
154  //----------------------------------------------//
155  // Momentum constraint //
156  //----------------------------------------------//
157 
158 const Tbl& Binaire::mom_constr() const {
159 
160  using namespace Unites ;
161 
162  if (p_mom_constr == 0x0) { // A new computation is required
163 
164  Tenseur divk1( star1.get_mp(), 1, CON, ref_triad ) ;
165  Tenseur divk2( star2.get_mp(), 1, CON, ref_triad ) ;
166 
167  Tenseur source1( star1.get_mp(), 1, CON, ref_triad ) ;
168  Tenseur source2( star2.get_mp(), 1, CON, ref_triad ) ;
169 
170  Tenseur* p_divk[2] ;
171  Tenseur* p_source[2] ;
172  p_divk[0] = &divk1 ;
173  p_divk[1] = &divk2 ;
174  p_source[0] = &source1 ;
175  p_source[1] = &source2 ;
176 
177 
178  // Computation of the l.h.s. and r.h.s. of the momentum
179  // constraint in each star.
180  // -------------------------------------------------------
181 
182  double somx = 0 ;
183  double somy = 0 ;
184  double somz = 0 ;
185 
186  for (int i=0; i<2; i++) {
187 
188  // (flat space) divergence of K^{ij}
189  // ---------------------------------
190 
191  const Tenseur& a_car = et[i]->get_a_car() ;
192  Tenseur kij_auto = et[i]->get_tkij_auto() / a_car ;
193 
194  kij_auto.dec2_dzpuis() ; // dzpuis : 2 --> 0
195  // so that in the external domain, kij_auto
196  // contains now exactly K^{ij}
197 
198  // The gradient of K^{ij} is computed on the local triad:
199  kij_auto.change_triad( (et[i]->get_mp()).get_bvect_cart() ) ;
200 
201  *(p_divk[i]) = contract( kij_auto.gradient(), 0, 1) ;
202 
203  // Back to the Reference triad :
204  p_divk[i]->change_triad( ref_triad ) ;
205  kij_auto.change_triad( ref_triad ) ;
206 
207  // Right-hand-side of the momentum constraint
208  // ------------------------------------------
209 
210  const Tenseur& u_euler = et[i]->get_u_euler() ;
211  const Tenseur& ener_euler = et[i]->get_ener_euler() ;
212  const Tenseur& press = et[i]->get_press() ;
213 
214 
215  Tenseur d_alpha = et[i]->get_d_beta_auto()
216  - et[i]->get_d_logn_auto()
217  + et[i]->get_d_beta_comp()
218  - et[i]->get_d_logn_comp() ;
219 
220  *(p_source[i]) = 2 * qpig * (ener_euler + press) * u_euler
221  - 5 * contract(kij_auto, 1, d_alpha, 0) ;
222 
223  // Relative differences
224  // --------------------
225  Tbl diffx = diffrel( (*(p_divk[i]))(0), (*(p_source[i]))(0)) ;
226  Tbl diffy = diffrel( (*(p_divk[i]))(1), (*(p_source[i]))(1)) ;
227  Tbl diffz = diffrel( (*(p_divk[i]))(2), (*(p_source[i]))(2)) ;
228 
229  cout << "Binaire::mom_constr : norme div(K) : " << endl ;
230  cout << "X component : " << norme( (*(p_divk[i]))(0) ) << endl ;
231  cout << "Y component : " << norme( (*(p_divk[i]))(1) ) << endl ;
232  cout << "Z component : " << norme( (*(p_divk[i]))(2) ) << endl ;
233 
234  cout << "Binaire::mom_constr : norme source : " << endl ;
235  cout << "X component : " << norme( (*(p_source[i]))(0) ) << endl ;
236  cout << "Y component : " << norme( (*(p_source[i]))(1) ) << endl ;
237  cout << "Z component : " << norme( (*(p_source[i]))(2) ) << endl ;
238 
239 
240  cout <<
241  "Binaire::mom_constr : rel. diff. div(K) <-> source : "
242  << endl ;
243  cout << "X component : " << diffx << endl ;
244  cout << "Y component : " << diffy << endl ;
245  cout << "Z component : " << diffz << endl ;
246 
247 
248  somx += max( abs(diffx) ) ;
249  somy += max( abs(diffy) ) ;
250  somz += max( abs(diffz) ) ;
251  }
252 
253  // Total error
254  // -----------
255 
256  p_mom_constr = new Tbl(3) ;
258 
259  p_mom_constr->set(0) = 0.5 * somx ;
260  p_mom_constr->set(1) = 0.5 * somy ;
261  p_mom_constr->set(2) = 0.5 * somz ;
262 
263 
264  }
265 
266  return *p_mom_constr ;
267 
268 }
269 }
const Map & get_mp() const
Returns the mapping.
Definition: etoile.h:659
void dec2_dzpuis()
dzpuis -= 2 ;
Definition: tenseur.C:1130
const Tenseur & get_a_car() const
Returns the total conformal factor .
Definition: etoile.h:733
Etoile_bin star1
First star of the system.
Definition: binaire.h:115
Lorene prototypes.
Definition: app_hor.h:64
const Tenseur & get_akcar_comp() const
Returns the part of the scalar generated by shift_auto and shift_comp , i.e.
Definition: etoile.h:1210
Tbl * p_mom_constr
Relative error on the momentum constraint.
Definition: binaire.h:163
Standard units of space, time and mass.
const Tenseur_sym & get_tkij_auto() const
Returns the part of the extrinsic curvature tensor generated by shift_auto .
Definition: etoile.h:1192
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
const Tenseur & get_d_beta_auto() const
Returns the gradient of beta_auto (Cartesian components with respect to ref_triad ) ...
Definition: etoile.h:1141
double * p_ham_constr
Relative error on the Hamiltonian constraint.
Definition: binaire.h:160
Tenseur flat_scalar_prod(const Tenseur &t1, const Tenseur &t2)
Scalar product of two Tenseur when the metric is : performs the contraction of the last index of t1 w...
const Base_vect_cart ref_triad
Cartesian triad of the absolute reference frame.
Definition: binaire.h:112
const Tenseur & get_logn_auto() const
Returns the logarithm of the part of the lapse N generated principaly by the star.
Definition: etoile.h:701
const Tenseur & get_d_logn_auto() const
Returns the gradient of logn_auto (Cartesian components with respect to ref_triad ) ...
Definition: etoile.h:1121
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
const Tenseur & get_press() const
Returns the fluid pressure.
Definition: etoile.h:682
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
const Tbl & mom_constr() const
Estimates the relative error on the momentum constraint equation by comparing ${}_j K^{ij}$ with {equ...
Etoile_bin * et[2]
Array of the two stars (to perform loops on the stars): { et[0]} contains the address of { star1} and...
Definition: binaire.h:124
double ham_constr() const
Estimates the relative error on the Hamiltonian constraint equation by comparing $ A$ with {equation}...
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition: tenseur.C:668
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:481
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:435
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
const Tenseur & get_d_beta_comp() const
Returns the gradient of beta_comp (Cartesian components with respect to ref_triad ) ...
Definition: etoile.h:1146
const Tenseur & get_ener_euler() const
Returns the total energy density with respect to the Eulerian observer.
Definition: etoile.h:685
const Tenseur & get_d_logn_comp() const
Returns the gradient of logn_comp (Cartesian components with respect to ref_triad ) ...
Definition: etoile.h:1131
const Tenseur & get_akcar_auto() const
Returns the part of the scalar generated by shift_auto , i.e.
Definition: etoile.h:1204
Cmp abs(const Cmp &)
Absolute value.
Definition: cmp_math.C:410
Basic array class.
Definition: tbl.h:161
const Tenseur & get_u_euler() const
Returns the fluid 3-velocity with respect to the Eulerian observer.
Definition: etoile.h:694
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:298
const Tenseur & get_beta_auto() const
Returns the logarithm of the part of the product AN generated principaly by the star.
Definition: etoile.h:724
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition: tenseur.C:1542
Etoile_bin star2
Second star of the system.
Definition: binaire.h:118