29 char isol_hor_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Isol_hor/isol_hor.C,v 1.35 2014/10/13 08:53:01 j_novak Exp $" ;
161 #include "utilitaires.h"
162 #include "time_slice.h"
163 #include "isol_hor.h"
166 #include "evolution.h"
178 mp(mpi), nz(mpi.get_mg()->get_nzone()), radius ((mpi.get_alpha())[0]),
179 omega(0), boost_x(0), boost_z(0), regul(0),
180 n_auto_evol(depth_in), n_comp_evol(depth_in),
181 psi_auto_evol(depth_in), psi_comp_evol(depth_in),
182 dn_evol(depth_in), dpsi_evol(depth_in),
183 beta_auto_evol(depth_in), beta_comp_evol(depth_in),
184 aa_auto_evol(depth_in), aa_comp_evol(depth_in),
185 aa_nn(depth_in), aa_quad_evol(depth_in),
186 met_gamt(mpi.flat_met_spher()), gamt_point(mpi, CON, mpi.get_bvect_spher()),
187 trK(mpi), trK_point(mpi), decouple(mpi){
200 ff_in.con(), psi_in*psi_in*psi_in*psi_in*psi_in*psi_in*aa_in,
202 mp(mpi), nz(mpi.get_mg()->get_nzone()), radius ((mpi.get_alpha())[0]),
203 omega(0), boost_x(0), boost_z(0), regul(0),
204 n_auto_evol(depth_in), n_comp_evol(depth_in),
205 psi_auto_evol(depth_in), psi_comp_evol(depth_in),
206 dn_evol(depth_in), dpsi_evol(depth_in),
207 beta_auto_evol(depth_in), beta_comp_evol(depth_in),
208 aa_auto_evol(depth_in), aa_comp_evol(depth_in),
209 aa_nn(depth_in), aa_quad_evol(depth_in),
210 met_gamt(metgamt), gamt_point(gamt_point_in),
211 trK(trK_in), trK_point(trK_point_in), decouple(lapse_in.get_mp()){
226 radius(isolhor_in.radius),
227 omega(isolhor_in.omega),
228 boost_x(isolhor_in.boost_x),
229 boost_z(isolhor_in.boost_z),
230 regul(isolhor_in.regul),
231 n_auto_evol(isolhor_in.n_auto_evol),
232 n_comp_evol(isolhor_in.n_comp_evol),
233 psi_auto_evol(isolhor_in.psi_auto_evol),
234 psi_comp_evol(isolhor_in.psi_comp_evol),
235 dn_evol(isolhor_in.dn_evol),
236 dpsi_evol(isolhor_in.dpsi_evol),
237 beta_auto_evol(isolhor_in.beta_auto_evol),
238 beta_comp_evol(isolhor_in.beta_comp_evol),
239 aa_auto_evol(isolhor_in.aa_auto_evol),
240 aa_comp_evol(isolhor_in.aa_comp_evol),
241 aa_nn(isolhor_in.aa_nn),
242 aa_quad_evol(isolhor_in.aa_quad_evol),
243 met_gamt(isolhor_in.met_gamt),
244 gamt_point(isolhor_in.gamt_point),
246 trK_point(isolhor_in.trK_point),
247 decouple(isolhor_in.decouple){
254 bool partial_read,
int depth_in)
256 fich, partial_read, depth_in),
257 mp(mpi), nz(mpi.get_mg()->get_nzone()), radius ((mpi.get_alpha())[0]),
258 omega(0), boost_x(0), boost_z(0), regul(0),
259 n_auto_evol(depth_in), n_comp_evol(depth_in),
260 psi_auto_evol(depth_in), psi_comp_evol(depth_in),
261 dn_evol(depth_in), dpsi_evol(depth_in),
262 beta_auto_evol(depth_in), beta_comp_evol(depth_in),
263 aa_auto_evol(depth_in), aa_comp_evol(depth_in),
264 aa_nn(depth_in), aa_quad_evol(depth_in),
265 met_gamt(mpi.flat_met_spher()),
266 gamt_point(mpi, CON, mpi.get_bvect_spher()),
267 trK(mpi), trK_point(mpi), decouple(mpi){
277 for (
int j=jmin; j<=
jtime; j++) {
278 fread_be(&indicator,
sizeof(
int), 1, fich) ;
279 if (indicator == 1) {
286 for (
int j=jmin; j<=
jtime; j++) {
287 fread_be(&indicator,
sizeof(
int), 1, fich) ;
288 if (indicator == 1) {
295 for (
int j=jmin; j<=
jtime; j++) {
296 fread_be(&indicator,
sizeof(
int), 1, fich) ;
297 if (indicator == 1) {
304 for (
int j=jmin; j<=
jtime; j++) {
305 fread_be(&indicator,
sizeof(
int), 1, fich) ;
306 if (indicator == 1) {
382 flux <<
'\n' <<
"radius of the horizon : " <<
radius <<
'\n' ;
383 flux <<
"boost in x-direction : " <<
boost_x <<
'\n' ;
384 flux <<
"boost in z-direction : " <<
boost_z <<
'\n' ;
385 flux <<
"angular velocity omega : " <<
omega_hor() <<
'\n' ;
386 flux <<
"area of the horizon : " <<
area_hor() <<
'\n' ;
387 flux <<
"ang. mom. of horizon : " <<
ang_mom_hor() <<
'\n' ;
388 flux <<
"ADM ang. mom. : " <<
ang_mom_adm() <<
'\n' ;
389 flux <<
"Mass of the horizon : " <<
mass_hor() <<
'\n' ;
390 flux <<
"ADM Mass : " <<
adm_mass() <<
'\n' ;
419 for (
int j=jmin; j<=
jtime; j++) {
420 int indicator = (
psi_evol.is_known(j)) ? 1 : 0 ;
421 fwrite_be(&indicator,
sizeof(
int), 1, fich) ;
422 if (indicator == 1)
psi_evol[j].sauve(fich) ;
426 for (
int j=jmin; j<=
jtime; j++) {
427 int indicator = (
n_auto_evol.is_known(j)) ? 1 : 0 ;
428 fwrite_be(&indicator,
sizeof(
int), 1, fich) ;
433 for (
int j=jmin; j<=
jtime; j++) {
435 fwrite_be(&indicator,
sizeof(
int), 1, fich) ;
440 for (
int j=jmin; j<=
jtime; j++) {
442 fwrite_be(&indicator,
sizeof(
int), 1, fich) ;
592 for (
int i=1 ; i<=3 ; i++){
593 if (auxi(i).get_etat() != ETATZERO)
600 for (
int i=1 ; i<=3 ; i++){
684 for (
int i=1 ; i<=3 ; i++){
685 if (auxi(i).get_etat() != ETATZERO)
692 for (
int i=1 ; i<=3 ; i++){
740 auxi = 0.5 - 0.5/
mp.
r ;
773 temp_vect1.
set(2) = 0. ;
774 temp_vect1.
set(3) = 0. ;
850 for (
int k=0; k<nnp; k++)
851 for (
int j=0; j<nnt; j++){
852 if (
nn().val_grid_point(1, k, j , 0) < 1e-12){
863 cout <<
"regul = " <<
regul << endl ;
868 for (
int i=1 ; i<=3 ; i++)
869 for (
int j=i ; j<=3 ; j++) {
870 auxi = aa_new(i, j) ;
871 auxi = division_xpun (auxi, 0) ;
872 aa_new.
set(i,j) = auxi / nn_sxpun / 2. ;
903 Scalar psi_perturb (
pow(gamm(3,3), 0.25)) ;
918 void Isol_hor::aa_kerr_ww(
double mm,
double aaa) {
928 Scalar rbl = rr + mm + (mm*mm - aaa*aaa) / (4*rr) ;
931 Scalar sigma_inv = 1. / (rbl*rbl + aaa*aaa*cost*cost) ;
936 ww_pert = - 1*(mm*aaa*aaa*aaa*
pow(sint, 4.)*cost) * sigma_inv ;
937 ww_pert.set_spectral_va().set_base_r(0,
R_CHEBPIM_P) ;
938 for (
int l=1; l<
nz-1; l++)
939 ww_pert.set_spectral_va().set_base_r(l,
R_CHEB) ;
940 ww_pert.set_spectral_va().set_base_r(
nz-1,
R_CHEBU) ;
941 ww_pert.set_spectral_va().set_base_t(
T_COSSIN_CI) ;
942 ww_pert.set_spectral_va().set_base_p(
P_COSSIN) ;
951 dw_by.set(2) = 3 * aaa*mm*sint*sint*sint / rr ;
953 dw_by.set(2).set_spectral_va().set_base_r(0,
R_CHEBPIM_P) ;
954 for (
int l=1; l<
nz-1; l++)
955 dw_by.set(2).set_spectral_va().set_base_r(l,
R_CHEB) ;
956 dw_by.set(2).set_spectral_va().set_base_r(
nz-1,
R_CHEBU) ;
957 dw_by.set(2).set_spectral_va().set_base_t(
T_COSSIN_SI) ;
958 dw_by.set(2).set_spectral_va().set_base_p(
P_COSSIN) ;
968 Vector dw_pert(ww_pert.derive_con(
ff)) ;
987 Scalar aquad = aquad_1 + aquad_2 + aquad_3 ;
1012 Scalar s_r (ww_pert.derive_cov(
ff)(2)) ;
1023 Scalar s_t (ww_pert.derive_cov(
ff)(1)) ;
1029 temp.div_rsint_dzpuis(2) ;
1045 aij.set(3,1) = s_r ;
1046 aij.set(3,1).div_rsint() ;
1047 aij.set(3,2) = s_t ;
1048 aij.set(3,2).div_rsint() ;
1050 Base_val base_31 (aij(3,1).get_spectral_va().get_base()) ;
1051 Base_val base_32 (aij(3,2).get_spectral_va().get_base()) ;
1053 aij.set(3,1) = aij(3,1) *
pow(
gam_dd()(1,1), 5./3.)
1055 aij.set(3,1) = aij(3,1) *
pow(
psi(), -6.) ;
1056 aij.set(3,1).set_spectral_va().set_base(base_31) ;
1057 aij.set(3,2) = aij(3,2) *
pow(
gam_dd()(1,1), 5./3.)
1059 aij.set(3,2) = aij(3,2) *
pow(
psi(), -6.) ;
1060 aij.set(3,2).set_spectral_va().set_base(base_32) ;
1080 kij = kij *
pow(
gam().determinant(), -1./3.) ;
1081 kij.std_spectral_base() ;
1116 double fonc_expansion(
double rr,
const Param& par_expansion) {
1125 void Isol_hor::adapt_hor(
double c_min,
double c_max) {
1128 Scalar app_hor(
mp) ;
1129 app_hor.annule_hard() ;
1133 double precis = 1.e-13 ;
1138 for (
int nt=0; nt<
mp.
get_mg()->get_nt(1); nt++)
1139 for (
int np=0; np<
mp.
get_mg()->get_np(1); np++) {
1144 Param par_expansion ;
1146 par_expansion.add_double(theta, 0) ;
1147 par_expansion.add_double(phi, 1) ;
1148 double r_app_hor =
zerosec_b(fonc_expansion, par_expansion, c_min, c_max,
1149 precis, nitmax, nit) ;
1151 app_hor.set_grid_point(1, np, nt, 0) = r_app_hor ;
1160 Scalar trans_11 (
mp) ;
1162 r_new.annule_hard() ;
1164 for (
int l=1; l<
nz; l++)
1165 for (
int nr=0; nr<
mp.
get_mg()->get_nr(1); nr++)
1166 for (
int nt=0; nt<
mp.
get_mg()->get_nt(1); nt++)
1167 for (
int np=0; np<
mp.
get_mg()->get_np(1); np++) {
1168 r_new.set_grid_point(l, np, nt, nr) = rr.val_grid_point(l, np, nt, nr) -
1169 app_hor.val_grid_point(1, np, nt, 0) + 1 ;
1173 r_new.std_spectral_base() ;
1179 Scalar trans_12 (r_new.dsdt()) ;
1181 Scalar trans_13 (r_new.stdsdp()) ;
1183 for (
int nr=0; nr<
mp.
get_mg()->get_nr(1); nr++)
1184 for (
int nt=0; nt<
mp.
get_mg()->get_nt(1); nt++)
1185 for (
int np=0; np<
mp.
get_mg()->get_np(1); np++) {
1186 trans_12.set_grid_point(
nz-1, np, nt, nr) = trans_12.val_grid_point(1, np, nt, nr) ;
1187 trans_13.set_grid_point(
nz-1, np, nt, nr) = trans_13.val_grid_point(1, np, nt, nr) ;
1192 trans.set(1,1) = 1 ;
1195 trans.set(2,2) = 1. ;
1196 trans.set(3,3) = 1. ;
1197 trans.set(2,1) = 0. ;
1198 trans.set(3,1) = 0. ;
1199 trans.set(3,2) = 0. ;
1200 trans.set(2,3) = 0. ;
1201 trans.std_spectral_base() ;
1203 cout <<
"trans(1,3)" << endl <<
norme(trans(1,3)) << endl ;
1205 Sym_tensor gamma_uu (
gam().con()) ;
1206 Sym_tensor kk_uu (
k_uu()) ;
1208 gamma_uu =
contract(gamma_uu, 0, 1, trans * trans, 1, 3) ;
1209 kk_uu =
contract(kk_uu, 0, 1, trans * trans, 1, 3) ;
1211 Sym_tensor copie_gamma (gamma_uu) ;
1212 Sym_tensor copie_kk (kk_uu) ;
1215 kk_uu.dec_dzpuis(2) ;
1216 for(
int i=1; i<=3; i++)
1217 for(
int j=i; j<=3; j++){
1218 kk_uu.set(i,j).annule_hard() ;
1219 gamma_uu.set(i,j).annule_hard() ;
1222 copie_kk.dec_dzpuis(2) ;
1224 Scalar expa_trans(
mp) ;
1225 expa_trans.annule_hard() ;
1240 for(
int i=1; i<=3; i++)
1241 for(
int j=i; j<=3; j++)
1242 for (
int l=1; l<
nz; l++)
1243 for (
int nr=0; nr<
mp.
get_mg()->get_nr(1); nr++)
1244 for (
int nt=0; nt<
mp.
get_mg()->get_nt(1); nt++)
1245 for (
int np=0; np<
mp.
get_mg()->get_np(1); np++) {
1249 double r_inv = rr.val_grid_point(l, np, nt, nr) +
1250 app_hor.val_grid_point(1, np, nt, 0) - 1. ;
1252 gamma_uu.set(i,j).set_grid_point(l, np, nt, nr) =
1253 copie_gamma(i,j).val_point(r_inv, theta, phi) ;
1254 kk_uu.set(i,j).set_grid_point(l, np, nt, nr) =
1255 copie_kk(i,j).val_point(r_inv, theta, phi) ;
1257 expa_trans.set_grid_point(l, np, nt, nr) = expa.
val_point(r_inv, theta, phi) ;
1259 kk_uu.std_spectral_base() ;
1260 gamma_uu.std_spectral_base() ;
1261 expa_trans.std_spectral_base() ;
1263 for (
int l=1; l<
nz; l++)
1264 for (
int nr=0; nr<
mp.
get_mg()->get_nr(1); nr++)
1265 for (
int nt=0; nt<
mp.
get_mg()->get_nt(1); nt++)
1266 for (
int np=0; np<
mp.
get_mg()->get_np(1); np++) {
1267 gamma_uu.set(1,2).set_grid_point(l,np,nt,nr) = gamma_uu.set(1,2).val_grid_point(l,np,nt,nr)
1268 / (1 +app_hor.val_grid_point(1, np, nt, 0)/ rr.val_grid_point(l, np, nt, nr)
1269 - 1 / rr.val_grid_point(l, np, nt, nr) ) ;
1270 gamma_uu.set(1,3).set_grid_point(l,np,nt,nr) = gamma_uu.set(1,3).val_grid_point(l,np,nt,nr)
1271 / (1 +app_hor.val_grid_point(1, np, nt, 0)/ rr.val_grid_point(l, np, nt, nr)
1272 - 1 / rr.val_grid_point(l, np, nt, nr) ) ;
1273 gamma_uu.set(2,2).set_grid_point(l,np,nt,nr) = gamma_uu.set(2,2).val_grid_point(l,np,nt,nr)
1274 / (1 +app_hor.val_grid_point(1, np, nt, 0)/ rr.val_grid_point(l, np, nt, nr)
1275 - 1 / rr.val_grid_point(l, np, nt, nr) )
1276 / (1 +app_hor.val_grid_point(1, np, nt, 0)/ rr.val_grid_point(l, np, nt, nr)
1277 - 1 / rr.val_grid_point(l, np, nt, nr) ) ;
1278 gamma_uu.set(2,3).set_grid_point(l,np,nt,nr) = gamma_uu.set(2,3).val_grid_point(l,np,nt,nr)
1279 / (1 +app_hor.val_grid_point(1, np, nt, 0)/ rr.val_grid_point(l, np, nt, nr)
1280 - 1 / rr.val_grid_point(l, np, nt, nr) )
1281 / (1 +app_hor.val_grid_point(1, np, nt, 0)/ rr.val_grid_point(l, np, nt, nr)
1282 - 1 / rr.val_grid_point(l, np, nt, nr) ) ;
1283 gamma_uu.set(3,3).set_grid_point(l,np,nt,nr) = gamma_uu.set(3,3).val_grid_point(l,np,nt,nr)
1284 / (1 +app_hor.val_grid_point(1, np, nt, 0)/ rr.val_grid_point(l, np, nt, nr)
1285 - 1 / rr.val_grid_point(l, np, nt, nr) )
1286 / (1 +app_hor.val_grid_point(1, np, nt, 0)/ rr.val_grid_point(l, np, nt, nr)
1287 - 1 / rr.val_grid_point(l, np, nt, nr) ) ;
1303 cout <<
"gamma_uu(1,1)" << endl <<
norme(gamma_uu(1,1)) << endl ;
1304 cout <<
"gamma_uu(2,1)" << endl <<
norme(gamma_uu(2,1)) << endl ;
1305 cout <<
"gamma_uu(3,1)" << endl <<
norme(gamma_uu(3,1)) << endl ;
1306 cout <<
"gamma_uu(2,2)" << endl <<
norme(gamma_uu(2,2)) << endl ;
1307 cout <<
"gamma_uu(3,2)" << endl <<
norme(gamma_uu(3,2)) << endl ;
1308 cout <<
"gamma_uu(3,3)" << endl <<
norme(gamma_uu(3,3)) << endl ;
1310 kk_uu.inc_dzpuis(2) ;
1311 expa_trans.inc_dzpuis(2) ;
1313 Metric gamm (gamma_uu) ;