LORENE
strot_diff_faitomeg.C
1 /*
2  * Functions Star_rot_diff::fait_omega_field
3  * Star_rot_diff::fait_prim_field
4  *
5  * (see file star_rot_diff.h for documentation).
6  *
7  */
8 
9 /*
10  * Copyright (c) 2025 Santiago Jaraba
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 // Lorene headers
31 #include "star_rot_diff.h"
32 #include "utilitaires.h"
33 #include "param.h"
34 
35 namespace Lorene {
36  double strot_diff_fzero(double omeg, const Param& par) ;
37 
38  //-----------------------------------//
39  // fait_omega_field //
40  //-----------------------------------//
41 
42  void Star_rot_diff::fait_omega_field(double omeg_min, double omeg_max,
43  double precis, int nitermax) {
44 
45  const Mg3d& mg = *(mp.get_mg()) ;
46  int nz = mg.get_nzone() ;
47  int nzm1 = nz - 1 ;
48 
49  // Computation of B^2 r^2 sin^2(theta):
50  // -----------------------------------
51 
52  Scalar brst2 = bbb ;
53  brst2.annule_domain(nzm1) ;
54  brst2.std_spectral_base() ;
55  brst2.mult_rsint() ; // Multiplication by r sin(theta)
56  brst2 = brst2*brst2 ;
57 
58  Scalar nnn2 = nn * nn ;
59 
60  Param par ;
61  par.add_scalar(nnn2, 0) ;
62  par.add_scalar(brst2, 1) ;
63  par.add_scalar(nphi, 2) ;
64 
65  int l, k, j, i ;
66  par.add_int(l, 0) ;
67  par.add_int(k, 1) ;
68  par.add_int(j, 2) ;
69  par.add_int(i, 3) ;
70 
71  par.add_star(*this) ;
72 
73  // Loop on the grid points
74  // -----------------------
75 
76  int niter ;
77 
78  bool prev_zero = (omega_field.get_etat() == ETATZERO) ;
79 
81 
82  for (l=0; l<nzet+1; l++) {
83  Tbl& tom = omega_field.set_domain(l) ;
84  for (k=0; k<mg.get_np(l); k++) {
85  for (j=0; j<mg.get_nt(l); j++) {
86  for (i=0; i<mg.get_nr(l); i++) {
87 
88  double& omeg = tom.set(k, j, i) ;
89 
90  double omeg_min1, omeg_max1 ;
91  if ( prev_zero || omeg == double(0)) {
92  omeg_min1 = omeg_min ;
93  omeg_max1 = omeg_max ;
94  }
95  else{
96  omeg_min1 = 0.8 * omeg ;
97  omeg_max1 = 1.2 * omeg ;
98  }
99 
100  omeg = zerosec(strot_diff_fzero, par, omeg_min1,
101  omeg_max1, precis, nitermax, niter) ;
102 
103  }
104  }
105  }
106  }
107 
108  // omega_field is set to 0 far from the star:
109  // ---------------------------------------------------
110  for (l=nzet+1; l<nz; l++) {
111  omega_field.set_domain(l) = 0 ;
112  }
113 
115 
116  // Min and Max of Omega
117  // --------------------
118 
119  omega_min = min( omega_field.domain(0) ) ;
120  for (l=1; l<nzet; l++) {
121  double xx = min( omega_field.domain(l) ) ;
122  omega_min = (xx < omega_min) ? xx : omega_min ;
123  }
124 
125  omega_max = max( max( omega_field ) ) ;
126 
127  // Update of prim_field
128  // --------------------
129  fait_prim_field() ;
130 
131 }
132 
133 
134  //----------------------------------------//
135  // strot_diff_fzero //
136  //----------------------------------------//
137 
138 double strot_diff_fzero(double omeg, const Param& par) {
139 
140  const Scalar& nnn2 = par.get_scalar(0) ;
141  const Scalar& brst2 = par.get_scalar(1) ;
142  const Scalar& nphi = par.get_scalar(2) ;
143  int l = par.get_int(0) ;
144  int k = par.get_int(1) ;
145  int j = par.get_int(2) ;
146  int i = par.get_int(3) ;
147 
148  const Star_rot_diff* star =
149  dynamic_cast<const Star_rot_diff*>(&par.get_star()) ;
150 
151  double omnp = omeg - nphi.val_grid_point(l, k, j, i) ;
152  double fom = brst2.val_grid_point(l, k, j, i) * omnp /
153  ( nnn2.val_grid_point(l, k, j, i) -
154  brst2.val_grid_point(l, k, j, i) * omnp * omnp ) ;
155 
156  return omeg - star->omega_funct(fom) ;
157 
158 }
159 
160 
161  //-----------------------------------//
162  // fait_prim_field //
163  //-----------------------------------//
164 
165 
167 
168  const Mg3d& mg = *(mp.get_mg()) ;
169  int nz = mg.get_nzone() ;
170  int nzm1 = nz - 1 ;
171 
172  Scalar brst2 = bbb ;
173  brst2.annule_domain(nzm1) ;
174  brst2.std_spectral_base() ;
175  brst2.mult_rsint() ; // Multiplication by r sin(theta)
176  brst2 = brst2*brst2 ;
177 
178  Scalar nnn2 = nn * nn ;
179 
180  // Loop on the grid points in the vicinity of the star
181  // ---------------------------------------------------
182 
184 
185  for (int l=0; l<nzet+1; l++) {
186  Tbl& tprim = prim_field.set_domain(l) ;
187  for (int k=0; k<mg.get_np(l); k++) {
188  for (int j=0; j<mg.get_nt(l); j++) {
189  for (int i=0; i<mg.get_nr(l); i++) {
190  double omnp = omega_field.val_grid_point(l, k, j, i) - nphi.val_grid_point(l, k, j, i) ;
191  double F = brst2.val_grid_point(l, k, j, i) * omnp /
192  ( nnn2.val_grid_point(l, k, j, i) -
193  brst2.val_grid_point(l, k, j, i) * omnp * omnp ) ;
194 
195  tprim.set(k, j, i) =
196  primfrot( F, par_frot ) ;
197 
198  }
199  }
200  }
201  }
202 
203  // prim_field is set to 0 far from the star:
204  // -----------------------------------------
205  for (int l=nzet+1; l<nz; l++) {
206  prim_field.set_domain(l) = 0 ;
207  }
208 
210 
211 
212 }
213 }
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:676
Class for differentially rotating stars in quasi-isotropic gauge and maximal slicing.
Definition: star_rot_diff.h:39
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition: param.C:249
Scalar prim_field
Field .
Definition: star_rot_diff.h:79
Map & mp
Mapping associated with the star.
Definition: star.h:180
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:481
Tbl par_frot
Parameters of the function .
Definition: star_rot_diff.h:70
const Tbl & domain(int l) const
Read-only of the value in a given domain.
Definition: scalar.h:637
Scalar bbb
Metric factor B.
Definition: star_rot.h:119
Lorene prototypes.
Definition: app_hor.h:67
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:792
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:399
const Star & get_star(int position=0) const
Returns the reference of a Star stored in the list.
Definition: param.C:1738
void fait_prim_field()
Computes the member prim_field from omega_field .
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:790
double omega_funct(double F) const
Evaluates , where F is linked to the components of the fluid 4-velocity by .
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition: scalar.C:371
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:566
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition: cmp_math.C:461
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Definition: scalar.h:627
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:92
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition: scalar.h:649
int nzet
Number of domains of *mp occupied by the star.
Definition: star.h:183
Scalar nphi
Metric coefficient .
Definition: star_rot.h:125
void add_star(const Star &eti, int position=0)
Adds the address of a new Star to the list.
Definition: param.C:1693
const Scalar & get_scalar(int position=0) const
Returns the reference of a Scalar stored in the list.
Definition: param.C:1396
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition: param.C:295
Parameter storage.
Definition: param.h:125
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:467
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:438
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
double omega_max
Maximum value of .
Definition: star_rot_diff.h:76
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:471
Multi-domain grid.
Definition: grilles.h:276
void add_scalar(const Scalar &ti, int position=0)
Adds the address of a new Scalar to the list.
Definition: param.C:1351
Scalar nn
Lapse function N .
Definition: star.h:225
Basic array class.
Definition: tbl.h:164
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:476
double omega_min
Minimum value of .
Definition: star_rot_diff.h:75
double(* primfrot)(double, const Tbl &)
Primitive of the function , which vanishes at the stellar center.
Definition: star_rot_diff.h:61
void fait_omega_field(double omeg_min, double omeg_max, double precis, int nitermax)
Computes (member omega_field ).
Scalar omega_field
Field .
Definition: star_rot_diff.h:73