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