# 1. 2D Linear Convection

We consider the 1d linear Convection equation, under a constant velocity

$$
\partial_t u + \mathbf{a} \cdot \nabla u = 0
$$


```python
# needed imports
from numpy import zeros, ones, linspace, zeros_like
from matplotlib.pyplot import plot, contourf, show, colorbar
%matplotlib inline
```


```python
# Initial condition
import numpy as np
u0 = lambda x,y: np.exp(-(x-.3)**2/.05**2)*np.exp(-(y-.3)**2/.05**2)

ts = linspace(0., 1., 401)
x,y = np.meshgrid(ts,ts)
u = u0(x,y) 
```


```python
contourf(x,y, u); colorbar() ; show()
```

### Time scheme

$$\frac{u^{n+1}-u^n}{\Delta t} + \mathbf{a} \cdot \nabla u^{n+1} = 0 $$

$$ \left(I + \Delta t \mathbf{a} \cdot \nabla \right) u^{n+1} = u^n $$

### Weak formulation

$$
\langle v, u^{n+1} \rangle + \Delta t ~ \langle v, \mathbf{a} \cdot \nabla u^{n+1} \rangle = \langle v, u^n \rangle
$$
if we assume $\mathbf{a} = \left( a_1, a_2 \right)^T$ is a constant, then our weak formulation writes
$$
\langle v, u^{n+1} \rangle - \Delta t ~ \langle \mathbf{a} \cdot \nabla v ,  u^{n+1} \rangle = \langle v, u^n \rangle
$$

expending $u^n$ over the fem basis, we get the linear system

$$A U^{n+1} = M U^n$$

where

$$
M_{ij} = \langle b_i, b_j \rangle
$$

$$
A_{ij} = \langle b_i, b_j \rangle - \Delta t ~ \langle \mathbf{a} \cdot \nabla b_i, b_j \rangle
$$


## Abstract Model using SymPDE


```python
from sympde.core import Constant
from sympde.expr import BilinearForm, LinearForm, integral
from sympde.topology import ScalarFunctionSpace, Square, element_of
from sympde.calculus import grad, dot
from sympy import Tuple
```


```python
# ... abstract model
domain = Square()

V = ScalarFunctionSpace('V', domain)

x,y = domain.coordinates

u,v = [element_of(V, name=i) for i in ['u', 'v']]

a1 = Constant('a1')
a2 = Constant('a2')
dt = Constant('dt')

a = Tuple(a1,a2)

# bilinear form
expr = v*u + dt* dot(a, grad(u))*v
a = BilinearForm((u,v), integral(domain , expr))

# bilinear form for the mass matrix
expr = u*v
m  = BilinearForm((u,v), integral(domain , expr))

# linear form for initial condition
from sympy import exp
expr = exp(-(x-.3)**2/.05**2)*exp(-(y-.3)**2/.05**2)*v
l = LinearForm(v, integral(domain, expr))
```

## Discretization using Psydac


```python
from psydac.api.discretization import discretize
```


```python
a1 = 1. ; a2 = 0. # wavespeed
T = 0.25  # T final time
dt = 0.001
niter = int(T / dt)

degree = [3,3]    # spline degree
ncells = [64,64]   # number of elements
```


```python
# Create computational domain from topological domain
domain_h = discretize(domain, ncells=ncells, comm=None)

# Discrete spaces
Vh = discretize(V, domain_h, degree=degree)

# Discretize the bilinear forms
ah = discretize(a, domain_h, [Vh, Vh])
mh = discretize(m, domain_h, [Vh, Vh])

# Discretize the linear form for the initial condition
lh = discretize(l, domain_h, Vh)
```


```python
# assemble matrices and convert them to scipy
M = mh.assemble().tosparse()
A = ah.assemble(a1=a1, a2=a2, dt=dt).tosparse()

# assemble the rhs and convert it to numpy array
rhs = lh.assemble().toarray()
```


```python
from scipy.sparse.linalg import gmres
```


```python
# L2 projection of the initial condition
un, status = gmres(M, rhs, tol=1.e-8, maxiter=5000) 
```


```python
from utilities.plot import plot_field_2d

nbasis = [W.nbasis for W in Vh.spaces]
plot_field_2d(Vh.knots, Vh.degree, un.reshape(nbasis)) ; colorbar() ; show()
```


```python
for i in range(0, niter):
    b = M.dot(un)
    un, status = gmres(A, b, tol=1.e-8, maxiter=5000)
```


```python
nbasis = [W.nbasis for W in Vh.spaces]
plot_field_2d(Vh.knots, Vh.degree, un.reshape(nbasis)) ; colorbar() ; show()
```


```python

```


```python

```

def980788d984b837facc48701a8109d69405c4b
ipynb
Jupyter Notebook
lessons/Chapter2/04_convection_2d.ipynb
pyccel/IGA-Python
e3604ba3d76a20e3d30ed3c7c952dcd2dc8147bb
MIT
2022-01-21T08:51:30.000Z
2022-03-17T12:14:02.000Z
lessons/Chapter2/04_convection_2d.ipynb
pyccel/IGA-Python
e3604ba3d76a20e3d30ed3c7c952dcd2dc8147bb
MIT
lessons/Chapter2/04_convection_2d.ipynb
pyccel/IGA-Python
e3604ba3d76a20e3d30ed3c7c952dcd2dc8147bb
MIT
2022-03-01T06:41:54.000Z
2022-03-01T06:41:54.000Z
Qwen/Qwen-72B
1. YES
2. YES

