# Scalar equation in PETSc block matrix system

Dear FreeFem community!

I would like to solve a nonlinear problem using FreeFem and PETSc. The problem consists of 3 [P2,P2,P1] equations and a scalar equation that sets the value of one of the vectorial functions at a single point. The scalar equation and the effect of the scalar value appears as row and column vectors in the block-assembled Jacobian.

I am confused about how the scalar equation should be calculated in the block matrix structure. In the example following example, the rectangular matrix assembly is illustrated well; however, I am not sure how to handle the rows. Should I make vectors from the variational formulation and rearrange their numbering using the āChangeNumberingā command, and then add them to the block matrix? If not, how can I define an (n,1) sized Mat object? I am also confused about changing from PETSc to FreeFem numbering: can this also be done for each [P2,P2,P1] component individually, just like when switching from FreeFem to PETSc?

Andras

This example is probably more meaningful if you want to learn how to add a single constraint.
But yes, you are right. First assemble your vector using a `varf`. Then, go to PETSc numbering with `ChangeNumbering`, and put that in the block `Mat`. You can go from PETSc to FreeFEM no problem using `ChangeNumbering(.., inverse = true);`.

Let me know if you need some more details.

1 Like

Thanks for the help! Based on this info, I managed to improve the code, but it is still not working properly.

When I call the functions for the residual and the Jacobian in the code ,they seem to work. However, when calling SNESSolve, I get the following error message:
current line = 98 mpirank 3 / 4
Assertion fail : (in->n == 0 || in->n == last - first)
line :3570, in file PETSc-code.hpp
err code 6 , mpirank 1

Do you have an idea what causes this error? I am struggling to figure it out.

Also two more small questions, which might relate to the error above, but due to my shallow understanding of FreeFem-PETSc I am not sure:

1. why is the command MatConvert needed? I would think that when I use Mat-s to assemble the block-matrix, it is ready for use by PETSC.
2. What is the suggested way for the variational formulation of setting the value of the function at a certain point? I was thinking about using levelset in the integration for an extremely small circle, or integrating against a dirac-delta function.

Also, I was not precise: the scalar equation is not a constrain, without it, the system would be underdetermined.

Andras

Could you please share a MWE, surrounded by triple backward single quote ```, like this:

``````my running code;
``````
1. `MatConvert` is used to convert a `MatNest` into a `MatAIJ`. That is, an unassembled block operator into a āstandardā assembled distributed operator. That is useful for debugging purposes, because with `MatNest`, you can only use inexact preconditioners such as `PCFIELDSPLIT`, whereas with an assembled operator, you can just use an exact factorization `PCLU` or `PCCHOLESKY` to debug the code.
2. a Gaussian function is a good way to achieve what you want, I think?

Yes, so itās kind of like a constraint indeed, a constraint to make the system definite

I created an MWE. It is unfortunately pretty long, and does not produce the exact same error message; however, the error seems to appear the same way as in my own code. I would really appreciate if you could take a look at it!
BlockPETScMWE.edp (15.3 KB)

Thanks
Andras

Your script is far from being minimalist, I canāt run it myself to help you debug. One obvious issue is that you define `omega` on all ranks, whereas it should only be defined on the first rank of the communicator (usually `mpirank == 0`). You can see in the previous example, for example those lines, there is a `if` to augmente the vector only on a single process. Your line 440 is thus not correct (augmente with `omega` only on rank 0), same with lines 472, 395, 384, 365. If you need `omega` globally, you can just broadcast it from processor 0, `broadcast(processor(0), omega);`.

Thank you so much, what was indeed the main issue!