The steady-state solution to a diffusion equation in cylindrical geometry using FiPy is rather different from the solution obtained from another software, eg., Mathematica.
The equation is:
$0 = \frac{1}{r}\frac{d}{dr}\left(\frac{r}{T^{1/2}}\frac{dT}{dr}\right) + cte*T^{3/2}$
which means that, by using a CylindricalGrid1D mesh, we can write the equation as:
mesh = CylindricalGrid1D(nr=100, dr=0.01, origin=0.0)
T = CellVariable(name='temperature', mesh=mesh, hasOld=True)
r = mesh.cellCenters()
#BC's
T.constrain(0., mesh.facesRight)
T.faceGrad.constrain(0.,mesh.facesLeft)
#initial temperature profile
T.setValue( 1-r**2)
eq = 0 == DiffusionTerm( coeff=T**(-1/2), var=T) + 20*ImplicitSourceTerm(coeff=T**(1/2), var=T)
viewer = Viewer(vars=T)
eq.solve()
viewer.plot()
raw_input(" Press <enter> to proceed...")
In here I have set cte=20, but the problem remains whatever this value is. I get the solution on the left whereas the solution given by Mathematica is the one on the right:
I then tried to sweep as recommended for such non-linear equation like this. So instead of eq.solve(), I did:
current_residual = 1.0e100
desired_residual = 1e-5
while current_residual > desired_residual:
current_residual = eq.sweep()
T.updateOld()
But I get the error:
/home/antonio/.local/lib/python2.7/site-packages/fipy/solvers/scipy/linearLUSolver.py:66: RuntimeWarning: overflow encountered in square
error0 = numerix.sqrt(numerix.sum((L * x - b)**2))
/home/antonio/.local/lib/python2.7/site-packages/fipy/solvers/scipy/linearLUSolver.py:71: RuntimeWarning: overflow encountered in square
if (numerix.sqrt(numerix.sum(errorVector**2)) / error0) <= self.tolerance:
/home/antonio/.local/lib/python2.7/site-packages/fipy/solvers/scipy/linearLUSolver.py:71: RuntimeWarning: invalid value encountered in double_scalars
if (numerix.sqrt(numerix.sum(errorVector**2)) / error0) <= self.tolerance:
/home/antonio/.local/lib/python2.7/site-packages/fipy/tools/numerix.py:966: RuntimeWarning: overflow encountered in square
return sqrt(add.reduce(arr**2))
/home/antonio/.local/lib/python2.7/site-packages/fipy/solvers/scipy/linearLUSolver.py:58: RuntimeWarning: overflow encountered in multiply
b = b * (1 / maxdiag)
Traceback (most recent call last):
File "stack.py", line 26, in <module>
current_residual = eq.sweep()
File "/home/antonio/.local/lib/python2.7/site-packages/fipy/terms/term.py", line 254, in sweep
solver._solve()
File "/home/antonio/.local/lib/python2.7/site-packages/fipy/solvers/scipy/scipySolver.py", line 61, in _solve
self.var[:] = numerix.reshape(self._solve_(self.matrix, self.var.ravel(), numerix.array(self.RHSvector)), self.var.shape)
File "/home/antonio/.local/lib/python2.7/site-packages/fipy/solvers/scipy/linearLUSolver.py", line 64, in _solve_
permc_spec=3)
File "/usr/lib/python2.7/dist-packages/scipy/sparse/linalg/dsolve/linsolve.py", line 257, in splu
ilu=False, options=_options)
RuntimeError: Factor is exactly singular
Finally, I re-wrote the initial equation in an equivalente form by using the variable V=T^{1/2}. It's easy to see that with V the equation becomes
$0 = \frac{1}{r}\frac{d}{dr}\left(r\frac{dV}{dr}\right) + \frac{cte}{2}V^3$
So I then used the code:
mesh = CylindricalGrid1D(nr=100, dr=0.01, origin=0.0)
V = CellVariable(name='V', mesh=mesh, hasOld = True)
r = mesh.cellCenters()
#BC's
V.constrain(0., mesh.facesRight)
V.faceGrad.constrain(0.,mesh.facesLeft)
#initial V profile
V.setValue( 1-r**2)
eqV = 0 == DiffusionTerm( coeff=1., var=V) + 20*0.5*ImplicitSourceTerm(coeff=V*V, var=V)
T = V*V
viewer = Viewer(vars=T)
eqV.solve()
viewer.plot()
raw_input(" Press <enter> to proceed...")
and the profile obtained is quite similar, but the values in the y-axis are not the same either from the first FiPy solution or Mathematica one! Sweeping is giving the same error as before.
