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
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 $" ;
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
91namespace Lorene {
92
93 using namespace Unites_mag ;
94
95// Algo du papier de 1995
96
97void 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())()
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
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 }
206
207 // Calcul du courant j_phi
208 j_phi = omega * j_t + (ener() + press())*f_j(A_phi, a_j) ;
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) ;
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) ;
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 ;
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
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
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
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
642double 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
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)
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}
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
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...
Definition cmp.C:323
void div_r()
Division by r everywhere.
Definition cmp_r_manip.C:78
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.
Definition cmp.h:899
Valeur va
The numerical value of the Cmp
Definition cmp.h:464
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition cmp.C:644
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:348
const Cmp & laplacien(int zec_mult_r=4) const
Returns the Laplacian of *this.
Definition cmp_deriv.C:242
int get_dzpuis() const
Returns dzpuis.
Definition cmp.h:903
void mult_r()
Multiplication by r everywhere.
Definition cmp_r_manip.C:91
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.
Definition cmp.h:724
void set_dzpuis(int)
Set a value to dzpuis.
Definition cmp.C:654
double integrale() const
Computes the integral over all space of *this .
Definition cmp_integ.C:55
const Cmp & srdsdt() const
Returns of *this .
Definition cmp_deriv.C:105
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
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.
Definition cmp_pde.C:94
const Cmp & dsdr() const
Returns of *this .
Definition cmp_deriv.C:84
void fait() const
Computes, at each point of the grid, the value of the coordinate or mapping derivative represented by...
Definition coord.C:116
Mtbl * c
The coordinate values at each grid point.
Definition coord.h:97
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.
Definition et_rot_mag.h:630
Vector J_I
Interaction momentum density 3-vector.
Definition et_rot_mag.h:627
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.
Definition et_rot_mag.h:677
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.
Definition et_rot_mag.h:624
virtual double mass_g() const
Gravitational mass.
Tenseur Srr_em
rr component of the electromagnetic stress 3-tensor, as measured in the Eulerian frame....
Definition et_rot_mag.h:170
Cmp j_phi
-component of the current 4-vector
Definition et_rot_mag.h:159
Cmp A_phi
-component of the electromagnetic potential 1-form divided by .
Definition et_rot_mag.h:155
double a_j
Amplitude of the curent/charge function.
Definition et_rot_mag.h:180
Tenseur Spp_em
component of the electromagnetic stress 3-tensor, as measured in the Eulerian frame.
Definition et_rot_mag.h:173
Cmp A_t
t-component of the elecctromagnetic potential 1-form, divided by .
Definition et_rot_mag.h:150
Cmp j_t
t-component of the current 4-vector
Definition et_rot_mag.h:158
Tenseur E_em
electromagnetic energy density in the Eulerian frame
Definition et_rot_mag.h:161
Tenseur Jp_em
component of the electromagnetic momentum density 3-vector, as measured in the Eulerian frame.
Definition et_rot_mag.h:167
double Q
In the case of a perfect conductor, the requated baryonic charge.
Definition et_rot_mag.h:179
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.
Definition etoile.h:1518
double omega
Rotation angular velocity ([f_unit] )
Definition etoile.h:1501
Tenseur & logn
Metric potential = logn_auto.
Definition etoile.h:1521
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.
Definition etoile.h:1642
Tenseur nphi
Metric coefficient .
Definition etoile.h:1510
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.
Definition etoile.h:1504
Tenseur ak_car
Scalar .
Definition etoile.h:1586
Tenseur & dzeta
Metric potential = beta_auto.
Definition etoile.h:1534
double * p_grv3
Error on the virial identity GRV3.
Definition etoile.h:1634
double * p_grv2
Error on the virial identity GRV2.
Definition etoile.h:1633
double * p_mom_quad_Bo
Part of the quadrupole moment.
Definition etoile.h:1643
double * p_angu_mom
Angular momentum.
Definition etoile.h:1631
Tenseur b_car
Square of the metric factor B.
Definition etoile.h:1507
Tenseur tnphi
Component of the shift vector.
Definition etoile.h:1515
int nzet
Number of domains of *mp occupied by the star.
Definition etoile.h:432
double * p_mass_g
Gravitational mass.
Definition etoile.h:548
Tenseur nnn
Total lapse function.
Definition etoile.h:509
Tenseur nbar
Baryon density in the fluid frame.
Definition etoile.h:459
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition etoile.h:474
Tenseur gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition etoile.h:471
Map & mp
Mapping associated with the star.
Definition etoile.h:429
Tenseur ener
Total energy density in the fluid frame.
Definition etoile.h:460
Tenseur press
Fluid pressure.
Definition etoile.h:461
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one.
Definition etoile.h:437
Tenseur ener_euler
Total energy density in the Eulerian frame.
Definition etoile.h:465
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
Definition etoile.h:468
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 .
Definition etoile.h:515
Base class for pure radial mappings.
Definition map.h:1536
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1549
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,...
Definition map.h:791
Coord r
r coordinate centered on the grid
Definition map.h:718
Coord tet
coordinate centered on the grid
Definition map.h:719
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Metric for tensor calculation.
Definition metric.h:90
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Multi-domain array.
Definition mtbl.h:118
Parameter storage.
Definition param.h:125
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
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.
Definition scalar.C:808
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.
Definition sym_tensor.h:223
Basic array class.
Definition tbl.h:161
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:361
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
double * t
The array of double.
Definition tbl.h:173
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tenseur.C:636
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition tenseur.C:1170
const Tenseur & gradient_spher() const
Returns the gradient of *this (Spherical coordinates) (scalar field only).
Definition tenseur.C:1548
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
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.
Definition tenseur.h:707
Values and coefficients of a (real-value) function.
Definition valeur.h:287
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...
Definition valeur.C:900
const Valeur & sx() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR )
Definition valeur_sx.C:110
const Valeur & ssint() const
Returns of *this.
Tensor field of valence 1.
Definition vector.h:188
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.
Definition vector.C:296
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:94
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
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).
Definition tensor.C:654
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.
Lorene prototypes.
Definition app_hor.h:64
Standard electro-magnetic units.
Standard units of space, time and mass.