Skip to content

betterscientificsoftware/hello-heat-equation

Repository files navigation

The Heat Equation is the "Hello World" of Scientific Computing

"Hello World!" is a prototypical program that is universally used by programmers to illustrate the basic syntax and implementation difficulty of a programming language (see several Hello World implementations at WikiBooks). In scientific computing, software can be written in any number on languages and must use various programming models or math libraries to achieve high performance for any particular problem. The Heat Equation, a well known partial differential equation (PDE), is the perfect formula to demonstrate a "Hello World!" comparison for programming options.

What is the Heat Equation?

The Heat Equation, also known as the Diffusion Equation, was developed by Joseph Fourier in 1822. Primarily used in the fields of Mathematics and Physics, the equation shows how heat (or something similar) changes over time in a solid medium. While there are many adaptations for different dimensions and scenerios, new programmers are often first introduced to the 1-dimensional version of the Heat Equation as a tutorial problem.

The following description is adapted from a Trinity University lecture.

Suppose you have a thin metal rod along the x-axis starting at position 0 and ending at position L. The goal is to model the temperature distribution throughout the rod at a given time, with the following (ideal) assumptions:

  • An initial temperature distribution is given
  • Boundary conditions (established tempuratures at each end of the rod) are known
  • Heat only moves horizontally through the bar
  • There are no external heating/cooling sources

To model the temperature distribution of the rod, we can use the 1-dimensional Heat Equation given by the formula:

∂u / ∂t = α∇ ²u

Where:

  • u = Temperature at a position, x
  • t = time
  • α = thermal diffusity of the material
  • = The LaPlace operator

Assumptions and Discretization

To make the problem more approachable, we make some simplifying assumptions:

  • α is constant
  • The heat/cooling source only comes from initial/boundary conditions
  • Only Cartesian coordinates in 1 dimension are used

Since the heat equation takes the form of a PDE over a continuous space, most programs create a discretization of the formula over a finite set of sampling points.

First, we approximate the derivative ∂u / ∂t, to be

∂u / ∂t | tk+1 (uk+1i - uki) / Δt

Then, we approximate α∇ ²u to be

α∇ ²u | xi α (uki-1 - 2 uki + uki+1) / Δx²

By setting these equations equal to each other, we can find the temperature at time k+1 at spot i using the discretized formula:

uk+1i = ruki+1 + (1 - 2 r)uki + ruki-1

(More information about the full discretization process is available.)

Where r = α Δt² / Δ. This discretization is referred to as the Forward-Time, Central-Space (FTCS) scheme which is considered numerically stable as long as r is less than or equal to 0.5.

Algorithm

After the continuous formula is discretized into a programmable problem we can define our variables and calculate each iteration of the temperature distribution algorithm. Before each stage, we first need to check that r = α * ∂t / (∂x²) is greater than 0.5 to verify it is numerically stable.

Once r has been verified, we can cycle through the non-boundary locations of the mesh and update it according to the discretization uk+1i = ruki+1 + (1 - 2 r)uki + ruki-1, ensuring that the boundary conditions stay constant at each time iteration.

After this process has repeated some specified number of times, we can return the final distribution of our mesh.

Why the Heat Equation?

PDEs are widely used among the different sciences to model various phenomena such as electromagnetism, fluid mechanics, and heat transfers. While PDEs extend across multiple fields, they can still be quite tricky to work with because unlike ordinary differential equations, they are comprised of derivatives of more than just a single variable. A single example of a PDE is the Heat Equation, which is used calculate the distribution of heat on a region over time. While here we just focus on the 1-dimensional version of the Heat Equation, it can actually take a multitude of forms including the Fourier, LaPlace (also known as steady-state), 2D, and 3D heat equations.

As a consequence of the Heat Equation being so expansive and well-known, it is commonly used as an introductory example to new users of scientific software. This problem is also very easy for users to visualize, which can help them understand how different algorithms approach the solution. It allows users to become familiar with the structure and syntax of a particular programming language, parallel paradigm, or mathematical library. Because this tutorial exercise is so common, several implementations exist that are able to be referenced with just a quick web search.

Implementations

Similar to "Hello World!" programs, various implementations of the heat equation can demonstrate syntactic features of a programming language. However, a sequential implementation demonstrates little more than for-loop syntax. Scientific software demands parallel programming models and libraries. Thus, we can demonstrate some interesting implementations.

Variable Definition

Throughout these implementations, we use the following parameters:

n     = Number of temperature samples        # Integer
uk1   = New temperatures across x-axis       # Array of Doubles
uk0   = Old temperatures across x-axis       # Array of Doubles
alpha = Thermal Diffusity                    # Double
dx    = Spacing in space                     # Double
dt    = Spacing in time                      # Double
bc0   = Beginning boundary condition (x = 0) # Double
bc1   = End boundary condition (x = L_x)     # Double

Programming Model Example: OpenMP (C)

OpenMP is a portable way to add parallelism to a program. The OpenMP Heat Equation demonstrates the use of #pragma syntax. We see that this implementation is a 2-line change from the base implementation.

Library Example: MFEM (C++)

Scientific software programmers are able to take advantage of numerical libraries instead of hand-coding implementations, which can optimize performance and allow access to the additional features. Most of these implementations have multiple steps that are outlined to produce accurate outputs and visual depictations of the solution.

MFEM is an opensource C++ library for finite element methods. They provide a serial implementation of the Heat Equation, though they also produce a parallelized version of this code. This is simplified even further, in a (very) condensed version of the code, to demonstrate the core Heat Equation algorithm.

This implementation clearly demonstrates the complexity introduced by using a third-party library. However, users gain access to a number of numerical solvers, parallelism, and even GPU capability.

Library Example 2: FD1D_HEAT_EXPLICIT (C)

While highly developed numerical libraries such as MFEM have a wide range of functionality and support options, others are smaller. They still provide support for different use cases of function like the heat equation, but do not come with the obscurity and extra features that may intimidate users.

FD1D_HEAT_EXPLICIT is a simple 1-dimensional Heat Equation solver which uses the explicit version of the finite difference method to handle time integration. Versions of this library are written in C, C++, Fortran90, Matlab, and Python. A code example from the library demonstrates a sequential implementation, but can easily be adapted to become parallelized.

While this code is almost identical the hand-coded example above, this library does include some extra features that allows users to implement some additional functionality including writing results to different kinds of R8 files, adding a timestamp, and solving for the Courant-Friedrichs-Loewy coefficient.

This repository is a showcase of 1D Heat Equation implementations.

This file shows sequential implementations of a function that calculates a single iteration of the 1-dimensional heat equation using the forward in time, centered difference (FTCS) discretization written in several programming languages. These languages include C/C++, Fortran, Julia, Matlab, and Python. While some differences in syntax can be observed, the overall structure of the code stays consistent through the implementations.

These examples were created by @Heatherms27

This file currently only shows a parallel implementation of a function that calculates a single iteration of the 1-dimensional heat equation using the forward in time, centered difference (FTCS) discretization written using OpenMP with C as its base. This file will be expanded.

This example was created by @Heatherms27

While a user may find it more comfortable to hand code their own implementation of the heat equation, it would be more beneficial to take advantage of pre-existing codes. Several libraries have already provided a code base that saves users time and allows them to access additional features the library may provide. Below are a couple of library implementations of the heat equation.
For a full list of finite element software packages, see this list from Wikipedia.

As a result of the Heat Equation being so expansive and taking on several different forms, many examples of heat equation implementations are readily availible through a quick Google search. This file contains a table that includes links to user implementations of the 1-dimensional, 2-dimensional, and LaPlace forms of the Heat Equation in C, C++, Fortran, Julia, Matlab, and Python.

About

Various implementations of the FTCS form of the Heat Equation

Resources

Stars

Watchers

Forks

Releases

No releases published

Packages