Jacobi method converging then diverging

Question

I am working to solve Poisson's equation (in 2d axisymmetric cylindrical coordinates) using the Jacobi method. The L2 norm decreases from ~1E3 on the first iteration (I have a really bad guess) to ~0.2 very slowly. Then, the L2 norm begins to increase over many iterations.

My geometry is parallel plates with sharp points at r = 0 on both plates. (If that matters).

Is there some error in my code? Do I need to go to a different algorithm? (I have a not yet working DADI algorithm.)

Here is my Jacobi method algorithm. Then this just in wrapped in a while loop.

subroutine Jacobi(PoissonRHS, V, resid)
implicit none
real, dimension(0:,0:) :: PoissonRHS, V
REAL resid
integer :: i,j, lb, ub
real, dimension(0:size(V,1)-1, 0:size(V,2)-1) :: oldV

real :: dr = delta(1)
real :: dz = delta(2)

real :: dr2 = (delta(1))**(-2)
real :: dz2 = (delta(2))**(-2)

integer :: M = cells(1)
integer :: N = cells(2)


oldV = V
!Note: All of the equations are second order accurate

!If at r = 0 and in the computational domain
! This is the smoothness condition, dV(r=0)/dr = 0
    V(0,:) = (4.0*oldV(1,:)-oldV(2,:))/3.0

!If at r = rMax and in the computational domain
! This is an approximation and should be fixed to improve accuracy, it should be
! lim r->inf V' = 0, while this is V'(r = R) = 0
    V(M, 1:N-1) = 0.5 / (dr2 + dz2) * (     &
        (2.0*dr2)*oldV(M-1,1:N-1) +         &
        dz2 * (oldV(M,2:N) + oldV(M,0:N-2)) &
        - PoissonRHS(M,1:N-1))

do i = 1, M-1
    lb = max(0, nint(lowerBoundary(i * dr) / dz)) + 1
    ub = min(N, nint(upperBoundary(i * dr) / dz)) - 1

    V(i,lb:ub) =  0.5 / (dr2 + dz2) * ( &
        ((1.0 - 0.5/dble(i))*dr2)*oldV(i-1,lb:ub) + &
        ((1.0 + 0.5/dble(i))*dr2)*oldV(i+1,lb:ub) + &
        dz2 * (oldV(i,lb+1:ub+1) + oldV(i,lb-1:ub-1)) &
        - PoissonRHS(i,lb:ub))

    V(i, 0:lb-1) = V0
    V(i, ub+1:N) = VL
enddo

!compare to old V values to check for convergence
resid = sqrt(sum((oldV-V)**2))

return
end subroutine Jacobi

Answer 1

Based on additional readings it seems like it was a precision problem. Because (for example), I had the expression

    V(i,lb:ub) =  0.5 / (dr2 + dz2) * ( &
        ((1.0 - 0.5/dble(i))*dr2)*oldV(i-1,lb:ub) + &
        ((1.0 + 0.5/dble(i))*dr2)*oldV(i+1,lb:ub) + &
        dz2 * (oldV(i,lb+1:ub+1) + oldV(i,lb-1:ub-1)) &
        - PoissonRHS(i,lb:ub))

where dr2 and dz2 are very large. So by distributing these I got terms that were ~1 and the code converges (slowly, but that's a function of the mathematics).

So my new code is

subroutine Preconditioned_Jacobi(PoissonRHS, V, resid)
implicit none
real, dimension(0:,0:) :: PoissonRHS, V
REAL resid
integer :: i,j, lb, ub
real, dimension(0:size(V,1)-1, 0:size(V,2)-1) :: oldV

real :: dr = delta(1)
real :: dz = delta(2)

real :: dr2 = (delta(1))**(-2)
real :: dz2 = (delta(2))**(-2)

real :: b,c,d

integer :: M = cells(1)
integer :: N = cells(2)

    b = 0.5*(dr**2)/((dr**2) + (dz**2))
    c = 0.5*(dz**2)/((dr**2) + (dz**2))
    d = -0.5 / (dr2 + dz2)

oldV = V
!Note: All of the equations are second order accurate

!If at r = 0 and in the computational domain
! This is the smoothness condition, dV(r=0)/dr = 0
    V(0,:) = (4.0*oldV(1,:)-oldV(2,:))/3.0 !same as: oldV(0,:) - 2.0/3.0 * (1.5 * oldV(0,:) - 2.0 * oldV(1,:) + 0.5 * oldV(2,:) - 0)

!If at r = rMax and in the computational domain
! This is an approximation and should be fixed to improve accuracy, it should be
! lim r->inf V' = 0, while this is V'(r = R) = 0
V(M,1:N-1) = d*PoissonRHS(M,1:N-1) &
                + 2.0*c * oldV(M-1,1:N-1) &
                + b * ( oldV(M,0:N) + oldV(M,2:N) )

do i = 1, M-1
    lb = max(0, nint(lowerBoundary(i * dr) / dz)) + 1
    ub = min(N, nint(upperBoundary(i * dr) / dz)) - 1

    V(i,lb:ub) = d*PoissonRHS(i,lb:ub)  &
            + (c * (1.0-0.5/dble(i)) * oldV(i-1,lb:ub)) &
            + (c * (1.0+0.5/dble(i)) * oldV(i+1,lb:ub)) &
            + b * (oldV(i,lb-1:ub-1) + oldV(i,lb+1:ub+1))

    V(i, 0:lb-1) = V0
    V(i, ub+1:N) = VL
enddo

!compare to old V values to check for convergence
resid = sum(abs(oldV-V))

return
end subroutine Preconditioned_Jacobi