Actual source code: ex84.c
1: static const char help[] = "Solves a Q2-Q1 Navier-Stokes problem.\n\n";
3: #include <petscksp.h>
4: #include <petscdmda.h>
6: PetscErrorCode LSCLoadOperators(Mat *A, Mat *Q, Mat *L, Vec *rhs, IS *velocity, IS *pressure)
7: {
8: PetscViewer viewer;
9: char filename[PETSC_MAX_PATH_LEN];
10: PetscBool flg;
12: PetscFunctionBeginUser;
13: PetscCall(MatCreate(PETSC_COMM_WORLD, A));
14: PetscCall(MatCreate(PETSC_COMM_WORLD, Q));
15: if (L) PetscCall(MatCreate(PETSC_COMM_WORLD, L));
16: PetscCall(ISCreate(PETSC_COMM_WORLD, velocity));
17: PetscCall(ISCreate(PETSC_COMM_WORLD, pressure));
18: PetscCall(VecCreate(PETSC_COMM_WORLD, rhs));
19: /* Load matrices from a Q2-Q1 discretisation of Navier-Stokes. The data is packed into one file. */
20: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg));
21: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must provide a data file with -f");
22: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &viewer));
23: PetscCall(MatLoad(*A, viewer));
24: PetscCall(VecLoad(*rhs, viewer));
25: PetscCall(ISLoad(*velocity, viewer));
26: PetscCall(ISLoad(*pressure, viewer));
27: PetscCall(MatLoad(*Q, viewer));
28: if (L) PetscCall(MatLoad(*L, viewer));
29: PetscCall(PetscViewerDestroy(&viewer));
30: PetscFunctionReturn(PETSC_SUCCESS);
31: }
33: PetscErrorCode port_lsd_bfbt(void)
34: {
35: Mat A, Q, L = NULL;
36: Vec x, b;
37: KSP ksp_A;
38: PC pc_A;
39: IS isu, isp;
40: PetscBool commute_lsc = PETSC_FALSE;
41: KSP *subksp; // This will be length two, with the former being the A KSP and the latter being the
42: // Schur complement KSP
43: KSP schur_complement_ksp;
44: PC lsc_pc;
45: PetscInt num_splits;
46: Mat lsc_pc_pmat;
48: PetscFunctionBeginUser;
49: PetscCall(PetscOptionsGetBool(NULL, NULL, "-commute_lsc", &commute_lsc, NULL));
50: if (commute_lsc) PetscCall(LSCLoadOperators(&A, &Q, &L, &b, &isu, &isp));
51: else PetscCall(LSCLoadOperators(&A, &Q, NULL, &b, &isu, &isp));
52: PetscCall(VecDuplicate(b, &x));
54: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp_A));
55: PetscCall(KSPSetOptionsPrefix(ksp_A, "fc_"));
56: PetscCall(KSPSetOperators(ksp_A, A, A));
58: PetscCall(KSPSetFromOptions(ksp_A));
59: PetscCall(KSPGetPC(ksp_A, &pc_A));
60: PetscCall(PCFieldSplitSetBlockSize(pc_A, 2));
61: PetscCall(PCFieldSplitSetIS(pc_A, "velocity", isu));
62: PetscCall(PCFieldSplitSetIS(pc_A, "pressure", isp));
64: // Need to call this before getting sub ksps, etc.
65: PetscCall(PCSetUp(pc_A));
66: PetscCall(PCFieldSplitGetSubKSP(pc_A, &num_splits, &subksp));
67: schur_complement_ksp = subksp[1];
69: PetscCall(KSPGetPC(schur_complement_ksp, &lsc_pc));
70: PetscCall(PCGetOperators(lsc_pc, NULL, &lsc_pc_pmat));
71: PetscCall(PetscObjectCompose((PetscObject)lsc_pc_pmat, "LSC_Qscale", (PetscObject)Q));
72: if (commute_lsc) PetscCall(PetscObjectCompose((PetscObject)lsc_pc_pmat, "LSC_L", (PetscObject)L));
74: PetscCall(KSPSolve(ksp_A, b, x));
76: PetscCall(KSPDestroy(&ksp_A));
77: PetscCall(MatDestroy(&A));
78: PetscCall(MatDestroy(&Q));
79: if (L) PetscCall(MatDestroy(&L));
80: PetscCall(VecDestroy(&x));
81: PetscCall(VecDestroy(&b));
82: PetscCall(ISDestroy(&isu));
83: PetscCall(ISDestroy(&isp));
84: PetscCall(PetscFree(subksp));
85: PetscFunctionReturn(PETSC_SUCCESS);
86: }
88: int main(int argc, char **argv)
89: {
90: PetscFunctionBeginUser;
91: PetscCall(PetscInitialize(&argc, &argv, 0, help));
92: PetscCall(port_lsd_bfbt());
93: PetscCall(PetscFinalize());
94: return 0;
95: }
97: /*TEST
99: test:
100: suffix: elman
101: args: -f ${DATAFILESPATH}/matrices/elman -fc_ksp_monitor -fc_ksp_type fgmres -fc_ksp_max_it 100 -fc_pc_type fieldsplit -fc_pc_fieldsplit_type SCHUR -fc_pc_fieldsplit_schur_fact_type full -fc_fieldsplit_velocity_ksp_type gmres -fc_fieldsplit_velocity_pc_type lu -fc_fieldsplit_pressure_ksp_max_it 100 -fc_fieldsplit_pressure_ksp_monitor -fc_fieldsplit_pressure_pc_type lsc -fc_fieldsplit_pressure_lsc_ksp_type gmres -fc_fieldsplit_pressure_lsc_ksp_max_it 100 -fc_fieldsplit_pressure_lsc_pc_type hypre -fc_fieldsplit_pressure_lsc_pc_hypre_type boomeramg -commute_lsc 0 -fc_ksp_converged_reason -fc_fieldsplit_pressure_ksp_converged_reason -fc_fieldsplit_pressure_ksp_type fgmres -fc_fieldsplit_pressure_pc_lsc_scale_diag
102: requires: datafilespath double hypre !complex defined(PETSC_USE_64BIT_INDICES)
104: test:
105: suffix: olshanskii
106: args: -f ${DATAFILESPATH}/matrices/olshanskii -fc_ksp_monitor -fc_ksp_type fgmres -fc_ksp_max_it 100 -fc_pc_type fieldsplit -fc_pc_fieldsplit_type SCHUR -fc_pc_fieldsplit_schur_fact_type full -fc_fieldsplit_velocity_ksp_type gmres -fc_fieldsplit_velocity_pc_type lu -fc_fieldsplit_pressure_ksp_max_it 100 -fc_fieldsplit_pressure_ksp_monitor -fc_fieldsplit_pressure_pc_type lsc -fc_fieldsplit_pressure_lsc_ksp_type gmres -fc_fieldsplit_pressure_lsc_ksp_max_it 100 -fc_fieldsplit_pressure_lsc_pc_type hypre -fc_fieldsplit_pressure_lsc_pc_hypre_type boomeramg -commute_lsc 1 -fc_ksp_converged_reason -fc_fieldsplit_pressure_ksp_converged_reason -fc_fieldsplit_pressure_ksp_type fgmres -fc_fieldsplit_pressure_lsc_mass_pc_type lu -fc_fieldsplit_pressure_lsc_mass_ksp_type gmres -fc_fieldsplit_pressure_lsc_mass_ksp_pc_side right -fc_fieldsplit_pressure_pc_lsc_commute 1
107: requires: datafilespath double hypre !complex defined(PETSC_USE_64BIT_INDICES)
109: TEST*/