et_rot_diff_faitomeg.C
/*
* Functions Et_rot_diff::fait_omega_field
* Et_rot_diff::fait_prim_field
*
* (see file et_rot_diff.h for documentation).
*
*/
/*
* Copyright (c) 2001 Eric Gourgoulhon
*
* This file is part of LORENE.
*
* LORENE is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* LORENE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LORENE; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
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 $" ;
/*
* $Id: et_rot_diff_faitomeg.C,v 1.3 2014/10/13 08:52:57 j_novak Exp $
* $Log: et_rot_diff_faitomeg.C,v $
* Revision 1.3 2014/10/13 08:52:57 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.2 2003/05/14 20:07:00 e_gourgoulhon
* Suppressed the outputs (cout)
*
* Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
* LORENE
*
* Revision 1.3 2001/10/25 09:43:18 eric
* omega_min est determine dans l'etoile seulement (l<nzet).
*
* Revision 1.2 2001/10/25 09:21:29 eric
* Recherche de Omega dans un intervalle de 20% autour de la valeur precedente.
*
* Revision 1.1 2001/10/19 08:18:23 eric
* Initial revision
*
*
* $Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_diff_faitomeg.C,v 1.3 2014/10/13 08:52:57 j_novak Exp $
*
*/
// Headers Lorene
#include "et_rot_diff.h"
#include "utilitaires.h"
#include "param.h"
namespace Lorene {
double et_rot_diff_fzero(double omeg, const Param& par) ;
//-----------------------------------//
// fait_omega_field //
//-----------------------------------//
void Et_rot_diff::fait_omega_field(double omeg_min, double omeg_max,
double precis, int nitermax) {
const Mg3d& mg = *(mp.get_mg()) ;
int nz = mg.get_nzone() ;
int nzm1 = nz - 1 ;
// Computation of B^2 r^2 sin^2(theta):
// -----------------------------------
Cmp brst2 = bbb() ;
brst2.annule(nzm1) ;
brst2.std_base_scal() ;
brst2.mult_rsint() ; // Multiplication by r sin(theta)
brst2 = brst2*brst2 ;
Cmp nnn2 = nnn() * nnn() ;
Param par ;
par.add_cmp(nnn2, 0) ;
par.add_cmp(brst2, 1) ;
par.add_cmp(nphi(), 2) ;
int l, k, j, i ;
par.add_int(l, 0) ;
par.add_int(k, 1) ;
par.add_int(j, 2) ;
par.add_int(i, 3) ;
par.add_etoile(*this) ;
// Loop on the grid points
// -----------------------
int niter ;
bool prev_zero = (omega_field.get_etat() == ETATZERO) ;
omega_field.allocate_all() ;
// cout << "Et_rot_diff::fait_omega_field: niter : " << endl ;
for (l=0; l<nzet+1; l++) {
Tbl& tom = omega_field.set().set(l) ;
for (k=0; k<mg.get_np(l); k++) {
for (j=0; j<mg.get_nt(l); j++) {
for (i=0; i<mg.get_nr(l); i++) {
double& omeg = tom.set(k, j, i) ;
double omeg_min1, omeg_max1 ;
if ( prev_zero || omeg == double(0)) {
omeg_min1 = omeg_min ;
omeg_max1 = omeg_max ;
}
else{
omeg_min1 = 0.8 * omeg ;
omeg_max1 = 1.2 * omeg ;
}
omeg = zerosec(et_rot_diff_fzero, par, omeg_min1,
omeg_max1, precis, nitermax, niter) ;
// cout << " " << niter ;
}
}
}
// cout << endl ;
}
// omega_field is set to 0 far from the star:
// ---------------------------------------------------
for (l=nzet+1; l<nz; l++) {
omega_field.set().set(l) = 0 ;
}
omega_field.set_std_base() ;
// Min and Max of Omega
// --------------------
omega_min = min( omega_field()(0) ) ;
for (l=1; l<nzet; l++) {
double xx = min( omega_field()(l) ) ;
omega_min = (xx < omega_min) ? xx : omega_min ;
}
omega_max = max( max( omega_field() ) ) ;
// Update of prim_field
// --------------------
fait_prim_field() ;
}
//-----------------------------------//
// et_rot_diff_fzero //
//-----------------------------------//
double et_rot_diff_fzero(double omeg, const Param& par) {
const Cmp& nnn2 = par.get_cmp(0) ;
const Cmp& brst2 = par.get_cmp(1) ;
const Cmp& nphi = par.get_cmp(2) ;
int l = par.get_int(0) ;
int k = par.get_int(1) ;
int j = par.get_int(2) ;
int i = par.get_int(3) ;
const Et_rot_diff* et =
dynamic_cast<const Et_rot_diff*>(&par.get_etoile()) ;
double fom = et->funct_omega(omeg) ;
double omnp = omeg - nphi(l, k, j, i) ;
return fom - brst2(l, k, j, i) * omnp /
( nnn2(l, k, j, i) - brst2(l, k, j, i) * omnp*omnp ) ;
}
//-----------------------------------//
// fait_prim_field //
//-----------------------------------//
void Et_rot_diff::fait_prim_field() {
const Mg3d& mg = *(mp.get_mg()) ;
int nz = mg.get_nzone() ;
// Loop on the grid points in the vicinity of the star
// ---------------------------------------------------
prim_field.allocate_all() ;
for (int l=0; l<nzet+1; l++) {
Tbl& tprim = prim_field.set().set(l) ;
for (int k=0; k<mg.get_np(l); k++) {
for (int j=0; j<mg.get_nt(l); j++) {
for (int i=0; i<mg.get_nr(l); i++) {
tprim.set(k, j, i) =
primfrot( omega_field()(l, k, j, i), par_frot ) ;
}
}
}
}
// prim_field is set to 0 far from the star:
// -----------------------------------------
for (int l=nzet+1; l<nz; l++) {
prim_field.set().set(l) = 0 ;
}
prim_field.set_std_base() ;
}
}
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Lorene::Cmp::mult_rsint
void mult_rsint()
Multiplication by .
Definition cmp_r_manip.C:116
Lorene::Cmp::std_base_scal
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition cmp.C:644
Lorene::Cmp::annule
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:348
Lorene::Et_rot_diff
Class for differentially rotating stars.
Definition et_rot_diff.h:70
Lorene::Et_rot_diff::omega_min
double omega_min
Minimum value of .
Definition et_rot_diff.h:107
Lorene::Et_rot_diff::fait_omega_field
void fait_omega_field(double omeg_min, double omeg_max, double precis, int nitermax)
Computes (member omega_field ).
Definition et_rot_diff_faitomeg.C:73
Lorene::Et_rot_diff::prim_field
Tenseur prim_field
Field .
Definition et_rot_diff.h:111
Lorene::Et_rot_diff::omega_max
double omega_max
Maximum value of .
Definition et_rot_diff.h:108
Lorene::Et_rot_diff::omega_field
Tenseur omega_field
Field .
Definition et_rot_diff.h:105
Lorene::Et_rot_diff::fait_prim_field
void fait_prim_field()
Computes the member prim_field from omga_field .
Definition et_rot_diff_faitomeg.C:199
Lorene::Et_rot_diff::par_frot
Tbl par_frot
Parameters of the function .
Definition et_rot_diff.h:102
Lorene::Et_rot_diff::primfrot
double(* primfrot)(double, const Tbl &)
Primitive of the function , which vanishes at the stellar center.
Definition et_rot_diff.h:93
Lorene::Et_rot_diff::funct_omega
double funct_omega(double omeg) const
Evaluates , where F is the function defining the rotation profile.
Definition et_rot_diff.C:315
Lorene::Etoile_rot::nphi
Tenseur nphi
Metric coefficient .
Definition etoile.h:1510
Lorene::Etoile_rot::bbb
Tenseur bbb
Metric factor B.
Definition etoile.h:1504
Lorene::Etoile::nzet
int nzet
Number of domains of *mp occupied by the star.
Definition etoile.h:432
Lorene::Etoile::nnn
Tenseur nnn
Total lapse function.
Definition etoile.h:509
Lorene::Etoile::mp
Map & mp
Mapping associated with the star.
Definition etoile.h:429
Lorene::Mg3d
Multi-domain grid.
Definition grilles.h:273
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Lorene::Param
Parameter storage.
Definition param.h:125
Lorene::Param::get_int
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition param.C:292
Lorene::Param::add_etoile
void add_etoile(const Etoile &eti, int position=0)
Adds the address of a new Etoile to the list.
Definition param.C:1624
Lorene::Param::add_cmp
void add_cmp(const Cmp &ti, int position=0)
Adds the address of a new Cmp to the list.
Definition param.C:935
Lorene::Param::get_cmp
const Cmp & get_cmp(int position=0) const
Returns the reference of a Cmp stored in the list.
Definition param.C:980
Lorene::Param::add_int
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition param.C:246
Lorene::Param::get_etoile
const Etoile & get_etoile(int position=0) const
Returns the reference of a Etoile stored in the list.
Definition param.C:1669
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::set
double & set(int i)
Read/write of a particular element (index i) (1D case).
Definition tbl.h:281
Lorene::min
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition cmp_math.C:458
Lorene::max
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
Lorene::zerosec
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
Lorene prototypes.
Definition app_hor.h:64