Actual source code: ex86.c

      1: static char help[] = "Solves a one-dimensional steady upwind advection system with KSP.\n\n";

      3: /*
      4:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
      5:   automatically includes:
      6:      petscsys.h    - base PETSc routines   petscvec.h - vectors
      7:      petscmat.h    - matrices              petscpc.h  - preconditioners
      8:      petscis.h     - index sets
      9:      petscviewer.h - viewers
     10: */
     11: #include <petscksp.h>

     13: int main(int argc, char **args)
     14: {
     15:   Vec         x, b, work_vec; /* approx solution, RHS, work vector */
     16:   Mat         A;              /* linear system matrix */
     17:   KSP         ksp;            /* linear solver context */
     18:   PC          pc;             /* preconditioner context */
     19:   PetscInt    i, j, n = 10, col[2];
     20:   PetscScalar work_scalar, value[2];
     21:   PetscRandom r;

     23:   PetscFunctionBeginUser;
     24:   PetscCall(PetscInitialize(&argc, &args, NULL, help));

     26:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));

     28:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     29:          Compute the matrix and right-hand-side vector that define
     30:          the linear system, Ax = b.
     31:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     33:   /*
     34:      Create vectors.  Note that we form 1 vector from scratch and
     35:      then duplicate as needed.
     36:   */
     37:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
     38:   PetscCall(PetscObjectSetName((PetscObject)x, "Solution"));
     39:   PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
     40:   PetscCall(VecSetFromOptions(x));
     41:   PetscCall(VecDuplicate(x, &b));
     42:   PetscCall(VecDuplicate(x, &work_vec));

     44:   /*
     45:      Create matrix.  When using MatCreate(), the matrix format can
     46:      be specified at runtime.

     48:      Performance tuning note:  For problems of substantial size,
     49:      preallocation of matrix memory is crucial for attaining good
     50:      performance. See the matrix chapter of the users manual for details.
     51:   */
     52:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
     53:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
     54:   PetscCall(MatSetFromOptions(A));
     55:   PetscCall(MatSetUp(A));

     57:   /*
     58:      Assemble matrix
     59:   */
     60:   value[0] = -1.0;
     61:   value[1] = 1.0;
     62:   for (i = 1; i < n; i++) {
     63:     col[0] = i - 1;
     64:     col[1] = i;
     65:     PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
     66:   }
     67:   i           = 0;
     68:   j           = 0;
     69:   work_scalar = 1;
     70:   PetscCall(MatSetValues(A, 1, &i, 1, &j, &work_scalar, INSERT_VALUES));
     71:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
     72:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

     74:   /*
     75:      Create x, b
     76:   */
     77:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &r));
     78:   PetscCall(VecSetRandom(x, r));
     79:   PetscCall(VecSetRandom(work_vec, r));
     80:   PetscCall(MatMult(A, work_vec, b));

     82:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     83:                 Create the linear solver and set various options
     84:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     85:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));

     87:   /*
     88:      Set operators. Here the matrix that defines the linear system
     89:      also serves as the matrix that defines the preconditioner.
     90:   */
     91:   PetscCall(KSPSetOperators(ksp, A, A));

     93:   /*
     94:      Set linear solver defaults for this problem (optional).
     95:      - By extracting the KSP and PC contexts from the KSP context,
     96:        we can then directly call any KSP and PC routines to set
     97:        various options.
     98:      - The following four statements are optional; all of these
     99:        parameters could alternatively be specified at runtime via
    100:        KSPSetFromOptions();
    101:   */
    102:   PetscCall(KSPGetPC(ksp, &pc));
    103:   PetscCall(PCSetType(pc, PCBJACOBI));
    104:   PetscCall(KSPSetTolerances(ksp, 1.e-8, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT));

    106:   /*
    107:     Set runtime options, e.g.,
    108:         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
    109:     These options will override those specified above as long as
    110:     KSPSetFromOptions() is called _after_ any other customization
    111:     routines.
    112:   */
    113:   PetscCall(KSPSetFromOptions(ksp));

    115:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    116:                       Solve the linear system
    117:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    118:   PetscCall(KSPSolve(ksp, b, x));

    120:   /*
    121:      Free work space.  All PETSc objects should be destroyed when they
    122:      are no longer needed.
    123:   */
    124:   PetscCall(VecDestroy(&x));
    125:   PetscCall(VecDestroy(&b));
    126:   PetscCall(VecDestroy(&work_vec));
    127:   PetscCall(MatDestroy(&A));
    128:   PetscCall(KSPDestroy(&ksp));
    129:   PetscCall(PetscRandomDestroy(&r));

    131:   /*
    132:      Always call PetscFinalize() before exiting a program.  This routine
    133:        - finalizes the PETSc libraries as well as MPI
    134:        - provides summary and diagnostic information if certain runtime
    135:          options are chosen (e.g., -log_view).
    136:   */
    137:   PetscCall(PetscFinalize());
    138:   return 0;
    139: }

    141: /*TEST

    143:    testset:
    144:       requires: hypre !single !complex !defined(PETSC_HAVE_HYPRE_DEVICE)
    145:       args: -ksp_monitor -pc_type hypre -pc_hypre_type boomeramg -pc_hypre_boomeramg_grid_sweeps_down 0 -pc_hypre_boomeramg_grid_sweeps_up 1 -pc_hypre_boomeramg_grid_sweeps_coarse 2 -pc_hypre_boomeramg_max_levels 2 -ksp_rtol 1e-7 -pc_hypre_boomeramg_max_coarse_size 16 -n 33 -ksp_max_it 30 -pc_hypre_boomeramg_relax_type_all Jacobi
    146:       test:
    147:          suffix: hypre
    148:       test:
    149:          suffix: hypre_air
    150:          args: -pc_hypre_boomeramg_restriction_type 1 -pc_hypre_boomeramg_postrelax F -pc_hypre_boomeramg_grid_sweeps_down 1 -pc_hypre_boomeramg_prerelax C

    152: TEST*/