[petsc-users] Nullspaces
Matthew Knepley
knepley at gmail.com
Tue Jan 18 14:24:35 CST 2022
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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220118/ddec39a9/attachment-0001.html>
More information about the petsc-users
mailing list