Solving Poisson PDE in an arbitrary 2D geometry, under different types of both internal and external boundary conditions

Specific Domains Numerics
1

Fairly new to Julia, and I am trying to solve a 2D PDE problem of the form:
image
with finite differences where the 2D domain may not necessarily be a pure rectangle (e.g., a cross-like shape, or rectangle with internal rectangular holes [one or more]). In addition, the problem should consider Dirichlet (external and/or internal), Neumann (external and/or internal), and periodic (only external) boundary conditions. The exact specification of the geometry and the type of boundary conditions are typically provided with an integer matrix with following notation:

  • 1: simulation region (R)
  • 0: outside-of-simulation region (E)
  • -1: Dirichlet boundary (D)
  • -2: Neumann boundary (N)
  • -3: periodic boundary (P).

This matrix usually has 2 extra rows and 2 extra columns to include the ghost regions for NBC and/or PBC. For instance, when I have all DBC on the external region (fully rectangular, no holes):

                  E E E E E E E E E E
                  E D D D D D D D D E
                  E D R R R R R R D E
                  E D R R R R R R D E
lattice_matrix =  E D R R R R R R D E
                  E D R R R R R R D E
                  E D R R R R R R D E
                  E D D D D D D D D E
                  E E E E E E E E E E

In case of all NBC (or PBC) on the external region (fully rectangular, no holes):

                  N N N N N N N N N
                  N R R R R R R R N
                  N R R R R R R R N
                  N R R R R R R R N
                  N R R R R R R R N
lattice_matrix =  N R R R R R R R N
                  N R R R R R R R N
                  N R R R R R R R N
                  N R R R R R R R N
                  N R R R R R R R N
                  N N N N N N N N N

the actual value of the D or N boundary conditions are specified in a separate matrix and are subsequently included in the right-hand side of the equation.
Equation is only solved at the R points.

I already tried to write some code to do this, but I am getting somewhat strange results with it (see below).
Could somebody help me a bit with this approach and let me know what am I doing wrong here?

Thanks in advance,
Željko

const REG::Int64 = 1;
const EXT::Int64 = 0;
const DBC::Int64 = -1;
const NBC::Int64 = -2;
const PBC::Int64 = -3;


function rectangularLattice(N::Int64, M::Int64, BC::String)
    
    lattice = ones(Int64, N+2, M+2).*EXT;
    lattice[2:N+1, 2:M+1] .= REG;

    if BC[1] == 'D'
        lattice[2:end-1, 2] .= DBC;
    elseif BC[1] == 'N'
        lattice[1:end, 1] .= NBC;
    elseif BC[1] == 'P'
        lattice[1:end, 1] .= PBC;
    end

    if BC[2] == 'D'
        lattice[2:end-1, end-1] .= DBC;
    elseif BC[2] == 'N'
        lattice[1:end, end] .= NBC;
    elseif BC[2] == 'P'
        lattice[1:end, end] .= PBC;
    end

    if BC[3] == 'D'
        lattice[2, 2:end-1] .= DBC;
    elseif BC[3] == 'N'
        lattice[1, 1:end] .= NBC;
    elseif BC[3] == 'P'
        lattice[1, 1:end] .= PBC;
    end
    
    if BC[4] == 'D'
        lattice[end-1, 2:end-1] .= DBC;
    elseif BC[4] == 'N'
        lattice[end, 1:end] .= NBC;
    elseif BC[4] == 'P'
        lattice[end, 1:end] .= PBC;
    end

    return lattice

end;

function rectangularLatticeNeumannBoundaries(N::Int64, M::Int64, lattice::Matrix{Int64}, BCval::Vector{Float64})

    Nr = zeros(Float64, N, M);
    Nl = zeros(Float64, N, M);
    Nt = zeros(Float64, N, M);
    Nb = zeros(Float64, N, M);

    for i = 1:N
        Nl[i, 1] = (lattice[1+i-1, 1+1] == NBC) ? BCval[1] : 0.0;
        Nr[i, M] = (lattice[1+i+1, 1+M] == NBC) ? BCval[2] : 0.0;
    end
    for j = 1:M
        Nb[1, j] = (lattice[1+1, 1+j-1] == NBC) ? BCval[3] : 0.0;
        Nt[N, j] = (lattice[1+N, 1+j+1] == NBC) ? BCval[4] : 0.0;
    end

    return Nr, Nl, Nt, Nb;

end;

function rectangularLatticeDirichletBoundaries(N::Int64, M::Int64, lattice::Matrix{Int64}, BCval::Vector{Float64})

    Dr = zeros(Float64, N, M);
    Dl = zeros(Float64, N, M);
    Dt = zeros(Float64, N, M);
    Db = zeros(Float64, N, M);

    for i = 1:N
        Dl[i, 1] = (lattice[1+i, 1+1] == DBC) ? BCval[1] : 0.0;
        Dr[i, M] = (lattice[1+i, 1+M] == DBC) ? BCval[2] : 0.0;
    end
    for j = 1:M
        Db[1, j] = (lattice[1+1, 1+j] == DBC) ? BCval[3] : 0.0;
        Dt[N, j] = (lattice[1+N, 1+j] == DBC) ? BCval[4] : 0.0;
    end

    return Dr, Dl, Dt, Db;

end;

function initializeUniformSolution(N::Int64, M::Int64, ϕ0::Float64, lattice::Matrix{Int64}, Dr::Matrix{Float64}, Dl::Matrix{Float64}, Dt::Matrix{Float64}, Db::Matrix{Float64})

    ϕ = ones(Float64, N, M).*ϕ0;

    for i = 1:N
        ϕ[i, 1] = (lattice[1+i, 1+1] == DBC) ? Dl[i, 1] : ϕ[i, 1];
        ϕ[i, M] = (lattice[1+i, 1+M] == DBC) ? Dr[i, M] : ϕ[i, M];
    end

    for j = 1:M
        ϕ[1, j] = (lattice[1+1, 1+j] == DBC) ? Db[1, j] : ϕ[1, j];
        ϕ[N, j] = (lattice[1+N, 1+j] == DBC) ? Dt[N, j] : ϕ[N, j];
    end

    return ϕ;

end;

function getStencilComponents(N::Int64, M::Int64, a2x::Float64, a2y::Float64, ϵ::Matrix{Float64})

    stR = ϵ.*a2x;
    stR[2:end-1, :] .+= (ϵ[3:end, :].-ϵ[1:end-2, :]).*0.25.*a2x;
    stL = ϵ.*a2x;
    stL[2:end-1, :] .-= (ϵ[3:end, :].-ϵ[1:end-2, :]).*0.25.*a2x;
    stT = ϵ.*a2y;
    stT[:, 2:end-1] .+= (ϵ[:, 3:end].-ϵ[:, 1:end-2]).*0.25.*a2y;
    stB = ϵ.*a2y;
    stB[:, 2:end-1] .-= (ϵ[:, 3:end].-ϵ[:, 1:end-2]).*0.25.*a2y;
    stC = 2.0.*ϵ.*(a2x.+a2y);


    return stR, stL, stT, stB, stC;

end;

function getSource(N::Int64, M::Int64, dx::Float64, dy::Float64)

    ρ = zeros(Float64, N, M);

    for i = 1:N
        for j = 1:M
            ρ[i, j] = 0.0;#-8.0*π*π*sin(2.0*π*(dx*(i-1)+dy*(j-1)));
        end
    end

    return ρ

end;

function getLHS(N::Int64, M::Int64, lattice::Matrix{Int64}, r::Matrix{Float64}, l::Matrix{Float64}, t::Matrix{Float64}, b::Matrix{Float64}, c::Matrix{Float64})

    Nbands = 5;
    Dim = Nbands*M*N;
    A = zeros(Float64, Dim);
    I = zeros(Int64, Dim);
    J = zeros(Int64, Dim);
    
    index = 0;
    for i = 1:N
        for j = 1:M
            Ilat = i+1;
            Jlat = j+1;
            kc = j+(i-1)*M;
            if lattice[Ilat, Jlat] == REG
                L = ((lattice[Ilat-1, Jlat] == REG) ? l[i, j] : 0.0) + ((lattice[Ilat+1, Jlat] == NBC) ? r[i, j] : 0.0);
                Lp = (lattice[Ilat-1, Jlat] == PBC) ? l[i, j] : 0.0;
                if !isapprox(L, 0.0; atol=eps(Float64), rtol=0)
                    index += 1;
                    I[index] = kc-M;
                    J[index] = kc;
                    A[index] = L;
                elseif !isapprox(Lp, 0.0; atol=eps(Float64), rtol=0)
                    index += 1;
                    I[index] = j+(N-1)*M;
                    J[index] = kc;
                    A[index] = Lp;
                end

                B = ((lattice[Ilat, Jlat-1] == REG) ? b[i, j] : 0.0) + ((lattice[Ilat, Jlat+1] == NBC) ? t[i, j] : 0.0);
                Bp = (lattice[Ilat, Jlat-1] == PBC) ? b[i, j] : 0.0;
                if !isapprox(B, 0.0; atol=eps(Float64), rtol=0)
                    index += 1;
                    I[index] = kc-1;
                    J[index] = kc;
                    A[index] = B;
                elseif !isapprox(Bp, 0.0; atol=eps(Float64), rtol=0)
                    index += 1;
                    I[index] = M+(i-1)*M;
                    J[index] = kc;
                    A[index] = Bp;
                end

                C = -c[i, j];
                if !isapprox(C, 0.0; atol=eps(Float64), rtol=0)
                    index += 1;
                    I[index] = kc;
                    J[index] = kc;
                    A[index] = C;
                end

                T = ((lattice[Ilat, Jlat+1] == REG) ? t[i, j] : 0.0) + ((lattice[Ilat, Jlat-1] == NBC) ? b[i, j] : 0.0);
                Tp = (lattice[Ilat, Jlat+1] == PBC) ? t[i, j] : 0.0;
                if !isapprox(T, 0.0; atol=eps(Float64), rtol=0)
                    index += 1;
                    I[index] = kc+1;
                    J[index] = kc;
                    A[index] = T;
                elseif !isapprox(Tp, 0.0; atol=eps(Float64), rtol=0)
                    index += 1;
                    I[index] = 1+(i-1)*M;
                    J[index] = kc;
                    A[index] = Tp;
                end

                R = ((lattice[Ilat+1, Jlat] == REG) ? r[i, j] : 0.0) + ((lattice[Ilat-1, Jlat] == NBC) ? l[i, j] : 0.0);
                Rp = (lattice[Ilat+1, Jlat] == PBC) ? r[i, j] : 0.0;
                if !isapprox(R, 0.0; atol=eps(Float64), rtol=0)
                    index += 1;
                    I[index] = kc+M;
                    J[index] = kc;
                    A[index] = R;
                elseif !isapprox(Rp, 0.0; atol=eps(Float64), rtol=0)
                    index += 1;
                    I[index] = j+(1-1)*M;
                    J[index] = kc;
                    A[index] = Rp;
                end
            end
        end
    end

    return A, I, J, index;

end

function getRHS(N::Int64, M::Int64, dx::Float64, dy::Float64, lattice::Matrix{Int64}, ρ::Matrix{Float64}, r::Matrix{Float64}, l::Matrix{Float64}, t::Matrix{Float64}, b::Matrix{Float64}, Nr::Matrix{Float64}, Nl::Matrix{Float64}, Nt::Matrix{Float64}, Nb::Matrix{Float64}, Dr::Matrix{Float64}, Dl::Matrix{Float64}, Dt::Matrix{Float64}, Db::Matrix{Float64})

    B = zeros(Float64, N*M);
    K = zeros(Int64, N*M);

    index = 0;
    for i = 1:N
        for j = 1:M
            Ilat = i+1;
            Jlat = j+1;
            if lattice[Ilat, Jlat] == REG
                index += 1;
                B[index] = ρ[i, j];
                B[index] += 2.0*dx*(((lattice[Ilat-1, Jlat] == NBC) ? (l[i, j]*Nl[i, j]) : 0.0) - ((lattice[Ilat+1, Jlat] == NBC) ? (r[i, j]*Nr[i, j]) : 0.0));
                B[index] += 2.0*dy*(((lattice[Ilat, Jlat-1] == NBC) ? (b[i, j]*Nb[i, j]) : 0.0) - ((lattice[Ilat, Jlat+1] == NBC) ? (t[i, j]*Nt[i, j]) : 0.0));
                B[index] -= ((lattice[Ilat-1, Jlat] == DBC) ? (l[i, j]*Dl[i, j]) : 0.0)+((lattice[Ilat+1, Jlat] == DBC) ? (r[i, j]*Dr[i, j]) : 0.0);
                B[index] -= ((lattice[Ilat, Jlat-1] == DBC) ? (b[i, j]*Db[i, j]) : 0.0)+((lattice[Ilat, Jlat+1] == DBC) ? (t[i, j]*Dt[i, j]) : 0.0);
                K[index] = j+(i-1)*M;
            end
        end
    end

    return B, K, index;

end


Lx = 1.0;
Ly = 2.0;

N = 11;
M = 21;

dx = Lx/(N-1);
dy = Ly/(M-1);
d2x = dx*dx;
d2y = dy*dy;

ax = 1.0/dx;
ay = 1.0/dy;
a2x = ax*ax;
a2y = ay*ay;

BC = "NNNN";
BCval = [1.0, 1.0, 0.0, 0.0];

lattice = rectangularLattice(N, M, BC);

Nr, Nl, Nt, Nb = rectangularLatticeNeumannBoundaries(N, M, lattice, BCval);
Dr, Dl, Dt, Db = rectangularLatticeDirichletBoundaries(N, M, lattice, BCval);

ϵ = ones(Float64, N, M);
ρ = getSource(N, M, dx, dy);

ϕ0 = 0.0;
ϕ = initializeUniformSolution(N, M, ϕ0, lattice, Dr, Dl, Dt, Db);
ϕ1D = zeros(Float64, N*M);
for i = 1:N
    for j = 1:M
        ϕ1D[j+(i-1)*M] = ϕ[i, j];
    end
end

stR, stL, stT, stB, stC = getStencilComponents(N, M, a2x, a2y, ϵ);

A, I, J, nonzero_size1 = getLHS(N, M, lattice, stR, stL, stT, stB, stC);
b, K, nonzero_size2 = getRHS(N, M, dx, dy, lattice, ρ, stR, stL, stT, stB, Nr, Nl, Nt, Nb, Dr, Dl, Dt, Db);

As = sparse(I[1:nonzero_size1], J[1:nonzero_size1], A[1:nonzero_size1]);
bs = b[1:nonzero_size2];

x = IterativeSolvers.gmres!(ϕ1D, As, bs);
1 Like
2

I haven’t looked through your code, but with these kind of things I would suggest you get it working for one BC type (and test/validate on some simple problems), then add a second (and validate by testing on a problem where only the second is used), and then add the third (and validate by testing on a problem using only the third). This will help to ensure you are handling each BC correctly before you start trying to combine them.

3

I think I managed to solve the issue - as I was porting my code from C, some indexing was not done properly.
However, without going into the coding matters, what I would like to know is if there is a package or packages in Julia that could as an input take some arbitrary geometry (e.g., given as a text file) with arbitrarily distributed boundary conditions (external and internal) and then discretize those?
Or even simpler, is there a detailed tutorial which starts from scratch, shows a Poisson/Laplace 2D PDE, shows how to discretize it, to apply boundary conditions and then benchmarks different solvers?

I think for myself this would also be an interesting project to do, but as I am new to Julia, it would take me some time to get around things and to familiarize myself with the environment.

I really like the speed appeal of Julia, and as such it would be nice to have a fast-solver collection, as, I would say, lots of people in the scientific community would benefit from it.

In any case, I will continue slowly working on this to see if a satisfying solution can be achieved

4

Check out Gridap.

5

This seems very nice, but is this finite-elements-based, as it seems so?
At first instance I would be more interested in the finite differences, as for some simple geometries, they might be faster. But it’s also good to be aware of the FEM, as when the complexity of the geometry grows, FEM becomes inevitable.

6

FD may be fast, but the treatment of boundary conditions is always a pain in the neck. FE will be less work in the long run.

7

Possibly GitHub - ziolai/finite_element_electrical_engineering helps. Has both FDM and FEM.