2
$\begingroup$

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):

enter image description here

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.

$\endgroup$

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.