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*/