diff --git a/.mailmap b/.mailmap index 8332207b8..eab7608f0 100644 --- a/.mailmap +++ b/.mailmap @@ -96,7 +96,8 @@ Megha Sharma megha sharma Megha Sharma Megha Sharma Megha Sharma Megha Sharma Megha Sharma Megha Sharma -Mike Lau Mike Lau <55525335+themikelau@users.noreply.github.com> +Mike Lau Mike Lau <55525335+themikelau@users.noreply.github.com> +Mike Lau Mike Lau Elisabeth Borchert emborchert <69176538+emborchert@users.noreply.github.com> Ward Homan ward Ward Homan wardhoman <33419533+wardhoman@users.noreply.github.com> @@ -104,6 +105,7 @@ Benedetta Veronesi benedetta veronesi Phantom benchmark bot Ubuntu Stephane Michoulier StephaneMichoulier +Stephane Michoulier MICHOULIER Stephane Jolien Malfait Jolien128 <72729152+Jolien128@users.noreply.github.com> Martina Toscani Martina Toscani @@ -115,4 +117,11 @@ Amena Faruqi <42060670+amenafaruqi@users.noreply.gi Amena Faruqi Amena Faruqi Alison Young Alison Young Simone Ceppi Simone Ceppi +Madeline Overton Madeline Nicole Overton +Madeline Overton Madeline Overton <85810161+moverton000@users.noreply.github.com> Nicolás Cuello Nicolas Cuello +Rebecca Martin rebeccagmartin <74937128+rebeccagmartin@users.noreply.github.com> +Stephen Nielson s-neilson <36410751+s-neilson@users.noreply.github.com> +Stephen Nielson Stephen Neilson <36410751+s-neilson@users.noreply.github.com> +Yann Bernard Yrisch +David Bamba DavidBamba diff --git a/AUTHORS b/AUTHORS index 649e5eb20..b3a77dfde 100644 --- a/AUTHORS +++ b/AUTHORS @@ -6,14 +6,14 @@ # Edit .mailmap if your name or email are wrong # #-------------------------------------------------------# Daniel Price -Mike Lau +Mike Lau Conrad Chan James Wurster David Liptai Lionel Siess Fangyi (Fitz) Hu Daniel Mentiplay -Yrisch +Yann Bernard Megha Sharma Arnaud Vericel Mark Hutchison @@ -26,8 +26,8 @@ Terrence Tricco Stephane Michoulier Simone Ceppi Spencer Magnall -Caitlyn Hardiman Enrico Ragusa +Caitlyn Hardiman Cristiano Longarini Sergei Biriukov Giovanni Dipierro @@ -35,41 +35,41 @@ Roberto Iaconi Amena Faruqi Hauke Worpel Alison Young -Stephen Neilson <36410751+s-neilson@users.noreply.github.com> +Stephen Nielson Martina Toscani Benedetta Veronesi -Sahl Rowther Simon Glover +Sahl Rowther Thomas Reichardt Jean-François Gonzalez Christopher Russell -Alessia Franchini +Madeline Overton Alex Pettitt -Jolien Malfait Phantom benchmark bot -Kieran Hirsh -Mike Lau +Jolien Malfait +Alessia Franchini Nicole Rodrigues -David Trevascus +Kieran Hirsh Farzana Meru +David Trevascus Nicolás Cuello Chris Nixon Miguel Gonzalez-Bolivar -Benoit Commercon -Giulia Ballabio -Joe Fisher Maxime Lombart -Orsola De Marco Zachary Pellow -s-neilson <36410751+s-neilson@users.noreply.github.com> -Ariel Chitan -Chunliang Mu -Cox, Samuel +Orsola De Marco +Joe Fisher +Benoit Commercon +Giulia Ballabio +Rebecca Martin +Jorge Cuadra Hugh Griffiths Jeremy Smallwood -Jorge Cuadra -MICHOULIER Stephane +David Bamba +Cox, Samuel +Chunliang Mu +Shunquan Huang Steven Rieder Stéven Toupin Taj Jankovič -rebeccagmartin <74937128+rebeccagmartin@users.noreply.github.com> +Ariel Chitan diff --git a/build/Makefile b/build/Makefile index 5acba64f2..201fa7a45 100644 --- a/build/Makefile +++ b/build/Makefile @@ -488,7 +488,7 @@ ifdef METRIC else SRCMETRIC= metric_minkowski.f90 endif -SRCGR=inverse4x4.f90 einsteintk_utils.f90 $(SRCMETRIC) metric_tools.f90 utils_gr.f90 interpolate3D.f90 tmunu2grid.f90 +SRCGR=inverse4x4.f90 metric_et_utils.f90 einsteintk_utils.f90 $(SRCMETRIC) metric_tools.f90 utils_gr.f90 interpolate3D.f90 tmunu2grid.f90 # # chemistry and cooling # @@ -1230,6 +1230,24 @@ combinedustdumps: checksys checkparams $(OBJCDD) cleancombinedustdumps: rm -f $(BINDIR)/combinedustdumps +#---------------------------------------------------- +# these are the sources for the tabulate_metric tility +.PHONY: tabulate_metric +SRCTAB= io.F90 utils_infiles.f90 metric_${METRIC}.f90 metric_et_utils.f90 tabulate_metric.f90 #metric_tools.F90 +OBJTAB1= $(SRCTAB:.F90=.o) +OBJTAB= $(OBJTAB1:.f90=.o) + +tabulate_metric: checksys $(OBJTAB) + @echo "" + @echo "tabulate_metric: Because grids are great" + @echo "" + $(FC) $(FFLAGS) -o $(BINDIR)/$@ $(OBJTAB) + +cleantabulatemetric: + rm -f $(BINDIR)/tabulate_metric + + + include Makefile_qscripts diff --git a/build/Makefile_setups b/build/Makefile_setups index 94c9d7e0d..f6dd95b81 100644 --- a/build/Makefile_setups +++ b/build/Makefile_setups @@ -77,7 +77,17 @@ endif ifeq ($(SETUP), asteroidwind) # asteroid emitting a wind (Trevascus et al. 2021) SETUPFILE=setup_asteroidwind.f90 - SRCINJECT=utils_binary.f90 inject_asteroidwind.f90 + SRCINJECT=utils_binary.f90 inject_randomwind.f90 + IND_TIMESTEPS=yes + CONST_AV=yes + ISOTHERMAL=yes + KNOWN_SETUP=yes +endif + +ifeq ($(SETUP), randomwind) +# asteroid emitting a wind (Trevascus et al. 2021) + SETUPFILE=setup_disc.f90 + SRCINJECT=utils_binary.f90 inject_randomwind.f90 IND_TIMESTEPS=yes CONST_AV=yes ISOTHERMAL=yes diff --git a/docs/developer-guide/datafiles.rst b/docs/developer-guide/datafiles.rst index 5ceaa8f03..0dc469ac6 100644 --- a/docs/developer-guide/datafiles.rst +++ b/docs/developer-guide/datafiles.rst @@ -8,7 +8,7 @@ order for your modules to be portable. Small files ----------- -For *very small files* (under 1Mb), you can simply add these to the git +For *very small files* (under 100Kb), you can simply add these to the git repository in a subdirectory of the phantom/data directory: :: diff --git a/src/main/H2regions.f90 b/src/main/H2regions.f90 index 8d81b4a12..cac51a26e 100644 --- a/src/main/H2regions.f90 +++ b/src/main/H2regions.f90 @@ -6,21 +6,17 @@ !--------------------------------------------------------------------------! module HIIRegion ! -! HIIRegion +! Feedback from HII regions ! ! :References: Fujii et al. (2021), Hopkins et al. (2012) ! -! :Owner: Yrisch +! :Owner: Yann Bernard ! ! :Runtime parameters: None ! ! :Dependencies: dim, eos, infile_utils, io, linklist, part, physcon, ! sortutils, timing, units ! -! contains routines to model HII region expansion due to ionization and radiation pressure.. -! routine originally made by Hopkins et al. (2012),reused by Fujii et al. (2021) -! and adapted in Phantom by Yann Bernard - implicit none public :: update_ionrates,update_ionrate, HII_feedback,initialize_H2R,read_options_H2R,write_options_H2R @@ -52,17 +48,18 @@ module HIIRegion contains - !----------------------------------------------------------------------- - !+ - ! Initialise stellar feedbacks - !+ - !----------------------------------------------------------------------- +!----------------------------------------------------------------------- +!+ +! Initialise stellar feedbacks +!+ +!----------------------------------------------------------------------- subroutine initialize_H2R use io, only:iprint,iverbose,id,master use part, only:isionised use units, only:udist,umass,utime use physcon, only:mass_proton_cgs,kboltz,pc,eV,solarm - use eos , only:gmw,gamma + use eos, only:gmw,gamma + isionised(:)=.false. !calculate the useful constant in code units mH = gmw*mass_proton_cgs @@ -85,7 +82,7 @@ subroutine initialize_H2R write(iprint,"(a,es18.10,es18.10)") " Max strögrem radius (code/pc) : ", Rst_max, Rmax write(iprint,"(a,es18.10,es18.10)") " Min feedback mass (code/Msun) : ", Minmass, Mmin endif - return + end subroutine initialize_H2R !----------------------------------------------------------------------- @@ -93,7 +90,6 @@ end subroutine initialize_H2R ! Calculation of the the ionizing photon rate of all stars (Only for restart) !+ !----------------------------------------------------------------------- - subroutine update_ionrates(nptmass,xyzmh_ptmass,h_acc) use io, only:iprint,iverbose use units, only:umass @@ -104,6 +100,7 @@ subroutine update_ionrates(nptmass,xyzmh_ptmass,h_acc) real, intent(in) :: h_acc real :: logmi,log_Q,mi,hi integer :: i + nHIIsources = 0 !$omp parallel do default(none) & !$omp shared(xyzmh_ptmass,iprint,iverbose,umass)& @@ -122,7 +119,7 @@ subroutine update_ionrates(nptmass,xyzmh_ptmass,h_acc) xyzmh_ptmass(irstrom,i) = -1. nHIIsources = nHIIsources + 1 if (iverbose >= 0) then - write(iprint,"(/a,es18.10,es18.10/)")"Massive stars detected : Log Q, Mass : ",log_Q,mi + write(iprint,"(/a,es18.10,es18.10/)") "Massive stars detected : Log Q, Mass : ",log_Q,mi endif else xyzmh_ptmass(irateion,i) = -1. @@ -133,9 +130,14 @@ subroutine update_ionrates(nptmass,xyzmh_ptmass,h_acc) if (iverbose > 1) then write(iprint,"(/a,i8/)") "nb_feedback sources : ",nHIIsources endif - return + end subroutine update_ionrates +!----------------------------------------------------------------------- +!+ +! update the ionizing photon rate +!+ +!----------------------------------------------------------------------- subroutine update_ionrate(i,xyzmh_ptmass,h_acc) use io, only:iprint,iverbose use units, only:umass @@ -145,6 +147,7 @@ subroutine update_ionrate(i,xyzmh_ptmass,h_acc) real, intent(inout) :: xyzmh_ptmass(:,:) real, intent(in) :: h_acc real :: logmi,log_Q,mi,hi + mi = xyzmh_ptmass(4,i) hi = xyzmh_ptmass(ihacc,i) if (mi > Minmass .and. hi < h_acc) then @@ -156,7 +159,7 @@ subroutine update_ionrate(i,xyzmh_ptmass,h_acc) xyzmh_ptmass(irstrom,i) = -1. nHIIsources = nHIIsources + 1 if (iverbose >= 0) then - write(iprint,"(/a,es18.10,es18.10/)")"Massive stars detected : Log Q, Mass : ",log_Q,mi + write(iprint,"(/a,es18.10,es18.10/)") "Massive stars detected : Log Q, Mass : ",log_Q,mi endif else xyzmh_ptmass(irateion,i) = -1. @@ -166,15 +169,14 @@ subroutine update_ionrate(i,xyzmh_ptmass,h_acc) if (iverbose > 1) then write(iprint,"(/a,i8/)") "nb_feedback sources : ",nHIIsources endif - return -end subroutine update_ionrate - !----------------------------------------------------------------------- - !+ - ! Main subroutine : Application of the HII feedback using Hopkins's like prescription - !+ - !----------------------------------------------------------------------- +end subroutine update_ionrate +!----------------------------------------------------------------------- +!+ +! Main subroutine : Application of the HII feedback using Hopkins's like prescription +!+ +!----------------------------------------------------------------------- subroutine HII_feedback(nptmass,npart,xyzh,xyzmh_ptmass,vxyzu,isionised,dt) use part, only:rhoh,massoftype,ihsoft,igas,irateion,isdead_or_accreted,& irstrom @@ -268,10 +270,10 @@ subroutine HII_feedback(nptmass,npart,xyzh,xyzmh_ptmass,vxyzu,isionised,dt) r_in = sqrt((xi-xyzh(1,j))**2 + (yi-xyzh(2,j))**2 + (zi-xyzh(3,j))**2) mHII = ((4.*pi*(r**3-r_in**3)*rhoh(xyzh(4,j),pmass))/3) if (mHII>3*pmass) then -!$omp parallel do default(none) & -!$omp shared(mHII,listneigh,xyzh,sigd,dt) & -!$omp shared(mH,vxyzu,log_Qi,hv_on_c,npartin,pmass,xi,yi,zi) & -!$omp private(j,dx,dy,dz,vkx,vky,vkz,xj,yj,zj,r,taud) + !$omp parallel do default(none) & + !$omp shared(mHII,listneigh,xyzh,sigd,dt) & + !$omp shared(mH,vxyzu,log_Qi,hv_on_c,npartin,pmass,xi,yi,zi) & + !$omp private(j,dx,dy,dz,vkx,vky,vkz,xj,yj,zj,r,taud) do k=1,npartin j = listneigh(1) xj = xyzh(1,j) @@ -290,28 +292,40 @@ subroutine HII_feedback(nptmass,npart,xyzh,xyzmh_ptmass,vxyzu,isionised,dt) vxyzu(2,j) = vxyzu(2,j) + vky*dt vxyzu(3,j) = vxyzu(3,j) + vkz*dt enddo -!$omp end parallel do + !$omp end parallel do endif endif enddo endif call get_timings(t2,tcpu2) call increment_timer(itimer_HII,t2-t1,tcpu2-tcpu1) - return + end subroutine HII_feedback +!----------------------------------------------------------------------- +!+ +! write options to input file +!+ +!----------------------------------------------------------------------- subroutine write_options_H2R(iunit) use infile_utils, only:write_inopt use physcon, only:solarm integer, intent(in) :: iunit - write(iunit,"(/,a)") '# options controlling HII region expansion feedback' - if (iH2R>0) then + + if (iH2R > 0) then + write(iunit,"(/,a)") '# options controlling HII region expansion feedback' call write_inopt(iH2R, 'iH2R', "enable the HII region expansion feedback in star forming reigon", iunit) call write_inopt(Mmin, 'Mmin', "Minimum star mass to trigger HII region (MSun)", iunit) call write_inopt(Rmax, 'Rmax', "Maximum radius for HII region (pc)", iunit) endif + end subroutine write_options_H2R +!----------------------------------------------------------------------- +!+ +! read options from input file +!+ +!----------------------------------------------------------------------- subroutine read_options_H2R(name,valstring,imatch,igotall,ierr) use io, only:fatal character(len=*), intent(in) :: name,valstring @@ -319,6 +333,7 @@ subroutine read_options_H2R(name,valstring,imatch,igotall,ierr) integer, intent(out) :: ierr integer, save :: ngot = 0 character(len=30), parameter :: label = 'read_options_H2R' + imatch = .true. select case(trim(name)) case('iH2R') @@ -337,6 +352,7 @@ subroutine read_options_H2R(name,valstring,imatch,igotall,ierr) imatch = .true. end select igotall = (ngot >= 3) + end subroutine read_options_H2R end module HIIRegion diff --git a/src/main/cons2prim.f90 b/src/main/cons2prim.f90 index a8bb68ffc..56878217a 100644 --- a/src/main/cons2prim.f90 +++ b/src/main/cons2prim.f90 @@ -6,9 +6,19 @@ !--------------------------------------------------------------------------! module cons2prim ! -! None +! Subroutines to swap between primitive variables needed on RHS of fluid +! equations (density,velocity,internal energy) and conserved/evolved +! variables on the LHS of the fluid equations (rho*,momentum,entropy) ! -! :References: None +! This is complicated in the GR code but also useful to structure +! things this way in the non-GR code, e.g. B/rho is the evolved variable +! while B is the "primitive" variable for magnetic field. Similarly +! sqrt(rho_d) is the evolved variable for dust species but the primitive +! variable is the dust mass fraction. +! +! :References: +! Liptai & Price (2019), MNRAS 485, 819-842 +! Ballabio et al. (2018), MNRAS 477, 2766-2771 ! ! :Owner: Elisabeth Borchert ! @@ -26,12 +36,14 @@ module cons2prim contains -!------------------------------------- -! -! Primitive to conservative routines -! -!------------------------------------- - +!---------------------------------------------------------------------- +!+ +! Primitive to conservative transform (for GR): +! Construct conserved variables (rho*,momentum,entropy) +! from the primitive/fluid rest frame variables +! (density,velocity,internal energy), for ALL particles +!+ +!---------------------------------------------------------------------- subroutine prim2consall(npart,xyzh,metrics,vxyzu,dens,pxyzu,use_dens) use part, only:isdead_or_accreted,ien_type,eos_vars,igasP,igamma,itemp use eos, only:gamma,ieos @@ -75,6 +87,12 @@ subroutine prim2consall(npart,xyzh,metrics,vxyzu,dens,pxyzu,use_dens) end subroutine prim2consall +!---------------------------------------------------------------------- +!+ +! Primitive to conservative transform (for GR): +! for a single SPH particle +!+ +!---------------------------------------------------------------------- subroutine prim2consi(xyzhi,metrici,vxyzui,dens_i,pri,tempi,pxyzui,use_dens,ien_type) use cons2primsolver, only:primitive2conservative use utils_gr, only:h2dens @@ -112,12 +130,13 @@ subroutine prim2consi(xyzhi,metrici,vxyzui,dens_i,pri,tempi,pxyzui,use_dens,ien_ end subroutine prim2consi -!--------------------------------------------- -! -! Conservative to primitive routines (for GR) -! -!--------------------------------------------- - +!---------------------------------------------------------------------- +!+ +! Conservative to primitive routines (for GR): +! Solve for primitive variables (density,velocity,internal energy) +! from the evolved/conservative variables (rho*,momentum,entropy) +!+ +!---------------------------------------------------------------------- subroutine cons2primall(npart,xyzh,metrics,pxyzu,vxyzu,dens,eos_vars) use cons2primsolver, only:conservative2primitive use part, only:isdead_or_accreted,massoftype,igas,rhoh,igasP,ics,ien_type,& @@ -165,12 +184,15 @@ subroutine cons2primall(npart,xyzh,metrics,pxyzu,vxyzu,dens,eos_vars) end subroutine cons2primall -!------------------------------------- -! -! Primitive variables from conservative variables +!----------------------------------------------------------------------------- +!+ +! Solve for primitive variables (v,u,P,B,dustfrac) from evolved variables +! (v,energy variable,B/rho,sqrt(rho*eps)) in the non-relativistic code ! -!------------------------------------- - +! In this case no "solver" is required, but we do need to call the +! equation of state to get the pressure +!+ +!----------------------------------------------------------------------------- subroutine cons2prim_everything(npart,xyzh,vxyzu,dvdx,rad,eos_vars,radprop,& Bevol,Bxyz,dustevol,dustfrac,alphaind) use part, only:isdead_or_accreted,massoftype,igas,rhoh,igasP,iradP,iradxi,ics,imu,iX,iZ,& @@ -347,7 +369,6 @@ subroutine cons2prim_everything(npart,xyzh,vxyzu,dvdx,rad,eos_vars,radprop,& endif endif - end subroutine cons2prim_everything end module cons2prim diff --git a/src/main/cullendehnen.f90 b/src/main/cullendehnen.f90 index 5ebd2e7c9..b8c5d2946 100644 --- a/src/main/cullendehnen.f90 +++ b/src/main/cullendehnen.f90 @@ -42,6 +42,12 @@ pure real function get_alphaloc(divvdti,spsoundi,hi,xi_limiter,alphamin,alphamax end function get_alphaloc +!------------------------------------------------------------------------------- +!+ +! return the xi_limiter function used in the Cullen & Dehnen switch +! based on the spatial velocity gradients +!+ +!------------------------------------------------------------------------------- pure real function xi_limiter(dvdx) real(kind=4), intent(in) :: dvdx(9) real :: dvxdx,dvxdy,dvxdz,dvydx,dvydy,dvydz,dvzdx,dvzdy,dvzdz diff --git a/src/main/eos_HIIR.f90 b/src/main/eos_HIIR.f90 index 315d97734..a0ec8b61c 100644 --- a/src/main/eos_HIIR.f90 +++ b/src/main/eos_HIIR.f90 @@ -10,7 +10,7 @@ module eos_HIIR ! ! :References: None ! -! :Owner: Yrisch +! :Owner: Yann Bernard ! ! :Runtime parameters: None ! diff --git a/src/main/evolve.F90 b/src/main/evolve.F90 index a5bf47b1d..c7f63d6b1 100644 --- a/src/main/evolve.F90 +++ b/src/main/evolve.F90 @@ -306,7 +306,7 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag) if (iH2R > 0 .and. id==master) then istepHII = 1 - if(ind_timesteps) then + if (ind_timesteps) then istepHII = 2**nbinmax/HIIuprate if (istepHII==0) istepHII = 1 endif diff --git a/src/main/initial.F90 b/src/main/initial.F90 index 7b8f26032..d050a7f8c 100644 --- a/src/main/initial.F90 +++ b/src/main/initial.F90 @@ -19,11 +19,11 @@ module initial ! damping, densityforce, deriv, dim, dust, dust_formation, ! einsteintk_utils, energies, eos, evwrite, extern_gr, externalforces, ! fastmath, fileutils, forcing, growth, inject, io, io_summary, -! krome_interface, linklist, metric_tools, mf_write, mpibalance, -! mpidomain, mpimemory, mpitree, mpiutils, nicil, nicil_sup, omputils, -! options, part, partinject, porosity, ptmass, radiation_utils, -! readwrite_dumps, readwrite_infile, subgroup, timestep, timestep_ind, -! timestep_sts, timing, tmunu2grid, units, writeheader +! krome_interface, linklist, metric, metric_et_utils, metric_tools, +! mf_write, mpibalance, mpidomain, mpimemory, mpitree, mpiutils, nicil, +! nicil_sup, omputils, options, part, partinject, porosity, ptmass, +! radiation_utils, readwrite_dumps, readwrite_infile, subgroup, timestep, +! timestep_ind, timestep_sts, timing, tmunu2grid, units, writeheader ! implicit none @@ -41,7 +41,7 @@ module initial !+ !---------------------------------------------------------------- subroutine initialise() - use dim, only:mpi + use dim, only:mpi,gr use io, only:fatal,die,id,master,nprocs,ievfile #ifdef FINVSQRT use fastmath, only:testsqrt @@ -56,6 +56,8 @@ subroutine initialise() use cpuinfo, only:print_cpuinfo use checkoptions, only:check_compile_time_settings use readwrite_dumps, only:init_readwrite_dumps + use metric, only:metric_type + use metric_et_utils, only:read_tabulated_metric,gridinit integer :: ierr ! !--write 'PHANTOM' and code version @@ -99,6 +101,13 @@ subroutine initialise() !--initialise MPI domains ! call init_domains(nprocs) +! +!--initialise metric if tabulated +! + if (gr .and. metric_type=='et') then + call read_tabulated_metric('tabuled_metric.dat',ierr) + if (ierr == 0) gridinit = .true. + endif call init_readwrite_dumps() diff --git a/src/main/inject_asteroidwind.f90 b/src/main/inject_randomwind.f90 similarity index 66% rename from src/main/inject_asteroidwind.f90 rename to src/main/inject_randomwind.f90 index 5e61737fe..c55eea947 100644 --- a/src/main/inject_asteroidwind.f90 +++ b/src/main/inject_randomwind.f90 @@ -6,7 +6,7 @@ !--------------------------------------------------------------------------! module inject ! -! Injection module for wind from an orbiting asteroid, as used +! Injection module for wind from an orbiting body, as used ! in Trevascus et al. (2021) ! ! :References: @@ -15,9 +15,11 @@ module inject ! :Owner: David Liptai ! ! :Runtime parameters: -! - mdot : *mass injection rate in grams/second* -! - mdot_type : *injection rate (0=const, 1=cos(t), 2=r^(-2))* -! - vlag : *percentage lag in velocity of wind* +! - delta_theta : *standard deviation for the gaussion distribution (random_type=1)* +! - mdot : *mass injection rate with unit, e.g. 1e8*g/s, 1e-7M_s/yr* +! - mdot_type : *injection rate (0=const, 1=cos(t), 2=r^(-2))* +! - random_type : *random position on the surface, 0 for random, 1 for gaussian* +! - vlag : *percentage lag in velocity of wind* ! ! :Dependencies: binaryutils, externalforces, infile_utils, io, options, ! part, partinject, physcon, random, units @@ -25,8 +27,9 @@ module inject use io, only:error use physcon, only:pi implicit none - character(len=*), parameter, public :: inject_type = 'asteroidwind' - real, public :: mdot = 5.e8 ! mass injection rate in grams/second + character(len=*), parameter, public :: inject_type = 'randomwind' + character(len=20), public :: mdot_str = "5.e8*g/s" + real, public :: mdot = 1.e8 ! mass injection rate in grams/second public :: init_inject,inject_particles,write_options_inject,read_options_inject,& set_default_options_inject,update_injected_par @@ -35,7 +38,12 @@ module inject real :: npartperorbit = 1000. ! particle injection rate in particles per orbit real :: vlag = 0.0 ! percentage lag in velocity of wind - integer :: mdot_type = 2 ! injection rate (0=const, 1=cos(t), 2=r^(-2)) + integer :: mdot_type = 0 ! injection rate (0=const, 1=cos(t), 2=r^(-2)) + integer :: random_type = 0 ! random position on the surface, 0 for random, 1 for gaussian + real :: delta_theta = 0.5 ! standard deviation for the gaussion distribution (random_type=1) + real :: have_injected = 0. + real :: t_old = 0. + real :: r_ref = 1. contains !----------------------------------------------------------------------- @@ -61,8 +69,8 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& use part, only:nptmass,massoftype,igas,hfact,ihsoft use partinject, only:add_or_update_particle use physcon, only:twopi,gg,kboltz,mass_proton_cgs - use random, only:get_random_pos_on_sphere - use units, only:umass, utime + use random, only:get_random_pos_on_sphere, get_gaussian_pos_on_sphere + use units, only:in_code_units use options, only:iexternalforce use externalforces,only:mass1 use binaryutils, only:get_orbit_bits @@ -71,18 +79,17 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& integer, intent(inout) :: npart, npart_old integer, intent(inout) :: npartoftype(:) real, intent(out) :: dtinject + integer :: ierr real, dimension(3) :: xyz,vxyz,r1,r2,v2,vhat,v1 integer :: i,ipart,npinject,seed,pt - real :: dmdt,rasteroid,h,u,speed,inject_this_step + real :: dmdt,rbody,h,u,speed,inject_this_step real :: m1,m2,r real :: dt - real, save :: have_injected,t_old - real, save :: semia - if (nptmass < 2 .and. iexternalforce == 0) & - call fatal('inject_asteroidwind','not enough point masses for asteroid wind injection') + if (nptmass < 1 .and. iexternalforce == 0) & + call fatal('inject_randomwind','not enough point masses for random wind injection') if (nptmass > 2) & - call fatal('inject_asteroidwind','too many point masses for asteroid wind injection') + call fatal('inject_randomwind','too many point masses for random wind injection') if (nptmass == 2) then pt = 2 @@ -97,7 +104,7 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& endif r2 = xyzmh_ptmass(1:3,pt) - rasteroid = xyzmh_ptmass(ihsoft,pt) + rbody = xyzmh_ptmass(ihsoft,pt) m2 = xyzmh_ptmass(4,pt) v2 = vxyz_ptmass(1:3,pt) @@ -109,7 +116,8 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& ! ! Add any dependency on radius to mass injection rate (and convert to code units) ! - dmdt = mdot*mdot_func(r,semia)/(umass/utime) ! Use semi-major axis as r_ref + mdot = in_code_units(mdot_str,ierr) + dmdt = mdot*mdot_func(r,r_ref) ! r_ref is the radius for which mdot_fund = mdot ! !-- How many particles do we need to inject? @@ -120,7 +128,7 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& else ! Calculate how many extra particles from previous step to now dt = time - t_old - inject_this_step = dt*mdot/massoftype(igas)/(umass/utime) + inject_this_step = dt*dmdt/massoftype(igas) npinject = max(0, int(0.5 + have_injected + inject_this_step - npartoftype(igas) )) @@ -130,14 +138,13 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& endif ! - !-- Randomly inject particles around the asteroids outer 'radius'. - ! Only inject them on the side that is facing the central sink + !-- Randomly inject particles around the body's outer 'radius'. ! do i=1,npinject - xyz = r2 + rasteroid*get_random_pos_on_sphere(seed) + xyz = r2 + rbody*get_pos_on_sphere(seed, delta_theta) vxyz = (1.-vlag/100)*speed*vhat u = 0. ! setup is isothermal so utherm is not stored - h = hfact*(rasteroid/2.) + h = hfact*(rbody/2.) ipart = npart + 1 call add_or_update_particle(igas,xyz,vxyz,h,u,ipart,npart,npartoftype,xyzh,vxyzu) enddo @@ -173,6 +180,26 @@ real function mdot_func(r,r_ref) end function mdot_func +!----------------------------------------------------------------------- +!+ +! Returns a random location on a sperical surface +!+ +!----------------------------------------------------------------------- +function get_pos_on_sphere(iseed, delta_theta) result(dx) + use random, only:get_random_pos_on_sphere, get_gaussian_pos_on_sphere + integer, intent(inout) :: iseed + real, intent(inout) :: delta_theta + real :: dx(3) + + select case (random_type) + case(1) + dx = get_gaussian_pos_on_sphere(iseed, delta_theta) + case(0) + dx = get_random_pos_on_sphere(iseed) + end select + +end function get_pos_on_sphere + !----------------------------------------------------------------------- !+ ! Writes input options to the input file. @@ -182,11 +209,18 @@ subroutine write_options_inject(iunit) use infile_utils, only:write_inopt integer, intent(in) :: iunit - call write_inopt(mdot,'mdot','mass injection rate in grams/second',iunit) + call write_inopt(mdot_str,'mdot','mass injection rate with unit, e.g. 1e8*g/s, 1e-7M_s/yr',iunit) call write_inopt(npartperorbit,'npartperorbit',& 'particle injection rate in particles/binary orbit',iunit) call write_inopt(vlag,'vlag','percentage lag in velocity of wind',iunit) - call write_inopt(mdot_type,'mdot_type','injection rate (0=const, 1=cos(t), 2=r^(-2))',iunit) + call write_inopt(mdot_type,'mdot_type','injection rate (0=const, 2=r^(-2))',iunit) + if (mdot_type==2) then + call write_inopt(r_ref,'r_ref','radius at whieh Mdot=mdot for 1/r^2 injection type',iunit) + endif + call write_inopt(random_type, 'random_type', 'random position on the surface, 0 for random, 1 for gaussian', iunit) + if (random_type==1) then + call write_inopt(delta_theta, 'delta_theta', 'standard deviation for the gaussion distribution', iunit) + endif end subroutine write_options_inject @@ -206,9 +240,9 @@ subroutine read_options_inject(name,valstring,imatch,igotall,ierr) imatch = .true. select case(trim(name)) case('mdot') - read(valstring,*,iostat=ierr) mdot + read(valstring,'(A)',iostat=ierr) mdot_str ngot = ngot + 1 - if (mdot < 0.) call fatal(label,'mdot < 0 in input options') + ! if (mdot < 0.) call fatal(label,'mdot < 0 in input options') case('npartperorbit') read(valstring,*,iostat=ierr) npartperorbit ngot = ngot + 1 @@ -219,6 +253,15 @@ subroutine read_options_inject(name,valstring,imatch,igotall,ierr) case('mdot_type') read(valstring,*,iostat=ierr) mdot_type ngot = ngot + 1 + case('r_ref') + read(valstring,*,iostat=ierr) r_ref + ngot = ngot + 1 + case('random_type') + read(valstring,*,iostat=ierr) random_type + ngot = ngot + 1 + case('delta_theta') + read(valstring,*,iostat=ierr) delta_theta + ngot = ngot + 1 case default imatch = .false. end select diff --git a/src/main/metric_et.f90 b/src/main/metric_et.f90 index ce133ea83..cd93c9c3b 100644 --- a/src/main/metric_et.f90 +++ b/src/main/metric_et.f90 @@ -6,7 +6,7 @@ !--------------------------------------------------------------------------! module metric ! -! None +! Generic module for a tabulated metric, e.g. from the Einstein Toolkit ! ! :References: None ! @@ -14,7 +14,7 @@ module metric ! ! :Runtime parameters: None ! -! :Dependencies: einsteintk_utils, eos_shen, infile_utils +! :Dependencies: metric_et_utils, table_utils ! implicit none character(len=*), parameter :: metric_type = 'et' @@ -28,44 +28,54 @@ module metric !---------------------------------------------------------------- !+ ! Compute the metric tensor in both covariant (gcov) and -! contravariant (gcon) form +! contravariant (gcon) form. Here we merely interpolate +! these values from the global grid. !+ !---------------------------------------------------------------- pure subroutine get_metric_cartesian(position,gcov,gcon,sqrtg) - use einsteintk_utils, only:gridinit + use metric_et_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 + integer :: ierr ! 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 + ! Perform trilinear interpolation if ( .not. gridinit) then + ierr = 1 ! 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. + if (ierr /= 0) then + 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. 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 +!----------------------------------------------------------------------- +!+ +! dummy routine to get the metric in spherical coordinates (not used) +!+ +!----------------------------------------------------------------------- pure subroutine get_metric_spherical(position,gcov,gcon,sqrtg) real, intent(in) :: position(3) real, intent(out) :: gcov(0:3,0:3) @@ -95,19 +105,39 @@ pure subroutine get_metric_spherical(position,gcov,gcon,sqrtg) end subroutine get_metric_spherical +!----------------------------------------------------------------------- +!+ +! cartesian metric derivatives, interpolates the derivatives from +! the grid +!+ +!----------------------------------------------------------------------- pure subroutine metric_cartesian_derivatives(position,dgcovdx, dgcovdy, dgcovdz) - use einsteintk_utils, only:gridinit + use metric_et_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) + integer :: ierr + if (.not. gridinit) then - dgcovdx = 0. - dgcovdy = 0. - dgcovdz = 0. + ierr = 1 + if (ierr /= 0) then + dgcovdx = 0. + dgcovdy = 0. + dgcovdz = 0. + else + ! gridinit = .true. + call interpolate_metric_derivs(position,dgcovdx,dgcovdy,dgcovdz) + endif else call interpolate_metric_derivs(position,dgcovdx,dgcovdy,dgcovdz) endif + end subroutine metric_cartesian_derivatives +!----------------------------------------------------------------------- +!+ +! dummy routine for spherical metric derivatives, not used +!+ +!----------------------------------------------------------------------- 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 @@ -127,6 +157,11 @@ pure subroutine metric_spherical_derivatives(position,dgcovdr, dgcovdtheta, dgco end subroutine metric_spherical_derivatives +!----------------------------------------------------------------------- +!+ +! dummy routine to convert cartesian to spherical coordinates +!+ +!----------------------------------------------------------------------- pure subroutine cartesian2spherical(xcart,xspher) real, intent(in) :: xcart(3) real, intent(out) :: xspher(3) @@ -142,6 +177,7 @@ pure subroutine cartesian2spherical(xcart,xspher) phi = atan2(y,x) xspher = (/r,theta,phi/) + end subroutine cartesian2spherical !----------------------------------------------------------------------- @@ -150,10 +186,11 @@ end subroutine cartesian2spherical !+ !----------------------------------------------------------------------- subroutine write_options_metric(iunit) - use infile_utils, only:write_inopt + !use infile_utils, only:write_inopt integer, intent(in) :: iunit - write(iunit,"(/,a)") '# There are no options relating to the '//trim(metric_type)//' metric' + !write(iunit,"(/,a)") '# There are no options relating to the '//trim(metric_type)//' metric' + !call write_inopt(metric_file,'metric_file','file from which to read tabulated metric (blank if used with einsteintk)',iunit) end subroutine write_options_metric @@ -166,9 +203,17 @@ 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. + !integer, save :: ngot = 0 + + select case(trim(name)) + !case('metric_file') + ! read(valstring,*,iostat=ierr) metric_file + ! ngot = ngot + 1 + case default + imatch = .false. + end select + !igotall = (ngot >= 1) + igotall = .true. end subroutine read_options_metric @@ -177,12 +222,11 @@ 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 + use table_utils, only:linear_interpolator_one_d + use metric_et_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 @@ -197,19 +241,9 @@ pure subroutine interpolate_metric(position,gcov,gcon,sqrtg) ! 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) @@ -286,12 +320,16 @@ pure subroutine interpolate_metric(position,gcov,gcon,sqrtg) sqrtg = interptmp(7) endif - end subroutine interpolate_metric +!----------------------------------------------------------------------- +!+ +! Interpolates derivatives of the metric from the grid to the position +!+ +!----------------------------------------------------------------------- pure subroutine interpolate_metric_derivs(position,dgcovdx, dgcovdy, dgcovdz) - use eos_shen, only:linear_interpolator_one_d - use einsteintk_utils, only:metricderivsgrid, dxgrid,gridorigin + use table_utils, only:linear_interpolator_one_d + use metric_et_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 @@ -300,13 +338,6 @@ pure subroutine interpolate_metric_derivs(position,dgcovdx, dgcovdy, dgcovdz) 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) @@ -383,13 +414,15 @@ pure subroutine interpolate_metric_derivs(position,dgcovdx, dgcovdy, dgcovdz) enddo enddo - - - end subroutine interpolate_metric_derivs +!----------------------------------------------------------------------- +!+ +! Utility routine to get the lower grid neighbours of a position +!+ +!----------------------------------------------------------------------- pure subroutine get_grid_neighbours(position,dx,xlower,ylower,zlower) - use einsteintk_utils, only:gridorigin + use metric_et_utils, only:gridorigin real, intent(in) :: position(3) real, intent(in) :: dx(3) integer, intent(out) :: xlower,ylower,zlower @@ -409,8 +442,7 @@ pure subroutine get_grid_neighbours(position,dx,xlower,ylower,zlower) ylower = ylower + 1 zlower = zlower + 1 - end subroutine get_grid_neighbours - end module metric + diff --git a/src/main/metric_et_utils.f90 b/src/main/metric_et_utils.f90 new file mode 100644 index 000000000..a3c3bebf5 --- /dev/null +++ b/src/main/metric_et_utils.f90 @@ -0,0 +1,193 @@ +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2024 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.github.io/ ! +!--------------------------------------------------------------------------! +module metric_et_utils +! +! Utilities for handling tabulated metrics from the Einstein Toolkit +! +! :References: Magnall et al. (2023), Phys. Rev. D 108, 103534 +! +! :Owner: Daniel Price +! +! :Runtime parameters: None +! +! :Dependencies: None +! + implicit none + + real, allocatable :: gcovgrid(:,:,:,:,:) + real, allocatable :: gcongrid(:,:,:,:,:) + real, allocatable :: sqrtggrid(:,:,:) + real, allocatable :: metricderivsgrid(:,:,:,:,:,:) + real :: dxgrid(3), gridorigin(3) + integer :: gridsize(3) + logical :: gridinit = .false. + + ! Declaration of grid limits and dimensions + integer, public :: nx,ny,nz + real, parameter :: xmin = -10.0, xmax = 10.0 + real, parameter :: ymin = -10.0, ymax = 10.0 + real, parameter :: zmin = -10.0, zmax = 10.0 + real, parameter :: mass = 1.0 ! Mass of the central object + +contains + +!--------------------------------------------------------------- +!+ +! allocate memory for the metric grid +!+ +!--------------------------------------------------------------- +subroutine allocate_grid(nxin,nyin,nzin,dx,dy,dz,originx,originy,originz) + integer, intent(in) :: nxin,nyin,nzin + real, intent(in) :: dx,dy,dz,originx,originy,originz + + nx = nxin + ny = nyin + nz = nzin + 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)) + + !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)) + +end subroutine allocate_grid + +!--------------------------------------------------------------- +!+ +! initialise a metric grid with a uniform grid +! (currently size is hardwired but just for testing...) +!+ +!--------------------------------------------------------------- +subroutine initialize_grid() + ! Local variable declarations + real :: dx, dy, dz, x0(3) + + nx = 100 + ny = 100 + nz = 100 + + ! Calculate the step size in each direction + dx = (xmax - xmin) / (nx - 1) + dy = (ymax - ymin) / (ny - 1) + dz = (zmax - zmin) / (nz - 1) + + x0 = [0.,0.,0.] + call allocate_grid(nx,ny,nz,dx,dy,dz,x0(1),x0(2),x0(3)) + + gridinit = .true. + +end subroutine initialize_grid + +!--------------------------------------------------------------- +!+ +! print information about the metric grid +!+ +!--------------------------------------------------------------- +subroutine print_metric_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_metric_grid + +!--------------------------------------------------------------- +!+ +! write tabulated metric to file +!+ +!--------------------------------------------------------------- +subroutine write_tabulated_metric(metric_file, ierr) + character(len=*), intent(in) :: metric_file + integer, intent(out) :: ierr + integer :: iunit + + ! Open the file for writing + open(newunit=iunit,file=metric_file,status='replace',form='unformatted',action='write',iostat=ierr) + if (ierr /= 0) then + ierr = 1 + return + endif + + ! Write the dimensions of the grid + write(iunit) gridsize + + ! Write the grid origin and spacing + write(iunit) gridorigin + write(iunit) dxgrid + + ! Write the metric values to the file + write(iunit) gcovgrid + write(iunit) gcongrid + write(iunit) sqrtggrid + write(iunit) metricderivsgrid + + ! Close the file + close(iunit) + ierr = 0 + +end subroutine write_tabulated_metric + +!--------------------------------------------------------------- +!+ +! read tabulated metric from file +!+ +!--------------------------------------------------------------- +subroutine read_tabulated_metric(metric_file, ierr) + character(len=*), intent(in) :: metric_file + integer, intent(out) :: ierr + integer :: iunit + + ! Open the file for reading + open(newunit=iunit,file=metric_file,status='old',form='unformatted',action='read',iostat=ierr) + if (ierr /= 0) return + + ! Read the dimensions of the grid + read(iunit) gridsize + + ! Read the grid origin and spacing + read(iunit) gridorigin + read(iunit) dxgrid + + nx = gridsize(1) + ny = gridsize(2) + nz = gridsize(3) + + call allocate_grid(nx,ny,nz,& + dxgrid(1),dxgrid(2),dxgrid(3),& + gridorigin(1),gridorigin(2),gridorigin(3)) + + ! Read the metric values from the file + read(iunit) gcovgrid + read(iunit) gcongrid + read(iunit) sqrtggrid + read(iunit) metricderivsgrid + + gridinit = .true. + + ! Close the file + close(iunit) + ierr = 0 + +end subroutine read_tabulated_metric + +end module metric_et_utils diff --git a/src/main/metric_tools.F90 b/src/main/metric_tools.f90 similarity index 81% rename from src/main/metric_tools.F90 rename to src/main/metric_tools.f90 index 8fd54fdf0..b0cf56c3f 100644 --- a/src/main/metric_tools.F90 +++ b/src/main/metric_tools.f90 @@ -29,24 +29,22 @@ module metric_tools icoord_cartesian = 1, & ! Cartesian coordinates icoord_spherical = 2 ! Spherical coordinates -!--- List of metrics + !--- List of metrics integer, public, parameter :: & imet_minkowski = 1, & ! Minkowski metric imet_schwarzschild = 2, & ! Schwarzschild 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) + !--- Choice of coordinate system + ! (When using this with PHANTOM, it should always be set to cartesian) integer, public, parameter :: icoordinate = icoord_cartesian -!--- Choice for contravariant metric -! false -> use analytic contravariant metric -! true -> invert the covariant metric + !--- Choice for contravariant metric + ! false -> use analytic contravariant metric + ! true -> invert the covariant metric logical, private, parameter :: useinv4x4 = .true. -!------------------------------------------------------------------------------- - public :: get_metric, get_metric_derivs, print_metricinfo, init_metric, pack_metric, unpack_metric public :: pack_metricderivs public :: imetric @@ -63,8 +61,8 @@ module metric_tools !+ !------------------------------------------------------------------------------- pure subroutine get_metric(position,gcov,gcon,sqrtg) - use metric, only: get_metric_cartesian,get_metric_spherical,cartesian2spherical - use inverse4x4, only: inv4x4 + use metric, only:get_metric_cartesian,get_metric_spherical,cartesian2spherical + use inverse4x4, only:inv4x4 real, intent(in) :: position(3) real, intent(out) :: gcov(0:3,0:3), gcon(0:3,0:3), sqrtg real :: det @@ -95,8 +93,8 @@ end subroutine get_metric ! of metric. !+ !------------------------------------------------------------------------------- -subroutine get_metric_derivs(position,dgcovdx1, dgcovdx2, dgcovdx3) - use metric, only: metric_cartesian_derivatives, metric_spherical_derivatives, imetric +subroutine get_metric_derivs(position,dgcovdx1,dgcovdx2,dgcovdx3) + use metric, only:metric_cartesian_derivatives,metric_spherical_derivatives,imetric real, intent(in) :: position(3) real, intent(out) :: dgcovdx1(0:3,0:3), dgcovdx2(0:3,0:3), dgcovdx3(0:3,0:3) @@ -114,7 +112,7 @@ end subroutine get_metric_derivs ! Numerical derivatives of the covariant metric tensor !+ !------------------------------------------------------------------------------- -pure subroutine numerical_metric_derivs(position,dgcovdx, dgcovdy, dgcovdz) +pure subroutine numerical_metric_derivs(position,dgcovdx,dgcovdy,dgcovdz) real, intent(in) :: position(3) real, intent(out), dimension(0:3,0:3) :: dgcovdx,dgcovdy,dgcovdz real :: gblah(0:3,0:3), temp(3), gplus(0:3,0:3),gminus(0:3,0:3),dx,dy,dz,di,sqrtgblag @@ -150,29 +148,10 @@ pure subroutine numerical_metric_derivs(position,dgcovdx, dgcovdy, dgcovdz) end subroutine numerical_metric_derivs !------------------------------------------------------------------------------- - -! This is not being used at the moment... -!-- Do a coordinate transformation of a 4x4 rank-2 tensor with both indices down -! subroutine tensortransform_dd(position,T_old,T_new) -! use metric, only: get_jacobian -! real, intent(in), dimension(3) :: position -! real, intent(in), dimension(0:3,0:3) :: T_old -! real, intent(out), dimension(0:3,0:3) :: T_new -! real, dimension(0:3,0:3) :: dxdx -! integer :: i,j,k,l -! call get_jacobian(position,dxdx) -! T_new = 0. -! do i=0,3 -! do j=0,3 -! do k=0,3 -! do l=0,3 -! T_new(i,j) = T_new(i,j)+dxdx(k,i)*dxdx(l,j)*T_old(k,l) -! enddo -! enddo -! enddo -! enddo -! end subroutine tensortransform_dd - +!+ +! print the metric type +!+ +!------------------------------------------------------------------------------- subroutine print_metricinfo(iprint) use metric, only:metric_type integer, intent(in) :: iprint @@ -181,6 +160,11 @@ subroutine print_metricinfo(iprint) end subroutine print_metricinfo +!------------------------------------------------------------------------------- +!+ +! initialise arrays for the metric and metric derivatives +!+ +!------------------------------------------------------------------------------- subroutine init_metric(npart,xyzh,metrics,metricderivs) integer, intent(in) :: npart real, intent(in) :: xyzh(:,:) @@ -188,7 +172,6 @@ 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) @@ -209,9 +192,11 @@ subroutine init_metric(npart,xyzh,metrics,metricderivs) end subroutine init_metric -! -!--- Subroutine to pack the metric (cov and con) into a single array -! +!------------------------------------------------------------------------------- +!+ +! subroutine to pack the metric into a 4x4x2 array +!+ +!------------------------------------------------------------------------------- pure subroutine pack_metric(xyz,metrici) real, intent(in) :: xyz(3) real, intent(out) :: metrici(:,:,:) @@ -221,6 +206,11 @@ pure subroutine pack_metric(xyz,metrici) end subroutine pack_metric +!------------------------------------------------------------------------------- +!+ +! subroutine to pack the metric derivatives into a 4x4x3 array +!+ +!------------------------------------------------------------------------------- subroutine pack_metricderivs(xyzi,metricderivsi) real, intent(in) :: xyzi(3) real, intent(out) :: metricderivsi(0:3,0:3,3) @@ -229,24 +219,19 @@ subroutine pack_metricderivs(xyzi,metricderivsi) end subroutine pack_metricderivs -! -!--- Subroutine to return metric/components from metrici array -! +!------------------------------------------------------------------------------- +!+ +! Subroutine to return metric/components from metrici array +!+ +!------------------------------------------------------------------------------- pure subroutine unpack_metric(metrici,gcov,gcon,gammaijdown,gammaijUP,alpha,betadown,betaUP) -#ifdef FINVSQRT - use fastmath, only:finvsqrt -#endif real, intent(in), dimension(0:3,0:3,2) :: metrici real, intent(out), dimension(0:3,0:3), optional :: gcov,gcon real, intent(out), dimension(1:3,1:3), optional :: gammaijdown,gammaijUP real, intent(out), optional :: alpha,betadown(3),betaUP(3) integer :: i -#ifdef FINVSQRT - if (present(alpha)) alpha = finvsqrt(-metrici(0,0,2)) -#else if (present(alpha)) alpha = sqrt(-1./metrici(0,0,2)) -#endif if (present(betaUP)) betaUP = metrici(0,1:3,2) * (-1./metrici(0,0,2)) ! = gcon(0,1:3)*alpha**2 @@ -264,6 +249,4 @@ pure subroutine unpack_metric(metrici,gcov,gcon,gammaijdown,gammaijUP,alpha,beta end subroutine unpack_metric - - end module metric_tools diff --git a/src/main/ptmass.F90 b/src/main/ptmass.F90 index a5362c8a0..2f79f2738 100644 --- a/src/main/ptmass.F90 +++ b/src/main/ptmass.F90 @@ -1639,6 +1639,11 @@ subroutine ptmass_create(nptmass,npart,itest,xyzh,vxyzu,fxyzu,fext,divcurlv,pote end subroutine ptmass_create +!------------------------------------------------------------------------- +!+ +! subroutine to create a bundh of star "seeds" inside a sink particle +!+ +!------------------------------------------------------------------------- subroutine ptmass_create_seeds(nptmass,itest,xyzmh_ptmass,linklist_ptmass,time) use part, only:itbirth,ihacc use random, only:ran2 @@ -1675,7 +1680,13 @@ subroutine ptmass_create_seeds(nptmass,itest,xyzmh_ptmass,linklist_ptmass,time) end subroutine ptmass_create_seeds -subroutine ptmass_create_stars(nptmass,itest,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,fxyz_ptmass_sinksink,linklist_ptmass,time) +!------------------------------------------------------------------------- +!+ +! subroutine to create a bundh of stars inside a sink particle +!+ +!------------------------------------------------------------------------- +subroutine ptmass_create_stars(nptmass,itest,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,& + fxyz_ptmass_sinksink,linklist_ptmass,time) use dim, only:maxptmass use physcon, only:solarm,pi use io, only:iprint,iverbose @@ -1923,22 +1934,34 @@ subroutine merge_sinks(time,nptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,linklis end subroutine merge_sinks +!----------------------------------------------------------------------- +!+ +! helper routine for managing the sink particle linked list +!+ +!----------------------------------------------------------------------- subroutine ptmass_endsize_lklist(i,k,n,linklist_ptmass) integer, intent(in) :: linklist_ptmass(:) integer, intent(in) :: i integer, intent(out) :: k,n integer :: l,g - g=i + + g = i n = 0 do while (g>0) l = g g = linklist_ptmass(l) n = n + 1 enddo - k=l -end subroutine ptmass_endsize_lklist + k = l +end subroutine ptmass_endsize_lklist +!----------------------------------------------------------------------- +!+ +! Swap between leapfrog and 4th order forward sympletic integrator +! for evolving sink particles +!+ +!----------------------------------------------------------------------- subroutine set_integration_precision if (use_fourthorder) then @@ -2197,7 +2220,7 @@ subroutine write_options_ptmass(iunit) integer, intent(in) :: iunit write(iunit,"(/,a)") '# options controlling sink particles' - call write_inopt(isink_potential,'isink_potential','sink potential(0=1/r,1=surf)',iunit) + call write_inopt(isink_potential,'isink_potential','sink potential (0=1/r,1=surf)',iunit) if (gravity) then call write_inopt(icreate_sinks,'icreate_sinks','allow automatic sink particle creation',iunit) if (icreate_sinks > 0) then @@ -2219,8 +2242,10 @@ subroutine write_options_ptmass(iunit) endif call write_inopt(h_soft_sinksink,'h_soft_sinksink','softening length between sink particles',iunit) call write_inopt(f_acc,'f_acc','particles < f_acc*h_acc accreted without checks',iunit) - call write_inopt(r_merge_uncond,'r_merge_uncond','sinks will unconditionally merge within this separation',iunit) - call write_inopt(r_merge_cond,'r_merge_cond','sinks will merge if bound within this radius',iunit) + if (gravity .and. icreate_sinks > 0) then + call write_inopt(r_merge_uncond,'r_merge_uncond','sinks will unconditionally merge within this separation',iunit) + call write_inopt(r_merge_cond,'r_merge_cond','sinks will merge if bound within this radius',iunit) + endif if (use_regnbody) then call write_inopt(use_regnbody, 'use_regnbody', 'allow subgroup integration method', iunit) call write_inopt(r_neigh, 'r_neigh', 'searching radius to detect subgroups', iunit) @@ -2321,7 +2346,7 @@ subroutine read_options_ptmass(name,valstring,imatch,igotall,ierr) imatch = .false. end select - if (icreate_sinks ==2) store_ll_ptmass = .true. + if (icreate_sinks == 2) store_ll_ptmass = .true. !--make sure we have got all compulsory options (otherwise, rewrite input file) if (icreate_sinks > 0) then @@ -2331,5 +2356,5 @@ subroutine read_options_ptmass(name,valstring,imatch,igotall,ierr) endif end subroutine read_options_ptmass -!----------------------------------------------------------------------- + end module ptmass diff --git a/src/main/random.f90 b/src/main/random.f90 index b5fc3bd88..dd2ba97c1 100644 --- a/src/main/random.f90 +++ b/src/main/random.f90 @@ -19,7 +19,8 @@ module random ! implicit none public :: ran2,get_random,rayleigh_deviate - public :: get_random_pos_on_sphere,gauss_random,divide_unit_seg + public :: get_random_pos_on_sphere,get_gaussian_pos_on_sphere + public :: gauss_random,divide_unit_seg real, parameter :: pi = 4.*atan(1.) private @@ -149,6 +150,31 @@ function get_random_pos_on_sphere(iseed) result(dx) end function get_random_pos_on_sphere +!------------------------------------------------------------------------- +! +! get Gaussian random position on sphere +! +!------------------------------------------------------------------------- +function get_gaussian_pos_on_sphere(iseed, deltheta) result(dx) + integer, intent(inout) :: iseed + real, intent(in) :: deltheta + real :: phi,theta,sintheta,costheta,sinphi,cosphi,gauss_theta + real :: dx(3) + + phi = 2.*pi*(ran2(iseed) - 0.5) + gauss_theta = gauss_random(iseed) * deltheta + do while (abs(gauss_theta) > 1.0) + gauss_theta = gauss_random(iseed) * deltheta + enddo + theta = acos(gauss_theta) + sintheta = sin(theta) + costheta = cos(theta) + sinphi = sin(phi) + cosphi = cos(phi) + dx = (/sintheta*cosphi,sintheta*sinphi,costheta/) + +end function get_gaussian_pos_on_sphere + !------------------------------------------------------------------------- ! ! get random number from gaussian @@ -167,7 +193,11 @@ real function gauss_random(iseed) end function gauss_random - +!------------------------------------------------------------------------- +! +! divide a unit segment into nlengths segments of random length +! +!------------------------------------------------------------------------- subroutine divide_unit_seg(lengths,mindist,nlengths,iseed) use sortutils, only:indexx integer, intent(in) :: nlengths @@ -192,7 +222,6 @@ subroutine divide_unit_seg(lengths,mindist,nlengths,iseed) mindist = 0.01 ! we'll have stars less massive than 0.08 solarmasses but it will assure to never brick the sim... endif - do i=2,nlengths close = .true. do while (close) @@ -217,5 +246,4 @@ subroutine divide_unit_seg(lengths,mindist,nlengths,iseed) end subroutine divide_unit_seg - end module random diff --git a/src/main/readwrite_dumps_common.f90 b/src/main/readwrite_dumps_common.f90 index c3059a51f..f67ae3c05 100644 --- a/src/main/readwrite_dumps_common.f90 +++ b/src/main/readwrite_dumps_common.f90 @@ -16,7 +16,7 @@ module readwrite_dumps_common ! ! :Dependencies: boundary, boundary_dyn, checkconserved, dim, dump_utils, ! dust, dust_formation, eos, externalforces, fileutils, gitinfo, io, -! options, part, ptmass, setup_params, sphNGutils, timestep, units +! options, part, setup_params, sphNGutils, timestep, units ! use dump_utils, only:lenid implicit none diff --git a/src/main/readwrite_dumps_fortran.f90 b/src/main/readwrite_dumps_fortran.f90 index 40015c958..429c6ff32 100644 --- a/src/main/readwrite_dumps_fortran.f90 +++ b/src/main/readwrite_dumps_fortran.f90 @@ -19,7 +19,7 @@ module readwrite_dumps_fortran ! :Runtime parameters: None ! ! :Dependencies: boundary_dyn, dim, dump_utils, eos, io, memory, -! metric_tools, mpiutils, options, part, ptmass, readwrite_dumps_common, +! metric_tools, mpiutils, options, part, readwrite_dumps_common, ! sphNGutils, timestep ! use dump_utils, only:lenid,ndatatypes,i_int,i_int1,i_int2,i_int4,i_int8,& diff --git a/src/main/subgroup.f90 b/src/main/subgroup.f90 index 1e6e52e28..940e8e92a 100644 --- a/src/main/subgroup.f90 +++ b/src/main/subgroup.f90 @@ -11,7 +11,7 @@ module subgroup ! ! :References: Makkino et Aarseth 2002,Wang et al. 2020, Wang et al. 2021, Rantala et al. 2023 ! -! :Owner: Yrisch +! :Owner: Yann Bernard ! ! :Runtime parameters: None ! diff --git a/src/main/substepping.F90 b/src/main/substepping.F90 index 504e7b9d3..a0d8d4693 100644 --- a/src/main/substepping.F90 +++ b/src/main/substepping.F90 @@ -21,7 +21,7 @@ module substepping ! Tuckerman, Berne & Martyna (1992), J. Chem. Phys. 97, 1990-2001 ! Rantala + (2020) (2023),Chin (2007a) ! -! :Owner: Yrisch +! :Owner: Yann Bernard ! ! :Runtime parameters: None ! diff --git a/src/main/units.f90 b/src/main/units.f90 index 62cae18aa..8af9b5dc9 100644 --- a/src/main/units.f90 +++ b/src/main/units.f90 @@ -34,7 +34,7 @@ module units public :: set_units, set_units_extra, print_units public :: get_G_code, get_c_code, get_radconst_code, get_kbmh_code public :: c_is_unity, G_is_unity, in_geometric_units - public :: is_time_unit, is_length_unit + public :: is_time_unit, is_length_unit, is_mdot_unit public :: in_solarr, in_solarm, in_solarl contains @@ -226,6 +226,12 @@ subroutine select_unit(string,unit,ierr) unit = minutes case('s','sec','second','seconds') unit = seconds + case("g/s","grams/second","g/second","grams/s","g/sec","grams/sec") + unit = 1.d0/seconds + case("Ms/yr","M_s/yr","ms/yr","m_s/yr","Msun/yr","M_sun/yr","Msolar/yr",& + "M_solar/yr","Ms/year","M_s/year","ms/year","m_s/year","Msun/year",& + "M_sun/year","Msolar/year","M_solar/year") + unit = solarm/years case default ierr = 1 unit = 1.d0 @@ -285,6 +291,32 @@ logical function is_length_unit(string) end function is_length_unit +!------------------------------------------------------------------------------------ +!+ +! check if string is a unit of mdot +!+ +!------------------------------------------------------------------------------------ +logical function is_mdot_unit(string) + character(len=*), intent(in) :: string + character(len=len(string)) :: unitstr + real(kind=8) :: fac + integer :: ierr + + ierr = 0 + call get_unit_multiplier(string,unitstr,fac,ierr) + + select case(trim(unitstr)) + case("g/s","gram/second","g/second","gram/s","g/sec","gram/sec",& + "Ms/yr","M_s/yr","ms/yr","m_s/yr","Msun/yr","M_sun/yr","Msolar/yr",& + "M_solar/yr","Ms/year","M_s/year","ms/year","m_s/year","Msun/year",& + "M_sun/year","Msolar/year","M_solar/year") + is_mdot_unit = .true. + case default + is_mdot_unit = .false. + end select + +end function is_mdot_unit + !------------------------------------------------------------------------------------ !+ ! parse a string like "10.*days" or "10*au" and return the value in code units @@ -301,6 +333,8 @@ real function in_code_units(string,ierr) result(rval) rval = real(val/utime) elseif (is_length_unit(string) .and. ierr == 0) then rval = real(val/udist) + elseif (is_mdot_unit(string) .and. ierr == 0) then + rval = real(val/(umass/utime)) else rval = real(val) ! no unit conversion endif diff --git a/src/main/utils_gr.F90 b/src/main/utils_gr.f90 similarity index 80% rename from src/main/utils_gr.F90 rename to src/main/utils_gr.f90 index 479476ca6..1308aeef8 100644 --- a/src/main/utils_gr.F90 +++ b/src/main/utils_gr.f90 @@ -6,9 +6,9 @@ !--------------------------------------------------------------------------! module utils_gr ! -! None +! Utility routines for the GR code ! -! :References: None +! :References: Liptai & Price (2019), MNRAS 485, 819-842 ! ! :Owner: David Liptai ! @@ -42,15 +42,14 @@ pure real function dot_product_gr(vec1,vec2,gcov) dot_product_gr = dot_product_gr + dot_product(gcov(:,i),vec1(i)*vec2(:)) enddo - return end function dot_product_gr -!------------------------------------------------------------------------------- - +!---------------------------------------------------------------- +!+ +! Function to return U^0, the time component of the 4-velocity +!+ +!---------------------------------------------------------------- pure subroutine get_u0(gcov,v,U0,ierr) -#ifdef FINVSQRT - use fastmath, only:finvsqrt -#endif real, intent(in) :: gcov(0:3,0:3), v(1:3) real, intent(out) :: U0 integer, intent(out) :: ierr @@ -60,23 +59,19 @@ pure subroutine get_u0(gcov,v,U0,ierr) v4(0) = 1. v4(1:3) = v(1:3) vv = dot_product_gr(v4,v4,gcov) -#ifdef FINVSQRT - U0 = finvsqrt(-vv) -#else U0 = 1./sqrt(-vv) -#endif if (vv > 0.) ierr = 1 end subroutine get_u0 -!------------------------------------------------------------------------------- - +!---------------------------------------------------------------- +!+ +! Function to return V^i, the velocity of an Eulerian observer +!+ +!---------------------------------------------------------------- subroutine get_bigv(metrici,v,bigv,bigv2,alpha,lorentz) use metric_tools, only:unpack_metric use io, only:fatal -#ifdef FINVSQRT - use fastmath, only:finvsqrt -#endif real, intent(in) :: metrici(0:3,0:3,2),v(1:3) real, intent(out) :: bigv(1:3),bigv2,alpha,lorentz real :: betaUP(1:3),gammaijdown(1:3,1:3) @@ -85,16 +80,16 @@ subroutine get_bigv(metrici,v,bigv,bigv2,alpha,lorentz) bigv = (v + betaUP)/alpha bigv2 = dot_product_gr(bigv,bigv,gammaijdown) if (bigv2 > 1.) call fatal('get_bigv','velocity faster than speed of light -- bigv2',val=bigv2) -#ifdef FINVSQRT - lorentz = finvsqrt(1.-bigv2) -#else lorentz = 1./sqrt(1.-bigv2) -#endif end subroutine get_bigv -!------------------------------------------------------------------------------- - +!---------------------------------------------------------------- +!+ +! get density in the fluid rest frame (primitive dens) from +! the conserved density rho* (stored as the smoothing length) +!+ +!---------------------------------------------------------------- subroutine h2dens(dens,xyzh,metrici,v) use part, only: rhoh,massoftype,igas real, intent(in) :: xyzh(1:4),metrici(:,:,:),v(1:3) @@ -108,6 +103,12 @@ subroutine h2dens(dens,xyzh,metrici,v) end subroutine h2dens +!---------------------------------------------------------------- +!+ +! get density in the fluid rest frame (primitive dens) from +! the conserved density rho* +!+ +!---------------------------------------------------------------- subroutine rho2dens(dens,rho,position,metrici,v) use metric_tools, only:unpack_metric use io, only:error @@ -116,7 +117,6 @@ subroutine rho2dens(dens,rho,position,metrici,v) integer :: ierror real :: gcov(0:3,0:3), sqrtg, U0 - call unpack_metric(metrici,gcov=gcov) call get_sqrtg(gcov, sqrtg) call get_u0(gcov,v,U0,ierror) @@ -126,6 +126,12 @@ subroutine rho2dens(dens,rho,position,metrici,v) end subroutine rho2dens +!---------------------------------------------------------------- +!+ +! get terms required on the RHS of the geodesic equation +! in the form dp_i/dt = a_i, as described in Liptai & Price (2019) +!+ +!---------------------------------------------------------------- subroutine get_geodesic_accel(axyz,npart,vxyz,metrics,metricderivs) use metric_tools, only:unpack_metric integer, intent(in) :: npart @@ -157,6 +163,11 @@ subroutine get_geodesic_accel(axyz,npart,vxyz,metrics,metricderivs) end subroutine get_geodesic_accel +!---------------------------------------------------------------- +!+ +! get determininant of the 4-metric +!+ +!---------------------------------------------------------------- subroutine get_sqrtg(gcov, sqrtg) use metric, only: metric_type real, intent(in) :: gcov(0:3,0:3) @@ -205,6 +216,11 @@ subroutine get_sqrtg(gcov, sqrtg) end subroutine get_sqrtg +!---------------------------------------------------------------- +!+ +! get determininant of the 3-metric +!+ +!---------------------------------------------------------------- subroutine get_sqrt_gamma(gcov,sqrt_gamma) use metric, only: metric_type real, intent(in) :: gcov(0:3,0:3) @@ -235,18 +251,20 @@ subroutine get_sqrt_gamma(gcov,sqrt_gamma) else sqrt_gamma = 1. - endif - end subroutine get_sqrt_gamma +!---------------------------------------------------------------- +!+ +! add a Newtonian gravitational perturbation to the metric +!+ +!---------------------------------------------------------------- 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 @@ -257,29 +275,12 @@ subroutine perturb_metric(phi,gcovper,gcov) gcovper(3,3) = 1. endif - ! Set the pertubed metric based on the Bardeen formulation + ! Set the perturbed 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 -! real, intent(in) :: dens,position(1:3),v(1:3) -! real, intent(out):: rho -! real :: gcov(0:3,0:3), gcon(0:3,0:3), sqrtg, U0 -! -! call get_metric(position,gcov,gcon,sqrtg) -! call get_u0(gcov,v,U0) -! -! rho = sqrtg*U0*dens -! -! end subroutine dens2rho - -!------------------------------------------------------------------------------- - end module utils_gr diff --git a/src/main/utils_kepler.f90 b/src/main/utils_kepler.f90 index deb5de94b..071844eda 100644 --- a/src/main/utils_kepler.f90 +++ b/src/main/utils_kepler.f90 @@ -10,7 +10,7 @@ module utils_kepler ! ! :References: None ! -! :Owner: Yrisch +! :Owner: Yann Bernard ! ! :Runtime parameters: None ! diff --git a/src/main/utils_subgroup.f90 b/src/main/utils_subgroup.f90 index 913a57606..91f20f713 100644 --- a/src/main/utils_subgroup.f90 +++ b/src/main/utils_subgroup.f90 @@ -10,7 +10,7 @@ module utils_subgroup ! ! :References: None ! -! :Owner: Yrisch +! :Owner: Yann Bernard ! ! :Runtime parameters: None ! diff --git a/src/main/utils_tables.f90 b/src/main/utils_tables.f90 index 47320b69c..747b7663e 100644 --- a/src/main/utils_tables.f90 +++ b/src/main/utils_tables.f90 @@ -19,7 +19,7 @@ module table_utils implicit none public :: yinterp, linspace, logspace, diff, flip_array, interpolator - public :: find_nearest_index, interp_1d + public :: find_nearest_index, interp_1d, linear_interpolator_one_d private @@ -201,11 +201,25 @@ end subroutine find_nearest_index ! 1D linear interpolation routine !+ !----------------------------------------------------------------------- -real function interp_1d(x,x1,x2,y1,y2) +pure real function interp_1d(x,x1,x2,y1,y2) real, intent(in) :: x, x1, x2, y1, y2 interp_1d = y1 + (x-x1)*(y2-y1)/(x2-x1) end function interp_1d +!----------------------------------------------------------------------- +!+ +! similar but just interpolates between two values +! val0 and val1 where u is the fraction of the way +!+ +!----------------------------------------------------------------------- +pure subroutine linear_interpolator_one_d(val0,val1,u,val) + real, intent(out) :: val + real, intent(in) :: val0,val1,u + + val=(1.-u)*val0+u*val1 + +end subroutine linear_interpolator_one_d + end module table_utils diff --git a/src/main/viscosity.f90 b/src/main/viscosity.f90 index 114165a0e..c5d17fbba 100644 --- a/src/main/viscosity.f90 +++ b/src/main/viscosity.f90 @@ -37,7 +37,6 @@ subroutine set_defaults_viscosity shearparam = 0.1 ! alphadisc (if irealvisc=2) or nu if irealvisc=1 bulkvisc = 0.0 ! bulk viscosity parameter in code units - return end subroutine set_defaults_viscosity !---------------------------------------------------------------- diff --git a/src/setup/set_orbit.f90 b/src/setup/set_orbit.f90 index 38f3348b3..a44a77ff1 100644 --- a/src/setup/set_orbit.f90 +++ b/src/setup/set_orbit.f90 @@ -13,12 +13,12 @@ module setorbit ! 0) Campbell elements for bound or unbound orbit (aeiOwf) ! 1) Flyby parameters (periapsis, initial separation, argument of periapsis, inclination) ! 2) position and velocity for both bodies - -! While Campbell elements can be used for unbound orbits, they require -! specifying the true anomaly at the start of the simulation. This is -! not always easy to determine, so the flyby option is provided as an -! alternative. There one specifies the initial separation instead, however -! the choice of angles is more restricted +! +! While Campbell elements can be used for unbound orbits, they require +! specifying the true anomaly at the start of the simulation. This is +! not always easy to determine, so the flyby option is provided as an +! alternative. There one specifies the initial separation instead, however +! the choice of angles is more restricted ! ! :References: None ! diff --git a/src/setup/setup_asteroidwind.f90 b/src/setup/setup_asteroidwind.f90 index 8ccc75fb8..27ce8ce81 100644 --- a/src/setup/setup_asteroidwind.f90 +++ b/src/setup/setup_asteroidwind.f90 @@ -22,7 +22,7 @@ module setup ! - ipot : *wd modelled by 0=sink or 1=externalforce* ! - m1 : *mass of white dwarf (solar mass)* ! - m2 : *mass of asteroid (ceres mass)* -! - mdot : *mass injection rate (g/s)* +! - mdot : *mass injection rate with unit, e.g. 1e8*g/s, 1e-7M_s/yr (from setup)* ! - norbits : *number of orbits* ! - npart_at_end : *number of particles injected after norbits* ! - rasteroid : *radius of asteroid (km)* @@ -33,6 +33,7 @@ module setup ! timestep, units ! use inject, only:mdot + use inject, only:mdot_str implicit none public :: setpart @@ -47,7 +48,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, use part, only:nptmass,xyzmh_ptmass,vxyz_ptmass,ihacc,ihsoft,idust,set_particle_type,igas use setbinary, only:set_binary,get_a_from_period use spherical, only:set_sphere - use units, only:set_units,umass,udist,utime,unit_velocity + use units, only:set_units,umass,udist,unit_velocity,in_code_units use physcon, only:solarm,au,pi,solarr,ceresm,km,kboltz,mass_proton_cgs use externalforces, only:iext_binary, iext_einsteinprec, update_externalforce, & mass1,accradius1 @@ -84,7 +85,8 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, rasteroid = 2338.3 ! (km) gastemp = 5000. ! (K) norbits = 1000. - mdot = 5.e8 ! Mass injection rate (g/s) + mdot = 5.e8 ! Mass injection rate (will change later by the mdot_str) + mdot_str = "5.e8*g/s" ! Mass injection rate with unit, e.g. 1e8*g/s, 1e-7M_s/yr npart_at_end = 1.0e6 ! Number of particles after norbits dumpsperorbit = 1 @@ -180,7 +182,8 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, endif ! we use the estimated injection rate and the final time to set the particle mass - massoftype(igas) = tmax*mdot/(umass/utime)/npart_at_end + mdot = in_code_units(mdot_str,ierr) + massoftype(igas) = tmax*mdot/npart_at_end hfact = hfact_default !npart_old = npart !call inject_particles(time,0.,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,npart,npart_old,npartoftype,dtinj) @@ -219,7 +222,7 @@ subroutine write_setupfile(filename) call write_inopt(norbits, 'norbits', 'number of orbits', iunit) call write_inopt(dumpsperorbit,'dumpsperorbit','number of dumps per orbit', iunit) call write_inopt(npart_at_end,'npart_at_end','number of particles injected after norbits',iunit) - call write_inopt(mdot,'mdot','mass injection rate (g/s)',iunit) + call write_inopt(mdot_str, 'mdot', 'mass injection rate with unit, e.g. 1e8*g/s, 1e-7M_s/yr (from setup)',iunit) close(iunit) end subroutine write_setupfile @@ -227,6 +230,7 @@ end subroutine write_setupfile subroutine read_setupfile(filename,ierr) use infile_utils, only:open_db_from_file,inopts,read_inopt,close_db use io, only:error + use units, only:in_code_units character(len=*), intent(in) :: filename integer, intent(out) :: ierr integer, parameter :: iunit = 21 @@ -248,7 +252,7 @@ subroutine read_setupfile(filename,ierr) call read_inopt(norbits, 'norbits', db,min=0.,errcount=nerr) call read_inopt(dumpsperorbit,'dumpsperorbit',db,min=0 ,errcount=nerr) call read_inopt(npart_at_end, 'npart_at_end', db,min=0 ,errcount=nerr) - call read_inopt(mdot, 'mdot', db,min=0.,errcount=nerr) + call read_inopt(mdot_str, 'mdot', db,errcount=nerr) call close_db(db) if (nerr > 0) then print "(1x,i2,a)",nerr,' error(s) during read of setup file: re-writing...' diff --git a/src/setup/setup_cluster.f90 b/src/setup/setup_cluster.f90 index c1d1ae01a..157414cfa 100644 --- a/src/setup/setup_cluster.f90 +++ b/src/setup/setup_cluster.f90 @@ -22,6 +22,7 @@ module setup ! - mass_fac : *mass unit in Msun* ! - mu : *mean molecular mass* ! - n_particles : *number of particles in sphere* +! - relax : *relax the cloud ?* ! ! :Dependencies: HIIRegion, centreofmass, cooling, datafiles, dim, eos, ! infile_utils, io, kernel, mpidomain, options, part, physcon, prompting, diff --git a/src/setup/setup_solarsystem.f90 b/src/setup/setup_solarsystem.f90 index 5b06d37af..4dcb7a861 100644 --- a/src/setup/setup_solarsystem.f90 +++ b/src/setup/setup_solarsystem.f90 @@ -55,7 +55,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, character(len=20), intent(in) :: fileprefix real, intent(out) :: vxyzu(:,:) character(len=120) :: filename - integer :: ierr,nbodies,i + integer :: ierr,nbodies,i,j,n,nsample logical :: iexist real :: period,semia,mtot,hpart integer, parameter :: max_bodies = 2000000 @@ -100,7 +100,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, semia = 1. ! Earth mtot = solarm/umass - hpart = 100.*au/udist + hpart = 10.*au/udist period = 2.*pi*sqrt(semia**3/mtot) tmax = norbits*period @@ -110,23 +110,36 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, call read_mpc(filename,nbodies,dat=dat) print "(a,i0,a)",' read orbital data for ',nbodies,' minor planets' + n = 0 + nsample = 1 ! can place many particles evenly sampling the orbit if desired do i=1,nbodies ! ! for each solar system object get the xyz positions from the orbital parameters ! !print*,i,'aeiOwM=',dat(i)%a,dat(i)%ecc,dat(i)%inc,dat(i)%O,dat(i)%w,dat(i)%M - call set_binary(mtot,epsilon(0.),dat(i)%a,dat(i)%ecc,0.02,1.e-15,& + do j=1,nsample + n = n + 1 + if (nsample==1) then + call set_binary(mtot,epsilon(0.),dat(i)%a,dat(i)%ecc,0.02,1.e-15,& xyzmh_ptmass,vxyz_ptmass,nptmass,ierr,incl=dat(i)%inc,& arg_peri=dat(i)%w,posang_ascnode=dat(i)%O,& mean_anomaly=dat(i)%M,verbose=.false.) - ! - ! now delete the point masses but set a dust particle as the secondary - ! - nptmass = 0 - xyzh(1:3,i) = xyzmh_ptmass(1:3,2) - xyzh(4,i) = hpart ! give a random length scale as the smoothing length - vxyzu(1:3,i) = vxyz_ptmass(1:3,2) - call set_particle_type(i,idust) + else + call set_binary(mtot,epsilon(0.),dat(i)%a,dat(i)%ecc,0.02,1.e-15,& + xyzmh_ptmass,vxyz_ptmass,nptmass,ierr,incl=dat(i)%inc,& + arg_peri=dat(i)%w,posang_ascnode=dat(i)%O,& + f=360.*(n-1)/nsample,verbose=.false.) + endif + ! + ! now delete the point masses but set a dust particle as the secondary + ! + nptmass = 0 + xyzh(1:3,n) = xyzmh_ptmass(1:3,2) + xyzh(4,n) = hpart ! give a random length scale as the smoothing length + vxyzu(1:3,n) = vxyz_ptmass(1:3,2) + call set_particle_type(n,idust) + enddo + enddo ! ! restore the Sun @@ -135,10 +148,11 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, ! ! set mass of all the minor bodies equal ! - npart = nbodies + npart = nbodies*nsample + print*,' n = ',n,' npart = ',npart ndustlarge = 1 ndusttypes = 1 - npartoftype(idust) = nbodies + npartoftype(idust) = nbodies*nsample massoftype(idust) = 1.e-20 grainsize(1:ndustlarge) = km/udist ! assume km-sized bodies graindens(1:ndustlarge) = 2./unit_density ! 2 g/cm^3 diff --git a/src/setup/setup_starcluster.f90 b/src/setup/setup_starcluster.f90 index 429558843..4a11efd07 100644 --- a/src/setup/setup_starcluster.f90 +++ b/src/setup/setup_starcluster.f90 @@ -11,7 +11,7 @@ module setup ! ! :References: Paumard et al. (2006) ! -! :Owner: Yrisch +! :Owner: Yann Bernard ! ! :Runtime parameters: ! - datafile : *filename for star data (m,x,y,z,vx,vy,vz)* diff --git a/src/utils/einsteintk_utils.f90 b/src/utils/einsteintk_utils.f90 index 6ec6668ef..3268df976 100644 --- a/src/utils/einsteintk_utils.f90 +++ b/src/utils/einsteintk_utils.f90 @@ -14,24 +14,20 @@ module einsteintk_utils ! ! :Runtime parameters: None ! -! :Dependencies: part +! :Dependencies: metric_et_utils, part ! + use metric_et_utils, only:gridorigin,dxgrid,gridsize 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. + real :: boundsize(3) logical :: exact_rendering character(len=128) :: logfilestor,evfilestor,dumpfilestor,infilestor contains subroutine init_etgrid(nx,ny,nz,dx,dy,dz,originx,originy,originz) + use metric_et_utils, only:allocate_grid,gridsize,dxgrid,gridorigin integer, intent(in) :: nx,ny,nz real, intent(in) :: dx,dy,dz,originx,originy,originz @@ -47,39 +43,23 @@ subroutine init_etgrid(nx,ny,nz,dx,dy,dz,originx,originy,originz) 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)) + call allocate_grid(nx,ny,nz,dx,dy,dz,originx,originy,originz) ! 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 + use metric_et_utils, only:print_metric_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) + call print_metric_grid() end subroutine print_etgrid diff --git a/src/utils/struct_part.f90 b/src/utils/struct_part.f90 index 99640148d..781a3c2fd 100644 --- a/src/utils/struct_part.f90 +++ b/src/utils/struct_part.f90 @@ -149,10 +149,10 @@ subroutine get_structure_fn(sf,nbins,norder,distmin,distmax,xbins,ncount,npart,x !$omp reduction(+:sf) do ipt=1,npts !$ if (.false.) then - if (mod(ipt,100)==0) then - call cpu_time(tcpu2) - print*,' ipt = ',ipt,tcpu2-tcpu1 - endif + if (mod(ipt,100)==0) then + call cpu_time(tcpu2) + print*,' ipt = ',ipt,tcpu2-tcpu1 + endif !$ endif i = list(ipt) xpt(1) = xyz(1,i) diff --git a/src/utils/tabulate_metric.f90 b/src/utils/tabulate_metric.f90 new file mode 100644 index 000000000..c02d2a870 --- /dev/null +++ b/src/utils/tabulate_metric.f90 @@ -0,0 +1,84 @@ +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2024 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.github.io/ ! +!--------------------------------------------------------------------------! +program tabulate_metric +! +! tabulate_metric +! +! :References: None +! +! :Owner: Daniel Price +! +! :Usage: tabulate_metric [no arguments] +! +! :Dependencies: metric, metric_et_utils +! + use metric_et_utils + !use metric + + implicit none + + integer :: ierr + character(len=64) :: metric_file = 'tabuled_metric.dat' + + + ! Init grid and tabulated metric + call initialize_grid() + + ! Fill and interpolate metric in the grid + call fill_grid() + + ! Write Data in file + call write_tabulated_metric(metric_file, ierr) + + if (ierr /= 0) then + print *, 'Error writing metric data to file' + else + print *, 'Metric data successfully written to file' + endif + +contains + +subroutine fill_grid() + use metric + integer :: i, j, k + real :: dx, dy, dz + real :: position(3) + real :: gcov(0:3,0:3) + real :: gcon(0:3,0:3) + real :: sqrtg + real :: dgcovdx(0:3,0:3) + real :: dgcovdy(0:3,0:3) + real :: dgcovdz(0:3,0:3) + ! Triple loop to fill the grid + dx = (xmax - xmin) / (nx - 1) + dy = (ymax - ymin) / (ny - 1) + dz = (zmax - zmin) / (nz - 1) + + do i = 1, nx + do j = 1, ny + do k = 1, nz + ! Calculate the current position in the grid + position(1) = xmin + (i - 1) * dx + position(2) = ymin + (j - 1) * dy + position(3) = zmin + (k - 1) * dz + ! Store the calculated values in the grid arrays + call get_metric_cartesian(position,gcov,gcon,sqrtg) + !call get_metric_derivs(position,dgcovdx, dgcovdy, dgcovdz) + call metric_cartesian_derivatives(position,dgcovdx, dgcovdy, dgcovdz) + gcovgrid(:,:,i,j,k) = gcov + gcongrid(:,:,i,j,k) = gcon + sqrtggrid(i,j,k) = sqrtg + metricderivsgrid(:,:,1,i,j,k) = dgcovdx + metricderivsgrid(:,:,2,i,j,k) = dgcovdy + metricderivsgrid(:,:,3,i,j,k) = dgcovdz + enddo + enddo + enddo +end subroutine fill_grid + +end program tabulate_metric + diff --git a/src/utils/utils_ephemeris.f90 b/src/utils/utils_ephemeris.f90 index c6d0a689c..91ad77224 100644 --- a/src/utils/utils_ephemeris.f90 +++ b/src/utils/utils_ephemeris.f90 @@ -71,7 +71,7 @@ subroutine construct_horizons_api_url(object,url,ierr) integer, intent(out) :: ierr character(len=3) :: cmd character(len=10) :: start_epoch,end_epoch - integer(kind=8) :: values(8),year,month,day + integer :: values(8),year,month,day ierr = 0 select case(trim(adjustl(object)))