Actual source code: ex39.c

  1: /*
  2: mpiexec -n 8 ./ex39 -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 32 -n2 32 -n3 32

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

  7: static char help[] = "Solves the PDE (in 3D) - laplacian(u) + gamma x dot grad(u) + beta u = 1\n\
  8: with zero Dirichlet condition. The discretization is standard centered\n\
  9: difference. Input parameters include:\n\
 10:   -n1        : number of mesh points in 1st dimension (default 32)\n\
 11:   -n2        : number of mesh points in 2nd dimension (default 32)\n\
 12:   -n3        : number of mesh points in 3rd dimension (default 32)\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: #include <petscksp.h>
 18: int main(int argc, char **args)
 19: {
 20:   Vec         x, b, u;        /* approx solution, RHS, working vector */
 21:   Mat         A;              /* linear system matrix */
 22:   KSP         ksp;            /* linear solver context */
 23:   PetscInt    n1, n2, n3;     /* parameters */
 24:   PetscReal   h, gamma, beta; /* parameters */
 25:   PetscInt    i, j, k, Ii, J, Istart, Iend;
 26:   PetscScalar v, co1, co2;

 29:   PetscInitialize(&argc, &args, (char *)0, help);
 30:   n1 = 32;
 31:   n2 = 32;
 32:   n3 = 32;
 33:   PetscOptionsGetInt(NULL, NULL, "-n1", &n1, NULL);
 34:   PetscOptionsGetInt(NULL, NULL, "-n2", &n2, NULL);
 35:   PetscOptionsGetInt(NULL, NULL, "-n3", &n3, NULL);

 37:   h     = 1.0 / n1;
 38:   gamma = 4.0 / h;
 39:   beta  = 0.01 / (h * h);
 40:   PetscOptionsGetReal(NULL, NULL, "-h", &h, NULL);
 41:   PetscOptionsGetReal(NULL, NULL, "-gamma", &gamma, NULL);
 42:   PetscOptionsGetReal(NULL, NULL, "-beta", &beta, NULL);

 44:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 45:          Compute the matrix and set right-hand-side vector.
 46:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 47:   MatCreate(PETSC_COMM_WORLD, &A);
 48:   MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n1 * n2 * n3, n1 * n2 * n3);
 49:   MatSetFromOptions(A);
 50:   MatMPIAIJSetPreallocation(A, 7, NULL, 7, NULL);
 51:   MatSeqAIJSetPreallocation(A, 7, NULL);
 52:   MatSetUp(A);
 53:   MatGetOwnershipRange(A, &Istart, &Iend);

 55:   /*
 56:      Set matrix elements for the 3-D, seven-point stencil in parallel.
 57:       - Each processor needs to insert only elements that it owns
 58:         locally (but any non-local elements will be sent to the
 59:         appropriate processor during matrix assembly).
 60:       - Always specify global rows and columns of matrix entries.
 61:    */
 62:   co1 = gamma * h * h / 2.0;
 63:   co2 = beta * h * h;
 64:   for (Ii = Istart; Ii < Iend; Ii++) {
 65:     i = Ii / (n2 * n3);
 66:     j = (Ii - i * n2 * n3) / n3;
 67:     k = Ii - i * n2 * n3 - j * n3;
 68:     if (i > 0) {
 69:       J = Ii - n2 * n3;
 70:       v = -1.0 + co1 * (PetscScalar)i;
 71:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 72:     }
 73:     if (i < n1 - 1) {
 74:       J = Ii + n2 * n3;
 75:       v = -1.0 + co1 * (PetscScalar)i;
 76:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 77:     }
 78:     if (j > 0) {
 79:       J = Ii - n3;
 80:       v = -1.0 + co1 * (PetscScalar)j;
 81:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 82:     }
 83:     if (j < n2 - 1) {
 84:       J = Ii + n3;
 85:       v = -1.0 + co1 * (PetscScalar)j;
 86:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 87:     }
 88:     if (k > 0) {
 89:       J = Ii - 1;
 90:       v = -1.0 + co1 * (PetscScalar)k;
 91:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 92:     }
 93:     if (k < n3 - 1) {
 94:       J = Ii + 1;
 95:       v = -1.0 + co1 * (PetscScalar)k;
 96:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 97:     }
 98:     v = 6.0 + co2;
 99:     MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES);
100:   }
101:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
102:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

104:   /* Create parallel vectors and Set right-hand side. */
105:   VecCreate(PETSC_COMM_WORLD, &b);
106:   VecSetSizes(b, PETSC_DECIDE, n1 * n2 * n3);
107:   VecSetFromOptions(b);
108:   VecDuplicate(b, &x);
109:   VecDuplicate(b, &u);
110:   VecSet(b, 1.0);

112:   /* Create linear solver context */
113:   KSPCreate(PETSC_COMM_WORLD, &ksp);
114:   KSPSetOperators(ksp, A, A);
115:   KSPSetTolerances(ksp, 1.e-6, PETSC_DEFAULT, PETSC_DEFAULT, 200);
116:   KSPSetFromOptions(ksp);

118:   /* Solve the linear system */
119:   KSPSolve(ksp, b, x);

121:   /* Free work space.  */
122:   KSPDestroy(&ksp);
123:   VecDestroy(&u);
124:   VecDestroy(&x);
125:   VecDestroy(&b);
126:   MatDestroy(&A);
127:   PetscFinalize();
128:   return 0;
129: }

131: /*TEST

133:    test:
134:       nsize: 8
135:       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 32 -n2 32 -n3 32

137:    test:
138:       suffix: 2
139:       nsize: 8
140:       args: -ksp_type fbcgsr -ksp_rtol 1.e-6 -sub_ksp_type bcgs -sub_ksp_rtol 1.e-3 -pc_type bjacobi -ksp_converged_reason -n1 32 -n2 32 -n3 32
141:       output_file: output/ex39_1.out

143:    test:
144:       suffix: 3
145:       nsize: 8
146:       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 32 -n2 32 -n3 32
147:       output_file: output/ex39_1.out

149: TEST*/