From d6af633ca06ad14d50bfaa4b9acbd3229eb5c528 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Thu, 18 Jan 2024 14:18:09 +0100 Subject: [PATCH 1/8] (CE-analysis) add comments to potential energy calculation in calculate_energies --- src/utils/analysis_common_envelope.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index 6cd1c6c27..aba02cdbe 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -793,7 +793,7 @@ subroutine calculate_energies(time,npart,particlemass,xyzh,vxyzu) real :: rhopart,ponrhoi,spsoundi,tempi,r_ij,radvel real, dimension(3) :: rcrossmv character(len=17), allocatable :: columns(:) - integer :: i, j, ncols + integer :: i,j,ncols logical :: inearsink integer, parameter :: ie_tot = 1 integer, parameter :: ie_pot = ie_tot + 1 @@ -826,9 +826,9 @@ subroutine calculate_energies(time,npart,particlemass,xyzh,vxyzu) ' pot energy',& ' kin energy',& 'therm energy',& - ' sink pot',& + ' sink pot',& ! does not include sink-gas potential energy ' sink kin',& - ' sink orb',& + ' sink orb',& ! sink kin + sink pot ' comp orb',& ' env pot',& ' env energy',& @@ -914,7 +914,7 @@ subroutine calculate_energies(time,npart,particlemass,xyzh,vxyzu) do j=i+1,nptmass if (xyzmh_ptmass(4,j) > 0.) then r_ij = separation(xyzmh_ptmass(1:3,i),xyzmh_ptmass(1:3,j)) - encomp(ipot_sink) = encomp(ipot_sink) - xyzmh_ptmass(4,i) * xyzmh_ptmass(4,j) / r_ij + encomp(ipot_sink) = encomp(ipot_sink) - xyzmh_ptmass(4,i) * xyzmh_ptmass(4,j) / r_ij ! Newtonian expression is fine as long as rij > hsofti + hsoftj if (i==1 .and. j==2) encomp(iorb_comp) = encomp(iorb_comp) - xyzmh_ptmass(4,i) * xyzmh_ptmass(4,j) / r_ij endif enddo From aa546cda0c6f32ec8e318a5c49368312cfd398f5 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Thu, 18 Jan 2024 14:18:46 +0100 Subject: [PATCH 2/8] (CE-analysis) add sound speed to energy_profile --- src/utils/analysis_common_envelope.f90 | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index aba02cdbe..629645cb2 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -2064,11 +2064,12 @@ subroutine energy_profile(time,npart,particlemass,xyzh,vxyzu) if (dump_number==0) then iquantity = 1 use_mass_coord = .false. - print "(4(/,a))",'1. Energy',& + print "(5(/,a))",'1. Energy',& '2. Entropy',& '3. Bernoulli energy',& - '4. Ion fractions' - call prompt("Select quantity to calculate",iquantity,1,4) + '4. Ion fractions',& + '5. Sound speed' + call prompt("Select quantity to calculate",iquantity,1,5) call prompt("Bin in mass coordinates instead of radius?",use_mass_coord) endif @@ -2087,7 +2088,7 @@ subroutine energy_profile(time,npart,particlemass,xyzh,vxyzu) call compute_energies(time) ! Allocate arrays for single variable outputs - if ( (iquantity==1) .or. (iquantity==2) .or. (iquantity==3) ) then + if (iquantity==1 .or. iquantity==2 .or. iquantity==3 .or. iquantity==5) then nvars = 1 else nvars = 5 @@ -2127,6 +2128,9 @@ subroutine energy_profile(time,npart,particlemass,xyzh,vxyzu) ' # HeI', & ' # HeII', & ' # HeIII' /) + case(5) ! Sound speed + filename = ' grid_cs.ev' + headerline = '# cs profile ' end select allocate(iorder(npart)) @@ -2174,6 +2178,8 @@ subroutine energy_profile(time,npart,particlemass,xyzh,vxyzu) quant(i,3) = xhe0 quant(i,4) = xhe1 quant(i,5) = xhe2 + case(5) ! Sound speed + quant(i,1) = spsoundi end select enddo From 55f23d5a3add6517a53624fdc20479497feadebd Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Tue, 30 Apr 2024 13:21:34 +0200 Subject: [PATCH 3/8] (CE-analysis) add radiation energy to thermal energy when using radiation --- src/main/ionization.f90 | 7 ++- src/utils/analysis_common_envelope.f90 | 69 ++++++++++++-------------- 2 files changed, 38 insertions(+), 38 deletions(-) diff --git a/src/main/ionization.f90 b/src/main/ionization.f90 index ebc536639..88f155561 100644 --- a/src/main/ionization.f90 +++ b/src/main/ionization.f90 @@ -338,13 +338,15 @@ end subroutine get_erec_components ! gas particle. Inputs and outputs in code units !+ !---------------------------------------------------------------- -subroutine calc_thermal_energy(particlemass,ieos,xyzh,vxyzu,presi,tempi,gamma,ethi) - use part, only:rhoh +subroutine calc_thermal_energy(particlemass,ieos,xyzh,vxyzu,presi,tempi,gamma,ethi,radprop) + use dim, only:do_radiation + use part, only:rhoh,iradxi use eos_idealplusrad, only:get_idealgasplusrad_tempfrompres,get_idealplusrad_enfromtemp use physcon, only:radconst,Rg 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),gamma + real, intent(in), optional :: radprop real, intent(out) :: ethi real :: hi,densi_cgs,mui @@ -357,6 +359,7 @@ subroutine calc_thermal_energy(particlemass,ieos,xyzh,vxyzu,presi,tempi,gamma,et ethi = particlemass * ethi / unit_ergg case default ! assuming internal energy = thermal energy ethi = particlemass * vxyzu(4) + if (do_radiation) ethi = ethi + particlemass*radprop(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 629645cb2..4386063b3 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -20,24 +20,26 @@ module analysis ! sortutils, table_utils, units, vectorutils ! - use part, only:xyzmh_ptmass,vxyz_ptmass,nptmass,poten,ihsoft,ihacc,& - rhoh,nsinkproperties,maxvxyzu,maxptmass,isdead_or_accreted - use units, only:print_units,umass,utime,udist,unit_ergg,unit_density,& - unit_pressure,unit_velocity,unit_Bfield,unit_energ - use physcon, only:gg,pi,c,Rg - use io, only:fatal - use prompting, only:prompt - use centreofmass, only:get_centreofmass, reset_centreofmass - use energies, only:compute_energies,ekin,etherm,epot,etot - use ptmass, only:get_accel_sink_gas,get_accel_sink_sink - use kernel, only:kernel_softening,radkern,wkern,cnormk - use eos, only:equationofstate,ieos,init_eos,X_in,Z_in,gmw,get_spsound,done_init_eos - use eos_gasradrec,only:irecomb - use eos_mesa, only:get_eos_kappa_mesa,get_eos_pressure_temp_mesa,& - get_eos_various_mesa,get_eos_pressure_temp_gamma1_mesa - use setbinary, only:Rochelobe_estimate,L1_point - use sortutils, only:set_r2func_origin,r2func_origin,indexxfunc - use table_utils, only:logspace + use part, only:xyzmh_ptmass,vxyz_ptmass,nptmass,poten,ihsoft,ihacc,& + rhoh,nsinkproperties,maxvxyzu,maxptmass,isdead_or_accreted + 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 + use physcon, only:gg,pi,c,Rg + use io, only:fatal + use prompting, only:prompt + use centreofmass, only:get_centreofmass, reset_centreofmass + use energies, only:compute_energies,ekin,etherm,epot,etot + use ptmass, only:get_accel_sink_gas,get_accel_sink_sink + use kernel, only:kernel_softening,radkern,wkern,cnormk + use ionization_mod,only:calc_thermal_energy + use eos, only:equationofstate,ieos,init_eos,X_in,Z_in,gmw,get_spsound,done_init_eos + use eos_gasradrec, only:irecomb + use eos_mesa, only:get_eos_kappa_mesa,get_eos_pressure_temp_mesa,& + get_eos_various_mesa,get_eos_pressure_temp_gamma1_mesa + use setbinary, only:Rochelobe_estimate,L1_point + use sortutils, only:set_r2func_origin,r2func_origin,indexxfunc + use table_utils, only:logspace implicit none character(len=20), parameter, public :: analysistype = 'common_envelope' integer :: analysis_to_perform @@ -623,9 +625,8 @@ end subroutine m_vs_t !+ !---------------------------------------------------------------- subroutine bound_mass(time,npart,particlemass,xyzh,vxyzu) - use part, only:eos_vars,itemp + use part, only:eos_vars,itemp,radprop use ptmass, only:get_accel_sink_gas - use ionization_mod, only:calc_thermal_energy use vectorutils, only:cross_product3D integer, intent(in) :: npart real, intent(in) :: time,particlemass @@ -702,7 +703,7 @@ subroutine bound_mass(time,npart,particlemass,xyzh,vxyzu) 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,gamma,ethi) + call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,tempi,gamma,ethi,radprop) 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 @@ -1382,7 +1383,7 @@ subroutine output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu) 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:calc_thermal_energy,ionisation_fraction + 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 integer, intent(in) :: npart @@ -1663,7 +1664,7 @@ subroutine track_particle(time,particlemass,xyzh,vxyzu) use part, only:eos_vars,itemp use eos, only:entropy use mesa_microphysics, only:getvalue_mesa - use ionization_mod, only:calc_thermal_energy,ionisation_fraction + use ionization_mod, only:ionisation_fraction real, intent(in) :: time,particlemass real, intent(inout) :: xyzh(:,:),vxyzu(:,:) integer, parameter :: nparttotrack=10,ncols=17 @@ -1888,7 +1889,7 @@ end subroutine tconv_profile !---------------------------------------------------------------- subroutine recombination_tau(time,npart,particlemass,xyzh,vxyzu) use part, only:eos_vars,itemp - use ionization_mod, only:calc_thermal_energy,ionisation_fraction + use ionization_mod, only:ionisation_fraction integer, intent(in) :: npart real, intent(in) :: time,particlemass real, intent(inout) :: xyzh(:,:),vxyzu(:,:) @@ -1979,7 +1980,6 @@ end subroutine recombination_tau !---------------------------------------------------------------- subroutine energy_hist(time,npart,particlemass,xyzh,vxyzu) use part, only:eos_vars,itemp - use ionization_mod, only:calc_thermal_energy integer, intent(in) :: npart real, intent(in) :: time,particlemass real, intent(inout) :: xyzh(:,:),vxyzu(:,:) @@ -2044,7 +2044,7 @@ subroutine energy_profile(time,npart,particlemass,xyzh,vxyzu) use part, only:eos_vars,itemp use eos, only:entropy use mesa_microphysics, only:getvalue_mesa - use ionization_mod, only:calc_thermal_energy,ionisation_fraction + use ionization_mod, only:ionisation_fraction integer, intent(in) :: npart real, intent(in) :: time,particlemass real, intent(inout) :: xyzh(:,:),vxyzu(:,:) @@ -2288,7 +2288,6 @@ end subroutine rotation_profile !---------------------------------------------------------------- subroutine velocity_histogram(time,num,npart,particlemass,xyzh,vxyzu) use part, only:eos_vars,itemp - use ionization_mod, only:calc_thermal_energy real, intent(in) :: time,particlemass integer, intent(in) :: npart,num real, intent(inout) :: xyzh(:,:),vxyzu(:,:) @@ -2571,7 +2570,6 @@ end subroutine planet_profile !+ !---------------------------------------------------------------- subroutine unbound_profiles(time,num,npart,particlemass,xyzh,vxyzu) - use ionization_mod, only:calc_thermal_energy integer, intent(in) :: npart,num real, intent(in) :: time,particlemass real, intent(inout) :: xyzh(:,:),vxyzu(:,:) @@ -2690,7 +2688,7 @@ end subroutine unbound_profiles !+ !---------------------------------------------------------------- subroutine unbound_ionfrac(time,npart,particlemass,xyzh,vxyzu) - use ionization_mod, only:calc_thermal_energy,get_xion,ionisation_fraction + use ionization_mod, only:get_xion,ionisation_fraction integer, intent(in) :: npart real, intent(in) :: time,particlemass real, intent(inout) :: xyzh(:,:),vxyzu(:,:) @@ -2766,7 +2764,7 @@ end subroutine unbound_ionfrac !---------------------------------------------------------------- subroutine unbound_temp(time,npart,particlemass,xyzh,vxyzu) use part, only:eos_vars,itemp - use ionization_mod, only:calc_thermal_energy,get_xion + use ionization_mod, only:get_xion integer, intent(in) :: npart real, intent(in) :: time,particlemass real, intent(inout) :: xyzh(:,:),vxyzu(:,:) @@ -2840,7 +2838,7 @@ end subroutine unbound_temp !---------------------------------------------------------------- subroutine recombination_stats(time,num,npart,particlemass,xyzh,vxyzu) use part, only:eos_vars,itemp - use ionization_mod, only:calc_thermal_energy,ionisation_fraction + use ionization_mod, only:ionisation_fraction integer, intent(in) :: npart,num real, intent(in) :: time,particlemass real, intent(inout) :: xyzh(:,:),vxyzu(:,:) @@ -3038,7 +3036,6 @@ end subroutine sink_properties subroutine env_binding_ene(npart,particlemass,xyzh,vxyzu) use part, only:eos_vars,itemp - use ionization_mod, only:calc_thermal_energy integer, intent(in) :: npart real, intent(in) :: particlemass real, intent(inout) :: xyzh(:,:),vxyzu(:,:) @@ -3724,7 +3721,6 @@ end subroutine print_dump_numbers subroutine analyse_disk(num,npart,particlemass,xyzh,vxyzu) use part, only:eos_vars,itemp use extern_corotate, only:get_companion_force - use ionization_mod, only:calc_thermal_energy use vectorutils, only:cross_product3D integer, intent(in) :: num,npart real, intent(in) :: particlemass @@ -3855,14 +3851,14 @@ end subroutine get_gas_omega ! and internal energy of a gas particle. !+ !---------------------------------------------------------------- -subroutine calc_gas_energies(particlemass,poten,xyzh,vxyzu,xyzmh_ptmass,phii,epoti,ekini,einti,etoti) +subroutine calc_gas_energies(particlemass,poten,xyzh,vxyzu,radprop,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 + use part, only:nptmass,iradxi real, intent(in) :: particlemass real(4), intent(in) :: poten - real, dimension(4), intent(in) :: xyzh,vxyzu + real, intent(in) :: xyzh(:),vxyzu(:),radprop(:) real, dimension(5,nptmass), intent(in) :: xyzmh_ptmass real, intent(out) :: phii,epoti,ekini,einti,etoti real :: fxi,fyi,fzi @@ -3874,6 +3870,7 @@ subroutine calc_gas_energies(particlemass,poten,xyzh,vxyzu,xyzmh_ptmass,phii,epo 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) etoti = epoti + ekini + einti end subroutine calc_gas_energies From 11de37ef28b80213ee893a6e15681c9a9d00ad89 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Tue, 30 Apr 2024 13:36:51 +0200 Subject: [PATCH 4/8] (CE-analysis) add radprop to calc_gas_energy --- src/main/ionization.f90 | 2 +- src/utils/analysis_common_envelope.f90 | 44 ++++++++++++++------------ 2 files changed, 24 insertions(+), 22 deletions(-) diff --git a/src/main/ionization.f90 b/src/main/ionization.f90 index 88f155561..a2b914c5b 100644 --- a/src/main/ionization.f90 +++ b/src/main/ionization.f90 @@ -346,7 +346,7 @@ subroutine calc_thermal_energy(particlemass,ieos,xyzh,vxyzu,presi,tempi,gamma,et 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),gamma - real, intent(in), optional :: radprop + real, intent(in), optional :: radprop(:) real, intent(out) :: ethi real :: hi,densi_cgs,mui diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index 4386063b3..5a1018278 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -21,7 +21,8 @@ module analysis ! use part, only:xyzmh_ptmass,vxyz_ptmass,nptmass,poten,ihsoft,ihacc,& - rhoh,nsinkproperties,maxvxyzu,maxptmass,isdead_or_accreted + rhoh,nsinkproperties,maxvxyzu,maxptmass,isdead_or_accreted,& + 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 @@ -281,7 +282,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) xyz_a(1:3) = xyzh(1:3,i) - com_xyz(1:3) vxyz_a(1:3) = vxyzu(1:3,i) - com_vxyz(1:3) - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) histogram_data(1:3,i) = xyzh(1:3,i) histogram_data(4,i) = distance(xyz_a(1:3)) histogram_data(5,i) = epoti + ekini @@ -697,13 +698,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),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,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,gamma,ethi,radprop) + call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,tempi,gamma,ethi,radprop(:,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 @@ -859,7 +860,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),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) encomp(ipot_ps) = encomp(ipot_ps) + particlemass * phii @@ -1105,7 +1106,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),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,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)) @@ -1281,7 +1282,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),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,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 @@ -1519,7 +1520,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),xyzmh_ptmass,phii,epoti,ekini,einti,dum1) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,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),gamma,ethi) quant(k,i) = (ekini + epoti + ethi) / particlemass ! Specific energy @@ -1727,7 +1728,7 @@ subroutine track_particle(time,particlemass,xyzh,vxyzu) 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),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + 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),gamma,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) @@ -1927,7 +1928,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),xyzmh_ptmass,phii,epoti,ekini,einti,dum) ! Calculate total energy + 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_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rho_part(i),eos_vars(itemp,i),gamma,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 @@ -2006,7 +2007,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),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,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),gamma,ethi) else @@ -2153,7 +2154,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),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + 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),gamma,ethi) quant(i,1) = ekini + epoti + ethi case(2) ! Entropy @@ -2169,7 +2170,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),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,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) @@ -2301,7 +2302,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),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + 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),gamma,ethi) vr(i) = dot_product(xyzh(1:3,i),vxyzu(1:3,i)) / sqrt(dot_product(xyzh(1:3,i),xyzh(1:3,i))) @@ -2609,7 +2610,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),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + 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,tempi,gamma,ethi) etoti = ekini + epoti + ethi @@ -2717,7 +2718,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),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + 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,tempi,gamma,ethi) etoti = ekini + epoti + ethi @@ -2787,7 +2788,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),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + 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),gamma,ethi) etoti = ekini + epoti + ethi @@ -2857,7 +2858,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),xyzmh_ptmass,phii,epoti,ekini,einti,dum) + 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),gamma,ethi) etoti = ekini + epoti + ethi @@ -3102,7 +3103,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),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) rhopart = rhoh(xyzh(4,i), particlemass) @@ -3447,7 +3448,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),xyzmh_ptmass,dum1,dum2,dum3,dum4,etoti) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,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 @@ -3996,7 +3997,8 @@ 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),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) + call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,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 c24c077b1d317b680e7bda61a321496a4ed166ff Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Tue, 30 Apr 2024 13:46:08 +0200 Subject: [PATCH 5/8] (CE-analysis) get rid of gammas in calc_therm_energy --- src/utils/analysis_common_envelope.f90 | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index ff567ae41..7edc51774 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -704,7 +704,7 @@ subroutine bound_mass(time,npart,particlemass,xyzh,vxyzu) 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,gamma,ethi,radprop(:,i)) + call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,tempi,ethi,radprop(:,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 @@ -1734,7 +1734,7 @@ subroutine track_particle(time,particlemass,xyzh,vxyzu) ! 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),gamma,ethi) + 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) @@ -1934,7 +1934,7 @@ subroutine recombination_tau(time,npart,particlemass,xyzh,vxyzu) 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_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rho_part(i),eos_vars(itemp,i),gamma,ethi) + 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 j=j+1 @@ -2160,7 +2160,7 @@ subroutine energy_profile(time,npart,particlemass,xyzh,vxyzu) 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_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),gamma,ethi) + 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 if ((ieos==10) .and. (ientropy==2)) then @@ -2308,7 +2308,7 @@ subroutine velocity_histogram(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),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),gamma,ethi) + 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))) if (ekini+epoti > 0.) then @@ -2616,7 +2616,7 @@ subroutine unbound_profiles(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),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) - call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,tempi,gamma,ethi) + call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,tempi,ethi) etoti = ekini + epoti + ethi ! Ekin + Epot + Eth > 0 @@ -2724,7 +2724,7 @@ subroutine unbound_ionfrac(time,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),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,dum) - call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,tempi,gamma,ethi) + call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,tempi,ethi) etoti = ekini + epoti + ethi if ((etoti > 0.) .and. (.not. prev_unbound(i))) then @@ -2794,7 +2794,7 @@ subroutine unbound_temp(time,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),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_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),gamma,ethi) + call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),ethi) etoti = ekini + epoti + ethi if ((etoti > 0.) .and. (.not. prev_unbound(i))) then @@ -2864,7 +2864,7 @@ 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),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),gamma,ethi) + call calc_thermal_energy(particlemass,ieos,xyzh(:,i),vxyzu(:,i),ponrhoi*rhopart,eos_vars(itemp,i),ethi) 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 From 5a7532ef303e3b47775d8d5b68a076b213128417 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Mon, 6 May 2024 10:52:54 +0200 Subject: [PATCH 6/8] (CE-analysis) remove case 14 for saving particle therm. quantities into file --- src/utils/analysis_common_envelope.f90 | 54 +------------------------- 1 file changed, 2 insertions(+), 52 deletions(-) diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index 7edc51774..0a40d006f 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -62,9 +62,6 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) integer :: unitnum,i,ncols logical :: requires_eos_opts - !case 5 variables - real :: rhopart - !case 7 variables character(len=17), allocatable :: columns(:) @@ -76,17 +73,9 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) real, allocatable :: histogram_data(:,:) real :: ang_vel - real :: pres_1i, proint_1i, peint_1i, temp_1i - real :: troint_1i, teint_1i, entrop_1i, abad_1i, gamma1_1i, gam_1i - - !case 16 variables - real, allocatable :: thermodynamic_quantities(:,:) - real, allocatable :: radius_1i, dens_1i - - !chose analysis type if (dump_number==0) then - print "(41(a,/))", & + print "(40(a,/))", & ' 1) Sink separation', & ' 2) Bound and unbound quantities', & ' 3) Energies', & @@ -99,7 +88,6 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) '11) Profile of newly unbound particles', & '12) Sink properties', & '13) MESA EoS compute total entropy and other average td quantities', & - '14) MESA EoS save on file thermodynamical quantities for all particles', & '15) Gravitational drag on sinks', & '16) CoM of gas around primary core', & '17) Miscellaneous', & @@ -136,7 +124,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) xyzmh_ptmass,vxyz_ptmass,omega_corotate,dump_number) ! List of analysis options that require specifying EOS options - requires_eos_opts = any((/2,3,4,5,6,8,9,11,13,14,15,20,21,22,23,24,25,26,29,30,31,32,33,35,41/) == analysis_to_perform) + requires_eos_opts = any((/2,3,4,5,6,8,9,11,13,15,20,21,22,23,24,25,26,29,30,31,32,33,35,41/) == analysis_to_perform) if (dump_number == 0 .and. requires_eos_opts) call set_eos_options(analysis_to_perform) select case(analysis_to_perform) @@ -210,48 +198,10 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) call sink_properties(time,npart,particlemass,xyzh,vxyzu) case(13) !MESA EoS compute total entropy and other average thermodynamical quantities call bound_unbound_thermo(time,npart,particlemass,xyzh,vxyzu) - case(14) !MESA EoS save on file thermodynamical quantities for all particles - allocate(thermodynamic_quantities(5,npart)) - do i=1,npart - - !particle radius - radius_1i = distance(xyzh(1:3,i)) * udist - - !particles density in code units - rhopart = rhoh(xyzh(4,i), particlemass) - dens_1i = rhopart * unit_density - - !gets entropy for the current particle - call get_eos_various_mesa(rhopart*unit_density,vxyzu(4,i) * unit_ergg, & - pres_1i,proint_1i,peint_1i,temp_1i,troint_1i, & - teint_1i,entrop_1i,abad_1i,gamma1_1i,gam_1i) - - !stores everything in an array - thermodynamic_quantities(1,i) = radius_1i - thermodynamic_quantities(2,i) = dens_1i - thermodynamic_quantities(3,i) = pres_1i - thermodynamic_quantities(4,i) = temp_1i - thermodynamic_quantities(5,i) = entrop_1i - - enddo - ncols = 5 - allocate(columns(ncols)) - columns = (/' radius', & - ' density', & - ' pressure', & - ' temperature', & - ' entropy'/) - call write_file('td_quantities', 'thermodynamics', columns, thermodynamic_quantities, npart, ncols, num) - - unitnum = unitnum + 1 - deallocate(thermodynamic_quantities) - case(15) !Gravitational drag on sinks call gravitational_drag(time,npart,particlemass,xyzh,vxyzu) - case(16) call get_core_gas_com(time,npart,xyzh,vxyzu) - case(17) ncols = 6 allocate(columns(ncols)) From 4eb6a35ee11e340bae9ae6b6d10ac0b116d4502f Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Mon, 6 May 2024 10:58:40 +0200 Subject: [PATCH 7/8] (CE-analysis) remove redundant case 17, replaced by divv functionality --- src/utils/analysis_common_envelope.f90 | 54 -------------------------- 1 file changed, 54 deletions(-) diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index 0a40d006f..376048ce8 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -62,17 +62,6 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) integer :: unitnum,i,ncols logical :: requires_eos_opts - !case 7 variables - character(len=17), allocatable :: columns(:) - - !case 12 variables - real :: etoti, ekini, einti, epoti, phii - - real, dimension(3) :: com_xyz, com_vxyz - real, dimension(3) :: xyz_a, vxyz_a - real, allocatable :: histogram_data(:,:) - real :: ang_vel - !chose analysis type if (dump_number==0) then print "(40(a,/))", & @@ -90,7 +79,6 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) '13) MESA EoS compute total entropy and other average td quantities', & '15) Gravitational drag on sinks', & '16) CoM of gas around primary core', & - '17) Miscellaneous', & '18) J-E plane', & '19) Rotation profile', & '20) Energy profile', & @@ -202,48 +190,6 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) call gravitational_drag(time,npart,particlemass,xyzh,vxyzu) case(16) call get_core_gas_com(time,npart,xyzh,vxyzu) - case(17) - ncols = 6 - allocate(columns(ncols)) - columns = (/' x', & - ' y', & - ' z', & - ' r', & - 'spec. energy', & - ' omega ratio'/) - - call orbit_com(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptmass,com_xyz,com_vxyz) - - ang_vel = 0. - - do i=1,nptmass - if (xyzmh_ptmass(4,i) > 0.) then - 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) - ang_vel = ang_vel + (-xyz_a(2) * vxyz_a(1) + xyz_a(1) * vxyz_a(2)) / dot_product(xyz_a(1:2), xyz_a(1:2)) - endif - enddo - - ang_vel = ang_vel / 2. - - allocate(histogram_data(6,npart)) - - do i=1,npart - xyz_a(1:3) = xyzh(1:3,i) - com_xyz(1:3) - vxyz_a(1:3) = vxyzu(1:3,i) - com_vxyz(1:3) - - call calc_gas_energies(particlemass,poten(i),xyzh(:,i),vxyzu(:,i),radprop(:,i),xyzmh_ptmass,phii,epoti,ekini,einti,etoti) - histogram_data(1:3,i) = xyzh(1:3,i) - histogram_data(4,i) = distance(xyz_a(1:3)) - histogram_data(5,i) = epoti + ekini - histogram_data(6,i) = (-xyz_a(2) * vxyz_a(1) + xyz_a(1) * vxyz_a(2)) / dot_product(xyz_a(1:2), xyz_a(1:2)) - histogram_data(6,i) = (histogram_data(6,i) - ang_vel) / ang_vel - enddo - - call write_file('specific_energy_particles', 'histogram', columns, histogram_data, size(histogram_data(1,:)), ncols, num) - - deallocate(histogram_data) - case(18) call J_E_plane(num,npart,particlemass,xyzh,vxyzu) end select From 27c5e3508f6fd5c902cbfbd413991cc0f87d19c8 Mon Sep 17 00:00:00 2001 From: Mike Lau Date: Mon, 6 May 2024 11:06:51 +0200 Subject: [PATCH 8/8] (CE-analysis) clean up case numbers --- src/utils/analysis_common_envelope.f90 | 131 ++++++++++++------------- 1 file changed, 65 insertions(+), 66 deletions(-) diff --git a/src/utils/analysis_common_envelope.f90 b/src/utils/analysis_common_envelope.f90 index 376048ce8..a7ec86cff 100644 --- a/src/utils/analysis_common_envelope.f90 +++ b/src/utils/analysis_common_envelope.f90 @@ -59,7 +59,6 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) integer, intent(in) :: num,npart,iunit real, intent(inout) :: xyzh(:,:),vxyzu(:,:) real, intent(in) :: particlemass,time - integer :: unitnum,i,ncols logical :: requires_eos_opts !chose analysis type @@ -74,37 +73,37 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) ' 7) Simulation units and particle properties', & ' 8) Output .divv', & ' 9) EoS testing', & - '11) Profile of newly unbound particles', & - '12) Sink properties', & - '13) MESA EoS compute total entropy and other average td quantities', & - '15) Gravitational drag on sinks', & - '16) CoM of gas around primary core', & - '18) J-E plane', & - '19) Rotation profile', & - '20) Energy profile', & - '21) Recombination statistics', & - '22) Optical depth profile', & - '23) Particle tracker', & - '24) Unbound ion fraction', & - '25) Optical depth at recombination', & - '26) Envelope binding energy', & - '27) Print dumps number matching separation', & - '28) Companion mass coordinate vs. time', & - '29) Energy histogram',& - '30) Analyse disk',& - '31) Recombination energy vs time',& - '32) Binding energy profile',& - '33) planet_rvm',& - '34) Velocity histogram',& - '35) Unbound temperature',& - '36) Planet mass distribution',& - '37) Planet profile',& - '38) Velocity profile',& - '39) Angular momentum profile',& - '40) Keplerian velocity profile',& - '41) Total dust mass' + '10) Profile of newly unbound particles', & + '11) Sink properties', & + '12) MESA EoS compute total entropy and other average td quantities', & + '13) Gravitational drag on sinks', & + '14) CoM of gas around primary core', & + '15) J-E plane', & + '16) Rotation profile', & + '17) Energy profile', & + '18) Recombination statistics', & + '19) Optical depth profile', & + '20) Particle tracker', & + '21) Unbound ion fraction', & + '22) Optical depth at recombination', & + '23) Envelope binding energy', & + '24) Print dumps number matching separation', & + '25) Companion mass coordinate vs. time', & + '26) Energy histogram',& + '27) Analyse disk',& + '28) Recombination energy vs time',& + '29) Binding energy profile',& + '30) planet_rvm',& + '31) Velocity histogram',& + '32) Unbound temperature',& + '33) Planet mass distribution',& + '34) Planet profile',& + '35) Velocity profile',& + '36) Angular momentum profile',& + '37) Keplerian velocity profile',& + '38) Total dust mass' analysis_to_perform = 1 - call prompt('Choose analysis type ',analysis_to_perform,1,41) + call prompt('Choose analysis type ',analysis_to_perform,1,38) endif call reset_centreofmass(npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptmass) @@ -112,7 +111,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) xyzmh_ptmass,vxyz_ptmass,omega_corotate,dump_number) ! List of analysis options that require specifying EOS options - requires_eos_opts = any((/2,3,4,5,6,8,9,11,13,15,20,21,22,23,24,25,26,29,30,31,32,33,35,41/) == analysis_to_perform) + requires_eos_opts = any((/2,3,4,5,6,8,9,10,13,17,18,19,20,21,22,23,26,27,28,29,30,32,38/) == analysis_to_perform) if (dump_number == 0 .and. requires_eos_opts) call set_eos_options(analysis_to_perform) select case(analysis_to_perform) @@ -134,64 +133,64 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) call output_divv_files(time,dumpfile,npart,particlemass,xyzh,vxyzu) case(9) !EoS testing call eos_surfaces - case(11) !New unbound particle profiles in time + case(10) !New unbound particle profiles in time call unbound_profiles(time,num,npart,particlemass,xyzh,vxyzu) - case(19) ! Rotation profile + case(11) !sink properties + call sink_properties(time,npart,particlemass,xyzh,vxyzu) + case(12) !MESA EoS compute total entropy and other average thermodynamical quantities + call bound_unbound_thermo(time,npart,particlemass,xyzh,vxyzu) + case(13) !Gravitational drag on sinks + call gravitational_drag(time,npart,particlemass,xyzh,vxyzu) + case(14) + call get_core_gas_com(time,npart,xyzh,vxyzu) + case(15) + call J_E_plane(num,npart,particlemass,xyzh,vxyzu) + case(16) ! Rotation profile call rotation_profile(time,num,npart,xyzh,vxyzu) - case(20) ! Energy profile + case(17) ! Energy profile call energy_profile(time,npart,particlemass,xyzh,vxyzu) - case(21) ! Recombination statistics + case(18) ! Recombination statistics call recombination_stats(time,num,npart,particlemass,xyzh,vxyzu) - case(22) ! Optical depth profile + case(19) ! Optical depth profile call tau_profile(time,num,npart,particlemass,xyzh) - case(23) ! Particle tracker + case(20) ! Particle tracker call track_particle(time,particlemass,xyzh,vxyzu) - case(24) ! Unbound ion fractions + case(21) ! Unbound ion fractions call unbound_ionfrac(time,npart,particlemass,xyzh,vxyzu) - case(25) ! Optical depth at recombination + case(22) ! Optical depth at recombination call recombination_tau(time,npart,particlemass,xyzh,vxyzu) - case(26) ! Calculate binding energy outside core + case(23) ! Calculate binding energy outside core call env_binding_ene(npart,particlemass,xyzh,vxyzu) - case(27) ! Print dump number corresponding to given set of sink-sink separations + case(24) ! Print dump number corresponding to given set of sink-sink separations call print_dump_numbers(dumpfile) - case(28) ! Companion mass coordinate (spherical mass shells) vs. time + case(25) ! Companion mass coordinate (spherical mass shells) vs. time call m_vs_t(time,npart,particlemass,xyzh) - case(29) ! Energy histogram + case(26) ! Energy histogram call energy_hist(time,npart,particlemass,xyzh,vxyzu) - case(30) ! Analyse disk around companion + case(27) ! Analyse disk around companion call analyse_disk(num,npart,particlemass,xyzh,vxyzu) - case(31) ! Recombination energy vs. time + case(28) ! Recombination energy vs. time call erec_vs_t(time,npart,particlemass,xyzh) - case(32) ! Binding energy profile + case(29) ! Binding energy profile call create_bindingEnergy_profile(time,num,npart,particlemass,xyzh,vxyzu) - case(33) ! Planet coordinates and mass + case(30) ! Planet coordinates and mass call planet_rvm(time,particlemass,xyzh,vxyzu) - case(34) ! Velocity histogram + case(31) ! Velocity histogram call velocity_histogram(time,num,npart,particlemass,xyzh,vxyzu) - case(35) ! Unbound temperatures + case(32) ! Unbound temperatures call unbound_temp(time,npart,particlemass,xyzh,vxyzu) - case(36) ! Planet mass distribution + case(33) ! Planet mass distribution call planet_mass_distribution(time,num,npart,xyzh) - case(37) ! Calculate planet profile + case(34) ! Calculate planet profile call planet_profile(num,dumpfile,particlemass,xyzh,vxyzu) - case(38) ! Velocity profile + case(35) ! Velocity profile call velocity_profile(time,num,npart,particlemass,xyzh,vxyzu) - case(39) ! Angular momentum profile + case(36) ! Angular momentum profile call angular_momentum_profile(time,num,npart,particlemass,xyzh,vxyzu) - case(40) ! Keplerian velocity profile + case(37) ! Keplerian velocity profile call vkep_profile(time,num,npart,particlemass,xyzh,vxyzu) - case(41) !Total dust mass + case(38) !Total dust mass call total_dust_mass(time,npart,particlemass,xyzh) - case(12) !sink properties - call sink_properties(time,npart,particlemass,xyzh,vxyzu) - case(13) !MESA EoS compute total entropy and other average thermodynamical quantities - call bound_unbound_thermo(time,npart,particlemass,xyzh,vxyzu) - case(15) !Gravitational drag on sinks - call gravitational_drag(time,npart,particlemass,xyzh,vxyzu) - case(16) - call get_core_gas_com(time,npart,xyzh,vxyzu) - case(18) - call J_E_plane(num,npart,particlemass,xyzh,vxyzu) end select !increase dump number counter dump_number = dump_number + 1