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 
28 
29 
30 /*
31  * $Id: sym_ttt_poisson.C,v 1.6 2016/12/05 16:18:17 j_novak Exp $
32  * $Log: sym_ttt_poisson.C,v $
33  * Revision 1.6 2016/12/05 16:18:17 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.5 2014/10/13 08:53:44 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.4 2004/12/28 10:37:24 j_novak
40  * Better way of enforcing zero divergence.
41  *
42  * Revision 1.3 2004/12/27 14:33:12 j_novak
43  * New algorithm for the tensor Poisson eq.
44  *
45  * Revision 1.2 2004/03/04 09:50:41 e_gourgoulhon
46  * Method poisson: use of new methods khi() and set_khi_mu.
47  *
48  * Revision 1.1 2003/11/07 16:53:52 e_gourgoulhon
49  * First version
50  *
51  *
52  * $Header: /cvsroot/Lorene/C++/Source/Tensor/sym_ttt_poisson.C,v 1.6 2016/12/05 16:18:17 j_novak Exp $
53  *
54  */
55 
56 // C headers
57 //#include <>
58 
59 // Lorene headers
60 #include "tensor.h"
61 #include "param_elliptic.h"
62 
63 
64 namespace Lorene {
66 
67  // All this has a meaning only for spherical components...
68  assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
69  //## ... and affine mapping, for the moment!
70  assert(dynamic_cast<const Map_af*>(mp) != 0x0) ;
71  assert( (dzfin == 0) || (dzfin == 2) ) ;
72  Sym_tensor_tt resu(*mp, *triad, *met_div) ;
73 
74  // Solution for the rr-component
75  // ----------------------------
76 
77  const Scalar& source_rr = operator()(1,1) ;
78  Scalar h_rr(*mp) ;
79  int nz = mp->get_mg()->get_nzone() ;
80 
81  if (source_rr.get_etat() != ETATZERO) {
82 
83  //------------------------
84  // The elliptic operator
85  //------------------------
86 
87  Param_elliptic param_hr(source_rr) ;
88  for (int lz=0; lz<nz; lz++)
89  param_hr.set_poisson_tens_rr(lz) ;
90 
91  h_rr = source_rr.sol_elliptic(param_hr) ;
92  }
93  else
94  h_rr.set_etat_zero() ;
95 
96  h_rr.inc_dzpuis(dzfin) ; //## can we improve here?
97  resu.set(1,1) = h_rr ;
98 
99  // Solution for (eta / r)
100  //-----------------------
101 // Scalar source_eta = - source_rr ;
102 // source_eta.mult_r_dzpuis(3) ;
103 // source_eta.mult_r_dzpuis(2) ;
104 // h_rr.set_spectral_va().ylm() ;
105 // Scalar tmp = 2*h_rr + h_rr.lapang() ;
106 // if (dzfin == 0)
107 // tmp.inc_dzpuis(2) ;
108 // source_eta += tmp ;
109 // source_eta = source_eta.primr() ;
110 
111 // source_eta.div_r_dzpuis(dzfin) ;
112 
113 // Scalar etasurr = (h_rr+source_eta).poisson_angu() ;
114 
115  Scalar source_eta = -resu(1,1).dsdr() ;
116  source_eta.mult_r_dzpuis(dzfin) ;
117  source_eta -= 3.*resu(1,1) ;
118  Scalar etasurr = source_eta.poisson_angu() ;
119 
120  // Solution for mu
121  // ---------------
122 
123  Scalar musurr = mu().poisson() ;
124  musurr.div_r_dzpuis(dzfin) ;
125 
126  resu.set(1,1).set_spectral_va().ylm_i() ;
127 
128  Scalar** rcmp = resu.cmp ;
129 
130  Itbl idx(2) ;
131  idx.set(0) = 1 ; // r index
132 
133  // h^{r theta} :
134  // ------------
135  idx.set(1) = 2 ; // theta index
136  *rcmp[position(idx)] = etasurr.dsdt() - musurr.stdsdp() ;
137 
138  // h^{r phi} :
139  // ------------
140  idx.set(1) = 3 ; // phi index
141  *rcmp[position(idx)] = etasurr.stdsdp() + musurr.dsdt() ;
142 
143  // h^{theta phi} and h^{phi phi}
144  // -----------------------------
145 
146  //-------------- Computation of T^theta --> taut :
147 
148  Scalar tautst = resu(1,2).dsdr() ;
149 
150  // dhrr contains dh^{rt}/dr in all domains but the CED,
151  // in the CED: r^2 dh^{rt}/dr if dzfin = 0 (1)
152  // r^3 dh^{rt}/dr if dzfin = 2 (2)
153 
154  // Multiplication by r of dh^{rt}/dr (with the same dzpuis than h^{rt})
155  tautst.mult_r_dzpuis(dzfin) ;
156 
157  // Addition of the remaining parts :
158  tautst += 3 * resu(1,2) - resu(1,1).dsdt() ;
159  tautst.mult_sint() ;
160 
161  Scalar tmp = resu(1,1) ;
162  tmp.mult_cost() ; // h^{rr} cos(th)
163 
164  tautst -= tmp ; // T^th / sin(th)
165 
166  Scalar taut = tautst ;
167  taut.mult_sint() ; // T^th
168 
169 
170  //----------- Computation of T^phi --> taup :
171 
172  Scalar taupst = - resu(1,3).dsdr() ;
173 
174  // dhrr contains - dh^{rp}/dr in all domains but the CED,
175  // in the CED: - r^2 dh^{rp}/dr if dzfin = 0 (3)
176  // - r^3 dh^{rp}/dr if dzfin = 2 (4)
177 
178  // Multiplication by r of -dh^{rp}/dr (with the same dzpuis than h^{rp})
179  taupst.mult_r_dzpuis(dzfin) ;
180 
181  // Addition of the remaining part :
182 
183  taupst -= 3 * resu(1,3) ;
184  taupst.mult_sint() ; // T^ph / sin(th)
185 
186  Scalar taup = taupst ;
187  taup.mult_sint() ; // T^ph
188 
189  //------------------- Computation of F and h^[ph ph}
190 
191  tmp = tautst ;
192  tmp.mult_cost() ; // T^th / tan(th)
193 
194  // dT^th/dth + T^th / tan(th) + 1/sin(th) dT^ph/dph :
195  tmp = taut.dsdt() + tmp + taup.stdsdp() ;
196 
197  Scalar tmp2 (*mp) ;
198  tmp2 = tmp.poisson_angu() ; // F
199  tmp2.div_sint() ;
200  tmp2.div_sint() ; // h^{ph ph}
201 
202  idx.set(0) = 3 ; // phi index
203  idx.set(1) = 3 ; // phi index
204  *rcmp[position(idx)] = tmp2 ; // h^{ph ph} is updated
205 
206 
207  //------------------- Computation of G and h^[th ph}
208 
209  tmp = taupst ;
210  tmp.mult_cost() ; // T^ph / tan(th)
211 
212  // - 1/sin(th) dT^th/dph + dT^ph/dth + T^ph / tan(th) :
213  tmp = - taut.stdsdp() + taup.dsdt() + tmp ;
214 
215  tmp2 = tmp.poisson_angu() ; // G
216  tmp2.div_sint() ;
217  tmp2.div_sint() ; // h^{th ph}
218 
219  idx.set(0) = 2 ; // theta index
220  idx.set(1) = 3 ; // phi index
221  *rcmp[position(idx)] = tmp2 ; // h^{th ph} is updated
222 
223  // h^{th th} (from the trace-free condition)
224  // ---------
225  idx.set(1) = 2 ; // theta index
226  *rcmp[position(idx)] = - resu(1,1) - resu(3,3) ;
227 
228  return resu ;
229 }
230 }
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 ).
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Definition: itbl.h:247
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:134
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
const Scalar & dsdt() const
Returns of *this .
Definition: scalar_deriv.C:208
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:777
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition: scalar_pde.C:139
void mult_sint()
Multiplication by .
Basic integer array class.
Definition: itbl.h:122
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tensor.h:309
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:560
Scalar poisson_angu(double lambda=0) const
Solves the (generalized) angular Poisson equation with *this as source.
Definition: scalar_pde.C:203
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition: tensor.h:321
Scalar sol_elliptic(Param_elliptic &params) const
Resolution of a general elliptic equation, putting zero at infinity.
Definition: scalar_pde.C:237
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
void set_poisson_tens_rr(int zone)
Sets the operator to in all domains, for only.
const Scalar & stdsdp() const
Returns of *this .
Definition: scalar_deriv.C:238
void mult_cost()
Multiplication by .
This class contains the parameters needed to call the general elliptic solver.
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:663
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:248
const Scalar & operator()(const Itbl &ind) const
Returns the value of a component (read-only version).
Definition: tensor.C:807
void div_sint()
Division by .
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:610
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition: tensor.h:301
void div_r_dzpuis(int ced_mult_r)
Division by r everywhere but with the output flag dzpuis set to ced_mult_r .
void mult_r_dzpuis(int ced_mult_r)
Multiplication by r everywhere but with the output flag dzpuis set to ced_mult_r ...
Transverse and traceless symmetric tensors of rank 2.
Definition: sym_tensor.h:933