Actual source code: ex78.c
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>
30: PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
31: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
33: int main(int argc,char **argv)
34: {
35: SNES snes; /* SNES context */
36: Mat J; /* Jacobian matrix */
37: DM da;
38: Vec x,r; /* vectors */
39: PetscInt N = 5;
40: MatNullSpace constants;
42: PetscInitialize(&argc,&argv,(char*)0,help);
43: PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL);
45: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
46: Create nonlinear solver context
47: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
49: SNESCreate(PETSC_COMM_WORLD,&snes);
51: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52: Create vector data structures; set function evaluation routine
53: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
55: /*
56: Create distributed array (DMDA) to manage parallel grid and vectors
57: */
58: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,N,1,1,NULL,&da);
59: DMSetFromOptions(da);
60: DMSetUp(da);
62: /*
63: Extract global and local vectors from DMDA; then duplicate for remaining
64: vectors that are the same types
65: */
66: DMCreateGlobalVector(da,&x);
67: VecDuplicate(x,&r);
69: /*
70: Set function evaluation routine and vector. Whenever the nonlinear
71: solver needs to compute the nonlinear function, it will call this
72: routine.
73: - Note that the final routine argument is the user-defined
74: context that provides application-specific data for the
75: function evaluation routine.
76: */
77: SNESSetFunction(snes,r,FormFunction,da);
79: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
80: Create matrix data structure; set Jacobian evaluation routine
81: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
82: DMCreateMatrix(da,&J);
83: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&constants);
84: MatSetNullSpace(J,constants);
85: SNESSetJacobian(snes,J,J,FormJacobian,da);
87: SNESSetFromOptions(snes);
88: SNESSolve(snes,NULL,x);
90: VecDestroy(&x);
91: VecDestroy(&r);
92: MatDestroy(&J);
93: MatNullSpaceDestroy(&constants);
94: SNESDestroy(&snes);
95: DMDestroy(&da);
96: PetscFinalize();
97: return 0;
98: }
100: /*
101: FormFunction - Evaluates nonlinear function, F(x).
103: Input Parameters:
104: . snes - the SNES context
105: . x - input vector
106: . ctx - optional user-defined context, as set by SNESSetFunction()
108: Output Parameter:
109: . f - function vector
111: Note:
112: The user-defined context can contain any application-specific
113: data needed for the function evaluation.
114: */
115: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
116: {
117: DM da = (DM) ctx;
118: PetscScalar *xx,*ff;
119: PetscReal h;
120: PetscInt i,M,xs,xm;
121: Vec xlocal;
124: /* Get local work vector */
125: DMGetLocalVector(da,&xlocal);
127: /*
128: Scatter ghost points to local vector, using the 2-step process
129: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
130: By placing code between these two statements, computations can
131: be done while messages are in transition.
132: */
133: DMGlobalToLocalBegin(da,x,INSERT_VALUES,xlocal);
134: DMGlobalToLocalEnd(da,x,INSERT_VALUES,xlocal);
136: /*
137: Get pointers to vector data.
138: - The vector xlocal includes ghost point; the vectors x and f do
139: NOT include ghost points.
140: - Using DMDAVecGetArray() allows accessing the values using global ordering
141: */
142: DMDAVecGetArray(da,xlocal,&xx);
143: DMDAVecGetArray(da,f,&ff);
145: /*
146: Get local grid boundaries (for 1-dimensional DMDA):
147: xs, xm - starting grid index, width of local grid (no ghost points)
148: */
149: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
150: DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
152: /*
153: Compute function over locally owned part of the grid
154: Note the [i-1] and [i+1] will automatically access the ghost points from other processes or the periodic points.
155: */
156: h = 1.0/M;
157: 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);
159: /*
160: Restore vectors
161: */
162: DMDAVecRestoreArray(da,xlocal,&xx);
163: DMDAVecRestoreArray(da,f,&ff);
164: DMRestoreLocalVector(da,&xlocal);
165: return 0;
166: }
167: /* ------------------------------------------------------------------- */
168: /*
169: FormJacobian - Evaluates Jacobian matrix.
171: Input Parameters:
172: . snes - the SNES context
173: . x - input vector
174: . dummy - optional user-defined context (not used here)
176: Output Parameters:
177: . jac - Jacobian matrix
178: . B - optionally different preconditioning matrix
179: . flag - flag indicating matrix structure
180: */
181: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *ctx)
182: {
183: PetscScalar *xx,A[3];
184: PetscInt i,M,xs,xm;
185: DM da = (DM) ctx;
186: MatStencil row,cols[3];
187: PetscReal h;
190: /*
191: Get pointer to vector data
192: */
193: DMDAVecGetArrayRead(da,x,&xx);
194: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
196: /*
197: Get range of locally owned matrix
198: */
199: DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
201: MatZeroEntries(jac);
202: h = 1.0/M;
203: /* because of periodic boundary conditions we can simply loop over all local nodes and access to the left and right */
204: for (i=xs; i<xs+xm; i++) {
205: row.i = i;
206: cols[0].i = i - 1;
207: cols[1].i = i;
208: cols[2].i = i + 1;
209: A[0] = A[2] = 1.0/(h*h); A[1] = -2.0/(h*h);
210: MatSetValuesStencil(jac,1,&row,3,cols,A,ADD_VALUES);
211: }
213: DMDAVecRestoreArrayRead(da,x,&xx);
214: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
215: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
216: return 0;
217: }
219: /*TEST
221: test:
222: args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3
223: requires: !single
225: test:
226: suffix: 2
227: args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3 -snes_type newtontrdc
228: requires: !single
230: test:
231: suffix: 3
232: args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3 -snes_type newtontrdc -snes_trdc_use_cauchy false
233: requires: !single
235: TEST*/