Navier–Stokes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\frac{\\partial \\vec{v}}{\\partial t}+(\\vec{v}\\cdot\\nabla)\\vec{v}=-\\frac{1}{\\rho}\\nabla p + \\nu \\nabla^2\\vec{v}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import init_printing, simplify, expand\n",
    "from sympy import symbols, Function\n",
    "from sympy.vector import Del, CoordSys3D, laplacian\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "C = CoordSys3D(\"C\")\n",
    "delop = Del()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y, z = C.x, C.y, C.z\n",
    "t, rho, nu = symbols(r\"t \\rho \\nu\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = Function(\"vx\")(t, x, y) * C.i + Function(\"vy\")(t, x, y) * C.j\n",
    "p = Function(\"p\")(t, x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "eq = v.diff(t) + v.dot(delop)(v) + delop(p) / rho - nu * laplacian(v)\n",
    "eq = eq.doit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [],
   "source": [
    "eqx = eq.coeff(C.i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "h = symbols(\"h\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import calculus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import Derivative"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqsAAAAsBAMAAACwH+yLAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAMmYiu80QdonvRN2ZVKvu110NAAAIc0lEQVRoBe1aTWxUVRQ+M9P56d8wMXHBwvQZAi40tiFoYkJCAxp1Y8efkiAqzx9CCIuOiQR2nQ0aApFGRQWCHTEujImtISQYBCYaRBdIZWN002pEFxgphLahWMZz7r3nvv95w+ugbxJvMu+ec+453/3umTf3vZ5egP/bv5SBhBEy0d0h45AwwjyCx0PBfUITho/RbgoHvcPufmtyOLjAWxWGmi+EePggZB4fDwlSw6HgPjA+8zm9QkFT/c4A0l455LX5WULBZdB6gMxRFb+tzw+o6Ge02ewIbE5XWdJ9NHAdbgn2+aKBpk0AT+Rz1gxK8rgIe1g6hFO2DNA9rXAmVe/sdjpVt+ZA4MElEyxl+pQUCZxRbL1jvmigtCJP5ALP0QzG7bjqrqJCPMPIjn6rQ/MoDgQeHTNZ6i4pKRI4o9h6x3zRQN9GOHdkxyzP0QzGnyJYZ1khHhd9tlYr5WrTygaQntCin+BAYIc9LIAmGQlcw1iCY75IoB3jiCYjLdgcGUVzMYbh2hzsq/FoUDqSm96CyyeX14yemQLAXoDPBp4oUVDHb9f+EMHtfwM8pmEgV7VkIQUiWH7PXXitT2qKZMPgFoaWll838jcuEO3DaFsEY8S5PgG0IE1HT5I+d886qbgZw2gV2spyDK+edKiR9ilIA5w305hVWIsf3mZ4j83WzMy7aN9+4MvVBkCyiLKjBSKwV2amLz8lFf7uGwZnDFv/AUAROgXtRTE+RTgAbYbtecLTDD0Nf0nZzRiGZ2GrCZDdcOjrr/zSIcOys7ASsYuvknoSP7zNdKn8ZXATqAI8NAfdowYyGEcfRwtEYK/uK9DWLxUmGQZucqxPP2ymDMBJ76KxRTAeNhMTiNCJF6ajZ5sswMNScTOGoTn6iaRGS/DTVVs6TB0rhMxVKONb1TsGaQfxw9uM3mPPj9Bz4XQBoAedEtPkaG+BCJtPUHuSSLeXZASTDAFPHTPtMzjlISOHhOeI9qIYD5UQB6C9Yj1PMhuJ8bEC3Vvv0SDmTTHXjLuuAz7kOq9gKjCtnA4P40uC32SFMDCtmXkSsA1VRAewZuF+vOevoZY2LBw1SF0QArsgozX4pUBucHDj/sHBKooh4J9spMdzQMuPUzq+y4j1LoJx18gWmoHSqumQgdpxSGHefBmnahemAIb7cfgpKx0exmNH+tDjLG2f9JPqmAaTRBiDjOghf+Nz3CXwwQUpE8fHUXC2IAT2QtK7FSh/9/XB8VbczcHevmN+GRp7jpg0tAjGifkVhECbgKZDBmx4b+VG5PI9jOHyiQLAQBXdHtTp8DLu/BMdcuYjeKVHVq7YjVHYdsEPoodErYiTz0nF55EFQQgqAnoB3k+aQmOSIeAmf6UM4ehPVVFNy7eTRTCGZ8oES48sTYcM2JLT0F4WSdebgOXSM4Mep0vCUafDwzi5gA47YKiK3Zu4e09tgSHMI3yfKWUIAGDThLpbSemq0tXRghDY6SLAvNjGrJ2qcXAGsfW9yAe65Wa1CMYwWSHQ3JR9rWQByBdhSaUsRL4RLMZpSpi4W7H3SQckDBzIjAPsmYHemRIAvlynHi1B5yzaV+3A78REAe7Fj9hbu3A96QpqzhaEwF54R+0tS4VJ1gGnolDbYY519QkDDduEsSiu0RnjUgQC7WoWHWGiOzj3rRS9jJMYIPfWczIdxDhxyJT+ePVUfzr75NhHyqUD9T74lTR6E2g38UVLDQV0bgSXG5MkcwB4HiealD8TVzCqFuOUIUbd8wWAepHQgi+W1I7KjiJ9mptxCvIGutGbAODziNJBjPO1EkqyrWeB+2RVSoYyJPHHdjMxQtrSq2YCn3zwhhoK6NwILrcOwzIEgRfx5eOm5eaQmHHKyJliwD1fEKgDRSgr4XVpvE92GOnXXIxhcmopuSVGV8CdhkoHMs7WKmSmli2Lzn5ZJ5SMoWwv4e59c/OE0LYfeKGAWwZC1G0uhDq+QeA7kbScUge/rCTNODetbjVwzRcEqqG0sC9lSLmtT/QYGdbIpff354VbdsP+XzgdyBh+1sH0lu9q8keWZKuBwqVnWaNecbCbnLILwTno0AzU/MC3OryE8qMyacapsyVlcs1noNkPVHnbutX47BAt1S86Q1zrXsglf6xi85HpcDLG7d7d+Bt027W+TEsBQihCQByZJXh6AmquTYDTelsY8x5bh1jgkGY8Viuz014W4tVjBeKBgLTGl3HqtE7r2nilk9lgjayL0qpKZmTmuzW+jGGgzIxPEuP4NayRUVq5ZEYEOa3xZUxpVYwPIuEltVtvGKZLPdEQguecFtULSiuXzOxpJcbBsYEjCHq7GVNaFWMiGcOGtTZMqy6Zwa7BwS8GB8W/QuPLGNPKjGP6k8K/JjGtumRGX3zMNwH6+3egzIzj+wDAtOqSmT2t8WWMaWXGWP2JY8OiEN+tTI/v1vgy5rsVGfu8XPM6/ss+XcHH0rQumREVTmt8GfPeikU+rv745bDRQ0h+sT62ho9hYayqkXHJjNA4rbFl/CKmlYt8XP3xSQOAePL6jkQypqsNh6kaGZfMKO6ICo4t41phdEIX+dbVWepCnbEIQ9YxrLBgXSOTJTOne1wZ1z68gUQVY1n9cRJXmnUIyXf4lo1jZqMhdWtkcWX88YlvrPWlDEt2STl8FWtm29MwmCwKBbi3HmPnQqxDSE57VM06hhUVISyuJRhbh5DCltPQuO0YVkP+EZxagvGkPoQUYYXeENsxLO9gcywtwfiMPoTUlEV3WcewmoLnA9ISjNUhJB/6kUx8DCtScGNBrcDYdgipsUWFeA1Zx7BCPKMOtwTjpHUIKeo6HXG2Y1gOe/OUlmBsO4TUlJVftI5hNQXPC9ISjG2HkLwriGDBQikfw4oQ3UhIbBj/A4HJ1/8iSADDAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$(\\frac{\\operatorname{vx}{\\left (t,t,\\mathbf{{y}_{C}} \\right )}}{h} - \\frac{\\operatorname{vx}{\\left (t,- h + t,\\mathbf{{y}_{C}} \\right )}}{h})\\mathbf{\\hat{i}_{C}} + (\\frac{\\operatorname{vy}{\\left (t,t,\\mathbf{{y}_{C}} \\right )}}{h} - \\frac{\\operatorname{vy}{\\left (t,- h + t,\\mathbf{{y}_{C}} \\right )}}{h})\\mathbf{\\hat{j}_{C}}$$"
      ],
      "text/plain": [
       "(vx(t, t, C.y)/h - vx(t, -h + t, C.y)/h)*C.i + (vy(t, t, C.y)/h - vy(t, -h + t, C.y)/h)*C.j"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Derivative.as_finite_difference(v.diff(x), [t-h, t])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "ders = list(eqx.find(Derivative))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [],
   "source": [
    "o = ders[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhUAAAAtBAMAAADrfISkAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMiKZu6uJRO92VGZ6zyUAAAAHP0lEQVRoBe1aXYwTVRQ+059tp912JxhjzBq3FiUqASqISpTQ8OALiVsfhAcSqCYKJiLrD/5EwxYx6oORBiJoNKGJmugLu0CIEYxtiEZNSLaJJj4o7IBKfDBhdwGBZWE85947f+10Z6Y788ZJeuec757zzZkzM3dm5yzADfFagce8OsJfrp7euRyoXOljqkOUVyhbtHuu3fSyHUArPtQGdQKSxU4zAvfBZWOS3iig7Uq/zhbk3bjlMPm+bQuQKrC/ZkPQiCgAN5VaUWc75wwbqB8uI4iUSJnGHA2zyB0A0vJZ5jtOLaCZiG06qUB61Iag8TP+xlvBDjb5ziZ+uGw8fU0yXeh7qwCZCVvc7MbN+vQ1UijclHQBkhdMk2tHcPNUK9jB7u+A67AfLj2GbUcU2rjQp0sAPTnm7234Q7ilJkmRRq1Rkcn2WqQK6LGCe41p1zLaJWvEmDYD+zUDiTQN1UnxxWUj2MMsF/rP0ClaZZ69mtYc14YsHIg0ZG3CggDotZALDF5km0SuScgu+Q52rTyjqYNTRQC5DKl3L73P3bYW4CHFFlEvQ6JqIOhslzlw2YgWnPqxxHKxoS2p7gX4cvi1BnPpOQ/RHTbn9H8Ar9sQoxaRLf9uwxmMt8kIHli6QuvIViVSxKmEarkJI5flIcQsMjYJ/YphZ3OGKpTuuaxM0lQpWQFwod+OIcbKtqSyz8qAy4GmSAftkH5dDNwF7+HMGfssPIh27yRswCLkvqe5aNNyE8brL9GB9y48/PTzNAkwMMPKKZBMgaPm2D2XyYEn4wIkhnAzO/1KDDFWtvWX7iOG48ceWKuSAhLeI2XcmohxXYwX4VWc+QV/FunBPYJ0Hao4fq3SRLpm3IRojdCKG6834PfrqKD0XAVcEHUkNsFRc+yeC6T5y1CWFmkvOUg3AFzoj6KnWNkAUloTzfUzkKmrqKBsHaXF1Yro1wUW8Ft0GGRuxvAT03ZSLWC8RgbVYoBpZPVfVbA0+KyJiVrEtVMVE2lPFrrmot3pgivixqJTLWz0WAtpWg+R6nRDrMKoQZVjG6/9Y0V25/MP5/PszWIFxOnxaa9FpgxPIDhyAAsIz7G7i+6REZDQJrlneAhgDH+wmNkAu5YVTYQ9KMSE2HTNZaXBk/EJKOBCvxKvhgl0YxLZhw+8XvxBROVI8vz9LYi4LrCA8igeIT2HTHkU4Ae0ou/gICsv0kRCBdgNv6KG0lPpvwIwXEb1cbJRBqdw0JH2xa17LsYuhnMA32QVh7XTRo9rp5zLFHnMZqlehQQ+OyAuqhPTcng4VkTUIjsB6erplrUzvuLjzQWMztKqcAIGyriRKwCbpYZEhyztV2JaDVY1UMc7J0djhHwFAj1lgmzSNZeV5S2AaRlPRdkKMt1K/xGdrNt4XrdOwuA0RGeYE0seYEkTTzNHOI+oRTIHfbUqwDNW9oimaVQLCYc9U3BuqoFXHerrTuDhKoBvVdWkppWGyywIX0VQshQgEIjUCLJJ11xWFjzje6tY95oVZLqVHq/x+JsNekUCWbsCu7SL7CpAP0oe4E/86QjZ+nMkoYL8LJr3EjarLOezqZLhxdaLLQBf4L4hqSKuI+sNH2fFD5cDgwt9VKSIeQlh60VPEyj5EnyAqI4wB3FdCOf4qFA6b87yqazpQc8ROAig4rOmsoZwHfmUjFnED5cDjQt9tsxjVDOUniNpXGjw3eRyjB2rQJjLAdMRNblhM52MBK/27eZcrH4aVqsgqQDndiwiXEdyZJBIPIgb5uiDywwyNCmnqx3ot7F5ykuXNdeV2GIATD57eV6TUIHoDpbtPIveQRXfX1TLdO/CQ/jhg66U5NIawzkiDhWRTIPBrYMPrtZQtF3p17EgysuQ48fuLrIrGHaK5YAjhoOp8EqatpO2wQl0xP420A61oHf77sWVPq52Tw5Z50t5DowitFMt5s7MGEKmDyjLG7VoK2TIJy48+nn09+CyO9sOaA5AeMmypEKmp33ga+ecBWnkfH7+oXy+jGrfnPnsBBPh0+MegpWQT1zI9Ddq4bsC8qYnvcT4P3F7X1G8EHOfkOk9JvIh4N/P7uI7WbmEfzl7lpDpPeaxHAa9vJ/5TjZail30mAK6hUzvmgjv7x6FjTVXV/zqoXpwEi6ss5uuxae9h4RM755IpMx8xn3c1+6k5LGAuWXYpyFvEb68wqDva7IUHvGViBdn+kqI30AabBP8EAY97++yfkugCfPOblsHL6h9hELP+7tng8rR4OGdXf2jlAEHpYRCz/q7qWqqGFSWgod3dr+CzwPmDZGe93dXn/wt6LWTdXZji06+EE4twqDn/d26Zv63RUC5j1NnN621/G9EQOT4eToEetHfDSxHg0h0dg07YCUMetHfDThTpBOd3eCJOWMY9AO8vxt4ynpnN3BiThgKvejvBp6y3tkNnJgThkKPf59Sfzdw0Tu7gRNzwlDot/P+buApJ0RnN3BiTuiP/n9c6UdAKMPy1QAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$- \\frac{2 \\operatorname{vx}{\\left (t,\\mathbf{{x}_{C}},\\mathbf{{y}_{C}} \\right )}}{h^{2}} + \\frac{\\operatorname{vx}{\\left (t,\\mathbf{{x}_{C}},\\mathbf{{y}_{C}} - h \\right )}}{h^{2}} + \\frac{\\operatorname{vx}{\\left (t,\\mathbf{{x}_{C}},\\mathbf{{y}_{C}} + h \\right )}}{h^{2}}$$"
      ],
      "text/plain": [
       "  2⋅vx(t, C_x, C_y)   vx(t, C_x, C_y - h)   vx(t, C_x, C_y + h)\n",
       "- ───────────────── + ─────────────────── + ───────────────────\n",
       "           2                    2                     2        \n",
       "          h                    h                     h         "
      ]
     },
     "execution_count": 114,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "o.as_finite_difference([y-h, y, y+h])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAK4AAAAzBAMAAAD82i5AAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAZokiELtEzTLvmd1Uq3Z1ltcVAAAD9UlEQVRYCa1XW4gbVRj+kpnNTq4NWihUpbHbvhSKQeqD9WHHFilFMOsFBaE26ItFhAgiLWzd+LT0KQtqLxRJRIsoitUKSr1sKvggFjZCwYIPO2/i0+626Wq7bU6/fyaZzCTZbdg9P5xzvv87/3zzz7kmQMh2vWyHfE1ONJ8uaZIKySTz1vUQocmJF41lTVK+zH3b88Sppk/oAdG8+SKVko4eOV9lC3CCzi6f0AQ+BzZlEalpkgvKJGt4Fa8HGR34vXp8wpq8/IIOrYDGdP3RQi2u1GKA0wBjp4F5R4NQj8Q4NaeKPaQGd94GzmjQ6ZX4DzCv9ZIb981bQEb3FmZaZotbOLfx/PoUDgHn+kgNxBOb99buIWNl7xHgd0fKPsSv93fxYPTYYHoQe3UQuRr3NqfhwmqdYX407K7pZXI89Yfc5RI7rMV5nSRKw0WbM8PFSdQbLO2lmFGqPqcmAs+ScaKBM2sy0NcHt9zMplcuTS1Z6ij7eJO8efiQ40YlVpC8EoqPc2M912XWvnY+AkpINCGyeJZlrvPkVGOnwH/OXNybdbmMss0PugxvtTVswTb41EuZssT8xPKZALHHb3wndQupqqdrchxqXeYPqFVskc8VslHWByVT4CTLty5iFVN11rNloOLp4usZmdgOU2Hv6pY+L7op76cKdU0PkTOr/OjMDaLRtu747YcCzNq6sWXZhIl/bdYyDrFFuJBqO6k5wrmC0WbSKz8EGFk8/RY996lH/lxju2OkJB5HI1pKlQUCX5nVHJI8BH2zFKN8ZvC8jeEvL35bnd/myOEM7GbmjWdQkJfsa6Ky7OXLMeGdQJuqt/MV5wuWyPYxgQG7gEredZ9mjrOliHqXHj/NeMVBsklOXcO8uu6Ob6KOWVuCH2DpMHiHzhieKrIJ2EmMFwOuB5Peq/Bat0dmP24jxp48/hS+zRgzTLcEi03Y5uywT6/zEyvb7TmxZFtvsYdr5o6n0WaiDgeNokvdWA/90kvQf97lzIAu99vxMvAkpe8cqLvdHnPAxcDFdttpZOv0mXeuM7lekzf9LSPqm5uC9SA+XFAtPKJ8fquPAkD29HAWyTPOmjR23EK1hpFc56lYLlbu4HW3R8pIt7DQxEGbC+XY2U8+xvTl34g3ZgY3zegyCi05YI2qg9+XUFXdIVmveprrPD6DxE28z334P8elb3GsS3pTiedhDoa61AAWJqjBFanBCjXK2fzx+2MZOEwHe1g2bgXHO64rcnrM0tFkIw0k5OtHb7Ny89UjbHyzXw41RM6zcsf3Sz3CqYfzXGBIZykn6wG8YbTZXEP+bMKqbsa0vECXbbsy6Upljp26qktTdNLfF3XKAXcBfwMT1vPJHQwAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\frac{\\partial^{2}}{\\partial \\mathbf{{y}_{C}}^{2}} \\operatorname{vx}{\\left (t,\\mathbf{{x}_{C}},\\mathbf{{y}_{C}} \\right )}$$"
      ],
      "text/plain": [
       "   2                  \n",
       "  ∂                   \n",
       "─────(vx(t, C_x, C_y))\n",
       "    2                 \n",
       "∂C_y                  "
      ]
     },
     "execution_count": 113,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "o"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}