!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! MOVE_SURFACE Nov. 2007 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine move_surface (surface,surface0,flag0,ov,dt,params,istep,is)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! subroutine to move a set of particles on a surface from a velocity
! field known at the nodes of an octree
! surface is the surface on which the particles need to be moved
! surface0 is the surface in its initial geometry
! flag0 is 1 if in midpoint configuration and 0 if at the end of the timestep
! if flag0=0, then the surface0 arrays are deallocated!!
! octree is the velocity octree
! noctree is its size
! unode,vnode and wnode are the components of the velocity
! nnode is the number of nodes on the octree
! icon_octree is the connectivity matrix on the octree
! nleaves is the number of leaves in the octree
! dt is the time step
! niter_move is the number of iterations used to move the particles
! it is read in input.txt
! here we use a simple method that performs a given (niter_move)
! number of iterations in a mid_point algorithm
! we also advect the normals by using the velocity gradient computed
! from the finite element shape function derivatives
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use definitions
implicit none
type (sheet) surface
type (sheet) surface0
type (octreev) ov
integer flag0
double precision dt
type(parameters) params
integer istep,is
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer i,ip,iter,nremove,nedge,nrem
integer i1,i2,i3,ij,k,err,nside
integer iproc,nproc,ierr
integer,dimension(:),allocatable::new,old
type (sheet) surfacep,surfacep0
double precision u,v,w,rat,eps
double precision x0,y0,z0,xi,yi,zi,x,y,z
double precision,dimension(:),allocatable::xbuf,ybuf,zbuf
double precision,dimension(:),allocatable::xnbuf,ynbuf,znbuf
double precision,dimension(:),allocatable::pos,zzz
double precision ux,uy,uz,vx,vy,vz,wx,wy,wz,xyzn
double precision xn0,yn0,zn0,xyn,theta,phi
double precision xn1,yn1,zn1,xn2,yn2,zn2
double precision xln1,yln1,zln1,xln2,yln2,zln2
double precision x1,x2,x3,y1,y2,y3,z1,z2,z3,xne,yne,zne
logical remove
character*72 :: shift
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
INCLUDE 'mpif.h'
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
shift=' '
allocate (xbuf(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc xbuf in move_surface$')
allocate (ybuf(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ybuf in move_surface$')
allocate (zbuf(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc zbuf in move_surface$')
allocate (xnbuf(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc xnbuf in move_surface$')
allocate (ynbuf(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ynbuf in move_surface$')
allocate (znbuf(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc znbuf in move_surface$')
xbuf=0.d0
ybuf=0.d0
zbuf=0.d0
xnbuf=0.d0
ynbuf=0.d0
znbuf=0.d0
rat=1.d0-1.d-10
eps=1.d-10
!nside=2**surface%levelt+1
nside=2**surface%levelt+2 !opla
do i=iproc+1,surface%nsurface,nproc
u=0.d0
v=0.d0
w=0.d0
x0=surface%x(i)
y0=surface%y(i)
z0=surface%z(i)
x=x0
y=y0
z=z0
do iter=1,params%niter_move
xi=(x0+x)/2.d0
yi=(y0+y)/2.d0
zi=(z0+z)/2.d0
call octree_interpolate_three (3,ov%octree,ov%noctree,ov%icon,ov%nleaves, &
ov%nnode,xi,yi,zi,ov%unode,u,ov%vnode,v,ov%wnode,w)
if (i==1 .or. i==nside .or. i==2*nside-1 .or. i==3*nside-2) then !corner points
z=z0+dt*w
elseif (i>1 .and. i<nside) then !y=0 edge
x=x0+dt*u
z=z0+dt*w
elseif (i>nside .and. i<2*nside-1) then !x=1 edge
y=y0+dt*v
z=z0+dt*w
elseif (i>2*nside-1 .and. i<3*nside-2) then !y=1 edge
x=x0+dt*u
z=z0+dt*w
elseif (i>3*nside-2 .and. i<4*nside-3) then !x=0 edge
y=y0+dt*v
z=z0+dt*w
else
x=x0+dt*u
y=y0+dt*v
z=z0+dt*w
end if
enddo
xbuf(i)=x
ybuf(i)=y
zbuf(i)=z
if (params%normaladvect) then
call octree_interpolate_three_derivative (3,ov%octree,ov%noctree,ov%icon,ov%nleaves, &
ov%nnode,xi,yi,zi, &
ov%unode,u,ux,uy,uz, &
ov%vnode,v,vx,vy,vz, &
ov%wnode,w,wx,wy,wz)
xn0=surface%xn(i)
yn0=surface%yn(i)
zn0=surface%zn(i)
xyn=xn0**2+yn0**2
if (xyn.gt.0.d0) xyn=sqrt(xyn)
theta=atan2(xyn,zn0)
phi=0.d0
if (sin(theta).ne.0.d0) phi=atan2(yn0,xn0)
xn1=cos(theta)*cos(phi)
yn1=cos(theta)*sin(phi)
zn1=-sin(theta)
xn2=-sin(phi)
yn2=cos(phi)
zn2=0.d0
xln1=ux*xn1+uy*yn1+uz*zn1
yln1=vx*xn1+vy*yn1+vz*zn1
zln1=wx*xn1+wy*yn1+wz*zn1
xln2=ux*xn2+uy*yn2+uz*zn2
yln2=vx*xn2+vy*yn2+vz*zn2
zln2=wx*xn2+wy*yn2+wz*zn2
xnbuf(i)=xn0+(yn1*zln2-zn1*yln2+yln1*zn2-zln1*yn2)*dt+(yln1*zln2-zln1*yln2)*dt**2
ynbuf(i)=yn0+(zn1*xln2-xn1*zln2+zln1*xn2-xln1*zn2)*dt+(zln1*xln2-xln1*zln2)*dt**2
znbuf(i)=zn0+(xn1*yln2-yn1*xln2+xln1*yn2-yln1*xn2)*dt+(xln1*yln2-yln1*xln2)*dt**2
if (xnbuf(i).gt.1.d50) write(*,*) 'on proc ', iproc, ', xnbuf is large: ', xnbuf(i)
if (ynbuf(i).gt.1.d50) write(*,*) 'on proc ', iproc, ', ynbuf is large: ', ynbuf(i)
if (znbuf(i).gt.1.d50) write(*,*) 'on proc ', iproc, ', znbuf is large: ', znbuf(i)
xyzn=sqrt(xnbuf(i)**2+ynbuf(i)**2+znbuf(i)**2)
xnbuf(i)=xnbuf(i)/xyzn
ynbuf(i)=ynbuf(i)/xyzn
znbuf(i)=znbuf(i)/xyzn
endif
enddo
surface%x=0.d0
surface%y=0.d0
surface%z=0.d0
call mpi_allreduce (xbuf,surface%x,surface%nsurface,mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
call mpi_allreduce (ybuf,surface%y,surface%nsurface,mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
call mpi_allreduce (zbuf,surface%z,surface%nsurface,mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
if (params%normaladvect) then
surface%xn=0.d0
surface%yn=0.d0
surface%zn=0.d0
call mpi_allreduce (xnbuf,surface%xn,surface%nsurface,mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
call mpi_allreduce (ynbuf,surface%yn,surface%nsurface,mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
call mpi_allreduce (znbuf,surface%zn,surface%nsurface,mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
else
call compute_normals (surface%nsurface,surface%x,surface%y,surface%z, &
surface%nt,surface%icon,surface%xn,surface%yn,surface%zn)
endif
deallocate (xbuf)
deallocate (ybuf)
deallocate (zbuf)
deallocate (xnbuf)
deallocate (ynbuf)
deallocate (znbuf)
! -------------------------------------------------------------------------------
! At this point, all points of the surface (and possibly the normal they carry)
! have been advected. A little bit of extra work is needed for the points on the
! convex hull (those at the intersection of the surface with the cube). On each
! edge of the square, there are nside points of coordinates (x_i,z(x_i)).
! Because of advection, the x_i are no more equidistant, which after several time
! steps can lead to one these points to be advected outside the cube and removed,
! which in turn leads to other problems. In order to avoid this problem, and also
! to keep the triangulation somewhat regular near the edges, the curve (x_i,z(x_i))
! is reinterpolated at equidistant X_i by means of the 'resample' routine in
! /RESAMPLE/libresample.a, based on linear interpolation (k=1) or spline
! interpolation (k=3,5). This procedure is carried out only at the end of a time
! step, i.e. when flag0=0.
allocate(pos(nside))
allocate(zzz(nside))
pos=surface%x(1:nside)
zzz=surface%z(1:nside)
call resample (pos,zzz,nside,1)
surface%x(1:nside)=pos
surface%z(1:nside)=zzz
pos=surface%y(nside:2*nside-1)
zzz=surface%z(nside:2*nside-1)
call resample (pos,zzz,nside,1)
surface%y(nside:2*nside-1)=pos
surface%z(nside:2*nside-1)=zzz
pos=surface%x(3*nside-2:2*nside-1:-1)
zzz=surface%z(3*nside-2:2*nside-1:-1)
call resample (pos,zzz,nside,1)
surface%x(3*nside-2:2*nside-1:-1)=pos
surface%z(3*nside-2:2*nside-1:-1)=zzz
pos(1)=surface%y(1)
zzz(1)=surface%z(1)
pos(2:nside)=surface%y(4*nside-4:3*nside-2:-1)
zzz(2:nside)=surface%z(4*nside-4:3*nside-2:-1)
call resample (pos,zzz,nside,1)
surface%y(4*nside-4:3*nside-2:-1)=pos(2:nside)
surface%z(4*nside-4:3*nside-2:-1)=zzz(2:nside)
deallocate(pos)
deallocate(zzz)
!-------------------------------------------------------------------------------
if (surface%itype.eq.0) then
surface%r=surface%x
surface%s=surface%y
endif
!-------------------------------------------------------------------------------
! We now check whether points have been advected outside the simulation domain
! (unit cube), and if so, we need to remove them. However, because the triangulation
! is only computed at startup and further updated, the removal of a point requires
! the triangulation to be sown back by hand. This will most probably lead to a
! non-Delaunay triangulation (locally), but this will be taken care of in the
! check_delaunay subroutine. The first 4*nside-4 points in the surface
! array correspond to those on the convex hull so that we need not test them.
nrem=0
remove=.true.
do while (remove)
remove=.false.
do i=4*nside-3,surface%nsurface
if (surface%x(i).lt.-eps .or. surface%x(i).gt.1.d0+eps .or. &
surface%y(i).lt.-eps .or. surface%y(i).gt.1.d0+eps .or. &
surface%z(i).lt.-eps .or. surface%z(i).gt.1.d0+eps) then
call remove_point(i,surface,flag0,surface0)
nrem=nrem+1
remove=.true.
exit
endif
end do
end do
if (iproc.eq.0 .and. nrem/=0) write(*,'(a,i4,a)') shift//'removing',nrem,' particle from surface'
call check_delaunay (params,surface,is,istep,'delaunay_aft_remove',0)
if (.not.params%normaladvect) then
call compute_normals (surface%nsurface,surface%x,surface%y,surface%z, &
surface%nt,surface%icon,surface%xn,surface%yn,surface%zn)
end if
if (flag0==0) then
nedge=surface%nsurface+surface%nt-1
call DoRuRe_surf_stats (params%doDoRuRe,istep,50,is,surface%nt,surface%nsurface,nedge,-nrem)
end if
if (iproc.eq.0) write (8,*) 'Removal complete'
end subroutine move_surface