Skip to content

Commit 05df909

Browse files
committed
Set nullspace for both J and P operators
These can be different for instance when doing PJFNK. This now matters post https://gitlab.com/petsc/petsc/-/merge_requests/8247
1 parent 1fba560 commit 05df909

1 file changed

Lines changed: 46 additions & 35 deletions

File tree

src/solvers/petsc_nonlinear_solver.C

Lines changed: 46 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -1008,39 +1008,6 @@ PetscNonlinearSolver<T>::solve (SparseMatrix<T> & pre_in, // System Preconditi
10081008
if (this->jacobian || this->jacobian_object || this->residual_and_jacobian_object)
10091009
LibmeshPetscCall(SNESSetJacobian (_snes, pre->mat(), pre->mat(), libmesh_petsc_snes_jacobian, this));
10101010

1011-
1012-
// Only set the nullspace if we have a way of computing it and the result is non-empty.
1013-
if (this->nullspace || this->nullspace_object)
1014-
{
1015-
WrappedPetsc<MatNullSpace> msp;
1016-
this->build_mat_null_space(this->nullspace_object, this->nullspace, msp.get());
1017-
if (msp)
1018-
LibmeshPetscCall(MatSetNullSpace(pre->mat(), msp));
1019-
}
1020-
1021-
// Only set the transpose nullspace if we have a way of computing it and the result is non-empty.
1022-
if (this->transpose_nullspace || this->transpose_nullspace_object)
1023-
{
1024-
#if PETSC_VERSION_LESS_THAN(3,6,0)
1025-
libmesh_warning("MatSetTransposeNullSpace is only supported for PETSc >= 3.6, transpose nullspace will be ignored.");
1026-
#else
1027-
WrappedPetsc<MatNullSpace> msp;
1028-
this->build_mat_null_space(this->transpose_nullspace_object, this->transpose_nullspace, msp.get());
1029-
if (msp)
1030-
LibmeshPetscCall(MatSetTransposeNullSpace(pre->mat(), msp));
1031-
#endif
1032-
}
1033-
1034-
// Only set the nearnullspace if we have a way of computing it and the result is non-empty.
1035-
if (this->nearnullspace || this->nearnullspace_object)
1036-
{
1037-
WrappedPetsc<MatNullSpace> msp;
1038-
this->build_mat_null_space(this->nearnullspace_object, this->nearnullspace, msp.get());
1039-
1040-
if (msp)
1041-
LibmeshPetscCall(MatSetNearNullSpace(pre->mat(), msp));
1042-
}
1043-
10441011
// Have the Krylov subspace method use our good initial guess rather than 0
10451012
KSP ksp;
10461013
LibmeshPetscCall(SNESGetKSP (_snes, &ksp));
@@ -1101,8 +1068,8 @@ PetscNonlinearSolver<T>::solve (SparseMatrix<T> & pre_in, // System Preconditi
11011068
#endif
11021069
LibmeshPetscCall(SNESSetUp(_snes));
11031070

1104-
Mat J;
1105-
LibmeshPetscCall(SNESGetJacobian(_snes, &J, LIBMESH_PETSC_NULLPTR,
1071+
Mat J, P;
1072+
LibmeshPetscCall(SNESGetJacobian(_snes, &J, &P,
11061073
LIBMESH_PETSC_NULLPTR,
11071074
LIBMESH_PETSC_NULLPTR));
11081075
LibmeshPetscCall(MatMFFDSetFunction(J, libmesh_petsc_snes_mffd_interface, this));
@@ -1120,6 +1087,50 @@ PetscNonlinearSolver<T>::solve (SparseMatrix<T> & pre_in, // System Preconditi
11201087
#endif
11211088
#endif
11221089

1090+
// Only set the nullspace if we have a way of computing it and the result is non-empty.
1091+
if (this->nullspace || this->nullspace_object)
1092+
{
1093+
WrappedPetsc<MatNullSpace> msp;
1094+
this->build_mat_null_space(this->nullspace_object, this->nullspace, msp.get());
1095+
if (msp)
1096+
{
1097+
LibmeshPetscCall(MatSetNullSpace(J, msp));
1098+
if (P != J)
1099+
LibmeshPetscCall(MatSetNullSpace(P, msp));
1100+
}
1101+
}
1102+
1103+
// Only set the transpose nullspace if we have a way of computing it and the result is non-empty.
1104+
if (this->transpose_nullspace || this->transpose_nullspace_object)
1105+
{
1106+
#if PETSC_VERSION_LESS_THAN(3,6,0)
1107+
libmesh_warning("MatSetTransposeNullSpace is only supported for PETSc >= 3.6, transpose nullspace will be ignored.");
1108+
#else
1109+
WrappedPetsc<MatNullSpace> msp;
1110+
this->build_mat_null_space(this->transpose_nullspace_object, this->transpose_nullspace, msp.get());
1111+
if (msp)
1112+
{
1113+
LibmeshPetscCall(MatSetTransposeNullSpace(J, msp));
1114+
if (P != J)
1115+
LibmeshPetscCall(MatSetTransposeNullSpace(P, msp));
1116+
}
1117+
#endif
1118+
}
1119+
1120+
// Only set the nearnullspace if we have a way of computing it and the result is non-empty.
1121+
if (this->nearnullspace || this->nearnullspace_object)
1122+
{
1123+
WrappedPetsc<MatNullSpace> msp;
1124+
this->build_mat_null_space(this->nearnullspace_object, this->nearnullspace, msp.get());
1125+
1126+
if (msp)
1127+
{
1128+
LibmeshPetscCall(MatSetNearNullSpace(J, msp));
1129+
if (P != J)
1130+
LibmeshPetscCall(MatSetNearNullSpace(P, msp));
1131+
}
1132+
}
1133+
11231134
SNESLineSearch linesearch;
11241135
if (linesearch_object)
11251136
{

0 commit comments

Comments
 (0)