H.1.5 Add_to_Matrix_Equation_Monomial Procedure

The main documentation of the Add_to_Matrix_Equation_Monomial Procedure contains additional explanation of this code listing.

  subroutine Add_to_Matrix_Equation_Monomial (Monomial, ELLM, RHS_MV)

    ! Use association information.

    use Caesar_Numbers_Module, only: one, zero

    ! Input variables.

    type(Monomial_type), intent(inout) :: Monomial  ! Monomial to be added.

    ! Input/Output variables.

    ! ELL_Matrix to be incremented.
    type(ELL_Matrix_type), intent(inout) :: ELLM
    ! RHS mathematic vector.
    type(Mathematic_Vector_type), intent(inout) :: RHS_MV   

    ! Internal variables.

    ! Arrays to manipulate values of the matrix.
    type(real,1) :: Matrix_Values   
    type(integer,1) :: Matrix_Rows
    type(integer,1) :: Matrix_Columns
    type(integer) :: i                          ! Loop parameter.
    type(integer) :: shift                      ! Index shift.
    type(real,1) :: Volume_Cells                ! Volume of the cells.

    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    ! Verify requirements.

    VERIFY(Valid_State(Monomial),5) ! Monomial is valid.
    VERIFY(Valid_State(ELLM),5)     ! ELLM is valid.
    VERIFY(Valid_State(RHS_MV),5)   ! RHS_MV is valid.

    ! Create some working vectors.

    call Initialize (Matrix_Values, NCells_PE(Monomial%Mesh))
    call Initialize (Matrix_Rows, NCells_PE(Monomial%Mesh))
    call Initialize (Matrix_Columns, NCells_PE(Monomial%Mesh))
    call Initialize (Volume_Cells, NCells_PE(Monomial%Mesh))

    ! Set Matrix_Rows, Matrix_Columns.
    !
    ! To spread out the rows and columns from this equation into a multiple
    ! equation matrix, translate local row numbers via:
    !   i --> NEquations*(i - 1) + Equation
    ! and global column numbers via:
    !   i + shift --> NEquations*(i + shift - 1) + Variable

    shift = First_PE(Monomial%Structure) - 1
    do i = 1, Length_PE(Monomial%Structure)
      Matrix_Rows(i) = Monomial%NEquations*(i - 1) + Monomial%Equation
      Matrix_Columns(i) = Monomial%NEquations*(i + shift - 1) + &
                          Monomial%Variable
    end do

    ! Increment ELLM values.

    call Get_Volume_Cells (Volume_Cells, Monomial%Mesh)
    if (Monomial%Exponent /= zero) then
      if (Monomial%Derivative) then
        ! Time derivative version.
        Matrix_Values(:) = Monomial%Coefficient * Monomial%Exponent * &
                           Monomial%Phi_k ** (Monomial%Exponent - one) * &
                           Volume_Cells / Monomial%Delta_t
      else
        ! Standard version.
        Matrix_Values(:) = Monomial%Coefficient * Monomial%Exponent * &
                           Monomial%Phi_k ** (Monomial%Exponent - one) * &
                           Volume_Cells
      end if
    else
      ! Treat exponent=0 (monomial=constant) case separately to avoid 
      ! division by zero if Phi_k=0 and to save a little time.
      Matrix_Values(:) = zero
    end if
    call Add_Values (ELLM, Matrix_Values, Matrix_Rows, Matrix_Columns, &
                     Global=.false.)

    ! Increment RHS_MV values. Note that they have been multiplied by -1
    ! to move them to the RHS of the equation.

    if (Monomial%Derivative) then
      ! Time derivative version.
      Matrix_Values(:) = Monomial%Coefficient *  &
                         ((Monomial%Exponent - one) * &
                         Monomial%Phi_k ** Monomial%Exponent + &
                         Monomial%Phi_n ** Monomial%Exponent) * &
                         Volume_Cells / Monomial%Delta_t
    else
      ! Standard version.
      Matrix_Values(:) = Monomial%Coefficient * &
                         (Monomial%Exponent - one) * &
                         Monomial%Phi_k ** Monomial%Exponent * &
                         Volume_Cells
    end if
    call Add_Values (RHS_MV, Matrix_Values, Matrix_Rows, Global=.false.)

    ! Finalize working vectors.

    call Finalize (Matrix_Values)
    call Finalize (Matrix_Rows)
    call Finalize (Matrix_Columns)
    call Finalize (Volume_Cells)

    ! Verify guarantees.

    VERIFY(Valid_State(ELLM),5)     ! ELLM is valid.
    VERIFY(Valid_State(RHS_MV),5)   ! RHS_MV is valid.

    return
  end subroutine Add_to_Matrix_Equation_Monomial