LORENE
valeur_equipot.C
1/*
2 * Method of the class Valeur to compute an equipotential surface.
3 *
4 * (see file valeur.h for the documentation).
5 */
6
7/*
8 * Copyright (c) 1999-2001 Eric Gourgoulhon
9 *
10 * This file is part of LORENE.
11 *
12 * LORENE is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU General Public License as published by
14 * the Free Software Foundation; either version 2 of the License, or
15 * (at your option) any later version.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28
29char valeur_equipot_C[] = "$Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_equipot.C,v 1.5 2014/10/13 08:53:50 j_novak Exp $" ;
30
31/*
32 * $Id: valeur_equipot.C,v 1.5 2014/10/13 08:53:50 j_novak Exp $
33 * $Log: valeur_equipot.C,v $
34 * Revision 1.5 2014/10/13 08:53:50 j_novak
35 * Lorene classes and functions now belong to the namespace Lorene.
36 *
37 * Revision 1.4 2014/10/06 15:13:23 j_novak
38 * Modified #include directives to use c++ syntax.
39 *
40 * Revision 1.3 2003/12/19 15:05:56 j_novak
41 * Added some initializations
42 *
43 * Revision 1.2 2002/10/16 14:37:15 j_novak
44 * Reorganization of #include instructions of standard C++, in order to
45 * use experimental version 3 of gcc.
46 *
47 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
48 * LORENE
49 *
50 * Revision 1.3 2000/01/28 17:06:37 eric
51 * Changement du cas l2<nz_search-1.
52 *
53 * Revision 1.2 2000/01/03 15:07:50 eric
54 * *** empty log message ***
55 *
56 * Revision 1.1 1999/12/29 13:12:08 eric
57 * Initial revision
58 *
59 *
60 * $Header: /cvsroot/Lorene/C++/Source/Valeur/valeur_equipot.C,v 1.5 2014/10/13 08:53:50 j_novak Exp $
61 *
62 */
63
64// Headers C
65#include <cstdlib>
66#include <cmath>
67
68// Headers Lorene
69#include "valeur.h"
70#include "itbl.h"
71#include "param.h"
72#include "nbr_spx.h"
73#include "utilitaires.h"
74
75// Local prototypes
76namespace Lorene {
77double valeur_equipot_fonc(double, const Param&) ;
78
79//****************************************************************************
80
81void Valeur::equipot(double uu0, int nz_search, double precis, int nitermax,
82 int& niter, Itbl& l_iso, Tbl& xi_iso) const {
83
84
85 int nz = mg->get_nzone() ;
86 int nt = mg->get_nt(0) ;
87 int np = mg->get_np(0) ;
88
89 // Protections
90 // -----------
91 assert(etat == ETATQCQ) ;
92 assert(nz_search > 0) ;
93 assert(nz_search <= nz) ;
94 for (int l=1; l<nz_search; l++) {
95 assert( mg->get_nt(l) == nt ) ;
96 assert( mg->get_np(l) == np ) ;
97 }
98
99 coef() ; // the spectral coefficients are required
100
101 l_iso.set_etat_qcq() ;
102 xi_iso.set_etat_qcq() ;
103
104 Param parf ;
105 int j, k, l2 ;
106 parf.add_int_mod(j, 0) ;
107 parf.add_int_mod(k, 1) ;
108 parf.add_int_mod(l2, 2) ;
109 parf.add_double(uu0) ;
110 parf.add_mtbl_cf(*c_cf) ;
111
112
113 // Loop of phi and theta
114 // ---------------------
115
116 niter = 0 ;
117
118 for (k=0; k<np; k++) {
119
120 for (j=0; j<nt; j++) {
121
122
123// 1. Recherche du point de collocation en r (l2,i2) egal a ou situe
124// immediatement avant r_iso(phi,theta)
125
126 int i2 = 0;
127 l2 = nz ;
128 for (int l=nz_search-1; l>= 0; l--) { // On part de la zone la plus extreme
129 int nr = mg->get_nr(l) ;
130
131 for (int i=nr-1; i >= 0; i--) {
132 double uux = (*this)(l, k, j, i) ;
133 if ( ( (uux > uu0) || ( fabs(uux-uu0) < 1.e-14 ) ) &&
134 (uux != __infinity) ) {
135 l2 = l ; // on a trouve le point
136 i2 = i ; //
137 break ;
138 }
139 }
140
141 if (l2 == l) break ;
142 } // fin de la boucle sur les zones
143
144 if (l2==nz) {
145 cout << "Valeur::equipot: the point uu >= uu0 has not been found" << endl ;
146 cout << " for the phi index " << k << endl ;
147 cout << " and the theta index " << j << endl ;
148 cout << " uu0 = " << uu0 << endl ;
149 abort () ;
150 }
151
152// 2. Point suivant (l2,i3)
153 int i3 = 0 ;
154 if (i2 < mg->get_nr(l2) -1) { // on reste dans la meme zone
155 i3 = i2 + 1 ;
156 }
157 else {
158 if (l2<nz_search-1) { // on change de zone
159
160 double uux = (*this)(l2, k, j, i2) ;
161 if ( ( fabs(uux-uu0) < 1.e-14 ) ) {
162 // it is OK
163 cout <<
164 "Valeur::equipot: WARNING : fabs(uux-uu0) < 1.e-14" << endl ;
165 l_iso.set(k, j) = l2 ;
166 xi_iso.set(k, j) = (mg->get_grille3d(l2))->x[i2] ;
167 }
168 else{
169 cout << "Valeur::equipot: PROBLEM !!!" << endl ;
170 cout << " k, j : " << k << " " << j << endl ;
171 cout << " uu0 : " << uu0 << endl ;
172 abort () ;
173 }
174 }
175 else { // on est a l'extremite de la grille
176 l_iso.set(k, j) = nz_search-1 ;
177 xi_iso.set(k, j) = (mg->get_grille3d(l2))->x[i2] ;
178 continue ; // on passe au theta suivant
179 }
180 }
181
182 l_iso.set(k, j) = l2 ;
183
184 double x2 = (mg->get_grille3d(l2))->x[i2] ;
185 double x3 = (mg->get_grille3d(l2))->x[i3] ;
186
187 double y2 = (*this)(l2, k, j, i2) - uu0 ;
188
189// 3. Recherche de x0
190
191 int niter0 = 0 ;
192 if (fabs(y2) < 1.e-14) {
193 xi_iso.set(k, j) = x2 ;
194 }
195 else {
196 xi_iso.set(k, j) = zerosec(valeur_equipot_fonc, parf, x2, x3, precis,
197 nitermax, niter0 ) ;
198 }
199
200 niter = ( niter0 > niter ) ? niter0 : niter ;
201
202 } // fin de la boucle sur theta
203 } // fin de la boucle sur phi
204
205}
206
207
208
209//****************************************************************************
210
211double valeur_equipot_fonc(double xi, const Param& par) {
212
213 int j = par.get_int_mod(0) ;
214 int k = par.get_int_mod(1) ;
215 int l = par.get_int_mod(2) ;
216 double uu0 = par.get_double() ;
217 const Mtbl_cf& cuu = par.get_mtbl_cf() ;
218
219 return cuu.val_point_jk(l, xi, j, k) - uu0 ;
220}
221
222}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Basic integer array class.
Definition itbl.h:122
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
Parameter storage.
Definition param.h:125
Basic array class.
Definition tbl.h:161
void equipot(double uu0, int nz_search, double precis, int nitermax, int &niter, Itbl &l_iso, Tbl &xi_iso) const
Determines an equipotential surface of the field represented by *this (inward search).
const Mg3d * mg
Multi-grid Mgd3 on which this is defined.
Definition valeur.h:292
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
void coef() const
Computes the coeffcients of *this.
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
Definition valeur.h:295
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