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