LORENE
single_bound.C
1 /*
2  * Method of class single_hor to compute boundary conditions
3  *
4  * (see file isol_hor.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2004 Jose Luis Jaramillo
10  * Francois Limousin
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License version 2
16  * as published by the Free Software Foundation.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 char single_bound_C[] = "$Header: /cvsroot/Lorene/C++/Source/Isol_hor/single_bound.C,v 1.3 2014/10/13 08:53:01 j_novak Exp $" ;
30 
31 /*
32  * $Id: single_bound.C,v 1.3 2014/10/13 08:53:01 j_novak Exp $
33  * $Log: single_bound.C,v $
34  * Revision 1.3 2014/10/13 08:53:01 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.2 2014/10/06 15:13:11 j_novak
38  * Modified #include directives to use c++ syntax.
39  *
40  * Revision 1.1 2007/04/13 15:28:35 f_limousin
41  * Lots of improvements, generalisation to an arbitrary state of
42  * rotation, implementation of the spatial metric given by Samaya.
43  *
44  *
45  * $Header: /cvsroot/Lorene/C++/Source/Isol_hor/single_bound.C,v 1.3 2014/10/13 08:53:01 j_novak Exp $
46  *
47  */
48 
49 // C++ headers
50 #include "headcpp.h"
51 
52 // C headers
53 #include <cstdlib>
54 #include <cassert>
55 
56 // Lorene headers
57 #include "time_slice.h"
58 #include "isol_hor.h"
59 #include "metric.h"
60 #include "evolution.h"
61 #include "unites.h"
62 #include "graphique.h"
63 #include "utilitaires.h"
64 #include "param.h"
65 
66 
67 
68 namespace Lorene {
70 
71  Metric flat (mp.flat_met_spher()) ;
72  Vector temp (mp, CON, mp.get_bvect_spher()) ;
73  temp.set(1) = 1.;
74  temp.set(2) = 0.;
75  temp.set(3) = 0.;
76  temp.std_spectral_base() ;
77 
78  Scalar tmp = psi * psi * psi * trK
79  - contract(get_k_dd(),0, 1, tgam.radial_vect() * tgam.radial_vect(), 0, 1)
80  / psi
82  - 4 * ( tgam.radial_vect()(2) * psi.derive_cov(ff)(2)
83  + tgam.radial_vect()(3) * psi.derive_cov(ff)(3) ) ;
84 
85  tmp = tmp / (4 * tgam.radial_vect()(1)) ;
86 
87  // in this case you don't have to substract any value
88 
89  Valeur psi_bound (mp.get_mg()->get_angu() ) ;
90 
91  int nnp = mp.get_mg()->get_np(1) ;
92  int nnt = mp.get_mg()->get_nt(1) ;
93 
94  psi_bound = 1 ;
95 
96  for (int k=0 ; k<nnp ; k++)
97  for (int j=0 ; j<nnt ; j++)
98  psi_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
99 
100  psi_bound.std_base_scal() ;
101 
102  return psi_bound ;
103 
104 }
105 
106 const Valeur Single_hor::boundary_nn_Neu(double cc)const {
107 
108  double rho = 1. ; // 1 is the standart case;
109 
110  Scalar tmp = - cc * nn ;
111  // Scalar tmp = - nn()/psi()*psi().dsdr() ;
112 
113  // in this case you don't have to substract any value
114  tmp += (rho - 1) * tgam.radial_vect()(1) * dn(1) ;
115  tmp = tmp / (rho * tgam.radial_vect()(1)) ;
116 
117  int nnp = mp.get_mg()->get_np(1) ;
118  int nnt = mp.get_mg()->get_nt(1) ;
119 
120  Valeur nn_bound (mp.get_mg()->get_angu()) ;
121 
122  nn_bound = 1 ; // Juste pour affecter dans espace des configs ;
123 
124  for (int k=0 ; k<nnp ; k++)
125  for (int j=0 ; j<nnt ; j++)
126  nn_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
127 
128  nn_bound.std_base_scal() ;
129 
130  return nn_bound ;
131 
132 }
133 
134 
135 const Valeur Single_hor::boundary_nn_Dir(double cc)const {
136 
137  Scalar rho(mp);
138  rho = 0. ; // 0 is the standard case
139 
140  Scalar tmp(mp) ;
141  tmp = cc ;
142 
143 
144  //tmp = 1./(2*psi()) ;
145  // tmp = - psi() * nn().dsdr() / (psi().dsdr()) ;
146 
147  // We have substracted 1, since we solve for zero condition at infinity
148  //and then we add 1 to the solution
149 
150  tmp = (tmp + rho * nn)/(1 + rho) ;
151 
152  tmp = tmp - 1 ;
153 
154  int nnp = mp.get_mg()->get_np(1) ;
155  int nnt = mp.get_mg()->get_nt(1) ;
156 
157  Valeur nn_bound (mp.get_mg()->get_angu()) ;
158 
159  nn_bound = 1 ; // Juste pour affecter dans espace des configs ;
160 
161  for (int k=0 ; k<nnp ; k++)
162  for (int j=0 ; j<nnt ; j++)
163  nn_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
164 
165  nn_bound.std_base_scal() ;
166 
167  return nn_bound ;
168 
169 }
170 
171 // Component x of boundary value of beta
172 //--------------------------------------
173 
174 const Valeur Single_hor:: boundary_beta_x(double om_orb,
175  double om_loc)const {
176 
177  // Les alignemenents pour le signe des CL.
178  double orientation = mp.get_rot_phi() ;
179  assert ((orientation == 0) || (orientation == M_PI)) ;
180  int aligne = (orientation == 0) ? 1 : -1 ;
181 
182  int nnp = mp.get_mg()->get_np(1) ;
183  int nnt = mp.get_mg()->get_nt(1) ;
184 
185  Vector tmp_vect = nn * get_gam().radial_vect() ;
186  tmp_vect.change_triad(mp.get_bvect_cart() ) ;
187 
188  //Isol_hor boundary conditions
189 
190  Valeur lim_x (mp.get_mg()->get_angu()) ;
191 
192  lim_x = 1 ; // Juste pour affecter dans espace des configs ;
193 
194  Mtbl ya_mtbl (mp.get_mg()) ;
195  ya_mtbl.set_etat_qcq() ;
196  ya_mtbl = mp.ya ;
197 
198  Mtbl yy_mtbl (mp.get_mg()) ;
199  yy_mtbl.set_etat_qcq() ;
200  yy_mtbl = mp.y ;
201 
202  for (int k=0 ; k<nnp ; k++)
203  for (int j=0 ; j<nnt ; j++)
204  lim_x.set(0, k, j, 0) = aligne * om_orb * ya_mtbl(1, k, j, 0)
205  + (om_loc-om_orb)* yy_mtbl(1, k, j, 0)
206  + tmp_vect(1).val_grid_point(1, k, j, 0) ;
207 
208  lim_x.set_base(*(mp.get_mg()->std_base_vect_cart()[0])) ;
209 
210  return lim_x ;
211 }
212 
213 
214 // Component y of boundary value of beta
215 //--------------------------------------
216 
217 const Valeur Single_hor:: boundary_beta_y(double om_orb,
218  double om_loc)const {
219 
220  // Les alignemenents pour le signe des CL.
221  double orientation = mp.get_rot_phi() ;
222  assert ((orientation == 0) || (orientation == M_PI)) ;
223  int aligne = (orientation == 0) ? 1 : -1 ;
224 
225 
226  int nnp = mp.get_mg()->get_np(1) ;
227  int nnt = mp.get_mg()->get_nt(1) ;
228 
229  Vector tmp_vect = nn * get_gam().radial_vect() ;
230  tmp_vect.change_triad(mp.get_bvect_cart() ) ;
231 
232  //Isol_hor boundary conditions
233 
234  Valeur lim_y (mp.get_mg()->get_angu()) ;
235 
236  lim_y = 1 ; // Juste pour affecter dans espace des configs ;
237 
238  Mtbl xa_mtbl (mp.get_mg()) ;
239  xa_mtbl.set_etat_qcq() ;
240  xa_mtbl = mp.xa ;
241 
242  Mtbl xx_mtbl (mp.get_mg()) ;
243  xx_mtbl.set_etat_qcq() ;
244  xx_mtbl = mp.x ;
245 
246  for (int k=0 ; k<nnp ; k++)
247  for (int j=0 ; j<nnt ; j++)
248  lim_y.set(0, k, j, 0) = - aligne *om_orb * xa_mtbl(1, k, j, 0) -
249  (om_loc-om_orb)*xx_mtbl(1, k, j, 0)
250  + tmp_vect(2).val_grid_point(1, k, j, 0) ;
251 
252  lim_y.set_base(*(mp.get_mg()->std_base_vect_cart()[1])) ;
253 
254  return lim_y ;
255 }
256 
257 // Component z of boundary value of beta
258 //--------------------------------------
259 
260 const Valeur Single_hor:: boundary_beta_z()const {
261 
262  int nnp = mp.get_mg()->get_np(1) ;
263  int nnt = mp.get_mg()->get_nt(1) ;
264 
265  Vector tmp_vect = nn * get_gam().radial_vect() ;
266  tmp_vect.change_triad(mp.get_bvect_cart() ) ;
267 
268  //Isol_hor boundary conditions
269 
270  Valeur lim_z (mp.get_mg()->get_angu()) ;
271 
272  lim_z = 1 ; // Juste pour affecter dans espace des configs ;
273 
274  for (int k=0 ; k<nnp ; k++)
275  for (int j=0 ; j<nnt ; j++)
276  lim_z.set(0, k, j, 0) = tmp_vect(3).val_grid_point(1, k, j, 0) ;
277 
278  lim_z.set_base(*(mp.get_mg()->std_base_vect_cart()[2])) ;
279 
280  return lim_z ;
281 }
282 
283 }
Coord xa
Absolute x coordinate.
Definition: map.h:730
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
Metric for tensor calculation.
Definition: metric.h:90
const Metric & get_gam() const
metric
Definition: single_hor.C:339
Vector dn
Covariant derivative of the lapse with respect to the flat metric .
Definition: isol_hor.h:937
Metric tgam
3 metric tilde
Definition: isol_hor.h:977
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
Multi-domain array.
Definition: mtbl.h:118
const Sym_tensor & get_k_dd() const
k_dd
Definition: single_hor.C:348
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:783
Lorene prototypes.
Definition: app_hor.h:64
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
virtual const Vector & radial_vect() const
Returns the radial vector normal to a spherical slicing and pointing toward spatial infinity...
Definition: metric.C:362
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:784
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition: map.h:775
void set(const Map *mp, Mtbl *(*construct)(const Map *))
Semi-constructor from a mapping and a method.
Definition: coord.C:134
Tensor field of valence 1.
Definition: vector.h:188
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition: scalar.h:637
const Valeur boundary_psi_app_hor() const
Neumann boundary condition for.
Definition: single_bound.C:69
Scalar psi
Conformal factor .
Definition: isol_hor.h:930
Scalar trK
Trace of the extrinsic curvature.
Definition: isol_hor.h:989
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition: mg3d.C:473
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition: vector.C:381
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl.C:299
Coord ya
Absolute y coordinate.
Definition: map.h:731
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:791
Map_af & mp
Affine mapping.
Definition: isol_hor.h:900
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:654
Coord y
y coordinate centered on the grid
Definition: map.h:727
Coord x
x coordinate centered on the grid
Definition: map.h:726
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
Metric_flat ff
3 metric flat
Definition: isol_hor.h:980
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this.
Definition: scalar_deriv.C:390
const Metric_flat & flat_met_spher() const
Returns the flat metric associated with the spherical coordinates and with components expressed in th...
Definition: map.C:321
Scalar nn
Lapse function .
Definition: isol_hor.h:921
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition: valeur.h:363