Issue solving 1D Poisson Eq.

[markowkes]

I'm trying to use IterativeSolvers.jl but am struggling to get a very simple example working. This code solves L*P=R where L is a Laplacian operator (discretized with a finite difference) and R is zero. The first and last equations set the boundary condition of P=10 and P=100, respectively. I'm solving the system with conjugate gradient, BiCGStab, IDRS, and the "\". All the methods are not working except the base Julia operator "\". Any information on why this is not working would be very helpful.

using IterativeSolvers
using Plots

# Create grid
Nx=100
x=range(0,Nx,Nx)
dx=x[2]-x[1]

# Create Laplacian - eventually I want to do this with LinearMap.jl
L = zeros(Nx,Nx)
L[ 1, 1] = 1.0  # Left BC: Dirchlet
L[Nx,Nx] = 1.0  # Right BC: Dirchlet 
# Interior points: Laplacian 
for i=2:Nx-1
    L[i,i  ] = -2/dx^2
    L[i,i-1] =  1/dx^2
    L[i,i+1] =  1/dx^2
end

# Define RHS 
R=zeros(Nx)
R[ 1] = 10  # Left BC 
R[Nx] = 100 # Right BC

# Solve for p using various methods
P1 =        cg(L, R; abstol=1e-6,verbose=true)
P2 = bicgstabl(L, R; abstol=1e-6,verbose=true)
P3 =      idrs(L, R, abstol=1e-6,verbose=true)
P4 = L\R 

# Plot results
plt1 = plot(x,P1,title="Conjugate Gradient",legend=false)
plt2 = plot(x,P2,title="BiCGStab"          ,legend=false)
plt3 = plot(x,P3,title="IDRS"              ,legend=false)
plt4 = plot(x,P4,title="Julia"             ,legend=false)
plot(plt1,plt2,plt3,plt4,layout = 4)

[markowkes]

I realized the issue was that the initial guess did not include the Dirichlet boundary conditions. FYI, here's a code that implements the Poisson equation with Dirchlet boundary conditions and shows how it can be solved with a LinearMap, IterativeSolvers, and KrylovKit.

using LinearMaps
import IterativeSolvers
import KrylovKit
using Plots

# Create grid
Nx=100
x=range(0,10,Nx)
dx=x[2]-x[1]

# Solve L*p = R 

# # Create matrix - eventually I want to do this with LinearMap.jl
# L = zeros(Nx,Nx)
# L[ 1, 1] = 1.0  # Left BC: Dirchlet
# L[Nx,Nx] = 1.0  # Right BC: Dirchlet 
# # Interior points: Laplacian 
# for i=2:Nx-1
#     L[i,i  ] =  2/dx^2
#     L[i,i-1] = -1/dx^2
#     L[i,i+1] = -1/dx^2
# end

"""
Define L operator
    - Input:  P 
    - Output: Lp = L*p
"""  
function L!(Lp::AbstractVector,P::AbstractVector)

    #println("P=",P)
    # Left BC: Dirchlet
    Lp[1] = P[1]

    # Right BC: Dirchlet 
    Lp[Nx] = P[Nx] 

    # Interior points: Laplacian 
    for i=2:Nx-1
        Lp[i] = (P[i-1] - 2*P[i] + P[i+1])/dx^2
    end
    return nothing
end

# Create linear map to let L! act as a matrix 
L = LinearMap(L!,Nx,ismutating=true)

# Define RHS 
R=ones(Nx)
R[ 1] = 10  # Left BC 
R[Nx] = 100 # Right BC

# Solve for p using various methods
Po = zeros(Nx); Po[1]=R[1]; Po[Nx]=R[Nx]
P1=copy(Po); P1 =        IterativeSolvers.cg!(P1, L, R; abstol=1e-6,verbose=true)
P2=copy(Po); P2 = IterativeSolvers.bicgstabl!(P2, L, R; abstol=1e-6,verbose=true,max_mv_products=200)
P3=copy(Po); P3 =      IterativeSolvers.idrs!(P3, L, R, abstol=1e-6,verbose=true,maxiter=200)
#P4=copy(Po); P4 = L\R  # Doesn't work with LinearMap
P5=copy(Po); P5 = KrylovKit.linsolve(L,R)
P6=copy(Po); P6 = KrylovKit.linsolve(L,R,P6,KrylovKit.CG())
P7=copy(Po); P7 = KrylovKit.linsolve(L,R,P7,KrylovKit.BiCGStab(maxiter=1000))
P8=copy(Po); P8 = KrylovKit.linsolve(L,R,P8,KrylovKit.GMRES())

plt1 = plot(x,P1,title="Conjugate Gradient",legend=false); 
plt2 = plot(x,P2,title="BiCGStab"          ,legend=false); 
plt3 = plot(x,P3,title="IDRS"              ,legend=false); 
#plt4 = plot(x,P4,title="Julia"             ,legend=false);
plt5 = plot(x,P5,title="KrylovKit.linsolve",legend=false); 
plt6 = plot(x,P6,title="KrylovKit.CG"      ,legend=false); 
plt7 = plot(x,P7,title="KrylovKit.BiCGStab",legend=false); 
plt8 = plot(x,P8,title="KrylovKit.GMRES"   ,legend=false); 
plot(plt1,plt2,plt3,plt5,plt6,plt7,plt8)