From ae14e6761eb706af56610d5cbf02967716ef369c Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Thu, 16 May 2024 11:42:57 +0200 Subject: [PATCH 01/19] (CE-analysis) mistake---rad array should be used where radprop used --- src/main/ionization.f90 | 6 ++-- src/utils/analysis_common_envelope.f90 | 46 +++++++++++++------------- 2 files changed, 26 insertions(+), 26 deletions(-) diff --git a/src/main/ionization.f90 b/src/main/ionization.f90 index 2ebad8398..4823171f7 100644 --- a/src/main/ionization.f90 +++ b/src/main/ionization.f90 @@ -338,7 +338,7 @@ end subroutine get_erec_components ! gas particle. Inputs and outputs in code units !+ !---------------------------------------------------------------- -subroutine calc_thermal_energy(particlemass,ieos,xyzh,vxyzu,presi,tempi,ethi,radprop) +subroutine calc_thermal_energy(particlemass,ieos,xyzh,vxyzu,presi,tempi,ethi,rad) use dim, only:do_radiation use part, only:rhoh,iradxi use eos_idealplusrad, only:get_idealgasplusrad_tempfrompres,get_idealplusrad_enfromtemp @@ -346,7 +346,7 @@ subroutine calc_thermal_energy(particlemass,ieos,xyzh,vxyzu,presi,tempi,ethi,rad use units, only:unit_density,unit_pressure,unit_ergg,unit_pressure integer, intent(in) :: ieos real, intent(in) :: particlemass,presi,tempi,xyzh(4),vxyzu(4) - real, intent(in), optional :: radprop(:) + real, intent(in), optional :: rad(:) real, intent(out) :: ethi real :: hi,densi_cgs,mui @@ -359,7 +359,7 @@ subroutine calc_thermal_energy(particlemass,ieos,xyzh,vxyzu,presi,tempi,ethi,rad ethi = particlemass * ethi / unit_ergg case default ! assuming internal energy = thermal energy ethi = particlemass * vxyzu(4) - if (do_radiation) ethi = ethi + particlemass*radprop(iradxi) + if (do_radiation) ethi = ethi + particlemass*rad(iradxi) end select end subroutine calc_thermal_energy diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index a7ec86cff..f882e923a 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -22,7 +22,7 @@ module analysis use part, only:xyzmh_ptmass,vxyz_ptmass,nptmass,poten,ihsoft,ihacc,& rhoh,nsinkproperties,maxvxyzu,maxptmass,isdead_or_accreted,& - radprop + rad,radprop use dim, only:do_radiation use units, only:print_units,umass,utime,udist,unit_ergg,unit_density,& unit_pressure,unit_velocity,unit_Bfield,unit_energ @@ -521,7 +521,7 @@ end subroutine m_vs_t !+ !---------------------------------------------------------------- subroutine bound_mass(time,npart,particlemass,xyzh,vxyzu) - use part, only:eos_vars,itemp,radprop + use part, only:eos_vars,itemp use ptmass, only:get_accel_sink_gas use vectorutils, only:cross_product3D integer, intent(in) :: npart @@ -593,13 +593,13 @@ subroutine bound_mass(time,npart,particlemass,xyzh,vxyzu) do i = 1,npart if (.not. isdead_or_accreted(xyzh(4,i))) then - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) call get_accel_sink_gas(nptmass,xyzh(1,i),xyzh(2,i),xyzh(3,i),xyzh(4,i),xyzmh_ptmass,dum1,dum2,dum3,phii) rhopart = rhoh(xyzh(4,i), particlemass) tempi = eos_vars(itemp,i) call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) call cross_product3D(xyzh(1:3,i), particlemass * vxyzu(1:3,i), rcrossmv) ! Angular momentum w.r.t. CoM - call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,tempi,ethi,radprop(:,i)) + call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,tempi,ethi,rad(:,i)) etoti = ekini + epoti + ethi ! Overwrite etoti outputted by calc_gas_energies to use ethi instead of einti else ! Output 0 for quantities pertaining to accreted particles @@ -755,7 +755,7 @@ subroutine calculate_energies(time,npart,particlemass,xyzh,vxyzu) jz = rcrossmv(3) encomp(ijz_tot) = encomp(ijz_tot) + jz - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) encomp(ipot_ps) = encomp(ipot_ps) + particlemass * phii @@ -1001,7 +1001,7 @@ subroutine roche_lobe_values(time,npart,particlemass,xyzh,vxyzu) call orbit_com(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptmass,com_xyz,com_vxyz) do i=1,npart - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) sep1 = separation(xyzmh_ptmass(1:3,1),xyzh(1:3,i)) sep2 = separation(xyzmh_ptmass(1:3,2),xyzh(1:3,i)) @@ -1177,7 +1177,7 @@ subroutine star_stabilisation_suite(time,npart,particlemass,xyzh,vxyzu) rhopart = rhoh(xyzh(4,i), particlemass) totvol = totvol + particlemass / rhopart ! Sum "volume" of all particles virialpart = virialpart + particlemass * ( dot_product(fxyzu(1:3,i),xyzh(1:3,i)) + dot_product(vxyzu(1:3,i),vxyzu(1:3,i)) ) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) totekin = totekin + ekini totepot = totepot + 0.5*epoti ! Factor of 1/2 to correct for double counting if (rhopart > rho_surface) then @@ -1417,7 +1417,7 @@ subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu) case(1,9) ! Total energy (kin + pot + therm) rhopart = rhoh(xyzh(4,i), particlemass) call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum1) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum1) if (quantities_to_calculate(k)==1) then call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),ethi) quant(k,i) = (ekini + epoti + ethi) / particlemass ! Specific energy @@ -1828,7 +1828,7 @@ subroutine recombination_tau(time,npart,particlemass,xyzh,vxyzu) call get_eos_kappa_mesa(rho_part(i)*unit_density,eos_vars(itemp,i),kappa,kappat,kappar) kappa_part(i) = kappa ! In cgs units call ionisation_fraction(rho_part(i)*unit_density,eos_vars(itemp,i),X_in,1.-X_in-Z_in,xh0,xh1,xhe0,xhe1,xhe2) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) ! Calculate total energy + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) ! Calculate total energy call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rho_part(i),eos_vars(itemp,i),ethi) etoti = ekini + epoti + ethi if ((xh0 > recomb_th) .and. (.not. prev_recombined(i)) .and. (etoti < 0.)) then ! Recombination event and particle is still bound @@ -1907,7 +1907,7 @@ subroutine energy_hist(time,npart,particlemass,xyzh,vxyzu) do i=1,npart rhopart = rhoh(xyzh(4,i), particlemass) call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) if (ieos==10 .or. ieos==20) then call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),ethi) else @@ -2054,7 +2054,7 @@ subroutine energy_profile(time,npart,particlemass,xyzh,vxyzu) call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) select case (iquantity) case(1) ! Energy - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),ethi) quant(i,1) = ekini + epoti + ethi case(2) ! Entropy @@ -2070,7 +2070,7 @@ subroutine energy_profile(time,npart,particlemass,xyzh,vxyzu) quant(i,1) = entropy(rhopart*unit_density,ponrhoi*rhopart*unit_pressure,mu,ientropy,ierr=ierr) endif case(3) ! Bernoulli energy (per unit mass) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) quant(i,1) = 0.5*dot_product(vxyzu(1:3,i),vxyzu(1:3,i)) + ponrhoi + vxyzu(4,i) + epoti/particlemass ! 1/2 v^2 + P/rho + phi case(4) ! Ion fraction call ionisation_fraction(rhopart*unit_density,eos_vars(itemp,i),X_in,1.-X_in-Z_in,xh0,xh1,xhe0,xhe1,xhe2) @@ -2202,7 +2202,7 @@ subroutine velocity_histogram(time,num,npart,particlemass,xyzh,vxyzu) do i = 1,npart rhopart = rhoh(xyzh(4,i), particlemass) call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),ethi) vr(i) = dot_product(xyzh(1:3,i),vxyzu(1:3,i)) / sqrt(dot_product(xyzh(1:3,i),xyzh(1:3,i))) @@ -2510,7 +2510,7 @@ subroutine unbound_profiles(time,num,npart,particlemass,xyzh,vxyzu) if (.not. isdead_or_accreted(xyzh(4,i))) then rhopart = rhoh(xyzh(4,i), particlemass) call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,tempi,ethi) etoti = ekini + epoti + ethi @@ -2618,7 +2618,7 @@ subroutine unbound_ionfrac(time,npart,particlemass,xyzh,vxyzu) do i=1,npart rhopart = rhoh(xyzh(4,i), particlemass) call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,tempi,ethi) etoti = ekini + epoti + ethi @@ -2688,7 +2688,7 @@ subroutine unbound_temp(time,npart,particlemass,xyzh,vxyzu) do i=1,npart rhopart = rhoh(xyzh(4,i), particlemass) call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),eos_vars(itemp,i),vxyzu(4,i)) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),ethi) etoti = ekini + epoti + ethi @@ -2758,7 +2758,7 @@ subroutine recombination_stats(time,num,npart,particlemass,xyzh,vxyzu) ! Calculate total energy rhopart = rhoh(xyzh(4,i), particlemass) call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),ethi) etoti = ekini + epoti + ethi @@ -3004,7 +3004,7 @@ subroutine bound_unbound_thermo(time,npart,particlemass,xyzh,vxyzu) call compute_energies(time) do i=1,npart - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) rhopart = rhoh(xyzh(4,i), particlemass) @@ -3350,7 +3350,7 @@ subroutine J_E_plane(num,npart,particlemass,xyzh,vxyzu) call get_centreofmass(com_xyz,com_vxyz,npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptmass) do i=1,npart - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),xyzmh_ptmass,dum1,dum2,dum3,dum4,etoti) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,dum1,dum2,dum3,dum4,etoti) data(1,i) = etoti call cross_product3D(xyzh(1:3,i)-xyzmh_ptmass(1:3,1), vxyzu(1:3,i)-vxyz_ptmass(1:3,1), angmom_core) data(5:7,i) = angmom_core @@ -3754,14 +3754,14 @@ end subroutine get_gas_omega ! and internal energy of a gas particle. !+ !---------------------------------------------------------------- -subroutine calc_gas_energies(particlemass,poten,xyzh,vxyzu,radprop,xyzmh_ptmass,phii,epoti,ekini,einti,etoti) +subroutine calc_gas_energies(particlemass,poten,xyzh,vxyzu,rad,xyzmh_ptmass,phii,epoti,ekini,einti,etoti) ! Warning: Do not sum epoti or etoti as it is to obtain a total energy; this would not give the correct ! total energy due to complications related to double-counting. use ptmass, only:get_accel_sink_gas use part, only:nptmass,iradxi real, intent(in) :: particlemass real(4), intent(in) :: poten - real, intent(in) :: xyzh(:),vxyzu(:),radprop(:) + real, intent(in) :: xyzh(:),vxyzu(:),rad(:) real, dimension(5,nptmass), intent(in) :: xyzmh_ptmass real, intent(out) :: phii,epoti,ekini,einti,etoti real :: fxi,fyi,fzi @@ -3773,7 +3773,7 @@ subroutine calc_gas_energies(particlemass,poten,xyzh,vxyzu,radprop,xyzmh_ptmass, epoti = 2.*poten + particlemass * phii ! For individual particles, need to multiply 2 to poten to get \sum_j G*mi*mj/r ekini = particlemass * 0.5 * dot_product(vxyzu(1:3),vxyzu(1:3)) einti = particlemass * vxyzu(4) - if (do_radiation) einti = einti + particlemass * radprop(iradxi) + if (do_radiation) einti = einti + particlemass * rad(iradxi) etoti = epoti + ekini + einti end subroutine calc_gas_energies @@ -3900,7 +3900,7 @@ subroutine stellar_profile(time,ncols,particlemass,npart,xyzh,vxyzu,profile,simp kappa = 1. endif - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),& + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),& xyzmh_ptmass,phii,epoti,ekini,einti,etoti) call ionisation_fraction(rhopart*unit_density,temp,X_in,1.-X_in-Z_in,xh0,xh1,xhe0,xhe1,xhe2) From 45bbf52406ecc498eb4371080d738599c1dac70f Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Thu, 16 May 2024 11:44:20 +0200 Subject: [PATCH 02/19] (CE-analysis) particle tracking for radiation quantities --- src/utils/analysis_common_envelope.f90 | 64 +++++++++++++++++--------- 1 file changed, 41 insertions(+), 23 deletions(-) diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index f882e923a..1c933a570 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -154,7 +154,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) case(19) ! Optical depth profile call tau_profile(time,num,npart,particlemass,xyzh) case(20) ! Particle tracker - call track_particle(time,particlemass,xyzh,vxyzu) + call track_particle(time,npart,particlemass,xyzh,vxyzu) case(21) ! Unbound ion fractions call unbound_ionfrac(time,npart,particlemass,xyzh,vxyzu) case(22) ! Optical depth at recombination @@ -1561,31 +1561,40 @@ end subroutine eos_surfaces ! Particle tracker: Paint the life of a particle !+ !---------------------------------------------------------------- -subroutine track_particle(time,particlemass,xyzh,vxyzu) - use part, only:eos_vars,itemp - use eos, only:entropy +subroutine track_particle(time,npart,particlemass,xyzh,vxyzu) + use part, only:itemp,iradxi,ilambda + use eos, only:entropy + use radiation_utils, only:Trad_from_rhoxi use mesa_microphysics, only:getvalue_mesa - use ionization_mod, only:ionisation_fraction + use ionization_mod, only:ionisation_fraction real, intent(in) :: time,particlemass + integer, intent(in) :: npart real, intent(inout) :: xyzh(:,:),vxyzu(:,:) - integer, parameter :: nparttotrack=10,ncols=17 + integer, parameter :: nparttotrack=6,ncols=20 real :: r,v,rhopart,ponrhoi,Si,spsoundi,tempi,machi,xh0,xh1,xhe0,xhe1,xhe2,& - ekini,einti,epoti,ethi,etoti,dum,phii,pgas,mu + ekini,einti,epoti,ethi,eradi,etoti,dum,phii,pgas,mu,rho_cgs,Tradi,lambdai real, dimension(ncols) :: datatable character(len=17) :: filenames(nparttotrack),columns(ncols) integer :: i,k,partID(nparttotrack),ientropy,ierr - partID = (/ 1,2,3,4,5,6,7,8,9,10 /) + ! pid_orig is a map from current particle ID to original particle ID, if particles have been removed +! call initial_to_current_IDs(npart,pid_orig) + +! partID = (/ 1,2,3,4,5,6,7,8,9,10 /) + partID = (/ 359018, 1669237, 342811, 598910, 1690937, 285745 /) columns = (/ ' r',& ' v',& ' rho',& ' temp',& + ' Trad',& + ' lambda',& 'entropy',& 'spsound',& ' mach',& ' ekin',& ' epot',& ' eth',& + ' erad',& ' eint',& ' etot',& ' xHI',& @@ -1605,36 +1614,45 @@ subroutine track_particle(time,particlemass,xyzh,vxyzu) r = separation(xyzh(1:3,i),xyzmh_ptmass(1:3,1)) v = separation(vxyzu(1:3,i),vxyz_ptmass(1:3,1)) rhopart = rhoh(xyzh(4,i), particlemass) + rho_cgs = rhopart*unit_density call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) machi = v / spsoundi select case(ieos) case(2) - ientropy = 1 + if (do_radiation) then + ientropy = 2 + else + ientropy = 1 + endif case(10,12) ientropy = 2 case default ientropy = -1 end select - if (ieos==10) then - call getvalue_mesa(rhopart*unit_density,vxyzu(4,i)*unit_ergg,3,pgas,ierr) ! Get gas pressure - mu = rhopart*unit_density * Rg * eos_vars(itemp,i) / pgas + if (ieos==10) then ! get mu + call getvalue_mesa(rho_cgs,vxyzu(4,i)*unit_ergg,3,pgas,ierr) ! Get gas pressure + mu = rho_cgs * Rg * tempi / pgas else mu = gmw endif - ! MESA ENTROPY - Si = 0. - if (ieos==10) then - Si = entropy(rhopart*unit_density,ponrhoi*rhopart*unit_pressure,mu,3,vxyzu(4,i)*unit_ergg,ierr) + Tradi = 0. + lambdai = 0. + if (do_radiation) then + lambdai = radprop(ilambda,i) + Tradi = Trad_from_rhoxi(rhopart,rad(iradxi,i)) + Si = entropy(rho_cgs,ponrhoi*rhopart*unit_pressure,mu,ientropy,vxyzu(4,i)*unit_ergg,ierr,Trad_in=Tradi) + else + Si = entropy(rho_cgs,ponrhoi*rhopart*unit_pressure,mu,ientropy,vxyzu(4,i)*unit_ergg,ierr) endif - ! MESA ENTROPY - ! Si = entropy(rhopart*unit_density,ponrhoi*rhopart*unit_pressure,mu,ientropy,vxyzu(4,i)*unit_ergg,ierr) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) - call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),ethi) - etoti = ekini + epoti + ethi - call ionisation_fraction(rhopart*unit_density,eos_vars(itemp,i),X_in,1.-X_in-Z_in,xh0,xh1,xhe0,xhe1,xhe2) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,tempi,ethi,rad(:,i)) + etoti = ekini + epoti + ethi ! ethi includes radiation energy + eradi = particlemass*rad(iradxi,i) + call ionisation_fraction(rho_cgs,tempi,X_in,1.-X_in-Z_in,xh0,xh1,xhe0,xhe1,xhe2) ! Write file - datatable = (/ r,v,rhopart,eos_vars(itemp,i),Si,spsoundi,machi,ekini,epoti,ethi,einti,etoti,xh0,xh1,xhe0,xhe1,xhe2 /) + datatable = (/ r,v,rhopart,tempi,Tradi,lambdai,Si,spsoundi,machi,ekini,epoti,ethi,eradi,einti,etoti,& + xh0,xh1,xhe0,xhe1,xhe2 /) call write_time_file(trim(adjustl(filenames(k))),columns,time,datatable,ncols,dump_number) enddo From 0f691db9509b5612d0674fd37686dfd1d8fbfd2b Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Thu, 16 May 2024 11:46:17 +0200 Subject: [PATCH 03/19] (radiation) entropy function to retrieve radiation temperature for rad arrays --- src/main/eos.f90 | 31 +++++++++++++++++++------------ src/main/radiation_utils.f90 | 24 ++++++++++++------------ src/setup/set_star_utils.f90 | 4 ++-- src/tests/test_sedov.F90 | 2 +- 4 files changed, 34 insertions(+), 27 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index be0410681..67a791772 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -907,36 +907,43 @@ end subroutine calc_rho_from_PT ! up to an additive integration constant, from density and pressure. !+ !----------------------------------------------------------------------- -function entropy(rho,pres,mu_in,ientropy,eint_in,ierr) +function entropy(rho,pres,mu_in,ientropy,eint_in,ierr,T_in,Trad_in) use io, only:fatal,warning - use physcon, only:radconst,kb_on_mh + use physcon, only:radconst,kb_on_mh,Rg use eos_idealplusrad, only:get_idealgasplusrad_tempfrompres use eos_mesa, only:get_eos_eT_from_rhop_mesa use mesa_microphysics, only:getvalue_mesa real, intent(in) :: rho,pres,mu_in - real, intent(in), optional :: eint_in + real, intent(in), optional :: eint_in,T_in,Trad_in integer, intent(in) :: ientropy integer, intent(out), optional :: ierr - real :: mu,entropy,logentropy,temp,eint + real :: mu,entropy,logentropy,temp,Trad,eint if (present(ierr)) ierr=0 mu = mu_in - select case(ientropy) - case(1) ! Include only gas entropy (up to additive constants) - temp = pres * mu / (rho * kb_on_mh) - entropy = kb_on_mh / mu * log(temp**1.5/rho) + if (present(T_in)) then ! is gas temperature specified? + temp = T_in + else + temp = pres * mu / (rho * Rg) ! used as initial guess for case 2 + endif + select case(ientropy) + case(1) ! Only include gas contribution ! check temp if (temp < tiny(0.)) call warning('entropy','temperature = 0 will give minus infinity with s entropy') + entropy = kb_on_mh / mu * log(temp**1.5/rho) - case(2) ! Include both gas and radiation entropy (up to additive constants) - temp = pres * mu / (rho * kb_on_mh) ! Guess for temp + case(2) ! Include both gas and radiation contributions (up to additive constants) call get_idealgasplusrad_tempfrompres(pres,rho,mu,temp) ! First solve for temp from rho and pres - entropy = kb_on_mh / mu * log(temp**1.5/rho) + 4.*radconst*temp**3 / (3.*rho) - + if (present(Trad_in)) then + Trad = Trad_in + else + Trad = temp ! assume thermal equilibrium + endif ! check temp if (temp < tiny(0.)) call warning('entropy','temperature = 0 will give minus infinity with s entropy') + entropy = kb_on_mh / mu * log(temp**1.5/rho) + 4.*radconst*Trad**3 / (3.*rho) case(3) ! Get entropy from MESA tables if using MESA EoS if (ieos /= 10 .and. ieos /= 20) call fatal('eos','Using MESA tables to calculate S from rho and pres, but not using MESA EoS') diff --git a/src/main/radiation_utils.f90 b/src/main/radiation_utils.f90 index 0147e01c7..bb42061ea 100644 --- a/src/main/radiation_utils.f90 +++ b/src/main/radiation_utils.f90 @@ -23,8 +23,8 @@ module radiation_utils public :: get_rad_R public :: radiation_equation_of_state public :: T_from_Etot - public :: radE_from_Trad - public :: Trad_from_radE + public :: radxi_from_rhoT + public :: Trad_from_rhoxi public :: ugas_from_Tgas public :: Tgas_from_ugas public :: get_opacity @@ -136,29 +136,29 @@ end function T_from_Etot !--------------------------------------------------------- !+ -! get the radiation energy from the radiation temperature +! get specific radiation energy density from the radiation temperature !+ !--------------------------------------------------------- -real function radE_from_Trad(Trad) result(radE) +real function radxi_from_rhoT(rho,Trad) result(radxi) use units, only:get_radconst_code - real, intent(in) :: Trad + real, intent(in) :: Trad,rho - radE = Trad**4*get_radconst_code() + radxi = Trad**4*get_radconst_code()/rho -end function radE_from_Trad +end function radxi_from_rhoT !--------------------------------------------------------- !+ -! get the radiation temperature from the radiation energy per unit volume +! get the radiation temperature from specific radiation energy !+ !--------------------------------------------------------- -real function Trad_from_radE(radE) result(Trad) +real function Trad_from_rhoxi(rho,radxi) result(Trad) use units, only:get_radconst_code - real, intent(in) :: radE + real, intent(in) :: rho,radxi - Trad = (radE/get_radconst_code())**0.25 + Trad = (rho*radxi/get_radconst_code())**0.25 -end function Trad_from_radE +end function Trad_from_rhoxi !--------------------------------------------------------- !+ diff --git a/src/setup/set_star_utils.f90 b/src/setup/set_star_utils.f90 index 73a5d7017..bcdc754d6 100644 --- a/src/setup/set_star_utils.f90 +++ b/src/setup/set_star_utils.f90 @@ -382,7 +382,7 @@ subroutine set_star_thermalenergy(ieos,den,pres,r,npts,npart,xyzh,vxyzu,rad,eos_ relaxed,use_var_comp,initialtemp,npin) use part, only:do_radiation,rhoh,massoftype,igas,itemp,igasP,iX,iZ,imu,iradxi use eos, only:equationofstate,calc_temp_and_ene,gamma,gmw - use radiation_utils, only:ugas_from_Tgas,radE_from_Trad + use radiation_utils, only:ugas_from_Tgas,radxi_from_rhoT use table_utils, only:yinterp use units, only:unit_density,unit_ergg,unit_pressure integer, intent(in) :: ieos,npart,npts @@ -439,7 +439,7 @@ subroutine set_star_thermalenergy(ieos,den,pres,r,npts,npart,xyzh,vxyzu,rad,eos_ endif if (do_radiation) then vxyzu(4,i) = ugas_from_Tgas(tempi,gamma,gmw) - rad(iradxi,i) = radE_from_Trad(tempi)/densi + rad(iradxi,i) = radxi_from_rhoT(densi,tempi) else vxyzu(4,i) = eni / unit_ergg endif diff --git a/src/tests/test_sedov.F90 b/src/tests/test_sedov.F90 index d12efb34a..9e96630d7 100644 --- a/src/tests/test_sedov.F90 +++ b/src/tests/test_sedov.F90 @@ -54,7 +54,7 @@ subroutine test_sedov(ntests,npass) use mpidomain, only:i_belong use checkconserved, only:etot_in,angtot_in,totmom_in,mdust_in use radiation_utils, only:set_radiation_and_gas_temperature_equal,& - T_from_Etot,Tgas_from_ugas,ugas_from_Tgas,radE_from_Trad,Trad_from_radE + T_from_Etot,Tgas_from_ugas,ugas_from_Tgas use readwrite_dumps, only:write_fulldump use step_lf_global, only:init_step integer, intent(inout) :: ntests,npass From d857b4dd867235d5fa25e9e5a7294f9d62ddbe31 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Thu, 16 May 2024 14:32:50 +0200 Subject: [PATCH 04/19] (CE-analysis) revert changes to radE subroutines and correct radiation energy calculation --- src/main/ionization.f90 | 13 +++++++------ src/main/radiation_utils.f90 | 25 ++++++++++++------------- src/setup/set_star_utils.f90 | 4 ++-- src/utils/analysis_common_envelope.f90 | 13 ++++++++----- 4 files changed, 29 insertions(+), 26 deletions(-) diff --git a/src/main/ionization.f90 b/src/main/ionization.f90 index 4823171f7..d43c2fd10 100644 --- a/src/main/ionization.f90 +++ b/src/main/ionization.f90 @@ -348,18 +348,19 @@ subroutine calc_thermal_energy(particlemass,ieos,xyzh,vxyzu,presi,tempi,ethi,rad real, intent(in) :: particlemass,presi,tempi,xyzh(4),vxyzu(4) real, intent(in), optional :: rad(:) real, intent(out) :: ethi - real :: hi,densi_cgs,mui + real :: hi,densi_cgs,mui,rhoi + rhoi = rhoh(hi,particlemass) select case (ieos) case(10,20) ! calculate just gas + radiation thermal energy hi = xyzh(4) - densi_cgs = rhoh(hi,particlemass)*unit_density + densi_cgs = rhoi*unit_density mui = densi_cgs * Rg * tempi / (presi*unit_pressure - radconst * tempi**4 / 3.) ! Get mu from pres and temp call get_idealplusrad_enfromtemp(densi_cgs,tempi,mui,ethi) ethi = particlemass * ethi / unit_ergg case default ! assuming internal energy = thermal energy ethi = particlemass * vxyzu(4) - if (do_radiation) ethi = ethi + particlemass*rad(iradxi) + if (do_radiation) ethi = ethi + particlemass*rad(iradxi)/rhoi end select end subroutine calc_thermal_energy @@ -419,9 +420,9 @@ subroutine ionisation_fraction(dens,temp,X,Y,xh0,xh1,xhe0,xhe1,xhe2) xhe2g = xhe2g + dx(3) enddo - xh1 = xh1g * n / nh - xhe1 = xhe1g * n / nhe - xhe2 = xhe2g * n / nhe + xh1 = max(xh1g * n / nh,1.e-99) + xhe1 = max(xhe1g * n / nhe,1.e-99) + xhe2 = max(xhe2g * n / nhe,1.e-99) xh0 = ((nh/n) - xh1g) * n / nh xhe0 = ((nhe/n) - xhe1g - xhe2g) * n / nhe diff --git a/src/main/radiation_utils.f90 b/src/main/radiation_utils.f90 index bb42061ea..d54015bd6 100644 --- a/src/main/radiation_utils.f90 +++ b/src/main/radiation_utils.f90 @@ -23,8 +23,8 @@ module radiation_utils public :: get_rad_R public :: radiation_equation_of_state public :: T_from_Etot - public :: radxi_from_rhoT - public :: Trad_from_rhoxi + public :: radE_from_Trad + public :: Trad_from_radE public :: ugas_from_Tgas public :: Tgas_from_ugas public :: get_opacity @@ -136,29 +136,28 @@ end function T_from_Etot !--------------------------------------------------------- !+ -! get specific radiation energy density from the radiation temperature +! get the radiation energy from the radiation temperature !+ !--------------------------------------------------------- -real function radxi_from_rhoT(rho,Trad) result(radxi) +real function radE_from_Trad(Trad) result(radE) use units, only:get_radconst_code - real, intent(in) :: Trad,rho + real, intent(in) :: Trad - radxi = Trad**4*get_radconst_code()/rho + radE = Trad**4*get_radconst_code() -end function radxi_from_rhoT +end function radE_from_Trad !--------------------------------------------------------- !+ -! get the radiation temperature from specific radiation energy -!+ +! get the radiation temperature from the radiation energy per unit volume!+ !--------------------------------------------------------- -real function Trad_from_rhoxi(rho,radxi) result(Trad) +real function Trad_from_radE(radE) result(Trad) use units, only:get_radconst_code - real, intent(in) :: rho,radxi + real, intent(in) :: radE - Trad = (rho*radxi/get_radconst_code())**0.25 + Trad = (radE/get_radconst_code())**0.25 -end function Trad_from_rhoxi +end function Trad_from_radE !--------------------------------------------------------- !+ diff --git a/src/setup/set_star_utils.f90 b/src/setup/set_star_utils.f90 index bcdc754d6..73a5d7017 100644 --- a/src/setup/set_star_utils.f90 +++ b/src/setup/set_star_utils.f90 @@ -382,7 +382,7 @@ subroutine set_star_thermalenergy(ieos,den,pres,r,npts,npart,xyzh,vxyzu,rad,eos_ relaxed,use_var_comp,initialtemp,npin) use part, only:do_radiation,rhoh,massoftype,igas,itemp,igasP,iX,iZ,imu,iradxi use eos, only:equationofstate,calc_temp_and_ene,gamma,gmw - use radiation_utils, only:ugas_from_Tgas,radxi_from_rhoT + use radiation_utils, only:ugas_from_Tgas,radE_from_Trad use table_utils, only:yinterp use units, only:unit_density,unit_ergg,unit_pressure integer, intent(in) :: ieos,npart,npts @@ -439,7 +439,7 @@ subroutine set_star_thermalenergy(ieos,den,pres,r,npts,npart,xyzh,vxyzu,rad,eos_ endif if (do_radiation) then vxyzu(4,i) = ugas_from_Tgas(tempi,gamma,gmw) - rad(iradxi,i) = radxi_from_rhoT(densi,tempi) + rad(iradxi,i) = radE_from_Trad(tempi)/densi else vxyzu(4,i) = eni / unit_ergg endif diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index 1c933a570..77567ead8 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -1564,7 +1564,7 @@ end subroutine eos_surfaces subroutine track_particle(time,npart,particlemass,xyzh,vxyzu) use part, only:itemp,iradxi,ilambda use eos, only:entropy - use radiation_utils, only:Trad_from_rhoxi + use radiation_utils, only:Trad_from_radE use mesa_microphysics, only:getvalue_mesa use ionization_mod, only:ionisation_fraction real, intent(in) :: time,particlemass @@ -1639,7 +1639,7 @@ subroutine track_particle(time,npart,particlemass,xyzh,vxyzu) lambdai = 0. if (do_radiation) then lambdai = radprop(ilambda,i) - Tradi = Trad_from_rhoxi(rhopart,rad(iradxi,i)) + Tradi = Trad_from_radE(rad(iradxi,i)) Si = entropy(rho_cgs,ponrhoi*rhopart*unit_pressure,mu,ientropy,vxyzu(4,i)*unit_ergg,ierr,Trad_in=Tradi) else Si = entropy(rho_cgs,ponrhoi*rhopart*unit_pressure,mu,ientropy,vxyzu(4,i)*unit_ergg,ierr) @@ -1647,7 +1647,7 @@ subroutine track_particle(time,npart,particlemass,xyzh,vxyzu) call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,tempi,ethi,rad(:,i)) etoti = ekini + epoti + ethi ! ethi includes radiation energy - eradi = particlemass*rad(iradxi,i) + eradi = particlemass*rad(iradxi,i)/rhopart call ionisation_fraction(rho_cgs,tempi,X_in,1.-X_in-Z_in,xh0,xh1,xhe0,xhe1,xhe2) ! Write file @@ -3782,7 +3782,7 @@ subroutine calc_gas_energies(particlemass,poten,xyzh,vxyzu,rad,xyzmh_ptmass,phii real, intent(in) :: xyzh(:),vxyzu(:),rad(:) real, dimension(5,nptmass), intent(in) :: xyzmh_ptmass real, intent(out) :: phii,epoti,ekini,einti,etoti - real :: fxi,fyi,fzi + real :: fxi,fyi,fzi,rhopart phii = 0.0 @@ -3791,7 +3791,10 @@ subroutine calc_gas_energies(particlemass,poten,xyzh,vxyzu,rad,xyzmh_ptmass,phii epoti = 2.*poten + particlemass * phii ! For individual particles, need to multiply 2 to poten to get \sum_j G*mi*mj/r ekini = particlemass * 0.5 * dot_product(vxyzu(1:3),vxyzu(1:3)) einti = particlemass * vxyzu(4) - if (do_radiation) einti = einti + particlemass * rad(iradxi) + if (do_radiation) then + rhopart = rhoh(xyzh(4),particlemass) + einti = einti + particlemass * rad(iradxi)/rhopart + endif etoti = epoti + ekini + einti end subroutine calc_gas_energies From e2631612975f439049a79de31a002830ce3041e1 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Wed, 22 May 2024 19:37:04 +0200 Subject: [PATCH 05/19] (eos_idealplusrad) add functions to calculate egas and erad separately --- src/main/eos_idealplusrad.f90 | 30 ++++++++++++++++++++++++++++-- 1 file changed, 28 insertions(+), 2 deletions(-) diff --git a/src/main/eos_idealplusrad.f90 b/src/main/eos_idealplusrad.f90 index e2a6c10aa..df37e4902 100644 --- a/src/main/eos_idealplusrad.f90 +++ b/src/main/eos_idealplusrad.f90 @@ -24,7 +24,7 @@ module eos_idealplusrad public :: get_idealplusrad_temp,get_idealplusrad_pres,get_idealplusrad_spsoundi,& get_idealgasplusrad_tempfrompres,get_idealplusrad_enfromtemp,& - get_idealplusrad_rhofrompresT + get_idealplusrad_rhofrompresT,egas_from_rhoT,erad_from_rhoT private @@ -125,11 +125,37 @@ subroutine get_idealplusrad_enfromtemp(densi,tempi,mu,eni) real, intent(in) :: densi,tempi,mu real, intent(out) :: eni - eni = 1.5*Rg*tempi/mu + radconst*tempi**4/densi + eni = egas_from_rhoT(densi,tempi,mu) + erad_from_rhoT(densi,tempi,mu) end subroutine get_idealplusrad_enfromtemp +!---------------------------------------------------------------- +!+ +! Calculates specific gas energy from density and temperature +!+ +!---------------------------------------------------------------- +real function egas_from_rhoT(densi,tempi,mu) result(egasi) + real, intent(in) :: densi,tempi,mu + + egasi = 1.5*Rg*tempi/mu + +end function egas_from_rhoT + + +!---------------------------------------------------------------- +!+ +! Calculates specific radiation energy from density and temperature +!+ +!---------------------------------------------------------------- +real function erad_from_rhoT(densi,tempi,mu) result(eradi) + real, intent(in) :: densi,tempi,mu + + eradi = radconst*tempi**4/densi + +end function erad_from_rhoT + + !---------------------------------------------------------------- !+ ! Calculates density from pressure and temperature From d44b98184a3b3d48bf8f0741c71f4d569afbd72a Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Thu, 23 May 2024 15:42:10 +0200 Subject: [PATCH 06/19] big fixes for radiation utility functions --- src/main/eos.f90 | 6 ++++-- src/main/ionization.f90 | 8 +++----- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 67a791772..af2cf2f8a 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -935,11 +935,13 @@ function entropy(rho,pres,mu_in,ientropy,eint_in,ierr,T_in,Trad_in) entropy = kb_on_mh / mu * log(temp**1.5/rho) case(2) ! Include both gas and radiation contributions (up to additive constants) - call get_idealgasplusrad_tempfrompres(pres,rho,mu,temp) ! First solve for temp from rho and pres if (present(Trad_in)) then Trad = Trad_in else - Trad = temp ! assume thermal equilibrium + if (.not. present(T_in)) then + call get_idealgasplusrad_tempfrompres(pres,rho,mu,temp) ! First solve for temp from rho and pres + Trad = temp ! assume thermal equilibrium + endif endif ! check temp if (temp < tiny(0.)) call warning('entropy','temperature = 0 will give minus infinity with s entropy') diff --git a/src/main/ionization.f90 b/src/main/ionization.f90 index d43c2fd10..103d19020 100644 --- a/src/main/ionization.f90 +++ b/src/main/ionization.f90 @@ -348,19 +348,17 @@ subroutine calc_thermal_energy(particlemass,ieos,xyzh,vxyzu,presi,tempi,ethi,rad real, intent(in) :: particlemass,presi,tempi,xyzh(4),vxyzu(4) real, intent(in), optional :: rad(:) real, intent(out) :: ethi - real :: hi,densi_cgs,mui,rhoi + real :: densi_cgs,mui - rhoi = rhoh(hi,particlemass) select case (ieos) case(10,20) ! calculate just gas + radiation thermal energy - hi = xyzh(4) - densi_cgs = rhoi*unit_density + densi_cgs = rhoh(xyzh(4),particlemass)*unit_density mui = densi_cgs * Rg * tempi / (presi*unit_pressure - radconst * tempi**4 / 3.) ! Get mu from pres and temp call get_idealplusrad_enfromtemp(densi_cgs,tempi,mui,ethi) ethi = particlemass * ethi / unit_ergg case default ! assuming internal energy = thermal energy ethi = particlemass * vxyzu(4) - if (do_radiation) ethi = ethi + particlemass*rad(iradxi)/rhoi + if (do_radiation) ethi = ethi + particlemass*rad(iradxi) end select end subroutine calc_thermal_energy From d76ed843d7c4798e9f010e30e5f2b24a2361ec9f Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Thu, 23 May 2024 15:43:00 +0200 Subject: [PATCH 07/19] revert Trad and E conversion functions back to Trad and xi --- src/main/radiation_utils.f90 | 25 +++++++++++++------------ src/setup/set_star_utils.f90 | 4 ++-- 2 files changed, 15 insertions(+), 14 deletions(-) diff --git a/src/main/radiation_utils.f90 b/src/main/radiation_utils.f90 index d54015bd6..19b176db9 100644 --- a/src/main/radiation_utils.f90 +++ b/src/main/radiation_utils.f90 @@ -23,8 +23,8 @@ module radiation_utils public :: get_rad_R public :: radiation_equation_of_state public :: T_from_Etot - public :: radE_from_Trad - public :: Trad_from_radE + public :: radxi_from_Trad + public :: Trad_from_radxi public :: ugas_from_Tgas public :: Tgas_from_ugas public :: get_opacity @@ -136,28 +136,29 @@ end function T_from_Etot !--------------------------------------------------------- !+ -! get the radiation energy from the radiation temperature +! get specific radiation energy from radiation temperature !+ !--------------------------------------------------------- -real function radE_from_Trad(Trad) result(radE) +real function radxi_from_Trad(rho,Trad) result(radxi) use units, only:get_radconst_code - real, intent(in) :: Trad + real, intent(in) :: rho,Trad - radE = Trad**4*get_radconst_code() + radxi = Trad**4*get_radconst_code()/rho -end function radE_from_Trad +end function radxi_from_Trad !--------------------------------------------------------- !+ -! get the radiation temperature from the radiation energy per unit volume!+ +! get radiation temperature from the specific radiation energy +!+ !--------------------------------------------------------- -real function Trad_from_radE(radE) result(Trad) +real function Trad_from_radxi(rho,radxi) result(Trad) use units, only:get_radconst_code - real, intent(in) :: radE + real, intent(in) :: rho,radxi - Trad = (radE/get_radconst_code())**0.25 + Trad = (rho*radxi/get_radconst_code())**0.25 -end function Trad_from_radE +end function Trad_from_radxi !--------------------------------------------------------- !+ diff --git a/src/setup/set_star_utils.f90 b/src/setup/set_star_utils.f90 index 73a5d7017..12942d71a 100644 --- a/src/setup/set_star_utils.f90 +++ b/src/setup/set_star_utils.f90 @@ -382,7 +382,7 @@ subroutine set_star_thermalenergy(ieos,den,pres,r,npts,npart,xyzh,vxyzu,rad,eos_ relaxed,use_var_comp,initialtemp,npin) use part, only:do_radiation,rhoh,massoftype,igas,itemp,igasP,iX,iZ,imu,iradxi use eos, only:equationofstate,calc_temp_and_ene,gamma,gmw - use radiation_utils, only:ugas_from_Tgas,radE_from_Trad + use radiation_utils, only:ugas_from_Tgas,radxi_from_Trad use table_utils, only:yinterp use units, only:unit_density,unit_ergg,unit_pressure integer, intent(in) :: ieos,npart,npts @@ -439,7 +439,7 @@ subroutine set_star_thermalenergy(ieos,den,pres,r,npts,npart,xyzh,vxyzu,rad,eos_ endif if (do_radiation) then vxyzu(4,i) = ugas_from_Tgas(tempi,gamma,gmw) - rad(iradxi,i) = radE_from_Trad(tempi)/densi + rad(iradxi,i) = radxi_from_Trad(tempi) else vxyzu(4,i) = eni / unit_ergg endif From a86125cec4ce6b10a8d30f234a048a4c61c4ba60 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Thu, 23 May 2024 15:45:03 +0200 Subject: [PATCH 08/19] (CE-analysis) calc_gas_energies to output gas thermal and radiation energies --- src/utils/analysis_common_envelope.f90 | 161 +++++++++++++++---------- 1 file changed, 97 insertions(+), 64 deletions(-) diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index 77567ead8..ae2d46ad4 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -527,7 +527,7 @@ subroutine bound_mass(time,npart,particlemass,xyzh,vxyzu) integer, intent(in) :: npart real, intent(in) :: time,particlemass real, intent(inout) :: xyzh(:,:),vxyzu(:,:) - real :: etoti,ekini,epoti,phii,einti,ethi + real :: etoti,ekini,epoti,phii,ereci,egasi,eradi,ethi real :: E_H2,E_HI,E_HeI,E_HeII real, save :: Xfrac,Yfrac,Zfrac real :: rhopart,ponrhoi,spsoundi,tempi,dum1,dum2,dum3 @@ -593,7 +593,8 @@ subroutine bound_mass(time,npart,particlemass,xyzh,vxyzu) do i = 1,npart if (.not. isdead_or_accreted(xyzh(4,i))) then - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,& + egasi,eradi,ereci,etoti) call get_accel_sink_gas(nptmass,xyzh(1,i),xyzh(2,i),xyzh(3,i),xyzh(4,i),xyzmh_ptmass,dum1,dum2,dum3,phii) rhopart = rhoh(xyzh(4,i), particlemass) tempi = eos_vars(itemp,i) @@ -606,7 +607,6 @@ subroutine bound_mass(time,npart,particlemass,xyzh,vxyzu) etoti = 0. epoti = 0. ekini = 0. - einti = 0. ethi = 0. phii = 0. ponrhoi = 0. @@ -686,7 +686,7 @@ subroutine calculate_energies(time,npart,particlemass,xyzh,vxyzu) integer, intent(in) :: npart real, intent(in) :: time,particlemass real, intent(inout) :: xyzh(:,:),vxyzu(:,:) - real :: etoti,ekini,einti,epoti,phii,phii1,jz,fxi,fyi,fzi + real :: etoti,ekini,egasi,eradi,ereci,epoti,phii,phii1,jz,fxi,fyi,fzi real :: rhopart,ponrhoi,spsoundi,tempi,r_ij,radvel real, dimension(3) :: rcrossmv character(len=17), allocatable :: columns(:) @@ -755,7 +755,8 @@ subroutine calculate_energies(time,npart,particlemass,xyzh,vxyzu) jz = rcrossmv(3) encomp(ijz_tot) = encomp(ijz_tot) + jz - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,& + epoti,ekini,egasi,eradi,ereci,etoti) encomp(ipot_ps) = encomp(ipot_ps) + particlemass * phii @@ -927,7 +928,7 @@ subroutine roche_lobe_values(time,npart,particlemass,xyzh,vxyzu) integer, parameter :: iFBV = 20 integer, parameter :: iFBJz = 21 real, dimension(iFBJz) :: MRL - real :: etoti, ekini, einti, epoti, phii, jz + real :: etoti, ekini, ereci, egasi, eradi, epoti, phii, jz logical, dimension(:), allocatable, save:: transferred real, save :: m1, m2 real :: sep, sep1, sep2 @@ -1001,7 +1002,8 @@ subroutine roche_lobe_values(time,npart,particlemass,xyzh,vxyzu) call orbit_com(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptmass,com_xyz,com_vxyz) do i=1,npart - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,& + epoti,ekini,egasi,eradi,ereci,etoti) sep1 = separation(xyzmh_ptmass(1:3,1),xyzh(1:3,i)) sep2 = separation(xyzmh_ptmass(1:3,2),xyzh(1:3,i)) @@ -1132,7 +1134,8 @@ subroutine star_stabilisation_suite(time,npart,particlemass,xyzh,vxyzu) integer, allocatable :: iorder(:),iorder_a(:) real, allocatable :: star_stability(:) real :: total_mass,rhovol,totvol,rhopart,virialpart,virialfluid - real :: phii,ponrhoi,spsoundi,tempi,epoti,ekini,einti,etoti,totekin,totepot,virialintegral,gamma + real :: phii,ponrhoi,spsoundi,tempi,epoti,ekini,egasi,eradi,ereci,etoti + real :: totekin,totepot,virialintegral,gamma integer, parameter :: ivoleqrad = 1 integer, parameter :: idensrad = 2 integer, parameter :: imassout = 3 @@ -1177,7 +1180,8 @@ subroutine star_stabilisation_suite(time,npart,particlemass,xyzh,vxyzu) rhopart = rhoh(xyzh(4,i), particlemass) totvol = totvol + particlemass / rhopart ! Sum "volume" of all particles virialpart = virialpart + particlemass * ( dot_product(fxyzu(1:3,i),xyzh(1:3,i)) + dot_product(vxyzu(1:3,i),vxyzu(1:3,i)) ) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,& + epoti,ekini,egasi,eradi,ereci,etoti) totekin = totekin + ekini totepot = totepot + 0.5*epoti ! Factor of 1/2 to correct for double counting if (rhopart > rho_surface) then @@ -1289,7 +1293,7 @@ subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu) integer :: i,k,Nquantities,ierr,iu integer, save :: quantities_to_calculate(4) integer, allocatable :: iorder(:) - real :: ekini,einti,epoti,ethi,phii,rho_cgs,ponrhoi,spsoundi,tempi,& + real :: ekini,epoti,egasi,eradi,ereci,ethi,phii,rho_cgs,ponrhoi,spsoundi,tempi,& omega_orb,kappai,kappat,kappar,pgas,mu,entropyi,rhopart,& dum1,dum2,dum3,dum4,dum5 real, allocatable, save :: init_entropy(:) @@ -1417,7 +1421,8 @@ subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu) case(1,9) ! Total energy (kin + pot + therm) rhopart = rhoh(xyzh(4,i), particlemass) call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum1) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,& + epoti,ekini,egasi,eradi,ereci,dum1) if (quantities_to_calculate(k)==1) then call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),ethi) quant(k,i) = (ekini + epoti + ethi) / particlemass ! Specific energy @@ -1564,7 +1569,7 @@ end subroutine eos_surfaces subroutine track_particle(time,npart,particlemass,xyzh,vxyzu) use part, only:itemp,iradxi,ilambda use eos, only:entropy - use radiation_utils, only:Trad_from_radE + use radiation_utils, only:Trad_from_radxi use mesa_microphysics, only:getvalue_mesa use ionization_mod, only:ionisation_fraction real, intent(in) :: time,particlemass @@ -1572,7 +1577,7 @@ subroutine track_particle(time,npart,particlemass,xyzh,vxyzu) real, intent(inout) :: xyzh(:,:),vxyzu(:,:) integer, parameter :: nparttotrack=6,ncols=20 real :: r,v,rhopart,ponrhoi,Si,spsoundi,tempi,machi,xh0,xh1,xhe0,xhe1,xhe2,& - ekini,einti,epoti,ethi,eradi,etoti,dum,phii,pgas,mu,rho_cgs,Tradi,lambdai + ekini,egasi,eradi,epoti,ereci,etoti,phii,pgas,mu,rho_cgs,Tradi,lambdai real, dimension(ncols) :: datatable character(len=17) :: filenames(nparttotrack),columns(ncols) integer :: i,k,partID(nparttotrack),ientropy,ierr @@ -1581,7 +1586,7 @@ subroutine track_particle(time,npart,particlemass,xyzh,vxyzu) ! call initial_to_current_IDs(npart,pid_orig) ! partID = (/ 1,2,3,4,5,6,7,8,9,10 /) - partID = (/ 359018, 1669237, 342811, 598910, 1690937, 285745 /) + partID = (/ 193332, 966126, 1303771, 1466288, 142011, 36840 /) columns = (/ ' r',& ' v',& ' rho',& @@ -1593,9 +1598,9 @@ subroutine track_particle(time,npart,particlemass,xyzh,vxyzu) ' mach',& ' ekin',& ' epot',& - ' eth',& + ' egas',& ' erad',& - ' eint',& + ' erec',& ' etot',& ' xHI',& ' xHII',& @@ -1639,19 +1644,17 @@ subroutine track_particle(time,npart,particlemass,xyzh,vxyzu) lambdai = 0. if (do_radiation) then lambdai = radprop(ilambda,i) - Tradi = Trad_from_radE(rad(iradxi,i)) + Tradi = Trad_from_radxi(rhopart,rad(iradxi,i)) Si = entropy(rho_cgs,ponrhoi*rhopart*unit_pressure,mu,ientropy,vxyzu(4,i)*unit_ergg,ierr,Trad_in=Tradi) else Si = entropy(rho_cgs,ponrhoi*rhopart*unit_pressure,mu,ientropy,vxyzu(4,i)*unit_ergg,ierr) endif - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) - call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,tempi,ethi,rad(:,i)) - etoti = ekini + epoti + ethi ! ethi includes radiation energy - eradi = particlemass*rad(iradxi,i)/rhopart + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,& + epoti,ekini,egasi,eradi,ereci,etoti) call ionisation_fraction(rho_cgs,tempi,X_in,1.-X_in-Z_in,xh0,xh1,xhe0,xhe1,xhe2) ! Write file - datatable = (/ r,v,rhopart,tempi,Tradi,lambdai,Si,spsoundi,machi,ekini,epoti,ethi,eradi,einti,etoti,& + datatable = (/ r,v,rhopart,tempi,Tradi,lambdai,Si,spsoundi,machi,ekini,epoti,egasi,eradi,ereci,etoti,& xh0,xh1,xhe0,xhe1,xhe2 /) call write_time_file(trim(adjustl(filenames(k))),columns,time,datatable,ncols,dump_number) enddo @@ -1819,7 +1822,7 @@ subroutine recombination_tau(time,npart,particlemass,xyzh,vxyzu) real, allocatable :: kappa_hist(:),rho_hist(:),tau_r(:),sepbins(:),sepbins_cm(:) logical, allocatable, save :: prev_recombined(:) real :: maxloga,minloga,kappa,kappat,kappar,xh0,xh1,xhe0,xhe1,xhe2,& - ponrhoi,spsoundi,tempi,etoti,ekini,einti,epoti,ethi,phii,dum + ponrhoi,spsoundi,tempi,etoti,ekini,ereci,egasi,eradi,epoti,ethi,phii,dum real, parameter :: recomb_th=0.9 integer :: i,j,nrecombined,bin_ind @@ -1846,7 +1849,8 @@ subroutine recombination_tau(time,npart,particlemass,xyzh,vxyzu) call get_eos_kappa_mesa(rho_part(i)*unit_density,eos_vars(itemp,i),kappa,kappat,kappar) kappa_part(i) = kappa ! In cgs units call ionisation_fraction(rho_part(i)*unit_density,eos_vars(itemp,i),X_in,1.-X_in-Z_in,xh0,xh1,xhe0,xhe1,xhe2) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) ! Calculate total energy + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,& + epoti,ekini,egasi,eradi,ereci,dum) call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rho_part(i),eos_vars(itemp,i),ethi) etoti = ekini + epoti + ethi if ((xh0 > recomb_th) .and. (.not. prev_recombined(i)) .and. (etoti < 0.)) then ! Recombination event and particle is still bound @@ -1906,7 +1910,7 @@ subroutine energy_hist(time,npart,particlemass,xyzh,vxyzu) character(len=40) :: data_formatter integer :: nbins,nhists,i,unitnum real, allocatable :: hist(:),coord(:,:),Emin(:),Emax(:) - real :: rhopart,ponrhoi,spsoundi,tempi,phii,epoti,ekini,einti,ethi,dum + real :: rhopart,ponrhoi,spsoundi,tempi,phii,epoti,ekini,ereci,egasi,eradi,ethi,dum real, allocatable :: quant(:) logical :: ilogbins @@ -1925,11 +1929,12 @@ subroutine energy_hist(time,npart,particlemass,xyzh,vxyzu) do i=1,npart rhopart = rhoh(xyzh(4,i), particlemass) call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,& + epoti,ekini,egasi,eradi,ereci,dum) if (ieos==10 .or. ieos==20) then call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),ethi) else - ethi = einti + ethi = ethi+ereci endif coord(i,1) = (ekini + epoti)/particlemass coord(i,2) = vxyzu(4,i) - ethi/particlemass @@ -1970,7 +1975,7 @@ subroutine energy_profile(time,npart,particlemass,xyzh,vxyzu) integer :: nbins real, allocatable :: coord(:) real, allocatable :: hist(:),quant(:,:) - real :: ekini,einti,epoti,ethi,phii,pgas,mu,dum,rhopart,ponrhoi,spsoundi,tempi,& + real :: ekini,ereci,egasi,eradi,epoti,ethi,phii,pgas,mu,dum,rhopart,ponrhoi,spsoundi,tempi,& maxcoord,mincoord,xh0,xh1,xhe0,xhe1,xhe2 character(len=17), allocatable :: filename(:),headerline(:) character(len=40) :: data_formatter @@ -2072,7 +2077,8 @@ subroutine energy_profile(time,npart,particlemass,xyzh,vxyzu) call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) select case (iquantity) case(1) ! Energy - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,& + epoti,ekini,egasi,eradi,ereci,dum) call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),ethi) quant(i,1) = ekini + epoti + ethi case(2) ! Entropy @@ -2088,7 +2094,8 @@ subroutine energy_profile(time,npart,particlemass,xyzh,vxyzu) quant(i,1) = entropy(rhopart*unit_density,ponrhoi*rhopart*unit_pressure,mu,ientropy,ierr=ierr) endif case(3) ! Bernoulli energy (per unit mass) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,& + epoti,ekini,egasi,eradi,ereci,dum) quant(i,1) = 0.5*dot_product(vxyzu(1:3,i),vxyzu(1:3,i)) + ponrhoi + vxyzu(4,i) + epoti/particlemass ! 1/2 v^2 + P/rho + phi case(4) ! Ion fraction call ionisation_fraction(rhopart*unit_density,eos_vars(itemp,i),X_in,1.-X_in-Z_in,xh0,xh1,xhe0,xhe1,xhe2) @@ -2213,14 +2220,15 @@ subroutine velocity_histogram(time,num,npart,particlemass,xyzh,vxyzu) character(len=40) :: data_formatter character(len=40) :: file_name1,file_name2 integer :: i,iu1,iu2,ncols - real :: ponrhoi,rhopart,spsoundi,phii,epoti,ekini,einti,tempi,ethi,dum + real :: ponrhoi,rhopart,spsoundi,phii,epoti,ekini,ereci,egasi,eradi,tempi,ethi,dum real, allocatable :: vbound(:),vunbound(:),vr(:) allocate(vbound(npart),vunbound(npart),vr(npart)) do i = 1,npart rhopart = rhoh(xyzh(4,i), particlemass) call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,& + epoti,ekini,egasi,eradi,ereci,dum) call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),ethi) vr(i) = dot_product(xyzh(1:3,i),vxyzu(1:3,i)) / sqrt(dot_product(xyzh(1:3,i),xyzh(1:3,i))) @@ -2495,7 +2503,7 @@ subroutine unbound_profiles(time,num,npart,particlemass,xyzh,vxyzu) integer, dimension(4) :: npart_hist real, dimension(5,npart) :: dist_part,rad_part real, dimension(:), allocatable :: hist_var - real :: etoti,ekini,einti,epoti,ethi,phii,dum,rhopart,ponrhoi,spsoundi,tempi + real :: etoti,ekini,ereci,egasi,eradi,epoti,ethi,phii,dum,rhopart,ponrhoi,spsoundi,tempi real :: maxloga,minloga character(len=18), dimension(4) :: grid_file character(len=40) :: data_formatter @@ -2528,7 +2536,8 @@ subroutine unbound_profiles(time,num,npart,particlemass,xyzh,vxyzu) if (.not. isdead_or_accreted(xyzh(4,i))) then rhopart = rhoh(xyzh(4,i), particlemass) call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,& + epoti,ekini,egasi,eradi,ereci,dum) call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,tempi,ethi) etoti = ekini + epoti + ethi @@ -2613,7 +2622,7 @@ subroutine unbound_ionfrac(time,npart,particlemass,xyzh,vxyzu) real, intent(inout) :: xyzh(:,:),vxyzu(:,:) character(len=17) :: columns(5) integer :: i - real :: etoti,ekini,einti,epoti,ethi,phii,dum,rhopart,xion(1:4),& + real :: etoti,ekini,egasi,eradi,ereci,epoti,ethi,phii,dum,rhopart,xion(1:4),& ponrhoi,spsoundi,tempi,xh0,xh1,xhe0,xhe1,xhe2 logical, allocatable, save :: prev_unbound(:),prev_bound(:) real, allocatable, save :: ionfrac(:,:) @@ -2636,7 +2645,7 @@ subroutine unbound_ionfrac(time,npart,particlemass,xyzh,vxyzu) do i=1,npart rhopart = rhoh(xyzh(4,i), particlemass) call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,egasi,eradi,ereci,dum) call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,tempi,ethi) etoti = ekini + epoti + ethi @@ -2689,7 +2698,7 @@ subroutine unbound_temp(time,npart,particlemass,xyzh,vxyzu) real, intent(inout) :: xyzh(:,:),vxyzu(:,:) character(len=17) :: columns(1) integer :: i,final_count(7) - real :: etoti,ekini,einti,epoti,ethi,phii,dum,rhopart,& + real :: etoti,ekini,ereci,egasi,eradi,epoti,ethi,phii,dum,rhopart,& ponrhoi,spsoundi,temp_bins(7) logical, allocatable, save :: prev_unbound(:),prev_bound(:) real, allocatable, save :: temp_unbound(:) @@ -2706,7 +2715,7 @@ subroutine unbound_temp(time,npart,particlemass,xyzh,vxyzu) do i=1,npart rhopart = rhoh(xyzh(4,i), particlemass) call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),eos_vars(itemp,i),vxyzu(4,i)) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,egasi,eradi,ereci,dum) call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),ethi) etoti = ekini + epoti + ethi @@ -2761,7 +2770,7 @@ subroutine recombination_stats(time,num,npart,particlemass,xyzh,vxyzu) integer, intent(in) :: npart,num real, intent(in) :: time,particlemass real, intent(inout) :: xyzh(:,:),vxyzu(:,:) - real :: etoti,ekini,einti,epoti,ethi,phii,dum,rhopart,& + real :: etoti,ekini,egasi,eradi,ereci,epoti,ethi,phii,dum,rhopart,& ponrhoi,spsoundi,tempi,pressure,temperature,xh0,xh1,xhe0,xhe1,xhe2 character(len=40) :: data_formatter,logical_format logical, allocatable :: isbound(:) @@ -2776,7 +2785,7 @@ subroutine recombination_stats(time,num,npart,particlemass,xyzh,vxyzu) ! Calculate total energy rhopart = rhoh(xyzh(4,i), particlemass) call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,egasi,eradi,ereci,dum) call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),ethi) etoti = ekini + epoti + ethi @@ -2999,7 +3008,7 @@ subroutine bound_unbound_thermo(time,npart,particlemass,xyzh,vxyzu) character(len=17), allocatable :: columns(:) integer :: i, ncols real, dimension(8) :: entropy_array - real :: etoti, ekini, einti, epoti, phii, rhopart + real :: etoti, ekini, epoti, phii, rhopart,egasi,eradi,ereci real :: pres_1, proint_1, peint_1, temp_1 real :: troint_1, teint_1, entrop_1, abad_1, gamma1_1, gam_1 integer, parameter :: ient_b = 1 @@ -3022,7 +3031,8 @@ subroutine bound_unbound_thermo(time,npart,particlemass,xyzh,vxyzu) call compute_energies(time) do i=1,npart - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,& + epoti,ekini,egasi,eradi,ereci,etoti) rhopart = rhoh(xyzh(4,i), particlemass) @@ -3034,7 +3044,7 @@ subroutine bound_unbound_thermo(time,npart,particlemass,xyzh,vxyzu) !sums entropy and other quantities for bound particles and unbound particles if (.not. switch(1)) then - etoti = etoti - einti + etoti = etoti - egasi - eradi - ereci endif if (etoti < 0.0) then !bound @@ -3352,7 +3362,7 @@ subroutine J_E_plane(num,npart,particlemass,xyzh,vxyzu) real, intent(in) :: particlemass,xyzh(:,:),vxyzu(:,:) character(len=17), allocatable :: columns(:) integer :: ncols,i - real :: com_xyz(3),com_vxyz(3),dum1,dum2,dum3,dum4,etoti,angmom_com(3),angmom_core(3) + real :: com_xyz(3),com_vxyz(3),dum1,dum2,dum3,dum4,dum5,dum6,etoti,angmom_com(3),angmom_core(3) real, allocatable :: data(:,:) ncols = 7 @@ -3368,7 +3378,7 @@ subroutine J_E_plane(num,npart,particlemass,xyzh,vxyzu) call get_centreofmass(com_xyz,com_vxyz,npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptmass) do i=1,npart - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,dum1,dum2,dum3,dum4,etoti) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,dum1,dum2,dum3,dum4,dum5,dum6,etoti) data(1,i) = etoti call cross_product3D(xyzh(1:3,i)-xyzmh_ptmass(1:3,1), vxyzu(1:3,i)-vxyz_ptmass(1:3,1), angmom_core) data(5:7,i) = angmom_core @@ -3769,33 +3779,58 @@ end subroutine get_gas_omega !---------------------------------------------------------------- !+ ! Calculate kinetic, gravitational potential (gas-gas and sink-gas), -! and internal energy of a gas particle. +! and other energies of a gas particle. !+ !---------------------------------------------------------------- -subroutine calc_gas_energies(particlemass,poten,xyzh,vxyzu,rad,xyzmh_ptmass,phii,epoti,ekini,einti,etoti) +subroutine calc_gas_energies(particlemass,poten,xyzh,vxyzu,rad,xyzmh_ptmass,phii,epoti,ekini,egasi,eradi,ereci,etoti) ! Warning: Do not sum epoti or etoti as it is to obtain a total energy; this would not give the correct ! total energy due to complications related to double-counting. - use ptmass, only:get_accel_sink_gas - use part, only:nptmass,iradxi + use ptmass, only:get_accel_sink_gas + use part, only:nptmass,iradxi,itemp + use eos_idealplusrad, only:get_idealplusrad_temp,egas_from_rhoT,erad_from_rhoT real, intent(in) :: particlemass real(4), intent(in) :: poten real, intent(in) :: xyzh(:),vxyzu(:),rad(:) real, dimension(5,nptmass), intent(in) :: xyzmh_ptmass - real, intent(out) :: phii,epoti,ekini,einti,etoti - real :: fxi,fyi,fzi,rhopart - - phii = 0.0 + real, intent(out) :: phii,epoti,ekini,egasi,eradi,ereci,etoti + real :: fxi,fyi,fzi,rhoi,spsoundi,ponrhoi,presi,tempi,egasradi + integer :: ierr + rhoi = rhoh(xyzh(4),particlemass) + phii = 0. call get_accel_sink_gas(nptmass,xyzh(1),xyzh(2),xyzh(3),xyzh(4),xyzmh_ptmass,fxi,fyi,fzi,phii) - epoti = 2.*poten + particlemass * phii ! For individual particles, need to multiply 2 to poten to get \sum_j G*mi*mj/r ekini = particlemass * 0.5 * dot_product(vxyzu(1:3),vxyzu(1:3)) - einti = particlemass * vxyzu(4) + + egasi = 0. + ereci = 0. if (do_radiation) then - rhopart = rhoh(xyzh(4),particlemass) - einti = einti + particlemass * rad(iradxi)/rhopart + eradi = rad(iradxi)*particlemass + else + eradi = 0. endif - etoti = epoti + ekini + einti + + select case (ieos) + case(2) + egasi = vxyzu(4)*particlemass + egasradi = egasi + eradi + case(10) ! not tested + eradi = 0. ! not implemented + egasi = 0. ! not implemented + call equationofstate(ieos,ponrhoi,spsoundi,rhoi,xyzh(1),xyzh(2),xyzh(3),tempi,vxyzu(4)) + presi = ponrhoi*rhoi + call calc_thermal_energy(particlemass,10,xyzh,vxyzu,presi,tempi,egasradi,rad) + ereci = vxyzu(4)*particlemass - egasradi + case(12) + call get_idealplusrad_temp(rhoi,vxyzu(4)*unit_ergg,gmw,tempi,ierr) + egasi = egas_from_rhoT(rhoi,tempi,gmw)*particlemass + eradi = erad_from_rhoT(rhoi,tempi,gmw)*particlemass + egasradi = egasi + eradi + case default + call fatal('calc_gas_energies',"EOS type not supported (currently, only supporting ieos=2,10,12)") + end select + + etoti = epoti + ekini + ereci + egasradi end subroutine calc_gas_energies @@ -3869,7 +3904,7 @@ subroutine stellar_profile(time,ncols,particlemass,npart,xyzh,vxyzu,profile,simp real :: proj(3),orth(3),proj_mag,orth_dist,orth_ratio real :: rhopart,ponrhoi,spsoundi,tempi real :: temp,kappa,kappat,kappar,pres - real :: ekini,epoti,einti,etoti,phii + real :: ekini,epoti,egasi,eradi,ereci,etoti,phii real :: xh0, xh1, xhe0, xhe1, xhe2 real :: temp_profile(ncols,npart) logical :: criteria @@ -3922,7 +3957,7 @@ subroutine stellar_profile(time,ncols,particlemass,npart,xyzh,vxyzu,profile,simp endif call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),& - xyzmh_ptmass,phii,epoti,ekini,einti,etoti) + xyzmh_ptmass,phii,epoti,ekini,egasi,eradi,ereci,etoti) call ionisation_fraction(rhopart*unit_density,temp,X_in,1.-X_in-Z_in,xh0,xh1,xhe0,xhe1,xhe2) @@ -4466,10 +4501,8 @@ subroutine set_eos_options(analysis_to_perform) case(2,12) gamma = 5./3. call prompt('Enter gamma:',gamma,0.) - if (ieos==12) then - gmw = 0.618212823 - call prompt('Enter mean molecular weight for gas+rad EoS:',gmw,0.) - endif + gmw = 0.618212823 + call prompt('Enter mean molecular weight for gas+rad EoS:',gmw,0.) case(10,20) gamma = 5./3. X_in = 0.69843 From cfc2e49dc0f5e450df1c92754e42bda9eeaa6a91 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Fri, 7 Jun 2024 17:23:51 +0200 Subject: [PATCH 09/19] (CE-analysis) clean up unbound_profiles subroutine --- src/utils/analysis_common_envelope.f90 | 35 +++++++++++--------------- 1 file changed, 14 insertions(+), 21 deletions(-) diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index ae2d46ad4..cbaa2cf73 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -2497,18 +2497,17 @@ end subroutine planet_profile !+ !---------------------------------------------------------------- subroutine unbound_profiles(time,num,npart,particlemass,xyzh,vxyzu) - integer, intent(in) :: npart,num - real, intent(in) :: time,particlemass - real, intent(inout) :: xyzh(:,:),vxyzu(:,:) - integer, dimension(4) :: npart_hist - real, dimension(5,npart) :: dist_part,rad_part - real, dimension(:), allocatable :: hist_var - real :: etoti,ekini,ereci,egasi,eradi,epoti,ethi,phii,dum,rhopart,ponrhoi,spsoundi,tempi - real :: maxloga,minloga - character(len=18), dimension(4) :: grid_file - character(len=40) :: data_formatter - logical, allocatable, save :: prev_unbound(:,:),prev_bound(:,:) - integer :: i,unitnum,nbins + integer, intent(in) :: npart,num + real, intent(in) :: time,particlemass + real, intent(inout) :: xyzh(:,:),vxyzu(:,:) + integer, dimension(4) :: npart_hist + real, dimension(5,npart) :: dist_part,rad_part + real, dimension(:), allocatable :: hist_var + real :: etoti,ekini,ereci,egasi,eradi,epoti,phii,dum,maxloga,minloga + character(len=18), dimension(4) :: grid_file + character(len=40) :: data_formatter + logical, allocatable, save :: prev_unbound(:,:),prev_bound(:,:) + integer :: i,unitnum,nbins call compute_energies(time) npart_hist = 0 ! Stores number of particles fulfilling each of the four bound/unbound criterion @@ -2534,12 +2533,9 @@ subroutine unbound_profiles(time,num,npart,particlemass,xyzh,vxyzu) do i=1,npart if (.not. isdead_or_accreted(xyzh(4,i))) then - rhopart = rhoh(xyzh(4,i), particlemass) - call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,& epoti,ekini,egasi,eradi,ereci,dum) - call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,tempi,ethi) - etoti = ekini + epoti + ethi + etoti = ekini + epoti + egasi + eradi ! Ekin + Epot + Eth > 0 if ((etoti > 0.) .and. (.not. prev_unbound(1,i))) then @@ -2590,15 +2586,12 @@ subroutine unbound_profiles(time,num,npart,particlemass,xyzh,vxyzu) write(data_formatter, "(a,I5,a)") "(", nbins+1, "(3x,es18.10e3,1x))" ! Time column plus nbins columns if (num == 0) then ! Write header line - unitnum = 1000 - open(unit=unitnum,file=trim(adjustl(grid_file(i))),status='replace') + open(newunit=unitnum,file=trim(adjustl(grid_file(i))),status='replace') write(unitnum, "(a)") '# Newly bound/unbound particles' close(unit=unitnum) endif - unitnum=1001+i - - open(unit=unitnum,file=trim(adjustl(grid_file(i))), position='append') + open(newunit=unitnum,file=trim(adjustl(grid_file(i))), position='append') write(unitnum,"()") write(unitnum,data_formatter) time,hist_var(:) From 78cc02d6b80938ec6b5a747983b769a0bcad5881 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Mon, 8 Jul 2024 16:24:19 +0200 Subject: [PATCH 10/19] (CE-analysis) bug fix: feed rho in phys units to thermal energy calculation --- src/main/eos_idealplusrad.f90 | 6 +++--- src/utils/analysis_common_envelope.f90 | 9 +++++---- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/src/main/eos_idealplusrad.f90 b/src/main/eos_idealplusrad.f90 index df37e4902..72ecc22db 100644 --- a/src/main/eos_idealplusrad.f90 +++ b/src/main/eos_idealplusrad.f90 @@ -125,7 +125,7 @@ subroutine get_idealplusrad_enfromtemp(densi,tempi,mu,eni) real, intent(in) :: densi,tempi,mu real, intent(out) :: eni - eni = egas_from_rhoT(densi,tempi,mu) + erad_from_rhoT(densi,tempi,mu) + eni = egas_from_rhoT(tempi,mu) + erad_from_rhoT(densi,tempi,mu) end subroutine get_idealplusrad_enfromtemp @@ -135,8 +135,8 @@ end subroutine get_idealplusrad_enfromtemp ! Calculates specific gas energy from density and temperature !+ !---------------------------------------------------------------- -real function egas_from_rhoT(densi,tempi,mu) result(egasi) - real, intent(in) :: densi,tempi,mu +real function egas_from_rhoT(tempi,mu) result(egasi) + real, intent(in) :: tempi,mu egasi = 1.5*Rg*tempi/mu diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index cbaa2cf73..07ac58cf6 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -3786,10 +3786,11 @@ subroutine calc_gas_energies(particlemass,poten,xyzh,vxyzu,rad,xyzmh_ptmass,phii real, intent(in) :: xyzh(:),vxyzu(:),rad(:) real, dimension(5,nptmass), intent(in) :: xyzmh_ptmass real, intent(out) :: phii,epoti,ekini,egasi,eradi,ereci,etoti - real :: fxi,fyi,fzi,rhoi,spsoundi,ponrhoi,presi,tempi,egasradi + real :: fxi,fyi,fzi,rhoi,rho_cgs,spsoundi,ponrhoi,presi,tempi,egasradi integer :: ierr rhoi = rhoh(xyzh(4),particlemass) + rho_cgs = rhoi*unit_density phii = 0. call get_accel_sink_gas(nptmass,xyzh(1),xyzh(2),xyzh(3),xyzh(4),xyzmh_ptmass,fxi,fyi,fzi,phii) epoti = 2.*poten + particlemass * phii ! For individual particles, need to multiply 2 to poten to get \sum_j G*mi*mj/r @@ -3815,9 +3816,9 @@ subroutine calc_gas_energies(particlemass,poten,xyzh,vxyzu,rad,xyzmh_ptmass,phii call calc_thermal_energy(particlemass,10,xyzh,vxyzu,presi,tempi,egasradi,rad) ereci = vxyzu(4)*particlemass - egasradi case(12) - call get_idealplusrad_temp(rhoi,vxyzu(4)*unit_ergg,gmw,tempi,ierr) - egasi = egas_from_rhoT(rhoi,tempi,gmw)*particlemass - eradi = erad_from_rhoT(rhoi,tempi,gmw)*particlemass + call get_idealplusrad_temp(rho_cgs,vxyzu(4)*unit_ergg,gmw,tempi,ierr) + egasi = egas_from_rhoT(tempi,gmw)/unit_ergg*particlemass + eradi = erad_from_rhoT(rho_cgs,tempi,gmw)/unit_ergg*particlemass egasradi = egasi + eradi case default call fatal('calc_gas_energies',"EOS type not supported (currently, only supporting ieos=2,10,12)") From 2edd6e8449f6dceafcfb6abf05f66bdf82ba7588 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Wed, 4 Sep 2024 14:53:12 +0200 Subject: [PATCH 11/19] (CE-analysis) change writing of divv files to splash extra-column files --- src/main/inject_rochelobe.f90 | 1 - src/utils/analysis_common_envelope.f90 | 324 ++++++++++++------------- 2 files changed, 161 insertions(+), 164 deletions(-) diff --git a/src/main/inject_rochelobe.f90 b/src/main/inject_rochelobe.f90 index ffe84d4dd..5aa329cd2 100644 --- a/src/main/inject_rochelobe.f90 +++ b/src/main/inject_rochelobe.f90 @@ -165,7 +165,6 @@ subroutine inject_particles(time,dtlast,xyzh,vxyzu,xyzmh_ptmass,vxyz_ptmass,& h = hfact*sw_chi/udist !add the particle call add_or_update_particle(part_type, xyzi, vxyz, h, u, i_part, npart, npartoftype, xyzh, vxyzu) - enddo ! !-- no constraint on timestep diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index cbaa2cf73..3f2140ce3 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -71,7 +71,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) ' 5) Roche-lobe utils', & ' 6) Star stabilisation suite', & ' 7) Simulation units and particle properties', & - ' 8) Output .divv', & + ' 8) Output extra quantities', & ' 9) EoS testing', & '10) Profile of newly unbound particles', & '11) Sink properties', & @@ -129,8 +129,8 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) call star_stabilisation_suite(time,npart,particlemass,xyzh,vxyzu) case(7) !Units call print_simulation_parameters(npart,particlemass) - case(8) !Output .divv - call output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu) + case(8) ! output extra quantities + call output_extra_quantities(time,dumpfile,npart,particlemass,xyzh,vxyzu) case(9) !EoS testing call eos_surfaces case(10) !New unbound particle profiles in time @@ -1274,110 +1274,183 @@ end subroutine print_simulation_parameters !---------------------------------------------------------------- !+ -! Write quantities (up to four) to divv file +! Write extra quantities to .extras files +! These files can be read via splash --extracols=[label1],[label2],... !+ !---------------------------------------------------------------- -subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu) +subroutine output_extra_quantities(time,dumpfile,npart,particlemass,xyzh,vxyzu) use part, only:eos_vars,itemp,nucleation,idK0,idK1,idK2,idK3,idJstar,idmu,idgamma use eos, only:entropy use eos_mesa, only:get_eos_kappa_mesa use mesa_microphysics, only:getvalue_mesa use sortutils, only:set_r2func_origin,r2func_origin,indexxfunc use ionization_mod, only:ionisation_fraction - use dust_formation, only:psat_C,eps,set_abundances,mass_per_H, chemical_equilibrium_light, calc_nucleation!, Scrit - !use dim, only:nElements + use dust_formation, only:psat_C,eps,set_abundances,mass_per_H,chemical_equilibrium_light,calc_nucleation integer, intent(in) :: npart character(len=*), intent(in) :: dumpfile real, intent(in) :: time,particlemass real, intent(inout) :: xyzh(:,:),vxyzu(:,:) - integer :: i,k,Nquantities,ierr,iu - integer, save :: quantities_to_calculate(4) - integer, allocatable :: iorder(:) + character(len=30) :: msg + character(len=17), allocatable:: labels(:) + integer :: i,k,Noptions,ierr + integer, save :: Nquant + integer, save, allocatable :: quants(:) + integer, allocatable :: iorder(:),iu(:) real :: ekini,epoti,egasi,eradi,ereci,ethi,phii,rho_cgs,ponrhoi,spsoundi,tempi,& - omega_orb,kappai,kappat,kappar,pgas,mu,entropyi,rhopart,& + omega_orb,kappai,kappat,kappar,pgas,mu,entropyi,rhopart,v_esci,& dum1,dum2,dum3,dum4,dum5 + real :: pC,pC2,pC2H,pC2H2,nH_tot,epsC,S,taustar,taugr,JstarS real, allocatable, save :: init_entropy(:) - real, allocatable :: quant(:,:) - real, dimension(3) :: com_xyz,com_vxyz,xyz_a,vxyz_a - real :: pC, pC2, pC2H, pC2H2, nH_tot, epsC, S - real :: taustar, taugr, JstarS - real :: v_esci + real, allocatable :: arr(:,:) + real, dimension(3) :: com_xyz,com_vxyz,xyz_a,vxyz_a,sinkcom_xyz,sinkcom_vxyz real, parameter :: Scrit = 2. ! Critical saturation ratio - logical :: verbose = .false. + logical :: req_eos_call,req_gas_energy,req_thermal_energy,verbose=.false. + + Noptions = 13 + allocate(labels(Noptions)) + labels = (/ 'e_kpt ',& + 'e_kp ',& + 'erec ',& + 'mach ',& + 'kappa ',& + 'omega_sinkCM',& + 'omega_core ',& + 'delta_omega ',& + 'entropy_mesa',& + 'entropy_gain',& + 'm ',& + 'vesc ',& + 'JstarS '& + /) - allocate(quant(4,npart)) - Nquantities = 14 if (dump_number == 0) then - print "(14(a,/))",& + call prompt('Enter number of extra quantities to write out: ',Nquant,0) + allocate(quants(Nquant)) + + print "(13(a,/))",& '1) Total energy (kin + pot + therm)', & - '2) Mach number', & - '3) Opacity from MESA tables', & - '4) Gas omega w.r.t. effective CoM', & - '5) Fractional difference between gas and orbital omega', & - '6) MESA EoS specific entropy', & - '7) Fractional entropy gain', & - '8) Specific recombination energy', & - '9) Total energy (kin + pot)', & - '10) Mass coordinate', & - '11) Gas omega w.r.t. CoM', & - '12) Gas omega w.r.t. sink 1',& - '13) JstarS', & - '14) Escape velocity' - - quantities_to_calculate = (/1,2,4,5/) - call prompt('Choose first quantity to compute ',quantities_to_calculate(1),0,Nquantities) - call prompt('Choose second quantity to compute ',quantities_to_calculate(2),0,Nquantities) - call prompt('Choose third quantity to compute ',quantities_to_calculate(3),0,Nquantities) - call prompt('Choose fourth quantity to compute ',quantities_to_calculate(4),0,Nquantities) + '2) Total energy (kin + pot)', & + '3) Specific recombination energy', & + '4) Mach number', & + '5) Opacity from MESA tables', & + '6) Gas omega w.r.t. sink CoM', & + '7) Gas omega w.r.t. sink 1',& + '8) Fractional difference between gas and orbital omega w.r.t. sink CoM', & + '9) MESA EoS specific entropy', & + '10) Fractional entropy gain', & + '11) Mass coordinate', & + '12) Escape velocity', & + '13) JstarS' + + do i=1,Nquant + write(msg, '(a,i2,a)') 'Enter quantity ',i,':' + call prompt(msg,quants(i),0,Noptions) + enddo endif + allocate(arr(Nquant,npart),iu(Nquant)) + arr = 0. + ! Calculations performed outside loop over particles call compute_energies(time) omega_orb = 0. com_xyz = 0. com_vxyz = 0. - do k=1,4 - select case (quantities_to_calculate(k)) - case(0,1,2,3,6,8,9,13,14) ! Nothing to do - case(4,5,11,12) ! Fractional difference between gas and orbital omega - if (quantities_to_calculate(k) == 4 .or. quantities_to_calculate(k) == 5) then - com_xyz = (xyzmh_ptmass(1:3,1)*xyzmh_ptmass(4,1) + xyzmh_ptmass(1:3,2)*xyzmh_ptmass(4,2)) & - / (xyzmh_ptmass(4,1) + xyzmh_ptmass(4,2)) - com_vxyz = (vxyz_ptmass(1:3,1)*xyzmh_ptmass(4,1) + vxyz_ptmass(1:3,2)*xyzmh_ptmass(4,2)) & - / (xyzmh_ptmass(4,1) + xyzmh_ptmass(4,2)) - elseif (quantities_to_calculate(k) == 11 .or. quantities_to_calculate(k) == 12) then - com_xyz = xyzmh_ptmass(1:3,1) - com_vxyz = vxyz_ptmass(1:3,1) - endif - do i=1,nptmass - xyz_a(1:3) = xyzmh_ptmass(1:3,i) - com_xyz(1:3) - vxyz_a(1:3) = vxyz_ptmass(1:3,i) - com_vxyz(1:3) - omega_orb = omega_orb + 0.5 * (-xyz_a(2) * vxyz_a(1) + xyz_a(1) * vxyz_a(2)) / dot_product(xyz_a(1:2), xyz_a(1:2)) - enddo - case(7) - if (dump_number==0) allocate(init_entropy(npart)) - case(10) - call set_r2func_origin(0.,0.,0.) - allocate(iorder(npart)) - call indexxfunc(npart,r2func_origin,xyzh,iorder) - deallocate(iorder) - case default - print*,"Error: Requested quantity is invalid." - stop - end select - enddo + epoti = 0. + ekini = 0. + + req_eos_call = any(quants==1 .or. quants==2 .or. quants==4 .or. quants==6 .or. quants==7 & + .or. quants==9 .or. quants==10 .or. quants==13) + req_gas_energy = any(quants==1 .or. quants==2 .or. quants==3) + req_thermal_energy = any(quants==1 .or. quants==3) + + if (any(quants==6 .or. quants==8)) then + sinkcom_xyz = (xyzmh_ptmass(1:3,1)*xyzmh_ptmass(4,1) + xyzmh_ptmass(1:3,2)*xyzmh_ptmass(4,2)) & + / (xyzmh_ptmass(4,1) + xyzmh_ptmass(4,2)) + sinkcom_vxyz = (vxyz_ptmass(1:3,1)*xyzmh_ptmass(4,1) + vxyz_ptmass(1:3,2)*xyzmh_ptmass(4,2)) & + / (xyzmh_ptmass(4,1) + xyzmh_ptmass(4,2)) + endif + + if (any(quants==11)) then + call set_r2func_origin(0.,0.,0.) + allocate(iorder(npart)) + call indexxfunc(npart,r2func_origin,xyzh,iorder) + endif + + if (any(quants==10) .and. dump_number==0) allocate(init_entropy(npart)) + + if (any(quants==13)) call set_abundances ! set initial abundances to get mass_per_H + - !set initial abundances to get mass_per_H - call set_abundances - ! Calculations performed in loop over particles do i=1,npart - do k=1,4 - select case (quantities_to_calculate(k)) + rhopart = rhoh(xyzh(4,i),particlemass) + rho_cgs = rhopart*unit_density + if (req_eos_call) then + call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) + endif + + if (req_gas_energy) then + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,& + epoti,ekini,egasi,eradi,ereci,dum1) + endif + + if (req_thermal_energy) then + call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),ethi) + endif + + do k=1,Nquant + select case (quants(k)) + case(1) ! Total energy (kin + pot + therm) + arr(k,i) = (ekini + epoti + ethi) / particlemass + case(2) ! Total energy (kin + pot) + arr(k,i) = (ekini + epoti) / particlemass + case(3) ! Specific recombination energy + arr(k,i) = vxyzu(4,i) - ethi / particlemass + case(4) ! Mach number + arr(k,i) = distance(vxyzu(1:3,i)) / spsoundi + case(5) ! Opacity from MESA tables + call ionisation_fraction(rho_cgs,eos_vars(itemp,i),X_in,1.-X_in-Z_in,dum1,dum2,dum3,dum4,dum5) + if (ieos == 10) then + call get_eos_kappa_mesa(rho_cgs,eos_vars(itemp,i),kappai,kappat,kappar) + arr(k,i) = kappai + else + arr(k,i) = 0. + endif + case(6) ! Gas omega w.r.t. sink CoM + xyz_a = xyzh(1:3,i) - sinkcom_xyz(1:3) + vxyz_a = vxyzu(1:3,i) - sinkcom_vxyz(1:3) + arr(k,i) = (-xyz_a(2) * vxyz_a(1) + xyz_a(1) * vxyz_a(2)) / dot_product(xyz_a(1:2), xyz_a(1:2)) + case(7) ! Gas omega w.r.t. sink 1 + xyz_a = xyzh(1:3,i) - xyzmh_ptmass(1:3,1) + vxyz_a = vxyzu(1:3,i) - vxyz_ptmass(1:3,1) + arr(k,i) = (-xyz_a(2) * vxyz_a(1) + xyz_a(1) * vxyz_a(2)) / dot_product(xyz_a(1:2), xyz_a(1:2)) + case(8) ! Fractional difference between gas and orbital omega + xyz_a = xyzh(1:3,i) - sinkcom_xyz(1:3) + vxyz_a = vxyzu(1:3,i) - sinkcom_vxyz(1:3) + omega_orb = omega_orb + 0.5 * (-xyz_a(2) * vxyz_a(1) + xyz_a(1) * vxyz_a(2)) / dot_product(xyz_a(1:2), xyz_a(1:2)) + arr(k,i) = (-xyz_a(2) * vxyz_a(1) + xyz_a(1) * vxyz_a(2)) / dot_product(xyz_a(1:2), xyz_a(1:2)) + arr(k,i) = arr(k,i)/omega_orb - 1. + case(9,10) ! Calculate MESA EoS entropy + entropyi = 0. + if (ieos==10) then + call getvalue_mesa(rho_cgs,vxyzu(4,i)*unit_ergg,3,pgas,ierr) ! Get gas pressure + mu = rho_cgs * Rg * eos_vars(itemp,i) / pgas + entropyi = entropy(rho_cgs,ponrhoi*rhopart*unit_pressure,mu,3,vxyzu(4,i)*unit_ergg,ierr) + elseif (ieos==2) then + entropyi = entropy(rho_cgs,ponrhoi*rhopart*unit_pressure,gmw,1) + endif + if (quants(k) == 10) then + if (dump_number == 0) init_entropy(i) = entropyi ! Store initial entropy on each particle + arr(k,i) = entropyi/init_entropy(i) - 1. + elseif (quants(k) == 9) then + arr(k,i) = entropyi + endif + case(11) ! Mass coordinate + arr(k,iorder(i)) = real(i,kind=kind(time)) * particlemass + case(12) ! Escape_velocity + call calc_escape_velocities(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),xyzmh_ptmass,phii,epoti,v_esci) + arr(k,i) = v_esci case(13) !to calculate JstarS - rhopart = rhoh(xyzh(4,i), particlemass) - rho_cgs = rhopart*unit_density - !call equationofstate to obtain temperature and store it in tempi - call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) JstarS = 0. !nH_tot is needed to normalize JstarS nH_tot = rho_cgs/mass_per_H @@ -1413,97 +1486,22 @@ subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu) print *,'eps = ',eps print *,'JstarS = ',JstarS endif - quant(k,i) = JstarS - - case(0) ! Skip - quant(k,i) = 0. - - case(1,9) ! Total energy (kin + pot + therm) - rhopart = rhoh(xyzh(4,i), particlemass) - call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,& - epoti,ekini,egasi,eradi,ereci,dum1) - if (quantities_to_calculate(k)==1) then - call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),ethi) - quant(k,i) = (ekini + epoti + ethi) / particlemass ! Specific energy - elseif (quantities_to_calculate(k)==9) then - quant(k,i) = (ekini + epoti) / particlemass ! Specific energy - endif - - case(2) ! Mach number - rhopart = rhoh(xyzh(4,i), particlemass) - call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) - quant(k,i) = distance(vxyzu(1:3,i)) / spsoundi - - case(3) ! Opacity from MESA tables - rhopart = rhoh(xyzh(4,i), particlemass) - call ionisation_fraction(rhopart*unit_density,eos_vars(itemp,i),X_in,1.-X_in-Z_in,dum1,dum2,dum3,dum4,dum5) - if (ieos == 10) then - call get_eos_kappa_mesa(rhopart*unit_density,eos_vars(itemp,i),kappai,kappat,kappar) - quant(k,i) = kappai - else - quant(k,i) = 0. - endif - - case(4,11,12) ! Gas omega w.r.t. effective CoM - xyz_a = xyzh(1:3,i) - com_xyz(1:3) - vxyz_a = vxyzu(1:3,i) - com_vxyz(1:3) - quant(k,i) = (-xyz_a(2) * vxyz_a(1) + xyz_a(1) * vxyz_a(2)) / dot_product(xyz_a(1:2), xyz_a(1:2)) - - case(5) ! Fractional difference between gas and orbital omega - xyz_a = xyzh(1:3,i) - com_xyz(1:3) - vxyz_a = vxyzu(1:3,i) - com_vxyz(1:3) - quant(k,i) = (-xyz_a(2) * vxyz_a(1) + xyz_a(1) * vxyz_a(2)) / dot_product(xyz_a(1:2), xyz_a(1:2)) - quant(k,i) = (quant(k,i) - omega_orb) / omega_orb - - case(6,7) ! Calculate MESA EoS entropy - entropyi = 0. - rhopart = rhoh(xyzh(4,i), particlemass) - call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) - if (ieos==10) then - call getvalue_mesa(rhopart*unit_density,vxyzu(4,i)*unit_ergg,3,pgas,ierr) ! Get gas pressure - mu = rhopart*unit_density * Rg * eos_vars(itemp,i) / pgas - entropyi = entropy(rhopart*unit_density,ponrhoi*rhopart*unit_pressure,mu,3,vxyzu(4,i)*unit_ergg,ierr) - elseif (ieos==2) then - entropyi = entropy(rhopart*unit_density,ponrhoi*rhopart*unit_pressure,gmw,1) - endif - - if (quantities_to_calculate(k) == 7) then - if (dump_number == 0) then - init_entropy(i) = entropyi ! Store initial entropy on each particle - endif - quant(k,i) = entropyi/init_entropy(i) - 1. - elseif (quantities_to_calculate(k) == 6) then - quant(k,i) = entropyi - endif - - case(8) ! Specific recombination energy - rhopart = rhoh(xyzh(4,i), particlemass) - call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) - call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),ethi) - quant(k,i) = vxyzu(4,i) - ethi / particlemass ! Specific energy - - case(10) ! Mass coordinate - quant(k,iorder(i)) = real(i,kind=kind(time)) * particlemass - - case(14) ! Escape_velocity - call calc_escape_velocities(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),xyzmh_ptmass,phii,epoti,v_esci) - quant(k,i) = v_esci + arr(k,i) = JstarS case default - print*,"Error: Requested quantity is invalid." - stop + call fatal('analysis_common_envelope','Requested quantity is invalid.') end select enddo enddo - open(newunit=iu,file=trim(dumpfile)//".divv",status='replace',form='unformatted') - do k=1,4 - write(iu) (quant(k,i),i=1,npart) + ! Open files + do k=1,Nquant + open(newunit=iu(k),file=trim(dumpfile)//"."//trim(labels(quants(k)))//".extras",status='replace',form='unformatted') + write(iu(k)) (arr(k,i),i=1,npart) + close(iu(k)) enddo - close(iu) - deallocate(quant) + deallocate(arr) -end subroutine output_divv_files +end subroutine output_extra_quantities @@ -3794,7 +3792,7 @@ subroutine calc_gas_energies(particlemass,poten,xyzh,vxyzu,rad,xyzmh_ptmass,phii call get_accel_sink_gas(nptmass,xyzh(1),xyzh(2),xyzh(3),xyzh(4),xyzmh_ptmass,fxi,fyi,fzi,phii) epoti = 2.*poten + particlemass * phii ! For individual particles, need to multiply 2 to poten to get \sum_j G*mi*mj/r ekini = particlemass * 0.5 * dot_product(vxyzu(1:3),vxyzu(1:3)) - + egasradi = 0. egasi = 0. ereci = 0. if (do_radiation) then From d321d38fb452e0da09e5265144bc83dcab91a1e8 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Wed, 4 Sep 2024 16:37:07 +0200 Subject: [PATCH 12/19] bug fix: add forgotten argument to radxi_from_Trad --- src/setup/set_star_utils.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/setup/set_star_utils.f90 b/src/setup/set_star_utils.f90 index 26e32e4cc..9d4c8e7ab 100644 --- a/src/setup/set_star_utils.f90 +++ b/src/setup/set_star_utils.f90 @@ -440,7 +440,7 @@ subroutine set_star_thermalenergy(ieos,den,pres,r,npts,npart,xyzh,vxyzu,rad,eos_ endif if (do_radiation) then vxyzu(4,i) = ugas_from_Tgas(tempi,gamma,gmw) - rad(iradxi,i) = radxi_from_Trad(tempi) + rad(iradxi,i) = radxi_from_Trad(densi,tempi) else vxyzu(4,i) = eni / unit_ergg endif From 3483fac0dc9aa71a1dc5e133ca5d69c81ebc4c66 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Thu, 5 Sep 2024 01:15:10 +1000 Subject: [PATCH 13/19] merge changes --- src/utils/analysis_common_envelope.f90 | 34 +++++++++++--------------- 1 file changed, 14 insertions(+), 20 deletions(-) diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index 199198056..7a11023b9 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -1478,10 +1478,6 @@ subroutine output_extra_quantities(time,dumpfile,npart,particlemass,xyzh,vxyzu) print *,'epsC = ',epsC print *,'tempi = ',tempi print *,'S = ',S - print *,'pC =',pC - print *,'psat_C(tempi) = ',psat_C(tempi) - print *,'nucleation(idmu,i) = ',nucleation(idmu,i) - print *,'nucleation(idgamma,i) = ',nucleation(idgamma,i) print *,'taustar = ',taustar print *,'eps = ',eps print *,'JstarS = ',JstarS @@ -2764,8 +2760,7 @@ subroutine recombination_stats(time,num,npart,particlemass,xyzh,vxyzu) real :: etoti,ekini,egasi,eradi,ereci,epoti,ethi,phii,dum,rhopart,& ponrhoi,spsoundi,tempi,pressure,temperature,xh0,xh1,xhe0,xhe1,xhe2 character(len=40) :: data_formatter,logical_format - logical, allocatable :: isbound(:) - integer, allocatable :: H_state(:),He_state(:) + integer, allocatable :: H_state(:),He_state(:),isbound(:) integer :: i real, parameter :: recomb_th=0.05 @@ -2777,17 +2772,16 @@ subroutine recombination_stats(time,num,npart,particlemass,xyzh,vxyzu) rhopart = rhoh(xyzh(4,i), particlemass) call equationofstate(ieos,ponrhoi,spsoundi,rhopart,xyzh(1,i),xyzh(2,i),xyzh(3,i),tempi,vxyzu(4,i)) call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,epoti,ekini,egasi,eradi,ereci,dum) - call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),ethi) - etoti = ekini + epoti + ethi + call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),ethi,rad(:,i)) + etoti = ekini + epoti! + ethi - call get_eos_pressure_temp_mesa(rhopart*unit_density,vxyzu(4,i)*unit_ergg,pressure,temperature) ! This should depend on ieos - call ionisation_fraction(rhopart*unit_density,temperature,X_in,1.-X_in-Z_in,xh0,xh1,xhe0,xhe1,xhe2) + call ionisation_fraction(rhopart*unit_density,tempi,X_in,1.-X_in-Z_in,xh0,xh1,xhe0,xhe1,xhe2) ! Is unbound? if (etoti > 0.) then - isbound(i) = .false. + isbound(i) = 0 else - isbound(i) = .true. + isbound(i) = 1 endif ! H ionisation state @@ -2811,8 +2805,8 @@ subroutine recombination_stats(time,num,npart,particlemass,xyzh,vxyzu) endif enddo - write(data_formatter, "(a,I5,a)") "(es18.10e3,", npart, "(1x,i1))" ! Time column plus npart columns - write(logical_format, "(a,I5,a)") "(es18.10e3,", npart, "(1x,L))" ! Time column plus npart columns + write(data_formatter, "(a,I7,a)") "(es18.10e3,", npart, "(1x,i1))" ! Time column plus npart columns + write(logical_format, "(a,I7,a)") "(es18.10e3,", npart, "(1x,i1))" ! Time column plus npart columns if (num == 0) then ! Write header line open(unit=1000,file="H_state.ev",status='replace') @@ -2830,13 +2824,13 @@ subroutine recombination_stats(time,num,npart,particlemass,xyzh,vxyzu) write(1000,data_formatter) time,H_state(:) close(unit=1000) - open(unit=1000,file="He_state.ev", position='append') - write(1000,data_formatter) time,He_state(:) - close(unit=1000) + open(unit=1001,file="He_state.ev", position='append') + write(1001,data_formatter) time,He_state(:) + close(unit=1001) - open(unit=1000,file="isbound.ev", position='append') - write(1000,logical_format) time,isbound(:) - close(unit=1000) + open(unit=1002,file="isbound.ev", position='append') + write(1002,logical_format) time,isbound(:) + close(unit=1002) deallocate(isbound,H_state,He_state) From ae6d48dc391d8a4917bc5ff9ec22f97d6b39b699 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Mon, 9 Sep 2024 13:52:13 -0500 Subject: [PATCH 14/19] (eos) fix logic in entropy function --- src/main/eos.f90 | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index c0323a254..39d79b6cc 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -953,13 +953,16 @@ function entropy(rho,pres,mu_in,ientropy,eint_in,ierr,T_in,Trad_in) entropy = kb_on_mh / mu * log(temp**1.5/rho) case(2) ! Include both gas and radiation contributions (up to additive constants) - if (present(Trad_in)) then - Trad = Trad_in - else - if (.not. present(T_in)) then - call get_idealgasplusrad_tempfrompres(pres,rho,mu,temp) ! First solve for temp from rho and pres - Trad = temp ! assume thermal equilibrium + if (present(T_in)) then + temp = T_in + if (present(Trad_in)) then + Trad = Trad_in + else + Trad = temp endif + else + call get_idealgasplusrad_tempfrompres(pres,rho,mu,temp) ! First solve for temp from rho and pres + Trad = temp endif ! check temp if (temp < tiny(0.)) call warning('entropy','temperature = 0 will give minus infinity with s entropy') From 9daca12696c9752763a1069020b2eb503e973cee Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Mon, 9 Sep 2024 13:52:44 -0500 Subject: [PATCH 15/19] (eos) get_p_from_rho_s should use Rg not kb_on_mh in temperature solving --- src/main/eos.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 39d79b6cc..7fe20028d 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -1051,7 +1051,7 @@ end subroutine get_rho_from_p_s !+ !----------------------------------------------------------------------- subroutine get_p_from_rho_s(ieos,S,rho,mu,P,temp) - use physcon, only:kb_on_mh,radconst,rg,mass_proton_cgs,kboltz + use physcon, only:kb_on_mh,radconst,Rg,mass_proton_cgs,kboltz use io, only:fatal use eos_idealplusrad, only:get_idealgasplusrad_tempfrompres,get_idealplusrad_pres use units, only:unit_density,unit_pressure,unit_ergg @@ -1072,7 +1072,7 @@ subroutine get_p_from_rho_s(ieos,S,rho,mu,P,temp) select case (ieos) case (2,5,17) temp = (cgsrho * exp(mu*cgss*mass_proton_cgs))**(2./3.) - cgsP = cgsrho*kb_on_mh*temp / mu + cgsP = cgsrho*Rg*temp / mu case (12) corr = huge(corr) do while (abs(corr) > eoserr .and. niter < nitermax) From b10c7db073bf148750c3d57d939c0412a3573745 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Mon, 9 Sep 2024 14:10:25 -0500 Subject: [PATCH 16/19] delete or replace a bunch of kb_on_mh's --- src/main/eos.f90 | 10 +++++----- src/setup/set_star_utils.f90 | 4 ++-- src/setup/setup_star.f90 | 2 +- src/tests/test_gr.f90 | 2 +- 4 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/main/eos.f90 b/src/main/eos.f90 index 7fe20028d..5df8d12ca 100644 --- a/src/main/eos.f90 +++ b/src/main/eos.f90 @@ -840,7 +840,7 @@ end subroutine calc_rec_ene !+ !----------------------------------------------------------------------- subroutine calc_temp_and_ene(eos_type,rho,pres,ene,temp,ierr,guesseint,mu_local,X_local,Z_local) - use physcon, only:kb_on_mh + use physcon, only:Rg 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 @@ -861,7 +861,7 @@ subroutine calc_temp_and_ene(eos_type,rho,pres,ene,temp,ierr,guesseint,mu_local, if (present(Z_local)) Z = Z_local select case(eos_type) case(2,5,17) ! Ideal gas - temp = pres / (rho * kb_on_mh) * mu + temp = pres / (rho * Rg) * mu ene = pres / ( (gamma-1.) * rho) case(12) ! Ideal gas + radiation call get_idealgasplusrad_tempfrompres(pres,rho,mu,temp) @@ -889,7 +889,7 @@ end subroutine calc_temp_and_ene !+ !----------------------------------------------------------------------- subroutine calc_rho_from_PT(eos_type,pres,temp,rho,ierr,mu_local,X_local,Z_local) - use physcon, only:kb_on_mh + use physcon, only:Rg use eos_idealplusrad, only:get_idealplusrad_rhofrompresT use eos_mesa, only:get_eos_eT_from_rhop_mesa use eos_gasradrec, only:calc_uT_from_rhoP_gasradrec @@ -910,7 +910,7 @@ subroutine calc_rho_from_PT(eos_type,pres,temp,rho,ierr,mu_local,X_local,Z_local if (present(Z_local)) Z = Z_local select case(eos_type) case(2) ! Ideal gas - rho = pres / (temp * kb_on_mh) * mu + rho = pres / (temp * Rg) * mu case(12) ! Ideal gas + radiation call get_idealplusrad_rhofrompresT(pres,temp,mu,rho) case default @@ -1051,7 +1051,7 @@ end subroutine get_rho_from_p_s !+ !----------------------------------------------------------------------- subroutine get_p_from_rho_s(ieos,S,rho,mu,P,temp) - use physcon, only:kb_on_mh,radconst,Rg,mass_proton_cgs,kboltz + use physcon, only:radconst,Rg,mass_proton_cgs,kboltz use io, only:fatal use eos_idealplusrad, only:get_idealgasplusrad_tempfrompres,get_idealplusrad_pres use units, only:unit_density,unit_pressure,unit_ergg diff --git a/src/setup/set_star_utils.f90 b/src/setup/set_star_utils.f90 index 9d4c8e7ab..f616345e7 100644 --- a/src/setup/set_star_utils.f90 +++ b/src/setup/set_star_utils.f90 @@ -458,7 +458,7 @@ end subroutine set_star_thermalenergy !----------------------------------------------------------------------- subroutine solve_uT_profiles(eos_type,r,den,pres,Xfrac,Yfrac,regrid_core,temp,en,mu) use eos, only:get_mean_molecular_weight,calc_temp_and_ene - use physcon, only:radconst,kb_on_mh + use physcon, only:radconst,Rg integer, intent(in) :: eos_type real, intent(in) :: r(:),den(:),pres(:),Xfrac(:),Yfrac(:) logical, intent(in) :: regrid_core @@ -475,7 +475,7 @@ subroutine solve_uT_profiles(eos_type,r,den,pres,Xfrac,Yfrac,regrid_core,temp,en mu(i) = get_mean_molecular_weight(Xfrac(i),1.-Xfrac(i)-Yfrac(i)) ! only used in u, T calculation if ieos==2,12 if (i==1) then guessene = 1.5*pres(i)/den(i) ! initial guess - tempi = min((3.*pres(i)/radconst)**0.25, pres(i)*mu(i)/(den(i)*kb_on_mh)) ! guess for temperature + tempi = min((3.*pres(i)/radconst)**0.25, pres(i)*mu(i)/(den(i)*Rg)) ! guess for temperature else guessene = en(i-1) tempi = temp(i-1) diff --git a/src/setup/setup_star.f90 b/src/setup/setup_star.f90 index e6fa4a8eb..7f4ded993 100644 --- a/src/setup/setup_star.f90 +++ b/src/setup/setup_star.f90 @@ -33,7 +33,7 @@ module setup ! use io, only:fatal,error,warning,master use part, only:gravity,gr - use physcon, only:solarm,solarr,km,pi,c,kb_on_mh,radconst + use physcon, only:solarm,solarr,km,pi,c,radconst use options, only:nfulldump,iexternalforce,calc_erot,use_var_comp use timestep, only:tmax,dtmax use eos, only:ieos diff --git a/src/tests/test_gr.f90 b/src/tests/test_gr.f90 index 905a15a0d..bdc1ad0fa 100644 --- a/src/tests/test_gr.f90 +++ b/src/tests/test_gr.f90 @@ -468,7 +468,7 @@ subroutine test_cons2prim_i(x,v,dens,u,p,ncheck,nfail,errmax,tol) use part, only:ien_entropy,ien_etotal,ien_entropy_s use metric_tools, only:pack_metric,unpack_metric use eos, only:ieos,equationofstate,calc_temp_and_ene - use physcon, only:radconst,kb_on_mh + use physcon, only:radconst real, intent(in) :: x(1:3),v(1:3),dens,p,tol real, intent(inout) :: u From 300c3396afb42b4cc4639c1afbfa583bd34b64ec Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Mon, 9 Sep 2024 22:56:38 -0500 Subject: [PATCH 17/19] (CE-analysis) remove unused variables declared --- src/utils/analysis_common_envelope.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index 4b8a17d35..ba9fecc06 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -2758,7 +2758,7 @@ subroutine recombination_stats(time,num,npart,particlemass,xyzh,vxyzu) real, intent(in) :: time,particlemass real, intent(inout) :: xyzh(:,:),vxyzu(:,:) real :: etoti,ekini,egasi,eradi,ereci,epoti,ethi,phii,dum,rhopart,& - ponrhoi,spsoundi,tempi,pressure,temperature,xh0,xh1,xhe0,xhe1,xhe2 + ponrhoi,spsoundi,tempi,xh0,xh1,xhe0,xhe1,xhe2 character(len=40) :: data_formatter,logical_format integer, allocatable :: H_state(:),He_state(:),isbound(:) integer :: i From 99730c9e18a4699e4f376e7149e143b564cbe90e Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Tue, 10 Sep 2024 10:24:16 -0500 Subject: [PATCH 18/19] (CE-analysis) update test script --- scripts/test_analysis_ce.sh | 2 ++ src/utils/analysis_common_envelope.f90 | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/scripts/test_analysis_ce.sh b/scripts/test_analysis_ce.sh index a149860b8..cd718ce17 100755 --- a/scripts/test_analysis_ce.sh +++ b/scripts/test_analysis_ce.sh @@ -29,6 +29,7 @@ SEP no 2 1.667 +0.6182 0.6984 0.0142 BOUND @@ -38,4 +39,5 @@ BOUND no 2 1.667 +0.6182 ENERGIES diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index ba9fecc06..04e1ba8a7 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -4488,7 +4488,7 @@ subroutine set_eos_options(analysis_to_perform) gamma = 5./3. call prompt('Enter gamma:',gamma,0.) gmw = 0.618212823 - call prompt('Enter mean molecular weight for gas+rad EoS:',gmw,0.) + call prompt('Enter mean molecular weight:',gmw,0.) case(10,20) gamma = 5./3. X_in = 0.69843 From e0fdcf94510e46b442edb37fece2f30e6a59c057 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Thu, 12 Sep 2024 12:37:44 -0500 Subject: [PATCH 19/19] (CE-analysis) clean up newly unbound particles --- src/utils/analysis_common_envelope.f90 | 100 +++++++++---------------- 1 file changed, 34 insertions(+), 66 deletions(-) diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index 04e1ba8a7..692c4ef43 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -2494,89 +2494,59 @@ subroutine unbound_profiles(time,num,npart,particlemass,xyzh,vxyzu) integer, intent(in) :: npart,num real, intent(in) :: time,particlemass real, intent(inout) :: xyzh(:,:),vxyzu(:,:) - integer, dimension(4) :: npart_hist - real, dimension(5,npart) :: dist_part,rad_part + integer, dimension(2) :: nunbound + real, dimension(2,npart) :: dist_part,rad_part real, dimension(:), allocatable :: hist_var - real :: etoti,ekini,ereci,egasi,eradi,epoti,phii,dum,maxloga,minloga - character(len=18), dimension(4) :: grid_file + real :: e_kp,e_kpt,etoti,ekini,ereci,egasi,eradi,epoti,phii,sep,maxloga,minloga + character(len=18), dimension(2) :: grid_file character(len=40) :: data_formatter - logical, allocatable, save :: prev_unbound(:,:),prev_bound(:,:) - integer :: i,unitnum,nbins + logical, allocatable, save :: prev_bound(:,:) + integer :: i,j,unitnum,nbins,maxj call compute_energies(time) - npart_hist = 0 ! Stores number of particles fulfilling each of the four bound/unbound criterion + nunbound = 0 ! Stores number of particles that have become newly unbound in this dump according to e_kp or e_kpt criterion nbins = 500 - rad_part = 0. ! (4,npart_hist)-array storing separations of particles - dist_part = 0. + rad_part = 0. ! (2,npart_hist)-array storing separations of newly unbound particles + dist_part = 0. ! Array of ones with size of 2? minloga = 0.5 maxloga = 4.3 allocate(hist_var(nbins)) - grid_file = (/ 'grid_unbound_th.ev', & - 'grid_unbound_kp.ev', & - ' grid_bound_kpt.ev', & - ' grid_bound_kp.ev' /) + grid_file = (/ 'grid_unbound_th.ev', 'grid_unbound_kp.ev' /) if (dump_number == 0) then allocate(prev_bound(2,npart)) - allocate(prev_unbound(2,npart)) - prev_bound = .false. - prev_unbound = .false. + prev_bound = .true. ! all particles bound to begin with endif - do i=1,npart - if (.not. isdead_or_accreted(xyzh(4,i))) then - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,& - epoti,ekini,egasi,eradi,ereci,dum) - etoti = ekini + epoti + egasi + eradi - - ! Ekin + Epot + Eth > 0 - if ((etoti > 0.) .and. (.not. prev_unbound(1,i))) then - npart_hist(1) = npart_hist(1) + 1 ! Keeps track of number of particles that have become newly unbound in this dump - rad_part(1,npart_hist(1)) = separation(xyzh(1:3,i),xyzmh_ptmass(1:3,1)) - dist_part(1,npart_hist(1)) = 1. ! Array of ones with size of npart_hist(1)? - prev_unbound(1,i) = .true. - elseif (etoti < 0.) then - prev_unbound(1,i) = .false. - endif - - ! Ekin + Epot > 0 - if ((ekini + epoti > 0.) .and. (.not. prev_unbound(2,i))) then - npart_hist(2) = npart_hist(2) + 1 - rad_part(2,npart_hist(2)) = separation(xyzh(1:3,i),xyzmh_ptmass(1:3,1)) - dist_part(2,npart_hist(2)) = 1. - prev_unbound(2,i) = .true. - elseif (ekini + epoti < 0.) then - prev_unbound(2,i) = .false. - endif - - ! Ekin + Epot + Eth < 0 - if ((etoti < 0.) .and. (.not. prev_bound(1,i))) then - npart_hist(3) = npart_hist(3) + 1 - rad_part(3,npart_hist(3)) = separation(xyzh(1:3,i),xyzmh_ptmass(1:3,1)) - dist_part(3,npart_hist(3)) = 1. - prev_bound(1,i) = .true. - elseif (etoti > 0.) then - prev_bound(1,i) = .false. - endif - - ! Ekin + Epot < 0 - if ((ekini + epoti < 0.) .and. (.not. prev_bound(2,i))) then - npart_hist(4) = npart_hist(4) + 1 - rad_part(4,npart_hist(4)) = separation(xyzh(1:3,i),xyzmh_ptmass(1:3,1)) - dist_part(4,npart_hist(4)) = 1. - prev_bound(2,i) = .true. - elseif (ekini + epoti > 0.) then - prev_bound(2,i) = .false. - endif + if (isdead_or_accreted(xyzh(4,i))) cycle + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),rad(:,i),xyzmh_ptmass,phii,& + epoti,ekini,egasi,eradi,ereci,etoti) + e_kp = ekini + epoti + e_kpt = e_kp + egasi + eradi + + if (e_kp > 0. .and. prev_bound(2,i)) then ! newly bound by e_kp criterion + maxj = 2 + sep = separation(xyzh(1:3,i),xyzmh_ptmass(1:3,1)) + elseif (e_kpt > 0. .and. prev_bound(1,i)) then ! newly bound by e_kpt but not e_kp criterion + maxj = 1 + sep = separation(xyzh(1:3,i),xyzmh_ptmass(1:3,1)) + else ! particle state has not changed + cycle endif + + do j = 1,maxj + nunbound(j) = nunbound(j) + 1 + rad_part(j,nunbound(j)) = sep + dist_part(j,nunbound(j)) = 1. + prev_bound(j,i) = .false. + enddo enddo - do i=1,4 - call histogram_setup(rad_part(i,1:npart_hist(i)),dist_part(i,1:npart_hist(i)),hist_var,npart_hist(i),maxloga,minloga,nbins,& + do i=1,2 + call histogram_setup(rad_part(i,1:nunbound(i)),dist_part(i,1:nunbound(i)),hist_var,nunbound(i),maxloga,minloga,nbins,& .false.,.true.) - write(data_formatter, "(a,I5,a)") "(", nbins+1, "(3x,es18.10e3,1x))" ! Time column plus nbins columns if (num == 0) then ! Write header line @@ -2586,10 +2556,8 @@ subroutine unbound_profiles(time,num,npart,particlemass,xyzh,vxyzu) endif open(newunit=unitnum,file=trim(adjustl(grid_file(i))), position='append') - write(unitnum,"()") write(unitnum,data_formatter) time,hist_var(:) - close(unit=unitnum) enddo deallocate(hist_var)