/* * Computes the equilibrium configuration of a binary system. * */ /* * Copyright (c) 2000-2003 Eric Gourgoulhon * Copyright (c) 2001-2002 Keisuke Taniguchi * * 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 * */ /* * $Id: coal.C,v 1.17 2016/12/05 16:18:23 j_novak Exp $ * $Log: coal.C,v $ * Revision 1.17 2016/12/05 16:18:23 j_novak * Suppression of some global variables (file names, loch, ...) to prevent redefinitions * * Revision 1.16 2014/10/13 08:53:54 j_novak * Lorene classes and functions now belong to the namespace Lorene. * * Revision 1.15 2014/10/06 15:09:42 j_novak * Modified #include directives to use c++ syntax. * * Revision 1.14 2008/11/14 13:53:23 e_gourgoulhon * Introduced the arrays ent_limit to force the enthalpy values at the * boundaries between the domains inside the stars. * * Revision 1.13 2004/09/28 15:53:25 f_limousin * Improve the rescaling of the domains for nzone = 4 and nzone = 5. * * Revision 1.12 2004/03/25 12:35:36 j_novak * now using namespace Unites * * Revision 1.11 2003/10/24 13:28:51 e_gourgoulhon * Suppressed des_explorer. * * Revision 1.10 2003/10/24 13:26:39 e_gourgoulhon * Corrected error on reduce_shift (previously: it was read and not used !). * * Revision 1.9 2003/09/18 15:06:59 e_gourgoulhon * Added printing of date and host name in the header of the file calcul.d. * * Revision 1.8 2003/09/17 15:40:33 e_gourgoulhon * Test of good opening of the initial data file. * * Revision 1.7 2003/09/16 13:36:50 e_gourgoulhon * -- initial analytical shift decreased by a factor reduce_shift which * can be specified in the parcoal.d file (default value = 0.6) * -- replaced the fpar.getline(blabla, 120) by fpar.ignore(1000,'\n') * -- new output file "resformat.d" at the end of the computation for * formatted output of global quantities (to be read by external code) * * Revision 1.6 2003/09/08 12:06:20 e_gourgoulhon * Added the printing of the virial errors in the output file "calcul.d". * * Revision 1.5 2003/02/20 15:08:42 e_gourgoulhon * Suppressed the qualifier ios::nocreate in call to fstream::open * (not supported by gcc 3.2). * * Revision 1.4 2003/01/09 11:07:48 j_novak * headcpp.h is now compliant with C++ norm. * The include files have been ordered, as well as the local_settings_linux * * Revision 1.3 2002/12/11 13:53:01 k_taniguchi * Add "thres_adapt" to star 2, * and modify the part of "resize". * * Revision 1.2 2001/12/06 16:20:14 e_gourgoulhon * Return type of main changed from 'void' to 'int' * * Revision 1.1.1.1 2001/11/20 15:19:31 e_gourgoulhon * LORENE * * Revision 2.18 2001/08/07 09:51:33 keisuke * Change of the argument in Etoile_bin::equilibrium. * Addition of the procedure to calculate a factor for "resize". * * Revision 2.17 2001/03/20 16:03:22 keisuke * Translation of the origin of the absolute frame. * * Revision 2.16 2000/07/07 14:18:01 eric * Sortie des quantites en unites polytropiques dans calcul.d * (appel de la nouvelle fonction Binaire::display_poly). * * Revision 2.15 2000/07/06 10:06:35 eric * Ajout des parametres et de l'identification dans calcul.d * Creation du fichier identif.d pour l'identification du code. * * Revision 2.14 2000/05/25 15:19:55 eric * Adaptation du mapping gele (cusp). * * Revision 2.13 2000/03/29 09:13:35 eric * Suppression des appels aux routines de dessin (car probleme sur mesioc). * * Revision 2.12 2000/03/22 16:42:22 eric * Sortie des erreurs sur les equations de Poisson. * * Revision 2.11 2000/03/22 14:00:44 eric * Les sorties pour Explorer ne sont effectuees que si graph = 1. * * Revision 2.10 2000/03/22 13:56:41 eric * Changement de prototypage d'Etoile_bin:equilibrium: on passe desormais * le Tbl differ pour avoir en sortie des estimateurs d'erreur. * Reorganisation des fichiers log : introduction de resconv1 et resconv2 * * Revision 2.9 2000/03/22 11:08:00 eric * Reintroduction des messages Log RCS des versions 2.4-2.6 qui avaient * disparus. * / * * Revision 2.8 2000/03/20 14:25:33 eric * resviriel.d rebaptise resconver.d * Theoreme du viriel calcule a la fin et sorti dans calcul.d * Suppression de l'ecriture du mapping dans calcul.d * * Revision 2.7 2000/03/17 16:52:50 eric * Rajout des headers Id et Log qui avaient disparus. * * * Revision 2.6 * Mutliplication du shift analytique par un facteur 0.6 * Omega_initial donne par Binaire::analytical_omega(). * * Revision 2.5 * Introduction d'un shift analytique dans le cas ou le shift initial est nul * * Revision 2.4 * Sortie des log(diff_ent) et log(diff_masse). * * Revision 2.3 2000/02/22 16:58:15 eric * Ajout de la convergence vers une masse baryonique donnee. * * Revision 2.2 2000/02/21 16:47:27 eric * Appel de fait_d_psi avant hydro_euler. * Suppression dessin surface. * * Revision 2.1 2000/02/17 19:54:44 eric * Premiere version operationnelle ! * * Revision 2.0 2000/02/12 18:41:28 eric * *** empty log message *** * * * $Header: /cvsroot/Lorene/Codes/Bin_star/coal.C,v 1.17 2016/12/05 16:18:23 j_novak Exp $ * */ // headers C #include #include #include // headers Lorene #include "binaire.h" #include "eos.h" #include "utilitaires.h" #include "graphique.h" #include "unites.h" namespace Lorene { // Local prototype Cmp raccord_c1(const Cmp& uu, int l1) ; } //****************************************************************************** using namespace Lorene ; int main(){ using namespace Unites ; // Identification of all the subroutines called by the code : // system("ident coal > identif.d") ; // For the display : char display_bold[]="x[1m" ; display_bold[0] = 27 ; char display_normal[] = "x[0m" ; display_normal[0] = 27 ; // To avoid some compilation warnings if (display_bold == 0x0) { cout << qpig << f_unit << mevpfm3 << endl ; } //------------------------------------------------------------------ // Parameters of the computation //------------------------------------------------------------------ char nomini[80] ; int mermax, mermax_eqb, prompt, graph, fmer_stop, fmer_save, mermax_poisson ; int mermax_potvit, mer_masse, fmer_upd_met, ind_rel_met, nofsteps_radius, nofsteps_mass,nmass_delay ; int massscan_rstep; double seuil, relax_poisson, relax_potvit, relax, aexp_masse ; double mbar_voulue_init[2], mbar_voulue_final[2],fact_separ, relax_met, relax_omeg ; double fact_omeg_min, fact_omeg_max, thres_adapt[2], reduce_shift ; double start_dist, radius_step ; ifstream fpar("parcoal_seq.d") ; if ( !fpar.good() ) { cout << "Problem with opening the file parcoal.d ! " << endl ; abort() ; } fpar.ignore(1000, '\n') ; fpar.ignore(1000, '\n') ; fpar.getline(nomini, 80) ; // fpar >> fact_separ ; fpar.ignore(1000, '\n'); fpar >> start_dist ; fpar.ignore(1000, '\n'); fpar >> radius_step ; fpar.ignore(1000, '\n'); fpar >> nofsteps_radius ; fpar.ignore(1000, '\n'); fpar >> mbar_voulue_init[0] ; fpar.ignore(1000, '\n') ; fpar >> mbar_voulue_init[1] ; fpar.ignore(1000, '\n') ; fpar >> mbar_voulue_final[0] ; fpar.ignore(1000, '\n') ; fpar >> mbar_voulue_final[1] ; fpar.ignore(1000, '\n') ; mbar_voulue_init[0] *= msol ; mbar_voulue_init[1] *= msol ; mbar_voulue_final[0] *= msol ; mbar_voulue_final[1] *= msol ; fpar >> nofsteps_mass ; fpar.ignore(1000, '\n') ; fpar >> nmass_delay ; fpar.ignore(1000, '\n') ; fpar >> massscan_rstep ; fpar.ignore(1000, '\n') ; fpar.ignore(1000, '\n') ; fpar >> mermax ; fpar.ignore(1000, '\n') ; fpar >> relax ; fpar.ignore(1000, '\n') ; fpar >> mermax_eqb ; fpar.ignore(1000, '\n') ; fpar >> prompt ; fpar.ignore(1000, '\n') ; fpar >> graph ; fpar.ignore(1000, '\n') ; fpar >> seuil ; fpar.ignore(1000, '\n') ; fpar >> fmer_stop ; fpar.ignore(1000, '\n') ; fpar >> fmer_save ; fpar.ignore(1000, '\n') ; fpar >> mermax_poisson ; fpar.ignore(1000, '\n') ; fpar >> relax_poisson ; fpar.ignore(1000, '\n') ; fpar >> mermax_potvit ; fpar.ignore(1000, '\n') ; fpar >> relax_potvit ; fpar.ignore(1000, '\n') ; fpar >> mer_masse ; fpar.ignore(1000, '\n') ; fpar >> aexp_masse ; fpar.ignore(1000, '\n') ; fpar >> fmer_upd_met ; fpar.ignore(1000, '\n'); fpar >> ind_rel_met ; fpar.ignore(1000, '\n'); fpar >> relax_met ; fpar.ignore(1000, '\n'); if (ind_rel_met == 0) relax_met = 1. ; fpar >> relax_omeg ; fpar.ignore(1000, '\n'); fpar >> fact_omeg_min ; fpar.ignore(1000, '\n'); fpar >> fact_omeg_max ; fpar.ignore(1000, '\n'); fpar >> thres_adapt[0] ; fpar.ignore(1000, '\n'); fpar >> thres_adapt[1] ; fpar.ignore(1000, '\n'); fpar >> reduce_shift ; if ( ! fpar.good() ) { // to ensure compatibility with old reduce_shift = 0.6 ; // parcoal.d files which did not had } // the reduce_shift line fpar.close() ; cout << endl << "==========================================================" << endl << " Physical parameters " << endl << "==========================================================" << endl ; cout << endl << endl ; cout << "File containing the initial conditions : " << nomini << endl ; cout << "Factor by which the initial separation will be multiplied : " << start_dist << endl ; if ( abs(mer_masse) < mermax ) { cout << "Initial Baryon mass required for star 1 [M_sol] : " << mbar_voulue_init[0] / msol << endl ; cout << "Initial Baryon mass required for star 2 [M_sol] : " << mbar_voulue_init[1] / msol << endl ; cout << "Final Baryon mass required for star 1 [M_sol] : " << mbar_voulue_final[0] / msol << endl ; cout << "Final Baryon mass required for star 2 [M_sol] : " << mbar_voulue_final[1] / msol << endl ; } cout << endl << "==========================================================" << endl << " Parameters of the computation " << endl << "==========================================================" << endl ; cout << "Maximum number of steps in the main iteration : " << mermax << endl ; cout << "Relaxation factor in the main iteration : " << relax << endl ; cout << "Maximum number of steps in Etoile_bin::equilibrium : " << mermax_eqb << endl ; cout << "Threshold on the enthalpy relative change for ending the computation : " << seuil << endl ; cout << "Step interval between safeguards of the whole configuration : " << fmer_save << endl ; cout << "Maximum number of steps in Map_et::poisson : " << mermax_poisson << endl ; cout << "Relaxation factor in Map_et::poisson : " << relax_poisson << endl ; cout << "Maximum number of steps in Map_radial::poisson_compact : " << mermax_potvit << endl ; cout << "Relaxation factor in Map_radial::poisson_compact : " << relax_potvit << endl ; cout << "Step from which the baryon mass is forced to converge : " << mer_masse << endl ; cout << "Exponent for the increase factor of the central enthalpy : " << aexp_masse << endl ; cout << "Step interval between metric updates : " << fmer_upd_met << endl ; if (ind_rel_met == 1) { cout << "Relaxation factor of the metric : " << relax_met << endl ; } else { cout << "No relaxation on the metric" << endl ; } cout << "Relaxation factor on Omega (orbital angular velocity) : " << relax_omeg << endl ; cout << "Relative low bound in the omega search : " << fact_omeg_min << endl ; cout << "Relative high bound in the omega search : " << fact_omeg_max << endl ; cout << "Threshold on |dH/dr|_eq / |dH/dr|_pole for the adaptation of the mapping for star 1" << endl << thres_adapt[0] << endl ; cout << "Threshold on |dH/dr|_eq / |dH/dr|_pole for the adaptation of the mapping for star 2" << endl << thres_adapt[1] << endl ; cout << "Factor by which the initial analytical shift is reduced : " << reduce_shift << endl ; arrete(prompt) ; //------------------------------------------------------------------ // Read of the initial conditions //------------------------------------------------------------------ FILE* fich = fopen(nomini, "r") ; if (fich == 0x0) { cout << "Problem in opening the file " << nomini << " ! " << endl ; perror(" reason") ; abort() ; } int mer_ini ; fread(&mer_ini, sizeof(int), 1, fich) ; Mg3d mg1(fich) ; Map_et mp1(mg1, fich) ; Eos* peos1 = Eos::eos_from_file(fich) ; Mg3d mg2(fich) ; Map_et mp2(mg2, fich) ; Eos* peos2 = Eos::eos_from_file(fich) ; Binaire star(mp1, *peos1, mp2, *peos2, fich) ; fclose(fich) ; //Tables with values of enthalpy at the borders of domains int nzet1 = star(1).get_nzet() ; int nzet2 = star(2).get_nzet() ; Tbl ent_limit1(nzet1) ; Tbl ent_limit2(nzet2) ; ent_limit1.set_etat_qcq() ; ent_limit2.set_etat_qcq() ; for(int j=0; jsauve(fresu) ; star(1).get_mp().sauve(fresu) ; star(1).get_eos().sauve(fresu) ; star(2).get_mp().get_mg()->sauve(fresu) ; star(2).get_mp().sauve(fresu) ; star(2).get_eos().sauve(fresu) ; star.sauve(fresu) ; fclose(fresu) ; //## cout << endl << "==========================================================" << endl << " Initial conditions " << endl << "==========================================================" << endl ; cout << star << endl ; if ( star(1).is_relativistic() ) { cout << "========================" << endl ; cout << "Relativistic computation" << endl ; cout << "========================" << endl ; } else { cout << "=====================" << endl ; cout << "Newtonian computation" << endl ; cout << "=====================" << endl ; } arrete(prompt) ; //---------------------------------------------------------------- // Auxiliary quantities //---------------------------------------------------------------- double ent_c[2] ; // Central enthalpy in each star double dentdx[2] ; // Central d/dx(enthalpy) in each star // Resizing factor of the first shell // Computation is shown just befor "equilibrium" Tbl fact_resize[] = {Tbl(2), Tbl(2)} ; fact_resize[0].set_etat_qcq() ; fact_resize[1].set_etat_qcq() ; // Error indicators in each star Tbl differ[] = {Tbl(7), Tbl(7)} ; differ[0].set_etat_qcq() ; differ[1].set_etat_qcq() ; for (int i=1 ; i<=2 ; i++) { ent_c[i-1] = star(i).get_ent()()(0, 0, 0, 0) ; differ[i-1].set(0) = 1 ; // diff_ent = 1 differ[i-1].set(1) = 1 ; // err_psi = 1 differ[i-1].set(2) = 1 ; // err_logn = 1 differ[i-1].set(3) = 1 ; // err_beta = 1 differ[i-1].set(4) = 1 ; // err_shift_x = 1 differ[i-1].set(5) = 1 ; // err_shift_y = 1 differ[i-1].set(6) = 1 ; // err_shift_z = 1 } double relax_jm1 = 1. - relax ; double relax_omeg_jm1 = 1. - relax_omeg ; //---------------------------------------------------------------- // Binary system at the previous step (for the relaxation) //---------------------------------------------------------------- Binaire star_jm1 = star ; double omega_jm1 = star_jm1.get_omega() ; // Notes step 9, seems to help to only call this the first time // logn_comp and pot_centri are initialized to 0 on star_jm1 : // --------------------------------------------------------- if(loop1 == 1) { for (int i=1 ; i<=2 ; i++) { star_jm1.set(i).set_logn_comp() = 0 ; star_jm1.set(i).set_pot_centri() = 0 ; } } //---------------------------------------------------------------- // Openning of log files //---------------------------------------------------------------- sprintf(name, "resglob_%e_%e_%e.d", distance,massfilename[0],massfilename[1]) ; ofstream fichresu(name) ; fichresu.precision(16) ; sprintf(name, "resrota_%e_%e_%e.d", distance,massfilename[0],massfilename[1]) ; ofstream fichrota(name) ; fichrota.precision(16) ; sprintf(name, "resdiffm_%e_%e_%e.d", distance,massfilename[0],massfilename[1]) ; ofstream fichvir(name) ; fichvir.precision(16) ; ofstream fichconv[2] ; sprintf(name, "resconv1_%e_%e_%e.d", distance,massfilename[0],massfilename[1]) ; fichconv[0].open(name) ; fichconv[0].precision(16) ; sprintf(name, "resconv2_%e_%e_%e.d", distance,massfilename[0],massfilename[1]) ; fichconv[1].open(name) ; fichconv[1].precision(16) ; ofstream fichet[2] ; sprintf(name, "resstar1_%e_%e_%e.d", distance,massfilename[0],massfilename[1]) ; fichet[0].open(name) ; fichet[0].precision(16) ; sprintf(name, "resstar2_%e_%e_%e.d", distance,massfilename[0],massfilename[1]) ; fichet[1].open(name) ; fichet[1].precision(16) ; fichrota << "# Omega [rad/s] x_axe [km] x_g (et 0) [km] x_g (et 1) [km] M_grav [M_sol] J [ G M_sol^2 / c]" << endl ; fichvir << "# diff_mass" << endl ; for (int i=1; i<=2; i++) { fichconv[i-1] << "# diff_ent err_psi err_logn err_beta err_shift_x err_shift_y err_shift_z" << endl ; fichet[i-1] << "# ori_x [km] ent_c M_bar [M_sol] R(theta=0) [km] R(pi/2, 0) [km] R(pi/2, pi/2) [km] R(pi/2, pi) [km]" << endl ; } double omega_kep, diff_mass ; int mer ; //============================================================================ // Start of iteration //============================================================================ //Don't stop until we reach the final mass!!! for (mer=0; (differ[0](0) > seuil || ((dist == massscan_rstep-1) && (mer <= mer_masse + nofsteps_mass + nmass_delay))) && (mer < mermax); mer++) { double mbar_voulue[2]; // this is the number of steps from the point we start scanning int nmasscount = mer - mer_masse - nmass_delay; if(dist < massscan_rstep - 1) { mbar_voulue[0] = mbar_voulue_init[0]; mbar_voulue[1] = mbar_voulue_init[1]; } else if (dist > massscan_rstep ) { mbar_voulue[0] = mbar_voulue_final[0]; mbar_voulue[1] = mbar_voulue_final[1]; } else { //dist == massscan_rstep - 1; we need to scan the mass! if(nmasscount < 0) { mbar_voulue[0] = mbar_voulue_init[0]; mbar_voulue[1] = mbar_voulue_init[1]; } else if (nmasscount > nofsteps_mass) { mbar_voulue[0] = mbar_voulue_final[0]; mbar_voulue[1] = mbar_voulue_final[1]; } else { mbar_voulue[0] = mbar_voulue_init[0] + nmasscount * (mbar_voulue_final[0] - mbar_voulue_init[0]) / (1.0*nofsteps_mass); mbar_voulue[1] = mbar_voulue_init[1] + nmasscount * (mbar_voulue_final[1] - mbar_voulue_init[1]) / (1.0*nofsteps_mass); } } cout<<"Baryon masses:"<(2) : " << differ[0](0) << " " << differ[1](0) << " " << differ[0](0) - differ[1](0) << endl ; cout << "==========================================================================" << endl ; fichresu << mer; fichresu << " step" << endl ; fichrota << mer ; fichvir << mer ; fichconv[0] << mer ; fichconv[1] << mer ; fichet[0] << mer ; fichet[1] << mer ; //------------------------------------------------------------------ // Computation of the metric coefficients //------------------------------------------------------------------ if ( (mer % fmer_upd_met) == 0 ) { for (int i=1; i<=2; i++) { (star.set(i)).update_metric(star(3-i), star_jm1(i), relax_met) ; } for (int i=1; i<=2; i++) { (star.set(i)).update_metric_der_comp(star(3-i)) ; } } //------------------------------------------------------------------ // Computation of the orbital angular velocity Omega //------------------------------------------------------------------ double xgg[2] ; star.orbit(fact_omeg_min, fact_omeg_max, xgg[0], xgg[1]) ; // Translation of the stars in order to set the origin // of the absolute frame on the rotation axis //----------------------------------------------------- double x_rot = star.get_x_axe() ; for (int i=1 ; i<=2 ; i++) { double ori_x_old = (star(i).get_mp()).get_ori_x() ; double ori_x_new = ori_x_old - x_rot ; ((star.set(i)).set_mp()).set_ori(ori_x_new, 0., 0.) ; } star.set_x_axe() = 0. ; // Relaxation on the orbital velocity // ---------------------------------- double omega_j = star.get_omega() ; omega_j = relax_omeg * omega_j + relax_omeg_jm1 * omega_jm1 ; omega_jm1 = omega_j ; star.set_omega() = omega_j ; cout << display_bold << "New orbital velocity Omega : " << star.get_omega() * f_unit << " rad/s" << display_normal << endl ; // Keplerian velocity (for comparison only) // ---------------------------------------- omega_kep = sqrt( g_si/g_unit * (star(1).mass_g() + star(2).mass_g()) / pow( star.separation(), 3.) ) ; cout << "``Keplerian'' velocity (for comparison only) : " << omega_kep * f_unit << " rad/s" << endl ; cout << "New X coordinate of the rotation axis : " << star.get_x_axe() / km << " km" << endl ; arrete(prompt) ; fichresu << star.get_x_axe() / km ; fichresu << " abscidia of the rotation axis [km] " << endl ; fichresu << star.get_omega() * f_unit ; fichresu << " Orbital frequency Omega [rad/s] " << endl ; fichresu << xgg[0] / km ; fichresu << " Abscidia ``center of mass'' star 1 [km] " << endl ; fichresu << xgg[1] / km ; fichresu << " Abscidia ``center of mass'' star 2 [km] " << endl ; fichrota << " " << star.get_omega() * f_unit ; fichrota << " " << star.get_x_axe() / km ; fichrota << " " << xgg[0] / km ; fichrota << " " << xgg[1] / km ; //------------------------------------------------------------------ // Computation of B^i/N (bsn) and pot_centri in each star //------------------------------------------------------------------ for (int i=1; i<=2; i++) { (star.set(i)).kinematics( star.get_omega(), star.get_x_axe() ) ; } //------------------------------------------------------------------ // Computation of gam_euler, u_euler, ener_euler, s_euler, // wit_w and loggam in each star //------------------------------------------------------------------ for (int i=1; i<=2; i++) { (star.set(i)).fait_d_psi() ; (star.set(i)).hydro_euler() ; // Check of the Binaire::orbit computation //---------------------------------------- Cmp tmp = star(i).get_logn_auto()() + star(i).get_logn_comp()() + star(i).get_loggam()() ; double grad1 = tmp.dsdx()(0, 0, 0, 0) ; double grad2 = star(i).get_pot_centri()().dsdx()(0, 0, 0, 0) ; dentdx[i-1] = star(i).get_ent()().dsdx()(0, 0, 0, 0) ; double grad3 = star(i).get_loggam()().dsdx()(0, 0, 0, 0) ; cout << "Star " << i << " : " << endl ; cout << " central dH/dx : " << dentdx[i-1] << endl ; cout << " central d(log(Gam))/dx : " << grad3 << endl ; cout << " central d/dx(nu + log(Gam)) : " << grad1 << endl ; cout << " central d/dx(pot_centri) : " << grad2 << endl ; cout << " central d/dx(nu + log(Gam) + pot_centri) : " << grad1 + grad2 << endl ; } //------------------------------------------------------------------ // Computation of the stellar equilibrium configurations //------------------------------------------------------------------ for (int i=1; i<=2; i++) { // Computation of the resizing factor double ray_eq_auto = star(i).ray_eq() ; double ray_eq_comp = star(3-i).ray_eq() ; double ray_eq_pi_comp = star(3-i).ray_eq_pi() ; int num_resize ; if (mg1.get_nzone() > 3) { assert( mg2.get_nzone() == mg1.get_nzone() ) ; num_resize = mg1.get_nzone() - 3 ; } else { num_resize = star(i).get_nzet() ; } double lambda_resize = 0.95 * (star.separation() - ray_eq_comp)/ray_eq_auto ; fact_resize[i-1].set(0) = (lambda_resize < 2.*num_resize) ? lambda_resize : 2.*num_resize ; fact_resize[i-1].set(1) = 1.05 * (star.separation() + ray_eq_pi_comp)/ray_eq_auto ; } for (int i=1; i<=2; i++) { // Relaxation on logn_comp (only if it has not been done by // update_metric) // -------------------------------------------------------- if ( (ind_rel_met == 0) || ( (mer % fmer_upd_met) != 0 ) ) { star.set(i).set_logn_comp() = relax * star(i).get_logn_comp() + relax_jm1 * star_jm1(i).get_logn_comp() ; } // Relaxation on pot_centri // ------------------------ star.set(i).set_pot_centri() = relax * star(i).get_pot_centri() + relax_jm1 * star_jm1(i).get_pot_centri() ; // Call to Etoile_bin::equilibrium // -------------------------------- //mbtest (star.set(i)).equilibrium(ent_c[i-1], mermax_eqb, mermax_poisson, relax_poisson, mermax_potvit, relax_potvit, thres_adapt[i-1], fact_resize[i-1], differ[i-1], pent_limit[i-1]) ; } //------------------------------------------------------------------ // Relaxations //------------------------------------------------------------------ for (int i=1; i<=2; i++) { star.set(i).relaxation( star_jm1(i), relax, relax_met, mer, fmer_upd_met ) ; star.set(i).hydro_euler() ; } /* if (mer % 1000 == 0) { double r_max = 1.2 * star(1).ray_eq() ; des_profile(star(1).get_nbar()(), 0., r_max, M_PI/2, 0., "n", "Baryon density") ; des_profile(star(1).get_ener()(), 0., r_max, M_PI/2, 0., "e", "Energy density") ; des_profile(star(1).get_press()(), 0., r_max, M_PI/2, 0., "p", "Pressure") ; des_profile(star(1).get_ent()(), 0., r_max, M_PI/2, 0., "H", "Enthalpy") ; } */ //------------------------------------------------------------------ // Change in the central enthalpy to get a fixed baryon mass //------------------------------------------------------------------ if (mer >= mer_masse) { for (int i=1; i<=2; i++) { double xx = star(i).mass_b() / mbar_voulue[i-1] - 1. ; cout << "Discrepancy M_b / wanted M_b : " << xx << endl ; double xprog = ( mer > 2*mer_masse) ? 1. : double(mer-mer_masse)/double(mer_masse) ; xx *= xprog ; double ax = .5 * ( 2. + xx ) / (1. + xx ) ; double fact_ent = pow(ax, aexp_masse) ; cout << " xprog, xx, ax, fact : " << xprog << " " << xx << " " << ax << " " << fact_ent << endl ; ent_c[i-1] *= fact_ent ; } } // Updates for the next step // ------------------------- star_jm1 = star ; cout << star << endl ; // Graphical output // ---------------- if ( (graph==1) && (mer % 5 == 0) ) { double xdes_min = - star(1).ray_eq_pi() + star(1).get_mp().get_ori_x() ; xdes_min *= 1.5 ; double xdes_max = star(2).ray_eq_pi() + star(2).get_mp().get_ori_x() ; xdes_max *= 1.5 ; double ydes_min = - 2.5 * star(1).ray_eq_pis2() ; double ydes_max = 2.5 * star(2).ray_eq_pis2() ; Cmp surf1 = star(1).get_ent()() ; Cmp surf1_ext(mp1) ; surf1_ext = - 0.2 * surf1(0, 0, 0, 0) ; surf1_ext.annule(0, star(1).get_nzet()-1) ; surf1.annule(star(1).get_nzet(), mg1.get_nzone()-1) ; surf1 = surf1 + surf1_ext ; surf1 = raccord_c1(surf1, star(1).get_nzet()) ; Cmp surf2 = star(2).get_ent()() ; Cmp surf2_ext(mp2) ; surf2_ext = - 0.2 * surf2(0, 0, 0, 0) ; surf2_ext.annule(0, star(2).get_nzet()-1) ; surf2.annule(star(2).get_nzet(), mg2.get_nzone()-1) ; surf2 = surf2 + surf2_ext ; surf2 = raccord_c1(surf2, star(2).get_nzet()) ; des_coupe_bin_z(star(1).get_nbar()(), star(2).get_nbar()(), 0., xdes_min, xdes_max, ydes_min, ydes_max, "Baryon density (z=0)", &surf1, &surf2 ) ; des_coupe_bin_y(star(1).get_nbar()(), star(2).get_nbar()(), 0., xdes_min, xdes_max, ydes_min, ydes_max, "Baryon density (y=0)", &surf1, &surf2 ) ; } //----------------------------------------------------------------------- // The whole configuration is saved in a file //----------------------------------------------------------------------- if ( (mer % fmer_save) == 0 ) { sprintf(name, "resu_%e.d", distance) ; FILE* fresu2 = fopen(name, "w") ; fwrite(&mer, sizeof(int), 1, fresu2) ; // mer star(1).get_mp().get_mg()->sauve(fresu2) ; star(1).get_mp().sauve(fresu2) ; star(1).get_eos().sauve(fresu2) ; star(2).get_mp().get_mg()->sauve(fresu2) ; star(2).get_mp().sauve(fresu2) ; star(2).get_eos().sauve(fresu2) ; star.sauve(fresu2) ; fclose(fresu2) ; } //-------------------------------------------- // Writing of global quantities in log files //-------------------------------------------- for (int i=1 ; i<=2 ; i++) { fichresu << star(i).mass_b() / msol ; fichresu << " Baryon mass of star " << i << " [M_sol] " << endl ; fichresu << differ[i-1](0) ; fichresu << " relative variation enth. star " << i << endl ; fichresu << star(i).ray_pole() / km ; fichresu << " R(theta=0) [km] " << endl ; fichresu << star(i).ray_eq() / km ; fichresu << " R(theta=pi/2, phi=0) [km] " << endl ; fichresu << star(i).ray_eq_pis2() / km ; fichresu << " R(theta=pi/2, phi=pi/2) [km] " << endl ; fichresu << star(i).ray_eq_pi() / km ; fichresu << " R(theta=pi/2, phi=pi) [km] " << endl ; fichconv[i-1] << " " << log10( fabs(differ[i-1](0)) + 1.e-16 ) ; fichconv[i-1] << " " << log10( fabs(differ[i-1](1)) + 1.e-16 ) ; fichconv[i-1] << " " << log10( fabs(differ[i-1](2)) + 1.e-16 ) ; fichconv[i-1] << " " << log10( fabs(differ[i-1](3)) + 1.e-16 ) ; fichconv[i-1] << " " << log10( fabs(differ[i-1](4)) + 1.e-16 ) ; fichconv[i-1] << " " << log10( fabs(differ[i-1](5)) + 1.e-16 ) ; fichconv[i-1] << " " << log10( fabs(differ[i-1](6)) + 1.e-16 ) ; fichet[i-1] << " " << star(i).get_mp().get_ori_x() / km ; fichet[i-1] << " " << ent_c[i-1] ; fichet[i-1] << " " << star(i).mass_b() / msol ; fichet[i-1] << " " << star(i).ray_pole() / km ; fichet[i-1] << " " << star(i).ray_eq() / km ; fichet[i-1] << " " << star(i).ray_eq_pis2() / km ; fichet[i-1] << " " << star(i).ray_eq_pi() / km ; } // End of loop on the stars diff_mass = ( star(2).mass_b() - star(1).mass_b() ) / star(1).mass_b() ; cout << "Relative difference between the baryon masses: " << diff_mass << endl ; fichresu << diff_mass ; fichresu << " Relative difference between the baryon masses" << endl ; fichvir << " " << log10( fabs(diff_mass) + 1.e-16 ) ; fichresu << " " << endl ; fichresu.flush() ; fichrota << " " << endl ; fichrota.flush() ; fichvir << " " << endl ; fichvir.flush() ; fichconv[0] << " " << endl ; fichconv[0].flush() ; fichconv[1] << " " << endl ; fichconv[1].flush() ; fichet[0] << " " << endl ; fichet[0].flush() ; fichet[1] << " " << endl ; fichet[1].flush() ; loop1=0; } // End of the main loop (mer) //============================================================================ // End of iteration //============================================================================ fichresu.close() ; fichrota.close() ; fichvir.close() ; fichconv[0].close() ; fichconv[1].close() ; fichet[0].close() ; fichet[1].close() ; //----------------------------------------------- // General features of the final configuration // saved in a file //----------------------------------------------- ofstream fichfinal("calcul.d") ; fichfinal.precision(6) ; time_t rawtime = time(0x0) ; fichfinal << "Date: " << asctime( localtime( &rawtime ) ) << endl ; char* hostname = getenv("HOST") ; if (hostname != 0x0) { fichfinal << "Computer: " << hostname << endl ; } fichfinal << "===================================================================" << endl << endl ; if ( star(1).is_relativistic() ) { fichfinal << "Relativistic computation" ; } else { fichfinal << "Newtonian computation" ; } if ( star(1).is_irrotational() ) { fichfinal << " Irrotational" << endl ; } else { fichfinal << " Co-rotating" << endl ; } fichfinal << star(1).get_eos() << endl ; fichfinal << "Omega = " << star.get_omega() * f_unit << " rad/s" << " Orbital frequency f = " << star.get_omega() / (2*M_PI) * f_unit << " Hz" << endl ; fichfinal << "Omega_kepler = " << omega_kep * f_unit << " rad/s" << endl ; fichfinal << "Coordinate separation : " << star.separation() / km << " km" << endl ; fichfinal << "1/2 ADM mass : " << 0.5 * star.mass_adm() / msol << " Mo" << endl ; fichfinal << "Total angular momentum : " << star.angu_mom()(2)/ ( qpig / (4* M_PI) * msol*msol) << " G M_sol^2 / c" << endl ; cout << "1/2 ADM mass : " << 0.5 * star.mass_adm() / msol << " Mo" << endl ; cout << "Total angular momentum : " << star.angu_mom()(2)/ ( qpig / (4* M_PI) * msol*msol) << " G M_sol^2 / c" << endl ; fichfinal << endl << "Number of steps : " << mer << endl ; for (int i=1 ; i<=2; i++) { fichfinal << endl << "===================================================================" << endl ; fichfinal << " Star no. " << i << endl ; fichfinal << "===================================================================" << endl ; fichfinal << "Grid : " << endl ; fichfinal << "------ " << endl ; fichfinal << *(star(i).get_mp().get_mg()) << endl ; fichfinal << endl << "Physical characteristics : " << endl ; fichfinal << "-------------------------" << endl ; fichfinal << star(i) << endl ; } fichfinal << endl ; star.display_poly(fichfinal) ; // Reduced quantities for polytropic EOS fichfinal << endl << "===================================================================" << endl ; fichfinal << "Diff_ent : star 1 : " << differ[0](0) << " star 2 : " << differ[1](0) << endl ; fichfinal << "Relative difference between the baryon masses of the two stars : " << diff_mass << endl ; fichfinal << "dH/dx at r = 0 : star 1 : " << dentdx[0] << " star 2 : " << dentdx[1] << endl ; if ( star(1).is_relativistic() ) { fichfinal << "Relative error on the virial theorem : " << endl ; fichfinal << " VE(M)= " << star.virial() << " VE(GB)= "<< star.virial_gb() << " VE(FUS)= " << star.virial_fus() << endl ; } else { fichfinal << "Relative error on the virial theorem : " << star.virial() << endl ; } fichfinal << endl << "================================================================" << endl ; fichfinal << " PARAMETERS USED FOR THE COMPUTATION : " << endl ; fichfinal << "================================================================" << endl ; fichfinal.close() ; system("cat parcoal.d >> calcul.d") ; // Identification du code et de ses sous-routines (no. de version RCS) : fichfinal.open("calcul.d", ios::app ) ; fichfinal << endl << "================================================================" << endl ; fichfinal << " IDENTIFICATION OF THE CODE : " << endl ; fichfinal << "================================================================" << endl ; fichfinal.close() ; system("ident coal >> calcul.d") ; // Preparation for CPU infos printing : fichfinal.open("calcul.d", ios::app ) ; fichfinal << endl << "================================================================" << endl ; fichfinal << " CPU TIME and MEMORY infos : " << endl ; fichfinal << "================================================================" << endl ; fichfinal << endl ; fichfinal.close() ; //----------------------------------------------- // General features of the final configuration // saved in a file with scientific notation and // 14 digits for further reading by a code //----------------------------------------------- sprintf(name, "resformat_%e.d", distance) ; ofstream seqfich(name) ; if ( !seqfich.good() ) { cout << "coal : problem with opening the file resformat.d !" << endl ; abort() ; } star.write_global(seqfich) ; seqfich.close() ; if (graph == 1) { //## double xdes_min = - star(1).ray_eq_pi() + star(1).get_mp().get_ori_x() ; // xdes_min *= 1.5 ; // double xdes_max = star(2).ray_eq_pi() + star(2).get_mp().get_ori_x() ; // xdes_max *= 1.5 ; // double ydes_min = - 2.5 * star(1).ray_eq_pis2() ; // double ydes_max = 2.5 * star(2).ray_eq_pis2() ; // // des_coupe_bin_z(star(1).get_nbar()(), star(2).get_nbar()(), 0, // xdes_min, xdes_max, ydes_min, ydes_max, "density (z=0)", // &(star(1).get_ent()()), &(star(2).get_ent()()) ) ; // // des_coupe_bin_z(star(1).get_logn_auto()(), star(2).get_logn_auto()(), 0, // xdes_min, xdes_max, ydes_min, ydes_max, "ln(N) (z=0)", // &(star(1).get_ent()()), &(star(2).get_ent()()) ) ; // // des_coupe_z(star(1).get_logn_comp()(), 0., 1, "logn_comp (z=0)", //## &(star(1).get_ent()()) ) ; } // Cleaning // -------- } delete peos1 ; delete peos2 ; return EXIT_SUCCESS ; }