Actual source code: ex78.c
petsc-3.9.4 2018-09-11
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*/
12: /*
13: Compare this example to ex3.c that handles Dirichlet boundary conditions
15: Though this is a linear problem it is treated as a nonlinear problem in this example to demonstrate
16: handling periodic boundary conditions for nonlinear problems.
18: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
19: Include "petscsnes.h" so that we can use SNES solvers. Note that this
20: file automatically includes:
21: petscsys.h - base PETSc routines petscvec.h - vectors
22: petscmat.h - matrices
23: petscis.h - index sets petscksp.h - Krylov subspace methods
24: petscviewer.h - viewers petscpc.h - preconditioners
25: petscksp.h - linear solvers
26: */
28: #include <petscdm.h>
29: #include <petscdmda.h>
30: #include <petscsnes.h>
33: PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
34: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
36: int main(int argc,char **argv)
37: {
38: SNES snes; /* SNES context */
39: Mat J; /* Jacobian matrix */
40: DM da;
41: Vec x,r; /* vectors */
43: PetscInt N = 5;
44: MatNullSpace constants;
46: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
47: PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL);
49: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: Create nonlinear solver context
51: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
53: SNESCreate(PETSC_COMM_WORLD,&snes);
55: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56: Create vector data structures; set function evaluation routine
57: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59: /*
60: Create distributed array (DMDA) to manage parallel grid and vectors
61: */
62: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,N,1,1,NULL,&da);
63: DMSetFromOptions(da);
64: DMSetUp(da);
66: /*
67: Extract global and local vectors from DMDA; then duplicate for remaining
68: vectors that are the same types
69: */
70: DMCreateGlobalVector(da,&x);
71: VecDuplicate(x,&r);
73: /*
74: Set function evaluation routine and vector. Whenever the nonlinear
75: solver needs to compute the nonlinear function, it will call this
76: routine.
77: - Note that the final routine argument is the user-defined
78: context that provides application-specific data for the
79: function evaluation routine.
80: */
81: SNESSetFunction(snes,r,FormFunction,da);
83: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: Create matrix data structure; set Jacobian evaluation routine
85: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
86: DMCreateMatrix(da,&J);
87: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&constants);
88: MatSetNullSpace(J,constants);
89: SNESSetJacobian(snes,J,J,FormJacobian,da);
91: SNESSetFromOptions(snes);
92: SNESSolve(snes,NULL,x);
94: VecDestroy(&x);
95: VecDestroy(&r);
96: MatDestroy(&J);
97: MatNullSpaceDestroy(&constants);
98: SNESDestroy(&snes);
99: DMDestroy(&da);
100: PetscFinalize();
101: return ierr;
102: }
104: /*
105: FormFunction - Evaluates nonlinear function, F(x).
107: Input Parameters:
108: . snes - the SNES context
109: . x - input vector
110: . ctx - optional user-defined context, as set by SNESSetFunction()
112: Output Parameter:
113: . f - function vector
115: Note:
116: The user-defined context can contain any application-specific
117: data needed for the function evaluation.
118: */
119: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
120: {
121: DM da = (DM) ctx;
122: PetscScalar *xx,*ff;
123: PetscReal h;
125: PetscInt i,M,xs,xm;
126: Vec xlocal;
129: /* Get local work vector */
130: DMGetLocalVector(da,&xlocal);
132: /*
133: Scatter ghost points to local vector, using the 2-step process
134: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
135: By placing code between these two statements, computations can
136: be done while messages are in transition.
137: */
138: DMGlobalToLocalBegin(da,x,INSERT_VALUES,xlocal);
139: DMGlobalToLocalEnd(da,x,INSERT_VALUES,xlocal);
141: /*
142: Get pointers to vector data.
143: - The vector xlocal includes ghost point; the vectors x and f do
144: NOT include ghost points.
145: - Using DMDAVecGetArray() allows accessing the values using global ordering
146: */
147: DMDAVecGetArray(da,xlocal,&xx);
148: DMDAVecGetArray(da,f,&ff);
150: /*
151: Get local grid boundaries (for 1-dimensional DMDA):
152: xs, xm - starting grid index, width of local grid (no ghost points)
153: */
154: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
155: DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
157: /*
158: Compute function over locally owned part of the grid
159: Note the [i-1] and [i+1] will automatically access the ghost points from other processes or the periodic points.
160: */
161: h = 1.0/M;
162: 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);
164: /*
165: Restore vectors
166: */
167: DMDAVecRestoreArray(da,xlocal,&xx);
168: DMDAVecRestoreArray(da,f,&ff);
169: DMRestoreLocalVector(da,&xlocal);
170: return(0);
171: }
172: /* ------------------------------------------------------------------- */
173: /*
174: FormJacobian - Evaluates Jacobian matrix.
176: Input Parameters:
177: . snes - the SNES context
178: . x - input vector
179: . dummy - optional user-defined context (not used here)
181: Output Parameters:
182: . jac - Jacobian matrix
183: . B - optionally different preconditioning matrix
184: . flag - flag indicating matrix structure
185: */
186: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *ctx)
187: {
188: PetscScalar *xx,A[3];
190: PetscInt i,M,xs,xm;
191: DM da = (DM) ctx;
192: MatStencil row,cols[3];
193: PetscReal h;
196: /*
197: Get pointer to vector data
198: */
199: DMDAVecGetArrayRead(da,x,&xx);
200: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
202: /*
203: Get range of locally owned matrix
204: */
205: DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
207: MatZeroEntries(jac);
208: h = 1.0/M;
209: /* because of periodic boundary conditions we can simply loop over all local nodes and access to the left and right */
210: for (i=xs; i<xs+xm; i++) {
211: row.i = i;
212: cols[0].i = i - 1;
213: cols[1].i = i;
214: cols[2].i = i + 1;
215: A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
216: MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
217: }
219: DMDAVecRestoreArrayRead(da,x,&xx);
220: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
221: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
222: return(0);
223: }
226: /*TEST
228: test:
229: args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3
230: requires: !single
232: TEST*/