Actual source code: ex67.c
petsc-3.11.4 2019-09-28
2: static char help[] = "Krylov methods to solve u'' = f in parallel with periodic boundary conditions.\n\n";
4: /*T
5: Concepts: KSP^basic parallel example
6: Concepts: periodic boundary conditions
7: Processors: n
8: T*/
12: /*
14: This tests solving singular inconsistent systems with GMRES
16: Default: Solves a symmetric system
17: -symmetric false: Solves a non-symmetric system where nullspace(A) != nullspace(A')
19: -ksp_pc_side left or right
21: See the KSPSolve() for a discussion of when right preconditioning with nullspace(A) != nullspace(A') can fail to produce the
22: norm minimizing solution.
24: Note that though this example does solve the system with right preconditioning and nullspace(A) != nullspace(A') it does not produce the
25: norm minimizing solution, that is the computed solution is not orthogonal to the nullspace(A).
28: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
29: Include "petscksp.h" so that we can use KSP solvers. Note that this
30: file automatically includes:
31: petscsys.h - base PETSc routines petscvec.h - vectors
32: petscmat.h - matrices
33: petscis.h - index sets petscksp.h - Krylov subspace methods
34: petscviewer.h - viewers petscpc.h - preconditioners
35: petscksp.h - linear solvers
36: */
38: #include <petscdm.h>
39: #include <petscdmda.h>
40: #include <petscsnes.h>
41: #include <petsc/private/kspimpl.h>
43: PetscBool symmetric = PETSC_TRUE;
45: PetscErrorCode FormMatrix(Mat,void*);
46: PetscErrorCode FormRightHandSide(Vec,void*);
48: int main(int argc,char **argv)
49: {
50: KSP ksp;
51: Mat J;
52: DM da;
53: Vec x,r; /* vectors */
55: PetscInt M = 10;
56: MatNullSpace constants,nconstants;
58: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
59: PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
60: PetscOptionsGetBool(NULL,NULL,"-symmetric",&symmetric,NULL);
62: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63: Create linear solver context
64: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
66: KSPCreate(PETSC_COMM_WORLD,&ksp);
68: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: Create vector data structures; set function evaluation routine
70: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
72: /*
73: Create distributed array (DMDA) to manage parallel grid and vectors
74: */
75: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,M,1,2,NULL,&da);
76: DMSetFromOptions(da);
77: DMSetUp(da);
79: /*
80: Extract global and local vectors from DMDA; then duplicate for remaining
81: vectors that are the same types
82: */
83: DMCreateGlobalVector(da,&x);
84: VecDuplicate(x,&r);
86: /*
87: Set function evaluation routine and vector. Whenever the nonlinear
88: solver needs to compute the nonlinear function, it will call this
89: routine.
90: - Note that the final routine argument is the user-defined
91: context that provides application-specific data for the
92: function evaluation routine.
93: */
94: FormRightHandSide(r,da);
96: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: Create matrix data structure;
98: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99: DMCreateMatrix(da,&J);
100: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&constants);
101: if (symmetric) {
102: MatSetOption(J,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
103: MatSetOption(J,MAT_SYMMETRIC,PETSC_TRUE);
104: } else {
105: Vec n;
106: PetscInt zero = 0;
107: PetscScalar zeros = 0.0;
108: VecDuplicate(x,&n);
109: /* the nullspace(A') is the constant vector but with a zero in the very first entry; hence nullspace(A') != nullspace(A) */
110: VecSet(n,1.0);
111: VecSetValues(n,1,&zero,&zeros,INSERT_VALUES);
112: VecAssemblyBegin(n);
113: VecAssemblyEnd(n);
114: VecNormalize(n,NULL);
115: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_FALSE,1,&n,&nconstants);
116: MatSetTransposeNullSpace(J,nconstants);
117: MatNullSpaceDestroy(&nconstants);
118: VecDestroy(&n);
119: }
120: MatSetNullSpace(J,constants);
121: FormMatrix(J,da);
123: KSPSetOperators(ksp,J,J);
125: KSPSetFromOptions(ksp);
126: KSPSolve(ksp,r,x);
127: KSPSolveTranspose(ksp,r,x);
129: VecDestroy(&x);
130: VecDestroy(&r);
131: MatDestroy(&J);
132: MatNullSpaceDestroy(&constants);
133: KSPDestroy(&ksp);
134: DMDestroy(&da);
135: PetscFinalize();
136: return ierr;
137: }
139: /*
141: This intentionally includes something in the right hand side that is not in the range of the matrix A.
142: Since MatSetNullSpace() is called and the matrix is symmetric; the Krylov method automatically removes this
143: portion of the right hand side before solving the linear system.
144: */
145: PetscErrorCode FormRightHandSide(Vec f,void *ctx)
146: {
147: DM da = (DM) ctx;
148: PetscScalar *ff;
150: PetscInt i,M,xs,xm;
151: PetscReal h;
154: DMDAVecGetArray(da,f,&ff);
156: /*
157: Get local grid boundaries (for 1-dimensional DMDA):
158: xs, xm - starting grid index, width of local grid (no ghost points)
159: */
160: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
161: DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
163: /*
164: Compute function over locally owned part of the grid
165: Note the [i-1] and [i+1] will automatically access the ghost points from other processes or the periodic points.
166: */
167: h = 1.0/M;
168: for (i=xs; i<xs+xm; i++) ff[i] = - PetscSinReal(2.0*PETSC_PI*i*h) + 1.0;
170: /*
171: Restore vectors
172: */
173: DMDAVecRestoreArray(da,f,&ff);
174: return(0);
175: }
176: /* ------------------------------------------------------------------- */
177: PetscErrorCode FormMatrix(Mat jac,void *ctx)
178: {
179: PetscScalar A[3];
181: PetscInt i,M,xs,xm;
182: DM da = (DM) ctx;
183: MatStencil row,cols[3];
184: PetscReal h;
187: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
189: /*
190: Get range of locally owned matrix
191: */
192: DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
194: MatZeroEntries(jac);
195: h = 1.0/M;
196: /* because of periodic boundary conditions we can simply loop over all local nodes and access to the left and right */
197: if (symmetric) {
198: for (i=xs; i<xs+xm; i++) {
199: row.i = i;
200: cols[0].i = i - 1;
201: cols[1].i = i;
202: cols[2].i = i + 1;
203: A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
204: MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
205: }
206: } else {
207: MatStencil *acols;
208: PetscScalar *avals;
210: /* only works for one process */
211: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
212: row.i = 0;
213: PetscMalloc1(M,&acols);
214: PetscMalloc1(M,&avals);
215: for (i=0; i<M; i++) {
216: acols[i].i = i;
217: avals[i] = (i % 2) ? 1 : -1;
218: }
219: MatSetValuesStencil(jac,1,&row,M,acols,avals,ADD_VALUES);
220: PetscFree(acols);
221: PetscFree(avals);
222: row.i = 1;
223: cols[0].i = - 1;
224: cols[1].i = 1;
225: cols[2].i = 1 + 1;
226: A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
227: MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
228: for (i=2; i<xs+xm-1; i++) {
229: row.i = i;
230: cols[0].i = i - 1;
231: cols[1].i = i;
232: cols[2].i = i + 1;
233: A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
234: MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
235: }
236: row.i = M - 1 ;
237: cols[0].i = M-2;
238: cols[1].i = M-1;
239: cols[2].i = M+1;
240: A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
241: MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
242: }
243: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
244: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
245: return(0);
246: }
249: /*TEST
251: test:
252: suffix: nonsymmetric_left
253: args: -symmetric false -ksp_view -ksp_converged_reason -pc_type jacobi -mat_no_inode -ksp_monitor_true_residual -ksp_rtol 1.e-14 -ksp_max_it 12 -ksp_pc_side left
254: filter: sed 's/ATOL/RTOL/g'
255: requires: !single
257: test:
258: suffix: nonsymmetric_right
259: args: -symmetric false -ksp_view -ksp_converged_reason -pc_type jacobi -mat_no_inode -ksp_monitor_true_residual -ksp_rtol 1.e-14 -ksp_max_it 12 -ksp_pc_side right
260: filter: sed 's/ATOL/RTOL/g'
261: requires: !single
263: test:
264: suffix: symmetric_left
265: args: -ksp_view -ksp_converged_reason -pc_type sor -mat_no_inode -ksp_monitor_true_residual -ksp_rtol 1.e-14 -ksp_max_it 12 -ksp_pc_side left
266: requires: !single
268: test:
269: suffix: symmetric_right
270: args: -ksp_view -ksp_converged_reason -pc_type sor -mat_no_inode -ksp_monitor_true_residual -ksp_rtol 1.e-14 -ksp_max_it 12 -ksp_pc_side right
271: filter: sed 's/ATOL/RTOL/g'
272: requires: !single
274: test:
275: suffix: transpose_asm
276: args: -symmetric false -ksp_monitor -ksp_view -pc_type asm -sub_pc_type lu -sub_pc_factor_zeropivot 1.e-33 -ksp_converged_reason
277: filter: sed 's/ATOL/RTOL/g'
279: TEST*/