# Prevent single node spikes in a FEM-simulation (using continuous Galerkin)

I am trying to solve a non-linear time-dependent heat equation $$\partial_tT=\nabla \left(k_T(T)\nabla T\right) + f$$ (similar to question Solving a non-linear heat equation with the galerkin method gives negative values) with a FEM-approach, using the continuous Galerkin method for the FEM part, the Newton method for linearization of the non-linear part and the implicit Euler method for the time stepping.
Initially, I had problems with parts of the results becoming negative, which could be solved (more or less) by clipping all parts which are below a certain value. Now I got the additional problem that the solution at some nodes becomes unstable, and due to a positive feedback loop increases into unrealistic values. This leads to results as shown in the following image:

The heat source is located at the bottom, approximately at 500, which can be seen due to the darker yellow color. Still, in the middle of the image, the solution shows some nodes with values significantly higher values than the surrounding area, leading to a non-smooth result.
Thus I was wondering if there is a method to prevent those runaways, such as a slope limiter (and if yes, how large will the error be)? I can not use a lowpass here, after those values can be reached during regular calculations.
Or is my approach in general flawed for such an equation, and I should switch to another method?

To address the comment: In general I try to solve the equations from https://research.utwente.nl/en/publications/two-temperature-model-for-pulsed-laser-induced-subsurface-modific, with the equation for the electron temperature ($$T_E$$) resulting in node spikes, while the other two equations do not make any problems. The boundary conditions are Neumann conditions with $$dT/dt=0$$. The domain itself is defined as a square from (0, 0) to (1, 1), due to the problem being dimensionless, while the basis functions are polynomials.

• How non-linear is $k_T$? Is it highly oscillatory? Is it piecewise linear?
– Paul
Sep 29, 2018 at 20:11

I believe these issues arise because you are solving the transient form of the heat equation which, being parabolic, can lead to some instability for continuous galerkin method.

There are ways to circumvent that. GGLS approaches (Galerkin Gradient Least Square) introduce stabilization term that dampen these oscillations by minimizing the Error on temperature not only in the L2 space but in H1. You can see this article by Ilinca and Hétu which is very detailed: https://www.sciencedirect.com/science/article/pii/S0045782502002426

There are also alternatives like the source stabilized petrov-galerkin stabilization (SSPG). There is also very good work by Ilinca on this : https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.2324

Disclaimer : F. Ilinca is someone I know.

• The first article has a very strange numerical example. Clearly the initial condition of $T(x)=1$, for all $x \in (0,1)$, is not consistent with the Dirichlet boundary condition at $x=0$, i.e. $T(0)=0$.
– knl
Sep 29, 2018 at 7:52
• Absolutely, but this kind of BC exacerbates this type of instability so it's a good show case I guess.
– BlaB
Sep 29, 2018 at 18:35
• Could this behaviour also be related to instabilities of the Newton method for linearization of the non-linear part? Oct 22, 2018 at 14:26