-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathdirect_sum_poisson3D.f90
97 lines (89 loc) · 3.7 KB
/
direct_sum_poisson3D.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
!------------------------------------------------------------------------------!
! NDSPMHD: A Smoothed Particle (Magneto)Hydrodynamics code for (astrophysical) !
! fluid dynamics simulations in 1, 2 and 3 spatial dimensions. !
! !
! (c) 2002-2015 Daniel Price !
! !
! http://users.monash.edu.au/~dprice/ndspmhd !
! daniel.price@monash.edu -or- dprice@cantab.net (forwards to current address) !
! !
! NDSPMHD comes with ABSOLUTELY NO WARRANTY. !
! This is free software; and you are welcome to redistribute !
! it under the terms of the GNU General Public License !
! (see LICENSE file for details) and the provision that !
! this notice remains intact. If you modify this file, please !
! note section 2a) of the GPLv2 states that: !
! !
! a) You must cause the modified files to carry prominent notices !
! stating that you changed the files and the date of any change. !
! !
! ChangeLog: !
!------------------------------------------------------------------------------!
!!----------------------------------------------------------------------------
!! Calculates the solution to any Poisson equation
!!
!! \nabla^2 \phi = \eta
!!
!! by a direct summation over the particles.
!!
!! Use this to check the accuracy of the tree code
!!
!! Input:
!!
!! x(ndim,ntot) : co-ordinates of the particles
!! source(ntot) : quantity to be summed over
!! for an arbitrary quantity \eta, this is given by
!! source = particle mass * eta / (4.*pi*density)
!! ie. source = particle mass for gravity
!! psoft : softening length for plummer softening
!!
!! Output:
!!
!! phitot : total potential phi
!! gradphi(ndim,ntot) : gradient of the potential (for gravity this = force)
!!----------------------------------------------------------------------------
subroutine direct_sum_poisson(x,source,phitot,gradphi,psoft,ntot)
use dimen_mhd, only:ndim
use debug, only:trace
use loguns, only:iprint
implicit none
integer, intent(in) :: ntot
real, dimension(ndim,ntot), intent(in) :: x
real, dimension(ntot), intent(in) :: source
real, dimension(ntot) :: phi
real, dimension(ndim,ntot), intent(out) :: gradphi
real, intent(in) :: psoft
real, intent(out) :: phitot
integer :: i,j
real, dimension(ndim) :: dx,term
real :: rij,rij2,sourcei
!
!--allow for tracing flow
!
if (trace) write(iprint,*) ' Entering subroutine direct_sum_poisson'
!
!--reset potential (but not force) initially
!
phi = 0.
!
!--calculate gravitational force by direct summation
!
do i=1,ntot
sourcei = source(i)
do j=i+1,ntot
dx = x(:,i) - x(:,j)
rij2 = dot_product(dx,dx) + psoft**2
rij = sqrt(rij2)
term(:) = dx(:)/(rij*rij2)
phi(i) = phi(i) - source(j)/rij
phi(j) = phi(j) - sourcei/rij
gradphi(:,i) = gradphi(:,i) - source(j)*term(:)
gradphi(:,j) = gradphi(:,j) + sourcei*term(:)
enddo
enddo
phitot = 0.
do i=1,ntot
phitot = phitot + 0.5*source(i)*phi(i)
enddo
return
end subroutine direct_sum_poisson