Should the init method of PETScSNESSolver be virtual?

Dear FEniCS developers,

I am trying to use the PETScSNESSolver class to solve a problem where the residual and jacobian callback functions return assembled tensors. The most simple rationale that comes to mind is to copy tensors that are provided by my callbacks to internal PETSc data structures.

The problem that I face is that, in order to be able to copy the jacobian matrix, its sparsity pattern should have been already defined in PETScSNESSolver::init when calling SNESSetJacobian. I thought of working around that inheriting from PETScSNESSolver and overriding its init method. However, PETScSNESSolver::init is not virtual.

I attach four files. The last one 3_patch is a simple patch to make PETScSNESSolver::init virtual Files 0_reference.py and 1_residual_jacobian_as_forms.py are to be considered as reference solutions (adapted from the FEniCS example on nonlinear Poisson problem) Finally, file 2_residual_jacobian_as_tensors.py provides an example of jacobian and residual callbacks that return assembled tensors. PETScSNESSolver::init is overridden, although admittedly in a hacky way. Look for comments starting with XXX for more details. Is there a better way to handle this case? If not, should at least the init method of PETScSNESSolver be virtual?

Best regards, Francesco

Comments (2)