algorithm - 雅可比方法先收敛再发散

Question

我正在使用 Jacobi 方法求解泊松方程(在二维轴对称圆柱坐标系中)。 L2 范数从第一次迭代的 ~1E3(我的猜测非常糟糕)非常缓慢地减少到 ~0.2。然后，经过多次迭代，L2 范数开始增加。

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

我的代码有什么错误吗？我需要使用不同的算法吗？ (我有一个尚未运行的 DADI 算法。)

这是我的雅可比法算法。然后这只是包裹在一个 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