Diffusion_linear_solve.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9bdc387a-fc3e-496e-8fc9-c026d78754a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "## ME 3275 Linear solver Assignment 3\n",
    "## A list of physical parameters % L, k, q, TA, TB\n",
    "## A list of numerical paramters % N, dx, cell center location x(N)\n",
    "## Setting up A matrix: first interior points, then boundary points\n",
    "## Setting up b vector: first interior points, then boundary points\n",
    "## Solve for T\n",
    "## Plot T as a function of x \n",
    "## Define TDMA as a direct solver\n",
    "## Define Jacobi as the iterative solver using Jacobi algorithm\n",
    "## Define gSeidel as the iterative solver using Gauss Seidel algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8a6ff94-6ac8-47d2-b25c-d1dff5737ba4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68fc61b4-8273-4a9b-ac24-46d28924e7a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "L = 0.02; # m\n",
    "k = 0.5; # W/mK\n",
    "#q = 1000000; # W/m3\n",
    "TA = 100; #Celsius\n",
    "TB = 200; #Celsius \n",
    "## Numerical Paramteres\n",
    "N = 1000; \n",
    "dx = L/N; \n",
    "Xc = np.zeros(N);\n",
    "for i in range(0, N):\n",
    "        Xc[i] = i*dx + 0.5*dx\n",
    "# To facilitate method of manufectured Solution\n",
    "q=np.zeros(N);\n",
    "\n",
    "for i in range(0,N):\n",
    "    q[i] = 1000000; # uncomment for the original diffusion problem W/m3\n",
    "    #q[i] = k*(TB-TA)/np.sin(L)*np.sin(Xc[i])\n",
    "print(q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c05d61c-71b2-4de8-ae5f-b187e66d6648",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Set up the A matrix\n",
    "A = np.zeros((N,N));\n",
    "for i in range(1,N-1):\n",
    "    A[i,i] = -2*k/dx;\n",
    "    A[i,i+1] = k/dx;\n",
    "    A[i,i-1] = k/dx;\n",
    "# for left boundary\n",
    "A[0,0] = -3*k/dx;\n",
    "A[0,1] = k/dx;\n",
    "# for right boundary\n",
    "A[N-1,N-1] = -3*k/dx;\n",
    "A[N-1,N-2] = k/dx;\n",
    "#print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e538d840-41f7-4d23-942f-cabec42f6a64",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Set up b vector\n",
    "b=np.zeros(N);\n",
    "for i in range(1,N-1):\n",
    "    b[i] = -q[i]*dx;\n",
    "#left boundary\n",
    "    b[0] = -q[0]*dx - 2*k*TA/dx;\n",
    "#right boundary\n",
    "    b[N-1] = -q[N-1]*dx -2*k*TB/dx;\n",
    "#print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "753ac0ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Defining our function as Jacobi which takes 4 arguments \n",
    "# as A matrix, Solution x, relaxation factor alpha, and right-hand side b vector \n",
    "def Jacobi(A,x,alpha,b):\n",
    "    N = len(A) # find dimension of A matrix, assuming it is a square matrix\n",
    "    xold = x.copy(); # make sure to use old values from previous iteration\n",
    "    for i in range(0,N):\n",
    "        temp = b[i]\n",
    "        #for j in range(0,N):\n",
    "        for j in range(max(i-1,0),min(i+1,N)):\n",
    "            temp-=A[i][j]*xold[j]\n",
    "        x[i] = xold[i]+alpha*temp/A[i][i]\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "927c5044",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Defining our function as seidel which takes 4 arguments \n",
    "# as A matrix, Solution x, relaxation factor alpha, and right-hand side b vector \n",
    "def gSeidel(A,x,alpha,b):\n",
    "    N = len(A) # find dimension of A matrix, assuming it is a square matrix\n",
    "    for i in range(0,N):\n",
    "        temp = b[i]\n",
    "        #for j in range(max(i-1,0),min(i+1,N)): # you could do this; however, it is not a fair comparison when TDMA is utilizing the knowledge of tridiagnal matrix \n",
    "        for j in range(max(i-1,0),min(i+1,N)): # profit from the knowledg of A is a tridiagnal matrix\n",
    "            temp-=A[i][j]*x[j] # x is updated while i is looped through\n",
    "        x[i] = x[i]+alpha*temp/A[i][i]\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defining our function to take 3 arguments: A matrix, Solution x, and right-hand side b vector \n",
    "\n",
    "def TDMA(A,b):\n",
    "    N = len(A) # find dimension of A matrix, assuming it is a square matrix\n",
    "    P = np.zeros(N)\n",
    "    Q = np.zeros(N)\n",
    "    x = np.zeros(N)\n",
    "    P[0] = -A[0][1]/A[0][0]\n",
    "    Q[0] = b[0]/A[0][0] #assuming A(0,0) is non zero\n",
    "\n",
    "    for i in range(1,N-1):\n",
    "        P[i] = -A[i][i+1]/(A[i][i]+A[i][i-1]*P[i-1])\n",
    "        Q[i] = (-A[i][i-1]*Q[i-1]+b[i])/(A[i][i]+A[i][i-1]*P[i-1])\n",
    "    #backward substitute\n",
    "    x[N-1] = (-A[N-1][N-2]*Q[N-2]+b[N-1])/(A[N-1][N-1]+A[N-1][N-2]*P[N-2])\n",
    "  \n",
    "    for i in range(N-2,-1,-1):\n",
    "        x[i] = P[i]*x[i+1]+Q[i]\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8fa2f24-5c8f-4d35-91a9-092221ce3786",
   "metadata": {},
   "outputs": [],
   "source": [
    "# solve for Temperature\n",
    "# Initial guess\n",
    "T = 100*np.ones(N); # initial guess for temperature \n",
    "abstol=1;\n",
    "niter = 0;\n",
    "tic = time.process_time()\n",
    "#T = np.linalg.solve(A,b)\n",
    "#T=TDMA(A,b)\n",
    "\n",
    "while (abstol >1E-3):\n",
    "    Told = T.copy();\n",
    " #   T = Jacobi(A,T,1.0,b);\n",
    "    T = gSeidel(A,T,1.0,b);\n",
    "    abstol = max(abs(T-Told));\n",
    "    niter = niter + 1;\n",
    "#print(niter)\n",
    "toc = time.process_time()\n",
    "print(toc-tic)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e36641e2-0345-47f4-a577-e357c3b42cd7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# analytical solution for the original problem\n",
    "def analytical_solution(x,q,L,k,TA,TB):\n",
    "    return -q/2/k*x*x + (TB-TA+q*L*L/2/k)*x/L+TA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f03305f-1578-48f6-83eb-e3cb32ef3661",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the results \n",
    "analytical_results = analytical_solution(Xc,q,L,k,TA,TB) #uncomment for the original problem\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.plot(Xc, T, label='CFD results')\n",
    "plt.plot(Xc, analytical_results, label='Analytical') #uncomment for the original problem\n",
    "#plt.plot(Xc, (TB-TA)/np.sin(L)*np.sin(Xc)+ TA, label='Analytical')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('T')\n",
    "plt.title('Comparison of temperature numerical vs analytical')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "025bc292-bbb0-403c-b40b-1fcdf926ed75",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}