# Scalar equation in PETSc block matrix system

**URL:** <https://community.freefem.org/t/scalar-equation-in-petsc-block-matrix-system/796>\
**Category:** General Discussion\
**Created:** [February 11, 2021, 1:41pm UTC](https://community.freefem.org/t/scalar-equation-in-petsc-block-matrix-system/796 "2021-02-11T13:41:39Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![aszaboa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/aszaboa/32/2918_2.png) [@aszaboa](https://community.freefem.org/u/aszaboa)\
**Post date:** [February 11, 2021, 1:41pm UTC](https://community.freefem.org/t/scalar-equation-in-petsc-block-matrix-system/796/1 "2021-02-11T13:41:39Z")

</div>

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](https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/stokes-block-2d-PETSc.edp), 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?

Thanks in advance  
Andras

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [February 11, 2021, 2:33pm UTC](https://community.freefem.org/t/scalar-equation-in-petsc-block-matrix-system/796/2 "2021-02-11T14:33:44Z")

</div>

This [example](https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/laplace-lagrange-PETSc.edp) 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.

---

<div class="post-metadata">

**Author:** ![aszaboa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/aszaboa/32/2918_2.png) [@aszaboa](https://community.freefem.org/u/aszaboa)\
**Post date:** [February 12, 2021, 9:48am UTC](https://community.freefem.org/t/scalar-equation-in-petsc-block-matrix-system/796/3 "2021-02-12T09:48:36Z")

</div>

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.

Thanks in advance  
Andras

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [February 12, 2021, 9:55am UTC](https://community.freefem.org/t/scalar-equation-in-petsc-block-matrix-system/796/4 "2021-02-12T09:55:47Z")

</div>

Could you please share a [MWE](https://en.wikipedia.org/wiki/Minimal_working_example), surrounded by triple backward single quote ```, like this:

```auto
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 😃

---

<div class="post-metadata">

**Author:** ![aszaboa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/aszaboa/32/2918_2.png) [@aszaboa](https://community.freefem.org/u/aszaboa)\
**Post date:** [February 12, 2021, 10:38pm UTC](https://community.freefem.org/t/scalar-equation-in-petsc-block-matrix-system/796/5 "2021-02-12T22:38:38Z")

</div>

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](https://community.freefem.org/uploads/short-url/871eQLnYUX44jbHxEe9KCxysrRK.edp) (15.3 KB)

Thanks  
Andras

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [February 13, 2021, 8:53am UTC](https://community.freefem.org/t/scalar-equation-in-petsc-block-matrix-system/796/6 "2021-02-13T08:53:15Z")

</div>

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](https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/laplace-lagrange-PETSc.edp#L38-L39) 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);`.

---

<div class="post-metadata">

**Author:** ![aszaboa](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/aszaboa/32/2918_2.png) [@aszaboa](https://community.freefem.org/u/aszaboa)\
**Post date:** [February 18, 2021, 5:32pm UTC](https://community.freefem.org/t/scalar-equation-in-petsc-block-matrix-system/796/7 "2021-02-18T17:32:30Z")

</div>

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