comp_mean_fields0.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE comp_mean_fields0

! Modules
      USE bc
      USE compar
      USE constant
      USE cutcell
      USE derived_types, only: pic
      USE desmpi
      USE discretelement
      USE drag
      USE fldvar
      USE functions, only: fluid_at
      USE functions, only: funijk
      USE functions, only: is_on_mype_wobnd
      USE geometry
      USE indices
      USE interpolation
      USE mfix_pic
      USE mpi_node_des, only: des_addnodevalues_mean_fields
      USE mpi_utility
      USE parallel
      USE param
      USE param1
      USE particle_filter, only: des_report_mass_interp
      USE physprop, only: mmax
      USE run, only: solids_model
      USE sendrecv

      IMPLICIT NONE
! Local variables
!---------------------------------------------------------------------//
! general i, j, k indices
      INTEGER :: I, J, K, IJK, &
                 II, JJ, KK
      INTEGER :: I1, I2, J1, J2, K1, K2
      INTEGER :: IDIM, IJK2
      INTEGER :: CUR_IJK
! indices used for interpolation stencil (unclear why IE, JN, KTP are
! needed)
      INTEGER :: IW, IE, JS, JN, KB, KTP
! i,j,k indices of the fluid cell the particle resides in minus 1
! (e.g., shifted 1 in west, south, bottom direction)
      INTEGER, DIMENSION(3):: PCELL
! order of interpolation set in the call to set_interpolation_scheme
! unless it is re/set later through the call to
! set_interpolation_stencil
      INTEGER :: ONEW
! constant whose value depends on dimension of system
! avg_factor=0.250 (in 3D) or =0.50 (in 2D)
      DOUBLE PRECISION :: AVG_FACTOR
! index of solid phase that particle NP belongs to
      INTEGER :: M
! particle number index, used for looping
      INTEGER :: NP, NINDX
! Statistical weight of the particle. Equal to one for DEM
      DOUBLE PRECISION :: WTP

      DOUBLE PRECISION :: MASS_SOL1, MASS_SOL2
! sum of mass_sol1 and mass_sol2 across all processors
      DOUBLE PRECISION :: MASS_SOL1_ALL, MASS_SOL2_ALL

      DOUBLE PRECISION :: TEMP1, TEMP2

      DOUBLE PRECISION, DIMENSION(3) :: DES_VEL_DENSITY
      DOUBLE PRECISION :: DES_ROP_DENSITY

      INTEGER :: COUNT_NODES_OUTSIDE, COUNT_NODES_INSIDE, &
                 COUNT_NODES_INSIDE_MAX

      DOUBLE PRECISION :: RESID_ROPS(dimension_m)
      DOUBLE PRECISION :: RESID_VEL(3, dimension_m)
      DOUBLE PRECISION :: NORM_FACTOR

!Handan Liu added on Jan 17 2013
      DOUBLE PRECISION, DIMENSION(2,2,2,3) :: gst_tmp
      DOUBLE PRECISION, DIMENSION(2,2,2) :: weight_ft

! total number of 'solids' phases in simulation
      INTEGER :: MMAX_TOT
!......................................................................!

! initializing
      mass_sol1 = zero
      mass_sol2 = zero
      mass_sol1_all = zero
      mass_sol2_all = zero
! avg_factor=0.25 (in 3D) or =0.5 (in 2D)
      avg_factor = merge(0.5d0, 0.25d0, no_k)

! cartesian_grid related quantities
      count_nodes_inside_max = merge(4, 8, no_k)

      mmax_tot = des_mmax+mmax
! initialize only information related to the discrete 'phases' of
! these continuous variables
      rop_s(:,mmax+1:mmax_tot) = zero
      u_s(:,mmax+1:mmax_tot) = zero
      v_s(:,mmax+1:mmax_tot) = zero
      IF(do_k) w_s(:,mmax+1:mmax_tot) = zero
! these are not 'continuous' variables but used only within this routine...
      des_vel_node = zero
      des_rops_node = zero

! sets several quantities including interp_scheme, scheme, and
! order and allocates arrays necessary for interpolation
      CALL set_interpolation_scheme(2)

!$omp parallel default(shared) &
!$omp private(IJK, I, J, K, PCELL, IW, IE, JS, JN, KB, KTP, ONEW, &
!$omp         GST_TMP, COUNT_NODES_INSIDE, II, JJ, KK, CUR_IJK, &
!$omp         NINDX, NP, WTP, M, WEIGHT_FT, I1, I2, J1, J2, K1, K2, &
!$omp         IDIM, IJK2, NORM_FACTOR, RESID_ROPS, RESID_VEL, &
!$omp         COUNT_NODES_OUTSIDE, TEMP1)
!$omp do reduction(+:MASS_SOL1) reduction(+:DES_ROPS_NODE,DES_VEL_NODE)
      DO ijk = ijkstart3,ijkend3

! Cycle this cell if not in the fluid domain or if it contains no
! particle/parcel
         IF(.NOT.fluid_at(ijk)) cycle
         IF( pinc(ijk) == 0) cycle

         pcell(1) = i_of(ijk)-1
         pcell(2) = j_of(ijk)-1
         pcell(3) = merge(k_of(ijk)-1, 1, do_k)

! setup the stencil based on the order of interpolation and factoring in
! whether the system has any periodic boundaries. sets onew to order.
         CALL set_interpolation_stencil(pcell, iw, ie, js, jn, kb, ktp,&
            interp_scheme, dimn, ordernew=onew)

         count_nodes_outside = 0
! Computing/setting the geometric stencil
         DO k=1, merge(1, onew, no_k)
         DO j=1, onew
         DO i=1, onew
            ii = iw + i-1
            jj = js + j-1
            kk = kb + k-1
            cur_ijk = funijk_map_c(ii,jj,kk)
            gst_tmp(i,j,k,1) = xe(ii)
            gst_tmp(i,j,k,2) = yn(jj)
            gst_tmp(i,j,k,3) = merge(dz(1), zt(kk), no_k)
            IF(cartesian_grid) THEN
               IF(scalar_node_atwall(cur_ijk))                   &
                  count_nodes_outside = count_nodes_outside + 1
            ENDIF
         ENDDO
         ENDDO
         ENDDO


! Calculate des_rops_node so rop_s, and in turn, ep_g can be updated
!----------------------------------------------------------------->>>

! looping through particles in the cell
         DO nindx=1, pinc(ijk)
            np = pic(ijk)%P(nindx)
            call drag_weightfactor(gst_tmp,des_pos_new(np,:),weight_ft)
            m = pijk(np,5)
            wtp = one

            IF(mppic) wtp = des_stat_wt(np)

            mass_sol1 = mass_sol1 + pmass(np)*wtp
            temp2 = pmass(np)*wtp
            DO k = 1, merge(1, onew, no_k)
            DO j = 1, onew
            DO i = 1, onew
! shift loop index to new variables for manipulation
               ii = iw + i-1
               jj = js + j-1
               kk = kb + k-1
               cur_ijk = funijk_map_c(ii,jj,kk)
               temp1 = weight_ft(i,j,k)*temp2

               des_rops_node(cur_ijk,m) = des_rops_node(cur_ijk,m) +   &
                  temp1
               des_vel_node(cur_ijk,:,m) = des_vel_node(cur_ijk,:,m) + &
                  temp1*des_vel_new(np,:)
            ENDDO
            ENDDO
            ENDDO
         ENDDO   ! end do (nindx=1,pinc(ijk))
!-----------------------------------------------------------------<<<


! Only for cutcell cases may count_nodes_inside become less than its
! original set value. In such an event, the contribution of scalar nodes
! that do not reside in the domain is added to a residual array. This
! array is then redistribited equally to the nodes that are in the fluid
! domain. These steps are done to conserve mass.
!----------------------------------------------------------------->>>
! initializing
         resid_rops = zero
         resid_vel = zero
         IF (cartesian_grid) THEN

! only for cartesian_grid will count_nodes_outside be modified from zero
            count_nodes_inside = count_nodes_inside_max - &
                                 count_nodes_outside

            IF(count_nodes_inside.LT.count_nodes_inside_max) THEN

! Convention used to number node numbers
! i=1, j=2           i=2, j=2
!   _____________________
!   |                   |
!   |  I = 2, J = 2     |
!   |___________________|
! i=1, j=1           i=2, j=1
! setting indices based on convention
               i = i_of(ijk)
               j = j_of(ijk)
               k = k_of(ijk)
               i1 = i-1
               i2 = i
               j1 = j-1
               j2 = j
               k1 = merge(k, k-1, no_k)
               k2 = k

! first calculate the residual des_rops_node and des_vel_node that was
! computed on nodes that do not belong to the domain
               DO kk = k1, k2
               DO jj = j1, j2
               DO ii = i1, i2
                  ijk2 = funijk(ii, jj, kk)
                  IF(scalar_node_atwall(ijk2)) THEN
                     DO m = mmax+1,mmax_tot
                        resid_rops(m) = resid_rops(m) + &
                           des_rops_node(ijk2,m)
                        des_rops_node(ijk2,m) = zero

                        DO idim = 1, merge(2,3,no_k)
                           resid_vel(idim,m) = resid_vel(idim, m) + &
                              des_vel_node(ijk2,idim, m)
                           des_vel_node(ijk2,idim, m) = zero
                        ENDDO
                     ENDDO
                  ENDIF
               ENDDO
               ENDDO
               ENDDO

! now add this residual equally to the remaining nodes
               norm_factor = one/REAL(count_nodes_inside)
               DO kk = k1, k2
               DO jj = j1, j2
               DO ii = i1, i2
                  ijk2 = funijk(ii, jj, kk)
                  IF(.NOT.scalar_node_atwall(ijk2)) THEN
                     DO m = mmax+1,mmax_tot
                        des_rops_node(ijk2,m) = &
                           des_rops_node(ijk2,m) + &
                           resid_rops(m)*norm_factor
                        DO idim = 1, merge(2,3,no_k)
                           des_vel_node(ijk2,idim, m) = &
                              des_vel_node(ijk2,idim, m) + &
                              resid_vel(idim, m)*norm_factor
                        ENDDO
                     ENDDO
                  ENDIF
               ENDDO
               ENDDO
               ENDDO

            ENDIF
         ENDIF   ! end if (cartesian_grid)
      ENDDO
!$omp end parallel


! At the interface des_rops_node has to be added since particles
! across the processors will contribute to the same scalar node.
! sendrecv will be called and the node values will be added
! at the junction. des_rops_node is altered by the routine when
! periodic boundaries are invoked
      CALL des_addnodevalues_mean_fields


! Now go from node to scalar center. Same convention as sketched
! earlier
!----------------------------------------------------------------->>>
! Explanation by RG: 08/17/2012
! the approach used here is to make it general enough for cutcells to be
! included as well. The new changes do not alter earlier calculations
! but make the technique general as to include cartesian grid (cut-cell)
! simulations.
! Previously, the volume of the node (by array des_vol_node) was used to
! first scale the nodal values. Subsequently, these nodal values were
! equally added to compute the cell centered values for the scalar cell.

! Consider an internal node next to an edge node (a node adjacent to a
! boundary). In 2D, the volume of an edge node will be half that of an
! internal node. And an edge node will contribute double compared to
! an internal node to the value of the scalar cell they share. These
! calculations were previously accomplished via the variable volume of
! node.  Now this is accomplished by the ratio vol(ijk2)/vol_sur, where
! vol(ijk2) is the volume of the scalar cell in consideration and
! vol_sur is the sum of all the scalar cell volumes that have this node
! as the common node.

!$omp parallel do default(none) collapse (3) &
!$omp shared(KSTART2, KEND1, JSTART2, JEND1, ISTART2, IEND1, DO_K, &
!$omp        VOL, DEAD_CELL_AT, FUNIJK_MAP_C, VOL_SURR, MMAX_TOT, &
!$omp        MMAX, DES_ROPS_NODE, DES_VEL_NODE) &
!$omp private(I, J, K, IJK, M, II, JJ, KK, IJK2, DES_ROP_DENSITY, &
!$omp         DES_VEL_DENSITY) &
!$omp reduction(+:ROP_S, U_S, V_S, W_S)
      DO k = kstart2, kend1
      DO j = jstart2, jend1
      DO i = istart2, iend1
         IF (dead_cell_at(i,j,k)) cycle  ! skip dead cells
         ijk = funijk(i,j,k)
         if (vol_surr(ijk).eq.zero) cycle ! no FLUID_AT any of the stencil points have

! looping over stencil points (NODE VALUES)
         DO m = mmax+1, mmax_tot
            des_rop_density = des_rops_node(ijk, m)/vol_surr(ijk)
            des_vel_density(:) = des_vel_node(ijk, :, m)/vol_surr(ijk)

            DO kk = k, merge(k+1, k, do_k)
            DO jj = j, j+1
            DO ii = i, i+1
               IF (dead_cell_at(ii,jj,kk)) cycle  ! skip dead cells

               ijk2 = funijk_map_c(ii, jj, kk)
               IF(fluid_at(ijk2).and.(is_on_mype_wobnd(ii, jj, kk))) THEN
! Since the data in the ghost cells is spurious anyway and overwritten during
! subsequent send receives, do not compute any value here as this will
! mess up the total mass value that is computed below to ensure mass conservation
! between Lagrangian and continuum representations
                  rop_s(ijk2, m) = rop_s(ijk2, m) + des_rop_density*vol(ijk2)
                  u_s(ijk2, m) = u_s(ijk2, m) + des_vel_density(1)*vol(ijk2)
                  v_s(ijk2, m) = v_s(ijk2, m) + des_vel_density(2)*vol(ijk2)
                  IF(do_k) w_s(ijk2, m) = w_s(ijk2, m) + des_vel_density(3)*vol(ijk2)
               ENDIF
            ENDDO  ! end do (ii=i1,i2)
            ENDDO  ! end do (jj=j1,j2)
            ENDDO  ! end do (kk=k1,k2)
         ENDDO
      ENDDO   ! end do (i=istart2,iend1)
      ENDDO   ! end do (j=jstart2,jend1)
      ENDDO   ! end do (k=kstart2,kend1)
!omp end parallel do
!-----------------------------------------------------------------<<<


!$omp parallel do default(none) private(IJK, M) &
!$omp shared(IJKSTART3, IJKEND3, DO_K, MMAX_TOT, MMAX, ROP_s, U_S, &
!$omp        V_S, W_S, VOL)
      DO ijk = ijkstart3, ijkend3
         IF(.NOT.fluid_at(ijk)) cycle
         DO m = mmax+1, mmax_tot
            IF(rop_s(ijk, m).GT.zero) THEN
               u_s(ijk, m) = u_s(ijk,m)/rop_s(ijk, m)
               v_s(ijk, m) = v_s(ijk,m)/rop_s(ijk, m)
               IF(do_k) w_s(ijk, m) = w_s(ijk,m)/rop_s(ijk, m)
! Divide by scalar cell volume to obtain the bulk density
               rop_s(ijk, m) = rop_s(ijk, m)/vol(ijk)
            ENDIF
         ENDDO   ! end loop over M=1,MMAX_TOT
      ENDDO  ! end loop over IJK=ijkstart3,ijkend3
!omp end parallel do


      IF (mppic) CALL send_recv(rop_s,2)
      CALL calc_epg_des

      IF(mppic) THEN
! Now calculate Eulerian mean velocity fields for discrete phases
         CALL send_recv(u_s,2)
         CALL send_recv(v_s,2)
         IF(do_k) CALL send_recv(w_s,2)

! The Eulerian velocity field is used to set up the stencil to interpolate
! mean solid velocity at the parcel's location. U_S could have also been
! used, but that also would have require the communication at this stage.
! The final interpolated value does not change if the stencil is formed by
! first obtaining face centered Eulerian velocities (U_S, etc.)
! and then computing the node velocities from them or directly computing
! the node velocities from cell centered average velocity field (U_S,
! etc.). We are using the first approach as it is more natural to set
! BC's on solid velocity field in the face centered represenation (U_S,
! etc.)
         IF(.NOT.cartesian_grid) THEN
            CALL mppic_comp_eulerian_vels_non_cg
         ELSE
            CALL mppic_comp_eulerian_vels_cg
         ENDIF
      ENDIF   ! end if (.not.mppic)

! turn on the below statements to check if the mass is conserved
! between discrete and continuum representations. Should be turned to
! false for any production runs.
      IF(des_report_mass_interp) THEN
         DO ijk = ijkstart3, ijkend3
            IF(.NOT.fluid_at(ijk)) cycle
            i = i_of(ijk)
            j = j_of(ijk)
            k = k_of(ijk)
! It is important to check both FLUID_AT and IS_ON_MYPE_WOBND.
            IF(is_on_mype_wobnd(i,j,k)) THEN
               DO m=mmax+1,mmax_tot
                  mass_sol2 = mass_sol2 + &
                     rop_s(ijk,m)*vol(ijk)
               ENDDO
            ENDIF
         ENDDO
         CALL global_sum(mass_sol1, mass_sol1_all)
         CALL global_sum(mass_sol2, mass_sol2_all)
         if(mype.eq.pe_io) THEN
            WRITE(*,'(/,5x,A,4(2x,g17.8),/)') &
              'SOLIDS MASS DISCRETE AND CONTINUUM =  ', &
              mass_sol1_all, mass_sol2_all
         ENDIF
      ENDIF

      END SUBROUTINE comp_mean_fields0



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!  Subroutine: DRAG_WEIGHTFACTOR                                        C
!  Purpose: DES - Calculate the fluid velocity interpolated at the      C
!           particle's location and weights. Replace 'interpolator'     C
!                       interface for OpenMP implementation.            C
!                                                                       C
!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC

      SUBROUTINE drag_weightfactor(GSTEN,DESPOS,WEIGHTFACTOR)

      use geometry, only: no_k

        IMPLICIT NONE

!-----------------------------------------------
! Local Variables
!-----------------------------------------------
        DOUBLE PRECISION, DIMENSION(2,2,2,3), INTENT(IN):: GSTEN
        DOUBLE PRECISION, DIMENSION(3), INTENT(IN):: DESPOS
        DOUBLE PRECISION, DIMENSION(2,2,2), INTENT(OUT) :: WEIGHTFACTOR
        INTEGER :: II, JJ, KK

        DOUBLE PRECISION, DIMENSION(2) :: XXVAL, YYVAL, ZZVAL
        DOUBLE PRECISION :: DXX, DYY, DZZ
        DOUBLE PRECISION, DIMENSION(3) :: ZETAA

        dxx = gsten(2,1,1,1) - gsten(1,1,1,1)
        dyy = gsten(1,2,1,2) - gsten(1,1,1,2)

        zetaa(1:2) = despos(1:2) - gsten(1,1,1,1:2)

        zetaa(1) = zetaa(1)/dxx
        zetaa(2) = zetaa(2)/dyy

        xxval(1)=1-zetaa(1)
        yyval(1)=1-zetaa(2)
        xxval(2)=zetaa(1)
        yyval(2)=zetaa(2)

        IF(no_k) THEN
           DO jj=1,2
              DO ii=1,2
                 weightfactor(ii,jj,1) = xxval(ii)*yyval(jj)
              ENDDO
           ENDDO
        ELSE
           dzz = gsten(1,1,2,3) - gsten(1,1,1,3)
           zetaa(3) = despos(3) - gsten(1,1,1,3)
           zetaa(3) = zetaa(3)/dzz
           zzval(1)=1-zetaa(3)
           zzval(2)=zetaa(3)
           DO kk=1,2
              DO jj=1,2
                 DO ii=1,2
                    weightfactor(ii,jj,kk) = xxval(ii)*yyval(jj)*zzval(kk)
                 ENDDO
              ENDDO
           ENDDO
        ENDIF

      END SUBROUTINE drag_weightfactor