LORENE: binhor_coal.C Source File

00001 /*
00002  *   Copyright (c) 2004-2005 Francois Limousin
00003  *                           Jose-Luis Jaramillo
00004  *
00005  *   This file is part of LORENE.
00006  *
00007  *   LORENE is free software; you can redistribute it and/or modify
00008  *   it under the terms of the GNU General Public License as published by
00009  *   the Free Software Foundation; either version 2 of the License, or
00010  *   (at your option) any later version.
00011  *
00012  *   LORENE is distributed in the hope that it will be useful,
00013  *   but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  *   GNU General Public License for more details.
00016  *
00017  *   You should have received a copy of the GNU General Public License
00018  *   along with LORENE; if not, write to the Free Software
00019  *   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00020  *
00021  */
00022 
00023 
00024 char binhor_coal_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bin_hor/binhor_coal.C,v 1.13 2007/04/13 15:28:55 f_limousin Exp $" ;
00025 
00026 /*
00027  * $Id: binhor_coal.C,v 1.13 2007/04/13 15:28:55 f_limousin Exp $
00028  * $Log: binhor_coal.C,v $
00029  * Revision 1.13  2007/04/13 15:28:55  f_limousin
00030  * Lots of improvements, generalisation to an arbitrary state of
00031  * rotation, implementation of the spatial metric given by Samaya.
00032  *
00033  * Revision 1.12  2006/08/01 14:37:19  f_limousin
00034  * New version
00035  *
00036  * Revision 1.11  2006/06/28 13:36:09  f_limousin
00037  * Convergence to a given irreductible mass
00038  *
00039  * Revision 1.10  2006/05/24 16:56:37  f_limousin
00040  * Many small modifs.
00041  *
00042  * Revision 1.9  2005/09/13 18:33:15  f_limousin
00043  * New function vv_bound_cart_bin(double) for computing binaries with
00044  * berlin condition for the shift vector.
00045  * Suppress all the symy and asymy in the importations.
00046  *
00047  * Revision 1.8  2005/07/11 08:21:57  f_limousin
00048  * Implementation of a new boundary condition for the lapse in the binary
00049  * case : boundary_nn_Dir_lapl().
00050  *
00051  * Revision 1.7  2005/03/10 17:21:52  f_limousin
00052  * Add the Berlin boundary condition for the shift.
00053  * Some changes to avoid warnings.
00054  *
00055  * Revision 1.6  2005/03/10 17:09:05  f_limousin
00056  * Display the logarithm of viriel and convergence.
00057  *
00058  * Revision 1.5  2005/03/10 16:57:00  f_limousin
00059  * Improve the convergence of the code coal_bh.
00060  *
00061  * Revision 1.4  2005/02/24 17:24:26  f_limousin
00062  * The boundary conditions for psi, N and beta are now parameters in
00063  * par_init.d and par_coal.d.
00064  *
00065  * Revision 1.3  2005/02/07 10:43:36  f_limousin
00066  * Add the printing of the regularisation of the shift in the case N=0
00067  * on the horizon.
00068  *
00069  * Revision 1.2  2004/12/31 15:40:21  f_limousin
00070  * Improve the initialisation of several quantities in set_statiques().
00071  *
00072  * Revision 1.1  2004/12/29 16:11:19  f_limousin
00073  * First version
00074  *
00075  *
00076  * $Header: /cvsroot/Lorene/C++/Source/Bin_hor/binhor_coal.C,v 1.13 2007/04/13 15:28:55 f_limousin Exp $
00077  *
00078  */
00079 
00080 //standard
00081 #include <stdlib.h>
00082 
00083 // Lorene
00084 #include "tensor.h"
00085 #include "isol_hor.h"
00086 #include "graphique.h"
00087 
00088 
00089 void Bin_hor::set_statiques (double precis, double relax, int bound_nn,
00090                  double lim_nn, int bound_psi) {
00091     
00092   int nz = hole1.mp.get_mg()->get_nzone() ;
00093     
00094   set_omega(0) ;
00095   hole1.init_met_trK() ;
00096   hole2.init_met_trK() ;
00097   init_bin_hor() ;
00098   extrinsic_curvature() ;
00099       
00100   int indic = 1 ;
00101   int conte = 0 ;
00102  
00103   cout << "Static black holes : " << endl ;
00104   while (indic == 1) {
00105     Scalar lapse_un_old (hole1.n_auto) ;
00106 
00107     solve_psi (precis, relax, bound_psi) ;
00108     solve_lapse (precis, relax, bound_nn, lim_nn) ;
00109 
00110     //  des_profile(hole1.nn(), 0, 20, M_PI/2, M_PI) ;
00111 
00112     double erreur = 0 ;
00113     Tbl diff (diffrelmax (lapse_un_old, hole1.n_auto)) ;
00114     for (int i=1 ; i<nz ; i++)
00115       if (diff(i) > erreur)
00116     erreur = diff(i) ;
00117     
00118     cout << "Step : " << conte << " Difference : " << erreur << endl ;
00119     
00120     if (erreur < precis)
00121       indic = -1 ;
00122     conte ++ ;
00123   }
00124 }
00125 
00126 double Bin_hor::coal (double angu_vel, double relax, int nb_ome,
00127               int nb_it, int bound_nn, double lim_nn, 
00128               int bound_psi, int bound_beta, double omega_eff,
00129               double alpha,
00130               ostream& fich_iteration, ostream& fich_correction,
00131               ostream& fich_viriel, ostream& fich_kss, 
00132               int step, int search_mass, double mass_irr, 
00133               const int sortie) {
00134     
00135   int nz = hole1.mp.get_mg()->get_nzone() ;
00136 
00137   double precis = 1e-7 ;
00138     
00139   // LOOP INCREASING OMEGA  : 
00140   cout << "OMEGA INCREASED STEP BY STEP." << endl ;
00141   double homme = get_omega() ;
00142   double inc_homme = (angu_vel - homme)/nb_ome ;
00143   for (int pas = 0 ; pas <nb_ome ; pas ++) {
00144       
00145     bool verif = false ;
00146     if (omega_eff == alpha*homme ) verif = true ;
00147       
00148     homme += inc_homme ;
00149     set_omega (homme) ;
00150     if (verif)
00151       omega_eff = alpha*homme ;
00152     Scalar beta_un_old (hole1.beta_auto(1)) ;
00153       
00154     solve_shift (precis, relax, bound_beta, omega_eff) ;
00155     extrinsic_curvature() ;
00156 
00157     solve_psi (precis, relax, bound_psi) ;
00158     solve_lapse (precis, relax, bound_nn, lim_nn) ;
00159     
00160     // Convergence to the given irreductible mass 
00161     if (search_mass == 1 && step >= 30) {
00162       double mass_area = sqrt(hole1.area_hor()/16/M_PI) + 
00163     sqrt(hole2.area_hor()/16/M_PI) ;
00164       double error_m = (mass_irr - mass_area) / mass_irr ;
00165       double scaling_r = pow((2-error_m)/(2-2*error_m), 1.) ;
00166       hole1.mp.homothetie_interne(scaling_r) ;
00167       hole1.radius = hole1.radius *scaling_r ;
00168       hole2.mp.homothetie_interne(scaling_r) ;
00169       hole2.radius = hole2.radius *scaling_r ;
00170     
00171       // Update of the different metrics (another possibility would 
00172       // be to set all derived quantities to 0x0, especially
00173       // the connection p_connect
00174       hole1.ff = hole1.mp.flat_met_spher() ;
00175       hole1.tgam = hole1.mp.flat_met_spher() ;
00176       hole2.ff = hole2.mp.flat_met_spher() ;
00177       hole2.tgam = hole1.mp.flat_met_spher() ;
00178     
00179     }
00180       
00181     cout << "Angular momentum computed at the horizon : " << ang_mom_hor()
00182      << endl ;
00183       
00184     double erreur = 0 ;
00185     Tbl diff (diffrelmax (beta_un_old, hole1.beta_auto(1))) ;
00186     for (int i=1 ; i<nz ; i++)
00187       if (diff(i) > erreur)
00188     erreur = diff(i) ;
00189       
00190     // Saving ok K_{ij}s^is^j
00191     // -----------------------
00192     
00193     Scalar kkss (contract(hole1.get_k_dd(), 0, 1, 
00194               hole1.get_gam().radial_vect()*
00195               hole1.get_gam().radial_vect(), 0, 1)) ;
00196     double max_kss = kkss.val_grid_point(1, 0, 0, 0) ;
00197     double min_kss = kkss.val_grid_point(1, 0, 0, 0) ;
00198     int nnp = hole1.mp.get_mg()->get_np(1) ;
00199     int nnt = hole2.mp.get_mg()->get_nt(1) ;
00200     for (int k=0 ; k<nnp ; k++)
00201       for (int j=0 ; j<nnt ; j++){
00202     if (kkss.val_grid_point(1, k, j, 0) > max_kss)
00203       max_kss = kkss.val_grid_point(1, k, j, 0) ;
00204     if (kkss.val_grid_point(1, k, j, 0) < min_kss)
00205       min_kss = kkss.val_grid_point(1, k, j, 0) ;
00206       }
00207 
00208     if (sortie != 0) {
00209       fich_iteration << step << " " << log10(erreur) << " " << homme << endl ;
00210       fich_correction << step << " " << log10(hole1.regul) << " " << homme << endl ;
00211       //      fich_viriel << step << " " << log10(fabs(viriel())) << " " << homme << endl ;
00212       fich_viriel << step << " " << viriel() << " " << homme << " " << hole1.omega_hor() - alpha*homme << " " << omega_eff << endl ;
00213       fich_kss << step << " " << max_kss << " " << min_kss << endl ;
00214     }
00215         
00216     cout << "STEP : " << step << " DIFFERENCE : " << erreur << endl ;
00217     step ++ ;
00218   }
00219     
00220   // LOOP WITH FIXED OMEGA :
00221 
00222   if (nb_it !=0)
00223     cout << "OMEGA FIXED" << endl ;
00224   double erreur ;
00225 
00226   for (int pas = 0 ; pas <nb_it ; pas ++) {
00227     
00228     Scalar beta_un_old (hole1.beta_auto(1)) ;
00229 
00230     solve_shift (precis, relax, bound_beta, omega_eff) ;
00231     extrinsic_curvature() ;
00232 
00233     solve_psi (precis, relax, bound_psi) ;
00234     solve_lapse (precis, relax, bound_nn, lim_nn) ;
00235 
00236     // Convergence to the given irreductible mass 
00237     if (search_mass == 1 && step >= 30) {
00238       double mass_area = sqrt(hole1.area_hor()/16/M_PI) + 
00239     sqrt(hole2.area_hor()/16/M_PI) ;
00240       double error_m = (mass_irr - mass_area) / mass_irr ;
00241       double scaling_r = pow((2-error_m)/(2-2*error_m), 1.) ;
00242       
00243       hole1.mp.homothetie_interne(scaling_r) ;
00244       hole1.radius = hole1.radius *scaling_r ;
00245       hole2.mp.homothetie_interne(scaling_r) ;
00246       hole2.radius = hole2.radius *scaling_r ;
00247     }
00248 
00249     erreur = 0 ;
00250     Tbl diff (diffrelmax (beta_un_old, hole1.beta_auto(1))) ;
00251     for (int i=1 ; i<nz ; i++)
00252       if (diff(i) > erreur)
00253     erreur = diff(i) ;
00254 
00255     // Saving ok K_{ij}s^is^j
00256     // -----------------------
00257     
00258     Scalar kkss (contract(hole1.get_k_dd(), 0, 1, 
00259               hole1.get_gam().radial_vect()*
00260               hole1.get_gam().radial_vect(), 0, 1)) ;
00261     double max_kss = kkss.val_grid_point(1, 0, 0, 0) ;
00262     double min_kss = kkss.val_grid_point(1, 0, 0, 0) ;
00263     int nnp = hole1.mp.get_mg()->get_np(1) ;
00264     int nnt = hole2.mp.get_mg()->get_nt(1) ;
00265     for (int k=0 ; k<nnp ; k++)
00266       for (int j=0 ; j<nnt ; j++){
00267     if (kkss.val_grid_point(1, k, j, 0) > max_kss)
00268       max_kss = kkss.val_grid_point(1, k, j, 0) ;
00269     if (kkss.val_grid_point(1, k, j, 0) < min_kss)
00270       min_kss = kkss.val_grid_point(1, k, j, 0) ;
00271       }
00272 
00273     
00274     if (sortie != 0) {
00275       fich_iteration << step << " " << log10(erreur) << " " << homme << endl ;
00276       fich_correction << step << " " << log10(hole1.regul) << " " << homme << endl ;
00277       //      fich_viriel << step << " " << log10(fabs(viriel())) << " " << homme << endl ;
00278       fich_viriel << step << " " << viriel() << " " << homme << " " << hole1.omega_hor() - alpha*homme << " " << omega_eff << endl ;
00279       fich_kss << step << " " << max_kss << " " << min_kss << endl ;
00280     }
00281 
00282     cout << "STEP : " << step << " DIFFERENCE : " << erreur << endl ;
00283     step ++ ;
00284   }
00285 
00286   if (nb_it != 0){
00287     fich_iteration << "#----------------------------"  << endl ;
00288     fich_correction << "#-----------------------------" << endl ;
00289     fich_viriel << "#------------------------------"  << endl ;
00290   }
00291 
00292   return viriel() ;
00293 }