1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 2 ! C 3 ! Subroutine: SOLVE_CONTINUITY C 4 ! Purpose: Solve for solids bulk density (i.e., material density C 5 ! multiplied by volume fraction). C 6 ! C 7 ! Author: M. Syamlal Date: 2-JUL-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_CONTINUITY(M,IER) 19 20 !----------------------------------------------- 21 ! Modules 22 !----------------------------------------------- 23 USE param 24 USE param1 25 USE physprop 26 USE geometry 27 USE fldvar 28 USE indices 29 USE residual 30 USE cont 31 USE leqsol 32 Use ambm 33 USE ur_facs 34 use run 35 use ps 36 37 IMPLICIT NONE 38 !----------------------------------------------- 39 ! Dummy arguments 40 !----------------------------------------------- 41 ! phase index 42 INTEGER, INTENT(IN) :: M 43 ! error index 44 INTEGER, INTENT(INOUT) :: IER 45 !----------------------------------------------- 46 ! Local variables 47 !----------------------------------------------- 48 ! solids volume fraction residual 49 DOUBLE PRECISION :: RESs 50 ! linear equation solver method and iterations 51 INTEGER :: LEQM, LEQI 52 53 ! temporary use of global arrays: 54 ! Septadiagonal matrix A_m, vector B_m 55 ! DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M) 56 ! DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M) 57 !----------------------------------------------- 58 59 call lock_ambm 60 61 62 IF (M==0) THEN 63 ! solve gas phase continuity equation. note that this branch will never 64 ! be entered given that the calling subroutine (iterate) only calls 65 ! solve_continuity when M>1 66 67 ! initializing 68 CALL INIT_AB_M (A_M, B_M, IJKMAX2, 0, IER) 69 70 ! forming the matrix equation 71 CALL CONV_ROP_G (A_M, B_M, IER) 72 CALL SOURCE_ROP_G (A_M, B_M, IER) 73 74 ! calculating the residual 75 CALL CALC_RESID_C (ROP_G, A_M, B_M, 0, NUM_RESID(RESID_RO,0), & 76 DEN_RESID(RESID_RO,0), RESID(RESID_RO,0), MAX_RESID(& 77 RESID_RO,0), IJK_RESID(RESID_RO,0), IER) 78 79 ! call check_ab_m(a_m, b_m, 0, .true., ier) 80 ! call write_ab_m(a_m, b_m, ijkmax2, 0, ier) 81 ! write(*,*) resid(resid_ro, 0), max_resid(resid_ro, 0), & 82 ! ijk_resid(resid_ro, 0) 83 ! call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0),& 84 ! 1, DO_K, ier) 85 86 ! solving gas continuity 87 CALL ADJUST_LEQ (RESID(RESID_RO,0), LEQ_IT(2), LEQ_METHOD(2),& 88 LEQI, LEQM) 89 CALL SOLVE_LIN_EQ ('ROP_g', 2, ROP_G, A_M, B_M, 0, LEQI, & 90 LEQM, LEQ_SWEEP(2), LEQ_TOL(2), LEQ_PC(2), IER) 91 CALL ADJUST_ROP (ROP_G) 92 ! call out_array(ROP_g, 'rop_g') 93 94 ELSE 95 ! solve solids phase M continuity equation. 96 97 ! initializing 98 CALL INIT_AB_M (A_M, B_M, IJKMAX2, M, IER) 99 100 ! forming the matrix equation 101 CALL CONV_ROP_S (A_M, B_M, M, IER) 102 CALL SOURCE_ROP_S (A_M, B_M, M, IER) 103 IF(POINT_SOURCE) CALL POINT_SOURCE_ROP_S (B_M, M, IER) 104 105 CALL CALC_RESID_C (ROP_S(1,M), A_M, B_M, M, & 106 NUM_RESID(RESID_RO,M), DEN_RESID(RESID_RO,M), & 107 RESID(RESID_RO,M), MAX_RESID(RESID_RO,M), & 108 IJK_RESID(RESID_RO,M), IER) 109 RESS = RESID(RESID_RO,M) 110 111 ! call check_ab_m(a_m, b_m, m, .true., ier) 112 ! write(*,*)'solve_cont=',resid(resid_ro, m),max_resid(resid_ro, m),& 113 ! ijk_resid(resid_ro, m),m 114 ! call write_ab_m(a_m, b_m, ijkmax2, m, ier) 115 ! call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, M), 1, 116 ! & DO_K, ier) 117 ! 118 CALL ADJUST_LEQ (RESID(RESID_RO,M), LEQ_IT(2), LEQ_METHOD(2),& 119 LEQI, LEQM) 120 CALL SOLVE_LIN_EQ ('ROP_s', 2, ROP_S(1,M), A_M, B_M, M, LEQI,& 121 LEQM,LEQ_SWEEP(2), LEQ_TOL(2), LEQ_PC(2), IER) 122 CALL ADJUST_ROP (ROP_S(1,M)) 123 ! call out_array(rop_s(1,m), 'rop_s') 124 ENDIF 125 126 call unlock_ambm 127 128 129 RETURN 130 END SUBROUTINE SOLVE_CONTINUITY 131