Finite Element and CFL condition for the heat equation

[pepgma]

I am solving the heat equation in a non comercial C++ finite elements code with explicit euler stepping, and adaptive meshes (coarse in the boundaries and finer in the center). I am aware the CFL condition for the heat equation depends on dt/h**2 for the 1D, 2D, 3D case. When I solve the equation in 2D this principle is followed and I require smaller grids following dt<h**2.

But in 3D the problem seems to be requiring finner and finer grids as I decrease the timestep in what appears to be a dt/h**3 behaviour. Does anyone have an idea what could be happening? is the CFL no longer valid in FEM and 3D? What other factors could be influencing?

 

[LawrenceC]

What partial differential equation are you solving? Fourier or something with a moving fluid where you have first derivatives in space?

 

[AlephZero]

It is possible the analytic solution has a singularity in the space dimension(s) , so a naive attempt to solve it with automatic mesh refinement will "diverge" to an infinite number of elements and an infinitesimal time step.

The simplest fix for that is just to limit the minimum size of elements and/or timesteps to something physcially meaningful for the problem.

If you need small elements in space to capture the solution, an implicit difference scheme in time with variable times steps can be orders of magnitude faster than a conditionally stable implicit scheme.

 

[LawrenceC]

A fully implicit integration is the best approach. With too large a time step, all you lose is accuracy, not stability.

 

[pepgma]

Thanks LawrenceC and AlephZero. Finally, the issue was merely a meshing problem. I am ussing gmsh, which allows to provide a 1d parameter of "typical length". Although this parameter is well conserved in 2d meshes, 3d meshes seem to result with very deformed elements (e.g. tetrahedron with one of the edges near to zero). Optimizing the mesh with the two 'optimize' buttons in Gmsh, helped tough, but still some of the tetrahedrons where very small.

To corroborate, I solved the same equations in a regular 3D grid, where I was sure the element size was conserved, and indeed, the dt<h**2 was followed.

As for your questions, the equation was really just the heat equation div(sigma*grad(phi))=K*dphi/dt

No automatic (i assume, during run time) mesh refinement was used

>The simplest fix for that is just to limit the minimum size of elements and/or timesteps to
>something physcially meaningful for the problem.

I guess this is equivalent to what i found, i was not having real control over the minimum size

>an implicit difference scheme in time with variable times steps can be orders of magnitude faster
>than a conditionally stable implicit scheme.

I understand the first, and I agree, I have moved to an implicit scheme now. Could you explain " conditionally stable implicit scheme" did you mean 'explicit scheme' ?

>A fully implicit integration is the best approach. With too large a time step, all you lose is
>accuracy, not stability.

Thanks!, hadnt see it that way.

 

[LawrenceC]

To advance the first order differential equations in time, you can go from explicit to fully implicit. There is a weight factor from 0 to 1 that controls it. If it is 0, the advancement is explicit. If it is 1, the solution is fully implicit. If it is 1/2, it's called Crank-Nicolson. With the heat equation you wind up with a set of equations written as follows. The order of the matricies is the number of gridpoints in the system.

[B * theta * [K] + [CAP]]{t'} = [[CAP] - (1-B) * theta * [K]]{t''} + {b} * theta

where

[ ] denotes NXN matrix
{} denotes NX1 vector
theta is time step
B is weight 0<B<1 (should be weak inequality - can't type it)
[K] is stiffness matrix consisting of FE approximation for Laplacian operator, NXN
[CAP] is density-specific heat matrix, diagonal, NXN
{b} is a NX1 vector where the boundary conditions wind up
{t'} is the vector of temperatures at the new time time step
{t''} is the vector of temperatures from the previous time step

So you can readily see that if B is anything but zero, you must solve an NXN set of simultaneous equations for each time step. Fortunately, there are some very fast routines (Gauss elimination) that can do this if the matrix is symmetric and the grid numbering is such that bandwidth is minimized. If it is very large, I have used conjugate gradient methods to solve large systems of equations.

If B=0, you have 'explicit' expressions for the answer at the new time step because [CAP] is diagonal. But, as mentioned, it is unstable. If B is greater or equal to 1/2, the solution is stable but it can have damping oscillations if too large a time step is chosen. For B=1, there are no oscillations.