# Examples on Using PyModPDE
---
__Mokbel Karam, James C. Sutherland and Tony Saad__

__Department of Chemical Engineering, University of Utah__

Below are the necessary libraries used in the examples below.


```python
import matplotlib.pyplot as plt
plt.rc("font", family='serif')
import numpy as np
from sympy import symbols,latex
from matplotlib import cm, animation
from IPython.display import Math
```

##  Example Using FTUS on Advection PDE (1D)
---

Consider the following PDE
\begin{equation}
u_t=-a u_x ; a>0
\end{equation}
to be solved using forward Euler discretization in time and upwind discretization in space(FTUS). The difference equation is
\begin{equation}
\frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t}=-a\left(\frac{u_{i}^{n}-u_{i-1}^{n}}{\Delta x}\right).
\end{equation}


```python
from src.pymodpde import DifferentialEquation, i, k, n # import the differential equation and the indicies we need to use

a= symbols('a') # define positive a constant for the advection velocity

DE1 = DifferentialEquation(dependentVarName='u',independentVarsNames=['x']) # construct a one D differential equation with 'x' and 'u' as independent and dependent variables, respectively

# method I of constructing the rhs:
# construct the first upwind derivative of the dependent variabe 'u' 
# with respect to the independent variable 'x' using the stencil -1, and 0 evaluated at time tn = n
advectionTerm1 = DE1.expr(order=1, directionName='x', time=n, stencil=[-1, 0])
                                                                        

# set the rhs of the differential equation
DE1.set_rhs(- a * advectionTerm1 )

# compute the modified equation up to two terms

DE1.generate_modified_equation(nterms=4)
```


```python
# display the mofified equation
DE1.display_modified_equation()
```


$$\frac{\partial u }{\partial t}  =  - a \frac{\partial u }{\partial x}  + \frac{a}{2} \left(- \Delta{t} a + \Delta{x}\right) \frac{\partial^{2} u }{\partial x^{2}}   + a \left(- \frac{\Delta{t}^{2} a^{2}}{3} + \frac{\Delta{t} \Delta{x}}{2} a - \frac{\Delta{x}^{2}}{6}\right) \frac{\partial^{3} u }{\partial x^{3}}   + a \left(- \frac{\Delta{t}^{3} a^{3}}{4} + \frac{\Delta{t}^{2} \Delta{x}}{2} a^{2} - \frac{7 \Delta{t}}{24} \Delta{x}^{2} a + \frac{\Delta{x}^{3}}{24}\right) \frac{\partial^{4} u }{\partial x^{4}}  $$



```python
# display the amplification factor
DE1.display_amp_factor()
```


$$e^{\Delta{t} \alpha} = - \frac{\Delta{t} a}{\Delta{x}} + \frac{\Delta{t} a}{\Delta{x}} e^{- i \Delta{x} k_{1}} + 1.0$$



```python
# print the latex form of the generate_modified_equation
print(DE1.latex_modified_equation())
```

    \frac{\partial u }{\partial t}  =  - a \frac{\partial u }{\partial x}  + \frac{a}{2} \left(- \Delta{t} a + \Delta{x}\right) \frac{\partial^{2} u }{\partial x^{2}}   + a \left(- \frac{\Delta{t}^{2} a^{2}}{3} + \frac{\Delta{t} \Delta{x}}{2} a - \frac{\Delta{x}^{2}}{6}\right) \frac{\partial^{3} u }{\partial x^{3}}   + a \left(- \frac{\Delta{t}^{3} a^{3}}{4} + \frac{\Delta{t}^{2} \Delta{x}}{2} a^{2} - \frac{7 \Delta{t}}{24} \Delta{x}^{2} a + \frac{\Delta{x}^{3}}{24}\right) \frac{\partial^{4} u }{\partial x^{4}}  



```python
# print the latex form of the amplification factor
print(DE1.latex_amp_factor())
```

    e^{\Delta{t} \alpha} = - \frac{\Delta{t} a}{\Delta{x}} + \frac{\Delta{t} a}{\Delta{x}} e^{- i \Delta{x} k_{1}} + 1.0


Using the second method for constructing the rhs


```python
from src.pymodpde import DifferentialEquation, j, n # import the differential equation and the indicies we need to use

b = symbols('b') # define positive a constant for the advection velocity

DE2 = DifferentialEquation(dependentVarName='f',independentVarsNames=['y'], indices=[j])# construct the first upwind derivative of the dependent variabe 'f' 
                                                                        # with respect to the independent variable 'y' using the stencil -1, and 0 evaluated at time tn = n
    

# method II of constructing the rhs:
advectionTerm2 = (DE2.f(time=n, y=j) - DE2.f(time=n, y = j-1)) / DE2.dy 
    
DE2.set_rhs(- b * advectionTerm2 ) # set the rhs of the Differential equation DE

DE2.generate_amp_factor() # computing the amplification factor
```


```python
# display the amplification factor
DE2.display_amp_factor()
```


$$e^{\Delta{t} \alpha} = - \frac{\Delta{t} b}{\Delta{y}} + \frac{\Delta{t} b}{\Delta{y}} e^{- i \Delta{y} k_{1}} + 1.0$$


The difference equation for Backward in time and Upwind in space is
\begin{equation}
\frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t}=-a\left(\frac{u_{i}^{n+1}-u_{i-1}^{n+1}}{\Delta x}\right).
\end{equation}


```python
from src.pymodpde import DifferentialEquation
from sympy import *

# define the advection velocity
a= symbols("a")
# define the spatial and temporal indices
i, n = symbols("i n")

# construct a time dependent differential equation
DE = DifferentialEquation(dependentVarName="u",independentVarsNames=["x"], indices=[i], timeIndex=n)

# method I of constructing the rhs:
advectionTerm = -a * DE.expr(order=1, directionName="x", time=n+1,stencil=[-1, 0])

# setting the rhs of the differential equation
DE.set_rhs( advectionTerm )

# computing and displaying the modified equation up to two terms
DE.generate_modified_equation(nterms=2)
```


```python
# display the mofified equation
DE.display_modified_equation()
```


$$\frac{\partial u }{\partial t}  =  - a \frac{\partial u }{\partial x}  + \frac{a}{2} \left(\Delta{t} a + \Delta{x}\right) \frac{\partial^{2} u }{\partial x^{2}}  $$



```python
# display the amplification factor
DE.display_amp_factor()
```


$$e^{\Delta{t} \alpha} = \frac{\Delta{x} e^{i \Delta{x} k_{1}}}{\Delta{t} a e^{i \Delta{x} k_{1}} - \Delta{t} a + \Delta{x} e^{i \Delta{x} k_{1}}}$$


## Example Using BTCS on Advection-Diffusion PDEs (1D)
---

Consider the following PDE
\begin{equation}
u_t=-cu_x+\gamma u_{xx}
\end{equation} 
to be solved using backward Euler discretization in time and central discretization in space(BTCS). The difference equation is
$$
\frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t}=-c\left(\frac{u_{i+1}^{n+1}-u_{i-1}^{n+1}}{2\Delta x}\right)+ \gamma \left(\frac{u_{i+1}^{n+1}-2u_i^{n+1}+u_{i-1}^{n+1}}{\Delta x^2}\right).
$$


```python
from src.pymodpde import DifferentialEquation, i, n # import the differential equation and the indicies we need to use

c, g = symbols('c gamma') # define positive a constant for the advection velocity

DE3 = DifferentialEquation(dependentVarName='u',independentVarsNames=['x'], indices=[i])# construct the first upwind derivative of the dependent variabe 'f' 
                                                                        # with respect to the independent variable 'y' using the stencil -1, and 0 evaluated at time tn = n
    

# method II of constructing the rhs:
advectionTerm3 = (DE3.u(time=n+1, x=i+1) - DE3.u(time=n+1, x=i-1)) / DE3.dx /2

diffusionTerm3 = (DE3.u(time=n+1, x=i+1) - 2* DE3.u(time=n+1, x=i) + DE3.u(time=n+1, x=i-1)) / DE3.dx / DE3.dx
    
DE3.set_rhs(- c * advectionTerm3  ) # set the rhs of the Differential equation DE

DE3.generate_modified_equation(nterms=3)  # compute the modified equation up to 2 terms and display the latex form of the modified equation
```


```python
# display the mofified equation
DE3.display_modified_equation()
```


$$\frac{\partial u }{\partial t}  =  - c \frac{\partial u }{\partial x}  + \frac{\Delta{t} c^{2}}{2} \frac{\partial^{2} u }{\partial x^{2}}   - c \left(\frac{\Delta{t}^{2} c^{2}}{3} + \frac{\Delta{x}^{2}}{6}\right) \frac{\partial^{3} u }{\partial x^{3}}  $$



```python
# display the amplification factor
DE3.display_amp_factor()
```


$$e^{\Delta{t} \alpha} = \frac{2.0 \Delta{x} e^{i \Delta{x} k_{1}}}{\Delta{t} c e^{2.0 i \Delta{x} k_{1}} - \Delta{t} c + 2.0 \Delta{x} e^{i \Delta{x} k_{1}}}$$


## Example Using FTUS on Advection PDEs (2D)
---

Consider the advection equation in a two dimensional space
\begin{equation}
u_{t}=-bu_{x}-cu_{y};\quad(b,c)\in\mathbb{R}^{+}
\end{equation}
to be discretized using a forward Euler method in time and UPWIND in space (FTUS)
\begin{equation}
\frac{u_{i,j}^{n+1}-u_{i,j}^{n}}{\Delta t}=-b\frac{u_{i,j}^{n}-u_{i-1,j}^{n}}{\Delta x}-c\frac{u_{i,j}^{n}-u_{i,j-1}^{n}}{\Delta y},
\end{equation}

with a modified equation 


```python
from src.pymodpde import DifferentialEquation, i, j, n # import the differential equation and the indicies we need to use

b, c= symbols('b c') # define positive a constant for the advection velocity

DE1 = DifferentialEquation(dependentVarName='u',independentVarsNames=['x','y']) # construct a one D differential equation with 'x' and 'u' as independent and dependent variables, respectively

# method II of constructing the rhs:
advectionTerm1 = (DE1.u(time=n, x=i, y=j) - DE1.u(time=n, x=i-1, y=j))/DE1.dx

advectionTerm2 = (DE1.u(time=n, x=i, y=j) - DE1.u(time=n, x=i, y=j-1))/DE1.dy

# set the rhs of the differential equation
DE1.set_rhs( -b * advectionTerm1 - c * advectionTerm2 )

# compute and the modified equation up to two terms
DE1.generate_modified_equation(nterms=2)
```


```python
# display the mofified equation
DE1.display_modified_equation()
```


$$\frac{\partial u }{\partial t}  =  - c \frac{\partial u }{\partial y}  - b \frac{\partial u }{\partial x}  + \frac{c}{2} \left(- \Delta{t} c + \Delta{y}\right) \frac{\partial^{2} u }{\partial y^{2}}   - \Delta{t} b c \frac{\partial^{2} u }{\partial x\partial y}   + \frac{b}{2} \left(- \Delta{t} b + \Delta{x}\right) \frac{\partial^{2} u }{\partial x^{2}}  $$



```python
# display the amplification factor
DE1.display_amp_factor()
```


$$e^{\Delta{t} \alpha} = - \frac{\Delta{t} c}{\Delta{y}} + \frac{\Delta{t} c}{\Delta{y}} e^{- i \Delta{y} k_{2}} - \frac{\Delta{t} b}{\Delta{x}} + \frac{\Delta{t} b}{\Delta{x}} e^{- i \Delta{x} k_{1}} + 1.0$$


### Numerical Solution


```python
nx = 64
ny = 64
Lx = 1.0
Ly = 1.0
x = np.linspace(0,Lx,nx)
y = np.linspace(0,Ly,ny)
dx = Lx/(nx-1)
dy = Ly/(ny-1)
# create a grid of coordinates
xx,yy = np.meshgrid(x,y)

u0 = lambda x,y : np.exp(-(x-0.2)**2/0.01)*np.exp(-(y-0.2)**2/0.01)

cx = 1.0
cy = 1.0

dt = 1/128
tend = 0.5 #s
t = 0

cflx = cx * dt/dx
cfly = cy * dt/dy

# setup the initial condition
sol = []
u = np.zeros([ny+2,nx+2]) # we will ghost cells to simplify the implementation of periodic BCs
u[1:-1, 1:-1] = u0(xx,yy)
# set periodic boundaries
u[:,0] = u[:,-3] #x-minus face
u[:,-1] = u[:,2] #x-plus face
u[0,:] = u[-3,:] #y-minus face
u[-1,:] = u[2,:] #y-plus face
sol.append(u)
```


```python
while t < tend:
    un = sol[-1]
    unew = un.copy()
    unew[1:-1,1:-1] = un[1:-1,1:-1] - cflx * (un[1:-1,1:-1] - un[1:-1,:-2]) - cfly * (un[1:-1,1:-1] - un[:-2,1:-1])
    # set periodic boundaries
    unew[:,0] = unew[:,-3] #x-minus face
    unew[:,-1] = unew[:,2] #x-plus face
    unew[0,:] = unew[-3,:] #y-minus face
    unew[-1,:] = unew[2,:] #y-plus face

    sol.append(unew)
    t += dt
```


```python
levs = np.linspace(0,1,15)
plt.figure(figsize=(10,8))
plt.contourf(xx,yy,sol[-1][1:-1,1:-1]+sol[0][1:-1,1:-1],cmap=cm.viridis,levels=levs,vmax=1.0,vmin=0)
plt.text(0.1, 0.4, 'time=0.0s', fontsize=12,family='serif')
plt.text(0.6, 0.92,'time={}s'.format(tend), fontsize=12,family='serif')
plt.xlabel('x', fontsize=14,family='serif')
plt.ylabel('y', fontsize=14,family='serif')
cbar = plt.colorbar()
plt.clim(-1,1)
cbar.set_ticks(np.linspace(-1,1,5))
plt.savefig('./advection_2d.pdf',transparent=True)
plt.show()
```

## Numerical Solver
---


```python
%matplotlib notebook
from numerical.layout import *
plty=PlotStyling()
```


    Tab(children=(VBox(children=(HBox(children=(FloatText(value=0.01, description='dt:'), FloatText(value=10.0, de…



```python

```

9e509c7e2cc1e2acfcf2c0cbc5bee6013fd5a394
ipynb
Jupyter Notebook
examples.ipynb
saadgroup/PyModPDE
5e9631aaece79cedc6f507dc52239f2a9de20d59
MIT
2020-02-20T16:03:40.000Z
2020-11-13T16:03:48.000Z
examples.ipynb
saadgroup/PyModPDE
5e9631aaece79cedc6f507dc52239f2a9de20d59
MIT
2020-03-14T21:13:03.000Z
2022-03-10T21:48:20.000Z
examples.ipynb
saadgroup/PyModPDE
5e9631aaece79cedc6f507dc52239f2a9de20d59
MIT
2020-02-20T03:58:22.000Z
2021-12-11T13:47:42.000Z
Qwen/Qwen-72B
1. YES
2. YES