File: N:\mfix\model\des\pic\time_march_pic.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !  Subroutine: PIC_TIME_MARCH                                          !
3     !  Author: R. Garg                                                     !
4     !                                                                      !
5     !  Purpose: Main PIC driver routine.                                   !
6     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
7           SUBROUTINE PIC_TIME_MARCH
8     
9     ! Global variables
10     !---------------------------------------------------------------------//
11     ! Fluid time, simulation end time, time step size, number of time steps
12           use run, only: TIME, TSTOP, DT, NSTEP
13     ! Discrete particle time, time step size
14           use discretelement, only: S_TIME, DTSOLID
15     ! MPPIC model step-size bounds
16           use mfix_pic, only: DTPIC_MAX, DTPIC_CFL, DTPIC_TAUP
17     ! Local particle count
18           use discretelement, only: PIP
19     ! Flag: Coupled fluid-solids simulation
20           use discretelement, only: DES_CONTINUUM_COUPLED
21     ! Flag: Store _OLD arrays
22           use discretelement, only: DO_OLD
23     ! Flag: Call user defined subroutines
24           use run, only: CALL_USR
25     ! Flag: Explicitly coupled gas-solids drag
26           use discretelement, only: DES_EXPLICITLY_COUPLED
27     ! Number of mass outflows/inflows
28           use pic_bc, only: PIC_BCMO, PIC_BCMI
29     
30     ! Module procedures
31     !---------------------------------------------------------------------//
32           use desgrid, only: DESGRID_PIC
33           use error_manager
34           use mpi_funs_des, only: DES_PAR_EXCHANGE
35           use mpi_utility, only: GLOBAL_ALL_SUM
36           use output_man, only: output_manager
37     
38           IMPLICIT NONE
39     
40     ! Local variables
41     !---------------------------------------------------------------------//
42     ! time till which the PIC loop will be run
43           double precision :: TEND_PIC_LOOP
44     ! number of PIC time steps
45           Integer :: PIC_ITERS
46     ! Global number of parcels.
47           INTEGER :: gPIP
48     !......................................................................!
49     
50     
51     ! Set solids time to fluid time.
52           S_TIME = TIME
53     
54           IF(DES_CONTINUUM_COUPLED) THEN
55              TEND_PIC_LOOP = TIME+DT
56              DTSOLID = min(DTPIC_MAX, DT)
57           ELSE
58              TEND_PIC_LOOP = TSTOP
59              DTSOLID = DT
60           ENDIF
61           PIC_ITERS = 0
62     
63     
64           IF(CALL_USR) CALL USR0_DES
65     
66     ! Compute the gas-phase pressure gradient
67           IF(DES_CONTINUUM_COUPLED) THEN
68              IF(DES_EXPLICITLY_COUPLED) CALL DRAG_GS_DES1
69              CALL CALC_PG_GRAD
70           ENDIF
71     
72     
73     ! If the current time in the discrete loop exceeds the current time in
74     ! the continuum simulation, exit the lagrangian loop
75           DO WHILE(S_TIME.LT.TEND_PIC_LOOP)
76     
77              PIC_ITERS  = PIC_ITERS + 1
78     
79     ! Set the solids time step
80     !         DTSOLID = MERGE(MIN(DTPIC_MAX, DT), DTPIC_MAX,                &
81     !            DES_CONTINUUM_COUPLED)
82     
83     ! If next time step in the discrete loop will exceed the current time
84     ! in the continuum simulation, modify the discrete time step so final
85     ! time will match
86              IF(S_TIME + DTSOLID > TEND_PIC_LOOP) &
87                 DTSOLID = TEND_PIC_LOOP - S_TIME
88     
89     ! Calculate the solids pressure
90              CALL CALC_PS_PIC
91              CALL CALC_PS_GRAD_PIC
92              CALL INTERPOLATE_PIC
93     
94              IF(DES_CONTINUUM_COUPLED) CALL CALC_DRAG_DES
95     
96              IF (DO_OLD) CALL CFUPDATEOLD
97     
98              CALL INTEGRATE_TIME_PIC 
99     
100     !         CALL WRITE_PARTICLE(6010)
101     
102     ! Apply mass outflow/inflow boundary conditions
103              IF(PIC_BCMO > 0) CALL MASS_OUTFLOW_PIC
104              IF(PIC_BCMI > 0) CALL MASS_INFLOW_PIC
105     
106     ! Impose the wall-particle boundary condition
107              CALL APPLY_WALL_BC_PIC
108     
109     ! Exchange particle crossing processor boundaries
110              CALL DESGRID_PIC(.TRUE.)
111              CALL DES_PAR_EXCHANGE
112     
113              IF(S_TIME + DTSOLID < TEND_PIC_LOOP .OR. &
114                 .NOT.DES_EXPLICITLY_COUPLED ) THEN
115     ! Bin particles to the fluid grid
116                 CALL PARTICLES_IN_CELL
117     ! Calculate interpolation weights
118                 CALL CALC_INTERP_WEIGHTS
119     ! Calculate mean fields
120                 CALL COMP_MEAN_FIELDS
121              ENDIF
122     
123     ! This was moved from particles in cell and the passed variables should
124     ! be added to particles in cell or made global.
125              CALL REPORT_STATS_PIC
126     ! Update time to reflect changes
127              S_TIME = S_TIME + DTSOLID
128     
129              DTPIC_MAX = MIN(DTPIC_CFL, DTPIC_TAUP)
130     
131     ! When coupled, all write calls are made in time_march (the continuum
132     ! portion) according to user settings for spx_time and res_time.
133     ! The following section targets data writes for DEM only cases:
134              IF(.NOT.DES_CONTINUUM_COUPLED) THEN
135     ! Keep track of TIME for DEM simulations
136                 TIME = S_TIME
137                 NSTEP = NSTEP + 1
138     ! Call the output manager to write RES and SPx data.
139                 CALL OUTPUT_MANAGER(.FALSE., .FALSE.)
140              ENDIF  ! end if (.not.des_continuum_coupled)
141     
142           ENDDO
143     
144           CALL GLOBAL_ALL_SUM(PIP, gPIP)
145           WRITE(ERR_MSG, 3000) trim(iVal(PIC_ITERS)), trim(iVal(gPIP))
146           CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
147     
148      3000 FORMAT(/'PIC NITs: ',A,3x,'Total PIP: ', A)
149     
150           RETURN
151           END SUBROUTINE PIC_TIME_MARCH
152     
153     
154     
155     
156     !         !DTPIC_MAX = MIN( 1e-04, DTPIC_MAX)
157     !         IF(MOD(PIC_ITERS, 10).eq.0) then
158     !            IF(DES_CONTINUUM_COUPLED) then
159     !               WRITE(ERR_MSG, 2000) DTSOLID, DTPIC_CFL, DTPIC_TAUP, DT
160     !            ELSE
161     !               WRITE(ERR_MSG, 2001) S_TIME, DTSOLID, DTPIC_CFL, DTPIC_TAUP, DT
162     !            ENDIF
163     !            CALL FLUSH_ERR_MSG(HEADER = .FALSE., FOOTER = .FALSE.)
164     !         ENDIF
165     !
166     ! 2000 FORMAT(/5x,'DTSOLID CURRENT  = ',g17.8,/5x,'DTPIC_CFL',8x,'= ',  &
167     !         g17.8, /5x,'DTPIC TAUP',7x,'= ',g17.8,/5x,'DT FLOW',10x,'= ', &
168     !         g17.8)
169     !
170     ! 2001 FORMAT(/5x,'TIME',13X,'= ',g17.8,/5x,'DTSOLID CURRENT  = ',g17.8,&
171     !         /5x,'DTPIC_CFL',8X,'= ', g17.8,/5x,'DTPIC TAUP',7x,'= ',g17.8,&
172     !         /5x,'DT FLOW',10X,'= ', g17.8)
173     
174     
175     
176     
177           SUBROUTINE WRITE_PARTICLE(NP)
178     
179           Use usr
180           use compar
181           use discretelement
182     
183           IMPLICIT NONE
184     
185           INTEGER, INTENT(IN) :: NP
186     
187           INTEGER, SAVE :: CALLS = 0
188           CHARACTER(len=128) :: FNAME
189     
190           FNAME=''; WRITE(FNAME, 2000) NP, myPE, CALLS
191      2000 FORMAT('DBG/DBG_',I9.9,'_',I4.4,'_',I5.5,'.vtp')
192     
193           OPEN(UNIT=555, FILE=trim(FNAME), STATUS='UNKNOWN')
194     
195           write(*,"('Saving: ',A,' at ',F15.8)") trim(FNAME), S_TIME
196     
197     
198           WRITE(555, 3000)
199      3000 FORMAT('<?xml version="1.0"?>')
200     
201           WRITE(555, 3001)
202      3001 FORMAT('<VTKFile type="PolyData" ' &
203              'version="0.1" byte_order="LittleEndian">')
204     
205           WRITE(555,"('<PolyData>')")
206     
207           WRITE(555, 3002)
208      3002 FORMAT('<Piece NumberOfPoints="1" ',         &
209              'NumberOfVerts="0" NumberOfLines="0" ', &
210              'NumberOfStrips="0" ',                    &
211              'NumberOfPolys="0">')
212     
213           WRITE(555,"('<Points>')")
214     
215      3003 FORMAT('<DataArray type="Float32" Name="Position" ', &
216              'NumberOfComponents="3" format="ascii">')
217           WRITE(555, 3003)
218           WRITE(555,"(3(3x,F15.8))") DES_POS_NEW(NP,:)
219           WRITE(555,"('</DataArray>')")
220     
221           WRITE(555,"('</Points>')")
222     
223      3004 FORMAT('<PointData Scalars="Diameter" Vectors="Velocity">')
224           WRITE(555, 3004)
225     
226      3005 FORMAT('<DataArray type="Float32" ', &
227              'Name="Diameter" format="ascii">')
228           WRITE(555, 3005)
229           WRITE(555,"(3x,F15.8)") DES_RADIUS(NP)*2.0d0
230           WRITE(555,"('</DataArray>')")
231     
232      3006 FORMAT('<DataArray type="Float32" Name="Velocity" ',&
233              'NumberOfComponents="3" format="ascii">')
234           WRITE(555, 3006)
235           WRITE(555,"(3(3x,F15.8))") DES_VEL_NEW(NP,:)
236           WRITE(555,"('</DataArray>')")
237     
238           WRITE(555,"('</PointData>')")
239           WRITE(555,"('<CellData></CellData>')")
240           WRITE(555,"('<Verts></Verts>')")
241           WRITE(555,"('<Lines></Lines>')")
242           WRITE(555,"('<Strips></Strips>')")
243           WRITE(555,"('<Polys></Polys>')")
244           WRITE(555,"('</Piece>')")
245           WRITE(555,"('</PolyData>')")
246           WRITE(555,"('</VTKFile>')")
247     
248           close(555)
249     
250           CALLS = CALLS+1
251     
252           RETURN
253           END SUBROUTINE WRITE_PARTICLE
254     
255