Disc patch revived

Implements the changes discussed in #311 (closed).

• The disc patch example is now 400x400x800 in length with correct internal energy.
• The accelerations get truncated at |z-z_disc| > 300 with a cosine transition
• The accelerations are 0 at |z-z_disc| > 380.
• The calculation is significantly faster thanks to the extraction of constant terms.

Could you let me know whether GIZMO does better with this version ?

Added 1 commit:

• 33787d31 - Removed the IDL scripts for the disc patch test case. Updated the python script …

It's still running, but this version seems to have the same problem as the old one.

Everything looks okay after a few snapshots:

But after a while, the spurious velocities start appearing again, and the double peak shows up:

In the old test, the density profile never entirely makes it to what it should be, while now it does temporarily, so there is some improvement. But we still seem to have a steady-state rather than a hydrostatic solution in the outskirts, and (consequently?) the emergence of velocity noise.

Ok. So that's still not enough but going somewhere. Could you try with

 z_trunc:         250.
 z_max:           300.

in the YAML file ? That would leave a larger range without any gravitational forces at all. Should be more than one particle width with these parameters.

Looks better, but still noisy:

The velocity is very close to zero for z<100 and z>700 and there is a clear break in the density pressure profile at these positions, so the cut off seems to work. Something is still seeding the velocity noise in the centre though...

Alright, that's quite a bit better.

Random idea: what if instead of using the forces on the particle you used the potential ? You could then compute the potential gradient across the particle/cell using the same technique than for other quantities. That would allow you to compute the gravity forces at the faces, where the other momentum fluxes are also computed.

Looking at the plot more closely, it looks like something weird happens at the transition between the normal potential and the truncated one. It could be too sharp a transition.

Added 1 commit:

• 2855cd74 - Incorrect division py pi instead of multiplication by pi in the disc patch potential.

Sorry.... I had a typo in my potential. The argument of the cosine truncation was divided by pi instead of multiplied by pi. That implied a large jump in accelerations a z=100 and z=700.

Added 1 commit:

• 1334f2e2 - More soft transition for the truncation of the potential. Set some plotting limi…

The correct potential gives more or less the same result as above (the double peak grows a bit more slowly, but it is still present and is about as strong).

If we were to use the gradient of the potential at the faces to compute the forces, we would have two different values for the gravitational acceleration at each face position (one from the interpolated left state and one from the interpolated right state). How would we decide which value to use? Even then, I don't think the force would balance the momentum flux exactly.

I have tried running it for longer just to see whether things become better or worse.

It seems that things don't grow but the turbulence never settles.

If we were to use the gradient of the potential at the faces to compute the forces, we would have two different values for the gravitational acceleration at each face position (one from the interpolated left state and one from the interpolated right state). How would we decide which value to use? Even then, I don't think the force would balance the momentum flux exactly.

How about just taking the average ? Or can a 1D Riemann solver with gravity be written ?

Okay, then I stick with my theory that the turbulence is somehow driven by the vertical velocity, which should be zero in a hydrostatic solution, but is non zero because the solution is more steady state like. I'm currently trying to write a very basic 1D test code to explore this problem further (with full control over the boundary conditions); I will keep you posted. That would also allow me to experiment with forces evaluated at the faces.

Not sure if a Riemann solver with gravity is possible, I will think about it.

I agree with your diagnostic.

I'll merge the changes to the ICs into master. This is an improvement over the current situation as it works with various box lengths and gets rid of the IDL analysis scripts.

Status changed to merged

mentioned in commit ad8aa1fb

mentioned in issue #311 (closed)

I have done some experiments with a 1D fixed grid finite volume scheme, and it looks like the (original, non suppressed) disc patch solution is very sensitive to the type of boundary condition that is used.

The best solution is obtained when you impose the solution at the boundaries (i.e. zero velocity and the analytic density and pressure at that position):

If periodic boundaries are used, the solution converges to the right answer, but this takes very long (dynamical time is 48 in the same time units as used in the plot), and the velocity stays relatively high throughout:

Inflow/outflow boundaries simply don't work, as in that case you get an infinite amount of mass flowing into your system, and the central density just keeps growing.

I will rerun these tests for the new, suppressed potential, and check if that one works better.