Skip to content

Fmm prec plot #786

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
Yrisch wants to merge 10 commits into
base: master
Choose a base branch
from
Draft

Fmm prec plot #786

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
69fe25d
(test_grav) long_range force comparison
Yrisch Dec 10, 2025
8c457b2
(test_grav) agree with shamrock
Yrisch Dec 12, 2025
80ff94a
(test_grav) Plummer sphere before polish
Yrisch Dec 18, 2025
8478a88
(test_grav) bad clean up
Yrisch Dec 19, 2025
e458621
(grav) final clean up, FMM method switch added
Yrisch Jan 2, 2026
138eed0
Merge branch 'master' into FMM_prec_plot
Yrisch Jan 15, 2026
065173e
(test-grav) final cleanups before PR
Yrisch Jan 15, 2026
fd63ae3
(test_grav) fix mase test
Yrisch Jan 15, 2026
e560232
(test_grav) add plotplummer to grav test strings
Yrisch Jan 21, 2026
4f96a88
(grav) minor fixes and adjust toli
Yrisch Jan 27, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 0 additions & 1 deletion src/main/kdtree.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,6 @@ module kdtree
public :: maketreeglobal
public :: empty_tree
public :: compute_M2L,expand_fgrav_in_taylor_series

integer, public :: maxlevel_indexed, maxlevel

type kdbuildstack
Expand Down
25 changes: 17 additions & 8 deletions src/main/neigh_kdtree.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ module neighkdtree
integer(kind=8), public :: ncells
real, public :: dxcell
real, public :: dcellx = 0.,dcelly = 0.,dcellz = 0.
logical, public :: force_dual_walk
integer, public :: neigh_switch = 1
integer :: globallevel,refinelevels

public :: allocate_neigh, deallocate_neigh
Expand Down Expand Up @@ -279,13 +279,19 @@ subroutine get_neighbour_list(inode,mylistneigh,nneigh,xyzh,xyzcache,ixyzcachesi
endif

! Find neighbours of this cell on this node
if ((get_f .or. force_dual_walk) .and. (.not.mpi)) then
call getneigh_dual(node,xpos,xsizei,rcuti,mylistneigh,nneigh,xyzcache,ixyzcachesize,&
leaf_is_active,get_j,get_f,fgrav,inode)
else
select case(neigh_switch)
case(1)
if (get_f .and. .not.(mpi)) then
call getneigh_dual(node,xpos,xsizei,rcuti,mylistneigh,nneigh,xyzcache,ixyzcachesize,&
leaf_is_active,get_j,get_f,fgrav,inode)
else
call getneigh(node,xpos,xsizei,rcuti,mylistneigh,nneigh,xyzcache,ixyzcachesize,&
leaf_is_active,get_j,get_f,fgrav)
endif
case(2)
call getneigh(node,xpos,xsizei,rcuti,mylistneigh,nneigh,xyzcache,ixyzcachesize,&
leaf_is_active,get_j,get_f,fgrav)
endif
end select

if (get_f) f = fgrav + fgrav_global

Expand Down Expand Up @@ -328,6 +334,7 @@ subroutine write_options_tree(iunit)

if (gravity) then
call write_inopt(tree_accuracy,'tree_accuracy','tree opening criterion (0.0-1.0)',iunit)
!call write_inopt(neigh_switch,'neigh_switch','SFMM=1, FMM=2 (STAY with SFMM!!)',iunit)
endif

end subroutine write_options_tree
Expand All @@ -344,8 +351,10 @@ subroutine read_options_tree(db,nerr)
type(inopts), intent(inout) :: db(:)
integer, intent(inout) :: nerr

if (gravity) call read_inopt(tree_accuracy,'tree_accuracy',db,errcount=nerr,min=0.,max=1.)

if (gravity) then
call read_inopt(tree_accuracy,'tree_accuracy',db,errcount=nerr,min=0.,max=1.)
call read_inopt(neigh_switch,'neigh_switch',db,errcount=nerr,min=1,max=2)
endif
end subroutine read_options_tree

!-----------------------------------------------------------------------
Expand Down
11 changes: 6 additions & 5 deletions src/setup/set_sphere.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,8 @@ module spherical

public :: set_sphere,set_ellipse,rho_func

integer,public :: iseed_mc = -1978

integer, parameter :: &
ierr_notinrange = 1, &
ierr_not_converged = 2, &
Expand Down Expand Up @@ -162,7 +164,7 @@ subroutine set_sphere_mc(id,master,rmin,rmax,hfact,np_requested,np,xyzh, &
integer(kind=8), intent(inout) :: nptot
logical, intent(in) :: verbose
procedure(mask_prototype) :: mask
integer :: i,npin,iseed,maxp
integer :: i,npin,maxp
real :: vol_sphere,rr,phi,theta,mr,dir(3)
real :: sintheta,costheta,sinphi,cosphi,psep
integer(kind=8) :: iparttot
Expand All @@ -172,24 +174,23 @@ subroutine set_sphere_mc(id,master,rmin,rmax,hfact,np_requested,np,xyzh, &
vol_sphere = 4./3.*pi*(rmax**3 - rmin**3)
! use mean particle spacing to set initial smoothing lengths
psep = (vol_sphere/real(np_requested))**(1./3.)
iseed = -1978
maxp = size(xyzh(1,:))
ierr = 1

do i=npin+1,npin+np_requested,2
!
! get random mass coordinate i.e. m(r)
!
mr = ran2(iseed)
mr = ran2(iseed_mc)
!
! invert to get mass coordinate from radial coordinate, i.e. r(m)
!
rr = (rmin**3+mr*(rmax**3-rmin**3))**(1./3.) ! uniform density
!
! get a random position on sphere
!
phi = 2.*pi*(ran2(iseed) - 0.5)
costheta = 2.*ran2(iseed) - 1.
phi = 2.*pi*(ran2(iseed_mc) - 0.5)
costheta = 2.*ran2(iseed_mc) - 1.
theta = acos(costheta)
sintheta = sin(theta)
sinphi = sin(phi)
Expand Down
186 changes: 181 additions & 5 deletions src/tests/test_gravity.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,6 @@ module testgravity
use io, only:id,master
implicit none
public :: test_gravity

private

contains
Expand All @@ -37,12 +36,14 @@ subroutine test_gravity(ntests,npass,string)
integer, intent(inout) :: ntests,npass
character(len=*), intent(in) :: string
logical :: testdirectsum,test_mom,testtaylorseries,testall,test_plummer
logical :: plot_plummer

testdirectsum = .false.
testtaylorseries = .false.
test_mom = .false.
testall = .false.
test_plummer = .false.
plot_plummer = .false.
select case(string)
case('taylorseries')
testtaylorseries = .true.
Expand All @@ -52,6 +53,8 @@ subroutine test_gravity(ntests,npass,string)
test_mom = .true.
case('spheres','plummer','hernquist')
test_plummer = .true.
case('plotplummer')
plot_plummer = .true.
case default
testall = .true.
end select
Expand All @@ -73,6 +76,10 @@ subroutine test_gravity(ntests,npass,string)
!--unit tests of Plummer and Hernquist spheres
!
if (test_plummer .or. testall) call test_spheres(ntests,npass)
!
!--Plot routine of Plummer and Homogeneous sphere (store data to be plotted)
!
if (plot_plummer) call plot_SFMM()

if (id==master) write(*,"(/,a)") '<-- SELF-GRAVITY TESTS COMPLETE'
else
Expand Down Expand Up @@ -762,7 +769,8 @@ subroutine test_sphere(ntests,npass,iprofile)
use setup_params,only:npart_total
use testutils, only:checkval,update_test_scores
use setplummer, only:get_accel_profile,profile_label,radius_from_mass,density_profile
use spherical, only:set_sphere
use spherical, only:set_sphere,iseed_mc
use kdtree, only:tree_accuracy
use kernel, only:hfact_default
use table_utils, only:linspace
use mpidomain, only:i_belong
Expand Down Expand Up @@ -793,6 +801,7 @@ subroutine test_sphere(ntests,npass,iprofile)

call init_part()
hfact = hfact_default
tree_accuracy = 0.5
gamma = 5./3.
polyk = 0.
ieos = 11
Expand All @@ -816,6 +825,7 @@ subroutine test_sphere(ntests,npass,iprofile)
ref_sum = 0.

do ireal=1,nrealisations
iseed_mc = ireal
npart = 0
npart_total = 0
call set_sphere('random',id,master,rmin,rmax,psep,hfact,npart, &
Expand All @@ -841,25 +851,191 @@ subroutine test_sphere(ntests,npass,iprofile)
err_local = reduceall_mpi('+',err_local)
ref_local = reduceall_mpi('+',ref_local)
if (iverbose > 0 .and. id==master) then
print*,' realisation ',ireal,' mase_local = ',sqrt(err_local/ref_local)
print*,' realisation ',ireal,' mase_local = ',err_local/npart
endif
err_sum = err_sum + err_local
ref_sum = ref_sum + ref_local
enddo

if (ref_sum > tiny(0.)) then
mase = sqrt(err_sum/ref_sum)
mase = err_sum/(nrealisations*npart)
else
mase = 0.
endif

mase_tol = 1.2e-1
mase_tol = 8.5e-4
nfailed = 0
call checkval(mase,0.,mase_tol,nfailed(1),'MASE '//trim(label))
call update_test_scores(ntests,nfailed,npass)

end subroutine test_sphere


!-----------------------------------------------------------------------
!+
! Generate plot data showing the perf and accuracy of self-gravity solver
! in Phantom (see Bernard et al. 2026)
!+
!-----------------------------------------------------------------------
subroutine plot_SFMM()
use setplummer, only:iprofile_plummer
integer :: ntarg(7),i

ntarg = (/1000,3000,10000,30000,100000,300000,1000000/)

if (id==master) write(*,*) '--> Plot routine : Plummer sphere tests with different Npart'
do i=1,size(ntarg)
if (id==master) write(*,*) 'Test with Npart = ',ntarg(i)
call get_plummer_prec_perf(ntarg(i),iprofile_plummer)
enddo

end subroutine plot_SFMM

!-----------------------------------------------------------------------
!+
! test function to compare FMM,SFMM to directsum for Plummer and Homo
! sphere. It tests precision and perf for different theta max and Npart
!+
!-----------------------------------------------------------------------
subroutine get_plummer_prec_perf(npart_target,iprofile)
use dim, only:maxp
use deriv, only:get_derivs_global
use eos, only:gamma,polyk
use mpiutils, only:reduceall_mpi
use options, only:ieos,alpha,alphau,alphaB,tolh
use part, only:init_part,npart,xyzh,fxyzu,hfact,&
npartoftype,massoftype,istar,maxphase,iphase,isetphase,rhoh
use setup_params,only:npart_total
use testutils, only:checkval,update_test_scores
use setplummer, only:get_accel_profile,profile_label,radius_from_mass,density_profile
use spherical, only:set_sphere,iseed_mc
use random, only:ran2
use kdtree, only:tree_accuracy
use kernel, only:hfact_default
use table_utils, only:linspace
use mpidomain, only:i_belong
use io, only:id,master,iverbose
use neighkdtree, only:neigh_switch
use timing, only:get_timings
use sortutils, only:indexx
integer, intent(in) :: iprofile,npart_target
integer :: i,it,itest
integer, parameter :: niter=10
real, allocatable :: fxyz_dir(:,:),err_rel(:)
integer,allocatable :: erridx(:)
real :: rsoft,mass_total,cut_fraction,rmin,rmax,psep,theta_crit
character(len=64) :: label,filename_max,type
integer, parameter :: ntab = 1000
real :: rgrid(ntab),rhotab(ntab)
real :: maxerr(3,niter+1),minerr(3,niter+1),meanerr(3,niter+1)
integer :: iunit
real(kind=4) :: t1,t2,tcpu1,tcpu2,timings(3,niter+1)

label = profile_label(iprofile)

write(filename_max,'("data_plot_plummer_sphere_",i8.8,".ev")') npart_target
filename_max = adjustl(filename_max)

open(newunit=iunit,file=trim(filename_max),action='write',status='replace')
write(iunit,"(a)") '# \theta, emax_SFMM, emax_FMM, emin_SFMM, emin_FMM, &
&tcpu_SFMM, tcpu_FMM, tcpu_direct'



call init_part()
hfact = hfact_default
gamma = 5./3.
polyk = 0.
ieos = 11
alpha = 0.; alphau = 0.; alphaB = 0.
tolh = 1.e-5
rsoft = 1.0
mass_total = 1.0

! construct tables for radius and density
cut_fraction = 0.999
rmin = 0.
rmax = radius_from_mass(iprofile,cut_fraction,rsoft)
call linspace(rgrid,0.,rmax)

do i=1,ntab
rhotab(i) = density_profile(iprofile,rgrid(i),rsoft,mass_total)
enddo

psep = rmax/real(ntab) ! this is not used for random placement anyway
iverbose = 1

iseed_mc = 1
npart = 0
npart_total = 0

call set_sphere('random',id,master,rmin,rmax,psep,hfact,npart, &
xyzh,npart_total,rhotab=rhotab,rtab=rgrid,exactN=.true.,&
np_requested=npart_target,mask=i_belong,verbose=.false.)
!call set_sphere('random',id,master,rmin,rmax,psep,hfact,npart,xyzh,npart_total,np_requested=npart_target)

massoftype(istar) = mass_total/real(npart_total)
npartoftype(istar) = npart

if (maxphase==maxp) then
iphase(1:npart) = isetphase(istar,iactive=.true.)
endif

allocate(fxyz_dir(3,npart))
allocate(err_rel(npart))
allocate(erridx(npart))


call get_derivs_global(icall=1)


tree_acc: do it=0,niter
theta_crit = 0.1 + it*0.05
do itest=3,1,-1
if (itest==1) then
type = "SFMM"
neigh_switch = 1
tree_accuracy = theta_crit
elseif (itest==2) then
type = "FMM"
neigh_switch = 2
tree_accuracy = theta_crit
else
type = "direct"
neigh_switch = 2
tree_accuracy = 0.
endif

if (itest==3 .and. it>0) cycle
call get_timings(t1,tcpu1)
call get_derivs_global(icall=2)
call get_timings(t2,tcpu2)

if (itest==3) fxyz_dir = fxyzu(0:3,1:npart)

err_rel = norm2(fxyzu(1:3,1:npart)-fxyz_dir,1)/norm2(fxyz_dir,1)
call indexx(npart, err_rel, erridx)

maxerr(itest,it+1) = maxval(err_rel)
minerr(itest,it+1) = err_rel(erridx(npart/10))
meanerr(itest,it+1) = sum(err_rel)/npart

if (itest==3) then
timings(itest,1:niter+1) = tcpu2-tcpu1
else
timings(itest,it+1) = tcpu2-tcpu1
endif
enddo
enddo tree_acc

do it=0,niter
write(iunit,*) 0.1 + it*0.05, maxerr(1:2,it+1), minerr(1:2,it+1), meanerr(1:2,it+1), timings(1:3,it+1)
enddo

close(iunit)

end subroutine get_plummer_prec_perf

!-----------------------------------------------------------------------
!+
! Copy gas particles to sinks
Expand Down
3 changes: 1 addition & 2 deletions src/tests/test_neigh.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ subroutine test_neigh(ntests,npass)
use mpidomain, only:i_belong
use part, only:maxphase,iphase,isetphase,igas,iactive,isdead_or_accreted
use testutils, only:checkval,checkvalbuf_start,checkvalbuf,checkvalbuf_end,update_test_scores
use neighkdtree, only:build_tree,get_neighbour_list,ncells,leaf_is_active,force_dual_walk
use neighkdtree, only:build_tree,get_neighbour_list,ncells,leaf_is_active
use kdtree, only:inodeparts,inoderange,tree_accuracy
use mpidomain, only:i_belong
use boundary, only:xmin,xmax,ymin,ymax,zmin,zmax,dybound,dzbound
Expand Down Expand Up @@ -103,7 +103,6 @@ subroutine test_neigh(ntests,npass)
!print*,'thread ',id,' npart = ',npart
!iverbose = 3

force_dual_walk = .false.
tree_accuracy = 0.5
rhozero = 7.5
hfact = 1.2
Expand Down
2 changes: 1 addition & 1 deletion src/tests/testsuite.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -194,7 +194,7 @@ subroutine testsuite(string,first,last,ntests,npass,nfail)
dosedov = .true.
case('indtstep','ind')
doindtstep = .true.
case('gravity','grav','plummer','hernquist','fmm','taylorseries','directsum')
case('gravity','grav','plummer','hernquist','fmm','taylorseries','directsum','plotplummer')
dogravity = .true.
case('dump','rwdump','dumprw')
dorwdump = .true.
Expand Down
Loading