2011/6/2 Jonathan Guyer <[email protected]> > > > On Jun 1, 2011, at 10:55 AM, Jeremy Laforet wrote: > > > > > I've been looking to use external software to simulate the model of > > uterine muscle we developed in our team. > > FiPy really catch my eye for its high level approach. > > But even after reading the documentation and searching the mailing > > list archive, I'm unsure if FiPy is able to solve this type of problem. > > > > The actual formulation of the model is using 3 variables (per cell): > > - their temporal evolution is described by non-linear coupled ODEs > > (one for each variable) > > - anisotropic diffusion occurs on one of the variables. > > > > Each cell of the mesh would account for a muscle cell, so my variables > > would be CellVariables. I think the ODEs would appear as Source terms. > > But I still don't see how to write my model's equations in the right > > way for FiPy. > > > > Is FiPy able to simulate such things ? If yes, could you point me to > > some examples (regarding multivariables) ? > > As we've just been discussing in > > http://thread.gmane.org/gmane.comp.python.fipy/2208 > > FiPy should certainly be able to do this, although we do not have dedicated > ODE solvers. > > One option is to use the ODE solvers in SciPy for that part of the problem > and manually assign those solution variables to your FiPy sources. > > The other option, as Biswa has done in his script, is to treat the ODEs as > degenerate PDEs. A combination of TransientTerm and SourceTerms will solve > the same ODE at every point in your domain and you can then use that > solution as the coefficients to your diffusive PDE. > > Let's say you are solving > > \frac{\partial \phi}{\partial t} = \nabla\cdot(m \nabla \phi) > > and > > \frac{d m}{d t} = -\alpha m > > then for the second case, you would write something like > > phi = CellVariable(mesh=mesh, value=phi0, hasOld=True) > m = CellVariable(mesh=mesh, value=m0, hasOld=True) > > phiEq = TransientTerm() == DiffusionTerm(coeff=m) > mEq = TransientTerm() == ImplicitSourceTerm(coeff=-alpha) > > for step in range(steps): > phi.updateOld() > m.updateOld() > for sweep in range(sweeps): > mEq.sweep(var=m, dt=dt) > phiEq.sweep(var=phi, dt=dt) > > > If you were to use SciPy (which might be a very good idea, as it has a > range of specialized ODE solvers), then I'd be inclined to write > > m = Variable(value=m0) > > and then, whenever you solve the ODE, do something like > > soln = odeint(mFn, m0, ...) > > m.setValue(soln[-1, 0]) >
Actually, the solvers in scipy start to show their age if you need more advanced features like preconditioning for the ODEs. The solvers in scipy have evolved into the C++ Sundials package https://computation.llnl.gov/casc/sundials/main.html Python bindings are unfortunately not being updated at the moment http://sourceforge.net/projects/pysundials/ (stuck on sundials 2.3) I wrote a high level interface to pysundials which does what I need (which is a subset of what is possible): http://scikits.appspot.com/odes, advantage of that is that you don't need to learn how sundials works as it is all hidden by the interface. Anyway, odeint is the old lsode solver, if you use ode from scipy you use vode. If you use pysundials or my package, you use cvode or idas, so you have something like ode15s in matlab. I advise you to use an ode solver and not pde, as the ode solvers have adaptive time stepping and error control, things you don't get in fipy. We combine ode with fipy, where we solve 1D problems with ode, and use fipy for higher dimensions. Benny > > to assign the ODE solution back to the FiPy Variable m. > > I recommend you familiarize yourself with the ODE examples in > http://www.scipy.org/Cookbook and then ask specific questions about how to > proceed. > > > > >
