Heatsink simulation — Thursday, April 29

This program sort of works, but the relaxation and scaling are both off. See the debugged and tweaked version, also posted.

heatsink.py

# -*- coding: utf-8 -*-
"""
Created on Mon Apr 26 10:52:51 2021

@author: rmontant
"""
import numpy as np

def pdf(dx):
    return ( np.exp(-(dx**2)/2) / np.sqrt(2*np.pi) )
#--------

class HeatSink:
    surround = -1e100
    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 = nfins * finwidth + self.ngaps * gapwidth
        self.length = self.width
        self.height = int( 1.2 * self.width ) # 1.2 seems like a good number

        ###
        ### 2021-04-22 afternoon:
        ### The cube MUST be initialized to the "surround" value here!
        ### (because it never happens anywhere else...)
        ###
        self.cube = np.full(
            (self.width + 2, self.length + 2, self.height + 2), self.surround )

        self.form_base()
        self.form_fins()
        self.build_pointarray()
        self.cube[:,:,0] = 2* self.initial_temperature

        self.tmax = self.initial_temperature
    #--------

    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
            self.cube[ \
                fstart:fend, 1:self.length-1 , self.baseheight+1:self.height-1 \
            ] = self.initial_temperature
    #--------

    def build_pointarray(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:
                        coords.append( (x,y,z) )
        self.coords = np.array( coords )
    #--------

    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_once(self, fraction=(1/8)):
        ###
        ### 2021-04-29:
        ### Make the new cube be a copy of the previous cube - this gets the
        ### hot base set properly, along with the surrounding cells:
        new_cube = self.cube.copy()

        for idx, p in enumerate(self.coords):
            w, l, h = p
            total = []
            for x in (-1,0,1):
                for y in (-1,0,1):
                    for z in (-1,0,1):
                        if self.cube[w+x, l+y, h+z] != self.surround:
                            total.append(self.cube[w+x, l+y, h+z])
                        else:
                            ###
                            ### 2021-04-29:
                            ### Let the edge points leak heat
                            ### into the surrounding points.
                            ### "Fraction" is a heuristic parameter
                            ### that governs the rate of heat loss.
                            neighbor_temp = \
                                (fraction * self.initial_temperature) \
                                    + ((1 - fraction) * self.cube[w,l,z])
                            total.append(neighbor_temp)
                            #pass

            new_cube[w,l,h] = sum(total)/len(total)
        return new_cube
    #--------

#--------

model-heatsink.py

# -*- coding: utf-8 -*-
"""
Created on Tue Apr 27 09:29:17 2021

@author: rmontant
"""
import sys, time
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.colors import Normalize
from mpl_toolkits.mplot3d import Axes3D

import heatsink

def draw_contour(axis, surface, tmin=None, tmax=None):
    if not tmin:
        tmin = np.min(surface)
    if not tmax:
        tmax = np.max(surface)
    normal = Normalize(tmin, tmax)
    pcolormap = axis.pcolormesh(surface, cmap=cm.afmhot, norm=normal)
    axis.contour(surface, cmap=cm.hot)
    plt.colorbar(pcolormap, ax=axis)
#--------

def rgb(v):
    return (np.sqrt(v), v, 0.9*v, 0.05)
#--------

def draw_heatsink(axis, heatsink):
    xcoords = [ p[0] for p in heatsink.coords]
    ycoords = [ p[1] for p in heatsink.coords]
    zcoords = [ p[2] for p in heatsink.coords]
    colors = [ rgb(0.6 * z / heatsink.height) for z in zcoords]
    axis.scatter(xcoords, ycoords, zcoords, color=colors, marker=',')
    axis.set_xlabel('width')
    axis.set_ylabel('length')
    axis.set_zlabel('height')
#--------

def display_results(heatsink):
    fig = plt.figure(figsize=(10,10))
    for layer in range(9):
        layer_idx = round((layer / 9) * heatsink.height)
        surface = heatsink.cube[ : , : , layer_idx]
        axis = fig.add_subplot(3,3, (layer+1))
        draw_contour(axis, surface,
                     tmin=heatsink.initial_temperature,
                     tmax=heatsink.tmax)
    plt.savefig('relaxed.png')
    plt.show()
#--------


def main(argv=[__name__]):
    htsnk = heatsink.HeatSink(nfins=2, finwidth=6, gapwidth=8, baseheight=10)

    nhotspots = 3
    for i in range(nhotspots):
        temperature = 100 + 10 * np.random.random()
        w = np.random.randint(1, htsnk.width-1)
        l = np.random.randint(1, htsnk.length-1)
        scale = 0.07 + 0.06 * np.random.random()
        print('hotspot {:d} is {:.3} at ({:d}, {:d})' \
              .format(i, scale, w, l))
        htsnk.set_hotspot(w, l, temperature, scale)

    fig = plt.figure(figsize=(8,8))
    a1 = fig.add_subplot(2, 2, 1)
    draw_contour(a1, htsnk.cube[:,:,0], \
                 tmin=htsnk.initial_temperature,
                 tmax=htsnk.tmax
    )

    a2 = fig.add_subplot(2,2, 2)
    surface = htsnk.cube[:, htsnk.length//2, :]
    draw_contour(a2, surface,
                 tmin=htsnk.initial_temperature,
                 tmax=htsnk.tmax,
    )

    a4 = fig.add_subplot(2,2, 4, projection='3d')
    draw_heatsink(a4, htsnk)

#    prows, pcols = 4, 4
#    fig2, axes = plt.subplots(nrows=prows, ncols=pcols, figsize=(8,8))
#    for i in range(prows):
#        for j in range(pcols):
#            zval = int( htsnk.height * (i*pcols + j) / (prows*pcols) )
#            surface = htsnk.cube[ :, :, zval]
#            draw_contour(axes[i][j], surface,
#                          tmin=htsnk.initial_temperature,
#                          tmax=htsnk.tmax,
#            )
    display_results(htsnk)

    timestart = time.process_time()
    changed = True
    passes = 0
    while changed:
        passes += 1
        print('\r{:5d} passes'.format(passes),end='')
        new_cube = htsnk.relax_once()
        changed = np.any( new_cube != htsnk.cube )
        # or...
        # changed = not np.all( new_heatsink.cube == htsnk.cube )
        htsnk.cube = new_cube
    print('\r{:5d} passes total'.format(passes))

    timedone = time.process_time()
    print('{:.3f} seconds elapsed'.format(timedone - timestart))

    display_results(htsnk)

    #plt.show()
#--------

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