Choosing a solver for a non-linear reaction-diffusion equation

1

Hi

I’m not very familiar with FEniCS and have encountered a problem when trying to solve a system of advection-diffusion equations. I can’t seem to get the system to converge. I have tried changing the solver and providing a good initial guess (which I gauged from an image off of the paper I copied the equation from), but nothing seemed to work.

Here is the problem:

u2v = c_3*u_1*u_2*u_1
F1 = D_u*inner(grad(u_1), grad(v_1))*dx - (c_1*v_1 - c__1*u_1*v_1 + u2v*v_1)*dx
F2 = D_v*inner(grad(u_2), grad(v_2))*dx - (c_2*v_2 - u2v*v_2)*dx

F = F1 + F2

where c_ are all constants, u_ are the trial functions (though they come from Function) and v_ are the TestFunctions.

This is the solver at the moment:

J = derivative(F, u)
problem = NonlinearVariationalProblem(F_init, u, [], J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters

prm['nonlinear_solver'] = 'newton'
prm['newton_solver']['linear_solver'] = 'gmres'
prm['newton_solver']['preconditioner'] = 'jacobi'
prm['newton_solver']['relaxation_parameter'] = 1.0
prm['newton_solver']['krylov_solver']['maximum_iterations'] = int(1e3)
prm['newton_solver']['krylov_solver']['absolute_tolerance'] = 1e-20

and the initial guess:

u_init = Expression(('2.9*pow(sin(.25*x[0])*sin(.25*x[1]), 10)', '0.25*pow(cos(.25*x[0])*cos(.25*x[1]), 10)'), degree=3)

But all I manage to get is the error message:

*** Error:   Unable to solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 1000 iterations (PETSc reason DIVERGED_ITS, residual norm ||r|| = 7.644914e-01).

I would appreciate any help in understanding how I can debug my problem or known solutions. Thank you very much for your time.

2

Hi paulvo,

The error message implies that your implementation solves for some iterations. This means that your code is working, and it is a good sign!

I do not have much experience in solving pattern-forming reaction-diffusion equations. Though I can point out, usually, there is a time derivative term in such problems. Your implementation is now solving a steady-state problem, which may not converge.

So I would suggest you define a Function u_10 and u_20
u_10 = interpolate(your_Expression1, your_function_space)
u_20 = interpolate(your_Expression2, your_function_space).

Then add the time derivative term in your form, e.g.,
F1 += v_1*(u_1-u_10)/dt*dx.

You should look for other FEniCS tutorials about solving time-dependent problems. I hope this is helpful!

1 Like
3

Hi @Powei,

Thank you so much for your time.

Sorry for the late reply. I understand your thinking but my objective is to find the steady state (which I know to exist) so that I can study the behaviour of the time-dependent problem with an initial guess near to the steady state. That is why the term \partial_t u has been set to 0 for both terms.

In essence I would like to do what you suggest but setting u_10 and u_20 to be close to the steady state.

Do you know how I could achieve this goal?

4

Hi @paulvo,

I would suggest defining an initial condition (a random field, for example), then use fixed-point iterations to reach a good enough steady-state solution. I found a tutorial that implements such methods: Link.

You can also try to implement Anderson acceleration, which may speed up the rate of finding steady-state solutions.

I wish you a good day. Happy solving PDEs!

5

Hi @Powei,

thank you very much for your reply. I will look into the proposed methods, which will hopefully work.

6

Hi @Powei,

it turns out that the problem of divergence persists. I have tried altering a variety of parameters but nothing seems to help. Do you think that in this case the only reasons can be the type solver or the initial guess?

Thank you so much for your time and patience.