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