solve_continuity.f

Go to the documentation of this file.

! -*- f90 -*-
MODULE cont

! Indicates whether the continuity equation needs to be
! solved

      LOGICAL, DIMENSION(:), ALLOCATABLE ::  do_cont

CONTAINS

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SOLVE_CONTINUITY                                        C
!  Purpose: Solve for solids bulk density (i.e., material density      C
!           multiplied by volume fraction).                            C
!                                                                      C
!  Author: M. Syamlal                                 Date: 2-JUL-96   C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE solve_continuity(M,IER)

! Modules
!---------------------------------------------------------------------//
      use ambm, only: a_m, b_m, lock_ambm, unlock_ambm
      use fldvar, only: rop_g, rop_s
      use geometry, only: ijkmax2
      use leqsol, only: leq_it, leq_method, leq_sweep, leq_pc, leq_tol
      use param, only: dimension_3, dimension_m
      use param1, only: undefined_i, zero, one
      use ps, only: point_source
      use residual, only: resid, max_resid, ijk_resid
      use residual, only: num_resid, den_resid
      use residual, only: resid_ro
      use usr_src, only: call_usr_source, calc_usr_source
      use usr_src, only: gas_continuity, solids_continuity
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! phase index
      INTEGER, INTENT(IN) :: M
! error index
      INTEGER, INTENT(INOUT) :: IER

! Local variables
!---------------------------------------------------------------------//
! solids volume fraction residual
      DOUBLE PRECISION :: RESs
! linear equation solver method and iterations
      INTEGER :: LEQM, LEQI

! temporary use of global arrays:
! Septadiagonal matrix A_m, vector B_m
!      DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
!      DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
!---------------------------------------------------------------------//

      call lock_ambm

      IF (m==0) THEN
! solve gas phase continuity equation. note that this branch will never
! be entered given that the calling subroutine (iterate) only calls
! solve_continuity when M>1

! initializing
         CALL init_ab_m (a_m, b_m, ijkmax2, 0)

! forming the matrix equation
         CALL conv_rop_g (a_m, b_m)
         CALL source_rop_g (a_m, b_m)
         IF(call_usr_source(2)) CALL calc_usr_source (gas_continuity,&
                              a_m, b_m, lm=0)

! calculating the residual
         CALL calc_resid_c (rop_g, a_m, b_m, 0, num_resid(resid_ro,0), &
            den_resid(resid_ro,0), resid(resid_ro,0), max_resid(&
            resid_ro,0), ijk_resid(resid_ro,0))

!         call check_ab_m(a_m, b_m, 0, .true., ier)
!         call write_ab_m(a_m, b_m, ijkmax2, 0, ier)
!         write(*,*) resid(resid_ro, 0), max_resid(resid_ro, 0), &
!            ijk_resid(resid_ro, 0)
!         call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0),&
!            1, DO_K, ier)

! solving gas continuity
         CALL adjust_leq (resid(resid_ro,0), leq_it(2), leq_method(2),&
            leqi, leqm)
         CALL solve_lin_eq ('ROP_g', 2, rop_g, a_m, b_m, 0, leqi, &
            leqm, leq_sweep(2), leq_tol(2), leq_pc(2), ier)
         CALL adjust_rop (rop_g)
!        call out_array(ROP_g, 'rop_g')

      ELSE
! solve solids phase M continuity equation.

! initializing
         CALL init_ab_m (a_m, b_m, ijkmax2, m)

! forming the matrix equation
         CALL conv_rop_s (a_m, b_m, m)
         CALL source_rop_s (a_m, b_m, m)
         IF(point_source) CALL point_source_rop_s (b_m, m)
         IF(call_usr_source(2)) CALL calc_usr_source (solids_continuity, &
                              a_m, b_m, lm=m)

         CALL calc_resid_c (rop_s(1,m), a_m, b_m, m, &
            num_resid(resid_ro,m), den_resid(resid_ro,m), &
            resid(resid_ro,m), max_resid(resid_ro,m), &
            ijk_resid(resid_ro,m))
         ress = resid(resid_ro,m)

!         call check_ab_m(a_m, b_m, m, .true., ier)
!         write(*,*) 'solve_cont= ', resid(resid_ro, m),&
!                    max_resid(resid_ro, m), ijk_resid(resid_ro, m), m
!         call write_ab_m(a_m, b_m, ijkmax2, m, ier)
!         call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, M), 1, &
!                          DO_K, ier)
!
         CALL adjust_leq (resid(resid_ro,m), leq_it(2), leq_method(2),&
            leqi, leqm)
         CALL solve_lin_eq ('ROP_s', 2, rop_s(1,m), a_m, b_m, m, leqi,&
            leqm,leq_sweep(2), leq_tol(2), leq_pc(2), ier)
         CALL adjust_rop (rop_s(1,m))
!         call out_array(rop_s(1,m), 'rop_s')

      ENDIF

      call unlock_ambm
      RETURN

    CONTAINS

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: ADJUST_ROP                                              C
!  Purpose: Remove small negative values of density.                   C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE adjust_rop(ROP)

! Modules
!---------------------------------------------------------------------//
      use compar, only: ijkstart3, ijkend3
      USE functions, only: fluid_at
      USE param1, only: zero
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! density
      DOUBLE PRECISION, INTENT(INOUT) :: ROP(dimension_3)

! Local variables
!---------------------------------------------------------------------//
! Indices
      INTEGER :: IJK
!---------------------------------------------------------------------//

      DO ijk = ijkstart3, ijkend3
         IF (fluid_at(ijk)) rop(ijk) = dmax1(zero,rop(ijk))
      ENDDO

      RETURN
      END SUBROUTINE adjust_rop

      END SUBROUTINE solve_continuity

      END MODULE cont