strot_diff_cfc_equil.C

/*
 *  Function Star_rot_diff_CFC::equilibrium
 *
 *    (see file star_rot_diff_cfc.h for documentation).
 *
 */

/*
 *   Copyright (c) 2025 Santiago Jaraba
 *
 *   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 version 2
 *   as published by the Free Software Foundation.
 *
 *   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
 *
 */

// C headers
#include <cmath>

// Lorene headers
#include "star_rot_diff_cfc.h"
#include "param.h"
#include "graphique.h"
#include "utilitaires.h"
#include "unites.h"


namespace Lorene{
void Star_rot_diff_CFC::equilibrium(double ent_c, double omega_c0, double fact_omega, 
                int , const Tbl& ent_limit, 
                const Itbl& icontrol,
                const Tbl& control, double mbar_wanted, 
                double aexp_mass, Tbl& diff, Param*){

    // Fundamental constants and units
    // --------------------------------
    using namespace Unites ;

    // For the display 
    // ---------------
    char display_bold[]="x[1m" ; display_bold[0] = 27 ;
    char display_normal[] = "x[0m" ; display_normal[0] = 27 ;


    // Grid parameters
    // ----------------

    const Mg3d* mg = mp.get_mg() ;
    int nz = mg->get_nzone() ;    // total number of domains
    int nzm1 = nz - 1 ;

    // Index of the point at phi=0, theta=pi/2 at the surface of the star:
    int type_t = mg->get_type_t() ; 
    assert( ( type_t == SYM) || (type_t == NONSYM) ) ; 
    int l_b = nzet - 1 ; 
    int i_b = mg->get_nr(l_b) - 1 ; 
    int j_b = (type_t == SYM ? mg->get_nt(l_b) - 1 : mg->get_nt(l_b)/2 ) ; 
    int k_b = 0 ;
    
    // Value of the enthalpy defining the surface of the star
    double ent_b = ent_limit(nzet-1) ;
    
    // Parameters to control the iteration
    // -----------------------------------
    
    int mer_max = icontrol(0) ; 
    int mer_rot = icontrol(1) ;
    int mer_change_omega = icontrol(2) ; 
    int mer_fix_omega = icontrol(3) ; 
    int mer_mass = icontrol(4) ; 
    int mer_triax = icontrol(5) ; 
    int delta_mer_kep = icontrol(6) ;  

    // Protections:
    if (mer_change_omega < mer_rot) {
    cout << "Star_rot_diff_CFC::equilibrium: mer_change_omega < mer_rot !" << endl ;
    cout << " mer_change_omega = " << mer_change_omega << endl ; 
    cout << " mer_rot = " << mer_rot << endl ; 
    abort() ; 
    }
    if (mer_fix_omega < mer_change_omega) {
    cout << "Star_rot_diff_CFC::equilibrium: mer_fix_omega < mer_change_omega !" 
         << endl ;
    cout << " mer_fix_omega = " << mer_fix_omega << endl ; 
    cout << " mer_change_omega = " << mer_change_omega << endl ; 
    abort() ; 
    }

    // In order to converge to a given baryon mass, shall the central
    // enthalpy be varied or Omega ?
    bool change_ent = true ; 
    if (mer_mass < 0) {
        change_ent = false ; 
        mer_mass = abs(mer_mass) ;
    }
    
    double precis = control(0) ; 
    double omega_ini = control(1) ; 
    double relax = control(2) ;
    double relax_prev = double(1) - relax ;  
    double ampli_triax = control(3) ; 


    // Error indicators
    // ----------------
    
    diff.annule_hard() ;
    double& diff_ent = diff.set(0) ; 
    double& vit_triax = diff.set(7) ; 

    double alpha_r = 1 ;  

    // Initializations
    // ---------------

    // Initial central angular velocity 
    double omega_c = 0 ;
    
    double accrois_omega = (omega_c0 - omega_ini) / 
                double(mer_fix_omega - mer_change_omega) ; 
 

    update_metric() ;    //update of the metric quantities 

    equation_of_state() ;  // update of the density, pressure,...etc

    hydro_euler() ; //update of the hydro quantities relative to the 
                    //  Eulerian observer

    // Quantities at the previous step :
    Scalar ent_prev = ent ;
    Scalar logn_prev = logn ;
    Scalar psi_prev = psi ;
    //  Vector beta_prev = beta ;

    // Output files
    // -------------

    ofstream fichconv("convergence.d") ;    // Output file for diff_ent
    fichconv << "#     diff_ent    GRV2    max_triax   vit_triax" << endl ;   

    ofstream fichfreq("frequency.d") ;    // Output file for omega_c
    fichfreq << "#       f [Hz]" << endl ; 

    ofstream fichevol("evolution.d") ;    // Output file for various quantities
    fichevol << 
        "#       |dH/dr_eq/dH/dr_pole|      r_pole/r_eq ent_c" 
        << endl ; 

    diff_ent = 1 ;
    double err_grv2 = 1 ;
    double max_triax_prev = 0 ;     // Triaxial amplitude at previous step


    //=========================================================================
    //          Start of iteration
    //=========================================================================

    for(int mer=0 ; (diff_ent > precis) && (mer<mer_max) ; mer++ ) {

    cout << "-----------------------------------------------" << endl ;
    cout << "step: " << mer << endl ;
    cout << "diff_ent = " << display_bold << diff_ent << display_normal
        << endl ;    
    cout << "err_grv2 = " << err_grv2 << endl ;    
    fichconv << mer ;
    fichfreq << mer ;
    fichevol << mer ;

    // Switch on rotation    
    if (mer >= mer_rot) {
    
        if (mer < mer_change_omega) {
        omega_c = omega_ini ; 
        }
        else {
        if (mer <= mer_fix_omega) {
            omega_c = omega_ini + accrois_omega * 
                    (mer - mer_change_omega) ;
        }
        }


    }

 //---------------------------------------------//
 // Resolution of the Poisson equation for psi  //
 //---------------------------------------------//

 Scalar psi_new(mp) ;

 solve_psi( psi_new ) ;

 psi_new.std_spectral_base() ;


 //---------------------------------------------------//
 // Resolution of the Poisson equation for logn       //
 // Note: ln_f is due to the fluid part               //
 //       ln_q is due to the quadratic metric part    //
 //---------------------------------------------------//

 Scalar ln_f_new(mp) ;
 Scalar ln_q_new(mp) ;

 solve_logn_f( ln_f_new ) ;
 solve_logn_q( ln_q_new ) ;

 ln_f_new.std_spectral_base() ;
 ln_q_new.std_spectral_base() ;

    //---------------------------------------
    // Triaxial perturbation of ln_ff
    //---------------------------------------
    
    if (mer == mer_triax) {
    
        if ( mg->get_np(0) == 1 ) {
        cout << 
        "Star_rot_diff::equilibrium: np must be stricly greater than 1"
        << endl << " to set a triaxial perturbation !" << endl ; 
        abort() ; 
        }
    
        const Coord& phi = mp.phi ; 
        const Coord& sint = mp.sint ; 
        Scalar perturb(mp) ; 
        perturb = 1 + ampli_triax * sint*sint * cos(2*phi) ;
        ln_f_new = ln_f_new * perturb ;  
        
        ln_f_new.std_spectral_base() ;    // set the bases for spectral expansions
                    //  to be the standard ones for a 
                    //  scalar field

    }
    
    // Monitoring of the triaxial perturbation
    // ---------------------------------------
    
    const Valeur& va_nuf = ln_f_new.get_spectral_va() ;
    va_nuf.coef() ;     // Computes the spectral coefficients
    double max_triax = 0 ; 

    if ( mg->get_np(0) > 1 ) {

        for (int l=0; l<nz; l++) {  // loop on the domains
        for (int j=0; j<mg->get_nt(l); j++) {
            for (int i=0; i<mg->get_nr(l); i++) {
        
            // Coefficient of cos(2 phi) : 
            double xcos2p = (*(va_nuf.c_cf))(l, 2, j, i) ; 

            // Coefficient of sin(2 phi) : 
            double xsin2p = (*(va_nuf.c_cf))(l, 3, j, i) ; 
            
            double xx = sqrt( xcos2p*xcos2p + xsin2p*xsin2p ) ; 
            
            max_triax = ( xx > max_triax ) ? xx : max_triax ;           
            }
        } 
        }
        
    }
    
    cout << "Triaxial part of nuf : " << max_triax << endl ;
    

 //--------------------------------------------------//
 // Resolution of the Poisson equation for shift     //
 //--------------------------------------------------//

 Vector beta_new(mp, CON, mp.get_bvect_spher()) ;

 solve_shift( beta_new ) ;

    //-----------------------------------------
    // Determination of the fluid velociy U
    //-----------------------------------------
    
    if (mer > mer_fix_omega + delta_mer_kep) {
        
            omega_c *= fact_omega ;  // Increase of the angular velocity if 
    }              //  fact_omega != 1
    
    
    bool omega_trop_grand = false ; 
    bool kepler = true ; 

    while ( kepler ) {
    
        // Possible decrease of Omega to ensure a velocity < c 
    
        bool superlum = true ; 
    
        while ( superlum ) {
    
        // Computation of Omega(r,theta)

        if (omega_c == 0.) {
            omega_field = 0 ; 
        }
        else if (par_frot.get_taille() == 5) { // Uryu rotation profile
            if (par_frot(1) == 1. && par_frot(2) == 1.) { // Rigid rotation limit in Uryu profile
                omega_field = omega_c ;
                for (int l=nzet+1; l<nz; l++) { // omega_field set to 0 far from the star
                    omega_field.set_domain(l) = 0 ;
                }
                prim_field = 0 ;
                omega_min = omega_c ;
                omega_max = omega_c ;
            }
            else { // Usual case
            double lambda1 = par_frot(1) ;
            double lambda2 = par_frot(2) ;
            int p = 1 ;
            int q = 3 ;

            // First, compute F_e and F_max
            Scalar p4rst2 = psi4 ;
            p4rst2.annule_domain(nzm1) ; 
            p4rst2.std_spectral_base() ;
            p4rst2.mult_rsint() ; p4rst2.mult_rsint() ;     //  Multiplication by r^2 sin^2(theta)       
            Scalar nnn2 = nn * nn ;

            double omnp = omega_field.val_grid_point(l_b, k_b, j_b, i_b) + beta(3).val_grid_point(l_b, k_b, j_b, i_b) ;
            double F_e = p4rst2.val_grid_point(l_b, k_b, j_b, i_b) * omnp /
            ( nnn2.val_grid_point(l_b, k_b, j_b, i_b) - 
            p4rst2.val_grid_point(l_b, k_b, j_b, i_b) * omnp * omnp ) ;

            if (F_e != double(0)) { // Not the first iteration
                int l, k, j, i ;
                double om_max = 0 ;
                double om_lkji = 0 ;
                double F_max = 0 ;
                for (l=0; l<nzet+1; 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++) {
                        om_lkji = omega_field.val_grid_point(l, k, j, i) ;
                        if (om_lkji > om_max) {
                            om_max = om_lkji ;
                            omnp = om_max + beta(3).val_grid_point(l, k, j, i) ;
                            F_max = p4rst2.val_grid_point(l, k, j, i) * omnp /
                            ( nnn2.val_grid_point(l, k, j, i) - 
                            p4rst2.val_grid_point(l, k, j, i) * omnp * omnp ) ;
                        }
                    }
                    }
                }
                }

                // Now A^2 and B^2, parameters for omega_frot, can be computed

                if (lambda2 == -1.) { // Special flag to keep A^2 and B^2 constant
                    // Computing lambda1 from A^2 and B^2
                    double A2 = par_frot(3) ;
                    double B2 = par_frot(4) ;
                    lambda1 = (1 + pow(F_max/(B2*omega_c), p))/(1 + pow(F_max/(A2*omega_c), q+p)) ;                 
                }
                else if (F_e < 0) { // Work in progress: fixing issues with A^2 and B^2 parameters
                    cout << "Warning: F_e=" << F_e << " < 0. Aborting." << endl ;
                    abort() ;
                    //F_e = 1.001*pow(lambda1/lambda2, 1./q)*F_max ;
                    //cout << F_e << endl ;
                    par_frot.set(3) = 10000. ;
                    par_frot.set(4) = 10000. ;
                }
                /*else if (F_e < F_max) {
                    cout << "Warning: F_e=" << F_e << " < F_max=" << F_max << ". Impossible to fix by changing lambda1 and lambda2." << endl ;
                    cout << "Trying with the parameters of the previous iteration" << endl ;
                }*/
                else {
                /*while (F_e < pow(lambda1/lambda2, 1./q)*F_max) {
                    // This would lead to unphysical A^2 and B^2. Relaxing lambda1 and lambda2 only in this iteration
                    lambda1 *= .9 ;
                    lambda2 *= 1.1 ;
                    // lambda2 = 1.1*lambda1*pow(F_max/F_e, q) ;
                    cout << "Warning: unphysical A^2 and B^2. Relaxing lambda1=" << par_frot(1)
                         << "->" << lambda1 << ", lambda2=" << par_frot(2) << "->" << lambda2 << endl ;

                }*/
                //cout << lambda2*pow(F_e, q)-lambda1*pow(F_max, q) << endl ;
                par_frot.set(3) = pow(pow(F_e*F_max,p)/pow(omega_c, p+q) * (lambda2*pow(F_e, q)-lambda1*pow(F_max, q) )
                                / ((lambda1-1)*pow(F_e, p) - (lambda2-1)*pow(F_max, p)), 1./(p+q) ) ;
                par_frot.set(4) = F_e*F_max/omega_c * pow((lambda2*pow(F_e, q)-lambda1*pow(F_max, q))
                                / (lambda2*(lambda1-1)*pow(F_e, p+q) - lambda1*(lambda2-1)*pow(F_max, p+q)), 1./p) ;

                cout << F_e << " " << F_max << endl ;
                }
                
            }

            par_frot.set(0) = omega_c ;

            cout << "Parameters of Omega(F) : " << par_frot << endl ; 
            
            double omeg_min = 0 ;
            double omeg_max = lambda1 * omega_c ;
            double precis1 = 1.e-14 ;
            int nitermax1 = 200 ; 

            fait_omega_field(omeg_min, omeg_max, precis1, nitermax1) ;
            }
        }
        else { // Usual rotation profile
            par_frot.set(0) = omega_c ; 
            if (par_frot(2) != double(0)) {  // fixed a = R_eq / R_0
            par_frot.set(1) = ray_eq() / par_frot(2) ; 
            }
            double omeg_min = 0 ; 
            double omeg_max = omega_c ;
            double precis1 = 1.e-14 ;
            int nitermax1 = 100 ; 

            fait_omega_field(omeg_min, omeg_max, precis1, nitermax1) ;
        }

        // New fluid velocity :
        //
        
        u_euler.set(1).set_etat_zero() ;
        u_euler.set(2).set_etat_zero() ;

        u_euler.set(3) = omega_field ;
        u_euler.set(3).annule(nzet,nzm1) ;   // nzet is defined in class Star
        u_euler.set(3).std_spectral_base() ;
        u_euler.set(3).mult_rsint() ;
        u_euler.set(3) += beta(3) ;
        u_euler.set(3).annule(nzet,nzm1) ;
        
        u_euler = u_euler / nn ; 


        // v2 (square of norm of u_euler)
        // -------------------------------

        v2 = contract(contract(gamma.cov(), 0, u_euler, 0), 0, u_euler, 0) ;

        // Is the new velocity larger than c in the equatorial plane ?
    
        superlum = false ; 
    
        for (int l=0; l<nzet; l++) {
            for (int i=0; i<mg->get_nr(l); i++) {
        
            double u2 = v2.val_grid_point(l, 0, j_b, i) ; 
            if (u2 >= 1.) {     // superluminal velocity
                superlum = true ; 
                cout << "U > c  for l, i : " << l << "  " << i 
                 << "   U = " << sqrt(u2) << endl ;  
            }
            }
        }
        if ( superlum ) {
            cout << "**** VELOCITY OF LIGHT REACHED ****" << endl ; 
            omega_c /= fact_omega ;    // Decrease of Omega_c
            cout << "New central rotation frequency : " 
             << omega/(2.*M_PI) * f_unit <<  " Hz" << endl ; 
            omega_trop_grand = true ;  
        }
        }   // end of while ( superlum )

    
        // New computation of U (which this time is not superluminal)
        //  as well as of gam_euler, ener_euler, etc...
        // -----------------------------------

        hydro_euler() ; 
    

        //------------------------------------------------------
        //  First integral of motion 
        //------------------------------------------------------
    
        Scalar mlngamma(mp) ;  //  -log( gam_euler )

        mlngamma = - log( gam_euler ) ;
    
        // Equatorial values of various potentials :
        double ln_f_b = ln_f_new.val_grid_point(l_b, k_b, j_b, i_b) ; 
        double ln_q_b = ln_q_new.val_grid_point(l_b, k_b, j_b, i_b) ; 
        double mlngamma_b = mlngamma.val_grid_point(l_b, k_b, j_b, i_b) ; 
        double primf_b = prim_field.val_grid_point(l_b, k_b, j_b, i_b) ;
        
        
        // Central values of various potentials :
        double ln_f_c = ln_f_new.val_grid_point(0,0,0,0) ; 
        double ln_q_c = ln_q_new.val_grid_point(0,0,0,0) ; 
        double mlngamma_c = 0 ;
        double primf_c = prim_field.val_grid_point(0,0,0,0) ; 
    
        // Scale factor to ensure that the (log of) enthalpy is equal to 
        // ent_b at the equator
        double alpha_r2 = ( ent_c - ent_b + mlngamma_c - mlngamma_b
                + ln_q_c - ln_q_b + primf_c - primf_b) 
                / ( ln_f_b - ln_f_c  ) ;
        alpha_r = sqrt(alpha_r2) ;

        cout << "alpha_r = " << alpha_r << endl ; 

        // Rescaling of the grid (no adaptation!)
        //---------------------------------------
        mp.homothetie(alpha_r) ;
        
        // Readjustment of logn :
        // -----------------------

        logn = alpha_r2 * ln_f_new + ln_q_new ;

        double logn_c = logn.val_grid_point(0,0,0,0) ;

        // First integral of motion -> (log of) enthalpy in all space
        // ----------------------------------------------------------

        ent = (ent_c + logn_c + mlngamma_c) - logn - mlngamma - prim_field ;

        // Test: is the enthalpy negative somewhere in the equatorial plane
        //  inside the star?
        // --------------------------------------------------------

        kepler = false ; 
        for (int l=0; l<nzet; l++) {
        int imax = mg->get_nr(l) - 1 ;
        if (l == l_b) imax-- ;  // The surface point is skipped
        for (int i=0; i<imax; i++) { 
            if ( ent.val_grid_point(l, 0, j_b, i) < 0. ) {
            kepler = true ;
            cout << "ent < 0 for l, i : " << l << "  " << i 
                 << "   ent = " << ent.val_grid_point(l, 0, j_b, i) << endl ;  
            } 
        }
        }
    
        if ( kepler ) {
        cout << "**** KEPLERIAN VELOCITY REACHED ****" << endl ; 
        omega_c /= fact_omega ;    // Omega is decreased
        cout << "New central rotation frequency : " 
             << omega_c/(2.*M_PI) * f_unit << " Hz" << endl ; 
        omega_trop_grand = true ;  
        }

    }   // End of while ( kepler )
    
    if ( omega_trop_grand ) {   // fact_omega is decreased for the
                                //  next step 
        fact_omega = sqrt( fact_omega ) ; 
        cout << "**** New fact_omega : " << fact_omega << endl ; 
    }


    //----------------------------------------------------
    // Equation of state  
    //----------------------------------------------------
    
    equation_of_state() ;  // computes new values for nbar (n), ener (e) 
                           // and press (p) from the new ent (H)

    hydro_euler() ;

    beta = beta_new ;

    //---------------------------------------
    // Calculate error of the GRV2 identity 
    //---------------------------------------

    err_grv2 = grv2() ;


    //--------------------------------------
    // Relaxations on some quantities....?
    //
    //--------------------------------------

    if (mer >= 10) {
    logn = relax * logn + relax_prev * logn_prev ;

    psi = relax * psi_new + relax_prev * psi_prev ;

    }
    
    // Update of the metric quantities :

    update_metric() ; 
    
    //-----------------------
    //  Informations display
    // More to come later......
    //-----------------------

    partial_display(cout) ; 
    fichfreq << "  " << omega_c / (2*M_PI) * f_unit ; 
    fichevol << "  " << ray_pole() / ray_eq() ; 
    fichevol << "  " << ent_c ; 

    //-----------------------------------------
    // Convergence towards a given baryon mass 
    //-----------------------------------------

    if (mer > mer_mass) {
        
        double xx ; 
        if (mbar_wanted > 0.) {
        xx = mass_b() / mbar_wanted - 1. ;
        cout << "Discrep. baryon mass <-> wanted bar. mass : " << xx 
             << endl ; 
        }
        else{
        xx = mass_g() / fabs(mbar_wanted) - 1. ;
        cout << "Discrep. grav. mass <-> wanted grav. mass : " << xx 
             << endl ; 
        }
        double xprog = ( mer > 2*mer_mass) ? 1. : 
                 double(mer-mer_mass)/double(mer_mass) ; 
        xx *= xprog ; 
        double ax = .5 * ( 2. + xx ) / (1. + xx ) ; 
        double fact = pow(ax, aexp_mass) ; 
        cout << "  xprog, xx, ax, fact : " << xprog << "  " <<
            xx << "  " << ax << "  " << fact << endl ; 

        if ( change_ent ) {
        ent_c *= fact ; 
        }
        else {
        if (mer%4 == 0) omega_c *= fact ; 
        }
    }
    

    //------------------------------------------------------------
    //  Relative change in enthalpy with respect to previous step 
    // ** Check:  Is diffrel(ent, ent_prev) ok?  
    //------------------------------------------------------------

    Tbl diff_ent_tbl = diffrel( ent, ent_prev ) ; 
    diff_ent = diff_ent_tbl(0) ; 
    for (int l=1; l<nzet; l++) {
        diff_ent += diff_ent_tbl(l) ; 
    }
    diff_ent /= nzet ; 
    
    fichconv << " " << log10( fabs(diff_ent) + 1.e-16 ) ;
    fichconv << " " << log10( fabs(err_grv2) + 1.e-16 ) ;
    fichconv << "  " << log10( fabs(max_triax) + 1.e-16 ) ;

    vit_triax = 0 ;
    if ( (mer > mer_triax+1) && (max_triax_prev > 1e-13) ) {
        vit_triax = (max_triax - max_triax_prev) / max_triax_prev ; 
    }

    fichconv << "  " << vit_triax ;
    
    //------------------------------
    //  Recycling for the next step
    //------------------------------
    
    ent_prev = ent ; 
    logn_prev = logn ; 
    psi_prev = psi ;
    max_triax_prev = max_triax ; 
    
    fichconv << endl ;
    fichfreq << endl ;
    fichevol << endl ;
    fichconv.flush() ; 
    fichfreq.flush() ; 
    fichevol.flush() ; 


    } // End of main loop
    
    //=========================================================================
    //          End of iteration
    //=========================================================================

    fichconv.close() ; 
    fichfreq.close() ; 
    fichevol.close() ; 

}
}