From 549dbc407d26ebf81ae40008b2a7d68e50d500ca Mon Sep 17 00:00:00 2001 From: Stephen Neilson <36410751+s-neilson@users.noreply.github.com> Date: Wed, 28 Sep 2022 00:32:53 +1000 Subject: [PATCH 1/3] Made compiling with gfortran on aarch64 processors possible. --- build/Makefile | 1 + build/Makefile_systems | 4 ++++ 2 files changed, 5 insertions(+) diff --git a/build/Makefile b/build/Makefile index 6088790ae..babdc6d97 100644 --- a/build/Makefile +++ b/build/Makefile @@ -31,6 +31,7 @@ SHELL = /bin/bash VPATH = "${RUNDIR}" ../src/main ../src/utils ../src/setup ../src/tests ../src/lib/NICIL/src BINDIR= ../bin UNAME=${shell uname} +UNAMEP=${shell uname -p} #---------------------------------------------------------------- # Sensible defaults for phantom configuration #---------------------------------------------------------------- diff --git a/build/Makefile_systems b/build/Makefile_systems index b380d6a43..71c5aa56d 100644 --- a/build/Makefile_systems +++ b/build/Makefile_systems @@ -140,7 +140,11 @@ endif ifeq ($(SYSTEM), gfortran) include Makefile_defaults_gfortran ifneq ($(UNAME), Darwin) +ifeq ($(UNAMEP), aarch64) # -mcmodel=medium is not an option for aarch64 type processors with gfortran. + FFLAGS+= -mcmodel=small +else FFLAGS+= -mcmodel=medium +endif endif endif From 3d0e903dbbaf2c682d94b79ca4ec8bc3fc230bc2 Mon Sep 17 00:00:00 2001 From: Stephen Neilson <36410751+s-neilson@users.noreply.github.com> Date: Fri, 14 Oct 2022 13:48:29 +1100 Subject: [PATCH 2/3] Added options for planetary destruction related information and binding energy profile in analysis_common_envelope.F90 --- src/utils/analysis_common_envelope.F90 | 204 ++++++++++++++++++++++++- 1 file changed, 201 insertions(+), 3 deletions(-) diff --git a/src/utils/analysis_common_envelope.F90 b/src/utils/analysis_common_envelope.F90 index 1c5d5c01e..7ec0b7268 100644 --- a/src/utils/analysis_common_envelope.F90 +++ b/src/utils/analysis_common_envelope.F90 @@ -86,7 +86,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) !chose analysis type if (dump_number==0) then - print "(29(a,/))", & + print "(31(a,/))", & ' 1) Sink separation', & ' 2) Bound and unbound quantities', & ' 3) Energies', & @@ -96,6 +96,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) ' 7) Simulation units and particle properties', & ' 8) Output .divv', & ' 9) EoS testing', & + '10) Planet destruction', & '11) Profile of newly unbound particles', & '12) Sink properties', & '13) MESA EoS compute total entropy and other average td quantities', & @@ -115,11 +116,12 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) '27) Print dumps number matching separation', & '28) Companion mass coordinate vs. time', & '29) Energy histogram',& - '30) Analyse disk' + '30) Analyse disk', & + '31) Binding energy profile' analysis_to_perform = 1 - call prompt('Choose analysis type ',analysis_to_perform,1,30) + call prompt('Choose analysis type ',analysis_to_perform,1,31) endif @@ -171,6 +173,8 @@ 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(10) !Planet destruction + call planet_destruction(time,npart,particlemass,xyzh,vxyzu) case(11) !New unbound particle profiles in time call unbound_profiles(time,num,npart,particlemass,xyzh,vxyzu) case(19) ! Rotation profile @@ -197,6 +201,8 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) call energy_hist(time,npart,particlemass,xyzh,vxyzu) case(30) ! Analyse disk around companion call analyse_disk(num,npart,particlemass,xyzh,vxyzu) + case(31) !Binding energy profile + call create_bindingEnergy_profile(time,num,npart,particlemass,xyzh,vxyzu) case(12) !sink properties call sink_properties(time,npart,particlemass,xyzh,vxyzu) case(13) !MESA EoS compute total entropy and other average thermodynamical quantities @@ -2639,6 +2645,151 @@ subroutine J_E_plane(num,npart,particlemass,xyzh,vxyzu) end subroutine J_E_plane +!------------------------------------------------------------------- +!+ +! Planet destruction +!+ +!------------------------------------------------------------------- +subroutine planet_destruction(time,npart,particlemass,xyzh,vxyzu) + use kernel, only:wkern + integer, intent(in) :: npart + real, intent(in) :: time,particlemass + real, intent(inout) :: xyzh(:,:),vxyzu(:,:) + character(len=17), allocatable :: columns(:) + character(len=18) :: filename + real, allocatable :: planetDestruction(:) + integer :: ncols,i,j + real, allocatable, save :: time_old + real, allocatable, save :: particleRho(:) + character(len=50) :: planetRadiusPromptString + real, allocatable, save :: planetRadii(:) !In units of Rsun + + real, dimension(3) :: currentGasVel, currentVelContrast + real :: currentRho(1) !Is a one element array because sphInterpolation returns a 1 dimensional array. + real :: currentRhoScaled,currentVelContrastScaled,currentPlanetRhoScaled + real :: currentPlanetMassScaled,currentPlanetRadiusScaled + real, allocatable, save :: currentKhAblatedMass(:) + + ncols=5 + allocate(columns(ncols)) + allocate(planetDestruction(ncols)) + columns=(/" rhoGas", & + " kh_rhoCrit", & + " kh_lmax", & + " kh_mdot", & + " kh_ablatedM" /) + + !Kelvin-Helmholtz instability planet destruction as described in "On the survival of brown dwarfs + !and planets by their giant host star" (https://arxiv.org/abs/1210.0879). Description of columns: + !rhoGas: background gas density at sink. In units of g/cm^3. + !kh_rhoCrit: paper equation 5. In units of g/cm^3. + !kh_lmax: paper equation 6. In units of Jupiter radii. + !kh_mdot: paper equation 7. In units of Jupiter mass/year. + !kh_ablatedM: kh_mdot integrated over time. In units of Jupiter masses. + + do i=1,nptmass + if (i==1) cycle !The first sink is assumed to be the core. + + if ((dump_number==0) .and. (i==2)) then !This is only done once. + allocate(planetRadii(nptmass)) + planetRadii=0.1 + do j=2,nptmass + write(planetRadiusPromptString,"(A13,I0,A32)") "Enter planet ",j-1," radius in units of solar radii" + call prompt(planetRadiusPromptString,planetRadii(i),0.0,1.0) + enddo + + allocate(time_old) + allocate(particleRho(npart)) + allocate(currentKhAblatedMass(nptmass)) + + time_old=0.0 + particleRho=getParticleRho(xyzh(4,:),particlemass) + currentKhAblatedMass=0.0 + endif + + + currentRho=sphInterpolation(npart,particlemass,particleRho,xyzh,xyzmh_ptmass(1:3,i),reshape(particleRho,(/1,npart/))) + currentGasVel=sphInterpolation(npart,particlemass,particleRho,xyzh,xyzmh_ptmass(1:3,i),vxyzu(1:3,:)) + currentVelContrast=vxyz_ptmass(1:3,i)-currentGasVel + + currentPlanetRadiusScaled=planetRadii(i)/0.1 !In units of 0.1 Rsun. + currentPlanetMassScaled=xyzmh_ptmass(4,i)*104.74 !In units of 10 jupiter masses. + currentPlanetRhoScaled=(xyzmh_ptmass(4,i)/((4.0/3.0)*pi*(planetRadii(i)**3.0)))*0.44 !In units of 13.34 g/cm^3 + currentRhoScaled=currentRho(1)*59000.0 !In units of 10^-4 g/cm^3. + currentVelContrastScaled=distance(currentVelContrast)*4.37 !In units of 100 km/s. + + planetDestruction(1)=currentRho(1)*5.9 + planetDestruction(2)=3.82*(currentPlanetRhoScaled**(4.0/3.0))*(currentPlanetMassScaled**(2.0/3.0))& + *(currentVelContrastScaled**(-2.0)) + planetDestruction(3)=0.0000263*(currentVelContrastScaled**2.0)*currentRhoScaled*(currentPlanetRhoScaled**((-5.0)/3.0))& + *(currentPlanetMassScaled**((-1.0)/3.0)) + planetDestruction(4)=11.0*currentVelContrastScaled*currentRhoScaled*(currentPlanetRadiusScaled**2.0)& + *(planetDestruction(3)/(currentPlanetRadiusScaled*0.973)) + + currentKhAblatedMass(i)=currentKhAblatedMass(i)+((time-time_old)*planetDestruction(4)*0.0000505) + planetDestruction(5)=currentKhAblatedMass(i) + + + write(filename, "(A17,I0)") "sink_destruction_",i + call write_time_file(filename, columns, time, planetDestruction, ncols, dump_number) + enddo + + time_old=time +end subroutine planet_destruction + +!----------------------------------------------------------------------------------------- +!+ +!Binding energy profile +!+ +!----------------------------------------------------------------------------------------- +subroutine create_bindingEnergy_profile(time,num,npart,particlemass,xyzh,vxyzu) + real, intent(in) :: time,particlemass + integer, intent(in) :: num,npart + real, intent(in) :: xyzh(4,npart),vxyzu(4,npart) + + character(len=17), allocatable :: columns(:) + real, allocatable :: profile(:,:) + integer :: ncols,i,j,iorder(npart) + real :: currentInteriorMass,currentParticleGPE,currentCoreParticleSeparation + real :: previousBindingEnergy,previousBindingEnergyU + + ncols=3 + allocate(columns(ncols)) + allocate(profile(ncols,npart)) + columns=(/" radius",& + " bEnergy",& !Binding energy without internal energy. + " bEnergy (u)"/) !Binding energy with internal energy. + + + call set_r2func_origin(xyzmh_ptmass(1,1),xyzmh_ptmass(2,1),xyzmh_ptmass(3,1)) + call indexxfunc(npart,r2func_origin,xyzh,iorder) + currentInteriorMass=xyzmh_ptmass(4,1)+(npart*particlemass) !Initally set to the entire mass of the star. + + do i=npart,1,-1 !Loops over all particles from outer to inner. + j=iorder(i) + currentInteriorMass=currentInteriorMass-particlemass + currentCoreParticleSeparation=separation(xyzmh_ptmass(1:3,1),xyzh(1:3,j)) + currentParticleGPE=(currentInteriorMass*particlemass)/currentCoreParticleSeparation + + !The binding energy at a particular radius is the sum of the gravitational potential energies + !(and internal energies in the case of the third column) of all particles beyond that radius. + if (i==npart) then + previousBindingEnergy=0.0 + previousBindingEnergyU=0.0 + else + previousBindingEnergy=profile(2,i+1) + previousBindingEnergyU=profile(3,i+1) + endif + + profile(1,i)=currentCoreParticleSeparation + profile(2,i)=previousBindingEnergy+currentParticleGPE + profile(3,i)=previousBindingEnergyU+currentParticleGPE-(vxyzu(4,j)*particlemass) + enddo + + + call write_file('bEnergyProfile','bEnergyProfiles',columns,profile,npart,ncols,num) +end subroutine create_bindingEnergy_profile + subroutine get_core_gas_com(time,npart,xyzh,vxyzu) use sortutils, only:set_r2func_origin,r2func_origin,indexxfunc @@ -3459,6 +3610,53 @@ real function separation(a,b) separation = distance(a - b) end function separation +!Creates an array of SPH particle densities for each value of h. +elemental real function getParticleRho(h,particlemass) + real, intent(in) :: h,particlemass + getParticleRho=rhoh(h,particlemass) +end function getParticleRho + +!Performs SPH interpolation on the SPH particle property toInterpolate at the location interpolateXyz. +!The smoothing length used is the smoothing length of the closest SPH particle to interpolateXyz. +function sphInterpolation(npart,particlemass,particleRho,particleXyzh,interpolateXyz,toInterpolate) result(interpolatedData) + use kernel, only:wkern + integer, intent(in) :: npart + real, intent(in) :: particlemass + real, intent(in) :: particleRho(npart) + real, intent(in) :: particleXyzh(4,npart) + real, intent(in) :: interpolateXyz(3) + real, intent(in) :: toInterpolate(:,:) + real :: interpolatedData(size(toInterpolate,1)) + + integer :: i,j,iorder(npart) + real :: currentR,currentQ,currentQ2 + real :: nearestSphH + real :: currentParticleRho,currentSphSummandFactor + + interpolatedData=0.0 + call set_r2func_origin(interpolateXyz(1),interpolateXyz(2),interpolateXyz(3)) + call indexxfunc(npart,r2func_origin,particleXyzh,iorder) !Gets the order of SPH particles from the interpolation point. + nearestSphH=particleXyzh(4,iorder(1)) !The smoothing length of the nearest SPH particle to the ineterpolation point. + + do i=1,npart + j=iorder(i) + + currentR=separation(interpolateXyz,particleXyzh(1:3,j)) + currentQ=currentR/nearestSphH !currentR is scaled in units of nearestSphH + currentQ2=currentQ**2.0 + + !All SPH particles beyond 2 smoothing lengths are ignored. + if (currentQ>2) then + exit + endif + + !SPH interpolation is done below. + currentParticleRho=particleRho(j) + currentSphSummandFactor=(particlemass/currentParticleRho)*((1.0/((nearestSphH**3.0)*pi))*wkern(currentQ2,currentQ)) + interpolatedData=interpolatedData+(currentSphSummandFactor*toInterpolate(:,j)) + enddo +end function sphInterpolation + !Sorting routines recursive subroutine quicksort(a, first, last, ncols, sortcol) integer, intent(in) :: first, last, ncols, sortcol From ebbdaf53fcc4ba0297f59b5bb1fcfec97979c75b Mon Sep 17 00:00:00 2001 From: s-neilson <36410751+s-neilson@users.noreply.github.com> Date: Tue, 18 Oct 2022 01:36:23 +1100 Subject: [PATCH 3/3] Revert "Added options for planetary destruction related information and binding energy profile in analysis_common_envelope.F90" in order to put it in a separate branch. This reverts commit 3d0e903dbbaf2c682d94b79ca4ec8bc3fc230bc2. --- src/utils/analysis_common_envelope.F90 | 204 +------------------------ 1 file changed, 3 insertions(+), 201 deletions(-) diff --git a/src/utils/analysis_common_envelope.F90 b/src/utils/analysis_common_envelope.F90 index 7ec0b7268..1c5d5c01e 100644 --- a/src/utils/analysis_common_envelope.F90 +++ b/src/utils/analysis_common_envelope.F90 @@ -86,7 +86,7 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) !chose analysis type if (dump_number==0) then - print "(31(a,/))", & + print "(29(a,/))", & ' 1) Sink separation', & ' 2) Bound and unbound quantities', & ' 3) Energies', & @@ -96,7 +96,6 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) ' 7) Simulation units and particle properties', & ' 8) Output .divv', & ' 9) EoS testing', & - '10) Planet destruction', & '11) Profile of newly unbound particles', & '12) Sink properties', & '13) MESA EoS compute total entropy and other average td quantities', & @@ -116,12 +115,11 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) '27) Print dumps number matching separation', & '28) Companion mass coordinate vs. time', & '29) Energy histogram',& - '30) Analyse disk', & - '31) Binding energy profile' + '30) Analyse disk' analysis_to_perform = 1 - call prompt('Choose analysis type ',analysis_to_perform,1,31) + call prompt('Choose analysis type ',analysis_to_perform,1,30) endif @@ -173,8 +171,6 @@ 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(10) !Planet destruction - call planet_destruction(time,npart,particlemass,xyzh,vxyzu) case(11) !New unbound particle profiles in time call unbound_profiles(time,num,npart,particlemass,xyzh,vxyzu) case(19) ! Rotation profile @@ -201,8 +197,6 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit) call energy_hist(time,npart,particlemass,xyzh,vxyzu) case(30) ! Analyse disk around companion call analyse_disk(num,npart,particlemass,xyzh,vxyzu) - case(31) !Binding energy profile - call create_bindingEnergy_profile(time,num,npart,particlemass,xyzh,vxyzu) case(12) !sink properties call sink_properties(time,npart,particlemass,xyzh,vxyzu) case(13) !MESA EoS compute total entropy and other average thermodynamical quantities @@ -2645,151 +2639,6 @@ subroutine J_E_plane(num,npart,particlemass,xyzh,vxyzu) end subroutine J_E_plane -!------------------------------------------------------------------- -!+ -! Planet destruction -!+ -!------------------------------------------------------------------- -subroutine planet_destruction(time,npart,particlemass,xyzh,vxyzu) - use kernel, only:wkern - integer, intent(in) :: npart - real, intent(in) :: time,particlemass - real, intent(inout) :: xyzh(:,:),vxyzu(:,:) - character(len=17), allocatable :: columns(:) - character(len=18) :: filename - real, allocatable :: planetDestruction(:) - integer :: ncols,i,j - real, allocatable, save :: time_old - real, allocatable, save :: particleRho(:) - character(len=50) :: planetRadiusPromptString - real, allocatable, save :: planetRadii(:) !In units of Rsun - - real, dimension(3) :: currentGasVel, currentVelContrast - real :: currentRho(1) !Is a one element array because sphInterpolation returns a 1 dimensional array. - real :: currentRhoScaled,currentVelContrastScaled,currentPlanetRhoScaled - real :: currentPlanetMassScaled,currentPlanetRadiusScaled - real, allocatable, save :: currentKhAblatedMass(:) - - ncols=5 - allocate(columns(ncols)) - allocate(planetDestruction(ncols)) - columns=(/" rhoGas", & - " kh_rhoCrit", & - " kh_lmax", & - " kh_mdot", & - " kh_ablatedM" /) - - !Kelvin-Helmholtz instability planet destruction as described in "On the survival of brown dwarfs - !and planets by their giant host star" (https://arxiv.org/abs/1210.0879). Description of columns: - !rhoGas: background gas density at sink. In units of g/cm^3. - !kh_rhoCrit: paper equation 5. In units of g/cm^3. - !kh_lmax: paper equation 6. In units of Jupiter radii. - !kh_mdot: paper equation 7. In units of Jupiter mass/year. - !kh_ablatedM: kh_mdot integrated over time. In units of Jupiter masses. - - do i=1,nptmass - if (i==1) cycle !The first sink is assumed to be the core. - - if ((dump_number==0) .and. (i==2)) then !This is only done once. - allocate(planetRadii(nptmass)) - planetRadii=0.1 - do j=2,nptmass - write(planetRadiusPromptString,"(A13,I0,A32)") "Enter planet ",j-1," radius in units of solar radii" - call prompt(planetRadiusPromptString,planetRadii(i),0.0,1.0) - enddo - - allocate(time_old) - allocate(particleRho(npart)) - allocate(currentKhAblatedMass(nptmass)) - - time_old=0.0 - particleRho=getParticleRho(xyzh(4,:),particlemass) - currentKhAblatedMass=0.0 - endif - - - currentRho=sphInterpolation(npart,particlemass,particleRho,xyzh,xyzmh_ptmass(1:3,i),reshape(particleRho,(/1,npart/))) - currentGasVel=sphInterpolation(npart,particlemass,particleRho,xyzh,xyzmh_ptmass(1:3,i),vxyzu(1:3,:)) - currentVelContrast=vxyz_ptmass(1:3,i)-currentGasVel - - currentPlanetRadiusScaled=planetRadii(i)/0.1 !In units of 0.1 Rsun. - currentPlanetMassScaled=xyzmh_ptmass(4,i)*104.74 !In units of 10 jupiter masses. - currentPlanetRhoScaled=(xyzmh_ptmass(4,i)/((4.0/3.0)*pi*(planetRadii(i)**3.0)))*0.44 !In units of 13.34 g/cm^3 - currentRhoScaled=currentRho(1)*59000.0 !In units of 10^-4 g/cm^3. - currentVelContrastScaled=distance(currentVelContrast)*4.37 !In units of 100 km/s. - - planetDestruction(1)=currentRho(1)*5.9 - planetDestruction(2)=3.82*(currentPlanetRhoScaled**(4.0/3.0))*(currentPlanetMassScaled**(2.0/3.0))& - *(currentVelContrastScaled**(-2.0)) - planetDestruction(3)=0.0000263*(currentVelContrastScaled**2.0)*currentRhoScaled*(currentPlanetRhoScaled**((-5.0)/3.0))& - *(currentPlanetMassScaled**((-1.0)/3.0)) - planetDestruction(4)=11.0*currentVelContrastScaled*currentRhoScaled*(currentPlanetRadiusScaled**2.0)& - *(planetDestruction(3)/(currentPlanetRadiusScaled*0.973)) - - currentKhAblatedMass(i)=currentKhAblatedMass(i)+((time-time_old)*planetDestruction(4)*0.0000505) - planetDestruction(5)=currentKhAblatedMass(i) - - - write(filename, "(A17,I0)") "sink_destruction_",i - call write_time_file(filename, columns, time, planetDestruction, ncols, dump_number) - enddo - - time_old=time -end subroutine planet_destruction - -!----------------------------------------------------------------------------------------- -!+ -!Binding energy profile -!+ -!----------------------------------------------------------------------------------------- -subroutine create_bindingEnergy_profile(time,num,npart,particlemass,xyzh,vxyzu) - real, intent(in) :: time,particlemass - integer, intent(in) :: num,npart - real, intent(in) :: xyzh(4,npart),vxyzu(4,npart) - - character(len=17), allocatable :: columns(:) - real, allocatable :: profile(:,:) - integer :: ncols,i,j,iorder(npart) - real :: currentInteriorMass,currentParticleGPE,currentCoreParticleSeparation - real :: previousBindingEnergy,previousBindingEnergyU - - ncols=3 - allocate(columns(ncols)) - allocate(profile(ncols,npart)) - columns=(/" radius",& - " bEnergy",& !Binding energy without internal energy. - " bEnergy (u)"/) !Binding energy with internal energy. - - - call set_r2func_origin(xyzmh_ptmass(1,1),xyzmh_ptmass(2,1),xyzmh_ptmass(3,1)) - call indexxfunc(npart,r2func_origin,xyzh,iorder) - currentInteriorMass=xyzmh_ptmass(4,1)+(npart*particlemass) !Initally set to the entire mass of the star. - - do i=npart,1,-1 !Loops over all particles from outer to inner. - j=iorder(i) - currentInteriorMass=currentInteriorMass-particlemass - currentCoreParticleSeparation=separation(xyzmh_ptmass(1:3,1),xyzh(1:3,j)) - currentParticleGPE=(currentInteriorMass*particlemass)/currentCoreParticleSeparation - - !The binding energy at a particular radius is the sum of the gravitational potential energies - !(and internal energies in the case of the third column) of all particles beyond that radius. - if (i==npart) then - previousBindingEnergy=0.0 - previousBindingEnergyU=0.0 - else - previousBindingEnergy=profile(2,i+1) - previousBindingEnergyU=profile(3,i+1) - endif - - profile(1,i)=currentCoreParticleSeparation - profile(2,i)=previousBindingEnergy+currentParticleGPE - profile(3,i)=previousBindingEnergyU+currentParticleGPE-(vxyzu(4,j)*particlemass) - enddo - - - call write_file('bEnergyProfile','bEnergyProfiles',columns,profile,npart,ncols,num) -end subroutine create_bindingEnergy_profile - subroutine get_core_gas_com(time,npart,xyzh,vxyzu) use sortutils, only:set_r2func_origin,r2func_origin,indexxfunc @@ -3610,53 +3459,6 @@ real function separation(a,b) separation = distance(a - b) end function separation -!Creates an array of SPH particle densities for each value of h. -elemental real function getParticleRho(h,particlemass) - real, intent(in) :: h,particlemass - getParticleRho=rhoh(h,particlemass) -end function getParticleRho - -!Performs SPH interpolation on the SPH particle property toInterpolate at the location interpolateXyz. -!The smoothing length used is the smoothing length of the closest SPH particle to interpolateXyz. -function sphInterpolation(npart,particlemass,particleRho,particleXyzh,interpolateXyz,toInterpolate) result(interpolatedData) - use kernel, only:wkern - integer, intent(in) :: npart - real, intent(in) :: particlemass - real, intent(in) :: particleRho(npart) - real, intent(in) :: particleXyzh(4,npart) - real, intent(in) :: interpolateXyz(3) - real, intent(in) :: toInterpolate(:,:) - real :: interpolatedData(size(toInterpolate,1)) - - integer :: i,j,iorder(npart) - real :: currentR,currentQ,currentQ2 - real :: nearestSphH - real :: currentParticleRho,currentSphSummandFactor - - interpolatedData=0.0 - call set_r2func_origin(interpolateXyz(1),interpolateXyz(2),interpolateXyz(3)) - call indexxfunc(npart,r2func_origin,particleXyzh,iorder) !Gets the order of SPH particles from the interpolation point. - nearestSphH=particleXyzh(4,iorder(1)) !The smoothing length of the nearest SPH particle to the ineterpolation point. - - do i=1,npart - j=iorder(i) - - currentR=separation(interpolateXyz,particleXyzh(1:3,j)) - currentQ=currentR/nearestSphH !currentR is scaled in units of nearestSphH - currentQ2=currentQ**2.0 - - !All SPH particles beyond 2 smoothing lengths are ignored. - if (currentQ>2) then - exit - endif - - !SPH interpolation is done below. - currentParticleRho=particleRho(j) - currentSphSummandFactor=(particlemass/currentParticleRho)*((1.0/((nearestSphH**3.0)*pi))*wkern(currentQ2,currentQ)) - interpolatedData=interpolatedData+(currentSphSummandFactor*toInterpolate(:,j)) - enddo -end function sphInterpolation - !Sorting routines recursive subroutine quicksort(a, first, last, ncols, sortcol) integer, intent(in) :: first, last, ncols, sortcol