hypre_obj_Setup Subroutine

    impure subroutine hypre_obj_Setup(this)
      !> Sets up the hypre objects in preparation for solves.
      ! Note: Setting up HYPRE is an expensive operation
      ! so it's best to do this only once during a simulation
      ! as opposed to setting-up and destorying each time-step.
      implicit none
      class(hypre_obj), intent(inout) :: this                                  !! Hypre machinery
      ! Work variables
      integer :: i,j,k
      real(wp):: dxi,dxmi(2)
      real(wp):: dyi,dymi(2)
      real(wp):: dzi,dzmi(2)

      if (.not.allocated(this%mat)) then
        ! Use default operator: 2nd order Laplacian
        ! Set stencil extent
        this%stm = 1

        ! Set matrix coefficients
        associate (x=>this%block%x, y=>this%block%y, z=>this%block%z,   &
          xm=>this%block%xm, ym=>this%block%ym, zm=>this%block%zm,      &
          lo=>this%block%lo, hi =>this%block%hi)
          allocate(this%mat(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3),3,-this%stm:this%stm))

          do k=lo(3),hi(3)
            do j=lo(2),hi(2)
              do i=lo(1),hi(1)
                dxi     = 1.0_wp/(x(i+1)-x(i))
                dxmi(1) = 1.0_wp/(xm(i) -xm(i-1))
                dxmi(2) = 1.0_wp/(xm(i+1) -xm(i))

                dyi     = 1.0_wp/(y(j+1)-y(j))
                dymi(1) = 1.0_wp/(ym(j) -ym(j-1))
                dymi(2) = 1.0_wp/(ym(j+1) -ym(j))

                dzi     = 1.0_wp/(z(k+1)-z(k))
                dzmi(1) = 1.0_wp/(zm(k) -zm(k-1))
                dzmi(2) = 1.0_wp/(zm(k+1) -zm(k))
                ! Second derivative with respect to x1 (d2dx1): lower/diagonal/upper
                this%mat(i,j,k,1,:) = [dxi*dxmi(1),dxi*(-dxmi(1))-dxi*dxmi(2),dxi*dxmi(2)]
                ! Second derivative with respect to x2 (d2dx2): lower/diagonal/upper
                this%mat(i,j,k,2,:) = [dyi*dymi(1),dyi*(-dymi(1))-dyi*dymi(2),dyi*dymi(2)]
                ! Second derivative with respect to x3 (d2dx3): lower/diagonal/upper
                this%mat(i,j,k,3,:) = [dzi*dzmi(1),dzi*(-dzmi(1))-dzi*dzmi(2),dzi*dzmi(2)]
              end do
            end do
          end do
        end associate
      end if

      select case (this%interface)
      case (HYPRE_INT_Struct)
        ! Initialize the hypre components
        call hypre_obj_SetupGridS(this)
        call hypre_obj_SetupStencilS(this)
        call hypre_obj_SetupRHSS(this)
        call hypre_obj_SetupSolS(this)
        call hypre_obj_SetupMatrixS(this)

        ! Now setup the sovler
        call hypre_obj_SetupSolverS(this)
      case (HYPRE_INT_SStruct)
        ! Not yet implemented
      case (HYPRE_INT_IJ)
        call hypre_obj_SetupPointersIJ(this)
        call hypre_obj_SetupRowsIJ(this)
        call hypre_obj_SetupMatrixIJ(this)
        call hypre_obj_SetupRHSIJ(this)
        call hypre_obj_SetupSolIJ(this)

        ! Now setup the sovler
        call hypre_obj_SetupSolIJver(this)
      case default
        call this%parallel%Stop("Unknown or uninitialized HYPRE interface")
      end select

      ! Clear internal matrix
      deallocate(this%mat)

      return
    end subroutine hypre_obj_Setup