LORENE
strot_dirac_diff_faitomeg.C
1/*
2 * Functions Star_rot_Dirac_diff::fait_omega_field
3 * Star_rot_Dirac_diff::fait_prim_field
4 *
5 * (see file star_rot_dirac_diff.h for documentation).
6 *
7 */
8
9/*
10 * Copyright (c) 2005 Motoyuki Saijo
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 strot_dirac_diff_faitomeg_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star/strot_dirac_diff_faitomeg.C,v 1.2 2014/10/13 08:53:39 j_novak Exp $" ;
32
33/*
34 * $Header: /cvsroot/Lorene/C++/Source/Star/strot_dirac_diff_faitomeg.C,v 1.2 2014/10/13 08:53:39 j_novak Exp $
35 *
36 */
37
38
39// Headers Lorene
40#include "star_rot_dirac_diff.h"
41#include "utilitaires.h"
42#include "param.h"
43
44namespace Lorene {
45double strot_dirac_diff_fzero(double omeg, const Param& par) ;
46
47
48 //-----------------------------------//
49 // fait_omega_field //
50 //-----------------------------------//
51
53 double precis, int nitermax) {
54
55 const Mg3d& mg = *(mp.get_mg()) ;
56 int nz = mg.get_nzone() ;
57 int nzm1 = nz - 1 ;
58
59 Scalar brst2 = gamma.cov()(3,3) ;
60 brst2.annule(nzm1, nzm1) ;
61 brst2.mult_rsint() ;
62 brst2.mult_rsint() ;
63 Scalar nnn2 = qqq * qqq / psi4 ;
64 Scalar betp = beta(3) ;
65 betp.div_rsint() ;
66
67 Param par ;
68 par.add_scalar(nnn2, 0) ;
69 par.add_scalar(brst2, 1) ;
70 par.add_scalar(betp, 2) ;
71
72 int l, k, j, i ;
73 par.add_int(l, 0) ;
74 par.add_int(k, 1) ;
75 par.add_int(j, 2) ;
76 par.add_int(i, 3) ;
77
78 par.add_star(*this) ;
79
80 // Loop on the grid points
81 // -----------------------
82
83 int niter ;
84
85 bool prev_zero = (omega_field.get_etat() == ETATZERO) ;
86
88
89 // cout << "Star_rot_Dirac_diff::fait_omega_field: niter : " << '\n' ;
90 for (l=0; l<nzet+1; l++) {
92 for (k=0; k<mg.get_np(l); k++) {
93 for (j=0; j<mg.get_nt(l); j++) {
94 for (i=0; i<mg.get_nr(l); i++) {
95
96 double& omeg = tom.set(k, j, i) ;
97
98 double omeg_min1, omeg_max1 ;
99 if ( prev_zero || omeg == double(0)) {
102 }
103 else{
104 omeg_min1 = 0.8 * omeg ;
105 omeg_max1 = 1.2 * omeg ;
106 }
107
108 omeg = zerosec(strot_dirac_diff_fzero, par, omeg_min1,
109 omeg_max1, precis, nitermax, niter) ;
110
111 // cout << " " << niter ;
112
113 }
114 }
115 }
116 // cout << '\n' ;
117 }
118
119 // omega_field is set to 0 far from the star:
120 // ---------------------------------------------------
121 for (l=nzet+1; l<nz; l++) {
123 }
124
126
127 // Min and Max of Omega
128 // --------------------
129
131 for (l=1; l<nzet; l++) {
132 double xx = min( omega_field.domain(l) ) ;
133 omega_min = (xx < omega_min) ? xx : omega_min ;
134 }
135
136 omega_max = max( max( omega_field ) ) ;
137
138 // Update of prim_field
139 // --------------------
140 fait_prim_field() ;
141
142}
143
144 //----------------------------------------//
145 // strot_dirac_diff_fzero //
146 //----------------------------------------//
147
148double strot_dirac_diff_fzero(double omeg, const Param& par) {
149
150 const Scalar& nnn2 = par.get_scalar(0) ;
151 const Scalar& brst2 = par.get_scalar(1) ;
152 const Scalar& betp = par.get_scalar(2) ;
153 int l = par.get_int(0) ;
154 int k = par.get_int(1) ;
155 int j = par.get_int(2) ;
156 int i = par.get_int(3) ;
157
158 const Star_rot_Dirac_diff* star =
159 dynamic_cast<const Star_rot_Dirac_diff*>(&par.get_star()) ;
160
161 double fom = star->funct_omega(omeg) ;
162 double omnp = omeg + betp.val_grid_point(l, k, j, i) ;
163
164 return fom - brst2.val_grid_point(l, k, j, i) * omnp /
165 ( nnn2.val_grid_point(l, k, j, i) -
166 brst2.val_grid_point(l, k, j, i) * omnp * omnp ) ;
167
168}
169
170
171 //-----------------------------------//
172 // fait_prim_field //
173 //-----------------------------------//
174
175
177
178 const Mg3d& mg = *(mp.get_mg()) ;
179 int nz = mg.get_nzone() ;
180
181 // Loop on the grid points in the vicinity of the star
182 // ---------------------------------------------------
183
185
186 for (int l=0; l<nzet+1; l++) {
188 for (int k=0; k<mg.get_np(l); k++) {
189 for (int j=0; j<mg.get_nt(l); j++) {
190 for (int i=0; i<mg.get_nr(l); i++) {
191
192 tprim.set(k, j, i) =
194 par_frot ) ;
195
196 }
197 }
198 }
199 }
200
201 // prim_field is set to 0 far from the star:
202 // -----------------------------------------
203 for (int l=nzet+1; l<nz; l++) {
204 prim_field.set_domain(l) = 0 ;
205 }
206
208
209
210}
211}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition metric.C:280
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
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition scalar.C:784
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:637
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Definition scalar.h:615
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition scalar.C:365
const Tbl & domain(int l) const
Read-only of the value in a given domain.
Definition scalar.h:625
Class for relativistic differentially rotating stars in Dirac gauge and maximal slicing.
void fait_omega_field(double omeg_min, double omeg_max, double precis, int nitermax)
Computes (member omega_field ).
double funct_omega(double omeg) const
Evaluates , where F is the function defining the rotation profile.
void fait_prim_field()
Computes the member prim_field from omega_field .
double omega_min
Minimum value of .
double(* primfrot)(double, const Tbl &)
Primitive of the function , which vanishes at the stellar center.
Tbl par_frot
Parameters of the function .
double omega_max
Maximum value of .
Scalar psi4
Conformal factor .
Metric gamma
3-metric
Definition star.h:235
Map & mp
Mapping associated with the star.
Definition star.h:180
int nzet
Number of domains of *mp occupied by the star.
Definition star.h:183
Vector beta
Shift vector.
Definition star.h:228
Basic array class.
Definition tbl.h:161
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