# Solution of the Laplace Equation on a Polar grid

In this workbook we are going to look at the solution of the Laplace equations on a polar coordinate grid.  In this case we are going to use the Bi-Conjugate Gradient solver to solve the problem.  It should be noted that there is a periodic boundary condition at $\theta=0$ and $\theta=2\pi$ so the coeficient matrix needs to reflecxt this.

$$
\frac{\partial^2 u}{\partial r^2} + \frac{1}{r}\frac{\partial u}{\partial r} + \frac{1}{r^2}\frac{\partial^2 u}{\partial \theta^2} = 0
$$

The boundary conditions on the inside and outside of the pipe are given by
$$\begin{align}
    u(1, \theta) &= 0, \quad 0\le\theta\le 2\pi \label{eq:inner_BC}\\
    \frac{\partial u}{\partial r}(3, \theta) &= \begin{cases}
        0, \quad & 0\le\theta<\pi \\
        1, \quad & \pi\le\theta\le 2\pi
        \end{cases}
\end{align} $$


```python
import numpy as np
import scipy.sparse.linalg as LA
import scipy.sparse as sps
import matplotlib.pyplot as plt
import time

class Grid:
    '''Class defining a 2D computational grid for
    a polar cartesian grid.  We assume that theta
    goes from 0 to 2\pi and has periodic boundaires
    in the theta direction.'''
    
    def __init__(self,ni,nj):
        # set up information about the grid
        self.origin = 0.0 # unit circle  
        self.extent = 1.0
        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.r = np.zeros((nj, ni))
        self.theta = np.zeros((nj, ni))

    def set_origin(self, r0):
        self.origin = r0
    
    def set_extent(self,r1):
        self.extent = r1
        
    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.'''
        
        theta = np.linspace(0.0, 2*np.pi, self.Ni+1)
        theta_ord = np.delete(theta,[-1])
        r_ord = np.linspace(self.origin, self.extent, self.Nj)
        self.theta,self.r = np.meshgrid(theta_ord, r_ord)
        if not Quiet:
            print(self)

    def Delta_r(self):
        # calculate delta x
        return self.r[1,0]-self.r[0,0]
    
    def Delta_theta(self):
        # calculate delta y
        return self.theta[0,1]-self.theta[0,0]
    
    def find(self,point):
        '''find the i and j ordinates of the grid cell which contains 
        the point (theta,r).  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_theta()
        grid_y = (point[1] - self.origin[1])/self.Delta_r()
        return int(grid_x), int(grid_y)
    
    def plot(self,title):
        '''produce a contour plot of u(theta,r)'''
        # convert polar to Cartesian coordinates
        x = self.r*np.cos(self.theta)
        y = self.r*np.sin(self.theta)
        
        # produce the plot
        fig, ax = plt.subplots()
        cmap = plt.get_cmap('jet')
        ax.axis('equal')
        ax.set_aspect('equal', 'box')
        cf = ax.contourf(x,y,self.u,cmap=cmap,levels = 21)
        fig.colorbar(cf, ax=ax)
        ax.set_title(f'{title} ({self.Ni} x {self.Nj} grid)')
        return plt
    
    def __str__(self):
        # describe the object when asked to print it
        return 'Uniform {}x{} polar grid from r = {} to {}.'.format(self.Ni, self.Nj, self.origin, self.extent)
```


```python
#Set up the coursework 
def coursework_one(ni,nj):
    # set up a mesh
    mesh = Grid(ni,nj)
    mesh.set_origin(1.0)
    mesh.set_extent(3.0)
    mesh.generate()
    return mesh
```


```python
# lets set up the test problem
test = coursework_one(41,41)
print(test)
```

Now we can implement the solver.  We can use the solver from the BiCGStab workbook as a template, remembering that in this case we know that the r=1 boundary is a Dirichlet BC with $u=(1,\theta)=0$.  The $r=3$ boundary condition is a Neumann condition where

$$\begin{align}
    u(1, \theta) &= 0, \quad 0\le\theta\le 2\pi \label{eq:inner_BC}\\
    \frac{\partial u}{\partial r}(3, \theta) &= \begin{cases}
        0, \quad & 0\le\theta<\pi \\
        1, \quad & \pi\le\theta\le 2\pi
        \end{cases}
\end{align} $$

We also need to enfroce the periodic boundary condition $u(r,0)=u(r,2\pi)$ and to modify the coefficients in the coefficient matrix to include both the additional derivative and the $1/r$ and $1/r^2$ coefficients.


```python
def PolarLaplaceSolver(grid,tol=0.5e-7):
    '''Solve the two dimensional laplace equation on a polar coordinates
    using the bi-conjugate gradient stabilised matrix solver (BiCGStab).  
    This function assembles the coeficient matrix A and the RHS vector b
    taking account of the boundary conditions specified in the question.
    It should call the  BiCGStab solver from scipy.  The value tol is p
    assed to BiCGStab routine The solution vector x is then unpacked 
    into grid.u. It returns the info value from BiCGStab if this is 
    zero everything worked.'''
    
    # Create the A matrix using the lil format and the b vector
    # as numpy vector.
    N = (grid.Nj-2)*(grid.Ni-2)
    A_mat = sps.lil_matrix((N, N), dtype=np.float64)
    b_vec = np.zeros(N, dtype=np.float64)

    # calculate the coefficients
    beta = grid.Delta_x()/grid.Delta_y()
    beta_sq = beta**2
    R_x = - 1/(2*(1+beta_sq))
    R_y = beta_sq * R_x

    # -----
    # your code goes here
    # -----
    
    # call bicgstab
    x_vec, info = LA.bicgstab(A_mat,b_vec,tol=tol)
    
    if info==0:
        # unpack x_vec into u
        for j in range(1, grid.Nj-1):
            for i in range(1, grid.Ni-1):
                k = (i-1) + (grid.Ni-2)*(j-1)
                grid.u[j,i]=x_vec[k]
    
    return info
```


```python
# Test the solution on the 40x40 grid
info = PolarLaplaceSolver(test)
if info==0:
    plt = test.plot('Coursework 1')
    plt.show()
else:
    print('Error code ',info,' returned by BiCGStab')
```

719e24d4c1b7f0d4b4bba6968fc92feb060c582a
ipynb
Jupyter Notebook
coursework/Laplace-polar.ipynb
SolarisRe/PDE3
5dbda654487f1cadcb02e72162fd319160ab18a7
CC0-1.0
2022-01-19T15:08:44.000Z
2022-01-27T15:42:34.000Z
coursework/Laplace-polar.ipynb
SolarisRe/PDE3
5dbda654487f1cadcb02e72162fd319160ab18a7
CC0-1.0
coursework/Laplace-polar.ipynb
SolarisRe/PDE3
5dbda654487f1cadcb02e72162fd319160ab18a7
CC0-1.0
2022-01-24T07:03:07.000Z
2022-03-23T12:06:03.000Z
Qwen/Qwen-72B
1. YES
2. YES