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