cfnewvalues.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: CFNEWVALUES                                            C
!
!  Purpose: DES - Calculate the new values of particle velocity,
!           position, angular velocity etc
!
!                                                                      C
!  Comments: Implements Eqns 1, 2, 3, 4 & 5  from the following paper:
!    Tsuji Y., Kawaguchi T., and Tanak T., "Lagrangian numerical
!    simulation of plug glow of cohesionless particles in a
!    horizontal pipe", Powder technology, 71, 239-250, 1992
!
!  pradeep : changes for parallel processing
!          1. periodic boundaries might lie in different proc. so adjust
!             particle position for periodic removed
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE cfnewvalues

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE constant
      USE des_bc
      USE derived_types, only: multisap, boxhandle
      USE discretelement
      USE fldvar
      USE mfix_pic
      USE mpi_utility
      USE parallel
      USE run
      USE param
      USE param1
      USE physprop
      use geometry, only: do_k, no_k
      use multi_sweep_and_prune, only: aabb_t, multisap_sort, multisap_update

      IMPLICIT NONE
!-----------------------------------------------
! Local Variables
!-----------------------------------------------
      INTEGER :: L
      LOGICAL, SAVE :: FIRST_PASS = .true.
      DOUBLE PRECISION :: OMEGA_MAG,OMEGA_UNIT(3),ROT_ANGLE

      type(aabb_t) aabb

!-----------------------------------------------

! Adams-Bashforth defaults to Euler for the first time step.
      IF(first_pass .AND. intg_adams_bashforth) THEN
         DO l =1, max_pip
            IF(is_nonexistent(l)) cycle                       ! Only real particles
            IF(is_entering(l).or.is_entering_ghost(l)) cycle  ! Only non-entering
            IF(is_ghost(l)) cycle                             ! Skip ghost particles
            des_acc_old(l,:) = fc(l,:)/pmass(l) + grav(:)
            rot_acc_old(l,:) = tow(l,:)
         ENDDO
      ENDIF


!!$omp parallel default(none)                    &
!!$omp shared(MAX_PIP,INTG_EULER,INTG_ADAMS_BASHFORTH,fc,tow,do_nsearch,   &
!!$omp       omega_new,omega_old,pmass,grav,des_vel_new,des_pos_new,       &
!!$omp       des_vel_old,des_pos_old,dtsolid,omoi,des_acc_old,rot_acc_old, &
!!$omp       ppos,neighbor_search_rad_ratio,des_radius,DO_OLD, iGlobal_ID, &
!!$omp       particle_orientation,orientation,boxhandle,multisap,particle_state) &
!!$omp private(l,rot_angle,omega_mag,omega_unit,aabb)

! If a particle is classified as new, then forces are ignored.
! Classification from new to existing is performed in routine
! des_check_new_particle.f

! Advance particle position, velocity
! first-order method
      IF (intg_euler) THEN
!!$omp sections
!!$omp section
         WHERE(particle_state == normal_particle) des_vel_new(:,1) =   &
            des_vel_new(:,1) + dtsolid*(fc(:,1)/pmass(:) + grav(1))
         WHERE(particle_state < normal_ghost) des_pos_new(:,1) =       &
            des_pos_new(:,1) + des_vel_new(:,1)*dtsolid
         fc(:,1) = zero

!!$omp section
         WHERE(particle_state == normal_particle) des_vel_new(:,2) =   &
            des_vel_new(:,2) + dtsolid*(fc(:,2)/pmass(:) + grav(2))
         WHERE(particle_state < normal_ghost) des_pos_new(:,2) =       &
            des_pos_new(:,2) + des_vel_new(:,2)*dtsolid
         fc(:,2) = zero

!!$omp section
         WHERE(particle_state == normal_particle) des_vel_new(:,3) =   &
            des_vel_new(:,3) + dtsolid*(fc(:,3)/pmass(:) + grav(3))
         WHERE(particle_state < normal_ghost) des_pos_new(:,3) =       &
            des_pos_new(:,3) + des_vel_new(:,3)*dtsolid
         fc(:,3) = zero

!!$omp section
         WHERE(particle_state == normal_particle) omega_new(:,1) =     &
            omega_new(:,1) + tow(:,1)*omoi(:)*dtsolid
         tow(:,1) = zero

!!$omp section
         WHERE(particle_state == normal_particle)omega_new(:,2) =      &
            omega_new(:,2) + tow(:,2)*omoi(:)*dtsolid
         tow(:,2) = zero

!!$omp section
         WHERE(particle_state == normal_particle)omega_new(:,3) =      &
            omega_new(:,3) + tow(:,3)*omoi(:)*dtsolid
         tow(:,3) = zero
!!$omp end sections

! Second-order Adams-Bashforth/Trapezoidal scheme
      ELSEIF (intg_adams_bashforth) THEN

!!$omp sections
!!$omp section
         WHERE(particle_state == normal_particle) &
            fc(:max_pip,1) = fc(:max_pip,1)/pmass(:max_pip) + grav(1)
         WHERE(particle_state == normal_particle) &
            des_vel_new(:max_pip,1) = des_vel_old(:max_pip,1) + 0.5d0* &
               (3.d0*fc(:max_pip,1)-des_acc_old(:max_pip,1) )*dtsolid
            des_acc_old(:max_pip,1) = fc(:max_pip,1)

         WHERE(particle_state < normal_ghost)                          &
            des_pos_new(:max_pip,1) = des_pos_old(:max_pip,1) + 0.5d0* &
               (des_vel_old(:max_pip,1)+des_vel_new(:max_pip,1))*dtsolid
         fc(:max_pip,1) = zero

!!$omp section
         WHERE(particle_state == normal_particle) &
            fc(:max_pip,2) = fc(:max_pip,2)/pmass(:max_pip) + grav(2)
         WHERE(particle_state == normal_particle) &
            des_vel_new(:max_pip,2) = des_vel_old(:max_pip,2) + 0.5d0* &
               (3.d0*fc(:max_pip,2)-des_acc_old(:max_pip,2) )*dtsolid
            des_acc_old(:max_pip,2) = fc(:max_pip,2)

         WHERE(particle_state < normal_ghost)                          &
            des_pos_new(:max_pip,2) = des_pos_old(:max_pip,2) + 0.5d0* &
               (des_vel_old(:max_pip,2)+des_vel_new(:max_pip,2))*dtsolid
         fc(:max_pip,2) = zero

!!$omp section
         WHERE(particle_state == normal_particle) &
            fc(:max_pip,3) = fc(:max_pip,3)/pmass(:max_pip) + grav(3)
         WHERE(particle_state == normal_particle) &
            des_vel_new(:max_pip,3) = des_vel_old(:max_pip,3) + 0.5d0* &
                 (3.d0*fc(:max_pip,3)-des_acc_old(:max_pip,3) )*dtsolid
            des_acc_old(:max_pip,3) = fc(:max_pip,3)

         WHERE(particle_state < normal_ghost)                          &
            des_pos_new(:max_pip,3) = des_pos_old(:max_pip,3) + 0.5d0* &
               (des_vel_old(:max_pip,3)+des_vel_new(:max_pip,3))*dtsolid
         fc(:max_pip,3) = zero

!!$omp section
        WHERE(particle_state(:max_pip) == normal_particle) &
           omega_new(:max_pip,1) = omega_old(:max_pip,1) + 0.5d0*     &
               (3.d0*tow(:max_pip,1)*omoi(:max_pip) -                  &
               rot_acc_old(:max_pip,1))*dtsolid
        WHERE(particle_state(:max_pip) == normal_particle) &
            rot_acc_old(:max_pip,1) = tow(:max_pip,1)*omoi(:max_pip)
         tow(:max_pip,1) = zero

!!$omp section
         WHERE(particle_state == normal_particle) &
            omega_new(:max_pip,2) = omega_old(:max_pip,2) + 0.5d0*     &
                 (3.d0*tow(:max_pip,2)*omoi(:max_pip)-                 &
                 rot_acc_old(:max_pip,2) )*dtsolid
         WHERE(particle_state == normal_particle) &
            rot_acc_old(:max_pip,2) = tow(:max_pip,2)*omoi(:max_pip)
         tow(:max_pip,2) = zero

!!$omp section
        WHERE(particle_state == normal_particle) &
            omega_new(:max_pip,3) = omega_old(:max_pip,3) + 0.5d0*     &
               (3.d0*tow(:max_pip,3)*omoi(:max_pip)-                   &
               rot_acc_old(:max_pip,3) )*dtsolid
        WHERE(particle_state == normal_particle) &
            rot_acc_old(:max_pip,3) = tow(:max_pip,3)*omoi(:max_pip)
         tow(:max_pip,3) = zero

!!$omp end sections
      ENDIF

#ifdef do_sap
!!$omp single
         DO l = 1, max_pip
            aabb%minendpoint(:) = des_pos_new(l,:)-des_radius(l)
            aabb%maxendpoint(:) = des_pos_new(l,:)+des_radius(l)
            call multisap_update(multisap,aabb,boxhandle(l))
         ENDDO
!!$omp end single
#endif

! Update particle orientation - Always first order
! When omega is non-zero, compute the rotation angle, and apply the
! Rodrigues' rotation formula

      IF(particle_orientation) THEN
         DO l = 1, max_pip
            omega_mag = omega_new(l,1)**2 +omega_new(l,2)**2 + omega_new(l,3)**2

            IF(omega_mag>zero) THEN
               omega_mag=dsqrt(omega_mag)
               omega_unit(:) = omega_new(l,:)/omega_mag
               rot_angle = omega_mag * dtsolid

               orientation(:,l) = orientation(:,l)*dcos(rot_angle) + &
                  des_crossprdct(omega_unit,orientation(:,l))*dsin(rot_angle) + &
                  omega_unit(:)*dot_product(omega_unit,orientation(:,l))*&
                  (one-dcos(rot_angle))
            ENDIF
         ENDDO
      ENDIF

! Check if the particle has moved a distance greater than or equal to
! its radius since the last time a neighbor search was called. if so,
! make sure that neighbor is called in des_time_march
      IF(.NOT.do_nsearch) THEN
!!$omp do reduction (.or.:do_nsearch)
         DO l = 1, max_pip
            do_nsearch = do_nsearch .or. &
               (des_pos_new(l,1) - ppos(l,1))**2+              &
               (des_pos_new(l,2) - ppos(l,2))**2+              &
               (des_pos_new(l,3) - ppos(l,3))**2  >=           &
               (neighbor_search_rad_ratio*des_radius(l))**2
         ENDDO
      ENDIF

!!$omp end parallel

#ifdef do_sap
            call multisap_sort(multisap)
#endif

      first_pass = .false.

 1002 FORMAT(/1x,70('*')//&
         ' From: CFNEWVALUES -',/&
         ' Message: Particle ',i10, ' moved a distance ', es17.9, &
         ' during a',/10x, 'single solids time step, which is ',&
         ' greater than',/10x,'its radius: ', es17.9)
 1003 FORMAT(1x,70('*')/)

      RETURN

      contains

      include 'functions.inc'

      END SUBROUTINE cfnewvalues