Actual source code: ex78.c
1: static char help[] = "Newton methods to solve u'' = f in parallel with periodic boundary conditions.\n\n";
3: /*
4: Compare this example to ex3.c that handles Dirichlet boundary conditions
6: Though this is a linear problem it is treated as a nonlinear problem in this example to demonstrate
7: handling periodic boundary conditions for nonlinear problems.
9: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
10: Include "petscsnes.h" so that we can use SNES solvers. Note that this
11: file automatically includes:
12: petscsys.h - base PETSc routines petscvec.h - vectors
13: petscmat.h - matrices
14: petscis.h - index sets petscksp.h - Krylov subspace methods
15: petscviewer.h - viewers petscpc.h - preconditioners
16: petscksp.h - linear solvers
17: */
19: #include <petscdm.h>
20: #include <petscdmda.h>
21: #include <petscsnes.h>
23: PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
24: PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
26: int main(int argc, char **argv)
27: {
28: SNES snes; /* SNES context */
29: Mat J; /* Jacobian matrix */
30: DM da;
31: Vec x, r; /* vectors */
32: PetscInt N = 5;
33: MatNullSpace constants;
35: PetscFunctionBeginUser;
36: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
37: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &N, NULL));
39: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
40: Create nonlinear solver context
41: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
43: PetscCall(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: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, N, 1, 1, NULL, &da));
53: PetscCall(DMSetFromOptions(da));
54: PetscCall(DMSetUp(da));
56: /*
57: Extract global and local vectors from DMDA; then duplicate for remaining
58: vectors that are the same types
59: */
60: PetscCall(DMCreateGlobalVector(da, &x));
61: PetscCall(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: PetscCall(SNESSetFunction(snes, r, FormFunction, da));
73: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74: Create matrix data structure; set Jacobian evaluation routine
75: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
76: PetscCall(DMCreateMatrix(da, &J));
77: PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, NULL, &constants));
78: PetscCall(MatSetNullSpace(J, constants));
79: PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, da));
81: PetscCall(SNESSetFromOptions(snes));
82: PetscCall(SNESSolve(snes, NULL, x));
84: PetscCall(VecDestroy(&x));
85: PetscCall(VecDestroy(&r));
86: PetscCall(MatDestroy(&J));
87: PetscCall(MatNullSpaceDestroy(&constants));
88: PetscCall(SNESDestroy(&snes));
89: PetscCall(DMDestroy(&da));
90: PetscCall(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;
117: PetscFunctionBeginUser;
118: /* Get local work vector */
119: PetscCall(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: PetscCall(DMGlobalToLocalBegin(da, x, INSERT_VALUES, xlocal));
128: PetscCall(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: PetscCall(DMDAVecGetArray(da, xlocal, &xx));
137: PetscCall(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: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
144: PetscCall(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: PetscCall(DMDAVecRestoreArray(da, xlocal, &xx));
157: PetscCall(DMDAVecRestoreArray(da, f, &ff));
158: PetscCall(DMRestoreLocalVector(da, &xlocal));
159: PetscFunctionReturn(PETSC_SUCCESS);
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;
183: PetscFunctionBeginUser;
184: /*
185: Get pointer to vector data
186: */
187: PetscCall(DMDAVecGetArrayRead(da, x, &xx));
188: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
190: /*
191: Get range of locally owned matrix
192: */
193: PetscCall(DMDAGetInfo(da, NULL, &M, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
195: PetscCall(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);
204: A[1] = -2.0 / (h * h);
205: PetscCall(MatSetValuesStencil(jac, 1, &row, 3, cols, A, ADD_VALUES));
206: }
208: PetscCall(DMDAVecRestoreArrayRead(da, x, &xx));
209: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
210: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
211: PetscFunctionReturn(PETSC_SUCCESS);
212: }
214: /*TEST
216: test:
217: args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3
218: requires: !single
220: test:
221: suffix: 2
222: args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3 -snes_type newtontrdc
223: requires: !single
225: test:
226: suffix: 3
227: args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3 -snes_type newtontrdc -snes_trdc_use_cauchy false
228: requires: !single
230: TEST*/