File: N:\mfix\model\usr_src_mod.f

1     ! or put the actual routine calc_usr_source here? move logic variable
2     !  usr_source into run or something...
3           module usr_src
4           use param, only: dim_eqs
5     
6     ! If .TRUE. call user-defined source subroutines
7           LOGICAL :: CALL_USR_SOURCE(DIM_EQS)
8     
9           ENUM, BIND(C)
10              ENUMERATOR :: Pressure_Correction, Solids_Correction
11              ENUMERATOR :: Gas_Continuity, Solids_Continuity
12              ENUMERATOR :: Gas_U_Mom, Solids_U_Mom, Gas_V_Mom, Solids_V_Mom
13              ENUMERATOR :: Gas_W_Mom, Solids_W_Mom
14              ENUMERATOR :: Gas_Energy, Solids_Energy
15              ENUMERATOR :: Gas_Species, Solids_Species
16              ENUMERATOR :: Gran_Energy
17              ENUMERATOR :: Usr_Scalar, K_Epsilon_K, K_Epsilon_E
18              ENUMERATOR :: BLANK
19              ENUMERATOR :: DUMMY
20           END ENUM
21     
22           CONTAINS
23     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
24     !                                                                      C
25     !  Subroutine: CALC_USR_SOURCE                                         C
26     !  Purpose: Driver routine to calculate user defined source terms      C
27     !  for each of the various equations                                   C
28     !                                                                      C
29     !                                                                      C
30     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
31           SUBROUTINE CALC_USR_SOURCE(lEQ_NO, A_M, B_M, lB_MMAX, lM, lN)
32     
33     ! Modules
34     !-----------------------------------------------
35           use compar, only: ijkstart3, ijkend3
36           use fldvar, only: ep_s, ep_g
37           use fun_avg, only: avg_x, avg_y, avg_z
38           use functions, only: fluid_at, wall_at
39           use functions, only: ip_at_e, ip_at_n, ip_at_t
40           use functions, only: sip_at_e, sip_at_n, sip_at_t
41           use functions, only: east_of, north_of, top_of
42           use indices, only: i_of, j_of, k_of
43           use param, only: dimension_3, dimension_m
44           use param1, only: undefined_i, zero
45           use pgcor, only: phase_4_p_g
46           use physprop, only: smax, mmax
47           use pscor, only: phase_4_p_s
48           use run, only: kt_type_enum, ghd_2007
49           use run, only: momentum_x_eq, momentum_y_eq, momentum_z_eq
50           use toleranc, only: dil_ep_s
51           use error_manager
52           IMPLICIT NONE
53     ! Dummy arguments
54     !-----------------------------------------------
55     ! Equation number
56     ! Table relating the eq_no set in the mfix.dat to the corresponding
57     ! pde/variable of interest
58     !     1, PP_g (pressure correction eq)
59     !     2, EPP (solids correction eq)
60     !     2, rop_g (gas continuity eq)
61     !     2, rop_s (solids continuity eq)
62     !     3, u_g (gas u-momentum eq)
63     !     3, u_s (solids u-momentum eq)
64     !     4, v_g (gas v-momentum eq)
65     !     4, v_s (solids v-momentum eq)
66     !     5, w_g (gas w-momentum eq)
67     !     5, w_s (solids w-momentum eq)
68     !     6, T_g (gas energy eq)
69     !     6, T_s (solids energy eq)
70     !     7, x_g (gas species eq)
71     !     7, x_s (solids species eq)
72     !     8, theta_m (granular energy eq)
73     !     9, scalar_eq (scalar eq)
74     !     9, k_epsilon k (k_epsilon eqs)
75     !     9, k_epsilon e (k_epsilon_eqs)
76           INTEGER, INTENT(IN) :: lEQ_NO
77     
78     ! Septadiagonal matrix A_m
79           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
80     ! Vector b_m
81           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
82     ! Vector b_mmax
83           DOUBLE PRECISION, OPTIONAL, INTENT(INOUT) :: lB_mmax(DIMENSION_3, 0:DIMENSION_M)
84     ! Phase index
85           INTEGER, OPTIONAL, INTENT(IN) :: lM
86     ! Species index OR scalar equation number
87           INTEGER, OPTIONAL, INTENT(IN) :: lN
88     
89     ! Local variables
90     !-----------------------------------------------
91     ! source terms which appear appear in the
92     ! center coefficient (lhs) - part of a_m matrix
93     ! source vector (rhs) - part of b_m vector
94           DOUBLE PRECISION :: sourcelhs, sourcerhs
95     ! indices
96           INTEGER :: IJK, I, J, K, IJKE, IJKN, IJKT
97     ! tmp indices
98           INTEGER :: L, ll, M, N
99     ! volume fractions
100           DOUBLE PRECISION :: EPGA, EPSA, EPtmp
101     !-----------------------------------------------
102     
103     ! set local values for phase index and species or scalar index
104           IF (.NOT.present(lM)) THEN
105              M = UNDEFINED_I
106           ELSE
107              M = lM
108           ENDIF
109           IF (.NOT.present(lN)) THEN
110              N = UNDEFINED_I
111           ELSE
112              N = lN
113           ENDIF
114     
115           SELECT CASE(lEQ_NO)
116     
117     ! gas pressure correction equation
118           CASE(Pressure_Correction)
119              DO IJK=IJKSTART3,IJKEND3
120                 IF (FLUID_AT(IJK)) THEN
121                    CALL USR_SOURCES(lEQ_NO, IJK, sourcelhs, sourcerhs, M, N)
122                    A_M(IJK,0,M) = A_M(IJK,0,M) - sourcelhs
123                    B_M(IJK,M) = B_M(IJK,M) - sourcerhs
124                    lB_MMAX(IJK,M) = max(abs(lB_MMAX(IJK,M)), abs(B_M(IJK,M)))
125                 ENDIF
126              ENDDO
127     
128     ! solids correction currently only called when smax==1 and
129     ! mcp/=undefined_i
130           CASE(Solids_correction)
131              DO IJK=IJKSTART3,IJKEND3
132                 IF (FLUID_AT(IJK)) THEN
133                    CALL USR_SOURCES(lEQ_NO, IJK, sourcelhs, sourcerhs, M, N)
134                    A_M(IJK,0,0) = A_M(IJK,0,0) - sourcelhs
135                    B_M(IJK,0) = B_M(IJK,0) - sourcerhs
136                    lB_MMAX(IJK,0) = max(abs(lB_MMAX(IJK,0)), abs(B_M(IJK,0)))
137                 ENDIF
138              ENDDO
139     
140     ! gas continuity
141           CASE(Gas_continuity)
142              DO IJK=IJKSTART3,IJKEND3
143                 IF (FLUID_AT(IJK)) THEN
144     !                            .AND. PHASE_4_P_G(IJK)/=M) THEN
145     ! have not included phase_4_p_g logic...which would require unique
146     ! structure here
147                    CALL USR_SOURCES(lEQ_NO, IJK, sourcelhs, sourcerhs, M, N)
148                    A_M(IJK,0,M) = A_M(IJK,0,M) - sourcelhs
149                    B_M(IJK,M) = B_M(IJK,M) - sourcerhs
150                 ENDIF
151              ENDDO
152     
153     ! soldis continuity
154           CASE(Solids_continuity)
155     ! have not included phase_4_p_g/p_s logic...which would require unique
156     ! structure here
157              DO IJK=IJKSTART3,IJKEND3
158                 IF (FLUID_AT(IJK)) THEN
159     !               .AND. PHASE_4_P_G(IJK)/=M .AND. PHASE_4_P_S(IJK)/=M) THEN
160                    CALL USR_SOURCES(lEQ_NO, IJK, sourcelhs, sourcerhs, M, N)
161                    A_M(IJK,0,M) = A_M(IJK,0,M) - sourcelhs
162                    B_M(IJK,M) = B_M(IJK,M) - sourcerhs
163                 ENDIF
164              ENDDO
165     
166     ! u-momentum equations
167           CASE(Gas_U_Mom)
168              M=0
169              IF(.NOT.MOMENTUM_X_EQ(M)) RETURN
170              DO IJK=IJKSTART3,IJKEND3
171                 IF (.NOT.FLUID_AT(IJK)) CYCLE
172                 I = I_OF(IJK)
173                 IJKE = EAST_OF(IJK)
174                 EPGA = AVG_X(EP_G(IJK),EP_G(IJKE),I)
175                 IF (IP_AT_E(IJK) .OR. EPGA <= DIL_EP_S) CYCLE
176                 CALL USR_SOURCES(lEQ_NO, IJK, sourcelhs, sourcerhs, M, N)
177                 A_M(IJK,0,M) = A_M(IJK,0,M) - sourcelhs
178                 B_M(IJK,M) = B_M(IJK,M) - sourcerhs
179              ENDDO
180     
181     ! u-momentum equations
182           CASE(Solids_U_Mom)
183              DO L=1,MMAX
184                 IF (KT_TYPE_ENUM /= GHD_2007 .OR. &
185                    (KT_TYPE_ENUM == GHD_2007 .AND. L==MMAX)) THEN
186                    IF(.NOT.MOMENTUM_X_EQ(L)) RETURN
187                    DO IJK=IJKSTART3,IJKEND3
188                       IF (.NOT.FLUID_AT(IJK)) CYCLE
189                       I = I_OF(IJK)
190                       IJKE = EAST_OF(IJK)
191                       EPtmp = ZERO
192                       IF (KT_TYPE_ENUM == GHD_2007) THEN
193                          DO lL = 1, SMAX
194                             EPtmp = EPtmp + AVG_X(EP_S(IJK,lL),EP_S(IJKE,lL),I)
195                          ENDDO
196                          EPSA = EPtmp
197                       ELSE
198                          EPSA = AVG_X(EP_S(IJK,L),EP_S(IJKE,L),I)
199                       ENDIF
200                       IF (IP_AT_E(IJK) .OR. SIP_AT_E(IJK) .OR. &
201                           EPSA <= DIL_EP_S) CYCLE
202                       CALL USR_SOURCES(lEQ_NO, IJK, sourcelhs, sourcerhs, L, N)
203                       A_M(IJK,0,L) = A_M(IJK,0,L) - sourcelhs
204                       B_M(IJK,L) = B_M(IJK,L) - sourcerhs
205                    ENDDO    ! enddo ijk
206                 ENDIF
207              ENDDO   ! enddo mmax
208     
209     ! v-momentum equations
210           CASE(Gas_V_Mom)
211              M=0
212              IF(.NOT.MOMENTUM_Y_EQ(M)) RETURN
213              DO IJK=IJKSTART3,IJKEND3
214                 IF (.NOT.FLUID_AT(IJK)) CYCLE
215                 J = J_OF(IJK)
216                 IJKN = NORTH_OF(IJK)
217                 EPGA = AVG_Y(EP_G(IJK),EP_G(IJKN),J)
218                 IF (IP_AT_N(IJK) .OR. EPGA <= DIL_EP_S) CYCLE
219                 CALL USR_SOURCES(lEQ_NO, IJK, sourcelhs, sourcerhs, M, N)
220                 A_M(IJK,0,M) = A_M(IJK,0,M) - sourcelhs
221                 B_M(IJK,M) = B_M(IJK,M) - sourcerhs
222              ENDDO
223     
224     ! v-momentum equations
225           CASE(Solids_V_Mom)
226              DO L=1,MMAX
227                 IF (KT_TYPE_ENUM /= GHD_2007 .OR. &
228                     (KT_TYPE_ENUM == GHD_2007 .AND. L==MMAX)) THEN
229                    IF(.NOT.MOMENTUM_Y_EQ(L)) RETURN
230                    DO IJK=IJKSTART3,IJKEND3
231                       IF (.NOT.FLUID_AT(IJK)) CYCLE
232                       J = J_OF(IJK)
233                       IJKN = NORTH_OF(IJK)
234                       EPtmp = ZERO
235                       IF (KT_TYPE_ENUM == GHD_2007) THEN
236                          DO lL = 1, SMAX
237                             EPtmp = EPtmp + AVG_Y(EP_S(IJK,lL),EP_S(IJKN,lL),J)
238                          ENDDO
239                          EPSA = EPtmp
240                       ELSE
241                          EPSA = AVG_Y(EP_S(IJK,L),EP_S(IJKN,L),J)
242                       ENDIF
243                       IF (IP_AT_N(IJK) .OR. SIP_AT_N(IJK) .OR. &
244                           EPSA <= DIL_EP_S) CYCLE
245                       CALL USR_SOURCES(lEQ_NO, IJK, sourcelhs, sourcerhs, L, N)
246                       A_M(IJK,0,L) = A_M(IJK,0,L) - sourcelhs
247                       B_M(IJK,L) = B_M(IJK,L) - sourcerhs
248                    ENDDO   ! enddo ijk
249                 ENDIF
250              ENDDO   ! enddo mmax
251     
252     ! w-momentum equations
253           CASE (Gas_W_Mom)
254              M=0
255              IF(.NOT.MOMENTUM_Z_EQ(M)) RETURN
256              DO IJK=IJKSTART3,IJKEND3
257                 IF (.NOT.FLUID_AT(IJK)) CYCLE
258                 K = K_OF(IJK)
259                 IJKT = TOP_OF(IJK)
260                 EPGA = AVG_Z(EP_G(IJK),EP_G(IJKT),K)
261                 IF (IP_AT_T(IJK) .OR. EPGA <= DIL_EP_S) CYCLE
262                 CALL USR_SOURCES(lEQ_NO, IJK, sourcelhs, sourcerhs, M, N)
263                 A_M(IJK,0,M) = A_M(IJK,0,M) - sourcelhs
264                 B_M(IJK,M) = B_M(IJK,M) - sourcerhs
265              ENDDO
266     
267     ! w-momentum equations
268           CASE (Solids_W_Mom)
269              DO L=1,MMAX
270                 IF (KT_TYPE_ENUM /= GHD_2007 .OR. &
271                     (KT_TYPE_ENUM == GHD_2007 .AND. L==MMAX)) THEN
272                    IF(.NOT.MOMENTUM_Z_EQ(L)) RETURN
273                    DO IJK=IJKSTART3,IJKEND3
274                       IF (.NOT.FLUID_AT(IJK)) CYCLE
275                       K = K_OF(IJK)
276                       IJKT = TOP_OF(IJK)
277                       EPtmp = ZERO
278                       IF (KT_TYPE_ENUM == GHD_2007) THEN
279                          DO lL = 1, SMAX
280                             EPtmp = EPtmp + AVG_Z(EP_S(IJK,lL),EP_S(IJKT,lL),K)
281                          ENDDO
282                          EPSA = EPtmp
283                       ELSE
284                          EPSA = AVG_Z(EP_S(IJK,L),EP_S(IJKT,L),K)
285                       ENDIF
286                       IF (IP_AT_T(IJK) .OR. SIP_AT_T(IJK) .OR. &
287                           EPSA <= DIL_EP_S) CYCLE
288                       CALL USR_SOURCES(lEQ_NO, IJK, sourcelhs, sourcerhs, L, N)
289                       A_M(IJK,0,L) = A_M(IJK,0,L) - sourcelhs
290                       B_M(IJK,L) = B_M(IJK,L) - sourcerhs
291                    ENDDO   ! enddo ijk
292                 ENDIF
293              ENDDO   ! enddo mmax
294     
295     
296     ! gas and solids energy equations (6)
297     ! gas and solids species mass fractions (7)
298     ! scalars, k_epsilon k and e (9)
299           CASE (Gas_Energy, Solids_Energy, Gas_Species, Solids_Species,&
300                 Usr_Scalar, K_Epsilon_K, K_Epsilon_E )
301               DO IJK = IJKSTART3, IJKEND3
302                 IF (FLUID_AT(IJK)) THEN
303                    IF (M==0) THEN
304                       eptmp=ep_g(ijk)
305                    ELSE
306                       eptmp=ep_s(ijk,M)
307                    ENDIF
308                    IF (eptmp <= DIL_EP_S) CYCLE
309                    CALL USR_SOURCES(lEQ_NO, IJK, sourcelhs, sourcerhs, M, N)
310                    A_M(IJK,0,M) = A_M(IJK,0,M) - sourcelhs
311                    B_M(IJK,M) = B_M(IJK,M) - sourcerhs
312                 ENDIF
313              ENDDO
314     
315     
316           CASE (Gran_Energy)   ! unique due to ghd
317     ! granular temperature
318              DO IJK = IJKSTART3, IJKEND3
319                 IF (FLUID_AT(IJK)) THEN
320                    IF (KT_TYPE_ENUM == GHD_2007) THEN
321                       eptmp = zero
322                       DO lL=1,SMAX
323                          eptmp=eptmp+EP_S(IJK,lL)
324                       ENDDO
325                    ELSE
326                       eptmp = ep_s(IJK,M)
327                    ENDIF
328                 ENDIF
329                 IF (eptmp<=dil_ep_s) cycle
330                 CALL USR_SOURCES(lEQ_NO, IJK, sourcelhs, sourcerhs, M, N)
331                 A_M(IJK,0,M) = A_M(IJK,0,M) - sourcelhs
332                 B_M(IJK,M) = B_M(IJK,M) - sourcerhs
333              ENDDO
334     
335     
336     ! error out
337           CASE DEFAULT
338     ! should never hit this
339     ! Initialize the error manager.
340              CALL INIT_ERR_MSG("CALC_USR_SOURCE")
341              WRITE(ERR_MSG, 1001) ival(leq_no)
342              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
343      1001 FORMAT('Error 1101: Unknown Equation= ', A)
344              CALL FINL_ERR_MSG
345           END SELECT   ! end selection of usr_source equation
346     
347           RETURN
348           END SUBROUTINE CALC_USR_SOURCE
349     
350           end module usr_src
351     
352