test_lin_eq.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: TEST_LIN_EQ(A_m, LEQIT, LEQMETHOD, LEQSWEEP, LEQTOL, LEQPC, TEST, IER)
!  Purpose: Routine for testing the accuracy of linear equation solver C
!                                                                      C
!                                                                      C
!  Author: M. Syamlal                                 Date: 4-JUN-96   C
!  Reviewer:                                          Date:            C
!                                                                      C
!  **** See Solve_energy_eq.f for an example  ****                     C
!
!  Literature/Document References:                                     C
!                                                                      C
!  Variables referenced:                                               C
!  Variables modified:                                                 C
!                                                                      C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
!
!
      SUBROUTINE test_lin_eq(A_M, LEQIT, LEQMETHOD, LEQSWEEP, LEQTOL, LEQPC,  TEST, IER)
!...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
!...Switches: -xf
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE param
      USE param1
      USE geometry
      USE indices
      USE compar
      USE functions
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER TEST  !=0 use the passed A_m; =1 construct a random A_m
      INTEGER IER
      DOUBLE PRECISION, DIMENSION(DIMENSION_3,-3:3) :: A_M
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: IJKERR
      INTEGER IJK, IpJK, ImJK, IJpK, IJmK, IJKp, IJKm
      DOUBLE PRECISION, DIMENSION(DIMENSION_3,-3:3) :: Am
      DOUBLE PRECISION, DIMENSION(DIMENSION_3) :: Bm, X_ACT, X_SOL
      DOUBLE PRECISION :: ERR, ERRMAX, ERRSUM, XSUM
      CHARACTER(LEN=80), DIMENSION(8) :: LINE
!
!                      linear equation solver method and iterations
      INTEGER          LEQMETHOD, LEQIT
      CHARACTER(LEN=4) ::   LEQSWEEP
      DOUBLE PRECISION LEQTOL
      CHARACTER(LEN=4) ::   LEQPC
      REAL  :: Harvest

!  Initialize the random number generator
!
      CALL random_seed
!
!  Fill the A and x arrays with random numbers, but ensuring that
!  the matrix is diagonally dominant
!
      DO ijk = ijkstart3, ijkend3
         CALL random_number(harvest)
         x_act(ijk) = dble(harvest) + 1.e-5
         x_sol(ijk) = 0.0
         IF (test == 0) THEN
            am(ijk,-3) = a_m(ijk,-3)
            am(ijk,-2) = a_m(ijk,-2)
            am(ijk,-1) = a_m(ijk,-1)
            am(ijk,0) = a_m(ijk,0)
            am(ijk,1) = a_m(ijk,1)
            am(ijk,2) = a_m(ijk,2)
            am(ijk,3) = a_m(ijk,3)
         ELSE
            CALL random_number(harvest)
            am(ijk,-3) = dble(harvest)
            CALL random_number(harvest)
            am(ijk,-2) = dble(harvest)
            CALL random_number(harvest)
            am(ijk,-1) = dble(harvest)
            CALL random_number(harvest)
            am(ijk,0) = -dble(max(harvest,0.1))*70.
            CALL random_number(harvest)
            am(ijk,1) = dble(harvest)
            CALL random_number(harvest)
            am(ijk,2) = dble(harvest)
            CALL random_number(harvest)
            am(ijk,3) = dble(harvest)
         ENDIF
      END DO


!
!!!$omp  parallel do private( IJK, IpJK, ImJK, IJpK, IJmK, IJKp, IJKm)
      DO ijk = ijkstart3, ijkend3

            imjk = im_of(ijk)
            ijmk = jm_of(ijk)
            ipjk = ip_of(ijk)
            ijpk = jp_of(ijk)
            bm(ijk) = am(ijk,0)*x_act(ijk)
            IF(i_of(ijk) > 1) bm(ijk) = bm(ijk) + am(ijk,west)*x_act(imjk)
            IF(i_of(ijk) < imax2) bm(ijk) = bm(ijk) +am(ijk,east)*x_act(ipjk)
            IF(j_of(ijk) > 1) bm(ijk) = bm(ijk) + am(ijk,south)*x_act(ijmk)
            IF(j_of(ijk) < jmax2) bm(ijk) = bm(ijk) + am(ijk,north)*x_act(ijpk)
            IF (do_k) THEN
               ijkm = km_of(ijk)
               ijkp = kp_of(ijk)
               IF(k_of(ijk) > 1) bm(ijk) = bm(ijk) + am(ijk,bottom)*x_act(ijkm)
               IF(k_of(ijk) < kmax2) bm(ijk) = bm(ijk) + am(ijk,top)*x_act(ijkp)
            ENDIF
      END DO

!
!  Solve the linear equation
!
      CALL solve_lin_eq ('Test', 1, x_sol, am, bm, 0, leqit, leqmethod, leqsweep, leqtol, leqpc, ier)
!
!  Check the solution
!
      errsum = 0.0
      xsum = 0.0
      errmax = 0.0
      ijkerr = 0
      DO ijk = ijkstart3, ijkend3
         IF (x_act(ijk) /= 0.0) THEN
            err = abs(x_sol(ijk)-x_act(ijk))/x_act(ijk)
         ELSE IF (x_sol(ijk) == 0.0) THEN
            err = 0.0
         ELSE
            err = 1.0d32
         ENDIF
         IF (err > errmax) THEN
            errmax = err
            ijkerr = ijk
         ENDIF
         errsum = errsum + abs(x_sol(ijk)-x_act(ijk))
         xsum = xsum + abs(x_act(ijk))
!        print *, ijk, i_of(ijk), j_of(ijk), nint(err * 100.)
      END DO
      IF (xsum /= 0.0) THEN
         err = errsum/xsum
      ELSE IF (errsum == 0.0) THEN
         err = 0.0
      ELSE
         err = 1.0d32
      ENDIF
!
      IF (err < leqtol) THEN
         ier = 0
         line(1) = 'Message: Lin equation solution satisfies tolerance.'
      ELSE
         ier = 1
         line(1) = 'Error: Lin equation solution does not satisfy tolerance!'
      ENDIF
!
      WRITE (line(2), *) 'Average normalized error = ', err
      WRITE (line(3), *) 'Max normalized error = ', errmax
      WRITE (line(4), *) 'Location of max error = ', ijkerr
      WRITE (line(5), *) 'Xa and Xs @ max error = ', x_act(ijkerr), x_sol(ijkerr)
      WRITE (line(6), *) 'Sample values of actual (Xa) and solution (Xs):'
      WRITE (line(7), '(A,G12.5, A, I6, A, G12.5, A, I6, A, G12.5)') 'Xa(1)=', &
         x_act(1), '  Xa(', ijkmax2/2, ')=', x_act(ijkmax2/2), '  Xa(', ijkmax2, ')=', &
         x_act(ijkmax2)
      WRITE (line(8), '(A,G12.5, A, I6, A, G12.5, A, I6, A, G12.5)') 'Xs(1)=', &
         x_sol(1), '  Xs(', ijkmax2/2, ')=', x_sol(ijkmax2/2), '  Xs(', ijkmax2&
         , ')=', x_sol(ijkmax2)
!
      CALL write_error ('TEST_LIN_EQ', line, 8)
!
      RETURN
      END SUBROUTINE test_lin_eq