heat_equation_1d

• Introduction
• The Model
• Implementation
• Operation
• Stability
• Usage

Introduction

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.

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:

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 );

Implementation

The 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.

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

element.t4:

static
    int init
    element_type_enum element_type

sub
    change my_change

pub
    double y 
    double x

• static:
• init is a flag that indicates whether the system has just started up (discussed below)
• element_type controls whether the element is at the START, MIDDLE or END of the bar, and therefore what to do w.r.t. the ghost elements.
• sub:
• my_change is the change object to which the element subscribes
• pub:
• x is the x-value of the element. It is derived from the element's id in the calculate function. In the default distribution these are set from 0.1 to 1.0
• y is the y-value for the element. If the system has just been started this will be cos(2*x). If the system has run through once this will be y = y + change.

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.

change.t4:

sub
    data_server my_data_server
    element up_element
    element this_element
    element down_element

pub
    double change

• sub:
• my_data_server is the data_server to which every change subscribes
• up_element is the element to the right of this_element along the bar
• this_element it the element to which this change applies
• down_element is the element to the left of this_element along the bar
• pub:
• change is the change that should be applied to this_element

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.

data_server.t4:

static
    double k
    double alpha
    double beta

• static:
• k it the thermal conductivity of the bar
• alpha is the given Dirichlet boundary value for the left end of the bar
• beta is the given Dirichlet boundary value for the right end of the bar

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.

monitor.t4:

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

static
    int counter
    int print_each_cycle
    double converged_limit
    double diverged_limit

pub
    double start_time

• sub:
• my_elements is an array containing the ids of all the element objects
• my_changes is an array containing the ids of all the change objects
• my_data_server is there so that the monitor can report on data_server's data (e.g. k)
• static:
• counter allows for the monitor to keep track of how many cycles run before convergence/divergence
• print_each_cycle controls wether the monitor prints out the system each cycle
• converged_limit is the threshold every change_change value has to be less than for the system to be considered converged. In worked_example1 this is set to 1.0e-06
• diverged_limit is the threshold any change_change value has to be greater than for the system to be considered diverged. In worked_example1 this is set to 1.0
• pub:
• start_time allows the monitor to time the system

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.

Operation

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.

Start-up

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:

load_all_data_servers()
load_all_changes()
load_all_monitors()
        load_monitor_my_changes()
        load_monitor_my_elements()
load_all_elements()

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:

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:

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

Next iteration

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:

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

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

Stability

Convergence

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.

Divergence

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.

Usage

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