31char CTS_gen[] =
"$Header: /cvsroot/Lorene/C++/Source/Isol_hor/CTS_gen.C,v 1.7 2014/10/13 08:53:00 j_novak Exp $" ;
34#include "time_slice.h"
40#include "utilitaires.h"
44void Isol_hor::init_data_CTS_gen(
int,
double,
int,
int,
45 int solve_lapse,
int solve_psi,
int solve_shift,
46 double precis,
double relax_nn ,
47 double relax_psi,
double relax_beta ,
48 int niter,
double a,
double ) {
61 ofstream conv(
"resconv.d") ;
62 ofstream kss_out(
"kss.d") ;
63 conv <<
" # diff_nn diff_psi diff_beta " << endl ;
66 Scalar nntilde_j =
nn() *
pow(
psi(), a) ;
67 nntilde_j.std_spectral_base() ;
68 Scalar psi_j =
psi() ;
69 Vector beta_j =
beta() ;
75 for (
int mer=0; mer<niter; mer++) {
80 Scalar temp_scal_0 (
mp) ;
81 Scalar temp_scal_2 (
mp) ;
82 Scalar temp_scal_3 (
mp) ;
83 Scalar temp_scal_4 (
mp) ;
85 Vector temp_vect_2 (
mp, CON, triad) ;
86 Vector temp_vect_3 (
mp, CON, triad) ;
87 Vector temp_vect_4 (
mp, CON, triad) ;
89 Scalar psi_j_a =
pow(psi_j, a) ;
90 psi_j_a.std_spectral_base() ;
91 Scalar psi_j_ma =
pow(psi_j, -a) ;
92 psi_j_ma.std_spectral_base() ;
93 Vector beta_j_d = beta_j.down(0,
met_gamt) ;
94 Scalar nnpsia_j = nntilde_j * psi_j_a /( psi_j*psi_j*psi_j*psi_j*psi_j*psi_j) ;
105 cout <<
" relax_nn = " << relax_nn << endl ;
106 cout <<
" relax_psi = " << relax_psi << endl ;
107 cout <<
" relax_beta = " << relax_beta << endl ;
116 Scalar sour_psi (
mp) ;
122 temp_scal_4 = 1./12. *
trK*
trK * psi_j*psi_j*psi_j*psi_j*psi_j
123 - 1./32. * psi_j*psi_j*psi_j*psi_j*psi_j * psi_j_ma*psi_j_ma / (nntilde_j*nntilde_j) *
126 sour_psi = temp_scal_3 + temp_scal_4 ;
132 temp_scal_3.inc_dzpuis() ;
135 sour_psi += temp_scal_3 + temp_scal_4 ;
137 sour_psi.annule_domain(0) ;
142 Scalar sour_nntilde (
mp) ;
146 temp_scal_2 = - psi_j*psi_j*psi_j*psi_j * psi_j_ma *
trK_point ;
153 temp_scal_4 = nntilde_j * (4-a)/12. * psi_j*psi_j*psi_j*psi_j *
trK*
trK
154 - 2 * (a+1) *
contract(psi_j.derive_cov(
ff), 0, nntilde_j.derive_con(
met_gamt), 0) / psi_j
156 + (a+8)/32. * psi_j*psi_j*psi_j*psi_j * psi_j_ma*psi_j_ma/nntilde_j
159 sour_nntilde = temp_scal_2 + temp_scal_3 + temp_scal_4 ;
164 temp_scal_3.inc_dzpuis() ;
168 sour_nntilde += temp_scal_3 + temp_scal_4 ;
169 sour_nntilde.annule_domain(0) ;
174 Vector sour_beta(
mp, CON, triad) ;
182 temp_vect_4 =
contract(beta_j.ope_killing_conf(
met_gamt), 1, nnpsia_j.derive_cov(
ff), 0) / nnpsia_j ;
184 sour_beta = temp_vect_3 + temp_vect_4 ;
189 temp_vect_3 = (lambda - 1./3.) *
contract (beta_j.derive_cov(
ff).derive_con(
ff), 0, 1)
198 sour_beta += temp_vect_3 + temp_vect_4 ;
207 Scalar tmp = psi_j * psi_j * psi_j *
trK
209 - psi_j*psi_j*psi_j*psi_j_ma/(2. * nntilde_j)
227 for (
int k=0 ; k<nnp ; k++)
228 for (
int j=0 ; j<nnt ; j++)
229 psi_bound.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
231 psi_bound.std_base_scal() ;
237 const Coord& y =
mp.
y ;
238 const Coord& x =
mp.
x ;
242 xx.std_spectral_base() ;
246 yy.std_spectral_base() ;
257 for (
int k=0 ; k<nnp ; k++)
258 for (
int j=0 ; j<nnt ; j++)
259 lim_x.set(0, k, j, 0) =
omega * yy.val_grid_point(1, k, j, 0)
260 + tmp_vect(1).val_grid_point(1, k, j, 0) ;
271 for (
int k=0 ; k<nnp ; k++)
272 for (
int j=0 ; j<nnt ; j++)
273 lim_y.set(0, k, j, 0) = -
omega * xx.val_grid_point(1, k, j, 0)
274 + tmp_vect(2).val_grid_point(1, k, j, 0) ;
284 for (
int k=0 ; k<nnp ; k++)
285 for (
int j=0 ; j<nnt ; j++)
286 lim_z.set(0, k, j, 0) = tmp_vect(3).val_grid_point(1, k, j, 0) ;
298 Scalar btilde_j =
contract (beta_j_d, 0, stilde_j, 0) ;
303 hh_tilde.dec_dzpuis(2) ;
306 tmp_vect = btilde_j * stilde_j ;
307 Vector VV_j = btilde_j * stilde_j - beta_j ;
310 Scalar accel = 2 *
contract( VV_j, 0,
contract( stilde_j, 0, stilde_j.down(0,
315 Sym_tensor qq_spher =
met_gamt.
con() - stilde_j * stilde_j ;
325 angular.set(2) =
omega * xxa ;
326 angular.set(3).annule_hard() ;
329 angular.set(1).set_spectral_va()
331 angular.set(2).set_spectral_va()
333 angular.set(3).set_spectral_va()
340 Scalar corr_nn_kappa (
mp) ;
343 corr_nn_kappa =
sqrt(
sqrt (corr_nn_kappa*corr_nn_kappa)) ;
344 corr_nn_kappa.std_spectral_base() ;
347 Scalar kss = -
kappa_hor() * psi_j_ma / nntilde_j ;
349 kss +=
contract(stilde_j, 0, nntilde_j.derive_cov(
ff), 0)/ (psi_j*psi_j*nntilde_j)
350 + a *
contract(stilde_j, 0, psi_j.derive_cov(
ff), 0) / (psi_j*psi_j*psi_j)
358 Scalar beta_r_corr = (rho - 1) * btilde_j * hh_tilde;
360 Scalar nn_KK = nntilde_j * psi_j_a *
trK ;
363 beta_r_corr.set_dzpuis(2) ;
364 nn_KK.set_dzpuis(2) ;
365 accel.set_dzpuis(2) ;
366 div_VV.set_dzpuis(2) ;
368 Scalar b_tilde_new (
mp) ;
369 b_tilde_new = 2 *
contract(stilde_j, 0, btilde_j.derive_cov(
ff), 0)
371 - 3 * nntilde_j * psi_j_a * kss + nn_KK + accel + div_VV ;
373 b_tilde_new = b_tilde_new / (hh_tilde * rho) ;
377 tmp_vect.set(1) = b_tilde_new * stilde_j(1) + angular(1) ;
378 tmp_vect.set(2) = b_tilde_new * stilde_j(2) + angular(2) ;
379 tmp_vect.set(3) = b_tilde_new * stilde_j(3) + angular(3) ;
388 for (
int k=0 ; k<nnp ; k++)
389 for (
int j=0 ; j<nnt ; j++)
390 lim_x_AEI.set(0, k, j, 0) = tmp_vect(1).val_grid_point(1, k, j, 0) ;
400 for (
int k=0 ; k<nnp ; k++)
401 for (
int j=0 ; j<nnt ; j++)
402 lim_y_AEI.set(0, k, j, 0) = tmp_vect(2).val_grid_point(1, k, j, 0) ;
411 for (
int k=0 ; k<nnp ; k++)
412 for (
int j=0 ; j<nnt ; j++)
413 lim_z_AEI.set(0, k, j, 0) = tmp_vect(3).val_grid_point(1, k, j, 0) ;
430 Scalar darea = psi_j* psi_j* psi_j* psi_j
434 darea.std_spectral_base() ;
436 Scalar area_int (darea) ;
437 area_int.raccord(1) ;
441 double radius = area / (4. * M_PI);
447 tmp.std_spectral_base() ;
453 Scalar kksphi = psi_j*psi_j*psi_j_ma/(2.*nntilde_j) *
455 1./3. *
trK * psi_j*psi_j *
458 kksphi = kksphi / (8. * M_PI) ;
459 Scalar kkspsi_int = kksphi * darea ;
476 kappa.std_spectral_base() ;
477 kappa.inc_dzpuis(2) ;
481 tmp = - a / psi_j * nntilde_j
483 + 1./2. * psi_j * psi_j * psi_j_ma
486 + 1./3. * nntilde_j * psi_j * psi_j *
trK
488 + psi_j * psi_j * psi_j_ma * kappa
501 for (
int k=0 ; k<nnp ; k++)
502 for (
int j=0 ; j<nnt ; j++)
503 nn_bound_kappa.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
505 nn_bound_kappa.std_base_scal() ;
540 tmp = 0.2 * psi_j_ma ;
541 tmp = (tmp + rho * nntilde_j)/(1 + rho) ;
549 for (
int k=0 ; k<nnp ; k++)
550 for (
int j=0 ; j<nnt ; j++)
551 nn_bound_const.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
553 nn_bound_const.std_base_scal() ;
560 tmp = btilde_j * psi_j*psi_j*psi_j_ma ;
561 tmp = (tmp + rho * nntilde_j)/(1 + rho) ;
569 for (
int k=0 ; k<nnp ; k++)
570 for (
int j=0 ; j<nnt ; j++)
571 nn_bound_AEI.set(0, k, j, 0) = tmp.val_grid_point(1, k, j, 0) ;
573 nn_bound_AEI.std_base_scal() ;
582 Scalar nntilde_jp1(
mp) ;
583 Vector beta_jp1(beta_j) ;
586 if (solve_psi == 1) {
587 psi_jp1 = sour_psi.poisson_neumann(psi_bound, 0) + 1. ;
590 maxabs(psi_jp1.laplacian() - sour_psi,
591 "Absolute error in the resolution of the equation for Psi") ;
593 psi_jp1 = relax_psi * psi_jp1 + (1 - relax_psi) * psi_j ;
597 if (solve_lapse == 1) {
601 nntilde_jp1 = sour_nntilde.poisson_dirichlet(nn_bound_AEI, 0) + 1. ;
604 maxabs(nntilde_jp1.laplacian() - sour_nntilde,
605 "Absolute error in the resolution of the equation for nntilde") ;
607 nntilde_jp1 = relax_nn * nntilde_jp1 + (1 - relax_nn) * nntilde_j ;
610 if (solve_shift == 1) {
611 double precision = 1e-8 ;
614 poisson_vect_boundary(lambda, sour_beta, beta_jp1, lim_x_AEI,
615 lim_y_AEI, lim_z_AEI, 0, precision, 20) ;
619 sour_beta.dec_dzpuis() ;
621 + lambda * beta_jp1.divergence(
ff)
622 .derive_con(
ff) - sour_beta,
623 "Absolute error in the resolution of the equation for beta") ;
628 beta_jp1 = relax_beta * beta_jp1 + (1 - relax_beta) * beta_j ;
637 double diff_nn, diff_psi, diff_beta ;
641 if (solve_lapse == 1)
642 diff_nn =
max(
diffrel(nntilde_j, nntilde_jp1) ) ;
645 if (solve_shift == 1)
646 diff_beta =
max(
maxabs(beta_jp1 - beta_j) ) ;
648 cout <<
"step = " << mer <<
" : diff_psi = " << diff_psi
649 <<
" diff_nn = " << diff_nn
650 <<
" diff_beta = " << diff_beta << endl ;
651 cout <<
"----------------------------------------------" << endl ;
652 if ((diff_psi<precis) && (diff_nn<precis) && (diff_beta<precis))
655 if (mer>0) {conv << mer <<
" " <<
log10(diff_nn) <<
" " <<
log10(diff_psi)
656 <<
" " <<
log10(diff_beta) << endl ; } ;
665 nntilde_j = nntilde_jp1 ;
674 Scalar kkss ( psi_j_ma/(2. * nntilde_j )
680 double max_kss = kkss.val_grid_point(1, 0, 0, 0) ;
681 double min_kss = kkss.val_grid_point(1, 0, 0, 0) ;
684 double max_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
685 double min_hh_tilde = hh_tilde.val_grid_point(1, 0, 0, 0) ;
688 for (
int k=0 ; k<nnp ; k++)
689 for (
int j=0 ; j<nnt ; j++){
690 if (kkss.val_grid_point(1, k, j, 0) > max_kss)
691 max_kss = kkss.val_grid_point(1, k, j, 0) ;
692 if (kkss.val_grid_point(1, k, j, 0) < min_kss)
693 min_kss = kkss.val_grid_point(1, k, j, 0) ;
695 if (hh_tilde.val_grid_point(1, k, j, 0) > max_hh_tilde)
696 max_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
697 if (hh_tilde.val_grid_point(1, k, j, 0) < min_hh_tilde)
698 min_hh_tilde = hh_tilde.val_grid_point(1, k, j, 0) ;
700 kss_out << mer <<
" " << max_kss <<
" " << min_kss
701 <<
" " << -1 * max_hh_tilde <<
" " << -1 * min_hh_tilde << endl ;
706 if (solve_lapse == 1)
708 if (solve_shift == 1)
711 if (solve_shift == 1)
718 Scalar psi_j_a =
pow(psi_j, a) ;
719 psi_j_a.std_spectral_base() ;
724 if (solve_lapse == 1)
726 if (solve_shift == 1)
729 if (solve_shift == 1)
const Tensor_sym & get_delta() const
Returns the tensor which defines the connection with respect to the flat one: is the difference bet...
void set(const Map *mp, Mtbl *(*construct)(const Map *))
Semi-constructor from a mapping and a method.
double omega
Angular velocity in LORENE's units.
Scalar trK
Trace of the extrinsic curvature.
Metric met_gamt
3 metric tilde
void set_psi(const Scalar &psi_in)
Sets the conformal factor relating the physical metric to the conformal one: .
Scalar trK_point
Time derivative of the trace of the extrinsic curvature.
double kappa_hor() const
Surface gravity
void update_aa()
Conformal representation of the traceless part of the extrinsic curvature: .
double radius
Radius of the horizon in LORENE's units.
Map_af & mp
Affine mapping.
const Valeur boundary_nn_Neu_kk(int nn=1) const
Neumann boundary condition for N using the extrinsic curvature.
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Coord y
y coordinate centered on the grid
Coord ya
Absolute y coordinate.
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Coord x
x coordinate centered on the grid
Coord xa
Absolute x coordinate.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the metric is defined.
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
virtual const Vector & radial_vect() const
Returns the radial vector normal to a spherical slicing and pointing toward spatial infinity.
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
const Sym_tensor & ricci() const
Returns the Ricci tensor (given by the Connection p_connect )
virtual const Connection & connect() const
Returns the connection.
virtual const Scalar & ricci_scal() const
Returns the Ricci scalar.
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
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...
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
const Vector & divergence(const Metric &) const
Returns the divergence of this with respect to a Metric .
void inc_dzpuis()
dzpuis += 1 ;
void dec_dzpuis()
dzpuis -= 1 ;
const Tenseur & derive_con(const Metrique &) const
Returns the contravariant derivative of *this , with respect to met .
virtual const Sym_tensor & hh(Param *=0x0, Param *=0x0) const
Deviation of the conformal metric from the flat metric : .
virtual const Scalar & nn() const
Lapse function N at the current time step (jtime )
virtual const Vector & hdirac() const
Vector which vanishes in Dirac gauge.
virtual const Scalar & psi() const
Conformal factor relating the physical metric to the conformal one: .
const Metric_flat & ff
Pointer on the flat metric with respect to which the conformal decomposition is performed.
int jtime
Time step index of the latest slice.
virtual const Vector & beta() const
shift vector at the current time step (jtime )
Evolution_std< Scalar > n_evol
Values at successive time steps of the lapse function N.
Evolution_std< double > the_time
Time label of each slice.
Evolution_std< Vector > beta_evol
Values at successive time steps of the shift vector .
const Scalar & divergence(const Metric &) const
The divergence of this with respect to a Metric .
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Cmp sqrt(const Cmp &)
Square root.
Cmp log10(const Cmp &)
Basis 10 logarithm.
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Cmp pow(const Cmp &, int)
Power .
const Tensor_sym & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
const Tensor_sym & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Tbl maxabs(const Tensor &aa, const char *comment=0x0, ostream &ost=cout, bool verb=true)
Maxima in each domain of the absolute values of the tensor components.
Standard units of space, time and mass.