{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "If you are reading this introduction outside of the notebook, note that\n", "[python](./hdg.py)-only and [notebook](./notebook/hdg-step0.ipynb) versions of the code are available. Note that within\n", "the notebook, these previous links are not available, you are already all set to\n", "start this hands-on session. You may want to have a look at [this jupyter\n", "tutorial](https://www.dataquest.io/blog/jupyter-notebook-tutorial/) if you are not familiar with jupyter.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Imports\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import random\n", "import numpy.polynomial.polynomial as polynom\n", "import numpy.linalg as linalg\n", "from scipy import sparse\n", "from scipy.sparse.linalg import spsolve\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mesh Class\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "class Mesh:\n", " def __init__(self, omega, number_media, length_media, rho_media, mu_media,\n", " h_glob, order_glob):\n", " #omega : frequency\n", " #number_media : number of layers\n", " #length_media : a vector of size number_media containing the length of each layer\n", " #rho_media : a vector of size number_media containing the density of each layer\n", " #mu_media : a vector of size number_media containing the bulk modulus of each layer\n", " length = 0 #the total length\n", " rho = [] #a vector containing the density of each cell of the mesh\n", " mu = [] #a vector containing the bulk modulus of each cell of the mesh\n", " order = [\n", " ] #a vector containing the polynomial order of each cell of the mesh\n", " nbpoints = 1\n", " # the number of nodes of the mesh\n", " nbele = 0\n", " # the number of elements of the mesh\n", " coor = [\n", " 0\n", " ] #a vector of size nbpoints, containing the coordinates of each node of the mesh\n", " ele = [\n", " ] # 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\n", " L = 0\n", " for i in range(0, number_media):\n", " #we loop over all the layers\n", " L = L + length_media[\n", " i] #we increase the total length by the size of the current layer\n", " c = np.sqrt(\n", " mu_media[i] /\n", " rho_media[i]) #we compute the speed of the acoustic wave\n", " h = h_glob #the size of the cells in the layer #TBC (exo 8)\n", " order_cell = order_glob #the order of the cells in the layer\n", " N = int(np.floor(length_media[i] /\n", " h)) #we compute the number of cells of the layer\n", " for j in range(0, N - 1):\n", " #we construct the cells of the layer i\n", " nbele = nbele + 1 #we increase the number of cells\n", " nbpoints = nbpoints + 1 #we increase the number of nodes\n", " length = length + h #we increase the length by h\n", " coor.append(length) #this is the coordinate of the new node\n", " rho.append(\n", " rho_media[i]) #the density of the cell is rho_media[i]\n", " mu.append(\n", " mu_media[i]) #the bulk modulus of the cell is mu_media[i]\n", " order.append(\n", " order_cell) #the bulk modulus of the cell is order_global\n", " ele.append(\n", " nbpoints - 2\n", " ) #the two nodes of the elements are nbpoints-2 and nbpoints-1\n", " ele.append(nbpoints - 1)\n", " #now, we define the last cell of the layer\n", " nbele = nbele + 1\n", " nbpoints = nbpoints + 1\n", " length = L\n", " coor.append(length)\n", " rho.append(rho_media[i])\n", " mu.append(mu_media[i])\n", " order.append(order_cell)\n", " ele.append(nbpoints - 2)\n", " ele.append(nbpoints - 1)\n", "\n", " self.ele = np.reshape(np.array(ele), [nbele, 2])\n", " self.coor = np.array(coor)\n", " self.rho = np.array(rho)\n", " self.mu = np.array(mu)\n", " self.nbele = nbele\n", " self.nbpoints = nbpoints\n", " self.order = order" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Basis Class\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "class Basis:\n", " def __init__(self, p):\n", " #p is the order of the basis\n", " #we build the p+1 elementary polynomials of degree 1\n", " poly_elem = np.zeros(\n", " (p + 1, 2)\n", " ) #poly_elem is table (p+1)x2 : each line contains the coefficients a and b such that the polynom reads as a+bx.\n", " for i in range(0, p + 1):\n", " #the root of polynonimal i is i/p, so that a=-i/p and b=1;\n", " poly_elem[i, 0] = -(i) / p\n", " poly_elem[i, 1] = 1\n", " #we initialize the basis\n", " phi = []\n", " #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\n", " for i in range(0, p + 1):\n", " #the basis function i is the product of all elementary polynomials except polynomial i. It is normalized so that phi[i](i)=1\n", " #we start with the constant polynomial\n", " phi.append()\n", " for j in range(0, p + 1):\n", " if i != j:\n", " phi[i] = 000000 #TBC (ex 1)\n", " #we multiply phi[i] by all elementary polynomials except polynomial i\n", " phi[i] = phi[i] / 000000 #TBC (ex 1)\n", " #we normalize phi[i]\n", " self.phi = phi\n", " massmat = np.zeros((p + 1, p + 1)) #we initialize the mass matrix\n", " stiffmat = np.zeros(\n", " (p + 1, p + 1)) #we initialize the stiffness matrix\n", " for i in range(0, p + 1):\n", " for j in range(0, p + 1):\n", " massmat[i, j] = 000000 #TBC (ex 2)\n", " stiffmat[i, j] = 000000 #TBC (ex 2)\n", " self.massmat = massmat\n", " self.stiffmat = stiffmat" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Action 1: Construction of the basis functions (Class Basis)\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Multiplying two polynomials\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The python command for multiplying two polynomials q and r is\n", " polynom.polymul(q,r). Compute the multiplication of polynomials phi[i] and\n", " poly_elem[j,:] and store it in phi[i].\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Evaluating a polynomial\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The python command for evaluating a polynomial q at point x is\n", " polynom.polyval(x,q). Normalize phi[i] by its value at point i/p.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Action 2: Construction of the reference matrices (Class Basis)\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Primitive of a polynomial\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The python command for computing the primitive of a polynomial P is\n", " polynom.polyint(P). Hence, $\\int_0^{1}P(x)dx$ can be computed thanks to\n", " sum(polynom.polyint(P)). Compute: $$massmat(i,j) = \\int_0^1 \\phi_i(x) \\phi_j(x) dx$$.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Derivative of a polynomial\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The python command for computing the derivative of a polynomial is\n", " polynom.polyder(P). Compute: $$stiffmat(i,j)= \\int_0^1 \\phi_i(x)\n", " \\frac{\\partial\\phi_j(x)}{\\partial x} dx$$\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Elem Class\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "#This class contains the local matrix and rhs\n", "class Elem:\n", " def __init__(self, x0, x1, rho, mu, order, alpha, basis, omega, xs):\n", " #x0 : left coordinate of the cell\n", " #x1 : right coordinate of the cell\n", " #rho : density of the cell\n", " #mu : bulk modulus of the cell\n", " #order : order of the cell\n", " #alpha : penalization coefficient\n", " #basis : the basis functions\n", " #omega : the frequency\n", " #xs : the position of the source, if any\n", " n = order + 1\n", " #the number of degree of freedom in the cell\n", " h = x1 - x0 #the length of the cell\n", " mat = np.zeros(\n", " (2 * n, 2 * n), dtype=complex) #initialization of the matrix\n", " rhs = np.zeros((2 * n), dtype=complex) #initialization of the rhs\n", "\n", " #the four blocks of the matrix\n", " mat[:n, :n] = 000000 #TBC(exo3)\n", " mat[n:, :n] = 000000 #TBC(exo3)\n", " mat[:n, n:] = 000000 #TBC(exo3)\n", " mat[n:, n:] = 000000 #TBC(exo3)\n", "\n", " #we add the penalization coefficients\n", " mat[n, n] = 000000 #TBC(exo3)\n", " mat[2 * n - 1, 2 * n - 1] = 000000 #TBC(exo3)\n", "\n", " F = np.zeros((2 * n, 2)) #the right hand side\n", " F[0, 0] = 000000 #TBC(exo3)\n", " F[n - 1, 1] = 000000 #TBC(exo3)\n", " F[n, 0] = 000000 #TBC(exo3)\n", " F[2 * n - 1, 1] = 000000 #TBC(exo3)\n", "\n", " #computation of the source if it is inside the element\n", " if (xs >= x0 and xs <= x1):\n", " #we compute the coordinates of the source in the reference element\n", " xs_ref = (xs - x0) / h\n", " for i in range(n):\n", " #we compute the values of basis functions at point xs_ref\n", " rhs[i] = polynom.polyval(xs_ref, basis.phi[i])\n", " #the source is only on the pressure\n", " rhs[n:] = 0\n", " rhs = linalg.solve(mat, rhs)\n", " else:\n", " rhs[:] = 0\n", "\n", " self.mat = linalg.solve(mat, F)\n", " self.rhs = rhs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Action 3: Construction of the local matrices and rhs on each element (Class Elem)\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Compute the four blocks of the local matrix (without the penalization coefficients)\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have:\n", "\n", "1. $M(0:n-1, 0:n-1) = \\frac{i\\omega}{\\mu} * Mass * h$\n", "2. $M(n:2*n-1, 0:n-1) = -Stiff$\n", "3. $M(0:n-1, n:2*n-1) = Stiff^T$\n", "4. $M(n:2*n-1, n:2*n-1) = i\\omega * \\rho * Mass * h$\n", "\n", "Note that:\n", "\n", "- the $i$ imaginary number is defined by 1j in python;\n", "- the mass matrix Mass is in basis.massmat;\n", "- the stiffness matrix Stiff is in basis.stiffmat and its transpose in basis.stiffmat.T.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Add the penalization coefficients\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have:\n", "\n", "1. $M(n,n) = M(n,n)-\\alpha$\n", "2. $M(2 * n-1, 2*n-1) = M(2*n-1,2*n-1)-\\alpha$\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Construct the local rhs\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Class Sparse\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "class Sparse:\n", " def __init__(self, alpha, elem, nbpoints, order):\n", " #alpha : penalization coefficient\n", " #elem : the local matrices associated to each elements\n", " #nbpoints : the number of vertices of the mesh\n", " #order : a vector containing the order of each element\n", "\n", " #the global matrices is stored in coordinate format :\n", " A = []\n", " #the vector containing the non-zero values\n", " I = []\n", " # the vector containing the lines of the non-zero values\n", " J = []\n", " # the vector containing the columns of the non-zero values\n", " #we have Mat(I[i], J[i])=A[i]\n", " B = []\n", " # the Right Hand Side\n", "\n", " #Left point of the mesh\n", " p = order\n", "\n", " res = 000000 #TBC (exo4)\n", " A.append(res)\n", " I.append(0)\n", " J.append(0)\n", "\n", " res = 000000 #TBC (exo4)\n", " A.append(res)\n", " I.append(1)\n", " J.append(0)\n", "\n", " #right hand side\n", " res = 000000 #TBC (exo4)\n", " B.append(res)\n", "\n", " #now we run on all the vertices of the mesh (except the left and right ones)\n", " for i in range(1, nbpoints - 1):\n", " #we get the order of elem i-1 and i\n", " p1 = order[i - 1]\n", " p2 = order[i]\n", "\n", " res = 000000 #TBC (exo4)\n", " A.append(res)\n", " I.append(i - 1)\n", " J.append(i)\n", "\n", " res = 000000 #TBC (exo4)\n", " A.append(res)\n", " I.append(i)\n", " J.append(i)\n", "\n", " res = 000000 #TBC (exo4)\n", " A.append(res)\n", " I.append(i + 1)\n", " J.append(i)\n", "\n", " #right hand side\n", " res = 000000 #TBC (exo4)\n", " B.append(res)\n", "\n", " #right point of the mesh\n", " p = order[nbpoints - 2]\n", " res = 000000 #TBC (exo4)\n", " A.append(res)\n", " I.append(nbpoints - 2)\n", " J.append(nbpoints - 1)\n", "\n", " res = 000000 #TBC (exo4)\n", " A.append(res)\n", " I.append(nbpoints - 1)\n", " J.append(nbpoints - 1)\n", "\n", " res = 000000 #TBC (exo4)\n", " B.append(res)\n", "\n", " # Left and right boundary condition. Comment only if there is a point soucre\n", " B = 0\n", " B[nbpoints - 1] = -np.sin(omega * mesh.coor[nbpoints - 1])\n", " self.matrix = sparse.coo_matrix(\n", " (A, (I, J)), shape=(nbpoints, nbpoints))\n", " self.rhs = B" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Action 4: Construction of the global matrix (Class Elem) and rhs\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### For internal vertices\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The element left to vertex $i$ is $i-1$, the element right to vertex $i$ is $i$.\n", " We denote by $p_1$ the order of cell $i-1$ and by $p_2$ the order of cell $i$.\n", " Check that:\n", "\n", "1. $A(i-1,i)=-M_{i-1}(p_1, 0)-\\alpha M_{i-1}(2p_1+1,0)$\n", "2. $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$\n", "3. $A(i+1,i)=M_i(0,1)-\\alpha M_i(p_2+1,1)$\n", "4. $B(i)=F_{i-1}(0)-F_i(p_2)+\\alpha(F_i(2p_2+1)+F_{i-1}(p1+1))$\n", "\n", "Note that $M_{i}(j, k)$ corresponds to elem[i].mat[j,k] and that\n", "Fi(j) corresponds to elem[i].rhs[j] in the code.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### For left vertex\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The element right to vertex $0$ is $0$, we denote by $p$ its order.\n", " Check that:\n", "\n", "1. $A(0,0)=M_0(0,0)-2 \\alpha M_0(p+1,0)+2\\alpha$\n", "2. $A(1,0)=M_0(0,1)-2 \\alpha M_0(p+1,1)$\n", "3. $B(0)=-F_0(p)+2\\alpha F_0(2p+1)$\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### For right vertex\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The element left to vertex $nbpoints-1$ is $nbpoints-2$, we denote by $p$ its order.\n", " Check that:\n", "\n", "1. $A(nbpoints-2,nbpoints-1)=-M_{nbpoints-2}(p,0)-2 \\alpha M_{nbpoints-2}(2p+1,0)$\n", "2. $A(nbpoints-1,nbpoints-1)=-M_{nbpoints-2}(p,1)-2 \\alpha M_{nbpoints-2}(2p+1,1)+2\\alpha$\n", "3. $B((nbpoints-1)=F_{nbpoints-2}(0)+2\\alpha F_{nbpoints-2}(p+1)$\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Main Program\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "alpha = 0.5 #penalization parameter\n", "omega = 10 #frequency\n", "x_source = 15 #position of the source (if no source, put it outside of the domain)\n", "nb_media = 1 #number of layer\n", "points_plot = 10 #number of points on each cell for plotting the solution\n", "h_glob = 0.1 #Global size of the cells\n", "order_glob = 1 #Global order of the cells\n", "length_media =  #Length of the layers\n", "rho_media =  #density of the layers\n", "mu_media =  #bulk modulus of the layers\n", "mesh = Mesh(omega, nb_media, length_media, rho_media, mu_media, h_glob,\n", " order_glob) #we create the mesh\n", "basis = list(range(max(mesh.order) + 1))\n", "for p in set(mesh.order):\n", " print(\"building basis functions of order\", p)\n", " basis[p] = Basis(p) #we define all the basis functions we need\n", "\n", "elem = list(range(mesh.nbele + 1))\n", "for i in range(mesh.nbele):\n", " #we run over all the cells of the mesh\n", " node = mesh.ele[i]\n", " x0 = mesh.coor[node]\n", " x1 = mesh.coor[node]\n", " rho = mesh.rho[i]\n", " mu = mesh.mu[i]\n", " order = mesh.order[i]\n", " elem[i] = Elem(x0, x1, rho, mu, order, alpha, basis[order], omega,\n", " x_source) #we construct the elementary basis functions\n", "\n", "linsys = Sparse(alpha, elem, mesh.nbpoints,\n", " mesh.order) # we create the global linear system\n", "Lambda = 000000 #TBC (exo 5) #we solve the linear system\n", "#plt.plot(np.imag(Lambda))\n", "#plt.show()\n", "#Computation of the solution\n", "\n", "#Now, we reconstruct the solution on each element of the mesh\n", "P = []\n", "U = []\n", "X = []\n", "P_ex = []\n", "U_ex = []\n", "err_P = 0 #for computing the error between the approximate solution and the exact one\n", "err_U = 0\n", "L2_P = 0 #L2 norm of the exact solution\n", "L2_U = 0\n", "#we run over all the elements of the mesh\n", "for i in range(mesh.nbele):\n", " #we find the two vertices of the element\n", " node = mesh.ele[i, :]\n", " #we get the two lambda associated to these vertices\n", " lambda_loc = Lambda[node]\n", " #we can now compute locally the solution on the element\n", " W_K = 000000 #TBC (exo 6)\n", " p = mesh.order[i]\n", " #the first degrees of freedom are associated to the velocity\n", " U_K = W_K[p + 1:]\n", " #the last degrees of freedom are associated to the pressure\n", " P_K = W_K[:p + 1]\n", " #we compute the value of P and U on the points of the submesh\n", " for j in range(points_plot + 1):\n", " x = j / points_plot\n", " xcoor = mesh.coor[node] + x * (mesh.coor[node] -\n", " mesh.coor[node])\n", " X.append(xcoor)\n", " P_exact = np.sin(omega * xcoor)\n", " U_exact = -1j * np.cos(omega * xcoor)\n", " P_ex.append(P_exact)\n", " U_ex.append(U_exact)\n", " valP = 0\n", " valU = 0\n", " #We get the order of the cell\n", " for k in range(p + 1):\n", " val = 000000 #TBC (exo 6)\n", " valP = valP + val * P_K[k]\n", " valU = valU + val * U_K[k]\n", " P.append(valP)\n", " U.append(valU)\n", " err_P = err_P + abs(valP - P_exact) * abs(valP - P_exact)\n", " err_U = err_U + abs(valU - U_exact) * abs(valU - U_exact)\n", " L2_P = L2_P + abs(P_exact) * abs(P_exact)\n", " L2_U = L2_U + abs(U_exact) * abs(U_exact)\n", "#computation of relative L2 error\n", "print(\"relative L2 error on P :\", np.sqrt(err_P / L2_P))\n", "print(\"relative L2 error on U :\", np.sqrt(err_U / L2_U))\n", "\n", "fig, axs = plt.subplots(2, 2)\n", "ax1 = axs[0, 0]\n", "ax1.plot(X, np.real(P), 'g')\n", "ax1.plot(X, np.real(P_ex), 'b--')\n", "ax1.set_title('P_real')\n", "ax2 = axs[1, 0]\n", "ax2.plot(X, np.imag(P), 'g')\n", "ax2.plot(X, np.imag(P_ex), 'b--')\n", "ax2.set_title('P_imag')\n", "ax3 = axs[0, 1]\n", "ax3.plot(X, np.real(U), 'g')\n", "ax3.plot(X, np.real(U_ex), 'b--')\n", "ax3.set_title('U_real')\n", "ax4 = axs[1, 1]\n", "ax4.plot(X, np.imag(U), 'g')\n", "ax4.plot(X, np.imag(U_ex), 'b--')\n", "ax4.set_title('U_imag')\n", "plt.show()\n", "fig, axs = plt.subplots(2, 2)\n", "ax1 = axs[0, 0]\n", "ax1.plot(X, np.real(P) - np.real(P_ex), 'g')\n", "ax1.set_title('P_real')\n", "ax2 = axs[1, 0]\n", "ax2.plot(X, np.imag(P) - np.imag(P_ex), 'g')\n", "ax2.set_title('P_imag')\n", "ax3 = axs[0, 1]\n", "ax3.plot(X, np.real(U) - np.real(U_ex), 'g')\n", "ax3.set_title('U_real')\n", "ax4 = axs[1, 1]\n", "ax4.plot(X, np.imag(U) - np.imag(U_ex), 'g')\n", "ax4.set_title('U_imag')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Action 5: Solution to the global system\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The solution to a sparse linear system $Ax=b$ is given by x=spsolve(A,b). \n", " Compute Λ, solution to linsys.matrix \\Lambda= linsys.rhs.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Action 6: Reconstruction of the solution, element by element\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### The matrix vector product is given by np.dot(A,x)\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute $W_K=A_K^{-1} \\lambda_{loc}+F_K$ for each element $K$\n", " (we remind that AK-1=elem[i].mat and FK is elem[i].rhs)\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### In order to compute the solution at point x, compute \\phi_k(x) for each k\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We remind that the value of a polynom p at point x is given by polynom.polyval(x,p)\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Action 7: Validation and plotting\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Check that the plotted solution is close to the exact one\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Set the order to 1\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check that the relative error is divided by 4 when you divide the space step h by 2\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Set the order to p=2...5\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check that the relative error is divided by 2p+1 when you divide the space step h by 2\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Advanced testing\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Action 8: Role of h and p\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Set h to 0.1\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How does the error behaves when you increase the order, for p from 1 to 10 ?\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Choosing the number of point per wavelength at a given order\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The wavelength is $\\frac{2\\pi c}{\\omega}$, with $c=\\sqrt{\\frac{\\mu}{\\rho}}$. Determine\n", " (roughly) the minimum number of point pers wavelength to obtain a relative error\n", " of 0.001, for p from 1 to 5.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Choosing the number of point per wavelength at a given h\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Conversely, determine the minimal order to obtain a relative error of 0.001, for a number of point per wavelength from 15 to 1.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### h-adaptivity\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Modify the class Mesh in order to define h on each layer, depending on c, omega and orderglob\n", " How does the error behaves if you increase the frequency ?\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### p-adaptivity\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Modify the class Mesh in order to define p on each layer, depending on c, omega and hglob\n", " How does the error behaves if you increase the frequency ?\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Action 9: Point source\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Comment the two lines\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- B=0;\n", "- B[nbpoints-1]=-np.sin(omega*mesh.coor[nbpoints-1]);\n", "\n", "They are at the end of the Sparse class.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Set the source inside the domain (for instance x_source=0.3)\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Remove the plots of the exact solution\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Action 10: Multilayers\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set n_media=3, set length_media=[1,1,1 ], rho_media=[1,2,1] and\n", " mu_media=[1,1,2]. Run the code, using either the h- or the p-adaptivity.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Correction\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The corrections below are to be fetched outside from the notebook. If you are in\n", "the notebook, you can directly load them from /home/hpcs-readonly/notebook\n", "(you can still retrieve the below links at the bottom of [this page](https://agullo-teach.gitlabpages.inria.fr/school/school2019/tps/hdg-1d/README.html)).\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### python version\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[step 1-8b](./python/hdg_solution_step_1-8b.py), [step 8d](./python/hdg_solution_step_8d.py), [step 8e](./python/hdg_solution_step_8e.py), [step 9](./python/hdg_solution_step_9.py), [step 10](./python/hdg_solution_step_10.py)\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### notebook version\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[step 1-8b](./notebook/hdg_solution_step_1-8b.ipynb), [step 8d](./notebook/hdg_solution_step_8d.ipynb), [step 8e](./notebook/hdg_solution_step_8e.ipynb), [step 9](./notebook/hdg_solution_step_9.ipynb), [step 10](./notebook/hdg_solution_step_10.ipynb)\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that create_mesh and basis_function are (almost) self-contained.\n", "\n", "It could be interesting to develop / find a tool for mesh visualization.\n", "\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.5" }, "org": null }, "nbformat": 4, "nbformat_minor": 1 }