Решение эллиптического уравнения в частных производных с помощью FiPy

Я пытаюсь решить эллиптический pde с помощью FiPy, и у меня возникают некоторые проблемы с конвергенцией. Уравнение, которое я пытаюсь решить:

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

где L = f (x), и я использую кортеж значений dx, так как решение для альфа зависит от сетки.

Я сделал следующий сценарий для решения уравнения с помощью 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: используемые данные доступны по этой ссылке.


python pde fipy
person Raqato    schedule 31.03.2021    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Если L ^ 2 достаточно мало, решение нестабильно, даже когда L ^ 2 постоянно.

Кажется, работает переход на неявный источник, например,

eq = DiffusionTerm(coeff=1.0, var=alpha) == ImplicitSourceTerm(coeff=1./coeff_L, var=alpha) - 1.0 / coeff_L
person jeguyer    schedule 01.04.2021
comment
Действительно, вы правы, малое значение L ^ 2 приводит к нестабильному решению. Спасибо за ответ! - person Raqato; 01.04.2021
comment
Я потратил путь слишком долго, пытаясь доказать себе, что он должен быть нестабильным, но я ужасен как в анализе стабильности фон Неймана, так и в аналитическом решении ODE. Я должен был просто попробовать ImplicitSourceTerm в первую очередь! - person jeguyer; 01.04.2021
comment
Признаюсь, в спешке с решением этой проблемы я не потратил времени на анализ стабильности. Однако, если вам интересно, аналитическое решение уравнения - [\ alpha (x) = 1.0 - \ exp (- x / L)]. Спасибо за разработку этого замечательного инструмента! - person Raqato; 02.04.2021
comment
Маленький грязный секрет: я никогда не трачу время на анализ стабильности! FWIW, с условиями Дирихле на обоих концах (и для константы L), я получил [\ alpha (x) = 1 - \ cosh ((2x - x_R) / 2L) / \ cosh ((x_R) / 2L)]. Казалось, что это согласуется с численным результатом (даже достаточно хорошо для неоднородного L). - person jeguyer; 02.04.2021