現在、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 return self._BinaryOperatorVariable(lambda a,b: a<=b, other)
トレースバック (最新call last):
.. ファイル "SBM_sphere3.py", 行 59, in
....res0 = eq0.sweep(var=phi, dt=dt)
.. ファイル "C:\Python27\lib\site-packages\ fipy\terms\term.py"、207 行、sweep
....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 行目、factorize
....return self.factorizeFnc(*args, **kwargs)
RuntimeError: Factor は厳密に特異です
このエラーは、 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)
エラーはまだ残っています (線形化されたシステムを解いた直後に、スイープが係数の新しい行列を計算しようとするためですか?)。
編集： FiPy 3.2でPython 2.7.14を使用しています。クエリに関連すると思われるコードの部分を以下に共有しています。コード全体は非常に拡張されています。いくつかのコンテキスト: サスペンション フローのバランス方程式を解いています。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)