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

1     ! -*- f90 -*-
2     MODULE cont
3     
4     ! Indicates whether the continuity equation needs to be
5     ! solved
6     
7           LOGICAL, DIMENSION(:), ALLOCATABLE ::  DO_CONT
8     
9     CONTAINS
10     
11     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
12     !                                                                      C
13     !  Subroutine: SOLVE_CONTINUITY                                        C
14     !  Purpose: Solve for solids bulk density (i.e., material density      C
15     !           multiplied by volume fraction).                            C
16     !                                                                      C
17     !  Author: M. Syamlal                                 Date: 2-JUL-96   C
18     !                                                                      C
19     !                                                                      C
20     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
21           SUBROUTINE SOLVE_CONTINUITY(M,IER)
22     
23     ! Modules
24     !---------------------------------------------------------------------//
25           use ambm, only: a_m, b_m, lock_ambm, unlock_ambm
26           use fldvar, only: rop_g, rop_s
27           use geometry, only: ijkmax2
28           use leqsol, only: leq_it, leq_method, leq_sweep, leq_pc, leq_tol
29           use param, only: dimension_3, dimension_m
30           use param1, only: undefined_i, zero, one
31           use ps, only: point_source
32           use residual, only: resid, max_resid, ijk_resid
33           use residual, only: num_resid, den_resid
34           use residual, only: resid_ro
35           use usr_src, only: call_usr_source, calc_usr_source
36           use usr_src, only: gas_continuity, solids_continuity
37           IMPLICIT NONE
38     
39     ! Dummy arguments
40     !---------------------------------------------------------------------//
41     ! phase index
42           INTEGER, INTENT(IN) :: M
43     ! error index
44           INTEGER, INTENT(INOUT) :: IER
45     
46     ! Local variables
47     !---------------------------------------------------------------------//
48     ! solids volume fraction residual
49           DOUBLE PRECISION :: RESs
50     ! linear equation solver method and iterations
51           INTEGER :: LEQM, LEQI
52     
53     ! temporary use of global arrays:
54     ! Septadiagonal matrix A_m, vector B_m
55     !      DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
56     !      DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
57     !---------------------------------------------------------------------//
58     
59           call lock_ambm
60     
61           IF (M==0) THEN
62     ! solve gas phase continuity equation. note that this branch will never
63     ! be entered given that the calling subroutine (iterate) only calls
64     ! solve_continuity when M>1
65     
66     ! initializing
67              CALL INIT_AB_M (A_M, B_M, IJKMAX2, 0)
68     
69     ! forming the matrix equation
70              CALL CONV_ROP_G (A_M, B_M)
71              CALL SOURCE_ROP_G (A_M, B_M)
72              IF(CALL_USR_SOURCE(2)) CALL CALC_USR_SOURCE (GAS_CONTINUITY,&
73                                   A_M, B_M, lM=0)
74     
75     ! calculating the residual
76              CALL CALC_RESID_C (ROP_G, A_M, B_M, 0, NUM_RESID(RESID_RO,0), &
77                 DEN_RESID(RESID_RO,0), RESID(RESID_RO,0), MAX_RESID(&
78                 RESID_RO,0), IJK_RESID(RESID_RO,0))
79     
80     !         call check_ab_m(a_m, b_m, 0, .true., ier)
81     !         call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
82     !         write(*,*) resid(resid_ro, 0), max_resid(resid_ro, 0), &
83     !            ijk_resid(resid_ro, 0)
84     !         call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0),&
85     !            1, DO_K, ier)
86     
87     ! solving gas continuity
88              CALL ADJUST_LEQ (RESID(RESID_RO,0), LEQ_IT(2), LEQ_METHOD(2),&
89                 LEQI, LEQM)
90              CALL SOLVE_LIN_EQ ('ROP_g', 2, ROP_G, A_M, B_M, 0, LEQI, &
91                 LEQM, LEQ_SWEEP(2), LEQ_TOL(2), LEQ_PC(2), IER)
92              CALL ADJUST_ROP (ROP_G)
93     !        call out_array(ROP_g, 'rop_g')
94     
95           ELSE
96     ! solve solids phase M continuity equation.
97     
98     ! initializing
99              CALL INIT_AB_M (A_M, B_M, IJKMAX2, M)
100     
101     ! forming the matrix equation
102              CALL CONV_ROP_S (A_M, B_M, M)
103              CALL SOURCE_ROP_S (A_M, B_M, M)
104              IF(POINT_SOURCE) CALL POINT_SOURCE_ROP_S (B_M, M)
105              IF(CALL_USR_SOURCE(2)) CALL CALC_USR_SOURCE (SOLIDS_CONTINUITY, &
106                                   A_M, B_M, lM=M)
107     
108              CALL CALC_RESID_C (ROP_S(1,M), A_M, B_M, M, &
109                 NUM_RESID(RESID_RO,M), DEN_RESID(RESID_RO,M), &
110                 RESID(RESID_RO,M), MAX_RESID(RESID_RO,M), &
111                 IJK_RESID(RESID_RO,M))
112              RESS = RESID(RESID_RO,M)
113     
114     !         call check_ab_m(a_m, b_m, m, .true., ier)
115     !         write(*,*) 'solve_cont= ', resid(resid_ro, m),&
116     !                    max_resid(resid_ro, m), ijk_resid(resid_ro, m), m
117     !         call write_ab_m(a_m, b_m, ijkmax2, m, ier)
118     !         call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, M), 1, &
119     !                          DO_K, ier)
120     !
121              CALL ADJUST_LEQ (RESID(RESID_RO,M), LEQ_IT(2), LEQ_METHOD(2),&
122                 LEQI, LEQM)
123              CALL SOLVE_LIN_EQ ('ROP_s', 2, ROP_S(1,M), A_M, B_M, M, LEQI,&
124                 LEQM,LEQ_SWEEP(2), LEQ_TOL(2), LEQ_PC(2), IER)
125              CALL ADJUST_ROP (ROP_S(1,M))
126     !         call out_array(rop_s(1,m), 'rop_s')
127     
128           ENDIF
129     
130           call unlock_ambm
131           RETURN
132     
133         CONTAINS
134     
135     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
136     !                                                                      C
137     !  Subroutine: ADJUST_ROP                                              C
138     !  Purpose: Remove small negative values of density.                   C
139     !                                                                      C
140     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
141           SUBROUTINE ADJUST_ROP(ROP)
142     
143     ! Modules
144     !---------------------------------------------------------------------//
145           use compar, only: ijkstart3, ijkend3
146           USE functions, only: fluid_at
147           USE param1, only: zero
148           IMPLICIT NONE
149     
150     ! Dummy arguments
151     !---------------------------------------------------------------------//
152     ! density
153           DOUBLE PRECISION, INTENT(INOUT) :: ROP(DIMENSION_3)
154     
155     ! Local variables
156     !---------------------------------------------------------------------//
157     ! Indices
158           INTEGER :: IJK
159     !---------------------------------------------------------------------//
160     
161           DO IJK = ijkstart3, ijkend3
162              IF (FLUID_AT(IJK)) ROP(IJK) = DMAX1(ZERO,ROP(IJK))
163           ENDDO
164     
165           RETURN
166           END SUBROUTINE ADJUST_ROP
167     
168           END SUBROUTINE SOLVE_CONTINUITY
169     
170           END MODULE cont
171