#! /usr/bin/env python3 # def boundary ( nx, ny, x, y, A, rhs ): #*****************************************************************************80 # ## boundary() sets up the matrix and right hand side at boundary nodes. # # Discussion: # # For this simple problem, the boundary conditions specify that the solution # is 10 on the left size, 100 on the right side, and 0 on the top and bottom. # # Nodes are assigned a single index K, which increases as: # # (NY-1)*NX+1 (NY-1)*NX+2 ... NY * NX # .... .... ... ..... # NX+1 NX+2 ... 2 * NX # 1 2 ... NX # # The index K of a node on the lower boundary satisfies: # 1 <= K <= NX # The index K of a node on the upper boundary satisfies: # (NY-1)*NX+1 <= K <= NY * NX # The index K of a node on the left boundary satisfies: # mod ( K, NX ) = 1 # The index K of a node on the right boundary satisfies: # mod ( K, NX ) = 0 # # If we number rows from bottom I = 1 to top I = NY # and columns from left J = 1 to right J = NX, then the relationship # between the single index K and the row and column indices I and J is: # K = ( I - 1 ) * NX + J # and # J = 1 + mod ( K - 1, NX ) # I = 1 + ( K - J ) / NX # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 12 March 2017 # # Author: # # John Burkardt # # Input: # # integer NX, NY, the number of grid points in X and Y. # # real X(NX), Y(NY), the coordinates of grid lines. # # real sparse A(N,N), the system matrix, with the entries for the # interior nodes filled in. # # real RHS(N), the system right hand side, with the entries for the # interior nodes filled in. # # Output: # # real sparse A(N,N), the system matrix, with the entries for all # nodes filled in. # # real RHS(N), the system right hand side, with the entries for # all nodes filled in. # # # Left boundary. # j = 0 for i in range ( 0, ny ): kc = i * nx + j xc = x[j] yc = y[i] A[kc,kc] = 1.0 rhs[kc] = 10.0 # # Right boundary. # j = nx - 1 for i in range ( 0, ny ): kc = i * nx + j xc = x[j] yc = y[i] A[kc,kc] = 1.0 rhs[kc] = 100.0 # # Lower boundary. # i = 0 for j in range ( 0, nx ): kc = i * nx + j xc = x[j] yc = y[i] A[kc,kc] = 1.0 rhs[kc] = 0.0 # # Upper boundary. # i = ny - 1 for j in range ( 0, nx ): kc = i * nx + j xc = x[j] yc = y[i] A[kc,kc] = 1.0 rhs[kc] = 0.0 return A, rhs def fd2d_heat_steady ( nx, ny, x, y, d, f ): #*****************************************************************************80 # ## fd2d_heat_steady() solves the steady 2D heat equation. # # Discussion: # # Nodes are assigned a singled index K, which increases as: # # (NY-1)*NX+1 (NY-1)*NX+2 ... NY * NX # .... .... ... ..... # NX+1 NX+2 ... 2 * NX # 1 2 ... NX # # Therefore, the neighbors of an interior node numbered C are # # C+NY # | # C-1 --- C --- C+1 # | # C-NY # # Nodes on the lower boundary satisfy: # 1 <= K <= NX # Nodes on the upper boundary satisfy: # (NY-1)*NX+1 <= K <= NY * NX # Nodes on the left boundary satisfy: # mod ( K, NX ) = 1 # Nodes on the right boundary satisfy: # mod ( K, NX ) = 0 # # If we number rows from bottom I = 1 to top I = NY # and columns from left J = 1 to right J = NX, we have # K = ( I - 1 ) * NX + J # and # J = 1 + mod ( K - 1, NX ) # I = 1 + ( K - J ) / NX # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 12 March 2017 # # Author: # # John Burkardt # # Input: # # integer NX, NY, the number of grid points in X and Y. # # real X(NX), Y(NY), the coordinates of grid lines. # # function D(X,Y), evaluates the thermal conductivity. # # function F(X,Y), evaluates the heat source term. # # Output: # # real U(NX,NY), the approximation to the solution at the grid points. # import numpy as np # # Set the total number of unknowns. # n = nx * ny # # Allocate the matrix and right hand side. # A = np.zeros ( [ n, n ] ) rhs = np.zeros ( n ) # # Define the matrix at interior points. # A, rhs = interior ( nx, ny, x, y, d, f, A, rhs ) # # Handle boundary conditions. # A, rhs = boundary ( nx, ny, x, y, A, rhs ) # for i in range ( 0, 5 ): # for j in range ( 0, nx ): # kc = i * nx + j # xc = x[j] # yc = y[i] # print ( ' %d %d %d %f %f %g %g' % ( i, j, kc, xc, yc, A[kc,kc], rhs[kc] ) ) # # Solve the linear system. # u = np.linalg.solve ( A, rhs ) u.shape = ( ny, nx ) return u def fd2d_heat_steady_test ( ): #*****************************************************************************80 # ## fd2d_heat_steady_test() tests fd2d_heat_steady(). # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 27 August 2013 # # Author: # # John Burkardt # print ( '' ) print ( 'fd2d_heat_steady_test:' ) print ( ' Python version' ) print ( ' Test fd2d_heat_steady().' ) fd2d_heat_steady_test01 ( ) # # Terminate. # print ( '' ) print ( 'fd2d_heat_steady_test:' ) print ( ' Normal end of execution.' ) return def fd2d_heat_steady_test01 ( ): #*****************************************************************************80 # ## fd2d_heat_steady_test01() demonstrates the use of fd2d_heat_steady. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 12 March 2017 # # Author: # # John Burkardt # import numpy as np from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt from matplotlib import cm # # Specify the spatial grid. # nx = 21 ny = 11 xvec = np.linspace ( 0.0, 2.0, nx ) yvec = np.linspace ( 0.0, 1.0, ny ) # # Solve the finite difference approximation to the steady 2D heat equation. # umat = fd2d_heat_steady ( nx, ny, xvec, yvec, d, f ) # # Plotting. # xmat, ymat = np.meshgrid ( xvec, yvec ) fig = plt.figure() ax = fig.add_subplot ( 111, projection = '3d' ) ax.plot_surface ( xmat, ymat, umat, cmap = cm.coolwarm, linewidth = 0, antialiased = False ) ax.set_xlabel ( '<--- Y --->' ) ax.set_ylabel ( '<--- X --->' ) ax.set_zlabel ( '<---U(X,Y)--->' ) ax.set_title ( 'Solution of steady heat equation' ) plt.draw ( ) filename = 'fd2d_heat_steady_test01.png' fig.savefig ( filename ) plt.show ( block = False ) plt.close ( ) print ( '' ) print ( ' Plotfile saved as "%s".' % ( filename ) ) return def d ( x, y ): #*****************************************************************************80 # ## d() evaluates the heat conductivity coefficient. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 23 July 2013 # # Author: # # John Burkardt # # Input: # # real X, Y, the evaluation point. # # Output: # # real VALUE, the value of the heat conductivity at (X,Y). # value = 1.0 return value def f ( x, y ): #*****************************************************************************80 # ## f() evaluates the heat source term. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 23 July 2013 # # Author: # # John Burkardt # # Input: # # real X, Y, the evaluation point. # # Output: # # real VALUE, the value of the heat source term at (X,Y). # value = 0.0 return value def interior ( nx, ny, x, y, d, f, A, rhs ): #*****************************************************************************80 # ## interior() sets up the matrix and right hand side at interior nodes. # # Discussion: # # Nodes are assigned a single index K, which increases as: # # (NY-1)*NX+1 (NY-1)*NX+2 ... NY * NX # .... .... ... ..... # NX+1 NX+2 ... 2 * NX # 1 2 ... NX # # Therefore, the neighbors of an interior node numbered C are # # C+NY # | # C-1 --- C --- C+1 # | # C-NY # # If we number rows from bottom I = 1 to top I = NY # and columns from left J = 1 to right J = NX, then the relationship # between the single index K and the row and column indices I and J is: # K = ( I - 1 ) * NX + J # and # J = 1 + mod ( K - 1, NX ) # I = 1 + ( K - J ) / NX # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 12 March 2017 # # Author: # # John Burkardt # # Input: # # integer NX, NY, the number of grid points in X and Y. # # real X(NX), Y(NY), the coordinates of grid lines. # # function pointer @D(X,Y), evaluates the thermal conductivity. # # function pointer @F(X,Y), evaluates the heat source term. # # real sparse A(N,N), the system matrix, without any entries set. # # real RHS(N), the system right hand side, without any entries set. # # Output: # # real sparse A(N,N), the system matrix, with the entries for the # interior nodes filled in. # # real RHS(N), the system right hand side, with the entries for the # interior nodes filled in. # import numpy as np # # For now, assume X and Y are equally spaced. # dx = x[1] - x[0] dy = y[1] - y[0] for ic in range ( 1, ny - 1 ): for jc in range ( 1, nx - 1 ): ino = ic + 1 iso = ic - 1 je = jc + 1 jw = jc - 1 kc = ic * nx + jc ke = kc + 1 kw = kc - 1 kn = kc + nx ks = kc - nx dce = d ( 0.5 * ( x[jc] + x[je] ), y[ic] ) dcw = d ( 0.5 * ( x[jc] + x[jw] ), y[ic] ) dcn = d ( x[jc], 0.5 * ( y[ic] + y[ino] ) ) dcs = d ( x[jc], 0.5 * ( y[ic] + y[iso] ) ) A[kc,kc] = ( dce + dcw ) / dx / dx + ( dcn + dcs ) / dy / dy A[kc,ke] = - dce / dx / dx A[kc,kw] = - dcw / dx / dx A[kc,kn] = - dcn / dy / dy A[kc,ks] = - dcs / dy / dy rhs[kc] = f ( x[jc], y[ic] ) return A, rhs def timestamp ( ): #*****************************************************************************80 # ## timestamp() prints the date as a timestamp. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 06 April 2013 # # Author: # # John Burkardt # import time t = time.time ( ) print ( time.ctime ( t ) ) return None if ( __name__ == '__main__' ): timestamp ( ) fd2d_heat_steady_test ( ) timestamp ( )