Hi! I am currently working with the Cahn Hilliard equation in a fairly similar form to the demo provided in the docs. However, upon making a couple of modifications, the equations become quite stiff and the solver fails to converge.
Main
Currently, the 4th order equation is split into two coupled 2nd order equations and solved accordingly:
L0 = c*q*dx - c0*q*dx + dt*dot(grad(mu_mid), grad(q))*dx
L1 = mu*v*dx - dfdc*v*dx - lmbda*dot(grad(c), grad(v))*dx
L = L0 + L1
One strategy (even though more computationally expensive) I was thinking of trying was to decouple the equations and solve L0 first before inputting those results into L1 and iterating this in a loop to convergence. But I dont know how to implement this. Any help here would be much appreciated!
Extra stuff on numerics
I have also attempted to fiddle around with the choice of finite elements (DG vs CG vs CR of various orders). I am currently using the default LU solver:
problem = CahnHilliardEquation(a, L)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "residual"
solver.parameters["relative_tolerance"] = 1e-6
Might there be alternatives in these settings or the use of a preconditioner which might improve stability?
Thanks for your help!