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