Actual source code: ex38.c

  1: /*

  3: mpiexec -n 8 ./ex38 -ksp_type fbcgs -ksp_rtol 1.e-6 -sub_ksp_type bcgs -sub_ksp_rtol 1.e-3 -pc_type bjacobi -ksp_converged_reason -ksp_monitor -n1 64 -n2 64

  5:   Contributed by Jie Chen for testing flexible BiCGStab algorithm
  6: */

  8: static char help[] = "Solves the PDE (in 2D) -laplacian(u) + gamma x dot grad(u) + beta u = 1\n\
  9: with zero Dirichlet condition. The discretization is standard centered\n\
 10: difference. Input parameters include:\n\
 11:   -n1        : number of mesh points in 1st dimension (default 64)\n\
 12:   -n2        : number of mesh points in 2nd dimension (default 64)\n\
 13:   -h         : spacing between mesh points (default 1/n1)\n\
 14:   -gamma     : gamma (default 4/h)\n\
 15:   -beta      : beta (default 0.01/h^2)\n\n";

 17: /*
 18:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
 19:   automatically includes:
 20:      petscsys.h    - base PETSc routines   petscvec.h - vectors
 21:      petscmat.h    - matrices
 22:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 23:      petscviewer.h - viewers               petscpc.h  - preconditioners
 24: */
 25: #include <petscksp.h>

 27: int main(int argc, char **args)
 28: {
 29:   Vec         x, b, u;        /* approx solution, RHS, working vector */
 30:   Mat         A;              /* linear system matrix */
 31:   KSP         ksp;            /* linear solver context */
 32:   PetscInt    n1, n2;         /* parameters */
 33:   PetscReal   h, gamma, beta; /* parameters */
 34:   PetscInt    i, j, Ii, J, Istart, Iend;
 35:   PetscScalar v, co1, co2;
 36: #if defined(PETSC_USE_LOG)
 37:   PetscLogStage stage;
 38: #endif

 41:   PetscInitialize(&argc, &args, (char *)0, help);
 42:   n1 = 64;
 43:   n2 = 64;
 44:   PetscOptionsGetInt(NULL, NULL, "-n1", &n1, NULL);
 45:   PetscOptionsGetInt(NULL, NULL, "-n2", &n2, NULL);

 47:   h     = 1.0 / n1;
 48:   gamma = 4.0;
 49:   beta  = 0.01;
 50:   PetscOptionsGetReal(NULL, NULL, "-h", &h, NULL);
 51:   PetscOptionsGetReal(NULL, NULL, "-gamma", &gamma, NULL);
 52:   PetscOptionsGetReal(NULL, NULL, "-beta", &beta, NULL);
 53:   gamma = gamma / h;
 54:   beta  = beta / (h * h);

 56:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 57:          Compute the matrix and set right-hand-side vector.
 58:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 59:   /*
 60:      Create parallel matrix, specifying only its global dimensions.
 61:      When using MatCreate(), the matrix format can be specified at
 62:      runtime. Also, the parallel partitioning of the matrix is
 63:      determined by PETSc at runtime.

 65:      Performance tuning note:  For problems of substantial size,
 66:      preallocation of matrix memory is crucial for attaining good
 67:      performance. See the matrix chapter of the users manual for details.
 68:   */
 69:   MatCreate(PETSC_COMM_WORLD, &A);
 70:   MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n1 * n2, n1 * n2);
 71:   MatSetFromOptions(A);
 72:   MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL);
 73:   MatSeqAIJSetPreallocation(A, 5, NULL);
 74:   MatSetUp(A);

 76:   /*
 77:      Currently, all PETSc parallel matrix formats are partitioned by
 78:      contiguous chunks of rows across the processors.  Determine which
 79:      rows of the matrix are locally owned.
 80:   */
 81:   MatGetOwnershipRange(A, &Istart, &Iend);

 83:   /*
 84:      Set matrix elements for the 2-D, five-point stencil in parallel.
 85:       - Each processor needs to insert only elements that it owns
 86:         locally (but any non-local elements will be sent to the
 87:         appropriate processor during matrix assembly).
 88:       - Always specify global rows and columns of matrix entries.
 89:    */
 90:   PetscLogStageRegister("Assembly", &stage);
 91:   PetscLogStagePush(stage);
 92:   co1 = gamma * h * h / 2.0;
 93:   co2 = beta * h * h;
 94:   for (Ii = Istart; Ii < Iend; Ii++) {
 95:     i = Ii / n2;
 96:     j = Ii - i * n2;
 97:     if (i > 0) {
 98:       J = Ii - n2;
 99:       v = -1.0 + co1 * (PetscScalar)i;
100:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
101:     }
102:     if (i < n1 - 1) {
103:       J = Ii + n2;
104:       v = -1.0 + co1 * (PetscScalar)i;
105:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
106:     }
107:     if (j > 0) {
108:       J = Ii - 1;
109:       v = -1.0 + co1 * (PetscScalar)j;
110:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
111:     }
112:     if (j < n2 - 1) {
113:       J = Ii + 1;
114:       v = -1.0 + co1 * (PetscScalar)j;
115:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
116:     }
117:     v = 4.0 + co2;
118:     MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES);
119:   }

121:   /*
122:      Assemble matrix, using the 2-step process:
123:        MatAssemblyBegin(), MatAssemblyEnd()
124:      Computations can be done while messages are in transition
125:      by placing code between these two statements.
126:   */
127:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
128:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
129:   PetscLogStagePop();

131:   /*
132:      Create parallel vectors.
133:       - We form 1 vector from scratch and then duplicate as needed.
134:       - When using VecCreate(), VecSetSizes and VecSetFromOptions()
135:         in this example, we specify only the
136:         vector's global dimension; the parallel partitioning is determined
137:         at runtime.
138:       - When solving a linear system, the vectors and matrices MUST
139:         be partitioned accordingly.  PETSc automatically generates
140:         appropriately partitioned matrices and vectors when MatCreate()
141:         and VecCreate() are used with the same communicator.
142:       - The user can alternatively specify the local vector and matrix
143:         dimensions when more sophisticated partitioning is needed
144:         (replacing the PETSC_DECIDE argument in the VecSetSizes() statement
145:         below).
146:   */
147:   VecCreate(PETSC_COMM_WORLD, &b);
148:   VecSetSizes(b, PETSC_DECIDE, n1 * n2);
149:   VecSetFromOptions(b);
150:   VecDuplicate(b, &x);
151:   VecDuplicate(b, &u);

153:   /*
154:      Set right-hand side.
155:   */
156:   VecSet(b, 1.0);

158:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159:                 Create the linear solver and set various options
160:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161:   /*
162:      Create linear solver context
163:   */
164:   KSPCreate(PETSC_COMM_WORLD, &ksp);

166:   /*
167:      Set operators. Here the matrix that defines the linear system
168:      also serves as the preconditioning matrix.
169:   */
170:   KSPSetOperators(ksp, A, A);

172:   /*
173:      Set linear solver defaults for this problem (optional).
174:      - By extracting the KSP and PC contexts from the KSP context,
175:        we can then directly call any KSP and PC routines to set
176:        various options.
177:   */
178:   KSPSetTolerances(ksp, 1.e-6, PETSC_DEFAULT, PETSC_DEFAULT, 200);

180:   /*
181:     Set runtime options, e.g.,
182:         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
183:     These options will override those specified above as long as
184:     KSPSetFromOptions() is called _after_ any other customization
185:     routines.
186:   */
187:   KSPSetFromOptions(ksp);

189:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
190:                       Solve the linear system
191:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

193:   KSPSolve(ksp, b, x);

195:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196:                       Clean up
197:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
198:   /*
199:      Free work space.  All PETSc objects should be destroyed when they
200:      are no longer needed.
201:   */
202:   KSPDestroy(&ksp);
203:   VecDestroy(&u);
204:   VecDestroy(&x);
205:   VecDestroy(&b);
206:   MatDestroy(&A);

208:   /*
209:      Always call PetscFinalize() before exiting a program.  This routine
210:        - finalizes the PETSc libraries as well as MPI
211:        - provides summary and diagnostic information if certain runtime
212:          options are chosen (e.g., -log_view).
213:   */
214:   PetscFinalize();
215:   return 0;
216: }

218: /*TEST

220:    test:
221:       nsize: 8
222:       args: -ksp_type fbcgs -ksp_rtol 1.e-6 -sub_ksp_type bcgs -sub_ksp_rtol 1.e-3 -pc_type bjacobi -ksp_converged_reason -n1 64 -n2 64

224:    test:
225:       suffix: 2
226:       nsize: 8
227:       args: -ksp_type qmrcgs -ksp_rtol 1.e-6 -sub_ksp_type bcgs -sub_ksp_rtol 1.e-3 -pc_type bjacobi -ksp_converged_reason -n1 64 -n2 64
228:       output_file: output/ex38_1.out

230: TEST*/