difficulties with solving the Gray Scott model

Question

Hi, I'm slightly new to fenics and I've been trying to solve the Gray-Scott model for some time. I've tried to adapt the Cahn-Hillard python demo file (found here: http://fenicsproject.org/documentation/dolfin/dev/python/_downloads/demo_cahn-hilliard.py) for this purpose and now I have trouble locating where my problem is. I assume it has some something to do with the weak form when I'm passing it to the Newton solver, i.e., I'm not quite sure I'm using the Test- and TrialFunctions in the right way.

Here is some of my code:

# Split mixed functions, where W: solution n+1
#                              W0: solution n
u,  v  = split(W)
u0, v0 = split(W0)

# Weak statement of the equations
# Implicit Euler/ Implicit time stepping, with q,p TestFunctions on the domain V*V
F1 = -u*q*dx + u0*q*dx \
      + D_u*inner(grad(u), grad(q))*dt*dx \
      - u*v0*v0*q*dt*dx \
      + c*(1-u)*q*dt*dx
F2 = -v*p*dx + v0*p*dx \
       + D_v*inner(grad(v), grad(p))*dt*dx \
       + u0*v*v*p*dt*dx \
       - (c+k)*v*p*dt*dx
F = F1 + F2

# Compute directional derivative about W in the direction of dp (Jacobian)
a = derivative(F, W, dp)

# Create nonlinear problem and Newton solver
problem = GrayScottEquations(a, F)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6

# Step in time
t = 0.0
T = 50*dt
while (t < T):
    t += dt
    W0.vector()[:] = W.vector()
    solver.solve(problem, W.vector())

plot(W.split()[1])
plot(W.split()[0])
interactive()

Comments

I solved the Gray-Scott model with FEniCS a while ago, so it is definitely possible ;-). But it is impossible to help you without more details. Why do you use different spaces W and W0? What is W anyway? What is the class GrayScottEquations? What are the exact errors which you observe?

---

Thank you very much for your reply. I choose only to post a part of my code and see this can be a problem now. This is the full code (with the problem stated below the code):

from dolfin import *
import numpy as np

# Set parameters
D_u = 8.0e-05; D_v = 4.0e-05; c = 0.024; k = 0.06; dt = 0.1;

# Class representing the intial conditions
class InitialConditions(Expression):
    def eval(self, val ,x):
        if between(x[0], (0.25, 0.75)) and between(x[1], (0.25, 0.75)):
            val[1] = 0.25*np.power(np.sin(4*np.pi*x[0]), 2)\
                        *np.power(np.sin(4*np.pi*x[1]), 2)
            val[0] = 1 - 2*val[1]
        else:
            val[1] = 0
            val[0] = 1
    def value_shape(self):
        return (2,)

# Class for interfacing with the Newton solver
class GrayScottEquations(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 mesh and function spaces
p0 = Point(0.0, 0.0)
p1 = Point(1.0, 1.0)
mesh = RectangleMesh(p0, p1, 49, 49)
V = FunctionSpace(mesh, 'CG', 1)
ME = V*V

# Define trial and test functions
dp = TrialFunction(ME)
q , p = TestFunctions(ME)

# Define functions
W   = Function(ME)  # current solution
W0  = Function(ME)  # solution from previous converged step

#Create intial conditions and interpolate
W_init = InitialConditions()
W0.interpolate(W_init)

# Split mixed functions
#du, dv = split(dp)
u,  v  = split(W)
u0, v0 = split(W0)

# Weak statement of the equations
# Implicit Euler/ Implicit time stepping
F1 = -u*q*dx + u0*q*dx + D_u*inner(grad(u), grad(q))*dt*dx \
        - u*v0*v0*q*dt*dx + c*(1-u)*q*dt*dx
F2 = -v*p*dx + v0*p*dx + D_v*inner(grad(v), grad(p))*dt*dx \
        + u0*v*v*p*dt*dx - (c+k)*v*p*dt*dx
F = F1 + F2

# Compute directional derivative about W in the direction of dp (Jacobian)
a = derivative(F, W, dp)

# Create nonlinear problem and Newton solver
problem = GrayScottEquations(a, F)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6

# Step in time
t = 0.0
T = 50*dt
while (t < T):
    solver.solve(problem, W.vector())
    W0.vector()[:] = W.vector()
    t += dt

plot(W.split()[1])
plot(W.split()[0])
interactive() 

My exact problem is, when I'm trying to solve the model with the Newton solver. It seems like it isn't solving the problem properly (I either get wierd oscillations or the solver doesn't converge at all). This, in turn makes me think that the way I'm using the split functions (u, v, u0, v0) in the weak form is wrong i.e., the weak form isn't started properly in regards to the functions W and W0 defined on the functionspace ME

Regards Nicklas.

---

Did you try smaller timesteps?

---

When trying smaller timesteps the Newton solver converges, but sadly I get big ocillations in the solution.

---

I am quite busy now and cannot go into details of your implementation... Here is mine from a couple of years ago (might need to be adjusted to work with current versions)

n     = 128           # mesh parameter
p     = 1             # polynomial degree
dt    = 1.0           # time step
gamma = 0.024         # feed rate
kappa = 0.06          # reaction rate
D     = (8e-5, 4e-5)  # diffusivities

# create mesh
mesh = RectangleMesh(0, 0, 2.5, 2.5, n, n)

# initial conditions
class Initial(Expression):
  def eval(self, value, x):
    value[1] = 0
    if between(x[0], (1, 1.5)) and between(x[1], (1, 1.5)):
      value[1] = (sin(4*pi*x[0])**2)*(sin(4*pi*x[1])**2)/4
    value[0] = 1 - 2*value[1]
    return
  def value_shape(self):
    return (2,)
initial = Initial()

# diffusion terms
def a(u, v):
  return inner(D[0]*grad(u[0]), grad(v[0]))*dx \
       + inner(D[1]*grad(u[1]), grad(v[1]))*dx

# reaction/replenishment/diminishment terms
def f(u, v):
  return (u[0]*u[1]**2)*(v[1] - v[0])*dx \
       + gamma*(1 - u[0])*v[0]*dx \
       - (gamma + kappa)*u[1]*v[1]*dx

V = VectorFunctionSpace(mesh, 'Lagrange', p)
U = interpolate(initial, V)
u = Function(V)
v = TestFunction(V)

# variational problem
F = dot(u,v)*dx - dot(U,v)*dx + dt*(a(u,v) - f(u,v))

# set up nonlinear problem
du = TrialFunction(V)
dF = derivative(F, u, du);
problem = NonlinearVariationalProblem(F, u, None, dF)

# set up nonlinear and linear solver
solver = NonlinearVariationalSolver(problem)
solver.parameters["linear_solver"]  = "cg"  # conjugate gradient method
solver.parameters["preconditioner"] = "ilu" # incomplete LU preconditioner
solver.parameters["newton_solver"]["maximum_iterations"] = 25

T = 2000 # total simulation time

ufile = File("output/gray-scott.pvd", "compressed")
ufile << U

# time-stepping loop
t = 0.0
u.assign(U)
while t <= T:
  t += dt
  solver.solve()
  U.assign(u)
  print 'time =', t
  ufile << u

Hope this helps to identify your error.

---

Thank you very much :-). This helps a lot. I'll try to locate my error and post the result here for anyone else trying to adapt the Cahn-Hillard python demo to the Gray-Scott model.

---

Okay I forgot to post this yesterday. I found the error to be:
I was using the wrong function space and I had a simple error in my weak variational formulation. Here's the corrected code:

from dolfin import *
import numpy as np

# Set parameters
D_u = 8.0e-05; D_v = 4.0e-05; c = 0.024; k = 0.06;  dt = 1.0;

# Class representing the intial conditions
class InitialConditions(Expression):
    def eval(self, val ,x):
        if between(x[0], (0.25, 0.75)) and between(x[1], (0.25, 0.75)):
            val[1] = 0.25*np.power(np.sin(4*np.pi*x[0]), 2)*np.power(np.sin(4*np.pi*x[1]), 2)
            val[0] = 1 - 2*val[1]
        else:
            val[1] = 0
            val[0] = 1
    def value_shape(self):
        return (2,)

# Class for interfacing with the Newton solver
class GrayScottEquations(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 mesh and function space
p0 = Point(0.0, 0.0)
p1 = Point(1.0, 1.0)
mesh = RectangleMesh(p0, p1, 28, 28)
V = VectorFunctionSpace(mesh, 'CG', 1)

# Define functions 
W_init = InitialConditions()
phi    = TestFunction(V)
dp     = TrialFunction(V)
W0     = Function(V)
W      = Function(V)
# Interpolate initial conditions and split functions
W0.interpolate(W_init)
q, p   = split(phi)
u,  v  = split(W)
u0, v0 = split(W0)

# Weak statement of the equations
F1 = u*q*dx - u0*q*dx + D_u*inner(grad(u), grad(q))*dt*dx + u*v*v*q*dt*dx - c*(1-u)*q*dt*dx
F2 = v*p*dx - v0*p*dx + D_v*inner(grad(v), grad(p))*dt*dx - u*v*v*p*dt*dx + (c+k)*v*p*dt*dx
F = F1 + F2

# Compute directional derivative about W in the direction of dp (Jacobian)
a = derivative(F, W, dp)

# Create nonlinear problem and Newton solver
problem = GrayScottEquations(a, F)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-5

t = 0.0
T = 200
W.assign(W0)
while (t < T):
    t += dt
    solver.solve(problem, W.vector())
    W0.assign(W)
    print 'time =', t    

plot(W.split()[1])
plot(W.split()[0])
interactive()

(As I've pointed out ealier, this is an adapted version of the Cahn-Hillard python demo file)

Thanks a lot for the help Christian!