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