1 static const char help[] = "Solves a Q2-Q1 Navier-Stokes problem.\n\n";
2
3 #include <petscksp.h>
4 #include <petscdmda.h>
5
LSCLoadOperators(Mat * A,Mat * Q,Mat * L,Vec * rhs,IS * velocity,IS * pressure)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;
11
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 }
32
port_lsd_bfbt(void)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;
47
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));
53
54 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp_A));
55 PetscCall(KSPSetOptionsPrefix(ksp_A, "fc_"));
56 PetscCall(KSPSetOperators(ksp_A, A, A));
57
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));
63
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];
68
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));
73
74 PetscCall(KSPSolve(ksp_A, b, x));
75
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 }
87
main(int argc,char ** argv)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 }
96
97 /*TEST
98
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)
103
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)
108
109 TEST*/
110