我目前正在使用扫描循环在 python 中使用 FiPy来求解关于我的单元变量phi的微分方程（ eq0 ）。因为我的方程是非线性的，所以我使用了一个扫描循环，如下面的代码摘录所示。
while res0 > resphi_tol:
res0 = eq0.sweep(var=phi, dt=dt)
但我不断收到以下错误：
C:\Python27\lib\site-packages\fipy\variables\variable.py:1100：RuntimeWarning：power return self._BinaryOperatorVariable(lambda a,b: pow(a,b), other, value1mattersForUnit=True )
C:\Python27\lib\site-packages\fipy\variables\variable.py:1186: RuntimeWarning: 在less_equal 中遇到无效值返回 self._BinaryOperatorVariable(lambda a,b: a<=b, other)
Traceback (最近最后调用）：
.. 文件“SBM_sphere3.py”，第 59 行，在
....res0 = eq0.sweep(var=phi, dt=dt)
.. 文件“C:\Python27\lib\site-packages\ fipy\terms\term.py"，第 207 行，扫描
....solver._solve()
.. 文件 "C:\Python27\lib\site-packages\fipy\solvers\pysparse\pysparseSolver.py"，行68，在_solve
....self。solve (self.matrix, array, self.RHSvector)
.. 文件“C:\Python27\lib\site-packages\fipy\solvers\pysparse\linearLUSolver.py”，第 53 行，在 _solve__
....LU = superlu .factorize(L.matrix.to_csr())
.. 文件“C:\Python27\lib\site-packages\pysparse\misc__init__.py”，第 29 行，在 newFunc
....return func(*args, ** kwargs)
.. 文件“C:\Python27\lib\site-packages\pysparse__init__.py”，第 47 行，因式分解
....return self.factorizeFnc(*args, **kwargs)
RuntimeError: Factor is exactly single
我很确定这个错误是由于eq0中存在的术语 phi^(2/3) 造成的。如果我用 abs(phi)^(2/3) 替换这个术语，错误就会消失。
我假设扫描循环在某些时候为phi中的几个单元格返回负值，导致错误，因为我们不能用非整数指数为负值供电。
所以我的问题是：有没有办法强制扫描以避免负面解决方案？
我试图包含一条在扫描之前将所有负值设置为 0 的行：
while res0 > resphi_tol:
phi.setValue(0.,where=phi<0.)
res0 = eq0.sweep(var=phi, dt=dt)
错误仍然存在（因为扫描试图在求解线性化系统后立即计算新的系数矩阵？）。
编辑：我将 Python 2.7.14 与 FiPy 3.2 一起使用。我在下面分享我认为与查询相关的代码部分。整个代码非常广泛。一些背景：我正在求解悬浮流的平衡方程。eq0对应于粒子相的质量平衡方程，phi是粒子的体积分数。
from pylab import *
from fipy import *
from fipy.tools import numerix
from scipy import misc
import osmotic_pressure_functions as opf
kic=96.91
lic=0.049
dt=1.e-2
steps=10
tol=1.e-6
Nx=8
Ny=4
Lx=Nx/Ny
dL=1./Ny
mesh = PeriodicGrid2DTopBottom(nx=Nx, ny=Ny, dx=dL, dy=dL)
x, y = mesh.cellCenters
phi = CellVariable(mesh=mesh, hasOld=True, value=0.,name='Volume fraction')
phi.constrain(0.01, mesh.facesLeft)
phi.constrain(0., mesh.facesRight)
rad=0.1
var1 = DistanceVariable(name='distance to center', mesh=mesh, value=numerix.sqrt((x-Nx*dL/2.)**2+(y-Ny*dL/2.)**2))
pi_ci = CellVariable(mesh=mesh, value=0.,name='Colloid-interface energy map')
pi_ci.setValue(kic*exp(-1.*(var1-rad)/(1.*lic)), where=(var1 > rad))
pi_ci.setValue(kic, where=(var1 <= rad))
def pi_cc_entr(x):
return opf.vantHoff(x)
def pi_cc_vdw(x):
return opf.van_der_waals(x,0.74,0.1)
def pi_cc(x):
return pi_cc_entr(x) + pi_cc_vdw(x)
diffusioncoeff = misc.derivative(pi_cc,phi,dx=1.e-6)
eq0 = TransientTerm() + ConvectionTerm(-pi_ci.faceGrad) == DiffusionTerm(coeff=diffusioncoeff)
step=0
t=0.
for step in range(steps):
print 'Step ', step
phi.updateOld()
res0 = 1e+10
while res0 > tol :
phi.setValue(0., where=phi<0)
res0 = eq0.sweep(var=phi, dt=dt) #ERROR HAPPENS HERE
函数vantHoff和van_der_waals在单独的文件中定义。
def vantHoff(phi):
return phi
def van_der_waals(phi,phi_cp,nd_v):
return (nd_v*phi**3) / ((phi_cp-(phi_cp)**(1./3.)*(phi)**(2./3.))**2)