使用 FiPy 求解椭圆偏微分方程

问题描述

我正在尝试使用 FiPy 求解椭圆 pde，但遇到了一些收敛问题。我试图解决的方程是：

\[ \frac{\partial^2 \alpha}{\partial x^2} = (\alpha -1)/L^2 \]

其中，L = f(x) 并且我使用的是 dx 值的元组，因为 alpha 的解决方案取决于网格。

我编写了以下脚本来使用 FiPy 求解方程：

from fipy import *
import numpy as np

deltax = tuple(np.genfromtxt('delx_eps.dat')[:,0])

mesh = Grid1D(dx=deltax,nx=257)

# compute L^2
CL = 0.161
Ceta = 80.0
eps = np.genfromtxt('delx_eps.dat')[:,1]
nu = np.empty_like(eps)
nu[:] = 1.0/590.0
L_sq = np.empty_like(eps)
L_sq = (CL*Ceta*(nu**3/eps)**(1/4))**2

coeff_L = CellVariable(mesh=mesh,value=L_sq,name='Lsquare')

alpha = CellVariable(mesh=mesh,name='Solution',value=0.0)
# Boundary conditions
alpha.constrain(0.,where=mesh.facesLeft)
alpha.constrain(0.,where=mesh.facesRight)

eq = DiffusionTerm(coeff=1.0,var=alpha) == (alpha-1.0)/coeff_L

mySolver = LinearLUSolver(iterations=10,tolerance=5e-6)
res = 1e+100
while res > 1e-8:
    res = eq.sweep(var=alpha,solver=mySolver)
    print(res)

解决方案发散直到 res 的值是 "inf" 导致错误：

/usr/local/lib/python3.8/dist-packages/fipy/variables/variable.py:1143: RuntimeWarning: overflow encountered in true_divide
  return self._BinaryOperatorVariable(lambda a,b: a / b,other)
/usr/local/lib/python3.8/dist-packages/fipy/variables/variable.py:1122: RuntimeWarning: invalid value encountered in multiply
  return self._BinaryOperatorVariable(lambda a,b: a*b,other)
Traceback (most recent call last):
  File "elliptic_shielding.py",line 73,in <module>
    res = eq.sweep(var=alpha,solver=mySolver)
  File "/usr/local/lib/python3.8/dist-packages/fipy/terms/term.py",line 232,in sweep
    solver._solve()
  File "/usr/local/lib/python3.8/dist-packages/fipy/solvers/scipy/scipySolver.py",line 26,in _solve
    self.var[:] = numerix.reshape(self._solve_(self.matrix,self.var.ravel(),numerix.array(self.RHsvector)),self.var.shape)
  File "/usr/local/lib/python3.8/dist-packages/fipy/solvers/scipy/linearLUSolver.py",line 31,in _solve_
    LU = splu(L.matrix.asformat("csc"),diag_pivot_thresh=1.,File "/usr/local/lib/python3.8/dist-packages/scipy/sparse/linalg/dsolve/linsolve.py",line 337,in splu
    return _superlu.gstrf(N,A.nnz,A.data,A.indices,A.indptr,RuntimeError: Factor is exactly singular

我注意到的是，当 L^2 的值是常数时，解会很好地收敛。但是，我无法使其与不同的 L 一起使用。 我应该如何最好地解决这个问题？

非常感谢任何帮助或指导。

提前致谢。

PS：使用的数据可用 through this link。

解决方法

如果 L^2 足够小，即使 L^2 是常数，解也是不稳定的。

更改为隐式来源似乎有效，例如

eq = DiffusionTerm(coeff=1.0,var=alpha) == ImplicitSourceTerm(coeff=1./coeff_L,var=alpha) - 1.0 / coeff_L