Why am I getting NaNs modeling PDE system?

General
1

I’m modeling the following PDE system:

T_t = \nabla \cdot(k\nabla T)-\vec{v}\cdot\nabla T+A\left(Se^{-B/(T-T_a)}-C(T-T_a)\right) \\ S_t = -C_SSe^{-B/(T-T_a)}, \;\;T>T_a

I’ve derived the variational form as follows:

We apply a forward Euler finite difference scheme for the time derivative. Then we multiply by a test function, q, and integrate over our domain to get:

\begin{align*} \int_{\Omega} \frac{(T_{n+1}-T_n)}{dt}q\; d\vec{x} &= \int_{\Omega}\bigg(k\Delta T_nq-\vec{v}\cdot\nabla T_nq + A\big(S_ne^{-B/(T_n-T_a)}q-C(T_n-T_a)q\big)\bigg)\; d\vec{x} \\ &= \int_{\partial\Omega}k\nabla T_n\cdot\hat{\textbf{n}}q\; d\vec{s} + \int_{\Omega}\bigg(-k\nabla T_n\cdot\nabla q-\vec{v}\cdot\nabla T_nq + A\big(S_ne^{-B/(T_n-T_a)}q-C(T_n-T_a)q\big)\bigg)\; d\vec{x} \end{align*}

by Green’s first identity where \hat{\textbf{n}} is the outward pointing unit normal at each point on the boundary.

We apply pure Neumann boundary conditions, so the integral over the boundary is zero. Then, we have

\begin{align*} \int_{\Omega} (T_{n+1}-T_n)q + dt\bigg(k\nabla T_n\cdot\nabla q + \vec{v}\cdot\nabla T_nq - A\big(S_ne^{-B/(T_n-T_a)}q-C(T_n-T_a)q\big)\bigg)\; d\vec{x} &= 0 \end{align*}

For our fuel equation (S), we have

\begin{align*} \int_{\Omega} \frac{(S_{n+1}-S_n)}{dt}q\; d\vec{x} &= \int_{\Omega}-C_SS_ne^{-B/(T_n-T_a)}q\; d\vec{x} \end{align*}

which simplifies to

\begin{align*} \int_{\Omega} (S_{n+1}-S_n)q + dt\big(C_SS_ne^{-B/(T_n-T_a)}q\big)\; d\vec{x} &= 0 \end{align*}

I’ve elected to solve for T (temperature) at each step, then use the new values for T to solve for S (fuel) at each step, and use those new values for T at the next time step.
My fenics code is as follows:

# Define function space, trial and test functions.
mesh = fa.RectangleMesh(fenics.Point(x1, y1), fenics.Point(x2, y2), Nx, Ny)
V = fenics.FunctionSpace(mesh, 'CG', 1)
T = fenics.TrialFunction(V)
S = fenics.TrialFunction(V)
q = fenics.TestFunction(V)
p = fenics.TestFunction(V)

# Define initial values.
T_0 = fenics.Expression('1200/cosh(pow(x[0]-b,2)+pow(x[1]-b,2))', degree=2, b=4)
T_n = fa.interpolate(T_0, V)
S_0 = fenics.Expression('1/cosh(pow(x[0]-b,2)+pow(x[1]-b,2))', degree=2, b=5)
S_n = fa.interpolate(S_0, V)

# Define Boundary Conditions: BC = dot(grad(T_n),N) where N is the unit normal vector.  
BC = fa.Expression("0", degree=1) # Pure Neumann BCs.

# Define variational forms.
F_T = (T-T_n)*q*dx \
        + dt*(k*dot(grad(T_n),grad(q)) + dot(v,grad(T_n))*q - A*(S_n*exp(-B/(T_n-Ta))*q - C*(T_n-Ta)*q))*dx \
        + (k*BC*q)*ds # ds = integral over boundary.
F_S = (S-S_n)*p*dx + dt*Cs*S_n*exp(-B/(T_n-Ta))*dx

# Separate into bilinear and linear form.
a_T, L_T = lhs(F_T), rhs(F_T)
a_S, L_S = lhs(F_S), rhs(F_S)

# Step through time.
Cs = [T_n.compute_vertex_values(mesh)]
Ds = [S_n.compute_vertex_values(mesh)]
T = fa.Function(V)
S = fa.Function(V)
for _ in t:
    # Solve for T, update solution.
    fa.solve(a_T==L_T, T)
    T_n.assign(T)
    # Solve for S, update solution
    fa.solve(a_S==L_S, S)
    S_n.assign(S)
    
    # Save.
    C = T.compute_vertex_values(mesh)
    Cs.append(C)
    D = S.compute_vertex_values(mesh)
    Ds.append(D)

The first iteration seems to execute okay, but on step 2 I get NaNs for T. I know there are a variety of reasons this could be, but I’m new to FEniCS and I’m not sure if I’ve set up the code correctly and there’s an issue of stability? Or if I’m missing some important piece of the code itself? Or in my variational form? I’d appreciate any help.

2

I think it’s just because the diffusion isn’t stable with the finite difference scheme I selected. I’m doing a backward scheme with a linearization via Taylor Approximation and it seems to be working now.