Folyadk szimulci Folyadkok Folykony anyagok Fstszer jelensgek Felhk

Folyadék szimuláció

Folyadékok �Folyékony anyagok �Füstszerű jelenségek �Felhők �Festékek

Folyadék állapota �Sebesség (vektormező) x = (x, y) pozíció u = (u, v) sebesség T idő u(x, t) = (u(x, t), v(x, t)) Kocka rács

Navier-Stokes egyenletek (1822) �Claude Navier és George Gabriel Stokes Folyékony anyagok mozgása, áramlása �Alap feltevések Az anyagban fellépő feszültség két összetevője ▪ A sebesség gradiensével arányos diffúzió ▪ Nyomás összetevő

Navier-Stokes egyenletek (1822) �Számos fizikai jelenség leírására alkalmas Időjárás Folyadékok áramlása nem kör keresztmetszetű csatornákban Repülőgépek szárnya körül fellépő áramlás Szilárd testek folyékony anyagokon keresztüli mozgása (pl. a csillagok galaxisokon belül leírt mozgása) Összekapcsolható a Maxwell egyenletekkel (magnetohidrodinamika)

Navier-Stokes egyenletek (1822) �Tisztán elméleti értelemben is fontos Nincs bizonyítva a három dimenziós érvényesség A „létezési és simasági” probléma annyira fontos, hogy a Clay Mathematics Institute az évezred hét legfontosabb matematikai problémái között tartja számon. ▪ A megoldásra egymillió dolláros díjat tűztek ki

Navier-Stokes egyenletek (1822) sűrűség viszkozitás Külső erők Összenyomhatatlan, homogén folyadékok

Az egyenlet tagjai I. �Advekció �Előre haladás, szállítás �Bármilyen mennyiséget �Saját vektormezőt is

Az egyenlet tagjai II. �Nyomás �Az erő nem hirtelen áramlik végig a folyadékon �A molekulák ütköznek, nyomás keletkezik �Gyorsulást (sebességváltozást) eredményez

Az egyenlet tagjai III. �Diffúzió �A különböző folyadékok, különbözőképpen mozognak: vannak sűrűbbek és vannak folyékonyabbak �Viszkozitás: mennyire ellenálló a folyadék az áramlásra �Ez az ellenállás sebesség diffúziót okoz

Az egyenlet tagjai IV. �Külső erők �Lehetnek lokálisak vagy globálisak (pl gravitáció)

Operátorok Operátor Gradiens Divergencia Laplace Definíció Véges differencia alak

Az egyenletek megoldása � 3 egyenlet: u, v, p �Analitikus megoldás ritkán, és csak egyszerű esetekben található �Numerikus módszerek, inkrementális megoldás �Ha animációt szeretnénk, az idő inkrementálás még jól is jön �A problémát kisebb lépésekre bontjuk (Stam, J. 1999. "Stable Fluids. " In Proceedings of SIGGRAPH 1999)

Helmholtz-Hodge dekompozíció (projekciós lépés) �(Bármely vektor felbontható bázisvektorok súlyozott összegére) �Bármely vektormező felbontható vektormezők összegére :

Helmholtz-Hodge dekompozíció (projekciós lépés) W

Hogyan számítsuk ki a nyomást? Poisson egyenlet

Mit is kell tenni egy szimulációs lépésben? Projekció Külső erők Diffúzió Advekció

Advekció �Euler módszer, előrelépés: �Nem stabil (és shaderből nehezen végrehajtható) �A megoldás a visszalépés:

Advekció

Diffúzió �Explicit megoldás: �Nem stabil ! Implicit megoldás: Poisson egyenlet

Projekció Poisson egyenlet

Poisson egyenlet megoldása �Iteratív megoldás, kiindulunk egy kezdeti állapotból és folyamatosan finomítjuk � �Nálunk alakú egyenlet a Laplace operátor �A legegyszerűbb megoldás a Jacobi iteráció

Jacobi iteráció Diffúzió Nyomás x sebesség (u) nyomás(p) b sebesség(u) sebesség divergenciája -1. 0 0. 25

Határfeltételek �Véges tartományon számítunk, kellenek határfeltételek �Ha az anyagot a szimulált tartományba zárjuk (falakkal vesszük körül) a sebességre és a nyomásra a feltételek: Sebesség: a határokon a sebesség nulla (no-slip feltétel) Nyomás: a határokon a nyomás változása nulla (Neumann feltétel)

Kiegészítés: örvénylés �A szimuláció és a diszkretizálás numerikus hibája elmossa a mozgás bizonyos részleteit, a finom örvényeket �Ezeket csaljuk vissza:

Implementáció �A mennyiségeket 2 D tömbökben tároljuk �Mivel a számítások során szomszédossági információk kellenek, néhány mennyiséget dupla bufferben kell tárolni (PING-PONG) �A tömbök frissítését az Open. CL kernelek végzik �Az egyes számítási lépésekhez külön kernelek szükségesek �A megjelenítés egyszerű képernyőre rajzolás �Kernel függvények. . . (folyt)

Advekció __kernel void advection(const int grid. Resolution, __global float 2* input. Velocity. Buffer, __global float 2* output. Velocity. Buffer) { int 2 id = (int 2)(get_global_id(0), get_global_id(1)); if(id. x > 0 && id. x < grid. Resolution - 1 && id. y > 0 && id. y < grid. Resolution – 1){ float 2 velocity = input. Velocity. Buffer[id. x + id. y * grid. Resolution]; float 2 p = (float 2)((float)id. x - dt * velocity. x, (float)id. y - dt * velocity. y); output. Velocity. Buffer[id. x + id. y * grid. Resolution] = get. Bil(p, grid. Resolution, input. Velocity. Buffer); } else{ //határfeltételek if(id. x == 0) output. Velocity. Buffer[id. x + id. y * grid. Resolution] = - input. Velocity. Buffer[id. x + 1 + id. y * grid. Resolution]; . . . } }

Divergencia __kernel void divergence(const int grid. Resolution, __global float 2* velocity. Buffer, __global float* divergence. Buffer) { int 2 id = (int 2)(get_global_id(0), get_global_id(1)); if(id. x > 0 && id. x < grid. Resolution - 1 && id. y > 0 && id. y < grid. Resolution - 1){ float 2 v. L v. R v. B v. T = = velocity. Buffer[id. x + + + 1 + id. y * (id. y - 1) (id. y + 1) grid. Resolution]; * grid. Resolution]; divergence. Buffer[id. x + id. y * grid. Resolution] = 0. 5 f * ((v. R. x - v. L. x) + (v. T. y - v. B. y)); } else{ divergence. Buffer[id. x + id. y * grid. Resolution] = 0. 0 f; } }

Nyomás számítása, Jacobi iteráció __kernelvoid pressure. Jacobi(const int grid. Resolution, __global float* input. Pressure. Buffer, __global float* output. Pressure. Buffer, __global float* divergence. Buffer) { int 2 id = (int 2)(get_global_id(0), get_global_id(1)); if(id. x > 0 && id. x < grid. Resolution - 1 && id. y > 0 && id. y < grid. Resolution - 1){ float float alpha = -1. 0 f; beta = 0. 25 f; v. L = input. Pressure. Buffer[id. x - 1 + id. y * grid. Resolution]; v. R = input. Pressure. Buffer[id. x + 1 + id. y * grid. Resolution]; v. B = input. Pressure. Buffer[id. x + (id. y - 1) * grid. Resolution]; v. T = input. Pressure. Buffer[id. x + (id. y + 1) * grid. Resolution]; divergence = divergence. Buffer[id. x + id. y * grid. Resolution]; output. Pressure. Buffer[id. x + id. y * grid. Resolution] = (v. L + v. R + v. B + v. T + alpha * divergence) * beta; }else{ //határfeltételek if(id. x == 0) output. Pressure. Buffer[id. x + id. y * grid. Resolution] = input. Pressure. Buffer[id. x + 1 + id. y * grid. Resolution]; . . . } }

Projekció __kernel void projection(const int grid. Resolution, __global float 2* input. Velocity. Buffer, __global float* pressure. Buffer, __global float 2* output. Velocity. Buffer) { int 2 id = (int 2)(get_global_id(0), get_global_id(1)); if(id. x > 0 && id. x < grid. Resolution - 1 && id. y > 0 && id. y < grid. Resolution - 1){ float p. L = pressure. Buffer[id. x - 1 + id. y * grid. Resolution]; float p. R = pressure. Buffer[id. x + 1 + id. y * grid. Resolution]; float p. B = pressure. Buffer[id. x + (id. y - 1) * grid. Resolution]; float p. T = pressure. Buffer[id. x + (id. y + 1) * grid. Resolution]; float 2 velocity = input. Velocity. Buffer[id. x + id. y * grid. Resolution]; output. Velocity. Buffer[id. x + id. y * grid. Resolution] = velocity – (float 2)(p. R - p. L, p. T - p. B); } else {//határfeltételek if(id. x == 0) output. Velocity. Buffer[id. x + id. y * grid. Resolution] = -input. Velocity. Buffer[id. x + 1 + id. y * grid. Resolution]; . . . } }

Diffúzió __kernel void diffusion(const int grid. Resolution, __global float 2* input. Velocity. Buffer, __global float 2* output. Velocity. Buffer) { int 2 id = (int 2)(get_global_id(0), get_global_id(1)); float viscousity = 0. 01 f; float alpha = 1. 0 f / (viscousity * dt); float beta = 1. 0 f / (4. 0 f + alpha); if(id. x > 0 && id. x < grid. Resolution - 1 && id. y > 0 && id. y < grid. Resolution - 1){ float 2 float 2 v. L = input. Velocity. Buffer[id. x - 1 + id. y * v. R = input. Velocity. Buffer[id. x + 1 + id. y * v. B = input. Velocity. Buffer[id. x + (id. y - 1) v. T = input. Velocity. Buffer[id. x + (id. y + 1) velocity = input. Velocity. Buffer[id. x + id. y grid. Resolution]; * grid. Resolution]; output. Velocity. Buffer[id. x + id. y * grid. Resolution] = (v. L + v. R + v. B + v. T + alpha * velocity) * beta; } else { output. Velocity. Buffer[id. x + id. y * grid. Resolution] = input. Velocity. Buffer[id. x + id. y * grid. Resolution]; }}

Örvénylés __kernel void vorticity(const int grid. Resolution, __global float 2* velocity. Buffer, __global float* vorticity. Buffer) { int 2 id = (int 2)(get_global_id(0), get_global_id(1)); if(id. x > 0 && id. x < grid. Resolution - 1 && id. y > 0 && id. y < grid. Resolution – 1){ float 2 v. L v. R v. B v. T = = velocity. Buffer[id. x + + + 1 + id. y * (id. y - 1) (id. y + 1) grid. Resolution]; * grid. Resolution]; vorticity. Buffer[id. x + id. y * grid. Resolution] = (v. R. y - v. L. y) - (v. T. x - v. B. x); } else{ vorticity. Buffer[id. x + id. y * grid. Resolution] = 0. 0 f; } }

Sebesség az örvénylésből __kernel void add. Vorticity(const int grid. Resolution, __global float* vorticity. Buffer, __global float 2* velocity. Buffer) { int 2 id = (int 2)(get_global_id(0), get_global_id(1)); const float scale = 0. 2 f; if(id. x > 0 && id. x < grid. Resolution - 1 && id. y > 0 && id. y < grid. Resolution - 1){ float v. L = vorticity. Buffer[id. x - 1 + id. y * grid. Resolution]; float v. R = vorticity. Buffer[id. x + 1 + id. y * grid. Resolution]; float v. B = vorticity. Buffer[id. x + (id. y - 1) * grid. Resolution]; float v. T = vorticity. Buffer[id. x + (id. y + 1) * grid. Resolution]; float 4 grad. V = (float 4)(v. R - v. L, v. T - v. B, 0. 0 f); float 4 z = (float 4)(0. 0 f, 1. 0 f, 0. 0 f); if(dot(grad. V, grad. V)){ float 4 vorticity. Force = scale * cross(grad. V, z); velocity. Buffer[id. x + id. y * grid. Resolution] += vorticity. Force. xy * dt; } } }

Külső erők __kernel void add. Force(const float x, const float y, const float 2 force, const int grid. Resolution, __global float 2* velocity. Buffer, const float 4 density, __global float 4* density. Buffer) { int 2 id = (int 2)(get_global_id(0), get_global_id(1)); float dx = ((float)id. x / (float)grid. Resolution) - x; float dy = ((float)id. y / (float)grid. Resolution) - y; float radius = 0. 001 f; float c = exp( - (dx * dx + dy * dy) / radius ) * dt; velocity. Buffer[id. x + id. y * grid. Resolution] += c * force; density. Buffer[id. x + id. y * grid. Resolution] += c * density; }

Kiegészítési lehetőségek �Felhajtóerő és gravitáció �Termodinamikai szimuláció (felhők) � 3 dimenzióban �Más rács típusok: a vektormezőkre FCC �Tömör testekkel való interakció (voxelizálás, határfeltételek kezelése)

Együttműködés a grafikus API-val �Célja az átjárás megteremtése Open. GL és Direct. X támogatás Megoszthatóak ▪ Általános buffer objektumok (pl. vertex buffer) ▪ Textúrák ▪ Render bufferek A megosztandó objektumokat a grafikus API hozza létre ▪ Open. CL-beli használat előtt zárolni kell Az objektum használata kizárólagos!

Együttműködés a grafikus API-val �Open. GL és Open. CL kontextus megosztás GL_SHARING_EXTENSION Open. GL kontextus információk cl_int cl. Get. GLContext. Info. KHR(const cl_context_properties *props, cl_gl_context_info param_name, size_t param_value_size, void* param_value, size_t* param_value_size_ret) ▪ CL_CURRENT_DEVICE_FOR_GL_CONTEXT_KHR ▪ CL_DEVICES_FOR_GL_CONTEXT_KHR

Együttműködés a grafikus API-val �Open. GL és Open. CL kontextus megosztás Open. CL kontextus létrehozás cl_context cl. Create. Context(const cl_context_properties *props, cl_uint num_devices, const cl_device_id *devices, void (*pfn_notify)(. . . ), void *user_data, cl_int *errcode_ret) ▪ Tulajdonságok: ▪ CL_GL_CONTEXT_KHR: Open. GL kontextus ▪ CL_WGL_HDC_KHR: az Open. GL kontextus HDC-je ▪ CL_CONTEXT_PLATFORM: platform_id

Együttműködés a grafikus API-val �Kontextus megosztása Init. GL(); cl_platform = create. Platform(); cl_device_id = create. Device(platform, CL_DEVICE_TYPE_GPU); cl_context shared. Context = 0; if(Check. Sharing. Support(device_id)){ cl_context_properties props[] = { CL_GL_CONTEXT_KHR, (cl_context_properties)wgl. Get. Current. Context(), CL_WGL_HDC_KHR, (cl_context_properties)wgl. Get. Current. DC(), CL_CONTEXT_PLATFORM, (cl_context_properties)platform, 0 }; shared. Context = cl. Create. Context(props, 1, &device_id, NULL, &err); }

Együttműködés a grafikus API-val �Buffer objektumok megosztása cl_mem cl. Create. From. GLBuffer(cl_context, cl_mem_flags, GLuint bufobj, cl_int* errcode_ret) �Image objektumok megosztása cl_mem cl. Create. From. GLTexture 2 D(cl_context, cl_mem_flags, GLenum texture_target, GLint miplevel, GLuint texture, cl_int* errcode_ret)

Együttműködés a grafikus API-val �Render buffer megosztása cl_mem cl. Create. From. GLRender. Buffer(cl_context, cl_mem_flags, GLuint renderbuffer, cl_int* errcode_ret) �Az Open. CL objektumok tulajdonságai Létrehozáskor aktuális értékek alapján Nem követik az Open. GL objektum változásait! ▪ Amennyiben változik újra meg kell osztani!

Együttműködés a grafikus API-val �Buffer objektum megosztása Open. GL vertex buffer mint Open. CL memória objektum GLuint vbo; gl. Gen. Buffers(1, &vbo); gl. Bind. Buffer(GL_ARRAY_BUFFER, vbo); gl. Buffer. Data(GL_ARRAY_BUFFER, size, 0, GL_DYNAMIC_DRAW); cl_mem vbo. CL; vbo. CL = cl. Create. From. GLBuffer(shared. Context, CL_MEM_WRITE_ONLY, vbo, NULL);

Együttműködés a grafikus API-val �Objektum lefoglalása cl_int cl. Enqueue. Acquire. GLObjects(cl_command_queue command, cl_uint num_objects, const cl_mem* mem_objects, . . . ) �Objektum felszabadítása cl_mem cl. Enqueue. Release. GLObjects(cl_command_queue command, cl_uint num_objects, const cl_mem* mem_objects, . . . ) Minden használat előtt le kell foglalni Használat után fel kell szabadítani

Együttműködés a grafikus API-val �Szinkronizáció Open. GL és Open. CL között Nincs explicit szinkronizáció! ▪ Szüksége lenne mindkét API támogatására Mindkét API oldalán a csővezeték kiürítése ▪ Open. GL: gl. Finish() ▪ Open. CL: cl. Finish() Implementáció függően más megoldás is lehet ▪ gl. Flush() és cl. Enqueue. Barrier()

Együttműködés a grafikus API-val �Buffer objektum használata Open. GL vertex buffer mint Open. CL memória objektum // Open. GL hívások gl. Finish(); cl. Enqueue. Acquire. GLObjects(command, 1, &vbo. CL, 0, NULL); // Kernel paraméterek beállítása és kernel végrehajtás cl. Finish(); cl. Enqueue. Release. GLObjects(commands, 1, &vbo. CL, 0, NULL); // Open. GL hívások