If anyone is interested in playing around with it, I threw it up at JSFiddle here: http://jsfiddle.net/zyAzg/
Excellent demo.
All that matrix code you see in the source here is basically just the minimum prerequisite for doing 3D graphics, it would be the same if you just wanted to get a spinning cube with some camera control.
The real magic happens in the shaders which compute water mesh displacement, normals and shading directly on the GPU (21 simulation passes + 1 rendering pass + 1 initialization pass, using several floating point textures, seems like stateful simulation with ping-ponging, basically using graphics rendering pipeline to do GPGPU without a need for OpenCL / CUDA).
It would be really cool to see a project that would make that easier to do in JS more generally. Basically a similar idea to Parakeet or Numba, but with GLSL as a target.
> It would be really cool to see a project that would make that easier to do in JS more generally. Basically a similar idea to Parakeet or Numba, but with GLSL as a target.
The problem with doing this on WebGL (which is essentially GLES2) is that the internal precision of the pipeline is not required to be a full 32 bit float. There are 10 bit integers and 20 bit floating point numbers and all sorts of funny number formats in the hardware. When most scientific/numeric computation is done in 64 bit double precision, going to less than 32 bits is a problem.
Of course it still could be done in WebGL, the problem is that running it on different hardware would yield different results. If this is acceptable (e.g. image processing where bit-accurate results are not required), then why not.
I noticed that, which I thought was pretty cool. I was playing with it in Chrome on OS X when I realized that it wasn't doing too much to my CPU, even when I cranked the simulation up to max.
GPGPU is more about e.g. doing multi-precision integer arithmetic on the GPU - where the processors are used for something which doesn't fit the 'streaming processor being fed vectors' pattern which the architecture was made for (of which this demo is a very good example) - quite literally 'general purpose'.
actually this is exactly what everyone thought when this hardware was new 10 years ago... and even before that we were thinking about it. look at rendering demos from that time that leveraged the shader - this was the dumbass obvious thing that almost everyone did.
this is not physical simulation - not even close. the navier-stokes equations are still really quite difficult to do anything real-time with unless you make some serious compromises.
i will conceded that this is very loosely based on physical reality. but actually, you can do a lot, /lot/ better with clever hacks.
as an example of one of the much, much better solutions which approximate reality than this - the old school 2d 'water' effect created from feedback and a clever 3x3 convolution filter shaped as a ring with a hole in the middle, for instance is visually much more impressive and computationally simpler - and it actually is a solution to a discretised wave equation. render this old effect as greyscale into a buffer than use it as a height/displacement-map style texture - with this solution implementing real-time ripples and fluid like response to objects entering the liquid comes down to rendering the outline into that buffer and these kinds of waves simply fall out rather than requiring maintenance or trig calls...
http://www.cs.unm.edu/~mskarim/fluid_smoke_gpu.htm
I was doing this in college 10 years ago. Cmon now.
I was at one point perusing a PhD in Computer Graphics and I studied many papers ranging from haptic devices w/ liquid simulations (great paper on making pancakes which change properties as they cook [1]) to research from our lab on volume rendered simulations to prevent temperature shock in introducing cool water (relatively) to prevent a meltdown of a nuclear power plant.
[1]http://www.gmrv.es/Publications/2013/CMOL13/haptic_multistat...)
I once used WolframAlpha to get me some complicated formulas and wrote a custom WA-plaintext-output-to-C translator for that output to use it in my code. Here is the script: https://github.com/albertz/helpers/blob/master/wolframalpha_...
Example input: http://www.wolframalpha.com/input/?i=solve+a*x_1%5E3+%2B+b*x...
Output: s->a = ((x1 + x2 - 2.) / pow(x2 - x1, 3.)); s->b = ((- (((x1 + x2 - 2.) * pow(x1, 2.)) / pow(x2 - x1, 3.)) - ((4. * x2 * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) + ((6. * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) - ((7. * pow(x2, 2.) * (x1 + x2 - 2.)) / pow(x2 - x1, 3.)) + ((6. * x2 * (x1 + x2 - 2.)) / pow(x2 - x1, 3.)) - 1.) / (4. * x2 - 4.)); s->c = (1. / 2.) * ((((x1 + x2 - 2.) * pow(x1, 2.)) / pow(x2 - x1, 3.)) + ((4. * x2 * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) - ((6. * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) + ((pow(x2, 2.) * (x1 + x2 - 2.)) / pow(x2 - x1, 3.)) - ((6. * x2 * (x1 + x2 - 2.)) / pow(x2 - x1, 3.)) - ((4. * (- (((x1 + x2 - 2.) * pow(x1, 2.)) / pow(x2 - x1, 3.)) - ((4. * x2 * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) + ((6. * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) - ((7. * pow(x2, 2.) * (x1 + x2 - 2.)) / pow(x2 - x1, 3.)) + ((6. * x2 * (x1 + x2 - 2.)) / pow(x2 - x1, 3.)) - 1.)) / (4. * x2 - 4.)) + 1.); s->d = (1. / 4.) * ((((x1 + x2 - 2.) * pow(x1, 3.)) / pow(x2 - x1, 3.)) - ((4. * x2 * (x1 + x2 - 2.) * pow(x1, 2.)) / pow(x2 - x1, 3.)) - (((x1 + x2 - 2.) * pow(x1, 2.)) / pow(x2 - x1, 3.)) - ((pow(x2, 2.) * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) + ((2. * x2 * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) + ((6. * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) + x1 - ((pow(x2, 2.) * (x1 + x2 - 2.)) / pow(x2 - x1, 3.)) + ((6. * x2 * (x1 + x2 - 2.)) / pow(x2 - x1, 3.)) + ((4. * (- (((x1 + x2 - 2.) * pow(x1, 2.)) / pow(x2 - x1, 3.)) - ((4. * x2 * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) + ((6. * (x1 + x2 - 2.) * x1) / pow(x2 - x1, 3.)) - ((7. * pow(x2, 2.) * (x1 + x2 - 2.)) / pow(x2 - x1, 3.)) + ((6. * x2 * (x1 + x2 - 2.)) / pow(x2 - x1, 3.)) - 1.)) / (4. * x2 - 4.)) + 1.);
(https://github.com/albertz/music-player/blob/master/ffmpeg_p... SmoothClipCalc::setX)
For example, in your case you are solving Ax=b which you can easily solve by a simple LU decomposition. Plus if your A changes, your exact formula is wrong, so using an LU decomposition just makes sense. Plus your code will look way cleaner as well. And way easier to debug.
A bit confusing since normally x1,x2 denote the unknowns.
Jargon is fun, unless I want to communicate across fields :/.
Well, to be fair, in that example, if the compiler does those `pow(x2 - x1, 3.)` calls multiple times, this would not be optimal, but otherwise, it should be ok.
I did this because I will probably never ever change the matrix A here and I wanted to make that code very fast. Otherwise, it's of course a good idea to use some more generic solution.
LU decomposition is numerically stable, for the most part. You can run into problems with floating point arithmetic if you're not careful. This just has to do with the way you solve the problem. By breaking up the matrix into an upper and lower triangular part and solving the resulting triangular system, you often avoid things like powers or square roots and just reduce everything to basic addition and multiplication which tend to have more stable numerical properties.
As for speed, well in this small case of N=4 there is probably very little speed increase since LU is O(N^{3}), although this can be improved depending on symmetries of the matrix. Maybe for N=5 I might write out the exact formulas but beyond that I would just a LU solver since the amount of code for an LU solver is fairly compact.
But, as you correctly point out, if you are never changing the matrix A then you are probably safe with writing it out this way.
I just come from a numerical fluids background so seeing exact formulas makes me uneasy. And in my work since the number of grid points N tends to be variable so general formulas are not available. Plus I've noticed a lot of programmers tend to not know numerical algorithms. In fairness, I don't blame them. All the numerical computations courses I've taken and TAed are very boring and don't make you do anything fun and numerical analysis has this reputation of being dry. If you actually made them write stuff like a Navier-Stokes solver, you would get them interested.
By the way if you want to learn more about numerical linear algebra, which is the cornerstone of most scientific and high performance computing, I personally enjoy Trefethen and Bau's book. Although it's aimed at a mathematical audience and assumes such.
There's plenty of cases in graphics and simulation where solving even a 4x4 matrix needs careful numeric consideration (see this paper about tetrahedral mesh simplification, section III.D, page 7: http://www.sci.utah.edu/~hvo/papers/tetstream.pdf)
If you care about accuracy at all, do not use the straightforward "exact formulas" like Cramer's rule. That's just asking for trouble.
And for numerical linear algebra, I'd go start with Strang and chapter 9: http://math.mit.edu/linearalgebra/
And good recommendation by Strang. I believe he has some fantastic MIT OpenCourse lectures on numerical linear algebra as well.
It is most often used as a hand rolled 4x4 matrix inversion function, not a generic NxN matrix inversion.
Yes, it probably is hand written. Or more likely, copied from some well known resource such as MESA.
Computer graphics deals extensively with 4x4 matrices, so there are hand written implementations of elementary operations such as matrix inverse, all unrolled and precomputed so that there's no loops or anything. (NOTE: this isn't really basic loop unrolling, these are not simple loops to begin with).
http://stackoverflow.com/questions/1148309/inverting-a-4x4-m... http://stackoverflow.com/questions/2624422/efficient-4x4-mat...
just because it's viewable, does not mean you are granted the rights to do whatever you want with it.
Ability to do something does not mean it is ok to do it.
As a musician, if I posted a piece of music to a site that permitted it to be played in their Flash player on the site, this does NOT mean I want it to be distributed does it? Unless I explicitly stated that it was alright to distribute it.
EDIT: I didn't mean this to come across as snappily as it did btw!
http://www.keithlantz.net/2011/10/ocean-simulation-part-one-...
There is a larger upfront cost to learning the right language for the right job. However, many people opt to try to fit a square peg into a round hole.
Then again, I might be biased there since I've spent a decent chunk of time porting MATLAB/Octave code to Python/Numpy so that we could get it into production in a robust way..
I thought the only consideration would be speed? The algebra is mainly looping. Isn't that way most languages try to handball algebra to LAPACK/BLAS?
How do you know the author is ok with that?