HDG 1D

1. Imports
2. Mesh Class
3. Basis Class
- 3.1. Action 1: Construction of the basis functions (Class Basis)
  - 3.1.1. Multiplying two polynomials
  - 3.1.2. Evaluating a polynomial
- 3.2. Action 2: Construction of the reference matrices (Class Basis)
  - 3.2.1. Primitive of a polynomial
  - 3.2.2. Derivative of a polynomial
4. Elem Class
- 4.1. Action 3: Construction of the local matrices and rhs on each element (Class Elem)
5. Class Sparse
- 5.1. Action 4: Construction of the global matrix (Class Elem) and rhs
6. Main Program
7. Advanced testing
8. Correction
- 8.1. python version
- 8.2. notebook version

If you are reading this introduction outside of the notebook, note that python-only and notebook versions of the code are available. Note that within the notebook, these previous links are not available, you are already all set to start this hands-on session. You may want to have a look at this jupyter tutorial if you are not familiar with jupyter.

1. Imports

import numpy as np
import random
import numpy.polynomial.polynomial as polynom
import numpy.linalg as linalg
from scipy import sparse
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt

2. Mesh Class

class Mesh:
    def __init__(self, omega, number_media, length_media, rho_media, mu_media,
		 h_glob, order_glob):
	#omega : frequency
	#number_media : number of layers
	#length_media : a vector of size number_media containing the length of each layer
	#rho_media : a vector of size number_media containing the density of each layer
	#mu_media : a vector of size number_media containing the bulk modulus of each layer
	length = 0  #the total length
	rho = []  #a vector containing the density of each cell of the mesh
	mu = []  #a vector containing the bulk modulus of each cell of the mesh
	order = [
	]  #a vector containing the polynomial order of each cell of the mesh
	nbpoints = 1
	# the number of nodes of the mesh
	nbele = 0
	# the number of elements of the mesh
	coor = [
	    0
	]  #a vector of size nbpoints, containing the coordinates of each node of the mesh
	ele = [
	]  # a table of size nbele x 2 : the first columns gives the first node of the celle, the second gives the second nodes. It is useless in 1D, but we set it in order to mimic 2D and 3D
	L = 0
	for i in range(0, number_media):
	    #we loop over all the layers
	    L = L + length_media[
		i]  #we increase the total length by the size of the current layer
	    c = np.sqrt(
		mu_media[i] /
		rho_media[i])  #we compute the speed of the acoustic wave
	    h = h_glob  #the size of the cells in the layer #TBC (exo 8)
	    order_cell = order_glob  #the order of the cells in the layer
	    N = int(np.floor(length_media[i] /
			     h))  #we compute the number of cells of the layer
	    for j in range(0, N - 1):
		#we construct the cells of the layer i
		nbele = nbele + 1  #we increase the number  of cells
		nbpoints = nbpoints + 1  #we increase the number  of nodes
		length = length + h  #we increase the length by h
		coor.append(length)  #this is the coordinate of the new node
		rho.append(
		    rho_media[i])  #the density of the cell is rho_media[i]
		mu.append(
		    mu_media[i])  #the bulk modulus of the cell is mu_media[i]
		order.append(
		    order_cell)  #the bulk modulus of the cell is order_global
		ele.append(
		    nbpoints - 2
		)  #the two nodes of the elements are nbpoints-2 and nbpoints-1
		ele.append(nbpoints - 1)
	    #now, we define the last cell of the layer
	    nbele = nbele + 1
	    nbpoints = nbpoints + 1
	    length = L
	    coor.append(length)
	    rho.append(rho_media[i])
	    mu.append(mu_media[i])
	    order.append(order_cell)
	    ele.append(nbpoints - 2)
	    ele.append(nbpoints - 1)

	self.ele = np.reshape(np.array(ele), [nbele, 2])
	self.coor = np.array(coor)
	self.rho = np.array(rho)
	self.mu = np.array(mu)
	self.nbele = nbele
	self.nbpoints = nbpoints
	self.order = order

3. Basis Class

class Basis:
    def __init__(self, p):
	#p is the order of the basis
	#we build the p+1 elementary polynomials of degree 1
	poly_elem = np.zeros(
	    (p + 1, 2)
	)  #poly_elem is table (p+1)x2 : each line contains the coefficients a and b such that the polynom reads as a+bx.
	for i in range(0, p + 1):
	    #the root of polynonimal i is i/p, so that a=-i/p and b=1;
	    poly_elem[i, 0] = -(i) / p
	    poly_elem[i, 1] = 1
	    #we initialize the basis
	    phi = []
	    #phi is a vector of vector : phi[i] is the basis function i which is a polynomial of degree p, it represented by a vector of size p+1 containing the coefficients of the polynomial fom the lowest to the highest
	for i in range(0, p + 1):
	    #the basis function i is the product of all elementary polynomials except polynomial i. It is normalized so that phi[i](i)=1
	    #we start with the constant polynomial
	    phi.append([1])
	    for j in range(0, p + 1):
		if i != j:
		    phi[i] = 000000  #TBC (ex 1)
		    #we multiply phi[i] by all elementary polynomials except polynomial i
	    phi[i] = phi[i] / 000000  #TBC (ex 1)
	    #we normalize phi[i]
	self.phi = phi
	massmat = np.zeros((p + 1, p + 1))  #we initialize the mass matrix
	stiffmat = np.zeros(
	    (p + 1, p + 1))  #we initialize the stiffness matrix
	for i in range(0, p + 1):
	    for j in range(0, p + 1):
		massmat[i, j] = 000000  #TBC (ex 2)
		stiffmat[i, j] = 000000  #TBC (ex 2)
	self.massmat = massmat
	self.stiffmat = stiffmat

3.1. Action 1: Construction of the basis functions (Class Basis)

3.1.1. Multiplying two polynomials

The python command for multiplying two polynomials q and r is polynom.polymul(q,r). Compute the multiplication of polynomials phi[i] and poly_elem[j,:] and store it in phi[i].

3.1.2. Evaluating a polynomial

The python command for evaluating a polynomial q at point x is polynom.polyval(x,q). Normalize phi[i] by its value at point i/p.

3.2. Action 2: Construction of the reference matrices (Class Basis)

3.2.1. Primitive of a polynomial

The python command for computing the primitive of a polynomial P is polynom.polyint(P). Hence, \(\int_0^{1}P(x)dx\) can be computed thanks to sum(polynom.polyint(P)). Compute: \[massmat(i,j) = \int_0^1 \phi_i(x) \phi_j(x) dx\].

3.2.2. Derivative of a polynomial

The python command for computing the derivative of a polynomial is polynom.polyder(P). Compute: \[stiffmat(i,j)= \int_0^1 \phi_i(x) \frac{\partial\phi_j(x)}{\partial x} dx\]

4. Elem Class

#This class contains the local matrix and rhs
class Elem:
    def __init__(self, x0, x1, rho, mu, order, alpha, basis, omega, xs):
	#x0 : left coordinate of the cell
	#x1 :  right coordinate of the cell
	#rho : density of the cell
	#mu : bulk modulus of the cell
	#order : order of the cell
	#alpha : penalization coefficient
	#basis : the basis functions
	#omega : the frequency
	#xs : the position of the source, if any
	n = order + 1
	#the number of degree of freedom in the cell
	h = x1 - x0  #the length of the cell
	mat = np.zeros(
	    (2 * n, 2 * n), dtype=complex)  #initialization of the matrix
	rhs = np.zeros((2 * n), dtype=complex)  #initialization of the rhs

	#the four blocks of the matrix
	mat[:n, :n] = 000000  #TBC(exo3)
	mat[n:, :n] = 000000  #TBC(exo3)
	mat[:n, n:] = 000000  #TBC(exo3)
	mat[n:, n:] = 000000  #TBC(exo3)

	#we add the penalization coefficients
	mat[n, n] = 000000  #TBC(exo3)
	mat[2 * n - 1, 2 * n - 1] = 000000  #TBC(exo3)

	F = np.zeros((2 * n, 2))  #the right hand side
	F[0, 0] = 000000  #TBC(exo3)
	F[n - 1, 1] = 000000  #TBC(exo3)
	F[n, 0] = 000000  #TBC(exo3)
	F[2 * n - 1, 1] = 000000  #TBC(exo3)

	#computation of the source if it is inside the element
	if (xs >= x0 and xs <= x1):
	    #we compute the coordinates of the source in the reference element
	    xs_ref = (xs - x0) / h
	    for i in range(n):
		#we compute the values of basis functions at point xs_ref
		rhs[i] = polynom.polyval(xs_ref, basis.phi[i])
	    #the source is only on the pressure
	    rhs[n:] = 0
	    rhs = linalg.solve(mat, rhs)
	else:
	    rhs[:] = 0

	self.mat = linalg.solve(mat, F)
	self.rhs = rhs

4.1. Action 3: Construction of the local matrices and rhs on each element (Class Elem)

4.1.1. Compute the four blocks of the local matrix (without the penalization coefficients)

We have:

\(M(0:n-1, 0:n-1) = \frac{i\omega}{\mu} * Mass * h\)
\(M(n:2*n-1, 0:n-1) = -Stiff\)
\(M(0:n-1, n:2*n-1) = Stiff^T\)
\(M(n:2*n-1, n:2*n-1) = i\omega * \rho * Mass * h\)

Note that:

the \(i\) imaginary number is defined by 1j in python;
the mass matrix Mass is in basis.massmat;
the stiffness matrix Stiff is in basis.stiffmat and its transpose in basis.stiffmat.T.

4.1.2. Add the penalization coefficients

We have:

\(M(n,n) = M(n,n)-\alpha\)
\(M(2 * n-1, 2*n-1) = M(2*n-1,2*n-1)-\alpha\)

4.1.3. Construct the local rhs

5. Class Sparse

class Sparse:
    def __init__(self, alpha, elem, nbpoints, order):
	#alpha : penalization coefficient
	#elem : the local matrices associated to each elements
	#nbpoints : the number of vertices of the mesh
	#order : a vector containing the order of each element

	#the global matrices is stored in coordinate format :
	A = []
	#the vector containing the non-zero values
	I = []
	# the vector containing the lines of the non-zero values
	J = []
	# the vector containing the columns of the non-zero values
	#we have Mat(I[i], J[i])=A[i]
	B = []
	# the Right Hand Side

	#Left point of the mesh
	p = order[0]

	res = 000000  #TBC (exo4)
	A.append(res)
	I.append(0)
	J.append(0)

	res = 000000  #TBC (exo4)
	A.append(res)
	I.append(1)
	J.append(0)

	#right hand side
	res = 000000  #TBC (exo4)
	B.append(res)

	#now we run on all the vertices of the mesh (except the left and right ones)
	for i in range(1, nbpoints - 1):
	    #we get the order of elem i-1 and i
	    p1 = order[i - 1]
	    p2 = order[i]

	    res = 000000  #TBC (exo4)
	    A.append(res)
	    I.append(i - 1)
	    J.append(i)

	    res = 000000  #TBC (exo4)
	    A.append(res)
	    I.append(i)
	    J.append(i)

	    res = 000000  #TBC (exo4)
	    A.append(res)
	    I.append(i + 1)
	    J.append(i)

	    #right hand side
	    res = 000000  #TBC (exo4)
	    B.append(res)

	#right point of the mesh
	p = order[nbpoints - 2]
	res = 000000  #TBC (exo4)
	A.append(res)
	I.append(nbpoints - 2)
	J.append(nbpoints - 1)

	res = 000000  #TBC (exo4)
	A.append(res)
	I.append(nbpoints - 1)
	J.append(nbpoints - 1)

	res = 000000  #TBC (exo4)
	B.append(res)

	# Left and right boundary condition. Comment only if there is a point soucre
	B[0] = 0
	B[nbpoints - 1] = -np.sin(omega * mesh.coor[nbpoints - 1])
	self.matrix = sparse.coo_matrix(
	    (A, (I, J)), shape=(nbpoints, nbpoints))
	self.rhs = B

5.1. Action 4: Construction of the global matrix (Class Elem) and rhs

5.1.1. For internal vertices

The element left to vertex \(i\) is \(i-1\), the element right to vertex \(i\) is \(i\). We denote by \(p_1\) the order of cell \(i-1\) and by \(p_2\) the order of cell \(i\). Check that:

\(A(i-1,i)=-M_{i-1}(p_1, 0)-\alpha M_{i-1}(2p_1+1,0)\)
\(A(i,i)=-M_{i-1}(p_1,1)-\alpha M_{i-1}(2p_1+1,1)+M_i(0,0)-\alpha M_i(p_2+1,0)+2\alpha\)
\(A(i+1,i)=M_i(0,1)-\alpha M_i(p_2+1,1)\)
\(B(i)=F_{i-1}(0)-F_i(p_2)+\alpha(F_i(2p_2+1)+F_{i-1}(p1+1))\)

Note that \(M_{i}(j, k)\) corresponds to elem[i].mat[j,k] and that F_i(j) corresponds to elem[i].rhs[j] in the code.

5.1.2. For left vertex

The element right to vertex \(0\) is \(0\), we denote by \(p\) its order. Check that:

\(A(0,0)=M_0(0,0)-2 \alpha M_0(p+1,0)+2\alpha\)
\(A(1,0)=M_0(0,1)-2 \alpha M_0(p+1,1)\)
\(B(0)=-F_0(p)+2\alpha F_0(2p+1)\)

5.1.3. For right vertex

The element left to vertex \(nbpoints-1\) is \(nbpoints-2\), we denote by \(p\) its order. Check that:

\(A(nbpoints-2,nbpoints-1)=-M_{nbpoints-2}(p,0)-2 \alpha M_{nbpoints-2}(2p+1,0)\)
\(A(nbpoints-1,nbpoints-1)=-M_{nbpoints-2}(p,1)-2 \alpha M_{nbpoints-2}(2p+1,1)+2\alpha\)
\(B((nbpoints-1)=F_{nbpoints-2}(0)+2\alpha F_{nbpoints-2}(p+1)\)

6. Main Program

alpha = 0.5  #penalization parameter
omega = 10  #frequency
x_source = 15  #position of the source (if no source, put it outside of the domain)
nb_media = 1  #number of layer
points_plot = 10  #number of points on each cell for plotting the solution
h_glob = 0.1  #Global size of the cells
order_glob = 1  #Global order of the cells
length_media = [1]  #Length of the layers
rho_media = [1]  #density of the layers
mu_media = [1]  #bulk modulus of the layers
mesh = Mesh(omega, nb_media, length_media, rho_media, mu_media, h_glob,
	    order_glob)  #we create the mesh
basis = list(range(max(mesh.order) + 1))
for p in set(mesh.order):
    print("building basis functions of order", p)
    basis[p] = Basis(p)  #we define all the basis functions we need

elem = list(range(mesh.nbele + 1))
for i in range(mesh.nbele):
    #we run over all the cells of the mesh
    node = mesh.ele[i]
    x0 = mesh.coor[node[0]]
    x1 = mesh.coor[node[1]]
    rho = mesh.rho[i]
    mu = mesh.mu[i]
    order = mesh.order[i]
    elem[i] = Elem(x0, x1, rho, mu, order, alpha, basis[order], omega,
		   x_source)  #we construct the elementary basis functions

linsys = Sparse(alpha, elem, mesh.nbpoints,
		mesh.order)  # we create the global linear system
Lambda = 000000  #TBC (exo 5) #we solve the linear system
#plt.plot(np.imag(Lambda))
#plt.show()
#Computation of the solution

#Now, we reconstruct the solution on each element of the mesh
P = []
U = []
X = []
P_ex = []
U_ex = []
err_P = 0  #for computing the error between the approximate solution and the exact one
err_U = 0
L2_P = 0  #L2 norm of the exact solution
L2_U = 0
#we run over all the elements of the mesh
for i in range(mesh.nbele):
    #we find the two vertices of the element
    node = mesh.ele[i, :]
    #we get the two lambda associated to these vertices
    lambda_loc = Lambda[node]
    #we can now compute locally the solution on the element
    W_K = 000000  #TBC (exo 6)
    p = mesh.order[i]
    #the first degrees of freedom are associated to the velocity
    U_K = W_K[p + 1:]
    #the last degrees of freedom are associated to the pressure
    P_K = W_K[:p + 1]
    #we compute the value of P and U on the points of the submesh
    for j in range(points_plot + 1):
	x = j / points_plot
	xcoor = mesh.coor[node[0]] + x * (mesh.coor[node[1]] -
					  mesh.coor[node[0]])
	X.append(xcoor)
	P_exact = np.sin(omega * xcoor)
	U_exact = -1j * np.cos(omega * xcoor)
	P_ex.append(P_exact)
	U_ex.append(U_exact)
	valP = 0
	valU = 0
	#We get the order of the cell
	for k in range(p + 1):
	    val = 000000  #TBC (exo 6)
	    valP = valP + val * P_K[k]
	    valU = valU + val * U_K[k]
	P.append(valP)
	U.append(valU)
	err_P = err_P + abs(valP - P_exact) * abs(valP - P_exact)
	err_U = err_U + abs(valU - U_exact) * abs(valU - U_exact)
	L2_P = L2_P + abs(P_exact) * abs(P_exact)
	L2_U = L2_U + abs(U_exact) * abs(U_exact)
#computation of relative  L2 error
print("relative L2 error on P :", np.sqrt(err_P / L2_P))
print("relative L2 error on U :", np.sqrt(err_U / L2_U))

fig, axs = plt.subplots(2, 2)
ax1 = axs[0, 0]
ax1.plot(X, np.real(P), 'g')
ax1.plot(X, np.real(P_ex), 'b--')
ax1.set_title('P_real')
ax2 = axs[1, 0]
ax2.plot(X, np.imag(P), 'g')
ax2.plot(X, np.imag(P_ex), 'b--')
ax2.set_title('P_imag')
ax3 = axs[0, 1]
ax3.plot(X, np.real(U), 'g')
ax3.plot(X, np.real(U_ex), 'b--')
ax3.set_title('U_real')
ax4 = axs[1, 1]
ax4.plot(X, np.imag(U), 'g')
ax4.plot(X, np.imag(U_ex), 'b--')
ax4.set_title('U_imag')
plt.show()
fig, axs = plt.subplots(2, 2)
ax1 = axs[0, 0]
ax1.plot(X, np.real(P) - np.real(P_ex), 'g')
ax1.set_title('P_real')
ax2 = axs[1, 0]
ax2.plot(X, np.imag(P) - np.imag(P_ex), 'g')
ax2.set_title('P_imag')
ax3 = axs[0, 1]
ax3.plot(X, np.real(U) - np.real(U_ex), 'g')
ax3.set_title('U_real')
ax4 = axs[1, 1]
ax4.plot(X, np.imag(U) - np.imag(U_ex), 'g')
ax4.set_title('U_imag')
plt.show()

6.1. Action 5: Solution to the global system

The solution to a sparse linear system \(Ax=b\) is given by x=spsolve(A,b). Compute Λ, solution to linsys.matrix \Lambda= linsys.rhs.

6.2. Action 6: Reconstruction of the solution, element by element

6.2.1. The matrix vector product is given by np.dot(A,x)

Compute \(W_K=A_K^{-1} \lambda_{loc}+F_K\) for each element \(K\) (we remind that A_K^-1=elem[i].mat and F_K is elem[i].rhs)

6.2.2. In order to compute the solution at point x, compute φ_k(x) for each k

We remind that the value of a polynom p at point x is given by polynom.polyval(x,p)

6.3. Action 7: Validation and plotting

6.3.1. Check that the plotted solution is close to the exact one

6.3.2. Set the order to 1

Check that the relative error is divided by 4 when you divide the space step h by 2

6.3.3. Set the order to p=2…5

Check that the relative error is divided by 2^p+1 when you divide the space step h by 2

7. Advanced testing

7.1. Action 8: Role of h and p

7.1.1. Set h to 0.1

How does the error behaves when you increase the order, for p from 1 to 10 ?

7.1.2. Choosing the number of point per wavelength at a given order

The wavelength is \(\frac{2\pi c}{\omega}\), with \(c=\sqrt{\frac{\mu}{\rho}}\). Determine (roughly) the minimum number of point pers wavelength to obtain a relative error of 0.001, for p from 1 to 5.

7.1.3. Choosing the number of point per wavelength at a given h

Conversely, determine the minimal order to obtain a relative error of 0.001, for a number of point per wavelength from 15 to 1.

7.1.4. h-adaptivity

Modify the class Mesh in order to define h on each layer, depending on c, omega and order_glob How does the error behaves if you increase the frequency ?

7.1.5. p-adaptivity

Modify the class Mesh in order to define p on each layer, depending on c, omega and h_glob How does the error behaves if you increase the frequency ?

7.2. Action 9: Point source

7.2.1. Comment the two lines

B[0]=0;
B[nbpoints-1]=-np.sin(omega*mesh.coor[nbpoints-1]);

They are at the end of the Sparse class.

7.2.2. Set the source inside the domain (for instance x_source=0.3)

7.2.3. Remove the plots of the exact solution

7.3. Action 10: Multilayers

Set n_media=3, set length_media=[1,1,1 ], rho_media=[1,2,1] and mu_media=[1,1,2]. Run the code, using either the h- or the p-adaptivity.

8. Correction

The corrections below are to be fetched outside from the notebook. If you are in the notebook, you can directly load them from /home/hpcs-readonly/notebook (you can still retrieve the below links at the bottom of this page).

8.1. python version

step 1-8b, step 8d, step 8e, step 9, step 10

8.2. notebook version