使用 FiPy 求解椭圆 PDE
Solving an elliptic PDE using 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
我正在尝试使用 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