Jacobi方法收敛然后发散

Question

我正在用Jacobi方法求解泊松方程(在2d轴对称圆柱坐标系中)。L2范数从第一次迭代时的~1E3 (我有一个非常糟糕的猜测)下降到~0.2，非常慢。然后，L2范数开始在多次迭代中增加。

我的几何学是平行的板，在两个板上r=0处有尖锐的点。(如果这重要的话)。

我的代码中有错误吗？我需要使用不同的算法吗？(我有一个尚未工作的DADI算法。)

这是我的Jacobi方法算法。然后把它封装在一个while循环中。

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

根据额外的读数，这似乎是一个精确的问题。因为(例如)，我有一个表达式

    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))

其中dr2和dz2非常大。因此，通过分发这些项，我得到了~1的项和代码的收敛(缓慢，但这是数学的函数)。

所以我的新代码是

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