bug in nullspace removal?

[amneetb]

I think there is a bug in removing the nullspace component in the solution vector here. Specifically, we are calculating L2_norm of the null vector, whereas it should be L2_norm^2

https://github.com/IBAMR/IBAMR/blob/master/src/navier_stokes/INSStaggeredHierarchyIntegrator.cpp#L1592

Vectors in d_nul_vecs are not normalized. Thus, in the final result there should be two L2_norm factors in the denominator

Remove component of vector $\vec{b}$ from $\vec{a}$

$\vec{a} = \vec{a} - \left(\frac{\vec{a} \cdot \vec{b}}{|| b||}\right) \frac{\vec{b}}{|| b ||} $

If we take a dot product of the above equation with $\vec{b}$ or $\vec{b}/|| b||$, we will get zero.

The fix is just to remove the std::sqrt on line 1592.

[boyceg]

Why should it be the norm squared?

[amneetb]

My comment explains it :-)

[amneetb]

Basically it is a Gram-Schmidt process. We work with orthonormal vectors q. The null vectors are not orthonormal. We need two normalization factors: first when taking a dot product, and second when using the null vector in the AXPY process.

[boyceg]

Ah, I missed that this was not just normalizing the vector. (I should have waited until I got off my phone to check. 🙄)

Could / should the nullspace basis vectors be stored as normalized vectors in the first place?

[amneetb]

Yes, that would be the right thing to do. We also pass the vectors to other libraries like petsc. Maybe petsc does normalization before using them. If not then it is also not using the correct vectors (our fault not petsc’s).

--Amneet

On Sun, Jan 21, 2024 at 8:01 AM Boyce Griffith @.***> wrote:

Ah, I missed that this was not just normalizing the vector.

Could / should the nullspace basis vectors be stored as normalized vectors in the first place?

— Reply to this email directly, view it on GitHub https://github.com/IBAMR/IBAMR/issues/1673#issuecomment-1902679558, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABYJJ2BSR5AQCERNSO6ECF3YPU3UJAVCNFSM6AAAAABCECKMP6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSMBSGY3TSNJVHA . You are receiving this because you authored the thread.Message ID: @.***>

[boyceg]

Agree that the current code is wrong. It should normalize by nul_dot_nul. But I think the better fix actually is to store the nullspace vectors as orthonormal vectors using modified Gram-Schmidt.

BTW, MatNullSpaceCreate does indeed request orthonormal nullspace vectors — which we are not currently providing.

[boyceg]

It looks like the only place we set up a PETSc nullspace is in PETScKrylovLinearSolver, and that normalizes the vectors but doesn't orthogonalize them.