Crank-Nicolson instability in 3D

[WastedMeat]

I have a perfectly working 2-D finite-difference thermal solver, using an alternating-direction implicit scheme (painfully written in Fortran 90 at my adviser’s insistence), that I have recently extended to 3-D. It is essentially a three step approximation to the Crank-Nicolson equation, which is shown to be unconditionally stable in homogeneous media with constant boundary conditions in the old textbook in which I found it.

Who the hell writes a thermal solver to simulate a single material with constant temperatures at the boundaries?

Anyway, I have adapted the method to an inhomogeneous material structure, and, regardless of boundary conditions, I have exponential error growth (in magnitude...the actual error oscillates), which starts at the material interfaces.

For almost prohibitively low time-steps, though still much larger than what would be required for an explicit method, the thermal solver actually appears stable and produces the correct results as far as I can tell, so I don't think there is any simple error in my implementation. It appears to me that the material inhomogeneities have broken the unconditional stability of the 3-D alternating direction implicit scheme, and I am having a very hard time finding relevant resources.

Posting this here is probably a bit of a fishing expedition, but if this sounds familiar to anyone else, please elaborate on your experiences.

Thanks in advance.

 

[AlephZero]

The first thing to check is if it works in 3D with a homogeneous material.

The "point" of ADI is that you are taking steps in different directiions which are unstable, except (at least in 2D) the instabilities in the two directions cancel out.

I don't know how that would extend either to 3D or inhomogeneous materials, but (leaving aside bugs in the code) I guess that's one place to start researching.

FWIW explicit forward-difference finite difference schemes (with low time steps for accuracy) were used for 3D problems back in the 1970s, but things have moved on a bit since then. For an explicit scheme, you can do the calculations point-by-point, so you don't need anything like "ADI" to avoid solving sets of linear equations.

You can generalize Crank Nicholson into a family of methods with a parameter (say $\theta$) where $\theta = 0$ is forward difference, $\theta = 1/2$ is central difference, and $\theta = 1$ is backward difference. Sometimes the best compromise of accuracy against stability is about $\theta = 0.6$ or $\theta = 2/3$. (It's possible to invent an arm-waving argument that says 2/3 is a "good" value to choose, but I'm not personally very convinced by it.)

 

[WastedMeat]

It works just fine with homogeneous media and simple boundary conditions, but that also would bypass a lot of potential errors in the code, e.g. typos in interpolating the staggered material mesh onto the temperature mesh.

For the system I need to simulate, which contains both glass and silver, and has 20 micron cylindrical substructures, I would need time steps on the order of femtoseconds for an explicit algorithm. It is starting to look like I can get away with 10's of nanoseconds with the ADI, but it certainly isn't unconditionally stable.

I have a fully functioning inhomogeneous θ=1 1-D straight C-N, and a 2-D ADI code that are stable in every scenario in which I have tested them, even coupled to radiative boundary conditions with source inhomogeneities (i.e, well outside the formal proof of unconditional stability), so I am not a complete novice on Crank-Nicolson implementations.

There is just not much literature on real three dimensional ADI implementations. While the source for my algorithm (I can post it and a copy of the relevant paragraph tomorrow) claims it is unconditionally stable for the example given for homogeneous media, I am not familiar enough with numerical analysis to determine easily if the structural inhomogeneities should break the unconditional stability, as they appear to. Then again, maybe nobody is skilled enough to do it "easily."

If in principal ADI's are stable due to canceling errors in each explicit part of the update, perhaps it is not general inhomogeneity that breaks the stability, but anisotropy. Some fully symmetric 3-D test cases do seem to require a less refined time-step than the extrusion style test cases.

 

[AlephZero]

Googling around, the Douglas-Gunn method seems to be the favorite 3D ADI method.

Re your test problems and interpolation etc, some problems where the steady state is a uniform heat flux in different directions, (including heat flux not parallel to the finite difference grid) might flush out some errors.

Most of my heat transfer work has been with FE codes rather than FD, so unfortunately I've forgottem most of what I used to know about FD algorithms

 

[RedBranchKnight]

I have a perfectly working 2-D finite-difference thermal solver, using an alternating-direction implicit scheme (painfully written in Fortran 90 at my adviser’s insistence), that I have recently extended to 3-D. It is essentially a three step approximation to the Crank-Nicolson equation, which is shown to be unconditionally stable in homogeneous media with constant boundary conditions in the old textbook in which I found it.

Who the hell writes a thermal solver to simulate a single material with constant temperatures at the boundaries?

Anyway, I have adapted the method to an inhomogeneous material structure, and, regardless of boundary conditions, I have exponential error growth (in magnitude...the actual error oscillates), which starts at the material interfaces.

For almost prohibitively low time-steps, though still much larger than what would be required for an explicit method, the thermal solver actually appears stable and produces the correct results as far as I can tell, so I don't think there is any simple error in my implementation. It appears to me that the material inhomogeneities have broken the unconditional stability of the 3-D alternating direction implicit scheme, and I am having a very hard time finding relevant resources.

Posting this here is probably a bit of a fishing expedition, but if this sounds familiar to anyone else, please elaborate on your experiences.

Thanks in advance.

Are you using the 3d equivalent of the 2d ADI? If so, then it is unstable. In 2d, the legs are separately unstable but they knock each other out.

The 'fix' is the Peaceman-Rachford method. BTW do you have mixed derivatives like u_xy etc.?

A better approach imo is the Method of Fractional Steps (Marchuk/Janenko). And the Alternating Direction Explicit (ADE) is used as well. Google will help.

 

[vagantes]

I have a perfectly working 2-D finite-difference thermal solver, using an alternating-direction implicit scheme (painfully written in Fortran 90 at my adviser’s insistence), that I have recently extended to 3-D. It is essentially a three step approximation to the Crank-Nicolson equation, which is shown to be unconditionally stable in homogeneous media with constant boundary conditions in the old textbook in which I found it.

Who the hell writes a thermal solver to simulate a single material with constant temperatures at the boundaries?

Anyway, I have adapted the method to an inhomogeneous material structure, and, regardless of boundary conditions, I have exponential error growth (in magnitude...the actual error oscillates), which starts at the material interfaces.

For almost prohibitively low time-steps, though still much larger than what would be required for an explicit method, the thermal solver actually appears stable and produces the correct results as far as I can tell, so I don't think there is any simple error in my implementation. It appears to me that the material inhomogeneities have broken the unconditional stability of the 3-D alternating direction implicit scheme, and I am having a very hard time finding relevant resources.

Posting this here is probably a bit of a fishing expedition, but if this sounds familiar to anyone else, please elaborate on your experiences.

Thanks in advance.

Hello!
1. I suppose that the proof of the unconditional stability of the scheme that you mention is based on the von Neuman analysis... It is important to understand that this analysis does NOT take into account any boundary conditions! It gives the stability criteria only for the FD scheme itself.

2. I know that the generalizations of a 2D ADI scheme to 3D case, known as the Peaceman-Rachford scheme, is NOT unconditionally stable. One must consider the schemes known as Douglas-Rachford or Brian-Douglas - these schemes are unconditionaly stable in 3D ADI.

3. Concerning the boundary conditions in ADI methods: there is a problem with the so-called intermediate boundary conditions - in a correct implementation one must specify boundary conditions at the intermediate time steps. For details see the book by Strikwerda "Finite difference schemes and partial differential equations" - there is a discussion of this problem for the 2D ADI scheme. The same is valid for 3D. The ill-posed boundary conditions can, in principle, produce instabilities. But if you have a non-linear heat equation (e.g., thermal conductivity depends on temperature) the situation with imtermediate boundary conditions becomes much more complicated...

4. An anticipated question: is there a method of stability analysis that could take into account the boundary conditions? I do not know. But I suppose that the method of energy inequalities could solve this problem. Unfortunately I am not familiar with this method and I would be very grateful if someone could explain it to me.
Probably one could apply a Matrix stability analysis.

Please, let me know if you have advanced in solving your problem with inhomogeneties.
Actually, I work on a similar problem - I have to pose 3-rd order boundary conditions in some internal points of the modelling region, i.e. I also have inhomogeneties.

good luck!