{ "cells": [ { "cell_type": "markdown", "id": "428e8785", "metadata": {}, "source": [ "## Résolution par différence finies de l'équation de Poisson sur un domaine 2D rectangulaire" ] }, { "cell_type": "markdown", "id": "5739a759", "metadata": {}, "source": [ "Position du Problème. On veut résoudre le problème de thermique suivant : \n", "\n", "(1) $ - \\lambda \\Delta u = f $, pour $(x,y) \\in \\Omega = ]0,L[ \\times ]0,H[$\n", "\n", "(2) $ u = g_d $, pour $(x,y) \\in [0,L] \\times \\{ 0 \\}$\n", "\n", "(3) $ u = g_u $, pour $(x,y) \\in [0,L] \\times \\{ H \\}$\n", "\n", "(4) $ u = g_l $, pour $(x,y) \\in \\{ 0 \\} \\times [0,H]$\n", "\n", "(5) $ u = g_r $, pour $(x,y) \\in \\{ L \\} \\times [0,H]$\n", "\n", "où $f, g_d, g_u, g_l, g_r$ sont des fonctions régulières données. \n", "\n", "Travail demandé : \n", " - Compléter les fonctions \"Matrice\" et \"Second_membre\" \n", " - Calculer la solution pour différentes conditions aux limites et différents termes source & vérifier si les résultats vous semblent physiquement cohérents.\n", " - Etudier l'influence du pas de maillage sur la solution\n", " - Modifier le programme pour résoudre le problème : \n", " \n", " (6) $ - \\epsilon^2 \\Delta u + u = f $, pour $(x,y) \\in \\Omega = ]0,L[ \\times ]0,H[$, $f$ fonction nulle sur $\\partial \\Omega$ \n", " \n", " (7) $u = 0$, pour $(x,y) \\in \\partial \\Omega$ \n", " \n", " \n", " - Que représente $\\epsilon$ ? Expliquer pourquoi la solution $u$ du problème est une version filtrée de $f$ (sans les plus hautes fréquences) ? \n", " " ] }, { "cell_type": "code", "execution_count": 2, "id": "4fdb55ec", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import numpy.linalg as npla\n", "import scipy as sp\n", "import scipy.sparse as spsp\n", "import scipy.sparse.linalg\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt\n", "from math import * " ] }, { "cell_type": "code", "execution_count": 103, "id": "46b480cc", "metadata": {}, "outputs": [], "source": [ "# Définition de la fonction f \n", "def funcf(x,y,L,H):\n", " xc = 0.8 * L \n", " yc = 0.8 * H\n", " delta_x = 0.1 * L \n", " delta_y = 0.1 * H\n", " fmax = 50000.0 \n", " f = fmax * exp(-((x-xc)/delta_x)**2 - ((y-yc)/delta_y)**2)\n", " return f" ] }, { "cell_type": "code", "execution_count": 104, "id": "97deda26", "metadata": {}, "outputs": [], "source": [ "# Définition des CL de Dirichlet \n", "def funcgb(x,L):\n", " # Bottom side\n", " gb = 0.0 # 20.0*sin(pi*x/L) # tempéraure imposée en °C \n", " return gb\n", "\n", "def funcgt(x,L):\n", " # Top side\n", " gt = 0.0 # tempéraure imposée en °C \n", " return gt\n", "\n", "def funcgl(y,H):\n", " # Left side\n", " gl = 0.0 # tempéraure imposée en °C \n", " return gl\n", "\n", "def funcgr(y,H):\n", " # Right side \n", " gr = 0.0 # tempéraure imposée en °C \n", " return gr" ] }, { "cell_type": "code", "execution_count": 105, "id": "925ad3ed", "metadata": {}, "outputs": [], "source": [ "def indk(i,j,I) :\n", " k = i + j*I\n", " return k " ] }, { "cell_type": "code", "execution_count": 106, "id": "870c0841", "metadata": {}, "outputs": [], "source": [ "def Mat_Laplacien_Dirichlet(I,J,dx,dy,lamda):\n", " K = I*J\n", " # On déclare une matrice creuse de dimension K x K\n", " A = spsp.csr_matrix((K, K))\n", " # On boucle sur les noeuds internes uniquement pour mettre à jour les coefficients A(k,k')\n", " for i in range(1,I-1):\n", " for j in range(1,J-1):\n", " k = indk(i,j,I)\n", " ke = indk(i+1,j,I)\n", " ko = indk(i-1,j,I)\n", " kn = indk(i,j+1,I)\n", " ks = indk(i,j-1,I)\n", " \n", " # A compléter \n", " \n", " # On met à 1 les coefficients diagonaux de A associés aux noeuds frontière pour gérer la CL de Dirichlet \n", " # Boucle sur les noeuds de la frontière x=0 \n", " for j in range(J):\n", " k = indk(0,j,I)\n", " A[k,k] = 1.0\n", " # Boucle sur les noeuds de la frontière x=L \n", " for j in range(J):\n", " k = indk(I-1,j,I)\n", " A[k,k] = 1.0\n", " # Boucle sur les noeuds de la frontière y=0 \n", " for i in range(I):\n", " k = indk(i,0,I)\n", " A[k,k] = 1.0\n", " # Boucle sur les noeuds de la frontière y=H \n", " for i in range(I):\n", " k = indk(i,J-1,I)\n", " A[k,k] = 1.0\n", " return A" ] }, { "cell_type": "code", "execution_count": 107, "id": "1173ef1d", "metadata": {}, "outputs": [], "source": [ "def Second_membre(I,J,dx,dy,L,H):\n", " K = I*J\n", " # Initialisation du second membre S \n", " S = np.zeros((K,1))\n", " # Boucle sur les noeuds internes pour la contribution de la source f à S \n", " for i in range(1,I-1):\n", " for j in range(1,J-1):\n", " k = indk(i,j,I)\n", " \n", " # A compléter \n", " \n", " # Boucle sur les noeuds de la frontière x=0 pour la contribution de gl à S\n", " for j in range(J):\n", " k = indk(0,j,I)\n", " # A compléter \n", " # Boucle sur les noeuds de la frontière x=L pour la contribution de gr à S \n", " for j in range(J):\n", " k = indk(I-1,j,I)\n", " # A compléter \n", " # Boucle sur les noeuds de la frontière y=0 pour la contribution de gb à S \n", " for i in range(I):\n", " k = indk(i,0,I)\n", " # A compléter \n", " # Boucle sur les noeuds de la frontière y=H pour la contribution de gt à S \n", " for i in range(I):\n", " k = indk(i,J-1,I)\n", " # A compléter \n", " return S " ] }, { "cell_type": "code", "execution_count": 108, "id": "c4f6f629", "metadata": {}, "outputs": [], "source": [ "def convert_u1d2d(u,I,J):\n", " u2d = np.zeros((I,J))\n", " for i in range(I):\n", " for j in range(J):\n", " k = indk(i,j,I)\n", " u2d[i,j] = u[k]\n", " return np.transpose(u2d)\n", "\n", "def calculf2d(I,J,dx,dy,L,K):\n", " f2d = np.zeros((I,J))\n", " for i in range(I):\n", " for j in range(J):\n", " f2d[i,j] = funcf(i*dx,j*dy,L,H)\n", " return np.transpose(f2d)" ] }, { "cell_type": "code", "execution_count": 109, "id": "13b341cd", "metadata": {}, "outputs": [], "source": [ "def flux_chaleur(u2d,lamda,I,J,dx,dy):\n", " # On doit prendre la transposée de u2d pour revenir à x = i et y = j\n", " u2d = np.transpose(u2d)\n", " Phi2dx = np.zeros((I,J))\n", " Phi2dy = np.zeros((I,J))\n", " # Boucle sur les noeuds internes \n", " for i in range(1,I-1):\n", " for j in range(1,J-1):\n", " Phi2dx[i,j] = - lamda * (u2d[i+1,j] - u2d[i-1,j])/(2.0 * dx)\n", " Phi2dy[i,j] = - lamda * (u2d[i,j+1] - u2d[i,j-1])/(2.0 * dy)\n", " # Boucle sur les noeuds de la frontière x=0 \n", " for j in range(1,J-1):\n", " i = 0\n", " Phi2dx[i,j] = - lamda * (u2d[i+1,j] - u2d[i,j])/dx\n", " Phi2dy[i,j] = - lamda * (u2d[i,j+1] - u2d[i,j-1])/(2.0 * dy)\n", " # Boucle sur les noeuds de la frontière x=L \n", " for j in range(1,J-1):\n", " i = I-1\n", " Phi2dx[i,j] = - lamda * (u2d[i,j] - u2d[i-1,j])/dx\n", " Phi2dy[i,j] = - lamda * (u2d[i,j+1] - u2d[i,j-1])/(2.0 * dy)\n", " # Boucle sur les noeuds de la frontière y=0 \n", " for i in range(1,I-1):\n", " j = 0\n", " Phi2dx[i,j] = - lamda * (u2d[i+1,j] - u2d[i-1,j])/(2.0 * dx)\n", " Phi2dy[i,j] = - lamda * (u2d[i,j+1] - u2d[i,j])/dy\n", " # Boucle sur les noeuds de la frontière y=H \n", " for i in range(1,I-1):\n", " j = J-1\n", " Phi2dx[i,j] = - lamda * (u2d[i+1,j] - u2d[i-1,j])/(2.0 * dx)\n", " Phi2dy[i,j] = - lamda * (u2d[i,j] - u2d[i,j-1])/dy\n", " # Cas des coins \n", " Phi2dx[0,0] = - lamda * (u2d[1,0] - u2d[0,0])/dx\n", " Phi2dy[0,0] = - lamda * (u2d[0,1] - u2d[0,0])/dy \n", " Phi2dx[0,J-1] = - lamda * (u2d[1,J-1] - u2d[0,J-1])/dx\n", " Phi2dy[0,J-1] = - lamda * (u2d[0,J-1] - u2d[0,J-2])/dy \n", " Phi2dx[I-1,0] = - lamda * (u2d[I-1,0] - u2d[I-2,0])/dx\n", " Phi2dy[I-1,0] = - lamda * (u2d[I-1,1] - u2d[I-1,0])/dy \n", " Phi2dx[I-1,J-1] = - lamda * (u2d[I-1,J-1] - u2d[I-2,J-1])/dx\n", " Phi2dy[I-1,J-1] = - lamda * (u2d[I-1,J-1] - u2d[I-1,J-2])/dy \n", " return np.transpose(Phi2dx), np.transpose(Phi2dy)" ] }, { "cell_type": "code", "execution_count": null, "id": "dcf469a7", "metadata": {}, "outputs": [], "source": [ "###################################### Programme principal #################################\n", "# Resolution de - lamda Lap(u) = f sur Omega, u = g sur dOmega \n", "\n", "# Données géométriques \n", "L = 1.0 \n", "H = 1.0 \n", "\n", "# Données physique\n", "lamda = 16.0 # Conductivité acier en W m^-1 K^-1 \n", "\n", "# Paramètres Numériques\n", "I = 51 # Nombre de noeuds selon x\n", "J = 51 # Nombre de noeuds selon y\n", "dx = L/(I-1) # Pas de maillage en x \n", "dy = H/(J-1) # Pas de maillage en y\n", "xn = np.linspace(0., L, I) # Abscisses des noeuds \n", "yn = np.linspace(0., H, J) # Ordonnées des noeuds \n", "X,Y = np.meshgrid(xn, yn) # Coordonnées des noeuds de la grille DF \n", "\n", "# Assemblage de la matrice associée au problème \n", "A = Matrice(I,J,dx,dy,lamda)\n", "\n", "# Assemblage du second membre associé au problème \n", "S = Second_membre(I,J,dx,dy,L,H)\n", "\n", "# Calcul de la solution du problème \n", "u = spsp.linalg.spsolve(A, S, permc_spec=None, use_umfpack=True)\n", "\n", "# Calculs auxillaires pour affichage de la solution \n", "u2d = convert_u1d2d(u,I,J) # solution en fonction de i,j\n", "Phi2dx, Phi2dy = flux_chaleur(u2d,lamda,I,J,dx,dy) # flux de chaleur \n", "f2d = calculf2d(I,J,dx,dy,L,H) # source de chaleur en fonction de i,j \n", "\n", "# Affichage \n", "fig, axs = plt.subplots(1, 2, figsize=(20, 7*H/L))\n", "axs = axs.flat\n", "Nblevels = 31\n", "\n", "## Iso-contour tempéraure u\n", "cs = axs[0].contourf(X, Y, u2d, Nblevels,cmap='hot' )\n", "fig.colorbar(cs, ax=axs[0], shrink=0.9)\n", "axs[0].set_title('Temperature')\n", "axs[0].set_xticks(np.linspace(0,L,5,endpoint=True))\n", "axs[0].set_yticks(np.linspace(0,H,5,endpoint=True))\n", "\n", "## Iso-niveaux source de chaleur f (source volumique)\n", "cs = axs[1].contourf(X, Y, f2d, Nblevels, cmap='hot')\n", "fig.colorbar(cs, ax=axs[1], shrink=0.9)\n", "axs[1].set_title('Heat Source')\n", "axs[1].set_xticks(np.linspace(0,L,5,endpoint=True))\n", "axs[1].set_yticks(np.linspace(0,H,5,endpoint=True))\n", " \n", "for ax in axs:\n", " ax.label_outer()\n", " fig.subplots_adjust(wspace=0.5, hspace=0.4)\n", " plt.show()\n", " \n", "fig, axs = plt.subplots(1, 2, figsize=(20, 7*H/L))\n", "axs = axs.flat\n", "Nblevels = 31\n", "\n", "## Iso-niveaux norme du flux de chaleur\n", "Phi2d = np.zeros((I,J))\n", "Phi2d[:,:] = (Phi2dx[:,:]**2 + Phi2dy[:,:]**2)**0.5\n", "cs = axs[0].contourf(X, Y, Phi2d, Nblevels, cmap='hot')\n", "fig.colorbar(cs, ax=axs[0], shrink=0.9)\n", "axs[0].set_title('Heat flux norm')\n", "axs[0].set_xticks(np.linspace(0,L,5,endpoint=True))\n", "axs[0].set_yticks(np.linspace(0,H,5,endpoint=True))\n", "\n", "## Lignes de courant associées au flux de chaleur \n", "cs = axs[1].streamplot(X, Y, Phi2dx, Phi2dy, linewidth=2)\n", "axs[1].set_title('Heat flux streamlines')\n", "\n", "for ax in axs:\n", " ax.label_outer()\n", " fig.subplots_adjust(wspace=0.5, hspace=0.4)\n", " plt.show()" ] }, { "cell_type": "markdown", "id": "dc45bb5d", "metadata": {}, "source": [ "A partir de cette cellule, on passe au problème de filtrage linéaire : \n", " \n", " $ - \\epsilon^2 \\Delta u + u = f $, pour $(x,y) \\in \\Omega = ]0,L[ \\times ]0,H[$, $f$ fonction nulle sur $\\partial \\Omega$\n", " \n", " \n", " $u = 0$, pour $(x,y) \\in \\partial \\Omega$ \n", " \n", " Il vous faut compléter les fonctions Matrice_bis et Second_membre_bis\n", " " ] }, { "cell_type": "code", "execution_count": null, "id": "4b0dddd3", "metadata": {}, "outputs": [], "source": [ "# Nouvelle définition de f pour le problème de filtrage linéaire \n", "\n", "def funcf_bis(x,y,L,H,kb): \n", " alpha = 0.05 \n", " f = sin(pi*x/L)*sin(pi*y/H) + alpha * bruit(x,y,L,H,kb)\n", " return f\n", "\n", "def bruit(x,y,L,H,kb):\n", " b = sin(kb*pi*x/L)*sin(kb*pi*y/H) \n", " return b" ] }, { "cell_type": "code", "execution_count": null, "id": "3e0400f5", "metadata": {}, "outputs": [], "source": [ "# Assemblage de la matrice de l'opérateur I - eps Laplacien \n", "def Matrice_bis(I,J,dx,dy,epsilon):\n", " K = I*J\n", " # On déclare une matrice creuse de dimension K x K\n", " A = spsp.csr_matrix((K, K))\n", " # On boucle sur les noeuds internes uniquement pour mettre à jour les coefficients A(k,k')\n", " for i in range(1,I-1):\n", " for j in range(1,J-1):\n", " k = indk(i,j,I)\n", " ke = indk(i+1,j,I)\n", " ko = indk(i-1,j,I)\n", " kn = indk(i,j+1,I)\n", " ks = indk(i,j-1,I)\n", " # A compléter \n", " # On rajoute l'identité à A \n", " # A compléter \n", " return A" ] }, { "cell_type": "code", "execution_count": null, "id": "b743c88a", "metadata": {}, "outputs": [], "source": [ "def Second_membre_bis(I,J,dx,dy,L,H,kb):\n", " K = I*J\n", " # Initialisation du second membre S \n", " S = np.zeros((K,1))\n", " # Boucle sur les noeuds internes pour la contribution de f à S \n", " # Pas besoin d'une boucle sur les noeuds frontière puisque dans ce pb u = 0 sur la frontière\n", " for i in range(1,I-1):\n", " for j in range(1,J-1):\n", " # A compléter \n", " return S " ] }, { "cell_type": "code", "execution_count": null, "id": "b9eca7e4", "metadata": {}, "outputs": [], "source": [ "def calculf2d_bis(I,J,dx,dy,L,K,kb):\n", " f2d = np.zeros((I,J))\n", " for i in range(I):\n", " for j in range(J):\n", " f2d[i,j] = funcf_bis(i*dx,j*dy,L,H,kb)\n", " return np.transpose(f2d)\n", "\n", "def calculuref2d(I,J,dx,dy,L,H):\n", " u2d = np.zeros((I,J))\n", " for i in range(I):\n", " for j in range(J):\n", " u2d[i,j] = sin(pi*i*dx/L)*sin(pi*j*dy/H)\n", " return np.transpose(u2d)" ] }, { "cell_type": "code", "execution_count": null, "id": "51a7c9f7", "metadata": {}, "outputs": [], "source": [ "###################################### Programme principal #################################\n", "# Resolution de - epsilon Lap(u) + u = f sur Omega, u = 0 sur dOmega \n", "\n", "\n", "# Données géométriques \n", "L = 1.0 \n", "H = 1.0 \n", "\n", "# Longueur de coupure \n", "epsilon = 0.1 \n", "# fréquence spatiale du bruit\n", "kb = 10 \n", "\n", "# Paramètres Numériques\n", "I = 51 # Nombre de noeuds selon x\n", "J = 51 # Nombre de noeuds selon y\n", "dx = L/(I-1) # Pas de maillage en x \n", "dy = H/(J-1) # Pas de maillage en y\n", "xn = np.linspace(0., L, I) # Abscisses des noeuds \n", "yn = np.linspace(0., H, J) # Ordonnées des noeuds \n", "X,Y = np.meshgrid(xn, yn) # Coordonnées des noeuds de la grille DF \n", "\n", "# Assemblage de la matrice associée au problème \n", "A = Matrice_bis(I,J,dx,dy,epsilon)\n", "\n", "# Assemblage du second membre associé au problème \n", "S = Second_membre_bis(I,J,dx,dy,L,H,kb)\n", "\n", "# Calcul de la solution du problème \n", "u = spsp.linalg.spsolve(A, S, permc_spec=None, use_umfpack=True)\n", "\n", "# Calculs auxillaires pour affichage de la solution \n", "u2d = convert_u1d2d(u,I,J) # solution u en fonction de i,j\n", "f2d = calculf2d_bis(I,J,dx,dy,L,H,kb) # fonction fbis en fonction de i,j \n", "uref2d = calculuref2d(I,J,dx,dy,L,H) # solution de reférence en fonction de i,j \n", "\n", "# Affichage \n", "fig, axs = plt.subplots(1, 3, figsize=(25, 6*H/L))\n", "axs = axs.flat\n", "Nblevels = 31\n", "\n", "## Iso-contour solution filtrée \n", "cs = axs[0].contourf(X, Y, u2d, Nblevels,cmap='hot' )\n", "fig.colorbar(cs, ax=axs[0], shrink=0.9)\n", "axs[0].set_title('Solution filtrée')\n", "axs[0].set_xticks(np.linspace(0,L,5,endpoint=True))\n", "axs[0].set_yticks(np.linspace(0,H,5,endpoint=True))\n", "\n", "## Iso-niveaux solution bruitée\n", "cs = axs[1].contourf(X, Y, f2d, Nblevels, cmap='hot')\n", "fig.colorbar(cs, ax=axs[1], shrink=0.9)\n", "axs[1].set_title('Solution bruitée')\n", "axs[1].set_xticks(np.linspace(0,L,5,endpoint=True))\n", "axs[1].set_yticks(np.linspace(0,H,5,endpoint=True))\n", "\n", "## Iso-niveaux solution idéale\n", "cs = axs[2].contourf(X, Y, uref2d, Nblevels, cmap='hot')\n", "fig.colorbar(cs, ax=axs[2], shrink=0.9)\n", "axs[2].set_title('Solution idéale')\n", "axs[2].set_xticks(np.linspace(0,L,5,endpoint=True))\n", "axs[2].set_yticks(np.linspace(0,H,5,endpoint=True))\n", " \n", "for ax in axs:\n", " ax.label_outer()\n", " fig.subplots_adjust(wspace=0.5, hspace=0.4)\n", " plt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.9.7" } }, "nbformat": 4, "nbformat_minor": 5 }