Thanks Barry.
I already did it, I mean a constant vec per domain with 1 on rows of that
domain and zero on rows of the other one, and exactly I did this:
Vec *constants = nullptr;
PetscInt nConstants = nActiveDomains;
VecDuplicateVecs(m_solution, nConstants, &constants);
for (PetscInt i = 0; i < nConstants; ++i) {
VecSet(constants[i],0.0);
}
std::vector<PetscScalar *> rawConstants(nConstants);
for (PetscInt i = 0; i < nConstants; ++i) {
VecGetArray(constants[i], &rawConstants[i]);
}
long nRows = getRowCount();
const CBaseLinearSolverMapping *cellRowMapping = getMapping();
for (long cellRow = 0; cellRow < nRows; ++cellRow) {
std::size_t cellRawId = cellRowMapping->getRowRawCell(cellRow);
int cellDomain =
m_grid->getCellStorage().getDomains().rawAt(cellRawId);
if (m_grid->isCellDomainActive(cellDomain)) {
std::size_t activeDomainId = activeDomainIds[cellDomain];
rawConstants[activeDomainId][cellRow] = 1.;
}
}
for (PetscInt i = 0; i < nConstants; ++i) {
VecRestoreArray(constants[i], &rawConstants[i]);
VecAssemblyBegin(constants[i]);
VecAssemblyEnd(constants[i]);
VecNormalize(constants[i], nullptr);
}
MatNullSpaceCreate(getCommunicator(), PETSC_FALSE, nConstants,
constants, &nullspace);
But this does not pass the test if I divide by cell volume every row of the
linear system.
It works if I do not perform the division
Thanks
Marco Cisternino, PhD
[email protected]<mailto:[email protected]>
______________________
Optimad Engineering Srl
Via Bligny 5, Torino, Italia.
+3901119719782
www.optimad.it<http://www.optimad.it/>
From: Barry Smith <[email protected]>
Sent: martedì 7 dicembre 2021 17:53
To: Marco Cisternino <[email protected]>
Cc: petsc-users <[email protected]>
Subject: Re: [petsc-users] Nullspaces
A side note: The MatNullSpaceTest tells you that the null space you provided
is in the null space of the operator, it does not say if you have found the
entire null space. In your case with two subdomains the null space is actually
two dimensional; all constant on one domain (0 on the other) and 0 on the first
domain and all constant on the second. So you need to pass two vectors into the
MatNullSpaceCreate(). But regardless the test should work for all constant on
both domains.
On Dec 7, 2021, at 11:19 AM, Marco Cisternino
<[email protected]<mailto:[email protected]>> 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
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
Probably, I have wrong expectations, but should not it be the same?
Thanks
Marco Cisternino, PhD
[email protected]<mailto:[email protected]>
______________________
Optimad Engineering Srl
Via Bligny 5, Torino, Italia.
+3901119719782
www.optimad.it<http://www.optimad.it/>
