How to implement this PDE system in FiPy?

Question

I'm solving a system of PDEs using FiPy which involves integration terms. However, the integration term is not allowed in FiPy. Can someone help me with this system?System of PDE

Answer 1

I'm not sure about this, but I think you'll have to keep track of the integration yourself.

You should be able to implement trapezoidal integration something like:

import fipy as fp

sigma = fp.CellVariable(..., hasOld=True)
phi = fp.CellVariable(..., hasOld=True)
integral = fp.CellVariable(..., value=0., hasOld=True)

sigma_eq = (fp.DiffusionTerm(var=sigma) - fp.ImplicitSourceTerm(coeff=phi, var=sigma) 
            == 0)
phi_eq = (fp.TransientTerm(var=phi)
          == fp.ImplicitSourceTerm(coeff=-integral + alpha, var=phi)

for step in steps:
    sigma.updateOld()
    phi.updateOld()
    integral.updateOld()
    for sweep in sweeps:
        sigma_res = sigma_eq.sweep(dt=dt)
        phi_res = phi_eq.sweep(dt=dt)
        integral.value = (integral.old 
                          + (sigma * phi + sigma.old * phi.old) * dt / 2).value

Instead of for sweep in sweeps, you'll actually want some sort of test of the convergence of sigma_res and phi_res. You'll need to experiment. This may be wildly unstable.