I’m trying to implement a time-dependent Landau model for ferroelectric materials on a 2D cartesian domain, where the time-dependent Landau equation
\frac{dP}{dt} = -\eta (\alpha P+\beta P^3 + \gamma P^5 - k_x \frac{\partial^2 P}{\partial x^2} - k_z \frac{\partial^2 P}{\partial z^2} + P\frac{dV}{dz}) is coupled with the Poisson equation
\nabla^2 V = \frac{1}{\epsilon} \frac{\partial P}{\partial z}
In this system, V is not explicitly dependent on time, and the typical solution routine solves for V for a given P, then uses that solution for V to time step P, and repeats.
I’m wondering which Julia package would be most suited for this task. My first attempt was with MethodOfLines.jl, but the following code froze upon execution (>5min on the discretize step, no output). Any help is appreciated!
mutable struct Ferroelectric
# Parameters
ϵ::Float64 # Permittivity
# Landau coefficients
α::Float64
β::Float64
γ::Float64
# Domain wall coupling coefficients
k_x::Float64
k_z::Float64
# dynamic coefficient
η::Float64
end
using MethodOfLines, ModelingToolkit, DomainSets, Plots, OrdinaryDiffEq
@parameters t x z
@variables P(..) V(..)
Dt= Differential(t); Dz = Differential(z); Dx = Differential(x)
Dxx = Differential(x)^2; Dzz = Differential(z)^2; Dtt = Differential(t)^2
function simulate(fe::Ferroelectric, lims::Tuple{Float64, Float64, Float64}, ΔV::Float64, Δx::Float64, Δz::Float64, Δt::Float64, order::Int)
x_lim, z_lim, t_lim = lims
# Define the PDE System
eqs = [Dt(P(x, z, t)) ~ -fe.η * (fe.α * P(x, z, t) + fe.β * P(x, z, t)^3 + fe.γ * P(x, z, t)^5 - fe.k_x * Dxx(P(x, z, t)) - fe.k_z * Dzz(P(x, z, t)) + P(x, z, t) * Dz(V(x, z, t))),
Dxx(V(x, z, t)) + Dzz(V(x, z, t)) ~ P(x, z, t) / fe.ϵ]
print("System defined \n")
# Define the domain
domain = [x ∈ Interval(0.0, x_lim),
z ∈ Interval(0.0, z_lim),
t ∈ Interval(0.0, t_lim)]
# boundary conditions
bcs = [P(0, z, t) ~ P(x_lim, z, t), # PBC for P
P(x, z, 0) ~ 5, # initial polarization
V(0, z, t) ~ V(x_lim, z, t), # PBC for V
V(x, 0, t) ~ 0,
V(x, z_lim, t) ~ ΔV, # voltage applied at the top
V(x, z, 0) ~ ΔV * (z / z_lim), # initial potential distribution
]
@named pdesys = PDESystem(eqs, bcs, domain, [x, z, t], [P(x, z, t), V(x, z, t)])
# Discretize the PDE System
discretization = MOLFiniteDifference([x => Δx, z => Δz], t, approx_order=order)
prob = discretize(pdesys, discretization)
print("System discretized \n")
# Solve the PDE System
sol = solve(prob, Tsit5(), saveat=Δt)
return sol
end
test_fe = Ferroelectric(16., 1., 1., 1., 0.1, 0.1, 0.1) #arbitrary test values
simulate(test_fe, (1.0, 1.0, 1.0), 1.0, 0.1, 0.1, 0.1, 2)