Hand Coded Heat

Why use numerical packages...

At a Glance

To begin this lesson

• Open the Answers Form
• Go to the directory for the hand-coded heat application

  cd HandsOnLessons/hand_coded_heat
  

A Simple Science Question

Lets say you live in a house with exterior walls made of a single material of thickness, \(L_x\). Inside the walls are some water pipes as pictured below.

You keep the inside temperature of the house always at 70 degrees F. But, there is an overnight storm coming. The temperature is expected to drop to -40 degrees F. Will your pipes freeze before the storm is over?

Governing Equations

In this lesson, we demonstrate the implementation and use of a hand-coded (e.g., does not use any numerical packages) C-language application to model the one dimensional heat conduction equation through a wall as pictured here …

In general, heat conduction is governed by the partial differential (PDE)… \[\frac{\partial u}{\partial t} - \nabla \cdot \alpha \nabla u = 0\]

where u is the temperature within the wall at spatial positions, x, and times, t, \( \alpha \), is the thermal diffusivity of the material(s) comprising the wall. This equation is known as the Diffusion Equation and also the Heat Equation.

Simplifying Assumptions

To make the problem tractable for this short lesson, we make some simplifying assumptions…

1. The thermal diffusivity, \( \alpha \) is constant for all space and time.
2. The only heat source is from the initial and/or boundary conditions.
3. We will deal only with the one dimensional problem in Cartesian coordinates.
4. We will discretize with constant spacing in both space, \(\Delta x\) and time, \(\Delta t\).

In this case, the PDE we need to develop an application to solve simplifies to… \[\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}\]

Discretization

The equation above is a continuous, partial differential equation (PDE) In order to write a computer program to solve this equation, numerically, the first thing we need to consider is how to discretize the equation into a form suitable for numerical computation.

Consider discretizing, independently, the left- and right-hand sides of equation 2. For the left-hand side, we can approximate the first derivative of u with respect to time, t, by the equation… \[\frac{\partial u}{\partial t} \Bigr\vert_{t_{k+1}} \approx \frac{u_i^{k+1}-u_i^k}{\Delta t}\]

We can approximate the right-hand side of equation 2 with the second derivative of u with respect to space, x, by the equation… \[\alpha \frac{\partial^2 u}{\partial x^2}\Bigr\vert_{x_i} \approx \alpha \frac{u_{i-1}^k-2u_i^k+u_{i+1}^k}{\Delta x^2}\]

Setting equations 3 and 4 equal to each other and re-arranging terms, we arrive at the following update scheme for producing the temperatures at the next time, k+1, from temperatures at the current time, k, as \[u_i^{k+1} = ru_{i+1}^k+(1-2r)u_i^k+ru_{i-1}^k\]

where \( r=\alpha\frac{\Delta t}{\Delta x^2} \)

Note that this equation now defines the solution at spatial position i and time k+1 in terms of values of u at time k . This is an explicit numerical method known as the forward in time, centered difference (FTCS) algorithm. As an explicit method, it has some nice properties:

• They are easy to implement.
• They typically require minimal memory.
• They are easy to parallelize.

Exercise #1: Implement the FTCS Algorithm (2 mins)

The function, solution_update_ftcs, is defined below without its body.

static void
solution_update_ftcs(
    int n,              // # of temperature samples in space
    Double *uk1,        // new temperatures @ t = k+1
    Double const *uk0,  // old/last temperatures @ t = k
    Double alpha,       // thermal diffusivity
    Double dx,          // spacing in space, x
    Double dt,          // spacing in time, t
    Double bc_0,        // boundary condition @ x=0
    Double bc_1         // boundary condition @ x=Lx
)
{

    Double const r = alpha * dt / (dx * dx);

    // Sanity check for stability
    if (r > 0.5) return false;

    // Update the solution using FTCS algorithm
    for (int i = 1; i < n-1; i++)
        uk1[i] = r*uk0[i+1] + (1-2*r)*uk0[i] + r*uk0[i-1];

    // Impose boundary conditions for solution indices i==0 and i==n-1
    curr[0  ] = bc_0;
    curr[n-1] = bc_1;

    return true;

}

Open ftcs.C and implement the FTCS numerical algorithm by coding the body of this function.

Exercise #2: Build and Test the Application (1 min)

To compile the code you have just written…

make

A Simple Sanity Check Run

As a sanity check, lets just run the application with no arguments and see what happens…

% ./heat
    runame="heat_results"
    prec="double"
    alpha=0.2
    lenx=1
    dx=0.1
    dt=0.004
    maxt=2
    bc0=0
    bc1=1
    ic="const(1)"
    alg="ftcs"
    savi=0
    save=0
    outi=100
    noout=0
Iteration 0000: last change l2=0.0909091
Iteration 0100: last change l2=2.42918e-06
Iteration 0200: last change l2=4.86446e-07
Iteration 0300: last change l2=1.00929e-07
Iteration 0400: last change l2=2.09483e-08
Iteration 0500: last change l2=4.41684e-09
Counts: Adds:24500, Mults:25001, Divs:1005, Bytes:176

Before running, the application dumps its command-line arguments so the user can see what parameters it was passed to run. In this case, you are seeing the default values. It then runs the problem as defined by the command-line arguments and saves result files to the directory specified by the runame= argument.

% ls -1t | head -n 1
heat_results
% file heat_results/*.*
heat_results/clargs.out:                    ASCII text
heat_results/heat_results_soln_00000.curve: ASCII text
heat_results/heat_results_soln_final.curve: ASCII text

For this simple application, the results are uncomplicated. They are simple ascii text files containing x/y pairs of the computed numerical results.

% cat heat_results/heat_results_soln_final.curve
# Temperature
       0        0
     0.1   0.1039
     0.2   0.2073
     0.3   0.3101
     0.4   0.4119
     0.5   0.5125
     0.6   0.6119
     0.7   0.7101
     0.8   0.8073
     0.9   0.9039
       1        1

Getting Help

Now that we have built the application, at any point, we can get help regarding various options for the heat application like so…

Usage: ./heat <arg>=<value> <arg>=<value>...
    runame="heat"                       name to give run and results dir (char*)
    alpha=0.2           material thermal diffusivity (sq-meters/second) (double)
    lenx=1                                     material length (meters) (double)
    dx=0.1                  x-incriment. Best if lenx/dx==int. (meters) (double)
    dt=0.004                                      t-increment (seconds) (double)
    maxt=2         >0:max sim time (seconds) | <0:min l2 change in soln (double)
    bc0=0                     boundary condition @ x=0: u(0,t) (Kelvin) (double)
    bc1=1               boundary condition @ x=lenx: u(lenx,t) (Kelvin) (double)
    ic="const(1)"               initial condition @ t=0: u(x,0) (Kelvin) (char*)
    alg="ftcs"                            algorithm ftcs|upwind15|crankn (char*)
    savi=0                                   save every i-th solution step (int)
    save=0                              save error in every saved solution (int)
    outi=100                      output progress every i-th solution step (int)
    noout=0                                       disable all file outputs (int)
    prec="double"                       precision half|float|double|quad (char*)
Examples...
    ./heat dx=0.01 dt=0.0002 alg=ftcs
    ./heat dx=0.1 bc0=273 bc1=273 ic="spikes(273,5,373)"

See the note below regarding more information on the ic argument to specify a variety of initial conditions.

When the heat application runs, by default it will store three files in a directory named runame

Testing The heat Application

Before we use our new application to solve our simple science question, how can we assure ourselves that the code we have written is correct or, at the very least, sanity check ourselves that there isn’t anything glaringly incorrect?

We know, maybe even intuitively, that if we maintain constant temperatures at \(A @ x=0\) and \(B @ x=L_x\), then after a long time (e.g. when the solution reaches steady state), we expect it to be a simple linear variation between temperatures A and B. For example, observe what happens after a long time in the one dimensional example below.

% ./heat dx=0.25 maxt=100 ic="rand(125489,100,50)" runame=test

% cat test/test_soln_00000.curve
# Temperature
       0    69.09
    0.25    143.6
     0.5     96.3
    0.75    52.61
       1    131.6
% cat test/test_soln_final.curve
# Temperature
       0        0
    0.25     0.25
     0.5      0.5
    0.75     0.75
       1        1

Exercise #3: Use Applicaton to Model Science Problem of Interest

Lets now use our heat application to model our simple science question.

Additional Information / Assumptions

• Outside temp has been same as inside temp for a long time, 70 degrees F
• Night/Storm will last 15.5 hours @ -40 degrees F
• Walls are 0.25 meters thick wood and pipe is 0.1 meters in diameter
• Pipe will freeze if center point drops below freezing.

Note: An all too common issue in simulation applications is being sure data is input in the correct units. Take care!

Exercise #4: Analyze Results and Do Some Science

Its time to use the results from our simulation to answer the science question of interest. Below we plot results

make plot PTOOL=gnuplot RUNAME=wall

Depending on your situation, the above command may or may not produce a plot looking like below.

Custom Coded Solutions are a Slippery Slope

We all like to write code and build useful tools. However, it is all too easy to see the unfamiliar as an impediment rather than enabler in reaching our sience goals. However, this is a slippery slope. We often start with relatively simple goals and over time wish to evolve our software solutions to ever more challenging science problems.

Examine the lines of code of the complete application here

$ wc -l *.[Ch]
  125 args.C     # User interface
   94 crankn.C   # Alternative solver
   88 Double.H   # Performance tracking
   38 exact.C    # Testing support
   24 ftcs.C     # FTCS solver
  222 heat.C     # Main
   19 heat.H     # Modularization
   27 upwind15.C # Alternative solver
  151 utils.C    # Utilities, I/O, Data Formats
  788 total

Developing generally useful science applications involves many considerations and software engineering challenges.

• More than just one spatial dimension
• More complex geometric shapes and non-cartesian coordinate systems
• Heat sources and radiation
• Laminated, anisotropic or non-linear materials
• Much larger objects involving billions of discretization points and requiring scalability in all phases of the solution.
• Parallelism such as MPI, and/or, GPU and/or many core and/or various parallel runtimes
• Alternative and interoperable discretizations, solvers, time-integrators, optimizers
• An agile and sustainable software design and implementation addressing understandability of code, with encapsulation of complexities, robustness, efficiency, scalability, portability, reproducibility, rigorous testing, etc.

When we employ numerical packages, many of these issues are addressed for us allowing us to focus more of our effort on the software engineering involved in developing applications that address our science questions of interest.

Evening Hands On Session

Short / Quick Follow-on Questions

Determine Optimum Wall Thicknesses

What are the minimum thicknesses of walls of Wood, Adobe and Common brick to prevent the pipes from freezing?

When you are done, go to Intro->Submit A Show Your Work using the hands-on activity name Optimized Walls and upload evidence of your completed solution.

Compare FTCS, Crank-Nicholson and Upwind15 Algorithms (5 points)

Crank-Nicolson Discretization

Using the Crank-Nicolson discretization, we arrive at the following discretization of equation 2… \[-ru_{i+1}^{k+1}+(1+2r)u_i^{k+1}-ru_{i-1}^{k+1} = ru_{i+1}^k+(1-2r)u_i^k+ru_{i-1}^k\]

where \[r= \alpha \frac{\Delta t}{2 \Delta x^2}\]

In equation 7, the solution at spatial position i and time k+1 now depends not only on values of u at time k but also on other values of u at time k+1. This means each time we advance the solution in time we must solve a linear system; in other words we must solve for all of the values at time k+1 in one step. This is an example of an implicit method. In this case, the system of equations is tri-diagonal – since each update for u at i only uses u at i-1 , i and i+1 – so it is easier to implement than a general matrix solve but is still more complicated than an explicit update.

The code to implement this method is more involved because it involves doing a tri-diagonal solve. It is in crankn.C. It involves code that sets up and LU factors the initial matrix. Then, the LU factored matrix is used on each solution timestep to solve for the new temperatures.

Run the same problems using each of these algorithms and observe total memory usage and operation counts (printed at the end) and provide your explanations for them and in comparison with the FTCS method.

When you are done, go to Intro->Submit A Show Your Work using the hands-on activity name Crank-Nicholson and upload evidence of your completed solution.

Use The Application to Solve The Pipeline Problem (5 points)

An pipeline carrying Ethenol-85 (E85) runs between a manure processing facility and a kerosene production factory. In the unlikely event that both facilities experience catastrophic explosion (burning methane at the manure facility and burning kerosene at the kerosene facility), that briefly increases the local air temperature on both sides of the pipe to the burning temperature of the respective materials, determine the minimum thermal diffusivity of the material used to coat/insulate the pipe to prevent the E-85 from exploding. Assume the pipe is 36 inches in diameter.

When you are done, go to Intro->Submit A Show Your Work using the hands-on activity name Pipeline and upload evidence of your completed solution.

Modify the Application to Support Two Materials (10 points)

Using other research, modify the application to work for a composite wall composed of two materials.

When you are done, go to Intro->Submit A Show Your Work using the hands-on activity name Composite Wall and upload evidence of your completed solution.

A note about the ic= argument to heat

The initial condition argument, ic, handles a few interesting cases

Constant, ic="const(V)"

Set initial condition to constant value, V

Ramp, ic="ramp(L,R)"

Set initial condition to a linear ramp having value L @ x=0 and R @ x=\(L_x\).

Step, ic="step(L,Mx,R)"

Set initial condition to a step function having value L for all x<Mx and value R for all x>=Mx.

Random, ic="rand(S,B,A)"

Set initial condition to random values in the range [B-A,B+A] using seed value S.

Sin, ic="sin(Pi*x)"

Set initial condition to \(sin(\pi x)\).

Spikes, ic="spikes(C,A0,X0,A1,X1,...)"

Set initial condition to a constant value, C with any number of spikes where each spike is the pair, Ai specifying the spike amplitude and Xi specifying its position in, x.

Back to all ATPESC 2018 Lessons