Kuramoto-Sivashinsky equation Pattern idea

[danwills]

More of a question than an issue, How doable might a Kuramoto-Sivashinsky Formula/Pattern be to implement in Ready? If it's simple enough I'm hoping that I could even do it myself with a few pointers?

https://en.wikipedia.org/wiki/Kuramoto%E2%80%93Sivashinsky_equation

https://encyclopediaofmath.org/wiki/Kuramoto-Sivashinsky_equation

From the looks of the formula (and example output on youtube)

https://www.youtube.com/watch?v=OWYz3bVKl6k https://www.youtube.com/watch?v=FO3mYe8zkrM

it seems like something that Ready could probably do, laplacian, bilaplacian and so on.

I have seen some solutions use really complicated-sounding integration schemes like "2nd order Runge-Kutta Exponential Time Differencing" (ETDRK2) or even ETDRK4! is this kind of thing needed for this formula for some reason or can it be converted into the usual reagent-deltas world of a Ready Formula-mode pattern?

Any help or suggestions would be enormously appreciated!

[timhutton]

Hi Dan,

That looks great.

My reading of it for a 2D system is:

du/dt = - laplacian(u) - bilaplacian(u) - ( (du/dx)^2 + (du/dy)^2 ) / 2

but there seems to be some confusion about the sign and magnitude of the gradient term.

Start with Pennybacker2013/spots.vti which already has these pieces.

It may be very 'stiff' and thus require tiny timesteps but I think it should work in Ready.

Go for it!

Tim

On Thu, 5 Nov 2020 at 10:15, Dan Wills notifications@github.com wrote:

More of a question than an issue, How doable might a Kuramoto-Sivashinsky Formula/Pattern be to implement in Ready? If it's simple enough I'm hoping that I could even do it myself with a few pointers?

https://en.wikipedia.org/wiki/Kuramoto%E2%80%93Sivashinsky_equation

https://encyclopediaofmath.org/wiki/Kuramoto-Sivashinsky_equation

From the looks of the formula (and example output on youtube)

https://www.youtube.com/watch?v=OWYz3bVKl6k https://www.youtube.com/watch?v=FO3mYe8zkrM

it seems like something that Ready could probably do, laplacian, bilaplacian and so on.

I have seen some solutions use really complicated-sounding integration schemes like "2nd order Runge-Kutta Exponential Time Differencing" (ETDRK2) or even ETDRK4! is this kind of thing needed for this formula for some reason or can it be converted into the usual reagent-deltas world of a Ready Formula-mode pattern?

Any help or suggestions would be enormously appreciated!

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/GollyGang/ready/issues/85, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAE57ND7UC4LWWAN7YFLMD3SOJ3L7ANCNFSM4TLFEDZA .

-- Tim Hutton - https://twitter.com/_tim_hutton_ - https://github.com/timhutton

[danwills]

I made a pattern exactly as you suggested and what an amazing suggestion for where to start! as you said all the bits were there and easy to snap together.

And It looks like it's working!! hooray! But currently the values seem to keep on getting more and more negative forever..(at least as far as I've gone, which led to setting the low range value as far as -1000! Ready's range-fitting with negative values seems a bit odd though.

I had a play with adding a constant in to kinda 'fight against the flow' but that didn't seem to help, yet. Would love to hear your intuition about what changes might allow it to be kinda DC-stabilized around some fixed value?

(zipped vti is attached)

Kuramoto-Sivashinsky.tar.gz

[timhutton]

Oh sweet. Nice job Dan.

In the first youtube it says they are plotting h-h_Bar, where presumably h_Bar is the mean of all the h values. I'm not sure how to do that in Ready but I'm sure we can think of a way.

What do you mean about Ready's range-fitting? We don't try to do any kind of automatic range estimation. (Maybe we should.)

[danwills]

Cheers Tim, so cool that it worked! Aah it was 'h minus hbar', makes perfect sense then! Yeah it would be awesome if we can work out how to do that in Ready, and automatic range-fitting too, would be hugely helpful in lots of cases (especially when the ranges of each reagent are wildly different).

Will describe what I found odd about the fitting in a bit.

[danwills]

I had some luck (after a bit of mucking around with things that would crash or insta-NaN-explode) with adding an additional reagent to control the general direction of evolution of the KSE.

It's a bit touchy, but the heavy diffusion of 'b' in this version seems to reduce the chances of it degenerating into a tight-spirals dynamic that seems to occur when the KSE evolution is able to outrun the 'b' diffusion.

This does regular KSE up to a threshold (in this case, 50.0), and then once 50 is exceeded, b starts dropping, and when b becomes negative, the evolution of 'a' swaps to downward too, until the threshold is exceeded on the negative side, at which point it reverses and repeats the cycle again.

Kuramoto-Sivashinsky-pocket-c.tar.gz

This is pretty cool and all but I still feel like there might be some term that we could add to the main formula to make it stay centered on some value rather than always going up or down.

I will contact one of the authors of the youtube videos: Steffen Richters and see if they have any ideas : )

[danwills]

Couple of screenshots of the 'pocket' variant.

[danwills]

Steffen hasn't responded yet.. I hope I have the right address.

I also wanted to take the chance refer to this 'Kuramoto-Sivashinsky Equation' thing as "KSE" specifically, to see if it might eventually show up on a google search for "KSE" ; )

[danwills]

I got a message back from Steffen (I think the 1st one was just lost) and awesomely he was able to give me advice that led to a naturally-stabilized version of KSE!

I also turned up the timestep a bit and turned down the number of steps so that it's a bit smoother but still has a good pace:

Kuramoto-Sivashinsky-stabilize2.vti.tar.gz

I also asked Steffen about what might give more control or differently-behaved dynamics and he suggested potentially adding in some amount of the tri-laplacian (as in the 'Nikolaevskiy Equation': https://www.youtube.com/watch?v=fmljvzy77KI this video is also by Steffen).

I'd love to experiment with adding this in but I'm not sure how to write it. I get that the bi-laplacian is the laplacian-of-the-laplacian, so this sounds like the laplacian-of-the-laplacian-of-the-laplacian, which sounds only slightly crazy but I am still not confident enough with the math/etc to do anything other than play, and I think this most likely needs an exactly-correct solution.

If anyone knows of any example code (for trilaplacian) anywhere that I could refer to, that would be utterly choice! : )

[timhutton]

Hi Dan. This looks great! Will you make a pull request to get it in?

[timhutton]

A tri-laplacian should be straightforward. The key thing to understand is convolution. Here how it works:

If you look in the code for this pattern you will see the 5-point stencil for the laplacian: [ 0,1,0; 1,-4,1; 0,1,0 ]. This means:

0  1  0
1 -4  1
0  1  0

To apply this to an image, place the stencil on top of the image and center it on the pixel you want to compute. Then multiply every pixel and add up the sum.

Example:

Our image is: 2  3  2
              4  1  0
              1  5  3

Our stencil is: 0  1  0
                1 -4  1 (and then divide by h^2)
                0  1  0

To find the new image we convolve it with the stencil. Place the stencil on top and multiply:

2 * 0   3 * 1   2 * 0              0   3  0
4 * 1   1 *-4   0 * 1     gives    4  -4  0
1 * 0   5 * 1   3 * 0              0   5  0

Add all the entries and that's the new pixel value in the center:

0 + 3 + 0 + 4 - 4 + 0 + 0 + 5 + 0 = 8  (and then divide by h^2)

Move the stencil onto a new pixel on the original image and repeat. When you have done every pixel you will have applied the laplacian to the whole image.

For a bi-laplacian we can compute the stencil by convolving the laplacian stencil with itself. You'll need to add zeros around the edge so that there's enough room. Try it. Check that you get the same result as you can see in the code:

0  0  1  0  0
0  2 -8  2  0
1 -8 20 -8  1 (and then divide by h^4)
0  2 -8  2  0
0  0  1  0  0 

To get the tri-laplacian stencil, convolve the bi-laplacian stencil with the laplacian stencil. You should have a 7x7 grid of numbers, all divided by h^6.

To get the tri-laplacian running in Ready, you'll need to collect all the pixels to apply it to. That's the xm1, xm2, ..., a_w2 code. You'll need to extend to xm3, ..., a_w3 to cover the 7x7 area.

[timhutton]

(There are better stencils than this but this is the simplest.)

[danwills]

G'day Tim,

Thanks heaps for taking a look! (yes i'll absolutely make a pull request) and for the great explanation of convolving the kernel with itself! I had wondered whether that was how it worked, so it was very nice to hear it confirmed.

I'll get the existing 3 into a pull request: Straight KSE, Stabilised KSE and the 'Pocket' experiment with the extra reagent. Do you think they should go in an author-named folder (KuramotoSivashinsky/), or perhaps in a subdir of Experiments/DanWills?

I'll also do some experimentation with adding higher order Laplacian to the formula, maybe with that in there we'll be able to add an example for "Soft Mode Turbulence" / Nikolaevskiy to Ready too! that would be cool.

I'm also so keen for the idea in #87 (providing some of these extra convolutions/derivatives in Formula mode). It'd be a really nice simplification to be able to convert the KSE formula to just the main line! : ) I'll see if I can contribute any progress in that direction.

[danwills]

Actually I think I misunderstood, reading wikipedia (https://en.wikipedia.org/wiki/Del#Laplacian) which says that nabla^2 (also known as del^2) is laplacian. I'm also watching a bunch of Khan Academy to try to understand better what this means, over on Khan laplacian is shown as the divergence of the gradient of F (and seeing that visualised gave a very nice intuitive geometrical interpretation. too)

I reckon I'll implement a thing for making and combining/convolving kernels with each other in Houdini, sounds like it could be useful (and fun to play with) beyond this nabla^6 thing.

[timhutton]

I agree that #87 is an important thing to do. I'm working on that now.

Update: See PR #90 - KSE can now be a 1-line formula!

[danwills]

Ah bloody awesome! Great news Tim, was trying to checkout the branch yesterday, but no need now! I'll crack on with #89 now.

[timhutton]

@danwills KSE works nicely in 3D and 1D. I've added suggested values of stabilize for these.

There is some really interesting behavior at small sizes with toroidal wrap-around (as is default):

1D: At 32x1x1 or 64x1x1 the whole pattern stabilizes pretty quickly into fixed waves. At 128x1x1 it can take a long time but happens very occasionally.

2D: At 32x32x1 it stabilizes but then flips between several states with hexagonal or square symmetries. I've never seen anything like this before - is this known behavior? Video: https://www.youtube.com/watch?v=32shO3hRpV0

3D: Hints of multi-stability at 32x16x16 but it seems more chaotic.

I'll try to verify at different values of dx, to make sure it's not an artefact of our stencils or something. Update: yes! Still works (64x64x1, dx=0.25, timestep=0.00025, stabilize=-0.05). This is a real phenomenon.

I think it's resonance, like standing waves in a pool - but bistable. What is this!?

If this was a bottle that you were blowing across, it would make one note for a while then switch to another note, then switch back after a while.

The resonance hypothesis seems to be supported by trying it in a rectangle - the phenomenon stops working.

Possibly this paper covers it: Secondary instabilities in the stabilized Kuramoto-Sivashinsky equation Chaouqi Misbah and Alexandre Valance Phys. Rev. E 49, 166 – Published 1 January 1994

[danwills]

Crikey Tim that is quite an interesting discovery! Makes me wonder whether KSE can perhaps somehow be made locally-resonant (in a larger-field (ie not just these tiny ones)) to create more interesting/larger patterns inspired by the same phenomena?

I took a look at your patterns and it's fascinating! (I was just now lucky enough to see westward-traveling hexagons form from a random initial state! The transitions are so beautiful! especially with the 3d view there as well!) Just, wicked!!! I cannot wait to couple this thing to some other RDs as well that is gonna be well beyond wave-equation possibilities.

We/I should crack on (especially now that you've already added in the trilaplacian keyword) with trying to get a 'Soft Mode Turbulence' / nikolaevskiy (Spelling probably wrong) formula going too!

Such awesome work lately @timhutton ! huge respect and thanks!

I would be super interested to hear whether that paper is talking about the same thing. Maybe I should see if Steffen might have any thoughts about it?