diff --git a/AUTHORS b/AUTHORS index 73d749b5b..d0b207025 100644 --- a/AUTHORS +++ b/AUTHORS @@ -29,6 +29,7 @@ Simone Ceppi MatsEsseldeurs Caitlyn Hardiman Enrico Ragusa +Spencer Magnall fhu Sergei Biriukov Cristiano Longarini @@ -54,11 +55,11 @@ Kieran Hirsh Nicole Rodrigues Amena Faruqi David Trevascus +Farzana Meru Chris Nixon Megha Sharma Nicolas Cuello Benoit Commercon -Farzana Meru Giulia Ballabio Joe Fisher Maxime Lombart @@ -70,6 +71,7 @@ s-neilson <36410751+s-neilson@users.noreply.github.com> Alison Young Cox, Samuel Jorge Cuadra +Miguel Gonzalez-Bolivar Nicolás Cuello Steven Rieder Stéven Toupin diff --git a/build/Makefile b/build/Makefile index 6da744192..1b928d74b 100644 --- a/build/Makefile +++ b/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) diff --git a/build/Makefile_setups b/build/Makefile_setups index 613f126e7..7d41d92be 100644 --- a/build/Makefile_setups +++ b/build/Makefile_setups @@ -1014,6 +1014,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 diff --git a/docs/phantomNR.rst b/docs/phantomNR.rst new file mode 100644 index 000000000..00bc1d662 --- /dev/null +++ b/docs/phantomNR.rst @@ -0,0 +1,79 @@ +PhantomNR +========= + +Using PhantomNR to simulate general relativistic hydrodynamics on dynamical spacetimes +-------------------------------------------------------------------------------------- + +About phantomNR +~~~~~~~~~~~~~~~ + +`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 `__). +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 +------------------------------- + + diff --git a/src/main/config.F90 b/src/main/config.F90 index bb548a994..561adf30e 100644 --- a/src/main/config.F90 +++ b/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 !-------------------- diff --git a/src/main/cons2primsolver.f90 b/src/main/cons2primsolver.f90 index 6d29e0727..ee101a69b 100644 --- a/src/main/cons2primsolver.f90 +++ b/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 diff --git a/src/main/deriv.F90 b/src/main/deriv.F90 index 7fcd13155..3eb0e9f92 100644 --- a/src/main/deriv.F90 +++ b/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) diff --git a/src/main/eos_shen.f90 b/src/main/eos_shen.f90 index 8212addb1..a62ddb51d 100644 --- a/src/main/eos_shen.f90 +++ b/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 diff --git a/src/main/extern_gr.F90 b/src/main/extern_gr.F90 index a3b0ae350..939d7b301 100644 --- a/src/main/extern_gr.F90 +++ b/src/main/extern_gr.F90 @@ -10,7 +10,7 @@ module extern_gr ! ! :References: None ! -! :Owner: David Liptai +! :Owner: Spencer Magnall ! ! :Runtime parameters: None ! @@ -19,7 +19,7 @@ module extern_gr ! implicit none - public :: get_grforce, get_grforce_all, update_grforce_leapfrog + public :: get_grforce, get_grforce_all, update_grforce_leapfrog, get_tmunu_all, get_tmunu_all_exact, get_tmunu private @@ -81,30 +81,28 @@ end subroutine get_grforce_all !--- Subroutine to calculate the timestep constraint from the 'external force' ! this is multiplied by the safety factor C_force elsewhere subroutine dt_grforce(xyzh,fext,dtf) -#ifdef FINVSQRT - use fastmath, only:finvsqrt -#endif use physcon, only:pi - use metric_tools, only:imet_minkowski,imetric + use metric_tools, only:imetric,imet_schwarzschild,imet_kerr real, intent(in) :: xyzh(4),fext(3) real, intent(out) :: dtf real :: r,r2,dtf1,dtf2,f2i integer, parameter :: steps_per_orbit = 100 - + f2i = fext(1)*fext(1) + fext(2)*fext(2) + fext(3)*fext(3) -#ifdef FINVSQRT - dtf1 = sqrt(xyzh(4)*finvsqrt(f2i)) -#else - dtf1 = sqrt(xyzh(4)/sqrt(f2i)) ! This is not really accurate since fi is a component of dp/dt, not da/dt -#endif + if (f2i > 0.) then + dtf1 = sqrt(xyzh(4)/sqrt(f2i)) ! This is not really accurate since fi is a component of dp/dt, not da/dt + else + dtf1 = huge(dtf1) + endif - if (imetric /= imet_minkowski) then + select case (imetric) + case (imet_schwarzschild,imet_kerr) r2 = xyzh(1)*xyzh(1) + xyzh(2)*xyzh(2) + xyzh(3)*xyzh(3) r = sqrt(r2) dtf2 = (2.*pi*sqrt(r*r2))/steps_per_orbit - else + case default dtf2 = huge(dtf2) - endif + end select dtf = min(dtf1,dtf2) @@ -223,4 +221,242 @@ subroutine update_grforce_leapfrog(vhalfx,vhalfy,vhalfz,fxi,fyi,fzi,fexti,dt,xi, end subroutine update_grforce_leapfrog +subroutine get_tmunu_all(npart,xyzh,metrics,vxyzu,metricderivs,dens,tmunus) + use eos, only:ieos,get_pressure + use part, only:isdead_or_accreted + integer, intent(in) :: npart + real, intent(in) :: xyzh(:,:), metrics(:,:,:,:), metricderivs(:,:,:,:), dens(:) + real, intent(inout) :: vxyzu(:,:),tmunus(:,:,:) + real :: pi + integer :: i + logical :: verbose + + verbose = .false. + ! TODO write openmp parallel code + !$omp parallel do default(none) & + !$omp shared(npart,xyzh,metrics,vxyzu,dens,ieos,tmunus) & + !$omp private(i,pi,verbose) + do i=1, npart + !print*, "i: ", i + if (i==1) then + verbose = .true. + else + verbose = .false. + endif + if (.not.isdead_or_accreted(xyzh(4,i))) then + pi = get_pressure(ieos,xyzh(:,i),dens(i),vxyzu(:,i)) + call get_tmunu(xyzh(:,i),metrics(:,:,:,i),& + vxyzu(1:3,i),dens(i),vxyzu(4,i),pi,tmunus(:,:,i),verbose) + endif + enddo + !$omp end parallel do + !print*, "tmunu calc val is: ", tmunus(0,0,5) +end subroutine get_tmunu_all + +subroutine get_tmunu_all_exact(npart,xyzh,metrics,vxyzu,metricderivs,dens,tmunus) + use eos, only:ieos,get_pressure + use part, only:isdead_or_accreted + integer, intent(in) :: npart + real, intent(in) :: xyzh(:,:), metrics(:,:,:,:), metricderivs(:,:,:,:), dens(:) + real, intent(inout) :: vxyzu(:,:),tmunus(:,:,:) + real :: pi + integer :: i + logical :: firstpart + real :: tmunu(4,4) + !print*, "entered get tmunu_all_exact" + tmunu = 0. + firstpart = .true. + ! TODO write openmp parallel code + do i=1, npart + if (.not.isdead_or_accreted(xyzh(4,i)) .and. firstpart) then + pi = get_pressure(ieos,xyzh(:,i),dens(i),vxyzu(:,i)) + call get_tmunu_exact(xyzh(:,i),metrics(:,:,:,i), metricderivs(:,:,:,i), & + vxyzu(1:3,i),dens(i),vxyzu(4,i),pi,tmunus(:,:,i)) + !print*, "finished get_tmunu call!" + firstpart = .false. + !print*, "tmunu: ", tmunu + !print*, "tmunus: ", tmunus(:,:,i) + tmunu(:,:) = tmunus(:,:,i) + !print*, "Got tmunu val: ", tmunu + !stop + else + !print*, "setting tmunu for part: ", i + tmunus(:,:,i) = tmunu(:,:) + endif + + enddo + !print*, "tmunu calc val is: ", tmunus(0,0,5) +end subroutine get_tmunu_all_exact + + +! Subroutine to calculate the covariant form of the stress energy tensor +! For a particle at position p +subroutine get_tmunu(x,metrici,v,dens,u,p,tmunu,verbose) + use metric_tools, only:unpack_metric + use utils_gr, only:get_u0 + real, intent(in) :: x(3),metrici(:,:,:),v(3),dens,u,p + real, intent(out) :: tmunu(0:3,0:3) + logical, optional, intent(in) :: verbose + real :: w,v4(0:3),uzero,u_upper(0:3),u_lower(0:3) + real :: gcov(0:3,0:3), gcon(0:3,0:3) + real :: gammaijdown(1:3,1:3),betadown(3),alpha + integer :: ierr,mu,nu + + ! Reference for all the variables used in this routine: + ! w - the enthalpy + ! gcov - the covariant form of the metric tensor + ! gcon - the contravariant form of the metric tensor + ! gammaijdown - the covariant form of the spatial metric + ! alpha - the lapse + ! betadown - the covariant component of the shift + ! v4 - the uppercase 4 velocity in covariant form + ! v - the fluid velocity v^x + ! vcov - the covariant form of big V_i + ! bigV - the uppercase contravariant V^i + + ! Calculate the enthalpy + w = 1 + u + p/dens + + ! Get cov and con versions of the metric + spatial metric and lapse and shift + ! Not entirely convinced that the lapse and shift calculations are acccurate for the general case!! + !print*, "Before unpack metric " + call unpack_metric(metrici,gcov=gcov,gcon=gcon,gammaijdown=gammaijdown,alpha=alpha,betadown=betadown) + !print*, "After unpack metric" + +! if (present(verbose) .and. verbose) then +! ! Do we get sensible values +! print*, "Unpacked metric quantities..." +! print*, "gcov: ", gcov +! print*, "gcon: ", gcon +! print*, "gammaijdown: ", gammaijdown +! print* , "alpha: ", alpha +! print*, "betadown: ", betadown +! print*, "v4: ", v4 +! endif + + ! ! Need to change Betadown to betaup + ! ! Won't matter at this point as it is allways zero + ! ! get big V + ! bigV(:) = (v(:) + betadown)/alpha + + ! ! We need the covariant version of the 3 velocity + ! ! gamma_ij v^j = v_i where gamma_ij is the spatial metric + ! do i=1, 3 + ! vcov(i) = gammaijdown(i,1)*bigv(1) + gammaijdown(i,2)*bigv(2) + gammaijdown(i,3)*bigv(3) + ! enddo + + + ! ! Calculate the lorentz factor + ! lorentz = (1. - (vcov(1)*bigv(1) + vcov(2)*bigv(2) + vcov(3)*bigv(3)))**(-0.5) + + ! ! Calculate the 4-velocity + ! velshiftterm = vcov(1)*betadown(1) + vcov(2)*betadown(2) + vcov(3)*betadown(3) + ! v4(0) = lorentz*(-alpha + velshiftterm) + ! ! This should be vcov not v + ! v4(1:3) = lorentz*vcov(1:3) + + + ! We are going to use the same Tmunu calc as force GR + ! And then lower it using the metric + ! i.e calc T^{\mu\nu} and then lower it using the metric + ! tensor + ! lower-case 4-velocity (contravariant) + v4(0) = 1. + v4(1:3) = v(:) + + + ! first component of the upper-case 4-velocity (contravariant) + call get_u0(gcov,v,uzero,ierr) + + u_upper = uzero*v4 + do mu=0,3 + u_lower(mu) = gcov(mu,0)*u_upper(0) + gcov(mu,1)*u_upper(1) & + + gcov(mu,2)*u_upper(2) + gcov(mu,3)*u_upper(3) + enddo + + ! Stress energy tensor in contravariant form + do nu=0,3 + do mu=0,3 + tmunu(mu,nu) = w*dens*u_lower(mu)*u_lower(nu) + p*gcov(mu,nu) + enddo + enddo + + +! if (present(verbose) .and. verbose) then +! ! Do we get sensible values +! print*, "Unpacked metric quantities..." +! print*, "gcov: ", gcov +! print*, "gcon: ", gcon +! print*, "gammaijdown: ", gammaijdown +! print* , "alpha: ", alpha +! print*, "betadown: ", betadown +! print*, "v4: ", v4 +! endif + +! if (verbose) then +! print*, "tmunu part: ", tmunu +! print*, "dens: ", dens +! print*, "w: ", w +! print*, "p: ", p +! print*, "gcov: ", gcov +! endif + + ! print*, "tmunu part: ", tmunu + ! print*, "dens: ", dens + ! print*, "w: ", w + ! print*, "p: ", p + ! print*, "gcov: ", gcov + ! stop +end subroutine get_tmunu + +subroutine get_tmunu_exact(x,metrici,metricderivsi,v,dens,u,p,tmunu) + use metric_tools, only:unpack_metric + use utils_gr, only:get_sqrtg + real, intent(in) :: x(3),metrici(:,:,:),metricderivsi(0:3,0:3,3),v(3),dens,u,p + real, intent(out) :: tmunu(0:3,0:3) + real :: w,v4(0:3),vcov(3),lorentz + real :: gcov(0:3,0:3), gcon(0:3,0:3) + real :: gammaijdown(1:3,1:3),betadown(3),alpha + real :: velshiftterm + real :: rhostar,rhoprim,negsqrtg + integer :: i,j + + ! Calculate the enthalpy + ! enthalpy should be 1 as we have zero pressure + ! or should have zero pressure + w = 1 + ! Calculate the exact value of density from conserved density + + call unpack_metric(metrici,gcov=gcov,gcon=gcon,gammaijdown=gammaijdown,alpha=alpha,betadown=betadown) + ! We need the covariant version of the 3 velocity + ! gamma_ij v^j = v_i where gamma_ij is the spatial metric + do i=1, 3 + vcov(i) = gammaijdown(i,1)*v(1) + gammaijdown(i,2)*v(2) + gammaijdown(i,3)*v(3) + enddo + + ! Calculate the lorentz factor + lorentz = (1. - (vcov(1)*v(1) + vcov(2)*v(2) + vcov(3)*v(3)))**(-0.5) + + ! Calculate the 4-velocity + velshiftterm = vcov(1)*betadown(1) + vcov(2)*betadown(2) + vcov(3)*betadown(3) + v4(0) = lorentz*(-alpha + velshiftterm) + v4(1:3) = lorentz*v(1:3) + + rhostar = 13.294563008157013D0 + call get_sqrtg(gcov,negsqrtg) + ! Set/Calculate primitive density using rhostar exactly + rhoprim = rhostar/(negsqrtg/alpha) + + + ! Stress energy tensor + do j=0,3 + do i=0,3 + tmunu(i,j) = rhoprim*w*v4(i)*v4(j) ! + p*gcov(i,j) neglect the pressure term as we don't care + enddo + enddo + + + +end subroutine get_tmunu_exact + end module extern_gr diff --git a/src/main/initial.F90 b/src/main/initial.F90 index a63657cb8..3784e7431 100644 --- a/src/main/initial.F90 +++ b/src/main/initial.F90 @@ -138,8 +138,10 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread) use part, only:metricderivs use cons2prim, only:prim2consall use eos, only:ieos - use extern_gr, only:get_grforce_all + use extern_gr, only:get_grforce_all,get_tmunu_all,get_tmunu_all_exact use metric_tools, only:init_metric,imet_minkowski,imetric + use einsteintk_utils + use tmunu2grid #endif use units, only:utime,umass,unit_Bfield use eos, only:gmw @@ -422,7 +424,6 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread) fxyzu,fext,alphaind,gradh,rad,radprop,dvdx) endif #ifndef PRIM2CONS_FIRST -! COMPUTE METRIC HERE call init_metric(npart,xyzh,metrics,metricderivs) call prim2consall(npart,xyzh,metrics,vxyzu,dens,pxyzu,use_dens=.false.) #endif diff --git a/src/main/interp_metric.F90 b/src/main/interp_metric.F90 new file mode 100644 index 000000000..fc4dd62bf --- /dev/null +++ b/src/main/interp_metric.F90 @@ -0,0 +1,60 @@ +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.bitbucket.io/ ! +!--------------------------------------------------------------------------! +module metric_interp +! +! metric_interp +! +! :References: None +! +! :Owner: Spencer Magnall +! +! :Runtime parameters: None +! +! :Dependencies: einsteintk_utils +! + + interface trilinear_interp + module procedure interp_g, interp_sqrtg, interp_gderiv + end interface trilinear_interp +contains + +subroutine interp_g() +end subroutine interp_g + +subroutine interp_sqrtg() + +end subroutine interp_sqrtg + +subroutine interp_gderiv() + +end subroutine interp_gderiv + +pure subroutine get_grid_neighbours(position,dx,xlower,ylower,zlower) + use einsteintk_utils, only:gridorigin + real, intent(in) :: position(3) + real, intent(in) :: dx(3) + integer, intent(out) :: xlower,ylower,zlower + + ! Get the lower grid neighbours of the position + ! If this is broken change from floor to int + ! How are we handling the edge case of a particle being + ! in exactly the same position as the grid? + ! Hopefully having different grid sizes in each direction + ! Doesn't break the lininterp + xlower = floor((position(1)-gridorigin(1))/dx(1)) + ylower = floor((position(2)-gridorigin(2))/dx(2)) + zlower = floor((position(3)-gridorigin(3))/dx(3)) + + ! +1 because fortran + xlower = xlower + 1 + ylower = ylower + 1 + zlower = zlower + 1 + + +end subroutine get_grid_neighbours + +end module metric_interp diff --git a/src/main/metric_et.f90 b/src/main/metric_et.f90 new file mode 100644 index 000000000..ca792fc92 --- /dev/null +++ b/src/main/metric_et.f90 @@ -0,0 +1,416 @@ +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.bitbucket.io/ ! +!--------------------------------------------------------------------------! +module metric +! +! None +! +! :References: None +! +! :Owner: Spencer Magnall +! +! :Runtime parameters: None +! +! :Dependencies: einsteintk_utils, eos_shen, infile_utils +! + implicit none + character(len=*), parameter :: metric_type = 'et' + integer, parameter :: imetric = 6 + ! This are dummy parameters to stop the compiler complaing + ! 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 +contains + +!---------------------------------------------------------------- +!+ +! Compute the metric tensor in both covariant (gcov) and +! contravariant (gcon) form +!+ +!---------------------------------------------------------------- +pure subroutine get_metric_cartesian(position,gcov,gcon,sqrtg) + use einsteintk_utils, only:gridinit + 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 + + ! The subroutine that computes the metric tensor for a given position + ! In this case it is interpolated from the global grid values + + ! Perform trilenar interpolation + if ( .not. gridinit) then + ! This is required for phantomsetup + ! As no grid information has been passed to phantom from ET + ! So interpolation cannot be performed + gcov = 0. + gcov(0,0) = -1. + gcov(1,1) = 1. + gcov(2,2) = 1. + gcov(3,3) = 1. + if (present(gcon)) then + gcon = 0. + gcon(0,0) = -1. + gcon(1,1) = 1. + gcon(2,2) = 1. + gcon(3,3) = 1. + endif + if (present(sqrtg)) sqrtg = -1. + elseif (present(gcon) .and. present(sqrtg)) then + call interpolate_metric(position,gcov,gcon,sqrtg) + else + call interpolate_metric(position,gcov) + 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 :: r2,sintheta + + gcov = 0. + + r2 = position(1)**2 + sintheta = sin(position(2)) + + gcov(0,0) = -1. + gcov(1,1) = 1. + gcov(2,2) = r2 + gcov(3,3) = r2*sintheta**2 + + if (present(gcon)) then + gcon = 0. + gcon(0,0) = -1. + gcon(1,1) = 1. + gcon(2,2) = 1./r2 + gcov(3,3) = 1./gcov(3,3) + endif + + if (present(sqrtg)) sqrtg = r2*sintheta + +end subroutine get_metric_spherical + +pure subroutine metric_cartesian_derivatives(position,dgcovdx, dgcovdy, dgcovdz) + use einsteintk_utils, only:gridinit + real, intent(in) :: position(3) + real, intent(out) :: dgcovdx(0:3,0:3), dgcovdy(0:3,0:3), dgcovdz(0:3,0:3) + if (.not. gridinit) then + dgcovdx = 0. + dgcovdy = 0. + dgcovdz = 0. + else + call interpolate_metric_derivs(position,dgcovdx,dgcovdy,dgcovdz) + endif +end subroutine metric_cartesian_derivatives + +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 + + r = position(1) + theta = position(2) + + dgcovdr = 0. + dgcovdtheta = 0. + dgcovdphi = 0. + + 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 + +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 + +!----------------------------------------------------------------------- +!+ +! 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)") '# There are no options relating to the '//trim(metric_type)//' metric' + +end subroutine write_options_metric + +!----------------------------------------------------------------------- +!+ +! reads metric options from the input file +!+ +!----------------------------------------------------------------------- +subroutine read_options_metric(name,valstring,imatch,igotall,ierr) + character(len=*), intent(in) :: name,valstring + logical, intent(out) :: imatch,igotall + integer, intent(out) :: ierr + + ! imatch = .true. + ! igotall = .true. + +end subroutine read_options_metric + +!----------------------------------------------------------------------- +!+ +! Interpolates value from grid to position +!+ +!----------------------------------------------------------------------- + +pure subroutine interpolate_metric(position,gcov,gcon,sqrtg) + ! linear and cubic interpolators should be moved to their own subroutine + ! away from eos_shen + use eos_shen, only:linear_interpolator_one_d + use einsteintk_utils, only:gcovgrid,gcongrid,sqrtggrid,dxgrid,gridorigin!,gridsize + real, intent(in) :: position(3) + real, intent(out) :: gcov(0:3,0:3) + real, intent(out), optional :: gcon(0:3,0:3), sqrtg + integer :: xlower,ylower,zlower!,xupper,yupper,zupper + real :: xlowerpos,ylowerpos,zlowerpos + real :: xd,yd,zd + real :: interptmp(7) + integer :: i,j + + ! If the issue is that the metric vals are undefined on + ! Setup since we have not recieved anything about the metric + ! from ET during phantomsetup + ! Then simply set gcov and gcon to 0 + ! as these values will be overwritten during the run anyway + !print*, "Calling interp metric!" + ! Get neighbours + call get_grid_neighbours(position, dxgrid, xlower, ylower, zlower) + !print*,"Neighbours: ", xlower,ylower,zlower + ! This is not true as upper neighbours on the boundary will be on the side + ! take a mod of grid size +! xupper = mod(xlower + 1, gridsize(1)) +! yupper = mod(ylower + 1, gridsize(2)) +! zupper = mod(zlower + 1, gridsize(3)) + ! xupper - xlower should always just be dx provided we are using a uniform grid + ! xd = (position(1) - xlower)/(xupper - xlower) + ! yd = (position(2) - ylower)/(yupper - ylower) + ! zd = (position(3) - zlower)/(zupper - zlower) + xlowerpos = gridorigin(1) + (xlower-1)*dxgrid(1) + ylowerpos = gridorigin(2) + (ylower-1)*dxgrid(2) + zlowerpos = gridorigin(3) + (zlower-1)*dxgrid(3) + + xd = (position(1) - xlowerpos)/(dxgrid(1)) + yd = (position(2) - ylowerpos)/(dxgrid(2)) + zd = (position(3) - zlowerpos)/(dxgrid(3)) + + interptmp = 0. + ! All the interpolation should go into an interface, then you should just call trilinear_interp + ! interpolate for gcov + do i=0, 3 + do j=0, 3 + ! Interpolate along x + call linear_interpolator_one_d(gcovgrid(i,j,xlower,ylower,zlower), & + gcovgrid(i,j,xlower+1,ylower,zlower),xd,interptmp(1)) + call linear_interpolator_one_d(gcovgrid(i,j,xlower,ylower,zlower+1), & + gcovgrid(i,j,xlower+1,ylower,zlower+1),xd,interptmp(2)) + call linear_interpolator_one_d(gcovgrid(i,j,xlower,ylower+1,zlower), & + gcovgrid(i,j,xlower+1,ylower+1,zlower),xd,interptmp(3)) + call linear_interpolator_one_d(gcovgrid(i,j,xlower,ylower+1,zlower+1), & + gcovgrid(i,j,xlower+1,ylower+1,zlower+1),xd,interptmp(4)) + ! Interpolate along y + call linear_interpolator_one_d(interptmp(1),interptmp(3),yd,interptmp(5)) + call linear_interpolator_one_d(interptmp(2),interptmp(4),yd,interptmp(6)) + ! Interpolate along z + call linear_interpolator_one_d(interptmp(5),interptmp(6),zd,interptmp(7)) + + gcov(i,j) = interptmp(7) + enddo + enddo + + if (present(gcon)) then + ! interpolate for gcon + do i=0, 3 + do j=0, 3 + ! Interpolate along x + call linear_interpolator_one_d(gcongrid(i,j,xlower,ylower,zlower), & + gcongrid(i,j,xlower+1,ylower,zlower),xd,interptmp(1)) + call linear_interpolator_one_d(gcongrid(i,j,xlower,ylower,zlower+1), & + gcongrid(i,j,xlower+1,ylower,zlower+1),xd,interptmp(2)) + call linear_interpolator_one_d(gcongrid(i,j,xlower,ylower+1,zlower), & + gcongrid(i,j,xlower+1,ylower+1,zlower),xd,interptmp(3)) + call linear_interpolator_one_d(gcongrid(i,j,xlower,ylower+1,zlower+1), & + gcongrid(i,j,xlower+1,ylower+1,zlower+1),xd,interptmp(4)) + ! Interpolate along y + call linear_interpolator_one_d(interptmp(1),interptmp(3),yd,interptmp(5)) + call linear_interpolator_one_d(interptmp(2),interptmp(4),yd,interptmp(6)) + ! Interpolate along z + call linear_interpolator_one_d(interptmp(5),interptmp(6),zd,interptmp(7)) + + gcon(i,j) = interptmp(7) + enddo + enddo + endif + + if (present(sqrtg)) then + ! Interpolate for sqrtg + ! Interpolate along x + call linear_interpolator_one_d(sqrtggrid(xlower,ylower,zlower), & + sqrtggrid(xlower+1,ylower,zlower),xd,interptmp(1)) + call linear_interpolator_one_d(sqrtggrid(xlower,ylower,zlower+1), & + sqrtggrid(xlower+1,ylower,zlower+1),xd,interptmp(2)) + call linear_interpolator_one_d(sqrtggrid(xlower,ylower+1,zlower), & + sqrtggrid(xlower+1,ylower+1,zlower),xd,interptmp(3)) + call linear_interpolator_one_d(sqrtggrid(xlower,ylower+1,zlower+1), & + sqrtggrid(xlower+1,ylower+1,zlower+1),xd,interptmp(4)) + ! Interpolate along y + call linear_interpolator_one_d(interptmp(1),interptmp(3),yd,interptmp(5)) + call linear_interpolator_one_d(interptmp(2),interptmp(4),yd,interptmp(6)) + ! Interpolate along z + call linear_interpolator_one_d(interptmp(5),interptmp(6),zd,interptmp(7)) + + sqrtg = interptmp(7) + endif + + +end subroutine interpolate_metric + +pure subroutine interpolate_metric_derivs(position,dgcovdx, dgcovdy, dgcovdz) + use eos_shen, only:linear_interpolator_one_d + use einsteintk_utils, only:metricderivsgrid, dxgrid,gridorigin + real, intent(out) :: dgcovdx(0:3,0:3), dgcovdy(0:3,0:3),dgcovdz(0:3,0:3) + real, intent(in) :: position(3) + integer :: xlower,ylower,zlower!,xupper,yupper,zupper + real :: xd,yd,zd,xlowerpos, ylowerpos,zlowerpos + real :: interptmp(7) + integer :: i,j + + call get_grid_neighbours(position, dxgrid, xlower, ylower, zlower) + !print*,"Neighbours: ", xlower,ylower,zlower +! xupper = xlower + 1 +! yupper = yupper + 1 +! zupper = zupper + 1 + ! xd = (position(1) - xlower)/(xupper - xlower) + ! yd = (position(2) - ylower)/(yupper - ylower) + ! zd = (position(3) - zlower)/(zupper - zlower) + + xlowerpos = gridorigin(1) + (xlower-1)*dxgrid(1) + ylowerpos = gridorigin(2) + (ylower-1)*dxgrid(2) + zlowerpos = gridorigin(3) + (zlower-1)*dxgrid(3) + + xd = (position(1) - xlowerpos)/(dxgrid(1)) + yd = (position(2) - ylowerpos)/(dxgrid(2)) + zd = (position(3) - zlowerpos)/(dxgrid(3)) + + interptmp = 0. + + ! Interpolate for dx + do i=0, 3 + do j=0, 3 + ! Interpolate along x + call linear_interpolator_one_d(metricderivsgrid(i,j,1,xlower,ylower,zlower), & + metricderivsgrid(i,j,1,xlower+1,ylower,zlower),xd,interptmp(1)) + call linear_interpolator_one_d(metricderivsgrid(i,j,1,xlower,ylower,zlower+1), & + metricderivsgrid(i,j,1,xlower+1,ylower,zlower+1),xd,interptmp(2)) + call linear_interpolator_one_d(metricderivsgrid(i,j,1,xlower,ylower+1,zlower), & + metricderivsgrid(i,j,1,xlower+1,ylower+1,zlower),xd,interptmp(3)) + call linear_interpolator_one_d(metricderivsgrid(i,j,1,xlower,ylower+1,zlower+1), & + metricderivsgrid(i,j,1,xlower+1,ylower+1,zlower+1),xd,interptmp(4)) + ! Interpolate along y + call linear_interpolator_one_d(interptmp(1),interptmp(3),yd,interptmp(5)) + call linear_interpolator_one_d(interptmp(2),interptmp(4),yd,interptmp(6)) + ! Interpolate along z + call linear_interpolator_one_d(interptmp(5),interptmp(6),zd,interptmp(7)) + + dgcovdx(i,j) = interptmp(7) + enddo + enddo + ! Interpolate for dy + do i=0, 3 + do j=0, 3 + ! Interpolate along x + call linear_interpolator_one_d(metricderivsgrid(i,j,2,xlower,ylower,zlower), & + metricderivsgrid(i,j,2,xlower+1,ylower,zlower),xd,interptmp(1)) + call linear_interpolator_one_d(metricderivsgrid(i,j,2,xlower,ylower,zlower+1), & + metricderivsgrid(i,j,2,xlower+1,ylower,zlower+1),xd,interptmp(2)) + call linear_interpolator_one_d(metricderivsgrid(i,j,2,xlower,ylower+1,zlower), & + metricderivsgrid(i,j,2,xlower+1,ylower+1,zlower),xd,interptmp(3)) + call linear_interpolator_one_d(metricderivsgrid(i,j,2,xlower,ylower+1,zlower+1), & + metricderivsgrid(i,j,2,xlower+1,ylower+1,zlower+1),xd,interptmp(4)) + ! Interpolate along y + call linear_interpolator_one_d(interptmp(1),interptmp(3),yd,interptmp(5)) + call linear_interpolator_one_d(interptmp(2),interptmp(4),yd,interptmp(6)) + ! Interpolate along z + call linear_interpolator_one_d(interptmp(5),interptmp(6),zd,interptmp(7)) + + dgcovdy(i,j) = interptmp(7) + enddo + enddo + + ! Interpolate for dz + do i=0, 3 + do j=0, 3 + ! Interpolate along x + call linear_interpolator_one_d(metricderivsgrid(i,j,3,xlower,ylower,zlower), & + metricderivsgrid(i,j,3,xlower+1,ylower,zlower),xd,interptmp(1)) + call linear_interpolator_one_d(metricderivsgrid(i,j,3,xlower,ylower,zlower+1), & + metricderivsgrid(i,j,3,xlower+1,ylower,zlower+1),xd,interptmp(2)) + call linear_interpolator_one_d(metricderivsgrid(i,j,3,xlower,ylower+1,zlower), & + metricderivsgrid(i,j,3,xlower+1,ylower+1,zlower),xd,interptmp(3)) + call linear_interpolator_one_d(metricderivsgrid(i,j,3,xlower,ylower+1,zlower+1), & + metricderivsgrid(i,j,3,xlower+1,ylower+1,zlower+1),xd,interptmp(4)) + ! Interpolate along y + call linear_interpolator_one_d(interptmp(1),interptmp(3),yd,interptmp(5)) + call linear_interpolator_one_d(interptmp(2),interptmp(4),yd,interptmp(6)) + ! Interpolate along z + call linear_interpolator_one_d(interptmp(5),interptmp(6),zd,interptmp(7)) + + dgcovdz(i,j) = interptmp(7) + enddo + enddo + + + + +end subroutine interpolate_metric_derivs + +pure subroutine get_grid_neighbours(position,dx,xlower,ylower,zlower) + use einsteintk_utils, only:gridorigin + real, intent(in) :: position(3) + real, intent(in) :: dx(3) + integer, intent(out) :: xlower,ylower,zlower + + ! Get the lower grid neighbours of the position + ! If this is broken change from floor to int + ! How are we handling the edge case of a particle being + ! in exactly the same position as the grid? + ! Hopefully having different grid sizes in each direction + ! Doesn't break the lininterp + xlower = floor((position(1)-gridorigin(1))/dx(1)) + ylower = floor((position(2)-gridorigin(2))/dx(2)) + zlower = floor((position(3)-gridorigin(3))/dx(3)) + + ! +1 because fortran + xlower = xlower + 1 + ylower = ylower + 1 + zlower = zlower + 1 + + +end subroutine get_grid_neighbours + + +end module metric diff --git a/src/main/metric_flrw.f90 b/src/main/metric_flrw.f90 new file mode 100644 index 000000000..68152b86d --- /dev/null +++ b/src/main/metric_flrw.f90 @@ -0,0 +1,239 @@ +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.bitbucket.io/ ! +!--------------------------------------------------------------------------! +module metric +! +! None +! +! :References: None +! +! :Owner: Spencer Magnall +! +! :Runtime parameters: None +! +! :Dependencies: infile_utils, timestep +! + + + use timestep, only: time + implicit none + character(len=*), parameter :: metric_type = 'flrw' + integer, parameter :: imetric = 5 + +contains + +!---------------------------------------------------------------- +!+ +! Compute the metric tensor in both covariant (gcov) and +! contravariant (gcon) 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 :: a,t + + t = time + gcov = 0. + ! Get the scale factor for the current time + call get_scale_factor(t,a) + gcov(0,0) = -1. + gcov(1,1) = a + gcov(2,2) = a + gcov(3,3) = a + + if (present(gcon)) then + gcon = 0. + gcon(0,0) = -1. + gcon(1,1) = 1./a + gcon(2,2) = 1./a + gcon(3,3) = 1./a + endif + if (present(sqrtg)) sqrtg = a*a*a + +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 :: r2,sintheta + real :: t,a + + t = time + ! Get the scale factor for the current time + call get_scale_factor(t,a) + + gcov = 0. + + r2 = position(1)**2 + sintheta = sin(position(2)) + + gcov(0,0) = -1. + gcov(1,1) = a + gcov(2,2) = r2*a + gcov(3,3) = a*r2*sintheta**2 + + if (present(gcon)) then + gcon = 0. + gcon(0,0) = -1. + gcon(1,1) = 1./a + gcon(2,2) = 1./(r2*a) + gcov(3,3) = 1./gcov(3,3) + endif + + if (present(sqrtg)) sqrtg = a*a*a + +end subroutine get_metric_spherical + +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) + dgcovdx = 1. + dgcovdy = 1. + dgcovdz = 1. +end subroutine metric_cartesian_derivatives + +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 + real :: t, a + + t = time + ! Get the scale factor for the current time + call get_scale_factor(t,a) + + + r = position(1) + theta = position(2) + + dgcovdr = 0. + dgcovdtheta = 0. + dgcovdphi = 0. + + dgcovdr(2,2) = 2*r*a + dgcovdr(3,3) = 2*r*sin(theta)**2 + + dgcovdtheta(3,3) = 2*a**r**2*cos(theta)*sin(theta) + +end subroutine metric_spherical_derivatives + +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 + +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 + +!----------------------------------------------------------------------- +!+ +! 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)") '# There are no options relating to the '//trim(metric_type)//' metric' + +end subroutine write_options_metric + +!----------------------------------------------------------------------- +!+ +! reads metric options from the input file +!+ +!----------------------------------------------------------------------- +subroutine read_options_metric(name,valstring,imatch,igotall,ierr) + character(len=*), intent(in) :: name,valstring + logical, intent(out) :: imatch,igotall + integer, intent(out) :: ierr + + ! imatch = .true. + ! igotall = .true. + +end subroutine read_options_metric + +pure subroutine get_scale_factor(t,a) + real, intent(in) :: t + real, intent(out) :: a + + a = t*(0.5) + 1 + +end subroutine get_scale_factor + +end module metric diff --git a/src/main/metric_tools.F90 b/src/main/metric_tools.F90 index 17cd33f10..c4acb5c4d 100644 --- a/src/main/metric_tools.F90 +++ b/src/main/metric_tools.F90 @@ -33,7 +33,8 @@ module metric_tools integer, public, parameter :: & imet_minkowski = 1, & ! Minkowski metric imet_schwarzschild = 2, & ! Schwarzschild metric - imet_kerr = 3 ! Kerr metric + imet_kerr = 3, & ! Kerr metric + imet_et = 6 ! Tabulated metric from Einstein toolkit !--- Choice of coordinate system ! (When using this with PHANTOM, it should always be set to cartesian) @@ -187,6 +188,7 @@ subroutine init_metric(npart,xyzh,metrics,metricderivs) real, optional, intent(out) :: metricderivs(:,:,:,:) integer :: i + !$omp parallel do default(none) & !$omp shared(npart,xyzh,metrics) & !$omp private(i) @@ -262,4 +264,6 @@ pure subroutine unpack_metric(metrici,gcov,gcon,gammaijdown,gammaijUP,alpha,beta end subroutine unpack_metric + + end module metric_tools diff --git a/src/main/part.F90 b/src/main/part.F90 index 9a95e47f5..c77b72cc0 100644 --- a/src/main/part.F90 +++ b/src/main/part.F90 @@ -166,6 +166,8 @@ module part real, allocatable :: dens(:) !dens(maxgr) real, allocatable :: metrics(:,:,:,:) !metrics(0:3,0:3,2,maxgr) real, allocatable :: metricderivs(:,:,:,:) !metricderivs(0:3,0:3,3,maxgr) + real, allocatable :: tmunus(:,:,:) !tmunus(0:3,0:3,maxgr) + real, allocatable :: sqrtgs(:) ! sqrtg(maxgr) ! !--sink particles ! @@ -400,6 +402,8 @@ subroutine allocate_part call allocate_array('dens', dens, maxgr) call allocate_array('metrics', metrics, 4, 4, 2, maxgr) call allocate_array('metricderivs', metricderivs, 4, 4, 3, maxgr) + call allocate_array('tmunus', tmunus, 4, 4, maxgr) + call allocate_array('sqrtgs', sqrtgs, maxgr) call allocate_array('xyzmh_ptmass', xyzmh_ptmass, nsinkproperties, maxptmass) call allocate_array('vxyz_ptmass', vxyz_ptmass, 3, maxptmass) call allocate_array('fxyz_ptmass', fxyz_ptmass, 4, maxptmass) @@ -479,6 +483,8 @@ subroutine deallocate_part if (allocated(dens)) deallocate(dens) if (allocated(metrics)) deallocate(metrics) if (allocated(metricderivs)) deallocate(metricderivs) + if (allocated(tmunus)) deallocate(tmunus) + if (allocated(sqrtgs)) deallocate(sqrtgs) if (allocated(xyzmh_ptmass)) deallocate(xyzmh_ptmass) if (allocated(vxyz_ptmass)) deallocate(vxyz_ptmass) if (allocated(fxyz_ptmass)) deallocate(fxyz_ptmass) diff --git a/src/main/readwrite_dumps_fortran.F90 b/src/main/readwrite_dumps_fortran.F90 index 10181a8e4..cbe7212ad 100644 --- a/src/main/readwrite_dumps_fortran.F90 +++ b/src/main/readwrite_dumps_fortran.F90 @@ -218,6 +218,7 @@ subroutine write_fulldump_fortran(t,dumpfile,ntotal,iorder,sphNG) rad,rad_label,radprop,radprop_label,do_radiation,maxirad,maxradprop,itemp,igasP,igamma,& iorig,iX,iZ,imu,nucleation,nucleation_label,n_nucleation,tau,itau_alloc,tau_lucy,itauL_alloc,& luminosity,eta_nimhd,eta_nimhd_label + use part, only:metrics,metricderivs,tmunus use options, only:use_dustfrac,use_var_comp,icooling use dump_utils, only:tag,open_dumpfile_w,allocate_header,& free_header,write_header,write_array,write_block_header @@ -227,6 +228,7 @@ subroutine write_fulldump_fortran(t,dumpfile,ntotal,iorder,sphNG) #ifdef PRDRAG use lumin_nsdisc, only:beta #endif + use metric_tools, only:imetric, imet_et real, intent(in) :: t character(len=*), intent(in) :: dumpfile integer, intent(in), optional :: iorder(:) @@ -363,7 +365,25 @@ subroutine write_fulldump_fortran(t,dumpfile,ntotal,iorder,sphNG) endif if (gr) then call write_array(1,pxyzu,pxyzu_label,maxvxyzu,npart,k,ipass,idump,nums,ierrs(8)) - call write_array(1,dens,'dens prim',npart,k,ipass,idump,nums,ierrs(8)) + call write_array(1,dens,'dens prim',npart,k,ipass,idump,nums,ierrs(8)) + if (imetric==imet_et) then + ! Output metric if imetric=iet + call write_array(1,metrics(1,1,1,:), 'gtt (covariant)',npart,k,ipass,idump,nums,ierrs(8)) + ! call write_array(1,metrics(1,2,1,:), 'gtx (covariant)',npart,k,ipass,idump,nums,ierrs(8)) + ! call write_array(1,metrics(1,3,1,:), 'gty (covariant)',npart,k,ipass,idump,nums,ierrs(8)) + ! call write_array(1,metrics(1,2,1,:), 'gtz (covariant)',npart,k,ipass,idump,nums,ierrs(8)) + ! call write_array(1,metrics(1,2,1,:), 'gtx (covariant)',npart,k,ipass,idump,nums,ierrs(8)) + call write_array(1,metrics(2,2,1,:), 'gxx (covariant)',npart,k,ipass,idump,nums,ierrs(8)) + call write_array(1,metrics(3,3,1,:), 'gyy (covariant)',npart,k,ipass,idump,nums,ierrs(8)) + call write_array(1,metrics(4,4,1,:), 'gzz (covariant)',npart,k,ipass,idump,nums,ierrs(8)) + + call write_array(1,metricderivs(1,1,1,:), 'dxgtt (covariant)',npart,k,ipass,idump,nums,ierrs(8)) + call write_array(1,metricderivs(2,2,1,:), 'dxgxx (covariant)',npart,k,ipass,idump,nums,ierrs(8)) + call write_array(1,metricderivs(3,3,1,:), 'dxgyy (covariant)',npart,k,ipass,idump,nums,ierrs(8)) + call write_array(1,metricderivs(4,4,1,:), 'dxgzz (covariant)',npart,k,ipass,idump,nums,ierrs(8)) + + call write_array(1,tmunus(1,1,:), 'tmunutt (covariant)',npart,k,ipass,idump,nums,ierrs(8)) + endif endif if (eos_is_non_ideal(ieos) .or. (.not.store_dust_temperature .and. icooling > 0)) then call write_array(1,eos_vars(itemp,:),eos_vars_label(itemp),npart,k,ipass,idump,nums,ierrs(12)) diff --git a/src/main/step_leapfrog.F90 b/src/main/step_leapfrog.F90 index 5e4ca24b0..172ff8340 100644 --- a/src/main/step_leapfrog.F90 +++ b/src/main/step_leapfrog.F90 @@ -151,7 +151,6 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) timei = t hdtsph = 0.5*dtsph dterr = bignumber - ! determine twas for each ibin if (ind_timesteps .and. sts_it_n) then time_now = timei + dtsph @@ -359,6 +358,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) endif enddo predict_sph !$omp end parallel do + if (use_dustgrowth) call check_dustprop(npart,dustproppred(1,:)) ! @@ -370,14 +370,17 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) if (npart > 0) then if (gr) vpred = vxyzu ! Need primitive utherm as a guess in cons2prim dt_too_small = .false. + call derivs(1,npart,nactive,xyzh,vpred,fxyzu,fext,divcurlv,& divcurlB,Bpred,dBevol,radpred,drad,radprop,dustproppred,ddustprop,& dustpred,ddustevol,dustfrac,eos_vars,timei,dtsph,dtnew,& ppred,dens,metrics) + if (do_radiation .and. implicit_radiation) then rad = radpred vxyzu(4,1:npart) = vpred(4,1:npart) endif + if (gr) vxyzu = vpred ! May need primitive variables elsewhere? if (dt_too_small) then ! dt < dtmax/2**nbinmax and exit @@ -656,7 +659,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew) endif endif enddo iterations - + ! MPI reduce summary variables nwake = int(reduceall_mpi('+', nwake)) nvfloorp = int(reduceall_mpi('+', nvfloorp)) @@ -804,7 +807,6 @@ subroutine step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,me nsubsteps = 0 dtextforce_min = huge(dt) done = .false. - substeps: do while (timei <= t_end_step .and. .not.done) hdt = 0.5*dt timei = timei + dt @@ -816,7 +818,6 @@ subroutine step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,me if (.not.last_step .and. iverbose > 1 .and. id==master) then write(iprint,"(a,f14.6)") '> external forces only : t=',timei endif - !--------------------------- ! predictor during substeps !--------------------------- @@ -946,7 +947,6 @@ subroutine step_extern_gr(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,pxyzu,dens,me naccreted = 0 nlive = 0 dtextforce_min = bignumber - !$omp parallel default(none) & !$omp shared(npart,xyzh,metrics,metricderivs,vxyzu,fext,iphase,ntypes,massoftype,hdt,timei) & !$omp shared(maxphase,maxp) & diff --git a/src/main/tmunu2grid.f90 b/src/main/tmunu2grid.f90 new file mode 100644 index 000000000..c2ff7ab27 --- /dev/null +++ b/src/main/tmunu2grid.f90 @@ -0,0 +1,455 @@ +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.bitbucket.io/ ! +!--------------------------------------------------------------------------! +module tmunu2grid +! +! tmunu2grid +! +! :References: None +! +! :Owner: Spencer Magnall +! +! :Runtime parameters: None +! +! :Dependencies: boundary, einsteintk_utils, interpolations3D, part +! + implicit none + +contains +subroutine get_tmunugrid_all(npart,xyzh,vxyzu,tmunus,calc_cfac) + use einsteintk_utils, only: dxgrid, gridorigin,gridsize,tmunugrid,rhostargrid + use interpolations3D, only: interpolate3D,interpolate3D_vecexact + use boundary, only: xmin,ymin,zmin,xmax,ymax,zmax + use part, only: massoftype,igas,rhoh + integer, intent(in) :: npart + real, intent(in) :: vxyzu(:,:), tmunus(:,:,:) + real, intent(inout) :: xyzh(:,:) + logical, intent(in), optional :: calc_cfac + real :: weight,h,rho,pmass + real :: weights(npart) + real, save :: cfac + real :: xmininterp(3) + integer :: ngrid(3) + real,allocatable :: datsmooth(:,:,:,:), dat(:,:) + integer :: nnodes,i,k,j, ilower, iupper, jlower, jupper, klower, kupper + logical :: normalise, vertexcen,periodicx,periodicy,periodicz + real :: totalmass + integer :: itype(npart),ilendat + + + ! total mass of the particles + totalmass = npart*massoftype(igas) + + !print*, "totalmass(part): ", totalmass + + ! Density interpolated to the grid + rhostargrid = 0. + if (.not. allocated(datsmooth)) allocate (datsmooth(16,gridsize(1),gridsize(2),gridsize(3))) + if (.not. allocated(dat)) allocate (dat(npart,16)) + ! All particles have equal weighting in the interp + ! Here we calculate the weight for the first particle + ! Get the smoothing length + h = xyzh(4,1) + ! Get pmass + pmass = massoftype(igas) + ! Get density + rho = rhoh(h,pmass) + call get_weight(pmass,h,rho,weight) + ! Correct for Kernel Bias, find correction factor + ! Wrap this into it's own subroutine + if (present(calc_cfac)) then + if (calc_cfac) call get_cfac(cfac,rho) + endif + + weights = weight + itype = 1 + !call get_cfac(cfac,rho) + !print*, "Weighting for particle smoothing is: ", weight + !weight = 1. + ! For now we can set this to the origin, but it might need to be + ! set to the grid origin of the CCTK_grid since we have boundary points + ! TODO This should also be the proper phantom values and not a magic number + !xmin(:) = gridorigin(:) - 0.5*dxgrid(:) ! We move the origin back by 0.5*dx to make a pseudo cell-centered grid + xmininterp(1) = xmin -dxgrid(1) !- 0.5*dxgrid(1) + xmininterp(2) = ymin -dxgrid(2) !- 0.5*dxgrid(2) + xmininterp(3) = zmin-dxgrid(3) !- 0.5*dxgrid(3) + + call get_particle_domain(gridorigin(1),xmin,xmax,dxgrid(1),ilower,iupper) + call get_particle_domain(gridorigin(2),ymin,ymax,dxgrid(2),jlower,jupper) + call get_particle_domain(gridorigin(3),zmin,zmax,dxgrid(3),klower,kupper) + ! nnodes is just the size of the mesh + ! might not be needed + ! We note that this is not actually the size of the einstein toolkit grid + ! As we want our periodic boundary to be on the particle domain not the + ! ET grid domain + ngrid(1) = (iupper-ilower) + 1 + ngrid(2) = (jupper-jlower) + 1 + ngrid(3) = (kupper-klower) + 1 + nnodes = (iupper-ilower)*(jupper-jlower)*(kupper-klower) + ! Do we want to normalise interpolations? + normalise = .true. + ! Is our NR GRID vertex centered? + vertexcen = .false. + periodicx = .true. + periodicy = .true. + periodicz = .true. + + + + ! tt component + + tmunugrid = 0. + datsmooth = 0. + + ! Vectorized tmunu calculation + + ! Put tmunu into an array of form + ! tmunu(npart,16) + do k=1, 4 + do j=1,4 + do i=1,npart + ! Check that this is correct!!! + ! print*,"i j is: ", k, j + ! print*, "Index in array is: ", (k-1)*4 + j + ! print*,tmunus(k,j,1) + dat(i, (k-1)*4 + j) = tmunus(k,j,i) + enddo + enddo +enddo +!stop +ilendat = 16 + +call interpolate3D_vecexact(xyzh,weights,dat,ilendat,itype,npart,& + xmininterp(1),xmininterp(2),xmininterp(3), & + datsmooth(:,ilower:iupper,jlower:jupper,klower:kupper),& + ngrid(1),ngrid(2),ngrid(3),dxgrid(1),dxgrid(2),dxgrid(3),& + normalise,periodicx,periodicy,periodicz) + +! Put the smoothed array into tmunugrid +do i=1,4 + do j=1,4 + ! Check this is correct too! + !print*,"i j is: ", i, j + !print*, "Index in array is: ", (i-1)*4 + j + tmunugrid(i-1,j-1,:,:,:) = datsmooth((i-1)*4 + j, :,:,:) + !print*, "tmunugrid: ", tmunugrid(i-1,j-1,10,10,10) + !print*, datsmooth((i-1)*4 + j, 10,10,10) + enddo +enddo +!stop +! do k=1,4 +! do j=1,4 +! do i=1,4 +! print*, "Lock index is: ", (k-1)*16+ (j-1)*4 + i +! enddo +! enddo +! enddo + +! tmunugrid(0,0,:,:,:) = datsmooth(1,:,:,:) + + ! TODO Unroll this loop for speed + using symmetries + ! Possiblly cleanup the messy indexing +! do k=1,4 +! do j=1,4 +! do i=1, npart +! dat(i) = tmunus(k,j,i) +! enddo + +! ! Get the position of the first grid cell x,y,z +! ! Call to interpolate 3D +! ! COMMENTED OUT AS NOT USED BY NEW INTERPOLATE ROUTINE +! ! call interpolate3D(xyzh,weight,npart, & +! ! xmininterp,tmunugrid(k-1,j-1,ilower:iupper,jlower:jupper,klower:kupper), & +! ! nnodes,dxgrid,normalise,dat,ngrid,vertexcen) + +! !print*, "Interpolated grid values are: ", datsmooth(4:38,4:38,4:38) +! !stop +! ! NEW INTERPOLATION ROUTINE +! call interpolate3D(xyzh,weights,dat,itype,npart,& +! xmininterp(1),xmininterp(2),xmininterp(3), & +! tmunugrid(k-1,j-1,ilower:iupper,jlower:jupper,klower:kupper),& +! ngrid(1),ngrid(2),ngrid(3),dxgrid(1),dxgrid(2),dxgrid(3),& +! normalise,periodicx,periodicy,periodicz) +! enddo +! enddo + + ! RHOSTARGRID CALCULATION IS NOW HANDLED BY AN EXTERNAL ROUTINE + ! THIS IS COMMENTED OUT IN CASE I BREAK EVERYTHING AND NEED TO GO BACK + ! Get the conserved density on the particles + ! dat = 0. + ! do i=1, npart + ! ! Get the smoothing length + ! h = xyzh(4,i) + ! ! Get pmass + ! pmass = massoftype(igas) + ! rho = rhoh(h,pmass) + ! dat(i) = rho + ! enddo + + ! Commented out as not used by new interpolate routine + ! call interpolate3D(xyzh,weight,npart, & + ! xmininterp,rhostargrid(ilower:iupper,jlower:jupper,klower:kupper), & + ! nnodes,dxgrid,.true.,dat,ngrid,vertexcen) + + + ! Calculate the total mass on the grid + !totalmassgrid = 0. + ! do i=ilower,iupper + ! do j=jlower,jupper + ! do k=klower, kupper + ! totalmassgrid = totalmassgrid + dxgrid(1)*dxgrid(2)*dxgrid(3)*rhostargrid(i,j,k) + + ! enddo + ! enddo + ! enddo + ! Explicitly set pressure to be 0 + ! Need to do this in the phantom setup file later + ! tmunugrid(1,0:3,:,:,:) = 0. + ! tmunugrid(2,0:3,:,:,:) = 0. + ! tmunugrid(3,0:3,:,:,:) = 0. + !tmunugrid(0,0,:,:,:) = tmunus(1,1,1) + ! Correction for kernel bias code + ! Hardcoded values for the cubic spline computed using + ! a constant density flrw universe. + ! Ideally this should be in a more general form + ! cfac = totalmass/totalmassgrid + ! ! Output total mass on grid, total mass on particles + ! ! and the residuals + ! !cfac = 0.99917535781746514D0 + ! tmunugrid = tmunugrid*cfac + ! if (iteration==0) then + ! write(666,*) "iteration ", "Mass(Grid) ", "Mass(Particles) ", "Mass(Grid-Particles)" + ! endif + ! write(666,*) iteration, totalmassgrid, totalmass, abs(totalmassgrid-totalmass) + ! close(unit=666) + ! iteration = iteration + 1 + + ! New rho/smoothing length calc based on correction?? + ! not sure that this is a valid thing to do + ! do i=1, npart + ! rho = rhoh(xyzh(i,4),pmass) + ! rho = rho*cfac + ! xyzh(i,4) = hfact*(pmass/rho)**(1./3.) + + ! enddo + + ! Correct rhostargrid using cfac + !rhostargrid = cfac*rhostargrid + + ! Calculate rho(prim), P and e on the grid + ! Apply kernel correction to primatives?? + ! Then calculate a stress energy tensor per grid and fill tmunu + ! A good consistency check would be to do it both ways and compare values + + ! Primative density + + +end subroutine get_tmunugrid_all + +subroutine get_weight(pmass,h,rhoi,weight) + real, intent(in) :: pmass,h,rhoi + real, intent(out) :: weight + + weight = (pmass)/(rhoi*h**3) + +end subroutine get_weight + +subroutine get_dat(tmunus,dat) + real, intent(in) :: tmunus + real, intent(out) :: dat + +end subroutine get_dat + + ! subroutine get_primdens(dens,dat) + ! real, intent(in) :: dens + ! real, intent(out) :: dat + ! integer :: i, npart + + ! ! Get the primative density on the particles + ! dat = 0. + ! do i=1, npart + ! dat(i) = dens(i) + ! enddo + + ! end subroutine get_primdens + + ! subroutine get_4velocity(vxyzu,dat) + ! real, intent(in) :: vxyzu(:,:) + ! real, intent(out) :: dat(:,:) + ! integer :: i,npart + + ! ! Get the primative density on the particles + ! dat = 0. + ! do i=1, npart + ! dat(:,i) = vxyzu(1:3,i) + ! enddo + + ! end subroutine get_4velocity + +subroutine get_particle_domain(gridorigin,xmin,xmax,dxgrid,ilower,iupper) + real, intent(in) :: gridorigin, xmin,xmax, dxgrid + integer, intent(out) :: ilower, iupper + + ! Changed from int to nint + ! to fix a bug + ilower = nint((xmin - gridorigin)/dxgrid) + 1 ! +1 since our arrays start at 1 not 0 + iupper = nint((xmax - gridorigin)/dxgrid) ! Removed the +1 as this was also a bug + ! The lower boundary is in the physical + ! domain but the upper is not; can't have both? +end subroutine get_particle_domain + +subroutine get_cfac(cfac,rho) + real, intent(in) :: rho + real, intent(out) :: cfac + real :: rhoexact + rhoexact = 13.294563008157013D0 + cfac = rhoexact/rho + +end subroutine get_cfac + +subroutine interpolate_to_grid(gridarray,dat) + use einsteintk_utils, only: dxgrid, gridorigin + use interpolations3D, only: interpolate3D + use boundary, only: xmin,ymin,zmin,xmax,ymax,zmax + use part, only:npart,xyzh,massoftype,igas,rhoh + real :: weight,h,rho,pmass + !real, save :: cfac + !integer, save :: iteration = 0 + real :: xmininterp(3) + integer :: ngrid(3) + integer :: nnodes,i, ilower, iupper, jlower, jupper, klower, kupper + logical :: normalise, vertexcen,periodicx, periodicy, periodicz + real, dimension(npart) :: weights + integer, dimension(npart) :: itype + real, intent(out) :: gridarray(:,:,:) ! Grid array to interpolate a quantity to + ! GRID MUST BE RESTRICTED WITH UPPER AND LOWER INDICIES + real, intent(in) :: dat(:) ! The particle data to interpolate to grid + real, allocatable :: interparray(:,:,:) + + + xmininterp(1) = xmin - dxgrid(1)!- 0.5*dxgrid(1) + xmininterp(2) = ymin - dxgrid(2) !- 0.5*dxgrid(2) + xmininterp(3) = zmin - dxgrid(3) !- 0.5*dxgrid(3) + !print*, "xminiterp: ", xmininterp + call get_particle_domain(gridorigin(1),xmin,xmax,dxgrid(1),ilower,iupper) + call get_particle_domain(gridorigin(2),ymin,ymax,dxgrid(2),jlower,jupper) + call get_particle_domain(gridorigin(3),zmin,zmax,dxgrid(3),klower,kupper) + + ! We note that this is not actually the size of the einstein toolkit grid + ! As we want our periodic boundary to be on the particle domain not the + ! ET grid domain + ngrid(1) = (iupper-ilower) + 1 + ngrid(2) = (jupper-jlower) + 1 + ngrid(3) = (kupper-klower) + 1 + allocate(interparray(ngrid(1),ngrid(2),ngrid(3))) + interparray = 0. + nnodes = (iupper-ilower)*(jupper-jlower)*(kupper-klower) + ! Do we want to normalise interpolations? + normalise = .true. + ! Is our NR GRID vertex centered? + vertexcen = .false. + periodicx = .true. + periodicy = .true. + periodicz = .true. + + + + do i=1, npart + h = xyzh(4,i) + ! Get pmass + pmass = massoftype(igas) + ! Get density + rho = rhoh(h,pmass) + call get_weight(pmass,h,rho,weight) + weights(i) = weight + enddo + itype = igas + ! call interpolate3D(xyzh,weight,npart, & + ! xmininterp,gridarray(ilower:iupper,jlower:jupper,klower:kupper), & + ! nnodes,dxgrid,normalise,dat,ngrid,vertexcen) + call interpolate3D(xyzh,weights,dat,itype,npart,& + xmininterp(1),xmininterp(2),xmininterp(3), & + !interparray, & + gridarray(ilower:iupper,jlower:jupper,klower:kupper),& + ngrid(1),ngrid(2),ngrid(3),dxgrid(1),dxgrid(2),dxgrid(3),& + normalise,periodicx,periodicy,periodicz) + + + + +end subroutine interpolate_to_grid + +subroutine check_conserved_dens(rhostargrid,cfac) + use part, only:npart,massoftype,igas + use einsteintk_utils, only: dxgrid, gridorigin + use boundary, only:xmin,xmax,ymin,ymax,zmin,zmax + real, intent(in) :: rhostargrid(:,:,:) + real(kind=16), intent(out) :: cfac + real :: totalmassgrid,totalmasspart + integer :: i,j,k,ilower,iupper,jlower,jupper,klower,kupper + + + call get_particle_domain(gridorigin(1),xmin,xmax,dxgrid(1),ilower,iupper) + call get_particle_domain(gridorigin(2),ymin,ymax,dxgrid(2),jlower,jupper) + call get_particle_domain(gridorigin(3),zmin,zmax,dxgrid(3),klower,kupper) + + totalmassgrid = 0. + do i=ilower,iupper + do j=jlower,jupper + do k=klower, kupper + totalmassgrid = totalmassgrid + dxgrid(1)*dxgrid(2)*dxgrid(3)*rhostargrid(i,j,k) + + enddo + enddo + enddo + + ! total mass of the particles + totalmasspart = npart*massoftype(igas) + + !print*, "Total mass grid: ", totalmassgrid + !print*, "Total mass part: ", totalmasspart + ! Calculate cfac + cfac = totalmasspart/totalmassgrid + + !print*, "cfac mass: ", cfac + +end subroutine check_conserved_dens + +subroutine check_conserved_p(pgrid,cfac) + use part, only:npart,massoftype,igas + use einsteintk_utils, only: dxgrid, gridorigin + use boundary, only:xmin,xmax,ymin,ymax,zmin,zmax + real, intent(in) :: pgrid(:,:,:) + real(kind=16), intent(out) :: cfac + real :: totalmomentumgrid,totalmomentumpart + integer :: i,j,k,ilower,iupper,jlower,jupper,klower,kupper + + call get_particle_domain(gridorigin(1),xmin,xmax,dxgrid(1),ilower,iupper) + call get_particle_domain(gridorigin(2),ymin,ymax,dxgrid(2),jlower,jupper) + call get_particle_domain(gridorigin(3),zmin,zmax,dxgrid(3),klower,kupper) + + ! I'm still a bit unsure what this conserved quantity is actually meant to be?? + totalmomentumgrid = 0. + do i=ilower,iupper + do j=jlower,jupper + do k=klower, kupper + !totalmomentumgrid = totalmomentumgrid + dxgrid(1)*dxgrid(2)*dxgrid(3)*rhostargrid(i,j,k) + + enddo + enddo + enddo + + ! total cons(momentum) of the particles + totalmomentumpart = npart*massoftype(igas) + + ! Calculate cfac + cfac = totalmomentumpart/totalmomentumgrid + + !print*, "cfac mass: ", cfac + +end subroutine check_conserved_p + +end module tmunu2grid diff --git a/src/main/utils_cpuinfo.f90 b/src/main/utils_cpuinfo.f90 index 4fa898b6a..317a6c18b 100644 --- a/src/main/utils_cpuinfo.f90 +++ b/src/main/utils_cpuinfo.f90 @@ -83,6 +83,8 @@ subroutine get_cpuinfo(ncpu,ncpureal,cpuspeed,cpumodel,cachesize,ierr) ncpu = 0 ncpureal = 0 cpuspeed = 0. + cachesizel2 = 0. + cachesizel3 = 0. cpumodel = '' cachesize = '' inquire(file='/proc/cpuinfo',exist=iexist) diff --git a/src/main/utils_gr.F90 b/src/main/utils_gr.F90 index 59e128afd..550b340ec 100644 --- a/src/main/utils_gr.F90 +++ b/src/main/utils_gr.F90 @@ -14,11 +14,12 @@ module utils_gr ! ! :Runtime parameters: None ! -! :Dependencies: fastmath, io, metric_tools, part +! :Dependencies: fastmath, io, metric, metric_tools, part ! implicit none - public :: dot_product_gr, get_u0, get_bigv, rho2dens, h2dens, get_geodesic_accel + public :: dot_product_gr, get_u0, get_bigv, rho2dens, h2dens, get_geodesic_accel, get_sqrtg, get_sqrt_gamma + public :: perturb_metric private @@ -115,9 +116,9 @@ subroutine rho2dens(dens,rho,position,metrici,v) integer :: ierror real :: gcov(0:3,0:3), sqrtg, U0 - ! Hard coded sqrtg=1 since phantom is always in cartesian coordinates - sqrtg = 1. + call unpack_metric(metrici,gcov=gcov) + call get_sqrtg(gcov, sqrtg) call get_u0(gcov,v,U0,ierror) dens = rho/(sqrtg*U0) @@ -156,6 +157,115 @@ subroutine get_geodesic_accel(axyz,npart,vxyz,metrics,metricderivs) end subroutine get_geodesic_accel +subroutine get_sqrtg(gcov, sqrtg) + use metric, only: metric_type + real, intent(in) :: gcov(0:3,0:3) + real, intent(out) :: sqrtg + real :: det + real :: a11,a12,a13,a14 + real :: a21,a22,a23,a24 + real :: a31,a32,a33,a34 + real :: a41,a42,a43,a44 + + + if (metric_type == 'et') then + + a11 = gcov(0,0) + a21 = gcov(1,0) + a31 = gcov(2,0) + a41 = gcov(3,0) + a12 = gcov(0,1) + a22 = gcov(1,1) + a32 = gcov(2,1) + a42 = gcov(3,1) + a13 = gcov(0,2) + a23 = gcov(1,2) + a33 = gcov(2,2) + a43 = gcov(3,2) + a14 = gcov(0,3) + a24 = gcov(1,3) + a34 = gcov(2,3) + a44 = gcov(3,3) + + ! Calculate the determinant + det = a14*a23*a32*a41 - a13*a24*a32*a41 - a14*a22*a33*a41 + a12*a24*a33*a41 + & + a13*a22*a34*a41 - a12*a23*a34*a41 - a14*a23*a31*a42 + a13*a24*a31*a42 + & + a14*a21*a33*a42 - a11*a24*a33*a42 - a13*a21*a34*a42 + a11*a23*a34*a42 + & + a14*a22*a31*a43 - a12*a24*a31*a43 - a14*a21*a32*a43 + a11*a24*a32*a43 + & + a12*a21*a34*a43 - a11*a22*a34*a43 - a13*a22*a31*a44 + a12*a23*a31*a44 + & + a13*a21*a32*a44 - a11*a23*a32*a44 - a12*a21*a33*a44 + a11*a22*a33*a44 + + sqrtg = sqrt(-det) + else + ! If we are not using an evolving metric then + ! Sqrtg = 1 + sqrtg = 1. + endif + + +end subroutine get_sqrtg + +subroutine get_sqrt_gamma(gcov,sqrt_gamma) + use metric, only: metric_type + real, intent(in) :: gcov(0:3,0:3) + real, intent(out) :: sqrt_gamma + real :: a11,a12,a13 + real :: a21,a22,a23 + real :: a31,a32,a33 + !real :: a41,a42,a43 + real :: det + + if (metric_type == 'et') then + ! Calculate the determinant of a 3x3 matrix + ! Spatial metric is just the physical metric + ! without the tt component + + a11 = gcov(1,1) + a12 = gcov(1,2) + a13 = gcov(1,3) + a21 = gcov(2,1) + a22 = gcov(2,2) + a23 = gcov(2,3) + a31 = gcov(3,1) + a32 = gcov(3,2) + a33 = gcov(3,3) + + det = a11*(a22*a33 - a23*a32) - a12*(a21*a33 - a23*a31) + a13*(a21*a32-a22*a31) + sqrt_gamma = sqrt(det) + + else + sqrt_gamma = 1. + + endif + + +end subroutine get_sqrt_gamma + +subroutine perturb_metric(phi,gcovper,gcov) + real, intent(in) :: phi + real, intent(out) :: gcovper(0:3,0:3) + real, optional, intent(in) :: gcov(0:3,0:3) + + + if (present(gcov)) then + gcovper = gcov + else + gcovper = 0. + gcovper(0,0) = -1. + gcovper(1,1) = 1. + gcovper(2,2) = 1. + gcovper(3,3) = 1. + endif + + ! Set the pertubed metric based on the Bardeen formulation + gcovper(0,0) = gcovper(0,0) - 2.*phi + gcovper(1,1) = gcovper(1,1) - 2.*phi + gcovper(2,2) = gcovper(2,2) - 2.*phi + gcovper(3,3) = gcovper(3,3) - 2.*phi + + +end subroutine perturb_metric + ! This is not being used at the moment. ! subroutine dens2rho(rho,dens,position,v) ! use metric_tools, only: get_metric diff --git a/src/main/utils_infiles.f90 b/src/main/utils_infiles.f90 index ec49fd6f4..2d47ad151 100644 --- a/src/main/utils_infiles.f90 +++ b/src/main/utils_infiles.f90 @@ -44,7 +44,7 @@ module infile_utils ! maximum length for input strings ! (if you change this, must also change format statements below) ! - integer, parameter, private :: maxlen = 20 ! max length of string containing variable + integer, parameter, private :: maxlen = 100 ! max length of string containing variable integer, parameter, private :: maxlenval = 100 ! max length of string containing value integer, parameter, private :: maxlenstring = 120 ! max length of string variable integer, parameter, private :: maxlenline = 120 ! maximum line length @@ -177,6 +177,7 @@ subroutine write_inopt_real8(rval,name,descript,iunit,ierr,exp,time) logical :: doexp,dotime integer :: nhr,nmin !,nsec character(len=16) :: tmpstring + character(len=3) :: fmts real(kind=8) :: trem integer :: ierror @@ -189,6 +190,9 @@ subroutine write_inopt_real8(rval,name,descript,iunit,ierr,exp,time) if (time) dotime = .true. endif + fmts = "a20" + if (len_trim(name) > 20) fmts = "a" + if (dotime) then trem = rval nhr = int(trem/3600.d0) @@ -197,12 +201,12 @@ subroutine write_inopt_real8(rval,name,descript,iunit,ierr,exp,time) if (nmin > 0) trem = trem - nmin*60.d0 !nsec = int(trem) - write(iunit,"(a20,' = ',5x,i3.3,':',i2.2,4x,'! ',a)",iostat=ierror) & + write(iunit,"("//trim(fmts)//",' = ',5x,i3.3,':',i2.2,4x,'! ',a)",iostat=ierror) & name,nhr,nmin,descript else if (doexp .or. (abs(rval) < 1.e-3 .and. abs(rval) > tiny(rval)) & .or. (abs(rval) >= 1.e4)) then - write(iunit,"(a20,' = ',1x,es10.3,4x,'! ',a)",iostat=ierror) & + write(iunit,"("//trim(fmts)//",' = ',1x,es10.3,4x,'! ',a)",iostat=ierror) & name,rval,descript else if (abs(rval) <= 1.e-1) then @@ -215,10 +219,11 @@ subroutine write_inopt_real8(rval,name,descript,iunit,ierr,exp,time) write(tmpstring,"(g16.9)",iostat=ierror) rval tmpstring = adjustl(strip_zeros(tmpstring,3)) endif + if (len_trim(tmpstring) > 10) then - write(iunit,"(a20,' = ',1x,a,2x,'! ',a)",iostat=ierror) name,adjustr(trim(tmpstring)),descript + write(iunit,"("//trim(fmts)//",' = ',1x,a,2x,'! ',a)",iostat=ierror) name,adjustr(trim(tmpstring)),descript else - write(iunit,"(a20,' = ',1x,a10,4x,'! ',a)",iostat=ierror) name,adjustr(trim(tmpstring)),descript + write(iunit,"("//trim(fmts)//",' = ',1x,a10,4x,'! ',a)",iostat=ierror) name,adjustr(trim(tmpstring)),descript endif endif endif @@ -268,12 +273,16 @@ subroutine write_inopt_string(sval,name,descript,iunit,ierr) integer, intent(in) :: iunit integer, intent(out), optional :: ierr character(len=40) :: fmtstring + character(len=3) :: fmts integer :: ierror + fmts = "a20" + if (len_trim(name) > 20) fmts = "a" + if (len_trim(sval) > 10) then - fmtstring = '(a20,'' = '',1x,a,3x,''! '',a)' + fmtstring = '('//fmts//','' = '',1x,a,3x,''! '',a)' else - fmtstring = '(a20,'' = '',1x,a10,4x,''! '',a)' + fmtstring = '('//fmts//','' = '',1x,a10,4x,''! '',a)' endif write(iunit,fmtstring,iostat=ierror) name,trim(sval),trim(descript) @@ -517,7 +526,6 @@ subroutine read_inopt_string(valstring,tag,db,err,errcount) ierr = 0 if (.not.match_inopt_in_db(db,tag,valstring)) ierr = -1 - if (present(err)) then err = ierr elseif (ierr /= 0) then diff --git a/src/setup/phantomsetup.F90 b/src/setup/phantomsetup.F90 index 5d0a13cb9..3c44b2de5 100644 --- a/src/setup/phantomsetup.F90 +++ b/src/setup/phantomsetup.F90 @@ -119,7 +119,6 @@ program phantomsetup call init_domains(nprocs) id = 0 endif - do myid=0,nprocsfake-1 myid1 = myid diff --git a/src/setup/set_star.f90 b/src/setup/set_star.f90 index 24071d154..66a5d476c 100644 --- a/src/setup/set_star.f90 +++ b/src/setup/set_star.f90 @@ -15,7 +15,19 @@ module setstar ! ! :Owner: Daniel Price ! -! :Runtime parameters: None +! :Runtime parameters: +! - Mstar : *mass of star* +! - Rstar : *radius of star* +! - hsoft : *Softening length of sink particle stellar core* +! - input_profile : *Path to input profile* +! - isinkcore : *Add a sink particle stellar core* +! - isoftcore : *0=no core softening, 1=cubic core, 2=constant entropy core* +! - isofteningopt : *1=supply rcore, 2=supply mcore, 3=supply both* +! - mcore : *Mass of sink particle stellar core* +! - np : *number of particles* +! - outputfilename : *Output path for softened MESA profile* +! - rcore : *Radius of core softening* +! - ui_coef : *specific internal energy (units of GM/R)* ! ! :Dependencies: centreofmass, dim, eos, extern_densprofile, infile_utils, ! io, mpiutils, part, physcon, prompting, radiation_utils, relaxstar, diff --git a/src/setup/set_unifdis.f90 b/src/setup/set_unifdis.f90 index 64ccd8bc5..e85096dc6 100644 --- a/src/setup/set_unifdis.f90 +++ b/src/setup/set_unifdis.f90 @@ -16,7 +16,7 @@ module unifdis ! ! :Dependencies: random, stretchmap ! - use stretchmap, only:rho_func + use stretchmap, only:rho_func, mass_func implicit none public :: set_unifdis, get_ny_nz_closepacked, get_xyzmin_xyzmax_exact public :: is_valid_lattice, is_closepacked, latticetype @@ -53,7 +53,7 @@ end function mask_prototype subroutine set_unifdis(lattice,id,master,xmin,xmax,ymin,ymax, & zmin,zmax,delta,hfact,np,xyzh,periodic, & rmin,rmax,rcylmin,rcylmax,rellipsoid,in_ellipsoid, & - nptot,npy,npz,npnew_in,rhofunc,inputiseed,verbose,centre,dir,geom,mask,err) + nptot,npy,npz,npnew_in,rhofunc,massfunc,inputiseed,verbose,centre,dir,geom,mask,err) use random, only:ran2 use stretchmap, only:set_density_profile !use mpidomain, only:i_belong @@ -70,6 +70,7 @@ subroutine set_unifdis(lattice,id,master,xmin,xmax,ymin,ymax, & integer(kind=8), intent(inout), optional :: nptot integer, intent(in), optional :: npy,npz,npnew_in,dir,geom procedure(rho_func), pointer, optional :: rhofunc + procedure(mass_func), pointer, optional :: massfunc integer, intent(in), optional :: inputiseed logical, intent(in), optional :: verbose,centre,in_ellipsoid integer, intent(out), optional :: err @@ -587,7 +588,7 @@ subroutine set_unifdis(lattice,id,master,xmin,xmax,ymin,ymax, & endif endif call set_density_profile(np,xyzh,min=xmins,max=xmaxs,rhofunc=rhofunc,& - start=npin,geom=igeom,coord=icoord,verbose=(id==master .and. is_verbose),err=ierr) + start=npin,geom=igeom,coord=icoord,verbose=(id==master .and. is_verbose),err=ierr)!,massfunc=massfunc) if (ierr > 0) then if (present(err)) err = ierr return diff --git a/src/setup/setup_flrw.f90 b/src/setup/setup_flrw.f90 new file mode 100644 index 000000000..7cdc8c868 --- /dev/null +++ b/src/setup/setup_flrw.f90 @@ -0,0 +1,604 @@ +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.bitbucket.io/ ! +!--------------------------------------------------------------------------! +module setup +! +! Setup routine for a constant density + petrubtations FLRW universe +! as described in Magnall et al. 2023 +! +! :References: None +! +! :Owner: Spencer Magnall +! +! :Runtime parameters: +! - cs0 : *initial sound speed in code units* +! - dist_unit : *distance unit (e.g. au)* +! - ilattice : *lattice type (1=cubic, 2=closepacked)* +! - mass_unit : *mass unit (e.g. solarm)* +! - nx : *number of particles in x direction* +! - radiation_dominated : *Radiation dominated universe (yes/no)* +! - rhozero : *initial density in code units* +! +! :Dependencies: boundary, dim, infile_utils, io, mpidomain, mpiutils, +! options, part, physcon, prompting, setup_params, stretchmap, unifdis, +! units, utils_gr +! + use dim, only:use_dust + use setup_params, only:rhozero + use physcon, only:radconst + implicit none + public :: setpart + + integer :: npartx,ilattice + real :: cs0,xmini,xmaxi,ymini,ymaxi,zmini,zmaxi,ampl,phaseoffset + character(len=20) :: dist_unit,mass_unit,perturb_direction,perturb,radiation_dominated + real :: perturb_wavelength + real :: rho_matter + real(kind=8) :: udist,umass + + private + +contains + +!---------------------------------------------------------------- +!+ +! setup for uniform particle distributions +!+ +!---------------------------------------------------------------- +subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,time,fileprefix) + use dim, only:maxvxyzu,gr + use setup_params, only:npart_total + use io, only:master + use unifdis, only:set_unifdis,rho_func!,mass_func + use boundary, only:xmin,ymin,zmin,xmax,ymax,zmax,dxbound,dybound,dzbound,set_boundary + use part, only:periodic + use physcon, only:years,pc,solarm,pi + use units, only:set_units + use mpidomain, only:i_belong + use stretchmap, only:set_density_profile + use utils_gr, only:perturb_metric, get_u0, get_sqrtg + !use cons2primsolver, only:primative2conservative + + integer, intent(in) :: id + integer, intent(inout) :: npart + integer, intent(out) :: npartoftype(:) + real, intent(out) :: xyzh(:,:) + real, intent(out) :: massoftype(:) + real, intent(out) :: polyk,gamma + real, intent(inout) :: hfact + real, intent(inout) :: time + character(len=20), intent(in) :: fileprefix + real, intent(out) :: vxyzu(:,:) + character(len=40) :: filename,lattice + real :: totmass,deltax + integer :: i,ierr + logical :: iexist + real :: kwave,denom,length, c1,c3,lambda + real :: xval + real :: Vup(0:3),phi,sqrtg,gcov(0:3,0:3),alpha,hub + real :: last_scattering_temp + procedure(rho_func), pointer :: density_func + + density_func => rhofunc ! desired density function + + + ! + !--general parameters + ! + perturb_wavelength = 1.0 + time = 0. + if (maxvxyzu < 4) then + gamma = 1. + else + ! 4/3 for radiation dominated case + ! irrelevant for + gamma = 4./3. + endif + + ! + ! default units + ! + mass_unit = 'solarm' + dist_unit = 'mpc' + ! + ! set boundaries to default values + ! + xmini = xmin; xmaxi = xmax + ymini = ymin; ymaxi = ymax + zmini = zmin; zmaxi = zmax + ! + ! set default values for input parameters + ! + npartx = 64 + ilattice = 1 + perturb = '"no"' + perturb_direction = '"none"' + radiation_dominated = '"no"' + + ! Ideally this should read the values of the box length + ! and initial Hubble parameter from the par file. + ! Then it should be set using the Friedmann equation: + !!!!!! rhozero = (3H^2)/(8*pi*a*a) + + hub = 10.553495658357338!/10. + !hub = 23.588901903912664 + !hub = 0.06472086375185665 + rhozero = 3. * hub**2 / (8. * pi) + phaseoffset = 0. + ampl = 0. + + ! Approx Temp of the CMB in Kelvins + !last_scattering_temp = 3000 + !last_scattering_temp = (rhozero/radconst)**(1./4.)*0.999999999999999d0 + last_scattering_temp = 0. + + ! Define some parameters for Linear pertubations + ! We assume ainit = 1, but this may not always be the case + c1 = 1./(4.*pi*rhozero) + !c2 = We set g(x^i) = 0 as we only want to extract the growing mode + c3 = - sqrt(1./(6.*pi*rhozero)) + !c3 = hub/(4.d0*PI*rhozero) + + + if (gr) then + ! 0 Because dust? + cs0 = 0. + else + cs0 = 1. + endif + + ! get disc setup parameters from file or interactive setup + ! + filename=trim(fileprefix)//'.setup' + inquire(file=filename,exist=iexist) + if (iexist) then + !--read from setup file + call read_setupfile(filename,ierr) + if (id==master) call write_setupfile(filename) + if (ierr /= 0) then + stop + endif + elseif (id==master) then + call setup_interactive(id,polyk) + call write_setupfile(filename) + stop 'rerun phantomsetup after editing .setup file' + else + stop + endif + ! + ! set units and boundaries + ! + if (gr) then + call set_units(dist=udist,c=1.d0,G=1.d0) + else + call set_units(dist=udist,mass=umass,G=1.d0) + endif + call set_boundary(xmini,xmaxi,ymini,ymaxi,zmini,zmaxi) + ! + ! setup particles + ! + + npart = 0 + npart_total = 0 + length = xmaxi - xmini + deltax = length/npartx +! +! general parameters +! +! time should be read in from the par file + !time = 0.08478563386065302 + time = 0.18951066686763596 ! z~1000 + lambda = perturb_wavelength*length + kwave = (2.d0*pi)/lambda + denom = length - ampl/kwave*(cos(kwave*length)-1.0) + rhozero = 3.d0 * hub**2 / (8.d0 * pi) + print*, rhozero + + select case(radiation_dominated) + case('"yes"') + + ! Set a value of rho_matter + rho_matter = 1.e-40 + !rhozero = rhozero - radconst*last_scattering_temp**4 + ! Solve for temperature + last_scattering_temp = ((rhozero-rho_matter)/radconst)**(1./4.) + rhozero = rho_matter + end select + + xval = density_func(0.75) + xval = density_func(0.5) + + select case(ilattice) + case(2) + lattice = 'closepacked' + case default + if (ilattice /= 1) print*,' error: chosen lattice not available, using cubic' + lattice = 'cubic' + end select + + select case(perturb) + case('"yes"') + select case(perturb_direction) + !TODO Z AND Y LINEAR PERTURBATIONS + case('"x"') + call set_unifdis(lattice,id,master,xmin,xmax,ymin,ymax,zmin,zmax,deltax,hfact,& + npart,xyzh,periodic,nptot=npart_total,mask=i_belong,rhofunc=density_func)!,massfunc=mass_function) + case('"y"') + call set_unifdis(lattice,id,master,xmin,xmax,ymin,ymax,zmin,zmax,deltax,hfact,& + npart,xyzh,periodic,nptot=npart_total,mask=i_belong) + call set_density_profile(npart,xyzh,min=ymin,max=ymax,rhofunc=density_func,& + geom=1,coord=2) + case('"all"') + call set_unifdis(lattice,id,master,xmin,xmax,ymin,ymax,zmin,zmax,deltax,hfact,& + npart,xyzh,periodic,nptot=npart_total,mask=i_belong,rhofunc=density_func) + call set_density_profile(npart,xyzh,min=ymin,max=ymax,rhofunc=density_func,& + geom=1,coord=2) + call set_density_profile(npart,xyzh,min=zmin,max=zmax,rhofunc=density_func,& + geom=1,coord=3) + end select + case('"no"') + call set_unifdis(lattice,id,master,xmin,xmax,ymin,ymax,zmin,zmax,deltax,hfact,& + npart,xyzh,periodic,nptot=npart_total,mask=i_belong) + end select + + npartoftype(:) = 0 + npartoftype(1) = npart + print*,' npart = ',npart,npart_total + + + totmass = rhozero*dxbound*dybound*dzbound + massoftype(1) = totmass/npart_total + if (id==master) print*,' particle mass = ',massoftype(1) + if (id==master) print*,' initial sound speed = ',cs0,' pressure = ',cs0**2/gamma + + + + if (maxvxyzu < 4 .or. gamma <= 1.) then + polyk = cs0**2 + else + polyk = 0. + endif + do i=1,npart + + select case(perturb_direction) + case ('"x"') + ! should not be zero, for a perturbed wave + vxyzu(1,i) = kwave*c3*ampl*cos((2.*pi*xyzh(1,i))/lambda - phaseoffset) + phi = ampl*sin(kwave*xyzh(1,i)-phaseoffset) + Vup(1) = kwave*c3*ampl*cos(2.*pi*xyzh(1,i) - phaseoffset) + Vup(2:3) = 0. + call perturb_metric(phi,gcov) + call get_sqrtg(gcov,sqrtg) + + alpha = sqrt(-gcov(0,0)) + vxyzu(1,i) = Vup(1)*alpha + vxyzu(2:3,i) = 0. + case ('"y"') + vxyzu(2,i) = kwave*c3*ampl*cos((2.*pi*xyzh(2,i))/lambda - phaseoffset) + phi = ampl*sin(kwave*xyzh(2,i)-phaseoffset) + Vup = 0. + Vup(2) = kwave*c3*ampl*cos(2.*pi*xyzh(2,i) - phaseoffset) + + call perturb_metric(phi,gcov) + call get_sqrtg(gcov,sqrtg) + + alpha = sqrt(-gcov(0,0)) + vxyzu(:,i) = 0. + vxyzu(2,i) = Vup(2)*alpha + + case ('"all"') + phi = ampl*(sin(kwave*xyzh(1,i)-phaseoffset) - sin(kwave*xyzh(2,i)-phaseoffset) - sin(kwave*xyzh(3,i)-phaseoffset)) + Vup(1) = kwave*c3*ampl*cos((2.*pi*xyzh(1,i))/lambda - phaseoffset) + Vup(2) = kwave*c3*ampl*cos((2.*pi*xyzh(2,i))/lambda - phaseoffset) + Vup(3) = kwave*c3*ampl*cos((2.*pi*xyzh(3,i))/lambda - phaseoffset) + + call perturb_metric(phi,gcov) + call get_sqrtg(gcov,sqrtg) + + alpha = sqrt(-gcov(0,0)) + + ! perturb the y and z velocities + vxyzu(1,i) = Vup(1)*alpha + vxyzu(2,i) = Vup(2)*alpha + vxyzu(3,i) = Vup(3)*alpha + end select + ! Setup the intial internal energy here? + ! This should be u = aT^4/\rho + ! Choose an initial temp of the cmb ~ 3000K + ! Set a=1 for now + ! Asssuming that this is constant density/pressure for now so I'm making sure that + ! Note that rhozero != rho + ! rhozero = rho + rho*u as this is the energy density + select case(radiation_dominated) + case('"yes"') + if (maxvxyzu >= 4 .and. gamma > 1.) vxyzu(4,i) = (radconst*(last_scattering_temp**4))/rhozero !vxyzu(4,i) = cs0**2/(gamma*(gamma-1.)) + ! Check that the pressure is correct + print*, "Pressure: ", (gamma-1)*rhozero*vxyzu(4,i) + print*, "Pressure from energy density: ", 3. * hub**2 / (8. * pi)/3. + print*, "Pressure 1/3 \rho u: ",radconst*(last_scattering_temp**4)/3. + print*, "particle mass: ", massoftype + end select + enddo + + +contains +!---------------------------------------------------- +!+ +! callback function giving desired density profile +!+ +!---------------------------------------------------- +real function rhofunc(x) + use utils_gr, only:perturb_metric, get_u0, get_sqrtg + real, intent(in) :: x + real :: const, phi, rhoprim, gcov(0:3,0:3), sqrtg,u0,v(3),Vup(3) + real :: alpha + integer :: ierr + + !rhofunc = 1.d0 + ampl*sin(kwave*(x-xmin)) + !rhofunc = ampl*sin(kwave*(x-xmin)) + ! Eq 28. in Macpherson+ 2017 + ! Although it is missing a negative sign + const = -kwave*kwave*c1 - 2. + phi = ampl*sin(kwave*x-phaseoffset) + !rhofunc = rhozero*(1.d0 + const*ampl*sin(kwave*x)) + ! Get the primative density from the linear perb + rhoprim = rhozero*(1.d0+const*phi) + + ! Get the perturbed 4-metric + call perturb_metric(phi,gcov) + ! Get sqrt(-det(g)) + call get_sqrtg(gcov,sqrtg) + ! Define the 3 velocities to calculate u0 + ! Three velocity will need to be converted from big V to small v + ! + Vup(1) = kwave*c3*ampl*cos((2.*pi*x)/lambda-phaseoffset) + Vup(2:3) = 0. + alpha = sqrt(-gcov(0,0)) + v(1) = Vup(1)*alpha + v(2:3) = 0. + ! calculate u0 + ! TODO Should probably handle this error at some point + call get_u0(gcov,v,u0,ierr) + ! Perform a prim2cons + rhofunc = rhoprim*u0*sqrtg + +end function rhofunc + +real function massfunc(x,xmin) + use utils_gr, only:perturb_metric, get_u0, get_sqrtg,dot_product_gr + real, intent(in) :: x,xmin + real :: const, gcov(0:3,0:3), sqrtg,u0,v(3),Vup(3) + real :: massprimx,massprimmin,massprim + real :: lorrentz, bigv2 + + ! The value inside the bracket + const = -kwave*kwave*c1 - 2. + phi = ampl*sin(kwave*x-phaseoffset) + !expr = ampl*(-(1./kwave))*cos(phaseoffset - (2.d0*pi*x)/lambda) + !exprmin = ampl*(-(1./kwave))*cos(phaseoffset - (2.d0*pi*xmin)/lambda) + massprimx = (x+deltaint(x)) + massprimmin = (xmin+deltaint(xmin)) + ! Evalutation of the integral + ! rho0[x-Acos(kx)]^x_0 + massprim = rhozero*(massprimx - massprimmin) + print*, massprim + + ! Get the perturbed 4-metric + call perturb_metric(phi,gcov) + ! Get sqrt(-det(g)) + call get_sqrtg(gcov,sqrtg) + ! Define the 3 velocities to calculate u0 + ! Three velocity will need to be converted from big V to small v + ! + Vup(1) = kwave*c3*ampl*cos((2.*pi*x)/lambda-phaseoffset) + Vup(2:3) = 0. + alpha = sqrt(-gcov(0,0)) + !v(0) = 1 + v(1) = Vup(1)*alpha + v(2:3) = 0. + bigv2 = dot_product_gr(Vup,Vup,gcov) + lorrentz = 1./sqrt(1.-bigv2) + call get_u0(gcov,v,u0,ierr) + massfunc = (massprim)!*lorrentz + massfunc = massprim!*sqrtg*u0 + +end function massfunc + +real function deltaint(x) + real, intent(in) :: x + + deltaint = (1./kwave)*(kwave*kwave*c1 - 2)*ampl*cos(2*pi*x/lambda) + +end function deltaint + +end subroutine setpart + +!------------------------------------------------------------------------ +! +! interactive setup +! +!------------------------------------------------------------------------ +subroutine setup_interactive(id,polyk) + use io, only:master + use mpiutils, only:bcast_mpi + use dim, only:maxp,maxvxyzu + use prompting, only:prompt + use units, only:select_unit + integer, intent(in) :: id + real, intent(out) :: polyk + integer :: ierr + + if (id==master) then + ierr = 1 + do while (ierr /= 0) + call prompt('Enter mass unit (e.g. solarm,jupiterm,earthm)',mass_unit) + call select_unit(mass_unit,umass,ierr) + if (ierr /= 0) print "(a)",' ERROR: mass unit not recognised' + enddo + ierr = 1 + do while (ierr /= 0) + call prompt('Enter distance unit (e.g. au,pc,kpc,0.1pc)',dist_unit) + call select_unit(dist_unit,udist,ierr) + if (ierr /= 0) print "(a)",' ERROR: length unit not recognised' + enddo + + call prompt('enter xmin boundary',xmini) + call prompt('enter xmax boundary',xmaxi,xmini) + call prompt('enter ymin boundary',ymini) + call prompt('enter ymax boundary',ymaxi,ymini) + call prompt('enter zmin boundary',zmini) + call prompt('enter zmax boundary',zmaxi,zmini) + endif + ! + ! number of particles + ! + if (id==master) then + print*,' uniform setup... (max = ',nint((maxp)**(1/3.)),')' + call prompt('enter number of particles in x direction ',npartx,1) + endif + call bcast_mpi(npartx) + ! + ! mean density + ! + if (id==master) call prompt(' enter density (gives particle mass)',rhozero,0.) + call bcast_mpi(rhozero) + ! + ! sound speed in code units + ! + if (id==master) then + call prompt(' enter sound speed in code units (sets polyk)',cs0,0.) + endif + call bcast_mpi(cs0) + ! + ! type of lattice + ! + if (id==master) then + call prompt(' select lattice type (1=cubic, 2=closepacked)',ilattice,1) + endif + call bcast_mpi(ilattice) +end subroutine setup_interactive + +!------------------------------------------------------------------------ +! +! write setup file +! +!------------------------------------------------------------------------ +subroutine write_setupfile(filename) + use infile_utils, only:write_inopt + character(len=*), intent(in) :: filename + integer :: iunit + + print "(/,a)",' writing setup options file '//trim(filename) + open(newunit=iunit,file=filename,status='replace',form='formatted') + write(iunit,"(a)") '# input file for uniform setup routine' + + write(iunit,"(/,a)") '# units' + call write_inopt(dist_unit,'dist_unit','distance unit (e.g. au)',iunit) + call write_inopt(mass_unit,'mass_unit','mass unit (e.g. solarm)',iunit) + ! + ! boundaries + ! + write(iunit,"(/,a)") '# boundaries' + call write_inopt(xmini,'CoordBase::xmin','xmin boundary',iunit) + call write_inopt(xmaxi,'CoordBase::xmax','xmax boundary',iunit) + call write_inopt(ymini,'CoordBase::ymin','ymin boundary',iunit) + call write_inopt(ymaxi,'CoordBase::ymax','ymax boundary',iunit) + call write_inopt(zmini,'CoordBase::zmin','zmin boundary',iunit) + call write_inopt(zmaxi,'CoordBase::zmax','zmax boundary',iunit) + + + + ! + ! other parameters + ! + write(iunit,"(/,a)") '# setup' + call write_inopt(npartx,'nx','number of particles in x direction',iunit) + call write_inopt(rhozero,'rhozero','initial density in code units',iunit) + call write_inopt(cs0,'cs0','initial sound speed in code units',iunit) + call write_inopt(perturb,'FLRWSolver::FLRW_perturb','Pertrubations of FLRW?',iunit) + call write_inopt(ampl,'FLRWSolver::phi_amplitude','Pertubation amplitude',iunit) + call write_inopt(phaseoffset,'FLRWSolver::phi_phase_offset','Pertubation phase offset',iunit) + call write_inopt(perturb_direction, 'FLRWSolver::FLRW_perturb_direction','Pertubation direction',iunit) + call write_inopt(radiation_dominated, 'radiation_dominated','Radiation dominated universe (yes/no)',iunit) + call write_inopt(perturb_wavelength,'FLRWSolver::single_perturb_wavelength','Perturbation wavelength',iunit) + call write_inopt(ilattice,'ilattice','lattice type (1=cubic, 2=closepacked)',iunit) + close(iunit) + +end subroutine write_setupfile + +!------------------------------------------------------------------------ +! +! read setup file +! +!------------------------------------------------------------------------ +subroutine read_setupfile(filename,ierr) + use infile_utils, only:open_db_from_file,inopts,read_inopt,close_db + use units, only:select_unit + use io, only:error + character(len=*), intent(in) :: filename + integer, intent(out) :: ierr + integer, parameter :: iunit = 21 + integer :: nerr + type(inopts), allocatable :: db(:) + + print "(a)",' reading setup options from '//trim(filename) + nerr = 0 + ierr = 0 + call open_db_from_file(db,filename,iunit,ierr) + ! + ! units + ! + call read_inopt(mass_unit,'mass_unit',db,errcount=nerr) + call read_inopt(dist_unit,'dist_unit',db,errcount=nerr) + ! + ! boundaries + ! + call read_inopt(xmini,'CoordBase::xmin',db,errcount=nerr) + call read_inopt(xmaxi,'CoordBase::xmax',db,min=xmini,errcount=nerr) + call read_inopt(ymini,'CoordBase::ymin',db,errcount=nerr) + call read_inopt(ymaxi,'CoordBase::ymax',db,min=ymini,errcount=nerr) + call read_inopt(zmini,'CoordBase::zmin',db,errcount=nerr) + call read_inopt(zmaxi,'CoordBase::zmax',db,min=zmini,errcount=nerr) + ! + ! other parameters + ! + call read_inopt(npartx,'nx',db,min=8,errcount=nerr) + call read_inopt(rhozero,'rhozero',db,min=0.,errcount=nerr) + call read_inopt(cs0,'cs0',db,min=0.,errcount=nerr) + + call read_inopt(perturb_direction,'FLRWSolver::FLRW_perturb_direction',db,errcount=nerr) + call read_inopt(ampl, 'FLRWSolver::phi_amplitude',db,errcount=nerr) + call read_inopt(phaseoffset,'FLRWSolver::phi_phase_offset',db,errcount=nerr) + call read_inopt(ilattice,'ilattice',db,min=1,max=2,errcount=nerr) + ! TODO Work out why this doesn't read in correctly + call read_inopt(perturb,'FLRWSolver::FLRW_perturb',db,errcount=nerr) + call read_inopt(radiation_dominated,'radiation_dominated',db,errcount=nerr) + call read_inopt(perturb_wavelength,'FLRWSolver::single_perturb_wavelength',db,errcount=nerr) + !print*, db + call close_db(db) + + if (nerr > 0) then + print "(1x,i2,a)",nerr,' error(s) during read of setup file: re-writing...' + ierr = nerr + endif + ! + ! parse units + ! + call select_unit(mass_unit,umass,nerr) + if (nerr /= 0) then + call error('setup_unifdis','mass unit not recognised') + ierr = ierr + 1 + endif + call select_unit(dist_unit,udist,nerr) + if (nerr /= 0) then + call error('setup_unifdis','length unit not recognised') + ierr = ierr + 1 + endif + + +end subroutine read_setupfile + +end module setup diff --git a/src/setup/setup_flrwpspec.f90 b/src/setup/setup_flrwpspec.f90 new file mode 100644 index 000000000..eef12efc8 --- /dev/null +++ b/src/setup/setup_flrwpspec.f90 @@ -0,0 +1,607 @@ +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.bitbucket.io/ ! +!--------------------------------------------------------------------------! +module setup +! +! Setup routine for realistic cosmological initial conditions based +! on the Zeldovich approximation. +! Requries velocity files generated from a powerspectrum. +! +! :References: None +! +! :Owner: Spencer Magnall +! +! :Runtime parameters: +! - cs0 : *initial sound speed in code units* +! - dist_unit : *distance unit (e.g. au)* +! - ilattice : *lattice type (1=cubic, 2=closepacked)* +! - mass_unit : *mass unit (e.g. solarm)* +! - nx : *number of particles in x direction* +! - radiation_dominated : *Radiation dominated universe (yes/no)* +! - rhozero : *initial density in code units* +! +! :Dependencies: boundary, dim, eos_shen, infile_utils, io, mpidomain, +! mpiutils, options, part, physcon, prompting, setup_params, stretchmap, +! unifdis, units, utils_gr +! + use dim, only:use_dust + use setup_params, only:rhozero + use physcon, only:radconst + implicit none + public :: setpart + + integer :: npartx,ilattice + real :: cs0,xmini,xmaxi,ymini,ymaxi,zmini,zmaxi,ampl,phaseoffset + character(len=20) :: dist_unit,mass_unit,perturb_direction,perturb,radiation_dominated + real :: perturb_wavelength + real(kind=8) :: udist,umass + + private + +contains + +!---------------------------------------------------------------- +!+ +! setup for uniform particle distributions +!+ +!---------------------------------------------------------------- +subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,time,fileprefix) + use dim, only:maxvxyzu,gr + use setup_params, only:npart_total + use io, only:master + use unifdis, only:set_unifdis,rho_func!,mass_func + use boundary, only:xmin,ymin,zmin,xmax,ymax,zmax,dxbound,dybound,dzbound,set_boundary,cross_boundary + use part, only:periodic + use physcon, only:years,pc,solarm + use units, only:set_units + use mpidomain, only:i_belong + use stretchmap, only:set_density_profile + use utils_gr, only:perturb_metric, get_u0, get_sqrtg + !use cons2primsolver, only:primative2conservative + + integer, intent(in) :: id + integer, intent(inout) :: npart + integer, intent(out) :: npartoftype(:) + real, intent(out) :: xyzh(:,:) + real, intent(out) :: massoftype(:) + real, intent(out) :: polyk,gamma + real, intent(inout) :: hfact + real, intent(inout) :: time + character(len=20), intent(in) :: fileprefix + real, intent(out) :: vxyzu(:,:) + character(len=40) :: filename,lattice,pspec_filename1,pspec_filename2,pspec_filename3 + real :: totmass,deltax,pi + integer :: i,ierr,ncross + logical :: iexist,isperiodic(3) + real :: length, c1,c3 + real :: hub + real :: last_scattering_temp + real :: scale_factor,gradphi(3),vxyz(3),dxgrid,gridorigin + integer :: nghost, gridres, gridsize + real, allocatable :: vxgrid(:,:,:),vygrid(:,:,:),vzgrid(:,:,:) + + ! + !--general parameters + ! + perturb_wavelength = 0. + time = 0. + if (maxvxyzu < 4) then + gamma = 1. + else + ! 4/3 for radiation dominated case + ! irrelevant for + gamma = 4./3. + endif + ! Redefinition of pi to fix numerical error + pi = 4.*atan(1.) + ! + ! default units + ! + mass_unit = 'solarm' + dist_unit = 'mpc' + ! + ! set boundaries to default values + ! + xmini = xmin; xmaxi = xmax + ymini = ymin; ymaxi = ymax + zmini = zmin; zmaxi = zmax + ! + ! set default values for input parameters + ! + npartx = 64 + length = 0. + ilattice = 1 + perturb = '"no"' + perturb_direction = '"none"' + radiation_dominated = '"no"' + ampl = 0. + + ! Ideally this should read the values of the box length + ! and initial Hubble parameter from the par file. + ! Then it should be set using the Friedmann equation: + !!!!!! rhozero = (3H^2)/(8*pi*a*a) + + hub = 10.553495658357338 + rhozero = 3. * hub**2 / (8. * pi) + phaseoffset = 0. + + ! Set some default values for the grid + nghost = 6 + gridres = 64 + + gridsize = nghost + gridres + gridorigin = 0. + xmax = 1. + dxgrid = xmax/gridres + gridorigin = gridorigin-3*dxgrid + + isperiodic = .true. + ncross = 0 + + + + ! Approx Temp of the CMB in Kelvins + last_scattering_temp = 3000 + last_scattering_temp = (rhozero/radconst)**(1./4.)*0.99999 + + ! Define some parameters for Linear pertubations + ! We assume ainit = 1, but this may not always be the case + c1 = 1./(4.*PI*rhozero) + !c2 = We set g(x^i) = 0 as we only want to extract the growing mode + c3 = - sqrt(1./(6.*PI*rhozero)) + + + if (gr) then + ! 0 Because dust? + cs0 = 0. + else + cs0 = 1. + endif + + ! get disc setup parameters from file or interactive setup + ! + filename=trim(fileprefix)//'.setup' + inquire(file=filename,exist=iexist) + if (iexist) then + !--read from setup file + call read_setupfile(filename,ierr) + if (id==master) call write_setupfile(filename) + if (ierr /= 0) then + stop + endif + elseif (id==master) then + call setup_interactive(id,polyk) + call write_setupfile(filename) + stop 'rerun phantomsetup after editing .setup file' + else + stop + endif + ! + ! set units and boundaries + ! + if (gr) then + call set_units(dist=udist,c=1.,G=1.) + else + call set_units(dist=udist,mass=umass,G=1.) + endif + call set_boundary(xmini,xmaxi,ymini,ymaxi,zmini,zmaxi) + + + allocate(vxgrid(gridsize,gridsize,gridsize)) + allocate(vygrid(gridsize,gridsize,gridsize)) + allocate(vzgrid(gridsize,gridsize,gridsize)) + ! + ! setup particles + ! + + npart = 0 + npart_total = 0 + length = xmaxi - xmini + deltax = length/npartx +! +! general parameters +! +! time should be read in from the par file + time = 0.18951066686763596 ! z~1000 +! lambda = perturb_wavelength*length +! kwave = (2.d0*pi)/lambda +! denom = length - ampl/kwave*(cos(kwave*length)-1.0) + rhozero = 3. * hub**2 / (8. * pi) + + + lattice = 'cubic' + + call set_unifdis(lattice,id,master,xmin,xmax,ymin,ymax,zmin,zmax,deltax,hfact,& + npart,xyzh,periodic,nptot=npart_total,mask=i_belong) + + npartoftype(:) = 0 + npartoftype(1) = npart + print*,' npart = ',npart,npart_total + + + totmass = rhozero*dxbound*dybound*dzbound + massoftype(1) = totmass/npart_total + if (id==master) print*,' particle mass = ',massoftype(1) + if (id==master) print*,' initial sound speed = ',cs0,' pressure = ',cs0**2/gamma + + + + if (maxvxyzu < 4 .or. gamma <= 1.) then + polyk = cs0**2 + else + polyk = 0. + endif + + pspec_filename1 = 'init_vel1_64.dat' + pspec_filename2 = 'init_vel2_64.dat' + pspec_filename3 = 'init_vel3_64.dat' + + ! Check if files exist otherwise skip and return flat space + if (.not. check_files(pspec_filename1,pspec_filename2,pspec_filename3)) then + print*, "Velocity files not found..." + print*, "Setting up flat space!" + return + endif + + + ! Read in velocities from vel file here + ! Should be made into a function at some point +! open(unit=444,file=pspec_filename,status='old') +! do k=1,gridsize +! do j=1,gridsize +! read(444,*) (vxgrid(i,j,k), i=1, 9) + +! enddo +! enddo +! close(444) + call read_veldata(vxgrid,pspec_filename1,gridsize) + call read_veldata(vygrid,pspec_filename2,gridsize) + call read_veldata(vzgrid,pspec_filename3,gridsize) + +! vxgrid = 1. +! vygrid = 2. +! vzgrid = 3. + !stop + do i=1,npart + ! Assign new particle possition + particle velocities here using the Zeldovich approximation: + ! Valid for Omega = 1 + ! x = q - a grad phi (1), where q is the non perturbed lattice point position + ! v = -aH grad phi (2) + ! Interpolate grid velocities to particles + ! big v vs small v? + ! Call interpolate from grid + !get_velocity_fromgrid(vxyz,pos) + ! CHECK THAT GRID ORIGIN IS CORRECT !!!!!!!!!!! + ! DO I NEED TO UPDATE THE GHOST CELLS?? + ! Get x velocity at particle position + call interpolate_val(xyzh(1:3,i),vxgrid,gridsize,gridorigin,dxgrid,vxyz(1)) + print*, "Finished x interp" + ! Get y velocity at particle position + call interpolate_val(xyzh(1:3,i),vygrid,gridsize,gridorigin,dxgrid,vxyz(2)) + ! Get z velocity at particle position + call interpolate_val(xyzh(1:3,i),vzgrid,gridsize,gridorigin,dxgrid,vxyz(3)) + + vxyzu(1:3,i) = vxyz + print*, vxyz + ! solve eqn (2) for grad phi + ! This is probally not constant?? + scale_factor = 1. + gradphi = -vxyz/(scale_factor*hub) + ! Set particle pos + xyzh(1:3,i) = xyzh(1:3,i) - scale_factor*gradphi + ! Apply periodic boundary conditions to particle position + call cross_boundary(isperiodic,xyzh(1:3,i),ncross) + + ! Calculate a new smoothing length?? Since the particle distrubtion has changed + + enddo + + + +end subroutine setpart + +!------------------------------------------------------------------------ +! +! interactive setup +! +!------------------------------------------------------------------------ +subroutine setup_interactive(id,polyk) + use io, only:master + use mpiutils, only:bcast_mpi + use dim, only:maxp,maxvxyzu + use prompting, only:prompt + use units, only:select_unit + integer, intent(in) :: id + real, intent(out) :: polyk + integer :: ierr + + if (id==master) then + ierr = 1 + do while (ierr /= 0) + call prompt('Enter mass unit (e.g. solarm,jupiterm,earthm)',mass_unit) + call select_unit(mass_unit,umass,ierr) + if (ierr /= 0) print "(a)",' ERROR: mass unit not recognised' + enddo + ierr = 1 + do while (ierr /= 0) + call prompt('Enter distance unit (e.g. au,pc,kpc,0.1pc)',dist_unit) + call select_unit(dist_unit,udist,ierr) + if (ierr /= 0) print "(a)",' ERROR: length unit not recognised' + enddo + + call prompt('enter xmin boundary',xmini) + call prompt('enter xmax boundary',xmaxi,xmini) + call prompt('enter ymin boundary',ymini) + call prompt('enter ymax boundary',ymaxi,ymini) + call prompt('enter zmin boundary',zmini) + call prompt('enter zmax boundary',zmaxi,zmini) + endif + ! + ! number of particles + ! + if (id==master) then + print*,' uniform setup... (max = ',nint((maxp)**(1/3.)),')' + call prompt('enter number of particles in x direction ',npartx,1) + endif + call bcast_mpi(npartx) + ! + ! mean density + ! + if (id==master) call prompt(' enter density (gives particle mass)',rhozero,0.) + call bcast_mpi(rhozero) + ! + ! sound speed in code units + ! + if (id==master) then + call prompt(' enter sound speed in code units (sets polyk)',cs0,0.) + endif + call bcast_mpi(cs0) + + ! + ! type of lattice + ! + if (id==master) then + call prompt(' select lattice type (1=cubic, 2=closepacked)',ilattice,1) + endif + call bcast_mpi(ilattice) +end subroutine setup_interactive + +!------------------------------------------------------------------------ +! +! write setup file +! +!------------------------------------------------------------------------ +subroutine write_setupfile(filename) + use infile_utils, only:write_inopt + character(len=*), intent(in) :: filename + integer :: iunit + + print "(/,a)",' writing setup options file '//trim(filename) + open(newunit=iunit,file=filename,status='replace',form='formatted') + write(iunit,"(a)") '# input file for uniform setup routine' + + write(iunit,"(/,a)") '# units' + call write_inopt(dist_unit,'dist_unit','distance unit (e.g. au)',iunit) + call write_inopt(mass_unit,'mass_unit','mass unit (e.g. solarm)',iunit) + ! + ! boundaries + ! + write(iunit,"(/,a)") '# boundaries' + call write_inopt(xmini,'CoordBase::xmin','xmin boundary',iunit) + call write_inopt(xmaxi,'CoordBase::xmax','xmax boundary',iunit) + call write_inopt(ymini,'CoordBase::ymin','ymin boundary',iunit) + call write_inopt(ymaxi,'CoordBase::ymax','ymax boundary',iunit) + call write_inopt(zmini,'CoordBase::zmin','zmin boundary',iunit) + call write_inopt(zmaxi,'CoordBase::zmax','zmax boundary',iunit) + + + + ! + ! other parameters + ! + write(iunit,"(/,a)") '# setup' + call write_inopt(npartx,'nx','number of particles in x direction',iunit) + call write_inopt(rhozero,'rhozero','initial density in code units',iunit) + call write_inopt(cs0,'cs0','initial sound speed in code units',iunit) + call write_inopt(perturb,'FLRWSolver::FLRW_perturb','Pertrubations of FLRW?',iunit) + call write_inopt(ampl,'FLRWSolver::phi_amplitude','Pertubation amplitude',iunit) + call write_inopt(phaseoffset,'FLRWSolver::phi_phase_offset','Pertubation phase offset',iunit) + call write_inopt(perturb_direction, 'FLRWSolver::FLRW_perturb_direction','Pertubation direction',iunit) + call write_inopt(radiation_dominated, 'radiation_dominated','Radiation dominated universe (yes/no)',iunit) + call write_inopt(perturb_wavelength,'FLRWSolver::single_perturb_wavelength','Perturbation wavelength',iunit) + call write_inopt(ilattice,'ilattice','lattice type (1=cubic, 2=closepacked)',iunit) + close(iunit) + +end subroutine write_setupfile + +!------------------------------------------------------------------------ +! +! read setup file +! +!------------------------------------------------------------------------ +subroutine read_setupfile(filename,ierr) + use infile_utils, only:open_db_from_file,inopts,read_inopt,close_db + use units, only:select_unit + use io, only:error + character(len=*), intent(in) :: filename + integer, intent(out) :: ierr + integer, parameter :: iunit = 21 + integer :: nerr + type(inopts), allocatable :: db(:) + + print "(a)",' reading setup options from '//trim(filename) + nerr = 0 + ierr = 0 + call open_db_from_file(db,filename,iunit,ierr) + ! + ! units + ! + call read_inopt(mass_unit,'mass_unit',db,errcount=nerr) + call read_inopt(dist_unit,'dist_unit',db,errcount=nerr) + ! + ! boundaries + ! + call read_inopt(xmini,'CoordBase::xmin',db,errcount=nerr) + call read_inopt(xmaxi,'CoordBase::xmax',db,min=xmini,errcount=nerr) + call read_inopt(ymini,'CoordBase::ymin',db,errcount=nerr) + call read_inopt(ymaxi,'CoordBase::ymax',db,min=ymini,errcount=nerr) + call read_inopt(zmini,'CoordBase::zmin',db,errcount=nerr) + call read_inopt(zmaxi,'CoordBase::zmax',db,min=zmini,errcount=nerr) + ! + ! other parameters + ! + call read_inopt(npartx,'nx',db,min=8,errcount=nerr) + call read_inopt(rhozero,'rhozero',db,min=0.,errcount=nerr) + call read_inopt(cs0,'cs0',db,min=0.,errcount=nerr) + + call read_inopt(perturb_direction,'FLRWSolver::FLRW_perturb_direction',db,errcount=nerr) + call read_inopt(ampl, 'FLRWSolver::phi_amplitude',db,errcount=nerr) + call read_inopt(phaseoffset,'FLRWSolver::phi_phase_offset',db,errcount=nerr) + call read_inopt(ilattice,'ilattice',db,min=1,max=2,errcount=nerr) + ! TODO Work out why this doesn't read in correctly + call read_inopt(perturb,'FLRWSolver::FLRW_perturb',db,errcount=nerr) + call read_inopt(radiation_dominated,'radiation_dominated',db,errcount=nerr) + call read_inopt(perturb_wavelength,'FLRWSolver::single_perturb_wavelength',db,errcount=nerr) + !print*, db + call close_db(db) + + if (nerr > 0) then + print "(1x,i2,a)",nerr,' error(s) during read of setup file: re-writing...' + ierr = nerr + endif + ! + ! parse units + ! + call select_unit(mass_unit,umass,nerr) + if (nerr /= 0) then + call error('setup_unifdis','mass unit not recognised') + ierr = ierr + 1 + endif + call select_unit(dist_unit,udist,nerr) + if (nerr /= 0) then + call error('setup_unifdis','length unit not recognised') + ierr = ierr + 1 + endif + + +end subroutine read_setupfile + +subroutine read_veldata(velarray,vfile,gridsize) + integer, intent(in) :: gridsize + character(len=20),intent(in) :: vfile + real,intent(out) :: velarray(:,:,:) + integer :: i,j,k,iu + + open(newunit=iu,file=vfile,status='old') + do k=1,gridsize + do j=1,gridsize + read(iu,*) (velarray(i,j,k), i=1, gridsize) + enddo + enddo + close(iu) + print*, "Finished reading ", vfile + +end subroutine read_veldata + +subroutine interpolate_val(position,valgrid,gridsize,gridorigin,dxgrid,val) + ! Subroutine to interpolate quanities to particle positions given a cube + ! Note we have assumed that the grid will always be cubic!!!! + use eos_shen, only:linear_interpolator_one_d + real, intent(in) :: valgrid(:,:,:) + real, intent(inout) :: position(3) + real, intent(inout) :: dxgrid,gridorigin + integer, intent(in) :: gridsize + real, intent(out) :: val + integer :: xupper,yupper,zupper,xlower,ylower,zlower + real :: xlowerpos,ylowerpos,zlowerpos!,xupperpos,yupperpos,zupperpos + real :: interptmp(7) + real :: xd,yd,zd + + + + call get_grid_neighbours(position,gridorigin,dxgrid,xlower,ylower,zlower) + + print*,"Neighbours: ", xlower,ylower,zlower + print*,"Position: ", position + ! This is not true as upper neighbours on the boundary will be on the side + ! take a mod of grid size + xupper = mod(xlower + 1, gridsize) + yupper = mod(ylower + 1, gridsize) + zupper = mod(zlower + 1, gridsize) + ! xupper - xlower should always just be dx provided we are using a uniform grid + ! xd = (position(1) - xlower)/(xupper - xlower) + ! yd = (position(2) - ylower)/(yupper - ylower) + ! zd = (position(3) - zlower)/(zupper - zlower) + xlowerpos = gridorigin + (xlower-1)*dxgrid + ylowerpos = gridorigin + (ylower-1)*dxgrid + zlowerpos = gridorigin + (zlower-1)*dxgrid + + xd = (position(1) - xlowerpos)/(dxgrid) + yd = (position(2) - ylowerpos)/(dxgrid) + zd = (position(3) - zlowerpos)/(dxgrid) + + interptmp = 0. + + call linear_interpolator_one_d(valgrid(xlower,ylower,zlower), & + valgrid(xlower+1,ylower,zlower),xd,interptmp(1)) + call linear_interpolator_one_d(valgrid(xlower,ylower,zlower+1), & + valgrid(xlower+1,ylower,zlower+1),xd,interptmp(2)) + call linear_interpolator_one_d(valgrid(xlower,ylower+1,zlower), & + valgrid(xlower+1,ylower+1,zlower),xd,interptmp(3)) + call linear_interpolator_one_d(valgrid(xlower,ylower+1,zlower+1), & + valgrid(xlower+1,ylower+1,zlower+1),xd,interptmp(4)) + ! Interpolate along y + call linear_interpolator_one_d(interptmp(1),interptmp(3),yd,interptmp(5)) + call linear_interpolator_one_d(interptmp(2),interptmp(4),yd,interptmp(6)) + ! Interpolate along z + call linear_interpolator_one_d(interptmp(5),interptmp(6),zd,interptmp(7)) + + val = interptmp(7) + +end subroutine interpolate_val + +subroutine get_grid_neighbours(position,gridorigin,dx,xlower,ylower,zlower) + ! TODO IDEALLY THIS SHOULDN'T BE HERE AND SHOULD BE IN A UTILS MODULE + ! WITH THE VERSION USED IN METRIC_ET + real, intent(in) :: position(3), gridorigin + real, intent(in) :: dx + integer, intent(out) :: xlower,ylower,zlower + + ! Get the lower grid neighbours of the position + ! If this is broken change from floor to int + ! How are we handling the edge case of a particle being + ! in exactly the same position as the grid? + ! Hopefully having different grid sizes in each direction + ! Doesn't break the lininterp + xlower = floor((position(1)-gridorigin)/dx) + print*, "pos x: ", position(1) + print*, "gridorigin: ", gridorigin + print*, "dx: ", dx + ylower = floor((position(2)-gridorigin)/dx) + zlower = floor((position(3)-gridorigin)/dx) + + ! +1 because fortran + xlower = xlower + 1 + ylower = ylower + 1 + zlower = zlower + 1 + + +end subroutine get_grid_neighbours + +logical function check_files(file1,file2,file3) + character(len=*), intent(in) :: file1,file2,file3 + logical :: file1_exist, file2_exist, file3_exist + + inquire(file=file1,exist=file1_exist) + inquire(file=file2,exist=file2_exist) + inquire(file=file3,exist=file3_exist) + + if ((.not. file1_exist) .or. (.not. file2_exist) .or. (.not. file3_exist)) then + check_files = .false. + endif +end function check_files + +end module setup diff --git a/src/setup/stretchmap.f90 b/src/setup/stretchmap.f90 index 94aac53c6..733c14497 100644 --- a/src/setup/stretchmap.f90 +++ b/src/setup/stretchmap.f90 @@ -30,11 +30,12 @@ module stretchmap public :: set_density_profile public :: get_mass_r public :: rho_func + public :: mass_func - integer, private :: ngrid = 1024 ! number of points used when integrating rho to get mass + integer, private :: ngrid = 8192 ! number of points used when integrating rho to get mass integer, parameter, private :: maxits = 100 ! max number of iterations integer, parameter, private :: maxits_nr = 30 ! max iterations with Newton-Raphson - real, parameter, private :: tol = 1.e-9 ! tolerance on iterations + real, parameter, private :: tol = 1.e-10 ! tolerance on iterations integer, parameter, public :: ierr_zero_size_density_table = 1, & ! error code ierr_memory_allocation = 2, & ! error code ierr_table_size_differs = 3, & ! error code @@ -45,11 +46,17 @@ real function rho_func(x) end function rho_func end interface + abstract interface + real function mass_func(x,xmin) + real, intent(in) :: x, xmin + end function mass_func + end interface + private contains -subroutine set_density_profile(np,xyzh,min,max,rhofunc,rhotab,xtab,start,geom,coord,verbose,err) +subroutine set_density_profile(np,xyzh,min,max,rhofunc,massfunc,rhotab,xtab,start,geom,coord,verbose,err) ! ! Subroutine to implement the stretch mapping procedure ! @@ -91,6 +98,7 @@ subroutine set_density_profile(np,xyzh,min,max,rhofunc,rhotab,xtab,start,geom,co real, intent(inout) :: xyzh(:,:) real, intent(in) :: min,max procedure(rho_func), pointer, optional :: rhofunc + procedure(mass_func), pointer, optional :: massfunc real, intent(in), optional :: rhotab(:),xtab(:) integer, intent(in), optional :: start, geom, coord logical, intent(in), optional :: verbose @@ -101,13 +109,16 @@ subroutine set_density_profile(np,xyzh,min,max,rhofunc,rhotab,xtab,start,geom,co real, allocatable :: xtable(:),masstab(:) integer :: i,its,igeom,icoord,istart,nt,nerr,ierr logical :: is_r, is_rcyl, bisect, isverbose - logical :: use_rhotab + logical :: use_rhotab, use_massfunc isverbose = .true. use_rhotab = .false. + use_massfunc = .false. + if (present(verbose)) isverbose = verbose if (present(rhotab)) use_rhotab = .true. - + if (present(massfunc)) use_massfunc = .true. + if (use_massfunc) print "(a)", 'Using massfunc rather than numerically-integrated table' if (present(rhofunc) .or. present(rhotab)) then if (isverbose) print "(a)",' >>>>>> s t r e t c h m a p p i n g <<<<<<' ! @@ -176,6 +187,8 @@ subroutine set_density_profile(np,xyzh,min,max,rhofunc,rhotab,xtab,start,geom,co totmass = get_mass_r(rhofunc,xmax,xmin) elseif (is_rcyl) then totmass = get_mass_rcyl(rhofunc,xmax,xmin) + elseif (use_massfunc) then + totmass = massfunc(xmax,min) else totmass = get_mass(rhofunc,xmax,xmin) endif @@ -203,8 +216,8 @@ subroutine set_density_profile(np,xyzh,min,max,rhofunc,rhotab,xtab,start,geom,co nerr = 0 !$omp parallel do default(none) & - !$omp shared(np,xyzh,rhozero,igeom,use_rhotab,rhotab,xtable,masstab,nt) & - !$omp shared(xmin,xmax,totmass,icoord,is_r,is_rcyl,istart,rhofunc) & + !$omp shared(np,xyzh,rhozero,igeom,use_rhotab,use_massfunc,rhotab,xtable,masstab,nt) & + !$omp shared(xmin,xmax,totmass,icoord,is_r,is_rcyl,istart,rhofunc,massfunc) & !$omp private(x,xold,xt,fracmassold,its,xprev,xi,hi,rhoi) & !$omp private(func,dfunc,xminbisect,xmaxbisect,bisect) & !$omp reduction(+:nerr) @@ -239,6 +252,8 @@ subroutine set_density_profile(np,xyzh,min,max,rhofunc,rhotab,xtab,start,geom,co func = get_mass_r(rhofunc,xi,xmin) elseif (is_rcyl) then func = get_mass_rcyl(rhofunc,xi,xmin) + elseif (use_massfunc) then + func = massfunc(xi,xmin) else func = get_mass(rhofunc,xi,xmin) endif @@ -266,6 +281,9 @@ subroutine set_density_profile(np,xyzh,min,max,rhofunc,rhotab,xtab,start,geom,co elseif (is_rcyl) then func = get_mass_rcyl(rhofunc,xi,xmin) - fracmassold dfunc = 2.*pi*xi*rhofunc(xi) + elseif (use_massfunc) then + func = massfunc(xi,xmin) - fracmassold + dfunc = rhofunc(xi) else func = get_mass(rhofunc,xi,xmin) - fracmassold dfunc = rhofunc(xi) diff --git a/src/utils/einsteintk_utils.f90 b/src/utils/einsteintk_utils.f90 new file mode 100644 index 000000000..880ac3096 --- /dev/null +++ b/src/utils/einsteintk_utils.f90 @@ -0,0 +1,172 @@ +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.bitbucket.io/ ! +!--------------------------------------------------------------------------! +module einsteintk_utils +! +! einsteintk_utils +! +! :References: None +! +! :Owner: Spencer Magnall +! +! :Runtime parameters: None +! +! :Dependencies: part +! + implicit none + real, allocatable :: gcovgrid(:,:,:,:,:) + real, allocatable :: gcongrid(:,:,:,:,:) + real, allocatable :: sqrtggrid(:,:,:) + real, allocatable :: tmunugrid(:,:,:,:,:) + real, allocatable :: rhostargrid(:,:,:) + real, allocatable :: pxgrid(:,:,:,:) + real, allocatable :: entropygrid(:,:,:) + real, allocatable :: metricderivsgrid(:,:,:,:,:,:) + real :: dxgrid(3), gridorigin(3), boundsize(3) + integer :: gridsize(3) + logical :: gridinit = .false. + logical :: exact_rendering + character(len=128) :: logfilestor,evfilestor,dumpfilestor,infilestor +contains +subroutine init_etgrid(nx,ny,nz,dx,dy,dz,originx,originy,originz) + integer, intent(in) :: nx,ny,nz + real, intent(in) :: dx,dy,dz,originx,originy,originz + + gridsize(1) = nx + gridsize(2) = ny + gridsize(3) = nz + + dxgrid(1) = dx + dxgrid(2) = dy + dxgrid(3) = dz + + gridorigin(1) = originx + gridorigin(2) = originy + gridorigin(3) = originz + + + allocate(gcovgrid(0:3,0:3,nx,ny,nz)) + allocate(gcongrid(0:3,0:3,nx,ny,nz)) + allocate(sqrtggrid(nx,ny,nz)) + + ! Will need to delete this at somepoint + ! For now it is the simplest way + allocate(tmunugrid(0:3,0:3,nx,ny,nz)) + + allocate(pxgrid(3,nx,ny,nz)) + + allocate(rhostargrid(nx,ny,nz)) + + allocate(entropygrid(nx,ny,nz)) + + ! metric derivs are stored in the form + ! mu comp, nu comp, deriv, gridx,gridy,gridz + ! Note that this is only the spatial derivs of + ! the metric and we will need an additional array + ! for time derivs + allocate(metricderivsgrid(0:3,0:3,3,nx,ny,nz)) + + gridinit = .true. + !exact_rendering = exact + +end subroutine init_etgrid + +subroutine print_etgrid() + ! Subroutine for printing quantities of the ET grid + + print*, "Grid spacing (x,y,z) is : ", dxgrid + print*, "Grid origin (x,y,z) is: ", gridorigin + print*, "Covariant metric tensor of the grid is: ", gcovgrid(:,:,1,1,1) + +end subroutine print_etgrid + +subroutine get_particle_rhs(i,vx,vy,vz,fx,fy,fz,e_rhs) + use part, only: vxyzu,fext!,fxyzu + integer, intent(in) :: i + real, intent(out) :: vx,vy,vz,fx,fy,fz,e_rhs + + !vxyz + vx = vxyzu(1,i) + vy = vxyzu(2,i) + vz = vxyzu(3,i) + + fx = fext(1,i) + fy = fext(2,i) + fz = fext(3,i) + + + ! de/dt + e_rhs = 0. + +end subroutine get_particle_rhs + +subroutine get_particle_val(i,x,y,z,px,py,pz,e) + use part, only: xyzh, pxyzu + integer, intent(in) :: i + real, intent(out) :: x,y,z,px,py,pz,e + + !xyz + x = xyzh(1,i) + y = xyzh(2,i) + z = xyzh(3,i) + + ! p + px = pxyzu(1,i) + py = pxyzu(2,i) + pz = pxyzu(3,i) + + ! e + ! ??? + e = pxyzu(4,i) + +end subroutine get_particle_val + +subroutine set_particle_val(i,x,y,z,px,py,pz,e) + use part, only: xyzh, pxyzu + integer, intent(in) :: i + real, intent(in) :: x,y,z,px,py,pz,e + ! Subroutine for setting the particle values in phantom + ! using the values stored in einstein toolkit before a dump + + !xyz + xyzh(1,i) = x + xyzh(2,i) = y + xyzh(3,i) = z + + ! p + pxyzu(1,i) = px + pxyzu(2,i) = py + pxyzu(3,i) = pz + pxyzu(4,i) = e + + +end subroutine set_particle_val + +subroutine get_phantom_dt(dtout) + use part, only:xyzh + real, intent(out) :: dtout + real, parameter :: safety_fac = 0.2 + real :: minh + + ! Get the smallest smoothing length + minh = minval(xyzh(4,:)) + + ! Courant esque condition from Rosswog 2021+ + ! Since c is allways one in our units + dtout = safety_fac*minh + print*, "dtout phantom: ", dtout + + +end subroutine get_phantom_dt + +subroutine set_rendering(flag) + logical, intent(in) :: flag + + exact_rendering = flag + +end subroutine set_rendering + +end module einsteintk_utils diff --git a/src/utils/einsteintk_wrapper.f90 b/src/utils/einsteintk_wrapper.f90 new file mode 100644 index 000000000..f7b5282e2 --- /dev/null +++ b/src/utils/einsteintk_wrapper.f90 @@ -0,0 +1,542 @@ +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.bitbucket.io/ ! +!--------------------------------------------------------------------------! +module einsteintk_wrapper +! +! einsteintk_wrapper +! +! :References: None +! +! :Owner: Spencer Magnall +! +! :Runtime parameters: None +! +! :Dependencies: cons2prim, densityforce, deriv, einsteintk_utils, evwrite, +! extern_gr, fileutils, initial, io, linklist, metric, metric_tools, +! mpiutils, part, readwrite_dumps, timestep, tmunu2grid +! + implicit none +contains + +subroutine init_et2phantom(infilestart,dt_et,nophantompart,dtout) + ! Wrapper that intialises phantom + ! Intended to hide all of the inner works of phantom from ET + ! Majority of the code from HelloHydro_init has been moved here + + use io, only:id,master,nprocs,set_io_unit_numbers,die + use mpiutils, only:init_mpi,finalise_mpi + use initial, only:initialise,finalise,startrun,endrun + !use evolve, only:evol_init + use tmunu2grid + use einsteintk_utils + use extern_gr + use metric + use part, only:npart!, tmunus + + + implicit none + character(len=*), intent(in) :: infilestart + real, intent(in) :: dt_et + integer, intent(inout) :: nophantompart + real, intent(out) :: dtout + !character(len=500) :: logfile,evfile,dumpfile,path + !integer :: i,j,k,pathstringlength + + ! For now we just hardcode the infile, to see if startrun actually works! + ! I'm not sure what the best way to actually do this is? + ! Do we store the phantom.in file in par and have it read from there? + !infile = "/Users/spencer/phantomET/phantom/test/flrw.in" + !infile = trim(infile)//'.in' + !print*, "phantom_path: ", phantom_path + !infile = phantom_path // "flrw.in" + !infile = trim(path) // "flrw.in" + !infile = 'flrw.in' + !infile = trim(infile) + !print*, "Phantom path is: ", path + !print*, "Infile is: ", infile + ! Use system call to copy phantom files to simulation directory + ! This is a digusting temporary fix + !call SYSTEM('cp ~/phantomET/phantom/test/flrw* ./') + + ! The infile from ET + infilestor = infilestart + + ! We should do everything that is done in phantom.f90 + + ! Setup mpi + id=0 + call init_mpi(id,nprocs) + ! setup io + call set_io_unit_numbers + ! routine that starts a phantom run + print*, "Start run called!" + ! Do we want to pass dt in here?? + call startrun(infilestor,logfilestor,evfilestor,dumpfilestor) + print*, "Start run finished!" + !print*, "tmunugrid: ", tmunugrid(1,1,6,6,6) + !stop + ! Intialises values for the evol routine: t, dt, etc.. + !call evol_init(infilestor,logfilestor,evfilestor,dumpfilestor,dt_et,nophantompart) + !print*, "Evolve init finished!" + nophantompart = npart + ! Calculate the stress energy tensor for each particle + ! Might be better to do this in evolve init + !call get_tmunugrid_all + ! Calculate the stress energy tensor + call get_metricderivs_all(dtout,dt_et) ! commented out to try and fix prim2cons + !call get_tmunu_all(npart,xyzh,metrics,vxyzu,metricderivs,dens,tmunus) ! commented out to try and fix prim2cons + !call get_tmunu_all_exact(npart,xyzh,metrics,vxyzu,metricderivs,dens,tmunus) + ! Interpolate stress energy tensor from particles back + ! to grid + !call get_tmunugrid_all(npart,xyzh,vxyzu,tmunus,calc_cfac=.true.) ! commented out to try and fix cons2prim + + call get_phantom_dt(dtout) + +end subroutine init_et2phantom + +subroutine init_et2phantomgrid(nx,ny,nz,originx,originy,originz,dx,dy,dz) + use einsteintk_utils + integer, intent(in) :: nx,ny,nz ! The maximum values of the grid in each dimension + real(8), intent(in) :: originx, originy, originz ! The origin of grid + real(8), intent(in) :: dx, dy, dz ! Grid spacing in each dimension + !integer, intent(in) :: boundsizex, boundsizey, boundsizez + + ! Setup metric grid + call init_etgrid(nx,ny,nz,originx,originy,originz,dx,dy,dz) + +end subroutine init_et2phantomgrid + +subroutine init_phantom2et() + ! Subroutine +end subroutine init_phantom2et + +subroutine et2phantom(rho,nx,ny,nz) + integer, intent(in) :: nx, ny, nz + real, intent(in) :: rho(nx,ny,nz) + + print*, "Grid limits: ", nx, ny, nz + ! get mpi thread number + ! send grid limits +end subroutine et2phantom + + ! DONT THINK THIS IS USED ANYWHERE!!! + ! subroutine step_et2phantom(infile,dt_et) + ! use einsteintk_utils + ! use evolve, only:evol_step + ! use tmunu2grid + ! character(len=*), intent(in) :: infile + ! real, intent(inout) :: dt_et + ! character(len=500) :: logfile,evfile,dumpfile,path + + + ! ! Print the values of logfile, evfile, dumpfile to check they are sensible + ! !print*, "logfile, evfile, dumpfile: ", logfile, evfile, dumpfile + ! print*, "stored values of logfile, evfile, dumpfile: ", logfilestor, evfilestor, dumpfilestor + + ! ! Interpolation stuff + ! ! Call et2phantom (construct global grid, metric, metric derivs, determinant) + ! ! Run phantom for a step + ! call evol_step(infile,logfilestor,evfilestor,dumpfilestor,dt_et) + ! ! Interpolation stuff back to et + ! !call get_tmunugrid_all() + ! ! call phantom2et (Tmunu_grid) + + ! end subroutine step_et2phantom + +subroutine phantom2et() + ! should take in the cctk_array for tmunu?? + ! Is it better if this routine is just + ! Calculate stress energy tensor for each particle + + ! Perform kernel interpolation from particles to grid positions + +end subroutine phantom2et + +subroutine step_et2phantom_MoL(infile,dt_et,dtout) + use part, only:xyzh,vxyzu,pxyzu,dens,metrics, npart, eos_vars + use cons2prim, only: cons2primall + use deriv + use extern_gr + use tmunu2grid + use einsteintk_utils, only: get_phantom_dt + character(len=*), intent(in) :: infile + real, intent(inout) :: dt_et + real, intent(out) :: dtout + real :: vbefore,vafter + + ! Metric should have already been passed in + ! and interpolated + ! Call get_derivs global + call get_derivs_global + + ! Get metric derivs + call get_metricderivs_all(dtout,dt_et) + ! Store our particle quantities somewhere / send them to ET + ! Cons2prim after moving the particles with the external force + vbefore = vxyzu(1,1) + call cons2primall(npart,xyzh,metrics,pxyzu,vxyzu,dens,eos_vars) + vafter = vxyzu(1,1) + + ! Does get_derivs_global perform a stress energy calc?? + ! If not do that here + + ! Perform the calculation of the stress energy tensor + ! Interpolate the stress energy tensor back to the ET grid! + ! Calculate the stress energy tensor + ! Interpolate stress energy tensor from particles back + ! to grid + call get_phantom_dt(dtout) + + +end subroutine step_et2phantom_MoL + +subroutine et2phantom_tmunu() + use part, only:npart,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,& + Bevol,rad,radprop,eos_vars,pxyzu,dens,metrics,tmunus,metricderivs,& + igas,rhoh,alphaind,dvdx,gradh + !use part, only:xyzh,vxyzu,fxyzu,pxyzu,dens,metricderivs, metrics, npart, tmunus,eos_vars + use cons2prim, only: cons2primall + use deriv + use extern_gr + use tmunu2grid + use einsteintk_utils, only: get_phantom_dt,rhostargrid,tmunugrid + use metric_tools, only:init_metric + use densityforce, only:densityiterate + use linklist, only:set_linklist + + real :: stressmax + real(kind=16) :: cfac + + stressmax = 0. + + ! Also probably need to pack the metric before I call things + call init_metric(npart,xyzh,metrics) + ! Might be better to just do this in get derivs global with a number 2 call? + ! Rebuild the tree + call set_linklist(npart,npart,xyzh,vxyzu) + ! Apparently init metric needs to be called again??? + !call init_metric(npart,xyzh,metrics) + ! Calculate the cons density + call densityiterate(1,npart,npart,xyzh,vxyzu,divcurlv,divcurlB,Bevol,& + stressmax,fxyzu,fext,alphaind,gradh,rad,radprop,dvdx) + ! Get primative variables for tmunu + call cons2primall(npart,xyzh,metrics,pxyzu,vxyzu,dens,eos_vars) + + call get_tmunu_all(npart,xyzh,metrics,vxyzu,metricderivs,dens,tmunus) + ! Interpolate stress energy tensor from particles back + ! to grid + call get_tmunugrid_all(npart,xyzh,vxyzu,tmunus) + + ! Interpolate density to grid + call phantom2et_rhostar + + ! Density check vs particles + call check_conserved_dens(rhostargrid,cfac) + + ! Correct Tmunu + ! Convert to 8byte real to stop compiler warning + tmunugrid = real(cfac)*tmunugrid + + +end subroutine et2phantom_tmunu + +subroutine phantom2et_consvar() + use part, only:npart,xyzh,vxyzu,fxyzu,fext,divcurlv,divcurlB,& + Bevol,rad,radprop,metrics,igas,rhoh,alphaind,dvdx,gradh + use densityforce, only:densityiterate + use metric_tools, only:init_metric + use linklist, only:set_linklist + use einsteintk_utils, only:rhostargrid,pxgrid,entropygrid + use tmunu2grid, only:check_conserved_dens + + real :: stressmax + real(kind=16) :: cfac + + ! Init metric + call init_metric(npart,xyzh,metrics) + + ! Might be better to just do this in get derivs global with a number 2 call? + ! Rebuild the tree + call set_linklist(npart,npart,xyzh,vxyzu) + ! Apparently init metric needs to be called again??? + call init_metric(npart,xyzh,metrics) + ! Calculate the cons density + call densityiterate(1,npart,npart,xyzh,vxyzu,divcurlv,divcurlB,Bevol,& + stressmax,fxyzu,fext,alphaind,gradh,rad,radprop,dvdx) + + ! Interpolate density to grid + call phantom2et_rhostar + + ! Interpolate momentum to grid + call phantom2et_momentum + + ! Interpolate entropy to grid + call phantom2et_entropy + + + ! Conserved quantity checks + corrections + + ! Density check vs particles + call check_conserved_dens(rhostargrid,cfac) + + ! Momentum check vs particles + + ! Correct momentum and Density + ! Conversion of cfac to 8byte real to avoid + ! compiler warning + rhostargrid = real(cfac)*rhostargrid + pxgrid = real(cfac)*pxgrid + entropygrid = real(cfac)*entropygrid + + +end subroutine phantom2et_consvar + +subroutine phantom2et_rhostar() + use part, only:xyzh,npart,& + igas, massoftype,rhoh + use cons2prim, only: cons2primall + use deriv + use extern_gr + use tmunu2grid + use einsteintk_utils, only: get_phantom_dt,rhostargrid + use metric_tools, only:init_metric + real :: dat(npart), h, pmass,rho + integer :: i + + + ! Get new cons density from new particle positions somehow (maybe)? + ! Set linklist to update the tree for neighbour finding + ! Calculate the density for the new particle positions + ! Call density iterate + + ! Interpolate from particles to grid + ! This can all go into its own function as it will essentially + ! be the same thing for all quantites + ! get particle data + ! get rho from xyzh and rhoh + ! Get the conserved density on the particles + dat = 0. + pmass = massoftype(igas) + ! $omp parallel do default(none) & + ! $omp shared(npart,xyzh,dat,pmass) & + ! $omp private(i,h,rho) + do i=1, npart + ! Get the smoothing length + h = xyzh(4,i) + ! Get pmass + + rho = rhoh(h,pmass) + dat(i) = rho + enddo + ! $omp end parallel do + rhostargrid = 0. + call interpolate_to_grid(rhostargrid,dat) + +end subroutine phantom2et_rhostar + +subroutine phantom2et_entropy() + use part, only:pxyzu,npart + use cons2prim, only: cons2primall + use deriv + use extern_gr + use tmunu2grid + use einsteintk_utils, only: get_phantom_dt,entropygrid + use metric_tools, only:init_metric + real :: dat(npart) + integer :: i + + + ! Get new cons density from new particle positions somehow (maybe)? + ! Set linklist to update the tree for neighbour finding + ! Calculate the density for the new particle positions + ! Call density iterate + + ! Interpolate from particles to grid + ! This can all go into its own function as it will essentially + ! be the same thing for all quantites + ! get particle data + ! get rho from xyzh and rhoh + ! Get the conserved density on the particles + dat = 0. + !$omp parallel do default(none) & + !$omp shared(npart,pxyzu,dat) & + !$omp private(i) + do i=1, npart + ! Entropy is the u component of pxyzu + dat(i) = pxyzu(4,i) + enddo + !$omp end parallel do + entropygrid = 0. + call interpolate_to_grid(entropygrid,dat) + +end subroutine phantom2et_entropy + +subroutine phantom2et_momentum() + use part, only:pxyzu, npart + use cons2prim, only: cons2primall + use deriv + use extern_gr + use tmunu2grid + use einsteintk_utils, only: get_phantom_dt,pxgrid + use metric_tools, only:init_metric + real :: dat(3,npart) + integer :: i + + + ! Pi is directly updated at the end of each MoL add + + ! Interpolate from particles to grid + ! get particle data for the x component of momentum + dat = 0. + !$omp parallel do default(none) & + !$omp shared(npart,pxyzu,dat) & + !$omp private(i) + do i=1, npart + dat(1,i) = pxyzu(1,i) + dat(2,i) = pxyzu(2,i) + dat(3,i) = pxyzu(3,i) + enddo + !$omp end parallel do + pxgrid = 0. + ! call interpolate 3d + ! In this case call it 3 times one for each vector component + ! px component + call interpolate_to_grid(pxgrid(1,:,:,:), dat(1,:)) + ! py component + call interpolate_to_grid(pxgrid(2,:,:,:), dat(2,:)) + ! pz component + call interpolate_to_grid(pxgrid(3,:,:,:),dat(3,:)) + + + +end subroutine phantom2et_momentum + + + + ! Subroutine for performing a phantom dump from einstein toolkit +subroutine et2phantom_dumphydro(time,dt_et,checkpointfile) + use einsteintk_utils + use evwrite, only:write_evfile,write_evlog + use readwrite_dumps, only:write_smalldump,write_fulldump + use fileutils, only:getnextfilename + use tmunu2grid, only:check_conserved_dens + real, intent(in) :: time, dt_et + !real(kind=16) :: cfac + !logical, intent(in), optional :: checkpoint + !integer, intent(in) :: checkpointno + character(*),optional, intent(in) :: checkpointfile + logical :: createcheckpoint + + if (present(checkpointfile)) then + createcheckpoint = .true. + else + createcheckpoint = .false. + endif + + ! Write EV_file + if (.not. createcheckpoint) then + call write_evfile(time,dt_et) + + evfilestor = getnextfilename(evfilestor) + logfilestor = getnextfilename(logfilestor) + dumpfilestor = getnextfilename(dumpfilestor) + call write_fulldump(time,dumpfilestor) + endif + + ! Write full dump + if (createcheckpoint) then + call write_fulldump(time,checkpointfile) + endif + + ! Quick and dirty write cfac to txtfile + + ! Density check vs particles +! call check_conserved_dens(rhostargrid,cfac) +! open(unit=777, file="cfac.txt", action='write', position='append') +! print*, time, cfac +! write(777,*) time, cfac +! close(unit=777) + +end subroutine et2phantom_dumphydro + + ! Provides the RHS derivs for a particle at index i +subroutine phantom2et_rhs(index, vx,vy,vz,fx,fy,fz,e_rhs) + use einsteintk_utils + real, intent(inout) :: vx,vy,vz,fx,fy,fz, e_rhs + integer, intent(in) :: index + + call get_particle_rhs(index,vx,vy,vz,fx,fy,fz,e_rhs) + +end subroutine phantom2et_rhs + +subroutine phantom2et_initial(index,x,y,z,px,py,pz,e) + use einsteintk_utils + real, intent(inout) :: x,y,z,px,py,pz,e + integer, intent(in) :: index + + call get_particle_val(index,x,y,z,px,py,pz,e) + +end subroutine phantom2et_initial + +subroutine et2phantom_setparticlevars(index,x,y,z,px,py,pz,e) + use einsteintk_utils + real, intent(inout) :: x,y,z,px,py,pz,e + integer, intent(in) :: index + + call set_particle_val(index,x,y,z,px,py,pz,e) + +end subroutine et2phantom_setparticlevars + + ! I really HATE this routine being here but it needs to be to fix dependency issues. +subroutine get_metricderivs_all(dtextforce_min,dt_et) + !use einsteintk_utils, only: metricderivsgrid + use part, only:npart,xyzh,vxyzu,dens,metrics,metricderivs,fext!,fxyzu + use timestep, only:bignumber,C_force + use extern_gr, only:get_grforce + use metric_tools, only:pack_metricderivs + real, intent(out) :: dtextforce_min + real, intent(in) :: dt_et + integer :: i + real :: pri,dtf + + pri = 0. + dtextforce_min = bignumber + + !$omp parallel do default(none) & + !$omp shared(npart, xyzh,metrics,metricderivs,vxyzu,dens,C_force,fext) & + !$omp firstprivate(pri) & + !$omp private(i,dtf) & + !$omp reduction(min:dtextforce_min) + do i=1, npart + call pack_metricderivs(xyzh(1:3,i),metricderivs(:,:,:,i)) + call get_grforce(xyzh(:,i),metrics(:,:,:,i),metricderivs(:,:,:,i), & + vxyzu(1:3,i),dens(i),vxyzu(4,i),pri,fext(1:3,i),dtf) + dtextforce_min = min(dtextforce_min,C_force*dtf) + enddo + !$omp end parallel do + ! manually add v contribution from gr + ! do i=1, npart + ! !fxyzu(:,i) = fxyzu(:,i) + fext(:,i) + ! vxyzu(1:3,i) = vxyzu(1:3,i) + fext(:,i)*dt_et + ! enddo +end subroutine get_metricderivs_all + +subroutine get_eos_quantities(densi,en) + use cons2prim, only:cons2primall + use part, only:dens,vxyzu,npart,metrics,xyzh,pxyzu,eos_vars + real, intent(out) :: densi,en + + !call h2dens(densi,xyzhi,metrici,vi) ! Compute dens from h + densi = dens(1) ! Feed the newly computed dens back out of the routine + !call cons2primall(npart,xyzh,metrics,vxyzu,dens,pxyzu,.true.) + call cons2primall(npart,xyzh,metrics,pxyzu,vxyzu,dens,eos_vars) + ! print*,"pxyzu: ",pxyzu(:,1) + ! print*, "vxyzu: ",vxyzu(:,1) + en = vxyzu(4,1) +end subroutine get_eos_quantities + + +end module einsteintk_wrapper diff --git a/src/utils/interpolate3D.F90 b/src/utils/interpolate3D.F90 new file mode 100644 index 000000000..5e1196284 --- /dev/null +++ b/src/utils/interpolate3D.F90 @@ -0,0 +1,1325 @@ +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.bitbucket.io/ ! +!--------------------------------------------------------------------------! +module interpolations3D +! +! interpolations3D +! +! :References: None +! +! :Owner: Spencer Magnall +! +! :Runtime parameters: None +! +! :Dependencies: einsteintk_utils, kernel +! + +!---------------------------------------------------------------------- +! +! Module containing all of the routines required for interpolation +! from 3D data to a 3D grid (SLOW!) +! +!---------------------------------------------------------------------- + + use einsteintk_utils, only:exact_rendering + use kernel, only:radkern2,radkern,cnormk,wkern!,wallint ! Moved to this module + !use interpolation, only:iroll ! Moved to this module + + !use timing, only:wall_time,print_time ! Using cpu_time for now + implicit none + integer, parameter :: doub_prec = kind(0.d0) + real :: cnormk3D = cnormk + public :: interpolate3D,interpolate3D_vecexact + +contains + !-------------------------------------------------------------------------- + ! subroutine to interpolate from particle data to even grid of pixels + ! + ! The data is interpolated according to the formula + ! + ! datsmooth(pixel) = sum_b weight_b dat_b W(r-r_b, h_b) + ! + ! where _b is the quantity at the neighbouring particle b and + ! W is the smoothing kernel, for which we use the usual cubic spline. + ! + ! For a standard SPH smoothing the weight function for each particle should be + ! + ! weight = pmass/(rho*h^3) + ! + ! this version is written for slices through a rectangular volume, ie. + ! assumes a uniform pixel size in x,y, whilst the number of pixels + ! in the z direction can be set to the number of cross-section slices. + ! + ! Input: particle coordinates : x,y,z (npart) + ! smoothing lengths : hh (npart) + ! weight for each particle : weight (npart) + ! scalar data to smooth : dat (npart) + ! + ! Output: smoothed data : datsmooth (npixx,npixy,npixz) + ! + ! Daniel Price, Institute of Astronomy, Cambridge 16/7/03 + ! Revised for "splash to grid", Monash University 02/11/09 + ! Maya Petkova contributed exact subgrid interpolation, April 2019 + !-------------------------------------------------------------------------- +subroutine interpolate3D(xyzh,weight,dat,itype,npart,& + xmin,ymin,zmin,datsmooth,npixx,npixy,npixz,pixwidthx,pixwidthy,pixwidthz,& + normalise,periodicx,periodicy,periodicz) + +integer, intent(in) :: npart,npixx,npixy,npixz +real, intent(in) :: xyzh(4,npart) +!real, intent(in), dimension(npart) :: x,y,z,hh ! change to xyzh() +real, intent(in), dimension(npart) :: weight,dat +integer, intent(in), dimension(npart) :: itype +real, intent(in) :: xmin,ymin,zmin,pixwidthx,pixwidthy,pixwidthz +real, intent(out), dimension(npixx,npixy,npixz) :: datsmooth +logical, intent(in) :: normalise,periodicx,periodicy,periodicz +!logical, intent(in), exact_rendering +real, allocatable :: datnorm(:,:,:) + +integer :: i,ipix,jpix,kpix +integer :: iprintinterval,iprintnext +integer :: ipixmin,ipixmax,jpixmin,jpixmax,kpixmin,kpixmax +integer :: ipixi,jpixi,kpixi,nxpix,nwarn,threadid +real :: xminpix,yminpix,zminpix,hmin !,dhmin3 +real, dimension(npixx) :: dx2i +real :: xi,yi,zi,hi,hi1,hi21,wab,q2,const,dyz2,dz2 +real :: term,termnorm,dy,dz,ypix,zpix,xpixi,pixwidthmax,dfac +real :: t_start,t_end,t_used +logical :: iprintprogress +real, dimension(npart) :: x,y,z,hh +real :: radkernel, radkernel2, radkernh + +! Exact rendering +real :: pixint, wint +!logical, parameter :: exact_rendering = .true. ! use exact rendering y/n +integer :: usedpart, negflag + + +!$ integer :: omp_get_num_threads,omp_get_thread_num +integer(kind=selected_int_kind(10)) :: iprogress,j ! up to 10 digits + +! Fill the particle data with xyzh +x(:) = xyzh(1,:) +y(:) = xyzh(2,:) +z(:) = xyzh(3,:) +hh(:) = xyzh(4,:) +!print*, "smoothing length: ", hh(1:10) +! cnormk3D set the value from the kernel routine +cnormk3D = cnormk +radkernel = radkern +radkernel2 = radkern2 +! print*, "radkern: ", radkern +! print*, "radkernel: ",radkernel +! print*, "radkern2: ", radkern2 + +! print*, "npix: ", npixx, npixy,npixz + +if (exact_rendering) then +print "(1x,a)",'interpolating to 3D grid (exact/Petkova+2018 on subgrid) ...' +elseif (normalise) then +print "(1x,a)",'interpolating to 3D grid (normalised) ...' +else +print "(1x,a)",'interpolating to 3D grid (non-normalised) ...' +endif +if (pixwidthx <= 0. .or. pixwidthy <= 0 .or. pixwidthz <= 0) then +print "(1x,a)",'interpolate3D: error: pixel width <= 0' +return +endif +if (any(hh(1:npart) <= tiny(hh))) then +print*,'interpolate3D: WARNING: ignoring some or all particles with h < 0' +endif + +!call wall_time(t_start) + +datsmooth = 0. +if (normalise) then +allocate(datnorm(npixx,npixy,npixz)) +datnorm = 0. +endif +! +!--print a progress report if it is going to take a long time +! (a "long time" is, however, somewhat system dependent) +! +iprintprogress = (npart >= 100000) .or. (npixx*npixy > 100000) !.or. exact_rendering +! +!--loop over particles +! +iprintinterval = 25 +if (npart >= 1e6) iprintinterval = 10 +iprintnext = iprintinterval +! +!--get starting CPU time +! +call cpu_time(t_start) + +usedpart = 0 + +xminpix = xmin !- 0.5*pixwidthx +yminpix = ymin !- 0.5*pixwidthy +zminpix = zmin !- 0.5*pixwidthz +! print*, "xminpix: ", xminpix +! print*, "yminpix: ", yminpix +! print*, "zminpix: ", zminpix +! print*, "dat: ", dat(1:10) +! print*, "weights: ", weight(1:10) +pixwidthmax = max(pixwidthx,pixwidthy,pixwidthz) +! +!--use a minimum smoothing length on the grid to make +! sure that particles contribute to at least one pixel +! +hmin = 0.5*pixwidthmax +!dhmin3 = 1./(hmin*hmin*hmin) + +const = cnormk3D ! normalisation constant (3D) +!print*, "const: ", const +nwarn = 0 +j = 0_8 +threadid = 1 +! +!--loop over particles +! +!$omp parallel default(none) & +!$omp shared(hh,z,x,y,weight,dat,itype,datsmooth,npart) & +!$omp shared(xmin,ymin,zmin,radkernel,radkernel2) & +!$omp shared(xminpix,yminpix,zminpix,pixwidthx,pixwidthy,pixwidthz) & +!$omp shared(npixx,npixy,npixz,const) & +!$omp shared(datnorm,normalise,periodicx,periodicy,periodicz,exact_rendering) & +!$omp shared(hmin,pixwidthmax) & +!$omp shared(iprintprogress,iprintinterval,j) & +!$omp private(hi,xi,yi,zi,radkernh,hi1,hi21) & +!$omp private(term,termnorm,xpixi,iprogress) & +!$omp private(ipixmin,ipixmax,jpixmin,jpixmax,kpixmin,kpixmax) & +!$omp private(ipix,jpix,kpix,ipixi,jpixi,kpixi) & +!$omp private(dx2i,nxpix,zpix,dz,dz2,dyz2,dy,ypix,q2,wab) & +!$omp private(pixint,wint,negflag,dfac,threadid) & +!$omp firstprivate(iprintnext) & +!$omp reduction(+:nwarn,usedpart) +!$omp master +!$ print "(1x,a,i3,a)",'Using ',omp_get_num_threads(),' cpus' +!$omp end master + +!$omp do schedule (guided, 2) +over_parts: do i=1,npart +! +!--report on progress +! +if (iprintprogress) then + !$omp atomic + j=j+1_8 +!$ threadid = omp_get_thread_num() + iprogress = 100*j/npart + if (iprogress >= iprintnext .and. threadid==1) then + write(*,"(i3,'%.')",advance='no') iprogress + iprintnext = iprintnext + iprintinterval + endif +endif +! +!--skip particles with itype < 0 +! +if (itype(i) < 0 .or. weight(i) < tiny(0.)) cycle over_parts + +hi = hh(i) +if (hi <= 0.) then + cycle over_parts +elseif (hi < hmin) then + ! + !--use minimum h to capture subgrid particles + ! (get better results *without* adjusting weights) + ! + termnorm = const*weight(i) !*(hi*hi*hi)*dhmin3 + if (.not.exact_rendering) hi = hmin +else + termnorm = const*weight(i) +endif + +! +!--set kernel related quantities +! +xi = x(i) +yi = y(i) +zi = z(i) + +hi1 = 1./hi +hi21 = hi1*hi1 +radkernh = radkernel*hi ! radius of the smoothing kernel +!termnorm = const*weight(i) +term = termnorm*dat(i) +dfac = hi**3/(pixwidthx*pixwidthy*pixwidthz*const) +!dfac = hi**3/(pixwidthx*pixwidthy*const) +! +!--for each particle work out which pixels it contributes to +! +ipixmin = int((xi - radkernh - xmin)/pixwidthx) +jpixmin = int((yi - radkernh - ymin)/pixwidthy) +kpixmin = int((zi - radkernh - zmin)/pixwidthz) +ipixmax = int((xi + radkernh - xmin)/pixwidthx) + 1 +jpixmax = int((yi + radkernh - ymin)/pixwidthy) + 1 +kpixmax = int((zi + radkernh - zmin)/pixwidthz) + 1 + +if (.not.periodicx) then + if (ipixmin < 1) ipixmin = 1 ! make sure they only contribute + if (ipixmax > npixx) ipixmax = npixx ! to pixels in the image +endif +if (.not.periodicy) then + if (jpixmin < 1) jpixmin = 1 + if (jpixmax > npixy) jpixmax = npixy +endif +if (.not.periodicz) then + if (kpixmin < 1) kpixmin = 1 + if (kpixmax > npixz) kpixmax = npixz +endif + +negflag = 0 + +! +!--precalculate an array of dx2 for this particle (optimisation) +! +! Check the x position of the grid cells +!open(unit=677,file="posxgrid.txt",action='write',position='append') +nxpix = 0 +do ipix=ipixmin,ipixmax + nxpix = nxpix + 1 + ipixi = ipix + if (periodicx) ipixi = iroll(ipix,npixx) + xpixi = xminpix + ipix*pixwidthx + !write(677,*) ipix, xpixi + !--watch out for errors with periodic wrapping... + if (nxpix <= size(dx2i)) then + dx2i(nxpix) = ((xpixi - xi)**2)*hi21 + endif +enddo + +!--if particle contributes to more than npixx pixels +! (i.e. periodic boundaries wrap more than once) +! truncate the contribution and give warning +if (nxpix > npixx) then + nwarn = nwarn + 1 + ipixmax = ipixmin + npixx - 1 +endif +! +!--loop over pixels, adding the contribution from this particle +! +do kpix = kpixmin,kpixmax + kpixi = kpix + if (periodicz) kpixi = iroll(kpix,npixz) + + zpix = zminpix + kpix*pixwidthz + dz = zpix - zi + dz2 = dz*dz*hi21 + + do jpix = jpixmin,jpixmax + jpixi = jpix + if (periodicy) jpixi = iroll(jpix,npixy) + + ypix = yminpix + jpix*pixwidthy + dy = ypix - yi + dyz2 = dy*dy*hi21 + dz2 + + nxpix = 0 + do ipix = ipixmin,ipixmax + if ((kpix==kpixmin).and.(jpix==jpixmin).and.(ipix==ipixmin)) then + usedpart = usedpart + 1 + endif + + nxpix = nxpix + 1 + ipixi = ipix + if (periodicx) ipixi = iroll(ipix,npixx) + + q2 = dx2i(nxpix) + dyz2 ! dx2 pre-calculated; dy2 pre-multiplied by hi21 + + if (exact_rendering .and. ipixmax-ipixmin <= 4) then + if (q2 < radkernel2 + 3.*pixwidthmax**2*hi21) then + xpixi = xminpix + ipix*pixwidthx + + ! Contribution of the cell walls in the xy-plane + pixint = 0.0 + wint = wallint(zpix-zi+0.5*pixwidthz,xi,yi,xpixi,ypix,pixwidthx,pixwidthy,hi) + pixint = pixint + wint + + wint = wallint(zi-zpix+0.5*pixwidthz,xi,yi,xpixi,ypix,pixwidthx,pixwidthy,hi) + pixint = pixint + wint + + ! Contribution of the cell walls in the xz-plane + wint = wallint(ypix-yi+0.5*pixwidthy,xi,zi,xpixi,zpix,pixwidthx,pixwidthz,hi) + pixint = pixint + wint + + wint = wallint(yi-ypix+0.5*pixwidthy,xi,zi,xpixi,zpix,pixwidthx,pixwidthz,hi) + pixint = pixint + wint + + ! Contribution of the cell walls in the yz-plane + wint = wallint(xpixi-xi+0.5*pixwidthx,zi,yi,zpix,ypix,pixwidthz,pixwidthy,hi) + pixint = pixint + wint + + wint = wallint(xi-xpixi+0.5*pixwidthx,zi,yi,zpix,ypix,pixwidthz,pixwidthy,hi) + pixint = pixint + wint + + wab = pixint*dfac ! /(pixwidthx*pixwidthy*pixwidthz*const)*hi**3 + + if (pixint < -0.01d0) then + print*, "Error: (",ipixi,jpixi,kpixi,") -> ", pixint, term*wab + endif + + ! + !--calculate data value at this pixel using the summation interpolant + ! + !$omp atomic + datsmooth(ipixi,jpixi,kpixi) = datsmooth(ipixi,jpixi,kpixi) + term*wab + if (normalise) then + !$omp atomic + datnorm(ipixi,jpixi,kpixi) = datnorm(ipixi,jpixi,kpixi) + termnorm*wab + endif + endif + else + if (q2 < radkernel2) then + + ! + !--SPH kernel - standard cubic spline + ! + wab = wkernel(q2) + ! + !--calculate data value at this pixel using the summation interpolant + ! + !$omp atomic + datsmooth(ipixi,jpixi,kpixi) = datsmooth(ipixi,jpixi,kpixi) + term*wab + if (normalise) then + !$omp atomic + datnorm(ipixi,jpixi,kpixi) = datnorm(ipixi,jpixi,kpixi) + termnorm*wab + endif + endif + endif + enddo + enddo +enddo +enddo over_parts +!$omp enddo +!$omp end parallel + +if (nwarn > 0) then +print "(a,i11,a,/,a)",' interpolate3D: WARNING: contributions truncated from ',nwarn,' particles',& + ' that wrap periodic boundaries more than once' +endif +! +!--normalise dat array +! +if (normalise) then +where (datnorm > tiny(datnorm)) + datsmooth = datsmooth/datnorm +end where +endif +if (allocated(datnorm)) deallocate(datnorm) + +!call wall_time(t_end) +call cpu_time(t_end) +t_used = t_end - t_start +print*, 'Interpolate3D completed in ',t_end-t_start,'s' +!if (t_used > 10.) call print_time(t_used) + +!print*, 'Number of particles in the volume: ', usedpart +! datsmooth(1,1,1) = 3.14159 +! datsmooth(32,32,32) = 3.145159 +! datsmooth(11,11,11) = 3.14159 +! datsmooth(10,10,10) = 3.145159 + +end subroutine interpolate3D + +subroutine interpolate3D_vecexact(xyzh,weight,dat,ilendat,itype,npart,& + xmin,ymin,zmin,datsmooth,npixx,npixy,npixz,pixwidthx,pixwidthy,pixwidthz,& + normalise,periodicx,periodicy,periodicz) + + integer, intent(in) :: npart,npixx,npixy,npixz,ilendat + real, intent(in) :: xyzh(4,npart) + !real, intent(in), dimension(npart) :: x,y,z,hh ! change to xyzh() + real, intent(in), dimension(npart) :: weight + real, intent(in),dimension(npart,ilendat) :: dat + integer, intent(in), dimension(npart) :: itype + real, intent(in) :: xmin,ymin,zmin,pixwidthx,pixwidthy,pixwidthz + real, intent(out), dimension(ilendat,npixx,npixy,npixz) :: datsmooth + logical, intent(in) :: normalise,periodicx,periodicy,periodicz + !logical, intent(in), exact_rendering + real, allocatable :: datnorm(:,:,:) + + integer :: i,ipix,jpix,kpix,lockindex,smoothindex + integer :: iprintinterval,iprintnext + integer :: ipixmin,ipixmax,jpixmin,jpixmax,kpixmin,kpixmax + integer :: ipixi,jpixi,kpixi,nxpix,nwarn,threadid + real :: xminpix,yminpix,zminpix,hmin !,dhmin3 + real, dimension(npixx) :: dx2i + real :: xi,yi,zi,hi,hi1,hi21,wab,q2,const,dyz2,dz2 + real :: term(ilendat),termnorm,dy,dz,ypix,zpix,xpixi,pixwidthmax,dfac + real :: t_start,t_end,t_used + logical :: iprintprogress + real, dimension(npart) :: x,y,z,hh + real :: radkernel, radkernel2, radkernh + + ! Exact rendering + real :: pixint, wint + !logical, parameter :: exact_rendering = .true. ! use exact rendering y/n + integer :: usedpart, negflag + + +!$ integer :: omp_get_num_threads,omp_get_thread_num + integer(kind=selected_int_kind(10)) :: iprogress,j ! up to 10 digits + + ! Fill the particle data with xyzh + x(:) = xyzh(1,:) + y(:) = xyzh(2,:) + z(:) = xyzh(3,:) + hh(:) = xyzh(4,:) + !print*, "smoothing length: ", hh(1:10) + ! cnormk3D set the value from the kernel routine + cnormk3D = cnormk + radkernel = radkern + radkernel2 = radkern2 +! print*, "radkern: ", radkern +! print*, "radkernel: ",radkernel +! print*, "radkern2: ", radkern2 + + !print*, "npix: ", npixx, npixy,npixz + + if (exact_rendering) then + print "(1x,a)",'interpolating to 3D grid (exact/Petkova+2018 on subgrid) ...' + elseif (normalise) then + print "(1x,a)",'interpolating to 3D grid (normalised) ...' + else + print "(1x,a)",'interpolating to 3D grid (non-normalised) ...' + endif + if (pixwidthx <= 0. .or. pixwidthy <= 0 .or. pixwidthz <= 0) then + print "(1x,a)",'interpolate3D: error: pixel width <= 0' + return + endif + if (any(hh(1:npart) <= tiny(hh))) then + print*,'interpolate3D: WARNING: ignoring some or all particles with h < 0' + endif + + !call wall_time(t_start) + +!! $ allocate(ilock(npixx*npixy*npixz)) +!! $ do i=1,npixx*npixy*npixz +!! $ call omp_init_lock(ilock(i)) +!! $ enddo + + datsmooth = 0. + if (normalise) then + allocate(datnorm(npixx,npixy,npixz)) + datnorm = 0. + endif + ! + !--print a progress report if it is going to take a long time + ! (a "long time" is, however, somewhat system dependent) + ! + iprintprogress = (npart >= 100000) .or. (npixx*npixy > 100000) !.or. exact_rendering + ! + !--loop over particles + ! + iprintinterval = 25 + if (npart >= 1e6) iprintinterval = 10 + iprintnext = iprintinterval + ! + !--get starting CPU time + ! + call cpu_time(t_start) + + usedpart = 0 + + xminpix = xmin !- 0.5*pixwidthx + yminpix = ymin !- 0.5*pixwidthy + zminpix = zmin !- 0.5*pixwidthz +! print*, "xminpix: ", xminpix +! print*, "yminpix: ", yminpix +! print*, "zminpix: ", zminpix +! print*, "dat: ", dat(1:10) +! print*, "weights: ", weight(1:10) + pixwidthmax = max(pixwidthx,pixwidthy,pixwidthz) + ! + !--use a minimum smoothing length on the grid to make + ! sure that particles contribute to at least one pixel + ! + hmin = 0.5*pixwidthmax + !dhmin3 = 1./(hmin*hmin*hmin) + + const = cnormk3D ! normalisation constant (3D) + !print*, "const: ", const + nwarn = 0 + j = 0_8 + threadid = 1 + ! + !--loop over particles + ! + !$omp parallel default(none) & + !$omp shared(hh,z,x,y,weight,dat,itype,datsmooth,npart) & + !$omp shared(xmin,ymin,zmin,radkernel,radkernel2) & + !$omp shared(xminpix,yminpix,zminpix,pixwidthx,pixwidthy,pixwidthz) & + !$omp shared(npixx,npixy,npixz,const,ilendat) & + !$omp shared(datnorm,normalise,periodicx,periodicy,periodicz,exact_rendering) & + !$omp shared(hmin,pixwidthmax) & + !$omp shared(iprintprogress,iprintinterval,j) & + !$omp private(hi,xi,yi,zi,radkernh,hi1,hi21) & + !$omp private(term,termnorm,xpixi,iprogress) & + !$omp private(ipixmin,ipixmax,jpixmin,jpixmax,kpixmin,kpixmax) & + !$omp private(ipix,jpix,kpix,ipixi,jpixi,kpixi) & + !$omp private(dx2i,nxpix,zpix,dz,dz2,dyz2,dy,ypix,q2,wab) & + !$omp private(pixint,wint,negflag,dfac,threadid,lockindex,smoothindex) & + !$omp firstprivate(iprintnext) & + !$omp reduction(+:nwarn,usedpart) + !$omp master +!$ print "(1x,a,i3,a)",'Using ',omp_get_num_threads(),' cpus' + !$omp end master + + !$omp do schedule (guided, 2) + over_parts: do i=1,npart + ! + !--report on progress + ! + if (iprintprogress) then + !$omp atomic + j=j+1_8 +!$ threadid = omp_get_thread_num() + iprogress = 100*j/npart + if (iprogress >= iprintnext .and. threadid==1) then + write(*,"(i3,'%.')",advance='no') iprogress + iprintnext = iprintnext + iprintinterval + endif + endif + ! + !--skip particles with itype < 0 + ! + if (itype(i) < 0 .or. weight(i) < tiny(0.)) cycle over_parts + + hi = hh(i) + if (hi <= 0.) then + cycle over_parts + elseif (hi < hmin) then + ! + !--use minimum h to capture subgrid particles + ! (get better results *without* adjusting weights) + ! + termnorm = const*weight(i) !*(hi*hi*hi)*dhmin3 + if (.not.exact_rendering) hi = hmin + else + termnorm = const*weight(i) + endif + + ! + !--set kernel related quantities + ! + xi = x(i) + yi = y(i) + zi = z(i) + + hi1 = 1./hi + hi21 = hi1*hi1 + radkernh = radkernel*hi ! radius of the smoothing kernel + !termnorm = const*weight(i) + term(:) = termnorm*dat(i,:) + dfac = hi**3/(pixwidthx*pixwidthy*pixwidthz*const) + !dfac = hi**3/(pixwidthx*pixwidthy*const) + ! + !--for each particle work out which pixels it contributes to + ! + ipixmin = int((xi - radkernh - xmin)/pixwidthx) + jpixmin = int((yi - radkernh - ymin)/pixwidthy) + kpixmin = int((zi - radkernh - zmin)/pixwidthz) + ipixmax = int((xi + radkernh - xmin)/pixwidthx) + 1 + jpixmax = int((yi + radkernh - ymin)/pixwidthy) + 1 + kpixmax = int((zi + radkernh - zmin)/pixwidthz) + 1 + + if (.not.periodicx) then + if (ipixmin < 1) ipixmin = 1 ! make sure they only contribute + if (ipixmax > npixx) ipixmax = npixx ! to pixels in the image + endif + if (.not.periodicy) then + if (jpixmin < 1) jpixmin = 1 + if (jpixmax > npixy) jpixmax = npixy + endif + if (.not.periodicz) then + if (kpixmin < 1) kpixmin = 1 + if (kpixmax > npixz) kpixmax = npixz + endif + + negflag = 0 + + ! + !--precalculate an array of dx2 for this particle (optimisation) + ! + ! Check the x position of the grid cells + !open(unit=677,file="posxgrid.txt",action='write',position='append') + nxpix = 0 + do ipix=ipixmin,ipixmax + nxpix = nxpix + 1 + ipixi = ipix + if (periodicx) ipixi = iroll(ipix,npixx) + xpixi = xminpix + ipix*pixwidthx + !write(677,*) ipix, xpixi + !--watch out for errors with periodic wrapping... + if (nxpix <= size(dx2i)) then + dx2i(nxpix) = ((xpixi - xi)**2)*hi21 + endif + enddo + + !--if particle contributes to more than npixx pixels + ! (i.e. periodic boundaries wrap more than once) + ! truncate the contribution and give warning + if (nxpix > npixx) then + nwarn = nwarn + 1 + ipixmax = ipixmin + npixx - 1 + endif + ! + !--loop over pixels, adding the contribution from this particle + ! + do kpix = kpixmin,kpixmax + kpixi = kpix + if (periodicz) kpixi = iroll(kpix,npixz) + + zpix = zminpix + kpix*pixwidthz + dz = zpix - zi + dz2 = dz*dz*hi21 + + do jpix = jpixmin,jpixmax + jpixi = jpix + if (periodicy) jpixi = iroll(jpix,npixy) + + ypix = yminpix + jpix*pixwidthy + dy = ypix - yi + dyz2 = dy*dy*hi21 + dz2 + + nxpix = 0 + do ipix = ipixmin,ipixmax + if ((kpix==kpixmin).and.(jpix==jpixmin).and.(ipix==ipixmin)) then + usedpart = usedpart + 1 + endif + + nxpix = nxpix + 1 + ipixi = ipix + if (periodicx) ipixi = iroll(ipix,npixx) + + q2 = dx2i(nxpix) + dyz2 ! dx2 pre-calculated; dy2 pre-multiplied by hi21 + + if (exact_rendering .and. ipixmax-ipixmin <= 4) then + if (q2 < radkernel2 + 3.*pixwidthmax**2*hi21) then + xpixi = xminpix + ipix*pixwidthx + + ! Contribution of the cell walls in the xy-plane + pixint = 0.0 + wint = wallint(zpix-zi+0.5*pixwidthz,xi,yi,xpixi,ypix,pixwidthx,pixwidthy,hi) + pixint = pixint + wint + + wint = wallint(zi-zpix+0.5*pixwidthz,xi,yi,xpixi,ypix,pixwidthx,pixwidthy,hi) + pixint = pixint + wint + + ! Contribution of the cell walls in the xz-plane + wint = wallint(ypix-yi+0.5*pixwidthy,xi,zi,xpixi,zpix,pixwidthx,pixwidthz,hi) + pixint = pixint + wint + + wint = wallint(yi-ypix+0.5*pixwidthy,xi,zi,xpixi,zpix,pixwidthx,pixwidthz,hi) + pixint = pixint + wint + + ! Contribution of the cell walls in the yz-plane + wint = wallint(xpixi-xi+0.5*pixwidthx,zi,yi,zpix,ypix,pixwidthz,pixwidthy,hi) + pixint = pixint + wint + + wint = wallint(xi-xpixi+0.5*pixwidthx,zi,yi,zpix,ypix,pixwidthz,pixwidthy,hi) + pixint = pixint + wint + + wab = pixint*dfac ! /(pixwidthx*pixwidthy*pixwidthz*const)*hi**3 + + if (pixint < -0.01d0) then + print*, "Error: (",ipixi,jpixi,kpixi,") -> ", pixint, term*wab + endif + + ! + !--calculate data value at this pixel using the summation interpolant + ! + ! Find out where this pixel sits in the lock array + ! lockindex = (k-1)*nx*ny + (j-1)*nx + i + !lockindex = (kpixi-1)*npixx*npixy + (jpixi-1)*npixx + ipixi + !!$call omp_set_lock(ilock(lockindex)) + !!$omp critical (datsmooth) + do smoothindex=1, ilendat + !$omp atomic + datsmooth(smoothindex,ipixi,jpixi,kpixi) = datsmooth(smoothindex,ipixi,jpixi,kpixi) + term(smoothindex)*wab + enddo + !!$omp end critical (datsmooth) + if (normalise) then + !$omp atomic + !!$omp critical (datnorm) + datnorm(ipixi,jpixi,kpixi) = datnorm(ipixi,jpixi,kpixi) + termnorm*wab + !!$omp end critical (datnorm) + endif + + !!$call omp_unset_lock(ilock(lockindex)) + endif + else + if (q2 < radkernel2) then + + ! + !--SPH kernel - standard cubic spline + ! + wab = wkernel(q2) + ! + !--calculate data value at this pixel using the summation interpolant + ! + !!$omp atomic ! Atomic statmements only work with scalars + !!$omp set lock ! Does this work with an array? + ! Find out where this pixel sits in the lock array + ! lockindex = (k-1)*nx*ny + (j-1)*nx + i + !lockindex = (kpixi-1)*npixx*npixy + (jpixi-1)*npixx + ipixi + !!$call omp_set_lock(ilock(lockindex)) + !!$omp critical (datsmooth) + do smoothindex=1,ilendat + !$omp atomic + datsmooth(smoothindex,ipixi,jpixi,kpixi) = datsmooth(smoothindex,ipixi,jpixi,kpixi) + term(smoothindex)*wab + enddo + !!$omp end critical (datsmooth) + if (normalise) then + !$omp atomic + !!$omp critical (datnorm) + datnorm(ipixi,jpixi,kpixi) = datnorm(ipixi,jpixi,kpixi) + termnorm*wab + !!$omp end critical (datnorm) + endif + !!$call omp_unset_lock(ilock(lockindex)) + + endif + endif + enddo + enddo + enddo + enddo over_parts + !$omp enddo + !$omp end parallel + +!!$ do i=1,npixx*npixy*npixz +!!$ call omp_destroy_lock(ilock(i)) +!!$ enddo +!!$ if (allocated(ilock)) deallocate(ilock) + + if (nwarn > 0) then + print "(a,i11,a,/,a)",' interpolate3D: WARNING: contributions truncated from ',nwarn,' particles',& + ' that wrap periodic boundaries more than once' + endif + ! + !--normalise dat array + ! + if (normalise) then + do i=1, ilendat + where (datnorm > tiny(datnorm)) + + datsmooth(i,:,:,:) = datsmooth(i,:,:,:)/datnorm(:,:,:) + end where + enddo + endif + if (allocated(datnorm)) deallocate(datnorm) + + !call wall_time(t_end) + call cpu_time(t_end) + t_used = t_end - t_start + print*, 'Interpolate3DVec completed in ',t_end-t_start,'s' + !if (t_used > 10.) call print_time(t_used) + + !print*, 'Number of particles in the volume: ', usedpart + ! datsmooth(1,1,1) = 3.14159 + ! datsmooth(32,32,32) = 3.145159 + ! datsmooth(11,11,11) = 3.14159 + ! datsmooth(10,10,10) = 3.145159 + +end subroutine interpolate3D_vecexact + + ! subroutine interpolate3D_vec(x,y,z,hh,weight,datvec,itype,npart,& + ! xmin,ymin,zmin,datsmooth,npixx,npixy,npixz,pixwidthx,pixwidthy,pixwidthz,& + ! normalise,periodicx,periodicy,periodicz) + + ! integer, intent(in) :: npart,npixx,npixy,npixz + ! real, intent(in), dimension(npart) :: x,y,z,hh,weight + ! real, intent(in), dimension(npart,3) :: datvec + ! integer, intent(in), dimension(npart) :: itype + ! real, intent(in) :: xmin,ymin,zmin,pixwidthx,pixwidthy,pixwidthz + ! real(doub_prec), intent(out), dimension(3,npixx,npixy,npixz) :: datsmooth + ! logical, intent(in) :: normalise,periodicx,periodicy,periodicz + ! real(doub_prec), dimension(npixx,npixy,npixz) :: datnorm + + ! integer :: i,ipix,jpix,kpix + ! integer :: iprintinterval,iprintnext + ! integer :: ipixmin,ipixmax,jpixmin,jpixmax,kpixmin,kpixmax + ! integer :: ipixi,jpixi,kpixi,nxpix,nwarn + ! real :: xminpix,yminpix,zminpix + ! real, dimension(npixx) :: dx2i + ! real :: xi,yi,zi,hi,hi1,hi21,radkern,wab,q2,const,dyz2,dz2 + ! real :: termnorm,dy,dz,ypix,zpix,xpixi,ddatnorm + ! real, dimension(3) :: term + ! !real :: t_start,t_end + ! logical :: iprintprogress + ! !$ integer :: omp_get_num_threads + ! integer(kind=selected_int_kind(10)) :: iprogress ! up to 10 digits + + ! datsmooth = 0. + ! datnorm = 0. + ! if (normalise) then + ! print "(1x,a)",'interpolating to 3D grid (normalised) ...' + ! else + ! print "(1x,a)",'interpolating to 3D grid (non-normalised) ...' + ! endif + ! if (pixwidthx <= 0. .or. pixwidthy <= 0. .or. pixwidthz <= 0.) then + ! print "(1x,a)",'interpolate3D: error: pixel width <= 0' + ! return + ! endif + ! if (any(hh(1:npart) <= tiny(hh))) then + ! print*,'interpolate3D: WARNING: ignoring some or all particles with h < 0' + ! endif + + ! ! + ! !--print a progress report if it is going to take a long time + ! ! (a "long time" is, however, somewhat system dependent) + ! ! + ! iprintprogress = (npart >= 100000) .or. (npixx*npixy > 100000) + ! !$ iprintprogress = .false. + ! ! + ! !--loop over particles + ! ! + ! iprintinterval = 25 + ! if (npart >= 1e6) iprintinterval = 10 + ! iprintnext = iprintinterval + ! ! + ! !--get starting CPU time + ! ! + ! !call cpu_time(t_start) + + ! xminpix = xmin - 0.5*pixwidthx + ! yminpix = ymin - 0.5*pixwidthy + ! zminpix = zmin - 0.5*pixwidthz + + ! const = cnormk3D ! normalisation constant (3D) + ! nwarn = 0 + + ! !$omp parallel default(none) & + ! !$omp shared(hh,z,x,y,weight,datvec,itype,datsmooth,npart) & + ! !$omp shared(xmin,ymin,zmin,radkernel,radkernel2) & + ! !$omp shared(xminpix,yminpix,zminpix,pixwidthx,pixwidthy,pixwidthz) & + ! !$omp shared(npixx,npixy,npixz,const) & + ! !$omp shared(iprintprogress,iprintinterval) & + ! !$omp shared(datnorm,normalise,periodicx,periodicy,periodicz) & + ! !$omp private(hi,xi,yi,zi,radkern,hi1,hi21) & + ! !$omp private(term,termnorm,xpixi) & + ! !$omp private(iprogress,iprintnext) & + ! !$omp private(ipixmin,ipixmax,jpixmin,jpixmax,kpixmin,kpixmax) & + ! !$omp private(ipix,jpix,kpix,ipixi,jpixi,kpixi) & + ! !$omp private(dx2i,nxpix,zpix,dz,dz2,dyz2,dy,ypix,q2,wab) & + ! !$omp reduction(+:nwarn) + ! !$omp master + ! !$ print "(1x,a,i3,a)",'Using ',omp_get_num_threads(),' cpus' + ! !$omp end master + ! ! + ! !--loop over particles + ! ! + ! !$omp do schedule (guided, 2) + ! over_parts: do i=1,npart + ! ! + ! !--report on progress + ! ! + ! if (iprintprogress) then + ! iprogress = 100*i/npart + ! if (iprogress >= iprintnext) then + ! write(*,"('(',i3,'% -',i12,' particles done)')") iprogress,i + ! iprintnext = iprintnext + iprintinterval + ! endif + ! endif + ! ! + ! !--skip particles with itype < 0 + ! ! + ! if (itype(i) < 0 .or. weight(i) < tiny(0.)) cycle over_parts + + ! hi = hh(i) + ! if (hi <= 0.) cycle over_parts + + ! ! + ! !--set kernel related quantities + ! ! + ! xi = x(i) + ! yi = y(i) + ! zi = z(i) + + ! hi1 = 1./hi + ! hi21 = hi1*hi1 + ! radkern = radkernel*hi ! radius of the smoothing kernel + ! termnorm = const*weight(i) + ! term(:) = termnorm*datvec(i,:) + ! ! + ! !--for each particle work out which pixels it contributes to + ! ! + ! ipixmin = int((xi - radkern - xmin)/pixwidthx) + ! jpixmin = int((yi - radkern - ymin)/pixwidthy) + ! kpixmin = int((zi - radkern - zmin)/pixwidthz) + ! ipixmax = int((xi + radkern - xmin)/pixwidthx) + 1 + ! jpixmax = int((yi + radkern - ymin)/pixwidthy) + 1 + ! kpixmax = int((zi + radkern - zmin)/pixwidthz) + 1 + + ! if (.not.periodicx) then + ! if (ipixmin < 1) ipixmin = 1 ! make sure they only contribute + ! if (ipixmax > npixx) ipixmax = npixx ! to pixels in the image + ! endif + ! if (.not.periodicy) then + ! if (jpixmin < 1) jpixmin = 1 + ! if (jpixmax > npixy) jpixmax = npixy + ! endif + ! if (.not.periodicz) then + ! if (kpixmin < 1) kpixmin = 1 + ! if (kpixmax > npixz) kpixmax = npixz + ! endif + ! ! + ! !--precalculate an array of dx2 for this particle (optimisation) + ! ! + ! nxpix = 0 + ! do ipix=ipixmin,ipixmax + ! nxpix = nxpix + 1 + ! ipixi = ipix + ! if (periodicx) ipixi = iroll(ipix,npixx) + ! xpixi = xminpix + ipix*pixwidthx + ! !--watch out for errors with perioic wrapping... + ! if (nxpix <= size(dx2i)) then + ! dx2i(nxpix) = ((xpixi - xi)**2)*hi21 + ! endif + ! enddo + + ! !--if particle contributes to more than npixx pixels + ! ! (i.e. periodic boundaries wrap more than once) + ! ! truncate the contribution and give warning + ! if (nxpix > npixx) then + ! nwarn = nwarn + 1 + ! ipixmax = ipixmin + npixx - 1 + ! endif + ! ! + ! !--loop over pixels, adding the contribution from this particle + ! ! + ! do kpix = kpixmin,kpixmax + ! kpixi = kpix + ! if (periodicz) kpixi = iroll(kpix,npixz) + ! zpix = zminpix + kpix*pixwidthz + ! dz = zpix - zi + ! dz2 = dz*dz*hi21 + + ! do jpix = jpixmin,jpixmax + ! jpixi = jpix + ! if (periodicy) jpixi = iroll(jpix,npixy) + ! ypix = yminpix + jpix*pixwidthy + ! dy = ypix - yi + ! dyz2 = dy*dy*hi21 + dz2 + + ! nxpix = 0 + ! do ipix = ipixmin,ipixmax + ! ipixi = ipix + ! if (periodicx) ipixi = iroll(ipix,npixx) + ! nxpix = nxpix + 1 + ! q2 = dx2i(nxpix) + dyz2 ! dx2 pre-calculated; dy2 pre-multiplied by hi21 + ! ! + ! !--SPH kernel - standard cubic spline + ! ! + ! if (q2 < radkernel2) then + ! wab = wkernel(q2) + ! ! + ! !--calculate data value at this pixel using the summation interpolant + ! ! + ! !$omp atomic + ! datsmooth(1,ipixi,jpixi,kpixi) = datsmooth(1,ipixi,jpixi,kpixi) + term(1)*wab + ! !$omp atomic + ! datsmooth(2,ipixi,jpixi,kpixi) = datsmooth(2,ipixi,jpixi,kpixi) + term(2)*wab + ! !$omp atomic + ! datsmooth(3,ipixi,jpixi,kpixi) = datsmooth(3,ipixi,jpixi,kpixi) + term(3)*wab + ! if (normalise) then + ! !$omp atomic + ! datnorm(ipixi,jpixi,kpixi) = datnorm(ipixi,jpixi,kpixi) + termnorm*wab + ! endif + ! endif + ! enddo + ! enddo + ! enddo + ! enddo over_parts + ! !$omp enddo + ! !$omp end parallel + + ! if (nwarn > 0) then + ! print "(a,i11,a,/,a)",' interpolate3D: WARNING: contributions truncated from ',nwarn,' particles',& + ! ' that wrap periodic boundaries more than once' + ! endif + ! ! + ! !--normalise dat array + ! ! + ! if (normalise) then + ! !$omp parallel do default(none) schedule(static) & + ! !$omp shared(datsmooth,datnorm,npixz,npixy,npixx) & + ! !$omp private(kpix,jpix,ipix,ddatnorm) + ! do kpix=1,npixz + ! do jpix=1,npixy + ! do ipix=1,npixx + ! if (datnorm(ipix,jpix,kpix) > tiny(datnorm)) then + ! ddatnorm = 1./datnorm(ipix,jpix,kpix) + ! datsmooth(1,ipix,jpix,kpix) = datsmooth(1,ipix,jpix,kpix)*ddatnorm + ! datsmooth(2,ipix,jpix,kpix) = datsmooth(2,ipix,jpix,kpix)*ddatnorm + ! datsmooth(3,ipix,jpix,kpix) = datsmooth(3,ipix,jpix,kpix)*ddatnorm + ! endif + ! enddo + ! enddo + ! enddo + ! !$omp end parallel do + ! endif + + ! return + + ! end subroutine interpolate3D_vec + + !------------------------------------------------------------ + ! interface to kernel routine to avoid problems with openMP + !----------------------------------------------------------- +real function wkernel(q2) + use kernel, only:wkern + real, intent(in) :: q2 + real :: q + q = sqrt(q2) + wkernel = wkern(q2,q) + +end function wkernel + + !------------------------------------------------------------ + ! 3D functions to evaluate exact overlap of kernel with wall boundaries + ! see Petkova, Laibe & Bonnell (2018), J. Comp. Phys + !------------------------------------------------------------ +real function wallint(r0, xp, yp, xc, yc, pixwidthx, pixwidthy, hi) + real, intent(in) :: r0, xp, yp, xc, yc, pixwidthx, pixwidthy, hi + real(doub_prec) :: R_0, d1, d2, dx, dy, h + + wallint = 0.0 + dx = xc - xp + dy = yc - yp + h = hi + + ! + ! Contributions from each of the 4 sides of a cell wall + ! + R_0 = 0.5*pixwidthy + dy + d1 = 0.5*pixwidthx - dx + d2 = 0.5*pixwidthx + dx + wallint = wallint + pint3D(r0, R_0, d1, d2, h) + + R_0 = 0.5*pixwidthy - dy + d1 = 0.5*pixwidthx + dx + d2 = 0.5*pixwidthx - dx + wallint = wallint + pint3D(r0, R_0, d1, d2, h) + + R_0 = 0.5*pixwidthx + dx + d1 = 0.5*pixwidthy + dy + d2 = 0.5*pixwidthy - dy + wallint = wallint + pint3D(r0, R_0, d1, d2, h) + + R_0 = 0.5*pixwidthx - dx + d1 = 0.5*pixwidthy - dy + d2 = 0.5*pixwidthy + dy + wallint = wallint + pint3D(r0, R_0, d1, d2, h) + +end function wallint + + +real function pint3D(r0, R_0, d1, d2, hi) + + real(doub_prec), intent(in) :: R_0, d1, d2, hi + real, intent(in) :: r0 + real(doub_prec) :: ar0, aR_0 + real(doub_prec) :: int1, int2 + !integer :: fflag = 0 + + if (abs(r0) < tiny(0.)) then + pint3D = 0.d0 + return + endif + + if (r0 > 0.d0) then + pint3D = 1.d0 + ar0 = r0 + else + pint3D = -1.d0 + ar0 = -r0 + endif + + if (R_0 > 0.d0) then + aR_0 = R_0 + else + pint3D = -pint3D + aR_0 = -R_0 + endif + + int1 = full_integral_3D(d1, ar0, aR_0, hi) + int2 = full_integral_3D(d2, ar0, aR_0, hi) + + if (int1 < 0.d0) int1 = 0.d0 + if (int2 < 0.d0) int2 = 0.d0 + + if (d1*d2 >= 0) then + pint3D = pint3D*(int1 + int2) + if (int1 + int2 < 0.d0) print*, 'Error: int1 + int2 < 0' + elseif (abs(d1) < abs(d2)) then + pint3D = pint3D*(int2 - int1) + if (int2 - int1 < 0.d0) print*, 'Error: int2 - int1 < 0: ', int1, int2, '(', d1, d2,')' + else + pint3D = pint3D*(int1 - int2) + if (int1 - int2 < 0.d0) print*, 'Error: int1 - int2 < 0: ', int1, int2, '(', d1, d2,')' + endif + +end function pint3D + +real(doub_prec) function full_integral_3D(d, r0, R_0, h) + + real(doub_prec), intent(in) :: d, r0, R_0, h + real(doub_prec) :: B1, B2, B3, a, h2 + real(doub_prec), parameter :: pi = 4.*atan(1.) + real(doub_prec) :: tanphi, phi, a2, cosp, r0h, r03, r0h2, r0h3, r0h_2, r0h_3 + real(doub_prec) :: r2, R_, linedist2, cosphi + real(doub_prec) :: I0, I1, I_2, I_3, I_4, I_5 + real(doub_prec) :: D2, D3 + + r0h = r0/h + tanphi = abs(d)/R_0 + phi = atan(tanphi) + + if (abs(r0h) < tiny(0.) .or. abs(R_0/h) < tiny(0.) .or. abs(phi) < tiny(0.)) then + full_integral_3D = 0.0 + return + endif + + h2 = h*h + r03 = r0*r0*r0 + r0h2 = r0h*r0h + r0h3 = r0h2*r0h + r0h_2 = 1./r0h2 + r0h_3 = 1./r0h3 + + ! Avoid Compiler warnings + B1 = 0. + B2 = 0. + + if (r0 >= 2.0*h) then + B3 = 0.25*h2*h + elseif (r0 > h) then + B3 = 0.25*r03 *(-4./3. + (r0h) - 0.3*r0h2 + 1./30.*r0h3 - 1./15. *r0h_3+ 8./5.*r0h_2) + B2 = 0.25*r03 *(-4./3. + (r0h) - 0.3*r0h2 + 1./30.*r0h3 - 1./15. *r0h_3) + else + B3 = 0.25*r03 *(-2./3. + 0.3*r0h2 - 0.1*r0h3 + 7./5.*r0h_2) + B2 = 0.25*r03 *(-2./3. + 0.3*r0h2 - 0.1*r0h3 - 1./5.*r0h_2) + B1 = 0.25*r03 *(-2./3. + 0.3*r0h2 - 0.1*r0h3) + endif + + a = R_0/r0 + a2 = a*a + + linedist2 = (r0*r0 + R_0*R_0) + cosphi = cos(phi) + R_ = R_0/cosphi + r2 = (r0*r0 + R_*R_) + + D2 = 0.0 + D3 = 0.0 + + if (linedist2 < h2) then + !////// phi1 business ///// + cosp = R_0/sqrt(h2-r0*r0) + call get_I_terms(cosp,a2,a,I0,I1,I_2,I_3,I_4,I_5) + + D2 = -1./6.*I_2 + 0.25*(r0h) *I_3 - 0.15*r0h2 *I_4 + 1./30.*r0h3 *I_5 - 1./60. *r0h_3 *I1 + (B1-B2)/r03 *I0 + endif + if (linedist2 < 4.*h2) then + !////// phi2 business ///// + cosp = R_0/sqrt(4.0*h2-r0*r0) + call get_I_terms(cosp,a2,a,I0,I1,I_2,I_3,I_4,I_5) + + D3 = 1./3.*I_2 - 0.25*(r0h) *I_3 + 3./40.*r0h2 *I_4 - 1./120.*r0h3 *I_5 + 4./15. *r0h_3 *I1 + (B2-B3)/r03 *I0 + D2 + endif + + !////////////////////////////// + call get_I_terms(cosphi,a2,a,I0,I1,I_2,I_3,I_4,I_5,phi=phi,tanphi=tanphi) + + if (r2 < h2) then + full_integral_3D = r0h3/pi * (1./6. *I_2 - 3./40.*r0h2 *I_4 + 1./40.*r0h3 *I_5 + B1/r03 *I0) + elseif (r2 < 4.*h2) then + full_integral_3D= r0h3/pi * (0.25 * (4./3. *I_2 - (r0/h) *I_3 + 0.3*r0h2 *I_4 - & + & 1./30.*r0h3 *I_5 + 1./15. *r0h_3 *I1) + B2/r03 *I0 + D2) + else + full_integral_3D = r0h3/pi * (-0.25*r0h_3 *I1 + B3/r03 *I0 + D3) + endif + +end function full_integral_3D + +subroutine get_I_terms(cosp,a2,a,I0,I1,I_2,I_3,I_4,I_5,phi,tanphi) + real(doub_prec), intent(in) :: cosp,a2,a + real(doub_prec), intent(out) :: I0,I1,I_2,I_3,I_4,I_5 + real(doub_prec), intent(in), optional :: phi,tanphi + real(doub_prec) :: cosp2,p,tanp,u2,u,logs,I_1,mu2_1,fac + + cosp2 = cosp*cosp + if (present(phi)) then + p = phi + tanp = tanphi + else + p = acos(cosp) + tanp = sqrt(1.-cosp2)/cosp ! tan(p) + endif + + mu2_1 = 1. / (1. + cosp2/a2) + I0 = p + I_2 = p + a2 * tanp + I_4 = p + 2.*a2 * tanp + 1./3.*a2*a2 * tanp*(2. + 1./cosp2) + + u2 = (1.-cosp2)*mu2_1 + u = sqrt(u2) + logs = log((1.+u)/(1.-u)) + I1 = atan2(u,a) + + fac = 1./(1.-u2) + I_1 = 0.5*a*logs + I1 + I_3 = I_1 + a*0.25*(1.+a2)*(2.*u*fac + logs) + I_5 = I_3 + a*(1.+a2)*(1.+a2)/16. *( (10.*u - 6.*u*u2)*fac*fac + 3.*logs) + +end subroutine get_I_terms + + !------------------------------------------------------------ + ! function to return a soft maximum for 1/x with no bias + ! for x >> eps using the cubic spline kernel softening + ! i.e. something equivalent to 1/sqrt(x**2 + eps**2) but + ! with compact support, i.e. f=1/x when x > 2*eps + !------------------------------------------------------------ +pure elemental real function soft_func(x,eps) result(f) + real, intent(in) :: x,eps + real :: q,q2,q4 + + q = x/eps + q2 = q*q + if (q < 1.) then + q4 = q2*q2 + f = (1./eps)*(q4*q/10. - 3.*q4/10. + 2.*q2/3. - 7./5.) + elseif (q < 2.) then + q4 = q2*q2 + f = (1./eps)*(q*(-q4*q + 9.*q4 - 30.*q2*q + 40.*q2 - 48.) + 2.)/(30.*q) + else + f = -1./x + endif + f = -f + +end function soft_func + + !-------------------------------------------------------------------------- + ! + ! utility to wrap pixel index around periodic domain + ! indices that roll beyond the last position are re-introduced at the first + ! + !-------------------------------------------------------------------------- +pure integer function iroll(i,n) + integer, intent(in) :: i,n + + if (i > n) then + iroll = mod(i-1,n) + 1 + elseif (i < 1) then + iroll = n + mod(i,n) ! mod is negative + else + iroll = i + endif + +end function iroll +end module interpolations3D + diff --git a/src/utils/interpolate3Dold.F90 b/src/utils/interpolate3Dold.F90 new file mode 100644 index 000000000..d1344fd96 --- /dev/null +++ b/src/utils/interpolate3Dold.F90 @@ -0,0 +1,367 @@ +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2023 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.bitbucket.io/ ! +!--------------------------------------------------------------------------! +module interpolations3D +! +! Module containing routine for interpolation from PHANTOM data +! to 3D adaptive mesh +! +! Requires adaptivemesh.f90 module +! +! :References: None +! +! :Owner: Spencer Magnall +! +! :Runtime parameters: None +! +! :Dependencies: kernel +! + + implicit none + real, parameter, private :: dpi = 1./3.1415926536d0 + public :: interpolate3D +!$ integer(kind=8), dimension(:), private, allocatable :: ilock + +contains +!-------------------------------------------------------------------------- +! subroutine to interpolate from particle data to even grid of pixels +! +! The data is interpolated according to the formula +! +! datsmooth(pixel) = sum_b weight_b dat_b W(r-r_b, h_b) +! +! where _b is the quantity at the neighbouring particle b and +! W is the smoothing kernel, for which we use the usual cubic spline. +! +! For a standard SPH smoothing the weight function for each particle should be +! +! weight = pmass/(rho*h^3) +! +! this version is written for slices through a rectangular volume, ie. +! assumes a uniform pixel size in x,y, whilst the number of pixels +! in the z direction can be set to the number of cross-section slices. +! +! Input: particle coordinates and h : xyzh(4,npart) +! weight for each particle : weight [ same on all parts in PHANTOM ] +! scalar data to smooth : dat (npart) +! +! Output: smoothed data : datsmooth (npixx,npixy,npixz) +! +! Daniel Price, Monash University 2010 +! daniel.price@monash.edu +!-------------------------------------------------------------------------- + +subroutine interpolate3D(xyzh,weight,npart, & + xmin,datsmooth,nnodes,dxgrid,normalise,dat,ngrid,vertexcen) + use kernel, only:wkern, radkern, radkern2, cnormk + !use adaptivemesh, only:ifirstlevel,nsub,ndim,gridnodes + integer, intent(in) :: npart,nnodes,ngrid(3) + real, intent(in) :: xyzh(:,:)! ,vxyzu(:,:) + real, intent(in) :: weight !,pmass + real, intent(in) :: xmin(3),dxgrid(3) + real, intent(out) :: datsmooth(:,:,:) + logical, intent(in) :: normalise, vertexcen + real, intent(in), optional :: dat(:) + real, allocatable :: datnorm(:,:,:) +! real, dimension(nsub**ndim,nnodes) :: datnorm + integer, parameter :: ndim = 3, nsub=1 + integer :: i,ipix,jpix,kpix,isubmesh,imesh,level,icell + integer :: iprintinterval,iprintnext + integer :: ipixmin,ipixmax,jpixmin,jpixmax,kpixmin,kpixmax + integer :: ipixi,jpixi,kpixi,npixx,npixy,npixz + real :: xi,yi,zi,hi,hi1,hi21,radkernh,qq,wab,q2,const,dyz2,dz2 + real :: xorigi,yorigi,zorigi,xpix,ypix,zpix,dx,dy,dz + real :: dxcell(ndim),xminnew(ndim), dxmax(ndim) + real :: t_start,t_end + real :: termnorm + real :: term + logical :: iprintprogress +!$ integer :: omp_get_num_threads,j +#ifndef _OPENMP + integer(kind=8) :: iprogress +#endif + + print*, "size: ", size(datsmooth) + print*, "datsmooth out of bounds: ", datsmooth(35,1,1) + datsmooth = 0. + dxmax(:) = dxgrid(:) + !datnorm = 0. + if (normalise) then + print "(1x,a)",'interpolating from particles to Einstein toolkit grid (normalised) ...' + else + print "(1x,a)",'interpolating from particles to Einstein toolkit grid (non-normalised) ...' + endif +! if (any(dxmax(:) <= 0.)) then +! print "(1x,a)",'interpolate3D: error: grid size <= 0' +! return +! endif +! if (ilendat /= 0) then +! print "(1x,a)",'interpolate3D: error in interface: dat has non-zero length but is not present' +! return +! endif + if (normalise) then + allocate(datnorm(ngrid(1),ngrid(2),ngrid(3))) + datnorm = 0. + endif + +!$ allocate(ilock(0:nnodes)) +!$ do i=0,nnodes +!$ call omp_init_lock(ilock(i)) +!$ enddo + + ! + !--print a progress report if it is going to take a long time + ! (a "long time" is, however, somewhat system dependent) + ! + iprintprogress = (npart >= 100000) .or. (nnodes > 10000) + ! + !--loop over particles + ! + iprintinterval = 25 + if (npart >= 1e6) iprintinterval = 10 + iprintnext = iprintinterval + ! + !--get starting CPU time + ! + call cpu_time(t_start) + + imesh = 1 + level = 1 + dxcell(:) = dxgrid(:)/real(nsub**level) +! xminpix(:) = xmin(:) - 0.5*dxcell(:) + npixx = ngrid(1) + npixy = ngrid(2) + npixz = ngrid(3) + print "(3(a,i4))",' root grid: ',npixx,' x ',npixy,' x ',npixz + print*, "position of i cell is: ", 1*dxcell(1) + xmin(1) + print*, "npart: ", npart + + const = cnormk ! kernel normalisation constant (3D) + print*,"const: ", const + !stop + + ! + !--loop over particles + ! + !$omp parallel default(none) & + !$omp shared(npart,xyzh,dat,datsmooth,datnorm,vertexcen,const,weight) & + !$omp shared(xmin,imesh,nnodes,level) & + !$omp shared(npixx,npixy,npixz,dxmax,dxcell,normalise) & + !$omp private(i,j,hi,hi1,hi21,termnorm,term) & + !$omp private(xpix,ypix,zpix,dx,dy,dz,dz2,dyz2,qq,q2,wab,radkernh) & + !$omp private(xi,yi,zi,xorigi,yorigi,zorigi,xminnew) & + !$omp private(ipix,jpix,kpix,ipixi,jpixi,kpixi,icell,isubmesh) & + !$omp private(ipixmin,ipixmax,jpixmin,jpixmax,kpixmin,kpixmax) + !$omp master +!$ print "(1x,a,i3,a)",'Using ',omp_get_num_threads(),' cpus' + !$omp end master + !$omp do schedule(guided,10) + over_parts: do i=1,npart + ! + !--report on progress + ! + !print*, i +#ifndef _OPENMP + if (iprintprogress) then + iprogress = nint(100.*i/npart) + if (iprogress >= iprintnext) then + write(*,"('(',i3,'% -',i12,' particles done)')") iprogress,i + iprintnext = iprintnext + iprintinterval + endif + endif +#endif + ! + !--set kernel related quantities + ! + xi = xyzh(1,i); xorigi = xi + yi = xyzh(2,i); yorigi = yi + zi = xyzh(3,i); zorigi = zi + hi = xyzh(4,i) + radkernh = radkern*hi + !print*, "hi: ", hi + if (hi <= 0.) cycle over_parts + hi1 = 1./hi; hi21 = hi1*hi1 + termnorm = const*weight + ! print*, "const: ", const + ! print*, "weight: ", weight + ! print*, "termnorm: ", termnorm + + !radkern = 2.*hi ! radius of the smoothing kernel + !print*, "radkern: ", radkern + !print*, "part pos: ", xi,yi,zi + term = termnorm*dat(i) ! weight for density calculation + ! I don't understand why this doesnt involve any actual smoothing? + !dfac = hi**3/(dxcell(1)*dxcell(2)*dxcell(3)*const) + ! + !--for each particle work out which pixels it contributes to + ! + !print*, "radkern: ", radkern + ipixmin = int((xi - radkernh - xmin(1))/dxcell(1)) + jpixmin = int((yi - radkernh - xmin(2))/dxcell(2)) + kpixmin = int((zi - radkernh - xmin(3))/dxcell(3)) + + ipixmax = int((xi + radkernh - xmin(1))/dxcell(1)) + 1 + jpixmax = int((yi + radkernh - xmin(2))/dxcell(2)) + 1 + kpixmax = nint((zi + radkernh - xmin(3))/dxcell(3)) + 1 + + !if (ipixmax == 33) stop + + + !if (ipixmin == 4 .and. jpixmin == 30 .and. kpixmin == 33) print*, "particle (min): ", i + !if (ipixmax == 4 .and. jpixmax == 30 .and. kpixmax == 33) print*, "particle (max): ", i +#ifndef PERIODIC + if (ipixmin < 1) ipixmin = 1 ! make sure they only contribute + if (jpixmin < 1) jpixmin = 1 ! to pixels in the image + if (kpixmin < 1) kpixmin = 1 + if (ipixmax > npixx) ipixmax = npixx + if (jpixmax > npixy) jpixmax = npixy + if (kpixmax > npixz) kpixmax = npixz + !print*, "ipixmin: ", ipixmin + !print*, "ipixmax: ", ipixmax + !print*, "jpixmin: ", jpixmin + !print*, "jpixmax: ", jpixmax + !print*, "kpixmin: ", kpixmin + !print*, "kpixmax: ", kpixmax +#endif + !print*,' part ',i,' lims = ',ipixmin,ipixmax,jpixmin,jpixmax,kpixmin,kpixmax + ! + !--loop over pixels, adding the contribution from this particle + ! (note that we handle the periodic boundary conditions + ! entirely on the root grid) + ! + do kpix = kpixmin,kpixmax + kpixi = kpix +#ifdef PERIODIC + if (kpixi < 1) then + kpixi = kpixi + npixz + zi = zorigi !+ dxmax(3) + elseif (kpixi > npixz) then + kpixi = kpixi - npixz + zi = zorigi !- dxmax(3) + else + zi = zorigi + endif +#endif + if (vertexcen) then + zpix = xmin(3) + (kpixi-1)*dxcell(3) + else + zpix = xmin(3) + (kpixi-0.5)*dxcell(3) + endif + dz = zpix - zi + dz2 = dz*dz*hi21 + + do jpix = jpixmin,jpixmax + jpixi = jpix +#ifdef PERIODIC + if (jpixi < 1) then + jpixi = jpixi + npixy + yi = yorigi !+ dxmax(2) + elseif (jpixi > npixy) then + jpixi = jpixi - npixy + yi = yorigi !- dxmax(2) + else + yi = yorigi + endif +#endif + if (vertexcen) then + ypix = xmin(2) + (jpixi-1)*dxcell(2) + else + ypix = xmin(2) + (jpixi-0.5)*dxcell(2) + endif + dy = ypix - yi + dyz2 = dy*dy*hi21 + dz2 + + do ipix = ipixmin,ipixmax + ipixi = ipix +#ifdef PERIODIC + if (ipixi < 1) then + ipixi = ipixi + npixx + xi = xorigi !+ dxmax(1) + elseif (ipixi > npixx) then + if (ipixi == 33) then + print*,"xi old: ", xorigi + print*, "xi new: ", xorigi-dxmax(1) + print*, "ipixi new: ", ipixi - npixx + endif + ipixi = ipixi - npixx + xi = xorigi !- dxmax(1) + else + xi = xorigi + endif +#endif + icell = ((kpixi-1)*nsub + (jpixi-1))*nsub + ipixi + ! + !--particle interpolates directly onto the root grid + ! + !print*,'onto root grid ',ipixi,jpixi,kpixi + if (vertexcen) then + xpix = xmin(1) + (ipixi-1)*dxcell(1) + else + xpix = xmin(1) + (ipixi-0.5)*dxcell(1) + endif + !print*, "xpix: ", xpix + !xpix = xmin(1) + (ipixi-1)*dxcell(1) ! Since we are vertex centered from Et + dx = xpix - xi + q2 = dx*dx*hi21 + dyz2 ! dx2 pre-calculated; dy2 pre-multiplied by hi21 + ! + !--SPH kernel - standard cubic spline + ! + if (q2 < radkern2) then + ! if (q2 < 1.0) then + ! qq = sqrt(q2) + ! wab = 1.-1.5*q2 + 0.75*q2*qq + ! else + ! qq = sqrt(q2) + ! wab = 0.25*(2.-qq)**3 + ! endif + ! Call the kernel routine + qq = sqrt(q2) + wab = wkern(q2,qq) + ! + !--calculate data value at this pixel using the summation interpolant + ! + ! Change this to the access the pixel coords x,y,z + !$omp critical + datsmooth(ipixi,jpixi,kpixi) = datsmooth(ipixi,jpixi,kpixi) + term*wab + + !if (ipixi==1 .and. jpixi==1 .and. kpixi==1) print*, "x position of 1,1,1", xi,yi,zi + if (normalise) then + datnorm(ipixi,jpixi,kpixi) = datnorm(ipixi,jpixi,kpixi) + termnorm*wab + endif + !$omp end critical + endif + enddo + enddo + enddo + enddo over_parts + !$omp enddo + !$omp end parallel + +!$ do i=0,nnodes +!$ call omp_destroy_lock(ilock(i)) +!$ enddo +!$ if (allocated(ilock)) deallocate(ilock) + + ! + !--normalise dat array + ! + if (normalise) then + where (datnorm > tiny(datnorm)) + datsmooth = datsmooth/datnorm + end where + endif + if (allocated(datnorm)) deallocate(datnorm) + ! + !--get ending CPU time + ! + call cpu_time(t_end) + print*,'completed in ',t_end-t_start,'s' + + return + +end subroutine interpolate3D + +end module interpolations3D