[petsc-users] Nullspaces

Wed Jan 19 08:15:31 CST 2022

ILU is LU for this 1D problem and its singular.
ILU might have some logic to deal with a singular system. Not sure. LU
should fail.
You might try -ksp_type preonly and -pc_type lu
And -pc_type jacobi should not have any numerical problems. Try that.

On Wed, Jan 19, 2022 at 7:19 AM Matthew Knepley <knepley at gmail.com> wrote:

> On Wed, Jan 19, 2022 at 4:52 AM Marco Cisternino <
> marco.cisternino at optimad.it> wrote:
>
>> Thank you Matthew.
>>
>> But I cannot get the point. I got the point about the test but to try to
>> explain my doubt I’m going to prepare another toy code.
>>
>>
>>
>> By words…
>>
>> I usually have a finite volume discretization of the Laplace operator
>> with homogeneous Neumann BC on an octree mesh and it reads
>>
>> Aij * xj = bi,
>>
>> being the discretization of
>>
>> Int|Vi(nabla^2 pi dV) = Int|Vi(nabla dot ui)
>>
>> (I omit constants), where Int|Vi(…dV) is a volume integral on the I cell,
>> pi is cell pressure, ui is the cell velocity.
>>
>> The computational domain contains 2 separated sub-domains.
>>
>> Let’s consider 4 cells into the whole domain and 2 cells for each
>> sub-domain.
>>
>> I would expect a null space made of 2 vectors and from your patch they
>> should look like [1/sqrt(2) 1/sqrt(2) 0 0] and [0 0 1/sqrt(2) 1/sqrt(2)],
>> i.e. norm2 = 1 for both.
>>
>> With MatNullSpaceCreate(getCommunicator(), PETSC_TRUE, 0, nullptr,
>> &nullspace) I’m saying that [1/sqrt(4) 1/sqrt(4) 1/sqrt(4) 1/sqrt(4)] is
>> the null space, which is not but it is in the null space.
>>
>> But this is not the case I sent you, not exactly.
>>
>>
>>
>> The case I sent is 1/Vi * Aij * xj = 1/Vi bi, where Vi is the volume of
>> the cell i.
>>
>> Let’s admit that yj is in the null space of of Aij, it should be in the
>> null space of 1/Vi * Aij, therefore Aij*yj = 0 and 1/Vi * Aij*yj = 0 too.
>>
>> But in the framework of the test, this is true with infinite precision.
>>
>> What happens if norm2(Aij*yj) = 10^-15 and Vi = 10^-5? Norm2(1/Vi * Aij
>> * yj) = 10^-10!!! Is yi still in the null space numerically?
>>
>> Let’s say yi is the constant vector over the whole domain, i.e.
>> [1/sqrt(4) 1/sqrt(4) 1/sqrt(4) 1/sqrt(4)]. Should this be in the null space
>> of 1/Vi * Aij, shouldn’t it?
>>
>> An analogous argument should be for the compatibility condition that
>> concerns bi.
>>
>>
>>
>> My current problem is that properly creating the null space for Aij, i.e.
>> [1/sqrt(2) 1/sqrt(2) 0 0] and [0 0 1/sqrt(2) 1/sqrt(2)], allows me to solve
>> and find xi, but multiplying by 1/Vi, I cannot get any solution using both
>> FGMRES+ILU and FGMRE+GAMG.
>>
>>
>>
>> The tiny problem will load Aij, Vi and bi and show the problem by testing
>> the proper null space and trying to solve. I will include the patch to my
>> PETSc version. I hope to come back to you very soon.
>>
>
> This sounds like a dimensionalization problem to me. It is best to choose
> length (and other) units that make the matrix entries about 1. It seems
> like you are far from this, and it is causing a loss of accuracy (your
> null vector has a residual of about 1e-7).
>
>   Thanks,
>
>       Matt
>
>
>> Thank you very much for your support!
>>
>>
>>
>>
>>
>> Marco Cisternino
>>
>>
>>
>> *From:* Matthew Knepley <knepley at gmail.com>
>> *Sent:* martedì 18 gennaio 2022 21:25
>> *To:* Marco Cisternino <marco.cisternino at optimad.it>
>> *Cc:* petsc-users <petsc-users at mcs.anl.gov>
>> *Subject:* Re: [petsc-users] Nullspaces
>>
>>
>>
>> I made a fix for this:
>>
>>
>>
>>   https://gitlab.com/petsc/petsc/-/merge_requests/4729
>>
>>
>>
>>   Thanks,
>>
>>
>>
>>      Matt
>>
>>
>>
>> On Tue, Jan 18, 2022 at 3:20 PM Matthew Knepley <knepley at gmail.com>
>> wrote:
>>
>> On Thu, Dec 16, 2021 at 11:09 AM Marco Cisternino <
>> marco.cisternino at optimad.it> wrote:
>>
>> Hello Matthew,
>>
>> as promised I prepared a minimal (112960 rows. I’m not able to produce
>> anything smaller than this and triggering the issue) example of the
>> behavior I was talking about some days ago.
>>
>> What I did is to produce matrix, right hand side and initial solution of
>> the linear system.
>>
>>
>>
>> As I told you before, this linear system is the discretization of the
>> pressure equation of a predictor-corrector method for NS equations in the
>> framework of finite volume method.
>>
>> This case has homogeneous Neumann boundary conditions. Computational
>> domain has two independent and separated sub-domains.
>>
>> I discretize the weak formulation and I divide every row of the linear
>> system by the volume of the relative cell.
>>
>> The underlying mesh is not uniform, therefore cells have different
>> volumes.
>>
>> The issue I’m going to explain does not show up if the mesh is uniform,
>> same volume for all the cells.
>>
>>
>>
>> I usually build the null space sub-domain by sub-domain with
>>
>> MatNullSpaceCreate(getCommunicator(), PETSC_FALSE, nConstants, constants,
>> &nullspace);
>>
>> Where nConstants = 2 and constants contains two normalized arrays with
>> constant values on degrees of freedom relative to the associated sub-domain
>> and zeros elsewhere.
>>
>>
>>
>> However, as a test I tried the constant over the whole domain using 2
>> alternatives that should produce the same null space:
>>
>>    1. MatNullSpaceCreate(getCommunicator(), PETSC_TRUE, 0, nullptr,
>>    &nullspace);
>>    2. Vec* nsp;
>>
>> VecDuplicateVecs(solution, 1, &nsp);
>>
>> VecSet(nsp[0],1.0);
>>
>> VecNormalize(nsp[0], nullptr);
>>
>> MatNullSpaceCreate(getCommunicator(), PETSC_FALSE, 1, nsp, &nullspace);
>>
>>
>>
>> Once I created the null space I test it using:
>>
>> MatNullSpaceTest(nullspace, m_A, &isNullSpaceValid);
>>
>>
>>
>> The case 1 pass the test while case 2 don’t.
>>
>>
>>
>> I have a small code for matrix loading, null spaces creation and testing.
>>
>> Unfortunately I cannot implement a small code able to produce that linear
>> system.
>>
>>
>>
>> As attachment you can find an archive containing the matrix, the initial
>> solution (used to manually build the null space) and the rhs (not used in
>> the test code) in binary format.
>>
>> You can also find the testing code in the same archive.
>>
>> I used petsc 3.12(gcc+openMPI) and petsc 3.15.2(intelOneAPI) same results.
>>
>> If the attachment is not delivered, I can share a link to it.
>>
>>
>>
>> Marco, please forgive me for taking so long to get to your issue. It has
>> been crazy.
>>
>>
>>
>> You are correct, we had a bug. it is in MatNullSpaceTest. The
>> normalization for the constant vector was wrong:
>>
>>
>>
>> diff --git a/src/mat/interface/matnull.c b/src/mat/interface/matnull.c
>> index f8ab2925988..0c4c3855be0 100644
>> --- a/src/mat/interface/matnull.c
>> +++ b/src/mat/interface/matnull.c
>> @@ -429,7 +429,7 @@ PetscErrorCode  MatNullSpaceTest(MatNullSpace sp,Mat
>> mat,PetscBool  *isNull)
>>    if (sp->has_cnst) {
>>      ierr = VecDuplicate(l,&r);CHKERRQ(ierr);
>>      ierr = VecGetSize(l,&N);CHKERRQ(ierr);
>> -    sum  = 1.0/N;
>>
>> +    sum  = 1.0/PetscSqrtReal(N);
>>      ierr = VecSet(l,sum);CHKERRQ(ierr);
>>      ierr = MatMult(mat,l,r);CHKERRQ(ierr);
>>      ierr = VecNorm(r,NORM_2,&nrm);CHKERRQ(ierr);
>>
>>
>>
>> With this fix, your two cases give the same answer, namely that the
>> constant vector is not a null vector of your
>>
>> operator, but it is close, as your can see
>> using -mat_null_space_test_view.
>>
>>
>>
>>   Thanks,
>>
>>
>>
>>      Matt
>>
>>
>>
>> Thanks for any help.
>>
>>
>>
>> Marco Cisternino
>>
>>
>>
>>
>>
>> Marco Cisternino, PhD
>> marco.cisternino at optimad.it
>>
>> ______________________
>>
>> Optimad Engineering Srl
>>
>> Via Bligny 5, Torino, Italia.
>> +3901119719782
>> www.optimad.it
>>
>>
>>
>> *From:* Marco Cisternino <marco.cisternino at optimad.it>
>> *Sent:* martedì 7 dicembre 2021 19:36
>> *To:* Matthew Knepley <knepley at gmail.com>
>> *Cc:* petsc-users <petsc-users at mcs.anl.gov>
>> *Subject:* Re: [petsc-users] Nullspaces
>>
>>
>>
>> I will, as soon as possible...
>>
>>
>>
>> Scarica Outlook per Android <https://aka.ms/AAb9ysg>
>> ------------------------------
>>
>> *From:* Matthew Knepley <knepley at gmail.com>
>> *Sent:* Tuesday, December 7, 2021 7:25:43 PM
>> *To:* Marco Cisternino <marco.cisternino at optimad.it>
>> *Cc:* petsc-users <petsc-users at mcs.anl.gov>
>> *Subject:* Re: [petsc-users] Nullspaces
>>
>>
>>
>> On Tue, Dec 7, 2021 at 11:19 AM Marco Cisternino <
>> marco.cisternino at optimad.it> wrote:
>>
>> Good morning,
>>
>> I’m still struggling with the Poisson equation with Neumann BCs.
>>
>> I discretize the equation by finite volume method and I divide every line
>> of the linear system by the volume of the cell. I could avoid this
>> division, but I’m trying to understand.
>>
>> My mesh is not uniform, i.e. cells have different volumes (it is an
>> octree mesh).
>>
>> Moreover, in my computational domain there are 2 separated sub-domains.
>>
>> I build the null space and then I use MatNullSpaceTest to check it.
>>
>>
>>
>> If I do this:
>>
>> MatNullSpaceCreate(getCommunicator(), PETSC_TRUE, 0, nullptr, &nullspace);
>>
>> It works
>>
>>
>>
>> This produces the normalized constant vector.
>>
>>
>>
>> If I do this:
>>
>> Vec nsp;
>>
>> VecDuplicate(m_rhs, &nsp);
>>
>> VecSet(nsp,1.0);
>>
>> VecNormalize(nsp, nullptr);
>>
>> MatNullSpaceCreate(getCommunicator(), PETSC_FALSE, 1, &nsp, &nullspace);
>>
>> It does not work
>>
>>
>>
>> This is also the normalized constant vector.
>>
>>
>>
>> So you are saying that these two vectors give different results with
>> MatNullSpaceTest()?
>>
>> Something must be wrong in the code. Can you send a minimal example of
>> this? I will go
>>
>> through and debug it.
>>
>>
>>
>>   Thanks,
>>
>>
>>
>>      Matt
>>
>>
>>
>> Probably, I have wrong expectations, but should not it be the same?
>>
>>
>>
>> Thanks
>>
>>
>>
>> Marco Cisternino, PhD
>> marco.cisternino at optimad.it
>>
>> ______________________
>>
>> Optimad Engineering Srl
>>
>> Via Bligny 5, Torino, Italia.
>> +3901119719782
>> www.optimad.it
>>
>>
>>
>>
>>
>>
>> --
>>
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which their
>> experiments lead.
>> -- Norbert Wiener
>>
>>
>>
>> https://www.cse.buffalo.edu/~knepley/
>> <http://www.cse.buffalo.edu/~knepley/>
>>
>>
>>
>>
>> --
>>
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which their
>> experiments lead.
>> -- Norbert Wiener
>>
>>
>>
>> https://www.cse.buffalo.edu/~knepley/
>> <http://www.cse.buffalo.edu/~knepley/>
>>
>>
>>
>>
>> --
>>
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which their
>> experiments lead.
>> -- Norbert Wiener
>>
>>
>>
>> https://www.cse.buffalo.edu/~knepley/
>> <http://www.cse.buffalo.edu/~knepley/>
>>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> <http://www.cse.buffalo.edu/~knepley/>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220119/744ecb00/attachment-0001.html>