`CreateBlockSmoother` fails Intermittently on well-posed local problems

General
1

Hi,

I’m trying to solve some element-wise nonlinear problems. In NGSolve, I believe there are no trivial ways to define local problems (please let me know if there are any).

So, I decided to make a block solver using the block Jacobi preconditioner, a.mat.CreateBlockSmoother.
But the block Jacobi shows unreliable behavior. It fails even if local problems are non-singular.

I created a simple test case, L2 projection of x^2 + y^2 for DG space with order 2.
The test fails roughly 1 out of 10 runs.
I believe there are some memory leaks or synchronization issues.

from math import sqrt
from ngsolve import BaseMatrix, BilinearForm, L2, Integrate, x, y, SymbolicBFI, GridFunction, CF
from ngsolve.nonlinearsolvers import NewtonSolver
from ngsolve.meshes import MakeStructured2DMesh

class BlockSolver(BaseMatrix):
    def _createBlock(self):
        self.blocks = []
        fes = self.a.space
        for elem in fes.mesh.Elements():
            self.blocks.append(set(d for d in fes.GetDofNrs(elem)))

    def __init__(self, a:BilinearForm):
        super(BlockSolver, self).__init__()
        self.a = a
        self._createBlock()
        self.solver = None # At each Newton iteration, Update will be called

    def Mult(self, x, y):
        self.solver.Smooth(y, x)

    def Update(self):
        # not sure whether this is necessary after the first call.
        self.solver = self.a.mat.CreateBlockSmoother(self.blocks)

if __name__ == "__main__":
    mesh = MakeStructured2DMesh(nx=2, ny=2)

    fes = L2(mesh, order=2, all_dofs_together=True)
    uex = x**2 + y**2 # expect uh == uex

    u, v = fes.TnT()
    a = BilinearForm(fes)
    a += SymbolicBFI((u-uex)*v)

    for i in range(4): # Mesh refinement
        mesh.ngmesh.Refine()
        fes.Update()
        uh = GridFunction(fes)
        uh.Set(CF(0.0))

        local_solver = BlockSolver(a)
        solver = NewtonSolver(a=a, u=uh, solver=local_solver)
        success, it = solver.Solve()
        print(f"Refinement {i}: Newton {it}, error {sqrt(Integrate((uh-uex)**2, mesh))}")
2

You have to initialize the y vector in the Mult function. With

def Mult(self, x, y):
    y[:] = 0
    self.solver.Smooth(y, x)

I get consistently 2 iterations.

3

Ah! Thank you — it now works well on my machine, too.

I found that setting condense=True essentially solves the local problems when there are no degrees of freedom associated with element interfaces. So, if fes consists of element-wise finite element spaces without inter-element continuity, the condensed system becomes a 0x0 matrix, and NewtonSolver handles it gracefully.

A quick performance test shows a 10× speed-up for the L^2-projection compared to the block Jacobi preconditioner approach.