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
29char 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
68namespace 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
106const 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
135const 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
174const 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
217const 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
260const 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}
void set(const Map *mp, Mtbl *(*construct)(const Map *))
Semi-constructor from a mapping and a method.
Definition coord.C:134
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Definition map.h:791
Coord y
y coordinate centered on the grid
Definition map.h:727
Coord ya
Absolute y coordinate.
Definition map.h:731
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
Coord x
x coordinate centered on the grid
Definition map.h:726
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition map.h:775
Coord xa
Absolute x coordinate.
Definition map.h:730
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
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
Metric for tensor calculation.
Definition metric.h:90
virtual const Vector & radial_vect() const
Returns the radial vector normal to a spherical slicing and pointing toward spatial infinity.
Definition metric.C:362
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition mg3d.C:473
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
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 Metric & get_gam() const
metric
Definition single_hor.C:339
const Valeur boundary_psi_app_hor() const
Neumann boundary condition for.
const Sym_tensor & get_k_dd() const
k_dd
Definition single_hor.C:348
Metric_flat ff
3 metric flat
Definition isol_hor.h:980
Scalar psi
Conformal factor .
Definition isol_hor.h:930
Map_af & mp
Affine mapping.
Definition isol_hor.h:900
Vector dn
Covariant derivative of the lapse with respect to the flat metric .
Definition isol_hor.h:937
Scalar nn
Lapse function .
Definition isol_hor.h:921
Scalar trK
Trace of the extrinsic curvature.
Definition isol_hor.h:989
Metric tgam
3 metric tilde
Definition isol_hor.h:977
Values and coefficients of a (real-value) function.
Definition valeur.h:287
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition valeur.h:363
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
Definition valeur.C:824
Tensor field of valence 1.
Definition vector.h:188
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:316
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
Definition vector.C:381
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition app_hor.h:64