How to use the Thomas-Algorithm to the Heat-diffusion-equation correctly

Question

My post is structured in four parts:

1. I give you some information about the context my principal questions refer to.

2. I will tell you what I believe to know about the Thomas Algorithm. If I am wrong in any regards, correct me please.

3. I suggest two approaches on how to implement time stepping for the heating process.

4. I will state my questions

Part I: Background Information

I use the finite-difference-method to simulate Joule-heating.

1. The first PDE is the Laplace-equation. As a solution it delivers the electrical potential field. This potential field is used to calculate the joule-heating, which is the source term in the heat-diffusion equation. Thus, for subseqent analysis we can regard the potential field/joule heating as given.

2. The second PDE is the heat-diffusion equation for an instationary/time-dependent heating. I use the Peaceman-Rachford scheme to split the heat-diffusion-equation in two half-time steps. Below you can see the heat-diffusion equation casted in the time-splitting scheme of Peaceman-Rachford:

This term $\mu_x=\frac{\lambda dt}{\rho c_p \Delta x^2}$ is called Fourier-number, where $\lambda$ is the thermal conductivity, $\rho$ is the densitiy, $c_p$ is the specific heat, $dt$ is the time step, $\Delta x$ is the spacial step in x-direction. $\delta_x^2$ is the differential operator for the second derivative. Subsequently, $\delta_x^2T_{ij}=T_{i+1,j}-2T_{i,j}+T_{i-1,j}$. The same applies analogously to $\mu_y$. $F_{i,j}$ is the source term representing the joule heating.

3. I apply the Thomas algorithm to each of both equations of the Peaceman-Rachford-scheme to solve for the temperature field over a geometry in the x-y-plane.

4. Based on the Joule-heating, material parameters, an initial value for the temperature field and the boundary conditions the algorithm is supposed to compute the instationary/transient temperature field. Of course, it also should be capable of computing the heating up to the stationary state.

Part II: What do I believe to know about the Thomas Algorithm

1. It is a direct solver for linear equation systems with tridiagonal coefficient matrices. Direct means that - depending on the number of linear equations of the system - you can solve the linear equation systems with a fixed number of basis operations or rather you don't need to perform iterations over the solution vector. In case of the Thomas Algorithm the number of basic operations is of Order 0(N4), where N stand for the number of equations.
2. It is an exact solver, you can obtain exact solution provided that the personal computer does not have any round off errors.
3. You should be able to get a result after running your Thomas-Algorithm over your linear equation systems only ONCE, i.e. you don't need to iterate your solution over and over again and therefore you don't need a convergence criterion.

Part III: time stepping

1. approach

After finishing the first time step, I take the solved temperature field as the input for the second time step. After finishing the second time step, I take the solved temperature field as the input for the third time step and so forth for all subsequent time steps. The pseudo-code for this method can be seen below

---

PSEUDO-CODE START

number_iteration = 100;
dt = 0.01 s;

for i = 1:number_iteration
do{

Use temperature field of time step n as input for first equation of PR-Scheme, i.e. apply following equation:

$(1-\frac{\lambda dt}{\rho c_p \Delta x^2}\delta_x^2)T_{i,j}^{n+0.5}=(1+\frac{\lambda dt}{\rho c_p \Delta y^2}\delta_y^2)T_{i,j}^{n}+\frac{\ dt}{2}F_{i,j}^{n+0.5}$

Apply Thomas Algorithm to first equation of PR-Scheme (directly above) to solve for temperature field of time step n+0.5.

Use temperature field of time step n+0.5 as input for second equation of PR-Scheme, i.e. apply subsequent equation:

$(1-\frac{\lambda dt}{\rho c_p \Delta y^2}\delta_y^2)T_{i,j}^{n+1}=(1+\frac{\lambda dt}{\rho c_p \Delta x^2}\delta_x^2)T_{i,j}^{n+0.5}+\frac{\ dt}{2}F_{i,j}^{n+0.5}$

Apply Thomas Algorithm to second equation of PR-Scheme (directly above) to solve for temperature field of time step n+1.

}

---

PSEUDO-CODE END

For each time step dt = 0.01 s is used. The overall simulation time amounts to 1 s (100 runs x 0.01 s).

2. approach

How does the approach 2 differs from approach 1? Approach 2 accumulates the time steps by dt = dt + dt_inc. Approach 1 always uses a constant time step dt. This is the only difference. Subsequently, in approach 2 an increased dt is used for calculating the temperature field for the next time step. The pseudo-code for this approach is stated below:

---

PSEUDO-CODE START

number_iteration = 100
dt = 0 s;
dt_inc = 0.01 s;

for i=1:number_iteration
do {
dt = dt + dt_inc

Use temperature field of time step n as input for first equation of PR-Scheme, i.e. apply following equation:

$(1-\frac{\lambda dt}{\rho c_p \Delta x^2}\delta_x^2)T_{i,j}^{n+0.5}=(1+\frac{\lambda dt}{\rho c_p \Delta y^2}\delta_y^2)T_{i,j}^{n}+\frac{\ dt}{2}F_{i,j}^{n+0.5}$

Apply Thomas Algorithm to first equation of PR-Scheme (directly above) to solve for temperature field of time step n+0.5.

Use temperature field of time step n+0.5 as input for second equation of PR-Scheme, i.e. apply subsequent equation:

$(1-\frac{\lambda dt}{\rho c_p \Delta y^2}\delta_y^2)T_{i,j}^{n+1}=(1+\frac{\lambda dt}{\rho c_p \Delta x^2}\delta_x^2)T_{i,j}^{n+0.5}+\frac{\ dt}{2}F_{i,j}^{n+0.5}$

Apply Thomas Algorithm to second equation of PR-Scheme (directly above) to solve for temperature field of time step n+1.

}

---

PSEUDO-CODE END

Part III: Questions

Which approach is the right one?

Answer 1

For a linear PDE, like the Laplace equation, when you discretize it you should get a linear system. Since you're 1D, the Thomas algorithm should be able to solve the system, and it's executed by running over the system once; Thomas algorithm is a direct solver, not an iterative one. If I understand your question correctly, you're asking what happens if you run Thomas algorithm multiple times over the same linear system. It should always return the same result (to within machine error) as it is not iterative. If you're trying to run a time stepping problem, you need to break it up into time steps and solve the linear system at each time step. If you're running the Thomas algorithm always at the same condition, it will not timestep, and you'll always get the same result, one step in the future. This is the best answer I can give you given the information that you gave. If you explained what a PC is and showed some math as to your actual PDEs and the splitting that would be appreciated. Also, LES is typically an abbreviation for Large Eddy Simulation, not linear equation system.

EDIT:: I asked for more information to be able to answer the question and it was provided.

The Peaceman-Rachford is your timestepping equation. You begin with the condition at step n then you solve a linear system with thomas algorithm to get the solution at $n + \frac{1}{2}$. Then you take state $n + \frac{1}{2}$ and plug that state into the right hand side of the second Peaceman-Rachford equation to get the solution at step n. You repeat this until you reach the desired time of the simulation (so sum your $\Delta t$ at each step until $t = t_f$). both approaches that you put up should work if you proper implement them, but your 2nd approach pseudo code has some pretty bad errors, like not being consistent as to whether dt is change in time or temperature...

Answer 2

With PC I mean personal computer. Sorry for not mentioning this in advance. LES is a typical abbreviation for Large Eddy Simulation? I did not know that before. Ok, I will memorize it and regard it in future posts.

Here comes the heat-diffusion equation casted in the time-spitting scheme of Peaceman-Rachford:

$\mu_x=\frac{\lambda dt}{\rho c_p \Delta x^2}$, where $\lambda$ is the thermal conductivity, $\rho$ is the densitiy, $c_p$ is the specific heat, $dt$ is the time step, $\Delta x$ is the spacial step. $\delta_x^2$ is the differential operator for the second derivative. Subsequently, $\delta_x^2T_{ij}=T_{i+1,j}-2T_{i,j}+T_{i-1,j}$. The same applies analogously to $\mu_y$. $F_{i,j}$ is the source term representing the joule heating. I apply the Thomas-Algorithm to both equations successively.

What do I do when I apply the Thomas Algorithm onto my Heat-Diffusion equation iteratively? After finishing the first run, I take the solved temperature field as the input for the second run. After finishing the second run, I take the solved temperature field as the input for the third run and so forth for all subsequent runs. Exactly, I just enforce the Thomas-Algorithm to behave like an iterative solver. In this case, the Thomas-Algorithm doesn’t start with the same conditions from run to the next run because the temperature field as the input changes from run to run. Consequently, I don’t get the same result for all runs. What confuses me, is, that I get a trustworthy result even if I enforce the Thomas-Algorithm to behave like an iterative solver. This described method is described as an pseudocode below by approach 1.

You are writing: “If you’re trying to run a time stepping problem, you need to break it up into time steps and solve the linear system at each time step.” Yes, I am trying to run a time stepping problem. This leads me to the next question. How should I implement the time stepping itself?

I have two suggestions for how to implement time stepping:

approach 1

number_iteration = 100;
dt = 0.01 s; 

for i = 1:number_iteration
do{
- use temperature field of last time step as input for next time step
- solve both heat-diffusion-equations (see figure at top of this post) by means of Thomas-Algorithm
-> result: temperature field of current time step 
}

For each (loop) run the time step dt = 0.01 s is used. Does this approach integrate time by means of the for-loop? If so, the overall simulation time would amount to 1 s (100 runs x 0.01 s). This also means that the time step dt in the above mentioned equations (the figure at the top this post) would always remain the same from run to run. I am not sure about the correctness of this approach 1. Is that right?

The second approach follows:

approach 2

simulation_time = 1 s;
dt = 0;
dt_inc = 0.010 s;

while (simulation_time > dt)
do {
dt = dt + dt_inc
- use temperature field of last time step as input for next time step
- solve both heat-diffusion-equation (see figure at top of this post) by means of Thomas-Algorithm
-> result: temperature field of current time step
}

The approach 2 calculates the temperature field by dt = dt + dt_inc, i. e. the temperature field for the time step dt is larger by dt_inc compared to the temperature field of the previous time step. This differs from the approach 1, where the time step always equals dt. The resulting temperature field (dt) in approach 2 serves as the input for calculating the temperature field of the next time step (dt = dt + dt_inc).

The question is: which approach is correct, the first or second approach? If none of these are right, what is the right implementation for the time stepping problem? I am not sure, because both approaches lead to results that look right.

You also write: “If you’re running the Thomas algorithm always at the same condition, it will not timestep, and you’ll always get the same result, one step in the future.” I don’t understand it, especially the last part “….it will not timestep, and you’ll always get the same result, one step in the future.” Is this what my suggested approach 1 of the time stepping implementation does? Could you explain me your statement with different words again, please?