From 56cacf707b73cd4367ca1b3392908b99d786e3aa Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 21 Oct 2024 08:17:06 +1100 Subject: [PATCH 01/36] (set_star) pass npartoftype to shift_star via interface instead of via module; adjust subroutine calls to account for this --- src/setup/set_star.f90 | 15 ++++++++------- src/setup/setup_binary.f90 | 2 +- src/setup/setup_grdisc.F90 | 3 ++- src/setup/setup_grtde.f90 | 2 +- 4 files changed, 12 insertions(+), 10 deletions(-) diff --git a/src/setup/set_star.f90 b/src/setup/set_star.f90 index c8a5f5ec5..94e719333 100644 --- a/src/setup/set_star.f90 +++ b/src/setup/set_star.f90 @@ -400,10 +400,11 @@ end subroutine set_stars ! shift star to the desired position and velocity !+ !----------------------------------------------------------------------- -subroutine shift_star(npart,xyz,vxyz,x0,v0,itype,corotate) - use part, only:get_particle_type,set_particle_type,igas,npartoftype +subroutine shift_star(npart,npartoftype,xyz,vxyz,x0,v0,itype,corotate) + use part, only:get_particle_type,set_particle_type,igas use vectorutils, only:cross_product3D - integer, intent(in) :: npart + integer, intent(in) :: npart + integer, intent(inout) :: npartoftype(:) real, intent(inout) :: xyz(:,:),vxyz(:,:) real, intent(in) :: x0(3),v0(3) integer, intent(in), optional :: itype @@ -434,7 +435,7 @@ subroutine shift_star(npart,xyz,vxyz,x0,v0,itype,corotate) if (mytype /= itype+istar_offset) cycle over_parts ! reset type back to gas call set_particle_type(i,igas) - npartoftype(itype+istar_offset) = npartoftype(itype+istar_offset) - 1 + npartoftype(mytype) = npartoftype(mytype) - 1 npartoftype(igas) = npartoftype(igas) + 1 endif xyz(1:3,i) = xyz(1:3,i) + x0(:) @@ -453,13 +454,13 @@ end subroutine shift_star !+ !----------------------------------------------------------------------- subroutine shift_stars(nstar,star,xyzmh_ptmass_in,vxyz_ptmass_in,& - xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,npart,nptmass,corotate) + xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,npart,npartoftype,nptmass,corotate) integer, intent(in) :: nstar,npart type(star_t), intent(in) :: star(nstar) real, intent(inout) :: xyzh(:,:),vxyzu(:,:) real, intent(in) :: xyzmh_ptmass_in(:,:),vxyz_ptmass_in(:,:) real, intent(inout) :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:) - integer, intent(inout) :: nptmass + integer, intent(inout) :: nptmass,npartoftype(:) logical, intent(in), optional :: corotate integer :: i logical :: do_corotate @@ -469,7 +470,7 @@ subroutine shift_stars(nstar,star,xyzmh_ptmass_in,vxyz_ptmass_in,& do i=1,min(nstar,size(xyzmh_ptmass_in(1,:))) if (star(i)%iprofile > 0) then - call shift_star(npart,xyzh,vxyzu,x0=xyzmh_ptmass_in(1:3,i),& + call shift_star(npart,npartoftype,xyzh,vxyzu,x0=xyzmh_ptmass_in(1:3,i),& v0=vxyz_ptmass_in(1:3,i),itype=i,corotate=do_corotate) else nptmass = nptmass + 1 diff --git a/src/setup/setup_binary.f90 b/src/setup/setup_binary.f90 index 0b6c87bd5..f66acfae9 100644 --- a/src/setup/setup_binary.f90 +++ b/src/setup/setup_binary.f90 @@ -136,7 +136,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,& !--place stars into orbit, or add real sink particles if iprofile=0 ! call shift_stars(nstar,star,xyzmh_ptmass_in,vxyz_ptmass_in,xyzh,vxyzu,& - xyzmh_ptmass,vxyz_ptmass,npart,nptmass,corotate=add_spin) + xyzmh_ptmass,vxyz_ptmass,npart,npartoftype,nptmass,corotate=add_spin) ! !--restore options ! diff --git a/src/setup/setup_grdisc.F90 b/src/setup/setup_grdisc.F90 index 2ae3e5427..73d794714 100644 --- a/src/setup/setup_grdisc.F90 +++ b/src/setup/setup_grdisc.F90 @@ -227,7 +227,8 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, ! shift the star to the position of the second body if (star(i)%iprofile > 0) then - call shift_star(npart,xyzh,vxyzu,x0=xyzmh_ptmass_in(:,2),v0=vxyz_ptmass_in(:,2),itype=i) + call shift_star(npart,npartoftype,xyzh,vxyzu,& + x0=xyzmh_ptmass_in(:,2),v0=vxyz_ptmass_in(:,2),itype=i) else nptmass = nptmass + 1 xyzmh_ptmass(:,nptmass) = xyzmh_ptmass_in(:,2) diff --git a/src/setup/setup_grtde.f90 b/src/setup/setup_grtde.f90 index 65df28533..475989343 100644 --- a/src/setup/setup_grtde.f90 +++ b/src/setup/setup_grtde.f90 @@ -216,7 +216,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, print "(a,3f10.3,/)",' Pericentre = ',rp endif - call shift_star(npart,xyzh,vxyzu,x0=xyzstar,v0=vxyzstar) + call shift_star(npart,npartoftype,xyzh,vxyzu,x0=xyzstar,v0=vxyzstar) if (id==master) print "(/,a,i10,/)",' Number of particles setup = ',npart From dc360c8ec92ec2351fd461248f78cbb5909f6e1c Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 30 Oct 2024 08:48:02 +1100 Subject: [PATCH 02/36] (extern_densprofile) add header to density profile output for ease of plotting --- src/main/extern_densprofile.f90 | 93 ++++++--------------------------- 1 file changed, 17 insertions(+), 76 deletions(-) diff --git a/src/main/extern_densprofile.f90 b/src/main/extern_densprofile.f90 index d8fc59c21..2a78573a6 100644 --- a/src/main/extern_densprofile.f90 +++ b/src/main/extern_densprofile.f90 @@ -41,6 +41,7 @@ module extern_densprofile !+ !---------------------------------------------- subroutine densityprofile_force(xi,yi,zi,fxi,fyi,fzi,phi) + use table_utils, only:yinterp real, intent(in) :: xi, yi, zi real, intent(out) :: fxi, fyi, fzi, phi @@ -58,7 +59,7 @@ subroutine densityprofile_force(xi,yi,zi,fxi,fyi,fzi,phi) if (ri2 >= r2surf) then f = fsurf * (r2surf / ri2)**1.5 else - f = yinterp(ri2, r2tab, ftab, ntab) + f = yinterp(ftab(1:ntab),r2tab(1:ntab),ri2) endif fxi = f * xi @@ -136,13 +137,15 @@ subroutine read_rhotab(filename, rsize, rtab, rhotab, nread, polyk, gamma, rhoc, return endif - ! First line: # K gamma rhoc + ! first line, skip header + read(iunit, *,iostat=ierr) + ! second line: # K gamma rhoc read(iunit, *,iostat=ierr) hash,polyk, gamma, rhoc if (ierr /= 0) then call error('extern_densityprofile','Error reading first line of header from '//trim(filename)) return endif - ! Second line: # nentries (number of r density entries in file) + ! third line: # nentries (number of r density entries in file) read(iunit,*,iostat=ierr) hash,nread if (ierr /= 0) then call error('extern_densityprofile','Error reading second line of header from '//trim(filename)) @@ -180,11 +183,14 @@ subroutine write_rhotab(filename, rtab, rhotab, ntab, polyk, gamma, rhoc, ierr) ierr = 0 open(newunit=iunit,file=filename,action='write',status='replace') + ! write header '# r,density' + write(iunit,"(a)") '# r, density' + ! First line: # K gamma rhoc - write(iunit,*) '# ', polyk, gamma, rhoc + write(iunit,"(a,3(g0,1x))") '# ', polyk, gamma, rhoc ! Second line: # nentries (number of r density entries in file) - write(iunit,*) '# ', ntab + write(iunit,"(a,i0)") '# ', ntab ! Loop over 'n' lines: r and density separated by space do i = 1,ntab @@ -207,82 +213,17 @@ subroutine calc_menc(n, r, rho, menc_out, totmass) r2 = r(1:n)**2 r2rho = r2(1:n) * rho(1:n) - if (.false.) then - ! NB: Ensure that this mass calculation is correct if it is to be used (J. Wurster) - ! Use trapezoid, Simpson's and Simpson's 3/8 for first entries then Simpson's for remaining - ! (trapezoidal term has largest order error: avoid using it as part of sum for later terms) - menc(1) = 0. - menc(2) = (r(2)-r(1)) * (r2rho(1) + r2(2) * rho(2)) / 2. - menc(3) = (r(3)-r(1)) * (r2rho(1) + 4.*r2rho(2) + r2rho(3)) / 6. - menc(4) = (r(4)-r(1)) * (r2rho(1) + 3.*r2rho(2) + 3.*r2rho(3) + r2rho(4)) / 8. - totalmass = menc(1) + menc(2) + menc(3) + menc(4) - do i = 5, n - menc(i) = menc(i-2) + (r2rho(i-2) + 4.*r2rho(i-1) + r2rho(i)) * (r(i) - r(i-2)) / 6. - totalmass = totalmass + menc(i) - enddo - menc(:) = 4.0 * pi * menc(:) - totalmass = 4.0 * pi * totalmass - else - menc(1) = 4.0/3.0*pi*r(1)**3 * rho(1) - do i = 2,n - menc(i) = menc(i-1) + 4.0/3.0*pi*(r(i)**3 - r(i-1)**3) * rho(i) - enddo - totalmass = menc(n) - endif + menc(1) = 4.0/3.0*pi*r(1)**3 * rho(1) + do i = 2,n + menc(i) = menc(i-1) + 4.0/3.0*pi*(r(i)**3 - r(i-1)**3) * rho(i) + enddo + totalmass = menc(n) + if (present(menc_out)) menc_out = menc if (present(totmass)) totmass = totalmass end subroutine calc_menc -! Linear 1D interpolation -real function yinterp(x, xtab, ytab, ntab) - real, intent(in) :: x - real, intent(in) :: ytab(:),xtab(:) - integer, intent(in) :: ntab - - integer :: ibelow, iabove - real :: slope - - yinterp = 0. - if (x <= xtab(1)) then - yinterp = ytab(1) - return - elseif (x >= xtab(ntab)) then - yinterp = ytab(ntab) - return - endif - - ibelow = indexbelow(x, xtab, ntab) - iabove = ibelow + 1 - - slope = (ytab(iabove) - ytab(ibelow)) / (xtab(iabove) - xtab(ibelow)) - yinterp = ytab(ibelow) + (x - xtab(ibelow)) * slope - -end function yinterp - -! Find index below value x in monotomic array xtab -pure integer function indexbelow(x, xtab, ntab) - real, intent(in) :: x - real, intent(in) :: xtab(:) - integer, intent(in) :: ntab - - integer :: ibelow, imid, iabove - - ibelow = 1 - iabove = ntab - do while ((iabove - ibelow) > 1) - imid = (iabove + ibelow) / 2 - if ((x > xtab(imid)) .eqv. (xtab(ntab) > xtab(1))) then - ibelow = imid - else - iabove = imid - endif - enddo - - indexbelow = ibelow - -end function indexbelow - !---------------------------------------------- !+ ! Wrapper to get the density profile (J.Wurster) From 9cb7869dfe26079b30e6448e78d5d6f92a4bfedb Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 30 Oct 2024 08:48:08 +1100 Subject: [PATCH 03/36] (extern_densprofile) add header to density profile output for ease of plotting --- data/neutronstar/ns-rdensity.tab | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/data/neutronstar/ns-rdensity.tab b/data/neutronstar/ns-rdensity.tab index 9e736f9e6..61ab9ef2e 100644 --- a/data/neutronstar/ns-rdensity.tab +++ b/data/neutronstar/ns-rdensity.tab @@ -1,5 +1,6 @@ - # 0.42442000000000002 1.6666700000000001 1.4249734000000001 - # 914 +# r, density +# 0.42442000000000002 1.6666700000000001 1.4249734000000001 +# 914 0.0000000000000000 1.4281200250341688 1.0952902519167577E-003 1.4281143125578766 2.1905805038335154E-003 1.4280971752432505 From 530a6b7292425e5390c5280bf949375677d00c4a Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 30 Oct 2024 08:48:35 +1100 Subject: [PATCH 04/36] (sort) fix warning in sort-by-radius function --- src/main/utils_sort.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/main/utils_sort.f90 b/src/main/utils_sort.f90 index 0ba6429a8..4ffbf6bd3 100644 --- a/src/main/utils_sort.f90 +++ b/src/main/utils_sort.f90 @@ -650,9 +650,9 @@ end subroutine find_rank ! coordinates as input !+ !---------------------------------------------------------------- -subroutine sort_by_radius(n,xyz,iorder,x0) +subroutine sort_by_radius(n,xyzh,iorder,x0) integer, intent(in) :: n - real, intent(in) :: xyz(3,n) + real, intent(in) :: xyzh(4,n) integer, intent(out) :: iorder(n) real, intent(in), optional :: x0(3) @@ -660,7 +660,7 @@ subroutine sort_by_radius(n,xyz,iorder,x0) if (present(x0)) call set_r2func_origin(x0(1),x0(2),x0(3)) ! sort by r^2 using the r2func function - call indexxfunc(n,r2func,xyz,iorder) + call indexxfunc(n,r2func,xyzh,iorder) end subroutine sort_by_radius From 87a099bde77d3e14858d2dad250e52446b15efff Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 30 Oct 2024 09:58:49 +1100 Subject: [PATCH 05/36] (utils_infiles) allow optional default= argument for read_inopt_string --- src/main/utils_infiles.f90 | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/main/utils_infiles.f90 b/src/main/utils_infiles.f90 index 1a5259b6c..8b788ff80 100644 --- a/src/main/utils_infiles.f90 +++ b/src/main/utils_infiles.f90 @@ -528,12 +528,13 @@ end subroutine read_inopt_real ! read a string variable from an input options database !+ !----------------------------------------------------------------- -subroutine read_inopt_string(valstring,tag,db,err,errcount) +subroutine read_inopt_string(valstring,tag,db,err,errcount,default) character(len=*), intent(out) :: valstring character(len=*), intent(in) :: tag type(inopts), allocatable, intent(inout) :: db(:) integer, intent(out), optional :: err integer, intent(inout), optional :: errcount + character(len=*), intent(in), optional :: default integer :: ierr ierr = 0 @@ -546,6 +547,12 @@ subroutine read_inopt_string(valstring,tag,db,err,errcount) if (present(errcount)) then if (ierr /= 0) errcount = errcount + 1 endif + ! default string to use if the string read is blank + if (present(default)) then + if (len_trim(valstring) <= 0) then + valstring = default + endif + endif end subroutine read_inopt_string From 3d5b3b9e59a1562db8bd6da4ee6674f0b64177b9 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 30 Oct 2024 10:01:44 +1100 Subject: [PATCH 06/36] (density_profiles) try to handle the case where stellar radius is unresolved by the grid; adjust dr until resolved --- src/setup/density_profiles.f90 | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/setup/density_profiles.f90 b/src/setup/density_profiles.f90 index b06e9b2f1..6189b73bc 100644 --- a/src/setup/density_profiles.f90 +++ b/src/setup/density_profiles.f90 @@ -163,6 +163,7 @@ subroutine rho_piecewise_polytrope(rtab,rhotab,rhocentre,mstar_in,get_dPdrho,npt lastsign = 1 iterate = .true. bisect = .false. + rtab = 0. ! !--Iterate to get the correct density profile do while ( iterate ) @@ -171,6 +172,10 @@ subroutine rho_piecewise_polytrope(rtab,rhotab,rhocentre,mstar_in,get_dPdrho,npt !--did not complete the profile; reset dr dr = 2.0*dr ierr = 0 + elseif (npts < size(rtab)/4) then + !--profile is unresolved by radial grid, take smaller dr + dr = 0.5*dr + ierr = 0 else call calc_mass_enc(npts,rtab,rhotab,mstar=mstar) !--iterate to get the correct mass @@ -234,6 +239,8 @@ subroutine integrate_rho_profile(rtab,rhotab,rhocentre,get_dPdrho,dr,npts,ierr) i = i + 1 rhotab(i) = rhotab(i-1) + dr*drhodr rtab(i) = rtab(i-1) + dr + if (rhotab(i) < 0.0) exit + dPdrho = get_dPdrho(rhotab(i)) if (i==2) then drhodr = drhodr - fourpi*rhotab(i-1)**2*dr/dPdrho @@ -243,7 +250,6 @@ subroutine integrate_rho_profile(rtab,rhotab,rhocentre,get_dPdrho,dr,npts,ierr) - (dPdrho-dPdrho_prev)/(dr*dPdrho)*drhodr - 2.0*drhodr/rtab(i) ) endif dPdrho_prev = dPdrho - if (rhotab(i) < 0.0) iterate = .false. if (i >=size(rtab)) then ierr = 1 iterate = .false. From b0f68b315ed042adb718704779b260b711b47294 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 30 Oct 2024 10:02:12 +1100 Subject: [PATCH 07/36] (infile) allow eos=9 even if maxvxyzu=4 --- src/main/readwrite_infile.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/readwrite_infile.F90 b/src/main/readwrite_infile.F90 index 1094e6a6a..e173b7fe7 100644 --- a/src/main/readwrite_infile.F90 +++ b/src/main/readwrite_infile.F90 @@ -702,7 +702,7 @@ subroutine read_infile(infile,logfile,evfile,dumpfile) #ifndef MCFOST if (maxvxyzu >= 4 .and. (ieos /= 2 .and. ieos /= 5 .and. ieos /= 4 .and. ieos /= 10 .and. & ieos /=11 .and. ieos /=12 .and. ieos /= 15 .and. ieos /= 16 .and. ieos /= 17 .and. & - ieos /= 20 .and. ieos/=22)) & + ieos /= 20 .and. ieos/=22 .and. ieos /= 9)) & call fatal(label,'only ieos=2 makes sense if storing thermal energy') #endif if (irealvisc < 0 .or. irealvisc > 12) call fatal(label,'invalid setting for physical viscosity') From b78643c9d17ac9014b2b9b2f5d67bde3ccbebdf8 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 30 Oct 2024 10:02:45 +1100 Subject: [PATCH 08/36] (eos) allow import of irecomb from eos module --- src/main/eos.f90 | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 10528d8ca..07e0f7ad9 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -46,8 +46,9 @@ module eos ! eos_stratified, infile_utils, io, mesa_microphysics, part, physcon, ! units ! - use part, only:ien_etotal,ien_entropy,ien_type - use dim, only:gr + use part, only:ien_etotal,ien_entropy,ien_type + use dim, only:gr + use eos_gasradrec, only:irecomb implicit none integer, parameter, public :: maxeos = 22 real, public :: polyk, polyk2, gamma @@ -64,6 +65,8 @@ module eos public :: init_eos,finish_eos,write_options_eos,read_options_eos public :: write_headeropts_eos, read_headeropts_eos + public :: irecomb ! propagated from eos_gasradrec + private integer, public :: ieos = 1 From 77cacb8c4202f27950975fb1781f1b43e98b40dd Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Wed, 30 Oct 2024 10:04:03 +1100 Subject: [PATCH 09/36] (relax_star) bug fixes for relaxation with ieos=9 --- src/setup/relax_star.f90 | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/src/setup/relax_star.f90 b/src/setup/relax_star.f90 index f081f1f4f..c05c2deaa 100644 --- a/src/setup/relax_star.f90 +++ b/src/setup/relax_star.f90 @@ -70,7 +70,7 @@ subroutine relax_star(nt,rho,pr,r,npart,xyzh,use_var_comp,Xfrac,Yfrac,mu,ierr,np use io, only:error,warning,fatal,id,master use fileutils, only:getnextfilename use readwrite_dumps, only:write_fulldump,init_readwrite_dumps - use eos, only:gamma,eos_outputs_mu + use eos, only:gamma,eos_outputs_mu,ieos use physcon, only:pi use options, only:iexternalforce use io_summary, only:summary_initialise @@ -165,9 +165,14 @@ subroutine relax_star(nt,rho,pr,r,npart,xyzh,use_var_comp,Xfrac,Yfrac,mu,ierr,np ! define utherm(r) based on P(r) and rho(r) ! and use this to set the thermal energy of all particles ! - entrop = pr/rho**gamma - utherm = pr/(rho*(gamma-1.)) - if (any(utherm <= 0.)) then + where (rho > 0) + entrop = pr/rho**gamma + utherm = pr/(rho*(gamma-1.)) + elsewhere + entrop = 0. + utherm = 0. + end where + if (any(utherm(1:nt-1) <= 0.)) then call error('relax_star','relax-o-matic needs non-zero pressure array set in order to work') call restore_original_options(i1,npart) ierr = ierr_no_pressure @@ -187,7 +192,7 @@ subroutine relax_star(nt,rho,pr,r,npart,xyzh,use_var_comp,Xfrac,Yfrac,mu,ierr,np ! ! perform sanity checks ! - if (etherm > abs(epot)) then + if (etherm > abs(epot) .and. ieos /= 9) then call error('relax_star','cannot relax star because it is unbound (etherm > epot)') if (id==master) print*,' Etherm = ',etherm,' Epot = ',Epot if (maxvxyzu < 4) print "(/,a,/)",' *** Try compiling with ISOTHERMAL=no instead... ***' @@ -463,7 +468,7 @@ subroutine set_options_for_relaxation(tdyn) ! ! turn on settings appropriate to relaxation ! - if (maxvxyzu >= 4) ieos = 2 + if (maxvxyzu >= 4 .and. ieos /= 9) ieos = 2 if (tdyn > 0.) then idamp = 2 tdyn_s = tdyn*utime From 911ad654483474bf718e4cef92b21a83f5a89162 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Thu, 28 Nov 2024 11:39:12 +1100 Subject: [PATCH 10/36] (set_star) allow arbitrary units in the .setup file for radii, masses and orbital properties when setting up stars and binary systems; shift equation of state options for stars inside set_stars options routines --- src/main/units.f90 | 220 +++++++++++++------ src/setup/set_orbit.f90 | 81 +++---- src/setup/set_star.f90 | 394 +++++++++++++++++++++++------------ src/setup/set_star_utils.f90 | 18 +- src/setup/set_units.f90 | 4 +- src/setup/setup_binary.f90 | 41 ++-- src/setup/setup_star.f90 | 204 +++++++++--------- 7 files changed, 597 insertions(+), 365 deletions(-) diff --git a/src/main/units.f90 b/src/main/units.f90 index 8af9b5dc9..9aa59d101 100644 --- a/src/main/units.f90 +++ b/src/main/units.f90 @@ -31,12 +31,16 @@ module units real(kind=8), public :: unit_ergg, unit_energ, unit_opacity, unit_luminosity real(kind=8), public :: unit_angmom - public :: set_units, set_units_extra, print_units + public :: set_units, set_units_extra, print_units, select_unit public :: get_G_code, get_c_code, get_radconst_code, get_kbmh_code - public :: c_is_unity, G_is_unity, in_geometric_units + public :: c_is_unity, G_is_unity, in_geometric_units, in_code_units, in_units public :: is_time_unit, is_length_unit, is_mdot_unit public :: in_solarr, in_solarm, in_solarl + integer, parameter :: len_utype = 10 + + private + contains !------------------------------------------------------------------------------------ @@ -178,12 +182,14 @@ end subroutine print_units ! Subroutine to recognise mass, length and time units from a string !+ !------------------------------------------------------------------------------------ -subroutine select_unit(string,unit,ierr) +subroutine select_unit(string,unit,ierr,unit_type) use physcon - character(len=*), intent(in) :: string - real(kind=8), intent(out) :: unit - integer, intent(out) :: ierr + character(len=*), intent(in) :: string + real(kind=8), intent(out) :: unit + integer, intent(out) :: ierr + character(len=len_utype), intent(out), optional :: unit_type character(len=len(string)) :: unitstr + character(len=len_utype) :: utype real(kind=8) :: fac ierr = 0 @@ -191,53 +197,105 @@ subroutine select_unit(string,unit,ierr) select case(trim(unitstr)) case('solarr','rsun') - unit = solarr + unit = solarr + utype = 'length' + case('jupiterr','rjup','rjupiter') + unit = jupiterr + utype = 'length' case('au') - unit = au + unit = au + utype = 'length' case('ly','lightyear') - unit = ly + unit = ly + utype = 'length' case('pc','parsec') - unit = pc + unit = pc + utype = 'length' case('kpc','kiloparsec') - unit = kpc + unit = kpc + utype = 'length' case('mpc','megaparsec') - unit = mpc + unit = mpc + utype = 'length' case('km','kilometres','kilometers') - unit = km + unit = km + utype = 'length' case('cm','centimetres','centimeters') unit = 1.d0 + utype = 'length' case('solarm','msun') unit = solarm + utype = 'mass' case('earthm','mearth') unit = earthm + utype = 'mass' case('jupiterm','mjup','mjupiter') unit = jupiterm + utype = 'mass' + case('ceresm','mceres') + unit = ceresm + utype = 'mass' case('g','grams') - unit = 1.d0 + unit = gram + utype = 'mass' case('days','day') unit = days + utype = 'time' case('Myr') unit = 1.d6*years + utype = 'time' case('yr','year','yrs','years') unit = years + utype = 'time' case('hr','hour','hrs','hours') unit = hours + utype = 'time' case('min','minute','mins','minutes') unit = minutes + utype = 'time' 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") + utype = 'time' + case('g/s','grams/second','g/second','grams/s','g/sec','grams/sec') + unit = gram/seconds + utype = 'mdot' + 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','msun/yr') unit = solarm/years + utype = 'mdot' + case('lsun','solarl','Lsun') + unit = solarl + utype = 'luminosity' + case('erg/s') + unit = 1.d0 + utype = 'luminosity' + case('cm/s') + unit = cm/seconds + utype = 'velocity' + case('m/s') + unit = 1.d2*cm/seconds + utype = 'velocity' + case('km/s') + unit = km/seconds + utype = 'velocity' + case('km/h') + unit = km/hours + utype = 'velocity' + case('au/yr') + unit = au/years + utype = 'velocity' + case('c') + unit = c + utype = 'velocity' case default - ierr = 1 + if (len_trim(unitstr) > 0) ierr = 1 unit = 1.d0 + utype = 'none' end select unit = unit*fac + if (present(unit_type)) unit_type = utype end subroutine select_unit @@ -248,21 +306,14 @@ end subroutine select_unit !------------------------------------------------------------------------------------ logical function is_time_unit(string) character(len=*), intent(in) :: string - character(len=len(string)) :: unitstr - real(kind=8) :: fac + character(len=len_utype) :: unit_type + real :: val integer :: ierr ierr = 0 - call get_unit_multiplier(string,unitstr,fac,ierr) + call select_unit(string,val,ierr,unit_type) - select case(trim(unitstr)) - case('days','day','Myr','yr','year','yrs','years',& - 'hr','hour','hrs','hours','min','minute','mins','minutes',& - 's','sec','second','seconds') - is_time_unit = .true. - case default - is_time_unit = .false. - end select + is_time_unit = (trim(unit_type) == 'time') end function is_time_unit @@ -273,21 +324,14 @@ end function is_time_unit !------------------------------------------------------------------------------------ logical function is_length_unit(string) character(len=*), intent(in) :: string - character(len=len(string)) :: unitstr - real(kind=8) :: fac + character(len=len_utype) :: unit_type + real :: val integer :: ierr ierr = 0 - call get_unit_multiplier(string,unitstr,fac,ierr) + call select_unit(string,val,ierr,unit_type) - select case(trim(unitstr)) - case('solarr','rsun','au','ly','lightyear','pc','parsec',& - 'kpc','kiloparsec','mpc','megaparsec','km','kilometres',& - 'kilometers','cm','centimetres','centimeters') - is_length_unit = .true. - case default - is_length_unit = .false. - end select + is_length_unit = (trim(unit_type) == 'length') end function is_length_unit @@ -298,45 +342,59 @@ end function is_length_unit !------------------------------------------------------------------------------------ logical function is_mdot_unit(string) character(len=*), intent(in) :: string - character(len=len(string)) :: unitstr - real(kind=8) :: fac + character(len=len_utype) :: unit_type + real :: val integer :: ierr ierr = 0 - call get_unit_multiplier(string,unitstr,fac,ierr) + call select_unit(string,val,ierr,unit_type) - 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 + is_mdot_unit = (trim(unit_type) == 'mdot') end function is_mdot_unit !------------------------------------------------------------------------------------ !+ -! parse a string like "10.*days" or "10*au" and return the value in code units +! parse a string like '10.*days' or '10*au' and return the value in code units ! if there is no recognisable units, the value is returned unscaled !+ !------------------------------------------------------------------------------------ -real function in_code_units(string,ierr) result(rval) +real function in_code_units(string,ierr,unit_type) result(rval) character(len=*), intent(in) :: string integer, intent(out) :: ierr + character(len=*), intent(in), optional :: unit_type real(kind=8) :: val + character(len=len_utype) :: utype + + call select_unit(string,val,ierr,unit_type=utype) + + ! return an error if incorrect dimensions (e.g. mass instead of length) + if (present(unit_type)) then + if ((trim(utype) /= 'none') .and. trim(utype) /= trim(unit_type)) then + ierr = 2 + rval = real(val) + return + endif + endif - call select_unit(string,val,ierr) - if (is_time_unit(string) .and. ierr == 0) then - 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)) + if (ierr /= 0) then + rval = real(val) + return else - rval = real(val) ! no unit conversion + select case(trim(utype)) + case('time') + rval = real(val/utime) + case('length') + rval = real(val/udist) + case('mass') + rval = real(val/umass) + case('mdot') + rval = real(val/(umass/utime)) + case('luminosity') + rval = real(val/unit_luminosity) + case default + rval = real(val) ! no unit conversion + end select endif end function in_code_units @@ -390,7 +448,7 @@ end subroutine get_unit_multiplier pure logical function is_digit(ch) character(len=1), intent(in) :: ch - is_digit = (iachar(ch) >= iachar('0') .and. iachar(ch) <= iachar('9')) + is_digit = (iachar(ch) >= iachar('0') .and. iachar(ch) <= iachar('9')) .or. (ch=='.') end function is_digit @@ -511,4 +569,36 @@ real(kind=8) function in_solarl(val) result(rval) end function in_solarl +!------------------------------------------------------------------------------------ +!+ +! print a value in physical units, e.g. give code value of mass and +! call this routine print*,in_units(mass,'solarm') +!+ +!------------------------------------------------------------------------------------ +real(kind=8) function in_units(val,unitstring) result(rval) + real, intent(in) :: val + character(len=*), intent(in) :: unitstring + character(len=len_utype) :: utype + integer :: ierr + real :: fac + + call select_unit(unitstring,fac,ierr,unit_type=utype) ! handle errors silently by ignoring ierr + + select case(trim(utype)) + case('time') + rval = fac*val/utime + case('length') + rval = fac*val/udist + case('mass') + rval = fac*val/umass + case('mdot') + rval = fac*val/(umass/utime) + case('luminosity') + rval = fac*val/unit_luminosity + case default + rval = real(fac*val) ! no unit conversion + end select + +end function in_units + end module units diff --git a/src/setup/set_orbit.f90 b/src/setup/set_orbit.f90 index a44a77ff1..605a8a467 100644 --- a/src/setup/set_orbit.f90 +++ b/src/setup/set_orbit.f90 @@ -36,8 +36,10 @@ module setorbit ! define data types with options needed ! to setup an orbit ! + integer, parameter :: len_str = 20 + type campbell_elems - character(len=20) :: semi_major_axis ! string because can specific units + character(len=len_str) :: semi_major_axis ! string because can specific units real :: e ! eccentricity real :: i ! inclination real :: O ! position angle of the ascending node @@ -46,15 +48,15 @@ module setorbit end type campbell_elems type posvel_elems - real :: x1(3) ! position of body 1 - real :: v1(3) ! velocity of body 1 - real :: x2(3) ! position of body 2 - real :: v2(3) ! velocity of body 2 + character(len=len_str) :: x1(3) ! position of body 1 + character(len=len_str) :: v1(3) ! velocity of body 1 + character(len=len_str) :: x2(3) ! position of body 2 + character(len=len_str) :: v2(3) ! velocity of body 2 end type posvel_elems type flyby_elems character(len=20) :: rp ! pericentre distance in arbitrary units - real :: d ! initial separation + character(len=20) :: d ! initial separation in arbitrary units real :: O ! position angle of the ascending node real :: i ! inclination end type flyby_elems @@ -90,18 +92,18 @@ subroutine set_defaults_orbit(orbit) orbit%elems%f = 180. ! start orbit at apocentre orbit%flyby%rp = '10.' - orbit%flyby%d = 100.0 + orbit%flyby%d = '100.0' orbit%flyby%O = 0.0 orbit%flyby%i = 0.0 - orbit%posvel%x1 = 0.0 - orbit%posvel%v1 = 0.0 - orbit%posvel%x2 = 0.0 - orbit%posvel%v2 = 0.0 - orbit%posvel%x1(1) = 10.0 - orbit%posvel%x2(1) = -10.0 - orbit%posvel%v1(2) = 1.0 - orbit%posvel%v2(2) = -1.0 + orbit%posvel%x1 = '0.0' + orbit%posvel%v1 = '0.0' + orbit%posvel%x2 = '0.0' + orbit%posvel%v2 = '0.0' + orbit%posvel%x1(1) = '10.0' + orbit%posvel%x2(1) = '-10.0' + orbit%posvel%v1(2) = '1.0' + orbit%posvel%v2(2) = '-1.0' end subroutine set_defaults_orbit @@ -122,25 +124,33 @@ subroutine set_orbit(orbit,m1,m2,hacc1,hacc2,xyzmh_ptmass,vxyz_ptmass,nptmass,ve logical, intent(in) :: verbose integer, intent(out) :: ierr real, intent(out), optional :: omega_corotate - real :: rp,a + real :: rp,d,a,x1(3),x2(3),v1(3),v2(3) + integer :: i ierr = 0 select case(orbit%itype) case(2) + do i=1,3 + x1(i) = in_code_units(orbit%posvel%x1(i),ierr,unit_type='length') + x2(i) = in_code_units(orbit%posvel%x2(i),ierr,unit_type='length') + v1(i) = in_code_units(orbit%posvel%v1(i),ierr,unit_type='velocity') + v2(i) = in_code_units(orbit%posvel%v2(i),ierr,unit_type='velocity') + enddo ! body 1 - xyzmh_ptmass(1:3,nptmass+1) = orbit%posvel%x1(1:3) + xyzmh_ptmass(1:3,nptmass+1) = x1 xyzmh_ptmass(4,nptmass+1) = m1 xyzmh_ptmass(5,nptmass+1) = hacc1 - vxyz_ptmass(1:3,nptmass+1) = orbit%posvel%v1(1:3) + vxyz_ptmass(1:3,nptmass+1) = v1 ! body 2 - xyzmh_ptmass(1:3,nptmass+2) = orbit%posvel%x2(1:3) + xyzmh_ptmass(1:3,nptmass+2) = x2 xyzmh_ptmass(4,nptmass+2) = m2 xyzmh_ptmass(5,nptmass+2) = hacc2 - vxyz_ptmass(1:3,nptmass+2) = orbit%posvel%v2(1:3) + vxyz_ptmass(1:3,nptmass+2) = v2 case(1) rp = in_code_units(orbit%flyby%rp,ierr) + d = in_code_units(orbit%flyby%d,ierr) - call set_flyby(m1,m2,rp,orbit%flyby%d,hacc1,hacc2,xyzmh_ptmass, & + call set_flyby(m1,m2,rp,d,hacc1,hacc2,xyzmh_ptmass, & vxyz_ptmass,nptmass,ierr,orbit%flyby%O,orbit%flyby%i,verbose=verbose) case default ! @@ -190,21 +200,21 @@ subroutine write_options_orbit(orbit,iunit,label) call write_inopt(orbit%itype,'itype'//trim(c),'type of orbital elements (0=aeiOwf,1=flyby,2=posvel)',iunit) select case(orbit%itype) case(2) - call write_inopt(orbit%posvel%x1(1),'x1'//trim(c),'x position body 1',iunit) - call write_inopt(orbit%posvel%x1(2),'y1'//trim(c),'y position body 1',iunit) - call write_inopt(orbit%posvel%x1(3),'z1'//trim(c),'z position body 1',iunit) - call write_inopt(orbit%posvel%v1(1),'vx1'//trim(c),'x velocity body 1',iunit) - call write_inopt(orbit%posvel%v1(2),'vy1'//trim(c),'y velocity body 1',iunit) - call write_inopt(orbit%posvel%v1(3),'vz1'//trim(c),'z velocity body 1',iunit) - call write_inopt(orbit%posvel%x2(1),'x2'//trim(c),'x position body 2',iunit) - call write_inopt(orbit%posvel%x2(2),'y2'//trim(c),'y position body 2',iunit) - call write_inopt(orbit%posvel%x2(3),'z2'//trim(c),'z position body 2',iunit) - call write_inopt(orbit%posvel%v2(1),'vx2'//trim(c),'x velocity body 2',iunit) - call write_inopt(orbit%posvel%v2(2),'vy2'//trim(c),'y velocity body 2',iunit) - call write_inopt(orbit%posvel%v2(3),'vz2'//trim(c),'z velocity body 2',iunit) + call write_inopt(orbit%posvel%x1(1),'x1'//trim(c),'x position body 1 (code units or e.g. 1*au)',iunit) + call write_inopt(orbit%posvel%x1(2),'y1'//trim(c),'y position body 1 (code units or e.g. 1*au)',iunit) + call write_inopt(orbit%posvel%x1(3),'z1'//trim(c),'z position body 1 (code units or e.g. 1*au)',iunit) + call write_inopt(orbit%posvel%v1(1),'vx1'//trim(c),'x velocity body 1 (code units or e.g. 1*km/s)',iunit) + call write_inopt(orbit%posvel%v1(2),'vy1'//trim(c),'y velocity body 1 (code units or e.g. 1*km/s)',iunit) + call write_inopt(orbit%posvel%v1(3),'vz1'//trim(c),'z velocity body 1 (code units or e.g. 1*km/s)',iunit) + call write_inopt(orbit%posvel%x2(1),'x2'//trim(c),'x position body 2 (code units or e.g. 1*au)',iunit) + call write_inopt(orbit%posvel%x2(2),'y2'//trim(c),'y position body 2 (code units or e.g. 1*au)',iunit) + call write_inopt(orbit%posvel%x2(3),'z2'//trim(c),'z position body 2 (code units or e.g. 1*au)',iunit) + call write_inopt(orbit%posvel%v2(1),'vx2'//trim(c),'x velocity body 2 (code units or e.g. 1*km/s)',iunit) + call write_inopt(orbit%posvel%v2(2),'vy2'//trim(c),'y velocity body 2 (code units or e.g. 1*km/s)',iunit) + call write_inopt(orbit%posvel%v2(3),'vz2'//trim(c),'z velocity body 2 (code units or e.g. 1*km/s)',iunit) case(1) - call write_inopt(orbit%flyby%rp,'rp'//trim(c),'pericentre distance',iunit) - call write_inopt(orbit%flyby%d,'d'//trim(c),'initial separation [same units as rp]',iunit) + call write_inopt(orbit%flyby%rp,'rp'//trim(c),'pericentre distance (code units or e.g. 1*au)',iunit) + call write_inopt(orbit%flyby%d,'d'//trim(c),'initial separation (code units or e.g. 1*au)',iunit) call write_inopt(orbit%flyby%O,'O'//trim(c),'position angle of the ascending node',iunit) call write_inopt(orbit%flyby%i,'i'//trim(c),'inclination',iunit) case default @@ -236,6 +246,7 @@ subroutine read_options_orbit(orbit,db,nerr,label) c = '' if (present(label)) c = trim(adjustl(label)) + call set_defaults_orbit(orbit) call read_inopt(orbit%itype,'itype'//trim(c),db,errcount=nerr,min=0,max=2) select case(orbit%itype) case(2) diff --git a/src/setup/set_star.f90 b/src/setup/set_star.f90 index 94e719333..8997a9cb0 100644 --- a/src/setup/set_star.f90 +++ b/src/setup/set_star.f90 @@ -8,8 +8,8 @@ module setstar ! ! General routine for setting up a 3D star from a 1D profile ! This is the main functionality from setup_star but in a single routine -! that can also be called from other setups. In principle this -! could also be used to setup multiple stars +! that can also be called from other setups. Also contains +! routines to setup multiple stars ! ! :References: None ! @@ -24,7 +24,7 @@ module setstar ! setstar_utils, unifdis, units, vectorutils ! use setstar_utils, only:ikepler,imesa,ibpwpoly,ipoly,iuniform,ifromfile,ievrard,& - need_polyk + need_polyk,need_mu implicit none ! @@ -37,15 +37,10 @@ module setstar logical :: isinkcore integer :: isofteningopt integer :: np - real :: Rstar - real :: Mstar + character(len=20) :: m,r,rcore,mcore,lcore,hsoft,hacc real :: ui_coef + real :: polyk real :: initialtemp - real :: rcore - real :: mcore - real :: lcore - real :: hsoft - real :: hacc ! accretion radius if star is a sink particle character(len=120) :: input_profile,dens_profile character(len=120) :: outputfilename ! outputfilename is the path to the cored profile character(len=2) :: label ! used to rename relax_star snapshots to relax1, relax2 etc. @@ -64,6 +59,8 @@ module setstar integer, parameter :: istar_offset = 3 ! offset for particle type to distinguish particles ! placed in stars from other particles in the simulation + integer, private :: EOSopt = 1 + private contains @@ -75,22 +72,21 @@ module setstar !+ !-------------------------------------------------------------------------- subroutine set_defaults_star(star) - use units, only:udist,umass - use physcon, only:solarm,solarr type(star_t), intent(out) :: star star%iprofile = 2 - star%rstar = 1.0*real(solarr/udist) - star%mstar = 1.0*real(solarm/umass) + star%r = '1.0' + star%m = '1.0' star%ui_coef = 0.05 + star%polyk = 0. star%initialtemp = 1.0e7 star%isoftcore = 0 star%isinkcore = .false. - star%hsoft = 0. - star%hacc = 1. - star%rcore = 0. - star%mcore = 0. - star%lcore = 0. + star%hsoft = '0.0' + star%hacc = '1.0' + star%rcore = '0.0' + star%mcore = '0.0' + star%lcore = '0.*lsun' star%isofteningopt = 1 ! By default, specify rcore star%np = 1000 star%input_profile = 'P12_Phantom_Profile.data' @@ -115,6 +111,48 @@ subroutine set_defaults_stars(stars) end subroutine set_defaults_stars +!-------------------------------------------------------------------------- +!+ +! utility routine to convert properties to code units while checking +! for errors +!+ +!-------------------------------------------------------------------------- +subroutine check_and_convert(var,desc,unit_type,val,nerr) + use units, only:in_code_units + use io, only:error + character(len=*), intent(in) :: var,desc,unit_type + real, intent(out) :: val + integer, intent(inout) :: nerr + integer :: ierr + + val = in_code_units(var,ierr,unit_type=unit_type) + if (ierr /= 0) then + call error('set_star','error parsing units for '//trim(desc),var=var) + nerr = nerr + 1 + endif + +end subroutine check_and_convert + +!-------------------------------------------------------------------------- +!+ +! utility routine to convert the stellar properties to code units +!+ +!-------------------------------------------------------------------------- +subroutine get_star_properties_in_code_units(star,rstar,mstar,rcore,mcore,hsoft,lcore,hacc,nerr) + type(star_t), intent(in) :: star + real, intent(out) :: rstar,mstar,rcore,mcore,hsoft,lcore,hacc + integer, intent(inout) :: nerr + + call check_and_convert(star%r,'stellar radius','length',rstar,nerr) + call check_and_convert(star%m,'stellar mass','mass',mstar,nerr) + call check_and_convert(star%rcore,'rcore','length',rcore,nerr) + call check_and_convert(star%mcore,'mcore','mass',mcore,nerr) + call check_and_convert(star%hsoft,'core softening','length',hsoft,nerr) + call check_and_convert(star%lcore,'core luminosity','luminosity',lcore,nerr) + call check_and_convert(star%hacc,'accretion radius','length',hacc,nerr) + +end subroutine get_star_properties_in_code_units + !-------------------------------------------------------------------------- !+ ! Master routine to setup a star from a specified file or density profile @@ -122,7 +160,7 @@ end subroutine set_defaults_stars !-------------------------------------------------------------------------- subroutine set_star(id,master,star,xyzh,vxyzu,eos_vars,rad,& npart,npartoftype,massoftype,hfact,& - xyzmh_ptmass,vxyz_ptmass,nptmass,ieos,polyk,gamma,X_in,Z_in,& + xyzmh_ptmass,vxyz_ptmass,nptmass,ieos,gamma,X_in,Z_in,& relax,use_var_comp,write_rho_to_file,& rhozero,npart_total,mask,ierr,x0,v0,itype,& write_files,density_error,energy_error) @@ -139,7 +177,7 @@ subroutine set_star(id,master,star,xyzh,vxyzu,eos_vars,rad,& use extern_densprofile, only:write_rhotab use unifdis, only:mask_prototype use physcon, only:pi - use units, only:utime,udist,umass,unit_density + use units, only:umass,udist,utime,unit_density use mpiutils, only:reduceall_mpi type(star_t), intent(inout) :: star integer, intent(in) :: id,master @@ -150,7 +188,7 @@ subroutine set_star(id,master,star,xyzh,vxyzu,eos_vars,rad,& real, intent(in) :: hfact logical, intent(in) :: relax,use_var_comp,write_rho_to_file integer, intent(in) :: ieos - real, intent(inout) :: polyk,gamma + real, intent(inout) :: gamma real, intent(in) :: X_in,Z_in real, intent(out) :: rhozero integer(kind=8), intent(out) :: npart_total @@ -160,11 +198,12 @@ subroutine set_star(id,master,star,xyzh,vxyzu,eos_vars,rad,& logical, intent(in), optional :: write_files real, intent(out), optional :: density_error,energy_error procedure(mask_prototype) :: mask - integer :: npts,ierr_relax + integer :: npts,ierr_relax,nerr integer :: ncols_compo,npart_old,i real, allocatable :: r(:),den(:),pres(:),temp(:),en(:),mtab(:),Xfrac(:),Yfrac(:),mu(:) real, allocatable :: composition(:,:) real :: rmin,rhocentre,rmserr,en_err + real :: rstar,mstar,rcore,mcore,hsoft,lcore,hacc character(len=20), allocatable :: comp_label(:) character(len=30) :: lattice ! The lattice type if stretchmap is used logical :: use_exactN,composition_exists,write_dumps @@ -185,14 +224,23 @@ subroutine set_star(id,master,star,xyzh,vxyzu,eos_vars,rad,& return endif ! + ! perform unit conversion on input quantities (if necessary) + ! + nerr = 0 + call get_star_properties_in_code_units(star,rstar,mstar,rcore,mcore,hsoft,lcore,hacc,nerr) + if (nerr /= 0) then + ierr = 2 + return + endif + ! ! get the desired tables of density, pressure, temperature and composition ! as a function of radius / mass fraction ! - call read_star_profile(star%iprofile,ieos,star%input_profile,gamma,polyk,& + call read_star_profile(star%iprofile,ieos,star%input_profile,gamma,star%polyk,& star%ui_coef,r,den,pres,temp,en,mtab,X_in,Z_in,Xfrac,Yfrac,mu,& - npts,rmin,star%rstar,star%mstar,rhocentre,& - star%isoftcore,star%isofteningopt,star%rcore,star%mcore,& - star%hsoft,star%outputfilename,composition,& + npts,rmin,rstar,mstar,rhocentre,& + star%isoftcore,star%isofteningopt,rcore,mcore,& + hsoft,star%outputfilename,composition,& comp_label,ncols_compo) ! @@ -205,26 +253,26 @@ subroutine set_star(id,master,star,xyzh,vxyzu,eos_vars,rad,& ierr = 2 return endif - if (star%mstar < 0.) then + if (mstar < 0.) then call fatal('set_star','cannot set up a star with negative mass!') ierr = 2 return endif - call set_star_density(lattice,id,master,rmin,star%rstar,star%mstar,hfact,& + call set_star_density(lattice,id,master,rmin,rstar,mstar,hfact,& npts,den,r,npart,npartoftype,massoftype,xyzh,use_exactN,& star%np,rhozero,npart_total,mask) ! ! die if stupid things done with GR ! if (gr) then - if (star%rstar < 6.*star%mstar) call fatal('set_star','R < 6GM/c^2 for star in GR violates weak field assumption') + if (rstar < 6.*mstar) call fatal('set_star','R < 6GM/c^2 for star in GR violates weak field assumption') endif ! ! add sink particle stellar core ! if (star%isinkcore) call set_stellar_core(nptmass,xyzmh_ptmass,vxyz_ptmass,ihsoft,& - star%mcore,star%hsoft,ilum,star%lcore,ierr) + mcore,hsoft,ilum,lcore,ierr) if (ierr==1) call fatal('set_stellar_core','mcore <= 0') if (ierr==2) call fatal('set_stellar_core','hsoft <= 0') if (ierr==3) call fatal('set_stellar_core','lcore < 0') @@ -232,7 +280,7 @@ subroutine set_star(id,master,star,xyzh,vxyzu,eos_vars,rad,& ! Write the desired profile to file (do this before relaxation) ! if (write_rho_to_file) call write_rhotab(star%dens_profile,& - r,den,npts,polyk,gamma,rhocentre,ierr) + r,den,npts,star%polyk,gamma,rhocentre,ierr) ! ! mask any existing particles as accreted so they are ! excluded from the centre of mass and relax_star calculations @@ -243,7 +291,7 @@ subroutine set_star(id,master,star,xyzh,vxyzu,eos_vars,rad,& ! if (relax) then if (reduceall_mpi('+',npart)==npart) then - polyk_eos = polyk + polyk_eos = star%polyk call relax_star(npts,den,pres,r,npart,xyzh,use_var_comp,Xfrac,Yfrac,& mu,ierr_relax,npin=npart_old,label=star%label,& write_dumps=write_dumps,density_error=rmserr,energy_error=en_err) @@ -269,7 +317,7 @@ subroutine set_star(id,master,star,xyzh,vxyzu,eos_vars,rad,& ! if (use_var_comp .or. eos_outputs_mu(ieos)) then call set_star_composition(use_var_comp,eos_outputs_mu(ieos),npart,& - xyzh,Xfrac,Yfrac,mu,mtab,star%mstar,eos_vars,npin=npart_old) + xyzh,Xfrac,Yfrac,mu,mtab,mstar,eos_vars,npin=npart_old) endif ! ! Write .comp file containing composition of each particle after interpolation @@ -327,12 +375,12 @@ subroutine set_star(id,master,star,xyzh,vxyzu,eos_vars,rad,& write(*,'(1x,a,f12.5)') 'gamma = ', gamma endif if (maxvxyzu <= 4 .and. need_polyk(star%iprofile)) then - write(*,'(1x,a,f12.5)') 'polyk = ', polyk - write(*,'(1x,a,f12.6,a)') 'specific int. energ = ', polyk*star%rstar/star%mstar,' GM/R' + write(*,'(1x,a,f12.5)') 'polyk = ', star%polyk + write(*,'(1x,a,f12.6,a)') 'specific int. energ = ', star%polyk*rstar/mstar,' GM/R' endif call write_mass('particle mass = ',massoftype(igas),umass) - call write_dist('Radius = ',star%rstar,udist) - call write_mass('Mass = ',star%mstar,umass) + call write_dist('Radius = ',rstar,udist) + call write_mass('Mass = ',mstar,umass) if (star%iprofile==ipoly) then write(*,'(1x,a,g0,a)') 'rho_central = ', rhocentre*unit_density,' g/cm^3' endif @@ -341,9 +389,7 @@ subroutine set_star(id,master,star,xyzh,vxyzu,eos_vars,rad,& , rhozero, ' code units' write(*,'(1x,a,es12.5,a)') 'free fall time = ', sqrt(3.*pi/(32.*rhozero))*utime,' s' - if (composition_exists) then - write(*,'(a)') 'Composition written to kepler.comp file.' - endif + if (composition_exists) write(*,'(a)') 'Composition written to .comp file.' write(*,"(70('='))") endif @@ -362,10 +408,13 @@ end subroutine set_star !-------------------------------------------------------------------------- subroutine set_stars(id,master,nstars,star,xyzh,vxyzu,eos_vars,rad,& npart,npartoftype,massoftype,hfact,& - xyzmh_ptmass,vxyz_ptmass,nptmass,ieos,polyk,gamma,X_in,Z_in,& + xyzmh_ptmass,vxyz_ptmass,nptmass,ieos,gamma,X_in,Z_in,& relax,use_var_comp,write_rho_to_file,& rhozero,npart_total,mask,ierr) - use unifdis, only:mask_prototype + use unifdis, only:mask_prototype + use eos, only:init_eos,finish_eos + use eos_piecewise, only:init_eos_piecewise_preset + use io, only:error type(star_t), intent(inout) :: star(:) integer, intent(in) :: id,master,nstars integer, intent(inout) :: npart,npartoftype(:),nptmass @@ -375,7 +424,7 @@ subroutine set_stars(id,master,nstars,star,xyzh,vxyzu,eos_vars,rad,& real, intent(in) :: hfact logical, intent(in) :: relax,use_var_comp,write_rho_to_file integer, intent(in) :: ieos - real, intent(inout) :: polyk,gamma + real, intent(inout) :: gamma real, intent(in) :: X_in,Z_in real, intent(out) :: rhozero integer(kind=8), intent(out) :: npart_total @@ -383,16 +432,26 @@ subroutine set_stars(id,master,nstars,star,xyzh,vxyzu,eos_vars,rad,& procedure(mask_prototype) :: mask integer :: i + print*,' HERE ieos=',ieos,EOSopt + if (ieos==9) call init_eos_piecewise_preset(EOSopt) + call init_eos(ieos,ierr) + if (ierr /= 0) then + call error('setup','could not initialise equation of state') + return + endif + do i=1,min(nstars,size(star)) if (star(i)%iprofile > 0) then print "(/,a,i0,a)",' --- STAR ',i,' ---' call set_star(id,master,star(i),xyzh,vxyzu,eos_vars,rad,npart,npartoftype,& - massoftype,hfact,xyzmh_ptmass,vxyz_ptmass,nptmass,ieos,polyk,gamma,& + massoftype,hfact,xyzmh_ptmass,vxyz_ptmass,nptmass,ieos,gamma,& X_in,Z_in,relax,use_var_comp,write_rho_to_file,& rhozero,npart_total,mask,ierr,itype=i) endif enddo + call finish_eos(ieos,ierr) + end subroutine set_stars !----------------------------------------------------------------------- @@ -525,7 +584,8 @@ subroutine set_defaults_given_profile(iprofile,filename,need_iso,ieos,mstar,poly character(len=120), intent(out) :: filename integer, intent(out) :: need_iso integer, intent(inout) :: ieos - real, intent(inout) :: mstar,polyk + real, intent(inout) :: polyk + character(len=20), intent(out) :: mstar need_iso = 0 select case(iprofile) @@ -546,8 +606,7 @@ subroutine set_defaults_given_profile(iprofile,filename,need_iso,ieos,mstar,poly ! Original Author: Madeline Marshall & Bernard Field ! Supervisors: James Wurster & Paul Lasky ieos = 9 - !dist_unit = 'km' - Mstar = 1.35 + Mstar = '1.35' polyk = 144. case(ievrard) need_iso = -1 @@ -560,28 +619,19 @@ end subroutine set_defaults_given_profile ! interactive prompting for setting up a star !+ !----------------------------------------------------------------------- -subroutine set_star_interactive(id,master,star,need_iso,use_var_comp,ieos,polyk) +subroutine set_star_interactive(id,master,star,need_iso,use_var_comp,ieos) use prompting, only:prompt use setstar_utils, only:nprofile_opts,profile_opt,need_inputprofile,need_rstar - use units, only:in_solarm,in_solarr,in_solarl,udist,umass,unit_luminosity - use physcon, only:solarr,solarm,solarl + use units, only:in_code_units type(star_t), intent(out) :: star integer, intent(in) :: id,master logical, intent(out) :: use_var_comp integer, intent(out) :: need_iso integer, intent(inout) :: ieos - real, intent(inout) :: polyk integer :: i - real :: mstar_msun,rstar_rsun,rcore_rsun,mcore_msun,lcore_lsun,hsoft_rsun ! set defaults call set_defaults_star(star) - mstar_msun = real(in_solarm(star%mstar)) - rstar_rsun = real(in_solarr(star%rstar)) - mcore_msun = real(in_solarm(star%mcore)) - lcore_lsun = real(in_solarl(star%lcore)) - rcore_rsun = real(in_solarr(star%rcore)) - hsoft_rsun = real(in_solarr(star%hsoft)) ! Select sphere & set default values do i = 1, nprofile_opts @@ -594,7 +644,7 @@ subroutine set_star_interactive(id,master,star,need_iso,use_var_comp,ieos,polyk) ! if (id==master) write(*,"('Setting up ',a)") trim(profile_opt(star%iprofile)) call set_defaults_given_profile(star%iprofile,star%input_profile,& - need_iso,ieos,star%mstar,polyk) + need_iso,ieos,star%m,star%polyk) ! resolution if (star%iprofile > 0) then @@ -606,11 +656,9 @@ subroutine set_star_interactive(id,master,star,need_iso,use_var_comp,ieos,polyk) if (need_inputprofile(star%iprofile)) then call prompt('Enter file name containing input profile',star%input_profile) else - call prompt('Enter the mass of the star (Msun)',mstar_msun,0.) - star%mstar = mstar_msun*real(solarm/umass) + call prompt('Enter the mass of the star (e.g. 1*msun)',star%m) if (need_rstar(star%iprofile)) then - call prompt('Enter the radius of the star (Rsun)',rstar_rsun,0.) - star%rstar = rstar_rsun*real(solarr/udist) + call prompt('Enter the radius of the star (e.g. 1*rsun)',star%r) endif endif @@ -628,12 +676,9 @@ subroutine set_star_interactive(id,master,star,need_iso,use_var_comp,ieos,polyk) case(0) call prompt('Add a sink particle stellar core?',star%isinkcore) if (star%isinkcore) then - call prompt('Enter mass of the created sink particle core [Msun]',mcore_msun,0.) - call prompt('Enter softening length of the sink particle core [Rsun]',hsoft_rsun,0.) - call prompt('Enter sink particle luminosity [Lsun]',lcore_lsun,0.) - star%mcore = mcore_msun*real(solarm/umass) - star%hsoft = hsoft_rsun*real(solarr/udist) - star%lcore = lcore_lsun*real(solarl/unit_luminosity) + call prompt('Enter mass of the created sink particle core (e.g. 0.1*Msun)',star%mcore) + call prompt('Enter softening length of the sink particle core (e.g. 0.1*Rsun)',star%hsoft) + call prompt('Enter sink particle luminosity (e.g. 1*Lsun)',star%lcore) endif case(1) star%isinkcore = .true. ! Create sink particle core automatically @@ -647,35 +692,27 @@ subroutine set_star_interactive(id,master,star,need_iso,use_var_comp,ieos,polyk) select case(star%isofteningopt) case(1) - call prompt('Enter core radius [Rsun]',rcore_rsun,0.) - star%rcore = rcore_rsun*real(solarr/udist) + call prompt('Enter core radius (e.g. 0.1*rsun)',star%rcore) case(2) - call prompt('Enter mass of the created sink particle core [Msun]',mcore_msun,0.) - star%mcore = mcore_msun*real(solarm/umass) + call prompt('Enter mass of the created sink particle core (e.g. 0.1*msun)',star%mcore) case(3) - call prompt('Enter mass of the created sink particle core [Msun]',mcore_msun,0.) - call prompt('Enter core radius [Rsun]',rcore_rsun,0.) - star%mcore = mcore_msun*real(solarm/umass) - star%rcore = rcore_rsun*real(solarr/udist) + call prompt('Enter mass of the created sink particle core (e.g. 0.1*msun)',star%mcore) + call prompt('Enter core radius (e.g. 1*rsun)',star%rcore) end select - call prompt('Enter sink particle luminosity [Lsun]',lcore_lsun,0.) - star%lcore = lcore_lsun*real(solarl/unit_luminosity) + call prompt('Enter sink particle luminosity [e.g. 1*Lsun]',star%lcore) case(2) star%isinkcore = .true. ! Create sink particle core automatically print*,'Specify core radius and initial guess for mass of sink particle core' - call prompt('Enter core radius in Rsun : ',rcore_rsun,0.) - call prompt('Enter guess for core mass in Msun : ',mcore_msun,0.) - call prompt('Enter sink particle luminosity [Lsun]',lcore_lsun,0.) + call prompt('Enter core radius (e.g. 0.1*rsun): ',star%rcore) + call prompt('Enter guess for core mass (e.g. 0.1*Msun): ',star%mcore) + call prompt('Enter sink particle luminosity (e.g. 1.*lsun',star%lcore) call prompt('Enter output file name of cored stellar profile:',star%outputfilename) - star%mcore = mcore_msun*real(solarm/umass) - star%rcore = rcore_rsun*real(solarr/udist) - star%lcore = lcore_lsun*real(solarl/unit_luminosity) end select case(ievrard) call prompt('Enter the specific internal energy (units of GM/R) ',star%ui_coef,0.) case(:0) - call prompt('Enter the accretion radius in code units',star%hacc,0.) + call prompt('Enter the accretion radius (e.g. 1.0)',star%hacc) end select end subroutine set_star_interactive @@ -685,12 +722,11 @@ end subroutine set_star_interactive ! write setupfile options needed for a star !+ !----------------------------------------------------------------------- -subroutine write_options_star(star,iunit,label) +subroutine write_options_star(star,iunit,ieos,label) use infile_utils, only:write_inopt,get_optstring - use setstar_utils, only:nprofile_opts,profile_opt,need_inputprofile,need_rstar - use units, only:in_solarm,in_solarr,in_solarl + use setstar_utils, only:nprofile_opts,profile_opt,need_inputprofile,need_rstar,need_polyk type(star_t), intent(in) :: star - integer, intent(in) :: iunit + integer, intent(in) :: iunit,ieos character(len=*), intent(in), optional :: label character(len=120) :: string character(len=10) :: c @@ -708,9 +744,9 @@ subroutine write_options_star(star,iunit,label) call write_inopt(star%input_profile,'input_profile'//trim(c),& 'Path to input profile',iunit) else - call write_inopt(in_solarm(star%Mstar),'Mstar'//trim(c),'mass of star '//trim(c)//' [Msun]',iunit) + call write_inopt(star%m,'Mstar'//trim(c),'mass of star '//trim(c)//' (code units or e.g. 1*msun)',iunit) if (need_rstar(star%iprofile)) & - call write_inopt(in_solarr(star%Rstar),'Rstar'//trim(c),'radius of star'//trim(c)//' [Rsun]',iunit) + call write_inopt(star%r,'Rstar'//trim(c),'radius of star'//trim(c)//' (code units or e.g. 1*rsun)',iunit) endif endif @@ -729,31 +765,28 @@ subroutine write_options_star(star,iunit,label) call write_inopt(star%isofteningopt,'isofteningopt'//trim(c),& '1=supply rcore, 2=supply mcore, 3=supply both',iunit) if ((star%isofteningopt == 1) .or. (star%isofteningopt == 3)) then - call write_inopt(in_solarr(star%rcore),'rcore'//trim(c),'Radius of core softening [Rsun]',iunit) + call write_inopt(star%rcore,'rcore'//trim(c),'Radius of core softening',iunit) endif if ((star%isofteningopt == 2) .or. (star%isofteningopt == 3)) then - call write_inopt(in_solarm(star%mcore),'mcore'//trim(c),& - 'Mass of point mass stellar core [Msun]',iunit) + call write_inopt(star%mcore,'mcore'//trim(c),'Mass of point mass stellar core',iunit) endif elseif (star%isoftcore == 2) then - call write_inopt(in_solarr(star%rcore),'rcore'//trim(c),& - 'Radius of core softening [Rsun]',iunit) - call write_inopt(in_solarm(star%mcore),'mcore'//trim(c),& - 'Initial guess for mass of sink particle stellar core [Msun]',iunit) + call write_inopt(star%rcore,'rcore'//trim(c),'Radius of core softening',iunit) + call write_inopt(star%mcore,'mcore'//trim(c),& + 'Initial guess for mass of sink particle stellar core',iunit) endif - call write_inopt(in_solarl(star%lcore),'lcore'//trim(c),& - 'Luminosity of point mass stellar core [Lsun]',iunit) + call write_inopt(star%lcore,'lcore'//trim(c),& + 'Luminosity of point mass stellar core',iunit) else call write_inopt(star%isinkcore,'isinkcore'//trim(c),& 'Add a sink particle stellar core',iunit) if (star%isinkcore) then - call write_inopt(in_solarm(star%mcore),'mcore'//trim(c),& + call write_inopt(star%mcore,'mcore'//trim(c),& 'Mass of sink particle stellar core',iunit) - call write_inopt(in_solarr(star%hsoft),'hsoft'//trim(c),& - 'Softening length of sink particle stellar core [Rsun]',iunit) + call write_inopt(star%hsoft,'hsoft'//trim(c),& + 'Softening length of sink particle stellar core',iunit) endif - call write_inopt(in_solarl(star%lcore),'lcore'//trim(c),& - 'Luminosity of sink core particle [Lsun]',iunit) + call write_inopt(star%lcore,'lcore'//trim(c),'Luminosity of sink core particle',iunit) endif case (ievrard) call write_inopt(star%ui_coef,'ui_coef'//trim(c),& @@ -762,9 +795,16 @@ subroutine write_options_star(star,iunit,label) call write_inopt(star%hacc,'hacc'//trim(c),'accretion radius for sink'//trim(c),iunit) end select + if (need_polyk(star%iprofile)) call write_inopt(star%polyk,'polyk'//trim(c),'polytropic constant (cs^2 if isothermal)',iunit) + + ! options for setting initial thermal energy (e.g. if degenerate matter eos) + select case(ieos) + case(15) + call write_inopt(star%initialtemp,'initialtemp'//trim(c),'initial temperature of star (e.g. if degenerate matter eos)',iunit) + end select + if (star%iprofile > 0 .and. (len_trim(c)==0 .or. c(1:1)=='1')) then call write_inopt(star%np,'np'//trim(c),'number of particles',iunit) - !call write_inopt(use_exactN,'use_exactN','find closest particle number to np',iunit) endif end subroutine write_options_star @@ -774,21 +814,18 @@ end subroutine write_options_star ! read setupfile options needed for a star !+ !----------------------------------------------------------------------- -subroutine read_options_star(star,need_iso,ieos,polyk,db,nerr,label) +subroutine read_options_star(star,need_iso,ieos,db,nerr,label) use infile_utils, only:inopts,read_inopt use setstar_utils, only:need_inputprofile,need_rstar,nprofile_opts - use units, only:umass,udist,unit_luminosity - use physcon, only:solarm,solarr,solarl type(star_t), intent(out) :: star type(inopts), allocatable, intent(inout) :: db(:) integer, intent(out) :: need_iso integer, intent(inout) :: ieos - real, intent(inout) :: polyk integer, intent(inout) :: nerr character(len=*), intent(in), optional :: label character(len=10) :: c - real :: mcore_msun,rcore_rsun,lcore_lsun,mstar_msun,rstar_rsun,hsoft_rsun integer :: ierr + real :: rstar,mstar,rcore,mcore,hsoft,lcore,hacc ! set defaults call set_defaults_star(star) @@ -800,7 +837,7 @@ subroutine read_options_star(star,need_iso,ieos,polyk,db,nerr,label) call read_inopt(star%iprofile,'iprofile'//trim(c),db,errcount=nerr,min=0,max=nprofile_opts) call set_defaults_given_profile(star%iprofile,star%input_profile,& - need_iso,ieos,star%mstar,polyk) + need_iso,ieos,star%m,star%polyk) if (need_inputprofile(star%iprofile)) then call read_inopt(star%input_profile,'input_profile'//trim(c),db,errcount=nerr) @@ -819,10 +856,8 @@ subroutine read_options_star(star,need_iso,ieos,polyk,db,nerr,label) if (star%isoftcore <= 0) then ! sink particle core without softening call read_inopt(star%isinkcore,'isinkcore'//trim(c),db,errcount=nerr) if (star%isinkcore) then - call read_inopt(mcore_msun,'mcore'//trim(c),db,errcount=nerr,min=0.,err=ierr) - if (ierr==0) star%mcore = mcore_msun*real(solarm/umass) - call read_inopt(hsoft_rsun,'hsoft'//trim(c),db,errcount=nerr,min=0.,err=ierr) - if (ierr==0) star%hsoft = hsoft_rsun*real(solarr/udist) + call read_inopt(star%mcore,'mcore'//trim(c),db,errcount=nerr,err=ierr) + call read_inopt(star%hsoft,'hsoft'//trim(c),db,errcount=nerr,err=ierr) endif else star%isinkcore = .true. @@ -834,24 +869,29 @@ subroutine read_options_star(star,need_iso,ieos,polyk,db,nerr,label) endif if ((star%isofteningopt==1) .or. (star%isofteningopt==3)) then - call read_inopt(rcore_rsun,'rcore'//trim(c),db,errcount=nerr,min=0.,err=ierr) - if (ierr==0) star%rcore = rcore_rsun*real(solarr/udist) + call read_inopt(star%rcore,'rcore'//trim(c),db,errcount=nerr,err=ierr) endif if ((star%isofteningopt==2) .or. (star%isofteningopt==3) & .or. (star%isoftcore==2)) then - call read_inopt(mcore_msun,'mcore'//trim(c),db,errcount=nerr,min=0.,err=ierr) - if (ierr==0) star%mcore = mcore_msun*real(solarm/umass) + call read_inopt(star%mcore,'mcore'//trim(c),db,errcount=nerr,err=ierr) endif endif if (star%isinkcore) then - call read_inopt(lcore_lsun,'lcore'//trim(c),db,errcount=nerr,min=0.,err=ierr) - if (ierr==0) star%lcore = lcore_lsun*real(solarl/unit_luminosity) + call read_inopt(star%lcore,'lcore'//trim(c),db,errcount=nerr,err=ierr) endif case(ievrard) - call read_inopt(star%ui_coef,'ui_coef'//trim(c),db,errcount=nerr,min=0.) + call read_inopt(star%ui_coef,'ui_coef'//trim(c),db,errcount=nerr) case(:0) - call read_inopt(star%hacc,'hacc'//trim(c),db,errcount=nerr,min=0.) + call read_inopt(star%hacc,'hacc'//trim(c),db,errcount=nerr) + end select + + if (need_polyk(star%iprofile)) call read_inopt(star%polyk,'polyk'//trim(c),db,errcount=nerr) + + ! options for setting initial thermal energy (e.g. if degenerate matter eos) + select case(ieos) + case(15) + call read_inopt(star%initialtemp,'initialtemp'//trim(c),db,errcount=nerr,min=0.,max=1e12) end select ! star properties @@ -859,15 +899,16 @@ subroutine read_options_star(star,need_iso,ieos,polyk,db,nerr,label) if (need_inputprofile(star%iprofile)) then call read_inopt(star%input_profile,'input_profile'//trim(c),db,errcount=nerr) else - call read_inopt(mstar_msun,'Mstar'//trim(c),db,errcount=nerr,min=0.,err=ierr) - if (ierr==0) star%mstar = mstar_msun*real(solarm/umass) + call read_inopt(star%m,'Mstar'//trim(c),db,errcount=nerr,err=ierr) if (need_rstar(star%iprofile)) then - call read_inopt(rstar_rsun,'Rstar'//trim(c),db,errcount=nerr,min=0.,err=ierr) - if (ierr==0) star%rstar = rstar_rsun*real(solarr/udist) + call read_inopt(star%r,'Rstar'//trim(c),db,errcount=nerr,err=ierr) endif endif endif + ! perform a unit conversion, just to check that there are no errors parsing the .setup file + call get_star_properties_in_code_units(star,rstar,mstar,rcore,mcore,hsoft,lcore,hacc,nerr) + end subroutine read_options_star !----------------------------------------------------------------------- @@ -875,11 +916,11 @@ end subroutine read_options_star ! write_options routine that writes options for multiple stars !+ !----------------------------------------------------------------------- -subroutine write_options_stars(star,relax,iunit,nstar) +subroutine write_options_stars(star,relax,ieos,iunit,nstar) use relaxstar, only:write_options_relax use infile_utils, only:write_inopt type(star_t), intent(in) :: star(:) - integer, intent(in) :: iunit + integer, intent(in) :: ieos,iunit logical, intent(in) :: relax integer, intent(in), optional :: nstar integer :: i,nstars @@ -894,9 +935,15 @@ subroutine write_options_stars(star,relax,iunit,nstar) ! write options for each star do i=1,nstars - call write_options_star(star(i),iunit,label=achar(i+48)) + call write_options_star(star(i),iunit,ieos,label=achar(i+48)) enddo + ! write equation of state options if any stars made of gas + if (nstars > 0) then + if (any(star(1:nstars)%iprofile > 0)) & + call write_options_stars_eos(star(1:nstars),ieos,iunit) + endif + ! write relaxation options if any stars are made of gas if (nstars > 0) then if (any(star(1:nstars)%iprofile > 0)) then @@ -913,14 +960,13 @@ end subroutine write_options_stars ! read_options routine that reads options for multiple stars !+ !----------------------------------------------------------------------- -subroutine read_options_stars(star,need_iso,ieos,polyk,relax,db,nerr,nstar) +subroutine read_options_stars(star,need_iso,ieos,relax,db,nerr,nstar) use relaxstar, only:read_options_relax use infile_utils, only:inopts,read_inopt type(star_t), intent(out) :: star(:) type(inopts), allocatable, intent(inout) :: db(:) integer, intent(out) :: need_iso integer, intent(inout) :: ieos - real, intent(inout) :: polyk logical, intent(out) :: relax integer, intent(inout) :: nerr integer, intent(out), optional :: nstar @@ -928,15 +974,20 @@ subroutine read_options_stars(star,need_iso,ieos,polyk,relax,db,nerr,nstar) ! optionally ask for number of stars if (present(nstar)) then - call read_inopt(nstar,'nstars',db,nerr,min=0,max=size(star)) + call read_inopt(nstar,'nstars',db,errcount=nerr,min=0,max=size(star)) nstars = nstar else nstars = size(star) endif + ! write equation of state options if any stars made of gas + if (nstars > 0) then + if (any(star(1:nstars)%iprofile > 0)) call read_options_stars_eos(star,ieos,db,nerr) + endif + ! read options for each star do i=1,nstars - call read_options_star(star(i),need_iso,ieos,polyk,db,nerr,label=achar(i+48)) + call read_options_star(star(i),need_iso,ieos,db,nerr,label=achar(i+48)) enddo ! read relaxation options if any stars are made of gas @@ -949,4 +1000,73 @@ subroutine read_options_stars(star,need_iso,ieos,polyk,relax,db,nerr,nstar) end subroutine read_options_stars +!----------------------------------------------------------------------- +!+ +! write equation of state options needed to setup stars +!+ +!----------------------------------------------------------------------- +subroutine write_options_stars_eos(star,ieos,iunit) + use eos, only:use_var_comp,X_in,Z_in,irecomb,gmw,gamma + use infile_utils, only:write_inopt + integer, intent(in) :: ieos,iunit + type(star_t), intent(in) :: star(:) + + write(iunit,"(/,a)") '# equation of state used to set the thermal energy profile' + call write_inopt(ieos,'ieos','1=isothermal,2=adiabatic,10=MESA,12=idealplusrad',iunit) + + if (any(star(:)%iprofile==imesa)) then + call write_inopt(use_var_comp,'use_var_comp','Use variable composition (X, Z, mu)',iunit) + endif + + select case(ieos) + case(9) + write(iunit,"(/,a)") '# Piecewise Polytrope default options' + call write_inopt(EOSopt,'EOSopt','EOS: 1=APR3,2=SLy,3=MS1,4=ENG (from Read et al 2009)',iunit) + case(2,12) + call write_inopt(gamma,'gamma','Adiabatic index',iunit) + if (any(need_mu(star(:)%isoftcore)) .and. (.not. use_var_comp)) call write_inopt(gmw,'mu','mean molecular weight',iunit) + case(10,20) + if (ieos==20) call write_inopt(irecomb,'irecomb','Species to include in recombination (0: H2+H+He, 1:H+He, 2:He',iunit) + if ( (.not. use_var_comp) .and. any(need_mu(star(:)%isoftcore))) then + call write_inopt(X_in,'X','hydrogen mass fraction',iunit) + call write_inopt(Z_in,'Z','metallicity',iunit) + endif + end select + +end subroutine write_options_stars_eos + +!----------------------------------------------------------------------- +!+ +! read equation of state options needed to setup stars +!+ +!----------------------------------------------------------------------- +subroutine read_options_stars_eos(star,ieos,db,nerr) + use eos, only:use_var_comp,X_in,Z_in,irecomb,gamma,gmw + use infile_utils, only:inopts,read_inopt + type(star_t), intent(out) :: star(:) + type(inopts), allocatable, intent(inout) :: db(:) + integer, intent(inout) :: ieos + integer, intent(inout) :: nerr + + ! equation of state + call read_inopt(ieos,'ieos',db,errcount=nerr) + if (any(star(:)%iprofile==imesa)) call read_inopt(use_var_comp,'use_var_comp',db,errcount=nerr) + + select case(ieos) + case(9) + call read_inopt(EOSopt,'EOSopt',db,min=0,max=4,errcount=nerr) + case(2,12) + call read_inopt(gamma,'gamma',db,min=1.,max=7.,errcount=nerr) + if ( (.not. use_var_comp) .and. any(need_mu(star(:)%isoftcore))) call read_inopt(gmw,'mu',db,min=0.,errcount=nerr) + case(10,20) + if (ieos==20) call read_inopt(irecomb,'irecomb',db,errcount=nerr) + ! if softening stellar core, composition is automatically determined at R/2 + if ( (.not. use_var_comp) .and. any(need_mu(star(:)%isoftcore))) then + call read_inopt(X_in,'X',db,min=0.,max=1.,errcount=nerr) + call read_inopt(Z_in,'Z',db,min=0.,max=1.,errcount=nerr) + endif + end select + +end subroutine read_options_stars_eos + end module setstar diff --git a/src/setup/set_star_utils.f90 b/src/setup/set_star_utils.f90 index a82a1be51..449c4eecc 100644 --- a/src/setup/set_star_utils.f90 +++ b/src/setup/set_star_utils.f90 @@ -49,7 +49,7 @@ module setstar_utils public :: set_star_thermalenergy public :: set_stellar_core public :: write_kepler_comp - public :: need_inputprofile,need_polyk,need_rstar + public :: need_inputprofile,need_polyk,need_rstar,need_mu public :: get_mass_coord private @@ -199,7 +199,7 @@ end function need_inputprofile ! polytropic constant !+ !------------------------------------------------------------------------------- -logical function need_polyk(iprofile) +logical elemental function need_polyk(iprofile) integer, intent(in) :: iprofile select case(iprofile) @@ -211,6 +211,18 @@ logical function need_polyk(iprofile) end function need_polyk +!------------------------------------------------------------------------------- +!+ +! query function for whether mean molecular weight is needed +!+ +!------------------------------------------------------------------------------- +logical elemental function need_mu(isoftcore) + integer, intent(in) :: isoftcore + + need_mu = (isoftcore <= 0) + +end function need_mu + !------------------------------------------------------------------------------- !+ ! query function for whether a particular profile needs to read the @@ -359,7 +371,7 @@ subroutine get_mass_coord(i1,npart,xyzh,mass_enclosed_r) allocate(mass_enclosed_r(npart-i1),iorder(npart-i1)) ! sort particles by radius - call sort_by_radius(npart-i1,xyzh(1:3,i1+1:npart),iorder) + call sort_by_radius(npart-i1,xyzh(:,i1+1:npart),iorder) ! calculate cumulative mass massri = 0. diff --git a/src/setup/set_units.f90 b/src/setup/set_units.f90 index 5c6de9e7e..218eedae7 100644 --- a/src/setup/set_units.f90 +++ b/src/setup/set_units.f90 @@ -106,8 +106,8 @@ subroutine read_options_units(db,umass,udist,nerr,gr) if (present(gr)) nogr = .not.gr ! units - call read_inopt(mass_unit,'mass_unit',db,errcount=nerr) - if (nogr) call read_inopt(dist_unit,'dist_unit',db,errcount=nerr) + call read_inopt(mass_unit,'mass_unit',db,errcount=nerr,default=trim(mass_unit)) + if (nogr) call read_inopt(dist_unit,'dist_unit',db,errcount=nerr,default=trim(dist_unit)) ! ! parse units diff --git a/src/setup/setup_binary.f90 b/src/setup/setup_binary.f90 index f66acfae9..1291b6805 100644 --- a/src/setup/setup_binary.f90 +++ b/src/setup/setup_binary.f90 @@ -54,6 +54,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,& use setunits, only:mass_unit,dist_unit use physcon, only:deg_to_rad use kernel, only:hfact_default + use units, only:in_code_units integer, intent(in) :: id integer, intent(inout) :: npart integer, intent(out) :: npartoftype(:) @@ -67,6 +68,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,& integer :: ierr,nstar,nptmass_in,iextern_prev logical :: iexist,write_profile,use_var_comp,add_spin real :: xyzmh_ptmass_in(nsinkproperties,2),vxyz_ptmass_in(3,2),angle + real :: m1,m2,hacc1,hacc2 logical, parameter :: set_oblateness = .false. ! !--general parameters @@ -75,7 +77,8 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,& mass_unit = 'solarm' time = 0. polyk = 0. - gamma = 1. + gamma = 5./3. + ieos = 2 hfact = hfact_default ! !--space available for injected gas particles @@ -93,17 +96,18 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,& call set_defaults_orbit(orbit) relax = .true. corotate = .false. - ieos = 2 + use_var_comp = .false. + write_profile = .false. if (id==master) print "(/,65('-'),1(/,a),/,65('-'),/)",& ' Welcome to the Ultimate Binary Setup' filename = trim(fileprefix)//'.setup' inquire(file=filename,exist=iexist) - if (iexist) call read_setupfile(filename,ieos,polyk,ierr) + if (iexist) call read_setupfile(filename,ieos,ierr) if (.not. iexist .or. ierr /= 0) then if (id==master) then - call write_setupfile(filename) + call write_setupfile(filename,ieos) print*,' Edit '//trim(filename)//' and rerun phantomsetup' endif stop @@ -111,24 +115,27 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,& ! !--setup and relax stars as needed ! - use_var_comp = .false. - write_profile = .false. iextern_prev = iexternalforce iexternalforce = 0 - gamma = 5./3. call set_stars(id,master,nstar,star,xyzh,vxyzu,eos_vars,rad,npart,npartoftype,& - massoftype,hfact,xyzmh_ptmass,vxyz_ptmass,nptmass,ieos,polyk,gamma,& + massoftype,hfact,xyzmh_ptmass,vxyz_ptmass,nptmass,ieos,gamma,& X_in,Z_in,relax,use_var_comp,write_profile,& rhozero,npart_total,i_belong,ierr) nptmass_in = 0 + ! convert mass and accretion radii to code units + m1 = in_code_units(star(1)%m,ierr) + m2 = in_code_units(star(2)%m,ierr) + hacc1 = in_code_units(star(1)%hacc,ierr) + hacc2 = in_code_units(star(2)%hacc,ierr) + if (iextern_prev==iext_corotate) then - call set_orbit(orbit,star(1)%mstar,star(2)%mstar,star(1)%hacc,star(2)%hacc,& - xyzmh_ptmass_in,vxyz_ptmass_in,nptmass_in,(id==master),ierr,omega_corotate) + call set_orbit(orbit,m1,m2,hacc1,hacc2,xyzmh_ptmass_in,vxyz_ptmass_in,& + nptmass_in,(id==master),ierr,omega_corotate) add_spin = .false. else - call set_orbit(orbit,star(1)%mstar,star(2)%mstar,star(1)%hacc,star(2)%hacc,& - xyzmh_ptmass_in,vxyz_ptmass_in,nptmass_in,(id==master),ierr) + call set_orbit(orbit,m1,m2,hacc1,hacc2,xyzmh_ptmass_in,vxyz_ptmass_in,& + nptmass_in,(id==master),ierr) add_spin = corotate endif if (ierr /= 0) call fatal ('setup_binary','error in call to set_orbit') @@ -167,12 +174,13 @@ end subroutine setpart ! write options to .setup file !+ !---------------------------------------------------------------- -subroutine write_setupfile(filename) +subroutine write_setupfile(filename,ieos) use infile_utils, only:write_inopt use setstar, only:write_options_stars use setorbit, only:write_options_orbit use setunits, only:write_options_units character(len=*), intent(in) :: filename + integer, intent(in) :: ieos integer :: iunit print "(a)",' writing setup options file '//trim(filename) @@ -180,7 +188,7 @@ subroutine write_setupfile(filename) write(iunit,"(a)") '# input file for binary setup routines' call write_options_units(iunit,gr) - call write_options_stars(star,relax,iunit) + call write_options_stars(star,relax,ieos,iunit) call write_inopt(corotate,'corotate','set stars in corotation',iunit) call write_options_orbit(orbit,iunit) close(iunit) @@ -192,7 +200,7 @@ end subroutine write_setupfile ! read options from .setup file !+ !---------------------------------------------------------------- -subroutine read_setupfile(filename,ieos,polyk,ierr) +subroutine read_setupfile(filename,ieos,ierr) use infile_utils, only:open_db_from_file,inopts,read_inopt,close_db use io, only:error,fatal use setstar, only:read_options_stars @@ -200,7 +208,6 @@ subroutine read_setupfile(filename,ieos,polyk,ierr) use setunits, only:read_options_and_set_units character(len=*), intent(in) :: filename integer, intent(inout) :: ieos - real, intent(inout) :: polyk integer, intent(out) :: ierr integer, parameter :: iunit = 21 integer :: nerr,need_iso @@ -210,7 +217,7 @@ subroutine read_setupfile(filename,ieos,polyk,ierr) ierr = 0 call open_db_from_file(db,filename,iunit,ierr) call read_options_and_set_units(db,nerr,gr) - call read_options_stars(star,need_iso,ieos,polyk,relax,db,nerr) + call read_options_stars(star,need_iso,ieos,relax,db,nerr) if (need_iso==1) call fatal('setup_binary','incompatible setup for eos') call read_inopt(corotate,'corotate',db,errcount=nerr) call read_options_orbit(orbit,db,nerr) diff --git a/src/setup/setup_star.f90 b/src/setup/setup_star.f90 index 26c310081..1a11292f3 100644 --- a/src/setup/setup_star.f90 +++ b/src/setup/setup_star.f90 @@ -50,7 +50,7 @@ module setup real :: maxvxyzu logical :: iexist logical :: relax_star_in_setup,write_rho_to_file - type(star_t) :: star + type(star_t) :: star(1) public :: setpart private @@ -125,7 +125,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, setupfile = trim(fileprefix)//'.setup' inquire(file=setupfile,exist=setexists) if (setexists) then - call read_setupfile(setupfile,gamma,polyk,need_iso,ierr) + call read_setupfile(setupfile,gamma,need_iso,ierr) if (ierr /= 0) then if (id==master) call write_setupfile(setupfile,gamma,polyk) stop 'please rerun phantomsetup with revised .setup file' @@ -133,7 +133,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, !--Prompt to get inputs and write to file elseif (id==master) then print "(a,/)",trim(setupfile)//' not found: using interactive setup' - call setup_interactive(polyk,gamma,iexist,id,master,ierr) + call setup_interactive(gamma,iexist,id,master,ierr) call write_setupfile(setupfile,gamma,polyk) stop 'please check and edit .setup file and rerun phantomsetup' endif @@ -161,8 +161,8 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, npart = 0 nptmass = 0 vxyzu = 0.0 - call set_star(id,master,star,xyzh,vxyzu,eos_vars,rad,npart,npartoftype,& - massoftype,hfact,xyzmh_ptmass,vxyz_ptmass,nptmass,ieos,polyk,gamma,& + call set_star(id,master,star(1),xyzh,vxyzu,eos_vars,rad,npart,npartoftype,& + massoftype,hfact,xyzmh_ptmass,vxyz_ptmass,nptmass,ieos,gamma,& X_in,Z_in,relax_star_in_setup,use_var_comp,write_rho_to_file,& rhozero,npart_total,i_belong,ierr) ! @@ -173,7 +173,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, ! ! override some default settings in the .in file for some cases ! - select case(star%iprofile) + select case(star(1)%iprofile) case(ibpwpoly) ! piecewise polytrope calc_erot = .true. case(ievrard) ! Evrard Collapse @@ -191,14 +191,14 @@ end subroutine setpart ! Ask questions of the user to determine which setup to use !+ !----------------------------------------------------------------------- -subroutine setup_interactive(polyk,gamma,iexist,id,master,ierr) +subroutine setup_interactive(gamma,iexist,id,master,ierr) use prompting, only:prompt use units, only:select_unit use eos, only:X_in,Z_in,gmw use eos_gasradrec, only:irecomb use setstar, only:set_star_interactive use setunits, only:set_units_interactive - real, intent(out) :: polyk,gamma + real, intent(out) :: gamma logical, intent(in) :: iexist integer, intent(in) :: id,master integer, intent(out) :: ierr @@ -209,35 +209,29 @@ subroutine setup_interactive(polyk,gamma,iexist,id,master,ierr) call set_units_interactive(gr) ! star - call set_star_interactive(id,master,star,need_iso,use_var_comp,ieos,polyk) - - ! equation of state - call prompt('Enter the desired EoS (1=isothermal,2=adiabatic,10=MESA,12=idealplusrad)',ieos) - select case(ieos) - case(15) ! Helmholtz - call prompt('Enter temperature',star%initialtemp,1.0e3,1.0e11) - case(9) - write(*,'(a)') 'EOS options: 1=APR3,2=SLy,3=MS1,4=ENG (from Read et al 2009)' - call prompt('Enter equation of state type',EOSopt,1,4) - case(2) - call prompt('Enter gamma (adiabatic index)',gamma,1.,7.) - case(20) - call prompt('Enter irecomb (0: H2+H+He, 1:H+He, 2:He)',irecomb,0) - end select - - if (need_polyk(star%iprofile)) then - call prompt('Enter polytropic constant (cs^2 if isothermal)',polyk,0.) - endif - - if ((.not. use_var_comp) .and. (star%isoftcore<=0)) then - if ( (ieos==12) .or. (ieos==2) ) call prompt('Enter mean molecular weight',gmw,0.) - if ( (ieos==10) .or. (ieos==20) ) then - call prompt('Enter hydrogen mass fraction (X)',X_in,0.,1.) - call prompt('Enter metals mass fraction (Z)',Z_in,0.,1.) - endif - endif - - call prompt('Relax star automatically during setup?',relax_star_in_setup) + !call set_star_interactive(id,master,star,need_iso,use_var_comp,ieos) + +! ! equation of state +! call prompt('Enter the desired EoS (1=isothermal,2=adiabatic,10=MESA,12=idealplusrad)',ieos) +! select case(ieos) +! case(9) +! write(*,'(a)') 'EOS options: 1=APR3,2=SLy,3=MS1,4=ENG (from Read et al 2009)' +! call prompt('Enter equation of state type',EOSopt,1,4) +! case(2) +! call prompt('Enter gamma (adiabatic index)',gamma,1.,7.) +! case(20) +! call prompt('Enter irecomb (0: H2+H+He, 1:H+He, 2:He)',irecomb,0) +! end select + +! if ((.not. use_var_comp) .and. (star%isoftcore<=0)) then +! if ( (ieos==12) .or. (ieos==2) ) call prompt('Enter mean molecular weight',gmw,0.) +! if ( (ieos==10) .or. (ieos==20) ) then +! call prompt('Enter hydrogen mass fraction (X)',X_in,0.,1.) +! call prompt('Enter metals mass fraction (Z)',Z_in,0.,1.) +! endif +! endif + +! call prompt('Relax star automatically during setup?',relax_star_in_setup) end subroutine setup_interactive @@ -250,9 +244,9 @@ subroutine write_setupfile(filename,gamma,polyk) use infile_utils, only:write_inopt use dim, only:tagline,use_apr use relaxstar, only:write_options_relax - use eos, only:X_in,Z_in,gmw + use eos, only:ieos use eos_gasradrec, only:irecomb - use setstar, only:write_options_star,need_polyk + use setstar, only:write_options_stars,need_polyk use setunits, only:write_options_units use apr, only:write_options_apr real, intent(in) :: gamma,polyk @@ -265,40 +259,38 @@ subroutine write_setupfile(filename,gamma,polyk) write(iunit,"(a)") '# input file for Phantom star setup' call write_options_units(iunit,gr) - call write_options_star(star,iunit) - - write(iunit,"(/,a)") '# equation of state' - call write_inopt(ieos,'ieos','1=isothermal,2=adiabatic,10=MESA,12=idealplusrad',iunit) - - if (star%iprofile==imesa) then - call write_inopt(use_var_comp,'use_var_comp','Use variable composition (X, Z, mu)',iunit) - endif - - select case(ieos) - case(15) ! Helmholtz - call write_inopt(star%initialtemp,'initialtemp','initial temperature of the star',iunit) - case(9) - write(iunit,"(/,a)") '# Piecewise Polytrope default options' - call write_inopt(EOSopt,'EOSopt','EOS: 1=APR3,2=SLy,3=MS1,4=ENG (from Read et al 2009)',iunit) - case(2) - call write_inopt(gamma,'gamma','Adiabatic index',iunit) - if ((star%isoftcore<=0) .and. (.not. use_var_comp)) call write_inopt(gmw,'mu','mean molecular weight',iunit) - case(10,20) - if (ieos==20) call write_inopt(irecomb,'irecomb','Species to include in recombination (0: H2+H+He, 1:H+He, 2:He',iunit) - if ( (.not. use_var_comp) .and. (star%isoftcore <= 0) ) then - call write_inopt(X_in,'X','hydrogen mass fraction',iunit) - call write_inopt(Z_in,'Z','metallicity',iunit) - endif - case(12) - call write_inopt(gamma,'gamma','Adiabatic index',iunit) - if ((star%isoftcore<=0) .and. (.not. use_var_comp)) call write_inopt(gmw,'mu','mean molecular weight',iunit) - end select - - if (need_polyk(star%iprofile)) call write_inopt(polyk,'polyk','polytropic constant (cs^2 if isothermal)',iunit) - - write(iunit,"(/,a)") '# relaxation options' - call write_inopt(relax_star_in_setup,'relax_star','relax star(s) automatically during setup',iunit) - if (relax_star_in_setup) call write_options_relax(iunit) + call write_options_stars(star,relax_star_in_setup,ieos,iunit) + + !write(iunit,"(/,a)") '# equation of state' + !call write_inopt(ieos,'ieos','1=isothermal,2=adiabatic,10=MESA,12=idealplusrad',iunit) + + !if (star%iprofile==imesa) then + ! call write_inopt(use_var_comp,'use_var_comp','Use variable composition (X, Z, mu)',iunit) + !endif + + !select case(ieos) + !case(15) ! Helmholtz + ! call write_inopt(star%initialtemp,'initialtemp','initial temperature of the star',iunit) + !case(9) + ! write(iunit,"(/,a)") '# Piecewise Polytrope default options' + ! call write_inopt(EOSopt,'EOSopt','EOS: 1=APR3,2=SLy,3=MS1,4=ENG (from Read et al 2009)',iunit) + !case(2) + ! call write_inopt(gamma,'gamma','Adiabatic index',iunit) + ! if ((star%isoftcore<=0) .and. (.not. use_var_comp)) call write_inopt(gmw,'mu','mean molecular weight',iunit) + !case(10,20) + ! if (ieos==20) call write_inopt(irecomb,'irecomb','Species to include in recombination (0: H2+H+He, 1:H+He, 2:He',iunit) + ! if ( (.not. use_var_comp) .and. (star%isoftcore <= 0) ) then + ! call write_inopt(X_in,'X','hydrogen mass fraction',iunit) + ! call write_inopt(Z_in,'Z','metallicity',iunit) + ! endif + !case(12) + ! call write_inopt(gamma,'gamma','Adiabatic index',iunit) + ! if ((star%isoftcore<=0) .and. (.not. use_var_comp)) call write_inopt(gmw,'mu','mean molecular weight',iunit) + !end select + + !write(iunit,"(/,a)") '# relaxation options' + !call write_inopt(relax_star_in_setup,'relax_star','relax star(s) automatically during setup',iunit) + !if (relax_star_in_setup) call write_options_relax(iunit) call write_inopt(write_rho_to_file,'write_rho_to_file','write density profile(s) to file',iunit) @@ -312,21 +304,21 @@ end subroutine write_setupfile ! Read setup parameters from input file !+ !----------------------------------------------------------------------- -subroutine read_setupfile(filename,gamma,polyk,need_iso,ierr) +subroutine read_setupfile(filename,gamma,need_iso,ierr) use infile_utils, only:open_db_from_file,inopts,close_db,read_inopt use io, only:error use units, only:select_unit use relaxstar, only:read_options_relax - use eos, only:X_in,Z_in,gmw + use eos, only:X_in,Z_in,gmw,ieos use eos_gasradrec, only:irecomb - use setstar, only:read_options_star + use setstar, only:read_options_stars use setunits, only:read_options_and_set_units use apr, only:apr_max_in,ref_dir,apr_type,apr_rad,apr_drad use dim, only:use_apr character(len=*), intent(in) :: filename integer, parameter :: lu = 21 integer, intent(out) :: need_iso,ierr - real, intent(out) :: gamma,polyk + real, intent(out) :: gamma integer :: nerr type(inopts), allocatable :: db(:) @@ -339,39 +331,39 @@ subroutine read_setupfile(filename,gamma,polyk,need_iso,ierr) call read_options_and_set_units(db,nerr,gr) ! star options - call read_options_star(star,need_iso,ieos,polyk,db,nerr) - - ! equation of state - call read_inopt(ieos,'ieos',db,errcount=nerr) - if (star%iprofile==imesa) call read_inopt(use_var_comp,'use_var_comp',db,errcount=nerr) - - select case(ieos) - case(15) ! Helmholtz - call read_inopt(star%initialtemp,'initialtemp',db,errcount=nerr) - case(9) - call read_inopt(EOSopt,'EOSopt',db,errcount=nerr) - case(2) - call read_inopt(gamma,'gamma',db,errcount=nerr) - if ( (.not. use_var_comp) .and. (star%isoftcore <= 0)) call read_inopt(gmw,'mu',db,errcount=nerr) - case(10,20) - if (ieos==20) call read_inopt(irecomb,'irecomb',db,errcount=nerr) - ! if softening stellar core, composition is automatically determined at R/2 - if ( (.not. use_var_comp) .and. (star%isoftcore <= 0)) then - call read_inopt(X_in,'X',db,errcount=nerr) - call read_inopt(Z_in,'Z',db,errcount=nerr) - endif - case(12) + call read_options_stars(star,need_iso,ieos,relax_star_in_setup,db,nerr) + + !! equation of state + !call read_inopt(ieos,'ieos',db,errcount=nerr) + !if (star%iprofile==imesa) call read_inopt(use_var_comp,'use_var_comp',db,errcount=nerr) + + !select case(ieos) + !case(15) ! Helmholtz + ! call read_inopt(star%initialtemp,'initialtemp',db,errcount=nerr) + !case(9) + ! call read_inopt(EOSopt,'EOSopt',db,errcount=nerr) + !case(2) + ! call read_inopt(gamma,'gamma',db,errcount=nerr) + ! if ( (.not. use_var_comp) .and. (star%isoftcore <= 0)) call read_inopt(gmw,'mu',db,errcount=nerr) + !case(10,20) + ! if (ieos==20) call read_inopt(irecomb,'irecomb',db,errcount=nerr) + ! ! if softening stellar core, composition is automatically determined at R/2 + ! if ( (.not. use_var_comp) .and. (star%isoftcore <= 0)) then + ! call read_inopt(X_in,'X',db,errcount=nerr) + ! call read_inopt(Z_in,'Z',db,errcount=nerr) + ! endif + !case(12) ! if softening stellar core, mu is automatically determined at R/2 - call read_inopt(gamma,'gamma',db,errcount=nerr) - if ( (.not. use_var_comp) .and. (star%isoftcore <= 0)) call read_inopt(gmw,'mu',db,errcount=nerr) - end select + ! call read_inopt(gamma,'gamma',db,errcount=nerr) + ! if ( (.not. use_var_comp) .and. (star%isoftcore <= 0)) call read_inopt(gmw,'mu',db,errcount=nerr) +! end select - if (need_polyk(star%iprofile)) call read_inopt(polyk,'polyk',db,errcount=nerr) + !if (need_polyk(star%iprofile)) call read_inopt(polyk,'polyk',db,errcount=nerr) ! relax star options - call read_inopt(relax_star_in_setup,'relax_star',db,errcount=nerr) - if (relax_star_in_setup) call read_options_relax(db,nerr) - if (nerr /= 0) ierr = ierr + 1 + !call read_inopt(relax_star_in_setup,'relax_star',db,errcount=nerr) + !if (relax_star_in_setup) call read_options_relax(db,nerr)! + !if (nerr /= 0) ierr = ierr + 1 ! option to write density profile to file call read_inopt(write_rho_to_file,'write_rho_to_file',db) From 7ed2c3f77096ff7eafd14b3740f471e382b68eb1 Mon Sep 17 00:00:00 2001 From: Yrisch Date: Mon, 2 Dec 2024 17:31:15 +0100 Subject: [PATCH 11/36] (substepping) reorganize substep and get_force routine for clarity --- src/main/evolve.F90 | 7 +- src/main/initial.F90 | 21 ++-- src/main/ptmass.F90 | 46 ++++----- src/main/substepping.F90 | 217 +++++++++++++++------------------------ 4 files changed, 107 insertions(+), 184 deletions(-) diff --git a/src/main/evolve.F90 b/src/main/evolve.F90 index 16d45adc2..b899b6d5a 100644 --- a/src/main/evolve.F90 +++ b/src/main/evolve.F90 @@ -333,12 +333,9 @@ subroutine evol(infile,logfile,evfile,dumpfile,flag) if (use_regnbody) then call group_identify(nptmass,n_group,n_ingroup,n_sing,xyzmh_ptmass,vxyz_ptmass,group_info,bin_info,nmatrix,& new_ptmass=.true.,dtext=dtextforce) - call get_force(nptmass,npart,0,1,time,dtextforce,xyzh,vxyzu,fext,xyzmh_ptmass,vxyz_ptmass,& - fxyz_ptmass,dsdt_ptmass,0.,0.,dummy,.false.,linklist_ptmass,bin_info,group_info=group_info) - else - call get_force(nptmass,npart,0,1,time,dtextforce,xyzh,vxyzu,fext,xyzmh_ptmass,vxyz_ptmass,& - fxyz_ptmass,dsdt_ptmass,0.,0.,dummy,.false.,linklist_ptmass,bin_info) endif + call get_force(nptmass,npart,0,1,time,dtextforce,xyzh,vxyzu,fext,xyzmh_ptmass,vxyz_ptmass,& + fxyz_ptmass,dsdt_ptmass,0.,0.,dummy,.false.,linklist_ptmass,bin_info,group_info) if (ipart_createseeds /= 0) ipart_createseeds = 0 ! reset pointer to zero if (ipart_createstars /= 0) ipart_createstars = 0 ! reset pointer to zero dummy = 0 diff --git a/src/main/initial.F90 b/src/main/initial.F90 index aaeb33633..67b480e15 100644 --- a/src/main/initial.F90 +++ b/src/main/initial.F90 @@ -529,14 +529,10 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread) if (use_regnbody) then call init_subgroup call group_identify(nptmass,n_group,n_ingroup,n_sing,xyzmh_ptmass,vxyz_ptmass,group_info,bin_info,nmatrix) - call get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,epot_sinksink,dtsinksink,& - iexternalforce,time,merge_ij,merge_n,dsdt_ptmass,& - group_info=group_info,bin_info=bin_info) - - else - call get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,epot_sinksink,dtsinksink,& - iexternalforce,time,merge_ij,merge_n,dsdt_ptmass) endif + call get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,epot_sinksink,dtsinksink,& + iexternalforce,time,merge_ij,merge_n,dsdt_ptmass,& + group_info,bin_info) dtsinksink = C_force*dtsinksink if (id==master) write(iprint,*) 'dt(sink-sink) = ',dtsinksink dtextforce = min(dtextforce,dtsinksink) @@ -555,14 +551,9 @@ subroutine startrun(infile,logfile,evfile,dumpfile,noread) elseif (use_apr) then pmassi = aprmassoftype(igas,apr_level(i)) endif - if (use_regnbody) then - call get_accel_sink_gas(nptmass,xyzh(1,i),xyzh(2,i),xyzh(3,i),xyzh(4,i),xyzmh_ptmass, & - fext(1,i),fext(2,i),fext(3,i),poti,pmassi,fxyz_ptmass,& - dsdt_ptmass,fonrmax,dtphi2,bin_info=bin_info) - else - call get_accel_sink_gas(nptmass,xyzh(1,i),xyzh(2,i),xyzh(3,i),xyzh(4,i),xyzmh_ptmass, & - fext(1,i),fext(2,i),fext(3,i),poti,pmassi,fxyz_ptmass,dsdt_ptmass,fonrmax,dtphi2) - endif + call get_accel_sink_gas(nptmass,xyzh(1,i),xyzh(2,i),xyzh(3,i),xyzh(4,i),xyzmh_ptmass, & + fext(1,i),fext(2,i),fext(3,i),poti,pmassi,fxyz_ptmass,& + dsdt_ptmass,fonrmax,dtphi2,bin_info) dtsinkgas = min(dtsinkgas,C_force*1./sqrt(fonrmax),C_force*sqrt(dtphi2)) endif enddo diff --git a/src/main/ptmass.F90 b/src/main/ptmass.F90 index b0a148042..4f4c9f995 100644 --- a/src/main/ptmass.F90 +++ b/src/main/ptmass.F90 @@ -150,8 +150,7 @@ module ptmass !+ !---------------------------------------------------------------- subroutine get_accel_sink_gas(nptmass,xi,yi,zi,hi,xyzmh_ptmass,fxi,fyi,fzi,phi, & - pmassi,fxyz_ptmass,dsdt_ptmass,fonrmax, & - dtphi2,extrapfac,fsink_old,bin_info) + pmassi,fxyz_ptmass,dsdt_ptmass,fonrmax,dtphi2,bin_info,extrapfac,fsink_old) #ifdef FINVSQRT use fastmath, only:finvsqrt #endif @@ -165,15 +164,15 @@ subroutine get_accel_sink_gas(nptmass,xi,yi,zi,hi,xyzmh_ptmass,fxi,fyi,fzi,phi, real, intent(in) :: xyzmh_ptmass(nsinkproperties,nptmass) real, optional, intent(in) :: pmassi,extrapfac real, optional, intent(inout) :: fxyz_ptmass(4,nptmass),dsdt_ptmass(3,nptmass) - real, optional, intent(inout) :: bin_info(6,nptmass) real, optional, intent(in) :: fsink_old(4,nptmass) real, optional, intent(out) :: fonrmax,dtphi2 + real, optional, intent(inout) :: bin_info(6,nptmass) real :: ftmpxi,ftmpyi,ftmpzi real :: dx,dy,dz,rr2,ddr,dr3,f1,f2,pmassj,J2,shat(3),Rsink real :: hsoft,hsoft1,hsoft21,q2i,qi,psoft,fsoft real :: fxj,fyj,fzj,dsx,dsy,dsz,fac,r integer :: j - logical :: tofrom,extrap,kappa + logical :: tofrom,extrap ! ! Determine if acceleration is from/to gas, or to gas ! @@ -191,12 +190,6 @@ subroutine get_accel_sink_gas(nptmass,xi,yi,zi,hi,xyzmh_ptmass,fxi,fyi,fzi,phi, extrap = .false. endif - if (present(bin_info)) then - kappa = .true. - else - kappa = .false. - endif - ftmpxi = 0. ! use temporary summation variable ftmpyi = 0. ! (better for round-off, plus we need this bit of @@ -303,7 +296,7 @@ subroutine get_accel_sink_gas(nptmass,xi,yi,zi,hi,xyzmh_ptmass,fxi,fyi,fzi,phi, ! timestep is sqrt(separation/force) fonrmax = max(f1,f2,fonrmax) - if (kappa) then + if (use_regnbody) then if (abs(bin_info(isemi,j))>tiny(f2)) then bin_info(ipert,j) = bin_info(ipert,j) + f2 endif @@ -338,8 +331,8 @@ end subroutine get_accel_sink_gas !+ !---------------------------------------------------------------- subroutine get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,phitot,dtsinksink,& - iexternalforce,ti,merge_ij,merge_n,dsdt_ptmass,extrapfac,fsink_old,& - group_info,bin_info) + iexternalforce,ti,merge_ij,merge_n,dsdt_ptmass,group_info,bin_info,& + extrapfac,fsink_old) #ifdef FINVSQRT use fastmath, only:finvsqrt #endif @@ -356,10 +349,10 @@ subroutine get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,phitot,dtsinksin real, intent(in) :: ti integer, intent(out) :: merge_ij(:),merge_n real, intent(out) :: dsdt_ptmass(3,nptmass) + integer, optional, intent(in) :: group_info(4,nptmass) + real, optional, intent(out) :: bin_info(6,nptmass) real, optional, intent(in) :: extrapfac real, optional, intent(in) :: fsink_old(4,nptmass) - real, optional, intent(out) :: bin_info(6,nptmass) - integer, optional, intent(in) :: group_info(4,nptmass) real :: xi,yi,zi,pmassi,pmassj,hacci,haccj,fxi,fyi,fzi,phii real :: ddr,dx,dy,dz,rr2,rr2j,dr3,f1,f2 real :: hsoft1,hsoft21,q2i,qi,psoft,fsoft @@ -368,7 +361,7 @@ subroutine get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,phitot,dtsinksin real :: J2i,rsinki,shati(3) real :: J2j,rsinkj,shatj(3) integer :: k,l,i,j,gidi,gidj,compi - logical :: extrap,subsys + logical :: extrap dtsinksink = huge(dtsinksink) fxyz_ptmass(:,:) = 0. @@ -376,6 +369,8 @@ subroutine get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,phitot,dtsinksin phitot = 0. merge_n = 0 merge_ij = 0 + gidi = 0 + gidj = 0 if (nptmass <= 0) return ! check if it is a force computed using Omelyan extrapolation method for FSI if (present(extrapfac) .and. present(fsink_old)) then @@ -384,11 +379,7 @@ subroutine get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,phitot,dtsinksin extrap = .false. endif - if (present(group_info) .and. present(bin_info)) then - subsys = .true. - else - subsys = .false. - endif + ! !--get self-contribution to the potential if sink-sink softening is used ! @@ -408,18 +399,19 @@ subroutine get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,phitot,dtsinksin !$omp shared(nptmass,xyzmh_ptmass,fxyz_ptmass,merge_ij,r_merge2,dsdt_ptmass) & !$omp shared(iexternalforce,ti,h_soft_sinksink,potensoft0,hsoft1,hsoft21) & !$omp shared(extrapfac,extrap,fsink_old,h_acc,icreate_sinks) & - !$omp shared(group_info,bin_info,subsys) & + !$omp shared(group_info,bin_info,use_regnbody) & !$omp private(i,j,xi,yi,zi,pmassi,pmassj,hacci,haccj) & - !$omp private(gidi,gidj,compi,pert_out) & + !$omp private(compi,pert_out) & !$omp private(dx,dy,dz,rr2,rr2j,ddr,dr3,f1,f2) & !$omp private(fxi,fyi,fzi,phii,dsx,dsy,dsz) & !$omp private(fextx,fexty,fextz,phiext) & !$omp private(q2i,qi,psoft,fsoft) & !$omp private(fterm,pterm,J2i,J2j,shati,shatj,rsinki,rsinkj) & + !$omp firstprivate(gidi,gidj)& !$omp reduction(min:dtsinksink) & !$omp reduction(+:phitot,merge_n) do k=1,nptmass - if (subsys) then + if (use_regnbody) then pert_out = 0. i = group_info(igarg,k) ! new id order when using group info gidi = group_info(igid,k) ! id of the group to identify which ptmasses are in the same group @@ -451,7 +443,7 @@ subroutine get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,phitot,dtsinksin dsy = 0. dsz = 0. do l=1,nptmass - if (subsys) then + if (use_regnbody) then j = group_info(igarg,l) gidj = group_info(igid,l) if (gidi==gidj) cycle @@ -551,13 +543,13 @@ subroutine get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,phitot,dtsinksin endif endif endif - if (subsys) then + if (use_regnbody) then if (compi /= i) pert_out = pert_out + f1 endif enddo phitot = phitot + 0.5*pmassi*phii ! total potential (G M_1 M_2/r) - if (subsys) bin_info(ipert,i) = pert_out + if (use_regnbody) bin_info(ipert,i) = pert_out ! !--apply external forces diff --git a/src/main/substepping.F90 b/src/main/substepping.F90 index f16a842d6..08f5d88a5 100644 --- a/src/main/substepping.F90 +++ b/src/main/substepping.F90 @@ -448,7 +448,7 @@ subroutine substep(npart,ntypes,nptmass,dtsph,dtextforce,time,xyzh,vxyzu,fext, & use io_summary, only:summary_variable,iosumextr,iosumextt use externalforces, only:is_velocity_dependent use ptmass, only:use_fourthorder,use_regnbody,ck,dk,ptmass_check_stars,icreate_sinks - use subgroup, only:group_identify,evolve_groups + use subgroup, only:group_identify integer, intent(in) :: npart,ntypes,nptmass integer, intent(inout) :: n_group,n_ingroup,n_sing integer, intent(inout) :: group_info(:,:) @@ -500,56 +500,33 @@ subroutine substep(npart,ntypes,nptmass,dtsph,dtextforce,time,xyzh,vxyzu,fext, & call kick(dk(1),dt,npart,nptmass,ntypes,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass, & fext,fxyz_ptmass,dsdt_ptmass,dptmass) - if (use_regnbody) then - call evolve_groups(n_group,nptmass,time_par,time_par+ck(1)*dt,group_info,bin_info, & - xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,gtgrad) + call drift(ck(1),dt,time_par,npart,nptmass,ntypes,xyzh,xyzmh_ptmass,& + vxyzu,vxyz_ptmass,fxyz_ptmass,gtgrad,n_group,n_ingroup,& + group_info,bin_info) - call drift(ck(1),dt,time_par,npart,nptmass,ntypes,xyzh,xyzmh_ptmass,vxyzu,vxyz_ptmass,n_ingroup,group_info) - - call get_force(nptmass,npart,nsubsteps,ntypes,time_par,dtextforce,xyzh,vxyzu,fext,xyzmh_ptmass, & - vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,dt,dk(2),force_count,extf_vdep_flag,linklist_ptmass,& - bin_info,group_info=group_info) - else - call drift(ck(1),dt,time_par,npart,nptmass,ntypes,xyzh,xyzmh_ptmass,vxyzu,vxyz_ptmass) - - call get_force(nptmass,npart,nsubsteps,ntypes,time_par,dtextforce,xyzh,vxyzu,fext,xyzmh_ptmass, & - vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,dt,dk(2),force_count,extf_vdep_flag,linklist_ptmass,& - bin_info,isionised=isionised) - endif + call get_force(nptmass,npart,nsubsteps,ntypes,time_par,dtextforce,xyzh,vxyzu,fext,xyzmh_ptmass, & + vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,dt,dk(2),force_count,extf_vdep_flag,linklist_ptmass,& + bin_info,group_info,isionised=isionised) if (use_fourthorder) then !! FSI 4th order scheme ! FSI extrapolation method (Omelyan 2006) - if (use_regnbody) then - call get_force(nptmass,npart,nsubsteps,ntypes,time_par,dtextforce,xyzh,vxyzu,fext,xyzmh_ptmass, & - vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,dt,dk(2),force_count,extf_vdep_flag,linklist_ptmass, & - bin_info,fsink_old,group_info=group_info) - - call kick(dk(2),dt,npart,nptmass,ntypes,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& - fext,fxyz_ptmass,dsdt_ptmass,dptmass) + call get_force(nptmass,npart,nsubsteps,ntypes,time_par,dtextforce,xyzh,vxyzu,fext,xyzmh_ptmass, & + vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,dt,dk(2),force_count,extf_vdep_flag,linklist_ptmass, & + bin_info,group_info,fsink_old) - call evolve_groups(n_group,nptmass,time_par,time_par+ck(2)*dt,group_info,bin_info, & - xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,gtgrad) + call kick(dk(2),dt,npart,nptmass,ntypes,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& + fext,fxyz_ptmass,dsdt_ptmass,dptmass) - call drift(ck(2),dt,time_par,npart,nptmass,ntypes,xyzh,xyzmh_ptmass,vxyzu,vxyz_ptmass,n_ingroup,group_info) + call drift(ck(2),dt,time_par,npart,nptmass,ntypes,xyzh,xyzmh_ptmass,& + vxyzu,vxyz_ptmass,fxyz_ptmass,gtgrad,n_group,n_ingroup,& + group_info,bin_info) - call get_force(nptmass,npart,nsubsteps,ntypes,time_par,dtextforce,xyzh,vxyzu,fext,xyzmh_ptmass, & - vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,dt,dk(3),force_count,extf_vdep_flag,linklist_ptmass, & - bin_info,group_info=group_info,isionised=isionised) - else - call get_force(nptmass,npart,nsubsteps,ntypes,time_par,dtextforce,xyzh,vxyzu,fext,xyzmh_ptmass, & - vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,dt,dk(2),force_count,extf_vdep_flag,linklist_ptmass,& - bin_info,fsink_old) - call kick(dk(2),dt,npart,nptmass,ntypes,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& - fext,fxyz_ptmass,dsdt_ptmass,dptmass) - - call drift(ck(2),dt,time_par,npart,nptmass,ntypes,xyzh,xyzmh_ptmass,vxyzu,vxyz_ptmass) + call get_force(nptmass,npart,nsubsteps,ntypes,time_par,dtextforce,xyzh,vxyzu,fext,xyzmh_ptmass, & + vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,dt,dk(3),force_count,extf_vdep_flag,linklist_ptmass, & + bin_info,group_info,isionised=isionised) - call get_force(nptmass,npart,nsubsteps,ntypes,time_par,dtextforce,xyzh,vxyzu,fext,xyzmh_ptmass, & - vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,dt,dk(3),force_count,extf_vdep_flag,linklist_ptmass,& - bin_info,isionised=isionised) - ! the last kick phase of the scheme will perform the accretion loop after velocity update - endif + ! the last kick phase of the scheme will perform the accretion loop after velocity update call kick(dk(3),dt,npart,nptmass,ntypes,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,fext, & fxyz_ptmass,dsdt_ptmass,dptmass,ibin_wake,nbinmax,timei, & @@ -560,11 +537,11 @@ subroutine substep(npart,ntypes,nptmass,dtsph,dtextforce,time,xyzh,vxyzu,fext, & dtext=dtextforce) call get_force(nptmass,npart,nsubsteps,ntypes,time_par,dtextforce,xyzh,vxyzu,fext,xyzmh_ptmass, & vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,dt,dk(3),force_count,extf_vdep_flag,linklist_ptmass, & - bin_info,group_info=group_info) + bin_info,group_info) elseif (accreted) then call get_force(nptmass,npart,nsubsteps,ntypes,time_par,dtextforce,xyzh,vxyzu,fext,xyzmh_ptmass, & vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,dt,dk(3),force_count,extf_vdep_flag,linklist_ptmass,& - bin_info) + bin_info,group_info) endif else !! standard leapfrog scheme ! the last kick phase of the scheme will perform the accretion loop after velocity update @@ -574,7 +551,7 @@ subroutine substep(npart,ntypes,nptmass,dtsph,dtextforce,time,xyzh,vxyzu,fext, & if (accreted) then call get_force(nptmass,npart,nsubsteps,ntypes,time_par,dtextforce,xyzh,vxyzu,fext,xyzmh_ptmass, & vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,dt,dk(2),force_count,extf_vdep_flag,linklist_ptmass,& - bin_info) + bin_info,group_info) endif endif @@ -610,18 +587,22 @@ end subroutine substep !+ !---------------------------------------------------------------- -subroutine drift(cki,dt,time_par,npart,nptmass,ntypes,xyzh,xyzmh_ptmass,vxyzu,vxyz_ptmass,n_ingroup,group_info) +subroutine drift(cki,dt,time_par,npart,nptmass,ntypes,xyzh,xyzmh_ptmass,vxyzu, & + vxyz_ptmass,fxyz_ptmass,gtgrad,n_group,n_ingroup,group_info, & + bin_info) use part, only: isdead_or_accreted,ispinx,ispiny,ispinz,igarg - use ptmass, only:ptmass_drift + use ptmass, only:ptmass_drift,use_regnbody + use subgroup, only:evolve_groups use io , only:id,master use mpiutils, only:bcast_mpi - real, intent(in) :: dt,cki - integer, intent(in) :: npart,nptmass,ntypes - real, intent(inout) :: time_par - real, intent(inout) :: xyzh(:,:),vxyzu(:,:) - real, intent(inout) :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:) - integer, optional, intent(in) :: n_ingroup - integer, optional, intent(in) :: group_info(:,:) + real, intent(in) :: dt,cki + integer, intent(in) :: npart,nptmass,ntypes + real, intent(inout) :: time_par + real, intent(inout) :: xyzh(:,:),vxyzu(:,:) + real, intent(inout) :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:) + real, intent(inout) :: fxyz_ptmass(:,:),gtgrad(:,:),bin_info(:,:) + integer, intent(in) :: n_ingroup,n_group + integer, intent(inout) :: group_info(:,:) integer :: i real :: ckdt @@ -644,7 +625,7 @@ subroutine drift(cki,dt,time_par,npart,nptmass,ntypes,xyzh,xyzmh_ptmass,vxyzu,vx ! Drift sink particles if (nptmass>0) then if (id==master) then - if (present(n_ingroup)) then + if (use_regnbody) then call ptmass_drift(nptmass,ckdt,xyzmh_ptmass,vxyz_ptmass,group_info,n_ingroup) else call ptmass_drift(nptmass,ckdt,xyzmh_ptmass,vxyz_ptmass) @@ -653,6 +634,11 @@ subroutine drift(cki,dt,time_par,npart,nptmass,ntypes,xyzh,xyzmh_ptmass,vxyzu,vx call bcast_mpi(xyzmh_ptmass(:,1:nptmass)) endif + if (use_regnbody) then + call evolve_groups(n_group,nptmass,time_par,time_par+cki*dt,group_info,bin_info, & + xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,gtgrad) + endif + time_par = time_par + ckdt !! update time for external potential in force routine end subroutine drift @@ -860,8 +846,8 @@ end subroutine kick !---------------------------------------------------------------- subroutine get_force(nptmass,npart,nsubsteps,ntypes,timei,dtextforce,xyzh,vxyzu, & fext,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,dsdt_ptmass,dt,dki, & - force_count,extf_vdep_flag,linklist_ptmass,bin_info,fsink_old,& - group_info,isionised) + force_count,extf_vdep_flag,linklist_ptmass,bin_info,group_info,& + fsink_old,isionised) use io, only:iverbose,master,id,iprint,warning,fatal use dim, only:maxp,maxvxyzu,itau_alloc,use_apr use ptmass, only:get_accel_sink_gas,get_accel_sink_sink,merge_sinks, & @@ -886,8 +872,8 @@ subroutine get_force(nptmass,npart,nsubsteps,ntypes,timei,dtextforce,xyzh,vxyzu, real, intent(in) :: timei,dki,dt logical, intent(in) :: extf_vdep_flag real, intent(inout) :: bin_info(:,:) + integer, intent(in) :: group_info(:,:) real, optional, intent(inout) :: fsink_old(4,nptmass) - integer, optional, intent(in) :: group_info(:,:) logical, optional, intent(in) :: isionised(:) integer :: merge_ij(nptmass) integer :: merge_n @@ -897,7 +883,7 @@ subroutine get_force(nptmass,npart,nsubsteps,ntypes,timei,dtextforce,xyzh,vxyzu, real :: fextx,fexty,fextz,xi,yi,zi,pmassi,damp_fac real :: fonrmaxi,phii,dtphi2i real :: dkdt,extrapfac - logical :: extrap,last,wsub + logical :: extrap,last if (present(fsink_old)) then fsink_old = fxyz_ptmass @@ -906,13 +892,6 @@ subroutine get_force(nptmass,npart,nsubsteps,ntypes,timei,dtextforce,xyzh,vxyzu, extrap = .false. endif - if (present(group_info)) then - wsub = .true. - else - wsub = .false. - endif - - force_count = force_count + 1 extrapfac = (1./24.)*dt**2 dkdt = dki*dt @@ -935,54 +914,31 @@ subroutine get_force(nptmass,npart,nsubsteps,ntypes,timei,dtextforce,xyzh,vxyzu, if (nptmass > 0) then if (id==master) then if (extrap) then - if (wsub) then + call get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,epot_sinksink,& + dtf,iexternalforce,timei,merge_ij,merge_n,dsdt_ptmass, & + group_info,bin_info,extrapfac,fsink_old) + if (merge_n > 0) then + call merge_sinks(timei,nptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,linklist_ptmass,merge_ij) call get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,epot_sinksink,& - dtf,iexternalforce,timei,merge_ij,merge_n,dsdt_ptmass, & - extrapfac,fsink_old,group_info,bin_info) - if (merge_n > 0) then - call merge_sinks(timei,nptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,linklist_ptmass,merge_ij) - call get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,epot_sinksink,& - dtf,iexternalforce,timei,merge_ij,merge_n,dsdt_ptmass, & - extrapfac,fsink_old,group_info,bin_info) - endif - else - call get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,epot_sinksink,& - dtf,iexternalforce,timei,merge_ij,merge_n, & - dsdt_ptmass,extrapfac,fsink_old) - if (merge_n > 0) then - call merge_sinks(timei,nptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,linklist_ptmass,merge_ij) - call get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,epot_sinksink,& - dtf,iexternalforce,timei,merge_ij,merge_n, & - dsdt_ptmass,extrapfac,fsink_old) - endif + dtf,iexternalforce,timei,merge_ij,merge_n,dsdt_ptmass, & + group_info,bin_info,extrapfac,fsink_old) endif else - if (wsub) then + call get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,epot_sinksink,& + dtf,iexternalforce,timei,merge_ij,merge_n,dsdt_ptmass, & + group_info,bin_info) + if (merge_n > 0) then + call merge_sinks(timei,nptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,linklist_ptmass,merge_ij) call get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,epot_sinksink,& - dtf,iexternalforce,timei,merge_ij,merge_n,dsdt_ptmass, & - group_info=group_info,bin_info=bin_info) - if (merge_n > 0) then - call merge_sinks(timei,nptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,linklist_ptmass,merge_ij) - call get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,epot_sinksink,& - dtf,iexternalforce,timei,merge_ij,merge_n,dsdt_ptmass, & - group_info=group_info,bin_info=bin_info) - endif - if (iverbose >= 2) write(iprint,*) 'dt(sink-sink) = ',C_force*dtf - fxyz_ptmass_sinksink(:,1:nptmass) = fxyz_ptmass (:,1:nptmass) - dsdt_ptmass_sinksink(:,1:nptmass) = dsdt_ptmass (:,1:nptmass) - else - call get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,epot_sinksink,& - dtf,iexternalforce,timei,merge_ij,merge_n,dsdt_ptmass) - if (merge_n > 0) then - call merge_sinks(timei,nptmass,xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,linklist_ptmass,merge_ij) - call get_accel_sink_sink(nptmass,xyzmh_ptmass,fxyz_ptmass,epot_sinksink,& - dtf,iexternalforce,timei,merge_ij,merge_n,dsdt_ptmass) - endif - if (iverbose >= 2) write(iprint,*) 'dt(sink-sink) = ',C_force*dtf - fxyz_ptmass_sinksink(:,1:nptmass) = fxyz_ptmass (:,1:nptmass) - dsdt_ptmass_sinksink(:,1:nptmass) = dsdt_ptmass (:,1:nptmass) + dtf,iexternalforce,timei,merge_ij,merge_n,dsdt_ptmass, & + group_info,bin_info) endif endif + if (iverbose >= 2) write(iprint,*) 'dt(sink-sink) = ',C_force*dtf + if (last) then + fxyz_ptmass_sinksink(:,1:nptmass) = fxyz_ptmass (:,1:nptmass) + dsdt_ptmass_sinksink(:,1:nptmass) = dsdt_ptmass (:,1:nptmass) + endif else fxyz_ptmass(:,1:nptmass) = 0. dsdt_ptmass(:,1:nptmass) = 0. @@ -997,7 +953,7 @@ subroutine get_force(nptmass,npart,nsubsteps,ntypes,timei,dtextforce,xyzh,vxyzu, ! !$omp parallel default(none) & - !$omp shared(maxp,maxphase,wsub) & + !$omp shared(maxp,maxphase) & !$omp shared(npart,nptmass,xyzh,vxyzu,xyzmh_ptmass,fext) & !$omp shared(eos_vars,dust_temp,idamp,damp_fac,abundance,iphase,ntypes,massoftype) & !$omp shared(dkdt,dt,timei,iexternalforce,extf_vdep_flag,last,aprmassoftype,apr_level) & @@ -1034,30 +990,17 @@ subroutine get_force(nptmass,npart,nsubsteps,ntypes,timei,dtextforce,xyzh,vxyzu, zi = xyzh(3,i) endif if (nptmass > 0) then - if (wsub) then - if (extrap) then - call get_accel_sink_gas(nptmass,xi,yi,zi,xyzh(4,i),xyzmh_ptmass,& + if (extrap) then + call get_accel_sink_gas(nptmass,xi,yi,zi,xyzh(4,i),xyzmh_ptmass,& fextx,fexty,fextz,phii,pmassi,fxyz_ptmass, & - dsdt_ptmass,fonrmaxi,dtphi2i,extrapfac,fsink_old,& - bin_info=bin_info) - else - call get_accel_sink_gas(nptmass,xi,yi,zi,xyzh(4,i),xyzmh_ptmass,& - fextx,fexty,fextz,phii,pmassi,fxyz_ptmass,dsdt_ptmass,fonrmaxi,dtphi2i,& - bin_info=bin_info) - fonrmax = max(fonrmax,fonrmaxi) - dtphi2 = min(dtphi2,dtphi2i) - endif + dsdt_ptmass,fonrmaxi,dtphi2i,bin_info,& + extrapfac,fsink_old) else - if (extrap) then - call get_accel_sink_gas(nptmass,xi,yi,zi,xyzh(4,i),xyzmh_ptmass,& - fextx,fexty,fextz,phii,pmassi,fxyz_ptmass, & - dsdt_ptmass,fonrmaxi,dtphi2i,extrapfac,fsink_old) - else - call get_accel_sink_gas(nptmass,xi,yi,zi,xyzh(4,i),xyzmh_ptmass,& - fextx,fexty,fextz,phii,pmassi,fxyz_ptmass,dsdt_ptmass,fonrmaxi,dtphi2i) - fonrmax = max(fonrmax,fonrmaxi) - dtphi2 = min(dtphi2,dtphi2i) - endif + call get_accel_sink_gas(nptmass,xi,yi,zi,xyzh(4,i),xyzmh_ptmass,& + fextx,fexty,fextz,phii,pmassi,fxyz_ptmass,& + dsdt_ptmass,fonrmaxi,dtphi2i,bin_info) + fonrmax = max(fonrmax,fonrmaxi) + dtphi2 = min(dtphi2,dtphi2i) endif endif @@ -1066,9 +1009,9 @@ subroutine get_force(nptmass,npart,nsubsteps,ntypes,timei,dtextforce,xyzh,vxyzu, ! if (iexternalforce > 0) then call get_external_force_gas(xi,yi,zi,xyzh(4,i),vxyzu(1,i), & - vxyzu(2,i),vxyzu(3,i),timei,i, & - dtextforcenew,dtf,dkdt,fextx,fexty,fextz, & - extf_vdep_flag,iexternalforce) + vxyzu(2,i),vxyzu(3,i),timei,i, & + dtextforcenew,dtf,dkdt,fextx,fexty, & + fextz,extf_vdep_flag,iexternalforce) endif ! ! damping @@ -1083,10 +1026,10 @@ subroutine get_force(nptmass,npart,nsubsteps,ntypes,timei,dtextforce,xyzh,vxyzu, if (extrap) then if (itau_alloc == 1) then call get_rad_accel_from_ptmass(nptmass,npart,i,xi,yi,zi,xyzmh_ptmass,fextx,fexty,fextz, & - tau=tau,fsink_old=fsink_old,extrapfac=extrapfac) + tau=tau,fsink_old=fsink_old,extrapfac=extrapfac) else call get_rad_accel_from_ptmass(nptmass,npart,i,xi,yi,zi,xyzmh_ptmass,fextx,fexty,fextz, & - fsink_old=fsink_old,extrapfac=extrapfac) + fsink_old=fsink_old,extrapfac=extrapfac) endif else if (itau_alloc == 1) then @@ -1105,7 +1048,7 @@ subroutine get_force(nptmass,npart,nsubsteps,ntypes,timei,dtextforce,xyzh,vxyzu, ! if (maxvxyzu >= 4 .and. itype==igas .and. last) then call cooling_abundances_update(i,pmassi,xyzh,vxyzu,eos_vars,abundance,nucleation,dust_temp, & - divcurlv,abundc,abunde,abundo,abundsi,dt,dphot0,isionised(i)) + divcurlv,abundc,abunde,abundo,abundsi,dt,dphot0,isionised(i)) endif endif enddo From 0b65f47b7c0878d2004130645892d24cfcf5538e Mon Sep 17 00:00:00 2001 From: Yrisch Date: Mon, 2 Dec 2024 17:32:29 +0100 Subject: [PATCH 12/36] (test_ptmass) fix eos initialisation in test_createsink if used with test_all --- src/tests/test_ptmass.f90 | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/tests/test_ptmass.f90 b/src/tests/test_ptmass.f90 index ce6751ab8..da7f7dac6 100644 --- a/src/tests/test_ptmass.f90 +++ b/src/tests/test_ptmass.f90 @@ -780,7 +780,9 @@ subroutine test_createsink(ntests,npass) use dim, only:gravity,maxp,maxphase use boundary, only:set_boundary use deriv, only:get_derivs_global + use eos, only:ieos,polyk use kdtree, only:tree_accuracy + use units, only:set_units use io, only:id,master,iverbose use part, only:init_part,npart,npartoftype,igas,xyzh,massoftype,hfact,rhoh,& iphase,isetphase,fext,divcurlv,vxyzu,fxyzu,poten, & @@ -801,10 +803,13 @@ subroutine test_createsink(ntests,npass) real :: etotin,angmomin,totmomin,rhomax,rhomax_test procedure(rho_func), pointer :: density_func + call set_units(mass=1.d0,dist=1.d0,G=1.d0) density_func => gaussianr t = 0. iverbose = 1 rho_crit = rho_crit_cgs + ieos = 1 + polyk = 0. do itest=1,3 select case(itest) @@ -822,6 +827,7 @@ subroutine test_createsink(ntests,npass) vxyzu(:,:) = 0. fxyzu(:,:) = 0. fext(:,:) = 0. + ! ! set a boundary that is larger than the sphere size, so test still works with periodic boundaries ! @@ -1318,6 +1324,7 @@ subroutine test_SDAR(ntests,npass) real :: fxyz_sinksink(4,3),dsdt_sinksink(3,3) ! we only use 3 sink particles in the tests here real :: xsec(3),vsec(3) real(kind=4) :: t1 + if (id==master) write(*,"(/,a)") '--> testing SDAR module : Kozai-Lidov effect' ! !--no gas particles ! From 245f2d2be81ea0acad73198bed752af7c70bb295 Mon Sep 17 00:00:00 2001 From: Kateryna Andrych Date: Wed, 4 Dec 2024 15:22:53 +1100 Subject: [PATCH 13/36] set isink=1for ieos=6 in case where isink was not set --- src/main/eos.f90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 5df8d12ca..59c4326c5 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -481,9 +481,8 @@ subroutine init_eos(eos_type,ierr) !--Check that if using ieos=6, then isink is set properly ! if (isink==0) then - call error('eos','ieos=6, but isink is not set') - ierr = ierr_isink_not_set - return + call error('eos','ieos=6, but isink is not set, setting to 1') + isink = 1 endif case(8) From 1a847eb3b93f7ace6120059d6d13298b478b252b Mon Sep 17 00:00:00 2001 From: Kateryna Andrych Date: Wed, 4 Dec 2024 15:50:34 +1100 Subject: [PATCH 14/36] does not override ieos to 3 for single disc in case where ieos==6 --- src/setup/setup_disc.f90 | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/src/setup/setup_disc.f90 b/src/setup/setup_disc.f90 index 590f8a2b4..171b7e404 100644 --- a/src/setup/setup_disc.f90 +++ b/src/setup/setup_disc.f90 @@ -687,6 +687,7 @@ subroutine equation_of_state(gamma) ieos = 6 print "(/,a)",' setting ieos=6 for locally isothermal disc around sink' else + isink = 0 if (discstrat > 0) then ieos = 7 print "(/,a)",' setting ieos=7 for locally isothermal disc with stratification' @@ -695,10 +696,15 @@ subroutine equation_of_state(gamma) polyk2 = (cs*(1./R_ref(onlydisc))**(-qfacdisc2))**2 z0 = z0_ref/R_ref(onlydisc)**beta_z else - ieos = 3 - print "(/,a)",' setting ieos=3 for locally isothermal disc around origin' + if (ieos == 6) then + ! handle the case where ieos=6 is already set in the .in file; do not override this + isink = 1 + print "(/,a)",' keeping ieos=6 for locally isothermal disc with bright primary' + else + ieos = 3 + print "(/,a)",' setting ieos=3 for locally isothermal disc around origin' + endif endif - isink = 0 ! In the case isink==3, to be generalized endif qfacdisc = qindex(onlydisc) endif From 269861aee0bcf86daf901c0818a4b31ec986538a Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Fri, 6 Dec 2024 14:36:20 +1100 Subject: [PATCH 15/36] (set_star) shift all functionality from setup_star into set_star; now includes eos, relax and apr options needed for one or more stars --- src/main/apr.f90 | 14 ++- src/setup/set_star.f90 | 137 ++++++++++++++++++++----- src/setup/setup_star.f90 | 204 +++++-------------------------------- src/tests/test_setstar.f90 | 3 +- 4 files changed, 145 insertions(+), 213 deletions(-) diff --git a/src/main/apr.f90 b/src/main/apr.f90 index 3adf55d70..b645ade72 100644 --- a/src/main/apr.f90 +++ b/src/main/apr.f90 @@ -24,10 +24,12 @@ module apr ! mpiforce, part, physcon, ptmass, quitdump, random, relaxem, ! timestep_ind, vectorutils ! + use dim, only:use_apr implicit none public :: init_apr,update_apr,read_options_apr,write_options_apr public :: create_or_update_apr_clump + public :: use_apr ! default values for runtime parameters integer, public :: apr_max_in = 3 @@ -56,7 +58,7 @@ module apr !+ !----------------------------------------------------------------------- subroutine init_apr(apr_level,ierr) - use dim, only:maxp_hard + use dim, only:maxp_hard use part, only:npart,massoftype,aprmassoftype use apr_region, only:set_apr_centre,set_apr_regions integer, intent(inout) :: ierr @@ -610,22 +612,26 @@ subroutine write_options_apr(iunit) call write_inopt(apr_max_in,'apr_max','number of additional refinement levels (3 -> 2x resolution)',iunit) call write_inopt(ref_dir,'ref_dir','increase (1) or decrease (-1) resolution',iunit) call write_inopt(apr_type,'apr_type','1: static, 2: moving sink, 3: create clumps',iunit) - select case (apr_type) + select case (apr_type) case(2) call write_inopt(track_part,'track_part','number of sink to track',iunit) - case default call write_inopt(apr_centre(1),'apr_centre(1)','centre of region x position',iunit) call write_inopt(apr_centre(2),'apr_centre(2)','centre of region y position',iunit) call write_inopt(apr_centre(3),'apr_centre(3)','centre of region z position',iunit) - end select + call write_inopt(apr_rad,'apr_rad','radius of innermost region',iunit) call write_inopt(apr_drad,'apr_drad','size of step to next region',iunit) end subroutine write_options_apr +!----------------------------------------------------------------------- +!+ +! Find the closest neighbour to a particle (needs replacing) +!+ +!----------------------------------------------------------------------- subroutine closest_neigh(i,next_door,rmin) use part, only:xyzh,npart integer, intent(in) :: i diff --git a/src/setup/set_star.f90 b/src/setup/set_star.f90 index 8997a9cb0..e0ea77a77 100644 --- a/src/setup/set_star.f90 +++ b/src/setup/set_star.f90 @@ -52,7 +52,7 @@ module setstar public :: shift_star,shift_stars public :: write_options_star,write_options_stars public :: read_options_star,read_options_stars - public :: set_star_interactive + public :: set_stars_interactive public :: ikepler,imesa,ibpwpoly,ipoly,iuniform,ifromfile,ievrard public :: need_polyk @@ -102,9 +102,15 @@ end subroutine set_defaults_star !+ !-------------------------------------------------------------------------- subroutine set_defaults_stars(stars) + use eos, only:use_var_comp,gmw,X_in,Z_in type(star_t), intent(out) :: stars(:) integer :: i + EOSopt = 1 + gmw = 0.5988 + X_in = 0.74 + Z_in = 0.02 + use_var_comp = .false. do i=1,size(stars) call set_defaults_star(stars(i)) enddo @@ -432,7 +438,6 @@ subroutine set_stars(id,master,nstars,star,xyzh,vxyzu,eos_vars,rad,& procedure(mask_prototype) :: mask integer :: i - print*,' HERE ieos=',ieos,EOSopt if (ieos==9) call init_eos_piecewise_preset(EOSopt) call init_eos(ieos,ierr) if (ierr /= 0) then @@ -579,15 +584,13 @@ end subroutine write_mass ! This routine should not do ANY prompting !+ !----------------------------------------------------------------------- -subroutine set_defaults_given_profile(iprofile,filename,need_iso,ieos,mstar,polyk) +subroutine set_defaults_given_profile(iprofile,filename,ieos,mstar,polyk) integer, intent(in) :: iprofile character(len=120), intent(out) :: filename - integer, intent(out) :: need_iso integer, intent(inout) :: ieos real, intent(inout) :: polyk character(len=20), intent(out) :: mstar - need_iso = 0 select case(iprofile) case(ifromfile) ! Read the density profile from file (e.g. for neutron star) @@ -605,11 +608,9 @@ subroutine set_defaults_given_profile(iprofile,filename,need_iso,ieos,mstar,poly ! piecewise polytrope ! Original Author: Madeline Marshall & Bernard Field ! Supervisors: James Wurster & Paul Lasky - ieos = 9 - Mstar = '1.35' - polyk = 144. - case(ievrard) - need_iso = -1 + ieos = 9 + Mstar = '1.35' + polyk = 144. end select end subroutine set_defaults_given_profile @@ -619,14 +620,11 @@ end subroutine set_defaults_given_profile ! interactive prompting for setting up a star !+ !----------------------------------------------------------------------- -subroutine set_star_interactive(id,master,star,need_iso,use_var_comp,ieos) +subroutine set_star_interactive(star,ieos) use prompting, only:prompt use setstar_utils, only:nprofile_opts,profile_opt,need_inputprofile,need_rstar use units, only:in_code_units type(star_t), intent(out) :: star - integer, intent(in) :: id,master - logical, intent(out) :: use_var_comp - integer, intent(out) :: need_iso integer, intent(inout) :: ieos integer :: i @@ -642,9 +640,9 @@ subroutine set_star_interactive(id,master,star,need_iso,use_var_comp,ieos) ! ! set default file output parameters ! - if (id==master) write(*,"('Setting up ',a)") trim(profile_opt(star%iprofile)) + write(*,"('Setting up ',a)") trim(profile_opt(star%iprofile)) call set_defaults_given_profile(star%iprofile,star%input_profile,& - need_iso,ieos,star%m,star%polyk) + ieos,star%m,star%polyk) ! resolution if (star%iprofile > 0) then @@ -664,8 +662,6 @@ subroutine set_star_interactive(id,master,star,need_iso,use_var_comp,ieos) select case (star%iprofile) case(imesa) - call prompt('Use variable composition?',use_var_comp) - print*,'Soften the core density profile and add a sink particle core?' print "(3(/,a))",'0: Do not soften profile', & '1: Use cubic softened density profile', & @@ -717,6 +713,44 @@ subroutine set_star_interactive(id,master,star,need_iso,use_var_comp,ieos) end subroutine set_star_interactive +!----------------------------------------------------------------------- +!+ +! as above but wrapper routine for multiple stars +!+ +!----------------------------------------------------------------------- +subroutine set_stars_interactive(star,ieos,relax,nstar) + use prompting, only:prompt + type(star_t), intent(out) :: star(:) + integer, intent(inout) :: ieos + logical, intent(out) :: relax + integer, intent(out), optional :: nstar + integer :: i,nstars + + call set_defaults_stars(star) + + ! optionally ask for number of stars, otherwise fix nstars to the input array size + if (present(nstar) .and. size(star) > 1) then + call prompt('how many stars to set up (0-'//achar(size(star)+48)//')',nstar,0,size(star)) + nstars = nstar + else + nstars = size(star) + endif + + do i=1,nstars + print "(/,'------------- STAR ',i0,'-------------')",i + call set_star_interactive(star(i),ieos) + enddo + + ! prompt for equation of state and relaxation options if any stars made of gas + if (nstars > 0) then + if (any(star(1:nstars)%iprofile > 0)) then + call set_star_eos_interactive(ieos,star) + call prompt('Relax stars automatically during setup?',relax) + endif + endif + +end subroutine set_stars_interactive + !----------------------------------------------------------------------- !+ ! write setupfile options needed for a star @@ -814,12 +848,11 @@ end subroutine write_options_star ! read setupfile options needed for a star !+ !----------------------------------------------------------------------- -subroutine read_options_star(star,need_iso,ieos,db,nerr,label) +subroutine read_options_star(star,ieos,db,nerr,label) use infile_utils, only:inopts,read_inopt use setstar_utils, only:need_inputprofile,need_rstar,nprofile_opts type(star_t), intent(out) :: star type(inopts), allocatable, intent(inout) :: db(:) - integer, intent(out) :: need_iso integer, intent(inout) :: ieos integer, intent(inout) :: nerr character(len=*), intent(in), optional :: label @@ -837,7 +870,7 @@ subroutine read_options_star(star,need_iso,ieos,db,nerr,label) call read_inopt(star%iprofile,'iprofile'//trim(c),db,errcount=nerr,min=0,max=nprofile_opts) call set_defaults_given_profile(star%iprofile,star%input_profile,& - need_iso,ieos,star%m,star%polyk) + ieos,star%m,star%polyk) if (need_inputprofile(star%iprofile)) then call read_inopt(star%input_profile,'input_profile'//trim(c),db,errcount=nerr) @@ -916,12 +949,13 @@ end subroutine read_options_star ! write_options routine that writes options for multiple stars !+ !----------------------------------------------------------------------- -subroutine write_options_stars(star,relax,ieos,iunit,nstar) +subroutine write_options_stars(star,relax,write_rho_to_file,ieos,iunit,nstar) use relaxstar, only:write_options_relax use infile_utils, only:write_inopt + use apr, only:use_apr,write_options_apr type(star_t), intent(in) :: star(:) integer, intent(in) :: ieos,iunit - logical, intent(in) :: relax + logical, intent(in) :: relax,write_rho_to_file integer, intent(in), optional :: nstar integer :: i,nstars @@ -950,9 +984,12 @@ subroutine write_options_stars(star,relax,ieos,iunit,nstar) write(iunit,"(/,a)") '# relaxation options' call write_inopt(relax,'relax','relax stars into equilibrium',iunit) call write_options_relax(iunit) + call write_inopt(write_rho_to_file,'write_rho_to_file','write density profile(s) to file',iunit) endif endif + if (use_apr) call write_options_apr(iunit) + end subroutine write_options_stars !----------------------------------------------------------------------- @@ -960,14 +997,14 @@ end subroutine write_options_stars ! read_options routine that reads options for multiple stars !+ !----------------------------------------------------------------------- -subroutine read_options_stars(star,need_iso,ieos,relax,db,nerr,nstar) +subroutine read_options_stars(star,ieos,relax,write_rho_to_file,db,nerr,nstar) use relaxstar, only:read_options_relax use infile_utils, only:inopts,read_inopt + use apr, only:use_apr,apr_max_in,ref_dir,apr_type,apr_rad,apr_drad type(star_t), intent(out) :: star(:) type(inopts), allocatable, intent(inout) :: db(:) - integer, intent(out) :: need_iso integer, intent(inout) :: ieos - logical, intent(out) :: relax + logical, intent(out) :: relax,write_rho_to_file integer, intent(inout) :: nerr integer, intent(out), optional :: nstar integer :: i,nstars @@ -987,17 +1024,27 @@ subroutine read_options_stars(star,need_iso,ieos,relax,db,nerr,nstar) ! read options for each star do i=1,nstars - call read_options_star(star(i),need_iso,ieos,db,nerr,label=achar(i+48)) + call read_options_star(star(i),ieos,db,nerr,label=achar(i+48)) enddo ! read relaxation options if any stars are made of gas if (nstars > 0) then if (any(star(1:nstars)%iprofile > 0)) then call read_inopt(relax,'relax',db,errcount=nerr) - call read_options_relax(db,nerr) + call read_options_relax(db,nerr) + ! option to write density profile to file + call read_inopt(write_rho_to_file,'write_rho_to_file',db,errcount=nerr) endif endif + if (use_apr) then + call read_inopt(apr_max_in,'apr_max',db,errcount=nerr) + call read_inopt(ref_dir,'ref_dir',db,errcount=nerr) + call read_inopt(apr_type,'apr_type',db,errcount=nerr) + call read_inopt(apr_rad,'apr_rad',db,errcount=nerr) + call read_inopt(apr_drad,'apr_drad',db,errcount=nerr) + endif + end subroutine read_options_stars !----------------------------------------------------------------------- @@ -1069,4 +1116,38 @@ subroutine read_options_stars_eos(star,ieos,db,nerr) end subroutine read_options_stars_eos +!------------------------------------------------------------------------ +!+ +! interactive prompt for equation of state options needed to setup stars +!+ +!------------------------------------------------------------------------ +subroutine set_star_eos_interactive(ieos,star) + use prompting, only:prompt + use eos, only:use_var_comp,X_in,Z_in,irecomb,gamma,gmw + integer, intent(inout) :: ieos + type(star_t), intent(in) :: star(:) + + ! equation of state + call prompt('Enter the desired EoS (1=isothermal,2=adiabatic,10=MESA,12=idealplusrad)',ieos) + if (any(star(:)%iprofile==imesa)) call prompt('Use variable composition?',use_var_comp) + + select case(ieos) + case(9) + write(*,'(a)') 'EOS options: 1=APR3,2=SLy,3=MS1,4=ENG (from Read et al 2009)' + call prompt('Enter equation of state type',EOSopt,1,4) + case(2,12) + call prompt('Enter gamma (adiabatic index)',gamma,1.,7.) + if ( (.not. use_var_comp) .and. any(need_mu(star(:)%isoftcore))) then + call prompt('Enter mean molecular weight',gmw,0.) + endif + case(10,20) + if (ieos==20) call prompt('Enter irecomb (0: H2+H+He, 1:H+He, 2:He)',irecomb,0) + if ( (.not. use_var_comp) .and. any(need_mu(star(:)%isoftcore))) then + call prompt('Enter hydrogen mass fraction (X)',X_in,0.,1.) + call prompt('Enter metals mass fraction (Z)',Z_in,0.,1.) + endif + end select + +end subroutine set_star_eos_interactive + end module setstar diff --git a/src/setup/setup_star.f90 b/src/setup/setup_star.f90 index 1a11292f3..21c8db2d8 100644 --- a/src/setup/setup_star.f90 +++ b/src/setup/setup_star.f90 @@ -38,15 +38,12 @@ module setup use timestep, only:tmax,dtmax use eos, only:ieos use externalforces, only:iext_densprofile - use extern_densprofile, only:nrhotab - use setstar, only:ibpwpoly,ievrard,imesa,star_t,need_polyk + use setstar, only:star_t use setunits, only:dist_unit,mass_unit implicit none ! ! Input parameters ! - integer :: EOSopt - integer :: need_iso real :: maxvxyzu logical :: iexist logical :: relax_star_in_setup,write_rho_to_file @@ -63,15 +60,12 @@ module setup !+ !----------------------------------------------------------------------- subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,time,fileprefix) - use units, only:set_units,select_unit use kernel, only:hfact_default - use eos, only:init_eos,finish_eos,gmw,X_in,Z_in - use eos_piecewise, only:init_eos_piecewise_preset use part, only:nptmass,xyzmh_ptmass,vxyz_ptmass,eos_vars,rad - use mpiutils, only:reduceall_mpi + use eos, only:X_in,Z_in use mpidomain, only:i_belong use setup_params, only:rhozero,npart_total - use setstar, only:set_star + use setstar, only:set_stars,ibpwpoly,ievrard integer, intent(in) :: id integer, intent(inout) :: npart integer, intent(out) :: npartoftype(:) @@ -92,23 +86,11 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, gamma = 5./3. hfact = hfact_default maxvxyzu = size(vxyzu(:,1)) - relax_star_in_setup = .false. - write_rho_to_file = .false. - ! ! set default options ! dist_unit = 'solarr' mass_unit = 'solarm' - EOSopt = 1 - gmw = 0.5988 - X_in = 0.74 - Z_in = 0.02 - use_var_comp = .false. - ! - ! defaults needed for error checking - ! - need_iso = 0 ! -1 = no; 0 = doesn't matter; 1 = yes ! ! determine if the .in file exists ! @@ -125,16 +107,16 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, setupfile = trim(fileprefix)//'.setup' inquire(file=setupfile,exist=setexists) if (setexists) then - call read_setupfile(setupfile,gamma,need_iso,ierr) + call read_setupfile(setupfile,ierr) if (ierr /= 0) then - if (id==master) call write_setupfile(setupfile,gamma,polyk) + if (id==master) call write_setupfile(setupfile) stop 'please rerun phantomsetup with revised .setup file' endif !--Prompt to get inputs and write to file elseif (id==master) then print "(a,/)",trim(setupfile)//' not found: using interactive setup' - call setup_interactive(gamma,iexist,id,master,ierr) - call write_setupfile(setupfile,gamma,polyk) + call setup_interactive(ieos) + call write_setupfile(setupfile) stop 'please check and edit .setup file and rerun phantomsetup' endif @@ -145,31 +127,17 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, iexternalforce = iext_densprofile write_rho_to_file = .true. endif - - if (maxvxyzu > 3 .and. need_iso == 1) call fatal('setup','require ISOTHERMAL=yes') - if (maxvxyzu < 4 .and. need_iso ==-1) call fatal('setup','require ISOTHERMAL=no') - ! - ! initialise the equation of state - ! - if (ieos==9) call init_eos_piecewise_preset(EOSopt) - call init_eos(ieos,ierr) - if (ierr /= 0) call fatal('setup','could not initialise equation of state') ! ! set up particles ! npartoftype(:) = 0 npart = 0 nptmass = 0 - vxyzu = 0.0 - call set_star(id,master,star(1),xyzh,vxyzu,eos_vars,rad,npart,npartoftype,& - massoftype,hfact,xyzmh_ptmass,vxyz_ptmass,nptmass,ieos,gamma,& - X_in,Z_in,relax_star_in_setup,use_var_comp,write_rho_to_file,& - rhozero,npart_total,i_belong,ierr) - ! - ! finish/deallocate equation of state tables - ! - call finish_eos(ieos,ierr) - + call set_stars(id,master,1,star,xyzh,vxyzu,eos_vars,rad,& + npart,npartoftype,massoftype,hfact,& + xyzmh_ptmass,vxyz_ptmass,nptmass,ieos,gamma,X_in,Z_in,& + relax_star_in_setup,use_var_comp,write_rho_to_file,& + rhozero,npart_total,i_belong,ierr) ! ! override some default settings in the .in file for some cases ! @@ -191,47 +159,16 @@ end subroutine setpart ! Ask questions of the user to determine which setup to use !+ !----------------------------------------------------------------------- -subroutine setup_interactive(gamma,iexist,id,master,ierr) - use prompting, only:prompt - use units, only:select_unit - use eos, only:X_in,Z_in,gmw - use eos_gasradrec, only:irecomb - use setstar, only:set_star_interactive +subroutine setup_interactive(ieos) + use setstar, only:set_stars_interactive use setunits, only:set_units_interactive - real, intent(out) :: gamma - logical, intent(in) :: iexist - integer, intent(in) :: id,master - integer, intent(out) :: ierr - - ierr = 0 + integer, intent(inout) :: ieos ! units call set_units_interactive(gr) ! star - !call set_star_interactive(id,master,star,need_iso,use_var_comp,ieos) - -! ! equation of state -! call prompt('Enter the desired EoS (1=isothermal,2=adiabatic,10=MESA,12=idealplusrad)',ieos) -! select case(ieos) -! case(9) -! write(*,'(a)') 'EOS options: 1=APR3,2=SLy,3=MS1,4=ENG (from Read et al 2009)' -! call prompt('Enter equation of state type',EOSopt,1,4) -! case(2) -! call prompt('Enter gamma (adiabatic index)',gamma,1.,7.) -! case(20) -! call prompt('Enter irecomb (0: H2+H+He, 1:H+He, 2:He)',irecomb,0) -! end select - -! if ((.not. use_var_comp) .and. (star%isoftcore<=0)) then -! if ( (ieos==12) .or. (ieos==2) ) call prompt('Enter mean molecular weight',gmw,0.) -! if ( (ieos==10) .or. (ieos==20) ) then -! call prompt('Enter hydrogen mass fraction (X)',X_in,0.,1.) -! call prompt('Enter metals mass fraction (Z)',Z_in,0.,1.) -! endif -! endif - -! call prompt('Relax star automatically during setup?',relax_star_in_setup) + call set_stars_interactive(star,ieos,relax_star_in_setup) end subroutine setup_interactive @@ -240,16 +177,10 @@ end subroutine setup_interactive ! Write setup parameters to input file !+ !----------------------------------------------------------------------- -subroutine write_setupfile(filename,gamma,polyk) - use infile_utils, only:write_inopt - use dim, only:tagline,use_apr - use relaxstar, only:write_options_relax - use eos, only:ieos - use eos_gasradrec, only:irecomb - use setstar, only:write_options_stars,need_polyk +subroutine write_setupfile(filename) + use dim, only:tagline + use setstar, only:write_options_stars use setunits, only:write_options_units - use apr, only:write_options_apr - real, intent(in) :: gamma,polyk character(len=*), intent(in) :: filename integer, parameter :: iunit = 20 @@ -259,43 +190,7 @@ subroutine write_setupfile(filename,gamma,polyk) write(iunit,"(a)") '# input file for Phantom star setup' call write_options_units(iunit,gr) - call write_options_stars(star,relax_star_in_setup,ieos,iunit) - - !write(iunit,"(/,a)") '# equation of state' - !call write_inopt(ieos,'ieos','1=isothermal,2=adiabatic,10=MESA,12=idealplusrad',iunit) - - !if (star%iprofile==imesa) then - ! call write_inopt(use_var_comp,'use_var_comp','Use variable composition (X, Z, mu)',iunit) - !endif - - !select case(ieos) - !case(15) ! Helmholtz - ! call write_inopt(star%initialtemp,'initialtemp','initial temperature of the star',iunit) - !case(9) - ! write(iunit,"(/,a)") '# Piecewise Polytrope default options' - ! call write_inopt(EOSopt,'EOSopt','EOS: 1=APR3,2=SLy,3=MS1,4=ENG (from Read et al 2009)',iunit) - !case(2) - ! call write_inopt(gamma,'gamma','Adiabatic index',iunit) - ! if ((star%isoftcore<=0) .and. (.not. use_var_comp)) call write_inopt(gmw,'mu','mean molecular weight',iunit) - !case(10,20) - ! if (ieos==20) call write_inopt(irecomb,'irecomb','Species to include in recombination (0: H2+H+He, 1:H+He, 2:He',iunit) - ! if ( (.not. use_var_comp) .and. (star%isoftcore <= 0) ) then - ! call write_inopt(X_in,'X','hydrogen mass fraction',iunit) - ! call write_inopt(Z_in,'Z','metallicity',iunit) - ! endif - !case(12) - ! call write_inopt(gamma,'gamma','Adiabatic index',iunit) - ! if ((star%isoftcore<=0) .and. (.not. use_var_comp)) call write_inopt(gmw,'mu','mean molecular weight',iunit) - !end select - - !write(iunit,"(/,a)") '# relaxation options' - !call write_inopt(relax_star_in_setup,'relax_star','relax star(s) automatically during setup',iunit) - !if (relax_star_in_setup) call write_options_relax(iunit) - - call write_inopt(write_rho_to_file,'write_rho_to_file','write density profile(s) to file',iunit) - - if (use_apr) call write_options_apr(iunit) - + call write_options_stars(star,relax_star_in_setup,write_rho_to_file,ieos,iunit) close(iunit) end subroutine write_setupfile @@ -304,21 +199,13 @@ end subroutine write_setupfile ! Read setup parameters from input file !+ !----------------------------------------------------------------------- -subroutine read_setupfile(filename,gamma,need_iso,ierr) - use infile_utils, only:open_db_from_file,inopts,close_db,read_inopt - use io, only:error - use units, only:select_unit - use relaxstar, only:read_options_relax - use eos, only:X_in,Z_in,gmw,ieos - use eos_gasradrec, only:irecomb +subroutine read_setupfile(filename,ierr) + use infile_utils, only:open_db_from_file,inopts,close_db use setstar, only:read_options_stars use setunits, only:read_options_and_set_units - use apr, only:apr_max_in,ref_dir,apr_type,apr_rad,apr_drad - use dim, only:use_apr character(len=*), intent(in) :: filename integer, parameter :: lu = 21 - integer, intent(out) :: need_iso,ierr - real, intent(out) :: gamma + integer, intent(out) :: ierr integer :: nerr type(inopts), allocatable :: db(:) @@ -331,50 +218,7 @@ subroutine read_setupfile(filename,gamma,need_iso,ierr) call read_options_and_set_units(db,nerr,gr) ! star options - call read_options_stars(star,need_iso,ieos,relax_star_in_setup,db,nerr) - - !! equation of state - !call read_inopt(ieos,'ieos',db,errcount=nerr) - !if (star%iprofile==imesa) call read_inopt(use_var_comp,'use_var_comp',db,errcount=nerr) - - !select case(ieos) - !case(15) ! Helmholtz - ! call read_inopt(star%initialtemp,'initialtemp',db,errcount=nerr) - !case(9) - ! call read_inopt(EOSopt,'EOSopt',db,errcount=nerr) - !case(2) - ! call read_inopt(gamma,'gamma',db,errcount=nerr) - ! if ( (.not. use_var_comp) .and. (star%isoftcore <= 0)) call read_inopt(gmw,'mu',db,errcount=nerr) - !case(10,20) - ! if (ieos==20) call read_inopt(irecomb,'irecomb',db,errcount=nerr) - ! ! if softening stellar core, composition is automatically determined at R/2 - ! if ( (.not. use_var_comp) .and. (star%isoftcore <= 0)) then - ! call read_inopt(X_in,'X',db,errcount=nerr) - ! call read_inopt(Z_in,'Z',db,errcount=nerr) - ! endif - !case(12) - ! if softening stellar core, mu is automatically determined at R/2 - ! call read_inopt(gamma,'gamma',db,errcount=nerr) - ! if ( (.not. use_var_comp) .and. (star%isoftcore <= 0)) call read_inopt(gmw,'mu',db,errcount=nerr) -! end select - - !if (need_polyk(star%iprofile)) call read_inopt(polyk,'polyk',db,errcount=nerr) - - ! relax star options - !call read_inopt(relax_star_in_setup,'relax_star',db,errcount=nerr) - !if (relax_star_in_setup) call read_options_relax(db,nerr)! - !if (nerr /= 0) ierr = ierr + 1 - - ! option to write density profile to file - call read_inopt(write_rho_to_file,'write_rho_to_file',db) - - if (use_apr) then - call read_inopt(apr_max_in,'apr_max',db,errcount=nerr) - call read_inopt(ref_dir,'ref_dir',db,errcount=nerr) - call read_inopt(apr_type,'apr_type',db,errcount=nerr) - call read_inopt(apr_rad,'apr_rad',db,errcount=nerr) - call read_inopt(apr_drad,'apr_drad',db,errcount=nerr) - endif + call read_options_stars(star,ieos,relax_star_in_setup,write_rho_to_file,db,nerr) if (nerr > 0) then print "(1x,a,i2,a)",'setup_star: ',nerr,' error(s) during read of setup file' diff --git a/src/tests/test_setstar.f90 b/src/tests/test_setstar.f90 index f9474ca05..a5202e827 100644 --- a/src/tests/test_setstar.f90 +++ b/src/tests/test_setstar.f90 @@ -156,9 +156,10 @@ subroutine test_polytrope(ntests,npass) ! do this test twice, to check the second star relaxes... do i=1,2 if (i==2) x0 = [3.,0.,0.] + call set_star(id,master,star,xyzh,vxyzu,eos_vars,rad,& npart,npartoftype,massoftype,hfact,& - xyzmh_ptmass,vxyz_ptmass,nptmass,ieos,polyk,gamma,X_in,Z_in,& + xyzmh_ptmass,vxyz_ptmass,nptmass,ieos,gamma,X_in,Z_in,& relax=.true.,use_var_comp=.false.,write_rho_to_file=.false.,& rhozero=rhozero,npart_total=ntot,mask=i_belong,ierr=ierr,& write_files=.false.,density_error=rmserr,energy_error=ekin,x0=x0) From fc30fd35762b40fb0d9d52efd6295725a039101a Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Fri, 6 Dec 2024 16:51:58 +1100 Subject: [PATCH 16/36] (set_star) further issues with new set_star calls fixed in binary, grtde and star setups --- src/setup/set_star.f90 | 49 +++++++++++++++++++++++--------------- src/setup/setup_binary.f90 | 11 ++++----- src/setup/setup_grtde.f90 | 38 ++++++++++++++++------------- src/setup/setup_star.f90 | 11 ++++++++- 4 files changed, 67 insertions(+), 42 deletions(-) diff --git a/src/setup/set_star.f90 b/src/setup/set_star.f90 index e0ea77a77..b8bc2b47c 100644 --- a/src/setup/set_star.f90 +++ b/src/setup/set_star.f90 @@ -91,7 +91,7 @@ subroutine set_defaults_star(star) star%np = 1000 star%input_profile = 'P12_Phantom_Profile.data' star%outputfilename = 'mysoftenedstar.dat' - star%dens_profile = 'density-profile.tab' + star%dens_profile = 'density.profile' star%label = '' end subroutine set_defaults_star @@ -112,7 +112,7 @@ subroutine set_defaults_stars(stars) Z_in = 0.02 use_var_comp = .false. do i=1,size(stars) - call set_defaults_star(stars(i)) + if (len_trim(stars(i)%m)==0) call set_defaults_star(stars(i)) enddo end subroutine set_defaults_stars @@ -514,32 +514,40 @@ end subroutine shift_star !----------------------------------------------------------------------- !+ -! As above but shifts all stars to desired positions and velocities +! Shifts all stars to desired positions and velocities !+ !----------------------------------------------------------------------- -subroutine shift_stars(nstar,star,xyzmh_ptmass_in,vxyz_ptmass_in,& - xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,npart,npartoftype,nptmass,corotate) +subroutine shift_stars(nstar,star,x0,v0,& + xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& + npart,npartoftype,nptmass,corotate) + use part, only:ihacc,ihsoft integer, intent(in) :: nstar,npart type(star_t), intent(in) :: star(nstar) + real, intent(in) :: x0(3,nstar),v0(3,nstar) real, intent(inout) :: xyzh(:,:),vxyzu(:,:) - real, intent(in) :: xyzmh_ptmass_in(:,:),vxyz_ptmass_in(:,:) real, intent(inout) :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:) integer, intent(inout) :: nptmass,npartoftype(:) logical, intent(in), optional :: corotate - integer :: i + integer :: i,ierr logical :: do_corotate + real :: rstar,mstar,rcore,mcore,hsoft,lcore,hacc do_corotate = .false. if (present(corotate)) do_corotate = corotate - do i=1,min(nstar,size(xyzmh_ptmass_in(1,:))) + do i=1,nstar if (star(i)%iprofile > 0) then - call shift_star(npart,npartoftype,xyzh,vxyzu,x0=xyzmh_ptmass_in(1:3,i),& - v0=vxyz_ptmass_in(1:3,i),itype=i,corotate=do_corotate) + call shift_star(npart,npartoftype,xyzh,vxyzu,x0=x0(1:3,i),& + v0=v0(1:3,i),itype=i,corotate=do_corotate) else + call get_star_properties_in_code_units(star(i),rstar,mstar,rcore,mcore,hsoft,lcore,hacc,ierr) + nptmass = nptmass + 1 - xyzmh_ptmass(:,nptmass) = xyzmh_ptmass_in(:,i) - vxyz_ptmass(:,nptmass) = vxyz_ptmass_in(:,i) + xyzmh_ptmass(1:3,nptmass) = x0(1:3,i) + xyzmh_ptmass(4,nptmass) = mstar + xyzmh_ptmass(ihsoft,nptmass) = hsoft + xyzmh_ptmass(ihacc, nptmass) = hacc + vxyz_ptmass(1:3,nptmass) = v0(1:3,i) endif enddo @@ -629,7 +637,7 @@ subroutine set_star_interactive(star,ieos) integer :: i ! set defaults - call set_defaults_star(star) + if (len_trim(star%m)==0) call set_defaults_star(star) ! Select sphere & set default values do i = 1, nprofile_opts @@ -654,9 +662,9 @@ subroutine set_star_interactive(star,ieos) if (need_inputprofile(star%iprofile)) then call prompt('Enter file name containing input profile',star%input_profile) else - call prompt('Enter the mass of the star (e.g. 1*msun)',star%m) + call prompt('Enter the mass of the star (e.g. 1*msun)',star%m,noblank=.true.) if (need_rstar(star%iprofile)) then - call prompt('Enter the radius of the star (e.g. 1*rsun)',star%r) + call prompt('Enter the radius of the star (e.g. 1*rsun)',star%r,noblank=.true.) endif endif @@ -756,11 +764,11 @@ end subroutine set_stars_interactive ! write setupfile options needed for a star !+ !----------------------------------------------------------------------- -subroutine write_options_star(star,iunit,ieos,label) +subroutine write_options_star(star,ieos,iunit,label) use infile_utils, only:write_inopt,get_optstring use setstar_utils, only:nprofile_opts,profile_opt,need_inputprofile,need_rstar,need_polyk type(star_t), intent(in) :: star - integer, intent(in) :: iunit,ieos + integer, intent(in) :: ieos,iunit character(len=*), intent(in), optional :: label character(len=120) :: string character(len=10) :: c @@ -861,12 +869,13 @@ subroutine read_options_star(star,ieos,db,nerr,label) real :: rstar,mstar,rcore,mcore,hsoft,lcore,hacc ! set defaults - call set_defaults_star(star) + if (len_trim(star%m)==0) call set_defaults_star(star) ! append optional label e.g. '1', '2' c = '' if (present(label)) c = trim(adjustl(label)) star%label = trim(c) + star%dens_profile = 'relax'//trim(c)//'.profile' call read_inopt(star%iprofile,'iprofile'//trim(c),db,errcount=nerr,min=0,max=nprofile_opts) call set_defaults_given_profile(star%iprofile,star%input_profile,& @@ -969,7 +978,7 @@ subroutine write_options_stars(star,relax,write_rho_to_file,ieos,iunit,nstar) ! write options for each star do i=1,nstars - call write_options_star(star(i),iunit,ieos,label=achar(i+48)) + call write_options_star(star(i),ieos,iunit,label=achar(i+48)) enddo ! write equation of state options if any stars made of gas @@ -1009,6 +1018,8 @@ subroutine read_options_stars(star,ieos,relax,write_rho_to_file,db,nerr,nstar) integer, intent(out), optional :: nstar integer :: i,nstars + call set_defaults_stars(star) + ! optionally ask for number of stars if (present(nstar)) then call read_inopt(nstar,'nstars',db,errcount=nerr,min=0,max=size(star)) diff --git a/src/setup/setup_binary.f90 b/src/setup/setup_binary.f90 index 1291b6805..5bebb3c3f 100644 --- a/src/setup/setup_binary.f90 +++ b/src/setup/setup_binary.f90 @@ -25,7 +25,7 @@ module setup implicit none public :: setpart - logical :: relax,corotate + logical :: relax,write_rho_to_file,corotate type(star_t) :: star(2) type(orbit_t) :: orbit @@ -142,7 +142,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,& ! !--place stars into orbit, or add real sink particles if iprofile=0 ! - call shift_stars(nstar,star,xyzmh_ptmass_in,vxyz_ptmass_in,xyzh,vxyzu,& + call shift_stars(nstar,star,xyzmh_ptmass_in(1:3,:),vxyz_ptmass_in(1:3,:),xyzh,vxyzu,& xyzmh_ptmass,vxyz_ptmass,npart,npartoftype,nptmass,corotate=add_spin) ! !--restore options @@ -188,7 +188,7 @@ subroutine write_setupfile(filename,ieos) write(iunit,"(a)") '# input file for binary setup routines' call write_options_units(iunit,gr) - call write_options_stars(star,relax,ieos,iunit) + call write_options_stars(star,relax,write_rho_to_file,ieos,iunit) call write_inopt(corotate,'corotate','set stars in corotation',iunit) call write_options_orbit(orbit,iunit) close(iunit) @@ -210,15 +210,14 @@ subroutine read_setupfile(filename,ieos,ierr) integer, intent(inout) :: ieos integer, intent(out) :: ierr integer, parameter :: iunit = 21 - integer :: nerr,need_iso + integer :: nerr type(inopts), allocatable :: db(:) nerr = 0 ierr = 0 call open_db_from_file(db,filename,iunit,ierr) call read_options_and_set_units(db,nerr,gr) - call read_options_stars(star,need_iso,ieos,relax,db,nerr) - if (need_iso==1) call fatal('setup_binary','incompatible setup for eos') + call read_options_stars(star,ieos,relax,write_rho_to_file,db,nerr) call read_inopt(corotate,'corotate',db,errcount=nerr) call read_options_orbit(orbit,db,nerr) call close_db(db) diff --git a/src/setup/setup_grtde.f90 b/src/setup/setup_grtde.f90 index 475989343..bb25f7e19 100644 --- a/src/setup/setup_grtde.f90 +++ b/src/setup/setup_grtde.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,igas,& gravity,eos_vars,rad,gr use setbinary, only:set_binary - use setstar, only:set_star,shift_star + use setstar, only:set_star,shift_star,set_defaults_star use units, only:set_units,umass,udist use physcon, only:solarm,pi,solarr use io, only:master,fatal,warning @@ -63,6 +63,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, use setup_params, only:rhozero,npart_total use systemutils, only:get_command_option use options, only:iexternalforce + use units, only:in_code_units integer, intent(in) :: id integer, intent(inout) :: npart integer, intent(out) :: npartoftype(:) @@ -77,7 +78,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, logical :: iexist,write_profile,use_var_comp real :: rtidal,rp,semia,period,hacc1,hacc2 real :: vxyzstar(3),xyzstar(3) - real :: r0,vel,lorentz + real :: r0,vel,lorentz,mstar,rstar real :: vhat(3),x0,y0 ! !-- general parameters @@ -101,8 +102,9 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, ! mhole = 1.e6 ! (solar masses) call set_units(mass=mhole*solarm,c=1.d0,G=1.d0) !--Set central mass to M=1 in code units - star%mstar = 1.*solarm/umass - star%rstar = 1.*solarr/udist + call set_defaults_star(star) + star%m = '1.*msun' + star%r = '1.*solarr' np_default = 1e6 star%np = int(get_command_option('np',default=np_default)) ! can set default value with --np=1e5 flag (mainly for testsuite) star%iprofile = 2 @@ -120,7 +122,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, if (id==master) print "(/,65('-'),1(/,a),/,65('-'),/)",' Tidal disruption in GR' filename = trim(fileprefix)//'.setup' inquire(file=filename,exist=iexist) - if (iexist) call read_setupfile(filename,ieos,polyk,ierr) + if (iexist) call read_setupfile(filename,ieos,ierr) if (.not. iexist .or. ierr /= 0) then if (id==master) then call write_setupfile(filename) @@ -133,24 +135,28 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, !--set up and relax a star ! call set_star(id,master,star,xyzh,vxyzu,eos_vars,rad,npart,npartoftype,& - massoftype,hfact,xyzmh_ptmass,vxyz_ptmass,nptmass,ieos,polyk,gamma,& + massoftype,hfact,xyzmh_ptmass,vxyz_ptmass,nptmass,ieos,gamma,& X_in,Z_in,relax,use_var_comp,write_profile,& rhozero,npart_total,i_belong,ierr) if (ierr /= 0) call fatal('setup','errors in set_star') + rstar = in_code_units(star%r,ierr,unit_type='length') + if (ierr /= 0) call fatal('setup','could not convert rstar to code units') + mstar = in_code_units(star%m,ierr,unit_type='mass') + if (ierr /= 0) call fatal('setup','could not convert mstar to code units') ! !--place star into orbit ! - rtidal = star%rstar*(mass1/star%mstar)**(1./3.) + rtidal = rstar*(mass1/mstar)**(1./3.) rp = rtidal/beta accradius1_hard = 5.*mass1 accradius1 = accradius1_hard a = 0. theta = theta*pi/180. - print*, 'mstar', star%mstar - print*, 'rstar', star%rstar + print*, 'mstar', mstar + print*, 'rstar', rstar print*, 'umass', umass print*, 'udist', udist print*, 'mass1', mass1 @@ -168,11 +174,11 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, semia = rp/(1.-ecc) period = 2.*pi*sqrt(semia**3/mass1) print*, 'period', period - hacc1 = star%rstar/1.e8 ! Something small so that set_binary doesnt warn about Roche lobe + hacc1 = rstar/1.e8 ! Something small so that set_binary doesnt warn about Roche lobe hacc2 = hacc1 ! apocentre = rp*(1.+ecc)/(1.-ecc) ! trueanom = acos((rp*(1.+ecc)/r0 - 1.)/ecc)*180./pi - call set_binary(mass1,star%mstar,semia,ecc,hacc1,hacc2,xyzmh_ptmass,vxyz_ptmass,nptmass,ierr,& + call set_binary(mass1,mstar,semia,ecc,hacc1,hacc2,xyzmh_ptmass,vxyz_ptmass,nptmass,ierr,& posang_ascnode=0.,arg_peri=90.,incl=0.,f=-180.) vxyzstar = vxyz_ptmass(1:3,2) xyzstar = xyzmh_ptmass(1:3,2) @@ -244,13 +250,14 @@ subroutine write_setupfile(filename) use infile_utils, only:write_inopt use setstar, only:write_options_star use relaxstar, only:write_options_relax + use eos, only:ieos character(len=*), intent(in) :: filename integer :: iunit print "(a)",' writing setup options file '//trim(filename) open(newunit=iunit,file=filename,status='replace',form='formatted') write(iunit,"(a)") '# input file for tidal disruption setup' - call write_options_star(star,iunit) + call write_options_star(star,ieos,iunit) call write_inopt(relax,'relax','relax star into hydrostatic equilibrium',iunit) if (relax) call write_options_relax(iunit) @@ -265,7 +272,7 @@ subroutine write_setupfile(filename) end subroutine write_setupfile -subroutine read_setupfile(filename,ieos,polyk,ierr) +subroutine read_setupfile(filename,ieos,ierr) use infile_utils, only:open_db_from_file,inopts,read_inopt,close_db use io, only:error use setstar, only:read_options_star @@ -274,10 +281,9 @@ subroutine read_setupfile(filename,ieos,polyk,ierr) use units, only:set_units character(len=*), intent(in) :: filename integer, intent(inout) :: ieos - real, intent(inout) :: polyk integer, intent(out) :: ierr integer, parameter :: iunit = 21 - integer :: nerr,need_iso + integer :: nerr type(inopts), allocatable :: db(:) print "(a)",'reading setup options from '//trim(filename) @@ -292,7 +298,7 @@ subroutine read_setupfile(filename,ieos,polyk,ierr) ! !--read star options and convert to code units ! - call read_options_star(star,need_iso,ieos,polyk,db,nerr) + call read_options_star(star,ieos,db,nerr) call read_inopt(relax,'relax',db,errcount=nerr) if (relax) call read_options_relax(db,nerr) diff --git a/src/setup/setup_star.f90 b/src/setup/setup_star.f90 index 21c8db2d8..563c2be8f 100644 --- a/src/setup/setup_star.f90 +++ b/src/setup/setup_star.f90 @@ -65,7 +65,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, use eos, only:X_in,Z_in use mpidomain, only:i_belong use setup_params, only:rhozero,npart_total - use setstar, only:set_stars,ibpwpoly,ievrard + use setstar, only:set_stars,shift_stars,ibpwpoly,ievrard integer, intent(in) :: id integer, intent(inout) :: npart integer, intent(out) :: npartoftype(:) @@ -78,6 +78,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, integer :: ierr logical :: setexists character(len=120) :: setupfile,inname + real :: x0(3,1),v0(3,1) ! ! Initialise parameters, including those that will not be included in *.setup ! @@ -139,6 +140,14 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, relax_star_in_setup,use_var_comp,write_rho_to_file,& rhozero,npart_total,i_belong,ierr) ! + ! put the star at the origin with zero velocity, + ! or replace with sink particle + ! + x0 = 0. + v0 = 0. + call shift_stars(1,star,x0,v0,xyzh,vxyzu,& + xyzmh_ptmass,vxyz_ptmass,npart,npartoftype,nptmass) + ! ! override some default settings in the .in file for some cases ! select case(star(1)%iprofile) From d92de15efff25d04c0dcf9c0af037b718b14cd36 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Fri, 6 Dec 2024 18:00:22 +1100 Subject: [PATCH 17/36] (set_star) bug fix with default values for mstar and rstar; now give defaults in msun and rsun; also fix issues with adding stars to grdisc --- src/setup/set_star.f90 | 22 +++----- src/setup/setup_grdisc.F90 | 103 ++++++++++++++++++++----------------- src/setup/setup_grtde.f90 | 2 +- src/setup/setup_star.f90 | 1 + 4 files changed, 63 insertions(+), 65 deletions(-) diff --git a/src/setup/set_star.f90 b/src/setup/set_star.f90 index b8bc2b47c..2b601342d 100644 --- a/src/setup/set_star.f90 +++ b/src/setup/set_star.f90 @@ -75,8 +75,8 @@ subroutine set_defaults_star(star) type(star_t), intent(out) :: star star%iprofile = 2 - star%r = '1.0' - star%m = '1.0' + star%r = '1.0*rsun' + star%m = '1.0*msun' star%ui_coef = 0.05 star%polyk = 0. star%initialtemp = 1.0e7 @@ -112,7 +112,7 @@ subroutine set_defaults_stars(stars) Z_in = 0.02 use_var_comp = .false. do i=1,size(stars) - if (len_trim(stars(i)%m)==0) call set_defaults_star(stars(i)) + call set_defaults_star(stars(i)) enddo end subroutine set_defaults_stars @@ -636,9 +636,6 @@ subroutine set_star_interactive(star,ieos) integer, intent(inout) :: ieos integer :: i - ! set defaults - if (len_trim(star%m)==0) call set_defaults_star(star) - ! Select sphere & set default values do i = 1, nprofile_opts write(*,"(i2,')',1x,a)") i, profile_opt(i) @@ -734,8 +731,6 @@ subroutine set_stars_interactive(star,ieos,relax,nstar) integer, intent(out), optional :: nstar integer :: i,nstars - call set_defaults_stars(star) - ! optionally ask for number of stars, otherwise fix nstars to the input array size if (present(nstar) .and. size(star) > 1) then call prompt('how many stars to set up (0-'//achar(size(star)+48)//')',nstar,0,size(star)) @@ -859,7 +854,7 @@ end subroutine write_options_star subroutine read_options_star(star,ieos,db,nerr,label) use infile_utils, only:inopts,read_inopt use setstar_utils, only:need_inputprofile,need_rstar,nprofile_opts - type(star_t), intent(out) :: star + type(star_t), intent(inout) :: star type(inopts), allocatable, intent(inout) :: db(:) integer, intent(inout) :: ieos integer, intent(inout) :: nerr @@ -868,9 +863,6 @@ subroutine read_options_star(star,ieos,db,nerr,label) integer :: ierr real :: rstar,mstar,rcore,mcore,hsoft,lcore,hacc - ! set defaults - if (len_trim(star%m)==0) call set_defaults_star(star) - ! append optional label e.g. '1', '2' c = '' if (present(label)) c = trim(adjustl(label)) @@ -949,7 +941,7 @@ subroutine read_options_star(star,ieos,db,nerr,label) endif ! perform a unit conversion, just to check that there are no errors parsing the .setup file - call get_star_properties_in_code_units(star,rstar,mstar,rcore,mcore,hsoft,lcore,hacc,nerr) + if (nerr==0) call get_star_properties_in_code_units(star,rstar,mstar,rcore,mcore,hsoft,lcore,hacc,nerr) end subroutine read_options_star @@ -1010,7 +1002,7 @@ subroutine read_options_stars(star,ieos,relax,write_rho_to_file,db,nerr,nstar) use relaxstar, only:read_options_relax use infile_utils, only:inopts,read_inopt use apr, only:use_apr,apr_max_in,ref_dir,apr_type,apr_rad,apr_drad - type(star_t), intent(out) :: star(:) + type(star_t), intent(inout) :: star(:) ! inout because can set default options manually in calling routine type(inopts), allocatable, intent(inout) :: db(:) integer, intent(inout) :: ieos logical, intent(out) :: relax,write_rho_to_file @@ -1018,8 +1010,6 @@ subroutine read_options_stars(star,ieos,relax,write_rho_to_file,db,nerr,nstar) integer, intent(out), optional :: nstar integer :: i,nstars - call set_defaults_stars(star) - ! optionally ask for number of stars if (present(nstar)) then call read_inopt(nstar,'nstars',db,errcount=nerr,min=0,max=size(star)) diff --git a/src/setup/setup_grdisc.F90 b/src/setup/setup_grdisc.F90 index 73d794714..00ef2c139 100644 --- a/src/setup/setup_grdisc.F90 +++ b/src/setup/setup_grdisc.F90 @@ -41,7 +41,7 @@ module setup real, private :: mhole,mdisc,r_in,r_out,r_ref,spin,honr,theta,p_index,q_index,accrad,gamma_ad integer, private :: np,nstars - logical, private :: ismooth,relax + logical, private :: ismooth,relax,write_rho_to_file integer, parameter :: max_stars = 10 type(star_t), private :: star(max_stars) type(orbit_t),private :: orbit(max_stars) @@ -61,7 +61,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, use io, only:master use externalforces, only:accradius1,accradius1_hard use options, only:iexternalforce,alphau,iexternalforce,ipdv_heating,ishock_heating - use units, only:set_units,umass + use units, only:set_units,umass,in_code_units use physcon, only:solarm,pi #ifdef GR use metric, only:a @@ -73,7 +73,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, use timestep, only:tmax,dtmax use eos, only:ieos,use_var_comp,X_in,Z_in use kernel, only:hfact_default - use setstar, only:shift_star,set_stars + use setstar, only:shift_star,set_stars,set_defaults_stars use setorbit, only:set_defaults_orbit,set_orbit use setunits, only:mass_unit use mpidomain, only:i_belong @@ -91,20 +91,15 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, integer :: ierr,nptmass_in,i integer(kind=8) :: npart_total logical :: iexist,write_profile - real :: cs2 + real :: cs2,mstar,rstar real :: xyzmh_ptmass_in(nsinkproperties,2),vxyz_ptmass_in(3,2) time = 0. alphau = 0.0 npartoftype(:) = 0 nptmass = 0 - iexternalforce = 1 hfact = hfact_default -#ifndef GR - iexternalforce = iext_einsteinprec -#endif - tmax = 2.e4 dtmax = 100. @@ -135,6 +130,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, ! stars nstars = 0 + call set_defaults_stars(star) do i=1,size(orbit) call set_defaults_orbit(orbit(i)) enddo @@ -146,16 +142,15 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, if (id==master) print "(/,65('-'),(/,1x,a),/,65('-'),/)",'General relativistic disc setup' filename = trim(fileprefix)//'.setup' inquire(file=filename,exist=iexist) - if (iexist) call read_setupfile(filename,ierr) + if (iexist) call read_setupfile(filename,ieos,ierr) if (.not. iexist .or. ierr /= 0) then if (id==master) then - call write_setupfile(filename) + call write_setupfile(filename,ieos) print*,' Edit '//trim(filename)//' and rerun phantomsetup' endif stop endif accradius1 = accrad - npart = np !-- Set gamma from the option read from .setup file gamma = gamma_ad @@ -167,14 +162,50 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, call set_units(G=1.,c=1.,mass=mhole) ! Set central mass to M=1 in code units mdisc = mdisc*solarm/umass accradius1_hard = accradius1 + massoftype(igas) = mdisc/np ! set particle mass from the disc mass + ! + ! add stars on desired orbits around the black hole, these could be + ! either sink particles or balls of gas + ! + if (nstars > 0) then + write_profile = .false. + iexternalforce = 0 + call set_stars(id,master,nstars,star,xyzh,vxyzu,eos_vars,rad,npart,npartoftype,& + massoftype,hfact,xyzmh_ptmass,vxyz_ptmass,nptmass,ieos,gamma,& + X_in,Z_in,relax,use_var_comp,write_profile,& + rhozero,npart_total,i_belong,ierr) + do i=1,nstars + nptmass_in = 0 + ! convert stellar mass and radius to code units + mstar = in_code_units(star(i)%m,ierr) + rstar = in_code_units(star(i)%r,ierr) + call set_orbit(orbit(i),mhole/umass,mstar,r_in,rstar, & + xyzmh_ptmass_in,vxyz_ptmass_in,nptmass_in,(id==master),ierr) + + ! shift the star to the position of the second body + if (star(i)%iprofile > 0) then + call shift_star(npart,npartoftype,xyzh,vxyzu,& + x0=xyzmh_ptmass_in(:,2),v0=vxyz_ptmass_in(:,2),itype=i) + else + nptmass = nptmass + 1 + xyzmh_ptmass(:,nptmass) = xyzmh_ptmass_in(:,2) + vxyz_ptmass(:,nptmass) = vxyz_ptmass_in(:,2) + endif + enddo +endif + +#ifndef GR + iexternalforce = iext_einsteinprec +#endif ! ! Convert to radians ! theta = theta/180. * pi call set_disc(id,master,& - npart = npart, & + npart = np, & + npart_start = npart+1, & rmin = r_in, & rmax = r_out, & rref = r_ref, & @@ -194,6 +225,8 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, bh_spin = spin, & prefix = fileprefix) + npart = npart + np + #ifdef GR a = spin ! Overwrite thermal energies to be correct for GR @@ -210,33 +243,6 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, npartoftype(1) = npart - ! - ! add stars on desired orbits around the black hole, these could be - ! either sink particles or balls of gas - ! - if (nstars > 0) then - write_profile = .false. - call set_stars(id,master,nstars,star,xyzh,vxyzu,eos_vars,rad,npart,npartoftype,& - massoftype,hfact,xyzmh_ptmass,vxyz_ptmass,nptmass,ieos,polyk,gamma,& - X_in,Z_in,relax,use_var_comp,write_profile,& - rhozero,npart_total,i_belong,ierr) - do i=1,nstars - nptmass_in = 0 - call set_orbit(orbit(i),mhole/umass,star(i)%mstar,r_in,star(i)%rstar, & - xyzmh_ptmass_in,vxyz_ptmass_in,nptmass_in,(id==master),ierr) - - ! shift the star to the position of the second body - if (star(i)%iprofile > 0) then - call shift_star(npart,npartoftype,xyzh,vxyzu,& - x0=xyzmh_ptmass_in(:,2),v0=vxyz_ptmass_in(:,2),itype=i) - else - nptmass = nptmass + 1 - xyzmh_ptmass(:,nptmass) = xyzmh_ptmass_in(:,2) - vxyz_ptmass(:,nptmass) = vxyz_ptmass_in(:,2) - endif - enddo - endif - ipdv_heating = 0 ishock_heating = 0 if (id==master) print "(/,a,/)",' ** SETTING ipdv_heating=0 and ishock_heating=0 for grdisc setup **' @@ -247,12 +253,13 @@ end subroutine setpart ! !---Read/write setup file-------------------------------------------------- ! -subroutine write_setupfile(filename) +subroutine write_setupfile(filename,ieos) use infile_utils, only:write_inopt use setstar, only:write_options_stars use setorbit, only:write_options_orbit use setunits, only:write_options_units character(len=*), intent(in) :: filename + integer, intent(in) :: ieos integer, parameter :: iunit = 20 integer :: i @@ -278,7 +285,7 @@ subroutine write_setupfile(filename) call write_inopt(np ,'np' ,'number of particles in disc' , iunit) write(iunit,"(/,a)") '# stars' - call write_options_stars(star,relax,iunit,nstar=nstars) + call write_options_stars(star,relax,write_rho_to_file,ieos,iunit,nstar=nstars) do i=1,nstars call write_options_orbit(orbit(i),iunit,label=achar(i+48)) enddo @@ -286,17 +293,17 @@ subroutine write_setupfile(filename) end subroutine write_setupfile -subroutine read_setupfile(filename,ierr) +subroutine read_setupfile(filename,ieos,ierr) use infile_utils, only:open_db_from_file,inopts,read_inopt,close_db use io, only:error use setstar, only:read_options_stars use setorbit, only:read_options_orbit - use eos, only:ieos,polyk use setunits, only:read_options_and_set_units - character(len=*), intent(in) :: filename - integer, intent(out) :: ierr + character(len=*), intent(in) :: filename + integer, intent(inout) :: ieos + integer, intent(out) :: ierr integer, parameter :: iunit = 21 - integer :: nerr,need_iso,i + integer :: nerr,i type(inopts), allocatable :: db(:) print "(a)",'reading setup options from '//trim(filename) @@ -319,7 +326,7 @@ subroutine read_setupfile(filename,ierr) call read_inopt(gamma_ad,'gamma' ,db,min=1.,errcount=nerr) call read_inopt(accrad ,'accrad' ,db,min=0.,errcount=nerr) call read_inopt(np ,'np ' ,db,min=0 ,errcount=nerr) - call read_options_stars(star,need_iso,ieos,polyk,relax,db,nerr,nstars) + call read_options_stars(star,ieos,relax,write_rho_to_file,db,nerr,nstars) do i=1,nstars call read_options_orbit(orbit(i),db,nerr,label=achar(i+48)) enddo diff --git a/src/setup/setup_grtde.f90 b/src/setup/setup_grtde.f90 index bb25f7e19..c0d251cb1 100644 --- a/src/setup/setup_grtde.f90 +++ b/src/setup/setup_grtde.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,igas,& gravity,eos_vars,rad,gr use setbinary, only:set_binary - use setstar, only:set_star,shift_star,set_defaults_star + use setstar, only:set_star,shift_star,set_defaults_stars use units, only:set_units,umass,udist use physcon, only:solarm,pi,solarr use io, only:master,fatal,warning diff --git a/src/setup/setup_star.f90 b/src/setup/setup_star.f90 index 563c2be8f..ce9236b90 100644 --- a/src/setup/setup_star.f90 +++ b/src/setup/setup_star.f90 @@ -92,6 +92,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, ! dist_unit = 'solarr' mass_unit = 'solarm' + call set_defaults_stars(star) ! ! determine if the .in file exists ! From 85a3d9018533c84ca6a30abeead6ae84f6c9e6af Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Fri, 6 Dec 2024 20:44:13 +1100 Subject: [PATCH 18/36] (setstar) build and test failures fixed --- src/setup/setup_star.f90 | 2 +- src/tests/test_setstar.f90 | 8 +++++--- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/setup/setup_star.f90 b/src/setup/setup_star.f90 index ce9236b90..288b8b119 100644 --- a/src/setup/setup_star.f90 +++ b/src/setup/setup_star.f90 @@ -65,7 +65,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, use eos, only:X_in,Z_in use mpidomain, only:i_belong use setup_params, only:rhozero,npart_total - use setstar, only:set_stars,shift_stars,ibpwpoly,ievrard + use setstar, only:set_defaults_stars,set_stars,shift_stars,ibpwpoly,ievrard integer, intent(in) :: id integer, intent(inout) :: npart integer, intent(out) :: npartoftype(:) diff --git a/src/tests/test_setstar.f90 b/src/tests/test_setstar.f90 index a5202e827..43f6b3468 100644 --- a/src/tests/test_setstar.f90 +++ b/src/tests/test_setstar.f90 @@ -131,13 +131,13 @@ subroutine test_polytrope(ntests,npass) use mpidomain, only:i_belong use options, only:ieos use physcon, only:solarr,solarm,pi - use eos, only:gamma,X_in,Z_in + use eos, only:gamma,X_in,Z_in,polyk use setstar, only:star_t,set_star,set_defaults_star,ipoly use units, only:set_units use checksetup, only:check_setup integer, intent(inout) :: ntests,npass type(star_t) :: star - real :: polyk,rhozero,rmserr,ekin,x0(3) + real :: rhozero,rmserr,ekin,x0(3) integer(kind=8) :: ntot integer :: ierr,nfail(1),i,nerror,nwarn @@ -148,7 +148,6 @@ subroutine test_polytrope(ntests,npass) call set_units(dist=solarr,mass=solarm,G=1.d0) ieos = 2 gamma = 5./3. - polyk = 1. call set_defaults_star(star) star%iprofile = ipoly ! a polytrope star%np = 1000 @@ -174,6 +173,9 @@ subroutine test_polytrope(ntests,npass) call checkval(rhozero,1./(4./3.*pi),1e-6,nfail(1),'mean density') call update_test_scores(ntests,nfail,npass) + call checkval(star%polyk,0.424304,1e-6,nfail(1),'polyk value for M=1,R=1') + call update_test_scores(ntests,nfail,npass) + call checkval(polyk,0.424304,1e-6,nfail(1),'polyk value for M=1,R=1') call update_test_scores(ntests,nfail,npass) From 4ceb2a5f5bb933d002b742efac802cd9b21a16eb Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Fri, 6 Dec 2024 20:58:47 +1100 Subject: [PATCH 19/36] (setstar) build and test failures fixed --- src/setup/setup_grtde.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/setup/setup_grtde.f90 b/src/setup/setup_grtde.f90 index c0d251cb1..bb25f7e19 100644 --- a/src/setup/setup_grtde.f90 +++ b/src/setup/setup_grtde.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,igas,& gravity,eos_vars,rad,gr use setbinary, only:set_binary - use setstar, only:set_star,shift_star,set_defaults_stars + use setstar, only:set_star,shift_star,set_defaults_star use units, only:set_units,umass,udist use physcon, only:solarm,pi,solarr use io, only:master,fatal,warning From 51f4661c5dbbb6927431cffe6cb751237c098813 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 9 Dec 2024 09:47:38 +1100 Subject: [PATCH 20/36] (wd) test failures fixed; also fixed issues with iprofile=0 in star setup --- src/main/utils_infiles.f90 | 49 ++++++++++++-- src/setup/set_star.f90 | 127 ++++++++++++++++++----------------- src/setup/set_star_utils.f90 | 6 +- 3 files changed, 115 insertions(+), 67 deletions(-) diff --git a/src/main/utils_infiles.f90 b/src/main/utils_infiles.f90 index 8b788ff80..4134a8114 100644 --- a/src/main/utils_infiles.f90 +++ b/src/main/utils_infiles.f90 @@ -21,6 +21,7 @@ module infile_utils public :: write_inopt, read_inopt public :: read_next_inopt, get_inopt public :: write_infile_series, check_infile, contains_loop, get_optstring + public :: int_to_string ! ! generic interface write_inopt to write an input option of any type ! @@ -1259,13 +1260,14 @@ end subroutine write_infile_lines ! Creates a string out of a list of options ! !--------------------------------------------------------------------------- -subroutine get_optstring(nopts,optstring,string,maxlen) +subroutine get_optstring(nopts,optstring,string,maxlen,from_zero) integer, intent(in) :: nopts character(len=*), intent(in) :: optstring(nopts) character(len=*), intent(out) :: string integer, intent(in), optional :: maxlen + logical, intent(in), optional :: from_zero character(len=len(string)) :: temp - integer :: i,maxl,ierr + integer :: i,maxl,ierr,ioffset if (present(maxlen)) then maxl = max(maxlen,1) @@ -1274,15 +1276,54 @@ subroutine get_optstring(nopts,optstring,string,maxlen) endif string = '' + !--allow for enumeration that starts from 0 instead of 1 + ioffset = 0 + if (present(from_zero)) then + if (from_zero) ioffset = 1 + endif + do i=1,nopts temp = adjustl(optstring(i)) if (i==nopts) then - write(string(len_trim(string)+1:),"(i0,'=',a)",iostat=ierr) i,trim(temp(1:maxl)) + write(string(len_trim(string)+1:),"(i0,'=',a)",iostat=ierr) i-ioffset,trim(temp(1:maxl)) else - write(string(len_trim(string)+1:),"(i0,'=',a,',')",iostat=ierr) i,trim(temp(1:maxl)) + write(string(len_trim(string)+1:),"(i0,'=',a,',')",iostat=ierr) i-ioffset,trim(temp(1:maxl)) endif enddo end subroutine get_optstring +!--------------------------------------------------------------------------- +! +! convert an integer to a string without using write statements +! so the function itself can be used in a print or write statement +! +!--------------------------------------------------------------------------- +function int_to_string(num) result(str) + integer, intent(in) :: num + character(len=20) :: str + integer :: i, n + character(len=1) :: digit + + n = abs(num) ! Get the absolute value of the number + str = '' ! Initialize the string + + ! Convert integer to string + do while (n > 0) + i = mod(n, 10) ! Get the last digit + digit = char(i + ichar('0')) ! Convert digit to character + str = trim(adjustl(digit)) // str ! Prepend digit to string + n = n / 10 ! Remove the last digit + enddo + + if (num < 0) then + str = '-' // trim(str) ! Add negative sign if necessary + endif + + if (num == 0) then + str = '0' ! Handle the case for zero + endif + +end function int_to_string + end module infile_utils diff --git a/src/setup/set_star.f90 b/src/setup/set_star.f90 index 2b601342d..a514bb1c1 100644 --- a/src/setup/set_star.f90 +++ b/src/setup/set_star.f90 @@ -592,12 +592,11 @@ end subroutine write_mass ! This routine should not do ANY prompting !+ !----------------------------------------------------------------------- -subroutine set_defaults_given_profile(iprofile,filename,ieos,mstar,polyk) +subroutine set_defaults_given_profile(iprofile,filename,mstar,polyk) integer, intent(in) :: iprofile character(len=120), intent(out) :: filename - integer, intent(inout) :: ieos real, intent(inout) :: polyk - character(len=20), intent(out) :: mstar + character(len=*), intent(inout) :: mstar select case(iprofile) case(ifromfile) @@ -616,8 +615,7 @@ subroutine set_defaults_given_profile(iprofile,filename,ieos,mstar,polyk) ! piecewise polytrope ! Original Author: Madeline Marshall & Bernard Field ! Supervisors: James Wurster & Paul Lasky - ieos = 9 - Mstar = '1.35' + Mstar = '1.35*msun' polyk = 144. end select @@ -628,26 +626,24 @@ end subroutine set_defaults_given_profile ! interactive prompting for setting up a star !+ !----------------------------------------------------------------------- -subroutine set_star_interactive(star,ieos) +subroutine set_star_interactive(star) use prompting, only:prompt use setstar_utils, only:nprofile_opts,profile_opt,need_inputprofile,need_rstar use units, only:in_code_units - type(star_t), intent(out) :: star - integer, intent(inout) :: ieos + type(star_t), intent(inout) :: star integer :: i ! Select sphere & set default values - do i = 1, nprofile_opts + do i = 0, nprofile_opts write(*,"(i2,')',1x,a)") i, profile_opt(i) enddo - call prompt('Enter which density profile to use',star%iprofile,1,nprofile_opts) + call prompt('Enter which density profile to use',star%iprofile,0,nprofile_opts) ! ! set default file output parameters ! write(*,"('Setting up ',a)") trim(profile_opt(star%iprofile)) - call set_defaults_given_profile(star%iprofile,star%input_profile,& - ieos,star%m,star%polyk) + call set_defaults_given_profile(star%iprofile,star%input_profile,star%m,star%polyk) ! resolution if (star%iprofile > 0) then @@ -724,8 +720,9 @@ end subroutine set_star_interactive !+ !----------------------------------------------------------------------- subroutine set_stars_interactive(star,ieos,relax,nstar) - use prompting, only:prompt - type(star_t), intent(out) :: star(:) + use prompting, only:prompt + use infile_utils, only:int_to_string + type(star_t), intent(inout) :: star(:) integer, intent(inout) :: ieos logical, intent(out) :: relax integer, intent(out), optional :: nstar @@ -733,7 +730,7 @@ subroutine set_stars_interactive(star,ieos,relax,nstar) ! optionally ask for number of stars, otherwise fix nstars to the input array size if (present(nstar) .and. size(star) > 1) then - call prompt('how many stars to set up (0-'//achar(size(star)+48)//')',nstar,0,size(star)) + call prompt('how many stars to set up (0-'//int_to_string(size(star))//')',nstar,0,size(star)) nstars = nstar else nstars = size(star) @@ -741,12 +738,13 @@ subroutine set_stars_interactive(star,ieos,relax,nstar) do i=1,nstars print "(/,'------------- STAR ',i0,'-------------')",i - call set_star_interactive(star(i),ieos) + call set_star_interactive(star(i)) enddo ! prompt for equation of state and relaxation options if any stars made of gas if (nstars > 0) then if (any(star(1:nstars)%iprofile > 0)) then + if (any(star(1:nstars)%iprofile==ibpwpoly)) ieos = 9 ! set default eos for piecewise polytropes call set_star_eos_interactive(ieos,star) call prompt('Relax stars automatically during setup?',relax) endif @@ -759,11 +757,11 @@ end subroutine set_stars_interactive ! write setupfile options needed for a star !+ !----------------------------------------------------------------------- -subroutine write_options_star(star,ieos,iunit,label) +subroutine write_options_star(star,iunit,label) use infile_utils, only:write_inopt,get_optstring use setstar_utils, only:nprofile_opts,profile_opt,need_inputprofile,need_rstar,need_polyk type(star_t), intent(in) :: star - integer, intent(in) :: ieos,iunit + integer, intent(in) :: iunit character(len=*), intent(in), optional :: label character(len=120) :: string character(len=10) :: c @@ -773,8 +771,8 @@ subroutine write_options_star(star,ieos,iunit,label) if (present(label)) c = trim(adjustl(label)) write(iunit,"(/,a)") '# options for star '//trim(c) - call get_optstring(nprofile_opts,profile_opt,string,4) - call write_inopt(star%iprofile,'iprofile'//trim(c),'0=Sink,'//trim(string(1:40)),iunit) + call get_optstring(nprofile_opts+1,profile_opt,string,4,from_zero=.true.) + call write_inopt(star%iprofile,'iprofile'//trim(c),trim(string(1:48)),iunit) if (star%isoftcore <= 0) then if (need_inputprofile(star%iprofile)) then @@ -834,12 +832,6 @@ subroutine write_options_star(star,ieos,iunit,label) if (need_polyk(star%iprofile)) call write_inopt(star%polyk,'polyk'//trim(c),'polytropic constant (cs^2 if isothermal)',iunit) - ! options for setting initial thermal energy (e.g. if degenerate matter eos) - select case(ieos) - case(15) - call write_inopt(star%initialtemp,'initialtemp'//trim(c),'initial temperature of star (e.g. if degenerate matter eos)',iunit) - end select - if (star%iprofile > 0 .and. (len_trim(c)==0 .or. c(1:1)=='1')) then call write_inopt(star%np,'np'//trim(c),'number of particles',iunit) endif @@ -851,12 +843,11 @@ end subroutine write_options_star ! read setupfile options needed for a star !+ !----------------------------------------------------------------------- -subroutine read_options_star(star,ieos,db,nerr,label) +subroutine read_options_star(star,db,nerr,label) use infile_utils, only:inopts,read_inopt use setstar_utils, only:need_inputprofile,need_rstar,nprofile_opts type(star_t), intent(inout) :: star type(inopts), allocatable, intent(inout) :: db(:) - integer, intent(inout) :: ieos integer, intent(inout) :: nerr character(len=*), intent(in), optional :: label character(len=10) :: c @@ -870,8 +861,7 @@ subroutine read_options_star(star,ieos,db,nerr,label) star%dens_profile = 'relax'//trim(c)//'.profile' call read_inopt(star%iprofile,'iprofile'//trim(c),db,errcount=nerr,min=0,max=nprofile_opts) - call set_defaults_given_profile(star%iprofile,star%input_profile,& - ieos,star%m,star%polyk) + call set_defaults_given_profile(star%iprofile,star%input_profile,star%m,star%polyk) if (need_inputprofile(star%iprofile)) then call read_inopt(star%input_profile,'input_profile'//trim(c),db,errcount=nerr) @@ -922,12 +912,6 @@ subroutine read_options_star(star,ieos,db,nerr,label) if (need_polyk(star%iprofile)) call read_inopt(star%polyk,'polyk'//trim(c),db,errcount=nerr) - ! options for setting initial thermal energy (e.g. if degenerate matter eos) - select case(ieos) - case(15) - call read_inopt(star%initialtemp,'initialtemp'//trim(c),db,errcount=nerr,min=0.,max=1e12) - end select - ! star properties if (star%isoftcore <= 0) then if (need_inputprofile(star%iprofile)) then @@ -952,17 +936,18 @@ end subroutine read_options_star !----------------------------------------------------------------------- subroutine write_options_stars(star,relax,write_rho_to_file,ieos,iunit,nstar) use relaxstar, only:write_options_relax - use infile_utils, only:write_inopt + use infile_utils, only:write_inopt,int_to_string use apr, only:use_apr,write_options_apr type(star_t), intent(in) :: star(:) integer, intent(in) :: ieos,iunit logical, intent(in) :: relax,write_rho_to_file integer, intent(in), optional :: nstar integer :: i,nstars + character(len=3) :: label(size(star)) ! optionally ask for number of stars, otherwise fix nstars to the input array size if (present(nstar)) then - call write_inopt(nstar,'nstars','number of stars to add (0-'//achar(size(star)+48)//')',iunit) + call write_inopt(nstar,'nstars','number of stars to add (0-'//int_to_string(size(star))//')',iunit) nstars = nstar else nstars = size(star) @@ -970,18 +955,15 @@ subroutine write_options_stars(star,relax,write_rho_to_file,ieos,iunit,nstar) ! write options for each star do i=1,nstars - call write_options_star(star(i),ieos,iunit,label=achar(i+48)) + label(i) = trim(int_to_string(i)) + call write_options_star(star(i),iunit,label=label(i)) enddo - ! write equation of state options if any stars made of gas - if (nstars > 0) then - if (any(star(1:nstars)%iprofile > 0)) & - call write_options_stars_eos(star(1:nstars),ieos,iunit) - endif - - ! write relaxation options if any stars are made of gas + ! write equation of state and relaxation options if any stars made of gas if (nstars > 0) then if (any(star(1:nstars)%iprofile > 0)) then + call write_options_stars_eos(nstars,star(1:nstars),label(1:nstars),ieos,iunit) + write(iunit,"(/,a)") '# relaxation options' call write_inopt(relax,'relax','relax stars into equilibrium',iunit) call write_options_relax(iunit) @@ -1000,7 +982,7 @@ end subroutine write_options_stars !----------------------------------------------------------------------- subroutine read_options_stars(star,ieos,relax,write_rho_to_file,db,nerr,nstar) use relaxstar, only:read_options_relax - use infile_utils, only:inopts,read_inopt + use infile_utils, only:inopts,read_inopt,int_to_string use apr, only:use_apr,apr_max_in,ref_dir,apr_type,apr_rad,apr_drad type(star_t), intent(inout) :: star(:) ! inout because can set default options manually in calling routine type(inopts), allocatable, intent(inout) :: db(:) @@ -1009,6 +991,7 @@ subroutine read_options_stars(star,ieos,relax,write_rho_to_file,db,nerr,nstar) integer, intent(inout) :: nerr integer, intent(out), optional :: nstar integer :: i,nstars + character(len=3) :: label(size(star)) ! optionally ask for number of stars if (present(nstar)) then @@ -1018,19 +1001,19 @@ subroutine read_options_stars(star,ieos,relax,write_rho_to_file,db,nerr,nstar) nstars = size(star) endif - ! write equation of state options if any stars made of gas - if (nstars > 0) then - if (any(star(1:nstars)%iprofile > 0)) call read_options_stars_eos(star,ieos,db,nerr) - endif - ! read options for each star do i=1,nstars - call read_options_star(star(i),ieos,db,nerr,label=achar(i+48)) + label(i) = trim(int_to_string(i)) + call read_options_star(star(i),db,nerr,label=label(i)) enddo - ! read relaxation options if any stars are made of gas + ! equation of state and relaxation options if any stars made of gas if (nstars > 0) then if (any(star(1:nstars)%iprofile > 0)) then + if (any(star(1:nstars)%iprofile==ibpwpoly)) ieos = 9 ! set default eos for piecewise polytropes + ! equation of state options + call read_options_stars_eos(nstars,star(1:nstars),label(1:nstars),ieos,db,nerr) + ! relaxation options call read_inopt(relax,'relax',db,errcount=nerr) call read_options_relax(db,nerr) ! option to write density profile to file @@ -1053,11 +1036,13 @@ end subroutine read_options_stars ! write equation of state options needed to setup stars !+ !----------------------------------------------------------------------- -subroutine write_options_stars_eos(star,ieos,iunit) +subroutine write_options_stars_eos(nstars,star,label,ieos,iunit) use eos, only:use_var_comp,X_in,Z_in,irecomb,gmw,gamma use infile_utils, only:write_inopt - integer, intent(in) :: ieos,iunit - type(star_t), intent(in) :: star(:) + integer, intent(in) :: nstars,ieos,iunit + type(star_t), intent(in) :: star(nstars) + character(len=*), intent(in) :: label(nstars) + integer :: i write(iunit,"(/,a)") '# equation of state used to set the thermal energy profile' call write_inopt(ieos,'ieos','1=isothermal,2=adiabatic,10=MESA,12=idealplusrad',iunit) @@ -1079,6 +1064,14 @@ subroutine write_options_stars_eos(star,ieos,iunit) call write_inopt(X_in,'X','hydrogen mass fraction',iunit) call write_inopt(Z_in,'Z','metallicity',iunit) endif + case(15) + ! options for setting initial thermal energy (e.g. if degenerate matter eos) + do i=1,nstars + if (star(i)%iprofile > 0) then + call write_inopt(star(i)%initialtemp,'initialtemp'//trim(label(i)),& + 'initial temperature of star (e.g. if degenerate matter eos)',iunit) + endif + enddo end select end subroutine write_options_stars_eos @@ -1088,13 +1081,16 @@ end subroutine write_options_stars_eos ! read equation of state options needed to setup stars !+ !----------------------------------------------------------------------- -subroutine read_options_stars_eos(star,ieos,db,nerr) +subroutine read_options_stars_eos(nstars,star,label,ieos,db,nerr) use eos, only:use_var_comp,X_in,Z_in,irecomb,gamma,gmw use infile_utils, only:inopts,read_inopt - type(star_t), intent(out) :: star(:) + integer, intent(in) :: nstars + type(star_t), intent(inout) :: star(nstars) + character(len=*), intent(in) :: label(nstars) type(inopts), allocatable, intent(inout) :: db(:) integer, intent(inout) :: ieos integer, intent(inout) :: nerr + integer :: i ! equation of state call read_inopt(ieos,'ieos',db,errcount=nerr) @@ -1105,14 +1101,23 @@ subroutine read_options_stars_eos(star,ieos,db,nerr) call read_inopt(EOSopt,'EOSopt',db,min=0,max=4,errcount=nerr) case(2,12) call read_inopt(gamma,'gamma',db,min=1.,max=7.,errcount=nerr) - if ( (.not. use_var_comp) .and. any(need_mu(star(:)%isoftcore))) call read_inopt(gmw,'mu',db,min=0.,errcount=nerr) + if ((.not. use_var_comp) .and. any(need_mu(star(:)%isoftcore))) then + call read_inopt(gmw,'mu',db,min=0.,errcount=nerr) + endif case(10,20) if (ieos==20) call read_inopt(irecomb,'irecomb',db,errcount=nerr) ! if softening stellar core, composition is automatically determined at R/2 - if ( (.not. use_var_comp) .and. any(need_mu(star(:)%isoftcore))) then + if ((.not. use_var_comp) .and. any(need_mu(star(:)%isoftcore))) then call read_inopt(X_in,'X',db,min=0.,max=1.,errcount=nerr) call read_inopt(Z_in,'Z',db,min=0.,max=1.,errcount=nerr) endif + case(15) + do i=1,nstars + if (star(i)%iprofile > 0) then + call read_inopt(star(i)%initialtemp,'initialtemp'//trim(label(i)),& + db,min=0.,max=1e12,errcount=nerr) + endif + enddo end select end subroutine read_options_stars_eos diff --git a/src/setup/set_star_utils.f90 b/src/setup/set_star_utils.f90 index 449c4eecc..074d0473d 100644 --- a/src/setup/set_star_utils.f90 +++ b/src/setup/set_star_utils.f90 @@ -26,6 +26,7 @@ module setstar_utils ! Index of setup options ! integer, parameter, public :: nprofile_opts = 7 ! maximum number of initial configurations + integer, parameter, public :: ipointmass = 0 integer, parameter, public :: iuniform = 1 integer, parameter, public :: ipoly = 2 integer, parameter, public :: ifromfile = 3 @@ -34,8 +35,9 @@ module setstar_utils integer, parameter, public :: ibpwpoly = 6 integer, parameter, public :: ievrard = 7 - character(len=*), parameter, public :: profile_opt(nprofile_opts) = & - (/'Uniform density profile ', & + character(len=*), parameter, public :: profile_opt(0:nprofile_opts) = & + (/'Sink particle/point mass ', & + 'Uniform density sphere ', & 'Polytrope ', & 'Density vs r from ascii file', & 'KEPLER star from file ', & From 32b4b15b6927d926087426efd4465828d8569e96 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 9 Dec 2024 10:53:58 +1100 Subject: [PATCH 21/36] (build) buildbot now checks that all iprofile options work in SETUP=star; issues with piecewise polytrope fixed --- scripts/buildbot.sh | 42 +++++++++++++++++++++++++++--------------- src/setup/set_star.f90 | 5 ++++- 2 files changed, 31 insertions(+), 16 deletions(-) diff --git a/scripts/buildbot.sh b/scripts/buildbot.sh index 09a1c8e19..2d72f04b2 100755 --- a/scripts/buildbot.sh +++ b/scripts/buildbot.sh @@ -72,11 +72,12 @@ echo "url = $url"; pwd=$PWD; phantomdir="$pwd/../"; listofcomponents='main setup analysis utils'; -#listofcomponents='setup' +listofcomponents='setup' # # get list of targets, components and setups to check # allsetups=`grep 'ifeq ($(SETUP)' $phantomdir/build/Makefile_setups | grep -v skip | cut -d, -f 2 | cut -d')' -f 1` +allsetups='star' setuparr=($allsetups) batchsize=$(( ${#setuparr[@]} / $nbatch + 1 )) offset=$(( ($batch-1) * $batchsize )) @@ -197,11 +198,11 @@ check_phantomsetup () # run phantomsetup up to 3 times to successfully create/rewrite the .setup file # infile="${prefix}.in" - ./phantomsetup $prefix $flags < myinput.txt > /dev/null; - ./phantomsetup $prefix $flags < myinput.txt > /dev/null; + ./phantomsetup $prefix $flags < mycleanin.txt > /dev/null; + ./phantomsetup $prefix $flags < mycleanin.txt > /dev/null; if [ -e "$prefix.setup" ]; then print_result "creates .setup file" $pass; - #test_setupfile_options "$prefix" "$prefix.setup" $infile; + test_setupfile_options "$prefix.setup" "$flags" "$setup" $infile; else print_result "no .setup file" $warn; fi @@ -255,22 +256,33 @@ check_phantomsetup () test_setupfile_options() { myfail=0; - setup=$1; - setupfile=$2; - infile=$3; + #"$prefix.setup" "$flags" "$setup" $infile + setupfile=$1; + flags=$2; + setup=$3; + infile=$4; range='' - if [ "X$setup"=="Xstar" ]; then - param='iprofile' - range='1 2 3 4 5 6 7' + if [ "X$setup" == "Xstar" ]; then + param='iprofile1' + range='0 1 2 3 4 5 6 7' fi for x in $range; do valstring="$param = $x" - echo "checking $valstring" + echo "checking $valstring for SETUP=$setup" + rm $setupfile; + ./phantomsetup $setupfile $flags < mycleanin.txt > /dev/null; sed "s/$param.*=.*$/$valstring/" $setupfile > ${setupfile}.tmp - cp ${setupfile}.tmp $setupfile + mv ${setupfile}.tmp $setupfile + if [ "X$x" == "X6" ]; then + #sed "s/ieos.*=.*$/ieos = 9/" $setupfile > ${setupfile}.tmp + #mv ${setupfile}.tmp $setupfile + sed "s/dist_unit.*=.*$/dist_unit = km/" $setupfile > ${setupfile}.tmp + mv ${setupfile}.tmp $setupfile + fi + echo $setupfile rm $infile - ./phantomsetup $setupfile < /dev/null > /dev/null; - ./phantomsetup $setupfile < /dev/null; + ./phantomsetup $setupfile $flags < /dev/null > /dev/null; + ./phantomsetup $setupfile $flags < /dev/null > /dev/null; if [ -e $infile ]; then print_result "successful phantomsetup with $valstring" $pass; @@ -424,7 +436,7 @@ for setup in $listofsetups; do echo $setup >> $faillog; fi if [ -e $errorlogold ]; then - diff --unchanged-line-format="" --old-line-format="" --new-line-format="%L" $errorlogold $errorlog | tail -20 > warnings.tmp + diff $errorlogold $errorlog | tail -20 > warnings.tmp if [ -s warnings.tmp ]; then newwarn=1; else diff --git a/src/setup/set_star.f90 b/src/setup/set_star.f90 index a514bb1c1..ed726506d 100644 --- a/src/setup/set_star.f90 +++ b/src/setup/set_star.f90 @@ -438,7 +438,10 @@ subroutine set_stars(id,master,nstars,star,xyzh,vxyzu,eos_vars,rad,& procedure(mask_prototype) :: mask integer :: i - if (ieos==9) call init_eos_piecewise_preset(EOSopt) + ! initialise piecewise polytropic equation of state if piecewise polytrope used + if (ieos==9 .or. any(star(:)%iprofile==ibpwpoly)) call init_eos_piecewise_preset(EOSopt) + if (any(star(:)%iprofile==ibpwpoly)) call init_eos(9,ierr) + call init_eos(ieos,ierr) if (ierr /= 0) then call error('setup','could not initialise equation of state') From c2b0c670963cbcede6e32183997a18683d495943 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 9 Dec 2024 11:15:21 +1100 Subject: [PATCH 22/36] (grtde) build failures fixed --- src/setup/set_star.f90 | 3 +++ src/setup/setup_grtde.f90 | 12 +++++------- src/setup/setup_star.f90 | 2 +- 3 files changed, 9 insertions(+), 8 deletions(-) diff --git a/src/setup/set_star.f90 b/src/setup/set_star.f90 index ed726506d..f9b410d9b 100644 --- a/src/setup/set_star.f90 +++ b/src/setup/set_star.f90 @@ -733,6 +733,7 @@ subroutine set_stars_interactive(star,ieos,relax,nstar) ! optionally ask for number of stars, otherwise fix nstars to the input array size if (present(nstar) .and. size(star) > 1) then + nstar = 1 call prompt('how many stars to set up (0-'//int_to_string(size(star))//')',nstar,0,size(star)) nstars = nstar else @@ -749,6 +750,8 @@ subroutine set_stars_interactive(star,ieos,relax,nstar) if (any(star(1:nstars)%iprofile > 0)) then if (any(star(1:nstars)%iprofile==ibpwpoly)) ieos = 9 ! set default eos for piecewise polytropes call set_star_eos_interactive(ieos,star) + + relax = .false. call prompt('Relax stars automatically during setup?',relax) endif endif diff --git a/src/setup/setup_grtde.f90 b/src/setup/setup_grtde.f90 index bb25f7e19..1fb776d80 100644 --- a/src/setup/setup_grtde.f90 +++ b/src/setup/setup_grtde.f90 @@ -113,7 +113,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, norbits = 5. dumpsperorbit = 100 theta = 0. - write_profile = .false. + write_profile = .true. use_var_comp = .false. relax = .true. ! @@ -122,7 +122,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, if (id==master) print "(/,65('-'),1(/,a),/,65('-'),/)",' Tidal disruption in GR' filename = trim(fileprefix)//'.setup' inquire(file=filename,exist=iexist) - if (iexist) call read_setupfile(filename,ieos,ierr) + if (iexist) call read_setupfile(filename,ierr) if (.not. iexist .or. ierr /= 0) then if (id==master) then call write_setupfile(filename) @@ -250,14 +250,13 @@ subroutine write_setupfile(filename) use infile_utils, only:write_inopt use setstar, only:write_options_star use relaxstar, only:write_options_relax - use eos, only:ieos character(len=*), intent(in) :: filename integer :: iunit print "(a)",' writing setup options file '//trim(filename) open(newunit=iunit,file=filename,status='replace',form='formatted') write(iunit,"(a)") '# input file for tidal disruption setup' - call write_options_star(star,ieos,iunit) + call write_options_star(star,iunit) call write_inopt(relax,'relax','relax star into hydrostatic equilibrium',iunit) if (relax) call write_options_relax(iunit) @@ -272,7 +271,7 @@ subroutine write_setupfile(filename) end subroutine write_setupfile -subroutine read_setupfile(filename,ieos,ierr) +subroutine read_setupfile(filename,ierr) use infile_utils, only:open_db_from_file,inopts,read_inopt,close_db use io, only:error use setstar, only:read_options_star @@ -280,7 +279,6 @@ subroutine read_setupfile(filename,ieos,ierr) use physcon, only:solarm,solarr use units, only:set_units character(len=*), intent(in) :: filename - integer, intent(inout) :: ieos integer, intent(out) :: ierr integer, parameter :: iunit = 21 integer :: nerr @@ -298,7 +296,7 @@ subroutine read_setupfile(filename,ieos,ierr) ! !--read star options and convert to code units ! - call read_options_star(star,ieos,db,nerr) + call read_options_star(star,db,nerr) call read_inopt(relax,'relax',db,errcount=nerr) if (relax) call read_options_relax(db,nerr) diff --git a/src/setup/setup_star.f90 b/src/setup/setup_star.f90 index 288b8b119..67cdc057f 100644 --- a/src/setup/setup_star.f90 +++ b/src/setup/setup_star.f90 @@ -127,8 +127,8 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, ! if (.not.gravity) then iexternalforce = iext_densprofile - write_rho_to_file = .true. endif + write_rho_to_file = .true. ! ! set up particles ! From 304ab17de528906a687ad162f8d58b8dea483c1b Mon Sep 17 00:00:00 2001 From: ryosuke-hirai Date: Mon, 9 Dec 2024 17:54:57 +1100 Subject: [PATCH 23/36] (readwrite_mesa) Add a few more options for column names --- src/setup/readwrite_mesa.f90 | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/src/setup/readwrite_mesa.f90 b/src/setup/readwrite_mesa.f90 index d69a8ff72..e72c313e8 100644 --- a/src/setup/readwrite_mesa.f90 +++ b/src/setup/readwrite_mesa.f90 @@ -101,8 +101,9 @@ subroutine read_mesa(filepath,rho,r,pres,m,ene,temp,X_in,Z_in,Xfrac,Yfrac,Mstar, call read_column_labels(iu,nheaderlines,ncols,nlabels,header) if (nlabels /= ncols) print*,' WARNING: different number of labels compared to columns' - allocate(m(lines),r(lines),pres(lines),rho(lines),ene(lines), & - temp(lines),Xfrac(lines),Yfrac(lines)) + allocate(m(lines)) + m = -1d0 + allocate(r,pres,rho,ene,temp,Xfrac,Yfrac,source=m) over_directions: do idir=1,2 ! try backwards, then forwards if (idir==1) then @@ -152,12 +153,18 @@ subroutine read_mesa(filepath,rho,r,pres,m,ene,temp,X_in,Z_in,Xfrac,Yfrac,Mstar, r = (10**dat(1:lines,i)) * solarr case('pressure','p') pres = dat(1:lines,i) + case('logP') + pres = 10**dat(1:lines,i) case('temperature','t') temp = dat(1:lines,i) - case('x_mass_fraction_h','xfrac') + case('x_mass_fraction_h','x_mass_fraction_H','x','xfrac','h1') Xfrac = dat(1:lines,i) - case('y_mass_fraction_he','yfrac') + case('log_x') + Xfrac = 10**dat(1:lines,i) + case('y_mass_fraction_he','y_mass_fraction_He','y','yfrac','he4') Yfrac = dat(1:lines,i) + case('log_y') + Yfrac = 10**dat(1:lines,i) case default got_column = .false. end select @@ -176,6 +183,13 @@ subroutine read_mesa(filepath,rho,r,pres,m,ene,temp,X_in,Z_in,Xfrac,Yfrac,Mstar, enddo over_directions close(iu) + if(min(minval(m),minval(r),minval(pres),minval(rho),minval(ene))<0d0)ierr = 1 + + if (ierr /= 0) then + print "(a,/)",' ERROR reading MESA file [missing required columns]' + return + endif + if (.not. usecgs) then m = m / umass r = r / udist From 9e0a59ab56fa7beb72ffae465a5936eb88752fdf Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 9 Dec 2024 18:09:01 +1100 Subject: [PATCH 24/36] BUG FIX: fixes #604 where star would fail to relax if Mstar /= 1; incorrect divide-by-mstar in relaxation routine --- src/setup/relax_star.f90 | 2 +- src/tests/test_setstar.f90 | 10 ++++++---- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/setup/relax_star.f90 b/src/setup/relax_star.f90 index c05c2deaa..b8acc853b 100644 --- a/src/setup/relax_star.f90 +++ b/src/setup/relax_star.f90 @@ -421,7 +421,7 @@ subroutine reset_u_and_get_errors(i1,npart,xyzh,vxyzu,rad,nt,mr,rho,& do i = i1+1,npart ri = sqrt(dot_product(xyzh(1:3,i),xyzh(1:3,i))) - massri = mass_enclosed_r(i-i1)/mstar + massri = mass_enclosed_r(i-i1) rhor = yinterp(rho,mr,massri) ! analytic rho(r) if (use_apr) then diff --git a/src/tests/test_setstar.f90 b/src/tests/test_setstar.f90 index 43f6b3468..a847e9e2e 100644 --- a/src/tests/test_setstar.f90 +++ b/src/tests/test_setstar.f90 @@ -145,10 +145,12 @@ subroutine test_polytrope(ntests,npass) npartoftype = 0 massoftype = 0. iverbose = 0 - call set_units(dist=solarr,mass=solarm,G=1.d0) + call set_units(dist=10.*solarr,mass=10.*solarm,G=1.d0) ieos = 2 gamma = 5./3. call set_defaults_star(star) + star%m = '1*msun' + star%r = '1*rsun' star%iprofile = ipoly ! a polytrope star%np = 1000 x0 = 0. @@ -170,13 +172,13 @@ subroutine test_polytrope(ntests,npass) call checkval(nerror+nwarn,0,0,nfail(1),'no errors or warnings') call update_test_scores(ntests,nfail,npass) - call checkval(rhozero,1./(4./3.*pi),1e-6,nfail(1),'mean density') + call checkval(rhozero,100./(4./3.*pi),1e-6,nfail(1),'mean density') call update_test_scores(ntests,nfail,npass) - call checkval(star%polyk,0.424304,1e-6,nfail(1),'polyk value for M=1,R=1') + call checkval(star%polyk,1.9694457e-2,1e-6,nfail(1),'polyk value for M=0.1,R=0.1') call update_test_scores(ntests,nfail,npass) - call checkval(polyk,0.424304,1e-6,nfail(1),'polyk value for M=1,R=1') + call checkval(polyk,1.9694457e-2,1e-6,nfail(1),'polyk value for M=0.1,R=0.1') call update_test_scores(ntests,nfail,npass) call checkval(rmserr,0.0,0.04,nfail(1),'error in density profile') From 121b48f015129ff481357ac61860173e469b5d97 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 9 Dec 2024 18:15:34 +1100 Subject: [PATCH 25/36] [header-bot] updated file headers --- src/main/extern_densprofile.f90 | 2 +- src/main/relaxem.f90 | 2 +- src/setup/relax_star.f90 | 3 +-- src/setup/set_star.f90 | 20 +++++++++++++++----- src/setup/set_star_utils.f90 | 4 ++-- src/setup/setup_binary.f90 | 2 +- src/setup/setup_star.f90 | 21 ++++----------------- src/tests/test_apr.f90 | 4 ++-- src/tests/test_setstar.f90 | 4 +++- src/tests/testsuite.F90 | 2 +- 10 files changed, 31 insertions(+), 33 deletions(-) diff --git a/src/main/extern_densprofile.f90 b/src/main/extern_densprofile.f90 index 2a78573a6..b759a311a 100644 --- a/src/main/extern_densprofile.f90 +++ b/src/main/extern_densprofile.f90 @@ -16,7 +16,7 @@ module extern_densprofile ! ! :Runtime parameters: None ! -! :Dependencies: datafiles, io, physcon, units +! :Dependencies: datafiles, io, physcon, table_utils, units ! implicit none diff --git a/src/main/relaxem.f90 b/src/main/relaxem.f90 index 417896bf1..f7cfe238f 100644 --- a/src/main/relaxem.f90 +++ b/src/main/relaxem.f90 @@ -14,7 +14,7 @@ module relaxem ! ! :Runtime parameters: None ! -! :Dependencies: boundary, deriv, dim, eos, kernel, mpidomain, options, +! :Dependencies: boundary, deriv, dim, eos, io, kernel, mpidomain, options, ! part ! implicit none diff --git a/src/setup/relax_star.f90 b/src/setup/relax_star.f90 index b8acc853b..3512837cc 100644 --- a/src/setup/relax_star.f90 +++ b/src/setup/relax_star.f90 @@ -20,8 +20,7 @@ module relaxstar ! :Dependencies: apr, checksetup, damping, deriv, dim, dump_utils, ! energies, eos, externalforces, fileutils, infile_utils, initial, io, ! io_summary, linklist, memory, options, part, physcon, ptmass, -! readwrite_dumps, setstar_utils, sortutils, step_lf_global, table_utils, -! units +! readwrite_dumps, setstar_utils, step_lf_global, table_utils, units ! implicit none public :: relax_star,write_options_relax,read_options_relax diff --git a/src/setup/set_star.f90 b/src/setup/set_star.f90 index f9b410d9b..b44e124f9 100644 --- a/src/setup/set_star.f90 +++ b/src/setup/set_star.f90 @@ -16,12 +16,22 @@ module setstar ! :Owner: Daniel Price ! ! :Runtime parameters: -! - nstars : *number of stars to add (0-'//achar(size(star)+48)//')* -! - relax : *relax stars into equilibrium* +! - EOSopt : *EOS: 1=APR3,2=SLy,3=MS1,4=ENG (from Read et al 2009)* +! - X : *hydrogen mass fraction* +! - gamma : *Adiabatic index* +! - ieos : *1=isothermal,2=adiabatic,10=MESA,12=idealplusrad* +! - irecomb : *Species to include in recombination (0: H2+H+He, 1:H+He, 2:He* +! - metallicity : *metallicity* +! - mu : *mean molecular weight* +! - nstars : *number of stars to add (0-'//int_to_string(size(star))//')* +! - relax : *relax stars into equilibrium* +! - use_var_comp : *Use variable composition (X, Z, mu)* +! - write_rho_to_file : *write density profile(s) to file* ! -! :Dependencies: centreofmass, dim, eos, extern_densprofile, infile_utils, -! io, mpiutils, part, physcon, prompting, radiation_utils, relaxstar, -! setstar_utils, unifdis, units, vectorutils +! :Dependencies: apr, centreofmass, dim, eos, eos_piecewise, +! extern_densprofile, infile_utils, io, mpiutils, part, physcon, +! prompting, radiation_utils, relaxstar, setstar_utils, unifdis, units, +! vectorutils ! use setstar_utils, only:ikepler,imesa,ibpwpoly,ipoly,iuniform,ifromfile,ievrard,& need_polyk,need_mu diff --git a/src/setup/set_star_utils.f90 b/src/setup/set_star_utils.f90 index 074d0473d..7ac2f2e18 100644 --- a/src/setup/set_star_utils.f90 +++ b/src/setup/set_star_utils.f90 @@ -14,8 +14,8 @@ module setstar_utils ! ! :Runtime parameters: None ! -! :Dependencies: eos, eos_piecewise, extern_densprofile, io, kernel, part, -! physcon, radiation_utils, readwrite_kepler, readwrite_mesa, +! :Dependencies: dim, eos, eos_piecewise, extern_densprofile, io, kernel, +! part, physcon, radiation_utils, readwrite_kepler, readwrite_mesa, ! rho_profile, setsoftenedcore, sortutils, spherical, table_utils, ! unifdis, units ! diff --git a/src/setup/setup_binary.f90 b/src/setup/setup_binary.f90 index 5bebb3c3f..37d450fbf 100644 --- a/src/setup/setup_binary.f90 +++ b/src/setup/setup_binary.f90 @@ -17,7 +17,7 @@ module setup ! ! :Dependencies: centreofmass, dim, eos, externalforces, infile_utils, io, ! kernel, mpidomain, options, part, physcon, setorbit, setstar, setunits, -! setup_params +! setup_params, units ! use setstar, only:star_t use setorbit, only:orbit_t diff --git a/src/setup/setup_star.f90 b/src/setup/setup_star.f90 index 67cdc057f..7fd11ef99 100644 --- a/src/setup/setup_star.f90 +++ b/src/setup/setup_star.f90 @@ -12,24 +12,11 @@ module setup ! ! :Owner: Daniel Price ! -! :Runtime parameters: -! - EOSopt : *EOS: 1=APR3,2=SLy,3=MS1,4=ENG (from Read et al 2009)* -! - X : *hydrogen mass fraction* -! - gamma : *Adiabatic index* -! - ieos : *1=isothermal,2=adiabatic,10=MESA,12=idealplusrad* -! - initialtemp : *initial temperature of the star* -! - irecomb : *Species to include in recombination (0: H2+H+He, 1:H+He, 2:He* -! - metallicity : *metallicity* -! - mu : *mean molecular weight* -! - polyk : *polytropic constant (cs^2 if isothermal)* -! - relax_star : *relax star(s) automatically during setup* -! - use_var_comp : *Use variable composition (X, Z, mu)* -! - write_rho_to_file : *write density profile(s) to file* +! :Runtime parameters: None ! -! :Dependencies: apr, dim, eos, eos_gasradrec, eos_piecewise, -! extern_densprofile, externalforces, infile_utils, io, kernel, -! mpidomain, mpiutils, options, part, physcon, prompting, relaxstar, -! setstar, setunits, setup_params, timestep, units +! :Dependencies: dim, eos, externalforces, infile_utils, io, kernel, +! mpidomain, options, part, physcon, setstar, setunits, setup_params, +! timestep ! use io, only:fatal,error,warning,master use part, only:gravity,gr diff --git a/src/tests/test_apr.f90 b/src/tests/test_apr.f90 index 002768959..eebb5c0f5 100644 --- a/src/tests/test_apr.f90 +++ b/src/tests/test_apr.f90 @@ -14,8 +14,8 @@ module testapr ! ! :Runtime parameters: None ! -! :Dependencies: apr, boundary, dim, io, linklist, mpidomain, mpiutils, -! part, physcon, testutils, unifdis, units +! :Dependencies: apr, boundary, dim, io, mpidomain, mpiutils, part, +! testutils, unifdis ! use testutils, only:checkval,update_test_scores use io, only:id,master,fatal diff --git a/src/tests/test_setstar.f90 b/src/tests/test_setstar.f90 index a847e9e2e..a68ffb84a 100644 --- a/src/tests/test_setstar.f90 +++ b/src/tests/test_setstar.f90 @@ -14,7 +14,9 @@ module testsetstar ! ! :Runtime parameters: None ! -! :Dependencies: +! :Dependencies: checksetup, dim, eos, io, mpidomain, options, part, +! physcon, setstar, setstar_utils, sortutils, table_utils, testutils, +! units ! use testutils, only:checkval,update_test_scores implicit none diff --git a/src/tests/testsuite.F90 b/src/tests/testsuite.F90 index e66c73210..86db01be7 100644 --- a/src/tests/testsuite.F90 +++ b/src/tests/testsuite.F90 @@ -21,7 +21,7 @@ module test ! testgrowth, testindtstep, testiorig, testkdtree, testkernel, testlink, ! testmath, testmpi, testnimhd, testpart, testpoly, testptmass, ! testradiation, testrwdump, testsedov, testsetdisc, testsethier, -! testsmol, teststep, testwind, timing +! testsetstar, testsmol, teststep, testwind, timing ! implicit none public :: testsuite From 1b9a59b1e6c53b697eca6f82bbf06e6abda1f897 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 9 Dec 2024 18:15:52 +1100 Subject: [PATCH 26/36] [space-bot] whitespace at end of lines removed --- src/setup/set_star.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/setup/set_star.f90 b/src/setup/set_star.f90 index b44e124f9..a73b411ed 100644 --- a/src/setup/set_star.f90 +++ b/src/setup/set_star.f90 @@ -754,7 +754,7 @@ subroutine set_stars_interactive(star,ieos,relax,nstar) print "(/,'------------- STAR ',i0,'-------------')",i call set_star_interactive(star(i)) enddo - + ! prompt for equation of state and relaxation options if any stars made of gas if (nstars > 0) then if (any(star(1:nstars)%iprofile > 0)) then @@ -765,7 +765,7 @@ subroutine set_stars_interactive(star,ieos,relax,nstar) call prompt('Relax stars automatically during setup?',relax) endif endif - + end subroutine set_stars_interactive !----------------------------------------------------------------------- @@ -1031,9 +1031,9 @@ subroutine read_options_stars(star,ieos,relax,write_rho_to_file,db,nerr,nstar) call read_options_stars_eos(nstars,star(1:nstars),label(1:nstars),ieos,db,nerr) ! relaxation options call read_inopt(relax,'relax',db,errcount=nerr) - call read_options_relax(db,nerr) + call read_options_relax(db,nerr) ! option to write density profile to file - call read_inopt(write_rho_to_file,'write_rho_to_file',db,errcount=nerr) + call read_inopt(write_rho_to_file,'write_rho_to_file',db,errcount=nerr) endif endif From 3c9315e2f3aac764f38bd6e625173592d0add81c Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 9 Dec 2024 18:16:34 +1100 Subject: [PATCH 27/36] [indent-bot] standardised indentation --- src/setup/set_star.f90 | 24 ++++++++++++------------ src/setup/setup_grdisc.F90 | 38 +++++++++++++++++++------------------- 2 files changed, 31 insertions(+), 31 deletions(-) diff --git a/src/setup/set_star.f90 b/src/setup/set_star.f90 index a73b411ed..52a412be1 100644 --- a/src/setup/set_star.f90 +++ b/src/setup/set_star.f90 @@ -1038,11 +1038,11 @@ subroutine read_options_stars(star,ieos,relax,write_rho_to_file,db,nerr,nstar) endif if (use_apr) then - call read_inopt(apr_max_in,'apr_max',db,errcount=nerr) - call read_inopt(ref_dir,'ref_dir',db,errcount=nerr) - call read_inopt(apr_type,'apr_type',db,errcount=nerr) - call read_inopt(apr_rad,'apr_rad',db,errcount=nerr) - call read_inopt(apr_drad,'apr_drad',db,errcount=nerr) + call read_inopt(apr_max_in,'apr_max',db,errcount=nerr) + call read_inopt(ref_dir,'ref_dir',db,errcount=nerr) + call read_inopt(apr_type,'apr_type',db,errcount=nerr) + call read_inopt(apr_rad,'apr_rad',db,errcount=nerr) + call read_inopt(apr_drad,'apr_drad',db,errcount=nerr) endif end subroutine read_options_stars @@ -1081,13 +1081,13 @@ subroutine write_options_stars_eos(nstars,star,label,ieos,iunit) call write_inopt(Z_in,'Z','metallicity',iunit) endif case(15) - ! options for setting initial thermal energy (e.g. if degenerate matter eos) - do i=1,nstars - if (star(i)%iprofile > 0) then - call write_inopt(star(i)%initialtemp,'initialtemp'//trim(label(i)),& + ! options for setting initial thermal energy (e.g. if degenerate matter eos) + do i=1,nstars + if (star(i)%iprofile > 0) then + call write_inopt(star(i)%initialtemp,'initialtemp'//trim(label(i)),& 'initial temperature of star (e.g. if degenerate matter eos)',iunit) - endif - enddo + endif + enddo end select end subroutine write_options_stars_eos @@ -1167,7 +1167,7 @@ subroutine set_star_eos_interactive(ieos,star) if ( (.not. use_var_comp) .and. any(need_mu(star(:)%isoftcore))) then call prompt('Enter hydrogen mass fraction (X)',X_in,0.,1.) call prompt('Enter metals mass fraction (Z)',Z_in,0.,1.) - endif + endif end select end subroutine set_star_eos_interactive diff --git a/src/setup/setup_grdisc.F90 b/src/setup/setup_grdisc.F90 index 00ef2c139..605d9daa3 100644 --- a/src/setup/setup_grdisc.F90 +++ b/src/setup/setup_grdisc.F90 @@ -169,31 +169,31 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, ! either sink particles or balls of gas ! if (nstars > 0) then - write_profile = .false. - iexternalforce = 0 - call set_stars(id,master,nstars,star,xyzh,vxyzu,eos_vars,rad,npart,npartoftype,& + write_profile = .false. + iexternalforce = 0 + call set_stars(id,master,nstars,star,xyzh,vxyzu,eos_vars,rad,npart,npartoftype,& massoftype,hfact,xyzmh_ptmass,vxyz_ptmass,nptmass,ieos,gamma,& X_in,Z_in,relax,use_var_comp,write_profile,& rhozero,npart_total,i_belong,ierr) - do i=1,nstars - nptmass_in = 0 - ! convert stellar mass and radius to code units - mstar = in_code_units(star(i)%m,ierr) - rstar = in_code_units(star(i)%r,ierr) - call set_orbit(orbit(i),mhole/umass,mstar,r_in,rstar, & + do i=1,nstars + nptmass_in = 0 + ! convert stellar mass and radius to code units + mstar = in_code_units(star(i)%m,ierr) + rstar = in_code_units(star(i)%r,ierr) + call set_orbit(orbit(i),mhole/umass,mstar,r_in,rstar, & xyzmh_ptmass_in,vxyz_ptmass_in,nptmass_in,(id==master),ierr) - ! shift the star to the position of the second body - if (star(i)%iprofile > 0) then - call shift_star(npart,npartoftype,xyzh,vxyzu,& + ! shift the star to the position of the second body + if (star(i)%iprofile > 0) then + call shift_star(npart,npartoftype,xyzh,vxyzu,& x0=xyzmh_ptmass_in(:,2),v0=vxyz_ptmass_in(:,2),itype=i) - else - nptmass = nptmass + 1 - xyzmh_ptmass(:,nptmass) = xyzmh_ptmass_in(:,2) - vxyz_ptmass(:,nptmass) = vxyz_ptmass_in(:,2) - endif - enddo -endif + else + nptmass = nptmass + 1 + xyzmh_ptmass(:,nptmass) = xyzmh_ptmass_in(:,2) + vxyz_ptmass(:,nptmass) = vxyz_ptmass_in(:,2) + endif + enddo + endif #ifndef GR iexternalforce = iext_einsteinprec From 9c5df94dad81a2ad17e4360327e156295f74da7d Mon Sep 17 00:00:00 2001 From: ryosuke-hirai Date: Mon, 9 Dec 2024 18:18:13 +1100 Subject: [PATCH 28/36] (readwrite_mesa) Add a few more options for column names --- src/setup/readwrite_mesa.f90 | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/src/setup/readwrite_mesa.f90 b/src/setup/readwrite_mesa.f90 index e72c313e8..8e02843a5 100644 --- a/src/setup/readwrite_mesa.f90 +++ b/src/setup/readwrite_mesa.f90 @@ -138,25 +138,36 @@ subroutine read_mesa(filepath,rho,r,pres,m,ene,temp,X_in,Z_in,Xfrac,Yfrac,Mstar, case('mass','#mass','m') m = dat(1:lines,i) if (ismesafile .or. maxval(m) < 1.e-10*solarm) m = m * solarm ! If reading MESA profile, 'mass' is in units of Msun + case('logM','log_mass') + m = 10**dat(1:lines,i) + if (ismesafile .or. maxval(m) < 1.e-10*solarm) m = m * solarm ! If reading MESA profile, 'mass' is in units of Msun case('rho','density') rho = dat(1:lines,i) - case('logrho') + case('logrho','logRho') rho = 10**(dat(1:lines,i)) case('energy','e_int','e_internal') ene = dat(1:lines,i) + case('logE') + ene = 10**dat(1:lines,i) case('radius_cm') r = dat(1:lines,i) + case('radius_km') + r = dat(1:lines,i) * 1e5 case('radius','r') r = dat(1:lines,i) if (ismesafile .or. maxval(r) < 1e-10*solarr) r = r * solarr - case('logr') + case('logr','logR') r = (10**dat(1:lines,i)) * solarr + case('logR_cm') + r = 10**dat(1:lines,i) case('pressure','p') pres = dat(1:lines,i) case('logP') pres = 10**dat(1:lines,i) case('temperature','t') temp = dat(1:lines,i) + case('logT') + temp = 10**dat(1:lines,i) case('x_mass_fraction_h','x_mass_fraction_H','x','xfrac','h1') Xfrac = dat(1:lines,i) case('log_x') From 60af9b289c9930c00a61f2ed5a31362fc6dab76b Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 9 Dec 2024 20:45:07 +1100 Subject: [PATCH 29/36] (build) check all setups; remove debugging line --- scripts/buildbot.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/buildbot.sh b/scripts/buildbot.sh index 2d72f04b2..9b8b7df51 100755 --- a/scripts/buildbot.sh +++ b/scripts/buildbot.sh @@ -77,7 +77,7 @@ listofcomponents='setup' # get list of targets, components and setups to check # allsetups=`grep 'ifeq ($(SETUP)' $phantomdir/build/Makefile_setups | grep -v skip | cut -d, -f 2 | cut -d')' -f 1` -allsetups='star' +#allsetups='star' setuparr=($allsetups) batchsize=$(( ${#setuparr[@]} / $nbatch + 1 )) offset=$(( ($batch-1) * $batchsize )) From bab66c9c6c7ee9610a6a62cf37399efb688afcf5 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Mon, 9 Dec 2024 20:45:34 +1100 Subject: [PATCH 30/36] (build) check all setups+components; remove debugging line --- scripts/buildbot.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/buildbot.sh b/scripts/buildbot.sh index 9b8b7df51..3a51adc0e 100755 --- a/scripts/buildbot.sh +++ b/scripts/buildbot.sh @@ -72,7 +72,7 @@ echo "url = $url"; pwd=$PWD; phantomdir="$pwd/../"; listofcomponents='main setup analysis utils'; -listofcomponents='setup' +#listofcomponents='setup' # # get list of targets, components and setups to check # From 59824d339191c5bc0000aaab7b286fae04daa2da Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Tue, 10 Dec 2024 08:22:03 +1100 Subject: [PATCH 31/36] (blob) build failure in single precision fixed --- src/main/units.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/main/units.f90 b/src/main/units.f90 index 9aa59d101..a859c768d 100644 --- a/src/main/units.f90 +++ b/src/main/units.f90 @@ -307,7 +307,7 @@ end subroutine select_unit logical function is_time_unit(string) character(len=*), intent(in) :: string character(len=len_utype) :: unit_type - real :: val + real(kind=8) :: val integer :: ierr ierr = 0 @@ -325,7 +325,7 @@ end function is_time_unit logical function is_length_unit(string) character(len=*), intent(in) :: string character(len=len_utype) :: unit_type - real :: val + real(kind=8) :: val integer :: ierr ierr = 0 @@ -343,7 +343,7 @@ end function is_length_unit logical function is_mdot_unit(string) character(len=*), intent(in) :: string character(len=len_utype) :: unit_type - real :: val + real(kind=8) :: val integer :: ierr ierr = 0 @@ -580,7 +580,7 @@ real(kind=8) function in_units(val,unitstring) result(rval) character(len=*), intent(in) :: unitstring character(len=len_utype) :: utype integer :: ierr - real :: fac + real(kind=8) :: fac call select_unit(unitstring,fac,ierr,unit_type=utype) ! handle errors silently by ignoring ierr From 6af713194881de68cbe5eee247ce7dc19ddd7d20 Mon Sep 17 00:00:00 2001 From: Daniel Price Date: Tue, 10 Dec 2024 08:57:46 +1100 Subject: [PATCH 32/36] build failure with default filename in extern_densprofile fixed --- src/main/extern_densprofile.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/extern_densprofile.f90 b/src/main/extern_densprofile.f90 index b759a311a..c106668f8 100644 --- a/src/main/extern_densprofile.f90 +++ b/src/main/extern_densprofile.f90 @@ -25,7 +25,7 @@ module extern_densprofile integer, public :: ntab ! *** Add option to .in file to specify density profile / mass enclosed filename? *** - character(len=*), public, parameter :: rhotabfile = 'density-profile.tab' + character(len=*), public, parameter :: rhotabfile = 'relax1.profile' integer, public, parameter :: nrhotab = 5000 ! maximum allowed size of r rho tabulated arrays public :: densityprofile_force, load_extern_densityprofile, read_rhotab, write_rhotab, calc_menc From 559c3cec7594094bd3dc065f8d9d8e39edbd382b Mon Sep 17 00:00:00 2001 From: ryosuke-hirai Date: Tue, 10 Dec 2024 12:29:09 +1100 Subject: [PATCH 33/36] (readwrite_mesa) loosen error criterion --- src/setup/readwrite_mesa.f90 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/setup/readwrite_mesa.f90 b/src/setup/readwrite_mesa.f90 index 8e02843a5..22fdf8733 100644 --- a/src/setup/readwrite_mesa.f90 +++ b/src/setup/readwrite_mesa.f90 @@ -102,7 +102,7 @@ subroutine read_mesa(filepath,rho,r,pres,m,ene,temp,X_in,Z_in,Xfrac,Yfrac,Mstar, if (nlabels /= ncols) print*,' WARNING: different number of labels compared to columns' allocate(m(lines)) - m = -1d0 + m = -1. allocate(r,pres,rho,ene,temp,Xfrac,Yfrac,source=m) over_directions: do idir=1,2 ! try backwards, then forwards @@ -143,11 +143,11 @@ subroutine read_mesa(filepath,rho,r,pres,m,ene,temp,X_in,Z_in,Xfrac,Yfrac,Mstar, if (ismesafile .or. maxval(m) < 1.e-10*solarm) m = m * solarm ! If reading MESA profile, 'mass' is in units of Msun case('rho','density') rho = dat(1:lines,i) - case('logrho','logRho') + case('logrho') rho = 10**(dat(1:lines,i)) case('energy','e_int','e_internal') ene = dat(1:lines,i) - case('logE') + case('logr') ene = 10**dat(1:lines,i) case('radius_cm') r = dat(1:lines,i) @@ -156,23 +156,23 @@ subroutine read_mesa(filepath,rho,r,pres,m,ene,temp,X_in,Z_in,Xfrac,Yfrac,Mstar, case('radius','r') r = dat(1:lines,i) if (ismesafile .or. maxval(r) < 1e-10*solarr) r = r * solarr - case('logr','logR') + case('logr') r = (10**dat(1:lines,i)) * solarr - case('logR_cm') + case('logr_cm') r = 10**dat(1:lines,i) case('pressure','p') pres = dat(1:lines,i) - case('logP') + case('logp') pres = 10**dat(1:lines,i) case('temperature','t') temp = dat(1:lines,i) - case('logT') + case('logt') temp = 10**dat(1:lines,i) - case('x_mass_fraction_h','x_mass_fraction_H','x','xfrac','h1') + case('x_mass_fraction_h','x','xfrac','h1') Xfrac = dat(1:lines,i) case('log_x') Xfrac = 10**dat(1:lines,i) - case('y_mass_fraction_he','y_mass_fraction_He','y','yfrac','he4') + case('y_mass_fraction_he','y','yfrac','he4') Yfrac = dat(1:lines,i) case('log_y') Yfrac = 10**dat(1:lines,i) @@ -194,7 +194,7 @@ subroutine read_mesa(filepath,rho,r,pres,m,ene,temp,X_in,Z_in,Xfrac,Yfrac,Mstar, enddo over_directions close(iu) - if(min(minval(m),minval(r),minval(pres),minval(rho),minval(ene))<0d0)ierr = 1 + if(min(minval(m),minval(r),minval(pres),minval(rho))<0d0)ierr = 1 if (ierr /= 0) then print "(a,/)",' ERROR reading MESA file [missing required columns]' From dcbf82af608b16d699df4966a32efdb3f5efd5ea Mon Sep 17 00:00:00 2001 From: ryosuke-hirai Date: Tue, 10 Dec 2024 12:30:47 +1100 Subject: [PATCH 34/36] (readwrite_mesa) fix typo --- src/setup/readwrite_mesa.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/setup/readwrite_mesa.f90 b/src/setup/readwrite_mesa.f90 index 22fdf8733..78de7ffca 100644 --- a/src/setup/readwrite_mesa.f90 +++ b/src/setup/readwrite_mesa.f90 @@ -147,7 +147,7 @@ subroutine read_mesa(filepath,rho,r,pres,m,ene,temp,X_in,Z_in,Xfrac,Yfrac,Mstar, rho = 10**(dat(1:lines,i)) case('energy','e_int','e_internal') ene = dat(1:lines,i) - case('logr') + case('loge') ene = 10**dat(1:lines,i) case('radius_cm') r = dat(1:lines,i) From dd9251f81873d8fb88a329c3957754d77959d3b4 Mon Sep 17 00:00:00 2001 From: ryosuke-hirai Date: Tue, 10 Dec 2024 12:36:32 +1100 Subject: [PATCH 35/36] (readwrite_mesa) loosen error criterion --- src/setup/readwrite_mesa.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/setup/readwrite_mesa.f90 b/src/setup/readwrite_mesa.f90 index 78de7ffca..3329edef4 100644 --- a/src/setup/readwrite_mesa.f90 +++ b/src/setup/readwrite_mesa.f90 @@ -194,7 +194,7 @@ subroutine read_mesa(filepath,rho,r,pres,m,ene,temp,X_in,Z_in,Xfrac,Yfrac,Mstar, enddo over_directions close(iu) - if(min(minval(m),minval(r),minval(pres),minval(rho))<0d0)ierr = 1 + if(min(minval(r),minval(pres),minval(rho))<0d0)ierr = 1 if (ierr /= 0) then print "(a,/)",' ERROR reading MESA file [missing required columns]' From ae882611681812e7507388a8e85bc9b0faaffc1c Mon Sep 17 00:00:00 2001 From: ryosuke-hirai Date: Tue, 10 Dec 2024 12:37:44 +1100 Subject: [PATCH 36/36] (readwrite_mesa) loosen error criterion --- src/setup/readwrite_mesa.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/setup/readwrite_mesa.f90 b/src/setup/readwrite_mesa.f90 index 3329edef4..70bb69dbd 100644 --- a/src/setup/readwrite_mesa.f90 +++ b/src/setup/readwrite_mesa.f90 @@ -194,7 +194,7 @@ subroutine read_mesa(filepath,rho,r,pres,m,ene,temp,X_in,Z_in,Xfrac,Yfrac,Mstar, enddo over_directions close(iu) - if(min(minval(r),minval(pres),minval(rho))<0d0)ierr = 1 + if(min(minval(pres),minval(rho))<0d0)ierr = 1 if (ierr /= 0) then print "(a,/)",' ERROR reading MESA file [missing required columns]'