Using GEKKO to solve 2-D Heat Equations

Question

I am basically trying to solve a dynamic second order partial differential equation using GEKKO. This is the equation for reference: 2-D Heat transfer equation.

Here, t is time, T is temperature, (k, rho, Cp,l e, sigma and Z) are all constants.

T has been inputted as an array (I think that is where the problem lies).

To solve the equation numerically, boundary conditions are required which were inputted as other equations.

This is the code:

import numpy as np
from gekko import GEKKO

m = GEKKO()

tf = 5*60*60
dt = int(tf/300)+1
m.time = np.linspace(0,tf,dt)

# Number of nodes
n = 41

# Length of domain
Lx = 1
Ly = Lx # square domain

x_div = np.linspace(0,Lx,n)
y_div = np.linspace(Ly,0,n)

[X, Y] = np.meshgrid(x_div, y_div)

# step size
dx = x_div[1] - x_div[0]
dy = y_div[1] - y_div[0]

# Temp. initialization
T = m.Var(np.ones((n,n))*290)

# Equation set-up

# Middle segments
for i in range(1,n-1):
    for j in range(1,n-1):
        m.Equations(T[i,j].dt() == (k/(rho*Cp)*((T[i+1,j]-2*T[i,j]+T[i-1,j])/dx**2 + (T[i,j+1]-2*T[i,j]+T[i,j-1])/dy**2))
                    + (G_total - em*sigma*(T[i,j]**4-T_surr**4))/(rho*Cp*Z))

# Boundary Conditions
m.Equations(T[0,:]==310,
           T[-1,:]==310,
           T[1:-2,0]==315,
           T[1:-2,-1]==315,
           
           T[0,0]==312,
           T[n,0]==312,
           T[0,n]==312,
           T[n,n]==312)

Basically, I am trying to solve this meshgrid consisting of temperatures. I get the following error: 'numpy.float64' object has no attribute 'dt'

If I just write T instead of T[i,j], I get this error: 'int' object is not subscriptable

My questions:

1. Is GEKKO able to solve such equations, that are 2 dimensional in nature? How do I go about it?
2. Are there any other cool libraries out there for this purpose? I need to be able to draw a contour plot having temperatures of the plate as time progresses; for which I need to solve the equations.

Thank you for your time.

Answer 1

Here is a solution with Lutz Lehmann's suggestion to use m.Array to define T.

import numpy as np
from gekko import GEKKO

m = GEKKO(remote=False)

tf = 5*60*60
dt = int(tf/300)+1
#m.time = np.linspace(0,tf,dt)
# for testing
m.time = [0,0.01,0.02,0.03]

# Number of nodes
n = 41

# Length of domain
Lx = 1
Ly = Lx # square domain

# Define for material
k = 1; rho = 8000; Cp = 500
G_total=1; em=1; sigma=1
T_surr=298; Z=1

x_div = np.linspace(0,Lx,n)
y_div = np.linspace(Ly,0,n)

[X, Y] = np.meshgrid(x_div, y_div)

# step size
dx = x_div[1] - x_div[0]
dy = y_div[1] - y_div[0]

# Temp. initialization
T = m.Array(m.Var,(n,n),value=290)

# Equation set-up

# Middle segments
for i in range(1,n-1):
    for j in range(1,n-1):
        m.Equation(rho*Cp*T[i,j].dt() == (k*\
                                    ((T[i+1,j]-2*T[i,j]+T[i-1,j])/dx**2 \
                                     + (T[i,j+1]-2*T[i,j]+T[i,j-1])/dy**2))
                    + (G_total - em*sigma*(T[i,j]**4-T_surr**4))/Z)

# Boundary Conditions
m.Equations([T[0,i]==310 for i in range(1,n-1)])
m.Equations([T[-1,i]==310 for i in range(1,n-1)])
m.Equations([T[i,0]==315 for i in range(1,n-1)])
m.Equations([T[i,-1]==315 for i in range(1,n-1)])
m.Equations([T[0,0]==312, T[n-1,0]==312, \
             T[0,n-1]==312, T[n-1,n-1]==312])

m.options.IMODE=7
m.solve(disp=False)

import matplotlib.pyplot as plt

plt.figure(figsize=(8,8))
for i in range(0,4):
    for j in range(0,4):
        plt.subplot(4,4,i*4+j+1)
        plt.plot(m.time,T[i,j].value)
plt.savefig('heat.png',dpi=600)
plt.show()

There are additional examples of hyperbolic and parabolic PDEs with Gekko. The strength of Gekko is with optimization. For simulation, it may be better to use standard simulation methods for PDEs. Also, you may get faster solutions for simulation using IMODE=7.