5
$\begingroup$

I am using backward Euler in a FEM scheme for a convection-diffusion problem. On a given mesh, I can take arbitrarily large time steps, as expected. But if I decrease time step, at some point it will generate oscillations in the solution (spikes). Is this a known behavior? It's not mentioned in text books, at least the ones I know. Would Crank–Nicolson scheme eliminate the issue?

The above mentioned problem is only the simplest one to demonstrate the phenomenon. I attach an image of what is happening in my real problem: transient (here only 2D) incompressible flow. $\rho=1$, dynamic $\eta=1\times 10^{-3}$, flat inflow 1, cylinder diameter 0.1 ($R_e=100$). If timestep is 0.0001, such oscillations occur. If timestep is 0.005, solution becomes smooth and I reproduce von Karman vortices with proper frequency, so it is highly unlikely I have a bug in the code.

This is standard Galerkin FEM with no stabilization solved with a direct solver.

Any other thoughts how this is possible and how to know how small timestep is "too small".

enter image description here

Thanks for any hints.

Dominik

Improve this question
$\endgroup$
7
  • 1
    $\begingroup$ Are you keeping the spatial grid spacing constant as you change the time step? $\endgroup$
    – David Ketcheson
    Dec 31, 2011 at 14:04
  • $\begingroup$ Also, are you solving a steady-state or time-dependent problem? $\endgroup$
    – David Ketcheson
    Jan 1, 2012 at 8:37
  • $\begingroup$ Yes, I keep the spatial grid size constant and I am doing a time dependent problem. When I solve for steady state, the oscillations in the solutions are gone. $\endgroup$
    – Dominik Lark
    Jan 2, 2012 at 8:43
  • $\begingroup$ I have re-edited the original question to reflect my progress. $\endgroup$
    – Dominik Lark
    Jan 7, 2012 at 12:31
  • 1
    $\begingroup$ This is normal. Write the equations down for 2 linear elements. If you initial solution is 1 everywhere and you then impose a value of 2 on the first node and 1 on the last (3rd) node, you will see that there exists a small timestep for which you can have oscillations on node 2. $\endgroup$
    – mfnx
    Mar 11, 2020 at 11:32

2 Answers 2

Reset to default
3
$\begingroup$

Assuming that you keep the spatial grid (mesh) size constant and decrease the time step, one thing that could be happening is that as you decrease the time step, you move into the small, unstable region of the backward Euler method in the right half-plane of the stability diagram. (See Ascher and Petzold, page 52.) The region of instability for backward Euler is $\{z \in \mathbb{C}: |1 - z| \leq 1\}$, where $z = h \lambda$, $h$ being your step size and $\lambda$ being the eigenvalue of the test equation $\dot{y}(t) = \lambda y(t)$.

Since the true solution operator for convection-diffusion only has eigenvalues with negative real part, if your semi-discretization yields eigenvalues with positive real part then it's a very bad semi-discretization!

Improve this answer
$\endgroup$
6
  • 3
    $\begingroup$ That could only be true if the semi-discretization has positive eigenvalues, which wouldn't make much sense for advection-diffusion. $\endgroup$
    – David Ketcheson
    Jan 1, 2012 at 8:34
  • 1
    $\begingroup$ Without a source term, I don't see why you would have positive eigenvalues either, which is why I said that I doubted that it was the problem. Assuming the operators are discretized correctly, I don't see why decreasing the time step should yield oscillations, but the oscillations could be a symptom of an issue not directly related to time-stepping (say, a bug in the implementation of the discretization of the operators, yielding an ODE system whose right-hand side has a Jacobian matrix with positive eigenvalues). $\endgroup$
    – Geoff Oxberry
    Jan 1, 2012 at 20:37
  • 1
    $\begingroup$ Also, the implicit solve may not be accurate enough. Depending on the solution algorithm, it's easy for the small time step limit to become poorly scaled. $\endgroup$
    – Jed Brown
    Jan 2, 2012 at 0:35
  • $\begingroup$ @GeoffOxberry I edited the answer to reflect our discussion. Hope you don't mind. $\endgroup$
    – David Ketcheson
    Jan 2, 2012 at 13:23
  • $\begingroup$ @ Jed: nice hint, but I checked, and no, the scaling seems OK, gets actually better (since smaller dt increases max value in A due to 1/dt term). I also checked with better accuracy of the implicit solve, the problem persists. @Geoff: thanks for pointing out stability condition for the backward Euler, I always believed it is unconditionally stable. But also, I can't see how I could be getting positive eigenvalues... using standard Galerkin FEM on second-order tets here. $\endgroup$
    – Dominik Lark
    Jan 2, 2012 at 22:38
0
$\begingroup$

Did you do convergence tests of your code?

From an absolute stability point of view (think necessary but not sufficient conditions of stability for nonlinear PDEs) it should be safe to reduce the time step while maintaining the spatial discretization fixed in a method of lines type of discretization. Thus, it could be seen as surprising that the stability appear to worsen when time step is refined.

How is the backward euler method implemented and solved? Except for linear PDEs where direct solvers can be used immediately, it is typical to use either a fix-point iteration scheme or a Newton-Raphson method for solving the implicit scheme to advance the solution one time step. In both cases it should be favorable to decrease the time step.

How did you impose boundary conditions?

A more detailed description of your numerical scheme would be helpful to understand exactly what you have done.

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ Thanks for the answer. What do you mean by convergence tests? The inner linear solves are direct, there is no convergence criterion. As for the outer solves, yes, I tightened the tolerance and the problem persists. As for the scheme, I am really doing standard stuff here: Bubnov-Galerkin 2nd order momentum 1st order pressure, fully implicit time integration, fixed point iteration. Boundary conditions: fixed flat inflow profile, gradient condition on outflow, no-slip elsewhere. The matrices and rhs's are well scaled. I am really surprised to find this effect and am determined to find an answer. $\endgroup$
    – Dominik Lark
    Jan 8, 2012 at 21:42
  • $\begingroup$ @Dominik: Could you explain where fixed point iteration is used? $\endgroup$
    – faleichik
    Jan 30, 2012 at 15:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.