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