LORENE
sym_ttt_poisson.C
1/*
2 * Resolution of the TT tensor Poisson equation
3 *
4 * (see file sym_tensor.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2003 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 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_ttt_poisson_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/sym_ttt_poisson.C,v 1.5 2014/10/13 08:53:44 j_novak Exp $" ;
29
30/*
31 * $Id: sym_ttt_poisson.C,v 1.5 2014/10/13 08:53:44 j_novak Exp $
32 * $Log: sym_ttt_poisson.C,v $
33 * Revision 1.5 2014/10/13 08:53:44 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.4 2004/12/28 10:37:24 j_novak
37 * Better way of enforcing zero divergence.
38 *
39 * Revision 1.3 2004/12/27 14:33:12 j_novak
40 * New algorithm for the tensor Poisson eq.
41 *
42 * Revision 1.2 2004/03/04 09:50:41 e_gourgoulhon
43 * Method poisson: use of new methods khi() and set_khi_mu.
44 *
45 * Revision 1.1 2003/11/07 16:53:52 e_gourgoulhon
46 * First version
47 *
48 *
49 * $Header: /cvsroot/Lorene/C++/Source/Tensor/sym_ttt_poisson.C,v 1.5 2014/10/13 08:53:44 j_novak Exp $
50 *
51 */
52
53// C headers
54//#include <>
55
56// Lorene headers
57#include "tensor.h"
58#include "param_elliptic.h"
59
60
61namespace Lorene {
63
64 // All this has a meaning only for spherical components...
65 assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
66 //## ... and affine mapping, for the moment!
67 assert(dynamic_cast<const Map_af*>(mp) != 0x0) ;
68 assert( (dzfin == 0) || (dzfin == 2) ) ;
70
71 // Solution for the rr-component
72 // ----------------------------
73
74 const Scalar& source_rr = operator()(1,1) ;
75 Scalar h_rr(*mp) ;
76 int nz = mp->get_mg()->get_nzone() ;
77
78 if (source_rr.get_etat() != ETATZERO) {
79
80 //------------------------
81 // The elliptic operator
82 //------------------------
83
85 for (int lz=0; lz<nz; lz++)
86 param_hr.set_poisson_tens_rr(lz) ;
87
88 h_rr = source_rr.sol_elliptic(param_hr) ;
89 }
90 else
91 h_rr.set_etat_zero() ;
92
93 h_rr.inc_dzpuis(dzfin) ; //## can we improve here?
94 resu.set(1,1) = h_rr ;
95
96 // Solution for (eta / r)
97 //-----------------------
98// Scalar source_eta = - source_rr ;
99// source_eta.mult_r_dzpuis(3) ;
100// source_eta.mult_r_dzpuis(2) ;
101// h_rr.set_spectral_va().ylm() ;
102// Scalar tmp = 2*h_rr + h_rr.lapang() ;
103// if (dzfin == 0)
104// tmp.inc_dzpuis(2) ;
105// source_eta += tmp ;
106// source_eta = source_eta.primr() ;
107
108// source_eta.div_r_dzpuis(dzfin) ;
109
110// Scalar etasurr = (h_rr+source_eta).poisson_angu() ;
111
112 Scalar source_eta = -resu(1,1).dsdr() ;
113 source_eta.mult_r_dzpuis(dzfin) ;
114 source_eta -= 3.*resu(1,1) ;
115 Scalar etasurr = source_eta.poisson_angu() ;
116
117 // Solution for mu
118 // ---------------
119
120 Scalar musurr = mu().poisson() ;
121 musurr.div_r_dzpuis(dzfin) ;
122
123 resu.set(1,1).set_spectral_va().ylm_i() ;
124
125 Scalar** rcmp = resu.cmp ;
126
127 Itbl idx(2) ;
128 idx.set(0) = 1 ; // r index
129
130 // h^{r theta} :
131 // ------------
132 idx.set(1) = 2 ; // theta index
133 *rcmp[position(idx)] = etasurr.dsdt() - musurr.stdsdp() ;
134
135 // h^{r phi} :
136 // ------------
137 idx.set(1) = 3 ; // phi index
138 *rcmp[position(idx)] = etasurr.stdsdp() + musurr.dsdt() ;
139
140 // h^{theta phi} and h^{phi phi}
141 // -----------------------------
142
143 //-------------- Computation of T^theta --> taut :
144
145 Scalar tautst = resu(1,2).dsdr() ;
146
147 // dhrr contains dh^{rt}/dr in all domains but the CED,
148 // in the CED: r^2 dh^{rt}/dr if dzfin = 0 (1)
149 // r^3 dh^{rt}/dr if dzfin = 2 (2)
150
151 // Multiplication by r of dh^{rt}/dr (with the same dzpuis than h^{rt})
152 tautst.mult_r_dzpuis(dzfin) ;
153
154 // Addition of the remaining parts :
155 tautst += 3 * resu(1,2) - resu(1,1).dsdt() ;
156 tautst.mult_sint() ;
157
158 Scalar tmp = resu(1,1) ;
159 tmp.mult_cost() ; // h^{rr} cos(th)
160
161 tautst -= tmp ; // T^th / sin(th)
162
163 Scalar taut = tautst ;
164 taut.mult_sint() ; // T^th
165
166
167 //----------- Computation of T^phi --> taup :
168
169 Scalar taupst = - resu(1,3).dsdr() ;
170
171 // dhrr contains - dh^{rp}/dr in all domains but the CED,
172 // in the CED: - r^2 dh^{rp}/dr if dzfin = 0 (3)
173 // - r^3 dh^{rp}/dr if dzfin = 2 (4)
174
175 // Multiplication by r of -dh^{rp}/dr (with the same dzpuis than h^{rp})
176 taupst.mult_r_dzpuis(dzfin) ;
177
178 // Addition of the remaining part :
179
180 taupst -= 3 * resu(1,3) ;
181 taupst.mult_sint() ; // T^ph / sin(th)
182
183 Scalar taup = taupst ;
184 taup.mult_sint() ; // T^ph
185
186 //------------------- Computation of F and h^[ph ph}
187
188 tmp = tautst ;
189 tmp.mult_cost() ; // T^th / tan(th)
190
191 // dT^th/dth + T^th / tan(th) + 1/sin(th) dT^ph/dph :
192 tmp = taut.dsdt() + tmp + taup.stdsdp() ;
193
194 Scalar tmp2 (*mp) ;
195 tmp2 = tmp.poisson_angu() ; // F
196 tmp2.div_sint() ;
197 tmp2.div_sint() ; // h^{ph ph}
198
199 idx.set(0) = 3 ; // phi index
200 idx.set(1) = 3 ; // phi index
201 *rcmp[position(idx)] = tmp2 ; // h^{ph ph} is updated
202
203
204 //------------------- Computation of G and h^[th ph}
205
206 tmp = taupst ;
207 tmp.mult_cost() ; // T^ph / tan(th)
208
209 // - 1/sin(th) dT^th/dph + dT^ph/dth + T^ph / tan(th) :
210 tmp = - taut.stdsdp() + taup.dsdt() + tmp ;
211
212 tmp2 = tmp.poisson_angu() ; // G
213 tmp2.div_sint() ;
214 tmp2.div_sint() ; // h^{th ph}
215
216 idx.set(0) = 2 ; // theta index
217 idx.set(1) = 3 ; // phi index
218 *rcmp[position(idx)] = tmp2 ; // h^{th ph} is updated
219
220 // h^{th th} (from the trace-free condition)
221 // ---------
222 idx.set(1) = 2 ; // theta index
223 *rcmp[position(idx)] = - resu(1,1) - resu(3,3) ;
224
225 return resu ;
226}
227}
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
This class contains the parameters needed to call the general elliptic solver.
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
const Metric *const met_div
Metric with respect to which the divergence and the trace are defined.
Definition sym_tensor.h:614
Transverse and traceless symmetric tensors of rank 2.
Definition sym_tensor.h:938
Sym_tensor_tt poisson(int dzfin=2) const
Computes the solution of a tensorial TT Poisson equation with *this as a source:
const Scalar & mu(Param *par=0x0) const
Gives the field (see member p_mu ).
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
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
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