A problem about the example 'mgpcg_advanced'

[timxu5]

in the example mgpcg_advanced,
the function prolongate is like this:

    @ti.kernel
    def prolongate(self, l: ti.template()):
        for I in ti.grouped(self.z[l]):
            self.z[l][I] = self.z[l + 1][I // 2]

which means: z_l = prolongate(z_(l+1)),
but refer to the paper, “A parallel multigrid Poisson solver for fluids simulation on large grids”,

I think the code should be like this:
z_l += prolongate(z_(l+1)).

    @ti.kernel
    def prolongate(self, l: ti.template()):
        for I in ti.grouped(self.z[l]):
            self.z[l][I] += self.z[l + 1][I // 2]

So I changed the code and the new mgpcg can converge.

I’m very confused which implementation is correct. Why both of them can converge?

[yuanming]

Interesting… Does one of them converge faster than the other?

[timxu5]

Two versions of them take nearly the same time to converge.
In the following images, the upper terminal shows the result using the code:

    @ti.kernel
    def prolongate(self, l: ti.template()):
        for I in ti.grouped(self.z[l]):
            self.z[l][I] += self.z[l + 1][I // 2]

and the bottom terminal shows the results using code

    @ti.kernel
    def prolongate(self, l: ti.template()):
        for I in ti.grouped(self.z[l]):
            self.z[l][I] = self.z[l + 1][I // 2]

When self.N =128, the run time is:

When self.N = 256, the run time is:

While I’m posting the results, I think I can try to explain the reasoning(but all with my imagination).
z is just a preconditioner, so the ‘slight’ difference of z won’t affect the overall algorithm convergence… I’m babbling, plz don’t take it seriously.

[citadel]

The paper one is faster for free surface.

BTW, the factor in restrict() should be 1.0 for 2D.