Alright, so I think we’re pretty much on the same page. Here are my definitions just to make sure:
x \in \mathbb{R}^n is basically the vector of all the free variables defining your system. That would be the pressure P, temperature T, the fraction of constituents y, \aleph, etc.
f : \mathbb{R}^n \to \mathbb{R} is the Gibbs energy you’re trying to minimize (I guess you denote it by G in your paper).
Ax = b simply regroups all of the linear equality constraints (e.g. \sum y = 1, etc).
Cx \geq d are all the linear inequality constraints in your problem (it could be for instance the non-negativity of the y’s, etc).
Putting aside the inequality constraints for now, my understanding is that you essentially use some form of Newton’s method to iteratively solve your problem. Given an initial guess x_0 (and assuming A x_0 = b, i.e. your initial guess satisfies your charge conservation constraint etc), you solve a sequence of problems of the form
\begin{bmatrix}
\mathcal{H}(x_t) & A^\top \\
A & 0
\end{bmatrix}
\begin{bmatrix}
\delta x \\
\delta v
\end{bmatrix}
=
\begin{bmatrix}
- \nabla f (x_t) \\
0
\end{bmatrix}
where subscript t denotes the iteration number, \delta x is the linear correction, and \delta v is related to the Lagrange multipliers required to enforce the desired constraints. Once you solved this linear system, you then set your new estimate as x_{t+1} = x_t + \delta x. Here, \mathcal{H}(x_t) is the Hessian (i.e. matrix of second derivatives) of the Gibbs energy evaluated at your current estimate of the solution.
I just want to make sure I get the mathematical problem right to see if there is anything (relatively easy) you could do from an optimization theoretic point of view to accelerate the solver since this appears to be the main complaint you here about the code.
Again, I am quite satisfied with OC but I can understand that someone who does not understand what is needed to find the equilibrium does not like to spend most of the CPU time solving a “subproblem” in his/her simulation.
Otherwise, I think you understand the algorithm OC uses. It was not invented by me, I think it goes back all the time to the initial thermodynamic calculations when working with ideal gases, for such a case the second derivatives of G are just RT/y.
When one has complex interactions as in solids or liquids this was replaced by algorithms using nonlinear minimization but 1981 Mats Hillert (Physica 103B (1981) pp 31-40) published an algorithm trying to apply this also to more complicated G functions and Bo Jansson implemented this in the Thermo-Calc software as part of his PhD thesis 1984. I simply implemented the same in OC which is free and managed to publish (in my view) a more simple derivation.
Inga problem, I understand that. Eventhough thermodynamics ain’t my forte, a handful of colleagues of mine are working on high-performance simulations of shock wave turbulence interactions for gases with non standard thermodynamics (e.g. dense gases, or systems with out-of-equilibrium chemistry and thermodynamics).
Since parallelization is not the answer to your problem, working on the minimization solver is the only way to speed-up the calculations (although it may lead you down the rabbit hole and that may not necessarily be something you want). I was simply trying to understand how you are currently solving the minimization problem to see if there is anything relatively easy you can do to speed-up your computations. With what you have and what I understand, there are a few things you could try which wouldn’t require a lot of development and could still lead to improved computational performances without sacrificing robustness.
I suppose that the computational bottleneck in your problem is not so much solving the linear system as much as computing all the n (n-1)/2 second derivatives needed to construct the Hessian matrix at every iteration. A simple trick would be to keep the Hessian matrix constant for a handful of iterations (even say just two). This is known as Picard iteration I believe. It may require a few more iterations to converge, but the overall cost of each iteration would be drastically smaller as you do not re-compute the second derivatives every time. That is probably the first thing I’d try if I were you.
Another possibility, since second-order methods (e.g. Newton) require computing the Hessian matrix which in itself is probably more costly than actually solving the resulting linear system, would be to use first order methods such as variants of gradient descent. From my understanding, the Gibbs energy has some form of coercivity (i.e. f(x) \to +\infty as \| x \| \to +\infty), so variants of gradient descent (such as the Heavy ball method for instance) should perform pretty well without careful tuning. The only thing required essentially is computing the gradient which you already do since it forms the right-hand side vector of your linear system.
Quasi-Newton methods such as BFGS should also perform remarkably well for a fraction of the cost of your current method eventhough it may require a bit more development time if you don’t have access to an existing library implementing it. It ain’t too complicated though and I teach to 2nd or 3rd year students in mechanical engineering.
I do have other ideas but that would quite require a bit more development and may be code restructuring which I guess you’re not interested in. Hope that helps anyway and give you some potential ideas
Thank you again. Skipping calculating the Hessian at every iteration is an interesting idea but as each phase can have different thermodynamic models it will not be trivial to implement. As I wrote using the Hessian allows more flexible set of conditions (specifying stable phases and phase constitution rather than just the overall composition) I am not sure if I can avoid it totally. BFSG I know nothing about, is there is free library for that (my software is free). The algorithm I use is anyway much faster than the normal Calphad algorithms which suffer from instabilities and cannot detect things such as miscibility gaps.
