Hello,
I wonder what block preconditioners build using CreateBlockSmoother
really do.
In my code, I want to use block Jacobi smoother and I observed that
rho.data = smoother * b
rho.data *= weight
gives different results than
b *= weight
rho.data = smoother * b
Do I miss something? I use quite particular construction of blocks for block smoother (I consider the dofs corresponding to each patch but those on the boundary.). Nevertheless, I would not expect such behavior for any block smoother.
Many thanks for any help,
Jan Papez
BlockSmoother used with the matrix-vector multiplication is Block-Jacobi, so linear.
BlockSmoother used multiplicatively is Block-Gauss-Seidel.
If rho is 0,
smoother.Smooth(rho, b)
and also the symmtric version
smoother.Smooth(rho, b)
smoother.SmoothBack(rho, b)
are linear operations.
Best,
Lukas
Thank you , Lukas. This is exactly what I learnt from the tutorials. However, then I would expect to have the same results (for the codes given in my initial post), which is not the case.
You should get the same result.
Could you upload a failing minimal example?
I enclose an example, which is not very minimal but provides different results for two variants (switch on the line 8 ) of computations on lines 101-106.
( I use a multigrid solver with adaptive stepsize computed on lines 157-160. )
Many thanks for your help!
Jan
https://ngsolve.org/media/kunena/attachments/1395/minimal_example.py
Attachment: minimal_example.py
In variant 2 you modify b, which you seem to use again at another point:
if variant == 1:
rho.data = self.smoothers[level] * b
rho.data *= (1./self.w1) #weightening
else:
b.data *= (1./self.w1)
rho.data = self.smoothers[level] * b
Changing this to:
if variant == 1:
rho.data = self.smoothers[level] * b
rho.data *= (1./self.w1) #weightening
else:
b.data *= (1./self.w1)
rho.data = self.smoothers[level] * b
b.data *= self.w1
solves the problem.
Another thing: It seems that you are only using blocks consisting of one DOF each. Using “.CreateSmoother()” instead of “.CreateBlockSmoother()” might be cheaper.
And one final remark: They way your Multigrid is programmed, you are only solving on the lowest order DOFs. Prolongations only transfer lowest order DOFs, so you need your smoother on the finest level to take care of the high order part by itself. But, the smoother blocks you are using are just one block for each vertex, consisting of the single lowest order DOF for that vertex.
Hope this helps,
Best,
Lukas
Thank you very much, Lukas.
Now I understand the problem. The way the multiplication is coded,
rho.data = pre * res
changes not only rho but also res !
As for your second remark, I used on purpose only lowest order smoother in the minimal example. In the code, I have a higher order smoother on the finest mesh.
Best and thank you once more,
Jan