LORENE
sym_tensor_tt_etamu.C
1/*
2 * Methods of class Sym_tensor_tt related to eta and mu
3 *
4 * (see file sym_tensor.h for documentation)
5 *
6 */
7
8/*
9 * Copyright (c) 2003-2004 Eric Gourgoulhon & Jerome Novak
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 sym_tensor_tt_etamu_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_tt_etamu.C,v 1.18 2014/10/13 08:53:44 j_novak Exp $" ;
31
32/*
33 * $Id: sym_tensor_tt_etamu.C,v 1.18 2014/10/13 08:53:44 j_novak Exp $
34 * $Log: sym_tensor_tt_etamu.C,v $
35 * Revision 1.18 2014/10/13 08:53:44 j_novak
36 * Lorene classes and functions now belong to the namespace Lorene.
37 *
38 * Revision 1.17 2014/10/06 15:13:19 j_novak
39 * Modified #include directives to use c++ syntax.
40 *
41 * Revision 1.16 2006/10/24 13:03:19 j_novak
42 * New methods for the solution of the tensor wave equation. Perhaps, first
43 * operational version...
44 *
45 * Revision 1.15 2005/04/01 14:28:32 j_novak
46 * Members p_eta and p_mu are now defined in class Sym_tensor.
47 *
48 * Revision 1.14 2004/06/04 09:25:58 e_gourgoulhon
49 * Method eta(): eta is no longer computed from h^rr but from khi (in the
50 * case where khi is known).
51 *
52 * Revision 1.13 2004/05/25 15:08:44 f_limousin
53 * Add parameters in argument of the functions update, eta and mu for
54 * the case of a Map_et.
55 *
56 * Revision 1.12 2004/05/24 13:45:29 e_gourgoulhon
57 * Added parameter dzp to method Sym_tensor_tt::update.
58 *
59 * Revision 1.11 2004/05/05 14:24:54 e_gourgoulhon
60 * Corrected a bug in method set_khi_mu: the division of khi by r^2
61 * was ommitted in the case dzp=2 !!!
62 *
63 * Revision 1.10 2004/04/08 16:38:43 e_gourgoulhon
64 * Sym_tensor_tt::set_khi_mu: added argument dzp (dzpuis of resulting h^{ij}).
65 *
66 * Revision 1.9 2004/03/04 09:53:04 e_gourgoulhon
67 * Methods eta(), mu() and upate(): use of Scalar::mult_r_dzpuis and
68 * change of dzpuis behavior of eta and mu.
69 *
70 * Revision 1.8 2004/03/03 13:16:21 j_novak
71 * New potential khi (p_khi) and the functions manipulating it.
72 *
73 * Revision 1.7 2004/02/05 13:44:50 e_gourgoulhon
74 * Major modif. of methods eta(), mu() and update() to treat
75 * any value of dzpuis, thanks to the new definitions of
76 * Scalar::mult_r(), Scalar::dsdr(), etc...
77 *
78 * Revision 1.6 2004/01/28 13:25:41 j_novak
79 * The ced_mult_r arguments have been suppressed from the Scalar::*dsd* methods.
80 * In the div/mult _r_dzpuis, there is no more default value.
81 *
82 * Revision 1.5 2003/11/05 15:28:31 e_gourgoulhon
83 * Corrected error in update.
84 *
85 * Revision 1.4 2003/11/04 23:03:34 e_gourgoulhon
86 * First full version of method update().
87 * Add method set_rr_mu.
88 * Method set_eta_mu ---> set_rr_eta_mu.
89 *
90 * Revision 1.3 2003/11/04 09:35:27 e_gourgoulhon
91 * First operational version of update_tp().
92 *
93 * Revision 1.2 2003/11/03 22:33:36 e_gourgoulhon
94 * Added methods update_tp and set_eta_mu.
95 *
96 * Revision 1.1 2003/11/03 17:08:37 e_gourgoulhon
97 * First version
98 *
99 *
100 * $Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_tt_etamu.C,v 1.18 2014/10/13 08:53:44 j_novak Exp $
101 *
102 */
103
104// Headers C
105#include <cstdlib>
106#include <cassert>
107
108// Headers Lorene
109#include "tensor.h"
110
111 //--------------//
112 // khi //
113 //--------------//
114
115namespace Lorene {
117
118 if (p_khi == 0x0) { // a new computation is necessary
119
120 // All this has a meaning only for spherical components:
121 assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
122
123 // khi is computed from $h^{rr}$ component
124
125 p_khi = new Scalar(operator()(1,1)) ;
126 p_khi->mult_r() ;
127 p_khi->mult_r() ;
128 }
129
130 return *p_khi ;
131
132}
133
134
135 //--------------//
136 // eta //
137 //--------------//
138
139
141
142
143 if (p_eta == 0x0) { // a new computation is necessary
144
145
146 // All this has a meaning only for spherical components:
147 assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
148
149 // eta is computed from the divergence-free condition:
150
151 int dzp = operator()(1,1).get_dzpuis() ;
152 int dzp_resu = ((dzp == 0) ? 0 : dzp-1) ;
153
154 Scalar source_eta(*mp) ;
155
156 if (p_khi == 0x0) { // eta is computed from h^rr
157 // -------------------------
158
159 source_eta = - operator()(1,1).dsdr() ;
160
161 // dhrr contains - dh^{rr}/dr in all domains but the CED,
162 // in the CED: - r^2 dh^{rr}/dr if dzp = 0 (1)
163 // - r^(dzp+1) dh^{rr}/dr if dzp > 0 (2)
164
165
166 // Multiplication by r of (-d h^{rr}/dr) (with the same dzpuis as h^{rr})
167 source_eta.mult_r_dzpuis( dzp ) ;
168
169 // Substraction of the h^rr part and multiplication by r :
170 source_eta -= 3. * operator()(1,1) ;
171
172 source_eta.mult_r_dzpuis(dzp_resu) ;
173 }
174 else { // eta is computed from khi
175 // ------------------------
176
177 source_eta = - p_khi->dsdr() ;
178 int diff_dzp = source_eta.get_dzpuis() - dzp_resu ;
179 assert( diff_dzp >= 0 ) ;
180 source_eta.dec_dzpuis(diff_dzp) ;
181
182 Scalar tmp(*p_khi) ;
183 tmp.div_r_dzpuis(dzp_resu) ;
184
185 source_eta -= tmp ;
186
187 }
188
189
190 // Resolution of the angular Poisson equation for eta
191 // --------------------------------------------------
192 if (dynamic_cast<const Map_af*>(mp) != 0x0) {
193 p_eta = new Scalar( source_eta.poisson_angu() ) ;
194 }
195 else {
196 Scalar resu (*mp) ;
197 resu = 0. ;
198 mp->poisson_angu(source_eta, *par, resu) ;
199 p_eta = new Scalar( resu ) ;
200 }
201
202 }
203
204 return *p_eta ;
205
206}
207
208
209
210 //-------------------//
211 // set_rr_eta_mu //
212 //-------------------//
213
214
216 const Scalar& mu_i) {
217
218 // All this has a meaning only for spherical components:
219 assert( dynamic_cast<const Base_vect_spher*>(triad) != 0x0 ) ;
220
221 set(1,1) = hrr ; // h^{rr}
222 // calls del_deriv() and therefore delete previous
223 // p_eta and p_mu
224
225 p_eta = new Scalar( eta_i ) ; // eta
226
227 p_mu = new Scalar( mu_i ) ; // mu
228
229 update( hrr.get_dzpuis() ) ; // all h^{ij}, except for h^{rr}
230
231}
232
233 //---------------//
234 // set_rr_mu //
235 //---------------//
236
237
239
240 // All this has a meaning only for spherical components:
241 assert( dynamic_cast<const Base_vect_spher*>(triad) != 0x0 ) ;
242
243 set(1,1) = hrr ; // h^{rr}
244 // calls del_deriv() and therefore delete previous
245 // p_eta and p_mu
246
247 p_mu = new Scalar( mu_i ) ; // mu
248
249 eta() ; // computes eta form the divergence-free condition
250
251 update( hrr.get_dzpuis() ) ; // all h^{ij}, except for h^{rr}
252
253}
254
255 //-------------------//
256 // set_khi_eta_mu //
257 //-------------------//
258
259
261 const Scalar& mu_i) {
262
263 // All this has a meaning only for spherical components:
264 assert( dynamic_cast<const Base_vect_spher*>(triad) != 0x0 ) ;
265
266 set(1,1) = khi_i ;
267 set(1,1).div_r() ;
268 set(1,1).div_r() ; // h^{rr}
269
270 // calls del_deriv() and therefore delete previous
271 // p_khi, p_eta and p_mu
272
273 p_khi = new Scalar( khi_i ) ; // khi
274
275 p_eta = new Scalar( eta_i ) ; // eta
276
277 p_mu = new Scalar( mu_i ) ; // mu
278
279 update( khi_i.get_dzpuis() ) ; // all h^{ij}, except for h^{rr}
280
281}
282
283 //---------------//
284 // set_khi_mu //
285 //---------------//
286
287
289 int dzp, Param* par1, Param* par2, Param* par3) {
290
291 // All this has a meaning only for spherical components:
292 assert( dynamic_cast<const Base_vect_spher*>(triad) != 0x0 ) ;
293
294 set(1,1) = khi_i ;
295 // calls del_deriv() and therefore delete previous
296 // p_eta and p_mu
297
298 assert( khi_i.check_dzpuis(0) ) ;
299 if (dzp == 0) {
300 set(1,1).div_r() ;
301 set(1,1).div_r() ; // h^{rr}
302 }
303 else {
304 assert(dzp == 2) ; //## temporary: the other cases are not treated yet
305 set(1,1).div_r_dzpuis(1) ;
306 set(1,1).div_r_dzpuis(2) ;
307 }
308
309 p_khi = new Scalar ( khi_i ) ; // khi
310
311 p_mu = new Scalar( mu_i ) ; // mu
312
313 if (dynamic_cast<const Map_af*>(mp) != 0x0) {
314 eta() ; // computes eta form the divergence-free condition
315 // dzp = 0 ==> eta.dzpuis = 0
316 // dzp = 2 ==> eta.dzpuis = 1
317
318 update(dzp) ; // all h^{ij}, except for h^{rr}
319 }
320 else {
321 eta(par1) ; // computes eta form the divergence-free condition
322 // dzp = 0 ==> eta.dzpuis = 0
323 // dzp = 2 ==> eta.dzpuis = 1
324
325 update(dzp, par2, par3) ; // all h^{ij}, except for h^{rr}
326 }
327
328}
329
330
331
332 //-------------//
333 // update //
334 //-------------//
335
336
337
339
340 // All this has a meaning only for spherical components:
341 assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
342
343 assert( (p_eta != 0x0) && (p_mu != 0x0) ) ;
344
345 Itbl idx(2) ;
346 idx.set(0) = 1 ; // r index
347
348 // h^{r theta} :
349 // ------------
350 idx.set(1) = 2 ; // theta index
351 *cmp[position(idx)] = p_eta->srdsdt() - p_mu->srstdsdp() ;
352
353 if (dzp == 0) {
354 assert( cmp[position(idx)]->check_dzpuis(2) ) ;
355 cmp[position(idx)]->dec_dzpuis(2) ;
356 }
357
358 assert( cmp[position(idx)]->check_dzpuis(dzp) ) ;
359
360 // h^{r phi} :
361 // ------------
362 idx.set(1) = 3 ; // phi index
363 *cmp[position(idx)] = p_eta->srstdsdp() + p_mu->srdsdt() ;
364
365 if (dzp == 0) {
366 assert( cmp[position(idx)]->check_dzpuis(2) ) ;
367 cmp[position(idx)]->dec_dzpuis(2) ;
368 }
369
370 assert( cmp[position(idx)]->check_dzpuis(dzp) ) ;
371
372 // h^{theta phi} and h^{phi phi}
373 // -----------------------------
374
375 //-------------- Computation of T^theta --> taut :
376
377 Scalar tautst = operator()(1,2).dsdr() ;
378
379 // dhrr contains dh^{rt}/dr in all domains but the CED,
380 // in the CED: r^2 dh^{rt}/dr if dzp = 0 (1)
381 // r^(dzp+1) dh^{rt}/dr if dzp > 0 (2)
382
383 // Multiplication by r of dh^{rt}/dr (with the same dzpuis than h^{rt})
384 tautst.mult_r_dzpuis( operator()(1,2).get_dzpuis() ) ;
385
386
387 // Addition of the remaining parts :
388 tautst += 3 * operator()(1,2) - operator()(1,1).dsdt() ;
389 tautst.mult_sint() ;
390
391 Scalar tmp = operator()(1,1) ;
392 tmp.mult_cost() ; // h^{rr} cos(th)
393
394 tautst -= tmp ; // T^th / sin(th)
395
396 Scalar taut = tautst ;
397 taut.mult_sint() ; // T^th
398
399
400 //----------- Computation of T^phi --> taup :
401
402 Scalar taupst = - operator()(1,3).dsdr() ;
403
404 // dhrr contains - dh^{rp}/dr in all domains but the CED,
405 // in the CED: - r^2 dh^{rp}/dr if dzp = 0 (3)
406 // - r^(dzp+1) dh^{rp}/dr if dzp > 0 (4)
407
408 // Multiplication by r of -dh^{rp}/dr (with the same dzpuis than h^{rp})
409 taupst.mult_r_dzpuis( operator()(1,3).get_dzpuis() ) ;
410
411
412 // Addition of the remaining part :
413
414 taupst -= 3 * operator()(1,3) ;
415 taupst.mult_sint() ; // T^ph / sin(th)
416
417 Scalar taup = taupst ;
418 taup.mult_sint() ; // T^ph
419
420
421 //------------------- Computation of F and h^[ph ph}
422
423 tmp = tautst ;
424 tmp.mult_cost() ; // T^th / tan(th)
425
426 // dT^th/dth + T^th / tan(th) + 1/sin(th) dT^ph/dph :
427 tmp = taut.dsdt() + tmp + taup.stdsdp() ;
428
429 Scalar tmp2 (*mp) ;
430 if (dynamic_cast<const Map_af*>(mp) != 0x0) {
431 tmp2 = tmp.poisson_angu() ; // F
432 }
433 else {
434 tmp2 = 0. ;
435 mp->poisson_angu(tmp, *par1, tmp2) ; // F
436 }
437
438
439 tmp2.div_sint() ;
440 tmp2.div_sint() ; // h^{ph ph}
441
442 idx.set(0) = 3 ; // phi index
443 idx.set(1) = 3 ; // phi index
444 *cmp[position(idx)] = tmp2 ; // h^{ph ph} is updated
445
446
447 //------------------- Computation of G and h^[th ph}
448
449 tmp = taupst ;
450 tmp.mult_cost() ; // T^ph / tan(th)
451
452 // - 1/sin(th) dT^th/dph + dT^ph/dth + T^ph / tan(th) :
453 tmp = - taut.stdsdp() + taup.dsdt() + tmp ;
454
455 if (dynamic_cast<const Map_af*>(mp) != 0x0) {
456 tmp2 = tmp.poisson_angu() ; // G
457 }
458 else {
459 tmp2 = 0. ;
460 mp->poisson_angu(tmp, *par2, tmp2) ; // G
461 }
462
463 tmp2.div_sint() ;
464 tmp2.div_sint() ; // h^{th ph}
465
466 idx.set(0) = 2 ; // theta index
467 idx.set(1) = 3 ; // phi index
468 *cmp[position(idx)] = tmp2 ; // h^{th ph} is updated
469
470 // h^{th th} (from the trace-free condition)
471 // ---------
472 idx.set(1) = 2 ; // theta index
473 *cmp[position(idx)] = - operator()(1,1) - operator()(3,3) ;
474
475
476 Sym_tensor_trans::del_deriv() ; //## in order not to delete p_eta and p_mu
477
478
479
480}
481
482
483 //-----------------//
484 // set_A_tilde_B //
485 //-----------------//
486
489
490 Scalar zero(*mp) ;
491 zero.set_etat_zero() ;
493 return ;
494}
495
496
497
498
499
500
501}
Spherical orthonormal vectorial bases (triads).
Definition base_vect.h:308
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Basic integer array class.
Definition itbl.h:122
Affine radial mapping.
Definition map.h:2027
Parameter storage.
Definition param.h:125
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
const Scalar & srdsdt() const
Returns of *this .
void div_r_dzpuis(int ced_mult_r)
Division by r everywhere but with the output flag dzpuis set to ced_mult_r .
const Scalar & srstdsdp() const
Returns of *this .
void div_r()
Division by r everywhere; dzpuis is not changed.
const Scalar & dsdr() const
Returns of *this .
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values of the Scalar in the co...
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
void set_AtB_trace(const Scalar &a_in, const Scalar &tb_in, const Scalar &trace, Param *par_bc=0x0, Param *par_mat=0x0)
Assigns the derived members A , and the trace.
virtual void del_deriv() const
Deletes the derived quantities.
void set_khi_mu(const Scalar &khi_i, const Scalar &mu_i, int dzp=0, Param *par1=0x0, Param *par2=0x0, Param *par3=0x0)
Sets the component , as well as the angular potential (see member p_khi and p_mu ).
void set_rr_mu(const Scalar &hrr, const Scalar &mu_i)
Sets the component , as well as the angular potential (see member p_mu ).
void set_A_tildeB(const Scalar &a_in, const Scalar &tb_in, Param *par_bc=0x0, Param *par_mat=0x0)
Assigns the derived members A and .
Scalar * p_khi
Field such that the component .
Definition sym_tensor.h:946
void set_khi_eta_mu(const Scalar &khi_i, const Scalar &eta_i, const Scalar &mu_i)
Sets the component , as well as the angular potentials and (see members p_khi , p_eta and p_mu ).
void update(int dzp, Param *par1=0x0, Param *par2=0x0)
Computes the components , , , and , from and the potentials and .
virtual const Scalar & eta(Param *par=0x0) const
Gives the field (see member p_eta ).
void set_rr_eta_mu(const Scalar &hrr, const Scalar &eta_i, const Scalar &mu_i)
Sets the component , as well as the angular potentials and (see members p_eta and p_mu ).
const Scalar & khi() const
Gives the field such that the component .
Scalar * p_mu
Field such that the components of the tensor are written (has only meaning with spherical component...
Definition sym_tensor.h:274
Scalar * p_eta
Field such that the components of the tensor are written (has only meaning with spherical component...
Definition sym_tensor.h:260
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
const Scalar & operator()(const Itbl &ind) const
Returns the value of a component (read-only version).
Definition tensor.C:798
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition tensor.h:315
virtual int position(const Itbl &ind) const
Returns the position in the array cmp of a component given by its indices.
Definition tensor_sym.C:245
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:654
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition tensor.h:303
Lorene prototypes.
Definition app_hor.h:64