<p>也许您可以尝试一种更受测试驱动的方法。一个帮助是Python<strong>slice</strong>-命令。它在numpy数组访问中速度很快,并且在有限差分中使用时不太容易出错。试一试:</p>
<pre><code>import numpy as np
def get_slices(mx, my, mz):
jxw, jxp, jxe = slice(0,mx-2), slice(1,mx-1), slice(2,mx) # west, center, east
jys, jyp, jyn = slice(0,my-2), slice(1,my-1), slice(2,my) # south, center, north
jzb, jzp, jzt = slice(0,mz-2), slice(1,mz-1), slice(2,mz) # botoom,center, top
return jxw, jxp, jxe, jys, jyp, jyn, jzb, jzp, jzt
nx, ny, nz = 15, 10, 15
jxw, jxp, jxe, \
jys, jyp, jyn, \
jzb, jzpc, jzt = get_slices(nx, ny, nz)
ax, ay, az = np.arange(nx), np.arange(ny), np.arange(nz)
print('axW', ax[jxw])
print('axP', ax[jxp])
print('axE', ax[jxe])
axW [ 0 1 2 3 4 5 6 7 8 9 10 11 12]
axP [ 1 2 3 4 5 6 7 8 9 10 11 12 13]
axE [ 2 3 4 5 6 7 8 9 10 11 12 13 14]
</code></pre>
<p>我不太清楚舒尔<strong>您代码中的图形显示了什么。因此,通过使用不同的T(或w2)初始化和不完全相同的<code>kx, ky, kz</code>值仔细检查它。注意<code>np.meshgrid()</code>命令中的<code>indexing="ij"</code>选项。使用<code>np.nan</code>在某种程度上是危险的,因为它改变了变量globaly</p>
<pre><code>import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
def get_grid(mx, my, mz, Lx,Ly,Lz):
ix, iy, iz = Lx*np.linspace(0,1,mx), Ly*np.linspace(0,1,my), Lz*np.linspace(0,1,mz)
x, y, z = np.meshgrid(ix,iy,iz, indexing='ij')
print('ix', ix), print('iy', iy), print('iz', iz)
return x,y,z
def plot_grid(x,y,z,T):
def plot_boundary_only(x,y,z,T):
mx, my, mz = x.shape
x[1:-1, 1:-1, 1:-1],y[1:-1, 1:-1, 1:-1],z[1:-1, 1:-1, 1:-1],T[1:-1, 1:-1, 1:-1] = np.nan, np.nan, np.nan, np.nan # remove interior
x[1:-1, 1:,0], y[1:-1, 1:,0], z[1:-1, 1:,0], T[1:-1, 1:,0] = np.nan, np.nan, np.nan, np.nan # remove bottom
x[1:-1, my-1, 1:-1], y[1:-1, my-1, 1:-1], z[1:-1, my-1, 1:-1], T[1:-1, my-1, 1:-1] = np.nan, np.nan, np.nan, np.nan # remove north face
x[0, 1:, :-1], y[0, 1:, :-1], z[0, 1:, :-1], T[0, 1:, :-1] = np.nan, np.nan, np.nan, np.nan # remove west face
return x,y,z,T
x,y,z,T = plot_boundary_only(x,y,z,T)
fig = plt.figure(figsize=(15,15))
ax = fig.add_subplot(111, projection='3d')
img = ax.scatter(x,y,z, c=T.reshape(-1), s=150, cmap=plt.jet())
fig.colorbar(img)
#ax.zaxis.set_rotate_label(False) # To disable automatic label rotation
ax.set_ylabel('Y')
ax.set_xlabel('X', rotation=0)
ax.set_zlabel('Z', rotation=0)
plt.savefig("cube.png")
plt.show()
def init_T(x,y,z):
T = np.zeros_like(x)
case = 'xz'
if case == 'x': T = x
if case == 'y': T = y
if case == 'z': T = z
if case == 'xz': T = x+z
return T
nx, ny, nz = 35, 20, 30
Lx, Ly, Lz = nx-1, ny-1, nz-1
x,y,z = get_grid(nx, ny, nz, Lx,Ly,Lz) # generate a grid with mesh size Δx = Δy = Δz = 1
T = init_T(x,y,z)
plot_grid(x,y,z,T)
</code></pre>
<p><a href="https://i.stack.imgur.com/oE5Fp.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/oE5Fp.png" alt="enter image description here"/></a></p>
<p>这是运行热方程求解器的主程序。请注意,我们不能使用在曲面上显示点的图形来查看任何结果!Neumann条件位于我们看不到的底部,面上的值是不变的定值边界条件</p>
<pre><code>def set_GradBC(u):
# - Neuman BC at bottom boundary -
u[:, :, 0] = u[:, :, 1]
return u
def dT(u,DU):
DU[jxp,jyp,jzp] = (u[jxe,jyp,jzp] - 2.0*u[jxp,jyp,jzp] + u[jxw,jyp,jzp])*Dx + \
(u[jxp,jyn,jzp] - 2.0*u[jxp,jyp,jzp] + u[jxp,jys,jzp])*Dy + \
(u[jxp,jyp,jzt] - 2.0*u[jxp,jyp,jzp] + u[jxp,jyp,jzb])*Dz
return DU
#mmmmmmmmmmmmm main mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm
# set the parameters
dt = 0.01 # time step size
nIterations = 5000 # number of iterations
nx, ny, nz = 40,20,30 # number of grid points
Lx, Ly, Lz = nx-1, ny-1, nz-1 # physical grid size
Kx, Ky, Kz = 1, 1, 1 # diffusion coeficients
dx, dy, dz = Lx/(nx-1), Ly/(ny-1), Lz/(nz-1) # mesh size
Dx, Dy, Dz = Kx/dx**2, Ky/dy**2, Kz/dz**2 # equation coefficients
# initialize the field variables -
x,y,z = get_grid(nx, ny, nz, Lx,Ly,Lz) # generate the grid
T = init_T(x,y,z) # initialize the temperature
DT = np.zeros_like(T) # initialize the change of temperature
jxw, jxp, jxe, \
jys, jyp, jyn, \
jzb, jzp, jzt = get_slices(nx, ny, nz) # slices for finite differencing
# run the solving algorithm
for jter in range(nIterations):
T = set_GradBC(T) # set Neumann boundary condition
T = T + dt*dT(T,DT) # update
</code></pre>
<p>下面是一个切片,沿JY1=常数显示结果:</p>
<pre><code>def plot_y_slice(JY1, x,y,z,T):
T2 = T[:,JY1,:]
X = x[:,JY1,:]
Z = z[:,JY1,:]
fig = plt.figure(figsize=(15,15))
ax = fig.add_subplot(111)
ax.set_aspect('equal')
plt.contourf(X,Z,T2, 40,alpha=0.99)
plt.savefig("y_slice.png")
plt.show()
JY1 = 10
plot_y_slice(JY1, x,y,z,T)
</code></pre>
<p><a href="https://i.stack.imgur.com/YdSFq.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/YdSFq.png" alt="enter image description here"/></a></p>