Implementation of Reissner-Nordstrom metric in gr.

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
 
 !------------------------------------------------------------------

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))

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

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)

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
 !----------------------------------------------------------------

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
 

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
 

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