Efficient solution strategies for 3D Cahn-Hilliard-Navier-Stokes linear systems (Newton iteration) with slow MUMPS direct solver in FEniCS

I am solving the 3D Cahn-Hilliard-Navier-Stokes (CHNS) coupled equations in FEniCS with a finite element discretization scheme as follows:

  • A mixed finite element method is applied to discretize the Navier-Stokes equations;
  • Standard P1 Lagrange finite elements are used for the remaining Cahn-Hilliard equations.

This spatial discretization yields a large-scale nonlinear discrete system, which I solve using the Newton iterative method. At each Newton update step, I have to solve the corresponding Jacobian linear system. Currently, I am using the MUMPS direct solver for this linear solve, but the computational speed is extremely slow for the 3D case, and the overall computational complexity is prohibitively high, which has become a bottleneck for my entire simulation.
Any relevant literature references, FEniCS code snippets, or personal experience with this type of problem would be of great help to me. I would be very grateful for any assistance or insights you can provide.

Without having the actual code, and a breakdown of timings, it’s hard to give guidance.
Some of the things that might speed up your code is covered in: