# Geol 5470, Humphrey 2022 # Our first, very straightforward 2D Finite Difference problem # uses constant coeffiecients and is steady state # Rectangular problem domain. Simplest approach # full (not sparse) matrices, with BCnodes included in matrix # Solve temperature field in a rectangular piece of lithosphere, assuming # we know all the bounding temperatures. So it is steady state, fixed BCs # the only difference here is the plotting in 3D import matplotlib.pyplot as plt import numpy as np from mpl_toolkits.mplot3d import Axes3D from matplotlib import cm #***** our problem assumes the block is 1500m wide, 1000m deep, with a sfc # temperature of 1, base temperature of 7 degrees, right side of 3, and # left side of 6 degrees. # Note, since K is constant, we can divide it out of the problem, # also, less obvious, we could divide the dimensions out of the problem # if we make delx and delz the same!! (this requires that the ratio of the # number of nodes in x over the nodes in z be the same as the ratio of the # side lengths. For this first problem, to make it trivial, we will do that # and divide delx or delz out of the problem. # discretize the domain lengthx = 1500 # width lengthz = 1000 # depth nodes_x = 35 # in 2D, you rapidly run out of memory!, try 200 nodes!! delx = lengthx/(nodes_x-1) pdelx = lengthx/(nodes_x) # note this is for plotting only!! # the solution to the memory problem is to use SPARSE arrays in scipy nodes_z = int(nodes_x * (lengthz/lengthx)) # BCs BCtop = 1 BCbottom = 7 BCright = 3 BCleft = 6 nodes = nodes_z*nodes_x # nodes is the total length of the Temperature field vector # including all the unknowns and the BC nodes # Our mesh is a grid, our numbering scheme is that the nodes are numbered # from 0 starting in the top left corner. Numbers increase downward # (column-wise) since the problem is wider than it is deep. In large # problems it is best to number in the shortest dimension of the grid. # Therefore the left hand side numbers from 0 to nodes_z - 1 # all BCs at boundaries are fixed # steady state and constant K # values of the diagonal and off diagonal coeff (but in 2D these appear in # 'bands' in the matrix, not merely on the main diagonal alpha = -4 # value of main diagonal beta = 1 # value of both off diagonals # better in large problems would be to use sparse matrices below # for 2D or 3D problems, most of the work is assembling the A matrix # It is necessary that you have a clear picture of the indexing of your nodes! # construct the non-trivial bands in the matrix, using a 5 point FD kernel centered on each node:- # beta*node to left, zeros, beta*node above, alpha*node, beta*node below, zeros, beta*node to right # You can assemble the A matrix using fancy 'slicing', but for this first time we will use a very # explicit assembly by looping over all the nodes to check if they are boundary nodes or not. # make an A matrix with '1's on the diagonal, zeros elsewhere A = np.identity(nodes) for i in np.arange(nodes_z,nodes-nodes_z): #this loops over all the nodes, # except the left and right hand sides (we are counting column-wise) # the ith row of the A matrix is the equation for the ith node # check if we are at a boundary node, if not then do: if (((i+1)%nodes_z) !=0 ) and ((i%nodes_z) !=0): #leave BC nodes as '1's #boundary bottom nodes are divisible by nodes_z-1, while top nodes are 1 more A[i,i]=alpha A[i,i-1]=beta A[i,i+1]=beta A[i,i-nodes_z]=beta #the node to the left of a node is i-nodes_z A[i,i+nodes_z]=beta # make the known vector, with BC's in correct locations C = np.zeros(nodes) # the following makes good use of python 'slicing' # (see if you can figure out what is being done!) C[:nodes_z] =BCleft # C[nodes_z:nodes-nodes_z-1:nodes_z] =BCtop C[2*nodes_z-1::nodes_z] =BCbottom #bottom BC C[nodes-nodes_z:] =BCright #C[-1] =BCbottom #improves symetry of plot (nothing else) # plot a SPY diagram, very useful to find out if your indexing is incorect! fig1 = plt.figure() ax1 = fig1.add_subplot(111) ax1.spy(A) ax1.set_title("SPY plot of A Array") plt.show() #solve Tc = np.linalg.solve(A,C) # direct solver in numpy, Tc now contains the solution vector of temperatures # we improve on the simple plotting below, showing several plotting possibilities fig2 = plt.figure() ax2 = fig2.add_subplot(1,1,1) # Tc is a temperature vector, we need to reshape it into an array with # nodes_z rows and nodes_x cols # NOTE in an array the 1st index is the row number, the 2nd the col (in 3D the 3rd is the col, # the 2nd is the row, and the 1st is the plane) # Thinking in 2 and 3D you need to be very careful, so below we want to make an array # since we numbered down the rows as our index, we are numbering in the WRONG direction B = Tc.reshape(nodes_x,nodes_z) # reshaping is very useful for plotting, # this produces an array with the col numbers going down and the rows across, in other words # it is the transpose of the array we want Bt = B.T # arrays have a '.T' method to transpose # pcolor plots a psuedo color image of the array, using the array values as gradations in color # you can give it x,y position arrays, or more easily give it x,y array eadge position vectors and # it will figure out the x,y coords (this is called braodcasting) # A subtlty is that the color is a patch to the right and down of the nodes, therefore yoiu need to # plot an array that is 1 larger in each direction to see the full psuedo plot f=ax2.pcolor(np.arange(nodes_x+1)*pdelx,np.arange(nodes_z+1)*pdelx,Bt,cmap=cm.magma,shading='auto') # the '.T' is the transpose # as a result of the expanded plot, the psuedo color image is too large, but we partially fix with # 'pdelx' , but if you look closely the interior nodes are not quite in the correct places #ax2.grid() ax2.invert_yaxis() # always good in Geoplots to make the surface at the top bar = plt.colorbar(f) # colorbar places a scale to the colors of the psuedo plot bar.ax.set_ylabel("Temperature in C") ax2.set_ylabel("Depth") ax2.set_xlabel("Width in Meters") ax2.set_title("Temperature field") # demonstrate a 3D surface plot. Use meshgrid to make x and z location arrays # surface plots do not have the offset position problems of pcolor plots, the B values # colored at the grid node points, but we need to give it the x,z coords xx,zz = np.meshgrid(np.arange(nodes_x)*delx,np.arange(nodes_z)*delx) #produces a complete set of x,z coords fig3 = plt.figure() ax3 = fig3.add_subplot(111, projection='3d') # a surface plot is a 3D projection ax3.plot_surface(xx,zz,Bt,cmap=cm.magma, rstride = 1, cstride = 1,edgecolor='k') #note using a 'colormap' ax3.set_ylabel("Depth") ax3.set_xlabel("Width in Meters") ax3.set_title("Temperature field as a 3D surface") # This is another useful figure, a contour plot. Here we make the contours stand out as black fig4 = plt.figure() ax4 = fig4.add_subplot(111) cs = ax4.contour(xx,zz,Bt,colors ='k') # contour just plots the contours ax4.contourf(xx,zz,Bt) # contourf fills the contours with colors ax4.clabel(cs,inline=1,fontsize=12) # clabel puts labels on the contours ax4.invert_yaxis() ax4.set_ylabel("Depth") ax4.set_xlabel("Width in Meters") ax4.set_title("Temperature field as a contour") # finally a plot showing the directions of heat flow, overlayed over the psuedo color plot fig5 = plt.figure() ax5 = fig5.add_subplot(111) # gradient produces the difference arrays in the 2 axis directions U,V = np.gradient(Bt) # this gives the positive gradient in x,z, returns 2 arrays with these gradients ax5.invert_yaxis() ax5.pcolor(np.arange(nodes_x+1)*pdelx,np.arange(nodes_z+1)*pdelx,Bt) # We want the heat to flow from hot to cold, so we need negative, also again the directions are reversed ax5.quiver(xx,zz,-V,U) ax5.set_ylabel("Depth") ax5.set_xlabel("Width in Meters") ax5.set_title("Temperature field with heat flow directions and Gradients")