Help: I'm stuck on implementation of Solidification (Stefan) problem

Dear Fenics Team,

I am stuck on the Fenics implementation of a solidification (Stefan) problem, maybe you can help.

The PDE and B.C's are described below:

We can rewrite the standard Stefan problem for the heat $T$
in terms of the enthalpy $\Theta$

$\begin{equation}
\frac{\partial\Theta}{\partial t}=\frac{\partial^{2}T}{\partial x^{2}},
\end{equation}$

where $T$ and $\Theta$ are related by,

$T=\ldots$

$
\begin{equation}
\frac{\kappa}{St}\Theta :\quad\Theta<0
\end{equation}
$

$
\begin{equation}
0 :\quad0<\Theta<\frac{1}{K}
\end{equation}
$

$
\begin{equation}
\frac{K}{St}\left(\Theta-\frac{1}{K}\right) :\quad\frac{1}{K}<\Theta
\end{equation}
$

Typical values for the non-dimensional parameters are,
$
[
K=1.8,\quad\kappa=1.8,\quad St=0.75
]$

The boundary conditions are,

$
\begin{equation}
\begin{array}{cc}
\frac{\partial T}{\partial x}=h_{left}\left(1+T\right) & :\quad x=0\
\
\frac{\partial T}{\partial x}=-h_{right}\left(1+T\right) & :\quad x=1
\end{array}\label{eq:bcs}
\end{equation}$

Typical values for the heat transfer coefficients are,
$
[
h_{left}=1.5,\quad h_{right}=0.1
]
$

Thus the weak formulation for the problem is,
$
\begin{equation}
\int_{0}^{1}\frac{\partial\Theta}{\partial t}v+\frac{\partial T}{\partial x}\frac{\partial v}{\partial x}\mathrm{d}x=\left[h_{left}\left(1+T\right)v\right]_{x=0}-\left[h_{right}\left(1+T\right)v\right]_{x=1},
\end{equation}$

where $v$ is the test function.}

Also attached is a python script with my first attempt at solving - It works up until a point when the enthalpy Theta reaches 1/K after which you get a runtime error.

I would ideally like to be able to solve this with an implicit scheme (which makes the problem non-linear but much faster)

Many thanks in advance for your time.

Graham (PhD Student in Applied Mathematics)

from dolfin import *

set_log_level(WARNING)

N = 100
mesh = UnitIntervalMesh(N)
V = FunctionSpace(mesh, "CG", 1)
delta_x = 1.0 / N

dt = Constant(delta_x)

# Define boundaries
class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0) and on_boundary

# Define right boundary
class RightBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0) and on_boundary


# Mark boundaries
left_boundary = LeftBoundary()
right_boundary = RightBoundary()
boundaries = FacetFunction("uint", mesh)
boundaries.set_all(0)
left_boundary.mark(boundaries, 1)
right_boundary.mark(boundaries, 2)

# Redefine boundary measure
ds = ds[boundaries]

Theta_old = Function(V)
Theta = Function(V)
v = TestFunction(V)

K = Constant(1.8)
kappa = Constant(1.8)
St = Constant(0.75)

Tfun_small = kappa / St * Theta
Tfun_mid = Constant(0.0)
Tfun_large = K / St * (Theta - 1.0 / K)

Tfun = min_value(Tfun_small, Tfun_mid) \
       + max_value(Tfun_large, Tfun_mid)

# Robin equations to substitute for inner(grad(T), n)
h_left = Constant(1.9)
leftBC = h_left * (1.0 + Tfun)

h_right = Constant(0.01)
rightBC = -h_right * (1.0 + Tfun)

# Assemble nonlinear functional
F = inner((Theta-Theta_old)/dt, v) * dx \
    + inner(grad(Tfun), grad(v))*dx \
    - (rightBC*v*ds(2) - leftBC*v*ds(1))

t = 0.0
tend = 1.0

# Initial conditions
TempIC = Constant(1.0)
ThetaIC = project((TempIC*St + 1.0) / K, V)

Theta_old.assign(ThetaIC)
Theta.assign(Theta_old)

J = derivative(F, Theta)
problem = NonlinearVariationalProblem(F, Theta, [], J)
solver  = NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['report'] = False

Thetafile = File('Theta.xdmf')
Thetafile << (Theta, t)
while t < tend:
    if t % 0.1 == 0:
        print "t = %.2f" % t
    solver.solve()
    t += float(dt)
    Thetafile << (Theta, t)
    Theta_old.assign(Theta)