0
$\begingroup$

I'm trying to write a Python code that is a numerical solver for 1-d heat conduction (using FVM) with a temperature dependent thermal conductivity.

The solver has three functions I need to iterate until convergence. This is basically the algorithm from Patankar's Numerical Heat Transfer and Fluid Flow book.

Conductivity function: takes an initial guess for temperature field and outputs nodal conductivities, dependent on past value of nodal temperatures.

Nodal coefficients: calculates nodal coefficients based on current nodal conductivities.

Nodal temperatures: uses a TDMA algorithm to calculate the nodal temperatures based on current coefficients.

Once the temperatures are calculated, the values are to be cycled through to the conductivity function to continue the loop until convergence.

I can make this work just fine when the conductivity is constant, but start getting issues when I try to implement this pseudocode. I get a maximum recursion depth exceeded error. My code is below:

import numpy as np
import matplotlib.pyplot as plt
import copy

# SOLUTION ALGORITHM

# Step 1: DATA  
def data(n, hw, he, TinfW, TinfE):    
    n = 5
    hw = 10**10
    he = 10**-10
    TinfW = 100
    TinfE = 10**12

    return n, hw, he, TinfW, TinfE    

# Step 2: GRID    
def grid(x1, x2, data):
    x1 = 0 
    x2 = 0.5
    xe = (x2-x1) / data(n, hw, he, TinfW, TinfE)[0]

    return xe

# Step 3: GAMSOR     
def gamsor(k, q, TDMA):
    k0 = 100
    beta = 20
    T0 = 25   
    for i in range(data(n, hw, he, TinfW, TinfE)[0]):
        k[i] = k0 + beta#*(TDMA(data, coeff)[i]-T0)
        q[i] = 10**5

    return k, q

 #Step 4: COEFF     
def coeff(Ae, Aw, Ap, data, gamsor, grid):
    for i in range(data(n, hw, he, TinfW, TinfE)[0]):       
        Ae[i] = Aw[i] = gamsor(k,q, TDMA(data, coeff))[0][i] /     grid(x1, x2, data)
        Ap[i] = Ae[i] + Aw[i]  

# Define "b" array - including source term and boundary conditions
        b = np.zeros((data(n, hw, he, TinfW, TinfE)[0], 1))
        b.fill(gamsor(k,q, TDMA(data, coeff))[1][i])    

        b[0] = (gamsor(k,q, TDMA(data, coeff))[1][i]*(grid(x1,x2,data)/2) + 
           data(n, hw, he, TinfW, TinfE)[1]*
           data(n, hw, he, TinfW, TinfE)[3])

        b[-1] = (gamsor(k,q, TDMA(data, coeff))[1][i]*(grid(x1,x2,data)/2)+ 
            data(n, hw, he, TinfW, TinfE)[2]*
            data(n, hw, he, TinfW, TinfE)[4])

# Change values in first and last coefficients to reflect specified BC's
        Aw[0] = 0
        Ap[0] = Ae[0] + Aw[0] + data(n, hw, he, TinfW, TinfE)[1]

        Ae[-1] = 0
        Ap[-1] = Ae[0] + Aw[0] + data(n, hw, he, TinfW, TinfE)[2]

    return Ae, Aw, Ap, b

 #Step 5: SOLVE
def TDMA(data, coeff):
# Initialize "T" array - the solution array
    T = np.zeros(data(n, hw, he, TinfW, TinfE)[0]);

# Initialize recursion functions as arrays
    P = np.zeros(data(n, hw, he, TinfW, TinfE)[0]);
    Q = np.zeros(data(n, hw, he, TinfW, TinfE)[0]);

## Step 1: evaluate at node 1
    P[0] = (coeff(Ae, Aw, Ap, data, gamsor, grid)[0][0] / 
       coeff(Ae, Aw, Ap, data, gamsor, grid)[2][0])

    Q[0] = (coeff(Ae, Aw, Ap, data, gamsor, grid)[3][0] / 
       coeff(Ae, Aw, Ap, data, gamsor, grid)[2][0])

## Step 2: sweep from node 2 to node (n-1) (python
## cell 1->(n-2) ) evaluating P and Q

    for i in range(1, data(n, hw, he, TinfW, TinfE)[0]-1):
        P[i] = ((coeff(Ae, Aw, Ap, data, gamsor, grid)[0][i]) / 
           (coeff(Ae, Aw, Ap, data, gamsor, grid)[2][i] - 
           coeff(Ae, Aw, Ap, data, gamsor, grid)[1][i]*P[i-1]))

        Q[i] = ((coeff(Ae, Aw, Ap, data, gamsor, grid)[3][i] + 
           coeff(Ae, Aw, Ap, data, gamsor, grid)[1][i]*Q[i-1]) / 
           (coeff(Ae, Aw, Ap, data, gamsor, grid)[2][i] - 
           coeff(Ae, Aw, Ap, data, gamsor, grid)[1][i]*P[i-1]))

## Step 3: calculate for node n

    P[data(n, hw, he, TinfW, TinfE)[0]-1] = 0

    Q[data(n, hw, he, TinfW, TinfE)[0]-1] = ( 
(coeff(Ae, Aw, Ap, data, gamsor, grid)[3][data(n, hw, he, 
TinfW, TinfE)[0]-1] + coeff(Ae, Aw, Ap, data, gamsor, grid)
[1][data(n, hw, he, TinfW, TinfE)[0]-1]*Q[data(n, hw, he, TinfW, 
TinfE)[0]-2]) / (coeff(Ae, Aw, Ap, data, gamsor, grid)[2][data(n, 
hw, he, TinfW, TinfE)[0]-1] - coeff(Ae, Aw, Ap, data, gamsor, 
grid)[1][data(n, hw, he, TinfW, TinfE)[0]-1]*P[data(n, hw, he, 
TinfW, TinfE)[0]-2]))

    T[data(n, hw, he, TinfW, TinfE)[0]-1] = ( 
    Q[data(n, hw, he, TinfW, TinfE)[0]-1] )

## Step 4: back fill, giving the temperatures
    for i in range(data(n, hw, he, TinfW, TinfE)[0] - 2, -1, -1):
        T[i] = P[i]*T[i+1] + Q[i]

    return T

def Temps(TDMA):
    Temps = copy.copy(TDMA(data,coeff)); Temps.fill(100)
    return Temps

if __name__ == '__main__': 

# Initialize variables
    n = 0
    hw = 0
    he = 0
    TinfW = 0
    TinfE = 0

    k = np.zeros(data(n, hw, he, TinfW, TinfE)[0])
    q = np.zeros(data(n, hw, he, TinfW, TinfE)[0])
    k0 = 0
    beta = 20
    T0 = 0

    x1 = 0
    x2 = 0

    Ae = np.zeros(data(n, hw, he, TinfW, TinfE)[0])
    Aw = np.zeros(data(n, hw, he, TinfW, TinfE)[0])
    Ap = np.zeros(data(n, hw, he, TinfW, TinfE)[0])
    b = 0

    data(n, hw, he, TinfW, TinfE)
    grid(x1, x2, data)

    Temps(TDMA)
    for _ in range(1):
        gamsor(k, q)
        coeff(Ae, Aw, Ap, data, gamsor, grid)
        TDMA(data, coeff)

#print Temps(TDMA) 



    print "Nodal Temperature Distribution (C):", TDMA(data, coeff) 
$\endgroup$
  • 3
    $\begingroup$ You should get rid of the recursion. Python doesn't do tail recursion optimization and has a maximum call depth. Use a loop. $\endgroup$ – Bill Barth Feb 3 '15 at 18:51
  • 1
    $\begingroup$ Your code is distinctly unpythonic, so I strongly recommend codereview.stackexchange.com $\endgroup$ – Kirill Feb 4 '15 at 10:13
  • $\begingroup$ Ok thanks for the code review link. I'm still really new at this. I'll try to figure out how to get rid of the recursion and use a loop. $\endgroup$ – Josh Feb 4 '15 at 12:46
2
$\begingroup$

If you are interested in a pre-existing solver you should take a look a FiPy. It's a Finite-Volume method solver implemented in Python.

http://www.ctcms.nist.gov/fipy/

This is a 1D example: (http://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.mesh1D.html)

The example shows how to work with time-dependent b.c.'s and varying diffusion coefficients. It might be useful for you.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.