File: /nfs/home/0/users/jenkins/mfix.git/model/cartesian_grid/get_delh.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name:  GET_DEL_H                                             C
4     !  Purpose: Finds the normal distance and unit vector from a cut face  C
5     !  to any point (x0,y0,z0)                                             C
6     !  The unit normal vector points away from the boundary,               C
7     !  towards the fluid.                                                  C
8     !  This subroutine must be called from a cut-cell                      C
9     !                                                                      C
10     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
11     !  Reviewer:                                          Date:            C
12     !                                                                      C
13     !  Revision Number #                                  Date: ##-###-##  C
14     !  Author: #                                                           C
15     !  Purpose: #                                                          C
16     !                                                                      C
17     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
18       SUBROUTINE GET_DEL_H(IJK,TYPE_OF_CELL,X0,Y0,Z0,Del_H,Nx,Ny,Nz)
19     
20           USE param
21           USE param1
22           USE parallel
23           USE constant
24           USE run
25           USE toleranc
26           USE geometry
27           USE indices
28           USE compar
29           USE sendrecv
30           USE quadric
31           USE cutcell
32           USE functions
33     
34           IMPLICIT NONE
35           CHARACTER (LEN=*) :: TYPE_OF_CELL
36           DOUBLE PRECISION:: X0,Y0,Z0,XREF,YREF,ZREF
37           INTEGER :: IJK,I,J,K
38           DOUBLE PRECISION :: Del_H,Diagonal
39           DOUBLE PRECISION :: Nx,Ny,Nz
40     
41           SELECT CASE (TYPE_OF_CELL)
42              CASE('SCALAR')
43     
44                 IF(.NOT.CUT_CELL_AT(IJK)) THEN
45                    WRITE(*,*)' EROR IN SUBROUTINE GET_DEL_H:'
46                    WRITE(*,*)' SCALAR CELL',IJK,' IS NOT A CUT CELL'
47                    WRITE(*,*)' MFiX will exit now.'
48                    CALL MFIX_EXIT(myPE)
49                 ENDIF
50     
51                 Nx = NORMAL_S(IJK,1)
52                 Ny = NORMAL_S(IJK,2)
53                 Nz = NORMAL_S(IJK,3)
54     
55                 Xref = REFP_S(IJK,1)
56                 Yref = REFP_S(IJK,2)
57                 Zref = REFP_S(IJK,3)
58     
59              CASE('U_MOMENTUM')
60     
61                 IF(.NOT.CUT_U_CELL_AT(IJK)) THEN
62                    WRITE(*,*)' EROR IN SUBROUTINE GET_DEL_H:'
63                    WRITE(*,*)' U-MOMENTUM CELL',IJK,' IS NOT A CUT CELL'
64                    WRITE(*,*)' MFiX will exit now.'
65                    CALL MFIX_EXIT(myPE)
66                 ENDIF
67     
68                 Nx = NORMAL_U(IJK,1)
69                 Ny = NORMAL_U(IJK,2)
70                 Nz = NORMAL_U(IJK,3)
71     
72                 Xref = REFP_U(IJK,1)
73                 Yref = REFP_U(IJK,2)
74                 Zref = REFP_U(IJK,3)
75     
76                 IF(WALL_U_AT(IJK)) THEN
77                    Nx = ZERO
78                    Ny = ZERO
79                    Nz = ZERO
80                    DEL_H = UNDEFINED
81                    RETURN
82                 ENDIF
83     
84              CASE('V_MOMENTUM')
85     
86                 IF(.NOT.CUT_V_CELL_AT(IJK)) THEN
87                    WRITE(*,*)' EROR IN SUBROUTINE GET_DEL_H:'
88                    WRITE(*,*)' V-MOMENTUM CELL',IJK,' IS NOT A CUT CELL'
89                    WRITE(*,*)' MFiX will exit now.'
90                    CALL MFIX_EXIT(myPE)
91                 ENDIF
92     
93                 Nx = NORMAL_V(IJK,1)
94                 Ny = NORMAL_V(IJK,2)
95                 Nz = NORMAL_V(IJK,3)
96     
97                 Xref = REFP_V(IJK,1)
98                 Yref = REFP_V(IJK,2)
99                 Zref = REFP_V(IJK,3)
100     
101                 IF(WALL_V_AT(IJK)) THEN
102                    Nx = ZERO
103                    Ny = ZERO
104                    Nz = ZERO
105                    DEL_H = UNDEFINED
106                    RETURN
107                 ENDIF
108     
109              CASE('W_MOMENTUM')
110     
111                 IF(.NOT.CUT_W_CELL_AT(IJK)) THEN
112                    WRITE(*,*)' EROR IN SUBROUTINE GET_DEL_H:'
113                    WRITE(*,*)' W-MOMENTUM CELL',IJK,' IS NOT A CUT CELL'
114                    WRITE(*,*)' MFiX will exit now.'
115                    CALL MFIX_EXIT(myPE)
116                 ENDIF
117     
118                 Nx = NORMAL_W(IJK,1)
119                 Ny = NORMAL_W(IJK,2)
120                 Nz = NORMAL_W(IJK,3)
121     
122                 Xref = REFP_W(IJK,1)
123                 Yref = REFP_W(IJK,2)
124                 Zref = REFP_W(IJK,3)
125     
126                 IF(WALL_W_AT(IJK)) THEN
127                    Nx = ZERO
128                    Ny = ZERO
129                    Nz = ZERO
130                    DEL_H = UNDEFINED
131                    RETURN
132                 ENDIF
133     
134              CASE DEFAULT
135                 WRITE(*,*)' EROR IN SUBROUTINE GET_DEL_H:'
136                 WRITE(*,*)'UNKNOWN TYPE OF CELL:',TYPE_OF_CELL
137                 WRITE(*,*)'ACCEPTABLE TYPES ARE:'
138                 WRITE(*,*)'SCALAR'
139                 WRITE(*,*)'U_MOMENTUM'
140                 WRITE(*,*)'V_MOMENTUM'
141                 WRITE(*,*)'W_MOMENTUM'
142                 CALL MFIX_EXIT(myPE)
143           END SELECT
144     
145     
146           DEL_H = Nx * (X0 - Xref) + Ny * (Y0 - Yref) + Nz * (Z0 - Zref)
147     
148     !======================================================================
149     ! Negative values of DEL_H are not accepted
150     !======================================================================
151     
152            I = I_OF(IJK)
153            J = J_OF(IJK)
154            K = K_OF(IJK)
155     
156            IF(NO_K) THEN
157               Diagonal = sqrt(DX(I)**2 + DY(J)**2 )
158            ELSE
159               Diagonal = sqrt(DX(I)**2 + DY(J)**2 + DZ(K)**2)
160            ENDIF
161     
162           IF (DEL_H <= TOL_DELH * Diagonal) THEN
163     
164              DEL_H = UNDEFINED
165              Nx = ZERO
166              Ny = ZERO
167              Nz = ZERO
168     
169           ENDIF
170     
171           RETURN
172     
173     
174           END SUBROUTINE GET_DEL_H
175     
176       SUBROUTINE GET_DEL_H_DES(IJK,TYPE_OF_CELL,X0,Y0,Z0,Del_H,Nx,Ny,Nz, allow_neg_dist)
177     
178           USE param
179           USE param1
180           USE parallel
181           USE constant
182           USE run
183           USE toleranc
184           USE geometry
185           USE indices
186           USE compar
187           USE sendrecv
188           USE quadric
189           USE cutcell
190           USE functions
191     
192           IMPLICIT NONE
193           CHARACTER (LEN=*) :: TYPE_OF_CELL
194           DOUBLE PRECISION:: X0,Y0,Z0,XREF,YREF,ZREF
195           INTEGER :: IJK,I,J,K
196           DOUBLE PRECISION :: Del_H,Diagonal
197           DOUBLE PRECISION :: Nx,Ny,Nz
198     
199           LOGICAL :: ALLOW_NEG_DIST
200     
201           SELECT CASE (TYPE_OF_CELL)
202              CASE('SCALAR')
203     
204                 IF(.NOT.CUT_CELL_AT(IJK)) THEN
205                    WRITE(*,*)' EROR IN SUBROUTINE GET_DEL_H_DES:'
206                    WRITE(*,*)' SCALAR CELL',IJK,' IS NOT A CUT CELL'
207                    WRITE(*,*)' I, J, K =',I_OF(IJK), J_OF(IJK), K_OF(IJK)
208                    WRITE(*,*)' MFiX will exit now.'
209                    CALL MFIX_EXIT(myPE)
210                 ENDIF
211     
212                 Nx = NORMAL_S(IJK,1)
213                 Ny = NORMAL_S(IJK,2)
214                 Nz = NORMAL_S(IJK,3)
215     
216                 Xref = REFP_S(IJK,1)
217                 Yref = REFP_S(IJK,2)
218                 Zref = REFP_S(IJK,3)
219     
220              CASE('U_MOMENTUM')
221     
222                 IF(.NOT.CUT_U_CELL_AT(IJK)) THEN
223                    WRITE(*,*)' EROR IN SUBROUTINE GET_DEL_H_DES:'
224                    WRITE(*,*)' U-MOMENTUM CELL',IJK,' IS NOT A CUT CELL'
225                    WRITE(*,*)' MFiX will exit now.'
226                    CALL MFIX_EXIT(myPE)
227                 ENDIF
228     
229                 Nx = NORMAL_U(IJK,1)
230                 Ny = NORMAL_U(IJK,2)
231                 Nz = NORMAL_U(IJK,3)
232     
233                 Xref = REFP_U(IJK,1)
234                 Yref = REFP_U(IJK,2)
235                 Zref = REFP_U(IJK,3)
236     
237                 IF(WALL_U_AT(IJK)) THEN
238                    Nx = ZERO
239                    Ny = ZERO
240                    Nz = ZERO
241                    DEL_H = UNDEFINED
242                    RETURN
243                 ENDIF
244     
245              CASE('V_MOMENTUM')
246     
247                 IF(.NOT.CUT_V_CELL_AT(IJK)) THEN
248                    WRITE(*,*)' EROR IN SUBROUTINE GET_DEL_H_DES:'
249                    WRITE(*,*)' V-MOMENTUM CELL',IJK,' IS NOT A CUT CELL'
250                    WRITE(*,*)' MFiX will exit now.'
251                    CALL MFIX_EXIT(myPE)
252                 ENDIF
253     
254                 Nx = NORMAL_V(IJK,1)
255                 Ny = NORMAL_V(IJK,2)
256                 Nz = NORMAL_V(IJK,3)
257     
258                 Xref = REFP_V(IJK,1)
259                 Yref = REFP_V(IJK,2)
260                 Zref = REFP_V(IJK,3)
261     
262                 IF(WALL_V_AT(IJK)) THEN
263                    Nx = ZERO
264                    Ny = ZERO
265                    Nz = ZERO
266                    DEL_H = UNDEFINED
267                    RETURN
268                 ENDIF
269     
270              CASE('W_MOMENTUM')
271     
272                 IF(.NOT.CUT_W_CELL_AT(IJK)) THEN
273                    WRITE(*,*)' EROR IN SUBROUTINE GET_DEL_H_DES:'
274                    WRITE(*,*)' W-MOMENTUM CELL',IJK,' IS NOT A CUT CELL'
275                    WRITE(*,*)' MFiX will exit now.'
276                    CALL MFIX_EXIT(myPE)
277                 ENDIF
278     
279                 Nx = NORMAL_W(IJK,1)
280                 Ny = NORMAL_W(IJK,2)
281                 Nz = NORMAL_W(IJK,3)
282     
283                 Xref = REFP_W(IJK,1)
284                 Yref = REFP_W(IJK,2)
285                 Zref = REFP_W(IJK,3)
286     
287                 IF(WALL_W_AT(IJK)) THEN
288                    Nx = ZERO
289                    Ny = ZERO
290                    Nz = ZERO
291                    DEL_H = UNDEFINED
292                    RETURN
293                 ENDIF
294     
295              CASE DEFAULT
296                 WRITE(*,*)' EROR IN SUBROUTINE GET_DEL_H_DES:'
297                 WRITE(*,*)'UNKNOWN TYPE OF CELL:',TYPE_OF_CELL
298                 WRITE(*,*)'ACCEPTABLE TYPES ARE:'
299                 WRITE(*,*)'SCALAR'
300                 WRITE(*,*)'U_MOMENTUM'
301                 WRITE(*,*)'V_MOMENTUM'
302                 WRITE(*,*)'W_MOMENTUM'
303                 CALL MFIX_EXIT(myPE)
304           END SELECT
305     
306           DEL_H = Nx * (X0 - Xref) + Ny * (Y0 - Yref) + Nz * (Z0 - Zref)
307     
308     !======================================================================
309     ! Negative values of DEL_H are not accepted
310     !======================================================================
311     
312            I = I_OF(IJK)
313            J = J_OF(IJK)
314            K = K_OF(IJK)
315     
316            IF(.NOT.ALLOW_NEG_DIST) THEN
317               Diagonal = sqrt(DX(I)**2 + DY(J)**2 + DZ(K)**2)
318     
319               IF (DEL_H <= TOL_DELH * Diagonal) THEN
320                  DEL_H = UNDEFINED
321                  Nx = ZERO
322                  Ny = ZERO
323                  Nz = ZERO
324     
325               ENDIF
326            ENDIF
327     
328           RETURN
329     
330           END SUBROUTINE GET_DEL_H_DES
331     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
332     !                                                                      C
333     !  Module name: STORE_CUT_FACE_INFO                                    C
334     !  Purpose: Compute and store unit normal vector and reference point   C
335     !           Defining a cut face                                        C
336     !                                                                      C
337     !  Author: Jeff Dietiker                              Date: 21-Feb-08  C
338     !  Reviewer:                                          Date:            C
339     !                                                                      C
340     !  Revision Number #                                  Date: ##-###-##  C
341     !  Author: #                                                           C
342     !  Purpose: #                                                          C
343     !                                                                      C
344     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
345       SUBROUTINE STORE_CUT_FACE_INFO(IJK,TYPE_OF_CELL,N_CUT_FACE_NODES,COORD_CUT_FACE_NODES,X_MEAN,Y_MEAN,Z_MEAN)
346     
347           USE param
348           USE param1
349           USE parallel
350           USE constant
351           USE run
352           USE toleranc
353           USE geometry
354           USE indices
355           USE compar
356           USE sendrecv
357           USE quadric
358           USE cutcell
359           USE functions
360     
361           IMPLICIT NONE
362           CHARACTER (LEN=*) :: TYPE_OF_CELL
363           INTEGER :: IJK
364           INTEGER :: N_CUT_FACE_NODES
365           DOUBLE PRECISION, DIMENSION(3,15) :: COORD_CUT_FACE_NODES
366           DOUBLE PRECISION :: X_MEAN,Y_MEAN,Z_MEAN
367           DOUBLE PRECISION, DIMENSION(3) :: N,V,N1,N2
368           DOUBLE PRECISION :: NORM,N1_dot_N2
369           DOUBLE PRECISION, DIMENSION(3) :: VECTEMP,VECTEMP2
370     
371          IF(NO_K) THEN  ! 2D case
372     
373             N(1) = COORD_CUT_FACE_NODES(2,1) - COORD_CUT_FACE_NODES(2,2) ! y1-y2
374             N(2) = COORD_CUT_FACE_NODES(1,2) - COORD_CUT_FACE_NODES(1,1) ! x2-x1
375             N(3) = ZERO
376     
377          ELSE
378     
379     !======================================================================
380     ! Make sure there are a least three points along the plane
381     !======================================================================
382     
383             IF(N_CUT_FACE_NODES < 3) THEN
384                WRITE(*,*)' ERROR IN SUBROUTINE STORE_CUT_FACE_INFO:'
385                WRITE(*,*)' CUT FACE HAS LESS THAN 3 NODES.'
386                WRITE(*,*)' MFIX WILL EXIT NOW.'
387                CALL MFIX_EXIT(myPE)
388             END IF
389     
390     !======================================================================
391     !  Find tentative unit normal vector
392     !  and reverse sign if necessary
393     !  (unit vector must be pointing towards the fluid)
394     !======================================================================
395     
396             VECTEMP  = COORD_CUT_FACE_NODES(:,2)-COORD_CUT_FACE_NODES(:,1)
397             VECTEMP2 = COORD_CUT_FACE_NODES(:,3)-COORD_CUT_FACE_NODES(:,1)
398             N = CROSS_PRODUCT(VECTEMP,VECTEMP2)
399     
400          ENDIF
401     
402           NORM = sqrt(dot_product(n(:),n(:)))
403           N = N / NORM
404     
405           V(1) = X_MEAN - COORD_CUT_FACE_NODES(1,1)
406           V(2) = Y_MEAN - COORD_CUT_FACE_NODES(1,2)
407           V(3) = Z_MEAN - COORD_CUT_FACE_NODES(1,3)
408     
409           IF (DOT_PRODUCT(N,V) < ZERO) N = - N
410     
411           IF(N_CUT_FACE_NODES > 3) THEN     ! FOR 3D geometry, check normal of plane defined by nodes 1,2, and 4
412     
413              N1 = N  ! Keep copy of previous N (nodes 1,2,3)
414     
415             VECTEMP  = COORD_CUT_FACE_NODES(:,2)-COORD_CUT_FACE_NODES(:,1)
416             VECTEMP2 = COORD_CUT_FACE_NODES(:,4)-COORD_CUT_FACE_NODES(:,1)
417     
418              N2 = CROSS_PRODUCT(VECTEMP,VECTEMP2)
419     
420              NORM = sqrt(dot_product(n2(:),n2(:)))
421              N2 = N2 / NORM
422     
423              V(1) = X_MEAN - COORD_CUT_FACE_NODES(1,1)
424              V(2) = Y_MEAN - COORD_CUT_FACE_NODES(2,1)
425              V(3) = Z_MEAN - COORD_CUT_FACE_NODES(3,1)
426     
427              IF (DOT_PRODUCT(N2,V) < ZERO) N2 = - N2
428     
429           ENDIF
430     
431           N1_dot_N2 = DOT_PRODUCT(N1,N2)
432           DEBUG_CG(IJK,1)=N1_dot_N2
433     
434           IF(N1_dot_N2<0.99) THEN
435     
436     !         What should be done when the unit vectors are different ?
437     
438           ENDIF
439     
440     !======================================================================
441     ! Store unit normal vector and reference point
442     !======================================================================
443     
444           SELECT CASE (TYPE_OF_CELL)
445              CASE('SCALAR')
446     
447                 NORMAL_S(IJK,:) = N
448                 REFP_S(IJK,:)   = COORD_CUT_FACE_NODES(:,1)
449     
450                 IF(DO_K) CALL TEST_DEL_H(IJK,'SCALAR') ! test for negative del_H
451     
452              CASE('U_MOMENTUM')
453     
454                 NORMAL_U(IJK,:) = N
455                 REFP_U(IJK,:)   = COORD_CUT_FACE_NODES(:,1)
456     
457              CASE('V_MOMENTUM')
458     
459                 NORMAL_V(IJK,:) = N
460                 REFP_V(IJK,:)   = COORD_CUT_FACE_NODES(:,1)
461     
462              CASE('W_MOMENTUM')
463     
464                 NORMAL_W(IJK,:) = N
465                 REFP_W(IJK,:)   = COORD_CUT_FACE_NODES(:,1)
466     
467              CASE DEFAULT
468                 WRITE(*,*)'SUBROUTINE: STORE_CUT_FACE_INFO'
469                 WRITE(*,*)'UNKNOWN TYPE OF CELL:',TYPE_OF_CELL
470                 WRITE(*,*)'ACCEPTABLE TYPES ARE:'
471                 WRITE(*,*)'SCALAR'
472                 WRITE(*,*)'U_MOMENTUM'
473                 WRITE(*,*)'V_MOMENTUM'
474                 WRITE(*,*)'W_MOMENTUM'
475                 CALL MFIX_EXIT(myPE)
476           END SELECT
477     
478           RETURN
479     
480           END SUBROUTINE STORE_CUT_FACE_INFO
481     
482     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
483     !                                                                      C
484     !  Module name: TEST_DEL_H                                             C
485     !  Purpose: tests the computation of wall distance                     C
486     !           If a negative distance is detected, the normal vector      C
487     !           is inverted                                                C
488     !                                                                      C
489     !  Author: Jeff Dietiker                              Date: 22-Feb-12  C
490     !  Reviewer:                                          Date:            C
491     !                                                                      C
492     !  Revision Number #                                  Date: ##-###-##  C
493     !  Author: #                                                           C
494     !  Purpose: #                                                          C
495     !                                                                      C
496     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
497       SUBROUTINE TEST_DEL_H(IJK,TYPE_OF_CELL)
498     
499           USE param
500           USE param1
501           USE parallel
502           USE constant
503           USE run
504           USE toleranc
505           USE geometry
506           USE indices
507           USE compar
508           USE sendrecv
509           USE quadric
510           USE cutcell
511     
512           IMPLICIT NONE
513           CHARACTER (LEN=*) :: TYPE_OF_CELL
514           DOUBLE PRECISION:: X_COPY,Y_COPY,Z_COPY
515           INTEGER :: IJK
516           INTEGER :: NODE,N_N1,N_N2,N
517           DOUBLE PRECISION :: Del_H
518           DOUBLE PRECISION :: Nx,Ny,Nz
519     
520           LOGICAL :: ALLOW_NEG_DIST = .TRUE.  ! forces GET_DEL_H_DES to output negative delh
521                                                ! i.e. do not let the subroutine overwrite negative values
522     
523     ! This subroutine tests values of del_H for nodes defining a cut cell.
524     ! Only nodes that are in the fluid region are tested.
525     ! Nodes belonging to the cut face should return zero (or near zero) values and are not tested.
526     ! If a negative del_H is detected, the unit normal vector is reversed.
527     
528           IF(NO_K) THEN
529              N_N1 = 5
530              N_N2 = 8
531           ELSE
532              N_N1 = 1
533              N_N2 = 8
534           ENDIF
535     
536           CALL GET_CELL_NODE_COORDINATES(IJK,TYPE_OF_CELL)
537     
538           DO NODE = 1,NUMBER_OF_NODES(IJK)
539              IF(CONNECTIVITY(IJK,NODE)<=IJKEND3) THEN  ! node does not belong to the cut-face
540                                                         ! i.e. is in the fluid region
541     
542                 DO N = N_N1,N_N2                     ! get node coordinate
543                    IF(CONNECTIVITY(IJK,NODE) == IJK_OF_NODE(N)) THEN
544                       X_COPY = X_NODE(N)
545                       Y_COPY = Y_NODE(N)
546                       Z_COPY = Z_NODE(N)
547                       EXIT
548                    ENDIF
549                 ENDDO
550     
551     !           Compute del_H
552                 CALL GET_DEL_H_DES(IJK,TYPE_OF_CELL,X_COPY,Y_COPY,Z_COPY,Del_H,Nx,Ny,Nz, ALLOW_NEG_DIST)
553     
554     
555                 IF(DEL_H<ZERO) THEN
556     
557     
558                    IF(PRINT_WARNINGS.AND.MyPE==PE_IO) THEN
559                       WRITE(*,*) ' Warning: Negative delh detected in scalar cell :',IJK
560                       WRITE(*,*) ' Location (X,Y,Z) = ',X_COPY,Y_COPY,Z_COPY
561                       WRITE(*,*) ' Reverting unit normal vector.'
562                    ENDIF
563     
564                    NORMAL_S(IJK,1) = -NORMAL_S(IJK,1)
565                    NORMAL_S(IJK,2) = -NORMAL_S(IJK,2)
566                    NORMAL_S(IJK,3) = -NORMAL_S(IJK,3)
567     
568                 ENDIF
569     
570              ENDIF
571           ENDDO
572     
573           RETURN
574     
575           END SUBROUTINE TEST_DEL_H
576     
577     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
578     !                                                                      C
579     !  Module name:  GET_DISTANCE_TO_WALL                                  C
580     !  Purpose: Finds the distance fraom any scalar cell to the closest    C
581     !  wall                                                                C
582     !                                                                      C
583     !  Author: Jeff Dietiker                              Date: 16-May-13  C
584     !  Reviewer:                                          Date:            C
585     !                                                                      C
586     !  Revision Number #                                  Date: ##-###-##  C
587     !  Author: #                                                           C
588     !  Purpose: #                                                          C
589     !                                                                      C
590     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
591       SUBROUTINE GET_DISTANCE_TO_WALL
592     
593           USE param
594           USE param1
595           USE parallel
596           USE constant
597           USE run
598           USE toleranc
599           USE geometry
600           USE indices
601           USE compar
602           USE sendrecv
603           USE quadric
604           USE cutcell
605           USE functions
606     
607           USE mpi_utility
608     
609           IMPLICIT NONE
610     !      CHARACTER (LEN=*) :: TYPE_OF_CELL
611           DOUBLE PRECISION:: X0,Y0,Z0,XREF,YREF,ZREF
612           INTEGER :: IJK,N
613     
614           DOUBLE PRECISION :: D_TO_CUT, D_TO_PE_REF
615     
616           INTEGER :: N_CUT_CELLS
617           INTEGER :: LIST_OF_CUT_CELLS(DIMENSION_3)
618     
619           INTEGER :: iproc,IERR,IJK_OFFSET,nb,n1,n2
620           INTEGER :: GLOBAL_N_CUT_CELLS
621           INTEGER, DIMENSION(0:numPEs-1) :: disp,rcount
622           LOGICAL, DIMENSION(0:numPEs-1) :: ALREADY_VISITED
623           DOUBLE PRECISION, DIMENSION(0:numPEs-1,3) :: PE_REFP,ALL_PE_REFP
624           DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: LOCAL_REFP_S,GLOBAL_REFP_S
625     
626           IF(MyPE==PE_IO) WRITE(*,*)'COMPUTING WALL DISTANCE...'
627     
628     ! Count the number of cut cells
629           N_CUT_CELLS = 0
630           DO IJK = IJKSTART3, IJKEND3
631              IF(INTERIOR_CELL_AT(IJK).AND.CUT_CELL_AT(IJK)) THEN
632                 N_CUT_CELLS = N_CUT_CELLS + 1
633                 LIST_OF_CUT_CELLS(N_CUT_CELLS) = IJK
634              ENDIF
635           ENDDO
636     
637     ! Store cut cell reference points (centroid of cut cell face)
638     ! and save a reference point for the entire processor
639           ALLOCATE (LOCAL_REFP_S(N_CUT_CELLS,3))
640     
641           N_CUT_CELLS = 0
642           PE_REFP(MyPE,1) = ZERO
643           PE_REFP(MyPE,2) = ZERO
644           PE_REFP(MyPE,3) = ZERO
645     
646           DO IJK = IJKSTART3, IJKEND3
647              IF(INTERIOR_CELL_AT(IJK).AND.CUT_CELL_AT(IJK)) THEN
648                 N_CUT_CELLS = N_CUT_CELLS + 1
649                 LOCAL_REFP_S(N_CUT_CELLS,1) = REFP_S(IJK,1)
650                 LOCAL_REFP_S(N_CUT_CELLS,2) = REFP_S(IJK,2)
651                 LOCAL_REFP_S(N_CUT_CELLS,3) = REFP_S(IJK,3)
652     
653                 PE_REFP(MyPE,1) = PE_REFP(MyPE,1) + REFP_S(IJK,1)
654                 PE_REFP(MyPE,2) = PE_REFP(MyPE,2) + REFP_S(IJK,2)
655                 PE_REFP(MyPE,3) = PE_REFP(MyPE,3) + REFP_S(IJK,3)
656              ENDIF
657           ENDDO
658     
659           IF(N_CUT_CELLS>0) THEN
660              PE_REFP(MyPE,1) = PE_REFP(MyPE,1) / N_CUT_CELLS
661              PE_REFP(MyPE,2) = PE_REFP(MyPE,2) / N_CUT_CELLS
662              PE_REFP(MyPE,3) = PE_REFP(MyPE,3) / N_CUT_CELLS
663           ELSE
664              PE_REFP(MyPE,1) = UNDEFINED
665              PE_REFP(MyPE,2) = UNDEFINED
666              PE_REFP(MyPE,3) = UNDEFINED
667           ENDIF
668     
669     !======================================================================
670     ! Now, gather the local reference points to head node
671     ! to get a global list, and broadcast it to each processor.
672     !======================================================================
673     
674     !======================================================================
675     ! First get the offset and build the rcount and disp arrays.
676     ! rcount is the number of elements to be gathered.
677     ! disp is the displacement of the variable size gather,
678     ! i.e. the cumulative sum at a given procesor.
679     !======================================================================
680     
681           CALL allgather_1i (N_CUT_CELLS,rcount,IERR)
682     
683           IF (myPE == 0) THEN
684              IJK_OFFSET = 0
685           ELSE
686              IJK_OFFSET = 0
687              DO iproc=0,myPE-1
688                 IJK_OFFSET = IJK_OFFSET + rcount(iproc)
689              ENDDO
690           ENDIF
691     
692           CALL allgather_1i (IJK_OFFSET,disp,IERR)
693     
694           CALL allgather_1d (PE_REFP(MyPE,1),ALL_PE_REFP(:,1),IERR)
695           CALL allgather_1d (PE_REFP(MyPE,2),ALL_PE_REFP(:,2),IERR)
696           CALL allgather_1d (PE_REFP(MyPE,3),ALL_PE_REFP(:,3),IERR)
697     
698     !======================================================================
699     ! Get the global number of cut cells and allocate the
700     ! global reference point array. Each processor gets its own
701     ! copy of all cut cell reference points !!
702     !======================================================================
703           CALL GLOBAL_ALL_SUM(N_CUT_CELLS, GLOBAL_N_CUT_CELLS,  PE_IO, ierr )
704           ALLOCATE (GLOBAL_REFP_S(GLOBAL_N_CUT_CELLS,3))
705     
706     !======================================================================
707     ! For a serial run, the global and local arrays are the same.
708     ! for a parallel run, first gather on head node, then broadcast to all.
709     !======================================================================
710     
711           IF(numPEs==1) THEN  ! Serial run
712              GLOBAL_REFP_S =  LOCAL_REFP_S
713           ELSE !Parallel run
714              call gatherv_1d( LOCAL_REFP_S(:,1), N_CUT_CELLS, GLOBAL_REFP_S(:,1), rcount, disp, PE_IO, ierr )
715              call gatherv_1d( LOCAL_REFP_S(:,2), N_CUT_CELLS, GLOBAL_REFP_S(:,2), rcount, disp, PE_IO, ierr )
716              call gatherv_1d( LOCAL_REFP_S(:,3), N_CUT_CELLS, GLOBAL_REFP_S(:,3), rcount, disp, PE_IO, ierr )
717     
718              call bcast(GLOBAL_REFP_S(:,1))
719              call bcast(GLOBAL_REFP_S(:,2))
720              call bcast(GLOBAL_REFP_S(:,3))
721           ENDIF
722     
723     
724           ALREADY_VISITED(:) = .FALSE.
725     !======================================================================
726     !  Loop through all scalar cells, grouped by processor.
727     !  Compute the distance to each cut cell and take the minimum
728     !  This is doe in three passes, in order of likelyhood to find a wall:
729     !  1) Within local processor.
730     !  2) Within neighboring processors (from send2 schedule)
731     !  3) Within all other processors
732     !
733     !  The idea is that it is very likely that cut cells located
734     !  on current or neighbor processors will contribute to the
735     !  wall distance calculation, whereas cut cells located
736     !  on the other processors will not contribute
737     !  During the third pass, the entire processor is skipped
738     !  if a reference point (average of all cut faces reference
739     !  points) is farther than the current wall distance.
740     !  This should speed up the brute calculation of going
741     !  explicitely through all cut cells.
742     !  However, it is not guaranteed that dwall will be
743     !  independent of the grid partitionning. This is considered
744     !  very unlikely at the moment
745     !  Change DWALL_BRUTE_FORCE to .TRUE. to force brute-force
746     !  calculation.
747     !======================================================================
748     
749           DO IJK = IJKSTART3, IJKEND3
750     
751              CALL WRITE_PROGRESS_BAR(IJK,IJKEND3 - IJKSTART3 + 1,'C')
752     
753              IF(INTERIOR_CELL_AT(IJK)) THEN
754                 DWALL(IJK) = UNDEFINED
755     
756                 IF(.NOT.BLOCKED_CELL_AT(IJK)) THEN
757     
758     ! Get coordinates of cell center
759                    CALL GET_CELL_NODE_COORDINATES(IJK,'SCALAR')
760     
761                    X0 = X_NODE(0)
762                    Y0 = Y_NODE(0)
763                    Z0 = Z_NODE(0)
764     
765     !======================================================================
766     ! First pass: Loop through local cut cells
767     !======================================================================
768     
769                    ALREADY_VISITED(MyPE) = .TRUE.
770     
771                    DO N = 1,N_CUT_CELLS
772     
773                       Xref = LOCAL_REFP_S(N,1)
774                       Yref = LOCAL_REFP_S(N,2)
775                       Zref = LOCAL_REFP_S(N,3)
776     
777                       D_TO_CUT = sqrt((X0 - Xref)**2 + (Y0 - Yref)**2 + (Z0 - Zref)**2)
778     
779                       DWALL(IJK) = MIN(DWALL(IJK),D_TO_CUT)
780     
781                    ENDDO
782     
783     
784     !======================================================================
785     ! Second pass: Loop through neighbor processors (use send2 schedule)
786     !======================================================================
787     
788                    DO nb=1,nsend2
789                       iproc = sendproc2(nb)
790                       ALREADY_VISITED(iproc) = .TRUE.
791                       n1 = disp(iproc)+1
792                       n2 = n1 + rcount(iproc) - 1
793     
794                       DO N = n1,n2
795     
796                          Xref = GLOBAL_REFP_S(N,1)
797                          Yref = GLOBAL_REFP_S(N,2)
798                          Zref = GLOBAL_REFP_S(N,3)
799     
800                          D_TO_CUT = sqrt((X0 - Xref)**2 + (Y0 - Yref)**2 + (Z0 - Zref)**2)
801     
802                          DWALL(IJK) = MIN(DWALL(IJK),D_TO_CUT)
803     
804                       ENDDO
805     
806                    ENDDO
807     
808     !======================================================================
809     ! Third pass: Loop through all other processors
810     ! skip already visited processors (myPE and its neigbhors)
811     ! First test if the average reference point is farther than dwall
812     ! if this is true, then skip the entire processor cut cells.
813     !======================================================================
814     
815     
816                    DO iproc=0,numPEs-1
817                       IF(ALREADY_VISITED(iproc)) CYCLE
818     
819                       Xref = ALL_PE_REFP(iproc,1)
820                       Yref = ALL_PE_REFP(iproc,2)
821                       Zref = ALL_PE_REFP(iproc,3)
822     
823                       D_TO_PE_REF = sqrt((X0 - Xref)**2 + (Y0 - Yref)**2 + (Z0 - Zref)**2)
824     
825                       IF((DWALL(IJK) < D_TO_PE_REF).AND. &
826                          (.NOT.DWALL_BRUTE_FORCE  )) THEN
827                          CYCLE
828                       ELSE
829                          n1 = disp(iproc)+1
830                          n2 = n1 + rcount(iproc) - 1
831     
832                          DO N = n1,n2
833     
834                             Xref = GLOBAL_REFP_S(N,1)
835                             Yref = GLOBAL_REFP_S(N,2)
836                             Zref = GLOBAL_REFP_S(N,3)
837     
838                             D_TO_CUT = sqrt((X0 - Xref)**2 + (Y0 - Yref)**2 + (Z0 - Zref)**2)
839     
840                             DWALL(IJK) = MIN(DWALL(IJK),D_TO_CUT)
841     
842                          ENDDO
843                       ENDIF
844                    ENDDO
845     
846                 ENDIF
847     
848              ENDIF
849     
850           ENDDO
851     
852           call send_recv(DWALL,2)
853     !======================================================================
854     !  Deallocate arrays before leaving
855     !======================================================================
856           DEALLOCATE (LOCAL_REFP_S)
857           DEALLOCATE (GLOBAL_REFP_S)
858           RETURN
859     
860           END SUBROUTINE GET_DISTANCE_TO_WALL
861