Fysiikan lakeja, joita virtaava neste tottelee

Kirjoitetaan fysiikan lakeja kuvaavien matemaattisten yhtälöiden perusteella python-funktio, jota tarvitaan numeerisessa simuloinnissa.

Jalo tarkoitukseni oli kirjoittaa tämä juttu niin, että se antaa yleisvaikutelman tämäntapaisesta ohjelmoinnista, vaikka koodiosuudet ja yhtälöt ohittaisi pikaisella silmäyksellä.

(Tällaiset ohjelmakoodit ovat yleensä vaikeita lukea. Luulenpa, että jos vuoden päästä palaan tämän koodin pariin, joudun pitkään hämmästelemään, mitä tässä tapahtuu ja miksi olen sen näin koodannut.)

Ohjelmoinnin perusteita olen yrittänyt esitellä toisaalla Ohjelmointia, matematiikkaa, fysiikka

---

Alla luetellaan valmiit python-modulit, joita tässä ohjelmassa käytetään.

Olennainen osa ohjelmointia on perehtyä, mitä moduleja on tarjolla ja miten niitä käytetään. Kokemuksen myötä oppii, millaisilla hauilla löytää webistä parhaat ohjeet ja esimerkit modulien soveltamiseen.

import datetime
import sympy as sp  # symbolisen laskennan moduli
import numpy as np # numeerisen laskennan moduli
from IPython.display import display  # Yhtälöiden tulostaminen
import re  # regular expression module, tarvitaan koodin generoinnissa
from textwrap import fill, wrap  # samaan tarkoitukseen

import params as par  # Nesteen ja järven ominaisuuksia, itse kirjoittamani 'moduli'

Allaolevia määrittelyjä tarvitaan symbolisessa laskennassa. Toisin kuin tavalliset muuttujat, symbolisen laskennan muuttujat pitää näin esitellä.

Aina kun alkaa ratkaista uudenlaista ongelmaa ohjelmoimalla, joutuu perehtymään uusiin asioihin. Uuden opettelu on olennainen osa ohjelmointia. Uuden opettelu toki edellyttää, että omaa perusymmärryksen ohjelmoinnista. Samalla lailla kuin matemaattinen ongelmanratkaisu edellyttää, että on omaksunut 'matemaattisen ajattelun', ohjelmoinnissa tarvitaan 'ohjelmoinnillista ajattelua'. Sitä on yhtä vaikea määritellä kuin matemaattista ajattelua.

---

Syrjähyppy menneeseen

Tietokoneeni käyttöliittymä on hyvin erilainen kuin se, mihin totuin -70 ja -80 -luvuilla.

Tämä on ohjelmoitu python-ohjelmointikielellä, jonka ensimmäinen versio julkaistiin 1991 ja, jolla opettelin ohjelmoimaan muutama vuosi sitten. Symboliseen laskentaan käyttämäni sympy-kirjaston ensimmäinen versio julkaistiin 2007. Tämä dokumentti on Jupyter-notebook, jonka on kehittänyt vuonna 2015 perustettu project Jupyter. Numeeriseen laskentaan käyttämäni numpy-kirjaston ensimmäinen versio on julkaistu 1995.

Kaikkia ylläolevia ja paljon muuta olen opetellut soveltamaan vasta viime vuosina harrastusmielessä. Kovin vaikeaa se ei ole ollut, koska olen tehnyt samanlaisia hommia -70 ja -80 -luvuilla. Minulle on kai silloin kehittynyt 'ohjelmoinnillinen ajattelu'.

Web-selaimen ja web-palvelimen välisen kommunikoinnin ohjelmointia tuskin jaksaisin opetella. Tietokannat, käyttäjien autentikointi web-sovelluksissa yms. ovat mysteereitä, joihin perehtymiseen en tunne houkutusta. Tuskin oppisin, koska siinä ei nuorena opittu auttaisi.

Tekoäly kuulostaa mielenkiintoiselta, mutta pelkkä ajatus opetella sen soveltamista alkaa ahdistaa. Meille opetettiin neuroverkkojen periaatteita -70 -luvulla, mutta niillä eväin ei pitkälle pötkitä.

x, j, y, i, t, n, c_p, mu, c_fr,  D_pnt, rho_0 = \
        sp.symbols('x j y i t n c_p \mu c_fr D_pnt rho_0', real=True)
ds = sp.Symbol(r'\Delta s', real=True)
dt = sp.Symbol(r'\Delta t', real=True)
rho = sp.Function('rho', real=True, nargs=3)
u = sp.Function('u', real=True, nargs=3)
v = sp.Function('v', real=True, nargs=3)
Fe_v = sp.Function('Fe_v', real=True, nargs=2)
Fe_u = sp.Function('Fe_u', real=True, nargs=2)
paint = sp.Function('paint', real=True, nargs=3)

rho_d = sp.Function('rho_d', real=True, nargs=2)
u_d = sp.Function('u_d', real=True, nargs=2)
v_d = sp.Function('v_d', real=True, nargs=2)
paint_d = sp.Function('paint_d', real=True, nargs=2)

D_u = sp.Function('D_u', real=True, nargs=2)
D_v = sp.Function('D_v', real=True, nargs=2)
D_rho = sp.Function('D_rho', real=True, nargs=2)
D_paint = sp.Function('D_paint', real=True, nargs=2)

# display-komennolla muuttujat saa näytettyä matemaatikoille tutummassa muodossa

display(x, mu, rho_0, Fe_v)

$\displaystyle x$

$\displaystyle \mu$

$\displaystyle \rho_{0}$

Fe_v

Symbolinen laskenta

Yhtälöitä on vaikea syöttää näppäimistöltä niin, että ne miellyttäisivät matemaatikon silmää. Siksi on kehitetty tapoja kirjoittaa yhtälöt helposti näppäiltävässä muodossa ja sen jälkeen erillisellä komennolla tulostaa ne muotoon, josta näkee, mitä tuli näppäiltyä.

Näissä kommenteissa käytän 80-luvulla kehitettyä LaTeX syntaksia. Esimerkiksi $$\nabla \cdot \overrightarrow{V} = \frac{\partial}{\partial x} u(x,y) + \frac{\partial}{\partial y} v(x,y)$$

kirjoitetaan

\nabla \cdot \overrightarrow{V} = \frac{\partial}{\partial x} u(x,y) + \frac{\partial}{\partial y} v(x,y)

Lähes kaikki matemaattisluonnontieteelliset tutkimusartikkelit kirjoitetaan LaTeXilla

LaTeX kuitenkin osaa pelkästään sommitella matemaattisia lausekkeita helposti luettavaan muotoon. Jotta yhtälöitä voisi käsitellä matematiikan sääntöjen mukaisesti, kirjoitan ne sympy-modulin vaatimaan muotoon.

(On olemassa muitakin symbolisen laskennan työkaluja. Sympyn esitystapa ei ehkä ole matemaattikkojen mielestä kaikista kaunein, mutta sympyn etuna on, että se on kiinteä osa python-ohjelmointikieltä.)

Tarvittaessa sympy-lausekkeet voi tulostaa helposti luettavassa muodossa display-komennolla, joka perustuu LaTeXiin:

esimerkki = sp.diff(u(x,y,t),x) + sp.diff(v(x,y,t),y)
display(esimerkki)

$\displaystyle \frac{\partial}{\partial x} u{\left(x,y,t \right)} + \frac{\partial}{\partial y} v{\left(x,y,t \right)}$

Kaasujen ja nesteiden (= fluidien) virtausten simuloinnin perusidea

Tarvittavien yhtälöiden johtamiseksi tarkastellaan pienen pientä fluidin mukana kulkevaa fluidi-kuutiota. Kuution tilan kertovat sen nopeus ja siinä olevan fluidin tiheys ja paine. (http://www.eng.auburn.edu/~tplacek/courses/fluidsreview-1.pdf)

Säilymislait

Tarpeeksi pitkään tuollaista kuution kuvaa tuijotettuaan pari tyyppiä päätyi kirjoittamaan yhtälöt kuution massan, liikemäärän ja energian säilymiselle ( navier-Stokes yhtälöt )

Koska on jokseenkin mahdoton tehdä ohjelmaa, jossa tuollaiset fluidi-elementit vaeltelevat ympäriinsä, matemaattiselle tempulla liikkuvan fluidin yhtälöt (non-conservation form) voidaan vääntää yhtälöiksi, jotka kertovat, mitä paikallaan olevassa kuutiossa tapahtuu (conservation form).

(Saattaa olla, etten ole kovin syvällisesti ymmärtänyt eri muotojen perimmäistä eroa. Sori siitä)

Navier-Stokes yhtälöt näyttävät niin luotaantyöntäviltä, että arvelin helpommaksi suoraan fysiikan peruslaeista lähtien kirjoittaa tarvittavat yhtälöt kuutioille, joista järveni rakennan. Muutaman päivän ja muutaman kymmenen yritys-erehdys kierroksen jälkeen nöyrryin selvittämään, miten Navier-Stokes yhtälöt on tapana ratkaista. Jälleen kerran sain muistutuksen siitä, ettei pahankaan näköinen yhtälö ole pohjimmiltaan ilkeä. Niiden kanssa tulee mukavasti juttuun, kunhan varaa kylliksi aikaa rauhalliseen tutustumiseen.

Massan säilyminen

Järvessäni massan säilyminen tarkoittaa, että vesi ei lisäänny eikä vähene, ellen laita puroja tuomaan tai viemään siitä vettä.

Miten noinkin pieni yhtälö voi näyttää noin ilkeältä? Minun matemaattinen intuitioni ei riitä niin pitkälle, että suoraan ymmärtäisin, mitä tuo tarkoittaa. Helpottaa, kun ymmärtää, että yhtälö on tiivistelmä pitkän johdattelun lopputuloksesta.

Toisaalta, näin taiteen tekemisessä ei haittaa, vaikkei intuitio kertoisikaan yhtälön merkitystä. Sen voi kopioida järveäni kuvaavien yhtälöiden sekaan ja antaa tietokoneen laskea, mitä siitä seuraa.

$\rho$: fluidin (= kaasun tai nesteen) tiheys, (järvessäni tämän voi tulkita vesikerroksen paksuudeksi.)

$\overrightarrow{V}$: fluidin nopeus. Nuoli tarkoittaa, että nopeudella on sekä suunta että suuruus

$\nabla: \overrightarrow{i}\frac{\partial}{\partial x} + \overrightarrow{j}\frac{\partial}{\partial y}$, lyhennysmerkintä derivoinnille paikan suhteen. Jos vuorenrinteellä korkeus merenpinnasta on h(x,y) niin $\nabla h$ on rinteen jyrkkyys

Jos $\overrightarrow{V} = u(x,y)\overrightarrow{i} + v(x,y)\overrightarrow{j}$ niin $\nabla \cdot \overrightarrow{V} = \frac{\partial}{\partial x} u(x,y) + \frac{\partial}{\partial y} v(x,y)$

Massan säilyminen sympy-muodossa

Kirjoitetaan yo. kuvan Eq. (2.27) sympy-muotoon

rho(x,y,t) ($\rho(,x,y,t$: Pinnankorkeus

u(x,y,t): veden nopeuden itä-länsi -suuntainen komponentti

v(x,y,t): veden nopeuden etelä-pohjois -suuntainen komponentti

# diffus_lqd = D_lqd*(sp.diff(rho(x,y,t),x,2) + sp.diff(rho(x,y,t),y,2))
#  Ylläolevan lisäämällä saisi mukaan nesteen diffuusion. 
# Kokeilin: Lisäsi simulointia hidastavaa laskentaa, 
# ei tuottanut jännittävämpiä tuloksia eikä tasannut algoritmin epätarkkuuden aiheuttamia kupruja
# kommentoin pois

eq_mass = sp.Eq(sp.diff(rho(x,y,t),t) + 
                sp.diff(rho(x,y,t)*u(x,y,t),x) + 
                sp.diff(rho(x,y,t)*v(x,y,t),y), 0)  #  - diffus_lqd, 0)

display(eq_mass)

$\displaystyle \rho{\left(x,y,t \right)} \frac{\partial}{\partial x} u{\left(x,y,t \right)} + \rho{\left(x,y,t \right)} \frac{\partial}{\partial y} v{\left(x,y,t \right)} + u{\left(x,y,t \right)} \frac{\partial}{\partial x} \rho{\left(x,y,t \right)} + v{\left(x,y,t \right)} \frac{\partial}{\partial y} \rho{\left(x,y,t \right)} + \frac{\partial}{\partial t} \rho{\left(x,y,t \right)} = 0$

'Väriaineiden' leviäminen nesteeseen

Simuloinnin alussa järveen roiskitaan kahta väriainetta. Lisäksi rannoilta liukenee veteen samoja väriaineita. Diffuusio levittää niitä veteen, vaikka vesi pysyisi paikallaan, mutta lisäksi värinaineet kulkevat nesteen virtauksen mukana.

Väriaineen massa on niin vähäinen, ettei se vaikuta nesteen virtaukseen.

Sekoittuessaan väriaineet neutraloivat toisiaan.

Diffuusio

Diffuusio on esimerkiksi sitä, että vaikka miten varovasti tiputtaisit tilkan maitoa teehesi, ajan myötä maito värjää teen kauttaaltaan. Jos tipautat maidon korkealta tai sekoitat teetä, maito sekoittuu nopeammin, koska diffuusion lisäksi nesteen virtaus levittää maitoa kuppiin.

---

Fick's second law tells how diffusion causes the concentration to change with respect to time. It is a partial differential equation which in one dimension reads:

$$ \frac {\partial \varphi }{\partial t} = D\,{\frac {\partial ^{2}\varphi }{\partial x^{2}}} $$

where

$ \varphi $ is the concentration;

$ \varphi = \varphi (x,t) $ is a function that depends on location x and time t

D is the diffusion coefficient

In two or more dimensions we must use the Laplacian Δ = ∇2, which generalises the second derivative, obtaining the equation

$$ \frac {\partial \varphi }{\partial t} = D\Delta \varphi $$

https://en.wikipedia.org/wiki/Fick%27s_laws_of_diffusion#Fick's_second_law

D_pnt: diffuusiokerroin paint(x,y,t): väriaineen konsentraatio

Simulaatiossa on vain yhtä väriainetta. Jos sen konsentraatio on negatiivinen, se piirretään vihertävänsävyisenä, jos positiivinen, punertavana. Kun punainen ja vihertävä kohtaavat ja sekoittuvat, summaksi tulee lähellä nollaa oleva luku, jonka väri on vaaleahko.

Alla Ficksin laki 2-ulotteisessa tapauksessa kirjoitettuna sympy-muodossa

# diffuusio
diffus_pnt = D_pnt*(sp.diff(paint(x,y,t),x,2) + 
                    sp.diff(paint(x,y,t),y,2))
display(diffus_pnt)

$\displaystyle D_{pnt} \left(\frac{\partial^{2}}{\partial x^{2}} \operatorname{paint}{\left(x,y,t \right)} + \frac{\partial^{2}}{\partial y^{2}} \operatorname{paint}{\left(x,y,t \right)}\right)$

Väriaineelle pätee sama massan säilymisen yhtälö kuin vedellekin lisättynä diffuusion vaikutuksella.

Olen koodannut tämän muutama viikko sitten. Miksi olen laittanut miinus-merkin diffuusion vaikutuksen eteen? Nyt allaoleva väriaineen massan säilymisen yhtälö kertoo, että pisteessä (x,y) virta vie ja diffuusio tuo väriainetta. Ehkä olen kopioinut yhtälöt mitään miettimättä. Mitä suotta miettimään, valmiita yhtälöitä ;-)

Entä jos vaihtaisi diffuusion merkin? Väriaine ei leviäisi itsestään vaan päinvastoin kerääntyisi pienten tiivistymien ympärille. Eli virta veisi, mutta diffuusio taistelisi vastaan. Pakko kokeilla!

(Kokeilinkin. Ei näyttänyt hyvältä)

# Väriaineen massan säilyminen 
eq_paint = sp.Eq(sp.diff(paint(x,y,t),t) + 
                 sp.diff(paint(x,y,t)*u(x,y,t),x) + 
                 sp.diff(paint(x,y,t)*v(x,y,t),y) - 
                 diffus_pnt, 0)
display(eq_paint)

$\displaystyle - D_{pnt} \left(\frac{\partial^{2}}{\partial x^{2}} \operatorname{paint}{\left(x,y,t \right)} + \frac{\partial^{2}}{\partial y^{2}} \operatorname{paint}{\left(x,y,t \right)}\right) + \operatorname{paint}{\left(x,y,t \right)} \frac{\partial}{\partial x} u{\left(x,y,t \right)} + \operatorname{paint}{\left(x,y,t \right)} \frac{\partial}{\partial y} v{\left(x,y,t \right)} + u{\left(x,y,t \right)} \frac{\partial}{\partial x} \operatorname{paint}{\left(x,y,t \right)} + v{\left(x,y,t \right)} \frac{\partial}{\partial y} \operatorname{paint}{\left(x,y,t \right)} + \frac{\partial}{\partial t} \operatorname{paint}{\left(x,y,t \right)} = 0$

Liikemäärän säilyminen

Kiinteälle kappaleelle pätee: $voima = massa\cdot kiihtyvyys$ eli kappaleen vauhti kiihtyy, kun siihen kohdistuu voima. Mitä suurempi kappaleen massa, sitä enemmän pitää työntää, että sen saa kunnolla vauhtiin.

Sama pätee fluideihin, mutta liikeyhtälöitä johdettaessa tarkastellaan pienen fluidikuution liikettä.

Fluidikuution kohdistuvat voimat

Seuraavat yhtälöt kertovat, miten viereiset kuutiot työntävät ja vetävät toisiaan ja millaisia voimia viereisten kuutioiden välinen 'kitka' aiheuttaa.

Jos ketä kiinnostaa, voi kaivaa webistä selityksen kuvalle ja yhtälöille

Pääosin ymmärrän selityksen, mutta silti tuntuu oudolta mallintaa kaasua tai nestettä kuutioina. Itse en olisi rohjennut muista ehdottaa. Vähän helpottaa, kun ajattelee, että kuutio on 'infinitesimally small'. Siis pienempi kuin pieni.

Fluidin liikemäärän säilymisen yhtälöt

Ensivilkaisulla allaolevan yhtälöt näyttivät niin pahoilta, että kuvittelin helpommaksi miettiä itse, mitä kuutioille tapahtuu. En kuitenkaan onnistunut, joten jouduin nöyrtymään ja perehtymään, mitä nuo yhtälöt tarkoittavat.

Eipä tuossa loppujen lopuksi paljoa tarvinnut ymmärtää. Kopioin yhtälöt ja ratkaisin ne :-)

Liike-energian säilyminen

Kaasut ovat kokoonpuristuvia eli massaa varastoituu tiheyden vaihteluihin. Kaasun paine, tiheys ja lämpötila riippuvat toisistaan, joten kaasuja simuloidessa pitää osata laskea, miten energiaa varastoituu lämmöksi. Sitä varten olisi pitänyt ottaa mukaan allaolevat yhtälöt.

En jaksanut perehtyä ja olisi ollut tyhmää tehdä 'järvivedestä' kokoonpuristuvaa,koska nesteet ovat yleensä käytännössä kokoonpuristumattomia, eli massaa ei varastoidu mihinkään.

Aloin kuitenkin epäillä, että taiteellisen vaikutelman kannalta ajatellen kokoonpuristumaton liemi saattaa käyttäytyä kovin tylsästi. Jos fluidi puristuu kokoon, muutokset etenevät aaltona. Kokoonpuristumattomassa nesteessä muutoksen vaikutus leviää silmänräpäyksessä koko systeemiin.

Ehkä ihan väärin luulin. Olisi kannattanut kokeilla yksinkertaisempaa ilmiötä, mutta halusin heti ekalla yrityksellä mukaan kaiken mahdollisen.

Keksin, että järvessä nestettä voi ajatella varastoituvan pinnankorkeuden vaihteluihin jolloin yhtälöissä esiintyvä paine voisi olla suoraan verrannollinen pinnankorkeuteen. Minun ei tarvitsisi opetella ymmärtämään ylläolevia energiayhtälöitä ja kuitenkin saisin simulointiin varastoitumisen!

Aluksi olin huolissani, toimiiko vesistöni vallan luonnottomasti, ellen vaadi energiaa säilymään. Koska kuitenkin yhtälöitä syntyi yhtä monta kuin tuntemattomia, matemaattisesti ajatellen tilanne oli ok. Toisinaan, kun luottaa matematiikan lakeihin, ei tarvitse ajatella eikä ymmärtää.

Virtausvastus, kitka

Jos kappaleeseen kohdistuu voima, sen vauhti kiihtyy äärettömäksi, ellei mikään voima sitä jarruta. Järvessäni vesi virtaa ikuisesti alamäkeen, joten ylinopeus oli vaarana. Aluksi veden virtausnopeus kasvoikin niin suureksi, että algoritmini sekosi.

Koska järveni on 2-dimensioinen, yhtälöissä ei esiinny ylä- eikä alapuolen vaikutusta. Nehän olisivat kolmatta dimensiota, jota ei kerta kaikkiaan ole olemassa. Lisäsin silti virtausnopeuteen verrannollisen jarruttavan voiman kuvaamaan järven pohjaa. Nopeus ei enää kasvanut liian suureksi.

Myöhemmin oivalsin 'loiventaa alamäkiä' ja tehdä vedestä 'paksumpaa' niin, että rantojen vastus riitti pitämään virtausnopeuden aisoissa.

Liikemäärän säilyminen sympy-muodossa

# momentum conservation, x direction

# 'paine' pinnankorkeuden funktiona
p = c_p*rho(x,y,t)

# viereisten kuutioiden aiheuttamat voimat
vb = -2*mu/3
tau_xx = vb*(sp.diff(u(x,y,t),x) + sp.diff(v(x,y,t),y)) + \
            2*mu*sp.diff(u(x,y,t),x)
tau_yy = vb*(sp.diff(u(x,y,t),x) + sp.diff(v(x,y,t),y)) + \
                2*mu*sp.diff(v(x,y,t),y)
tau_xy = mu*(sp.diff(v(x,y,t),x) + sp.diff(u(x,y,t),y))

# Järven pohjan aiheuttama kitka
# fr_x = -c_fr*u(x, y, t)

# Liikemäärän yhtälö, Fe_u on nesteeseen kohdistuva ulkoinen voima, esim. maan vetovoiman aiheuttama voima koskessa
eq_mom_x = sp.Eq(sp.diff(rho(x,y,t)*u(x,y,t),t) + 
                 sp.diff(rho(x,y,t)*u(x,y,t)**2,x) + 
                 sp.diff(rho(x,y,t)*u(x,y,t)*v(x,y,t),y), 
                 -sp.diff(p,x) + sp.diff(tau_xx,x) + 
                 sp.diff(tau_xy,y) + rho(x, y, t)*Fe_u(x,y))  # + fr_x 
display(eq_mom_x)

$\displaystyle 2 \rho{\left(x,y,t \right)} u{\left(x,y,t \right)} \frac{\partial}{\partial x} u{\left(x,y,t \right)} + \rho{\left(x,y,t \right)} u{\left(x,y,t \right)} \frac{\partial}{\partial y} v{\left(x,y,t \right)} + \rho{\left(x,y,t \right)} v{\left(x,y,t \right)} \frac{\partial}{\partial y} u{\left(x,y,t \right)} + \rho{\left(x,y,t \right)} \frac{\partial}{\partial t} u{\left(x,y,t \right)} + u^{2}{\left(x,y,t \right)} \frac{\partial}{\partial x} \rho{\left(x,y,t \right)} + u{\left(x,y,t \right)} v{\left(x,y,t \right)} \frac{\partial}{\partial y} \rho{\left(x,y,t \right)} + u{\left(x,y,t \right)} \frac{\partial}{\partial t} \rho{\left(x,y,t \right)} = - \frac{2 \mu \left(\frac{\partial^{2}}{\partial x^{2}} u{\left(x,y,t \right)} + \frac{\partial^{2}}{\partial y\partial x} v{\left(x,y,t \right)}\right)}{3} + \mu \left(\frac{\partial^{2}}{\partial y^{2}} u{\left(x,y,t \right)} + \frac{\partial^{2}}{\partial y\partial x} v{\left(x,y,t \right)}\right) + 2 \mu \frac{\partial^{2}}{\partial x^{2}} u{\left(x,y,t \right)} - c_{p} \frac{\partial}{\partial x} \rho{\left(x,y,t \right)} + \operatorname{Fe_{u}}{\left(x,y \right)} \rho{\left(x,y,t \right)}$

# momentum conservation, y direction

# fr_y = -c_fr*v(x, y, t)

eq_mom_y = sp.Eq(sp.diff(rho(x,y,t)*v(x,y,t),t) + 
                 sp.diff(rho(x,y,t)*u(x,y,t)*v(x,y,t),x) + 
                 sp.diff(rho(x,y,t)*v(x,y,t)**2,y), 
                -sp.diff(p,y) + sp.diff(tau_yy,y) + sp.diff(tau_xy,x) 
                 + rho(x, y, t)*Fe_v(x,y))  # + fr_y 
display(eq_mom_y)

$\displaystyle \rho{\left(x,y,t \right)} u{\left(x,y,t \right)} \frac{\partial}{\partial x} v{\left(x,y,t \right)} + \rho{\left(x,y,t \right)} v{\left(x,y,t \right)} \frac{\partial}{\partial x} u{\left(x,y,t \right)} + 2 \rho{\left(x,y,t \right)} v{\left(x,y,t \right)} \frac{\partial}{\partial y} v{\left(x,y,t \right)} + \rho{\left(x,y,t \right)} \frac{\partial}{\partial t} v{\left(x,y,t \right)} + u{\left(x,y,t \right)} v{\left(x,y,t \right)} \frac{\partial}{\partial x} \rho{\left(x,y,t \right)} + v^{2}{\left(x,y,t \right)} \frac{\partial}{\partial y} \rho{\left(x,y,t \right)} + v{\left(x,y,t \right)} \frac{\partial}{\partial t} \rho{\left(x,y,t \right)} = \mu \left(\frac{\partial^{2}}{\partial x^{2}} v{\left(x,y,t \right)} + \frac{\partial^{2}}{\partial y\partial x} u{\left(x,y,t \right)}\right) - \frac{2 \mu \left(\frac{\partial^{2}}{\partial y^{2}} v{\left(x,y,t \right)} + \frac{\partial^{2}}{\partial y\partial x} u{\left(x,y,t \right)}\right)}{3} + 2 \mu \frac{\partial^{2}}{\partial y^{2}} v{\left(x,y,t \right)} - c_{p} \frac{\partial}{\partial y} \rho{\left(x,y,t \right)} + \operatorname{Fe_{v}}{\left(x,y \right)} \rho{\left(x,y,t \right)}$

Differentiaaliyhtälöistä differenssiyhtälöiksi

Ylläoleva yhtälöt ovat differentiaaliyhtälöitä, eli ne kertovat, mitä äärettömän pienet muutokset ajassa tai paikassa vaikuttavat. Tietokoneeni muistiin ei sovi ääretön määrä isoja eikä pieniä muutoksia eikä tietokoneeni erota äärettömän pientä nollasta.

Tietokoneella simuloitaessa maailma jaetaan pieniksi kuutioiksi ja lasketaan, mitä millekin kuutiolle tapahtuu. Pienin mahdollinen siirtymä on hyppy kuutiosta viereiseen tai $\Delta T$:n suuruinen hyppy ajassa. Yhtälöt pitää muuttaa kuvaamaan tuollaisia hyppyjä.

Wikipediasta löytyi tähänkin ongelmaan allaolevan mukainen valmis ratkaisu.

Esimerkiksi allaolevista ylimmän yhtälön mukaisesti pinnankorkeuden osittaisderivaatta itä-länsi suunnassa voidaan approksimoida laskemalla yhteen pinnankorkeus kahdessa viereisessä kuutiossa ja jakamalla kuutioiden keskipisteiden etäisyydellä.

Yksinkertaisimmat allaolevista kaavoista johdetaan kaiketi jo lukiomatikassa, mutta pahimpia en ehkä olisi itse osannut johtaa — siis saada oikeaa tulosta.

---

Finite differences can be considered in more than one variable. They are analogous to partial derivatives in several variables.

Some partial derivative approximations are:

$$ \begin{aligned}f_{x}(x,y)&\approx {\frac {f(x+h,y)-f(x-h,y)}{2h}}\\f_{y}(x,y)&\approx {\frac {f(x,y+k)-f(x,y-k)}{2k}}\\f_{xx}(x,y)&\approx {\frac {f(x+h,y)-2f(x,y)+f(x-h,y)}{h^{2}}}\\f_{yy}(x,y)&\approx {\frac {f(x,y+k)-2f(x,y)+f(x,y-k)}{k^{2}}}\\f_{xy}(x,y)&\approx {\frac {f(x+h,y+k)-f(x+h,y-k)-f(x-h,y+k)+f(x-h,y-k)}{4hk}}\end{aligned} $$

/https://en.wikipedia.org/wiki/Finite_difference#multivariate_finite_differences/

---

Alla ensin yksi esimerkki diskretoinnista sympyn avulla ja sitten kaikki yo. säännöt sympy-muodossa

kitkavoima = c_fr*u(x,y,t)
display('kitkavoima:', kitkavoima) 
diskretointiSaanto = [(sp.diff(u(x,y,t),x), (u_d(i,j+1)-u_d(i,j-1))/ds/2)]
jatkuva = sp.diff(kitkavoima, x)
diskretoitu = jatkuva.subs(diskretointiSaanto)
print('kitkavoiman derivaatta = diskretoitu')
display(sp.Eq(jatkuva, diskretoitu))

'kitkavoima:'

$\displaystyle c_{fr} u{\left(x,y,t \right)}$

kitkavoiman derivaatta = diskretoitu

$\displaystyle c_{fr} \frac{\partial}{\partial x} u{\left(x,y,t \right)} = \frac{c_{fr} \left(- \operatorname{u_{d}}{\left(i,j - 1 \right)} + \operatorname{u_{d}}{\left(i,j + 1 \right)}\right)}{2 \Delta s}$

discr_subs = [
      (sp.diff(v(x,y,t),x,2), 
       (v_d(i,j+1)-2*v_d(i,j)+v_d(i,j-1))/ds**2), 
      (sp.diff(v(x,y,t),y,2), 
       (v_d(i+1,j)-2*v_d(i,j)+v_d(i-1,j))/ds**2),
      (sp.diff(v(x,y,t),x,y), 
       (v_d(i+1,j+1) - v_d(i+1,j-1) - v_d(i-1,j+1) + v_d(i-1,j-1))/4/ds**2),
      (sp.diff(v(x,y,t),x), 
       (v_d(i,j+1)-v_d(i,j-1))/ds/2), 
      (sp.diff(v(x,y,t),y), 
       (v_d(i+1,j)-v_d(i-1,j))/ds/2), 

      (sp.diff(u(x,y,t),x,2), 
       (u_d(i,j+1)-2*u_d(i,j)+u_d(i,j-1))/ds**2),
      (sp.diff(u(x,y,t),y,2), 
       (u_d(i+1,j)-2*u_d(i,j)+u_d(i-1,j))/ds**2),
      (sp.diff(u(x,y,t),x,y), 
       (u_d(i+1,j+1) - u_d(i+1,j-1) - u_d(i-1,j+1) + u_d(i-1,j-1))/4/ds**2),
      (sp.diff(u(x,y,t),x), 
       (u_d(i,j+1)-u_d(i,j-1))/ds/2),
      (sp.diff(u(x,y,t),y), 
       (u_d(i+1,j)-u_d(i-1,j))/ds/2),
    
      (sp.diff(rho(x,y,t),x,2), 
       (rho_d(i,j+1)-2*rho_d(i,j)+rho_d(i,j-1))/ds**2), 
      (sp.diff(rho(x,y,t),y,2), 
       (rho_d(i+1,j)-2*rho_d(i,j)+rho_d(i-1,j))/ds**2), 
      (sp.diff(rho(x,y,t),x,y), 
       (rho_d(i+1,j+1) - rho_d(i+1,j-1) - rho_d(i-1,j+1) + rho_d(i-1,j-1))/4/ds**2),
      (sp.diff(rho(x,y,t),x), 
       (rho_d(i,j+1)-rho_d(i,j-1))/ds/2), 
      (sp.diff(rho(x,y,t),y), 
       (rho_d(i+1,j)-rho_d(i-1,j))/ds/2),
    
      (sp.diff(paint(x,y,t),x,2), 
       (paint_d(i,j+1)-2*paint_d(i,j)+paint_d(i,j-1))/ds**2), 
      (sp.diff(paint(x,y,t),y,2), 
       (paint_d(i+1,j)-2*paint_d(i,j)+paint_d(i-1,j))/ds**2), 
      (sp.diff(paint(x,y,t),x,y), 
       (paint_d(i+1,j+1) - paint_d(i+1,j-1) - paint_d(i-1,j+1) + paint_d(i-1,j-1))/4/ds**2),
      (sp.diff(paint(x,y,t),x), 
       (paint_d(i,j+1)-paint_d(i,j-1))/ds/2), 
      (sp.diff(paint(x,y,t),y), 
       (paint_d(i+1,j)-paint_d(i-1,j))/ds/2),
    
    (sp.diff(u(x,y,t),t), D_u(i,j)),
    (sp.diff(v(x,y,t),t), D_v(i,j)),
    (sp.diff(rho(x,y,t),t), D_rho(i,j)), 
    (sp.diff(paint(x,y,t),t), D_paint(i,j)),
    
      (rho(x,y,t), rho_d(i,j)), 
    (paint(x,y,t), paint_d(i,j)), 
      (u(x,y,t), u_d(i,j)),
    (v(x,y,t), v_d(i,j)), 
      (Fe_u(x,y), Fe_u(i,j)), 
    (Fe_v(x,y), Fe_v(i,j))
    ]

Koska ylläolevia korvaussääntöjä lienee vaikea lukea, tein esimerkin niiden soveltamisesta. Kirjoitin listallisen derivaattojen lausekkeita ja ...

difs = [
    sp.diff(u(x,y,t),t), sp.diff(u(x,y,t),x), sp.diff(u(x,y,t),y), sp.diff(u(x,y,t),x,2), sp.diff(u(x,y,t),y,2), sp.diff(u(x,y,t),x,y),
    sp.diff(v(x,y,t),t), sp.diff(v(x,y,t),x), sp.diff(v(x,y,t),y), sp.diff(v(x,y,t),x,2), sp.diff(v(x,y,t),y,2), sp.diff(v(x,y,t),x,y),
    sp.diff(rho(x,y,t),t), sp.diff(rho(x,y,t),x), sp.diff(rho(x,y,t),y), sp.diff(rho(x,y,t),x,2), 
    sp.diff(rho(x,y,t),y,2), sp.diff(rho(x,y,t),x,y),
    sp.diff(paint(x,y,t),t), sp.diff(paint(x,y,t),x), sp.diff(paint(x,y,t),y), sp.diff(paint(x,y,t),x,2), 
    sp.diff(paint(x,y,t),y,2), sp.diff(paint(x,y,t),x,y)
            ]

for dif in difs:
    display(sp.Eq(dif, dif.subs(discr_subs)))
    print('--------------------------------------')

$\displaystyle \frac{\partial}{\partial t} u{\left(x,y,t \right)} = \operatorname{D_{u}}{\left(i,j \right)}$

--------------------------------------

$\displaystyle \frac{\partial}{\partial x} u{\left(x,y,t \right)} = \frac{- \operatorname{u_{d}}{\left(i,j - 1 \right)} + \operatorname{u_{d}}{\left(i,j + 1 \right)}}{2 \Delta s}$

--------------------------------------

$\displaystyle \frac{\partial}{\partial y} u{\left(x,y,t \right)} = \frac{- \operatorname{u_{d}}{\left(i - 1,j \right)} + \operatorname{u_{d}}{\left(i + 1,j \right)}}{2 \Delta s}$

--------------------------------------

$\displaystyle \frac{\partial^{2}}{\partial x^{2}} u{\left(x,y,t \right)} = \frac{- 2 \operatorname{u_{d}}{\left(i,j \right)} + \operatorname{u_{d}}{\left(i,j - 1 \right)} + \operatorname{u_{d}}{\left(i,j + 1 \right)}}{\Delta s^{2}}$

--------------------------------------

$\displaystyle \frac{\partial^{2}}{\partial y^{2}} u{\left(x,y,t \right)} = \frac{- 2 \operatorname{u_{d}}{\left(i,j \right)} + \operatorname{u_{d}}{\left(i - 1,j \right)} + \operatorname{u_{d}}{\left(i + 1,j \right)}}{\Delta s^{2}}$

--------------------------------------

$\displaystyle \frac{\partial^{2}}{\partial y\partial x} u{\left(x,y,t \right)} = \frac{\frac{\operatorname{u_{d}}{\left(i - 1,j - 1 \right)}}{4} - \frac{\operatorname{u_{d}}{\left(i - 1,j + 1 \right)}}{4} - \frac{\operatorname{u_{d}}{\left(i + 1,j - 1 \right)}}{4} + \frac{\operatorname{u_{d}}{\left(i + 1,j + 1 \right)}}{4}}{\Delta s^{2}}$

--------------------------------------

$\displaystyle \frac{\partial}{\partial t} v{\left(x,y,t \right)} = \operatorname{D_{v}}{\left(i,j \right)}$

--------------------------------------

$\displaystyle \frac{\partial}{\partial x} v{\left(x,y,t \right)} = \frac{- \operatorname{v_{d}}{\left(i,j - 1 \right)} + \operatorname{v_{d}}{\left(i,j + 1 \right)}}{2 \Delta s}$

--------------------------------------

$\displaystyle \frac{\partial}{\partial y} v{\left(x,y,t \right)} = \frac{- \operatorname{v_{d}}{\left(i - 1,j \right)} + \operatorname{v_{d}}{\left(i + 1,j \right)}}{2 \Delta s}$

--------------------------------------

$\displaystyle \frac{\partial^{2}}{\partial x^{2}} v{\left(x,y,t \right)} = \frac{- 2 \operatorname{v_{d}}{\left(i,j \right)} + \operatorname{v_{d}}{\left(i,j - 1 \right)} + \operatorname{v_{d}}{\left(i,j + 1 \right)}}{\Delta s^{2}}$

--------------------------------------

$\displaystyle \frac{\partial^{2}}{\partial y^{2}} v{\left(x,y,t \right)} = \frac{- 2 \operatorname{v_{d}}{\left(i,j \right)} + \operatorname{v_{d}}{\left(i - 1,j \right)} + \operatorname{v_{d}}{\left(i + 1,j \right)}}{\Delta s^{2}}$

--------------------------------------

$\displaystyle \frac{\partial^{2}}{\partial y\partial x} v{\left(x,y,t \right)} = \frac{\frac{\operatorname{v_{d}}{\left(i - 1,j - 1 \right)}}{4} - \frac{\operatorname{v_{d}}{\left(i - 1,j + 1 \right)}}{4} - \frac{\operatorname{v_{d}}{\left(i + 1,j - 1 \right)}}{4} + \frac{\operatorname{v_{d}}{\left(i + 1,j + 1 \right)}}{4}}{\Delta s^{2}}$

--------------------------------------

$\displaystyle \frac{\partial}{\partial t} \rho{\left(x,y,t \right)} = \operatorname{D_{\rho}}{\left(i,j \right)}$

--------------------------------------

$\displaystyle \frac{\partial}{\partial x} \rho{\left(x,y,t \right)} = \frac{- \rho_{d}{\left(i,j - 1 \right)} + \rho_{d}{\left(i,j + 1 \right)}}{2 \Delta s}$

--------------------------------------

$\displaystyle \frac{\partial}{\partial y} \rho{\left(x,y,t \right)} = \frac{- \rho_{d}{\left(i - 1,j \right)} + \rho_{d}{\left(i + 1,j \right)}}{2 \Delta s}$

--------------------------------------

$\displaystyle \frac{\partial^{2}}{\partial x^{2}} \rho{\left(x,y,t \right)} = \frac{- 2 \rho_{d}{\left(i,j \right)} + \rho_{d}{\left(i,j - 1 \right)} + \rho_{d}{\left(i,j + 1 \right)}}{\Delta s^{2}}$

--------------------------------------

$\displaystyle \frac{\partial^{2}}{\partial y^{2}} \rho{\left(x,y,t \right)} = \frac{- 2 \rho_{d}{\left(i,j \right)} + \rho_{d}{\left(i - 1,j \right)} + \rho_{d}{\left(i + 1,j \right)}}{\Delta s^{2}}$

--------------------------------------

$\displaystyle \frac{\partial^{2}}{\partial y\partial x} \rho{\left(x,y,t \right)} = \frac{\frac{\rho_{d}{\left(i - 1,j - 1 \right)}}{4} - \frac{\rho_{d}{\left(i - 1,j + 1 \right)}}{4} - \frac{\rho_{d}{\left(i + 1,j - 1 \right)}}{4} + \frac{\rho_{d}{\left(i + 1,j + 1 \right)}}{4}}{\Delta s^{2}}$

--------------------------------------

$\displaystyle \frac{\partial}{\partial t} \operatorname{paint}{\left(x,y,t \right)} = \operatorname{D_{paint}}{\left(i,j \right)}$

--------------------------------------

$\displaystyle \frac{\partial}{\partial x} \operatorname{paint}{\left(x,y,t \right)} = \frac{- \operatorname{paint_{d}}{\left(i,j - 1 \right)} + \operatorname{paint_{d}}{\left(i,j + 1 \right)}}{2 \Delta s}$

--------------------------------------

$\displaystyle \frac{\partial}{\partial y} \operatorname{paint}{\left(x,y,t \right)} = \frac{- \operatorname{paint_{d}}{\left(i - 1,j \right)} + \operatorname{paint_{d}}{\left(i + 1,j \right)}}{2 \Delta s}$

--------------------------------------

$\displaystyle \frac{\partial^{2}}{\partial x^{2}} \operatorname{paint}{\left(x,y,t \right)} = \frac{- 2 \operatorname{paint_{d}}{\left(i,j \right)} + \operatorname{paint_{d}}{\left(i,j - 1 \right)} + \operatorname{paint_{d}}{\left(i,j + 1 \right)}}{\Delta s^{2}}$

--------------------------------------

$\displaystyle \frac{\partial^{2}}{\partial y^{2}} \operatorname{paint}{\left(x,y,t \right)} = \frac{- 2 \operatorname{paint_{d}}{\left(i,j \right)} + \operatorname{paint_{d}}{\left(i - 1,j \right)} + \operatorname{paint_{d}}{\left(i + 1,j \right)}}{\Delta s^{2}}$

--------------------------------------

$\displaystyle \frac{\partial^{2}}{\partial y\partial x} \operatorname{paint}{\left(x,y,t \right)} = \frac{\frac{\operatorname{paint_{d}}{\left(i - 1,j - 1 \right)}}{4} - \frac{\operatorname{paint_{d}}{\left(i - 1,j + 1 \right)}}{4} - \frac{\operatorname{paint_{d}}{\left(i + 1,j - 1 \right)}}{4} + \frac{\operatorname{paint_{d}}{\left(i + 1,j + 1 \right)}}{4}}{\Delta s^{2}}$

--------------------------------------

---

---

Järvi kuutioiksi

Allaolevan kuvan on tarkoitus auttaa ymmärtämään, miten sovelsin ylläolevia diskretointisääntöjä.

Pari kappaletta numeerisesta integroinnista. Ellei numeerinen integrointi kiinnosta, tämän kappaleen voi ohittaa.

Ensimmäisellä yrityksellä tein ylläolevan kaltaisen korvauksen myös derivaatoille ajan suhteen. Sain yhtälöt, joista saattoi ratkaista muuttujien arvot hetkellä $t + \Delta t$ kun tiesi muuttujien arvot hetkellä $t$. Menetelmää kutsutaan Eulerin menetelmäksi. Tiesin sen epätarkaksi, mutten pitänyt sitä ongelmana. Enhän ole suunnittelemassa lentokoneita.

Eulerin menetelmä on kuitenkin epätarkka tavalla, jonka takia muuttujien arvot karkasivat äärettömyyksiin, joten jouduin turvautumaan pikkuisen monimutkaisempaan Heunin algoritmiin. Siitä lisää myöhemmin numeerista laskentaa esittelevässä luvussa.

Heunin algoritmi varten tässä ratkaistaan muuttujien aikaderivaatat.

Reunaehdot

Ylläolevat säilymislait pätevät keskellä 'järveä', mutta rannoilla pätee omat ehtonsa. Pienellä kaistaleella rannan tuntumassa vesi ei virtaa — ei rannan suuntaisesti, eikä kohtisuoraan rantaan nähden. (Oikeissa järvissä vesi loiskuaa edes takaisin rannalla, mutta eipä siinäkään taida mitään nettovirtausta olla.)

Eli mitä tehdään, kun ylläolevissa yhtälöissä viereinen ruutu onkin maalla?

Simulointiohjelmassani järvi on iso taulukko. Simuloinnin alussa nopeudet laitetaan nolliksi koko taulukossa, myös 'maalla', eikä 'maalla' olevia arvoja käydä sen jälkeen muuttamassa, joten nopeuksien reunaehto toteutuu itsestään. Sama pätee väriaineen konsetraatioon.

Pieni syrjähyppy: Muuttamalla virtauksen nopeutta rantaviivalla simulaatioon voi lisätä puron, josta vesi laskee järveen.

'Reunoilla' ei ole paine-eroa kerrottiin eräässä artikkelissa. Kuulosti järkevältä, vaan miten nollata paine-ero?

Ohjelmani laskee uusia arvoja vain järvellä oleville ruuduille, muttei erota, ovatko laskennassa käytettävät naapuriruudut maalla vai järvellä. Niinpä kopioin järveltä pinnankorkeuden arvon viereiselle rantaviivalle jokaisella aika-askelella. Jos maalla ja rannassa pinnakorkeus on sama, paine-ero on nolla. Selitän tarkemmin myöhemmin.

Ratkaistaan aikaderivaattojen lausekkeet

Ratkaistaan pinnankorkeuden, itä-länsi ja pohjois-etelä suuntaisten nopeuksien ja väriaineen määrän derivaatat ajan suhteen

# ratkaistavat muuttujat
sys_vars = [D_rho(i,j), D_u(i,j), D_v(i,j), D_paint(i,j)]

# Säilymislait differentiaaliyhtälöryhmänä
eqs = [eq_mass, eq_mom_x, eq_mom_y, eq_paint]

# joista saadaan säilymislait differenssiyhtälöinä
deqs = [eq.subs(discr_subs) for eq in eqs]

for eq in (deqs):
    display(eq)
    print('\n-------------------------\n')

$\displaystyle \operatorname{D_{\rho}}{\left(i,j \right)} + \frac{\left(- \rho_{d}{\left(i,j - 1 \right)} + \rho_{d}{\left(i,j + 1 \right)}\right) \operatorname{u_{d}}{\left(i,j \right)}}{2 \Delta s} + \frac{\left(- \rho_{d}{\left(i - 1,j \right)} + \rho_{d}{\left(i + 1,j \right)}\right) \operatorname{v_{d}}{\left(i,j \right)}}{2 \Delta s} + \frac{\left(- \operatorname{u_{d}}{\left(i,j - 1 \right)} + \operatorname{u_{d}}{\left(i,j + 1 \right)}\right) \rho_{d}{\left(i,j \right)}}{2 \Delta s} + \frac{\left(- \operatorname{v_{d}}{\left(i - 1,j \right)} + \operatorname{v_{d}}{\left(i + 1,j \right)}\right) \rho_{d}{\left(i,j \right)}}{2 \Delta s} = 0$

-------------------------

$\displaystyle \operatorname{D_{\rho}}{\left(i,j \right)} \operatorname{u_{d}}{\left(i,j \right)} + \operatorname{D_{u}}{\left(i,j \right)} \rho_{d}{\left(i,j \right)} + \frac{\left(- \rho_{d}{\left(i,j - 1 \right)} + \rho_{d}{\left(i,j + 1 \right)}\right) \operatorname{u_{d}}^{2}{\left(i,j \right)}}{2 \Delta s} + \frac{\left(- \rho_{d}{\left(i - 1,j \right)} + \rho_{d}{\left(i + 1,j \right)}\right) \operatorname{u_{d}}{\left(i,j \right)} \operatorname{v_{d}}{\left(i,j \right)}}{2 \Delta s} + \frac{\left(- \operatorname{u_{d}}{\left(i,j - 1 \right)} + \operatorname{u_{d}}{\left(i,j + 1 \right)}\right) \rho_{d}{\left(i,j \right)} \operatorname{u_{d}}{\left(i,j \right)}}{\Delta s} + \frac{\left(- \operatorname{u_{d}}{\left(i - 1,j \right)} + \operatorname{u_{d}}{\left(i + 1,j \right)}\right) \rho_{d}{\left(i,j \right)} \operatorname{v_{d}}{\left(i,j \right)}}{2 \Delta s} + \frac{\left(- \operatorname{v_{d}}{\left(i - 1,j \right)} + \operatorname{v_{d}}{\left(i + 1,j \right)}\right) \rho_{d}{\left(i,j \right)} \operatorname{u_{d}}{\left(i,j \right)}}{2 \Delta s} = - \frac{2 \mu \left(\frac{- 2 \operatorname{u_{d}}{\left(i,j \right)} + \operatorname{u_{d}}{\left(i,j - 1 \right)} + \operatorname{u_{d}}{\left(i,j + 1 \right)}}{\Delta s^{2}} + \frac{\frac{\operatorname{v_{d}}{\left(i - 1,j - 1 \right)}}{4} - \frac{\operatorname{v_{d}}{\left(i - 1,j + 1 \right)}}{4} - \frac{\operatorname{v_{d}}{\left(i + 1,j - 1 \right)}}{4} + \frac{\operatorname{v_{d}}{\left(i + 1,j + 1 \right)}}{4}}{\Delta s^{2}}\right)}{3} + \mu \left(\frac{- 2 \operatorname{u_{d}}{\left(i,j \right)} + \operatorname{u_{d}}{\left(i - 1,j \right)} + \operatorname{u_{d}}{\left(i + 1,j \right)}}{\Delta s^{2}} + \frac{\frac{\operatorname{v_{d}}{\left(i - 1,j - 1 \right)}}{4} - \frac{\operatorname{v_{d}}{\left(i - 1,j + 1 \right)}}{4} - \frac{\operatorname{v_{d}}{\left(i + 1,j - 1 \right)}}{4} + \frac{\operatorname{v_{d}}{\left(i + 1,j + 1 \right)}}{4}}{\Delta s^{2}}\right) + \operatorname{Fe_{u}}{\left(i,j \right)} \rho_{d}{\left(i,j \right)} - \frac{c_{p} \left(- \rho_{d}{\left(i,j - 1 \right)} + \rho_{d}{\left(i,j + 1 \right)}\right)}{2 \Delta s} + \frac{2 \mu \left(- 2 \operatorname{u_{d}}{\left(i,j \right)} + \operatorname{u_{d}}{\left(i,j - 1 \right)} + \operatorname{u_{d}}{\left(i,j + 1 \right)}\right)}{\Delta s^{2}}$

-------------------------

$\displaystyle \operatorname{D_{\rho}}{\left(i,j \right)} \operatorname{v_{d}}{\left(i,j \right)} + \operatorname{D_{v}}{\left(i,j \right)} \rho_{d}{\left(i,j \right)} + \frac{\left(- \rho_{d}{\left(i,j - 1 \right)} + \rho_{d}{\left(i,j + 1 \right)}\right) \operatorname{u_{d}}{\left(i,j \right)} \operatorname{v_{d}}{\left(i,j \right)}}{2 \Delta s} + \frac{\left(- \rho_{d}{\left(i - 1,j \right)} + \rho_{d}{\left(i + 1,j \right)}\right) \operatorname{v_{d}}^{2}{\left(i,j \right)}}{2 \Delta s} + \frac{\left(- \operatorname{u_{d}}{\left(i,j - 1 \right)} + \operatorname{u_{d}}{\left(i,j + 1 \right)}\right) \rho_{d}{\left(i,j \right)} \operatorname{v_{d}}{\left(i,j \right)}}{2 \Delta s} + \frac{\left(- \operatorname{v_{d}}{\left(i,j - 1 \right)} + \operatorname{v_{d}}{\left(i,j + 1 \right)}\right) \rho_{d}{\left(i,j \right)} \operatorname{u_{d}}{\left(i,j \right)}}{2 \Delta s} + \frac{\left(- \operatorname{v_{d}}{\left(i - 1,j \right)} + \operatorname{v_{d}}{\left(i + 1,j \right)}\right) \rho_{d}{\left(i,j \right)} \operatorname{v_{d}}{\left(i,j \right)}}{\Delta s} = \mu \left(\frac{- 2 \operatorname{v_{d}}{\left(i,j \right)} + \operatorname{v_{d}}{\left(i,j - 1 \right)} + \operatorname{v_{d}}{\left(i,j + 1 \right)}}{\Delta s^{2}} + \frac{\frac{\operatorname{u_{d}}{\left(i - 1,j - 1 \right)}}{4} - \frac{\operatorname{u_{d}}{\left(i - 1,j + 1 \right)}}{4} - \frac{\operatorname{u_{d}}{\left(i + 1,j - 1 \right)}}{4} + \frac{\operatorname{u_{d}}{\left(i + 1,j + 1 \right)}}{4}}{\Delta s^{2}}\right) - \frac{2 \mu \left(\frac{- 2 \operatorname{v_{d}}{\left(i,j \right)} + \operatorname{v_{d}}{\left(i - 1,j \right)} + \operatorname{v_{d}}{\left(i + 1,j \right)}}{\Delta s^{2}} + \frac{\frac{\operatorname{u_{d}}{\left(i - 1,j - 1 \right)}}{4} - \frac{\operatorname{u_{d}}{\left(i - 1,j + 1 \right)}}{4} - \frac{\operatorname{u_{d}}{\left(i + 1,j - 1 \right)}}{4} + \frac{\operatorname{u_{d}}{\left(i + 1,j + 1 \right)}}{4}}{\Delta s^{2}}\right)}{3} + \operatorname{Fe_{v}}{\left(i,j \right)} \rho_{d}{\left(i,j \right)} - \frac{c_{p} \left(- \rho_{d}{\left(i - 1,j \right)} + \rho_{d}{\left(i + 1,j \right)}\right)}{2 \Delta s} + \frac{2 \mu \left(- 2 \operatorname{v_{d}}{\left(i,j \right)} + \operatorname{v_{d}}{\left(i - 1,j \right)} + \operatorname{v_{d}}{\left(i + 1,j \right)}\right)}{\Delta s^{2}}$

-------------------------

$\displaystyle - D_{pnt} \left(\frac{- 2 \operatorname{paint_{d}}{\left(i,j \right)} + \operatorname{paint_{d}}{\left(i,j - 1 \right)} + \operatorname{paint_{d}}{\left(i,j + 1 \right)}}{\Delta s^{2}} + \frac{- 2 \operatorname{paint_{d}}{\left(i,j \right)} + \operatorname{paint_{d}}{\left(i - 1,j \right)} + \operatorname{paint_{d}}{\left(i + 1,j \right)}}{\Delta s^{2}}\right) + \operatorname{D_{paint}}{\left(i,j \right)} + \frac{\left(- \operatorname{paint_{d}}{\left(i,j - 1 \right)} + \operatorname{paint_{d}}{\left(i,j + 1 \right)}\right) \operatorname{u_{d}}{\left(i,j \right)}}{2 \Delta s} + \frac{\left(- \operatorname{paint_{d}}{\left(i - 1,j \right)} + \operatorname{paint_{d}}{\left(i + 1,j \right)}\right) \operatorname{v_{d}}{\left(i,j \right)}}{2 \Delta s} + \frac{\left(- \operatorname{u_{d}}{\left(i,j - 1 \right)} + \operatorname{u_{d}}{\left(i,j + 1 \right)}\right) \operatorname{paint_{d}}{\left(i,j \right)}}{2 \Delta s} + \frac{\left(- \operatorname{v_{d}}{\left(i - 1,j \right)} + \operatorname{v_{d}}{\left(i + 1,j \right)}\right) \operatorname{paint_{d}}{\left(i,j \right)}}{2 \Delta s} = 0$

-------------------------

Ratkaistaan ylläolevasta derivaatat ajan suhteen (Ei onnistuisi ainakaan minulta käsityönä)

sys_eq = sp.solve(deqs, sys_vars)

for eq in (sys_eq):
    display(sp.Eq(eq, sys_eq[eq]))
    print('\n-------------------------\n')

$\displaystyle \operatorname{D_{\rho}}{\left(i,j \right)} = \frac{\rho_{d}{\left(i,j \right)} \operatorname{u_{d}}{\left(i,j - 1 \right)}}{2 \Delta s} - \frac{\rho_{d}{\left(i,j \right)} \operatorname{u_{d}}{\left(i,j + 1 \right)}}{2 \Delta s} + \frac{\rho_{d}{\left(i,j \right)} \operatorname{v_{d}}{\left(i - 1,j \right)}}{2 \Delta s} - \frac{\rho_{d}{\left(i,j \right)} \operatorname{v_{d}}{\left(i + 1,j \right)}}{2 \Delta s} + \frac{\rho_{d}{\left(i,j - 1 \right)} \operatorname{u_{d}}{\left(i,j \right)}}{2 \Delta s} - \frac{\rho_{d}{\left(i,j + 1 \right)} \operatorname{u_{d}}{\left(i,j \right)}}{2 \Delta s} + \frac{\rho_{d}{\left(i - 1,j \right)} \operatorname{v_{d}}{\left(i,j \right)}}{2 \Delta s} - \frac{\rho_{d}{\left(i + 1,j \right)} \operatorname{v_{d}}{\left(i,j \right)}}{2 \Delta s}$

-------------------------

$\displaystyle \operatorname{D_{u}}{\left(i,j \right)} = \operatorname{Fe_{u}}{\left(i,j \right)} + \frac{c_{p} \rho_{d}{\left(i,j - 1 \right)}}{2 \Delta s \rho_{d}{\left(i,j \right)}} - \frac{c_{p} \rho_{d}{\left(i,j + 1 \right)}}{2 \Delta s \rho_{d}{\left(i,j \right)}} + \frac{\operatorname{u_{d}}{\left(i,j \right)} \operatorname{u_{d}}{\left(i,j - 1 \right)}}{2 \Delta s} - \frac{\operatorname{u_{d}}{\left(i,j \right)} \operatorname{u_{d}}{\left(i,j + 1 \right)}}{2 \Delta s} + \frac{\operatorname{u_{d}}{\left(i - 1,j \right)} \operatorname{v_{d}}{\left(i,j \right)}}{2 \Delta s} - \frac{\operatorname{u_{d}}{\left(i + 1,j \right)} \operatorname{v_{d}}{\left(i,j \right)}}{2 \Delta s} - \frac{14 \mu \operatorname{u_{d}}{\left(i,j \right)}}{3 \Delta s^{2} \rho_{d}{\left(i,j \right)}} + \frac{4 \mu \operatorname{u_{d}}{\left(i,j - 1 \right)}}{3 \Delta s^{2} \rho_{d}{\left(i,j \right)}} + \frac{4 \mu \operatorname{u_{d}}{\left(i,j + 1 \right)}}{3 \Delta s^{2} \rho_{d}{\left(i,j \right)}} + \frac{\mu \operatorname{u_{d}}{\left(i - 1,j \right)}}{\Delta s^{2} \rho_{d}{\left(i,j \right)}} + \frac{\mu \operatorname{u_{d}}{\left(i + 1,j \right)}}{\Delta s^{2} \rho_{d}{\left(i,j \right)}} + \frac{\mu \operatorname{v_{d}}{\left(i - 1,j - 1 \right)}}{12 \Delta s^{2} \rho_{d}{\left(i,j \right)}} - \frac{\mu \operatorname{v_{d}}{\left(i - 1,j + 1 \right)}}{12 \Delta s^{2} \rho_{d}{\left(i,j \right)}} - \frac{\mu \operatorname{v_{d}}{\left(i + 1,j - 1 \right)}}{12 \Delta s^{2} \rho_{d}{\left(i,j \right)}} + \frac{\mu \operatorname{v_{d}}{\left(i + 1,j + 1 \right)}}{12 \Delta s^{2} \rho_{d}{\left(i,j \right)}}$

-------------------------

$\displaystyle \operatorname{D_{v}}{\left(i,j \right)} = \operatorname{Fe_{v}}{\left(i,j \right)} + \frac{c_{p} \rho_{d}{\left(i - 1,j \right)}}{2 \Delta s \rho_{d}{\left(i,j \right)}} - \frac{c_{p} \rho_{d}{\left(i + 1,j \right)}}{2 \Delta s \rho_{d}{\left(i,j \right)}} + \frac{\operatorname{u_{d}}{\left(i,j \right)} \operatorname{v_{d}}{\left(i,j - 1 \right)}}{2 \Delta s} - \frac{\operatorname{u_{d}}{\left(i,j \right)} \operatorname{v_{d}}{\left(i,j + 1 \right)}}{2 \Delta s} + \frac{\operatorname{v_{d}}{\left(i,j \right)} \operatorname{v_{d}}{\left(i - 1,j \right)}}{2 \Delta s} - \frac{\operatorname{v_{d}}{\left(i,j \right)} \operatorname{v_{d}}{\left(i + 1,j \right)}}{2 \Delta s} + \frac{\mu \operatorname{u_{d}}{\left(i - 1,j - 1 \right)}}{12 \Delta s^{2} \rho_{d}{\left(i,j \right)}} - \frac{\mu \operatorname{u_{d}}{\left(i - 1,j + 1 \right)}}{12 \Delta s^{2} \rho_{d}{\left(i,j \right)}} - \frac{\mu \operatorname{u_{d}}{\left(i + 1,j - 1 \right)}}{12 \Delta s^{2} \rho_{d}{\left(i,j \right)}} + \frac{\mu \operatorname{u_{d}}{\left(i + 1,j + 1 \right)}}{12 \Delta s^{2} \rho_{d}{\left(i,j \right)}} - \frac{14 \mu \operatorname{v_{d}}{\left(i,j \right)}}{3 \Delta s^{2} \rho_{d}{\left(i,j \right)}} + \frac{\mu \operatorname{v_{d}}{\left(i,j - 1 \right)}}{\Delta s^{2} \rho_{d}{\left(i,j \right)}} + \frac{\mu \operatorname{v_{d}}{\left(i,j + 1 \right)}}{\Delta s^{2} \rho_{d}{\left(i,j \right)}} + \frac{4 \mu \operatorname{v_{d}}{\left(i - 1,j \right)}}{3 \Delta s^{2} \rho_{d}{\left(i,j \right)}} + \frac{4 \mu \operatorname{v_{d}}{\left(i + 1,j \right)}}{3 \Delta s^{2} \rho_{d}{\left(i,j \right)}}$

-------------------------

$\displaystyle \operatorname{D_{paint}}{\left(i,j \right)} = - \frac{4 D_{pnt} \operatorname{paint_{d}}{\left(i,j \right)}}{\Delta s^{2}} + \frac{D_{pnt} \operatorname{paint_{d}}{\left(i,j - 1 \right)}}{\Delta s^{2}} + \frac{D_{pnt} \operatorname{paint_{d}}{\left(i,j + 1 \right)}}{\Delta s^{2}} + \frac{D_{pnt} \operatorname{paint_{d}}{\left(i - 1,j \right)}}{\Delta s^{2}} + \frac{D_{pnt} \operatorname{paint_{d}}{\left(i + 1,j \right)}}{\Delta s^{2}} + \frac{\operatorname{paint_{d}}{\left(i,j \right)} \operatorname{u_{d}}{\left(i,j - 1 \right)}}{2 \Delta s} - \frac{\operatorname{paint_{d}}{\left(i,j \right)} \operatorname{u_{d}}{\left(i,j + 1 \right)}}{2 \Delta s} + \frac{\operatorname{paint_{d}}{\left(i,j \right)} \operatorname{v_{d}}{\left(i - 1,j \right)}}{2 \Delta s} - \frac{\operatorname{paint_{d}}{\left(i,j \right)} \operatorname{v_{d}}{\left(i + 1,j \right)}}{2 \Delta s} + \frac{\operatorname{paint_{d}}{\left(i,j - 1 \right)} \operatorname{u_{d}}{\left(i,j \right)}}{2 \Delta s} - \frac{\operatorname{paint_{d}}{\left(i,j + 1 \right)} \operatorname{u_{d}}{\left(i,j \right)}}{2 \Delta s} + \frac{\operatorname{paint_{d}}{\left(i - 1,j \right)} \operatorname{v_{d}}{\left(i,j \right)}}{2 \Delta s} - \frac{\operatorname{paint_{d}}{\left(i + 1,j \right)} \operatorname{v_{d}}{\left(i,j \right)}}{2 \Delta s}$

-------------------------

Koodin generointia

Muutetaan symboliset lausekkeet python-koodiksi ja kirjoitetaan python-funktio, joka laskee pinnankorkeuden ja nopeuksien derivaattojen arvot kullekin järveäni edustavalla taulukon ruudulle.

Tähän en keksinyt mitään eleganttia ratkaisua, joten tuskin ketään kiinnostaa, miten tämän tein.

# Muunnetaan yhtälön oikean puoli sympy-lausekkeesta merkkijonoksi
def conv_rhs(rhs0):

    # Parametrien arvoja joutuu muuttelemaan, kun yrittää saada halutun taiteellisen vaikutelman.
    # Numeeriset arvot parametreille ovat moduulissa par.
    # Jos parametrille antaa arvon tässä vaiheessa, sen muuttaminen vaatii koodin generoinnin uudestaan,
    # mutta saattaa nopeuttaa numeerista laskentaa.
    # Päätin, että parametrin ds arvoa en muuta, joten annan sille numeerisen arvon tässä vaiheessa.
    para_subs = [(ds, par.ds)]
    # sympy-lauseke merkkijonoksi
    rhs = str(sp.N(sp.simplify(rhs0.subs(para_subs))))

    # 
    code_subs = [
        (r'\\mu', 'par.mu'), ('c_p', 'par.c_p'), ('c_fr', 'par.c_fr'), ('D_pnt', 'par.D_pnt'), 
        ('rho_0', 'par.rho_0'), ('roo', 'rho_d'),
        (' \+ ', '+'), (' - ', '-'), 
        ('\(i', '[i'), ('j\)', 'j]'), ('j\+1\)', 'j+1]'), ('j-1\)', 'j-1]')
            ]
    for sub in code_subs:
        rhs = re.sub(sub[0], sub[1], rhs)
    return rhs

def conv_lhs(lhs0):
    code_subs = [('D_rho(i, j)', 'D_rho[i, j]'), 
                 ('D_u(i, j)', 'D_u[i, j]'), 
                 ('D_v(i, j)', 'D_v[i, j]'), 
                 ('D_paint(i, j)', 'D_paint[i, j]')]
    lhs = str(lhs0)
    for sub in code_subs:
        if lhs == sub[0]:
            return sub[1]
    return lhs    

Tehdään ylläesitellyt korvaukset systeemiyhtälöihin.

Onnistuu pythonilla aika näppärästi. Toisaalta pythonia tuntemattoman lienee vaikea ymmärtää, mitä alla tapahtuu.

(sympy.Solve antaa ratkaisun yhtälöryhmälle dictionaryna, jonka kelaan lävitse)

py_eqs = [(conv_lhs(eq), conv_rhs(sys_eq[eq])) for eq in sys_eq]

def gen_syseq():
    sis = "    "
    sis2 = sis + sis
    sis3 = sis2 + sis
    sis4 = sis2 + sis2
    sis5 = sis4 + ' '
    sis6 = sis5 + ' '

    with open('dfdt_virtaukset.py', 'w') as fil:
        fil.write("# -*- coding: utf-8 -*-\n\n")
        fil.write("import params as par\n")
        fil.write("import saaret as ly\n\n\n")
        fil.write("def dfdt(kartta, rho_d, u_d, v_d, paint_d, Fe_u, Fe_v, D_rho, D_u, D_v, D_paint):\n")

        fil.write(sis + "for i in range(1, ly.DIM_i-1):\n")
        fil.write(sis2 + "for j in range(1, ly.DIM_j-1):\n")

        fil.write(sis3 + "if kartta[i, j] == 'jarvi':\n")
        for eq in py_eqs:
            fil.write(sis4 + eq[0]+" = \\\n")
            fil.write(fill(eq[1], initial_indent=sis5 + "(", subsequent_indent=sis6, break_long_words=False))
            fil.write(")\n\n")

gen_syseq()

with open('dfdt_virtaukset.py', 'r') as code:
    print(code.read())

# -*- coding: utf-8 -*-

import params as par
import saaret as ly


def dfdt(kartta, rho_d, u_d, v_d, paint_d, Fe_u, Fe_v, D_rho, D_u, D_v, D_paint):
    for i in range(1, ly.DIM_i-1):
        for j in range(1, ly.DIM_j-1):
            if kartta[i, j] == 'jarvi':
                D_rho[i, j] = \
                 (0.5*rho_d[i, j]*u_d[i, j-1]-0.5*rho_d[i, j]*u_d[i,
                  j+1]+0.5*rho_d[i, j]*v_d[i-1, j]-0.5*rho_d[i,
                  j]*v_d[i+1, j]+0.5*rho_d[i, j-1]*u_d[i,
                  j]-0.5*rho_d[i, j+1]*u_d[i, j]+0.5*rho_d[i-1,
                  j]*v_d[i, j]-0.5*rho_d[i+1, j]*v_d[i, j])

                D_u[i, j] = \
                 ((-4.66666666666667*par.mu*u_d[i,
                  j]+1.33333333333333*par.mu*u_d[i,
                  j-1]+1.33333333333333*par.mu*u_d[i,
                  j+1]+1.0*par.mu*u_d[i-1, j]+1.0*par.mu*u_d[i+1,
                  j]+0.0833333333333333*par.mu*v_d[i-1,
                  j-1]-0.0833333333333333*par.mu*v_d[i-1,
                  j+1]-0.0833333333333333*par.mu*v_d[i+1,
                  j-1]+0.0833333333333333*par.mu*v_d[i+1,
                  j+1]+0.5*par.c_p*rho_d[i, j-1]-0.5*par.c_p*rho_d[i,
                  j+1]+(Fe_u[i, j]+0.5*u_d[i, j]*u_d[i,
                  j-1]-0.5*u_d[i, j]*u_d[i, j+1]+0.5*u_d[i-1,
                  j]*v_d[i, j]-0.5*u_d[i+1, j]*v_d[i, j])*rho_d[i,
                  j])/rho_d[i, j])

                D_v[i, j] = \
                 ((0.0833333333333333*par.mu*u_d[i-1,
                  j-1]-0.0833333333333333*par.mu*u_d[i-1,
                  j+1]-0.0833333333333333*par.mu*u_d[i+1,
                  j-1]+0.0833333333333333*par.mu*u_d[i+1,
                  j+1]-4.66666666666667*par.mu*v_d[i,
                  j]+1.0*par.mu*v_d[i, j-1]+1.0*par.mu*v_d[i,
                  j+1]+1.33333333333333*par.mu*v_d[i-1,
                  j]+1.33333333333333*par.mu*v_d[i+1,
                  j]+0.5*par.c_p*rho_d[i-1, j]-0.5*par.c_p*rho_d[i+1,
                  j]+(Fe_v[i, j]+0.5*u_d[i, j]*v_d[i, j-1]-0.5*u_d[i,
                  j]*v_d[i, j+1]+0.5*v_d[i, j]*v_d[i-1, j]-0.5*v_d[i,
                  j]*v_d[i+1, j])*rho_d[i, j])/rho_d[i, j])

                D_paint[i, j] = \
                 (-4.0*par.D_pnt*paint_d[i,
                  j]+1.0*par.D_pnt*paint_d[i,
                  j-1]+1.0*par.D_pnt*paint_d[i,
                  j+1]+1.0*par.D_pnt*paint_d[i-1,
                  j]+1.0*par.D_pnt*paint_d[i+1, j]+0.5*paint_d[i,
                  j]*u_d[i, j-1]-0.5*paint_d[i, j]*u_d[i,
                  j+1]+0.5*paint_d[i, j]*v_d[i-1, j]-0.5*paint_d[i,
                  j]*v_d[i+1, j]+0.5*paint_d[i, j-1]*u_d[i,
                  j]-0.5*paint_d[i, j+1]*u_d[i, j]+0.5*paint_d[i-1,
                  j]*v_d[i, j]-0.5*paint_d[i+1, j]*v_d[i, j])

Simuloinnin kompleksisuus

Videon tekemistä varten joudun simuloimaan ehkä 10000 aika-askelta. Myöhemmin esiteltävä Heunin algoritmi edellyttää, että ylläolevaa funktiolla lasketaan 'järven' jokaisen ruudun pinnankorkeudelle, virtausnopeudelle ja maalin määrälle uusi arvo kaksi kertaa jokaisella aika-askeleella. Taulukossa on 240x200 ruutua, joista enintään puolet on 'maalla'. Uudet arvot lasketaan siis noin N kertaa.

ruudut = 240*200
heun = 2
jarvella = 0.5
aika_askelet = 10000
N = heun*ruudut*aika_askelet*jarvella
print('Uudet arvot lasketaan lasketaan noin {:,} kertaa'.format(N))

Uudet arvot lasketaan lasketaan noin 480,000,000.0 kertaa

Herää kysymys, eikö funktiota dfdt(...) voisi formuloida vektorioperaatioiksi, jotka voisi toteuttaa numpyn avulla. Laskenta nopeutuisi ties kuinka paljon. 10-kertaiseksi? 100-kertaiseksi?

Kysymys herää, mutta onneksi ei sentään herää into koodata. Olisi monen päivän työ päästä selville edes siitä, onnistuuko mokoma.

print(datetime.datetime.now())
print('valmis')

2021-01-14 17:54:27.619693
valmis