1 static char help[] = "Newton methods to solve u'' = f in parallel with periodic boundary conditions.\n\n";
2
3 /*
4 Compare this example to ex3.c that handles Dirichlet boundary conditions
5
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.
8
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 */
18
19 #include <petscdm.h>
20 #include <petscdmda.h>
21 #include <petscsnes.h>
22
23 PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
24 PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
25
main(int argc,char ** argv)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;
34
35 PetscFunctionBeginUser;
36 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
37 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &N, NULL));
38
39 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
40 Create nonlinear solver context
41 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
42
43 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
44
45 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
46 Create vector data structures; set function evaluation routine
47 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
48
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));
55
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));
62
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));
72
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));
80
81 PetscCall(SNESSetFromOptions(snes));
82 PetscCall(SNESSolve(snes, NULL, x));
83
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 }
93
94 /*
95 FormFunction - Evaluates nonlinear function, F(x).
96
97 Input Parameters:
98 . snes - the SNES context
99 . x - input vector
100 . ctx - optional user-defined context, as set by SNESSetFunction()
101
102 Output Parameter:
103 . f - function vector
104
105 Note:
106 The user-defined context can contain any application-specific
107 data needed for the function evaluation.
108 */
FormFunction(SNES snes,Vec x,Vec f,PetscCtx ctx)109 PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, PetscCtx ctx)
110 {
111 DM da = (DM)ctx;
112 PetscScalar *xx, *ff;
113 PetscReal h;
114 PetscInt i, M, xs, xm;
115 Vec xlocal;
116
117 PetscFunctionBeginUser;
118 /* Get local work vector */
119 PetscCall(DMGetLocalVector(da, &xlocal));
120
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));
129
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));
138
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));
145
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);
152
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.
164
165 Input Parameters:
166 . snes - the SNES context
167 . x - input vector
168 . dummy - optional user-defined context (not used here)
169
170 Output Parameters:
171 . jac - Jacobian matrix
172 . B - optionally different matrix used to construct the preconditioner
173 */
FormJacobian(SNES snes,Vec x,Mat jac,Mat B,PetscCtx ctx)174 PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, PetscCtx ctx)
175 {
176 PetscScalar *xx, A[3];
177 PetscInt i, M, xs, xm;
178 DM da = (DM)ctx;
179 MatStencil row, cols[3];
180 PetscReal h;
181
182 PetscFunctionBeginUser;
183 /*
184 Get pointer to vector data
185 */
186 PetscCall(DMDAVecGetArrayRead(da, x, &xx));
187 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
188
189 /*
190 Get range of locally owned matrix
191 */
192 PetscCall(DMDAGetInfo(da, NULL, &M, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
193
194 PetscCall(MatZeroEntries(jac));
195 h = 1.0 / M;
196 /* because of periodic boundary conditions we can simply loop over all local nodes and access to the left and right */
197 for (i = xs; i < xs + xm; i++) {
198 row.i = i;
199 cols[0].i = i - 1;
200 cols[1].i = i;
201 cols[2].i = i + 1;
202 A[0] = A[2] = 1.0 / (h * h);
203 A[1] = -2.0 / (h * h);
204 PetscCall(MatSetValuesStencil(jac, 1, &row, 3, cols, A, ADD_VALUES));
205 }
206
207 PetscCall(DMDAVecRestoreArrayRead(da, x, &xx));
208 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
209 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
210 PetscFunctionReturn(PETSC_SUCCESS);
211 }
212
213 /*TEST
214
215 test:
216 args: -snes_monitor_short -ksp_monitor_short -pc_type sor -snes_converged_reason -da_refine 3
217 requires: !single
218
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
223
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
228
229 TEST*/
230