Steady State Solution of Drift Diffusion Equation

[apekshasingh]

Hi there, thanks for your work on these packages. I had originally submitted an issue on the DiffEqOperators package and was suggested to try out MethodOfLines instead (https://github.com/SciML/DiffEqOperators.jl/issues/518). I'm trying to implement the same equation as described in my previous issue and solve for steady state solutions. Following the steady state heat equation example under the tutorial, I implemented the following. As you can see, my PDE equation relies on additional defined functions and I'm not sure if I'm implementing them appropriately. Any guidance would be greatly appreciated!

@parameters x y
@variables c(..)
Dx = Differential(x)
Dy = Differential(y)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

∇²(c) = Dxx(c) + Dyy(c)

D = 0.01

function V(x, y)
    [x^2-x*y, y^2-x*y]
end

function R(x, y)
    r = sqrt(x^2 + y^2)
    (r^4/(0.1 + r^4)) + 0.01
end

function S(x, y)
    if x==y==0
        1
    else
        0
    end
end

@register V(x, y)
@register R(x, y)
@register S(x, y)

x_min = y_min = 0.0
x_max = y_max = 1.0

eq = 0.5*D*∇²(c(x, y)) - Dx(c(x,y)*V(x,y)) - Dy(c(x,y)*V(x,y)) + R(x, y)*c(x, y) + S(x, y) ~ 0

bcs = [c(x_max, y) ~ Dx(c(x_max, y)),
c(x, y_max) ~ Dy(c(x, y_max)),
Dx(c(x_min, y)) ~ 0,
Dy(c(x, y_min)) ~ 0]

domains = [x ∈ Interval(0.0, 1.0),
           y ∈ Interval(0.0, 1.0)]

@named pdesys = PDESystem([eq], bcs, domains, [x, y], [c(x, y)])

dx = 0.1
dy = 0.1

discretization = MOLFiniteDifference([x => dx, y => dy], nothing, approx_order=2)

prob = discretize(pdesys, discretization)
sol = NonlinearSolve.solve(prob, NewtonRaphson())

grid = get_discrete(pdesys, discretization)

c_sol = map(d -> sol[d], grid[c(x, y)])

using Plots

heatmap(grid[x], grid[y], c_sol, xlabel="x values", ylabel="y values",
        title="Steady State Solution") 

[apekshasingh]

Rereading the documentation, I see that I cannot have terms like Dx(c(x,y)*V(x,y)). I have rewritten it in the following way:

@parameters x y
@variables c(..)
Dx = Differential(x)
Dy = Differential(y)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

∇²(c) = Dxx(c) + Dyy(c)

D = 0.01

function Vx(x, y)
    x^2-x*y
end

function Vy(x, y)
    y^2-x*y
end

function dVx(x, y)
    2x - y
end

function dVy(x, y)
    2y - x
end

function R(x, y)
    r = sqrt(x^2 + y^2)
    (r^4/(0.1 + r^4)) + 0.01
end

function S(x, y)
    if x==y==0
        1
    else
        0
    end
end

@register Vx(x, y)
@register Vy(x, y)
@register dVx(x, y)
@register dVy(x, y)
@register R(x, y)
@register S(x, y)

x_min = y_min = 0.0
x_max = y_max = 1.0

eq = 0.5*D*∇²(c(x, y)) - Dx(c(x,y))*Vx(x,y)- c(x,y)*dVx(x,y) - Dy(c(x,y))*Vy(x,y)- c(x,y)*dVy(x,y) + R(x, y)*c(x, y) + S(x, y) ~ 0

bcs = [ Dx(c(x_min, y)) ~ 0,
Dy(c(x, y_min)) ~ 0, 
c(x_max, y) ~ 0,
c(x, y_max) ~ 0]

domains = [x ∈ Interval(0.0, 1.0),
           y ∈ Interval(0.0, 1.0)]

@named pdesys = PDESystem([eq], bcs, domains, [x, y], [c(x, y)])

dx = 0.1
dy = 0.1

discretization = MOLFiniteDifference([x => dx, y => dy], nothing, approx_order=2)

prob = discretize(pdesys, discretization)
sol = NonlinearSolve.solve(prob, NewtonRaphson())

grid = get_discrete(pdesys, discretization)

c_sol = map(d -> sol[d], grid[c(x, y)])

using Plots

heatmap(grid[x], grid[y], c_sol, xlabel="x values", ylabel="y values",
        title="Steady State Solution")

This does generate results, however they do not look like how I would expect. I believe everything in my model is symmetric in space (x, y), however when I plot solutions they are not symmetric about the origin. I'm not sure if it is an issue with the way I have specified the boundary conditions (at the boundary x=0 or y=0, I will most likely have reflective boundaries; at x=1 or y=1 I will likely have absorptive boundaries; it's a bit unclear to me how corner points are treated) or with one of the terms in the equation.

[xtalax]

Have you updated to v0.3.1? This fixed an issue with lhs advection terms like those appearing in your equation.

[apekshasingh]

Thanks for your quick reply. I just started with the package yesterday for the first time and the version is listed as MethodOfLines v0.3.1. Simplifying the equation to just ∇²(c(x, y)) - Dx(c(x,y))- Dy(c(x,y)) +S(x,y) ~ 0 still results in solutions not symmetric about x=y as might be expected.

[xtalax]

dVx and dVy and S(x, y) don't look symmetrical, are you sure the result should be?

[apekshasingh]

That was my expectation, but I could be wrong. Working with the simpler ∇²(c(x, y)) - Dx(c(x,y))- Dy(c(x,y)) +S(x,y) ~ 0 and the bcs = [c(x_min, y) ~ 1, c(x, y_min) ~ 1, c(x_max, y) ~ 0, c(x, y_max) ~ 0] seem to give an ok result. Do you have any examples of drift-diffusion steady state solutions with non-Dirichlet bcs? Also could you clarify how boundary conditions at the corner points of the space are treated? Thanks so much for your help.

[xtalax]

I don't actually have any drift diffusion examples, I'd like to get this working so it can be added. Corner points are currently dirichlet 0, but they are not used in the computation.

[xtalax]

@apekshasingh Were you satisfied that this was working in the end?

[apekshasingh]

Thank you for following up. I actually did not end up pursuing this further at this time, but greatly appreciate your help.

On Fri, Oct 7, 2022 at 11:23 AM Alex Jones @.***> wrote:

@apekshasingh https://github.com/apekshasingh Were you satisfied that this was working in the end?

— Reply to this email directly, view it on GitHub https://github.com/SciML/MethodOfLines.jl/issues/111#issuecomment-1271919737, or unsubscribe https://github.com/notifications/unsubscribe-auth/AL5QHLMM6IT3N4EUFU5QS43WCBTA3ANCNFSM5YAUIAWA . You are receiving this because you were mentioned.Message ID: @.***>