HeatSink Simulation

HeatSink class

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# 2021-04-20 v0.2
"""
This is a MODULE containing the HeatSink and Framed_heatsink classes.
For inclusion by other programs.

2021-04-16 v0.1
@author: bobmon
"""
import numpy as np

#-----------------------------------------------------------
# functions to generate hotspots on the heatsink's base.

def pdf(dx):
    '''Gaussian distribution about a central point.'''
    return ( np.exp(-(dx**2)/2) / np.sqrt(2*np.pi) )
#--------

class HeatSink:
    surround = -1e3    # low enough
    initial_temperature = 20

    def __init__(self, nfins=4, finwidth=4, gapwidth=6, baseheight=5):
        self.nfins = nfins
        self.finwidth = finwidth
        self.gapwidth = gapwidth
        self.baseheight = baseheight
        self.ngaps = self.nfins - 1

        self.width = int(2 + self.nfins * self.finwidth + self.ngaps * self.gapwidth)
        self.length = int(1.8 * self.width)
        self.height = int(2 + 3 * self.baseheight)  # keep it short

        self.cube = np.full((self.width, self.length, self.height), self.surround)

        self.form_base()
        self.form_fins()
        self.build_pointarray()

        self.tmin = self.initial_temperature
        self.tmax = self.initial_temperature    # changes when hotspots are added
    #--------

    def form_base(self):
        self.cube[1:self.width-1, 1:self.length-1, 1:self.baseheight+1] \
            = self.initial_temperature

    def form_fins(self):
        for f in range(self.nfins):
            fstart = 1 + f * (self.gapwidth + self.finwidth)
            fend = fstart + self.finwidth
            print('fin {:d} - {:d}-{:d}'.format(f, fstart, fend))

            self.cube[fstart:fend, 1:self.length-1, self.baseheight+1:self.height-1] \
                = self.initial_temperature
    #--------

    def build_pointarray(self):
        self.coords = []
        for x in range(self.width):
            for y in range(self.length):
                for z in range(self.height):
                    if self.cube[x,y,z] > self.surround:
                        # Keep track of this point's coordinates:
                        self.coords.append( (x,y,z) )
    #--------

    def set_hotspot(self, w, l, t, scale):
        rmax = ((self.width**2 + self.length**2)**0.5) * scale
        for x in range(1, self.width-1):
            for y in range(1, self.length-1):
                r = ( (w-x)**2 + (l-y)**2 )**0.5
                s = pdf(r/rmax) / pdf(0)
                self.cube[x,y,0] = (s * t) + ((1-s) * self.cube[x,y,0])
        if t > self.tmax:
            self.tmax = t
    #--------

    def relax_precise(self, cool_ratio=0.1):
        '''
        "Relax" the heatsink by averaging the temperature of every point
        that's part of the heatsink, with its neighbors' temperatures.
        This version operates with maximum precision, at a cost in speed.
        '''
        warm_ratio = 1.0 - cool_ratio

        nextcube = self.cube.copy()
        for idx, p in enumerate(self.coords):
            x, y, z = tuple(p)
            neighborhood = []
            for i in (-1, 0, 1):
                for j in (-1, 0, 1):
                    for k in (-1, 0, 1):
                        vn = self.cube[x + i, y + j, z + k]
                        if vn > self.surround:
                            neighborhood.append(vn)
                        else:
                            # primitive modelling of cooling:
                            # make it proportional to the delta-t
                            cool_temp = (cool_ratio * self.initial_temperature) \
                                    + (warm_ratio * self.cube[x,y,z])
                            neighborhood.append(cool_temp)

            nextcube[x,y,z] = (sum(neighborhood) / len(neighborhood))

        changed = np.any( self.cube != nextcube )
        if changed:
            self.cube = nextcube
            #self.cube[1:self.width-1, 1:self.length-1, 1:self.height-1] \
            #    = nextcube[1:self.width-1, 1:self.length-1, 1:self.height-1]
        return changed
    #--------

    def relax_rounded(self, cool_ratio=0.1, resolution=2):
        '''
        "Relax" the heatsink by averaging the temperature of every point
        that's part of the heatsink, with its neighbors' temperatures.
        This version rounds calculations to a few decimal places, to trade
        precision for faster convergence.
        '''
        warm_ratio = 1.0 - cool_ratio

        nextcube = self.cube.copy()
        for idx, p in enumerate(self.coords):
            x, y, z = tuple(p)
            neighborhood = []
            for i in (-1, 0, 1):
                for j in (-1, 0, 1):
                    for k in (-1, 0, 1):
                        vn = self.cube[x + i, y + j, z + k]
                        if vn > self.surround:
                            neighborhood.append(vn)
                        else:
                            # primitive modelling of cooling:
                            # make it proportional to the delta-t
                            cool_temp = (cool_ratio * self.initial_temperature) \
                                    + (warm_ratio * self.cube[x,y,z])
                            neighborhood.append(cool_temp)

            nextcube[x,y,z] = sum(neighborhood) / len(neighborhood)

        changed = np.any(
            np.round(self.cube, resolution) != np.round(nextcube, resolution)
        )
        if changed:
            self.cube = nextcube
        return changed
    #--------
#--------

class FramedHeatSink(HeatSink):
    '''
    This class adds a wireframe of the heatsink volume.
    Also accepts an optional initial value.
    '''
    def __init__(self, \
        nfins=4, finwidth=4, gapwidth=6, baseheight=5,
        initial_value=None
    ):
        HeatSink.__init__(self,
            nfins=nfins, finwidth=finwidth, gapwidth=gapwidth, baseheight=baseheight
        )
        if initial_value:
            self.cube[ : , : , : ] = initial_value

        self.build_wireframe()

        #print('Framed_HeatSink __init__() done.')
    #--------

    def finframe(self, xa):
        '''Return a list of line segments outlining one fin.'''
        xb = xa + self.finwidth
        l = self.length # local name-of-convenience
        return [ \
            # verticals
            ((xa,xa), (0,0), (self.baseheight, self.height)),
            ((xb,xb), (0,0), (self.baseheight, self.height)),
            ((xa,xa), (l,l), (self.baseheight, self.height)),
            ((xb,xb), (l,l), (self.baseheight, self.height)),
            # y-laterals
            ((xb,xb), (0,l), (self.height, self.height)),
            ((xa,xa), (0,l), (self.height, self.height)),
            # x-laterals
            ((xa,xb), (0,0), (self.height, self.height)),
            ((xa,xb), (l,l), (self.height, self.height)),
        ]
    #--------

    def build_wireframe(self):
        '''Add a list of line segments that outline the heatsink,
        for each object.'''
        w, l, d = self.width, self.length, self.height   # local names-of-convenience
        gap_spacing = self.finwidth + self.gapwidth
        n_gap_positions = self.ngaps * gap_spacing

        wireframe = [ \
            # base
            ((0,w), (0,0), (0,0)),
            ((w,w), (0,l), (0,0)),
            ((0,w), (l,l), (0,0)),
            ((0,0), (0,l), (0,0)),
            # base + baseheight
            #((0,0), (0,l), (self.baseheight, self.baseheight)),
            #((w,w), (0,l), (self.baseheight, self.baseheight)),
            # join-base
            ((0,0), (l,l), (0, self.baseheight)),
            ((0,0), (0,0), (0, self.baseheight)),
            ((w,w), (0,0), (0, self.baseheight)),
            ((w,w), (l,l), (0, self.baseheight)),
            ((0,0), (l,l), (0, self.baseheight)),
        ]
        # gaps
        for gap_posn in range(self.finwidth, n_gap_positions, gap_spacing):
            gap_extent = min(self.width, gap_posn + self.gapwidth)
            gap_x_coords = (gap_posn, gap_extent)
            gap_z_coords = (self.baseheight, self.baseheight)
            wireframe.append( (gap_x_coords, (0,0), gap_z_coords) )
            wireframe.append( (gap_x_coords, (l,l), gap_z_coords) )

            wireframe.append( ((gap_posn, gap_posn), (0,l), gap_z_coords) )
            wireframe.append( ((gap_extent, gap_extent), (0,l), gap_z_coords) )
        # fins
        for f in range(self.nfins):
            wireframe += self.finframe(f * gap_spacing)

        self.wireframe = np.array(wireframe)
    #--------
#--------

model-heatsink.py

#!/usr/bin/env python3
# next-gen model of the heatsink
# 2021-04-24 Rearrange heatsink views, partially optimize for batch use
# 2021-04-22 Add data-file writing
# 2021-04-20
import sys
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.colors import Normalize
from mpl_toolkits.mplot3d import Axes3D   # needed for '3d" projection?

import heatsink

#-----------------------------------------------------------

def reddish(frac, transparency=0.1):
    return (np.sqrt(frac), 0.95*frac, 0.9*frac, transparency)
#--------

def bluish(frac, transparency=0.1):
    return (0.9*frac, 0.95*frac, np.sqrt(frac), transparency)
#--------

def draw_colormesh(axis, surface, tmin, tmax):
    normal = Normalize(tmin, tmax)
    pcolormap = axis.pcolormesh(surface, cmap=cm.hot, norm=normal)
    plt.colorbar(pcolormap, ax=axis)
#--------

# Produce a "scatter" plot of the coordinates of any object, viz. a heatsink:
def heatsink_hot_view(axis, heatsink):
    #
    # Define each layer's color based on its temperature.
    # This can be based on the initialized element values, as done in class,
    # or directly on the z-coordinate, as done here:
    xcoords = [p[0] for p in heatsink.coords]
    ycoords = [p[1] for p in heatsink.coords]
    zcoords = [p[2] for p in heatsink.coords]
    values = [heatsink.cube[p[0],p[1],p[2]] for p in heatsink.coords]
    color = [reddish(v/heatsink.tmax) for v in values]
    axis.scatter( xcoords, ycoords, zcoords,
        color=color,    # temperature-dependent color
        marker=',',     # use a small dot for each point
    )
    axis.set_title('Lighter is hotter')
    axis.set_xlabel('width')
    axis.set_ylabel('length')
    axis.set_zlabel('height')
#--------

# Produce a "scatter" plot of the coordinates of any object, viz. a heatsink:
def heatsink_cold_view(axis, heatsink):
    #
    # Define each layer's color based on its height (z-value).
    # This can be based on the initialized element values, as done in class,
    # or directly on the z-coordinate, as done here:
    xcoords = [p[0] for p in heatsink.coords]
    ycoords = [p[1] for p in heatsink.coords]
    zcoords = [p[2] for p in heatsink.coords]
    color = [bluish(0.80*z/heatsink.height) for z in zcoords]
    axis.scatter( xcoords, ycoords, zcoords,
        color=color,    # z-dependent color
        marker=',',     # use a small dot for each point
    )
    axis.set_xlabel('width')
    axis.set_ylabel('length')
    axis.set_zlabel('height')
    axis.set_title('Lighter is higher')
#--------

def display_heatsink(htsnk, cubeviewFtn, tmin=None, tmax=None):
    if not tmin:
        tmin = htsnk.tmin
    if not tmax:
        tmax = htsnk.tmax
    hsfig = plt.figure(figsize=(9, 3))  # figsize: (width X height)

    aF = hsfig.add_subplot(1, 3, 1)
    surface = htsnk.cube[ : , htsnk.length//2, : ]
    draw_colormesh(aF, surface, tmin=tmin, tmax=tmax)
    aF.contour(surface, cmap=cm.copper)
    aF.set_title('halfway through the width')

    aC = hsfig.add_subplot(1, 3, 2, projection='3d')
    cubeviewFtn(aC, htsnk)

    aB = hsfig.add_subplot(1, 3, 3)
    wslice = int(htsnk.finwidth + 0.5*htsnk.gapwidth)
    surface = htsnk.cube[ wslice, : , : ]
    draw_colormesh(aB, surface, tmin=tmin, tmax=tmax)
    aB.contour(surface, cmap=cm.copper)
    aB.set_title('through the length, between fins')
#--------

def display_cross_sections(htsnk, title, tmin=None, tmax=None):
    if not tmin:
        tmin = htsnk.tmin
    if not tmax:
        tmax = htsnk.tmax
    plotrows, plotcols = 4, 4

    hsfig, axes = plt.subplots(nrows=plotrows, ncols=plotcols, figsize=(10,10))
    #hsfig.tight_layout()

    for i in range(plotrows):
        for j in range(plotcols):
            z_value = int(htsnk.height * (i * plotcols + j) / (plotrows*plotcols))
            surface = htsnk.cube[ : , : , z_value]
            draw_colormesh(axes[i][j], surface, tmin=tmin, tmax=tmax)
            axes[i][j].contour(surface, cmap=cm.copper)
            axes[i][j].set_title(str(z_value))

    hsfig.suptitle(title)
#--------


def main(argv=[__name__]):
    if '-i' in argv:
        interactive = True
    else:
        interactive = False
    if '-S' in argv:
        saves = True
    else:
        saves = False

    if '-f' in argv:
        idx = argv.index('-f')
        nfins = int(argv[idx+1])
    else:
        nfins = 3
    if '-w' in argv:
        idx = argv.index('-w')
        finwidth = int(argv[idx+1])
    else:
        finwidth = 8
    if '-g' in argv:
        idx = argv.index('-g')
        gapwidth = int(argv[idx+1])
    else:
        gapwidth = 10
    if '-b' in argv:
        idx = argv.index('-b')
        baseheight = int(argv[idx+1])
    else:
        baseheight = 10

    if '-R' in argv:
        rounding = True
    else:
        rounding = False
    if '-r' in argv:
        idx = argv.index('-r')
        rounding_decimals = int(argv[idx+1])
    else:
        rounding_decimals = 2

    if '-o' in argv:
        idx = argv.index('-o')
        outfile = argv[idx+1]
    else:
        outfile = 'heatsink-relaxed.data'


    htsnk = heatsink.HeatSink(
        nfins=nfins, finwidth=finwidth, gapwidth=gapwidth, baseheight=baseheight)
    w, l, z = htsnk.cube.shape
    n_cpts = w * l * z
    n_hspts = len(htsnk.coords)
    print( '{:,} heatsink points in a volume of {:,} points, fraction is {:.3f}' \
        .format(n_hspts, n_cpts, n_hspts/n_cpts) )


    hotside = 100
    htsnk.cube[1:htsnk.width-1, 1:htsnk.length-1, 0] = hotside
    nhotspots = 3
    for i in range(nhotspots):
        temperature = hotside + 100 * np.random.random()
        w = np.random.randint(1, htsnk.width-1)
        l = np.random.randint(1, htsnk.length-1)
        scale = 0.12 + 0.03 * np.random.random()
        print('hotspot {:d} is {:.3} at ({:d}, {:d})' .format(i, scale, w, l))
        htsnk.set_hotspot(w, l, temperature, scale)

    tmin_initial = hotside - 15
    tmin_display = htsnk.initial_temperature + 10
    tmax_display = htsnk.tmax


    display_heatsink(htsnk, heatsink_cold_view)
    if saves:
        plt.savefig('heatsink-views-initial.png')

    count = 0
    display_cross_sections(htsnk, 'iteration {:4d}'.format(count),
        tmin=tmin_display, tmax=tmax_display)
    if saves:
        plt.savefig('{:04d}.png'.format(count))
    if interactive:
        plt.show()

    sequence = 1
    changed = True
    if rounding:
        relax_once = lambda cool_ratio: \
                        htsnk.relax_rounded(cool_ratio=cool_ratio,
                        resolution=rounding_decimals)
    else:
        relax_once = lambda cool_ratio: \
                        htsnk.relax_precise(cool_ratio=cool_ratio)

    while changed:
        count += 1
        changed = relax_once(cool_ratio=0.05)
        if count == sequence:
            display_cross_sections(htsnk, 'iteration {:d}'.format(count),
                tmin=tmin_display, tmax=tmax_display)
            if saves:
                plt.savefig('{:04d}.png'.format(count))
            if interactive:
                plt.show()

            sequence *= 2
        print('\r{:5d} relaxations'.format(count), end='')
    print()


    display_cross_sections(htsnk, 'iteration {:d} (Final)'.format(count),
        tmin=tmin_display, tmax=tmax_display)
    if saves:
        plt.savefig('{:04d}.png'.format(count))

    display_heatsink(htsnk, heatsink_hot_view)
    if saves:
        plt.savefig('heatsink-views.png')
    if interactive:
        plt.show()

    data = [
        '# After {:d} iterations:'.format(count),
        '# {:3d} {:3d} {:3d}'.format(htsnk.width, htsnk.length, htsnk.height),
        '# {:3d} {:3d} {:3d}'.format(htsnk.baseheight, htsnk.nfins, htsnk.finwidth),
        '# {:e}  {:e}'.format(htsnk.initial_temperature, htsnk.surround),
    ]
    data += [ \
        '{:3d} {:3d} {:3d}  {:e}'.format(p[0], p[1], p[2], htsnk.cube[p[0],p[1],p[2]]) \
        for p in htsnk.coords ]
    data.append('')

    with open(outfile, 'w') as h:
        h.write('\n'.join(data))
    print('Data saved as {:s}'.format(outfile))
#--------

if __name__ == '__main__':
    sys.exit(main(sys.argv))
#--------