LORENE
et_rot_mag_mag.C
1/*
2 * Computes magnetic fields and derived quantities for rotating equilibrium
3 *
4 * (see file et_rot_mag.h for documentation)
5 *
6 */
7
8/*
9 * Copyright (c) 2002 Emmanuel Marcq
10 * Copyright (c) 2002 Jerome Novak
11 *
12 * This file is part of LORENE.
13 *
14 * LORENE is free software; you can redistribute it and/or modify
15 * it under the terms of the GNU General Public License as published by
16 * the Free Software Foundation; either version 2 of the License, or
17 * (at your option) any later version.
18 *
19 * LORENE is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License
25 * along with LORENE; if not, write to the Free Software
26 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27 *
28 */
29
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 $" ;
31
32/*
33 * $Id: et_rot_mag_mag.C,v 1.17 2014/10/13 08:52:58 j_novak Exp $
34 * $Log: et_rot_mag_mag.C,v $
35 * Revision 1.17 2014/10/13 08:52:58 j_novak
36 * Lorene classes and functions now belong to the namespace Lorene.
37 *
38 * Revision 1.16 2014/09/03 15:33:42 j_novak
39 * Filtering of Maxwell sources is now optional.
40 *
41 * Revision 1.15 2014/07/04 12:15:12 j_novak
42 * Added filtering.
43 *
44 * Revision 1.14 2005/06/03 15:31:56 j_novak
45 * Better computation when more than one point in phi.
46 *
47 * Revision 1.13 2003/10/03 15:58:47 j_novak
48 * Cleaning of some headers
49 *
50 * Revision 1.12 2002/09/09 13:00:39 e_gourgoulhon
51 * Modification of declaration of Fortran 77 prototypes for
52 * a better portability (in particular on IBM AIX systems):
53 * All Fortran subroutine names are now written F77_* and are
54 * defined in the new file C++/Include/proto_f77.h.
55 *
56 * Revision 1.11 2002/06/05 15:15:59 j_novak
57 * The case of non-adapted mapping is treated.
58 * parmag.d and parrot.d have been merged.
59 *
60 * Revision 1.10 2002/06/03 13:23:16 j_novak
61 * The case when the mapping is not adapted is now treated
62 *
63 * Revision 1.9 2002/06/03 13:00:45 e_marcq
64 *
65 * conduc parameter read in parmag.d
66 *
67 * Revision 1.7 2002/05/20 10:31:59 j_novak
68 * *** empty log message ***
69 *
70 * Revision 1.6 2002/05/17 15:08:01 e_marcq
71 *
72 * Rotation progressive plug-in, units corrected, Q and a_j new member data
73 *
74 * Revision 1.5 2002/05/16 10:02:09 j_novak
75 * Errors in stress energy tensor corrected
76 *
77 * Revision 1.4 2002/05/15 09:54:00 j_novak
78 * First operational version
79 *
80 * Revision 1.3 2002/05/14 13:38:36 e_marcq
81 *
82 *
83 * Unit update, new outputs
84 *
85 * Revision 1.1 2002/05/10 09:26:52 j_novak
86 * Added new class Et_rot_mag for magnetized rotating neutron stars (under development)
87 *
88 *
89 * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_mag_mag.C,v 1.17 2014/10/13 08:52:58 j_novak Exp $
90 *
91 */
92
93// Headers C
94#include <cstdlib>
95#include <cmath>
96
97// Headers Lorene
98#include "et_rot_mag.h"
99#include "utilitaires.h"
100#include "param.h"
101#include "proto_f77.h"
102#include "graphique.h"
103#include "tensor.h"
104
105namespace Lorene {
106// Local prototype (for drawings only)
107Cmp raccord_c1(const Cmp& uu, int l1) ;
108
109// Algo du papier de 1995
110
111void Et_rot_mag::magnet_comput(const int adapt_flag,
112 Cmp (*f_j)(const Cmp&, const double),
113 Param& par_poisson_At,
114 Param& par_poisson_Avect){
115 double relax_mag = 0.5 ;
116
117 int mag_filter = 0 ;
118 if (par_poisson_At.get_n_int() > 1)
119 mag_filter = par_poisson_At.get_int(1) ;
120
121 int Z = mp.get_mg()->get_nzone();
122
123 if(is_conduct()) {
124 bool adapt(adapt_flag) ;
125 /****************************************************************
126 * Assertion that all zones have same number of points in theta
127 ****************************************************************/
128 int nt = mp.get_mg()->get_nt(nzet-1) ;
129 for (int l=0; l<Z; l++) assert(mp.get_mg()->get_nt(l) == nt) ;
130
131 Tbl Rsurf(nt) ;
132 Rsurf.set_etat_qcq() ;
133 mp.r.fait() ;
134 mp.tet.fait() ;
135 Mtbl* theta = mp.tet.c ;
136 const Map_radial* mpr = dynamic_cast<const Map_radial*>(&mp) ;
137 assert (mpr != 0x0) ;
138 for (int j=0; j<nt; j++)
139 Rsurf.set(j) = mpr->val_r_jk(l_surf()(0,j), xi_surf()(0,j), j, 0) ;
140
141
142 // Calcul de A_0t dans l'etoile (conducteur parfait)
143
144 Cmp A_0t(- omega * A_phi) ;
145 //A_0t.annule(nzet,Z-1) ;
146
147 Tenseur ATTENS(A_t) ;
148 Tenseur APTENS(A_phi) ;
149 Tenseur BMN(-logn) ;
150 BMN = BMN + log(bbb) ;
151 BMN.set_std_base() ;
152
153
155 nphi.gradient_spher())());
157 nphi.gradient_spher())()) ;
159 BMN.gradient_spher())()
161 BMN.gradient_spher())()) ;
162
163 Cmp ATANT(A_phi.srdsdt()); // Constrction par copie pour mapping
164
165 ATANT.va = ATANT.va.mult_ct().ssint() ;
166
167 Cmp ttnphi(tnphi()) ;
168 ttnphi.mult_rsint() ;
169 Cmp BLAH(- b_car()/(nnn()*nnn())*ttnphi*grad1) ;
170 BLAH -= (1+b_car()/(nnn()*nnn())*tnphi()*tnphi())*grad2 ;
171 Cmp nphisr(nphi()) ;
172 nphisr.div_r() ;
173 Cmp npgrada(2*nphisr*(A_phi.dsdr()+ATANT )) ;
174 npgrada.inc2_dzpuis() ;
175 BLAH -= grad3 + npgrada ;
176 Cmp gtt(-nnn()*nnn()+b_car()*tnphi()*tnphi()) ;
177 Cmp gtphi( - b_car()*ttnphi) ;
178
179 // Calcul de j_t grace a Maxwell-Gauss
180 Cmp tmp(((BLAH - A_0t.laplacien())/a_car() - gtphi*j_phi)
181 / gtt);
182 tmp.annule(nzet, Z-1) ;
183 if (adapt) {
184 j_t = tmp ;
185 }
186 else {
187 j_t.allocate_all() ;
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) ) ;
193 j_t.annule(nzet,Z-1) ;
194 }
196
197 // Calcul du courant j_phi
198 j_phi = omega * j_t + (ener() + press())*f_j(A_phi, a_j) ;
200
201 // Resolution de Maxwell Ampere (-> A_phi)
202 // Calcul des termes sources avec A-t du pas precedent.
203
205 BMN.gradient_spher())());
206
207 Tenseur source_tAphi(mp, 1, CON, mp.get_bvect_spher()) ;
208
209 source_tAphi.set_etat_qcq() ;
210 Cmp tjphi(j_phi) ;
211 tjphi.mult_rsint() ;
212 Cmp tgrad1(grad1) ;
213 tgrad1.mult_rsint() ;
214 Cmp d_grad4(grad4) ;
215 d_grad4.div_rsint() ;
216 source_tAphi.set(0)=0 ;
217 source_tAphi.set(1)=0 ;
218 source_tAphi.set(2)= -b_car()*a_car()*(tjphi-tnphi()*j_t)
219 + b_car()/(nnn()*nnn())*(tgrad1+tnphi()*grad2)+d_grad4 ;
220
221 source_tAphi.change_triad(mp.get_bvect_cart());
222 if (mag_filter == 1) {
223 for (int i=0; i<3; i++) {
224 Scalar tmp_filter = source_tAphi(i) ;
225 tmp_filter.exponential_filter_r(0, 2, 1) ;
226 tmp_filter.exponential_filter_ylm(0, 2, 1) ;
227 source_tAphi.set(i) = tmp_filter ;
228 }
229 }
230
231
232 Tenseur WORK_VECT(mp, 1, CON, mp.get_bvect_cart()) ;
233 WORK_VECT.set_etat_qcq() ;
234 for (int i=0; i<3; i++) {
235 WORK_VECT.set(i) = 0 ;
236 }
237 Tenseur WORK_SCAL(mp) ;
238 WORK_SCAL.set_etat_qcq() ;
239 WORK_SCAL.set() = 0 ;
240
241 double lambda_mag = 0. ; // No 3D version !
242
243 Tenseur AVECT(source_tAphi) ;
244 if (source_tAphi.get_etat() != ETATZERO) {
245
246 for (int i=0; i<3; i++) {
247 if(source_tAphi(i).dz_nonzero()) {
248 assert( source_tAphi(i).get_dzpuis() == 4 ) ;
249 }
250 else{
251 (source_tAphi.set(i)).set_dzpuis(4) ;
252 }
253 }
254
255 }
256 source_tAphi.poisson_vect(lambda_mag, par_poisson_Avect, AVECT, WORK_VECT,
257 WORK_SCAL) ;
259 Cmp A_phi_n(AVECT(2));
260 A_phi_n.mult_rsint() ;
261
262 // Resolution de Maxwell-Ampere : A_1
263
264 Cmp source_A_1t(-a_car()*(j_t*gtt + j_phi*gtphi) + BLAH);
265 if (mag_filter == 1) {
266 Scalar tmp_filter = source_A_1t ;
267 tmp_filter.exponential_filter_r(0, 2, 1) ;
268 tmp_filter.exponential_filter_ylm(0, 2, 1) ;
269 source_A_1t = tmp_filter ;
270 }
271
272 Cmp A_1t(mp);
273 A_1t = 0 ;
274 source_A_1t.poisson(par_poisson_At, A_1t) ;
275
276 int L = mp.get_mg()->get_nt(0) ;
277
278 Tbl MAT(L,L) ;
279 Tbl MAT_PHI(L,L);
280 Tbl VEC(L) ;
281
282 MAT.set_etat_qcq() ;
283 VEC.set_etat_qcq() ;
284 MAT_PHI.set_etat_qcq() ;
285
286 Tbl leg(L,2*L) ;
287 leg.set_etat_qcq() ;
288
289 Cmp psi(mp);
290 Cmp psi2(mp);
291 Cmp psi3(mp);
292 psi.allocate_all() ;
293 psi2.allocate_all() ;
294 psi3.allocate_all() ;
295
296 Tbl VEC3(L) ;
297 VEC3.set_etat_qcq() ;
298 for (int i=0; i<L; i++)
299 VEC3.set(i) = 1. / double(i+1) ;
300
301 for (int p=0; p<mp.get_mg()->get_np(0); p++) {
302 // leg[k,l] : legendre_l(cos(theta_k))
303 // Construction par recurrence de degre 2
304 for(int k=0;k<L;k++){
305 for(int l=0;l<2*L;l++){
306
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)
310 * cos((*theta)(l_surf()(p,k),p,k,0))
311 * leg(k,l-1)-double(l-1)/double(l)*leg(k,l-2) ;
312 }
313 }
314
315 for(int k=0;k<L;k++){
316
317 // Valeurs a la surface trouvees via va.val_point_jk(l,xisurf,k,p)
318
319 VEC.set(k) = A_0t.va.val_point_jk(l_surf()(p,k), xi_surf()(p,k), k, p)
320 -A_1t.va.val_point_jk(l_surf()(p,k), xi_surf()(p,k), k, p);
321
322 for(int l=0;l<L;l++) MAT.set(l,k) = leg(k,2*l)/pow(Rsurf(k),2*l+1);
323
324 }
325 // appel fortran :
326
327 int* IPIV=new int[L] ;
328 int INFO ;
329
330 Tbl MAT_SAVE(MAT) ;
331 Tbl VEC2(L) ;
332 VEC2.set_etat_qcq() ;
333 int un = 1 ;
334
335 F77_dgesv(&L, &un, MAT.t, &L, IPIV, VEC.t, &L, &INFO) ;
336
337 // coeffs a_l dans VEC
338
339 for(int k=0;k<L;k++) {VEC2.set(k)=1. ; }
340
341 F77_dgesv(&L, &un, MAT_SAVE.t, &L, IPIV, VEC2.t, &L, &INFO) ;
342
343 delete [] IPIV ;
344
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)) ;
358 }
359 }
360 }
361 }
362 }
363 psi.std_base_scal() ;
364 psi2.std_base_scal() ;
365
366 assert(psi.get_dzpuis() == 0) ;
367 int dif = A_1t.get_dzpuis() ;
368 if (dif > 0) {
369 for (int d=0; d<dif; d++) A_1t.dec_dzpuis() ;
370 }
371
372 if (adapt) {
373 Cmp A_t_ext(A_1t + psi) ;
374 A_t_ext.annule(0,nzet-1) ;
375 A_0t += A_t_ext ;
376 }
377 else {
378 tmp = A_0t ;
379 A_0t.allocate_all() ;
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) ) ;
385 }
386 A_0t.std_base_scal() ;
387
388 if (mag_filter == 1) {
389 Scalar tmp_filter = A_0t ;
390 tmp_filter.exponential_filter_r(0, 2, 1) ;
391 tmp_filter.exponential_filter_ylm(0, 2, 1) ;
392 A_0t = tmp_filter ;
393 }
394
395 Valeur** asymp = A_0t.asymptot(1) ;
396
397 double Q_0 = -4*M_PI*(*asymp[1])(Z-1,0,0,0) ; // utilise A_0t plutot que E
398 delete asymp[0] ;
399 delete asymp[1] ;
400
401 delete [] asymp ;
402
403 asymp = psi2.asymptot(1) ;
404
405 double Q_2 = -4*M_PI*(*asymp[1])(Z-1,0,0,0) ; // A_2t = psi2 a l'infini
406 delete asymp[0] ;
407 delete asymp[1] ;
408
409 delete [] asymp ;
410
411 // solution definitive de A_t:
412
413 double C = (Q-Q_0)/Q_2 ;
414
415 assert(psi2.get_dzpuis() == 0) ;
416 dif = A_0t.get_dzpuis() ;
417 if (dif > 0) {
418 for (int d=0; d<dif; d++) A_0t.dec_dzpuis() ;
419 }
420 Cmp A_t_n(mp) ;
421 if (adapt) {
422 A_t_n = A_0t + C ;
423 Cmp A_t_ext(A_0t + C*psi2) ;
424 A_t_ext.annule(0,nzet-1) ;
425 A_t_n.annule(nzet,Z-1) ;
426 A_t_n += A_t_ext ;
427 }
428 else {
429 A_t_n.allocate_all() ;
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 ) ;
436 }
437 }
438 A_t_n.std_base_scal() ;
439 if (mag_filter == 1) {
440 Scalar tmp_filter = A_t_n ;
441 tmp_filter.exponential_filter_r(0, 2, 1) ;
442 tmp_filter.exponential_filter_ylm(0, 2, 1) ;
443 A_t_n = tmp_filter ;
444 }
445
446 asymp = A_t_n.asymptot(1) ;
447
448 delete asymp[0] ;
449 delete asymp[1] ;
450
451 delete [] asymp ;
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 ;
454
455 } // End of perfect conductor case
456
457 else
458 {
459
460 /***************
461 * CAS ISOLANT *
462 ***************/
463
464 // Calcul de j_t
465 j_t = Q*nbar() + (ener()+press())*f_j(omega* A_phi - A_t,a_j) ;
466 j_t.annule(nzet,Z-1) ;
468
469 // Calcul de j_phi
470 j_phi = omega * j_t ;
472
473 // Resolution de A_t
474
475 Tenseur ATTENS(A_t) ;
476 Tenseur APTENS(A_phi) ;
477 Tenseur BMN(-logn) ;
478 BMN = BMN + log(bbb) ;
479 BMN.set_std_base() ;
480
481
483 nphi.gradient_spher())());
485 nphi.gradient_spher())()) ;
487 BMN.gradient_spher())()
489 BMN.gradient_spher())()) ;
490
491 Cmp ATANT(A_phi.srdsdt());
492
493 ATANT.va = ATANT.va.mult_ct().ssint() ;
494
495 Cmp ttnphi(tnphi()) ;
496 ttnphi.mult_rsint() ;
497 Cmp BLAH(- b_car()/(nnn()*nnn())*ttnphi*grad1) ;
498 BLAH -= (1+b_car()/(nnn()*nnn())*tnphi()*tnphi())*grad2 ;
499 Cmp nphisr(nphi()) ;
500 nphisr.div_r() ;
501 Cmp npgrada(2*nphisr*(A_phi.dsdr()+ATANT )) ;
502 npgrada.inc2_dzpuis() ;
503 BLAH -= grad3 + npgrada ;
504 Cmp gtt(-nnn()*nnn()+b_car()*tnphi()*tnphi()) ;
505 Cmp gtphi( - b_car()*ttnphi) ;
506
507 Cmp source_A_t_n(mp);
508 if (relativistic) {
509 source_A_t_n = (-a_car()*(j_t*gtt + j_phi*gtphi) + BLAH);
510 source_A_t_n.std_base_scal();}
511 else{
512 source_A_t_n = j_t;}
513
514 Cmp A_t_n(A_t) ;
515 A_t_n = 0 ;
516 A_t_n.std_base_scal() ;
517
518 source_A_t_n.poisson(par_poisson_At, A_t_n) ;
519
520 // Resolution de A_phi
521
523 BMN.gradient_spher())());
524
525 Tenseur source_tAphi(mp, 1, CON, mp.get_bvect_spher()) ;
526
527 source_tAphi.set_etat_qcq() ;
528
529 Cmp tjphi(j_phi) ;
530 tjphi.mult_rsint() ;
531 Cmp tgrad1(grad1) ;
532 tgrad1.mult_rsint() ;
533 Cmp d_grad4(grad4) ;
534 d_grad4.div_rsint() ;
535 source_tAphi.set(0)=0 ;
536 source_tAphi.set(1)=0 ;
537
538 if (relativistic) {
539 source_tAphi.set(2)= -b_car()*a_car()*(tjphi-tnphi()*j_t)
540 + b_car()/(nnn()*nnn())*(tgrad1+tnphi()*grad2)+d_grad4 ;}
541 else{
542 source_tAphi.set(2)= - tjphi ;}
543
544 source_tAphi.change_triad(mp.get_bvect_cart());
545
546 Tenseur WORK_VECT(mp, 1, CON, mp.get_bvect_cart()) ;
547 WORK_VECT.set_etat_qcq() ;
548 for (int i=0; i<3; i++) {
549 WORK_VECT.set(i) = 0 ;
550 }
551 Tenseur WORK_SCAL(mp) ;
552 WORK_SCAL.set_etat_qcq() ;
553 WORK_SCAL.set() = 0 ;
554
555 double lambda_mag = 0. ; // No 3D version !
556
557 Tenseur AVECT(source_tAphi) ;
558 if (source_tAphi.get_etat() != ETATZERO) {
559
560 for (int i=0; i<3; i++) {
561 if(source_tAphi(i).dz_nonzero()) {
562 assert( source_tAphi(i).get_dzpuis() == 4 ) ;
563 }
564 else{
565 (source_tAphi.set(i)).set_dzpuis(4) ;
566 }
567 }
568
569 }
570 source_tAphi.poisson_vect(lambda_mag, par_poisson_Avect, AVECT, WORK_VECT,
571 WORK_SCAL) ;
573 Cmp A_phi_n(AVECT(2));
574 A_phi_n.mult_rsint() ;
575
576 // Relaxation
577
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 ;
580
581 }
582
583
584}
585
586
587
588
589
590
591
592
593
594
595
596
597
598}
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...
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 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
const Cmp & srdsdt() const
Returns of *this .
Definition cmp_deriv.C:105
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
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
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 .
Definition et_rot_mag.h:150
bool is_conduct() const
Tells if the star is made of conducting or isolating material.
Definition et_rot_mag.h:241
Cmp j_t
t-component of the current 4-vector
Definition et_rot_mag.h:158
double Q
In the case of a perfect conductor, the requated baryonic charge.
Definition et_rot_mag.h:179
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...
Tenseur nphi
Metric coefficient .
Definition etoile.h:1510
Tenseur bbb
Metric factor B.
Definition etoile.h:1504
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
Tenseur nnn
Total lapse function.
Definition etoile.h:509
Tenseur nbar
Baryon density in the fluid frame.
Definition etoile.h:459
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
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
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
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
Multi-domain array.
Definition mtbl.h:118
Parameter storage.
Definition param.h:125
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition param.C:292
int get_n_int() const
Returns the number of stored int 's addresses.
Definition param.C:239
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.
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.
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 & ssint() const
Returns of *this.
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
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