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 char 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 
61 namespace 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) ) ;
69  Sym_tensor_tt resu(*mp, *triad, *met_div) ;
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 
84  Param_elliptic param_hr(source_rr) ;
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 }
Lorene::Scalar::inc_dzpuis
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...
Definition: scalar_r_manip.C:470
Lorene::Itbl::set
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Definition: itbl.h:247
Lorene::Scalar::sol_elliptic
Scalar sol_elliptic(Param_elliptic &params) const
Resolution of a general elliptic equation, putting zero at infinity.
Definition: scalar_pde.C:234
Lorene::Sym_tensor_tt
Transverse and traceless symmetric tensors of rank 2.
Definition: sym_tensor.h:938
Lorene::Tensor::operator()
const Scalar & operator()(const Itbl &ind) const
Returns the value of a component (read-only version).
Definition: tensor.C:798
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Tensor::cmp
Scalar ** cmp
Array of size n_comp of pointers onto the components.
Definition: tensor.h:315
Lorene::Scalar::mult_r_dzpuis
void mult_r_dzpuis(int ced_mult_r)
Multiplication by r everywhere but with the output flag dzpuis set to ced_mult_r .
Definition: scalar_r_manip.C:221
Lorene::Itbl
Basic integer array class.
Definition: itbl.h:122
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::Param_elliptic
This class contains the parameters needed to call the general elliptic solver.
Definition: param_elliptic.h:135
Lorene::Scalar::div_sint
void div_sint()
Division by .
Definition: scalar_th_manip.C:98
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
Lorene::Tensor_sym::position
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
Lorene::Tensor::set
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:654
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_tt::poisson
Sym_tensor_tt poisson(int dzfin=2) const
Computes the solution of a tensorial TT Poisson equation with *this as a source:
Definition: sym_ttt_poisson.C:62
Lorene::Base_vect_spher
Spherical orthonormal vectorial bases (triads).
Definition: base_vect.h:308
Lorene::Scalar::div_r_dzpuis
void div_r_dzpuis(int ced_mult_r)
Division by r everywhere but with the output flag dzpuis set to ced_mult_r .
Definition: scalar_r_manip.C:147
Lorene::Scalar::mult_sint
void mult_sint()
Multiplication by .
Definition: scalar_th_manip.C:84
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::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Tensor::triad
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition: tensor.h:303
Lorene::Map_af
Affine radial mapping.
Definition: map.h:2027
Lorene::Scalar::set_etat_zero
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:324
Lorene::Scalar::dsdt
const Scalar & dsdt() const
Returns of *this .
Definition: scalar_deriv.C:208
Lorene::Scalar::get_etat
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
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::Scalar::poisson
Scalar poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition: scalar_pde.C:136
Lorene::Scalar::poisson_angu
Scalar poisson_angu(double lambda=0) const
Solves the (generalized) angular Poisson equation with *this as source.
Definition: scalar_pde.C:200
Lorene::Scalar::stdsdp
const Scalar & stdsdp() const
Returns of *this .
Definition: scalar_deriv.C:238
Lorene::Param_elliptic::set_poisson_tens_rr
void set_poisson_tens_rr(int zone)
Sets the operator to in all domains, for only.
Definition: param_elliptic.C:560
Lorene::Scalar::mult_cost
void mult_cost()
Multiplication by .
Definition: scalar_th_manip.C:58
Lorene::Valeur::ylm_i
void ylm_i()
Inverse of ylm()
Definition: valeur_ylm_i.C:131