1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 2 ! C 3 ! Subroutine: SOLVE_Epp C 4 ! Purpose: Solve solids volume fraction correction equation. C 5 ! C 6 ! Notes: MCP must be defined to call this routine. C 7 ! C 8 ! Author: M. Syamlal Date: 25-SEP-96 C 9 ! C 10 ! C 11 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 12 SUBROUTINE SOLVE_EPP(NORMS, RESS, IER) 13 14 ! Modules 15 !---------------------------------------------------------------------// 16 use ambm, only: a_m, b_m, lock_ambm, unlock_ambm 17 use geometry, only: ijkmax2 18 use leqsol, only: leq_it, leq_method, leq_sweep, leq_pc, leq_tol 19 use param, only: dimension_3, dimension_m 20 use param1, only: undefined_i, zero, one 21 use ps, only: point_source 22 use pscor, only: epp, mcp 23 use residual, only: resid, max_resid, ijk_resid 24 use residual, only: num_resid, den_resid 25 use residual, only: resid_p 26 use run, only: momentum_x_eq, momentum_y_eq 27 use usr_src, only: call_usr_source, calc_usr_source 28 use usr_src, only: solids_correction 29 IMPLICIT NONE 30 31 ! Local parameters 32 !---------------------------------------------------------------------// 33 ! Parameter to make tolerance for residual scaled with max value 34 ! compatible with residual scaled with first iteration residual. 35 ! Increase it to tighten convergence. 36 DOUBLE PRECISION, PARAMETER :: DEN = 1.0D1 !5.0D2 37 38 ! Dummy arguments 39 !---------------------------------------------------------------------// 40 ! Normalization factor for solids volume fraction correction residual. 41 ! At start of the iterate loop norms will either be 1 (i.e. not 42 ! normalized) or a user defined value given by norm_s. If norm_s 43 ! was set to zero then the normalization is based on dominate 44 ! term in the equation 45 DOUBLE PRECISION, INTENT(IN) :: NORMs 46 ! solids volume fraction correction residual 47 DOUBLE PRECISION, INTENT(OUT) :: RESs 48 ! Error index 49 INTEGER, INTENT(INOUT) :: IER 50 51 ! Local variables 52 !----------------------------------------------- 53 ! solids phase index locally assigned to mcp 54 ! mcp is the lowest index of those solids phases that are close_packed 55 ! and of the solids phase that is used for the solids correction 56 ! equation. 57 INTEGER :: M 58 ! Normalization factor for solids volume fraction correction residual 59 DOUBLE PRECISION :: NORMSloc 60 ! linear equation solver method and iterations 61 INTEGER :: LEQM, LEQI 62 63 ! temporary use of global arrays: 64 ! arraym1 (locally b_mmax) 65 ! vector B_M based on dominate term in correction equation 66 DOUBLE PRECISION :: B_MMAX(DIMENSION_3, DIMENSION_M) 67 ! Septadiagonal matrix A_m, vector B_m 68 ! DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M) 69 ! DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M) 70 !---------------------------------------------------------------------// 71 72 ! Note that currently this subroutine is only called when MMAX=1 AND 73 ! MCP is defined. This combintion effectively means solids phase 1 74 ! can close pack. 75 IF (MCP == UNDEFINED_I) THEN 76 ! this error should be caught earlier in the routines so that this 77 ! branch should never be entered 78 RETURN 79 ELSE 80 ! the lowest solids phase index of those solids phases that can close 81 ! pack (i.e. close_packed=T) and the index of the solids phase that is 82 ! used to form the solids correction equation. 83 M = MCP 84 ENDIF 85 call lock_ambm 86 87 ! Form the sparse matrix equation. Note that the index 0 is explicitly 88 ! used throughout this routine for creating the matrix equation. 89 ! However, the equation is based on the index of MCP. 90 91 ! initializing 92 CALL INIT_AB_M (A_M, B_M, IJKMAX2, 0) 93 EPP(:) = ZERO 94 95 CALL CONV_SOURCE_EPP (A_M, B_M, B_mmax, M) 96 97 ! Add point source contributions. 98 IF(POINT_SOURCE) CALL POINT_SOURCE_EPP (B_M, B_mmax, M) 99 100 ! Add usr source contributions 101 IF(CALL_USR_SOURCE(2)) CALL CALC_USR_SOURCE(SOLIDS_CORRECTION, & 102 A_M, B_M, lB_MMAX=B_MMAX, lM=M) 103 104 ! call check_ab_m(a_m, b_m, 0, .false., ier) 105 ! call write_ab_m(a_m, b_m, ijkmax2, 0, ier) 106 ! call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), 1, & 107 ! DO_K, ier) 108 109 110 ! Find average residual, maximum residual and location 111 NORMSloc = NORMS 112 IF(NORMS == ZERO) THEN 113 ! calculate the residual based on dominate term in correction equation 114 ! and use this to form normalization factor 115 CALL CALC_RESID_PP (B_MMAX, ONE, NUM_RESID(RESID_P,M), & 116 DEN_RESID(RESID_P,M), RESID(RESID_P,M), & 117 MAX_RESID(RESID_P,M), IJK_RESID(RESID_P,M)) 118 NORMSloc = RESID(RESID_P,M)/DEN 119 ENDIF 120 121 CALL CALC_RESID_PP (B_M, NORMSloc, NUM_RESID(RESID_P,M), & 122 DEN_RESID(RESID_P,M), RESID(RESID_P,M), MAX_RESID(RESID_P,M), & 123 IJK_RESID(RESID_P,M)) 124 RESS = RESID(RESID_P,M) 125 ! write(*,*) resid(resid_p, 1), max_resid(resid_p, 1), & 126 ! ijk_resid(resid_p, 1) 127 128 129 ! Solve EP_s_prime equation 130 CALL ADJUST_LEQ(RESID(RESID_P,M), LEQ_IT(2), LEQ_METHOD(2),& 131 LEQI, LEQM) 132 ! note index 0 is used here since that is the index that was used for 133 ! creating this matrix equation 134 CALL SOLVE_LIN_EQ ('EPp', 2, EPP, A_M, B_M, 0, LEQI, LEQM, & 135 LEQ_SWEEP(2), LEQ_TOL(2), LEQ_PC(2), IER) 136 137 ! call out_array(EPp, 'EPp') 138 139 call unlock_ambm 140 141 RETURN 142 END SUBROUTINE SOLVE_EPP 143