LORENE
et_rot_diff_faitomeg.C
1/*
2 * Functions Et_rot_diff::fait_omega_field
3 * Et_rot_diff::fait_prim_field
4 *
5 * (see file et_rot_diff.h for documentation).
6 *
7 */
8
9/*
10 * Copyright (c) 2001 Eric Gourgoulhon
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
30
31char et_rot_diff_faitomeg_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_diff_faitomeg.C,v 1.3 2014/10/13 08:52:57 j_novak Exp $" ;
32
33/*
34 * $Id: et_rot_diff_faitomeg.C,v 1.3 2014/10/13 08:52:57 j_novak Exp $
35 * $Log: et_rot_diff_faitomeg.C,v $
36 * Revision 1.3 2014/10/13 08:52:57 j_novak
37 * Lorene classes and functions now belong to the namespace Lorene.
38 *
39 * Revision 1.2 2003/05/14 20:07:00 e_gourgoulhon
40 * Suppressed the outputs (cout)
41 *
42 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
43 * LORENE
44 *
45 * Revision 1.3 2001/10/25 09:43:18 eric
46 * omega_min est determine dans l'etoile seulement (l<nzet).
47 *
48 * Revision 1.2 2001/10/25 09:21:29 eric
49 * Recherche de Omega dans un intervalle de 20% autour de la valeur precedente.
50 *
51 * Revision 1.1 2001/10/19 08:18:23 eric
52 * Initial revision
53 *
54 *
55 * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_diff_faitomeg.C,v 1.3 2014/10/13 08:52:57 j_novak Exp $
56 *
57 */
58
59
60// Headers Lorene
61#include "et_rot_diff.h"
62#include "utilitaires.h"
63#include "param.h"
64
65namespace Lorene {
66double et_rot_diff_fzero(double omeg, const Param& par) ;
67
68
69 //-----------------------------------//
70 // fait_omega_field //
71 //-----------------------------------//
72
73void Et_rot_diff::fait_omega_field(double omeg_min, double omeg_max,
74 double precis, int nitermax) {
75
76 const Mg3d& mg = *(mp.get_mg()) ;
77 int nz = mg.get_nzone() ;
78 int nzm1 = nz - 1 ;
79
80 // Computation of B^2 r^2 sin^2(theta):
81 // -----------------------------------
82
83 Cmp brst2 = bbb() ;
84 brst2.annule(nzm1) ;
85 brst2.std_base_scal() ;
86 brst2.mult_rsint() ; // Multiplication by r sin(theta)
87 brst2 = brst2*brst2 ;
88
89 Cmp nnn2 = nnn() * nnn() ;
90
91 Param par ;
92 par.add_cmp(nnn2, 0) ;
93 par.add_cmp(brst2, 1) ;
94 par.add_cmp(nphi(), 2) ;
95
96 int l, k, j, i ;
97 par.add_int(l, 0) ;
98 par.add_int(k, 1) ;
99 par.add_int(j, 2) ;
100 par.add_int(i, 3) ;
101
102 par.add_etoile(*this) ;
103
104 // Loop on the grid points
105 // -----------------------
106
107 int niter ;
108
109 bool prev_zero = (omega_field.get_etat() == ETATZERO) ;
110
112
113 // cout << "Et_rot_diff::fait_omega_field: niter : " << endl ;
114 for (l=0; l<nzet+1; l++) {
115 Tbl& tom = omega_field.set().set(l) ;
116 for (k=0; k<mg.get_np(l); k++) {
117 for (j=0; j<mg.get_nt(l); j++) {
118 for (i=0; i<mg.get_nr(l); i++) {
119
120 double& omeg = tom.set(k, j, i) ;
121
122 double omeg_min1, omeg_max1 ;
123 if ( prev_zero || omeg == double(0)) {
124 omeg_min1 = omeg_min ;
125 omeg_max1 = omeg_max ;
126 }
127 else{
128 omeg_min1 = 0.8 * omeg ;
129 omeg_max1 = 1.2 * omeg ;
130 }
131
132 omeg = zerosec(et_rot_diff_fzero, par, omeg_min1,
133 omeg_max1, precis, nitermax, niter) ;
134
135 // cout << " " << niter ;
136
137 }
138 }
139 }
140 // cout << endl ;
141 }
142
143 // omega_field is set to 0 far from the star:
144 // ---------------------------------------------------
145 for (l=nzet+1; l<nz; l++) {
146 omega_field.set().set(l) = 0 ;
147 }
148
150
151 // Min and Max of Omega
152 // --------------------
153
154 omega_min = min( omega_field()(0) ) ;
155 for (l=1; l<nzet; l++) {
156 double xx = min( omega_field()(l) ) ;
157 omega_min = (xx < omega_min) ? xx : omega_min ;
158 }
159
160 omega_max = max( max( omega_field() ) ) ;
161
162 // Update of prim_field
163 // --------------------
164 fait_prim_field() ;
165
166}
167
168 //-----------------------------------//
169 // et_rot_diff_fzero //
170 //-----------------------------------//
171
172double et_rot_diff_fzero(double omeg, const Param& par) {
173
174 const Cmp& nnn2 = par.get_cmp(0) ;
175 const Cmp& brst2 = par.get_cmp(1) ;
176 const Cmp& nphi = par.get_cmp(2) ;
177 int l = par.get_int(0) ;
178 int k = par.get_int(1) ;
179 int j = par.get_int(2) ;
180 int i = par.get_int(3) ;
181
182 const Et_rot_diff* et =
183 dynamic_cast<const Et_rot_diff*>(&par.get_etoile()) ;
184
185 double fom = et->funct_omega(omeg) ;
186 double omnp = omeg - nphi(l, k, j, i) ;
187
188 return fom - brst2(l, k, j, i) * omnp /
189 ( nnn2(l, k, j, i) - brst2(l, k, j, i) * omnp*omnp ) ;
190
191}
192
193
194 //-----------------------------------//
195 // fait_prim_field //
196 //-----------------------------------//
197
198
200
201 const Mg3d& mg = *(mp.get_mg()) ;
202 int nz = mg.get_nzone() ;
203
204 // Loop on the grid points in the vicinity of the star
205 // ---------------------------------------------------
206
208
209 for (int l=0; l<nzet+1; l++) {
210 Tbl& tprim = prim_field.set().set(l) ;
211 for (int k=0; k<mg.get_np(l); k++) {
212 for (int j=0; j<mg.get_nt(l); j++) {
213 for (int i=0; i<mg.get_nr(l); i++) {
214
215 tprim.set(k, j, i) =
216 primfrot( omega_field()(l, k, j, i), par_frot ) ;
217
218 }
219 }
220 }
221 }
222
223 // prim_field is set to 0 far from the star:
224 // -----------------------------------------
225 for (int l=nzet+1; l<nz; l++) {
226 prim_field.set().set(l) = 0 ;
227 }
228
230
231
232}
233}
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
void mult_rsint()
Multiplication by .
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
Tbl & set(int l)
Read/write of the value in a given domain.
Definition cmp.h:724
Class for differentially rotating stars.
Definition et_rot_diff.h:70
double omega_min
Minimum value of .
void fait_omega_field(double omeg_min, double omeg_max, double precis, int nitermax)
Computes (member omega_field ).
Tenseur prim_field
Field .
double omega_max
Maximum value of .
Tenseur omega_field
Field .
void fait_prim_field()
Computes the member prim_field from omga_field .
Tbl par_frot
Parameters of the function .
double(* primfrot)(double, const Tbl &)
Primitive of the function , which vanishes at the stellar center.
Definition et_rot_diff.h:93
double funct_omega(double omeg) const
Evaluates , where F is the function defining the rotation profile.
Tenseur nphi
Metric coefficient .
Definition etoile.h:1510
Tenseur bbb
Metric factor B.
Definition etoile.h:1504
int nzet
Number of domains of *mp occupied by the star.
Definition etoile.h:432
Tenseur nnn
Total lapse function.
Definition etoile.h:509
Map & mp
Mapping associated with the star.
Definition etoile.h:429
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Multi-domain grid.
Definition grilles.h:273
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
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
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
void add_etoile(const Etoile &eti, int position=0)
Adds the address of a new Etoile to the list.
Definition param.C:1624
void add_cmp(const Cmp &ti, int position=0)
Adds the address of a new Cmp to the list.
Definition param.C:935
const Cmp & get_cmp(int position=0) const
Returns the reference of a Cmp stored in the list.
Definition param.C:980
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition param.C:246
const Etoile & get_etoile(int position=0) const
Returns the reference of a Etoile stored in the list.
Definition param.C:1669
Basic array class.
Definition tbl.h:161
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition tenseur.C:657
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition tenseur.C:1170
int get_etat() const
Returns the logical state.
Definition tenseur.h:707
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition cmp_math.C:458
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
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