et_bin_hydro.C

/*
 * Methods of the class Etoile_bin for computing hydro quantities
 *
 * (see file etoile.h for documentation)
 */

/*
 *   Copyright (c) 2000-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
 *
 */


 

/*
 * $Id: et_bin_hydro.C,v 1.7 2016/12/05 16:17:53 j_novak Exp $
 * $Log: et_bin_hydro.C,v $
 * Revision 1.7  2016/12/05 16:17:53  j_novak
 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
 *
 * Revision 1.6  2014/10/13 08:52:55  j_novak
 * Lorene classes and functions now belong to the namespace Lorene.
 *
 * Revision 1.5  2005/08/29 15:10:16  p_grandclement
 * Addition of things needed :
 *   1) For BBH with different masses
 *   2) Provisory files for the mixted binaries (Bh and NS) : THIS IS NOT
 *   WORKING YET !!!
 *
 * Revision 1.4  2003/03/03 19:46:09  f_limousin
 * Set standard bases for s_euler.
 *
 * Revision 1.3  2003/01/17 13:34:56  f_limousin
 * Replace A^2*flat_scalar_prod by sprod
 *
 * Revision 1.2  2002/12/10 15:29:29  k_taniguchi
 * Change the multiplication "*" to "%"
 *   and flat_scalar_prod to flat_scalar_prod_desal.
 *
 * Revision 1.1.1.1  2001/11/20 15:19:28  e_gourgoulhon
 * LORENE
 *
 * Revision 2.11  2000/12/22  13:10:32  eric
 * Modif des annulations en dehors de l'etoile.
 *
 * Revision 2.10  2000/03/29  11:47:53  eric
 * Suppression affichage.
 *
 * Revision 2.9  2000/03/01  16:12:09  eric
 * Appel de set_std_base sur u_euler dans le cas irrotationnel.
 *
 * Revision 2.8  2000/02/24  09:34:10  keisuke
 * Ajout de l'appel a set_std_base() sur wit_w.
 *
 * Revision 2.7  2000/02/21  13:59:09  eric
 * Appel de set_dzpuis(0) sur loggam.
 *
 * Revision 2.6  2000/02/21  08:43:17  keisuke
 * Ajout de l'appel a set_std_base() sur loggam.
 *
 * Revision 2.5  2000/02/17  14:42:45  eric
 * Ajout de l'appel a set_std_base() sur gam_euler.
 *
 * Revision 2.4  2000/02/14  11:06:15  eric
 * Suppression de l'appel a annule(nzet.nzm1) sur gam_euler dans le cas en
 *   corotation.
 * Ajout de l'appel a annule(nzet,nzm1) sur wit_w.
 *
 * Revision 2.3  2000/02/12  18:36:23  eric
 * Appel de set_std_base sur loggam.
 *
 * Revision 2.2  2000/02/08  19:29:03  eric
 * Calcul sur les tenseurs.
 * wit_w est calcule.
 *
 * Revision 2.1  2000/02/04  16:37:28  eric
 * Premiere version implementee (non testee !)/
 *
 * Revision 2.0  2000/01/31  15:57:30  eric
 * *** empty log message ***
 *
 *
 * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_hydro.C,v 1.7 2016/12/05 16:17:53 j_novak Exp $
 *
 */

// Headers C

// Headers Lorene
#include "etoile.h"

namespace Lorene {
void Etoile_bin::hydro_euler(){

    int nz = mp.get_mg()->get_nzone() ; 
    int nzm1 = nz - 1 ; 

    //----------------------------------
    // Specific relativistic enthalpy           ---> hhh
    //----------------------------------
    
    Tenseur hhh = exp(unsurc2 * ent) ;  // = 1 at the Newtonian limit
    hhh.set_std_base() ;

    //---------------------------------------------------
    // Lorentz factor between the co-orbiting             ---> gam0
    // observer and the Eulerian one
    // See Eq (23) and (24) from Gourgoulhon et al. (2001)
    //---------------------------------------------------

    Tenseur gam0 = 1/sqrt(1-unsurc2*sprod(bsn,bsn)) ;
    gam0.set_std_base() ;

    //------------------------------------------
    // Lorentz factor and 3-velocity of the fluid 
    //  with respect to the Eulerian observer
    //------------------------------------------
    
    if (irrotational) { 

//##    cout << "Etoile_bin::hydro_euler : ### warning : " << endl ; 
//  cout << "   d_psi.d_psi would be better computed in spher. coord. !"
//##         << endl ; 

        d_psi.set_std_base() ;

    // See Eq (32) from Gourgoulhon et al. (2001)
    gam_euler = sqrt( 1 + unsurc2 * sprod(d_psi, d_psi)
                    / (hhh%hhh) ) ; 

    gam_euler.set_std_base() ;  // sets the standard spectral bases for
                    // a scalar field


    // See Eq (31) from Gourgoulhon et al. (2001)
    Tenseur vtmp = d_psi / ( hhh % gam_euler % a_car ) ; 

    // The assignment of u_euler is performed component by component 
    //  because u_euler is contravariant and d_psi is covariant
    if (vtmp.get_etat() == ETATZERO) {
        u_euler.set_etat_zero() ; 
    }
    else{
        assert(vtmp.get_etat() == ETATQCQ) ; 
        u_euler.set_etat_qcq() ; 
        for (int i=0; i<3; i++) {
        u_euler.set(i) = vtmp(i) ;
        }
        u_euler.set_triad( *(vtmp.get_triad()) ) ; 
    }
    
    u_euler.set_std_base() ; 

    }
    else {      // Rigid rotation
            // --------------

    gam_euler = gam0 ; 
    gam_euler.set_std_base() ;  // sets the standard spectral bases for
                    // a scalar field

    u_euler = -bsn ; 

    }
    
    //------------------------------------
    //  Energy density E with respect to the Eulerian observer
    // See Eq (53) from Gourgoulhon et al. (2001)  
    //------------------------------------
    
    ener_euler = gam_euler % gam_euler % ( ener + press ) - press ; 
    
    //------------------------------------
    // Trace of the stress tensor with respect to the Eulerian observer
    // See Eq (54) from Gourgoulhon et al. (2001)  
    //------------------------------------

    //Tenseur tmp00 = sprod(u_euler, u_euler) ; 
    //cout << hex << u_euler(0).va.base.b[0] << endl ; 
    //cout << hex << u_euler(1).va.base.b[0] << endl ; 
    //cout << hex << u_euler(2).va.base.b[0] << endl ; 
    //cout << hex << tmp00().va.base.b[0] << endl ; 

    s_euler = 3 * press  +  ( ener_euler + press ) *
      sprod(u_euler, u_euler) ;
    s_euler.set_std_base() ; 


    //-------------------------------------------
    //  Lorentz factor between the fluid and        ---> gam
    //  co-orbiting observers
    // See Eq (58) from Gourgoulhon et al. (2001)  
    //--------------------------------------------
    
    Tenseur tmp = ( 1+unsurc2*sprod(bsn,u_euler) ) ;
    tmp.set_std_base() ;
    Tenseur gam = gam0 % gam_euler % tmp ;

    //-------------------------------------------
    //  Spatial projection of the fluid 3-velocity
    //  with respect to the co-orbiting observer
    //--------------------------------------------

    wit_w = gam_euler / gam % u_euler + gam0 % bsn ; 

    wit_w.set_std_base() ;  // set the bases for spectral expansions

    wit_w.annule(nzm1) ;    // zero in the ZEC


    //-------------------------------------------
    //  Logarithm of the Lorentz factor between 
    //  the fluid and co-orbiting observers
    //--------------------------------------------

    if (relativistic) {

    loggam = log( gam ) ;
    
        loggam.set_std_base() ;   // set the bases for spectral expansions
    }
    else {
  
    loggam = 0.5 * flat_scalar_prod_desal(wit_w, wit_w) ;

        loggam.set_std_base() ;   // set the bases for spectral expansions

//### Forcage a zero des termes en sin(m*phi) : 
//  loggam.coef() ; 
//  int np = mgrille->np[0] ; 
//  int nt = mgrille->nt[0] ;
//  int nr = mgrille->nr[0] ;
//  int ntnr = nt * nr ; 
//  for (int k = 1; k < np+2; k+=2) {
//      for (int j = 0; j < nt; j++) {
//      for (int i = 0; i<nr; i++) {
//          loggam.c_cf[0]->t[0]->t[ntnr*k + nr*j + i] = 0 ;
//      }
//      }
//  }
//  loggam.c_ajx[0] = 0 ; 
//  loggam.c_aj = 0 ; 
//  loggam.coef_i() ; 
//###
    }


    //-------------------------------------------------
    // Velocity fields set to zero in external domains
    //-------------------------------------------------

    loggam.annule(nzm1) ;       // zero in the ZEC only

    wit_w.annule(nzm1) ;        // zero outside the star     

    u_euler.annule(nzm1) ;  // zero outside the star     


    if (loggam.get_etat() != ETATZERO) {
    (loggam.set()).set_dzpuis(0) ; 
    }
    
    //################
    // verification: test on gam_euler

    // if (irrotational) {
    
    // Tenseur gam_test = 1. / sqrt( 1 - unsurc2 * sprod(u_euler, u_euler) ) ;
    
    // cout << "Etoile_bin::hydro_euler: test of gam_euler : " << endl ; 
    // cout << "  maximum error : " << endl ;
    // cout << max(gam_test() - gam_euler()) << endl ; 
    //cout << "  relative error : " << endl ; 
    // cout << diffrel(gam_test(), gam_euler()) << endl ; 
    
    // }
        
    //##################


    //### Test

    if (!irrotational) {

    //  Tenseur diff = gam - 1 ; 
    //  cout << "Etoile_bin::hydro_euler: rigid motion: difference gam <-> 1 : " 
    //       << endl ; 
    //  cout << norme(diff()) / norme(gam()) << endl ; 
    //
    //  cout << "Etoile_bin::hydro_euler: rigid motion: difference wit_w <-> 0 : " 
    //       << endl ; 
    //  cout << "Component x : " << endl << norme(wit_w(0)) << endl ; 
    //  cout << "Component y : " << endl << norme(wit_w(1)) << endl ; 
    //  cout << "Component z : " << endl << norme(wit_w(2)) << endl ; 

//####
    gam = 1 ; 
    loggam = 0 ; 
    wit_w = 0 ; 
    }
    
    // The derived quantities are obsolete
    // -----------------------------------
    
    del_deriv() ;                
    
    
}
}