[Bug]: `IDAKLUSolver` with `num_threads > 1` randomly crashes the process (double free / access violation) on Windows

[remoscream]

PyBaMM Version

26.3.1

Python Version

3.12.10

Describe the bug

Description

When using IDAKLUSolver with OpenMP multi-threading enabled (options={"num_threads": N}, N > 1), the Python process sometimes crashes silently — no Python exception, no traceback, the process just disappears. This happens randomly, not on every run, which makes it hard to reproduce on demand.

I confirmed the root cause using Windows crash dumps and WinDbg: it is a data race (thread-safety bug) inside the OpenMP parallel code of the solver.

Steps to reproduce

1. Create an IDAKLUSolver with options={"num_threads": N} where N > 1.
2. Run parameter identification / repeated simulations for an extended period (the crash is intermittent, so it may take many runs to appear).
3. The process is killed without any Python-level error. Windows Event Viewer shows an Application Error (Event ID 1000) pointing to idaklu.cp312-win_amd64.pyd.

Root cause

I analyzed two separate crash dump files with WinDbg (!analyze -v, then .ecxr + kb). Both crashes happened on OpenMP worker threads (frames from vcomp140.dll: vcomp_fork → ParallelRegion::HandlerThreadFunc), inside the idaklu module.

Crash 1 — double free:

Failure.Bucket: HEAP_CORRUPTION_ACTIONABLE_BlockNotBusy_DOUBLE_FREE_c0000374_ntdll.dll!RtlpFreeHeapInternal

A memory block is freed twice, from inside the OpenMP parallel region.

Crash 2 — invalid pointer read:

Failure.Bucket: INVALID_POINTER_READ_c0000005_idaklu.cp312-win_amd64.pyd!Unknown

Faulting instruction:

mov rcx, qword ptr [rax+rcx*8]   ; rcx holds a garbage value, not a valid index

Also on an OpenMP worker thread. The value used as an array index looks corrupted — consistent with a shared variable being read while another thread is still writing to it.

Both point to the same conclusion: some shared state inside the OpenMP parallel region in the IDAKLU solver is not properly synchronized between threads.

Scope

num_threads	Crash observed
1	No (never observed)
> 1	Yes (intermittent, not every run)

Expected behaviour

Running the solver with num_threads > 1 should not crash the process, regardless of thread scheduling timing.

Actual behaviour

The process is randomly killed with no Python-level error, only visible in the OS-level crash log (Windows Application Error) and confirmed as a double free / invalid memory read via WinDbg.

Environment

• OS: Windows 11
• Python: 3.12
• PyBaMM version: 26.3.1
• Used via: PyBOP (pybop-team/PyBOP), calling IDAKLUSolver with num_threads > 1

Additional context

We first started noticing this crash after we began identifying "function-type" parameters — i.e. adding an extra pybamm.InputParameter knob on top of an existing FunctionParameter (such as electrolyte_diffusivity_Nyman2008), via functools.partial-wrapped Python callables. Before that, we only identified plain scalar parameters and never observed this crash.

We are not certain whether this is a direct trigger or just a coincidence of timing — our best guess is that identifying function-type parameters adds extra sensitivity equations for IDAKLU to solve (since each new InputParameter needs a sensitivity w.r.t. it for the optimiser), which increases the parallel workload in the OpenMP region and may make the pre-existing race condition easier to hit. We are reporting this correlation in case it helps narrow down the affected code path, not as a confirmed root cause.

Steps to Reproduce

1. Build a PyBaMM model whose parameter set includes a FunctionParameter (e.g. electrolyte_diffusivity_Nyman2008(c_e, T)), and wrap it with an extra pybamm.InputParameter "knob" so the knob can be identified/optimised. Example:

   import pybamm
   
   def wrapped_diffusivity(original_function, knob_input, c_e, T):
       return knob_input * original_function(c_e, T)
   
   knob = pybamm.InputParameter("diffusivity_scale")
   parameter_values["Electrolyte diffusivity [m2.s-1]"] = functools.partial(
       wrapped_diffusivity, original_diffusivity_function, knob
   )

2. Create an IDAKLUSolver with multi-threading enabled:

   solver = pybamm.IDAKLUSolver(options={
       "num_threads": 8,   # any value > 1
       "num_solvers": 8,
   })

3. Run this model repeatedly inside an optimisation loop (e.g. via PyBOP), so the solver runs many times in a row (hundreds to thousands of solves), each solve computing sensitivities with respect to the knob input in parallel across threads.
4. Let it run for an extended period (in our case, anywhere from a few minutes to a few hours).
5. The process is killed silently — no Python exception, no traceback. Windows Event Viewer logs an Application Error (Event ID 1000) pointing to idaklu.cp312-win_amd64.pyd.

Notes:

• This is intermittent, not deterministic. It does not happen on every run, and we have never observed it with num_threads=1.
• We only started seeing this crash after we began identifying "function-type" parameters this way (step 1). Before that, we only identified plain scalar parameters and never saw the crash. We are not sure if this is a direct cause or just a coincidence — our guess is that the extra InputParameter knob adds more sensitivity equations for IDAKLU to compute, increasing the parallel workload and making the (likely pre-existing) race condition easier to trigger. See "Additional context" below.
• Because this is a data race (see "Root cause"), we could not build a small, always-reproducing minimal example. It only shows up in our real workload after repeated solves under load.

Relevant log output

Crash A — double free (python.exe.15484.dmp)

Failure.Bucket: HEAP_CORRUPTION_ACTIONABLE_BlockNotBusy_DOUBLE_FREE_c0000374_ntdll.dll!RtlpFreeHeapInternal
Failure.ProblemClass.Primary: HEAP_CORRUPTION
Failure.Exception.Code: 0xc0000374

EXCEPTION_RECORD: ntdll!RtlReportFatalFailure+0x9
NumberParameters: 1
Parameter[0]: 00007ff931ab48b0

Real call stack (.ecxr + kb):
00 ntdll!RtlReportFatalFailure+0x9
01 ntdll!RtlReportCriticalFailure+0x97
02 ntdll!RtlpHeapHandleError+0x12
03 ntdll!RtlpHpHeapHandleError+0x7a
04 ntdll!RtlpLogHeapFailure+0x45
05 ntdll!RtlpFreeHeapInternal+0x84c
06 ntdll!RtlFreeHeap+0x51
07 ucrtbase!_free_base+0x1b
08 idaklu_cp312_win_amd64!PyInit_idaklu+0xf86c8
09 idaklu_cp312_win_amd64!PyInit_idaklu+0x106e28
0a idaklu_cp312_win_amd64!PyInit_idaklu+0x4d018
0b idaklu_cp312_win_amd64!PyInit_idaklu+0x16f9e5
0c idaklu_cp312_win_amd64!PyInit_idaklu+0x110837
0d idaklu_cp312_win_amd64!PyInit_idaklu+0x4d107
0e idaklu_cp312_win_amd64!PyInit_idaklu+0x40f60
0f idaklu_cp312_win_amd64!PyInit_idaklu+0x1ee09
10 idaklu_cp312_win_amd64!PyInit_idaklu+0x2baa65
11 idaklu_cp312_win_amd64!PyInit_idaklu+0x2c882b
12 idaklu_cp312_win_amd64!PyInit_idaklu+0x2c7758
13 idaklu_cp312_win_amd64!PyInit_idaklu+0x2051a
14 idaklu_cp312_win_amd64!PyInit_idaklu+0x19c9d
15 vcomp140_31af29c03643f8396a6f26bcd601c636!_vcomp_fork_helper+0x6a
16 vcomp140_31af29c03643f8396a6f26bcd601c636!_vcomp::fork_helper_wrapper+0x9
17 vcomp140_31af29c03643f8396a6f26bcd601c636!_vcomp::ParallelRegion::HandlerThreadFunc+0x72
18 vcomp140_31af29c03643f8396a6f26bcd601c636!InvokeThreadTeam+0x115
19 vcomp140_31af29c03643f8396a6f26bcd601c636!vcomp_fork+0x1b1
1b idaklu_cp312_win_amd64!PyInit_idaklu+0x36e23
1c idaklu_cp312_win_amd64!PyInit_idaklu+0xc70b
1d idaklu_cp312_win_amd64!PyInit_idaklu+0x16ecb
1e idaklu_cp312_win_amd64!PyInit_idaklu+0x16f24
1f idaklu_cp312_win_amd64+0x21219
20 python312!PyEval_EvalFrameDefault+... (normal Python call chain continues)

Crash B — invalid pointer read (python.exe.9172.dmp)

Failure.Bucket: INVALID_POINTER_READ_c0000005_idaklu.cp312-win_amd64.pyd!Unknown
Failure.ProblemClass.Primary: INVALID_POINTER_READ
AV.Type: Read

EXCEPTION_RECORD: idaklu_cp312_win_amd64!PyInit_idaklu+0x16f920
ExceptionCode: c0000005 (Access violation)
Attempt to read from address 000014e12dfb4a30

Faulting instruction:
idaklu_cp312_win_amd64!PyInit_idaklu+0x16f920:
mov rcx, qword ptr [rax+rcx*8]

Real call stack (.ecxr + kb):
00 idaklu_cp312_win_amd64!PyInit_idaklu+0x16f920
01 idaklu_cp312_win_amd64!PyInit_idaklu+0x110837
02 idaklu_cp312_win_amd64!PyInit_idaklu+0x4d107
03 idaklu_cp312_win_amd64!PyInit_idaklu+0x4d03e
04 idaklu_cp312_win_amd64!PyInit_idaklu+0x16f9e5
05 idaklu_cp312_win_amd64!PyInit_idaklu+0x110837
06 idaklu_cp312_win_amd64!PyInit_idaklu+0x4d107
07 idaklu_cp312_win_amd64!PyInit_idaklu+0x40f60
08 idaklu_cp312_win_amd64!PyInit_idaklu+0x1ee09
09 idaklu_cp312_win_amd64!PyInit_idaklu+0x2baa65
0a idaklu_cp312_win_amd64!PyInit_idaklu+0x2c882b
0b idaklu_cp312_win_amd64!PyInit_idaklu+0x2c7758
0c idaklu_cp312_win_amd64!PyInit_idaklu+0x2051a
0d idaklu_cp312_win_amd64!PyInit_idaklu+0x19c9d
0e idaklu_cp312_win_amd64!PyInit_idaklu+0x2e854e
0f vcomp140_31af29c03643f8396a6f26bcd601c636!_vcomp_fork_helper+0x6a
10 vcomp140_31af29c03643f8396a6f26bcd601c636!_vcomp::fork_helper_wrapper+0x9
11 vcomp140_31af29c03643f8396a6f26bcd601c636!_vcomp::ParallelRegion::HandlerThreadFunc+0x72
12 vcomp140_31af29c03643f8396a6f26bcd601c636!_vcomp::PersistentThreadFunc+0x5f
13 kernel32!BaseThreadInitThunk+0x1d
14 ntdll!RtlUserThreadStart+0x28