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
28char 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
116namespace Lorene {
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
128 tens_tt.set_rr_mu(hrr, mu_in) ;
129 tens_tt.inc_dzpuis(2) ;
131 dec_dzpuis(2) ;
132
133 return ;
134
135}
136
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) ;
147 assert( (precis > 0.) && (it_max > 0) ) ;
148
149 //Computation of mu and X from A
150 //-------------------------------
152 Scalar x_new(*mp) ;
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) ;
168 Scalar w_tt(*mp) ;
171 hijtt.set_auxiliary(hrr_tt, eta_sr_tt, mu_over_r, w_tt, x_new, zero) ;
172
173 Scalar hrr_new(*mp) ;
175 Scalar w_new(*mp) ;
176
177 for (int it=0; it<=it_max; it++) {
180
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) {
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 ;
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
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) ;
255 Scalar w_new(*mp) ;
256
257 for (int it=0; it<=it_max; it++) {
260
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) {
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
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 //-------------------------------
316 Scalar x_new(*mp) ;
318
319 // Computation of the other potentials
320 //------------------------------------
321 Scalar hrr_new(*mp) ;
323 Scalar w_new(*mp) ;
324
326
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}
Spherical orthonormal vectorial bases (triads).
Definition base_vect.h:308
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Parameter storage.
Definition param.h:125
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
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
Transverse symmetric tensors of rank 2.
Definition sym_tensor.h:608
Sym_tensor_tt * p_tt
Traceless part with respect to the metric *met_div
Definition sym_tensor.h:620
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 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...
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...
const Metric *const met_div
Metric with respect to which the divergence and the trace are defined.
Definition sym_tensor.h:614
Scalar * p_trace
Trace with respect to the metric *met_div
Definition sym_tensor.h:617
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_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.
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...
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 .
Transverse and traceless symmetric tensors of rank 2.
Definition sym_tensor.h:938
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
Scalar * p_tilde_b
Field defined from and h insensitive to the longitudinal part of the Sym_tensor.
Definition sym_tensor.h:334
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 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.
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 ,...
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
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
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
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
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