Python differential analysis of heat loss across a pipe

I made a differential equation solver to forecast heat loss through a pipe with heat loss coefficient contours. It does what I want it to do but I had to use global variables to change the outputs of my simulation function - is there a better way to deal with this issue? Another problem I have is that my code is very slow - when I increase dx the answer seems to become less accurate (the graphs seem to approach a single result past 0.0001)

import scipy
# R_w = 0.028# W/mk air heat loss ocefficient
R_w = 0.654# W/mk water heat loss ocefficient
R_f = 0.042 # W/mk pipe heat loss coefficient
r = 0.2 # m pipe radius
T_0 = 50 # C initial temperature
T_a = 20 # C ambiant temperature
dx = 0.0001 # m differential distance
L = 25 # m length of pipe

def heat_balance(T_1,T_2):
heat_in = (T_1-T_2) * scipy.pi * r**2 * R_w
heat_out = ( (T_1 + T_2) / 2 - T_a ) * 2 * r * scipy.pi * dx * R_f
return abs(heat_in-heat_out)

def simulation():
T = [T_0]
X = [0]
print 'r is ',R_f
import numpy as np
from scipy.optimize import minimize
for x in np.arange(dx,L,dx):
X.append(x)
fun = lambda s: heat_balance(T[-1],s)
T.append(float(minimize(fun,T[-1]-dx,method="SLSQP")['x']))
return([T,X])

from matplotlib import pyplot as plt
typs = ['+','x','|','_','*']
#typs = ['D','s','8','p','^']
cols = ['red','blue','black','purple','brown']

for i, rf in enumerate([7,14,50,100,200]):
R_f = 1.0/rf
ans = simulation()
plt.scatter(ans[1],ans[0],c=cols[i],marker=typs[i],s = 100,label=rf,lw=3)

plt.legend(fontsize=20,title='Pipe Insolation R Value')
plt.xlabel('Distance (m)',fontsize=20)
plt.ylabel('Average Temperature (C)',fontsize=20)
plt.show()


Globals in Python are not necessarily a bad thing, especially for what you're using them for. However, they can be prevented. There are a couple of major flaws I found so I'm glad you came over to get a review. The point you're most worried about will be handled at the end.

Style

Clean Python code adheres to PEP8, the general style guide for Python. Your code violates this guide in a couple of ways. I'll explain the violations I deem most important and fixed it up for you to show the difference. I do recommend you read PEP8 yourself as well.

To verify whether your code adheres to the guide, you can either use the pep8 tool inside Python itself or pep8online. This will validate your code against most of the guide.

Import

Import statements should be on top of your file for maintainability purposes. This is a common practice in many languages, including C, C++, C# and Java. If I now skim your code, I'll see unexpected functions like plt.xlabel() because I didn't notice the following line:

from matplotlib import pyplot as plt


All imports should be on the top.

Globals

Imports should be followed by your global variables. Often those globals are supposed to be constant (note that Python does not actually enforce them being constant). The naming convention for constants is UPPERCASE_WITH_UNDERSCORES. This makes it clear to anyone reading the code the value should not be changed and is probably defined at module level.

This also means that other variables should not use the same casing. All relevant casing styles are in the PEP8.

Variable naming

Now we're handling variables anyway, please use more descriptive variable names. Single-letter variables and functions make for hard-to-read code. Code should be self descriptive. With good names, many of your comments would become redundant. Functions with names are fun can definitely use something more descriptive.

Whitespace

Give your operators and commas some breathing room. foo, bar(woof + 1), fizz is a lot easier on the eyes than foo,bar(woof+1),fizz.

Indentation is done with spaces, not with tabs. Multiples of 4 are preferred, but multiples of 2 are also common.

import scipy
import numpy as np
from scipy.optimize import minimize
from matplotlib import pyplot as plt

R_w = 0.654  # W/mk water heat loss ocefficient
R_f = 0.042  # W/mk pipe heat loss coefficient
PIPE_LENGTH = 25  # m length of pipe
TEMPERATURE_INITIAL = 50  # C initial temperature
TEMPERATURE_AMBIANT = 20  # C ambiant temperature
dx = 0.0001  # m differential distance

def heat_balance(T_1, T_2):
heat_in = (T_1-T_2) * scipy.pi * PIPE_RADIUS**2 * R_w
heat_out = ( (T_1 + T_2) / 2 - TEMPERATURE_AMBIANT ) * 2 * PIPE_RADIUS * scipy.pi * dx * R_f
return abs(heat_in-heat_out)

def simulation():
T = [TEMPERATURE_INITIAL]
X = [0]
print 'r is ', R_f
for x in np.arange(dx, PIPE_LENGTH, dx):
X.append(x)
fun = lambda s: heat_balance(T[-1], s)
T.append(float(minimize(fun, T[-1] - dx, method="SLSQP")['x']))
return([T, X])

typs = ['+', 'x', '|', '_', '*']
cols = ['red', 'blue', 'black', 'purple', 'brown']

for i, rf in enumerate([7, 14, 50, 100, 200]):
R_f = 1.0/rf
ans = simulation()
plt.scatter(ans[1], ans[0], c=cols[i], marker=typs[i], s = 100, label=rf, lw=3)

plt.legend(fontsize=20, title='Pipe Insolation R Value')
plt.xlabel('Distance (m)', fontsize=20)
plt.ylabel('Average Temperature (C)', fontsize=20)
plt.show()


Objective Oriented Programming

This is one of the best solutions against global variables: objective oriented programming.

Did you know you can make classes in Python?

class PipeCalculator:
self.R_w = 0.654  # W/mk water heat loss ocefficient
self.R_f = 0.042  # W/mk pipe heat loss coefficient
self.PIPE_LENGTH = 25  # m length of pipe
self.TEMPERATURE_INITIAL = 50  # C initial temperature
self.TEMPERATURE_AMBIANT = 20  # C ambiant temperature
self.dx = 0.0001  # m differential distance

def heat_balance(self, T_1, T_2):
heat_in = (T_1-T_2) * scipy.pi * self.PIPE_RADIUS**2 * R_w
heat_out = ( (T_1 + T_2) / 2 - self.TEMPERATURE_AMBIANT ) * 2 * self.PIPE_RADIUS * scipy.pi * dx * R_f
return abs(heat_in-heat_out)

def simulation(self):
T = [self.TEMPERATURE_INITIAL]
X = [0]
print 'r is ', self.R_f
for x in np.arange(self.dx, self.PIPE_LENGTH, self.dx):
X.append(x)
fun = lambda s: self.heat_balance(T[-1], s)
T.append(float(minimize(fun, T[-1] - self.dx, method="SLSQP")['x']))
return([T, X])

def plot_calculations(self):
typs = ['+', 'x', '|', '_', '*']
cols = ['red', 'blue', 'black', 'purple', 'brown']

for i, rf in enumerate([7, 14, 50, 100, 200]):
self.R_f = 1.0/rf
ans = simulation()
plt.scatter(ans[1], ans[0], c=cols[i], marker=typs[i], s = 100, label=rf, lw=3)

plt.legend(fontsize=20, title='Pipe Insolation R Value')
plt.xlabel('Distance (m)', fontsize=20)
plt.ylabel('Average Temperature (C)', fontsize=20)
plt.show()

calculator = PipeCalculator()
calculator.plot_calculations()


And poof, all the globals are gone.

Of-course, a straight port from what you have to object oriented still looks ugly. Pass some variables around, use the power of arguments and when in doubt don't hesitate to split up to more functions. But at least you now have a start and it will be easier to re-use this code later on in different projects.

As two final notes:

• Your lambda could be re-written to a normal function for improved readability. This may help you figure out how to optimize the calculations.
• Pipes usually have an inside and outside diameter. You may want to explicitly state which one you're using.
• @ferada Yes, there are many more additions to be made here. However, there is such a thing as overkill and I didn't want to overdo it. If OP would post a second iteration in a new question later, I'd definitely remark such. – Mast Nov 13 '15 at 12:13
• This is a superficial use of a class. It's just moving variables and functions from the module namespace to the class namespace. – hpaulj Nov 13 '15 at 17:22
• @hpaulj Isn't that exactly what he wanted, to move variables away from the module namespace? – Mast Nov 13 '15 at 17:43

Accuracy

I don't really see much of an accuracy problem by increasing $\Delta x$. Compare your original:

… to the result using my code below, with $\Delta x=0.01$ (and some clarifying visual tweaks):

I also tried my code using $\Delta x=0.001$, and the results are visually indistinguishable from the bottom plot, except that the dots in the initial temperature drop are less sparse. (If you are concerned about connectedness, then just change your plot configuration.)

If anything, it's your choice of scatter plot markers that hurts the precision of the output.

The heat-transfer equations should, after all, lead to stable results. If the temperature at point $x$ is estimated to be a bit higher than the true value, then the temperature drop over the interval between $x$ and $x + \Delta x$ would increase as well.

Increasing $\Delta x$ by a factor of 100 naturally cuts the computation by a factor of 100.

Preamble

For ease of maintenance, imports should go at the top, unless you have a good reason for hiding them.

from collections import namedtuple
from matplotlib import pyplot as plt
from numpy import arange
import scipy
from scipy.optimize import minimize
import functools


The global variables would be tidier if you placed them in a dictionary. The tradeoff is that the functions that use them become more verbose. (Watch your spelling in the comments.)

model = {
'R_w': 0.654,       # W/mk water heat loss coefficient
'R_f': None,        # W/mk pipe heat loss coefficient
'r':   0.2,         # m pipe radius
'T_0': 50,          # C initial temperature
'T_a': 20,          # C ambient temperature
'L':   25,          # m length of pipe
}


Instead of defining separate lists for typs, cols, and rf values, it would be clearer to unify them. I would also move these definitions to the top of the program.

Series = namedtuple('Series', 'rf marker color')
R_values = [
Series(  7, '.', 'red'),
Series( 14, '.', 'blue'),
Series( 50, '.', 'black'),
Series(100, '.', 'purple'),
Series(200, '.', 'brown'),
]


Optimization

The heat_balance function has a lot of multiplications that could be precomputed for each scenario.

def make_heat_balance(model, dx):
heat_in_mult = scipy.pi * model['r']**2 * model['R_w']
heat_out_mult = 2 * model['r'] * scipy.pi * dx * model['R_f']

def heat_balance(T_1, T_2):
heat_in = heat_in_mult * (T_1 - T_2)
heat_out = heat_out_mult * ((T_1 + T_2) / 2 - model['T_a'])
return abs(heat_in - heat_out)

return heat_balance


The simulation() function should take parameters. I've made dx a parameter to simulation() because it's not really part of the model.

Mixing temperature and distance units in T[-1]-dx feels unclean to me, even if it works. It's also unnecessary, since it's merely an initial guess for the minimization procedure.

I prefer to define X all at once.

The lambda can be replaced with functools.partial().

If you just tweak the print call, you can make the whole program compatible with Python 3. Calling it "r" instead of "R_f" is confusing — it looks like you are varying the pipe radius.

def simulation(model, dx=0.01):
print('R_f is {0:f}'.format(model['R_f']))
X = arange(0, model['L'], dx).tolist()
T = [model['T_0']]
heat_balance = make_heat_balance(model, dx)
for x in X[1:]:
heat_balance_T1 = functools.partial(heat_balance, T[-1])
T.append(minimize(heat_balance_T1, x0=T[-1], method="SLSQP")['x'])
return X, T


Plotting

In simulation(), I've switched the order of the return values, which prettifies the call to plt.scatter().

"Insolation" should be spelled "Insulation". I'm not sure what you mean by "Average Temperature" — isn't it just "Temperature"?

def main():
for series in R_values:
plt.scatter(*simulation(dict(model, R_f=1.0 / series.rf)),
label=series.rf,
c=series.color, marker=series.marker, edgecolors='face',
)

plt.legend(fontsize=20, title='Pipe Insulation R Value')
plt.xlabel('Distance (m)', fontsize=20)
plt.ylabel('Temperature (C)', fontsize=20)
plt.show()

if __name__ == '__main__':
main()

• In my timing tests, a make_heat_balance kind of precalculation (I used a class), shaves about 30% off a simulation() call. – hpaulj Nov 14 '15 at 0:10

scipy.optimize functions take an args parameter, as a way of passing extra arguments to the function. You could use that to pass a number of the parameters to heat_balance.

For example, let's assume r, T_a and dx are more 'variable' than R_w; that is, more likely to vary from run to run or test case. I'll also switch the order of T_1 and T_2, since the later is, apparently, the variable you want to minimize.

def heat_balance(T_2,T_1, T_a, r, dx):
heat_in = (T_1-T_2) * scipy.pi * r**2 * R_w
heat_out = ( (T_1 + T_2) / 2 - T_a ) * 2 * r * scipy.pi * dx * R_f
return abs(heat_in-heat_out)


Then in simulation minimize could be called with:

result = optimize.minimize(heat_balance, T[-1]-dx, (T[-1], T_a, r, dx),
method="SLSQP")
newT = result['x']  # not result.x?
# what's the purpose of the float() call?
T.append(newT)


The simulation definition could also include these 'global' variables.

If we are going to define object classes, I'd think more in terms of 'pipe' objects rather than 'simulation' objects. A 'pipe' would have certain properties - dimensions, conductivity etc, that could differ from one pipe to another. So it could be convenient to define several such objects, and pass them to the simulation.

I'd rewrite the simulation iteration as something like:

def simulation(...):
X = np.arange(dx, L, dx)
T = np.zeros_like(X)
for i, x in enumerate(X):
result = optimize.minimize(heat_balance, T[-1]-dx, (T[-1], T_a, r, dx), method="SLSQP")
newT = result['x']
T[i] = newT
return T,X   # return a simple tuple


And later call:

T,X = simulation(<parameters>)
plot(T, X,...)


I haven't tested these changes, but I think they give an adequate idea of how the code and calling structure could be cleaned up.

These changes don't address the speed issue. Repeatedly calling minimize over the length of pipe, does sound expensive.

Here's a version of your script that could make it easier to play with some of the parameters, such as the dx. I can import it into a shell (e.g. ipython session) and perform timings on pieces. I'm experimenting with 2 forms of the func. Precalculating values may say 20-30% time, not a change. dx has most effect on simulation times.

I also collect the nfev statistic, incase that gives any ideas of how costly the minimize call is. On average it calls the balance function 25 times per call, or 6000 times for a dx=.1.

Fine tuning the heat_balance equation helps some with speed, but I think being smarter with dx and the use of minimize help more.

import numpy as np
from scipy.optimize import minimize
from matplotlib import pyplot as plt

Rw = 0.654 # W/mk water heat loss ocefficient
Rf = 0.042 # W/mk pipe heat loss coefficient
r = 0.2 # m pipe radius
T0 = 50 # C initial temperature
Ta = 20 # C ambiant temperature
dx = 0.01 # m differential distance
L = 25 # m length of pipe
pi = np.pi
dt = 0   # tweak to last T

def heat_balance(T2, T1, Rf, r, Ta):
heat_in = (T1-T2) * pi * r**2 * Rw
heat_out = ( (T1 + T2) / 2 - Ta ) * 2 * r * pi * dx * Rf
return abs(heat_in-heat_out)

class Balance(object):
# heat calc with precalc
# seems to save 30% on time
def __init__(self, Rf, r, Ta):
self.coeff1 = pi * r**2 * Rw
self.coeff2 = 2 * r * pi * dx * Rf
self.Ta = Ta
def __call__(self, s, T1):
heat_in = (T1-s) * self.coeff1
heat_out = ( (T1+s)/2 - self.Ta) * self.coeff2
return abs(heat_in - heat_out)

def simulation1(Rf, T0, r, Ta):
# simulation with plain function
X = np.arange(0, L, dx)
T = np.zeros_like(X)
T[0] = lastT = T0
cnt = np.zeros_like(X)
func = heat_balance
for i,x in enumerate(X[1:]):
args = (lastT, Rf, r, Ta)
result = minimize(func, lastT-dt, args, method="SLSQP")
T[i+1] = lastT = result.x
cnt[i+1] = result.nfev
cnt = cnt[1:]
return X, T, cnt

def simulation2(Rf, T0, r, Ta):
# simulation with Balance class
X = np.arange(0, L, dx)
T = np.zeros_like(X)
T[0] = lastT = T0
cnt = np.zeros_like(X)
func = Balance(Rf, r, Ta)
for i,x in enumerate(X[1:]):
args = (lastT,)
result = minimize(func, lastT-dt, args, method="SLSQP")
T[i+1] = lastT = result.x
cnt[i+1] = result.nfev
cnt = cnt[1:]
return X, T, cnt

simulation = simulation2

if __name__ == '__main__':
import sys
if sys.argv[1:]:
dx = float(sys.argv[1])
if sys.argv[2:]:
dt = float(sys.argv[2]) # temp offset
Rfs = [7,14,50,100,200]
typs = ['+','x','|','_','*']
cols = ['red','blue','black','purple','brown']

for rf, typ, col in zip(Rfs, typs, cols):
Rf = 1.0/rf
print(Rf)
X, T, cnt = simulation(Rf, T0, r, Ta)
print('nfev count: ',len(cnt), cnt.min(), cnt.max(), cnt.sum(), cnt.mean())
plt.scatter(X, T, c=col, marker=typ, s = 100, label=rf, lw=3)

plt.legend(fontsize=20,title='Pipe Insolation R Value')
plt.xlabel('Distance (m)',fontsize=20)
plt.ylabel('Average Temperature (C)',fontsize=20)
plt.show()