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?
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);.
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:
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.
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.
Could you please share a MWE, surrounded by triple backward single quote ```, like this:
my running code;
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.
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)
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);.
