Commit 9225abd2 authored by m.tavelli's avatar m.tavelli

Bugfix in TBB implementation of Metric

parent 95af9357
......@@ -177,7 +177,7 @@ END SUBROUTINE NSTOV_x
SUBROUTINE NSTOV_rbar(r,NSTOV_q)
RECURSIVE SUBROUTINE NSTOV_rbar(r,NSTOV_q)
USE MainVariables, ONLY : EQN,myrank,NSTOVVar,NSTOVVar_barNew,NSTOVVar_bar, NSTOV_nODE, NSTOV_rho_c, NSTOV_rho_atmo, NSTOV_kappa, NSTOV_p_atmo, NSTOV_nODE_p
IMPLICIT NONE
! input argument
......
......@@ -5,11 +5,12 @@
!#include "expintegrator_type.f90"
RECURSIVE SUBROUTINE InitParameters(STRLEN,PARSETUP)
USE MainVariables, ONLY: nVar , nDim, EQN, ICType
USE MainVariables, ONLY: nVar , nDim, EQN, ICType, NSTOV_rho_atmo, NSTOV_kappa,NSTOV_p_atmo
USE NSTOV_mod
IMPLICIT NONE
integer :: STRLEN
character(len=STRLEN) :: PARSETUP
real :: igamma
ICType=trim(parsetup)
print *, "****************************************************************"
......@@ -74,7 +75,9 @@ RECURSIVE SUBROUTINE InitParameters(STRLEN,PARSETUP)
EQN%CCZ4LapseType = 0 ! harmonic lapse
EQN%EinsteinAutoAux = 0
igamma = 1.0/EQN%gamma
NSTOV_rho_atmo = (NSTOV_p_atmo/NSTOV_kappa)**igamma
CALL NSTOV_Main
#else
......@@ -149,6 +152,8 @@ RECURSIVE SUBROUTINE InitialData(xGP, tGp, Q)
!
case('CCZ4TOV')
V0 = 0.0
! Compute the metric and its derivatives
CALL METRIC( xGP, alpha, gp, gm, beta, Kex, g_cov, g_contr, phi )
CALL DMETRIC( xGP, AA, BB, DD, PP )
......@@ -304,6 +309,10 @@ RECURSIVE SUBROUTINE InitialData(xGP, tGp, Q)
! The trace of the extrinsic curvature at the initial time
V0(59) = K0
!
if(V0(55)>1.0 .or. V0(55)<0.0) then
print *, V0
pause
end if
#ifdef RNSTOV
......@@ -333,7 +342,7 @@ RECURSIVE SUBROUTINE InitialData(xGP, tGp, Q)
#if defined(CCZ4EINSTEIN) || defined(CCZ4GRHD)
V0(60:64) = (/ rho, VV_cov(1:3), p /)
!print *, V0
!pause
#endif
......
SUBROUTINE METRIC ( xc, lapse, gp, gm, shift, Kex, g_cov, g_contr, phi )
USE MainVariables, ONLY: ICType, aom, Mbh, P_eps,NSTOVVar
USE MainVariables, ONLY: ICType, aom, Mbh, P_eps,NSTOVVar,NSTOV_nODE
#if defined(RNSTOV)
USE NSTOV_mod
#endif
......@@ -20,6 +20,9 @@ SUBROUTINE METRIC ( xc, lapse, gp, gm, shift, Kex, g_cov, g_contr, phi )
REAL :: lx, ly, lz, HH, SS, detg, rho, r_hz, sig, theta, a, M
REAL :: st, st2, delta, rho2, sigma, zz, V0Z4(54),V070(70)
REAL :: transl(3), xGP(3),Mbh_loc,rbar
#if defined(RNSTOV)
REAL :: qloc(NSTOV_nODE)
#endif
REAL, DIMENSION(3) :: xGP_loc,xGP_sph
!
Kex = 0.0
......@@ -326,16 +329,16 @@ SUBROUTINE METRIC ( xc, lapse, gp, gm, shift, Kex, g_cov, g_contr, phi )
r=XGP_sph(1)
!
#ifdef SPHERICAL ! then we are in the original coordinate system
CALL NSTOV_x(r,NSTOVVar%qloc)
CALL NSTOV_x(r,qloc)
#elif CYLINDRICAL
PRINT *, 'CYLINDRICAL COORDINATES NOT TESTED FOR RNSTOV'
#else
CALL NSTOV_rbar(r,NSTOVVar%qloc)
CALL NSTOV_rbar(r,qloc)
!print *, r,NSTOVVar%qloc
!pause
#endif
!
lapse = EXP(NSTOVVar%qloc(2))
lapse = EXP(qloc(2))
shift(1:3) = 0.
!gammaij(1:6) = 0.
!gammaij(1) = 1.0
......@@ -343,16 +346,16 @@ SUBROUTINE METRIC ( xc, lapse, gp, gm, shift, Kex, g_cov, g_contr, phi )
!gammaij(6) = 1.0
g_cov = 0.
#ifdef SPHERICAL
g_cov(1,1) = 1.0/(1.0-2.0*NSTOVVar%qloc(1)/r)
g_cov(1,1) = 1.0/(1.0-2.0*qloc(1)/r)
g_cov(2,2) = r**2
g_cov(3,3) = SIN(xGP_sph(2))**2*r**2
!
#else
!rbar = r !*NSTOVVar%C*exp(NSTOVVar%q(4,n))
!rbar = 0.5*r/NSTOVVar%radius*(SQRT(NSTOVVar%radius**2-2*NSTOVVar%q(1,NSTOVVar%iradius)*NSTOVVar%radius)+NSTOVVar%radius-NSTOVVar%q(1,NSTOVVar%iradius))*EXP(NSTOVVar%qloc(4)-NSTOVVar%q(4,NSTOVVar%iradius))
g_cov(1,1) = 1.0/(NSTOVVar%C**2*exp(2.0*NSTOVVar%qloc(4))) ! rbar**2/(r+1e-14)**2
g_cov(2,2) = 1.0/(NSTOVVar%C**2*exp(2.0*NSTOVVar%qloc(4))) ! rbar**2/(r+1e-14)**2
g_cov(3,3) = 1.0/(NSTOVVar%C**2*exp(2.0*NSTOVVar%qloc(4))) ! rbar**2/(r+1e-14)**2
g_cov(1,1) = 1.0/(NSTOVVar%C**2*exp(2.0*qloc(4))) ! rbar**2/(r+1e-14)**2
g_cov(2,2) = 1.0/(NSTOVVar%C**2*exp(2.0*qloc(4))) ! rbar**2/(r+1e-14)**2
g_cov(3,3) = 1.0/(NSTOVVar%C**2*exp(2.0*qloc(4))) ! rbar**2/(r+1e-14)**2
!
#endif
g_contr=0.
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment