Tensile instability or something else?

[Qcellaris]

Dear PySPH community,

I was trying to perform some simulations but I kept running into instabilities. I think the simulations that I tried were relatively long and that the problem shows up only if simulations took a certain amount of simulation time. In order to test this hypothesis, I tried running a preinstalled example (couette.py) that used the same scheme (TVFScheme) as I was using for a very long simulation time to see if these instabilities also would show up there. So, the only thing I changed compared to the original preinstalled example was --tf 1920000.0 instead of --tf 100.0.

The results indeed showed the instability that I also encountered in my own simulations, so I know I didn't make a mistake somewhere. Nevertheless, I'm really wondering what this instability might be and how it is caused. At first, to me it seemed to be a tensile instability but tensile instabilities should be suppressed in the TVF scheme. I also tried to increase the background pressure of the TVF scheme to make sure the the effect is strong enough, but this also didn't help. Therefore, I suppose it is not a tensile instability. Furthermore, I also tested the scheme with an angular conserving viscous term in the momentum equation instead of the one used by Adami et al., but unfortunately this didn't help either.

I noticed that he instabilities start near the corners of the periodic domain so might it have something to do with the implementation of periodic boundary conditions in PySPH or with the solid wall boundary conditions as prescribed by Adami et al.? I would like to hear your views on this and am looking forward to a discussion.

You can find the full output at every 800 iterations here: couette_output.zip. To illustrate the issue I have attached two screenshots of the output at different times (one while it is still stable and one where particles are flying around vigorously).

[pawansnegi]

Yes, the problem is due to boundary condition implementation along with periodic BC. The Neumann pressure and no-slip boundary conditions require an extrapolation step just before acceleration computation. Ideally, one should copy the extrapolated properties to the corresponding ones as they have the pressure from the last update only. This does not affect the values for shorter time, but it accumulates for longer time.

An easy fix is do nnps_update=True for the group where solid BC implementation is performed.

A computationally faster fix is to write an equation and copy appropriate properties to the corresponding periodic particle.

[pawansnegi]

Also, the time integrator and the SPH discretizations must be accurate to minimize errors after many iterations. The TVF scheme does not use convergent SPH discretizations.

[Qcellaris]

Thank you for your suggestion. I did a test with the proposed solution and indeed it remains stable for longer. With the nnps_update=False it remained stable up to approximately t=2.3e5 while with nnps_update=True it remains stable up to t=3.2e5 more or less.

I am still wondering if I incorporated your suggestion correctly, see code below. I copied the exact equations as used in the class TVFScheme(scheme) and added the update_nnps=True to the third group. I was wondering if I also have to incorporate it in the second group, because in the second group we extrapolate the fluid velocity to the wall particles and in the third group we extrapolate the fluid pressure to the wall particles. I thought the most important thing is that the nnps is updated before we compute the SolidWallNoSlipBC in the last group.

def create_equations(self):
        equations = [
            Group(
                equations=[
                    SummationDensity(dest='fluid', sources=['fluid', 'channel'])
                ],
                real=False
            ),
            Group(
                equations=[
                    StateEquation(dest='fluid', sources=None, p0=p0, rho0=rho0, b=1.0),
                    SetWallVelocity(dest='channel', sources=['fluid'])
                ],
                real=False
            ),
            Group(
                equations=[
                    SolidWallPressureBC(dest='channel', sources=['fluid'], b=1.0, rho0=rho0, p0=p0, gx=gx, gy=gy, gz=gz)
                ],
                real=False,
                update_nnps=True
            ),
            Group(
                equations=[
                    MomentumEquationPressureGradient(dest='fluid', sources=['fluid', 'channel'], pb=pb, gx=gx, gy=gy, gz=gz, tdamp=tdamp),
                    MomentumEquationViscosity(dest='fluid', sources=['fluid'], nu=nu),
                    SolidWallNoSlipBC(dest='fluid', sources=['channel'], nu=nu),
                    MomentumEquationArtificialStress(dest='fluid', sources=['fluid'])
                ],
                real=True
            )
        ]
        return equations

[pawansnegi]

Hi,

So your time-step dt is 3e-2, and spacing dx, I assume, is around 1e-2. So even if we consider a second-order accurate time-stepping and second-order accurate space discretization, the error will be of order dt^2 * dx^2 = 9e-8. So after 5.3e7 iterations the error accumulation will be of order 4.5e-1 which is huge. Therefore, after so many iterations you should not expect the error to be less.

These were hand wavy calculations. However, I hope you get the point.

If you want to run a lot of iterations, your time integrator and space discretization must be highly accurate.

[Qcellaris]

Yes, I definitely get your point. I'm not planning on running such long simulations, but I was just trying to understand where the instabilities came from and what the issues might be just to get a better understanding. Why is update_nnps=True actually solving this issue? I still don't completely understand what this does. In the integration scheme we have:

class PECIntegrator(Integrator):

    def one_timestep(self, t, dt):
        self.initialize()

        self.stage1()
        self.update_domain()

        self.do_post_stage(0.5*dt, 1)

        self.compute_accelerations()

        self.stage2()
        self.update_domain()

        self.do_post_stage(dt, 2)

So to my understanding, the self.update_domain() is the step that actually creates the ghost particles right? It basically copies particles close to de periodic boundary over to the other side of the domain as ghost particles, including all their properties right? I figured out that in the self.compute_accelerations() step there is also a nnps.update() call. What does this nnps.update() do with the particles right before the accelerations are evaluated? And finally, what does update_nnps=True in a group of equations do to the particles? I can't seem to find a clear explanation in the documentation about this.

P.S. did I incorporate the update_nnps=True correctly in the code block I showed you or should it also be included in the second group?

[pawansnegi]

It is the same thing. Notice that before the acceleration computation, you compute some properties for the solid particle arrays. The ghost particles, which are a copy of actual particles, do not have those updated values. So when you do nnps_update=True, it recreates the particles and thus copies all properties.

Therefore, I said a faster fix would be to write an equation that copies properties from parent to ghost.

The nnps_update=True is required whenever you find that the periodic properties do not have the updated value from the previous group. For the example you showed, I think only one nnps_update is enough.

[pawansnegi]

Yes, I definitely get your point. I'm not planning on running such long simulations, but I was just trying to understand where the instabilities came from and what the issues might be just to get a better understanding. Why is update_nnps=True actually solving this issue?

It provides the ghost particles the updated values from previous group. So, it perform nnps_update().

I still don't completely understand what this does. In the integration scheme we have:

class PECIntegrator(Integrator):

    def one_timestep(self, t, dt):
        self.initialize()

        self.stage1()
        self.update_domain()

        self.do_post_stage(0.5*dt, 1)

        self.compute_accelerations()

        self.stage2()
        self.update_domain()

        self.do_post_stage(dt, 2)

So to my understanding, the self.update_domain() is the step that actually creates the ghost particles right? It basically copies particles close to de periodic boundary over to the other side of the domain as ghost particles, including all their properties right?

Yes.

I figured out that in the self.compute_accelerations() step there is also a nnps.update() call. What does this nnps.update() do with the particles right before the accelerations are evaluated?

Sometimes, you may have some particle shifting to other calculations before actual self.compute_accelerations() in do_post_stage(). Therefore, to ensure updated ghost, nnps_update() is there. you can customize it if you want and make it faster.

And finally, what does update_nnps=True in a group of equations do to the particles? I can't seem to find a clear explanation in the documentation about this.

P.S. did I incorporate the update_nnps=True correctly in the code block I showed you or should it also be included in the second group?

Find in the original reply.

[Qcellaris]

Thank you for your reply. I think I understand it now.

Sometimes, you may have some particle shifting to other calculations before actual self.compute_accelerations() in do_post_stage(). Therefore, to ensure updated ghost, nnps_update() is there. you can customize it if you want and make it faster.

I was indeed thinking about speeding up the simulation by using self.compute_accelerations(update_nnps=False) since we already have self.update_domain() right before the compute step and we don't do anything in the self.do_post_stage() step. This would save at least one domain update.

[Qcellaris]

In order to speed up the simulation I tried using self.compute_accelerations(update_nnps=False) which gave very inaccurate results. Next, I decided to try removing the self.update_domain(), which is just after stage 1 in the integrator, while using self.compute_accelerations(update_nnps=True). Unfortunately, this also gave instabilities. So apparently, both are required together. This indicates to me that either the default self.do_post_stage() does change the particles in some way, or that the update_domain() and nnps_update() are different operations.