29 char et_bin_nsbh_equilibrium_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_nsbh_equilibrium.C,v 1.13 2014/10/13 08:52:56 j_novak Exp $" ;
89 #include "et_bin_nsbh.h"
92 #include "graphique.h"
93 #include "utilitaires.h"
98 int mermax_poisson,
double relax_poisson,
99 int mermax_potvit,
double relax_potvit,
113 Map_et& mp_et =
dynamic_cast<Map_et&
>(
mp) ;
117 double& diff_ent = diff.set(0) ;
118 double& diff_vel_pot = diff.set(1) ;
119 double& diff_lapse = diff.set(2) ;
120 double& diff_confpsi = diff.set(3) ;
121 double& diff_shift_x = diff.set(4) ;
122 double& diff_shift_y = diff.set(5) ;
123 double& diff_shift_z = diff.set(6) ;
128 double precis_poisson = 1.e-16 ;
131 par_poisson1.add_int(mermax_poisson, 0) ;
132 par_poisson1.add_double(relax_poisson, 0) ;
133 par_poisson1.add_double(precis_poisson, 1) ;
134 par_poisson1.add_int_mod(niter, 0) ;
141 par_poisson2.add_int(mermax_poisson, 0) ;
142 par_poisson2.add_double(relax_poisson, 0) ;
143 par_poisson2.add_double(precis_poisson, 1) ;
144 par_poisson2.add_int_mod(niter, 0) ;
150 Param par_poisson_vect ;
151 par_poisson_vect.add_int(mermax_poisson, 0) ;
153 par_poisson_vect.add_double(relax_poisson, 0) ;
154 par_poisson_vect.add_double(precis_poisson, 1) ;
155 par_poisson_vect.add_cmp_mod(
ssjm1_khi ) ;
157 par_poisson_vect.add_int_mod(niter, 0) ;
163 int adapt_flag = (adapt) ? 1 : 0 ;
164 int nz_search =
nzet + 1 ;
165 double precis_secant = 1.e-14 ;
167 double reg_map = 1. ;
170 Tbl ent_limit(
nzet) ;
172 par_adapt.add_int(nitermax, 0) ;
173 par_adapt.add_int(
nzet, 1) ;
174 par_adapt.add_int(nz_search, 2) ;
175 par_adapt.add_int(adapt_flag, 3) ;
176 par_adapt.add_int(j_b, 4) ;
177 par_adapt.add_int(k_b, 5) ;
178 par_adapt.add_int_mod(niter_adapt, 0) ;
179 par_adapt.add_double(precis_secant, 0) ;
180 par_adapt.add_double(reg_map, 1) ;
181 par_adapt.add_double(alpha_r, 2) ;
182 par_adapt.add_tbl(ent_limit, 0) ;
189 Tenseur ent_jm1 =
ent ;
199 for(
int mer=0 ; mer<mermax ; mer++ ) {
201 cout <<
"-----------------------------------------------" << endl ;
202 cout <<
"step: " << mer << endl ;
203 cout <<
"diff_ent = " << diff_ent << endl ;
218 int nt = mg->get_nt(
nzet-1) ;
219 int np = mg->get_np(
nzet-1) ;
220 int nr = mg->get_nr(
nzet-1) ;
223 double hc =
exp(ent_c) ;
224 double gamma_c =
exp(
loggam())(0,0,0,0) ;
226 double n_auto_c =
n_auto()(0,0,0,0) ;
227 double n_comp_c =
n_comp()(0,0,0,0) ;
229 double alpha_square = 0 ;
230 double constante = 0;
231 for (
int k=0; k<np; k++) {
232 for (
int j=0; j<nt; j++) {
241 double alpha_square_courant = (gamma_0_c*gamma_b*n_comp_b - hc*gamma_c*gamma_0_b*n_comp_c) /
242 (hc*gamma_c*gamma_0_b*n_auto_c-gamma_0_c*gamma_b*n_auto_b) ;
243 double constante_courant = gamma_b*(n_comp_b+alpha_square_courant*n_auto_b)/gamma_0_b ;
245 if (alpha_square_courant > alpha_square) {
246 alpha_square = alpha_square_courant ;
249 constante = constante_courant ;
254 alpha_r =
sqrt(alpha_square) ;
255 cout <<
"Adaptation : " << k_b <<
" " << j_b <<
" " << alpha_r << endl ;
259 potentiel.set_std_base() ;
260 for (
int l=
nzet+1 ; l<nz ; l++)
261 potentiel.set().va.set(l) = 1 ;
263 Map_et mp_prev = mp_et ;
269 for (
int l=0; l<
nzet; l++) {
270 ent_limit.set(l) =
ent()(l, k_b, j_b, nr-1) ;
275 mp_prev.homothetie(alpha_r) ;
277 for (
int l=
nzet ; l<nz-1 ; l++)
299 Tbl diff_ent_tbl =
diffrel(
ent(), ent_jm1() ) ;
300 diff_ent = diff_ent_tbl(0) ;
301 for (
int l=1; l<
nzet; l++) {
302 diff_ent += diff_ent_tbl(l) ;
325 tmp2.set_etat_qcq() ;
326 for (
int i=0 ; i<3 ; i++) {
327 tmp2.set() = tmp(i, i) ;
332 +
nnn * confpsi_q * kk
338 "WARNING : Et_bin_nsbh is for the relativistic calculation"
343 source.set_std_base() ;
347 Cmp n_auto_old (
n_auto()) ;
348 source().poisson(par_poisson1,
n_auto.
set()) ;
356 "Relative difference on n_auto : "
358 for (
int l=0; l<nz; l++) {
359 cout << tdiff_lapse(l) <<
" " ;
362 diff_lapse =
max(
abs(tdiff_lapse)) ;
379 tmp2.set_etat_qcq() ;
380 for (
int i=0 ; i<3 ; i++) {
381 tmp2.set() = tmp(i, i) ;
386 - 0.125 * confpsi_c * kk ;
402 "Relative difference on confpsi_auto : "
404 for (
int l=0; l<nz; l++) {
405 cout << tdiff_confpsi(l) <<
" " ;
408 diff_confpsi =
max(
abs(tdiff_confpsi)) ;
426 for (
int i=0 ; i<3 ; i++)
427 if (source_shift(i).get_etat() != ETATZERO)
428 source_shift.set(i).va.coef_i() ;
430 for (
int i=0; i<3; i++)
431 if ((source_shift(i).get_etat() != ETATZERO) && (source_shift(i).va.c->t[nz-1]->get_etat() != ETATZERO))
432 source_shift.set(i).filtre(4) ;
433 for (
int i=0; i<3; i++) {
434 if(source_shift(i).dz_nonzero()) {
435 assert( source_shift(i).get_dzpuis() == 4 ) ;
438 (source_shift.set(i)).set_dzpuis(4) ;
444 double lambda_shift = double(1) / double(3) ;
448 source_shift.poisson_vect_oohara(lambda_shift, par_poisson_vect,
462 "Relative difference on shift_auto : "
464 cout <<
"x component : " ;
465 for (
int l=0; l<nz; l++) {
466 cout << tdiff_shift_x(l) <<
" " ;
469 cout <<
"y component : " ;
470 for (
int l=0; l<nz; l++) {
471 cout << tdiff_shift_y(l) <<
" " ;
474 cout <<
"z component : " ;
475 for (
int l=0; l<nz; l++) {
476 cout << tdiff_shift_z(l) <<
" " ;
480 diff_shift_x =
max(
abs(tdiff_shift_x)) ;
481 diff_shift_y =
max(
abs(tdiff_shift_y)) ;
482 diff_shift_z =
max(
abs(tdiff_shift_z)) ;
495 int,
double,
double,
const Tbl&,
Tbl&) {
497 cout <<
"Not implemented !" << endl ;