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