1

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} = (\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.

Raqato
  • 13
  • 3

1 Answers1

0

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
jeguyer
  • 2,379
  • 1
  • 11
  • 15
  • Indeed you are right, a small value of L^2 leads to an unstable solution. Thank you for the answer! – Raqato Apr 01 '21 at 11:38
  • I spent *way* too long trying to prove to myself that it should be unstable, but I'm terrible at both von Neumann stability analysis and analytical solution of ODEs. I should have just tried the `ImplicitSourceTerm` in the first place! – jeguyer Apr 01 '21 at 14:27
  • I admit, in a hurry to solve this issue I didn't spend the time to do the stability analysis. However, if you're curious, the analytical solution to the equation is \[ \alpha(x) = 1.0 - \exp( - x / L) \]. Thanks for developing this amazing tool! – Raqato Apr 02 '21 at 07:42
  • Dirtly little secret: I don't *ever* spend the time to do the stability analysis! FWIW, with Dirichlet conditions at both ends (and for constant L), I got \[ \alpha(x) = 1 - \cosh((2x - x_R)/2L) / \cosh((x_R)/2L) \]. That seemed to agree with the numerical result (even reasonably well for non-uniform L). – jeguyer Apr 02 '21 at 13:02