Laplace's Equation with Internal Dirichlet Conditions

In this tutorial, we consider Laplace's equation with some additional complexity put into the problem via internal Dirichlet conditions:

\[\begin{equation} \begin{aligned} \grad^2 u &= 0 & \vb x \in [0, 1]^2, \\ u(0, y) &= 100y & 0 \leq y \leq 1, \\ u(1, y) &= 100y & 0 \leq y \leq 1, \\ u(x, 0) &= 0 & 0 \leq x \leq 1, \\ u(x, 1) &= 100 & 0 \leq x \leq 1, \\ u(1/2, y) &= 0 & 0 \leq y \leq 2/5. \end{aligned} \end{equation}\]

To start with solving this problem, let us define an initial mesh.

using DelaunayTriangulation, FiniteVolumeMethod
tri = triangulate_rectangle(0, 1, 0, 1, 50, 50, single_boundary=false)
Delaunay Triangulation.
    Constrained: true
    Has ghost triangles: true
    Number of points: 2500
    Number of triangles: 4998
    Number of edges: 7501

In this mesh, we don't have any points that lie exactly on the line $\{x = 1/2, 0 \leq y \leq 2/5\}$, so we cannot enforce this constraint exactly.[1] Instead, we need to add these points into tri. We do not need to add any constrained edges in this case, since these internal conditions are enforced only at points.

Let us now add in the points.

using CairoMakie
new_points = LinRange(0, 2 / 5, 250)
for y in new_points
    add_point!(tri, 1 / 2, y)
end
fig, ax, sc = triplot(tri)
fig

It may also help to refine the mesh slightly.

refine!(tri, max_area=1e-4)
fig, ax, sc = triplot(tri)
fig

mesh = FVMGeometry(tri)
FVMGeometry with 10223 control volumes, 20042 triangles, and 30264 edges

Now that we have the mesh, we can define the boundary conditions. Remember that the order of the boundary indices is the bottom wall, right wall, top wall, and then the left wall.

bc_bot = (x, y, t, u, p) -> zero(u)
bc_right = (x, y, t, u, p) -> oftype(u, 100y) # helpful to have each bc return the same type
bc_top = (x, y, t, u, p) -> oftype(u, 100)
bc_left = (x, y, t, u, p) -> oftype(u, 100y)
bcs = (bc_bot, bc_right, bc_top, bc_left)
types = (Dirichlet, Dirichlet, Dirichlet, Dirichlet)
BCs = BoundaryConditions(mesh, bcs, types)
BoundaryConditions with 4 boundary conditions with types (Dirichlet, Dirichlet, Dirichlet, Dirichlet)

We now need to define the internal conditions. This is done using InternalConditions. First, we need to find all the vertices that lie on the line $\{x = 1/2, 0 \leq y \leq 2/5\}$. We could compute these manually, but let's find them programmatically instead for the sake of demonstration.

function find_all_points_on_line(tri)
    vertices = Int[]
    for i in each_solid_vertex(tri)
        x, y = get_point(tri, i)
        if x == 1 / 2 && 0 ≤ y ≤ 2 / 5
            push!(vertices, i)
        end
    end
    return vertices
end
vertices = find_all_points_on_line(tri)
fig, ax, sc = triplot(tri)
points = [get_point(tri, i) for i in vertices]
scatter!(ax, points, color=:red, markersize=10)
fig

Now that we have the vertices, we can define the internal conditions. We need to provide InternalConditions with a Dict that maps each vertex in vertices to a function index that corresponds to the condition for that vertex. In this case, that function index is 1 as we only have a single function.

ICs = InternalConditions((x, y, t, u, p) -> zero(u),
    dirichlet_nodes=Dict(vertices .=> 1))
InternalConditions with 250 Dirichlet nodes and 0 Dudt nodes

Now we can define the problem. As discussed in the Helmholtz tutorial, we are looking to define a steady state problem, and so the initial condition needs to be a suitable initial guess of what the solution could be. Looking to the boundary and internal conditions, one suitable guess is $u(x, y) = 100y$ with $u(1/2, y) = 0$ for $0 \leq y \leq 2/5$; in fact, $u(x, y) = 100y$ is the solution of the problem without the internal condition. Let us now use this to define our initial condition.

initial_condition = zeros(DelaunayTriangulation.num_solid_vertices(tri))
for i in each_solid_vertex(tri)
    x, y = get_point(tri, i)
    initial_condition[i] = ifelse(x == 1 / 2 && 0 ≤ y ≤ 2 / 5, 0, 100y)
end

Now let's define the problem. The internal conditions are provided as the third argument of FVMProblem.

diffusion_function = (x, y, t, u, p) -> one(u) # ∇²u = ∇⋅[D∇u], D = 1
final_time = Inf
prob = FVMProblem(mesh, BCs, ICs;
    diffusion_function,
    initial_condition,
    final_time)
FVMProblem with 10223 nodes and time span (0.0, Inf)
steady_prob = SteadyFVMProblem(prob)
SteadyFVMProblem with 10223 nodes

Now let's solve the problem.

using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq
sol = solve(steady_prob, DynamicSS(TRBDF2(linsolve=KLUFactorization())))
u: 10223-element Vector{Float64}:
  0.0
  0.0
  0.0
  0.0
  0.0
  ⋮
 52.99056010002759
 60.20408163265305
 82.65306122448979
  0.0
 64.28571428571428
fig, ax, sc = tricontourf(tri, sol.u, levels=LinRange(0, 100, 28))
tightlimits!(ax)
fig

Just the code

An uncommented version of this example is given below. You can view the source code for this file here.

using DelaunayTriangulation, FiniteVolumeMethod
tri = triangulate_rectangle(0, 1, 0, 1, 50, 50, single_boundary=false)

using CairoMakie
new_points = LinRange(0, 2 / 5, 250)
for y in new_points
    add_point!(tri, 1 / 2, y)
end
fig, ax, sc = triplot(tri)
fig

refine!(tri, max_area=1e-4)
fig, ax, sc = triplot(tri)
fig

mesh = FVMGeometry(tri)

bc_bot = (x, y, t, u, p) -> zero(u)
bc_right = (x, y, t, u, p) -> oftype(u, 100y) # helpful to have each bc return the same type
bc_top = (x, y, t, u, p) -> oftype(u, 100)
bc_left = (x, y, t, u, p) -> oftype(u, 100y)
bcs = (bc_bot, bc_right, bc_top, bc_left)
types = (Dirichlet, Dirichlet, Dirichlet, Dirichlet)
BCs = BoundaryConditions(mesh, bcs, types)

function find_all_points_on_line(tri)
    vertices = Int[]
    for i in each_solid_vertex(tri)
        x, y = get_point(tri, i)
        if x == 1 / 2 && 0 ≤ y ≤ 2 / 5
            push!(vertices, i)
        end
    end
    return vertices
end
vertices = find_all_points_on_line(tri)
fig, ax, sc = triplot(tri)
points = [get_point(tri, i) for i in vertices]
scatter!(ax, points, color=:red, markersize=10)
fig

ICs = InternalConditions((x, y, t, u, p) -> zero(u),
    dirichlet_nodes=Dict(vertices .=> 1))

initial_condition = zeros(DelaunayTriangulation.num_solid_vertices(tri))
for i in each_solid_vertex(tri)
    x, y = get_point(tri, i)
    initial_condition[i] = ifelse(x == 1 / 2 && 0 ≤ y ≤ 2 / 5, 0, 100y)
end

diffusion_function = (x, y, t, u, p) -> one(u) # ∇²u = ∇⋅[D∇u], D = 1
final_time = Inf
prob = FVMProblem(mesh, BCs, ICs;
    diffusion_function,
    initial_condition,
    final_time)

steady_prob = SteadyFVMProblem(prob)

using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq
sol = solve(steady_prob, DynamicSS(TRBDF2(linsolve=KLUFactorization())))

fig, ax, sc = tricontourf(tri, sol.u, levels=LinRange(0, 100, 28))
tightlimits!(ax)
fig

This page was generated using Literate.jl.

  • 1Of course, by defining the grid spacing appropriately we could have such points, but we just want to show here how we can add these points in if needed.