力变量在Fipy中为非负

Question

我目前正在使用一个sweep循环来解决一个微分方程（eq0），这个微分方程是关于我在python中使用FiPy的单元格变量phi。因为我的方程是非线性的，所以我使用了一个sweep循环，如下面代码的摘录所示。你知道吗

while  res0 > resphi_tol:    
        res0 = eq0.sweep(var=phi, dt=dt)

但我不断得到以下错误：

C:\Python27\lib\site-packages\fipy\variables\variable.py:1100: RuntimeWarning: invalid value encountered in power return self._BinaryOperatorVariable(lambda a,b: pow(a,b), other, value1mattersForUnit=True)
C:\Python27\lib\site-packages\fipy\variables\variable.py:1186: RuntimeWarning: invalid value encountered in less_equal return self._BinaryOperatorVariable(lambda a,b: a<=b, other)
Traceback (most recent call last):
.. File "SBM_sphere3.py", line 59, in
....res0 = eq0.sweep(var=phi, dt=dt)
.. File "C:\Python27\lib\site-packages\fipy\terms\term.py", line 207, in sweep
....solver._solve()
.. File "C:\Python27\lib\site-packages\fipy\solvers\pysparse\pysparseSolver.py", line 68, in _solve
....self.solve(self.matrix, array, self.RHSvector)
.. File "C:\Python27\lib\site-packages\fipy\solvers\pysparse\linearLUSolver.py", line 53, in _solve__
....LU = superlu.factorize(L.matrix.to_csr())
.. File "C:\Python27\lib\site-packages\pysparse\misc__init__.py", line 29, in newFunc
....return func(*args, **kwargs)
.. File "C:\Python27\lib\site-packages\pysparse__init__.py", line 47, in factorize
....return self.factorizeFnc(*args, **kwargs)
RuntimeError: Factor is exactly singular

我很确定这个错误是由于等式0中的φ（2/3）项引起的。如果我用abs（phi）^（2/3）替换这个项，错误就消失了。你知道吗

我假设sweep循环在某个点为phi中的几个单元格返回负值，这会导致错误，因为我们不能用非整数指数对负值进行幂运算。你知道吗

所以我的问题是：有没有一种方法可以强制扫描以避免负面的解决方案？你知道吗

我尝试在扫掠之前包含一条将所有负值设置为0的线：

while  res0 > resphi_tol:    
        phi.setValue(0.,where=phi<0.)
        res0 = eq0.sweep(var=phi, dt=dt)

误差仍然存在（因为sweep试图在解线性化系统之后计算新的系数矩阵？）。你知道吗

编辑：我正在使用python2.7.14和fipy3.2。我在下面分享我认为与查询相关的代码部分。整个代码相当扩展。 背景：我在解悬浮流的平衡方程。等式0对应于颗粒相的质量平衡方程，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)

Answer 1

产生错误的原因是DiffusionTerm的系数都是nan。这又是因为扩散系数定义为

(((((((Volume fraction + -1e-06) + (((pow((Volume fraction + -1e-06), 3)) * 0.1) / (pow((0.74 - ((pow((Volume fraction + -1e-06), 0.6666666666666666)) * 0.9045041696510275)), 2)))) * -0.5) + 0.0) + (((Volume fraction + 0.0) + (((pow((Volume fraction + 0.0), 3)) * 0.1) / (pow((0.74 - ((pow((Volume fraction + 0.0), 0.6666666666666666)) * 0.9045041696510275)), 2)))) * 0.0)) + (((Volume fraction + 1e-06) + (((pow((Volume fraction + 1e-06), 3)) * 0.1) / (pow((0.74 - ((pow((Volume fraction + 1e-06), 0.6666666666666666)) * 0.9045041696510275)), 2)))) * 0.5)) / 1e-06)

而且Volume fraction（phi）都是零，所以-1e-06被提升为分数次幂，这是未定义的。你知道吗

-1e-06的因素来自于你使用scipy.misc.derivative()来计算符号导数？我不相信这是故意的。使用SymPy可能会有更好的运气。你知道吗