!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! UPDATE_CLOUD_STRUCTURE Feb. 2007 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine update_cloud_structure (cl,os,params,ninject,nremove,istep)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this file contains the operations done on a 3D cloud of particles
! to check the density of particles, accordingly add or remove particles
! and transfer the fields carried by the cloud (here the original
! position of the particles)
! First it initializes the fields to the current (thus original) position
! then we also interpolate the fields (original position) onto the
! corresponding fields of the octree solve in order for the information
! to be available during matrix building operation
! cl is the cloud
! os is the octree solve
! npmin is the minimum number of particles in any leaf of the octree solve
! npmax is the maximum number of particles in the smallest leaves of the
! octree solve
! ninject is number of new particles injected in the cloud
! nremove is number of particles removed from the cloud
! levelmax_oct is maximum level of octree solve
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use definitions
!use mpi
implicit none
include 'mpif.h'
type (cloud) cl
type (octreesolve) os
type (parameters) params
integer ninject
integer nremove
type (octreev) ov
integer istep
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer ierr,iproc,nproc,i,j,k,ic,jc,kk,ii
integer np,ne,nn,npnew,ie,npswap,ntest,ninjectp,kinject,kkinject,nchoice
integer,dimension(:),allocatable::loc,nloc,ntot,remove,lev,done
integer err,leaf,level,locc,ndist,nisolated,ntodo,nremovep
integer, dimension(:,:), allocatable :: imemswap
double precision x0,y0,z0,dxyz,distmax,dist,distmin,xx,yy,zz,h(8),ww,r,s,t
double precision,dimension(:),allocatable::xtest,ytest,ztest
double precision,dimension(:),allocatable::xinject,yinject,zinject
double precision,dimension(:),allocatable::straininject,x0inject,y0inject,z0inject
double precision,dimension(:),allocatable::e2dpinject
double precision,dimension(:),allocatable::strain,w,e2dp
double precision,dimension(:,:),allocatable::memswap
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
np=cl%np
ne=os%nleaves
nn=os%nnode
ntest=10
! check number of particles to inject/remove to adjust the size of arrays
! in main program
allocate (ntot(ne),stat=err) ; if (err.ne.0) call stop_run ('error in allocating ntot in update_cloud$')
allocate (lev(ne),stat=err) ; if (err.ne.0) call stop_run ('error in allocating lev in update_cloud$')
ntot=0
do i=1,np
call octree_find_leaf (os%octree,os%noctree,cl%x(i),cl%y(i),cl%z(i), &
leaf,level,locc,x0,y0,z0,dxyz)
ntot(leaf)=ntot(leaf)+1
enddo
call octree_find_element_level (os%octree,os%noctree,lev,ne)
ninject=0
nremove=0
do ie=1,ne
if (ntot(ie).lt.params%npmin) ninject=ninject+(params%npmin-ntot(ie))
if (ntot(ie).gt.params%npmax .and. lev(ie).ge.params%levelmax_oct) &
nremove=nremove+(ntot(ie)-params%npmax)
enddo
deallocate (ntot)
deallocate (lev)
if (ninject.gt.0.and.iproc.eq.0) write (8,*) 'Injecting ',ninject,' particles'
if (nremove.gt.0.and.iproc.eq.0) write (8,*) 'Removing ',nremove,' particles'
npswap=cl%np
if (npswap.ne.0) then
allocate (memswap(11,npswap),stat=err) ; if (err.ne.0) call stop_run ('error in alloc memswap in main$')
allocate (imemswap(1,npswap),stat=err) ; if (err.ne.0) call stop_run ('error in alloc imemswap in main$')
memswap(1,:)=cl%x
memswap(2,:)=cl%y
memswap(3,:)=cl%z
memswap(4,:)=cl%x0
memswap(5,:)=cl%y0
memswap(6,:)=cl%z0
memswap(7,:)=cl%strain
memswap(8,:)=cl%lsf0
memswap(9,:)=cl%temp
memswap(10,:)=cl%press
memswap(11,:)=cl%e2dp
imemswap(1,:)=cl%tag
endif
deallocate (cl%x)
deallocate (cl%y)
deallocate (cl%z)
deallocate (cl%x0)
deallocate (cl%y0)
deallocate (cl%z0)
deallocate (cl%strain)
deallocate (cl%temp)
deallocate (cl%press)
deallocate (cl%lsf0)
deallocate (cl%e2dp)
deallocate (cl%tag)
npnew=cl%np+ninject-nremove
npnew=max(npnew,npswap)
allocate (cl%x(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%x in main$')
allocate (cl%y(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%y in main$')
allocate (cl%z(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%z in main$')
allocate (cl%x0(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%x0 in main$')
allocate (cl%y0(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%y0 in main$')
allocate (cl%z0(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%z0 in main$')
allocate (cl%strain(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%strain in main$')
allocate (cl%temp(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%temp in main$')
allocate (cl%press(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%press in main$')
allocate (cl%lsf0(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%lsf0 in main$')
allocate (cl%e2dp(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%e2dp in main$')
allocate (cl%tag(npnew),stat=err) ; if (err.ne.0) call stop_run ('error in alloc cl%tag in main$')
if (npswap.ne.0) then
cl%x(1:npswap)=memswap(1,:)
cl%y(1:npswap)=memswap(2,:)
cl%z(1:npswap)=memswap(3,:)
cl%x0(1:npswap)=memswap(4,:)
cl%y0(1:npswap)=memswap(5,:)
cl%z0(1:npswap)=memswap(6,:)
cl%strain(1:npswap)=memswap(7,:)
cl%lsf0(1:npswap)=memswap(8,:)
cl%temp(1:npswap)=memswap(9,:)
cl%press(1:npswap)=memswap(10,:)
cl%e2dp(1:npswap)=memswap(11,:)
cl%tag(1:npswap)=imemswap(1,:)
deallocate (memswap)
deallocate (imemswap)
else
cl%strain=0.d0
cl%lsf0=0.d0
cl%temp=0.d0
cl%press=0.d0
cl%e2dp=0.d0
cl%tag=1
endif
! first build location array
! ntot(ie) is number of particles in leaf ie
! nloc(ie) is location in array loc of where the particles belonging to ie are stored
! loc(nloc(ie):nloc(ie)+ntot(ie)-1) is list of particles belonging to ie
allocate (loc(np),stat=err) ; if (err.ne.0) call stop_run ('error in alloc loc in update_cloud$')
allocate (nloc(ne),stat=err) ; if (err.ne.0) call stop_run ('error in alloc nloc in update_cloud$')
allocate (ntot(ne),stat=err) ; if (err.ne.0) call stop_run ('error in alloc ntot in update_cloud$')
allocate (lev(ne),stat=err) ; if (err.ne.0) call stop_run ('error in alloc lev in update_cloud$')
allocate (remove(np),stat=err) ; if (err.ne.0) call stop_run ('error in alloc remove in update_cloud$')
ntot=0
do i=1,np
call octree_find_leaf (os%octree,os%noctree,cl%x(i),cl%y(i),cl%z(i), &
leaf,level,locc,x0,y0,z0,dxyz)
ntot(leaf)=ntot(leaf)+1
enddo
nloc(1)=1
do ie=2,ne
nloc(ie)=nloc(ie-1)+ntot(ie-1)
enddo
ntot=0
loc=0
do i=1,np
call octree_find_leaf (os%octree,os%noctree,cl%x(i),cl%y(i),cl%z(i), &
leaf,level,locc,x0,y0,z0,dxyz)
loc(nloc(leaf)+ntot(leaf))=i
ntot(leaf)=ntot(leaf)+1
enddo
call octree_find_element_level (os%octree,os%noctree,lev,ne)
! first check for injection
ninjectp=ninject
ninject=0
do ie=1,ne
if (ntot(ie).lt.params%npmin) then
ninject=ninject+(params%npmin-ntot(ie))
endif
enddo
allocate (xinject(ninject),stat=err) ; if (err.ne.0) call stop_run ('error in alloc xinject in update_cloud$')
allocate (yinject(ninject),stat=err) ; if (err.ne.0) call stop_run ('error in alloc yinject in update_cloud$')
allocate (zinject(ninject),stat=err) ; if (err.ne.0) call stop_run ('error in alloc zinject in update_cloud$')
allocate (straininject(ninject),stat=err) ; if (err.ne.0) call stop_run ('error in alloc straininject in update_cloud$')
allocate (x0inject(ninject),stat=err) ; if (err.ne.0) call stop_run ('error in alloc x0inject in update_cloud$')
allocate (y0inject(ninject),stat=err) ; if (err.ne.0) call stop_run ('error in alloc y0inject in update_cloud$')
allocate (z0inject(ninject),stat=err) ; if (err.ne.0) call stop_run ('error in alloc z0inject in update_cloud$')
allocate (e2dpinject(ninject),stat=err) ; if (err.ne.0) call stop_run ('error in alloc e2dpinject in update_cloud$')
if (ninject.ne.ninjectp) call stop_run ('error 1 in update_cloud$')
kinject=0
do ie=1,ne
if (ntot(ie).lt.params%npmin) then
kkinject=kinject+1
do ii=1,params%npmin-ntot(ie)
allocate (xtest(ntest),stat=err) ; if (err.ne.0) call stop_run ('error in alloc xtest in update_cloud$')
allocate (ytest(ntest),stat=err) ; if (err.ne.0) call stop_run ('error in alloc ytest in update_cloud$')
allocate (ztest(ntest),stat=err) ; if (err.ne.0) call stop_run ('error in alloc ztest in update_cloud$')
call random_number (xtest)
call random_number (ytest)
call random_number (ztest)
xtest=os%x(os%icon(1,ie))+xtest*(os%x(os%icon(2,ie))-os%x(os%icon(1,ie)))
ytest=os%y(os%icon(1,ie))+ytest*(os%y(os%icon(3,ie))-os%y(os%icon(1,ie)))
ztest=os%z(os%icon(1,ie))+ztest*(os%z(os%icon(5,ie))-os%z(os%icon(1,ie)))
distmax=0.d0
nchoice=0
do i=1,ntest
dist=0.d0
ndist=0
do k=1,ntot(ie)
dist=dist+(xtest(i)-cl%x(loc(nloc(ie)+k-1)))**2+ &
(ytest(i)-cl%y(loc(nloc(ie)+k-1)))**2+ &
(ztest(i)-cl%z(loc(nloc(ie)+k-1)))**2
ndist=ndist+1
enddo
do k=1,8
dist=dist+(xtest(i)-os%x(os%icon(k,ie)))**2+ &
(ytest(i)-os%y(os%icon(k,ie)))**2+ &
(ztest(i)-os%z(os%icon(k,ie)))**2
ndist=ndist+1
enddo
do k=kkinject,kinject
dist=dist+(xtest(i)-xinject(k))**2+ &
(ytest(i)-yinject(k))**2+ &
(ztest(i)-zinject(k))**2
ndist=ndist+1
enddo
dist=dist/ndist
if (dist.gt.distmax) then
distmax=dist
nchoice=i
endif
enddo
if (nchoice.eq.0) call stop_run ('error in finding inject part loc$')
! found the random particle that is the furthest away from the others and the
! 8 nodes of the leaf and previously injected ones; now we inject a particle
kinject=kinject+1
xinject(kinject)=xtest(nchoice)
yinject(kinject)=ytest(nchoice)
zinject(kinject)=ztest(nchoice)
deallocate (xtest)
deallocate (ytest)
deallocate (ztest)
enddo
endif
enddo
if (kinject.ne.ninjectp) call stop_run ('error 2 in update_cloud$')
! second check for removal
remove=0
do ie=1,ne
if (ntot(ie).gt.params%npmax.and.lev(ie).ge.params%levelmax_oct) then
do ii=1,ntot(ie)-params%npmax
distmin=1.d2
nchoice=0
do i=1,ntot(ie)
if (remove(loc(nloc(ie)+i-1)).eq.0) then
dist=0.d0
ndist=0
xx=cl%x(loc(nloc(ie)+i-1))
yy=cl%y(loc(nloc(ie)+i-1))
zz=cl%z(loc(nloc(ie)+i-1))
do k=1,ntot(ie)
if (remove(loc(nloc(ie)+k-1)).eq.0 .and. k.ne.i) then
dist=dist+(xx-cl%x(loc(nloc(ie)+k-1)))**2+ &
(yy-cl%y(loc(nloc(ie)+k-1)))**2+ &
(zz-cl%z(loc(nloc(ie)+k-1)))**2
ndist=ndist+1
endif
enddo
do k=1,8
dist=dist+(xx-os%x(os%icon(k,ie)))**2+ &
(yy-os%y(os%icon(k,ie)))**2+ &
(zz-os%z(os%icon(k,ie)))**2
ndist=ndist+1
enddo
dist=dist/ndist
if (dist.lt.distmin) then
distmin=dist
nchoice=loc(nloc(ie)+i-1)
endif
endif
enddo
if (nchoice.eq.0) call stop_run ('error in finding inject part loc$')
! found the particle that is the closest to the others and the
! 8 nodes of the leaf; we remove it
remove(nchoice)=1
enddo
endif
enddo
! now interpolate onto the nodes of the octree solve
os%strain=0.d0
os%e2dp=0.d0
allocate (strain(nn),stat=err) ; if (err.ne.0) call stop_run ('error in alloc strain in update_cloud$')
allocate (e2dp(nn),stat=err) ; if (err.ne.0) call stop_run ('error in alloc e2dp in update_cloud$')
allocate (w(nn),stat=err) ; if (err.ne.0) call stop_run ('error in alloc w in update_cloud')
strain=0.d0
e2dp=0.d0
w=0.d0
do i=1,np
call octree_find_leaf (os%octree,os%noctree,cl%x(i),cl%y(i),cl%z(i), &
leaf,level,locc,x0,y0,z0,dxyz)
r=(cl%x(i)-x0)/dxyz*2.d0-1.d0
s=(cl%y(i)-y0)/dxyz*2.d0-1.d0
t=(cl%z(i)-z0)/dxyz*2.d0-1.d0
h(1)=(1.d0-r)*(1.d0-s)*(1.d0-t)/8.d0
h(2)=(1.d0+r)*(1.d0-s)*(1.d0-t)/8.d0
h(3)=(1.d0-r)*(1.d0+s)*(1.d0-t)/8.d0
h(4)=(1.d0+r)*(1.d0+s)*(1.d0-t)/8.d0
h(5)=(1.d0-r)*(1.d0-s)*(1.d0+t)/8.d0
h(6)=(1.d0+r)*(1.d0-s)*(1.d0+t)/8.d0
h(7)=(1.d0-r)*(1.d0+s)*(1.d0+t)/8.d0
h(8)=(1.d0+r)*(1.d0+s)*(1.d0+t)/8.d0
do k=1,8
ic=os%icon(k,leaf)
ww=h(k)/dxyz**3
w(ic)=w(ic)+ww
strain(ic)=strain(ic)+ww*cl%strain(i)
e2dp(ic)=e2dp(ic)+ww*cl%e2dp(i)
enddo
enddo
if (np.gt.0) then
allocate (done(os%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc done in update_cloud$')
nisolated=0
done=0
do i=1,os%nnode
if (w(i).gt.0.d0) then
os%strain(i)=strain(i)/w(i)
os%e2dp(i)=e2dp(i)/w(i)
else
done(i)=1
nisolated=nisolated+1
endif
enddo
if (iproc.eq.0) write (8,*) 'There were ',nisolated,' octree nodes'
111 continue
strain=0.d0
e2dp=0.d0
w=0.d0
do ie=1,os%nleaves
do k=1,8
ic=os%icon(k,ie)
if (done(ic).eq.1) then
do kk=1,8
jc=os%icon(kk,ie)
if (done(jc).eq.0) then
ww=1.d0/dxyz**3
w(ic)=w(ic)+ww
strain(ic)=strain(ic)+os%strain(jc)*ww
e2dp(ic)=e2dp(ic)+os%e2dp(jc)*ww
endif
enddo
endif
enddo
enddo
ntodo=0
do i=1,os%nnode
if (done(i).eq.1) then
if (w(i).gt.0.d0) then
done(i)=0
os%strain(i)=strain(i)/w(i)
os%e2dp(i)=e2dp(i)/w(i)
else
ntodo=ntodo+1
endif
endif
enddo
if (iproc.eq.0) write (8,*) ntodo,' Isolated left '
if (ntodo.ne.0) goto 111
deallocate (done)
endif
deallocate (strain)
deallocate (e2dp)
deallocate (w)
! now calculate the field value at the new (injected points)
straininject=0.d0
e2dpinject=0.d0
x0inject=0.d0
y0inject=0.d0
z0inject=0.d0
do i=1,ninject
call octree_find_leaf (os%octree,os%noctree,xinject(i),yinject(i),zinject(i), &
leaf,level,locc,x0,y0,z0,dxyz)
r=(xinject(i)-x0)/dxyz*2.d0-1.d0
s=(yinject(i)-y0)/dxyz*2.d0-1.d0
t=(zinject(i)-z0)/dxyz*2.d0-1.d0
h(1)=(1.d0-r)*(1.d0-s)*(1.d0-t)/8.d0
h(2)=(1.d0+r)*(1.d0-s)*(1.d0-t)/8.d0
h(3)=(1.d0-r)*(1.d0+s)*(1.d0-t)/8.d0
h(4)=(1.d0+r)*(1.d0+s)*(1.d0-t)/8.d0
h(5)=(1.d0-r)*(1.d0-s)*(1.d0+t)/8.d0
h(6)=(1.d0+r)*(1.d0-s)*(1.d0+t)/8.d0
h(7)=(1.d0-r)*(1.d0+s)*(1.d0+t)/8.d0
h(8)=(1.d0+r)*(1.d0+s)*(1.d0+t)/8.d0
do k=1,8
straininject(i)=straininject(i)+os%strain(os%icon(k,leaf))*h(k)
e2dpinject(i)=e2dpinject(i)+os%e2dp(os%icon(k,leaf))*h(k)
x0inject(i)=x0inject(i)+os%x(os%icon(k,leaf))*h(k)
y0inject(i)=y0inject(i)+os%y(os%icon(k,leaf))*h(k)
z0inject(i)=z0inject(i)+os%z(os%icon(k,leaf))*h(k)
enddo
enddo
! modify the particle cloud for removal/injection of particles
nremovep=0
do i=1,np
if (remove(i).eq.1) then
nremovep=nremovep+1
else
cl%x(i-nremovep)=cl%x(i)
cl%y(i-nremovep)=cl%y(i)
cl%z(i-nremovep)=cl%z(i)
cl%x0(i-nremovep)=cl%x0(i)
cl%y0(i-nremovep)=cl%y0(i)
cl%z0(i-nremovep)=cl%z0(i)
cl%strain(i-nremovep)=cl%strain(i)
cl%lsf0(i-nremovep)=cl%lsf0(i)
cl%temp(i-nremovep)=cl%temp(i)
cl%press(i-nremovep)=cl%press(i)
cl%e2dp(i-nremovep)=cl%e2dp(i)
cl%tag(i-nremovep)=cl%tag(i)
endif
enddo
if (nremovep.ne.nremove) call stop_run ('error in update_cloud/remove$')
np=np-nremovep
do i=1,ninject
np=np+1
cl%x(np)=xinject(i)
cl%y(np)=yinject(i)
cl%z(np)=zinject(i)
cl%x0(np)=x0inject(i)
cl%y0(np)=y0inject(i)
cl%z0(np)=z0inject(i)
cl%strain(np)=straininject(i)
cl%lsf0(np)=0.d0
cl%temp(np)=0.d0
cl%press(np)=0.d0
cl%e2dp(np)=e2dpinject(i)
cl%ntag=cl%ntag+1
cl%tag(np)=cl%ntag
enddo
! in first step, we set the original particle positions to their
! current position
if (cl%np.eq.0) then
do i=1,np
cl%x0(i)=cl%x(i)
cl%y0(i)=cl%y(i)
cl%z0(i)=cl%z(i)
cl%strain(i)=0.d0
cl%temp(i)=0.d0
cl%press(i)=0.d0
cl%lsf0(i)=0.d0
cl%e2dp(i)=0.d0
cl%tag(i)=i
enddo
cl%ntag=np
endif
cl%np=np
deallocate (xinject)
deallocate (yinject)
deallocate (zinject)
deallocate (straininject)
deallocate (e2dpinject)
deallocate (x0inject,y0inject,z0inject)
deallocate (loc)
deallocate (nloc)
deallocate (ntot)
deallocate (remove)
deallocate (lev)
if (iproc.eq.0) then
write (8,*) 'Strain on osolve ',minval(os%strain),maxval(os%strain)
write (8,*) cl%np,' particles in 3D cloud'
end if
return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|