<jupyter_start><jupyter_text>1. 2D Linear ConvectionWe consider the 1d linear Convection equation, under a constant velocity$$\partial_t u + \mathbf{a} \cdot \nabla u - \nu \nabla^2 u = 0$$<jupyter_code># needed imports
from numpy import zeros, ones, linspace, zeros_like
from matplotlib.pyplot import plot, contourf, show, colorbar
%matplotlib inline
# 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) 
contourf(x,y, u); colorbar() ; show()<jupyter_output><empty_output><jupyter_text>Time scheme$$\frac{u^{n+1}-u^n}{\Delta t} + \mathbf{a} \cdot \nabla u^{n+1} - \nu \nabla^2 u_{n+1} = 0 $$$$ \left(I + \Delta t \mathbf{a} \cdot \nabla - \nu \nabla^2 \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 + \nu ~\Delta t~ \langle \nabla v, \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 + \nu ~ \Delta t~\langle \nabla v, \nabla 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 + \nu ~\Delta t~ \langle \nabla b_i , \nabla b_j \rangle$$  Abstract Model using SymPDE<jupyter_code>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
# ... 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')
nu = Constant('nu')

a = Tuple(a1,a2)

# bilinear form
expr = v*u + dt* dot(a, grad(u))*v + nu*dt*dot(grad(u), grad(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))<jupyter_output><empty_output><jupyter_text>Discretization using Psydac<jupyter_code>from psydac.api.discretization import discretize
a1 = 1. ; a2 = 0. # wavespeed
nu = .03   # viscosity
T = 0.2  # T final time
dt = 0.001
niter = int(T / dt)

degree = [3,3]    # spline degree
ncells = [64,64]   # number of elements
# 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)
# assemble matrices and convert them to scipy
M = mh.assemble().tosparse()
A = ah.assemble(a1=a1, a2=a2, nu=nu, dt=dt).tosparse()

# assemble the rhs and convert it to numpy array
rhs = lh.assemble().toarray()
from scipy.sparse.linalg import gmres
# L2 projection of the initial condition
un, status = gmres(M, rhs, tol=1.e-8, maxiter=5000) 
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()
for i in range(0, niter):
    b = M.dot(un)
    un, status = gmres(A, b, tol=1.e-8, maxiter=5000)
nbasis = [W.nbasis for W in Vh.spaces]
plot_field_2d(Vh.knots, Vh.degree, un.reshape(nbasis)) ; colorbar() ; show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16653796653796654, 'fraction_numerical': 0.023938223938223938, 'mean_word_length': 4.4198047419804745, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '20159', 'license': 'MIT', 'n_tokens_mistral': 1584, 'n_tokens_neox': 1518, 'n_words': 438, 'path': 'lessons/Chapter2/06_convection_diffusion_2d.ipynb', 'repo_name': 'pyccel/IGA-Python', 'size': '3885'}