Actual source code: lostnullspace.c
1: static char help[] = "Losing nullspaces in PCFIELDSPLIT after zeroing rows.\n";
3: // Contributed by Jeremy Theler
5: #include <petscksp.h>
7: int main(int argc, char **args)
8: {
9: KSP ksp, *sub_ksp;
10: Vec x, b, rigid_mode[6];
11: PetscViewer viewer;
12: PetscInt rows, cols, size, bs, n_splits = 0;
13: PetscBool has_columns = PETSC_FALSE;
14: Mat A, K;
15: MatNullSpace nullsp, near_null_space;
16: IS is_thermal, is_mech;
17: PC pc, pc_thermal, pc_mech;
18: const PetscInt *bc_thermal_indexes, *bc_mech_indexes;
19: char datafilespath[PETSC_MAX_PATH_LEN], datafile[PETSC_MAX_PATH_LEN];
21: PetscFunctionBeginUser;
22: PetscCall(PetscInitialize(&argc, &args, NULL, help));
23: PetscCall(PetscOptionsGetString(NULL, NULL, "-datafilespath", datafilespath, sizeof(datafilespath), NULL));
25: PetscCall(PetscStrcpy(datafile, datafilespath));
26: PetscCall(PetscStrcat(datafile, "/lostnullspace/"));
27: PetscCall(PetscStrcat(datafile, "A.bin"));
28: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, datafile, FILE_MODE_READ, &viewer));
29: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
30: PetscCall(MatLoad(A, viewer));
31: PetscCall(PetscViewerDestroy(&viewer));
33: PetscCall(MatGetNearNullSpace(A, &nullsp));
34: PetscCall(MatGetSize(A, &rows, &cols));
35: PetscCall(MatGetBlockSize(A, &bs));
36: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "A has rows = %" PetscInt_FMT ", cols = %" PetscInt_FMT ", bs = %" PetscInt_FMT ", nearnullsp = %p\n", rows, cols, bs, nullsp));
38: PetscCall(PetscStrcpy(datafile, datafilespath));
39: PetscCall(PetscStrcat(datafile, "/lostnullspace/"));
40: PetscCall(PetscStrcat(datafile, "is.bin"));
41: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, datafile, FILE_MODE_READ, &viewer));
43: PetscCall(ISCreate(PETSC_COMM_WORLD, &is_thermal));
44: PetscCall(ISLoad(is_thermal, viewer));
45: PetscCall(ISGetSize(is_thermal, &size));
46: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "thermal field size = %" PetscInt_FMT " \n", size));
48: PetscCall(ISCreate(PETSC_COMM_WORLD, &is_mech));
49: PetscCall(ISLoad(is_mech, viewer));
50: PetscCall(ISGetSize(is_mech, &size));
51: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mechanical field size = %" PetscInt_FMT " \n", size));
52: PetscCall(PetscViewerDestroy(&viewer));
54: PetscCall(MatCreateVecs(A, &x, &b));
55: PetscCall(VecZeroEntries(x));
56: PetscCall(VecZeroEntries(b));
58: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
59: PetscCall(KSPSetOperators(ksp, A, A));
61: PetscCall(KSPSetType(ksp, KSPPREONLY));
63: PetscCall(KSPGetPC(ksp, &pc));
64: PetscCall(PCSetType(pc, PCFIELDSPLIT));
65: PetscCall(PCFieldSplitSetIS(pc, "thermal", is_thermal));
66: PetscCall(PCFieldSplitSetIS(pc, "mechanical", is_mech));
67: PetscCall(PCSetUp(pc));
68: PetscCall(PCFieldSplitGetSubKSP(pc, &n_splits, &sub_ksp));
69: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "n_splits = %" PetscInt_FMT " \n", n_splits));
70: PetscCall(KSPSetType(sub_ksp[0], KSPGMRES));
72: PetscCall(KSPGetPC(sub_ksp[0], &pc_thermal));
73: PetscCall(PCSetType(pc_thermal, PCJACOBI));
74: PetscCall(KSPSetFromOptions(sub_ksp[0]));
76: PetscCall(KSPSetType(sub_ksp[1], KSPCG));
77: PetscCall(KSPGetPC(sub_ksp[1], &pc_mech));
78: PetscCall(PCSetType(pc_mech, PCGAMG));
79: PetscCall(KSPSetFromOptions(sub_ksp[1]));
81: PetscCall(PetscStrcpy(datafile, datafilespath));
82: PetscCall(PetscStrcat(datafile, "/lostnullspace/"));
83: PetscCall(PetscStrcat(datafile, "rigid-modes.bin"));
84: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, datafile, FILE_MODE_READ, &viewer));
85: for (PetscInt i = 0; i < 6; i++) {
86: PetscCall(VecCreate(PETSC_COMM_WORLD, &rigid_mode[i]));
87: PetscCall(VecLoad(rigid_mode[i], viewer));
88: }
89: PetscCall(PetscViewerDestroy(&viewer));
91: PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_FALSE, 6, rigid_mode, &near_null_space));
92: PetscCall(KSPGetOperators(sub_ksp[1], &K, PETSC_NULLPTR));
93: PetscCall(MatSetNearNullSpace(K, near_null_space));
94: PetscCall(MatSetBlockSize(K, 3));
96: PetscCall(MatGetSize(K, &rows, &cols));
97: PetscCall(MatGetBlockSize(K, &bs));
98: PetscCall(MatGetNearNullSpace(K, &nullsp));
99: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "K has rows = %" PetscInt_FMT ", cols = %" PetscInt_FMT ", bs = %" PetscInt_FMT ", nearnullsp = %p\n", rows, cols, bs, nullsp));
101: PetscCall(ISGetIndices(is_thermal, &bc_thermal_indexes));
102: PetscCall(ISGetIndices(is_mech, &bc_mech_indexes));
104: // check that MatZeroRows() without MAT_KEEP_NONZERO_PATTERN does not remove the near null spaces attached to the submatrices
105: PetscCall(PetscOptionsHasName(PETSC_NULLPTR, PETSC_NULLPTR, "-columns", &has_columns));
106: if (has_columns == PETSC_TRUE) {
107: PetscCall(MatZeroRowsColumns(A, 3, bc_mech_indexes, 1, PETSC_NULLPTR, PETSC_NULLPTR));
108: PetscCall(MatZeroRowsColumns(A, 1, bc_thermal_indexes, 1, PETSC_NULLPTR, PETSC_NULLPTR));
109: } else {
110: PetscCall(MatZeroRows(A, 3, bc_mech_indexes, 1, PETSC_NULLPTR, PETSC_NULLPTR));
111: PetscCall(MatZeroRows(A, 1, bc_thermal_indexes, 1, PETSC_NULLPTR, PETSC_NULLPTR));
112: }
113: PetscCall(ISRestoreIndices(is_mech, &bc_mech_indexes));
114: PetscCall(ISRestoreIndices(is_thermal, &bc_thermal_indexes));
116: PetscCall(KSPSetFromOptions(ksp));
117: PetscCall(KSPSetUp(ksp));
118: PetscCall(KSPSolve(ksp, b, x));
120: PetscCall(MatDestroy(&A));
121: PetscCall(KSPDestroy(&ksp));
122: PetscCall(ISDestroy(&is_mech));
123: PetscCall(ISDestroy(&is_thermal));
124: PetscCall(VecDestroy(&x));
125: PetscCall(VecDestroy(&b));
126: PetscCall(MatNullSpaceDestroy(&near_null_space));
127: for (PetscInt i = 0; i < 6; i++) PetscCall(VecDestroy(&rigid_mode[i]));
128: PetscCall(PetscFree(sub_ksp));
129: PetscCall(PetscFinalize());
130: return 0;
131: }
133: /*TEST
135: test:
136: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
137: args: -datafilespath ${DATAFILESPATH} -ksp_view
138: filter: grep "near null"
140: TEST*/