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:   PetscLogStage stage;

 38:   PetscFunctionBeginUser;
 39:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 40:   n1 = 64;
 41:   n2 = 64;
 42:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n1", &n1, NULL));
 43:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n2", &n2, NULL));

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

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

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

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

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

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

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

151:   /*
152:      Set right-hand side.
153:   */
154:   PetscCall(VecSet(b, 1.0));

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

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

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

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

187:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188:                       Solve the linear system
189:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

216: /*TEST

218:    test:
219:       nsize: 8
220:       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

222:    test:
223:       suffix: 2
224:       nsize: 8
225:       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
226:       output_file: output/ex38_1.out

228: TEST*/