From 4303043a1a4e924bf00a3f09d796f4bbc2eee327 Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Wed, 7 Aug 2024 11:26:24 +1000 Subject: [PATCH 01/29] (tillotson) added module and subroutines for tillotson eos --- build/Makefile | 2 +- src/main/eos.f90 | 15 +++++++ src/main/eos_tillotson.f90 | 88 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 104 insertions(+), 1 deletion(-) create mode 100644 src/main/eos_tillotson.f90 diff --git a/build/Makefile b/build/Makefile index 219f2eb0b..3387bc040 100644 --- a/build/Makefile +++ b/build/Makefile @@ -509,7 +509,7 @@ SRCCHEM= fs_data.f90 mol_data.f90 utils_spline.f90 \ # equations of state # SRCMESA= eos_mesa_microphysics.f90 eos_mesa.f90 -SRCEOS = eos_barotropic.f90 eos_stratified.f90 eos_piecewise.f90 ${SRCMESA} eos_shen.f90 eos_helmholtz.f90 eos_idealplusrad.f90 ionization.F90 eos_gasradrec.f90 eos.f90 +SRCEOS = eos_tillotson.f90 eos_barotropic.f90 eos_stratified.f90 eos_piecewise.f90 ${SRCMESA} eos_shen.f90 eos_helmholtz.f90 eos_idealplusrad.f90 ionization.F90 eos_gasradrec.f90 eos.f90 ifeq ($(HDF5), yes) SRCREADWRITE_DUMPS= utils_hdf5.f90 utils_dumpfiles_hdf5.f90 readwrite_dumps_common.F90 readwrite_dumps_fortran.F90 readwrite_dumps_hdf5.F90 readwrite_dumps.F90 diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 45776d30a..bb886a207 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -115,6 +115,7 @@ subroutine equationofstate(eos_type,ponrhoi,spsoundi,rhoi,xi,yi,zi,tempi,eni,gam use eos_stratified, only:get_eos_stratified use eos_barotropic, only:get_eos_barotropic use eos_piecewise, only:get_eos_piecewise + use eos_tillotson, only:equationofstate_tillotson integer, intent(in) :: eos_type real, intent(in) :: rhoi,xi,yi,zi real, intent(out) :: ponrhoi,spsoundi @@ -423,6 +424,19 @@ subroutine equationofstate(eos_type,ponrhoi,spsoundi,rhoi,xi,yi,zi,tempi,eni,gam if (present(mu_local)) mu_local = 1./imui if (present(gamma_local)) gamma_local = gammai + case(21) + ! + !-- Tillotson EOS + ! + ! from Tillotson (1962) + ! + call equationofstate_tillotson(cgsrhoi,cgseni*cgsrhoi,temperaturei,mui,X_i,1.-X_i-Z_i,cgspresi,cgsspsoundi,gammai) + ponrhoi = real(cgspresi / (unit_pressure * rhoi)) + spsoundi = real(cgsspsoundi / unit_velocity) + tempi = temperaturei + if (present(mu_local)) mu_local = mui + if (present(gamma_local)) gamma_local = gammai + case default spsoundi = 0. ! avoids compiler warnings ponrhoi = 0. @@ -1433,6 +1447,7 @@ subroutine write_options_eos(iunit) use eos_barotropic, only:write_options_eos_barotropic use eos_piecewise, only:write_options_eos_piecewise use eos_gasradrec, only:write_options_eos_gasradrec + use eos_tillotson, only:write_options_eos_tillotson integer, intent(in) :: iunit write(iunit,"(/,a)") '# options controlling equation of state' diff --git a/src/main/eos_tillotson.f90 b/src/main/eos_tillotson.f90 new file mode 100644 index 000000000..5b5946810 --- /dev/null +++ b/src/main/eos_tillotson.f90 @@ -0,0 +1,88 @@ +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2024 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.github.io/ ! +!--------------------------------------------------------------------------! +module eos_tillotson +! +! EoS from Tillotson et al. (1962) +! +! :References: https://apps.dtic.mil/sti/pdfs/AD0486711.pdf +! +! Implementation from Brundage, A. 2013 +! :DOI: 10.1016/j.proeng.2013.05.053 +! +! :Owner: +! +! :Runtime parameters: None +! +! :Dependencies: infile_utils, io +! + implicit none + +contains +!----------------------------------------------------------------------- +!+ +! EoS from Tillotson et al. (1962) +!+ +!----------------------------------------------------------------------- +subroutine equationofstate_tillotson(rho,energy,temperature,mu,X,Z,pressure,spsound,gamma) + real, intent(in) :: rho,energy,temperature,mu,X,Z + real, intent(out) :: pressure,spsound,gamma + +end subroutine equationofstate_tillotson + +!---------------------------------------------------------------- +!+ +! print eos information +!+ +!---------------------------------------------------------------- +subroutine eos_info_tillotson(iprint) + integer, intent(in) :: iprint + + write(iprint,"(/,a)") ' Tillotson EoS' + +end subroutine eos_info_tillotson + +!----------------------------------------------------------------------- +!+ +! reads equation of state options from the input file +!+ +!----------------------------------------------------------------------- + subroutine read_options_eos_tillotson(name,valstring,imatch,igotall,ierr) + use io, only:fatal + character(len=*), intent(in) :: name,valstring + logical, intent(out) :: imatch,igotall + integer, intent(out) :: ierr + integer, save :: ngot = 0 + character(len=30), parameter :: label = 'eos_tillotson' + + imatch = .true. + select case(trim(name)) +! case('irecomb') +! read(valstring,*,iostat=ierr) irecomb +! if ((irecomb < 0) .or. (irecomb > 3)) call fatal(label,'irecomb = 0,1,2,3') +! ngot = ngot + 1 + case default + imatch = .false. + end select + + igotall = (ngot >= 1) + + end subroutine read_options_eos_tillotson + +!----------------------------------------------------------------------- +!+ +! writes equation of state options to the input file +!+ +!----------------------------------------------------------------------- +subroutine write_options_eos_tillotson(iunit) + use infile_utils, only:write_inopt + integer, intent(in) :: iunit + +! call write_inopt(irecomb,'irecomb','recombination energy to include. 0=H2+H+He, 1=H+He, 2=He, 3=none',iunit) + +end subroutine write_options_eos_tillotson + +end module eos_tillotson \ No newline at end of file From ea43d24908d072952896ece7e5e828b09304a00d Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Mon, 12 Aug 2024 12:27:11 +1000 Subject: [PATCH 02/29] (Tillotson) wrote psuedocode for EOS --- src/setup/setup_solarsystem.f90 | 10 ++++++---- src/utils/utils_ephemeris.f90 | 8 ++++++-- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/src/setup/setup_solarsystem.f90 b/src/setup/setup_solarsystem.f90 index 4dcb7a861..14d8b2568 100644 --- a/src/setup/setup_solarsystem.f90 +++ b/src/setup/setup_solarsystem.f90 @@ -36,10 +36,10 @@ module setup !---------------------------------------------------------------- subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,time,fileprefix) use part, only:nptmass,xyzmh_ptmass,vxyz_ptmass,idust,set_particle_type,& - grainsize,graindens,ndustlarge,ndusttypes + grainsize,graindens,ndustlarge,ndusttypes,ihacc use setbinary, only:set_binary use units, only:set_units,umass,udist,unit_density - use physcon, only:solarm,au,pi,km + use physcon, only:solarm,au,pi,km,solarr use io, only:master,fatal use timestep, only:tmax,dtmax use mpc, only:read_mpc,mpc_entry @@ -106,8 +106,9 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, tmax = norbits*period dtmax = period/dumpsperorbit - filename = find_datafile('Distant.txt',url='https://www.minorplanetcenter.net/iau/MPCORB/') - call read_mpc(filename,nbodies,dat=dat) + nbodies = 0 +! filename = find_datafile('Distant.txt',url='https://www.minorplanetcenter.net/iau/MPCORB/') +! call read_mpc(filename,nbodies,dat=dat) print "(a,i0,a)",' read orbital data for ',nbodies,' minor planets' n = 0 @@ -145,6 +146,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, ! restore the Sun ! nptmass = 1 + xyzmh_ptmass(ihacc,1) = solarr/udist ! ! set mass of all the minor bodies equal ! diff --git a/src/utils/utils_ephemeris.f90 b/src/utils/utils_ephemeris.f90 index 91ad77224..81efb8325 100644 --- a/src/utils/utils_ephemeris.f90 +++ b/src/utils/utils_ephemeris.f90 @@ -88,13 +88,17 @@ subroutine construct_horizons_api_url(object,url,ierr) case('mars') cmd = '499' ! mars barycentre case('earth') - cmd = '399' ! earth-moon barycentre + cmd = '399' ! earth case('venus') cmd = '299' ! venus barycentre case('mercury') cmd = '199' ! mercury barycentre + case('moon') + cmd = '301' ! moon + case('apophis') + cmd = '2004 MN4' ! 99942 Apophis case default - ierr = 1 + cmd = trim(object) end select !if (present(epoch)) then From 91c9f0db6c01bcabb47b8a03780d6de275fb1bc4 Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Mon, 12 Aug 2024 14:18:56 +1000 Subject: [PATCH 03/29] (units) added earth radii to unit options --- src/main/eos.f90 | 1 + src/main/eos_tillotson.f90 | 54 +++++++++++++++++++++++++++++++++++--- 2 files changed, 52 insertions(+), 3 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 85aa8eca9..801afcec4 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -27,6 +27,7 @@ module eos ! 16 = Shen eos ! 17 = polytropic EOS with varying mu (depending on H2 formation) ! 20 = Ideal gas + radiation + various forms of recombination energy from HORMONE (Hirai et al., 2020) +! 23 = Hypervelocity Impact of solids-fluids from Tillotson EOS (Tillotson 1962 - implemented by Brundage A. 2013 ! ! :References: ! Lodato & Pringle (2007) diff --git a/src/main/eos_tillotson.f90 b/src/main/eos_tillotson.f90 index 5b5946810..f7b1a3c91 100644 --- a/src/main/eos_tillotson.f90 +++ b/src/main/eos_tillotson.f90 @@ -17,20 +17,68 @@ module eos_tillotson ! ! :Runtime parameters: None ! -! :Dependencies: infile_utils, io +! :Dependencies: ! implicit none contains !----------------------------------------------------------------------- !+ -! EoS from Tillotson et al. (1962) +! EoS from Tillotson et al. (1962) ; Implementation from Brundage (2013) and Ahrens and O'Keefe (1977) !+ !----------------------------------------------------------------------- subroutine equationofstate_tillotson(rho,energy,temperature,mu,X,Z,pressure,spsound,gamma) real, intent(in) :: rho,energy,temperature,mu,X,Z real, intent(out) :: pressure,spsound,gamma - +! +! Define: rho_0, E_0, V_0 (STP (std T and P) specific volume) +! Calc: rho (g/cm3), A (bulk modulus in kbar) +! +! Internal Energy, E (erg/g) +! E(rho,T) = E_C(rho) + E_T(rho,T) + E_e(rho,T) +! E_T = C_V * T , where E_e ~ 0 (small relative to the other terms) +! +! --- Compressed state --- +! if rho >= rho_0 and E >= 0. +! Define: +! a (Tillotson param), a = 0.5 (polytropic constant - 1, at high Temp) +! b (Tillotson param in kbar), b: defined sth (a+b) = STP Gruneisen param +! b, B, E_0 obtained by fitting to Thomas-Fermi calcs of EOS in 10^2 Mbar range and Hugoniot data. +! gamma = V(del P/ del E)_v , gamma_e = 2/3 (free electron gas), 1/2 (real gas) +! eta (compression, rho/rho_0, and/or V_0/V), mu (strain, eta-1) +! +! pressure(rho,E) = [ a + [ b / ( E / ( E_0 * eta**2 ) + 1. ) ] ] * rho*E + A*mu + B*mu**2 +! +! Cold curve: compression (4th order Runge-Kutta) +! d E_C / d rho = P1(rho,E_C) / rho**2 , (IC: rho = rho_0, E_C = 0) , P1 = pressure @ compressed state (C) +! +! --- Expanded state --- [4x regions] +! Cold curve: expansion same method to solve as compression except RHS piecewise function of eqs 1-4 below +! +! 1: Cold expanded state: +! Define: rho_IV, E_IV (IV = incipient vaporisation) +! if rho_0 > rho >= rho_IV and E <= E_IV +! pressure(rho,E) = (same eq as compressed) +! +! 2: Hot expanded state: +! Define: E_CV (CV = completely vaporised) +! if rho_0 > rho and E >= E_CV +! Define: alpha, beta, a, b, A, B +! pressure(rho,E) = a*rho*E + [ (b*rho*E)/(E/(E_0 * eta**2)+1) + A*mu*exp( -beta*[ (rho_0/rho) -1 ] )] * exp( -alpha*[ (rho_0/rho) -1 ]**2 ) +! As density approaches 0, second term in expression disappears and the material approaches classical Thomas-Fermi limit (high pressures and energy during shock compression) +! P_e = rho*gamma_e *E_e [rho_0 >> rho, E >> E_CV], gamma_e = 2/3 (free electron gas), 1/2 (real gas) +! Pressure in 2. chosen to match Pressure in 1. at rho_0 (init) +! +! 3: Mixing region between gas and liquid (solid): +! Define: P3 = pressure @ 2. , P2 = pressure @ 1. +! if rho_0 > rho > rho_IV and E_IV < E < E_CV +! pressure(rho,E) = [(E-E_IV)*P3 + (E_CV - E)*P2] / (E_CV - E_IV) +! +! 4: Low energy expansion state: +! if rho < rho_IV and E < E_CV +! Define: a, b, A +! pressure(rho,E) = [ a + [ b / ( E / ( E_0 * eta**2 ) + 1 ) ] ] * rho*E + A*mu (same eq as compressed but without B*mu**2 term) +! end subroutine equationofstate_tillotson !---------------------------------------------------------------- From 9db149cf961a7cc4e9fb7286cd25116f405244f1 Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Mon, 12 Aug 2024 14:27:41 +1000 Subject: [PATCH 04/29] (units) added earth radii to unit options --- src/main/units.f90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/main/units.f90 b/src/main/units.f90 index 8af9b5dc9..0f588ff97 100644 --- a/src/main/units.f90 +++ b/src/main/units.f90 @@ -192,6 +192,8 @@ subroutine select_unit(string,unit,ierr) select case(trim(unitstr)) case('solarr','rsun') unit = solarr + case('earthr','rearth') + unit = earthr case('au') unit = au case('ly','lightyear') @@ -281,9 +283,9 @@ logical function is_length_unit(string) call get_unit_multiplier(string,unitstr,fac,ierr) select case(trim(unitstr)) - case('solarr','rsun','au','ly','lightyear','pc','parsec',& - 'kpc','kiloparsec','mpc','megaparsec','km','kilometres',& - 'kilometers','cm','centimetres','centimeters') + case('solarr','rsun','rearth','earthr','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. From c8c7280123efea7a3458609336e997810157cfee Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Tue, 13 Aug 2024 10:43:07 +1000 Subject: [PATCH 05/29] (tillotson) added psuedocode - tillotson eos parameters for basalt and ice --- src/main/eos_tillotson.f90 | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/main/eos_tillotson.f90 b/src/main/eos_tillotson.f90 index f7b1a3c91..4f0532684 100644 --- a/src/main/eos_tillotson.f90 +++ b/src/main/eos_tillotson.f90 @@ -79,6 +79,18 @@ subroutine equationofstate_tillotson(rho,energy,temperature,mu,X,Z,pressure,spso ! Define: a, b, A ! pressure(rho,E) = [ a + [ b / ( E / ( E_0 * eta**2 ) + 1 ) ] ] * rho*E + A*mu (same eq as compressed but without B*mu**2 term) ! +! ---Tillotson EOS params [new subroutine?]--- +! From Benz and Asphaug 1999 [Table 2] DOI: 10.1006/icar.1999.6204 +! ---Basalt--- +! rho_0 = 2.7 [g/cm3]; A = 2.67e11 [erg/cm3]; B = 2.67e11[erg/cm3]; +! E_0 = 4.87e12 [erg/g]; E_IV = 4.72e10 [erg/g]; E_CV = 1.82e11 [erg/g]; +! a = 0.5; b = 1.5; alpha = 5.0; beta = 5.0 +! +! ---Ice--- +! rho_0 = 0.917 [g/cm3]; A = 9.47e10 [erg/cm3]; B = 9.47e10[erg/cm3]; +! E_0 = 1.00e11 [erg/g]; E_IV = 7.73e9 [erg/g]; E_CV = 3.04e10 [erg/g]; +! a = 0.3; b = 0.1; alpha = 10.0; beta = 5.0 +! end subroutine equationofstate_tillotson !---------------------------------------------------------------- From 93466fa1af0ab79a59b98e9815b32656731bd788 Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Wed, 14 Aug 2024 16:35:33 +1000 Subject: [PATCH 06/29] (tillotson) compiled code and params for basalt from benz & asphaug 1999 --- src/main/eos.f90 | 6 ++- src/main/eos_tillotson.f90 | 99 +++++++++++--------------------------- 2 files changed, 33 insertions(+), 72 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 801afcec4..56908fd04 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -442,10 +442,12 @@ subroutine equationofstate(eos_type,ponrhoi,spsoundi,rhoi,xi,yi,zi,tempi,eni,gam ! ! from Tillotson (1962) ! - call equationofstate_tillotson(cgsrhoi,cgseni*cgsrhoi,temperaturei,mui,X_i,1.-X_i-Z_i,cgspresi,cgsspsoundi,gammai) + cgsrhoi = rhoi * unit_density + cgseni = eni * unit_ergg + call equationofstate_tillotson(cgsrhoi,cgseni,cgspresi,cgsspsoundi,gammai) ponrhoi = real(cgspresi / (unit_pressure * rhoi)) spsoundi = real(cgsspsoundi / unit_velocity) - tempi = temperaturei + tempi = 0. !temperaturei if (present(mu_local)) mu_local = mui if (present(gamma_local)) gamma_local = gammai diff --git a/src/main/eos_tillotson.f90 b/src/main/eos_tillotson.f90 index 4f0532684..1d49bfb9f 100644 --- a/src/main/eos_tillotson.f90 +++ b/src/main/eos_tillotson.f90 @@ -10,8 +10,7 @@ module eos_tillotson ! ! :References: https://apps.dtic.mil/sti/pdfs/AD0486711.pdf ! -! Implementation from Brundage, A. 2013 -! :DOI: 10.1016/j.proeng.2013.05.053 +! Implementation from Benz and Asphaug (1999) ! ! :Owner: ! @@ -20,77 +19,37 @@ module eos_tillotson ! :Dependencies: ! implicit none + real :: rho0 = 2.7 ! g/cm^3 zero-pressure density (Basalt) from Benz & Asphaug 1999 +! aparam, bparam, A, B, energy0 material-dependent Tillotson parameters (Basalt) + real :: aparam = 0.5 , bparam = 1.5 , alpha = 5. , beta = 5. + real :: A = 2.67e11 , B = 2.67e11 ! erg/cm^3 + real :: energy0 = 4.87e12 ! erg/g + real :: eta, mu + private :: rho0, aparam, bparam, A, B, alpha, beta, energy0, eta, mu contains !----------------------------------------------------------------------- !+ -! EoS from Tillotson et al. (1962) ; Implementation from Brundage (2013) and Ahrens and O'Keefe (1977) +! EoS from Tillotson et al. (1962) ; Implementation from Benz and Asphaug (1999) !+ !----------------------------------------------------------------------- -subroutine equationofstate_tillotson(rho,energy,temperature,mu,X,Z,pressure,spsound,gamma) - real, intent(in) :: rho,energy,temperature,mu,X,Z +subroutine equationofstate_tillotson(rho,energy,pressure,spsound,gamma) + real, intent(in) :: rho,energy real, intent(out) :: pressure,spsound,gamma -! -! Define: rho_0, E_0, V_0 (STP (std T and P) specific volume) -! Calc: rho (g/cm3), A (bulk modulus in kbar) ! -! Internal Energy, E (erg/g) -! E(rho,T) = E_C(rho) + E_T(rho,T) + E_e(rho,T) -! E_T = C_V * T , where E_e ~ 0 (small relative to the other terms) -! -! --- Compressed state --- -! if rho >= rho_0 and E >= 0. -! Define: -! a (Tillotson param), a = 0.5 (polytropic constant - 1, at high Temp) -! b (Tillotson param in kbar), b: defined sth (a+b) = STP Gruneisen param -! b, B, E_0 obtained by fitting to Thomas-Fermi calcs of EOS in 10^2 Mbar range and Hugoniot data. -! gamma = V(del P/ del E)_v , gamma_e = 2/3 (free electron gas), 1/2 (real gas) -! eta (compression, rho/rho_0, and/or V_0/V), mu (strain, eta-1) -! -! pressure(rho,E) = [ a + [ b / ( E / ( E_0 * eta**2 ) + 1. ) ] ] * rho*E + A*mu + B*mu**2 -! -! Cold curve: compression (4th order Runge-Kutta) -! d E_C / d rho = P1(rho,E_C) / rho**2 , (IC: rho = rho_0, E_C = 0) , P1 = pressure @ compressed state (C) -! -! --- Expanded state --- [4x regions] -! Cold curve: expansion same method to solve as compression except RHS piecewise function of eqs 1-4 below -! -! 1: Cold expanded state: -! Define: rho_IV, E_IV (IV = incipient vaporisation) -! if rho_0 > rho >= rho_IV and E <= E_IV -! pressure(rho,E) = (same eq as compressed) -! -! 2: Hot expanded state: -! Define: E_CV (CV = completely vaporised) -! if rho_0 > rho and E >= E_CV -! Define: alpha, beta, a, b, A, B -! pressure(rho,E) = a*rho*E + [ (b*rho*E)/(E/(E_0 * eta**2)+1) + A*mu*exp( -beta*[ (rho_0/rho) -1 ] )] * exp( -alpha*[ (rho_0/rho) -1 ]**2 ) -! As density approaches 0, second term in expression disappears and the material approaches classical Thomas-Fermi limit (high pressures and energy during shock compression) -! P_e = rho*gamma_e *E_e [rho_0 >> rho, E >> E_CV], gamma_e = 2/3 (free electron gas), 1/2 (real gas) -! Pressure in 2. chosen to match Pressure in 1. at rho_0 (init) -! -! 3: Mixing region between gas and liquid (solid): -! Define: P3 = pressure @ 2. , P2 = pressure @ 1. -! if rho_0 > rho > rho_IV and E_IV < E < E_CV -! pressure(rho,E) = [(E-E_IV)*P3 + (E_CV - E)*P2] / (E_CV - E_IV) -! -! 4: Low energy expansion state: -! if rho < rho_IV and E < E_CV -! Define: a, b, A -! pressure(rho,E) = [ a + [ b / ( E / ( E_0 * eta**2 ) + 1 ) ] ] * rho*E + A*mu (same eq as compressed but without B*mu**2 term) -! -! ---Tillotson EOS params [new subroutine?]--- -! From Benz and Asphaug 1999 [Table 2] DOI: 10.1006/icar.1999.6204 -! ---Basalt--- -! rho_0 = 2.7 [g/cm3]; A = 2.67e11 [erg/cm3]; B = 2.67e11[erg/cm3]; -! E_0 = 4.87e12 [erg/g]; E_IV = 4.72e10 [erg/g]; E_CV = 1.82e11 [erg/g]; -! a = 0.5; b = 1.5; alpha = 5.0; beta = 5.0 -! -! ---Ice--- -! rho_0 = 0.917 [g/cm3]; A = 9.47e10 [erg/cm3]; B = 9.47e10[erg/cm3]; -! E_0 = 1.00e11 [erg/g]; E_IV = 7.73e9 [erg/g]; E_CV = 3.04e10 [erg/g]; -! a = 0.3; b = 0.1; alpha = 10.0; beta = 5.0 -! +eta = rho / rho0 +mu = eta - 1. + +if (rho >= rho0 .and. energy >= 0.) then + pressure = ( aparam + ( bparam / ( energy / ( energy0 * eta**2 ) + 1. ) ) ) * rho*energy + A*mu + B*mu**2 +else + pressure = aparam*rho*energy + ( (bparam*rho*energy)/(energy/(energy0 * eta**2) + 1.) + & + A*mu*exp(-beta*( (rho0/rho) - 1.))) * exp(-alpha*( (rho0/rho) - 1.)**2) +endif + +gamma = 0.5 ! gamma = V(del P/ del E)_v , gamma_e = 2/3 (free electron gas), 1/2 (real gas) +spsound = sqrt(A / rho) ! wave speed estimate + end subroutine equationofstate_tillotson !---------------------------------------------------------------- @@ -120,10 +79,10 @@ subroutine read_options_eos_tillotson(name,valstring,imatch,igotall,ierr) imatch = .true. select case(trim(name)) -! case('irecomb') -! read(valstring,*,iostat=ierr) irecomb -! if ((irecomb < 0) .or. (irecomb > 3)) call fatal(label,'irecomb = 0,1,2,3') -! ngot = ngot + 1 + case('rho0') + read(valstring,*,iostat=ierr) rho0 + if ((rho0 < 0.)) call fatal(label,'rho0 < 0') + ngot = ngot + 1 case default imatch = .false. end select @@ -141,7 +100,7 @@ subroutine write_options_eos_tillotson(iunit) use infile_utils, only:write_inopt integer, intent(in) :: iunit -! call write_inopt(irecomb,'irecomb','recombination energy to include. 0=H2+H+He, 1=H+He, 2=He, 3=none',iunit) + call write_inopt(rho0,'rho0','reference density',iunit) end subroutine write_options_eos_tillotson From 665fcb3039e54acbf32a85879fce93586b5a0dca Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Thu, 15 Aug 2024 10:30:29 +1000 Subject: [PATCH 07/29] (tillotson) added initialise eos case 23 subroutine --- src/main/eos.f90 | 4 +++- src/main/eos_tillotson.f90 | 4 ++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 56908fd04..c89ec799b 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -477,6 +477,7 @@ subroutine init_eos(eos_type,ierr) use eos_gasradrec, only:init_eos_gasradrec use eos_HIIR, only:init_eos_HIIR use dim, only:maxvxyzu,do_radiation + use eos_tillotson, only:init_eos_tillotson integer, intent(in) :: eos_type integer, intent(out) :: ierr integer :: ierr_mesakapp @@ -556,7 +557,8 @@ subroutine init_eos(eos_type,ierr) case(21,22) call init_eos_HIIR() - + case(23) + call init_eos_tillotson(ierr) end select done_init_eos = .true. diff --git a/src/main/eos_tillotson.f90 b/src/main/eos_tillotson.f90 index 1d49bfb9f..0d3e59b04 100644 --- a/src/main/eos_tillotson.f90 +++ b/src/main/eos_tillotson.f90 @@ -33,6 +33,10 @@ module eos_tillotson ! EoS from Tillotson et al. (1962) ; Implementation from Benz and Asphaug (1999) !+ !----------------------------------------------------------------------- +subroutine init_eos_tillotson(ierr) + integer, intent(out) :: ierr +end subroutine init_eos_tillotson + subroutine equationofstate_tillotson(rho,energy,pressure,spsound,gamma) real, intent(in) :: rho,energy real, intent(out) :: pressure,spsound,gamma From 00c42fe73435428aacad3e670028823c1702acc1 Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Thu, 15 Aug 2024 11:39:53 +1000 Subject: [PATCH 08/29] (tillotson) update init and read eos subroutines for case 23 --- src/main/eos.f90 | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index c89ec799b..cc492f8b4 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -438,18 +438,16 @@ subroutine equationofstate(eos_type,ponrhoi,spsoundi,rhoi,xi,yi,zi,tempi,eni,gam case(23) ! - !-- Tillotson EOS + !-- Tillotson EOS from Tillotson (1962); Implementation from Benz and Asphaug (1999) ! - ! from Tillotson (1962) - ! - cgsrhoi = rhoi * unit_density - cgseni = eni * unit_ergg - call equationofstate_tillotson(cgsrhoi,cgseni,cgspresi,cgsspsoundi,gammai) - ponrhoi = real(cgspresi / (unit_pressure * rhoi)) - spsoundi = real(cgsspsoundi / unit_velocity) - tempi = 0. !temperaturei - if (present(mu_local)) mu_local = mui - if (present(gamma_local)) gamma_local = gammai + ! cgsrhoi = rhoi * unit_density + ! cgseni = eni * unit_ergg + call equationofstate_tillotson(rhoi,eni,presi,spsoundi,gammai) !cgsrhoi,cgseni,cgspresi,cgsspsoundi,gammai + ! ponrhoi = real(cgspresi / (unit_pressure * rhoi)) + ! spsoundi = real(cgsspsoundi / unit_velocity) + ! tempi = 0. !temperaturei + ! if (present(mu_local)) mu_local = mui + ! if (present(gamma_local)) gamma_local = gammai case default spsoundi = 0. ! avoids compiler warnings @@ -557,8 +555,11 @@ subroutine init_eos(eos_type,ierr) case(21,22) call init_eos_HIIR() + case(23) + call init_eos_tillotson(ierr) + end select done_init_eos = .true. @@ -1509,6 +1510,7 @@ subroutine read_options_eos(name,valstring,imatch,igotall,ierr) use eos_barotropic, only:read_options_eos_barotropic use eos_piecewise, only:read_options_eos_piecewise use eos_gasradrec, only:read_options_eos_gasradrec + use eos_tillotson, only:read_options_eos_tillotson character(len=*), intent(in) :: name,valstring logical, intent(out) :: imatch,igotall integer, intent(out) :: ierr @@ -1548,6 +1550,7 @@ subroutine read_options_eos(name,valstring,imatch,igotall,ierr) if (.not.imatch .and. ieos== 8) call read_options_eos_barotropic(name,valstring,imatch,igotall_barotropic,ierr) if (.not.imatch .and. ieos== 9) call read_options_eos_piecewise( name,valstring,imatch,igotall_piecewise, ierr) if (.not.imatch .and. ieos==20) call read_options_eos_gasradrec( name,valstring,imatch,igotall_gasradrec, ierr) + if (.not.imatch .and. ieos==23) call read_options_eos_tillotson(name,valstring,imatch,igotall,ierr) !--make sure we have got all compulsory options (otherwise, rewrite input file) igotall = (ngot >= 1) .and. igotall_piecewise .and. igotall_barotropic .and. igotall_gasradrec From bf14df6bb388a3d8f4cd236844592dbf949d56f0 Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Thu, 15 Aug 2024 16:41:26 +1000 Subject: [PATCH 09/29] (tillotson) updated eos and readwrite infile to include case 23 eos --- src/main/eos.f90 | 2 +- src/main/eos_tillotson.f90 | 26 +++++++++++++++++++++++++- src/main/readwrite_infile.F90 | 2 +- 3 files changed, 27 insertions(+), 3 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index cc492f8b4..56bfee1d2 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -50,7 +50,7 @@ module eos use part, only:ien_etotal,ien_entropy,ien_type use dim, only:gr implicit none - integer, parameter, public :: maxeos = 22 + integer, parameter, public :: maxeos = 23 real, public :: polyk, polyk2, gamma real, public :: qfacdisc = 0.75, qfacdisc2 = 0.75 logical, public :: extract_eos_from_hdr = .false. diff --git a/src/main/eos_tillotson.f90 b/src/main/eos_tillotson.f90 index 0d3e59b04..0b6e0216d 100644 --- a/src/main/eos_tillotson.f90 +++ b/src/main/eos_tillotson.f90 @@ -18,6 +18,7 @@ module eos_tillotson ! ! :Dependencies: ! + use units, only: unit_density, unit_energ, umass, udist implicit none real :: rho0 = 2.7 ! g/cm^3 zero-pressure density (Basalt) from Benz & Asphaug 1999 ! aparam, bparam, A, B, energy0 material-dependent Tillotson parameters (Basalt) @@ -35,6 +36,20 @@ module eos_tillotson !----------------------------------------------------------------------- subroutine init_eos_tillotson(ierr) integer, intent(out) :: ierr + +! call write_options_eos_tillotson(iunit) +! +! Convert to code units +! +rho0 = real(rho0/unit_density) +A = real((A/unit_energ)*(udist**3)) ! erg/cm^3 +B = real((B/unit_energ)*(udist**3)) ! erg/cm^3 +energy0 = real((energy0/unit_energ) *umass) ! erg/g + +! call read_options_eos_tillotson(name,valstring,imatch,igotall,ierr) +! call equationofstate_tillotson(rho,energy,pressure,spsound,gamma) + + end subroutine init_eos_tillotson subroutine equationofstate_tillotson(rho,energy,pressure,spsound,gamma) @@ -104,7 +119,16 @@ subroutine write_options_eos_tillotson(iunit) use infile_utils, only:write_inopt integer, intent(in) :: iunit - call write_inopt(rho0,'rho0','reference density',iunit) + call write_inopt(rho0,'rho0','reference density g/cm^3',iunit) + call write_inopt(aparam,'aparam','material-dependent Tillotson parameter, unitless',iunit) + call write_inopt(bparam,'bparam','material-dependent Tillotson parameter, unitless',iunit) + call write_inopt(A,'A','material-dependent Tillotson parameter, erg/cm^3',iunit) + call write_inopt(B,'B','material-dependent Tillotson parameter, erg/cm^3',iunit) + call write_inopt(alpha,'alpha','material-dependent Tillotson parameter, unitless',iunit) + call write_inopt(beta,'beta','material-dependent Tillotson parameter, unitless',iunit) + call write_inopt(energy0,'energy0','material-dependent Tillotson parameter, erg/g',iunit) + call write_inopt(eta,'eta','compression, (rho/rho0), unitless',iunit) + call write_inopt(mu,'mu','strain, (eta-1), unitless',iunit) end subroutine write_options_eos_tillotson diff --git a/src/main/readwrite_infile.F90 b/src/main/readwrite_infile.F90 index 450b42f88..29e0235ee 100644 --- a/src/main/readwrite_infile.F90 +++ b/src/main/readwrite_infile.F90 @@ -693,7 +693,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/=23)) & 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 03f849317d79f58636ed77db247176d2b2d8e0e9 Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Sat, 17 Aug 2024 15:48:15 +1000 Subject: [PATCH 10/29] (tillotson) remove unit conversion from write options eos --- src/main/eos_tillotson.f90 | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/src/main/eos_tillotson.f90 b/src/main/eos_tillotson.f90 index 0b6e0216d..b7131dc71 100644 --- a/src/main/eos_tillotson.f90 +++ b/src/main/eos_tillotson.f90 @@ -119,17 +119,6 @@ subroutine write_options_eos_tillotson(iunit) use infile_utils, only:write_inopt integer, intent(in) :: iunit - call write_inopt(rho0,'rho0','reference density g/cm^3',iunit) - call write_inopt(aparam,'aparam','material-dependent Tillotson parameter, unitless',iunit) - call write_inopt(bparam,'bparam','material-dependent Tillotson parameter, unitless',iunit) - call write_inopt(A,'A','material-dependent Tillotson parameter, erg/cm^3',iunit) - call write_inopt(B,'B','material-dependent Tillotson parameter, erg/cm^3',iunit) - call write_inopt(alpha,'alpha','material-dependent Tillotson parameter, unitless',iunit) - call write_inopt(beta,'beta','material-dependent Tillotson parameter, unitless',iunit) - call write_inopt(energy0,'energy0','material-dependent Tillotson parameter, erg/g',iunit) - call write_inopt(eta,'eta','compression, (rho/rho0), unitless',iunit) - call write_inopt(mu,'mu','strain, (eta-1), unitless',iunit) - end subroutine write_options_eos_tillotson end module eos_tillotson \ No newline at end of file From cd8df39eff02fad9323e9226a371bf41a393203f Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Tue, 20 Aug 2024 13:12:45 +1000 Subject: [PATCH 11/29] (setup_solarsystem) added mass of sun in binary system setup --- src/setup/setup_solarsystem.f90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/setup/setup_solarsystem.f90 b/src/setup/setup_solarsystem.f90 index 14d8b2568..19cf321bc 100644 --- a/src/setup/setup_solarsystem.f90 +++ b/src/setup/setup_solarsystem.f90 @@ -146,7 +146,9 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, ! restore the Sun ! nptmass = 1 + xyzmh_ptmass(1:3,1) = 0. xyzmh_ptmass(ihacc,1) = solarr/udist + xyzmh_ptmass(4,1) = solarm/umass ! ! set mass of all the minor bodies equal ! From 0a5f9013092d40caeb991781719e209641530900 Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Tue, 20 Aug 2024 13:14:00 +1000 Subject: [PATCH 12/29] (setup_binary) removed hardcoded eos in subroutine --- src/setup/setup_binary.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/setup/setup_binary.f90 b/src/setup/setup_binary.f90 index 0b6c87bd5..7bcba69d9 100644 --- a/src/setup/setup_binary.f90 +++ b/src/setup/setup_binary.f90 @@ -93,7 +93,6 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,& call set_defaults_orbit(orbit) relax = .true. corotate = .false. - ieos = 2 if (id==master) print "(/,65('-'),1(/,a),/,65('-'),/)",& ' Welcome to the Ultimate Binary Setup' From 19a988be2358dec43b2b9a7d271b6bf84a1bb6fd Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Wed, 21 Aug 2024 10:47:55 +1000 Subject: [PATCH 13/29] (tillotson) fixed units read-write functions and updated condensed-evaporated cases --- src/main/eos.f90 | 21 ++++--- src/main/eos_tillotson.f90 | 115 +++++++++++++++++++++++++++---------- 2 files changed, 98 insertions(+), 38 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 56bfee1d2..924eba8dd 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -440,14 +440,12 @@ subroutine equationofstate(eos_type,ponrhoi,spsoundi,rhoi,xi,yi,zi,tempi,eni,gam ! !-- Tillotson EOS from Tillotson (1962); Implementation from Benz and Asphaug (1999) ! - ! cgsrhoi = rhoi * unit_density - ! cgseni = eni * unit_ergg - call equationofstate_tillotson(rhoi,eni,presi,spsoundi,gammai) !cgsrhoi,cgseni,cgspresi,cgsspsoundi,gammai - ! ponrhoi = real(cgspresi / (unit_pressure * rhoi)) - ! spsoundi = real(cgsspsoundi / unit_velocity) + cgsrhoi = rhoi * unit_density + cgseni = eni * unit_ergg + call equationofstate_tillotson(cgsrhoi,cgseni,cgspresi,cgsspsoundi,gammai) + ponrhoi = real(cgspresi / (unit_pressure * rhoi)) + spsoundi = real(cgsspsoundi / unit_velocity) ! tempi = 0. !temperaturei - ! if (present(mu_local)) mu_local = mui - ! if (present(gamma_local)) gamma_local = gammai case default spsoundi = 0. ! avoids compiler warnings @@ -1495,6 +1493,8 @@ subroutine write_options_eos(iunit) call write_inopt(X_in,'X','H mass fraction (ignored if variable composition)',iunit) call write_inopt(Z_in,'Z','metallicity (ignored if variable composition)',iunit) endif + case(23) + call write_options_eos_tillotson(iunit) end select end subroutine write_options_eos @@ -1517,11 +1517,13 @@ subroutine read_options_eos(name,valstring,imatch,igotall,ierr) integer, save :: ngot = 0 character(len=30), parameter :: label = 'read_options_eos' logical :: igotall_barotropic,igotall_piecewise,igotall_gasradrec + logical :: igotall_tillotson imatch = .true. igotall_barotropic = .true. igotall_piecewise = .true. igotall_gasradrec = .true. + igotall_tillotson = .true. select case(trim(name)) case('ieos') @@ -1550,10 +1552,11 @@ subroutine read_options_eos(name,valstring,imatch,igotall,ierr) if (.not.imatch .and. ieos== 8) call read_options_eos_barotropic(name,valstring,imatch,igotall_barotropic,ierr) if (.not.imatch .and. ieos== 9) call read_options_eos_piecewise( name,valstring,imatch,igotall_piecewise, ierr) if (.not.imatch .and. ieos==20) call read_options_eos_gasradrec( name,valstring,imatch,igotall_gasradrec, ierr) - if (.not.imatch .and. ieos==23) call read_options_eos_tillotson(name,valstring,imatch,igotall,ierr) + if (.not.imatch .and. ieos==23) call read_options_eos_tillotson(name,valstring,imatch,igotall_tillotson,ierr) !--make sure we have got all compulsory options (otherwise, rewrite input file) - igotall = (ngot >= 1) .and. igotall_piecewise .and. igotall_barotropic .and. igotall_gasradrec + igotall = (ngot >= 1) .and. igotall_piecewise .and. igotall_barotropic .and. igotall_gasradrec & + .and. igotall_tillotson end subroutine read_options_eos diff --git a/src/main/eos_tillotson.f90 b/src/main/eos_tillotson.f90 index b7131dc71..dcddc0291 100644 --- a/src/main/eos_tillotson.f90 +++ b/src/main/eos_tillotson.f90 @@ -18,15 +18,13 @@ module eos_tillotson ! ! :Dependencies: ! - use units, only: unit_density, unit_energ, umass, udist implicit none real :: rho0 = 2.7 ! g/cm^3 zero-pressure density (Basalt) from Benz & Asphaug 1999 ! aparam, bparam, A, B, energy0 material-dependent Tillotson parameters (Basalt) real :: aparam = 0.5 , bparam = 1.5 , alpha = 5. , beta = 5. real :: A = 2.67e11 , B = 2.67e11 ! erg/cm^3 - real :: energy0 = 4.87e12 ! erg/g - real :: eta, mu - private :: rho0, aparam, bparam, A, B, alpha, beta, energy0, eta, mu + real :: u0 = 4.87e12 ! erg/g + private :: rho0, aparam, bparam, A, B, alpha, beta, u0 contains !----------------------------------------------------------------------- @@ -37,37 +35,59 @@ module eos_tillotson subroutine init_eos_tillotson(ierr) integer, intent(out) :: ierr -! call write_options_eos_tillotson(iunit) -! -! Convert to code units -! -rho0 = real(rho0/unit_density) -A = real((A/unit_energ)*(udist**3)) ! erg/cm^3 -B = real((B/unit_energ)*(udist**3)) ! erg/cm^3 -energy0 = real((energy0/unit_energ) *umass) ! erg/g - -! call read_options_eos_tillotson(name,valstring,imatch,igotall,ierr) -! call equationofstate_tillotson(rho,energy,pressure,spsound,gamma) - + ierr = 0 end subroutine init_eos_tillotson -subroutine equationofstate_tillotson(rho,energy,pressure,spsound,gamma) - real, intent(in) :: rho,energy +subroutine equationofstate_tillotson(rho,u,pressure,spsound,gamma) + real, intent(in) :: rho,u real, intent(out) :: pressure,spsound,gamma -! + real :: eta, mu, neu, omega, spsoundmin + eta = rho / rho0 mu = eta - 1. - -if (rho >= rho0 .and. energy >= 0.) then - pressure = ( aparam + ( bparam / ( energy / ( energy0 * eta**2 ) + 1. ) ) ) * rho*energy + A*mu + B*mu**2 -else - pressure = aparam*rho*energy + ( (bparam*rho*energy)/(energy/(energy0 * eta**2) + 1.) + & - A*mu*exp(-beta*( (rho0/rho) - 1.))) * exp(-alpha*( (rho0/rho) - 1.)**2) +neu = (1. / eta) - 1. ! Kegerreis et al. 2020 +omega = (u / (u0*eta**2) + 1.) + +if (rho >= rho0 .and. u >= 0.) then ! Condensed state + ! pressure = ( aparam + ( bparam / ( u / ( u0 * eta**2 ) + 1. ) ) ) * rho*u + A*mu + B*mu**2 + pressure = ( aparam + ( bparam / omega ) ) * rho*u + A*mu + B*mu**2 + if (pressure <= 0.) then + pressure = 0. + endif + ! sound speed from ! Kegereis et al. 2020 + spsound = sqrt(pressure / rho * (1. + aparam + (bparam / omega )) + & + (bparam*(omega - 1.) / omega**2) * (2.*u - pressure/rho)) + & + ((1./rho)*(A + B*((eta**2) -1.))) + spsoundmin = sqrt(A / rho) ! wave speed estimate + if (spsound < spsoundmin) then + spsound = spsoundmin + endif +else ! Expanded state + ! pressure = aparam*rho*u + ( (bparam*rho*u)/(u/(u0 * eta**2) + 1.) + & + ! A*mu*exp(-beta*( (rho0/rho) - 1.))) * exp(-alpha*( (rho0/rho) - 1.)**2) + pressure = aparam*rho*u + ( ((bparam*rho*u)/omega) + & + A*mu*exp(-beta*neu)) * exp(-alpha*(neu**2)) + if (pressure <= 0.) then + pressure = 0. + endif + spsound = sqrt( (pressure/rho)*( 1. + aparam + ((bparam / omega )*exp(-alpha*neu**2))) + & + ( ((bparam*rho*u)/((omega**2) * (eta**2))) * ((1/(u0*rho))*(2.*u - pressure/rho) & + + ((2.*alpha*eta*omega)/rho0)) + ((A/rho0)*(1.+ (mu/eta**2)*(beta*2.*alpha - eta))) & + * exp(-beta*neu) )*exp(-alpha*neu**2) ) + spsoundmin = sqrt(A / rho) ! wave speed estimate + if (spsound < spsoundmin) then + spsound = spsoundmin + endif endif - -gamma = 0.5 ! gamma = V(del P/ del E)_v , gamma_e = 2/3 (free electron gas), 1/2 (real gas) -spsound = sqrt(A / rho) ! wave speed estimate +! Hybrid Pressure (inbetween) +! uiv = +! ucv = +! frac = ((u / uiv) - 1.) / ((ucv/uiv) - 1.) +! pressure = +! ce = spsound @ expanded +! cc = spsound @ condensed +! spsound = sqrt(( ((u - uiv)* ce**2) + ((ucv - u)*cc**2))/(ucv - uiv)) end subroutine equationofstate_tillotson @@ -102,11 +122,39 @@ subroutine read_options_eos_tillotson(name,valstring,imatch,igotall,ierr) read(valstring,*,iostat=ierr) rho0 if ((rho0 < 0.)) call fatal(label,'rho0 < 0') ngot = ngot + 1 + case('aparam') + read(valstring,*,iostat=ierr) aparam + if ((aparam < 0.)) call fatal(label,'aparam < 0') + ngot = ngot + 1 + case('bparam') + read(valstring,*,iostat=ierr) bparam + if ((bparam < 0.)) call fatal(label,'bparam < 0') + ngot = ngot + 1 + case('A') + read(valstring,*,iostat=ierr) A + if ((A < 0.)) call fatal(label,'A < 0') + ngot = ngot + 1 + case('B') + read(valstring,*,iostat=ierr) B + if ((B < 0.)) call fatal(label,'B < 0') + ngot = ngot + 1 + case('alpha_t') + read(valstring,*,iostat=ierr) alpha + if ((alpha < 0.)) call fatal(label,'alpha < 0') + ngot = ngot + 1 + case('beta_t') + read(valstring,*,iostat=ierr) beta + if ((beta < 0.)) call fatal(label,'beta < 0') + ngot = ngot + 1 + case('u0') + read(valstring,*,iostat=ierr) u0 + if ((u0 < 0.)) call fatal(label,'u0 < 0') + ngot = ngot + 1 case default imatch = .false. end select - igotall = (ngot >= 1) + igotall = (ngot >= 8) end subroutine read_options_eos_tillotson @@ -118,6 +166,15 @@ end subroutine read_options_eos_tillotson subroutine write_options_eos_tillotson(iunit) use infile_utils, only:write_inopt integer, intent(in) :: iunit + + call write_inopt(rho0,'rho0','reference density g/cm^3',iunit) + call write_inopt(aparam,'aparam','material-dependent Tillotson parameter, unitless',iunit) + call write_inopt(bparam,'bparam','material-dependent Tillotson parameter, unitless',iunit) + call write_inopt(A,'A','material-dependent Tillotson parameter, erg/cm^3',iunit) + call write_inopt(B,'B','material-dependent Tillotson parameter, erg/cm^3',iunit) + call write_inopt(alpha,'alpha_t','material-dependent Tillotson parameter, unitless',iunit) + call write_inopt(beta,'beta_t','material-dependent Tillotson parameter, unitless',iunit) + call write_inopt(u0,'u0','material-dependent Tillotson parameter, erg/g',iunit) end subroutine write_options_eos_tillotson From 14b9bc2f2a6e43c0fa92f11c2769165a3dd4f5f9 Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Tue, 3 Sep 2024 14:03:18 +1000 Subject: [PATCH 14/29] (tillotson) added calculation of temp to internal energy to eos subroutine call --- src/main/eos.f90 | 9 ++++++++- src/main/eos_tillotson.f90 | 31 +++++++++++++++++++++++++++++++ 2 files changed, 39 insertions(+), 1 deletion(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 924eba8dd..a89a5535d 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -860,13 +860,15 @@ subroutine calc_temp_and_ene(eos_type,rho,pres,ene,temp,ierr,guesseint,mu_local, use eos_idealplusrad, only:get_idealgasplusrad_tempfrompres,get_idealplusrad_enfromtemp use eos_mesa, only:get_eos_eT_from_rhop_mesa use eos_gasradrec, only:calc_uT_from_rhoP_gasradrec + use eos_tillotson, only:calc_uT_from_rhoP_tillotson + use units, only:unit_ergg,unit_pressure,unit_density integer, intent(in) :: eos_type real, intent(in) :: rho,pres real, intent(inout) :: ene,temp real, intent(in), optional :: guesseint,X_local,Z_local real, intent(inout), optional :: mu_local integer, intent(out) :: ierr - real :: mu,X,Z + real :: mu,X,Z,cgseni,cgsrhoi,cgspresi ierr = 0 mu = gmw @@ -887,6 +889,11 @@ subroutine calc_temp_and_ene(eos_type,rho,pres,ene,temp,ierr,guesseint,mu_local, case(20) ! Ideal gas + radiation + recombination (from HORMONE, Hirai et al., 2020) call calc_uT_from_rhoP_gasradrec(rho,pres,X,1.-X-Z,temp,ene,mu,ierr) if (present(mu_local)) mu_local = mu + case(23) ! tillotson + cgsrhoi = rho * unit_density + cgspresi = pres * unit_pressure + call calc_uT_from_rhoP_tillotson(cgsrhoi,cgspresi,temp,cgseni,ierr) + ene = cgseni / unit_ergg case default ierr = 1 end select diff --git a/src/main/eos_tillotson.f90 b/src/main/eos_tillotson.f90 index dcddc0291..46ac7a42c 100644 --- a/src/main/eos_tillotson.f90 +++ b/src/main/eos_tillotson.f90 @@ -103,6 +103,37 @@ subroutine eos_info_tillotson(iprint) end subroutine eos_info_tillotson +!----------------------------------------------------------------------- +!+ +! calc internal energy from pressure and density +!+ +!----------------------------------------------------------------------- +subroutine calc_uT_from_rhoP_tillotson(rho,pres,temp,u,ierr) + real, intent(in) :: rho,pres + real, intent(out) :: temp,u + integer, intent(out) :: ierr + real :: eta, mu, neu, omega, expterm + + eta = rho / rho0 + mu = eta - 1. + neu = (1. / eta) - 1. ! Kegerreis et al. 2020 + omega = (u / (u0*eta**2) + 1.) + + if (rho >= rho0) then ! Condensed state + u = (pres - A*mu - B*mu**2) / ( aparam + ( bparam / omega ) * rho) + else + expterm = exp(-alpha*(neu**2)) + u = (pres - A*mu*exp(-beta*neu)*expterm ) & + / (aparam*rho + ((bparam*rho)/omega)*expterm) + endif + temp = 300. + ierr = 0 + if (u<0.) then + u = 0. + ierr = 1 + endif + +end subroutine calc_uT_from_rhoP_tillotson !----------------------------------------------------------------------- !+ ! reads equation of state options from the input file From 5063782ad9314d0478919c4a5207dfa1556b678d Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Tue, 3 Sep 2024 14:04:36 +1000 Subject: [PATCH 15/29] (solarsystem) included asteroid apophis into ephemeris and solarsystem setup --- src/setup/setup_solarsystem.f90 | 2 +- src/utils/utils_ephemeris.f90 | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/setup/setup_solarsystem.f90 b/src/setup/setup_solarsystem.f90 index 19cf321bc..d42611364 100644 --- a/src/setup/setup_solarsystem.f90 +++ b/src/setup/setup_solarsystem.f90 @@ -196,7 +196,7 @@ subroutine set_solarsystem_planets(nptmass,xyzmh_ptmass,vxyz_ptmass) 'saturn ', & 'uranus ', & 'neptune', & - 'pluto '/) ! for nostalgia's sake + 'apophis'/) ! 'moon ','apophis', 'pluto ' for nostalgia's sake real :: elems(nelem),xyz_tmp(size(xyzmh_ptmass(:,1)),2),vxyz_tmp(3,2),gm_cgs real :: msun,mplanet,a,e,inc,O,w,f integer :: i,ierr,ntmp diff --git a/src/utils/utils_ephemeris.f90 b/src/utils/utils_ephemeris.f90 index 81efb8325..ddf76634b 100644 --- a/src/utils/utils_ephemeris.f90 +++ b/src/utils/utils_ephemeris.f90 @@ -69,7 +69,7 @@ subroutine construct_horizons_api_url(object,url,ierr) character(len=*), intent(in) :: object ! name of the solar system object character(len=*), intent(out) :: url ! url for query integer, intent(out) :: ierr - character(len=3) :: cmd + character(len=6) :: cmd character(len=10) :: start_epoch,end_epoch integer :: values(8),year,month,day @@ -96,7 +96,7 @@ subroutine construct_horizons_api_url(object,url,ierr) case('moon') cmd = '301' ! moon case('apophis') - cmd = '2004 MN4' ! 99942 Apophis + cmd = '99942' ! 99942 Apophis case default cmd = trim(object) end select From 03b50a90b4e208795f83fda87da66bb78434aeed Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Tue, 10 Sep 2024 11:10:05 +1000 Subject: [PATCH 16/29] (tillotson) fix windows crap --- src/main/eos_tillotson.f90 | 422 ++++++++++++++++++------------------- 1 file changed, 211 insertions(+), 211 deletions(-) diff --git a/src/main/eos_tillotson.f90 b/src/main/eos_tillotson.f90 index 46ac7a42c..f2d6b0f71 100644 --- a/src/main/eos_tillotson.f90 +++ b/src/main/eos_tillotson.f90 @@ -1,212 +1,212 @@ -!--------------------------------------------------------------------------! -! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! -! Copyright (c) 2007-2024 The Authors (see AUTHORS) ! -! See LICENCE file for usage and distribution conditions ! -! http://phantomsph.github.io/ ! -!--------------------------------------------------------------------------! -module eos_tillotson -! -! EoS from Tillotson et al. (1962) -! -! :References: https://apps.dtic.mil/sti/pdfs/AD0486711.pdf -! -! Implementation from Benz and Asphaug (1999) -! -! :Owner: -! -! :Runtime parameters: None -! -! :Dependencies: -! - implicit none - real :: rho0 = 2.7 ! g/cm^3 zero-pressure density (Basalt) from Benz & Asphaug 1999 -! aparam, bparam, A, B, energy0 material-dependent Tillotson parameters (Basalt) - real :: aparam = 0.5 , bparam = 1.5 , alpha = 5. , beta = 5. - real :: A = 2.67e11 , B = 2.67e11 ! erg/cm^3 - real :: u0 = 4.87e12 ! erg/g - private :: rho0, aparam, bparam, A, B, alpha, beta, u0 - -contains -!----------------------------------------------------------------------- -!+ -! EoS from Tillotson et al. (1962) ; Implementation from Benz and Asphaug (1999) -!+ -!----------------------------------------------------------------------- -subroutine init_eos_tillotson(ierr) - integer, intent(out) :: ierr - - ierr = 0 - -end subroutine init_eos_tillotson - -subroutine equationofstate_tillotson(rho,u,pressure,spsound,gamma) - real, intent(in) :: rho,u - real, intent(out) :: pressure,spsound,gamma - real :: eta, mu, neu, omega, spsoundmin - -eta = rho / rho0 -mu = eta - 1. -neu = (1. / eta) - 1. ! Kegerreis et al. 2020 -omega = (u / (u0*eta**2) + 1.) - -if (rho >= rho0 .and. u >= 0.) then ! Condensed state - ! pressure = ( aparam + ( bparam / ( u / ( u0 * eta**2 ) + 1. ) ) ) * rho*u + A*mu + B*mu**2 - pressure = ( aparam + ( bparam / omega ) ) * rho*u + A*mu + B*mu**2 - if (pressure <= 0.) then - pressure = 0. - endif - ! sound speed from ! Kegereis et al. 2020 - spsound = sqrt(pressure / rho * (1. + aparam + (bparam / omega )) + & - (bparam*(omega - 1.) / omega**2) * (2.*u - pressure/rho)) + & - ((1./rho)*(A + B*((eta**2) -1.))) - spsoundmin = sqrt(A / rho) ! wave speed estimate - if (spsound < spsoundmin) then - spsound = spsoundmin - endif -else ! Expanded state - ! pressure = aparam*rho*u + ( (bparam*rho*u)/(u/(u0 * eta**2) + 1.) + & - ! A*mu*exp(-beta*( (rho0/rho) - 1.))) * exp(-alpha*( (rho0/rho) - 1.)**2) - pressure = aparam*rho*u + ( ((bparam*rho*u)/omega) + & - A*mu*exp(-beta*neu)) * exp(-alpha*(neu**2)) - if (pressure <= 0.) then - pressure = 0. - endif - spsound = sqrt( (pressure/rho)*( 1. + aparam + ((bparam / omega )*exp(-alpha*neu**2))) + & - ( ((bparam*rho*u)/((omega**2) * (eta**2))) * ((1/(u0*rho))*(2.*u - pressure/rho) & - + ((2.*alpha*eta*omega)/rho0)) + ((A/rho0)*(1.+ (mu/eta**2)*(beta*2.*alpha - eta))) & - * exp(-beta*neu) )*exp(-alpha*neu**2) ) - spsoundmin = sqrt(A / rho) ! wave speed estimate - if (spsound < spsoundmin) then - spsound = spsoundmin - endif -endif -! Hybrid Pressure (inbetween) -! uiv = -! ucv = -! frac = ((u / uiv) - 1.) / ((ucv/uiv) - 1.) -! pressure = -! ce = spsound @ expanded -! cc = spsound @ condensed -! spsound = sqrt(( ((u - uiv)* ce**2) + ((ucv - u)*cc**2))/(ucv - uiv)) - -end subroutine equationofstate_tillotson - -!---------------------------------------------------------------- -!+ -! print eos information -!+ -!---------------------------------------------------------------- -subroutine eos_info_tillotson(iprint) - integer, intent(in) :: iprint - - write(iprint,"(/,a)") ' Tillotson EoS' - -end subroutine eos_info_tillotson - -!----------------------------------------------------------------------- -!+ -! calc internal energy from pressure and density -!+ -!----------------------------------------------------------------------- -subroutine calc_uT_from_rhoP_tillotson(rho,pres,temp,u,ierr) - real, intent(in) :: rho,pres - real, intent(out) :: temp,u - integer, intent(out) :: ierr - real :: eta, mu, neu, omega, expterm - - eta = rho / rho0 - mu = eta - 1. - neu = (1. / eta) - 1. ! Kegerreis et al. 2020 - omega = (u / (u0*eta**2) + 1.) - - if (rho >= rho0) then ! Condensed state - u = (pres - A*mu - B*mu**2) / ( aparam + ( bparam / omega ) * rho) - else - expterm = exp(-alpha*(neu**2)) - u = (pres - A*mu*exp(-beta*neu)*expterm ) & - / (aparam*rho + ((bparam*rho)/omega)*expterm) - endif - temp = 300. - ierr = 0 - if (u<0.) then - u = 0. - ierr = 1 - endif - -end subroutine calc_uT_from_rhoP_tillotson -!----------------------------------------------------------------------- -!+ -! reads equation of state options from the input file -!+ -!----------------------------------------------------------------------- - subroutine read_options_eos_tillotson(name,valstring,imatch,igotall,ierr) - use io, only:fatal - character(len=*), intent(in) :: name,valstring - logical, intent(out) :: imatch,igotall - integer, intent(out) :: ierr - integer, save :: ngot = 0 - character(len=30), parameter :: label = 'eos_tillotson' - - imatch = .true. - select case(trim(name)) - case('rho0') - read(valstring,*,iostat=ierr) rho0 - if ((rho0 < 0.)) call fatal(label,'rho0 < 0') - ngot = ngot + 1 - case('aparam') - read(valstring,*,iostat=ierr) aparam - if ((aparam < 0.)) call fatal(label,'aparam < 0') - ngot = ngot + 1 - case('bparam') - read(valstring,*,iostat=ierr) bparam - if ((bparam < 0.)) call fatal(label,'bparam < 0') - ngot = ngot + 1 - case('A') - read(valstring,*,iostat=ierr) A - if ((A < 0.)) call fatal(label,'A < 0') - ngot = ngot + 1 - case('B') - read(valstring,*,iostat=ierr) B - if ((B < 0.)) call fatal(label,'B < 0') - ngot = ngot + 1 - case('alpha_t') - read(valstring,*,iostat=ierr) alpha - if ((alpha < 0.)) call fatal(label,'alpha < 0') - ngot = ngot + 1 - case('beta_t') - read(valstring,*,iostat=ierr) beta - if ((beta < 0.)) call fatal(label,'beta < 0') - ngot = ngot + 1 - case('u0') - read(valstring,*,iostat=ierr) u0 - if ((u0 < 0.)) call fatal(label,'u0 < 0') - ngot = ngot + 1 - case default - imatch = .false. - end select - - igotall = (ngot >= 8) - - end subroutine read_options_eos_tillotson - -!----------------------------------------------------------------------- -!+ -! writes equation of state options to the input file -!+ -!----------------------------------------------------------------------- -subroutine write_options_eos_tillotson(iunit) - use infile_utils, only:write_inopt - integer, intent(in) :: iunit - - call write_inopt(rho0,'rho0','reference density g/cm^3',iunit) - call write_inopt(aparam,'aparam','material-dependent Tillotson parameter, unitless',iunit) - call write_inopt(bparam,'bparam','material-dependent Tillotson parameter, unitless',iunit) - call write_inopt(A,'A','material-dependent Tillotson parameter, erg/cm^3',iunit) - call write_inopt(B,'B','material-dependent Tillotson parameter, erg/cm^3',iunit) - call write_inopt(alpha,'alpha_t','material-dependent Tillotson parameter, unitless',iunit) - call write_inopt(beta,'beta_t','material-dependent Tillotson parameter, unitless',iunit) - call write_inopt(u0,'u0','material-dependent Tillotson parameter, erg/g',iunit) - -end subroutine write_options_eos_tillotson - +!--------------------------------------------------------------------------! +! The Phantom Smoothed Particle Hydrodynamics code, by Daniel Price et al. ! +! Copyright (c) 2007-2024 The Authors (see AUTHORS) ! +! See LICENCE file for usage and distribution conditions ! +! http://phantomsph.github.io/ ! +!--------------------------------------------------------------------------! +module eos_tillotson +! +! EoS from Tillotson et al. (1962) +! +! :References: https://apps.dtic.mil/sti/pdfs/AD0486711.pdf +! +! Implementation from Benz and Asphaug (1999) +! +! :Owner: +! +! :Runtime parameters: None +! +! :Dependencies: +! + implicit none + real :: rho0 = 2.7 ! g/cm^3 zero-pressure density (Basalt) from Benz & Asphaug 1999 +! aparam, bparam, A, B, energy0 material-dependent Tillotson parameters (Basalt) + real :: aparam = 0.5 , bparam = 1.5 , alpha = 5. , beta = 5. + real :: A = 2.67e11 , B = 2.67e11 ! erg/cm^3 + real :: u0 = 4.87e12 ! erg/g + private :: rho0, aparam, bparam, A, B, alpha, beta, u0 + +contains +!----------------------------------------------------------------------- +!+ +! EoS from Tillotson et al. (1962) ; Implementation from Benz and Asphaug (1999) +!+ +!----------------------------------------------------------------------- +subroutine init_eos_tillotson(ierr) + integer, intent(out) :: ierr + + ierr = 0 + +end subroutine init_eos_tillotson + +subroutine equationofstate_tillotson(rho,u,pressure,spsound,gamma) + real, intent(in) :: rho,u + real, intent(out) :: pressure,spsound,gamma + real :: eta, mu, neu, omega, spsoundmin + +eta = rho / rho0 +mu = eta - 1. +neu = (1. / eta) - 1. ! Kegerreis et al. 2020 +omega = (u / (u0*eta**2) + 1.) + +if (rho >= rho0 .and. u >= 0.) then ! Condensed state + ! pressure = ( aparam + ( bparam / ( u / ( u0 * eta**2 ) + 1. ) ) ) * rho*u + A*mu + B*mu**2 + pressure = ( aparam + ( bparam / omega ) ) * rho*u + A*mu + B*mu**2 + if (pressure <= 0.) then + pressure = 0. + endif + ! sound speed from ! Kegereis et al. 2020 + spsound = sqrt(pressure / rho * (1. + aparam + (bparam / omega )) + & + (bparam*(omega - 1.) / omega**2) * (2.*u - pressure/rho)) + & + ((1./rho)*(A + B*((eta**2) -1.))) + spsoundmin = sqrt(A / rho) ! wave speed estimate + if (spsound < spsoundmin) then + spsound = spsoundmin + endif +else ! Expanded state + ! pressure = aparam*rho*u + ( (bparam*rho*u)/(u/(u0 * eta**2) + 1.) + & + ! A*mu*exp(-beta*( (rho0/rho) - 1.))) * exp(-alpha*( (rho0/rho) - 1.)**2) + pressure = aparam*rho*u + ( ((bparam*rho*u)/omega) + & + A*mu*exp(-beta*neu)) * exp(-alpha*(neu**2)) + if (pressure <= 0.) then + pressure = 0. + endif + spsound = sqrt( (pressure/rho)*( 1. + aparam + ((bparam / omega )*exp(-alpha*neu**2))) + & + ( ((bparam*rho*u)/((omega**2) * (eta**2))) * ((1/(u0*rho))*(2.*u - pressure/rho) & + + ((2.*alpha*eta*omega)/rho0)) + ((A/rho0)*(1.+ (mu/eta**2)*(beta*2.*alpha - eta))) & + * exp(-beta*neu) )*exp(-alpha*neu**2) ) + spsoundmin = sqrt(A / rho) ! wave speed estimate + if (spsound < spsoundmin) then + spsound = spsoundmin + endif +endif +! Hybrid Pressure (inbetween) +! uiv = +! ucv = +! frac = ((u / uiv) - 1.) / ((ucv/uiv) - 1.) +! pressure = +! ce = spsound @ expanded +! cc = spsound @ condensed +! spsound = sqrt(( ((u - uiv)* ce**2) + ((ucv - u)*cc**2))/(ucv - uiv)) + +end subroutine equationofstate_tillotson + +!---------------------------------------------------------------- +!+ +! print eos information +!+ +!---------------------------------------------------------------- +subroutine eos_info_tillotson(iprint) + integer, intent(in) :: iprint + + write(iprint,"(/,a)") ' Tillotson EoS' + +end subroutine eos_info_tillotson + +!----------------------------------------------------------------------- +!+ +! calc internal energy from pressure and density +!+ +!----------------------------------------------------------------------- +subroutine calc_uT_from_rhoP_tillotson(rho,pres,temp,u,ierr) + real, intent(in) :: rho,pres + real, intent(out) :: temp,u + integer, intent(out) :: ierr + real :: eta, mu, neu, omega, expterm + + eta = rho / rho0 + mu = eta - 1. + neu = (1. / eta) - 1. ! Kegerreis et al. 2020 + omega = (u / (u0*eta**2) + 1.) + + if (rho >= rho0) then ! Condensed state + u = (pres - A*mu - B*mu**2) / ( aparam + ( bparam / omega ) * rho) + else + expterm = exp(-alpha*(neu**2)) + u = (pres - A*mu*exp(-beta*neu)*expterm ) & + / (aparam*rho + ((bparam*rho)/omega)*expterm) + endif + temp = 300. + ierr = 0 + if (u<0.) then + u = 0. + ierr = 1 + endif + +end subroutine calc_uT_from_rhoP_tillotson +!----------------------------------------------------------------------- +!+ +! reads equation of state options from the input file +!+ +!----------------------------------------------------------------------- + subroutine read_options_eos_tillotson(name,valstring,imatch,igotall,ierr) + use io, only:fatal + character(len=*), intent(in) :: name,valstring + logical, intent(out) :: imatch,igotall + integer, intent(out) :: ierr + integer, save :: ngot = 0 + character(len=30), parameter :: label = 'eos_tillotson' + + imatch = .true. + select case(trim(name)) + case('rho0') + read(valstring,*,iostat=ierr) rho0 + if ((rho0 < 0.)) call fatal(label,'rho0 < 0') + ngot = ngot + 1 + case('aparam') + read(valstring,*,iostat=ierr) aparam + if ((aparam < 0.)) call fatal(label,'aparam < 0') + ngot = ngot + 1 + case('bparam') + read(valstring,*,iostat=ierr) bparam + if ((bparam < 0.)) call fatal(label,'bparam < 0') + ngot = ngot + 1 + case('A') + read(valstring,*,iostat=ierr) A + if ((A < 0.)) call fatal(label,'A < 0') + ngot = ngot + 1 + case('B') + read(valstring,*,iostat=ierr) B + if ((B < 0.)) call fatal(label,'B < 0') + ngot = ngot + 1 + case('alpha_t') + read(valstring,*,iostat=ierr) alpha + if ((alpha < 0.)) call fatal(label,'alpha < 0') + ngot = ngot + 1 + case('beta_t') + read(valstring,*,iostat=ierr) beta + if ((beta < 0.)) call fatal(label,'beta < 0') + ngot = ngot + 1 + case('u0') + read(valstring,*,iostat=ierr) u0 + if ((u0 < 0.)) call fatal(label,'u0 < 0') + ngot = ngot + 1 + case default + imatch = .false. + end select + + igotall = (ngot >= 8) + + end subroutine read_options_eos_tillotson + +!----------------------------------------------------------------------- +!+ +! writes equation of state options to the input file +!+ +!----------------------------------------------------------------------- +subroutine write_options_eos_tillotson(iunit) + use infile_utils, only:write_inopt + integer, intent(in) :: iunit + + call write_inopt(rho0,'rho0','reference density g/cm^3',iunit) + call write_inopt(aparam,'aparam','material-dependent Tillotson parameter, unitless',iunit) + call write_inopt(bparam,'bparam','material-dependent Tillotson parameter, unitless',iunit) + call write_inopt(A,'A','material-dependent Tillotson parameter, erg/cm^3',iunit) + call write_inopt(B,'B','material-dependent Tillotson parameter, erg/cm^3',iunit) + call write_inopt(alpha,'alpha_t','material-dependent Tillotson parameter, unitless',iunit) + call write_inopt(beta,'beta_t','material-dependent Tillotson parameter, unitless',iunit) + call write_inopt(u0,'u0','material-dependent Tillotson parameter, erg/g',iunit) + +end subroutine write_options_eos_tillotson + end module eos_tillotson \ No newline at end of file From db69327930641a16edd897573d1df55656d05918 Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Wed, 11 Sep 2024 16:26:52 +1000 Subject: [PATCH 17/29] (solarsystem) initial change from multiple binary setups to single star with particles setup --- src/setup/setup_solarsystem.f90 | 81 +++++++++++++++++++++++++++------ 1 file changed, 67 insertions(+), 14 deletions(-) diff --git a/src/setup/setup_solarsystem.f90 b/src/setup/setup_solarsystem.f90 index d42611364..683def261 100644 --- a/src/setup/setup_solarsystem.f90 +++ b/src/setup/setup_solarsystem.f90 @@ -19,32 +19,43 @@ module setup ! :Dependencies: centreofmass, datautils, ephemeris, infile_utils, io, mpc, ! part, physcon, setbinary, timestep, units ! + use setstar, only:star_t +! use setorbit, only:orbit_t +! use dim, only:gr implicit none public :: setpart real :: norbits integer :: dumpsperorbit + logical :: relax !,corotate + type(star_t) :: star +! type(orbit_t) :: orbit + private contains - !---------------------------------------------------------------- !+ ! setup for solar system orbits !+ !---------------------------------------------------------------- subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,time,fileprefix) - use part, only:nptmass,xyzmh_ptmass,vxyz_ptmass,idust,set_particle_type,& - grainsize,graindens,ndustlarge,ndusttypes,ihacc - use setbinary, only:set_binary - use units, only:set_units,umass,udist,unit_density - use physcon, only:solarm,au,pi,km,solarr - use io, only:master,fatal - use timestep, only:tmax,dtmax - use mpc, only:read_mpc,mpc_entry - use datautils, only:find_datafile - use centreofmass, only:reset_centreofmass + use part, only:nptmass,xyzmh_ptmass,vxyz_ptmass,idust,set_particle_type,& + grainsize,graindens,ndustlarge,ndusttypes,ihacc + use setorbit, only:set_defaults_orbit,set_orbit + use setstar, only:set_defaults_star,shift_star,set_star + use setbinary, only:set_binary + use options, only:iexternalforce + use externalforces, only:iext_corotate,iext_geopot,iext_star,omega_corotate,mass1,accradius1 + use units, only:set_units,umass,udist,unit_density + use physcon, only:solarm,au,pi,km,solarr + use io, only:master,fatal + use timestep, only:tmax,dtmax + use mpc, only:read_mpc,mpc_entry + use datautils, only:find_datafile + use centreofmass, only:reset_centreofmass +! use kernel, only:hfact_default integer, intent(in) :: id integer, intent(inout) :: npart integer, intent(out) :: npartoftype(:) @@ -61,6 +72,31 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, integer, parameter :: max_bodies = 2000000 type(mpc_entry), allocatable :: dat(:) ! +!--general parameters +! + dist_unit = 'solarr' + mass_unit = 'solarm' + time = 0. + polyk = 0. + gamma = 1. +! hfact = hfact_default +! +!--space available for injected gas particles +! in case only sink particles are used +! + npart = 0 + npartoftype(:) = 0 + massoftype = 0. + + xyzh(:,:) = 0. + vxyzu(:,:) = 0. + nptmass = 0 + nstar = 1 + call set_defaults_star(star) + call set_defaults_orbit(orbit) + relax = .true. + corotate = .false. +! ! default runtime parameters ! norbits = 1000. @@ -68,7 +104,9 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, ! ! read runtime parameters from setup file ! - if (id==master) print "(/,65('-'),1(/,a),/,65('-'),/)",' solar system' + if (id==master) print "(/,65('-'),1(/,a),/,65('-'),/)",& + ' Welcome to the Superb Solar System Setup' + filename = trim(fileprefix)//'.setup' inquire(file=filename,exist=iexist) if (iexist) call read_setupfile(filename,ierr) @@ -79,6 +117,20 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact, endif stop endif + ! + !--setup and relax stars as needed + ! + use_var_comp = .false. + write_profile = .false. + iextern_prev = iexternalforce + iexternalforce = 0 + gamma = 5./3. + call set_star(id,master,nstar,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) + + call reset_centreofmass(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptmass) ! ! set units ! @@ -186,17 +238,18 @@ subroutine set_solarsystem_planets(nptmass,xyzmh_ptmass,vxyz_ptmass) use setbinary, only:set_binary integer, intent(inout) :: nptmass real, intent(inout) :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:) - integer, parameter :: nplanets = 9 + integer, parameter :: nplanets = 10 ! including moon character(len=*), parameter :: planet_name(nplanets) = & (/'mercury', & 'venus ', & 'earth ', & + 'moon ', & 'mars ', & 'jupiter', & 'saturn ', & 'uranus ', & 'neptune', & - 'apophis'/) ! 'moon ','apophis', 'pluto ' for nostalgia's sake + 'pluto '/) ! for nostalgia's sake real :: elems(nelem),xyz_tmp(size(xyzmh_ptmass(:,1)),2),vxyz_tmp(3,2),gm_cgs real :: msun,mplanet,a,e,inc,O,w,f integer :: i,ierr,ntmp From a7ad58f8f525462117a8a9032c73c495cca08846 Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Fri, 13 Sep 2024 11:52:41 +1000 Subject: [PATCH 18/29] (tillotson) update eos to include hybrid phase and fix greek parameter errors --- src/main/eos_tillotson.f90 | 107 +++++++++++++++++++++++-------------- 1 file changed, 67 insertions(+), 40 deletions(-) diff --git a/src/main/eos_tillotson.f90 b/src/main/eos_tillotson.f90 index f2d6b0f71..fa1ccb98a 100644 --- a/src/main/eos_tillotson.f90 +++ b/src/main/eos_tillotson.f90 @@ -24,7 +24,11 @@ module eos_tillotson real :: aparam = 0.5 , bparam = 1.5 , alpha = 5. , beta = 5. real :: A = 2.67e11 , B = 2.67e11 ! erg/cm^3 real :: u0 = 4.87e12 ! erg/g - private :: rho0, aparam, bparam, A, B, alpha, beta, u0 + real :: u_iv = 4.72e10 ! erg/g + real :: u_cv = 1.82e11 ! erg/g + real :: rho_iv = 2.57 ! g/cm^3 incipient vaporisation density of gabbroic anorthosite lpp from Ahrens and Okeefe (1977) table 1 +! real :: rho_cv = 8.47 ! g/cm^3 complete vaporisation density of gabbroic anorthosite lpp from Ahrens and Okeefe (1977) table 2 + private :: rho0, aparam, bparam, A, B, alpha, beta, u0, rho_iv, u_iv, u_cv contains !----------------------------------------------------------------------- @@ -41,53 +45,76 @@ end subroutine init_eos_tillotson subroutine equationofstate_tillotson(rho,u,pressure,spsound,gamma) real, intent(in) :: rho,u - real, intent(out) :: pressure,spsound,gamma - real :: eta, mu, neu, omega, spsoundmin + real, intent(out) :: pressure, spsound, gamma + real :: eta, mu, neu, omega, spsoundmin, spsound2, spsound3, pressure2, pressure3 eta = rho / rho0 mu = eta - 1. neu = (1. / eta) - 1. ! Kegerreis et al. 2020 omega = (u / (u0*eta**2) + 1.) -if (rho >= rho0 .and. u >= 0.) then ! Condensed state - ! pressure = ( aparam + ( bparam / ( u / ( u0 * eta**2 ) + 1. ) ) ) * rho*u + A*mu + B*mu**2 - pressure = ( aparam + ( bparam / omega ) ) * rho*u + A*mu + B*mu**2 - if (pressure <= 0.) then - pressure = 0. - endif - ! sound speed from ! Kegereis et al. 2020 - spsound = sqrt(pressure / rho * (1. + aparam + (bparam / omega )) + & - (bparam*(omega - 1.) / omega**2) * (2.*u - pressure/rho)) + & - ((1./rho)*(A + B*((eta**2) -1.))) - spsoundmin = sqrt(A / rho) ! wave speed estimate - if (spsound < spsoundmin) then - spsound = spsoundmin - endif -else ! Expanded state - ! pressure = aparam*rho*u + ( (bparam*rho*u)/(u/(u0 * eta**2) + 1.) + & - ! A*mu*exp(-beta*( (rho0/rho) - 1.))) * exp(-alpha*( (rho0/rho) - 1.)**2) - pressure = aparam*rho*u + ( ((bparam*rho*u)/omega) + & - A*mu*exp(-beta*neu)) * exp(-alpha*(neu**2)) - if (pressure <= 0.) then - pressure = 0. - endif - spsound = sqrt( (pressure/rho)*( 1. + aparam + ((bparam / omega )*exp(-alpha*neu**2))) + & - ( ((bparam*rho*u)/((omega**2) * (eta**2))) * ((1/(u0*rho))*(2.*u - pressure/rho) & - + ((2.*alpha*eta*omega)/rho0)) + ((A/rho0)*(1.+ (mu/eta**2)*(beta*2.*alpha - eta))) & - * exp(-beta*neu) )*exp(-alpha*neu**2) ) - spsoundmin = sqrt(A / rho) ! wave speed estimate - if (spsound < spsoundmin) then - spsound = spsoundmin +! Brundage 2013 eq (2) +pressure2 = (( aparam + ( bparam / omega ) ) * rho*u) + A*mu + B*mu**2 + +! Brundage 2013 eq (3) +pressure3 = (aparam*rho*u) + ( ((bparam*rho*u)/omega) + & + A*mu*exp(-beta*neu)) * exp(-alpha*(neu**2)) + +! sound speed from Kegereis et al. 2020 +spsound2 = sqrt((pressure / rho) * (1. + aparam + (bparam / omega )) + & + ((bparam*(omega - 1.)) / omega**2) * (2.*u - pressure/rho)) + & + ((1./rho)*(A + B*((eta**2) - 1.))) + +spsound3 = sqrt( (pressure/rho)*( 1. + aparam + ((bparam / omega )*exp(-alpha*(neu**2)))) + & + ( ((bparam*rho*u)/((omega**2)*(eta**2))) * ((1. /(u0*rho))*((2.*u) - (pressure/rho)) & + + ((2.*alpha*neu*omega)/rho0)) + ((A/rho0)*(1.+ (mu/eta**2)*((beta*2.*alpha*neu) - eta))) & + * exp(-beta*neu) )*exp(-alpha*(neu**2)) ) + +if (rho < rho0) then + if (u >= u_cv) then ! completely vaporised + pressure = pressure3 + if (pressure <= 0.) then + pressure = 0. + endif + spsound = spsound3 + spsoundmin = sqrt(A / rho) ! wave speed estimate + if (spsound < spsoundmin) then + spsound = spsoundmin + endif + elseif (rho >= rho_iv) then + if (u <= u_iv) then ! cold expanded + pressure = pressure2 + if (pressure <= 0.) then + pressure = 0. + endif + spsound = spsound2 + spsoundmin = sqrt(A / rho) ! wave speed estimate + if (spsound < spsoundmin) then + spsound = spsoundmin + endif + elseif (u < u_cv) then ! hybrid + pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) + if (pressure <= 0.) then + pressure = 0. + endif + spsound = sqrt(( ((u-u_iv)*(spsound3)**2) + ((u_cv-u)*(spsound2)**2) ) / (u_cv - u_iv)) + spsoundmin = sqrt(A / rho) ! wave speed estimate + if (spsound < spsoundmin) then + spsound = spsoundmin + endif + endif endif +else ! rho >= rho0 compressed state + pressure = pressure3 + if (pressure <= 0.) then + pressure = 0. + endif + spsound = spsound3 + spsoundmin = sqrt(A / rho) ! wave speed estimate + if (spsound < spsoundmin) then + spsound = spsoundmin + endif endif -! Hybrid Pressure (inbetween) -! uiv = -! ucv = -! frac = ((u / uiv) - 1.) / ((ucv/uiv) - 1.) -! pressure = -! ce = spsound @ expanded -! cc = spsound @ condensed -! spsound = sqrt(( ((u - uiv)* ce**2) + ((ucv - u)*cc**2))/(ucv - uiv)) end subroutine equationofstate_tillotson From 40d7cb2847abf45a7c5f2932775afc5446dd313c Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Fri, 13 Sep 2024 15:51:49 +1000 Subject: [PATCH 19/29] (tillotson) updated internal energy calculation --- src/main/eos_tillotson.f90 | 155 +++++++++++++++++++------------------ 1 file changed, 79 insertions(+), 76 deletions(-) diff --git a/src/main/eos_tillotson.f90 b/src/main/eos_tillotson.f90 index fa1ccb98a..173ac54b1 100644 --- a/src/main/eos_tillotson.f90 +++ b/src/main/eos_tillotson.f90 @@ -29,6 +29,8 @@ module eos_tillotson real :: rho_iv = 2.57 ! g/cm^3 incipient vaporisation density of gabbroic anorthosite lpp from Ahrens and Okeefe (1977) table 1 ! real :: rho_cv = 8.47 ! g/cm^3 complete vaporisation density of gabbroic anorthosite lpp from Ahrens and Okeefe (1977) table 2 private :: rho0, aparam, bparam, A, B, alpha, beta, u0, rho_iv, u_iv, u_cv +! private :: spsoundmin, spsound2, spsound3, pressure2, pressure3 +! private :: expterm, expterm2, sqrtterm, sqrtterm1, sqrtterm2 contains !----------------------------------------------------------------------- @@ -48,73 +50,62 @@ subroutine equationofstate_tillotson(rho,u,pressure,spsound,gamma) real, intent(out) :: pressure, spsound, gamma real :: eta, mu, neu, omega, spsoundmin, spsound2, spsound3, pressure2, pressure3 -eta = rho / rho0 -mu = eta - 1. -neu = (1. / eta) - 1. ! Kegerreis et al. 2020 -omega = (u / (u0*eta**2) + 1.) + eta = rho / rho0 + mu = eta - 1. + neu = (1. / eta) - 1. ! Kegerreis et al. 2020 + omega = (u / (u0*eta**2) + 1.) + spsoundmin = sqrt(A / rho) ! wave speed estimate +! gamma = 2./3. -! Brundage 2013 eq (2) -pressure2 = (( aparam + ( bparam / omega ) ) * rho*u) + A*mu + B*mu**2 +! Brundage 2013 eq (2) cold expanded and compressed states + pressure2 = (( aparam + ( bparam / omega ) ) * rho*u) + A*mu + B*mu**2 + if (pressure2 <= 0.) then + pressure2 = 0. + endif -! Brundage 2013 eq (3) -pressure3 = (aparam*rho*u) + ( ((bparam*rho*u)/omega) + & +! Brundage 2013 eq (3) completely vaporised + pressure3 = (aparam*rho*u) + ( ((bparam*rho*u)/omega) + & A*mu*exp(-beta*neu)) * exp(-alpha*(neu**2)) + if (pressure3 <= 0.) then + pressure3 = 0. + endif ! sound speed from Kegereis et al. 2020 -spsound2 = sqrt((pressure / rho) * (1. + aparam + (bparam / omega )) + & - ((bparam*(omega - 1.)) / omega**2) * (2.*u - pressure/rho)) + & + spsound2 = sqrt((pressure2 / rho) * (1. + aparam + (bparam / omega )) + & + ((bparam*(omega - 1.)) / omega**2) * (2.*u - pressure2/rho)) + & ((1./rho)*(A + B*((eta**2) - 1.))) + if (spsound2 < spsoundmin) then + spsound2 = spsoundmin + endif -spsound3 = sqrt( (pressure/rho)*( 1. + aparam + ((bparam / omega )*exp(-alpha*(neu**2)))) + & - ( ((bparam*rho*u)/((omega**2)*(eta**2))) * ((1. /(u0*rho))*((2.*u) - (pressure/rho)) & - + ((2.*alpha*neu*omega)/rho0)) + ((A/rho0)*(1.+ (mu/eta**2)*((beta*2.*alpha*neu) - eta))) & + spsound3 = sqrt( (pressure3/rho)*( 1. + aparam + ((bparam / omega )*exp(-alpha*(neu**2)))) + & + ( ((bparam*rho*u)/((omega**2)*(eta**2))) * ((1. /(u0*rho))*((2.*u) - (pressure3/rho)) & + + ((2.*alpha*neu*omega)/rho0)) + ((A/rho0)*(1. + (mu/eta**2)*((beta*2.*alpha*neu) - eta))) & * exp(-beta*neu) )*exp(-alpha*(neu**2)) ) - -if (rho < rho0) then - if (u >= u_cv) then ! completely vaporised - pressure = pressure3 - if (pressure <= 0.) then - pressure = 0. - endif - spsound = spsound3 - spsoundmin = sqrt(A / rho) ! wave speed estimate - if (spsound < spsoundmin) then - spsound = spsoundmin - endif - elseif (rho >= rho_iv) then - if (u <= u_iv) then ! cold expanded - pressure = pressure2 - if (pressure <= 0.) then - pressure = 0. - endif - spsound = spsound2 - spsoundmin = sqrt(A / rho) ! wave speed estimate - if (spsound < spsoundmin) then - spsound = spsoundmin - endif - elseif (u < u_cv) then ! hybrid - pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) - if (pressure <= 0.) then - pressure = 0. - endif - spsound = sqrt(( ((u-u_iv)*(spsound3)**2) + ((u_cv-u)*(spsound2)**2) ) / (u_cv - u_iv)) - spsoundmin = sqrt(A / rho) ! wave speed estimate - if (spsound < spsoundmin) then - spsound = spsoundmin - endif - endif - endif -else ! rho >= rho0 compressed state - pressure = pressure3 - if (pressure <= 0.) then - pressure = 0. - endif - spsound = spsound3 - spsoundmin = sqrt(A / rho) ! wave speed estimate - if (spsound < spsoundmin) then - spsound = spsoundmin - endif -endif + if (spsound3 < spsoundmin) then + spsound3 = spsoundmin + endif + + + if (rho < rho0) then + if (rho <= rho_iv) then ! u <= u_iv cold expanded + pressure = pressure2 + spsound = spsound2 + ! print*, "cold expanded, pressure = ", pressure + elseif (rho > rho_iv) then ! u_iv < u < u_cv hybrid + pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) + spsound = sqrt(( ((u-u_iv)*(spsound3)**2) + ((u_cv-u)*(spsound2)**2) ) / (u_cv - u_iv)) + ! print*, "hybrid, pressure = ", pressure + else ! u >= u_cv completely vaporised + pressure = pressure3 + spsound = spsound3 + ! print*, "vaporised, pressure = ", pressure + endif + else ! rho >= rho0 ; u >= 0. compressed state + pressure = pressure2 + spsound = spsound2 + ! print*, "compressed, pressure = ", pressure + endif end subroutine equationofstate_tillotson @@ -136,29 +127,41 @@ end subroutine eos_info_tillotson !+ !----------------------------------------------------------------------- subroutine calc_uT_from_rhoP_tillotson(rho,pres,temp,u,ierr) - real, intent(in) :: rho,pres - real, intent(out) :: temp,u - integer, intent(out) :: ierr - real :: eta, mu, neu, omega, expterm + real, intent(in) :: rho,pres + real, intent(out) :: temp,u + integer, intent(out) :: ierr + real :: eta, mu, neu, omega, expterm, expterm2, sqrtterm1, sqrtterm2, sqrtterm eta = rho / rho0 mu = eta - 1. neu = (1. / eta) - 1. ! Kegerreis et al. 2020 omega = (u / (u0*eta**2) + 1.) - if (rho >= rho0) then ! Condensed state - u = (pres - A*mu - B*mu**2) / ( aparam + ( bparam / omega ) * rho) - else - expterm = exp(-alpha*(neu**2)) - u = (pres - A*mu*exp(-beta*neu)*expterm ) & - / (aparam*rho + ((bparam*rho)/omega)*expterm) - endif - temp = 300. - ierr = 0 - if (u<0.) then - u = 0. - ierr = 1 - endif + if (rho >= rho0) then ! Condensed state + ! u = (pres - A*mu - B*mu**2) / ( aparam + ( bparam / omega ) * rho) + sqrtterm1 = ((aparam*(eta**2)*u0 + bparam*(eta**2)*u0 & + + (A*mu + B*(mu**2) - pres )/rho )**2 )/(4.*aparam**2) + sqrtterm2 = ((A*mu*(eta**2)*u0 + B*mu*(eta**2)*u0 & + - (eta**2)*u0*pres) / rho) / aparam + sqrtterm = sqrtterm1 - sqrtterm2 + u = abs(sqrt(sqrtterm)) - ( (aparam*(eta**2)*u0 + bparam*(eta**2)*u0 & + + ((A*mu + B*(mu**2) - pres) / rho)) / 2.*aparam) + else + expterm = exp(-alpha*(neu**2)) + expterm2 = exp(-beta*neu) + ! u = (pres - A*mu*exp(-beta*neu)*expterm ) & + ! / (aparam*rho + ((bparam*rho)/omega)*expterm) + sqrtterm = u0*(eta**2)*((pres - ((A*mu*expterm2)*expterm))/rho) + u = abs(sqrt(sqrtterm)) - (((aparam*(eta**2)*u0 + bparam*(eta**2)*u0*expterm) & + - (((pres/rho) + (A*mu*expterm2)/rho)*expterm)) / 2.*aparam) + endif + + temp = 300. + ierr = 0 + if (u < 0.) then + u = 0. + ierr = 1 + endif end subroutine calc_uT_from_rhoP_tillotson !----------------------------------------------------------------------- From fa0b709fb4711d38617428e56b8243c35749c168 Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Mon, 16 Sep 2024 11:56:06 +1000 Subject: [PATCH 20/29] (tillotson) updated pressure and energy conditional calculations --- src/main/eos_tillotson.f90 | 92 +++++++++++++++++++++----------------- 1 file changed, 50 insertions(+), 42 deletions(-) diff --git a/src/main/eos_tillotson.f90 b/src/main/eos_tillotson.f90 index 173ac54b1..ac06fd135 100644 --- a/src/main/eos_tillotson.f90 +++ b/src/main/eos_tillotson.f90 @@ -35,7 +35,7 @@ module eos_tillotson contains !----------------------------------------------------------------------- !+ -! EoS from Tillotson et al. (1962) ; Implementation from Benz and Asphaug (1999) +! EoS from Tillotson et al. (1962) ; Implementation from Benz and Asphaug (1999) and Brundage (2013) !+ !----------------------------------------------------------------------- subroutine init_eos_tillotson(ierr) @@ -48,7 +48,7 @@ end subroutine init_eos_tillotson subroutine equationofstate_tillotson(rho,u,pressure,spsound,gamma) real, intent(in) :: rho,u real, intent(out) :: pressure, spsound, gamma - real :: eta, mu, neu, omega, spsoundmin, spsound2, spsound3, pressure2, pressure3 + real :: eta, mu, neu, omega, spsoundmin, spsound2, spsound2h, spsound3, pressure2, pressure3 eta = rho / rho0 mu = eta - 1. @@ -71,9 +71,9 @@ subroutine equationofstate_tillotson(rho,u,pressure,spsound,gamma) endif ! sound speed from Kegereis et al. 2020 - spsound2 = sqrt((pressure2 / rho) * (1. + aparam + (bparam / omega )) + & - ((bparam*(omega - 1.)) / omega**2) * (2.*u - pressure2/rho)) + & - ((1./rho)*(A + B*((eta**2) - 1.))) + spsound2 = sqrt( ((pressure2 / rho) * (1. + aparam + (bparam / omega ))) + & + (((bparam*(omega - 1.)) / (omega**2)) * (2.*u - pressure2/rho)) + & + ((1./rho)*(A + B*((eta**2) - 1.)))) if (spsound2 < spsoundmin) then spsound2 = spsoundmin endif @@ -86,25 +86,35 @@ subroutine equationofstate_tillotson(rho,u,pressure,spsound,gamma) spsound3 = spsoundmin endif - - if (rho < rho0) then - if (rho <= rho_iv) then ! u <= u_iv cold expanded - pressure = pressure2 - spsound = spsound2 - ! print*, "cold expanded, pressure = ", pressure - elseif (rho > rho_iv) then ! u_iv < u < u_cv hybrid - pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) - spsound = sqrt(( ((u-u_iv)*(spsound3)**2) + ((u_cv-u)*(spsound2)**2) ) / (u_cv - u_iv)) - ! print*, "hybrid, pressure = ", pressure - else ! u >= u_cv completely vaporised - pressure = pressure3 - spsound = spsound3 - ! print*, "vaporised, pressure = ", pressure - endif - else ! rho >= rho0 ; u >= 0. compressed state - pressure = pressure2 - spsound = spsound2 - ! print*, "compressed, pressure = ", pressure + spsound2h = sqrt( ((pressure2 / rho) * (1. + aparam + (bparam / omega ))) + & + (((bparam*(omega - 1.)) / (omega**2)) * (2.*u - pressure2/rho)) + & + ((1./rho)*(A + B*((eta**2) - 1.))) - ((((2.*B)/rho0**2)*rho) - (2.*B/rho0))) + if (spsound2h < spsoundmin) then + spsound2h = spsoundmin + endif + + if (rho >= rho0) then ! compressed state + pressure = pressure2 + spsound = spsound2 + else ! rho < rho0 + if (u >= u_cv) then ! vaporised expanded state + pressure = pressure3 + spsound = spsound3 + ! need to check rho ~0 and u >> u_cv (zero density becomes pressure = rho*gamma*u fermi, gamma = 2/3) + else ! u < u_cv + if (rho >= rho_iv) then + if (u > u_iv) then ! Hybrid state + pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) + spsound = sqrt(( ((u-u_iv)*(spsound3)**2) + ((u_cv-u)*(spsound2)**2) ) / (u_cv - u_iv)) + else ! u <= u_iv cold expanded + pressure = pressure2 + spsound = spsound2 + endif + else ! rho < rho_iv Low energy state + pressure = pressure2 - (B*mu**2) + spsound = spsound2h + endif + endif endif end subroutine equationofstate_tillotson @@ -136,24 +146,22 @@ subroutine calc_uT_from_rhoP_tillotson(rho,pres,temp,u,ierr) mu = eta - 1. neu = (1. / eta) - 1. ! Kegerreis et al. 2020 omega = (u / (u0*eta**2) + 1.) - - if (rho >= rho0) then ! Condensed state - ! u = (pres - A*mu - B*mu**2) / ( aparam + ( bparam / omega ) * rho) - sqrtterm1 = ((aparam*(eta**2)*u0 + bparam*(eta**2)*u0 & - + (A*mu + B*(mu**2) - pres )/rho )**2 )/(4.*aparam**2) - sqrtterm2 = ((A*mu*(eta**2)*u0 + B*mu*(eta**2)*u0 & - - (eta**2)*u0*pres) / rho) / aparam - sqrtterm = sqrtterm1 - sqrtterm2 - u = abs(sqrt(sqrtterm)) - ( (aparam*(eta**2)*u0 + bparam*(eta**2)*u0 & - + ((A*mu + B*(mu**2) - pres) / rho)) / 2.*aparam) - else - expterm = exp(-alpha*(neu**2)) - expterm2 = exp(-beta*neu) - ! u = (pres - A*mu*exp(-beta*neu)*expterm ) & - ! / (aparam*rho + ((bparam*rho)/omega)*expterm) - sqrtterm = u0*(eta**2)*((pres - ((A*mu*expterm2)*expterm))/rho) - u = abs(sqrt(sqrtterm)) - (((aparam*(eta**2)*u0 + bparam*(eta**2)*u0*expterm) & - - (((pres/rho) + (A*mu*expterm2)/rho)*expterm)) / 2.*aparam) + expterm2 = exp(-alpha*(neu**2)) + expterm = exp(-beta*neu) + + if (rho >= rho0) then ! compressed state + u = (mu*(eta**2)*u0*(aparam-bparam)*(A+(B*mu))) / & + (aparam*(A*mu+bparam*(eta**2)*u0*rho+B*(mu**2)-pres)) + elseif (rho >= rho_iv) then ! cold expanded + u = (mu*(eta**2)*u0*(aparam-bparam)*(A+(B*mu))) / & + (aparam*(A*mu+bparam*(eta**2)*u0*rho+B*(mu**2)-pres)) + elseif (rho < rho_iv) then ! rho < rho_iv Low energy state + u = (A*mu*(eta**2)*u0*(aparam-bparam)) / & + (aparam*(A*mu+bparam*(eta**2)*u0*rho-pres)) + else ! rho < rho0 vaporised expanded state + u = ((eta**2)*u0*(bparam-aparam)*(pres-A*mu*expterm*expterm2) ) / & + (aparam*rho*(A*mu*expterm + bparam*(eta**2)*u0 - pres)) + ! need to check rho ~0 and u >> u_cv (zero density becomes u = pres / (rho*gamma) fermi, gamma = 2/3) endif temp = 300. From 4eec9e74354a2108a2caafd55b7a059124a01414 Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Thu, 19 Sep 2024 08:31:30 +1000 Subject: [PATCH 21/29] (tillotson) quadratic equation fix for internal energy subroutine --- src/main/eos_tillotson.f90 | 115 ++++++++++++++++++++++++++++++++----- 1 file changed, 102 insertions(+), 13 deletions(-) diff --git a/src/main/eos_tillotson.f90 b/src/main/eos_tillotson.f90 index ac06fd135..551bff2df 100644 --- a/src/main/eos_tillotson.f90 +++ b/src/main/eos_tillotson.f90 @@ -29,6 +29,7 @@ module eos_tillotson real :: rho_iv = 2.57 ! g/cm^3 incipient vaporisation density of gabbroic anorthosite lpp from Ahrens and Okeefe (1977) table 1 ! real :: rho_cv = 8.47 ! g/cm^3 complete vaporisation density of gabbroic anorthosite lpp from Ahrens and Okeefe (1977) table 2 private :: rho0, aparam, bparam, A, B, alpha, beta, u0, rho_iv, u_iv, u_cv +! private :: quadratic_eq ! private :: spsoundmin, spsound2, spsound3, pressure2, pressure3 ! private :: expterm, expterm2, sqrtterm, sqrtterm1, sqrtterm2 @@ -53,7 +54,7 @@ subroutine equationofstate_tillotson(rho,u,pressure,spsound,gamma) eta = rho / rho0 mu = eta - 1. neu = (1. / eta) - 1. ! Kegerreis et al. 2020 - omega = (u / (u0*eta**2) + 1.) + omega = (u / (u0*eta**2)) + 1. spsoundmin = sqrt(A / rho) ! wave speed estimate ! gamma = 2./3. @@ -137,10 +138,13 @@ end subroutine eos_info_tillotson !+ !----------------------------------------------------------------------- subroutine calc_uT_from_rhoP_tillotson(rho,pres,temp,u,ierr) + use cubic, only:cubicsolve real, intent(in) :: rho,pres real, intent(out) :: temp,u integer, intent(out) :: ierr - real :: eta, mu, neu, omega, expterm, expterm2, sqrtterm1, sqrtterm2, sqrtterm + real :: eta, mu, neu, omega, expterm, expterm2, u1, u2, X, Y, quada, quadb, quadc + real :: usol(3) + integer :: nsol eta = rho / rho0 mu = eta - 1. @@ -150,28 +154,113 @@ subroutine calc_uT_from_rhoP_tillotson(rho,pres,temp,u,ierr) expterm = exp(-beta*neu) if (rho >= rho0) then ! compressed state - u = (mu*(eta**2)*u0*(aparam-bparam)*(A+(B*mu))) / & - (aparam*(A*mu+bparam*(eta**2)*u0*rho+B*(mu**2)-pres)) + quada = 1. + X = (pres - (A*mu) - (B*(mu**2)))/rho + quadb = -((X-(aparam*u0*(eta**2))-(bparam*u0*(eta**2))) / aparam) + quadc = -((u0*(eta**2)*X)/aparam) + call cubicsolve(0.,quada,quadb,quadc,usol,nsol) + ! call quadratic_eq(quada,quadb,quadc,u1,u2) + ! if (u2 <= 0.) then + ! u2 = 0. + ! elseif (u2 == u1) then + ! u = u1 + ! elseif (u2 > u1) then + ! u = u1 + ! endif + + ! u = (((1./aparam)-(1./bparam))*(-A*mu - B*(mu**2))/ & + ! (rho*(1. - ((-A*mu - B*(mu**2)+pres)/(bparam*(eta**2)*u0*rho))))) + ! u = (mu*(eta**2)*u0*(aparam-bparam)*(A+(B*mu))) / & + ! (aparam*(A*mu+bparam*(eta**2)*u0*rho+B*(mu**2)-pres)) elseif (rho >= rho_iv) then ! cold expanded - u = (mu*(eta**2)*u0*(aparam-bparam)*(A+(B*mu))) / & - (aparam*(A*mu+bparam*(eta**2)*u0*rho+B*(mu**2)-pres)) + quada = 1. + X = (pres - (A*mu) - (B*(mu**2)))/rho + quadb = -((X-(aparam*u0*(eta**2))-(bparam*u0*(eta**2))) / aparam) + quadc = -((u0*(eta**2)*X)/aparam) + call cubicsolve(0.,quada,quadb,quadc,usol,nsol) + ! call quadratic_eq(quada,quadb,quadc,u1,u2) + ! if (u2 <= 0.) then + ! u2 = 0. + ! elseif (u2 == u1) then + ! u = u1 + ! elseif (u2 > u1) then + ! u = u1 + ! endif + ! u = (((1./aparam)-(1./bparam))*(-A*mu - B*(mu**2))/ & + ! (rho*(1. - ((-A*mu - B*(mu**2)+pres)/(bparam*(eta**2)*u0*rho))))) + ! u = (mu*(eta**2)*u0*(aparam-bparam)*(A+(B*mu))) / & + ! (aparam*(A*mu+bparam*(eta**2)*u0*rho+B*(mu**2)-pres)) elseif (rho < rho_iv) then ! rho < rho_iv Low energy state - u = (A*mu*(eta**2)*u0*(aparam-bparam)) / & - (aparam*(A*mu+bparam*(eta**2)*u0*rho-pres)) + quada = 1. + X = (pres - (A*mu))/rho + quadb = -((X-(aparam*u0*(eta**2))-(bparam*u0*(eta**2))) / aparam) + quadc = -((u0*(eta**2)*X)/aparam) + call cubicsolve(0.,quada,quadb,quadc,usol,nsol) + ! call quadratic_eq(quada,quadb,quadc,u1,u2) + ! if (u2 <= 0.) then + ! u2 = 0. + ! elseif (u2 == u1) then + ! u = u1 + ! elseif (u2 > u1) then + ! u = u1 + ! endif + ! u = ((-A*mu*((1./aparam)-(1./bparam))) /& + ! (rho*(1.-((pres-A*mu)/(bparam*(eta**2)*u0*rho))))) + ! u = (A*mu*(eta**2)*u0*(aparam-bparam)) / & + ! (aparam*(A*mu+bparam*(eta**2)*u0*rho-pres)) else ! rho < rho0 vaporised expanded state - u = ((eta**2)*u0*(bparam-aparam)*(pres-A*mu*expterm*expterm2) ) / & - (aparam*rho*(A*mu*expterm + bparam*(eta**2)*u0 - pres)) + quada = 1. + X = (aparam*rho)/(expterm2) + quadb = ((bparam*rho*u0*(eta**2))+((aparam*rho*u0*(eta**2))/expterm2)+(A*mu*expterm)-(pres/(expterm2)))/X + quadc = ((A*mu*expterm*u0*(eta**2)) -((pres*u0*(eta**2))/expterm2))/X + call cubicsolve(0.,quada,quadb,quadc,usol,nsol) + !call quadratic_eq(quada,quadb,quadc,u1,u2) + u = maxval(usol(1:nsol)) + ! u = ((eta**2)*u0*(bparam-aparam)*(pres-A*mu*expterm*expterm2) ) / & + ! (aparam*rho*(A*mu*expterm + bparam*(eta**2)*u0 - pres)) ! need to check rho ~0 and u >> u_cv (zero density becomes u = pres / (rho*gamma) fermi, gamma = 2/3) endif temp = 300. ierr = 0 - if (u < 0.) then - u = 0. - ierr = 1 + if (u <= 0.) then + u = 0. + ierr = 1 endif end subroutine calc_uT_from_rhoP_tillotson + +!----------------------------------------------------------------------- +!+ +! calc quadratic for internal energy +!+ +!----------------------------------------------------------------------- +! function quadratic_eq(a,b,c,x1,x2) +! real, intent(in) :: a,b,c +! real, intent(out) :: x1,x2 +! real :: dis, sqrt_val + +! dis = (b**2) - (4. * a * c) ! calc discriminant using formula +! sqrt_val = sqrt(abs(dis)) + +! if (dis > 0) then ! check discriminant +! ! real and different roots +! x1 = (-b + sqrt_val)/(2 * a) +! x2 = (-b - sqrt_val)/(2 * a) + +! elseif (dis == 0) then +! ! real and same roots +! x1 = -b / (2. * a) +! x2 = x1 + +! else ! if discriminant < 0 + ! Complex Roots + ! print(- b / (2 * a), + i, sqrt_val / (2 * a)) + ! print(- b / (2 * a), - i, sqrt_val / (2 * a)) +! endif + +! end function quadratic_eq + !----------------------------------------------------------------------- !+ ! reads equation of state options from the input file From 211829a0fa8786b9b5b35a5814959ae307ad6a58 Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Thu, 19 Sep 2024 08:34:32 +1000 Subject: [PATCH 22/29] (tillotson) in Makefile changed order for reading in cubicsolve --- build/Makefile | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/build/Makefile b/build/Makefile index 8de3ebe81..ba1a86987 100644 --- a/build/Makefile +++ b/build/Makefile @@ -524,7 +524,7 @@ endif SOURCES= physcon.f90 ${CONFIG} ${SRCKERNEL} io.F90 units.f90 \ mpi_dens.F90 mpi_force.F90 mpi_utils.F90 dtype_kdtree.F90 utils_omp.F90 utils_cpuinfo.f90 \ - utils_infiles.f90 utils_allocate.f90 icosahedron.f90 quartic.f90 \ + utils_infiles.f90 utils_allocate.f90 icosahedron.f90 cubicsolve.f90 quartic.f90 \ utils_system.f90 utils_mathfunc.f90 part.F90 ${DOMAIN} boundary.f90 boundary_dynamic.f90 utils_timing.f90 mpi_balance.F90 \ setup_params.f90 timestep.f90 utils_dumpfiles.f90 utils_indtimesteps.F90 \ utils_sort.f90 utils_supertimestep.F90 utils_tables.f90 utils_gravwave.f90 \ @@ -614,7 +614,7 @@ SRCDUMP= physcon.f90 ${CONFIG} ${SRCKERNEL} utils_omp.F90 io.F90 units.f90 mpi_u utils_timing.f90 utils_infiles.f90 dtype_kdtree.f90 utils_allocate.f90 part.F90 \ ${DOMAIN} mpi_dens.F90 mpi_force.F90 boundary.f90 boundary_dynamic.f90 \ mpi_balance.F90 mpi_memory.f90 mpi_derivs.F90 mpi_tree.F90 kdtree.F90 linklist_kdtree.F90 \ - utils_dumpfiles.f90 utils_vectors.f90 utils_mathfunc.f90 \ + utils_dumpfiles.f90 utils_vectors.f90 utils_mathfunc.f90 cubicsolve.f90 \ utils_datafiles.f90 utils_filenames.f90 utils_system.f90 utils_tables.f90 datafiles.f90 gitinfo.f90 \ centreofmass.f90 \ timestep.f90 ${SRCEOS} cullendehnen.f90 dust_formation.f90 \ From 3bf2555d9e1448bd42e7f498dc5d58eff33ceb2c Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Fri, 20 Sep 2024 16:00:39 +1000 Subject: [PATCH 23/29] (tillotson) update piecewise pressure and internal energy --- src/main/eos_tillotson.f90 | 178 +++++++++++++++---------------------- 1 file changed, 70 insertions(+), 108 deletions(-) diff --git a/src/main/eos_tillotson.f90 b/src/main/eos_tillotson.f90 index 551bff2df..a275ce694 100644 --- a/src/main/eos_tillotson.f90 +++ b/src/main/eos_tillotson.f90 @@ -21,6 +21,7 @@ module eos_tillotson implicit none real :: rho0 = 2.7 ! g/cm^3 zero-pressure density (Basalt) from Benz & Asphaug 1999 ! aparam, bparam, A, B, energy0 material-dependent Tillotson parameters (Basalt) + real :: pressure = 0. ! THIS COULD WORK real :: aparam = 0.5 , bparam = 1.5 , alpha = 5. , beta = 5. real :: A = 2.67e11 , B = 2.67e11 ! erg/cm^3 real :: u0 = 4.87e12 ! erg/g @@ -29,7 +30,6 @@ module eos_tillotson real :: rho_iv = 2.57 ! g/cm^3 incipient vaporisation density of gabbroic anorthosite lpp from Ahrens and Okeefe (1977) table 1 ! real :: rho_cv = 8.47 ! g/cm^3 complete vaporisation density of gabbroic anorthosite lpp from Ahrens and Okeefe (1977) table 2 private :: rho0, aparam, bparam, A, B, alpha, beta, u0, rho_iv, u_iv, u_cv -! private :: quadratic_eq ! private :: spsoundmin, spsound2, spsound3, pressure2, pressure3 ! private :: expterm, expterm2, sqrtterm, sqrtterm1, sqrtterm2 @@ -49,24 +49,24 @@ end subroutine init_eos_tillotson subroutine equationofstate_tillotson(rho,u,pressure,spsound,gamma) real, intent(in) :: rho,u real, intent(out) :: pressure, spsound, gamma - real :: eta, mu, neu, omega, spsoundmin, spsound2, spsound2h, spsound3, pressure2, pressure3 + real :: eta, mu_t, neu, omega, spsoundmin, spsound2, spsoundh, spsound2h, spsound3, pressure2, pressure3 eta = rho / rho0 - mu = eta - 1. + mu_t = eta - 1. neu = (1. / eta) - 1. ! Kegerreis et al. 2020 omega = (u / (u0*eta**2)) + 1. spsoundmin = sqrt(A / rho) ! wave speed estimate ! gamma = 2./3. ! Brundage 2013 eq (2) cold expanded and compressed states - pressure2 = (( aparam + ( bparam / omega ) ) * rho*u) + A*mu + B*mu**2 + pressure2 = (( aparam + ( bparam / omega ) ) * rho*u) + A*mu_t + B*mu_t**2 if (pressure2 <= 0.) then pressure2 = 0. endif ! Brundage 2013 eq (3) completely vaporised pressure3 = (aparam*rho*u) + ( ((bparam*rho*u)/omega) + & - A*mu*exp(-beta*neu)) * exp(-alpha*(neu**2)) + A*mu_t*exp(-beta*neu)) * exp(-alpha*(neu**2)) if (pressure3 <= 0.) then pressure3 = 0. endif @@ -81,7 +81,7 @@ subroutine equationofstate_tillotson(rho,u,pressure,spsound,gamma) spsound3 = sqrt( (pressure3/rho)*( 1. + aparam + ((bparam / omega )*exp(-alpha*(neu**2)))) + & ( ((bparam*rho*u)/((omega**2)*(eta**2))) * ((1. /(u0*rho))*((2.*u) - (pressure3/rho)) & - + ((2.*alpha*neu*omega)/rho0)) + ((A/rho0)*(1. + (mu/eta**2)*((beta*2.*alpha*neu) - eta))) & + + ((2.*alpha*neu*omega)/rho0)) + ((A/rho0)*(1. + (mu_t/eta**2)*((beta*2.*alpha*neu) - eta))) & * exp(-beta*neu) )*exp(-alpha*(neu**2)) ) if (spsound3 < spsoundmin) then spsound3 = spsoundmin @@ -94,29 +94,59 @@ subroutine equationofstate_tillotson(rho,u,pressure,spsound,gamma) spsound2h = spsoundmin endif - if (rho >= rho0) then ! compressed state - pressure = pressure2 - spsound = spsound2 - else ! rho < rho0 - if (u >= u_cv) then ! vaporised expanded state - pressure = pressure3 - spsound = spsound3 - ! need to check rho ~0 and u >> u_cv (zero density becomes pressure = rho*gamma*u fermi, gamma = 2/3) - else ! u < u_cv - if (rho >= rho_iv) then - if (u > u_iv) then ! Hybrid state - pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) - spsound = sqrt(( ((u-u_iv)*(spsound3)**2) + ((u_cv-u)*(spsound2)**2) ) / (u_cv - u_iv)) - else ! u <= u_iv cold expanded + if (rho < rho0) then + if (rho >= rho_iv) then ! u <= u_iv cold expanded + if (u < u_cv) then ! hybrid state pressure = pressure2 spsound = spsound2 + elseif (u > u_iv) then ! hybrid state + pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) + spsoundh = sqrt(( ((u-u_iv)*(spsound3)**2) + ((u_cv-u)*(spsound2)**2) ) / (u_cv - u_iv)) + if (spsoundh < spsoundmin) then + spsound = spsoundmin endif - else ! rho < rho_iv Low energy state - pressure = pressure2 - (B*mu**2) - spsound = spsound2h + else ! (u >= u_cv) vaporised expanded state + pressure = pressure3 + spsound = spsound3 endif + else ! (rho < rho_iv) then (u = rho0) compressed state + pressure = pressure2 + spsound = spsound2 +endif + +! else ! (rho >= rho0) compressed state +! pressure = pressure2 +! spsound = spsound2 +! endif + +! if (rho >= rho0) then ! compressed state +! pressure = pressure2 +! spsound = spsound2 +! else ! rho < rho0 +! if (u >= u_cv) then ! vaporised expanded state +! pressure = pressure3 +! spsound = spsound3 +! ! need to check rho ~0 and u >> u_cv (zero density becomes pressure = rho*gamma*u fermi, gamma = 2/3) +! else ! u < u_cv +! ! if (rho >= rho_iv) then +! if (u > u_iv) then ! Hybrid state +! pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) +! spsound = sqrt(( ((u-u_iv)*(spsound3)**2) + ((u_cv-u)*(spsound2)**2) ) / (u_cv - u_iv)) +! else ! u <= u_iv cold expanded +! pressure = pressure2 +! spsound = spsound2 +! endif +! ! else ! rho < rho_iv Low energy state +! ! pressure = pressure2 - (B*mu_t**2) +! ! spsound = spsound2h +! ! endif +! endif +! endif +! ! print*,' INSIDE EOS: rho ',rho,' u ',u,' OUT pres ',pressure end subroutine equationofstate_tillotson @@ -128,7 +158,7 @@ end subroutine equationofstate_tillotson subroutine eos_info_tillotson(iprint) integer, intent(in) :: iprint - write(iprint,"(/,a)") ' Tillotson EoS' + write(iprint,"(/,a)") ' Tillotson EoS',pressure end subroutine eos_info_tillotson @@ -142,12 +172,12 @@ subroutine calc_uT_from_rhoP_tillotson(rho,pres,temp,u,ierr) real, intent(in) :: rho,pres real, intent(out) :: temp,u integer, intent(out) :: ierr - real :: eta, mu, neu, omega, expterm, expterm2, u1, u2, X, Y, quada, quadb, quadc + real :: eta, mu_t, neu, omega, expterm, expterm2, X, quada, quadb, quadc real :: usol(3) integer :: nsol eta = rho / rho0 - mu = eta - 1. + mu_t = eta - 1. neu = (1. / eta) - 1. ! Kegerreis et al. 2020 omega = (u / (u0*eta**2) + 1.) expterm2 = exp(-alpha*(neu**2)) @@ -155,70 +185,32 @@ subroutine calc_uT_from_rhoP_tillotson(rho,pres,temp,u,ierr) if (rho >= rho0) then ! compressed state quada = 1. - X = (pres - (A*mu) - (B*(mu**2)))/rho - quadb = -((X-(aparam*u0*(eta**2))-(bparam*u0*(eta**2))) / aparam) + X = (pres - (A*mu_t) - (B*(mu_t**2)))/rho + quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) quadc = -((u0*(eta**2)*X)/aparam) call cubicsolve(0.,quada,quadb,quadc,usol,nsol) - ! call quadratic_eq(quada,quadb,quadc,u1,u2) - ! if (u2 <= 0.) then - ! u2 = 0. - ! elseif (u2 == u1) then - ! u = u1 - ! elseif (u2 > u1) then - ! u = u1 - ! endif - - ! u = (((1./aparam)-(1./bparam))*(-A*mu - B*(mu**2))/ & - ! (rho*(1. - ((-A*mu - B*(mu**2)+pres)/(bparam*(eta**2)*u0*rho))))) - ! u = (mu*(eta**2)*u0*(aparam-bparam)*(A+(B*mu))) / & - ! (aparam*(A*mu+bparam*(eta**2)*u0*rho+B*(mu**2)-pres)) + u = maxval(usol(1:nsol)) elseif (rho >= rho_iv) then ! cold expanded quada = 1. - X = (pres - (A*mu) - (B*(mu**2)))/rho - quadb = -((X-(aparam*u0*(eta**2))-(bparam*u0*(eta**2))) / aparam) + X = (pres - (A*mu_t) - (B*(mu_t**2)))/rho + quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) quadc = -((u0*(eta**2)*X)/aparam) call cubicsolve(0.,quada,quadb,quadc,usol,nsol) - ! call quadratic_eq(quada,quadb,quadc,u1,u2) - ! if (u2 <= 0.) then - ! u2 = 0. - ! elseif (u2 == u1) then - ! u = u1 - ! elseif (u2 > u1) then - ! u = u1 - ! endif - ! u = (((1./aparam)-(1./bparam))*(-A*mu - B*(mu**2))/ & - ! (rho*(1. - ((-A*mu - B*(mu**2)+pres)/(bparam*(eta**2)*u0*rho))))) - ! u = (mu*(eta**2)*u0*(aparam-bparam)*(A+(B*mu))) / & - ! (aparam*(A*mu+bparam*(eta**2)*u0*rho+B*(mu**2)-pres)) + u = maxval(usol(1:nsol)) elseif (rho < rho_iv) then ! rho < rho_iv Low energy state quada = 1. - X = (pres - (A*mu))/rho - quadb = -((X-(aparam*u0*(eta**2))-(bparam*u0*(eta**2))) / aparam) + X = (pres - (A*mu_t))/rho + quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) quadc = -((u0*(eta**2)*X)/aparam) call cubicsolve(0.,quada,quadb,quadc,usol,nsol) - ! call quadratic_eq(quada,quadb,quadc,u1,u2) - ! if (u2 <= 0.) then - ! u2 = 0. - ! elseif (u2 == u1) then - ! u = u1 - ! elseif (u2 > u1) then - ! u = u1 - ! endif - ! u = ((-A*mu*((1./aparam)-(1./bparam))) /& - ! (rho*(1.-((pres-A*mu)/(bparam*(eta**2)*u0*rho))))) - ! u = (A*mu*(eta**2)*u0*(aparam-bparam)) / & - ! (aparam*(A*mu+bparam*(eta**2)*u0*rho-pres)) + u = maxval(usol(1:nsol)) else ! rho < rho0 vaporised expanded state quada = 1. - X = (aparam*rho)/(expterm2) - quadb = ((bparam*rho*u0*(eta**2))+((aparam*rho*u0*(eta**2))/expterm2)+(A*mu*expterm)-(pres/(expterm2)))/X - quadc = ((A*mu*expterm*u0*(eta**2)) -((pres*u0*(eta**2))/expterm2))/X + X = (pres - (A*mu_t*expterm*expterm2))/rho + quadb = ((bparam*u0*(eta**2)*expterm2)/aparam)+((u0*(eta**2)))-(X/aparam) + quadc = -((X*u0*(eta**2))/aparam) call cubicsolve(0.,quada,quadb,quadc,usol,nsol) - !call quadratic_eq(quada,quadb,quadc,u1,u2) u = maxval(usol(1:nsol)) - ! u = ((eta**2)*u0*(bparam-aparam)*(pres-A*mu*expterm*expterm2) ) / & - ! (aparam*rho*(A*mu*expterm + bparam*(eta**2)*u0 - pres)) - ! need to check rho ~0 and u >> u_cv (zero density becomes u = pres / (rho*gamma) fermi, gamma = 2/3) endif temp = 300. @@ -227,40 +219,10 @@ subroutine calc_uT_from_rhoP_tillotson(rho,pres,temp,u,ierr) u = 0. ierr = 1 endif +! print*,' INSIDE EOS: rho ',rho,' pres ',pres,' OUT u ',u end subroutine calc_uT_from_rhoP_tillotson -!----------------------------------------------------------------------- -!+ -! calc quadratic for internal energy -!+ -!----------------------------------------------------------------------- -! function quadratic_eq(a,b,c,x1,x2) -! real, intent(in) :: a,b,c -! real, intent(out) :: x1,x2 -! real :: dis, sqrt_val - -! dis = (b**2) - (4. * a * c) ! calc discriminant using formula -! sqrt_val = sqrt(abs(dis)) - -! if (dis > 0) then ! check discriminant -! ! real and different roots -! x1 = (-b + sqrt_val)/(2 * a) -! x2 = (-b - sqrt_val)/(2 * a) - -! elseif (dis == 0) then -! ! real and same roots -! x1 = -b / (2. * a) -! x2 = x1 - -! else ! if discriminant < 0 - ! Complex Roots - ! print(- b / (2 * a), + i, sqrt_val / (2 * a)) - ! print(- b / (2 * a), - i, sqrt_val / (2 * a)) -! endif - -! end function quadratic_eq - !----------------------------------------------------------------------- !+ ! reads equation of state options from the input file From 47d566360ab567b71e08c6ed6f1f1456d8fc8211 Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Tue, 24 Sep 2024 14:15:18 +1000 Subject: [PATCH 24/29] (tillotson) add calc u from prho to testing suite --- src/tests/test_eos.f90 | 63 ++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 60 insertions(+), 3 deletions(-) diff --git a/src/tests/test_eos.f90 b/src/tests/test_eos.f90 index 44f09afe8..fbeacadcd 100644 --- a/src/tests/test_eos.f90 +++ b/src/tests/test_eos.f90 @@ -44,13 +44,15 @@ subroutine test_eos(ntests,npass) call set_units(mass=solarm,dist=1.d16,G=1.d0) call test_init(ntests, npass) + + call test_u_from_Prho(ntests,npass) call test_barotropic(ntests, npass) ! call test_helmholtz(ntests, npass) call test_idealplusrad(ntests, npass) - do irecomb = 0,3 - call test_hormone(ntests,npass) - enddo +do irecomb = 0,3 + call test_hormone(ntests,npass) +enddo call test_eos_stratified(ntests,npass) @@ -104,7 +106,62 @@ subroutine test_init(ntests, npass) end subroutine test_init +!---------------------------------------------------------------------------- +!+ +! test that the routine to solve for u from pressure and density works +!+ +!---------------------------------------------------------------------------- +subroutine test_u_from_Prho(ntests, npass) + use io, only:id,master,stdout + use eos, only:init_eos,equationofstate,calc_temp_and_ene,gamma + use testutils, only:checkval,checkvalbuf_start,checkvalbuf,checkvalbuf_end,update_test_scores + use units, only:unit_density,unit_pressure,unit_ergg + use physcon, only:Rg,kb_on_mh + use table_utils, only:logspace + integer, intent(inout) :: ntests,npass + integer :: npts,ieos,ierr,i,j,nfail(1),ncheck(1) + real :: rhocodei,presi,pres2,dum,csound,eni,temp,ponrhoi,mu,tol,errmax(1),code_eni,en_back + real, allocatable :: rhogrid(:),ugrid(:),Pgrid(:) + + if (id==master) write(*,"(/,a)") '--> testing retrieval of u from P and rho' + + ieos = 23 + gamma = 5./3. + npts = 30 +! call get_rhoT_grid(npts,rhogrid,ugrid) + allocate(rhogrid(npts),ugrid(npts)) + call logspace(rhogrid,1e01,1e05) ! cgs + call logspace(ugrid,1e09,1e15) ! cgs +! print*,'rhogrid ',rhogrid,'ugrid ',ugrid + + dum = 0. + tol = 1.e-09 + nfail = 0; ncheck = 0; errmax = 0. + call init_eos(ieos,ierr) + do i=1,npts + do j=1,npts + + ! get u from P, rho + rhocodei = rhogrid(i)/unit_density + code_eni = ugrid(i)/unit_ergg + ! print*,' rhogridi ',rhogrid(i),' ugridi ',ugrid(i) + ! print*,' rhocodei ',rhocodei,' code_eni ',code_eni + call equationofstate(ieos,ponrhoi,csound,rhocodei,dum,dum,dum,temp,code_eni) + presi = ponrhoi * rhocodei + ! print*,' ponrhoi ',ponrhoi + ! print*,' sending en ',code_eni,'sending rho ',rhocodei,' got pres ',presi + + call calc_temp_and_ene(ieos,rhocodei,presi,en_back,temp,ierr) ! out = energy and temp + ! print*,' sending rho ',rhocodei,'sending pres ',presi,' got en ',en_back + ! read* + + call checkvalbuf(code_eni,en_back,tol,'Check recovery of u from rho, P',nfail(1),ncheck(1),errmax(1),use_rel_tol) + enddo + enddo + call checkvalbuf_end('Check recovery of u from rho, P',ncheck(1),nfail(1),errmax(1),tol) + call update_test_scores(ntests,nfail,npass) +end subroutine test_u_from_Prho !---------------------------------------------------------------------------- !+ ! test ideal gas plus radiation eos: Check that P, T calculated from rho, u gives back From e3091aa1f713035dd8c0aaa57b9df4830e514d57 Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Tue, 24 Sep 2024 14:43:34 +1000 Subject: [PATCH 25/29] (tillotson) removed unnecessary code --- src/main/eos_tillotson.f90 | 30 ------------------------------ 1 file changed, 30 deletions(-) diff --git a/src/main/eos_tillotson.f90 b/src/main/eos_tillotson.f90 index a275ce694..8eedd28a3 100644 --- a/src/main/eos_tillotson.f90 +++ b/src/main/eos_tillotson.f90 @@ -118,36 +118,6 @@ subroutine equationofstate_tillotson(rho,u,pressure,spsound,gamma) spsound = spsound2 endif -! else ! (rho >= rho0) compressed state -! pressure = pressure2 -! spsound = spsound2 -! endif - -! if (rho >= rho0) then ! compressed state -! pressure = pressure2 -! spsound = spsound2 -! else ! rho < rho0 -! if (u >= u_cv) then ! vaporised expanded state -! pressure = pressure3 -! spsound = spsound3 -! ! need to check rho ~0 and u >> u_cv (zero density becomes pressure = rho*gamma*u fermi, gamma = 2/3) -! else ! u < u_cv -! ! if (rho >= rho_iv) then -! if (u > u_iv) then ! Hybrid state -! pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) -! spsound = sqrt(( ((u-u_iv)*(spsound3)**2) + ((u_cv-u)*(spsound2)**2) ) / (u_cv - u_iv)) -! else ! u <= u_iv cold expanded -! pressure = pressure2 -! spsound = spsound2 -! endif -! ! else ! rho < rho_iv Low energy state -! ! pressure = pressure2 - (B*mu_t**2) -! ! spsound = spsound2h -! ! endif -! endif -! endif -! ! print*,' INSIDE EOS: rho ',rho,' u ',u,' OUT pres ',pressure - end subroutine equationofstate_tillotson !---------------------------------------------------------------- From 6e4764fb9fc4fd761b9396c8d7b8aa5d8c7d62e4 Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Tue, 1 Oct 2024 14:48:33 +1000 Subject: [PATCH 26/29] (Tillotson) updated eos piecewise and testing of eos --- src/main/eos_tillotson.f90 | 167 +++++++++++++++++++++++++++---------- src/tests/test_eos.f90 | 30 +++---- 2 files changed, 139 insertions(+), 58 deletions(-) diff --git a/src/main/eos_tillotson.f90 b/src/main/eos_tillotson.f90 index 8eedd28a3..a38815d7d 100644 --- a/src/main/eos_tillotson.f90 +++ b/src/main/eos_tillotson.f90 @@ -94,29 +94,62 @@ subroutine equationofstate_tillotson(rho,u,pressure,spsound,gamma) spsound2h = spsoundmin endif - if (rho < rho0) then - if (rho >= rho_iv) then ! u <= u_iv cold expanded - if (u < u_cv) then ! hybrid state - pressure = pressure2 - spsound = spsound2 - elseif (u > u_iv) then ! hybrid state - pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) - spsoundh = sqrt(( ((u-u_iv)*(spsound3)**2) + ((u_cv-u)*(spsound2)**2) ) / (u_cv - u_iv)) - if (spsoundh < spsoundmin) then - spsound = spsoundmin - endif - else ! (u >= u_cv) vaporised expanded state + + if (rho < rho_iv) then ! low energy expanded + pressure = pressure2 - (B*mu_t**2) + spsound = spsound2h + print*,'in low energy state, p = ',pressure + else ! (rho >= rho_iv) + if (rho >= rho0) then ! compressed + pressure = pressure2 + spsound = spsound2 + print*,'in compressed state, p = ',pressure + else ! (rho_iv <= rho < rho0) + if (u >= u_cv) then ! hot / vapourised expanded pressure = pressure3 spsound = spsound3 + print*,'in hot expanded state, p = ',pressure + else ! (u < u_cv) + ! if (u <= u_iv) then ! cold expanded + pressure = pressure2 + spsound = spsound2 + print*,'in cold expanded state, p = ',pressure + ! else ! (u > u_iv) hybrid + ! pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) + ! spsoundh = sqrt(( ((u-u_iv)*(spsound3)**2) + ((u_cv-u)*(spsound2)**2) ) / (u_cv - u_iv)) + ! if (spsoundh < spsoundmin) then + ! spsound = spsoundmin + ! endif + ! print*,'in hybrid state, p = ',pressure + ! endif endif - else ! (rho < rho_iv) then (u = rho0) compressed state - pressure = pressure2 - spsound = spsound2 -endif + endif + +! if (rho < rho0) then +! if (rho >= rho_iv) then ! u <= u_iv cold expanded +! if (u < u_cv) then ! hybrid state +! pressure = pressure2 +! spsound = spsound2 +! elseif (u > u_iv) then ! hybrid state +! pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) +! spsoundh = sqrt(( ((u-u_iv)*(spsound3)**2) + ((u_cv-u)*(spsound2)**2) ) / (u_cv - u_iv)) +! if (spsoundh < spsoundmin) then +! spsound = spsoundmin +! endif +! print*,'in hybrid state, p = ',pressure +! else ! (u >= u_cv) vaporised expanded state +! pressure = pressure3 +! spsound = spsound3 +! endif +! else ! (rho < rho_iv) then (u = rho0) compressed state +! pressure = pressure2 +! spsound = spsound2 +! endif end subroutine equationofstate_tillotson @@ -142,7 +175,7 @@ subroutine calc_uT_from_rhoP_tillotson(rho,pres,temp,u,ierr) real, intent(in) :: rho,pres real, intent(out) :: temp,u integer, intent(out) :: ierr - real :: eta, mu_t, neu, omega, expterm, expterm2, X, quada, quadb, quadc + real :: eta, mu_t, neu, omega, expterm, expterm2, X, quada, quadb, quadc, pres_iv real :: usol(3) integer :: nsol @@ -153,36 +186,84 @@ subroutine calc_uT_from_rhoP_tillotson(rho,pres,temp,u,ierr) expterm2 = exp(-alpha*(neu**2)) expterm = exp(-beta*neu) - if (rho >= rho0) then ! compressed state - quada = 1. - X = (pres - (A*mu_t) - (B*(mu_t**2)))/rho - quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) - quadc = -((u0*(eta**2)*X)/aparam) - call cubicsolve(0.,quada,quadb,quadc,usol,nsol) - u = maxval(usol(1:nsol)) - elseif (rho >= rho_iv) then ! cold expanded - quada = 1. - X = (pres - (A*mu_t) - (B*(mu_t**2)))/rho - quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) - quadc = -((u0*(eta**2)*X)/aparam) - call cubicsolve(0.,quada,quadb,quadc,usol,nsol) - u = maxval(usol(1:nsol)) - elseif (rho < rho_iv) then ! rho < rho_iv Low energy state - quada = 1. + pres_iv = (( aparam + ( bparam / ((u_iv / (u0*(rho_iv/rho0)**2) + 1.)) ) ) * rho_iv*u_iv) & + + A*((rho_iv/rho0)-1.) + B*((rho_iv/rho0)-1.)**2 + + quada = 1. + + if (rho < rho_iv) then ! low energy expanded + ! quada = 1. X = (pres - (A*mu_t))/rho quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) quadc = -((u0*(eta**2)*X)/aparam) call cubicsolve(0.,quada,quadb,quadc,usol,nsol) u = maxval(usol(1:nsol)) - else ! rho < rho0 vaporised expanded state - quada = 1. - X = (pres - (A*mu_t*expterm*expterm2))/rho - quadb = ((bparam*u0*(eta**2)*expterm2)/aparam)+((u0*(eta**2)))-(X/aparam) - quadc = -((X*u0*(eta**2))/aparam) - call cubicsolve(0.,quada,quadb,quadc,usol,nsol) - u = maxval(usol(1:nsol)) + print*,'in low energy state, u = ',u + else ! (rho >= rho_iv) + if (rho >= rho0) then ! compressed + ! quada = 1. + X = (pres - (A*mu_t) - (B*(mu_t**2)))/rho + quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) + quadc = -((u0*(eta**2)*X)/aparam) + call cubicsolve(0.,quada,quadb,quadc,usol,nsol) + u = maxval(usol(1:nsol)) + print*,'in compressed state, u = ',u + else ! (rho_iv <= rho < rho0) + if (pres <= pres_iv) then ! cold expanded + ! quada = 1. + X = (pres - (A*mu_t) - (B*(mu_t**2)))/rho + quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) + quadc = -((u0*(eta**2)*X)/aparam) + call cubicsolve(0.,quada,quadb,quadc,usol,nsol) + u = maxval(usol(1:nsol)) + print*,'in cold expanded state, u = ',u + else ! (pres > pres_iv) hot / vapourised expanded + ! quada = 1. + X = (pres - (A*mu_t*expterm*expterm2))/rho + quadb = ((bparam*u0*(eta**2)*expterm2)/aparam)+((u0*(eta**2)))-(X/aparam) + quadc = -((X*u0*(eta**2))/aparam) + call cubicsolve(0.,quada,quadb,quadc,usol,nsol) + u = maxval(usol(1:nsol)) + print*,'in hot expanded state, u = ',u + ! else ! (u > u_iv) hybrid + ! pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) + ! spsoundh = sqrt(( ((u-u_iv)*(spsound3)**2) + ((u_cv-u)*(spsound2)**2) ) / (u_cv - u_iv)) + ! print*,'in hybrid state, u = ',u + ! endif + endif + endif endif +! if (rho >= rho0) then ! compressed state +! quada = 1. +! X = (pres - (A*mu_t) - (B*(mu_t**2)))/rho +! quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) +! quadc = -((u0*(eta**2)*X)/aparam) +! call cubicsolve(0.,quada,quadb,quadc,usol,nsol) +! u = maxval(usol(1:nsol)) +! elseif (rho >= rho_iv) then ! cold expanded +! quada = 1. +! X = (pres - (A*mu_t) - (B*(mu_t**2)))/rho +! quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) +! quadc = -((u0*(eta**2)*X)/aparam) +! call cubicsolve(0.,quada,quadb,quadc,usol,nsol) +! u = maxval(usol(1:nsol)) +! elseif (rho < rho_iv) then ! rho < rho_iv Low energy state +! quada = 1. +! X = (pres - (A*mu_t))/rho +! quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) +! quadc = -((u0*(eta**2)*X)/aparam) +! call cubicsolve(0.,quada,quadb,quadc,usol,nsol) +! u = maxval(usol(1:nsol)) +! else ! rho < rho0 vaporised expanded state +! quada = 1. +! X = (pres - (A*mu_t*expterm*expterm2))/rho +! quadb = ((bparam*u0*(eta**2)*expterm2)/aparam)+((u0*(eta**2)))-(X/aparam) +! quadc = -((X*u0*(eta**2))/aparam) +! call cubicsolve(0.,quada,quadb,quadc,usol,nsol) +! u = maxval(usol(1:nsol)) +! endif + temp = 300. ierr = 0 if (u <= 0.) then diff --git a/src/tests/test_eos.f90 b/src/tests/test_eos.f90 index fbeacadcd..4b8c3f92b 100644 --- a/src/tests/test_eos.f90 +++ b/src/tests/test_eos.f90 @@ -46,15 +46,15 @@ subroutine test_eos(ntests,npass) call test_init(ntests, npass) call test_u_from_Prho(ntests,npass) - call test_barotropic(ntests, npass) -! call test_helmholtz(ntests, npass) - call test_idealplusrad(ntests, npass) +! call test_barotropic(ntests, npass) +! ! call test_helmholtz(ntests, npass) +! call test_idealplusrad(ntests, npass) -do irecomb = 0,3 - call test_hormone(ntests,npass) -enddo +! do irecomb = 0,3 +! call test_hormone(ntests,npass) +! enddo - call test_eos_stratified(ntests,npass) +! call test_eos_stratified(ntests,npass) if (id==master) write(*,"(/,a)") '<-- EQUATION OF STATE TEST COMPLETE' @@ -130,8 +130,8 @@ subroutine test_u_from_Prho(ntests, npass) npts = 30 ! call get_rhoT_grid(npts,rhogrid,ugrid) allocate(rhogrid(npts),ugrid(npts)) - call logspace(rhogrid,1e01,1e05) ! cgs - call logspace(ugrid,1e09,1e15) ! cgs + call logspace(rhogrid,1e-5,1e15) ! cgs + call logspace(ugrid,1e-5,1e25) ! cgs ! print*,'rhogrid ',rhogrid,'ugrid ',ugrid dum = 0. @@ -144,16 +144,16 @@ subroutine test_u_from_Prho(ntests, npass) ! get u from P, rho rhocodei = rhogrid(i)/unit_density code_eni = ugrid(i)/unit_ergg - ! print*,' rhogridi ',rhogrid(i),' ugridi ',ugrid(i) - ! print*,' rhocodei ',rhocodei,' code_eni ',code_eni + print*,' rhogridi ',rhogrid(i),' ugridi ',ugrid(i) + print*,' rhocodei ',rhocodei,' code_eni ',code_eni call equationofstate(ieos,ponrhoi,csound,rhocodei,dum,dum,dum,temp,code_eni) presi = ponrhoi * rhocodei - ! print*,' ponrhoi ',ponrhoi - ! print*,' sending en ',code_eni,'sending rho ',rhocodei,' got pres ',presi + print*,' ponrhoi ',ponrhoi + print*,' sending en ',code_eni,'sending rho ',rhocodei,' got pres ',presi call calc_temp_and_ene(ieos,rhocodei,presi,en_back,temp,ierr) ! out = energy and temp - ! print*,' sending rho ',rhocodei,'sending pres ',presi,' got en ',en_back - ! read* + print*,' sending rho ',rhocodei,'sending pres ',presi,' got en ',en_back + read* call checkvalbuf(code_eni,en_back,tol,'Check recovery of u from rho, P',nfail(1),ncheck(1),errmax(1),use_rel_tol) enddo From 5034e52e5cc2bc4b3ad88ad2425468458a21bcc0 Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Fri, 4 Oct 2024 17:17:57 +1000 Subject: [PATCH 27/29] (tillotson) fixed piecewise eos to within 1.E-12 tol --- src/main/eos_tillotson.f90 | 240 +++++++++++++++++++++++-------------- 1 file changed, 150 insertions(+), 90 deletions(-) diff --git a/src/main/eos_tillotson.f90 b/src/main/eos_tillotson.f90 index a38815d7d..08fc61b5d 100644 --- a/src/main/eos_tillotson.f90 +++ b/src/main/eos_tillotson.f90 @@ -28,7 +28,7 @@ module eos_tillotson real :: u_iv = 4.72e10 ! erg/g real :: u_cv = 1.82e11 ! erg/g real :: rho_iv = 2.57 ! g/cm^3 incipient vaporisation density of gabbroic anorthosite lpp from Ahrens and Okeefe (1977) table 1 -! real :: rho_cv = 8.47 ! g/cm^3 complete vaporisation density of gabbroic anorthosite lpp from Ahrens and Okeefe (1977) table 2 + real :: rho_cv = 8.47 ! g/cm^3 complete vaporisation density of gabbroic anorthosite lpp from Ahrens and Okeefe (1977) table 2 private :: rho0, aparam, bparam, A, B, alpha, beta, u0, rho_iv, u_iv, u_cv ! private :: spsoundmin, spsound2, spsound3, pressure2, pressure3 ! private :: expterm, expterm2, sqrtterm, sqrtterm1, sqrtterm2 @@ -47,84 +47,112 @@ subroutine init_eos_tillotson(ierr) end subroutine init_eos_tillotson subroutine equationofstate_tillotson(rho,u,pressure,spsound,gamma) - real, intent(in) :: rho,u + use units, only: unit_density, unit_pressure, unit_ergg + real, intent(inout) :: rho,u real, intent(out) :: pressure, spsound, gamma real :: eta, mu_t, neu, omega, spsoundmin, spsound2, spsoundh, spsound2h, spsound3, pressure2, pressure3 + +! rho = rho*unit_density +! u = u*unit_ergg - eta = rho / rho0 + eta = rho/ rho0 mu_t = eta - 1. neu = (1. / eta) - 1. ! Kegerreis et al. 2020 - omega = (u / (u0*eta**2)) + 1. - spsoundmin = sqrt(A / rho) ! wave speed estimate + omega = (u / (u0*(eta**2))) + 1. + spsoundmin = sqrt(A / rho) ! wave speed estimate Benz and Asphaug ! gamma = 2./3. ! Brundage 2013 eq (2) cold expanded and compressed states - pressure2 = (( aparam + ( bparam / omega ) ) * rho*u) + A*mu_t + B*mu_t**2 + pressure2 = (( aparam + ( bparam / omega ) ) * rho*u) + (A*mu_t) + (B*(mu_t**2)) if (pressure2 <= 0.) then pressure2 = 0. endif ! Brundage 2013 eq (3) completely vaporised pressure3 = (aparam*rho*u) + ( ((bparam*rho*u)/omega) + & - A*mu_t*exp(-beta*neu)) * exp(-alpha*(neu**2)) + A*mu_t*exp((-1.*beta)*neu)) * exp((-1.*alpha)*(neu**2)) if (pressure3 <= 0.) then pressure3 = 0. endif ! sound speed from Kegereis et al. 2020 spsound2 = sqrt( ((pressure2 / rho) * (1. + aparam + (bparam / omega ))) + & - (((bparam*(omega - 1.)) / (omega**2)) * (2.*u - pressure2/rho)) + & - ((1./rho)*(A + B*((eta**2) - 1.)))) + (((bparam*(omega - 1.)) / (omega**2)) * ((2.*u) - (pressure2/rho))) + & + ((1./rho)*(A + (B*((eta**2) - 1.))))) if (spsound2 < spsoundmin) then spsound2 = spsoundmin endif - spsound3 = sqrt( (pressure3/rho)*( 1. + aparam + ((bparam / omega )*exp(-alpha*(neu**2)))) + & - ( ((bparam*rho*u)/((omega**2)*(eta**2))) * ((1. /(u0*rho))*((2.*u) - (pressure3/rho)) & - + ((2.*alpha*neu*omega)/rho0)) + ((A/rho0)*(1. + (mu_t/eta**2)*((beta*2.*alpha*neu) - eta))) & - * exp(-beta*neu) )*exp(-alpha*(neu**2)) ) + spsound3 = sqrt( ((pressure3/rho)*( 1. + aparam + ((bparam / omega )*exp((-1.*alpha)*(neu**2))))) + & + ( (((bparam*rho*u)/((omega**2)*(eta**2))) * ((1. /(u0*rho)) * ((2.*u) - (pressure3/rho)) & + + ((2.*alpha*neu*omega)/rho0))) + ((A/rho0)*(1. + ((mu_t/(eta**2))*((beta+(2.*alpha*neu) - eta))) & + * exp((-1.*beta)*neu) )*exp((-1.*alpha)*(neu**2)) ))) if (spsound3 < spsoundmin) then spsound3 = spsoundmin endif spsound2h = sqrt( ((pressure2 / rho) * (1. + aparam + (bparam / omega ))) + & - (((bparam*(omega - 1.)) / (omega**2)) * (2.*u - pressure2/rho)) + & - ((1./rho)*(A + B*((eta**2) - 1.))) - ((((2.*B)/rho0**2)*rho) - (2.*B/rho0))) + (((bparam*(omega - 1.)) / (omega**2)) * ((2.*u) - (pressure2/rho))) + & + ((1./rho)*(A + (B*((eta**2) - 1.)))) - (((B*(mu_t**2))/rho)*(1+aparam+(bparam/(omega**2))))) +! spsound2h = sqrt( ((pressure2 / rho) * (1. + aparam + (bparam / omega ))) + & +! (((bparam*(omega - 1.)) / (omega**2)) * ((2.*u) - pressure2/rho)) + & +! ((1./rho)*(A + B*((eta**2) - 1.))) - ((((2.*B)/(rho0**2))*rho) - (2.*B/rho0))) if (spsound2h < spsoundmin) then spsound2h = spsoundmin endif - - - if (rho < rho_iv) then ! low energy expanded - pressure = pressure2 - (B*mu_t**2) - spsound = spsound2h - print*,'in low energy state, p = ',pressure - else ! (rho >= rho_iv) - if (rho >= rho0) then ! compressed - pressure = pressure2 - spsound = spsound2 - print*,'in compressed state, p = ',pressure - else ! (rho_iv <= rho < rho0) - if (u >= u_cv) then ! hot / vapourised expanded - pressure = pressure3 - spsound = spsound3 - print*,'in hot expanded state, p = ',pressure - else ! (u < u_cv) - ! if (u <= u_iv) then ! cold expanded - pressure = pressure2 - spsound = spsound2 - print*,'in cold expanded state, p = ',pressure - ! else ! (u > u_iv) hybrid - ! pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) - ! spsoundh = sqrt(( ((u-u_iv)*(spsound3)**2) + ((u_cv-u)*(spsound2)**2) ) / (u_cv - u_iv)) - ! if (spsoundh < spsoundmin) then - ! spsound = spsoundmin - ! endif - ! print*,'in hybrid state, p = ',pressure - ! endif - endif - endif + +! print*,' in eos - calc pres from rho and u' +! print*,' rho in = ',rho,' u in = ',u +! print*,'in low energy state, p = ',pressure,' pres in code units = ',pressure/unit_pressure ! rho < rho_iv ; u < u_cv = ISSUES +! print*,'in hybrid state, p = ',pressure,' pres in code units = ',pressure/unit_pressure ! ISSUES +! print*,'in compressed state, p = ',pressure,' pres in code units = ',pressure/unit_pressure ! rho > rho0 ; u >= u0 +! print*,'in hot expanded state, p = ',pressure,' pres in code units = ',pressure/unit_pressure ! rho_iv <= rho < rho0 ; u >= u_cv + + +! if (rho < rho_iv) then ! low energy expanded +! pressure = pressure2 - (B*(mu_t**2)) +! spsound = spsound2h +! print*,' in eos - calc pres from rho and u' +! print*,' rho in = ',rho,' u in = ',u +! print*,'in low energy state, p = ',pressure,' pres in code units = ',pressure/unit_pressure +! else ! (rho >= rho_iv) + if (rho >= rho0) then ! compressed and cold expanded + pressure = pressure2 + spsound = spsound2 + ! print*,' in eos - calc pres from rho and u' + ! print*,' rho in = ',rho,' u in = ',u + ! print*,' spsound = ',spsound + ! print*,'in compressed state, p = ',pressure,' pres in code units = ',pressure/unit_pressure + else ! (rho_iv <= rho < rho0) +! if (u >= u_cv) then ! hot / vapourised expanded + pressure = pressure3 + spsound = spsound3 + ! print*,' in eos - calc pres from rho and u' + ! print*,' rho in = ',rho,' u in = ',u + ! print*,' spsound = ',spsound + ! print*,'in hot expanded state, p = ',pressure,' pres in code units = ',pressure/unit_pressure endif +! print*,'leaving eos (pres calc) ' +! else ! (u < u_cv) +! if (u <= u_iv) then ! cold expanded +! pressure = pressure2 +! spsound = spsound2 +! print*,' in eos - calc pres from rho and u' +! print*,' rho in = ',rho,' u in = ',u +! print*,'in cold expanded state, p = ',pressure,' pres in code units = ',pressure/unit_pressure +! else ! (u > u_iv) hybrid +! pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) +! spsoundh = sqrt(( ((u-u_iv)*(spsound3)**2) + ((u_cv-u)*(spsound2)**2) ) / (u_cv - u_iv)) +! if (spsoundh < spsoundmin) then +! spsound = spsoundmin +! endif +! print*,' in eos - calc pres from rho and u' +! print*,' rho in = ',rho,' u in = ',u +! print*,'in hybrid state, p = ',pressure,' pres in code units = ',pressure/unit_pressure +! endif +! endif +! endif +! endif ! if (rho < rho0) then ! if (rho >= rho_iv) then ! u <= u_iv cold expanded @@ -172,6 +200,7 @@ end subroutine eos_info_tillotson !----------------------------------------------------------------------- subroutine calc_uT_from_rhoP_tillotson(rho,pres,temp,u,ierr) use cubic, only:cubicsolve + use units, only:unit_density,unit_pressure,unit_ergg real, intent(in) :: rho,pres real, intent(out) :: temp,u integer, intent(out) :: ierr @@ -179,61 +208,92 @@ subroutine calc_uT_from_rhoP_tillotson(rho,pres,temp,u,ierr) real :: usol(3) integer :: nsol +! rho = rho*unit_density +! pres = pres*unit_pressure + eta = rho / rho0 mu_t = eta - 1. neu = (1. / eta) - 1. ! Kegerreis et al. 2020 - omega = (u / (u0*eta**2) + 1.) + omega = (u / (u0*(eta**2)) + 1.) expterm2 = exp(-alpha*(neu**2)) expterm = exp(-beta*neu) - pres_iv = (( aparam + ( bparam / ((u_iv / (u0*(rho_iv/rho0)**2) + 1.)) ) ) * rho_iv*u_iv) & - + A*((rho_iv/rho0)-1.) + B*((rho_iv/rho0)-1.)**2 +! pres_iv = (( aparam + ( bparam / ((u_iv / (u0*(rho_iv/rho0)**2) + 1.)) ) ) * rho_iv*u_iv) & + ! + A*((rho_iv/rho0)-1.) + B*((rho_iv/rho0)-1.)**2 +! pres_cv = (aparam*rho_cv*u_cv) + ( ((bparam*rho_cv*u_cv)/((u_cv / (u0*(rho_cv/rho0)**2) + 1.))) & + ! + A*((rho_cv/rho0) - 1.)*exp((-1.*beta)*((rho0/rho_cv) - 1.))) * exp((-1.*alpha)*(((rho0/rho_cv) - 1.)**2)) - quada = 1. + +! quada = 1. +! X = (pres - (A*mu_t))/rho +! quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) +! quadc = -((u0*(eta**2)*X)/aparam) +! print*,' X = ',X,' a = ',quada,' b = ',quadb,' c = ',quadc +! call cubicsolve(0.,quada,quadb,quadc,usol,nsol) +! print*,' Solutions : ',usol,' # of sols : ',nsol +! u = maxval(usol(1:nsol)) +! print*,' in eos - calc u from rho and pres' +! print*,' rho in = ',rho,' pres in = ',pres +! print*,'in low energy state, u = ',u,' u in code units = ',u/unit_ergg ! rho < rho_iv ; pre < pres_cv = ISSUES +! print*,'in hot expanded state, u = ',u,' u in code units = ',u/unit_ergg ! rho_iv <= rho < rho0 ; pres >= pres_cv +! print*,'in compressed state, u = ',u,' u in code units = ',u/unit_ergg ! rho > rho0 ; pres >= pres0 +! print*,'leaving eos (u calc) ' - if (rho < rho_iv) then ! low energy expanded - ! quada = 1. - X = (pres - (A*mu_t))/rho - quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) - quadc = -((u0*(eta**2)*X)/aparam) +! if (rho < rho_iv) then ! low energy expanded +! quada = 1. +! X = (pres - (A*mu_t))/rho +! quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) +! quadc = -((u0*(eta**2)*X)/aparam) +! call cubicsolve(0.,quada,quadb,quadc,usol,nsol) +! u = maxval(usol(1:nsol)) +! print*,' in eos - calc u from rho and pres' +! print*,' rho in = ',rho,' pres in = ',pres +! print*,'in low energy state, u = ',u,' u in code units = ',u/unit_ergg +! else ! (rho >= rho_iv) + if (rho >= rho0) then ! compressed and cold expanded + quada = 1. + X = pres - (A*mu_t) - (B*(mu_t**2)) + quadb = (((u0*(eta**2)*rho)*(aparam + bparam)) - X)/(aparam*rho) + quadc = -1.*((X*u0*(eta**2))/(rho*aparam)) + call cubicsolve(0.,quada,quadb,quadc,usol,nsol) + ! u = maxval(usol(1:nsol)) + ! print*,' in eos - calc u from rho and pres' + ! print*,' rho in = ',rho,' pres in = ',pres + ! print*,'in compressed state, u = ',u,' u in code units = ',u/unit_ergg + else ! (rho_iv <= rho < rho0) +! if (pres >= pres_cv) then ! hot / vapourised expanded + quada = 1. + X = (pres - (A*mu_t*expterm*expterm2)) + quadb = ((aparam*rho*u0*(eta**2)) + (bparam*rho*u0*(eta**2)*expterm2) -X)/(rho*aparam) ! (u0*(eta**2))+((bparam*u0*(eta**2)*expterm2)/aparam) - (X/aparam) + quadc = -1.*(X*u0*(eta**2))/(aparam*rho) ! -1.*((X*u0*(eta**2))/aparam) call cubicsolve(0.,quada,quadb,quadc,usol,nsol) u = maxval(usol(1:nsol)) - print*,'in low energy state, u = ',u - else ! (rho >= rho_iv) - if (rho >= rho0) then ! compressed - ! quada = 1. - X = (pres - (A*mu_t) - (B*(mu_t**2)))/rho - quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) - quadc = -((u0*(eta**2)*X)/aparam) - call cubicsolve(0.,quada,quadb,quadc,usol,nsol) - u = maxval(usol(1:nsol)) - print*,'in compressed state, u = ',u - else ! (rho_iv <= rho < rho0) - if (pres <= pres_iv) then ! cold expanded - ! quada = 1. - X = (pres - (A*mu_t) - (B*(mu_t**2)))/rho - quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) - quadc = -((u0*(eta**2)*X)/aparam) - call cubicsolve(0.,quada,quadb,quadc,usol,nsol) - u = maxval(usol(1:nsol)) - print*,'in cold expanded state, u = ',u - else ! (pres > pres_iv) hot / vapourised expanded - ! quada = 1. - X = (pres - (A*mu_t*expterm*expterm2))/rho - quadb = ((bparam*u0*(eta**2)*expterm2)/aparam)+((u0*(eta**2)))-(X/aparam) - quadc = -((X*u0*(eta**2))/aparam) - call cubicsolve(0.,quada,quadb,quadc,usol,nsol) - u = maxval(usol(1:nsol)) - print*,'in hot expanded state, u = ',u - ! else ! (u > u_iv) hybrid - ! pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) - ! spsoundh = sqrt(( ((u-u_iv)*(spsound3)**2) + ((u_cv-u)*(spsound2)**2) ) / (u_cv - u_iv)) - ! print*,'in hybrid state, u = ',u - ! endif - endif - endif + ! print*,' in eos - calc u from rho and pres' + ! print*,' rho in = ',rho,' pres in = ',pres + ! print*,'in hot expanded state, u = ',u,' u in code units = ',u/unit_ergg endif - +! print*,'leaving eos (pres calc) ' +! else ! (pres < pres_cv) +! if (pres <= pres_iv) then ! cold expanded +! quada = 1. +! X = (pres - (A*mu_t) - (B*(mu_t**2)))/rho +! quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) +! quadc = -((u0*(eta**2)*X)/aparam) +! call cubicsolve(0.,quada,quadb,quadc,usol,nsol) +! u = maxval(usol(1:nsol)) +! print*,' in eos - calc u from rho and pres' +! print*,' rho in = ',rho,' pres in = ',pres +! print*,'in cold expanded state, u = ',u,' u in code units = ',u/unit_ergg +! else ! (u > u_iv) hybrid +! ! pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) +! ! spsoundh = sqrt(( ((u-u_iv)*(spsound3)**2) + ((u_cv-u)*(spsound2)**2) ) / (u_cv - u_iv)) +! ! print*,' in eos - calc u from rho and pres' +! ! print*,' rho in = ',rho,' pres in = ',pres +! ! print*,'in hybrid state, u = ',u,' u in code units = ',u/unit_ergg +! endif +! endif +! endif +! print*,'leaving eos ' ! if (rho >= rho0) then ! compressed state ! quada = 1. ! X = (pres - (A*mu_t) - (B*(mu_t**2)))/rho From f59bfd19cbe1f6060234c26f6048d50f00d7fe1c Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Wed, 9 Oct 2024 16:44:29 +1100 Subject: [PATCH 28/29] (tillotson) removed low energy cold expansion state as negative pressures ensue --- src/main/eos_tillotson.f90 | 157 ++++++++----------------------------- 1 file changed, 33 insertions(+), 124 deletions(-) diff --git a/src/main/eos_tillotson.f90 b/src/main/eos_tillotson.f90 index 08fc61b5d..83137308c 100644 --- a/src/main/eos_tillotson.f90 +++ b/src/main/eos_tillotson.f90 @@ -21,7 +21,7 @@ module eos_tillotson implicit none real :: rho0 = 2.7 ! g/cm^3 zero-pressure density (Basalt) from Benz & Asphaug 1999 ! aparam, bparam, A, B, energy0 material-dependent Tillotson parameters (Basalt) - real :: pressure = 0. ! THIS COULD WORK + real :: pressure = 0. ! THIS COULD WORK [P] = erg/cm^3 = 1E-5 Pa (1 microbar) real :: aparam = 0.5 , bparam = 1.5 , alpha = 5. , beta = 5. real :: A = 2.67e11 , B = 2.67e11 ! erg/cm^3 real :: u0 = 4.87e12 ! erg/g @@ -59,7 +59,7 @@ subroutine equationofstate_tillotson(rho,u,pressure,spsound,gamma) mu_t = eta - 1. neu = (1. / eta) - 1. ! Kegerreis et al. 2020 omega = (u / (u0*(eta**2))) + 1. - spsoundmin = sqrt(A / rho) ! wave speed estimate Benz and Asphaug + spsoundmin = sqrt(A / rho0) ! wave speed estimate Benz and Asphaug ! gamma = 2./3. ! Brundage 2013 eq (2) cold expanded and compressed states @@ -84,100 +84,59 @@ subroutine equationofstate_tillotson(rho,u,pressure,spsound,gamma) endif spsound3 = sqrt( ((pressure3/rho)*( 1. + aparam + ((bparam / omega )*exp((-1.*alpha)*(neu**2))))) + & - ( (((bparam*rho*u)/((omega**2)*(eta**2))) * ((1. /(u0*rho)) * ((2.*u) - (pressure3/rho)) & + (( (((bparam*rho*u)/((omega**2)*(eta**2))) * ((1. /(u0*rho)) * ((2.*u) - (pressure3/rho)) & + ((2.*alpha*neu*omega)/rho0))) + ((A/rho0)*(1. + ((mu_t/(eta**2))*((beta+(2.*alpha*neu) - eta))) & - * exp((-1.*beta)*neu) )*exp((-1.*alpha)*(neu**2)) ))) + * exp((-1.*beta)*neu) )))*exp((-1.*alpha)*(neu**2)) )) if (spsound3 < spsoundmin) then spsound3 = spsoundmin endif spsound2h = sqrt( ((pressure2 / rho) * (1. + aparam + (bparam / omega ))) + & (((bparam*(omega - 1.)) / (omega**2)) * ((2.*u) - (pressure2/rho))) + & - ((1./rho)*(A + (B*((eta**2) - 1.)))) - (((B*(mu_t**2))/rho)*(1+aparam+(bparam/(omega**2))))) + (A/rho)) ! spsound2h = sqrt( ((pressure2 / rho) * (1. + aparam + (bparam / omega ))) + & ! (((bparam*(omega - 1.)) / (omega**2)) * ((2.*u) - pressure2/rho)) + & ! ((1./rho)*(A + B*((eta**2) - 1.))) - ((((2.*B)/(rho0**2))*rho) - (2.*B/rho0))) if (spsound2h < spsoundmin) then - spsound2h = spsoundmin + spsound2h = spsoundmin endif - -! print*,' in eos - calc pres from rho and u' -! print*,' rho in = ',rho,' u in = ',u -! print*,'in low energy state, p = ',pressure,' pres in code units = ',pressure/unit_pressure ! rho < rho_iv ; u < u_cv = ISSUES -! print*,'in hybrid state, p = ',pressure,' pres in code units = ',pressure/unit_pressure ! ISSUES -! print*,'in compressed state, p = ',pressure,' pres in code units = ',pressure/unit_pressure ! rho > rho0 ; u >= u0 -! print*,'in hot expanded state, p = ',pressure,' pres in code units = ',pressure/unit_pressure ! rho_iv <= rho < rho0 ; u >= u_cv - -! if (rho < rho_iv) then ! low energy expanded + spsoundh = sqrt(( ((u-u_iv)*(spsound3)**2) + ((u_cv-u)*(spsound2)**2) ) / (u_cv - u_iv)) + if (spsoundh < spsoundmin) then + spsoundh = spsoundmin + endif + +! if (rho < rho_iv) then ! low energy expanded = [ISSUES!!] ! pressure = pressure2 - (B*(mu_t**2)) ! spsound = spsound2h ! print*,' in eos - calc pres from rho and u' ! print*,' rho in = ',rho,' u in = ',u ! print*,'in low energy state, p = ',pressure,' pres in code units = ',pressure/unit_pressure ! else ! (rho >= rho_iv) - if (rho >= rho0) then ! compressed and cold expanded +if (rho >= rho0) then ! compressed and cold expanded pressure = pressure2 spsound = spsound2 ! print*,' in eos - calc pres from rho and u' ! print*,' rho in = ',rho,' u in = ',u ! print*,' spsound = ',spsound ! print*,'in compressed state, p = ',pressure,' pres in code units = ',pressure/unit_pressure - else ! (rho_iv <= rho < rho0) -! if (u >= u_cv) then ! hot / vapourised expanded - pressure = pressure3 - spsound = spsound3 - ! print*,' in eos - calc pres from rho and u' - ! print*,' rho in = ',rho,' u in = ',u - ! print*,' spsound = ',spsound - ! print*,'in hot expanded state, p = ',pressure,' pres in code units = ',pressure/unit_pressure - endif +else ! (rho_iv <= rho < rho0) + if (u < u_cv) then ! (u > u_iv) hybrid + pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) + spsound = spsoundh + ! print*,' in eos - calc pres from rho and u' + ! print*,' rho in = ',rho,' u in = ',u + ! print*,'in hybrid state, p = ',pressure,' pres in code units = ',pressure/unit_pressure + else ! if (u >= u_cv) then ! hot / vapourised expanded + pressure = pressure3 + spsound = spsound3 + ! print*,' in eos - calc pres from rho and u' + ! print*,' rho in = ',rho,' u in = ',u + ! print*,' spsound = ',spsound + ! print*,'in hot expanded state, p = ',pressure,' pres in code units = ',pressure/unit_pressure + endif +endif ! print*,'leaving eos (pres calc) ' -! else ! (u < u_cv) -! if (u <= u_iv) then ! cold expanded -! pressure = pressure2 -! spsound = spsound2 -! print*,' in eos - calc pres from rho and u' -! print*,' rho in = ',rho,' u in = ',u -! print*,'in cold expanded state, p = ',pressure,' pres in code units = ',pressure/unit_pressure -! else ! (u > u_iv) hybrid -! pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) -! spsoundh = sqrt(( ((u-u_iv)*(spsound3)**2) + ((u_cv-u)*(spsound2)**2) ) / (u_cv - u_iv)) -! if (spsoundh < spsoundmin) then -! spsound = spsoundmin -! endif -! print*,' in eos - calc pres from rho and u' -! print*,' rho in = ',rho,' u in = ',u -! print*,'in hybrid state, p = ',pressure,' pres in code units = ',pressure/unit_pressure -! endif -! endif -! endif -! endif - -! if (rho < rho0) then -! if (rho >= rho_iv) then ! u <= u_iv cold expanded -! if (u < u_cv) then ! hybrid state -! pressure = pressure2 -! spsound = spsound2 -! elseif (u > u_iv) then ! hybrid state -! pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) -! spsoundh = sqrt(( ((u-u_iv)*(spsound3)**2) + ((u_cv-u)*(spsound2)**2) ) / (u_cv - u_iv)) -! if (spsoundh < spsoundmin) then -! spsound = spsoundmin -! endif -! print*,'in hybrid state, p = ',pressure -! else ! (u >= u_cv) vaporised expanded state -! pressure = pressure3 -! spsound = spsound3 -! endif -! else ! (rho < rho_iv) then (u = rho0) compressed state -! pressure = pressure2 -! spsound = spsound2 -! endif end subroutine equationofstate_tillotson @@ -218,28 +177,8 @@ subroutine calc_uT_from_rhoP_tillotson(rho,pres,temp,u,ierr) expterm2 = exp(-alpha*(neu**2)) expterm = exp(-beta*neu) -! pres_iv = (( aparam + ( bparam / ((u_iv / (u0*(rho_iv/rho0)**2) + 1.)) ) ) * rho_iv*u_iv) & - ! + A*((rho_iv/rho0)-1.) + B*((rho_iv/rho0)-1.)**2 -! pres_cv = (aparam*rho_cv*u_cv) + ( ((bparam*rho_cv*u_cv)/((u_cv / (u0*(rho_cv/rho0)**2) + 1.))) & - ! + A*((rho_cv/rho0) - 1.)*exp((-1.*beta)*((rho0/rho_cv) - 1.))) * exp((-1.*alpha)*(((rho0/rho_cv) - 1.)**2)) - -! quada = 1. -! X = (pres - (A*mu_t))/rho -! quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) -! quadc = -((u0*(eta**2)*X)/aparam) -! print*,' X = ',X,' a = ',quada,' b = ',quadb,' c = ',quadc -! call cubicsolve(0.,quada,quadb,quadc,usol,nsol) -! print*,' Solutions : ',usol,' # of sols : ',nsol -! u = maxval(usol(1:nsol)) -! print*,' in eos - calc u from rho and pres' -! print*,' rho in = ',rho,' pres in = ',pres -! print*,'in low energy state, u = ',u,' u in code units = ',u/unit_ergg ! rho < rho_iv ; pre < pres_cv = ISSUES -! print*,'in hot expanded state, u = ',u,' u in code units = ',u/unit_ergg ! rho_iv <= rho < rho0 ; pres >= pres_cv -! print*,'in compressed state, u = ',u,' u in code units = ',u/unit_ergg ! rho > rho0 ; pres >= pres0 -! print*,'leaving eos (u calc) ' - -! if (rho < rho_iv) then ! low energy expanded +! if (rho < rho_iv) then ! low energy expanded = [ISSUES!!] ! quada = 1. ! X = (pres - (A*mu_t))/rho ! quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) @@ -249,14 +188,14 @@ subroutine calc_uT_from_rhoP_tillotson(rho,pres,temp,u,ierr) ! print*,' in eos - calc u from rho and pres' ! print*,' rho in = ',rho,' pres in = ',pres ! print*,'in low energy state, u = ',u,' u in code units = ',u/unit_ergg -! else ! (rho >= rho_iv) + if (rho >= rho0) then ! compressed and cold expanded quada = 1. X = pres - (A*mu_t) - (B*(mu_t**2)) quadb = (((u0*(eta**2)*rho)*(aparam + bparam)) - X)/(aparam*rho) quadc = -1.*((X*u0*(eta**2))/(rho*aparam)) call cubicsolve(0.,quada,quadb,quadc,usol,nsol) - ! u = maxval(usol(1:nsol)) + u = maxval(usol(1:nsol)) ! print*,' in eos - calc u from rho and pres' ! print*,' rho in = ',rho,' pres in = ',pres ! print*,'in compressed state, u = ',u,' u in code units = ',u/unit_ergg @@ -272,7 +211,7 @@ subroutine calc_uT_from_rhoP_tillotson(rho,pres,temp,u,ierr) ! print*,' rho in = ',rho,' pres in = ',pres ! print*,'in hot expanded state, u = ',u,' u in code units = ',u/unit_ergg endif -! print*,'leaving eos (pres calc) ' +! print*,'leaving eos (u calc) ' ! else ! (pres < pres_cv) ! if (pres <= pres_iv) then ! cold expanded ! quada = 1. @@ -293,36 +232,6 @@ subroutine calc_uT_from_rhoP_tillotson(rho,pres,temp,u,ierr) ! endif ! endif ! endif -! print*,'leaving eos ' -! if (rho >= rho0) then ! compressed state -! quada = 1. -! X = (pres - (A*mu_t) - (B*(mu_t**2)))/rho -! quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) -! quadc = -((u0*(eta**2)*X)/aparam) -! call cubicsolve(0.,quada,quadb,quadc,usol,nsol) -! u = maxval(usol(1:nsol)) -! elseif (rho >= rho_iv) then ! cold expanded -! quada = 1. -! X = (pres - (A*mu_t) - (B*(mu_t**2)))/rho -! quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) -! quadc = -((u0*(eta**2)*X)/aparam) -! call cubicsolve(0.,quada,quadb,quadc,usol,nsol) -! u = maxval(usol(1:nsol)) -! elseif (rho < rho_iv) then ! rho < rho_iv Low energy state -! quada = 1. -! X = (pres - (A*mu_t))/rho -! quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) -! quadc = -((u0*(eta**2)*X)/aparam) -! call cubicsolve(0.,quada,quadb,quadc,usol,nsol) -! u = maxval(usol(1:nsol)) -! else ! rho < rho0 vaporised expanded state -! quada = 1. -! X = (pres - (A*mu_t*expterm*expterm2))/rho -! quadb = ((bparam*u0*(eta**2)*expterm2)/aparam)+((u0*(eta**2)))-(X/aparam) -! quadc = -((X*u0*(eta**2))/aparam) -! call cubicsolve(0.,quada,quadb,quadc,usol,nsol) -! u = maxval(usol(1:nsol)) -! endif temp = 300. ierr = 0 From 8318ad223fee9b67d28637c1999fb13e36df87cf Mon Sep 17 00:00:00 2001 From: Amber Tilly Date: Mon, 14 Oct 2024 08:09:26 +1100 Subject: [PATCH 29/29] (tillotson) update set star utils and eos to ensure zero initial internal energy --- src/main/eos_tillotson.f90 | 108 +++++++++++++++++------------------ src/setup/set_star_utils.f90 | 2 + 2 files changed, 56 insertions(+), 54 deletions(-) diff --git a/src/main/eos_tillotson.f90 b/src/main/eos_tillotson.f90 index 83137308c..5960f97a1 100644 --- a/src/main/eos_tillotson.f90 +++ b/src/main/eos_tillotson.f90 @@ -94,9 +94,6 @@ subroutine equationofstate_tillotson(rho,u,pressure,spsound,gamma) spsound2h = sqrt( ((pressure2 / rho) * (1. + aparam + (bparam / omega ))) + & (((bparam*(omega - 1.)) / (omega**2)) * ((2.*u) - (pressure2/rho))) + & (A/rho)) -! spsound2h = sqrt( ((pressure2 / rho) * (1. + aparam + (bparam / omega ))) + & -! (((bparam*(omega - 1.)) / (omega**2)) * ((2.*u) - pressure2/rho)) + & -! ((1./rho)*(A + B*((eta**2) - 1.))) - ((((2.*B)/(rho0**2))*rho) - (2.*B/rho0))) if (spsound2h < spsoundmin) then spsound2h = spsoundmin endif @@ -106,34 +103,36 @@ subroutine equationofstate_tillotson(rho,u,pressure,spsound,gamma) spsoundh = spsoundmin endif -! if (rho < rho_iv) then ! low energy expanded = [ISSUES!!] + if (rho < rho_iv) then ! low energy expanded = [ISSUES!!] + pressure = 1.0 ! 1.E-16 GPa [10^7 erg = 1Pa ; Giga = 10^9] + spsound = spsoundmin ! pressure = pressure2 - (B*(mu_t**2)) -! spsound = spsound2h ! print*,' in eos - calc pres from rho and u' ! print*,' rho in = ',rho,' u in = ',u ! print*,'in low energy state, p = ',pressure,' pres in code units = ',pressure/unit_pressure -! else ! (rho >= rho_iv) -if (rho >= rho0) then ! compressed and cold expanded - pressure = pressure2 - spsound = spsound2 - ! print*,' in eos - calc pres from rho and u' - ! print*,' rho in = ',rho,' u in = ',u - ! print*,' spsound = ',spsound - ! print*,'in compressed state, p = ',pressure,' pres in code units = ',pressure/unit_pressure -else ! (rho_iv <= rho < rho0) - if (u < u_cv) then ! (u > u_iv) hybrid - pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) - spsound = spsoundh - ! print*,' in eos - calc pres from rho and u' - ! print*,' rho in = ',rho,' u in = ',u - ! print*,'in hybrid state, p = ',pressure,' pres in code units = ',pressure/unit_pressure - else ! if (u >= u_cv) then ! hot / vapourised expanded - pressure = pressure3 - spsound = spsound3 + else ! (rho >= rho_iv) + if (rho >= rho0) then ! compressed and cold expanded + pressure = pressure2 + spsound = spsound2 ! print*,' in eos - calc pres from rho and u' ! print*,' rho in = ',rho,' u in = ',u ! print*,' spsound = ',spsound - ! print*,'in hot expanded state, p = ',pressure,' pres in code units = ',pressure/unit_pressure + ! print*,'in compressed state, p = ',pressure,' pres in code units = ',pressure/unit_pressure + else ! (rho_iv <= rho < rho0) + if (u < u_cv) then ! (u > u_iv) hybrid + pressure = ( (u - u_iv)*pressure3 + (u_cv - u)*pressure2 ) / (u_cv - u_iv) + spsound = spsoundh + ! print*,' in eos - calc pres from rho and u' + ! print*,' rho in = ',rho,' u in = ',u + ! print*,'in hybrid state, p = ',pressure,' pres in code units = ',pressure/unit_pressure + else ! if (u >= u_cv) then ! hot / vapourised expanded + pressure = pressure3 + spsound = spsound3 + ! print*,' in eos - calc pres from rho and u' + ! print*,' rho in = ',rho,' u in = ',u + ! print*,' spsound = ',spsound + ! print*,'in hot expanded state, p = ',pressure,' pres in code units = ',pressure/unit_pressure + endif endif endif ! print*,'leaving eos (pres calc) ' @@ -177,40 +176,41 @@ subroutine calc_uT_from_rhoP_tillotson(rho,pres,temp,u,ierr) expterm2 = exp(-alpha*(neu**2)) expterm = exp(-beta*neu) - -! if (rho < rho_iv) then ! low energy expanded = [ISSUES!!] -! quada = 1. -! X = (pres - (A*mu_t))/rho -! quadb = ((-X/aparam)+(u0*(eta**2))+((bparam*u0*(eta**2))/aparam)) -! quadc = -((u0*(eta**2)*X)/aparam) -! call cubicsolve(0.,quada,quadb,quadc,usol,nsol) -! u = maxval(usol(1:nsol)) -! print*,' in eos - calc u from rho and pres' -! print*,' rho in = ',rho,' pres in = ',pres -! print*,'in low energy state, u = ',u,' u in code units = ',u/unit_ergg - - if (rho >= rho0) then ! compressed and cold expanded - quada = 1. - X = pres - (A*mu_t) - (B*(mu_t**2)) + + if (rho < rho_iv) then ! low energy expanded = [ISSUES!!] + quada = 1. ! min pressure = 1. -> min u + X = 1. - (A*mu_t) quadb = (((u0*(eta**2)*rho)*(aparam + bparam)) - X)/(aparam*rho) quadc = -1.*((X*u0*(eta**2))/(rho*aparam)) call cubicsolve(0.,quada,quadb,quadc,usol,nsol) u = maxval(usol(1:nsol)) - ! print*,' in eos - calc u from rho and pres' - ! print*,' rho in = ',rho,' pres in = ',pres - ! print*,'in compressed state, u = ',u,' u in code units = ',u/unit_ergg - else ! (rho_iv <= rho < rho0) -! if (pres >= pres_cv) then ! hot / vapourised expanded - quada = 1. - X = (pres - (A*mu_t*expterm*expterm2)) - quadb = ((aparam*rho*u0*(eta**2)) + (bparam*rho*u0*(eta**2)*expterm2) -X)/(rho*aparam) ! (u0*(eta**2))+((bparam*u0*(eta**2)*expterm2)/aparam) - (X/aparam) - quadc = -1.*(X*u0*(eta**2))/(aparam*rho) ! -1.*((X*u0*(eta**2))/aparam) - call cubicsolve(0.,quada,quadb,quadc,usol,nsol) - u = maxval(usol(1:nsol)) - ! print*,' in eos - calc u from rho and pres' - ! print*,' rho in = ',rho,' pres in = ',pres - ! print*,'in hot expanded state, u = ',u,' u in code units = ',u/unit_ergg - endif +! print*,' in eos - calc u from rho and pres' +! print*,' rho in = ',rho,' pres in = ',pres +! print*,'in low energy state, u = ',u,' u in code units = ',u/unit_ergg + else + if (rho >= rho0) then ! compressed and cold expanded + quada = 1. + X = pres - (A*mu_t) - (B*(mu_t**2)) + quadb = (((u0*(eta**2)*rho)*(aparam + bparam)) - X)/(aparam*rho) + quadc = -1.*((X*u0*(eta**2))/(rho*aparam)) + call cubicsolve(0.,quada,quadb,quadc,usol,nsol) + u = maxval(usol(1:nsol)) + ! print*,' in eos - calc u from rho and pres' + ! print*,' rho in = ',rho,' pres in = ',pres + ! print*,'in compressed state, u = ',u,' u in code units = ',u/unit_ergg + else ! (rho_iv <= rho < rho0) + ! if (pres >= pres_cv) then ! hot / vapourised expanded + quada = 1. + X = (pres - (A*mu_t*expterm*expterm2)) + quadb = ((aparam*rho*u0*(eta**2)) + (bparam*rho*u0*(eta**2)*expterm2) -X)/(rho*aparam) ! (u0*(eta**2))+((bparam*u0*(eta**2)*expterm2)/aparam) - (X/aparam) + quadc = -1.*(X*u0*(eta**2))/(aparam*rho) ! -1.*((X*u0*(eta**2))/aparam) + call cubicsolve(0.,quada,quadb,quadc,usol,nsol) + u = maxval(usol(1:nsol)) + ! print*,' in eos - calc u from rho and pres' + ! print*,' rho in = ',rho,' pres in = ',pres + ! print*,'in hot expanded state, u = ',u,' u in code units = ',u/unit_ergg + endif +endif ! print*,'leaving eos (u calc) ' ! else ! (pres < pres_cv) ! if (pres <= pres_iv) then ! cold expanded diff --git a/src/setup/set_star_utils.f90 b/src/setup/set_star_utils.f90 index 4be4dbe7e..8ca396330 100644 --- a/src/setup/set_star_utils.f90 +++ b/src/setup/set_star_utils.f90 @@ -419,6 +419,8 @@ subroutine set_star_thermalenergy(ieos,den,pres,r,npts,npart,xyzh,vxyzu,rad,eos_ endif select case(ieos) + case(23) ! Tillotson + vxyzu(4,i) = 0. case(16) ! Shen EoS vxyzu(4,i) = initialtemp case(15) ! Helmholtz EoS