How to combine an ODE system with a FEM system

Question

I have a dynamic model set up as a (stiff) system of ODEs. I currently solve this with CVODE (from the SUNDIALS package in the Assimulo python package) and all is good.

I now want to add a new 3D heat sink (with temperature-dependent thermal parameters) to the problem. Instead of writing out all the equations from scratch for the 3D heat equation, my idea is to use an existing FEM or FVM framework to provide to me an interface that will allow me to easily provide the (t, y) for the 3D block to a routine, and get back the residuals y'. The principle is to use the equations from the FEM system but not the solver. CVODE can exploit sparsity, but the combined system is expected to solve slower than the FEM system would solve on its own, being tailored for such.

# pseudocode of a residuals function for CVODE
def residual(t, y):

    # ODE system of n equations 
    res[0] = <function of t,y>;
    res[1] = <function of t,y>;
    ...
    res[n] = <function of t,y>;

    # Here we add the FEM/FVM residuals
    for i in range(FEMcount):
        res[n+1+i] = FEMequations[FEMcount](t,y)

    return res

My question is whether (a) this approach is sane, and (b) is there a FEM or FVM library that will easily let me treat it as a system of equations, such that I can "tack it on" to my existing set of ODE equations.

If can't let the two systems share the same time axis, then I will have to run them in a stepping mode, where I run the one model for a short time, update the boundary conditions for the other, run that one, update the first model's BCs, and so on.

I have some experience with the wonderful library FiPy, and I am expecting to eventually end up using that library in the manner described above. But I want to know about experience with other systems in problems of this nature, and also other approaches that I have missed.

---

Edit: I now have some example code that appears to be working, showing how the FiPy mesh diffusion residuals can be solved with CVODE. However, this is only one approach (using FiPy) and the remainder of my other questions and concerns still stand. Any suggestions welcome.

from fipy import *
from fipy.solvers.scipy import DefaultSolver
solverFIPY = DefaultSolver()

from assimulo.solvers import CVode as solverASSIMULO
from assimulo.problem import Explicit_Problem as Problem

# FiPy Setup - Using params from the Mesh1D example
###################################################
nx = 50; dx = 1.; D = 1.
mesh = Grid1D(nx = nx, dx = dx)
phi = CellVariable(name="solution variable", mesh=mesh, value=0.)
valueLeft, valueRight = 1., 0.

phi.constrain(valueRight, mesh.facesRight)
phi.constrain(valueLeft, mesh.facesLeft)

# Instead of eqX = TransientTerm() == ExplicitDiffusionTerm(coeff=D),
# Rather just operate on the diffusion term. CVODE will calculate the
# Transient side
edt = ExplicitDiffusionTerm(coeff=D)

timeStepDuration = 0.9 * dx**2 / (2 * D)
steps = 100

# For comparison with an analytical solution - again,
# taken from the Mesh1D.py example
phiAnalytical = CellVariable(name="analytical value", mesh=mesh)
x = mesh.cellCenters[0]
t = timeStepDuration * steps
from scipy.special import erf
phiAnalytical.setValue(1 - erf(x / (2 * numerix.sqrt(D * t))))
if __name__ == '__main__':
    viewer = Viewer(vars=(phi, phiAnalytical))#, datamin=0., datamax=1.)
    viewer.plot()

raw_input('Press a key...')

# Now for the Assimulo/Sundials solver setup
############################################

def residual(t, X):
    # Pretty straightforward, phi is the unknown
    phi.value = X # This is a vector, 50 elements
    # Can immediately return the residuals, CVODE sees this vector
    # of 50 elements as X'(t), which is like TransientTerm() from FiPy
    return edt.justResidualVector(var=phi, solver=solverFIPY)

x0 = phi.value
t0 = 0.
model = Problem(residual, x0, t0)
simulation = solverASSIMULO(model)
tfinal = steps * timeStepDuration # s,

cell_tol = [1.0e-8]*50
simulation.atol = cell_tol
simulation.rtol = 1e-6
simulation.iter = 'Newton'

t, x = simulation.simulate(tfinal, 0)

print x[-1]
# Write back the answer to compare
phi.value = x[-1]
viewer.plot()
raw_input('Press a key...')

This will produce a graph showing a perfect match:

Answer 1

An ODE is a differential equation in one dimension.

An FEM model is for problems that are spacial ie, problems in higher dimensions. You want a finite difference method. It's easier to solve and understand from the perspective someone coming from ODE world.

The question I think you should really be asking is how do you take your ODE, and transfer it to a 3D problem space.

Multidimensional partial differential equations are difficult to solve, yet I'll refer you to a CFD method for doing just that however only in 2D. http://lorenabarba.com/blog/cfd-python-12-steps-to-navier-stokes/

It should take you a solid afternoon to get through that! Good Luck!