30char et_rot_mag_mag_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_mag_mag.C,v 1.17 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"
107Cmp raccord_c1(
const Cmp& uu,
int l1) ;
112 Cmp (*f_j)(
const Cmp&,
const double),
113 Param& par_poisson_At,
114 Param& par_poisson_Avect){
115 double relax_mag = 0.5 ;
119 mag_filter = par_poisson_At.
get_int(1) ;
124 bool adapt(adapt_flag) ;
137 assert (mpr != 0x0) ;
138 for (
int j=0; j<nt; j++)
175 BLAH -= grad3 + npgrada ;
188 for (
int j=0; j<nt; j++)
189 for (
int l=0; l<
nzet; l++)
190 for (
int i=0; i<
mp.
get_mg()->get_nr(l); i++)
191 j_t.
set(l,0,j,i) = ( (*
mp.
r.
c)(l,0,j,i) > Rsurf(j) ?
192 0. : tmp(l,0,j,i) ) ;
215 d_grad4.div_rsint() ;
216 source_tAphi.
set(0)=0 ;
217 source_tAphi.
set(1)=0 ;
222 if (mag_filter == 1) {
223 for (
int i=0; i<3; i++) {
224 Scalar tmp_filter = source_tAphi(i) ;
227 source_tAphi.
set(i) = tmp_filter ;
234 for (
int i=0; i<3; i++) {
235 WORK_VECT.
set(i) = 0 ;
239 WORK_SCAL.
set() = 0 ;
241 double lambda_mag = 0. ;
244 if (source_tAphi.
get_etat() != ETATZERO) {
246 for (
int i=0; i<3; i++) {
247 if(source_tAphi(i).dz_nonzero()) {
248 assert( source_tAphi(i).get_dzpuis() == 4 ) ;
251 (source_tAphi.
set(i)).set_dzpuis(4) ;
256 source_tAphi.
poisson_vect(lambda_mag, par_poisson_Avect, AVECT, WORK_VECT,
259 Cmp A_phi_n(AVECT(2));
265 if (mag_filter == 1) {
266 Scalar tmp_filter = source_A_1t ;
269 source_A_1t = tmp_filter ;
274 source_A_1t.
poisson(par_poisson_At, A_1t) ;
298 for (
int i=0; i<L; i++)
299 VEC3.
set(i) = 1. / double(i+1) ;
301 for (
int p=0; p<
mp.
get_mg()->get_np(0); p++) {
304 for(
int k=0;k<L;k++){
305 for(
int l=0;l<2*L;l++){
307 if(l==0) leg.
set(k,l)=1. ;
308 if(l==1) leg.
set(k,l)=
cos((*theta)(
l_surf()(p,k),p,k,0)) ;
309 if(l>=2) leg.
set(k,l) = double(2*l-1)/double(l)
311 * leg(k,l-1)-double(l-1)/double(l)*leg(k,l-2) ;
315 for(
int k=0;k<L;k++){
322 for(
int l=0;l<L;l++) MAT.
set(l,k) = leg(k,2*l)/
pow(Rsurf(k),2*l+1);
327 int* IPIV=
new int[L] ;
335 F77_dgesv(&L, &un, MAT.
t, &L, IPIV, VEC.
t, &L, &INFO) ;
339 for(
int k=0;k<L;k++) {VEC2.
set(k)=1. ; }
341 F77_dgesv(&L, &un, MAT_SAVE.
t, &L, IPIV, VEC2.
t, &L, &INFO) ;
345 for(
int nz=0;nz < Z; nz++){
346 for(
int i=0;i<
mp.
get_mg()->get_nr(nz);i++){
347 for(
int k=0;k<L;k++){
348 psi.
set(nz,p,k,i) = 0. ;
349 psi2.
set(nz,p,k,i) = 0. ;
350 psi3.
set(nz,p,k,i) = 0. ;
351 for(
int l=0;l<L;l++){
352 psi.
set(nz,p,k,i) += VEC(l)*leg(k,2*l) /
353 pow((*
mp.
r.
c)(nz,p,k,i),2*l+1);
354 psi2.
set(nz,p,k,i) += VEC2(l)*leg(k,2*l)/
355 pow((*
mp.
r.
c)(nz, p, k,i),2*l+1);
356 psi3.
set(nz,p,k,i) += VEC3(l)*leg(k,2*l)/
357 (
pow((*
mp.
r.
c)(nz, p, k,i),2*l+1)) ;
373 Cmp A_t_ext(A_1t + psi) ;
380 for (
int j=0; j<nt; j++)
381 for (
int l=0; l<Z; l++)
382 for (
int i=0; i<
mp.
get_mg()->get_nr(l); i++)
383 A_0t.
set(l,0,j,i) = ( (*
mp.
r.
c)(l,0,j,i) > Rsurf(j) ?
384 A_1t(l,0,j,i) + psi(l,0,j,i) : tmp(l,0,j,i) ) ;
388 if (mag_filter == 1) {
389 Scalar tmp_filter = A_0t ;
397 double Q_0 = -4*M_PI*(*asymp[1])(Z-1,0,0,0) ;
405 double Q_2 = -4*M_PI*(*asymp[1])(Z-1,0,0,0) ;
413 double C = (
Q-Q_0)/Q_2 ;
423 Cmp A_t_ext(A_0t + C*psi2) ;
430 for (
int j=0; j<nt; j++)
431 for (
int l=0; l<Z; l++)
432 for (
int i=0; i<
mp.
get_mg()->get_nr(l); i++) {
433 A_t_n.
set(l,0,j,i) = ( (*
mp.
r.
c)(l,0,j,i) > Rsurf(j) ?
434 A_0t(l,0,j,i) + C*psi2(l,0,j,i) :
435 A_0t(l,0,j,i) + C ) ;
439 if (mag_filter == 1) {
440 Scalar tmp_filter = A_t_n ;
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 ;
503 BLAH -= grad3 + npgrada ;
507 Cmp source_A_t_n(
mp);
518 source_A_t_n.
poisson(par_poisson_At, A_t_n) ;
534 d_grad4.div_rsint() ;
535 source_tAphi.
set(0)=0 ;
536 source_tAphi.
set(1)=0 ;
542 source_tAphi.
set(2)= - tjphi ;}
548 for (
int i=0; i<3; i++) {
549 WORK_VECT.
set(i) = 0 ;
553 WORK_SCAL.
set() = 0 ;
555 double lambda_mag = 0. ;
558 if (source_tAphi.
get_etat() != ETATZERO) {
560 for (
int i=0; i<3; i++) {
561 if(source_tAphi(i).dz_nonzero()) {
562 assert( source_tAphi(i).get_dzpuis() == 4 ) ;
565 (source_tAphi.
set(i)).set_dzpuis(4) ;
570 source_tAphi.
poisson_vect(lambda_mag, par_poisson_Avect, AVECT, WORK_VECT,
573 Cmp A_phi_n(AVECT(2));
578 A_t = relax_mag*A_t_n + (1.-relax_mag)*
A_t ;
579 A_phi = relax_mag*A_phi_n + (1. - relax_mag)*
A_phi ;
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...
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 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.
const Cmp & srdsdt() const
Returns of *this .
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.
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.
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,...
Cmp A_t
t-component of the elecctromagnetic potential 1-form, divided by .
bool is_conduct() const
Tells if the star is made of conducting or isolating material.
Cmp j_t
t-component of the current 4-vector
double Q
In the case of a perfect conductor, the requated baryonic charge.
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...
Tenseur nphi
Metric coefficient .
Tenseur bbb
Metric factor B.
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.
Tenseur nnn
Total lapse function.
Tenseur nbar
Baryon density in the fluid frame.
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.
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.
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.
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.
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
int get_n_int() const
Returns the number of stored int 's addresses.
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.
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 & ssint() const
Returns of *this.
Cmp pow(const Cmp &, int)
Power .
Cmp cos(const Cmp &)
Cosine.
Cmp log(const Cmp &)
Neperian logarithm.
Tenseur flat_scalar_prod_desal(const Tenseur &t1, const Tenseur &t2)
Same as flat_scalar_prod but with desaliasing.