Best method to solve this system of PDEs?

Question

I have a system of PDEs constituting an initial value problem (IVP) consisting of three coupled PDEs:

\begin{align} \partial_t \rho + \partial_x(\rho v) &= \left(k_A (1-\phi) + k_B \phi \right)\rho\ \\ \partial_t \phi + v \partial_x{\phi} &= D \partial^2_x \phi + D \frac{\partial_x \rho \partial_x \phi}{\rho} + \left(k_B - k_A \right) \phi(1-\phi) + k_D (1-\phi) \\ \eta \nabla^2 v - \gamma v &= E \frac{1}{\rho} \partial_x \rho \end{align} Here,
\begin{align} E &= \left( E_A + (E_B - E_A) \phi \right) \rho \\ k_D (1-\phi) &= \alpha (E_B - E_A) \phi(1 - \phi) \\ k_i(\rho) &= \frac{1}{\tau} \frac{\rho^c_i-\rho}{\rho^c_i} \hspace{2cm} \text{(i=A, B)}, \end{align} and the rest of the variables are constants.

The equations come from continuum mechanics applied to biological tissues, in this case composed of two different cell types A and B. The unknown variables to solve for are $\rho(x, t), \phi(x,t)$ and $v(x,t)$. $\rho$ represents a density field, $0 \leq \phi \leq 1$ represents a fraction of cells belonging one of the two cell types and $v$ is the velocity of the medium. $k_A, k_B, k_D$ are reaction rates, $D$ is a diffusion constant, and $\eta, \gamma, E$ are material constants. $\rho_i^{c}$ represents a homeostatic density and $\alpha, \tau$ are proportionality constants.

I take Dirichlet boundary conditions for all variables, $\rho(x=-L) = \rho(x=L) = \rho_0$, $\phi(x=-L) \approx 1, \phi(x=L) \approx 0$ (such that it is consistent with $\phi(x, t=0)$) and $V(x=-L)=v(x=L)=0$. Initial conditions are $\rho(x, t=0)= \rho_0$, $\phi(x, t=0)$ is a sigmoid function going from 1 to 0 and $V(x, t=0)=0$.

I've been looking into some high-level methods to solve this IVP, but none of them work well enough. The best so far has been Mathematica (link to mathematica.stackexchange post), which works to some extent. However, the solutions are highly unstable. For some parameters I get a nice solution, but when I change some parameters Mathematica gets stuck in some infinite loop and never finishes the computation.

Now I'm learning a finite element solver (https://fenicsproject.org/), but I'm not even sure whether in principle this system can be solved using finite elements. Could someone confirm whether this is the right approach?

Or should I be looking at a more low-level approach and program the numerics myself?

My guess is that some methods for solving the Navier-Stokes equations or those for elasticity might also be applicable here, but I'm totally inexperienced with numerically solving PDEs so correct me if I'm wrong.

Answer 1

On a cursory look it seems you are solving a advection diffusion reaction equations. The first one seems to be conservation of mass while the second one is AD equation for variable phi. You can always discretize using finite differences and use upwind schemes for advection terms,central differences for diffusion and a second order implicit or a semi implicit scheme like crank Nicolson for transient term should work. Since we have cross terms and it looks nonlinear a stability analysis of linearized version can help you derive a CFL like condition. I have written a second order accurate axi symmetric code for some industrial heterogeneous catalytic systems and the above said method worked okay though the source terms were a tad simpler. My understanding of FE is that it has some issues in handling advection terms and hence nearly all commerical fluid codes are erstwhile of FD type and now of Finite volume which uses a conservative form of naiver stokes equations.