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: PetscInitialize(&argc,&args,(char*)0,help);
29: n1 = 32;
30: n2 = 32;
31: n3 = 32;
32: PetscOptionsGetInt(NULL,NULL,"-n1",&n1,NULL);
33: PetscOptionsGetInt(NULL,NULL,"-n2",&n2,NULL);
34: PetscOptionsGetInt(NULL,NULL,"-n3",&n3,NULL);
36: h = 1.0/n1;
37: gamma = 4.0/h;
38: beta = 0.01/(h*h);
39: PetscOptionsGetReal(NULL,NULL,"-h",&h,NULL);
40: PetscOptionsGetReal(NULL,NULL,"-gamma",&gamma,NULL);
41: PetscOptionsGetReal(NULL,NULL,"-beta",&beta,NULL);
43: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
44: Compute the matrix and set right-hand-side vector.
45: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
46: MatCreate(PETSC_COMM_WORLD,&A);
47: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n1*n2*n3,n1*n2*n3);
48: MatSetFromOptions(A);
49: MatMPIAIJSetPreallocation(A,7,NULL,7,NULL);
50: MatSeqAIJSetPreallocation(A,7,NULL);
51: MatSetUp(A);
52: MatGetOwnershipRange(A,&Istart,&Iend);
54: /*
55: Set matrix elements for the 3-D, seven-point stencil in parallel.
56: - Each processor needs to insert only elements that it owns
57: locally (but any non-local elements will be sent to the
58: appropriate processor during matrix assembly).
59: - Always specify global rows and columns of matrix entries.
60: */
61: co1 = gamma * h * h / 2.0;
62: co2 = beta * h * h;
63: for (Ii=Istart; Ii<Iend; Ii++) {
64: i = Ii/(n2*n3); j = (Ii - i*n2*n3)/n3; k = Ii - i*n2*n3 - j*n3;
65: if (i>0) {
66: J = Ii - n2*n3; v = -1.0 + co1*(PetscScalar)i;
67: MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);
68: }
69: if (i<n1-1) {
70: J = Ii + n2*n3; v = -1.0 + co1*(PetscScalar)i;
71: MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);
72: }
73: if (j>0) {
74: J = Ii - n3; v = -1.0 + co1*(PetscScalar)j;
75: MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);
76: }
77: if (j<n2-1) {
78: J = Ii + n3; v = -1.0 + co1*(PetscScalar)j;
79: MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);
80: }
81: if (k>0) {
82: J = Ii - 1; v = -1.0 + co1*(PetscScalar)k;
83: MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);
84: }
85: if (k<n3-1) {
86: J = Ii + 1; v = -1.0 + co1*(PetscScalar)k;
87: MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);
88: }
89: v = 6.0 + co2;
90: MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);
91: }
92: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
93: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
95: /* Create parallel vectors and Set right-hand side. */
96: VecCreate(PETSC_COMM_WORLD,&b);
97: VecSetSizes(b,PETSC_DECIDE,n1*n2*n3);
98: VecSetFromOptions(b);
99: VecDuplicate(b,&x);
100: VecDuplicate(b,&u);
101: VecSet(b,1.0);
103: /* Create linear solver context */
104: KSPCreate(PETSC_COMM_WORLD,&ksp);
105: KSPSetOperators(ksp,A,A);
106: KSPSetTolerances(ksp,1.e-6,PETSC_DEFAULT,PETSC_DEFAULT,200);
107: KSPSetFromOptions(ksp);
109: /* Solve the linear system */
110: KSPSolve(ksp,b,x);
112: /* Free work space. */
113: KSPDestroy(&ksp);
114: VecDestroy(&u)); PetscCall(VecDestroy(&x);
115: VecDestroy(&b)); PetscCall(MatDestroy(&A);
116: PetscFinalize();
117: return 0;
118: }
120: /*TEST
122: test:
123: nsize: 8
124: 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
126: test:
127: suffix: 2
128: nsize: 8
129: 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
130: output_file: output/ex39_1.out
132: test:
133: suffix: 3
134: nsize: 8
135: 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
136: output_file: output/ex39_1.out
138: TEST*/