Solver did not converge on Steady-State Reaction-Diffusion equation

1

Hi all,
I’m a bit new to Fenics and I’m trying to solve this system of two coupled 2D stationary PDE:

\large\frac{\partial a(x,y,t)}{\partial t} = D_{a}\nabla^{2}a + -ab^{2} + f(1-a)

\large\frac{\partial b(x,y,t)}{\partial t} = D_{b}\nabla^{2}b + ab^{2} - (f+k)b

with \frac{\partial a(x,y,t)}{\partial t}=\frac{\partial b(x,y,t)}{\partial t}=0.

I’ve implemented the resolution in Fenics considering the vector (a,b) using the Newton Solver:

k = 0.05
f = 0.04

mesh = RectangleMesh(Point(0, 0), Point(100,100), 30, 30)

V = VectorFunctionSpace(mesh, 'P', 2)

X = Function(V)
Y = TestFunction(V)

#Initial Conditions
class IC(UserExpression):
        def eval(self,values,x):
            if (x[0]-50)**2+(x[1]-50)**2<=100:  #Small perturbation in the middle
                values[0] = 0.5 + 0.01*np.random.randn()
                values[1] = 0.25 + 0.01*np.random.randn()
            else:
                values[0] = 1.0
                values[1] = 0.0
        def value_shape(self):
            return(2,)
        
    X.interpolate(IC())
    (a,b) = split(X)
    (u,v) = split(Y)

    F = -Da*dot(grad(a), grad(u))*dx - a*b*b*u*dx + f*(1-a)*u*dx -Db*dot(grad(b),grad(v))*dx + a*b*b*v*dx - (f+k)*b*v*dx

    solve(F == 0, X, solver_parameters={"newton_solver":{"relative_tolerance":1e-3},"newton_solver":{"maximum_iterations":700}})

The thing is that I want to print in a few seconds those beautiful Turing’s patterns, but when I run this script, I usually got the error (sometime I can have a result btw):

RuntimeError                              Traceback (most recent call last)
<ipython-input-434-647364710ae2> in <module>

---> solve(F == 0, X, solver_parameters={"newton_solver":{"relative_tolerance":1e-2},"newton_solver":{"maximum_iterations":700}})

    ...  

RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to solve nonlinear system with NewtonSolver.
*** Reason:  Newton solver did not converge because maximum number of iterations reached.
*** Where:   This error was encountered inside NewtonSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

I’ve seen that many people got the same error, but by applying their solutions I couldn’t find mine (reducing the solver tolerance for example).

So I’m wondering if there are parameters to change or an other way to implement the problem that can be more efficient in this situation.

(I’m using the latest version of Fenics, and I’m running it on Windows 10 via WSL)

Thank you in advance for your precious help! :blush:

2

Fisher’s equation is pretty horrific at the best of times. And this looks even more horrific.

Are you sure there is a steady state solution? and that your initial guess is in an attractor region?

It may be worth just paying the computational cost and solving the time dependent problem until it reaches steady state to some satisfactory tolerance.

3

Thank you for your reply,

Indeed, I rushed to the steady-state solution, but considering the time dependant equation with t tends to infinity seemed to be quite efficient.

Thanks again for opening my mind!