Newbie questions around solving the brusselator equations

[vasilvas]

Hi, I’ve gone through some of the examples for FEniCSx and I’m exploring solving the Brusselator equations as they does not deviate too much from the heat equation/Cahn-Hilliard examples, yet have some nonlinear source terms that I thought would be a good way to learn more about FEniCSx.

As a reference I am using this:
2.14. Brusselator - Using the PDE class — py-pde unknown documentation

Long story short, my code is mostly taken from the Cahn-Hilliard equation example for Dolfinx:

#!/usr/bin/env python3

import numpy as np
import ufl
from dolfinx import fem, io, mesh
from dolfinx.fem import Function, FunctionSpace, Expression
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from ufl import dot, dx, grad

# Define temporal parameters
t = 0  # Start time
T = 20  # Final time
num_steps = 200
dt = T / num_steps  # time step size

A = 1
B = 3

Dc = 1  # diffusion coefficient for c
Dmu = 0.1  # diffusion coefficient for mu

msh = mesh.create_unit_square(MPI.COMM_WORLD, 128, 128, mesh.CellType.triangle)
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
ME = FunctionSpace(msh, P1 * P1)


q, v = ufl.TestFunctions(ME)

u = Function(ME)  # current solution
u0 = Function(ME)  # solution from previous converged step

c, mu = ufl.split(u)  # split vector into components u = (c, mu)
c0, mu0 = ufl.split(u0)

# zero-out array
u.x.array[:] = 0.0

# i.c.
u.sub(0).interpolate(lambda x: A + x[0]*0.0)
u.sub(1).interpolate(lambda x: B / A + A*np.random.rand(x.shape[1]))


u.x.scatter_forward()

react_c = A - (1 + B) * c + mu * c * c
react_mu = B * c - mu * c * c


F0 = (
    c * q * dx
    + Dc * dt * dot(grad(c), grad(q)) * dx
    - c0 * q * dx
    - dt * react_c * q * dx
)
F1 = (
    mu * v * dx
    + Dmu * dt * dot(grad(mu), grad(v)) * dx
    - mu0 * v * dx
    - dt * react_mu * v * dx
)
F = F0 + F1

problem = NonlinearProblem(F, u)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
ksp.setFromOptions()

file = io.XDMFFile(MPI.COMM_WORLD, "bruss_output.xdmf", "w")

file.write_mesh(msh)
t = 0.0

c = u.sub(0)
mu = u.sub(1)
u0.x.array[:] = u.x.array

while t < T:
    t += dt
    r = solver.solve(u)
    print(f"Step {int(t/dt)}: num iterations: {r[0]}")
    u0.x.array[:] = u.x.array
    file.write_function(c, t)
    file.write_function(mu, t)

file.close()

Based on it I have stumbled onto the following problems:

1. All I ever see is the diffusion of the initial conditions, no oscillatory behaviour, even though the i.c. is close to the equilbirium point + some random perturbance. I think my variational formulation is alright as it’s a simple reaction-diffusion equation?

2. If I do not reference x in the initial condition lambda-function i get a seg-fault. That’s why i have that nasty 0.0*x[0] hack there. Why is this the case?

3. Is this really the best to handle this system of equations with FEniCSx? I have done an implementation of the brusselator “manually” before and the system matrix is constant, all that ever changes is the “load vector”. This way it feels that it would assemble the whole system matrix on each iteration step (=> slow). What can I do better?

[dokken]

@vasilvas

Is this really the best to handle this system of equations with FEniCSx? I have done an implementation of the brusselator “manually” before and the system matrix is constant, all that ever changes is the “load vector”. This way it feels that it would assemble the whole system matrix on each iteration step (=> slow). What can I do better?

Surely, if the system is non-linear, the matrix (jacobian) cannot be constant, as it will depend on the previous state. When you solve a non-linear problem, the most common thing is to use a Newton-solver, as described in:
https://jsdokken.com/dolfinx-tutorial/chapter4/newton-solver.html
Then you need the jacobian evaluated at the previous newton iteration.

[vasilvas]

I started writing-out the weak formulation to better explain myself, and then I realized where’s my mistake - that the non-linear source depends on the concentration of the components…
Blah, I feel so dumb now. Anyway, thanks for poiting me in the right direction!

Do you perhaps have an idea about the seg-fault issue?

[dokken]

Im not at a laptop, so i cant run your code at the moment. I will have a look tomorrow if I find time for it.

Note that in general; convergence of newton methods are local, so if your initial guess is 0, and this is far away from the solution, it might diverge.

@nate has written a good guide on non linear problems:

Default absolute tolerance and relative tolerance

I write the following in the context of FE modelling rather than numerical analysis. I.e. this should be interpreted as tips and tricks. Perhaps someone else can contribute something more rigorous, or their own experience. Do you know what to expect the residual of your system to be to machine precision? Adjust your absolute tolerance accordingly. E.g: A Poisson problem with a diffusion coefficient of 1 may yield a residual to machine precision around 10^{-12}. Linear elasticity with a Young’…

[vasilvas]

Thank you for the support!
Edit: Found the issue, i was defining an unit square, not a 128x128 square with create_unit_square (of course ). Fixed that, changed to a quad element and made it break if it fails to converge on the solver. Now this produces the expected behaviour:

#!/usr/bin/env python3

import ufl

import numpy as np

from dolfinx import io, mesh
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.fem import Function, FunctionSpace

from mpi4py import MPI
from petsc4py import PETSc
from ufl import dot, dx, grad

# Define temporal parameters
t = 0  # Start time
T = 25  # Final time
num_steps = 3000
dt = T / num_steps  # time step size

A = 1
B = 3

Dc = 1  # diffusion coefficient for c
Dmu = 0.1  # diffusion coefficient for mu

msh = mesh.create_rectangle(
    comm=MPI.COMM_WORLD,
    points=[(0, 0), (64, 64)],
    n=[64, 64],
    cell_type=mesh.CellType.quadrilateral,
)
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
ME = FunctionSpace(msh, P1 * P1)


q, v = ufl.TestFunctions(ME)

u = Function(ME)  # current solution
u0 = Function(ME)  # solution from previous converged step

c, mu = ufl.split(u)  # split vector into components u = (c, mu)
c0, mu0 = ufl.split(u0)

# zero-out array
u.x.array[:] = 0.0

# i.c.
u.sub(0).interpolate(lambda x: A + 0.0 * x[0])
u.sub(1).interpolate(
    lambda x: B / A + 0.1 * A * np.random.standard_normal(size=x.shape[1])
)


u.x.scatter_forward()

react_c = A - (1 + B) * c + mu * c * c
react_mu = B * c - mu * c * c


F0 = (
    c * q * dx
    + Dc * dt * dot(grad(c), grad(q)) * dx
    - c0 * q * dx
    - dt * react_c * q * dx
)
F1 = (
    mu * v * dx
    + Dmu * dt * dot(grad(mu), grad(v)) * dx
    - mu0 * v * dx
    - dt * react_mu * v * dx
)
F = F0 + F1

problem = NonlinearProblem(F, u)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-10

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
ksp.setFromOptions()

file = io.XDMFFile(MPI.COMM_WORLD, "bruss_output.xdmf", "w")

file.write_mesh(msh)
t = 0.0

c = u.sub(0)
mu = u.sub(1)
u0.x.array[:] = u.x.array

file.write_function(c, t)
file.write_function(mu, t)

while t < T:
    t += dt
    r = solver.solve(u)
    print(f"Step {int(t/dt)}: num iterations: {r[0]}")
    if not r[1]:
        print("Solver failed to converge! Exiting")
        break
    u0.x.array[:] = u.x.array
    file.write_function(c, t)
    file.write_function(mu, t)

file.close()

Only issue left is the segfault