29 char et_magnetisation_comp_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Etoile/et_magnetisation_comp.C,v 1.13 2015/06/12 12:38:25 j_novak Exp $" ;
81 #include "et_rot_mag.h" 83 #include "utilitaires.h" 85 #include "proto_f77.h" 95 Cmp (*f_j)(
const Cmp&,
const double),
96 Param& par_poisson_At,
97 Param& par_poisson_Avect){
98 double relax_mag = 0.5 ;
100 int Z = mp.get_mg()->get_nzone();
102 bool adapt(adapt_flag) ;
106 int nt = mp.get_mg()->get_nt(nzet-1) ;
107 for (
int l=0; l<Z; l++) assert(mp.get_mg()->get_nt(l) == nt) ;
113 Mtbl* theta = mp.tet.c ;
115 assert (mpr != 0x0) ;
116 for (
int j=0; j<nt; j++)
117 Rsurf.
set(j) = mpr->
val_r_jk(l_surf()(0,j), xi_surf()(0,j), j, 0) ;
122 Cmp A_0t(- omega * A_phi) ;
128 BMN = BMN +
log(bbb) ;
141 Cmp ATANT(A_phi.srdsdt());
145 Cmp ttnphi(tnphi()) ;
147 Cmp BLAH(- b_car()/(nnn()*nnn())*ttnphi*grad1) ;
148 BLAH -= (1+b_car()/(nnn()*nnn())*tnphi()*tnphi())*grad2 ;
151 Cmp npgrada(2*nphisr*(A_phi.dsdr()+ATANT )) ;
153 BLAH -= grad3 + npgrada ;
154 Cmp gtt(-nnn()*nnn()+b_car()*tnphi()*tnphi()) ;
155 Cmp gtphi( - b_car()*ttnphi) ;
160 Cmp F01 = 1/(a_car()*nnn()*nnn())*A_0t.
dsdr()
161 + 1/(a_car()*nnn()*nnn())*nphi()*A_phi.
dsdr() ;
163 Cmp F02 = 1/(a_car()*nnn()*nnn())*A_0t.
srdsdt()
164 + 1/(a_car()*nnn()*nnn())*nphi()*A_phi.
srdsdt() ;
166 Cmp tmp = A_phi.
dsdr() / (bbb() * bbb() * a_car() );
169 Cmp F31 = 1/(a_car()*nnn()*nnn())*nphi()*nphi()*A_phi.
dsdr()
170 + 1/(a_car()*nnn()*nnn())*nphi()*A_0t.
dsdr()
173 tmp = A_phi.
srdsdt() / (bbb() * bbb() * a_car() );
176 Cmp F32 = 1/(a_car()*nnn()*nnn())*nphi()*nphi()*A_phi.
srdsdt()
177 + 1/(a_car()*nnn()*nnn())*nphi()*A_0t.
srdsdt()
180 Cmp x = get_magnetisation();
181 Cmp one_minus_x = 1 - x ;
184 tmp = ((BLAH - A_0t.
laplacien())*one_minus_x/a_car()
187 - gtphi*(F31*x.
dsdr()+F32*x.
srdsdt()) ) / gtt ;
195 for (
int j=0; j<nt; j++)
196 for (
int l=0; l<nzet; l++)
197 for (
int i=0; i<mp.get_mg()->get_nr(l); i++)
198 j_t.set(l,0,j,i) = ( (*mp.r.c)(l,0,j,i) > Rsurf(j) ?
199 0. : tmp(l,0,j,i) ) ;
200 j_t.annule(nzet,Z-1) ;
202 j_t.std_base_scal() ;
205 j_phi = omega * j_t + (ener() + press())*f_j(A_phi, a_j) ;
206 j_phi.std_base_scal() ;
214 Tenseur source_tAphi(mp, 1, CON, mp.get_bvect_spher()) ;
222 d_grad4.div_rsint() ;
223 source_tAphi.set(0)=0 ;
224 source_tAphi.set(1)=0 ;
227 Cmp phifac = (F31-nphi()*F01)*x.
dsdr()
228 + (F32-nphi()*F02)*x.
srdsdt() ;
230 source_tAphi.set(2)= -b_car()*a_car()/one_minus_x
231 *(tjphi-tnphi()*j_t + phifac)
232 + b_car()/(nnn()*nnn())*(tgrad1+tnphi()*grad2)
235 source_tAphi.change_triad(mp.get_bvect_cart());
238 for (
int i=0; i<3; i++) {
239 Scalar tmp_filter = source_tAphi(i) ;
242 source_tAphi.set(i) = tmp_filter ;
245 Tenseur WORK_VECT(mp, 1, CON, mp.get_bvect_cart()) ;
247 for (
int i=0; i<3; i++) {
248 WORK_VECT.set(i) = 0 ;
252 WORK_SCAL.
set() = 0 ;
254 double lambda_mag = 0. ;
257 if (source_tAphi.get_etat() != ETATZERO) {
259 for (
int i=0; i<3; i++) {
260 if(source_tAphi(i).dz_nonzero()) {
261 assert( source_tAphi(i).get_dzpuis() == 4 ) ;
264 (source_tAphi.set(i)).set_dzpuis(4) ;
270 source_tAphi.poisson_vect(lambda_mag, par_poisson_Avect, AVECT, WORK_VECT,
273 Cmp A_phi_n(AVECT(2));
278 Cmp source_A_1t(-a_car()*( j_t*gtt + j_phi*gtphi
280 + gtphi*(F31*x.
dsdr()+F32*x.
srdsdt()) )/one_minus_x
282 Scalar tmp_filter = source_A_1t ;
284 tmp_filter.exponential_filter_ylm(0, 2, 1) ;
285 source_A_1t = tmp_filter ;
289 source_A_1t.
poisson(par_poisson_At, A_1t) ;
291 int L = mp.get_mg()->get_nt(0);
309 for (
int p=0; p<mp.get_mg()->get_np(0); p++) {
312 for(
int k=0;k<L;k++){
313 for(
int l=0;l<2*L;l++){
315 if(l==0) leg.
set(k,l)=1. ;
316 if(l==1) leg.
set(k,l)=
cos((*theta)(l_surf()(p,k),p,k,0)) ;
317 if(l>=2) leg.
set(k,l) = double(2*l-1)/double(l)
318 *
cos((*theta)(l_surf()(p,k),p,k,0))
319 * leg(k,l-1)-double(l-1)/double(l)*leg(k,l-2) ;
323 for(
int k=0;k<L;k++){
330 for(
int l=0;l<L;l++) MAT.
set(l,k) = leg(k,2*l)/
pow(Rsurf(k),2*l+1);
335 int* IPIV=
new int[L] ;
343 F77_dgesv(&L, &un, MAT.
t, &L, IPIV, VEC.
t, &L, &INFO) ;
347 for(
int k=0;k<L;k++) {VEC2.
set(k)=1. ; }
349 F77_dgesv(&L, &un, MAT_SAVE.
t, &L, IPIV, VEC2.
t, &L, &INFO) ;
353 for(
int nz=0;nz < Z; nz++){
354 for(
int i=0;i< mp.get_mg()->get_nr(nz);i++){
355 for(
int k=0;k<L;k++){
356 psi.
set(nz,p,k,i) = 0. ;
357 psi2.
set(nz,p,k,i) = 0. ;
358 for(
int l=0;l<L;l++){
359 psi.
set(nz,p,k,i) += VEC(l)*leg(k,2*l) /
360 pow((*mp.r.c)(nz,p,k,i),2*l+1);
361 psi2.
set(nz,p,k,i) += VEC2(l)*leg(k,2*l)/
362 pow((*mp.r.c)(nz, p, k,i),2*l+1);
378 Cmp A_t_ext(A_1t + psi) ;
379 A_t_ext.
annule(0,nzet-1) ;
385 for (
int j=0; j<nt; j++)
386 for (
int l=0; l<Z; l++)
387 for (
int i=0; i<mp.get_mg()->get_nr(l); i++)
388 A_0t.
set(l,0,j,i) = ( (*mp.r.c)(l,0,j,i) > Rsurf(j) ?
389 A_1t(l,0,j,i) + psi(l,0,j,i) : tmp(l,0,j,i) ) ;
394 tmp_filter.exponential_filter_r(0, 2, 1) ;
395 tmp_filter.exponential_filter_ylm(0, 2, 1) ;
400 double Q_0 = -4*M_PI*(*asymp[1])(Z-1,0,0,0) ;
408 double Q_2 = -4*M_PI*(*asymp[1])(Z-1,0,0,0) ;
416 double C = (Q-Q_0)/Q_2 ;
426 Cmp A_t_ext(A_0t + C*psi2) ;
427 A_t_ext.
annule(0,nzet-1) ;
433 for (
int j=0; j<nt; j++)
434 for (
int l=0; l<Z; l++)
435 for (
int i=0; i<mp.get_mg()->get_nr(l); i++)
436 A_t_n.
set(l,0,j,i) = ( (*mp.r.c)(l,0,j,i) > Rsurf(j) ?
437 A_0t(l,0,j,i) + C*psi2(l,0,j,i) :
438 A_0t(l,0,j,i) + C ) ;
442 tmp_filter.exponential_filter_r(0, 2, 1) ;
443 tmp_filter.exponential_filter_ylm(0, 2, 1) ;
452 A_t = relax_mag*A_t_n + (1.-relax_mag)*A_t ;
453 A_phi = relax_mag*A_phi_n + (1. - relax_mag)*A_phi ;
472 if (ApAp.get_etat() != ETATZERO) {
473 ApAp.
set().div_rsint() ;
474 ApAp.set().div_rsint() ;
476 if (ApAt.get_etat() != ETATZERO)
477 ApAt.set().div_rsint() ;
479 E_em = 0.5*mu0 * ( 1/(a_car*nnn*nnn) * (AtAt + 2*tnphi*ApAt)
480 + ( (tnphi*tnphi/(a_car*nnn*nnn)) + 1/(a_car*b_car) )*ApAp );
481 Jp_em = -mu0 * (ApAt + tnphi*ApAp) /(a_car*nnn) ;
482 if (Jp_em.get_etat() != ETATZERO) Jp_em.set().mult_rsint() ;
494 Vector U_up(mp, CON, mp.get_bvect_cart()) ;
495 for (
int i=1; i<=3; i++)
496 U_up.
set(i) = u_euler(i-1) ;
497 U_up.change_triad(mp.get_bvect_spher()) ;
499 Sym_tensor gamij(mp, COV, mp.get_bvect_spher()) ;
500 for (
int i=1; i<=3; i++)
501 for (
int j=1; j<i; j++) {
504 gamij.
set(1,1) = a_car() ;
505 gamij.set(2,2) = a_car() ;
506 gamij.set(3,3) = b_car() ;
511 Vector B_up(mp, CON, mp.get_bvect_spher()) ;
517 fac =
Scalar(gam_euler()*gam_euler()) ;
519 E_I = get_magnetisation() * EiEi / mu0 ;
521 J_I = get_magnetisation() * BiBi * Ui / mu0 ;
522 Sij_I = get_magnetisation()
523 * ( (BiBi / fac) * gamij + BiBi*Ui*Ui - Bi*Bi / fac ) / mu0 ;
525 for (
int i=1; i<=3; i++)
526 for (
int j=i; j<=3; j++)
527 Sij_I.set(i,j).set_dzpuis(0) ;
537 if (p_mass_g == 0x0) {
542 Tenseur SrrplusStt(
Cmp(Sij_I(1, 1) + Sij_I(2, 2)) ) ;
543 SrrplusStt = SrrplusStt / a_car ;
552 Tenseur source = nnn * (ener_euler + E_em + E_i
553 + s_euler + Spp_em + SrrplusStt + Spp) +
555 + 2 * bbb * (ener_euler + press) * tnphi * uuu ;
557 source = a_car * bbb * source ;
561 p_mass_g =
new double( source().integrale() ) ;
566 p_mass_g =
new double( mass_b() ) ;
581 if (p_angu_mom == 0x0) {
586 dens.
va = (dens.
va).mult_st() ;
589 dens = a_car() * (b_car() * (ener_euler() + press())
590 * dens + bbb() * (Jp_em() +
Cmp(J_I(3)) ) ) ;
593 dens = nbar() * dens ;
598 p_angu_mom =
new double( dens.
integrale() ) ;
620 Tenseur sou_m = 2 * qpig * a_car * (press + (ener_euler+press)
623 Tenseur sou_q = 2 * qpig * a_car * Spp_em + 1.5 * ak_car
626 p_grv2 =
new double(
double(1) - lambda_grv2(sou_m(), sou_q()) ) ;
656 source = 0.75 * ak_car
658 logn.gradient_spher() )
662 Cmp aa = alpha() - 0.5 * beta() ;
668 cout <<
"Etoile_rot::grv3: the mapping does not belong" 669 <<
" to the class Map_radial !" << endl ;
676 vdaadt = vdaadt.
ssint() ;
681 temp = ( bbb() - a_car()/bbb() ) * temp ;
689 vtemp = (mpr->
xsr) * vtemp ;
697 source = bbb() * source() + 0.5 * temp ;
702 logn.gradient_spher() ) ;
707 double int_grav = source().integrale() ;
715 Tenseur SrrplusStt(
Cmp(Sij_I(1, 1) + Sij_I(2, 2)) ) ;
716 SrrplusStt = SrrplusStt / a_car ;
721 source = qpig * a_car * bbb * ( s_euler + Spp_em + SrrplusStt + Spp ) ;
724 source = qpig * ( 3 * press + nbar * uuu * uuu ) ;
729 double int_mat = source().integrale() ;
734 *ost <<
"Et_magnetisation::grv3 : gravitational term : " << int_grav
736 *ost <<
"Et_magnetisation::grv3 : matter term : " << int_mat
740 p_grv3 =
new double( (int_grav + int_mat) / int_mat ) ;
754 if (p_mom_quad_old == 0x0) {
766 Tenseur SrrplusStt(
Cmp(Sij_I(1, 1) + Sij_I(2, 2)) ) ;
767 SrrplusStt = SrrplusStt / a_car ;
777 source = qpig * a_car *( ener_euler + E_em + E_i
778 + s_euler + Spp_em + SrrplusStt + Spp)
783 source = qpig * nbar ;
793 Cmp& csource = source.
set() ;
805 assert( dynamic_cast<const Map_radial*>(&mp) != 0x0 ) ;
815 source = 0.5 * source() - 1.5 * temp ;
819 p_mom_quad_old =
new double( source().integrale() / qpig ) ;
821 return *p_mom_quad_old ;
829 if (p_mom_quad_Bo == 0x0) {
832 Tenseur SrrplusStt(
Cmp(Sij_I(1, 1) + Sij_I(2, 2)) ) ;
834 Cmp dens = a_car() * press() ;
835 dens = bbb() * nnn() * (SrrplusStt() + 2*dens) ;
839 p_mom_quad_Bo =
new double( - 16. * dens.
integrale() / qpig ) ;
841 return *p_mom_quad_Bo ;
const Cmp & dsdr() const
Returns of *this .
Metric for tensor calculation.
virtual double mom_quad_Bo() const
Part of the quadrupole moment.
Cmp log(const Cmp &)
Neperian logarithm.
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
const Tenseur & gradient_spher() const
Returns the gradient of *this (Spherical coordinates) (scalar field only).
Cmp sqrt(const Cmp &)
Square root.
void dec_dzpuis()
Decreases by 1 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
void annule(int l)
Sets the Cmp to zero in a given domain.
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Standard units of space, time and mass.
int get_etat() const
Returns the logical state.
double & set(int i)
Read/write of a particular element (index i) (1D case)
Tensor field of valence 0 (or component of a tensorial field).
Tenseur flat_scalar_prod(const Tenseur &t1, const Tenseur &t2)
Scalar product of two Tenseur when the metric is : performs the contraction of the last index of t1 w...
const Cmp & srdsdt() const
Returns of *this .
const Valeur & sx() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR ) ...
virtual void exponential_filter_r(int lzmin, int lzmax, int p, double alpha=-16.)
Applies an exponential filter to the spectral coefficients in the radial direction.
Values and coefficients of a (real-value) function.
Tensor field of valence 1.
Cmp cos(const Cmp &)
Cosine.
void div_r()
Division by r everywhere.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
virtual double angu_mom() const
Angular momentum.
void mult_r()
Multiplication by r everywhere.
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Tenseur flat_scalar_prod_desal(const Tenseur &t1, const Tenseur &t2)
Same as flat_scalar_prod but with desaliasing.
const Valeur & ssint() const
Returns of *this.
double * t
The array of double.
Base class for pure radial mappings.
void mult_rsint()
Multiplication by .
virtual double val_r_jk(int l, double xi, int j, int k) const =0
Returns the value of the radial coordinate r for a given and a given collocation point in in a give...
double integrale() const
Computes the integral over all space of *this .
double val_point_jk(int l, double x, int j, int k) const
Computes the value of the field represented by *this at an arbitrary point in , but collocation point...
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain...
Cmp pow(const Cmp &, int)
Power .
void inc2_dzpuis()
Increases by 2 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
virtual double mom_quad_old() const
Part of the quadrupole moment.
Cmp poisson() const
Solves the scalar Poisson equation with *this as a source.
virtual double grv2() const
Error on the virial identity GRV2.
void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Tbl & set(int l)
Read/write of the value in a given domain.
int get_dzpuis() const
Returns dzpuis.
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
virtual void magnet_comput(const int adapt_flag, Cmp(*f_j)(const Cmp &x, const double), Param &par_poisson_At, Param &par_poisson_Avect)
Computes the electromagnetic quantities solving the Maxwell equations (6) and (7) of [Bocquet...
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
virtual double mass_g() const
Gravitational mass.
Valeur ** asymptot(int n, const int flag=0) const
Asymptotic expansion at r = infinity.
const Valeur & mult_ct() const
Returns applied to *this.
virtual double grv3(ostream *ost=0x0) const
Error on the virial identity GRV3.
void set_dzpuis(int)
Set a value to dzpuis.
Scalar & set(int)
Read/write access to a component.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
virtual void exponential_filter_ylm(int lzmin, int lzmax, int p, double alpha=-16.)
Applies an exponential filter to the spectral coefficients in the angular directions.
const Cmp & laplacien(int zec_mult_r=4) const
Returns the Laplacian of *this.
Valeur va
The numerical value of the Cmp.
Class intended to describe valence-2 symmetric tensors.
Standard electro-magnetic units.
Tensor handling *** DEPRECATED : use class Tensor instead ***.
virtual void MHD_comput()
Computes the electromagnetic part of the stress-energy tensor.
Tensor down(int ind, const Metric &gam) const
Computes a new tensor by lowering an index of *this.