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 <knep...@gmail.com> wrote: > On Thu, Dec 16, 2021 at 11:09 AM Marco Cisternino < > 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. >> > > 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.cistern...@optimad.it >> >> ______________________ >> >> Optimad Engineering Srl >> >> Via Bligny 5, Torino, Italia. >> +3901119719782 >> www.optimad.it >> >> >> >> *From:* Marco Cisternino <marco.cistern...@optimad.it> >> *Sent:* martedì 7 dicembre 2021 19:36 >> *To:* Matthew Knepley <knep...@gmail.com> >> *Cc:* petsc-users <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> >> *Sent:* Tuesday, December 7, 2021 7:25:43 PM >> *To:* Marco Cisternino <marco.cistern...@optimad.it> >> *Cc:* petsc-users <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> 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 >> >> ______________________ >> >> 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/>