File: N:\mfix\model\des\diffuse_mean_fields.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: DIFFUSE_MEAN_FIELDS                                     !
4     !  Author: J.Musser                                   Date: 11-NOV-14  !
5     !                                                                      !
6     !  Purpose: Given the field variable PHI (e.g., volume fraction),      !
7     !  diffuse it across the Eulerian grid such that the Full Width at     !
8     !  Half Maximum (FWHM) equals the user specified diffusion width.      !
9     !                                                                      !
10     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
11           SUBROUTINE DIFFUSE_MEAN_FIELD(PHI, VNAME)
12     
13     ! Global Variables:
14     !---------------------------------------------------------------------//
15     ! Max bound for array sizes.
16           use geometry, only: IJKMAX2
17     ! Coefficient matrix and force vector.
18           use ambm, only: A_M, B_M
19     ! Method to solve linear system and max iterations
20           use leqsol, only: LEQ_METHOD, LEQ_IT
21     ! Preconditioner, sweep method, convergence tolerance
22           use leqsol, only: LEQ_PC, LEQ_SWEEP, LEQ_TOL
23     ! Size of fluid variable arrays.
24           use param, only: DIMENSION_3
25     
26     ! Module procedures
27     !---------------------------------------------------------------------//
28     ! Routines to manage messages to user.
29           use error_manager
30           use machine, only: wall_time
31     
32           IMPLICIT NONE
33     
34     ! Dummy Arguments:
35     !---------------------------------------------------------------------//
36     ! Variable to diffuse
37           DOUBLE PRECISION, INTENT(INOUT) :: PHI(DIMENSION_3)
38     ! Name of variable to diffuse
39           CHARACTER(LEN=*), INTENT(IN) :: VNAME
40     
41     ! Local Variables:
42     !---------------------------------------------------------------------//
43     ! Integer error flag
44           INTEGER :: IER
45     ! Linear equation solver method and iterations
46           INTEGER :: LEQM, LEQI
47     ! Start, stop and step size of diffusion time march
48           DOUBLE PRECISION :: DIF_TIME, DIF_STOP, DIF_DT
49     ! wall time at start
50           DOUBLE PRECISION :: WALL_START
51     ! Local flag to print debug messages
52           LOGICAL, PARAMETER :: setDBG = .TRUE.
53           DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: DIF
54     !......................................................................!
55     
56           IF(setDBG) THEN
57              WRITE(ERR_MSG, "(/3x,'Diffusing Variable: ',A)") VNAME
58              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
59              WALL_START = WALL_TIME()
60           ENDIF
61     
62     ! Populate the diffusion coefficients
63           CALL CALC_DIF_DES(DIF, setDBG, IER)
64     
65           DIF_STOP = 1.0d0
66           DIF_TIME = 0.0d0
67           DIF_DT = DIF_STOP/5.0
68     
69     ! Integrate the diffusion equation (time, space)
70           DO WHILE(DIF_TIME < DIF_STOP)
71     ! Initialize the coefficient matrix and force vector
72              CALL INIT_AB_M (A_M, B_M, IJKMAX2, 0)
73     ! Calculate the coefficients
74              CALL DIF_PHI_DES(0, DIF, A_M, B_M)
75     ! Apply zero-flux BC at all walls
76              CALL DIF_PHI_BC_DES(PHI, 0, A_M, B_M)
77     ! Collect the center coefficient and force vector
78              CALL DIF_PHI_SOURCE_DES(PHI, 0, A_M, B_M, DIF_DT)
79     ! Set the local method and iterations.
80              CALL ADJUST_LEQ(0.0d0, LEQ_IT(10), LEQ_METHOD(10), LEQI, LEQM)
81     ! Solve the linear system.
82              CALL SOLVE_LIN_EQ (VNAME, 10, PHI, A_M, B_M, 0, LEQI, LEQM, &
83                 LEQ_SWEEP(10), LEQ_TOL(10), LEQ_PC(10), IER)
84     ! Advance time.
85              DIF_TIME = DIF_TIME + DIF_DT
86           ENDDO
87     
88     ! Debugging information
89           IF(setDBG) THEN
90              WRITE(ERR_MSG, 9001) WALL_TIME() - WALL_START
91              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
92           ENDIF
93     
94      9001 FORMAT(5x,'Wall Time: ',g11.4)
95     
96           RETURN
97           END SUBROUTINE DIFFUSE_MEAN_FIELD
98     
99     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
100     !                                                                      !
101     !  Subroutine: CALC_DIF_DES                                            !
102     !  Author: J.Musser                                   Date: 11-NOV-14  !
103     !                                                                      !
104     !  Purpose: Calculate the diffusion coefficient for diffusing mean     !
105     !  fields. Presently the diffusion coefficient is constant.            !
106     !                                                                      !
107     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
108           SUBROUTINE CALC_DIF_DES(DIF, lDBG, IER)
109     
110           use param
111           use param1
112     
113           use compar, only: IJKStart3, IJKEnd3
114           use particle_filter, only: DES_DIFFUSE_WIDTH
115           use functions, only: FLUID_AT
116     
117           use error_manager
118     
119           IMPLICIT NONE
120     
121           DOUBLE PRECISION, INTENT(INOUT) :: DIF(DIMENSION_3)
122           LOGICAL, INTENT(IN) :: lDBG
123           INTEGER, INTENT(INOUT) :: IER
124     
125     ! Fluid Cell indices
126           INTEGER :: IJK
127     
128           DOUBLE PRECISION :: lDIF
129     
130           IER = 0
131     
132     ! The diffusion coefficient is set so that over one second, the
133     ! quantity diffuses such that the Full Width at Half Maximum (FWHM)
134     ! equals what the user specified as the "filter wideth."
135           lDIF = ((0.5*DES_DIFFUSE_WIDTH)**2) / &
136              (2.0*sqrt(2.0*log(2.0)))
137     
138     ! Store the diffusion coefficient in all fluid cells.
139           DO IJK = IJKStart3, IJKEnd3
140              DIF(IJK) = ZERO
141              IF(FLUID_AT(IJK)) DIF(IJK) = lDIF
142           ENDDO
143     
144     ! Information included for debugging.
145           IF(lDBG) THEN
146              WRITE(ERR_MSG, 9100) iVal(lDIF)
147              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
148           ENDIF
149     
150      9100 FORMAT(/3x,'Diffusion Coefficient: ',A)
151     
152           RETURN
153           END SUBROUTINE CALC_DIF_DES
154