PETSc SNES implementation in deal.ii

[Lorenzo Mingazzi]

Dear all,

I am currently trying to implement the PETSc SNES solver  (https://petsc.org/release/docs/manual/snes/)  into my code.
As a starting point, I am trying to implement it into the step-15 from the tutorials.

Attached can be found a MWE of the current version of the code which follows an older implementation (https://github.com/gujans/dealii-snes/blob/master/step-15.cc).
The code performs the following steps:
1 - create the mesh;
2 - initialize the PETScWrappers::MPI matrix/vectors;
3 - set the boundary values;
4 - initialize the SNES solver;
5 - assemble the FormFunction and FormJacobian using the assemble_rhs and assemble_system as functions to set the values to the system matrix and rhs;
6 - solve the problem using the SNES solver.

At the moment the code runs but it produces a solution which is constant zero, moreover the SNES iteration count returns zero as well.
Am I missing something while applying the boundary conditions?

How can I pass the correct assembled matrices to the SNES Solver?

To be noted that if I pass to the solver the PETScWrappers::MPI::Vector present_solution, which is initialized during the setup_system(), instead of a Vec test_x which i created just to store the solution from the solver, it returns the following error:
[0]PETSC ERROR: ----- Error Message --------

[0]PETSC ERROR: Object is in wrong state
[0]PETSC ERROR: Not for unassembled vector

Lastly, I also found an old post in the group (https://groups.google.com/g/dealii/c/pg9bj_z8dbk/m/OTtfdCY0AAAJ) which hinted another possible approach to the one i am using, however i am not sure how to extract PETSc Vec object from PETScWrappers::MPI::Vector without passing from a PETScWrapper::VectorBase.

Thanks in advance,
Lorenzo

[Timo Heister]

Hi Lorenzo,

I didn't look at your code, but I can help with this:

> however i am not sure how to extract PETSc Vec object from PETScWrappers::MPI::Vector without passing from a PETScWrapper::VectorBase.

You should be able to use the conversion operator
https://www.dealii.org/developer/doxygen/deal.II/classPETScWrappers_1_1VectorBase.html#ade5e08d408ee79da10cceb7642bfa1da

PETScWrappers::MPI::Vector x;
Vec px = static_cast<Vec>(x);

[Lorenzo Mingazzi]

Hi,

thank you very much for the reply, I completely missed that cast operator into a Vec.

I will try to see if i can make it work properly.

Best regards, 

Lorenzo