30 char et_rot_mag_mag_plus_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_mag_mag_plus.C,v 1.3 2014/10/13 08:52:58 j_novak Exp $" ;
98 #include "et_rot_mag.h"
99 #include "utilitaires.h"
101 #include "proto_f77.h"
102 #include "graphique.h"
118 Param& par_poisson_At,
119 Param& par_poisson_Avect){
120 double relax_mag = 1.0 ;
122 int Z = mp.get_mg()->get_nzone();
125 bool adapt(adapt_flag) ;
129 int nt = mp.get_mg()->get_nt(nzet-1) ;
130 for (
int l=0; l<Z; l++) assert(mp.get_mg()->get_nt(l) == nt) ;
136 Mtbl* theta = mp.tet.c ;
138 assert (mpr != 0x0) ;
139 for (
int j=0; j<nt; j++)
140 Rsurf.
set(j) = mpr->
val_r_jk(l_surf()(0,j), xi_surf()(0,j), j, 0) ;
145 Cmp A_0t(- omega * A_phi) ;
151 BMN = BMN +
log(bbb) ;
164 Cmp ATANT(A_phi.srdsdt());
168 Cmp ttnphi(tnphi()) ;
170 Cmp BLAH(- b_car()/(nnn()*nnn())*ttnphi*grad1) ;
171 BLAH -= (1+b_car()/(nnn()*nnn())*tnphi()*tnphi())*grad2 ;
174 Cmp npgrada(2*nphisr*(A_phi.dsdr()+ATANT )) ;
176 BLAH -= grad3 + npgrada ;
177 Cmp gtt(-nnn()*nnn()+b_car()*tnphi()*tnphi()) ;
178 Cmp gtphi( - b_car()*ttnphi) ;
181 Cmp tmp(((BLAH - A_0t.
laplacien())/a_car() - gtphi*j_phi)
189 for (
int j=0; j<nt; j++)
190 for (
int l=0; l<nzet; l++)
191 for (
int i=0; i<mp.get_mg()->get_nr(l); i++)
192 j_t.set(l,0,j,i) = ( (*mp.r.c)(l,0,j,i) > Rsurf(j) ?
193 0. : tmp(l,0,j,i) ) ;
194 j_t.annule(nzet,Z-1) ;
196 j_t.std_base_scal() ;
202 if (maxA_phi.
set(0) == 0) {
204 cout <<
"Initializing j_phi" << endl;
208 for (
int i=0; i<nd - 4;i++){
209 aini = aini + fabs(an_j(i));
211 aini = aini + fabs (an_j(3)) + fabs (an_j(5));
214 for (
int i=0; i<nd;i++){
215 bini = bini + fabs(bn_j(i));
220 j_phi = (ener() + press())*aini + bini ;
221 j_phi.std_base_scal() ;
222 j_phi.annule(nzet,Z-1) ;
223 j_phi.std_base_scal() ;
227 j_phi = (ener() + press())*aini + bini;
228 j_phi.std_base_scal() ;
230 j_phi.annule(nzet,Z-1) ;
231 j_phi.std_base_scal() ;
235 j_phi = (ener() + press())*aini + bini ;
236 j_phi.std_base_scal() ;
241 j_phi.annule(nzet,Z-1) ;
242 j_phi.std_base_scal() ;
246 cout <<
"ERROR" << endl;
247 cout <<
"initial_j = " << initial_j << endl;
256 double maxA_phi_surf = 0.;
257 for (
int j = 0; j < nt ; j++ ) {
259 double maxA_phi_surf_tmp = fabs(A_phi(nzet,0,j,0));
260 maxA_phi_surf=
max(maxA_phi_surf, maxA_phi_surf_tmp);
263 Cmp A_phi_scaled = A_phi / maxA_phi.
set(0);
264 Cmp A_phi_scaled2 = (A_phi - maxA_phi_surf)/ maxA_phi.
set(0);
267 for (
int j=0; j<nt; j++) {
268 for (
int l=0; l<Z; l++) {
269 for (
int i=0; i<mp.get_mg()->get_nr(l); i++) {
270 if (A_phi_scaled2(l,0,j,i) < 0.) {
271 A_phi_scaled2.
set(l,0,j,i) = 0 ;
279 Cmp j_phi_ff = g_j (A_phi_scaled2, bn_j)/maxA_phi.
set(0);
282 j_phi_ff.div_rsint();
283 j_phi_ff.div_rsint();
284 j_phi_ff = j_phi_ff / b_car() / nnn() / nnn();
287 + (ener() + press())*f_j(A_phi_scaled, an_j)/maxA_phi.
set(0)/g_si
291 j_phi.std_base_scal() ;
295 B_phi = N_j (A_phi_scaled2, bn_j);
296 B_phi.std_base_scal() ;
297 B_phi = B_phi / nnn();
311 Tenseur source_tAphi(mp, 1, CON, mp.get_bvect_spher()) ;
319 d_grad4.div_rsint() ;
320 source_tAphi.
set(0)=0 ;
321 source_tAphi.
set(1)=0 ;
322 source_tAphi.
set(2)= -b_car()*a_car()*(tjphi-tnphi()*j_t)
323 + b_car()/(nnn()*nnn())*(tgrad1+tnphi()*grad2)+d_grad4 ;
327 Tenseur WORK_VECT(mp, 1, CON, mp.get_bvect_cart()) ;
329 for (
int i=0; i<3; i++) {
330 WORK_VECT.
set(i) = 0 ;
334 WORK_SCAL.
set() = 0 ;
336 double lambda_mag = 0. ;
339 if (source_tAphi.
get_etat() != ETATZERO) {
341 for (
int i=0; i<3; i++) {
342 if(source_tAphi(i).dz_nonzero()) {
343 assert( source_tAphi(i).get_dzpuis() == 4 ) ;
346 (source_tAphi.
set(i)).set_dzpuis(4) ;
351 source_tAphi.
poisson_vect(lambda_mag, par_poisson_Avect, AVECT, WORK_VECT,
354 Cmp A_phi_n(AVECT(2));
359 Cmp source_A_1t(-a_car()*(j_t*gtt + j_phi*gtphi) + BLAH);
363 source_A_1t.
poisson(par_poisson_At, A_1t) ;
365 int L = mp.get_mg()->get_nt(0);
383 for (
int p=0; p<mp.get_mg()->get_np(0); p++) {
386 for(
int k=0;k<L;k++){
387 for(
int l=0;l<2*L;l++){
389 if(l==0) leg.
set(k,l)=1. ;
390 if(l==1) leg.
set(k,l)=
cos((*theta)(l_surf()(p,k),p,k,0)) ;
391 if(l>=2) leg.
set(k,l) = double(2*l-1)/double(l)
392 *
cos((*theta)(l_surf()(p,k),p,k,0))
393 * leg(k,l-1)-double(l-1)/double(l)*leg(k,l-2) ;
397 for(
int k=0;k<L;k++){
404 for(
int l=0;l<L;l++) MAT.
set(l,k) = leg(k,2*l)/
pow(Rsurf(k),2*l+1);
409 int* IPIV=
new int[L] ;
417 F77_dgesv(&L, &un, MAT.
t, &L, IPIV, VEC.
t, &L, &INFO) ;
421 for(
int k=0;k<L;k++) {VEC2.
set(k)=1. ; }
423 F77_dgesv(&L, &un, MAT_SAVE.
t, &L, IPIV, VEC2.
t, &L, &INFO) ;
427 for(
int nz=0;nz < Z; nz++){
428 for(
int i=0;i< mp.get_mg()->get_nr(nz);i++){
429 for(
int k=0;k<L;k++){
430 psi.
set(nz,p,k,i) = 0. ;
431 psi2.
set(nz,p,k,i) = 0. ;
432 for(
int l=0;l<L;l++){
433 psi.
set(nz,p,k,i) += VEC(l)*leg(k,2*l) /
434 pow((*mp.r.c)(nz,p,k,i),2*l+1);
435 psi2.
set(nz,p,k,i) += VEC2(l)*leg(k,2*l)/
436 pow((*mp.r.c)(nz, p, k,i),2*l+1);
452 Cmp A_t_ext(A_1t + psi) ;
453 A_t_ext.
annule(0,nzet-1) ;
459 for (
int j=0; j<nt; j++)
460 for (
int l=0; l<Z; l++)
461 for (
int i=0; i<mp.get_mg()->get_nr(l); i++)
462 A_0t.
set(l,0,j,i) = ( (*mp.r.c)(l,0,j,i) > Rsurf(j) ?
463 A_1t(l,0,j,i) + psi(l,0,j,i) : tmp(l,0,j,i) ) ;
469 double Q_0 = -4*M_PI*(*asymp[1])(Z-1,0,0,0) ;
477 double Q_2 = -4*M_PI*(*asymp[1])(Z-1,0,0,0) ;
485 double C = (Q-Q_0)/Q_2 ;
495 Cmp A_t_ext(A_0t + C*psi2) ;
496 A_t_ext.
annule(0,nzet-1) ;
502 for (
int j=0; j<nt; j++)
503 for (
int l=0; l<Z; l++)
504 for (
int i=0; i<mp.get_mg()->get_nr(l); i++)
505 A_t_n.
set(l,0,j,i) = ( (*mp.r.c)(l,0,j,i) > Rsurf(j) ?
506 A_0t(l,0,j,i) + C*psi2(l,0,j,i) :
507 A_0t(l,0,j,i) + C ) ;
517 A_t = relax_mag*A_t_n + (1.-relax_mag)*A_t ;
518 A_phi = relax_mag*A_phi_n + (1. - relax_mag)*A_phi ;
522 cout <<
"not implemented" << endl;