python - solve - Cómo combinar un sistema ODE con un sistema FEM

python differential equation library

Question

Tengo un modelo dinámico configurado como un sistema (rígido) de EDO. Actualmente resuelvo esto con CVODE (del paquete SUNDIALS en el paquete Assimulo python) y todo está bien.

Ahora quiero agregar un nuevo disipador de calor 3D (con parámetros térmicos dependientes de la temperatura) al problema. En lugar de escribir todas las ecuaciones desde cero para la ecuación de calor 3D, mi idea es utilizar un marco FEM o FVM existente para proporcionarme una interfaz que me permita proporcionar fácilmente (t, y) el bloque 3D para una rutina, y recuperar los residuos y'' . El principio es usar las ecuaciones del sistema FEM pero no el solucionador. CVODE puede explotar la escasez, pero se espera que el sistema combinado se resuelva más lentamente de lo que el sistema FEM resolvería por sí solo, y está diseñado para eso.

# 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

Mi pregunta es si (a) este enfoque es sensato, y (b) existe una biblioteca FEM o FVM que me permita tratarla como un sistema de ecuaciones, de manera que pueda "ajustarlo" a mi conjunto existente de Ecuaciones ODE.

Si no puedo permitir que los dos sistemas compartan el mismo eje de tiempo, entonces tendré que ejecutarlos en un modo paso a paso, donde ejecuto el modelo por un corto tiempo, actualizo las condiciones de contorno para el otro, ejecuto ese, actualice El primer modelo de BC, y así sucesivamente.

Tengo algo de experiencia con la maravillosa biblioteca FiPy, y estoy esperando que termine usando esa biblioteca de la manera descrita anteriormente. Pero quiero saber acerca de la experiencia con otros sistemas en problemas de esta naturaleza, y también otros enfoques que he omitido.

Edición: ahora tengo un código de ejemplo que parece estar funcionando, mostrando cómo se pueden resolver los residuos de difusión de la malla FiPy con CVODE. Sin embargo, este es solo un enfoque (usando FiPy) y el resto de mis otras preguntas y preocupaciones siguen en pie. Cualquier sugerencia de bienvenida.

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...'')

Esto producirá un gráfico que muestra una coincidencia perfecta:

Answer 1

Una EDO es una ecuación diferencial en una dimensión.

Un modelo FEM es para problemas que son espaciales, es decir, problemas en dimensiones más altas. Quieres un método de diferencia finita. Es más fácil de resolver y entender desde la perspectiva de alguien que viene del mundo ODE.

La pregunta que creo que realmente deberías estar preguntando es cómo tomas tu EDO y la transfieres a un espacio de problemas 3D.

Las ecuaciones diferenciales parciales multidimensionales son difíciles de resolver, sin embargo, te referiré a un método CFD para hacerlo solo en 2D. http://lorenabarba.com/blog/cfd-python-12-steps-to-navier-stokes/

¡Debería llevarte una tarde sólida para superar eso! ¡Buena suerte!