From 236548f6239f5f9cd8e155a77e88c615bc26958e Mon Sep 17 00:00:00 2001 From: Madeline Nicole Overton Date: Tue, 9 Jul 2024 11:34:27 -0700 Subject: [PATCH 01/31] add inject_randomwind.f90 name change asteroidwind --- src/main/inject_randomwind.f90 | 222 +++++++++++++++++++++++++++++++++ 1 file changed, 222 insertions(+) create mode 100644 src/main/inject_randomwind.f90 diff --git a/src/main/inject_randomwind.f90 b/src/main/inject_randomwind.f90 new file mode 100644 index 000000000..8bf8e9355 --- /dev/null +++ b/src/main/inject_randomwind.f90 @@ -0,0 +1,222 @@ +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2022 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.bitbucket.io/ ! +!--------------------------------------------------------------------------! +module inject +! +! None +! +! :References: None +! +! :Owner: David Liptai +! +! :Runtime parameters: +! - mdot : *mass injection rate in grams/second* +! - mdot_type : *injection rate (0=const, 1=cos(t), 2=r^(-2))* +! - npartperorbit : *particle injection rate in particles/binary orbit* +! - vlag : *percentage lag in velocity of wind* +! +! :Dependencies: binaryutils, externalforces, infile_utils, io, options, +! part, partinject, physcon, random, units +! + use io, only:error + use physcon, only:pi + implicit none + character(len=*), parameter, public :: inject_type = 'randomwind' + real, public :: mdot = 5.e8 ! mass injection rate in grams/second + real,save :: dndt_scaling ! scaling to get ninject correct + + public :: init_inject,inject_particles,write_options_inject,read_options_inject + + private + + 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)) + logical,save :: scaling_set ! has the scaling been set (initially false) + +contains +!----------------------------------------------------------------------- +!+ +! Initialize global variables or arrays needed for injection routine +!+ +!----------------------------------------------------------------------- +subroutine init_inject(ierr) + integer, intent(inout) :: ierr + + scaling_set = .false. + + ierr = 0 + +end subroutine init_inject + +!----------------------------------------------------------------------- +!+ +! Inject particles +!+ +!----------------------------------------------------------------------- +subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& + npart,npartoftype,dtinject) + use io, only:fatal + 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 options, only:iexternalforce + use externalforces,only:mass1 + use binaryutils, only:get_orbit_bits + real, intent(in) :: time, dtlast + real, intent(inout) :: xyzh(:,:), vxyzu(:,:), xyzmh_ptmass(:,:), vxyz_ptmass(:,:) + integer, intent(inout) :: npart + integer, intent(inout) :: npartoftype(:) + real, intent(out) :: dtinject + real, dimension(3) :: xyz,vxyz,r1,r2,v2,vhat,v1 + integer :: i,ipart,npinject,seed,pt + real :: dmdt,robject,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_randomwind','not enough point masses for random wind injection') + if (nptmass > 2) call fatal('inject_randomwind','too many point masses for random wind injection') + + if (nptmass == 2) then + pt = 2 + r1 = xyzmh_ptmass(1:3,1) + m1 = xyzmh_ptmass(4,1) + v1 = vxyz_ptmass(1:3,1) + else + pt = 1 + r1 = 0. + m1 = mass1 + v1 = 0. + endif + + r2 = xyzmh_ptmass(1:3,pt) + robject = xyzmh_ptmass(ihsoft,pt) + m2 = xyzmh_ptmass(4,pt) + v2 = vxyz_ptmass(1:3,pt) + + speed = sqrt(dot_product(v2,v2)) + vhat = v2/speed + + r = sqrt(dot_product(r1-r2,r1-r2)) + +! +! 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 + +!-- How many particles do we need to inject? +! (Seems to need at least eight gas particles to not crash) <-- This statement may or may not be true... +! + if (npartoftype(igas)<8) then + npinject = 8-npartoftype(igas) + else + ! Calculate how many extra particles from previous step to now + dt = time - t_old + inject_this_step = dt*mdot/massoftype(igas)/(umass/utime) + + npinject = max(0, int(0.5 + have_injected + inject_this_step - npartoftype(igas) )) + + ! Save for next step (faster than integrating the whole thing each time) + t_old = time + have_injected = have_injected + inject_this_step + endif + +!-- Randomly inject particles around the asteroids outer 'radius' +! + do i=1,npinject + xyz = r2 + robject*get_random_pos_on_sphere(seed) + vxyz = (1.-vlag/100)*speed*vhat + u = 0. ! setup is isothermal so utherm is not stored + h = hfact*(robject/2.) + ipart = npart + 1 + call add_or_update_particle(igas,xyz,vxyz,h,u,ipart,npart,npartoftype,xyzh,vxyzu) + enddo + + ! + !-- no constraint on timestep + ! + dtinject = huge(dtinject) + +end subroutine inject_particles + +!----------------------------------------------------------------------- +!+ +! Returns dndt(t) depending on which function is chosen +! Note that time in this function is strictly the fraction +! of the orbit, not absolute time +!+ +!----------------------------------------------------------------------- + +real function mdot_func(r,r_ref) + real, intent(in) :: r,r_ref + + select case (mdot_type) + case (2) + mdot_func = (r_ref/r)**2 + case default + mdot_func = 1.0 + end select + +end function mdot_func + +!----------------------------------------------------------------------- +!+ +! Writes input options to the input file. +!+ +!----------------------------------------------------------------------- +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(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) + +end subroutine write_options_inject + +!----------------------------------------------------------------------- +!+ +! Reads input options from the input file. +!+ +!----------------------------------------------------------------------- +subroutine read_options_inject(name,valstring,imatch,igotall,ierr) + use io, only:fatal + character(len=*), intent(in) :: name,valstring + logical, intent(out) :: imatch,igotall + integer, intent(out) :: ierr + integer, save :: ngot = 0 + character(len=30), parameter :: label = 'read_options_inject' + + imatch = .true. + select case(trim(name)) + case('mdot') + read(valstring,*,iostat=ierr) mdot + ngot = ngot + 1 + if (mdot < 0.) call fatal(label,'mdot < 0 in input options') + case('npartperorbit') + read(valstring,*,iostat=ierr) npartperorbit + ngot = ngot + 1 + if (npartperorbit < 0.) call fatal(label,'npartperorbit < 0 in input options') + case('vlag') + read(valstring,*,iostat=ierr) vlag + ngot = ngot + 1 + case('mdot_type') + read(valstring,*,iostat=ierr) mdot_type + ngot = ngot + 1 + case default + imatch = .false. + end select + + igotall = (ngot >= 1) + +end subroutine read_options_inject + +end module inject From 25f0c6043b631f00644e205dd2bd49f9f1749f0a Mon Sep 17 00:00:00 2001 From: Madeline Nicole Overton Date: Tue, 9 Jul 2024 12:00:04 -0700 Subject: [PATCH 02/31] remove *asteroidwind* from Makefile_setups --- build/Makefile_setups | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/build/Makefile_setups b/build/Makefile_setups index d5f34b9a5..5b875708a 100644 --- a/build/Makefile_setups +++ b/build/Makefile_setups @@ -74,10 +74,10 @@ ifeq ($(SETUP), wddisc) DUST=yes endif -ifeq ($(SETUP), asteroidwind) -# asteroid emitting a wind (Trevascus et al. 2021) +ifeq ($(SETUP), randomwind) +# object 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 From 70025fad989dcf4244b552ec2ab057a4e6486366 Mon Sep 17 00:00:00 2001 From: Madeline Overton Date: Thu, 11 Jul 2024 10:12:41 -0700 Subject: [PATCH 03/31] fixed name change inject_asteroidwind to inject_randomwind and added setup randomwind --- build/Makefile_setups | 16 ++- src/main/inject_asteroidwind.f90 | 235 ------------------------------- src/main/inject_randomwind.f90 | 84 ++++++----- 3 files changed, 62 insertions(+), 273 deletions(-) delete mode 100644 src/main/inject_asteroidwind.f90 diff --git a/build/Makefile_setups b/build/Makefile_setups index 5b875708a..e236c8b44 100644 --- a/build/Makefile_setups +++ b/build/Makefile_setups @@ -74,8 +74,8 @@ ifeq ($(SETUP), wddisc) DUST=yes endif -ifeq ($(SETUP), randomwind) -# object emitting a wind (Trevascus et al. 2021) +ifeq ($(SETUP), asteroidwind) +# asteroid emitting a wind (Trevascus et al. 2021) SETUPFILE=setup_asteroidwind.f90 SRCINJECT=utils_binary.f90 inject_randomwind.f90 IND_TIMESTEPS=yes @@ -165,6 +165,18 @@ ifeq ($(SETUP), disc) IND_TIMESTEPS=yes endif +ifeq ($(SETUP), randomwind) +# one component of binary emitting a wind + DISC_VISCOSITY=yes + SETUPFILE=setup_disc.f90 + ANALYSIS= analysis_disc.f90 + SRCINJECT=utils_binary.f90 inject_randomwind.f90 + IND_TIMESTEPS=yes + CONST_AV=yes + ISOTHERMAL=yes + KNOWN_SETUP=yes +endif + ifeq ($(SETUP), grtde) # tidal disruption event in general relativity SETUPFILE= setup_grtde.f90 diff --git a/src/main/inject_asteroidwind.f90 b/src/main/inject_asteroidwind.f90 deleted file mode 100644 index 37072a760..000000000 --- a/src/main/inject_asteroidwind.f90 +++ /dev/null @@ -1,235 +0,0 @@ -!--------------------------------------------------------------------------! -! 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 inject -! -! Injection module for wind from an orbiting asteroid, as used -! in Trevascus et al. (2021) -! -! :References: -! Trevascus et al. (2021), MNRAS 505, L21-L25 -! -! :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* -! -! :Dependencies: binaryutils, externalforces, infile_utils, io, options, -! part, partinject, physcon, random, units -! - 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 - - public :: init_inject,inject_particles,write_options_inject,read_options_inject,& - set_default_options_inject,update_injected_par - - private - - 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)) - -contains -!----------------------------------------------------------------------- -!+ -! Initialize global variables or arrays needed for injection routine -!+ -!----------------------------------------------------------------------- -subroutine init_inject(ierr) - integer, intent(inout) :: ierr - - ierr = 0 - -end subroutine init_inject - -!----------------------------------------------------------------------- -!+ -! Inject particles -!+ -!----------------------------------------------------------------------- -subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& - npart,npart_old,npartoftype,dtinject) - use io, only:fatal - 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 options, only:iexternalforce - use externalforces,only:mass1 - use binaryutils, only:get_orbit_bits - real, intent(in) :: time, dtlast - real, intent(inout) :: xyzh(:,:), vxyzu(:,:), xyzmh_ptmass(:,:), vxyz_ptmass(:,:) - integer, intent(inout) :: npart, npart_old - integer, intent(inout) :: npartoftype(:) - real, intent(out) :: dtinject - 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 :: 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 > 2) & - call fatal('inject_asteroidwind','too many point masses for asteroid wind injection') - - if (nptmass == 2) then - pt = 2 - r1 = xyzmh_ptmass(1:3,1) - m1 = xyzmh_ptmass(4,1) - v1 = vxyz_ptmass(1:3,1) - else - pt = 1 - r1 = 0. - m1 = mass1 - v1 = 0. - endif - - r2 = xyzmh_ptmass(1:3,pt) - rasteroid = xyzmh_ptmass(ihsoft,pt) - m2 = xyzmh_ptmass(4,pt) - v2 = vxyz_ptmass(1:3,pt) - - speed = sqrt(dot_product(v2,v2)) - vhat = v2/speed - - r = sqrt(dot_product(r1-r2,r1-r2)) - - ! - ! 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 - - ! - !-- How many particles do we need to inject? - ! (Seems to need at least eight gas particles to not crash) <-- This statement may or may not be true... - ! - if (npartoftype(igas) < 8) then - npinject = 8-npartoftype(igas) - else - ! Calculate how many extra particles from previous step to now - dt = time - t_old - inject_this_step = dt*mdot/massoftype(igas)/(umass/utime) - - npinject = max(0, int(0.5 + have_injected + inject_this_step - npartoftype(igas) )) - - ! Save for next step (faster than integrating the whole thing each time) - t_old = time - have_injected = have_injected + inject_this_step - endif - - ! - !-- Randomly inject particles around the asteroids outer 'radius'. - ! Only inject them on the side that is facing the central sink - ! - do i=1,npinject - xyz = r2 + rasteroid*get_random_pos_on_sphere(seed) - vxyz = (1.-vlag/100)*speed*vhat - u = 0. ! setup is isothermal so utherm is not stored - h = hfact*(rasteroid/2.) - ipart = npart + 1 - call add_or_update_particle(igas,xyz,vxyz,h,u,ipart,npart,npartoftype,xyzh,vxyzu) - enddo - - ! - !-- no constraint on timestep - ! - dtinject = huge(dtinject) - -end subroutine inject_particles - -subroutine update_injected_par - ! -- placeholder function - ! -- does not do anything and will never be used -end subroutine - -!----------------------------------------------------------------------- -!+ -! Returns dndt(t) depending on which function is chosen -! Note that time in this function is strictly the fraction -! of the orbit, not absolute time -!+ -!----------------------------------------------------------------------- -real function mdot_func(r,r_ref) - real, intent(in) :: r,r_ref - - select case (mdot_type) - case (2) - mdot_func = (r_ref/r)**2 - case default - mdot_func = 1.0 - end select - -end function mdot_func - -!----------------------------------------------------------------------- -!+ -! Writes input options to the input file. -!+ -!----------------------------------------------------------------------- -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(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) - -end subroutine write_options_inject - -!----------------------------------------------------------------------- -!+ -! Reads input options from the input file. -!+ -!----------------------------------------------------------------------- -subroutine read_options_inject(name,valstring,imatch,igotall,ierr) - use io, only:fatal - character(len=*), intent(in) :: name,valstring - logical, intent(out) :: imatch,igotall - integer, intent(out) :: ierr - integer, save :: ngot = 0 - character(len=30), parameter :: label = 'read_options_inject' - - imatch = .true. - select case(trim(name)) - case('mdot') - read(valstring,*,iostat=ierr) mdot - ngot = ngot + 1 - if (mdot < 0.) call fatal(label,'mdot < 0 in input options') - case('npartperorbit') - read(valstring,*,iostat=ierr) npartperorbit - ngot = ngot + 1 - if (npartperorbit < 0.) call fatal(label,'npartperorbit < 0 in input options') - case('vlag') - read(valstring,*,iostat=ierr) vlag - ngot = ngot + 1 - case('mdot_type') - read(valstring,*,iostat=ierr) mdot_type - ngot = ngot + 1 - case default - imatch = .false. - end select - - igotall = (ngot >= 1) - -end subroutine read_options_inject - -subroutine set_default_options_inject(flag) - integer, optional, intent(in) :: flag - -end subroutine set_default_options_inject - -end module inject diff --git a/src/main/inject_randomwind.f90 b/src/main/inject_randomwind.f90 index 8bf8e9355..fc91f058b 100644 --- a/src/main/inject_randomwind.f90 +++ b/src/main/inject_randomwind.f90 @@ -1,22 +1,23 @@ !--------------------------------------------------------------------------! ! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! -! Copyright (c) 2007-2022 The Authors (see AUTHORS) ! +! Copyright (c) 2007-2024 The Authors (see AUTHORS) ! ! See LICENCE file for usage and distribution conditions ! -! http://phantomsph.bitbucket.io/ ! +! http://phantomsph.github.io/ ! !--------------------------------------------------------------------------! module inject ! -! None +! Injection module for wind from an orbiting body, as used +! in Trevascus et al. (2021) ! -! :References: None +! :References: +! Trevascus et al. (2021), MNRAS 505, L21-L25 ! ! :Owner: David Liptai ! ! :Runtime parameters: -! - mdot : *mass injection rate in grams/second* -! - mdot_type : *injection rate (0=const, 1=cos(t), 2=r^(-2))* -! - npartperorbit : *particle injection rate in particles/binary orbit* -! - vlag : *percentage lag in velocity of wind* +! - 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* ! ! :Dependencies: binaryutils, externalforces, infile_utils, io, options, ! part, partinject, physcon, random, units @@ -25,17 +26,16 @@ module inject use physcon, only:pi implicit none character(len=*), parameter, public :: inject_type = 'randomwind' - real, public :: mdot = 5.e8 ! mass injection rate in grams/second - real,save :: dndt_scaling ! scaling to get ninject correct + real, public :: mdot = 5.e8 ! mass injection rate in grams/second - public :: init_inject,inject_particles,write_options_inject,read_options_inject + public :: init_inject,inject_particles,write_options_inject,read_options_inject,& + set_default_options_inject,update_injected_par private 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)) - logical,save :: scaling_set ! has the scaling been set (initially false) contains !----------------------------------------------------------------------- @@ -46,8 +46,6 @@ module inject subroutine init_inject(ierr) integer, intent(inout) :: ierr - scaling_set = .false. - ierr = 0 end subroutine init_inject @@ -58,7 +56,7 @@ end subroutine init_inject !+ !----------------------------------------------------------------------- subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& - npart,npartoftype,dtinject) + npart,npart_old,npartoftype,dtinject) use io, only:fatal use part, only:nptmass,massoftype,igas,hfact,ihsoft use partinject, only:add_or_update_particle @@ -70,19 +68,21 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& use binaryutils, only:get_orbit_bits real, intent(in) :: time, dtlast real, intent(inout) :: xyzh(:,:), vxyzu(:,:), xyzmh_ptmass(:,:), vxyz_ptmass(:,:) - integer, intent(inout) :: npart + integer, intent(inout) :: npart, npart_old integer, intent(inout) :: npartoftype(:) real, intent(out) :: dtinject real, dimension(3) :: xyz,vxyz,r1,r2,v2,vhat,v1 integer :: i,ipart,npinject,seed,pt - real :: dmdt,robject,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_randomwind','not enough point masses for random wind injection') - if (nptmass > 2) call fatal('inject_randomwind','too many point masses for random wind injection') + if (nptmass < 2 .and. iexternalforce == 0) & + call fatal('inject_randomwind','not enough point masses for random wind injection') + if (nptmass > 2) & + call fatal('inject_randomwind','too many point masses for random wind injection') if (nptmass == 2) then pt = 2 @@ -97,7 +97,7 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& endif r2 = xyzmh_ptmass(1:3,pt) - robject = xyzmh_ptmass(ihsoft,pt) + rbody = xyzmh_ptmass(ihsoft,pt) m2 = xyzmh_ptmass(4,pt) v2 = vxyz_ptmass(1:3,pt) @@ -106,15 +106,16 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& r = sqrt(dot_product(r1-r2,r1-r2)) -! -! Add any dependency on radius to mass injection rate (and convert to code units) -! + ! + ! 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 -!-- How many particles do we need to inject? -! (Seems to need at least eight gas particles to not crash) <-- This statement may or may not be true... -! - if (npartoftype(igas)<8) then + ! + !-- How many particles do we need to inject? + ! (Seems to need at least eight gas particles to not crash) <-- This statement may or may not be true... + ! + if (npartoftype(igas) < 8) then npinject = 8-npartoftype(igas) else ! Calculate how many extra particles from previous step to now @@ -128,13 +129,14 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& have_injected = have_injected + inject_this_step endif -!-- Randomly inject particles around the asteroids outer 'radius' -! + ! + !-- Randomly inject particles around the body's outer 'radius'. + ! do i=1,npinject - xyz = r2 + robject*get_random_pos_on_sphere(seed) + xyz = r2 + rbody*get_random_pos_on_sphere(seed) vxyz = (1.-vlag/100)*speed*vhat u = 0. ! setup is isothermal so utherm is not stored - h = hfact*(robject/2.) + h = hfact*(rbody/2.) ipart = npart + 1 call add_or_update_particle(igas,xyz,vxyz,h,u,ipart,npart,npartoftype,xyzh,vxyzu) enddo @@ -146,6 +148,11 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& end subroutine inject_particles +subroutine update_injected_par + ! -- placeholder function + ! -- does not do anything and will never be used +end subroutine + !----------------------------------------------------------------------- !+ ! Returns dndt(t) depending on which function is chosen @@ -153,7 +160,6 @@ end subroutine inject_particles ! of the orbit, not absolute time !+ !----------------------------------------------------------------------- - real function mdot_func(r,r_ref) real, intent(in) :: r,r_ref @@ -175,10 +181,11 @@ 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(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,'mdot','mass injection rate in grams/second',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) end subroutine write_options_inject @@ -219,4 +226,9 @@ subroutine read_options_inject(name,valstring,imatch,igotall,ierr) end subroutine read_options_inject +subroutine set_default_options_inject(flag) + integer, optional, intent(in) :: flag + +end subroutine set_default_options_inject + end module inject From a3a7a0a745716098a156b140db18425bd0d5c0dc Mon Sep 17 00:00:00 2001 From: Madeline Overton Date: Thu, 11 Jul 2024 10:15:52 -0700 Subject: [PATCH 04/31] fixed email --- .mailmap | 1 + 1 file changed, 1 insertion(+) diff --git a/.mailmap b/.mailmap index 8332207b8..c8f50cfe7 100644 --- a/.mailmap +++ b/.mailmap @@ -115,4 +115,5 @@ Amena Faruqi <42060670+amenafaruqi@users.noreply.gi Amena Faruqi Amena Faruqi Alison Young Alison Young Simone Ceppi Simone Ceppi +Madeline Overton Madeline Nicole Overton Nicolás Cuello Nicolas Cuello From 043922d998df9e0b374378829be61f986fa32022 Mon Sep 17 00:00:00 2001 From: Madeline Overton Date: Fri, 12 Jul 2024 05:07:19 -0700 Subject: [PATCH 05/31] fixed compile time choices for setup randomwind --- build/Makefile_setups | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/build/Makefile_setups b/build/Makefile_setups index e236c8b44..1eead4b84 100644 --- a/build/Makefile_setups +++ b/build/Makefile_setups @@ -171,10 +171,10 @@ ifeq ($(SETUP), randomwind) SETUPFILE=setup_disc.f90 ANALYSIS= analysis_disc.f90 SRCINJECT=utils_binary.f90 inject_randomwind.f90 - IND_TIMESTEPS=yes - CONST_AV=yes ISOTHERMAL=yes KNOWN_SETUP=yes + MULTIRUNFILE= multirun.f90 + IND_TIMESTEPS=yes endif ifeq ($(SETUP), grtde) From f7a98c17fe0924ba7c24d187675e237334ff2d89 Mon Sep 17 00:00:00 2001 From: DavidBamba Date: Fri, 12 Jul 2024 14:20:16 -0230 Subject: [PATCH 06/31] (gr) can use tabulated metric independently of einstein toolkit --- build/Makefile | 20 ++++- src/main/initial.F90 | 11 ++- src/main/metric_et.f90 | 83 +++++++++++------- src/main/metric_et_utils.f90 | 150 +++++++++++++++++++++++++++++++++ src/utils/einsteintk_utils.f90 | 36 ++------ src/utils/tabulate_metric.f90 | 67 +++++++++++++++ 6 files changed, 307 insertions(+), 60 deletions(-) create mode 100644 src/main/metric_et_utils.f90 create mode 100644 src/utils/tabulate_metric.f90 diff --git a/build/Makefile b/build/Makefile index a6a52c554..4a609420c 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/src/main/initial.F90 b/src/main/initial.F90 index 4c9a97aa2..6f81d142f 100644 --- a/src/main/initial.F90 +++ b/src/main/initial.F90 @@ -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) + gridinit = .true. + endif call init_readwrite_dumps() diff --git a/src/main/metric_et.f90 b/src/main/metric_et.f90 index ce133ea83..97dceb66d 100644 --- a/src/main/metric_et.f90 +++ b/src/main/metric_et.f90 @@ -32,33 +32,37 @@ module metric !+ !---------------------------------------------------------------- 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 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 @@ -95,17 +99,26 @@ pure subroutine get_metric_spherical(position,gcov,gcon,sqrtg) 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 + use metric_et_utils, only:gridinit +! use grid, only:read_tabulated_metric + 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 + 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 pure subroutine metric_spherical_derivatives(position,dgcovdr, dgcovdtheta, dgcovdphi) @@ -154,6 +167,7 @@ subroutine write_options_metric(iunit) integer, intent(in) :: iunit 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 +180,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 @@ -181,8 +203,8 @@ end subroutine read_options_metric 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 eos_shen, 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 @@ -291,7 +313,7 @@ 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 + 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 @@ -389,7 +411,7 @@ pure subroutine interpolate_metric_derivs(position,dgcovdx, dgcovdy, dgcovdz) end subroutine interpolate_metric_derivs 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 @@ -414,3 +436,4 @@ 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..2a8a89493 --- /dev/null +++ b/src/main/metric_et_utils.f90 @@ -0,0 +1,150 @@ +module metric_et_utils + 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 + + 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 + + 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 + + subroutine print_metric_grid() + ! 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_metric_grid + + 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 + + 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/utils/einsteintk_utils.f90 b/src/utils/einsteintk_utils.f90 index 6ec6668ef..c3fe29ae7 100644 --- a/src/utils/einsteintk_utils.f90 +++ b/src/utils/einsteintk_utils.f90 @@ -16,22 +16,18 @@ module einsteintk_utils ! ! :Dependencies: 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,40 +43,24 @@ 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 - - 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) - + use metric_et_utils, only:print_metric_grid + + call print_metric_grid() + end subroutine print_etgrid subroutine get_particle_rhs(i,vx,vy,vz,fx,fy,fz,e_rhs) diff --git a/src/utils/tabulate_metric.f90 b/src/utils/tabulate_metric.f90 new file mode 100644 index 000000000..1e8232a6c --- /dev/null +++ b/src/utils/tabulate_metric.f90 @@ -0,0 +1,67 @@ +program tabulate_metric + 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 + end do + end do + end do +end subroutine fill_grid + +end program tabulate_metric + From 5714d7bc8d6179f04fd9cf84350644af4ad2eecd Mon Sep 17 00:00:00 2001 From: Shunquan Huang Date: Wed, 24 Jul 2024 17:26:51 -0700 Subject: [PATCH 07/31] an unpadte to randomwind --- build/Makefile_setups | 12 +++- ...asteroidwind.f90 => inject_randomwind.f90} | 55 +++++++++++++++---- src/main/random.f90 | 29 +++++++++- src/main/units.f90 | 36 +++++++++++- src/setup/setup_asteroidwind.f90 | 17 ++++-- 5 files changed, 129 insertions(+), 20 deletions(-) rename src/main/{inject_asteroidwind.f90 => inject_randomwind.f90} (76%) 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/src/main/inject_asteroidwind.f90 b/src/main/inject_randomwind.f90 similarity index 76% rename from src/main/inject_asteroidwind.f90 rename to src/main/inject_randomwind.f90 index 37072a760..e0f78cf44 100644 --- a/src/main/inject_asteroidwind.f90 +++ b/src/main/inject_randomwind.f90 @@ -15,9 +15,12 @@ module inject ! :Owner: David Liptai ! ! :Runtime parameters: +! - mdot_str : *mdot with unit* ! - 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* +! - random_type : random position on the surface, 0 for random, 1 for gaussian +! - delta_theta : standard deviation for the gaussion distribution (random_type=1) ! ! :Dependencies: binaryutils, externalforces, infile_utils, io, options, ! part, partinject, physcon, random, units @@ -26,7 +29,8 @@ module inject 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=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 +39,9 @@ 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) contains !----------------------------------------------------------------------- @@ -61,8 +67,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:umass, utime, in_code_units use options, only:iexternalforce use externalforces,only:mass1 use binaryutils, only:get_orbit_bits @@ -71,6 +77,7 @@ 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 @@ -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,semia) ! Use semi-major axis as r_ref ! !-- 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) )) @@ -134,7 +142,7 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& ! Only inject them on the side that is facing the central sink ! do i=1,npinject - xyz = r2 + rasteroid*get_random_pos_on_sphere(seed) + xyz = r2 + rasteroid*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.) @@ -173,6 +181,25 @@ 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,13 @@ 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(random_type, 'random_type', 'random position on the surface, 0 for random, 1 for gaussian', iunit) + call write_inopt(delta_theta, 'delta_theta', 'standard deviation for the gaussion distribution (random_type=1)', iunit) end subroutine write_options_inject @@ -206,9 +235,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 +248,12 @@ subroutine read_options_inject(name,valstring,imatch,igotall,ierr) case('mdot_type') read(valstring,*,iostat=ierr) mdot_type 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/random.f90 b/src/main/random.f90 index e77444401..4984c9886 100644 --- a/src/main/random.f90 +++ b/src/main/random.f90 @@ -19,7 +19,7 @@ module random ! implicit none public :: ran2,get_random,rayleigh_deviate - public :: get_random_pos_on_sphere,gauss_random + public :: get_random_pos_on_sphere,gauss_random,get_gaussian_pos_on_sphere real, parameter :: pi = 4.*atan(1.) private @@ -167,4 +167,29 @@ real function gauss_random(iseed) end function gauss_random -end module random +!------------------------------------------------------------------------- +! +! get 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 + end do + 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 + +end module random \ No newline at end of file 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/setup/setup_asteroidwind.f90 b/src/setup/setup_asteroidwind.f90 index 8ccc75fb8..0a2af47e0 100644 --- a/src/setup/setup_asteroidwind.f90 +++ b/src/setup/setup_asteroidwind.f90 @@ -22,7 +22,8 @@ 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 +! - mdot_str : *mdot with unit* ! - norbits : *number of orbits* ! - npart_at_end : *number of particles injected after norbits* ! - rasteroid : *radius of asteroid (km)* @@ -33,6 +34,7 @@ module setup ! timestep, units ! use inject, only:mdot + use inject, only:mdot_str implicit none public :: setpart @@ -47,7 +49,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,utime,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 +86,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 +183,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 +223,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 +231,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 +253,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...' From 0a42f84b8cfa533e37d73ae8d67ea64ee6731905 Mon Sep 17 00:00:00 2001 From: Madeline Overton Date: Tue, 30 Jul 2024 11:01:30 -0700 Subject: [PATCH 08/31] Unnecessary Setup 'randomwind' is removed --- build/Makefile_setups | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/build/Makefile_setups b/build/Makefile_setups index 8578aa771..faee3ef76 100644 --- a/build/Makefile_setups +++ b/build/Makefile_setups @@ -165,18 +165,6 @@ ifeq ($(SETUP), disc) IND_TIMESTEPS=yes endif -ifeq ($(SETUP), randomwind) -# one component of binary emitting a wind - DISC_VISCOSITY=yes - SETUPFILE=setup_disc.f90 - ANALYSIS= analysis_disc.f90 - SRCINJECT=utils_binary.f90 inject_randomwind.f90 - ISOTHERMAL=yes - KNOWN_SETUP=yes - MULTIRUNFILE= multirun.f90 - IND_TIMESTEPS=yes -endif - ifeq ($(SETUP), grtde) # tidal disruption event in general relativity SETUPFILE= setup_grtde.f90 From 6895b7fb9154e42458a65db9271e1f931ef99028 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 09:39:07 +1000 Subject: [PATCH 09/31] (solarsystem) build failure on gfortran v14 fixed; allow sampling of minor body orbits --- src/setup/setup_solarsystem.f90 | 40 ++++++++++++++++++++++----------- src/utils/utils_ephemeris.f90 | 2 +- 2 files changed, 28 insertions(+), 14 deletions(-) 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/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))) From 7df0b70137435610b94ca0be4fd748f26c00c010 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 09:59:21 +1000 Subject: [PATCH 10/31] (metric_et) remove dependency on eos_shen --- src/main/metric_et.f90 | 121 ++++++++++++++++++++------------------ src/main/utils_tables.f90 | 18 +++++- 2 files changed, 81 insertions(+), 58 deletions(-) diff --git a/src/main/metric_et.f90 b/src/main/metric_et.f90 index 97dceb66d..a4a50b22c 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: einsteintk_utils, infile_utils, table_utils ! implicit none character(len=*), parameter :: metric_type = 'et' @@ -28,7 +28,8 @@ 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) @@ -41,8 +42,7 @@ pure subroutine get_metric_cartesian(position,gcov,gcon,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 + ! Perform trilinear interpolation if ( .not. gridinit) then ierr = 1 ! This is required for phantomsetup @@ -68,8 +68,14 @@ pure subroutine get_metric_cartesian(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) @@ -99,28 +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 metric_et_utils, only:gridinit -! use grid, only:read_tabulated_metric - 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 - 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 + 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 + 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 @@ -140,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) @@ -155,6 +177,7 @@ pure subroutine cartesian2spherical(xcart,xspher) phi = atan2(y,x) xspher = (/r,theta,phi/) + end subroutine cartesian2spherical !----------------------------------------------------------------------- @@ -163,7 +186,7 @@ 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' @@ -180,12 +203,12 @@ subroutine read_options_metric(name,valstring,imatch,igotall,ierr) character(len=*), intent(in) :: name,valstring logical, intent(out) :: imatch,igotall integer, intent(out) :: ierr - integer, save :: ngot = 0 + !integer, save :: ngot = 0 select case(trim(name)) - !case('metric_file') - ! read(valstring,*,iostat=ierr) metric_file - ! ngot = ngot + 1 + !case('metric_file') + ! read(valstring,*,iostat=ierr) metric_file + ! ngot = ngot + 1 case default imatch = .false. end select @@ -199,11 +222,10 @@ 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 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) @@ -219,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) @@ -308,11 +320,15 @@ 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 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) @@ -322,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) @@ -405,11 +414,13 @@ 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 metric_et_utils, only:gridorigin real, intent(in) :: position(3) @@ -431,9 +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/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 From 3122cbd3fe8d7f3e627dd23b1bb060bd0a9ef99f Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 10:12:42 +1000 Subject: [PATCH 11/31] (randomwind) fixed warnings --- src/main/inject_randomwind.f90 | 2 +- src/main/random.f90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/main/inject_randomwind.f90 b/src/main/inject_randomwind.f90 index 22dfe316a..3e9d7edc8 100644 --- a/src/main/inject_randomwind.f90 +++ b/src/main/inject_randomwind.f90 @@ -68,7 +68,7 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& use partinject, only:add_or_update_particle use physcon, only:twopi,gg,kboltz,mass_proton_cgs use random, only:get_random_pos_on_sphere, get_gaussian_pos_on_sphere - use units, only:umass, utime, in_code_units + use units, only:in_code_units use options, only:iexternalforce use externalforces,only:mass1 use binaryutils, only:get_orbit_bits diff --git a/src/main/random.f90 b/src/main/random.f90 index 058f29842..e4adeaba0 100644 --- a/src/main/random.f90 +++ b/src/main/random.f90 @@ -20,7 +20,7 @@ module random implicit none public :: ran2,get_random,rayleigh_deviate public :: get_random_pos_on_sphere,get_gaussian_pos_on_sphere - public :: gauss_random,divide_unit_seq + public :: gauss_random,divide_unit_seg real, parameter :: pi = 4.*atan(1.) private From c9e1e4e9f33e4d84751010f82c03b0d7a284003f Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 10:14:53 +1000 Subject: [PATCH 12/31] (randomwind) test failure fixed --- src/main/inject_randomwind.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/inject_randomwind.f90 b/src/main/inject_randomwind.f90 index 3e9d7edc8..96bc565f3 100644 --- a/src/main/inject_randomwind.f90 +++ b/src/main/inject_randomwind.f90 @@ -86,7 +86,7 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& real, save :: have_injected,t_old real, save :: semia - if (nptmass < 2 .and. iexternalforce == 0) & + if (nptmass < 1 .and. iexternalforce == 0) & call fatal('inject_randomwind','not enough point masses for random wind injection') if (nptmass > 2) & call fatal('inject_randomwind','too many point masses for random wind injection') From 61838a3e7ff17649cc3c8698526e93a9ac7d6462 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 10:18:07 +1000 Subject: [PATCH 13/31] [header-bot] updated file headers --- src/main/H2regions.f90 | 3 --- src/main/initial.F90 | 10 +++++----- src/main/inject_randomwind.f90 | 11 +++++------ src/main/metric_et.f90 | 2 +- src/main/metric_et_utils.f90 | 17 +++++++++++++++++ src/main/readwrite_dumps_common.f90 | 2 +- src/main/readwrite_dumps_fortran.f90 | 2 +- src/setup/set_orbit.f90 | 15 +++++++++++++++ src/setup/setup_asteroidwind.f90 | 3 +-- src/setup/setup_cluster.f90 | 1 + src/utils/einsteintk_utils.f90 | 2 +- src/utils/tabulate_metric.f90 | 17 +++++++++++++++++ 12 files changed, 65 insertions(+), 20 deletions(-) diff --git a/src/main/H2regions.f90 b/src/main/H2regions.f90 index 8d81b4a12..5ca0ad90d 100644 --- a/src/main/H2regions.f90 +++ b/src/main/H2regions.f90 @@ -17,9 +17,6 @@ module HIIRegion ! :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 diff --git a/src/main/initial.F90 b/src/main/initial.F90 index 626f4e9a1..c1f0b6842 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 diff --git a/src/main/inject_randomwind.f90 b/src/main/inject_randomwind.f90 index 96bc565f3..f2f8648e8 100644 --- a/src/main/inject_randomwind.f90 +++ b/src/main/inject_randomwind.f90 @@ -15,12 +15,11 @@ module inject ! :Owner: David Liptai ! ! :Runtime parameters: -! - mdot_str : *mdot with unit* -! - 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* -! - random_type : random position on the surface, 0 for random, 1 for gaussian -! - delta_theta : standard deviation for the gaussion distribution (random_type=1) +! - 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 diff --git a/src/main/metric_et.f90 b/src/main/metric_et.f90 index a4a50b22c..a4c757a73 100644 --- a/src/main/metric_et.f90 +++ b/src/main/metric_et.f90 @@ -14,7 +14,7 @@ module metric ! ! :Runtime parameters: None ! -! :Dependencies: einsteintk_utils, infile_utils, table_utils +! :Dependencies: metric_et_utils, table_utils ! implicit none character(len=*), parameter :: metric_type = 'et' diff --git a/src/main/metric_et_utils.f90 b/src/main/metric_et_utils.f90 index 2a8a89493..d48ec5930 100644 --- a/src/main/metric_et_utils.f90 +++ b/src/main/metric_et_utils.f90 @@ -1,4 +1,21 @@ +!--------------------------------------------------------------------------! +! 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 +! +! metric_et_utils +! +! :References: None +! +! :Owner: DavidBamba +! +! :Runtime parameters: None +! +! :Dependencies: None +! implicit none real, allocatable :: gcovgrid(:,:,:,:,:) 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/setup/set_orbit.f90 b/src/setup/set_orbit.f90 index 38f3348b3..9b05809b5 100644 --- a/src/setup/set_orbit.f90 +++ b/src/setup/set_orbit.f90 @@ -14,6 +14,21 @@ module setorbit ! 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 +! +! :References: None +! +! :Owner: Daniel Price +! +! :Runtime parameters: None +! +! :Dependencies: infile_utils, physcon, setbinary, setflyby, units +! + ! 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 diff --git a/src/setup/setup_asteroidwind.f90 b/src/setup/setup_asteroidwind.f90 index 0a2af47e0..289b8329d 100644 --- a/src/setup/setup_asteroidwind.f90 +++ b/src/setup/setup_asteroidwind.f90 @@ -22,8 +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 -! - mdot_str : *mdot with unit* +! - 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)* 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/utils/einsteintk_utils.f90 b/src/utils/einsteintk_utils.f90 index c3fe29ae7..45f56c13c 100644 --- a/src/utils/einsteintk_utils.f90 +++ b/src/utils/einsteintk_utils.f90 @@ -14,7 +14,7 @@ module einsteintk_utils ! ! :Runtime parameters: None ! -! :Dependencies: part +! :Dependencies: metric_et_utils, part ! use metric_et_utils, only:gridorigin,dxgrid,gridsize implicit none diff --git a/src/utils/tabulate_metric.f90 b/src/utils/tabulate_metric.f90 index 1e8232a6c..5d6a37869 100644 --- a/src/utils/tabulate_metric.f90 +++ b/src/utils/tabulate_metric.f90 @@ -1,4 +1,21 @@ +!--------------------------------------------------------------------------! +! 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: DavidBamba +! +! :Usage: tabulate_metric [no arguments] +! +! :Dependencies: metric, metric_et_utils +! use metric_et_utils !use metric From d5a13a042b98faab79a77f2be419e2c1577e9fb0 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 10:18:23 +1000 Subject: [PATCH 14/31] [space-bot] whitespace at end of lines removed --- src/main/metric_et_utils.f90 | 18 +++++++++--------- src/utils/einsteintk_utils.f90 | 4 ++-- src/utils/tabulate_metric.f90 | 6 +++--- 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/main/metric_et_utils.f90 b/src/main/metric_et_utils.f90 index d48ec5930..60eca7854 100644 --- a/src/main/metric_et_utils.f90 +++ b/src/main/metric_et_utils.f90 @@ -32,9 +32,9 @@ module metric_et_utils 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 - + 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 @@ -45,11 +45,11 @@ subroutine allocate_grid(nxin,nyin,nzin,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 @@ -57,14 +57,14 @@ subroutine allocate_grid(nxin,nyin,nzin,dx,dy,dz,originx,originy,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 subroutine initialize_grid() @@ -89,11 +89,11 @@ end subroutine initialize_grid subroutine print_metric_grid() ! 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_metric_grid subroutine write_tabulated_metric(metric_file, ierr) @@ -130,7 +130,7 @@ 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) diff --git a/src/utils/einsteintk_utils.f90 b/src/utils/einsteintk_utils.f90 index 45f56c13c..3268df976 100644 --- a/src/utils/einsteintk_utils.f90 +++ b/src/utils/einsteintk_utils.f90 @@ -58,9 +58,9 @@ end subroutine init_etgrid subroutine print_etgrid() use metric_et_utils, only:print_metric_grid - + call print_metric_grid() - + end subroutine print_etgrid subroutine get_particle_rhs(i,vx,vy,vz,fx,fy,fz,e_rhs) diff --git a/src/utils/tabulate_metric.f90 b/src/utils/tabulate_metric.f90 index 5d6a37869..a2ddf0537 100644 --- a/src/utils/tabulate_metric.f90 +++ b/src/utils/tabulate_metric.f90 @@ -20,14 +20,14 @@ program tabulate_metric !use metric implicit none - - integer :: ierr + + 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() From c294323bd70184eb315a0c75b58c6671b5cf3226 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 10:18:23 +1000 Subject: [PATCH 15/31] [author-bot] updated AUTHORS file --- AUTHORS | 34 +++++++++++++++++++--------------- 1 file changed, 19 insertions(+), 15 deletions(-) diff --git a/AUTHORS b/AUTHORS index 649e5eb20..b4d2b8428 100644 --- a/AUTHORS +++ b/AUTHORS @@ -26,50 +26,54 @@ Terrence Tricco Stephane Michoulier Simone Ceppi Spencer Magnall -Caitlyn Hardiman Enrico Ragusa +Caitlyn Hardiman Cristiano Longarini Sergei Biriukov Giovanni Dipierro Roberto Iaconi -Amena Faruqi Hauke Worpel +Amena Faruqi Alison Young Stephen Neilson <36410751+s-neilson@users.noreply.github.com> Martina Toscani Benedetta Veronesi -Sahl Rowther Simon Glover +Sahl Rowther Thomas Reichardt Jean-François Gonzalez Christopher Russell Alessia Franchini Alex Pettitt Jolien Malfait +Madeline Overton Phantom benchmark bot -Kieran Hirsh Mike Lau +Kieran Hirsh Nicole Rodrigues +Nicolás Cuello David Trevascus Farzana Meru -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 -Hugh Griffiths -Jeremy Smallwood +Joe Fisher +Giulia Ballabio +Benoit Commercon +Zachary Pellow +Madeline Overton <85810161+moverton000@users.noreply.github.com> +DavidBamba Jorge Cuadra -MICHOULIER Stephane +Cox, Samuel Steven Rieder Stéven Toupin Taj Jankovič +Chunliang Mu +MICHOULIER Stephane rebeccagmartin <74937128+rebeccagmartin@users.noreply.github.com> +Ariel Chitan +Hugh Griffiths +Shunquan Huang +Jeremy Smallwood From 0b8ea0931be320c41749892d1a062637735108b5 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 10:18:43 +1000 Subject: [PATCH 16/31] [format-bot] end if -> endif; end do -> enddo; if( -> if ( --- src/main/evolve.F90 | 2 +- src/main/metric_et_utils.f90 | 4 ++-- src/main/random.f90 | 2 +- src/utils/tabulate_metric.f90 | 6 +++--- 4 files changed, 7 insertions(+), 7 deletions(-) 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/metric_et_utils.f90 b/src/main/metric_et_utils.f90 index 60eca7854..05d5b75e5 100644 --- a/src/main/metric_et_utils.f90 +++ b/src/main/metric_et_utils.f90 @@ -102,7 +102,7 @@ subroutine write_tabulated_metric(metric_file, ierr) integer :: iunit ! Open the file for writing - open(newunit=iunit, file=metric_file, status='replace', form='unformatted',action='write', iostat=ierr) + open(newunit=iunit,file=metric_file,status='replace',form='unformatted',action='write',iostat=ierr) if (ierr /= 0) then ierr = 1 return @@ -133,7 +133,7 @@ subroutine read_tabulated_metric(metric_file, ierr) ! Open the file for reading - open(newunit=iunit, file=metric_file, status='old', form='unformatted', action='read', iostat=ierr) + open(newunit=iunit,file=metric_file,status='old', form='unformatted',action='read',iostat=ierr) if (ierr /= 0) return ! Read the dimensions of the grid diff --git a/src/main/random.f90 b/src/main/random.f90 index e4adeaba0..dd2ba97c1 100644 --- a/src/main/random.f90 +++ b/src/main/random.f90 @@ -165,7 +165,7 @@ function get_gaussian_pos_on_sphere(iseed, deltheta) result(dx) gauss_theta = gauss_random(iseed) * deltheta do while (abs(gauss_theta) > 1.0) gauss_theta = gauss_random(iseed) * deltheta - end do + enddo theta = acos(gauss_theta) sintheta = sin(theta) costheta = cos(theta) diff --git a/src/utils/tabulate_metric.f90 b/src/utils/tabulate_metric.f90 index a2ddf0537..7a73bcac3 100644 --- a/src/utils/tabulate_metric.f90 +++ b/src/utils/tabulate_metric.f90 @@ -75,9 +75,9 @@ subroutine fill_grid() metricderivsgrid(:,:,1,i,j,k) = dgcovdx metricderivsgrid(:,:,2,i,j,k) = dgcovdy metricderivsgrid(:,:,3,i,j,k) = dgcovdz - end do - end do - end do + enddo + enddo + enddo end subroutine fill_grid end program tabulate_metric From f4ffc6a238ec573eb07924e9b62dae3ea4826af9 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 10:19:04 +1000 Subject: [PATCH 17/31] [indent-bot] standardised indentation --- src/main/metric_et_utils.f90 | 268 +++++++++++++++++----------------- src/utils/struct_part.f90 | 8 +- src/utils/tabulate_metric.f90 | 100 ++++++------- 3 files changed, 188 insertions(+), 188 deletions(-) diff --git a/src/main/metric_et_utils.f90 b/src/main/metric_et_utils.f90 index 05d5b75e5..05b106208 100644 --- a/src/main/metric_et_utils.f90 +++ b/src/main/metric_et_utils.f90 @@ -16,152 +16,152 @@ module metric_et_utils ! ! :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 - - 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 + 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 + +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 - - 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 - - subroutine print_metric_grid() - ! 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_metric_grid - - 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 + 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 + +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 + +subroutine print_metric_grid() + ! 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_metric_grid + +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 - subroutine read_tabulated_metric(metric_file, ierr) - character(len=*), intent(in) :: metric_file - integer, intent(out) :: ierr - integer :: iunit +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 + ! 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 dimensions of the grid + read(iunit) gridsize - ! Read the grid origin and spacing - read(iunit) gridorigin - read(iunit) dxgrid + ! Read the grid origin and spacing + read(iunit) gridorigin + read(iunit) dxgrid - nx = gridsize(1) - ny = gridsize(2) - nz = gridsize(3) + nx = gridsize(1) + ny = gridsize(2) + nz = gridsize(3) - call allocate_grid(nx,ny,nz,& + 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 + ! Read the metric values from the file + read(iunit) gcovgrid + read(iunit) gcongrid + read(iunit) sqrtggrid + read(iunit) metricderivsgrid - gridinit = .true. + gridinit = .true. - ! Close the file - close(iunit) - ierr = 0 - end subroutine read_tabulated_metric + ! Close the file + close(iunit) + ierr = 0 +end subroutine read_tabulated_metric end module metric_et_utils 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 index 7a73bcac3..9216df362 100644 --- a/src/utils/tabulate_metric.f90 +++ b/src/utils/tabulate_metric.f90 @@ -16,68 +16,68 @@ program tabulate_metric ! ! :Dependencies: metric, metric_et_utils ! - use metric_et_utils - !use metric + use metric_et_utils + !use metric - implicit none + implicit none - integer :: ierr - character(len=64) :: metric_file = 'tabuled_metric.dat' + integer :: ierr + character(len=64) :: metric_file = 'tabuled_metric.dat' - ! Init grid and tabulated metric - call initialize_grid() + ! Init grid and tabulated metric + call initialize_grid() - ! Fill and interpolate metric in the grid - call fill_grid() + ! Fill and interpolate metric in the grid + call fill_grid() - ! Write Data in file - call write_tabulated_metric(metric_file, ierr) + ! 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 + 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) + 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 + 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 From bf0ecdbb965c628990358c8de3c7aa34903d6281 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 10:35:10 +1000 Subject: [PATCH 18/31] update .mailmap --- .mailmap | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/.mailmap b/.mailmap index c8f50cfe7..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 @@ -116,4 +118,10 @@ 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 From 82fe89a338abbd971a28aa2e227d3829d6725c51 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 10:38:46 +1000 Subject: [PATCH 19/31] [header-bot] updated file headers --- src/main/H2regions.f90 | 2 +- src/main/eos_HIIR.f90 | 2 +- src/main/metric_et_utils.f90 | 2 +- src/main/subgroup.f90 | 2 +- src/main/substepping.F90 | 2 +- src/main/utils_kepler.f90 | 2 +- src/main/utils_subgroup.f90 | 2 +- src/setup/set_orbit.f90 | 15 +++++++++++++++ src/setup/setup_starcluster.f90 | 2 +- src/utils/tabulate_metric.f90 | 2 +- 10 files changed, 24 insertions(+), 9 deletions(-) diff --git a/src/main/H2regions.f90 b/src/main/H2regions.f90 index 5ca0ad90d..507d42a7c 100644 --- a/src/main/H2regions.f90 +++ b/src/main/H2regions.f90 @@ -10,7 +10,7 @@ module HIIRegion ! ! :References: Fujii et al. (2021), Hopkins et al. (2012) ! -! :Owner: Yrisch +! :Owner: Yann Bernard ! ! :Runtime parameters: None ! 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/metric_et_utils.f90 b/src/main/metric_et_utils.f90 index 05b106208..288520cf5 100644 --- a/src/main/metric_et_utils.f90 +++ b/src/main/metric_et_utils.f90 @@ -10,7 +10,7 @@ module metric_et_utils ! ! :References: None ! -! :Owner: DavidBamba +! :Owner: Daniel Price ! ! :Runtime parameters: None ! 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/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/setup/set_orbit.f90 b/src/setup/set_orbit.f90 index 9b05809b5..f4ebbe1c8 100644 --- a/src/setup/set_orbit.f90 +++ b/src/setup/set_orbit.f90 @@ -29,6 +29,21 @@ module setorbit ! :Dependencies: infile_utils, physcon, setbinary, setflyby, units ! +! 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 +! +! :Owner: Daniel Price +! +! :Runtime parameters: None +! +! :Dependencies: infile_utils, physcon, setbinary, setflyby, units +! + ! 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 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/tabulate_metric.f90 b/src/utils/tabulate_metric.f90 index 9216df362..c02d2a870 100644 --- a/src/utils/tabulate_metric.f90 +++ b/src/utils/tabulate_metric.f90 @@ -10,7 +10,7 @@ program tabulate_metric ! ! :References: None ! -! :Owner: DavidBamba +! :Owner: Daniel Price ! ! :Usage: tabulate_metric [no arguments] ! From 2f7c230b93e1fa6b0bad871b507082893a52edd4 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 10:39:04 +1000 Subject: [PATCH 20/31] [author-bot] updated AUTHORS file --- AUTHORS | 40 ++++++++++++++++++---------------------- 1 file changed, 18 insertions(+), 22 deletions(-) diff --git a/AUTHORS b/AUTHORS index b4d2b8428..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 @@ -32,10 +32,10 @@ Cristiano Longarini Sergei Biriukov Giovanni Dipierro Roberto Iaconi -Hauke Worpel Amena Faruqi +Hauke Worpel Alison Young -Stephen Neilson <36410751+s-neilson@users.noreply.github.com> +Stephen Nielson Martina Toscani Benedetta Veronesi Simon Glover @@ -43,37 +43,33 @@ Sahl Rowther Thomas Reichardt Jean-François Gonzalez Christopher Russell -Alessia Franchini -Alex Pettitt -Jolien Malfait Madeline Overton +Alex Pettitt Phantom benchmark bot -Mike Lau -Kieran Hirsh +Jolien Malfait +Alessia Franchini Nicole Rodrigues -Nicolás Cuello -David Trevascus +Kieran Hirsh Farzana Meru +David Trevascus +Nicolás Cuello Chris Nixon Miguel Gonzalez-Bolivar Maxime Lombart +Zachary Pellow Orsola De Marco -s-neilson <36410751+s-neilson@users.noreply.github.com> Joe Fisher -Giulia Ballabio Benoit Commercon -Zachary Pellow -Madeline Overton <85810161+moverton000@users.noreply.github.com> -DavidBamba +Giulia Ballabio +Rebecca Martin Jorge Cuadra +Hugh Griffiths +Jeremy Smallwood +David Bamba Cox, Samuel +Chunliang Mu +Shunquan Huang Steven Rieder Stéven Toupin Taj Jankovič -Chunliang Mu -MICHOULIER Stephane -rebeccagmartin <74937128+rebeccagmartin@users.noreply.github.com> Ariel Chitan -Hugh Griffiths -Shunquan Huang -Jeremy Smallwood From c1ee6a9ffa34aeb6777afca482f5d0ba50d0ef92 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 17:47:02 +1000 Subject: [PATCH 21/31] (#55) remove ifdefs, fix test failure with SETUP=flrw and properly comment subroutines --- src/main/initial.F90 | 2 +- src/main/metric_et.f90 | 2 +- src/main/metric_et_utils.f90 | 36 +++++++++++-- src/main/metric_tools.F90 | 87 ++++++++++++------------------ src/main/utils_gr.F90 | 91 ++++++++++++++++---------------- src/setup/setup_asteroidwind.f90 | 2 +- 6 files changed, 115 insertions(+), 105 deletions(-) diff --git a/src/main/initial.F90 b/src/main/initial.F90 index c1f0b6842..d050a7f8c 100644 --- a/src/main/initial.F90 +++ b/src/main/initial.F90 @@ -106,7 +106,7 @@ subroutine initialise() ! if (gr .and. metric_type=='et') then call read_tabulated_metric('tabuled_metric.dat',ierr) - gridinit = .true. + if (ierr == 0) gridinit = .true. endif call init_readwrite_dumps() diff --git a/src/main/metric_et.f90 b/src/main/metric_et.f90 index a4c757a73..cd93c9c3b 100644 --- a/src/main/metric_et.f90 +++ b/src/main/metric_et.f90 @@ -189,7 +189,7 @@ 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' + !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 diff --git a/src/main/metric_et_utils.f90 b/src/main/metric_et_utils.f90 index 288520cf5..a3c3bebf5 100644 --- a/src/main/metric_et_utils.f90 +++ b/src/main/metric_et_utils.f90 @@ -6,9 +6,9 @@ !--------------------------------------------------------------------------! module metric_et_utils ! -! metric_et_utils +! Utilities for handling tabulated metrics from the Einstein Toolkit ! -! :References: None +! :References: Magnall et al. (2023), Phys. Rev. D 108, 103534 ! ! :Owner: Daniel Price ! @@ -35,6 +35,11 @@ module metric_et_utils 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 @@ -67,6 +72,12 @@ subroutine allocate_grid(nxin,nyin,nzin,dx,dy,dz,originx,originy,originz) 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) @@ -87,8 +98,12 @@ subroutine initialize_grid() end subroutine initialize_grid +!--------------------------------------------------------------- +!+ +! print information about the metric grid +!+ +!--------------------------------------------------------------- subroutine print_metric_grid() - ! Subroutine for printing quantities of the ET grid print*, "Grid spacing (x,y,z) is : ", dxgrid print*, "Grid origin (x,y,z) is: ", gridorigin @@ -96,6 +111,11 @@ subroutine print_metric_grid() 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 @@ -124,16 +144,21 @@ subroutine write_tabulated_metric(metric_file, ierr) ! 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) + open(newunit=iunit,file=metric_file,status='old',form='unformatted',action='read',iostat=ierr) if (ierr /= 0) return ! Read the dimensions of the grid @@ -162,6 +187,7 @@ subroutine read_tabulated_metric(metric_file, ierr) ! 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 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/utils_gr.F90 b/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/setup/setup_asteroidwind.f90 b/src/setup/setup_asteroidwind.f90 index 289b8329d..27ce8ce81 100644 --- a/src/setup/setup_asteroidwind.f90 +++ b/src/setup/setup_asteroidwind.f90 @@ -48,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,in_code_units + 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 From 728aafad1ceadfc7a9b2ad66cc7064a81bc002db Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 17:49:17 +1000 Subject: [PATCH 22/31] (#55) .F90->.f90 for metric utility routines --- src/main/{metric_tools.F90 => metric_tools.f90} | 0 src/main/{utils_gr.F90 => utils_gr.f90} | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename src/main/{metric_tools.F90 => metric_tools.f90} (100%) rename src/main/{utils_gr.F90 => utils_gr.f90} (100%) diff --git a/src/main/metric_tools.F90 b/src/main/metric_tools.f90 similarity index 100% rename from src/main/metric_tools.F90 rename to src/main/metric_tools.f90 diff --git a/src/main/utils_gr.F90 b/src/main/utils_gr.f90 similarity index 100% rename from src/main/utils_gr.F90 rename to src/main/utils_gr.f90 From 6d780da71bb3b4ee55fbaa008f780a1a401deeaf Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 17:52:11 +1000 Subject: [PATCH 23/31] (set_orbit) fix problem with header bot giving repeated splats of same paragraph in header --- src/setup/set_orbit.f90 | 40 +++++----------------------------------- 1 file changed, 5 insertions(+), 35 deletions(-) diff --git a/src/setup/set_orbit.f90 b/src/setup/set_orbit.f90 index f4ebbe1c8..a44a77ff1 100644 --- a/src/setup/set_orbit.f90 +++ b/src/setup/set_orbit.f90 @@ -13,42 +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 -! -! :References: None -! -! :Owner: Daniel Price -! -! :Runtime parameters: None -! -! :Dependencies: infile_utils, physcon, setbinary, setflyby, units -! - -! 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 -! -! :Owner: Daniel Price ! -! :Runtime parameters: None -! -! :Dependencies: infile_utils, physcon, setbinary, setflyby, units -! - -! 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 ! From d0b8020774c801d1cab5e5b51318c7cfd217ad18 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 18:55:31 +1000 Subject: [PATCH 24/31] test failure fixed --- src/main/inject_randomwind.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/main/inject_randomwind.f90 b/src/main/inject_randomwind.f90 index f2f8648e8..d1b7d694c 100644 --- a/src/main/inject_randomwind.f90 +++ b/src/main/inject_randomwind.f90 @@ -41,6 +41,8 @@ module inject 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. contains !----------------------------------------------------------------------- @@ -82,8 +84,6 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& 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 < 1 .and. iexternalforce == 0) & call fatal('inject_randomwind','not enough point masses for random wind injection') @@ -116,7 +116,7 @@ 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) ! mdot = in_code_units(mdot_str,ierr) - dmdt = mdot*mdot_func(r,semia) ! Use semi-major axis as r_ref + dmdt = mdot*mdot_func(r,rbody) ! Use rbody as r_ref, currently the softening length of the body ! !-- How many particles do we need to inject? From ac7156f1cf1d0ebfd71f4fc79bb1171864268a29 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 19:58:15 +1000 Subject: [PATCH 25/31] (HII regions) minor formatting issues fixed; do not print useless lines to .in file --- src/main/H2regions.f90 | 75 ++++++++++++++++++++++++++---------------- 1 file changed, 47 insertions(+), 28 deletions(-) diff --git a/src/main/H2regions.f90 b/src/main/H2regions.f90 index 507d42a7c..cac51a26e 100644 --- a/src/main/H2regions.f90 +++ b/src/main/H2regions.f90 @@ -6,7 +6,7 @@ !--------------------------------------------------------------------------! module HIIRegion ! -! HIIRegion +! Feedback from HII regions ! ! :References: Fujii et al. (2021), Hopkins et al. (2012) ! @@ -17,7 +17,6 @@ module HIIRegion ! :Dependencies: dim, eos, infile_utils, io, linklist, part, physcon, ! sortutils, timing, units ! - implicit none public :: update_ionrates,update_ionrate, HII_feedback,initialize_H2R,read_options_H2R,write_options_H2R @@ -49,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 @@ -82,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 !----------------------------------------------------------------------- @@ -90,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 @@ -101,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)& @@ -119,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. @@ -130,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 @@ -142,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 @@ -153,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. @@ -163,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 @@ -265,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) @@ -287,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 @@ -316,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') @@ -334,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 From 6d079a7a898ccb237305d3fe870aa10440361a26 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 19:58:37 +1000 Subject: [PATCH 26/31] (ptmass) avoid unneccessary cluttering of .in file --- src/main/ptmass.F90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/main/ptmass.F90 b/src/main/ptmass.F90 index a5362c8a0..30964e48d 100644 --- a/src/main/ptmass.F90 +++ b/src/main/ptmass.F90 @@ -2197,7 +2197,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 +2219,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) 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 +2323,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 From 6833c4c9adfd957df2739e008f4a819e0e86ec04 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 19:59:11 +1000 Subject: [PATCH 27/31] (inject_randomwind) fix mdot_type = 2; read/write r_ref from .in file for this case --- src/main/inject_randomwind.f90 | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/src/main/inject_randomwind.f90 b/src/main/inject_randomwind.f90 index d1b7d694c..c55eea947 100644 --- a/src/main/inject_randomwind.f90 +++ b/src/main/inject_randomwind.f90 @@ -43,6 +43,7 @@ module inject 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 !----------------------------------------------------------------------- @@ -116,7 +117,7 @@ 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) ! mdot = in_code_units(mdot_str,ierr) - dmdt = mdot*mdot_func(r,rbody) ! Use rbody as r_ref, currently the softening length of the body + 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? @@ -212,9 +213,14 @@ subroutine write_options_inject(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) - call write_inopt(delta_theta, 'delta_theta', 'standard deviation for the gaussion distribution (random_type=1)', 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 @@ -247,6 +253,9 @@ 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 From fe823e9de188e004e9c04b1d5425c046d43e9c2c Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 20:14:33 +1000 Subject: [PATCH 28/31] (ptmass/viscosity) minor formatting; subroutine purposes added --- src/main/ptmass.F90 | 35 +++++++++++++++++++++++++++++------ src/main/viscosity.f90 | 1 - 2 files changed, 29 insertions(+), 7 deletions(-) diff --git a/src/main/ptmass.F90 b/src/main/ptmass.F90 index 30964e48d..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 @@ -2219,7 +2242,7 @@ 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) - if (gravity .and. icreate_sinks) then + 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 @@ -2333,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/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 !---------------------------------------------------------------- From 53756bd72b9567be1b1c4e047d8a2b87c5abb1c3 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 20:34:27 +1000 Subject: [PATCH 29/31] (cons2prim) improved comments/documentation --- src/main/cons2prim.f90 | 61 ++++++++++++++++++++++++++++-------------- 1 file changed, 41 insertions(+), 20 deletions(-) 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 From 2aa5493a58c262cf1482e2c28e62859dbb599002 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 31 Jul 2024 20:36:48 +1000 Subject: [PATCH 30/31] (cullendehnen) improved comments/documentation --- src/main/cullendehnen.f90 | 6 ++++++ 1 file changed, 6 insertions(+) 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 From c4a2d9f61044b75474b13dae27fe49bf4bcef0d7 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 5 Aug 2024 08:55:03 +1000 Subject: [PATCH 31/31] (datafiles) update docs --- docs/developer-guide/datafiles.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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: ::