time_slice.h

/*
 *  Definition of Lorene class Time_slice
 *
 */

/*
 *   Copyright (c) 2004 Eric Gourgoulhon, Jose Luis Jaramillo & Jerome Novak
 *
 *   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
 *
 */

#ifndef __TIME_SLICE_H_ 
#define __TIME_SLICE_H_ 

/*
 * $Id: time_slice.h,v 1.33 2019/12/06 13:56:54 j_novak Exp $
 * $Log: time_slice.h,v $
 * Revision 1.33  2019/12/06 13:56:54  j_novak
 * Corrected mistake in a comment
 *
 * Revision 1.32  2014/10/13 08:52:37  j_novak
 * Lorene classes and functions now belong to the namespace Lorene.
 *
 * Revision 1.31  2012/02/06 12:59:07  j_novak
 * Correction of some errors.
 *
 * Revision 1.30  2010/10/20 07:58:09  j_novak
 * Better implementation of the explicit time-integration. Not fully-tested yet.
 *
 * Revision 1.29  2008/12/04 18:22:49  j_novak
 * Enhancement of the dzpuis treatment + various bug fixes.
 *
 * Revision 1.28  2008/12/02 15:02:21  j_novak
 * Implementation of the new constrained formalism, following Cordero et al. 2009
 * paper. The evolution eqs. are solved as a first-order system. Not tested yet!
 *
 * Revision 1.27  2008/08/19 06:41:59  j_novak
 * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
 * cast-type operations, and constant strings that must be defined as const char*
 *
 * Revision 1.26  2007/11/06 14:47:06  j_novak
 * New constructor from a rotating star in Dirac gauge (class Star_rot_Dirac).
 * Evolution can take into account matter terms.
 *
 * Revision 1.25  2007/04/25 15:20:59  j_novak
 * Corrected an error in the initialization of tildeB in
 * Tslice_dirac_max::initial_dat_cts. + New method for solve_hij_AB.
 *
 * Revision 1.24  2007/03/21 14:51:48  j_novak
 * Introduction of potentials A and tilde(B) of h^{ij} into Tslice_dirac_max.
 *
 * Revision 1.23  2005/03/28 19:44:00  f_limousin
 * Function tgam() is now virtual.
 *
 * Revision 1.22  2004/06/24 20:36:07  e_gourgoulhon
 * Class Time_slice_conf: added method check_psi_dot.
 *
 * Revision 1.21  2004/06/14 20:46:35  e_gourgoulhon
 * Added argument method_poisson to Tslice_dirac_max::solve_hij.
 *
 * Revision 1.20  2004/05/31 20:28:20  e_gourgoulhon
 * -- Class Time_slice : added inline functions get_latest_j() and
 *    get_time()
 * -- Class Tslice_dirac_max: method hh_det_one takes now a time step
 *    as argument, to compute h^{ij} from khi and mu at an arbitrary
 *    time step and not only the latest one.
 *
 * Revision 1.19  2004/05/27 15:22:28  e_gourgoulhon
 * Added functions save and sauve, as well as constructors from file.
 *
 * Revision 1.18  2004/05/20 20:30:37  e_gourgoulhon
 * Added arguments check_mod and save_mod to method Tsclice_dirac_max::evolve.
 *
 * Revision 1.17  2004/05/17 19:52:16  e_gourgoulhon
 * -- Method initial_data_cts: added arguments graph_device and
 *    method_poisson_vect.
 * -- Method Tslice_dirac_max::solve_beta : added argument method
 * -- Method Tslice_dirac_max::solve_hij : added argument graph_device
 * -- Method Tslice_dirac_max::evolve : added arguments
 *    method_poisson_vect, nopause and graph_device.
 *
 * Revision 1.16  2004/05/12 15:16:25  e_gourgoulhon
 * Added #include "metric.h" before #include "evolution.h".
 *
 * Revision 1.15  2004/05/09 20:56:09  e_gourgoulhon
 * Added member adm_mass_evol and corresponding virtual method adm_mass().
 *
 * Revision 1.14  2004/05/06 15:23:10  e_gourgoulhon
 * initial_data_cts is know a virtual function of class Time_slice_conf
 * and is implemented also for class Tslice_dirac_max.
 *
 * Revision 1.13  2004/05/03 14:46:11  e_gourgoulhon
 * Class Tslice_dirac_max: -- changed prototype of method solve_hij
 *                         -- added new method evolve
 *
 * Revision 1.12  2004/04/30 14:36:15  j_novak
 * Added the method Tslice_dirac_max::solve_hij(...)
 * NOT READY YET!!!
 *
 * Revision 1.11  2004/04/30 10:51:38  e_gourgoulhon
 * Class Tslice_dirac_max: added methods solve_n, solve_q and solve_beta
 * for resolution of the elliptic part of Einstein equations.
 *
 * Revision 1.10  2004/04/29 17:07:27  e_gourgoulhon
 * Added argument pdt to Time_slice_conf::initial_data_cts.
 *
 * Revision 1.9  2004/04/08 16:42:11  e_gourgoulhon
 * Many changes:
 * -- class Time_slice_conf: added methods set_*, changed argument list of
 *    method initial_data_cts.
 * -- class Tslice_dirac_max: added methods set_* and  hh_det_one().
 *
 * Revision 1.8  2004/04/05 21:21:51  e_gourgoulhon
 * class Time_slice_conf: added method initial_data_cts (computation of
 *  initial data from conformally thin sandwich method).
 * classes  Time_slice_conf and Tslice_dirac_max: added constructor as
 *  standard time slice of Minkowski spacetime.
 *
 * Revision 1.7  2004/04/01 16:09:01  j_novak
 * Trace of K_ij is now member of Time_slice (it was member of Time_slice_conf).
 * Added new methods for checking 3+1 Einstein equations (preliminary).
 *
 * Revision 1.6  2004/03/30 14:00:30  j_novak
 * New class Tslide_dirac_max (first version).
 *
 * Revision 1.5  2004/03/29 11:58:53  e_gourgoulhon
 * Many modif. to class Time_slice_conf.
 * Minor modif. to class Time_slice.
 *
 * Revision 1.4  2004/03/28 21:33:14  e_gourgoulhon
 * Constructor  Time_slice::Time_slice(int depth_in) declared "explicit".
 *
 * Revision 1.3  2004/03/28 21:27:57  e_gourgoulhon
 * Class Time_slice: - renamed the Evolution_std with suffix "_evol".
 *                   - added protected constructor for derived classes
 * Added class Time_slice_conf.
 *
 * Revision 1.2  2004/03/26 13:33:02  j_novak
 * New methods for accessing/updating members (nn(), beta(), gam_uu(), k_uu(), ...)
 *
 * Revision 1.1  2004/03/24 14:56:18  e_gourgoulhon
 * First version
 *
 *
 * $Header: /cvsroot/Lorene/C++/Include/time_slice.h,v 1.33 2019/12/06 13:56:54 j_novak Exp $
 *
 */

#include "star_rot_dirac.h"
#include "evolution.h"

                    //---------------------------//
                    //      class Time_slice     //
                    //---------------------------//

namespace Lorene {
class Time_slice {

    // Data : 
    // -----
    protected:
        int depth ; 
        
        int scheme_order ; 
        
        int jtime ; 
        
    Evolution_std<double> the_time ;
        
    mutable Evolution_std<Sym_tensor> gam_dd_evol ; 
        
    mutable Evolution_std<Sym_tensor> gam_uu_evol ; 

    mutable Evolution_std<Sym_tensor> k_dd_evol ; 

    mutable Evolution_std<Sym_tensor> k_uu_evol ; 

    mutable Evolution_std<Scalar> n_evol ; 
        
    mutable Evolution_std<Vector> beta_evol ; 
        
    mutable Evolution_std<Scalar> trk_evol ; 

        mutable Evolution_full<Tbl> adm_mass_evol ; 

    // Derived data : 
    // ------------
    protected:
    mutable Metric* p_gamma ;   

    // Constructors - Destructor
    // -------------------------
    public:
    
    Time_slice(const Scalar& lapse_in, const Vector& shift_in,
           const Sym_tensor& gamma_in, const Sym_tensor& kk_in,
           int depth_in = 3) ; 
    
    Time_slice(const Scalar& lapse_in, const Vector& shift_in,
           const Evolution_std<Sym_tensor>& gamma_in) ; 
    
    Time_slice(const Map& mp, const Base_vect& triad, int depth_in = 3) ; 
    
    Time_slice(const Map& mp, const Base_vect& triad, FILE* fich, 
           bool partial_read, int depth_in = 3) ; 
    
    Time_slice(const Time_slice& ) ;   
    
    protected:
    explicit Time_slice(int depth_in) ; 
    
    public:
    virtual ~Time_slice() ;         
 

    // Memory management
    // -----------------
    protected:
        
    virtual void del_deriv() const ; 
    
    void set_der_0x0() const ; 


    // Mutators / assignment
    // ---------------------
    public:
    void operator=(const Time_slice&) ;

    void set_scheme_order(int ord) { 
        assert ((0<= ord)&&(ord < 4)) ;
      scheme_order = ord ; } ; 
    
    // Accessors
    // ---------
    public:

    int get_scheme_order() const { return scheme_order ; } ;
        
        int get_latest_j() const {return jtime; } ;
        
        const Evolution_std<double>& get_time() const {return the_time; } ; 
    
    virtual const Scalar& nn() const ;
    
    virtual const Vector& beta() const ;
    
    const Metric& gam() const ;
    
    virtual const Sym_tensor& gam_dd() const ;
    
    virtual const Sym_tensor& gam_uu() const ;
    
    virtual const Sym_tensor& k_dd() const ;
    
    virtual const Sym_tensor& k_uu() const ;
    
        virtual const Scalar& trk() const ; 
        
    
    // Computational functions
    // -----------------------
    public:
     Tbl check_hamiltonian_constraint(const Scalar* energy_density = 0x0,
                      ostream& ost = cout, bool verb=true) const ;
    
     Tbl check_momentum_constraint(const Vector* momentum_density = 0x0,
                       ostream& ost = cout, bool verb=true) const ;
    
     Tbl check_dynamical_equations(const Sym_tensor* strain_tensor = 0x0,
                       const Scalar* energy_density = 0x0,
                       ostream& ost = cout, bool verb=true) const  ; 
    
        virtual double adm_mass() const ; 
        
    // Outputs
    // -------
    protected:
    virtual ostream& operator>>(ostream& ) const ; 
    
    friend ostream& operator<<(ostream& , const Time_slice& ) ; 

    public:
    void save(const char* rootname) const ; 
    
    protected:
    virtual void sauve(FILE* fich, bool partial_save) const ; 

};

ostream& operator<<(ostream& , const Time_slice& ) ;    



                    //---------------------------//
                    //   class Time_slice_conf   //
                    //---------------------------//

class Time_slice_conf : public Time_slice {

    // Data : 
    // -----
    protected: 
    
        const Metric_flat& ff ;  

    mutable Evolution_std<Scalar> psi_evol ; 
        
    mutable Evolution_std<Scalar> npsi_evol ; 
        
        
    mutable Evolution_std<Sym_tensor> hh_evol ; 

    mutable Evolution_std<Sym_tensor> hata_evol ; 

    mutable Evolution_std<Scalar> A_hata_evol ;
        
    mutable Evolution_std<Scalar> B_hata_evol ;
        
    // Derived data : 
    // ------------
    protected:
        mutable Metric* p_tgamma ; 
        
        mutable Scalar* p_psi4 ; 
        
        mutable Scalar* p_ln_psi ; 
        
        mutable Vector* p_hdirac ; 

    mutable Vector* p_vec_X ;

    // Constructors - Destructor
    // -------------------------
    public:
    
    Time_slice_conf(const Scalar& lapse_in, const Vector& shift_in,
            const Metric_flat& ff_in, const Scalar& psi_in, 
            const Sym_tensor& hh_in, const Sym_tensor& hata_in, 
            const Scalar& trk_in, int depth_in = 3) ; 
    
    
    Time_slice_conf(const Scalar& lapse_in, const Vector& shift_in,
               const Sym_tensor& gamma_in, const Sym_tensor& kk_in,
               const Metric_flat& ff_in, int depth_in = 3) ; 
               
    Time_slice_conf(const Map& mp, const Base_vect& triad, 
                    const Metric_flat& ff_in, int depth_in = 3) ; 
    
    Time_slice_conf(const Map& mp, const Base_vect& triad, 
                    const Metric_flat& ff_in, FILE* fich, 
                    bool partial_read, int depth_in = 3) ; 

    Time_slice_conf(const Time_slice_conf& ) ;  

    virtual ~Time_slice_conf() ;            
 

    // Memory management
    // -----------------
    protected:
        
    virtual void del_deriv() const ; 
    
    void set_der_0x0() const ; 


    // Mutators / assignment
    // ---------------------
    public:
    void operator=(const Time_slice_conf&) ;

    void operator=(const Time_slice&) ;
        
        virtual void set_psi_del_npsi(const Scalar& psi_in) ; 
        
        virtual void set_psi_del_n(const Scalar& psi_in) ; 
        
        virtual void set_npsi_del_psi(const Scalar& npsi_in) ; 
        
        virtual void set_npsi_del_n(const Scalar& npsi_in) ; 
        
        virtual void set_hh(const Sym_tensor& hh_in) ; 

        virtual void set_hata(const Sym_tensor& hata_in) ; 

    virtual void set_hata_TT(const Sym_tensor_tt& hata_tt) ;

    virtual void set_hata_from_XAB(Param* par_bc=0x0, Param* par_mat=0x0) ;

    // Accessors
    // ---------
    public:

        // Virtual functions from base class Time_slice:
        // ---------------------------------------------

    virtual const Scalar& nn() const ;
    
    virtual const Sym_tensor& gam_dd() const ;
    
    virtual const Sym_tensor& gam_uu() const ;
    
    virtual const Sym_tensor& k_dd() const ;
    
    virtual const Sym_tensor& k_uu() const ;
    
        // Virtual functions from this class:
        // ----------------------------------

    virtual const Scalar& A_hata() const ; 
    
    virtual const Scalar& B_hata() const ;
    
        virtual const Scalar& psi() const ; 
        
        const Scalar& psi4() const ; 
        
        const Scalar& ln_psi() const ; 
        
        virtual const Scalar& npsi() const ; 
        
        virtual const Metric& tgam() const ; 

        virtual const Sym_tensor& hh(Param* = 0x0, Param* = 0x0) const ; 

        virtual const Sym_tensor& hata() const ; 

        virtual Sym_tensor aa() const ; 

        virtual const Scalar& trk() const ; 
        
        virtual const Vector& hdirac() const ;

        virtual const Vector& vec_X(int method_poisson=6) const ;
        
    // Computational methods
    // ---------------------
    public:
    void compute_X_from_momentum_constraint
        (const Vector& hat_S, const Sym_tensor_tt& hata_tt, 
          int iter_max = 200, double precis = 1.e-12,
         double relax = 0.8, int methode_poisson = 6) ;

    virtual void set_AB_hata(const Scalar& A_in, const Scalar& B_in) ; 

         virtual void initial_data_cts(const Sym_tensor& uu, const Scalar& trk_in, 
                const Scalar& trk_point, double pdt, double precis = 1.e-12,
                int method_poisson_vect = 6, const char* graph_device = 0x0, 
                const Scalar* ener_dens = 0x0, const Vector* mom_dens = 0x0, 
                const Scalar* trace_stress = 0x0 ) ; 
        
        virtual double adm_mass() const ; 
        
        void check_psi_dot(Tbl& tlnpsi_dot, Tbl& tdiff, Tbl& tdiff_rel) const ; 
        
    // Outputs
    // -------
    protected:
    virtual ostream& operator>>(ostream& ) const ; 
    
    virtual void sauve(FILE* fich, bool partial_save) const ; 

} ; 
                    //----------------------------//
                    //   class Tslice_dirac_max   //
                    //----------------------------//

class Tslice_dirac_max : public Time_slice_conf {

  // Data : 
  // -----
 protected: 
  mutable Evolution_std<Scalar> A_hh_evol ;
  
  mutable Evolution_std<Scalar> B_hh_evol ;

  mutable Evolution_std<Scalar> source_A_hh_evol ;
  
  mutable Evolution_std<Scalar> source_B_hh_evol ;

  mutable Evolution_std<Scalar> source_A_hata_evol ;
  
  mutable Evolution_std<Scalar> source_B_hata_evol ;

  mutable Evolution_std<Scalar> trh_evol ;


    // Constructors - Destructor
    // -------------------------
    public:
    Tslice_dirac_max(const Scalar& lapse_in, const Vector& shift_in,
            const Metric_flat& ff_in, const Scalar& psi_in, 
            const Sym_tensor_trans& hh_in, const Sym_tensor& hata_in, 
            int depth_in = 3) ; 
    
    Tslice_dirac_max(const Map& mp, const Base_vect& triad, 
                     const Metric_flat& ff_in, int depth_in = 3) ; 
    
    Tslice_dirac_max(const Map& mp, const Base_vect& triad, 
                     const Metric_flat& ff_in, FILE* fich, 
                     bool partial_read, int depth_in = 3) ; 

    Tslice_dirac_max(const Star_rot_Dirac& star, double pdt, int depth_in = 3) ;

    Tslice_dirac_max(const Tslice_dirac_max& ) ;   

    virtual ~Tslice_dirac_max() ;           
 

    // Mutators / assignment
    // ---------------------
    public:
    void operator=(const Tslice_dirac_max&) ;   
    
        // Virtual functions from base class Time_slice_conf:
        // -------------------------------------------------

        virtual void set_hh(const Sym_tensor& hh_in) ; 

         virtual void initial_data_cts(const Sym_tensor& uu, const Scalar& trk_in, 
                const Scalar& trk_point, double pdt, double precis = 1.e-12,
                int method_poisson_vect = 6, const char* graph_device = 0x0, 
                const Scalar* ener_dens = 0x0, const Vector* mom_dens = 0x0, 
                const Scalar* trace_stress = 0x0 ) ; 
        
        // Virtual functions from this class:
        // ----------------------------------

    virtual void set_khi_mu(const Scalar& khi_in, const Scalar& mu_in) ; 

    virtual void set_AB_hh(const Scalar& A_in, const Scalar& B_in) ; 

    virtual void set_trh(const Scalar& trh_in) ;
    
    virtual Scalar solve_psi(const Scalar* ener_dens=0x0) const ; 
        
    virtual Scalar solve_npsi(const Scalar* ener_dens=0x0,
                  const Scalar* trace_stress=0x0) const ; 
        
    virtual Vector solve_beta(int method = 6) const ; 
        
    void evolve(double pdt, int nb_time_steps, int niter_elliptic,
                double relax_elliptic, int check_mod, int save_mod,
                int method_poisson_vect = 6, int nopause = 1, 
                const char* graph_device = 0x0, bool verbose=true, 
        const Scalar* ener_euler = 0x0,
        const Vector* mom_euler = 0x0, const Scalar* s_euler = 0x0,
        const Sym_tensor* strain_euler = 0x0) ; 
        
    virtual double adm_mass() const ; 
        
    protected:
    void compute_sources(const Sym_tensor* strain_tensor = 0x0) const ;

    void initialize_sources_copy() const ;

    void hh_det_one(int j, Param* par_bc = 0x0, Param* par_mat = 0x0) const ; 
    
    void hh_det_one(const Sym_tensor_tt& hijtt, Param* par_mat = 0x0) const ;

    // Accessors
    // ---------
    public:
    // Virtual functions from base class Time_slice_conf:
    // -------------------------------------------------
    
    virtual const Sym_tensor& hh(Param* par_bc = 0x0, Param* par_mat = 0x0) const ; 
    
    virtual const Scalar& trk() const ; 
    
    virtual const Vector& hdirac() const ; 
    
    // Virtual functions from this class:
    // ----------------------------------
    
    virtual const Scalar& A_hh() const ; 
    
    virtual const Scalar& B_hh() const ;
    
    virtual const Scalar& trh() const ;
    
    
    // Outputs
    // -------
    protected:
    virtual ostream& operator>>(ostream& ) const ;  
    
    virtual void sauve(FILE* fich, bool partial_save) const ; 
  
};
}
#endif