# Solving coupled PDEs with implicit time-dependence

**URL:** https://discourse.julialang.org/t/solving-coupled-pdes-with-implicit-time-dependence/115191
**Category:** Modelling & Simulations
**Tags:** physics, differentialequation, methodoflines
**Created:** [June 4, 2024, 5:36pm UTC](https://discourse.julialang.org/t/solving-coupled-pdes-with-implicit-time-dependence/115191 "2024-06-04T17:36:46Z")
**Posts on this page:** 2
**Page:** 1

<div class="post-metadata">

### Author: ![Danny\_Han](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danny_han/32/209679_2.png) [@Danny\_Han](https://discourse.julialang.org/u/Danny_Han)
#### Post date: [June 4, 2024, 5:36pm UTC](https://discourse.julialang.org/t/solving-coupled-pdes-with-implicit-time-dependence/115191/1 "2024-06-04T17:36:46Z")

</div>

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!

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

```

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [July 5, 2024, 6:26pm UTC](https://discourse.julialang.org/t/solving-coupled-pdes-with-implicit-time-dependence/115191/2 "2024-07-05T18:26:57Z")

</div>

Can you open an issue on MethodOfLines.jl? I just saw this and it seems to have gotten lost because it was never reported.