MFIX  2016-1
check_data_cartesian.f
Go to the documentation of this file.
1 MODULE check_data_cg
2  USE exit, only: mfix_exit
3 
4  CONTAINS
5 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
6 ! C
7 ! Module name: CHECK_DATA_CARTESIAN C
8 ! Purpose: check the data related to cartesian grid implementation C
9 ! C
10 ! Author: Jeff Dietiker Date: 21-Feb-08 C
11 ! Reviewer: Date: C
12 ! C
13 ! Revision Number # Date: ##-###-## C
14 ! Author: # C
15 ! Purpose: # C
16 ! C
17 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
18 
19  SUBROUTINE check_data_cartesian
20 
21 !-----------------------------------------------
22 ! Modules
23 !-----------------------------------------------
24  USE bc
25  USE compar
26  USE constant
27  USE cutcell
28  USE dashboard
29  USE discretelement
30  USE funits
31  USE indices
32  USE iterate, ONLY: max_nit
33  USE leqsol
34  USE mpi_utility
35  USE param
36  USE param1
37  USE physprop
38  USE polygon
39  USE quadric
40  USE run
41  USE scalars
42  USE stl
43  USE vtk
44  IMPLICIT NONE
45 !-----------------------------------------------
46 ! Local variables
47 !-----------------------------------------------
48  INTEGER :: I,J,Q,BCV
49  DOUBLE PRECISION :: norm, tan_half_angle
50  CHARACTER(LEN=9) :: GR
51 !-----------------------------------------------
52 
53  IF(.NOT.cartesian_grid) RETURN
54 
55  IF(discrete_element) THEN
56  IF(mype == pe_io) THEN
57 
58  WRITE(*,10)'######################################################################'
59  WRITE(*,10)'## ##'
60  WRITE(*,10)'## ===> WARNING: RUNNING CARTESIAN GRID WITH DISCRETE ELEMENT. ##'
61  WRITE(*,10)'## THIS HAS NOT BEEN FULLY TESTED. ##'
62  WRITE(*,10)'## PLEASE USE WITH CAUTION. ##'
63  WRITE(*,10)'## ##'
64  WRITE(*,10)'######################################################################'
65 
66  ENDIF
67  ENDIF
68 
69 10 FORMAT(1x,a)
70 
71  IF(coordinates=='CYLINDRICAL') THEN
72  IF(mype == pe_io) THEN
73  WRITE(*,*)'INPUT ERROR: CARTESIAN GRID OPTION NOT AVAILABLE'
74  WRITE(*,*)'WITH CYLINDRICAL COORDINATE SYSTEM.'
75  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
76  ENDIF
77  CALL mfix_exit(mype)
78  ENDIF
79 
80  IF(use_stl.AND.(.NOT.use_msh)) THEN
81  IF(do_k) THEN
82  CALL get_stl_data
83  ELSE
84  IF(mype == pe_io) WRITE(*,*) &
85  'ERROR: STL METHOD VALID ONLY IN 3D.'
86  CALL mfix_exit(mype)
87  ENDIF
88  IF(n_quadric > 0) THEN
89  IF(mype == pe_io) WRITE(*,*) 'ERROR: BOTH QUADRIC(S) AND ',&
90  'STL INPUT ARE SPECIFIED.'
91  IF(mype == pe_io) WRITE(*,*) 'MFIX HANDLES ONLY ONE TYPE ',&
92  'OF SURFACE INPUT.'
93  CALL mfix_exit(mype)
94  ENDIF
95  ENDIF
96 
97  IF(use_msh.AND.(.NOT.use_stl)) THEN
98  IF(do_k) THEN
99  CALL get_msh_data
100  ELSE
101  IF(mype == pe_io) WRITE(*,*) &
102  'ERROR: MSH METHOD VALID ONLY IN 3D.'
103  CALL mfix_exit(mype)
104  ENDIF
105  IF(n_quadric > 0) THEN
106  IF(mype == pe_io) WRITE(*,*) 'ERROR: BOTH QUADRIC(S) AND ',&
107  'MSH INPUT ARE SPECIFIED.'
108  IF(mype == pe_io) WRITE(*,*) 'MFIX HANDLES ONLY ONE TYPE ',&
109  'OF SURFACE INPUT.'
110  CALL mfix_exit(mype)
111  ENDIF
112  ENDIF
113 
114  IF(use_polygon) THEN
115  IF(do_k) THEN
116  IF(mype == pe_io) WRITE(*,*) 'ERROR: POLYGON METHOD ',&
117  'VALID ONLY IN 2D.'
118  CALL mfix_exit(mype)
119  ELSE
120  CALL get_poly_data
121  ENDIF
122  ENDIF
123 
124  IF(n_quadric > 0) THEN
125  IF(n_polygon > 0) THEN
126  IF(mype == pe_io) THEN
127  WRITE(*,*) 'ERROR: BOTH QUADRIC(S) AND POLYGON(S) ',&
128  'DEFINED.'
129  WRITE(*,*) 'MFIX HANDLES ONLY ONE TYPE OF SURFACE INPUT.'
130  ENDIF
131  CALL mfix_exit(mype)
132  ENDIF
133  IF(n_usr_def > 0) THEN
134  IF(mype == pe_io) THEN
135  WRITE(*,*) 'ERROR: BOTH QUADRIC(S) AND USER-DEFINED ',&
136  'FUNCTION DEFINED.'
137  WRITE(*,*) 'MFIX HANDLES ONLY ONE TYPE OF SURFACE.'
138  ENDIF
139  CALL mfix_exit(mype)
140  ENDIF
141  IF(quadric_scale <= zero) THEN
142  IF(mype == pe_io) THEN
143  WRITE(*,*) 'ERROR: QUADRIC_SCALE MUST BE POSITIVE.'
144  ENDIF
145  CALL mfix_exit(mype)
146  ELSEIF(quadric_scale /= one) THEN
147  DO q = 1, n_quadric
148  lambda_x(q) = lambda_x(q) * quadric_scale**2
149  lambda_y(q) = lambda_y(q) * quadric_scale**2
150  lambda_z(q) = lambda_z(q) * quadric_scale**2
151  radius(q) = radius(q) * quadric_scale
152  t_x(q) = t_x(q) * quadric_scale
153  t_y(q) = t_y(q) * quadric_scale
154  t_z(q) = t_z(q) * quadric_scale
155  clip_xmin(q) = clip_xmin(q) * quadric_scale
156  clip_xmax(q) = clip_xmax(q) * quadric_scale
157  clip_ymin(q) = clip_ymin(q) * quadric_scale
158  clip_ymax(q) = clip_ymax(q) * quadric_scale
159  clip_zmin(q) = clip_zmin(q) * quadric_scale
160  clip_zmax(q) = clip_zmax(q) * quadric_scale
161  piece_xmin(q) = piece_xmin(q) * quadric_scale
162  piece_xmax(q) = piece_xmax(q) * quadric_scale
163  piece_ymin(q) = piece_ymin(q) * quadric_scale
164  piece_ymax(q) = piece_ymax(q) * quadric_scale
165  piece_zmin(q) = piece_zmin(q) * quadric_scale
166  piece_zmax(q) = piece_zmax(q) * quadric_scale
167  ENDDO
168  ENDIF
169  ELSE
170  IF((n_polygon > 0).AND.(n_usr_def > 0)) THEN
171  IF(mype == pe_io) THEN
172  WRITE(*,*) 'ERROR: POLYGON(S) AND USER-DEFINED ',&
173  'FUNCTION DEFINED.'
174  WRITE(*,*) 'MFIX HANDLES ONLY ONE TYPE OF SURFACE.'
175  ENDIF
176  CALL mfix_exit(mype)
177  ENDIF
178  ENDIF
179 
180 
181  IF(n_quadric > dim_quadric) THEN
182  IF(mype == pe_io) THEN
183  WRITE(*,*)'INPUT ERROR: INVALID VALUE OF N_QUADRIC =', &
184  n_quadric
185  WRITE(*,*)'MAXIMUM ACCEPTABLE VALUE IS DIM_QUADRIC =', &
186  dim_quadric
187  WRITE(*,*)'CHANGE MAXIMUM VALUE IN QUADRIC_MOD.F ',&
188  'IF NECESSARY.'
189  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
190  ENDIF
191  CALL mfix_exit(mype)
192  ENDIF
193 
194 
195  DO q = 1, n_quadric
196 
197 
198 
199  SELECT CASE (trim(quadric_form(q)))
200 
201  CASE ('NORMAL')
202 
203  lambda_x(q) = lambda_x(q)
204  lambda_y(q) = lambda_y(q)
205  lambda_z(q) = lambda_z(q)
206 
207  norm = dsqrt(lambda_x(q)**2 + lambda_y(q)**2 + &
208  lambda_z(q)**2)
209 
210  IF(norm < tol_f) THEN
211  IF(mype == pe_io) THEN
212  WRITE(*,*)'INPUT ERROR: QUADRIC:', q, &
213  ' HAS ZERO COEFFICIENTS.'
214  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
215  ENDIF
216  CALL mfix_exit(mype)
217  ENDIF
218 
219  CASE ('PLANE') ! The quadric is predefined as a plane
220 
221  lambda_x(q) = n_x(q)
222  lambda_y(q) = n_y(q)
223  lambda_z(q) = n_z(q)
224 
225  norm = dsqrt(lambda_x(q)**2 + lambda_y(q)**2 + &
226  lambda_z(q)**2)
227 
228  IF( norm > tol_f) THEN
229  lambda_x(q) = lambda_x(q) / norm
230  lambda_y(q) = lambda_y(q) / norm
231  lambda_z(q) = lambda_z(q) / norm
232  ELSE
233  IF(mype == pe_io) THEN
234  WRITE(*,*)'INPUT ERROR: PLANE:', q, &
235  ' HAS ZERO NORMAL VECTOR.'
236  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
237  ENDIF
238  CALL mfix_exit(mype)
239  ENDIF
240 
241  dquadric(q) = - (lambda_x(q)*t_x(q) + lambda_y(q)*t_y(q) + &
242  lambda_z(q)*t_z(q))
243 
244  CASE ('X_CYL_INT') ! The quadric is predefined as a cylinder, along x-axis
245  ! Internal flow
246 
247  IF( radius(q) <= zero) THEN
248  IF(mype == pe_io) THEN
249  WRITE(*,*)'INPUT ERROR: CYLINDER:', q, &
250  ' HAS ZERO RADIUS.'
251  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
252  ENDIF
253  CALL mfix_exit(mype)
254  ELSE
255  lambda_x(q) = zero
256  lambda_y(q) = one
257  lambda_z(q) = one
258  dquadric(q) = -radius(q)**2
259  ENDIF
260 
261  CASE ('Y_CYL_INT') ! The quadric is predefined as a cylinder, along y-axis
262  ! Internal flow
263 
264  IF( radius(q) <= zero) THEN
265  IF(mype == pe_io) THEN
266  WRITE(*,*)'INPUT ERROR: CYLINDER:', q, &
267  ' HAS ZERO RADIUS.'
268  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
269  ENDIF
270  CALL mfix_exit(mype)
271  ELSE
272  lambda_x(q) = one
273  lambda_y(q) = zero
274  lambda_z(q) = one
275  dquadric(q) = -radius(q)**2
276  ENDIF
277 
278  CASE ('Z_CYL_INT') ! The quadric is predefined as a cylinder, along z-axis
279  ! Internal flow
280 
281  IF( radius(q) <= zero) THEN
282  IF(mype == pe_io) THEN
283  WRITE(*,*)'INPUT ERROR: CYLINDER:', q, &
284  ' HAS ZERO RADIUS.'
285  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
286  ENDIF
287  CALL mfix_exit(mype)
288  ELSE
289  lambda_x(q) = one
290  lambda_y(q) = one
291  lambda_z(q) = zero
292  dquadric(q) = -radius(q)**2
293  ENDIF
294 
295 
296  CASE ('X_CYL_EXT') ! The quadric is predefined as a cylinder, along x-axis
297  ! External flow
298 
299  IF( radius(q) <= zero) THEN
300  IF(mype == pe_io) THEN
301  WRITE(*,*)'INPUT ERROR: CYLINDER:', q, &
302  ' HAS ZERO RADIUS.'
303  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
304  ENDIF
305  CALL mfix_exit(mype)
306  ELSE
307  lambda_x(q) = zero
308  lambda_y(q) = -one
309  lambda_z(q) = -one
310  dquadric(q) = radius(q)**2
311  ENDIF
312 
313  CASE ('Y_CYL_EXT') ! The quadric is predefined as a cylinder, along y-axis
314  ! External flow
315 
316  IF( radius(q) <= zero) THEN
317  IF(mype == pe_io) THEN
318  WRITE(*,*)'INPUT ERROR: CYLINDER:', q, &
319  ' HAS ZERO RADIUS.'
320  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
321  ENDIF
322  CALL mfix_exit(mype)
323  ELSE
324  lambda_x(q) = -one
325  lambda_y(q) = zero
326  lambda_z(q) = -one
327  dquadric(q) = radius(q)**2
328  ENDIF
329 
330  CASE ('Z_CYL_EXT') ! The quadric is predefined as a cylinder, along z-axis
331  ! External flow
332 
333  IF( radius(q) <= zero) THEN
334  IF(mype == pe_io) THEN
335  WRITE(*,*)'INPUT ERROR: CYLINDER:', q, &
336  ' HAS ZERO RADIUS.'
337  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
338  ENDIF
339  CALL mfix_exit(mype)
340  ELSE
341  lambda_x(q) = -one
342  lambda_y(q) = -one
343  lambda_z(q) = zero
344  dquadric(q) = radius(q)**2
345  ENDIF
346 
347  CASE ('SPHERE_INT') ! The quadric is predefined as a sphere
348  ! Internal flow
349 
350  IF( radius(q) <= zero) THEN
351  WRITE(*,*)'INPUT ERROR: SPHERE:', q, &
352  ' HAS INVALID RADIUS.'
353  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
354  CALL mfix_exit(mype)
355  ELSE
356  lambda_x(q) = one
357  lambda_y(q) = one
358  lambda_z(q) = one
359  dquadric(q) = -radius(q)**2
360  ENDIF
361 
362  CASE ('SPHERE_EXT') ! The quadric is predefined as a sphere
363  ! External flow
364 
365  IF( radius(q) <= zero) THEN
366  WRITE(*,*)'INPUT ERROR: SPHERE:', q, &
367  ' HAS INVALID RADIUS.'
368  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
369  CALL mfix_exit(mype)
370  ELSE
371  lambda_x(q) = -one
372  lambda_y(q) = -one
373  lambda_z(q) = -one
374  dquadric(q) = radius(q)**2
375  ENDIF
376 
377 
378  CASE ('X_CONE') ! The quadric is predefined as a cone, along x-axis
379  ! Internal flow
380 
381  IF(half_angle(q) <= zero .OR. half_angle(q) >= 90.0) THEN
382  IF(mype == pe_io) THEN
383  WRITE(*,*)'INPUT ERROR: CONE:', q, &
384  ' HAS INCORRECT HALF-ANGLE.'
385  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
386  ENDIF
387  CALL mfix_exit(mype)
388  ELSE
389  tan_half_angle = dtan(half_angle(q)/180.0*pi)
390  lambda_x(q) = -one
391  lambda_y(q) = one/(tan_half_angle)**2
392  lambda_z(q) = one/(tan_half_angle)**2
393  dquadric(q) = zero
394  ENDIF
395 
396  CASE ('Y_CONE') ! The quadric is predefined as a cone, along y-axis
397  ! Internal flow
398 
399  IF(half_angle(q) <= zero .OR. half_angle(q) >= 90.0) THEN
400  IF(mype == pe_io) THEN
401  WRITE(*,*)'INPUT ERROR: CONE:', q, &
402  ' HAS INCORRECT HALF-ANGLE.'
403  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
404  ENDIF
405  CALL mfix_exit(mype)
406  ELSE
407  tan_half_angle = dtan(half_angle(q)/180.0*pi)
408  lambda_x(q) = one/(tan_half_angle)**2
409  lambda_y(q) = -one
410  lambda_z(q) = one/(tan_half_angle)**2
411  dquadric(q) = zero
412  ENDIF
413 
414  CASE ('Z_CONE') ! The quadric is predefined as a cone, along z-axis
415  ! Internal flow
416 
417  IF(half_angle(q) <= zero .OR. half_angle(q) >= 90.0) THEN
418  IF(mype == pe_io) THEN
419  WRITE(*,*)'INPUT ERROR: CONE:', q, &
420  ' HAS INCORRECT HALF-ANGLE.'
421  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
422  ENDIF
423  CALL mfix_exit(mype)
424  ELSE
425  tan_half_angle = dtan(half_angle(q)/180.0*pi)
426  lambda_x(q) = one/(tan_half_angle)**2
427  lambda_y(q) = one/(tan_half_angle)**2
428  lambda_z(q) = -one
429  dquadric(q) = zero
430  ENDIF
431 
432  CASE ('C2C') ! Cylinder to cylinder junction using cone
433  ! Internal flow
434 
435  CALL build_cone_for_c2c(q)
436 
437 
438  CASE ('TORUS_INT','TORUS_EXT') ! Torus - Hard coded in define_quadrics.f
439  IF((torus_r1(q) <= zero).OR.(torus_r1(q)==undefined)) THEN
440  IF(mype == pe_io) THEN
441  WRITE(*,*)'INPUT ERROR: TORUS:', q, &
442  ' HAS INVALID RADIUS R1:',torus_r1(q)
443  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
444  ENDIF
445  CALL mfix_exit(mype)
446  ENDIF
447  IF((torus_r2(q) <= zero).OR.(torus_r2(q)==undefined)) THEN
448  IF(mype == pe_io) THEN
449  WRITE(*,*)'INPUT ERROR: TORUS:', q, &
450  ' HAS INVALID RADIUS R2:',torus_r2(q)
451  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
452  ENDIF
453  CALL mfix_exit(mype)
454  ENDIF
455 
456  CASE ('Y_UCOIL_EXT') ! UCOIL - Hard coded in define_quadrics.f
457  IF((ucoil_r1(q) <= zero).OR.(ucoil_r1(q)==undefined)) THEN
458  IF(mype == pe_io) THEN
459  WRITE(*,*)'INPUT ERROR: Y_UCOIL_EXT:', q, &
460  ' HAS INVALID RADIUS R1:',ucoil_r1(q)
461  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
462  ENDIF
463  CALL mfix_exit(mype)
464  ENDIF
465  IF((ucoil_r2(q) <= zero).OR.(ucoil_r2(q)==undefined)) THEN
466  IF(mype == pe_io) THEN
467  WRITE(*,*)'INPUT ERROR: Y_UCOIL_EXT:', q, &
468  ' HAS INVALID RADIUS R2:',ucoil_r2(q)
469  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
470  ENDIF
471  CALL mfix_exit(mype)
472  ENDIF
473 
474  IF(ucoil_y2(q)<ucoil_y1(q)) THEN
475  IF(mype == pe_io) THEN
476  WRITE(*,*)'INPUT ERROR: Y_UCOIL_EXT:', q, &
477  ' COIL_Y2 < COIL_Y1: Y2,Y1=',ucoil_y2(q),ucoil_y1(q)
478  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
479  ENDIF
480  CALL mfix_exit(mype)
481  ENDIF
482 
483  CASE ('Y_UCOIL2_EXT') ! UCOIL - Hard coded in define_quadrics.f
484  IF((ucoil_r1(q) <= zero).OR.(ucoil_r1(q)==undefined)) THEN
485  IF(mype == pe_io) THEN
486  WRITE(*,*)'INPUT ERROR: Y_UCOIL2_EXT:', q, &
487  ' HAS INVALID RADIUS R1:',ucoil_r1(q)
488  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
489  ENDIF
490  CALL mfix_exit(mype)
491  ENDIF
492  IF((ucoil_r2(q) <= zero).OR.(ucoil_r2(q)==undefined)) THEN
493  IF(mype == pe_io) THEN
494  WRITE(*,*)'INPUT ERROR: Y_UCOIL2_EXT:', q, &
495  ' HAS INVALID RADIUS R2:',ucoil_r2(q)
496  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
497  ENDIF
498  CALL mfix_exit(mype)
499  ENDIF
500 
501  IF(ucoil_y2(q)<ucoil_y1(q)) THEN
502  IF(mype == pe_io) THEN
503  WRITE(*,*)'INPUT ERROR: Y_UCOIL2_EXT:', q, &
504  ' COIL_Y2 < COIL_Y1: Y2,Y1=',ucoil_y2(q),ucoil_y1(q)
505  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
506  ENDIF
507  CALL mfix_exit(mype)
508  ENDIF
509 
510  CASE ('XY_BEND_INT') ! Bend - Hard coded in define_quadrics.f
511  IF((bend_r1(q) <= zero).OR.(bend_r1(q)==undefined)) THEN
512  IF(mype == pe_io) THEN
513  WRITE(*,*)'INPUT ERROR: XY_BEND_INT:', q, &
514  ' HAS INVALID RADIUS R1:',bend_r1(q)
515  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
516  ENDIF
517  CALL mfix_exit(mype)
518  ENDIF
519 
520  IF((bend_r2(q) <= zero).OR.(bend_r2(q)==undefined)) THEN
521  IF(mype == pe_io) THEN
522  WRITE(*,*)'INPUT ERROR: XY_BEND_INT:', q, &
523  ' HAS INVALID RADIUS R2:',bend_r2(q)
524  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
525  ENDIF
526  CALL mfix_exit(mype)
527  ENDIF
528 
529  IF((bend_theta1(q) < zero).OR.(bend_theta1(q)>360.0)) THEN
530  IF(mype == pe_io) THEN
531  WRITE(*,*)'INPUT ERROR: XY_BEND_INT:', q, &
532  ' HAS INVALID ANGLE THETA1:',bend_theta1(q)
533  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
534  ENDIF
535  CALL mfix_exit(mype)
536  ENDIF
537 
538  IF((bend_theta2(q) < zero).OR.(bend_theta2(q)>360.0)) THEN
539  IF(mype == pe_io) THEN
540  WRITE(*,*)'INPUT ERROR: XY_BEND_INT:', q, &
541  ' HAS INVALID ANGLE THETA2:',bend_theta2(q)
542  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
543  ENDIF
544  CALL mfix_exit(mype)
545  ENDIF
546 
547  CASE ('Y_C2C_INT') ! Cylinder-cone-cylinder - Hard coded in define_quadrics.f
548  IF((c2c_r1(q) <= zero).OR.(c2c_r1(q)==undefined)) THEN
549  IF(mype == pe_io) THEN
550  WRITE(*,*)'INPUT ERROR: Y_C2C_INT:', q, &
551  ' HAS INVALID RADIUS R1:',c2c_r1(q)
552  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
553  ENDIF
554  CALL mfix_exit(mype)
555  ENDIF
556 
557  IF((c2c_r2(q) <= zero).OR.(c2c_r2(q)==undefined)) THEN
558  IF(mype == pe_io) THEN
559  WRITE(*,*)'INPUT ERROR: C2C_XY_INT:', q, &
560  ' HAS INVALID RADIUS R2:',c2c_r2(q)
561  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
562  ENDIF
563  CALL mfix_exit(mype)
564  ENDIF
565 
566  IF(c2c_y2(q) < c2c_y1(q)) THEN
567  IF(mype == pe_io) THEN
568  WRITE(*,*)'INPUT ERROR: Y_C2C_INT:', q
569  WRITE(*,*)'MUST HAVE C2C_Y2 >= C2C_Y1.'
570  WRITE(*,*)'C2C_Y1,C2C_Y2 =', c2c_y1(q),c2c_y2(q)
571  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
572  ENDIF
573  CALL mfix_exit(mype)
574  ENDIF
575 
576  IF((c2c_y1(q) == c2c_y2(q)).AND.(c2c_r1(q)/=c2c_r2(q))) THEN
577  IF(mype == pe_io) THEN
578  WRITE(*,*)'INPUT ERROR: Y_C2C_INT:', q, &
579  ' C2C_Y1=C2C_Y2 BUT C2C_R1/=C2C_R2:'
580  WRITE(*,*)'C2C_Y1,C2C_Y2 =', c2c_y1(q),c2c_y2(q)
581  WRITE(*,*)'C2C_R1,C2C_R2 =', c2c_r1(q),c2c_r2(q)
582  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
583  ENDIF
584  CALL mfix_exit(mype)
585  ENDIF
586 
587  CASE ('REACTOR1') ! Cylinder-cone-cylinder - Hard coded in define_quadrics.f
588 
589  IF(reactor1_y2(q) < reactor1_y1(q)) THEN
590  IF(mype == pe_io) THEN
591  WRITE(*,*)'INPUT ERROR: REACTOR1:', q
592  WRITE(*,*)'MUST HAVE REACTOR1_Y2 >= REACTOR1_Y1.'
593  WRITE(*,*)'REACTOR1_Y1,REACTOR1_Y2 =', reactor1_y1(q),reactor1_y2(q)
594  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
595  ENDIF
596  CALL mfix_exit(mype)
597  ENDIF
598 
599  IF((reactor1_y1(q) == reactor1_y2(q)).AND.(reactor1_r1(q)/=reactor1_r2(q))) THEN
600  IF(mype == pe_io) THEN
601  WRITE(*,*)'INPUT ERROR: REACTOR1:', q, &
602  ' REACTOR1_Y1=REACTOR1_Y2 BUT REACTOR1_R1/=REACTOR1_R2:'
603  WRITE(*,*)'REACTOR1_Y1,REACTOR1_Y2 =', reactor1_y1(q),reactor1_y2(q)
604  WRITE(*,*)'REACTOR1_R1,REACTOR1_R2 =', reactor1_r1(q),reactor1_r2(q)
605  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
606  ENDIF
607  CALL mfix_exit(mype)
608  ENDIF
609 
610 
611  IF(reactor1_yr2(q) <= reactor1_y2(q)) THEN
612  IF(mype == pe_io) THEN
613  WRITE(*,*)'INPUT ERROR: REACTOR1:', q
614  WRITE(*,*)'MUST HAVE REACTOR1_YR2 > REACTOR1_Y2.'
615  WRITE(*,*)'REACTOR1_YR2,REACTOR1_Y2 =', reactor1_yr2(q),reactor1_y2(q)
616  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
617  ENDIF
618  CALL mfix_exit(mype)
619  ENDIF
620 
621  IF(reactor1_yr1(q) >= reactor1_y1(q)) THEN
622  IF(mype == pe_io) THEN
623  WRITE(*,*)'INPUT ERROR: REACTOR1:', q
624  WRITE(*,*)'MUST HAVE REACTOR1_YR1 < REACTOR1_Y1.'
625  WRITE(*,*)'REACTOR1_YR1,REACTOR1_Y1 =', reactor1_yr1(q),reactor1_y1(q)
626  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
627  ENDIF
628  CALL mfix_exit(mype)
629  ENDIF
630 
631  IF(reactor1_theta1(q) <= zero.OR.reactor1_theta1(q) > 90.0d0) THEN
632  IF(mype == pe_io) THEN
633  WRITE(*,*)'INPUT ERROR: REACTOR1:', q
634  WRITE(*,*)'MUST HAVE 0.0 < REACTOR1_THETA1 <= 90 DEGREES.'
635  WRITE(*,*)'REACTOR1_THETA1 =', reactor1_theta1(q)
636  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
637  ENDIF
638  CALL mfix_exit(mype)
639  ENDIF
640 
641  IF(reactor1_theta2(q) <= zero.OR.reactor1_theta2(q) > 90.0d0) THEN
642  IF(mype == pe_io) THEN
643  WRITE(*,*)'INPUT ERROR: REACTOR1:', q
644  WRITE(*,*)'MUST HAVE 0.0 < REACTOR1_THETA2 <= 90 DEGREES.'
645  WRITE(*,*)'REACTOR1_THETA2 =', reactor1_theta2(q)
646  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
647  ENDIF
648  CALL mfix_exit(mype)
649  ENDIF
650 ! Convert angles from degrees to radians
651  reactor1_theta1(q) = reactor1_theta1(q)/180.0d0*pi
652  reactor1_theta2(q) = reactor1_theta2(q)/180.0d0*pi
653 
654 
655  CASE DEFAULT
656  IF(mype == pe_io) THEN
657  WRITE(*,*)'INPUT ERROR: QUADRIC:', q, &
658  ' HAS INCORRECT FORM: ',quadric_form(q)
659  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
660  ENDIF
661  CALL mfix_exit(mype)
662 
663  END SELECT
664 
665  IF(bc_id_q(q) == undefined_i) THEN
666  IF(mype == pe_io) THEN
667  WRITE(*,*)'INPUT ERROR: QUADRIC:', q, &
668  ' HAS NO ASSIGNED BC ID.'
669  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
670  ENDIF
671  CALL mfix_exit(mype)
672  ENDIF
673 
674  ENDDO
675 
676 
677  IF(n_quadric>0) THEN
678 
679 
680  IF(n_group > dim_group) THEN
681  IF(mype == pe_io) THEN
682  WRITE(*,*)'INPUT ERROR: INVALID VALUE OF N_GROUP =', n_group
683  WRITE(*,*)'MAXIMUM ACCEPTABLE VALUE IS DIM_GROUP =', dim_group
684  WRITE(*,*)'CHANGE MAXIMUM VALUE IN QUADRIC_MOD.F IF NECESSARY.'
685  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
686  ENDIF
687  CALL mfix_exit(mype)
688  ENDIF
689 
690 
691  DO i = 1,n_group
692 
693  IF(group_size(i) < 1 .OR. group_size(i) > n_quadric) THEN
694  IF(mype == pe_io) THEN
695  WRITE(*,*)'INPUT ERROR: GROUP:', i, ' HAS INCORRECT SIZE:', group_size(i)
696  WRITE(*,*)'VALID GROUP SIZE RANGE IS:', 1, ' TO ', n_quadric
697  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
698  ENDIF
699  CALL mfix_exit(mype)
700  ENDIF
701 
702  DO j = 1,group_size(i)
703  IF(group_q(i,j) < 1 .OR. group_q(i,j) > n_quadric) THEN
704  IF(mype == pe_io) THEN
705  WRITE(*,*)'INPUT ERROR: GROUP_Q(', i,',',j, ') HAS INCORRECT VALUE:', group_q(i,j)
706  WRITE(*,*)'VALID GROUP_Q RANGE IS:', 1, ' TO ', n_quadric
707  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
708  ENDIF
709  CALL mfix_exit(mype)
710  ENDIF
711  ENDDO
712 
713  gr = trim(group_relation(i))
714 
715  IF(gr/='OR'.AND.gr/='AND'.AND.gr/='PIECEWISE') THEN
716  IF(mype == pe_io) THEN
717  WRITE(*,*)'INPUT ERROR: GROUP:', i, ' HAS INCORRECT GROUP RELATION: ', gr
718  WRITE(*,*)'VALID GROUP RELATIONS ARE ''OR'',''AND'', AND ''PIECEWISE''. '
719  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
720  ENDIF
721  CALL mfix_exit(mype)
722  ENDIF
723 
724  ENDDO
725 
726  DO i = 2,n_group
727 
728  gr = trim(relation_with_previous(i))
729 
730  IF(gr/='OR'.AND.gr/='AND') THEN
731  IF(mype == pe_io) THEN
732  WRITE(*,*)'INPUT ERROR: GROUP:', i, ' HAS INCORRECT RELATION WITH PREVIOUS: ', gr
733  WRITE(*,*)'VALID GROUP RELATIONS ARE ''OR'', AND ''AND''. '
734  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
735  ENDIF
736  CALL mfix_exit(mype)
737  ENDIF
738 
739  ENDDO
740 
741  ENDIF
742 
743 
744  IF(tol_snap(1)<zero.OR.tol_snap(1)>half) THEN
745  IF(mype == pe_io) THEN
746  WRITE(*,*)'INPUT ERROR: INVALID VALUE OF TOL_SNAP IN X-DIRECTION =', tol_snap(1)
747  WRITE(*,*)'ACCEPTABLE VALUES ARE BETWEEN 0.0 AND 0.5.'
748  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
749  ENDIF
750  CALL mfix_exit(mype)
751  ENDIF
752 
753  IF(tol_snap(2)==undefined) tol_snap(2)=tol_snap(1)
754 
755  IF(tol_snap(2)<zero.OR.tol_snap(2)>half) THEN
756  IF(mype == pe_io) THEN
757  WRITE(*,*)'INPUT ERROR: INVALID VALUE OF TOL_SNAP IN Y-DIRECTION =', tol_snap(2)
758  WRITE(*,*)'ACCEPTABLE VALUES ARE BETWEEN 0.0 AND 0.5.'
759  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
760  ENDIF
761  CALL mfix_exit(mype)
762  ENDIF
763 
764  IF(tol_snap(3)==undefined) tol_snap(3)=tol_snap(1)
765 
766  IF(tol_snap(3)<zero.OR.tol_snap(3)>half) THEN
767  IF(mype == pe_io) THEN
768  WRITE(*,*)'INPUT ERROR: INVALID VALUE OF TOL_SNAP IN Z-DIRECTION =', tol_snap(3)
769  WRITE(*,*)'ACCEPTABLE VALUES ARE BETWEEN 0.0 AND 0.5.'
770  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
771  ENDIF
772  CALL mfix_exit(mype)
773  ENDIF
774 
775 
776  IF(tol_delh<zero.OR.tol_delh>one) THEN
777  IF(mype == pe_io) THEN
778  WRITE(*,*)'INPUT ERROR: INVALID VALUE OF TOL_DELH =', tol_delh
779  WRITE(*,*)'ACCEPTABLE VALUES ARE BETWEEN 0.0 AND 1.0.'
780  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
781  ENDIF
782  CALL mfix_exit(mype)
783  ENDIF
784 
785  IF(tol_small_cell<zero.OR.tol_small_cell>one) THEN
786  IF(mype == pe_io) THEN
787  WRITE(*,*)'INPUT ERROR: INVALID VALUE OF TOL_SMALL_CELL =', tol_small_cell
788  WRITE(*,*)'ACCEPTABLE VALUES ARE BETWEEN 0.0 AND 1.0.'
789  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
790  ENDIF
791  CALL mfix_exit(mype)
792  ENDIF
793 
794  IF(tol_small_area<zero.OR.tol_small_area>one) THEN
795  IF(mype == pe_io) THEN
796  WRITE(*,*)'INPUT ERROR: INVALID VALUE OF TOL_SMALL_AREA =', tol_small_area
797  WRITE(*,*)'ACCEPTABLE VALUES ARE BETWEEN 0.0 AND 1.0.'
798  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
799  ENDIF
800  CALL mfix_exit(mype)
801  ENDIF
802 
803  IF(alpha_max<zero) THEN
804  IF(mype == pe_io) THEN
805  WRITE(*,*)'INPUT ERROR: NEGATIVE VALUE OF ALPHA_MAX =', alpha_max
806  WRITE(*,*)'ACCEPTABLE VALUES ARE POSITIVE NUMBERS (E.G. 1.0).'
807  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
808  ENDIF
809  CALL mfix_exit(mype)
810  ENDIF
811 
812 
813  IF(tol_f<zero) THEN
814  IF(mype == pe_io) THEN
815  WRITE(*,*)'INPUT ERROR: NEGATIVE VALUE OF TOL_F =', tol_f
816  WRITE(*,*)'ACCEPTABLE VALUES ARE SMALL POSITIVE NUMBERS (E.G. 1.0E-9).'
817  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
818  ENDIF
819  CALL mfix_exit(mype)
820  ENDIF
821 
822  IF(tol_poly<zero) THEN
823  IF(mype == pe_io) THEN
824  WRITE(*,*)'INPUT ERROR: NEGATIVE VALUE OF TOL_POLY =', tol_poly
825  WRITE(*,*)'ACCEPTABLE VALUES ARE SMALL POSITIVE NUMBERS (E.G. 1.0E-9).'
826  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
827  ENDIF
828  CALL mfix_exit(mype)
829  ENDIF
830 
831  IF(itermax_int<0) THEN
832  IF(mype == pe_io) THEN
833  WRITE(*,*)'INPUT ERROR: NEGATIVE VALUE OF ITERMAX_INT =', itermax_int
834  WRITE(*,*)'ACCEPTABLE VALUES ARE LARGE POSITIVE INTEGERS (E.G. 10000).'
835  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
836  ENDIF
837  CALL mfix_exit(mype)
838  ENDIF
839 
840  IF(fac_dim_max_cut_cell<0.05.OR.fac_dim_max_cut_cell>5.0) THEN
841  IF(mype == pe_io) THEN
842  WRITE(*,*)'INPUT ERROR: NEGATIVE VALUE OF FAC_DIM_MAX_CUT_CELL =', fac_dim_max_cut_cell
843  WRITE(*,*)'ACCEPTABLE VALUES ARE BETWEEN 0.05 AND 5.0.'
844  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
845  ENDIF
846  CALL mfix_exit(mype)
847  ENDIF
848 
849 
850  IF((cg_safe_mode(1)==1).AND.(pg_option/=0)) THEN
851  pg_option = 0
852  IF(mype == pe_io) WRITE(*,*)'WARNING: SAFE_MODE ACTIVATED FOR GAS PRESSURE, REVERTING TO PG_OPTION = 0'
853  ENDIF
854 
855  IF(pg_option <0 .OR. pg_option>2) THEN
856  IF(mype == pe_io) THEN
857  WRITE(*,*)'INPUT ERROR: INVALID VALUE OF PG_OPTION =', pg_option
858  WRITE(*,*)'ACCEPTABLE VALUES ARE 0,1,AND 2.'
859  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
860  ENDIF
861  CALL mfix_exit(mype)
862  ENDIF
863 
864  IF(cg_ur_fac(2)<zero.OR.cg_ur_fac(2)>one) THEN
865  IF(mype == pe_io) THEN
866  WRITE(*,*)'INPUT ERROR: INVALID VALUE OF CG_UR_FAC(2) =', cg_ur_fac(2)
867  WRITE(*,*)'ACCEPTABLE VALUES ARE BETWEEN 0.0 AND 1.0.'
868  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
869  ENDIF
870  CALL mfix_exit(mype)
871  ENDIF
872 
873  IF(bar_width<10.OR.bar_width>80) THEN
874  IF(mype == pe_io) THEN
875  WRITE(*,*)'INPUT ERROR: INVALID VALUE OF BAR_WIDTH =', bar_width
876  WRITE(*,*)'ACCEPTABLE VALUES ARE BETWEEN 10 AND 80.'
877  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
878  ENDIF
879  CALL mfix_exit(mype)
880  ENDIF
881 
882  IF(bar_resolution<one.OR.bar_resolution>100.0) THEN
883  IF(mype == pe_io) THEN
884  WRITE(*,*)'INPUT ERROR: INVALID VALUE OF BAR_RESOLUTION =', bar_resolution
885  WRITE(*,*)'ACCEPTABLE VALUES ARE BETWEEN 0.0 AND 100.0.'
886  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
887  ENDIF
888  CALL mfix_exit(mype)
889  ENDIF
890 
891  IF(f_dashboard<1) THEN
892  IF(mype == pe_io) THEN
893  WRITE(*,*)'INPUT ERROR: INVALID VALUE OF F_DASHBOARD =', f_dashboard
894  WRITE(*,*)'ACCEPTABLE VALUES ARE INTEGERS >= 1.'
895  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
896  ENDIF
897  CALL mfix_exit(mype)
898  ENDIF
899 
900 
901 
902 !======================================================================
903 ! Data initialization for Dashboard
904 !======================================================================
905  init_time = time
906  smmin = large_number
907  smmax = -large_number
908 
909  dtmin = large_number
910  dtmax = -large_number
911 
912  nit_min = max_nit
913  nit_max = 0
914 
915  n_dashboard = 0
916 
917 
918 
919  cg_mi_converted_to_ps = .false.
920 
921  DO bcv = 1, dimension_bc
922 
923  IF (bc_type_enum(bcv) == cg_mi) THEN
924  bc_type_enum(bcv) = cg_nsw
925  cg_mi_converted_to_ps(bcv) = .true.
926  if(mype==0) print*,'From check_data_cartesian: Converted CG_MI to CG_NSW for BC#',bcv
927  ENDIF
928 
929  ENDDO
930 
931 
932  IF(re_indexing) THEN
933 
934  IF(mype==0) THEN
935  WRITE(*,*)' From check_data_cartesian: RE_INDEXING is turned on.'
936  WRITE(*,*)' The preconditionner will be turned off for all equations'
937  WRITE(*,*)' regardless of the mfix.dat setting.'
938  leq_pc = 'NONE'
939  ENDIF
940  ENDIF
941 
942  RETURN
943  END SUBROUTINE check_data_cartesian
944 
945 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
946 ! C
947 ! Module name: CHECK_BC_FLAGS C
948 ! Purpose: check the boundary conditions flags C
949 ! C
950 ! Author: Jeff Dietiker Date: 21-Feb-08 C
951 ! Reviewer: Date: C
952 ! C
953 ! Revision Number # Date: ##-###-## C
954 ! Author: # C
955 ! Purpose: # C
956 ! C
957 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
958 
959  SUBROUTINE check_bc_flags
960 
961 !-----------------------------------------------
962 ! Modules
963 !-----------------------------------------------
964  USE bc
965  USE compar
966  USE constant
967  USE cutcell
968  USE dashboard
969  USE fldvar
970  USE functions
971  USE funits
972  USE indices
973  USE leqsol
974  USE mpi_utility
975  USE param
976  USE param1
977  USE physprop
978  USE polygon
979  USE quadric
980  USE run
981  USE scalars
982  USE toleranc
983  USE vtk
984  use error_manager
985 
986  IMPLICIT NONE
987 !-----------------------------------------------
988 ! Local variables
989 !-----------------------------------------------
990  INTEGER :: IJK,IJKW,IJKS,IJKB,M,NN
991  INTEGER :: IJKWW,IJKSS,IJKBB
992  INTEGER :: BCV,BCV_U,BCV_V,BCV_W
993 !-----------------------------------------------
994  DOUBLE PRECISION SUM, SUM_EP
995 !-----------------------------------------------
996 !======================================================================
997 ! Boundary conditions
998 !======================================================================
999 
1000  CALL init_err_msg("CHECK_BC_FLAGS")
1001 
1002  DO bcv = 1, dimension_bc
1003  IF(cg_mi_converted_to_ps(bcv)) THEN
1004  IF(bc_massflow_g(bcv)==undefined) THEN
1005  WRITE(err_msg, 1710) bcv
1006  CALL flush_err_msg(abort=.true.)
1007  ENDIF
1008  DO m = 1, mmax
1009  IF(bc_massflow_s(bcv,m)==undefined) THEN
1010  WRITE(err_msg, 1711) bcv, m
1011  CALL flush_err_msg(abort=.true.)
1012  ENDIF
1013  ENDDO
1014  ENDIF
1015  ENDDO
1016 
1017 1710 FORMAT('Error 1110: BC :',i3,'. When using CG_MI, the gas mass flow rate',/1x, &
1018  'must be specified, including when it is zero.',/1x, &
1019  ' Please correct the mfix.dat file.')
1020 
1021 1711 FORMAT('Error 1111: BC :',i3,'. When using CG_MI, the solids mass flow rate',/1x, &
1022  'for M=',i4,' must be specified, including when it is zero.',/1x, &
1023  ' Please correct the mfix.dat file.')
1024 
1025  DO ijk = ijkstart3, ijkend3
1026  bcv = bc_id(ijk)
1027  IF(bcv>0) THEN
1028 
1029  IF(bc_type_enum(bcv) == cg_mi) THEN
1030 
1031  print*,'CG_MI at', ijk ! This should not be printed on the screen anymore after conversion to point source.
1032 
1033 ! FLAG(IJK) = 20
1034 ! FLAG_E(IJK) = UNDEFINED_I
1035 ! FLAG_N(IJK) = UNDEFINED_I
1036 ! FLAG_T(IJK) = UNDEFINED_I
1037 
1038  ELSEIF(bc_type_enum(bcv) == cg_po) THEN
1039 
1040  flag(ijk) = 11
1041  flag_e(ijk) = undefined_i
1042  flag_n(ijk) = undefined_i
1043  flag_t(ijk) = undefined_i
1044 
1045  ijkw = west_of(ijk)
1046  bcv_u = bc_u_id(ijkw)
1047  IF(bcv_u>0) THEN
1048  IF(bc_type_enum(bcv_u) == cg_po) THEN
1049  flag(ijkw) = 11
1050  flag_e(ijkw) = undefined_i
1051  flag_n(ijkw) = undefined_i
1052  flag_t(ijkw) = undefined_i
1053  ENDIF
1054  ENDIF
1055 
1056  ijks = south_of(ijk)
1057  bcv_v = bc_v_id(ijks)
1058  IF(bcv_v>0) THEN
1059  IF(bc_type_enum(bcv_v) == cg_po) THEN
1060  flag(ijks) = 11
1061  flag_e(ijks) = undefined_i
1062  flag_n(ijks) = undefined_i
1063  flag_t(ijks) = undefined_i
1064  ENDIF
1065  ENDIF
1066 
1067  IF(do_k) THEN
1068  ijkb = bottom_of(ijk)
1069  bcv_w = bc_w_id(ijkb)
1070  IF(bcv_w>0) THEN
1071  IF(bc_type_enum(bcv_w) == cg_po) THEN
1072  flag(ijkb) = 11
1073  flag_e(ijkb) = undefined_i
1074  flag_n(ijkb) = undefined_i
1075  flag_t(ijkb) = undefined_i
1076  ENDIF
1077  ENDIF
1078  ENDIF
1079 
1080  ENDIF
1081  ENDIF
1082  ENDDO
1083 
1084 
1085  DO ijk = ijkstart3, ijkend3
1086  bcv = bc_id(ijk)
1087  IF(bcv>0) THEN
1088  IF(bc_type_enum(bcv) == cg_mi) THEN
1089 
1090 ! IJKW = WEST_OF(IJK)
1091 ! IF(FLUID_AT(IJKW)) THEN
1092 ! FLAG_E(IJKW) = 2020
1093 ! ENDIF
1094 
1095 ! IJKS = SOUTH_OF(IJK)
1096 ! IF(FLUID_AT(IJKS)) THEN
1097 ! FLAG_N(IJKS) = 2020
1098 ! ENDIF
1099 
1100 ! IJKB = BOTTOM_OF(IJK)
1101 ! IF(FLUID_AT(IJKB)) THEN
1102 ! FLAG_T(IJKB) = 2020
1103 ! ENDIF
1104 
1105  IF (bc_u_g(bcv) == undefined) THEN
1106  IF (no_i) THEN
1107  bc_u_g(bcv) = zero
1108  ELSEIF(bc_volflow_g(bcv)==undefined.AND. &
1109  bc_massflow_g(bcv)==undefined.AND.&
1110  bc_velmag_g(bcv)==undefined) THEN
1111  IF(dmp_log)WRITE (unit_log, 900) 'BC_U_g', bcv
1112  call mfix_exit(mype)
1113  ENDIF
1114  ENDIF
1115  IF (bc_v_g(bcv) == undefined) THEN
1116  IF (no_j) THEN
1117  bc_v_g(bcv) = zero
1118  ELSEIF(bc_volflow_g(bcv)==undefined.AND. &
1119  bc_massflow_g(bcv)==undefined.AND.&
1120  bc_velmag_g(bcv)==undefined) THEN
1121  IF(dmp_log)WRITE (unit_log, 900) 'BC_V_g', bcv
1122  call mfix_exit(mype)
1123  ENDIF
1124  ENDIF
1125  IF (bc_w_g(bcv) == undefined) THEN
1126  IF (no_k) THEN
1127  bc_w_g(bcv) = zero
1128  ELSEIF(bc_volflow_g(bcv)==undefined.AND. &
1129  bc_massflow_g(bcv)==undefined.AND.&
1130  bc_velmag_g(bcv)==undefined) THEN
1131  IF(dmp_log)WRITE (unit_log, 900) 'BC_W_g', bcv
1132  call mfix_exit(mype)
1133  ENDIF
1134  ENDIF
1135  IF (k_epsilon .AND. bc_k_turb_g(bcv) == undefined) THEN
1136  IF(dmp_log)WRITE (unit_log, 1000) 'BC_K_Turb_G', bcv
1137  call mfix_exit(mype)
1138  ENDIF
1139  IF (k_epsilon .AND. bc_e_turb_g(bcv) == undefined) THEN
1140  IF(dmp_log)WRITE (unit_log, 1000) 'BC_E_Turb_G', bcv
1141  call mfix_exit(mype)
1142  ENDIF
1143 
1144 ! Check whether the bc velocity components have the correct sign
1145 ! SELECT CASE (BC_PLANE(BCV))
1146 ! CASE ('W')
1147 ! IF (BC_U_G(BCV) > ZERO) THEN
1148 ! IF(DMP_LOG)WRITE (UNIT_LOG, 1050) BCV, 'BC_U_g', '<'
1149 ! CALL MFIX_EXIT(myPE)
1150 ! ENDIF
1151 ! CASE ('E')
1152 ! IF (BC_U_G(BCV) < ZERO) THEN
1153 ! IF(DMP_LOG)WRITE (UNIT_LOG, 1050) BCV, 'BC_U_g', '>'
1154 ! CALL MFIX_EXIT(myPE)
1155 ! ENDIF
1156 ! CASE ('S')
1157 ! IF (BC_V_G(BCV) > ZERO) THEN
1158 ! IF(DMP_LOG)WRITE (UNIT_LOG, 1050) BCV, 'BC_V_g', '<'
1159 ! CALL MFIX_EXIT(myPE)
1160 ! ENDIF
1161 ! CASE ('N')
1162 ! IF (BC_V_G(BCV) < ZERO) THEN
1163 ! IF(DMP_LOG)WRITE (UNIT_LOG, 1050) BCV, 'BC_V_g', '>'
1164 ! CALL MFIX_EXIT(myPE)
1165 ! ENDIF
1166 ! CASE ('B')
1167 ! IF (BC_W_G(BCV) > ZERO) THEN
1168 ! IF(DMP_LOG)WRITE (UNIT_LOG, 1050) BCV, 'BC_W_g', '<'
1169 ! CALL MFIX_EXIT(myPE)
1170 ! ENDIF
1171 ! CASE ('T')
1172 ! IF (BC_W_G(BCV) < ZERO) THEN
1173 ! IF(DMP_LOG)WRITE (UNIT_LOG, 1050) BCV, 'BC_W_g', '>'
1174 ! CALL MFIX_EXIT(myPE)
1175 ! ENDIF
1176 ! END SELECT
1177 
1178  sum_ep = bc_ep_g(bcv)
1179  DO m = 1, mmax
1180  IF (bc_rop_s(bcv,m) == undefined) THEN
1181  IF (bc_ep_g(bcv) == one) THEN
1182  bc_rop_s(bcv,m) = zero
1183  ELSEIF (mmax == 1) THEN
1184  bc_rop_s(bcv,m) = (one - bc_ep_g(bcv))*ro_s0(m)
1185  ELSE
1186  IF(dmp_log)WRITE (unit_log, 1100) 'BC_ROP_s', bcv, m
1187  call mfix_exit(mype)
1188  ENDIF
1189  ENDIF
1190 
1191  sum_ep = sum_ep + bc_rop_s(bcv,m)/ro_s0(m)
1192  IF (species_eq(m)) THEN
1193  sum = zero
1194  DO nn = 1, nmax(m)
1195  IF(bc_x_s(bcv,m,nn)/=undefined)sum=sum+bc_x_s(bcv,m,nn)
1196  ENDDO
1197 
1198  IF (bc_rop_s(bcv,m)==zero .AND. sum==zero) THEN
1199  bc_x_s(bcv,m,1) = one
1200  sum = one
1201  ENDIF
1202 
1203  DO nn = 1, nmax(m)
1204  IF (bc_x_s(bcv,m,nn) == undefined) THEN
1205  IF(.NOT.compare(one,sum) .AND. dmp_log)WRITE (unit_log,1110)bcv,m,nn
1206  bc_x_s(bcv,m,nn) = zero
1207  ENDIF
1208  ENDDO
1209 
1210  IF (.NOT.compare(one,sum)) THEN
1211  IF(dmp_log)WRITE (unit_log, 1120) bcv, m
1212  call mfix_exit(mype)
1213  ENDIF
1214  ENDIF
1215 
1216  IF (bc_u_s(bcv,m) == undefined) THEN
1217  IF (bc_rop_s(bcv,m)==zero .OR. no_i) THEN
1218  bc_u_s(bcv,m) = zero
1219  ELSEIF(bc_volflow_s(bcv,m)==undefined.AND. &
1220  bc_massflow_s(bcv,m)==undefined.AND.&
1221  bc_velmag_s(bcv,m)==undefined) THEN
1222  IF(dmp_log)WRITE (unit_log, 910) 'BC_U_s', bcv, m
1223  call mfix_exit(mype)
1224  ENDIF
1225  ENDIF
1226 
1227  IF (bc_v_s(bcv,m) == undefined) THEN
1228  IF (bc_rop_s(bcv,m)==zero .OR. no_j) THEN
1229  bc_v_s(bcv,m) = zero
1230  ELSEIF(bc_volflow_s(bcv,m)==undefined.AND. &
1231  bc_massflow_s(bcv,m)==undefined.AND.&
1232  bc_velmag_s(bcv,m)==undefined) THEN
1233  IF(dmp_log)WRITE (unit_log, 910) 'BC_V_s', bcv, m
1234  call mfix_exit(mype)
1235  ENDIF
1236  ENDIF
1237 
1238  IF (bc_w_s(bcv,m) == undefined) THEN
1239  IF (bc_rop_s(bcv,m)==zero .OR. no_k) THEN
1240  bc_w_s(bcv,m) = zero
1241  ELSEIF(bc_volflow_s(bcv,m)==undefined.AND. &
1242  bc_massflow_s(bcv,m)==undefined.AND.&
1243  bc_velmag_s(bcv,m)==undefined) THEN
1244  IF(dmp_log)WRITE (unit_log, 910) 'BC_W_s', bcv, m
1245  call mfix_exit(mype)
1246  ENDIF
1247  ENDIF
1248 
1249  IF (energy_eq .AND. bc_t_s(bcv,m)==undefined) THEN
1250  IF (bc_rop_s(bcv,m) == zero) THEN
1251  bc_t_s(bcv,m) = bc_t_g(bcv)
1252  ELSE
1253  IF(dmp_log)WRITE (unit_log, 1100) 'BC_T_s', bcv, m
1254  call mfix_exit(mype)
1255  ENDIF
1256  ENDIF
1257 
1258  IF (granular_energy .AND. bc_theta_m(bcv,m)==undefined) THEN
1259  IF (bc_rop_s(bcv,m) == zero) THEN
1260  bc_theta_m(bcv,m) = zero
1261  ELSE
1262  IF(dmp_log)WRITE (unit_log, 1100) 'BC_Theta_m', bcv, m
1263  call mfix_exit(mype)
1264  ENDIF
1265  ENDIF
1266 
1267 ! Check whether the bc velocity components have the correct sign
1268 ! SELECT CASE (TRIM(BC_PLANE(BCV)))
1269 ! CASE ('W')
1270 ! IF (BC_U_S(BCV,M) > ZERO) THEN
1271 ! IF(DMP_LOG)WRITE (UNIT_LOG, 1150) BCV, 'BC_U_s', M, '<'
1272 ! CALL MFIX_EXIT(myPE)
1273 ! ENDIF
1274 ! CASE ('E')
1275 ! IF (BC_U_S(BCV,M) < ZERO) THEN
1276 ! IF(DMP_LOG)WRITE (UNIT_LOG, 1150) BCV, 'BC_U_s', M, '>'
1277 ! CALL MFIX_EXIT(myPE)
1278 ! ENDIF
1279 ! CASE ('S')
1280 ! IF (BC_V_S(BCV,M) > ZERO) THEN
1281 ! IF(DMP_LOG)WRITE (UNIT_LOG, 1150) BCV, 'BC_V_s', M, '<'
1282 ! CALL MFIX_EXIT(myPE)
1283 ! ENDIF
1284 ! CASE ('N')
1285 ! IF (BC_V_S(BCV,M) < ZERO) THEN
1286 ! IF(DMP_LOG)WRITE (UNIT_LOG, 1150) BCV, 'BC_V_s', M, '>'
1287 ! CALL MFIX_EXIT(myPE)
1288 ! ENDIF
1289 ! CASE ('B')
1290 ! IF (BC_W_S(BCV,M) > ZERO) THEN
1291 ! IF(DMP_LOG)WRITE (UNIT_LOG, 1150) BCV, 'BC_W_s', M, '<'
1292 ! CALL MFIX_EXIT(myPE)
1293 ! ENDIF
1294 ! CASE ('T')
1295 ! IF (BC_W_S(BCV,M) < ZERO) THEN
1296 ! IF(DMP_LOG)WRITE (UNIT_LOG, 1150) BCV, 'BC_W_s', M, '>'
1297 ! CALL MFIX_EXIT(myPE)
1298 ! ENDIF
1299 ! END SELECT
1300 
1301 
1302  ENDDO
1303 
1304  IF (.NOT.compare(one,sum_ep)) THEN
1305  IF(dmp_log)WRITE (unit_log, 1125) bcv
1306  call mfix_exit(mype)
1307  ENDIF
1308 
1309  DO nn = 1, nscalar
1310  IF (bc_scalar(bcv,nn) == undefined) THEN
1311  IF(dmp_log)WRITE (unit_log, 1004) 'BC_Scalar', bcv, nn
1312  CALL mfix_exit(mype)
1313  ENDIF
1314  ENDDO
1315 
1316 
1317  ELSEIF(bc_type_enum(bcv) == cg_po) THEN
1318 
1319  ijkw = west_of(ijk)
1320  IF(fluid_at(ijkw)) THEN
1321  flag_e(ijkw) = 2011
1322  ENDIF
1323 
1324  bcv_u = bc_u_id(ijkw)
1325  IF(bcv_u>0) THEN
1326  IF(bc_type_enum(bcv_u) == cg_po) THEN
1327  ijkww = west_of(ijkw)
1328  IF(fluid_at(ijkww)) THEN
1329  flag_e(ijkww) = 2011
1330  ENDIF
1331  ENDIF
1332  ENDIF
1333 
1334  ijks = south_of(ijk)
1335  IF(fluid_at(ijks)) THEN
1336  flag_n(ijks) = 2011
1337  ENDIF
1338 
1339  bcv_v = bc_v_id(ijks)
1340  IF(bcv_v>0) THEN
1341  IF(bc_type_enum(bcv_v) == cg_po) THEN
1342  ijkss = south_of(ijks)
1343  IF(fluid_at(ijkss)) THEN
1344  flag_n(ijkss) = 2011
1345  ENDIF
1346  ENDIF
1347  ENDIF
1348 
1349 
1350  IF(do_k) THEN
1351  ijkb = bottom_of(ijk)
1352  IF(fluid_at(ijkb)) THEN
1353  flag_t(ijkb) = 2011
1354  ENDIF
1355 
1356  bcv_w = bc_w_id(ijkb)
1357  IF(bcv_w>0) THEN
1358  IF(bc_type_enum(bcv_w) == cg_po) THEN
1359  ijkbb = bottom_of(ijkb)
1360  IF(fluid_at(ijkbb)) THEN
1361  flag_t(ijkbb) = 2011
1362  ENDIF
1363  ENDIF
1364  ENDIF
1365 
1366  ENDIF
1367 
1368  IF (bc_p_g(bcv) == undefined) THEN
1369  IF(dmp_log)WRITE (unit_log, 1000) 'BC_P_g', bcv
1370  call mfix_exit(mype)
1371  ELSEIF (bc_p_g(bcv)<=zero .AND. ro_g0==undefined) THEN
1372  IF(dmp_log)WRITE (unit_log, 1010) bcv, bc_p_g(bcv)
1373  call mfix_exit(mype)
1374  ENDIF
1375 
1376  ENDIF
1377 
1378  ENDIF
1379 
1380  ENDDO
1381 
1382  RETURN
1383 
1384 
1385  900 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',a,'(',i2,&
1386  ') not specified',/1x,'One of the following must be specified:',/1x,&
1387  'BC_VOLFLOW_g, BC_MASSFLOW_g or BC_VELMAG_g',/1x,70('*')/)
1388 
1389  910 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',a,'(',i2,',',i1,&
1390  ') not specified',/1x,'One of the following must be specified:',/1x,&
1391  'BC_VOLFLOW_g, BC_MASSFLOW_g or BC_VELMAG_g',/1x,70('*')/)
1392 
1393  1000 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',a,'(',i2,&
1394  ') not specified',/1x,70('*')/)
1395  1001 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/&
1396  ' Message: Illegal BC_TYPE for BC # = ',i2,/' BC_TYPE = ',a,/&
1397  ' Valid BC_TYPE are: ')
1398  1002 FORMAT(5x,a16)
1399  1003 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',a,'(',i2,&
1400  ') value is unphysical',/1x,70('*')/)
1401  1004 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',a,'(',i2,',',i2,&
1402  ') not specified',/1x,70('*')/)
1403  1005 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',a,'(',i2,',',i2,&
1404  ') value is unphysical',/1x,70('*')/)
1405  1010 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/' Message: BC_P_g( ',i2,&
1406  ') = ',g12.5,/&
1407  ' Pressure should be greater than zero for compressible flow',/1x,70(&
1408  '*')/)
1409  1050 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/' Message: BC number:',i2,&
1410  ' - ',a,' should be ',a,' zero.',/1x,70('*')/)
1411  1060 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/' Message: BC_X_g(',i2,',',i2&
1412  ,') not specified',/1x,70('*')/)
1413  1065 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/' Message: BC number:',i2,&
1414  ' - Sum of gas mass fractions is NOT equal to one',/1x,70('*')/)
1415  1100 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',a,'(',i2,',',i1,&
1416  ') not specified',/1x,70('*')/)
1417  1103 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',a,'(',i2,',',i1,&
1418  ') value is unphysical',/1x,70('*')/)
1419  1104 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',a,'(',i2,',',i2,&
1420  ',',i2,') not specified',/1x,70('*')/)
1421  1105 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',a,'(',i2,',',i2,&
1422  ',',i2,') value is unphysical',/1x,70('*')/)
1423  1110 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/' Message: BC_X_s(',i2,',',i2&
1424  ,',',i2,') not specified',/1x,70('*')/)
1425  1120 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/' Message: BC number:',i2,&
1426  ' - Sum of solids-',i1,' mass fractions is NOT equal to one',/1x,70(&
1427  '*')/)
1428  1125 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/' Message: BC number:',i2,&
1429  ' - Sum of volume fractions is NOT equal to one',/1x,70('*')/)
1430  1150 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/' Message: BC number:',i2,&
1431  ' - ',a,i1,' should be ',a,' zero.',/1x,70('*')/)
1432  1160 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/&
1433  ' Message: Boundary condition no', &
1434  i2,' is a second outflow condition.',/1x,&
1435  ' Only one outflow is allowed. Consider using P_OUTFLOW.',/1x, 70('*')/)
1436  1200 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',a,'(',i2,&
1437  ') specified',' for an undefined BC location',/1x,70('*')/)
1438  1300 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',a,'(',i2,',',i1,&
1439  ') specified',' for an undefined BC location',/1x,70('*')/)
1440  1400 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/&
1441  ' Message: No initial or boundary condition specified',/&
1442  ' I J K')
1443  1410 FORMAT(i5,3x,i5,3x,i5)
1444  1420 FORMAT(/1x,70('*')/)
1445 
1446  1500 FORMAT(/1x,70('*')//' From: CHECK_BC_FLAGS',/&
1447  ' Message: No initial or boundary condition specified',/&
1448  ' I J K')
1449 
1450 
1451  END SUBROUTINE check_bc_flags
1452 
1453 
1454 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1455 ! C
1456 ! Module name: CALL BUILD_CONE_FOR_C2C C
1457 ! Purpose: Define cone parameters for Cylider to Cylinder junction. C
1458 ! The C2C quadric ID must be between the two cylinders ID C
1459 ! (e.g., if Quadric 4 is a C2C, then Quadrics 3 and 5 C
1460 ! must be cylinders). The two cylinders must be aligned C
1461 ! in the same direction and be clipped to define the extent C
1462 ! of the conical junction. C
1463 ! This method is currentl available for internal flow only. C
1464 ! C
1465 ! C
1466 ! Author: Jeff Dietiker Date: 02-Dec-10 C
1467 ! Reviewer: Date: C
1468 ! C
1469 ! Revision Number # Date: ##-###-## C
1470 ! Author: # C
1471 ! Purpose: # C
1472 ! C
1473 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1474 !
1475  SUBROUTINE build_cone_for_c2c(Q)
1476 !
1477 !-----------------------------------------------
1478 ! M o d u l e s
1479 !-----------------------------------------------
1480  USE param
1481  USE param1
1482  USE constant
1483  USE run
1484  USE physprop
1485  USE indices
1486  USE scalars
1487  USE funits
1488  USE leqsol
1489  USE compar
1490  USE mpi_utility
1491  USE bc
1492  USE discretelement
1493 
1494  USE cutcell
1495  USE quadric
1496  USE vtk
1497  USE polygon
1498  USE dashboard
1499  USE stl
1500 
1501 
1502  IMPLICIT NONE
1503 !-----------------------------------------------
1504 ! G l o b a l P a r a m e t e r s
1505 !-----------------------------------------------
1506 !-----------------------------------------------
1507 ! L o c a l P a r a m e t e r s
1508 !-----------------------------------------------
1509 !-----------------------------------------------
1510 ! L o c a l V a r i a b l e s
1511 !-----------------------------------------------
1512  INTEGER :: Q,QM1,QP1
1513  DOUBLE PRECISION :: x1,x2,y1,y2,z1,z2,R1,R2
1514  DOUBLE PRECISION :: tan_half_angle
1515  LOGICAL :: aligned
1516 !-----------------------------------------------
1517 !
1518 
1519  qm1 = q-1
1520  qp1 = q+1
1521 
1522  IF(mype == pe_io) THEN
1523  WRITE(*,*)' INFO FOR QUADRIC', q
1524  WRITE(*,*)' Defining Cone for Cylinder to Cylinder junction'
1525  WRITE(*,*)' Between Quadrics ',qm1,' AND ', qp1
1526  ENDIF
1527 
1528 
1529  IF((trim(quadric_form(qm1))=='X_CYL_INT').AND. &
1530  (trim(quadric_form(qp1))=='X_CYL_INT')) THEN !Internal flow x-direction
1531 
1532  quadric_form(q) = 'X_CONE'
1533 
1534  aligned = (t_y(qm1)==t_y(qp1)).AND.(t_z(qm1)==t_z(qp1))
1535  IF(.NOT.aligned) THEN
1536  IF(mype == pe_io) THEN
1537  WRITE(*,*)' ERROR: CYLINDERS ',qm1, ' AND ', qp1, ' ARE NOT ALIGNED'
1538  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1539  ENDIF
1540  call mfix_exit(mype)
1541  ENDIF
1542 
1543  r1 = radius(qm1)
1544  r2 = radius(qp1)
1545  IF(r1==r2) THEN
1546  IF(mype == pe_io) THEN
1547  WRITE(*,*)' ERROR: CYLINDERS ',qm1, ' AND ', qp1, ' HAVE THE SAME RADIUS'
1548  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1549  ENDIF
1550  call mfix_exit(mype)
1551  ENDIF
1552 
1553  x1 = piece_xmax(qm1)
1554  x2 = piece_xmin(qp1)
1555  IF(x2<=x1) THEN
1556  IF(mype == pe_io) THEN
1557  WRITE(*,*)' ERROR: CYLINDERS ',qm1, ' AND ', qp1, ' ARE NOT PIECED PROPERLY'
1558  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1559  ENDIF
1560  call mfix_exit(mype)
1561  ENDIF
1562 
1563  tan_half_angle = (r2-r1)/(x2-x1)
1564 
1565  half_angle(q) = datan(tan_half_angle)/pi*180.0d0
1566  lambda_x(q) = -one
1567  lambda_y(q) = one/(tan_half_angle)**2
1568  lambda_z(q) = one/(tan_half_angle)**2
1569  dquadric(q) = zero
1570 
1571  piece_xmin(q) = x1
1572  piece_xmax(q) = x2
1573 
1574  t_x(q) = x1 - r1/tan_half_angle
1575  t_y(q) = t_y(qm1)
1576  t_z(q) = t_z(qm1)
1577 
1578  IF(mype == pe_io) THEN
1579  WRITE(*,*) ' QUADRIC:',q, ' WAS DEFINED AS ', trim(quadric_form(q))
1580  WRITE(*,*) ' WITH AN HALF-ANGLE OF ', half_angle(q), 'DEG.'
1581  ENDIF
1582 
1583 
1584  ELSEIF((trim(quadric_form(qm1))=='X_CYL_EXT').AND. &
1585  (trim(quadric_form(qp1))=='X_CYL_EXT')) THEN !External flow x-direction
1586 
1587  quadric_form(q) = 'X_CONE'
1588 
1589  aligned = (t_y(qm1)==t_y(qp1)).AND.(t_z(qm1)==t_z(qp1))
1590  IF(.NOT.aligned) THEN
1591  IF(mype == pe_io) THEN
1592  WRITE(*,*)' ERROR: CYLINDERS ',qm1, ' AND ', qp1, ' ARE NOT ALIGNED'
1593  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1594  ENDIF
1595  call mfix_exit(mype)
1596  ENDIF
1597 
1598  r1 = radius(qm1)
1599  r2 = radius(qp1)
1600  IF(r1==r2) THEN
1601  IF(mype == pe_io) THEN
1602  WRITE(*,*)' ERROR: CYLINDERS ',qm1, ' AND ', qp1, ' HAVE THE SAME RADIUS'
1603  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1604  ENDIF
1605  call mfix_exit(mype)
1606  ENDIF
1607 
1608  x1 = piece_xmax(qm1)
1609  x2 = piece_xmin(qp1)
1610  IF(x2<=x1) THEN
1611  IF(mype == pe_io) THEN
1612  WRITE(*,*)' ERROR: CYLINDERS ',qm1, ' AND ', qp1, ' ARE NOT PIECED PROPERLY'
1613  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1614  ENDIF
1615  call mfix_exit(mype)
1616  ENDIF
1617 
1618  tan_half_angle = (r2-r1)/(x2-x1)
1619 
1620  half_angle(q) = datan(tan_half_angle)/pi*180.0d0
1621  lambda_x(q) = one
1622  lambda_y(q) = -one/(tan_half_angle)**2
1623  lambda_z(q) = -one/(tan_half_angle)**2
1624  dquadric(q) = zero
1625 
1626  piece_xmin(q) = x1
1627  piece_xmax(q) = x2
1628 
1629  t_x(q) = x1 - r1/tan_half_angle
1630  t_y(q) = t_y(qm1)
1631  t_z(q) = t_z(qm1)
1632 
1633  IF(mype == pe_io) THEN
1634  WRITE(*,*) ' QUADRIC:',q, ' WAS DEFINED AS ', trim(quadric_form(q))
1635  WRITE(*,*) ' WITH AN HALF-ANGLE OF ', half_angle(q), 'DEG.'
1636  ENDIF
1637 
1638 
1639 
1640  ELSEIF((trim(quadric_form(qm1))=='Y_CYL_INT').AND. &
1641  (trim(quadric_form(qp1))=='Y_CYL_INT')) THEN !Internal flow y-direction
1642 
1643  quadric_form(q) = 'Y_CONE'
1644 
1645  aligned = (t_x(qm1)==t_x(qp1)).AND.(t_z(qm1)==t_z(qp1))
1646  IF(.NOT.aligned) THEN
1647  IF(mype == pe_io) THEN
1648  WRITE(*,*)' ERROR: CYLINDERS ',qm1, ' AND ', qp1, ' ARE NOT ALIGNED'
1649  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1650  ENDIF
1651  call mfix_exit(mype)
1652  ENDIF
1653 
1654  r1 = radius(qm1)
1655  r2 = radius(qp1)
1656  IF(r1==r2) THEN
1657  IF(mype == pe_io) THEN
1658  WRITE(*,*)' ERROR: CYLINDERS ',qm1, ' AND ', qp1, ' HAVE THE SAME RADIUS'
1659  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1660  ENDIF
1661  call mfix_exit(mype)
1662  ENDIF
1663 
1664  y1 = piece_ymax(qm1)
1665  y2 = piece_ymin(qp1)
1666  IF(y2<=y1) THEN
1667  IF(mype == pe_io) THEN
1668  WRITE(*,*)' ERROR: CYLINDERS ',qm1, ' AND ', qp1, ' ARE NOT PIECED PROPERLY'
1669  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1670  ENDIF
1671  call mfix_exit(mype)
1672  ENDIF
1673 
1674  tan_half_angle = (r2-r1)/(y2-y1)
1675 
1676  half_angle(q) = datan(tan_half_angle)/pi*180.0d0
1677  lambda_x(q) = one/(tan_half_angle)**2
1678  lambda_y(q) = -one
1679  lambda_z(q) = one/(tan_half_angle)**2
1680  dquadric(q) = zero
1681 
1682  piece_ymin(q) = y1
1683  piece_ymax(q) = y2
1684 
1685  t_x(q) = t_x(qm1)
1686  t_y(q) = y1 - r1/tan_half_angle
1687  t_z(q) = t_z(qm1)
1688 
1689  IF(mype == pe_io) THEN
1690  WRITE(*,*) ' QUADRIC:',q, ' WAS DEFINED AS ', trim(quadric_form(q))
1691  WRITE(*,*) ' WITH AN HALF-ANGLE OF ', half_angle(q), 'DEG.'
1692  ENDIF
1693 
1694  ELSEIF((trim(quadric_form(qm1))=='Y_CYL_EXT').AND. &
1695  (trim(quadric_form(qp1))=='Y_CYL_EXT')) THEN !External flow y-direction
1696 
1697  quadric_form(q) = 'Y_CONE'
1698 
1699  aligned = (t_x(qm1)==t_x(qp1)).AND.(t_z(qm1)==t_z(qp1))
1700  IF(.NOT.aligned) THEN
1701  IF(mype == pe_io) THEN
1702  WRITE(*,*)' ERROR: CYLINDERS ',qm1, ' AND ', qp1, ' ARE NOT ALIGNED'
1703  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1704  ENDIF
1705  call mfix_exit(mype)
1706  ENDIF
1707 
1708  r1 = radius(qm1)
1709  r2 = radius(qp1)
1710  IF(r1==r2) THEN
1711  IF(mype == pe_io) THEN
1712  WRITE(*,*)' ERROR: CYLINDERS ',qm1, ' AND ', qp1, ' HAVE THE SAME RADIUS'
1713  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1714  ENDIF
1715  call mfix_exit(mype)
1716  ENDIF
1717 
1718  y1 = piece_ymax(qm1)
1719  y2 = piece_ymin(qp1)
1720  IF(y2<=y1) THEN
1721  IF(mype == pe_io) THEN
1722  WRITE(*,*)' ERROR: CYLINDERS ',qm1, ' AND ', qp1, ' ARE NOT PIECED PROPERLY'
1723  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1724  ENDIF
1725  call mfix_exit(mype)
1726  ENDIF
1727 
1728  tan_half_angle = (r2-r1)/(y2-y1)
1729 
1730  half_angle(q) = datan(tan_half_angle)/pi*180.0d0
1731  lambda_x(q) = -one/(tan_half_angle)**2
1732  lambda_y(q) = one
1733  lambda_z(q) = -one/(tan_half_angle)**2
1734  dquadric(q) = zero
1735 
1736  piece_ymin(q) = y1
1737  piece_ymax(q) = y2
1738 
1739  t_x(q) = t_x(qm1)
1740  t_y(q) = y1 - r1/tan_half_angle
1741  t_z(q) = t_z(qm1)
1742 
1743  IF(mype == pe_io) THEN
1744  WRITE(*,*) ' QUADRIC:',q, ' WAS DEFINED AS ', trim(quadric_form(q))
1745  WRITE(*,*) ' WITH AN HALF-ANGLE OF ', half_angle(q), 'DEG.'
1746  ENDIF
1747 
1748 
1749 
1750  ELSEIF((trim(quadric_form(qm1))=='Z_CYL_INT').AND. &
1751  (trim(quadric_form(qp1))=='Z_CYL_INT')) THEN !Internal flow z-direction
1752 
1753  quadric_form(q) = 'Z_CONE'
1754 
1755  aligned = (t_x(qm1)==t_x(qp1)).AND.(t_y(qm1)==t_y(qp1))
1756  IF(.NOT.aligned) THEN
1757  IF(mype == pe_io) THEN
1758  WRITE(*,*)' ERROR: CYLINDERS ',qm1, ' AND ', qp1, ' ARE NOT ALIGNED'
1759  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1760  ENDIF
1761  call mfix_exit(mype)
1762  ENDIF
1763 
1764  r1 = radius(qm1)
1765  r2 = radius(qp1)
1766  IF(r1==r2) THEN
1767  IF(mype == pe_io) THEN
1768  WRITE(*,*)' ERROR: CYLINDERS ',qm1, ' AND ', qp1, ' HAVE THE SAME RADIUS'
1769  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1770  ENDIF
1771  call mfix_exit(mype)
1772  ENDIF
1773 
1774  z1 = piece_zmax(qm1)
1775  z2 = piece_zmin(qp1)
1776  IF(z2<=z1) THEN
1777  IF(mype == pe_io) THEN
1778  WRITE(*,*)' ERROR: CYLINDERS ',qm1, ' AND ', qp1, ' ARE NOT PIECED PROPERLY'
1779  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1780  ENDIF
1781  call mfix_exit(mype)
1782  ENDIF
1783 
1784  tan_half_angle = (r2-r1)/(z2-z1)
1785 
1786  half_angle(q) = datan(tan_half_angle)/pi*180.0d0
1787  lambda_x(q) = one/(tan_half_angle)**2
1788  lambda_y(q) = one/(tan_half_angle)**2
1789  lambda_z(q) = -one
1790  dquadric(q) = zero
1791 
1792  piece_zmin(q) = z1
1793  piece_zmax(q) = z2
1794 
1795  t_x(q) = t_x(qm1)
1796  t_y(q) = t_y(qm1)
1797  t_z(q) = z1 - r1/tan_half_angle
1798 
1799  IF(mype == pe_io) THEN
1800  WRITE(*,*) ' QUADRIC:',q, ' WAS DEFINED AS ', trim(quadric_form(q))
1801  WRITE(*,*) ' WITH AN HALF-ANGLE OF ', half_angle(q), 'DEG.'
1802  ENDIF
1803 
1804  ELSEIF((trim(quadric_form(qm1))=='Z_CYL_EXT').AND. &
1805  (trim(quadric_form(qp1))=='Z_CYL_EXT')) THEN !External flow z-direction
1806 
1807  quadric_form(q) = 'Z_CONE'
1808 
1809  aligned = (t_x(qm1)==t_x(qp1)).AND.(t_y(qm1)==t_y(qp1))
1810  IF(.NOT.aligned) THEN
1811  IF(mype == pe_io) THEN
1812  WRITE(*,*)' ERROR: CYLINDERS ',qm1, ' AND ', qp1, ' ARE NOT ALIGNED'
1813  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1814  ENDIF
1815  call mfix_exit(mype)
1816  ENDIF
1817 
1818  r1 = radius(qm1)
1819  r2 = radius(qp1)
1820  IF(r1==r2) THEN
1821  IF(mype == pe_io) THEN
1822  WRITE(*,*)' ERROR: CYLINDERS ',qm1, ' AND ', qp1, ' HAVE THE SAME RADIUS'
1823  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1824  ENDIF
1825  call mfix_exit(mype)
1826  ENDIF
1827 
1828  z1 = piece_zmax(qm1)
1829  z2 = piece_zmin(qp1)
1830  IF(z2<=z1) THEN
1831  IF(mype == pe_io) THEN
1832  WRITE(*,*)' ERROR: CYLINDERS ',qm1, ' AND ', qp1, ' ARE NOT PIECED PROPERLY'
1833  WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1834  ENDIF
1835  call mfix_exit(mype)
1836  ENDIF
1837 
1838  tan_half_angle = (r2-r1)/(z2-z1)
1839 
1840  half_angle(q) = datan(tan_half_angle)/pi*180.0d0
1841  lambda_x(q) = -one/(tan_half_angle)**2
1842  lambda_y(q) = -one/(tan_half_angle)**2
1843  lambda_z(q) = one
1844  dquadric(q) = zero
1845 
1846  piece_zmin(q) = z1
1847  piece_zmax(q) = z2
1848 
1849  t_x(q) = t_x(qm1)
1850  t_y(q) = t_y(qm1)
1851  t_z(q) = z1 - r1/tan_half_angle
1852 
1853  IF(mype == pe_io) THEN
1854  WRITE(*,*) ' QUADRIC:',q, ' WAS DEFINED AS ', trim(quadric_form(q))
1855  WRITE(*,*) ' WITH AN HALF-ANGLE OF ', half_angle(q), 'DEG.'
1856  ENDIF
1857 
1858 
1859  ELSE
1860  IF(mype == pe_io) THEN
1861  WRITE(*,*) ' ERROR: C2C MUST BE DEFINED BETWEEN 2 CYLINDERS'
1862  WRITE(*,*) ' QUADRIC:',qm1, ' IS ', trim(quadric_form(qm1))
1863  WRITE(*,*) ' QUADRIC:',qp1, ' IS ', trim(quadric_form(qp1))
1864  ENDIF
1865  call mfix_exit(mype)
1866 
1867  ENDIF
1868 
1869  RETURN
1870  END SUBROUTINE build_cone_for_c2c
1871 
1872 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1873 ! C
1874 ! Module name: GET_DXYZ_FROM_CONTROL_POINTS C
1875 ! Purpose: Define DX, DY, and DZ using control points C
1876 ! C
1877 ! C
1878 ! Author: Jeff Dietiker Date: 02-Dec-10 C
1879 ! Reviewer: Date: C
1880 ! C
1881 ! Revision Number # Date: ##-###-## C
1882 ! Author: # C
1883 ! Purpose: # C
1884 ! C
1885 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1886 !
1887  SUBROUTINE get_dxyz_from_control_points
1888 !
1889 !-----------------------------------------------
1890 ! M o d u l e s
1891 !-----------------------------------------------
1892  USE param
1893  USE param1
1894  USE constant
1895  USE run
1896  USE physprop
1897  USE indices
1898  USE scalars
1899  USE funits
1900  USE leqsol
1901  USE compar
1902  USE mpi_utility
1903  USE bc
1904  USE discretelement
1905 
1906  USE cutcell
1907  USE quadric
1908  USE vtk
1909  USE polygon
1910  USE dashboard
1911  USE stl
1912 
1913 
1914  IMPLICIT NONE
1915 !-----------------------------------------------
1916 ! G l o b a l P a r a m e t e r s
1917 !-----------------------------------------------
1918 !-----------------------------------------------
1919 ! L o c a l P a r a m e t e r s
1920 !-----------------------------------------------
1921 !-----------------------------------------------
1922 ! L o c a l V a r i a b l e s
1923 !-----------------------------------------------
1924 
1925  INTEGER :: NN,NX,NY,NZ
1926  INTEGER :: I,I1,I2,J,J1,J2,K,K1,K2
1927  DOUBLE PRECISION :: L,CELL_RATIO
1928 
1929  LOGICAL,DIMENSION(MAX_CP) :: INDEPENDENT_SEGMENT
1930 
1931 !-----------------------------------------------
1932 !
1933 
1934 !======================================================================
1935 ! X-DIRECTION
1936 !======================================================================
1937 
1938 ! Step 1. Input verification
1939 ! 1.1 Shift control points arrays such that the user only needs to enter
1940 ! CPX(1) and above, and CPX(0) is automatically set to zero.
1941 
1942  DO nn = max_cp,1,-1
1943  cpx(nn) = cpx(nn-1)
1944  ENDDO
1945 
1946  cpx(0) = zero
1947 
1948 ! 1.2. Last control point must match domain length.
1949 
1950  nx = 0
1951  DO nn = 1,max_cp
1952  IF(cpx(nn)>zero) nx = nx + 1
1953  ENDDO
1954 
1955  IF(nx>0) THEN
1956  IF(mype==0) WRITE(*,*)' INFO: DEFINING GRID SPACING IN X-DIRECTION... '
1957  IF(mype==0) WRITE(*,*)' INFO: NUMBER OF CONTROL POINTS IN X-DIRECTION = ',nx
1958  IF(cpx(nx)/=xlength) THEN
1959  IF(mype==0) WRITE(*,*)' ERROR: LAST CONTROL POINT MUST BE EQUAL TO XLENGTH.'
1960  IF(mype==0) WRITE(*,*)' XLENGTH = ',xlength
1961  IF(mype==0) WRITE(*,*)' LAST CONTROL POINT = ',cpx(nx)
1962  call mfix_exit(mype)
1963  ENDIF
1964  ENDIF
1965 
1966 ! 1.3. Check for acceptable values, and identify independent segments. If
1967 ! the first or last cell dimension is given, it is converted into an
1968 ! expansion ratio.
1969 
1970  independent_segment = .true.
1971 
1972  DO nn = 1,nx ! For each segment
1973 
1974  IF(cpx(nn) <= cpx(nn-1)) THEN
1975  IF(mype==0) WRITE(*,*)' ERROR: CONTROL POINTS ALONG X MUST BE SORTED IN ASCENDING ORDER.'
1976  IF(mype==0) WRITE(*,*)' CPX = ',cpx(0:nx)
1977  call mfix_exit(mype)
1978  ENDIF
1979 
1980  IF(ncx(nn) <= 1) THEN
1981  IF(mype==0) WRITE(*,*)' ERROR: NUMBER OF CELLS MUST BE LARGER THAN 1 IN X-SEGMENT :',nn
1982  IF(mype==0) WRITE(*,*)' NCX = ',ncx(nn)
1983  call mfix_exit(mype)
1984  ENDIF
1985 
1986  IF(erx(nn) <= zero) THEN
1987  IF(mype==0) WRITE(*,*)' ERROR: EXPANSION RATIO MUST BE POSITIVE IN X-SEGMENT :',nn
1988  IF(mype==0) WRITE(*,*)' ERX = ',erx(nn)
1989  call mfix_exit(mype)
1990  ENDIF
1991 
1992  ENDDO
1993 
1994  DO nn = 1,nx ! For each segment
1995 
1996  IF(first_dx(nn)/=zero.AND.last_dx(nn)/=zero) THEN
1997  IF(mype==0) WRITE(*,*)' ERROR: FIRST AND LAST DX ARE DEFINED, WHICH IS NOT ALLOWED IN X-SEGMENT :',nn
1998  IF(mype==0) WRITE(*,*)' FIRST DX = ',first_dx(nn)
1999  IF(mype==0) WRITE(*,*)' LAST DX = ',last_dx(nn)
2000  call mfix_exit(mype)
2001  ELSEIF(first_dx(nn)>zero) THEN
2002  IF(mype==0) WRITE(*,*)' INFO: FIRST DX DEFINED IN X-SEGMENT :',nn
2003  IF(mype==0) WRITE(*,*)' FIRST DX = ',first_dx(nn)
2004  l = cpx(nn) - cpx(nn-1) ! Size of the current segment
2005  IF(l<=first_dx(nn)+tol_f) THEN
2006  IF(mype==0) WRITE(*,*)' ERROR: FIRST DX IS NOT SMALLER THAN SEGMENT LENGTH IN X-SEGMENT :',nn
2007  IF(mype==0) WRITE(*,*)' FIRST DX = ',first_dx(nn)
2008  IF(mype==0) WRITE(*,*)' SEGMENT LENGTH = ',l
2009  call mfix_exit(mype)
2010  ENDIF
2011  CALL find_cell_ratio('FIRST',first_dx(nn),l,ncx(nn),cell_ratio)
2012  erx(nn) = cell_ratio**(ncx(nn)-1)
2013  IF(mype==0) WRITE(*,*)' CORRESPONDING EXPANSION RATIO = ',erx(nn)
2014  ELSEIF(last_dx(nn)>zero) THEN
2015  IF(mype==0) WRITE(*,*)' INFO: LAST DX DEFINED IN X-SEGMENT :',nn
2016  IF(mype==0) WRITE(*,*)' LAST DX = ',last_dx(nn)
2017  l = cpx(nn) - cpx(nn-1) ! Size of the current segment
2018  IF(l<=last_dx(nn)+tol_f) THEN
2019  IF(mype==0) WRITE(*,*)' ERROR: LAST DX IS NOT SMALLER THAN SEGMENT LENGTH IN X-SEGMENT :',nn
2020  IF(mype==0) WRITE(*,*)' LAST DX = ',last_dx(nn)
2021  IF(mype==0) WRITE(*,*)' SEGMENT LENGTH = ',l
2022  call mfix_exit(mype)
2023  ENDIF
2024  CALL find_cell_ratio('LAST ',last_dx(nn),l,ncx(nn),cell_ratio)
2025  erx(nn) = cell_ratio**(ncx(nn)-1)
2026  IF(mype==0) WRITE(*,*)' CORRESPONDING EXPANSION RATIO = ',erx(nn)
2027  ELSEIF(first_dx(nn)<zero) THEN
2028  IF(nn==1) THEN
2029  IF(mype==0) WRITE(*,*)' ERROR: FIRST DX CANNOT MATCH PREVIOUS DX FOR FIRST SEGMENT.'
2030  call mfix_exit(mype)
2031  ELSE
2032  IF(mype==0) WRITE(*,*)' INFO: FIRST DX WILL ATTEMPT TO MATCH PREVIOUS DX FOR SEGMENT :',nn
2033  independent_segment(nn) = .false.
2034  ENDIF
2035  ELSEIF(last_dx(nn)<zero) THEN
2036  IF(nn==nx) THEN
2037  IF(mype==0) WRITE(*,*)' ERROR: LAST DX CANNOT MATCH NEXT DX FOR LAST SEGMENT.'
2038  call mfix_exit(mype)
2039  ELSE
2040  IF(mype==0) WRITE(*,*)' INFO: LAST DX WILL ATTEMPT TO MATCH NEXT DX FOR SEGMENT :',nn
2041  independent_segment(nn) = .false.
2042  ENDIF
2043  ENDIF
2044 
2045  ENDDO
2046 
2047 ! Step 3. Computation of cell sizes.
2048 
2049 ! 3.1 First pass: Set-up all independent segments
2050 
2051 
2052  i1 = 0 ! First index of segment
2053  i2 = 0 ! Last index of segment
2054 
2055  DO nn = 1,nx ! For each segment
2056 
2057  i2 = i1 + ncx(nn) - 1
2058 
2059  IF(independent_segment(nn)) THEN
2060 
2061  l = cpx(nn) - cpx(nn-1) ! Size of the current segment
2062 
2063  IF(erx(nn)/=one) THEN
2064  cell_ratio = erx(nn)**(one/dble(ncx(nn)-1)) ! Ratio between two consecutive cells
2065  dx(i1) = l * (one - cell_ratio) / (one - cell_ratio**ncx(nn)) ! First cell size
2066 
2067  DO i = i1+1,i2 ! All other cell sizes, geometric series
2068  dx(i) = dx(i-1) * cell_ratio
2069  ENDDO
2070 
2071  ELSE
2072  dx(i1:i2) = l / ncx(nn) ! Uniform size if expansion ratio is unity.
2073  ENDIF
2074 
2075  ENDIF
2076 
2077  i1 = i2 + 1 ! Prepare First index for next segment
2078 
2079  ENDDO
2080 
2081 ! 3.2 Second pass: Set-up all dependent segments
2082 
2083 
2084  i1 = 0 ! First index of segment
2085  i2 = 0 ! Last index of segment
2086 
2087  DO nn = 1,nx ! For each segment
2088 
2089  i2 = i1 + ncx(nn) - 1
2090 
2091  IF(.NOT.independent_segment(nn)) THEN
2092 
2093  l = cpx(nn) - cpx(nn-1) ! Size of the current segment
2094 
2095  IF(first_dx(nn)<zero) THEN
2096  dx(i1) = dx(i1-1) ! First cell size
2097  CALL find_cell_ratio('FIRST',dx(i1),l,ncx(nn),cell_ratio)
2098  DO i = i1+1,i2 ! All other cell sizes, geometric series
2099  dx(i) = dx(i-1) * cell_ratio
2100  ENDDO
2101  ELSEIF(last_dx(nn)<zero) THEN
2102  dx(i2) = dx(i2+1) ! Last cell size
2103  CALL find_cell_ratio('LAST ',dx(i2),l,ncx(nn),cell_ratio)
2104  DO i = i2-1,i1,-1 ! All other cell sizes, geometric series
2105  dx(i) = dx(i+1) / cell_ratio
2106  ENDDO
2107  ENDIF
2108 
2109  ENDIF
2110 
2111  i1 = i2 + 1 ! Prepare First index for next segment
2112 
2113  ENDDO
2114 
2115 
2116 ! Step 4. Verify that the sum of cells among all segment matches the total number of cells
2117 
2118  IF(i1>0.AND.i1/=imax) THEN
2119  IF(mype==0) WRITE(*,*)' ERROR: IMAX MUST BE EQUAL TO THE SUM OF NCX.'
2120  IF(mype==0) WRITE(*,*)' IMAX = ', imax
2121  IF(mype==0) WRITE(*,*)' SUM OF NCX = ', i1
2122  call mfix_exit(mype)
2123  ENDIF
2124 
2125 
2126 !======================================================================
2127 ! Y-DIRECTION
2128 !======================================================================
2129 
2130 ! Step 1. Input verification
2131 ! 1.1 Shift control points arrays such that the user only needs to enter
2132 ! CPY(1) and above, and CPY(0) is automatically set to zero.
2133 
2134  DO nn = max_cp,1,-1
2135  cpy(nn) = cpy(nn-1)
2136  ENDDO
2137 
2138  cpy(0) = zero
2139 
2140 ! 1.2. Last control point must match domain length.
2141 
2142  ny = 0
2143  DO nn = 1,max_cp
2144  IF(cpy(nn)>zero) ny = ny + 1
2145  ENDDO
2146 
2147  IF(ny>0) THEN
2148  IF(mype==0) WRITE(*,*)' INFO: DEFINING GRID SPACING IN Y-DIRECTION... '
2149  IF(mype==0) WRITE(*,*)' INFO: NUMBER OF CONTROL POINTS IN Y-DIRECTION = ',ny
2150  IF(cpy(ny)/=ylength) THEN
2151  IF(mype==0) WRITE(*,*)' ERROR: LAST CONTROL POINT MUST BE EQUAL TO YLENGTH.'
2152  IF(mype==0) WRITE(*,*)' YLENGTH = ',ylength
2153  IF(mype==0) WRITE(*,*)' LAST CONTROL POINT = ',cpy(ny)
2154  call mfix_exit(mype)
2155  ENDIF
2156  ENDIF
2157 
2158 ! 1.3. Check for acceptable values, and identify independent segments. If
2159 ! the first or last cell dimension is given, it is converted into an
2160 ! expansion ratio.
2161 
2162  independent_segment = .true.
2163 
2164  DO nn = 1,ny ! For each segment
2165 
2166  IF(cpy(nn) <= cpy(nn-1)) THEN
2167  IF(mype==0) WRITE(*,*)' ERROR: CONTROL POINTS ALONG Y MUST BE SORTED IN ASCENDING ORDER.'
2168  IF(mype==0) WRITE(*,*)' CPY = ',cpy(0:ny)
2169  call mfix_exit(mype)
2170  ENDIF
2171 
2172  IF(ncy(nn) <= 1) THEN
2173  IF(mype==0) WRITE(*,*)' ERROR: NUMBER OF CELLS MUST BE LARGER THAN 1 IN Y-SEGMENT :',nn
2174  IF(mype==0) WRITE(*,*)' NCY = ',ncy(nn)
2175  call mfix_exit(mype)
2176  ENDIF
2177 
2178  IF(ery(nn) <= zero) THEN
2179  IF(mype==0) WRITE(*,*)' ERROR: EXPANSION RATIO MUST BE POSITIVE IN Y-SEGMENT :',nn
2180  IF(mype==0) WRITE(*,*)' ERY = ',ery(nn)
2181  call mfix_exit(mype)
2182  ENDIF
2183 
2184  ENDDO
2185 
2186  DO nn = 1,ny ! For each segment
2187 
2188  IF(first_dy(nn)/=zero.AND.last_dy(nn)/=zero) THEN
2189  IF(mype==0) WRITE(*,*)' ERROR: FIRST AND LAST DY ARE DEFINED, WHICH IS NOT ALLOWED IN Y-SEGMENT :',nn
2190  IF(mype==0) WRITE(*,*)' FIRST DY = ',first_dy(nn)
2191  IF(mype==0) WRITE(*,*)' LAST DY = ',last_dy(nn)
2192  call mfix_exit(mype)
2193  ELSEIF(first_dy(nn)>zero) THEN
2194  IF(mype==0) WRITE(*,*)' INFO: FIRST DY DEFINED IN Y-SEGMENT :',nn
2195  IF(mype==0) WRITE(*,*)' FIRST DY = ',first_dy(nn)
2196  l = cpy(nn) - cpy(nn-1) ! Size of the current segment
2197  IF(l<=first_dy(nn)+tol_f) THEN
2198  IF(mype==0) WRITE(*,*)' ERROR: FIRST DY IS NOT SMALLER THAN SEGMENT LENGTH IN Y-SEGMENT :',nn
2199  IF(mype==0) WRITE(*,*)' FIRST DY = ',first_dy(nn)
2200  IF(mype==0) WRITE(*,*)' SEGMENT LENGTH = ',l
2201  call mfix_exit(mype)
2202  ENDIF
2203  CALL find_cell_ratio('FIRST',first_dy(nn),l,ncy(nn),cell_ratio)
2204  ery(nn) = cell_ratio**(ncy(nn)-1)
2205  IF(mype==0) WRITE(*,*)' CORRESPONDING EXPANSION RATIO = ',ery(nn)
2206  ELSEIF(last_dy(nn)>zero) THEN
2207  IF(mype==0) WRITE(*,*)' INFO: LAST DY DEFINED IN Y-SEGMENT :',nn
2208  IF(mype==0) WRITE(*,*)' LAST DY = ',last_dy(nn)
2209  l = cpy(nn) - cpy(nn-1) ! Size of the current segment
2210  IF(l<=last_dy(nn)+tol_f) THEN
2211  IF(mype==0) WRITE(*,*)' ERROR: LAST DY IS NOT SMALLER THAN SEGMENT LENGTH IN Y-SEGMENT :',nn
2212  IF(mype==0) WRITE(*,*)' LAST DY = ',last_dy(nn)
2213  IF(mype==0) WRITE(*,*)' SEGMENT LENGTH = ',l
2214  call mfix_exit(mype)
2215  ENDIF
2216  CALL find_cell_ratio('LAST ',last_dy(nn),l,ncy(nn),cell_ratio)
2217  ery(nn) = cell_ratio**(ncy(nn)-1)
2218  IF(mype==0) WRITE(*,*)' CORRESPONDING EXPANSION RATIO = ',ery(nn)
2219  ELSEIF(first_dy(nn)<zero) THEN
2220  IF(nn==1) THEN
2221  IF(mype==0) WRITE(*,*)' ERROR: FIRST DY CANNOT MATCH PREVIOUS DY FOR FIRST SEGMENT.'
2222  call mfix_exit(mype)
2223  ELSE
2224  IF(mype==0) WRITE(*,*)' INFO: FIRST DY WILL ATTEMPT TO MATCH PREVIOUS DY FOR SEGMENT :',nn
2225  independent_segment(nn) = .false.
2226  ENDIF
2227  ELSEIF(last_dy(nn)<zero) THEN
2228  IF(nn==ny) THEN
2229  IF(mype==0) WRITE(*,*)' ERROR: LAST DY CANNOT MATCH NEXT DY FOR LAST SEGMENT.'
2230  call mfix_exit(mype)
2231  ELSE
2232  IF(mype==0) WRITE(*,*)' INFO: LAST DY WILL ATTEMPT TO MATCH NEXT DY FOR SEGMENT :',nn
2233  independent_segment(nn) = .false.
2234  ENDIF
2235  ENDIF
2236 
2237  ENDDO
2238 
2239 ! Step 3. Computation of cell sizes.
2240 
2241 ! 3.1 First pass: Set-up all independent segments
2242 
2243 
2244  j1 = 0 ! First index of segment
2245  j2 = 0 ! Last index of segment
2246 
2247  DO nn = 1,ny ! For each segment
2248 
2249  j2 = j1 + ncy(nn) - 1
2250 
2251  IF(independent_segment(nn)) THEN
2252 
2253  l = cpy(nn) - cpy(nn-1) ! Size of the current segment
2254 
2255  IF(ery(nn)/=one) THEN
2256  cell_ratio = ery(nn)**(one/dble(ncy(nn)-1)) ! Ratio between two consecutive cells
2257  dy(j1) = l * (one - cell_ratio) / (one - cell_ratio**ncy(nn)) ! First cell size
2258 
2259  DO j = j1+1,j2 ! All other cell sizes, geometric series
2260  dy(j) = dy(j-1) * cell_ratio
2261  ENDDO
2262 
2263  ELSE
2264  dy(j1:j2) = l / ncy(nn) ! Uniform size if expansion ratio is unity.
2265  ENDIF
2266 
2267  ENDIF
2268 
2269  j1 = j2 + 1 ! Prepare First index for next segment
2270 
2271  ENDDO
2272 
2273 ! 3.2 Second pass: Set-up all dependent segments
2274 
2275 
2276  j1 = 0 ! First index of segment
2277  j2 = 0 ! Last index of segment
2278 
2279  DO nn = 1,ny ! For each segment
2280 
2281  j2 = j1 + ncy(nn) - 1
2282 
2283  IF(.NOT.independent_segment(nn)) THEN
2284 
2285  l = cpy(nn) - cpy(nn-1) ! Size of the current segment
2286 
2287  IF(first_dy(nn)<zero) THEN
2288  dy(j1) = dy(j1-1) ! First cell size
2289  CALL find_cell_ratio('FIRST',dy(j1),l,ncy(nn),cell_ratio)
2290  DO j = j1+1,j2 ! All other cell sizes, geometric series
2291  dy(j) = dy(j-1) * cell_ratio
2292  ENDDO
2293  ELSEIF(last_dy(nn)<zero) THEN
2294  dy(j2) = dy(j2+1) ! Last cell size
2295  CALL find_cell_ratio('LAST ',dy(j2),l,ncy(nn),cell_ratio)
2296  DO j = j2-1,j1,-1 ! All other cell sizes, geometric series
2297  dy(j) = dy(j+1) / cell_ratio
2298  ENDDO
2299  ENDIF
2300 
2301  ENDIF
2302 
2303  j1 = j2 + 1 ! Prepare First index for next segment
2304 
2305  ENDDO
2306 
2307 
2308 ! Step 4. Verify that the sum of cells among all segment matches the total number of cells
2309 
2310  IF(j1>0.AND.j1/=jmax) THEN
2311  IF(mype==0) WRITE(*,*)' ERROR: JMAX MUST BE EQUAL TO THE SUM OF NCY.'
2312  IF(mype==0) WRITE(*,*)' JMAX = ', jmax
2313  IF(mype==0) WRITE(*,*)' SUM OF NCY = ', j1
2314  call mfix_exit(mype)
2315  ENDIF
2316 
2317 
2318 !======================================================================
2319 ! Z-DIRECTION
2320 !======================================================================
2321 
2322  IF(no_k) RETURN
2323 
2324 ! Step 1. Input verification
2325 ! 1.1 Shift control points arrays such that the user only needs to enter
2326 ! CPZ(1) and above, and CPZ(0) is automatically set to zero.
2327 
2328  DO nn = max_cp,1,-1
2329  cpz(nn) = cpz(nn-1)
2330  ENDDO
2331 
2332  cpz(0) = zero
2333 
2334 ! 1.2. Last control point must match domain length.
2335 
2336  nz = 0
2337  DO nn = 1,max_cp
2338  IF(cpz(nn)>zero) nz = nz + 1
2339  ENDDO
2340 
2341  IF(nz>0) THEN
2342  IF(mype==0) WRITE(*,*)' INFO: DEFINING GRID SPACING IN Z-DIRECTION... '
2343  IF(mype==0) WRITE(*,*)' INFO: NUMBER OF CONTROL POINTS IN Z-DIRECTION = ',nz
2344  IF(cpz(nz)/=zlength) THEN
2345  IF(mype==0) WRITE(*,*)' ERROR: LAST CONTROL POINT MUST BE EQUAL TO ZLENGTH.'
2346  IF(mype==0) WRITE(*,*)' ZLENGTH = ',zlength
2347  IF(mype==0) WRITE(*,*)' LAST CONTROL POINT = ',cpz(nz)
2348  call mfix_exit(mype)
2349  ENDIF
2350  ENDIF
2351 
2352 ! 1.3. Check for acceptable values, and identify independent segments. If
2353 ! the first or last cell dimension is given, it is converted into an
2354 ! expansion ratio.
2355 
2356  independent_segment = .true.
2357 
2358  DO nn = 1,nz ! For each segment
2359 
2360  IF(cpz(nn) <= cpz(nn-1)) THEN
2361  IF(mype==0) WRITE(*,*)' ERROR: CONTROL POINTS ALONG Z MUST BE SORTED IN ASCENDING ORDER.'
2362  IF(mype==0) WRITE(*,*)' CPZ = ',cpz(0:nz)
2363  call mfix_exit(mype)
2364  ENDIF
2365 
2366  IF(ncz(nn) <= 1) THEN
2367  IF(mype==0) WRITE(*,*)' ERROR: NUMBER OF CELLS MUST BE LARGER THAN 1 IN Z-SEGMENT :',nn
2368  IF(mype==0) WRITE(*,*)' NCZ = ',ncz(nn)
2369  call mfix_exit(mype)
2370  ENDIF
2371 
2372  IF(erz(nn) <= zero) THEN
2373  IF(mype==0) WRITE(*,*)' ERROR: EXPANSION RATIO MUST BE POSITIVE IN Z-SEGMENT :',nn
2374  IF(mype==0) WRITE(*,*)' ERZ = ',erz(nn)
2375  call mfix_exit(mype)
2376  ENDIF
2377 
2378  ENDDO
2379 
2380  DO nn = 1,nz ! For each segment
2381 
2382  IF(first_dz(nn)/=zero.AND.last_dz(nn)/=zero) THEN
2383  IF(mype==0) WRITE(*,*)' ERROR: FIRST AND LAST DZ ARE DEFINED, WHICH IS NOT ALLOWED IN Z-SEGMENT :',nn
2384  IF(mype==0) WRITE(*,*)' FIRST DZ = ',first_dz(nn)
2385  IF(mype==0) WRITE(*,*)' LAST DZ = ',last_dz(nn)
2386  call mfix_exit(mype)
2387  ELSEIF(first_dz(nn)>zero) THEN
2388  IF(mype==0) WRITE(*,*)' INFO: FIRST DZ DEFINED IN Z-SEGMENT :',nn
2389  IF(mype==0) WRITE(*,*)' FIRST DZ = ',first_dz(nn)
2390  l = cpz(nn) - cpz(nn-1) ! Size of the current segment
2391  IF(l<=first_dz(nn)+tol_f) THEN
2392  IF(mype==0) WRITE(*,*)' ERROR: FIRST DZ IS NOT SMALLER THAN SEGMENT LENGTH IN Z-SEGMENT :',nn
2393  IF(mype==0) WRITE(*,*)' FIRST DZ = ',first_dz(nn)
2394  IF(mype==0) WRITE(*,*)' SEGMENT LENGTH = ',l
2395  call mfix_exit(mype)
2396  ENDIF
2397  CALL find_cell_ratio('FIRST',first_dz(nn),l,ncz(nn),cell_ratio)
2398  erz(nn) = cell_ratio**(ncz(nn)-1)
2399  IF(mype==0) WRITE(*,*)' CORRESPONDING EXPANSION RATIO = ',erz(nn)
2400  ELSEIF(last_dz(nn)>zero) THEN
2401  IF(mype==0) WRITE(*,*)' INFO: LAST DZ DEFINED IN Z-SEGMENT :',nn
2402  IF(mype==0) WRITE(*,*)' LAST DZ = ',last_dz(nn)
2403  l = cpz(nn) - cpz(nn-1) ! Size of the current segment
2404  IF(l<=last_dz(nn)+tol_f) THEN
2405  IF(mype==0) WRITE(*,*)' ERROR: LAST DZ IS NOT SMALLER THAN SEGMENT LENGTH IN Z-SEGMENT :',nn
2406  IF(mype==0) WRITE(*,*)' LAST DZ = ',last_dz(nn)
2407  IF(mype==0) WRITE(*,*)' SEGMENT LENGTH = ',l
2408  call mfix_exit(mype)
2409  ENDIF
2410  CALL find_cell_ratio('LAST ',last_dz(nn),l,ncz(nn),cell_ratio)
2411  erz(nn) = cell_ratio**(ncz(nn)-1)
2412  IF(mype==0) WRITE(*,*)' CORRESPONDING EXPANSION RATIO = ',erz(nn)
2413  ELSEIF(first_dz(nn)<zero) THEN
2414  IF(nn==1) THEN
2415  IF(mype==0) WRITE(*,*)' ERROR: FIRST DZ CANNOT MATCH PREVIOUS DZ FOR FIRST SEGMENT.'
2416  call mfix_exit(mype)
2417  ELSE
2418  IF(mype==0) WRITE(*,*)' INFO: FIRST DZ WILL ATTEMPT TO MATCH PREVIOUS DZ FOR SEGMENT :',nn
2419  independent_segment(nn) = .false.
2420  ENDIF
2421  ELSEIF(last_dz(nn)<zero) THEN
2422  IF(nn==nz) THEN
2423  IF(mype==0) WRITE(*,*)' ERROR: LAST DZ CANNOT MATCH NEXT DZ FOR LAST SEGMENT.'
2424  call mfix_exit(mype)
2425  ELSE
2426  IF(mype==0) WRITE(*,*)' INFO: LAST DZ WILL ATTEMPT TO MATCH NEXT DZ FOR SEGMENT :',nn
2427  independent_segment(nn) = .false.
2428  ENDIF
2429  ENDIF
2430 
2431  ENDDO
2432 
2433 ! Step 3. Computation of cell sizes.
2434 
2435 ! 3.1 First pass: Set-up all independent segments
2436 
2437 
2438  k1 = 0 ! First index of segment
2439  k2 = 0 ! Last index of segment
2440 
2441  DO nn = 1,nz ! For each segment
2442 
2443  k2 = k1 + ncz(nn) - 1
2444 
2445  IF(independent_segment(nn)) THEN
2446 
2447  l = cpz(nn) - cpz(nn-1) ! Size of the current segment
2448 
2449  IF(erz(nn)/=one) THEN
2450  cell_ratio = erz(nn)**(one/dble(ncz(nn)-1)) ! Ratio between two consecutive cells
2451  dz(k1) = l * (one - cell_ratio) / (one - cell_ratio**ncz(nn)) ! First cell size
2452 
2453  DO k = k1+1,k2 ! All other cell sizes, geometric series
2454  dz(k) = dz(k-1) * cell_ratio
2455  ENDDO
2456 
2457  ELSE
2458  dz(k1:k2) = l / ncz(nn) ! Uniform size if expansion ratio is unity.
2459  ENDIF
2460 
2461  ENDIF
2462 
2463  k1 = k2 + 1 ! Prepare First index for next segment
2464 
2465  ENDDO
2466 
2467 ! 3.2 Second pass: Set-up all dependent segments
2468 
2469 
2470  k1 = 0 ! First index of segment
2471  k2 = 0 ! Last index of segment
2472 
2473  DO nn = 1,nz ! For each segment
2474 
2475  k2 = k1 + ncz(nn) - 1
2476 
2477  IF(.NOT.independent_segment(nn)) THEN
2478 
2479  l = cpz(nn) - cpz(nn-1) ! Size of the current segment
2480 
2481  IF(first_dz(nn)<zero) THEN
2482  dz(k1) = dz(k1-1) ! First cell size
2483  CALL find_cell_ratio('FIRST',dz(k1),l,ncz(nn),cell_ratio)
2484  DO k = k1+1,k2 ! All other cell sizes, geometric series
2485  dz(k) = dz(k-1) * cell_ratio
2486  ENDDO
2487  ELSEIF(last_dz(nn)<zero) THEN
2488  dz(k2) = dz(k2+1) ! Last cell size
2489  CALL find_cell_ratio('LAST ',dz(k2),l,ncz(nn),cell_ratio)
2490  DO k = k2-1,k1,-1 ! All other cell sizes, geometric series
2491  dz(k) = dz(k+1) / cell_ratio
2492  ENDDO
2493  ENDIF
2494 
2495  ENDIF
2496 
2497  k1 = k2 + 1 ! Prepare First index for next segment
2498 
2499  ENDDO
2500 
2501 
2502 ! Step 4. Verify that the sum of cells among all segment matches the total number of cells
2503 
2504  IF(k1>0.AND.k1/=kmax) THEN
2505  IF(mype==0) WRITE(*,*)' ERROR: KMAX MUST BE EQUAL TO THE SUM OF NCZ.'
2506  IF(mype==0) WRITE(*,*)' KMAX = ', kmax
2507  IF(mype==0) WRITE(*,*)' SUM OF NCZ = ', k1
2508  call mfix_exit(mype)
2509  ENDIF
2510 
2511 
2512 
2513  RETURN
2514 
2515  END SUBROUTINE get_dxyz_from_control_points
2516 
2517 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2518 ! C
2519 ! Module name: FIND_CELL_RATIO C
2520 ! Purpose: Given the interval length L, number of cells N, and the C
2521 ! target value of D_target, find the cell ratio alpha3 C
2522 ! such that D(POS) matches D_Target. POS can be either C
2523 ! FIRST or LAST cell in the segment. C
2524 ! C
2525 ! Author: Jeff Dietiker Date: 21-Feb-08 C
2526 ! Reviewer: Date: C
2527 ! C
2528 ! Revision Number # Date: ##-###-## C
2529 ! Author: # C
2530 ! Purpose: # C
2531 ! C
2532 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
2533  SUBROUTINE find_cell_ratio(POS,D_Target,L,NN,ALPHA3)
2534 
2535  USE param
2536  USE param1
2537  USE parallel
2538  USE constant
2539  USE run
2540  USE toleranc
2541  USE geometry
2542  USE indices
2543  USE compar
2544  USE sendrecv
2545  USE quadric
2546 
2547  IMPLICIT NONE
2548  LOGICAL :: SOLUTION_FOUND
2549 
2550  DOUBLE PRECISION :: f1,f2,f3
2551  DOUBLE PRECISION :: ALPHA1,ALPHA2,ALPHA3,D_Target,L,DU
2552  DOUBLE PRECISION, PARAMETER :: ALPHAMAX = 10.0d0 ! maximum value of cell ratio
2553  INTEGER :: NN,niter
2554  CHARACTER (LEN=5) :: POS
2555 
2556  du = l / dble(nn) ! Cell size if uniform distribution
2557 
2558  IF(du==d_target) THEN
2559  alpha3 = 1.0
2560  solution_found = .true.
2561  RETURN
2562  ELSE
2563 
2564  IF(trim(pos)=='FIRST') THEN ! Determine two initial guesses
2565  IF(d_target<du) THEN
2566  alpha1 = one
2567  alpha2 = alphamax
2568  ELSE
2569  alpha1 = one/alphamax
2570  alpha2 = one
2571  ENDIF
2572  ELSEIF(trim(pos)=='LAST') THEN
2573  IF(d_target>du) THEN
2574  alpha1 = one
2575  alpha2 = alphamax
2576  ELSE
2577  alpha1 = one/alphamax
2578  alpha2 = one
2579  ENDIF
2580  ELSE
2581  IF(mype==0) WRITE(*,*)' ERROR, IN FUNCTION F: POS MUST BE FIRST OR LAST.'
2582  call mfix_exit(mype)
2583  ENDIF
2584 
2585  ENDIF
2586 
2587  f1 = f(pos,alpha1,d_target,l,nn)
2588  f2 = f(pos,alpha2,d_target,l,nn)
2589 
2590 !======================================================================
2591 ! The cell ratio is solution of F(alpha) = zero. The root is found by
2592 ! the secant method, based on two inital guesses.
2593 !======================================================================
2594 
2595  niter = 1
2596  solution_found = .false.
2597 
2598  if(dabs(f1)<tol_f) then ! First guess is solution
2599  solution_found = .true.
2600  alpha3 = alpha1
2601  elseif(dabs(f2)<tol_f) then ! Second guess is solution
2602  solution_found = .true.
2603  alpha3 = alpha2
2604  elseif(f1*f2 < zero) then ! Solution is between two guesses
2605  niter = 0
2606  f3 = 2.0d0*tol_f
2607  do while ( (abs(f3) > tol_f) .AND. (niter<itermax_int) )
2608 
2609  alpha3 = alpha1 - f1*(alpha2-alpha1)/(f2-f1) ! secant point
2610 
2611  f3 = f(pos,alpha3,d_target,l,nn)
2612 
2613  if(f1*f3<0) then ! Reduce size of interval
2614  alpha2 = alpha3
2615  f2 = f3
2616  else
2617  alpha1 = alpha3
2618  f1 = f3
2619  endif
2620  niter = niter + 1
2621 
2622  end do
2623  if (niter < itermax_int) then
2624  solution_found = .true.
2625  else
2626  WRITE(*,*) 'Unable to find a solution'
2627  WRITE(*,1000)'between ALPHA1 = ', alpha1
2628  WRITE(*,1000)' and ALPHA2 = ', alpha2
2629  WRITE(*,1000)'Current value of ALPHA3 = ', alpha3
2630  WRITE(*,1000)'Current value of abs(f) = ', dabs(f3)
2631  WRITE(*,1000)'Tolerance = ', tol_f
2632  WRITE(*,*)'Maximum number of iterations = ', itermax_int
2633  WRITE(*,*) 'Please increase the intersection tolerance, '
2634  WRITE(*,*) 'or the maximum number of iterations, and try again.'
2635  WRITE(*,*) 'MFiX will exit now.'
2636  CALL mfix_exit(mype)
2637  solution_found = .false.
2638  endif
2639  else
2640  WRITE(*,*) 'Unable to find a solution'
2641  WRITE(*,*) 'MFiX will exit now.'
2642  CALL mfix_exit(mype)
2643  solution_found = .false.
2644  endif
2645 
2646  1000 FORMAT(a,3(2x,g12.5))
2647 
2648  RETURN
2649 
2650  contains
2651 
2652  DOUBLE PRECISION Function f(POS,ALPHAC,D_Target,L,N)
2653  use param1, only: one
2654  USE constant
2655  USE mpi_utility
2656 
2657  IMPLICIT NONE
2658  DOUBLE PRECISION:: ALPHAC,D,D_Target,DU,L
2659  INTEGER:: N
2660  CHARACTER (LEN=5) :: POS
2661 
2662  du = l / dble(nn) ! Cell size if uniform distribution
2663 
2664  IF(alphac==one) THEN
2665  d = du
2666  ELSE
2667  IF(trim(pos)=='FIRST') THEN
2668  d = l * (one - alphac) / (one -alphac**n)
2669  ELSEIF(trim(pos)=='LAST') THEN
2670  d = l * (one - alphac) / (one -alphac**n) * alphac**(n-1)
2671  ELSE
2672  IF(mype==0) WRITE(*,*)' ERROR, IN FUNCTION F: POS MUST BE FIRST OR LAST.'
2673  call mfix_exit(mype)
2674  ENDIF
2675  ENDIF
2676 
2677  f = d - d_target
2678 
2679  RETURN
2680 
2681  END FUNCTION f
2682 
2683  END SUBROUTINE find_cell_ratio
2684 
2685 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2686 ! C
2687 ! Module name: ADJUST_IJK_SIZE C
2688 ! Purpose: Adjust domain size of each processor, based on an C
2689 ! estimate on the total number of useful cells. C
2690 ! The objective is to reduce load imbalance. C
2691 ! This is the parallel version C
2692 ! C
2693 ! C
2694 ! Author: Jeff Dietiker Date: 02-Dec-10 C
2695 ! Reviewer: Date: C
2696 ! C
2697 ! Revision Number # Date: ##-###-## C
2698 ! Author: # C
2699 ! Purpose: # C
2700 ! C
2701 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
2702 !
2703  SUBROUTINE adjust_ijk_size
2704 !
2705 !-----------------------------------------------
2706 ! M o d u l e s
2707 !-----------------------------------------------
2708  USE bc
2709  USE compar
2710  USE constant
2711  USE cut_cell_preproc, only: eval_f, print_cg_header
2712  USE cutcell
2713  USE dashboard
2714  USE discretelement
2715  USE funits
2716  USE gridmap
2717  USE indices
2718  USE leqsol
2719  USE mpi_utility
2720  USE param
2721  USE param1
2722  USE physprop
2723  USE polygon
2724  USE quadric
2725  USE run
2726  USE scalars
2727  USE stl
2728  USE vtk
2729 
2730  IMPLICIT NONE
2731 !-----------------------------------------------
2732 ! G l o b a l P a r a m e t e r s
2733 !-----------------------------------------------
2734 !-----------------------------------------------
2735 ! L o c a l P a r a m e t e r s
2736 !-----------------------------------------------
2737 !-----------------------------------------------
2738 ! L o c a l V a r i a b l e s
2739 !-----------------------------------------------
2740  INTEGER :: LC,I,J,K,Q_ID,TOTAL_NUC,IDEAL_NCPP
2741  DOUBLE PRECISION :: X_COPY,Y_COPY,Z_COPY,F_COPY
2742  LOGICAL :: SHIFT,CLIP_FLAG,PRESENT
2743  DOUBLE PRECISION, DIMENSION(0:DIM_I) :: DXT ,XCC
2744  DOUBLE PRECISION, DIMENSION(0:DIM_J) :: DYT ,YCC
2745  DOUBLE PRECISION, DIMENSION(0:DIM_K) :: DZT, ZCC
2746 
2747  INTEGER, DIMENSION(0:DIM_I) :: NUC_I,GLOBAL_NUC_I
2748  INTEGER, DIMENSION(0:DIM_J) :: NUC_J,GLOBAL_NUC_J
2749  INTEGER, DIMENSION(0:DIM_K) :: NUC_K,GLOBAL_NUC_K
2750 
2751  INTEGER :: IPROC,PSUM
2752 
2753  INTEGER, DIMENSION(0:numPEs-1) :: NCPP_OLD,NCPP,NCPP_WITH_GHOST
2754 
2755  INTEGER, DIMENSION(0:NODESJ-1) :: JSIZE_OLD
2756 
2757  INTEGER :: JSIZE, JREMAIN
2758  INTEGER :: MAXVAL_NCPP_OLD,MINVAL_NCPP_OLD,MAXVAL_NCPP,MINVAL_NCPP
2759  INTEGER :: AVG_NCPP_OLD,AVG_NCPP
2760  DOUBLE PRECISION :: LIP_OLD,LIP,MAXSPEEDUP_OLD,P
2761 
2762  INTEGER :: I_OFFSET,J_OFFSET,K_OFFSET,IERR
2763 
2764  INTEGER, DIMENSION(0:numPEs-1) :: disp,rcount
2765 
2766  LOGICAL :: PRINT_STATISTICS
2767 
2768 
2769 !-----------------------------------------------
2770 !
2771 
2772  IF(.NOT.cartesian_grid) RETURN ! Perform adjustement only when both CG
2773  IF(.NOT.adjust_proc_domain_size) RETURN ! and domain adjustment
2774  IF(nodesi*nodesj*nodesk==1) RETURN ! and parallel run are active
2775 
2776 
2777  INQUIRE(file='gridmap.dat',exist=present)
2778  IF(present) THEN
2779  IF (mype == pe_io) WRITE(*,*)'gridmap was assigned from grimap.dat. Skipping the adjustment.'
2780  RETURN
2781 
2782  ENDIF
2783 
2784  IF (mype == pe_io) CALL print_cg_header
2785 
2786  allocate( ncpp_uniform(0:numpes-1))
2787 
2788  short_gridmap_init = .true.
2789  CALL gridmap_init
2790  short_gridmap_init = .false.
2791 
2792  allocate( isize_all(0:nodesi-1))
2793  allocate( jsize_all(0:nodesj-1))
2794  allocate( ksize_all(0:nodesk-1))
2795 
2796 
2797 
2798  ncpp_uniform = ijksize3_all
2799 
2800 ! print*,'------> MyPE,NCPP_UNIFORM=',MyPE,NCPP_UNIFORM
2801 
2802  dimension_i = imax3
2803  dimension_j = jmax3
2804  dimension_k = kmax3
2805 
2806  partial_check_03 = .true.
2807 ! CALL CHECK_DATA_03(SHIFT)
2808  partial_check_03 = .false.
2809  CALL check_data_cartesian ! Make sure CG data is valid
2810 
2811  CALL define_quadrics
2812 
2813  shift = .true.
2814 
2815  IF(shift) THEN ! Shift DX,DY,DZ and store it into temporary DXT,DYT,DZT
2816 
2817  IF (do_i) THEN
2818  dxt(imax3) = dx(imax-1)
2819  dxt(imax2) = dx(imax-1)
2820  DO lc = imax1, imin1, -1
2821  dxt(lc) = dx(lc-2)
2822  ENDDO
2823  dxt(imin2) = dx(imin1)
2824  dxt(imin3) =dx(imin2)
2825 
2826  xcc(imin1) = half*dxt(imin1)
2827  DO i=imin1+1,imax1
2828  xcc(i) = xcc(i-1) + half*(dxt(i-1) + dxt(i))
2829  ENDDO
2830 
2831  ENDIF
2832  !
2833  IF (do_j) THEN
2834  dyt(jmax3) = dy(jmax-1)
2835  dyt(jmax2) = dy(jmax-1)
2836  DO lc = jmax1, jmin1, -1
2837  dyt(lc) = dy(lc-2)
2838  ENDDO
2839  dyt(jmin2) = dy(jmin1)
2840  dyt(jmin3) =dy(jmin2)
2841 
2842  ycc(jmin1) = half*dyt(jmin1)
2843  DO j=jmin1+1,jmax1
2844  ycc(j) = ycc(j-1) + half*(dyt(j-1) + dyt(j))
2845  ENDDO
2846 
2847 
2848  ENDIF
2849  !
2850  IF (do_k) THEN
2851  dzt(kmax3) = dz(kmax-1)
2852  dzt(kmax2) = dz(kmax-1)
2853  DO lc = kmax1, kmin1, -1
2854  dzt(lc) = dz(lc-2)
2855  ENDDO
2856  dzt(kmin2) = dz(kmin1)
2857  dzt(kmin3) =dz(kmin2)
2858 
2859  zcc(kmin1) = half*dzt(kmin1)
2860  DO k=kmin1+1,kmax1
2861  zcc(k) = zcc(k-1) + half*(dzt(k-1) + dzt(k))
2862  ENDDO
2863 
2864  ENDIF
2865 
2866  ENDIF ! SHIFT
2867 
2868 
2869 
2870  IF(nodesi>1) THEN ! DOMAIN DECOMPOSITION IN I-DIRECTION
2871 
2872  IF(mype == 0.AND.nodesj*nodesk/=1) THEN
2873  WRITE(*,*)'ERROR IN SUBROUTINE ADJUST_IJK_SIZE.'
2874  WRITE(*,*)'ADJUSTMENT POSSIBLE ONLY FOR DOMAIN DECOMPOSITION In ONE DIRECTION.'
2875  WRITE(*,*)'NODESI,NODESJ,NODESK =',nodesi,nodesj,nodesk
2876  WRITE(*,*)'MFIX WILL EXIT NOW.'
2877  CALL mfix_exit(mype)
2878  ENDIF
2879 
2880 
2881  jsize_all(0:nodesj-1) = jmax1-jmin1+1 ! Assign size in J and K-direction
2882  ksize_all(0:nodesk-1) = kmax1-kmin1+1
2883 
2884 ! Assign size in I-direction
2885  IF (mype == 0) THEN
2886  WRITE(*,1000)'INFO FROM ADJUST_IJK_SIZE:'
2887  WRITE(*,1000)'ATTEMPTING TO ADJUST DOMAIN SIZE IN I-DIRECTION ...'
2888  WRITE(*,1000)'THIS IS BASED ON AN ESTIMATED (NOT EXACT) NUMBER OF USEFUL CELLS,'
2889  WRITE(*,1000)'AND IT INCLUDES GHOST LAYERS.'
2890 
2891  ENDIF
2892 
2893  DO i = istart1,iend1
2894 
2895  nuc_i(i) = 0 ! NUC : Number of Useful Cells
2896 
2897  DO j = jstart1,jend1
2898  DO k = kstart1,kend1
2899  q_id = 1
2900  clip_flag = .false.
2901 
2902  x_copy = xcc(i)
2903  y_copy = ycc(j)
2904  z_copy = zcc(k)
2905 
2906 
2907  CALL eval_f('QUADRIC',x_copy,y_copy,z_copy,q_id,f_copy,clip_flag)
2908 
2909  ! CALL EVAL_F('POLYGON',X_COPY,Y_COPY,Z_COPY,N_POLYGON,F_COPY,CLIP_FLAG)
2910 
2911  ! CALL EVAL_F('USR_DEF',X_COPY,Y_COPY,Z_COPY,N_USR_DEF,F_COPY,CLIP_FLAG)
2912 
2913  ! CALL EVAL_STL_FCT(x1,x2,x3,Q,f,CLIP_FLAG,BCID)
2914 
2915 
2916  IF (f_copy < tol_f ) THEN ! Interior point, counted as useful
2917  nuc_i(i) = nuc_i(i) + 1
2918  ENDIF
2919 
2920  ENDDO
2921  ENDDO
2922 
2923  ENDDO
2924 
2925 ! Gather NUC onto the head node
2926 
2927  CALL allgather_1i (iend1-istart1+1,rcount,ierr)
2928 
2929  IF (mype == 0) THEN
2930  i_offset = 0
2931  ELSE
2932  i_offset = 0
2933  DO iproc=0,mype-1
2934  i_offset = i_offset + rcount(iproc)
2935  ENDDO
2936  ENDIF
2937 
2938  CALL allgather_1i (i_offset,disp,ierr)
2939 
2940  call gatherv_1i( nuc_i(istart1:iend1), iend1-istart1+1, global_nuc_i(imin1:imax1), rcount, disp, pe_io, ierr )
2941 
2942 
2943 
2944  ELSEIF(nodesj>1) THEN ! DOMAIN DECOMPOSITION IN J-DIRECTION
2945 
2946 
2947  IF(mype == 0.AND.nodesi*nodesk/=1) THEN
2948  WRITE(*,*)'ERROR IN SUBROUTINE ADJUST_IJK_SIZE.'
2949  WRITE(*,*)'ADJUSTMENT POSSIBLE ONLY FOR DOMAIN DECOMPOSITION In ONE DIRECTION.'
2950  WRITE(*,*)'NODESI,NODESJ,NODESK =',nodesi,nodesj,nodesk
2951  WRITE(*,*)'MFIX WILL EXIT NOW.'
2952  CALL mfix_exit(mype)
2953  ENDIF
2954 
2955 
2956  isize_all(0:nodesi-1) = imax1-imin1+1 ! Assign size in I and K-direction
2957  ksize_all(0:nodesk-1) = kmax1-kmin1+1
2958 
2959 
2960 
2961 ! WRITE(*,*)'Attempting to read gridmap from grimap.dat...',MyPE
2962 ! INQUIRE(FILE='gridmap.dat',EXIST=PRESENT)
2963 ! IF(PRESENT) THEN
2964 ! WRITE(*,*)'Reading gridmap from grimap.dat...'
2965 ! OPEN(CONVERT='BIG_ENDIAN',UNIT=777, FILE='gridmap.dat', STATUS='OLD')
2966 ! DO IPROC = 0,NumPEs-1
2967 ! READ(777,*) jsize_all(IPROC)
2968 ! ENDDO
2969 ! CLOSE(777)
2970 ! GOTO 9999
2971 ! ENDIF
2972 
2973 
2974 ! Assign size in J-direction
2975  IF (mype == 0) THEN
2976  WRITE(*,1000)'INFO FROM ADJUST_IJK_SIZE:'
2977  WRITE(*,1000)'ATTEMPTING TO ADJUST DOMAIN SIZE IN J-DIRECTION ...'
2978  WRITE(*,1000)'THIS IS BASED ON AN ESTIMATED (NOT EXACT) NUMBER OF USEFUL CELLS,'
2979  WRITE(*,1000)'AND IT INCLUDES GHOST LAYERS.'
2980 
2981  ENDIF
2982 
2983  DO j = jstart1,jend1
2984 
2985  nuc_j(j) = 0 ! NUC : Number of Useful Cells
2986 
2987  DO i = istart1,iend1
2988  DO k = kstart1,kend1
2989  q_id = 1
2990  clip_flag = .false.
2991 
2992  x_copy = xcc(i)
2993  y_copy = ycc(j)
2994  z_copy = zcc(k)
2995 
2996 
2997  CALL eval_f('QUADRIC',x_copy,y_copy,z_copy,q_id,f_copy,clip_flag)
2998 
2999  ! CALL EVAL_F('POLYGON',X_COPY,Y_COPY,Z_COPY,N_POLYGON,F_COPY,CLIP_FLAG)
3000 
3001  ! CALL EVAL_F('USR_DEF',X_COPY,Y_COPY,Z_COPY,N_USR_DEF,F_COPY,CLIP_FLAG)
3002 
3003  ! CALL EVAL_STL_FCT(x1,x2,x3,Q,f,CLIP_FLAG,BCID)
3004 
3005 
3006  IF (f_copy < tol_f ) THEN ! Interior point, counted as useful
3007  nuc_j(j) = nuc_j(j) + 1
3008  ENDIF
3009 
3010  ENDDO
3011  ENDDO
3012 
3013  ENDDO
3014 
3015 ! Gather NUC onto the head node
3016 
3017  CALL allgather_1i (jend1-jstart1+1,rcount,ierr)
3018 
3019  IF (mype == 0) THEN
3020  j_offset = 0
3021  ELSE
3022  j_offset = 0
3023  DO iproc=0,mype-1
3024  j_offset = j_offset + rcount(iproc)
3025  ENDDO
3026  ENDIF
3027 
3028  CALL allgather_1i (j_offset,disp,ierr)
3029 
3030  call gatherv_1i( nuc_j(jstart1:jend1), jend1-jstart1+1, global_nuc_j(jmin1:jmax1), rcount, disp, pe_io, ierr )
3031 
3032 
3033  ELSEIF(nodesk>1) THEN ! DOMAIN DECOMPOSITION IN K-DIRECTION
3034 
3035  IF(mype == 0.AND.nodesi*nodesj/=1) THEN
3036  WRITE(*,*)'ERROR IN SUBROUTINE ADJUST_IJK_SIZE.'
3037  WRITE(*,*)'ADJUSTMENT POSSIBLE ONLY FOR DOMAIN DECOMPOSITION In ONE DIRECTION.'
3038  WRITE(*,*)'NODESI,NODESJ,NODESK =',nodesi,nodesj,nodesk
3039  WRITE(*,*)'MFIX WILL EXIT NOW.'
3040  CALL mfix_exit(mype)
3041  ENDIF
3042 
3043 
3044  isize_all(0:nodesi-1) = imax1-imin1+1 ! Assign size in I and J-direction
3045  jsize_all(0:nodesj-1) = jmax1-jmin1+1
3046 
3047 ! Assign size in K-direction
3048  IF (mype == 0) THEN
3049  WRITE(*,1000)'INFO FROM ADJUST_IJK_SIZE:'
3050  WRITE(*,1000)'ATTEMPTING TO ADJUST DOMAIN SIZE IN K-DIRECTION ...'
3051  WRITE(*,1000)'THIS IS BASED ON AN ESTIMATED (NOT EXACT) NUMBER OF USEFUL CELLS,'
3052  WRITE(*,1000)'AND IT INCLUDES GHOST LAYERS.'
3053  ENDIF
3054 
3055  DO k = kstart1,kend1
3056 
3057  nuc_k(k) = 0 ! NUC : Number of Useful Cells
3058 
3059  DO i = istart1,iend1
3060  DO j = jstart1,jend1
3061  q_id = 1
3062  clip_flag = .false.
3063 
3064  x_copy = xcc(i)
3065  y_copy = ycc(j)
3066  z_copy = zcc(k)
3067 
3068 
3069  CALL eval_f('QUADRIC',x_copy,y_copy,z_copy,q_id,f_copy,clip_flag)
3070 
3071  ! CALL EVAL_F('POLYGON',X_COPY,Y_COPY,Z_COPY,N_POLYGON,F_COPY,CLIP_FLAG)
3072 
3073  ! CALL EVAL_F('USR_DEF',X_COPY,Y_COPY,Z_COPY,N_USR_DEF,F_COPY,CLIP_FLAG)
3074 
3075  ! CALL EVAL_STL_FCT(x1,x2,x3,Q,f,CLIP_FLAG,BCID)
3076 
3077 
3078  IF (f_copy < tol_f ) THEN ! Interior point, counted as useful
3079  nuc_k(k) = nuc_k(k) + 1
3080  ENDIF
3081 
3082  ENDDO
3083  ENDDO
3084 
3085  ENDDO
3086 
3087 ! Gather NUC onto the head node
3088 
3089  CALL allgather_1i (kend1-kstart1+1,rcount,ierr)
3090 
3091  IF (mype == 0) THEN
3092  k_offset = 0
3093  ELSE
3094  k_offset = 0
3095  DO iproc=0,mype-1
3096  k_offset = k_offset + rcount(iproc)
3097  ENDDO
3098  ENDIF
3099 
3100  CALL allgather_1i (k_offset,disp,ierr)
3101 
3102  call gatherv_1i( nuc_k(kstart1:kend1), kend1-kstart1+1, global_nuc_k(kmin1:kmax1), rcount, disp, pe_io, ierr )
3103 
3104  ENDIF !(NODESI,NODESJ, OR NODESK)
3105 
3106  IF (mype == pe_io) THEN ! DETERMINE BEST DOMAIN DECOMPOSITION FROM HEAD NODE
3107 
3108  IF(nodesi>1) THEN ! DOMAIN DECOMPOSITION IN I-DIRECTION
3109 
3110  ELSEIF(nodesj>1) THEN ! DOMAIN DECOMPOSITION IN J-DIRECTION
3111 
3112  ! For comparison with old grid size, determine the size in j direction (without adjustment) and add the remainder sequentially
3113 
3114  jsize = (jmax1-jmin1+1)/nodesj
3115  jsize_old(0:nodesj-1) = jsize
3116 
3117  jremain = (jmax1-jmin1+1) - nodesj*jsize
3118  IF (jremain.GE.1) THEN
3119  jsize_old( 0:(jremain-1) ) = jsize + 1
3120  ENDIF
3121 
3122 
3123 ! DO IPROC = 0 ,numPEs-1
3124 ! NCPP_UNIFORM(IPROC) = (imax3-imin3+1)*(JSIZE_OLD(IPROC)+4)
3125 ! IF(DO_K) NCPP_UNIFORM(IPROC) = NCPP_UNIFORM(IPROC)*(kmax3-kmin3+1)
3126 ! ENDDO
3127 
3128 
3129 
3130  jsize_all = jsize_old
3131 
3132  CALL minimize_load_imbalance(nodesj,global_nuc_j(jmin1:jmax1),jmin1,jmax1,jsize_all,ncpp,ncpp_with_ghost)
3133 
3134 
3135  total_nuc = sum(ncpp_with_ghost(0:numpes-1))
3136  ideal_ncpp = total_nuc / numpes
3137 
3138  WRITE (*, 1000) 'AFTER OPTIMIZATION:'
3139  WRITE (*, 1010) 'TOTAL NUMBER OF USEFUL CELLS = ',total_nuc
3140  WRITE (*, 1010) 'IDEAL NUMBER OF CELLS/PROC. = ',ideal_ncpp
3141  WRITE (*, 1000) 'ACTUALL CELL COUNT:'
3142  WRITE (*, 1000) '================================================='
3143  WRITE (*, 1000) ' PROCESSOR J-SIZE CELLS/PROC. DIFF. (%)'
3144  WRITE (*, 1000) '================================================='
3145  DO iproc = 0,numpes-1
3146  WRITE (*, 1020) iproc,jsize_all(iproc),ncpp_with_ghost(iproc), &
3147  dble(ncpp_with_ghost(iproc)-ideal_ncpp)/dble(ideal_ncpp)*100.0d0
3148  ENDDO
3149 
3150  DO iproc = 0 ,numpes-1
3151  ncpp_old(iproc) = (imax3-imin3+1)*(jsize_all(iproc)+4)
3152  IF(do_k) ncpp_old(iproc) = ncpp_old(iproc)*(kmax3-kmin3+1)
3153  ENDDO
3154 
3155 
3156  ! Verify that the sum of all JSIZE_ALL matches JMAX
3157 
3158  psum = 0
3159  DO iproc = 0,numpes-1
3160  psum = psum + jsize_all(iproc)
3161  IF(jsize_all(iproc)<5) THEN
3162  WRITE (*, 1010) 'ERROR: J-SIZE TOO SMALL FOR PROCESSOR:',iproc
3163  WRITE (*, 1010) 'J-SIZE = ',jsize_all(iproc)
3164  CALL mfix_exit(mype)
3165  ENDIF
3166 
3167  ENDDO
3168 
3169  IF(psum/=jmax) THEN
3170  WRITE (*, 1000) 'ERROR IN ADJUST_IJK_SIZE: UNABLE TO ASSIGN JSIZE TO PROCESSORS.'
3171  WRITE (*, 1000) 'SUM OF JSIZE_ALL DOES NOT MATCH JMAX:'
3172  WRITE (*, 1010) 'SUM OF JSIZE_ALL = ',psum
3173  WRITE (*, 1010) 'JMAX1 = ',jmax
3174  CALL mfix_exit(mype)
3175  ENDIF
3176 
3177  ELSEIF(nodesk>1) THEN ! DOMAIN DECOMPOSITION IN K-DIRECTION
3178 
3179  ENDIF !(NODESI,NODESJ, OR NODESK)
3180 
3181  print_statistics = .true.
3182 
3183  IF(print_statistics) THEN
3184 
3185 ! MAXVAL_NCPP_OLD = MAXVAL(NCPP_OLD)
3186 ! MINVAL_NCPP_OLD = MINVAL(NCPP_OLD)
3187 ! AVG_NCPP_OLD = SUM(NCPP_OLD)/NUMPES
3188 
3189 ! LIP_OLD = DBLE(MAXVAL_NCPP_OLD-AVG_NCPP_OLD)/DBLE(AVG_NCPP_OLD)*100.0D0
3190 
3191 ! P = DBLE(MAXVAL_NCPP_OLD)/DBLE(AVG_NCPP_OLD)
3192 
3193 ! MAXSPEEDUP_OLD = DBLE(NumPes)*(ONE-LIP_OLD/100.0D0)
3194 
3195 
3196 
3197  maxval_ncpp_old = maxval(ncpp_uniform)
3198  minval_ncpp_old = minval(ncpp_uniform)
3199  avg_ncpp_old = sum(ncpp_uniform)/numpes
3200 
3201  lip_old = dble(maxval_ncpp_old-avg_ncpp_old)/dble(avg_ncpp_old)*100.0d0
3202 
3203  p = dble(maxval_ncpp_old)/dble(avg_ncpp_old)
3204 
3205  maxspeedup_old = dble(numpes)*(one-lip_old/100.0d0)
3206 
3207 
3208 
3209 
3210  maxval_ncpp = maxval(ncpp_with_ghost)
3211  minval_ncpp = minval(ncpp_with_ghost)
3212  avg_ncpp = sum(ncpp_with_ghost(0:numpes-1))/numpes
3213 
3214  lip = dble(maxval_ncpp-avg_ncpp)/dble(avg_ncpp)*100.0d0
3215 
3216  p = dble(maxval_ncpp)/dble(avg_ncpp)
3217 
3218 ! MAXSPEEDUP = DBLE(NumPes)*(ONE-LIP/100.0D0)
3219 
3220 
3221  WRITE (*, 1000) '================================================='
3222  WRITE (*, 1000) 'ESTIMATED PARALLEL LOAD BALANCING STATISTICS'
3223  WRITE (*, 1000) 'COMPARISION BETWEEN UNIFORM SIZE (OLD)'
3224  WRITE (*, 1000) 'AND ADJUSTED SIZE (NEW)'
3225  WRITE (*, 1000) '================================================='
3226  WRITE (*, 1000) ' OLD NEW'
3227  WRITE (*, 1010) 'MAX CELL COUNT : ',maxval_ncpp_old,maxval_ncpp
3228  WRITE (*, 1010) 'AT PROCESSOR : ',maxloc(ncpp_old)-1,maxloc(ncpp_with_ghost)-1
3229  WRITE (*, 1010) 'MIN CELL COUNT : ',minval_ncpp_old,minval_ncpp
3230  WRITE (*, 1010) 'AT PROCESSOR : ',minloc(ncpp_old)-1,minloc(ncpp_with_ghost)-1
3231  WRITE (*, 1010) 'AVG CELL COUNT : ',avg_ncpp_old,avg_ncpp
3232  WRITE (*, 1000) ''
3233  WRITE (*, 1030) 'LOAD IMBALANCE (%) : ',lip_old,lip
3234  WRITE (*, 1000) ''
3235 ! WRITE (*, 1030) 'IDEAL SPEEDUP : ',DBLE(NumPEs),DBLE(NumPEs)
3236 ! WRITE (*, 1030) 'MAX SPEEDUP : ',MAXSPEEDUP_OLD,MAXSPEEDUP
3237 ! WRITE (*, 1030) 'MAX EFFICIENCY (%) : ',100.0D0 - LIP_OLD,100.0D0 - LIP
3238 
3239  WRITE (*, 1000) '================================================='
3240 
3241  WRITE (*, 1000) 'NOTE: ACTUAL LOAD BALANCING WILL BE COMPUTED AFTER PRE_PROCESSING.'
3242 
3243 
3244  ENDIF !(PRINT_STATISTICS)
3245 
3246  ENDIF ! (MyPE==PE_IO)
3247 
3248 9999 CONTINUE
3249 
3250 ! Broadcast Domain sizes to all processors
3251 
3252  CALL bcast(isize_all)
3253  CALL bcast(jsize_all)
3254  CALL bcast(ksize_all)
3255 
3256  domain_size_adjusted = .true.
3257 
3258 
3259 1000 FORMAT(1x,a)
3260 1010 FORMAT(1x,a,i10,i10)
3261 1020 FORMAT(1x,i8,2(i12),f12.2)
3262 1030 FORMAT(1x,a,2(f10.1))
3263 1040 FORMAT(f10.1)
3264 1050 FORMAT(1x,3(a))
3265 
3266 
3267  RETURN
3268 
3269  END SUBROUTINE adjust_ijk_size
3270 
3271 
3272 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
3273 ! C
3274 ! Module name: REPORT_BEST_PROCESSOR_SIZE C
3275 ! Purpose: Adjust domain size of each processor, based on an C
3276 ! estimate on the total number of useful cells. C
3277 ! The objective is to reduce load imbalance. C
3278 ! This subroutine can be called in serial and parallel mode C
3279 ! C
3280 ! C
3281 ! Author: Jeff Dietiker Date: 02-Dec-10 C
3282 ! Reviewer: Date: C
3283 ! C
3284 ! Revision Number # Date: ##-###-## C
3285 ! Author: # C
3286 ! Purpose: # C
3287 ! C
3288 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
3289 !
3290  SUBROUTINE report_best_processor_size
3291 !
3292 !-----------------------------------------------
3293 ! M o d u l e s
3294 !-----------------------------------------------
3295  USE gridmap
3296  USE param
3297  USE param1
3298  USE constant
3299  USE run
3300  USE physprop
3301  USE indices
3302  USE scalars
3303  USE funits
3304  USE leqsol
3305  USE compar
3306  USE mpi_utility
3307  USE bc
3308  USE discretelement
3309 
3310  USE parallel
3311 
3312  USE cutcell
3313  USE quadric
3314  USE vtk
3315  USE polygon
3316  USE dashboard
3317  USE stl
3318 
3319 
3320  IMPLICIT NONE
3321 !-----------------------------------------------
3322 ! G l o b a l P a r a m e t e r s
3323 !-----------------------------------------------
3324 !-----------------------------------------------
3325 ! L o c a l P a r a m e t e r s
3326 !-----------------------------------------------
3327 !-----------------------------------------------
3328 !
3329 
3330  is_serial=(numpes==1)
3331 
3332  IF(is_serial) THEN ! Temporarily mimick a parallel run
3333 
3334  nodesi = nodesi_report
3335  nodesj = nodesj_report
3336  nodesk = nodesk_report
3337  numpes = nodesi*nodesj*nodesk
3338 
3339  IF(numpes>1.AND.mype==0) THEN
3340  WRITE(*,1000)'TEMPORARILY SETTING:'
3341  WRITE(*,1010)'NODESI = ',nodesi
3342  WRITE(*,1010)'NODESJ = ',nodesj
3343  WRITE(*,1010)'NODESK = ',nodesk
3344  WRITE(*,1000)'TO REPORT BEST DOMAIN SIZE FOR PARALLEL RUN'
3345  ENDIF
3346 
3347 
3348  ENDIF
3349 
3350 
3351 
3352  IF(numpes>1) CALL report_best_ijk_size
3353 
3354 
3355  IF(is_serial) THEN ! Revert to serial values
3356  ! These values were changed to allow reporting
3357  ! optimized sizes even for a serial run
3358 
3359  nodesi = 1
3360  nodesj = 1
3361  nodesk = 1
3362  numpes = 1
3363 
3364  ENDIF
3365 
3366 1000 FORMAT(1x,a)
3367 1010 FORMAT(1x,a,i10)
3368 
3369  RETURN
3370 
3371  END SUBROUTINE report_best_processor_size
3372 
3373 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
3374 ! C
3375 ! Module name: REPORT_BEST_IJK_SIZE C
3376 ! Purpose: Adjust domain size of each processor, based on an C
3377 ! estimate on the total number of useful cells. C
3378 ! The objective is to reduce load imbalance. C
3379 ! This is the parallel version C
3380 ! C
3381 ! C
3382 ! Author: Jeff Dietiker Date: 02-Dec-10 C
3383 ! Reviewer: Date: C
3384 ! C
3385 ! Revision Number # Date: ##-###-## C
3386 ! Author: # C
3387 ! Purpose: # C
3388 ! C
3389 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
3390 !
3391  SUBROUTINE report_best_ijk_size
3392 !
3393 !-----------------------------------------------
3394 ! M o d u l e s
3395 !-----------------------------------------------
3396  USE gridmap
3397  USE param
3398  USE param1
3399  USE constant
3400  USE run
3401  USE physprop
3402  USE indices
3403  USE scalars
3404  USE funits
3405  USE leqsol
3406  USE compar
3407  USE mpi_utility
3408  USE bc
3409  USE discretelement
3410 
3411  USE parallel
3412 
3413  USE cutcell
3414  USE quadric
3415  USE vtk
3416  USE polygon
3417  USE dashboard
3418  USE stl
3419 
3420 
3421  IMPLICIT NONE
3422 !-----------------------------------------------
3423 ! G l o b a l P a r a m e t e r s
3424 !-----------------------------------------------
3425 !-----------------------------------------------
3426 ! L o c a l P a r a m e t e r s
3427 !-----------------------------------------------
3428 !-----------------------------------------------
3429 ! L o c a l V a r i a b l e s
3430 !-----------------------------------------------
3431  INTEGER :: I,J,K,TOTAL_NUC,IDEAL_NCPP
3432 
3433  INTEGER, DIMENSION(0:DIM_I) :: NUC_I
3434  INTEGER, DIMENSION(0:DIM_J) :: NUC_J
3435  INTEGER, DIMENSION(0:DIM_K) :: NUC_K
3436 
3437 
3438  INTEGER :: ilistsize,jlistsize,klistsize ! size of list of cells
3439  INTEGER, ALLOCATABLE, DIMENSION(:) :: ALL_NUC_I ! Number of Useful Cells at I for all processors
3440  ! (I will repeat if decomposing in J or K direction)
3441  INTEGER, ALLOCATABLE, DIMENSION(:) :: ALL_LIST_I ! List of I for all processors
3442  INTEGER, ALLOCATABLE, DIMENSION(:) :: GLOBAL_NUC_I ! Number of Useful Cells at Global I
3443 
3444  INTEGER, ALLOCATABLE, DIMENSION(:) :: ALL_NUC_J ! Number of Useful Cells at J for all processors
3445  ! (J will repeat if decomposing in I or K direction)
3446  INTEGER, ALLOCATABLE, DIMENSION(:) :: ALL_LIST_J ! List of J for all processors
3447  INTEGER, ALLOCATABLE, DIMENSION(:) :: GLOBAL_NUC_J ! Number of Useful Cells at Global J
3448 
3449  INTEGER, ALLOCATABLE, DIMENSION(:) :: ALL_NUC_K ! Number of Useful Cells at K for all processors
3450  ! (K will repeat if decomposing in I or J direction)
3451  INTEGER, ALLOCATABLE, DIMENSION(:) :: ALL_LIST_K ! List of K for all processors
3452  INTEGER, ALLOCATABLE, DIMENSION(:) :: GLOBAL_NUC_K ! Number of Useful Cells at Global K
3453 
3454 
3455  INTEGER :: IPROC,PSUM
3456 
3457  INTEGER, DIMENSION(0:numPEs-1) :: NCPP_OLD,NCPP,NCPP_OLD_WITH_GHOST,NCPP_WITH_GHOST
3458 
3459  INTEGER, DIMENSION(0:NODESI-1) :: ISIZE_OLD
3460  INTEGER, DIMENSION(0:NODESJ-1) :: JSIZE_OLD
3461  INTEGER, DIMENSION(0:NODESK-1) :: KSIZE_OLD
3462 
3463  INTEGER :: JSIZE, IREMAIN,ISIZE, JREMAIN,KSIZE, KREMAIN
3464  DOUBLE PRECISION :: LIP_OLD
3465 
3466  INTEGER :: I_OFFSET,J_OFFSET,K_OFFSET,IERR
3467 
3468  INTEGER, DIMENSION(0:numPEs-1) :: disp,rcount
3469 
3470  INTEGER :: IPROC_OF_MAX_OLD,IPROC_OF_MIN_OLD
3471 
3472 !-----------------------------------------------
3473 !
3474  IF(.NOT.report_best_domain_size) RETURN
3475 
3476  IF(.NOT.cartesian_grid) RETURN ! Perform adjustement only when both CG
3477 ! IF(.NOT.ADJUST_PROC_DOMAIN_SIZE) RETURN ! and domain adjustment
3478  IF(nodesi*nodesj*nodesk==1) RETURN ! and parallel run are active
3479 
3480  IF(.not.allocated(isize_all)) allocate( isize_all(0:nodesi-1))
3481  IF(.not.allocated(jsize_all)) allocate( jsize_all(0:nodesj-1))
3482  IF(.not.allocated(ksize_all)) allocate( ksize_all(0:nodesk-1))
3483 
3484  isize_all(0:nodesi-1) = imax1-imin1+1 ! Assign default sizes in I, J and K-direction
3485  jsize_all(0:nodesj-1) = jmax1-jmin1+1
3486  ksize_all(0:nodesk-1) = kmax1-kmin1+1
3487 
3488 
3489  IF(nodesi>1) THEN ! DOMAIN DECOMPOSITION IN I-DIRECTION
3490 
3491 ! Assign size in I-direction
3492 
3493  DO i = istart1,iend1
3494  nuc_i(i) = 0 ! NUC : Number of Useful Cells
3495  DO j = jstart3,jend3
3496  DO k = kstart3,kend3
3497  IF( .NOT.dead_cell_at(i,j,k)) THEN ! Count number of useful cells
3498  nuc_i(i) = nuc_i(i) + 1
3499  ENDIF
3500  ENDDO
3501  ENDDO
3502  ENDDO
3503 
3504 ! Gather NUC onto the head node
3505 ! rcount is the size for each processor
3506 ! disp is the cumulative sum for each processor
3507 
3508  CALL allgather_1i (iend1-istart1+1,rcount,ierr)
3509 
3510  IF (mype == 0) THEN
3511  i_offset = 0
3512  ELSE
3513  i_offset = 0
3514  DO iproc=0,mype-1
3515  i_offset = i_offset + rcount(iproc)
3516  ENDDO
3517  ENDIF
3518 
3519  CALL allgather_1i (i_offset,disp,ierr)
3520 
3521  ilistsize=sum(rcount)
3522 
3523  allocate( all_nuc_i(ilistsize))
3524  allocate( all_list_i(ilistsize))
3525  allocate( global_nuc_i(imin1:imax1))
3526 
3527 ! Gather list of I and NUC, each processor has its own list
3528  call gatherv_1i( (/(i,i=istart1,iend1)/), iend1-istart1+1, all_list_i(:), rcount, disp, pe_io, ierr )
3529  call gatherv_1i( nuc_i(istart1:iend1), iend1-istart1+1, all_nuc_i(:), rcount, disp, pe_io, ierr )
3530 
3531 ! Get the glocal NUC for each unique value of I
3532  IF (mype == 0) THEN
3533  global_nuc_i = 0
3534  DO i=1,ilistsize
3535  global_nuc_i(all_list_i(i)) = global_nuc_i(all_list_i(i)) + all_nuc_i(i)
3536  ENDDO
3537  ENDIF
3538 
3539  ENDIF ! NODESI
3540 
3541  IF(nodesj>1) THEN ! DOMAIN DECOMPOSITION IN J-DIRECTION
3542 
3543 ! Assign size in J-direction
3544 
3545  DO j = jstart1,jend1
3546  nuc_j(j) = 0 ! NUC : Number of Useful Cells
3547  DO i = istart3,iend3
3548  DO k = kstart3,kend3
3549  IF( .NOT.dead_cell_at(i,j,k)) THEN ! Count number of useful cells
3550  nuc_j(j) = nuc_j(j) + 1
3551  ENDIF
3552  ENDDO
3553  ENDDO
3554  ENDDO
3555 
3556 ! Gather NUC onto the head node
3557 ! rcount is the size for each processor
3558 ! disp is the cumulative sum for each processor
3559 
3560  CALL allgather_1i (jend1-jstart1+1,rcount,ierr)
3561 
3562  IF (mype == 0) THEN
3563  j_offset = 0
3564  ELSE
3565  j_offset = 0
3566  DO iproc=0,mype-1
3567  j_offset = j_offset + rcount(iproc)
3568  ENDDO
3569  ENDIF
3570 
3571  CALL allgather_1i (j_offset,disp,ierr)
3572 
3573  jlistsize=sum(rcount)
3574 
3575  allocate( all_nuc_j(jlistsize))
3576  allocate( all_list_j(jlistsize))
3577  allocate( global_nuc_j(jmin1:jmax1))
3578 
3579 ! Gather list of J and NUC, each processor has its own list
3580  call gatherv_1i( (/(j,j=jstart1,jend1)/), jend1-jstart1+1, all_list_j(:), rcount, disp, pe_io, ierr )
3581  call gatherv_1i( nuc_j(jstart1:jend1), jend1-jstart1+1, all_nuc_j(:), rcount, disp, pe_io, ierr )
3582 
3583 ! Get the glocal NUC for each unique value of J
3584  IF (mype == 0) THEN
3585  global_nuc_j = 0
3586  DO j=1,jlistsize
3587  global_nuc_j(all_list_j(j)) = global_nuc_j(all_list_j(j)) + all_nuc_j(j)
3588  ENDDO
3589  ENDIF
3590 
3591  ENDIF ! NODESJ
3592 
3593  IF(nodesk>1) THEN ! DOMAIN DECOMPOSITION IN K-DIRECTION
3594 
3595 ! Assign size in K-direction
3596 
3597  DO k = kstart1,kend1
3598  nuc_k(k) = 0 ! NUC : Number of Useful Cells
3599  DO i = istart3,iend3
3600  DO j = jstart3,jend3
3601  IF( .NOT.dead_cell_at(i,j,k)) THEN ! Count number of useful cells
3602  nuc_k(k) = nuc_k(k) + 1
3603  ENDIF
3604  ENDDO
3605  ENDDO
3606  ENDDO
3607 
3608 ! Gather NUC onto the head node
3609 ! rcount is the size for each processor
3610 ! disp is the cumulative sum for each processor
3611 
3612  CALL allgather_1i (kend1-kstart1+1,rcount,ierr)
3613 
3614  IF (mype == 0) THEN
3615  k_offset = 0
3616  ELSE
3617  k_offset = 0
3618  DO iproc=0,mype-1
3619  k_offset = k_offset + rcount(iproc)
3620  ENDDO
3621  ENDIF
3622 
3623  CALL allgather_1i (k_offset,disp,ierr)
3624 
3625  klistsize=sum(rcount)
3626 
3627  allocate( all_nuc_k(klistsize))
3628  allocate( all_list_k(klistsize))
3629  allocate( global_nuc_k(kmin1:kmax1))
3630 
3631 ! Gather list of K and NUC, each processor has its own list
3632  call gatherv_1i( (/(k,k=kstart1,kend1)/), kend1-kstart1+1, all_list_k(:), rcount, disp, pe_io, ierr )
3633  call gatherv_1i( nuc_k(kstart1:kend1), kend1-kstart1+1, all_nuc_k(:), rcount, disp, pe_io, ierr )
3634 
3635 ! Get the glocal NUC for each unique value of K
3636  IF (mype == 0) THEN
3637  global_nuc_k = 0
3638  DO k=1,klistsize
3639  global_nuc_k(all_list_k(k)) = global_nuc_k(all_list_k(k)) + all_nuc_k(k)
3640  ENDDO
3641  ENDIF
3642 
3643  ENDIF ! NODESK
3644 
3645 
3646 
3647 ! DETERMINE BEST DOMAIN DECOMPOSITION FROM HEAD NODE
3648 
3649  IF (mype == pe_io) THEN
3650 
3651  IF(nodesi>1) THEN ! DOMAIN DECOMPOSITION IN I-DIRECTION
3652 
3653 ! For comparison with old grid size, determine the size in i direction (without adjustment) and add the remainder sequentially
3654 
3655  isize = (imax1-imin1+1)/nodesi
3656  isize_old(0:nodesi-1) = isize
3657 
3658  iremain = (imax1-imin1+1) - nodesi*isize
3659  IF (iremain.GE.1) THEN
3660  isize_old( 0:(iremain-1) ) = isize + 1
3661  ENDIF
3662 
3663  isize_all = isize_old
3664 
3665 ! Get load imbalance before optimization
3666  CALL get_lip_with_ghost_layers(nodesi,global_nuc_i(imin1:imax1),imin1, &
3667  imax1,isize_all,ncpp_old,ncpp_old_with_ghost, &
3668  lip_old,iproc_of_max_old,iproc_of_min_old)
3669 
3670  total_nuc = sum(ncpp_old_with_ghost(0:nodesi-1))
3671  ideal_ncpp = total_nuc / nodesi
3672  WRITE(*,1000)'INFO FROM REPORT_BEST_IJK_SIZE:'
3673  WRITE(*,1000)'ATTEMPTING TO REPORT BEST DOMAIN SIZE IN I-DIRECTION ...'
3674  WRITE (*, 1010) 'TOTAL NUMBER OF USEFUL CELLS = ',total_nuc
3675  WRITE (*, 1010) 'IDEAL NUMBER OF CELLS/I-NODE = ',ideal_ncpp
3676  WRITE (*, 1000) 'BEFORE OPTIMIZATION:'
3677  WRITE (*, 1000) '================================================='
3678  WRITE (*, 1000) ' I-NODE I-SIZE CELLS/NODE DIFF. (%)'
3679  WRITE (*, 1000) '================================================='
3680  DO iproc = 0,nodesi-1
3681  WRITE (*, 1020) iproc,isize_all(iproc),ncpp_old_with_ghost(iproc), &
3682  dble(ncpp_old_with_ghost(iproc)-ideal_ncpp)/dble(ideal_ncpp)*100.0d0
3683  ENDDO
3684  WRITE (*, 1000) '================================================='
3685 
3686 ! Optimize load imbalance
3687  CALL minimize_load_imbalance(nodesi,global_nuc_i(imin1:imax1),imin1,imax1,isize_all,ncpp,ncpp_with_ghost)
3688 
3689  total_nuc = sum(ncpp_with_ghost(0:nodesi-1))
3690  ideal_ncpp = total_nuc / nodesi
3691  WRITE (*, 1010) 'TOTAL NUMBER OF USEFUL CELLS = ',total_nuc
3692  WRITE (*, 1010) 'IDEAL NUMBER OF CELLS/I-NODE = ',ideal_ncpp
3693  WRITE (*, 1010) 'SINCE GHOST CELLS ARE INCLUDED, THE TOTALS'
3694  WRITE (*, 1010) 'BEFORE AND AFTER OPTIMIZATION MAY NOT MATCH.'
3695  WRITE (*, 1000) 'AFTER OPTIMIZATION:'
3696  WRITE (*, 1000) '================================================='
3697  WRITE (*, 1000) ' I-NODE I-SIZE CELLS/NODE DIFF. (%)'
3698  WRITE (*, 1000) '================================================='
3699  DO iproc = 0,nodesi-1
3700  WRITE (*, 1020) iproc,isize_all(iproc),ncpp_with_ghost(iproc), &
3701  dble(ncpp_with_ghost(iproc)-ideal_ncpp)/dble(ideal_ncpp)*100.0d0
3702  ENDDO
3703  WRITE (*, 1000) '================================================='
3704 
3705 ! Verify that the sum of all ISIZE_ALL matches IMAX
3706  psum = 0
3707  DO iproc = 0,nodesi-1
3708  psum = psum + isize_all(iproc)
3709  IF(isize_all(iproc)<5) THEN
3710  WRITE (*, 1010) 'ERROR: I-SIZE TOO SMALL FOR I-NODE:',iproc
3711  WRITE (*, 1010) 'I-SIZE = ',isize_all(iproc)
3712  CALL mfix_exit(mype)
3713  ENDIF
3714  ENDDO
3715 
3716  IF(psum/=imax) THEN
3717  WRITE (*, 1000) 'ERROR IN ADJUST_IJK_SIZE: UNABLE TO ASSIGN ISIZE TO PROCESSORS.'
3718  WRITE (*, 1000) 'SUM OF ISIZE_ALL DOES NOT MATCH IMAX:'
3719  WRITE (*, 1010) 'SUM OF ISIZE_ALL = ',psum
3720  WRITE (*, 1010) 'IMAX = ',imax
3721  CALL mfix_exit(mype)
3722  ENDIF
3723 
3724  ENDIF ! DOMAIN DECOMPOSITION IN I-DIRECTION
3725 
3726 
3727  IF(nodesj>1) THEN ! DOMAIN DECOMPOSITION IN J-DIRECTION
3728 ! For comparison with old grid size, determine the size in i direction (without adjustment) and add the remainder sequentially
3729 
3730  jsize = (jmax1-jmin1+1)/nodesj
3731  jsize_old(0:nodesj-1) = jsize
3732 
3733  jremain = (jmax1-jmin1+1) - nodesj*jsize
3734  IF (jremain.GE.1) THEN
3735  jsize_old( 0:(jremain-1) ) = jsize + 1
3736  ENDIF
3737 
3738  jsize_all = jsize_old
3739 
3740 ! Get load imbalance before optimization
3741  CALL get_lip_with_ghost_layers(nodesj,global_nuc_j(jmin1:jmax1),jmin1,jmax1,jsize_all, &
3742  ncpp_old,ncpp_old_with_ghost,lip_old,iproc_of_max_old,iproc_of_min_old)
3743 
3744  total_nuc = sum(ncpp_old_with_ghost(0:nodesj-1))
3745  ideal_ncpp = total_nuc / nodesj
3746  WRITE(*,1000)'INFO FROM REPORT_BEST_IJK_SIZE:'
3747  WRITE(*,1000)'ATTEMPTING TO REPORT BEST DOMAIN SIZE IN J-DIRECTION ...'
3748  WRITE (*, 1010) 'TOTAL NUMBER OF USEFUL CELLS = ',total_nuc
3749  WRITE (*, 1010) 'IDEAL NUMBER OF CELLS/J-NODE = ',ideal_ncpp
3750  WRITE (*, 1000) 'BEFORE OPTIMIZATION:'
3751  WRITE (*, 1000) '================================================='
3752  WRITE (*, 1000) ' J-NODE J-SIZE CELLS/NODE DIFF. (%)'
3753  WRITE (*, 1000) '================================================='
3754  DO iproc = 0,nodesj-1
3755  WRITE (*, 1020) iproc,jsize_all(iproc),ncpp_old_with_ghost(iproc), &
3756  dble(ncpp_old_with_ghost(iproc)-ideal_ncpp)/dble(ideal_ncpp)*100.0d0
3757  ENDDO
3758  WRITE (*, 1000) '================================================='
3759 
3760 ! Optimize load imbalance
3761  CALL minimize_load_imbalance(nodesj,global_nuc_j(jmin1:jmax1),jmin1,jmax1,jsize_all,ncpp,ncpp_with_ghost)
3762 
3763  total_nuc = sum(ncpp_with_ghost(0:nodesj-1))
3764  ideal_ncpp = total_nuc / nodesj
3765  WRITE (*, 1010) 'TOTAL NUMBER OF USEFUL CELLS = ',total_nuc
3766  WRITE (*, 1010) 'IDEAL NUMBER OF CELLS/J-NODE = ',ideal_ncpp
3767  WRITE (*, 1010) 'SINCE GHOST CELLS ARE INCLUDED, THE TOTALS'
3768  WRITE (*, 1010) 'BEFORE AND AFTER OPTIMIZATION MAY NOT MATCH.'
3769  WRITE (*, 1000) 'AFTER OPTIMIZATION:'
3770  WRITE (*, 1000) '================================================='
3771  WRITE (*, 1000) ' J-NODE J-SIZE CELLS/NODE DIFF. (%)'
3772  WRITE (*, 1000) '================================================='
3773  DO iproc = 0,nodesj-1
3774  WRITE (*, 1020) iproc,jsize_all(iproc),ncpp_with_ghost(iproc), &
3775  dble(ncpp_with_ghost(iproc)-ideal_ncpp)/dble(ideal_ncpp)*100.0d0
3776  ENDDO
3777  WRITE (*, 1000) '================================================='
3778 
3779 ! Verify that the sum of all JSIZE_ALL matches JMAX
3780  psum = 0
3781  DO iproc = 0,nodesj-1
3782  psum = psum + jsize_all(iproc)
3783  IF(jsize_all(iproc)<5) THEN
3784  WRITE (*, 1010) 'ERROR: J-SIZE TOO SMALL FOR J-NODE:',iproc
3785  WRITE (*, 1010) 'J-SIZE = ',jsize_all(iproc)
3786  CALL mfix_exit(mype)
3787  ENDIF
3788  ENDDO
3789 
3790  IF(psum/=jmax) THEN
3791  WRITE (*, 1000) 'ERROR IN ADJUST_IJK_SIZE: UNABLE TO ASSIGN JSIZE TO PROCESSORS.'
3792  WRITE (*, 1000) 'SUM OF JSIZE_ALL DOES NOT MATCH JMAX:'
3793  WRITE (*, 1010) 'SUM OF JSIZE_ALL = ',psum
3794  WRITE (*, 1010) 'JMAX = ',jmax
3795  CALL mfix_exit(mype)
3796  ENDIF
3797 
3798  ENDIF ! DOMAIN DECOMPOSITION IN J-DIRECTION
3799 
3800  IF(nodesk>1) THEN ! DOMAIN DECOMPOSITION IN K-DIRECTION
3801 ! For comparison with old grid size, determine the size in i direction (without adjustment) and add the remainder sequentially
3802 
3803  ksize = (kmax1-kmin1+1)/nodesk
3804  ksize_old(0:nodesk-1) = ksize
3805 
3806  kremain = (kmax1-kmin1+1) - nodesk*ksize
3807  IF (kremain.GE.1) THEN
3808  ksize_old( 0:(kremain-1) ) = ksize + 1
3809  ENDIF
3810 
3811  ksize_all = ksize_old
3812 
3813 ! Get load imbalance before optimization
3814  CALL get_lip_with_ghost_layers(nodesk,global_nuc_k(kmin1:kmax1),kmin1,kmax1,ksize_all, &
3815  ncpp_old,ncpp_old_with_ghost,lip_old,iproc_of_max_old,iproc_of_min_old)
3816 
3817  total_nuc = sum(ncpp_old_with_ghost(0:nodesk-1))
3818  ideal_ncpp = total_nuc / nodesk
3819  WRITE(*,1000)'INFO FROM REPORT_BEST_IJK_SIZE:'
3820  WRITE(*,1000)'ATTEMPTING TO REPORT BEST DOMAIN SIZE IN K-DIRECTION ...'
3821  WRITE (*, 1010) 'TOTAL NUMBER OF USEFUL CELLS = ',total_nuc
3822  WRITE (*, 1010) 'IDEAL NUMBER OF CELLS/K_NODE = ',ideal_ncpp
3823  WRITE (*, 1000) 'BEFORE OPTIMIZATION:'
3824  WRITE (*, 1000) '================================================='
3825  WRITE (*, 1000) ' K-NODE K-SIZE CELLS/NODE DIFF. (%)'
3826  WRITE (*, 1000) '================================================='
3827  DO iproc = 0,nodesk-1
3828  WRITE (*, 1020) iproc,ksize_all(iproc),ncpp_old_with_ghost(iproc), &
3829  dble(ncpp_old_with_ghost(iproc)-ideal_ncpp)/dble(ideal_ncpp)*100.0d0
3830  ENDDO
3831  WRITE (*, 1000) '================================================='
3832 
3833 ! Optimize load imbalance
3834  CALL minimize_load_imbalance(nodesk,global_nuc_k(kmin1:kmax1),kmin1,kmax1,ksize_all,ncpp,ncpp_with_ghost)
3835 
3836  total_nuc = sum(ncpp_with_ghost(0:nodesk-1))
3837  ideal_ncpp = total_nuc / nodesk
3838  WRITE (*, 1010) 'TOTAL NUMBER OF USEFUL CELLS = ',total_nuc
3839  WRITE (*, 1010) 'IDEAL NUMBER OF CELLS/K_NODE = ',ideal_ncpp
3840  WRITE (*, 1010) 'SINCE GHOST CELLS ARE INCLUDED, THE TOTALS'
3841  WRITE (*, 1010) 'BEFORE AND AFTER OPTIMIZATION MAY NOT MATCH.'
3842  WRITE (*, 1000) 'AFTER OPTIMIZATION:'
3843  WRITE (*, 1000) '================================================='
3844  WRITE (*, 1000) ' K-NODE K-SIZE CELLS/NODE DIFF. (%)'
3845  WRITE (*, 1000) '================================================='
3846  DO iproc = 0,nodesk-1
3847  WRITE (*, 1020) iproc,ksize_all(iproc),ncpp_with_ghost(iproc), &
3848  dble(ncpp_with_ghost(iproc)-ideal_ncpp)/dble(ideal_ncpp)*100.0d0
3849  ENDDO
3850  WRITE (*, 1000) '================================================='
3851 
3852 ! Verify that the sum of all KSIZE_ALL matches KMAX
3853  psum = 0
3854  DO iproc = 0,nodesk-1
3855  psum = psum + ksize_all(iproc)
3856  IF(ksize_all(iproc)<5) THEN
3857  WRITE (*, 1010) 'ERROR: K-SIZE TOO SMALL FOR K-NODE:',iproc
3858  WRITE (*, 1010) 'K-SIZE = ',ksize_all(iproc)
3859  CALL mfix_exit(mype)
3860  ENDIF
3861  ENDDO
3862 
3863  IF(psum/=kmax) THEN
3864  WRITE (*, 1000) 'ERROR IN ADJUST_IJK_SIZE: UNABLE TO ASSIGN KSIZE TO PROCESSORS.'
3865  WRITE (*, 1000) 'SUM OF KSIZE_ALL DOES NOT MATCH KMAX:'
3866  WRITE (*, 1010) 'SUM OF KSIZE_ALL = ',psum
3867  WRITE (*, 1010) 'KMAX = ',kmax
3868  CALL mfix_exit(mype)
3869  ENDIF
3870 
3871 
3872  ENDIF ! DOMAIN DECOMPOSITION IN K-DIRECTION
3873 
3874 
3875  OPEN(convert='BIG_ENDIAN',unit=777, file='suggested_gridmap.dat')
3876  WRITE (777, 1005) nodesi,nodesj,nodesk, ' ! NODESI, NODESJ, NODESK'
3877  DO iproc = 0,nodesi-1
3878  WRITE(777,1060) iproc,isize_all(iproc)
3879  ENDDO
3880  DO iproc = 0,nodesj-1
3881  WRITE(777,1060) iproc,jsize_all(iproc)
3882  ENDDO
3883  DO iproc = 0,nodesk-1
3884  WRITE(777,1060) iproc,ksize_all(iproc)
3885  ENDDO
3886 
3887  CLOSE(777)
3888  WRITE (*, 1000) '================================================='
3889  WRITE (*, 1000) 'GRID PARTITION SAVED IN FILE: suggested_gridmap.dat'
3890  WRITE (*, 1000) 'TO USE THIS DISTRIBUTION, RENAME THE FILE AS: gridmap.dat'
3891  WRITE (*, 1000) 'AND RUN MFIX AGAIN.'
3892 ! WRITE (*, 1000) 'MFIX WILL STOP NOW.'
3893  WRITE (*, 1000) '================================================='
3894 
3895 
3896  ENDIF ! (MyPE==PE_IO)
3897 
3898 ! Finalize and terminate MPI
3899  call parallel_fin
3900 
3901  stop
3902 
3903 ! Broadcast Domain sizes to all processors
3904 
3905 ! CALL BCAST(ISIZE_ALL)
3906 ! CALL BCAST(JSIZE_ALL)
3907 ! CALL BCAST(KSIZE_ALL)
3908 
3909 ! DOMAIN_SIZE_ADJUSTED = .TRUE.
3910 
3911 1000 FORMAT(1x,a)
3912 1005 FORMAT(1x,i10,i10,i10,a)
3913 1010 FORMAT(1x,a,i10,i10)
3914 1020 FORMAT(1x,i8,2(i12),f12.2)
3915 1030 FORMAT(1x,a,2(f10.1))
3916 1040 FORMAT(f10.1)
3917 1050 FORMAT(1x,3(a))
3918 1060 FORMAT(1x,i10,i10)
3919 
3920  RETURN
3921 
3922  END SUBROUTINE report_best_ijk_size
3923 
3924 
3925 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
3926 ! C
3927 ! Module name: GET_LIP_WITH_GHOST_LAYERS C
3928 ! Purpose: Compute Load Imbalance percentage C
3929 ! by including size of ghost layers C
3930 ! C
3931 ! C
3932 ! Author: Jeff Dietiker Date: 02-Dec-10 C
3933 ! Reviewer: Date: C
3934 ! C
3935 ! Revision Number # Date: ##-###-## C
3936 ! Author: # C
3937 ! Purpose: # C
3938 ! C
3939 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
3940 !
3941  SUBROUTINE get_lip_with_ghost_layers(NODESL,NUC_L,LMIN1,LMAX1,L_SIZE,NCPP,NCPP_WITH_GHOST,LIP,IPROC_OF_MAX,IPROC_OF_MIN)
3942 !
3943 !-----------------------------------------------
3944 ! M o d u l e s
3945 !-----------------------------------------------
3946  USE gridmap
3947  USE param
3948  USE param1
3949  USE constant
3950  USE run
3951  USE physprop
3952  USE indices
3953  USE scalars
3954  USE funits
3955  USE leqsol
3956  USE compar
3957  USE mpi_utility
3958  USE bc
3959  USE discretelement
3960 
3961  USE cutcell
3962  USE quadric
3963  USE vtk
3964  USE polygon
3965  USE dashboard
3966  USE stl
3967 
3968 
3969  IMPLICIT NONE
3970 !-----------------------------------------------
3971 ! G l o b a l P a r a m e t e r s
3972 !-----------------------------------------------
3973 !-----------------------------------------------
3974 ! L o c a l P a r a m e t e r s
3975 !-----------------------------------------------
3976 !-----------------------------------------------
3977 ! L o c a l V a r i a b l e s
3978 !-----------------------------------------------
3979  INTEGER :: NODESL,L,LMIN1,LMAX1,TOTAL_NUC,TOTAL_NUC_WITH_GHOST,IPROC_OF_MAX,IPROC_OF_MIN
3980 
3981  INTEGER :: LCOUNT1,LCOUNT2,MINVAL_NCPP,MAXVAL_NCPP,IDEAL_NCPP
3982  INTEGER, DIMENSION(LMIN1:LMAX1) :: NUC_L
3983 
3984  INTEGER :: IPROC
3985 
3986  INTEGER, DIMENSION(0:NODESL-1) :: NCPP,NCPP_WITH_GHOST,L_SIZE,L1,L2
3987 
3988 
3989  DOUBLE PRECISION :: LIP
3990 
3991 !-----------------------------------------------
3992 
3993  lcount1 = lmax1 - lmin1 + 1
3994  lcount2 = sum(l_size(0:nodesl-1))
3995 
3996  IF(lcount1/=lcount2) THEN
3997  WRITE(*,*)' ERROR: SUM OF CELLS DO NOT MATCH:',lcount1,lcount2
3998  CALL mfix_exit(mype)
3999  ENDIF
4000 
4001  l1(0) = lmin1
4002  l2(0) = l1(0) + l_size(0) - 1
4003 
4004  DO iproc = 1,nodesl-1
4005  l1(iproc) = l2(iproc-1) + 1
4006  l2(iproc) = l1(iproc) + l_size(iproc) - 1
4007  ENDDO
4008 
4009  DO iproc = 0,nodesl-1
4010  ncpp(iproc) = sum(nuc_l(l1(iproc):l2(iproc)))
4011 ! print*,'NUC=',NUC_L(L1(IPROC):L2(IPROC))
4012 ! print*,'L1,L2=',IPROC,L1(IPROC),L2(IPROC),NCPP(IPROC)
4013  ENDDO
4014 
4015  total_nuc = 0
4016 
4017  DO l=lmin1,lmax1
4018  total_nuc = total_nuc + nuc_l(l)
4019  ENDDO
4020 
4021  ncpp_with_ghost(0) = ncpp(0) + 2*nuc_l(l1(0)) + nuc_l(l1(1)) + nuc_l(l1(1)+1)
4022 
4023  DO iproc = 1,nodesl-2
4024  ncpp_with_ghost(iproc) = ncpp(iproc) &
4025  + nuc_l(l2(iproc-1)) + nuc_l(l2(iproc-1)-1) &
4026  + nuc_l(l1(iproc+1)) + nuc_l(l1(iproc+1)+1)
4027  ENDDO
4028 
4029  ncpp_with_ghost(nodesl-1) = ncpp(nodesl-1) + 2*nuc_l(l2(nodesl-1)) + nuc_l(l2(nodesl-2)) + nuc_l(l2(nodesl-2)-1)
4030 
4031  total_nuc_with_ghost = 0
4032  DO iproc = 0,nodesl-1
4033 ! print*,'NCPP_WITH_GHOST=',IPROC,L_SIZE(IPROC),NCPP(IPROC),NCPP_WITH_GHOST(IPROC)
4034  total_nuc_with_ghost = total_nuc_with_ghost + ncpp_with_ghost(iproc)
4035  ENDDO
4036 
4037 
4038  ideal_ncpp = total_nuc_with_ghost / numpes
4039 
4040  maxval_ncpp = maxval(ncpp_with_ghost)
4041  minval_ncpp = minval(ncpp_with_ghost)
4042 
4043  lip = dble(maxval_ncpp-ideal_ncpp)/dble(ideal_ncpp)*100.0d0
4044 ! LIP = DBLE(MAXVAL_NCPP-MINVAL_NCPP)/DBLE(MINVAL_NCPP)*100.0D0
4045 
4046 
4047  iproc_of_max = maxloc(ncpp_with_ghost,1)-1
4048  iproc_of_min = minloc(ncpp_with_ghost,1)-1
4049 
4050 
4051 ! print*,'IPROC_OF_MIN=',IPROC_OF_MIN,MINVAL_NCPP
4052 ! print*,'IPROC_OF_MAX=',IPROC_OF_MAX,MAXVAL_NCPP
4053 
4054  RETURN
4055  END SUBROUTINE get_lip_with_ghost_layers
4056 
4057 
4058 
4059 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
4060 ! C
4061 ! Module name: MINIMIZE_LOAD_IMBALANCE C
4062 ! Purpose: Rearrange L_SIZE to minimize load imbalance C
4063 ! by including size of ghost layers C
4064 ! C
4065 ! C
4066 ! Author: Jeff Dietiker Date: 02-Dec-10 C
4067 ! Reviewer: Date: C
4068 ! C
4069 ! Revision Number # Date: ##-###-## C
4070 ! Author: # C
4071 ! Purpose: # C
4072 ! C
4073 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
4074 !
4075  SUBROUTINE minimize_load_imbalance(NODESL,NUC_L,LMIN1,LMAX1,L_SIZE,NCPP,NCPP_WITH_GHOST)
4076 !
4077 !-----------------------------------------------
4078 ! M o d u l e s
4079 !-----------------------------------------------
4080  USE gridmap
4081  USE param
4082  USE param1
4083  USE constant
4084  USE run
4085  USE physprop
4086  USE indices
4087  USE scalars
4088  USE funits
4089  USE leqsol
4090  USE compar
4091  USE mpi_utility
4092  USE bc
4093  USE discretelement
4094 
4095  USE cutcell
4096  USE quadric
4097  USE vtk
4098  USE polygon
4099  USE dashboard
4100  USE stl
4101 
4102 
4103  IMPLICIT NONE
4104 !-----------------------------------------------
4105 ! G l o b a l P a r a m e t e r s
4106 !-----------------------------------------------
4107 !-----------------------------------------------
4108 ! L o c a l P a r a m e t e r s
4109 !-----------------------------------------------
4110 !-----------------------------------------------
4111 ! L o c a l V a r i a b l e s
4112 !-----------------------------------------------
4113  INTEGER :: NODESL,LMIN1,LMAX1,IPROC_OF_MAX,IPROC_OF_MIN
4114 
4115  INTEGER, DIMENSION(LMIN1:LMAX1) :: NUC_L
4116 
4117  INTEGER :: NN,NOIMPROVEMENT
4118 
4119  INTEGER,PARAMETER :: NAMAX=10000 ! maximum number of adjustments, increase if optimized load is not reached
4120 
4121  INTEGER, DIMENSION(0:numPEs-1) :: NCPP,NCPP_WITH_GHOST,L_SIZE,BEST_L_SIZE,BEST_NCPP,BEST_NCPP_WITH_GHOST
4122 
4123  DOUBLE PRECISION :: LIP,BEST_LIP
4124 
4125 !-----------------------------------------------
4126 
4127 
4128 ! NA = NAMAX ! Number of adjustments
4129 
4130 ! Initial estimate of LIP
4131 
4132  CALL get_lip_with_ghost_layers(nodesl,nuc_l,lmin1,lmax1,l_size,ncpp,ncpp_with_ghost,lip,iproc_of_max,iproc_of_min)
4133 
4134  best_lip = lip
4135  best_l_size = l_size
4136 
4137 ! print*,'INITIAL ESTIMATE OF LIP:',LIP
4138 ! WRITE (*, 1000) '================================================='
4139 ! WRITE (*, 1010) 'AFTER STEP:',N
4140 ! WRITE (*, 1010) 'NOIMPROVEMENT=',NOIMPROVEMENT
4141 ! WRITE (*, 1000) ' PROCESSOR J-SIZE CELLS/PROC. ERROR (%)'
4142 ! WRITE (*, 1000) '================================================='
4143 ! DO IPROC = 0,numPEs-1
4144 ! WRITE (*, 1020) IPROC,L_SIZE(IPROC),NCPP_WITH_GHOST(IPROC)
4145 ! ENDDO
4146 ! WRITE (*, 1000) '================================================='
4147 
4148 ! print*,'MIN ',IPROC_OF_MIN,NCPP_WITH_GHOST(IPROC_OF_MIN)
4149 ! print*,'MAX ',IPROC_OF_MAX,NCPP_WITH_GHOST(IPROC_OF_MAX)
4150 
4151  print*,'ATTEMPTING TO OPTIMIZE LOAD BALANCE...'
4152 
4153  noimprovement=0
4154 
4155  DO nn = 1,namax
4156 
4157  l_size(iproc_of_max) = l_size(iproc_of_max) - 1
4158  l_size(iproc_of_min) = l_size(iproc_of_min) + 1
4159 
4160  CALL get_lip_with_ghost_layers(nodesl,nuc_l,lmin1,lmax1,l_size,ncpp,ncpp_with_ghost,lip,iproc_of_max,iproc_of_min)
4161 
4162 ! print*,'After adjustment of LIP:',N, LIP
4163 ! print*,'MIN ',IPROC_OF_MIN,NCPP_WITH_GHOST(IPROC_OF_MIN)
4164 ! print*,'MAX ',IPROC_OF_MAX,NCPP_WITH_GHOST(IPROC_OF_MAX)
4165 
4166  IF(lip<best_lip) THEN
4167  best_lip = lip
4168  best_l_size = l_size
4169  best_ncpp = ncpp
4170  best_ncpp_with_ghost = ncpp_with_ghost
4171  noimprovement=0
4172  ELSE
4173  noimprovement = noimprovement + 1
4174  ENDIF
4175 
4176 ! WRITE (*, 1000) '================================================='
4177 ! WRITE (*, 1010) 'AFTER STEP:',N
4178 ! WRITE (*, 1010) 'NOIMPROVEMENT=',NOIMPROVEMENT
4179 ! WRITE (*, 1000) ' PROCESSOR J-SIZE CELLS/PROC. ERROR (%)'
4180 ! WRITE (*, 1000) '================================================='
4181 ! DO IPROC = 0,numPEs-1
4182 ! WRITE (*, 1020) IPROC,L_SIZE(IPROC),NCPP_WITH_GHOST(IPROC)
4183 ! ENDDO
4184 
4185  IF(noimprovement==10) THEN
4186  WRITE (*, 1000) 'OPTIMIZED LOAD BALANCE REACHED.'
4187  EXIT
4188  ENDIF
4189 
4190  ENDDO
4191 
4192 ! print*,'Best LIP = ',BEST_LIP
4193  l_size = best_l_size
4194  ncpp = best_ncpp
4195  ncpp_with_ghost = best_ncpp_with_ghost
4196 
4197 
4198 ! WRITE (*, 1000) '================================================='
4199 ! WRITE (*, 1000) 'AFTER OPTIMIZATION'
4200 ! WRITE (*, 1000) ' PROCESSOR J-SIZE CELLS/PROC. ERROR (%)'
4201 ! WRITE (*, 1000) '================================================='
4202 ! DO IPROC = 0,numPEs-1
4203 ! WRITE (*, 1020) IPROC,L_SIZE(IPROC),NCPP_WITH_GHOST(IPROC)
4204 ! ENDDO
4205 
4206 1000 FORMAT(1x,a)
4207 1010 FORMAT(1x,a,i10,i10)
4208 1020 FORMAT(1x,i8,2(i12),f12.2)
4209 
4210  RETURN
4211  END SUBROUTINE minimize_load_imbalance
4212 
4213 
4214 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
4215 ! C
4216 ! Module name: REPORT_BEST_IJK_SIZE0 C
4217 ! Purpose: Adjust domain size of each processor, based on an C
4218 ! estimate on the total number of useful cells. C
4219 ! The objective is to reduce load imbalance. C
4220 ! This is the parallel version C
4221 ! C
4222 ! C
4223 ! Author: Jeff Dietiker Date: 02-Dec-10 C
4224 ! Reviewer: Date: C
4225 ! C
4226 ! Revision Number # Date: ##-###-## C
4227 ! Author: # C
4228 ! Purpose: # C
4229 ! C
4230 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
4231 !
4232  SUBROUTINE report_best_ijk_size0
4233 !
4234 !-----------------------------------------------
4235 ! M o d u l e s
4236 !-----------------------------------------------
4237  USE gridmap
4238  USE param
4239  USE param1
4240  USE constant
4241  USE run
4242  USE physprop
4243  USE indices
4244  USE scalars
4245  USE funits
4246  USE leqsol
4247  USE compar
4248  USE mpi_utility
4249  USE bc
4250  USE discretelement
4251 
4252  USE parallel
4253 
4254  USE cutcell
4255  USE quadric
4256  USE vtk
4257  USE polygon
4258  USE dashboard
4259  USE stl
4260 
4261 
4262  IMPLICIT NONE
4263 !-----------------------------------------------
4264 ! G l o b a l P a r a m e t e r s
4265 !-----------------------------------------------
4266 !-----------------------------------------------
4267 ! L o c a l P a r a m e t e r s
4268 !-----------------------------------------------
4269 !-----------------------------------------------
4270 ! L o c a l V a r i a b l e s
4271 !-----------------------------------------------
4272  INTEGER :: I,J,K,TOTAL_NUC,IDEAL_NCPP
4273 
4274  INTEGER, DIMENSION(0:DIM_I) :: NUC_I,GLOBAL_NUC_I
4275  INTEGER, DIMENSION(0:DIM_J) :: NUC_J,GLOBAL_NUC_J
4276  INTEGER, DIMENSION(0:DIM_K) :: NUC_K,GLOBAL_NUC_K
4277 
4278  INTEGER :: IPROC,PSUM
4279 
4280  INTEGER, DIMENSION(0:numPEs-1) :: NCPP_OLD,NCPP,NCPP_OLD_WITH_GHOST,NCPP_WITH_GHOST
4281 
4282  INTEGER, DIMENSION(0:NODESJ-1) :: JSIZE_OLD
4283 
4284  INTEGER :: JSIZE, JREMAIN
4285  INTEGER :: MAXVAL_NCPP_OLD,MINVAL_NCPP_OLD,MAXVAL_NCPP,MINVAL_NCPP
4286  INTEGER :: AVG_NCPP_OLD,AVG_NCPP
4287  DOUBLE PRECISION :: LIP_OLD,LIP,MAXSPEEDUP_OLD,MAXSPEEDUP,P
4288 
4289  INTEGER :: I_OFFSET,J_OFFSET,K_OFFSET,IERR
4290 
4291  INTEGER, DIMENSION(0:numPEs-1) :: disp,rcount
4292 
4293  INTEGER :: IPROC_OF_MAX_OLD,IPROC_OF_MIN_OLD
4294 
4295 !-----------------------------------------------
4296 !
4297  RETURN
4298 
4299  IF(.NOT.cartesian_grid) RETURN ! Perform adjustement only when both CG
4300 ! IF(.NOT.ADJUST_PROC_DOMAIN_SIZE) RETURN ! and domain adjustment
4301  IF(nodesi*nodesj*nodesk==1) RETURN ! and parallel run are active
4302 
4303 
4304 
4305  IF(.not.allocated(isize_all)) allocate( isize_all(0:nodesi-1))
4306  IF(.not.allocated(jsize_all)) allocate( jsize_all(0:nodesj-1))
4307  IF(.not.allocated(ksize_all)) allocate( ksize_all(0:nodesk-1))
4308 
4309 ! allocate( ISIZE_ALL(0:NODESI-1))
4310 ! allocate( JSIZE_ALL(0:NODESJ-1))
4311 ! allocate( KSIZE_ALL(0:NODESK-1))
4312 
4313 
4314 
4315 
4316  IF(nodesi>1) THEN ! DOMAIN DECOMPOSITION IN I-DIRECTION
4317 
4318  IF(mype == 0.AND.nodesj*nodesk/=1) THEN
4319  WRITE(*,*)'ERROR IN SUBROUTINE REPORT_BEST_IJK_SIZE.'
4320  WRITE(*,*)'ADJUSTMENT POSSIBLE ONLY FOR DOMAIN DECOMPOSITION In ONE DIRECTION.'
4321  WRITE(*,*)'NODESI,NODESJ,NODESK =',nodesi,nodesj,nodesk
4322  WRITE(*,*)'MFIX WILL EXIT NOW.'
4323  CALL mfix_exit(mype)
4324  ENDIF
4325 
4326 
4327  jsize_all(0:nodesj-1) = jmax1-jmin1+1 ! Assign size in J and K-direction
4328  ksize_all(0:nodesk-1) = kmax1-kmin1+1
4329 
4330 ! Assign size in I-direction
4331  IF (mype == 0) THEN
4332  WRITE(*,1000)'INFO FROM REPORT_BEST_IJK_SIZE:'
4333  WRITE(*,1000)'ATTEMPTING TO REPORT BEST DOMAIN SIZE IN I-DIRECTION ...'
4334  ENDIF
4335 
4336  DO i = istart1,iend1
4337 
4338  nuc_i(i) = 0 ! NUC : Number of Useful Cells
4339 
4340  DO j = jstart3,jend3
4341  DO k = kstart3,kend3
4342 
4343  IF( .NOT.dead_cell_at(i,j,k)) THEN ! Count number of useful cells
4344  nuc_i(i) = nuc_i(i) + 1
4345  ENDIF
4346 
4347  ENDDO
4348  ENDDO
4349 
4350  ENDDO
4351 
4352 ! Gather NUC onto the head node
4353 
4354  CALL allgather_1i (iend1-istart1+1,rcount,ierr)
4355 
4356  IF (mype == 0) THEN
4357  i_offset = 0
4358  ELSE
4359  i_offset = 0
4360  DO iproc=0,mype-1
4361  i_offset = i_offset + rcount(iproc)
4362  ENDDO
4363  ENDIF
4364 
4365  CALL allgather_1i (i_offset,disp,ierr)
4366 
4367  call gatherv_1i( nuc_i(istart1:iend1), iend1-istart1+1, global_nuc_i(imin1:imax1), rcount, disp, pe_io, ierr )
4368 
4369 
4370 
4371  ELSEIF(nodesj>1) THEN ! DOMAIN DECOMPOSITION IN J-DIRECTION
4372 
4373 
4374  IF(mype == 0.AND.nodesi*nodesk/=1) THEN
4375  WRITE(*,*)'ERROR IN SUBROUTINE REPORT_BEST_IJK_SIZE.'
4376  WRITE(*,*)'ADJUSTMENT POSSIBLE ONLY FOR DOMAIN DECOMPOSITION In ONE DIRECTION.'
4377  WRITE(*,*)'NODESI,NODESJ,NODESK =',nodesi,nodesj,nodesk
4378  WRITE(*,*)'MFIX WILL EXIT NOW.'
4379  CALL mfix_exit(mype)
4380  ENDIF
4381 
4382 
4383  isize_all(0:nodesi-1) = imax1-imin1+1 ! Assign size in I and K-direction
4384  ksize_all(0:nodesk-1) = kmax1-kmin1+1
4385 
4386 ! Assign size in J-direction
4387  IF (mype == 0) THEN
4388  WRITE(*,1000)'INFO FROM REPORT_BEST_IJK_SIZE:'
4389  WRITE(*,1000)'ATTEMPTING TO REPORT BEST DOMAIN SIZE IN J-DIRECTION ...'
4390  ENDIF
4391 
4392  DO j = jstart1,jend1
4393 
4394  nuc_j(j) = 0 ! NUC : Number of Useful Cells
4395 
4396  DO i = istart3,iend3
4397  DO k = kstart3,kend3
4398 
4399  IF( .NOT.dead_cell_at(i,j,k)) THEN ! Count number of useful cells
4400  nuc_j(j) = nuc_j(j) + 1
4401  ENDIF
4402 
4403  ENDDO
4404  ENDDO
4405 
4406  ENDDO
4407 
4408 ! Gather NUC onto the head node
4409 
4410  IF(is_serial) THEN
4411 
4412  DO j=jmin1,jmax1
4413  global_nuc_j(j) = nuc_j(j)
4414  ENDDO
4415 
4416 
4417  ELSE
4418 
4419  CALL allgather_1i (jend1-jstart1+1,rcount,ierr)
4420 
4421  IF (mype == 0) THEN
4422  j_offset = 0
4423  ELSE
4424  j_offset = 0
4425  DO iproc=0,mype-1
4426  j_offset = j_offset + rcount(iproc)
4427  ENDDO
4428  ENDIF
4429 
4430  CALL allgather_1i (j_offset,disp,ierr)
4431 
4432  call gatherv_1i( nuc_j(jstart1:jend1), jend1-jstart1+1, global_nuc_j(jmin1:jmax1), rcount, disp, pe_io, ierr )
4433 
4434  ENDIF
4435 
4436 ! IF (myPE == 0) THEN
4437 ! DO J=JMIN1,JMAX1
4438 ! print*,'J,NUC=',J,GLOBAL_NUC_J(J)
4439 ! ENDDO
4440 ! ENDIF
4441 
4442 
4443 
4444  ELSEIF(nodesk>1) THEN ! DOMAIN DECOMPOSITION IN K-DIRECTION
4445 
4446  IF(mype == 0.AND.nodesi*nodesj/=1) THEN
4447  WRITE(*,*)'ERROR IN SUBROUTINE REPORT_BEST_IJK_SIZE.'
4448  WRITE(*,*)'ADJUSTMENT POSSIBLE ONLY FOR DOMAIN DECOMPOSITION In ONE DIRECTION.'
4449  WRITE(*,*)'NODESI,NODESJ,NODESK =',nodesi,nodesj,nodesk
4450  WRITE(*,*)'MFIX WILL EXIT NOW.'
4451  CALL mfix_exit(mype)
4452  ENDIF
4453 
4454 
4455  isize_all(0:nodesi-1) = imax1-imin1+1 ! Assign size in I and J-direction
4456  jsize_all(0:nodesj-1) = jmax1-jmin1+1
4457 
4458 ! Assign size in K-direction
4459  IF (mype == 0) THEN
4460  WRITE(*,1000)'INFO FROM REPORT_BEST_IJK_SIZE:'
4461  WRITE(*,1000)'ATTEMPTING TO REPORT BEST DOMAIN SIZE IN K-DIRECTION ...'
4462  ENDIF
4463 
4464  DO k = kstart1,kend1
4465 
4466  nuc_k(k) = 0 ! NUC : Number of Useful Cells
4467 
4468  DO i = istart1,iend1
4469  DO j = jstart1,jend1
4470 
4471  IF( .NOT.dead_cell_at(i,j,k)) THEN ! Count number of useful cells
4472  nuc_k(k) = nuc_k(k) + 1
4473  ENDIF
4474 
4475  ENDDO
4476  ENDDO
4477 
4478  ENDDO
4479 
4480 ! Gather NUC onto the head node
4481 
4482  CALL allgather_1i (kend1-kstart1+1,rcount,ierr)
4483 
4484  IF (mype == 0) THEN
4485  k_offset = 0
4486  ELSE
4487  k_offset = 0
4488  DO iproc=0,mype-1
4489  k_offset = k_offset + rcount(iproc)
4490  ENDDO
4491  ENDIF
4492 
4493  CALL allgather_1i (k_offset,disp,ierr)
4494 
4495  call gatherv_1i( nuc_k(kstart1:kend1), kend1-kstart1+1, global_nuc_k(kmin1:kmax1), rcount, disp, pe_io, ierr )
4496 
4497 
4498 
4499  ENDIF !(NODESI,NODESJ, OR NODESK)
4500 
4501 
4502 
4503 
4504  IF (mype == pe_io) THEN ! DETERMINE BEST DOMAIN DECOMPOSITION FROM HEAD NODE
4505 
4506  IF(nodesi>1) THEN ! DOMAIN DECOMPOSITION IN I-DIRECTION
4507 
4508 
4509 
4510  ELSEIF(nodesj>1) THEN ! DOMAIN DECOMPOSITION IN J-DIRECTION
4511 
4512  ! For comparison with old grid size, determine the size in j direction (without adjustment) and add the remainder sequentially
4513 
4514  jsize = (jmax1-jmin1+1)/nodesj
4515  jsize_old(0:nodesj-1) = jsize
4516 
4517  jremain = (jmax1-jmin1+1) - nodesj*jsize
4518  IF (jremain.GE.1) THEN
4519  jsize_old( 0:(jremain-1) ) = jsize + 1
4520  ENDIF
4521 
4522 
4523 
4524  jsize_all = jsize_old
4525 
4526  CALL get_lip_with_ghost_layers(nodesj,global_nuc_j(jmin1:jmax1),jmin1,jmax1,jsize_all, &
4527  ncpp_old,ncpp_old_with_ghost,lip_old,iproc_of_max_old,iproc_of_min_old)
4528 
4529 
4530 ! print*,'INITIAL ESTIMATE OF LIP, before minimizing LIP:',LIP_OLD
4531 ! WRITE (*, 1000) '================================================='
4532 ! WRITE (*, 1000) ' PROCESSOR J-SIZE CELLS/PROC. ERROR (%)'
4533 ! WRITE (*, 1000) '================================================='
4534 ! DO IPROC = 0,numPEs-1
4535 ! WRITE (*, 1020) IPROC,JSIZE_ALL(IPROC),NCPP_OLD_WITH_GHOST(IPROC)
4536 ! ENDDO
4537 ! WRITE (*, 1000) '================================================='
4538 
4539  CALL minimize_load_imbalance(nodesj,global_nuc_j(jmin1:jmax1),jmin1,jmax1,jsize_all,ncpp,ncpp_with_ghost)
4540 
4541 
4542  total_nuc = sum(ncpp_with_ghost(0:numpes-1))
4543  ideal_ncpp = total_nuc / numpes
4544 
4545  WRITE (*, 1000) 'AFTER OPTIMIZATION:'
4546  WRITE (*, 1010) 'TOTAL NUMBER OF USEFUL CELLS = ',total_nuc
4547  WRITE (*, 1010) 'IDEAL NUMBER OF CELLS/PROC. = ',ideal_ncpp
4548  WRITE (*, 1000) 'ACTUALL CELL COUNT:'
4549  WRITE (*, 1000) '================================================='
4550  WRITE (*, 1000) ' PROCESSOR J-SIZE CELLS/PROC. DIFF. (%)'
4551  WRITE (*, 1000) '================================================='
4552  DO iproc = 0,numpes-1
4553  WRITE (*, 1020) iproc,jsize_all(iproc),ncpp_with_ghost(iproc), &
4554  dble(ncpp_with_ghost(iproc)-ideal_ncpp)/dble(ideal_ncpp)*100.0d0
4555  ENDDO
4556 
4557 
4558 ! DO IPROC = 0 ,numPEs-1
4559 ! NCPP_OLD(IPROC) = (imax3-imin3+1)*(JSIZE_ALL(IPROC)+4)
4560 ! IF(DO_K) NCPP_OLD(IPROC) = NCPP_OLD(IPROC)*(kmax3-kmin3+1)
4561 ! ENDDO
4562 
4563 
4564  ! Verify that the sum of all JSIZE_ALL matches JMAX
4565 
4566  psum = 0
4567  DO iproc = 0,numpes-1
4568  psum = psum + jsize_all(iproc)
4569  IF(jsize_all(iproc)<5) THEN
4570  WRITE (*, 1010) 'ERROR: J-SIZE TOO SMALL FOR PROCESSOR:',iproc
4571  WRITE (*, 1010) 'J-SIZE = ',jsize_all(iproc)
4572  CALL mfix_exit(mype)
4573  ENDIF
4574 
4575  ENDDO
4576 
4577  IF(psum/=jmax) THEN
4578  WRITE (*, 1000) 'ERROR IN ADJUST_IJK_SIZE: UNABLE TO ASSIGN JSIZE TO PROCESSORS.'
4579  WRITE (*, 1000) 'SUM OF JSIZE_ALL DOES NOT MATCH JMAX:'
4580  WRITE (*, 1010) 'SUM OF JSIZE_ALL = ',psum
4581  WRITE (*, 1010) 'JMAX1 = ',jmax
4582  CALL mfix_exit(mype)
4583  ENDIF
4584 
4585 
4586  OPEN(convert='BIG_ENDIAN',unit=777, file='suggested_gridmap.dat')
4587  WRITE (777, 1000) 'J-SIZE DISTRIBUTION'
4588  WRITE (777, 1010) 'NUMBER OF PROCESSORS = ',numpes
4589  WRITE (777, 1000) '================================================='
4590  WRITE (777, 1000) ' PROCESSOR J-SIZE'
4591  WRITE (777, 1000) '================================================='
4592 
4593  DO iproc = 0,numpes-1
4594  WRITE(777,1060) iproc,jsize_all(iproc)
4595  ENDDO
4596  CLOSE(777)
4597  WRITE (*, 1000) '================================================='
4598  WRITE (*, 1000) 'J-SIZE DISTRIBUTION SAVED IN FILE: suggested_gridmap.dat'
4599  WRITE (*, 1000) 'TO USE THIS DISTRIBUTION, RENAME THE FILE AS: gridmap.dat'
4600  WRITE (*, 1000) 'AND RUN MFIX AGAIN.'
4601  WRITE (*, 1000) '================================================='
4602 
4603 
4604  ELSEIF(nodesk>1) THEN ! DOMAIN DECOMPOSITION IN K-DIRECTION
4605 
4606 
4607 
4608 
4609  ENDIF !(NODESI,NODESJ, OR NODESK)
4610 
4611 
4612  maxval_ncpp_old = maxval(ncpp_old_with_ghost)
4613  minval_ncpp_old = minval(ncpp_old_with_ghost)
4614  avg_ncpp_old = sum(ncpp_old_with_ghost)/numpes
4615 
4616 ! LIP_OLD = DBLE(MAXVAL_NCPP_OLD-AVG_NCPP_OLD)/DBLE(AVG_NCPP_OLD)*100.0D0
4617  lip_old = dble(maxval_ncpp_old-minval_ncpp_old)/dble(minval_ncpp_old)*100.0d0
4618 
4619 ! P = DBLE(MAXVAL_NCPP_OLD)/DBLE(AVG_NCPP_OLD)
4620  p = dble(minval_ncpp_old)/dble(maxval_ncpp_old)
4621 
4622 ! MAXSPEEDUP_OLD = DBLE(NumPes)*(ONE-LIP_OLD/100.0D0)
4623  maxspeedup_old = one / ((one-p) + p/numpes)
4624 
4625  maxval_ncpp = maxval(ncpp_with_ghost)
4626  minval_ncpp = minval(ncpp_with_ghost)
4627  avg_ncpp = sum(ncpp_with_ghost(0:numpes-1))/numpes
4628 
4629 ! LIP = DBLE(MAXVAL_NCPP-AVG_NCPP)/DBLE(AVG_NCPP)*100.0D0
4630  lip = dble(maxval_ncpp-minval_ncpp)/dble(minval_ncpp)*100.0d0
4631 
4632 ! P = DBLE(MAXVAL_NCPP)/DBLE(AVG_NCPP)
4633  p = dble(minval_ncpp)/dble(maxval_ncpp)
4634 
4635 ! MAXSPEEDUP = DBLE(NumPes)*(ONE-LIP/100.0D0)
4636  maxspeedup = one / ((one-p) + p/numpes)
4637 
4638  WRITE (*, 1000) '================================================='
4639  WRITE (*, 1000) 'ESTIMATED PARALLEL LOAD BALANCING STATISTICS'
4640  WRITE (*, 1000) 'COMPARISION BETWEEN UNIFORM SIZE (OLD)'
4641  WRITE (*, 1000) 'AND SUGGESTED SIZE (NEW)'
4642  WRITE (*, 1000) '================================================='
4643  WRITE (*, 1000) ' OLD NEW'
4644  WRITE (*, 1010) 'MAX CELL COUNT : ',maxval_ncpp_old,maxval_ncpp
4645  WRITE (*, 1010) 'AT PROCESSOR : ',maxloc(ncpp_old_with_ghost)-1,maxloc(ncpp_with_ghost)-1
4646  WRITE (*, 1010) 'MIN CELL COUNT : ',minval_ncpp_old,minval_ncpp
4647  WRITE (*, 1010) 'AT PROCESSOR : ',minloc(ncpp_old_with_ghost)-1,minloc(ncpp_with_ghost)-1
4648  WRITE (*, 1010) 'AVG CELL COUNT : ',avg_ncpp_old,avg_ncpp
4649  WRITE (*, 1000) ''
4650  WRITE (*, 1030) 'LOAD IMBALANCE (%) : ',lip_old,lip
4651  WRITE (*, 1000) ''
4652 ! WRITE (*, 1030) 'IDEAL SPEEDUP : ',DBLE(NumPEs),DBLE(NumPEs)
4653 ! WRITE (*, 1030) 'MAX SPEEDUP : ',MAXSPEEDUP_OLD,MAXSPEEDUP
4654 ! WRITE (*, 1030) 'MAX EFFICIENCY (%) : ',MAXSPEEDUP_OLD/DBLE(NumPEs)*100.0,MAXSPEEDUP/DBLE(NumPEs)*100.0
4655 
4656  WRITE (*, 1000) '================================================='
4657 
4658 
4659 
4660  ENDIF ! (MyPE==PE_IO)
4661 
4662 
4663 
4664 ! Broadcast Domain sizes to all processors
4665 
4666 ! CALL BCAST(ISIZE_ALL)
4667 ! CALL BCAST(JSIZE_ALL)
4668 ! CALL BCAST(KSIZE_ALL)
4669 
4670 ! DOMAIN_SIZE_ADJUSTED = .TRUE.
4671 
4672 
4673 
4674 
4675 1000 FORMAT(1x,a)
4676 1010 FORMAT(1x,a,i10,i10)
4677 1020 FORMAT(1x,i8,2(i12),f12.2)
4678 1030 FORMAT(1x,a,2(f10.1))
4679 1040 FORMAT(f10.1)
4680 1050 FORMAT(1x,3(a))
4681 1060 FORMAT(1x,i8,i12)
4682 
4683  RETURN
4684 
4685  END SUBROUTINE report_best_ijk_size0
4686 
4687 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
4688 ! C
4689 ! Module name: GET_LIP_WITH_GHOST_LAYERS0 C
4690 ! Purpose: Compute Load Imbalance percentage C
4691 ! by including size of ghost layers C
4692 ! C
4693 ! C
4694 ! Author: Jeff Dietiker Date: 02-Dec-10 C
4695 ! Reviewer: Date: C
4696 ! C
4697 ! Revision Number # Date: ##-###-## C
4698 ! Author: # C
4699 ! Purpose: # C
4700 ! C
4701 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
4702 !
4703  SUBROUTINE get_lip_with_ghost_layers0(NODESL,NUC_L,LMIN1,LMAX1,L_SIZE,NCPP,NCPP_WITH_GHOST,LIP,IPROC_OF_MAX,IPROC_OF_MIN)
4704 !
4705 !-----------------------------------------------
4706 ! M o d u l e s
4707 !-----------------------------------------------
4708  USE gridmap
4709  USE param
4710  USE param1
4711  USE constant
4712  USE run
4713  USE physprop
4714  USE indices
4715  USE scalars
4716  USE funits
4717  USE leqsol
4718  USE compar
4719  USE mpi_utility
4720  USE bc
4721  USE discretelement
4722 
4723  USE cutcell
4724  USE quadric
4725  USE vtk
4726  USE polygon
4727  USE dashboard
4728  USE stl
4729 
4730 
4731  IMPLICIT NONE
4732 !-----------------------------------------------
4733 ! G l o b a l P a r a m e t e r s
4734 !-----------------------------------------------
4735 !-----------------------------------------------
4736 ! L o c a l P a r a m e t e r s
4737 !-----------------------------------------------
4738 !-----------------------------------------------
4739 ! L o c a l V a r i a b l e s
4740 !-----------------------------------------------
4741  INTEGER :: NODESL,L,LMIN1,LMAX1,TOTAL_NUC,TOTAL_NUC_WITH_GHOST,IPROC_OF_MAX,IPROC_OF_MIN
4742 
4743  INTEGER :: LCOUNT1,LCOUNT2,MINVAL_NCPP,MAXVAL_NCPP,IDEAL_NCPP
4744  INTEGER, DIMENSION(LMIN1:LMAX1) :: NUC_L
4745 
4746  INTEGER :: IPROC
4747 
4748  INTEGER, DIMENSION(0:NODESL-1) :: NCPP,NCPP_WITH_GHOST,L_SIZE,L1,L2
4749 
4750 
4751  DOUBLE PRECISION :: LIP
4752 
4753 !-----------------------------------------------
4754 
4755  lcount1 = lmax1 - lmin1 + 1
4756  lcount2 = sum(l_size(0:nodesl-1))
4757 
4758  IF(lcount1/=lcount2) THEN
4759  WRITE(*,*)' ERROR: SUM OF CELLS DO NOT MATCH:',lcount1,lcount2
4760  CALL mfix_exit(mype)
4761  ENDIF
4762 
4763  l1(0) = lmin1
4764  l2(0) = l1(0) + l_size(0) - 1
4765 
4766  DO iproc = 1,nodesl-1
4767  l1(iproc) = l2(iproc-1) + 1
4768  l2(iproc) = l1(iproc) + l_size(iproc) - 1
4769  ENDDO
4770 
4771  DO iproc = 0,nodesl-1
4772  ncpp(iproc) = sum(nuc_l(l1(iproc):l2(iproc)))
4773 ! print*,'NUC=',NUC_L(L1(IPROC):L2(IPROC))
4774 ! print*,'L1,L2=',IPROC,L1(IPROC),L2(IPROC),NCPP(IPROC)
4775  ENDDO
4776 
4777  total_nuc = 0
4778 
4779  DO l=lmin1,lmax1
4780  total_nuc = total_nuc + nuc_l(l)
4781  ENDDO
4782 
4783  ncpp_with_ghost(0) = ncpp(0) + 2*nuc_l(l1(0)) + nuc_l(l1(1)) + nuc_l(l1(1)+1)
4784 
4785  DO iproc = 1,nodesl-2
4786  ncpp_with_ghost(iproc) = ncpp(iproc) &
4787  + nuc_l(l2(iproc-1)) + nuc_l(l2(iproc-1)-1) &
4788  + nuc_l(l1(iproc+1)) + nuc_l(l1(iproc+1)+1)
4789  ENDDO
4790 
4791  ncpp_with_ghost(nodesl-1) = ncpp(nodesl-1) + 2*nuc_l(l2(nodesl-1)) + nuc_l(l2(nodesl-2)) + nuc_l(l2(nodesl-2)-1)
4792 
4793  total_nuc_with_ghost = 0
4794  DO iproc = 0,nodesl-1
4795 ! print*,'NCPP_WITH_GHOST=',IPROC,L_SIZE(IPROC),NCPP(IPROC),NCPP_WITH_GHOST(IPROC)
4796  total_nuc_with_ghost = total_nuc_with_ghost + ncpp_with_ghost(iproc)
4797  ENDDO
4798 
4799 
4800  ideal_ncpp = total_nuc_with_ghost / numpes
4801 
4802  maxval_ncpp = maxval(ncpp_with_ghost)
4803  minval_ncpp = minval(ncpp_with_ghost)
4804 
4805 ! LIP = DBLE(MAXVAL_NCPP-IDEAL_NCPP)/DBLE(IDEAL_NCPP)*100.0D0
4806  lip = dble(maxval_ncpp-minval_ncpp)/dble(minval_ncpp)*100.0d0
4807 
4808 
4809  iproc_of_max = maxloc(ncpp_with_ghost,1)-1
4810  iproc_of_min = minloc(ncpp_with_ghost,1)-1
4811 
4812 
4813 ! print*,'IPROC_OF_MIN=',IPROC_OF_MIN,MINVAL_NCPP
4814 ! print*,'IPROC_OF_MAX=',IPROC_OF_MAX,MAXVAL_NCPP
4815 
4816  RETURN
4817  END SUBROUTINE get_lip_with_ghost_layers0
4818 
4819 
4820 
4821 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
4822 ! C
4823 ! Module name: MINIMIZE_LOAD_IMBALANCE0 C
4824 ! Purpose: Rearrange L_SIZE to minimize load imbalance C
4825 ! by including size of ghost layers C
4826 ! C
4827 ! C
4828 ! Author: Jeff Dietiker Date: 02-Dec-10 C
4829 ! Reviewer: Date: C
4830 ! C
4831 ! Revision Number # Date: ##-###-## C
4832 ! Author: # C
4833 ! Purpose: # C
4834 ! C
4835 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
4836 !
4837  SUBROUTINE minimize_load_imbalance0(NODESL,NUC_L,LMIN1,LMAX1,L_SIZE,NCPP,NCPP_WITH_GHOST)
4838 !
4839 !-----------------------------------------------
4840 ! M o d u l e s
4841 !-----------------------------------------------
4842  USE gridmap
4843  USE param
4844  USE param1
4845  USE constant
4846  USE run
4847  USE physprop
4848  USE indices
4849  USE scalars
4850  USE funits
4851  USE leqsol
4852  USE compar
4853  USE mpi_utility
4854  USE bc
4855  USE discretelement
4856 
4857  USE cutcell
4858  USE quadric
4859  USE vtk
4860  USE polygon
4861  USE dashboard
4862  USE stl
4863 
4864 
4865  IMPLICIT NONE
4866 !-----------------------------------------------
4867 ! G l o b a l P a r a m e t e r s
4868 !-----------------------------------------------
4869 !-----------------------------------------------
4870 ! L o c a l P a r a m e t e r s
4871 !-----------------------------------------------
4872 !-----------------------------------------------
4873 ! L o c a l V a r i a b l e s
4874 !-----------------------------------------------
4875  INTEGER :: NODESL,LMIN1,LMAX1,IPROC_OF_MAX,IPROC_OF_MIN
4876 
4877  INTEGER, DIMENSION(LMIN1:LMAX1) :: NUC_L
4878 
4879  INTEGER :: NN,NOIMPROVEMENT
4880 
4881  INTEGER,PARAMETER :: NAMAX=10000 ! maximum number of adjustments, increase if optimized load is not reached
4882 
4883  INTEGER, DIMENSION(0:numPEs-1) :: NCPP,NCPP_WITH_GHOST,L_SIZE,BEST_L_SIZE,BEST_NCPP,BEST_NCPP_WITH_GHOST
4884 
4885 
4886  DOUBLE PRECISION :: LIP,BEST_LIP
4887 
4888 !-----------------------------------------------
4889 
4890 
4891 ! NA = NAMAX ! Number of adjustments
4892 
4893 ! Initial estimate of LIP
4894 
4895  CALL get_lip_with_ghost_layers0(nodesl,nuc_l,lmin1,lmax1,l_size,ncpp,ncpp_with_ghost,lip,iproc_of_max,iproc_of_min)
4896 
4897  best_lip = lip
4898  best_l_size = l_size
4899 
4900 ! print*,'INITIAL ESTIMATE OF LIP:',LIP
4901 ! WRITE (*, 1000) '================================================='
4902 ! WRITE (*, 1010) 'AFTER STEP:',N
4903 ! WRITE (*, 1010) 'NOIMPROVEMENT=',NOIMPROVEMENT
4904 ! WRITE (*, 1000) ' PROCESSOR J-SIZE CELLS/PROC. ERROR (%)'
4905 ! WRITE (*, 1000) '================================================='
4906 ! DO IPROC = 0,numPEs-1
4907 ! WRITE (*, 1020) IPROC,L_SIZE(IPROC),NCPP_WITH_GHOST(IPROC)
4908 ! ENDDO
4909 ! WRITE (*, 1000) '================================================='
4910 
4911 ! print*,'MIN ',IPROC_OF_MIN,NCPP_WITH_GHOST(IPROC_OF_MIN)
4912 ! print*,'MAX ',IPROC_OF_MAX,NCPP_WITH_GHOST(IPROC_OF_MAX)
4913 
4914  print*,'ATTEMPTING TO OPTIMIZE LOAD BALANCE...'
4915 
4916  noimprovement=0
4917 
4918  DO nn = 1,namax
4919 
4920  l_size(iproc_of_max) = l_size(iproc_of_max) - 1
4921  l_size(iproc_of_min) = l_size(iproc_of_min) + 1
4922 
4923  CALL get_lip_with_ghost_layers(nodesl,nuc_l,lmin1,lmax1,l_size,ncpp,ncpp_with_ghost,lip,iproc_of_max,iproc_of_min)
4924 
4925 ! print*,'After adjustment of LIP:',N, LIP
4926 ! print*,'MIN ',IPROC_OF_MIN,NCPP_WITH_GHOST(IPROC_OF_MIN)
4927 ! print*,'MAX ',IPROC_OF_MAX,NCPP_WITH_GHOST(IPROC_OF_MAX)
4928 
4929  IF(lip<best_lip) THEN
4930  best_lip = lip
4931  best_l_size = l_size
4932  best_ncpp = ncpp
4933  best_ncpp_with_ghost = ncpp_with_ghost
4934  noimprovement=0
4935  ELSE
4936  noimprovement = noimprovement + 1
4937  ENDIF
4938 
4939 ! WRITE (*, 1000) '================================================='
4940 ! WRITE (*, 1010) 'AFTER STEP:',N
4941 ! WRITE (*, 1010) 'NOIMPROVEMENT=',NOIMPROVEMENT
4942 ! WRITE (*, 1000) ' PROCESSOR J-SIZE CELLS/PROC. ERROR (%)'
4943 ! WRITE (*, 1000) '================================================='
4944 ! DO IPROC = 0,numPEs-1
4945 ! WRITE (*, 1020) IPROC,L_SIZE(IPROC),NCPP_WITH_GHOST(IPROC)
4946 ! ENDDO
4947 
4948  IF(noimprovement==10) THEN
4949  WRITE (*, 1000) 'OPTIMIZED LOAD BALANCE REACHED.'
4950  EXIT
4951  ENDIF
4952 
4953  ENDDO
4954 
4955 ! print*,'Best LIP = ',BEST_LIP
4956  l_size = best_l_size
4957  ncpp = best_ncpp
4958  ncpp_with_ghost = best_ncpp_with_ghost
4959 
4960 
4961 ! WRITE (*, 1000) '================================================='
4962 ! WRITE (*, 1000) 'AFTER OPTIMIZATION'
4963 ! WRITE (*, 1000) ' PROCESSOR J-SIZE CELLS/PROC. ERROR (%)'
4964 ! WRITE (*, 1000) '================================================='
4965 ! DO IPROC = 0,numPEs-1
4966 ! WRITE (*, 1020) IPROC,L_SIZE(IPROC),NCPP_WITH_GHOST(IPROC)
4967 ! ENDDO
4968 
4969 1000 FORMAT(1x,a)
4970 1010 FORMAT(1x,a,i10,i10)
4971 1020 FORMAT(1x,i8,2(i12),f12.2)
4972 
4973  RETURN
4974  END SUBROUTINE minimize_load_imbalance0
4975 END MODULE check_data_cg
cutcell::last_dy
double precision, dimension(max_cp) last_dy
Definition: cutcell_mod.f:482
mpi_utility
Definition: mpi_utility_mod.f:6
param::dimension_i
integer dimension_i
Definition: param_mod.f:7
compar::iend3
integer iend3
Definition: compar_mod.f:80
dashboard::dtmin
double precision dtmin
Definition: dashboard_mod.f:6
quadric::n_quadric
integer n_quadric
Definition: quadric_mod.f:6
cutcell::re_indexing
logical re_indexing
Definition: cutcell_mod.f:16
cutcell::last_dz
double precision, dimension(max_cp) last_dz
Definition: cutcell_mod.f:483
dashboard::n_dashboard
integer n_dashboard
Definition: dashboard_mod.f:4
funits::dmp_log
logical dmp_log
Definition: funits_mod.f:6
cutcell::tol_small_area
double precision tol_small_area
Definition: cutcell_mod.f:384
bc::bc_volflow_g
double precision, dimension(dimension_bc) bc_volflow_g
Definition: bc_mod.f:195
check_data_cg::report_best_processor_size
subroutine report_best_processor_size
Definition: check_data_cartesian.f:3291
compar::jstart3
integer jstart3
Definition: compar_mod.f:80
cutcell::ncz
integer, dimension(max_cp) ncz
Definition: cutcell_mod.f:480
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
compar::kend1
integer kend1
Definition: compar_mod.f:80
compar::istart1
integer istart1
Definition: compar_mod.f:80
bc::bc_t_g
double precision, dimension(dimension_bc) bc_t_g
Definition: bc_mod.f:97
cutcell::partial_check_03
logical partial_check_03
Definition: cutcell_mod.f:462
toleranc::compare
logical function compare(V1, V2)
Definition: toleranc_mod.f:94
dashboard::f_dashboard
integer f_dashboard
Definition: dashboard_mod.f:4
quadric::reactor1_y1
double precision, dimension(dim_quadric) reactor1_y1
Definition: quadric_mod.f:42
compar::iend1
integer iend1
Definition: compar_mod.f:80
quadric::itermax_int
integer itermax_int
Definition: quadric_mod.f:78
quadric::group_size
integer, dimension(dim_group) group_size
Definition: quadric_mod.f:68
param1::one
double precision, parameter one
Definition: param1_mod.f:29
cutcell::alpha_max
double precision alpha_max
Definition: cutcell_mod.f:386
compar::short_gridmap_init
logical short_gridmap_init
Definition: compar_mod.f:164
check_data_cg::get_lip_with_ghost_layers
subroutine get_lip_with_ghost_layers(NODESL, NUC_L, LMIN1, LMAX1, L_SIZE, NCPP, NCPP_WITH_GHOST, LIP, IPROC_OF_MAX, IPROC_OF_MIN)
Definition: check_data_cartesian.f:3942
quadric::c2c_r2
double precision, dimension(dim_quadric) c2c_r2
Definition: quadric_mod.f:36
param::dimension_k
integer dimension_k
Definition: param_mod.f:9
quadric::c2c_r1
double precision, dimension(dim_quadric) c2c_r1
Definition: quadric_mod.f:36
quadric::lambda_x
double precision, dimension(dim_quadric) lambda_x
Definition: quadric_mod.f:14
cutcell::fac_dim_max_cut_cell
double precision fac_dim_max_cut_cell
Definition: cutcell_mod.f:28
quadric::torus_r2
double precision, dimension(dim_quadric) torus_r2
Definition: quadric_mod.f:24
cutcell::ncx
integer, dimension(max_cp) ncx
Definition: cutcell_mod.f:480
quadric::radius
double precision, dimension(dim_quadric) radius
Definition: quadric_mod.f:22
run::species_eq
logical, dimension(0:dim_m) species_eq
Definition: run_mod.f:115
bc::bc_w_s
double precision, dimension(dimension_bc, dim_m) bc_w_s
Definition: bc_mod.f:129
cutcell::cg_safe_mode
integer, dimension(10) cg_safe_mode
Definition: cutcell_mod.f:415
cutcell::adjust_proc_domain_size
logical adjust_proc_domain_size
Definition: cutcell_mod.f:19
bc::bc_x_s
double precision, dimension(dimension_bc, dim_m, dim_n_s) bc_x_s
Definition: bc_mod.f:254
quadric::clip_xmax
double precision, dimension(dim_quadric) clip_xmax
Definition: quadric_mod.f:56
cutcell::cpz
double precision, dimension(0:max_cp) cpz
Definition: cutcell_mod.f:478
cutcell::first_dz
double precision, dimension(max_cp) first_dz
Definition: cutcell_mod.f:483
check_data_cg::check_data_cartesian
subroutine check_data_cartesian
Definition: check_data_cartesian.f:20
check_data_cg::minimize_load_imbalance
subroutine minimize_load_imbalance(NODESL, NUC_L, LMIN1, LMAX1, L_SIZE, NCPP, NCPP_WITH_GHOST)
Definition: check_data_cartesian.f:4076
compar::kstart3
integer kstart3
Definition: compar_mod.f:80
mpi_utility::allgather_1i
subroutine allgather_1i(lbuf, gbuf, idebug)
Definition: mpi_utility_mod.f:146
param::dimension_bc
integer, parameter dimension_bc
Definition: param_mod.f:61
vtk
Definition: vtk_mod.f:1
bc::bc_type_enum
integer, dimension(dimension_bc) bc_type_enum
Definition: bc_mod.f:146
quadric::t_x
double precision, dimension(dim_quadric) t_x
Definition: quadric_mod.f:18
cutcell::nodesj_report
integer nodesj_report
Definition: cutcell_mod.f:490
dashboard::nit_max
integer nit_max
Definition: dashboard_mod.f:5
param1::undefined
double precision, parameter undefined
Definition: param1_mod.f:18
quadric::relation_with_previous
character(len=9), dimension(dim_group) relation_with_previous
Definition: quadric_mod.f:74
quadric::ucoil_y1
double precision, dimension(dim_quadric) ucoil_y1
Definition: quadric_mod.f:28
cutcell::cpy
double precision, dimension(0:max_cp) cpy
Definition: cutcell_mod.f:478
scalars::a
double precision, dimension(:), allocatable a
Definition: scalars_mod.f:29
bc::bc_v_g
double precision, dimension(dimension_bc) bc_v_g
Definition: bc_mod.f:117
bc::bc_velmag_s
double precision, dimension(dimension_bc, dim_m) bc_velmag_s
Definition: bc_mod.f:140
compar::kstart1
integer kstart1
Definition: compar_mod.f:80
compar::kend3
integer kend3
Definition: compar_mod.f:80
quadric::ucoil_r2
double precision, dimension(dim_quadric) ucoil_r2
Definition: quadric_mod.f:26
quadric::n_z
double precision, dimension(dim_quadric) n_z
Definition: quadric_mod.f:50
cutcell::cg_ur_fac
double precision, dimension(9) cg_ur_fac
Definition: cutcell_mod.f:441
quadric::quadric_scale
double precision quadric_scale
Definition: quadric_mod.f:12
quadric::dim_group
integer, parameter dim_group
Definition: quadric_mod.f:64
compar::numpes
integer numpes
Definition: compar_mod.f:24
dashboard::smmax
double precision smmax
Definition: dashboard_mod.f:7
quadric::n_group
integer n_group
Definition: quadric_mod.f:66
quadric::reactor1_r1
double precision, dimension(dim_quadric) reactor1_r1
Definition: quadric_mod.f:40
error_manager::init_err_msg
subroutine init_err_msg(CALLER)
Definition: error_manager_mod.f:171
quadric::c2c_y2
double precision, dimension(dim_quadric) c2c_y2
Definition: quadric_mod.f:34
dashboard::smmin
double precision smmin
Definition: dashboard_mod.f:7
compar::pe_io
integer pe_io
Definition: compar_mod.f:30
dashboard::init_time
double precision init_time
Definition: dashboard_mod.f:7
cutcell::erz
double precision, dimension(max_cp) erz
Definition: cutcell_mod.f:479
bc::bc_volflow_s
double precision, dimension(dimension_bc, dim_m) bc_volflow_s
Definition: bc_mod.f:198
dashboard::dtmax
double precision dtmax
Definition: dashboard_mod.f:6
quadric::torus_r1
double precision, dimension(dim_quadric) torus_r1
Definition: quadric_mod.f:24
quadric::bc_id_q
integer, dimension(dim_quadric) bc_id_q
Definition: quadric_mod.f:62
check_data_cg::get_dxyz_from_control_points
subroutine get_dxyz_from_control_points
Definition: check_data_cartesian.f:1888
cutcell::nodesi_report
integer nodesi_report
Definition: cutcell_mod.f:490
cutcell::bc_u_id
integer, dimension(:), allocatable bc_u_id
Definition: cutcell_mod.f:434
cutcell::cpx
double precision, dimension(0:max_cp) cpx
Definition: cutcell_mod.f:478
exit::mfix_exit
subroutine mfix_exit(myID, normal_termination)
Definition: exit.f:5
cutcell::use_stl
logical use_stl
Definition: cutcell_mod.f:428
physprop::mmax
integer mmax
Definition: physprop_mod.f:19
cutcell::tol_small_cell
double precision tol_small_cell
Definition: cutcell_mod.f:383
compar::dead_cell_at
logical, dimension(:,:,:), allocatable dead_cell_at
Definition: compar_mod.f:127
bc::bc_t_s
double precision, dimension(dimension_bc, dim_m) bc_t_s
Definition: bc_mod.f:101
physprop::ro_g0
double precision ro_g0
Definition: physprop_mod.f:59
get_poly_data
subroutine get_poly_data
Definition: get_poly_data.f:26
compar::jend3
integer jend3
Definition: compar_mod.f:80
stl
Definition: stl_mod.f:1
cutcell::bc_w_id
integer, dimension(:), allocatable bc_w_id
Definition: cutcell_mod.f:436
exit
Definition: exit.f:2
cutcell::nodesk_report
integer nodesk_report
Definition: cutcell_mod.f:490
quadric::ucoil_r1
double precision, dimension(dim_quadric) ucoil_r1
Definition: quadric_mod.f:26
quadric::bend_r2
double precision, dimension(dim_quadric) bend_r2
Definition: quadric_mod.f:30
check_data_cg::adjust_ijk_size
subroutine adjust_ijk_size
Definition: check_data_cartesian.f:2704
bc::bc_scalar
double precision, dimension(dimension_bc, dim_scalar) bc_scalar
Definition: bc_mod.f:384
f
double precision function f(POS, ALPHAC, D_Target, L, N)
Definition: check_data_cartesian.f:2653
quadric::bend_r1
double precision, dimension(dim_quadric) bend_r1
Definition: quadric_mod.f:30
cut_cell_preproc::print_cg_header
subroutine print_cg_header
Definition: cut_cell_preprocessing.f:781
quadric::tol_f
double precision tol_f
Definition: quadric_mod.f:76
cutcell::pg_option
integer pg_option
Definition: cutcell_mod.f:393
bc::bc_p_g
double precision, dimension(dimension_bc) bc_p_g
Definition: bc_mod.f:80
gridmap::gridmap_init
subroutine gridmap_init
Definition: gridmap_mod.f:170
cutcell::n_usr_def
integer n_usr_def
Definition: cutcell_mod.f:424
quadric::reactor1_theta1
double precision, dimension(dim_quadric) reactor1_theta1
Definition: quadric_mod.f:48
cutcell::tol_snap
double precision, dimension(3) tol_snap
Definition: cutcell_mod.f:374
quadric::half_angle
double precision, dimension(dim_quadric) half_angle
Definition: quadric_mod.f:38
parallel::is_serial
logical is_serial
Definition: parallel_mod.f:11
quadric::dquadric
double precision, dimension(dim_quadric) dquadric
Definition: quadric_mod.f:16
cutcell::ncy
integer, dimension(max_cp) ncy
Definition: cutcell_mod.f:480
param1::half
double precision, parameter half
Definition: param1_mod.f:28
funits::unit_log
integer, parameter unit_log
Definition: funits_mod.f:21
run
Definition: run_mod.f:13
cutcell::erx
double precision, dimension(max_cp) erx
Definition: cutcell_mod.f:479
quadric::ucoil_y2
double precision, dimension(dim_quadric) ucoil_y2
Definition: quadric_mod.f:28
compar::jsize_all
integer, dimension(:), allocatable jsize_all
Definition: compar_mod.f:153
cutcell::ery
double precision, dimension(max_cp) ery
Definition: cutcell_mod.f:479
param1::large_number
double precision, parameter large_number
Definition: param1_mod.f:23
bc::bc_velmag_g
double precision, dimension(dimension_bc) bc_velmag_g
Definition: bc_mod.f:136
cutcell::bc_v_id
integer, dimension(:), allocatable bc_v_id
Definition: cutcell_mod.f:435
param
Definition: param_mod.f:2
quadric::lambda_y
double precision, dimension(dim_quadric) lambda_y
Definition: quadric_mod.f:14
define_quadrics
subroutine define_quadrics
Definition: define_quadrics.f:15
cutcell::cartesian_grid
logical cartesian_grid
Definition: cutcell_mod.f:13
bc::bc_v_s
double precision, dimension(dimension_bc, dim_m) bc_v_s
Definition: bc_mod.f:121
bc::bc_massflow_g
double precision, dimension(dimension_bc) bc_massflow_g
Definition: bc_mod.f:201
check_data_cg::find_cell_ratio
subroutine find_cell_ratio(POS, D_Target, L, NN, ALPHA3)
Definition: check_data_cartesian.f:2534
physprop::nmax
integer, dimension(0:dim_m) nmax
Definition: physprop_mod.f:119
bc::bc_massflow_s
double precision, dimension(dimension_bc, dim_m) bc_massflow_s
Definition: bc_mod.f:204
compar::isize_all
integer, dimension(:), allocatable isize_all
Definition: compar_mod.f:153
run::k_epsilon
logical k_epsilon
Definition: run_mod.f:97
compar::mype
integer mype
Definition: compar_mod.f:24
bc::bc_u_g
double precision, dimension(dimension_bc) bc_u_g
Definition: bc_mod.f:109
check_data_cg::report_best_ijk_size0
subroutine report_best_ijk_size0
Definition: check_data_cartesian.f:4233
compar::nodesj
integer nodesj
Definition: compar_mod.f:37
run::energy_eq
logical energy_eq
Definition: run_mod.f:100
quadric::quadric_form
character(len=12), dimension(dim_quadric) quadric_form
Definition: quadric_mod.f:10
compar::domain_size_adjusted
logical domain_size_adjusted
Definition: compar_mod.f:155
compar::ijkstart3
integer ijkstart3
Definition: compar_mod.f:80
cutcell::use_msh
logical use_msh
Definition: cutcell_mod.f:430
quadric::n_y
double precision, dimension(dim_quadric) n_y
Definition: quadric_mod.f:50
param1::undefined_i
integer, parameter undefined_i
Definition: param1_mod.f:19
cutcell::last_dx
double precision, dimension(max_cp) last_dx
Definition: cutcell_mod.f:481
bc::bc_u_s
double precision, dimension(dimension_bc, dim_m) bc_u_s
Definition: bc_mod.f:113
polygon::n_polygon
integer n_polygon
Definition: polygon_mod.f:6
error_manager::err_msg
character(len=line_length), dimension(line_count) err_msg
Definition: error_manager_mod.f:29
quadric::group_relation
character(len=9), dimension(dim_group) group_relation
Definition: quadric_mod.f:72
iterate
Definition: iterate.f:2
cutcell::use_polygon
logical use_polygon
Definition: cutcell_mod.f:426
bc::bc_e_turb_g
double precision, dimension(dimension_bc) bc_e_turb_g
Definition: bc_mod.f:404
cutcell::bc_id
integer, dimension(:), allocatable bc_id
Definition: cutcell_mod.f:433
compar::nodesk
integer nodesk
Definition: compar_mod.f:37
scalars::nscalar
integer nscalar
Definition: scalars_mod.f:7
check_data_cg::get_lip_with_ghost_layers0
subroutine get_lip_with_ghost_layers0(NODESL, NUC_L, LMIN1, LMAX1, L_SIZE, NCPP, NCPP_WITH_GHOST, LIP, IPROC_OF_MAX, IPROC_OF_MIN)
Definition: check_data_cartesian.f:4704
quadric::clip_xmin
double precision, dimension(dim_quadric) clip_xmin
Definition: quadric_mod.f:56
cutcell::first_dx
double precision, dimension(max_cp) first_dx
Definition: cutcell_mod.f:481
bc::bc_theta_m
double precision, dimension(dimension_bc, dim_m) bc_theta_m
Definition: bc_mod.f:105
cutcell::tol_delh
double precision tol_delh
Definition: cutcell_mod.f:377
iterate::max_nit
integer max_nit
Definition: iterate.f:15
bc::bc_ep_g
double precision, dimension(dimension_bc) bc_ep_g
Definition: bc_mod.f:77
check_data_cg::minimize_load_imbalance0
subroutine minimize_load_imbalance0(NODESL, NUC_L, LMIN1, LMAX1, L_SIZE, NCPP, NCPP_WITH_GHOST)
Definition: check_data_cartesian.f:4838
bc::bc_k_turb_g
double precision, dimension(dimension_bc) bc_k_turb_g
Definition: bc_mod.f:403
compar::nodesi
integer nodesi
Definition: compar_mod.f:37
check_data_cg::check_bc_flags
subroutine check_bc_flags
Definition: check_data_cartesian.f:960
error_manager
Definition: error_manager_mod.f:8
physprop::ro_s0
double precision, dimension(dim_m) ro_s0
Definition: physprop_mod.f:28
quadric::dim_quadric
integer, parameter dim_quadric
Definition: quadric_mod.f:4
mpi_utility::bcast
Definition: mpi_utility_mod.f:53
quadric::reactor1_yr1
double precision, dimension(dim_quadric) reactor1_yr1
Definition: quadric_mod.f:44
cutcell::report_best_domain_size
logical report_best_domain_size
Definition: cutcell_mod.f:22
cut_cell_preproc
Definition: cut_cell_preprocessing.f:1
quadric::t_y
double precision, dimension(dim_quadric) t_y
Definition: quadric_mod.f:18
quadric::piece_xmin
double precision, dimension(dim_quadric) piece_xmin
Definition: quadric_mod.f:58
cut_cell_preproc::eval_f
subroutine eval_f(METHOD, x1, x2, x3, Q, f, CLIP_FLAG)
Definition: cut_cell_preprocessing.f:855
constant::pi
double precision, parameter pi
Definition: constant_mod.f:158
polygon::tol_poly
double precision tol_poly
Definition: polygon_mod.f:16
quadric::bend_theta1
double precision, dimension(dim_quadric) bend_theta1
Definition: quadric_mod.f:32
run::time
double precision time
Definition: run_mod.f:45
cutcell::max_cp
integer, parameter max_cp
Definition: cutcell_mod.f:476
quadric::t_z
double precision, dimension(dim_quadric) t_z
Definition: quadric_mod.f:18
compar::istart3
integer istart3
Definition: compar_mod.f:80
run::granular_energy
logical granular_energy
Definition: run_mod.f:112
bc::bc_w_g
double precision, dimension(dimension_bc) bc_w_g
Definition: bc_mod.f:125
bc::cg_mi_converted_to_ps
logical, dimension(dimension_bc) cg_mi_converted_to_ps
Definition: bc_mod.f:406
compar::ijksize3_all
integer, dimension(:), allocatable ijksize3_all
Definition: compar_mod.f:65
check_data_cg::build_cone_for_c2c
subroutine build_cone_for_c2c(Q)
Definition: check_data_cartesian.f:1476
quadric::n_x
double precision, dimension(dim_quadric) n_x
Definition: quadric_mod.f:50
check_data_cg
Definition: check_data_cartesian.f:1
quadric::c2c_y1
double precision, dimension(dim_quadric) c2c_y1
Definition: quadric_mod.f:34
quadric::group_q
integer, dimension(dim_group, dim_quadric) group_q
Definition: quadric_mod.f:70
compar::ncpp_uniform
integer, dimension(:), allocatable ncpp_uniform
Definition: compar_mod.f:157
dashboard::nit_min
integer nit_min
Definition: dashboard_mod.f:5
geometry::x
double precision, dimension(:), allocatable x
Definition: geometry_mod.f:129
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
compar::jend1
integer jend1
Definition: compar_mod.f:80
error_manager::flush_err_msg
subroutine flush_err_msg(DEBUG, HEADER, FOOTER, ABORT, LOG, CALL_TREE)
Definition: error_manager_mod.f:305
compar::jstart1
integer jstart1
Definition: compar_mod.f:80
bc::bc_rop_s
double precision, dimension(dimension_bc, dim_m) bc_rop_s
Definition: bc_mod.f:92
mpi_utility::gatherv_1i
subroutine gatherv_1i(lbuf, sendcnt, gbuf, rcount, disp, mroot, idebug)
Definition: mpi_utility_mod.f:204
bc
Definition: bc_mod.f:23
get_msh_data
subroutine get_msh_data
Definition: get_stl_data.f:28
param::dimension_j
integer dimension_j
Definition: param_mod.f:8
quadric::piece_xmax
double precision, dimension(dim_quadric) piece_xmax
Definition: quadric_mod.f:58
leqsol::leq_pc
character(len=4), dimension(dim_eqs) leq_pc
Definition: leqsol_mod.f:26
get_stl_data
subroutine get_stl_data
Definition: get_stl_data.f:657
cutcell::first_dy
double precision, dimension(max_cp) first_dy
Definition: cutcell_mod.f:482
check_data_cg::report_best_ijk_size
subroutine report_best_ijk_size
Definition: check_data_cartesian.f:3392