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;

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

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

 44:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 45:          Compute the matrix and set right-hand-side vector.
 46:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 47:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 48:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n1 * n2 * n3, n1 * n2 * n3));
 49:   PetscCall(MatSetFromOptions(A));
 50:   PetscCall(MatMPIAIJSetPreallocation(A, 7, NULL, 7, NULL));
 51:   PetscCall(MatSeqAIJSetPreallocation(A, 7, NULL));
 52:   PetscCall(MatSetUp(A));
 53:   PetscCall(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:       PetscCall(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:       PetscCall(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:       PetscCall(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:       PetscCall(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:       PetscCall(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:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
 97:     }
 98:     v = 6.0 + co2;
 99:     PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
100:   }
101:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
102:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

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

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

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

121:   /* Free work space.  */
122:   PetscCall(KSPDestroy(&ksp));
123:   PetscCall(VecDestroy(&u));
124:   PetscCall(VecDestroy(&x));
125:   PetscCall(VecDestroy(&b));
126:   PetscCall(MatDestroy(&A));
127:   PetscCall(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*/