Preconditioning a VI restricted system with the full unrestricted operator

[colinjcotter]

I am interested in using iterative methods with the VI SNES methods in PETSc. These methods maintain an "active set", a subset of the DOFs which can change from one SNES iteration to the next, and require that the Jacobian operator is restricted to the active set for the Newton update.

I'm solving a mixed system with MG+PatchPC which is executed for the full system. I think what is needed is my PC P^{-1} to be replaced with AP^{-1}A^T where A is the restriction operator from the full DOFs set to the active set.

Is this the right way to do it? What is the right way in Firedrake?

[ioannisPApapadopoulos]

Hi @colinjcotter! @pefarrell and I used MG+PatchPC with vinewtonrsls in this paper: https://doi.org/10.1137/22M1478598

Unfortunately I never merged the code required into Firedrake and as far as I can tell, it does not work with the current version of Firedrake.

A is the restriction operator from the full DOFs set to the active set

I think one would want the restriction to the inactive set, right?

In general I think there are the following steps to get this working:

1. If required, implement a matrixfree method for vinewtonrsls. I needed this since I had a real block in my system and Firedrake would not assemble it. If you can skip this step then great. At the time I did it by modifying vinewtonrsls.c and leveraging the DirichletBC interface, but Lawrence suggested an alternative here: matfree vinewtonrsls #2549

2. Decide how to pick your active set on the different MG levels. One option is directly by the active set definition with your injected functions. However, we found that definition resulted in poor iteration counts. We had a piecewise constant DG0 discretization and found that a better choice turned out to be the following: a coarse level cell is in the active set if half or more of its parent cells were also in the active set. I imagine this would also work for CG1 dofs as well. Note that vinewtonrsls does not work higher order discretizations anyway.

Here are the PCs I created for the injection of the active set and functions (specific to my application): https://github.com/ioannisPApapadopoulos/fir3dab/blob/main/fir3dab/dabMGPC.py

[colinjcotter]

Thanks John! This is super useful.

[ioannisPApapadopoulos]

I am happy to talk more about it if you like. I've been meaning to invest some time to update fir3dab. Let me know if that's something you would be interested in.

[colinjcotter]

Yes I would be

[colinjcotter]

A couple of (probably related) questions:

1. How are you getting the active set? I read the code but couldn't see it.
2. What is the numbering from the full set of local DOFs to the active set DOFs?

[ioannisPApapadopoulos]

The active set is defined in PETSc and then is passed, via a call to my matrix-free implementation of zeroRowsColumns, to an object I made called ZeroRowsColumnsBC. This contains the indices of the active set that I want to pretend are additional DirichletBC values.

It might be easiest if I write down the flow of the code. Recall that I have modified operators.py in https://github.com/ioannisPApapadopoulos/fir3dab/blob/main/firedrake-install/firedrake.diff and virs.c in https://github.com/ioannisPApapadopoulos/fir3dab/blob/main/firedrake-install/petsc.diff.

So in virs.c, the active set is kept in IS_act. Then the following call is made ierr = MatZeroRowsColumnsIS(jac_inact_inact, IS_act, 1.0, NULL, F_inact);CHKERRQ(ierr); which in turn, via petsc4py, calls its matrix free implementation in operators.py which is zeroRowsColumns(self, mat, active_rows, diag=1.0, x=None, b=None). So the active_rows argument is precisely the PETSc IS_act object. Now this needs to converted into the numbering that a DirichletBC requires and this is achieved via a call to find_element_of_which_sub_block(rows, ises) which again is something in https://github.com/ioannisPApapadopoulos/fir3dab/blob/main/firedrake-install/firedrake.diff

The following lines of code are probably the most relevant:

+        The way we zero rows and columns of unassembled matrices is by
+        constructing a DirichetBC corresponding to the rows and columns.
+        By nature of how bcs are implemented, DirichletBC is equivalent to
+        zeroing the rows and columns and adding a 1 to the diagonal
+        """
+        #print("Inside zeroRowsColumns")
+        if not numpy.allclose(diag, 1.0):
+            # DirichletBC adds a 1 onto the diagonal, this is part of the implementation and is not easy to change
+            raise NotImplementedError("We do not know how to implement matrix-free zeroRowsColumns with diag not equal to 1")
+        comm  = self._y.function_space().dof_dset.field_ises[0].comm.tompi4py()
+        ises = self._y.function_space().dof_dset.local_ises
+ 
+        lgmap = self._y.function_space().dof_dset.lgmap
+        brows = lgmap.applyInverse(active_rows)
+
+        # Unsure why but for parallel we need to gather global indices
+        active_rows = numpy.concatenate(comm.allgather(active_rows))
+        # active_rows is global indices, need to do a global to local mapping for DirichletBC
+        active_rows = lgmap.applyInverse(active_rows)
+ 
+        # Find the blocks which the rows are a part of and find the row shift
+        # since DirichletBC starts from 0 for each block
+        
+        if len(ises) > 1:
+            # This is slow, ideally we do not do it, but how else to know which field is where?
+            (block, shift) = find_element_of_which_sub_block(active_rows, ises)
+        else:
+            block = {0: list(active_rows)}
+            shift = [0]

Then when defining things for MG, the active rows are grabbed from the ZeroRowsColumnsBC class which was created when the call to zeroRowsColumns was made. E.g. here https://github.com/ioannisPApapadopoulos/fir3dab/blob/fc81dd7d7ef2cd07a3a29f6014115e170f61fd96/fir3dab/dabMGPC.py#L167-L170

[rckirby]

I'm keenly interested in this issue, as @danshapero and are are using VIs to enforce bounds constraints: https://link.springer.com/article/10.1007/s00211-024-01405-y.

We are currently just using vinewtonrsls with LU, which has its limitations. However, we are finding that this strategy does work for quadratic element. For cubics, it worked sometimes but failed to converge in some examples as well.

Another issue to discuss as part of this would be composition with the P1 preconditioner somehow.

[colinjcotter]

nice paper! I hadn’t spotted that one