# *FEniCS tutorial:* Diffusion of a Gaussian hill in a two-dimensional disc

In this demo, we solve the two-dimensional diffusion equation with Dirichlet
bounadary conditions $u_D$ and source term $f = 0$. The initial condition,
$u_0$, we define as a Gaussian hill.

$$
\begin{align}
  u'   &= \nabla^2 u + f \quad\text{in the unit circle} \\
  u    &= u_D  \hphantom{u+f}\quad\text{on the boundary} \\
  u    &= u_0  \hphantom{u+f}\quad\;\text{at $t = 0$} \\
\end{align}
$$
with
$$
\begin{align}
  u_D  &= 0 \\
  u_0  &= \exp\bigl(-a(x^2 + y^2)\bigr) \\
  f    &= 0
\end{align}
$$


```python
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import time as tm

# import os
# import sys
# import re
# from IPython.display import Image
# from IPython.core.interactiveshell import InteractiveShell
# InteractiveShell.ast_node_interactivity = "all"
```


```python
from fenics import *
from mshr import *
```


```python
## set up the problem and define the main simulation function evolve()

# simulation parameters
# -- time step
T = 0.25          # total simulation time
n_steps = 15      # number of time steps
dt = T / n_steps  # size of time step
nip = 3           # number of intervals between plots
# -- mesh density
md = 24           # mesh density within disc

# create mesh and define function space
domain = Circle(Point(0, 0), 1)
mesh = generate_mesh(domain, md)
V = FunctionSpace(mesh, 'P', 2)

# define boundary condition
def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, Constant(0), boundary)

# define initial value
u_0 = Expression('exp(-a*pow(x[0] - x0, 2) - a*pow(x[1] - y0, 2))',
                 degree=2, x0=0., y0=0.30, a=5)
u_n = interpolate(u_0, V)

u_min = 0.
u_max = 1.

# define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
a, L = lhs(F), rhs(F)

# define time-evolution function
def evolve():

    # report initial state
    yield u_n

    # time-stepping
    u = Function(V)
    t = 0
    for n in range(1, n_steps + 1):

        # update current time
        t += dt

        # compute solution
        solve(a == L, u, bc)

        # report solution at nip-step intervals
        if n % nip == 0:
            yield u

        # update previous solution
        u_n.assign(u)
```


```python
n_rows = 3
n_cols = 5
fig_wd = 15
# default sizing here yields unit aspect ratio
plt.figure(figsize = (fig_wd, fig_wd * n_rows // n_cols))
plt.subplots_adjust(wspace=0.3, hspace=0.3)

idx = 0
for u in evolve():
    idx += 1
    plt.subplot(n_rows, n_cols, idx)
    plot(u, vmin=u_min, vmax=u_max)
```


```python
mesh
```


```python
u
```


```python
plot(u)
```


```python
matplotlib.tri.tricontour.TriContourSet?
```


```python
??plot
```


```python
??_plot_matplotlib
```


```python

```

d3e41ce386431fe1c8a91450aa94e53599bd40fb
ipynb
Jupyter Notebook
examples/notebooks/thermal_transport/xt03_heat_2D_circle.ipynb
radiasoft/rslaser
096ebf41f56a1f34cea91a597f804d427bea34bc
Apache-2.0
examples/notebooks/thermal_transport/xt03_heat_2D_circle.ipynb
radiasoft/rslaser
096ebf41f56a1f34cea91a597f804d427bea34bc
Apache-2.0
2020-06-29T17:04:49.000Z
2022-03-28T21:39:28.000Z
examples/notebooks/thermal_transport/xt03_heat_2D_circle.ipynb
radiasoft/rslaser
096ebf41f56a1f34cea91a597f804d427bea34bc
Apache-2.0
2021-02-25T13:28:30.000Z
2021-02-25T13:28:30.000Z
Qwen/Qwen-72B
1. YES
2. YES