source: trunk/00FlowSolve_KW/src/immersed_boundary.f90 @ 2

Last change on this file since 2 was 2, checked in by xlvlod, 18 years ago

initial import from /home2/xlvlod/IDRIS/SVN_BASE_TRUNK
File size: 13.4 KB
Line 
1
2
3!!=====================================================================
4  module immersed_boundary 
5!!====================================================================
6   implicit none
7   
8   public  :: InitializeImmersedBoundary
9   public  :: IB_SourceStrength
10   public  :: IB_Update
11   public  :: ForcingFilter
12   public  :: ij_from_count
13   public  :: count_from_ij
14   public  :: unit_normal
15   
16   real(kind=8)                  :: eps(3),d_inf(3),delta_min
17   real(kind=8),allocatable      :: x_IB(:,:),y_IB(:,:),h_IB(:,:)
18   real(kind=8),allocatable      :: nhat(:,:,:),SS(:,:),f_IB(:,:,:)
19   real(kind=8),allocatable      :: IB_Filter(:,:,:)
20   integer                       :: n_i,n_j,n_k
21   integer,allocatable           :: labels(:,:),kplus(:,:)
22
23  Contains 
24 
25!===================================================================
26!===================================================================
27subroutine InitializeImmersedBoundary
28 use class_DataDecomposition
29 use independent_variables
30 use dependent_variables
31 use intermediate_variables
32 use user_data
33 use user_routines
34 use etc
35 use mpi_params
36 
37 implicit none
38 integer          :: idim,locnz
39 real(kind=8)     :: x(3),ds(3),h_x,h_y
40 real(kind=8)     :: s1,s2,s3,ans(3)
41 real(kind=8)     :: x_s,y_s,z_s
42 real(kind=8)     :: dx,dy,dz,xx,acosh2
43 real(kind=8),allocatable  :: work(:)
44 external            acosh2
45 
46 if( .not. do_immersed_boundary ) return
47       
48  !! get 'ds' in each of the x y z directions ==> ds(1-3)
49  if(ny==1) then
50   ds(2)=0.d0
51  else
52   ds(2)=Grid(0)%SpaceDims(2)%s(2)-Grid(0)%SpaceDims(2)%s(1)
53  endif
54  ds(1)=Grid(0)%SpaceDims(1)%s(2)-Grid(0)%SpaceDims(1)%s(1) 
55  ds(3)=Grid(0)%SpaceDims(3)%s(2)-Grid(0)%SpaceDims(3)%s(1)
56  locnz = ubound(velocity(1)%SlabVals,3)
57 
58  allocate( x_IB(nx,ny),y_IB(nx,ny),h_IB(nx,ny) )
59  allocate( nhat(nx,ny,3),IB_Filter(nx,ny,locnz) )
60  allocate( f_IB(nx,ny,3) )
61  allocate( SS(nx,ny))
62  allocate( labels(nx,ny),kplus(nx,ny) )
63  labels(:,:)=1
64 
65 
66 
67  do j=1,ny
68   do i=1,nx
69   
70   x(1)=Grid(0)%SpaceDims(1)%x(i)
71   x(2)=Grid(0)%SpaceDims(2)%x(j)
72   
73   call eval_surface(x,x(3),h_x,h_y)
74   call unit_normal(x,ans)
75   
76   !! x,y,z on IB surface
77   x_IB(i,j)=x(1)
78   y_IB(i,j)=x(2)
79   h_IB(i,j)=x(3)
80
81   !! unit normal vector
82   nhat(i,j,1)=ans(1)
83   nhat(i,j,2)=ans(2)
84   nhat(i,j,3)=ans(3)
85     
86   !! find first globally indexed k value ABOVE surface
87   do k=1,nz
88    kplus(i,j)=k
89    if( Grid(0)%SpaceDims(3)%x(k) > h_IB(i,j) ) goto 99
90   enddo
9199 continue
92   
93   
94   if( i==1 .and. j==1 ) then
95   
96   !! Get the param values s in each direction -> (s1,s2,s3)
97   !! corresponding too the surface point x(1:3).
98   !! Remember, this user routine has DIMENSIONAL inputs
99   call point_map_x2s(x(1)*Lz,1,xdim_periodic,Lx,nx,s1)
100   call point_map_x2s(x(2)*Lz,2,ydim_periodic,Ly,ny,s2)
101   call point_map_x2s(x(3)*Lz,3,zdim_periodic,Lz,nz,s3)
102   
103   !! now get the corresponding dx/ds,dy/ds,dz/ds values -> (x_s,y_s,z_s)
104   !! Remember, this user routine has DIMENSIONAL outputs
105   idim=1
106   call point_map_s2x(s1,idim,xdim_periodic,Lx,nx,xx,x_s)
107   
108   idim=2
109   call point_map_s2x(s2,idim,ydim_periodic,Ly,ny,xx,y_s)
110   
111   idim=3
112   call point_map_s2x(s3,idim,zdim_periodic,Lz,nz,xx,z_s)
113   
114   dx = (x_s/Lz) * ds(1)   !! local dxds * ds
115   dy = (y_s/Lz) * ds(2)   !! local dyds * ds
116   dz = (z_s/Lz) * ds(3)   !! local dzds * ds
117   if(ny==1) dy=1.d0
118     
119    delta_min = 1.e-6
120    eps(1) = 3.*dx
121    eps(2) = 3.*dy
122    eps(3) = 3.*dz
123   
124    d_inf(:) =  4.  !acosh2( 1./sqrt(2.*eps(:)*delta_min) )
125    n_i = eps(1)*d_inf(1)/dx
126    n_j = eps(2)*d_inf(2)/dy
127    n_k = eps(3)*d_inf(3)/dz
128   endif   
129             
130  enddo
131 enddo
132 
133 !! construct spatial filter so user forcing can be specified
134 !! more easily, i.e. globally w/o worrying about IB location
135 call ForcingFilter
136     
137 return
138end subroutine InitializeImmersedBoundary
139!=================================================================
140!=================================================================
141
142
143 
144  !!==============================================================
145  !!==============================================================
146  subroutine IB_SourceStrength
147   use user_data,  only: nx,ny,dt
148   use dependent_variables
149   use intermediate_variables
150   use etc, only: istep,N,MM0,M_oldest
151   implicit none
152   integer           :: i,j,k,m(3),npts
153   real(kind=8)      :: dtinv
154
155   npts=nx*ny
156   m(:)=0
157   dtinv=1.d0/dt
158   
159                     
160   !!====================================================
161   !! do the 3 interpolations
162   !!====================================================
163   do k=1,3                                         
164                                         
165    call spline_interp_value( x_IB,                 &    !x'
166                              y_IB,                 &    !y'
167                              h_IB,                 &    !z' on surface
168                              npts,velocity(k),m,   &    !field to interp
169                              f_IB(1,1,k),          &    ! ans
170                              labels,               &    !==1
171                              scratch(2)%SlabVals,  &    !work
172                              scratch(3)%SlabVals,  &    !work
173                              scratch(4)%SlabVals )      !work
174   enddo
175   
176   !! now fill in the SS w/ dot product with nhat
177   SS(:,:) = (f_IB(:,:,1)*nhat(:,:,1)  &
178           +  f_IB(:,:,2)*nhat(:,:,2)  &
179           +  f_IB(:,:,3)*nhat(:,:,3))*dtinv
180   
181   return
182  end subroutine IB_SourceStrength
183  !!====================================================================
184  !!====================================================================
185 
186 
187 
188  !!====================================================================
189  !!====================================================================
190   subroutine IB_Update
191    use user_data,              only: nx,ny,dt
192    use etc,                    only: numzplanes_below,N,istep                                 
193    use mpi_params,             only: myid
194    use independent_variables
195    use dependent_variables
196    implicit none
197    integer           :: i,j,k,ii,jj,kk,idim
198    integer,save      :: locnz,npts,k0
199    integer           :: istart,iend
200    integer           :: jstart,jend
201    real(kind=8),save :: ds(3)
202    real(kind=8)      :: x(3),dxds,dyds,dzds
203    real(kind=8)      :: dx,dy,dz,delta
204    real(kind=8)      :: d1,d2,d3,xx
205    logical,save      :: first_entry=.TRUE.
206    external          delta
207   
208  !!====================================================================
209  !!====================================================================
210   if( first_entry ) then
211    locnz=ubound(velocity(1)%SlabVals,3)
212    npts=nx*ny
213    k0=numzplanes_below(myid)
214    if(ny==1) then
215     ds(2)=0.d0
216    else
217     ds(2)=Grid(0)%SpaceDims(2)%s(2)-Grid(0)%SpaceDims(2)%s(1)
218    endif
219    ds(1)=Grid(0)%SpaceDims(1)%s(2)-Grid(0)%SpaceDims(1)%s(1) 
220    ds(3)=Grid(0)%SpaceDims(3)%s(2)-Grid(0)%SpaceDims(3)%s(1)
221    first_entry=.FALSE.
222   endif
223  !!====================================================================
224  !!====================================================================
225   
226    !! loop over "all" the grid points on myid (ii,jj,kk)
227    do kk=1,locnz
228     x(3)=Grid(0)%SpaceDims(3)%x(k0+kk) !; write(0,*) 'K = ',kk
229      do jj=1,ny
230       x(2)=Grid(0)%SpaceDims(2)%x(jj) 
231        do ii=1,nx       
232         x(1)=Grid(0)%SpaceDims(1)%x(ii)
233   
234         !!========================================================
235         !! if x(ii,jj,kk) is 'near' IB surface, compute integral,
236         !! otherwise, if x is too far from surface, do nothing
237         !!========================================================
238         if( abs( ( x(3) - h_IB(ii,jj) ) ) < d_inf(3)*eps(3) ) then
239   
240          !!================================================
241          !! above/below x(ii,jj,kk), define a range of i,j
242          !! values where d1 and d2 could be nonzero
243          !!================================================
244          istart = maxval( (/1,  ii - n_i/) )
245          iend   = minval( (/nx, ii + n_i/) )
246   
247          jstart = maxval( (/1,  jj - n_j/) )
248          jend   = minval( (/ny, jj + n_j/) )
249   
250          !!================================================
251          !! loop over this portion of the surface
252          !! and accumulate any contributions to F(x)
253          !!================================================
254          do j=jstart,jend
255           dyds= Grid(0)%SpaceDims(2)%dxds(j) 
256            dy = dyds*ds(2)
257            if(ny==1) dy=1.
258     
259            do i=istart,iend
260             dxds= Grid(0)%SpaceDims(1)%dxds(i) 
261             dx = dxds*ds(1)
262       
263              k=kplus(i,j)
264              dzds= Grid(0)%SpaceDims(3)%dxds(k) 
265              dz = dzds*ds(3)
266       
267              d1 = delta( (x(1)-x_IB(i,j)),eps(1) )
268              d2 = delta( (x(2)-y_IB(i,j)),eps(2) )
269              d3 = delta( (x(3)-h_IB(i,j)),eps(3) )
270              xx = SS(i,j)*d1*d2*d3*dx*dy*dz
271       
272              !! acumulate terms for the "idim" velocity component
273              do idim=1,3
274               velocity(idim)%SlabVals(ii,jj,kk)                         &
275                                 =  velocity(idim)%SlabVals(ii,jj,kk)    &
276                                 -  nhat(i,j,idim)*xx*dt
277              enddo
278       
279            enddo
280          enddo     !! end loop over section of surface for x(i,j)
281         
282         endif  !!  end process or no process for x(ii,jj,kk)
283         
284       enddo    !! next ii
285      enddo     !! next jj     
286     enddo      !! next kk, i.e. end loop over grid points on myid
287       
288   return
289  end subroutine IB_Update
290  !!====================================================================
291  !!====================================================================
292 
293  !!====================================================================
294  !!====================================================================
295   subroutine ForcingFilter
296    use user_data,              only: nx,ny
297    use etc,                    only: numzplanes_below                                 
298    use mpi_params,             only: myid
299    use intermediate_variables, only: scratch
300    use independent_variables,  only: Grid
301    implicit none
302    integer           :: i,j,k
303    integer,save      :: locnz,k0
304    real(kind=8)      :: x(3),step,delta,XX,YY,ZZ
305    external             step,delta
306   
307  !!====================================================================
308  !!====================================================================
309    locnz=ubound(scratch(1)%SlabVals,3)
310    k0=numzplanes_below(myid)
311   
312    !! loop over "all" the grid points on myid (i,j,k)
313    do k=1,locnz
314     x(3)=Grid(0)%SpaceDims(3)%x(k0+k) 
315      do j=1,ny
316       x(2)=Grid(0)%SpaceDims(2)%x(j) 
317        do i=1,nx       
318         x(1)=Grid(0)%SpaceDims(1)%x(i)
319   
320         
321          !! assuming the vector from x(:), that intersects
322          !! the IB surface normally, is not too far from
323          !! the surface point just above or below x(:).
324          !! Use that point and it's unit normal to to get
325          !! the x,y and z distances from x(:) to the
326          !! actual closest point on IB to x(:).
327          XX = (x(3) - h_IB(i,j))*nhat(i,j,1)/(3.*eps(1))
328          YY = (x(3) - h_IB(i,j))*nhat(i,j,2)/(3.*eps(2))
329          ZZ = (x(3) - h_IB(i,j))*nhat(i,j,3)/(3.*eps(3))
330         
331          if(ZZ>=0.d0)  then
332           IB_Filter(i,j,k) = 1.d0 
333          else
334           IB_Filter(i,j,k) = exp(-XX**2) * exp(-YY**2) * exp(-ZZ**2)
335          endif
336       
337                           
338       enddo   
339      enddo       
340     enddo   
341     
342   
343       
344   return
345  end subroutine ForcingFilter
346  !!====================================================================
347  !!====================================================================
348   
349
350  !!====================================================================
351  !!====================================================================
352   subroutine ij_from_count(i,j,count)
353    use user_data, only: nx
354    implicit none
355    integer   :: i,j,count
356    j = (count-1)/nx + 1
357    i = count - (j-1)/nx
358    return
359   end subroutine ij_from_count
360  !!====================================================================
361  !!====================================================================
362   
363  !!====================================================================
364  !!====================================================================   
365   subroutine count_from_ij(i,j,count)
366    use user_data, only: nx
367    implicit none
368    integer   :: i,j,count
369    count = (j-1)/nx + i
370    return
371   end subroutine count_from_ij
372  !!====================================================================
373  !!====================================================================
374 
375 
376!======================================================================
377!======================================================================
378subroutine unit_normal(x,n)
379 use user_routines, only: eval_surface
380 implicit none
381 real(kind=8) :: x(2),n(3),h,h_x,h_y
382 real(kind=8) :: mag
383 
384   call eval_surface(x,h,h_x,h_y)
385   mag = sqrt( h_x**2 + h_y**2 + 1. )
386   n(1)=-h_x/mag
387   n(2)=-h_y/mag
388   n(3)=1.d0/mag
389   
390 return
391end subroutine unit_normal  ! end of function unit_normal
392!======================================================================
393!======================================================================
394
395
396
397 
398!!==================================================================================== 
399  end module immersed_boundary 
400!!====================================================================================
Note: See TracBrowser for help on using the repository browser.