morgan
/
blog
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
blog/public/posts/jeonsanyuceyeoghag-cfd-with.../index.html

776 lines
263 KiB
HTML
Raw Normal View History

Changed bunch
2023-10-21 16:05:44 +09:00
<!doctype html><html lang=en dir=auto><head><meta charset=utf-8><meta http-equiv=x-ua-compatible content="IE=edge"><meta name=viewport content="width=device-width,initial-scale=1,shrink-to-fit=no"><meta name=robots content="index, follow"><title>[전산유체역학] CFD with Python (Navier-Stokes Equation) | Morgan's Blog</title><meta name=keywords content><meta name=description content="1-D Linear Convection 1차원 선형 열전도 방정식은 가장 심플하면서도 가장 기초적인 방정식입니다.
$$ \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0 $$
이 식을 오일러 방정식으로 변환하여 수치해석적으로 해를 구할 수 있도록 변환을 해줍니다.
$$ u_i^{n+1} = u_i^n - c \frac{\Delta t}{\Delta x}(u_i^n-u_{i-1}^n) $$
이제 이 오일러 방정식을 파이썬으로 구현해봅니다.
import numpy from matplotlib import pyplot import time, sys %matplotlib inline nx = 41 # try changing this number from 41 to 81 and Run All ."><meta name=author content="Me"><link rel=canonical href=http://blog.morgan.kr/posts/jeonsanyuceyeoghag-cfd-with-python-navier-stokes-equation/><meta name=google-site-verification content="XYZabc"><meta name=yandex-verification content="XYZabc"><meta name=msvalidate.01 content="XYZabc"><link crossorigin=anonymous href=/assets/css/stylesheet.31527a12923607f33c1cac9636a2fa755f6ade7c55866bdb96e44c6bcaf6cfbb.css integrity="sha256-MVJ6EpI2B/M8HKyWNqL6dV9q3nxVhmvbluRMa8r2z7s=" rel="preload stylesheet" as=style><script defer crossorigin=anonymous src=/assets/js/highlight.f413e19d0714851f6474e7ee9632408e58ac146fbdbe62747134bea2fa3415e0.js integrity="sha256-9BPhnQcUhR9kdOfuljJAjlisFG+9vmJ0cTS+ovo0FeA=" onload=hljs.initHighlightingOnLoad()></script>
<link rel=icon href=https://blog.morgan.kr/favicon.ico><link rel=icon type=image/png sizes=16x16 href=http://blog.morgan.kr/favicon-16x16.png><link rel=icon type=image/png sizes=32x32 href=http://blog.morgan.kr/favicon-32x32.png><link rel=apple-touch-icon href=https://blog.morgan.kr/favicon.ico><link rel=mask-icon href=https://blog.morgan.kr/favicon.ico><meta name=theme-color content="#2e2e33"><meta name=msapplication-TileColor content="#2e2e33"><noscript><style>#theme-toggle,.top-link{display:none}</style><style>@media(prefers-color-scheme:dark){:root{--theme:rgb(29, 30, 32);--entry:rgb(46, 46, 51);--primary:rgb(218, 218, 219);--secondary:rgb(155, 156, 157);--tertiary:rgb(65, 66, 68);--content:rgb(196, 196, 197);--hljs-bg:rgb(46, 46, 51);--code-bg:rgb(55, 56, 62);--border:rgb(51, 51, 51)}.list{background:var(--theme)}.list:not(.dark)::-webkit-scrollbar-track{background:0 0}.list:not(.dark)::-webkit-scrollbar-thumb{border-color:var(--theme)}}</style></noscript><script type=application/javascript>var doNotTrack=!1;doNotTrack||(function(e,t,n,s,o,i,a){e.GoogleAnalyticsObject=o,e[o]=e[o]||function(){(e[o].q=e[o].q||[]).push(arguments)},e[o].l=1*new Date,i=t.createElement(n),a=t.getElementsByTagName(n)[0],i.async=1,i.src=s,a.parentNode.insertBefore(i,a)}(window,document,"script","https://www.google-analytics.com/analytics.js","ga"),ga("create","UA-123-45","auto"),ga("send","pageview"))</script><meta property="og:title" content="[전산유체역학] CFD with Python (Navier-Stokes Equation)"><meta property="og:description" content="1-D Linear Convection 1차원 선형 열전도 방정식은 가장 심플하면서도 가장 기초적인 방정식입니다.
$$ \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0 $$
이 식을 오일러 방정식으로 변환하여 수치해석적으로 해를 구할 수 있도록 변환을 해줍니다.
$$ u_i^{n+1} = u_i^n - c \frac{\Delta t}{\Delta x}(u_i^n-u_{i-1}^n) $$
이제 이 오일러 방정식을 파이썬으로 구현해봅니다.
import numpy from matplotlib import pyplot import time, sys %matplotlib inline nx = 41 # try changing this number from 41 to 81 and Run All ."><meta property="og:type" content="article"><meta property="og:url" content="http://blog.morgan.kr/posts/jeonsanyuceyeoghag-cfd-with-python-navier-stokes-equation/"><meta property="og:image" content="http://blog.morgan.kr"><meta property="article:section" content="posts"><meta property="article:published_time" content="2021-07-10T08:23:47+00:00"><meta property="article:modified_time" content="2021-07-10T08:23:47+00:00"><meta property="og:site_name" content="Morgan's Blog"><meta name=twitter:card content="summary_large_image"><meta name=twitter:image content="http://blog.morgan.kr"><meta name=twitter:title content="[전산유체역학] CFD with Python (Navier-Stokes Equation)"><meta name=twitter:description content="1-D Linear Convection 1차원 선형 열전도 방정식은 가장 심플하면서도 가장 기초적인 방정식입니다.
$$ \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0 $$
이 식을 오일러 방정식으로 변환하여 수치해석적으로 해를 구할 수 있도록 변환을 해줍니다.
$$ u_i^{n+1} = u_i^n - c \frac{\Delta t}{\Delta x}(u_i^n-u_{i-1}^n) $$
이제 이 오일러 방정식을 파이썬으로 구현해봅니다.
import numpy from matplotlib import pyplot import time, sys %matplotlib inline nx = 41 # try changing this number from 41 to 81 and Run All ."><script type=application/ld+json>{"@context":"https://schema.org","@type":"BreadcrumbList","itemListElement":[{"@type":"ListItem","position":2,"name":"Posts","item":"http://blog.morgan.kr/posts/"},{"@type":"ListItem","position":3,"name":"[전산유체역학] CFD with Python (Navier-Stokes Equation)","item":"http://blog.morgan.kr/posts/jeonsanyuceyeoghag-cfd-with-python-navier-stokes-equation/"}]}</script><script type=application/ld+json>{"@context":"https://schema.org","@type":"BlogPosting","headline":"[전산유체역학] CFD with Python (Navier-Stokes Equation)","name":"[전산유체역학] CFD with Python (Navier-Stokes Equation)","description":"1-D Linear Convection 1차원 선형 열전도 방정식은 가장 심플하면서도 가장 기초적인 방정식입니다.\n$$ \\frac{\\partial u}{\\partial t} + c \\frac{\\partial u}{\\partial x} = 0 $$\n이 식을 오일러 방정식으로 변환하여 수치해석적으로 해를 구할 수 있도록 변환을 해줍니다.\n$$ u_i^{n+1} = u_i^n - c \\frac{\\Delta t}{\\Delta x}(u_i^n-u_{i-1}^n) $$\n이제 이 오일러 방정식을 파이썬으로 구현해봅니다.\nimport numpy from matplotlib import pyplot import time, sys %matplotlib inline nx = 41 # try changing this number from 41 to 81 and Run All .","keywords":[],"articleBody":"1-D Linear Convection 1차원 선형 열전도 방정식은 가장 심플하면서도 가장 기초적인 방정식입니다.\n$$ \\frac{\\partial u}{\\partial t} + c \\frac{\\partial u}{\\partial x} = 0 $$\n이 식을 오일러 방정식으로 변환하여 수치해석적으로 해를 구할 수 있도록 변환을 해줍니다.\n$$ u_i^{n+1} = u_i^n - c \\frac{\\Delta t}{\\Delta x}(u_i^n-u_{i-1}^n) $$\n이제 이 오일러 방정식을 파이썬으로 구현해봅니다.\nimport numpy from matplotlib import pyplot import time, sys %matplotlib inline nx = 41 # try changing this number from 41 to 81 and Run All ... what happens? dx = 2 / (nx-1) nt = 25 #nt is the number of timesteps we want to calculate dt = .025 #dt is the amount of time each timestep covers (delta t) c = 1 #assume wavespeed of c = 1 u = numpy.ones(nx) #numpy function ones() u[int(.5 / dx):int(1 / dx + 1)] = 2 #setting u = 2 between 0.5 and 1 as per our I.C.s un = numpy.ones(nx) #initialize a temporary array for n in range(nt): #loop for values of n from 0 to nt, so it will run nt times un = u.copy() ##copy the existing values of u into un for i in range(1, nx): ## you can try commenting this line and... #for i in range(nx): ## ... uncommenting this line and see what happens! u[i] = un[i] - c * dt / dx * (un[i] - un[i-1]) pyplot.plot(numpy.linspace(0, 2, nx), u); 1-D Convection Equation (Non-Linear) $$ \\frac{\\partial u}{\\partial t} + u \\frac{\\partial u}{\\partial x} = 0 $$\n$$ u_i^{n+1} = u_i^n - u_i^n \\frac{\\Delta t}{\\Delta x} (u_i^n - u_{i-1}^n) $$\nimport numpy # we're importing numpy from matplotlib import pyplot # and our 2D plotting library %matplotlib inline nx = 41 dx = 2 / (nx - 1) nt = 20 #nt is the number of timesteps we want to calculate dt = .025 #dt is the amount of time each timestep covers (delta t) u = numpy.ones(nx) #as before, we initialize u with every value equal to 1. u[int(.5 / dx) : int(1 / dx + 1)] = 2 #then set u = 2 between 0.5 and 1 as per our I.C.s un = numpy.ones(nx) #initialize our placeholder array un, to hold the time-stepped solution for n in range(nt): #iterate through time un = u.copy() ##copy the existing values of u into un for i in range(1, nx): ##now we'll iterate through the u array ###This is the line from Step 1, copied exactly. Edit it for our new equation. ###then uncomment it and run the cell to evaluate Step 2 ###u[i] = un[i] - c * dt / dx * (un[i] - un[i-1]) pyplot.plot(numpy.linspace(0, 2, nx), u) ##Plot the results 1-D Diffusion Equation $$ \\frac{\\partial u}{\\partial t}= \\nu \\frac{\\partial^2 u}{\\partial x^2} $$\n$$ u_{i}^{n+1}=u_{i}^{n}+\\frac{\\nu\\Delta t}{\\Delta
MathJax.Hub.Config({
tex2jax: {
inlineMath: [['$','$'], ['\\(','\\)']],
displayMath: [['$$','$$']],
},
"HTML-CSS": {
scale: 80
},
});
</script><script src="//cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script><header class=header><nav class=nav><div class=logo><div class=logo-switches><button id=theme-toggle accesskey=t title="(Alt + T)"><svg id="moon" xmlns="http://www.w3.org/2000/svg" width="24" height="18" viewBox="0 0 24 24" fill="none" stroke="currentcolor" stroke-width="2" stroke-linecap="round" stroke-linejoin="round"><path d="M21 12.79A9 9 0 1111.21 3 7 7 0 0021 12.79z"/></svg><svg id="sun" xmlns="http://www.w3.org/2000/svg" width="24" height="18" viewBox="0 0 24 24" fill="none" stroke="currentcolor" stroke-width="2" stroke-linecap="round" stroke-linejoin="round"><circle cx="12" cy="12" r="5"/><line x1="12" y1="1" x2="12" y2="3"/><line x1="12" y1="21" x2="12" y2="23"/><line x1="4.22" y1="4.22" x2="5.64" y2="5.64"/><line x1="18.36" y1="18.36" x2="19.78" y2="19.78"/><line x1="1" y1="12" x2="3" y2="12"/><line x1="21" y1="12" x2="23" y2="12"/><line x1="4.22" y1="19.78" x2="5.64" y2="18.36"/><line x1="18.36" y1="5.64" x2="19.78" y2="4.22"/></svg></button></div></div><ul id=menu><li><a href=http://blog.morgan.kr/categories/ title=Categories><span>Categories</span></a></li><li><a href=http://blog.morgan.kr/tags/ title=Tags><span>Tags</span></a></li><li><a href=http://blog.morgan.kr/posts/ title=Posts><span>Posts</span></a></li></ul></nav></header><main class=main><article class=post-single><header class=post-header><div class=breadcrumbs><a href=http://blog.morgan.kr>Home</a>&nbsp;»&nbsp;<a href=http://blog.morgan.kr/posts/>Posts</a></div><h1 class=post-title>[전산유체역학] CFD with Python (Navier-Stokes Equation)</h1><div class=post-meta><span title='2021-07-10 08:23:47 +0000 UTC'>July 10, 10000</span>&nbsp;·&nbsp;4414 words&nbsp;·&nbsp;Me</div></header><div class=post-content><h2 id=1-d-linear-convection>1-D Linear Convection<a hidden class=anchor aria-hidden=true href=#1-d-linear-convection>#</a></h2><p>1차원 선형 열전도 방정식은 가장 심플하면서도 가장 기초적인 방정식입니다.</p><p>$$ \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0 $$</p><p>이 식을 오일러 방정식으로 변환하여 수치해석적으로 해를 구할 수 있도록 변환을 해줍니다.</p><p>$$ u_i^{n+1} = u_i^n - c \frac{\Delta t}{\Delta x}(u_i^n-u_{i-1}^n) $$</p><p>이제 이 오일러 방정식을 파이썬으로 구현해봅니다.</p><div class=highlight><pre tabindex=0 class=chroma><code class=language-python data-lang=python><span class=line><span class=cl><span class=kn>import</span> <span class=nn>numpy</span>
</span></span><span class=line><span class=cl><span class=kn>from</span> <span class=nn>matplotlib</span> <span class=kn>import</span> <span class=n>pyplot</span>
</span></span><span class=line><span class=cl><span class=kn>import</span> <span class=nn>time</span><span class=o>,</span> <span class=nn>sys</span>
</span></span><span class=line><span class=cl><span class=o>%</span><span class=n>matplotlib</span> <span class=n>inline</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>nx</span> <span class=o>=</span> <span class=mi>41</span> <span class=c1># try changing this number from 41 to 81 and Run All ... what happens?</span>
</span></span><span class=line><span class=cl><span class=n>dx</span> <span class=o>=</span> <span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=n>nx</span><span class=o>-</span><span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>nt</span> <span class=o>=</span> <span class=mi>25</span> <span class=c1>#nt is the number of timesteps we want to calculate</span>
</span></span><span class=line><span class=cl><span class=n>dt</span> <span class=o>=</span> <span class=mf>.025</span> <span class=c1>#dt is the amount of time each timestep covers (delta t)</span>
</span></span><span class=line><span class=cl><span class=n>c</span> <span class=o>=</span> <span class=mi>1</span> <span class=c1>#assume wavespeed of c = 1</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>u</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>ones</span><span class=p>(</span><span class=n>nx</span><span class=p>)</span> <span class=c1>#numpy function ones()</span>
</span></span><span class=line><span class=cl><span class=n>u</span><span class=p>[</span><span class=nb>int</span><span class=p>(</span><span class=mf>.5</span> <span class=o>/</span> <span class=n>dx</span><span class=p>):</span><span class=nb>int</span><span class=p>(</span><span class=mi>1</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>+</span> <span class=mi>1</span><span class=p>)]</span> <span class=o>=</span> <span class=mi>2</span> <span class=c1>#setting u = 2 between 0.5 and 1 as per our I.C.s</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>un</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>ones</span><span class=p>(</span><span class=n>nx</span><span class=p>)</span> <span class=c1>#initialize a temporary array</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=k>for</span> <span class=n>n</span> <span class=ow>in</span> <span class=nb>range</span><span class=p>(</span><span class=n>nt</span><span class=p>):</span> <span class=c1>#loop for values of n from 0 to nt, so it will run nt times</span>
</span></span><span class=line><span class=cl> <span class=n>un</span> <span class=o>=</span> <span class=n>u</span><span class=o>.</span><span class=n>copy</span><span class=p>()</span> <span class=c1>##copy the existing values of u into un</span>
</span></span><span class=line><span class=cl> <span class=k>for</span> <span class=n>i</span> <span class=ow>in</span> <span class=nb>range</span><span class=p>(</span><span class=mi>1</span><span class=p>,</span> <span class=n>nx</span><span class=p>):</span> <span class=c1>## you can try commenting this line and...</span>
</span></span><span class=line><span class=cl> <span class=c1>#for i in range(nx): ## ... uncommenting this line and see what happens!</span>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[</span><span class=n>i</span><span class=p>]</span> <span class=o>=</span> <span class=n>un</span><span class=p>[</span><span class=n>i</span><span class=p>]</span> <span class=o>-</span> <span class=n>c</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>*</span> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=n>i</span><span class=p>]</span> <span class=o>-</span> <span class=n>un</span><span class=p>[</span><span class=n>i</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>plot</span><span class=p>(</span><span class=n>numpy</span><span class=o>.</span><span class=n>linspace</span><span class=p>(</span><span class=mi>0</span><span class=p>,</span> <span class=mi>2</span><span class=p>,</span> <span class=n>nx</span><span class=p>),</span> <span class=n>u</span><span class=p>);</span>
</span></span></code></pre></div><figure><img loading=lazy src=https://blog.kakaocdn.net/dn/ZbM8j/btq9fWovXzY/D1HOqkCqgkw9YLDpyMFxb1/img.png></figure><h2 id=1-d-convection-equation-non-linear>1-D Convection Equation (Non-Linear)<a hidden class=anchor aria-hidden=true href=#1-d-convection-equation-non-linear>#</a></h2><p>$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = 0 $$</p><p>$$ u_i^{n+1} = u_i^n - u_i^n \frac{\Delta t}{\Delta x} (u_i^n - u_{i-1}^n) $$</p><div class=highlight><pre tabindex=0 class=chroma><code class=language-python data-lang=python><span class=line><span class=cl><span class=kn>import</span> <span class=nn>numpy</span> <span class=c1># we&#39;re importing numpy </span>
</span></span><span class=line><span class=cl><span class=kn>from</span> <span class=nn>matplotlib</span> <span class=kn>import</span> <span class=n>pyplot</span> <span class=c1># and our 2D plotting library</span>
</span></span><span class=line><span class=cl><span class=o>%</span><span class=n>matplotlib</span> <span class=n>inline</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>nx</span> <span class=o>=</span> <span class=mi>41</span>
</span></span><span class=line><span class=cl><span class=n>dx</span> <span class=o>=</span> <span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=n>nx</span> <span class=o>-</span> <span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>nt</span> <span class=o>=</span> <span class=mi>20</span> <span class=c1>#nt is the number of timesteps we want to calculate</span>
</span></span><span class=line><span class=cl><span class=n>dt</span> <span class=o>=</span> <span class=mf>.025</span> <span class=c1>#dt is the amount of time each timestep covers (delta t)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>u</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>ones</span><span class=p>(</span><span class=n>nx</span><span class=p>)</span> <span class=c1>#as before, we initialize u with every value equal to 1.</span>
</span></span><span class=line><span class=cl><span class=n>u</span><span class=p>[</span><span class=nb>int</span><span class=p>(</span><span class=mf>.5</span> <span class=o>/</span> <span class=n>dx</span><span class=p>)</span> <span class=p>:</span> <span class=nb>int</span><span class=p>(</span><span class=mi>1</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>+</span> <span class=mi>1</span><span class=p>)]</span> <span class=o>=</span> <span class=mi>2</span> <span class=c1>#then set u = 2 between 0.5 and 1 as per our I.C.s</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>un</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>ones</span><span class=p>(</span><span class=n>nx</span><span class=p>)</span> <span class=c1>#initialize our placeholder array un, to hold the time-stepped solution</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=k>for</span> <span class=n>n</span> <span class=ow>in</span> <span class=nb>range</span><span class=p>(</span><span class=n>nt</span><span class=p>):</span> <span class=c1>#iterate through time</span>
</span></span><span class=line><span class=cl> <span class=n>un</span> <span class=o>=</span> <span class=n>u</span><span class=o>.</span><span class=n>copy</span><span class=p>()</span> <span class=c1>##copy the existing values of u into un</span>
</span></span><span class=line><span class=cl> <span class=k>for</span> <span class=n>i</span> <span class=ow>in</span> <span class=nb>range</span><span class=p>(</span><span class=mi>1</span><span class=p>,</span> <span class=n>nx</span><span class=p>):</span> <span class=c1>##now we&#39;ll iterate through the u array</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=c1>###This is the line from Step 1, copied exactly. Edit it for our new equation.</span>
</span></span><span class=line><span class=cl> <span class=c1>###then uncomment it and run the cell to evaluate Step 2 </span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=c1>###u[i] = un[i] - c * dt / dx * (un[i] - un[i-1]) </span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>plot</span><span class=p>(</span><span class=n>numpy</span><span class=o>.</span><span class=n>linspace</span><span class=p>(</span><span class=mi>0</span><span class=p>,</span> <span class=mi>2</span><span class=p>,</span> <span class=n>nx</span><span class=p>),</span> <span class=n>u</span><span class=p>)</span> <span class=c1>##Plot the results</span>
</span></span></code></pre></div><h2 id=1-d-diffusion-equation>1-D Diffusion Equation<a hidden class=anchor aria-hidden=true href=#1-d-diffusion-equation>#</a></h2><p>$$ \frac{\partial u}{\partial t}= \nu \frac{\partial^2 u}{\partial x^2} $$</p><p>$$ u_{i}^{n+1}=u_{i}^{n}+\frac{\nu\Delta t}{\Delta x^2}(u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}) $$</p><div class=highlight><pre tabindex=0 class=chroma><code class=language-python data-lang=python><span class=line><span class=cl><span class=kn>import</span> <span class=nn>numpy</span> <span class=c1>#loading our favorite library</span>
</span></span><span class=line><span class=cl><span class=kn>from</span> <span class=nn>matplotlib</span> <span class=kn>import</span> <span class=n>pyplot</span> <span class=c1>#and the useful plotting library</span>
</span></span><span class=line><span class=cl><span class=o>%</span><span class=n>matplotlib</span> <span class=n>inline</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>nx</span> <span class=o>=</span> <span class=mi>41</span>
</span></span><span class=line><span class=cl><span class=n>dx</span> <span class=o>=</span> <span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=n>nx</span> <span class=o>-</span> <span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>nt</span> <span class=o>=</span> <span class=mi>20</span> <span class=c1>#the number of timesteps we want to calculate</span>
</span></span><span class=line><span class=cl><span class=n>nu</span> <span class=o>=</span> <span class=mf>0.3</span> <span class=c1>#the value of viscosity</span>
</span></span><span class=line><span class=cl><span class=n>sigma</span> <span class=o>=</span> <span class=mf>.2</span> <span class=c1>#sigma is a parameter, we&#39;ll learn more about it later</span>
</span></span><span class=line><span class=cl><span class=n>dt</span> <span class=o>=</span> <span class=n>sigma</span> <span class=o>*</span> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>/</span> <span class=n>nu</span> <span class=c1>#dt is defined using sigma ... more later!</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>u</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>ones</span><span class=p>(</span><span class=n>nx</span><span class=p>)</span> <span class=c1>#a numpy array with nx elements all equal to 1.</span>
</span></span><span class=line><span class=cl><span class=n>u</span><span class=p>[</span><span class=nb>int</span><span class=p>(</span><span class=mf>.5</span> <span class=o>/</span> <span class=n>dx</span><span class=p>):</span><span class=nb>int</span><span class=p>(</span><span class=mi>1</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>+</span> <span class=mi>1</span><span class=p>)]</span> <span class=o>=</span> <span class=mi>2</span> <span class=c1>#setting u = 2 between 0.5 and 1 as per our I.C.s</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>un</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>ones</span><span class=p>(</span><span class=n>nx</span><span class=p>)</span> <span class=c1>#our placeholder array, un, to advance the solution in time</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=k>for</span> <span class=n>n</span> <span class=ow>in</span> <span class=nb>range</span><span class=p>(</span><span class=n>nt</span><span class=p>):</span> <span class=c1>#iterate through time</span>
</span></span><span class=line><span class=cl> <span class=n>un</span> <span class=o>=</span> <span class=n>u</span><span class=o>.</span><span class=n>copy</span><span class=p>()</span> <span class=c1>##copy the existing values of u into un</span>
</span></span><span class=line><span class=cl> <span class=k>for</span> <span class=n>i</span> <span class=ow>in</span> <span class=nb>range</span><span class=p>(</span><span class=mi>1</span><span class=p>,</span> <span class=n>nx</span> <span class=o>-</span> <span class=mi>1</span><span class=p>):</span>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[</span><span class=n>i</span><span class=p>]</span> <span class=o>=</span> <span class=n>un</span><span class=p>[</span><span class=n>i</span><span class=p>]</span> <span class=o>+</span> <span class=n>nu</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=n>i</span><span class=o>+</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>un</span><span class=p>[</span><span class=n>i</span><span class=p>]</span> <span class=o>+</span> <span class=n>un</span><span class=p>[</span><span class=n>i</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>plot</span><span class=p>(</span><span class=n>numpy</span><span class=o>.</span><span class=n>linspace</span><span class=p>(</span><span class=mi>0</span><span class=p>,</span> <span class=mi>2</span><span class=p>,</span> <span class=n>nx</span><span class=p>),</span> <span class=n>u</span><span class=p>);</span>
</span></span></code></pre></div><figure><img loading=lazy src=https://blog.kakaocdn.net/dn/oaKgN/btq9iA57vBH/weCFyYoImjFkasFiDMir3k/img.png></figure><h2 id=burgers-equation>Burger&rsquo;s Equation<a hidden class=anchor aria-hidden=true href=#burgers-equation>#</a></h2><p>$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = \nu \frac{\partial ^2u}{\partial x^2} $$</p><p>$$ u_i^{n+1} = u_i^n - u_i^n \frac{\Delta t}{\Delta x} (u_i^n - u_{i-1}^n) + \nu \frac{\Delta t}{\Delta x^2}(u_{i+1}^n - 2u_i^n + u_{i-1}^n) $$</p><div class=highlight><pre tabindex=0 class=chroma><code class=language-python data-lang=python><span class=line><span class=cl><span class=kn>import</span> <span class=nn>numpy</span>
</span></span><span class=line><span class=cl><span class=kn>import</span> <span class=nn>sympy</span>
</span></span><span class=line><span class=cl><span class=kn>from</span> <span class=nn>sympy</span> <span class=kn>import</span> <span class=n>init_printing</span>
</span></span><span class=line><span class=cl><span class=kn>from</span> <span class=nn>matplotlib</span> <span class=kn>import</span> <span class=n>pyplot</span>
</span></span><span class=line><span class=cl><span class=kn>from</span> <span class=nn>sympy.utilities.lambdify</span> <span class=kn>import</span> <span class=n>lambdify</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=o>%</span><span class=n>matplotlib</span> <span class=n>inline</span>
</span></span><span class=line><span class=cl><span class=n>init_printing</span><span class=p>(</span><span class=n>use_latex</span><span class=o>=</span><span class=kc>True</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>x</span><span class=p>,</span> <span class=n>nu</span><span class=p>,</span> <span class=n>t</span> <span class=o>=</span> <span class=n>sympy</span><span class=o>.</span><span class=n>symbols</span><span class=p>(</span><span class=s1>&#39;x nu t&#39;</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>phi</span> <span class=o>=</span> <span class=p>(</span><span class=n>sympy</span><span class=o>.</span><span class=n>exp</span><span class=p>(</span><span class=o>-</span><span class=p>(</span><span class=n>x</span> <span class=o>-</span> <span class=mi>4</span> <span class=o>*</span> <span class=n>t</span><span class=p>)</span><span class=o>**</span><span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=mi>4</span> <span class=o>*</span> <span class=n>nu</span> <span class=o>*</span> <span class=p>(</span><span class=n>t</span> <span class=o>+</span> <span class=mi>1</span><span class=p>)))</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>sympy</span><span class=o>.</span><span class=n>exp</span><span class=p>(</span><span class=o>-</span><span class=p>(</span><span class=n>x</span> <span class=o>-</span> <span class=mi>4</span> <span class=o>*</span> <span class=n>t</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>sympy</span><span class=o>.</span><span class=n>pi</span><span class=p>)</span><span class=o>**</span><span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=mi>4</span> <span class=o>*</span> <span class=n>nu</span> <span class=o>*</span> <span class=p>(</span><span class=n>t</span> <span class=o>+</span> <span class=mi>1</span><span class=p>))))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>phiprime</span> <span class=o>=</span> <span class=n>phi</span><span class=o>.</span><span class=n>diff</span><span class=p>(</span><span class=n>x</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>u</span> <span class=o>=</span> <span class=o>-</span><span class=mi>2</span> <span class=o>*</span> <span class=n>nu</span> <span class=o>*</span> <span class=p>(</span><span class=n>phiprime</span> <span class=o>/</span> <span class=n>phi</span><span class=p>)</span> <span class=o>+</span> <span class=mi>4</span>
</span></span><span class=line><span class=cl><span class=n>ufunc</span> <span class=o>=</span> <span class=n>lambdify</span><span class=p>((</span><span class=n>t</span><span class=p>,</span> <span class=n>x</span><span class=p>,</span> <span class=n>nu</span><span class=p>),</span> <span class=n>u</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=c1>###variable declarations</span>
</span></span><span class=line><span class=cl><span class=n>nx</span> <span class=o>=</span> <span class=mi>101</span>
</span></span><span class=line><span class=cl><span class=n>nt</span> <span class=o>=</span> <span class=mi>100</span>
</span></span><span class=line><span class=cl><span class=n>dx</span> <span class=o>=</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>numpy</span><span class=o>.</span><span class=n>pi</span> <span class=o>/</span> <span class=p>(</span><span class=n>nx</span> <span class=o>-</span> <span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>nu</span> <span class=o>=</span> <span class=mf>.07</span>
</span></span><span class=line><span class=cl><span class=n>dt</span> <span class=o>=</span> <span class=n>dx</span> <span class=o>*</span> <span class=n>nu</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>x</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>linspace</span><span class=p>(</span><span class=mi>0</span><span class=p>,</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>numpy</span><span class=o>.</span><span class=n>pi</span><span class=p>,</span> <span class=n>nx</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>un</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>empty</span><span class=p>(</span><span class=n>nx</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>t</span> <span class=o>=</span> <span class=mi>0</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>u</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>asarray</span><span class=p>([</span><span class=n>ufunc</span><span class=p>(</span><span class=n>t</span><span class=p>,</span> <span class=n>x0</span><span class=p>,</span> <span class=n>nu</span><span class=p>)</span> <span class=k>for</span> <span class=n>x0</span> <span class=ow>in</span> <span class=n>x</span><span class=p>])</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=k>for</span> <span class=n>n</span> <span class=ow>in</span> <span class=nb>range</span><span class=p>(</span><span class=n>nt</span><span class=p>):</span>
</span></span><span class=line><span class=cl> <span class=n>un</span> <span class=o>=</span> <span class=n>u</span><span class=o>.</span><span class=n>copy</span><span class=p>()</span>
</span></span><span class=line><span class=cl> <span class=k>for</span> <span class=n>i</span> <span class=ow>in</span> <span class=nb>range</span><span class=p>(</span><span class=mi>1</span><span class=p>,</span> <span class=n>nx</span><span class=o>-</span><span class=mi>1</span><span class=p>):</span>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[</span><span class=n>i</span><span class=p>]</span> <span class=o>=</span> <span class=n>un</span><span class=p>[</span><span class=n>i</span><span class=p>]</span> <span class=o>-</span> <span class=n>un</span><span class=p>[</span><span class=n>i</span><span class=p>]</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>*</span><span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=n>i</span><span class=p>]</span> <span class=o>-</span> <span class=n>un</span><span class=p>[</span><span class=n>i</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>+</span> <span class=n>nu</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>\
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=n>i</span><span class=o>+</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>un</span><span class=p>[</span><span class=n>i</span><span class=p>]</span> <span class=o>+</span> <span class=n>un</span><span class=p>[</span><span class=n>i</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[</span><span class=mi>0</span><span class=p>]</span> <span class=o>=</span> <span class=n>un</span><span class=p>[</span><span class=mi>0</span><span class=p>]</span> <span class=o>-</span> <span class=n>un</span><span class=p>[</span><span class=mi>0</span><span class=p>]</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>*</span> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>0</span><span class=p>]</span> <span class=o>-</span> <span class=n>un</span><span class=p>[</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>+</span> <span class=n>nu</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>\
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>un</span><span class=p>[</span><span class=mi>0</span><span class=p>]</span> <span class=o>+</span> <span class=n>un</span><span class=p>[</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=n>u</span><span class=p>[</span><span class=mi>0</span><span class=p>]</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>u_analytical</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>asarray</span><span class=p>([</span><span class=n>ufunc</span><span class=p>(</span><span class=n>nt</span> <span class=o>*</span> <span class=n>dt</span><span class=p>,</span> <span class=n>xi</span><span class=p>,</span> <span class=n>nu</span><span class=p>)</span> <span class=k>for</span> <span class=n>xi</span> <span class=ow>in</span> <span class=n>x</span><span class=p>])</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>figure</span><span class=p>(</span><span class=n>figsize</span><span class=o>=</span><span class=p>(</span><span class=mi>11</span><span class=p>,</span> <span class=mi>7</span><span class=p>),</span> <span class=n>dpi</span><span class=o>=</span><span class=mi>100</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>plot</span><span class=p>(</span><span class=n>x</span><span class=p>,</span><span class=n>u</span><span class=p>,</span> <span class=n>marker</span><span class=o>=</span><span class=s1>&#39;o&#39;</span><span class=p>,</span> <span class=n>lw</span><span class=o>=</span><span class=mi>2</span><span class=p>,</span> <span class=n>label</span><span class=o>=</span><span class=s1>&#39;Computational&#39;</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>plot</span><span class=p>(</span><span class=n>x</span><span class=p>,</span> <span class=n>u_analytical</span><span class=p>,</span> <span class=n>label</span><span class=o>=</span><span class=s1>&#39;Analytical&#39;</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>xlim</span><span class=p>([</span><span class=mi>0</span><span class=p>,</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>numpy</span><span class=o>.</span><span class=n>pi</span><span class=p>])</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>ylim</span><span class=p>([</span><span class=mi>0</span><span class=p>,</span> <span class=mi>10</span><span class=p>])</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>legend</span><span class=p>();</span>
</span></span></code></pre></div><figure><img loading=lazy src=https://blog.kakaocdn.net/dn/cE7P8B/btq9dF2BVpt/J6GbNhRT4dX1nfB2GPkurK/img.png></figure><h2 id=2-d-linear-convection>2-D Linear Convection<a hidden class=anchor aria-hidden=true href=#2-d-linear-convection>#</a></h2><p>$$ \frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x} + c\frac{\partial u}{\partial y} = 0 $$</p><p>$$ u_{i,j}^{n+1} = u_{i,j}^n-c \frac{\Delta t}{\Delta x}(u_{i,j}^n-u_{i-1,j}^n)-c \frac{\Delta t}{\Delta y}(u_{i,j}^n-u_{i,j-1}^n) $$</p><div class=highlight><pre tabindex=0 class=chroma><code class=language-python data-lang=python><span class=line><span class=cl><span class=kn>from</span> <span class=nn>mpl_toolkits.mplot3d</span> <span class=kn>import</span> <span class=n>Axes3D</span> <span class=c1>##New Library required for projected 3d plots</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=kn>import</span> <span class=nn>numpy</span>
</span></span><span class=line><span class=cl><span class=kn>from</span> <span class=nn>matplotlib</span> <span class=kn>import</span> <span class=n>pyplot</span><span class=p>,</span> <span class=n>cm</span>
</span></span><span class=line><span class=cl><span class=o>%</span><span class=n>matplotlib</span> <span class=n>inline</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=c1>###variable declarations</span>
</span></span><span class=line><span class=cl><span class=n>nx</span> <span class=o>=</span> <span class=mi>81</span>
</span></span><span class=line><span class=cl><span class=n>ny</span> <span class=o>=</span> <span class=mi>81</span>
</span></span><span class=line><span class=cl><span class=n>nt</span> <span class=o>=</span> <span class=mi>100</span>
</span></span><span class=line><span class=cl><span class=n>c</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl><span class=n>dx</span> <span class=o>=</span> <span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=n>nx</span> <span class=o>-</span> <span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>dy</span> <span class=o>=</span> <span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=n>ny</span> <span class=o>-</span> <span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>sigma</span> <span class=o>=</span> <span class=mf>.2</span>
</span></span><span class=line><span class=cl><span class=n>dt</span> <span class=o>=</span> <span class=n>sigma</span> <span class=o>*</span> <span class=n>dx</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>x</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>linspace</span><span class=p>(</span><span class=mi>0</span><span class=p>,</span> <span class=mi>2</span><span class=p>,</span> <span class=n>nx</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>y</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>linspace</span><span class=p>(</span><span class=mi>0</span><span class=p>,</span> <span class=mi>2</span><span class=p>,</span> <span class=n>ny</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>u</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>ones</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span> <span class=c1>##create a 1xn vector of 1&#39;s</span>
</span></span><span class=line><span class=cl><span class=n>un</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>ones</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span> <span class=c1>##</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=c1>###Assign initial conditions</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=c1>##set hat function I.C. : u(.5&lt;=x&lt;=1 &amp;&amp; .5&lt;=y&lt;=1 ) is 2</span>
</span></span><span class=line><span class=cl><span class=n>u</span><span class=p>[</span><span class=nb>int</span><span class=p>(</span><span class=mf>.5</span> <span class=o>/</span> <span class=n>dy</span><span class=p>):</span><span class=nb>int</span><span class=p>(</span><span class=mi>1</span> <span class=o>/</span> <span class=n>dy</span> <span class=o>+</span> <span class=mi>1</span><span class=p>),</span><span class=nb>int</span><span class=p>(</span><span class=mf>.5</span> <span class=o>/</span> <span class=n>dx</span><span class=p>):</span><span class=nb>int</span><span class=p>(</span><span class=mi>1</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>+</span> <span class=mi>1</span><span class=p>)]</span> <span class=o>=</span> <span class=mi>2</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=c1>###Plot Initial Condition</span>
</span></span><span class=line><span class=cl><span class=c1>##the figsize parameter can be used to produce different sized images</span>
</span></span><span class=line><span class=cl><span class=n>fig</span> <span class=o>=</span> <span class=n>pyplot</span><span class=o>.</span><span class=n>figure</span><span class=p>(</span><span class=n>figsize</span><span class=o>=</span><span class=p>(</span><span class=mi>11</span><span class=p>,</span> <span class=mi>7</span><span class=p>),</span> <span class=n>dpi</span><span class=o>=</span><span class=mi>100</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>ax</span> <span class=o>=</span> <span class=n>fig</span><span class=o>.</span><span class=n>gca</span><span class=p>(</span><span class=n>projection</span><span class=o>=</span><span class=s1>&#39;3d&#39;</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>X</span><span class=p>,</span> <span class=n>Y</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>meshgrid</span><span class=p>(</span><span class=n>x</span><span class=p>,</span> <span class=n>y</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>surf</span> <span class=o>=</span> <span class=n>ax</span><span class=o>.</span><span class=n>plot_surface</span><span class=p>(</span><span class=n>X</span><span class=p>,</span> <span class=n>Y</span><span class=p>,</span> <span class=n>u</span><span class=p>[:],</span> <span class=n>cmap</span><span class=o>=</span><span class=n>cm</span><span class=o>.</span><span class=n>viridis</span><span class=p>)</span>
</span></span></code></pre></div><figure><img loading=lazy src=https://blog.kakaocdn.net/dn/MH3sO/btq9fmIfvXs/Tbral2sgJxUQHfJgs1hG61/img.png></figure><h2 id=2-d-convection>2-D Convection<a hidden class=anchor aria-hidden=true href=#2-d-convection>#</a></h2><p>$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = 0 $$$$ \frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} = 0 $$$$ u_{i,j}^{n+1} = u_{i,j}^n - u_{i,j} \frac{\Delta t}{\Delta x} (u_{i,j}^n-u_{i-1,j}^n) - v_{i,j}^n \frac{\Delta t}{\Delta y} (u_{i,j}^n-u_{i,j-1}^n) $$$$ v_{i,j}^{n+1} = v_{i,j}^n - u_{i,j} \frac{\Delta t}{\Delta x} (v_{i,j}^n-v_{i-1,j}^n) - v_{i,j}^n \frac{\Delta t}{\Delta y} (v_{i,j}^n-v_{i,j-1}^n) $$</p><div class=highlight><pre tabindex=0 class=chroma><code class=language-python data-lang=python><span class=line><span class=cl><span class=kn>from</span> <span class=nn>mpl_toolkits.mplot3d</span> <span class=kn>import</span> <span class=n>Axes3D</span>
</span></span><span class=line><span class=cl><span class=kn>from</span> <span class=nn>matplotlib</span> <span class=kn>import</span> <span class=n>pyplot</span><span class=p>,</span> <span class=n>cm</span>
</span></span><span class=line><span class=cl><span class=kn>import</span> <span class=nn>numpy</span>
</span></span><span class=line><span class=cl><span class=o>%</span><span class=n>matplotlib</span> <span class=n>inline</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=c1>###variable declarations</span>
</span></span><span class=line><span class=cl><span class=n>nx</span> <span class=o>=</span> <span class=mi>101</span>
</span></span><span class=line><span class=cl><span class=n>ny</span> <span class=o>=</span> <span class=mi>101</span>
</span></span><span class=line><span class=cl><span class=n>nt</span> <span class=o>=</span> <span class=mi>80</span>
</span></span><span class=line><span class=cl><span class=n>c</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl><span class=n>dx</span> <span class=o>=</span> <span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=n>nx</span> <span class=o>-</span> <span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>dy</span> <span class=o>=</span> <span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=n>ny</span> <span class=o>-</span> <span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>sigma</span> <span class=o>=</span> <span class=mf>.2</span>
</span></span><span class=line><span class=cl><span class=n>dt</span> <span class=o>=</span> <span class=n>sigma</span> <span class=o>*</span> <span class=n>dx</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>x</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>linspace</span><span class=p>(</span><span class=mi>0</span><span class=p>,</span> <span class=mi>2</span><span class=p>,</span> <span class=n>nx</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>y</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>linspace</span><span class=p>(</span><span class=mi>0</span><span class=p>,</span> <span class=mi>2</span><span class=p>,</span> <span class=n>ny</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>u</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>ones</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span> <span class=c1>##create a 1xn vector of 1&#39;s</span>
</span></span><span class=line><span class=cl><span class=n>v</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>ones</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>un</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>ones</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>vn</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>ones</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=c1>###Assign initial conditions</span>
</span></span><span class=line><span class=cl><span class=c1>##set hat function I.C. : u(.5&lt;=x&lt;=1 &amp;&amp; .5&lt;=y&lt;=1 ) is 2</span>
</span></span><span class=line><span class=cl><span class=n>u</span><span class=p>[</span><span class=nb>int</span><span class=p>(</span><span class=mf>.5</span> <span class=o>/</span> <span class=n>dy</span><span class=p>):</span><span class=nb>int</span><span class=p>(</span><span class=mi>1</span> <span class=o>/</span> <span class=n>dy</span> <span class=o>+</span> <span class=mi>1</span><span class=p>),</span> <span class=nb>int</span><span class=p>(</span><span class=mf>.5</span> <span class=o>/</span> <span class=n>dx</span><span class=p>):</span><span class=nb>int</span><span class=p>(</span><span class=mi>1</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>+</span> <span class=mi>1</span><span class=p>)]</span> <span class=o>=</span> <span class=mi>2</span>
</span></span><span class=line><span class=cl><span class=c1>##set hat function I.C. : v(.5&lt;=x&lt;=1 &amp;&amp; .5&lt;=y&lt;=1 ) is 2</span>
</span></span><span class=line><span class=cl><span class=n>v</span><span class=p>[</span><span class=nb>int</span><span class=p>(</span><span class=mf>.5</span> <span class=o>/</span> <span class=n>dy</span><span class=p>):</span><span class=nb>int</span><span class=p>(</span><span class=mi>1</span> <span class=o>/</span> <span class=n>dy</span> <span class=o>+</span> <span class=mi>1</span><span class=p>),</span> <span class=nb>int</span><span class=p>(</span><span class=mf>.5</span> <span class=o>/</span> <span class=n>dx</span><span class=p>):</span><span class=nb>int</span><span class=p>(</span><span class=mi>1</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>+</span> <span class=mi>1</span><span class=p>)]</span> <span class=o>=</span> <span class=mi>2</span>
</span></span><span class=line><span class=cl><span class=n>fig</span> <span class=o>=</span> <span class=n>pyplot</span><span class=o>.</span><span class=n>figure</span><span class=p>(</span><span class=n>figsize</span><span class=o>=</span><span class=p>(</span><span class=mi>11</span><span class=p>,</span> <span class=mi>7</span><span class=p>),</span> <span class=n>dpi</span><span class=o>=</span><span class=mi>100</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>ax</span> <span class=o>=</span> <span class=n>fig</span><span class=o>.</span><span class=n>gca</span><span class=p>(</span><span class=n>projection</span><span class=o>=</span><span class=s1>&#39;3d&#39;</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>X</span><span class=p>,</span> <span class=n>Y</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>meshgrid</span><span class=p>(</span><span class=n>x</span><span class=p>,</span> <span class=n>y</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>ax</span><span class=o>.</span><span class=n>plot_surface</span><span class=p>(</span><span class=n>X</span><span class=p>,</span> <span class=n>Y</span><span class=p>,</span> <span class=n>u</span><span class=p>,</span> <span class=n>cmap</span><span class=o>=</span><span class=n>cm</span><span class=o>.</span><span class=n>viridis</span><span class=p>,</span> <span class=n>rstride</span><span class=o>=</span><span class=mi>2</span><span class=p>,</span> <span class=n>cstride</span><span class=o>=</span><span class=mi>2</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>ax</span><span class=o>.</span><span class=n>set_xlabel</span><span class=p>(</span><span class=s1>&#39;$x$&#39;</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>ax</span><span class=o>.</span><span class=n>set_ylabel</span><span class=p>(</span><span class=s1>&#39;$y$&#39;</span><span class=p>);</span>
</span></span></code></pre></div><figure><img loading=lazy src=https://blog.kakaocdn.net/dn/bbiKmO/btq9hyAHR99/3KeBvxXPvCzXXYqrTlSj9k/img.png></figure><h2 id=2d-diffusion>2D Diffusion<a hidden class=anchor aria-hidden=true href=#2d-diffusion>#</a></h2><p>$$ \frac{\partial u}{\partial t} = \nu \frac{\partial ^2 u}{\partial x^2} + \nu \frac{\partial ^2 u}{\partial y^2} $$$$ \begin{split}u_{i,j}^{n+1} = u_{i,j}^n &+ \frac{\nu \Delta t}{\Delta x^2}(u_{i+1,j}^n - 2 u_{i,j}^n + u_{i-1,j}^n) \&+ \frac{\nu \Delta t}{\Delta y^2}(u_{i,j+1}^n-2 u_{i,j}^n + u_{i,j-1}^n)\end{split} $$</p><div class=highlight><pre tabindex=0 class=chroma><code class=language-python data-lang=python><span class=line><span class=cl><span class=kn>import</span> <span class=nn>numpy</span>
</span></span><span class=line><span class=cl><span class=kn>from</span> <span class=nn>matplotlib</span> <span class=kn>import</span> <span class=n>pyplot</span><span class=p>,</span> <span class=n>cm</span>
</span></span><span class=line><span class=cl><span class=kn>from</span> <span class=nn>mpl_toolkits.mplot3d</span> <span class=kn>import</span> <span class=n>Axes3D</span> <span class=c1>##library for 3d projection plots</span>
</span></span><span class=line><span class=cl><span class=o>%</span><span class=n>matplotlib</span> <span class=n>inline</span>
</span></span><span class=line><span class=cl><span class=c1>###variable declarations</span>
</span></span><span class=line><span class=cl><span class=n>nx</span> <span class=o>=</span> <span class=mi>31</span>
</span></span><span class=line><span class=cl><span class=n>ny</span> <span class=o>=</span> <span class=mi>31</span>
</span></span><span class=line><span class=cl><span class=n>nt</span> <span class=o>=</span> <span class=mi>17</span>
</span></span><span class=line><span class=cl><span class=n>nu</span> <span class=o>=</span> <span class=mf>.05</span>
</span></span><span class=line><span class=cl><span class=n>dx</span> <span class=o>=</span> <span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=n>nx</span> <span class=o>-</span> <span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>dy</span> <span class=o>=</span> <span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=n>ny</span> <span class=o>-</span> <span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>sigma</span> <span class=o>=</span> <span class=mf>.25</span>
</span></span><span class=line><span class=cl><span class=n>dt</span> <span class=o>=</span> <span class=n>sigma</span> <span class=o>*</span> <span class=n>dx</span> <span class=o>*</span> <span class=n>dy</span> <span class=o>/</span> <span class=n>nu</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>x</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>linspace</span><span class=p>(</span><span class=mi>0</span><span class=p>,</span> <span class=mi>2</span><span class=p>,</span> <span class=n>nx</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>y</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>linspace</span><span class=p>(</span><span class=mi>0</span><span class=p>,</span> <span class=mi>2</span><span class=p>,</span> <span class=n>ny</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>u</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>ones</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span> <span class=c1># create a 1xn vector of 1&#39;s</span>
</span></span><span class=line><span class=cl><span class=n>un</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>ones</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=c1>###Assign initial conditions</span>
</span></span><span class=line><span class=cl><span class=c1># set hat function I.C. : u(.5&lt;=x&lt;=1 &amp;&amp; .5&lt;=y&lt;=1 ) is 2</span>
</span></span><span class=line><span class=cl><span class=n>u</span><span class=p>[</span><span class=nb>int</span><span class=p>(</span><span class=mf>.5</span> <span class=o>/</span> <span class=n>dy</span><span class=p>):</span><span class=nb>int</span><span class=p>(</span><span class=mi>1</span> <span class=o>/</span> <span class=n>dy</span> <span class=o>+</span> <span class=mi>1</span><span class=p>),</span><span class=nb>int</span><span class=p>(</span><span class=mf>.5</span> <span class=o>/</span> <span class=n>dx</span><span class=p>):</span><span class=nb>int</span><span class=p>(</span><span class=mi>1</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>+</span> <span class=mi>1</span><span class=p>)]</span> <span class=o>=</span> <span class=mi>2</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=c1>###Run through nt timesteps</span>
</span></span><span class=line><span class=cl><span class=k>def</span> <span class=nf>diffuse</span><span class=p>(</span><span class=n>nt</span><span class=p>):</span>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[</span><span class=nb>int</span><span class=p>(</span><span class=mf>.5</span> <span class=o>/</span> <span class=n>dy</span><span class=p>):</span><span class=nb>int</span><span class=p>(</span><span class=mi>1</span> <span class=o>/</span> <span class=n>dy</span> <span class=o>+</span> <span class=mi>1</span><span class=p>),</span><span class=nb>int</span><span class=p>(</span><span class=mf>.5</span> <span class=o>/</span> <span class=n>dx</span><span class=p>):</span><span class=nb>int</span><span class=p>(</span><span class=mi>1</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>+</span> <span class=mi>1</span><span class=p>)]</span> <span class=o>=</span> <span class=mi>2</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=k>for</span> <span class=n>n</span> <span class=ow>in</span> <span class=nb>range</span><span class=p>(</span><span class=n>nt</span> <span class=o>+</span> <span class=mi>1</span><span class=p>):</span>
</span></span><span class=line><span class=cl> <span class=n>un</span> <span class=o>=</span> <span class=n>u</span><span class=o>.</span><span class=n>copy</span><span class=p>()</span>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>nu</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>2</span><span class=p>:]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>nu</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span><span class=mi>1</span><span class=p>:</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>un</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]))</span>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[</span><span class=mi>0</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[:,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[:,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=n>fig</span> <span class=o>=</span> <span class=n>pyplot</span><span class=o>.</span><span class=n>figure</span><span class=p>()</span>
</span></span><span class=line><span class=cl> <span class=n>ax</span> <span class=o>=</span> <span class=n>fig</span><span class=o>.</span><span class=n>gca</span><span class=p>(</span><span class=n>projection</span><span class=o>=</span><span class=s1>&#39;3d&#39;</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>surf</span> <span class=o>=</span> <span class=n>ax</span><span class=o>.</span><span class=n>plot_surface</span><span class=p>(</span><span class=n>X</span><span class=p>,</span> <span class=n>Y</span><span class=p>,</span> <span class=n>u</span><span class=p>[:],</span> <span class=n>rstride</span><span class=o>=</span><span class=mi>1</span><span class=p>,</span> <span class=n>cstride</span><span class=o>=</span><span class=mi>1</span><span class=p>,</span> <span class=n>cmap</span><span class=o>=</span><span class=n>cm</span><span class=o>.</span><span class=n>viridis</span><span class=p>,</span>
</span></span><span class=line><span class=cl> <span class=n>linewidth</span><span class=o>=</span><span class=mi>0</span><span class=p>,</span> <span class=n>antialiased</span><span class=o>=</span><span class=kc>True</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>ax</span><span class=o>.</span><span class=n>set_zlim</span><span class=p>(</span><span class=mi>1</span><span class=p>,</span> <span class=mf>2.5</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>ax</span><span class=o>.</span><span class=n>set_xlabel</span><span class=p>(</span><span class=s1>&#39;$x$&#39;</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>ax</span><span class=o>.</span><span class=n>set_ylabel</span><span class=p>(</span><span class=s1>&#39;$y$&#39;</span><span class=p>);</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>diffuse</span><span class=p>(</span><span class=mi>14</span><span class=p>)</span>
</span></span></code></pre></div><figure><img loading=lazy src=https://blog.kakaocdn.net/dn/eLwQEW/btq9e0ysgnx/YVHruuNPlMa6pJODaIGJdK/img.png></figure><h2 id=burgers-equation-in-2d>Burgers&rsquo; Equation in 2D<a hidden class=anchor aria-hidden=true href=#burgers-equation-in-2d>#</a></h2><p>$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = \nu ; \left(\frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{\partial y^2}\right) $$$$ \frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} = \nu ; \left(\frac{\partial ^2 v}{\partial x^2} + \frac{\partial ^2 v}{\partial y^2}\right) $$$$ \begin{split}v_{i,j}^{n+1} = & v_{i,j}^n - \frac{\Delta t}{\Delta x} u_{i,j}^n (v_{i,j}^n - v_{i-1,j}^n) - \frac{\Delta t}{\Delta y} v_{i,j}^n (v_{i,j}^n - v_{i,j-1}^n) \&+ \frac{\nu \Delta t}{\Delta x^2}(v_{i+1,j}^n-2v_{i,j}^n+v_{i-1,j}^n) + \frac{\nu \Delta t}{\Delta y^2} (v_{i,j+1}^n - 2v_{i,j}^n + v_{i,j-1}^n)\end{split} $$</p><div class=highlight><pre tabindex=0 class=chroma><code class=language-python data-lang=python><span class=line><span class=cl><span class=kn>import</span> <span class=nn>numpy</span>
</span></span><span class=line><span class=cl><span class=kn>from</span> <span class=nn>matplotlib</span> <span class=kn>import</span> <span class=n>pyplot</span><span class=p>,</span> <span class=n>cm</span>
</span></span><span class=line><span class=cl><span class=kn>from</span> <span class=nn>mpl_toolkits.mplot3d</span> <span class=kn>import</span> <span class=n>Axes3D</span>
</span></span><span class=line><span class=cl><span class=o>%</span><span class=n>matplotlib</span> <span class=n>inline</span>
</span></span><span class=line><span class=cl><span class=c1>###variable declarations</span>
</span></span><span class=line><span class=cl><span class=n>nx</span> <span class=o>=</span> <span class=mi>41</span>
</span></span><span class=line><span class=cl><span class=n>ny</span> <span class=o>=</span> <span class=mi>41</span>
</span></span><span class=line><span class=cl><span class=n>nt</span> <span class=o>=</span> <span class=mi>120</span>
</span></span><span class=line><span class=cl><span class=n>c</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl><span class=n>dx</span> <span class=o>=</span> <span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=n>nx</span> <span class=o>-</span> <span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>dy</span> <span class=o>=</span> <span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=n>ny</span> <span class=o>-</span> <span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>sigma</span> <span class=o>=</span> <span class=mf>.0009</span>
</span></span><span class=line><span class=cl><span class=n>nu</span> <span class=o>=</span> <span class=mf>0.01</span>
</span></span><span class=line><span class=cl><span class=n>dt</span> <span class=o>=</span> <span class=n>sigma</span> <span class=o>*</span> <span class=n>dx</span> <span class=o>*</span> <span class=n>dy</span> <span class=o>/</span> <span class=n>nu</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>x</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>linspace</span><span class=p>(</span><span class=mi>0</span><span class=p>,</span> <span class=mi>2</span><span class=p>,</span> <span class=n>nx</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>y</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>linspace</span><span class=p>(</span><span class=mi>0</span><span class=p>,</span> <span class=mi>2</span><span class=p>,</span> <span class=n>ny</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>u</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>ones</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span> <span class=c1># create a 1xn vector of 1&#39;s</span>
</span></span><span class=line><span class=cl><span class=n>v</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>ones</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>un</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>ones</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>vn</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>ones</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>comb</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>ones</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=c1>###Assign initial conditions</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=c1>##set hat function I.C. : u(.5&lt;=x&lt;=1 &amp;&amp; .5&lt;=y&lt;=1 ) is 2</span>
</span></span><span class=line><span class=cl><span class=n>u</span><span class=p>[</span><span class=nb>int</span><span class=p>(</span><span class=mf>.5</span> <span class=o>/</span> <span class=n>dy</span><span class=p>):</span><span class=nb>int</span><span class=p>(</span><span class=mi>1</span> <span class=o>/</span> <span class=n>dy</span> <span class=o>+</span> <span class=mi>1</span><span class=p>),</span><span class=nb>int</span><span class=p>(</span><span class=mf>.5</span> <span class=o>/</span> <span class=n>dx</span><span class=p>):</span><span class=nb>int</span><span class=p>(</span><span class=mi>1</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>+</span> <span class=mi>1</span><span class=p>)]</span> <span class=o>=</span> <span class=mi>2</span>
</span></span><span class=line><span class=cl><span class=c1>##set hat function I.C. : u(.5&lt;=x&lt;=1 &amp;&amp; .5&lt;=y&lt;=1 ) is 2</span>
</span></span><span class=line><span class=cl><span class=n>v</span><span class=p>[</span><span class=nb>int</span><span class=p>(</span><span class=mf>.5</span> <span class=o>/</span> <span class=n>dy</span><span class=p>):</span><span class=nb>int</span><span class=p>(</span><span class=mi>1</span> <span class=o>/</span> <span class=n>dy</span> <span class=o>+</span> <span class=mi>1</span><span class=p>),</span><span class=nb>int</span><span class=p>(</span><span class=mf>.5</span> <span class=o>/</span> <span class=n>dx</span><span class=p>):</span><span class=nb>int</span><span class=p>(</span><span class=mi>1</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>+</span> <span class=mi>1</span><span class=p>)]</span> <span class=o>=</span> <span class=mi>2</span>
</span></span><span class=line><span class=cl><span class=c1>###(plot ICs)</span>
</span></span><span class=line><span class=cl><span class=k>for</span> <span class=n>n</span> <span class=ow>in</span> <span class=nb>range</span><span class=p>(</span><span class=n>nt</span> <span class=o>+</span> <span class=mi>1</span><span class=p>):</span> <span class=c1>##loop across number of time steps</span>
</span></span><span class=line><span class=cl> <span class=n>un</span> <span class=o>=</span> <span class=n>u</span><span class=o>.</span><span class=n>copy</span><span class=p>()</span>
</span></span><span class=line><span class=cl> <span class=n>vn</span> <span class=o>=</span> <span class=n>v</span><span class=o>.</span><span class=n>copy</span><span class=p>()</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>*</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>dt</span> <span class=o>/</span> <span class=n>dy</span> <span class=o>*</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>un</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>nu</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span><span class=mi>2</span><span class=p>:]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>nu</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>un</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=n>v</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>*</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>dt</span> <span class=o>/</span> <span class=n>dy</span> <span class=o>*</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>vn</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>nu</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>2</span><span class=p>:]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>nu</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>vn</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[</span><span class=mi>0</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[:,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[:,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=n>v</span><span class=p>[</span><span class=mi>0</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl> <span class=n>v</span><span class=p>[</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl> <span class=n>v</span><span class=p>[:,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl> <span class=n>v</span><span class=p>[:,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl><span class=n>fig</span> <span class=o>=</span> <span class=n>pyplot</span><span class=o>.</span><span class=n>figure</span><span class=p>(</span><span class=n>figsize</span><span class=o>=</span><span class=p>(</span><span class=mi>11</span><span class=p>,</span> <span class=mi>7</span><span class=p>),</span> <span class=n>dpi</span><span class=o>=</span><span class=mi>100</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>ax</span> <span class=o>=</span> <span class=n>fig</span><span class=o>.</span><span class=n>gca</span><span class=p>(</span><span class=n>projection</span><span class=o>=</span><span class=s1>&#39;3d&#39;</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>X</span><span class=p>,</span> <span class=n>Y</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>meshgrid</span><span class=p>(</span><span class=n>x</span><span class=p>,</span> <span class=n>y</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>ax</span><span class=o>.</span><span class=n>plot_surface</span><span class=p>(</span><span class=n>X</span><span class=p>,</span> <span class=n>Y</span><span class=p>,</span> <span class=n>u</span><span class=p>,</span> <span class=n>cmap</span><span class=o>=</span><span class=n>cm</span><span class=o>.</span><span class=n>viridis</span><span class=p>,</span> <span class=n>rstride</span><span class=o>=</span><span class=mi>1</span><span class=p>,</span> <span class=n>cstride</span><span class=o>=</span><span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>ax</span><span class=o>.</span><span class=n>plot_surface</span><span class=p>(</span><span class=n>X</span><span class=p>,</span> <span class=n>Y</span><span class=p>,</span> <span class=n>v</span><span class=p>,</span> <span class=n>cmap</span><span class=o>=</span><span class=n>cm</span><span class=o>.</span><span class=n>viridis</span><span class=p>,</span> <span class=n>rstride</span><span class=o>=</span><span class=mi>1</span><span class=p>,</span> <span class=n>cstride</span><span class=o>=</span><span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>ax</span><span class=o>.</span><span class=n>set_xlabel</span><span class=p>(</span><span class=s1>&#39;$x$&#39;</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>ax</span><span class=o>.</span><span class=n>set_ylabel</span><span class=p>(</span><span class=s1>&#39;$y$&#39;</span><span class=p>);</span>
</span></span></code></pre></div><figure><img loading=lazy src=https://blog.kakaocdn.net/dn/PL2CD/btq9fcyC1VV/MC1B8I2YedaaCFr5Lr06KK/img.png></figure><h2 id=2d-laplace-equation>2D Laplace Equation<a hidden class=anchor aria-hidden=true href=#2d-laplace-equation>#</a></h2><p>$$ \frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = 0 $$$$ p_{i,j}^n = \frac{\Delta y^2(p_{i+1,j}^n+p_{i-1,j}^n)+\Delta x^2(p_{i,j+1}^n + p_{i,j-1}^n)}{2(\Delta x^2 + \Delta y^2)} $$</p><div class=highlight><pre tabindex=0 class=chroma><code class=language-python data-lang=python><span class=line><span class=cl><span class=kn>import</span> <span class=nn>numpy</span>
</span></span><span class=line><span class=cl><span class=kn>from</span> <span class=nn>matplotlib</span> <span class=kn>import</span> <span class=n>pyplot</span><span class=p>,</span> <span class=n>cm</span>
</span></span><span class=line><span class=cl><span class=kn>from</span> <span class=nn>mpl_toolkits.mplot3d</span> <span class=kn>import</span> <span class=n>Axes3D</span>
</span></span><span class=line><span class=cl><span class=o>%</span><span class=n>matplotlib</span> <span class=n>inline</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=k>def</span> <span class=nf>plot2D</span><span class=p>(</span><span class=n>x</span><span class=p>,</span> <span class=n>y</span><span class=p>,</span> <span class=n>p</span><span class=p>):</span>
</span></span><span class=line><span class=cl> <span class=n>fig</span> <span class=o>=</span> <span class=n>pyplot</span><span class=o>.</span><span class=n>figure</span><span class=p>(</span><span class=n>figsize</span><span class=o>=</span><span class=p>(</span><span class=mi>11</span><span class=p>,</span> <span class=mi>7</span><span class=p>),</span> <span class=n>dpi</span><span class=o>=</span><span class=mi>100</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>ax</span> <span class=o>=</span> <span class=n>fig</span><span class=o>.</span><span class=n>gca</span><span class=p>(</span><span class=n>projection</span><span class=o>=</span><span class=s1>&#39;3d&#39;</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>X</span><span class=p>,</span> <span class=n>Y</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>meshgrid</span><span class=p>(</span><span class=n>x</span><span class=p>,</span> <span class=n>y</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>surf</span> <span class=o>=</span> <span class=n>ax</span><span class=o>.</span><span class=n>plot_surface</span><span class=p>(</span><span class=n>X</span><span class=p>,</span> <span class=n>Y</span><span class=p>,</span> <span class=n>p</span><span class=p>[:],</span> <span class=n>rstride</span><span class=o>=</span><span class=mi>1</span><span class=p>,</span> <span class=n>cstride</span><span class=o>=</span><span class=mi>1</span><span class=p>,</span> <span class=n>cmap</span><span class=o>=</span><span class=n>cm</span><span class=o>.</span><span class=n>viridis</span><span class=p>,</span>
</span></span><span class=line><span class=cl> <span class=n>linewidth</span><span class=o>=</span><span class=mi>0</span><span class=p>,</span> <span class=n>antialiased</span><span class=o>=</span><span class=kc>False</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>ax</span><span class=o>.</span><span class=n>set_xlim</span><span class=p>(</span><span class=mi>0</span><span class=p>,</span> <span class=mi>2</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>ax</span><span class=o>.</span><span class=n>set_ylim</span><span class=p>(</span><span class=mi>0</span><span class=p>,</span> <span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>ax</span><span class=o>.</span><span class=n>view_init</span><span class=p>(</span><span class=mi>30</span><span class=p>,</span> <span class=mi>225</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>ax</span><span class=o>.</span><span class=n>set_xlabel</span><span class=p>(</span><span class=s1>&#39;$x$&#39;</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>ax</span><span class=o>.</span><span class=n>set_ylabel</span><span class=p>(</span><span class=s1>&#39;$y$&#39;</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=k>def</span> <span class=nf>laplace2d</span><span class=p>(</span><span class=n>p</span><span class=p>,</span> <span class=n>y</span><span class=p>,</span> <span class=n>dx</span><span class=p>,</span> <span class=n>dy</span><span class=p>,</span> <span class=n>l1norm_target</span><span class=p>):</span>
</span></span><span class=line><span class=cl> <span class=n>l1norm</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl> <span class=n>pn</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>empty_like</span><span class=p>(</span><span class=n>p</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=k>while</span> <span class=n>l1norm</span> <span class=o>&gt;</span> <span class=n>l1norm_target</span><span class=p>:</span>
</span></span><span class=line><span class=cl> <span class=n>pn</span> <span class=o>=</span> <span class=n>p</span><span class=o>.</span><span class=n>copy</span><span class=p>()</span>
</span></span><span class=line><span class=cl> <span class=n>p</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=p>((</span><span class=n>dy</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span> <span class=p>(</span><span class=n>pn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>2</span><span class=p>:]</span> <span class=o>+</span> <span class=n>pn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span> <span class=p>(</span><span class=n>pn</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>pn</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]))</span> <span class=o>/</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=p>(</span><span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>+</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span><span class=p>)))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=n>p</span><span class=p>[:,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>=</span> <span class=mi>0</span> <span class=c1># p = 0 @ x = 0</span>
</span></span><span class=line><span class=cl> <span class=n>p</span><span class=p>[:,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=n>y</span> <span class=c1># p = y @ x = 2</span>
</span></span><span class=line><span class=cl> <span class=n>p</span><span class=p>[</span><span class=mi>0</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=n>p</span><span class=p>[</span><span class=mi>1</span><span class=p>,</span> <span class=p>:]</span> <span class=c1># dp/dy = 0 @ y = 0</span>
</span></span><span class=line><span class=cl> <span class=n>p</span><span class=p>[</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=n>p</span><span class=p>[</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=p>:]</span> <span class=c1># dp/dy = 0 @ y = 1</span>
</span></span><span class=line><span class=cl> <span class=n>l1norm</span> <span class=o>=</span> <span class=p>(</span><span class=n>numpy</span><span class=o>.</span><span class=n>sum</span><span class=p>(</span><span class=n>numpy</span><span class=o>.</span><span class=n>abs</span><span class=p>(</span><span class=n>p</span><span class=p>[:])</span> <span class=o>-</span> <span class=n>numpy</span><span class=o>.</span><span class=n>abs</span><span class=p>(</span><span class=n>pn</span><span class=p>[:]))</span> <span class=o>/</span>
</span></span><span class=line><span class=cl> <span class=n>numpy</span><span class=o>.</span><span class=n>sum</span><span class=p>(</span><span class=n>numpy</span><span class=o>.</span><span class=n>abs</span><span class=p>(</span><span class=n>pn</span><span class=p>[:])))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=k>return</span> <span class=n>p</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>nx</span> <span class=o>=</span> <span class=mi>31</span>
</span></span><span class=line><span class=cl><span class=n>ny</span> <span class=o>=</span> <span class=mi>31</span>
</span></span><span class=line><span class=cl><span class=n>c</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl><span class=n>dx</span> <span class=o>=</span> <span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=n>nx</span> <span class=o>-</span> <span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>dy</span> <span class=o>=</span> <span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=n>ny</span> <span class=o>-</span> <span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=c1>##initial conditions</span>
</span></span><span class=line><span class=cl><span class=n>p</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span> <span class=c1># create a XxY vector of 0&#39;s</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=c1>##plotting aids</span>
</span></span><span class=line><span class=cl><span class=n>x</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>linspace</span><span class=p>(</span><span class=mi>0</span><span class=p>,</span> <span class=mi>2</span><span class=p>,</span> <span class=n>nx</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>y</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>linspace</span><span class=p>(</span><span class=mi>0</span><span class=p>,</span> <span class=mi>1</span><span class=p>,</span> <span class=n>ny</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=c1>##boundary conditions</span>
</span></span><span class=line><span class=cl><span class=n>p</span><span class=p>[:,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>=</span> <span class=mi>0</span> <span class=c1># p = 0 @ x = 0</span>
</span></span><span class=line><span class=cl><span class=n>p</span><span class=p>[:,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=n>y</span> <span class=c1># p = y @ x = 2</span>
</span></span><span class=line><span class=cl><span class=n>p</span><span class=p>[</span><span class=mi>0</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=n>p</span><span class=p>[</span><span class=mi>1</span><span class=p>,</span> <span class=p>:]</span> <span class=c1># dp/dy = 0 @ y = 0</span>
</span></span><span class=line><span class=cl><span class=n>p</span><span class=p>[</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=n>p</span><span class=p>[</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=p>:]</span> <span class=c1># dp/dy = 0 @ y = 1</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>p</span> <span class=o>=</span> <span class=n>laplace2d</span><span class=p>(</span><span class=n>p</span><span class=p>,</span> <span class=n>y</span><span class=p>,</span> <span class=n>dx</span><span class=p>,</span> <span class=n>dy</span><span class=p>,</span> <span class=mf>1e-4</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>plot2D</span><span class=p>(</span><span class=n>x</span><span class=p>,</span> <span class=n>y</span><span class=p>,</span> <span class=n>p</span><span class=p>)</span>
</span></span></code></pre></div><figure><img loading=lazy src=https://blog.kakaocdn.net/dn/bxvdGX/btq9goyjEC1/YkjgRDKkIZuAe2isKNsv60/img.png></figure><h2 id=2d-poisson-equation>2D Poisson Equation<a hidden class=anchor aria-hidden=true href=#2d-poisson-equation>#</a></h2><p>$$ \frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = b $$$$ p_{i,j}^{n}=\frac{(p_{i+1,j}^{n}+p_{i-1,j}^{n})\Delta y^2+(p_{i,j+1}^{n}+p_{i,j-1}^{n})\Delta x^2-b_{i,j}^{n}\Delta x^2\Delta y^2}{2(\Delta x^2+\Delta y^2)} $$</p><div class=highlight><pre tabindex=0 class=chroma><code class=language-python data-lang=python><span class=line><span class=cl><span class=kn>import</span> <span class=nn>numpy</span>
</span></span><span class=line><span class=cl><span class=kn>from</span> <span class=nn>matplotlib</span> <span class=kn>import</span> <span class=n>pyplot</span><span class=p>,</span> <span class=n>cm</span>
</span></span><span class=line><span class=cl><span class=kn>from</span> <span class=nn>mpl_toolkits.mplot3d</span> <span class=kn>import</span> <span class=n>Axes3D</span>
</span></span><span class=line><span class=cl><span class=o>%</span><span class=n>matplotlib</span> <span class=n>inline</span>
</span></span><span class=line><span class=cl><span class=c1># Parameters</span>
</span></span><span class=line><span class=cl><span class=n>nx</span> <span class=o>=</span> <span class=mi>50</span>
</span></span><span class=line><span class=cl><span class=n>ny</span> <span class=o>=</span> <span class=mi>50</span>
</span></span><span class=line><span class=cl><span class=n>nt</span> <span class=o>=</span> <span class=mi>100</span>
</span></span><span class=line><span class=cl><span class=n>xmin</span> <span class=o>=</span> <span class=mi>0</span>
</span></span><span class=line><span class=cl><span class=n>xmax</span> <span class=o>=</span> <span class=mi>2</span>
</span></span><span class=line><span class=cl><span class=n>ymin</span> <span class=o>=</span> <span class=mi>0</span>
</span></span><span class=line><span class=cl><span class=n>ymax</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>dx</span> <span class=o>=</span> <span class=p>(</span><span class=n>xmax</span> <span class=o>-</span> <span class=n>xmin</span><span class=p>)</span> <span class=o>/</span> <span class=p>(</span><span class=n>nx</span> <span class=o>-</span> <span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>dy</span> <span class=o>=</span> <span class=p>(</span><span class=n>ymax</span> <span class=o>-</span> <span class=n>ymin</span><span class=p>)</span> <span class=o>/</span> <span class=p>(</span><span class=n>ny</span> <span class=o>-</span> <span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=c1># Initialization</span>
</span></span><span class=line><span class=cl><span class=n>p</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>pd</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>b</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>x</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>linspace</span><span class=p>(</span><span class=n>xmin</span><span class=p>,</span> <span class=n>xmax</span><span class=p>,</span> <span class=n>nx</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>y</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>linspace</span><span class=p>(</span><span class=n>xmin</span><span class=p>,</span> <span class=n>xmax</span><span class=p>,</span> <span class=n>ny</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=c1># Source</span>
</span></span><span class=line><span class=cl><span class=n>b</span><span class=p>[</span><span class=nb>int</span><span class=p>(</span><span class=n>ny</span> <span class=o>/</span> <span class=mi>4</span><span class=p>),</span> <span class=nb>int</span><span class=p>(</span><span class=n>nx</span> <span class=o>/</span> <span class=mi>4</span><span class=p>)]</span> <span class=o>=</span> <span class=mi>100</span>
</span></span><span class=line><span class=cl><span class=n>b</span><span class=p>[</span><span class=nb>int</span><span class=p>(</span><span class=mi>3</span> <span class=o>*</span> <span class=n>ny</span> <span class=o>/</span> <span class=mi>4</span><span class=p>),</span> <span class=nb>int</span><span class=p>(</span><span class=mi>3</span> <span class=o>*</span> <span class=n>nx</span> <span class=o>/</span> <span class=mi>4</span><span class=p>)]</span> <span class=o>=</span> <span class=o>-</span><span class=mi>100</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=k>for</span> <span class=n>it</span> <span class=ow>in</span> <span class=nb>range</span><span class=p>(</span><span class=n>nt</span><span class=p>):</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=n>pd</span> <span class=o>=</span> <span class=n>p</span><span class=o>.</span><span class=n>copy</span><span class=p>()</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=n>p</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=p>(((</span><span class=n>pd</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>2</span><span class=p>:]</span> <span class=o>+</span> <span class=n>pd</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>*</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>pd</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>pd</span><span class=p>[:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>*</span> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>b</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>*</span> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span><span class=p>)</span> <span class=o>/</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=p>(</span><span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>+</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span><span class=p>)))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=n>p</span><span class=p>[</span><span class=mi>0</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=mi>0</span>
</span></span><span class=line><span class=cl> <span class=n>p</span><span class=p>[</span><span class=n>ny</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=mi>0</span>
</span></span><span class=line><span class=cl> <span class=n>p</span><span class=p>[:,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>=</span> <span class=mi>0</span>
</span></span><span class=line><span class=cl> <span class=n>p</span><span class=p>[:,</span> <span class=n>nx</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=mi>0</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=k>def</span> <span class=nf>plot2D</span><span class=p>(</span><span class=n>x</span><span class=p>,</span> <span class=n>y</span><span class=p>,</span> <span class=n>p</span><span class=p>):</span>
</span></span><span class=line><span class=cl> <span class=n>fig</span> <span class=o>=</span> <span class=n>pyplot</span><span class=o>.</span><span class=n>figure</span><span class=p>(</span><span class=n>figsize</span><span class=o>=</span><span class=p>(</span><span class=mi>11</span><span class=p>,</span> <span class=mi>7</span><span class=p>),</span> <span class=n>dpi</span><span class=o>=</span><span class=mi>100</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>ax</span> <span class=o>=</span> <span class=n>fig</span><span class=o>.</span><span class=n>gca</span><span class=p>(</span><span class=n>projection</span><span class=o>=</span><span class=s1>&#39;3d&#39;</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>X</span><span class=p>,</span> <span class=n>Y</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>meshgrid</span><span class=p>(</span><span class=n>x</span><span class=p>,</span> <span class=n>y</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>surf</span> <span class=o>=</span> <span class=n>ax</span><span class=o>.</span><span class=n>plot_surface</span><span class=p>(</span><span class=n>X</span><span class=p>,</span> <span class=n>Y</span><span class=p>,</span> <span class=n>p</span><span class=p>[:],</span> <span class=n>rstride</span><span class=o>=</span><span class=mi>1</span><span class=p>,</span> <span class=n>cstride</span><span class=o>=</span><span class=mi>1</span><span class=p>,</span> <span class=n>cmap</span><span class=o>=</span><span class=n>cm</span><span class=o>.</span><span class=n>viridis</span><span class=p>,</span>
</span></span><span class=line><span class=cl> <span class=n>linewidth</span><span class=o>=</span><span class=mi>0</span><span class=p>,</span> <span class=n>antialiased</span><span class=o>=</span><span class=kc>False</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>ax</span><span class=o>.</span><span class=n>view_init</span><span class=p>(</span><span class=mi>30</span><span class=p>,</span> <span class=mi>225</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>ax</span><span class=o>.</span><span class=n>set_xlabel</span><span class=p>(</span><span class=s1>&#39;$x$&#39;</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>ax</span><span class=o>.</span><span class=n>set_ylabel</span><span class=p>(</span><span class=s1>&#39;$y$&#39;</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>plot2D</span><span class=p>(</span><span class=n>x</span><span class=p>,</span> <span class=n>y</span><span class=p>,</span> <span class=n>p</span><span class=p>)</span>
</span></span></code></pre></div><h2 id=heading><a hidden class=anchor aria-hidden=true href=#heading>#</a></h2><h2 id=cavity-flow-with-navierstokes>Cavity Flow with Navier–Stokes<a hidden class=anchor aria-hidden=true href=#cavity-flow-with-navierstokes>#</a></h2><p>$$ \frac{\partial \vec{v}}{\partial t}+(\vec{v}\cdot\nabla)\vec{v}=-\frac{1}{\rho}\nabla p + \nu \nabla^2\vec{v} $$$$ \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial x}+\nu \left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2} \right) $$$$ \frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2} = -\rho\left(\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+2\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+\frac{\partial v}{\partial y}\frac{\partial v}{\partial y} \right) $$$$ \begin{split}p_{i,j}^{n} = & \frac{\left(p_{i+1,j}^{n}+p_{i-1,j}^{n}\right) \Delta y^2 + \left(p_{i,j+1}^{n}+p_{i,j-1}^{n}\right) \Delta x^2}{2\left(\Delta x^2+\Delta y^2\right)} \& -\frac{\rho\Delta x^2\Delta y^2}{2\left(\Delta x^2+\Delta y^2\right)} \& \times \left[\frac{1}{\Delta t}\left(\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}+\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right)-\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x} -2\frac{u_{i,j+1}-u_{i,j-1}}{2\Delta y}\frac{v_{i+1,j}-v_{i-1,j}}{2\Delta x}-\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right]\end{split} $$</p><div class=highlight><pre tabindex=0 class=chroma><code class=language-python data-lang=python><span class=line><span class=cl><span class=kn>import</span> <span class=nn>numpy</span>
</span></span><span class=line><span class=cl><span class=kn>from</span> <span class=nn>matplotlib</span> <span class=kn>import</span> <span class=n>pyplot</span><span class=p>,</span> <span class=n>cm</span>
</span></span><span class=line><span class=cl><span class=kn>from</span> <span class=nn>mpl_toolkits.mplot3d</span> <span class=kn>import</span> <span class=n>Axes3D</span>
</span></span><span class=line><span class=cl><span class=o>%</span><span class=n>matplotlib</span> <span class=n>inline</span>
</span></span><span class=line><span class=cl><span class=n>nx</span> <span class=o>=</span> <span class=mi>41</span>
</span></span><span class=line><span class=cl><span class=n>ny</span> <span class=o>=</span> <span class=mi>41</span>
</span></span><span class=line><span class=cl><span class=n>nt</span> <span class=o>=</span> <span class=mi>500</span>
</span></span><span class=line><span class=cl><span class=n>nit</span> <span class=o>=</span> <span class=mi>50</span>
</span></span><span class=line><span class=cl><span class=n>c</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl><span class=n>dx</span> <span class=o>=</span> <span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=n>nx</span> <span class=o>-</span> <span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>dy</span> <span class=o>=</span> <span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=n>ny</span> <span class=o>-</span> <span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>x</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>linspace</span><span class=p>(</span><span class=mi>0</span><span class=p>,</span> <span class=mi>2</span><span class=p>,</span> <span class=n>nx</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>y</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>linspace</span><span class=p>(</span><span class=mi>0</span><span class=p>,</span> <span class=mi>2</span><span class=p>,</span> <span class=n>ny</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>X</span><span class=p>,</span> <span class=n>Y</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>meshgrid</span><span class=p>(</span><span class=n>x</span><span class=p>,</span> <span class=n>y</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>rho</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl><span class=n>nu</span> <span class=o>=</span> <span class=mf>.1</span>
</span></span><span class=line><span class=cl><span class=n>dt</span> <span class=o>=</span> <span class=mf>.001</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>u</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>v</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>p</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>b</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=k>def</span> <span class=nf>build_up_b</span><span class=p>(</span><span class=n>b</span><span class=p>,</span> <span class=n>rho</span><span class=p>,</span> <span class=n>dt</span><span class=p>,</span> <span class=n>u</span><span class=p>,</span> <span class=n>v</span><span class=p>,</span> <span class=n>dx</span><span class=p>,</span> <span class=n>dy</span><span class=p>):</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=n>b</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=p>(</span><span class=n>rho</span> <span class=o>*</span> <span class=p>(</span><span class=mi>1</span> <span class=o>/</span> <span class=n>dt</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>((</span><span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>2</span><span class=p>:]</span> <span class=o>-</span> <span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>/</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dx</span><span class=p>)</span> <span class=o>+</span> <span class=p>(</span><span class=n>v</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>v</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dy</span><span class=p>))</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=p>((</span><span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>2</span><span class=p>:]</span> <span class=o>-</span> <span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dx</span><span class=p>))</span><span class=o>**</span><span class=mi>2</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=mi>2</span> <span class=o>*</span> <span class=p>((</span><span class=n>u</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>u</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dy</span><span class=p>)</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>v</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>2</span><span class=p>:]</span> <span class=o>-</span> <span class=n>v</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dx</span><span class=p>))</span><span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=p>((</span><span class=n>v</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>v</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dy</span><span class=p>))</span><span class=o>**</span><span class=mi>2</span><span class=p>))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=k>return</span> <span class=n>b</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=k>def</span> <span class=nf>pressure_poisson</span><span class=p>(</span><span class=n>p</span><span class=p>,</span> <span class=n>dx</span><span class=p>,</span> <span class=n>dy</span><span class=p>,</span> <span class=n>b</span><span class=p>):</span>
</span></span><span class=line><span class=cl> <span class=n>pn</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>empty_like</span><span class=p>(</span><span class=n>p</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>pn</span> <span class=o>=</span> <span class=n>p</span><span class=o>.</span><span class=n>copy</span><span class=p>()</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=k>for</span> <span class=n>q</span> <span class=ow>in</span> <span class=nb>range</span><span class=p>(</span><span class=n>nit</span><span class=p>):</span>
</span></span><span class=line><span class=cl> <span class=n>pn</span> <span class=o>=</span> <span class=n>p</span><span class=o>.</span><span class=n>copy</span><span class=p>()</span>
</span></span><span class=line><span class=cl> <span class=n>p</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=p>(((</span><span class=n>pn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>2</span><span class=p>:]</span> <span class=o>+</span> <span class=n>pn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>*</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>pn</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>pn</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>*</span> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span><span class=p>)</span> <span class=o>/</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=p>(</span><span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>+</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span><span class=p>))</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=p>(</span><span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>+</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span><span class=p>))</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=n>b</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=n>p</span><span class=p>[:,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=n>p</span><span class=p>[:,</span> <span class=o>-</span><span class=mi>2</span><span class=p>]</span> <span class=c1># dp/dx = 0 at x = 2</span>
</span></span><span class=line><span class=cl> <span class=n>p</span><span class=p>[</span><span class=mi>0</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=n>p</span><span class=p>[</span><span class=mi>1</span><span class=p>,</span> <span class=p>:]</span> <span class=c1># dp/dy = 0 at y = 0</span>
</span></span><span class=line><span class=cl> <span class=n>p</span><span class=p>[:,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>=</span> <span class=n>p</span><span class=p>[:,</span> <span class=mi>1</span><span class=p>]</span> <span class=c1># dp/dx = 0 at x = 0</span>
</span></span><span class=line><span class=cl> <span class=n>p</span><span class=p>[</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=mi>0</span> <span class=c1># p = 0 at y = 2</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=k>return</span> <span class=n>p</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=k>def</span> <span class=nf>cavity_flow</span><span class=p>(</span><span class=n>nt</span><span class=p>,</span> <span class=n>u</span><span class=p>,</span> <span class=n>v</span><span class=p>,</span> <span class=n>dt</span><span class=p>,</span> <span class=n>dx</span><span class=p>,</span> <span class=n>dy</span><span class=p>,</span> <span class=n>p</span><span class=p>,</span> <span class=n>rho</span><span class=p>,</span> <span class=n>nu</span><span class=p>):</span>
</span></span><span class=line><span class=cl> <span class=n>un</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>empty_like</span><span class=p>(</span><span class=n>u</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>vn</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>empty_like</span><span class=p>(</span><span class=n>v</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>b</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=k>for</span> <span class=n>n</span> <span class=ow>in</span> <span class=nb>range</span><span class=p>(</span><span class=n>nt</span><span class=p>):</span>
</span></span><span class=line><span class=cl> <span class=n>un</span> <span class=o>=</span> <span class=n>u</span><span class=o>.</span><span class=n>copy</span><span class=p>()</span>
</span></span><span class=line><span class=cl> <span class=n>vn</span> <span class=o>=</span> <span class=n>v</span><span class=o>.</span><span class=n>copy</span><span class=p>()</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=n>b</span> <span class=o>=</span> <span class=n>build_up_b</span><span class=p>(</span><span class=n>b</span><span class=p>,</span> <span class=n>rho</span><span class=p>,</span> <span class=n>dt</span><span class=p>,</span> <span class=n>u</span><span class=p>,</span> <span class=n>v</span><span class=p>,</span> <span class=n>dx</span><span class=p>,</span> <span class=n>dy</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>p</span> <span class=o>=</span> <span class=n>pressure_poisson</span><span class=p>(</span><span class=n>p</span><span class=p>,</span> <span class=n>dx</span><span class=p>,</span> <span class=n>dy</span><span class=p>,</span> <span class=n>b</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span><span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dy</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>un</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>dt</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>rho</span> <span class=o>*</span> <span class=n>dx</span><span class=p>)</span> <span class=o>*</span> <span class=p>(</span><span class=n>p</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>2</span><span class=p>:]</span> <span class=o>-</span> <span class=n>p</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>nu</span> <span class=o>*</span> <span class=p>(</span><span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>2</span><span class=p>:]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>dt</span> <span class=o>/</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>un</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=n>v</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dy</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>vn</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>dt</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>rho</span> <span class=o>*</span> <span class=n>dy</span><span class=p>)</span> <span class=o>*</span> <span class=p>(</span><span class=n>p</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>p</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>nu</span> <span class=o>*</span> <span class=p>(</span><span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>2</span><span class=p>:]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>dt</span> <span class=o>/</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>vn</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[</span><span class=mi>0</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=mi>0</span>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[:,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>=</span> <span class=mi>0</span>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[:,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=mi>0</span>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=mi>1</span> <span class=c1># set velocity on cavity lid equal to 1</span>
</span></span><span class=line><span class=cl> <span class=n>v</span><span class=p>[</span><span class=mi>0</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=mi>0</span>
</span></span><span class=line><span class=cl> <span class=n>v</span><span class=p>[</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=mi>0</span>
</span></span><span class=line><span class=cl> <span class=n>v</span><span class=p>[:,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>=</span> <span class=mi>0</span>
</span></span><span class=line><span class=cl> <span class=n>v</span><span class=p>[:,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=mi>0</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=k>return</span> <span class=n>u</span><span class=p>,</span> <span class=n>v</span><span class=p>,</span> <span class=n>p</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>u</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>v</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>p</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>b</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>nt</span> <span class=o>=</span> <span class=mi>100</span>
</span></span><span class=line><span class=cl><span class=n>u</span><span class=p>,</span> <span class=n>v</span><span class=p>,</span> <span class=n>p</span> <span class=o>=</span> <span class=n>cavity_flow</span><span class=p>(</span><span class=n>nt</span><span class=p>,</span> <span class=n>u</span><span class=p>,</span> <span class=n>v</span><span class=p>,</span> <span class=n>dt</span><span class=p>,</span> <span class=n>dx</span><span class=p>,</span> <span class=n>dy</span><span class=p>,</span> <span class=n>p</span><span class=p>,</span> <span class=n>rho</span><span class=p>,</span> <span class=n>nu</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>fig</span> <span class=o>=</span> <span class=n>pyplot</span><span class=o>.</span><span class=n>figure</span><span class=p>(</span><span class=n>figsize</span><span class=o>=</span><span class=p>(</span><span class=mi>11</span><span class=p>,</span><span class=mi>7</span><span class=p>),</span> <span class=n>dpi</span><span class=o>=</span><span class=mi>100</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=c1># plotting the pressure field as a contour</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>contourf</span><span class=p>(</span><span class=n>X</span><span class=p>,</span> <span class=n>Y</span><span class=p>,</span> <span class=n>p</span><span class=p>,</span> <span class=n>alpha</span><span class=o>=</span><span class=mf>0.5</span><span class=p>,</span> <span class=n>cmap</span><span class=o>=</span><span class=n>cm</span><span class=o>.</span><span class=n>viridis</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>colorbar</span><span class=p>()</span>
</span></span><span class=line><span class=cl><span class=c1># plotting the pressure field outlines</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>contour</span><span class=p>(</span><span class=n>X</span><span class=p>,</span> <span class=n>Y</span><span class=p>,</span> <span class=n>p</span><span class=p>,</span> <span class=n>cmap</span><span class=o>=</span><span class=n>cm</span><span class=o>.</span><span class=n>viridis</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=c1># plotting velocity field</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>quiver</span><span class=p>(</span><span class=n>X</span><span class=p>[::</span><span class=mi>2</span><span class=p>,</span> <span class=p>::</span><span class=mi>2</span><span class=p>],</span> <span class=n>Y</span><span class=p>[::</span><span class=mi>2</span><span class=p>,</span> <span class=p>::</span><span class=mi>2</span><span class=p>],</span> <span class=n>u</span><span class=p>[::</span><span class=mi>2</span><span class=p>,</span> <span class=p>::</span><span class=mi>2</span><span class=p>],</span> <span class=n>v</span><span class=p>[::</span><span class=mi>2</span><span class=p>,</span> <span class=p>::</span><span class=mi>2</span><span class=p>])</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>xlabel</span><span class=p>(</span><span class=s1>&#39;X&#39;</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>ylabel</span><span class=p>(</span><span class=s1>&#39;Y&#39;</span><span class=p>);</span>
</span></span></code></pre></div><figure><img loading=lazy src=https://blog.kakaocdn.net/dn/VgtZK/btq9flvNFsA/4d03urU7VcLPRqzS5g40m1/img.png></figure><div class=highlight><pre tabindex=0 class=chroma><code class=language-python data-lang=python><span class=line><span class=cl><span class=n>u</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>v</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>p</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>b</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>nt</span> <span class=o>=</span> <span class=mi>700</span>
</span></span><span class=line><span class=cl><span class=n>u</span><span class=p>,</span> <span class=n>v</span><span class=p>,</span> <span class=n>p</span> <span class=o>=</span> <span class=n>cavity_flow</span><span class=p>(</span><span class=n>nt</span><span class=p>,</span> <span class=n>u</span><span class=p>,</span> <span class=n>v</span><span class=p>,</span> <span class=n>dt</span><span class=p>,</span> <span class=n>dx</span><span class=p>,</span> <span class=n>dy</span><span class=p>,</span> <span class=n>p</span><span class=p>,</span> <span class=n>rho</span><span class=p>,</span> <span class=n>nu</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>fig</span> <span class=o>=</span> <span class=n>pyplot</span><span class=o>.</span><span class=n>figure</span><span class=p>(</span><span class=n>figsize</span><span class=o>=</span><span class=p>(</span><span class=mi>11</span><span class=p>,</span> <span class=mi>7</span><span class=p>),</span> <span class=n>dpi</span><span class=o>=</span><span class=mi>100</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>contourf</span><span class=p>(</span><span class=n>X</span><span class=p>,</span> <span class=n>Y</span><span class=p>,</span> <span class=n>p</span><span class=p>,</span> <span class=n>alpha</span><span class=o>=</span><span class=mf>0.5</span><span class=p>,</span> <span class=n>cmap</span><span class=o>=</span><span class=n>cm</span><span class=o>.</span><span class=n>viridis</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>colorbar</span><span class=p>()</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>contour</span><span class=p>(</span><span class=n>X</span><span class=p>,</span> <span class=n>Y</span><span class=p>,</span> <span class=n>p</span><span class=p>,</span> <span class=n>cmap</span><span class=o>=</span><span class=n>cm</span><span class=o>.</span><span class=n>viridis</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>quiver</span><span class=p>(</span><span class=n>X</span><span class=p>[::</span><span class=mi>2</span><span class=p>,</span> <span class=p>::</span><span class=mi>2</span><span class=p>],</span> <span class=n>Y</span><span class=p>[::</span><span class=mi>2</span><span class=p>,</span> <span class=p>::</span><span class=mi>2</span><span class=p>],</span> <span class=n>u</span><span class=p>[::</span><span class=mi>2</span><span class=p>,</span> <span class=p>::</span><span class=mi>2</span><span class=p>],</span> <span class=n>v</span><span class=p>[::</span><span class=mi>2</span><span class=p>,</span> <span class=p>::</span><span class=mi>2</span><span class=p>])</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>xlabel</span><span class=p>(</span><span class=s1>&#39;X&#39;</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>ylabel</span><span class=p>(</span><span class=s1>&#39;Y&#39;</span><span class=p>);</span>
</span></span></code></pre></div><figure><img loading=lazy src=https://blog.kakaocdn.net/dn/boYcRB/btq9fmasDP2/6w1UPPVU4mG7frDcjQtBIK/img.png></figure><div class=highlight><pre tabindex=0 class=chroma><code class=language-python data-lang=python><span class=line><span class=cl><span class=n>fig</span> <span class=o>=</span> <span class=n>pyplot</span><span class=o>.</span><span class=n>figure</span><span class=p>(</span><span class=n>figsize</span><span class=o>=</span><span class=p>(</span><span class=mi>11</span><span class=p>,</span> <span class=mi>7</span><span class=p>),</span> <span class=n>dpi</span><span class=o>=</span><span class=mi>100</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>contourf</span><span class=p>(</span><span class=n>X</span><span class=p>,</span> <span class=n>Y</span><span class=p>,</span> <span class=n>p</span><span class=p>,</span> <span class=n>alpha</span><span class=o>=</span><span class=mf>0.5</span><span class=p>,</span> <span class=n>cmap</span><span class=o>=</span><span class=n>cm</span><span class=o>.</span><span class=n>viridis</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>colorbar</span><span class=p>()</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>contour</span><span class=p>(</span><span class=n>X</span><span class=p>,</span> <span class=n>Y</span><span class=p>,</span> <span class=n>p</span><span class=p>,</span> <span class=n>cmap</span><span class=o>=</span><span class=n>cm</span><span class=o>.</span><span class=n>viridis</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>streamplot</span><span class=p>(</span><span class=n>X</span><span class=p>,</span> <span class=n>Y</span><span class=p>,</span> <span class=n>u</span><span class=p>,</span> <span class=n>v</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>xlabel</span><span class=p>(</span><span class=s1>&#39;X&#39;</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>ylabel</span><span class=p>(</span><span class=s1>&#39;Y&#39;</span><span class=p>);</span>
</span></span></code></pre></div><figure><img loading=lazy src=https://blog.kakaocdn.net/dn/0J2aJ/btq9e57PmFB/UqArsnX9hzJ84H4rW5AtB1/img.png></figure><h2 id=channel-flow-with-navierstokes>Channel Flow with Navier–Stokes<a hidden class=anchor aria-hidden=true href=#channel-flow-with-navierstokes>#</a></h2><p>$$ \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=-\frac{1}{\rho}\frac{\partial p}{\partial x}+\nu\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}\right)+F $$$$ \frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2}=-\rho\left(\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+2\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+\frac{\partial v}{\partial y}\frac{\partial v}{\partial y}\right) $$$$ \begin{split}p_{i,j}^{n} = & \frac{\left(p_{i+1,j}^{n}+p_{i-1,j}^{n}\right) \Delta y^2 + \left(p_{i,j+1}^{n}+p_{i,j-1}^{n}\right) \Delta x^2}{2(\Delta x^2+\Delta y^2)} \& -\frac{\rho\Delta x^2\Delta y^2}{2\left(\Delta x^2+\Delta y^2\right)} \& \times \left[\frac{1}{\Delta t} \left(\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x} + \frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right) - \frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x} - 2\frac{u_{i,j+1}-u_{i,j-1}}{2\Delta y}\frac{v_{i+1,j}-v_{i-1,j}}{2\Delta x} - \frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right]\end{split} $$</p><div class=highlight><pre tabindex=0 class=chroma><code class=language-python data-lang=python><span class=line><span class=cl><span class=kn>import</span> <span class=nn>numpy</span>
</span></span><span class=line><span class=cl><span class=kn>from</span> <span class=nn>matplotlib</span> <span class=kn>import</span> <span class=n>pyplot</span><span class=p>,</span> <span class=n>cm</span>
</span></span><span class=line><span class=cl><span class=kn>from</span> <span class=nn>mpl_toolkits.mplot3d</span> <span class=kn>import</span> <span class=n>Axes3D</span>
</span></span><span class=line><span class=cl><span class=o>%</span><span class=n>matplotlib</span> <span class=n>inline</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=k>def</span> <span class=nf>build_up_b</span><span class=p>(</span><span class=n>rho</span><span class=p>,</span> <span class=n>dt</span><span class=p>,</span> <span class=n>dx</span><span class=p>,</span> <span class=n>dy</span><span class=p>,</span> <span class=n>u</span><span class=p>,</span> <span class=n>v</span><span class=p>):</span>
</span></span><span class=line><span class=cl> <span class=n>b</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros_like</span><span class=p>(</span><span class=n>u</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>b</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=p>(</span><span class=n>rho</span> <span class=o>*</span> <span class=p>(</span><span class=mi>1</span> <span class=o>/</span> <span class=n>dt</span> <span class=o>*</span> <span class=p>((</span><span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>2</span><span class=p>:]</span> <span class=o>-</span> <span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dx</span><span class=p>)</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>v</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>v</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dy</span><span class=p>))</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=p>((</span><span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>2</span><span class=p>:]</span> <span class=o>-</span> <span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dx</span><span class=p>))</span><span class=o>**</span><span class=mi>2</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=mi>2</span> <span class=o>*</span> <span class=p>((</span><span class=n>u</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>u</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dy</span><span class=p>)</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>v</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>2</span><span class=p>:]</span> <span class=o>-</span> <span class=n>v</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dx</span><span class=p>))</span><span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=p>((</span><span class=n>v</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>v</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dy</span><span class=p>))</span><span class=o>**</span><span class=mi>2</span><span class=p>))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=c1># Periodic BC Pressure @ x = 2</span>
</span></span><span class=line><span class=cl> <span class=n>b</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=p>(</span><span class=n>rho</span> <span class=o>*</span> <span class=p>(</span><span class=mi>1</span> <span class=o>/</span> <span class=n>dt</span> <span class=o>*</span> <span class=p>((</span><span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>-</span> <span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dx</span><span class=p>)</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>v</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>v</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dy</span><span class=p>))</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=p>((</span><span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>-</span> <span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dx</span><span class=p>))</span><span class=o>**</span><span class=mi>2</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=mi>2</span> <span class=o>*</span> <span class=p>((</span><span class=n>u</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>u</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dy</span><span class=p>)</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>v</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>-</span> <span class=n>v</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dx</span><span class=p>))</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=p>((</span><span class=n>v</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>v</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dy</span><span class=p>))</span><span class=o>**</span><span class=mi>2</span><span class=p>))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=c1># Periodic BC Pressure @ x = 0</span>
</span></span><span class=line><span class=cl> <span class=n>b</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>=</span> <span class=p>(</span><span class=n>rho</span> <span class=o>*</span> <span class=p>(</span><span class=mi>1</span> <span class=o>/</span> <span class=n>dt</span> <span class=o>*</span> <span class=p>((</span><span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dx</span><span class=p>)</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>v</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>-</span> <span class=n>v</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>0</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dy</span><span class=p>))</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=p>((</span><span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dx</span><span class=p>))</span><span class=o>**</span><span class=mi>2</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=mi>2</span> <span class=o>*</span> <span class=p>((</span><span class=n>u</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>-</span> <span class=n>u</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>0</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dy</span><span class=p>)</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>v</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>v</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dx</span><span class=p>))</span><span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=p>((</span><span class=n>v</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>-</span> <span class=n>v</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>0</span><span class=p>])</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dy</span><span class=p>))</span><span class=o>**</span><span class=mi>2</span><span class=p>))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=k>return</span> <span class=n>b</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=k>def</span> <span class=nf>pressure_poisson_periodic</span><span class=p>(</span><span class=n>p</span><span class=p>,</span> <span class=n>dx</span><span class=p>,</span> <span class=n>dy</span><span class=p>):</span>
</span></span><span class=line><span class=cl> <span class=n>pn</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>empty_like</span><span class=p>(</span><span class=n>p</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=k>for</span> <span class=n>q</span> <span class=ow>in</span> <span class=nb>range</span><span class=p>(</span><span class=n>nit</span><span class=p>):</span>
</span></span><span class=line><span class=cl> <span class=n>pn</span> <span class=o>=</span> <span class=n>p</span><span class=o>.</span><span class=n>copy</span><span class=p>()</span>
</span></span><span class=line><span class=cl> <span class=n>p</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=p>(((</span><span class=n>pn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>2</span><span class=p>:]</span> <span class=o>+</span> <span class=n>pn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>*</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>pn</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>pn</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>*</span> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span><span class=p>)</span> <span class=o>/</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=p>(</span><span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>+</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span><span class=p>))</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=p>(</span><span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>+</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span><span class=p>))</span> <span class=o>*</span> <span class=n>b</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=c1># Periodic BC Pressure @ x = 2</span>
</span></span><span class=line><span class=cl> <span class=n>p</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=p>(((</span><span class=n>pn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>+</span> <span class=n>pn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>2</span><span class=p>])</span><span class=o>*</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>pn</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>pn</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>*</span> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span><span class=p>)</span> <span class=o>/</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=p>(</span><span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>+</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span><span class=p>))</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=p>(</span><span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>+</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span><span class=p>))</span> <span class=o>*</span> <span class=n>b</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>])</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=c1># Periodic BC Pressure @ x = 0</span>
</span></span><span class=line><span class=cl> <span class=n>p</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>=</span> <span class=p>(((</span><span class=n>pn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>pn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>])</span><span class=o>*</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>pn</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>+</span> <span class=n>pn</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>0</span><span class=p>])</span> <span class=o>*</span> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span><span class=p>)</span> <span class=o>/</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=p>(</span><span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>+</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span><span class=p>))</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=p>(</span><span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>+</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span><span class=p>))</span> <span class=o>*</span> <span class=n>b</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>])</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=c1># Wall boundary conditions, pressure</span>
</span></span><span class=line><span class=cl> <span class=n>p</span><span class=p>[</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span><span class=n>p</span><span class=p>[</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=p>:]</span> <span class=c1># dp/dy = 0 at y = 2</span>
</span></span><span class=line><span class=cl> <span class=n>p</span><span class=p>[</span><span class=mi>0</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=n>p</span><span class=p>[</span><span class=mi>1</span><span class=p>,</span> <span class=p>:]</span> <span class=c1># dp/dy = 0 at y = 0</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=k>return</span> <span class=n>p</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=c1>##variable declarations</span>
</span></span><span class=line><span class=cl><span class=n>nx</span> <span class=o>=</span> <span class=mi>41</span>
</span></span><span class=line><span class=cl><span class=n>ny</span> <span class=o>=</span> <span class=mi>41</span>
</span></span><span class=line><span class=cl><span class=n>nt</span> <span class=o>=</span> <span class=mi>10</span>
</span></span><span class=line><span class=cl><span class=n>nit</span> <span class=o>=</span> <span class=mi>50</span>
</span></span><span class=line><span class=cl><span class=n>c</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl><span class=n>dx</span> <span class=o>=</span> <span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=n>nx</span> <span class=o>-</span> <span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>dy</span> <span class=o>=</span> <span class=mi>2</span> <span class=o>/</span> <span class=p>(</span><span class=n>ny</span> <span class=o>-</span> <span class=mi>1</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>x</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>linspace</span><span class=p>(</span><span class=mi>0</span><span class=p>,</span> <span class=mi>2</span><span class=p>,</span> <span class=n>nx</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>y</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>linspace</span><span class=p>(</span><span class=mi>0</span><span class=p>,</span> <span class=mi>2</span><span class=p>,</span> <span class=n>ny</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>X</span><span class=p>,</span> <span class=n>Y</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>meshgrid</span><span class=p>(</span><span class=n>x</span><span class=p>,</span> <span class=n>y</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=c1>##physical variables</span>
</span></span><span class=line><span class=cl><span class=n>rho</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl><span class=n>nu</span> <span class=o>=</span> <span class=mf>.1</span>
</span></span><span class=line><span class=cl><span class=n>F</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl><span class=n>dt</span> <span class=o>=</span> <span class=mf>.01</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=c1>#initial conditions</span>
</span></span><span class=line><span class=cl><span class=n>u</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>un</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>v</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>vn</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>p</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>ones</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl><span class=n>pn</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>ones</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>b</span> <span class=o>=</span> <span class=n>numpy</span><span class=o>.</span><span class=n>zeros</span><span class=p>((</span><span class=n>ny</span><span class=p>,</span> <span class=n>nx</span><span class=p>))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>udiff</span> <span class=o>=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl><span class=n>stepcount</span> <span class=o>=</span> <span class=mi>0</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=k>while</span> <span class=n>udiff</span> <span class=o>&gt;</span> <span class=mf>.001</span><span class=p>:</span>
</span></span><span class=line><span class=cl> <span class=n>un</span> <span class=o>=</span> <span class=n>u</span><span class=o>.</span><span class=n>copy</span><span class=p>()</span>
</span></span><span class=line><span class=cl> <span class=n>vn</span> <span class=o>=</span> <span class=n>v</span><span class=o>.</span><span class=n>copy</span><span class=p>()</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=n>b</span> <span class=o>=</span> <span class=n>build_up_b</span><span class=p>(</span><span class=n>rho</span><span class=p>,</span> <span class=n>dt</span><span class=p>,</span> <span class=n>dx</span><span class=p>,</span> <span class=n>dy</span><span class=p>,</span> <span class=n>u</span><span class=p>,</span> <span class=n>v</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>p</span> <span class=o>=</span> <span class=n>pressure_poisson_periodic</span><span class=p>(</span><span class=n>p</span><span class=p>,</span> <span class=n>dx</span><span class=p>,</span> <span class=n>dy</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dy</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>un</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>dt</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>rho</span> <span class=o>*</span> <span class=n>dx</span><span class=p>)</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>p</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>2</span><span class=p>:]</span> <span class=o>-</span> <span class=n>p</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>nu</span> <span class=o>*</span> <span class=p>(</span><span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>2</span><span class=p>:]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>dt</span> <span class=o>/</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>un</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]))</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>F</span> <span class=o>*</span> <span class=n>dt</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=n>v</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dy</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>vn</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>dt</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>rho</span> <span class=o>*</span> <span class=n>dy</span><span class=p>)</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>p</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>p</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>nu</span> <span class=o>*</span> <span class=p>(</span><span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>2</span><span class=p>:]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>dt</span> <span class=o>/</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>vn</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>])))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=c1># Periodic BC u @ x = 2 </span>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dy</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>un</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>dt</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>rho</span> <span class=o>*</span> <span class=n>dx</span><span class=p>)</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>p</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>-</span> <span class=n>p</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>nu</span> <span class=o>*</span> <span class=p>(</span><span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span><span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>dt</span> <span class=o>/</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>un</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]))</span> <span class=o>+</span> <span class=n>F</span> <span class=o>*</span> <span class=n>dt</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=c1># Periodic BC u @ x = 0</span>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>=</span> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>-</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>-</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dy</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>-</span> <span class=n>un</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>0</span><span class=p>])</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>dt</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>rho</span> <span class=o>*</span> <span class=n>dx</span><span class=p>)</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>p</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>p</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>nu</span> <span class=o>*</span> <span class=p>(</span><span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>+</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>dt</span> <span class=o>/</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>un</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>+</span> <span class=n>un</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>0</span><span class=p>]))</span> <span class=o>+</span> <span class=n>F</span> <span class=o>*</span> <span class=n>dt</span><span class=p>)</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=c1># Periodic BC v @ x = 2</span>
</span></span><span class=line><span class=cl> <span class=n>v</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>=</span> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dy</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>vn</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>dt</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>rho</span> <span class=o>*</span> <span class=n>dy</span><span class=p>)</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>p</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=n>p</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>nu</span> <span class=o>*</span> <span class=p>(</span><span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>2</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>dt</span> <span class=o>/</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>]</span> <span class=o>+</span> <span class=n>vn</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>])))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=c1># Periodic BC v @ x = 0</span>
</span></span><span class=line><span class=cl> <span class=n>v</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>=</span> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>-</span> <span class=n>un</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>-</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>*</span> <span class=n>dt</span> <span class=o>/</span> <span class=n>dy</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>-</span> <span class=n>vn</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>0</span><span class=p>])</span> <span class=o>-</span>
</span></span><span class=line><span class=cl> <span class=n>dt</span> <span class=o>/</span> <span class=p>(</span><span class=mi>2</span> <span class=o>*</span> <span class=n>rho</span> <span class=o>*</span> <span class=n>dy</span><span class=p>)</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>p</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>-</span> <span class=n>p</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>0</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>nu</span> <span class=o>*</span> <span class=p>(</span><span class=n>dt</span> <span class=o>/</span> <span class=n>dx</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>1</span><span class=p>]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>+</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=o>-</span><span class=mi>1</span><span class=p>])</span> <span class=o>+</span>
</span></span><span class=line><span class=cl> <span class=n>dt</span> <span class=o>/</span> <span class=n>dy</span><span class=o>**</span><span class=mi>2</span> <span class=o>*</span>
</span></span><span class=line><span class=cl> <span class=p>(</span><span class=n>vn</span><span class=p>[</span><span class=mi>2</span><span class=p>:,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>-</span> <span class=mi>2</span> <span class=o>*</span> <span class=n>vn</span><span class=p>[</span><span class=mi>1</span><span class=p>:</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=mi>0</span><span class=p>]</span> <span class=o>+</span> <span class=n>vn</span><span class=p>[</span><span class=mi>0</span><span class=p>:</span><span class=o>-</span><span class=mi>2</span><span class=p>,</span> <span class=mi>0</span><span class=p>])))</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=c1># Wall BC: u,v = 0 @ y = 0,2</span>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[</span><span class=mi>0</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=mi>0</span>
</span></span><span class=line><span class=cl> <span class=n>u</span><span class=p>[</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=mi>0</span>
</span></span><span class=line><span class=cl> <span class=n>v</span><span class=p>[</span><span class=mi>0</span><span class=p>,</span> <span class=p>:]</span> <span class=o>=</span> <span class=mi>0</span>
</span></span><span class=line><span class=cl> <span class=n>v</span><span class=p>[</span><span class=o>-</span><span class=mi>1</span><span class=p>,</span> <span class=p>:]</span><span class=o>=</span><span class=mi>0</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl> <span class=n>udiff</span> <span class=o>=</span> <span class=p>(</span><span class=n>numpy</span><span class=o>.</span><span class=n>sum</span><span class=p>(</span><span class=n>u</span><span class=p>)</span> <span class=o>-</span> <span class=n>numpy</span><span class=o>.</span><span class=n>sum</span><span class=p>(</span><span class=n>un</span><span class=p>))</span> <span class=o>/</span> <span class=n>numpy</span><span class=o>.</span><span class=n>sum</span><span class=p>(</span><span class=n>u</span><span class=p>)</span>
</span></span><span class=line><span class=cl> <span class=n>stepcount</span> <span class=o>+=</span> <span class=mi>1</span>
</span></span><span class=line><span class=cl>
</span></span><span class=line><span class=cl><span class=n>fig</span> <span class=o>=</span> <span class=n>pyplot</span><span class=o>.</span><span class=n>figure</span><span class=p>(</span><span class=n>figsize</span> <span class=o>=</span> <span class=p>(</span><span class=mi>11</span><span class=p>,</span><span class=mi>7</span><span class=p>),</span> <span class=n>dpi</span><span class=o>=</span><span class=mi>100</span><span class=p>)</span>
</span></span><span class=line><span class=cl><span class=n>pyplot</span><span class=o>.</span><span class=n>quiver</span><span class=p>(</span><span class=n>X</span><span class=p>[::</span><span class=mi>3</span><span class=p>,</span> <span class=p>::</span><span class=mi>3</span><span class=p>],</span> <span class=n>Y</span><span class=p>[::</span><span class=mi>3</span><span class=p>,</span> <span class=p>::</span><span class=mi>3</span><span class=p>],</span> <span class=n>u</span><span class=p>[::</span><span class=mi>3</span><span class=p>,</span> <span class=p>::</span><span class=mi>3</span><span class=p>],</span> <span class=n>v</span><span class=p>[::</span><span class=mi>3</span><span class=p>,</span> <span class=p>::</span><span class=mi>3</span><span class=p>]);</span>
</span></span></code></pre></div><figure><img loading=lazy src=https://blog.kakaocdn.net/dn/du6hla/btq9fdKZP6o/ifKi67Tsr8khMmReNSHn5K/img.png></figure><p>출처> <a href=https://lorenabarba.com/blog/cfd-python-12-steps-to-navier-stokes/>CFD Python: 12 steps to Navier-Stokes :: Lorena A. Barba Group (lorenabarba.com)</a></p></div><footer class=post-footer><ul class=post-tags></ul><nav class=paginav><a class=prev href=http://blog.morgan.kr/posts/hardware-security/><span class=title>« Prev</span><br><span>Hardware Security</span></a>
<a class=next href=http://blog.morgan.kr/posts/bandoceyi-weonri-jongryu-soja-yeogsa-saneob-mosfetbuteo-ram-flash-geurigo-intel-4004ggaji/><span class=title>Next »</span><br><span>반도체의 원리, 종류, 소자, 역사, 산업. (MOSFET부터 RAM, FLASH, 그리고 Intel 4004까지.</span></a></nav><br></footer></article></main><footer class=footer><span>&copy; 2023 <a href=http://blog.morgan.kr>Morgan's Blog</a></span>
<span>Powered by
<a href=https://gohugo.io/ rel="noopener noreferrer" target=_blank>Hugo</a> &
<a href=https://github.com/adityatelange/hugo-PaperMod/ rel=noopener target=_blank>PaperMod</a></span></footer><a href=#top aria-label="go to top" title="Go to Top (Alt + G)" class=top-link id=top-link accesskey=g><svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 12 6" fill="currentcolor"><path d="M12 6H0l6-6z"/></svg></a><script>let menu=document.getElementById("menu");menu&&(menu.scrollLeft=localStorage.getItem("menu-scroll-position"),menu.onscroll=function(){localStorage.setItem("menu-scroll-position",menu.scrollLeft)}),document.querySelectorAll('a[href^="#"]').forEach(e=>{e.addEventListener("click",function(e){e.preventDefault();var t=this.getAttribute("href").substr(1);window.matchMedia("(prefers-reduced-motion: reduce)").matches?document.querySelector(`[id='${decodeURIComponent(t)}']`).scrollIntoView():document.querySelector(`[id='${decodeURIComponent(t)}']`).scrollIntoView({behavior:"smooth"}),t==="top"?history.replaceState(null,null," "):history.pushState(null,null,`#${t}`)})})</script><script>var mybutton=document.getElementById("top-link");window.onscroll=function(){document.body.scrollTop>800||document.documentElement.scrollTop>800?(mybutton.style.visibility="visible",mybutton.style.opacity="1"):(mybutton.style.visibility="hidden",mybutton.style.opacity="0")}</script><script>document.getElementById("theme-toggle").addEventListener("click",()=>{document.body.className.includes("dark")?(document.body.classList.remove("dark"),localStorage.setItem("pref-theme","light")):(document.body.classList.add("dark"),localStorage.setItem("pref-theme","dark"))})</script></body></html>