Solve steady state reaction-diffusion/Helmholtz equation numerically

Question

I am solving a problem of the form:

$\dfrac{\partial u(x,y,t)}{\partial t} = \nabla^2 u(x,y,t) - f(x,y,t)u(x,y,t) - \kappa(x,y,t)$

At the moment, I am solving this at each time step by assuming a quasi-steady-state:

$\nabla^2 u(x,y,t) = f(x,y,t)u(x,y,t) + \kappa(x,y,t)$

To do this, I use finite differences, which means solving the following (forgetting boundary conditions for now, but for reference they are Neumann) at each time step:

$\dfrac{u^t_{i, j+1} + u^t_{i+1, j} + u^t_{i, j-1} + u^t_{i-1, j} - 4u^t_{i,j}}{h^2} = f^t_{i,j} u^t_{i,j} + \kappa^t_{i,j}$

This means I have to solve a linear system of the form:

$Ax = b$

Where $A$ is a matrix made of the laplacian, as well as the $f^t_{i,j} u^t_{i,j}$ term. For my 100x100 system (resulting in a laplacian of size 10000x10000), this takes about 0.4 seconds for each time step using Eigen's C++ library solver:SimplicialLLT. However, I would like to speed up solving by any means possible. The issue is that I am limited in the pre-conditioning that I can do since the matrix $A$ changes each time step.

Does anyone have any ideas? Finite element, spectral methods, non-steady state approximations all welcome.

Before anyone asks, the terms $f^t_{i,j}$ and $\kappa^t_{i,j}$ cannot be forecast ahead of time, as they depend on $u^{t-1}_{i,j}$. Hence parallelisation is not going to be possible I think.

Best,

Ben

Answer 1

One approach to acceleration would go as follows: Assume that $A^t=A(u^t)$ is the matrix you try to solve with, i.e., you are looking to solve the linear systems $$ A^t x^t = b^t. $$ Let me assume for a moment that $f^t \ge 0$, then $A^t$ is a symmetric and positive definite matrix. (If my assumption should be wrong, then it is still symmetric but may no longer be positive definite.) Given that, you can solve the linear system with the Conjugate Gradients (CG) method. The question of course is how to precondition it.

My approach would be to compute a Cholesky decomposition of $A^0$ in the first time step and use this as the preconditioner for the next time steps. You will be able to solve time step zero in exactly one CG iteration since the preconditioner is perfect. In subsequent time steps, $A^t$ will start to deviate from the preconditioner and you will need more iterations. Once the number of iterations becomes too large, you re-compute the Cholesky decomposition in the next time step. This works because the expensive step is to compute the Cholesky decomposition, but you only do that every few time steps.