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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: CHECK_DATA_09                                          C
4     !  Purpose: Check point source specifications.                         C
5     !                                                                      C
6     !  Author: J. Musser                                  Date: 10-JUN-13  C
7     !                                                                      C
8     !  Literature/Document References:                                     C
9     !                                                                      C
10     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
11           SUBROUTINE SET_PS
12     
13           use param
14           use run
15           use physprop
16           use ps
17           use compar
18           use geometry
19           use mpi_utility
20           use functions
21     
22           implicit none
23     
24           INTEGER :: IJK, I, J, K, M, N
25     
26           INTEGER PSV
27     
28           CHARACTER(LEN=64) :: eMsg
29     
30           INTEGER :: PS_SIZE
31     
32           DOUBLE PRECISION, allocatable :: lData_dp(:)
33           DOUBLE PRECISION, allocatable :: gData_dp(:)
34     
35           LOGICAL, parameter :: dbg_PS = .FALSE.
36     
37           if(.NOT.POINT_SOURCE) return
38     
39     ! DETERMINE WHICH BOUNDARY CONDITION INDICES HAVE VALUES
40           L50: do PSV = 1, DIMENSION_PS
41     
42              IF(.NOT.PS_DEFINED(PSV)) cycle L50
43     
44     
45     ! Calculate the velocity magnitude and normalize the axial components.
46              CALL CALC_PS_VEL_MAG(PS_VEL_MAG_g(PSV), PS_U_g(PSV),          &
47                 PS_V_g(PSV), PS_W_g(PSV))
48     
49              CALL CALC_PS_CpxMFLOW(PS_CpxMFLOW_g(PSV), PS_MASSFLOW_g(PSV), &
50                 PS_T_g(PSV), PS_X_g(PSV,:), 0, C_PG0, DIM_N_g, MW_g)
51     
52              do M=1, MMAX
53                 CALL CALC_PS_VEL_MAG(PS_VEL_MAG_s(PSV,M), PS_U_s(PSV,M),   &
54                    PS_V_s(PSV,M), PS_W_s(PSV,M))
55     
56                 CALL CALC_PS_CpxMFLOW(PS_CpxMFLOW_s(PSV,M),                &
57                    PS_MASSFLOW_s(PSV,M), PS_T_s(PSV,M), PS_X_s(PSV,M,:), M,&
58                    C_PS0(M), DIM_N_s, MW_s(M,:))
59              enddo
60     
61     
62     ! Calculate the number of cells comprising the point source. This
63     ! information is used to allocate some temp arrays.
64     !---------------------------------------------------------------------//
65              PS_SIZE = (PS_I_E(PSV) - PS_I_W(PSV) + 1) * &
66                        (PS_J_N(PSV) - PS_J_S(PSV) + 1)
67              if(DO_K) PS_SIZE = PS_SIZE * (PS_K_T(PSV) - PS_K_B(PSV) + 1)
68     
69              if(PS_SIZE < 1) then
70                  eMsg = ''; write(eMsg,"('Invalid PS size: ', I4)")PS_SIZE
71                  goto 500
72              endif
73     
74     
75              allocate(lData_dp( 0:numPEs-1)); lData_dp = ZERO
76              allocate(gData_dp( 0:numPEs-1)); gData_dp = ZERO
77     
78     
79     ! Calculate the volume of the PointSource cells
80     !---------------------------------------------------------------------//
81     ! Initialize the loop counter
82              PS_VOLUME(PSV) = ZERO
83     
84              do k = PS_K_B(PSV), PS_K_T(PSV)
85              do j = PS_J_S(PSV), PS_J_N(PSV)
86              do i = PS_I_W(PSV), PS_I_E(PSV)
87     
88                 if(IS_ON_myPE_owns(I, J, K)) then
89                 IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
90                    ijk = funijk(i,j,k)
91     
92                    if(fluid_at(ijk)) &
93                       lData_dp(myPE) = lData_dp(myPE) + VOL(IJK)
94                 endif
95     
96              enddo
97              enddo
98              enddo
99     
100     ! Each process (in DMP) only knows about the volume of the point source
101     ! it sees. Invoke send/recv so that all process can calculate the total
102     ! volume.
103              CALL global_all_sum(lData_dp, gData_dp)
104     
105              PS_VOLUME(PSV) = sum(gData_dp)
106              if(PS_VOLUME(PSV) == ZERO) then
107                 eMsg = 'No PS_VOLUME == ZERO'
108                 CALL DEBUG_PS(PSV, PS_SIZE)
109                 goto 501
110              endif
111     
112              if(allocated(lData_dp)) deallocate(lData_dp)
113              if(allocated(gData_dp)) deallocate(gData_dp)
114     
115     
116              IF(dbg_PS) CALL DEBUG_PS(PSV, PS_SIZE)
117     
118     
119           enddo L50
120     
121           return
122     
123       500 continue
124           if(myPE == PE_IO) &
125              write(*,"('PointSource setup Error: ',A)") trim(eMsg)
126     
127           call mfix_exit(myPE)
128     
129     
130       501 continue
131           write(*,"('PointSource setup Error: ',I3,2x,A)") myPE, trim(eMsg)
132     
133           call mfix_exit(myPE)
134     
135     
136           RETURN
137     
138           contains
139     
140     
141     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
142     !                                                                      C
143     !  Module name: CALC_PS_CpxMFLOW                                       C
144     !  Purpose: Calculate velocity magnitude and normalize U/V/W.          C
145     !                                                                      C
146     !  Author: J. Musser                                  Date: 10-JUN-13  C
147     !                                                                      C
148     !  Literature/Document References:                                     C
149     !                                                                      C
150     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
151           SUBROUTINE CALC_PS_CpxMFLOW(CpxMFLOW, PS_MFLOW, PS_T, PS_X, lM, &
152              Cp0, lDIM_N, lMW)
153     
154           use constant, only: GAS_CONST_cal
155           use read_thermochemical, only: calc_CpoR
156     
157           DOUBLE PRECISION, intent(out) :: CpxMFLOW
158     
159           INTEGER, intent(in) :: lM
160           INTEGER, intent(in) :: lDIM_N
161     
162           DOUBLE PRECISION, intent(in) :: Cp0
163           DOUBLE PRECISION, intent(in) :: PS_MFLOW
164           DOUBLE PRECISION, intent(in) :: PS_T          ! Temperature
165           DOUBLE PRECISION, intent(in) :: PS_X(lDIM_N)  ! Species Mass Frac
166     
167           DOUBLE PRECISION, intent(in) :: lMW(lDIM_N)
168     
169           INTEGER :: IER
170     
171     ! If there is no mass flow for this phase, then there is no need to
172     ! calculate a CPxMFLUX. Set it to zero and return.
173           if(.NOT.ENERGY_EQ .OR. PS_MFLOW == ZERO) then
174              CpxMFLOW = ZERO
175              return
176           endif
177     
178     ! Calculate the average specific heat.
179           if(Cp0 == UNDEFINED) then
180              IF(.NOT.database_read) call read_database0(IER)
181              CpxMFLOW = ZERO
182              do N = 1, NMAX(lM)
183                 CpxMFLOW = CpxMFLOW + PS_X(N) * (GAS_CONST_cal / lMW(N)) * &
184                   calc_CpoR(PS_T, lM, N, IER)
185              enddo
186           else
187              CpxMFLOW = Cp0
188           endif
189     
190     ! Unit conversion (if needed)
191           if(UNITS == 'SI') CpxMFLOW = 4.183925d3*CpxMFLOW
192           CpxMFLOW = CpxMFLOW*PS_MFLOW
193     
194           END SUBROUTINE CALC_PS_CpxMFLOW
195     
196     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
197     !                                                                      C
198     !  Module name: CALC_PS_VEL_MAG                                        C
199     !  Purpose: Calculate velocity magnitude and normalize U/V/W.          C
200     !                                                                      C
201     !  Author: J. Musser                                  Date: 10-JUN-13  C
202     !                                                                      C
203     !  Literature/Document References:                                     C
204     !                                                                      C
205     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
206           SUBROUTINE CALC_PS_VEL_MAG(VEL_MAG, lU, lV, lW)
207     
208           DOUBLE PRECISION, intent(inout) :: VEL_MAG
209           DOUBLE PRECISION, intent(inout) :: lU, lV, lW
210     
211     ! Normalize velocities:
212           VEL_MAG = lU**2 + lV**2
213           if(DO_K) VEL_MAG = VEL_MAG + lW**2
214     
215           VEL_MAG = sqrt(VEL_MAG)
216     
217           if(VEL_MAG > small_number) then
218              lU = lU/VEL_MAG
219              lV = lV/VEL_MAG
220              lW = lW/VEL_MAG
221           else
222              VEL_MAG = ZERO
223              lU = ZERO
224              lV = ZERO
225              lW = ZERO
226           endif
227     
228     
229           END SUBROUTINE CALC_PS_VEL_MAG
230     
231     
232           END SUBROUTINE SET_PS
233     
234     
235     
236     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
237     !                                                                      C
238     !  Module name: DEBUG_PS                                               C
239     !                                                                      C
240     !  Purpose: Write out some information about that point source that    C
241     !  may be useful in debugging.                                         C
242     !                                                                      C
243     !  Author: J. Musser                                  Date: 24-JUN-13  C
244     !                                                                      C
245     !  Literature/Document References:                                     C
246     !                                                                      C
247     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
248           SUBROUTINE DEBUG_PS(lPSV, lPS_SIZE)
249     
250           use bc
251           use compar
252           use constant
253           use cutcell
254           use fldvar
255           use geometry
256           use ic
257           use indices
258           use mflux
259           use parallel
260           use param
261           use param1
262           use physprop
263           use run
264           use sendrecv
265           use toleranc
266           use usr
267           use rxns
268           use ps
269           use mpi_utility
270           use functions
271     
272           IMPLICIT NONE
273     
274     ! Index of PS source to debug.
275           INTEGER, intent(in) :: lPSV
276     ! Number of cells comprising the point source.
277           INTEGER, intent(in) :: lPS_SIZE
278     
279           INTEGER :: IJK, I, J, K, M, N
280     
281           INTEGER :: lc1
282     
283           INTEGER, allocatable :: lFlags_i(:,:)
284           INTEGER, allocatable :: gFlags_i(:,:)
285     
286           if(myPE == PE_IO) then
287              write(*,"(3/,3x,'Debug Point Source Index: ',I3)") lPSV
288              write(*,"(/3x,'Size: ',I4)") lPS_SIZE
289           endif
290     
291           allocate(lFlags_i(lPS_SIZE,1:2) );   lFlags_i = 0
292           allocate(gFlags_i(lPS_SIZE,1:2) );   gFlags_i = 0
293     
294           lc1 = 0
295     
296           do k = PS_K_B(lPSV), PS_K_T(lPSV)
297           do j = PS_J_S(lPSV), PS_J_N(lPSV)
298           do i = PS_I_W(lPSV), PS_I_E(lPSV)
299     
300              lc1 = lc1 + 1
301              if(IS_ON_myPE_owns(I, J, K)) then
302              IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
303                 ijk = funijk(i,j,k)
304                 if(fluid_at(ijk)) then
305                    lFlags_i(lc1,1) = myPE
306                    lFlags_i(lc1,2) = FLAG(IJK)
307                 endif
308              endif
309           enddo
310           enddo
311           enddo
312     
313     ! Collect flag information on root.
314           CALL global_sum(lFlags_i, gFlags_i)
315     
316     ! Write some information to the screen.
317           if(myPE == PE_IO) then
318              write(*,"(/5x,'Location:')")
319              write(*,"( 5x,'X:',2(2x,g12.5),' :: ',2(2x,I4))")&
320                 PS_X_w(lPSV), PS_X_e(lPSV), PS_I_w(lPSV), PS_I_e(lPSV)
321              write(*,"( 5x,'Y:',2(2x,g12.5),' :: ',2(2x,I4))")&
322                 PS_Y_s(lPSV), PS_Y_n(lPSV), PS_J_s(lPSV), PS_J_n(lPSV)
323              if(DO_K)write(*,"( 5x,'Z:',2(2x,g12.5),' :: ',2(2x,I4))")&
324                 PS_Z_b(lPSV), PS_Z_t(lPSV), PS_K_b(lPSV), PS_K_t(lPSV)
325     
326              write(*,"(/5x,'Volume: ',g12.5)") PS_VOLUME(lPSV)
327     
328     
329              if(PS_MASSFLOW_G(lPSV) > small_number) then
330                 write(*,"(//5x,'Point Source Gas Phase:')")
331                 write(*,"(7x,'Mass Flow Rate: ',g12.5)")PS_MASSFLOW_G(lPSV)
332                 write(*,"(7x,'Velocity Magnitude: ',g12.5)") PS_VEL_MAG_g(lPSV)
333                 write(*,"(7x,'Normal:')")
334                 write(*,"(9x,'x-Axis: ',g12.5)")PS_U_g(lPSV)
335                 write(*,"(9x,'y-Axis: ',g12.5)")PS_V_g(lPSV)
336                 if(DO_K) write(*,"(9x,'z-Axis: ',g12.5)")PS_W_g(lPSV)
337                 if(energy_eq) &
338                    write(*,"(7x,'Temperature: ',g12.5)")PS_T_g(lPSV)
339                 if(species_eq(0)) then
340                    write(*,"(7x,'Species Composition:')")
341                    do n=1,nmax(0)
342                       write(*,"(9x,A,': ',g12.5)") &
343                          trim(SPECIES_ALIAS_g(n)), PS_X_g(lPSV,N)
344                    enddo
345                 endif
346              else
347                 write(*,"(//5x,'No gas phase point source.')")
348              endif
349     
350     
351              do m=1,mmax
352                 if(PS_MASSFLOW_S(lPSV,M) > small_number) then
353                    write(*,"(//5x,'Point Source Solids Phase ',I1,':')")M
354                    write(*,"(7x,'Mass Flow Rate: ',g12.5)")PS_MASSFLOW_S(lPSV,M)
355                    write(*,"(7x,'Velocity Magnitude: ',g12.5)") PS_VEL_MAG_s(lPSV,M)
356                    write(*,"(7x,'Normal:')")
357                    write(*,"(9x,'x-Axis: ',g12.5)")PS_U_s(lPSV,M)
358                    write(*,"(9x,'y-Axis: ',g12.5)")PS_V_s(lPSV,M)
359                    if(DO_K) write(*,"(9x,'z-Axis: ',g12.5)")PS_W_s(lPSV,M)
360                    if(energy_eq) &
361                       write(*,"(7x,'Temperature: ',g12.5)")PS_T_s(lPSV,M)
362                    if(species_eq(m)) then
363                       write(*,"(7x,'Species Composition:')")
364                       do n=1,nmax(m)
365                          write(*,"(9x,A,': ',g12.5)") &
366                             trim(SPECIES_ALIAS_s(m,n)), PS_X_s(lPSV,m,N)
367                       enddo
368                    endif
369                 else
370                    write(*,"(//5x,'No solids phase ',I1,' point source.')") m
371                 endif
372              enddo
373     
374              write(*,"(//5x,'Point Source Cells:')")
375              write(*,"(9x,'IJK',3(6x,A1),3x,'OWNS',3x,'FLAG')") 'I','J','K'
376     
377              lc1 = 0
378              do k = PS_K_B(lPSV), PS_K_T(lPSV)
379              do j = PS_J_S(lPSV), PS_J_N(lPSV)
380              do i = PS_I_W(lPSV), PS_I_E(lPSV)
381                 IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
382                 lc1 = lc1 + 1
383                 write(*,"(4x,I8,5(3x,I4))") IJK, I, J, K,  gFlags_i(lc1,:)
384              enddo
385              enddo
386              enddo
387     
388           endif
389     
390           if(allocated(lFlags_i)) deallocate(lFlags_i)
391           if(allocated(gFlags_i)) deallocate(gFlags_i)
392     
393           END SUBROUTINE DEBUG_PS
394     
395