MFIX  2016-1
partial_elim.f
Go to the documentation of this file.
1 ! TODO:
2 ! 1. Account for solids-solids transfer terms for scalars. Will need
3 ! to pass the exchange coefficients to this routine (unlike momentum
4 ! eqs)
5 ! 2. repetitive code may be eliminated by defining and using functions
6 
7 
8 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
9 ! C
10 ! Subroutine: PARTIAL_ELIM_S C
11 ! Purpose: Do partial elimination for scalar quantities C
12 ! C
13 ! C
14 ! Author: M. Syamlal Date: 21-MAY-96 C
15 ! Reviewer: Date: C
16 ! C
17 ! Comments: C
18 ! Currently does not evaluate partial elimination terms between C
19 ! unlike solids phases. The exchange coefficient must be hard C
20 ! coded explicitly for this case - it is not a passed dummy C
21 ! argument C
22 ! C
23 ! Literature/Document References: C
24 ! Spalding, D.B., 1980, Numerical Computation of Multi-phase C
25 ! fluid flow and heat transfer, Recent Advances in Numerical C
26 ! Methods in Fluids, C. Taylor et al., eds, Pineridge Press C
27 ! Syamlal, M., 1998, MFIX Documentation: Numerical Technique, C
28 ! Technical Note, DOE/MC-31346-5824, NTIS/DE98002029, EG&G C
29 ! Technical Services of West Virginia, Inc., Morgantown, WV. C
30 ! C
31 ! C
32 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
33 
34  SUBROUTINE partial_elim_s(VAR_G, VAR_S, VXF, A_M, B_M)
35 
36 !-----------------------------------------------
37 ! Modules
38 !-----------------------------------------------
39  USE param
40  USE param1
41  USE parallel
42  USE geometry
43  USE physprop
44  USE indices
45  USE compar
46  USE drag
47  USE fldvar
48  USE functions
49  IMPLICIT NONE
50 !-----------------------------------------------
51 ! Dummy arguments
52 !-----------------------------------------------
53 ! gas phase variable
54  DOUBLE PRECISION, INTENT(IN) :: Var_g(dimension_3)
55 ! solids phase variable
56  DOUBLE PRECISION, INTENT(IN) :: Var_s(dimension_3, dimension_m)
57 ! Volume x gas-solids transfer coefficient
58  DOUBLE PRECISION, INTENT(IN) :: VxF(dimension_3, dimension_m)
59 ! Septadiagonal matrix A_m
60  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
61 ! Vector b_m
62  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
63 !-----------------------------------------------
64 ! Local variables
65 !-----------------------------------------------
66 ! Indices
67  INTEGER :: IJK, IJKW, IJKS, IJKB, IJKE, IJKN, IJKT
68  INTEGER :: L, M, LP
69 !
70  DOUBLE PRECISION :: SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME, DEN
71 ! a0, b0 etc.
72  DOUBLE PRECISION :: a(0:dimension_m), BB(0:dimension_m),&
73  F(0:dimension_m,0:dimension_m),&
74  Saxf(0:dimension_m)
75 !-----------------------------------------------
76 
77 !!$omp parallel do private( IJKW, IJKS, IJKB, IJKE, IJKN, IJKT, &
78 !!$omp& a, bb,F, Saxf,SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME,L, M,LP, den) &
79 !!$omp& schedule(static)
80  DO ijk = ijkstart3, ijkend3
81  IF (fluid_at(ijk)) THEN
82  ijkw = west_of(ijk)
83  ijks = south_of(ijk)
84  ijke = east_of(ijk)
85  ijkn = north_of(ijk)
86  DO m=0, mmax
87  a(m)=a_m(ijk,0,m)
88  bb(m)=b_m(ijk,m)
89 
90  if (m .ne. 0) then
91  f(m,0)=-vxf(ijk,m)
92  f(0,m)=f(m,0)
93  else
94  f(0,0) = zero
95  endif
96 
97  DO l =1, mmax
98  IF ((l .NE. m) .AND. (m .NE. 0)) THEN
99  f(m,l) = zero !insert solids-solids exchange coefficients here
100  ELSE
101  f(m,l) = zero
102  ENDIF
103  f(l,m) = zero
104  ENDDO
105 
106  IF (m == 0 ) THEN
107  saxf(m) = -(a_m(ijk,east,m)*var_g(ijke)+a_m(ijk,west,m)*var_g(ijkw&
108  )+a_m(ijk,north,m)*var_g(ijkn)+a_m(ijk,south,m)*var_g(ijks))
109  ELSE
110  saxf(m) = -(a_m(ijk,east,m)*var_s(ijke,m)+a_m(ijk,west,m)*var_s(&
111  ijkw,m)+a_m(ijk,north,m)*var_s(ijkn,m)+a_m(ijk,south,m)*var_s(&
112  ijks,m))
113  ENDIF
114 
115  IF (do_k) THEN
116  ijkb = bottom_of(ijk)
117  ijkt = top_of(ijk)
118  IF ( m ==0) THEN
119  saxf(m) = saxf(m) - (a_m(ijk,top,m)*var_g(ijkt)+a_m(ijk,bottom,m)*&
120  var_g(ijkb))
121  ELSE
122  saxf(m) = saxf(m) - (a_m(ijk,top,m)*var_s(ijkt,m)+a_m(ijk,bottom,m)*&
123  var_s(ijkb,m))
124  ENDIF
125  ENDIF
126  ENDDO
127 
128  DO m=0,mmax
129  sum_a = zero
130  sum_b = zero
131 
132  DO l =0,mmax
133  sum_a_lprime = zero
134  sum_b_lprime = zero
135 
136  DO lp=0,mmax
137  IF ( lp .NE. m) THEN
138  sum_a_lprime=sum_a_lprime+f(l,lp)
139  IF (lp == 0) THEN
140  sum_b_lprime=sum_b_lprime+f(l,lp)*var_g(ijk)
141  ELSE
142  sum_b_lprime=sum_b_lprime+f(l,lp)*var_s(ijk,lp)
143  ENDIF
144  ENDIF
145  ENDDO
146 
147  den = a(l) + sum_a_lprime + f(l,m)
148  IF ( den .NE. zero) THEN
149  sum_a = sum_a + ((f(l,m)*(a(l)+sum_a_lprime))/den)
150  sum_b = sum_b + f(l,m)*(saxf(l) + bb(l)+sum_b_lprime &
151  )/den
152  ENDIF
153  ENDDO
154 
155  a_m(ijk,0,m)= sum_a+a(m)
156  b_m(ijk,m) = sum_b+bb(m)
157  ENDDO
158  ENDIF
159  ENDDO
160 
161  RETURN
162  END SUBROUTINE partial_elim_s
163 
164 
165 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
166 ! C
167 ! Subroutine: PARTIAL_ELIM_IA C
168 ! Purpose: Do partial elimination for granular temperature coupling C
169 ! C
170 ! C
171 ! Author: J. Galvin Date: 21-MAY-05 C
172 ! Reviewer: S. Benyahia Date: C
173 ! C
174 ! Literature/Document References: C
175 ! Iddir, Y.H., Modeling of the multiphase mixture of particles C
176 ! using the kinetic theory approach, PhD Thesis, Illinois C
177 ! Institute of Technology, Chicago, Illinois, 2004 C
178 ! Iddir, Y.H., & H. Arastoopour, Modeling of multitype particle C
179 ! flow using the kinetic theory approach, AIChE J., Vol 51, C
180 ! No 6, June 2005 C
181 ! C
182 ! C
183 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
184 
185  SUBROUTINE partial_elim_ia(VAR_S, VXTCSS, A_M, B_M)
186 
187 !-----------------------------------------------
188 ! Modules
189 !-----------------------------------------------
190  USE param
191  USE param1
192  USE parallel
193  USE geometry
194  USE physprop
195  USE indices
196  USE compar
197  USE drag
198  USE fldvar
199  USE functions
200  IMPLICIT NONE
201 !-----------------------------------------------
202 ! Dummy arguments
203 !-----------------------------------------------
204 ! solids phase variable
205  DOUBLE PRECISION, INTENT(IN) :: Var_s(dimension_3, dimension_m)
206 ! Volume x solids-solids transfer coefficient
207  DOUBLE PRECISION, INTENT(IN) :: VxTCss (dimension_3, dimension_lm)
208 ! Septadiagonal matrix A_m
209  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
210 ! Vector b_m
211  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
212 !-----------------------------------------------
213 ! Local variables
214 !-----------------------------------------------
215 ! Indices
216  INTEGER :: IJK, IJKW, IJKS, IJKB, IJKE, IJKN, IJKT
217  INTEGER :: L, M, LP, LM
218 !
219  DOUBLE PRECISION :: SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME, DEN
220 ! a0, b0 etc.
221  DOUBLE PRECISION :: a(dimension_m), BB(dimension_m),&
222  F(dimension_m,dimension_m),&
223  Saxf(dimension_m)
224 !-----------------------------------------------
225 
226 !!$omp parallel do private( IJKW, IJKS, IJKB, IJKE, IJKN, IJKT, &
227 !!$omp& a, bb,F, Saxf,SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME,L, M,LP, den) &
228 !!$omp& schedule(static)
229 
230  DO ijk = ijkstart3, ijkend3
231 
232  IF (fluid_at(ijk)) THEN
233  ijkw = west_of(ijk)
234  ijks = south_of(ijk)
235  ijke = east_of(ijk)
236  ijkn = north_of(ijk)
237 
238  f(:,:) = zero
239 
240  DO m = 1, mmax
241  a(m)=a_m(ijk,0,m)
242  bb(m)=b_m(ijk,m)
243 
244  DO l = 1, mmax
245  IF ( l .NE. m ) THEN
246  lm = funlm(l,m)
247  f(m,l)=-vxtcss(ijk,lm)
248  ENDIF
249  f(l,m) = f(m,l)
250  ENDDO
251 
252  saxf(m) = -(a_m(ijk,east,m)*var_s(ijke,m)+&
253  a_m(ijk,west,m)*var_s(ijkw,m)+&
254  a_m(ijk,north,m)*var_s(ijkn,m)+&
255  a_m(ijk,south,m)*var_s(ijks,m))
256 
257  IF (do_k) THEN
258  ijkb = bottom_of(ijk)
259  ijkt = top_of(ijk)
260  saxf(m) = saxf(m) - (a_m(ijk,top,m)*var_s(ijkt,m)+&
261  a_m(ijk,bottom,m)*var_s(ijkb,m))
262  ENDIF
263  ENDDO
264 
265  DO m = 1, mmax
266  sum_a = zero
267  sum_b = zero
268 
269  DO l = 1, mmax
270  sum_a_lprime = zero
271  sum_b_lprime = zero
272 
273  DO lp = 1, mmax
274  IF ( lp .NE. m) THEN
275  sum_a_lprime=sum_a_lprime+f(l,lp)
276  sum_b_lprime=sum_b_lprime+f(l,lp)*var_s(ijk,lp)
277  ENDIF
278  ENDDO
279 
280  den = a(l) + sum_a_lprime + f(l,m)
281 
282  IF ( den .NE. zero) THEN
283  sum_a = sum_a + ((f(l,m)*(a(l)+sum_a_lprime))/den)
284  sum_b = sum_b + f(l,m)*(saxf(l) + bb(l)+&
285  sum_b_lprime)/den
286  ENDIF
287  ENDDO
288 
289  a_m(ijk,0,m)= sum_a+a(m)
290  b_m(ijk,m) = sum_b+bb(m)
291  ENDDO
292 
293  ENDIF
294  ENDDO
295 
296  RETURN
297  END SUBROUTINE partial_elim_ia
298 
299 
300 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
301 ! C
302 ! Subroutine: PARTIAL_ELIM_U C
303 ! Purpose: Do partial elimination for X vector quantities C
304 ! C
305 ! C
306 ! Author: M. Syamlal Date: 21-MAY-96 C
307 ! Reviewer: Date: C
308 ! C
309 ! Comments: C
310 ! Unlike gas-solids exchange coefficient (multiplied by volume), C
311 ! the solids-solids exchange coefficient is not as a passed dummy C
312 ! argument to this routine but simply re-evaluated here. C
313 ! C
314 ! Literature/Document References: C
315 ! Spalding, D.B., 1980, Numerical Computation of Multi-phase C
316 ! fluid flow and heat transfer, Recent Advances in Numerical C
317 ! Methods in Fluids, C. Taylor et al., eds, Pineridge Press C
318 ! Syamlal, M., 1998, MFIX Documentation: Numerical Technique, C
319 ! Technical Note, DOE/MC-31346-5824, NTIS/DE98002029, EG&G C
320 ! Technical Services of West Virginia, Inc., Morgantown, WV. C
321 ! C
322 ! C
323 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
324 
325  SUBROUTINE partial_elim_u(VAR_G, VAR_S, VXF, A_M, B_M)
326 
327 !-----------------------------------------------
328 ! Modules
329 !-----------------------------------------------
330  USE param
331  USE param1
332  USE parallel
333  USE geometry
334  USE physprop
335  USE indices
336  USE run
337  USE compar
338  USE drag
339  USE fldvar
340  USE fun_avg
341  USE functions
342  IMPLICIT NONE
343 !-----------------------------------------------
344 ! Dummy arguments
345 !-----------------------------------------------
346 ! gas phase variable
347  DOUBLE PRECISION, INTENT(IN) :: Var_g(dimension_3)
348 ! solids phase variable
349  DOUBLE PRECISION, INTENT(IN) :: Var_s(dimension_3, dimension_m)
350 ! Volume x gas-solids transfer coefficient
351  DOUBLE PRECISION, INTENT(IN) :: VxF(dimension_3, dimension_m)
352 ! Septadiagonal matrix A_m
353  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
354 ! Vector b_m
355  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
356 !-----------------------------------------------
357 ! Local variables
358 !-----------------------------------------------
359 ! Indices
360  INTEGER :: IJK, IMJK, IJMK, IJKM, IPJK, IJPK, IJKP, I, IJKE
361  INTEGER :: L, M, LP, LM
362 !
363  DOUBLE PRECISION :: SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME, DEN
364 ! a0, b0 etc.
365  DOUBLE PRECISION :: a(0:dimension_m), BB(0:dimension_m), &
366  F(0:dimension_m,0:dimension_m),&
367  Saxf(0:dimension_m)
368 !-----------------------------------------------
369 
370 !!$omp parallel do private( IMJK, IJMK, IJKM, IPJK, IJPK, IJKP, &
371 !!$omp& a, bb,F, Saxf,SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME,L, M,LP, den) &
372 !!$omp& schedule(static)
373  DO ijk = ijkstart3, ijkend3
374  IF (flow_at_e(ijk)) THEN
375  imjk = im_of(ijk)
376  ijmk = jm_of(ijk)
377  ipjk = ip_of(ijk)
378  ijpk = jp_of(ijk)
379 
380  f = zero
381  DO m=0, mmax
382  a(m) = a_m(ijk,0,m)
383  bb(m)= b_m(ijk,m)
384 
385 ! locally storing vol x avg gas-solids drag coefficient at momentum cell
386 ! center
387  IF (m .NE. 0) THEN
388  f(m,0)=-vxf(ijk,m)
389  f(0,m) = f(m,0)
390  ENDIF
391 
392 ! locally storing vol x avg solids-solids drag coefficient at momentum cell
393 ! center
394  DO l =1, mmax
395  IF ((l .NE. m) .AND. (m .NE. 0)) THEN
396  lm = funlm(l,m)
397  IF (.NOT.ip_at_e(ijk)) THEN
398  i = i_of(ijk)
399  ijke = east_of(ijk)
400  f(m,l) = -avg_x(f_ss(ijk,lm),f_ss(ijke,lm),i)*vol_u(ijk)
401  ENDIF
402  ENDIF
403  f(l,m)=f(m,l)
404  ENDDO
405 
406  IF (m == 0) THEN
407  saxf(m) = -(a_m(ijk,east,m)*var_g(ipjk)+a_m(ijk,west,m)*var_g(imjk&
408  )+a_m(ijk,north,m)*var_g(ijpk)+a_m(ijk,south,m)*var_g(ijmk))
409  ELSE
410  saxf(m) = -(a_m(ijk,east,m)*var_s(ipjk,m)+a_m(ijk,west,m)*var_s(&
411  imjk,m)+a_m(ijk,north,m)*var_s(ijpk,m)+a_m(ijk,south,m)*var_s(&
412  ijmk,m))
413  ENDIF
414 
415  IF (do_k) THEN
416  ijkm = km_of(ijk)
417  ijkp = kp_of(ijk)
418  IF (m ==0) THEN
419  saxf(m) = saxf(m) - (a_m(ijk,top,m)*var_g(ijkp)+a_m(ijk,bottom,m)*&
420  var_g(ijkm))
421  ELSE
422  saxf(m) = saxf(m) - (a_m(ijk,top,m)*var_s(ijkp,m)+a_m(ijk,bottom,m)*&
423  var_s(ijkm,m))
424  ENDIF
425  ENDIF
426  ENDDO ! end do loop for M over 0 to MMAX
427 
428  DO m=0,mmax
429  IF (momentum_x_eq(m)) THEN
430  sum_a = zero
431  sum_b = zero
432 
433  DO l =0,mmax
434 ! - should only need to loop for L/=M. since the resulting term should
435 ! evaluate to zero anyway, don't bother to add an explicit conditional.
436 ! - if phase L's momentum equation is not solved then do not evaluate the
437 ! outer sum for phase L (as phase L has no momentum equation). however,
438 ! the inner sum may include terms from phase L as all phases still see
439 ! phase L and therefore experience drag with phase L
440  IF (momentum_x_eq(l)) THEN
441  sum_a_lprime = zero
442  sum_b_lprime = zero
443 
444  DO lp=0,mmax ! only need to loop for LP/=L; but
445  ! evaluates to zero so don't bother
446  ! with an explicit conditional
447  IF ( lp .NE. m) THEN
448  sum_a_lprime=sum_a_lprime+f(l,lp)
449  IF (lp == 0) THEN
450  sum_b_lprime=sum_b_lprime+f(l,lp)*var_g(ijk)
451  ELSE
452  sum_b_lprime=sum_b_lprime+f(l,lp)*var_s(ijk,lp)
453  ENDIF
454  ENDIF
455  ENDDO
456 
457  den = a(l) + sum_a_lprime + f(l,m)
458  IF ( den .NE. zero) THEN
459  sum_a = sum_a + ((f(l,m)*(a(l)+sum_a_lprime))/den)
460  sum_b = sum_b + f(l,m)*(saxf(l) + bb(l)+sum_b_lprime &
461  )/den
462  ENDIF
463  ELSE
464  sum_a = sum_a + f(l,m)
465  IF (l == 0) THEN
466  sum_b = sum_b + f(l,m)*var_g(ijk)
467  ELSE
468  sum_b = sum_b + f(l,m)*var_s(ijk,l)
469  ENDIF
470  ENDIF ! end if momentum_x_eq(L)
471  ENDDO !end do loop for L over 0 to mmax
472 
473  a_m(ijk,0,m) = sum_a+a(m)
474  b_m(ijk,m) = sum_b+bb(m)
475  ENDIF ! end if momentum_x_eq(M)
476  ENDDO ! end do loop for M over 0 to MMAX
477 
478  ENDIF ! end if flow_at_e(ijk)
479  ENDDO ! end do loop for ijk over ijkstart3,ijkend3
480 
481  RETURN
482  END SUBROUTINE partial_elim_u
483 
484 
485 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
486 ! C
487 ! Module name: PARTIAL_ELIM_V C
488 ! Purpose: Do partial elimination for Y vector quantities C
489 ! C
490 ! C
491 ! Author: M. Syamlal Date: 21-MAY-96 C
492 ! Reviewer: Date: C
493 ! C
494 ! Comments: C
495 ! Unlike gas-solids exchange coefficient (multiplied by volume), C
496 ! the solids-solids exchange coefficient is not as a passed dummy C
497 ! argument to this routine but simply re-evaluated here. C
498 ! C
499 ! Literature/Document References: C
500 ! Spalding, D.B., 1980, Numerical Computation of Multi-phase C
501 ! fluid flow and heat transfer, Recent Advances in Numerical C
502 ! Methods in Fluids, C. Taylor et al., eds, Pineridge Press C
503 ! Syamlal, M., 1998, MFIX Documentation: Numerical Technique, C
504 ! Technical Note, DOE/MC-31346-5824, NTIS/DE98002029, EG&G C
505 ! Technical Services of West Virginia, Inc., Morgantown, WV. C
506 ! C
507 ! C
508 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
509 
510  SUBROUTINE partial_elim_v(VAR_G, VAR_S, VXF, A_M, B_M)
511 
512 !-----------------------------------------------
513 ! Modules
514 !-----------------------------------------------
515  USE param
516  USE param1
517  USE parallel
518  USE geometry
519  USE physprop
520  USE indices
521  USE run
522  USE compar
523  USE drag
524  USE fldvar
525  USE fun_avg
526  USE functions
527  IMPLICIT NONE
528 !-----------------------------------------------
529 ! Dummy arguments
530 !-----------------------------------------------
531 ! gas phase variable
532  DOUBLE PRECISION, INTENT(IN) :: Var_g(dimension_3)
533 ! solids phase variable
534  DOUBLE PRECISION, INTENT(IN) :: Var_s(dimension_3, dimension_m)
535 ! Volume x gas-solids transfer coefficient
536  DOUBLE PRECISION, INTENT(IN) :: VxF(dimension_3, dimension_m)
537 ! Septadiagonal matrix A_m
538  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
539 ! Vector b_m
540  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
541 !-----------------------------------------------
542 ! Local variables
543 !-----------------------------------------------
544 ! Indices
545  INTEGER :: IJK, IMJK, IJMK, IJKM, IPJK, IJPK, IJKP, IJKN, J
546  INTEGER :: L, M, LP, LM
547 !
548  DOUBLE PRECISION :: SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME, DEN
549 ! a0, b0 etc.
550  DOUBLE PRECISION :: a(0:dimension_m), BB(0:dimension_m),&
551  F(0:dimension_m,0:dimension_m),&
552  Saxf(0:dimension_m)
553 !-----------------------------------------------
554 
555 !!$omp parallel do private( IMJK, IJMK, IJKM, IPJK, IJPK, IJKP, &
556 !!$omp& a, bb,F, Saxf,SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME,L, M,LP, DEN) &
557 !!$omp& schedule(static)
558  DO ijk = ijkstart3, ijkend3
559  IF (flow_at_n(ijk)) THEN
560  imjk = im_of(ijk)
561  ijmk = jm_of(ijk)
562  ipjk = ip_of(ijk)
563  ijpk = jp_of(ijk)
564 
565  f = zero
566  DO m = 0, mmax
567  a(m)=a_m(ijk,0,m)
568  bb(m)=b_m(ijk,m)
569 
570 ! locally storing vol x avg gas-solids drag coefficient at momentum cell
571 ! center
572  IF (m .NE. 0) THEN
573  f(m,0)=-vxf(ijk,m)
574  f(0,m)=f(m,0)
575  ENDIF
576 
577 ! locally storing vol x avg solids-solids drag coefficient at momentum cell
578 ! center
579  DO l =1, mmax
580  IF ((l .NE. m) .AND. (m .NE. 0)) THEN
581  lm = funlm(l,m)
582  IF (.NOT.ip_at_n(ijk)) THEN
583  j = j_of(ijk)
584  ijkn = north_of(ijk)
585  f(m,l)=-avg_y(f_ss(ijk,lm),f_ss(ijkn,lm),j)*vol_v(ijk)
586  ENDIF
587  ENDIF
588  f(l,m)=f(m,l)
589  ENDDO
590 
591  IF (m == 0) THEN
592  saxf(m) = -(a_m(ijk,east,m)*var_g(ipjk)+a_m(ijk,west,m)*var_g(imjk&
593  )+a_m(ijk,north,m)*var_g(ijpk)+a_m(ijk,south,m)*var_g(ijmk))
594  ELSE
595  saxf(m) = -(a_m(ijk,east,m)*var_s(ipjk,m)+a_m(ijk,west,m)*var_s(&
596  imjk,m)+a_m(ijk,north,m)*var_s(ijpk,m)+a_m(ijk,south,m)*var_s(&
597  ijmk,m))
598  ENDIF
599 
600  IF (do_k) THEN
601  ijkm = km_of(ijk)
602  ijkp = kp_of(ijk)
603  IF (m == 0) THEN
604  saxf(m) = saxf(m) - (a_m(ijk,top,m)*var_g(ijkp)+a_m(ijk,bottom,m)*&
605  var_g(ijkm))
606  ELSE
607  saxf(m) = saxf(m) - (a_m(ijk,top,m)*var_s(ijkp,m)+a_m(ijk,bottom,m)*&
608  var_s(ijkm,m))
609  ENDIF
610  ENDIF
611  ENDDO ! end do loop for M over 0 to MMAX
612 
613 
614  DO m=0,mmax
615  IF (momentum_y_eq(m)) THEN
616  sum_a = zero
617  sum_b = zero
618 
619  DO l =0,mmax
620  IF(momentum_y_eq(l)) THEN
621  sum_a_lprime = zero
622  sum_b_lprime = zero
623 
624  DO lp=0,mmax
625  IF ( lp .NE. m) THEN
626  sum_a_lprime=sum_a_lprime+f(l,lp)
627  IF (lp == 0) THEN
628  sum_b_lprime=sum_b_lprime+f(l,lp)*var_g(ijk)
629  ELSE
630  sum_b_lprime=sum_b_lprime+f(l,lp)*var_s(ijk,lp)
631  ENDIF
632  ENDIF
633  ENDDO
634 
635  den = a(l) + sum_a_lprime + f(l,m)
636  IF ( den .NE. zero) THEN
637  sum_a = sum_a + ((f(l,m)*(a(l)+sum_a_lprime))/den)
638  sum_b = sum_b + f(l,m)*(saxf(l) + bb(l)+sum_b_lprime &
639  )/den
640  ENDIF
641  ELSE
642  sum_a = sum_a + f(l,m)
643  IF (l == 0) THEN
644  sum_b = sum_b + f(l,m)*var_g(ijk)
645  ELSE
646  sum_b = sum_b + f(l,m)*var_s(ijk,l)
647  ENDIF
648  ENDIF ! end if momentum_y_eq(l)
649  ENDDO !end do loop for L over 0 to mmax
650 
651  a_m(ijk,0,m)=sum_a+a(m)
652  b_m(ijk,m) = sum_b+bb(m)
653  ENDIF ! end if momentum_y_eq(M)
654  ENDDO ! end do loop for M over 0 to MMAX
655 
656  ENDIF ! end if flow_at_n(ijk)
657  ENDDO ! end do loop for ijk over ijkstart3,ijkend3
658 
659  RETURN
660  END SUBROUTINE partial_elim_v
661 
662 
663 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
664 ! C
665 ! Subroutine: PARTIAL_ELIM_W C
666 ! Purpose: Do partial elimination for Z vector quantities
667 ! C
668 ! C
669 ! Author: M. Syamlal Date: 21-MAY-96 C
670 ! Reviewer: Date: C
671 ! C
672 ! Comments: C
673 ! Unlike gas-solids exchange coefficient (multiplied by volume), C
674 ! the solids-solids exchange coefficient is not as a passed dummy C
675 ! argument to this routine but simply re-evaluated here. C
676 ! C
677 ! Literature/Document References: C
678 ! Spalding, D.B., 1980, Numerical Computation of Multi-phase C
679 ! fluid flow and heat transfer, Recent Advances in Numerical C
680 ! Methods in Fluids, C. Taylor et al., eds, Pineridge Press C
681 ! Syamlal, M., 1998, MFIX Documentation: Numerical Technique, C
682 ! Technical Note, DOE/MC-31346-5824, NTIS/DE98002029, EG&G C
683 ! Technical Services of West Virginia, Inc., Morgantown, WV. C
684 ! C
685 ! C
686 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
687 
688  SUBROUTINE partial_elim_w(VAR_G, VAR_S, VXF, A_M, B_M)
689 
690 !-----------------------------------------------
691 ! Modules
692 !-----------------------------------------------
693  USE param
694  USE param1
695  USE parallel
696  USE geometry
697  USE physprop
698  USE indices
699  USE run
700  USE compar
701  USE drag
702  USE fldvar
703  USE fun_avg
704  USE functions
705  IMPLICIT NONE
706 !-----------------------------------------------
707 ! Dummy arguments
708 !-----------------------------------------------
709 ! gas phase variable
710  DOUBLE PRECISION, INTENT(IN) :: Var_g(dimension_3)
711 ! solids phase variable
712  DOUBLE PRECISION, INTENT(IN) :: Var_s(dimension_3, dimension_m)
713 ! Volume x gas-solids transfer coefficient
714  DOUBLE PRECISION, INTENT(IN) :: VxF(dimension_3, dimension_m)
715 ! Septadiagonal matrix A_m
716  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
717 ! Vector b_m
718  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
719 !-----------------------------------------------
720 ! Local variables
721 !-----------------------------------------------
722 ! Indices
723  INTEGER :: IJK, IMJK, IJMK, IJKM, IPJK, IJPK, IJKP, K, IJKT
724  INTEGER :: L, M, LP, LM
725 !
726  DOUBLE PRECISION :: SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME, DEN
727 ! a0, b0 etc.
728  DOUBLE PRECISION :: a(0:dimension_m), BB(0:dimension_m), &
729  F(0:dimension_m,0:dimension_m),&
730  Saxf(0:dimension_m)
731 !-----------------------------------------------
732 
733 !!$omp parallel do private( IMJK, IJMK, IJKM, IPJK, IJPK, IJKP, &
734 !!$omp& a, bb, F, Saxf,SUM_A, SUM_B, SUM_A_LPRIME,SUM_B_LPRIME,L, M,LP, DEN) &
735 !!$omp& schedule(static)
736  DO ijk = ijkstart3, ijkend3
737  IF (flow_at_t(ijk)) THEN
738  imjk = im_of(ijk)
739  ijmk = jm_of(ijk)
740  ipjk = ip_of(ijk)
741  ijpk = jp_of(ijk)
742 
743  f = zero
744  DO m=0, mmax
745  a(m)=a_m(ijk,0,m)
746  bb(m)=b_m(ijk,m)
747 
748 ! locally storing vol x avg gas-solids drag coefficient at momentum cell
749 ! center
750  IF (m .NE. 0) THEN
751  f(m,0)=-vxf(ijk,m)
752  f(0,m)=f(m,0)
753  ENDIF
754 
755 ! locally storing vol x avg solids-solids drag coefficient at momentum cell
756 ! center
757  DO l =1, mmax
758  IF ((l .NE. m) .AND. (m .NE. 0)) THEN
759  lm = funlm(l,m)
760  IF (.NOT.ip_at_t(ijk)) THEN
761  k = k_of(ijk)
762  ijkt = top_of(ijk)
763  f(m,l) = -avg_z(f_ss(ijk,lm),f_ss(ijkt,lm),k)*vol_w(ijk)
764  ENDIF
765  ENDIF
766  f(l,m)=f(m,l)
767  ENDDO
768 
769  IF (m == 0) THEN
770  saxf(m) = -(a_m(ijk,east,m)*var_g(ipjk)+a_m(ijk,west,m)*var_g(imjk&
771  )+a_m(ijk,north,m)*var_g(ijpk)+a_m(ijk,south,m)*var_g(ijmk))
772  ELSE
773  saxf(m) = -(a_m(ijk,east,m)*var_s(ipjk,m)+a_m(ijk,west,m)*var_s(&
774  imjk,m)+a_m(ijk,north,m)*var_s(ijpk,m)+a_m(ijk,south,m)*var_s(&
775  ijmk,m))
776  ENDIF
777 
778  IF (do_k) THEN
779  ijkm = km_of(ijk)
780  ijkp = kp_of(ijk)
781  IF (m == 0) THEN
782  saxf(m) = saxf(m) - (a_m(ijk,top,m)*var_g(ijkp)+a_m(ijk,bottom,m)*&
783  var_g(ijkm))
784  ELSE
785  saxf(m) = saxf(m) - (a_m(ijk,top,m)*var_s(ijkp,m)+a_m(ijk,bottom,m)*&
786  var_s(ijkm,m))
787  ENDIF
788  ENDIF
789  ENDDO ! end do loop for M over 0 to MMAX
790 
791 
792  DO m=0,mmax
793  IF (momentum_z_eq(m)) THEN
794  sum_a = zero
795  sum_b = zero
796 
797  DO l =0,mmax
798  IF (momentum_z_eq(l)) THEN
799  sum_a_lprime = zero
800  sum_b_lprime = zero
801 
802  DO lp=0,mmax
803  IF ( lp .NE. m) THEN
804  sum_a_lprime=sum_a_lprime+f(l,lp)
805  IF (lp == 0) THEN
806  sum_b_lprime=sum_b_lprime+f(l,lp)*var_g(ijk)
807  ELSE
808  sum_b_lprime=sum_b_lprime+f(l,lp)*var_s(ijk,lp)
809  ENDIF
810  ENDIF
811  ENDDO
812 
813  den = a(l) + sum_a_lprime + f(l,m)
814  IF ( den .NE. zero) THEN
815  sum_a = sum_a + ((f(l,m)*(a(l)+sum_a_lprime))/den)
816  sum_b = sum_b + f(l,m)*(saxf(l) + bb(l)+sum_b_lprime &
817  )/den
818  ENDIF
819  ELSE
820  sum_a = sum_a + f(l,m)
821  IF (l == 0) THEN
822  sum_b = sum_b + f(l,m)*var_g(ijk)
823  ELSE
824  sum_b = sum_b + f(l,m)*var_s(ijk,l)
825  ENDIF
826  ENDIF ! end if momentum_z_eq(L)
827  ENDDO ! end do loop for L over 0 to mmax
828 
829  a_m(ijk,0,m)=sum_a+a(m)
830  b_m(ijk,m) = sum_b+bb(m)
831  ENDIF ! end if momentum_Z_eq(M)
832  ENDDO ! end do loop for M over 0 to MMAX
833 
834  ENDIF ! end if flow_at_t(ijk)
835  ENDDO ! end do loop for ijk over ijkstart3,ijkend3
836 
837  RETURN
838  END SUBROUTINE partial_elim_w
839 
840 
partial_elim_u
subroutine partial_elim_u(VAR_G, VAR_S, VXF, A_M, B_M)
Definition: partial_elim.f:326
run::momentum_y_eq
logical, dimension(0:dim_m) momentum_y_eq
Definition: run_mod.f:77
geometry::vol_w
double precision, dimension(:), allocatable vol_w
Definition: geometry_mod.f:242
indices::i_of
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
param1::dimension_lm
integer dimension_lm
Definition: param1_mod.f:13
param::dimension_3
integer dimension_3
Definition: param_mod.f:11
drag::f_ss
double precision, dimension(:,:), allocatable f_ss
Definition: drag_mod.f:17
run::momentum_x_eq
logical, dimension(0:dim_m) momentum_x_eq
Definition: run_mod.f:74
drag
Definition: drag_mod.f:11
run::momentum_z_eq
logical, dimension(0:dim_m) momentum_z_eq
Definition: run_mod.f:80
param::east
integer east
Definition: param_mod.f:29
partial_elim_w
subroutine partial_elim_w(VAR_G, VAR_S, VXF, A_M, B_M)
Definition: partial_elim.f:689
indices::k_of
integer, dimension(:), allocatable k_of
Definition: indices_mod.f:47
fun_avg
Definition: fun_avg_mod.f:1
partial_elim_ia
subroutine partial_elim_ia(VAR_S, VXTCSS, A_M, B_M)
Definition: partial_elim.f:186
physprop::mmax
integer mmax
Definition: physprop_mod.f:19
indices::j_of
integer, dimension(:), allocatable j_of
Definition: indices_mod.f:46
param::north
integer north
Definition: param_mod.f:37
param::south
integer south
Definition: param_mod.f:41
run
Definition: run_mod.f:13
param
Definition: param_mod.f:2
geometry::do_k
logical do_k
Definition: geometry_mod.f:30
compar::ijkstart3
integer ijkstart3
Definition: compar_mod.f:80
param::west
integer west
Definition: param_mod.f:33
param::top
integer top
Definition: param_mod.f:45
geometry::vol_u
double precision, dimension(:), allocatable vol_u
Definition: geometry_mod.f:224
partial_elim_v
subroutine partial_elim_v(VAR_G, VAR_S, VXF, A_M, B_M)
Definition: partial_elim.f:511
param::dimension_m
integer dimension_m
Definition: param_mod.f:18
param::bottom
integer bottom
Definition: param_mod.f:49
partial_elim_s
subroutine partial_elim_s(VAR_G, VAR_S, VXF, A_M, B_M)
Definition: partial_elim.f:35
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
geometry::vol_v
double precision, dimension(:), allocatable vol_v
Definition: geometry_mod.f:233