Actual source code: ex87.c

      1: #include <petscksp.h>
      2: #include <petsc/private/petscimpl.h>

      4: static char help[] = "Solves a saddle-point linear system using PCHPDDM.\n\n";

      6: static PetscErrorCode MatAndISLoad(const char *prefix, const char *identifier, Mat A, IS is, Mat N, PetscMPIInt rank, PetscMPIInt size);

      8: int main(int argc, char **args)
      9: {
     10:   Vec               b, x;            /* computed solution and RHS */
     11:   Mat               A[4], aux[2], S; /* linear system matrix */
     12:   KSP               ksp, *subksp;    /* linear solver context */
     13:   PC                pc;
     14:   IS                is[2];
     15:   PetscMPIInt       rank, size;
     16:   PetscInt          m, M, n, N, id = 0;
     17:   PetscViewer       viewer;
     18:   const char *const system[] = {"elasticity", "stokes"};
     19:   /* "elasticity":
     20:    *    2D linear elasticity with rubber-like and steel-like material coefficients, i.e., Poisson's ratio \in {0.4999, 0.35} and Young's modulus \in {0.01 GPa, 200.0 GPa}
     21:    *      discretized by order 2 (resp. 0) Lagrange finite elements in displacements (resp. pressure) on a triangle mesh
     22:    * "stokes":
     23:    *    2D lid-driven cavity with constant viscosity
     24:    *      discretized by order 2 (resp. 1) Lagrange finite elements, i.e., lowest-order Taylor--Hood finite elements, in velocities (resp. pressure) on a triangle mesh
     25:    *      if the option -empty_A11 is not set (or set to false), a pressure with a zero mean-value is computed
     26:    */
     27:   char      dir[PETSC_MAX_PATH_LEN], prefix[PETSC_MAX_PATH_LEN];
     28:   PetscBool flg[4] = {PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE};

     30:   PetscFunctionBeginUser;
     31:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
     32:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
     33:   PetscCheck(size == 4, PETSC_COMM_WORLD, PETSC_ERR_USER, "This example requires 4 processes");
     34:   PetscCall(PetscOptionsGetEList(NULL, NULL, "-system", system, PETSC_STATIC_ARRAY_LENGTH(system), &id, NULL));
     35:   if (id == 1) PetscCall(PetscOptionsGetBool(NULL, NULL, "-empty_A11", flg, NULL));
     36:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
     37:   for (PetscInt i = 0; i < 2; ++i) {
     38:     PetscCall(MatCreate(PETSC_COMM_WORLD, A + (i ? 3 : 0)));
     39:     PetscCall(ISCreate(PETSC_COMM_SELF, is + i));
     40:     PetscCall(MatCreate(PETSC_COMM_SELF, aux + i));
     41:   }
     42:   PetscCall(PetscStrncpy(dir, ".", sizeof(dir)));
     43:   PetscCall(PetscOptionsGetString(NULL, NULL, "-load_dir", dir, sizeof(dir), NULL));
     44:   /* loading matrices and auxiliary data for the diagonal blocks */
     45:   PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%s/%s", dir, id == 1 ? "B" : "A"));
     46:   PetscCall(MatAndISLoad(prefix, "00", A[0], is[0], aux[0], rank, size));
     47:   PetscCall(MatAndISLoad(prefix, "11", A[3], is[1], aux[1], rank, size));
     48:   /* loading the off-diagonal block with a coherent row/column layout */
     49:   PetscCall(MatCreate(PETSC_COMM_WORLD, A + 2));
     50:   PetscCall(MatGetLocalSize(A[0], &n, NULL));
     51:   PetscCall(MatGetSize(A[0], &N, NULL));
     52:   PetscCall(MatGetLocalSize(A[3], &m, NULL));
     53:   PetscCall(MatGetSize(A[3], &M, NULL));
     54:   PetscCall(MatSetSizes(A[2], m, n, M, N));
     55:   PetscCall(MatSetUp(A[2]));
     56:   PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%s/%s10.dat", dir, id == 1 ? "B" : "A"));
     57:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, prefix, FILE_MODE_READ, &viewer));
     58:   PetscCall(MatLoad(A[2], viewer));
     59:   PetscCall(PetscViewerDestroy(&viewer));
     60:   /* transposing the off-diagonal block */
     61:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-transpose", flg + 1, NULL));
     62:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-permute", flg + 2, NULL));
     63:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-explicit", flg + 3, NULL));
     64:   if (flg[1]) {
     65:     if (flg[2]) {
     66:       PetscCall(MatTranspose(A[2], MAT_INITIAL_MATRIX, A + 1));
     67:       PetscCall(MatDestroy(A + 2));
     68:     }
     69:     if (!flg[3]) PetscCall(MatCreateTranspose(A[2 - flg[2]], A + 1 + flg[2]));
     70:     else PetscCall(MatTranspose(A[2 - flg[2]], MAT_INITIAL_MATRIX, A + 1 + flg[2]));
     71:   } else {
     72:     if (flg[2]) {
     73:       PetscCall(MatHermitianTranspose(A[2], MAT_INITIAL_MATRIX, A + 1));
     74:       PetscCall(MatDestroy(A + 2));
     75:     }
     76:     if (!flg[3]) PetscCall(MatCreateHermitianTranspose(A[2 - flg[2]], A + 1 + flg[2]));
     77:     else PetscCall(MatHermitianTranspose(A[2 - flg[2]], MAT_INITIAL_MATRIX, A + 1 + flg[2]));
     78:   }
     79:   if (flg[0]) PetscCall(MatDestroy(A + 3));
     80:   else {
     81:     PetscCall(PetscOptionsGetBool(NULL, NULL, "-diagonal_A11", flg, NULL));
     82:     if (flg[0]) {
     83:       PetscCall(MatDestroy(A + 3));
     84:       PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD, m, m, M, M, PETSC_SMALL, A + 3));
     85:     }
     86:   }
     87:   flg[1] = PETSC_FALSE;
     88:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-all_transpose", flg + 1, NULL));
     89:   if (flg[1] && flg[2]) {
     90:     PetscCall(MatTranspose(A[1], MAT_INITIAL_MATRIX, &S));
     91:     PetscCall(MatDestroy(A + 1));
     92:     PetscCall(MatCreateHermitianTranspose(S, A + 1));
     93:     PetscCall(MatDestroy(&S));
     94:   }
     95:   /* global coefficient matrix */
     96:   PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, A, &S));
     97:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
     98:   PetscCall(KSPSetOperators(ksp, S, S));
     99:   PetscCall(KSPGetPC(ksp, &pc));
    100:   /* outer preconditioner */
    101:   PetscCall(PCSetType(pc, PCFIELDSPLIT));
    102:   PetscCall(PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR));
    103:   PetscCall(PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_SELF, NULL));
    104:   PetscCall(PCSetFromOptions(pc));
    105:   PetscCall(PCSetUp(pc));
    106:   PetscCall(PCFieldSplitGetSubKSP(pc, &n, &subksp));
    107:   PetscCall(KSPGetPC(subksp[0], &pc));
    108:   /* inner preconditioner associated to top-left block */
    109: #if defined(PETSC_HAVE_HPDDM) && defined(PETSC_HAVE_DYNAMIC_LIBRARIES) && defined(PETSC_USE_SHARED_LIBRARIES)
    110:   PetscCall(PCSetType(pc, PCHPDDM));
    111:   PetscCall(PCHPDDMSetAuxiliaryMat(pc, is[0], aux[0], NULL, NULL));
    112: #endif
    113:   PetscCall(PCSetFromOptions(pc));
    114:   PetscCall(KSPGetPC(subksp[1], &pc));
    115:   /* inner preconditioner associated to Schur complement, which will be set internally to a PCKSP */
    116: #if defined(PETSC_HAVE_HPDDM) && defined(PETSC_HAVE_DYNAMIC_LIBRARIES) && defined(PETSC_USE_SHARED_LIBRARIES)
    117:   PetscCall(PCSetType(pc, PCHPDDM));
    118:   if (!flg[0]) PetscCall(PCHPDDMSetAuxiliaryMat(pc, is[1], aux[1], NULL, NULL));
    119: #endif
    120:   PetscCall(PCSetFromOptions(pc));
    121:   PetscCall(PetscFree(subksp));
    122:   PetscCall(KSPSetFromOptions(ksp));
    123:   PetscCall(MatCreateVecs(S, &b, &x));
    124:   PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%s/rhs_%s.dat", dir, id == 1 ? "B" : "A"));
    125:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, prefix, FILE_MODE_READ, &viewer));
    126:   PetscCall(VecLoad(b, viewer));
    127:   PetscCall(PetscViewerDestroy(&viewer));
    128:   PetscCall(KSPSolve(ksp, b, x));
    129:   flg[1] = PETSC_FALSE;
    130:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-viewer", flg + 1, NULL));
    131:   if (flg[1]) PetscCall(PCView(pc, PETSC_VIEWER_STDOUT_WORLD));
    132:   flg[1] = PETSC_FALSE;
    133:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-successive_solves", flg + 1, NULL));
    134:   if (flg[1]) {
    135:     KSPConvergedReason reason[2];
    136:     PetscInt           iterations[2];
    137:     PetscCall(KSPGetConvergedReason(ksp, reason));
    138:     PetscCall(KSPGetTotalIterations(ksp, iterations));
    139:     PetscCall(KSPMonitorCancel(ksp));
    140:     PetscCall(PetscOptionsClearValue(NULL, "-ksp_monitor"));
    141:     PetscCall(PetscObjectStateIncrease((PetscObject)S));
    142:     PetscCall(KSPSetUp(ksp));
    143:     PetscCall(KSPGetPC(ksp, &pc));
    144:     PetscCall(PCFieldSplitGetSubKSP(pc, &n, &subksp));
    145:     PetscCall(KSPGetPC(subksp[0], &pc));
    146: #if defined(PETSC_HAVE_HPDDM) && defined(PETSC_HAVE_DYNAMIC_LIBRARIES) && defined(PETSC_USE_SHARED_LIBRARIES)
    147:     PetscCall(PCHPDDMSetAuxiliaryMat(pc, is[0], aux[0], NULL, NULL));
    148: #endif
    149:     PetscCall(PCSetFromOptions(pc));
    150:     PetscCall(KSPGetPC(subksp[1], &pc));
    151: #if defined(PETSC_HAVE_HPDDM) && defined(PETSC_HAVE_DYNAMIC_LIBRARIES) && defined(PETSC_USE_SHARED_LIBRARIES)
    152:     if (!flg[0]) PetscCall(PCHPDDMSetAuxiliaryMat(pc, is[1], aux[1], NULL, NULL));
    153: #endif
    154:     PetscCall(PCSetFromOptions(pc));
    155:     PetscCall(PetscFree(subksp));
    156:     PetscCall(KSPSolve(ksp, b, x));
    157:     PetscCall(KSPGetConvergedReason(ksp, reason + 1));
    158:     PetscCall(KSPGetTotalIterations(ksp, iterations + 1));
    159:     iterations[1] -= iterations[0];
    160:     PetscCheck(reason[0] == reason[1] && PetscAbs(iterations[0] - iterations[1]) <= 3, PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "Successive calls to KSPSolve() did not converge for the same reason (%s v. %s) or with the same number of iterations (+/- 3, %" PetscInt_FMT " v. %" PetscInt_FMT ")", KSPConvergedReasons[reason[0]], KSPConvergedReasons[reason[1]], iterations[0], iterations[1]);
    161:   }
    162:   PetscCall(VecDestroy(&x));
    163:   PetscCall(VecDestroy(&b));
    164:   PetscCall(KSPDestroy(&ksp));
    165:   PetscCall(MatDestroy(&S));
    166:   PetscCall(MatDestroy(A + 1));
    167:   PetscCall(MatDestroy(A + 2));
    168:   for (PetscInt i = 0; i < 2; ++i) {
    169:     PetscCall(MatDestroy(A + (i ? 3 : 0)));
    170:     PetscCall(MatDestroy(aux + i));
    171:     PetscCall(ISDestroy(is + i));
    172:   }
    173:   PetscCall(PetscFinalize());
    174:   return 0;
    175: }

    177: PetscErrorCode MatAndISLoad(const char *prefix, const char *identifier, Mat A, IS is, Mat aux, PetscMPIInt rank, PetscMPIInt size)
    178: {
    179:   IS              sizes;
    180:   const PetscInt *idx;
    181:   PetscViewer     viewer;
    182:   char            name[PETSC_MAX_PATH_LEN];

    184:   PetscFunctionBeginUser;
    185:   PetscCall(PetscSNPrintf(name, sizeof(name), "%s%s_sizes_%d_%d.dat", prefix, identifier, rank, size));
    186:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, name, FILE_MODE_READ, &viewer));
    187:   PetscCall(ISCreate(PETSC_COMM_SELF, &sizes));
    188:   PetscCall(ISLoad(sizes, viewer));
    189:   PetscCall(ISGetIndices(sizes, &idx));
    190:   PetscCall(MatSetSizes(A, idx[0], idx[1], idx[2], idx[3]));
    191:   PetscCall(MatSetUp(A));
    192:   PetscCall(ISRestoreIndices(sizes, &idx));
    193:   PetscCall(ISDestroy(&sizes));
    194:   PetscCall(PetscViewerDestroy(&viewer));
    195:   PetscCall(PetscSNPrintf(name, sizeof(name), "%s%s.dat", prefix, identifier));
    196:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, name, FILE_MODE_READ, &viewer));
    197:   PetscCall(MatLoad(A, viewer));
    198:   PetscCall(PetscViewerDestroy(&viewer));
    199:   PetscCall(PetscSNPrintf(name, sizeof(name), "%s%s_is_%d_%d.dat", prefix, identifier, rank, size));
    200:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, name, FILE_MODE_READ, &viewer));
    201:   PetscCall(ISLoad(is, viewer));
    202:   PetscCall(PetscViewerDestroy(&viewer));
    203:   PetscCall(PetscSNPrintf(name, sizeof(name), "%s%s_aux_%d_%d.dat", prefix, identifier, rank, size));
    204:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, name, FILE_MODE_READ, &viewer));
    205:   PetscCall(MatLoad(aux, viewer));
    206:   PetscCall(PetscViewerDestroy(&viewer));
    207:   PetscFunctionReturn(PETSC_SUCCESS);
    208: }

    210: /*TEST

    212:    testset:
    213:       requires: datafilespath hpddm slepc double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
    214:       nsize: 4
    215:       args: -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO -ksp_monitor -ksp_rtol 1e-4 -fieldsplit_ksp_max_it 100 -fieldsplit_pc_hpddm_levels_1_eps_nev 10 -fieldsplit_pc_hpddm_levels_1_st_share_sub_ksp -fieldsplit_pc_hpddm_has_neumann -fieldsplit_pc_hpddm_define_subdomains -fieldsplit_1_pc_hpddm_schur_precondition geneo -fieldsplit_pc_hpddm_coarse_pc_type redundant -fieldsplit_pc_hpddm_coarse_redundant_pc_type cholesky -fieldsplit_pc_hpddm_levels_1_sub_pc_type lu -fieldsplit_ksp_type fgmres -ksp_type fgmres -ksp_max_it 10 -fieldsplit_1_pc_hpddm_coarse_correction balanced -fieldsplit_1_pc_hpddm_levels_1_eps_gen_non_hermitian -fieldsplit_1_pc_hpddm_coarse_p 2
    216:       test:
    217:         requires: mumps
    218:         suffix: 1
    219:         args: -viewer -system {{elasticity stokes}separate output} -fieldsplit_1_pc_hpddm_ksp_pc_side left -fieldsplit_1_pc_hpddm_levels_1_sub_mat_mumps_icntl_26 1
    220:         filter: grep -v -e "action of " -e "                            " -e "block size" -e "total: nonzeros=" -e "using I-node" -e "aij" -e "transpose" -e "diagonal" -e "total number of" -e "                rows="
    221:       test:
    222:         requires: mumps
    223:         suffix: 2
    224:         output_file: output/ex87_1_system-stokes.out
    225:         args: -viewer -system stokes -empty_A11 -transpose {{false true}shared output} -permute {{false true}shared output} -fieldsplit_1_pc_hpddm_ksp_pc_side right -fieldsplit_1_pc_hpddm_coarse_mat_type baij -fieldsplit_1_pc_hpddm_levels_1_sub_mat_mumps_icntl_26 1 -explicit {{false true}shared output}
    226:         filter: grep -v -e "action of " -e "                            " -e "block size" -e "total: nonzeros=" -e "using I-node" -e "aij" -e "transpose" -e "diagonal" -e "total number of" -e "                rows=" | sed -e "s/      right preconditioning/      left preconditioning/g" -e "s/      using UNPRECONDITIONED/      using PRECONDITIONED/g"
    227:       test:
    228:         suffix: 1_petsc
    229:         args: -system {{elasticity stokes}separate output} -fieldsplit_1_pc_hpddm_ksp_pc_side left -fieldsplit_1_pc_hpddm_levels_1_sub_pc_factor_mat_solver_type petsc -fieldsplit_1_pc_hpddm_levels_1_eps_threshold 0.3 -permute
    230:       test:
    231:         suffix: 2_petsc
    232:         output_file: output/ex87_1_petsc_system-stokes.out
    233:         args: -system stokes -empty_A11 -transpose -fieldsplit_1_pc_hpddm_ksp_pc_side right -fieldsplit_1_pc_hpddm_levels_1_sub_pc_factor_mat_solver_type petsc -fieldsplit_1_pc_hpddm_coarse_mat_type baij -fieldsplit_1_pc_hpddm_levels_1_eps_threshold 0.3 -fieldsplit_1_pc_hpddm_levels_1_sub_pc_factor_shift_type inblocks -successive_solves
    234:         filter: sed -e "s/type: transpose/type: hermitiantranspose/g"
    235:       test:
    236:         suffix: threshold
    237:         requires: !defined(PETSC_HAVE_MKL_SPARSE_SP2M_FEATURE)
    238:         output_file: output/ex87_1_petsc_system-elasticity.out
    239:         args: -fieldsplit_1_pc_hpddm_ksp_pc_side left -fieldsplit_1_pc_hpddm_levels_1_eps_threshold 0.2 -fieldsplit_1_pc_hpddm_coarse_mat_type {{baij sbaij}shared output} -successive_solves
    240:    testset:
    241:       requires: datafilespath hpddm slepc double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
    242:       nsize: 4
    243:       args: -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO -ksp_monitor -ksp_rtol 1e-4 -fieldsplit_ksp_max_it 100 -fieldsplit_pc_hpddm_levels_1_st_share_sub_ksp -fieldsplit_pc_hpddm_define_subdomains -fieldsplit_1_pc_hpddm_schur_precondition geneo -fieldsplit_pc_hpddm_coarse_pc_type redundant -fieldsplit_pc_hpddm_coarse_redundant_pc_type cholesky -fieldsplit_pc_hpddm_levels_1_sub_pc_type lu -fieldsplit_ksp_type fgmres -ksp_type fgmres -ksp_max_it 10 -fieldsplit_1_pc_hpddm_coarse_correction balanced -fieldsplit_1_pc_hpddm_levels_1_eps_gen_non_hermitian -fieldsplit_1_pc_hpddm_coarse_p 2 -system stokes -fieldsplit_1_pc_hpddm_ksp_pc_side left -fieldsplit_1_pc_hpddm_levels_1_sub_pc_factor_mat_solver_type petsc -fieldsplit_1_pc_hpddm_levels_1_eps_threshold 0.3
    244:       test:
    245:         suffix: diagonal
    246:         output_file: output/ex87_1_petsc_system-stokes.out
    247:         args: -fieldsplit_pc_hpddm_levels_1_eps_nev 10 -fieldsplit_0_pc_hpddm_has_neumann -diagonal_A11 {{false true}shared output}
    248:       test:
    249:         suffix: harmonic_overlap_2
    250:         output_file: output/ex87_1_petsc_system-stokes.out
    251:         args: -fieldsplit_0_pc_hpddm_harmonic_overlap 2 -fieldsplit_0_pc_hpddm_levels_1_svd_nsv 20 -diagonal_A11 -permute {{false true}shared output} -all_transpose

    253:    test:
    254:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) !hpddm !memkind
    255:       nsize: 4
    256:       suffix: selfp
    257:       output_file: output/ex41_1.out
    258:       filter: grep -v "CONVERGED_RTOL iterations"
    259:       args: -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO -system stokes -ksp_rtol 1e-4 -ksp_converged_reason -ksp_max_it 30 -pc_type fieldsplit -pc_fieldsplit_type schur -fieldsplit_ksp_type preonly -pc_fieldsplit_schur_precondition selfp -fieldsplit_pc_type bjacobi -fieldsplit_sub_pc_type lu -transpose {{false true}shared output} -fieldsplit_1_mat_schur_complement_ainv_type lump

    261: TEST*/