Re: [petsc-users] Nullspaces

[Marco Cisternino]

Hello Stefano and thank you for your support.
I never used SLEPc before but I did this:
right after the matrix loading from file I added the following lines to my 
shared tiny code
    MatLoad(matrix, v);

    EPS eps;
    EPSCreate(PETSC_COMM_WORLD, &eps);
    EPSSetOperators( eps, matrix, NULL );
    EPSSetWhichEigenpairs(eps, EPS_SMALLEST_MAGNITUDE);
    EPSSetProblemType( eps, EPS_NHEP );
    EPSSetTolerances(eps, 1.0e-12, 1000);
    EPSSetFromOptions( eps );
    EPSSolve( eps );

    Vec xr, xi; /* eigenvector, x */
    PetscScalar kr, ki; /* eigenvalue, k */
    PetscInt j, nconv;
    PetscReal error;
    EPSGetConverged( eps, &nconv );
    for (j=0; j<nconv; j++) {
        EPSGetEigenpair( eps, j, &kr, &ki, xr, xi );
        EPSComputeError( eps, j, EPS_ERROR_ABSOLUTE, &error );
        std::cout << "The smallest eigenvalue is = (" << kr << ", " << ki << ") 
with error = " << error << std::endl;
    }

I launched using
mpirun -n 1 ./testnullspace -eps_monitor

and the output is

  1 EPS nconv=0 first unconverged value (error) -1499.29 (6.57994794e+01)
  2 EPS nconv=0 first unconverged value (error) -647.468 (5.39939262e+01)
  3 EPS nconv=0 first unconverged value (error) -177.157 (9.49337698e+01)
  4 EPS nconv=0 first unconverged value (error) 59.6771 (1.62531943e+02)
  5 EPS nconv=0 first unconverged value (error) 41.755 (1.41965990e+02)
  6 EPS nconv=0 first unconverged value (error) -11.5462 (3.60453662e+02)
  7 EPS nconv=0 first unconverged value (error) -6.04493 (4.60890030e+02)
  8 EPS nconv=0 first unconverged value (error) -22.7362 (8.67630086e+01)
  9 EPS nconv=0 first unconverged value (error) -12.9637 (1.08507821e+02)
10 EPS nconv=0 first unconverged value (error) 7.7234 (1.53561979e+02)
…
111 EPS nconv=0 first unconverged value (error) -2.27e-08 (6.84762319e+00)
112 EPS nconv=0 first unconverged value (error) -2.60619e-08 (4.45245528e+00)
113 EPS nconv=0 first unconverged value (error) -5.49592e-09 (1.87798984e+01)
114 EPS nconv=0 first unconverged value (error) -9.9456e-09 (7.96711076e+00)
115 EPS nconv=0 first unconverged value (error) -1.89779e-08 (4.15471472e+00)
116 EPS nconv=0 first unconverged value (error) -2.05288e-08 (2.52953194e+00)
117 EPS nconv=0 first unconverged value (error) -2.02919e-08 (2.90090711e+00)
118 EPS nconv=0 first unconverged value (error) -3.8706e-08 (8.03595736e-01)
119 EPS nconv=1 first unconverged value (error) -61751.8 (9.58036571e-07)
Computed 1 pairs
The smallest eigenvalue is = (-3.8706e-08, 0) with error = 4.9707e-07

Am I using SLEPc in the right way at least for the first smallest eigenvalue? 
If I’m on the right way I can find out how to compute the second one.

Thanks a lot

Marco Cisternino

From: Stefano Zampini <stefano.zamp...@gmail.com>
Sent: martedì 4 gennaio 2022 17:39
To: Marco Cisternino <marco.cistern...@optimad.it>
Cc: Mark Adams <mfad...@lbl.gov>; petsc-users <petsc-users@mcs.anl.gov>
Subject: Re: [petsc-users] Nullspaces


Il Mar 4 Gen 2022, 16:56 Marco Cisternino 
<marco.cistern...@optimad.it<mailto:marco.cistern...@optimad.it>> ha scritto:
Hello Mark,
I analyzed the codes with valgrind, both the real code and the tiny one.
I obviously used memcheck tool but with full leak check compiling the codes 
with debug info.
Not considering OpenMPI events (I have no wrappers on the machine I used for 
the analysis), the real code gave zero errors and the tiny one gave this
==17911== 905,536 (1,552 direct, 903,984 indirect) bytes in 1 blocks are 
definitely lost in loss record 33 of 33
==17911==    at 0x483E340: memalign (in 
/usr/lib/x86_64-linux-gnu/valgrind/vgpreload_memcheck-amd64-linux.so)
==17911==    by 0x49CB672: PetscMallocAlign (in 
/usr/lib/x86_64-linux-gnu/libpetsc_real.so.3.12.4)
==17911==    by 0x49CBE1D: PetscMallocA (in 
/usr/lib/x86_64-linux-gnu/libpetsc_real.so.3.12.4)
==17911==    by 0x4B26187: VecCreate (in 
/usr/lib/x86_64-linux-gnu/libpetsc_real.so.3.12.4)
==17911==    by 0x10940D: main (testNullSpace.cpp:30)

due to the fact that I forgot to destroy the solution Vec (adding 
VecDestroy(&solution) at the end of the main, the error disappear).
For both the codes, I analyzed the two ways of passing the constant to the null 
space of the operator, no memory errors but still the same results from 
MatNullSpaceTest, i.e.

MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, nullptr, &nullspace);

passes the test while

Vec* nsp;
VecDuplicateVecs(solution, 1, &nsp);
VecSet(nsp[0],1.0);
VecNormalize(nsp[0], nullptr);
MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_FALSE, 1, nsp, &nullspace);
VecDestroyVecs(1,&nsp);
PetscFree(nsp);

does not.

I hope this can satisfy your doubt about the memory behavior, but please do not 
hesitate to ask for more analysis if it cannot.

As Matthew said some weeks ago, something should be wrong in the code, I would 
say in the matrix, that’s why I provided the matrix and the way I test it.
Unfortunately, it is hard (read impossible) for me to share the code producing 
the matrix. I hope the minimal code I provided is enough to understand 
something.

Try running slepc to find the smallest eigenvalues.  There should be two zero 
eigenvalues,  then inspect the eigenvectors

Thank you all.

Marco Cisternino

From: Marco Cisternino
Sent: lunedì 3 gennaio 2022 16:08
To: Mark Adams <mfad...@lbl.gov<mailto:mfad...@lbl.gov>>
Cc: Matthew Knepley <knep...@gmail.com<mailto:knep...@gmail.com>>; petsc-users 
<petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>>
Subject: RE: [petsc-users] Nullspaces

We usually analyze the code with valgrind, when important changes are 
implemented.
I have to admit that this analysis is still not automatic and the case we are 
talking about is not a test case for our workload.
The test cases we have give no errors in valgrind analysis.

However, I will analyze both the real code and the tiny one for this case with 
valgrind and report the results.

Thank you,


Marco Cisternino

From: Mark Adams <mfad...@lbl.gov<mailto:mfad...@lbl.gov>>
Sent: lunedì 3 gennaio 2022 15:50
To: Marco Cisternino 
<marco.cistern...@optimad.it<mailto:marco.cistern...@optimad.it>>
Cc: Matthew Knepley <knep...@gmail.com<mailto:knep...@gmail.com>>; petsc-users 
<petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] Nullspaces

I have not looked at your code, but as a general observation you want to have 
some sort of memory checker, like valgrid for CPUs, in your workflow.
It is the fastest way to find some classes of bugs.

On Mon, Jan 3, 2022 at 8:47 AM Marco Cisternino 
<marco.cistern...@optimad.it<mailto:marco.cistern...@optimad.it>> wrote:
Are you talking about the code that produce the linear system or about the tiny 
code that test the null space?
In the first case, it is absolutely possible, but I would expect no problem in 
the tiny code, do you agree?
It is important to remark that the real code and the tiny one behave in the 
same way when testing the null space of the operator. I can analyze with 
valgrind and I will, but I would not expect great insights.

Thanks,

Marco Cisternino, PhD
marco.cistern...@optimad.it<mailto:marco.cistern...@optimad.it>
______________________
Optimad Engineering Srl
Via Bligny 5, Torino, Italia.
+3901119719782
www.optimad.it<http://www.optimad.it/>

From: Mark Adams <mfad...@lbl.gov<mailto:mfad...@lbl.gov>>
Sent: lunedì 3 gennaio 2022 14:42
To: Marco Cisternino 
<marco.cistern...@optimad.it<mailto:marco.cistern...@optimad.it>>
Cc: Matthew Knepley <knep...@gmail.com<mailto:knep...@gmail.com>>; petsc-users 
<petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] Nullspaces

There could be a memory bug that does not cause a noticeable problem until it 
hits some vital data and valgrind might find it on a small problem.

However you might have a bug like a hardwired buffer size that overflows that 
is in fact not a bug until you get to this large size and in that case valgrid 
would need to be run on the large case and would have a good chance of finding 
it.


On Mon, Jan 3, 2022 at 4:42 AM Marco Cisternino 
<marco.cistern...@optimad.it<mailto:marco.cistern...@optimad.it>> wrote:
My comments are between the Mark’s lines and they starts with “#”

Marco Cisternino

From: Mark Adams <mfad...@lbl.gov<mailto:mfad...@lbl.gov>>
Sent: sabato 25 dicembre 2021 14:59
To: Marco Cisternino 
<marco.cistern...@optimad.it<mailto:marco.cistern...@optimad.it>>
Cc: Matthew Knepley <knep...@gmail.com<mailto:knep...@gmail.com>>; petsc-users 
<petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] Nullspaces

If  "triggering the issue" requires a substantial mesh, that makes me think 
there is a logic bug somewhere. Maybe use valgrind.

# Are you suggesting to use valgrind on this tiny toy code or on the original 
one? However, considering the purpose of the tiny code, i.e. testing the 
constant null space, why there should be a logical bug? Case 1 passes and case 
2 should be exactly the same, shouldn’t be it?

Also you say you divide by the cell volume. Maybe I am not understanding this 
but that is basically diagonal scaling and that will change the null space (ie, 
not a constant anymore)

# I agree on this, but it pushes a question: why the case 1 passes the test?
# Thank you, Mark.

On Thu, Dec 16, 2021 at 11:11 AM Marco Cisternino 
<marco.cistern...@optimad.it<mailto:marco.cistern...@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.

Thanks for any help.

Marco Cisternino


Marco Cisternino, PhD
marco.cistern...@optimad.it<mailto:marco.cistern...@optimad.it>
______________________
Optimad Engineering Srl
Via Bligny 5, Torino, Italia.
+3901119719782
www.optimad.it<http://www.optimad.it/>

From: Marco Cisternino 
<marco.cistern...@optimad.it<mailto:marco.cistern...@optimad.it>>
Sent: martedì 7 dicembre 2021 19:36
To: Matthew Knepley <knep...@gmail.com<mailto:knep...@gmail.com>>
Cc: petsc-users <petsc-users@mcs.anl.gov<mailto:petsc-users@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 <knep...@gmail.com<mailto:knep...@gmail.com>>
Sent: Tuesday, December 7, 2021 7:25:43 PM
To: Marco Cisternino 
<marco.cistern...@optimad.it<mailto:marco.cistern...@optimad.it>>
Cc: petsc-users <petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] Nullspaces

On Tue, Dec 7, 2021 at 11:19 AM Marco Cisternino 
<marco.cistern...@optimad.it<mailto:marco.cistern...@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.cistern...@optimad.it<mailto:marco.cistern...@optimad.it>

______________________

Optimad Engineering Srl

Via Bligny 5, Torino, Italia.
+3901119719782
www.optimad.it<http://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/>