# 1-D Diffusion equation

$$\frac{\partial u}{\partial t}= \nu \frac{\partial^2 u}{\partial x^2}$$



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


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

grid = linspace(0., 1., 401)
u = u0(grid) 
plot(grid, u) ; show()
```

### Time scheme

$$\frac{u^{n+1}-u^n}{\Delta t} - \nu \partial_{xx} u^{n+1} = 0 $$

$$ \left(I - \nu \Delta t \partial_{xx} \right) u^{n+1} = u^n $$


### Weak formulation

$$
\langle v, u^{n+1} \rangle + \nu \Delta t ~ \langle \partial_x v, \partial_x 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 + \nu \Delta t ~ \langle \partial_x b_i, \partial_x 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, Line, element_of, dx
from sympde.topology import dx1 # TODO: this is a bug right now
```


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

V = ScalarFunctionSpace('V', domain)

x = domain.coordinates

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

nu = Constant('nu')
dt = Constant('dt')

# bilinear form
# expr = v*u - c*dt*dx(v)*u # TODO BUG not working
expr = v*u + nu*dt*dx1(v)*dx1(u)
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-.5)**2/.05**2)*v
l = LinearForm(v, integral(domain, expr))
```

## Discretization using Psydac


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


```python
nu = 0.3   # viscosity
T = 0.02  # T final time
dt = 0.001
niter = int(T / dt)

degree = [3]    # spline degree
ncells = [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(nu=nu, dt=dt).tosparse()

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


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


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


```python
from simplines import plot_field_1d
plot_field_1d(Vh.knots[0], Vh.degree[0], un, nx=401)
```


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


```python
plot_field_1d(Vh.knots[0], Vh.degree[0], un, nx=401)
```


```python

```

ac9bf615b1c504710a40b3d9459c2e202ff86e43
ipynb
Jupyter Notebook
lessons/Chapter2/02_diffusion_1d.ipynb
ratnania/IGA-Python
a9d7aa9bd14d4b3f1b12cdfbc2f9bf3c0a68fff4
MIT
2018-04-27T15:40:17.000Z
2020-08-13T08:45:35.000Z
lessons/Chapter2/02_diffusion_1d.ipynb
GabrielJie/IGA-Python
a9d7aa9bd14d4b3f1b12cdfbc2f9bf3c0a68fff4
MIT
2021-06-08T22:59:19.000Z
2022-01-17T20:36:56.000Z
lessons/Chapter2/02_diffusion_1d.ipynb
GabrielJie/IGA-Python
a9d7aa9bd14d4b3f1b12cdfbc2f9bf3c0a68fff4
MIT
2018-10-06T01:30:20.000Z
2021-12-31T02:42:05.000Z
Qwen/Qwen-72B
1. YES
2. YES