Fairly new to Julia, and I am trying to solve a 2D PDE problem of the form:

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