<jupyter_start><jupyter_text>Week 4 worksheet: Matrix methods for the solution of Elliptic equations.This worksheet contains a number of exercises covering only the numerical aspects of the course. Some parts, however, still require you to solve the problem by hand, i.e. with pen and paper. The rest needs you to write pythob code. It should usually be obvious which parts require which. Suggested readingYou will see lists of links to further reading and resources throughout the worksheets, in sections titled **Learn more:**. These will include links to the Python documentation on the topic at hand, or links to relevant book sections or other online resources. Unless explicitly indicated, these are not mandatory reading, although of course we strongly recommend that you consult them! Displaying solutionsSolutions will be released after the workshop, as a new `.txt` file in the same GitHub repository. After pulling the file to Noteable, **run the following cell** to create clickable buttons under each exercise, which will allow you to reveal the solutions. Note:This workbook expects to find a diretory called figures in the same folder as well as the scripts folder.  Please make sure you download figures (and the files it contains) from the GitHub.<jupyter_code>%run scripts/create_widgets.py W04<jupyter_output><empty_output><jupyter_text>*How it works: You will see cells located below each exercise, each containing a command starting with `%run scripts/show_solutions.py`. You don't need to run those yourself; the command above runs a script which automatically runs these specific cells for you. The commands in each of these cells each create the button for the corresponding exercise. The Python code to achieve this is contained in `scripts/show_solutions.py`, and relies on [IPython widgets](https://ipywidgets.readthedocs.io/en/latest/examples/Widget%20Basics.html) --- feel free to take a look at the code if you are curious.*<jupyter_code>%%javascript
MathJax.Hub.Config({
    TeX: { equationNumbers: { autoNumber: "AMS" } }
});<jupyter_output><empty_output><jupyter_text>Exercise 1$$\newcommand{\vect}[1]{\bm 1}\newcommand{\grad}{\nabla}\newcommand{\pderiv}[2]{\frac{\partial 1}{\partial 2}}\newcommand{\pdderiv}[2]{\frac{\partial^2 1}{\partial 2^2}}$$In this laboratory we will solve$$0=T_{xx}+T_{yy}$$subject to the boundary conditions$$\begin{align*}T(x,0)&=0&0.0\le x\le 0.2 \\T(0,y)&=0&0.0\le y \le 0.2\\T(x,0.2)&=250x&0.0\le x\le 0.2 \\T(0.2,y)&=250y&0.0\le y\le 0.2.\end{align*}$$using the stabilised bi-conjugate gradient solver applied to a sparse matrix representing the coefficients of the linear system.  As explained in the lecture we are starting with the 2nd order, central finite difference approximation of the 2D Laplacian.$$T_{ij}=\frac{T_{i-1j}+T_{i+1j}+\beta^2\left(T_{ij-1}+T_{ij-1}\right)}{2(1+\beta^2)}$$where $\beta=\frac{\Delta x}{\Delta y}$. Part a)Re-arrange the finite difference approximation into the form $$r_y T_{ij-1} + r_x T_{i-1j} +  T_{ij}  + r_x T_{i+1j}  + r_y T_{ij+1} = 0$$and hence find the coefficients $r_x$ and $r_y$.<jupyter_code>%run scripts/show_solutions.py W04_ex1_parta<jupyter_output><empty_output><jupyter_text>Part b)The following diagram shows a uniformly spaced, $7\times7$, grid of $T$ values with boundary conditions stored in outer most entries.  Following Python notation the bottom left T value is T[0,0] and the top right T value is  T[6,6].   The grid also shows the computational molecule for the Laplacian operator and the cells in which the top & bottom and left & right boundary values must be applied.![](figures/LaplaceGrid.png)Starting with the first interior node (shown in white) write in the $k$ index values for each interior node, then write down the RHS values $b_k$ associated with the interior points and the points at which the boundary conditions must be applied.<jupyter_code>%run scripts/show_solutions.py W04_ex1_partb<jupyter_output><empty_output><jupyter_text>Part C)In this test we are **not** using the boundary conditions for our test problem but assuming that $T=1$ on all boundaries.On the following grid fill in the non-zero elements in the coefficient matrix, $A$, and the RHS vector $\vec b$.  You should use different colours or symbols to show entries in the coefficient matrix which are  $r_x$, $r_x$ and $1$ respectively.  As in the previous question _Python_ indexing has been used so the $25\times25$ matrix has elements $a_{ij}$ where $i,j = 0, 1, 2, 3, 5 \cdots 22, 23, 24.$ ![](figures/SparsityPattern.png)<jupyter_code>%run scripts/show_solutions.py W04_ex1_partc<jupyter_output><empty_output><jupyter_text>Part d)Using the Grid class defined below, write a program which assembles both the coefficient matrix and the RHS vector.  Produce a spy plot for the coefficient matrix and check it against the expected pattern (from part c)<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import scipy.sparse as sps
import time

class Grid:
    '''Class defining a 2D computational grid.  The grid object
    contains is a regular cartesian grid with a single variable, u.
    It stores information about the number of grid points in the i
    and j directions, the ordinates of these points and the bottom
    left corner of the gird (the origin) and the top right corner
    (the extent).
    
    Written by Prof <NAME>, School of Engineering
    (c) 2021 The University of Edinburgh
    Licensed under CC-BY-NC.'''
    
    DIRICHLET_BC = 0
    NEUMANN_BC = 1

    BC_NAME = ['left', 'right', 'top', 'bottom']
    
    def __init__(self,ni,nj):
        # set up information about the grid
        self.origin = (0.0, 0.0)  # bottom left
        self.extent = (1.0, 1.0)  # top right
        self.Ni = ni # grid points in i direction
        self.Nj = nj # grid points in j direction
        
        # initialse x,y and u arrays
        self.u = np.zeros((nj, ni))
        self.x = np.zeros((nj, ni))
        self.y = np.zeros((nj, ni))
        
        # boundary conditions (left right top and bottom)
        self.BC = [self.DIRICHLET_BC, self.DIRICHLET_BC, 
                   self.DIRICHLET_BC, self.DIRICHLET_BC]

    def set_origin(self,x0,y0):
        self.origin = (x0, y0)
    
    def set_extent(self,x1,y1):
        self.extent = (x1, y1)
        
    def generate(self,Quiet=True):
        '''generate a uniformly spaced grid covering the domain from the
        origin to the extent.  We are going to do this using linspace from
        numpy to create lists of x and y ordinates and then the meshgrid
        function to turn these into 2D arrays of grid point ordinates.'''
        x_ord = np.linspace(self.origin[0], self.origin[0] + self.extent[0], self.Ni)
        y_ord = np.linspace(self.origin[1], self.origin[1] + self.extent[1], self.Nj)
        self.x, self.y = np.meshgrid(x_ord,y_ord)
        if not Quiet:
            print(self)

    def Delta_x(self):
        # calculate delta x
        return self.x[0,1]-self.x[0,0]
    
    def Delta_y(self):
        # calculate delta y
        return self.y[1,0]-self.y[0,0]
    
    def find(self,point):
        '''find the i and j ordinates of the grid cell which contains 
        the point (x,y).  To do this we calculate the distance from
        the point to the origin in the x and y directions and then
        divide this by delta x and delta y.  The resulting real ordinates
        are converted to indices using the int() function.'''
        grid_x = (point[0] - self.origin[0])/self.Delta_x()
        grid_y = (point[1] - self.origin[1])/self.Delta_y()
        return int(grid_x), int(grid_y)
    
    def set_Neumann_bc(self,side):
        try:
             self.BC[self.BC_NAME.index(side)] = self.NEUMANN_BC
        except:
             print('error {} must be one of {}'.format(side,self.BC_NAME))
     
    def set_Dirichlet_bc(self,side):
        try:
            self.BC[self.BC_NAME.index(side)] = self.DIRICHLET_BC
        except:
            print('error {} must be one of {}'.format(side,self.BC_NAME))
    
    def report_BC(self):
        '''compile a string listing the boundary conditions on each side.
        We build up a string of four {side name}: {BC type} pairs and
        return it'''
        
        # initialise the string
        string = ''
        
        # loop over the sides
        for side in range(4):
            # add the side name
            string = string + self.BC_NAME[side]
            # and the bounday condition type
            if self.BC[side] == self.DIRICHLET_BC:
                string = string + ': Dirichlet, '
            elif self.BC[side] == self.NEUMANN_BC:
                string = string + ': Neumann, '
        return string[:-2] +'.' # lose the last comma and space.
    
    def __str__(self):
        # describe the object when asked to print it
        describe = 'Uniform {}x{} grid from {} to {}.'.format(self.Ni, self.Nj, self.origin, self.extent)
        boundaries = self.report_BC()
        return describe + '\nBoundaries conditions are - ' + boundaries
def assemble_matrix(mesh):
    '''assemble the matrix and the RHS vector'''
    

    
    return A_mat, b_vec
%run scripts/show_solutions.py W04_ex1_partd<jupyter_output><empty_output><jupyter_text>Part e)Use the stabilised bi-conjugate gradient solver to solve the system $A\vec x=\vec b$.  Use matrix algebra to check that $A\vec x-\vec b\approx 0$, and check that $T_{ij}=1$ for all internal points.You can use the _numpy_ function allclose() to check if two vectors or arrays are similar to each other.  You may also find the numpy functions zeros_like() and ones_like() which return the arrays of 0 and 1 the same shape as the array passed to them useful. Note:This is one of the famous _"Nothing should happen"_ tests which are extreemly useful in checking code. We know that the solution must be that T=1 everywhere but this will only happen if the coefficients are set up properly.  If the test fails its quite easy to find the errors as we know the values we should be calculating (and the grid is small).<jupyter_code>import scipy.sparse.linalg as LA



%run scripts/show_solutions.py W04_ex1_parte<jupyter_output><empty_output><jupyter_text>Part f)Finally, set up and solve the test problem.  You need to write a function unpack(mesh, x_vec) which unpacks the x vector into the u values of mesh.then plot the solution.<jupyter_code>%run scripts/show_solutions.py W04_ex1_partf<jupyter_output><empty_output><jupyter_text>Part g)Conduct a mesh refinement study using a sequence of 5 successively refine grids, starting with a $21\times 21$ grid.  Conduct a mesh refinement study using the integral quantity $Q=\int_{0.15}^{0.18} u(x,0.18,)\ dx$, calculated on each grid.   Use the following integrate function to calculate Q.  You should also use the ILU factorisation to accelerate the BiCGStab solver.ilu = LA.spilu(A.tocsc(), drop_tol=1e-6, fill_factor=100)M_mat = LA.LinearOperator(A_mat.shape, ilu.solve)x_vec, info = LA.bicgstab(A,b,atol=0.5e-12,M=M_mat)Use the refinement_analysis toolkit to perform a refinement analysis, comment on your findings.<jupyter_code>import scipy.integrate as integrate

def integrate_u_dx(mesh,x0,x1,y):
    '''Calculate U=\int_{1.0}^{1.95)u(x,0.25) dx  using the
    u value stored on the grid and simpsons rule'''
    
    # find the left and right grid points
    i0,j = mesh.find((x0,y))
    i1,j = mesh.find((x1,y))

    # add 1 to i1 as we are going to use it as an upper 
    # bound forarray slicing
    i1 = i1 +1
    
    # integrate
    return integrate.simps(mesh.u[j,i0:i1],mesh.x[j,i0:i1])

# run cases om the 321x321, 161x161, 81x81, 41x41 and 21x21 grids
U_val = []
dx_val = []
run_time = []
n_pts =[]

from refinement_analysis import refinement_analysis
# lets do the refinement analysis

%run scripts/show_solutions.py W04_ex1_partg<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09341861639454638, 'fraction_numerical': 0.01893620602592156, 'mean_word_length': 3.5062571103526734, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 38, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '310819', 'license': 'CC0-1.0', 'n_tokens_mistral': 3761, 'n_tokens_neox': 3568, 'n_words': 1400, 'path': 'laboratories/W04_PDE3_Matrix_methods_elliptic_PDEs.ipynb', 'repo_name': 'SolarisRe/PDE3', 'size': '11888'}