LORENE
et_magnetisation_comp.C
1 /*
2  * Computational functions for magnetized rotating equilibrium
3  *
4  * (see file et_rot_mag.h for documentation)
5  *
6  */
7 
8 /*
9  * Copyright (c) 2013 Debarati Chatterjee, Jerome Novak
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License as published by
15  * the Free Software Foundation; either version 2 of the License, or
16  * (at your option) any later version.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
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 $" ;
30 
31 /*
32  * $Id: et_magnetisation_comp.C,v 1.14 2016/11/01 09:12:59 j_novak Exp $
33  * $Log: et_magnetisation_comp.C,v $
34  * Revision 1.14 2016/11/01 09:12:59 j_novak
35  * Correction of a missing '-' in mom_quad_old().
36  *
37  * Revision 1.13 2015/06/12 12:38:25 j_novak
38  * Implementation of the corrected formula for the quadrupole momentum.
39  *
40  * Revision 1.12 2014/10/21 09:23:54 j_novak
41  * Addition of global functions mass_g(), angu_mom(), grv2/3() and mom_quad().
42  *
43  * Revision 1.11 2014/10/13 08:52:57 j_novak
44  * Lorene classes and functions now belong to the namespace Lorene.
45  *
46  * Revision 1.10 2014/07/04 12:08:02 j_novak
47  * Added some filtering.
48  *
49  * Revision 1.9 2014/05/14 15:19:05 j_novak
50  * The magnetisation field is now filtered.
51  *
52  * Revision 1.8 2014/05/13 15:37:12 j_novak
53  * Updated to new magnetic units.
54  *
55  * Revision 1.7 2014/05/01 13:07:16 j_novak
56  * Fixed two bugs: in the computation of F31,F32 and the triad of U_up.
57  *
58  * Revision 1.6 2014/04/29 13:46:07 j_novak
59  * Addition of switches 'use_B_in_eos' and 'include_magnetisation' to control the model.
60  *
61  * Revision 1.5 2014/04/28 14:53:29 j_novak
62  * Minor modif.
63  *
64  * Revision 1.4 2014/04/28 12:48:13 j_novak
65  * Minor modifications.
66  *
67  * Revision 1.2 2013/12/19 17:05:40 j_novak
68  * Corrected a dzpuis problem.
69  *
70  * Revision 1.1 2013/12/13 16:36:51 j_novak
71  * Addition and computation of magnetisation terms in the Einstein equations.
72  *
73  *
74  *
75  * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_magnetisation_comp.C,v 1.14 2016/11/01 09:12:59 j_novak Exp $
76  *
77  */
78 
79 // Headers C
80 #include <cstdlib>
81 #include <cmath>
82 
83 // Headers Lorene
84 #include "et_rot_mag.h"
85 #include "metric.h"
86 #include "utilitaires.h"
87 #include "param.h"
88 #include "proto_f77.h"
89 #include "unites.h"
90 
91 namespace Lorene {
92 
93  using namespace Unites_mag ;
94 
95 // Algo du papier de 1995
96 
97 void Et_magnetisation::magnet_comput(const int adapt_flag,
98  Cmp (*f_j)(const Cmp&, const double),
99  Param& par_poisson_At,
100  Param& par_poisson_Avect){
101  double relax_mag = 0.5 ;
102 
103  int Z = mp.get_mg()->get_nzone();
104 
105  bool adapt(adapt_flag) ;
106  /****************************************************************
107  * Assertion that all zones have same number of points in theta
108  ****************************************************************/
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) ;
111 
112  Tbl Rsurf(nt) ;
113  Rsurf.set_etat_qcq() ;
114  mp.r.fait() ;
115  mp.tet.fait() ;
116  Mtbl* theta = mp.tet.c ;
117  const Map_radial* mpr = dynamic_cast<const Map_radial*>(&mp) ;
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) ;
121 
122 
123  // Calcul de A_0t dans l'etoile (conducteur parfait)
124 
125  Cmp A_0t(- omega * A_phi) ;
126  A_0t.annule(nzet,Z-1) ;
127 
128  Tenseur ATTENS(A_t) ;
129  Tenseur APTENS(A_phi) ;
130  Tenseur BMN(-logn) ;
131  BMN = BMN + log(bbb) ;
132  BMN.set_std_base() ;
133 
134 
136  nphi.gradient_spher())());
138  nphi.gradient_spher())()) ;
140  BMN.gradient_spher())()
141  + 2*nphi()*flat_scalar_prod_desal(APTENS.gradient_spher(),
142  BMN.gradient_spher())()) ;
143 
144  Cmp ATANT(A_phi.srdsdt()); // Constrction par copie pour mapping
145 
146  ATANT.va = ATANT.va.mult_ct().ssint() ;
147 
148  Cmp ttnphi(tnphi()) ;
149  ttnphi.mult_rsint() ;
150  Cmp BLAH(- b_car()/(nnn()*nnn())*ttnphi*grad1) ;
151  BLAH -= (1+b_car()/(nnn()*nnn())*tnphi()*tnphi())*grad2 ;
152  Cmp nphisr(nphi()) ;
153  nphisr.div_r() ;
154  Cmp npgrada(2*nphisr*(A_phi.dsdr()+ATANT )) ;
155  npgrada.inc2_dzpuis() ;
156  BLAH -= grad3 + npgrada ;
157  Cmp gtt(-nnn()*nnn()+b_car()*tnphi()*tnphi()) ;
158  Cmp gtphi( - b_car()*ttnphi) ;
159 
160  // Computation of j_t thanks to Maxwell-Gauss
161  // modified to include Magnetisation
162  // components of F
163  Cmp F01 = 1/(a_car()*nnn()*nnn())*A_0t.dsdr()
164  + 1/(a_car()*nnn()*nnn())*nphi()*A_phi.dsdr() ;
165 
166  Cmp F02 = 1/(a_car()*nnn()*nnn())*A_0t.srdsdt()
167  + 1/(a_car()*nnn()*nnn())*nphi()*A_phi.srdsdt() ;
168 
169  Cmp tmp = A_phi.dsdr() / (bbb() * bbb() * a_car() );
170  tmp.div_rsint() ;
171  tmp.div_rsint() ;
172  Cmp F31 = 1/(a_car()*nnn()*nnn())*nphi()*nphi()*A_phi.dsdr()
173  + 1/(a_car()*nnn()*nnn())*nphi()*A_0t.dsdr()
174  + tmp ;
175 
176  tmp = A_phi.srdsdt() / (bbb() * bbb() * a_car() );
177  tmp.div_rsint() ;
178  tmp.div_rsint() ;
179  Cmp F32 = 1/(a_car()*nnn()*nnn())*nphi()*nphi()*A_phi.srdsdt()
180  + 1/(a_car()*nnn()*nnn())*nphi()*A_0t.srdsdt()
181  + tmp ;
182 
183  Cmp x = get_magnetisation();
184  Cmp one_minus_x = 1 - x ;
185  one_minus_x.std_base_scal() ;
186 
187  tmp = ((BLAH - A_0t.laplacien())*one_minus_x/a_car()
188  - gtphi*j_phi
189  - gtt*(F01*x.dsdr()+F02*x.srdsdt())
190  - gtphi*(F31*x.dsdr()+F32*x.srdsdt()) ) / gtt ;
191 
192  tmp.annule(nzet, Z-1) ;
193  if (adapt) {
194  j_t = tmp ;
195  }
196  else {
197  j_t.allocate_all() ;
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) ;
204  }
205  j_t.std_base_scal() ;
206 
207  // Calcul du courant j_phi
208  j_phi = omega * j_t + (ener() + press())*f_j(A_phi, a_j) ;
209  j_phi.std_base_scal() ;
210 
211  // Resolution de Maxwell Ampere (-> A_phi)
212  // Calcul des termes sources avec A-t du pas precedent.
213 
215  BMN.gradient_spher())());
216 
217  Tenseur source_tAphi(mp, 1, CON, mp.get_bvect_spher()) ;
218 
219  source_tAphi.set_etat_qcq() ;
220  Cmp tjphi(j_phi) ;
221  tjphi.mult_rsint() ;
222  Cmp tgrad1(grad1) ;
223  tgrad1.mult_rsint() ;
224  Cmp d_grad4(grad4) ;
225  d_grad4.div_rsint() ;
226  source_tAphi.set(0)=0 ;
227  source_tAphi.set(1)=0 ;
228 
229 // modified to include Magnetisation
230  Cmp phifac = (F31-nphi()*F01)*x.dsdr()
231  + (F32-nphi()*F02)*x.srdsdt() ;
232  phifac.mult_rsint();
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)
236  + d_grad4 ;
237 
238  source_tAphi.change_triad(mp.get_bvect_cart());
239 
240  // Filtering
241  for (int i=0; i<3; i++) {
242  Scalar tmp_filter = source_tAphi(i) ;
243  tmp_filter.exponential_filter_r(0, 2, 1) ;
244  tmp_filter.exponential_filter_ylm(0, 2, 1) ;
245  source_tAphi.set(i) = tmp_filter ;
246  }
247 
248  Tenseur WORK_VECT(mp, 1, CON, mp.get_bvect_cart()) ;
249  WORK_VECT.set_etat_qcq() ;
250  for (int i=0; i<3; i++) {
251  WORK_VECT.set(i) = 0 ;
252  }
253  Tenseur WORK_SCAL(mp) ;
254  WORK_SCAL.set_etat_qcq() ;
255  WORK_SCAL.set() = 0 ;
256 
257  double lambda_mag = 0. ; // No 3D version !
258 
259  Tenseur AVECT(source_tAphi) ;
260  if (source_tAphi.get_etat() != ETATZERO) {
261 
262  for (int i=0; i<3; i++) {
263  if(source_tAphi(i).dz_nonzero()) {
264  assert( source_tAphi(i).get_dzpuis() == 4 ) ;
265  }
266  else{
267  (source_tAphi.set(i)).set_dzpuis(4) ;
268  }
269  }
270 
271  }
272 
273  source_tAphi.poisson_vect(lambda_mag, par_poisson_Avect, AVECT, WORK_VECT,
274  WORK_SCAL) ;
275  AVECT.change_triad(mp.get_bvect_spher());
276  Cmp A_phi_n(AVECT(2));
277  A_phi_n.mult_rsint() ;
278 
279  // Solution to Maxwell-Ampere : A_1
280  // modified to include Magnetisation
281  Cmp source_A_1t(-a_car()*( j_t*gtt + j_phi*gtphi
282  + gtt*(F01*x.dsdr()+F02*x.srdsdt())
283  + gtphi*(F31*x.dsdr()+F32*x.srdsdt()) )/one_minus_x
284  + BLAH);
285  Scalar tmp_filter = source_A_1t ;
286  tmp_filter.exponential_filter_r(0, 2, 1) ;
287  tmp_filter.exponential_filter_ylm(0, 2, 1) ;
288  source_A_1t = tmp_filter ;
289 
290  Cmp A_1t(mp);
291  A_1t = 0 ;
292  source_A_1t.poisson(par_poisson_At, A_1t) ;
293 
294  int L = mp.get_mg()->get_nt(0);
295 
296  Tbl MAT(L,L) ;
297  Tbl MAT_PHI(L,L);
298  Tbl VEC(L) ;
299 
300  MAT.set_etat_qcq() ;
301  VEC.set_etat_qcq() ;
302  MAT_PHI.set_etat_qcq() ;
303 
304  Tbl leg(L,2*L) ;
305  leg.set_etat_qcq() ;
306 
307  Cmp psi(mp);
308  Cmp psi2(mp);
309  psi.allocate_all() ;
310  psi2.allocate_all() ;
311 
312  for (int p=0; p<mp.get_mg()->get_np(0); p++) {
313  // leg[k,l] : legendre_l(cos(theta_k))
314  // Construction par recurrence de degre 2
315  for(int k=0;k<L;k++){
316  for(int l=0;l<2*L;l++){
317 
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) ;
323  }
324  }
325 
326  for(int k=0;k<L;k++){
327 
328  // Valeurs a la surface trouvees via va.val_point_jk(l,xisurf,k,p)
329 
330  VEC.set(k) = A_0t.va.val_point_jk(l_surf()(p,k), xi_surf()(p,k), k, p)
331  -A_1t.va.val_point_jk(l_surf()(p,k), xi_surf()(p,k), k, p);
332 
333  for(int l=0;l<L;l++) MAT.set(l,k) = leg(k,2*l)/pow(Rsurf(k),2*l+1);
334 
335  }
336  // appel fortran :
337 
338  int* IPIV=new int[L] ;
339  int INFO ;
340 
341  Tbl MAT_SAVE(MAT) ;
342  Tbl VEC2(L) ;
343  VEC2.set_etat_qcq() ;
344  int un = 1 ;
345 
346  F77_dgesv(&L, &un, MAT.t, &L, IPIV, VEC.t, &L, &INFO) ;
347 
348  // coeffs a_l dans VEC
349 
350  for(int k=0;k<L;k++) {VEC2.set(k)=1. ; }
351 
352  F77_dgesv(&L, &un, MAT_SAVE.t, &L, IPIV, VEC2.t, &L, &INFO) ;
353 
354  delete [] IPIV ;
355 
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);
366  }
367  }
368  }
369  }
370  }
371  psi.std_base_scal() ;
372  psi2.std_base_scal() ;
373 
374  assert(psi.get_dzpuis() == 0) ;
375  int dif = A_1t.get_dzpuis() ;
376  if (dif > 0) {
377  for (int d=0; d<dif; d++) A_1t.dec_dzpuis() ;
378  }
379 
380  if (adapt) {
381  Cmp A_t_ext(A_1t + psi) ;
382  A_t_ext.annule(0,nzet-1) ;
383  A_0t += A_t_ext ;
384  }
385  else {
386  tmp = A_0t ;
387  A_0t.allocate_all() ;
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) ) ;
393  }
394  A_0t.std_base_scal() ;
395 
396  tmp_filter = A_0t ;
397  tmp_filter.exponential_filter_r(0, 2, 1) ;
398  tmp_filter.exponential_filter_ylm(0, 2, 1) ;
399  A_0t = tmp_filter ;
400 
401  Valeur** asymp = A_0t.asymptot(1) ;
402 
403  double Q_0 = -4*M_PI*(*asymp[1])(Z-1,0,0,0) ; // utilise A_0t plutot que E
404  delete asymp[0] ;
405  delete asymp[1] ;
406 
407  delete [] asymp ;
408 
409  asymp = psi2.asymptot(1) ;
410 
411  double Q_2 = -4*M_PI*(*asymp[1])(Z-1,0,0,0) ; // A_2t = psi2 a l'infini
412  delete asymp[0] ;
413  delete asymp[1] ;
414 
415  delete [] asymp ;
416 
417  // solution definitive de A_t:
418 
419  double C = (Q-Q_0)/Q_2 ;
420 
421  assert(psi2.get_dzpuis() == 0) ;
422  dif = A_0t.get_dzpuis() ;
423  if (dif > 0) {
424  for (int d=0; d<dif; d++) A_0t.dec_dzpuis() ;
425  }
426  Cmp A_t_n(mp) ;
427  if (adapt) {
428  A_t_n = A_0t + C ;
429  Cmp A_t_ext(A_0t + C*psi2) ;
430  A_t_ext.annule(0,nzet-1) ;
431  A_t_n.annule(nzet,Z-1) ;
432  A_t_n += A_t_ext ;
433  }
434  else {
435  A_t_n.allocate_all() ;
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 ) ;
442  }
443  A_t_n.std_base_scal() ;
444  tmp_filter = A_t_n ;
445  tmp_filter.exponential_filter_r(0, 2, 1) ;
446  tmp_filter.exponential_filter_ylm(0, 2, 1) ;
447  A_t_n = tmp_filter ;
448 
449  asymp = A_t_n.asymptot(1) ;
450 
451  delete asymp[0] ;
452  delete asymp[1] ;
453 
454  delete [] asymp ;
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 ;
457 
458 }
459 
460 
462  // Computes the E-M terms of the stress-energy tensor...
463 
464  Tenseur ATTENS(A_t) ;
465 
466  Tenseur APTENS(A_phi) ;
467 
469  APTENS.gradient_spher())() );
471  ATTENS.gradient_spher())() );
473  ATTENS.gradient_spher())() );
474 
475  if (ApAp.get_etat() != ETATZERO) {
476  ApAp.set().div_rsint() ;
477  ApAp.set().div_rsint() ;
478  }
479  if (ApAt.get_etat() != ETATZERO)
480  ApAt.set().div_rsint() ;
481 
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() ;
486  Srr_em = 0 ;
487  // Stt_em = -Srr_em
488  Spp_em = E_em ;
489 
490  // ... and those corresponding to the magnetization.
491  Tenseur Efield = Elec() ;
492  Tenseur Bfield = Magn() ;
493 
494  Scalar EiEi ( flat_scalar_prod(Efield, Efield)() ) ;
495  Scalar BiBi ( flat_scalar_prod(Bfield, Bfield)() ) ;
496 
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) ;
500  U_up.change_triad(mp.get_bvect_spher()) ;
501 
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++) {
505  gamij.set(i,j) = 0 ;
506  }
507  gamij.set(1,1) = a_car() ;
508  gamij.set(2,2) = a_car() ;
509  gamij.set(3,3) = b_car() ;
510  Metric met(gamij) ;
511  Vector Ui = U_up.down(0, met) ;
512 
513  Scalar fac = sqrt(a_car()) ;
514  Vector B_up(mp, CON, mp.get_bvect_spher()) ;
515  B_up.set(1) = Scalar(Bfield(0)) / fac ;
516  B_up.set(2) = Scalar(Bfield(1)) / fac ;
517  B_up.set(3) = 0 ;
518  Vector Bi = B_up.down(0, met) ;
519 
520  fac = Scalar(gam_euler()*gam_euler()) ;
521 
522  E_I = get_magnetisation() * EiEi / mu0 ;
523 
524  J_I = get_magnetisation() * BiBi * Ui / mu0 ;
525  Sij_I = get_magnetisation()
526  * ( (BiBi / fac) * gamij + BiBi*Ui*Ui - Bi*Bi / fac ) / mu0 ;
527 
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) ;
531 
532 }
533 
534  //----------------------------//
535  // Gravitational mass //
536  //----------------------------//
537 
538 double Et_magnetisation::mass_g() const {
539 
540  if (p_mass_g == 0x0) { // a new computation is required
541 
542  if (relativistic) {
543 
544  // Magnetisation: S_{rr} + S_{\theta\theta}
545  Tenseur SrrplusStt( Cmp(Sij_I(1, 1) + Sij_I(2, 2)) ) ;
546  SrrplusStt = SrrplusStt / a_car ; // S^r_r + S^\theta_\theta
547 
548  Tenseur Spp (Cmp(Sij_I(3, 3))) ; // Magnetisation: S_{\phi\phi}
549  Spp = Spp / b_car ; // S^\phi_\phi
550 
551  Cmp temp(E_I) ;
552  Tenseur E_i (temp) ;
553  Tenseur J_i (Cmp(J_I(3))) ;
554 
555  Tenseur source = nnn * (ener_euler + E_em + E_i
556  + s_euler + Spp_em + SrrplusStt + Spp) +
557  nphi * (Jp_em + J_i)
558  + 2 * bbb * (ener_euler + press) * tnphi * uuu ;
559 
560  source = a_car * bbb * source ;
561 
562  source.set_std_base() ;
563 
564  p_mass_g = new double( source().integrale() ) ;
565 
566 
567  }
568  else{ // Newtonian case
569  p_mass_g = new double( mass_b() ) ; // in the Newtonian case
570  // M_g = M_b
571  }
572  }
573 
574  return *p_mass_g ;
575 
576 }
577 
578  //----------------------------//
579  // Angular momentum //
580  //----------------------------//
581 
583 
584  if (p_angu_mom == 0x0) { // a new computation is required
585 
586  Cmp dens = uuu() ;
587 
588  dens.mult_r() ; // Multiplication by
589  dens.va = (dens.va).mult_st() ; // r sin(theta)
590 
591  if (relativistic) {
592  dens = a_car() * (b_car() * (ener_euler() + press())
593  * dens + bbb() * (Jp_em() + Cmp(J_I(3)) ) ) ;
594  }
595  else { // Newtonian case
596  dens = nbar() * dens ;
597  }
598 
599  dens.std_base_scal() ;
600 
601  p_angu_mom = new double( dens.integrale() ) ;
602 
603  }
604 
605  return *p_angu_mom ;
606 
607 }
608 
609  //----------------------------//
610  // GRV2 //
611  //----------------------------//
612 
613 double Et_magnetisation::grv2() const {
614 
615  if (p_grv2 == 0x0) { // a new computation is required
616 
617  // To get qpig:
618  using namespace Unites ;
619 
620  Tenseur Spp (Cmp(Sij_I(3, 3))) ; //S_{\phi\phi}
621  Spp = Spp / b_car ; // S^\phi_\phi
622 
623  Tenseur sou_m = 2 * qpig * a_car * (press + (ener_euler+press)
624  * uuu*uuu + Spp) ;
625 
626  Tenseur sou_q = 2 * qpig * a_car * Spp_em + 1.5 * ak_car
627  - flat_scalar_prod(logn.gradient_spher(), logn.gradient_spher() ) ;
628 
629  p_grv2 = new double( double(1) - lambda_grv2(sou_m(), sou_q()) ) ;
630 
631  }
632 
633  return *p_grv2 ;
634 
635 }
636 
637 
638  //----------------------------//
639  // GRV3 //
640  //----------------------------//
641 
642 double Et_magnetisation::grv3(ostream* ost) const {
643 
644  if (p_grv3 == 0x0) { // a new computation is required
645 
646  // To get qpig:
647  using namespace Unites ;
648 
649  Tenseur source(mp) ;
650 
651  // Gravitational term [cf. Eq. (43) of Gourgoulhon & Bonazzola
652  // ------------------ Class. Quantum Grav. 11, 443 (1994)]
653 
654  if (relativistic) {
655  Tenseur alpha = dzeta - logn ;
656  Tenseur beta = log( bbb ) ;
657  beta.set_std_base() ;
658 
659  source = 0.75 * ak_car
660  - flat_scalar_prod(logn.gradient_spher(),
661  logn.gradient_spher() )
662  + 0.5 * flat_scalar_prod(alpha.gradient_spher(),
663  beta.gradient_spher() ) ;
664 
665  Cmp aa = alpha() - 0.5 * beta() ;
666  Cmp daadt = aa.srdsdt() ; // 1/r d/dth
667 
668  // What follows is valid only for a mapping of class Map_radial :
669  const Map_radial* mpr = dynamic_cast<const Map_radial*>(&mp) ;
670  if (mpr == 0x0) {
671  cout << "Etoile_rot::grv3: the mapping does not belong"
672  << " to the class Map_radial !" << endl ;
673  abort() ;
674  }
675 
676  // Computation of 1/tan(theta) * 1/r daa/dtheta
677  if (daadt.get_etat() == ETATQCQ) {
678  Valeur& vdaadt = daadt.va ;
679  vdaadt = vdaadt.ssint() ; // division by sin(theta)
680  vdaadt = vdaadt.mult_ct() ; // multiplication by cos(theta)
681  }
682 
683  Cmp temp = aa.dsdr() + daadt ;
684  temp = ( bbb() - a_car()/bbb() ) * temp ;
685  temp.std_base_scal() ;
686 
687  // Division by r
688  Valeur& vtemp = temp.va ;
689  vtemp = vtemp.sx() ; // division by xi in the nucleus
690  // Id in the shells
691  // division by xi-1 in the ZEC
692  vtemp = (mpr->xsr) * vtemp ; // multiplication by xi/r in the nucleus
693  // by 1/r in the shells
694  // by r(xi-1) in the ZEC
695 
696  // In the ZEC, a multiplication by r has been performed instead
697  // of the division:
698  temp.set_dzpuis( temp.get_dzpuis() + 2 ) ;
699 
700  source = bbb() * source() + 0.5 * temp ;
701 
702  }
703  else{
704  source = - 0.5 * flat_scalar_prod(logn.gradient_spher(),
705  logn.gradient_spher() ) ;
706  }
707 
708  source.set_std_base() ;
709 
710  double int_grav = source().integrale() ;
711 
712  // Matter term
713  // -----------
714 
715  if (relativistic) {
716 
717  // S_{rr} + S_{\theta\theta}
718  Tenseur SrrplusStt( Cmp(Sij_I(1, 1) + Sij_I(2, 2)) ) ;
719  SrrplusStt = SrrplusStt / a_car ; // S^r_r + S^\theta_\theta
720 
721  Tenseur Spp (Cmp(Sij_I(3, 3))) ; //S_{\phi\phi}
722  Spp = Spp / b_car ; // S^\phi_\phi
723 
724  source = qpig * a_car * bbb * ( s_euler + Spp_em + SrrplusStt + Spp ) ;
725  }
726  else{
727  source = qpig * ( 3 * press + nbar * uuu * uuu ) ;
728  }
729 
730  source.set_std_base() ;
731 
732  double int_mat = source().integrale() ;
733 
734  // Virial error
735  // ------------
736  if (ost != 0x0) {
737  *ost << "Et_magnetisation::grv3 : gravitational term : " << int_grav
738  << endl ;
739  *ost << "Et_magnetisation::grv3 : matter term : " << int_mat
740  << endl ;
741  }
742 
743  p_grv3 = new double( (int_grav + int_mat) / int_mat ) ;
744 
745  }
746 
747  return *p_grv3 ;
748 
749 }
750 
751  //----------------------------//
752  // Quadrupole moment //
753  //----------------------------//
754 
756 
757  if (p_mom_quad_old == 0x0) { // a new computation is required
758 
759  // To get qpig:
760  using namespace Unites ;
761 
762  // Source for of the Poisson equation for nu
763  // -----------------------------------------
764 
765  Tenseur source(mp) ;
766 
767  if (relativistic) {
768  // S_{rr} + S_{\theta\theta}
769  Tenseur SrrplusStt( Cmp(Sij_I(1, 1) + Sij_I(2, 2)) ) ;
770  SrrplusStt = SrrplusStt / a_car ; // S^r_r + S^\theta_\theta
771 
772  Tenseur Spp (Cmp(Sij_I(3, 3))) ; //S_{\phi\phi}
773  Spp = Spp / b_car ; // S^\phi_\phi
774 
775  Cmp temp(E_I) ;
776  Tenseur E_i(temp) ;
777 
778  Tenseur beta = log(bbb) ;
779  beta.set_std_base() ;
780  source = qpig * a_car *( ener_euler + E_em + E_i
781  + s_euler + Spp_em + SrrplusStt + Spp)
782  + ak_car - flat_scalar_prod(logn.gradient_spher(),
783  logn.gradient_spher() + beta.gradient_spher()) ;
784  }
785  else {
786  source = qpig * nbar ;
787  }
788  source.set_std_base() ;
789 
790  // Multiplication by -r^2 P_2(cos(theta))
791  // [cf Eq.(7) of Salgado et al. Astron. Astrophys. 291, 155 (1994) ]
792  // ------------------------------------------------------------------
793 
794  // Multiplication by r^2 :
795  // ----------------------
796  Cmp& csource = source.set() ;
797  csource.mult_r() ;
798  csource.mult_r() ;
799  if (csource.check_dzpuis(2)) {
800  csource.inc2_dzpuis() ;
801  }
802 
803  // Muliplication by cos^2(theta) :
804  // -----------------------------
805  Cmp temp = csource ;
806 
807  // What follows is valid only for a mapping of class Map_radial :
808  assert( dynamic_cast<const Map_radial*>(&mp) != 0x0 ) ;
809 
810  if (temp.get_etat() == ETATQCQ) {
811  Valeur& vtemp = temp.va ;
812  vtemp = vtemp.mult_ct() ; // multiplication by cos(theta)
813  vtemp = vtemp.mult_ct() ; // multiplication by cos(theta)
814  }
815 
816  // Muliplication by -P_2(cos(theta)) :
817  // ----------------------------------
818  source = 0.5 * source() - 1.5 * temp ;
819 
820  // Final result
821  // ------------
822  p_mom_quad_old = new double( - source().integrale() / qpig ) ;
823  }
824  return *p_mom_quad_old ;
825  }
826 
827 
829 
830  using namespace Unites ;
831 
832  if (p_mom_quad_Bo == 0x0) { // a new computation is required
833 
834  // S_{rr} + S_{\theta\theta} = A^2*(S^r_r + S^\theta_\theta)
835  Tenseur SrrplusStt( Cmp(Sij_I(1, 1) + Sij_I(2, 2)) ) ;
836 
837  Cmp dens = a_car() * press() ;
838  dens = bbb() * nnn() * (SrrplusStt() + 2*dens) ;
839  dens.mult_rsint() ;
840  dens.std_base_scal() ;
841 
842  p_mom_quad_Bo = new double( - 16. * dens.integrale() / qpig ) ;
843  }
844  return *p_mom_quad_Bo ;
845  }
846 
847 
848 
849 }
const Cmp & dsdr() const
Returns of *this .
Definition: cmp_deriv.C:84
Metric for tensor calculation.
Definition: metric.h:90
virtual double mom_quad_Bo() const
Part of the quadrupole moment.
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:296
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
const Tenseur & gradient_spher() const
Returns the gradient of *this (Spherical coordinates) (scalar field only).
Definition: tenseur.C:1548
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
void dec_dzpuis()
Decreases by 1 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
Definition: cmp_r_manip.C:154
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:348
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1170
Multi-domain array.
Definition: mtbl.h:118
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
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 .
Definition: cmp_deriv.C:105
const Valeur & sx() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR ) ...
Definition: valeur_sx.C:110
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.
Definition: valeur.h:287
Tensor field of valence 1.
Definition: vector.h:188
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:94
void div_r()
Division by r everywhere.
Definition: cmp_r_manip.C:78
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
virtual double angu_mom() const
Angular momentum.
void mult_r()
Multiplication by r everywhere.
Definition: cmp_r_manip.C:91
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition: tenseur.C:668
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.
Definition: valeur_ssint.C:112
double * t
The array of double.
Definition: tbl.h:173
Parameter storage.
Definition: param.h:125
Base class for pure radial mappings.
Definition: map.h:1536
void mult_rsint()
Multiplication by .
Definition: cmp_r_manip.C:116
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 .
Definition: cmp_integ.C:55
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...
Definition: valeur.C:900
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain...
Definition: map.h:1549
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
void inc2_dzpuis()
Increases by 2 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
Definition: cmp_r_manip.C:192
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition: cmp.C:644
virtual double mom_quad_old() const
Part of the quadrupole moment.
Cmp poisson() const
Solves the scalar Poisson equation with *this as a source.
Definition: cmp_pde.C:94
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...
Definition: cmp.C:323
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:724
int get_dzpuis() const
Returns dzpuis.
Definition: cmp.h:903
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:654
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...
Definition: cmp.C:715
virtual double mass_g() const
Gravitational mass.
Valeur ** asymptot(int n, const int flag=0) const
Asymptotic expansion at r = infinity.
Definition: cmp_asymptot.C:71
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.
Definition: cmp.C:654
Basic array class.
Definition: tbl.h:161
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:636
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.
Definition: cmp_deriv.C:242
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:464
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:223
Standard electro-magnetic units.
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:301
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.