File: N:\mfix\model\cartesian_grid\check_data_cartesian.f
1 MODULE CHECK_DATA_CG
2 USE exit, only: mfix_exit
3
4 CONTAINS
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19 SUBROUTINE CHECK_DATA_CARTESIAN
20
21
22
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
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')
220
221 (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')
245
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')
262
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')
279
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')
297
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')
314
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')
331
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')
348
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')
363
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')
379
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')
397
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')
415
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')
433
434
435 CALL BUILD_CONE_FOR_C2C(Q)
436
437
438 CASE ('TORUS_INT','TORUS_EXT')
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')
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')
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')
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')
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')
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
651 (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
904
905 = 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
946
947
948
949
950
951
952
953
954
955
956
957
958
959 SUBROUTINE CHECK_BC_FLAGS
960
961
962
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
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
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
1032
1033
1034
1035
1036
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
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
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
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178 = 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
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
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
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475 SUBROUTINE BUILD_CONE_FOR_C2C(Q)
1476
1477
1478
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
1505
1506
1507
1508
1509
1510
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 = 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
1531
1532 (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
1586
1587 (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
1642
1643 (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
1696
1697 (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
1752
1753 (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
1806
1807 (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
1873
1874
1875
1876
1877
1878
1879
1880
1881
1882
1883
1884
1885
1886
1887 SUBROUTINE GET_DXYZ_FROM_CONTROL_POINTS
1888
1889
1890
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
1917
1918
1919
1920
1921
1922
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
1936
1937
1938
1939
1940
1941
1942 DO NN = MAX_CP,1,-1
1943 CPX(nn) = CPX(NN-1)
1944 ENDDO
1945
1946 CPX(0) = ZERO
1947
1948
1949
1950 = 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
1967
1968
1969
1970 = .TRUE.
1971
1972 DO NN = 1,NX
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
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)
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)
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
2048
2049
2050
2051
2052 = 0
2053 = 0
2054
2055 DO NN = 1,NX
2056
2057 = I1 + NCX(nn) - 1
2058
2059 IF(INDEPENDENT_SEGMENT(nn)) THEN
2060
2061 L = CPX(nn) - CPX(NN-1)
2062
2063 IF(ERX(nn)/=ONE) THEN
2064 CELL_RATIO = ERX(nn)**(ONE/DBLE(NCX(nn)-1))
2065 (I1) = L * (ONE - CELL_RATIO) / (ONE - CELL_RATIO**NCX(nn))
2066
2067 DO I = I1+1,I2
2068 (I) = DX(I-1) * CELL_RATIO
2069 ENDDO
2070
2071 ELSE
2072 DX(I1:I2) = L / NCX(nn)
2073 ENDIF
2074
2075 ENDIF
2076
2077 I1 = I2 + 1
2078
2079 ENDDO
2080
2081
2082
2083
2084 = 0
2085 = 0
2086
2087 DO NN = 1,NX
2088
2089 = I1 + NCX(nn) - 1
2090
2091 IF(.NOT.INDEPENDENT_SEGMENT(nn)) THEN
2092
2093 L = CPX(nn) - CPX(NN-1)
2094
2095 IF(FIRST_DX(nn)<ZERO) THEN
2096 DX(I1) = DX(I1-1)
2097 CALL FIND_CELL_RATIO('FIRST',DX(I1),L,NCX(NN),CELL_RATIO)
2098 DO I = I1+1,I2
2099 (I) = DX(I-1) * CELL_RATIO
2100 ENDDO
2101 ELSEIF(LAST_DX(nn)<ZERO) THEN
2102 DX(I2) = DX(I2+1)
2103 CALL FIND_CELL_RATIO('LAST ',DX(I2),L,NCX(nn),CELL_RATIO)
2104 DO I = I2-1,I1,-1
2105 (I) = DX(I+1) / CELL_RATIO
2106 ENDDO
2107 ENDIF
2108
2109 ENDIF
2110
2111 I1 = I2 + 1
2112
2113 ENDDO
2114
2115
2116
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
2128
2129
2130
2131
2132
2133
2134 DO NN = MAX_CP,1,-1
2135 CPY(nn) = CPY(NN-1)
2136 ENDDO
2137
2138 CPY(0) = ZERO
2139
2140
2141
2142 = 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
2159
2160
2161
2162 = .TRUE.
2163
2164 DO NN = 1,NY
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
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)
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)
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
2240
2241
2242
2243
2244 = 0
2245 = 0
2246
2247 DO NN = 1,NY
2248
2249 = J1 + NCY(nn) - 1
2250
2251 IF(INDEPENDENT_SEGMENT(nn)) THEN
2252
2253 L = CPY(nn) - CPY(NN-1)
2254
2255 IF(ERY(nn)/=ONE) THEN
2256 CELL_RATIO = ERY(nn)**(ONE/DBLE(NCY(nn)-1))
2257 (J1) = L * (ONE - CELL_RATIO) / (ONE - CELL_RATIO**NCY(nn))
2258
2259 DO J = J1+1,J2
2260 (J) = DY(J-1) * CELL_RATIO
2261 ENDDO
2262
2263 ELSE
2264 DY(J1:J2) = L / NCY(nn)
2265 ENDIF
2266
2267 ENDIF
2268
2269 J1 = J2 + 1
2270
2271 ENDDO
2272
2273
2274
2275
2276 = 0
2277 = 0
2278
2279 DO NN = 1,NY
2280
2281 = J1 + NCY(nn) - 1
2282
2283 IF(.NOT.INDEPENDENT_SEGMENT(nn)) THEN
2284
2285 L = CPY(nn) - CPY(NN-1)
2286
2287 IF(FIRST_DY(nn)<ZERO) THEN
2288 DY(J1) = DY(J1-1)
2289 CALL FIND_CELL_RATIO('FIRST',DY(J1),L,NCY(nn),CELL_RATIO)
2290 DO J = J1+1,J2
2291 (J) = DY(J-1) * CELL_RATIO
2292 ENDDO
2293 ELSEIF(LAST_DY(nn)<ZERO) THEN
2294 DY(J2) = DY(J2+1)
2295 CALL FIND_CELL_RATIO('LAST ',DY(J2),L,NCY(nn),CELL_RATIO)
2296 DO J = J2-1,J1,-1
2297 (J) = DY(J+1) / CELL_RATIO
2298 ENDDO
2299 ENDIF
2300
2301 ENDIF
2302
2303 J1 = J2 + 1
2304
2305 ENDDO
2306
2307
2308
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
2320
2321
2322 IF(NO_K) RETURN
2323
2324
2325
2326
2327
2328 DO NN = MAX_CP,1,-1
2329 CPZ(nn) = CPZ(NN-1)
2330 ENDDO
2331
2332 CPZ(0) = ZERO
2333
2334
2335
2336 = 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
2353
2354
2355
2356 = .TRUE.
2357
2358 DO NN = 1,NZ
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
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)
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)
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
2434
2435
2436
2437
2438 = 0
2439 = 0
2440
2441 DO NN = 1,NZ
2442
2443 = K1 + NCZ(nn) - 1
2444
2445 IF(INDEPENDENT_SEGMENT(nn)) THEN
2446
2447 L = CPZ(nn) - CPZ(NN-1)
2448
2449 IF(ERZ(nn)/=ONE) THEN
2450 CELL_RATIO = ERZ(nn)**(ONE/DBLE(NCZ(nn)-1))
2451 (K1) = L * (ONE - CELL_RATIO) / (ONE - CELL_RATIO**NCZ(nn))
2452
2453 DO K = K1+1,K2
2454 (K) = DZ(K-1) * CELL_RATIO
2455 ENDDO
2456
2457 ELSE
2458 DZ(K1:K2) = L / NCZ(nn)
2459 ENDIF
2460
2461 ENDIF
2462
2463 K1 = K2 + 1
2464
2465 ENDDO
2466
2467
2468
2469
2470 = 0
2471 = 0
2472
2473 DO NN = 1,NZ
2474
2475 = K1 + NCZ(nn) - 1
2476
2477 IF(.NOT.INDEPENDENT_SEGMENT(nn)) THEN
2478
2479 L = CPZ(nn) - CPZ(NN-1)
2480
2481 IF(FIRST_DZ(nn)<ZERO) THEN
2482 DZ(K1) = DZ(K1-1)
2483 CALL FIND_CELL_RATIO('FIRST',DZ(K1),L,NCZ(nn),CELL_RATIO)
2484 DO K = K1+1,K2
2485 (K) = DZ(K-1) * CELL_RATIO
2486 ENDDO
2487 ELSEIF(LAST_DZ(nn)<ZERO) THEN
2488 DZ(K2) = DZ(K2+1)
2489 CALL FIND_CELL_RATIO('LAST ',DZ(K2),L,NCZ(nn),CELL_RATIO)
2490 DO K = K2-1,K1,-1
2491 (K) = DZ(K+1) / CELL_RATIO
2492 ENDDO
2493 ENDIF
2494
2495 ENDIF
2496
2497 K1 = K2 + 1
2498
2499 ENDDO
2500
2501
2502
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
2518
2519
2520
2521
2522
2523
2524
2525
2526
2527
2528
2529
2530
2531
2532
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
2553 INTEGER :: NN,niter
2554 CHARACTER (LEN=5) :: POS
2555
2556 DU = L / DBLE(NN)
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
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
2592
2593
2594
2595 = 1
2596 SOLUTION_FOUND = .FALSE.
2597
2598 if(DABS(f1)<TOL_F) then
2599 = .TRUE.
2600 ALPHA3 = ALPHA1
2601 elseif(DABS(f2)<TOL_F) then
2602 = .TRUE.
2603 ALPHA3 = ALPHA2
2604 elseif(f1*f2 < ZERO) then
2605 = 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)
2610
2611 = F(POS,ALPHA3,D_Target,L,NN)
2612
2613 if(f1*f3<0) then
2614 = 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)
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
2686
2687
2688
2689
2690
2691
2692
2693
2694
2695
2696
2697
2698
2699
2700
2701
2702
2703 SUBROUTINE ADJUST_IJK_SIZE
2704
2705
2706
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
2733
2734
2735
2736
2737
2738
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
2773 IF(.NOT.ADJUST_PROC_DOMAIN_SIZE) RETURN
2774 IF(NODESI*NODESJ*NODESK==1) RETURN
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
2801
2802 = IMAX3
2803 DIMENSION_J = JMAX3
2804 DIMENSION_K = KMAX3
2805
2806 PARTIAL_CHECK_03 = .TRUE.
2807
2808 = .FALSE.
2809 CALL CHECK_DATA_CARTESIAN
2810
2811 CALL DEFINE_QUADRICS
2812
2813 SHIFT = .TRUE.
2814
2815 IF(SHIFT) THEN
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
2867
2868
2869
2870 IF(NODESI>1) THEN
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
2882 (0:NODESK-1) = kmax1-kmin1+1
2883
2884
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
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
2910
2911
2912
2913
2914
2915
2916 IF (F_COPY < TOL_F ) THEN
2917 (I) = NUC_I(I) + 1
2918 ENDIF
2919
2920 ENDDO
2921 ENDDO
2922
2923 ENDDO
2924
2925
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
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
2957 (0:NODESK-1) = kmax1-kmin1+1
2958
2959
2960
2961
2962
2963
2964
2965
2966
2967
2968
2969
2970
2971
2972
2973
2974
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
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
3000
3001
3002
3003
3004
3005
3006 IF (F_COPY < TOL_F ) THEN
3007 (J) = NUC_J(J) + 1
3008 ENDIF
3009
3010 ENDDO
3011 ENDDO
3012
3013 ENDDO
3014
3015
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
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
3045 (0:NODESJ-1) = jmax1-jmin1+1
3046
3047
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
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
3072
3073
3074
3075
3076
3077
3078 IF (F_COPY < TOL_F ) THEN
3079 (K) = NUC_K(K) + 1
3080 ENDIF
3081
3082 ENDDO
3083 ENDDO
3084
3085 ENDDO
3086
3087
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
3105
3106 IF (myPE == PE_IO) THEN
3107
3108 IF(NODESI>1) THEN
3109
3110 ELSEIF(NODESJ>1) THEN
3111
3112
3113
3114 = (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
3124
3125
3126
3127
3128
3129
3130 = 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
3157
3158 = 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
3178
3179 ENDIF
3180
3181 = .TRUE.
3182
3183 IF(PRINT_STATISTICS) THEN
3184
3185
3186
3187
3188
3189
3190
3191
3192
3193
3194
3195
3196
3197 = 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
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
3236
3237
3238
3239 WRITE (*, 1000) '================================================='
3240
3241 WRITE (*, 1000) 'NOTE: ACTUAL LOAD BALANCING WILL BE COMPUTED AFTER PRE_PROCESSING.'
3242
3243
3244 ENDIF
3245
3246 ENDIF
3247
3248 CONTINUE
3249
3250
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
3273
3274
3275
3276
3277
3278
3279
3280
3281
3282
3283
3284
3285
3286
3287
3288
3289
3290 SUBROUTINE REPORT_BEST_PROCESSOR_SIZE
3291
3292
3293
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
3323
3324
3325
3326
3327
3328
3329
3330 =(numPEs==1)
3331
3332 IF(IS_SERIAL) THEN
3333
3334 = 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
3356
3357
3358
3359 = 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
3374
3375
3376
3377
3378
3379
3380
3381
3382
3383
3384
3385
3386
3387
3388
3389
3390
3391 SUBROUTINE REPORT_BEST_IJK_SIZE
3392
3393
3394
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
3424
3425
3426
3427
3428
3429
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
3439 INTEGER, ALLOCATABLE, DIMENSION(:) :: ALL_NUC_I
3440
3441 INTEGER, ALLOCATABLE, DIMENSION(:) :: ALL_LIST_I
3442 INTEGER, ALLOCATABLE, DIMENSION(:) :: GLOBAL_NUC_I
3443
3444 INTEGER, ALLOCATABLE, DIMENSION(:) :: ALL_NUC_J
3445
3446 INTEGER, ALLOCATABLE, DIMENSION(:) :: ALL_LIST_J
3447 INTEGER, ALLOCATABLE, DIMENSION(:) :: GLOBAL_NUC_J
3448
3449 INTEGER, ALLOCATABLE, DIMENSION(:) :: ALL_NUC_K
3450
3451 INTEGER, ALLOCATABLE, DIMENSION(:) :: ALL_LIST_K
3452 INTEGER, ALLOCATABLE, DIMENSION(:) :: GLOBAL_NUC_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
3477
3478 IF(NODESI*NODESJ*NODESK==1) RETURN
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
3485 (0:NODESJ-1) = jmax1-jmin1+1
3486 KSIZE_ALL(0:NODESK-1) = kmax1-kmin1+1
3487
3488
3489 IF(NODESI>1) THEN
3490
3491
3492
3493 DO I = ISTART1,IEND1
3494 NUC_I(I) = 0
3495 DO J = JSTART3,JEND3
3496 DO K = KSTART3,KEND3
3497 IF( .NOT.DEAD_CELL_AT(I,J,K)) THEN
3498 (I) = NUC_I(I) + 1
3499 ENDIF
3500 ENDDO
3501 ENDDO
3502 ENDDO
3503
3504
3505
3506
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
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
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
3540
3541 IF(NODESJ>1) THEN
3542
3543
3544
3545 DO J = JSTART1,JEND1
3546 NUC_J(J) = 0
3547 DO I = ISTART3,IEND3
3548 DO K = KSTART3,KEND3
3549 IF( .NOT.DEAD_CELL_AT(I,J,K)) THEN
3550 (J) = NUC_J(J) + 1
3551 ENDIF
3552 ENDDO
3553 ENDDO
3554 ENDDO
3555
3556
3557
3558
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
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
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
3592
3593 IF(NODESK>1) THEN
3594
3595
3596
3597 DO K = KSTART1,KEND1
3598 NUC_K(K) = 0
3599 DO I = ISTART3,IEND3
3600 DO J = JSTART3,JEND3
3601 IF( .NOT.DEAD_CELL_AT(I,J,K)) THEN
3602 (K) = NUC_K(K) + 1
3603 ENDIF
3604 ENDDO
3605 ENDDO
3606 ENDDO
3607
3608
3609
3610
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
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
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
3644
3645
3646
3647
3648
3649 IF (myPE == PE_IO) THEN
3650
3651 IF(NODESI>1) THEN
3652
3653
3654
3655 = (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
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
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
3706 = 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
3725
3726
3727 IF(NODESJ>1) THEN
3728
3729
3730 = (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
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
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
3780 = 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
3799
3800 IF(NODESK>1) THEN
3801
3802
3803 = (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
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
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
3853 = 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
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
3893 WRITE (*, 1000) '================================================='
3894
3895
3896 ENDIF
3897
3898
3899 call parallel_fin
3900
3901 STOP
3902
3903
3904
3905
3906
3907
3908
3909
3910
3911 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
3926
3927
3928
3929
3930
3931
3932
3933
3934
3935
3936
3937
3938
3939
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
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
3972
3973
3974
3975
3976
3977
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 = 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
4012
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
4034 = 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
4045
4046
4047 = MAXLOC(NCPP_WITH_GHOST,1)-1
4048 IPROC_OF_MIN = MINLOC(NCPP_WITH_GHOST,1)-1
4049
4050
4051
4052
4053
4054 RETURN
4055 END SUBROUTINE GET_LIP_WITH_GHOST_LAYERS
4056
4057
4058
4059
4060
4061
4062
4063
4064
4065
4066
4067
4068
4069
4070
4071
4072
4073
4074
4075 SUBROUTINE MINIMIZE_LOAD_IMBALANCE(NODESL,NUC_L,LMIN1,LMAX1,L_SIZE,NCPP,NCPP_WITH_GHOST)
4076
4077
4078
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
4106
4107
4108
4109
4110
4111
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
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
4129
4130
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
4138
4139
4140
4141
4142
4143
4144
4145
4146
4147
4148
4149
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
4163
4164
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
4177
4178
4179
4180
4181
4182
4183
4184
4185 IF(NOIMPROVEMENT==10) THEN
4186 WRITE (*, 1000) 'OPTIMIZED LOAD BALANCE REACHED.'
4187 EXIT
4188 ENDIF
4189
4190 ENDDO
4191
4192
4193 = BEST_L_SIZE
4194 NCPP = BEST_NCPP
4195 NCPP_WITH_GHOST = BEST_NCPP_WITH_GHOST
4196
4197
4198
4199
4200
4201
4202
4203
4204
4205
4206 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
4215
4216
4217
4218
4219
4220
4221
4222
4223
4224
4225
4226
4227
4228
4229
4230
4231
4232 SUBROUTINE REPORT_BEST_IJK_SIZE0
4233
4234
4235
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
4265
4266
4267
4268
4269
4270
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
4300
4301 IF(NODESI*NODESJ*NODESK==1) RETURN
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
4310
4311
4312
4313
4314
4315
4316 IF(NODESI>1) THEN
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
4328 (0:NODESK-1) = kmax1-kmin1+1
4329
4330
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
4339
4340 DO J = JSTART3,JEND3
4341 DO K = KSTART3,KEND3
4342
4343 IF( .NOT.DEAD_CELL_AT(I,J,K)) THEN
4344 (I) = NUC_I(I) + 1
4345 ENDIF
4346
4347 ENDDO
4348 ENDDO
4349
4350 ENDDO
4351
4352
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
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
4384 (0:NODESK-1) = kmax1-kmin1+1
4385
4386
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
4395
4396 DO I = ISTART3,IEND3
4397 DO K = KSTART3,KEND3
4398
4399 IF( .NOT.DEAD_CELL_AT(I,J,K)) THEN
4400 (J) = NUC_J(J) + 1
4401 ENDIF
4402
4403 ENDDO
4404 ENDDO
4405
4406 ENDDO
4407
4408
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
4437
4438
4439
4440
4441
4442
4443
4444 ELSEIF(NODESK>1) THEN
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
4456 (0:NODESJ-1) = jmax1-jmin1+1
4457
4458
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
4467
4468 DO I = ISTART1,IEND1
4469 DO J = JSTART1,JEND1
4470
4471 IF( .NOT.DEAD_CELL_AT(I,J,K)) THEN
4472 (K) = NUC_K(K) + 1
4473 ENDIF
4474
4475 ENDDO
4476 ENDDO
4477
4478 ENDDO
4479
4480
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
4500
4501
4502
4503
4504 IF (myPE == PE_IO) THEN
4505
4506 IF(NODESI>1) THEN
4507
4508
4509
4510 ELSEIF(NODESJ>1) THEN
4511
4512
4513
4514 = (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
4531
4532
4533
4534
4535
4536
4537
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
4559
4560
4561
4562
4563
4564
4565
4566 = 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
4605
4606
4607
4608
4609 ENDIF
4610
4611
4612 = 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
4617 = DBLE(MAXVAL_NCPP_OLD-MINVAL_NCPP_OLD)/DBLE(MINVAL_NCPP_OLD)*100.0D0
4618
4619
4620 = DBLE(MINVAL_NCPP_OLD)/DBLE(MAXVAL_NCPP_OLD)
4621
4622
4623 = 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
4630 = DBLE(MAXVAL_NCPP-MINVAL_NCPP)/DBLE(MINVAL_NCPP)*100.0D0
4631
4632
4633 = DBLE(MINVAL_NCPP)/DBLE(MAXVAL_NCPP)
4634
4635
4636 = 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
4653
4654
4655
4656 WRITE (*, 1000) '================================================='
4657
4658
4659
4660 ENDIF
4661
4662
4663
4664
4665
4666
4667
4668
4669
4670
4671
4672
4673
4674
4675 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
4688
4689
4690
4691
4692
4693
4694
4695
4696
4697
4698
4699
4700
4701
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
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
4734
4735
4736
4737
4738
4739
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 = 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
4774
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
4796 = 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
4806 = 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
4814
4815
4816 RETURN
4817 END SUBROUTINE GET_LIP_WITH_GHOST_LAYERS0
4818
4819
4820
4821
4822
4823
4824
4825
4826
4827
4828
4829
4830
4831
4832
4833
4834
4835
4836
4837 SUBROUTINE MINIMIZE_LOAD_IMBALANCE0(NODESL,NUC_L,LMIN1,LMAX1,L_SIZE,NCPP,NCPP_WITH_GHOST)
4838
4839
4840
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
4868
4869
4870
4871
4872
4873
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
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
4892
4893
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
4901
4902
4903
4904
4905
4906
4907
4908
4909
4910
4911
4912
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
4926
4927
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
4940
4941
4942
4943
4944
4945
4946
4947
4948 IF(NOIMPROVEMENT==10) THEN
4949 WRITE (*, 1000) 'OPTIMIZED LOAD BALANCE REACHED.'
4950 EXIT
4951 ENDIF
4952
4953 ENDDO
4954
4955
4956 = BEST_L_SIZE
4957 NCPP = BEST_NCPP
4958 NCPP_WITH_GHOST = BEST_NCPP_WITH_GHOST
4959
4960
4961
4962
4963
4964
4965
4966
4967
4968
4969 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
4976