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, (char *)0, 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_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));

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