Problem with the solutions of a Cahn-Hilliard type problem.

Question

Hi,

I am a little desperate because it is quite a long time I work on a project and my code not gives me the good solutions, sometimes the solver diverges. Here are the equations :

$$ \frac{d~u_1}{dt} = D_n \cdot \Delta u_1 - k_1 \cdot \nabla (\frac{u_1}{1-u_3}\cdot \nabla u_3) - \rho_0 \nabla (u_1 \nabla u_2)$$
$$ \frac{d~u_2}{dt} = \omega u_1 - \mu u_1 u_2$$
$$ \frac{d~u_3}{dt} = - \lambda u_1 u_3$$

And here is my code :

#coding=utf8

##############################

from fenics import *
from numpy import ones
import random

##############################

mesh = RectangleMesh(Point(0.0, 0.0), Point(15.0, 10.0), 90, 60)
        # final time
num_steps = 50    # number of time steps
dt = 5.0e-03 # time step size

# Parameters

Dn       = 0.001
chi_0    = 1.
k_1      = 0.1  
rho_0    = 0.2
omega    = 0.1
mu       = 0.5
lambda_0 = 10.


# Class representing the intial conditions
class InitialConditions(Expression):
    def __init__(self, **kwargs):
        random.seed(2 + MPI.rank(mpi_comm_world()))
    def eval(self, values, x):
        values[0] = exp(-(pow(x[0],2)+pow(x[1]-5,2))/0.1)+0.02*(0.5 - random.random())
        values[1] = (-atan(x[0]/3-2.5)/2.6)+1
        values[2] = (atan(x[0]/3-2.5)/2.6)+1
    def value_shape(self):
        return (3,)

# Class for interfacing with the Newton solver
class CellEquation(NonlinearProblem):
    def __init__(self, a, L):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
    def F(self, b, x):
        assemble(self.L, tensor=b)
    def J(self, A, x):
        assemble(self.a, tensor=A)

# Form compiler options
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "quadrature"

# Define function space for system of concentrations

P1 = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
element = MixedElement([P1, P1, P1])
V = FunctionSpace(mesh, element)
P = FunctionSpace(mesh, P1)

# Define test functions

v_1, v_2, v_3 = TestFunctions(V)

# Define functions for pressure and concentrations
du = TrialFunction(V)
#w = Function(W)
u = Function(V)
u_n = Function(V)   

# Define initial condition
u_init = InitialConditions(degree=1)
u.interpolate(u_init)
u_n.interpolate(u_init)

# Split systeme functions to access components
u_1, u_2, u_3 = split(u)
u_n1, u_n2, u_n3 = split(u_n)


# Define expressions used in variational forms

k = dt
Dn = Constant(Dn)
chi_0 = Constant(chi_0)
k_1 = Constant(k_1)
rho_0 = Constant(rho_0)
omega = Constant(omega)
mu = Constant(mu)
lambda_0 = Constant(lambda_0)

# Define variational problem

F = (u_1 - u_n1)*v_1*dx \
  - k*Dn*dot(grad(u_1), grad(v_1))*dx \
  + k*k_1*(u_1 / (1 + u_3))*dot(grad(u_3), grad(v_1))*dx \
  + k*rho_0*u_1*dot(grad(u_2), grad(v_1))*dx \
  + (u_2 - u_n2)*v_2*dx \
  - k*omega*u_1*v_2*dx + k*mu*u_1*u_2*v_2*dx \
  + (u_3 - u_n3)*v_3*dx + k*lambda_0*u_1*u_3*v_3*dx

# Create VTK files for visualization output

vtkfile_u_1 = File('endothecial_cells_migration/u_1.pvd')
vtkfile_u_2 = File('endothecial_cells_migration/u_2.pvd')
vtkfile_u_3 = File('endothecial_cells_migration/u_3.pvd')

# Create progress bar
a = derivative(F, u, du)

problem = CellEquation(a, F)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6

# Time-stepping
t = 0
T = num_steps*dt
while (t < T):

    # Update current time
    t += dt

    u_n.vector()[:] = u.vector()
    solver.solve(problem, u.vector())
    vtkfile_u_1 << (u.split()[0], t)
    vtkfile_u_2 << (u.split()[1], t)
    vtkfile_u_3 << (u.split()[2], t)


# Hold plot
interactive ()

I don't understand what does not work. I tried multiple things...

Thank you for your time.