MFIX  2016-1
calc_grad_des.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2 ! !
3 ! Subroutine: CALC_PG_GRAD !
4 ! Purpose: Calculate cell centered pressure force exerted on the !
5 ! particles in the cell by the gas/fluid phase !
6 ! Note that P_force is evaluated as -dp/dx !
7 ! !
8 ! Notes: This pressure force only needs to be calculated once during !
9 ! the DEM loop (at the beginning) since the gas/fluid phase !
10 ! is essentially static at that point (i.e., gas field is not !
11 ! updated during DEM loop !
12 ! !
13 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
14  SUBROUTINE calc_grad_des(PHI, DEL_PHI)
15 
16 ! Flag for cut-cells
17  use cutcell, only: cartesian_grid
18 
19  use param, only: dimension_3
20 
21  IMPLICIT NONE
22 
23 ! Dummy Arguments:
24 !---------------------------------------------------------------------//
25  DOUBLE PRECISION, INTENT(IN) :: PHI(dimension_3)
26  DOUBLE PRECISION, INTENT(OUT) :: DEL_PHI(3,dimension_3)
27 !......................................................................!
28 
29  IF(cartesian_grid) THEN
30  CALL calc_grad_des_cg(phi, del_phi)
31  ELSE
32  CALL calc_grad_des_std(phi, del_phi)
33  ENDIF
34 
35  RETURN
36  END SUBROUTINE calc_grad_des
37 
38 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
39 ! !
40 ! Subroutine: CALC_GRAD_DES_STD !
41 ! !
42 ! Purpose: Calculate cell centered gradient (DEL_PHI) of scalar PHI. !
43 ! This routine calculates the average scalar value at the !
44 ! cell faces to calculate the graident across the cell. !
45 ! !
46 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
47  SUBROUTINE calc_grad_des_std(PHI, DEL_PHI)
48 
49 ! Modules
50 !-----------------------------------------------
51 
52  use geometry, only: do_k
53  use geometry, only: odx, ody, odz
54  USE geometry, only: imin1, jmin1, kmin1
55  USE geometry, only: imax1, jmax1, kmax1
56  USE indices, only: i_of, j_of, k_of
57  USE compar, only: ijkstart3, ijkend3
58 
59  use functions, only: fluid_at
60  use functions, only: im_of, jm_of, km_of
61  use functions, only: ip_of, jp_of, kp_of
62  use fun_avg, only: avg_x, avg_y, avg_z
63 
64  USE param, only: dimension_3
65  USE param1, only: zero
66 
67  IMPLICIT NONE
68 
69 ! Dummy Arguments:
70 !---------------------------------------------------------------------//
71  DOUBLE PRECISION, INTENT(IN) :: PHI(dimension_3)
72  DOUBLE PRECISION, INTENT(OUT) :: DEL_PHI(3,dimension_3)
73 
74 ! Local variables
75 !---------------------------------------------------------------------//
76 ! general i, j, k indices
77  INTEGER :: I, J, K, IJK
78  INTEGER :: IPJK, IJPK, IJKP, IMJK, IJMK, IJKM
79 !......................................................................!
80 
81 
82  DO ijk = ijkstart3, ijkend3
83  del_phi(:,ijk) = zero
84  IF(.NOT.fluid_at(ijk)) cycle
85 
86  i = i_of(ijk)
87  imjk = im_of(ijk)
88  ipjk = ip_of(ijk)
89 
90  IF(imin1 == imax1) THEN
91  ELSEIF((i>imin1).AND.(i<imax1)) THEN
92  del_phi(1,ijk) = odx(i)*(avg_x(phi(ijk),phi(ipjk),i) - &
93  avg_x(phi(imjk),phi(ijk),i-1))
94  ELSEIF(i == imin1) THEN
95  del_phi(1,ijk) = 2.0d0*odx(i) * &
96  (avg_x(phi(ijk),phi(ipjk),i) - phi(ijk))
97  ELSEIF(i == imax1) THEN
98  del_phi(1,ijk) = 2.0d0*odx(i) * &
99  (phi(ijk) - avg_x(phi(imjk), phi(ijk), i-1))
100  ELSE
101  del_phi(1,ijk) = zero
102  ENDIF
103 
104 
105  j = j_of(ijk)
106  ijmk = jm_of(ijk)
107  ijpk = jp_of(ijk)
108 
109  IF(jmin1 == jmax1) THEN
110  ELSEIF((j>jmin1) .AND. (j<jmax1)) THEN
111  del_phi(2,ijk) = ody(j)*(avg_y(phi(ijk),phi(ijpk),j) - &
112  avg_y(phi(ijmk),phi(ijk),j-1))
113  ELSEIF(j == jmin1) THEN
114  del_phi(2,ijk) = 2.0d0*ody(j) * &
115  (avg_y(phi(ijk),phi(ijpk),j) - phi(ijk))
116  ELSEIF(j == jmax1) THEN
117  del_phi(2,ijk) = 2.0d0*ody(j) * &
118  (phi(ijk)- avg_y(phi(ijmk),phi(ijk),j-1))
119  ELSE
120  del_phi(2,ijk) = zero
121  ENDIF
122 
123  IF(do_k) THEN
124 
125  k = k_of(ijk)
126  ijkm = km_of(ijk)
127  ijkp = kp_of(ijk)
128 
129  IF(kmin1 == kmax1) THEN
130  ELSEIF((k>kmin1) .AND. (k<kmax1)) THEN
131  del_phi(3,ijk) = odz(k)*(avg_z(phi(ijk),phi(ijkp),k) - &
132  avg_z(phi(ijkm),phi(ijk),k-1))
133  ELSEIF(k == kmin1) THEN
134  del_phi(3,ijk) = 2.0d0*odz(k) * &
135  (avg_z(phi(ijk),phi(ijkp),k) - phi(ijk))
136  ELSEIF(k == kmax1) THEN
137  del_phi(3,ijk) = 2.0d0*odz(k) * &
138  (phi(ijk) - avg_z(phi(ijkm),phi(ijk),k-1))
139  ENDIF
140  ENDIF
141  ENDDO
142 
143  RETURN
144  END SUBROUTINE calc_grad_des_std
145 
146 
147 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
148 ! !
149 ! Subroutine: CALC_GRAD_DES_CG !
150 ! !
151 ! Purpose: Calculate cell centered gradient (DEL_PHI) of scalar PHI. !
152 ! This routine calculates the graidents at the cell face !
153 ! and averages the face values to the cell center. !
154 ! !
155 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
156  SUBROUTINE calc_grad_des_cg(PHI, DEL_PHI)
157 
158 ! Modules
159 !---------------------------------------------------------------------//
160  use geometry, only: do_k
161  use geometry, only: dx, dy, dz
162  use indices, only: i_of, j_of, k_of
163  use compar, only: ijkstart3, ijkend3
164 
165  use functions, only: fluid_at
166  use functions, only: im_of, jm_of, km_of
167  use functions, only: ip_of, jp_of, kp_of
168 
169  use param, only: dimension_3
170  use param1, only: zero
171 
172  IMPLICIT NONE
173 
174 ! Dummy Arguments:
175 !---------------------------------------------------------------------//
176  DOUBLE PRECISION, INTENT(IN) :: PHI(dimension_3)
177  DOUBLE PRECISION, INTENT(OUT) :: DEL_PHI(3,dimension_3)
178 
179 ! Local variables
180 !---------------------------------------------------------------------//
181 ! General i, j, k indices
182  INTEGER :: I, J, K, IJK
183  INTEGER :: IJKE, IJKW, IJKN, IJKS, IJKT, IJKB
184 ! Counter
185  DOUBLE PRECISION :: dLC
186 !......................................................................!
187 
188 
189  DO ijk = ijkstart3, ijkend3
190 
191  i = i_of(ijk)
192  j = j_of(ijk)
193  k = k_of(ijk)
194  del_phi(:,ijk) = zero
195 
196  IF(.NOT.fluid_at(ijk)) cycle
197 
198  dlc = 0.0d0
199  ijke = ip_of(ijk)
200  ijkw = im_of(ijk)
201 
202  IF(fluid_at(ijke)) THEN
203  del_phi(1,ijk) = del_phi(1,ijk) + &
204  2.0d0*(phi(ijke) - phi(ijk))/(dx(i) + dx(i_of(ijke)))
205  dlc = dlc + 1.0d0
206  ENDIF
207  IF(fluid_at(ijkw)) THEN
208  del_phi(1,ijk) = del_phi(1,ijk) + &
209  2.0d0*(phi(ijk) - phi(ijkw))/(dx(i) + dx(i_of(ijkw)))
210  dlc = dlc + 1.0d0
211  ENDIF
212  del_phi(1,ijk) = del_phi(1,ijk)/max(1.0d0,dlc)
213 
214 
215  dlc = 0.0d0
216  ijkn = jp_of(ijk)
217  ijks = jm_of(ijk)
218 
219  IF(fluid_at(ijkn)) THEN
220  del_phi(2,ijk) = del_phi(2,ijk) + &
221  2.0d0*(phi(ijkn) - phi(ijk))/(dy(j) + dy(j_of(ijkn)))
222  dlc = dlc + 1.0d0
223  ENDIF
224 
225  IF(fluid_at(ijks)) THEN
226  del_phi(2,ijk) = del_phi(2,ijk) + &
227  2.d0*(phi(ijk) - phi(ijks))/(dy(j) + dy(j_of(ijks)))
228  dlc = dlc + 1.0d0
229  ENDIF
230  del_phi(2,ijk) = del_phi(2,ijk)/max(1.0d0, dlc)
231 
232  IF(do_k) THEN
233 
234  dlc = 0.0d0
235  ijkt = kp_of(ijk)
236  ijkb = km_of(ijk)
237 
238  IF(fluid_at(ijkt)) THEN
239  del_phi(3,ijk) = del_phi(3,ijk) + &
240  2.0d0*(phi(ijkt) - phi(ijk))/(dz(k) + dz(k_of(ijkt)))
241  dlc = dlc + 1.0d0
242  ENDIF
243  IF(fluid_at(ijkb)) THEN
244  del_phi(3,ijk) = del_phi(3,ijk) + &
245  2.0d0*(phi(ijk) - phi(ijkb))/(dz(k) + dz(k_of(ijkb)))
246  dlc = dlc + 1.0d0
247  ENDIF
248  del_phi(3,ijk) = del_phi(3,ijk)/max(1.0d0, dlc)
249  ENDIF
250  ENDDO
251 
252  RETURN
253  END SUBROUTINE calc_grad_des_cg
indices::i_of
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
calc_grad_des_cg
subroutine calc_grad_des_cg(PHI, DEL_PHI)
Definition: calc_grad_des.f:157
geometry::ody
double precision, dimension(:), allocatable ody
Definition: geometry_mod.f:116
geometry::odx
double precision, dimension(:), allocatable odx
Definition: geometry_mod.f:114
param::dimension_3
integer dimension_3
Definition: param_mod.f:11
geometry::dy
double precision, dimension(0:dim_j) dy
Definition: geometry_mod.f:70
geometry::dz
double precision, dimension(0:dim_k) dz
Definition: geometry_mod.f:72
indices::k_of
integer, dimension(:), allocatable k_of
Definition: indices_mod.f:47
fun_avg
Definition: fun_avg_mod.f:1
geometry::kmax1
integer kmax1
Definition: geometry_mod.f:58
indices::j_of
integer, dimension(:), allocatable j_of
Definition: indices_mod.f:46
geometry::imax1
integer imax1
Definition: geometry_mod.f:54
geometry::dx
double precision, dimension(0:dim_i) dx
Definition: geometry_mod.f:68
geometry::jmax1
integer jmax1
Definition: geometry_mod.f:56
param
Definition: param_mod.f:2
cutcell::cartesian_grid
logical cartesian_grid
Definition: cutcell_mod.f:13
geometry::jmin1
integer jmin1
Definition: geometry_mod.f:42
geometry::do_k
logical do_k
Definition: geometry_mod.f:30
geometry::odz
double precision, dimension(:), allocatable odz
Definition: geometry_mod.f:118
compar::ijkstart3
integer ijkstart3
Definition: compar_mod.f:80
calc_grad_des
subroutine calc_grad_des(PHI, DEL_PHI)
Definition: calc_grad_des.f:15
geometry::imin1
integer imin1
Definition: geometry_mod.f:40
geometry::kmin1
integer kmin1
Definition: geometry_mod.f:44
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
calc_grad_des_std
subroutine calc_grad_des_std(PHI, DEL_PHI)
Definition: calc_grad_des.f:48