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