Using stiff solvers for a large system of ode

General Usage Performance
, ,
1

Hi,

I am trying to solve the following set of partial differential equation:

\begin{cases} \partial_t u(x, y, t) = \partial_x f(x, y, t, u, \partial_x u) \\ \partial_t v(x, y, t) = \partial_y f(x, y, t, u, \partial_x u) \end{cases}

for some known function f. I want to use the a finite discretisation scheme of the shape

\begin{cases} \partial_t u(x, y, t) = D_x^f f(x, y, t, u, D_x^d u) \\ \partial_t v(x, y, t) = D_y^f f(x, y, t, u, D_x^d u) \end{cases}

where D^f stands for forward first order finite difference derivative and D^d for downward first order finite difference derivative, with the boundary condition that D^f=D^u on the right and D^d = 0 on the left.

Using the DifferentialEquations.jl package I arrived at the following working implementation:

using DifferentialEquations

function flux(i, j, t, f, p)
    T, μ, xl, yl = p

    m2π = f[i, j, 1]

    if i == 1
        m2σ = f[i, j, 1]
    else
        m2σ = f[i, j, 1] + xl[i] * (f[i-1, j, 1] - f[i, j, 1]) / (xl[i-1] / 2 - xl[i] / 2)
    end

    # details of the function `f`
    ϵπ = sqrt(t^2 + m2π)
    ϵσ = sqrt(t^2 + m2σ)
    ϵf = sqrt(t^2 + 3.25^2 * xl[i])
    Eplus = sqrt((ϵf + μ)^2 + 2 * 3.25^2 * yl[j])
    Eminus = sqrt((ϵf - μ)^2 + 2 * 3.25^2 * yl[j])

    return -t^4 / (3π^2 * ϵf) * (
        2(ϵf + μ) / Eplus * tanh(Eplus / 2T) +
        2(ϵf - μ) / Eminus * tanh(Eminus / 2T) +
        tanh((ϵf + μ) / 2T) +
        tanh((ϵf - μ) / 2T)
    ) + t^4 / 12π^2 * (3 / ϵπ * coth(ϵπ / 2T) + 1 / ϵσ * coth(ϵσ / 2T))
end


function rhs_pde!(dfdt, f, p, t)
    for I in CartesianIndices((length(p.xl), length(p.yl)))
        i, j = Tuple(I)

        # dudt
        if i == length(p.xl)
            dfdt[i, j, 1] =
                (flux(i - 1, j, t, f, p) - flux(i, j, t, f, p)) /
                (p.xl[i-1] / 2 - p.xl[i] / 2)
        else
            dfdt[i, j, 1] =
                (flux(i + 1, j, t, f, p) - flux(i, j, t, f, p)) /
                (p.xl[i+1] / 2 - p.xl[i] / 2)
        end

        # dvdt
        if j == length(p.yl)
            dfdt[i, j, 2] =
                (flux(i, j - 1, t, f, p) - flux(i, j, t, f, p)) /
                (p.yl[j-1] / 2 - p.yl[j] / 2)
        else
            dfdt[i, j, 2] =
                (flux(i, j + 1, t, f, p) - flux(i, j, t, f, p)) /
                (p.yl[j+1] / 2 - p.yl[j] / 2)
        end
    end
end

Nx = 80
Ny = 20

xl = range(0, 170, Nx) .^ 2
yl = range(0, 100, Ny) .^ 2

u_UV = [500359 + 1 / 2 * 21.0115 * x for x in xl, y in yl]
v_UV = [2 * 8.2e5 + 40 * y for x in xl, y in yl]

f0 = zeros(Nx, Ny, 2)
f0[:, :, 1] = u_UV
f0[:, :, 2] = v_UV

tspan = (1000, 200)

params = (T = 1, μ = 0, xl = xl, yl = yl)
prob = ODEProblem(rhs_pde!, f0, tspan, params)

sol = solve(
    prob,
    Tsit5(),
    save_everystep = false,
    reltol = 1e-12,
    abstol = 1e-12,
    dt = 1e-12,
    progress = true,
    progress_steps = 500,
)

I tried to adapt my code to a minimal working example in the context that it is used, hopefully the main structure is still understandable. This give the expected result in a reasonable amount of time using non-stiff algorithm such as DP5() or Tsit5().

My problem is that the equation get stiff for the value small values of t that I want to reach, as such I would like to use a stiff algorithm, but my code seems to be extremely slow with stiff solver. I checked that my function rhs_pde!() is non allocating but I am short on idea of how to further improve the performances. Any recommendation of what I could try, or I am doing something wrong?

2

I highly recommend reading Solving Large Stiff Equations · DifferentialEquations.jl

2 Likes
3

With great surprise the solution to my problem was on the documentation page carrying the name of my problem… I should have read it more carefully. Declaring a sparse Jacobian greatly improved the performances.

1 Like
4

One last thing worth mentioning is that the autoswitch methods may be valuable if your problem is only sometimes stiff.