2D reaction-diffusion system simulation
I am a complete beginner in numerical simulation and I am pretty lost about how to tackle this problem.
I have been trying for some time to find the steady state (or simulate), the following system
$$ \dot S = D_S \Delta S + \rho_S \frac{ F^n } { F^n + 1} - \delta_S S -\delta_C P^2 S $$
$$ \dot P = D_P \Delta P + \rho_P - \delta_P P + (v - 2 \delta_C) P ^ 2 S $$
$$ \dot F = \Delta + \rho_F \frac{ 1 } { (P^2S)^n + 1 } - \delta_F F $$
on a mesh that looks like this (with no boundary conditions given):
I realize that this systems contains the temporal derivatives, so I'm not explicitly going for the steady state, but when I took out the TransientTerm it complained that it was Exactly singular.
S = CellVariable(mesh = mesh, hasOld=True)
P = CellVariable(mesh = mesh, hasOld=True)
F = CellVariable(mesh = mesh, hasOld=True)
sdot = TransientTerm(var=S) == DiffusionTerm(coeff=Ds, var=S) + \
rs * pow(F,2) / (pow(F,2) + 1) - ds * S - dc * pow(P, 2) * S
pdot = TransientTerm(var=P) == DiffusionTerm(coeff=Dp, var=P) + \
rp - dp * P + (v - 2. * dc) * S * pow(P, 2)
fdot = TransientTerm(var=F) == DiffusionTerm(coeff=1, var=F) + \
rf / (pow(S * pow(P, 2), 2) + 1.0 ) - df * F
dt = 0.005
coupled.sweep(dt = dt)
for i in range(0, 100000):
v1.plot('plots/SHH_'+str(dt * i)+'.png')
v2.plot('plots/Ptc_' + str(dt * i)+'.png')
v3.plot('plots/FGF_'+str(dt * i)+'.png')
S.updateOld()
P.updateOld()
F.updateOld()
dt = 0.005
for i in range(0, 100000):
v1.plot('plots/SHH_'+str(dt * i)+'.png')
v2.plot('plots/Ptc_' + str(dt * i)+'.png')
v3.plot('plots/FGF_'+str(dt * i)+'.png')
S.updateOld()
P.updateOld()
F.updateOld()
coupled.sweep(dt = dt)
The simulation initially appears to be proceeding as normal, but eventually, regardless of the size of my time-step, the $S$ values and $P$ values drop below zero and subsequently $S$ blows up do its large production rate $\rho_s = 1600$, causing an overflow error.
I should mention that when I simulated this problem in 1D using pdepe, the simulation was find and converged to a steady state.
How do I debug this, and is there other problems with my simulation? Should I try a different package such as Firedrake? If so, how do I solve the convergence problems there?
Thanks very much in advance. 935