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, NULL, 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_CURRENT, PETSC_CURRENT, 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*/