I’m not aware of any changes in MUMPS that would cause this. But @amontoison is the expert, so he might have some ideas.
I suspect that the culprit is the new version of OpenBLAS32_jll.jl that can install with Julia 1.10.
New versions of OpenBLAS were recently compiled with different options (like the version of GCC).
You should try an older one with Julia 1.10.
Pkg> add OpenBLAS32_jll#0.3.21
Contrary to OpenBLAS_jll.jl, we can install different versions of OpenBLAS32_jll.jl for a specific version of Julia.
Thanks @tmigot for pointing towards OptimalControl.jl / controltoolbox
@A_C Interesting to see converging efforts for numerical optimal control in Julia About the points you raised:
 not sure about what you mean by multiple phases, but problems with several dynamics can easily be patched using the current framework
 mesh refinement: not a top priority but, yes, there have been collocation schemes with mesh refinement available for a while, would be nice to include this kind of solver, e.g. for BVPs
 not sure either about what you mean by scaling, but it is true that the direct part currently based on ADNLModels.jl  that very nicely fits our needs  needs performance improvement for large problems.
PS. @odow Great job with function tracing: wasn’t aware of that in JuMP, this is a welcome improvement to allow pure Julia code to be used
The problem with MUMPS is fixed. We should not compile the sequential version (without MPI) with OpenMP.
The release 1.6.1
of Ipopt.jl
provides the new release of Ipopt (v3.14.14) with upgraded linear solvers (MUMPS / SPRAL).
For large problems, the main bottleneck is the computation of the Hessian of the Lagrangian (if we use Ipopt).
The code for doing is very efficient in MathOptInterface.jl
without subexpressions.
One possibility is to create a MathOptNLPModel
with NLPModelsJuMP.jl for a part of the model and an ADNLPModel
for the other part that is more problematic with JuMP.
I have been referring to the book by Betts for Numerical Optimal Control.
 I think it is the same thing. A new phase starts at some event and then continues until some other event. For example, if one wants to first reach point x_1 and later the point x_2. Then first phase would be reaching point x_1. Travelling from x_1 to x_2 would be second phase. And these are linked in this example (though they need not always be) since initial state for second phase is the final state of first phase.
 One advantage of mesh refinement is that it allows you to start with small number of grid points. So it will be easier for the solver to converge and provide a solution which acts as initial guess for next refined grid which contains more points. Since the later grids have initial solution which is feasible it will be easier for solver to converge.
 Scaling is required to increase the robustness. For example, some input might have magnitude in
[0,1]
, while some other can be in between[0,10000]
. One can bring them to same range manually while formulating the problem. But this can also be done automatically.

OK. Multiphase problems can easily be recast as higher dimensional single phase one. Automating this procedure can be interesting within Julia that offers a powerful combination of symbolic / numeric to do so

Sure

Sure