Actual source code: ex78.c
petsc-3.8.4 2018-03-24
2: static char help[] = "Newton methods to solve u'' = f in parallel with periodic boundary conditions.\n\n";
4: /*T
5: Concepts: SNES^basic parallel example
6: Concepts: periodic boundary conditions
7: Processors: n
8: T*/
10: /*
11: Compare this example to ex3.c that handles Dirichlet boundary conditions
13: Though this is a linear problem it is treated as a nonlinear problem in this example to demonstrate
14: handling periodic boundary conditions for nonlinear problems.
16: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
17: Include "petscsnes.h" so that we can use SNES solvers. Note that this
18: file automatically includes:
19: petscsys.h - base PETSc routines petscvec.h - vectors
20: petscmat.h - matrices
21: petscis.h - index sets petscksp.h - Krylov subspace methods
22: petscviewer.h - viewers petscpc.h - preconditioners
23: petscksp.h - linear solvers
24: */
26: #include <petscdm.h>
27: #include <petscdmda.h>
28: #include <petscsnes.h>
31: PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
32: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
34: int main(int argc,char **argv)
35: {
36: SNES snes; /* SNES context */
37: Mat J; /* Jacobian matrix */
38: DM da;
39: Vec x,r; /* vectors */
41: PetscInt N = 5;
42: MatNullSpace constants;
44: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
45: PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL);
47: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
48: Create nonlinear solver context
49: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
51: SNESCreate(PETSC_COMM_WORLD,&snes);
53: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54: Create vector data structures; set function evaluation routine
55: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
57: /*
58: Create distributed array (DMDA) to manage parallel grid and vectors
59: */
60: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,N,1,1,NULL,&da);
61: DMSetFromOptions(da);
62: DMSetUp(da);
64: /*
65: Extract global and local vectors from DMDA; then duplicate for remaining
66: vectors that are the same types
67: */
68: DMCreateGlobalVector(da,&x);
69: VecDuplicate(x,&r);
71: /*
72: Set function evaluation routine and vector. Whenever the nonlinear
73: solver needs to compute the nonlinear function, it will call this
74: routine.
75: - Note that the final routine argument is the user-defined
76: context that provides application-specific data for the
77: function evaluation routine.
78: */
79: SNESSetFunction(snes,r,FormFunction,da);
81: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82: Create matrix data structure; set Jacobian evaluation routine
83: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
84: DMCreateMatrix(da,&J);
85: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&constants);
86: MatSetNullSpace(J,constants);
87: SNESSetJacobian(snes,J,J,FormJacobian,da);
89: SNESSetFromOptions(snes);
90: SNESSolve(snes,NULL,x);
92: VecDestroy(&x);
93: VecDestroy(&r);
94: MatDestroy(&J);
95: MatNullSpaceDestroy(&constants);
96: SNESDestroy(&snes);
97: DMDestroy(&da);
98: PetscFinalize();
99: return ierr;
100: }
102: /*
103: FormFunction - Evaluates nonlinear function, F(x).
105: Input Parameters:
106: . snes - the SNES context
107: . x - input vector
108: . ctx - optional user-defined context, as set by SNESSetFunction()
110: Output Parameter:
111: . f - function vector
113: Note:
114: The user-defined context can contain any application-specific
115: data needed for the function evaluation.
116: */
117: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
118: {
119: DM da = (DM) ctx;
120: PetscScalar *xx,*ff;
121: PetscReal h;
123: PetscInt i,M,xs,xm;
124: Vec xlocal;
127: /* Get local work vector */
128: DMGetLocalVector(da,&xlocal);
130: /*
131: Scatter ghost points to local vector, using the 2-step process
132: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
133: By placing code between these two statements, computations can
134: be done while messages are in transition.
135: */
136: DMGlobalToLocalBegin(da,x,INSERT_VALUES,xlocal);
137: DMGlobalToLocalEnd(da,x,INSERT_VALUES,xlocal);
139: /*
140: Get pointers to vector data.
141: - The vector xlocal includes ghost point; the vectors x and f do
142: NOT include ghost points.
143: - Using DMDAVecGetArray() allows accessing the values using global ordering
144: */
145: DMDAVecGetArray(da,xlocal,&xx);
146: DMDAVecGetArray(da,f,&ff);
148: /*
149: Get local grid boundaries (for 1-dimensional DMDA):
150: xs, xm - starting grid index, width of local grid (no ghost points)
151: */
152: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
153: DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
155: /*
156: Compute function over locally owned part of the grid
157: Note the [i-1] and [i+1] will automatically access the ghost points from other processes or the periodic points.
158: */
159: h = 1.0/M;
160: for (i=xs; i<xs+xm; i++) ff[i] = (xx[i-1] - 2.0*xx[i] + xx[i+1])/(h*h) - PetscSinReal(2.0*PETSC_PI*i*h);
162: /*
163: Restore vectors
164: */
165: DMDAVecRestoreArray(da,xlocal,&xx);
166: DMDAVecRestoreArray(da,f,&ff);
167: DMRestoreLocalVector(da,&xlocal);
168: return(0);
169: }
170: /* ------------------------------------------------------------------- */
171: /*
172: FormJacobian - Evaluates Jacobian matrix.
174: Input Parameters:
175: . snes - the SNES context
176: . x - input vector
177: . dummy - optional user-defined context (not used here)
179: Output Parameters:
180: . jac - Jacobian matrix
181: . B - optionally different preconditioning matrix
182: . flag - flag indicating matrix structure
183: */
184: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *ctx)
185: {
186: PetscScalar *xx,A[3];
188: PetscInt i,M,xs,xm;
189: DM da = (DM) ctx;
190: MatStencil row,cols[3];
191: PetscReal h;
194: /*
195: Get pointer to vector data
196: */
197: DMDAVecGetArrayRead(da,x,&xx);
198: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
200: /*
201: Get range of locally owned matrix
202: */
203: DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
205: MatZeroEntries(jac);
206: h = 1.0/M;
207: /* because of periodic boundary conditions we can simply loop over all local nodes and access to the left and right */
208: for (i=xs; i<xs+xm; i++) {
209: row.i = i;
210: cols[0].i = i - 1;
211: cols[1].i = i;
212: cols[2].i = i + 1;
213: A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
214: MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
215: }
217: DMDAVecRestoreArrayRead(da,x,&xx);
218: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
219: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
220: return(0);
221: }