From 5751f0c658f613bf352518c075fd368cf46237e2 Mon Sep 17 00:00:00 2001 From: Angelos Karakonstantakis Date: Wed, 8 Oct 2025 10:28:10 +0200 Subject: [PATCH] Implementation of Reissner-Nordstrom metric Implemented the Reissner-Nordstrom metric and its derivatives in Cartesian-like and spherical coordinates - Define a dummy charge parameter in all metrices So the compiler won't complain for unused variables. - Integration tests for RN: CO in various charges. --- src/main/checksetup.f90 | 13 +- src/main/extern_gr.f90 | 4 +- src/main/externalforces.f90 | 1 + src/main/externalforces_gr.f90 | 8 +- src/main/metric_et.f90 | 1 + src/main/metric_kerr-schild.f90 | 1 + src/main/metric_kerr.f90 | 3 +- src/main/metric_minkowski.f90 | 1 + src/main/metric_rn.f90 | 497 ++++++++++++++++++++++++++++++ src/main/metric_schwarzschild.f90 | 1 + src/main/metric_tools.f90 | 3 +- src/setup/setup_grtde.f90 | 12 +- src/setup/setup_testparticles.f90 | 2 +- src/tests/test_gr.f90 | 73 ++++- 14 files changed, 601 insertions(+), 19 deletions(-) create mode 100644 src/main/metric_rn.f90 diff --git a/src/main/checksetup.f90 b/src/main/checksetup.f90 index dffa91ff9..759d35116 100644 --- a/src/main/checksetup.f90 +++ b/src/main/checksetup.f90 @@ -932,7 +932,8 @@ end subroutine check_setup_dustfrac !+ !------------------------------------------------------------------ subroutine check_gr(npart,nerror,xyzh,vxyzu) - use metric_tools, only:pack_metric,unpack_metric + use metric_tools, only:pack_metric,unpack_metric,imet_rn,imetric + use metric, only:charge,mass1 use utils_gr, only:get_u0 use part, only:isdead_or_accreted,ien_type,ien_entropy,ien_etotal,ien_entropy_s use units, only:in_geometric_units,get_G_code,get_c_code @@ -977,6 +978,16 @@ subroutine check_gr(npart,nerror,xyzh,vxyzu) nerror = nerror + 1 endif + if (imetric==imet_rn) then + if (abs(mass1-1.) > tiny(0.)) then + print*, ' mass1 in code units shall be unity for proper interpretation' + nerror = nerror + 1 + endif + elseif (abs(charge) > 0.) then + print*,' charge should be zero for this metric' + nerror = nerror + 1 + endif + end subroutine check_gr !------------------------------------------------------------------ diff --git a/src/main/extern_gr.f90 b/src/main/extern_gr.f90 index e3c73efc3..1dd96abfd 100644 --- a/src/main/extern_gr.f90 +++ b/src/main/extern_gr.f90 @@ -118,7 +118,7 @@ end subroutine get_grforce_all !--------------------------------------------------------------------------- subroutine dt_grforce(xyzh,fext,dtf) use physcon, only:pi - use metric_tools, only:imetric,imet_schwarzschild,imet_kerr + use metric_tools, only:imetric,imet_schwarzschild,imet_kerr,imet_rn use metric, only:mass1 real, intent(in) :: xyzh(4),fext(3) real, intent(out) :: dtf @@ -133,7 +133,7 @@ subroutine dt_grforce(xyzh,fext,dtf) endif select case (imetric) - case (imet_schwarzschild,imet_kerr) + case (imet_schwarzschild,imet_kerr,imet_rn) r2 = xyzh(1)*xyzh(1) + xyzh(2)*xyzh(2) + xyzh(3)*xyzh(3) r = sqrt(r2) omega = sqrt(mass1/(r2*r)) diff --git a/src/main/externalforces.f90 b/src/main/externalforces.f90 index bf7600efb..4cf0f32af 100644 --- a/src/main/externalforces.f90 +++ b/src/main/externalforces.f90 @@ -46,6 +46,7 @@ module externalforces logical, public :: extract_iextern_from_hdr = .false. public :: mass1,a + real, public :: charge ! ! enumerated list of external forces diff --git a/src/main/externalforces_gr.f90 b/src/main/externalforces_gr.f90 index e22e02356..b54af912f 100644 --- a/src/main/externalforces_gr.f90 +++ b/src/main/externalforces_gr.f90 @@ -19,7 +19,7 @@ module externalforces ! :Dependencies: dump_utils, infile_utils, io, metric, metric_tools, part, ! units ! - use metric, only:mass1,a + use metric, only:mass1,a,charge implicit none private @@ -36,7 +36,7 @@ module externalforces ! integer, parameter, public :: iext_gr = 1 - public :: mass1,a ! exported from metric module + public :: mass1,a,charge ! exported from metric module real, public :: accradius1 = 0. real, public :: accradius1_hard = 0. real, public :: accretedmass1 = 0. @@ -186,7 +186,7 @@ end subroutine update_externalforce !+ !----------------------------------------------------------------------- subroutine accrete_particles(iexternalforce,xi,yi,zi,hi,mi,ti,accreted,i) - use metric_tools, only:imet_minkowski,imet_schwarzschild,imet_kerr,imetric,imet_binarybh + use metric_tools, only:imet_minkowski,imet_schwarzschild,imet_kerr,imet_rn,imetric,imet_binarybh use part, only:set_particle_type,iboundary,maxphase,maxp,igas,npartoftype use metric, only:accrete_particles_metric integer, intent(in) :: iexternalforce @@ -202,7 +202,7 @@ subroutine accrete_particles(iexternalforce,xi,yi,zi,hi,mi,ti,accreted,i) case(imet_minkowski) if (first) print*,"WARNING: Accrete particles: but Metric = Minkowski" - case(imet_schwarzschild,imet_kerr) + case(imet_schwarzschild,imet_kerr,imet_rn) r2 = xi*xi + yi*yi + zi*zi if (accradius1>accradius1_hard .and. r2 < accradius1**2 .and. maxphase==maxp .and. present(i)) then call set_particle_type(i,iboundary) diff --git a/src/main/metric_et.f90 b/src/main/metric_et.f90 index c54f5fb3a..29b561fa8 100644 --- a/src/main/metric_et.f90 +++ b/src/main/metric_et.f90 @@ -23,6 +23,7 @@ module metric ! Not used anywhere in the code - Needs a fix! real, public :: mass1 = 1. ! mass of central object real, public :: a = 0.0 ! spin of central object + real, public :: charge= 0. ! charge of central object contains !---------------------------------------------------------------- diff --git a/src/main/metric_kerr-schild.f90 b/src/main/metric_kerr-schild.f90 index 26932db2f..32e4ce536 100644 --- a/src/main/metric_kerr-schild.f90 +++ b/src/main/metric_kerr-schild.f90 @@ -25,6 +25,7 @@ module metric real, public :: mass1 = 1. ! mass of central object real, public :: a = 0. ! spin of central object + real, public :: charge = 0. ! charge of central object contains diff --git a/src/main/metric_kerr.f90 b/src/main/metric_kerr.f90 index 86c0ac569..d08c497ff 100644 --- a/src/main/metric_kerr.f90 +++ b/src/main/metric_kerr.f90 @@ -26,7 +26,8 @@ module metric integer, parameter :: imetric = 3 real, public :: mass1 = 1. ! mass of central object - real, public :: a = 0.9 ! spin of central object + real, public :: a = 0.9 ! spin of central object + real, public :: charge= 0. ! charge of central object contains diff --git a/src/main/metric_minkowski.f90 b/src/main/metric_minkowski.f90 index a9615aa2c..bf5cc6e50 100644 --- a/src/main/metric_minkowski.f90 +++ b/src/main/metric_minkowski.f90 @@ -22,6 +22,7 @@ module metric integer, parameter :: imetric = 1 real, public :: mass1 = 1. ! mass of central object real, public :: a = 0. ! spin of central object + real, public :: charge= 0. ! charge of central object contains diff --git a/src/main/metric_rn.f90 b/src/main/metric_rn.f90 new file mode 100644 index 000000000..8d6aa6757 --- /dev/null +++ b/src/main/metric_rn.f90 @@ -0,0 +1,497 @@ +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2025 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.github.io/ ! +!--------------------------------------------------------------------------! +module metric +! +! None +! +! :References: None +! +! :Owner: +! - David Liptai +! - modified for Reissner-Nordstrom metric by Angelos Karakonstantakis +! +! :Runtime parameters: +! - mass1 : *black hole mass in code units* +! - charge : *black hole charge* +! +! :Dependencies: infile_utils, io +! + implicit none + character(len=*), parameter :: metric_type = 'RN' + integer, parameter :: imetric = 8 + + real, public :: mass1 = 1. ! mass of central object + real, public :: a = 0. ! spin of central object + real, public :: charge= 0.0 ! charge of central object + +contains + +!---------------------------------------------------------------- +!+ +! Metric tensors +!+ +!---------------------------------------------------------------- + +!--- The metric tensor in 'CARTESIAN-like form' +pure subroutine get_metric_cartesian(position,gcov,gcon,sqrtg) + real, intent(in) :: position(3) + real, intent(out) :: gcov(0:3,0:3) + real, intent(out), optional :: gcon(0:3,0:3) + real, intent(out), optional :: sqrtg + real :: r,r2,r3,r4,rs_on_r3,charge2_on_r4,coeff,x,y,z,x2,y2,z2,term + real :: rs, charge2 + rs = 2.*mass1 + charge2 = charge**2 + + r2 = dot_product(position,position) + r = sqrt(r2) + r3 = r*r2 + r4 = r2*r2 + rs_on_r3 = rs/r3 + charge2_on_r4 = charge2/r4 + x = position(1) + y = position(2) + z = position(3) + x2 = x**2 + y2 = y**2 + z2 = z**2 + + !--- The Reissner-Nordstrom metric tensor in CARTESIAN-like form + if (present(sqrtg)) sqrtg = 1. + + term = (1. - rs/r + charge2/r2) + coeff = 1./term + + gcov = 0. + gcov(0,0) = -term + + gcov(1,1) = coeff*(1.-(rs_on_r3 - charge2_on_r4)*(y2+z2)) + gcov(2,1) = coeff*x*y*(rs_on_r3 - charge2_on_r4) + gcov(3,1) = coeff*x*z*(rs_on_r3 - charge2_on_r4) + + gcov(1,2) = gcov(2,1) + gcov(2,2) = coeff*(1.-(rs_on_r3 - charge2_on_r4)*(x2+z2)) + gcov(3,2) = coeff*y*z*(rs_on_r3 - charge2_on_r4) + + gcov(1,3) = gcov(3,1) + gcov(2,3) = gcov(3,2) + gcov(3,3) = coeff*(1.-(rs_on_r3 - charge2_on_r4)*(x2+y2)) + + if (present(gcon)) then + gcon = 0. + gcon(0,0) = -1./term + + gcon(1,1) = 1.-(rs_on_r3 - charge2_on_r4)*x2 + gcon(2,1) = -(rs_on_r3 - charge2_on_r4)*x*y + gcon(3,1) = -(rs_on_r3 - charge2_on_r4)*x*z + + gcon(1,2) = gcon(2,1) + gcon(2,2) = 1.-(rs_on_r3 - charge2_on_r4)*y2 + gcon(3,2) = -(rs_on_r3 - charge2_on_r4)*y*z + + gcon(1,3) = gcon(3,1) + gcon(2,3) = gcon(3,2) + gcon(3,3) = 1.-(rs_on_r3 - charge2_on_r4)*z2 + endif + +end subroutine get_metric_cartesian + +pure subroutine get_metric_spherical(position,gcov,gcon,sqrtg) + real, intent(in) :: position(3) + real, intent(out) :: gcov(0:3,0:3) + real, intent(out), optional :: gcon(0:3,0:3) + real, intent(out), optional :: sqrtg + real :: r,theta,sintheta,r2 + real :: rs, charge2 + rs = 2.*mass1 + charge2 = charge**2 + + r = position(1) + r2 = r**2 + theta = position(2) + sintheta = sin(theta) + + gcov(0,0) = -1. + rs/r - charge2 / r2 + gcov(1,0) = 0. + gcov(2,0) = 0. + gcov(3,0) = 0. + gcov(0,1) = 0. + gcov(1,1) = 1./(1. - rs/r + charge2/r2) + gcov(2,1) = 0. + gcov(3,1) = 0. + gcov(0,2) = 0. + gcov(1,2) = 0. + gcov(2,2) = r2 + gcov(3,2) = 0. + gcov(0,3) = 0. + gcov(1,3) = 0. + gcov(2,3) = 0. + gcov(3,3) = r2*sintheta**2 + + if (present(gcon)) then + gcon(0,0) = -(r2/(r2 - rs*r + charge2)) + gcon(1,0) = 0. + gcon(2,0) = 0. + gcon(3,0) = 0. + gcon(0,1) = 0. + gcon(1,1) = 1. - rs/r + charge2/r2 + gcon(2,1) = 0. + gcon(3,1) = 0. + gcon(0,2) = 0. + gcon(1,2) = 0. + gcon(2,2) = 1./r2 + gcon(3,2) = 0. + gcon(0,3) = 0. + gcon(1,3) = 0. + gcon(2,3) = 0. + gcon(3,3) = 1./(r2*sintheta**2) + endif + + if (present(sqrtg)) sqrtg = r2*sintheta + +end subroutine get_metric_spherical + +!---------------------------------------------------------------- +!+ +! Metric tensors derivatives +!+ +!---------------------------------------------------------------- + +!--- Derivatives of the covariant 'CARTEISAN' metric +pure subroutine metric_cartesian_derivatives(position,dgcovdx, dgcovdy, dgcovdz) + real, intent(in) :: position(3) + real, intent(out) :: dgcovdx(0:3,0:3), dgcovdy(0:3,0:3), dgcovdz(0:3,0:3) + real :: x,y,z,r,r2,r3,r4,r5,rs_on_r3,x2,y2,z2,rs2 + real :: rs + rs = 2.*mass1 + + dgcovdx = 0. + dgcovdy = 0. + dgcovdz = 0. + x = position(1) + y = position(2) + z = position(3) + x2= x**2 + y2= y**2 + z2= z**2 + r2 = dot_product(position,position) + r = sqrt(r2) + r3 = r*r2 + r4 = r2*r2 + r5 = r*r4 + + rs_on_r3 = rs/r3 + + ! dx + dgcovdx(0,0) = x*(2*charge**2 - r*rs)/r**4 + dgcovdx(1,1) = x*((2*charge**2 - r*rs)*(r**4 + (charge**2 - r*rs)* & + (y**2 + z**2)) - (4*charge**2 - 3*r*rs)*(y**2 + z**2)*(charge**2 & + + r**2 - r*rs))/(r**4*(charge**2 + r**2 - r*rs)**2) + dgcovdx(2,2) = x*((2*charge**2 - r*rs)*(r**4 + (charge**2 - r*rs)* & + (x**2 + z**2)) + (2*r**2*(charge**2 - r*rs) - (4*charge**2 - 3*r* & + rs)*(x**2 + z**2))*(charge**2 + r**2 - r*rs))/(r**4*(charge**2 + & + r**2 - r*rs)**2) + dgcovdx(3,3) = x*((2*charge**2 - r*rs)*(r**4 + (charge**2 - r*rs)* & + (x**2 + y**2)) + (2*r**2*(charge**2 - r*rs) - (4*charge**2 - 3*r* & + rs)*(x**2 + y**2))*(charge**2 + r**2 - r*rs))/(r**4*(charge**2 + & + r**2 - r*rs)**2) + dgcovdx(1,2) = y*(-r**2*(charge**2 - r*rs)*(charge**2 + r**2 - r* & + rs) - x**2*(charge**2 - r*rs)*(2*charge**2 - r*rs) + x**2*(4* & + charge**2 - 3*r*rs)*(charge**2 + r**2 - r*rs))/(r**4*(charge**2 + & + r**2 - r*rs)**2) + dgcovdx(1,3) = z*(-r**2*(charge**2 - r*rs)*(charge**2 + r**2 - r* & + rs) - x**2*(charge**2 - r*rs)*(2*charge**2 - r*rs) + x**2*(4* & + charge**2 - 3*r*rs)*(charge**2 + r**2 - r*rs))/(r**4*(charge**2 + & + r**2 - r*rs)**2) + dgcovdx(2,3) = x*y*z*(-(charge**2 - r*rs)*(2*charge**2 - r*rs) + ( & + 4*charge**2 - 3*r*rs)*(charge**2 + r**2 - r*rs))/(r**4*(charge**2 & + + r**2 - r*rs)**2) + + dgcovdx(1,0) = 0. + dgcovdx(2,0) = 0. + dgcovdx(3,0) = 0. + dgcovdx(0,1) = 0. + dgcovdx(0,2) = 0. + dgcovdx(0,3) = 0. + dgcovdx(2,1) = dgcovdx(1,2) + dgcovdx(3,1) = dgcovdx(1,3) + dgcovdx(3,2) = dgcovdx(2,3) + + ! dy + dgcovdy(1,0) = 0. + dgcovdy(2,0) = 0. + dgcovdy(3,0) = 0. + dgcovdy(0,1) = 0. + dgcovdy(0,2) = 0. + dgcovdy(0,3) = 0. + + dgcovdy(0,0) = y*(2*charge**2 - r*rs)/r**4 + dgcovdy(1,1) = y*((2*charge**2 - r*rs)*(r**4 + (charge**2 - r*rs)* & + (y**2 + z**2)) + (2*r**2*(charge**2 - r*rs) - (4*charge**2 - 3*r* & + rs)*(y**2 + z**2))*(charge**2 + r**2 - r*rs))/(r**4*(charge**2 + & + r**2 - r*rs)**2) + dgcovdy(2,2) = y*((2*charge**2 - r*rs)*(r**4 + (charge**2 - r*rs)* & + (x**2 + z**2)) - (4*charge**2 - 3*r*rs)*(x**2 + z**2)*(charge**2 & + + r**2 - r*rs))/(r**4*(charge**2 + r**2 - r*rs)**2) + dgcovdy(3,3) = y*((2*charge**2 - r*rs)*(r**4 + (charge**2 - r*rs)* & + (x**2 + y**2)) + (2*r**2*(charge**2 - r*rs) - (4*charge**2 - 3*r* & + rs)*(x**2 + y**2))*(charge**2 + r**2 - r*rs))/(r**4*(charge**2 + & + r**2 - r*rs)**2) + dgcovdy(1,2) = x*(-r**2*(charge**2 - r*rs)*(charge**2 + r**2 - r* & + rs) - y**2*(charge**2 - r*rs)*(2*charge**2 - r*rs) + y**2*(4* & + charge**2 - 3*r*rs)*(charge**2 + r**2 - r*rs))/(r**4*(charge**2 + & + r**2 - r*rs)**2) + dgcovdy(1,3) = x*y*z*(-(charge**2 - r*rs)*(2*charge**2 - r*rs) + ( & + 4*charge**2 - 3*r*rs)*(charge**2 + r**2 - r*rs))/(r**4*(charge**2 & + + r**2 - r*rs)**2) + dgcovdy(2,3) = z*(-r**2*(charge**2 - r*rs)*(charge**2 + r**2 - r* & + rs) - y**2*(charge**2 - r*rs)*(2*charge**2 - r*rs) + y**2*(4* & + charge**2 - 3*r*rs)*(charge**2 + r**2 - r*rs))/(r**4*(charge**2 + & + r**2 - r*rs)**2) + + dgcovdy(2,1) = dgcovdy(1,2) + dgcovdy(3,1) = dgcovdy(1,3) + dgcovdy(3,2) = dgcovdy(2,3) + + ! dz + dgcovdz(1,0) = 0. + dgcovdz(2,0) = 0. + dgcovdz(3,0) = 0. + dgcovdz(0,1) = 0. + dgcovdz(0,2) = 0. + dgcovdz(0,3) = 0. + + dgcovdz(0,0) = z*(2*charge**2 - r*rs)/r**4 + dgcovdz(1,1) = z*((2*charge**2 - r*rs)*(r**4 + (charge**2 - r*rs)* & + (y**2 + z**2)) + (2*r**2*(charge**2 - r*rs) - (4*charge**2 - 3*r* & + rs)*(y**2 + z**2))*(charge**2 + r**2 - r*rs))/(r**4*(charge**2 + & + r**2 - r*rs)**2) + dgcovdz(2,2) = z*((2*charge**2 - r*rs)*(r**4 + (charge**2 - r*rs)* & + (x**2 + z**2)) + (2*r**2*(charge**2 - r*rs) - (4*charge**2 - 3*r* & + rs)*(x**2 + z**2))*(charge**2 + r**2 - r*rs))/(r**4*(charge**2 + & + r**2 - r*rs)**2) + dgcovdz(3,3) = z*((2*charge**2 - r*rs)*(r**4 + (charge**2 - r*rs)* & + (x**2 + y**2)) - (4*charge**2 - 3*r*rs)*(x**2 + y**2)*(charge**2 & + + r**2 - r*rs))/(r**4*(charge**2 + r**2 - r*rs)**2) + dgcovdz(1,2) = x*y*z*(-(charge**2 - r*rs)*(2*charge**2 - r*rs) + ( & + 4*charge**2 - 3*r*rs)*(charge**2 + r**2 - r*rs))/(r**4*(charge**2 & + + r**2 - r*rs)**2) + dgcovdz(1,3) = x*(-r**2*(charge**2 - r*rs)*(charge**2 + r**2 - r* & + rs) - z**2*(charge**2 - r*rs)*(2*charge**2 - r*rs) + z**2*(4* & + charge**2 - 3*r*rs)*(charge**2 + r**2 - r*rs))/(r**4*(charge**2 + & + r**2 - r*rs)**2) + dgcovdz(2,3) = y*(-r**2*(charge**2 - r*rs)*(charge**2 + r**2 - r* & + rs) - z**2*(charge**2 - r*rs)*(2*charge**2 - r*rs) + z**2*(4* & + charge**2 - 3*r*rs)*(charge**2 + r**2 - r*rs))/(r**4*(charge**2 + & + r**2 - r*rs)**2) + + dgcovdz(2,1) = dgcovdz(1,2) + dgcovdz(3,1) = dgcovdz(1,3) + dgcovdz(3,2) = dgcovdz(2,3) + +end subroutine metric_cartesian_derivatives + +!--- Derivatives of the covariant 'SPHERICAL' metric +pure subroutine metric_spherical_derivatives(position,dgcovdr, dgcovdtheta, dgcovdphi) + real, intent(in) :: position(3) + real, intent(out), dimension(0:3,0:3) :: dgcovdr,dgcovdtheta,dgcovdphi + real :: r, theta, r2, r3 + real :: rs, charge2, charge2_on_r + rs = 2.*mass1 + charge2 = charge**2 + + r = position(1) + theta = position(2) + + r2 = r*r + r3 = r2*r + charge2_on_r = charge2 / r + + dgcovdr = 0. + dgcovdtheta = 0. + dgcovdphi = 0. + + dgcovdr(0,0) = -rs/r**2 + 2. * charge2 / r3 + dgcovdr(1,1) = (2.*charge2_on_r-rs)/(r-rs+charge2_on_r)**2 + dgcovdr(2,2) = 2.*r + dgcovdr(3,3) = 2.*r*sin(theta)**2 + + dgcovdtheta(3,3) = 2.*r**2*cos(theta)*sin(theta) + +end subroutine metric_spherical_derivatives + +!---------------------------------------------------------------- +!+ +! Coordinate transformations +!+ +!---------------------------------------------------------------- + +!--- (Jacobian tensor) Derivatives of Schwarzschild 'Spherical' with respect to 'Cartesian' coordinates +pure subroutine get_jacobian(position,dxdx) + real, intent(in), dimension(3) :: position + real, intent(out), dimension(0:3,0:3) :: dxdx + real, dimension(3) :: dSPHERICALdx,dSPHERICALdy,dSPHERICALdz + real :: drdx,drdy,drdz + real :: dthetadx,dthetady,dthetadz + real :: dphidx,dphidy,dphidz + real :: x,y,z,x2,y2,z2,r2,r,rcyl2,rcyl + + x = position(1) + y = position(2) + z = position(3) + x2 = x**2 + y2 = y**2 + z2 = z**2 + r2 = x2+y2+z2 + r = sqrt(r2) + rcyl2 = x2+y2 + rcyl = sqrt(x2+y2) + + drdx = x/r + drdy = y/r + drdz = z/r + + dthetadx = x*z/(r2*rcyl) + dthetady = y*z/(r2*rcyl) + dthetadz = -rcyl/r2 + + dphidx = -y/(x2+y2) + dphidy = x/(x2+y2) + dphidz = 0. + + dSPHERICALdx=(/drdx,dthetadx,dphidx/) + dSPHERICALdy=(/drdy,dthetady,dphidy/) + dSPHERICALdz=(/drdz,dthetadz,dphidz/) + + dxdx = 0. + dxdx(0,0) = 1. + dxdx(1:3,1) = dSPHERICALdx + dxdx(1:3,2) = dSPHERICALdy + dxdx(1:3,3) = dSPHERICALdz + +end subroutine get_jacobian + +pure subroutine cartesian2spherical(xcart,xspher) + real, intent(in) :: xcart(3) + real, intent(out) ::xspher(3) + real :: x,y,z + real :: r,theta,phi + + x = xcart(1) + y = xcart(2) + z = xcart(3) + + r = sqrt(x**2+y**2+z**2) + theta = acos(z/r) + phi = atan2(y,x) + + xspher = (/r,theta,phi/) + +end subroutine cartesian2spherical + +pure subroutine spherical2cartesian(xspher,xcart) + real, intent(in) :: xspher(3) + real, intent(out) :: xcart(3) + real :: x,y,z,r,theta,phi + + r = xspher(1) + theta = xspher(2) + phi = xspher(3) + x = r*sin(theta)*cos(phi) + y = r*sin(theta)*sin(phi) + z = r*cos(theta) + + xcart = (/x,y,z/) + +end subroutine spherical2cartesian + +!----------------------------------------------------------------------- +!+ +! writes metric options to the input file +!+ +!----------------------------------------------------------------------- +subroutine write_options_metric(iunit) + use infile_utils, only:write_inopt + integer, intent(in) :: iunit + + write(iunit,"(/,a)") '# options relating to the '//trim(metric_type)//' metric' + + call write_inopt(mass1,'mass1','black hole mass in code units',iunit) + call write_inopt(charge,'charge','charge parameter for Reissner-Nordstrom metric',iunit) + +end subroutine write_options_metric + +!----------------------------------------------------------------------- +!+ +! reads metric options from the input file +!+ +!----------------------------------------------------------------------- +subroutine read_options_metric(db,nerr) + use io, only:warn + use infile_utils, only:inopts,read_inopt + type(inopts), intent(inout) :: db(:) + integer, intent(inout) :: nerr + + call read_inopt(mass1,'mass1',db,errcount=nerr,min=0.,max=1.e12) + if (mass1 <= tiny(mass1)) call warn('metric','black hole mass: mass1 = 0') + call read_inopt(charge,'charge',db,errcount=nerr,min=0.,max=1.e12) + +end subroutine read_options_metric + +!------------------------------------------------------------------------------- +!+ +! Subroutine to update the metric inputs if time dependent +!+ +!------------------------------------------------------------------------------- +subroutine update_metric(time) + real, intent(in) :: time + +end subroutine update_metric + +!----------------------------------------------------------------------- +!+ +! Check if a particle should be accreted by the black hole +!+ +!----------------------------------------------------------------------- +subroutine accrete_particles_metric(xi,yi,zi,mi,ti,accradius,accreted) + real, intent(in) :: xi,yi,zi,mi,ti,accradius + logical, intent(out) :: accreted + + accreted = .false. + +end subroutine accrete_particles_metric + +!----------------------------------------------------------------------- +!+ +! writes relevant options to the header of the dump file +!+ +!----------------------------------------------------------------------- +subroutine write_headeropts_metric(hdr,time,accradius,ierr) + use dump_utils, only:dump_h + type(dump_h), intent(inout) :: hdr + real, intent(in) :: time,accradius + integer, intent(out) :: ierr + + ierr = 0 + +end subroutine write_headeropts_metric + +!----------------------------------------------------------------------- +!+ +! reads relevant options from the header of the dump file +!+ +!----------------------------------------------------------------------- +subroutine read_headeropts_metric(hdr,ierr) + use dump_utils, only:dump_h + type(dump_h), intent(in) :: hdr + integer, intent(out) :: ierr + + ierr = 0 + +end subroutine read_headeropts_metric + +end module metric diff --git a/src/main/metric_schwarzschild.f90 b/src/main/metric_schwarzschild.f90 index fcd1b8896..83c7e68e5 100644 --- a/src/main/metric_schwarzschild.f90 +++ b/src/main/metric_schwarzschild.f90 @@ -24,6 +24,7 @@ module metric real, public :: mass1 = 1. ! mass of central object real, public :: a = 0. ! spin of central object + real, public :: charge= 0. ! charge of central object contains diff --git a/src/main/metric_tools.f90 b/src/main/metric_tools.f90 index 256042954..1970023b4 100644 --- a/src/main/metric_tools.f90 +++ b/src/main/metric_tools.f90 @@ -37,7 +37,8 @@ module metric_tools imet_kerrschild = 4, & ! Kerr metric, Kerr-Schild coordinates imet_binarybh = 5, & ! Binary black hole metric imet_flrw = 6, & ! Friedmann-Lemaître-Robertson-Walker metric - imet_et = 7 ! Tabulated metric from Einstein toolkit + imet_et = 7, & ! Tabulated metric from Einstein toolkit + imet_rn = 8 ! Reissner-Nordstrom metric !--- Choice of coordinate system ! (When using this with PHANTOM, it should always be set to cartesian) diff --git a/src/setup/setup_grtde.f90 b/src/setup/setup_grtde.f90 index ed2ec19ae..dc94faca2 100644 --- a/src/setup/setup_grtde.f90 +++ b/src/setup/setup_grtde.f90 @@ -44,7 +44,7 @@ module setup use setstar, only:star_t use setorbit, only:orbit_t - use externalforces, only:mass1,a + use externalforces, only:mass1,a,charge implicit none public :: setpart @@ -99,7 +99,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, real, intent(out) :: massoftype(:) real, intent(out) :: polyk,gamma,hfact real, intent(inout) :: time - character(len=20), intent(in) :: fileprefix + character(len=*), intent(in) :: fileprefix real, intent(out) :: vxyzu(:,:) integer :: ierr,np_default integer :: nptmass_in @@ -139,6 +139,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, mass_unit = '1.e6*solarm' racc = '6.' mhole = 1.e6 ! (solar masses) + charge = 0. call set_units(mass=mhole*solarm,c=1.d0,G=1.d0) !--Set central mass to M=1 in code units call set_defaults_stars(star) call set_defaults_orbit(orbit) @@ -219,8 +220,9 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, if (ierr /= 0) call fatal('setup','could not convert racc to code units') accradius1 = accradius1_hard - a = 0. - theta_bh = theta_bh*pi/180. + a = 0. + charge = 0. + theta_bh = theta_bh*pi/180. if (id==master) then print "(4(/,1x,a,1f10.3))",' black hole mass: ',mass1, & @@ -381,6 +383,7 @@ subroutine write_setupfile(filename) write(iunit,"(/,a)") '# options for central object' call write_inopt(mhole, 'mhole', 'mass of black hole (solar mass)', iunit) + call write_inopt(charge, 'charge', 'charge of black hole', iunit) call write_inopt(racc, 'racc', 'accretion radius for the central object (code units or e.g. 1*km)', iunit) call write_options_stars(star,relax,write_profile,ieos,iunit,nstar) @@ -441,6 +444,7 @@ subroutine read_setupfile(filename,ierr) !--read black hole mass in solar masses ! call read_inopt(mhole,'mhole',db,min=0.,errcount=nerr) + call read_inopt(charge,'charge',db,errcount=nerr) call read_inopt(racc, 'racc', db,errcount=nerr) mass1 = mhole*solarm/umass ! diff --git a/src/setup/setup_testparticles.f90 b/src/setup/setup_testparticles.f90 index 216bc8d6a..fb3d8a3e3 100644 --- a/src/setup/setup_testparticles.f90 +++ b/src/setup/setup_testparticles.f90 @@ -36,7 +36,7 @@ module setup ! Module variables for setup parameters integer :: orbtype integer :: dumpsperorbit - real :: spin + real :: spin, charge real :: r real :: norbits real :: x0, y0, z0 diff --git a/src/tests/test_gr.f90 b/src/tests/test_gr.f90 index 1e69f8c95..c7b22870c 100644 --- a/src/tests/test_gr.f90 +++ b/src/tests/test_gr.f90 @@ -42,30 +42,32 @@ subroutine test_gr(ntests,npass) call test_combinations_all(ntests,npass) call test_precession(ntests,npass) call test_inccirc(ntests,npass) + call test_rn_charged(ntests,npass) if (id==master) write(*,"(/,a)") '<-- GR TESTS COMPLETE' end subroutine test_gr !----------------------------------------------------------------------- !+ -! Test of orbital precession in the Kerr metric +! Test of orbital precession in the Kerr/Schwarzschild/RN metric !+ !----------------------------------------------------------------------- subroutine test_precession(ntests,npass) - use metric_tools, only:imetric,imet_kerr,imet_schwarzschild - use metric, only:a + use metric_tools, only:imetric,imet_kerr,imet_schwarzschild,imet_rn + use metric, only:a,charge integer, intent(inout) :: ntests,npass integer :: nerr(6),norbits,nstepsperorbit,itest real :: dt,period,x0,vy0,tmax,angtol,postol real :: angmom(3),angmom0(3),xyz(3),vxyz(3) write(*,'(/,a)') '--> testing substep_gr (precession)' - if (imetric /= imet_kerr .and. imetric /= imet_schwarzschild) then - write(*,'(/,a)') ' Skipping test! Metric is not Kerr (or Schwarzschild).' + if (imetric /= imet_kerr .and. imetric /= imet_schwarzschild .and. imetric /= imet_rn) then + write(*,'(/,a)') ' Skipping test! Metric is not Kerr, Schwarzschild or RN.' return endif a = 0. + charge = 0. x0 = 90. vy0 = 0.0521157 period = 2390. ! approximate @@ -603,4 +605,65 @@ subroutine test_cons2prim_i(x,v,dens,u,p,ncheck,nfail,errmax,tol) end subroutine test_cons2prim_i +!----------------------------------------------------------------------- +!+ +! Test of circular orbits in the Reissner-Nordstrom metric +! with various charges (Q/M = 0, 0.5, 0.9, 1.1) +!+ +!----------------------------------------------------------------------- +subroutine test_rn_charged(ntests,npass) + use io, only:id,master + use metric_tools, only:imetric,imet_rn + use metric, only:mass1,charge + integer, intent(inout) :: ntests,npass + integer :: nerr(6),itest,iq + real :: dt,period,tmax,r0,v0,mass_val,charge_val + real :: angmom(3),angmom0(3),xyz(3),vxyz(3) + real :: charges_to_test(4) + real, parameter :: postol = 1.e-5, angtol = 1.e-10 + + if (id==master) write(*,'(/,a)') '--> testing RN charged orbits' + + if (imetric /= imet_rn) then + if (id==master) write(*,'(/,a)') ' Skipping test! Metric is not RN.' + return + endif + + ! Q/M ratios to test: Schwarzschild-like, Sub-extremal, Near-extremal, Naked Singularity + charges_to_test = (/0.0, 0.5, 0.9, 1.1/) + mass_val = 1.0 + r0 = 10.0 + + do iq=1,size(charges_to_test) + charge_val = charges_to_test(iq) + charge = charge_val + mass1 = mass_val + + ! Circular orbit coordinate velocity for RN metric: + ! v_circ = sqrt( M/r - Q^2/r^2 ) + v0 = sqrt( mass_val/r0 - charge_val**2/r0**2 ) + + period = 2.0 * 3.14159265358979 * r0 / v0 + tmax = 1.0 * period + dt = period / 1000.0 + + if (id==master) write(*,'(a,F4.1,a,F4.2)') ' testing Q/M = ', charge_val/mass_val, ' v0 = ', v0 + + do itest=1,2 + xyz = (/r0, 0., 0./) + vxyz = (/0., v0, 0./) + call integrate_geodesic(tmax,dt,xyz,vxyz,angmom0,angmom,use_sink=(itest==2)) + + nerr = 0 + call checkval(angmom(1),angmom0(1),angtol,nerr(1),'error in angmomx') + call checkval(angmom(2),angmom0(2),angtol,nerr(2),'error in angmomy') + call checkval(angmom(3),angmom0(3),angtol,nerr(3),'error in angmomz') + call checkval(sqrt(dot_product(xyz,xyz)), r0, postol, nerr(4),'error in final r position') + + call update_test_scores(ntests,nerr,npass) + enddo + enddo + +end subroutine test_rn_charged + end module testgr