29 char et_magnetisation_comp_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Etoile/et_magnetisation_comp.C,v 1.14 2016/11/01 09:12:59 j_novak Exp $" ;
84 #include "et_rot_mag.h"
86 #include "utilitaires.h"
88 #include "proto_f77.h"
98 Cmp (*f_j)(
const Cmp&,
const double),
99 Param& par_poisson_At,
100 Param& par_poisson_Avect){
101 double relax_mag = 0.5 ;
103 int Z = mp.get_mg()->get_nzone();
105 bool adapt(adapt_flag) ;
109 int nt = mp.get_mg()->get_nt(nzet-1) ;
110 for (
int l=0; l<Z; l++) assert(mp.get_mg()->get_nt(l) == nt) ;
116 Mtbl* theta = mp.tet.c ;
118 assert (mpr != 0x0) ;
119 for (
int j=0; j<nt; j++)
120 Rsurf.
set(j) = mpr->
val_r_jk(l_surf()(0,j), xi_surf()(0,j), j, 0) ;
125 Cmp A_0t(- omega * A_phi) ;
131 BMN = BMN +
log(bbb) ;
144 Cmp ATANT(A_phi.srdsdt());
148 Cmp ttnphi(tnphi()) ;
150 Cmp BLAH(- b_car()/(nnn()*nnn())*ttnphi*grad1) ;
151 BLAH -= (1+b_car()/(nnn()*nnn())*tnphi()*tnphi())*grad2 ;
154 Cmp npgrada(2*nphisr*(A_phi.dsdr()+ATANT )) ;
156 BLAH -= grad3 + npgrada ;
157 Cmp gtt(-nnn()*nnn()+b_car()*tnphi()*tnphi()) ;
158 Cmp gtphi( - b_car()*ttnphi) ;
163 Cmp F01 = 1/(a_car()*nnn()*nnn())*A_0t.
dsdr()
164 + 1/(a_car()*nnn()*nnn())*nphi()*A_phi.
dsdr() ;
166 Cmp F02 = 1/(a_car()*nnn()*nnn())*A_0t.
srdsdt()
167 + 1/(a_car()*nnn()*nnn())*nphi()*A_phi.
srdsdt() ;
169 Cmp tmp = A_phi.
dsdr() / (bbb() * bbb() * a_car() );
172 Cmp F31 = 1/(a_car()*nnn()*nnn())*nphi()*nphi()*A_phi.
dsdr()
173 + 1/(a_car()*nnn()*nnn())*nphi()*A_0t.
dsdr()
176 tmp = A_phi.
srdsdt() / (bbb() * bbb() * a_car() );
179 Cmp F32 = 1/(a_car()*nnn()*nnn())*nphi()*nphi()*A_phi.
srdsdt()
180 + 1/(a_car()*nnn()*nnn())*nphi()*A_0t.
srdsdt()
183 Cmp x = get_magnetisation();
184 Cmp one_minus_x = 1 - x ;
187 tmp = ((BLAH - A_0t.
laplacien())*one_minus_x/a_car()
190 - gtphi*(F31*x.
dsdr()+F32*x.
srdsdt()) ) / gtt ;
198 for (
int j=0; j<nt; j++)
199 for (
int l=0; l<nzet; l++)
200 for (
int i=0; i<mp.get_mg()->get_nr(l); i++)
201 j_t.set(l,0,j,i) = ( (*mp.r.c)(l,0,j,i) > Rsurf(j) ?
202 0. : tmp(l,0,j,i) ) ;
203 j_t.annule(nzet,Z-1) ;
205 j_t.std_base_scal() ;
208 j_phi = omega * j_t + (ener() + press())*f_j(A_phi, a_j) ;
209 j_phi.std_base_scal() ;
217 Tenseur source_tAphi(mp, 1, CON, mp.get_bvect_spher()) ;
225 d_grad4.div_rsint() ;
226 source_tAphi.
set(0)=0 ;
227 source_tAphi.
set(1)=0 ;
230 Cmp phifac = (F31-nphi()*F01)*x.
dsdr()
231 + (F32-nphi()*F02)*x.
srdsdt() ;
233 source_tAphi.
set(2)= -b_car()*a_car()/one_minus_x
234 *(tjphi-tnphi()*j_t + phifac)
235 + b_car()/(nnn()*nnn())*(tgrad1+tnphi()*grad2)
241 for (
int i=0; i<3; i++) {
242 Scalar tmp_filter = source_tAphi(i) ;
245 source_tAphi.
set(i) = tmp_filter ;
248 Tenseur WORK_VECT(mp, 1, CON, mp.get_bvect_cart()) ;
250 for (
int i=0; i<3; i++) {
251 WORK_VECT.
set(i) = 0 ;
255 WORK_SCAL.
set() = 0 ;
257 double lambda_mag = 0. ;
260 if (source_tAphi.
get_etat() != ETATZERO) {
262 for (
int i=0; i<3; i++) {
263 if(source_tAphi(i).dz_nonzero()) {
264 assert( source_tAphi(i).get_dzpuis() == 4 ) ;
267 (source_tAphi.
set(i)).set_dzpuis(4) ;
273 source_tAphi.
poisson_vect(lambda_mag, par_poisson_Avect, AVECT, WORK_VECT,
276 Cmp A_phi_n(AVECT(2));
281 Cmp source_A_1t(-a_car()*( j_t*gtt + j_phi*gtphi
283 + gtphi*(F31*x.
dsdr()+F32*x.
srdsdt()) )/one_minus_x
285 Scalar tmp_filter = source_A_1t ;
288 source_A_1t = tmp_filter ;
292 source_A_1t.
poisson(par_poisson_At, A_1t) ;
294 int L = mp.get_mg()->get_nt(0);
312 for (
int p=0; p<mp.get_mg()->get_np(0); p++) {
315 for(
int k=0;k<L;k++){
316 for(
int l=0;l<2*L;l++){
318 if(l==0) leg.
set(k,l)=1. ;
319 if(l==1) leg.
set(k,l)=
cos((*theta)(l_surf()(p,k),p,k,0)) ;
320 if(l>=2) leg.
set(k,l) = double(2*l-1)/double(l)
321 *
cos((*theta)(l_surf()(p,k),p,k,0))
322 * leg(k,l-1)-double(l-1)/double(l)*leg(k,l-2) ;
326 for(
int k=0;k<L;k++){
333 for(
int l=0;l<L;l++) MAT.
set(l,k) = leg(k,2*l)/
pow(Rsurf(k),2*l+1);
338 int* IPIV=
new int[L] ;
346 F77_dgesv(&L, &un, MAT.
t, &L, IPIV, VEC.
t, &L, &INFO) ;
350 for(
int k=0;k<L;k++) {VEC2.
set(k)=1. ; }
352 F77_dgesv(&L, &un, MAT_SAVE.
t, &L, IPIV, VEC2.
t, &L, &INFO) ;
356 for(
int nz=0;nz < Z; nz++){
357 for(
int i=0;i< mp.get_mg()->get_nr(nz);i++){
358 for(
int k=0;k<L;k++){
359 psi.
set(nz,p,k,i) = 0. ;
360 psi2.
set(nz,p,k,i) = 0. ;
361 for(
int l=0;l<L;l++){
362 psi.
set(nz,p,k,i) += VEC(l)*leg(k,2*l) /
363 pow((*mp.r.c)(nz,p,k,i),2*l+1);
364 psi2.
set(nz,p,k,i) += VEC2(l)*leg(k,2*l)/
365 pow((*mp.r.c)(nz, p, k,i),2*l+1);
381 Cmp A_t_ext(A_1t + psi) ;
382 A_t_ext.
annule(0,nzet-1) ;
388 for (
int j=0; j<nt; j++)
389 for (
int l=0; l<Z; l++)
390 for (
int i=0; i<mp.get_mg()->get_nr(l); i++)
391 A_0t.
set(l,0,j,i) = ( (*mp.r.c)(l,0,j,i) > Rsurf(j) ?
392 A_1t(l,0,j,i) + psi(l,0,j,i) : tmp(l,0,j,i) ) ;
403 double Q_0 = -4*M_PI*(*asymp[1])(Z-1,0,0,0) ;
411 double Q_2 = -4*M_PI*(*asymp[1])(Z-1,0,0,0) ;
419 double C = (Q-Q_0)/Q_2 ;
429 Cmp A_t_ext(A_0t + C*psi2) ;
430 A_t_ext.
annule(0,nzet-1) ;
436 for (
int j=0; j<nt; j++)
437 for (
int l=0; l<Z; l++)
438 for (
int i=0; i<mp.get_mg()->get_nr(l); i++)
439 A_t_n.
set(l,0,j,i) = ( (*mp.r.c)(l,0,j,i) > Rsurf(j) ?
440 A_0t(l,0,j,i) + C*psi2(l,0,j,i) :
441 A_0t(l,0,j,i) + C ) ;
455 A_t = relax_mag*A_t_n + (1.-relax_mag)*A_t ;
456 A_phi = relax_mag*A_phi_n + (1. - relax_mag)*A_phi ;
476 ApAp.
set().div_rsint() ;
477 ApAp.
set().div_rsint() ;
480 ApAt.
set().div_rsint() ;
482 E_em = 0.5*mu0 * ( 1/(a_car*nnn*nnn) * (AtAt + 2*tnphi*ApAt)
483 + ( (tnphi*tnphi/(a_car*nnn*nnn)) + 1/(a_car*b_car) )*ApAp );
484 Jp_em = -mu0 * (ApAt + tnphi*ApAp) /(a_car*nnn) ;
485 if (Jp_em.get_etat() != ETATZERO) Jp_em.set().mult_rsint() ;
497 Vector U_up(mp, CON, mp.get_bvect_cart()) ;
498 for (
int i=1; i<=3; i++)
499 U_up.
set(i) = u_euler(i-1) ;
502 Sym_tensor gamij(mp, COV, mp.get_bvect_spher()) ;
503 for (
int i=1; i<=3; i++)
504 for (
int j=1; j<i; j++) {
507 gamij.
set(1,1) = a_car() ;
508 gamij.
set(2,2) = a_car() ;
509 gamij.
set(3,3) = b_car() ;
514 Vector B_up(mp, CON, mp.get_bvect_spher()) ;
520 fac =
Scalar(gam_euler()*gam_euler()) ;
522 E_I = get_magnetisation() * EiEi / mu0 ;
524 J_I = get_magnetisation() * BiBi * Ui / mu0 ;
525 Sij_I = get_magnetisation()
526 * ( (BiBi / fac) * gamij + BiBi*Ui*Ui - Bi*Bi / fac ) / mu0 ;
528 for (
int i=1; i<=3; i++)
529 for (
int j=i; j<=3; j++)
530 Sij_I.set(i,j).set_dzpuis(0) ;
540 if (p_mass_g == 0x0) {
545 Tenseur SrrplusStt(
Cmp(Sij_I(1, 1) + Sij_I(2, 2)) ) ;
546 SrrplusStt = SrrplusStt / a_car ;
555 Tenseur source = nnn * (ener_euler + E_em + E_i
556 + s_euler + Spp_em + SrrplusStt + Spp) +
558 + 2 * bbb * (ener_euler + press) * tnphi * uuu ;
560 source = a_car * bbb * source ;
564 p_mass_g =
new double( source().integrale() ) ;
569 p_mass_g =
new double( mass_b() ) ;
584 if (p_angu_mom == 0x0) {
589 dens.
va = (dens.
va).mult_st() ;
592 dens = a_car() * (b_car() * (ener_euler() + press())
593 * dens + bbb() * (Jp_em() +
Cmp(J_I(3)) ) ) ;
596 dens = nbar() * dens ;
601 p_angu_mom =
new double( dens.
integrale() ) ;
623 Tenseur sou_m = 2 * qpig * a_car * (press + (ener_euler+press)
626 Tenseur sou_q = 2 * qpig * a_car * Spp_em + 1.5 * ak_car
629 p_grv2 =
new double(
double(1) - lambda_grv2(sou_m(), sou_q()) ) ;
659 source = 0.75 * ak_car
661 logn.gradient_spher() )
665 Cmp aa = alpha() - 0.5 * beta() ;
671 cout <<
"Etoile_rot::grv3: the mapping does not belong"
672 <<
" to the class Map_radial !" << endl ;
679 vdaadt = vdaadt.
ssint() ;
684 temp = ( bbb() - a_car()/bbb() ) * temp ;
692 vtemp = (mpr->
xsr) * vtemp ;
700 source = bbb() * source() + 0.5 * temp ;
705 logn.gradient_spher() ) ;
710 double int_grav = source().integrale() ;
718 Tenseur SrrplusStt(
Cmp(Sij_I(1, 1) + Sij_I(2, 2)) ) ;
719 SrrplusStt = SrrplusStt / a_car ;
724 source = qpig * a_car * bbb * ( s_euler + Spp_em + SrrplusStt + Spp ) ;
727 source = qpig * ( 3 * press + nbar * uuu * uuu ) ;
732 double int_mat = source().integrale() ;
737 *ost <<
"Et_magnetisation::grv3 : gravitational term : " << int_grav
739 *ost <<
"Et_magnetisation::grv3 : matter term : " << int_mat
743 p_grv3 =
new double( (int_grav + int_mat) / int_mat ) ;
757 if (p_mom_quad_old == 0x0) {
769 Tenseur SrrplusStt(
Cmp(Sij_I(1, 1) + Sij_I(2, 2)) ) ;
770 SrrplusStt = SrrplusStt / a_car ;
780 source = qpig * a_car *( ener_euler + E_em + E_i
781 + s_euler + Spp_em + SrrplusStt + Spp)
786 source = qpig * nbar ;
796 Cmp& csource = source.
set() ;
808 assert(
dynamic_cast<const Map_radial*
>(&mp) != 0x0 ) ;
818 source = 0.5 * source() - 1.5 * temp ;
822 p_mom_quad_old =
new double( - source().integrale() / qpig ) ;
824 return *p_mom_quad_old ;
832 if (p_mom_quad_Bo == 0x0) {
835 Tenseur SrrplusStt(
Cmp(Sij_I(1, 1) + Sij_I(2, 2)) ) ;
837 Cmp dens = a_car() * press() ;
838 dens = bbb() * nnn() * (SrrplusStt() + 2*dens) ;
842 p_mom_quad_Bo =
new double( - 16. * dens.
integrale() / qpig ) ;
844 return *p_mom_quad_Bo ;
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
void mult_rsint()
Multiplication by .
void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
void div_r()
Division by r everywhere.
void dec_dzpuis()
Decreases by 1 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
int get_etat() const
Returns the logical state.
Valeur va
The numerical value of the Cmp
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
void annule(int l)
Sets the Cmp to zero in a given domain.
const Cmp & laplacien(int zec_mult_r=4) const
Returns the Laplacian of *this.
int get_dzpuis() const
Returns dzpuis.
void mult_r()
Multiplication by r everywhere.
void inc2_dzpuis()
Increases by 2 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
void set_dzpuis(int)
Set a value to dzpuis.
double integrale() const
Computes the integral over all space of *this .
const Cmp & srdsdt() const
Returns of *this .
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...
Valeur ** asymptot(int n, const int flag=0) const
Asymptotic expansion at r = infinity.
Cmp poisson() const
Solves the scalar Poisson equation with *this as a source.
const Cmp & dsdr() const
Returns of *this .
Tbl & set(int l)
Read/write of the value in a given domain.
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,...
virtual double grv3(ostream *ost=0x0) const
Error on the virial identity GRV3.
virtual double mom_quad_Bo() const
Part of the quadrupole moment.
virtual double grv2() const
Error on the virial identity GRV2.
virtual double mom_quad_old() const
Part of the quadrupole moment.
virtual void MHD_comput()
Computes the electromagnetic part of the stress-energy tensor.
virtual double angu_mom() const
Angular momentum.
virtual double mass_g() const
Gravitational mass.
Base class for pure radial mappings.
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain.
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...
Metric for tensor calculation.
Tensor field of valence 0 (or component of a tensorial field).
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.
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.
Class intended to describe valence-2 symmetric tensors.
double & set(int i)
Read/write of a particular element (index i) (1D case)
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
double * t
The array of double.
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
void set_std_base()
Set the standard spectal basis of decomposition for each component.
const Tenseur & gradient_spher() const
Returns the gradient of *this (Spherical coordinates) (scalar field only).
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
void poisson_vect(double lambda, Param &par, Tenseur &shift, Tenseur &vect, Tenseur &scal) const
Solves the vectorial Poisson equation : .
int get_etat() const
Returns the logical state.
Values and coefficients of a (real-value) function.
const Valeur & mult_ct() const
Returns applied to *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...
const Valeur & sx() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR )
const Valeur & ssint() const
Returns of *this.
Tensor field of valence 1.
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Scalar & set(int)
Read/write access to a component.
Cmp sqrt(const Cmp &)
Square root.
Cmp pow(const Cmp &, int)
Power .
Cmp cos(const Cmp &)
Cosine.
Cmp log(const Cmp &)
Neperian logarithm.
Tensor down(int ind, const Metric &gam) const
Computes a new tensor by lowering an index of *this.
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
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...
Tenseur flat_scalar_prod_desal(const Tenseur &t1, const Tenseur &t2)
Same as flat_scalar_prod but with desaliasing.
Standard electro-magnetic units.
Standard units of space, time and mass.