solve_pp_g.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SOLVE_Pp_g                                              C
!  Purpose: Solve fluid pressure correction equation                   C
!                                                                      C
!  Author: M. Syamlal                                 Date: 19-JUN-96  C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE solve_pp_g(NORMG, RESG, IER)

! Modules
!---------------------------------------------------------------------//
      use ambm, only: a_m, b_m, lock_ambm, unlock_ambm
      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: zero, one, undefined
      use pgcor, only: pp_g
      use physprop, only: mmax, ro_g0
      use ps, only: point_source
      use residual, only: resid, max_resid, ijk_resid
      use residual, only: num_resid, den_resid
      use residual, only: resid_p
      use run, only: momentum_x_eq, momentum_y_eq
      use usr_src, only: call_usr_source, calc_usr_source
      use usr_src, only: pressure_correction
      IMPLICIT NONE

! Local parameters
!---------------------------------------------------------------------//
! Parameter to make tolerance for residual scaled with max value
! compatible with residual scaled with first iteration residual.
! Increase it to tighten convergence.
      DOUBLE PRECISION, PARAMETER :: DEN = 1.0d1   !5.0D2

! Dummy arguments
!---------------------------------------------------------------------//
! Normalization factor for gas pressure correction residual.
! At start of the iterate loop normg will either be 1 (i.e. not
! normalized) or a user defined value given by norm_g.  If norm_g
! was set to zero then the normalization is based on dominate
! term in the equation
      DOUBLE PRECISION, INTENT(IN) :: NORMg
! gas pressure correction residual
      DOUBLE PRECISION, INTENT(OUT) :: RESg
! Error index
      INTEGER, INTENT(INOUT) :: IER

! Local variables
!---------------------------------------------------------------------//
! phase index
      INTEGER :: M
! Normalization factor for gas pressure correction residual
      DOUBLE PRECISION :: NORMGloc
! linear equation solver method and iterations
      INTEGER :: LEQM, LEQI

! temporary use of global arrays:
! arraym1 (locally b_mmax)
! vector B_M based on dominate term in correction equation
      DOUBLE PRECISION :: B_MMAX(dimension_3, dimension_m)
! 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

! initializing
      pp_g(:) = zero
      DO m = 0, mmax
         CALL init_ab_m (a_m, b_m, ijkmax2, m)
      ENDDO

! If gas momentum equations in x and y directions are not solved return
      IF (.NOT.(momentum_x_eq(0) .OR. momentum_y_eq(0)) .AND.&
          ro_g0 .NE. undefined) THEN
        call unlock_ambm
        RETURN
      ENDIF

! Forming the sparse matrix equation.
      CALL conv_pp_g (a_m, b_m)
      CALL source_pp_g (a_m, b_m, b_mmax)
      IF(point_source) CALL point_source_pp_g (b_m, b_mmax)
      IF(call_usr_source(1)) CALL calc_usr_source(pressure_correction,&
                           a_m, b_m, lb_mmax=b_mmax, lm=0)

!      call check_ab_m(a_m, b_m, 0, .false., ier)
!      call write_ab_m(a_m, b_m, ijkmax2, 0, ier)


! Find average residual, maximum residual and location
      normgloc = normg
      IF(normg == zero) THEN
! calculating the residual based on dominate term in correction equation
! and use this to form normalization factor
        CALL calc_resid_pp (b_mmax, one, num_resid(resid_p,0), &
         den_resid(resid_p,0), resid(resid_p,0), max_resid(resid_p,0), &
         ijk_resid(resid_p,0))
         normgloc = resid(resid_p,0)/den
      ENDIF
      CALL calc_resid_pp (b_m, normgloc, num_resid(resid_p,0),  &
         den_resid(resid_p,0), resid(resid_p,0), max_resid(resid_p,0), &
         ijk_resid(resid_p,0))
      resg = resid(resid_p,0)
!      write(*,*) resid(resid_p, 0), max_resid(resid_p, 0), &
!         ijk_resid(resid_p, 0)


! Solve P_g_prime equation
       leqi = leq_it(1)
       leqm = leq_method(1)
!      CALL ADJUST_LEQ(RESID(RESID_P,0),LEQ_IT(1),LEQ_METHOD(1),LEQI,LEQM,IER)

!     call check_symmetry(A_m, 0, IER)
!     call test_lin_eq(A_M, LEQ_IT(1),LEQ_METHOD(1), LEQ_SWEEP(1), LEQ_TOL(1), LEQ_PC(1),0,IER)
      CALL solve_lin_eq ('Pp_g', 1, pp_g, a_m, b_m, 0, leqi, leqm, &
                         leq_sweep(1), leq_tol(1), leq_pc(1), ier)

!      call out_array(Pp_g, 'Pp_g')

      call unlock_ambm

      RETURN
      END SUBROUTINE solve_pp_g

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: POINT_SOURCE_Pp_g                                       C
!  Purpose: Adds point sources to the Pressure correction equation.    C
!                                                                      C
!  Notes: The off-diagonal coefficients are positive. The center       C
!         coefficient and the source vector are negative. See          C
!         conv_Pp_g                                                    C
!                                                                      C
!  Author: J. Musser                                  Date: 10-JUN-13  C
!  Reviewer:                                          Date:            C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE point_source_pp_g(B_M, B_mmax)

! Modules 
!-----------------------------------------------
      use compar, only: dead_cell_at
      use geometry, only: vol
      use functions, only: fluid_at, funijk
      use functions, only: is_on_mype_plus2layers
      use param, only: dimension_3, dimension_m
      use param1, only: small_number
      use ps, only: ps_defined, dimension_ps
      use ps, only: ps_massflow_g, ps_volume
      use ps, only: ps_k_b, ps_k_t
      use ps, only: ps_j_s, ps_j_n
      use ps, only: ps_i_w, ps_i_e
      IMPLICIT NONE

! Dummy arguments
!-----------------------------------------------
! Vector b_m
      DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
! maximum term in b_m expression
      DOUBLE PRECISION, INTENT(INOUT) :: B_mmax(dimension_3, 0:dimension_m)

! Local Variables
!-----------------------------------------------
! Indices
      INTEGER :: IJK, I, J, K
      INTEGER :: PSV

! terms of bm expression
      DOUBLE PRECISION pSource

!-----------------------------------------------
      ps_lp: do psv = 1, dimension_ps

         if(.NOT.ps_defined(psv)) cycle ps_lp
         if(ps_massflow_g(psv) < small_number) cycle ps_lp

         do k = ps_k_b(psv), ps_k_t(psv)
         do j = ps_j_s(psv), ps_j_n(psv)
         do i = ps_i_w(psv), ps_i_e(psv)

            if(.NOT.is_on_mype_plus2layers(i,j,k)) cycle
            IF (dead_cell_at(i,j,k)) cycle  ! skip dead cells

            ijk = funijk(i,j,k)
            if(fluid_at(ijk)) then
               psource = ps_massflow_g(psv) * (vol(ijk)/ps_volume(psv))

               b_m(ijk,0) = b_m(ijk,0) - psource
               b_mmax(ijk,0) = max(abs(b_mmax(ijk,0)), abs(b_m(ijk,0)))
            endif

         enddo
         enddo
         enddo

      enddo ps_lp

      RETURN
      END SUBROUTINE point_source_pp_g