Merge pull request #480 from spencermagnall/master

PhantomNR

AUTHORS

@@ -29,6 +29,7 @@ Simone Ceppi <simone.ceppi@unimi.it>
 MatsEsseldeurs <matsesseldeurs@yahoo.com>
 Caitlyn Hardiman <caitlyn.hardiman1@monash.edu>
 Enrico Ragusa <enr.ragusa@gmail.com>
+Spencer Magnall <spencermagnall@gmail.com>
 fhu <fhuu0005@student.monash.edu>
 Sergei Biriukov <svbiriukov@gmail.com>
 Cristiano Longarini <cristiano.longarini@unimi.it>
@@ -54,11 +55,11 @@ Kieran Hirsh <kieran.hirsh1@monash.edu>
 Nicole Rodrigues <nicole.rodrigues@monash.edu>
 Amena Faruqi <Amena.Faruqi@warwick.ac.uk>
 David Trevascus <dtre10@student.monash.edu>
+Farzana Meru <f.meru@warwick.ac.uk>
 Chris Nixon <cjn@leicester.ac.uk>
 Megha Sharma <megha@meghas-air.home>
 Nicolas Cuello <cuellonicolas@gmail.com>
 Benoit Commercon <benoit.commercon@gmail.com>
-Farzana Meru <f.meru@warwick.ac.uk>
 Giulia Ballabio <giulia.ballabio2@studenti.unimi.it>
 Joe Fisher <jwfis1@student.monash.edu>
 Maxime Lombart <maxime.lombart@ens-lyon.fr>
@@ -70,6 +71,7 @@ s-neilson <36410751+s-neilson@users.noreply.github.com>
 Alison Young <ayoung@astro.ex.ac.uk>
 Cox, Samuel <sc676@leicester.ac.uk>
 Jorge Cuadra <jcuadra@astro.puc.cl>
+Miguel Gonzalez-Bolivar <miguelgb@astro.unam.mx>
 Nicolás Cuello <cuello.nicolas@gmail.com>
 Steven Rieder <steven@rieder.nl>
 Stéven Toupin <steven.toupin@gmail.com>

build/Makefile

@@ -494,7 +494,7 @@ ifdef METRIC
 else
     SRCMETRIC= metric_minkowski.f90
 endif
-SRCGR=inverse4x4.f90 $(SRCMETRIC) metric_tools.f90 utils_gr.f90
+SRCGR=inverse4x4.f90 einsteintk_utils.f90 $(SRCMETRIC) metric_tools.f90 utils_gr.f90 interpolate3D.f90 tmunu2grid.f90 
 #
 # chemistry and cooling
 #
@@ -547,6 +547,10 @@ SOURCES= physcon.f90 ${CONFIG} ${SRCKERNEL} io.F90 units.f90 \
          mf_write.f90 evolve.F90 utils_orbits.f90 utils_linalg.f90 \
          checksetup.F90 initial.F90
 
+# Needed as einsteintk_wrapper depends on initial
+ifeq ($(GR),yes)
+	SOURCES+=einsteintk_wrapper.f90
+endif
 OBJECTS1 = $(SOURCES:.f90=.o)
 OBJECTS = $(OBJECTS1:.F90=.o)
 

build/Makefile_setups

@@ -1036,6 +1036,25 @@ ifeq ($(SETUP), testgr)
     SETUPFILE= setup_grdisc.f90
 endif
 
+ifeq ($(SETUP), flrw)
+#   constant density FLRW cosmology with perturbations  
+ GR=yes
+ KNOWN_SETUP=yes
+ IND_TIMESTEPS=no
+ METRIC=et
+ SETUPFILE= setup_flrw.f90
+ PERIODIC=yes
+endif
+
+ifeq ($(SETUP), flrwpspec)
+#  FLRW universe using a CMB powerspectrum and the Zeldovich approximation
+ GR=yes
+ KNOWN_SETUP=yes
+ IND_TIMESTEPS=no
+ METRIC=et
+ SETUPFILE= setup_flrwpspec.f90
+ PERIODIC=yes
+endif
 ifeq ($(SETUP), default)
 #   default setup, uniform box
     KNOWN_SETUP=yes

docs/phantomNR.rst

@@ -0,0 +1,79 @@
+PhantomNR
+=========
+
+Using PhantomNR to simulate general relativistic hydrodynamics on dynamical spacetimes 
+--------------------------------------------------------------------------------------
+
+About phantomNR
+~~~~~~~~~~~~~~~
+
+`phantomNR <https://github.com/spencermagnall/phantomNR>`__ is
+an extension to the General Relativistic Smoothed Particle Hydrodynamics code Phantom,
+that allows for the evolution of relativistic fluids with evolving spacetime metrics. 
+This is acomplished via coupling with the numerical relativity framework Einstein Toolkit (ET).
+phantomNR's current usage is as a fully relativistic  N-Body code for the simulation of inhomogenous
+cosmologies (see `Magnall et al. 2023 <https://ui.adsabs.harvard.edu/abs/2023arXiv230715194M/abstract>`__). 
+Einstein Toolkit acts as a "driver" for both the spacetime evolution, and the hydrodynamic evolution. 
+As a consquence, simulations are started and mointered entirely within ET, and are setup using a .par 
+parameter file which describes the parameters of the simulation. In addition, phantomNR also requires
+particle information, which is provided via the standard phantom dump file.   
+
+
+Compilation and linking
+~~~~~~~~~~~~~~~~~~~~~~~ 
+You will first need to compile phantom and phantomsetup 
+using the flrw setup
+
+::
+
+   scripts/writemake.sh flrw > Makefile 
+
+   make; make setup
+
+which compiles the libphantom.a static library which is 
+required for linking and the phantom and phantomsetup binaries. 
+
+You will also need to set the include directory of phantom in Einstein Toolkit
+e.g: 
+
+::
+
+   PHANTOM_DIR = /Users/smag0001/phantom/phantomET/bin
+
+Generating a phantom dump file from phantom setup
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Particles can be setup using phantomsetup in two ways:
+
+1. **Using a regular .setup file** 
+   e.g ./phantomsetup flrw.setup will produce a dump file and .in file using an interactive setup routine.
+
+
+3. **Using a .par file** By appending .setup to the end of an Einstein Toolkit parameter file, phantomsetup 
+   will automatically read in (most) relevant quantities about the simulation setup and generate an appropriate
+   distribution of particles
+
+
+Troubleshooting
+---------------
+
+**Issue**: Large Constraint Violations
+
+
+
+**Solution**: Generally, this is indicative of a mismatch between the spacetime setup by Einstein Toolkit
+and the particle distribution which is setup by Phantom. A large raw constraint violation, may not always be indicative
+of a poorly initialised setup however. It is important to check the relative constraint violations (TODO insert equations)
+
+In many cases, a poor initial constraint is simply a consquence of not setting spacetime and consistently (e.g phi=1e-4 for particles, but phi=1e-6 for spacetime).
+We reccomend that the .in and dumpfiles are generated using the .par file of Einstein Toolkit to alleviate this issue.
+
+Constraint violations may also occur due to a low particle and/or grid resolution 
+
+
+
+
+Using phantomNR on Ozstar/NT 
+-------------------------------
+
+

src/main/config.F90

@@ -270,6 +270,15 @@ module dim
  logical, parameter :: gr = .false.
 #endif
 
+!---------------------
+! Numerical relativity
+!---------------------
+#ifdef NR
+ logical, parameter :: nr = .true.
+#else
+ logical, parameter :: nr = .false.
+#endif  
+
 !--------------------
 ! Supertimestepping
 !--------------------

src/main/cons2primsolver.f90

@@ -75,7 +75,7 @@ end subroutine get_u
 !+
 !----------------------------------------------------------------
 subroutine primitive2conservative(x,metrici,v,dens,u,P,rho,pmom,en,ien_type)
- use utils_gr,     only:get_u0
+ use utils_gr,     only:get_u0,get_sqrtg
  use metric_tools, only:unpack_metric
  use io,           only:error
  use eos,          only:gmw,get_entropy
@@ -93,9 +93,9 @@ subroutine primitive2conservative(x,metrici,v,dens,u,P,rho,pmom,en,ien_type)
 
  enth = 1. + u + P/dens
 
- ! Hard coded sqrtg=1 since phantom is always in cartesian coordinates
- sqrtg = 1.
+ ! get determinant of metric
  call unpack_metric(metrici,gcov=gcov)
+ call get_sqrtg(gcov,sqrtg)
 
  call get_u0(gcov,v,U0,ierror)
  if (ierror > 0) call error('get_u0 in prim2cons','1/sqrt(-v_mu v^mu) ---> non-negative: v_mu v^mu')
@@ -134,6 +134,7 @@ end subroutine primitive2conservative
 !+
 !----------------------------------------------------------------
 subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ierr,ien_type)
+ use utils_gr,     only:get_sqrtg,get_sqrt_gamma
  use metric_tools, only:unpack_metric
  use eos,          only:ieos,gmw,get_entropy,get_p_from_rho_s,gamma_global=>gamma
  use io,           only:fatal
@@ -147,21 +148,21 @@ subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ie
  integer, intent(in)  :: ien_type
  real, dimension(1:3,1:3) :: gammaijUP
  real :: sqrtg,sqrtg_inv,lorentz_LEO,pmom2,alpha,betadown(1:3),betaUP(1:3),enth_old,v3d(1:3)
- real :: f,df,term,lorentz_LEO2,gamfac,pm_dot_b,sqrt_gamma_inv,enth,gamma1
+ real :: f,df,term,lorentz_LEO2,gamfac,pm_dot_b,sqrt_gamma_inv,enth,gamma1,sqrt_gamma
  real(kind=8) :: cgsdens,cgsu
  integer :: niter, i
  real, parameter :: tol = 1.e-12
  integer, parameter :: nitermax = 100
  logical :: converged
+ real    :: gcov(0:3,0:3)
  ierr = 0
 
- ! Hard coding sqrgt=1 since phantom is always in cartesian coordinates
- sqrtg = 1.
- sqrtg_inv = 1./sqrtg
-
  ! Get metric components from metric array
- call unpack_metric(metrici,gammaijUP=gammaijUP,alpha=alpha,betadown=betadown,betaUP=betaUP)
+ call unpack_metric(metrici,gcov=gcov,gammaijUP=gammaijUP,alpha=alpha,betadown=betadown,betaUP=betaUP)
 
+ ! Retrieve sqrt(g)
+ call get_sqrtg(gcov,sqrtg)
+ sqrtg_inv = 1./sqrtg
  pmom2 = 0.
  do i=1,3
     pmom2 = pmom2 + pmom(i)*dot_product(gammaijUP(:,i),pmom(:))
@@ -172,6 +173,7 @@ subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ie
 
  niter = 0
  converged = .false.
+ call get_sqrt_gamma(gcov,sqrt_gamma)
  sqrt_gamma_inv = alpha*sqrtg_inv ! get determinant of 3 spatial metric
  term = rho*sqrt_gamma_inv
  gamfac = gamma/(gamma-1.)
@@ -238,6 +240,7 @@ subroutine conservative2primitive(x,metrici,v,dens,u,P,temp,gamma,rho,pmom,en,ie
 
  if (.not.converged) ierr = 1
 
+
  lorentz_LEO = sqrt(1.+pmom2/enth**2)
  dens = term/lorentz_LEO
 

src/main/deriv.F90

@@ -113,7 +113,7 @@ subroutine derivs(icall,npart,nactive,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,&
     if (gr) then
        ! Recalculate the metric after moving particles to their new tasks
        call init_metric(npart,xyzh,metrics)
-       call prim2consall(npart,xyzh,metrics,vxyzu,dens,pxyzu,use_dens=.false.)
+       !call prim2consall(npart,xyzh,metrics,vxyzu,dens,pxyzu,use_dens=.false.)
     endif
 
     if (nptmass > 0 .and. periodic) call ptmass_boundary_crossing(nptmass,xyzmh_ptmass)

src/main/eos_shen.f90

@@ -249,7 +249,7 @@ end subroutine Cint
 !  Interpolate between values using linear interpolation in 1D
 !+
 !------------------------------------------------------------------------
-subroutine linear_interpolator_one_d(val0,val1,u,val)
+pure subroutine linear_interpolator_one_d(val0,val1,u,val)
  real, intent(out) :: val
  real, intent(in)  :: val0,val1,u