Newer
Older
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! MOVE_SURFACE Nov. 2007 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine move_surface (surface,surface0,flag0,flag1,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!!
! flag1 is 1 if surface is only being shifted vertically due to isostasy
! 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
!use mpi
include 'mpif.h'
type (sheet) surface
type (sheet) surface0
type (octreev) ov
integer flag0,flag1
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
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
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
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
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
if (flag1) then
! Only vertical velocity from isostasy
call octree_interpolate_three (3,ov%octree,ov%noctree,ov%icon,ov%nleaves, &
ov%nnode,xi,yi,zi,0.d0*ov%unode,u,0.d0*ov%vnode,v,ov%wnodeiso,w)
else
! Full velocity field
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/params%vex,w)
endif
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
if (flag1) then
! Only vertical velocity from isostasy
call octree_interpolate_three_derivative (3,ov%octree,ov%noctree,ov%icon,ov%nleaves, &
ov%nnode,xi,yi,zi, &
0.d0*ov%unode,u,ux,uy,uz, &
0.d0*ov%vnode,v,vx,vy,vz, &
ov%wnodeiso,w,wx,wy,wz)
else
! Full velocity field
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/params%vex,w,wx,wy,wz)
endif
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)
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
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