des_allocate_mod.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: DES_ALLOCATE                                           C
!                                                                      C
!  Purpose: subroutines to allocate all DEM arrays                     C
!                                                                      C
!  Author: Rahul Garg                               Date: 1-Dec-2013   C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

MODULE des_allocate

  PUBLIC:: des_allocate_arrays, add_pair, particle_grow, allocate_dem_mi

CONTAINS

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: DES_ALLOCATE_ARRAYS                                     C
!  Purpose: Original allocate arrays subroutines for DES               C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE des_allocate_arrays

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE compar
      USE constant
      USE cutcell
      USE derived_types, only: boxhandle, pic
      USE des_bc
      USE des_rxns
      USE des_thermo
      USE discretelement
      USE functions

      USE des_thermo_cond, only: des_qw_cond 

      USE funits
      USE geometry
      USE indices
      USE mfix_pic
      USE param
      USE param1
      USE physprop
      USE pic_bc, only: pic_bcmo, pic_bcmi

      USE run, only: energy_eq
      USE run, only: any_species_eq

      USE particle_filter, only: des_interp_scheme_enum
      USE particle_filter, only: des_interp_garg
      USE particle_filter, only: des_interp_dpvm
      USE particle_filter, only: des_interp_gauss
      USE particle_filter, only: des_interp_lhat
      USE particle_filter, only: filter_size
      USE particle_filter, only: filter_cell
      USE particle_filter, only: filter_weight

! Use the error manager for posting error messages.
!---------------------------------------------------------------------//
      USE error_manager

      IMPLICIT NONE
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! indices
      INTEGER :: IJK
!-----------------------------------------------

      CALL init_err_msg("DES_ALLOCATE_ARRAYS")

! For parallel processing the array size required should be either
! specified by the user or could be determined from total particles
! with some factor.
      max_pip = merge(0, particles/numpes, particles==undefined_i)
      max_pip = max(max_pip,4)

      WRITE(err_msg,1000) trim(ival(max_pip))
      CALL flush_err_msg(header = .false., footer = .false.)

 1000 FORMAT('Initial DES Particle array size: ',a)

! DES Allocatable arrays
!-----------------------------------------------
! Dynamic particle info including another index for parallel
! processing for ghost
      ALLOCATE( particle_state(max_pip) )
      ALLOCATE (iglobal_id(max_pip))

! R.Garg: Allocate necessary arrays for PIC mass inlet/outlet BCs
      IF(pic_bcmi /= 0 .OR. pic_bcmo /=0) CALL allocate_pic_mio

! Particle attributes
! Radius, density, mass, moment of inertia
      Allocate(  des_radius(max_pip) )
      Allocate(  ro_sol(max_pip) )
      Allocate(  pvol(max_pip) )
      Allocate(  pmass(max_pip) )
      Allocate(  omoi(max_pip) )

! Old and new particle positions, velocities (translational and
! rotational)
      Allocate(  des_pos_new(max_pip,dimn) )
      Allocate(  des_vel_new(max_pip,dimn) )
      Allocate(  omega_new(max_pip,dimn) )

      IF(particle_orientation) Allocate(  orientation(dimn,max_pip) )

      IF (do_old) THEN
         Allocate(  des_pos_old(max_pip,dimn) )
         Allocate(  des_vel_old(max_pip,dimn) )
         Allocate(  des_acc_old(max_pip,dimn) )
         Allocate(  omega_old(max_pip,dimn) )
         Allocate(  rot_acc_old(max_pip,dimn))
      ENDIF

! Allocating user defined array
      IF(des_usr_var_size > 0) &
         Allocate( des_usr_var(des_usr_var_size,max_pip) )

! Particle positions at the last call neighbor search algorithm call
      Allocate(  ppos(max_pip,dimn) )

! Total, normal and tangetial forces
      Allocate(  fc(max_pip,dimn) )

! Torque
      Allocate(  tow(max_pip,dimn) )


! allocate variable for des grid binning
      allocate(dg_pijk(max_pip)); dg_pijk=0
      allocate(dg_pijkprv(max_pip)); dg_pijkprv=0

! allocate variables related to ghost particles
      allocate(ighost_updated(max_pip))



      Allocate(  wall_collision_facet_id(collision_array_max, max_pip) )
      wall_collision_facet_id(:,:) = -1
      Allocate(  wall_collision_pft(dimn, collision_array_max, max_pip) )

! Temporary variables to store wall position, velocity and normal vector
      Allocate(  wall_normal(nwalls,dimn) )

      Allocate(  neighbor_index(max_pip) )
      Allocate(  neighbor_index_old(max_pip) )
      Allocate(  neighbors(max_pip) )
      neighbors(:) = 0

      Allocate(  neighbors_old(max_pip) )
      Allocate(  pft_neighbor(3,max_pip) )
      Allocate(  pft_neighbor_old(3,max_pip) )
#ifdef do_sap
      Allocate(  boxhandle(max_pip) )
#endif

! Variable that stores the particle in cell information (ID) on the
! computational fluid grid defined by imax, jmax and kmax in mfix.dat
      ALLOCATE(pic(dimension_3))
      DO ijk=1,dimension_3
        NULLIFY(pic(ijk)%p)
      ENDDO

! Particles in a computational fluid cell (for volume fraction)
      Allocate(  pinc(dimension_3) )

! For each particle track its i,j,k location on computational fluid grid
! defined by imax, jmax and kmax in mfix.dat and phase no.
      Allocate(  pijk(max_pip,5) )

      ALLOCATE(drag_am(dimension_3))
      ALLOCATE(drag_bm(dimension_3, dimn))
      ALLOCATE(f_gp(max_pip ))
      f_gp(1:max_pip)  = zero

! Explicit drag force acting on a particle.
      Allocate(drag_fc(max_pip,dimn) )

! force due to gas-pressure gradient
      ALLOCATE(p_force(dimn, dimension_3))

! Volume of nodes
      ALLOCATE(des_vol_node(dimension_3))

      ALLOCATE(f_gds(dimension_3))
      ALLOCATE(vxf_gds(dimension_3))

      SELECT CASE(des_interp_scheme_enum)
      CASE(des_interp_dpvm, des_interp_gauss, des_interp_lhat)
         ALLOCATE(filter_cell(filter_size, max_pip))
         ALLOCATE(filter_weight(filter_size, max_pip))
      CASE(des_interp_garg)
         ALLOCATE(des_rops_node(dimension_3, dimension_m))
         ALLOCATE(des_vel_node(dimension_3, dimn, dimension_m))
      END SELECT

! Variables for hybrid model
      IF (des_continuum_hybrid) THEN
         ALLOCATE(sdrag_am(dimension_3,dimension_m))
         ALLOCATE(sdrag_bm(dimension_3, dimn,dimension_m))

         ALLOCATE(f_sds(dimension_3,dimension_m))
         ALLOCATE(vxf_sds(dimension_3,dimension_m))
      ENDIF

! MP-PIC related
      IF(mppic) THEN
         Allocate(ps_force_pic(3, dimension_3))
         ALLOCATE(des_stat_wt(max_pip))
         ALLOCATE(des_vel_max(dimn))
         ALLOCATE(ps_grad(3,max_pip))
         ALLOCATE(avgsolvel_p(3, max_pip))
         ALLOCATE(epg_p(max_pip))

         Allocate(pic_u_s(dimension_3, dimension_m) )
         Allocate(pic_v_s(dimension_3, dimension_m) )
         Allocate(pic_w_s(dimension_3, dimension_m) )
         Allocate(pic_p_s(dimension_3, dimension_m) )
         pic_u_s = zero
         pic_v_s = zero
         pic_w_s = zero
         pic_p_s = zero
      ENDIF

! Averaged velocity obtained by averaging over all the particles
      ALLOCATE(des_vel_avg(dimn) )

! Global Granular Energy
      ALLOCATE(global_gran_energy(dimn) )
      ALLOCATE(global_gran_temp(dimn) )

! variable for bed height of solids phase M
      ALLOCATE(bed_height(dimension_m))

! ---------------------------------------------------------------->>>
! BEGIN COHESION
      IF(use_cohesion) THEN
! Matrix location of particle  (should be allocated in case user wishes
! to invoke routines in /cohesion subdirectory
         Allocate(  postcohesive(max_pip) )
      ENDIF
! END COHESION
! ----------------------------------------------------------------<<<

! ---------------------------------------------------------------->>>
! BEGIN Thermodynamic Allocation
      IF(energy_eq)THEN
! Particle temperature
         Allocate( des_t_s( max_pip ) )
! Specific heat
         Allocate( des_c_ps( max_pip ) )
! Species mass fractions comprising a particle. This array may not be
! needed for all thermo problems.
         Allocate( des_x_s( max_pip, dimension_n_s))
! Total rate of heat transfer to individual particles.
         Allocate( q_source( max_pip ) )
! Average solids temperature in fluid cell
         Allocate(avgdes_t_s(dimension_3) )
! Gas/Solids convective heat transfer coupling
         IF(calc_conv_des) THEN
! Fluid phase energy equation source terms
            Allocate(conv_sc(dimension_3) )
            Allocate(conv_sp(dimension_3) )
! Particle convection source term (explicit coupled)
            Allocate(conv_qs(max_pip))
! Gas-particle heat transfer coefficient TIMES surface area
            Allocate(gammaxsa(max_pip))
         ENDIF

! Allocate the history variables for Adams-Bashforth integration
         IF (intg_adams_bashforth) &
            Allocate( q_source0( max_pip ) )
! Allocate the array for storing particle-wall heat transfer per unit area
         IF (any(calc_cond_des)) &
            Allocate( des_qw_cond( dimension_3, dimension_m))
      ENDIF
! End Thermodynamic Allocation
! ----------------------------------------------------------------<<<


! ---------------------------------------------------------------->>>
! BEGIN Species Allocation
      IF(any_species_eq)THEN
! Rate of solids phase production/consumption for each species
         Allocate( des_r_s( max_pip, dimension_n_s) )


         Allocate( des_r_gp( dimension_3, dimension_n_g ) )
         Allocate( des_r_gc( dimension_3, dimension_n_g ) )
         Allocate( des_sum_r_g( dimension_3 ) )
         Allocate( des_r_phase( dimension_3, dimension_lm+dimension_m-1 ) )
         Allocate( des_hor_g( dimension_3 ) )


! Allocate the history variables for Adams-Bashforth integration
         IF (intg_adams_bashforth) THEN
! Rate of change of particle mass
            Allocate( dmdt_old( max_pip ) )
! Rate of change of particle mass percent species
            Allocate( dxdt_old( max_pip, dimension_n_s) )
         ENDIF

! Energy generation from reaction (cal/sec)
         Allocate( rxns_qs( max_pip ) )
      ENDIF
! End Species Allocation
! ----------------------------------------------------------------<<<

      CALL finl_err_msg

      RETURN
      END SUBROUTINE des_allocate_arrays

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: ALLOCATE_DEM_MIO                                        !
!                                                                      !
!  Purpose:                                                            !
!                                                                      !
!  Author: J.Musser                                   Date: 17-Aug-09  !
!                                                                      !
!  Comments:                                                           !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!

      SUBROUTINE allocate_dem_mi

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param1, only: undefined
      USE des_bc, only: dem_bcmi
      USE des_bc, only: pi_factor, pi_count
      use des_bc, only: numfrac_limit
      use des_bc, only: dem_mi_time, dem_bc_poly_layout
      use des_bc, only: dem_mi
      use des_bc, only: dem_bcmi_ijkstart, dem_bcmi_ijkend
      IMPLICIT NONE
!-----------------------------------------------

! Particle injection factor
      Allocate( pi_factor(dem_bcmi) )
! Particle injection count (injection number)
      Allocate( pi_count(dem_bcmi) )
! Particle injection time scale
      Allocate( dem_mi_time(dem_bcmi) )
! Array used for polydisperse inlets: stores the particle number
! distribution of an inlet scaled with numfrac_limit
      Allocate( dem_bc_poly_layout( dem_bcmi, numfrac_limit ) )
! Data structure for storing BC data.
      Allocate( dem_mi(dem_bcmi) )

! Initializiation
! Integer arrays
      pi_factor(:) = -1
      pi_count(:) = -1
      dem_bc_poly_layout(:,:) = -1
! Double precision arrays
      dem_mi_time(:) = undefined

      allocate( dem_bcmi_ijkstart(dem_bcmi) )
      allocate( dem_bcmi_ijkend(dem_bcmi) )

      dem_bcmi_ijkstart = -1
      dem_bcmi_ijkend   = -1

! Boundary classification
!         Allocate( PARTICLE_PLCMNT (DES_BCMI) )
! Character precision arrays
!         PARTICLE_PLCMNT(:) = UNDEFINED_C

      RETURN
      END SUBROUTINE allocate_dem_mi


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: ALLOCATE_PIC_MIO                                        !
!                                                                      !
!  Purpose:                                                            !
!                                                                      !
!  Author: R. Garg                                    Date: 11-Jun-14  !
!                                                                      !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE allocate_pic_mio

! Modules
!-----------------------------------------------
      USE pic_bc, only: pic_bcmi, pic_bcmo
      USE pic_bc, only: pic_bcmi_ijkstart, pic_bcmi_ijkend
      USE pic_bc, only: pic_bcmo_ijkstart, pic_bcmo_ijkend
      USE pic_bc, only: pic_bcmi_normdir
      USE pic_bc, only: pic_bcmi_offset
      USE pic_bc, only: pic_bcmi_incl_cutcell
      IMPLICIT NONE
!-----------------------------------------------

! Allocate/Initialize for inlets
      IF(pic_bcmi /= 0)THEN
         allocate( pic_bcmi_ijkstart(pic_bcmi) )
         allocate( pic_bcmi_ijkend(pic_bcmi) )
         allocate( pic_bcmi_normdir(pic_bcmi,3) )
         ALLOCATE( pic_bcmi_offset(pic_bcmi,3))
         ALLOCATE( pic_bcmi_incl_cutcell(pic_bcmi) )
         pic_bcmi_ijkstart = -1
         pic_bcmi_ijkend   = -1
      ENDIF  ! end if PIC_BCMI /= 0


      IF(pic_bcmo > 0)THEN
         allocate( pic_bcmo_ijkstart(pic_bcmo) )
         allocate( pic_bcmo_ijkend(pic_bcmo) )
         pic_bcmo_ijkstart = -1
         pic_bcmo_ijkend   = -1
      ENDIF

      RETURN
      END SUBROUTINE allocate_pic_mio



!``````````````````````````````````````````````````````````````````````!
! Subroutine: ADD_PAIR                                                 !
!                                                                      !
! Purpose: Adds a neighbor pair to the pairs array.                    !
!                                                                      !
!``````````````````````````````````````````````````````````````````````!
      DOUBLE PRECISION FUNCTION add_pair(ii,jj)
      USE discretelement
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: ii,jj

      CALL neighbor_grow(neighbor_index(ii))

      neighbors(neighbor_index(ii)) = jj
      neighbor_index(ii) = neighbor_index(ii) + 1
      add_pair = neighbor_index(ii)

      RETURN
      END FUNCTION add_pair

!``````````````````````````````````````````````````````````````````````!
! Subroutine: NEIGHBOR_GROW                                            !
!                                                                      !
! Purpose: Grow neighbors arrays to new_neigh_max. Note that neighbor      !
! max should be increased before calling this routine. Also, no        !
! assumption to the previous array size is made as needed for restarts.!
!``````````````````````````````````````````````````````````````````````!
      SUBROUTINE neighbor_grow(new_neigh_max)
        USE discretelement
        USE geometry
        IMPLICIT NONE

        integer, intent(in) :: new_neigh_max

        INTEGER :: lSIZE1
        INTEGER, DIMENSION(:), ALLOCATABLE :: neigh_tmp
        DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: pf_tmp
        INTEGER new_size

        lsize1 = size(neighbors,1)

        IF ( new_neigh_max .le. lsize1 ) RETURN

        new_size = lsize1

        DO WHILE(new_size < new_neigh_max)
           new_size = 2*new_size
        ENDDO

        allocate(neigh_tmp(new_size))
        neigh_tmp(1:lsize1) = neighbors(1:lsize1)
        neigh_tmp(lsize1+1:) = 0
        call move_alloc(neigh_tmp,neighbors)

        allocate(neigh_tmp(new_size))
        neigh_tmp(1:lsize1) = neighbors_old(1:lsize1)
        neigh_tmp(lsize1+1:) = 0
        call move_alloc(neigh_tmp,neighbors_old)

        allocate(pf_tmp(3,new_size))
        pf_tmp(:,1:lsize1) = pft_neighbor(:,1:lsize1)
        pf_tmp(:,lsize1+1:) = 0
        call move_alloc(pf_tmp,pft_neighbor)

        allocate(pf_tmp(3,new_size))
        pf_tmp(:,1:lsize1) = pft_neighbor_old(:,1:lsize1)
        pf_tmp(:,lsize1+1:) = 0
        call move_alloc(pf_tmp,pft_neighbor_old)


      END SUBROUTINE neighbor_grow

!``````````````````````````````````````````````````````````````````````!
! Subroutine: PARTICLE_GROW                                            !
!                                                                      !
! Purpose: Grow particle arrays to new_max_pip. Note that pair         !
! max should be increased before calling this routine. Also, no        !
! assumption to the previous array size is made as needed for restarts.!
!``````````````````````````````````````````````````````````````````````!
      SUBROUTINE particle_grow(new_max_pip)

        USE des_rxns
        USE des_thermo
        USE derived_types, only: boxhandle
        USE discretelement
        USE mfix_pic
        USE multi_sweep_and_prune, ONLY: boxhandle_grow
        USE particle_filter
        USE resize
        USE run

        IMPLICIT NONE

        integer, intent(in) :: new_max_pip
        integer :: old_size, new_size

        IF (new_max_pip .le. size(des_radius)) RETURN
        max_pip = max(max_pip, new_max_pip)

        old_size = size(des_radius)

        new_size = old_size

        DO WHILE (new_size < new_max_pip)
           new_size = 2*new_size
        ENDDO

#ifdef do_sap
        call boxhandle_grow(boxhandle,new_size)
#endif
        call real_grow(des_radius,new_size)
        call real_grow(ro_sol,new_size)
        call real_grow(pvol,new_size)
        call real_grow(pmass,new_size)
        call real_grow(omoi,new_size)
        call real_grow2_reverse(des_pos_new,new_size)
        call real_grow2_reverse(des_vel_new,new_size)
        call real_grow2_reverse(omega_new,new_size)
        call real_grow2_reverse(ppos,new_size)
        call byte_grow(particle_state,new_size)
        call integer_grow(iglobal_id,new_size)
        call integer_grow2_reverse(pijk,new_size)
        call integer_grow(dg_pijk,new_size)
        call integer_grow(dg_pijkprv,new_size)
        call logical_grow(ighost_updated,new_size)
        call real_grow2_reverse(fc,new_size)
        call real_grow2_reverse(tow,new_size)
        call real_grow(f_gp,new_size)
        call integer_grow2(wall_collision_facet_id,new_size)
        call real_grow3(wall_collision_pft,new_size)
        call real_grow2_reverse(drag_fc,new_size)

        call integer_grow(neighbor_index,new_size)
        call integer_grow(neighbor_index_old,new_size)

        IF(particle_orientation) call real_grow2(orientation,new_size)

        IF(filter_size > 0) THEN
           call integer_grow2(filter_cell,new_size)
           call real_grow2(filter_weight,new_size)
        ENDIF

        IF(mppic) THEN
           call real_grow(des_stat_wt,new_size)
           call real_grow2(ps_grad,new_size)
           call real_grow2(avgsolvel_p,new_size)
           call real_grow(epg_p,new_size)
        ENDIF

        IF(use_cohesion) THEN
           call real_grow(postcohesive,new_size)
        ENDIF

        IF (do_old) THEN
           call real_grow2_reverse(des_pos_old,new_size)
           call real_grow2_reverse(des_vel_old,new_size)
           call real_grow2_reverse(des_acc_old,new_size)
           call real_grow2_reverse(omega_old,new_size)
           call real_grow2_reverse(rot_acc_old,new_size)
        ENDIF

        IF(energy_eq)THEN
           call real_grow(des_t_s,new_size)
           call real_grow(des_c_ps,new_size)
           call real_grow2_reverse(des_x_s,new_size)
           call real_grow(q_source,new_size)
           IF(calc_conv_des) THEN
              call real_grow(conv_qs, new_size)
              call real_grow(gammaxsa, new_size)
           ENDIF
           IF(intg_adams_bashforth) &
              call real_grow(q_source0,new_size)
        ENDIF

        IF(any_species_eq)THEN
           call real_grow2_reverse( des_r_s, new_size )

           IF (intg_adams_bashforth) THEN
              call real_grow( dmdt_old, new_size )
              call real_grow2_reverse( dxdt_old, new_size )
           ENDIF

           call real_grow( rxns_qs, new_size )
        ENDIF

        IF(des_usr_var_size > 0) &
        call real_grow2(des_usr_var,new_size)

        CALL des_init_particle_arrays(old_size+1,new_size)

        RETURN

      END SUBROUTINE particle_grow

    END MODULE des_allocate