I am trying to solve the following set of partial differential equation:
for some known function f. I want to use the a finite discretisation scheme of the shape
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.
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
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?