# Examples on Using MOIRA
---
__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
# %config InlineBackend.figure_format = 'svg'
import matplotlib.pyplot as plt
# plt.rcParams['figure.figsize'] = [3.5, 4]
plt.rc("font", family='serif')
# plt.rc('text', usetex=True)
import numpy as np
from sympy import symbols, pretty_print
from IPython.display import display, Math, clear_output
from matplotlib import cm, animation
```

##  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.MOIRA import DifferentialEquation, i, 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(dependentVar='u',independentVars=['x']) # construct a one D differential equation with 'x' and 'u' as independent and dependent variables, respectively

# method I of constructing the rhs:
advectionTerm1 = DE1.expr(points=[1, 0], direction='x', order=1, time=n+1) # 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

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

# compute the modified equation up to two terms
DE1.modified_equation(nterms=3)

# display the modified equation in latex form
display(Math(DE1.latex()))

```


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



```python
from src.MOIRA 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(dependentVar='f',independentVars=['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.modified_equation(nterms=2)  # compute the modified equation up to 2 terms

display(Math(DE2.latex())) # display the latex form of the modified equation
```


$$\frac{\partial}{\partial t} f{\left (y,t \right )} =  - b \frac{\partial}{\partial y} f{\left (y,t \right )} + \frac{b}{2} \left(\Delta y - \Delta{t} b\right) \frac{\partial^{2}}{\partial y^{2}}  f{\left (y,t \right )}$$


## 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
\begin{equation}
\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).
\end{equation}


```python
from src.MOIRA 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(dependentVar='u',independentVars=['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 + g*diffusionTerm3 ) # set the rhs of the Differential equation DE

DE3.modified_equation(nterms=3)  # compute the modified equation up to 2 terms

display(Math(DE3.latex())) # display the latex form of the modified equation
```


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


## 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.MOIRA import DifferentialEquation, i, 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(dependentVar='u',independentVars=['x','y']) # construct a one D differential equation with 'x' and 'u' as independent and dependent variables, respectively

# method I of constructing the rhs:
advectionTerm1 = - b * DE1.expr(points=[-1, 0], direction='x', order=1, time=n)  - c * DE1.expr(points=[-1, 0], direction='y', order=1, time=n)

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

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

# display the modified equation in latex form
display(Math(DE1.latex()))
```


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


### 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(-1,1,20)
plt.figure(figsize=(10,8))
plt.contourf(xx,yy,sol[-1][1:-1,1:-1]+sol[0][1:-1,1:-1],cmap=cm.jet,levels=levs,vmax=1.0,vmin=-1.0)
plt.text(0.1, 0.4, 'time=0.0s', fontsize=12)
plt.text(0.6, 0.92,'time={}s'.format(tend), fontsize=12)
plt.xlabel('x', fontsize=14)
plt.ylabel('y', fontsize=14)
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…


70af0e0d7d5422bc351241ac468a294338e12009
ipynb
Jupyter Notebook
examples.ipynb
saadtony/modified_equation_code_test
551f40e88a6c0f011e8d34cd7d821139add42aaf
MIT
examples.ipynb
saadtony/modified_equation_code_test
551f40e88a6c0f011e8d34cd7d821139add42aaf
MIT
examples.ipynb
saadtony/modified_equation_code_test
551f40e88a6c0f011e8d34cd7d821139add42aaf
MIT
2021-12-28T19:05:55.000Z
2021-12-28T19:05:55.000Z
Qwen/Qwen-72B
1. YES
2. YES