29char 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"
98 Cmp (*f_j)(
const Cmp&,
const double),
99 Param& par_poisson_At,
100 Param& par_poisson_Avect){
101 double relax_mag = 0.5 ;
105 bool adapt(adapt_flag) ;
118 assert (mpr != 0x0) ;
119 for (
int j=0; j<nt; j++)
156 BLAH -= grad3 + npgrada ;
184 Cmp one_minus_x = 1 - x ;
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++)
201 j_t.
set(l,0,j,i) = ( (*
mp.
r.
c)(l,0,j,i) > Rsurf(j) ?
202 0. : tmp(l,0,j,i) ) ;
225 d_grad4.div_rsint() ;
226 source_tAphi.
set(0)=0 ;
227 source_tAphi.
set(1)=0 ;
241 for (
int i=0; i<3; i++) {
242 Scalar tmp_filter = source_tAphi(i) ;
245 source_tAphi.
set(i) = tmp_filter ;
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));
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) ;
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)
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) ;
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) ;
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() ;
498 for (
int i=1; i<=3; i++)
503 for (
int i=1; i<=3; i++)
504 for (
int j=1; j<i; j++) {
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++)
546 SrrplusStt = SrrplusStt /
a_car ;
564 p_mass_g =
new double( source().integrale() ) ;
589 dens.
va = (dens.
va).mult_st() ;
596 dens =
nbar() * dens ;
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() ;
692 vtemp = (mpr->
xsr) * vtemp ;
700 source =
bbb() * source() + 0.5 * temp ;
710 double int_grav = source().integrale() ;
719 SrrplusStt = SrrplusStt /
a_car ;
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 ) ;
770 SrrplusStt = SrrplusStt /
a_car ;
786 source = qpig *
nbar ;
796 Cmp& csource = source.
set() ;
818 source = 0.5 * source() - 1.5 * temp ;
838 dens =
bbb() *
nnn() * (SrrplusStt() + 2*dens) ;
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...
Tbl & set(int l)
Read/write of the value in a given domain.
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 .
void fait() const
Computes, at each point of the grid, the value of the coordinate or mapping derivative represented by...
Mtbl * c
The coordinate values at each grid point.
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.
Sym_tensor Sij_I
Interaction stress 3-tensor.
Vector J_I
Interaction momentum density 3-vector.
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.
const Scalar & get_magnetisation() const
Accessor to the magnetisation scalar field.
virtual void MHD_comput()
Computes the electromagnetic part of the stress-energy tensor.
virtual double angu_mom() const
Angular momentum.
Scalar E_I
Interaction (magnetisation) energy density.
virtual double mass_g() const
Gravitational mass.
Tenseur Srr_em
rr component of the electromagnetic stress 3-tensor, as measured in the Eulerian frame....
Cmp j_phi
-component of the current 4-vector
Cmp A_phi
-component of the electromagnetic potential 1-form divided by .
double a_j
Amplitude of the curent/charge function.
Tenseur Spp_em
component of the electromagnetic stress 3-tensor, as measured in the Eulerian frame.
Cmp A_t
t-component of the elecctromagnetic potential 1-form, divided by .
Cmp j_t
t-component of the current 4-vector
Tenseur E_em
electromagnetic energy density in the Eulerian frame
Tenseur Jp_em
component of the electromagnetic momentum density 3-vector, as measured in the Eulerian frame.
double Q
In the case of a perfect conductor, the requated baryonic charge.
Tenseur Elec() const
Computes the electric field spherical components in Lorene's units.
Tenseur Magn() const
Computes the magnetic field spherical components in Lorene's units.
Tenseur uuu
Norm of u_euler.
double omega
Rotation angular velocity ([f_unit] )
Tenseur & logn
Metric potential = logn_auto.
virtual const Itbl & l_surf() const
Description of the stellar surface: returns a 2-D Itbl containing the values of the domain index l on...
double * p_mom_quad_old
Part of the quadrupole moment.
Tenseur nphi
Metric coefficient .
virtual double mass_b() const
Baryon mass.
static double lambda_grv2(const Cmp &sou_m, const Cmp &sou_q)
Computes the coefficient which ensures that the GRV2 virial identity is satisfied.
Tenseur bbb
Metric factor B.
Tenseur & dzeta
Metric potential = beta_auto.
double * p_grv3
Error on the virial identity GRV3.
double * p_grv2
Error on the virial identity GRV2.
double * p_mom_quad_Bo
Part of the quadrupole moment.
double * p_angu_mom
Angular momentum.
Tenseur b_car
Square of the metric factor B.
Tenseur tnphi
Component of the shift vector.
int nzet
Number of domains of *mp occupied by the star.
double * p_mass_g
Gravitational mass.
Tenseur nnn
Total lapse function.
Tenseur nbar
Baryon density in the fluid frame.
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Tenseur gam_euler
Lorentz factor between the fluid and Eulerian observers.
Map & mp
Mapping associated with the star.
Tenseur ener
Total energy density in the fluid frame.
Tenseur press
Fluid pressure.
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one.
Tenseur ener_euler
Total energy density in the Eulerian frame.
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
const Tbl & xi_surf() const
Description of the stellar surface: returns a 2-D Tbl containing the values of the radial coordinat...
Tenseur a_car
Total conformal factor .
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...
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Coord r
r coordinate centered on the grid
Coord tet
coordinate centered on the grid
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Metric for tensor calculation.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
int get_nzone() const
Returns the number of domains.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
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.
void set_dzpuis(int)
Modifies the dzpuis flag.
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.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
double & set(int i)
Read/write of a particular element (index i) (1D case)
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.