[petsc-users] Nullspaces

Tue Jan 18 14:20:04 CST 2022

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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220118/48b045df/attachment.html>
