File: N:\mfix\model\des\cfassign.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CFASSIGN                                                C
4     !                                                                      C
5     !  Purpose: Assign the necessary values for DEM computation. For       C
6     !           example:                                                   C
7     !     - assigning DEM boundaries from the values entered for           C
8     !       MFIX input in mfix.datat                                       C
9     !     - assigning DEM gravity vector from MFIX input.                  C
10     !     - calculating DTSOLID based on particle properties: spring       C
11     !       coefficient, damping factor & mass                             C
12     !                                                                      C
13     !  Reviewer: Sreekanth Pannala                        Date: 09-Nov-06  C
14     !                                                                      C
15     !  Reviewer: Rahul Garg                               Date: 25-Mar-14  C
16     !  Comments: Breaking this subroutine into several subroutines for DEM C
17     !            and PIC models                                            C
18     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
19     
20           SUBROUTINE CFASSIGN
21     
22     !-----------------------------------------------
23     ! Modules
24     !-----------------------------------------------
25           USE param1
26           USE constant, only : GRAVITY_X, GRAVITY_Y, GRAVITY_Z
27           USE discretelement
28           USE mfix_pic
29           use error_manager
30     ! Flag: DEM solids present.
31           use run, only: DEM_SOLIDS
32     ! Runtime flag specifying MPPIC solids
33           use run, only: PIC_SOLIDS
34     
35           IMPLICIT NONE
36     !-----------------------------------------------
37     ! Local Variables
38     !-----------------------------------------------
39     
40           CALL INIT_ERR_MSG("CFASSIGN")
41     
42     ! The common assignments are done in this routine.
43     ! The model spcific assignmets are moved to the specific subroutines
44     
45           GRAV(1) = GRAVITY_X
46           GRAV(2) = GRAVITY_Y
47           GRAV(3) = GRAVITY_Z
48           GRAV_MAG = sqrt(dot_product(GRAV,GRAV))
49     
50           IF(DEM_SOLIDS) CALL CFASSIGN_DEM
51           IF(PIC_SOLIDS) CALL CFASSIGN_PIC
52     
53     ! Finalize the error manager.
54           CALL FINL_ERR_MSG
55     
56           RETURN
57           END SUBROUTINE CFASSIGN
58     
59     
60     
61     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
62     !                                                                      C
63     !  Subroutine: CFASSIGN_PIC                                            C
64     !                                                                      C
65     !  Purpose: PIC related cfassign source code moved here                C
66     !           example:                                                   C
67     !     - calculating DTSOLID based on particle response time            C
68     !                                                                      C
69     !                                                                      C
70     !  Reviewer: Rahul Garg                               Date: 25-Mar-14  C
71     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
72     
73           SUBROUTINE CFASSIGN_PIC
74     
75     
76     !-----------------------------------------------
77     ! Modules
78     !-----------------------------------------------
79           USE discretelement, only: des_mmax, dtsolid
80           use param1, only: large_number
81           USE physprop, only: MU_g0, mmax, d_p0, ro_s0
82           USE mfix_pic, only : dtpic_taup, des_tau_p
83           USE mfix_pic, only: MPPIC_PDRAG_IMPLICIT
84           use error_manager
85           IMPLICIT NONE
86     !-----------------------------------------------
87     ! Local Variables
88     !-----------------------------------------------
89           INTEGER :: M, MMAX_TOT
90     !-----------------------------------------------
91     
92           CALL INIT_ERR_MSG("CFASSIGN_PIC")
93     
94           MMAX_TOT = DES_MMAX+MMAX
95           DO M = MMAX+1, MMAX_TOT
96     ! account for a possible offset index when using d_p0 and ro_s.
97              DES_TAU_P(M) = RO_S0(M)*(D_P0(M)**2.d0)/(18.d0*MU_g0)
98              WRITE(err_msg,'(/A,I2,A,G17.8)') &
99              'TAU_P FOR ', M,'th SOLID PHASE= ', DES_TAU_P(M)
100              CALL FLUSH_ERR_MSG (Header = .false., Footer = .false.)
101     
102           ENDDO
103     
104           DTSOLID = MINVAL(DES_TAU_P(MMAX+1:MMAX_TOT))
105           DTPIC_TAUP = DTSOLID      !maximum dt for point-particles based on taup
106           if(MPPIC_PDRAG_IMPLICIT) DTPIC_TAUP = LARGE_NUMBER
107     
108           WRITE(ERR_MSG, 1000) DTSolid
109           CALL FLUSH_ERR_MSG(Header = .false.)
110     
111      1000 format('MPPIC: Point-particle ',&
112           'approximation for particle-particle and ', /, &
113           'particle-wall collisions', /, &
114           'DTSOLID based on particle time response Taup', /, &
115           'DTSOLID = ', 2x, E17.10)
116     
117     ! Finalize the error manager.
118           CALL FINL_ERR_MSG
119     
120           END SUBROUTINE CFASSIGN_PIC
121     
122     
123     
124     
125     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
126     !                                                                      C
127     !  Subroutine: CFASSIGN_DEM                                            C
128     !                                                                      C
129     !  Purpose: DEM related cfassign source code moved here                C
130     !           example:                                                   C
131     !     - calculating DTSOLID based on particle properties: spring       C
132     !       coefficient, damping factor & mass                             C
133     !                                                                      C
134     !                                                                      C
135     !  Reviewer: Rahul Garg                               Date: 25-Mar-14  C
136     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
137     
138           SUBROUTINE CFASSIGN_DEM
139     
140     
141     !-----------------------------------------------
142     ! Modules
143     !-----------------------------------------------
144           USE param1
145           USE discretelement
146           use error_manager
147           IMPLICIT NONE
148     !-----------------------------------------------
149     ! Local Variables
150     !-----------------------------------------------
151     
152           CALL INIT_ERR_MSG("CFASSIGN_DEM")
153     
154     !      WRITE(ERR_MSG,'(A)') 'Setting collision model parameters for DEM'
155     !      CALL FLUSH_ERR_MSG (Footer = .false.)
156     
157     ! Finalize the error manager.
158           CALL FINL_ERR_MSG
159           END SUBROUTINE CFASSIGN_DEM
160     
161     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
162     ! subroutine: compute_volume_of_nodes                                      C
163     ! Author: Rahul Garg                                                       C
164     ! Dare: Dec 20, 2011                                                       C
165     ! Purpose:  Due to the staggered grid, the interpolaiton of mean fields    C
166     ! is always first done at the nodes (of the scalar cell) for any quantity. C
167     ! For a quantity like drag force or ep_s, one needs to assign a geometric  C
168     ! volume to a node. In the past this was done on-the-fly in drag_fgs.f.    C
169     ! VCELL was the variable that was used and it used to be incorrecty set to C
170     ! the volume of the scalar cell that the particle belongs to.              C
171     ! This will be ok for uniform grid and will not hold for non-uniform grids.C
172     ! In addition, at the edges (in 2-D, the node volume is only half of       C
173     ! the standard cell volume. And for corner (in 2-D), it is only one fourth.C
174     ! This was also done on-the-fly earlier                                    C
175     ! But again the volume of the cell was used, which was not correct         C
176     ! also not extendable to cut-cell. So this routine computes the geoemetric C
177     ! volume of the nodes.                                                     C
178     !                                                                          C
179     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
180     
181           SUBROUTINE compute_volume_of_nodes
182     
183     !-----------------------------------------------
184     ! Modules
185     !-----------------------------------------------
186           USE param
187           USE param1
188           USE parallel
189           USE physprop
190           USE fldvar
191           USE run
192           USE geometry
193           USE indices
194           USE bc
195           USE compar
196           USE sendrecv
197           USE discretelement
198           USE cutcell
199           USE functions
200           implicit none
201     !-----------------------------------------------
202     ! Local Variables
203     !-----------------------------------------------
204     ! ijk index of fluid grid and corresponding i, j, k indices
205           integer :: ijk, iraw, jraw, kraw
206     ! i, j, k and (i+1, j+1, k+1) indices corrected for any
207     ! cyclic ghost boundaries on fluid grid
208           INTEGER :: I, J, K, ip, jp, kp
209           integer :: ipjk, ijpk, ipjpk, ijkp, ipjkp, ijpkp, ipjpkp
210     ! volume of indicated grid
211           double precision :: vol_ijk, vol_ipjk, vol_ijpk, vol_ipjpk
212           double precision :: vol_ijkp, vol_ipjkp, vol_ijpkp, vol_ipjpkp
213           double precision :: vol_node_count, vol_node_actual_count
214     ! weighting factor
215           double precision :: avg_factor
216     ! not used?
217           double precision :: vol_node_uncorr
218     !-----------------------------------------------
219     
220           avg_factor = merge(0.25d0, 0.125d0, NO_K)
221     
222     ! compute the volume at the grid nodes
223     ! grid nodes start from istart2 to iend1
224           vol_node_count = merge(4., 8., NO_K)
225     
226           DO ijk = ijkstart3,ijkend3
227              des_vol_node(ijk) = zero
228              iraw = i_of(ijk)
229              jraw = j_of(ijk)
230              kraw = k_of(ijk)
231     
232     
233     
234     ! start at 1 (ghost cell) and go to last fluid cell. why start at a
235     ! ghost cell and not a fluid cell?
236     ! See below
237     
238     ! Since we are interested in nodes making up the interpolation stencil,
239     ! their numbering goes from 1 to iend1.
240     ! Think of a case with IMAX = 3. Here the nodes where the interpolation will be
241     ! done will run from 1 (=istart2) to 4 (=iend1).
242              IF(iraw.LT.istart2 .OR. iraw.GT.iend1) CYCLE
243              IF(jraw.LT.jstart2 .OR. jraw.GT.jend1) CYCLE
244              IF(kraw.LT.kstart2 .OR. kraw.GT.kend1) CYCLE
245     
246     ! this will force indices of ghost cells on cyclic borders back to
247     ! the corresponding fluid cells. since we are using i, j, k indices and
248     ! not just a composite ijk index we need these to be shifted to account
249     ! for periodicity
250              I = imap_c(iraw)
251              J = jmap_c(jraw)
252              K = kmap_c(kraw)
253              IP = imap_c(iraw+1)
254              JP = jmap_c(jraw+1)
255     
256     ! using a function like ip_of(ijk) should work the same as getting funijk
257     ! of the shifted i, j, k indices.  however, small differences will
258     ! occur on the 'edges/corners'. so retaining the latter method at this
259     ! time. see j. galvin for discussion
260              ipjk = funijk(IP,J,K)
261              ijpk = funijk(I,JP,K)
262              ipjpk = funijk(IP,JP,K)
263     
264     ! the existing variable vol(ijk) is not used here for cut-cell reasons
265     ! see r. garg for discussion
266              vol_ijk   = dx(I) *dy(J) *dz(K)
267              vol_ipjk  = dx(IP)*dy(J) *dz(K)
268              vol_ijpk  = dx(I) *dy(JP)*dz(K)
269              vol_ipjpk = dx(IP)*dy(JP)*dz(K)
270     
271              vol_node_uncorr = avg_factor*(vol_ijk + vol_ipjk + vol_ijpk + vol_ipjpk)
272              vol_node_actual_count = vol_node_count
273     
274              IF(.NOT.FLUID_AT(ijk)) THEN
275                 vol_node_actual_count = vol_node_actual_count - 1
276                 vol_ijk  = zero
277              ENDIF
278     
279              IF(.NOT.FLUID_AT(ipjk)) THEN
280                 vol_node_actual_count = vol_node_actual_count - 1
281                 vol_ipjk  = zero
282              ENDIF
283     
284              IF(.NOT.FLUID_AT(ijpk)) THEN
285                 vol_node_actual_count = vol_node_actual_count - 1
286                 vol_ijpk  = zero
287              ENDIF
288     
289              IF(.NOT.FLUID_AT(ipjpk)) THEN
290                 vol_node_actual_count = vol_node_actual_count - 1
291                 vol_ipjpk = zero
292              ENDIF
293     
294     ! this will have non-zero values for non-fluid cells at the
295     ! west/south/bottom borders but not for east/north/top borders?
296              des_vol_node(ijk) = avg_factor*(vol_ijk + vol_ipjk + &
297                 vol_ijpk + vol_ipjpk)
298     
299              IF (DO_K) THEN
300                 KP     = kmap_c(kraw+1)
301                 ijkp   = funijk(I, J, KP)
302                 ipjkp  = funijk(IP,J, KP)
303                 ijpkp  = funijk(I, JP,KP)
304                 ipjpkp = funijk(IP,JP,KP)
305     
306                 vol_ijkp   = dx(I) *dy(J) *dz(KP)
307                 vol_ipjkp  = dx(IP)*dy(J) *dz(KP)
308                 vol_ijpkp  = dx(I) *dy(JP)*dz(KP)
309                 vol_ipjpkp = dx(IP)*dy(JP)*dz(KP)
310     
311                 vol_node_uncorr = avg_factor*(vol_node_uncorr + vol_ijkp + &
312                    vol_ipjkp + vol_ijpkp + vol_ipjpkp)
313     
314                 IF(.NOT.FLUID_AT(ijkp)) THEN
315                    vol_node_actual_count = vol_node_actual_count - 1
316                    vol_ijkp   = zero
317                 ENDIF
318     
319                 IF(.NOT.FLUID_AT(ipjkp)) THEN
320                    vol_node_actual_count = vol_node_actual_count - 1
321                    vol_ipjkp  = zero
322                 ENDIF
323     
324                 IF(.NOT.FLUID_AT(ijpkp)) THEN
325                    vol_node_actual_count = vol_node_actual_count - 1
326                    vol_ijpkp  = zero
327                 ENDIF
328     
329                 IF(.NOT.FLUID_AT(ipjpkp)) THEN
330                    vol_node_actual_count = vol_node_actual_count - 1
331                    vol_ipjpkp = zero
332                 ENDIF
333     
334                 des_vol_node(ijk) = des_vol_node(ijk) + avg_factor*&
335                    (vol_ijkp + vol_ipjpkp + vol_ijpkp + vol_ipjkp)
336     
337              ENDIF
338     
339     
340           ENDDO   ! end do ijk=ijkstart3,ijkend3
341     
342           RETURN
343           RETURN
344           END SUBROUTINE compute_volume_of_nodes
345     
346