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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SOLVE_Pp_g                                              C
4     !  Purpose: Solve fluid pressure correction equation                   C
5     !                                                                      C
6     !  Author: M. Syamlal                                 Date: 19-JUN-96  C
7     !                                                                      C
8     !                                                                      C
9     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
10           SUBROUTINE SOLVE_PP_G(NORMG, RESG, IER)
11     
12     ! Modules
13     !---------------------------------------------------------------------//
14           use ambm, only: a_m, b_m, lock_ambm, unlock_ambm
15           use geometry, only: ijkmax2
16           use leqsol, only: leq_it, leq_method, leq_sweep, leq_pc, leq_tol
17           use param, only: dimension_3, dimension_m
18           use param1, only: zero, one, undefined
19           use pgcor, only: pp_g
20           use physprop, only: mmax, ro_g0
21           use ps, only: point_source
22           use residual, only: resid, max_resid, ijk_resid
23           use residual, only: num_resid, den_resid
24           use residual, only: resid_p
25           use run, only: momentum_x_eq, momentum_y_eq
26           use usr_src, only: call_usr_source, calc_usr_source
27           use usr_src, only: pressure_correction
28           IMPLICIT NONE
29     
30     ! Local parameters
31     !---------------------------------------------------------------------//
32     ! Parameter to make tolerance for residual scaled with max value
33     ! compatible with residual scaled with first iteration residual.
34     ! Increase it to tighten convergence.
35           DOUBLE PRECISION, PARAMETER :: DEN = 1.0D1   !5.0D2
36     
37     ! Dummy arguments
38     !---------------------------------------------------------------------//
39     ! Normalization factor for gas pressure correction residual.
40     ! At start of the iterate loop normg will either be 1 (i.e. not
41     ! normalized) or a user defined value given by norm_g.  If norm_g
42     ! was set to zero then the normalization is based on dominate
43     ! term in the equation
44           DOUBLE PRECISION, INTENT(IN) :: NORMg
45     ! gas pressure correction residual
46           DOUBLE PRECISION, INTENT(OUT) :: RESg
47     ! Error index
48           INTEGER, INTENT(INOUT) :: IER
49     
50     ! Local variables
51     !---------------------------------------------------------------------//
52     ! phase index
53           INTEGER :: M
54     ! Normalization factor for gas pressure correction residual
55           DOUBLE PRECISION :: NORMGloc
56     ! linear equation solver method and iterations
57           INTEGER :: LEQM, LEQI
58     
59     ! temporary use of global arrays:
60     ! arraym1 (locally b_mmax)
61     ! vector B_M based on dominate term in correction equation
62           DOUBLE PRECISION :: B_MMAX(DIMENSION_3, DIMENSION_M)
63     ! Septadiagonal matrix A_m, vector B_m
64     !      DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
65     !      DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
66     !---------------------------------------------------------------------//
67     
68           call lock_ambm
69     
70     ! initializing
71           PP_G(:) = ZERO
72           DO M = 0, MMAX
73              CALL INIT_AB_M (A_M, B_M, IJKMAX2, M)
74           ENDDO
75     
76     ! If gas momentum equations in x and y directions are not solved return
77           IF (.NOT.(MOMENTUM_X_EQ(0) .OR. MOMENTUM_Y_EQ(0)) .AND.&
78               RO_G0 .NE. UNDEFINED) THEN
79             call unlock_ambm
80             RETURN
81           ENDIF
82     
83     ! Forming the sparse matrix equation.
84           CALL CONV_PP_G (A_M, B_M)
85           CALL SOURCE_PP_G (A_M, B_M, B_MMAX)
86           IF(POINT_SOURCE) CALL POINT_SOURCE_PP_G (B_M, B_MMAX)
87           IF(CALL_USR_SOURCE(1)) CALL CALC_USR_SOURCE(Pressure_correction,&
88                                A_M, B_M, lB_MMAX=B_MMAX, lM=0)
89     
90     !      call check_ab_m(a_m, b_m, 0, .false., ier)
91     !      call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
92     
93     
94     ! Find average residual, maximum residual and location
95           NORMGloc = NORMG
96           IF(NORMG == ZERO) THEN
97     ! calculating the residual based on dominate term in correction equation
98     ! and use this to form normalization factor
99             CALL CALC_RESID_PP (B_MMAX, ONE, NUM_RESID(RESID_P,0), &
100              DEN_RESID(RESID_P,0), RESID(RESID_P,0), MAX_RESID(RESID_P,0), &
101              IJK_RESID(RESID_P,0))
102              NORMGloc = RESID(RESID_P,0)/DEN
103           ENDIF
104           CALL CALC_RESID_PP (B_M, NORMGloc, NUM_RESID(RESID_P,0),  &
105              DEN_RESID(RESID_P,0), RESID(RESID_P,0), MAX_RESID(RESID_P,0), &
106              IJK_RESID(RESID_P,0))
107           RESG = RESID(RESID_P,0)
108     !      write(*,*) resid(resid_p, 0), max_resid(resid_p, 0), &
109     !         ijk_resid(resid_p, 0)
110     
111     
112     ! Solve P_g_prime equation
113            LEQI = LEQ_IT(1)
114            LEQM = LEQ_METHOD(1)
115     !      CALL ADJUST_LEQ(RESID(RESID_P,0),LEQ_IT(1),LEQ_METHOD(1),LEQI,LEQM,IER)
116     
117     !     call check_symmetry(A_m, 0, IER)
118     !     call test_lin_eq(A_M, LEQ_IT(1),LEQ_METHOD(1), LEQ_SWEEP(1), LEQ_TOL(1), LEQ_PC(1),0,IER)
119           CALL SOLVE_LIN_EQ ('Pp_g', 1, PP_G, A_M, B_M, 0, LEQI, LEQM, &
120                              LEQ_SWEEP(1), LEQ_TOL(1), LEQ_PC(1), IER)
121     
122     !      call out_array(Pp_g, 'Pp_g')
123     
124           call unlock_ambm
125     
126           RETURN
127           END SUBROUTINE SOLVE_PP_G
128     
129     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
130     !                                                                      C
131     !  Subroutine: POINT_SOURCE_Pp_g                                       C
132     !  Purpose: Adds point sources to the Pressure correction equation.    C
133     !                                                                      C
134     !  Notes: The off-diagonal coefficients are positive. The center       C
135     !         coefficient and the source vector are negative. See          C
136     !         conv_Pp_g                                                    C
137     !                                                                      C
138     !  Author: J. Musser                                  Date: 10-JUN-13  C
139     !  Reviewer:                                          Date:            C
140     !                                                                      C
141     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
142           SUBROUTINE POINT_SOURCE_PP_G(B_M, B_mmax)
143     
144     ! Modules 
145     !-----------------------------------------------
146           use compar, only: dead_cell_at
147           use geometry, only: vol
148           use functions, only: fluid_at, funijk
149           use functions, only: is_on_myPe_plus2layers
150           use param, only: dimension_3, dimension_m
151           use param1, only: small_number
152           use ps, only: ps_defined, dimension_ps
153           use ps, only: ps_massflow_g, ps_volume
154           use ps, only: ps_k_b, ps_k_t
155           use ps, only: ps_j_s, ps_j_n
156           use ps, only: ps_i_w, ps_i_e
157           IMPLICIT NONE
158     
159     ! Dummy arguments
160     !-----------------------------------------------
161     ! Vector b_m
162           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
163     ! maximum term in b_m expression
164           DOUBLE PRECISION, INTENT(INOUT) :: B_mmax(DIMENSION_3, 0:DIMENSION_M)
165     
166     ! Local Variables
167     !-----------------------------------------------
168     ! Indices
169           INTEGER :: IJK, I, J, K
170           INTEGER :: PSV
171     
172     ! terms of bm expression
173           DOUBLE PRECISION pSource
174     
175     !-----------------------------------------------
176           PS_LP: do PSV = 1, DIMENSION_PS
177     
178              if(.NOT.PS_DEFINED(PSV)) cycle PS_LP
179              if(PS_MASSFLOW_G(PSV) < small_number) cycle PS_LP
180     
181              do k = PS_K_B(PSV), PS_K_T(PSV)
182              do j = PS_J_S(PSV), PS_J_N(PSV)
183              do i = PS_I_W(PSV), PS_I_E(PSV)
184     
185                 if(.NOT.IS_ON_myPE_plus2layers(I,J,K)) cycle
186                 IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
187     
188                 ijk = funijk(i,j,k)
189                 if(fluid_at(ijk)) then
190                    pSource = PS_MASSFLOW_G(PSV) * (VOL(IJK)/PS_VOLUME(PSV))
191     
192                    B_M(IJK,0) = B_M(IJK,0) - pSource
193                    B_MMAX(IJK,0) = max(abs(B_MMAX(IJK,0)), abs(B_M(IJK,0)))
194                 endif
195     
196              enddo
197              enddo
198              enddo
199     
200           enddo PS_LP
201     
202           RETURN
203           END SUBROUTINE POINT_SOURCE_PP_G
204