advection_1D.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a02f1044-cd45-45d7-81fc-00259bfb67f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Данный код решает скалярное линейное уравнение переноса -- самую простую модель сплошной среды,\n",
    "которая подчиняется уравнению в частных производных U_t + vel * U_x = 0, \n",
    "где vel - постоянная скорость потока (например, массы)\n",
    "несмотря на простоту модели, данная задача уже может указать на большинство проблем, \n",
    "которые могут встретиться при численном моделировании самых разных сплошных сред.\n",
    "Точное решение уравнения известно -- профиль, задаваемый начальным условием, сдвинется на расстояние \n",
    "S = time * vel за время time без какого-либо изменения\n",
    "'''\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import display, clear_output  # добавляем импорт clear_output и display\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c8791b9-dfbe-4e6f-b792-794df16a1bd6",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "#для использования метода конечного объема зададим одномерную сетку на призме единичного сечения \n",
    "#     +-----+-----+-----+-----+-----+-----+-----+-----+\n",
    "#    /  1  /     /     /     /     /     /     /     /|\n",
    "#   /     /     /     /     /     /     /     /     / |\n",
    "#  +-----+-----+-----+-----+-----+-----+-----+-----+  |\n",
    "#  | 1   |     |     |     |     |     |     |     |  |\n",
    "#  |     |     |     |     |     |     |     |     |  +\n",
    "#  |     |     |     |     |     |     |     |     | /\n",
    "#  |     |     |     |     |     |     |     |     |/\n",
    "#  +--dx-+-----+-----+-----+-----+-----+-----+-----+\n",
    "'''\n",
    "\n",
    "def grid_setup(xmin, xmax, Nx):\n",
    "    #наша расчетная область - это призма с объемом x1*x2*x3 = (x1max - x1min)*1*1, которую мы \"нарезаем\" вдоль x1\n",
    "    #количество фиктивных (\"ghost\") ячеек\n",
    "    Ngc = 2\n",
    "    #разрешение сетки dx1\n",
    "    dx = (xmax - xmin)/Nx\n",
    "    #задаем площадь грани и объем ячейки: fSurf = 1*1, сVol = 1*1*dx1\n",
    "    fSurf = 1.0\n",
    "    cVol = dx\n",
    "    #координата центра ячейки\n",
    "    cx = np.linspace(xmin - (Ngc - 0.5) * dx, xmax + (Ngc - 0.5) * dx, Nx + 2 * Ngc, dtype=np.double)\n",
    "    #координата грани ячейки\n",
    "    fx = np.linspace(xmin - Ngc * dx, xmax + Ngc * dx, Nx + 2 * Ngc + 1, dtype=np.double)\n",
    "\n",
    "    return Ngc, dx, fSurf, cVol, cx, fx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98173054-bd9d-40fc-a200-b2bdc5cc1257",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "зададим решение в момент времени time (setup = 'disc' имитирует ударную волну/контактный разрыв в газе, \n",
    "a setup = 'smooth' имитирует гладкое течение в газе)\n",
    "'''\n",
    "def solution(xmin, xmax, cx, Nx, Ngc, time, setup):\n",
    "\n",
    "    #зададим фоновую скорость потока\n",
    "    vel = 1.0\n",
    "\n",
    "    Len = xmax - xmin\n",
    "    \n",
    "    if setup == 'disc':\n",
    "        \n",
    "        #зададим исходные координаты разрыва\n",
    "        x0_1 = 0.2\n",
    "        x0_2 = 0.4\n",
    "        \n",
    "        #предполагаем что наша область периодическая, \n",
    "        #т.е. вещество, что выходит с одной стороны, сразу входит в область с другой стороны с сохранением параметров\n",
    "        x_1 = x0_1 + vel * time - np.floor((x0_1 + vel * time)/Len) * Len\n",
    "        x_2 = x0_2 + vel * time - np.floor((x0_2 + vel * time)/Len) * Len\n",
    "    \n",
    "        #задаем массив начальных данных в реальных ячейках\n",
    "        for i in range(Ngc, Nx + Ngc):\n",
    "            if x_2 > x_1:\n",
    "                if cx[i] < x_1:\n",
    "                    U0[i] = 1.0\n",
    "                elif cx[i] < x_2:\n",
    "                    U0[i] = 2.0\n",
    "                else:\n",
    "                    U0[i] = 1.0    \n",
    "            else:\n",
    "                U0[i] = 1.0\n",
    "                if cx[i] <= x_2:\n",
    "                    U0[i] = 2.0\n",
    "                if cx[i] > x_1:\n",
    "                    U0[i] = 2.0\n",
    "                    \n",
    "    elif setup == 'smooth':\n",
    "        \n",
    "        #начальная координата гладкого профиля\n",
    "        x0 = 0.3\n",
    "        \n",
    "        #зададим исходную толщину профиля\n",
    "        delta = 0.1\n",
    "    \n",
    "        #предполагаем что наша область периодическая, \n",
    "        #т.е. вещество, что выходит с одной стороны, сразу входит в область с другой стороны с сохранением параметров\n",
    "        x = x0 + vel * time - np.floor((x0 + vel * time)/Len) * Len\n",
    "    \n",
    "        #задаем массив начальных данных в реальных ячейках\n",
    "        for i in range(Ngc, Nx + Ngc):\n",
    "            U0[i] = np.exp(-(cx[i] - x)**2/delta**2) + np.exp(-(cx[i] - x - np.sign(vel) * Len)**2/delta**2)\n",
    "\n",
    "    \n",
    "    #возвращаем решение во время time\n",
    "    return U0, vel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "845e4905-3da1-4676-b59a-1dba54255202",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "ограничитель наклона для схем со вторым порядком точности \n",
    "обеспечивает монотонность решения вблизи разрывов без ухудшения его качества в областях гладкости\n",
    "-- не вносит в решение нефизичного осцилляторного поведения\n",
    "'''\n",
    "def limiter_func(x, y):\n",
    "\n",
    "    #smoothness analyzer \n",
    "    r = (y + 1e-14) / (x + 1e-14)\n",
    "    \n",
    "    #van Leer limiter  function\n",
    "    xy = x * (np.abs(r) + r) / (1.0 + np.abs(r))\n",
    "    \n",
    "    return xy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b844c3fd-d75a-4110-9f99-49021280828d",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "#здесь начинается основной код программы\n",
    "#-------------------------------------------------------------------------------------------------------------------\n",
    "#-------------------------------------------------------------------------------------------------------------------\n",
    "#-------------------------------------------------------------------------------------------------------------------\n",
    "#-------------------------------------------------------------------------------------------------------------------\n",
    "#-------------------------------------------------------------------------------------------------------------------\n",
    "'''\n",
    "#задаем нашу сетку - координаты границ и число ячеек\n",
    "xmin = 0.0\n",
    "xmax = 1.0\n",
    "Nx = 200\n",
    "Ngc, dx, fSurf, cVol, cx, fx = grid_setup(xmin,xmax,Nx)\n",
    "\n",
    "#инициализируем массивы решения U0 (текущий шаг по времени) и U1 (следующий шаг по времени), а также точное решение Ua\n",
    "U0 = np.zeros(Nx + Ngc * 2, dtype=np.double)\n",
    "U1 = np.zeros(Nx + Ngc * 2, dtype=np.double)\n",
    "Ua = np.zeros(Nx + Ngc * 2, dtype=np.double)\n",
    "\n",
    "#инициализируем массив потоков в реальных гранях\n",
    "flux = np.zeros(Nx+1, dtype=np.double)\n",
    "print(np.size(flux))\n",
    "print(Nx)\n",
    "#зададим параметр Куранта-Фридрихса-Леви\n",
    "CFL = 0.5\n",
    "\n",
    "\n",
    "#зададим вид начального условия\n",
    "setup = 'smooth'\n",
    "\n",
    "#зададим начальные условия для времени time = 0, а также финальное время time_fin\n",
    "U0, vel = solution(xmin, xmax, cx, Nx, Ngc, 0.0, setup)\n",
    "time_fin = 2.0\n",
    "\n",
    "#задаем начальное время равным нулю\n",
    "time = 0.0\n",
    "#задаем начальное число шагов\n",
    "nts = 0\n",
    "\n",
    "#если не работает ipython - можно поставить runtimeplot = 0\n",
    "runtimeplot = 1\n",
    "if runtimeplot == 1:\n",
    "    # Создание фигуры для графика\n",
    "    fig, ax = plt.subplots()\n",
    "    line, = ax.plot(cx[Ngc:-Ngc], U0[Ngc:-Ngc])\n",
    "    ax.set_title('sol at time = ' + str(np.round(time, 3)))\n",
    "    ax.set_xlabel('x1')\n",
    "    ax.set_ylabel('solution')\n",
    "    plt.close()  # Закрываем фигуру, чтобы она не отображалась сразу\n",
    "\n",
    "\n",
    "#цикл по времени\n",
    "while time < time_fin - 1e-14:\n",
    "    \n",
    "    #задаем шаг по времени \n",
    "    dt = min(CFL * dx / np.abs(vel), time_fin - time)\n",
    "    time += dt\n",
    "    nts += 1\n",
    "\n",
    "    #задаем граничные условия с использованием фиктивных ячеек (в данном примере они периодические)\n",
    "    for i in range(0,Ngc):\n",
    "        #внутренняя граница\n",
    "        U0[i] = U0[Nx+i]\n",
    "        #внешняя граница\n",
    "        U0[Nx+Ngc+i] = U0[Ngc+i]\n",
    "    \n",
    "    #записываем разные виды потока\n",
    "\n",
    "    #поток от величины слева \n",
    "    #flux = vel * U0[Ngc-1:Nx+Ngc]\n",
    "    \n",
    "    #поток от величины справа \n",
    "    #flux = vel * U0[Ngc:Nx+Ngc+1]\n",
    "    \n",
    "    #полусумма из двух соседних ячеек\n",
    "    #flux = vel * (U0[Ngc-1:Nx+Ngc] + U0[Ngc:Nx+Ngc+1])/2.0\n",
    "    \n",
    "    #поток из задачи Римана (cхема с первым порядком по времени и пространству)\n",
    "    #flux = vel * (U0[Ngc-1:Nx+Ngc] + U0[Ngc:Nx+Ngc+1])/2.0 + np.abs(vel) * (U0[Ngc-1:Nx+Ngc] - U0[Ngc:Nx+Ngc+1])/2.0\n",
    "    \n",
    "    #Схема Лакса-Вендроффа (второй порядок по времени и пространству)\n",
    "    flux = vel * (U0[Ngc-1:Nx+Ngc] + U0[Ngc:Nx+Ngc+1])/2.0 + vel * (vel * dt / dx) * (U0[Ngc-1:Nx+Ngc] - U0[Ngc:Nx+Ngc+1])/2.0\n",
    "    \n",
    "    #схема с ограничителем наклона (второй порядок по времени и пространству везде кроме областей с экстремумами/разрывами/сильными градиентами)\n",
    "    #Ul = U0[Ngc-1:Nx+Ngc] + 1.0/2.0 * (1.0 - vel*dt/dx) * limiter_func(U0[Ngc-1:Nx+Ngc] - U0[Ngc-2:Nx+Ngc-1], U0[Ngc:Nx+Ngc+1] - U0[Ngc-1:Nx+Ngc])\n",
    "    #Ur = U0[Ngc:Nx+Ngc+1] - 1.0/2.0 * (1.0 + vel*dt/dx) * limiter_func(U0[Ngc:Nx+Ngc+1] - U0[Ngc-1:Nx+Ngc], U0[Ngc+1:Nx+Ngc+2] - U0[Ngc:Nx+Ngc+1])\n",
    "    #flux = vel * (Ul + Ur) / 2.0 + np.abs(vel) * (Ul - Ur) / 2.0\n",
    "    \n",
    "    #обновляем наше численное решение\n",
    "    U1[Ngc:-Ngc] = U0[Ngc:-Ngc] - dt * ( flux[1:] * fSurf - flux[:-1] * fSurf ) / cVol\n",
    "    \n",
    "    #перезаписываем наше численное решение для использования на следующем шаге по времени (в итоге мы храним только два решения, а не все слои по t)\n",
    "    U0[Ngc:-Ngc] = U1[Ngc:-Ngc]\n",
    "\n",
    "    #вывод решения на экран\n",
    "    if ((nts%20 == 0 or time > time_fin - 1e-13) & (runtimeplot == 1)):\n",
    "        print('num of timesteps = ', nts)\n",
    "        line.set_data(cx[Ngc:-Ngc], U1[Ngc:-Ngc])\n",
    "        ax.set_title('sol at time = '+ str(np.round(time, 2)))\n",
    "        \n",
    "        # Автоматическое масштабирование осей\n",
    "        ax.relim()\n",
    "        ax.autoscale_view()\n",
    "        \n",
    "        clear_output(wait=True)\n",
    "        # Обновление графика с задержкой\n",
    "        plt.pause(0.1)\n",
    "        display(fig)\n",
    "\n",
    "\n",
    "plt.ioff()\n",
    "plt.show()\n",
    "\n",
    "print('end of simulation')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a02b5542-2f6c-49ce-b383-9a5878aa9953",
   "metadata": {},
   "outputs": [],
   "source": [
    "#точное решение \n",
    "Ua, vel = solution(xmin, xmax, cx, Nx, Ngc, time, setup)\n",
    "\n",
    "\n",
    "# Создаем новый график\n",
    "plt.figure(figsize=(8, 6))\n",
    "\n",
    "# Построение первой кривой (синяя линия)\n",
    "plt.plot(cx[Ngc:-Ngc], U1[Ngc:-Ngc], label='numerical', color='blue')\n",
    "\n",
    "# Построение второй кривой (красная линия)\n",
    "plt.plot(cx[Ngc:-Ngc], Ua[Ngc:-Ngc], label='analytical', color='red')\n",
    "\n",
    "# Добавление заголовка и меток осей\n",
    "plt.title('analytical vs numerical, time = ' + str(np.round(time, 2)))\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('solution')\n",
    "plt.grid(True)  # добавляем сетку\n",
    "\n",
    "# Добавляем легенду\n",
    "plt.legend()\n",
    "\n",
    "# Отображаем график\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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}