Actual source code: ex67.c

petsc-3.11.3 2019-06-26
Report Typos and Errors

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