Actual source code: ex62.c

      1: static char help[] = "Illustrates use of PCGASM.\n\
      2: The Generalized Additive Schwarz Method for solving a linear system in parallel with KSP.  The\n\
      3: code indicates the procedure for setting user-defined subdomains.\n\
      4: See section 'ex62' below for command-line options.\n\
      5: Without -user_set_subdomains, the general PCGASM options are meaningful:\n\
      6:   -pc_gasm_total_subdomains\n\
      7:   -pc_gasm_print_subdomains\n\
      8: \n";

     10: /*
     11:    Note:  This example focuses on setting the subdomains for the GASM
     12:    preconditioner for a problem on a 2D rectangular grid.  See ex1.c
     13:    and ex2.c for more detailed comments on the basic usage of KSP
     14:    (including working with matrices and vectors).

     16:    The GASM preconditioner is fully parallel.  The user-space routine
     17:    CreateSubdomains2D that computes the domain decomposition is also parallel
     18:    and attempts to generate both subdomains straddling processors and multiple
     19:    domains per processor.

     21:    This matrix in this linear system arises from the discretized Laplacian,
     22:    and thus is not very interesting in terms of experimenting with variants
     23:    of the GASM preconditioner.
     24: */

     26: /*
     27:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
     28:   automatically includes:
     29:      petscsys.h    - base PETSc routines   petscvec.h - vectors
     30:      petscmat.h    - matrices
     31:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
     32:      petscviewer.h - viewers               petscpc.h  - preconditioners
     33: */
     34: #include <petscksp.h>

     36: PetscErrorCode AssembleMatrix(Mat, PetscInt m, PetscInt n);

     38: int main(int argc, char **args)
     39: {
     40:   Vec         x, b, u;           /* approx solution, RHS, exact solution */
     41:   Mat         A;                 /* linear system matrix */
     42:   KSP         ksp;               /* linear solver context */
     43:   PC          pc;                /* PC context */
     44:   IS         *inneris, *outeris; /* array of index sets that define the subdomains */
     45:   PetscInt    overlap;           /* width of subdomain overlap */
     46:   PetscInt    Nsub;              /* number of subdomains */
     47:   PetscInt    m, n;              /* mesh dimensions in x- and y- directions */
     48:   PetscInt    M, N;              /* number of subdomains in x- and y- directions */
     49:   PetscMPIInt size;
     50:   PetscBool   flg                 = PETSC_FALSE;
     51:   PetscBool   user_set_subdomains = PETSC_FALSE;
     52:   PetscReal   one, e;

     54:   PetscFunctionBeginUser;
     55:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
     56:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
     57:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex62", "PCGASM");
     58:   m = 15;
     59:   PetscCall(PetscOptionsInt("-M", "Number of mesh points in the x-direction", "PCGASMCreateSubdomains2D", m, &m, NULL));
     60:   n = 17;
     61:   PetscCall(PetscOptionsInt("-N", "Number of mesh points in the y-direction", "PCGASMCreateSubdomains2D", n, &n, NULL));
     62:   user_set_subdomains = PETSC_FALSE;
     63:   PetscCall(PetscOptionsBool("-user_set_subdomains", "Use the user-specified 2D tiling of mesh by subdomains", "PCGASMCreateSubdomains2D", user_set_subdomains, &user_set_subdomains, NULL));
     64:   M = 2;
     65:   PetscCall(PetscOptionsInt("-Mdomains", "Number of subdomain tiles in the x-direction", "PCGASMSetSubdomains2D", M, &M, NULL));
     66:   N = 1;
     67:   PetscCall(PetscOptionsInt("-Ndomains", "Number of subdomain tiles in the y-direction", "PCGASMSetSubdomains2D", N, &N, NULL));
     68:   overlap = 1;
     69:   PetscCall(PetscOptionsInt("-overlap", "Size of tile overlap.", "PCGASMSetSubdomains2D", overlap, &overlap, NULL));
     70:   PetscOptionsEnd();

     72:   /* -------------------------------------------------------------------
     73:          Compute the matrix and right-hand-side vector that define
     74:          the linear system, Ax = b.
     75:      ------------------------------------------------------------------- */

     77:   /*
     78:      Assemble the matrix for the five point stencil, YET AGAIN
     79:   */
     80:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
     81:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n));
     82:   PetscCall(MatSetFromOptions(A));
     83:   PetscCall(MatSetUp(A));
     84:   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE));
     85:   PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
     86:   PetscCall(AssembleMatrix(A, m, n));

     88:   /*
     89:      Create and set vectors
     90:   */
     91:   PetscCall(VecCreate(PETSC_COMM_WORLD, &b));
     92:   PetscCall(VecSetSizes(b, PETSC_DECIDE, m * n));
     93:   PetscCall(VecSetFromOptions(b));
     94:   PetscCall(VecDuplicate(b, &u));
     95:   PetscCall(VecDuplicate(b, &x));
     96:   one = 1.0;
     97:   PetscCall(VecSet(u, one));
     98:   PetscCall(MatMult(A, u, b));

    100:   /*
    101:      Create linear solver context
    102:   */
    103:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));

    105:   /*
    106:      Set operators. Here the matrix that defines the linear system
    107:      also serves as the matrix from which the preconditioner is constructed.
    108:   */
    109:   PetscCall(KSPSetOperators(ksp, A, A));

    111:   /*
    112:      Set the default preconditioner for this program to be GASM
    113:   */
    114:   PetscCall(KSPGetPC(ksp, &pc));
    115:   PetscCall(PCSetType(pc, PCGASM));

    117:   /* -------------------------------------------------------------------
    118:                   Define the problem decomposition
    119:      ------------------------------------------------------------------- */

    121:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    122:        Basic method, should be sufficient for the needs of many users.
    123:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    125:      Set the overlap, using the default PETSc decomposition via
    126:          PCGASMSetOverlap(pc,overlap);
    127:      Could instead use the option -pc_gasm_overlap <ovl>

    129:      Set the total number of blocks via -pc_gasm_blocks <blks>
    130:      Note:  The GASM default is to use 1 block per processor.  To
    131:      experiment on a single processor with various overlaps, you
    132:      must specify use of multiple blocks!
    133:   */

    135:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    136:        More advanced method, setting user-defined subdomains
    137:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    139:      Firstly, create index sets that define the subdomains.  The utility
    140:      routine PCGASMCreateSubdomains2D() is a simple example, which partitions
    141:      the 2D grid into MxN subdomains with an optional overlap.
    142:      More generally, the user should write a custom routine for a particular
    143:      problem geometry.

    145:      Then call PCGASMSetLocalSubdomains() with resulting index sets
    146:      to set the subdomains for the GASM preconditioner.
    147:   */

    149:   if (user_set_subdomains) { /* user-control version */
    150:     PetscCall(PCGASMCreateSubdomains2D(pc, m, n, M, N, 1, overlap, &Nsub, &inneris, &outeris));
    151:     PetscCall(PCGASMSetSubdomains(pc, Nsub, inneris, outeris));
    152:     PetscCall(PCGASMDestroySubdomains(Nsub, &inneris, &outeris));
    153:     flg = PETSC_FALSE;
    154:     PetscCall(PetscOptionsGetBool(NULL, NULL, "-subdomain_view", &flg, NULL));
    155:     if (flg) {
    156:       PetscInt i;
    157:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Nmesh points: %" PetscInt_FMT " x %" PetscInt_FMT "; subdomain partition: %" PetscInt_FMT " x %" PetscInt_FMT "; overlap: %" PetscInt_FMT "; Nsub: %" PetscInt_FMT "\n", m, n, M, N, overlap, Nsub));
    158:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Outer IS:\n"));
    159:       for (i = 0; i < Nsub; i++) {
    160:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  outer IS[%" PetscInt_FMT "]\n", i));
    161:         PetscCall(ISView(outeris[i], PETSC_VIEWER_STDOUT_SELF));
    162:       }
    163:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Inner IS:\n"));
    164:       for (i = 0; i < Nsub; i++) {
    165:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  inner IS[%" PetscInt_FMT "]\n", i));
    166:         PetscCall(ISView(inneris[i], PETSC_VIEWER_STDOUT_SELF));
    167:       }
    168:     }
    169:   } else { /* basic setup */
    170:     PetscCall(KSPSetFromOptions(ksp));
    171:   }

    173:   /* -------------------------------------------------------------------
    174:                 Set the linear solvers for the subblocks
    175:      ------------------------------------------------------------------- */

    177:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    178:        Basic method, should be sufficient for the needs of most users.
    179:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    181:      By default, the GASM preconditioner uses the same solver on each
    182:      block of the problem.  To set the same solver options on all blocks,
    183:      use the prefix -sub before the usual PC and KSP options, e.g.,
    184:           -sub_pc_type <pc> -sub_ksp_type <ksp> -sub_ksp_rtol 1.e-4

    186:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    187:         Advanced method, setting different solvers for various blocks.
    188:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    190:      Note that each block's KSP context is completely independent of
    191:      the others, and the full range of uniprocessor KSP options is
    192:      available for each block.

    194:      - Use PCGASMGetSubKSP() to extract the array of KSP contexts for
    195:        the local blocks.
    196:      - See ex7.c for a simple example of setting different linear solvers
    197:        for the individual blocks for the block Jacobi method (which is
    198:        equivalent to the GASM method with zero overlap).
    199:   */

    201:   flg = PETSC_FALSE;
    202:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-user_set_subdomain_solvers", &flg, NULL));
    203:   if (flg) {
    204:     KSP      *subksp;           /* array of KSP contexts for local subblocks */
    205:     PetscInt  i, nlocal, first; /* number of local subblocks, first local subblock */
    206:     PC        subpc;            /* PC context for subblock */
    207:     PetscBool isasm;

    209:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "User explicitly sets subdomain solvers.\n"));

    211:     /*
    212:        Set runtime options
    213:     */
    214:     PetscCall(KSPSetFromOptions(ksp));

    216:     /*
    217:        Flag an error if PCTYPE is changed from the runtime options
    218:      */
    219:     PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &isasm));
    220:     PetscCheck(isasm, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Cannot Change the PCTYPE when manually changing the subdomain solver settings");

    222:     /*
    223:        Call KSPSetUp() to set the block Jacobi data structures (including
    224:        creation of an internal KSP context for each block).

    226:        Note: KSPSetUp() MUST be called before PCGASMGetSubKSP().
    227:     */
    228:     PetscCall(KSPSetUp(ksp));

    230:     /*
    231:        Extract the array of KSP contexts for the local blocks
    232:     */
    233:     PetscCall(PCGASMGetSubKSP(pc, &nlocal, &first, &subksp));

    235:     /*
    236:        Loop over the local blocks, setting various KSP options
    237:        for each block.
    238:     */
    239:     for (i = 0; i < nlocal; i++) {
    240:       PetscCall(KSPGetPC(subksp[i], &subpc));
    241:       PetscCall(PCSetType(subpc, PCILU));
    242:       PetscCall(KSPSetType(subksp[i], KSPGMRES));
    243:       PetscCall(KSPSetTolerances(subksp[i], 1.e-7, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));
    244:     }
    245:   } else {
    246:     /*
    247:        Set runtime options
    248:     */
    249:     PetscCall(KSPSetFromOptions(ksp));
    250:   }

    252:   /* -------------------------------------------------------------------
    253:                       Solve the linear system
    254:      ------------------------------------------------------------------- */

    256:   PetscCall(KSPSolve(ksp, b, x));

    258:   /* -------------------------------------------------------------------
    259:         Assemble the matrix again to test repeated setup and solves.
    260:      ------------------------------------------------------------------- */

    262:   PetscCall(AssembleMatrix(A, m, n));
    263:   PetscCall(KSPSolve(ksp, b, x));

    265:   /* -------------------------------------------------------------------
    266:                       Compare result to the exact solution
    267:      ------------------------------------------------------------------- */
    268:   PetscCall(VecAXPY(x, -1.0, u));
    269:   PetscCall(VecNorm(x, NORM_INFINITY, &e));

    271:   flg = PETSC_FALSE;
    272:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-print_error", &flg, NULL));
    273:   if (flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Infinity norm of the error: %g\n", (double)e));

    275:   /*
    276:      Free work space.  All PETSc objects should be destroyed when they
    277:      are no longer needed.
    278:   */

    280:   PetscCall(KSPDestroy(&ksp));
    281:   PetscCall(VecDestroy(&u));
    282:   PetscCall(VecDestroy(&x));
    283:   PetscCall(VecDestroy(&b));
    284:   PetscCall(MatDestroy(&A));
    285:   PetscCall(PetscFinalize());
    286:   return 0;
    287: }

    289: PetscErrorCode AssembleMatrix(Mat A, PetscInt m, PetscInt n)
    290: {
    291:   PetscInt    i, j, Ii, J, Istart, Iend;
    292:   PetscScalar v;

    294:   PetscFunctionBegin;
    295:   PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));
    296:   for (Ii = Istart; Ii < Iend; Ii++) {
    297:     v = -1.0;
    298:     i = Ii / n;
    299:     j = Ii - i * n;
    300:     if (i > 0) {
    301:       J = Ii - n;
    302:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
    303:     }
    304:     if (i < m - 1) {
    305:       J = Ii + n;
    306:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
    307:     }
    308:     if (j > 0) {
    309:       J = Ii - 1;
    310:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
    311:     }
    312:     if (j < n - 1) {
    313:       J = Ii + 1;
    314:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
    315:     }
    316:     v = 4.0;
    317:     PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
    318:   }
    319:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    320:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    321:   PetscFunctionReturn(PETSC_SUCCESS);
    322: }

    324: /*TEST

    326:    test:
    327:       suffix: 2D_1
    328:       args: -M 7 -N 9 -user_set_subdomains -Mdomains 1 -Ndomains 3 -overlap 1 -print_error -pc_gasm_print_subdomains

    330:    test:
    331:       suffix: 2D_2
    332:       nsize: 2
    333:       args: -M 7 -N 9 -user_set_subdomains -Mdomains 1 -Ndomains 3 -overlap 1 -print_error -pc_gasm_print_subdomains

    335:    test:
    336:       suffix: 2D_3
    337:       nsize: 3
    338:       args: -M 7 -N 9 -user_set_subdomains -Mdomains 1 -Ndomains 3 -overlap 1 -print_error -pc_gasm_print_subdomains

    340:    test:
    341:       suffix: hp
    342:       nsize: 4
    343:       requires: superlu_dist
    344:       args: -M 7 -N 9 -pc_gasm_overlap 1 -sub_pc_type lu -sub_pc_factor_mat_solver_type superlu_dist -ksp_monitor -print_error -pc_gasm_total_subdomains 2 -pc_gasm_use_hierachical_partitioning 1
    345:       output_file: output/empty.out
    346:       TODO: bug, triggers New nonzero at (0,15) caused a malloc in MatCreateSubMatrices_MPIAIJ_SingleIS_Local

    348:    test:
    349:       suffix: superlu_dist_1
    350:       requires: superlu_dist
    351:       args: -M 7 -N 9 -print_error -pc_gasm_total_subdomains 1 -pc_gasm_print_subdomains -sub_pc_type lu -sub_pc_factor_mat_solver_type superlu_dist

    353:    test:
    354:       suffix: superlu_dist_2
    355:       nsize: 2
    356:       requires: superlu_dist
    357:       args: -M 7 -N 9 -print_error -pc_gasm_total_subdomains 1 -pc_gasm_print_subdomains -sub_pc_type lu -sub_pc_factor_mat_solver_type superlu_dist

    359:    test:
    360:       suffix: superlu_dist_3
    361:       nsize: 3
    362:       requires: superlu_dist
    363:       args: -M 7 -N 9 -print_error -pc_gasm_total_subdomains 2 -pc_gasm_print_subdomains -sub_pc_type lu -sub_pc_factor_mat_solver_type superlu_dist

    365:    test:
    366:       suffix: superlu_dist_4
    367:       nsize: 4
    368:       requires: superlu_dist
    369:       args: -M 7 -N 9 -print_error -pc_gasm_total_subdomains 2 -pc_gasm_print_subdomains -sub_pc_type lu -sub_pc_factor_mat_solver_type superlu_dist

    371:    test:
    372:       suffix: gasm_view
    373:       nsize: 8
    374:       args: -pc_type gasm -ksp_view -Mdomains 2 -Ndomains 2 -user_set_subdomains

    376: TEST*/