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) 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) 95 CALL SOURCE_PP_G (A_M, B_M, B_MMAX) 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)) 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)) 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 param1, only: small_number 159 use physprop 160 use ps 161 use run 162 use functions 163 164 IMPLICIT NONE 165 !----------------------------------------------- 166 ! Dummy arguments 167 !----------------------------------------------- 168 ! Vector b_m 169 DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M) 170 ! maximum term in b_m expression 171 DOUBLE PRECISION, INTENT(INOUT) :: B_mmax(DIMENSION_3, 0:DIMENSION_M) 172 !----------------------------------------------- 173 ! Local Variables 174 !----------------------------------------------- 175 176 ! Indices 177 INTEGER :: IJK, I, J, K 178 INTEGER :: PSV 179 180 ! terms of bm expression 181 DOUBLE PRECISION pSource 182 183 !----------------------------------------------- 184 PS_LP: do PSV = 1, DIMENSION_PS 185 186 if(.NOT.PS_DEFINED(PSV)) cycle PS_LP 187 if(PS_MASSFLOW_G(PSV) < small_number) cycle PS_LP 188 189 do k = PS_K_B(PSV), PS_K_T(PSV) 190 do j = PS_J_S(PSV), PS_J_N(PSV) 191 do i = PS_I_W(PSV), PS_I_E(PSV) 192 193 if(.NOT.IS_ON_myPE_plus2layers(I,J,K)) cycle 194 IF (DEAD_CELL_AT(I,J,K)) CYCLE ! skip dead cells 195 196 ijk = funijk(i,j,k) 197 if(fluid_at(ijk)) then 198 pSource = PS_MASSFLOW_G(PSV) * (VOL(IJK)/PS_VOLUME(PSV)) 199 200 B_M(IJK,0) = B_M(IJK,0) - pSource 201 B_MMAX(IJK,0) = max(abs(B_MMAX(IJK,0)), abs(B_M(IJK,0))) 202 endif 203 204 enddo 205 enddo 206 enddo 207 208 enddo PS_LP 209 210 RETURN 211 END SUBROUTINE POINT_SOURCE_PP_G 212