MFIX  2016-1
qmomk_time_march.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Module name: QMOMK_TIME_MARCH C
4 ! Purpose: Called in time_march.f to do QMOMK calculations C
5 ! C
6 ! Author: Alberto Passalacqua Date: C
7 ! Reviewer: Date: C
8 ! C
9 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
10 !
11 SUBROUTINE qmomk_time_march
12 
13  USE param
14  USE param1
15  USE constant
16  USE run
17  USE output
18  USE physprop
19  USE fldvar
20  USE geometry
21  USE cont
22  USE tau_g
23  USE tau_s
24  USE visc_g
25  USE visc_s
26  USE funits
27  USE vshear
28  USE scalars
29  USE drag
30  USE rxns
31  USE compar
32  USE time_cpu
33  USE is
34  USE indices
35  USE sendrecv
36  USE qmom_kinetic_equation
37  USE qmomk_fluxes
38  USE qmomk_quadrature
39  USE qmomk_collision
40  USE qmomk_parameters
41  USE drag
42  USE ur_facs
43  USE fun_avg
44  USE functions
45 
46  IMPLICIT NONE
47 
48  INTEGER :: I,J,K,M,M2,IN
49  INTEGER :: IJK, IPJK, IJPK, IJKP, IMJK, IJMK, IJKM
50  DOUBLE PRECISION :: Vmax, min_space_delta
51  DOUBLE PRECISION :: vrel, Rep, CD, beta_drag, min_tau_drag
52  DOUBLE PRECISION :: UGC, VGC, WGC, mom0, QMOMK_TCOL, drag_exp, QMOMK_OMEGA
53  DOUBLE PRECISION QMOMK_TIME, QMOMK_DT_TMP, TMP_DTS
54  INTEGER :: TIME_FACTOR, TIME_I
55  LOGICAL, SAVE :: FIRST_PASS = .true.
56 
57  DOUBLE PRECISION, DIMENSION(QMOMK_NN) :: Nlminus, Nlplus, Nrminus, Nrplus
58  DOUBLE PRECISION, DIMENSION(QMOMK_NN) :: Ulminus, Ulplus, Urminus, Urplus
59  DOUBLE PRECISION, DIMENSION(QMOMK_NN) :: Vlminus, Vlplus, Vrminus, Vrplus
60  DOUBLE PRECISION, DIMENSION(QMOMK_NN) :: Wlminus, Wlplus, Wrminus, Wrplus
61  DOUBLE PRECISION, DIMENSION(QMOMK_NMOM) :: F_x_left, F_x_right, F_y_left
62  DOUBLE PRECISION, DIMENSION(QMOMK_NMOM) :: F_y_right, F_z_left, F_z_right
63 
64  IF (first_pass) THEN
65  first_pass = .false.
66  ! Setting initial guess for drag time
67  qmomk_time = 0.d0
68  qmomk_tau_drag = 1.0d-5
69 
70  ! A.P. For restart runs, IC's are already set and would be overwritten
71  IF (run_type == 'NEW') THEN
72  CALL qmomk_initial_conditions
73  CALL qmomk_init_bc
74  ELSE
75  CALL qmomk_init_bc
76  ENDIF
77 
78  qmomk_m1 = qmomk_m0
79  qmomk_n1 = qmomk_n0
80  qmomk_u1 = qmomk_u0
81  qmomk_v1 = qmomk_v0
82  qmomk_w1 = qmomk_w0
83  END IF
84 
85  ! Finding initial QMOMK time step
86  ! CFL condition
87  vmax = 0.d0
88  DO m = 1, mmax
89  DO ijk = ijkstart3, ijkend3
90  DO i = 1, qmomk_nn
91  IF (vmax < max(abs(qmomk_u0(i,ijk,m)), abs(qmomk_v0(i,ijk,m)), abs(qmomk_w0(i,ijk,m)))) THEN
92  vmax = max(abs(qmomk_u0(i,ijk,m)), abs(qmomk_v0(i,ijk,m)), abs(qmomk_w0(i,ijk,m)))
93  END IF
94  END DO
95  CALL compute_collision_time(qmomk_m0(:,ijk,m), d_p0(m), theta_m(ijk,m), qmomk_collision_time(ijk,m), dt)
96  END DO
97  END DO
98 
99  min_space_delta = min(minval(dx), minval(dy), minval(dz))
100 
101  qmomk_dt = qmomk_cfl*min_space_delta/vmax
102 
103  ! Drag time and collision time
104  min_tau_drag = 1.d0
105  qmomk_tcol = 1.d10
106  DO ijk = ijkstart3, ijkend3
107  IF (fluid_at(ijk)) THEN
108  IF (min_tau_drag > minval(qmomk_tau_drag(:,ijk,:))) THEN
109  min_tau_drag = minval(qmomk_tau_drag(:,ijk,:))
110  END IF
111  IF (qmomk_tcol > minval(qmomk_collision_time(ijk,:))) THEN
112  qmomk_tcol = minval(qmomk_collision_time(ijk,:))
113  END IF
114  END IF
115  END DO
116 
117  IF ( qmomk_collisions /= 'NONE') THEN
118  qmomk_dt = min(qmomk_tcol,qmomk_dt)
119  END IF
120 
121  qmomk_dt = min(qmomk_dt,min_tau_drag/10.d0)
122 
123  ! Preparing for QMOMK internal time stepping, if required
124  qmomk_time = time
125  IF (dt >= qmomk_dt) THEN
126  time_factor = ceiling(real(dt/qmomk_dt)) + 1
127  qmomk_dt_tmp = zero
128  ELSE
129  time_factor = 1
130  qmomk_dt_tmp = qmomk_dt
131  qmomk_dt = dt
132  END IF
133 
134  tmp_dts = zero
135 
136  DO time_i = 1, time_factor ! Starting QMOMK internal time stepping
137 
138  ! Exit if time is greater than flow time
139  IF (qmomk_time .GE. (time + dt)) EXIT
140 
141  IF ((qmomk_time + qmomk_dt) .GT. (time + dt)) THEN
142  tmp_dts = qmomk_dt
143  qmomk_dt = time + dt - qmomk_time
144  END IF
145 
146  print *,'Flow time step: ',dt
147  print *,'QMOM DT = ',qmomk_dt
148  print *,'QMOMK internal time step: ',time_i
149  print *,'Time factor', time_factor
150 
151  qmomk_omega = (1.d0 + c_e)/2.d0
152 
153  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
154  ! First step of Runge-Kutta time integration !
155  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
156  ! Calculating kinetic based fluxes (Half time step)
157  DO m = 1, mmax
158  DO ijk = ijkstart3, ijkend3
159  IF (fluid_at(ijk)) THEN
160  i = i_of(ijk)
161  j = j_of(ijk)
162  k = k_of(ijk)
163 
164  imjk = im_of(ijk)
165  ijmk = jm_of(ijk)
166  ijkm = km_of(ijk)
167 
168  ipjk = ip_of(ijk)
169  ijpk = jp_of(ijk)
170  ijkp = kp_of(ijk)
171 
172  ! 2D Case
173  IF (no_k) THEN
174  ! x - direction
175  nlminus = qmomk_n0(:,imjk,m)
176  ulminus = qmomk_u0(:,imjk,m)
177  vlminus = qmomk_v0(:,imjk,m)
178  wlminus = qmomk_w0(:,imjk,m)
179 
180  nlplus = qmomk_n0(:,ijk,m)
181  ulplus = qmomk_u0(:,ijk,m)
182  vlplus = qmomk_v0(:,ijk,m)
183  wlplus = qmomk_w0(:,ijk,m)
184 
185  nrminus = qmomk_n0(:,ijk,m)
186  urminus = qmomk_u0(:,ijk,m)
187  vrminus = qmomk_v0(:,ijk,m)
188  wrminus = qmomk_w0(:,ijk,m)
189 
190  nrplus = qmomk_n0(:,ipjk,m)
191  urplus = qmomk_u0(:,ipjk,m)
192  vrplus = qmomk_v0(:,ipjk,m)
193  wrplus = qmomk_w0(:,ipjk,m)
194 
195  CALL kinetic_flux_x_twenty_eight_nodes (nlminus, ulminus, vlminus, wlminus, &
196  nlplus, ulplus, vlplus, wlplus, f_x_left)
197  CALL kinetic_flux_x_twenty_eight_nodes (nrminus, urminus, vrminus, wrminus, &
198  nrplus, urplus, vrplus, wrplus, f_x_right)
199 
200  ! y - direction
201  nlminus = qmomk_n0(:,ijmk,m)
202  ulminus = qmomk_u0(:,ijmk,m)
203  vlminus = qmomk_v0(:,ijmk,m)
204  wlminus = qmomk_w0(:,ijmk,m)
205 
206  nlplus = qmomk_n0(:,ijk,m)
207  ulplus = qmomk_u0(:,ijk,m)
208  vlplus = qmomk_v0(:,ijk,m)
209  wlplus = qmomk_w0(:,ijk,m)
210 
211  nrminus = qmomk_n0(:,ijk,m)
212  urminus = qmomk_u0(:,ijk,m)
213  vrminus = qmomk_v0(:,ijk,m)
214  wrminus = qmomk_w0(:,ijk,m)
215 
216  nrplus = qmomk_n0(:,ijpk,m)
217  urplus = qmomk_u0(:,ijpk,m)
218  vrplus = qmomk_v0(:,ijpk,m)
219  wrplus = qmomk_w0(:,ijpk,m)
220 
221  CALL kinetic_flux_y_twenty_eight_nodes (nlminus, ulminus, vlminus, wlminus, &
222  nlplus, ulplus, vlplus, wlplus, f_y_left)
223  CALL kinetic_flux_y_twenty_eight_nodes (nrminus, urminus, vrminus, wrminus, &
224  nrplus, urplus, vrplus, wrplus, f_y_right)
225 
226  qmomk_m1(:,ijk,m) = qmomk_m0(:,ijk,m) &
227  - 0.5*(qmomk_dt/minval(dx))*(f_x_right - f_x_left) - 0.5*(qmomk_dt/minval(dy))*(f_y_right - f_y_left)
228 
229 
230  ELSE ! 3D case
231  ! x - direction
232  nlminus = qmomk_n0(:,imjk,m)
233  ulminus = qmomk_u0(:,imjk,m)
234  vlminus = qmomk_v0(:,imjk,m)
235  wlminus = qmomk_w0(:,imjk,m)
236 
237  nlplus = qmomk_n0(:,ijk,m)
238  ulplus = qmomk_u0(:,ijk,m)
239  vlplus = qmomk_v0(:,ijk,m)
240  wlplus = qmomk_w0(:,ijk,m)
241 
242  nrminus = qmomk_n0(:,ijk,m)
243  urminus = qmomk_u0(:,ijk,m)
244  vrminus = qmomk_v0(:,ijk,m)
245  wrminus = qmomk_w0(:,ijk,m)
246 
247  nrplus = qmomk_n0(:,ipjk,m)
248  urplus = qmomk_u0(:,ipjk,m)
249  vrplus = qmomk_v0(:,ipjk,m)
250  wrplus = qmomk_w0(:,ipjk,m)
251 
252  CALL kinetic_flux_x_twenty_eight_nodes (nlminus, ulminus, vlminus, wlminus, &
253  nlplus, ulplus, vlplus, wlplus, f_x_left)
254  CALL kinetic_flux_x_twenty_eight_nodes (nrminus, urminus, vrminus, wrminus, &
255  nrplus, urplus, vrplus, wrplus, f_x_right)
256 
257  ! y - direction
258  nlminus = qmomk_n0(:,ijmk,m)
259  ulminus = qmomk_u0(:,ijmk,m)
260  vlminus = qmomk_v0(:,ijmk,m)
261  wlminus = qmomk_w0(:,ijmk,m)
262 
263  nlplus = qmomk_n0(:,ijk,m)
264  ulplus = qmomk_u0(:,ijk,m)
265  vlplus = qmomk_v0(:,ijk,m)
266  wlplus = qmomk_w0(:,ijk,m)
267 
268  nrminus = qmomk_n0(:,ijk,m)
269  urminus = qmomk_u0(:,ijk,m)
270  vrminus = qmomk_v0(:,ijk,m)
271  wrminus = qmomk_w0(:,ijk,m)
272 
273  nrplus = qmomk_n0(:,ijpk,m)
274  urplus = qmomk_u0(:,ijpk,m)
275  vrplus = qmomk_v0(:,ijpk,m)
276  wrplus = qmomk_w0(:,ijpk,m)
277 
278  CALL kinetic_flux_y_twenty_eight_nodes (nlminus, ulminus, vlminus, wlminus, &
279  nlplus, ulplus, vlplus, wlplus, f_y_left)
280  CALL kinetic_flux_y_twenty_eight_nodes (nrminus, urminus, vrminus, wrminus, &
281  nrplus, urplus, vrplus, wrplus, f_y_right)
282 
283  ! z - direction
284  nlminus = qmomk_n0(:,ijkm,m)
285  ulminus = qmomk_u0(:,ijkm,m)
286  vlminus = qmomk_v0(:,ijkm,m)
287  wlminus = qmomk_w0(:,ijkm,m)
288 
289  nlplus = qmomk_n0(:,ijk,m)
290  ulplus = qmomk_u0(:,ijk,m)
291  vlplus = qmomk_v0(:,ijk,m)
292  wlplus = qmomk_w0(:,ijk,m)
293 
294  nrminus = qmomk_n0(:,ijk,m)
295  urminus = qmomk_u0(:,ijk,m)
296  vrminus = qmomk_v0(:,ijk,m)
297  wrminus = qmomk_w0(:,ijk,m)
298 
299  nrplus = qmomk_n0(:,ijkp,m)
300  urplus = qmomk_u0(:,ijkp,m)
301  vrplus = qmomk_v0(:,ijkp,m)
302  wrplus = qmomk_w0(:,ijkp,m)
303 
304  CALL kinetic_flux_z_twenty_eight_nodes (nlminus, ulminus, vlminus, wlminus, &
305  nlplus, ulplus, vlplus, wlplus, f_z_left)
306  CALL kinetic_flux_z_twenty_eight_nodes (nrminus, urminus, vrminus, wrminus, &
307  nrplus, urplus, vrplus, wrplus, f_z_right)
308 
309  qmomk_m1(:,ijk,m) = qmomk_m0(:,ijk,m) - 0.5*(qmomk_dt/minval(dx))*(f_x_right - f_x_left) - &
310  0.5*(qmomk_dt/minval(dy))*(f_y_right - f_y_left) - 0.5*(qmomk_dt/minval(dz))*(f_z_right - f_z_left)
311  END IF
312  END IF
313  END DO
314  END DO
315  ! End of fluxes calcuation
316 
317  ! Update weights and abscissas
318  DO m = 1, mmax
319  DO ijk = ijkstart3, ijkend3
320  IF (fluid_at(ijk)) THEN
321  CALL eight_node_3d (qmomk_m1(:,ijk,m), qmomk_n1(:,ijk,m), &
322  qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m))
323  CALL moments_twenty_eight_nodes (qmomk_n1(:,ijk,m), &
324  qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m), qmomk_m1(:,ijk,m))
325  END IF
326  END DO
327  END DO
328 
329  ! RHS of moments equations
330  ! Gravity and drag contribution (Half time step)
331  DO m = 1, mmax
332  DO ijk = ijkstart3, ijkend3
333  IF (fluid_at(ijk)) THEN
334  i = i_of(ijk)
335  j = j_of(ijk)
336  k = k_of(ijk)
337  imjk = im_of(ijk)
338  ijmk = jm_of(ijk)
339  ijkm = km_of(ijk)
340 
341  ugc = avg_x_e(u_g(imjk),u_g(ijk),i)
342  vgc = avg_y_n(v_g(ijmk),v_g(ijk))
343  wgc = avg_z_t(w_g(ijkm),w_g(ijk))
344 
345  DO in = 1, qmomk_nn
346  vrel = sqrt(((qmomk_u1(in,ijk,m) - ugc)**2) + ((qmomk_v1(in,ijk,m) - vgc)**2) + ((qmomk_w1(in,ijk,m) - wgc)**2))
347 
348  ! Particle Reynolds number
349  rep = ro_g0*d_p0(m)*vrel/mu_g0
350  cd = 24.d0*(ep_g(ijk)**(-2.65))*(1.d0 + 0.15d0*((ep_g(ijk)*rep)**0.687))/(rep + small_number)
351  beta_drag = 3.d0*ro_g0*cd*vrel/(4.d0*d_p0(m)*ro_s(ijk,m))
352 
353  drag_exp = exp(-0.5*qmomk_dt*beta_drag)
354 
355  IF (beta_drag > small_number) THEN !CHECK HERE, it was /= 0
356  qmomk_tau_drag(in,ijk,m) = 1.d0/beta_drag
357  ELSE
358  qmomk_tau_drag(in,ijk,m) = large_number
359  END IF
360 
361  ! Drag calculation for QMOMK
362  qmomk_u1(in,ijk,m) = drag_exp*qmomk_u1(in,ijk,m) + (1.d0 - drag_exp)*ugc
363  qmomk_v1(in,ijk,m) = drag_exp*qmomk_v1(in,ijk,m) + (1.d0 - drag_exp)*(vgc - gravity*qmomk_tau_drag(in,ijk,m))
364  qmomk_w1(in,ijk,m) = drag_exp*qmomk_w1(in,ijk,m) + (1.d0 - drag_exp)*wgc
365  END DO
366  CALL moments_twenty_eight_nodes (qmomk_n1(:,ijk,m), &
367  qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m), qmomk_m1(:,ijk,m))
368  CALL eight_node_3d (qmomk_m1(:,ijk,m), qmomk_n1(:,ijk,m), &
369  qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m))
370  CALL moments_twenty_eight_nodes (qmomk_n1(:,ijk,m), &
371  qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m), qmomk_m1(:,ijk,m))
372  END IF
373  END DO
374  END DO
375 
376  DO m = 1, mmax
377  DO ijk = ijkstart3, ijkend3
378  IF (fluid_at(ijk)) THEN
379  CALL bind_theta(qmomk_m1(:,ijk,m), minval(qmomk_tau_drag(:,ijk,m)))
380  CALL eight_node_3d (qmomk_m1(:,ijk,m), qmomk_n1(:,ijk,m), &
381  qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m))
382  CALL moments_twenty_eight_nodes (qmomk_n1(:,ijk,m), &
383  qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m), qmomk_m1(:,ijk,m))
384  END IF
385  END DO
386  END DO
387 
388  ! End of gravity and drag contribution
389 
390  print *,'Before collisions'
391 
392  ! Collisions contribution (Half time step)
393  IF (qmomk_collisions == 'BGK' .AND. mmax == 1) THEN
394  DO m = 1, mmax
395  DO ijk = ijkstart3, ijkend3
396  IF (fluid_at(ijk)) THEN
397  mom0 = qmomk_m1(1,ijk,m)
398  IF (mom0 > epsn) THEN
399  CALL collisions_bgk(qmomk_m1(:,ijk,m), 0.5*qmomk_dt, qmomk_tcol, d_p0(m), c_e)
400  CALL eight_node_3d (qmomk_m1(:,ijk,m), qmomk_n1(:,ijk,m), &
401  qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m))
402  CALL moments_twenty_eight_nodes (qmomk_n1(:,ijk,m), &
403  qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m), qmomk_m1(:,ijk,m))
404  END IF
405  END IF
406  END DO
407  END DO
408  ELSE IF (qmomk_collisions == 'BOLTZMANN' .AND. mmax == 1) THEN
409  DO ijk = ijkstart3, ijkend3
410  IF (fluid_at(ijk)) THEN
411  CALL solve_boltzmann_collisions_one_specie(qmomk_m1(:,ijk,1), qmomk_n1(:,ijk,1), &
412  qmomk_u1(:,ijk,1), qmomk_v1(:,ijk,1), qmomk_w1(:,ijk,1), &
413  0.5*qmomk_dt, c_e, d_p0(1), qmomk_collisions_order)
414  CALL eight_node_3d (qmomk_m1(:,ijk,1), qmomk_n1(:,ijk,1), qmomk_u1(:,ijk,1), qmomk_v1(:,ijk,1), qmomk_w1(:,ijk,1))
415  CALL moments_twenty_eight_nodes (qmomk_n1(:,ijk,1), &
416  qmomk_u1(:,ijk,1), qmomk_v1(:,ijk,1), qmomk_w1(:,ijk,1), qmomk_m1(:,ijk,1))
417  END IF
418  END DO
419  ELSE
420  DO m = 1, mmax
421  DO ijk = ijkstart3, ijkend3
422  IF (fluid_at(ijk)) THEN
423  DO m2 = 1, mmax
424  CALL solve_boltzmann_collisions_two_species(qmomk_m1(:,ijk,m), qmomk_n1(:,ijk,m), &
425  qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m), &
426  qmomk_m1(:,ijk,m2), qmomk_n1(:,ijk,m2), qmomk_u1(:,ijk,m2), qmomk_v1(:,ijk,m2), qmomk_w1(:,ijk,m2), &
427  0.5*qmomk_dt, 1.d0/6.d0*pi*(d_p0(m)**3)*ro_s(ijk,m), 1.d0/6.d0*pi*(d_p0(m2)**3)*ro_s(ijk,m2), &
428  d_p0(m), d_p0(m2), c_e, c_e, qmomk_collisions_order)
429  END DO
430  CALL eight_node_3d (qmomk_m1(:,ijk,m), qmomk_n1(:,ijk,m), qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m))
431  CALL moments_twenty_eight_nodes (qmomk_n1(:,ijk,m), &
432  qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m), qmomk_m1(:,ijk,m))
433  END IF
434  END DO
435  END DO
436  END IF
437 
438  print *,'After collisions'
439 
440  ! End of collisions contribution
441 
442  ! Boundary conditions update
443  CALL qmomk_set_bc
444 
445  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
446  ! Second step of Runge-Kutta time integration !
447  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
448 
449  ! M0 -> old values
450  ! M1 -> first step RK values
451 
452  ! Calculating kinetic based fluxes (Full time step)
453  DO m = 1, mmax
454  DO ijk = ijkstart3, ijkend3
455  IF (fluid_at(ijk)) THEN
456  i = i_of(ijk)
457  j = j_of(ijk)
458  k = k_of(ijk)
459 
460  imjk = im_of(ijk)
461  ijmk = jm_of(ijk)
462  ijkm = km_of(ijk)
463 
464  ipjk = ip_of(ijk)
465  ijpk = jp_of(ijk)
466  ijkp = kp_of(ijk)
467 
468  ! 2D Case
469 
470  IF (no_k) THEN
471  ! x - direction
472  nlminus = qmomk_n1(:,imjk,m)
473  ulminus = qmomk_u1(:,imjk,m)
474  vlminus = qmomk_v1(:,imjk,m)
475  wlminus = qmomk_w1(:,imjk,m)
476 
477  nlplus = qmomk_n1(:,ijk,m)
478  ulplus = qmomk_u1(:,ijk,m)
479  vlplus = qmomk_v1(:,ijk,m)
480  wlplus = qmomk_w1(:,ijk,m)
481 
482  nrminus = qmomk_n1(:,ijk,m)
483  urminus = qmomk_u1(:,ijk,m)
484  vrminus = qmomk_v1(:,ijk,m)
485  wrminus = qmomk_w1(:,ijk,m)
486 
487  nrplus = qmomk_n1(:,ipjk,m)
488  urplus = qmomk_u1(:,ipjk,m)
489  vrplus = qmomk_v1(:,ipjk,m)
490  wrplus = qmomk_w1(:,ipjk,m)
491 
492  CALL kinetic_flux_x_twenty_eight_nodes (nlminus, ulminus, vlminus, wlminus, &
493  nlplus, ulplus, vlplus, wlplus, f_x_left)
494  CALL kinetic_flux_x_twenty_eight_nodes (nrminus, urminus, vrminus, wrminus, &
495  nrplus, urplus, vrplus, wrplus, f_x_right)
496 
497  ! y - direction
498  nlminus = qmomk_n1(:,ijmk,m)
499  ulminus = qmomk_u1(:,ijmk,m)
500  vlminus = qmomk_v1(:,ijmk,m)
501  wlminus = qmomk_w1(:,ijmk,m)
502 
503  nlplus = qmomk_n1(:,ijk,m)
504  ulplus = qmomk_u1(:,ijk,m)
505  vlplus = qmomk_v1(:,ijk,m)
506  wlplus = qmomk_w1(:,ijk,m)
507 
508  nrminus = qmomk_n1(:,ijk,m)
509  urminus = qmomk_u1(:,ijk,m)
510  vrminus = qmomk_v1(:,ijk,m)
511  wrminus = qmomk_w1(:,ijk,m)
512 
513  nrplus = qmomk_n1(:,ijpk,m)
514  urplus = qmomk_u1(:,ijpk,m)
515  vrplus = qmomk_v1(:,ijpk,m)
516  wrplus = qmomk_w1(:,ijpk,m)
517 
518  CALL kinetic_flux_y_twenty_eight_nodes (nlminus, ulminus, vlminus, wlminus, &
519  nlplus, ulplus, vlplus, wlplus, f_y_left)
520  CALL kinetic_flux_y_twenty_eight_nodes (nrminus, urminus, vrminus, wrminus, &
521  nrplus, urplus, vrplus, wrplus, f_y_right)
522 
523  qmomk_m1(:,ijk,m) = qmomk_m0(:,ijk,m) - (qmomk_dt/minval(dx))*(f_x_right - f_x_left) &
524  - (qmomk_dt/minval(dy))*(f_y_right - f_y_left)
525 
526  ! 3D case
527  ELSE
528  ! x - direction
529  nlminus = qmomk_n1(:,imjk,m)
530  ulminus = qmomk_u1(:,imjk,m)
531  vlminus = qmomk_v1(:,imjk,m)
532  wlminus = qmomk_w1(:,imjk,m)
533 
534  nlplus = qmomk_n1(:,ijk,m)
535  ulplus = qmomk_u1(:,ijk,m)
536  vlplus = qmomk_v1(:,ijk,m)
537  wlplus = qmomk_w1(:,ijk,m)
538 
539  nrminus = qmomk_n1(:,ijk,m)
540  urminus = qmomk_u1(:,ijk,m)
541  vrminus = qmomk_v1(:,ijk,m)
542  wrminus = qmomk_w1(:,ijk,m)
543 
544  nrplus = qmomk_n1(:,ipjk,m)
545  urplus = qmomk_u1(:,ipjk,m)
546  vrplus = qmomk_v1(:,ipjk,m)
547  wrplus = qmomk_w1(:,ipjk,m)
548 
549  CALL kinetic_flux_x_twenty_eight_nodes (nlminus, ulminus, vlminus, wlminus, &
550  nlplus, ulplus, vlplus, wlplus, f_x_left)
551  CALL kinetic_flux_x_twenty_eight_nodes (nrminus, urminus, vrminus, wrminus, &
552  nrplus, urplus, vrplus, wrplus, f_x_right)
553 
554  ! y - direction
555  nlminus = qmomk_n1(:,ijmk,m)
556  ulminus = qmomk_u1(:,ijmk,m)
557  vlminus = qmomk_v1(:,ijmk,m)
558  wlminus = qmomk_w1(:,ijmk,m)
559 
560  nlplus = qmomk_n1(:,ijk,m)
561  ulplus = qmomk_u1(:,ijk,m)
562  vlplus = qmomk_v1(:,ijk,m)
563  wlplus = qmomk_w1(:,ijk,m)
564 
565  nrminus = qmomk_n1(:,ijk,m)
566  urminus = qmomk_u1(:,ijk,m)
567  vrminus = qmomk_v1(:,ijk,m)
568  wrminus = qmomk_w1(:,ijk,m)
569 
570  nrplus = qmomk_n1(:,ijpk,m)
571  urplus = qmomk_u1(:,ijpk,m)
572  vrplus = qmomk_v1(:,ijpk,m)
573  wrplus = qmomk_w1(:,ijpk,m)
574 
575  CALL kinetic_flux_y_twenty_eight_nodes (nlminus, ulminus, vlminus, wlminus, &
576  nlplus, ulplus, vlplus, wlplus, f_y_left)
577  CALL kinetic_flux_y_twenty_eight_nodes (nrminus, urminus, vrminus, wrminus, &
578  nrplus, urplus, vrplus, wrplus, f_y_right)
579 
580  ! z - direction
581  nlminus = qmomk_n1(:,ijkm,m)
582  ulminus = qmomk_u1(:,ijkm,m)
583  vlminus = qmomk_v1(:,ijkm,m)
584  wlminus = qmomk_w1(:,ijkm,m)
585 
586  nlplus = qmomk_n1(:,ijk,m)
587  ulplus = qmomk_u1(:,ijk,m)
588  vlplus = qmomk_v1(:,ijk,m)
589  wlplus = qmomk_w1(:,ijk,m)
590 
591  nrminus = qmomk_n1(:,ijk,m)
592  urminus = qmomk_u1(:,ijk,m)
593  vrminus = qmomk_v1(:,ijk,m)
594  wrminus = qmomk_w1(:,ijk,m)
595 
596  nrplus = qmomk_n1(:,ijkp,m)
597  urplus = qmomk_u1(:,ijkp,m)
598  vrplus = qmomk_v1(:,ijkp,m)
599  wrplus = qmomk_w1(:,ijkp,m)
600 
601  CALL kinetic_flux_z_twenty_eight_nodes (nlminus, ulminus, vlminus, wlminus, &
602  nlplus, ulplus, vlplus, wlplus, f_z_left)
603  CALL kinetic_flux_z_twenty_eight_nodes (nrminus, urminus, vrminus, wrminus, &
604  nrplus, urplus, vrplus, wrplus, f_z_right)
605 
606  qmomk_m1(:,ijk,m) = qmomk_m0(:,ijk,m) - (qmomk_dt/minval(dx))*(f_x_right - f_x_left) - &
607  (qmomk_dt/minval(dy))*(f_y_right - f_y_left) - (qmomk_dt/minval(dz))*(f_z_right - f_z_left)
608  END IF
609  END IF
610  END DO
611  END DO
612 
613  ! Update weights and abscissas
614  DO m = 1, mmax
615  DO ijk = ijkstart3, ijkend3
616  IF (fluid_at(ijk)) THEN
617  CALL eight_node_3d (qmomk_m1(:,ijk,m), qmomk_n1(:,ijk,m), qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m))
618  CALL moments_twenty_eight_nodes (qmomk_n1(:,ijk,m), &
619  qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m), qmomk_m1(:,ijk,m))
620  END IF
621  END DO
622  END DO
623 
624  ! RHS of moments equations
625  ! Collisions contribution (Full time step)
626  IF (qmomk_collisions == 'BGK' .AND. mmax == 1) THEN
627  DO m = 1, mmax
628  DO ijk = ijkstart3, ijkend3
629  IF (fluid_at(ijk)) THEN
630  mom0 = qmomk_m1(1,ijk,m)
631  IF (mom0 > epsn) THEN
632  CALL collisions_bgk(qmomk_m1(:,ijk,m), qmomk_dt, qmomk_tcol, d_p0(m), c_e)
633  CALL eight_node_3d (qmomk_m1(:,ijk,m), qmomk_n1(:,ijk,m), &
634  qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m))
635  CALL moments_twenty_eight_nodes (qmomk_n1(:,ijk,m), &
636  qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m), qmomk_m1(:,ijk,m))
637  END IF
638  END IF
639  END DO
640  END DO
641  ELSE IF (qmomk_collisions == 'BOLTZMANN' .AND. mmax == 1) THEN
642  DO ijk = ijkstart3, ijkend3
643  IF (fluid_at(ijk)) THEN
644  CALL solve_boltzmann_collisions_one_specie(qmomk_m1(:,ijk,1), qmomk_n1(:,ijk,1), &
645  qmomk_u1(:,ijk,1), qmomk_v1(:,ijk,1), qmomk_w1(:,ijk,1), &
646  qmomk_dt, c_e, d_p0(1), qmomk_collisions_order)
647  CALL eight_node_3d (qmomk_m1(:,ijk,1), qmomk_n1(:,ijk,1), qmomk_u1(:,ijk,1), qmomk_v1(:,ijk,1), qmomk_w1(:,ijk,1))
648  CALL moments_twenty_eight_nodes (qmomk_n1(:,ijk,1), &
649  qmomk_u1(:,ijk,1), qmomk_v1(:,ijk,1), qmomk_w1(:,ijk,1), qmomk_m1(:,ijk,1))
650  END IF
651  END DO
652  ELSE
653  DO m = 1, mmax
654  DO ijk = ijkstart3, ijkend3
655  IF (fluid_at(ijk)) THEN
656  DO m2 = 1, mmax
657  CALL solve_boltzmann_collisions_two_species(qmomk_m1(:,ijk,m), qmomk_n1(:,ijk,m), &
658  qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m), &
659  qmomk_m1(:,ijk,m2), qmomk_n1(:,ijk,m2), &
660  qmomk_u1(:,ijk,m2), qmomk_v1(:,ijk,m2), qmomk_w1(:,ijk,m2), &
661  qmomk_dt, 1.d0/6.d0*pi*(d_p0(m)**3)*ro_s(ijk,m), &
662  1.d0/6.d0*pi*(d_p0(m2)**3)*ro_s(ijk,m2), d_p0(m), d_p0(m2), &
663  c_e, c_e, qmomk_collisions_order)
664  END DO
665  CALL eight_node_3d (qmomk_m1(:,ijk,m), qmomk_n1(:,ijk,m), qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m))
666  CALL moments_twenty_eight_nodes (qmomk_n1(:,ijk,m), &
667  qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m), qmomk_m1(:,ijk,m))
668  END IF
669  END DO
670  END DO
671  END IF
672  ! End of collisions contribution
673 
674  ! Gravity and drag contribution (Full time step)
675  DO m = 1, mmax
676  DO ijk = ijkstart3, ijkend3
677  IF (fluid_at(ijk)) THEN
678  i = i_of(ijk)
679 
680  imjk = im_of(ijk)
681  ijmk = jm_of(ijk)
682  ijkm = km_of(ijk)
683 
684  ugc = avg_x_e(u_g(imjk),u_g(ijk),i)
685  vgc = avg_y_n(v_g(ijmk),v_g(ijk))
686  wgc = avg_z_t(w_g(ijkm),w_g(ijk))
687 
688  DO in = 1, qmomk_nn
689  vrel = sqrt(((qmomk_u1(in,ijk,m) - ugc)**2) + ((qmomk_v1(in,ijk,m) - vgc)**2) + ((qmomk_w1(in,ijk,m) - wgc)**2))
690 
691  ! Particle Reynolds number
692  rep = ro_g0*d_p0(m)*vrel/mu_g0
693  cd = 24.d0*(ep_g(ijk)**(-2.65))*(1.d0 + 0.15*((ep_g(ijk)*rep)**0.687))/(rep+small_number)
694  beta_drag = 3.d0*ro_g0*cd*vrel/(4.d0*d_p0(m)*ro_s(ijk,m))
695 
696  IF (beta_drag > small_number) THEN ! CHECK HERE, it was /= 0.
697  qmomk_tau_drag(in,ijk,m) = 1.d0/beta_drag
698  ELSE
699  qmomk_tau_drag(in,ijk,m) = large_number
700  END IF
701 
702  drag_exp = exp(-qmomk_dt*beta_drag)
703 
704  ! Drag calculation for QMOMK
705  qmomk_u1(in,ijk,m) = drag_exp*qmomk_u1(in,ijk,m) + (1.d0 - drag_exp)*ugc
706  qmomk_v1(in,ijk,m) = drag_exp*qmomk_v1(in,ijk,m) + (1.d0 - drag_exp)*(vgc - gravity*qmomk_tau_drag(in,ijk,m))
707  qmomk_w1(in,ijk,m) = drag_exp*qmomk_w1(in,ijk,m) + (1.d0 - drag_exp)*wgc
708  END DO
709  CALL moments_twenty_eight_nodes (qmomk_n1(:,ijk,m), &
710  qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m), qmomk_m1(:,ijk,m))
711  CALL eight_node_3d (qmomk_m1(:,ijk,m), qmomk_n1(:,ijk,m), qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m))
712  CALL moments_twenty_eight_nodes (qmomk_n1(:,ijk,m), &
713  qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m), qmomk_m1(:,ijk,m))
714  END IF
715  END DO
716  END DO
717  ! End of gravity and drag contribution
718 
719  DO m = 1, mmax
720  DO ijk = ijkstart3, ijkend3
721  IF (fluid_at(ijk)) THEN
722  CALL bind_theta(qmomk_m1(:,ijk,m), minval(qmomk_tau_drag(:,ijk,m)))
723  CALL eight_node_3d (qmomk_m1(:,ijk,m), qmomk_n1(:,ijk,m), qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m))
724  CALL moments_twenty_eight_nodes (qmomk_n1(:,ijk,m), &
725  qmomk_u1(:,ijk,m), qmomk_v1(:,ijk,m), qmomk_w1(:,ijk,m), qmomk_m1(:,ijk,m))
726  END IF
727  END DO
728  END DO
729 
730  ! Boundary conditions update
731  CALL qmomk_set_bc
732 
733  ! Preparing data for next time step
734  qmomk_m0 = qmomk_m1
735  qmomk_n0 = qmomk_n1
736  qmomk_u0 = qmomk_u1
737  qmomk_v0 = qmomk_v1
738  qmomk_w0 = qmomk_w1
739 
740  print *,'QMOM DT 5 = ',qmomk_dt
741 
742  qmomk_time = qmomk_time + qmomk_dt
743 
744  print *,'QMOM DT 6 = ',qmomk_dt
745  END DO ! QMOMK internal time stepping
746 
747  IF (dt .LT. qmomk_dt_tmp) THEN
748  qmomk_dt = qmomk_dt_tmp
749  ENDIF
750 
751  IF (tmp_dts .NE. zero) THEN
752  qmomk_dt = tmp_dts
753  tmp_dts = zero
754  ENDIF
755 
756  print *,'Time = ',time
757  print *,'Time QMOMK = ',qmomk_time
758 
759  ! Resetting gas volume fraction to one, in order to update it
760  DO ijk = ijkstart3, ijkend3
761  IF (fluid_at(ijk)) THEN
762  ep_g(ijk) = one
763  END IF
764  END DO
765 
766  ! Passing values to the fluid solver and for storage
767  DO m = 1, mmax
768  DO ijk = ijkstart3, ijkend3
769  IF (fluid_at(ijk)) THEN
770  ! Mean particle velocities
771  u_s(ijk, m) = qmomk_m1(2,ijk,m)/qmomk_m1(1,ijk,m)
772  v_s(ijk, m) = qmomk_m1(3,ijk,m)/qmomk_m1(1,ijk,m)
773  w_s(ijk, m) = qmomk_m1(4,ijk,m)/qmomk_m1(1,ijk,m)
774 
775  ! Gas-phase volume fraction
776  ep_g(ijk) = ep_g(ijk) - qmomk_m1(1,ijk,m)
777 
778  ! Particle phases "densities"
779  rop_s(ijk, m) = qmomk_m1(1,ijk,m)*ro_s(ijk,m)
780 
781  ! Particle phases granular temperatures
782  theta_m(ijk,m) = ((qmomk_m1(5,ijk,m)/qmomk_m1(1,ijk,m) - &
783  (qmomk_m1(2,ijk,m)/qmomk_m1(1,ijk,m))*(qmomk_m1(2,ijk,m)/qmomk_m1(1,ijk,m))) + &
784  (qmomk_m1(8,ijk,m)/qmomk_m1(1,ijk,m) - &
785  (qmomk_m1(3,ijk,m)/qmomk_m1(1,ijk,m))*(qmomk_m1(3,ijk,m)/qmomk_m1(1,ijk,m))) + &
786  (qmomk_m1(10,ijk,m)/qmomk_m1(1,ijk,m) - &
787  (qmomk_m1(4,ijk,m)/qmomk_m1(1,ijk,m))*(qmomk_m1(4,ijk,m)/qmomk_m1(1,ijk,m))))/3.d0
788 
789  END IF
790  END DO
791  END DO
792 
793  ! Calculating momentum exchange terms
794  IF (qmomk_coupled) THEN
795  DO m = 1, mmax
796  DO ijk = ijkstart3, ijkend3
797  IF (fluid_at(ijk)) THEN
798  i = i_of(ijk)
799 
800  imjk = im_of(ijk)
801  ijmk = jm_of(ijk)
802  ijkm = km_of(ijk)
803 
804  ugc = avg_x_e(u_g(imjk),u_g(ijk),i)
805  vgc = avg_y_n(v_g(ijmk),v_g(ijk))
806  wgc = avg_z_t(w_g(ijkm),w_g(ijk))
807 
808  DO in = 1, qmomk_nn
809  vrel = sqrt(((qmomk_u1(in,ijk,m) - ugc)**2) + ((qmomk_v1(in,ijk,m) - vgc)**2) + ((qmomk_w1(in,ijk,m) - wgc)**2))
810 
811  rep = ro_g0*d_p0(m)*vrel/mu_g0
812  cd = 24.d0*(ep_g(ijk)**(-2.65))*(1.d0 + 0.15*((ep_g(ijk)*rep)**0.687))/(rep+small_number)
813  beta_drag = 3.d0*qmomk_n1(in,ijk,m)*ro_g0*cd*vrel/(4.d0*d_p0(m))
814 
815  qmomk_f_gs(in,ijk,m) = beta_drag
816  END DO
817  END IF
818  END DO
819  END DO
820  END IF
821 END SUBROUTINE qmomk_time_march
qmomk_time_march
subroutine qmomk_time_march
Definition: qmomk_time_march.f:12
constant::c_e
double precision c_e
Definition: constant_mod.f:105
fldvar::v_s
double precision, dimension(:,:), allocatable v_s
Definition: fldvar_mod.f:105
qmomk_quadrature::moments_twenty_eight_nodes
subroutine, public moments_twenty_eight_nodes(n, u, v, w, mom)
Definition: qmomk_quadrature_mod.f:91
qmomk_initial_conditions
subroutine qmomk_initial_conditions
Definition: qmomk_initial_conditions.f:10
indices::i_of
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
visc_s
Definition: visc_s_mod.f:1
physprop::d_p0
double precision, dimension(dim_m) d_p0
Definition: physprop_mod.f:25
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
fldvar::ep_g
double precision, dimension(:), allocatable ep_g
Definition: fldvar_mod.f:15
param1::one
double precision, parameter one
Definition: param1_mod.f:29
qmomk_parameters::epsn
double precision, parameter epsn
Definition: qmomk_parameters_mod.f:21
fldvar::w_s
double precision, dimension(:,:), allocatable w_s
Definition: fldvar_mod.f:117
rxns
Definition: rxns_mod.f:1
qmom_kinetic_equation::qmomk_w1
double precision, dimension(:,:,:), allocatable qmomk_w1
Definition: qmom_kinetic_equation_mod.f:33
qmomk_collision::solve_boltzmann_collisions_one_specie
subroutine, public solve_boltzmann_collisions_one_specie(M, N, U, V, W, dt, e, dp, order)
Definition: qmomk_collision_mod.f:692
qmom_kinetic_equation::qmomk_collisions_order
integer qmomk_collisions_order
Definition: qmom_kinetic_equation_mod.f:65
physprop::mu_g0
double precision mu_g0
Definition: physprop_mod.f:62
drag
Definition: drag_mod.f:11
geometry::dy
double precision, dimension(0:dim_j) dy
Definition: geometry_mod.f:70
qmomk_parameters
Definition: qmomk_parameters_mod.f:1
tau_s
Definition: tau_s_mod.f:1
run::dt
double precision dt
Definition: run_mod.f:51
qmom_kinetic_equation::qmomk_v0
double precision, dimension(:,:,:), allocatable qmomk_v0
Definition: qmom_kinetic_equation_mod.f:24
qmom_kinetic_equation::qmomk_w0
double precision, dimension(:,:,:), allocatable qmomk_w0
Definition: qmom_kinetic_equation_mod.f:25
qmomk_fluxes::kinetic_flux_y_twenty_eight_nodes
subroutine, public kinetic_flux_y_twenty_eight_nodes(Nl, Ul, Vl, Wl, Nr, Ur, Vr, Wr, f)
Definition: qmomk_fluxes_mod.f:123
geometry::dz
double precision, dimension(0:dim_k) dz
Definition: geometry_mod.f:72
is
Definition: is_mod.f:11
fldvar::u_s
double precision, dimension(:,:), allocatable u_s
Definition: fldvar_mod.f:93
qmomk_collision::collisions_bgk
subroutine, public collisions_bgk(M, Dt, tcol, dp, e)
Definition: qmomk_collision_mod.f:121
qmom_kinetic_equation::qmomk_collisions
character(64) qmomk_collisions
Definition: qmom_kinetic_equation_mod.f:58
qmomk_set_bc
subroutine qmomk_set_bc
Definition: qmomk_set_bc.f:10
indices::k_of
integer, dimension(:), allocatable k_of
Definition: indices_mod.f:47
fun_avg
Definition: fun_avg_mod.f:1
physprop::mmax
integer mmax
Definition: physprop_mod.f:19
qmom_kinetic_equation::qmomk_n1
double precision, dimension(:,:,:), allocatable qmomk_n1
Definition: qmom_kinetic_equation_mod.f:27
qmomk_quadrature::eight_node_3d
subroutine, public eight_node_3d(mom, n, u, v, w)
Definition: qmomk_quadrature_mod.f:220
indices::j_of
integer, dimension(:), allocatable j_of
Definition: indices_mod.f:46
tau_g
Definition: tau_g_mod.f:1
physprop::ro_g0
double precision ro_g0
Definition: physprop_mod.f:59
qmom_kinetic_equation::qmomk_collision_time
double precision, dimension(:,:), allocatable qmomk_collision_time
Definition: qmom_kinetic_equation_mod.f:44
param1::small_number
double precision, parameter small_number
Definition: param1_mod.f:24
qmomk_collision::compute_collision_time
subroutine, public compute_collision_time(M, dp, theta, tcol, dt)
Definition: qmomk_collision_mod.f:246
ur_facs
Definition: ur_facs_mod.f:1
qmom_kinetic_equation::qmomk_n0
double precision, dimension(:,:,:), allocatable qmomk_n0
Definition: qmom_kinetic_equation_mod.f:19
run::run_type
character(len=16) run_type
Definition: run_mod.f:33
qmom_kinetic_equation::qmomk_u1
double precision, dimension(:,:,:), allocatable qmomk_u1
Definition: qmom_kinetic_equation_mod.f:31
fldvar::theta_m
double precision, dimension(:,:), allocatable theta_m
Definition: fldvar_mod.f:149
qmom_kinetic_equation::qmomk_tau_drag
double precision, dimension(:,:,:), allocatable qmomk_tau_drag
Definition: qmom_kinetic_equation_mod.f:40
time_cpu
Definition: time_cpu_mod.f:1
qmomk_fluxes
Definition: qmomk_fluxes_mod.f:23
fldvar::v_g
double precision, dimension(:), allocatable v_g
Definition: fldvar_mod.f:99
geometry::dx
double precision, dimension(0:dim_i) dx
Definition: geometry_mod.f:68
fldvar::w_g
double precision, dimension(:), allocatable w_g
Definition: fldvar_mod.f:111
run
Definition: run_mod.f:13
qmomk_fluxes::kinetic_flux_x_twenty_eight_nodes
subroutine, public kinetic_flux_x_twenty_eight_nodes(Nl, Ul, Vl, Wl, Nr, Ur, Vr, Wr, f)
Definition: qmomk_fluxes_mod.f:41
param1::large_number
double precision, parameter large_number
Definition: param1_mod.f:23
param
Definition: param_mod.f:2
fldvar::ro_s
double precision, dimension(:,:), allocatable ro_s
Definition: fldvar_mod.f:45
qmom_kinetic_equation::qmomk_coupled
logical qmomk_coupled
Definition: qmom_kinetic_equation_mod.f:49
geometry::no_k
logical no_k
Definition: geometry_mod.f:28
qmom_kinetic_equation
Definition: qmom_kinetic_equation_mod.f:12
constant::gravity
double precision gravity
Definition: constant_mod.f:149
compar::ijkstart3
integer ijkstart3
Definition: compar_mod.f:80
qmomk_collision
Definition: qmomk_collision_mod.f:19
visc_g
Definition: visc_g_mod.f:1
fldvar::u_g
double precision, dimension(:), allocatable u_g
Definition: fldvar_mod.f:87
qmomk_quadrature
Definition: qmomk_quadrature_mod.f:16
fldvar::rop_s
double precision, dimension(:,:), allocatable rop_s
Definition: fldvar_mod.f:51
qmom_kinetic_equation::qmomk_v1
double precision, dimension(:,:,:), allocatable qmomk_v1
Definition: qmom_kinetic_equation_mod.f:32
qmomk_init_bc
subroutine qmomk_init_bc
Definition: qmomk_init_bc.f:10
qmomk_quadrature::bind_theta
subroutine, public bind_theta(mom, tau_p)
Definition: qmomk_quadrature_mod.f:552
qmomk_fluxes::kinetic_flux_z_twenty_eight_nodes
subroutine, public kinetic_flux_z_twenty_eight_nodes(Nl, Ul, Vl, Wl, Nr, Ur, Vr, Wr, f)
Definition: qmomk_fluxes_mod.f:206
constant::pi
double precision, parameter pi
Definition: constant_mod.f:158
run::time
double precision time
Definition: run_mod.f:45
qmom_kinetic_equation::qmomk_dt
double precision qmomk_dt
Definition: qmom_kinetic_equation_mod.f:68
qmom_kinetic_equation::qmomk_cfl
double precision qmomk_cfl
Definition: qmom_kinetic_equation_mod.f:71
qmom_kinetic_equation::qmomk_m0
double precision, dimension(:,:,:), allocatable qmomk_m0
Definition: qmom_kinetic_equation_mod.f:21
qmomk_parameters::qmomk_nn
integer, parameter qmomk_nn
Definition: qmomk_parameters_mod.f:11
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
qmom_kinetic_equation::qmomk_f_gs
double precision, dimension(:,:,:), allocatable qmomk_f_gs
Definition: qmom_kinetic_equation_mod.f:41
qmomk_collision::solve_boltzmann_collisions_two_species
subroutine, public solve_boltzmann_collisions_two_species(M1, N1, U1, V1, W1, M2, N2, U2, V2, W2, dt, m_1, m_2, dp1, dp2, e11, e12, order)
Definition: qmomk_collision_mod.f:733
qmom_kinetic_equation::qmomk_m1
double precision, dimension(:,:,:), allocatable qmomk_m1
Definition: qmom_kinetic_equation_mod.f:29
qmom_kinetic_equation::qmomk_u0
double precision, dimension(:,:,:), allocatable qmomk_u0
Definition: qmom_kinetic_equation_mod.f:23