LORENE
map_et_radius.C
1/*
2 * Methods of the class Map_et relative to the function
3 * r = R_l(xi, theta', phi')
4 */
5
6/*
7 * Copyright (c) 1999-2001 Eric Gourgoulhon
8 *
9 * This file is part of LORENE.
10 *
11 * LORENE is free software; you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License as published by
13 * the Free Software Foundation; either version 2 of the License, or
14 * (at your option) any later version.
15 *
16 * LORENE is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with LORENE; if not, write to the Free Software
23 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24 *
25 */
26
27
28char map_et_radius_C[] = "$Header: /cvsroot/Lorene/C++/Source/Map/map_et_radius.C,v 1.6 2014/10/13 08:53:05 j_novak Exp $" ;
29
30/*
31 * $Id: map_et_radius.C,v 1.6 2014/10/13 08:53:05 j_novak Exp $
32 * $Log: map_et_radius.C,v $
33 * Revision 1.6 2014/10/13 08:53:05 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.5 2013/06/05 15:10:42 j_novak
37 * Suppression of FINJAC sampling in r. This Jacobi(0,2) base is now
38 * available by setting colloc_r to BASE_JAC02 in the Mg3d constructor.
39 *
40 * Revision 1.4 2008/08/27 08:49:16 jl_cornou
41 * Added R_JACO02 case
42 *
43 * Revision 1.3 2004/01/26 16:58:35 j_novak
44 * Added initialization to avoid compiler warning.
45 *
46 * Revision 1.2 2003/12/19 16:21:43 j_novak
47 * Shadow hunt
48 *
49 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
50 * LORENE
51 *
52 * Revision 1.5 2000/10/20 13:14:47 eric
53 * Map_et::val_lx : la valeur par defaut de precis est desormais 1e-14
54 * (et non plus 1e-15).
55 *
56 * Revision 1.4 2000/01/04 13:05:47 eric
57 * \Corrections dans val_lx et val_lx_jk : initialisation de ftp/gtp
58 * par double(0) dans les cas ou il ne sont pas calcules.
59 *
60 * Revision 1.3 1999/12/17 11:03:36 eric
61 * Fonctions val_r_jk et val_lx_jk operationnelles.
62 *
63 * Revision 1.2 1999/12/17 09:29:22 eric
64 * val_lx : initialisation de niter a 0.
65 *
66 * Revision 1.1 1999/12/16 14:19:39 eric
67 * Initial revision
68 *
69 *
70 * $Header: /cvsroot/Lorene/C++/Source/Map/map_et_radius.C,v 1.6 2014/10/13 08:53:05 j_novak Exp $
71 *
72 */
73
74
75// Headers Lorene
76#include "map.h"
77#include "param.h"
78#include "nbr_spx.h"
79#include "utilitaires.h"
80
81// Local prototypes
82namespace Lorene {
83double fonc_invr_map_et_noyau(double, const Param&) ;
84double fonc_invr_map_et_coq(double, const Param&) ;
85double fonc_invr_map_et_zec(double, const Param&) ;
86
87 //------------------------------//
88 // val_r //
89 //------------------------------//
90
91
92double Map_et::val_r(int l, double xi, double theta, double pphi) const {
93
94 assert( l>=0 ) ;
95 assert( l<mg->get_nzone() ) ;
96
97 double resu ;
98 double ftp = ff.val_point(l, 0, theta, pphi) ; // value of F_l(theta,phi)
99
100 switch( mg->get_type_r(l) ) {
101
102 case RARE: {
103 double gtp = gg.val_point(l, 0, theta, pphi) ;
104 double xi_2 = xi * xi ;
105 double xi_3 = xi * xi_2 ;
106 double a = xi_2 * xi_2 * (3. - 2.*xi_2) ;
107 double b = ( 2.5 - 1.5 * xi_2 ) * xi_3 ;
108 resu = alpha[l] * ( xi + a * ftp + b * gtp ) + beta[l] ;
109 break ;
110 }
111
112 case FIN: {
113 double gtp = gg.val_point(l, 0, theta, pphi) ;
114 double xm1 = xi - 1. ;
115 double xp1 = xi + 1. ;
116 double a = 0.25* xm1 * xm1 * (xi + 2.) ;
117 double b = 0.25* xp1 * xp1 * (2. - xi) ;
118 resu = alpha[l] * ( xi + a * ftp + b * gtp ) + beta[l] ;
119 break ;
120 }
121
122 case UNSURR: {
123 double xm1 = xi - 1. ;
124 double a = 0.25* xm1 * xm1 * (xi + 2.) ;
125 resu = double(1) / ( alpha[l] * ( xi + a * ftp ) + beta[l] ) ;
126 break ;
127 }
128
129 default: {
130 cout << "Map_et::val_r: unknown type_r ! " << endl ;
131 abort () ;
132 }
133 }
134
135 return resu ;
136}
137
138 //------------------------------//
139 // val_lx //
140 //------------------------------//
141
142//-------------------------------
143// Version with default precision
144//-------------------------------
145
146void Map_et::val_lx(double rr, double theta, double pphi,
147 int& lz, double& xi) const {
148
149 int nitermax = 100 ; // Maximum number of iteration in the secant method
150 int niter ;
151 double precis = 1e-14 ; // Absolute precision in the secant method
152
153 Param par ;
154 par.add_int(nitermax) ;
155 par.add_int_mod(niter) ;
156 par.add_double(precis) ;
157
158 // Call of the version with precision parameters
159
160 val_lx(rr, theta, pphi, par, lz, xi) ;
161
162}
163
164
165//---------------------------------
166// Version with specified precision
167//---------------------------------
168
169void Map_et::val_lx(double rr, double theta, double pphi,
170 const Param& par, int& lz, double& xi) const {
171
172 int nz = mg->get_nzone() ;
173
174 // Precision in the secant method :
175 // ------------------------------
176
177 int nitermax = par.get_int() ;
178 int& niter = par.get_int_mod() ;
179 niter = 0 ; // initialisation
180 double precis = par.get_double() ;
181
182 // Particular case of r = infinity
183 //--------------------------------
184 if (rr == __infinity) {
185 if ( (mg->get_type_r(nz-1) == UNSURR) &&
186 (alpha[nz-1] == - beta[nz-1]) ) {
187 lz = nz - 1 ;
188 xi = 1 ;
189 }
190 else {
191 cout.precision(16);
192 cout.setf(ios::showpoint);
193 cout << "Map_et::val_lx: the domain containing r = " << rr <<
194 " has not been found ! "
195 << endl ;
196 abort () ;
197 }
198 return ;
199 }
200
201
202 // In which domain is located r ?
203 // ----------------------------
204 lz = - 1 ;
205 double rmax = 0. ;
206 double ftp = double(0) ;
207 double gtp = double(0) ;
208 for (int l=0; l<nz; l++) {
209
210 switch( mg->get_type_r(l) ) {
211
212 case RARE: {
213 ftp = ff.val_point(l, 0, theta, pphi) ; // value of F_l(theta,phi)
214 gtp = gg.val_point(l, 0, theta, pphi) ; // value of G_l(theta,phi)
215 rmax = alpha[l] * ( double(1) + ftp + gtp ) + beta[l] ;
216 break ;
217 }
218
219 case FIN: {
220 ftp = double(0) ;
221 gtp = gg.val_point(l, 0, theta, pphi) ;
222 rmax = alpha[l] * ( double(1) + gtp ) + beta[l] ;
223 break ;
224 }
225
226 case UNSURR: {
227 ftp = double(0) ;
228 gtp = double(0) ;
229 rmax = double(1) / ( alpha[l] + beta[l] ) ;
230 break ;
231 }
232
233 default: {
234 cout << "Map_et::val_lx: unknown type_r ! " << endl ;
235 abort() ;
236 }
237 } // end of switch on type_r
238
239
240 if ( rr <= rmax + 1.e-14 ) {
241 lz = l ;
242 if (ftp == double(0)) ftp = ff.val_point(l, 0, theta, pphi) ;
243 if (gtp == double(0)) gtp = gg.val_point(l, 0, theta, pphi) ;
244 break ;
245 }
246 } // End of loop onto the domains
247
248 if (lz == -1) { // The domain has not been found
249 cout.precision(16);
250 cout.setf(ios::showpoint);
251 cout << "Map_et::val_lx: the domain containing r = " << rr <<
252 " has not been found ! "
253 << endl ;
254 for (int l=0; l<nz; l++) {
255 ftp = ff.val_point(l, 0, theta, pphi) ;
256 gtp = gg.val_point(l, 0, theta, pphi) ;
257 switch( mg->get_type_r(l) ) {
258 case RARE: {
259 rmax = alpha[l] * ( double(1) + ftp + gtp ) + beta[l] ;
260 break ;
261 }
262 case FIN: {
263 rmax = alpha[l] * ( double(1) + gtp ) + beta[l] ;
264 break ;
265 }
266 case UNSURR: {
267 rmax = double(1) / ( alpha[l] + beta[l] ) ;
268 break ;
269 }
270 default: {
271 cout << "Map_et::val_lx: unknown type_r ! " << endl ;
272 abort () ;
273 }
274 } // end of switch on type_r
275
276 cout << "domain: " << l << " theta = " << theta << " phi = "
277 << pphi << " : rmax = " << rmax << endl ;
278 }
279 abort() ;
280 }
281
282 // Computation of xi
283 // -----------------
284
285 if ( (rr >= rmax) && ( rr <= rmax + 1.e-14) ) {
286 xi = double(1) ;
287 }
288 else {
289
290 // Search of xi by the secant method
291 // ---------------------------------
292
293 Param parzerosec ;
294 parzerosec.add_double(rr, 0) ;
295 parzerosec.add_double(ftp, 1) ;
296 parzerosec.add_double(gtp, 2) ;
297 parzerosec.add_double(alpha[lz], 3) ;
298 parzerosec.add_double(beta[lz], 4) ;
299
300
301 switch( mg->get_type_r(lz) ) {
302
303 case RARE: {
304 if ( (ff.get_etat()==ETATZERO) && (gg.get_etat()==ETATZERO) ) {
305 xi = ( rr - beta[lz] ) / alpha[lz] ;
306 }
307 else {
308 double xmin = 0 ;
309 double xmax = 1 ;
310 xi = zerosec(fonc_invr_map_et_noyau, parzerosec, xmin, xmax,
311 precis, nitermax, niter) ;
312 }
313 break ;
314 }
315
316 case FIN: {
317 if ( (ff.get_etat()==ETATZERO) && (gg.get_etat()==ETATZERO) ) {
318 xi = ( rr - beta[lz] ) / alpha[lz] ;
319 }
320 else {
321 double xmin = -1 ;
322 double xmax = 1 ;
323 xi = zerosec(fonc_invr_map_et_coq, parzerosec, xmin, xmax,
324 precis, nitermax, niter) ;
325 }
326 break ;
327 }
328
329 case UNSURR: {
330 if ( (ff.get_etat()==ETATZERO) ) {
331 xi = ( double(1)/rr - beta[lz] ) / alpha[lz] ;
332 }
333 else {
334 assert(ff.get_etat()==ETATQCQ) ;
335 if ( ff.c->t[lz]->get_etat() == ETATZERO) {
336 xi = ( double(1) / rr - beta[lz] ) / alpha[lz] ;
337 }
338 else {
339 double xmin = -1 ;
340 double xmax = 1 ;
341 xi = zerosec(fonc_invr_map_et_zec, parzerosec, xmin, xmax,
342 precis, nitermax, niter) ;
343 }
344 }
345 break ;
346 }
347
348 default: {
349 cout << "Map_et::val_lx: unknown type_r ! " << endl ;
350 abort () ;
351 }
352 }
353
354
355 } // End of search by the secant method
356
357}
358
359
360
361
362 //------------------------------//
363 // val_r_jk //
364 //------------------------------//
365
366
367double Map_et::val_r_jk(int l, double xi, int j, int k) const {
368
369 assert( l>=0 ) ;
370 assert( l<mg->get_nzone() ) ;
371
372 double resu ;
373 double ftp = ff(l, k, j, 0) ; // value of F_l(theta_j, phi_k)
374
375 switch( mg->get_type_r(l) ) {
376
377 case RARE: {
378 double gtp = gg(l, k, j, 0) ; // value of G_l(theta_j, phi_k)
379 double xi_2 = xi * xi ;
380 double xi_3 = xi * xi_2 ;
381 double a = xi_2 * xi_2 * (3. - 2.*xi_2) ;
382 double b = ( 2.5 - 1.5 * xi_2 ) * xi_3 ;
383 resu = alpha[l] * ( xi + a * ftp + b * gtp ) + beta[l] ;
384 break ;
385 }
386
387 case FIN: {
388 double gtp = gg(l, k, j, 0) ; // value of G_l(theta_j, phi_k)
389 double xm1 = xi - 1. ;
390 double xp1 = xi + 1. ;
391 double a = 0.25* xm1 * xm1 * (xi + 2.) ;
392 double b = 0.25* xp1 * xp1 * (2. - xi) ;
393 resu = alpha[l] * ( xi + a * ftp + b * gtp ) + beta[l] ;
394 break ;
395 }
396
397 case UNSURR: {
398 double xm1 = xi - 1. ;
399 double a = 0.25* xm1 * xm1 * (xi + 2.) ;
400 resu = double(1) / ( alpha[l] * ( xi + a * ftp ) + beta[l] ) ;
401 break ;
402 }
403
404 default: {
405 cout << "Map_et::val_r_jk: unknown type_r ! " << endl ;
406 abort () ;
407 }
408 }
409
410 return resu ;
411
412}
413
414 //------------------------------//
415 // val_lx_jk //
416 //------------------------------//
417
418void Map_et::val_lx_jk(double rr, int j, int k, const Param& par,
419 int& lz, double& xi) const {
420
421 int nz = mg->get_nzone() ;
422
423 // Precision in the secant method :
424 // ------------------------------
425
426 int nitermax = par.get_int() ;
427 int& niter = par.get_int_mod() ;
428 niter = 0 ; // initialisation
429 double precis = par.get_double() ;
430
431 // Particular case of r = infinity
432 //--------------------------------
433 if (rr == __infinity) {
434 if ( (mg->get_type_r(nz-1) == UNSURR) &&
435 (alpha[nz-1] == - beta[nz-1]) ) {
436 lz = nz - 1 ;
437 xi = 1 ;
438 }
439 else {
440 cout.precision(16);
441 cout.setf(ios::showpoint);
442 cout << "Map_et::val_lx_jk: the domain containing r = " << rr <<
443 " has not been found ! "
444 << endl ;
445 abort () ;
446 }
447 return ;
448 }
449
450
451 // In which domain is located r ?
452 // ----------------------------
453 lz = - 1 ;
454 double rmax = 0 ;
455 double ftp = double(0) ;
456 double gtp = double(0) ;
457 for (int l=0; l<nz; l++) {
458
459 switch( mg->get_type_r(l) ) {
460
461 case RARE: {
462 ftp = ff(l, k, j, 0) ; // value of F_l(theta_j, phi_k)
463 gtp = gg(l, k, j, 0) ; // value of G_l(theta_j, phi_k)
464 rmax = alpha[l] * ( double(1) + ftp + gtp ) + beta[l] ;
465 break ;
466 }
467
468 case FIN: {
469 ftp = double(0) ;
470 gtp = gg(l, k, j, 0) ; // value of G_l(theta_j, phi_k)
471 rmax = alpha[l] * ( double(1) + gtp ) + beta[l] ;
472 break ;
473 }
474
475 case UNSURR: {
476 ftp = double(0) ;
477 gtp = double(0) ;
478 rmax = double(1) / ( alpha[l] + beta[l] ) ;
479 break ;
480 }
481
482 default: {
483 cout << "Map_et::val_lx_jk: unknown type_r ! " << endl ;
484 abort() ;
485 }
486 } // end of switch on type_r
487
488
489 if ( rr <= rmax + 1.e-14 ) {
490 lz = l ;
491 if (ftp == double(0)) ftp = ff(l, k, j, 0) ;
492 if (gtp == double(0)) gtp = gg(l, k, j, 0) ;
493 break ;
494 }
495 } // End of loop onto the domains
496
497 if (lz == -1) { // The domain has not been found
498 cout.precision(16);
499 cout.setf(ios::showpoint);
500 cout << "Map_et::val_lx_jk: the domain containing r = " << rr <<
501 " has not been found ! "
502 << endl ;
503 for (int l=0; l<nz; l++) {
504 ftp = ff(l, k, j, 0) ;
505 gtp = gg(l, k, j, 0) ;
506 switch( mg->get_type_r(l) ) {
507 case RARE: {
508 rmax = alpha[l] * ( double(1) + ftp + gtp ) + beta[l] ;
509 break ;
510 }
511 case FIN: {
512 rmax = alpha[l] * ( double(1) + gtp ) + beta[l] ;
513 break ;
514 }
515 case UNSURR: {
516 rmax = double(1) / ( alpha[l] + beta[l] ) ;
517 break ;
518 }
519 default: {
520 cout << "Map_et::val_lx_jk: unknown type_r ! " << endl ;
521 abort () ;
522 }
523 } // end of switch on type_r
524
525 cout << "domain: " << l << " j = " << j << " k = " << k
526 << " : rmax = " << rmax << endl ;
527 }
528 abort() ;
529 }
530
531 // Computation of xi
532 // -----------------
533
534 if ( (rr >= rmax) && ( rr <= rmax + 1.e-14) ) {
535 xi = double(1) ;
536 }
537 else {
538
539 // Search of xi by the secant method
540 // ---------------------------------
541
542 Param parzerosec ;
543 parzerosec.add_double(rr, 0) ;
544 parzerosec.add_double(ftp, 1) ;
545 parzerosec.add_double(gtp, 2) ;
546 parzerosec.add_double(alpha[lz], 3) ;
547 parzerosec.add_double(beta[lz], 4) ;
548
549
550 switch( mg->get_type_r(lz) ) {
551
552 case RARE: {
553 if ( (ff.get_etat()==ETATZERO) && (gg.get_etat()==ETATZERO) ) {
554 xi = ( rr - beta[lz] ) / alpha[lz] ;
555 }
556 else {
557 double xmin = 0 ;
558 double xmax = 1 ;
559 xi = zerosec(fonc_invr_map_et_noyau, parzerosec, xmin, xmax,
560 precis, nitermax, niter) ;
561 }
562 break ;
563 }
564
565 case FIN: {
566 if ( (ff.get_etat()==ETATZERO) && (gg.get_etat()==ETATZERO) ) {
567 xi = ( rr - beta[lz] ) / alpha[lz] ;
568 }
569 else {
570 double xmin = -1 ;
571 double xmax = 1 ;
572 xi = zerosec(fonc_invr_map_et_coq, parzerosec, xmin, xmax,
573 precis, nitermax, niter) ;
574 }
575 break ;
576 }
577
578 case UNSURR: {
579 if ( (ff.get_etat()==ETATZERO) ) {
580 xi = ( double(1)/rr - beta[lz] ) / alpha[lz] ;
581 }
582 else {
583 assert(ff.get_etat()==ETATQCQ) ;
584 if ( ff.c->t[lz]->get_etat() == ETATZERO) {
585 xi = ( double(1) / rr - beta[lz] ) / alpha[lz] ;
586 }
587 else {
588 double xmin = -1 ;
589 double xmax = 1 ;
590 xi = zerosec(fonc_invr_map_et_zec, parzerosec, xmin, xmax,
591 precis, nitermax, niter) ;
592 }
593 }
594 break ;
595 }
596
597 default: {
598 cout << "Map_et::val_lx_jk: unknown type_r ! " << endl ;
599 abort () ;
600 }
601 }
602
603
604 } // End of search by the secant method
605}
606
607
608//=============================================================================//
609// Auxilary functions //
610//=============================================================================//
611
612double fonc_invr_map_et_noyau(double x, const Param& par) {
613
614 double r = par.get_double(0) ;
615 double f = par.get_double(1) ;
616 double g = par.get_double(2) ;
617 double alp = par.get_double(3) ;
618 double bet = par.get_double(4) ;
619 double x_2 = x * x ;
620 double x_3 = x_2 * x ;
621 double a = x_2 * x_2 * (3. - 2.*x_2) ;
622 double b = ( 2.5 - 1.5 * x_2 ) * x_3 ;
623
624 return alp * ( x + a * f + b * g ) + bet - r ;
625
626}
627
628//****************************************************************************
629
630double fonc_invr_map_et_coq(double x, const Param& par) {
631
632 double r = par.get_double(0) ;
633 double f = par.get_double(1) ;
634 double g = par.get_double(2) ;
635 double alp = par.get_double(3) ;
636 double bet = par.get_double(4) ;
637 double xm1 = x - 1. ;
638 double xp1 = x + 1. ;
639 double a = 0.25* xm1 * xm1 * (x + 2.) ;
640 double b = 0.25* xp1 * xp1 * (2. - x) ;
641
642 return alp * ( x + a * f + b * g ) + bet - r ;
643
644}
645
646//****************************************************************************
647
648double fonc_invr_map_et_zec(double x, const Param& par) {
649
650 double r = par.get_double(0) ;
651 double f = par.get_double(1) ;
652 double alp = par.get_double(3) ;
653 double bet = par.get_double(4) ;
654 double xm1 = x - 1. ;
655 double a = 0.25* xm1 * xm1 * (x + 2.) ;
656
657 return alp * ( x + a * f ) + bet - double(1) / r ;
658
659}
660
661}
double * beta
Array (size: mg->nzone ) of the values of in each domain.
Definition map.h:2760
virtual void val_lx_jk(double rr, int j, int k, const Param &par, int &l, double &xi) const
Computes the domain index l and the value of corresponding to a point of arbitrary r but collocation...
virtual double val_r(int l, double xi, double theta, double pphi) const
Returns the value of the radial coordinate r for a given in a given domain.
Valeur ff
Values of the function at the nt*np angular collocation points in each domain.
Definition map.h:2819
virtual void val_lx(double rr, double theta, double pphi, int &l, double &xi) const
Computes the domain index l and the value of corresponding to a point given by its physical coordina...
double * alpha
Array (size: mg->nzone ) of the values of in each domain.
Definition map.h:2758
Valeur gg
Values of the function at the nt*np angular collocation points in each domain.
Definition map.h:2826
virtual double val_r_jk(int l, double xi, int j, int k) const
Returns the value of the radial coordinate r for a given and a given collocation point in in a give...
const Mg3d * mg
Pointer on the multi-grid Mgd3 on which this is defined
Definition map.h:676
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition grilles.h:474
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
Definition mtbl.h:132
Parameter storage.
Definition param.h:125
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition param.C:315
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition param.C:292
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition param.C:361
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
Definition param.C:385
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
Definition param.C:430
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition param.C:246
int get_etat() const
Gives the logical state.
Definition tbl.h:394
double val_point(int l, double x, double theta, double phi) const
Computes the value of the field represented by *this at an arbitrary point, by means of the spectral ...
Definition valeur.C:882
int get_etat() const
Returns the logical state.
Definition valeur.h:726
Mtbl * c
Values of the function at the points of the multi-grid
Definition valeur.h:299
double zerosec(double(*f)(double, const Param &), const Param &par, double a, double b, double precis, int nitermax, int &niter, bool abort=true)
Finding the zero a function.
Definition zerosec.C:89
Lorene prototypes.
Definition app_hor.h:64