time_march_pic.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!  Subroutine: PIC_TIME_MARCH                                          !
!  Author: R. Garg                                                     !
!                                                                      !
!  Purpose: Main PIC driver routine.                                   !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE pic_time_march

! Global variables
!---------------------------------------------------------------------//
! Fluid time, simulation end time, time step size, number of time steps
      use run, only: time, tstop, dt, nstep
! Discrete particle time, time step size
      use discretelement, only: s_time, dtsolid
! MPPIC model step-size bounds
      use mfix_pic, only: dtpic_max, dtpic_cfl, dtpic_taup
! Local particle count
      use discretelement, only: pip
! Flag: Coupled fluid-solids simulation
      use discretelement, only: des_continuum_coupled
! Flag: Store _OLD arrays
      use discretelement, only: do_old
! Flag: Call user defined subroutines
      use run, only: call_usr
! Flag: Explicitly coupled gas-solids drag
      use discretelement, only: des_explicitly_coupled
! Number of mass outflows/inflows
      use pic_bc, only: pic_bcmo, pic_bcmi

! Module procedures
!---------------------------------------------------------------------//
      use desgrid, only: desgrid_pic
      use error_manager
      use mpi_funs_des, only: des_par_exchange
      use mpi_utility, only: global_all_sum
      use output_man, only: output_manager

      IMPLICIT NONE

! Local variables
!---------------------------------------------------------------------//
! time till which the PIC loop will be run
      double precision :: TEND_PIC_LOOP
! number of PIC time steps
      Integer :: PIC_ITERS
! Global number of parcels.
      INTEGER :: gPIP
!......................................................................!


! Set solids time to fluid time.
      s_time = time

      IF(des_continuum_coupled) THEN
         tend_pic_loop = time+dt
         dtsolid = min(dtpic_max, dt)
      ELSE
         tend_pic_loop = tstop
         dtsolid = dt
      ENDIF
      pic_iters = 0


      IF(call_usr) CALL usr0_des

! Compute the gas-phase pressure gradient
      IF(des_continuum_coupled) THEN
         IF(des_explicitly_coupled) CALL drag_gs_des1
         CALL calc_pg_grad
      ENDIF


! If the current time in the discrete loop exceeds the current time in
! the continuum simulation, exit the lagrangian loop
      DO WHILE(s_time.LT.tend_pic_loop)

         pic_iters  = pic_iters + 1

! Set the solids time step
!         DTSOLID = MERGE(MIN(DTPIC_MAX, DT), DTPIC_MAX,                &
!            DES_CONTINUUM_COUPLED)

! If next time step in the discrete loop will exceed the current time
! in the continuum simulation, modify the discrete time step so final
! time will match
         IF(s_time + dtsolid > tend_pic_loop) &
            dtsolid = tend_pic_loop - s_time

! Calculate the solids pressure
         CALL calc_ps_pic
         CALL calc_ps_grad_pic
         CALL interpolate_pic

         IF(des_continuum_coupled) CALL calc_drag_des

         IF (do_old) CALL cfupdateold

         CALL integrate_time_pic 

!         CALL WRITE_PARTICLE(6010)

! Apply mass outflow/inflow boundary conditions
         IF(pic_bcmo > 0) CALL mass_outflow_pic
         IF(pic_bcmi > 0) CALL mass_inflow_pic

! Impose the wall-particle boundary condition
         CALL apply_wall_bc_pic

! Exchange particle crossing processor boundaries
         CALL desgrid_pic(.true.)
         CALL des_par_exchange

         IF(s_time + dtsolid < tend_pic_loop .OR. &
            .NOT.des_explicitly_coupled ) THEN
! Bin particles to the fluid grid
            CALL particles_in_cell
! Calculate interpolation weights
            CALL calc_interp_weights
! Calculate mean fields
            CALL comp_mean_fields
         ENDIF

! This was moved from particles in cell and the passed variables should
! be added to particles in cell or made global.
         CALL report_stats_pic
! Update time to reflect changes
         s_time = s_time + dtsolid

         dtpic_max = min(dtpic_cfl, dtpic_taup)

! When coupled, all write calls are made in time_march (the continuum
! portion) according to user settings for spx_time and res_time.
! The following section targets data writes for DEM only cases:
         IF(.NOT.des_continuum_coupled) THEN
! Keep track of TIME for DEM simulations
            time = s_time
            nstep = nstep + 1
! Call the output manager to write RES and SPx data.
            CALL output_manager(.false., .false.)
         ENDIF  ! end if (.not.des_continuum_coupled)

      ENDDO

      CALL global_all_sum(pip, gpip)
      WRITE(err_msg, 3000) trim(ival(pic_iters)), trim(ival(gpip))
      CALL flush_err_msg(header=.false., footer=.false.)

 3000 FORMAT(/'PIC NITs: ',a,3x,'Total PIP: ', a)

      RETURN
      END SUBROUTINE pic_time_march




!         !DTPIC_MAX = MIN( 1e-04, DTPIC_MAX)
!         IF(MOD(PIC_ITERS, 10).eq.0) then
!            IF(DES_CONTINUUM_COUPLED) then
!               WRITE(ERR_MSG, 2000) DTSOLID, DTPIC_CFL, DTPIC_TAUP, DT
!            ELSE
!               WRITE(ERR_MSG, 2001) S_TIME, DTSOLID, DTPIC_CFL, DTPIC_TAUP, DT
!            ENDIF
!            CALL FLUSH_ERR_MSG(HEADER = .FALSE., FOOTER = .FALSE.)
!         ENDIF
!
! 2000 FORMAT(/5x,'DTSOLID CURRENT  = ',g17.8,/5x,'DTPIC_CFL',8x,'= ',  &
!         g17.8, /5x,'DTPIC TAUP',7x,'= ',g17.8,/5x,'DT FLOW',10x,'= ', &
!         g17.8)
!
! 2001 FORMAT(/5x,'TIME',13X,'= ',g17.8,/5x,'DTSOLID CURRENT  = ',g17.8,&
!         /5x,'DTPIC_CFL',8X,'= ', g17.8,/5x,'DTPIC TAUP',7x,'= ',g17.8,&
!         /5x,'DT FLOW',10X,'= ', g17.8)




      SUBROUTINE write_particle(NP)

      Use usr
      use compar
      use discretelement

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: NP

      INTEGER, SAVE :: CALLS = 0
      CHARACTER(len=128) :: FNAME

      fname=''; WRITE(fname, 2000) np, mype, calls
 2000 FORMAT('DBG/DBG_',i9.9,'_',i4.4,'_',i5.5,'.vtp')

      OPEN(unit=555, file=trim(fname), status='UNKNOWN')

      write(*,"('Saving: ',A,' at ',F15.8)") trim(fname), s_time


      WRITE(555, 3000)
 3000 FORMAT('<?xml version="1.0"?>')

      WRITE(555, 3001)
 3001 FORMAT('<VTKFile type="PolyData" ' &
         'version="0.1" byte_order="LittleEndian">')

      WRITE(555,"('<PolyData>')")

      WRITE(555, 3002)
 3002 FORMAT('<Piece NumberOfPoints="1" ',         &
         'NumberOfVerts="0" NumberOfLines="0" ', &
         'NumberOfStrips="0" ',                    &
         'NumberOfPolys="0">')

      WRITE(555,"('<Points>')")

 3003 FORMAT('<DataArray type="Float32" Name="Position" ', &
         'NumberOfComponents="3" format="ascii">')
      WRITE(555, 3003)
      WRITE(555,"(3(3x,F15.8))") des_pos_new(np,:)
      WRITE(555,"('</DataArray>')")

      WRITE(555,"('</Points>')")

 3004 FORMAT('<PointData Scalars="Diameter" Vectors="Velocity">')
      WRITE(555, 3004)

 3005 FORMAT('<DataArray type="Float32" ', &
         'Name="Diameter" format="ascii">')
      WRITE(555, 3005)
      WRITE(555,"(3x,F15.8)") des_radius(np)*2.0d0
      WRITE(555,"('</DataArray>')")

 3006 FORMAT('<DataArray type="Float32" Name="Velocity" ',&
         'NumberOfComponents="3" format="ascii">')
      WRITE(555, 3006)
      WRITE(555,"(3(3x,F15.8))") des_vel_new(np,:)
      WRITE(555,"('</DataArray>')")

      WRITE(555,"('</PointData>')")
      WRITE(555,"('<CellData></CellData>')")
      WRITE(555,"('<Verts></Verts>')")
      WRITE(555,"('<Lines></Lines>')")
      WRITE(555,"('<Strips></Strips>')")
      WRITE(555,"('<Polys></Polys>')")
      WRITE(555,"('</Piece>')")
      WRITE(555,"('</PolyData>')")
      WRITE(555,"('</VTKFile>')")

      close(555)

      calls = calls+1

      RETURN
      END SUBROUTINE write_particle