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