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 
31 char 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 
65 namespace Lorene {
66 double et_rot_diff_fzero(double omeg, const Param& par) ;
67 
68 
69  //-----------------------------------//
70  // fait_omega_field //
71  //-----------------------------------//
72 
73 void 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 
172 double 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 }
const Etoile & get_etoile(int position=0) const
Returns the reference of a Etoile stored in the list.
Definition: param.C:1669
const Cmp & get_cmp(int position=0) const
Returns the reference of a Cmp stored in the list.
Definition: param.C:980
Tenseur prim_field
Field .
Definition: et_rot_diff.h:111
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:443
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_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:348
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1170
Lorene prototypes.
Definition: app_hor.h:64
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:281
Tenseur nnn
Total lapse function.
Definition: etoile.h:509
Tenseur nphi
Metric coefficient .
Definition: etoile.h:1510
double(* primfrot)(double, const Tbl &)
Primitive of the function , which vanishes at the stellar center.
Definition: et_rot_diff.h:93
void fait_omega_field(double omeg_min, double omeg_max, double precis, int nitermax)
Computes (member omega_field ).
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition: cmp_math.C:458
double funct_omega(double omeg) const
Evaluates , where F is the function defining the rotation profile.
Definition: et_rot_diff.C:315
double omega_max
Maximum value of .
Definition: et_rot_diff.h:108
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
Class for differentially rotating stars.
Definition: et_rot_diff.h:70
void fait_prim_field()
Computes the member prim_field from omga_field .
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition: param.C:292
void add_cmp(const Cmp &ti, int position=0)
Adds the address of a new Cmp to the list.
Definition: param.C:935
Tbl par_frot
Parameters of the function .
Definition: et_rot_diff.h:102
Parameter storage.
Definition: param.h:125
void mult_rsint()
Multiplication by .
Definition: cmp_r_manip.C:116
Map & mp
Mapping associated with the star.
Definition: etoile.h:429
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:435
int get_etat() const
Returns the logical state.
Definition: tenseur.h:704
Tenseur bbb
Metric factor B.
Definition: etoile.h:1504
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition: cmp.C:644
void add_etoile(const Etoile &eti, int position=0)
Adds the address of a new Etoile to the list.
Definition: param.C:1624
int nzet
Number of domains of *mp occupied by the star.
Definition: etoile.h:432
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 grid.
Definition: grilles.h:273
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:721
void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition: tenseur.C:657
Basic array class.
Definition: tbl.h:161
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
double omega_min
Minimum value of .
Definition: et_rot_diff.h:107
Tenseur omega_field
Field .
Definition: et_rot_diff.h:105