LORENE
sym_tensor_trans_aux.C
1 /*
2  * Manipulation of auxiliary potentials for Sym_tensor_trans objects.
3  *
4  * (see file sym_tensor.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2005-2006 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 version 2
15  * as published by the Free Software Foundation.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 char sym_tensor_trans_aux_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_trans_aux.C,v 1.21 2014/10/13 08:53:43 j_novak Exp $" ;
29 
30 /*
31  * $Id: sym_tensor_trans_aux.C,v 1.21 2014/10/13 08:53:43 j_novak Exp $
32  * $Log: sym_tensor_trans_aux.C,v $
33  * Revision 1.21 2014/10/13 08:53:43 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.20 2014/10/06 15:13:19 j_novak
37  * Modified #include directives to use c++ syntax.
38  *
39  * Revision 1.19 2010/10/11 10:23:03 j_novak
40  * Removed methods Sym_tensor_trans::solve_hrr() and Sym_tensor_trans::set_WX_det_one(), as they are no longer relevant.
41  *
42  * Revision 1.18 2008/12/05 08:46:19 j_novak
43  * New method Sym_tensor_trans_aux::set_tt_part_det_one.
44  *
45  * Revision 1.17 2007/04/27 13:52:55 j_novak
46  * In method Sym_tensor_trans::set_AtBtt_det_one(), the member p_tt (the TT-part)
47  * is updated.
48  *
49  * Revision 1.16 2007/03/20 12:20:56 j_novak
50  * In Sym_tensor_trans::set_AtBtt_det_one(), the trace is stored in the resulting
51  * object.
52  *
53  * Revision 1.15 2006/10/24 13:03:19 j_novak
54  * New methods for the solution of the tensor wave equation. Perhaps, first
55  * operational version...
56  *
57  * Revision 1.14 2006/09/15 08:48:01 j_novak
58  * Suppression of some messages in the optimized version.
59  *
60  * Revision 1.13 2006/08/31 12:13:22 j_novak
61  * Added an argument of type Param to Sym_tensor_trans::sol_Dirac_A().
62  *
63  * Revision 1.12 2006/06/26 09:28:13 j_novak
64  * Added a forgotten initialisation in set_AtB_trace_zero().
65  *
66  * Revision 1.11 2006/06/20 12:07:15 j_novak
67  * Improved execution speed for sol_Dirac_tildeB...
68  *
69  * Revision 1.10 2006/06/14 10:04:21 j_novak
70  * New methods sol_Dirac_l01, set_AtB_det_one and set_AtB_trace_zero.
71  *
72  * Revision 1.9 2006/01/17 15:50:53 j_novak
73  * Slight re-arrangment of the det=1 formula.
74  *
75  * Revision 1.8 2005/11/28 14:45:17 j_novak
76  * Improved solution of the Poisson tensor equation in the case of a transverse
77  * tensor.
78  *
79  * Revision 1.7 2005/09/16 13:34:44 j_novak
80  * Back to dzpuis 1 for the source for mu. eta is computed the same way as hrr.
81  *
82  * Revision 1.6 2005/09/08 16:00:23 j_novak
83  * Change of dzpuis for source for mu (temporary?).
84  *
85  * Revision 1.5 2005/09/07 16:47:43 j_novak
86  * Removed method Sym_tensor_trans::T_from_det_one
87  * Modified Sym_tensor::set_auxiliary, so that it takes eta/r and mu/r as
88  * arguments.
89  * Modified Sym_tensor_trans::set_hrr_mu.
90  * Added new protected method Sym_tensor_trans::solve_hrr
91  *
92  * Revision 1.4 2005/04/08 08:22:04 j_novak
93  * New methods set_hrr_mu_det_one() and set_WX_det_one(). Not tested yet...
94  *
95  * Revision 1.3 2005/04/07 07:56:22 j_novak
96  * Better handling of dzpuis (first try).
97  *
98  * Revision 1.2 2005/04/06 15:49:46 j_novak
99  * Error corrected.
100  *
101  * Revision 1.1 2005/04/06 15:43:59 j_novak
102  * New method Sym_tensor_trans::T_from_det_one(...).
103  *
104  *
105  * $Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_trans_aux.C,v 1.21 2014/10/13 08:53:43 j_novak Exp $
106  *
107  */
108 
109 // C headers
110 #include <cassert>
111 
112 // Lorene headers
113 #include "metric.h"
114 #include "param.h"
115 
116 namespace Lorene {
117 void Sym_tensor_trans::set_hrr_mu_det_one(const Scalar& hrr, const Scalar& mu_in,
118  double precis, int it_max ) {
119 
120  // All this has a meaning only for spherical components:
121  assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
122  assert(hrr.check_dzpuis(0)) ;
123  assert(mu_in.check_dzpuis(0)) ;
124  assert(&mu_in != p_mu) ;
125  assert( (precis > 0.) && (it_max > 0) ) ;
126 
127  Sym_tensor_tt tens_tt(*mp, *triad, *met_div) ;
128  tens_tt.set_rr_mu(hrr, mu_in) ;
129  tens_tt.inc_dzpuis(2) ;
130  trace_from_det_one(tens_tt, precis, it_max) ;
131  dec_dzpuis(2) ;
132 
133  return ;
134 
135 }
136 
137 void Sym_tensor_trans::set_AtBtt_det_one(const Scalar& a_in, const Scalar& tbtt_in,
138  const Scalar* h_prev, Param* par_bc,
139  Param* par_mat, double precis,
140  int it_max ) {
141  // All this has a meaning only for spherical components:
142  assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
143  assert(a_in.check_dzpuis(2)) ;
144  assert(tbtt_in.check_dzpuis(2)) ;
145  assert(&a_in != p_aaa) ;
146  assert(&tbtt_in != p_tilde_b) ;
147  assert( (precis > 0.) && (it_max > 0) ) ;
148 
149  //Computation of mu and X from A
150  //-------------------------------
151  Scalar mu_over_r(*mp) ;
152  Scalar x_new(*mp) ;
153  sol_Dirac_A(a_in, mu_over_r, x_new, par_bc) ;
154 
155  // Preparation for the iteration
156  //------------------------------
157  Scalar h_old(*mp) ;
158  if (h_prev != 0x0)
159  h_old = *h_prev ;
160  else
161  h_old.set_etat_zero() ;
162  double lambda = 1. ;
163  Scalar zero(*mp) ;
164  zero.set_etat_zero() ;
165 
166  Scalar hrr_tt(*mp) ;
167  Scalar eta_sr_tt(*mp) ;
168  Scalar w_tt(*mp) ;
169  sol_Dirac_tilde_B(tbtt_in, zero, hrr_tt, eta_sr_tt, w_tt, par_bc, par_mat) ;
170  Sym_tensor_tt hijtt(*mp, *triad, *met_div) ;
171  hijtt.set_auxiliary(hrr_tt, eta_sr_tt, mu_over_r, w_tt, x_new, zero) ;
172 
173  Scalar hrr_new(*mp) ;
174  Scalar eta_over_r(*mp) ;
175  Scalar w_new(*mp) ;
176 
177  for (int it=0; it<=it_max; it++) {
178  Scalar tilde_B = get_tilde_B_from_TT_trace(zero, h_old) ;
179  sol_Dirac_tilde_B(tilde_B, h_old, hrr_new, eta_over_r, w_new, 0x0, par_mat) ;
180 
181  set_auxiliary(hrr_new+hrr_tt, eta_over_r+eta_sr_tt, mu_over_r,
182  w_new+w_tt, x_new, h_old - hrr_new-hrr_tt) ;
183 
184  const Sym_tensor_trans& hij = *this ;
185  Scalar h_new = (1 + hij(1,1))*( hij(2,3)*hij(2,3) - hij(2,2)*hij(3,3) )
186  + hij(1,2)*hij(1,2)*(1 + hij(3,3))
187  + hij(1,3)*hij(1,3)*(1 + hij(2,2))
188  - hij(1,1)*(hij(2,2) + hij(3,3)) - 2*hij(1,2)*hij(1,3)*hij(2,3) ;
189  h_new.set_spectral_base(hrr_tt.get_spectral_base()) ;
190 
191  double diff = max(max(abs(h_new - h_old))) ;
192 #ifndef NDEBUG
193  cout << "Sym_tensor_trans::set_AtB_det_one : "
194  << "iteration : " << it << " convergence on h: "
195  << diff << endl ;
196 #endif
197  if (diff < precis) {
198  set_auxiliary(hrr_new+hrr_tt, eta_over_r+eta_sr_tt, mu_over_r,
199  w_new+w_tt, x_new, h_old - hrr_new-hrr_tt) ;
200  if (p_aaa != 0x0) delete p_aaa ;
201  p_aaa = new Scalar(a_in) ;
202  if (p_tilde_b != 0x0) delete p_tilde_b ;
203  p_tilde_b = new Scalar(tilde_B + tbtt_in) ;
204  if (p_trace != 0x0) delete p_trace ;
205  p_trace = new Scalar(h_old) ;
206  if (p_tt != 0x0) delete p_tt ;
207  p_tt = new Sym_tensor_tt(hijtt) ;
208  p_tt->inc_dzpuis(2) ;
209  break ;
210  }
211  else {
212  h_old = lambda*h_new +(1-lambda)*h_old ;
213  }
214 
215  if (it == it_max) {
216  cout << "Sym_tensor_trans:::set_AtBtt_det_one : convergence not reached \n" ;
217  cout << " for the required accuracy (" << precis << ") ! "
218  << endl ;
219  abort() ;
220  }
221  }
222  return ;
223 
224 }
225 
227  Scalar* h_prev, Param* par_mat,
228  double precis, int it_max ) {
229  // All this has a meaning only for spherical components:
230  assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
231  assert( (precis > 0.) && (it_max > 0) ) ;
232 
233  Scalar mu_over_r = hijtt.mu() ;
234  mu_over_r.div_r() ;
235  const Scalar& x_new = hijtt.xxx() ;
236 
237  // Preparation for the iteration
238  //------------------------------
239  Scalar h_old(*mp) ;
240  if (h_prev != 0x0)
241  h_old = *h_prev ;
242  else
243  h_old.set_etat_zero() ;
244  double lambda = 1. ;
245  Scalar zero(*mp) ;
246  zero.set_etat_zero() ;
247 
248  const Scalar& hrr_tt = hijtt( 1, 1 ) ;
249  Scalar eta_sr_tt = hijtt.eta() ;
250  eta_sr_tt.div_r() ;
251  const Scalar w_tt = hijtt.www() ;
252 
253  Scalar hrr_new(*mp) ;
254  Scalar eta_over_r(*mp) ;
255  Scalar w_new(*mp) ;
256 
257  for (int it=0; it<=it_max; it++) {
258  Scalar tilde_B = get_tilde_B_from_TT_trace(zero, h_old) ;
259  sol_Dirac_tilde_B(tilde_B, h_old, hrr_new, eta_over_r, w_new, 0x0, par_mat) ;
260 
261  set_auxiliary(hrr_new+hrr_tt, eta_over_r+eta_sr_tt, mu_over_r,
262  w_new+w_tt, x_new, h_old - hrr_new-hrr_tt) ;
263 
264  const Sym_tensor_trans& hij = *this ;
265  Scalar h_new = (1 + hij(1,1))*( hij(2,3)*hij(2,3) - hij(2,2)*hij(3,3) )
266  + hij(1,2)*hij(1,2)*(1 + hij(3,3))
267  + hij(1,3)*hij(1,3)*(1 + hij(2,2))
268  - hij(1,1)*(hij(2,2) + hij(3,3)) - 2*hij(1,2)*hij(1,3)*hij(2,3) ;
269  h_new.set_spectral_base(hrr_tt.get_spectral_base()) ;
270 
271  double diff = max(max(abs(h_new - h_old))) ;
272 #ifndef NDEBUG
273  cout << "Sym_tensor_trans::set_tt_part_det_one : "
274  << "iteration : " << it << " convergence on h: "
275  << diff << endl ;
276 #endif
277  if (diff < precis) {
278  set_auxiliary(hrr_new+hrr_tt, eta_over_r+eta_sr_tt, mu_over_r,
279  w_new+w_tt, x_new, h_old - hrr_new-hrr_tt) ;
280  if (p_trace != 0x0) delete p_trace ;
281  p_trace = new Scalar(h_old) ;
282  if (p_tt != 0x0) delete p_tt ;
283  p_tt = new Sym_tensor_tt(hijtt) ;
284  p_tt->inc_dzpuis(2) ;
285  break ;
286  }
287  else {
288  h_old = lambda*h_new +(1-lambda)*h_old ;
289  }
290 
291  if (it == it_max) {
292  cout << "Sym_tensor_trans:::set_AtBtt_det_one : convergence not reached \n" ;
293  cout << " for the required accuracy (" << precis << ") ! "
294  << endl ;
295  abort() ;
296  }
297  }
298  return ;
299 
300 }
301 
302 void Sym_tensor_trans::set_AtB_trace(const Scalar& a_in, const Scalar& tb_in,
303  const Scalar& hh, Param* par_bc, Param* par_mat) {
304 
305  // All this has a meaning only for spherical components:
306  assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
307  assert(a_in.check_dzpuis(2)) ;
308  assert(tb_in.check_dzpuis(2)) ;
309  assert(hh.check_dzpuis(0)) ;
310  assert(&a_in != p_aaa) ;
311  assert(&tb_in != p_tilde_b) ;
312 
313  //Computation of mu and X from A
314  //-------------------------------
315  Scalar mu_over_r(*mp) ;
316  Scalar x_new(*mp) ;
317  sol_Dirac_A(a_in, mu_over_r, x_new, par_bc) ;
318 
319  // Computation of the other potentials
320  //------------------------------------
321  Scalar hrr_new(*mp) ;
322  Scalar eta_over_r(*mp) ;
323  Scalar w_new(*mp) ;
324 
325  sol_Dirac_tilde_B(tb_in, hh, hrr_new, eta_over_r, w_new, par_bc, par_mat) ;
326 
327  set_auxiliary(hrr_new, eta_over_r, mu_over_r, w_new, x_new, hh - hrr_new) ;
328  if (p_aaa != 0x0) delete p_aaa ;
329  p_aaa = new Scalar(a_in) ;
330  if (p_tilde_b != 0x0) delete p_tilde_b ;
331  p_tilde_b = new Scalar(tb_in) ;
332 
333  return ;
334 
335 }
336 }
Lorene::Sym_tensor_trans::set_hrr_mu_det_one
void set_hrr_mu_det_one(const Scalar &hrr, const Scalar &mu_in, double precis=1.e-14, int it_max=100)
Assigns the rr component and the derived member .
Definition: sym_tensor_trans_aux.C:117
Lorene::Sym_tensor_trans::sol_Dirac_A
void sol_Dirac_A(const Scalar &aaa, Scalar &tilde_mu, Scalar &xxx, const Param *par_bc=0x0) const
Solves a system of two coupled first-order PDEs obtained from the divergence-free condition (Dirac ga...
Definition: sym_tensor_trans_dirac.C:82
Lorene::Scalar::set_spectral_base
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
Definition: scalar.C:797
Lorene::Sym_tensor::get_tilde_B_from_TT_trace
Scalar get_tilde_B_from_TT_trace(const Scalar &tilde_B_tt_in, const Scalar &trace) const
Computes (see Sym_tensor::p_tilde_b ) from its transverse-traceless part and the trace.
Definition: sym_tensor_aux.C:531
Lorene::Sym_tensor_tt
Transverse and traceless symmetric tensors of rank 2.
Definition: sym_tensor.h:938
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Sym_tensor::www
const Scalar & www() const
Gives the field W (see member p_www ).
Definition: sym_tensor_aux.C:209
Lorene::Tensor::dec_dzpuis
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition: tensor.C:808
Lorene::Sym_tensor_trans::sol_Dirac_tilde_B
void sol_Dirac_tilde_B(const Scalar &tilde_b, const Scalar &hh, Scalar &hrr, Scalar &tilde_eta, Scalar &www, Param *par_bc=0x0, Param *par_mat=0x0) const
Solves a system of three coupled first-order PDEs obtained from divergence-free conditions (Dirac gau...
Definition: sym_tensor_trans_dirac.C:583
Lorene::Sym_tensor::set_auxiliary
void set_auxiliary(const Scalar &trr, const Scalar &eta_over_r, const Scalar &mu_over_r, const Scalar &www, const Scalar &xxx, const Scalar &ttt)
Assigns the component and the derived members p_eta , p_mu , p_www, p_xxx and p_ttt ,...
Definition: sym_tensor_aux.C:266
Lorene::Sym_tensor_tt::eta
virtual const Scalar & eta(Param *par=0x0) const
Gives the field (see member p_eta ).
Definition: sym_tensor_tt_etamu.C:140
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:435
Lorene::Sym_tensor::p_mu
Scalar * p_mu
Field such that the components of the tensor are written (has only meaning with spherical component...
Definition: sym_tensor.h:274
Lorene::Scalar::check_dzpuis
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition: scalar.C:873
Lorene::abs
Cmp abs(const Cmp &)
Absolute value.
Definition: cmp_math.C:410
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::Sym_tensor_trans::set_AtB_trace
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.
Definition: sym_tensor_trans_aux.C:302
Lorene::Scalar::get_spectral_base
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va
Definition: scalar.h:1294
Lorene::Sym_tensor::p_tilde_b
Scalar * p_tilde_b
Field defined from and h insensitive to the longitudinal part of the Sym_tensor.
Definition: sym_tensor.h:334
Lorene::Sym_tensor_trans::p_trace
Scalar * p_trace
Trace with respect to the metric *met_div
Definition: sym_tensor.h:617
Lorene::Tensor::mp
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:295
Lorene::Sym_tensor_trans::trace_from_det_one
void trace_from_det_one(const Sym_tensor_tt &htt, double precis=1.e-14, int it_max=100)
Assigns the derived member p_tt and computes the trace so that *this + the flat metric has a determin...
Definition: sym_tensor_trans.C:315
Lorene::Base_vect_spher
Spherical orthonormal vectorial bases (triads).
Definition: base_vect.h:308
Lorene::Sym_tensor::p_aaa
Scalar * p_aaa
Field A defined from X and insensitive to the longitudinal part of the Sym_tensor (only for ).
Definition: sym_tensor.h:322
Lorene::Sym_tensor_trans::met_div
const Metric *const met_div
Metric with respect to which the divergence and the trace are defined.
Definition: sym_tensor.h:614
Lorene::Tensor::triad
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tensor.h:303
Lorene::Sym_tensor_trans::p_tt
Sym_tensor_tt * p_tt
Traceless part with respect to the metric *met_div
Definition: sym_tensor.h:620
Lorene::Scalar::set_etat_zero
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:324
Lorene::Sym_tensor_trans::set_tt_part_det_one
void set_tt_part_det_one(const Sym_tensor_tt &hijtt, const Scalar *h_prev=0x0, Param *par_mat=0x0, double precis=1.e-14, int it_max=100)
Assignes the TT-part of the tensor.
Definition: sym_tensor_trans_aux.C:226
Lorene::Sym_tensor::mu
const Scalar & mu(Param *par=0x0) const
Gives the field (see member p_mu ).
Definition: sym_tensor_aux.C:151
Lorene::Tensor::inc_dzpuis
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition: tensor.C:816
Lorene::Param
Parameter storage.
Definition: param.h:125
Lorene::Sym_tensor_tt::set_rr_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 ).
Definition: sym_tensor_tt_etamu.C:238
Lorene::Sym_tensor::xxx
const Scalar & xxx() const
Gives the field X (see member p_xxx ).
Definition: sym_tensor_aux.C:240
Lorene::Sym_tensor_trans::set_AtBtt_det_one
void set_AtBtt_det_one(const Scalar &a_in, const Scalar &tbtt_in, const Scalar *h_prev=0x0, Param *par_bc=0x0, Param *par_mat=0x0, double precis=1.e-14, int it_max=100)
Assigns the derived member A and computes from its TT-part (see Sym_tensor::compute_tilde_B_tt() ).
Definition: sym_tensor_trans_aux.C:137
Lorene::Sym_tensor_trans
Transverse symmetric tensors of rank 2.
Definition: sym_tensor.h:608
Lorene::Scalar::div_r
void div_r()
Division by r everywhere; dzpuis is not changed.
Definition: scalar_r_manip.C:133