jacobian_block_get_tridiag Subroutine

  recursive subroutine jacobian_block_get_tridiag(this, diagonal, subdiagonal, &
                                                  superdiagonal)
    !* Author: Chris MacMackin
    !  Date: May 2017
    !
    ! Computes the tridiagonal matrix used to solve for this Jacobian block.
    !
    class(jacobian_block), intent(in)                :: this
    real(r8), dimension(:), allocatable, intent(out) :: diagonal
    real(r8), dimension(:), allocatable, intent(out) :: subdiagonal
    real(r8), dimension(:), allocatable, intent(out) :: superdiagonal

    integer                                   :: n, i, pos
    class(vector_field), pointer              :: grid
    logical                                   :: use_cached
    class(scalar_field), allocatable, save    :: cached_field_type
    real(r8), dimension(:), allocatable, save :: cached_dx_c,   &
                                                 cached_dx2_uc, &
                                                 cached_dx2_ud, &
                                                 cached_dx2_dc
    real(r8), save                            :: cached_upper_bound, &
                                                 cached_lower_bound
    real(r8), dimension(:), allocatable       :: cont, deriv, diag, &
                                                 subdiag, supdiag
    real(r8), dimension(:,:), allocatable     :: domain

    n = this%contents%raw_size()
    domain = this%contents%domain()
    ! Try to use cached copy of inverse grid spacing, if available and suitable 
    if (allocated(cached_field_type)) then
      use_cached = same_type_as(this%contents, cached_field_type) .and. &
                   abs(domain(1,1) - cached_lower_bound) < 1e-10 .and.  &
                   abs(domain(1,2) - cached_upper_bound) < 1e-10 .and.  &
                   n==size(cached_dx_c)
    else
      use_cached = .false.
    end if
    ! Construct an array containing the distance between grid
    ! points, if necessary
    if (.not. use_cached) then
      if (allocated(cached_field_type)) then   
        deallocate(cached_field_type)
        deallocate(cached_dx_c)
        deallocate(cached_dx2_uc)
        deallocate(cached_dx2_dc)
        deallocate(cached_dx2_ud)
      end if
#ifdef DEBUG
      call logger%debug('jacobian_block%get_tridiag','Calculating and caching '// &
                       'grid spacings.')
#endif
      cached_lower_bound = domain(1,1)
      cached_upper_bound = domain(1,2)
      allocate(cached_field_type, mold=this%contents)
      allocate(cached_dx_c(n))
      allocate(cached_dx2_uc(n-1))
      allocate(cached_dx2_dc(n-1))
      allocate(cached_dx2_ud(n))
      call this%contents%allocate_vector_field(grid)
      grid = this%contents%grid_spacing()
      ! Use this array to temporarily hold the grid spacings
      cached_dx_c = grid%raw()
      ! Work out the upwinded grid spacings
      cached_dx2_uc(1) = cached_dx_c(1)
      do i = 2, n-1
        cached_dx2_uc(i) = 2*cached_dx_c(i) - cached_dx2_uc(i-1)
      end do
      ! Work out the inverse of the downwinded grid spacings
      cached_dx2_dc(1:n-2) = 1._r8/cached_dx2_uc(1:n-2)
      cached_dx2_dc(n-1) = 1._r8/cached_dx_c(n)
      ! Work out the -2 times product of the inverse up- and
      ! down-winded grid spacings
      cached_dx2_ud(1) = cached_dx2_uc(1)**(-2)
      cached_dx2_ud(2:n-1) = -2 * cached_dx2_dc(2:n-1)/cached_dx2_uc(1:n-2)
      cached_dx2_ud(n) = cached_dx2_dc(n-1)**2
      ! Work out the inverse of product of the up-winded and
      ! centred grid spacings
      cached_dx2_uc(1) = -cached_dx2_uc(1)**(-2)
      cached_dx2_uc(2:n-1) = 1._r8/(cached_dx2_uc(2:n-1) * cached_dx_c(2:n-1))
      ! Work out the inverse of the product of the down-winded and
      ! centred grid spacings
      cached_dx2_dc(1:n-2) = cached_dx2_dc(1:n-2) / cached_dx_c(2:n-1)
      cached_dx2_dc(n-1) = -cached_dx2_dc(n-1)**2
      ! Work out half the inverse of the grid spacing
      cached_dx_c(2:n-1) = 0.5_r8/cached_dx_c(2:n-1)
      cached_dx_c(1) = 1._r8/cached_dx_c(1)
      cached_dx_c(n) = 1._r8/cached_dx_c(n)
    end if
    cont = this%contents%raw()
    if (this%extra_derivative == no_extra_derivative) then
      ! Create the tridiagonal matrix when there is no additional
      ! differentiation operator on the RHS
      diagonal = this%coef*this%derivative%raw()
      select case(this%has_increment)
      case(0)
        continue
      case(1)
        diagonal = diagonal + this%real_increment
      case(2)
        diagonal = diagonal + this%field_increment%raw()
      case default
        error stop ('Invalid increment has been added.')
      end select
      diagonal(1) = diagonal(1) - cont(1)*cached_dx_c(1)
      diagonal(n) = diagonal(n) + cont(n)*cached_dx_c(n)
      superdiagonal = cont(1:n-1) * cached_dx_c(1:n-1)
      subdiagonal = -cont(2:n) * cached_dx_c(2:n)
    else if (this%extra_derivative == this%direction) then
      ! Create the tridiagonal matrix when the additional
      ! differentiation operator on the RHS operates in the same
      ! direction as the derivative of the contents.
      deriv = this%coef*this%derivative%raw()
      allocate(diagonal(n))
      diagonal(2:n-1) = cont(2:n-1) * cached_dx2_ud(2:n-1)
      diagonal(1) = cont(1) * cached_dx2_ud(1) &
                       - deriv(1) * cached_dx_c(1)
      diagonal(n) = cont(n) * cached_dx2_ud(n) &
                       + deriv(n) * cached_dx_c(n)
      !print*,diagonal(10)
      select case(this%has_increment)
      case(0)
        continue
      case(1)
        diagonal = diagonal + this%real_increment
      case(2)
        diagonal = diagonal + this%field_increment%raw()
      case default
        error stop ('Invalid increment has been added.')
      end select
      allocate(superdiagonal(n-1))
      superdiagonal(2:n-1) = cont(2:n-1) * cached_dx2_uc(2:n-1) &
                           + deriv(2:n-1) * cached_dx_c(2:n-1)
      superdiagonal(1) = -2 * cont(1) * cached_dx2_uc(1) &
                       + deriv(1) * cached_dx_c(1)
      allocate(subdiagonal(n-1))
      subdiagonal(1:n-2) = cont(2:n-1) * cached_dx2_dc(1:n-2) &
                         - deriv(2:n-1) * cached_dx_c(2:n-1)
      subdiagonal(n-1) = -2 * cont(n) * cached_dx2_dc(n-1) &
                       - deriv(n) * cached_dx_c(n)
    else
      call logger%fatal('jacobian_block%solve_for', 'Currently no support '// &
                        'for extra derivatives other than in the 1-direction')
      error stop
    end if
    if (associated(this%block_increment)) then
      call this%block_increment%get_tridiag(diag, subdiag, supdiag)
      diagonal = diagonal + diag
      subdiagonal = subdiagonal + subdiag
      superdiagonal = superdiagonal + supdiag
    end if
    ! Set the boundary conditions for this problem
    do i = 1, size(this%boundary_locs)
      pos = this%boundary_locs(i)
      if (this%boundary_types(i) == dirichlet) then
        diagonal(pos) = 1._r8
        if (pos < n) superdiagonal(pos) = 0._r8
        if (pos > 1) subdiagonal(pos-1) = 0._r8
      else if (this%boundary_types(i) == neumann) then
        if (pos == n) then
          diagonal(n) = cached_dx_c(n)
          subdiagonal(n-1) = -cached_dx_c(n)
        else if (pos == 1) then
          diagonal(1) = -cached_dx_c(1)
          superdiagonal(1) = cached_dx_c(1)
        else
          superdiagonal(pos) = cached_dx_c(pos)
          subdiagonal(pos-1) = cached_dx_c(pos)
          diagonal(pos) = 0._r8
        end if
      else
        call logger%fatal('jacobian_block%solve_for','Boundary condition of '// &
                          'type other than Dirichlet or Neumann encountered.')
        error stop
      end if
    end do
#ifdef DEBUG
    call logger%debug('jacobian_block%get_tridiag','Constructed tridiagonal matrix.')
#endif
  end subroutine jacobian_block_get_tridiag