MFIX  2016-1
qmomk_quadrature_mod.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Module: qmomk_quadrature C
4 ! Purpose: Routines to perform the quadrature calculations C
5 ! Contains the following subroutines: C
6 ! quadrature_bounded, moments_twenty_eight_nodes, C
7 ! check_moments_twenty, eight_node_3d, bind_theta C
8 ! C
9 ! Author: A. Passalacqua Date: 18-Jun-2008 C
10 ! Reviewer: Date: dd-mmm-yy C
11 ! C
12 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
13 
14 #include "version.inc"
15 
16 MODULE qmomk_quadrature
17 
18  USE qmomk_tools
19  USE qmomk_parameters
20 
21  IMPLICIT NONE
22 
23  PRIVATE
24 
25  PUBLIC :: moments_twenty_eight_nodes
26  PUBLIC :: eight_node_3d
27  PUBLIC :: bind_theta
28 
29 CONTAINS
30 
31 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
32 ! C
33 ! Two-nodes quadrature C
34 ! C
35 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
36 
37  SUBROUTINE quadrature_bounded (m, w, u)
38 
39  IMPLICIT NONE
40 
41  DOUBLE PRECISION, INTENT(IN), DIMENSION(4) :: m
42  DOUBLE PRECISION, INTENT(OUT), DIMENSION(2) :: w
43  DOUBLE PRECISION, INTENT(OUT), DIMENSION(2) :: u
44 
45  DOUBLE PRECISION :: eps = 1.d-8
46  DOUBLE PRECISION :: smax = 1.d6
47 
48  DOUBLE PRECISION :: rho = 0.d0
49  DOUBLE PRECISION :: um, sigma, x, q, s
50 
51  rho = m(1)
52  um = m(2)/rho
53 
54  sigma = (m(3)*rho - m(2)*m(2))/(rho**2)
55 
56  IF (sigma <= 0.d0) THEN
57  print *,'QMOMK: Negative variance in quadrature'
58  sigma = 0.d0
59  x = 0
60  ELSE
61  sigma = sqrt(sigma)
62  q = m(4)/rho - um**3 - 3.d0*um*(sigma**2)
63  s = min(smax, abs(q)/(sigma**3))
64  x = 0.d0
65  IF ((s**2) < eps) THEN
66  x = 0.d0
67  ELSE
68  x = 0.5 * sign(1.d0, q)/sqrt(1.d0 + 4.d0/(s**2))
69  END IF
70  END IF
71 
72  w(1) = rho*(0.5 + x)
73  u(1) = um - sqrt((0.5-x)/(0.5+x))*sigma
74 
75  w(2) = rho*(0.5-x)
76  u(2) = um + sqrt((0.5+x)/(0.5-x))*sigma
77 
78  IF (abs(x) >= 0.5d0) THEN
79  print *,'QMOMK: Warning: x > 0.5'
80  print *,x
81  END IF
82  END SUBROUTINE quadrature_bounded
83 
84 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
85 ! C
86 ! Calculation of the twenty moments from weights and abscissas C
87 ! C
88 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
89 
90  SUBROUTINE moments_twenty_eight_nodes (n, u, v, w, mom)
91 
92  USE param1, only: small_number
93  IMPLICIT NONE
94 
95  INTEGER :: i
96 
97  DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: n
98  DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: u
99  DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: v
100  DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: w
101  DOUBLE PRECISION, INTENT(OUT), DIMENSION(QMOMK_NMOM) :: mom
102 
103  DOUBLE PRECISION, dimension(QMOMK_NMOM) :: M
104 
105  m=0.d0
106 
107  DO i = 1, qmomk_nn
108  m(1) = m(1) + n(i)
109  m(2) = m(2) + n(i)*u(i)
110  m(3) = m(3) + n(i)*v(i)
111  m(4) = m(4) + n(i)*w(i)
112  m(5) = m(5) + n(i)*(u(i)**2)
113  m(6) = m(6) + n(i)*u(i)*v(i)
114  m(7) = m(7) + n(i)*u(i)*w(i)
115  m(8) = m(8) + n(i)*(v(i)**2)
116  m(9) = m(9) + n(i)*v(i)*w(i)
117  m(10)= m(10)+ n(i)*w(i)**2
118  m(11)= m(11)+ n(i)*(u(i)**3)
119  m(12)= m(12)+ n(i)*(u(i)**2)*v(i)
120  m(13)= m(13)+ n(i)*(u(i)**2)*w(i)
121  m(14)= m(14)+ n(i)*u(i)*(v(i)**2)
122  m(15)= m(15)+ n(i)*u(i)*v(i)*w(i)
123  m(16)= m(16)+ n(i)*u(i)*(w(i)**2)
124  m(17)= m(17)+ n(i)*(v(i)**3)
125  m(18)= m(18)+ n(i)*(v(i)**2)*w(i)
126  m(19)= m(19)+ n(i)*v(i)*(w(i)**2)
127  m(20)= m(20)+ n(i)*(w(i)**3)
128  END DO
129 
130  IF (m(1) < small_number) THEN
131  print *,'QMOMK: Zero order moment is very small!'
132  END IF
133 
134  mom = m
135  END SUBROUTINE moments_twenty_eight_nodes
136 
137 
138 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
139 ! C
140 ! Moments check C
141 ! C
142 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
143 
144  SUBROUTINE check_moments_twenty(mom, n, u, v, w, Check)
145 
146  IMPLICIT NONE
147 
148  DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NMOM) :: mom
149  DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: n
150  DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: u
151  DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: v
152  DOUBLE PRECISION, INTENT(IN), DIMENSION(QMOMK_NN) :: w
153  DOUBLE PRECISION, INTENT(OUT), DIMENSION(QMOMK_NMOM) :: Check
154 
155  DOUBLE PRECISION, DIMENSION(QMOMK_NMOM) :: F
156 
157  INTEGER :: i
158 
159  f = 0.d0
160 
161  f(1) = - mom(1)
162  f(2) = - mom(2)
163  f(3) = - mom(3)
164  f(4) = - mom(4)
165  f(5) = - mom(5)
166  f(6) = - mom(6)
167  f(7) = - mom(7)
168  f(8) = - mom(8)
169  f(9) = - mom(9)
170  f(10)= - mom(10)
171  f(11)= - mom(11)
172  f(12)= - mom(12)
173  f(13)= - mom(13)
174  f(14)= - mom(14)
175  f(15)= - mom(15)
176  f(16)= - mom(16)
177  f(17)= - mom(17)
178  f(18)= - mom(18)
179  f(19)= - mom(19)
180  f(20)= - mom(20)
181 
182  DO i = 1, qmomk_nn
183  f(1) = f(1) + n(i)
184  f(2) = f(2) + n(i)*u(i)
185  f(3) = f(3) + n(i)*v(i)
186  f(4) = f(4) + n(i)*w(i)
187  f(5) = f(5) + n(i)*(u(i)**2)
188  f(6) = f(6) + n(i)*u(i)*v(i)
189  f(7) = f(7) + n(i)*u(i)*w(i)
190  f(8) = f(8) + n(i)*(v(i)**2)
191  f(9) = f(9) + n(i)*v(i)*w(i)
192  f(10)= f(10)+ n(i)*(w(i)**2)
193  f(11)= f(11)+ n(i)*(u(i)**3)
194  f(12)= f(12)+ n(i)*(u(i)**2)*v(i)
195  f(13)= f(13)+ n(i)*(u(i)**2)*w(i)
196  f(14)= f(14)+ n(i)*u(i)*(v(i)**2)
197  f(15)= f(15)+ n(i)*u(i)*v(i)*w(i)
198  f(16)= f(16)+ n(i)*u(i)*(w(i)**2)
199  f(17)= f(17)+ n(i)*(v(i)**3)
200  f(18)= f(18)+ n(i)*(v(i)**2)*w(i)
201  f(19)= f(19)+ n(i)*v(i)*(w(i)**2)
202  f(20)= f(20)+ n(i)*(w(i)**3)
203  END DO
204 
205  DO i = 1, qmomk_nmom
206  f(i) = abs(f(i))
207  END DO
208 
209  check = f
210  END SUBROUTINE check_moments_twenty
211 
212 
213 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
214 ! C
215 ! Calculation of weights and abscissas C
216 ! C
217 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
218 
219  SUBROUTINE eight_node_3d (mom, n, u, v, w)
220 
221  IMPLICIT NONE
222 
223  DOUBLE PRECISION, INTENT(INOUT), DIMENSION(QMOMK_NMOM) :: mom
224  DOUBLE PRECISION, INTENT(OUT), DIMENSION(QMOMK_NN) :: n
225  DOUBLE PRECISION, INTENT(OUT), DIMENSION(QMOMK_NN) :: u
226  DOUBLE PRECISION, INTENT(OUT), DIMENSION(QMOMK_NN) :: v
227  DOUBLE PRECISION, INTENT(OUT), DIMENSION(QMOMK_NN) :: w
228 
229  DOUBLE PRECISION :: um, vm, wm, s1, s2, s3, s12, s13, s23
230  DOUBLE PRECISION :: q000, q101, q110, q011, q200, q020, q002, q210, q300
231  DOUBLE PRECISION :: q201, q120, q111, q102, q030, q021, q012, q003
232 
233  DOUBLE PRECISION :: b11, b12, b13, b21, b22, b23, b31, b32, b33, Fmax
234 
235  DOUBLE PRECISION :: m000, m100, m010, m001, m200, m020, m002, m300, m030, m003, m111
236 
237  DOUBLE PRECISION, DIMENSION (QMOMK_NCOV, QMOMK_NCOV) :: sig
238  DOUBLE PRECISION, DIMENSION (QMOMK_NCOV, QMOMK_NCOV) :: chol
239  DOUBLE PRECISION, DIMENSION (QMOMK_NCOV, QMOMK_NCOV) :: scal, invScale
240  DOUBLE PRECISION, DIMENSION (QMOMK_NCOV, QMOMK_NCOV) :: Umat, Vmat
241 
242  DOUBLE PRECISION, DIMENSION (4) :: alpha1, alpha2, alpha3
243  DOUBLE PRECISION, DIMENSION (2) :: aq1, aq2, aq3, wq1, wq2, wq3
244 
245  DOUBLE PRECISION :: a, b, c, u1, u2, v1, v2, w1, w2, rho123
246 
247  DOUBLE PRECISION, DIMENSION(QMOMK_NN) :: nstar, ones
248 
249  DOUBLE PRECISION, DIMENSION(32) :: x
250 
251  DOUBLE PRECISION, DIMENSION(QMOMK_NMOM) :: F
252 
253  q000 = mom(1)
254  um = mom(2)/q000
255  vm = mom(3)/q000
256  wm = mom(4)/q000
257  q200 = mom(5)/q000
258  q110 = mom(6)/q000
259  q101 = mom(7)/q000
260  q020 = mom(8)/q000
261  q011 = mom(9)/q000
262  q002 = mom(10)/q000
263  q300 = mom(11)/q000
264  q210 = mom(12)/q000
265  q201 = mom(13)/q000
266  q120 = mom(14)/q000
267  q111 = mom(15)/q000
268  q102 = mom(16)/q000
269  q030 = mom(17)/q000
270  q021 = mom(18)/q000
271  q012 = mom(19)/q000
272  q003 = mom(20)/q000
273 
274  sig = 0.d0
275 
276  sig(1,1) = q200 - um**2
277  sig(2,2) = q020 - vm**2
278  sig(3,3) = q002 - wm**2
279 
280  IF ( sig(1,1) < 0.d0 ) THEN ! check that u variance is nonnegative
281  print *,'QMOMK: Negative u variance',sig
282  sig(1,1) = 0.d0
283  error_stop
284  END IF
285 
286  IF ( sig(2,2) < 0.d0 ) THEN ! check that v variance is nonnegative
287  print *,'QMOMK: Negative v variance',sig
288  sig(2,2) = 0.d0
289  error_stop
290  END IF
291 
292  IF ( sig(3,3) < 0.d0 ) THEN ! check that v variance is nonnegative
293  print *,'QMOMK: Negative w variance',sig
294  sig(3,3) = 0.d0
295  error_stop
296  END IF
297 
298  s1 = sqrt(sig(1,1))
299  s2 = sqrt(sig(2,2))
300  s3 = sqrt(sig(3,3))
301 
302  sig(1,2) = q110 - um*vm
303  sig(2,1) = sig(1,2)
304  s12 = sig(1,2)
305 
306  sig(1,3) = q101 - um*wm
307  sig(3,1) = sig(1,3)
308  s13 = sig(1,3)
309 
310  sig(2,3) = q011 - vm*wm
311  sig(3,2) = sig(2,3)
312  s23 = sig(2,3)
313 
314  IF ( s1*s2 < abs(s12) ) THEN ! check that covariance matrix is nonnegative
315  print *,'QMOMK: unphysical uv correlation', sig(1,2), s1, s2
316  sig(1,2) = sign(1.d0, s12)*s1*s2
317  sig(2,1) = sig(1,2)
318 
319  !q110 = 0.99*q110
320  !mom(6) = q000*q110
321  !sig(1,2) = q110 - um*vm
322  !sig(2,1) = sig(1,2)
323  print *,'Sigma = ',sig
324  print *,'M = ', mom
325  print *,'N = ', n
326  print *,'U = ', u
327  print *,'V = ', v
328  print *,'W = ', w
329  error_stop
330  END IF
331 
332  IF ( s1*s3 < abs(s13) ) THEN ! check that covariance matrix is nonnegative
333  sig(1,3) = sign(1.d0, s13)*s1*s3
334  sig(3,1) = sig(1,3)
335  print *,'QMOMK: unphysical uw correlation',sig(1,3)
336  error_stop
337  END IF
338 
339  IF ( s2*s3 < abs(s23) ) THEN ! check that covariance matrix is nonnegative
340  sig(2,3) = sign(1.d0, s23)*s2*s3
341  sig(3,2) = sig(2,3)
342  print *,'QMOMK: unphysical vw correlation'
343  error_stop
344  END IF
345 
346  CALL cholesky3(sig, chol) ! Bottom cholesky decomposition
347  CALL diag3(s1, s2, s3, scal)
348  CALL inv3(scal, invscale)
349  CALL multiplymatrix3(chol, invscale, umat)
350  CALL inv3(umat, vmat)
351 
352  b11 = vmat(1,1)
353  b12 = vmat(1,2)
354  b13 = vmat(1,3)
355 
356  b21 = vmat(2,1)
357  b22 = vmat(2,2)
358  b23 = vmat(2,3)
359 
360  b31 = vmat(3,1)
361  b32 = vmat(3,2)
362  b33 = vmat(3,3)
363 
364  m000 = 1
365  m100 = 0
366  m010 = 0
367  m001 = 0
368 
369  m200 = 2.d0*b11*b12*q110 - 2.d0*b11*um*b12*vm + 2.d0*b11*b13*q101 - &
370  2.d0*b11*um*b13*wm + 2.d0*b12*b13*q011 - 2.d0*b12*vm*b13*wm - &
371  (b11**2)*(um**2) - (b12**2)*(vm**2) - (b13**2)*(wm**2) + &
372  (b11**2)*q200 + (b12**2)*q020 + (b13**2)*q002
373 
374  m020 = - (b21**2)*(um**2) - (b22**2)*(vm**2) - (b23**2)*(wm**2) + &
375  2.d0*b21*b23*q101 - 2.d0*b21*um*b23*wm + 2.d0*b22*b23*q011 - &
376  2.d0*b22*vm*b23*wm - 2.d0*b21*um*b22*vm + 2.d0*b21*b22*q110 + &
377  (b21**2)*q200 + (b22**2)*q020 + (b23**2)*q002
378 
379  m002 = - (b33**2)*(wm**2) - (b31**2)*(um**2) - (b32**2)*(vm**2) - &
380  2.d0*b31*um*b32*vm + 2.d0*b31*b32*q110 + (b31**2)*q200 + &
381  (b32**2)*q020 + (b33**2)*q002 + 2.d0*b31*b33*q101 - 2.d0*b31*um*b33*wm + &
382  2.d0*b32*b33*q011 - 2.d0*b32*vm*b33*wm
383 
384  m300 = - 3.d0*(b13**3)*wm*q002 + 3.d0*b11*(b12**2)*q120 + 3.d0*(b12**2)*b13*q021 + &
385  3.d0*(b11**2)*b12*q210 + 3.d0*(b11**2)*b13*q201 + 3.d0*b11*(b13**2)*q102 + &
386  3.d0*b12*(b13**2)*q012 + 6.d0*b12*vm*(b13**2)*(wm**2) - 3.d0*(b12**3)*vm*q020 - &
387  3.d0*(b11**3)*um*q200 - 6.d0*b11*b12*b13*wm*q110 - 6.d0*b11*b12*vm*b13*q101 + &
388  (b11**3)*q300 + 2.d0*(b11**3)*(um**3) - 6.d0*b11*um*b12*b13*q011 + 2.d0*(b12**3)*(vm**3) + &
389  2.d0*(b13**3)*(wm**3) + (b13**3)*q003 + (b12**3)*q030 - 3.d0*b11*um*(b12**2)*q020 - &
390  3.d0*b11*um*(b13**2)*q002 - 3.d0*(b12**2)*b13*wm*q020 + 6.d0*(b11**2)*(um**2)*b12*vm + &
391  6.d0*(b12**2)*(vm**2)*b13*wm + 6.d0*(b11**2)*(um**2)*b13*wm + 6.d0*b11*um*(b12**2)*(vm**2) + &
392  6.d0*b11*um*(b13**2)*(wm**2) - 6.d0*b11*(b12**2)*vm*q110 - 6.d0*b11*(b13**2)*wm*q101 - &
393  6.d0*(b12**2)*vm*b13*q011 + 12.d0*b11*um*b12*vm*b13*wm - 6.d0*b12*(b13**2)*wm*q011 - &
394  6.d0*(b11**2)*um*b13*q101 + 6.d0*b11*b12*b13*q111 - 3.d0*b12*vm*(b13**2)*q002 - &
395  6.d0*(b11**2)*um*b12*q110 - 3.d0*(b11**2)*b13*wm*q200 - 3.d0*(b11**2)*b12*vm*q200
396 
397  m030 = 6.d0*(b21**2)*(um**2)*b22*vm + 6.d0*(b22**2)*(vm**2)*b23*wm - 3.d0*(b21**3)*um*q200 + &
398  2.d0*(b22**3)*(vm**3) + 3.d0*(b21**2)*b23*q201 + 2.d0*(b21**3)*(um**3) - 3.d0*(b23**3)*wm*q002 - &
399  3.d0*(b21**2)*b22*vm*q200 + 6.d0*b21*um*(b22**2)*(vm**2) + 3.d0*b21*(b22**2)*q120 + &
400  3.d0*(b21**2)*b22*q210 + 6.d0*(b21**2)*(um**2)*b23*wm - 3.d0*b21*um*(b22**2)*q020 + &
401  3.d0*b21*(b23**2)*q102 + 3.d0*(b22**2)*b23*q021 + 6.d0*b22*vm*(b23**2)*(wm**2) - &
402  3.d0*(b21**2)*b23*wm*q200 + (b22**3)*q030 + (b21**3)*q300 - 3.d0*(b22**3)*vm*q020 + &
403  6.d0*b21*b22*b23*q111 - 6.d0*(b21**2)*um*b22*q110 + (b23**3)*q003 - 3.d0*b22*vm*(b23**2)*q002 - &
404  6.d0*(b21**2)*um*b23*q101 + 3.d0*b22*(b23**2)*q012 - 6.d0*(b22**2)*vm*b23*q011 - &
405  3.d0*b21*um*(b23**2)*q002 - 3.d0*(b22**2)*b23*wm*q020 - 6.d0*b22*(b23**2)*wm*q011 - &
406  6.d0*b21*b22*b23*wm*q110 - 6.d0*b21*um*b22*b23*q011 + 6.d0*b21*um*(b23**2)*(wm**2) + &
407  2.d0*(b23**3)*(wm**3) - 6.d0*b21*(b23**2)*wm*q101 - 6.d0*b21*b22*vm*b23*q101 + &
408  12.d0*b21*um*b22*vm*b23*wm - 6.d0*b21*(b22**2)*vm*q110
409 
410  m003 = (b31**3)*q300 + 6.d0*b31*um*(b33**2)*(wm**2) + (b32**3)*q030 + (b33**3)*q003 - &
411  3.d0*b31*um*(b33**2)*q002 + 6.d0*b32*vm*(b33**2)*(wm**2) + 6.d0*b31*um*(b32**2)*(vm**2) + &
412  2.d0*(b31**3)*(um**3) + 2.d0*(b32**3)*(vm**3) + 3.d0*(b31**2)*b33*q201 - &
413  3.d0*b32*vm*(b33**2)*q002 - 6.d0*b31*(b32**2)*vm*q110 - 3.d0*(b31**2)*b32*vm*q200 - &
414  3.d0*(b31**2)*b33*wm*q200 - 6.d0*b31*(b33**2)*wm*q101 - 3.d0*(b33**3)*wm*q002 + &
415  6.d0*(b31**2)*(um**2)*b33*wm + 3.d0*b32*(b33**2)*q012 + 6.d0*b31*b32*b33*q111 - &
416  6.d0*(b31**2)*um*b33*q101 - 6.d0*b32*(b33**2)*wm*q011 + 6.d0*(b31**2)*(um**2)*b32*vm - &
417  3.d0*b31*um*(b32**2)*q020 - 6.d0*b31*b32*b33*wm*q110 + 3.d0*b31*(b33**2)*q102 - &
418  3.d0*(b32**2)*b33*wm*q020 + 12.d0*b31*um*b32*vm*b33*wm + 3.d0*(b31**2)*b32*q210 - &
419  6.d0*(b32**2)*vm*b33*q011 + 3.d0*b31*(b32**2)*q120 + 6.d0*(b32**2)*(vm**2)*b33*wm + &
420  3.d0*(b32**2)*b33*q021 - 3.d0*(b31**3)*um*q200 - 3.d0*(b32**3)*vm*q020 - &
421  6.d0*b31*um*b32*b33*q011 + 2.d0*(b33**3)*(wm**3) - 6.d0*b31*b32*vm*b33*q101 - &
422  6.d0*(b31**2)*um*b32*q110
423 
424  m111 = - b12*b23*wm*b31*q110 + 2.d0*b11*(um**2)*b22*vm*b31 - b12*b21*b33*wm*q110 - b12*vm*b23*b33*q002 + &
425  2.d0*b12*vm*b23*wm*b31*um - b12*vm*b21*b31*q200 - b13*b23*b31*um*q002 - b12*b22*b31*um*q020 + &
426  b12*b22*b31*q120 - b11*b21*b33*wm*q200 + 2.d0*b12*(vm**2)*b23*wm*b32 - b12*b23*wm*b32*q020 - &
427  b11*b21*b32*vm*q200 - b12*b21*um*b32*q020 - 2.d0*b12*b22*vm*b33*q011 - b11*um*b23*b33*q002 - &
428  2.d0*b12*b22*vm*b31*q110 - 3.d0*b11*b21*b31*um*q200 + b11*b21*b32*q210 + b11*b21*b31*q300 - &
429  2.d0*b12*b23*b33*wm*q011 + b13*b23*b33*q003 + b12*b23*b33*q012 + b12*b22*b33*q021 - &
430  b12*b23*b31*um*q011 - b12*b22*b33*wm*q020 + b13*b23*b31*q102 - 2.d0*b12*b21*b31*um*q110 + &
431  b12*b21*b32*q120 + b12*b23*b31*q111 - 2.d0*b12*b21*b32*vm*q110 + 2.d0*b12*(vm**2)*b21*um*b32 - &
432  3.d0*b12*b22*b32*vm*q020 + 2.d0*b12*(vm**2)*b22*b31*um + 2.d0*b12*vm*b21*(um**2)*b31 - &
433  2.d0*b12*b23*b32*vm*q011 - b12*b21*um*b33*q011 + 2.d0*b12*(vm**2)*b22*b33*wm + b12*b23*b32*q021 + &
434  b11*b21*b33*q201 - 3.d0*b13*b23*b33*wm*q002 - 2.d0*b11*b21*um*b32*q110 + 2.d0*b12*vm*b23*(wm**2)*b33 - &
435  2.d0*b13*b21*b33*wm*q101 + 2.d0*b13*(wm**2)*b22*vm*b33 - b13*b21*um*b32*q011 - b12*vm*b21*b33*q101 - &
436  2.d0*b11*b21*um*b33*q101 - b13*b22*vm*b33*q002 - b13*b21*um*b33*q002 + 2.d0*b11*(um**3)*b21*b31 + &
437  2.d0*b11*b21*(um**2)*b32*vm + b13*b21*b31*q201 + b13*b22*b33*q012 + b13*b22*b31*q111 + &
438  2.d0*b13*b23*(wm**2)*b31*um - b13*b23*b32*vm*q002 + b11*b23*b31*q201 - b11*b22*vm*b31*q200 + &
439  2.d0*b13*wm*b21*(um**2)*b31 - b13*b22*vm*b31*q101 + 2.d0*b11*b21*(um**2)*b33*wm + b13*b23*b32*q012 - &
440  2.d0*b11*b22*b31*um*q110 - 2.d0*b13*b22*b32*vm*q011 + b12*b21*b31*q210 - b13*wm*b22*b32*q020 - &
441  2.d0*b13*b23*wm*b31*q101 + b11*b22*b33*q111 + 2.d0*b13*(wm**2)*b21*um*b33 - b13*b21*b32*vm*q101 + &
442  b12*b22*b32*q030 + 2.d0*b13*b23*(wm**2)*b32*vm - b11*b22*b33*wm*q110 - b12*vm*b23*b31*q101 - &
443  b13*wm*b21*b32*q110 + 2.d0*b11*um*b23*wm*b32*vm + 2.d0*b11*um*b22*(vm**2)*b32 + 2*b13*(wm**3)*b23*b33 + &
444  2.d0*b12*(vm**3)*b22*b32 - 2.d0*b11*b23*b31*um*q101 - b11*b22*vm*b33*q101 - b11*um*b22*b32*q020 + &
445  2.d0*b13*wm*b22*(vm**2)*b32 + b13*b21*b33*q102 + b13*b21*b32*q111 - 2.d0*b13*b21*b31*um*q101 - &
446  b13*b22*b31*um*q011 + b11*b22*b32*q120 - b13*wm*b22*b31*q110 + 2.d0*b11*um*b23*(wm**2)*b33 + &
447  2.d0*b11*(um**2)*b23*wm*b31 + b11*b23*b32*q111 + b13*b22*b32*q021 + 2.d0*b12*vm*b21*um*b33*wm + &
448  2.d0*b13*wm*b22*vm*b31*um - 2.d0*b13*b23*wm*b32*q011 + b12*b21*b33*q111 - 2.d0*b11*b23*b33*wm*q101 + &
449  2.d0*b11*um*b22*vm*b33*wm - b11*um*b23*b32*q011 - b11*um*b22*b33*q011 - 2.d0*b11*b22*b32*vm*q110 + &
450  b11*b22*b31*q210 - b13*wm*b21*b31*q200 - 2.d0*b13*b22*b33*wm*q011 - b11*b23*wm*b32*q110 + &
451  b11*b23*b33*q102 - b11*b23*wm*b31*q200 + 2.d0*b13*wm*b21*um*b32*vm - b11*b23*b32*vm*q101
452 
453  alpha1(1) = m000
454  alpha1(2) = m100
455  alpha1(3) = m200
456  alpha1(4) = m300
457 
458  CALL quadrature_bounded(alpha1, wq1, aq1)
459 
460  alpha2(1) = m000
461  alpha2(2) = m010
462  alpha2(3) = m020
463  alpha2(4) = m030
464 
465  call quadrature_bounded(alpha2, wq2, aq2)
466 
467  alpha3(1) = m000
468  alpha3(2) = m001
469  alpha3(3) = m002
470  alpha3(4) = m003
471 
472  call quadrature_bounded(alpha3, wq3, aq3)
473 
474  a = wq1(1)
475  b = wq2(1)
476  c = wq3(1)
477 
478  u1 = aq1(1)
479  u2 = aq1(2)
480 
481  v1 = aq2(1)
482  v2 = aq2(2)
483 
484  w1 = aq3(1)
485  w2 = aq3(2)
486 
487  rho123 = sqrt(a*(1-a)*b*(1-b)*c*(1-c))*m111/sqrt(m200*m020*m002)
488 
489  IF (rho123 > 0.d0) THEN
490  rho123 = min(rho123 , a*b*c , (1-a)*(1-b)*c , (1-a)*b*(1-c) , a*(1-b)*(1-c))
491  ELSE
492  rho123 = -min(-rho123 , (1-a)*b*c , a*(1-b)*c , a*b*(1-c) , (1-a)*(1-b)*(1-c))
493  END IF
494 
495  IF (a > 1.) THEN
496  print *,a
497  END IF
498  IF (b > 1.) THEN
499  print *,b
500  END IF
501  IF (c > 1.) THEN
502  print *,c
503  END IF
504 
505  nstar(1) = a*b*c - rho123
506  nstar(2) = (1.d0-a)*b*c + rho123
507  nstar(3) = a*(1.d0-b)*c + rho123
508  nstar(4) = (1.d0-a)*(1.d0-b)*c - rho123
509  nstar(5) = a*b*(1.d0-c) + rho123
510  nstar(6) = (1.d0-a)*b*(1.d0-c) - rho123
511  nstar(7) = a*(1.d0-b)*(1.d0-c) - rho123
512  nstar(8) = (1.d0-a)*(1.d0-b)*(1.d0-c) + rho123
513 
514  IF (minval(nstar) < -1.d-16) THEN
515  print *,'QMOMK: Negative weight in quadrature'
516  print *,nstar
517  nstar = abs(nstar)
518  error_stop
519  END IF
520 
521  x(1:8) = nstar
522  x(9:16) = (/u1, u2, u1, u2, u1, u2, u1, u2/)
523  x(17:24) = (/v1, v1, v2, v2, v1, v1, v2, v2/)
524  x(25:32) = (/w1, w1, w1, w1, w2, w2, w2, w2/)
525 
526  ones = (/1.,1.,1.,1.,1.,1.,1.,1./)
527 
528  n = x(1:8)*q000
529  u = umat(1,1)*x(9:16) + umat(1,2)*x(17:24) + umat(1,3)*x(25:32) + um*ones(1:8)
530  v = umat(2,1)*x(9:16) + umat(2,2)*x(17:24) + umat(2,3)*x(25:32) + vm*ones(1:8)
531  w = umat(3,1)*x(9:16) + umat(3,2)*x(17:24) + umat(3,3)*x(25:32) + wm*ones(1:8)
532 
533  CALL check_moments_twenty(mom, n, u, v, w, f)
534  fmax = maxval(f(1:10))
535 
536  IF (fmax > 1.d-10) THEN
537  print *,'QMOMK: First 10 moments not satisfied.'
538  print *,'F(1:10) = ',f(1:10)
539  print *,'M(1:10) = ',mom(1:10)
540  error_stop
541  ENDIF
542  END SUBROUTINE eight_node_3d
543 
544 
545 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
546 ! C
547 ! Temperature limiter C
548 ! C
549 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
550 
551  SUBROUTINE bind_theta (mom, tau_p)
552 
553  DOUBLE PRECISION, INTENT(INOUT), DIMENSION(QMOMK_NMOM) :: mom
554  DOUBLE PRECISION, INTENT(IN) :: tau_p
555 
556  DOUBLE PRECISION :: sigma2
557 
558  sigma2 = min(minimum_theta/tau_p, maximum_sigma)
559 
560  ! Second order moments
561  mom(5) = mom(5) + 2.0*sigma2*mom(1)
562  mom(8) = mom(8) + 2.0*sigma2*mom(1)
563  mom(10) = mom(10) + 2.0*sigma2*mom(1)
564 
565  ! Third order moments
566  mom(11) = mom(11) + 6.0*sigma2*mom(2)
567  mom(12) = mom(12) + 2.0*sigma2*mom(3)
568  mom(13) = mom(13) + 2.0*sigma2*mom(4)
569  mom(14) = mom(14) + 2.0*sigma2*mom(2)
570 
571  mom(16) = mom(16) + 2.0*sigma2*mom(2)
572  mom(17) = mom(17) + 6.0*sigma2*mom(3)
573  mom(18) = mom(18) + 2.0*sigma2*mom(4)
574 
575  mom(19) = mom(19) + 2.0*sigma2*mom(3)
576  mom(20) = mom(20) + 6.0*sigma2*mom(4)
577 
578  END SUBROUTINE bind_theta
579 
580 END MODULE qmomk_quadrature
qmomk_quadrature::moments_twenty_eight_nodes
subroutine, public moments_twenty_eight_nodes(n, u, v, w, mom)
Definition: qmomk_quadrature_mod.f:91
qmomk_parameters::qmomk_nmom
integer, parameter qmomk_nmom
Definition: qmomk_parameters_mod.f:9
qmomk_quadrature::check_moments_twenty
subroutine check_moments_twenty(mom, n, u, v, w, Check)
Definition: qmomk_quadrature_mod.f:145
qmomk_parameters
Definition: qmomk_parameters_mod.f:1
qmomk_tools::diag3
subroutine, public diag3(a, b, c, mat)
Definition: qmomk_tools_mod.f:70
qmomk_tools::inv3
subroutine, public inv3(A, B)
Definition: qmomk_tools_mod.f:41
qmomk_parameters::maximum_sigma
double precision, parameter maximum_sigma
Definition: qmomk_parameters_mod.f:35
qmomk_quadrature::eight_node_3d
subroutine, public eight_node_3d(mom, n, u, v, w)
Definition: qmomk_quadrature_mod.f:220
qmomk_parameters::eps
double precision, parameter eps
Definition: qmomk_parameters_mod.f:16
param1::small_number
double precision, parameter small_number
Definition: param1_mod.f:24
qmomk_parameters::minimum_theta
double precision, parameter minimum_theta
Definition: qmomk_parameters_mod.f:34
qmomk_tools
Definition: qmomk_tools_mod.f:13
qmomk_tools::cholesky3
subroutine, public cholesky3(A, L)
Definition: qmomk_tools_mod.f:89
qmomk_quadrature
Definition: qmomk_quadrature_mod.f:16
qmomk_quadrature::bind_theta
subroutine, public bind_theta(mom, tau_p)
Definition: qmomk_quadrature_mod.f:552
qmomk_parameters::qmomk_nn
integer, parameter qmomk_nn
Definition: qmomk_parameters_mod.f:11
qmomk_tools::multiplymatrix3
subroutine, public multiplymatrix3(A, B, P)
Definition: qmomk_tools_mod.f:140