[SOLVED] Solving an elliptic PDE using FiPy

Solving an elliptic PDE using FiPy

I'm trying to solve for an elliptic pde using FiPy and I'm running into some convergence problems. The equation I'm trying to solve is:

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

where, L = f(x) and I'm using a tuple of dx values since the solution for alpha is dependent on the mesh.

I've made the following script to solve the equation using 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)

The solution diverges until the value of res is "inf" resulting in the error:

/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

What I have noticed is that the solution converges well when the value of L^2 is constant. However, I've not been able to make it work with varying L. How should I best go about solving this issue?

Any help or guidance is much appreciated.

Thanks in advance.

PS: The data used is available through this link.

Solution

If L^2 is small enough, the solution is unstable even when L^2 is constant.

Changing to an implicit source seems to work, e.g.,

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