У меня есть динамическая модель, настроенная как (жесткая) система ODE. В настоящее время я решаю это с помощью CVODE (из пакета SUNDIALS в пакете python от Assimulo), и все хорошо.
Теперь я хочу добавить новый 3D-радиатор (с температурно-зависимыми температурными параметрами) к проблеме. Вместо того, чтобы записывать все уравнения с нуля для трехмерного уравнения теплоты, моя идея состоит в том, чтобы использовать существующую инфраструктуру FEM или FVM, чтобы предоставить мне интерфейс, который позволит мне легко предоставить (t, y)
для 3D-блока для обычной, и вернуть остатки y'
. Принцип заключается в использовании уравнений из системы FEM, но не в решателе. CVODE может использовать разреженность, но ожидается, что объединенная система будет работать медленнее, чем система FEM будет решать сама по себе, будучи адаптированной для таких целей.
# 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
Мой вопрос в том, является ли (a) этот подход разумным, и (b) есть библиотека FEM или FVM, которая легко позволит мне рассматривать ее как систему уравнений, такую, что я могу "применить ее" к моему существующий набор уравнений ODE.
Если две системы не могут использовать одну и ту же ось времени, тогда мне придется запускать их в режиме шага, где я запускаю одну модель на короткое время, обновляю граничные условия для другого, запускаю один, обновить первую модель BC и т.д.
У меня есть некоторый опыт работы с замечательной библиотекой FiPy, и я ожидаю, что в конечном итоге вы сможете использовать эту библиотеку так, как описано выше. Но я хочу знать об опыте с другими системами в таких проблемах, а также о других подходах, которые я пропустил.
Изменить: теперь у меня есть примерный код, который работает, показывая, как остатки диффузии сетки FiPy могут быть решены с помощью CVODE. Однако это только один подход (с использованием FiPy), а остальные мои другие вопросы и проблемы по-прежнему остаются. Любые предложения приветствуются.
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...')
Это даст график, показывающий идеальное соответствие: