make_rot.cpp

#include "mpi.h"
#include "kadath_spheric.hpp"

using namespace Kadath ;

int main(int argc, char** argv) {
		
	int rc = MPI_Init (&argc, &argv) ;
	int rank = 0 ;
	MPI_Comm_rank (MPI_COMM_WORLD, &rank) ;
	
	double aa, ome, nhor, vhor ;
	FILE* ff = fopen(argv[1], "r") ;
	Space_spheric espace (ff) ;
	fread_be (&ome, sizeof(double), 1, ff) ;
	fread_be (&nhor, sizeof(double), 1, ff) ;
	fread_be (&vhor, sizeof(double), 1, ff) ;
	Scalar lapse(espace, ff) ;
	Vector shift(espace, ff) ;
	Metric_tensor gmet(espace, ff) ;
	fclose(ff) ;
	
	int ndom = espace.get_nbr_domains() ;
	
	// Number of point
	Base_tensor basis (shift.get_basis()) ;
	Metric_tensor gfixed (espace, COV, basis) ;
	for (int i=1 ; i<=3 ; i++)
		for (int j=1 ; j<=3 ; j++) 
			if (i==j)
				gfixed.set(i,j) = 1. ;
			else
				gfixed.set(i,j).annule_hard() ;
	gfixed.std_base() ;

  
    Vector scov (espace, COV, basis) ;
    scov = 0 ;
    scov.set(1) = 1. ;
    scov.std_base() ;

    Vector st (espace, COV, basis) ;
    st = 0 ;
    st.set(2) = 1. ;
    st.std_base() ;


    Vector er (espace, CON, basis) ;
    er = 0 ;
    er.set(1).set_domain(1) = 1. ;
    er.std_base() ;

    Vector mm (espace, CON, basis) ;
    for (int i=1 ; i<=3 ; i++)
	mm.set(i) = 0. ;
    Val_domain xx (espace.get_domain(1)->get_cart(1)) ;
    Val_domain yy (espace.get_domain(1)->get_cart(2)) ;
    mm.set(3).set_domain(1) = sqrt(xx*xx + yy*yy) ;
    mm.std_base() ;
    
    Scalar one (espace) ;
    one = 1 ; 
    one.std_base() ;
    
    List_comp used (6,2) ;
    used.set(0)->set(0) = 1 ; used.set(0)->set(1) = 1 ;
    used.set(1)->set(0) = 1 ; used.set(1)->set(1) = 2 ; 
    used.set(2)->set(0) = 1 ; used.set(2)->set(1) = 3 ;
    used.set(3)->set(0) = 2 ; used.set(3)->set(1) = 2 ;
    used.set(4)->set(0) = 2 ; used.set(4)->set(1) = 3 ; 
    used.set(5)->set(0) = 3 ; used.set(5)->set(1) = 3 ;

    List_comp dege(3,2) ;
    dege.set(0)->set(0) = 2 ; dege.set(0)->set(1) = 2 ;
    dege.set(1)->set(0) = 2 ; dege.set(1)->set(1) = 3 ; 
    dege.set(2)->set(0) = 3 ; dege.set(2)->set(1) = 3 ;

    List_comp nondege(2,2) ;
    nondege.set(0)->set(0) = 1 ; nondege.set(0)->set(1) = 2 ;
    nondege.set(1)->set(0) = 1 ; nondege.set(1)->set(1) = 3 ; 
 
    List_comp rrcomp (1,2) ;
    rrcomp.set(0)->set(0) = 1 ; rrcomp.set(0)->set(1) = 1 ;

    List_comp ozers(3,2) ;
    ozers.set(0)->set(0) = 1 ; ozers.set(0)->set(1) = 1 ;
    ozers.set(1)->set(0) = 1 ; ozers.set(1)->set(1) = 2 ; 
    ozers.set(2)->set(0) = 1 ; ozers.set(2)->set(1) = 3 ;

    double eightpi = 8.*M_PI ;
    Metric_general met (gmet) ;

	// Solve the equation in space outside the nucleus
	System_of_eqs syst (espace, 1, ndom-1) ;

	// Unknowns
	syst.add_var ("N", lapse) ;
	syst.add_var ("B", shift) ;
	met.set_system (syst, "g") ;
	syst.add_var ("ome", ome) ;

	// User defined constants
	syst.add_cst ("s", scov) ;	
	syst.add_cst ("gf", gfixed) ;
	syst.add_cst ("nhor", nhor) ;
	syst.add_cst ("vhor", vhor) ;
	
	syst.add_cst ("mm", mm) ;
	syst.add_cst ("er", er) ;
	syst.add_cst ("eightpi", eightpi) ; 
	syst.add_cst ("spin", spintarget) ; 

	// For speed
	syst.add_def ("DBS_i^j = D_i B^j") ;
	syst.add_def ("K_ij = (DBS_ij + DBS_ji) /2./N") ;

	syst.add_def (1, "sn^i = s^i / sqrt(s_i * s^i)") ;
	syst.add_def (1, "sigma = sn^i * sn^j * K_ij + D_i sn^i") ;

	// Quantities to compute the spin
	syst.add_def (1, "detq = (sn^k * s_k)^2 * determinant(g_ij)") ;
	syst.add_def (1, "intS = -sqrt(detq)*K_ij * mm^j *sn^i / eightpi") ;
	espace.add_eq_int_hor (syst, "integ(intS) = spin") ;

	syst.add_def ("V^i = g^kl * Gam_kl^i") ;
	syst.add_def ("Gauge_ij = 0.5*(D_i V_j + D_j V_i)") ;
	syst.add_def ("Ope_ij = R_ij - Gauge_ij") ;
	syst.add_def ("Hamilton = g^kl * Ope_kl - K_ij * K^ij") ;
	syst.add_def ("Momentum^i = D_j K^ij") ;
	syst.add_def ("Evol_ij = N * (Ope_ij-2*K_ik*K_j^k) - D_i D_j N + B^k * D_k K_ij + K_ik * DBS_j^k + K_jk * DBS_i^k ") ;
		
	// Inner BCs :
	syst.add_eq_bc (1, INNER_BC, "N = nhor") ; 
	syst.add_eq_bc (1, INNER_BC, "B^i = N * sn^i  + ome * mm^i") ;
	syst.add_eq_bc_exception (1, INNER_BC, "sigma=0", "g_ij * er^i * er^j = vhor") ;
	syst.add_eq_bc (1, INNER_BC, "g_ij = gf_ij", nondege) ;

	// Bulk
	for (int d=1 ; d<ndom ; d++) {
		syst.add_eq_inside (d, "Hamilton=0") ;
		syst.add_eq_inside (d, "Momentum^i=0") ;
		
		if (d==1) {
			syst.add_eq_inside (d, "Evol_ij=0", ozers) ;
			syst.add_eq_one_side (d, "Evol_ij=0", dege) ;
		}
		else
			syst.add_eq_inside (d, "Evol_ij=0", used) ;

		if (d!=ndom-1) {
		syst.add_eq_matching (d, OUTER_BC, "N") ;
		syst.add_eq_matching (d, OUTER_BC, "dn(N)") ;
		syst.add_eq_matching (d, OUTER_BC, "B^i") ;
		syst.add_eq_matching (d, OUTER_BC, "dn(B^i)") ;
		syst.add_eq_matching (d, OUTER_BC, "g_ij", used) ;
		syst.add_eq_matching (d, OUTER_BC, "dn(g_ij)", used) ;
		}	
	}

	// OUter BC
	syst.add_eq_bc (ndom-1, OUTER_BC, "N=1") ;
	syst.add_eq_bc (ndom-1, OUTER_BC, "B^i=0") ;
	syst.add_eq_bc (ndom-1, OUTER_BC, "g_ij = gf_ij", used) ;

	// Newton-Raphson
	double conv ;
	bool endloop = false ;
	int ite = 1 ;
	char name[100] ;
	
	
	while (!endloop) {
		endloop = syst.do_newton(1e-8, conv) ;
		if (rank==0) {
		    cout << "Newton iteration " << ite << " " << conv << endl ;
		    cout << "Omega = " << ome << endl ;
		}
		ite++ ;
	}
				
	sprintf (name, "kerr_%f_%f_%f_%d.dat", nhor, vhor, ome, espace.get_domain(1)->get_nbr_points()(0)) ;	
	//sprintf (name, "kerr.dat") ;	
	// Save
		if (rank==0) {
		FILE* ff = fopen(name, "w") ;
		espace.save(ff) ;
		fwrite_be (&ome, sizeof(double), 1, ff) ;
		fwrite_be (&nhor, sizeof(double), 1, ff) ;
		fwrite_be (&vhor, sizeof(double), 1, ff) ;
		lapse.save(ff) ;
		shift.save(ff) ;
		gmet.save(ff) ;
		fclose(ff) ;
		}
	
	MPI_Finalize() ;

	return EXIT_SUCCESS ;
}