Newer
Older
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
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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! 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_many (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_many_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
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