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.
I'm not convinced this problem has any solution other than T = 0. Further, that solution seems to be unstable against different values of the initial condition and/or of cte. That instability is not altogether surprising given that the equation in T form would have an infinite diffusivity when T = 0.
I suspect that Mathematica is doing roughly what FiPy is doing in your first set of figures, which is to show the first sweep of this non-linear problem. It's not the answer; just the first guess. I don't know anything about solving PDEs with Mma, though, either analytically or numerically.
The plot after one sweep after your V solution looks different, btw, because you didn't adjust the initial condition. It should be:
V.setValue( numerix.sqrt(1-r**2))