Actual source code: ex67.c
petsc-3.8.4 2018-03-24
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*/
10: /*
12: This tests solving singular inconsistent systems with GMRES
14: Default: Solves a symmetric system
15: -symmetric false: Solves a non-symmetric system where nullspace(A) != nullspace(A')
17: -ksp_pc_side left or right
19: See the KSPSolve() for a discussion of when right preconditioning with nullspace(A) != nullspace(A') can fail to produce the
20: norm minimizing solution.
22: Note that though this example does solve the system with right preconditioning and nullspace(A) != nullspace(A') it does not produce the
23: norm minimizing solution, that is the computed solution is not orthogonal to the nullspace(A).
26: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
27: Include "petscksp.h" so that we can use KSP solvers. Note that this
28: file automatically includes:
29: petscsys.h - base PETSc routines petscvec.h - vectors
30: petscmat.h - matrices
31: petscis.h - index sets petscksp.h - Krylov subspace methods
32: petscviewer.h - viewers petscpc.h - preconditioners
33: petscksp.h - linear solvers
34: */
36: #include <petscdm.h>
37: #include <petscdmda.h>
38: #include <petscsnes.h>
39: #include <petsc/private/kspimpl.h>
41: PetscBool symmetric = PETSC_TRUE;
43: PetscErrorCode FormMatrix(Mat,void*);
44: PetscErrorCode FormRightHandSide(Vec,void*);
46: int main(int argc,char **argv)
47: {
48: KSP ksp;
49: Mat J;
50: DM da;
51: Vec x,r; /* vectors */
53: PetscInt M = 10;
54: MatNullSpace constants,nconstants;
56: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
57: PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
58: PetscOptionsGetBool(NULL,NULL,"-symmetric",&symmetric,NULL);
60: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61: Create linear solver context
62: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64: KSPCreate(PETSC_COMM_WORLD,&ksp);
66: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67: Create vector data structures; set function evaluation routine
68: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
70: /*
71: Create distributed array (DMDA) to manage parallel grid and vectors
72: */
73: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,M,1,2,NULL,&da);
74: DMSetFromOptions(da);
75: DMSetUp(da);
77: /*
78: Extract global and local vectors from DMDA; then duplicate for remaining
79: vectors that are the same types
80: */
81: DMCreateGlobalVector(da,&x);
82: VecDuplicate(x,&r);
84: /*
85: Set function evaluation routine and vector. Whenever the nonlinear
86: solver needs to compute the nonlinear function, it will call this
87: routine.
88: - Note that the final routine argument is the user-defined
89: context that provides application-specific data for the
90: function evaluation routine.
91: */
92: FormRightHandSide(r,da);
94: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: Create matrix data structure;
96: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
97: DMCreateMatrix(da,&J);
98: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&constants);
99: if (symmetric) {
100: MatSetOption(J,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
101: MatSetOption(J,MAT_SYMMETRIC,PETSC_TRUE);
102: } else {
103: Vec n;
104: PetscInt zero = 0;
105: PetscScalar zeros = 0.0;
106: VecDuplicate(x,&n);
107: /* the nullspace(A') is the constant vector but with a zero in the very first entry; hence nullspace(A') != nullspace(A) */
108: VecSet(n,1.0);
109: VecSetValues(n,1,&zero,&zeros,INSERT_VALUES);
110: VecAssemblyBegin(n);
111: VecAssemblyEnd(n);
112: VecNormalize(n,NULL);
113: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_FALSE,1,&n,&nconstants);
114: MatSetTransposeNullSpace(J,nconstants);
115: MatNullSpaceDestroy(&nconstants);
116: VecDestroy(&n);
117: }
118: MatSetNullSpace(J,constants);
119: FormMatrix(J,da);
121: KSPSetOperators(ksp,J,J);
123: KSPSetFromOptions(ksp);
124: KSPSolve(ksp,r,x);
125: KSPSolveTranspose(ksp,r,x);
127: VecDestroy(&x);
128: VecDestroy(&r);
129: MatDestroy(&J);
130: MatNullSpaceDestroy(&constants);
131: KSPDestroy(&ksp);
132: DMDestroy(&da);
133: PetscFinalize();
134: return ierr;
135: }
137: /*
139: This intentionally includes something in the right hand side that is not in the range of the matrix A.
140: Since MatSetNullSpace() is called and the matrix is symmetric; the Krylov method automatically removes this
141: portion of the right hand side before solving the linear system.
142: */
143: PetscErrorCode FormRightHandSide(Vec f,void *ctx)
144: {
145: DM da = (DM) ctx;
146: PetscScalar *ff;
148: PetscInt i,M,xs,xm;
149: PetscReal h;
152: DMDAVecGetArray(da,f,&ff);
154: /*
155: Get local grid boundaries (for 1-dimensional DMDA):
156: xs, xm - starting grid index, width of local grid (no ghost points)
157: */
158: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
159: DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
161: /*
162: Compute function over locally owned part of the grid
163: Note the [i-1] and [i+1] will automatically access the ghost points from other processes or the periodic points.
164: */
165: h = 1.0/M;
166: for (i=xs; i<xs+xm; i++) ff[i] = - PetscSinReal(2.0*PETSC_PI*i*h) + 1.0;
168: /*
169: Restore vectors
170: */
171: DMDAVecRestoreArray(da,f,&ff);
172: return(0);
173: }
174: /* ------------------------------------------------------------------- */
175: PetscErrorCode FormMatrix(Mat jac,void *ctx)
176: {
177: PetscScalar A[3];
179: PetscInt i,M,xs,xm;
180: DM da = (DM) ctx;
181: MatStencil row,cols[3];
182: PetscReal h;
185: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
187: /*
188: Get range of locally owned matrix
189: */
190: DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
192: MatZeroEntries(jac);
193: h = 1.0/M;
194: /* because of periodic boundary conditions we can simply loop over all local nodes and access to the left and right */
195: if (symmetric) {
196: for (i=xs; i<xs+xm; i++) {
197: row.i = i;
198: cols[0].i = i - 1;
199: cols[1].i = i;
200: cols[2].i = i + 1;
201: A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
202: MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
203: }
204: } else {
205: MatStencil *acols;
206: PetscScalar *avals;
208: /* only works for one process */
209: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
210: row.i = 0;
211: PetscMalloc1(M,&acols);
212: PetscMalloc1(M,&avals);
213: for (i=0; i<M; i++) {
214: acols[i].i = i;
215: avals[i] = (i % 2) ? 1 : -1;
216: }
217: MatSetValuesStencil(jac,1,&row,M,acols,avals,ADD_VALUES);
218: PetscFree(acols);
219: PetscFree(avals);
220: row.i = 1;
221: cols[0].i = - 1;
222: cols[1].i = 1;
223: cols[2].i = 1 + 1;
224: A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
225: MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
226: for (i=2; i<xs+xm-1; i++) {
227: row.i = i;
228: cols[0].i = i - 1;
229: cols[1].i = i;
230: cols[2].i = i + 1;
231: A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
232: MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
233: }
234: row.i = M - 1 ;
235: cols[0].i = M-2;
236: cols[1].i = M-1;
237: cols[2].i = M+1;
238: A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
239: MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
240: }
241: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
242: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
243: return(0);
244: }