gpt4 book ai didi

python - 如何将 ODE 系统与 FEM 系统相结合

转载 作者:IT老高 更新时间:2023-10-28 20:29:22 25 4
gpt4 key购买 nike

我有一个动态模型设置为 ODE 的(刚性)系统。我目前用 CVODE(来自 Assimulo python 包中的 SUNDIALS 包)解决了这个问题,一切都很好。

我现在想为问题添加一个新的 3D 散热器(具有与温度相关的热参数)。我的想法是使用现有的 FEM 或 FVM 框架为我提供一个接口(interface),让我可以轻松地提供 (t, y) 将 3D block 转换为例程,并取回残差 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 有一些经验,我希望最终以上述方式使用该库。但我想了解其他系统在此类问题上的经验,以及我错过的其他方法。


编辑: 我现在有一些示例代码似乎正在运行,展示了如何使用 CVODE 解决 FiPy 网格扩散残差。但是,这只是一种方法(使用 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...')

这将产生一个显示完美匹配的图表:

Perfect match

最佳答案

ODE 是一维微分方程。

FEM 模型适用于空间问题,即更高维度的问题。你想要一个有限差分法。从来自 ODE 世界的人的角度来看,更容易解决和理解。

我认为您真正应该问的问题是如何获取 ODE,并将其转移到 3D 问题空间。

多维偏微分方程很难求解,但我会推荐您使用 CFD 方法来解决此问题,但仅限于二维。 http://lorenabarba.com/blog/cfd-python-12-steps-to-navier-stokes/

这需要你一个充实的下午才能完成!祝你好运!

关于python - 如何将 ODE 系统与 FEM 系统相结合,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/21565517/

25 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com