!
! This module implements gas/dust self-gravity using a Barnes-Hut octree
! algorithm. The main function, do_barneshut, takes the variable 'phi' as input,
! which is the gas/dust density grid from the selfgravity module, and modifies
! 'phi' in place, becoming the gravitational potential.
!
! To use this module:
! 1. In Makefile.local, set POISSON=experimental/barneshut (FOURIER is
! subsequently not needed)
! 2. Update start.in with &poisson_init_pars/appropriate parameters.
!
! The variable 'themap' contains a complete list of the region/point pairings
! for the BH map for each processor (where 'regions' are being integrated to
! calculate the potential at each 'point'). Each such pairing is given by ten
! integers: the indices of the local point coordinates (i,j,k), the upper/lower
! boundaries of the region in each dimension, and the processor ID of the
! region. The subroutine 'mkmap' is run twice: first, a "dummy" run of the
! subroutine calculates the total number of regions, after which 'themap' is
! allocated. A second pass then populates 'themap' with the needed information.
! The memory required per core for the map for a reasonably-sized run in
! spherical coordinates is of order 10^2 MB. The runtime then consists of
! MPI-based exchange of density fields between processors, then a single loop
! over the point-region pairings in 'themap' takes care of the integration.
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: lpoisson=.true.
!
! MVAR CONTRIBUTION 0
! MAUX CONTRIBUTION 0
!
!***************************************************************
module Poisson
!
use Cdata
use Cparam
use Messages
!
implicit none
!
real :: octree_theta = 0.5
logical :: lshowtime = .false. ! Print BH integration at each time step for
! root processor. Debugging/testing only.
logical :: lsquareregions = .false. ! .true. = make regions as square-shaped as
! possible. Results in fewer regions.
logical :: lprecalcdists = .false. ! Pre-calculate point-region distances.
! Not technically correct, but may work as an
! approximation for sufficiently small octree_theta,
! and provides a significant speedup.
logical :: lnorepeatsumming = .true. ! (very) slow to start, quick to run
logical :: lmakecartoon = .false. ! print some map data to stdout for tree illustration purposes
logical :: ltreestatus = .false.
logical :: lreadoctree = .false.
logical :: lwriteoctree = .false.
real :: octree_maxdist = 1e308
real :: octree_smoothdist = 0.15
real, dimension (nx,ny,nz) :: vols
integer, dimension (0:ncpus-1) :: sx, sy, sz
real, dimension (nx) :: xc
real, dimension (ny) :: yc
real, dimension (nz) :: zc
real, dimension (nx,ny,nz) :: xmesh, ymesh, zmesh
real, dimension (nx,0:ncpus-1) :: xrecv
real, dimension (ny,0:ncpus-1) :: yrecv
real, dimension (nz,0:ncpus-1) :: zrecv
integer, allocatable :: themap_group(:,:), themap_single(:,:)
real, allocatable :: regdist1_single(:), regsmooth_single(:)
real, allocatable :: regdist1_group(:), regsmooth_group(:)
logical, allocatable :: luseprevioussum(:)
integer :: nsingle, ngroup, iregion, irl, iru, jrl, jru, krl, kru, nprecalc=0
integer, parameter :: nlt = 1e7 ! 3*80MB for 1e7
real :: lkmin, lkmax, dxlt
real, dimension(nlt) :: xlt, sinlt, coslt
!
namelist /poisson_init_pars/ &
octree_theta, lshowtime, lsquareregions, lprecalcdists, octree_maxdist, &
octree_smoothdist, lnorepeatsumming, ltreestatus, lwriteoctree, lreadoctree, &
lmakecartoon
!
namelist /poisson_run_pars/ &
octree_theta, lshowtime, lsquareregions, lprecalcdists, octree_maxdist, &
octree_smoothdist, lnorepeatsumming, ltreestatus, lwriteoctree, lreadoctree
!
contains
!***********************************************************************
subroutine inverse_laplacian(phi)
!
use General, only: keep_compiler_quiet
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (nx,ny,nz) :: phi
!
integer :: pp, i, j, k, xs, ys, zs, ii, jj, kk
intent(inout) :: phi
!
if (modulo(log(real(nx))/log(2.0),1.0) .gt. tini .or. &
modulo(log(real(ny))/log(2.0),1.0) .gt. tini .or. &
modulo(log(real(nz))/log(2.0),1.0) .gt. tini) then
if (lroot) print*,'Grid zones per processor in each axis '//&
'must be a power of 2.'
call fatal_error("barneshut","")
endif
!
call do_barneshut(phi)
call keep_compiler_quiet(f)
!
endsubroutine inverse_laplacian
!***********************************************************************
subroutine initialize_poisson
!
use Mpicomm
!
real :: xw, yw, zw, dx1, dy1, dz1
integer, dimension (nx) :: xind
integer, dimension (ny) :: yind
integer, dimension (nz) :: zind
real, dimension (nx) :: dxc_1
real, dimension (ny) :: dyc_1
real, dimension (nz) :: dzc_1
integer :: pp, i, j, k, xs, ys, zs, ii, jj, kk
integer :: iprecalc, iregtmp, info, itmp
logical :: lprecalc = .false., cartooncond, doingsort
logical, allocatable :: tmpmask(:)
character :: treestatus
!
! Trimming ghost zones
xc = x(l1:l2)
yc = y(m1:m2)
zc = z(n1:n2)
dxc_1 = dx_1(l1:l2)
dyc_1 = dy_1(m1:m2)
dzc_1 = dz_1(n1:n2)
!
do i=1,nx
do j=1,ny
do k=1,nz
if (lspherical_coords) then
xmesh(i,j,k) = xc(i)*sin(yc(j))*cos(zc(k))
ymesh(i,j,k) = xc(i)*sin(yc(j))*sin(zc(k))
zmesh(i,j,k) = xc(i)*cos(yc(j))
elseif (lcylindrical_coords) then
xmesh(i,j,k) = xc(i)*cos(yc(j))
ymesh(i,j,k) = xc(i)*sin(yc(j))
zmesh(i,j,k) = zc(k)
elseif (lcartesian_coords) then
xmesh(i,j,k) = xc(i)
ymesh(i,j,k) = yc(j)
zmesh(i,j,k) = zc(k)
endif
enddo
enddo
enddo
!
do pp=0,ncpus-1
if (pp/=iproc) then
call mpisendrecv_real(xc,nx,pp,281,xrecv(:,pp),pp,281)
call mpisendrecv_real(yc,ny,pp,282,yrecv(:,pp),pp,282)
call mpisendrecv_real(zc,nz,pp,283,zrecv(:,pp),pp,283)
endif
enddo
!
xrecv(:,iproc) = xc
yrecv(:,iproc) = yc
zrecv(:,iproc) = zc
!
! Calculate volumes for converting density to mass
do i=1,nx
do j=1,ny
do k=1,nz
if (nxgrid ==1) then
dx1=1.0!/Lxyz(1)
else
dx1=dxc_1(i)
endif
if (nygrid ==1) then
dy1=1.0!/Lxyz(2)
else
dy1=dyc_1(j)
endif
if (nzgrid ==1) then
dz1=1.0!/Lxyz(3)
else
dz1=dzc_1(k)
endif
vols(i,j,k) = 1.0/(dx1*dy1*dz1)
if (lcylindrical_coords .and. nygrid/=1) then
vols(i,j,k) = vols(i,j,k)*xc(i)
elseif (lspherical_coords) then
if (nygrid/=1) vols(i,j,k) = vols(i,j,k)*xc(i)
if (nzgrid/=1) vols(i,j,k) = vols(i,j,k)*sin(yc(j))*xc(i)
endif
enddo
enddo
enddo
!
if (lsquareregions) then
do pp=0,ncpus-1
if (lspherical_coords) then
xw = abs(xrecv(nx,pp)-xrecv(1,pp))
yw = abs(yrecv(ny,pp)-yrecv(1,pp))*xrecv(nx/2,pp)
zw = abs(zrecv(nz,pp)-zrecv(1,pp))*xrecv(nx/2,pp)*sin(yrecv(ny/2,pp)) ! this needs to get fixed for 2D case (where ny/2=0)
elseif (lcylindrical_coords) then
xw = abs(xrecv(nx,pp)-xrecv(1,pp))
yw = abs(yrecv(ny,pp)-yrecv(1,pp))*xrecv(nx/2,pp)
zw = abs(zrecv(nz,pp)-zrecv(1,pp))
elseif (lcartesian_coords) then
xw = abs(xrecv(nx,pp)-xrecv(1,pp))
yw = abs(yrecv(ny,pp)-yrecv(1,pp))
zw = abs(zrecv(nz,pp)-zrecv(1,pp))
endif
! If a region is longer along one axis (or two), it will be split along
! that axis (or axes) to make the resulting sub-regions roughly cubical.
sx(pp) = max(nx/max(roundtwo(xw/yw),roundtwo(xw/zw),1),1)
sy(pp) = max(ny/max(roundtwo(yw/xw),roundtwo(yw/zw),1),1)
sz(pp) = max(nz/max(roundtwo(zw/xw),roundtwo(zw/yw),1),1)
enddo
endif
!
! Lookup table setup:
if (lspherical_coords .or. lcylindrical_coords) then
lkmin = -pi-dy-dz
lkmax = pi+dy+dz
dxlt = (lkmax-lkmin)/real(nlt)
do i=1,nlt
xlt(i) = lkmin+dxlt*real(i-1)
sinlt(i) = sin(xlt(i))
coslt(i) = cos(xlt(i))
enddo
endif
!
do i=1,nx
xind(i) = i
enddo
do j=1,ny
yind(j) = j
enddo
do k=1,nz
zind(k) = k
enddo
!
! First pass only counts regions, second pass actually populates 'themap'
do iprecalc=1,2
if (lprecalc) then
! if (lroot) print*,"barneshut: 1x1x1 ~regions/proc:",nsingle,"; >1x1x1 ~regions/proc:",ngroup
if (lroot) print '("barneshut: ",I0," 1x1x1 and ",I0," >1x1x1 region pairings on proc 0")', nsingle, ngroup
allocate(themap_group(10,ngroup))
allocate(themap_single(10,nsingle))
allocate(regsmooth_single(nsingle))
allocate(regsmooth_group(ngroup))
regsmooth_single = 1.0 ; regsmooth_group = 1.0
allocate(regdist1_single(nsingle))
if (lprecalcdists) allocate(regdist1_group(ngroup))
if (lreadoctree) then
call read_octree
exit
endif
endif
nsingle = 0 ! Advanced inside 'mkmap'
ngroup = 0 ! Advanced inside 'mkmap'
do pp=0,ncpus-1
if (lsquareregions) then
do kk=1,nz/sz(pp)
do jj=1,ny/sy(pp)
do ii=1,nx/sx(pp)
do k=1,nz
do j=1,ny
do i=1,nx
irl = (ii-1)*sx(pp)+1 ; iru = ii*sx(pp)
jrl = (jj-1)*sy(pp)+1 ; jru = jj*sy(pp)
krl = (kk-1)*sz(pp)+1 ; kru = kk*sz(pp)
call mkmap((/i,j,k/),(/sx(pp),sy(pp),sz(pp)/), &
xrecv(irl:iru,pp),yrecv(jrl:jru,pp),zrecv(krl:kru,pp), &
xind(irl:iru),yind(jrl:jru),zind(krl:kru),pp,lprecalc)
enddo !i
enddo !j
enddo !k
enddo !ii
enddo !jj
enddo !kk
else
do k=1,nz
do j=1,ny
do i=1,nx
call mkmap((/i,j,k/),(/nx,ny,nz/), xrecv(:,pp),yrecv(:,pp), &
zrecv(:,pp),xind,yind,zind,pp,lprecalc)
if (ltreestatus.and.lprecalc) then
write(*,char(nsingle+ngroup))
endif
enddo !i
enddo !j
enddo !k
endif !lsquareregions
enddo !pp
lprecalc = .true.
enddo !iprecalc
!
allocate(luseprevioussum(ngroup))
luseprevioussum = .false.
if (lnorepeatsumming) then
do iregion=2,ngroup
if (all(themap_group(4:,iregion).eq.themap_group(4:,iregion-1))) &
luseprevioussum(iregion) = .true.
enddo
endif
!
if (lwriteoctree) call write_octree
!
! if (lmakecartoon) then
! do iregion=1,nregions
! if (lcylindrical_coords.or.lcartesian_coords) then
! cartooncond = (themap(2,iregion).eq.(ny/2))
! else
! cartooncond = (themap(3,iregion).eq.(nz/2))
! endif
! if ((themap(1,iregion).eq.(nx/2)).and.cartooncond) then
! print*,themap(:,iregion)
! endif
! enddo
! endif
!
endsubroutine initialize_poisson
!***********************************************************************
subroutine do_barneshut(phi) ! 'phi' is density from selfgravity module,
! will be transformed into potential.
!
use Mpicomm
!
real, dimension (nx,ny,nz,0:ncpus-1) :: phirecv=0.0
real, dimension (nx,ny,nz) :: phi
real :: xreg, yreg, zreg, summreg, summreg1
real :: tstart_loop1, tstop_loop1, tstart_loop2, tstop_loop2, tstart_mpi, tstop_mpi
integer :: pp, i, j, k, xs, ys, zs, ii, jj, kk
!
phi = phi*vols ! 'phi' now in mass units
!
! Send masses to other processors
if (lroot .and. lshowtime) call cpu_time(tstart_mpi)
do pp=0,ncpus-1
if (pp/=iproc) then
call mpisendrecv_real(phi,(/nx,ny,nz/),pp,284,phirecv(:,:,:,pp),pp,284)
endif
enddo
!
if (lroot .and. lshowtime) call cpu_time(tstop_mpi)
!
phirecv(:,:,:,iproc) = phi
!
phi = 0.0
!
if (lroot .and. lshowtime) call cpu_time(tstart_loop1)
do iregion=1,nsingle
i = themap_single(1, iregion) ; irl = themap_single(5,iregion) ; iru = themap_single(6,iregion)
j = themap_single(2, iregion) ; jrl = themap_single(7,iregion) ; jru = themap_single(8,iregion)
k = themap_single(3, iregion) ; krl = themap_single(9,iregion) ; kru = themap_single(10,iregion)
pp = themap_single(4, iregion)
phi(i,j,k) = phi(i,j,k) - regsmooth_single(iregion)* &
sum(phirecv(irl:iru,jrl:jru,krl:kru,pp))*regdist1_single(iregion)
enddo
if (lroot .and. lshowtime) call cpu_time(tstop_loop1)
if (lroot .and. lshowtime) call cpu_time(tstart_loop2)
if (.not.lprecalcdists) then
do iregion=1,ngroup
i = themap_group(1,iregion)
j = themap_group(2,iregion)
k = themap_group(3,iregion)
if (.not.luseprevioussum(iregion)) then
pp = themap_group(4,iregion)
irl = themap_group(5,iregion) ; iru = themap_group(6,iregion)
jrl = themap_group(7,iregion) ; jru = themap_group(8,iregion)
krl = themap_group(9,iregion) ; kru = themap_group(10,iregion)
summreg = sum(phirecv(irl:iru,jrl:jru,krl:kru,pp))
summreg1 = 1.0/summreg
xreg = sum(xrecv(irl:iru,pp) &
*sum(sum(phirecv(irl:iru,jrl:jru,krl:kru,pp),3),2))*summreg1
yreg = sum(yrecv(jrl:jru,pp) &
*sum(sum(phirecv(irl:iru,jrl:jru,krl:kru,pp),3),1))*summreg1
zreg = sum(zrecv(krl:kru,pp) &
*sum(sum(phirecv(irl:iru,jrl:jru,krl:kru,pp),2),1))*summreg1
endif
phi(i,j,k) = phi(i,j,k) - summreg*regsmooth_group(iregion)/get_dist((/i,j,k/),(/xreg,yreg,zreg/))
enddo
else
do iregion=1,ngroup
i = themap_group(1,iregion)
j = themap_group(2,iregion)
k = themap_group(3,iregion)
pp = themap_group(4,iregion)
irl = themap_group(5,iregion) ; iru = themap_group(6,iregion)
jrl = themap_group(7,iregion) ; jru = themap_group(8,iregion)
krl = themap_group(9,iregion) ; kru = themap_group(10,iregion)
if (.not.luseprevioussum(iregion)) then
summreg = sum(phirecv(irl:iru,jrl:jru,krl:kru,pp))
endif
phi(i,j,k) = phi(i,j,k) - summreg*regsmooth_group(iregion)*regdist1_group(iregion)
enddo
endif
!
if (lroot .and. lshowtime) call cpu_time(tstop_loop2)
!
phi = phi/(4.0*pi) ! The selfgravity module is going to multiply by a factor
! of 4pi that we don't want (I think).
!
if (lroot .and. lshowtime) then
print '("barneshut: MPI time = ",f10.6," seconds.")',tstop_mpi-tstart_mpi
print '("barneshut: Loop1 time = ",f10.6," seconds.")',tstop_loop1-tstart_loop1
print '("barneshut: Loop2 time = ",f10.6," seconds.")',tstop_loop2-tstop_loop1
endif
!
endsubroutine do_barneshut
!***********************************************************************
recursive subroutine mkmap(ipos,dimi,xi,yi,zi,xsind,ysind,zsind,ppi,laddreg)
integer, dimension(3) :: ipos, dimi ! Local position; dimensions of region
real, dimension(dimi(1)) :: xi ! Region coordinates
real, dimension(dimi(2)) :: yi !
real, dimension(dimi(3)) :: zi !
integer, dimension(dimi(1)) :: xsind ! Region indices (for building map)
integer, dimension(dimi(2)) :: ysind !
integer, dimension(dimi(3)) :: zsind !
integer, dimension(10) :: newregion
integer :: sx, sy, sz, ii, ji, ki, il, iu, jl, ju, kl, ku, ppi, newindex
real :: lxi, lyi, lzi, dist, xwi, ywi, zwi, width
logical :: laddreg, lsplit, lregionmatched, lnotself
!
! If the point where the potential is being calculated is inside the
! region in question, then we MUST subdivide further (unless the region is
! a single cell-- this is taken care of in the 'if' statement below).
lsplit = ((ipos(1) .ge. xsind(1) .and. ipos(1) .le. xsind(dimi(1))) &
.and. (ipos(2) .ge. ysind(1) .and. ipos(2) .le. ysind(dimi(2))) &
.and. (ipos(3) .ge. zsind(1) .and. ipos(3) .le. zsind(dimi(3))) &
.and. (iproc .eq. ppi))
lnotself = (((ipos(1) .ne. xsind(1) .or. ipos(1) .ne. xsind(dimi(1))) &
.or. (ipos(2) .ne. ysind(1) .or. ipos(2) .ne. ysind(dimi(2))) &
.or. (ipos(3) .ne. zsind(1) .or. ipos(3) .ne. zsind(dimi(3))) &
.or. (iproc .ne. ppi)))
!
lxi = sum(xi)/real(dimi(1))
lyi = sum(yi)/real(dimi(2))
lzi = sum(zi)/real(dimi(3))
dist = get_dist(ipos,(/lxi,lyi,lzi/))
xwi = abs(xi(dimi(1))-xi(1))
ywi = abs(yi(dimi(2))-yi(1))
zwi = abs(zi(dimi(3))-zi(1))
if (lspherical_coords) then
width = max(xwi,lxi*ywi,lxi*zwi*sin(lyi))
elseif (lcylindrical_coords) then
width = max(xwi,lxi*ywi,zwi)
elseif (lcartesian_coords) then
width = max(xwi,ywi,zwi)
endif
if ((width/dist > octree_theta .or. lsplit) .and. product(dimi)/=1) then
sx = max(dimi(1)/2,1)
sy = max(dimi(2)/2,1)
sz = max(dimi(3)/2,1)
do ki=1,dimi(3)/sz
do ji=1,dimi(2)/sy
do ii=1,dimi(1)/sx
il = sx*(ii-1)+1; iu = sx*ii
jl = sy*(ji-1)+1; ju = sy*ji
kl = sz*(ki-1)+1; ku = sz*ki
call mkmap(ipos,(/sx,sy,sz/),xi(il:iu),yi(jl:ju),zi(kl:ku), &
xsind(il:iu),ysind(jl:ju),zsind(kl:ku),ppi,laddreg)
enddo
enddo
enddo
elseif (((.not. lsplit).or.(product(dimi).eq.1)) .and. (dist .lt. octree_maxdist) .and. lnotself) then
newregion = (/ipos(1),ipos(2),ipos(3), ppi, xsind(1), &
xsind(dimi(1)), ysind(1),ysind(dimi(2)),zsind(1),zsind(dimi(3))/)
if ((newregion(5).ne.newregion(6)).or.(newregion(7).ne.newregion(8)).or. &
(newregion(9).ne.newregion(10))) then
if (laddreg) then
lregionmatched = .false.
if (lnorepeatsumming) then
do iregion=1,ngroup
if (all(newregion(4:10).eq.themap_group(4:10,iregion))) then
themap_group(:,iregion+2:ngroup+1) = themap_group(:,iregion+1:ngroup)
themap_group(:,iregion+1) = newregion
lregionmatched = .true.
exit
endif
enddo
endif
if (.not.lregionmatched) themap_group(:,ngroup+1) = newregion
endif
ngroup = ngroup+1
if (lprecalcdists) regdist1_group(ngroup) = 1.0/dist
if (dist .gt. (octree_maxdist-octree_smoothdist)) then
regsmooth_group(ngroup) = 0.5*(cos(pi*((octree_maxdist-dist)/octree_smoothdist+1.0))+1.0)
endif
else
nsingle = nsingle+1
if (laddreg) then
themap_single(1:10,nsingle) = newregion(1:10)
regdist1_single(nsingle) = 1.0/dist
endif
if (dist .gt. (octree_maxdist-octree_smoothdist)) then
regsmooth_single(nsingle) = 0.5*(cos(pi*((octree_maxdist-dist)/octree_smoothdist+1.0))+1.0)
endif
endif
endif
!
endsubroutine mkmap
!***********************************************************************
function get_dist(p1,p2) result(rdist)
integer, dimension(3) :: p1
real, dimension(3) :: p2
real :: x2, y2, z2, rdist
real, dimension(2) :: sincosy, sincosz
!
if (lspherical_coords) then
sincosy = sincoslf(p2(2)) ! returns (/sin, cos/)
sincosz = sincoslf(p2(3))
x2 = p2(1)*sincosy(1)*sincosz(2)
y2 = p2(1)*sincosy(1)*sincosz(1)
z2 = p2(1)*sincosy(2)
elseif (lcylindrical_coords) then
sincosy = sincoslf(p2(2))
x2 = p2(1)*sincosy(2)
y2 = p2(1)*sincosy(1)
z2 = p2(3)
elseif (lcartesian_coords) then
x2 = p2(1)
y2 = p2(2)
z2 = p2(3)
endif
rdist = sqrt((xmesh(p1(1),p1(2),p1(3))-x2)**2+ &
(ymesh(p1(1),p1(2),p1(3))-y2)**2+(zmesh(p1(1),p1(2),p1(3))-z2)**2)
!
endfunction get_dist
!***********************************************************************
function sincoslf(ang) ! returns (/sin,cos/)
real,dimension(2) :: sincoslf
real :: ang, sinx, cosx !, a0, a1, a2
integer :: ilk
!
ilk = nint((ang-lkmin)/dxlt)+1
sinx = sinlt(ilk)
cosx = coslt(ilk)
sincoslf = (/sinx,cosx/)
!
endfunction sincoslf
!***********************************************************************
function roundtwo(rin) result(iout)
real :: rin
integer :: iout
!
iout = 2**nint(log(rin)/log(2.0))
endfunction roundtwo
!***********************************************************************
subroutine write_octree
call system('mkdir -p '//trim(directory_snap)//'/bhmap')
open(unit = 10, status='replace', &
file=trim(directory_snap)//trim('/bhmap/themap_single.dat'),form='unformatted')
write(10) themap_single
close(10)
open(unit = 10, status='replace', &
file=trim(directory_snap)//trim('/bhmap/themap_group.dat'),form='unformatted')
write(10) themap_group
close(10)
open(unit = 10, status='replace', &
file=trim(directory_snap)//trim('/bhmap/regsmooth_single.dat'),form='unformatted')
write(10) regsmooth_single
close(10)
open(unit = 10, status='replace', &
file=trim(directory_snap)//trim('/bhmap/regsmooth_group.dat'),form='unformatted')
write(10) regsmooth_group
close(10)
open(unit = 10, status='replace', &
file=trim(directory_snap)//trim('/bhmap/regdist1_single.dat'),form='unformatted')
write(10) regdist1_single
close(10)
if (lprecalcdists) then
open(unit = 10, status='replace', &
file=trim(directory_snap)//trim('/bhmap/regdist1_group.dat'),form='unformatted')
write(10) regdist1_group
close(10)
endif
endsubroutine write_octree
!***********************************************************************
subroutine read_octree
open(unit = 10, status='old', &
file=trim(directory_snap)//trim('/bhmap/themap_single.dat'),form='unformatted')
read(10) themap_single
close(10)
open(unit = 10, status='old', &
file=trim(directory_snap)//trim('/bhmap/themap_group.dat'),form='unformatted')
read(10) themap_group
close(10)
open(unit = 10, status='old', &
file=trim(directory_snap)//trim('/bhmap/regsmooth_single.dat'),form='unformatted')
read(10) regsmooth_single
close(10)
open(unit = 10, status='old', &
file=trim(directory_snap)//trim('/bhmap/regsmooth_group.dat'),form='unformatted')
read(10) regsmooth_single
close(10)
open(unit = 10, status='old', &
file=trim(directory_snap)//trim('/bhmap/regdist1_single.dat'),form='unformatted')
read(10) regdist1_single
close(10)
if (lprecalcdists) then
open(unit = 10, status='old', &
file=trim(directory_snap)//trim('/bhmap/regdist1_group.dat'),form='unformatted')
read(10) regdist1_group
close(10)
endif
endsubroutine read_octree
!***********************************************************************
subroutine read_poisson_init_pars(iostat)
!
use File_io, only: parallel_unit
!
integer, intent(out) :: iostat
!
read(parallel_unit, NML=poisson_init_pars, IOSTAT=iostat)
!
endsubroutine read_poisson_init_pars
!***********************************************************************
subroutine write_poisson_init_pars(unit)
!
integer, intent(in) :: unit
!
write(unit, NML=poisson_init_pars)
!
endsubroutine write_poisson_init_pars
!***********************************************************************
subroutine read_poisson_run_pars(iostat)
!
use File_io, only: parallel_unit
!
integer, intent(out) :: iostat
!
read(parallel_unit, NML=poisson_run_pars, IOSTAT=iostat)
!
endsubroutine read_poisson_run_pars
!***********************************************************************
subroutine write_poisson_run_pars(unit)
!
integer, intent(in) :: unit
!
write(unit, NML=poisson_run_pars)
!
endsubroutine write_poisson_run_pars
!***********************************************************************
subroutine inverse_laplacian_semispectral(f,phi)
!
! Solve the Poisson equation by Fourier transforming on a periodic grid.
!
! 15-may-2006/anders+jeff: dummy
!
use General, only: keep_compiler_quiet
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (nx,ny,nz) :: phi
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(phi)
call keep_compiler_quiet(f)
!
endsubroutine inverse_laplacian_semispectral
!***********************************************************************
subroutine get_acceleration(acceleration)
!
use General, only: keep_compiler_quiet
!
real, dimension(nx,ny,nz,3), intent(out) :: acceleration !should I (CAN I?) make this allocatable?
!
call keep_compiler_quiet(acceleration)
!
endsubroutine get_acceleration
!***********************************************************************
endmodule Poisson