4
$\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

$\endgroup$
  • 1
    $\begingroup$ Are you keeping the spatial grid spacing constant as you change the time step? $\endgroup$ – David Ketcheson Dec 31 '11 at 14:04
  • $\begingroup$ Also, are you solving a steady-state or time-dependent problem? $\endgroup$ – David Ketcheson Jan 1 '12 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 '12 at 8:43
  • $\begingroup$ I have re-edited the original question to reflect my progress. $\endgroup$ – Dominik Lark Jan 7 '12 at 12:31
  • $\begingroup$ It's always helpful to write down precisely the equations you want to solve and describe the discretization carefully. I take it we are talking about incompressible Navier-Stokes? How do you enforce incompressibility? Through a Poisson solve? What is eta? $\endgroup$ – David Ketcheson Jan 7 '12 at 18:06
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!

$\endgroup$
  • 2
    $\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 '12 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 '12 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 '12 at 0:35
  • $\begingroup$ @GeoffOxberry I edited the answer to reflect our discussion. Hope you don't mind. $\endgroup$ – David Ketcheson Jan 2 '12 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 '12 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.

$\endgroup$
  • $\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 '12 at 21:42
  • $\begingroup$ @Dominik: Could you explain where fixed point iteration is used? $\endgroup$ – faleichik Jan 30 '12 at 15:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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