Implement multiphase collisions

[uschille]

Preparations:

• [x] Implement calculation of hydrodynamic moments on the fly.
• [x] Implement the derivatives.
• [x] Implement interface force.

Separate sub-functions in lb_collisions(double *f) for:

• [x] bulk collisions
• [x] interface collisions
• [x] correction collisions

[uschille]

I wonder if we can try a first implementation without the correction collisions...

[uschille]

Since the hydrodynamic moments are needed for the derivatives, it makes sense to compute them on the fly upon lb_write() is called.

[niktre]

It is definitely possible to start without correction operator. As far as I understand, one should see then in a test a drop of "liquid" in "vapor" (or vice versa) with some spurious currents. The implementation of the interfacial force in Jasna's program was modified by me, it should be usable I think. However, the FD must be handled from scratch.

Do you think we may work on this in parallel?

[uschille]

Let's take a shot without correction collisions then. I think all the new collision functions should go into a new file, and perhaps the FD into a separate file. I guess we can work on this in parallel - I'll try to implement the calculation of the moments in the main loop tomorrow. At any rate feel free to go ahead, I can review the code and pull into my own fork as we work in parallel.

[uschille]

The moments are calculated during the update loop now. They reside in memory behind all the populations and write_profile() shows how they can be accessed.

[uschille]

How important do you think is the implicit algorithm? Turns out it is not so easy to integrate as each iteration requires a loop over the lattice.

[niktre]

It should be actually crutial to cancel the spurious currents (if not totally then to a great extent). One of the tests Burkhard mentioned was exactly turning on/off this iteration loop (as well as turning on and off the correction operator itself).

What I still do not now is if the loop may be done locally. I see at least no contradiction there. By this mean we would not need to loop over the whole lattice..

Does it make sense to you?

Am 21.07.15 um 18:09 schrieb Ulf Schiller:

How important do you think is the implicit algorithm? Turns out it is not so easy to integrate as each iteration requires a loop over the lattice.

— Reply to this email directly or view it on GitHub https://github.com/niktre/mlb_jasna/issues/6#issuecomment-123387144.

Dr. Nikita Tretyakov (Polymer Theory Group) Max Planck Institute for Polymer Research Ackermannweg 10 55128 Mainz, Germany Tel. +49-6131-379-473

[uschille]

I'm thinking that perhaps in the stationary state, when you initialize the correction current with the previous value, the implicit algorithm should converge very quickly. So just one step might be good enough an approximation.

I'm not sure it can be done locally. We need the gradient of the tensor in Suppl Eq. (7), which depends in turn on derivatives of the velocity. So this requires at least two neighbourhood stencils and some care with ordering of the data. But perhaps one can plug (6) into (7) into (8) and work out an expression directly with second derivatives.

[niktre]

In the end one has a big matrix where non-zero elements are located on diagonal and its vicinity only. The whole iterative procedure should inverse the whole thing. The question is if physics can help us to avoid this mathematical problem. Burkhard thinks it is unavoidable and one should loose a factor of 100 inside the iteration.

Von meinem iPhone gesendet

Am 21.07.2015 um 18:35 schrieb Ulf Schiller notifications@github.com:

I'm thinking that perhaps in the stationary state, when you initialize the correction current with the previous value, the implicit algorithm should converge very quickly. So just one step might be good enough an approximation.

I'm not sure it can be done locally. We need the gradient of the tensor in Suppl Eq. (7), which depends in turn on derivatives of the velocity. So this requires at least two neighbourhood stencils and some care with ordering of the data. But perhaps one can plug (6) into (7) into (8) and work out an expression directly with second derivatives.

— Reply to this email directly or view it on GitHub.

[niktre]

I think one should get rid of Jasna's code completely. I've put EOS into mlb.c however I am a bit confused with R1 and R0. In App. B they define c^2_s(\rho_1) = 0.6, but in Eq. B9 there is c^2_s(\rho_0).

I've branched with an EOS that I think corresponds to the paper 1 to 1, if I haven't done mistake upon integration (which happens). Can you access this branch?

[uschille]

I'll have a look at the EOS branch.

Did Jasna's code ever work in any way? There are some parts that I don't really understand...

[niktre]

I think it has never worked. Burkhard claims it kind of worked before she has added the implicit convergence for correction currents, but even in that case the change of the number of grid points killed everything. They basically saw a lot of (sound) waves traveling in x- and y-directions

[uschille]

Ok. I think the implicit algorithm can not have worked because I don't see periodic boundary conditions and the velocities are overwritten while the old values are still needed for gradients (although the latter may not be so critical). Anyhow, so we don't have a reference to compare to, which will make debugging a bit non-trivial.

[uschille]

I agree the EOS is confusing. Can it be that at some point they used a different function than the one in the paper?

[niktre]

Yes, they've changed it a couple of times (or at least once). Jasna's version is based not on sin but on straight lines.

Von meinem iPhone gesendet

Am 26.07.2015 um 00:14 schrieb Ulf Schiller notifications@github.com:

I agree the EOS is confusing. Can it be that at some point they used a different function than the one in the paper?

— Reply to this email directly or view it on GitHub.

[uschille]

Reopening for further discussion and testing.

[uschille]

I believe the EOS is correct now, but even without all the multiphase collisions and only the bulk collisions enabled the system becomes unstable. It looks like the spurious currents are rather larger, too large for LB. So probably there's a bug in the correction current.

[niktre]

Try with S1 0.5. It makes everything much more stable for me.

[uschille]

Thanks. So without the correction current it is running, but with the correction current it goes bonkers.

[uschille]

The problem comes from Xi[i][j][k], when I skip it in the correction collisions it works better. See comment in ccbc90965c36d79a0b19f74af2b741a7e0810c5f.

[uschille]

As of 481dd9b0b8fd45202412cd1f0de48456cd93a9cb I think we need to look into the implicit algorithm. Somehow after some iterations it becomes unstable. What is strange is that when I change the sign of the correction current it seems to be more stable.

[uschille]

So, I'm wondering again whether the jump in the second derivative of the pressure could cause problems. Essentially, if the density is just slightly different between two lattice sites, the correction tensor Sigma will reflect this jump. Not sure what that does then.

[niktre]

I was talking to Jasna at some point about this. She wasn't sure as well whether it is good to have such a jump. Maybe, it is better to ask Burkhard directly about this.

Am 29.07.15 um 18:26 schrieb Ulf Schiller:

So, I'm wondering again whether the jump in the second derivative of the pressure could cause problems. Essentially, if the density is just slightly different between two lattice sites, the correction tensor Sigma will reflect this jump. Not sure what that does then.

— Reply to this email directly or view it on GitHub https://github.com/niktre/mlb_jasna/issues/6#issuecomment-126007565.

[uschille]

As of 711e55708af257ab6d08fcbbba38b0efb4206d0d there is an alternative EOS with continuous second derivatives. I am currently experimenting with the parameters.