Application version 3.0
Application licence BSD
Trad4 version 3.0
Document version 3.0
Author schevans
Date 23-06-09


The heat equation is part of an important class of problems known as Partial Differential Equations (PDEs). The trad4 implementation uses Finite Difference Methods (FDMs) which along with the trad4 approach means that the system can evolve concurrently.

This approach will scale for 2D, 3D or nD models of the heat, wave and convection-diffusion equations, all of which are coming soon.

While the 1D heat equation has been discussed extensively elsewhere, I will briefly summarise the problem here.

Take a fully insulated metal bar with one dimension (think of a lagged copper wire). Provide an initial heat distribution across the bar, and hold each end at a specific temperature (i.e., using Dirichlet boundary conditions). Now let this bar settle into a steady state, taking snapshots of the temperature distribution across the bar as it evolves over time. We should see the temperature converge on a linear slope between the start and end temperatures.

There is a parameter k that gives the thermal conductivity of the bar, i.e. the rate at which heat dissipates through the bar.

The Model

As discussed, we'll be using FDMs for this model.

Firstly, we need to create our grid. For this example (worked_example1) we'll be using a very rough grid of 10 elements. Each element will have an x value which won't change through the life of the object, and an y value which will start with our initial conditions and then evolve into equilibrium.

Let's look at the starting values first. Each object has an id, an x value and a y value, where y = cos(2*x), and where y = 0.0 if y < 0.0.

id 1 2 3 4 5 6 7 8 9 10
x 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
y 0.980067 0.921061 0.825336 0.696707 0.540302 0.362358 0.169967 0.000000 0.000000 0.000000

Now we need to compute the next cycle of the system's evolution, and for this well be using the Central Difference Method, the derivation of which won't be given here.

We have the state at t=n. At t=n+1, the change in any given element is given by:

my_change_change = my_data_server_k * ( up_element_y - 2*this_element_y + down_element_y );          (1)

So that y at t+1 is given by:

element_y = element_y + my_change_change;          (2)

As discussed above, this first model uses the Dirichlet boundary conditions. This just means that the ends of the bar are held at a specific temperature. The way we visualise this is by imagining an extra ghost element at the edges of our bar, like so:

id ghost1 1 2 3 4 5 6 7 8 9 10 ghost2
x n/a 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 n/a
y 0.980067 0.980067 0.921061 0.825336 0.696707 0.540302 0.362358 0.169967 0.000000 0.000000 0.000000 0.000000

In this model the edges are given somewhat arbitrary values of alpha=0.980067 and beta=0.0 respectively.

It's important to stress that these ghost nodes don't actually exists - they just allow for the equation above to work at the edges. In this way the equation above becomes, for element id=1:

change_change = my_data_server_k * ( up_element_y - 2*this_element_y + my_data_server_alpha );


The abstract diagram

Heat 1D Abstract Diagram

The change object calculates the change that is applied to the object that subscribes to that change object. It uses the y-value of the object that subscribes to it, plus the y-values of it's neighbouring elements. The change object also subscribes to the singleton data_server object.

The concrete diagram

Only 5 nodes are shown below when there are ten nodes in worked_example1. The data_server is not shown on this to keep things simple, but remember every change object also subscribes to a single instance of the data_server.

Heat 1D Concrete Diagram

As you can see from the concrete diagram, a change object subscribes to three element objects so that it is able to calculate equation (1) above. This change gets fed back in to the element object. As such, this is the first application we've seen that has a feedback loop.

The t4 files


    int init
    element_type_enum element_type

    change my_change

    double y 
    double x

The element object holds the current value of the temperature at a particular x-point on the bar. The element object has an init flag in the static section that is set to 0 in the DB. When the object first fires, if init has not been set it will set the initial x and y values, in this case x = id/10 and y = cos(2*x). If init has been set that means the system has run through once and there are some valid change_change values to apply.


    data_server my_data_server
    element up_element
    element this_element
    element down_element

    double change

The change subscribes to three element objects - the one below and the two either side of the graph. It also subscribes to the data_server object.


    double k
    double alpha
    double beta

The data_server is a singleton object, meaning there is just one instance of the type. All change objects subscribe to this single instance. The advantage of this approach is that you don't get the duplication of data where every change object would carry a copy of alpha in it's static section and the resultant space and maintains issues. The disadvantage of this approach is that every change object has to keep track of and read from this data_server object which could have performance implications.


    change my_changes[NUM_NODES]
    element my_elements[NUM_NODES]
    data_server my_data_server

    int counter
    int print_each_cycle
    double converged_limit
    double diverged_limit

    double start_time

The monitor has two jobs in this system. The first is to subscribe to all the change objects and after each change tier (T2) is run, check to see if the system has converged, diverged or still running. The system is considered to have converged if all the change_change values are below converged_limit, and diverged if the change_change values are above diverged_limit.

The monitor is unusual in a few respects. Firstly, it is a singleton like data_server, but unlike data_server which is a T1 object along with all the elements, the monitor object is the only T3 object. This is because it needs to check all the published change_change values to test for convergence. As such it could be viewed as a concurrency bottleneck, but it need not be as the option exists to either a) create a bunch of T3 monitors all subscribing to a sub-set of the T2 objects with a T4 uber_monitor aggregating those monitors, or b) split the functionality of the monitor into two - one to check for convergence and one to check for divergence.

Secondly, it it the first object we've seen that uses the new (in trad4v2.0.1) sub vector mechanism which allows for an object to subscribe to many or all objects of a particular type. The list of change object to which the monitor object subscribes is given in the DB table 'monitor_my_changes'.

Lastly, the monitor object has the job of printing out the data in a consumable format. This is somewhat of a work-around for the lack of support for any feed out of trad4, but it works ok, so expect to see this again in other applications. There is a consequence of this that isn't covered in the above abstract and concrete diagrams and that is in order to print the element's data it must subscribe to them too. An indeed it does - as given in the DB table 'monitor_my_elements'. These links are not shown in the above two diagrams mainly because they make a hell of a mess of the diagrams, but also, arguably, because the act of printing the system is not relevent to it's function.


We'll now step through a couple of iterations, examining the log file as we go. The runtime logs are on-line in full -. there is one with print_each_cycle=0, one with print_each_cycle=1 and a full set of data in csv format.


As usual on start-up, the types are created and their shared objects loaded:

Creating new type data_server, type id: 1, tier: 1
Creating new type change, type id: 3, tier: 2
Creating new type monitor, type id: 4, tier: 3
Creating new type element, type id: 2, tier: 1

Then, the data is loaded:


Next, the system is checked and it's constituents reported:

Checking tier 1. Num objects this tier: 11
Checking tier 2. Num objects this tier: 10
Checking tier 3. Num objects this tier: 1
Checking tier 4. Num objects this tier: 0
Checking tier 5. Num objects this tier: 0

Next, the objects are validated:

Validating objects...

Then the threads start:

Starting thread 1
Starting thread 2
Starting thread 3
Starting thread 4

Then tier1 runs. At this stage, as it's the first time the system has been run, the init flag will not have been set so the element just calculate their y-values using the equation defined in their calculate method - in this case y = cos(2*x) and where y = 0.0 if y < 0.0. We've seen this before above:

id 1 2 3 4 5 6 7 8 9 10
x 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
y(t) 0.980067 0.921061 0.825336 0.696707 0.540302 0.362358 0.169967 0.000000 0.000000 0.000000

Then tier 2 runs. At this stage the changes are calculated. This is done using equation (1) above, but bear in mind the ghost elements caused by the boundary conditions we've selected. Running this tier will give:

id 1 2 3 4 5 6 7 8 9 10
x 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
y(t) 0.980067 0.921061 0.825336 0.696707 0.540302 0.362358 0.169967 0.000000 0.000000 0.000000
change -0.0195361 -0.0183599 -0.0164518 -0.0138877 -0.0107701 -0.00722303 0.0112117 0.0849836 0 0

And lastly the tier-3 monitor runs. This checks for convergence and prints out the data if configured to do so.

Next iteration

Due to the presence of the feedback loop, the system needs no external influence to continue to run. When the T1 element's tier is run, they will see that the objects they subscribe to - the T2 objects - have changed, and that they therefore need to fire.

The next thing they will see is that their init flag's been set, and that therefore they will not need to seed themselves, but rather to update their published values given the delta they pull up from the change objects to which they subscribe. This means that equation (2) above runs, which results in:

id 1 2 3 4 5 6 7 8 9 10
x 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
y(t) 0.980067 0.921061 0.825336 0.696707 0.540302 0.362358 0.169967 0 0 0
change(t) -0.0195361 -0.0183599 -0.0164518 -0.0138877 -0.0107701 -0.00722303 0.0112117 0.0849836 0 0
y(t+1) 0.96053 0.902701 0.808884 0.682819 0.529532 0.355135 0.181179 0.0849836 0 0

The change object then fires again as discussed above, which produces:

id 1 2 3 4 5 6 7 8 9 10
x 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
y(t) 0.980067 0.921061 0.825336 0.696707 0.540302 0.362358 0.169967 0 0 0
change -0.0195361 -0.0183599 -0.0164518 -0.0138877 -0.0107701 -0.00722303 0.0112117 0.0849836 0 0
y(t+1) 0.96053 0.902701 0.808884 0.682819 0.529532 0.355135 0.181179 0.0849836 0 0
change(t+1) -0.00917995 -0.0179939 -0.0161238 -0.0136109 -0.0105554 0.00022083 0.0388803 0.00560587 0.0424918 0

These changes continue to ripple round the system until the system reaches equilibrium (more on this below).



In the worked_example1 default data set, convergence is considered to have occurred if the value of every change_change is less than 1.0e-06. This is controlable at run-time through the monitor static, and the value you'd set would be a trade-off between your floating point precision and how much time you have on your hands.


In the worked_example1 default data set, divergence is considered to have occurred if the value of any change_change is greater than 1. This is a somewhat arbitrary value and could be set lower or, better, derived from the step_size, but in practise divergence is detected quite quickly.

Divergence is caused by the curve being steeper than the step_size can support, or the rate at which the heat dissipates though the bar - k - again being too fast for the granularity of this grid.


Database changes

Database values like k, alpha and beta are simple to change through a sql-prompt:

$ t4db
SQL> update data_server set k=0.3;

You can update the monitor static in the same way:

SQL> update monitor set print_each_cycle=1;

Functional changes

To change the shape of the initial heat distribution you will have to modify the calculate_element function (found in heat_equation_1d/objects/element.c). There you will be able to change both how the x and y values are derived from the object's id. So if you wanted x = 3-30 you would set x as:

element_x = id * 3.0;

And if you wanted y = sin(2*x) you would set y as:

element_y = sin( 2 * element_x );

Extending the grid

There is a 100-node grid checked in to worked_example2 but due to a limitation in the current trad4 infrastructure you will need to increase the value of NUM_NODES in heat_equation_1d/src/constants.t4s, run the precompiler and recompile:

heat_equation_1d$ t4p
heat_equation_1d$ make all

Get trad4 at SourceForge.net. Fast, secure and Free Open Source software downloads