File: RELATIVE:/../../../mfix.git/model/cartesian_grid/check_data_cartesian.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15 SUBROUTINE CHECK_DATA_CARTESIAN
16
17
18
19
20 USE param
21 USE param1
22 USE constant
23 USE run
24 USE physprop
25 USE indices
26 USE scalars
27 USE funits
28 USE leqsol
29 USE compar
30 USE mpi_utility
31 USE bc
32 USE DISCRETELEMENT
33
34 USE cutcell
35 USE quadric
36 USE vtk
37 USE polygon
38 USE dashboard
39 USE stl
40 IMPLICIT NONE
41
42
43
44 INTEGER :: I,J,Q,BCV
45 DOUBLE PRECISION :: norm, tan_half_angle
46 CHARACTER(LEN=9) :: GR
47
48
49
50 IF(.NOT.CARTESIAN_GRID) RETURN
51
52 IF(DISCRETE_ELEMENT) THEN
53 IF(MyPE == PE_IO) THEN
54
55 WRITE(*,10)'######################################################################'
56 WRITE(*,10)'## ##'
57 WRITE(*,10)'## ===> WARNING: RUNNING CARTESIAN GRID WITH DISCRETE ELEMENT. ##'
58 WRITE(*,10)'## THIS HAS NOT BEEN FULLY TESTED. ##'
59 WRITE(*,10)'## PLEASE USE WITH CAUTION. ##'
60 WRITE(*,10)'## ##'
61 WRITE(*,10)'######################################################################'
62
63 ENDIF
64 ENDIF
65
66
67 10 FORMAT(1X,A)
68
69
70 IF(COORDINATES=='CYLINDRICAL') THEN
71 IF(MyPE == PE_IO) THEN
72 WRITE(*,*)'INPUT ERROR: CARTESIAN GRID OPTION NOT AVAILABLE'
73 WRITE(*,*)'WITH CYLINDRICAL COORDINATE SYSTEM.'
74 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
75 ENDIF
76 CALL MFIX_EXIT(MYPE)
77 ENDIF
78
79 IF(USE_STL.AND.(.NOT.USE_MSH)) THEN
80 IF(DO_K) THEN
81 CALL GET_STL_DATA
82 ELSE
83 IF(MyPE == PE_IO) WRITE(*,*) &
84 'ERROR: STL METHOD VALID ONLY IN 3D.'
85 CALL MFIX_EXIT(MYPE)
86 ENDIF
87 IF(N_QUADRIC > 0) THEN
88 IF(MyPE == PE_IO) WRITE(*,*) 'ERROR: BOTH QUADRIC(S) AND ',&
89 'STL INPUT ARE SPECIFIED.'
90 IF(MyPE == PE_IO) WRITE(*,*) 'MFIX HANDLES ONLY ONE TYPE ',&
91 'OF SURFACE INPUT.'
92 CALL MFIX_EXIT(MYPE)
93 ENDIF
94 ENDIF
95
96 IF(USE_MSH.AND.(.NOT.USE_STL)) THEN
97 IF(DO_K) THEN
98 CALL GET_MSH_DATA
99 ELSE
100 IF(MyPE == PE_IO) WRITE(*,*) &
101 'ERROR: MSH METHOD VALID ONLY IN 3D.'
102 CALL MFIX_EXIT(MYPE)
103 ENDIF
104 IF(N_QUADRIC > 0) THEN
105 IF(MyPE == PE_IO) WRITE(*,*) 'ERROR: BOTH QUADRIC(S) AND ',&
106 'MSH INPUT ARE SPECIFIED.'
107 IF(MyPE == PE_IO) WRITE(*,*) 'MFIX HANDLES ONLY ONE TYPE ',&
108 'OF SURFACE INPUT.'
109 CALL MFIX_EXIT(MYPE)
110 ENDIF
111 ENDIF
112
113 IF(USE_POLYGON) THEN
114 IF(DO_K) THEN
115 IF(MyPE == PE_IO) WRITE(*,*) 'ERROR: POLYGON METHOD ',&
116 'VALID ONLY IN 2D.'
117 CALL MFIX_EXIT(MYPE)
118 ELSE
119 CALL GET_POLY_DATA
120 ENDIF
121 ENDIF
122
123 IF(N_QUADRIC > 0) THEN
124 IF(N_POLYGON > 0) THEN
125 IF(MyPE == PE_IO) THEN
126 WRITE(*,*) 'ERROR: BOTH QUADRIC(S) AND POLYGON(S) ',&
127 'DEFINED.'
128 WRITE(*,*) 'MFIX HANDLES ONLY ONE TYPE OF SURFACE INPUT.'
129 ENDIF
130 CALL MFIX_EXIT(MYPE)
131 ENDIF
132 IF(N_USR_DEF > 0) THEN
133 IF(MyPE == PE_IO) THEN
134 WRITE(*,*) 'ERROR: BOTH QUADRIC(S) AND USER-DEFINED ',&
135 'FUNCTION DEFINED.'
136 WRITE(*,*) 'MFIX HANDLES ONLY ONE TYPE OF SURFACE.'
137 ENDIF
138 CALL MFIX_EXIT(MYPE)
139 ENDIF
140 IF(QUADRIC_SCALE <= ZERO) THEN
141 IF(MyPE == PE_IO) THEN
142 WRITE(*,*) 'ERROR: QUADRIC_SCALE MUST BE POSITIVE.'
143 ENDIF
144 CALL MFIX_EXIT(MYPE)
145 ELSEIF(QUADRIC_SCALE /= ONE) THEN
146 DO Q = 1, N_QUADRIC
147 lambda_x(Q) = lambda_x(Q) * quadric_scale**2
148 lambda_y(Q) = lambda_y(Q) * quadric_scale**2
149 lambda_z(Q) = lambda_z(Q) * quadric_scale**2
150 Radius(Q) = Radius(Q) * quadric_scale
151 t_x(Q) = t_x(Q) * quadric_scale
152 t_y(Q) = t_y(Q) * quadric_scale
153 t_z(Q) = t_z(Q) * quadric_scale
154 clip_xmin(Q) = clip_xmin(Q) * quadric_scale
155 clip_xmax(Q) = clip_xmax(Q) * quadric_scale
156 clip_ymin(Q) = clip_ymin(Q) * quadric_scale
157 clip_ymax(Q) = clip_ymax(Q) * quadric_scale
158 clip_zmin(Q) = clip_zmin(Q) * quadric_scale
159 clip_zmax(Q) = clip_zmax(Q) * quadric_scale
160 piece_xmin(Q) = piece_xmin(Q) * quadric_scale
161 piece_xmax(Q) = piece_xmax(Q) * quadric_scale
162 piece_ymin(Q) = piece_ymin(Q) * quadric_scale
163 piece_ymax(Q) = piece_ymax(Q) * quadric_scale
164 piece_zmin(Q) = piece_zmin(Q) * quadric_scale
165 piece_zmax(Q) = piece_zmax(Q) * quadric_scale
166 ENDDO
167 ENDIF
168 ELSE
169 IF((N_POLYGON > 0).AND.(N_USR_DEF > 0)) THEN
170 IF(MyPE == PE_IO) THEN
171 WRITE(*,*) 'ERROR: POLYGON(S) AND USER-DEFINED ',&
172 'FUNCTION DEFINED.'
173 WRITE(*,*) 'MFIX HANDLES ONLY ONE TYPE OF SURFACE.'
174 ENDIF
175 CALL MFIX_EXIT(MYPE)
176 ENDIF
177 ENDIF
178
179
180 IF(N_QUADRIC > DIM_QUADRIC) THEN
181 IF(MyPE == PE_IO) THEN
182 WRITE(*,*)'INPUT ERROR: INVALID VALUE OF N_QUADRIC =', &
183 N_QUADRIC
184 WRITE(*,*)'MAXIMUM ACCEPTABLE VALUE IS DIM_QUADRIC =', &
185 DIM_QUADRIC
186 WRITE(*,*)'CHANGE MAXIMUM VALUE IN QUADRIC_MOD.F ',&
187 'IF NECESSARY.'
188 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
189 ENDIF
190 CALL MFIX_EXIT(MYPE)
191 ENDIF
192
193
194 DO Q = 1, N_QUADRIC
195
196
197
198 SELECT CASE (TRIM(QUADRIC_FORM(Q)))
199
200 CASE ('NORMAL')
201
202 lambda_x(Q) = lambda_x(Q)
203 lambda_y(Q) = lambda_y(Q)
204 lambda_z(Q) = lambda_z(Q)
205
206 norm = dsqrt(lambda_x(Q)**2 + lambda_y(Q)**2 + &
207 lambda_z(Q)**2)
208
209 IF(norm < TOL_F) THEN
210 IF(MyPE == PE_IO) THEN
211 WRITE(*,*)'INPUT ERROR: QUADRIC:', Q, &
212 ' HAS ZERO COEFFICIENTS.'
213 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
214 ENDIF
215 CALL MFIX_EXIT(MYPE)
216 ENDIF
217
218 CASE ('PLANE')
219
220 (Q) = n_x(Q)
221 lambda_y(Q) = n_y(Q)
222 lambda_z(Q) = n_z(Q)
223
224 norm = dsqrt(lambda_x(Q)**2 + lambda_y(Q)**2 + &
225 lambda_z(Q)**2)
226
227 IF( norm > TOL_F) THEN
228 lambda_x(Q) = lambda_x(Q) / norm
229 lambda_y(Q) = lambda_y(Q) / norm
230 lambda_z(Q) = lambda_z(Q) / norm
231 ELSE
232 IF(MyPE == PE_IO) THEN
233 WRITE(*,*)'INPUT ERROR: PLANE:', Q, &
234 ' HAS ZERO NORMAL VECTOR.'
235 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
236 ENDIF
237 CALL MFIX_EXIT(MYPE)
238 ENDIF
239
240 dquadric(Q) = - (lambda_x(Q)*t_x(Q) + lambda_y(Q)*t_y(Q) + &
241 lambda_z(Q)*t_z(Q))
242
243 CASE ('X_CYL_INT')
244
245
246 IF( Radius(Q) <= ZERO) THEN
247 IF(MyPE == PE_IO) THEN
248 WRITE(*,*)'INPUT ERROR: CYLINDER:', Q, &
249 ' HAS ZERO RADIUS.'
250 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
251 ENDIF
252 CALL MFIX_EXIT(MYPE)
253 ELSE
254 lambda_x(Q) = ZERO
255 lambda_y(Q) = ONE
256 lambda_z(Q) = ONE
257 dquadric(Q) = -Radius(Q)**2
258 ENDIF
259
260 CASE ('Y_CYL_INT')
261
262
263 IF( Radius(Q) <= ZERO) THEN
264 IF(MyPE == PE_IO) THEN
265 WRITE(*,*)'INPUT ERROR: CYLINDER:', Q, &
266 ' HAS ZERO RADIUS.'
267 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
268 ENDIF
269 CALL MFIX_EXIT(MYPE)
270 ELSE
271 lambda_x(Q) = ONE
272 lambda_y(Q) = ZERO
273 lambda_z(Q) = ONE
274 dquadric(Q) = -Radius(Q)**2
275 ENDIF
276
277 CASE ('Z_CYL_INT')
278
279
280 IF( Radius(Q) <= ZERO) THEN
281 IF(MyPE == PE_IO) THEN
282 WRITE(*,*)'INPUT ERROR: CYLINDER:', Q, &
283 ' HAS ZERO RADIUS.'
284 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
285 ENDIF
286 CALL MFIX_EXIT(MYPE)
287 ELSE
288 lambda_x(Q) = ONE
289 lambda_y(Q) = ONE
290 lambda_z(Q) = ZERO
291 dquadric(Q) = -Radius(Q)**2
292 ENDIF
293
294
295 CASE ('X_CYL_EXT')
296
297
298 IF( Radius(Q) <= ZERO) THEN
299 IF(MyPE == PE_IO) THEN
300 WRITE(*,*)'INPUT ERROR: CYLINDER:', Q, &
301 ' HAS ZERO RADIUS.'
302 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
303 ENDIF
304 CALL MFIX_EXIT(MYPE)
305 ELSE
306 lambda_x(Q) = ZERO
307 lambda_y(Q) = -ONE
308 lambda_z(Q) = -ONE
309 dquadric(Q) = Radius(Q)**2
310 ENDIF
311
312 CASE ('Y_CYL_EXT')
313
314
315 IF( Radius(Q) <= ZERO) THEN
316 IF(MyPE == PE_IO) THEN
317 WRITE(*,*)'INPUT ERROR: CYLINDER:', Q, &
318 ' HAS ZERO RADIUS.'
319 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
320 ENDIF
321 CALL MFIX_EXIT(MYPE)
322 ELSE
323 lambda_x(Q) = -ONE
324 lambda_y(Q) = ZERO
325 lambda_z(Q) = -ONE
326 dquadric(Q) = Radius(Q)**2
327 ENDIF
328
329 CASE ('Z_CYL_EXT')
330
331
332 IF( Radius(Q) <= ZERO) THEN
333 IF(MyPE == PE_IO) THEN
334 WRITE(*,*)'INPUT ERROR: CYLINDER:', Q, &
335 ' HAS ZERO RADIUS.'
336 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
337 ENDIF
338 CALL MFIX_EXIT(MYPE)
339 ELSE
340 lambda_x(Q) = -ONE
341 lambda_y(Q) = -ONE
342 lambda_z(Q) = ZERO
343 dquadric(Q) = Radius(Q)**2
344 ENDIF
345
346 CASE ('SPHERE_INT')
347
348
349 IF( Radius(Q) <= ZERO) THEN
350 WRITE(*,*)'INPUT ERROR: SPHERE:', Q, &
351 ' HAS INVALID RADIUS.'
352 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
353 CALL MFIX_EXIT(MYPE)
354 ELSE
355 lambda_x(Q) = ONE
356 lambda_y(Q) = ONE
357 lambda_z(Q) = ONE
358 dquadric(Q) = -Radius(Q)**2
359 ENDIF
360
361 CASE ('SPHERE_EXT')
362
363
364 IF( Radius(Q) <= ZERO) THEN
365 WRITE(*,*)'INPUT ERROR: SPHERE:', Q, &
366 ' HAS INVALID RADIUS.'
367 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
368 CALL MFIX_EXIT(MYPE)
369 ELSE
370 lambda_x(Q) = -ONE
371 lambda_y(Q) = -ONE
372 lambda_z(Q) = -ONE
373 dquadric(Q) = Radius(Q)**2
374 ENDIF
375
376
377 CASE ('X_CONE')
378
379
380 IF(HALF_ANGLE(Q) <= ZERO .OR. HALF_ANGLE(Q) >= 90.0) THEN
381 IF(MyPE == PE_IO) THEN
382 WRITE(*,*)'INPUT ERROR: CONE:', Q, &
383 ' HAS INCORRECT HALF-ANGLE.'
384 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
385 ENDIF
386 CALL MFIX_EXIT(MYPE)
387 ELSE
388 tan_half_angle = DTAN(HALF_ANGLE(Q)/180.0*PI)
389 lambda_x(Q) = -ONE
390 lambda_y(Q) = ONE/(tan_half_angle)**2
391 lambda_z(Q) = ONE/(tan_half_angle)**2
392 dquadric(Q) = ZERO
393 ENDIF
394
395 CASE ('Y_CONE')
396
397
398 IF(HALF_ANGLE(Q) <= ZERO .OR. HALF_ANGLE(Q) >= 90.0) THEN
399 IF(MyPE == PE_IO) THEN
400 WRITE(*,*)'INPUT ERROR: CONE:', Q, &
401 ' HAS INCORRECT HALF-ANGLE.'
402 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
403 ENDIF
404 CALL MFIX_EXIT(MYPE)
405 ELSE
406 tan_half_angle = DTAN(HALF_ANGLE(Q)/180.0*PI)
407 lambda_x(Q) = ONE/(tan_half_angle)**2
408 lambda_y(Q) = -ONE
409 lambda_z(Q) = ONE/(tan_half_angle)**2
410 dquadric(Q) = ZERO
411 ENDIF
412
413 CASE ('Z_CONE')
414
415
416 IF(HALF_ANGLE(Q) <= ZERO .OR. HALF_ANGLE(Q) >= 90.0) THEN
417 IF(MyPE == PE_IO) THEN
418 WRITE(*,*)'INPUT ERROR: CONE:', Q, &
419 ' HAS INCORRECT HALF-ANGLE.'
420 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
421 ENDIF
422 CALL MFIX_EXIT(MYPE)
423 ELSE
424 tan_half_angle = DTAN(HALF_ANGLE(Q)/180.0*PI)
425 lambda_x(Q) = ONE/(tan_half_angle)**2
426 lambda_y(Q) = ONE/(tan_half_angle)**2
427 lambda_z(Q) = -ONE
428 dquadric(Q) = ZERO
429 ENDIF
430
431 CASE ('C2C')
432
433
434 CALL BUILD_CONE_FOR_C2C(Q)
435
436
437 CASE ('TORUS_INT','TORUS_EXT')
438 IF((Torus_R1(Q) <= ZERO).OR.(Torus_R1(Q)==UNDEFINED)) THEN
439 IF(MyPE == PE_IO) THEN
440 WRITE(*,*)'INPUT ERROR: TORUS:', Q, &
441 ' HAS INVALID RADIUS R1:',Torus_R1(Q)
442 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
443 ENDIF
444 CALL MFIX_EXIT(MYPE)
445 ENDIF
446 IF((Torus_R2(Q) <= ZERO).OR.(Torus_R2(Q)==UNDEFINED)) THEN
447 IF(MyPE == PE_IO) THEN
448 WRITE(*,*)'INPUT ERROR: TORUS:', Q, &
449 ' HAS INVALID RADIUS R2:',Torus_R2(Q)
450 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
451 ENDIF
452 CALL MFIX_EXIT(MYPE)
453 ENDIF
454
455 CASE ('Y_UCOIL_EXT')
456 IF((UCOIL_R1(Q) <= ZERO).OR.(UCOIL_R1(Q)==UNDEFINED)) THEN
457 IF(MyPE == PE_IO) THEN
458 WRITE(*,*)'INPUT ERROR: Y_UCOIL_EXT:', Q, &
459 ' HAS INVALID RADIUS R1:',UCOIL_R1(Q)
460 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
461 ENDIF
462 CALL MFIX_EXIT(MYPE)
463 ENDIF
464 IF((UCOIL_R2(Q) <= ZERO).OR.(UCOIL_R2(Q)==UNDEFINED)) THEN
465 IF(MyPE == PE_IO) THEN
466 WRITE(*,*)'INPUT ERROR: Y_UCOIL_EXT:', Q, &
467 ' HAS INVALID RADIUS R2:',UCOIL_R2(Q)
468 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
469 ENDIF
470 CALL MFIX_EXIT(MYPE)
471 ENDIF
472
473 IF(UCOIL_Y2(Q)<UCOIL_Y1(Q)) THEN
474 IF(MyPE == PE_IO) THEN
475 WRITE(*,*)'INPUT ERROR: Y_UCOIL_EXT:', Q, &
476 ' COIL_Y2 < COIL_Y1: Y2,Y1=',UCOIL_Y2(Q),UCOIL_Y1(Q)
477 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
478 ENDIF
479 CALL MFIX_EXIT(MYPE)
480 ENDIF
481
482 CASE ('Y_UCOIL2_EXT')
483 IF((UCOIL_R1(Q) <= ZERO).OR.(UCOIL_R1(Q)==UNDEFINED)) THEN
484 IF(MyPE == PE_IO) THEN
485 WRITE(*,*)'INPUT ERROR: Y_UCOIL2_EXT:', Q, &
486 ' HAS INVALID RADIUS R1:',UCOIL_R1(Q)
487 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
488 ENDIF
489 CALL MFIX_EXIT(MYPE)
490 ENDIF
491 IF((UCOIL_R2(Q) <= ZERO).OR.(UCOIL_R2(Q)==UNDEFINED)) THEN
492 IF(MyPE == PE_IO) THEN
493 WRITE(*,*)'INPUT ERROR: Y_UCOIL2_EXT:', Q, &
494 ' HAS INVALID RADIUS R2:',UCOIL_R2(Q)
495 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
496 ENDIF
497 CALL MFIX_EXIT(MYPE)
498 ENDIF
499
500 IF(UCOIL_Y2(Q)<UCOIL_Y1(Q)) THEN
501 IF(MyPE == PE_IO) THEN
502 WRITE(*,*)'INPUT ERROR: Y_UCOIL2_EXT:', Q, &
503 ' COIL_Y2 < COIL_Y1: Y2,Y1=',UCOIL_Y2(Q),UCOIL_Y1(Q)
504 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
505 ENDIF
506 CALL MFIX_EXIT(MYPE)
507 ENDIF
508
509 CASE ('XY_BEND_INT')
510 IF((BEND_R1(Q) <= ZERO).OR.(BEND_R1(Q)==UNDEFINED)) THEN
511 IF(MyPE == PE_IO) THEN
512 WRITE(*,*)'INPUT ERROR: XY_BEND_INT:', Q, &
513 ' HAS INVALID RADIUS R1:',BEND_R1(Q)
514 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
515 ENDIF
516 CALL MFIX_EXIT(MYPE)
517 ENDIF
518
519 IF((BEND_R2(Q) <= ZERO).OR.(BEND_R2(Q)==UNDEFINED)) THEN
520 IF(MyPE == PE_IO) THEN
521 WRITE(*,*)'INPUT ERROR: XY_BEND_INT:', Q, &
522 ' HAS INVALID RADIUS R2:',BEND_R2(Q)
523 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
524 ENDIF
525 CALL MFIX_EXIT(MYPE)
526 ENDIF
527
528 IF((BEND_THETA1(Q) < ZERO).OR.(BEND_THETA1(Q)>360.0)) THEN
529 IF(MyPE == PE_IO) THEN
530 WRITE(*,*)'INPUT ERROR: XY_BEND_INT:', Q, &
531 ' HAS INVALID ANGLE THETA1:',BEND_THETA1(Q)
532 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
533 ENDIF
534 CALL MFIX_EXIT(MYPE)
535 ENDIF
536
537 IF((BEND_THETA2(Q) < ZERO).OR.(BEND_THETA2(Q)>360.0)) THEN
538 IF(MyPE == PE_IO) THEN
539 WRITE(*,*)'INPUT ERROR: XY_BEND_INT:', Q, &
540 ' HAS INVALID ANGLE THETA2:',BEND_THETA2(Q)
541 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
542 ENDIF
543 CALL MFIX_EXIT(MYPE)
544 ENDIF
545
546 CASE ('Y_C2C_INT')
547 IF((C2C_R1(Q) <= ZERO).OR.(C2C_R1(Q)==UNDEFINED)) THEN
548 IF(MyPE == PE_IO) THEN
549 WRITE(*,*)'INPUT ERROR: Y_C2C_INT:', Q, &
550 ' HAS INVALID RADIUS R1:',C2C_R1(Q)
551 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
552 ENDIF
553 CALL MFIX_EXIT(MYPE)
554 ENDIF
555
556 IF((C2C_R2(Q) <= ZERO).OR.(C2C_R2(Q)==UNDEFINED)) THEN
557 IF(MyPE == PE_IO) THEN
558 WRITE(*,*)'INPUT ERROR: C2C_XY_INT:', Q, &
559 ' HAS INVALID RADIUS R2:',C2C_R2(Q)
560 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
561 ENDIF
562 CALL MFIX_EXIT(MYPE)
563 ENDIF
564
565 IF(C2C_Y2(Q) < C2C_Y1(Q)) THEN
566 IF(MyPE == PE_IO) THEN
567 WRITE(*,*)'INPUT ERROR: Y_C2C_INT:', Q
568 WRITE(*,*)'MUST HAVE C2C_Y2 >= C2C_Y1.'
569 WRITE(*,*)'C2C_Y1,C2C_Y2 =', C2C_Y1(Q),C2C_Y2(Q)
570 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
571 ENDIF
572 CALL MFIX_EXIT(MYPE)
573 ENDIF
574
575 IF((C2C_Y1(Q) == C2C_Y2(Q)).AND.(C2C_R1(Q)/=C2C_R2(Q))) THEN
576 IF(MyPE == PE_IO) THEN
577 WRITE(*,*)'INPUT ERROR: Y_C2C_INT:', Q, &
578 ' C2C_Y1=C2C_Y2 BUT C2C_R1/=C2C_R2:'
579 WRITE(*,*)'C2C_Y1,C2C_Y2 =', C2C_Y1(Q),C2C_Y2(Q)
580 WRITE(*,*)'C2C_R1,C2C_R2 =', C2C_R1(Q),C2C_R2(Q)
581 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
582 ENDIF
583 CALL MFIX_EXIT(MYPE)
584 ENDIF
585
586 CASE ('REACTOR1')
587
588 IF(REACTOR1_Y2(Q) < REACTOR1_Y1(Q)) THEN
589 IF(MyPE == PE_IO) THEN
590 WRITE(*,*)'INPUT ERROR: REACTOR1:', Q
591 WRITE(*,*)'MUST HAVE REACTOR1_Y2 >= REACTOR1_Y1.'
592 WRITE(*,*)'REACTOR1_Y1,REACTOR1_Y2 =', REACTOR1_Y1(Q),REACTOR1_Y2(Q)
593 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
594 ENDIF
595 CALL MFIX_EXIT(MYPE)
596 ENDIF
597
598 IF((REACTOR1_Y1(Q) == REACTOR1_Y2(Q)).AND.(REACTOR1_R1(Q)/=REACTOR1_R2(Q))) THEN
599 IF(MyPE == PE_IO) THEN
600 WRITE(*,*)'INPUT ERROR: REACTOR1:', Q, &
601 ' REACTOR1_Y1=REACTOR1_Y2 BUT REACTOR1_R1/=REACTOR1_R2:'
602 WRITE(*,*)'REACTOR1_Y1,REACTOR1_Y2 =', REACTOR1_Y1(Q),REACTOR1_Y2(Q)
603 WRITE(*,*)'REACTOR1_R1,REACTOR1_R2 =', REACTOR1_R1(Q),REACTOR1_R2(Q)
604 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
605 ENDIF
606 CALL MFIX_EXIT(MYPE)
607 ENDIF
608
609
610 IF(REACTOR1_YR2(Q) <= REACTOR1_Y2(Q)) THEN
611 IF(MyPE == PE_IO) THEN
612 WRITE(*,*)'INPUT ERROR: REACTOR1:', Q
613 WRITE(*,*)'MUST HAVE REACTOR1_YR2 > REACTOR1_Y2.'
614 WRITE(*,*)'REACTOR1_YR2,REACTOR1_Y2 =', REACTOR1_YR2(Q),REACTOR1_Y2(Q)
615 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
616 ENDIF
617 CALL MFIX_EXIT(MYPE)
618 ENDIF
619
620 IF(REACTOR1_YR1(Q) >= REACTOR1_Y1(Q)) THEN
621 IF(MyPE == PE_IO) THEN
622 WRITE(*,*)'INPUT ERROR: REACTOR1:', Q
623 WRITE(*,*)'MUST HAVE REACTOR1_YR1 < REACTOR1_Y1.'
624 WRITE(*,*)'REACTOR1_YR1,REACTOR1_Y1 =', REACTOR1_YR1(Q),REACTOR1_Y1(Q)
625 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
626 ENDIF
627 CALL MFIX_EXIT(MYPE)
628 ENDIF
629
630 IF(REACTOR1_THETA1(Q) <= ZERO.OR.REACTOR1_THETA1(Q) > 90.0D0) THEN
631 IF(MyPE == PE_IO) THEN
632 WRITE(*,*)'INPUT ERROR: REACTOR1:', Q
633 WRITE(*,*)'MUST HAVE 0.0 < REACTOR1_THETA1 <= 90 DEGREES.'
634 WRITE(*,*)'REACTOR1_THETA1 =', REACTOR1_THETA1(Q)
635 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
636 ENDIF
637 CALL MFIX_EXIT(MYPE)
638 ENDIF
639
640 IF(REACTOR1_THETA2(Q) <= ZERO.OR.REACTOR1_THETA2(Q) > 90.0D0) THEN
641 IF(MyPE == PE_IO) THEN
642 WRITE(*,*)'INPUT ERROR: REACTOR1:', Q
643 WRITE(*,*)'MUST HAVE 0.0 < REACTOR1_THETA2 <= 90 DEGREES.'
644 WRITE(*,*)'REACTOR1_THETA2 =', REACTOR1_THETA2(Q)
645 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
646 ENDIF
647 CALL MFIX_EXIT(MYPE)
648 ENDIF
649
650 (Q) = REACTOR1_THETA1(Q)/180.0D0*PI
651 REACTOR1_THETA2(Q) = REACTOR1_THETA2(Q)/180.0D0*PI
652
653
654 CASE DEFAULT
655 IF(MyPE == PE_IO) THEN
656 WRITE(*,*)'INPUT ERROR: QUADRIC:', Q, &
657 ' HAS INCORRECT FORM: ',quadric_form(Q)
658 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
659 ENDIF
660 CALL MFIX_EXIT(MYPE)
661
662 END SELECT
663
664 IF(BC_ID_Q(Q) == UNDEFINED_I) THEN
665 IF(MyPE == PE_IO) THEN
666 WRITE(*,*)'INPUT ERROR: QUADRIC:', Q, &
667 ' HAS NO ASSIGNED BC ID.'
668 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
669 ENDIF
670 CALL MFIX_EXIT(MYPE)
671 ENDIF
672
673 ENDDO
674
675
676 IF(N_QUADRIC>0) THEN
677
678
679 IF(N_GROUP > DIM_GROUP) THEN
680 IF(MyPE == PE_IO) THEN
681 WRITE(*,*)'INPUT ERROR: INVALID VALUE OF N_GROUP =', N_GROUP
682 WRITE(*,*)'MAXIMUM ACCEPTABLE VALUE IS DIM_GROUP =', DIM_GROUP
683 WRITE(*,*)'CHANGE MAXIMUM VALUE IN QUADRIC_MOD.F IF NECESSARY.'
684 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
685 ENDIF
686 CALL MFIX_EXIT(MYPE)
687 ENDIF
688
689
690 DO I = 1,N_GROUP
691
692 IF(GROUP_SIZE(I) < 1 .OR. GROUP_SIZE(I) > N_QUADRIC) THEN
693 IF(MyPE == PE_IO) THEN
694 WRITE(*,*)'INPUT ERROR: GROUP:', I, ' HAS INCORRECT SIZE:', GROUP_SIZE(I)
695 WRITE(*,*)'VALID GROUP SIZE RANGE IS:', 1, ' TO ', N_QUADRIC
696 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
697 ENDIF
698 CALL MFIX_EXIT(MYPE)
699 ENDIF
700
701 DO J = 1,GROUP_SIZE(I)
702 IF(GROUP_Q(I,J) < 1 .OR. GROUP_Q(I,J) > N_QUADRIC) THEN
703 IF(MyPE == PE_IO) THEN
704 WRITE(*,*)'INPUT ERROR: GROUP_Q(', I,',',J, ') HAS INCORRECT VALUE:', GROUP_Q(I,J)
705 WRITE(*,*)'VALID GROUP_Q RANGE IS:', 1, ' TO ', N_QUADRIC
706 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
707 ENDIF
708 CALL MFIX_EXIT(MYPE)
709 ENDIF
710 ENDDO
711
712 GR = TRIM(GROUP_RELATION(I))
713
714 IF(GR/='OR'.AND.GR/='AND'.AND.GR/='PIECEWISE') THEN
715 IF(MyPE == PE_IO) THEN
716 WRITE(*,*)'INPUT ERROR: GROUP:', I, ' HAS INCORRECT GROUP RELATION: ', GR
717 WRITE(*,*)'VALID GROUP RELATIONS ARE ''OR'',''AND'', AND ''PIECEWISE''. '
718 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
719 ENDIF
720 CALL MFIX_EXIT(MYPE)
721 ENDIF
722
723 ENDDO
724
725 DO I = 2,N_GROUP
726
727 GR = TRIM(RELATION_WITH_PREVIOUS(I))
728
729 IF(GR/='OR'.AND.GR/='AND') THEN
730 IF(MyPE == PE_IO) THEN
731 WRITE(*,*)'INPUT ERROR: GROUP:', I, ' HAS INCORRECT RELATION WITH PREVIOUS: ', GR
732 WRITE(*,*)'VALID GROUP RELATIONS ARE ''OR'', AND ''AND''. '
733 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
734 ENDIF
735 CALL MFIX_EXIT(MYPE)
736 ENDIF
737
738 ENDDO
739
740 ENDIF
741
742
743 IF(TOL_SNAP(1)<ZERO.OR.TOL_SNAP(1)>HALF) THEN
744 IF(MyPE == PE_IO) THEN
745 WRITE(*,*)'INPUT ERROR: INVALID VALUE OF TOL_SNAP IN X-DIRECTION =', TOL_SNAP(1)
746 WRITE(*,*)'ACCEPTABLE VALUES ARE BETWEEN 0.0 AND 0.5.'
747 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
748 ENDIF
749 CALL MFIX_EXIT(MYPE)
750 ENDIF
751
752 IF(TOL_SNAP(2)==UNDEFINED) TOL_SNAP(2)=TOL_SNAP(1)
753
754 IF(TOL_SNAP(2)<ZERO.OR.TOL_SNAP(2)>HALF) THEN
755 IF(MyPE == PE_IO) THEN
756 WRITE(*,*)'INPUT ERROR: INVALID VALUE OF TOL_SNAP IN Y-DIRECTION =', TOL_SNAP(2)
757 WRITE(*,*)'ACCEPTABLE VALUES ARE BETWEEN 0.0 AND 0.5.'
758 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
759 ENDIF
760 CALL MFIX_EXIT(MYPE)
761 ENDIF
762
763 IF(TOL_SNAP(3)==UNDEFINED) TOL_SNAP(3)=TOL_SNAP(1)
764
765 IF(TOL_SNAP(3)<ZERO.OR.TOL_SNAP(3)>HALF) THEN
766 IF(MyPE == PE_IO) THEN
767 WRITE(*,*)'INPUT ERROR: INVALID VALUE OF TOL_SNAP IN Z-DIRECTION =', TOL_SNAP(3)
768 WRITE(*,*)'ACCEPTABLE VALUES ARE BETWEEN 0.0 AND 0.5.'
769 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
770 ENDIF
771 CALL MFIX_EXIT(MYPE)
772 ENDIF
773
774
775 IF(TOL_DELH<ZERO.OR.TOL_DELH>ONE) THEN
776 IF(MyPE == PE_IO) THEN
777 WRITE(*,*)'INPUT ERROR: INVALID VALUE OF TOL_DELH =', TOL_DELH
778 WRITE(*,*)'ACCEPTABLE VALUES ARE BETWEEN 0.0 AND 1.0.'
779 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
780 ENDIF
781 CALL MFIX_EXIT(MYPE)
782 ENDIF
783
784 IF(TOL_SMALL_CELL<ZERO.OR.TOL_SMALL_CELL>ONE) THEN
785 IF(MyPE == PE_IO) THEN
786 WRITE(*,*)'INPUT ERROR: INVALID VALUE OF TOL_SMALL_CELL =', TOL_SMALL_CELL
787 WRITE(*,*)'ACCEPTABLE VALUES ARE BETWEEN 0.0 AND 1.0.'
788 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
789 ENDIF
790 CALL MFIX_EXIT(MYPE)
791 ENDIF
792
793 IF(TOL_SMALL_AREA<ZERO.OR.TOL_SMALL_AREA>ONE) THEN
794 IF(MyPE == PE_IO) THEN
795 WRITE(*,*)'INPUT ERROR: INVALID VALUE OF TOL_SMALL_AREA =', TOL_SMALL_AREA
796 WRITE(*,*)'ACCEPTABLE VALUES ARE BETWEEN 0.0 AND 1.0.'
797 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
798 ENDIF
799 CALL MFIX_EXIT(MYPE)
800 ENDIF
801
802 IF(ALPHA_MAX<ZERO) THEN
803 IF(MyPE == PE_IO) THEN
804 WRITE(*,*)'INPUT ERROR: NEGATIVE VALUE OF ALPHA_MAX =', ALPHA_MAX
805 WRITE(*,*)'ACCEPTABLE VALUES ARE POSITIVE NUMBERS (E.G. 1.0).'
806 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
807 ENDIF
808 CALL MFIX_EXIT(MYPE)
809 ENDIF
810
811
812 IF(TOL_F<ZERO) THEN
813 IF(MyPE == PE_IO) THEN
814 WRITE(*,*)'INPUT ERROR: NEGATIVE VALUE OF TOL_F =', TOL_F
815 WRITE(*,*)'ACCEPTABLE VALUES ARE SMALL POSITIVE NUMBERS (E.G. 1.0E-9).'
816 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
817 ENDIF
818 CALL MFIX_EXIT(MYPE)
819 ENDIF
820
821 IF(TOL_POLY<ZERO) THEN
822 IF(MyPE == PE_IO) THEN
823 WRITE(*,*)'INPUT ERROR: NEGATIVE VALUE OF TOL_POLY =', TOL_POLY
824 WRITE(*,*)'ACCEPTABLE VALUES ARE SMALL POSITIVE NUMBERS (E.G. 1.0E-9).'
825 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
826 ENDIF
827 CALL MFIX_EXIT(MYPE)
828 ENDIF
829
830 IF(ITERMAX_INT<0) THEN
831 IF(MyPE == PE_IO) THEN
832 WRITE(*,*)'INPUT ERROR: NEGATIVE VALUE OF ITERMAX_INT =', ITERMAX_INT
833 WRITE(*,*)'ACCEPTABLE VALUES ARE LARGE POSITIVE INTEGERS (E.G. 10000).'
834 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
835 ENDIF
836 CALL MFIX_EXIT(MYPE)
837 ENDIF
838
839 IF(FAC_DIM_MAX_CUT_CELL<0.05.OR.FAC_DIM_MAX_CUT_CELL>ONE) THEN
840 IF(MyPE == PE_IO) THEN
841 WRITE(*,*)'INPUT ERROR: NEGATIVE VALUE OF FAC_DIM_MAX_CUT_CELL =', FAC_DIM_MAX_CUT_CELL
842 WRITE(*,*)'ACCEPTABLE VALUES ARE BETWEEN 0.05 AND 1.0.'
843 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
844 ENDIF
845 CALL MFIX_EXIT(MYPE)
846 ENDIF
847
848
849 IF((CG_SAFE_MODE(1)==1).AND.(PG_OPTION/=0)) THEN
850 PG_OPTION = 0
851 IF(MyPE == PE_IO) WRITE(*,*)'WARNING: SAFE_MODE ACTIVATED FOR GAS PRESSURE, REVERTING TO PG_OPTION = 0'
852 ENDIF
853
854 IF(PG_OPTION <0 .OR. PG_OPTION>2) THEN
855 IF(MyPE == PE_IO) THEN
856 WRITE(*,*)'INPUT ERROR: INVALID VALUE OF PG_OPTION =', PG_OPTION
857 WRITE(*,*)'ACCEPTABLE VALUES ARE 0,1,AND 2.'
858 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
859 ENDIF
860 CALL MFIX_EXIT(MYPE)
861 ENDIF
862
863 IF(CG_UR_FAC(2)<ZERO.OR.CG_UR_FAC(2)>ONE) THEN
864 IF(MyPE == PE_IO) THEN
865 WRITE(*,*)'INPUT ERROR: INVALID VALUE OF CG_UR_FAC(2) =', CG_UR_FAC(2)
866 WRITE(*,*)'ACCEPTABLE VALUES ARE BETWEEN 0.0 AND 1.0.'
867 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
868 ENDIF
869 CALL MFIX_EXIT(MYPE)
870 ENDIF
871
872 IF(BAR_WIDTH<10.OR.BAR_WIDTH>80) THEN
873 IF(MyPE == PE_IO) THEN
874 WRITE(*,*)'INPUT ERROR: INVALID VALUE OF BAR_WIDTH =', BAR_WIDTH
875 WRITE(*,*)'ACCEPTABLE VALUES ARE BETWEEN 10 AND 80.'
876 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
877 ENDIF
878 CALL MFIX_EXIT(MYPE)
879 ENDIF
880
881 IF(BAR_RESOLUTION<ONE.OR.BAR_RESOLUTION>100.0) THEN
882 IF(MyPE == PE_IO) THEN
883 WRITE(*,*)'INPUT ERROR: INVALID VALUE OF BAR_RESOLUTION =', BAR_RESOLUTION
884 WRITE(*,*)'ACCEPTABLE VALUES ARE BETWEEN 0.0 AND 100.0.'
885 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
886 ENDIF
887 CALL MFIX_EXIT(MYPE)
888 ENDIF
889
890 IF(F_DASHBOARD<1) THEN
891 IF(MyPE == PE_IO) THEN
892 WRITE(*,*)'INPUT ERROR: INVALID VALUE OF F_DASHBOARD =', F_DASHBOARD
893 WRITE(*,*)'ACCEPTABLE VALUES ARE INTEGERS >= 1.'
894 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
895 ENDIF
896 CALL MFIX_EXIT(MYPE)
897 ENDIF
898
899
900
901
902
903
904 = TIME
905 SMMIN = LARGE_NUMBER
906 SMMAX = -LARGE_NUMBER
907
908 DTMIN = LARGE_NUMBER
909 DTMAX = -LARGE_NUMBER
910
911 NIT_MIN = MAX_NIT
912 NIT_MAX = 0
913
914 N_DASHBOARD = 0
915
916
917
918 CG_MI_CONVERTED_TO_PS = .FALSE.
919
920 DO BCV = 1, DIMENSION_BC
921
922 IF (BC_TYPE(BCV) == 'CG_MI') THEN
923 BC_TYPE(BCV) = 'CG_NSW'
924 CG_MI_CONVERTED_TO_PS(BCV) = .TRUE.
925 if(MyPE==0) print*,'From check_data_cartesian: Converted CG_MI to CG_NSW for BC#',BCV
926 ENDIF
927
928 ENDDO
929
930
931 IF(RE_INDEXING) THEN
932
933 IF(MyPE==0) THEN
934 WRITE(*,*)' From check_data_cartesian: RE_INDEXING is turned on.'
935 WRITE(*,*)' The preconditionner will be turned off for all equations'
936 WRITE(*,*)' regardless of the mfix.dat setting.'
937 LEQ_PC = 'NONE'
938 ENDIF
939 ENDIF
940
941 RETURN
942 END SUBROUTINE CHECK_DATA_CARTESIAN
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958 SUBROUTINE CHECK_BC_FLAGS
959
960
961
962
963 USE bc
964 USE compar
965 USE constant
966 USE cutcell
967 USE dashboard
968 USE fldvar
969 USE functions
970 USE funits
971 USE indices
972 USE leqsol
973 USE mpi_utility
974 USE param
975 USE param1
976 USE physprop
977 USE polygon
978 USE quadric
979 USE run
980 USE scalars
981 USE toleranc
982 USE vtk
983 use error_manager
984
985 IMPLICIT NONE
986
987
988
989 INTEGER :: IJK,IJKW,IJKS,IJKB,M,N
990 INTEGER :: IJKWW,IJKSS,IJKBB
991 INTEGER :: BCV,BCV_U,BCV_V,BCV_W
992
993 DOUBLE PRECISION SUM, SUM_EP
994
995
996
997
998
999 CALL INIT_ERR_MSG("CHECK_BC_FLAGS")
1000
1001 DO BCV = 1, DIMENSION_BC
1002 IF(CG_MI_CONVERTED_TO_PS(BCV)) THEN
1003 IF(BC_MASSFLOW_g(BCV)==UNDEFINED) THEN
1004 WRITE(ERR_MSG, 1710) BCV
1005 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
1006 ENDIF
1007 DO M = 1, MMAX
1008 IF(BC_MASSFLOW_s(BCV,M)==UNDEFINED) THEN
1009 WRITE(ERR_MSG, 1711) BCV, M
1010 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
1011 ENDIF
1012 ENDDO
1013 ENDIF
1014 ENDDO
1015
1016 1710 FORMAT('Error 1110: BC :',I3,'. When using CG_MI, the gas mass flow rate',/1X, &
1017 'must be specified, including when it is zero.',/1X, &
1018 ' Please correct the mfix.dat file.')
1019
1020 1711 FORMAT('Error 1111: BC :',I3,'. When using CG_MI, the solids mass flow rate',/1X, &
1021 'for M=',I4,' must be specified, including when it is zero.',/1X, &
1022 ' Please correct the mfix.dat file.')
1023
1024 DO IJK = ijkstart3, ijkend3
1025 BCV = BC_ID(IJK)
1026 IF(BCV>0) THEN
1027
1028 IF(BC_TYPE(BCV) == 'CG_MI') THEN
1029
1030 print*,'CG_MI at', IJK
1031
1032
1033
1034
1035
1036
1037 ELSEIF(BC_TYPE(BCV) == 'CG_PO') THEN
1038
1039 FLAG(IJK) = 11
1040 FLAG_E(IJK) = UNDEFINED_I
1041 FLAG_N(IJK) = UNDEFINED_I
1042 FLAG_T(IJK) = UNDEFINED_I
1043
1044 IJKW = WEST_OF(IJK)
1045 BCV_U = BC_U_ID(IJKW)
1046 IF(BCV_U>0) THEN
1047 IF(BC_TYPE(BCV_U) == 'CG_PO') THEN
1048 FLAG(IJKW) = 11
1049 FLAG_E(IJKW) = UNDEFINED_I
1050 FLAG_N(IJKW) = UNDEFINED_I
1051 FLAG_T(IJKW) = UNDEFINED_I
1052 ENDIF
1053 ENDIF
1054
1055 IJKS = SOUTH_OF(IJK)
1056 BCV_V = BC_V_ID(IJKS)
1057 IF(BCV_V>0) THEN
1058 IF(BC_TYPE(BCV_V) == 'CG_PO') THEN
1059 FLAG(IJKS) = 11
1060 FLAG_E(IJKS) = UNDEFINED_I
1061 FLAG_N(IJKS) = UNDEFINED_I
1062 FLAG_T(IJKS) = UNDEFINED_I
1063 ENDIF
1064 ENDIF
1065
1066 IF(DO_K) THEN
1067 IJKB = BOTTOM_OF(IJK)
1068 BCV_W = BC_W_ID(IJKB)
1069 IF(BCV_W>0) THEN
1070 IF(BC_TYPE(BCV_W) == 'CG_PO') THEN
1071 FLAG(IJKB) = 11
1072 FLAG_E(IJKB) = UNDEFINED_I
1073 FLAG_N(IJKB) = UNDEFINED_I
1074 FLAG_T(IJKB) = UNDEFINED_I
1075 ENDIF
1076 ENDIF
1077 ENDIF
1078
1079 ENDIF
1080 ENDIF
1081 ENDDO
1082
1083
1084 DO IJK = ijkstart3, ijkend3
1085 BCV = BC_ID(IJK)
1086 IF(BCV>0) THEN
1087 IF(BC_TYPE(BCV) == 'CG_MI') THEN
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104 IF (BC_U_G(BCV) == UNDEFINED) THEN
1105 IF (NO_I) THEN
1106 BC_U_G(BCV) = ZERO
1107 ELSEIF(BC_VOLFLOW_g(BCV)==UNDEFINED.AND. &
1108 BC_MASSFLOW_g(BCV)==UNDEFINED.AND.&
1109 BC_VELMAG_g(BCV)==UNDEFINED) THEN
1110 IF(DMP_LOG)WRITE (UNIT_LOG, 900) 'BC_U_g', BCV
1111 call mfix_exit(myPE)
1112 ENDIF
1113 ENDIF
1114 IF (BC_V_G(BCV) == UNDEFINED) THEN
1115 IF (NO_J) THEN
1116 BC_V_G(BCV) = ZERO
1117 ELSEIF(BC_VOLFLOW_g(BCV)==UNDEFINED.AND. &
1118 BC_MASSFLOW_g(BCV)==UNDEFINED.AND.&
1119 BC_VELMAG_g(BCV)==UNDEFINED) THEN
1120 IF(DMP_LOG)WRITE (UNIT_LOG, 900) 'BC_V_g', BCV
1121 call mfix_exit(myPE)
1122 ENDIF
1123 ENDIF
1124 IF (BC_W_G(BCV) == UNDEFINED) THEN
1125 IF (NO_K) THEN
1126 BC_W_G(BCV) = ZERO
1127 ELSEIF(BC_VOLFLOW_g(BCV)==UNDEFINED.AND. &
1128 BC_MASSFLOW_g(BCV)==UNDEFINED.AND.&
1129 BC_VELMAG_g(BCV)==UNDEFINED) THEN
1130 IF(DMP_LOG)WRITE (UNIT_LOG, 900) 'BC_W_g', BCV
1131 call mfix_exit(myPE)
1132 ENDIF
1133 ENDIF
1134 IF (K_Epsilon .AND. BC_K_Turb_G(BCV) == UNDEFINED) THEN
1135 IF(DMP_LOG)WRITE (UNIT_LOG, 1000) 'BC_K_Turb_G', BCV
1136 call mfix_exit(myPE)
1137 ENDIF
1138 IF (K_Epsilon .AND. BC_E_Turb_G(BCV) == UNDEFINED) THEN
1139 IF(DMP_LOG)WRITE (UNIT_LOG, 1000) 'BC_E_Turb_G', BCV
1140 call mfix_exit(myPE)
1141 ENDIF
1142
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 = BC_EP_G(BCV)
1178 DO M = 1, MMAX
1179 IF (BC_ROP_S(BCV,M) == UNDEFINED) THEN
1180 IF (BC_EP_G(BCV) == ONE) THEN
1181 BC_ROP_S(BCV,M) = ZERO
1182 ELSEIF (MMAX == 1) THEN
1183 BC_ROP_S(BCV,M) = (ONE - BC_EP_G(BCV))*RO_S0(M)
1184 ELSE
1185 IF(DMP_LOG)WRITE (UNIT_LOG, 1100) 'BC_ROP_s', BCV, M
1186 call mfix_exit(myPE)
1187 ENDIF
1188 ENDIF
1189
1190 SUM_EP = SUM_EP + BC_ROP_S(BCV,M)/RO_S0(M)
1191 IF (SPECIES_EQ(M)) THEN
1192 SUM = ZERO
1193 DO N = 1, NMAX(M)
1194 IF(BC_X_S(BCV,M,N)/=UNDEFINED)SUM=SUM+BC_X_S(BCV,M,N)
1195 ENDDO
1196
1197 IF (BC_ROP_S(BCV,M)==ZERO .AND. SUM==ZERO) THEN
1198 BC_X_S(BCV,M,1) = ONE
1199 SUM = ONE
1200 ENDIF
1201
1202 DO N = 1, NMAX(M)
1203 IF (BC_X_S(BCV,M,N) == UNDEFINED) THEN
1204 IF(.NOT.COMPARE(ONE,SUM) .AND. DMP_LOG)WRITE (UNIT_LOG,1110)BCV,M,N
1205 BC_X_S(BCV,M,N) = ZERO
1206 ENDIF
1207 ENDDO
1208
1209 IF (.NOT.COMPARE(ONE,SUM)) THEN
1210 IF(DMP_LOG)WRITE (UNIT_LOG, 1120) BCV, M
1211 call mfix_exit(myPE)
1212 ENDIF
1213 ENDIF
1214
1215 IF (BC_U_S(BCV,M) == UNDEFINED) THEN
1216 IF (BC_ROP_S(BCV,M)==ZERO .OR. NO_I) THEN
1217 BC_U_S(BCV,M) = ZERO
1218 ELSEIF(BC_VOLFLOW_s(BCV,M)==UNDEFINED.AND. &
1219 BC_MASSFLOW_s(BCV,M)==UNDEFINED.AND.&
1220 BC_VELMAG_s(BCV,M)==UNDEFINED) THEN
1221 IF(DMP_LOG)WRITE (UNIT_LOG, 910) 'BC_U_s', BCV, M
1222 call mfix_exit(myPE)
1223 ENDIF
1224 ENDIF
1225
1226 IF (BC_V_S(BCV,M) == UNDEFINED) THEN
1227 IF (BC_ROP_S(BCV,M)==ZERO .OR. NO_J) THEN
1228 BC_V_S(BCV,M) = ZERO
1229 ELSEIF(BC_VOLFLOW_s(BCV,M)==UNDEFINED.AND. &
1230 BC_MASSFLOW_s(BCV,M)==UNDEFINED.AND.&
1231 BC_VELMAG_s(BCV,M)==UNDEFINED) THEN
1232 IF(DMP_LOG)WRITE (UNIT_LOG, 910) 'BC_V_s', BCV, M
1233 call mfix_exit(myPE)
1234 ENDIF
1235 ENDIF
1236
1237 IF (BC_W_S(BCV,M) == UNDEFINED) THEN
1238 IF (BC_ROP_S(BCV,M)==ZERO .OR. NO_K) THEN
1239 BC_W_S(BCV,M) = ZERO
1240 ELSEIF(BC_VOLFLOW_s(BCV,M)==UNDEFINED.AND. &
1241 BC_MASSFLOW_s(BCV,M)==UNDEFINED.AND.&
1242 BC_VELMAG_s(BCV,M)==UNDEFINED) THEN
1243 IF(DMP_LOG)WRITE (UNIT_LOG, 910) 'BC_W_s', BCV, M
1244 call mfix_exit(myPE)
1245 ENDIF
1246 ENDIF
1247
1248 IF (ENERGY_EQ .AND. BC_T_S(BCV,M)==UNDEFINED) THEN
1249 IF (BC_ROP_S(BCV,M) == ZERO) THEN
1250 BC_T_S(BCV,M) = BC_T_G(BCV)
1251 ELSE
1252 IF(DMP_LOG)WRITE (UNIT_LOG, 1100) 'BC_T_s', BCV, M
1253 call mfix_exit(myPE)
1254 ENDIF
1255 ENDIF
1256
1257 IF (GRANULAR_ENERGY .AND. BC_THETA_M(BCV,M)==UNDEFINED) THEN
1258 IF (BC_ROP_S(BCV,M) == ZERO) THEN
1259 BC_THETA_M(BCV,M) = ZERO
1260 ELSE
1261 IF(DMP_LOG)WRITE (UNIT_LOG, 1100) 'BC_Theta_m', BCV, M
1262 call mfix_exit(myPE)
1263 ENDIF
1264 ENDIF
1265
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 ENDDO
1302
1303 IF (.NOT.COMPARE(ONE,SUM_EP)) THEN
1304 IF(DMP_LOG)WRITE (UNIT_LOG, 1125) BCV
1305 call mfix_exit(myPE)
1306 ENDIF
1307
1308 DO N = 1, NScalar
1309 IF (BC_Scalar(BCV,N) == UNDEFINED) THEN
1310 IF(DMP_LOG)WRITE (UNIT_LOG, 1004) 'BC_Scalar', BCV, N
1311 CALL MFIX_EXIT(myPE)
1312 ENDIF
1313 ENDDO
1314
1315
1316 ELSEIF(BC_TYPE(BCV) == 'CG_PO') THEN
1317
1318 IJKW = WEST_OF(IJK)
1319 IF(FLUID_AT(IJKW)) THEN
1320 FLAG_E(IJKW) = 2011
1321 ENDIF
1322
1323 BCV_U = BC_U_ID(IJKW)
1324 IF(BCV_U>0) THEN
1325 IF(BC_TYPE(BCV_U) == 'CG_PO') THEN
1326 IJKWW = WEST_OF(IJKW)
1327 IF(FLUID_AT(IJKWW)) THEN
1328 FLAG_E(IJKWW) = 2011
1329 ENDIF
1330 ENDIF
1331 ENDIF
1332
1333 IJKS = SOUTH_OF(IJK)
1334 IF(FLUID_AT(IJKS)) THEN
1335 FLAG_N(IJKS) = 2011
1336 ENDIF
1337
1338 BCV_V = BC_V_ID(IJKS)
1339 IF(BCV_V>0) THEN
1340 IF(BC_TYPE(BCV_V) == 'CG_PO') THEN
1341 IJKSS = SOUTH_OF(IJKS)
1342 IF(FLUID_AT(IJKSS)) THEN
1343 FLAG_N(IJKSS) = 2011
1344 ENDIF
1345 ENDIF
1346 ENDIF
1347
1348
1349 IF(DO_K) THEN
1350 IJKB = BOTTOM_OF(IJK)
1351 IF(FLUID_AT(IJKB)) THEN
1352 FLAG_T(IJKB) = 2011
1353 ENDIF
1354
1355 BCV_W = BC_W_ID(IJKB)
1356 IF(BCV_W>0) THEN
1357 IF(BC_TYPE(BCV_W) == 'CG_PO') THEN
1358 IJKBB = BOTTOM_OF(IJKB)
1359 IF(FLUID_AT(IJKBB)) THEN
1360 FLAG_T(IJKBB) = 2011
1361 ENDIF
1362 ENDIF
1363 ENDIF
1364
1365 ENDIF
1366
1367 IF (BC_P_G(BCV) == UNDEFINED) THEN
1368 IF(DMP_LOG)WRITE (UNIT_LOG, 1000) 'BC_P_g', BCV
1369 call mfix_exit(myPE)
1370 ELSEIF (BC_P_G(BCV)<=ZERO .AND. RO_G0==UNDEFINED) THEN
1371 IF(DMP_LOG)WRITE (UNIT_LOG, 1010) BCV, BC_P_G(BCV)
1372 call mfix_exit(myPE)
1373 ENDIF
1374
1375 ENDIF
1376
1377 ENDIF
1378
1379 ENDDO
1380
1381 RETURN
1382
1383
1384 900 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',A,'(',I2,&
1385 ') not specified',/1X,'One of the following must be specified:',/1X,&
1386 'BC_VOLFLOW_g, BC_MASSFLOW_g or BC_VELMAG_g',/1X,70('*')/)
1387
1388 910 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',A,'(',I2,',',I1,&
1389 ') not specified',/1X,'One of the following must be specified:',/1X,&
1390 'BC_VOLFLOW_g, BC_MASSFLOW_g or BC_VELMAG_g',/1X,70('*')/)
1391
1392 1000 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',A,'(',I2,&
1393 ') not specified',/1X,70('*')/)
1394 1001 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/&
1395 ' Message: Illegal BC_TYPE for BC # = ',I2,/' BC_TYPE = ',A,/&
1396 ' Valid BC_TYPE are: ')
1397 1002 FORMAT(5X,A16)
1398 1003 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',A,'(',I2,&
1399 ') value is unphysical',/1X,70('*')/)
1400 1004 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',A,'(',I2,',',I2,&
1401 ') not specified',/1X,70('*')/)
1402 1005 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',A,'(',I2,',',I2,&
1403 ') value is unphysical',/1X,70('*')/)
1404 1010 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/' Message: BC_P_g( ',I2,&
1405 ') = ',G12.5,/&
1406 ' Pressure should be greater than zero for compressible flow',/1X,70(&
1407 '*')/)
1408 1050 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/' Message: BC number:',I2,&
1409 ' - ',A,' should be ',A,' zero.',/1X,70('*')/)
1410 1060 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/' Message: BC_X_g(',I2,',',I2&
1411 ,') not specified',/1X,70('*')/)
1412 1065 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/' Message: BC number:',I2,&
1413 ' - Sum of gas mass fractions is NOT equal to one',/1X,70('*')/)
1414 1100 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',A,'(',I2,',',I1,&
1415 ') not specified',/1X,70('*')/)
1416 1103 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',A,'(',I2,',',I1,&
1417 ') value is unphysical',/1X,70('*')/)
1418 1104 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',A,'(',I2,',',I2,&
1419 ',',I2,') not specified',/1X,70('*')/)
1420 1105 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',A,'(',I2,',',I2,&
1421 ',',I2,') value is unphysical',/1X,70('*')/)
1422 1110 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/' Message: BC_X_s(',I2,',',I2&
1423 ,',',I2,') not specified',/1X,70('*')/)
1424 1120 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/' Message: BC number:',I2,&
1425 ' - Sum of solids-',I1,' mass fractions is NOT equal to one',/1X,70(&
1426 '*')/)
1427 1125 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/' Message: BC number:',I2,&
1428 ' - Sum of volume fractions is NOT equal to one',/1X,70('*')/)
1429 1150 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/' Message: BC number:',I2,&
1430 ' - ',A,I1,' should be ',A,' zero.',/1X,70('*')/)
1431 1160 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/&
1432 ' Message: Boundary condition no', &
1433 I2,' is a second outflow condition.',/1X,&
1434 ' Only one outflow is allowed. Consider using P_OUTFLOW.',/1X, 70('*')/)
1435 1200 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',A,'(',I2,&
1436 ') specified',' for an undefined BC location',/1X,70('*')/)
1437 1300 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/' Message: ',A,'(',I2,',',I1,&
1438 ') specified',' for an undefined BC location',/1X,70('*')/)
1439 1400 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/&
1440 ' Message: No initial or boundary condition specified',/&
1441 ' I J K')
1442 1410 FORMAT(I5,3X,I5,3X,I5)
1443 1420 FORMAT(/1X,70('*')/)
1444
1445 1500 FORMAT(/1X,70('*')//' From: CHECK_BC_FLAGS',/&
1446 ' Message: No initial or boundary condition specified',/&
1447 ' I J K')
1448
1449
1450 END SUBROUTINE CHECK_BC_FLAGS
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474 SUBROUTINE BUILD_CONE_FOR_C2C(Q)
1475
1476
1477
1478
1479 USE param
1480 USE param1
1481 USE constant
1482 USE run
1483 USE physprop
1484 USE indices
1485 USE scalars
1486 USE funits
1487 USE leqsol
1488 USE compar
1489 USE mpi_utility
1490 USE bc
1491 USE DISCRETELEMENT
1492
1493 USE cutcell
1494 USE quadric
1495 USE vtk
1496 USE polygon
1497 USE dashboard
1498 USE stl
1499
1500
1501 IMPLICIT NONE
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511 INTEGER :: Q,QM1,QP1
1512 DOUBLE PRECISION :: x1,x2,y1,y2,z1,z2,R1,R2
1513 DOUBLE PRECISION :: tan_half_angle
1514 LOGICAL :: aligned
1515
1516
1517
1518 = Q-1
1519 QP1 = Q+1
1520
1521 IF(MyPE == PE_IO) THEN
1522 WRITE(*,*)' INFO FOR QUADRIC', Q
1523 WRITE(*,*)' Defining Cone for Cylinder to Cylinder junction'
1524 WRITE(*,*)' Between Quadrics ',QM1,' AND ', QP1
1525 ENDIF
1526
1527
1528 IF((TRIM(QUADRIC_FORM(QM1))=='X_CYL_INT').AND. &
1529 (TRIM(QUADRIC_FORM(QP1))=='X_CYL_INT')) THEN
1530
1531 (Q) = 'X_CONE'
1532
1533 aligned = (t_y(QM1)==t_y(QP1)).AND.(t_z(QM1)==t_z(QP1))
1534 IF(.NOT.aligned) THEN
1535 IF(MyPE == PE_IO) THEN
1536 WRITE(*,*)' ERROR: CYLINDERS ',QM1, ' AND ', QP1, ' ARE NOT ALIGNED'
1537 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1538 ENDIF
1539 call mfix_exit(myPE)
1540 ENDIF
1541
1542 R1 = RADIUS(QM1)
1543 R2 = RADIUS(QP1)
1544 IF(R1==R2) THEN
1545 IF(MyPE == PE_IO) THEN
1546 WRITE(*,*)' ERROR: CYLINDERS ',QM1, ' AND ', QP1, ' HAVE THE SAME RADIUS'
1547 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1548 ENDIF
1549 call mfix_exit(myPE)
1550 ENDIF
1551
1552 x1 = piece_xmax(QM1)
1553 x2 = piece_xmin(QP1)
1554 IF(x2<=x1) THEN
1555 IF(MyPE == PE_IO) THEN
1556 WRITE(*,*)' ERROR: CYLINDERS ',QM1, ' AND ', QP1, ' ARE NOT PIECED PROPERLY'
1557 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1558 ENDIF
1559 call mfix_exit(myPE)
1560 ENDIF
1561
1562 tan_half_angle = (R2-R1)/(x2-x1)
1563
1564 HALF_ANGLE(Q) = DATAN(tan_half_angle)/PI*180.0D0
1565 lambda_x(Q) = -ONE
1566 lambda_y(Q) = ONE/(tan_half_angle)**2
1567 lambda_z(Q) = ONE/(tan_half_angle)**2
1568 dquadric(Q) = ZERO
1569
1570 piece_xmin(Q) = x1
1571 piece_xmax(Q) = x2
1572
1573 t_x(Q) = x1 - R1/tan_half_angle
1574 t_y(Q) = t_y(QM1)
1575 t_z(Q) = t_z(QM1)
1576
1577 IF(MyPE == PE_IO) THEN
1578 WRITE(*,*) ' QUADRIC:',Q, ' WAS DEFINED AS ', TRIM(QUADRIC_FORM(Q))
1579 WRITE(*,*) ' WITH AN HALF-ANGLE OF ', HALF_ANGLE(Q), 'DEG.'
1580 ENDIF
1581
1582
1583 ELSEIF((TRIM(QUADRIC_FORM(QM1))=='X_CYL_EXT').AND. &
1584 (TRIM(QUADRIC_FORM(QP1))=='X_CYL_EXT')) THEN
1585
1586 (Q) = 'X_CONE'
1587
1588 aligned = (t_y(QM1)==t_y(QP1)).AND.(t_z(QM1)==t_z(QP1))
1589 IF(.NOT.aligned) THEN
1590 IF(MyPE == PE_IO) THEN
1591 WRITE(*,*)' ERROR: CYLINDERS ',QM1, ' AND ', QP1, ' ARE NOT ALIGNED'
1592 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1593 ENDIF
1594 call mfix_exit(myPE)
1595 ENDIF
1596
1597 R1 = RADIUS(QM1)
1598 R2 = RADIUS(QP1)
1599 IF(R1==R2) THEN
1600 IF(MyPE == PE_IO) THEN
1601 WRITE(*,*)' ERROR: CYLINDERS ',QM1, ' AND ', QP1, ' HAVE THE SAME RADIUS'
1602 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1603 ENDIF
1604 call mfix_exit(myPE)
1605 ENDIF
1606
1607 x1 = piece_xmax(QM1)
1608 x2 = piece_xmin(QP1)
1609 IF(x2<=x1) THEN
1610 IF(MyPE == PE_IO) THEN
1611 WRITE(*,*)' ERROR: CYLINDERS ',QM1, ' AND ', QP1, ' ARE NOT PIECED PROPERLY'
1612 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1613 ENDIF
1614 call mfix_exit(myPE)
1615 ENDIF
1616
1617 tan_half_angle = (R2-R1)/(x2-x1)
1618
1619 HALF_ANGLE(Q) = DATAN(tan_half_angle)/PI*180.0D0
1620 lambda_x(Q) = ONE
1621 lambda_y(Q) = -ONE/(tan_half_angle)**2
1622 lambda_z(Q) = -ONE/(tan_half_angle)**2
1623 dquadric(Q) = ZERO
1624
1625 piece_xmin(Q) = x1
1626 piece_xmax(Q) = x2
1627
1628 t_x(Q) = x1 - R1/tan_half_angle
1629 t_y(Q) = t_y(QM1)
1630 t_z(Q) = t_z(QM1)
1631
1632 IF(MyPE == PE_IO) THEN
1633 WRITE(*,*) ' QUADRIC:',Q, ' WAS DEFINED AS ', TRIM(QUADRIC_FORM(Q))
1634 WRITE(*,*) ' WITH AN HALF-ANGLE OF ', HALF_ANGLE(Q), 'DEG.'
1635 ENDIF
1636
1637
1638
1639 ELSEIF((TRIM(QUADRIC_FORM(QM1))=='Y_CYL_INT').AND. &
1640 (TRIM(QUADRIC_FORM(QP1))=='Y_CYL_INT')) THEN
1641
1642 (Q) = 'Y_CONE'
1643
1644 aligned = (t_x(QM1)==t_x(QP1)).AND.(t_z(QM1)==t_z(QP1))
1645 IF(.NOT.aligned) THEN
1646 IF(MyPE == PE_IO) THEN
1647 WRITE(*,*)' ERROR: CYLINDERS ',QM1, ' AND ', QP1, ' ARE NOT ALIGNED'
1648 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1649 ENDIF
1650 call mfix_exit(myPE)
1651 ENDIF
1652
1653 R1 = RADIUS(QM1)
1654 R2 = RADIUS(QP1)
1655 IF(R1==R2) THEN
1656 IF(MyPE == PE_IO) THEN
1657 WRITE(*,*)' ERROR: CYLINDERS ',QM1, ' AND ', QP1, ' HAVE THE SAME RADIUS'
1658 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1659 ENDIF
1660 call mfix_exit(myPE)
1661 ENDIF
1662
1663 y1 = piece_ymax(QM1)
1664 y2 = piece_ymin(QP1)
1665 IF(y2<=y1) THEN
1666 IF(MyPE == PE_IO) THEN
1667 WRITE(*,*)' ERROR: CYLINDERS ',QM1, ' AND ', QP1, ' ARE NOT PIECED PROPERLY'
1668 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1669 ENDIF
1670 call mfix_exit(myPE)
1671 ENDIF
1672
1673 tan_half_angle = (R2-R1)/(y2-y1)
1674
1675 HALF_ANGLE(Q) = DATAN(tan_half_angle)/PI*180.0D0
1676 lambda_x(Q) = ONE/(tan_half_angle)**2
1677 lambda_y(Q) = -ONE
1678 lambda_z(Q) = ONE/(tan_half_angle)**2
1679 dquadric(Q) = ZERO
1680
1681 piece_ymin(Q) = y1
1682 piece_ymax(Q) = y2
1683
1684 t_x(Q) = t_x(QM1)
1685 t_y(Q) = y1 - R1/tan_half_angle
1686 t_z(Q) = t_z(QM1)
1687
1688 IF(MyPE == PE_IO) THEN
1689 WRITE(*,*) ' QUADRIC:',Q, ' WAS DEFINED AS ', TRIM(QUADRIC_FORM(Q))
1690 WRITE(*,*) ' WITH AN HALF-ANGLE OF ', HALF_ANGLE(Q), 'DEG.'
1691 ENDIF
1692
1693 ELSEIF((TRIM(QUADRIC_FORM(QM1))=='Y_CYL_EXT').AND. &
1694 (TRIM(QUADRIC_FORM(QP1))=='Y_CYL_EXT')) THEN
1695
1696 (Q) = 'Y_CONE'
1697
1698 aligned = (t_x(QM1)==t_x(QP1)).AND.(t_z(QM1)==t_z(QP1))
1699 IF(.NOT.aligned) THEN
1700 IF(MyPE == PE_IO) THEN
1701 WRITE(*,*)' ERROR: CYLINDERS ',QM1, ' AND ', QP1, ' ARE NOT ALIGNED'
1702 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1703 ENDIF
1704 call mfix_exit(myPE)
1705 ENDIF
1706
1707 R1 = RADIUS(QM1)
1708 R2 = RADIUS(QP1)
1709 IF(R1==R2) THEN
1710 IF(MyPE == PE_IO) THEN
1711 WRITE(*,*)' ERROR: CYLINDERS ',QM1, ' AND ', QP1, ' HAVE THE SAME RADIUS'
1712 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1713 ENDIF
1714 call mfix_exit(myPE)
1715 ENDIF
1716
1717 y1 = piece_ymax(QM1)
1718 y2 = piece_ymin(QP1)
1719 IF(y2<=y1) THEN
1720 IF(MyPE == PE_IO) THEN
1721 WRITE(*,*)' ERROR: CYLINDERS ',QM1, ' AND ', QP1, ' ARE NOT PIECED PROPERLY'
1722 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1723 ENDIF
1724 call mfix_exit(myPE)
1725 ENDIF
1726
1727 tan_half_angle = (R2-R1)/(y2-y1)
1728
1729 HALF_ANGLE(Q) = DATAN(tan_half_angle)/PI*180.0D0
1730 lambda_x(Q) = -ONE/(tan_half_angle)**2
1731 lambda_y(Q) = ONE
1732 lambda_z(Q) = -ONE/(tan_half_angle)**2
1733 dquadric(Q) = ZERO
1734
1735 piece_ymin(Q) = y1
1736 piece_ymax(Q) = y2
1737
1738 t_x(Q) = t_x(QM1)
1739 t_y(Q) = y1 - R1/tan_half_angle
1740 t_z(Q) = t_z(QM1)
1741
1742 IF(MyPE == PE_IO) THEN
1743 WRITE(*,*) ' QUADRIC:',Q, ' WAS DEFINED AS ', TRIM(QUADRIC_FORM(Q))
1744 WRITE(*,*) ' WITH AN HALF-ANGLE OF ', HALF_ANGLE(Q), 'DEG.'
1745 ENDIF
1746
1747
1748
1749 ELSEIF((TRIM(QUADRIC_FORM(QM1))=='Z_CYL_INT').AND. &
1750 (TRIM(QUADRIC_FORM(QP1))=='Z_CYL_INT')) THEN
1751
1752 (Q) = 'Z_CONE'
1753
1754 aligned = (t_x(QM1)==t_x(QP1)).AND.(t_y(QM1)==t_y(QP1))
1755 IF(.NOT.aligned) THEN
1756 IF(MyPE == PE_IO) THEN
1757 WRITE(*,*)' ERROR: CYLINDERS ',QM1, ' AND ', QP1, ' ARE NOT ALIGNED'
1758 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1759 ENDIF
1760 call mfix_exit(myPE)
1761 ENDIF
1762
1763 R1 = RADIUS(QM1)
1764 R2 = RADIUS(QP1)
1765 IF(R1==R2) THEN
1766 IF(MyPE == PE_IO) THEN
1767 WRITE(*,*)' ERROR: CYLINDERS ',QM1, ' AND ', QP1, ' HAVE THE SAME RADIUS'
1768 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1769 ENDIF
1770 call mfix_exit(myPE)
1771 ENDIF
1772
1773 z1 = piece_zmax(QM1)
1774 z2 = piece_zmin(QP1)
1775 IF(z2<=z1) THEN
1776 IF(MyPE == PE_IO) THEN
1777 WRITE(*,*)' ERROR: CYLINDERS ',QM1, ' AND ', QP1, ' ARE NOT PIECED PROPERLY'
1778 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1779 ENDIF
1780 call mfix_exit(myPE)
1781 ENDIF
1782
1783 tan_half_angle = (R2-R1)/(z2-z1)
1784
1785 HALF_ANGLE(Q) = DATAN(tan_half_angle)/PI*180.0D0
1786 lambda_x(Q) = ONE/(tan_half_angle)**2
1787 lambda_y(Q) = ONE/(tan_half_angle)**2
1788 lambda_z(Q) = -ONE
1789 dquadric(Q) = ZERO
1790
1791 piece_zmin(Q) = z1
1792 piece_zmax(Q) = z2
1793
1794 t_x(Q) = t_x(QM1)
1795 t_y(Q) = t_y(QM1)
1796 t_z(Q) = z1 - R1/tan_half_angle
1797
1798 IF(MyPE == PE_IO) THEN
1799 WRITE(*,*) ' QUADRIC:',Q, ' WAS DEFINED AS ', TRIM(QUADRIC_FORM(Q))
1800 WRITE(*,*) ' WITH AN HALF-ANGLE OF ', HALF_ANGLE(Q), 'DEG.'
1801 ENDIF
1802
1803 ELSEIF((TRIM(QUADRIC_FORM(QM1))=='Z_CYL_EXT').AND. &
1804 (TRIM(QUADRIC_FORM(QP1))=='Z_CYL_EXT')) THEN
1805
1806 (Q) = 'Z_CONE'
1807
1808 aligned = (t_x(QM1)==t_x(QP1)).AND.(t_y(QM1)==t_y(QP1))
1809 IF(.NOT.aligned) THEN
1810 IF(MyPE == PE_IO) THEN
1811 WRITE(*,*)' ERROR: CYLINDERS ',QM1, ' AND ', QP1, ' ARE NOT ALIGNED'
1812 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1813 ENDIF
1814 call mfix_exit(myPE)
1815 ENDIF
1816
1817 R1 = RADIUS(QM1)
1818 R2 = RADIUS(QP1)
1819 IF(R1==R2) THEN
1820 IF(MyPE == PE_IO) THEN
1821 WRITE(*,*)' ERROR: CYLINDERS ',QM1, ' AND ', QP1, ' HAVE THE SAME RADIUS'
1822 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1823 ENDIF
1824 call mfix_exit(myPE)
1825 ENDIF
1826
1827 z1 = piece_zmax(QM1)
1828 z2 = piece_zmin(QP1)
1829 IF(z2<=z1) THEN
1830 IF(MyPE == PE_IO) THEN
1831 WRITE(*,*)' ERROR: CYLINDERS ',QM1, ' AND ', QP1, ' ARE NOT PIECED PROPERLY'
1832 WRITE(*,*)'PLEASE CORRECT MFIX.DAT AND TRY AGAIN.'
1833 ENDIF
1834 call mfix_exit(myPE)
1835 ENDIF
1836
1837 tan_half_angle = (R2-R1)/(z2-z1)
1838
1839 HALF_ANGLE(Q) = DATAN(tan_half_angle)/PI*180.0D0
1840 lambda_x(Q) = -ONE/(tan_half_angle)**2
1841 lambda_y(Q) = -ONE/(tan_half_angle)**2
1842 lambda_z(Q) = ONE
1843 dquadric(Q) = ZERO
1844
1845 piece_zmin(Q) = z1
1846 piece_zmax(Q) = z2
1847
1848 t_x(Q) = t_x(QM1)
1849 t_y(Q) = t_y(QM1)
1850 t_z(Q) = z1 - R1/tan_half_angle
1851
1852 IF(MyPE == PE_IO) THEN
1853 WRITE(*,*) ' QUADRIC:',Q, ' WAS DEFINED AS ', TRIM(QUADRIC_FORM(Q))
1854 WRITE(*,*) ' WITH AN HALF-ANGLE OF ', HALF_ANGLE(Q), 'DEG.'
1855 ENDIF
1856
1857
1858 ELSE
1859 IF(MyPE == PE_IO) THEN
1860 WRITE(*,*) ' ERROR: C2C MUST BE DEFINED BETWEEN 2 CYLINDERS'
1861 WRITE(*,*) ' QUADRIC:',QM1, ' IS ', TRIM(QUADRIC_FORM(QM1))
1862 WRITE(*,*) ' QUADRIC:',QP1, ' IS ', TRIM(QUADRIC_FORM(QP1))
1863 ENDIF
1864 call mfix_exit(myPE)
1865
1866 ENDIF
1867
1868 RETURN
1869 END SUBROUTINE BUILD_CONE_FOR_C2C
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880
1881
1882
1883
1884
1885 SUBROUTINE CG_FLOW_TO_VEL
1886
1887 USE bc
1888 USE compar
1889 USE constant
1890 USE cutcell
1891 USE eos, ONLY: EOSG
1892 USE fldvar
1893 USE functions
1894 USE funits
1895 USE geometry
1896 USE indices
1897 USE mpi_utility
1898 USE parallel
1899 USE param
1900 USE param1
1901 USE physprop
1902 USE quadric
1903 USE run
1904 USE scales
1905 USE sendrecv
1906 USE toleranc
1907 USE vtk
1908
1909 IMPLICIT NONE
1910
1911
1912
1913
1914
1915
1916
1917
1918
1919
1920
1921 INTEGER :: M, BCV
1922
1923 DOUBLE PRECISION :: VOLFLOW
1924
1925 DOUBLE PRECISION :: EPS
1926
1927 DOUBLE PRECISION :: MW
1928
1929
1930
1931 DO BCV = 1, DIMENSION_BC
1932
1933 IF (BC_TYPE(BCV)=='CG_MI') THEN
1934
1935 IF(BC_VELMAG_g(BCV)==UNDEFINED) THEN
1936
1937
1938
1939 IF (BC_MASSFLOW_G(BCV) /= UNDEFINED) THEN
1940 IF (RO_G0 /= UNDEFINED) THEN
1941 VOLFLOW = BC_MASSFLOW_G(BCV)/RO_G0
1942 ELSE
1943 IF (BC_P_G(BCV)/=UNDEFINED .AND. BC_T_G(BCV)/=UNDEFINED) &
1944 THEN
1945 IF (MW_AVG == UNDEFINED) THEN
1946 MW = CALC_MW(BC_X_G,DIMENSION_BC,BCV,NMAX(0),MW_G)
1947 ELSE
1948 MW = MW_AVG
1949 ENDIF
1950 VOLFLOW = BC_MASSFLOW_G(BCV)/EOSG(MW,(BC_P_G(BCV)-P_REF), &
1951 BC_T_G(BCV))
1952 ELSE
1953 IF (BC_TYPE(BCV) == 'CG_MO') THEN
1954 IF (BC_MASSFLOW_G(BCV) == ZERO) THEN
1955 VOLFLOW = ZERO
1956 ENDIF
1957 ELSE
1958 IF(DMP_LOG)WRITE (UNIT_LOG, 1020) BCV
1959 call mfix_exit(myPE)
1960 ENDIF
1961 ENDIF
1962 ENDIF
1963
1964
1965
1966 IF (BC_VOLFLOW_G(BCV) /= UNDEFINED) THEN
1967 IF (.NOT.COMPARE(VOLFLOW,BC_VOLFLOW_G(BCV))) THEN
1968 IF(DMP_LOG)WRITE (UNIT_LOG, 1000) BCV, VOLFLOW, BC_VOLFLOW_G(BCV)
1969 call mfix_exit(myPE)
1970 ENDIF
1971 ELSE
1972 BC_VOLFLOW_G(BCV) = VOLFLOW
1973 ENDIF
1974 ENDIF
1975
1976
1977
1978 IF (BC_VOLFLOW_G(BCV) /= UNDEFINED) THEN
1979 IF (BC_EP_G(BCV) /= UNDEFINED) THEN
1980 BC_VELMAG_g(BCV) = BC_VOLFLOW_G(BCV)/(BC_AREA(BCV)*BC_EP_G(BCV))
1981 ELSE
1982 RETURN
1983 ENDIF
1984 ENDIF
1985
1986 ENDIF
1987
1988
1989
1990
1991
1992 DO M = 1, MMAX
1993
1994 IF(BC_VELMAG_s(BCV,M)==UNDEFINED) THEN
1995
1996
1997
1998 IF (BC_MASSFLOW_S(BCV,M) /= UNDEFINED) THEN
1999 IF (RO_S0(M) /= UNDEFINED) THEN
2000 VOLFLOW = BC_MASSFLOW_S(BCV,M)/RO_S0(M)
2001 ELSE
2002 RETURN
2003 ENDIF
2004
2005
2006
2007 IF (BC_VOLFLOW_S(BCV,M) /= UNDEFINED) THEN
2008 IF (.NOT.COMPARE(VOLFLOW,BC_VOLFLOW_S(BCV,M))) THEN
2009 IF(DMP_LOG)WRITE(UNIT_LOG,1200)BCV,VOLFLOW,M,BC_VOLFLOW_S(BCV,M)
2010 call mfix_exit(myPE)
2011 ENDIF
2012 ELSE
2013 BC_VOLFLOW_S(BCV,M) = VOLFLOW
2014 ENDIF
2015 ENDIF
2016
2017 IF (BC_ROP_S(BCV,M)==UNDEFINED .AND. MMAX==1) BC_ROP_S(BCV,M)&
2018 = (ONE - BC_EP_G(BCV))*RO_S0(M)
2019 IF (BC_VOLFLOW_S(BCV,M) /= UNDEFINED) THEN
2020 IF (BC_ROP_S(BCV,M) /= UNDEFINED) THEN
2021 EPS = BC_ROP_S(BCV,M)/RO_S0(M)
2022 IF (EPS /= ZERO) THEN
2023 BC_VELMAG_s(BCV,M) = BC_VOLFLOW_S(BCV,M)/(BC_AREA(BCV)*EPS)
2024 ELSE
2025 IF (BC_VOLFLOW_S(BCV,M) == ZERO) THEN
2026 BC_VELMAG_s(BCV,M) = ZERO
2027 ELSE
2028 IF(DMP_LOG)WRITE (UNIT_LOG, 1250) BCV, M
2029 call mfix_exit(myPE)
2030 ENDIF
2031 ENDIF
2032 ELSE
2033 IF (BC_VOLFLOW_S(BCV,M) == ZERO) THEN
2034 BC_VELMAG_s(BCV,M) = ZERO
2035 ELSE
2036 IF(DMP_LOG)WRITE (UNIT_LOG, 1260) BCV, M
2037 call mfix_exit(myPE)
2038 ENDIF
2039 ENDIF
2040 ENDIF
2041
2042 ENDIF
2043 END DO
2044 ENDIF
2045 END DO
2046
2047 100 FORMAT(1X,A,I8)
2048 110 FORMAT(1X,A,A)
2049 120 FORMAT(1X,A,F14.8,/)
2050 130 FORMAT(1X,A,I8,F14.8,/)
2051
2052 1000 FORMAT(/1X,70('*')//' From: FLOW_TO_VEL',/' Message: BC No:',I2,/,&
2053 ' Computed volumetric flow is not equal to specified value',/,&
2054 ' Value computed from mass flow = ',G14.7,/,&
2055 ' Specified value (BC_VOLFLOW_g) = ',G14.7,/1X,70('*')/)
2056
2057
2058 1020 FORMAT(/1X,70('*')//' From: FLOW_TO_VEL',/' Message: BC No:',I2,&
2059 ' BC_P_g, BC_T_g, and BC_X_g',/' should be specified',/1X,70('*')/)
2060
2061
2062 1200 FORMAT(/1X,70('*')//' From: FLOW_TO_VEL',/' Message: BC No:',I2,/,&
2063 ' Computed volumetric flow is not equal to specified value',/,&
2064 ' Value computed from mass flow = ',G14.7,/,&
2065 ' Specified value (BC_VOLFLOW_s',I1,') = ',G14.7,/1X,70('*')/)
2066
2067 1250 FORMAT(/1X,70('*')//' From: FLOW_TO_VEL',/' Message: BC No:',I2,/,&
2068 ' Non-zero vol. or mass flow specified with BC_ROP_s',&
2069 I1,' = 0.',/1X,70('*')/)
2070 1260 FORMAT(/1X,70('*')//' From: FLOW_TO_VEL',/' Message: BC No:',I2,/,&
2071 ' BC_ROP_s',I1,' not specified',/1X,70('*')/)
2072 RETURN
2073
2074
2075 END SUBROUTINE CG_FLOW_TO_VEL
2076
2077
2078
2079
2080
2081
2082
2083
2084
2085
2086
2087
2088
2089
2090
2091 SUBROUTINE CONVERT_CG_MI_TO_PS
2092
2093 USE bc
2094 USE compar
2095 USE constant
2096 USE cutcell
2097 USE eos, only: EOSG
2098 USE fldvar
2099 USE functions
2100 USE funits
2101 USE geometry
2102 USE indices
2103 USE mpi_utility
2104 USE parallel
2105 USE param
2106 USE param1
2107 USE physprop
2108 USE ps
2109 USE quadric
2110 USE run
2111 USE scales
2112 USE sendrecv
2113 USE toleranc
2114 USE vtk
2115
2116 IMPLICIT NONE
2117
2118
2119
2120
2121
2122
2123
2124
2125
2126
2127
2128 INTEGER :: IJK, M, N, BCV
2129
2130 INTEGER :: iproc
2131 INTEGER :: NPS,PSV
2132
2133
2134
2135
2136
2137 CALL MPI_BARRIER(MPI_COMM_WORLD,mpierr)
2138
2139
2140
2141 do iproc = 0,NumPEs-1
2142 if (MyPE==iproc) Then
2143
2144
2145
2146 = 0
2147
2148 PS_LP: do PSV = 1, DIMENSION_PS
2149 if(.NOT.PS_DEFINED(PSV)) cycle PS_LP
2150 NPS = PSV
2151 enddo PS_LP
2152
2153
2154
2155
2156
2157 DO IJK = ijkstart3, ijkend3
2158 BCV = BC_ID(IJK)
2159 IF(BCV>0) THEN
2160 IF(CG_MI_CONVERTED_TO_PS(BCV).AND.INTERIOR_CELL_AT(IJK).AND.VOL(IJK)>ZERO) THEN
2161
2162 NPS = NPS + 1
2163
2164
2165
2166 (NPS) = .TRUE.
2167
2168 POINT_SOURCE = .TRUE.
2169
2170 PS_I_w(NPS) = I_OF(IJK)
2171 PS_I_e(NPS) = I_OF(IJK)
2172 PS_J_s(NPS) = J_OF(IJK)
2173 PS_J_n(NPS) = J_OF(IJK)
2174 PS_K_b(NPS) = K_OF(IJK)
2175 PS_K_t(NPS) = K_OF(IJK)
2176
2177 PS_VOLUME(NPS) = VOL(IJK)
2178
2179 PS_MASSFLOW_g(NPS) = BC_MASSFLOW_g(BCV) * VOL(IJK) / BC_VOL(BCV)
2180
2181 PS_T_g(NPS) = BC_T_g(BCV)
2182
2183 IF(BC_U_g(BCV)==UNDEFINED) THEN
2184 PS_U_g(NPS) = Normal_S(IJK,1)
2185 ELSE
2186 PS_U_g(NPS) = BC_U_g(BCV)
2187 ENDIF
2188
2189 IF(BC_V_g(BCV)==UNDEFINED) THEN
2190 PS_V_g(NPS) = Normal_S(IJK,2)
2191 ELSE
2192 PS_V_g(NPS) = BC_V_g(BCV)
2193 ENDIF
2194
2195 IF(BC_W_g(BCV)==UNDEFINED) THEN
2196 PS_W_g(NPS) = Normal_S(IJK,3)
2197 ELSE
2198 PS_W_g(NPS) = BC_W_g(BCV)
2199 ENDIF
2200
2201 DO N=1,NMAX(0)
2202 PS_X_g(NPS,N) = BC_X_g(BCV,N)
2203 ENDDO
2204
2205 DO M=1, MMAX
2206 PS_MASSFLOW_s(NPS,M) = BC_MASSFLOW_s(BCV,M) * VOL(IJK) / BC_VOL(BCV)
2207
2208 PS_T_s(NPS,1) = BC_T_s(BCV,M)
2209
2210 IF(BC_U_s(BCV,M)==UNDEFINED) THEN
2211 PS_U_s(NPS,M) = Normal_S(IJK,1)
2212 ELSE
2213 PS_U_s(NPS,M) = BC_U_s(BCV,M)
2214 ENDIF
2215
2216 IF(BC_V_s(BCV,M)==UNDEFINED) THEN
2217 PS_V_s(NPS,M) = Normal_S(IJK,2)
2218 ELSE
2219 PS_V_s(NPS,M) = BC_V_s(BCV,M)
2220 ENDIF
2221
2222 IF(BC_W_s(BCV,M)==UNDEFINED) THEN
2223 PS_W_s(NPS,M) = Normal_S(IJK,3)
2224 ELSE
2225 PS_W_s(NPS,M) = BC_W_s(BCV,M)
2226 ENDIF
2227
2228
2229 DO N=1,NMAX(M)
2230 PS_X_s(NPS,M,N) = BC_X_s(BCV,M,N)
2231 ENDDO
2232
2233 ENDDO
2234
2235
2236 ENDIF
2237 ENDIF
2238
2239 ENDDO
2240
2241 endif
2242
2243 CALL MPI_BARRIER(MPI_COMM_WORLD,mpierr)
2244 call bcast(POINT_SOURCE,iproc)
2245 call bcast(PS_DEFINED,iproc)
2246 call bcast(PS_I_w,iproc)
2247 call bcast(PS_I_e,iproc)
2248 call bcast(PS_J_s,iproc)
2249 call bcast(PS_J_n,iproc)
2250 call bcast(PS_K_b,iproc)
2251 call bcast(PS_K_t,iproc)
2252 call bcast(PS_MASSFLOW_g,iproc)
2253 call bcast(PS_U_g,iproc)
2254 call bcast(PS_V_g,iproc)
2255 call bcast(PS_W_g,iproc)
2256 call bcast(PS_X_g,iproc)
2257 call bcast(PS_T_g,iproc)
2258 call bcast(PS_MASSFLOW_s,iproc)
2259 call bcast(PS_U_s,iproc)
2260 call bcast(PS_V_s,iproc)
2261 call bcast(PS_W_s,iproc)
2262 call bcast(PS_X_s,iproc)
2263 call bcast(PS_T_s,iproc)
2264 call bcast(PS_VOLUME,iproc)
2265
2266 enddo
2267
2268 CALL MPI_BARRIER(MPI_COMM_WORLD,mpierr)
2269
2270
2271
2272 RETURN
2273
2274 END SUBROUTINE CONVERT_CG_MI_TO_PS
2275
2276
2277
2278
2279
2280
2281
2282
2283
2284
2285
2286
2287
2288
2289
2290 SUBROUTINE CONVERT_CG_MI_TO_PS_PE
2291
2292 USE bc
2293 USE compar
2294 USE constant
2295 USE cutcell
2296 USE eos, ONLY: EOSG
2297 USE fldvar
2298 USE functions
2299 USE funits
2300 USE geometry
2301 USE indices
2302 USE mpi_utility
2303 USE parallel
2304 USE param
2305 USE param1
2306 USE physprop
2307 USE ps
2308 USE quadric
2309 USE run
2310 USE scales
2311 USE sendrecv
2312 USE toleranc
2313 USE vtk
2314
2315 IMPLICIT NONE
2316
2317
2318
2319
2320
2321
2322
2323
2324
2325
2326
2327 INTEGER :: IJK, BCV
2328
2329 INTEGER :: NPS,PSV
2330
2331
2332
2333
2334
2335
2336
2337
2338
2339
2340
2341
2342
2343
2344
2345
2346 = 0
2347
2348 PS_LP: do PSV = 1, DIMENSION_PS
2349 if(.NOT.PS_DEFINED(PSV)) cycle PS_LP
2350 NPS = PSV
2351 enddo PS_LP
2352
2353 print *,'Last PS=',NPS
2354
2355
2356
2357
2358
2359 DO IJK = ijkstart3, ijkend3
2360 BCV = BC_ID(IJK)
2361 IF(BCV>0) THEN
2362 IF(CG_MI_CONVERTED_TO_PS(BCV).AND.INTERIOR_CELL_AT(IJK).AND.VOL(IJK)>ZERO) THEN
2363
2364 NPS = NPS + 1
2365
2366 PS_DEFINED(NPS) = .TRUE.
2367
2368 POINT_SOURCE = .TRUE.
2369
2370 PS_I_w(NPS) = I_OF(IJK)
2371 PS_I_e(NPS) = I_OF(IJK)
2372 PS_J_s(NPS) = J_OF(IJK)
2373 PS_J_n(NPS) = J_OF(IJK)
2374 PS_K_b(NPS) = K_OF(IJK)
2375 PS_K_t(NPS) = K_OF(IJK)
2376
2377 PS_MASSFLOW_g(NPS) = BC_MASSFLOW_g(BCV) * VOL(IJK) / BC_VOL(BCV)
2378
2379 PS_VOLUME(NPS) = VOL(IJK)
2380
2381 PS_T_g(NPS) = BC_T_g(BCV)
2382
2383 IF(BC_U_g(NPS)==UNDEFINED) THEN
2384 PS_U_g(NPS) = Normal_S(IJK,1)
2385 ELSE
2386 PS_U_g(NPS) = BC_U_g(NPS)
2387 ENDIF
2388
2389 IF(BC_V_g(NPS)==UNDEFINED) THEN
2390 PS_V_g(NPS) = Normal_S(IJK,2)
2391 ELSE
2392 PS_V_g(NPS) = BC_V_g(NPS)
2393 ENDIF
2394
2395 IF(BC_W_g(NPS)==UNDEFINED) THEN
2396 PS_W_g(NPS) = Normal_S(IJK,3)
2397 ELSE
2398 PS_W_g(NPS) = BC_W_g(NPS)
2399 ENDIF
2400
2401
2402 (NPS,1) = 0.0
2403
2404 PS_T_s(NPS,1) = 298.0
2405
2406 PS_U_s(NPS,1) = Normal_S(IJK,1)
2407 PS_V_s(NPS,1) = Normal_S(IJK,2)
2408 PS_W_s(NPS,1) = Normal_S(IJK,3)
2409
2410 PS_U_s(NPS,1) = ZERO
2411 PS_V_s(NPS,1) = ZERO
2412 PS_W_s(NPS,1) = ZERO
2413
2414
2415
2416
2417
2418
2419
2420
2421
2422
2423
2424
2425
2426
2427
2428
2429
2430
2431
2432
2433
2434
2435
2436
2437
2438
2439 print*,'PS created:',NPS,PS_MASSFLOW_g(NPS),PS_VOLUME(NPS),PS_I_w(NPS),PS_J_n(NPS),PS_K_b(NPS), &
2440 INTERIOR_CELL_AT(IJK),PS_U_g(NPS),PS_V_g(NPS),PS_W_g(NPS), &
2441 CUT_CELL_AT(IJK),Normal_S(IJK,1),Normal_S(IJK,2),Normal_S(IJK,3)
2442
2443
2444
2445 ENDIF
2446 ENDIF
2447 ENDDO
2448
2449
2450
2451
2452
2453
2454
2455
2456 FORMAT(1X,A,I8)
2457 110 FORMAT(1X,A,A)
2458 120 FORMAT(1X,A,F14.8,/)
2459 130 FORMAT(1X,A,I8,F14.8,/)
2460
2461
2462 1000 FORMAT(/1X,70('*')//' From: FLOW_TO_VEL',/' Message: BC No:',I2,/,&
2463 ' Computed volumetric flow is not equal to specified value',/,&
2464 ' Value computed from mass flow = ',G14.7,/,&
2465 ' Specified value (BC_VOLFLOW_g) = ',G14.7,/1X,70('*')/)
2466
2467
2468 1020 FORMAT(/1X,70('*')//' From: FLOW_TO_VEL',/' Message: BC No:',I2,&
2469 ' BC_P_g, BC_T_g, and BC_X_g',/' should be specified',/1X,70('*')/)
2470
2471
2472 1200 FORMAT(/1X,70('*')//' From: FLOW_TO_VEL',/' Message: BC No:',I2,/,&
2473 ' Computed volumetric flow is not equal to specified value',/,&
2474 ' Value computed from mass flow = ',G14.7,/,&
2475 ' Specified value (BC_VOLFLOW_s',I1,') = ',G14.7,/1X,70('*')/)
2476
2477 1250 FORMAT(/1X,70('*')//' From: FLOW_TO_VEL',/' Message: BC No:',I2,/,&
2478 ' Non-zero vol. or mass flow specified with BC_ROP_s',&
2479 I1,' = 0.',/1X,70('*')/)
2480 1260 FORMAT(/1X,70('*')//' From: FLOW_TO_VEL',/' Message: BC No:',I2,/,&
2481 ' BC_ROP_s',I1,' not specified',/1X,70('*')/)
2482 RETURN
2483
2484
2485 END SUBROUTINE CONVERT_CG_MI_TO_PS_PE
2486
2487
2488
2489
2490
2491
2492
2493
2494
2495
2496
2497
2498
2499
2500
2501
2502
2503 SUBROUTINE GET_DXYZ_FROM_CONTROL_POINTS
2504
2505
2506
2507
2508 USE param
2509 USE param1
2510 USE constant
2511 USE run
2512 USE physprop
2513 USE indices
2514 USE scalars
2515 USE funits
2516 USE leqsol
2517 USE compar
2518 USE mpi_utility
2519 USE bc
2520 USE DISCRETELEMENT
2521
2522 USE cutcell
2523 USE quadric
2524 USE vtk
2525 USE polygon
2526 USE dashboard
2527 USE stl
2528
2529
2530 IMPLICIT NONE
2531
2532
2533
2534
2535
2536
2537
2538
2539
2540
2541 INTEGER :: N,NX,NY,NZ
2542 INTEGER :: I,I1,I2,J,J1,J2,K,K1,K2
2543 DOUBLE PRECISION :: L,CELL_RATIO
2544
2545 LOGICAL,DIMENSION(MAX_CP) :: INDEPENDENT_SEGMENT
2546
2547
2548
2549
2550
2551
2552
2553
2554
2555
2556
2557
2558 DO N = MAX_CP,1,-1
2559 CPX(N) = CPX(N-1)
2560 ENDDO
2561
2562 CPX(0) = ZERO
2563
2564
2565
2566 = 0
2567 DO N = 1,MAX_CP
2568 IF(CPX(N)>ZERO) NX = NX + 1
2569 ENDDO
2570
2571 IF(NX>0) THEN
2572 IF(MyPE==0) WRITE(*,*)' INFO: DEFINING GRID SPACING IN X-DIRECTION... '
2573 IF(MyPE==0) WRITE(*,*)' INFO: NUMBER OF CONTROL POINTS IN X-DIRECTION = ',NX
2574 IF(CPX(NX)/=XLENGTH) THEN
2575 IF(MyPE==0) WRITE(*,*)' ERROR: LAST CONTROL POINT MUST BE EQUAL TO XLENGTH.'
2576 IF(MyPE==0) WRITE(*,*)' XLENGTH = ',XLENGTH
2577 IF(MyPE==0) WRITE(*,*)' LAST CONTROL POINT = ',CPX(NX)
2578 call mfix_exit(myPE)
2579 ENDIF
2580 ENDIF
2581
2582
2583
2584
2585
2586 = .TRUE.
2587
2588 DO N = 1,NX
2589
2590 IF(CPX(N) <= CPX(N-1)) THEN
2591 IF(MyPE==0) WRITE(*,*)' ERROR: CONTROL POINTS ALONG X MUST BE SORTED IN ASCENDING ORDER.'
2592 IF(MyPE==0) WRITE(*,*)' CPX = ',CPX(0:NX)
2593 call mfix_exit(myPE)
2594 ENDIF
2595
2596 IF(NCX(N) <= 1) THEN
2597 IF(MyPE==0) WRITE(*,*)' ERROR: NUMBER OF CELLS MUST BE LARGER THAN 1 IN X-SEGMENT :',N
2598 IF(MyPE==0) WRITE(*,*)' NCX = ',NCX(N)
2599 call mfix_exit(myPE)
2600 ENDIF
2601
2602 IF(ERX(N) <= ZERO) THEN
2603 IF(MyPE==0) WRITE(*,*)' ERROR: EXPANSION RATIO MUST BE POSITIVE IN X-SEGMENT :',N
2604 IF(MyPE==0) WRITE(*,*)' ERX = ',ERX(N)
2605 call mfix_exit(myPE)
2606 ENDIF
2607
2608 ENDDO
2609
2610 DO N = 1,NX
2611
2612 IF(FIRST_DX(N)/=ZERO.AND.LAST_DX(N)/=ZERO) THEN
2613 IF(MyPE==0) WRITE(*,*)' ERROR: FIRST AND LAST DX ARE DEFINED, WHICH IS NOT ALLOWED IN X-SEGMENT :',N
2614 IF(MyPE==0) WRITE(*,*)' FIRST DX = ',FIRST_DX(N)
2615 IF(MyPE==0) WRITE(*,*)' LAST DX = ',LAST_DX(N)
2616 call mfix_exit(myPE)
2617 ELSEIF(FIRST_DX(N)>ZERO) THEN
2618 IF(MyPE==0) WRITE(*,*)' INFO: FIRST DX DEFINED IN X-SEGMENT :',N
2619 IF(MyPE==0) WRITE(*,*)' FIRST DX = ',FIRST_DX(N)
2620 L = CPX(N) - CPX(N-1)
2621 IF(L<=FIRST_DX(N)+TOL_F) THEN
2622 IF(MyPE==0) WRITE(*,*)' ERROR: FIRST DX IS NOT SMALLER THAN SEGMENT LENGTH IN X-SEGMENT :',N
2623 IF(MyPE==0) WRITE(*,*)' FIRST DX = ',FIRST_DX(N)
2624 IF(MyPE==0) WRITE(*,*)' SEGMENT LENGTH = ',L
2625 call mfix_exit(myPE)
2626 ENDIF
2627 CALL FIND_CELL_RATIO('FIRST',FIRST_DX(N),L,NCX(N),CELL_RATIO)
2628 ERX(N) = CELL_RATIO**(NCX(N)-1)
2629 IF(MyPE==0) WRITE(*,*)' CORRESPONDING EXPANSION RATIO = ',ERX(N)
2630 ELSEIF(LAST_DX(N)>ZERO) THEN
2631 IF(MyPE==0) WRITE(*,*)' INFO: LAST DX DEFINED IN X-SEGMENT :',N
2632 IF(MyPE==0) WRITE(*,*)' LAST DX = ',LAST_DX(N)
2633 L = CPX(N) - CPX(N-1)
2634 IF(L<=LAST_DX(N)+TOL_F) THEN
2635 IF(MyPE==0) WRITE(*,*)' ERROR: LAST DX IS NOT SMALLER THAN SEGMENT LENGTH IN X-SEGMENT :',N
2636 IF(MyPE==0) WRITE(*,*)' LAST DX = ',LAST_DX(N)
2637 IF(MyPE==0) WRITE(*,*)' SEGMENT LENGTH = ',L
2638 call mfix_exit(myPE)
2639 ENDIF
2640 CALL FIND_CELL_RATIO('LAST ',LAST_DX(N),L,NCX(N),CELL_RATIO)
2641 ERX(N) = CELL_RATIO**(NCX(N)-1)
2642 IF(MyPE==0) WRITE(*,*)' CORRESPONDING EXPANSION RATIO = ',ERX(N)
2643 ELSEIF(FIRST_DX(N)<ZERO) THEN
2644 IF(N==1) THEN
2645 IF(MyPE==0) WRITE(*,*)' ERROR: FIRST DX CANNOT MATCH PREVIOUS DX FOR FIRST SEGMENT.'
2646 call mfix_exit(myPE)
2647 ELSE
2648 IF(MyPE==0) WRITE(*,*)' INFO: FIRST DX WILL ATTEMPT TO MATCH PREVIOUS DX FOR SEGMENT :',N
2649 INDEPENDENT_SEGMENT(N) = .FALSE.
2650 ENDIF
2651 ELSEIF(LAST_DX(N)<ZERO) THEN
2652 IF(N==NX) THEN
2653 IF(MyPE==0) WRITE(*,*)' ERROR: LAST DX CANNOT MATCH NEXT DX FOR LAST SEGMENT.'
2654 call mfix_exit(myPE)
2655 ELSE
2656 IF(MyPE==0) WRITE(*,*)' INFO: LAST DX WILL ATTEMPT TO MATCH NEXT DX FOR SEGMENT :',N
2657 INDEPENDENT_SEGMENT(N) = .FALSE.
2658 ENDIF
2659 ENDIF
2660
2661 ENDDO
2662
2663
2664
2665
2666
2667
2668 = 0
2669 = 0
2670
2671 DO N = 1,NX
2672
2673 = I1 + NCX(N) - 1
2674
2675 IF(INDEPENDENT_SEGMENT(N)) THEN
2676
2677 L = CPX(N) - CPX(N-1)
2678
2679 IF(ERX(N)/=ONE) THEN
2680 CELL_RATIO = ERX(N)**(ONE/DBLE(NCX(N)-1))
2681 (I1) = L * (ONE - CELL_RATIO) / (ONE - CELL_RATIO**NCX(N))
2682
2683 DO I = I1+1,I2
2684 (I) = DX(I-1) * CELL_RATIO
2685 ENDDO
2686
2687 ELSE
2688 DX(I1:I2) = L / NCX(N)
2689 ENDIF
2690
2691 ENDIF
2692
2693 I1 = I2 + 1
2694
2695 ENDDO
2696
2697
2698
2699
2700 = 0
2701 = 0
2702
2703 DO N = 1,NX
2704
2705 = I1 + NCX(N) - 1
2706
2707 IF(.NOT.INDEPENDENT_SEGMENT(N)) THEN
2708
2709 L = CPX(N) - CPX(N-1)
2710
2711 IF(FIRST_DX(N)<ZERO) THEN
2712 DX(I1) = DX(I1-1)
2713 CALL FIND_CELL_RATIO('FIRST',DX(I1),L,NCX(N),CELL_RATIO)
2714 DO I = I1+1,I2
2715 (I) = DX(I-1) * CELL_RATIO
2716 ENDDO
2717 ELSEIF(LAST_DX(N)<ZERO) THEN
2718 DX(I2) = DX(I2+1)
2719 CALL FIND_CELL_RATIO('LAST ',DX(I2),L,NCX(N),CELL_RATIO)
2720 DO I = I2-1,I1,-1
2721 (I) = DX(I+1) / CELL_RATIO
2722 ENDDO
2723 ENDIF
2724
2725 ENDIF
2726
2727 I1 = I2 + 1
2728
2729 ENDDO
2730
2731
2732
2733
2734 IF(I1>0.AND.I1/=IMAX) THEN
2735 IF(MyPE==0) WRITE(*,*)' ERROR: IMAX MUST BE EQUAL TO THE SUM OF NCX.'
2736 IF(MyPE==0) WRITE(*,*)' IMAX = ', IMAX
2737 IF(MyPE==0) WRITE(*,*)' SUM OF NCX = ', I1
2738 call mfix_exit(myPE)
2739 ENDIF
2740
2741
2742
2743
2744
2745
2746
2747
2748
2749
2750 DO N = MAX_CP,1,-1
2751 CPY(N) = CPY(N-1)
2752 ENDDO
2753
2754 CPY(0) = ZERO
2755
2756
2757
2758 = 0
2759 DO N = 1,MAX_CP
2760 IF(CPY(N)>ZERO) NY = NY + 1
2761 ENDDO
2762
2763 IF(NY>0) THEN
2764 IF(MyPE==0) WRITE(*,*)' INFO: DEFINING GRID SPACING IN Y-DIRECTION... '
2765 IF(MyPE==0) WRITE(*,*)' INFO: NUMBER OF CONTROL POINTS IN Y-DIRECTION = ',NY
2766 IF(CPY(NY)/=YLENGTH) THEN
2767 IF(MyPE==0) WRITE(*,*)' ERROR: LAST CONTROL POINT MUST BE EQUAL TO YLENGTH.'
2768 IF(MyPE==0) WRITE(*,*)' YLENGTH = ',YLENGTH
2769 IF(MyPE==0) WRITE(*,*)' LAST CONTROL POINT = ',CPY(NY)
2770 call mfix_exit(myPE)
2771 ENDIF
2772 ENDIF
2773
2774
2775
2776
2777
2778 = .TRUE.
2779
2780 DO N = 1,NY
2781
2782 IF(CPY(N) <= CPY(N-1)) THEN
2783 IF(MyPE==0) WRITE(*,*)' ERROR: CONTROL POINTS ALONG Y MUST BE SORTED IN ASCENDING ORDER.'
2784 IF(MyPE==0) WRITE(*,*)' CPY = ',CPY(0:NY)
2785 call mfix_exit(myPE)
2786 ENDIF
2787
2788 IF(NCY(N) <= 1) THEN
2789 IF(MyPE==0) WRITE(*,*)' ERROR: NUMBER OF CELLS MUST BE LARGER THAN 1 IN Y-SEGMENT :',N
2790 IF(MyPE==0) WRITE(*,*)' NCY = ',NCY(N)
2791 call mfix_exit(myPE)
2792 ENDIF
2793
2794 IF(ERY(N) <= ZERO) THEN
2795 IF(MyPE==0) WRITE(*,*)' ERROR: EXPANSION RATIO MUST BE POSITIVE IN Y-SEGMENT :',N
2796 IF(MyPE==0) WRITE(*,*)' ERY = ',ERY(N)
2797 call mfix_exit(myPE)
2798 ENDIF
2799
2800 ENDDO
2801
2802 DO N = 1,NY
2803
2804 IF(FIRST_DY(N)/=ZERO.AND.LAST_DY(N)/=ZERO) THEN
2805 IF(MyPE==0) WRITE(*,*)' ERROR: FIRST AND LAST DY ARE DEFINED, WHICH IS NOT ALLOWED IN Y-SEGMENT :',N
2806 IF(MyPE==0) WRITE(*,*)' FIRST DY = ',FIRST_DY(N)
2807 IF(MyPE==0) WRITE(*,*)' LAST DY = ',LAST_DY(N)
2808 call mfix_exit(myPE)
2809 ELSEIF(FIRST_DY(N)>ZERO) THEN
2810 IF(MyPE==0) WRITE(*,*)' INFO: FIRST DY DEFINED IN Y-SEGMENT :',N
2811 IF(MyPE==0) WRITE(*,*)' FIRST DY = ',FIRST_DY(N)
2812 L = CPY(N) - CPY(N-1)
2813 IF(L<=FIRST_DY(N)+TOL_F) THEN
2814 IF(MyPE==0) WRITE(*,*)' ERROR: FIRST DY IS NOT SMALLER THAN SEGMENT LENGTH IN Y-SEGMENT :',N
2815 IF(MyPE==0) WRITE(*,*)' FIRST DY = ',FIRST_DY(N)
2816 IF(MyPE==0) WRITE(*,*)' SEGMENT LENGTH = ',L
2817 call mfix_exit(myPE)
2818 ENDIF
2819 CALL FIND_CELL_RATIO('FIRST',FIRST_DY(N),L,NCY(N),CELL_RATIO)
2820 ERY(N) = CELL_RATIO**(NCY(N)-1)
2821 IF(MyPE==0) WRITE(*,*)' CORRESPONDING EXPANSION RATIO = ',ERY(N)
2822 ELSEIF(LAST_DY(N)>ZERO) THEN
2823 IF(MyPE==0) WRITE(*,*)' INFO: LAST DY DEFINED IN Y-SEGMENT :',N
2824 IF(MyPE==0) WRITE(*,*)' LAST DY = ',LAST_DY(N)
2825 L = CPY(N) - CPY(N-1)
2826 IF(L<=LAST_DY(N)+TOL_F) THEN
2827 IF(MyPE==0) WRITE(*,*)' ERROR: LAST DY IS NOT SMALLER THAN SEGMENT LENGTH IN Y-SEGMENT :',N
2828 IF(MyPE==0) WRITE(*,*)' LAST DY = ',LAST_DY(N)
2829 IF(MyPE==0) WRITE(*,*)' SEGMENT LENGTH = ',L
2830 call mfix_exit(myPE)
2831 ENDIF
2832 CALL FIND_CELL_RATIO('LAST ',LAST_DY(N),L,NCY(N),CELL_RATIO)
2833 ERY(N) = CELL_RATIO**(NCY(N)-1)
2834 IF(MyPE==0) WRITE(*,*)' CORRESPONDING EXPANSION RATIO = ',ERY(N)
2835 ELSEIF(FIRST_DY(N)<ZERO) THEN
2836 IF(N==1) THEN
2837 IF(MyPE==0) WRITE(*,*)' ERROR: FIRST DY CANNOT MATCH PREVIOUS DY FOR FIRST SEGMENT.'
2838 call mfix_exit(myPE)
2839 ELSE
2840 IF(MyPE==0) WRITE(*,*)' INFO: FIRST DY WILL ATTEMPT TO MATCH PREVIOUS DY FOR SEGMENT :',N
2841 INDEPENDENT_SEGMENT(N) = .FALSE.
2842 ENDIF
2843 ELSEIF(LAST_DY(N)<ZERO) THEN
2844 IF(N==NY) THEN
2845 IF(MyPE==0) WRITE(*,*)' ERROR: LAST DY CANNOT MATCH NEXT DY FOR LAST SEGMENT.'
2846 call mfix_exit(myPE)
2847 ELSE
2848 IF(MyPE==0) WRITE(*,*)' INFO: LAST DY WILL ATTEMPT TO MATCH NEXT DY FOR SEGMENT :',N
2849 INDEPENDENT_SEGMENT(N) = .FALSE.
2850 ENDIF
2851 ENDIF
2852
2853 ENDDO
2854
2855
2856
2857
2858
2859
2860 = 0
2861 = 0
2862
2863 DO N = 1,NY
2864
2865 = J1 + NCY(N) - 1
2866
2867 IF(INDEPENDENT_SEGMENT(N)) THEN
2868
2869 L = CPY(N) - CPY(N-1)
2870
2871 IF(ERY(N)/=ONE) THEN
2872 CELL_RATIO = ERY(N)**(ONE/DBLE(NCY(N)-1))
2873 (J1) = L * (ONE - CELL_RATIO) / (ONE - CELL_RATIO**NCY(N))
2874
2875 DO J = J1+1,J2
2876 (J) = DY(J-1) * CELL_RATIO
2877 ENDDO
2878
2879 ELSE
2880 DY(J1:J2) = L / NCY(N)
2881 ENDIF
2882
2883 ENDIF
2884
2885 J1 = J2 + 1
2886
2887 ENDDO
2888
2889
2890
2891
2892 = 0
2893 = 0
2894
2895 DO N = 1,NY
2896
2897 = J1 + NCY(N) - 1
2898
2899 IF(.NOT.INDEPENDENT_SEGMENT(N)) THEN
2900
2901 L = CPY(N) - CPY(N-1)
2902
2903 IF(FIRST_DY(N)<ZERO) THEN
2904 DY(J1) = DY(J1-1)
2905 CALL FIND_CELL_RATIO('FIRST',DY(J1),L,NCY(N),CELL_RATIO)
2906 DO J = J1+1,J2
2907 (J) = DY(J-1) * CELL_RATIO
2908 ENDDO
2909 ELSEIF(LAST_DY(N)<ZERO) THEN
2910 DY(J2) = DY(J2+1)
2911 CALL FIND_CELL_RATIO('LAST ',DY(J2),L,NCY(N),CELL_RATIO)
2912 DO J = J2-1,J1,-1
2913 (J) = DY(J+1) / CELL_RATIO
2914 ENDDO
2915 ENDIF
2916
2917 ENDIF
2918
2919 J1 = J2 + 1
2920
2921 ENDDO
2922
2923
2924
2925
2926 IF(J1>0.AND.J1/=JMAX) THEN
2927 IF(MyPE==0) WRITE(*,*)' ERROR: JMAX MUST BE EQUAL TO THE SUM OF NCY.'
2928 IF(MyPE==0) WRITE(*,*)' JMAX = ', JMAX
2929 IF(MyPE==0) WRITE(*,*)' SUM OF NCY = ', J1
2930 call mfix_exit(myPE)
2931 ENDIF
2932
2933
2934
2935
2936
2937
2938 IF(NO_K) RETURN
2939
2940
2941
2942
2943
2944 DO N = MAX_CP,1,-1
2945 CPZ(N) = CPZ(N-1)
2946 ENDDO
2947
2948 CPZ(0) = ZERO
2949
2950
2951
2952 = 0
2953 DO N = 1,MAX_CP
2954 IF(CPZ(N)>ZERO) NZ = NZ + 1
2955 ENDDO
2956
2957 IF(NZ>0) THEN
2958 IF(MyPE==0) WRITE(*,*)' INFO: DEFINING GRID SPACING IN Z-DIRECTION... '
2959 IF(MyPE==0) WRITE(*,*)' INFO: NUMBER OF CONTROL POINTS IN Z-DIRECTION = ',NZ
2960 IF(CPZ(NZ)/=ZLENGTH) THEN
2961 IF(MyPE==0) WRITE(*,*)' ERROR: LAST CONTROL POINT MUST BE EQUAL TO ZLENGTH.'
2962 IF(MyPE==0) WRITE(*,*)' ZLENGTH = ',ZLENGTH
2963 IF(MyPE==0) WRITE(*,*)' LAST CONTROL POINT = ',CPZ(NZ)
2964 call mfix_exit(myPE)
2965 ENDIF
2966 ENDIF
2967
2968
2969
2970
2971
2972 = .TRUE.
2973
2974 DO N = 1,NZ
2975
2976 IF(CPZ(N) <= CPZ(N-1)) THEN
2977 IF(MyPE==0) WRITE(*,*)' ERROR: CONTROL POINTS ALONG Z MUST BE SORTED IN ASCENDING ORDER.'
2978 IF(MyPE==0) WRITE(*,*)' CPZ = ',CPZ(0:NZ)
2979 call mfix_exit(myPE)
2980 ENDIF
2981
2982 IF(NCZ(N) <= 1) THEN
2983 IF(MyPE==0) WRITE(*,*)' ERROR: NUMBER OF CELLS MUST BE LARGER THAN 1 IN Z-SEGMENT :',N
2984 IF(MyPE==0) WRITE(*,*)' NCZ = ',NCZ(N)
2985 call mfix_exit(myPE)
2986 ENDIF
2987
2988 IF(ERZ(N) <= ZERO) THEN
2989 IF(MyPE==0) WRITE(*,*)' ERROR: EXPANSION RATIO MUST BE POSITIVE IN Z-SEGMENT :',N
2990 IF(MyPE==0) WRITE(*,*)' ERZ = ',ERZ(N)
2991 call mfix_exit(myPE)
2992 ENDIF
2993
2994 ENDDO
2995
2996 DO N = 1,NZ
2997
2998 IF(FIRST_DZ(N)/=ZERO.AND.LAST_DZ(N)/=ZERO) THEN
2999 IF(MyPE==0) WRITE(*,*)' ERROR: FIRST AND LAST DZ ARE DEFINED, WHICH IS NOT ALLOWED IN Z-SEGMENT :',N
3000 IF(MyPE==0) WRITE(*,*)' FIRST DZ = ',FIRST_DZ(N)
3001 IF(MyPE==0) WRITE(*,*)' LAST DZ = ',LAST_DZ(N)
3002 call mfix_exit(myPE)
3003 ELSEIF(FIRST_DZ(N)>ZERO) THEN
3004 IF(MyPE==0) WRITE(*,*)' INFO: FIRST DZ DEFINED IN Z-SEGMENT :',N
3005 IF(MyPE==0) WRITE(*,*)' FIRST DZ = ',FIRST_DZ(N)
3006 L = CPZ(N) - CPZ(N-1)
3007 IF(L<=FIRST_DZ(N)+TOL_F) THEN
3008 IF(MyPE==0) WRITE(*,*)' ERROR: FIRST DZ IS NOT SMALLER THAN SEGMENT LENGTH IN Z-SEGMENT :',N
3009 IF(MyPE==0) WRITE(*,*)' FIRST DZ = ',FIRST_DZ(N)
3010 IF(MyPE==0) WRITE(*,*)' SEGMENT LENGTH = ',L
3011 call mfix_exit(myPE)
3012 ENDIF
3013 CALL FIND_CELL_RATIO('FIRST',FIRST_DZ(N),L,NCZ(N),CELL_RATIO)
3014 ERZ(N) = CELL_RATIO**(NCZ(N)-1)
3015 IF(MyPE==0) WRITE(*,*)' CORRESPONDING EXPANSION RATIO = ',ERZ(N)
3016 ELSEIF(LAST_DZ(N)>ZERO) THEN
3017 IF(MyPE==0) WRITE(*,*)' INFO: LAST DZ DEFINED IN Z-SEGMENT :',N
3018 IF(MyPE==0) WRITE(*,*)' LAST DZ = ',LAST_DZ(N)
3019 L = CPZ(N) - CPZ(N-1)
3020 IF(L<=LAST_DZ(N)+TOL_F) THEN
3021 IF(MyPE==0) WRITE(*,*)' ERROR: LAST DZ IS NOT SMALLER THAN SEGMENT LENGTH IN Z-SEGMENT :',N
3022 IF(MyPE==0) WRITE(*,*)' LAST DZ = ',LAST_DZ(N)
3023 IF(MyPE==0) WRITE(*,*)' SEGMENT LENGTH = ',L
3024 call mfix_exit(myPE)
3025 ENDIF
3026 CALL FIND_CELL_RATIO('LAST ',LAST_DZ(N),L,NCZ(N),CELL_RATIO)
3027 ERZ(N) = CELL_RATIO**(NCZ(N)-1)
3028 IF(MyPE==0) WRITE(*,*)' CORRESPONDING EXPANSION RATIO = ',ERZ(N)
3029 ELSEIF(FIRST_DZ(N)<ZERO) THEN
3030 IF(N==1) THEN
3031 IF(MyPE==0) WRITE(*,*)' ERROR: FIRST DZ CANNOT MATCH PREVIOUS DZ FOR FIRST SEGMENT.'
3032 call mfix_exit(myPE)
3033 ELSE
3034 IF(MyPE==0) WRITE(*,*)' INFO: FIRST DZ WILL ATTEMPT TO MATCH PREVIOUS DZ FOR SEGMENT :',N
3035 INDEPENDENT_SEGMENT(N) = .FALSE.
3036 ENDIF
3037 ELSEIF(LAST_DZ(N)<ZERO) THEN
3038 IF(N==NZ) THEN
3039 IF(MyPE==0) WRITE(*,*)' ERROR: LAST DZ CANNOT MATCH NEXT DZ FOR LAST SEGMENT.'
3040 call mfix_exit(myPE)
3041 ELSE
3042 IF(MyPE==0) WRITE(*,*)' INFO: LAST DZ WILL ATTEMPT TO MATCH NEXT DZ FOR SEGMENT :',N
3043 INDEPENDENT_SEGMENT(N) = .FALSE.
3044 ENDIF
3045 ENDIF
3046
3047 ENDDO
3048
3049
3050
3051
3052
3053
3054 = 0
3055 = 0
3056
3057 DO N = 1,NZ
3058
3059 = K1 + NCZ(N) - 1
3060
3061 IF(INDEPENDENT_SEGMENT(N)) THEN
3062
3063 L = CPZ(N) - CPZ(N-1)
3064
3065 IF(ERZ(N)/=ONE) THEN
3066 CELL_RATIO = ERZ(N)**(ONE/DBLE(NCZ(N)-1))
3067 (K1) = L * (ONE - CELL_RATIO) / (ONE - CELL_RATIO**NCZ(N))
3068
3069 DO K = K1+1,K2
3070 (K) = DZ(K-1) * CELL_RATIO
3071 ENDDO
3072
3073 ELSE
3074 DZ(K1:K2) = L / NCZ(N)
3075 ENDIF
3076
3077 ENDIF
3078
3079 K1 = K2 + 1
3080
3081 ENDDO
3082
3083
3084
3085
3086 = 0
3087 = 0
3088
3089 DO N = 1,NZ
3090
3091 = K1 + NCZ(N) - 1
3092
3093 IF(.NOT.INDEPENDENT_SEGMENT(N)) THEN
3094
3095 L = CPZ(N) - CPZ(N-1)
3096
3097 IF(FIRST_DZ(N)<ZERO) THEN
3098 DZ(K1) = DZ(K1-1)
3099 CALL FIND_CELL_RATIO('FIRST',DZ(K1),L,NCZ(N),CELL_RATIO)
3100 DO K = K1+1,K2
3101 (K) = DZ(K-1) * CELL_RATIO
3102 ENDDO
3103 ELSEIF(LAST_DZ(N)<ZERO) THEN
3104 DZ(K2) = DZ(K2+1)
3105 CALL FIND_CELL_RATIO('LAST ',DZ(K2),L,NCZ(N),CELL_RATIO)
3106 DO K = K2-1,K1,-1
3107 (K) = DZ(K+1) / CELL_RATIO
3108 ENDDO
3109 ENDIF
3110
3111 ENDIF
3112
3113 K1 = K2 + 1
3114
3115 ENDDO
3116
3117
3118
3119
3120 IF(K1>0.AND.K1/=KMAX) THEN
3121 IF(MyPE==0) WRITE(*,*)' ERROR: KMAX MUST BE EQUAL TO THE SUM OF NCZ.'
3122 IF(MyPE==0) WRITE(*,*)' KMAX = ', KMAX
3123 IF(MyPE==0) WRITE(*,*)' SUM OF NCZ = ', K1
3124 call mfix_exit(myPE)
3125 ENDIF
3126
3127
3128
3129 RETURN
3130
3131 END SUBROUTINE GET_DXYZ_FROM_CONTROL_POINTS
3132
3133
3134
3135
3136
3137
3138
3139
3140
3141
3142
3143
3144
3145
3146
3147
3148
3149 SUBROUTINE FIND_CELL_RATIO(POS,D_Target,L,N,ALPHA3)
3150
3151 USE param
3152 USE param1
3153 USE parallel
3154 USE constant
3155 USE run
3156 USE toleranc
3157 USE geometry
3158 USE indices
3159 USE compar
3160 USE sendrecv
3161 USE quadric
3162
3163 IMPLICIT NONE
3164 LOGICAL :: SOLUTION_FOUND
3165
3166 DOUBLE PRECISION :: f1,f2,f3
3167 DOUBLE PRECISION :: ALPHA1,ALPHA2,ALPHA3,D_Target,L,DU
3168 DOUBLE PRECISION, PARAMETER :: ALPHAMAX = 10.0D0
3169 INTEGER :: N,niter
3170 CHARACTER (LEN=5) :: POS
3171
3172 DU = L / DBLE(N)
3173
3174 IF(DU==D_TARGET) THEN
3175 ALPHA3 = 1.0
3176 SOLUTION_FOUND = .TRUE.
3177 RETURN
3178 ELSE
3179
3180 IF(TRIM(POS)=='FIRST') THEN
3181 IF(D_TARGET<DU) THEN
3182 ALPHA1 = ONE
3183 ALPHA2 = ALPHAMAX
3184 ELSE
3185 ALPHA1 = ONE/ALPHAMAX
3186 ALPHA2 = ONE
3187 ENDIF
3188 ELSEIF(TRIM(POS)=='LAST') THEN
3189 IF(D_TARGET>DU) THEN
3190 ALPHA1 = ONE
3191 ALPHA2 = ALPHAMAX
3192 ELSE
3193 ALPHA1 = ONE/ALPHAMAX
3194 ALPHA2 = ONE
3195 ENDIF
3196 ELSE
3197 IF(MyPE==0) WRITE(*,*)' ERROR, IN FUNCTION F: POS MUST BE FIRST OR LAST.'
3198 call mfix_exit(myPE)
3199 ENDIF
3200
3201 ENDIF
3202
3203
3204 f1 = F(POS,ALPHA1,D_Target,L,N)
3205 f2 = F(POS,ALPHA2,D_Target,L,N)
3206
3207
3208
3209
3210
3211
3212
3213 = 1
3214 SOLUTION_FOUND = .FALSE.
3215
3216 if(DABS(f1)<TOL_F) then
3217 = .TRUE.
3218 ALPHA3 = ALPHA1
3219 elseif(DABS(f2)<TOL_F) then
3220 = .TRUE.
3221 ALPHA3 = ALPHA2
3222 elseif(f1*f2 < ZERO) then
3223 = 0
3224 f3 = 2.0d0*TOL_F
3225 do while ( (abs(f3) > TOL_F) .AND. (niter<ITERMAX_INT) )
3226
3227 ALPHA3 = ALPHA1 - f1*(ALPHA2-ALPHA1)/(f2-f1)
3228
3229 = F(POS,ALPHA3,D_Target,L,N)
3230
3231 if(f1*f3<0) then
3232 = ALPHA3
3233 f2 = f3
3234 else
3235 ALPHA1 = ALPHA3
3236 f1 = f3
3237 endif
3238 niter = niter + 1
3239
3240 end do
3241 if (niter < ITERMAX_INT) then
3242 SOLUTION_FOUND = .TRUE.
3243 else
3244 WRITE(*,*) 'Unable to find a solution'
3245 WRITE(*,1000)'between ALPHA1 = ', ALPHA1
3246 WRITE(*,1000)' and ALPHA2 = ', ALPHA2
3247 WRITE(*,1000)'Current value of ALPHA3 = ', ALPHA3
3248 WRITE(*,1000)'Current value of abs(f) = ', DABS(f3)
3249 WRITE(*,1000)'Tolerance = ', TOL_F
3250 WRITE(*,*)'Maximum number of iterations = ', ITERMAX_INT
3251 WRITE(*,*) 'Please increase the intersection tolerance, '
3252 WRITE(*,*) 'or the maximum number of iterations, and try again.'
3253 WRITE(*,*) 'MFiX will exit now.'
3254 CALL MFIX_EXIT(myPE)
3255 SOLUTION_FOUND = .FALSE.
3256 endif
3257 else
3258 WRITE(*,*) 'Unable to find a solution'
3259 WRITE(*,*) 'MFiX will exit now.'
3260 CALL MFIX_EXIT(myPE)
3261 SOLUTION_FOUND = .FALSE.
3262 endif
3263
3264 1000 FORMAT(A,3(2X,G12.5))
3265
3266 RETURN
3267
3268 contains
3269
3270 DOUBLE PRECISION Function F(POS,ALPHAC,D_Target,L,N)
3271 use param1, only: one
3272 USE constant
3273 USE mpi_utility
3274
3275 IMPLICIT NONE
3276 DOUBLE PRECISION:: ALPHAC,D,D_Target,DU,L
3277 INTEGER:: N
3278 CHARACTER (LEN=5) :: POS
3279
3280 DU = L / DBLE(N)
3281
3282 IF(ALPHAC==ONE) THEN
3283 D = DU
3284 ELSE
3285 IF(TRIM(POS)=='FIRST') THEN
3286 D = L * (ONE - ALPHAC) / (ONE -ALPHAC**N)
3287 ELSEIF(TRIM(POS)=='LAST') THEN
3288 D = L * (ONE - ALPHAC) / (ONE -ALPHAC**N) * ALPHAC**(N-1)
3289 ELSE
3290 IF(MyPE==0) WRITE(*,*)' ERROR, IN FUNCTION F: POS MUST BE FIRST OR LAST.'
3291 call mfix_exit(myPE)
3292 ENDIF
3293 ENDIF
3294
3295 F = D - D_Target
3296
3297 RETURN
3298
3299 END FUNCTION F
3300
3301 END SUBROUTINE FIND_CELL_RATIO
3302
3303
3304
3305
3306
3307
3308
3309
3310
3311
3312
3313
3314
3315
3316
3317
3318
3319
3320
3321 SUBROUTINE ADJUST_IJK_SIZE
3322
3323
3324
3325
3326 USE gridmap
3327 USE param
3328 USE param1
3329 USE constant
3330 USE run
3331 USE physprop
3332 USE indices
3333 USE scalars
3334 USE funits
3335 USE leqsol
3336 USE compar
3337 USE mpi_utility
3338 USE bc
3339 USE DISCRETELEMENT
3340
3341 USE cutcell
3342 USE quadric
3343 USE vtk
3344 USE polygon
3345 USE dashboard
3346 USE stl
3347
3348
3349 IMPLICIT NONE
3350
3351
3352
3353
3354
3355
3356
3357
3358
3359 INTEGER :: LC,I,J,K,Q_ID,TOTAL_NUC,IDEAL_NCPP
3360 DOUBLE PRECISION :: X_COPY,Y_COPY,Z_COPY,F_COPY
3361 LOGICAL :: SHIFT,CLIP_FLAG,PRESENT
3362 DOUBLE PRECISION, DIMENSION(0:DIM_I) :: DXT ,XCC
3363 DOUBLE PRECISION, DIMENSION(0:DIM_J) :: DYT ,YCC
3364 DOUBLE PRECISION, DIMENSION(0:DIM_K) :: DZT, ZCC
3365
3366 INTEGER, DIMENSION(0:DIM_I) :: NUC_I,GLOBAL_NUC_I
3367 INTEGER, DIMENSION(0:DIM_J) :: NUC_J,GLOBAL_NUC_J
3368 INTEGER, DIMENSION(0:DIM_K) :: NUC_K,GLOBAL_NUC_K
3369
3370 INTEGER :: IPROC,PSUM
3371
3372 INTEGER, DIMENSION(0:numPEs-1) :: NCPP_OLD,NCPP,NCPP_WITH_GHOST
3373
3374 INTEGER, DIMENSION(0:NODESJ-1) :: JSIZE_OLD
3375
3376 INTEGER :: JSIZE, JREMAIN
3377 INTEGER :: MAXVAL_NCPP_OLD,MINVAL_NCPP_OLD,MAXVAL_NCPP,MINVAL_NCPP
3378 INTEGER :: AVG_NCPP_OLD,AVG_NCPP
3379 DOUBLE PRECISION :: LIP_OLD,LIP,MAXSPEEDUP_OLD,P
3380
3381 INTEGER :: I_OFFSET,J_OFFSET,K_OFFSET,IERR
3382
3383 INTEGER, DIMENSION(0:numPEs-1) :: disp,rcount
3384
3385 LOGICAL :: PRINT_STATISTICS
3386
3387
3388
3389
3390
3391 IF(.NOT.CARTESIAN_GRID) RETURN
3392 IF(.NOT.ADJUST_PROC_DOMAIN_SIZE) RETURN
3393 IF(NODESI*NODESJ*NODESK==1) RETURN
3394
3395
3396 INQUIRE(FILE='gridmap.dat',EXIST=PRESENT)
3397 IF(PRESENT) THEN
3398 IF (myPE == PE_IO) WRITE(*,*)'gridmap was assigned from grimap.dat. Skipping the adjustment.'
3399 RETURN
3400
3401 ENDIF
3402
3403 IF (myPE == PE_IO) CALL PRINT_CG_HEADER
3404
3405 allocate( NCPP_UNIFORM(0:NumPEs-1))
3406
3407 SHORT_GRIDMAP_INIT = .TRUE.
3408 CALL gridmap_init
3409 SHORT_GRIDMAP_INIT = .FALSE.
3410
3411 allocate( ISIZE_ALL(0:NODESI-1))
3412 allocate( JSIZE_ALL(0:NODESJ-1))
3413 allocate( KSIZE_ALL(0:NODESK-1))
3414
3415
3416
3417 NCPP_UNIFORM = ijksize3_all
3418
3419
3420
3421 = IMAX3
3422 DIMENSION_J = JMAX3
3423 DIMENSION_K = KMAX3
3424
3425 PARTIAL_CHECK_03 = .TRUE.
3426
3427 = .FALSE.
3428 CALL CHECK_DATA_CARTESIAN
3429
3430 CALL DEFINE_QUADRICS
3431
3432 SHIFT = .TRUE.
3433
3434 IF(SHIFT) THEN
3435
3436 IF (DO_I) THEN
3437 DXT(IMAX3) = DX(IMAX-1)
3438 DXT(IMAX2) = DX(IMAX-1)
3439 DO LC = IMAX1, IMIN1, -1
3440 DXT(LC) = DX(LC-2)
3441 ENDDO
3442 DXT(IMIN2) = DX(IMIN1)
3443 DXT(IMIN3) =DX(IMIN2)
3444
3445 XCC(IMIN1) = HALF*DXT(IMIN1)
3446 DO I=IMIN1+1,IMAX1
3447 XCC(I) = XCC(I-1) + HALF*(DXT(I-1) + DXT(I))
3448 ENDDO
3449
3450 ENDIF
3451
3452 IF (DO_J) THEN
3453 DYT(JMAX3) = DY(JMAX-1)
3454 DYT(JMAX2) = DY(JMAX-1)
3455 DO LC = JMAX1, JMIN1, -1
3456 DYT(LC) = DY(LC-2)
3457 ENDDO
3458 DYT(JMIN2) = DY(JMIN1)
3459 DYT(JMIN3) =DY(JMIN2)
3460
3461 YCC(JMIN1) = HALF*DYT(JMIN1)
3462 DO J=JMIN1+1,JMAX1
3463 YCC(J) = YCC(J-1) + HALF*(DYT(J-1) + DYT(J))
3464 ENDDO
3465
3466
3467 ENDIF
3468
3469 IF (DO_K) THEN
3470 DZT(KMAX3) = DZ(KMAX-1)
3471 DZT(KMAX2) = DZ(KMAX-1)
3472 DO LC = KMAX1, KMIN1, -1
3473 DZT(LC) = DZ(LC-2)
3474 ENDDO
3475 DZT(KMIN2) = DZ(KMIN1)
3476 DZT(KMIN3) =DZ(KMIN2)
3477
3478 ZCC(KMIN1) = HALF*DZT(KMIN1)
3479 DO K=KMIN1+1,KMAX1
3480 ZCC(K) = ZCC(K-1) + HALF*(DZT(K-1) + DZT(K))
3481 ENDDO
3482
3483 ENDIF
3484
3485 ENDIF
3486
3487
3488
3489 IF(NODESI>1) THEN
3490
3491 IF(myPE == 0.AND.NODESJ*NODESK/=1) THEN
3492 WRITE(*,*)'ERROR IN SUBROUTINE ADJUST_IJK_SIZE.'
3493 WRITE(*,*)'ADJUSTMENT POSSIBLE ONLY FOR DOMAIN DECOMPOSITION In ONE DIRECTION.'
3494 WRITE(*,*)'NODESI,NODESJ,NODESK =',NODESI,NODESJ,NODESK
3495 WRITE(*,*)'MFIX WILL EXIT NOW.'
3496 CALL MFIX_EXIT(myPE)
3497 ENDIF
3498
3499
3500 JSIZE_ALL(0:NODESJ-1) = jmax1-jmin1+1
3501 (0:NODESK-1) = kmax1-kmin1+1
3502
3503
3504 IF (myPE == 0) THEN
3505 WRITE(*,1000)'INFO FROM ADJUST_IJK_SIZE:'
3506 WRITE(*,1000)'ATTEMPTING TO ADJUST DOMAIN SIZE IN I-DIRECTION ...'
3507 WRITE(*,1000)'THIS IS BASED ON AN ESTIMATED (NOT EXACT) NUMBER OF USEFUL CELLS,'
3508 WRITE(*,1000)'AND IT INCLUDES GHOST LAYERS.'
3509
3510 ENDIF
3511
3512 DO I = ISTART1,IEND1
3513
3514 NUC_I(I) = 0
3515
3516 DO J = JSTART1,JEND1
3517 DO K = KSTART1,KEND1
3518 Q_ID = 1
3519 CLIP_FLAG = .FALSE.
3520
3521 X_COPY = XCC(I)
3522 Y_COPY = YCC(J)
3523 Z_COPY = ZCC(K)
3524
3525
3526 CALL EVAL_F('QUADRIC',X_COPY,Y_COPY,Z_COPY,Q_ID,F_COPY,CLIP_FLAG)
3527
3528
3529
3530
3531
3532
3533
3534
3535 IF (F_COPY < TOL_F ) THEN
3536 (I) = NUC_I(I) + 1
3537 ENDIF
3538
3539 ENDDO
3540 ENDDO
3541
3542 ENDDO
3543
3544
3545
3546 CALL allgather_1i (IEND1-ISTART1+1,rcount,IERR)
3547
3548 IF (myPE == 0) THEN
3549 I_OFFSET = 0
3550 ELSE
3551 I_OFFSET = 0
3552 DO iproc=0,myPE-1
3553 I_OFFSET = I_OFFSET + rcount(iproc)
3554 ENDDO
3555 ENDIF
3556
3557 CALL allgather_1i (I_OFFSET,disp,IERR)
3558
3559 call gatherv_1i( NUC_I(ISTART1:IEND1), IEND1-ISTART1+1, GLOBAL_NUC_I(IMIN1:IMAX1), rcount, disp, PE_IO, ierr )
3560
3561
3562
3563 ELSEIF(NODESJ>1) THEN
3564
3565
3566 IF(myPE == 0.AND.NODESI*NODESK/=1) THEN
3567 WRITE(*,*)'ERROR IN SUBROUTINE ADJUST_IJK_SIZE.'
3568 WRITE(*,*)'ADJUSTMENT POSSIBLE ONLY FOR DOMAIN DECOMPOSITION In ONE DIRECTION.'
3569 WRITE(*,*)'NODESI,NODESJ,NODESK =',NODESI,NODESJ,NODESK
3570 WRITE(*,*)'MFIX WILL EXIT NOW.'
3571 CALL MFIX_EXIT(myPE)
3572 ENDIF
3573
3574
3575 ISIZE_ALL(0:NODESI-1) = imax1-imin1+1
3576 (0:NODESK-1) = kmax1-kmin1+1
3577
3578
3579
3580
3581
3582
3583
3584
3585
3586
3587
3588
3589
3590
3591
3592
3593
3594 IF (myPE == 0) THEN
3595 WRITE(*,1000)'INFO FROM ADJUST_IJK_SIZE:'
3596 WRITE(*,1000)'ATTEMPTING TO ADJUST DOMAIN SIZE IN J-DIRECTION ...'
3597 WRITE(*,1000)'THIS IS BASED ON AN ESTIMATED (NOT EXACT) NUMBER OF USEFUL CELLS,'
3598 WRITE(*,1000)'AND IT INCLUDES GHOST LAYERS.'
3599
3600 ENDIF
3601
3602 DO J = JSTART1,JEND1
3603
3604 NUC_J(J) = 0
3605
3606 DO I = ISTART1,IEND1
3607 DO K = KSTART1,KEND1
3608 Q_ID = 1
3609 CLIP_FLAG = .FALSE.
3610
3611 X_COPY = XCC(I)
3612 Y_COPY = YCC(J)
3613 Z_COPY = ZCC(K)
3614
3615
3616 CALL EVAL_F('QUADRIC',X_COPY,Y_COPY,Z_COPY,Q_ID,F_COPY,CLIP_FLAG)
3617
3618
3619
3620
3621
3622
3623
3624
3625 IF (F_COPY < TOL_F ) THEN
3626 (J) = NUC_J(J) + 1
3627 ENDIF
3628
3629 ENDDO
3630 ENDDO
3631
3632 ENDDO
3633
3634
3635
3636 CALL allgather_1i (JEND1-JSTART1+1,rcount,IERR)
3637
3638 IF (myPE == 0) THEN
3639 J_OFFSET = 0
3640 ELSE
3641 J_OFFSET = 0
3642 DO iproc=0,myPE-1
3643 J_OFFSET = J_OFFSET + rcount(iproc)
3644 ENDDO
3645 ENDIF
3646
3647 CALL allgather_1i (J_OFFSET,disp,IERR)
3648
3649 call gatherv_1i( NUC_J(JSTART1:JEND1), JEND1-JSTART1+1, GLOBAL_NUC_J(JMIN1:JMAX1), rcount, disp, PE_IO, ierr )
3650
3651
3652 ELSEIF(NODESK>1) THEN
3653
3654 IF(myPE == 0.AND.NODESI*NODESJ/=1) THEN
3655 WRITE(*,*)'ERROR IN SUBROUTINE ADJUST_IJK_SIZE.'
3656 WRITE(*,*)'ADJUSTMENT POSSIBLE ONLY FOR DOMAIN DECOMPOSITION In ONE DIRECTION.'
3657 WRITE(*,*)'NODESI,NODESJ,NODESK =',NODESI,NODESJ,NODESK
3658 WRITE(*,*)'MFIX WILL EXIT NOW.'
3659 CALL MFIX_EXIT(myPE)
3660 ENDIF
3661
3662
3663 ISIZE_ALL(0:NODESI-1) = imax1-imin1+1
3664 (0:NODESJ-1) = jmax1-jmin1+1
3665
3666
3667 IF (myPE == 0) THEN
3668 WRITE(*,1000)'INFO FROM ADJUST_IJK_SIZE:'
3669 WRITE(*,1000)'ATTEMPTING TO ADJUST DOMAIN SIZE IN K-DIRECTION ...'
3670 WRITE(*,1000)'THIS IS BASED ON AN ESTIMATED (NOT EXACT) NUMBER OF USEFUL CELLS,'
3671 WRITE(*,1000)'AND IT INCLUDES GHOST LAYERS.'
3672 ENDIF
3673
3674 DO K = KSTART1,KEND1
3675
3676 NUC_K(K) = 0
3677
3678 DO I = ISTART1,IEND1
3679 DO J = JSTART1,JEND1
3680 Q_ID = 1
3681 CLIP_FLAG = .FALSE.
3682
3683 X_COPY = XCC(I)
3684 Y_COPY = YCC(J)
3685 Z_COPY = ZCC(K)
3686
3687
3688 CALL EVAL_F('QUADRIC',X_COPY,Y_COPY,Z_COPY,Q_ID,F_COPY,CLIP_FLAG)
3689
3690
3691
3692
3693
3694
3695
3696
3697 IF (F_COPY < TOL_F ) THEN
3698 (K) = NUC_K(K) + 1
3699 ENDIF
3700
3701 ENDDO
3702 ENDDO
3703
3704 ENDDO
3705
3706
3707
3708 CALL allgather_1i (KEND1-KSTART1+1,rcount,IERR)
3709
3710 IF (myPE == 0) THEN
3711 K_OFFSET = 0
3712 ELSE
3713 K_OFFSET = 0
3714 DO iproc=0,myPE-1
3715 K_OFFSET = K_OFFSET + rcount(iproc)
3716 ENDDO
3717 ENDIF
3718
3719 CALL allgather_1i (K_OFFSET,disp,IERR)
3720
3721 call gatherv_1i( NUC_K(KSTART1:KEND1), KEND1-KSTART1+1, GLOBAL_NUC_K(KMIN1:KMAX1), rcount, disp, PE_IO, ierr )
3722
3723
3724
3725 ENDIF
3726
3727
3728
3729
3730
3731 IF (myPE == PE_IO) THEN
3732
3733 IF(NODESI>1) THEN
3734
3735
3736
3737
3738
3739
3740 ELSEIF(NODESJ>1) THEN
3741
3742
3743
3744 = (JMAX1-JMIN1+1)/NODESJ
3745 JSIZE_OLD(0:NODESJ-1) = JSIZE
3746
3747 JREMAIN = (JMAX1-JMIN1+1) - NODESJ*JSIZE
3748 IF (JREMAIN.GE.1) THEN
3749 JSIZE_OLD( 0:(JREMAIN-1) ) = JSIZE + 1
3750 ENDIF
3751
3752
3753
3754
3755
3756
3757
3758
3759
3760 = JSIZE_OLD
3761
3762 CALL MINIMIZE_LOAD_IMBALANCE(NODESJ,GLOBAL_NUC_J(JMIN1:JMAX1),JMIN1,JMAX1,JSIZE_ALL,NCPP,NCPP_WITH_GHOST)
3763
3764
3765 TOTAL_NUC = SUM(NCPP_WITH_GHOST(0:NumPEs-1))
3766 IDEAL_NCPP = TOTAL_NUC / NumPEs
3767
3768 WRITE (*, 1000) 'AFTER OPTIMIZATION:'
3769 WRITE (*, 1010) 'TOTAL NUMBER OF USEFUL CELLS = ',TOTAL_NUC
3770 WRITE (*, 1010) 'IDEAL NUMBER OF CELLS/PROC. = ',IDEAL_NCPP
3771 WRITE (*, 1000) 'ACTUALL CELL COUNT:'
3772 WRITE (*, 1000) '================================================='
3773 WRITE (*, 1000) ' PROCESSOR J-SIZE CELLS/PROC. DIFF. (%)'
3774 WRITE (*, 1000) '================================================='
3775 DO IPROC = 0,numPEs-1
3776 WRITE (*, 1020) IPROC,JSIZE_ALL(IPROC),NCPP_WITH_GHOST(IPROC), &
3777 DBLE(NCPP_WITH_GHOST(IPROC)-IDEAL_NCPP)/DBLE(IDEAL_NCPP)*100.0D0
3778 ENDDO
3779
3780 DO IPROC = 0 ,numPEs-1
3781 NCPP_OLD(IPROC) = (imax3-imin3+1)*(JSIZE_ALL(IPROC)+4)
3782 IF(DO_K) NCPP_OLD(IPROC) = NCPP_OLD(IPROC)*(kmax3-kmin3+1)
3783 ENDDO
3784
3785
3786
3787
3788 = 0
3789 DO IPROC = 0,numPEs-1
3790 PSUM = PSUM + JSIZE_ALL(IPROC)
3791 IF(JSIZE_ALL(IPROC)<5) THEN
3792 WRITE (*, 1010) 'ERROR: J-SIZE TOO SMALL FOR PROCESSOR:',IPROC
3793 WRITE (*, 1010) 'J-SIZE = ',JSIZE_ALL(IPROC)
3794 CALL MFIX_EXIT(myPE)
3795 ENDIF
3796
3797 ENDDO
3798
3799 IF(PSUM/=JMAX) THEN
3800 WRITE (*, 1000) 'ERROR IN ADJUST_IJK_SIZE: UNABLE TO ASSIGN JSIZE TO PROCESSORS.'
3801 WRITE (*, 1000) 'SUM OF JSIZE_ALL DOES NOT MATCH JMAX:'
3802 WRITE (*, 1010) 'SUM OF JSIZE_ALL = ',PSUM
3803 WRITE (*, 1010) 'JMAX1 = ',JMAX
3804 CALL MFIX_EXIT(myPE)
3805 ENDIF
3806
3807 ELSEIF(NODESK>1) THEN
3808
3809 ENDIF
3810
3811 = .TRUE.
3812
3813 IF(PRINT_STATISTICS) THEN
3814
3815
3816
3817
3818
3819
3820
3821
3822
3823
3824
3825
3826
3827 = MAXVAL(NCPP_UNIFORM)
3828 MINVAL_NCPP_OLD = MINVAL(NCPP_UNIFORM)
3829 AVG_NCPP_OLD = SUM(NCPP_UNIFORM)/NUMPES
3830
3831 LIP_OLD = DBLE(MAXVAL_NCPP_OLD-AVG_NCPP_OLD)/DBLE(AVG_NCPP_OLD)*100.0D0
3832
3833 P = DBLE(MAXVAL_NCPP_OLD)/DBLE(AVG_NCPP_OLD)
3834
3835 MAXSPEEDUP_OLD = DBLE(NumPes)*(ONE-LIP_OLD/100.0D0)
3836
3837
3838
3839
3840 MAXVAL_NCPP = MAXVAL(NCPP_WITH_GHOST)
3841 MINVAL_NCPP = MINVAL(NCPP_WITH_GHOST)
3842 AVG_NCPP = SUM(NCPP_WITH_GHOST(0:NumPEs-1))/NUMPES
3843
3844 LIP = DBLE(MAXVAL_NCPP-AVG_NCPP)/DBLE(AVG_NCPP)*100.0D0
3845
3846 P = DBLE(MAXVAL_NCPP)/DBLE(AVG_NCPP)
3847
3848
3849
3850
3851 WRITE (*, 1000) '================================================='
3852 WRITE (*, 1000) 'ESTIMATED PARALLEL LOAD BALANCING STATISTICS'
3853 WRITE (*, 1000) 'COMPARISION BETWEEN UNIFORM SIZE (OLD)'
3854 WRITE (*, 1000) 'AND ADJUSTED SIZE (NEW)'
3855 WRITE (*, 1000) '================================================='
3856 WRITE (*, 1000) ' OLD NEW'
3857 WRITE (*, 1010) 'MAX CELL COUNT : ',MAXVAL_NCPP_OLD,MAXVAL_NCPP
3858 WRITE (*, 1010) 'AT PROCESSOR : ',MAXLOC(NCPP_OLD)-1,MAXLOC(NCPP_WITH_GHOST)-1
3859 WRITE (*, 1010) 'MIN CELL COUNT : ',MINVAL_NCPP_OLD,MINVAL_NCPP
3860 WRITE (*, 1010) 'AT PROCESSOR : ',MINLOC(NCPP_OLD)-1,MINLOC(NCPP_WITH_GHOST)-1
3861 WRITE (*, 1010) 'AVG CELL COUNT : ',AVG_NCPP_OLD,AVG_NCPP
3862 WRITE (*, 1000) ''
3863 WRITE (*, 1030) 'LOAD IMBALANCE (%) : ',LIP_OLD,LIP
3864 WRITE (*, 1000) ''
3865
3866
3867
3868
3869 WRITE (*, 1000) '================================================='
3870
3871 WRITE (*, 1000) 'NOTE: ACTUAL LOAD BALANCING WILL BE COMPUTED AFTER PRE_PROCESSING.'
3872
3873
3874 ENDIF
3875
3876 ENDIF
3877
3878 CONTINUE
3879
3880
3881
3882 CALL BCAST(ISIZE_ALL)
3883 CALL BCAST(JSIZE_ALL)
3884 CALL BCAST(KSIZE_ALL)
3885
3886 DOMAIN_SIZE_ADJUSTED = .TRUE.
3887
3888
3889 1000 FORMAT(1x,A)
3890 1010 FORMAT(1x,A,I10,I10)
3891 1020 FORMAT(1X,I8,2(I12),F12.2)
3892 1030 FORMAT(1X,A,2(F10.1))
3893 1040 FORMAT(F10.1)
3894 1050 FORMAT(1X,3(A))
3895
3896
3897 RETURN
3898
3899 END SUBROUTINE ADJUST_IJK_SIZE
3900
3901
3902
3903
3904
3905
3906
3907
3908
3909
3910
3911
3912
3913
3914
3915
3916
3917
3918
3919
3920 SUBROUTINE REPORT_BEST_PROCESSOR_SIZE
3921
3922
3923
3924
3925 USE gridmap
3926 USE param
3927 USE param1
3928 USE constant
3929 USE run
3930 USE physprop
3931 USE indices
3932 USE scalars
3933 USE funits
3934 USE leqsol
3935 USE compar
3936 USE mpi_utility
3937 USE bc
3938 USE DISCRETELEMENT
3939
3940 USE parallel
3941
3942 USE cutcell
3943 USE quadric
3944 USE vtk
3945 USE polygon
3946 USE dashboard
3947 USE stl
3948
3949
3950 IMPLICIT NONE
3951
3952
3953
3954
3955
3956
3957
3958
3959
3960 =(numPEs==1)
3961
3962 IF(IS_SERIAL) THEN
3963
3964 = NODESI_REPORT
3965 NODESJ = NODESJ_REPORT
3966 NODESK = NODESK_REPORT
3967 numPEs = NODESI*NODESJ*NODESK
3968
3969 IF(numPEs>1.AND.myPE==0) THEN
3970 WRITE(*,1000)'TEMPORARILY SETTING:'
3971 WRITE(*,1010)'NODESI = ',NODESI
3972 WRITE(*,1010)'NODESJ = ',NODESJ
3973 WRITE(*,1010)'NODESK = ',NODESK
3974 WRITE(*,1000)'TO REPORT BEST DOMAIN SIZE FOR PARALLEL RUN'
3975 ENDIF
3976
3977
3978 ENDIF
3979
3980
3981
3982 IF(numPEs>1) CALL REPORT_BEST_IJK_SIZE
3983
3984
3985 IF(IS_SERIAL) THEN
3986
3987
3988
3989 = 1
3990 NODESJ = 1
3991 NODESK = 1
3992 numPEs = 1
3993
3994 ENDIF
3995
3996 1000 FORMAT(1x,A)
3997 1010 FORMAT(1x,A,I10)
3998
3999 RETURN
4000
4001 END SUBROUTINE REPORT_BEST_PROCESSOR_SIZE
4002
4003
4004
4005
4006
4007
4008
4009
4010
4011
4012
4013
4014
4015
4016
4017
4018
4019
4020
4021 SUBROUTINE REPORT_BEST_IJK_SIZE
4022
4023
4024
4025
4026 USE gridmap
4027 USE param
4028 USE param1
4029 USE constant
4030 USE run
4031 USE physprop
4032 USE indices
4033 USE scalars
4034 USE funits
4035 USE leqsol
4036 USE compar
4037 USE mpi_utility
4038 USE bc
4039 USE DISCRETELEMENT
4040
4041 USE parallel
4042
4043 USE cutcell
4044 USE quadric
4045 USE vtk
4046 USE polygon
4047 USE dashboard
4048 USE stl
4049
4050
4051 IMPLICIT NONE
4052
4053
4054
4055
4056
4057
4058
4059
4060
4061 INTEGER :: I,J,K,TOTAL_NUC,IDEAL_NCPP
4062
4063 INTEGER, DIMENSION(0:DIM_I) :: NUC_I
4064 INTEGER, DIMENSION(0:DIM_J) :: NUC_J
4065 INTEGER, DIMENSION(0:DIM_K) :: NUC_K
4066
4067
4068 INTEGER :: ilistsize,jlistsize,klistsize
4069 INTEGER, ALLOCATABLE, DIMENSION(:) :: ALL_NUC_I
4070
4071 INTEGER, ALLOCATABLE, DIMENSION(:) :: ALL_LIST_I
4072 INTEGER, ALLOCATABLE, DIMENSION(:) :: GLOBAL_NUC_I
4073
4074 INTEGER, ALLOCATABLE, DIMENSION(:) :: ALL_NUC_J
4075
4076 INTEGER, ALLOCATABLE, DIMENSION(:) :: ALL_LIST_J
4077 INTEGER, ALLOCATABLE, DIMENSION(:) :: GLOBAL_NUC_J
4078
4079 INTEGER, ALLOCATABLE, DIMENSION(:) :: ALL_NUC_K
4080
4081 INTEGER, ALLOCATABLE, DIMENSION(:) :: ALL_LIST_K
4082 INTEGER, ALLOCATABLE, DIMENSION(:) :: GLOBAL_NUC_K
4083
4084
4085 INTEGER :: IPROC,PSUM
4086
4087 INTEGER, DIMENSION(0:numPEs-1) :: NCPP_OLD,NCPP,NCPP_OLD_WITH_GHOST,NCPP_WITH_GHOST
4088
4089 INTEGER, DIMENSION(0:NODESI-1) :: ISIZE_OLD
4090 INTEGER, DIMENSION(0:NODESJ-1) :: JSIZE_OLD
4091 INTEGER, DIMENSION(0:NODESK-1) :: KSIZE_OLD
4092
4093 INTEGER :: JSIZE, IREMAIN,ISIZE, JREMAIN,KSIZE, KREMAIN
4094 DOUBLE PRECISION :: LIP_OLD
4095
4096 INTEGER :: I_OFFSET,J_OFFSET,K_OFFSET,IERR
4097
4098 INTEGER, DIMENSION(0:numPEs-1) :: disp,rcount
4099
4100 INTEGER :: IPROC_OF_MAX_OLD,IPROC_OF_MIN_OLD
4101
4102
4103
4104 IF(.NOT.REPORT_BEST_DOMAIN_SIZE) RETURN
4105
4106 IF(.NOT.CARTESIAN_GRID) RETURN
4107
4108 IF(NODESI*NODESJ*NODESK==1) RETURN
4109
4110 IF(.not.allocated(ISIZE_ALL)) allocate( ISIZE_ALL(0:NODESI-1))
4111 IF(.not.allocated(JSIZE_ALL)) allocate( JSIZE_ALL(0:NODESJ-1))
4112 IF(.not.allocated(KSIZE_ALL)) allocate( KSIZE_ALL(0:NODESK-1))
4113
4114 ISIZE_ALL(0:NODESI-1) = imax1-imin1+1
4115 (0:NODESJ-1) = jmax1-jmin1+1
4116 KSIZE_ALL(0:NODESK-1) = kmax1-kmin1+1
4117
4118
4119 IF(NODESI>1) THEN
4120
4121
4122
4123 DO I = ISTART1,IEND1
4124 NUC_I(I) = 0
4125 DO J = JSTART3,JEND3
4126 DO K = KSTART3,KEND3
4127 IF( .NOT.DEAD_CELL_AT(I,J,K)) THEN
4128 (I) = NUC_I(I) + 1
4129 ENDIF
4130 ENDDO
4131 ENDDO
4132 ENDDO
4133
4134
4135
4136
4137
4138 CALL allgather_1i (IEND1-ISTART1+1,rcount,IERR)
4139
4140 IF (myPE == 0) THEN
4141 I_OFFSET = 0
4142 ELSE
4143 I_OFFSET = 0
4144 DO iproc=0,myPE-1
4145 I_OFFSET = I_OFFSET + rcount(iproc)
4146 ENDDO
4147 ENDIF
4148
4149 CALL allgather_1i (I_OFFSET,disp,IERR)
4150
4151 ilistsize=SUM(rcount)
4152
4153 allocate( ALL_NUC_I(ilistsize))
4154 allocate( ALL_LIST_I(ilistsize))
4155 allocate( GLOBAL_NUC_I(IMIN1:IMAX1))
4156
4157
4158 call gatherv_1i( (/(I,I=ISTART1,IEND1)/), IEND1-ISTART1+1, ALL_LIST_I(:), rcount, disp, PE_IO, ierr )
4159 call gatherv_1i( NUC_I(ISTART1:IEND1), IEND1-ISTART1+1, ALL_NUC_I(:), rcount, disp, PE_IO, ierr )
4160
4161
4162 IF (myPE == 0) THEN
4163 GLOBAL_NUC_I = 0
4164 DO I=1,ilistsize
4165 GLOBAL_NUC_I(ALL_LIST_I(I)) = GLOBAL_NUC_I(ALL_LIST_I(I)) + ALL_NUC_I(I)
4166 ENDDO
4167 ENDIF
4168
4169 ENDIF
4170
4171 IF(NODESJ>1) THEN
4172
4173
4174
4175 DO J = JSTART1,JEND1
4176 NUC_J(J) = 0
4177 DO I = ISTART3,IEND3
4178 DO K = KSTART3,KEND3
4179 IF( .NOT.DEAD_CELL_AT(I,J,K)) THEN
4180 (J) = NUC_J(J) + 1
4181 ENDIF
4182 ENDDO
4183 ENDDO
4184 ENDDO
4185
4186
4187
4188
4189
4190 CALL allgather_1i (JEND1-JSTART1+1,rcount,IERR)
4191
4192 IF (myPE == 0) THEN
4193 J_OFFSET = 0
4194 ELSE
4195 J_OFFSET = 0
4196 DO iproc=0,myPE-1
4197 J_OFFSET = J_OFFSET + rcount(iproc)
4198 ENDDO
4199 ENDIF
4200
4201 CALL allgather_1i (J_OFFSET,disp,IERR)
4202
4203 Jlistsize=SUM(rcount)
4204
4205 allocate( ALL_NUC_J(Jlistsize))
4206 allocate( ALL_LIST_J(Jlistsize))
4207 allocate( GLOBAL_NUC_J(JMIN1:JMAX1))
4208
4209
4210 call gatherv_1i( (/(J,J=JSTART1,JEND1)/), JEND1-JSTART1+1, ALL_LIST_J(:), rcount, disp, PE_IO, ierr )
4211 call gatherv_1i( NUC_J(JSTART1:JEND1), JEND1-JSTART1+1, ALL_NUC_J(:), rcount, disp, PE_IO, ierr )
4212
4213
4214 IF (myPE == 0) THEN
4215 GLOBAL_NUC_J = 0
4216 DO J=1,Jlistsize
4217 GLOBAL_NUC_J(ALL_LIST_J(J)) = GLOBAL_NUC_J(ALL_LIST_J(J)) + ALL_NUC_J(J)
4218 ENDDO
4219 ENDIF
4220
4221 ENDIF
4222
4223 IF(NODESK>1) THEN
4224
4225
4226
4227 DO K = KSTART1,KEND1
4228 NUC_K(K) = 0
4229 DO I = ISTART3,IEND3
4230 DO J = JSTART3,JEND3
4231 IF( .NOT.DEAD_CELL_AT(I,J,K)) THEN
4232 (K) = NUC_K(K) + 1
4233 ENDIF
4234 ENDDO
4235 ENDDO
4236 ENDDO
4237
4238
4239
4240
4241
4242 CALL allgather_1i (KEND1-KSTART1+1,rcount,IERR)
4243
4244 IF (myPE == 0) THEN
4245 K_OFFSET = 0
4246 ELSE
4247 K_OFFSET = 0
4248 DO iproc=0,myPE-1
4249 K_OFFSET = K_OFFSET + rcount(iproc)
4250 ENDDO
4251 ENDIF
4252
4253 CALL allgather_1i (K_OFFSET,disp,IERR)
4254
4255 Klistsize=SUM(rcount)
4256
4257 allocate( ALL_NUC_K(Klistsize))
4258 allocate( ALL_LIST_K(Klistsize))
4259 allocate( GLOBAL_NUC_K(KMIN1:KMAX1))
4260
4261
4262 call gatherv_1i( (/(K,K=KSTART1,KEND1)/), KEND1-KSTART1+1, ALL_LIST_K(:), rcount, disp, PE_IO, ierr )
4263 call gatherv_1i( NUC_K(KSTART1:KEND1), KEND1-KSTART1+1, ALL_NUC_K(:), rcount, disp, PE_IO, ierr )
4264
4265
4266 IF (myPE == 0) THEN
4267 GLOBAL_NUC_K = 0
4268 DO K=1,Klistsize
4269 GLOBAL_NUC_K(ALL_LIST_K(K)) = GLOBAL_NUC_K(ALL_LIST_K(K)) + ALL_NUC_K(K)
4270 ENDDO
4271 ENDIF
4272
4273 ENDIF
4274
4275
4276
4277
4278
4279 IF (myPE == PE_IO) THEN
4280
4281 IF(NODESI>1) THEN
4282
4283
4284
4285 = (IMAX1-IMIN1+1)/NODESI
4286 ISIZE_OLD(0:NODESI-1) = ISIZE
4287
4288 IREMAIN = (IMAX1-IMIN1+1) - NODESI*ISIZE
4289 IF (IREMAIN.GE.1) THEN
4290 ISIZE_OLD( 0:(IREMAIN-1) ) = ISIZE + 1
4291 ENDIF
4292
4293 ISIZE_ALL = ISIZE_OLD
4294
4295
4296 CALL GET_LIP_WITH_GHOST_LAYERS(NODESI,GLOBAL_NUC_I(IMIN1:IMAX1),IMIN1, &
4297 IMAX1,ISIZE_ALL,NCPP_OLD,NCPP_OLD_WITH_GHOST, &
4298 LIP_OLD,IPROC_OF_MAX_OLD,IPROC_OF_MIN_OLD)
4299
4300 TOTAL_NUC = SUM(NCPP_OLD_WITH_GHOST(0:NODESI-1))
4301 IDEAL_NCPP = TOTAL_NUC / NODESI
4302 WRITE(*,1000)'INFO FROM REPORT_BEST_IJK_SIZE:'
4303 WRITE(*,1000)'ATTEMPTING TO REPORT BEST DOMAIN SIZE IN I-DIRECTION ...'
4304 WRITE (*, 1010) 'TOTAL NUMBER OF USEFUL CELLS = ',TOTAL_NUC
4305 WRITE (*, 1010) 'IDEAL NUMBER OF CELLS/I-NODE = ',IDEAL_NCPP
4306 WRITE (*, 1000) 'BEFORE OPTIMIZATION:'
4307 WRITE (*, 1000) '================================================='
4308 WRITE (*, 1000) ' I-NODE I-SIZE CELLS/NODE DIFF. (%)'
4309 WRITE (*, 1000) '================================================='
4310 DO IPROC = 0,NODESI-1
4311 WRITE (*, 1020) IPROC,ISIZE_ALL(IPROC),NCPP_OLD_WITH_GHOST(IPROC), &
4312 DBLE(NCPP_OLD_WITH_GHOST(IPROC)-IDEAL_NCPP)/DBLE(IDEAL_NCPP)*100.0D0
4313 ENDDO
4314 WRITE (*, 1000) '================================================='
4315
4316
4317 CALL MINIMIZE_LOAD_IMBALANCE(NODESI,GLOBAL_NUC_I(IMIN1:IMAX1),IMIN1,IMAX1,ISIZE_ALL,NCPP,NCPP_WITH_GHOST)
4318
4319 TOTAL_NUC = SUM(NCPP_WITH_GHOST(0:NODESI-1))
4320 IDEAL_NCPP = TOTAL_NUC / NODESI
4321 WRITE (*, 1010) 'TOTAL NUMBER OF USEFUL CELLS = ',TOTAL_NUC
4322 WRITE (*, 1010) 'IDEAL NUMBER OF CELLS/I-NODE = ',IDEAL_NCPP
4323 WRITE (*, 1010) 'SINCE GHOST CELLS ARE INCLUDED, THE TOTALS'
4324 WRITE (*, 1010) 'BEFORE AND AFTER OPTIMIZATION MAY NOT MATCH.'
4325 WRITE (*, 1000) 'AFTER OPTIMIZATION:'
4326 WRITE (*, 1000) '================================================='
4327 WRITE (*, 1000) ' I-NODE I-SIZE CELLS/NODE DIFF. (%)'
4328 WRITE (*, 1000) '================================================='
4329 DO IPROC = 0,NODESI-1
4330 WRITE (*, 1020) IPROC,ISIZE_ALL(IPROC),NCPP_WITH_GHOST(IPROC), &
4331 DBLE(NCPP_WITH_GHOST(IPROC)-IDEAL_NCPP)/DBLE(IDEAL_NCPP)*100.0D0
4332 ENDDO
4333 WRITE (*, 1000) '================================================='
4334
4335
4336 = 0
4337 DO IPROC = 0,NODESI-1
4338 PSUM = PSUM + ISIZE_ALL(IPROC)
4339 IF(ISIZE_ALL(IPROC)<5) THEN
4340 WRITE (*, 1010) 'ERROR: I-SIZE TOO SMALL FOR I-NODE:',IPROC
4341 WRITE (*, 1010) 'I-SIZE = ',ISIZE_ALL(IPROC)
4342 CALL MFIX_EXIT(myPE)
4343 ENDIF
4344 ENDDO
4345
4346 IF(PSUM/=IMAX) THEN
4347 WRITE (*, 1000) 'ERROR IN ADJUST_IJK_SIZE: UNABLE TO ASSIGN ISIZE TO PROCESSORS.'
4348 WRITE (*, 1000) 'SUM OF ISIZE_ALL DOES NOT MATCH IMAX:'
4349 WRITE (*, 1010) 'SUM OF ISIZE_ALL = ',PSUM
4350 WRITE (*, 1010) 'IMAX = ',IMAX
4351 CALL MFIX_EXIT(myPE)
4352 ENDIF
4353
4354 ENDIF
4355
4356
4357 IF(NODESJ>1) THEN
4358
4359
4360 = (JMAX1-JMIN1+1)/NODESJ
4361 JSIZE_OLD(0:NODESJ-1) = JSIZE
4362
4363 JREMAIN = (JMAX1-JMIN1+1) - NODESJ*JSIZE
4364 IF (JREMAIN.GE.1) THEN
4365 JSIZE_OLD( 0:(JREMAIN-1) ) = JSIZE + 1
4366 ENDIF
4367
4368 JSIZE_ALL = JSIZE_OLD
4369
4370
4371 CALL GET_LIP_WITH_GHOST_LAYERS(NODESJ,GLOBAL_NUC_J(JMIN1:JMAX1),JMIN1,JMAX1,JSIZE_ALL, &
4372 NCPP_OLD,NCPP_OLD_WITH_GHOST,LIP_OLD,IPROC_OF_MAX_OLD,IPROC_OF_MIN_OLD)
4373
4374 TOTAL_NUC = SUM(NCPP_OLD_WITH_GHOST(0:NODESJ-1))
4375 IDEAL_NCPP = TOTAL_NUC / NODESJ
4376 WRITE(*,1000)'INFO FROM REPORT_BEST_IJK_SIZE:'
4377 WRITE(*,1000)'ATTEMPTING TO REPORT BEST DOMAIN SIZE IN J-DIRECTION ...'
4378 WRITE (*, 1010) 'TOTAL NUMBER OF USEFUL CELLS = ',TOTAL_NUC
4379 WRITE (*, 1010) 'IDEAL NUMBER OF CELLS/J-NODE = ',IDEAL_NCPP
4380 WRITE (*, 1000) 'BEFORE OPTIMIZATION:'
4381 WRITE (*, 1000) '================================================='
4382 WRITE (*, 1000) ' J-NODE J-SIZE CELLS/NODE DIFF. (%)'
4383 WRITE (*, 1000) '================================================='
4384 DO IPROC = 0,NODESJ-1
4385 WRITE (*, 1020) IPROC,JSIZE_ALL(IPROC),NCPP_OLD_WITH_GHOST(IPROC), &
4386 DBLE(NCPP_OLD_WITH_GHOST(IPROC)-IDEAL_NCPP)/DBLE(IDEAL_NCPP)*100.0D0
4387 ENDDO
4388 WRITE (*, 1000) '================================================='
4389
4390
4391 CALL MINIMIZE_LOAD_IMBALANCE(NODESJ,GLOBAL_NUC_J(JMIN1:JMAX1),JMIN1,JMAX1,JSIZE_ALL,NCPP,NCPP_WITH_GHOST)
4392
4393 TOTAL_NUC = SUM(NCPP_WITH_GHOST(0:NODESJ-1))
4394 IDEAL_NCPP = TOTAL_NUC / NODESJ
4395 WRITE (*, 1010) 'TOTAL NUMBER OF USEFUL CELLS = ',TOTAL_NUC
4396 WRITE (*, 1010) 'IDEAL NUMBER OF CELLS/J-NODE = ',IDEAL_NCPP
4397 WRITE (*, 1010) 'SINCE GHOST CELLS ARE INCLUDED, THE TOTALS'
4398 WRITE (*, 1010) 'BEFORE AND AFTER OPTIMIZATION MAY NOT MATCH.'
4399 WRITE (*, 1000) 'AFTER OPTIMIZATION:'
4400 WRITE (*, 1000) '================================================='
4401 WRITE (*, 1000) ' J-NODE J-SIZE CELLS/NODE DIFF. (%)'
4402 WRITE (*, 1000) '================================================='
4403 DO IPROC = 0,NODESJ-1
4404 WRITE (*, 1020) IPROC,JSIZE_ALL(IPROC),NCPP_WITH_GHOST(IPROC), &
4405 DBLE(NCPP_WITH_GHOST(IPROC)-IDEAL_NCPP)/DBLE(IDEAL_NCPP)*100.0D0
4406 ENDDO
4407 WRITE (*, 1000) '================================================='
4408
4409
4410 = 0
4411 DO IPROC = 0,NODESJ-1
4412 PSUM = PSUM + JSIZE_ALL(IPROC)
4413 IF(JSIZE_ALL(IPROC)<5) THEN
4414 WRITE (*, 1010) 'ERROR: J-SIZE TOO SMALL FOR J-NODE:',IPROC
4415 WRITE (*, 1010) 'J-SIZE = ',JSIZE_ALL(IPROC)
4416 CALL MFIX_EXIT(myPE)
4417 ENDIF
4418 ENDDO
4419
4420 IF(PSUM/=JMAX) THEN
4421 WRITE (*, 1000) 'ERROR IN ADJUST_IJK_SIZE: UNABLE TO ASSIGN JSIZE TO PROCESSORS.'
4422 WRITE (*, 1000) 'SUM OF JSIZE_ALL DOES NOT MATCH JMAX:'
4423 WRITE (*, 1010) 'SUM OF JSIZE_ALL = ',PSUM
4424 WRITE (*, 1010) 'JMAX = ',JMAX
4425 CALL MFIX_EXIT(myPE)
4426 ENDIF
4427
4428 ENDIF
4429
4430 IF(NODESK>1) THEN
4431
4432
4433 = (KMAX1-KMIN1+1)/NODESK
4434 KSIZE_OLD(0:NODESK-1) = KSIZE
4435
4436 KREMAIN = (KMAX1-KMIN1+1) - NODESK*KSIZE
4437 IF (KREMAIN.GE.1) THEN
4438 KSIZE_OLD( 0:(KREMAIN-1) ) = KSIZE + 1
4439 ENDIF
4440
4441 KSIZE_ALL = KSIZE_OLD
4442
4443
4444 CALL GET_LIP_WITH_GHOST_LAYERS(NODESK,GLOBAL_NUC_K(KMIN1:KMAX1),KMIN1,KMAX1,KSIZE_ALL, &
4445 NCPP_OLD,NCPP_OLD_WITH_GHOST,LIP_OLD,IPROC_OF_MAX_OLD,IPROC_OF_MIN_OLD)
4446
4447 TOTAL_NUC = SUM(NCPP_OLD_WITH_GHOST(0:NODESK-1))
4448 IDEAL_NCPP = TOTAL_NUC / NODESK
4449 WRITE(*,1000)'INFO FROM REPORT_BEST_IJK_SIZE:'
4450 WRITE(*,1000)'ATTEMPTING TO REPORT BEST DOMAIN SIZE IN K-DIRECTION ...'
4451 WRITE (*, 1010) 'TOTAL NUMBER OF USEFUL CELLS = ',TOTAL_NUC
4452 WRITE (*, 1010) 'IDEAL NUMBER OF CELLS/K_NODE = ',IDEAL_NCPP
4453 WRITE (*, 1000) 'BEFORE OPTIMIZATION:'
4454 WRITE (*, 1000) '================================================='
4455 WRITE (*, 1000) ' K-NODE K-SIZE CELLS/NODE DIFF. (%)'
4456 WRITE (*, 1000) '================================================='
4457 DO IPROC = 0,NODESK-1
4458 WRITE (*, 1020) IPROC,KSIZE_ALL(IPROC),NCPP_OLD_WITH_GHOST(IPROC), &
4459 DBLE(NCPP_OLD_WITH_GHOST(IPROC)-IDEAL_NCPP)/DBLE(IDEAL_NCPP)*100.0D0
4460 ENDDO
4461 WRITE (*, 1000) '================================================='
4462
4463
4464 CALL MINIMIZE_LOAD_IMBALANCE(NODESK,GLOBAL_NUC_K(KMIN1:KMAX1),KMIN1,KMAX1,KSIZE_ALL,NCPP,NCPP_WITH_GHOST)
4465
4466 TOTAL_NUC = SUM(NCPP_WITH_GHOST(0:NODESK-1))
4467 IDEAL_NCPP = TOTAL_NUC / NODESK
4468 WRITE (*, 1010) 'TOTAL NUMBER OF USEFUL CELLS = ',TOTAL_NUC
4469 WRITE (*, 1010) 'IDEAL NUMBER OF CELLS/K_NODE = ',IDEAL_NCPP
4470 WRITE (*, 1010) 'SINCE GHOST CELLS ARE INCLUDED, THE TOTALS'
4471 WRITE (*, 1010) 'BEFORE AND AFTER OPTIMIZATION MAY NOT MATCH.'
4472 WRITE (*, 1000) 'AFTER OPTIMIZATION:'
4473 WRITE (*, 1000) '================================================='
4474 WRITE (*, 1000) ' K-NODE K-SIZE CELLS/NODE DIFF. (%)'
4475 WRITE (*, 1000) '================================================='
4476 DO IPROC = 0,NODESK-1
4477 WRITE (*, 1020) IPROC,KSIZE_ALL(IPROC),NCPP_WITH_GHOST(IPROC), &
4478 DBLE(NCPP_WITH_GHOST(IPROC)-IDEAL_NCPP)/DBLE(IDEAL_NCPP)*100.0D0
4479 ENDDO
4480 WRITE (*, 1000) '================================================='
4481
4482
4483 = 0
4484 DO IPROC = 0,NODESK-1
4485 PSUM = PSUM + KSIZE_ALL(IPROC)
4486 IF(KSIZE_ALL(IPROC)<5) THEN
4487 WRITE (*, 1010) 'ERROR: K-SIZE TOO SMALL FOR K-NODE:',IPROC
4488 WRITE (*, 1010) 'K-SIZE = ',KSIZE_ALL(IPROC)
4489 CALL MFIX_EXIT(myPE)
4490 ENDIF
4491 ENDDO
4492
4493 IF(PSUM/=KMAX) THEN
4494 WRITE (*, 1000) 'ERROR IN ADJUST_IJK_SIZE: UNABLE TO ASSIGN KSIZE TO PROCESSORS.'
4495 WRITE (*, 1000) 'SUM OF KSIZE_ALL DOES NOT MATCH KMAX:'
4496 WRITE (*, 1010) 'SUM OF KSIZE_ALL = ',PSUM
4497 WRITE (*, 1010) 'KMAX = ',KMAX
4498 CALL MFIX_EXIT(myPE)
4499 ENDIF
4500
4501
4502 ENDIF
4503
4504
4505 OPEN(CONVERT='BIG_ENDIAN',UNIT=777, FILE='suggested_gridmap.dat')
4506 WRITE (777, 1005) NODESI,NODESJ,NODESK, ' ! NODESI, NODESJ, NODESK'
4507 DO IPROC = 0,NODESI-1
4508 WRITE(777,1060) IPROC,Isize_all(IPROC)
4509 ENDDO
4510 DO IPROC = 0,NODESJ-1
4511 WRITE(777,1060) IPROC,Jsize_all(IPROC)
4512 ENDDO
4513 DO IPROC = 0,NODESK-1
4514 WRITE(777,1060) IPROC,Ksize_all(IPROC)
4515 ENDDO
4516
4517 CLOSE(777)
4518 WRITE (*, 1000) '================================================='
4519 WRITE (*, 1000) 'GRID PARTITION SAVED IN FILE: suggested_gridmap.dat'
4520 WRITE (*, 1000) 'TO USE THIS DISTRIBUTION, RENAME THE FILE AS: gridmap.dat'
4521 WRITE (*, 1000) 'AND RUN MFIX AGAIN.'
4522
4523 WRITE (*, 1000) '================================================='
4524
4525
4526 ENDIF
4527
4528
4529 call parallel_fin
4530
4531 STOP
4532
4533
4534
4535
4536
4537
4538
4539
4540
4541 FORMAT(1x,A)
4542 1005 FORMAT(1x,I10,I10,I10,A)
4543 1010 FORMAT(1x,A,I10,I10)
4544 1020 FORMAT(1X,I8,2(I12),F12.2)
4545 1030 FORMAT(1X,A,2(F10.1))
4546 1040 FORMAT(F10.1)
4547 1050 FORMAT(1X,3(A))
4548 1060 FORMAT(1x,I10,I10)
4549
4550 RETURN
4551
4552 END SUBROUTINE REPORT_BEST_IJK_SIZE
4553
4554
4555
4556
4557
4558
4559
4560
4561
4562
4563
4564
4565
4566
4567
4568
4569
4570
4571 SUBROUTINE GET_LIP_WITH_GHOST_LAYERS(NODESL,NUC_L,LMIN1,LMAX1,L_SIZE,NCPP,NCPP_WITH_GHOST,LIP,IPROC_OF_MAX,IPROC_OF_MIN)
4572
4573
4574
4575
4576 USE gridmap
4577 USE param
4578 USE param1
4579 USE constant
4580 USE run
4581 USE physprop
4582 USE indices
4583 USE scalars
4584 USE funits
4585 USE leqsol
4586 USE compar
4587 USE mpi_utility
4588 USE bc
4589 USE DISCRETELEMENT
4590
4591 USE cutcell
4592 USE quadric
4593 USE vtk
4594 USE polygon
4595 USE dashboard
4596 USE stl
4597
4598
4599 IMPLICIT NONE
4600
4601
4602
4603
4604
4605
4606
4607
4608
4609 INTEGER :: NODESL,L,LMIN1,LMAX1,TOTAL_NUC,TOTAL_NUC_WITH_GHOST,IPROC_OF_MAX,IPROC_OF_MIN
4610
4611 INTEGER :: LCOUNT1,LCOUNT2,MINVAL_NCPP,MAXVAL_NCPP,IDEAL_NCPP
4612 INTEGER, DIMENSION(LMIN1:LMAX1) :: NUC_L
4613
4614 INTEGER :: IPROC
4615
4616 INTEGER, DIMENSION(0:NODESL-1) :: NCPP,NCPP_WITH_GHOST,L_SIZE,L1,L2
4617
4618
4619 DOUBLE PRECISION :: LIP
4620
4621
4622
4623 = LMAX1 - LMIN1 + 1
4624 LCOUNT2 = SUM(L_SIZE(0:NODESL-1))
4625
4626 IF(LCOUNT1/=LCOUNT2) THEN
4627 WRITE(*,*)' ERROR: SUM OF CELLS DO NOT MATCH:',LCOUNT1,LCOUNT2
4628 CALL MFIX_EXIT(myPE)
4629 ENDIF
4630
4631 L1(0) = LMIN1
4632 L2(0) = L1(0) + L_SIZE(0) - 1
4633
4634 DO IPROC = 1,NODESL-1
4635 L1(IPROC) = L2(IPROC-1) + 1
4636 L2(IPROC) = L1(IPROC) + L_SIZE(IPROC) - 1
4637 ENDDO
4638
4639 DO IPROC = 0,NODESL-1
4640 NCPP(IPROC) = SUM(NUC_L(L1(IPROC):L2(IPROC)))
4641
4642
4643 ENDDO
4644
4645 TOTAL_NUC = 0
4646
4647 DO L=LMIN1,LMAX1
4648 TOTAL_NUC = TOTAL_NUC + NUC_L(L)
4649 ENDDO
4650
4651 NCPP_WITH_GHOST(0) = NCPP(0) + 2*NUC_L(L1(0)) + NUC_L(L1(1)) + NUC_L(L1(1)+1)
4652
4653 DO IPROC = 1,NODESL-2
4654 NCPP_WITH_GHOST(IPROC) = NCPP(IPROC) &
4655 + NUC_L(L2(IPROC-1)) + NUC_L(L2(IPROC-1)-1) &
4656 + NUC_L(L1(IPROC+1)) + NUC_L(L1(IPROC+1)+1)
4657 ENDDO
4658
4659 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)
4660
4661 TOTAL_NUC_WITH_GHOST = 0
4662 DO IPROC = 0,NODESL-1
4663
4664 = TOTAL_NUC_WITH_GHOST + NCPP_WITH_GHOST(IPROC)
4665 ENDDO
4666
4667
4668 IDEAL_NCPP = TOTAL_NUC_WITH_GHOST / NumPEs
4669
4670 MAXVAL_NCPP = MAXVAL(NCPP_WITH_GHOST)
4671 MINVAL_NCPP = MINVAL(NCPP_WITH_GHOST)
4672
4673 LIP = DBLE(MAXVAL_NCPP-IDEAL_NCPP)/DBLE(IDEAL_NCPP)*100.0D0
4674
4675
4676
4677 = MAXLOC(NCPP_WITH_GHOST,1)-1
4678 IPROC_OF_MIN = MINLOC(NCPP_WITH_GHOST,1)-1
4679
4680
4681
4682
4683
4684 RETURN
4685 END SUBROUTINE GET_LIP_WITH_GHOST_LAYERS
4686
4687
4688
4689
4690
4691
4692
4693
4694
4695
4696
4697
4698
4699
4700
4701
4702
4703
4704
4705 SUBROUTINE MINIMIZE_LOAD_IMBALANCE(NODESL,NUC_L,LMIN1,LMAX1,L_SIZE,NCPP,NCPP_WITH_GHOST)
4706
4707
4708
4709
4710 USE gridmap
4711 USE param
4712 USE param1
4713 USE constant
4714 USE run
4715 USE physprop
4716 USE indices
4717 USE scalars
4718 USE funits
4719 USE leqsol
4720 USE compar
4721 USE mpi_utility
4722 USE bc
4723 USE DISCRETELEMENT
4724
4725 USE cutcell
4726 USE quadric
4727 USE vtk
4728 USE polygon
4729 USE dashboard
4730 USE stl
4731
4732
4733 IMPLICIT NONE
4734
4735
4736
4737
4738
4739
4740
4741
4742
4743 INTEGER :: NODESL,LMIN1,LMAX1,IPROC_OF_MAX,IPROC_OF_MIN
4744
4745 INTEGER, DIMENSION(LMIN1:LMAX1) :: NUC_L
4746
4747 INTEGER :: N,NOIMPROVEMENT
4748
4749 INTEGER,PARAMETER :: NAMAX=10000
4750
4751 INTEGER, DIMENSION(0:numPEs-1) :: NCPP,NCPP_WITH_GHOST,L_SIZE,BEST_L_SIZE,BEST_NCPP,BEST_NCPP_WITH_GHOST
4752
4753 DOUBLE PRECISION :: LIP,BEST_LIP
4754
4755
4756
4757
4758
4759
4760
4761
4762 CALL GET_LIP_WITH_GHOST_LAYERS(NODESL,NUC_L,LMIN1,LMAX1,L_SIZE,NCPP,NCPP_WITH_GHOST,LIP,IPROC_OF_MAX,IPROC_OF_MIN)
4763
4764 BEST_LIP = LIP
4765 BEST_L_SIZE = L_SIZE
4766
4767
4768
4769
4770
4771
4772
4773
4774
4775
4776
4777
4778
4779
4780
4781 print*,'ATTEMPTING TO OPTIMIZE LOAD BALANCE...'
4782
4783 NOIMPROVEMENT=0
4784
4785 DO N = 1,NAMAX
4786
4787 L_SIZE(IPROC_OF_MAX) = L_SIZE(IPROC_OF_MAX) - 1
4788 L_SIZE(IPROC_OF_MIN) = L_SIZE(IPROC_OF_MIN) + 1
4789
4790 CALL GET_LIP_WITH_GHOST_LAYERS(NODESL,NUC_L,LMIN1,LMAX1,L_SIZE,NCPP,NCPP_WITH_GHOST,LIP,IPROC_OF_MAX,IPROC_OF_MIN)
4791
4792
4793
4794
4795
4796 IF(LIP<BEST_LIP) THEN
4797 BEST_LIP = LIP
4798 BEST_L_SIZE = L_SIZE
4799 BEST_NCPP = NCPP
4800 BEST_NCPP_WITH_GHOST = NCPP_WITH_GHOST
4801 NOIMPROVEMENT=0
4802 ELSE
4803 NOIMPROVEMENT = NOIMPROVEMENT + 1
4804 ENDIF
4805
4806
4807
4808
4809
4810
4811
4812
4813
4814
4815 IF(NOIMPROVEMENT==10) THEN
4816 WRITE (*, 1000) 'OPTIMIZED LOAD BALANCE REACHED.'
4817 EXIT
4818 ENDIF
4819
4820 ENDDO
4821
4822
4823 = BEST_L_SIZE
4824 NCPP = BEST_NCPP
4825 NCPP_WITH_GHOST = BEST_NCPP_WITH_GHOST
4826
4827
4828
4829
4830
4831
4832
4833
4834
4835
4836 FORMAT(1x,A)
4837 1010 FORMAT(1x,A,I10,I10)
4838 1020 FORMAT(1X,I8,2(I12),F12.2)
4839
4840 RETURN
4841 END SUBROUTINE MINIMIZE_LOAD_IMBALANCE
4842
4843
4844
4845
4846
4847
4848
4849
4850
4851
4852
4853
4854
4855
4856
4857
4858
4859
4860
4861
4862 SUBROUTINE REPORT_BEST_IJK_SIZE0
4863
4864
4865
4866
4867 USE gridmap
4868 USE param
4869 USE param1
4870 USE constant
4871 USE run
4872 USE physprop
4873 USE indices
4874 USE scalars
4875 USE funits
4876 USE leqsol
4877 USE compar
4878 USE mpi_utility
4879 USE bc
4880 USE DISCRETELEMENT
4881
4882 USE parallel
4883
4884 USE cutcell
4885 USE quadric
4886 USE vtk
4887 USE polygon
4888 USE dashboard
4889 USE stl
4890
4891
4892 IMPLICIT NONE
4893
4894
4895
4896
4897
4898
4899
4900
4901
4902 INTEGER :: I,J,K,TOTAL_NUC,IDEAL_NCPP
4903
4904 INTEGER, DIMENSION(0:DIM_I) :: NUC_I,GLOBAL_NUC_I
4905 INTEGER, DIMENSION(0:DIM_J) :: NUC_J,GLOBAL_NUC_J
4906 INTEGER, DIMENSION(0:DIM_K) :: NUC_K,GLOBAL_NUC_K
4907
4908 INTEGER :: IPROC,PSUM
4909
4910 INTEGER, DIMENSION(0:numPEs-1) :: NCPP_OLD,NCPP,NCPP_OLD_WITH_GHOST,NCPP_WITH_GHOST
4911
4912 INTEGER, DIMENSION(0:NODESJ-1) :: JSIZE_OLD
4913
4914 INTEGER :: JSIZE, JREMAIN
4915 INTEGER :: MAXVAL_NCPP_OLD,MINVAL_NCPP_OLD,MAXVAL_NCPP,MINVAL_NCPP
4916 INTEGER :: AVG_NCPP_OLD,AVG_NCPP
4917 DOUBLE PRECISION :: LIP_OLD,LIP,MAXSPEEDUP_OLD,MAXSPEEDUP,P
4918
4919 INTEGER :: I_OFFSET,J_OFFSET,K_OFFSET,IERR
4920
4921 INTEGER, DIMENSION(0:numPEs-1) :: disp,rcount
4922
4923 INTEGER :: IPROC_OF_MAX_OLD,IPROC_OF_MIN_OLD
4924
4925
4926
4927 RETURN
4928
4929 IF(.NOT.CARTESIAN_GRID) RETURN
4930
4931 IF(NODESI*NODESJ*NODESK==1) RETURN
4932
4933
4934
4935 IF(.not.allocated(ISIZE_ALL)) allocate( ISIZE_ALL(0:NODESI-1))
4936 IF(.not.allocated(JSIZE_ALL)) allocate( JSIZE_ALL(0:NODESJ-1))
4937 IF(.not.allocated(KSIZE_ALL)) allocate( KSIZE_ALL(0:NODESK-1))
4938
4939
4940
4941
4942
4943
4944
4945
4946 IF(NODESI>1) THEN
4947
4948 IF(myPE == 0.AND.NODESJ*NODESK/=1) THEN
4949 WRITE(*,*)'ERROR IN SUBROUTINE REPORT_BEST_IJK_SIZE.'
4950 WRITE(*,*)'ADJUSTMENT POSSIBLE ONLY FOR DOMAIN DECOMPOSITION In ONE DIRECTION.'
4951 WRITE(*,*)'NODESI,NODESJ,NODESK =',NODESI,NODESJ,NODESK
4952 WRITE(*,*)'MFIX WILL EXIT NOW.'
4953 CALL MFIX_EXIT(myPE)
4954 ENDIF
4955
4956
4957 JSIZE_ALL(0:NODESJ-1) = jmax1-jmin1+1
4958 (0:NODESK-1) = kmax1-kmin1+1
4959
4960
4961 IF (myPE == 0) THEN
4962 WRITE(*,1000)'INFO FROM REPORT_BEST_IJK_SIZE:'
4963 WRITE(*,1000)'ATTEMPTING TO REPORT BEST DOMAIN SIZE IN I-DIRECTION ...'
4964 ENDIF
4965
4966 DO I = ISTART1,IEND1
4967
4968 NUC_I(I) = 0
4969
4970 DO J = JSTART3,JEND3
4971 DO K = KSTART3,KEND3
4972
4973 IF( .NOT.DEAD_CELL_AT(I,J,K)) THEN
4974 (I) = NUC_I(I) + 1
4975 ENDIF
4976
4977 ENDDO
4978 ENDDO
4979
4980 ENDDO
4981
4982
4983
4984 CALL allgather_1i (IEND1-ISTART1+1,rcount,IERR)
4985
4986 IF (myPE == 0) THEN
4987 I_OFFSET = 0
4988 ELSE
4989 I_OFFSET = 0
4990 DO iproc=0,myPE-1
4991 I_OFFSET = I_OFFSET + rcount(iproc)
4992 ENDDO
4993 ENDIF
4994
4995 CALL allgather_1i (I_OFFSET,disp,IERR)
4996
4997 call gatherv_1i( NUC_I(ISTART1:IEND1), IEND1-ISTART1+1, GLOBAL_NUC_I(IMIN1:IMAX1), rcount, disp, PE_IO, ierr )
4998
4999
5000
5001 ELSEIF(NODESJ>1) THEN
5002
5003
5004 IF(myPE == 0.AND.NODESI*NODESK/=1) THEN
5005 WRITE(*,*)'ERROR IN SUBROUTINE REPORT_BEST_IJK_SIZE.'
5006 WRITE(*,*)'ADJUSTMENT POSSIBLE ONLY FOR DOMAIN DECOMPOSITION In ONE DIRECTION.'
5007 WRITE(*,*)'NODESI,NODESJ,NODESK =',NODESI,NODESJ,NODESK
5008 WRITE(*,*)'MFIX WILL EXIT NOW.'
5009 CALL MFIX_EXIT(myPE)
5010 ENDIF
5011
5012
5013 ISIZE_ALL(0:NODESI-1) = imax1-imin1+1
5014 (0:NODESK-1) = kmax1-kmin1+1
5015
5016
5017 IF (myPE == 0) THEN
5018 WRITE(*,1000)'INFO FROM REPORT_BEST_IJK_SIZE:'
5019 WRITE(*,1000)'ATTEMPTING TO REPORT BEST DOMAIN SIZE IN J-DIRECTION ...'
5020 ENDIF
5021
5022 DO J = JSTART1,JEND1
5023
5024 NUC_J(J) = 0
5025
5026 DO I = ISTART3,IEND3
5027 DO K = KSTART3,KEND3
5028
5029 IF( .NOT.DEAD_CELL_AT(I,J,K)) THEN
5030 (J) = NUC_J(J) + 1
5031 ENDIF
5032
5033 ENDDO
5034 ENDDO
5035
5036 ENDDO
5037
5038
5039
5040 IF(IS_SERIAL) THEN
5041
5042 DO J=JMIN1,JMAX1
5043 GLOBAL_NUC_J(J) = NUC_J(J)
5044 ENDDO
5045
5046
5047 ELSE
5048
5049 CALL allgather_1i (JEND1-JSTART1+1,rcount,IERR)
5050
5051 IF (myPE == 0) THEN
5052 J_OFFSET = 0
5053 ELSE
5054 J_OFFSET = 0
5055 DO iproc=0,myPE-1
5056 J_OFFSET = J_OFFSET + rcount(iproc)
5057 ENDDO
5058 ENDIF
5059
5060 CALL allgather_1i (J_OFFSET,disp,IERR)
5061
5062 call gatherv_1i( NUC_J(JSTART1:JEND1), JEND1-JSTART1+1, GLOBAL_NUC_J(JMIN1:JMAX1), rcount, disp, PE_IO, ierr )
5063
5064 ENDIF
5065
5066
5067
5068
5069
5070
5071
5072
5073
5074 ELSEIF(NODESK>1) THEN
5075
5076 IF(myPE == 0.AND.NODESI*NODESJ/=1) THEN
5077 WRITE(*,*)'ERROR IN SUBROUTINE REPORT_BEST_IJK_SIZE.'
5078 WRITE(*,*)'ADJUSTMENT POSSIBLE ONLY FOR DOMAIN DECOMPOSITION In ONE DIRECTION.'
5079 WRITE(*,*)'NODESI,NODESJ,NODESK =',NODESI,NODESJ,NODESK
5080 WRITE(*,*)'MFIX WILL EXIT NOW.'
5081 CALL MFIX_EXIT(myPE)
5082 ENDIF
5083
5084
5085 ISIZE_ALL(0:NODESI-1) = imax1-imin1+1
5086 (0:NODESJ-1) = jmax1-jmin1+1
5087
5088
5089 IF (myPE == 0) THEN
5090 WRITE(*,1000)'INFO FROM REPORT_BEST_IJK_SIZE:'
5091 WRITE(*,1000)'ATTEMPTING TO REPORT BEST DOMAIN SIZE IN K-DIRECTION ...'
5092 ENDIF
5093
5094 DO K = KSTART1,KEND1
5095
5096 NUC_K(K) = 0
5097
5098 DO I = ISTART1,IEND1
5099 DO J = JSTART1,JEND1
5100
5101 IF( .NOT.DEAD_CELL_AT(I,J,K)) THEN
5102 (K) = NUC_K(K) + 1
5103 ENDIF
5104
5105 ENDDO
5106 ENDDO
5107
5108 ENDDO
5109
5110
5111
5112 CALL allgather_1i (KEND1-KSTART1+1,rcount,IERR)
5113
5114 IF (myPE == 0) THEN
5115 K_OFFSET = 0
5116 ELSE
5117 K_OFFSET = 0
5118 DO iproc=0,myPE-1
5119 K_OFFSET = K_OFFSET + rcount(iproc)
5120 ENDDO
5121 ENDIF
5122
5123 CALL allgather_1i (K_OFFSET,disp,IERR)
5124
5125 call gatherv_1i( NUC_K(KSTART1:KEND1), KEND1-KSTART1+1, GLOBAL_NUC_K(KMIN1:KMAX1), rcount, disp, PE_IO, ierr )
5126
5127
5128
5129 ENDIF
5130
5131
5132
5133
5134 IF (myPE == PE_IO) THEN
5135
5136 IF(NODESI>1) THEN
5137
5138
5139
5140 ELSEIF(NODESJ>1) THEN
5141
5142
5143
5144 = (JMAX1-JMIN1+1)/NODESJ
5145 JSIZE_OLD(0:NODESJ-1) = JSIZE
5146
5147 JREMAIN = (JMAX1-JMIN1+1) - NODESJ*JSIZE
5148 IF (JREMAIN.GE.1) THEN
5149 JSIZE_OLD( 0:(JREMAIN-1) ) = JSIZE + 1
5150 ENDIF
5151
5152
5153
5154 JSIZE_ALL = JSIZE_OLD
5155
5156 CALL GET_LIP_WITH_GHOST_LAYERS(NODESJ,GLOBAL_NUC_J(JMIN1:JMAX1),JMIN1,JMAX1,JSIZE_ALL, &
5157 NCPP_OLD,NCPP_OLD_WITH_GHOST,LIP_OLD,IPROC_OF_MAX_OLD,IPROC_OF_MIN_OLD)
5158
5159
5160
5161
5162
5163
5164
5165
5166
5167
5168
5169 CALL MINIMIZE_LOAD_IMBALANCE(nodesj,GLOBAL_NUC_J(JMIN1:JMAX1),JMIN1,JMAX1,JSIZE_ALL,NCPP,NCPP_WITH_GHOST)
5170
5171
5172 TOTAL_NUC = SUM(NCPP_WITH_GHOST(0:NumPEs-1))
5173 IDEAL_NCPP = TOTAL_NUC / NumPEs
5174
5175 WRITE (*, 1000) 'AFTER OPTIMIZATION:'
5176 WRITE (*, 1010) 'TOTAL NUMBER OF USEFUL CELLS = ',TOTAL_NUC
5177 WRITE (*, 1010) 'IDEAL NUMBER OF CELLS/PROC. = ',IDEAL_NCPP
5178 WRITE (*, 1000) 'ACTUALL CELL COUNT:'
5179 WRITE (*, 1000) '================================================='
5180 WRITE (*, 1000) ' PROCESSOR J-SIZE CELLS/PROC. DIFF. (%)'
5181 WRITE (*, 1000) '================================================='
5182 DO IPROC = 0,numPEs-1
5183 WRITE (*, 1020) IPROC,JSIZE_ALL(IPROC),NCPP_WITH_GHOST(IPROC), &
5184 DBLE(NCPP_WITH_GHOST(IPROC)-IDEAL_NCPP)/DBLE(IDEAL_NCPP)*100.0D0
5185 ENDDO
5186
5187
5188
5189
5190
5191
5192
5193
5194
5195
5196 = 0
5197 DO IPROC = 0,numPEs-1
5198 PSUM = PSUM + JSIZE_ALL(IPROC)
5199 IF(JSIZE_ALL(IPROC)<5) THEN
5200 WRITE (*, 1010) 'ERROR: J-SIZE TOO SMALL FOR PROCESSOR:',IPROC
5201 WRITE (*, 1010) 'J-SIZE = ',JSIZE_ALL(IPROC)
5202 CALL MFIX_EXIT(myPE)
5203 ENDIF
5204
5205 ENDDO
5206
5207 IF(PSUM/=JMAX) THEN
5208 WRITE (*, 1000) 'ERROR IN ADJUST_IJK_SIZE: UNABLE TO ASSIGN JSIZE TO PROCESSORS.'
5209 WRITE (*, 1000) 'SUM OF JSIZE_ALL DOES NOT MATCH JMAX:'
5210 WRITE (*, 1010) 'SUM OF JSIZE_ALL = ',PSUM
5211 WRITE (*, 1010) 'JMAX1 = ',JMAX
5212 CALL MFIX_EXIT(myPE)
5213 ENDIF
5214
5215
5216 OPEN(CONVERT='BIG_ENDIAN',UNIT=777, FILE='suggested_gridmap.dat')
5217 WRITE (777, 1000) 'J-SIZE DISTRIBUTION'
5218 WRITE (777, 1010) 'NUMBER OF PROCESSORS = ',NumPEs
5219 WRITE (777, 1000) '================================================='
5220 WRITE (777, 1000) ' PROCESSOR J-SIZE'
5221 WRITE (777, 1000) '================================================='
5222
5223 DO IPROC = 0,NumPEs-1
5224 WRITE(777,1060) IPROC,jsize_all(IPROC)
5225 ENDDO
5226 CLOSE(777)
5227 WRITE (*, 1000) '================================================='
5228 WRITE (*, 1000) 'J-SIZE DISTRIBUTION SAVED IN FILE: suggested_gridmap.dat'
5229 WRITE (*, 1000) 'TO USE THIS DISTRIBUTION, RENAME THE FILE AS: gridmap.dat'
5230 WRITE (*, 1000) 'AND RUN MFIX AGAIN.'
5231 WRITE (*, 1000) '================================================='
5232
5233
5234 ELSEIF(NODESK>1) THEN
5235
5236
5237
5238
5239 ENDIF
5240
5241
5242 = MAXVAL(NCPP_OLD_WITH_GHOST)
5243 MINVAL_NCPP_OLD = MINVAL(NCPP_OLD_WITH_GHOST)
5244 AVG_NCPP_OLD = SUM(NCPP_OLD_WITH_GHOST)/NUMPES
5245
5246
5247 = DBLE(MAXVAL_NCPP_OLD-MINVAL_NCPP_OLD)/DBLE(MINVAL_NCPP_OLD)*100.0D0
5248
5249
5250 = DBLE(MINVAL_NCPP_OLD)/DBLE(MAXVAL_NCPP_OLD)
5251
5252
5253 = ONE / ((ONE-P) + P/NumPes)
5254
5255 MAXVAL_NCPP = MAXVAL(NCPP_WITH_GHOST)
5256 MINVAL_NCPP = MINVAL(NCPP_WITH_GHOST)
5257 AVG_NCPP = SUM(NCPP_WITH_GHOST(0:NumPEs-1))/NUMPES
5258
5259
5260 = DBLE(MAXVAL_NCPP-MINVAL_NCPP)/DBLE(MINVAL_NCPP)*100.0D0
5261
5262
5263 = DBLE(MINVAL_NCPP)/DBLE(MAXVAL_NCPP)
5264
5265
5266 = ONE / ((ONE-P) + P/NumPes)
5267
5268 WRITE (*, 1000) '================================================='
5269 WRITE (*, 1000) 'ESTIMATED PARALLEL LOAD BALANCING STATISTICS'
5270 WRITE (*, 1000) 'COMPARISION BETWEEN UNIFORM SIZE (OLD)'
5271 WRITE (*, 1000) 'AND SUGGESTED SIZE (NEW)'
5272 WRITE (*, 1000) '================================================='
5273 WRITE (*, 1000) ' OLD NEW'
5274 WRITE (*, 1010) 'MAX CELL COUNT : ',MAXVAL_NCPP_OLD,MAXVAL_NCPP
5275 WRITE (*, 1010) 'AT PROCESSOR : ',MAXLOC(NCPP_OLD_WITH_GHOST)-1,MAXLOC(NCPP_WITH_GHOST)-1
5276 WRITE (*, 1010) 'MIN CELL COUNT : ',MINVAL_NCPP_OLD,MINVAL_NCPP
5277 WRITE (*, 1010) 'AT PROCESSOR : ',MINLOC(NCPP_OLD_WITH_GHOST)-1,MINLOC(NCPP_WITH_GHOST)-1
5278 WRITE (*, 1010) 'AVG CELL COUNT : ',AVG_NCPP_OLD,AVG_NCPP
5279 WRITE (*, 1000) ''
5280 WRITE (*, 1030) 'LOAD IMBALANCE (%) : ',LIP_OLD,LIP
5281 WRITE (*, 1000) ''
5282
5283
5284
5285
5286 WRITE (*, 1000) '================================================='
5287
5288
5289
5290 ENDIF
5291
5292
5293
5294
5295
5296
5297
5298
5299
5300
5301
5302
5303
5304
5305 FORMAT(1x,A)
5306 1010 FORMAT(1x,A,I10,I10)
5307 1020 FORMAT(1X,I8,2(I12),F12.2)
5308 1030 FORMAT(1X,A,2(F10.1))
5309 1040 FORMAT(F10.1)
5310 1050 FORMAT(1X,3(A))
5311 1060 FORMAT(1x,I8,I12)
5312
5313 RETURN
5314
5315 END SUBROUTINE REPORT_BEST_IJK_SIZE0
5316
5317
5318
5319
5320
5321
5322
5323
5324
5325
5326
5327
5328
5329
5330
5331
5332
5333 SUBROUTINE GET_LIP_WITH_GHOST_LAYERS0(NODESL,NUC_L,LMIN1,LMAX1,L_SIZE,NCPP,NCPP_WITH_GHOST,LIP,IPROC_OF_MAX,IPROC_OF_MIN)
5334
5335
5336
5337
5338 USE gridmap
5339 USE param
5340 USE param1
5341 USE constant
5342 USE run
5343 USE physprop
5344 USE indices
5345 USE scalars
5346 USE funits
5347 USE leqsol
5348 USE compar
5349 USE mpi_utility
5350 USE bc
5351 USE DISCRETELEMENT
5352
5353 USE cutcell
5354 USE quadric
5355 USE vtk
5356 USE polygon
5357 USE dashboard
5358 USE stl
5359
5360
5361 IMPLICIT NONE
5362
5363
5364
5365
5366
5367
5368
5369
5370
5371 INTEGER :: NODESL,L,LMIN1,LMAX1,TOTAL_NUC,TOTAL_NUC_WITH_GHOST,IPROC_OF_MAX,IPROC_OF_MIN
5372
5373 INTEGER :: LCOUNT1,LCOUNT2,MINVAL_NCPP,MAXVAL_NCPP,IDEAL_NCPP
5374 INTEGER, DIMENSION(LMIN1:LMAX1) :: NUC_L
5375
5376 INTEGER :: IPROC
5377
5378 INTEGER, DIMENSION(0:NODESL-1) :: NCPP,NCPP_WITH_GHOST,L_SIZE,L1,L2
5379
5380
5381 DOUBLE PRECISION :: LIP
5382
5383
5384
5385 = LMAX1 - LMIN1 + 1
5386 LCOUNT2 = SUM(L_SIZE(0:NODESL-1))
5387
5388 IF(LCOUNT1/=LCOUNT2) THEN
5389 WRITE(*,*)' ERROR: SUM OF CELLS DO NOT MATCH:',LCOUNT1,LCOUNT2
5390 CALL MFIX_EXIT(myPE)
5391 ENDIF
5392
5393 L1(0) = LMIN1
5394 L2(0) = L1(0) + L_SIZE(0) - 1
5395
5396 DO IPROC = 1,NODESL-1
5397 L1(IPROC) = L2(IPROC-1) + 1
5398 L2(IPROC) = L1(IPROC) + L_SIZE(IPROC) - 1
5399 ENDDO
5400
5401 DO IPROC = 0,NODESL-1
5402 NCPP(IPROC) = SUM(NUC_L(L1(IPROC):L2(IPROC)))
5403
5404
5405 ENDDO
5406
5407 TOTAL_NUC = 0
5408
5409 DO L=LMIN1,LMAX1
5410 TOTAL_NUC = TOTAL_NUC + NUC_L(L)
5411 ENDDO
5412
5413 NCPP_WITH_GHOST(0) = NCPP(0) + 2*NUC_L(L1(0)) + NUC_L(L1(1)) + NUC_L(L1(1)+1)
5414
5415 DO IPROC = 1,NODESL-2
5416 NCPP_WITH_GHOST(IPROC) = NCPP(IPROC) &
5417 + NUC_L(L2(IPROC-1)) + NUC_L(L2(IPROC-1)-1) &
5418 + NUC_L(L1(IPROC+1)) + NUC_L(L1(IPROC+1)+1)
5419 ENDDO
5420
5421 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)
5422
5423 TOTAL_NUC_WITH_GHOST = 0
5424 DO IPROC = 0,NODESL-1
5425
5426 = TOTAL_NUC_WITH_GHOST + NCPP_WITH_GHOST(IPROC)
5427 ENDDO
5428
5429
5430 IDEAL_NCPP = TOTAL_NUC_WITH_GHOST / NumPEs
5431
5432 MAXVAL_NCPP = MAXVAL(NCPP_WITH_GHOST)
5433 MINVAL_NCPP = MINVAL(NCPP_WITH_GHOST)
5434
5435
5436 = DBLE(MAXVAL_NCPP-MINVAL_NCPP)/DBLE(MINVAL_NCPP)*100.0D0
5437
5438
5439 IPROC_OF_MAX = MAXLOC(NCPP_WITH_GHOST,1)-1
5440 IPROC_OF_MIN = MINLOC(NCPP_WITH_GHOST,1)-1
5441
5442
5443
5444
5445
5446 RETURN
5447 END SUBROUTINE GET_LIP_WITH_GHOST_LAYERS0
5448
5449
5450
5451
5452
5453
5454
5455
5456
5457
5458
5459
5460
5461
5462
5463
5464
5465
5466
5467 SUBROUTINE MINIMIZE_LOAD_IMBALANCE0(NODESL,NUC_L,LMIN1,LMAX1,L_SIZE,NCPP,NCPP_WITH_GHOST)
5468
5469
5470
5471
5472 USE gridmap
5473 USE param
5474 USE param1
5475 USE constant
5476 USE run
5477 USE physprop
5478 USE indices
5479 USE scalars
5480 USE funits
5481 USE leqsol
5482 USE compar
5483 USE mpi_utility
5484 USE bc
5485 USE DISCRETELEMENT
5486
5487 USE cutcell
5488 USE quadric
5489 USE vtk
5490 USE polygon
5491 USE dashboard
5492 USE stl
5493
5494
5495 IMPLICIT NONE
5496
5497
5498
5499
5500
5501
5502
5503
5504
5505 INTEGER :: NODESL,LMIN1,LMAX1,IPROC_OF_MAX,IPROC_OF_MIN
5506
5507 INTEGER, DIMENSION(LMIN1:LMAX1) :: NUC_L
5508
5509 INTEGER :: N,NOIMPROVEMENT
5510
5511 INTEGER,PARAMETER :: NAMAX=10000
5512
5513 INTEGER, DIMENSION(0:numPEs-1) :: NCPP,NCPP_WITH_GHOST,L_SIZE,BEST_L_SIZE,BEST_NCPP,BEST_NCPP_WITH_GHOST
5514
5515
5516 DOUBLE PRECISION :: LIP,BEST_LIP
5517
5518
5519
5520
5521
5522
5523
5524
5525 CALL GET_LIP_WITH_GHOST_LAYERS0(NODESL,NUC_L,LMIN1,LMAX1,L_SIZE,NCPP,NCPP_WITH_GHOST,LIP,IPROC_OF_MAX,IPROC_OF_MIN)
5526
5527 BEST_LIP = LIP
5528 BEST_L_SIZE = L_SIZE
5529
5530
5531
5532
5533
5534
5535
5536
5537
5538
5539
5540
5541
5542
5543
5544 print*,'ATTEMPTING TO OPTIMIZE LOAD BALANCE...'
5545
5546 NOIMPROVEMENT=0
5547
5548 DO N = 1,NAMAX
5549
5550 L_SIZE(IPROC_OF_MAX) = L_SIZE(IPROC_OF_MAX) - 1
5551 L_SIZE(IPROC_OF_MIN) = L_SIZE(IPROC_OF_MIN) + 1
5552
5553 CALL GET_LIP_WITH_GHOST_LAYERS(NODESL,NUC_L,LMIN1,LMAX1,L_SIZE,NCPP,NCPP_WITH_GHOST,LIP,IPROC_OF_MAX,IPROC_OF_MIN)
5554
5555
5556
5557
5558
5559 IF(LIP<BEST_LIP) THEN
5560 BEST_LIP = LIP
5561 BEST_L_SIZE = L_SIZE
5562 BEST_NCPP = NCPP
5563 BEST_NCPP_WITH_GHOST = NCPP_WITH_GHOST
5564 NOIMPROVEMENT=0
5565 ELSE
5566 NOIMPROVEMENT = NOIMPROVEMENT + 1
5567 ENDIF
5568
5569
5570
5571
5572
5573
5574
5575
5576
5577
5578 IF(NOIMPROVEMENT==10) THEN
5579 WRITE (*, 1000) 'OPTIMIZED LOAD BALANCE REACHED.'
5580 EXIT
5581 ENDIF
5582
5583 ENDDO
5584
5585
5586 = BEST_L_SIZE
5587 NCPP = BEST_NCPP
5588 NCPP_WITH_GHOST = BEST_NCPP_WITH_GHOST
5589
5590
5591
5592
5593
5594
5595
5596
5597
5598
5599 FORMAT(1x,A)
5600 1010 FORMAT(1x,A,I10,I10)
5601 1020 FORMAT(1X,I8,2(I12),F12.2)
5602
5603 RETURN
5604 END SUBROUTINE MINIMIZE_LOAD_IMBALANCE0
5605
5606
5607