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