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