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
30char 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
68namespace Lorene {
69double 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
158const 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}
double * p_ham_constr
Relative error on the Hamiltonian constraint.
Definition binaire.h:160
Etoile_bin * et[2]
Array of the two stars (to perform loops on the stars): {\tt et[0]} contains the address of {\tt star...
Definition binaire.h:124
Tbl * p_mom_constr
Relative error on the momentum constraint.
Definition binaire.h:163
Etoile_bin star2
Second star of the system.
Definition binaire.h:118
const Base_vect_cart ref_triad
Cartesian triad of the absolute reference frame.
Definition binaire.h:112
const Tbl & mom_constr() const
Estimates the relative error on the momentum constraint equation by comparing ${\overline\nabla}_j K^...
Etoile_bin star1
First star of the system.
Definition binaire.h:115
double ham_constr() const
Estimates the relative error on the Hamiltonian constraint equation by comparing $\underline\Delta\ln...
const Tenseur & get_d_logn_auto() const
Returns the gradient of logn_auto (Cartesian components with respect to ref_triad )
Definition etoile.h:1121
const Tenseur_sym & get_tkij_auto() const
Returns the part of the extrinsic curvature tensor generated by shift_auto .
Definition etoile.h:1192
const Tenseur & get_d_beta_auto() const
Returns the gradient of beta_auto (Cartesian components with respect to ref_triad )
Definition etoile.h:1141
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
const Tenseur & get_akcar_auto() const
Returns the part of the scalar generated by shift_auto , i.e.
Definition etoile.h:1204
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_d_beta_comp() const
Returns the gradient of beta_comp (Cartesian components with respect to ref_triad )
Definition etoile.h:1146
const Map & get_mp() const
Returns the mapping.
Definition etoile.h:659
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 & get_a_car() const
Returns the total conformal factor .
Definition etoile.h:733
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_u_euler() const
Returns the fluid 3-velocity with respect to the Eulerian observer.
Definition etoile.h:694
const Tenseur & get_press() const
Returns the fluid pressure.
Definition etoile.h:682
const Tenseur & get_ener_euler() const
Returns the total energy density with respect to the Eulerian observer.
Definition etoile.h:685
Basic array class.
Definition tbl.h:161
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
void dec2_dzpuis()
dzpuis -= 2 ;
Definition tenseur.C:1130
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition tenseur.C:1542
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 diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition cmp_math.C:504
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
Cmp abs(const Cmp &)
Absolute value.
Definition cmp_math.C:410
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
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...
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.