MKL issue for Brusselator on Windows+AMD

[ArnoStrouwen]

using OrdinaryDiffEq, LinearAlgebra, SparseArrays

const N = 32
const xyd_brusselator = range(0, stop = 1, length = N)
brusselator_f(x, y, t) = (((x - 0.3)^2 + (y - 0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.0
limit(a, N) = a == N + 1 ? 1 : a == 0 ? N : a
function brusselator_2d_loop(du, u, p, t)
    A, B, alpha, dx = p
    alpha = alpha / dx^2
    @inbounds for I in CartesianIndices((N, N))
        i, j = Tuple(I)
        x, y = xyd_brusselator[I[1]], xyd_brusselator[I[2]]
        ip1, im1, jp1, jm1 = limit(i + 1, N), limit(i - 1, N), limit(j + 1, N),
        limit(j - 1, N)
        du[i, j, 1] = alpha * (u[im1, j, 1] + u[ip1, j, 1] + u[i, jp1, 1] + u[i, jm1, 1] -
                       4u[i, j, 1]) +
                      B + u[i, j, 1]^2 * u[i, j, 2] - (A + 1) * u[i, j, 1] +
                      brusselator_f(x, y, t)
        du[i, j, 2] = alpha * (u[im1, j, 2] + u[ip1, j, 2] + u[i, jp1, 2] + u[i, jm1, 2] -
                       4u[i, j, 2]) +
                      A * u[i, j, 1] - u[i, j, 1]^2 * u[i, j, 2]
    end
end
p = (3.4, 1.0, 10.0, step(xyd_brusselator))

function init_brusselator_2d(xyd)
    N = length(xyd)
    u = zeros(N, N, 2)
    for I in CartesianIndices((N, N))
        x = xyd[I[1]]
        y = xyd[I[2]]
        u[I, 1] = 22 * (y * (1 - y))^(3 / 2)
        u[I, 2] = 27 * (x * (1 - x))^(3 / 2)
    end
    u
end
u0 = init_brusselator_2d(xyd_brusselator)
prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop, u0, (0.0, 11.5), p)
solve(prob_ode_brusselator_2d, TRBDF2(), save_everystep = false)

(tester) pkg> status
Status `C:\Users\arno\Desktop\tester\Project.toml`
  [1dea7af3] OrdinaryDiffEq v6.89.0 `https://github.com/SciML/OrdinaryDiffEq.jl.git#master`

julia> versioninfo()
Julia Version 1.10.5
Commit 6f3fdf7b36 (2024-08-27 14:19 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 24 × AMD Ryzen 9 5900X 12-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, znver3)
Threads: 12 default, 0 interactive, 6 GC (on 24 virtual cores)
Environment:
  JULIA_NUM_THREADS = 12

Linux works, it is Windows specific.
Rober works on both operating systems, so it is something specific to the Brusselator example.
The automatic solver choice works, so it is something specific to BDF maybe?

[ArnoStrouwen]

@SebastianM-C could not reproduce, so I did some more digging:

while integrator.t < 11.5
    println("integrator.t\n  ", integrator.t)
    println("computer time")
    @time step!(integrator)
    println()
end

integrator.t
  0.0
computer time
  1.438802 seconds (1.44 M allocations: 96.120 MiB, 1.20% gc time, 93.44% compilation time)

integrator.t
  1.0099846788902428e-7
computer time
  0.026778 seconds (28 allocations: 49.375 KiB)

integrator.t
  1.1109831467792668e-6
computer time
  0.028883 seconds (28 allocations: 49.375 KiB)

integrator.t
  7.85468002318867e-6
computer time

Then nothing after that.

[SebastianM-C]

Digging more into this, it looks like the issue might be more specific than Windows, it seems that AMD processors are affected.

On

julia> versioninfo()
Julia Version 1.10.5
Commit 6f3fdf7b36 (2024-08-27 14:19 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 8 × 11th Gen Intel(R) Core(TM) i7-1185G7 @ 3.00GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, tigerlake)
Threads: 1 default, 0 interactive, 1 GC (on 8 virtual cores)

it works, but on

julia> versioninfo()
Julia Version 1.10.5
Commit 6f3fdf7b36 (2024-08-27 14:19 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 16 × AMD Ryzen 7 5800H with Radeon Graphics
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, znver3)
Threads: 1 default, 0 interactive, 1 GC (on 16 virtual cores)

I can reproduce the issue.

Additionally, when I tried to interrupt the solver I got:

julia> solve(prob_ode_brusselator_2d, TRBDF2(), save_everystep = false)
WARNING: Force throwing a SIGINT
ERROR: InterruptException:
Stacktrace:
  [1] getrf!(A::Matrix{Float64}; ipiv::Vector{Int64}, info::Base.RefValue{Int64}, check::Bool)
    @ LinearSolve C:\Users\sebastian\.julia\packages\LinearSolve\0Q5jw\src\mkl.jl:63
  [2] getrf!
    @ C:\Users\sebastian\.julia\packages\LinearSolve\0Q5jw\src\mkl.jl:51 [inlined]
  [3] solve!(cache::LinearSolve.LinearCache{…}, alg::LinearSolve.MKLLUFactorization; kwargs::@Kwargs{…})
    @ LinearSolve C:\Users\sebastian\.julia\packages\LinearSolve\0Q5jw\src\mkl.jl:219
  [4] solve!
    @ C:\Users\sebastian\.julia\packages\LinearSolve\0Q5jw\src\mkl.jl:213 [inlined]
  [5] macro expansion
    @ C:\Users\sebastian\.julia\packages\LinearSolve\0Q5jw\src\default.jl:363 [inlined]
  [6] solve!(::LinearSolve.LinearCache{…}, ::LinearSolve.DefaultLinearSolver; assump::LinearSolve.OperatorAssumptions{…}, kwargs::@Kwargs{…})
    @ LinearSolve C:\Users\sebastian\.julia\packages\LinearSolve\0Q5jw\src\default.jl:356
  [7] solve!
    @ C:\Users\sebastian\.julia\packages\LinearSolve\0Q5jw\src\default.jl:356 [inlined]
  [8] #solve!#8
    @ C:\Users\sebastian\.julia\packages\LinearSolve\0Q5jw\src\common.jl:274 [inlined]
  [9] solve!
    @ C:\Users\sebastian\.julia\packages\LinearSolve\0Q5jw\src\common.jl:273 [inlined]
 [10] #dolinsolve#2
    @ C:\Users\sebastian\.julia\packages\OrdinaryDiffEqDifferentiation\wkzwo\src\linsolve_utils.jl:29 [inlined]
 [11] dolinsolve
    @ C:\Users\sebastian\.julia\packages\OrdinaryDiffEqDifferentiation\wkzwo\src\linsolve_utils.jl:5 [inlined]
 [12] compute_step!(nlsolver::OrdinaryDiffEqNonlinearSolve.NLSolver{…}, integrator::OrdinaryDiffEqCore.ODEIntegrator{…},
 γW::Float64)
    @ OrdinaryDiffEqNonlinearSolve C:\Users\sebastian\.julia\packages\OrdinaryDiffEqNonlinearSolve\EXjDj\src\newton.jl:224
 [13] nlsolve!(nlsolver::OrdinaryDiffEqNonlinearSolve.NLSolver{…}, integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, cache::OrdinaryDiffEqSDIRK.TRBDF2Cache{…}, repeat_step::Bool)
    @ OrdinaryDiffEqNonlinearSolve C:\Users\sebastian\.julia\packages\OrdinaryDiffEqNonlinearSolve\EXjDj\src\nlsolve.jl:52
 [14] perform_step!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, cache::OrdinaryDiffEqSDIRK.TRBDF2Cache{…}, repeat_step::Bool)
    @ OrdinaryDiffEqSDIRK C:\Users\sebastian\.julia\packages\OrdinaryDiffEqSDIRK\9xDev\src\sdirk_perform_step.jl:474
 [15] perform_step!
    @ C:\Users\sebastian\.julia\packages\OrdinaryDiffEqSDIRK\9xDev\src\sdirk_perform_step.jl:451 [inlined]
 [16] solve!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…})
    @ OrdinaryDiffEqCore C:\Users\sebastian\.julia\packages\OrdinaryDiffEqCore\55UVY\src\solve.jl:551
 [17] __solve(::ODEProblem{…}, ::TRBDF2{…}; kwargs::@Kwargs{…})
    @ OrdinaryDiffEqCore C:\Users\sebastian\.julia\packages\OrdinaryDiffEqCore\55UVY\src\solve.jl:7
 [18] __solve
    @ C:\Users\sebastian\.julia\packages\OrdinaryDiffEqCore\55UVY\src\solve.jl:1 [inlined]
 [19] #solve_call#44
    @ C:\Users\sebastian\.julia\packages\DiffEqBase\uqSeD\src\solve.jl:612 [inlined]
 [20] solve_call
    @ C:\Users\sebastian\.julia\packages\DiffEqBase\uqSeD\src\solve.jl:569 [inlined]
 [21] #solve_up#53
    @ C:\Users\sebastian\.julia\packages\DiffEqBase\uqSeD\src\solve.jl:1092 [inlined]
 [22] solve_up
    @ C:\Users\sebastian\.julia\packages\DiffEqBase\uqSeD\src\solve.jl:1078 [inlined]
 [23] #solve#51
    @ C:\Users\sebastian\.julia\packages\DiffEqBase\uqSeD\src\solve.jl:1015 [inlined]
 [24] top-level scope
    @ REPL[13]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> versioninfo()

Please submit a bug report with steps to reproduce this fault, and any error messages that follow (in their entirety). Thanks.
Exception:
Please submit a bug repor at 0x7ffa35961f9e --  at 0x7ffa35961f9e -- port with steps to reproduce this fault, and any error messages that follow (in their entirety). Thanks.
Exception: EXCEPTION_ACCESS_VIOLATION at 0x7ffa35961f9e -- IONrs\sebastian\.julia\artifacts\07d4e2a4c926ce5a99f7ff402617dc9caa2187a0\bin\mkl_intel_thread.2.dll (unknown line)EXCEPTION_ACCESS_VIOLATION at 0x7ffa35961f9e --  at 0x7ffa35961f9e -- mkl_lapack_dgetrf at C:\Users\sebastian\.julia\artifacts\07d4e2a4c926ce5a99f7ff402617dc9caa2187a0\bin\mkl_intel_thread.2.dll (unknown line)
in expression starting at REPL[14]:1
Allocations: 35014975 (Pool: 34986227; Big: 28748); GC: 38
Allocations: 35014975 (Pool: 34986227; Big: 28748); GC: 38
4e2a4c926ce5a99f7ff402617dc9caa2187a0\bin\mkl_intel_thread.2.dll (unknown line)
in expression starting at REPL[14]:1
Allocations: 35014975 (Pool: 34986227; Big: 28748); GC: 38
Allocations: 35014975 (Pool: 34986227; Big: 28748); GC: 38
4e2a4c926ce5a99f7ff402617dc9caa2187a0\bin\mkl_intel_thread.2.dll (unknown line)

Also, OrdinaryDiffEq is not needed, it happens on the latest release of OrdinaryDiffEqSDIRK too.

[ChrisRackauckas]

That looks like an MKL issue on AMD?

[ChrisRackauckas]

@staticfloat @giordano do you know if the latest MKL has a 64-bit version that is being vendored correctly?

[giordano]

I'm not aware of particular issues with MKL, we just ship whatever binary they distribute.