1 static char help[] = "Newton's method to solve a two-variable system that comes from the Rosenbrock function.\n\n";
2
3 /*
4 Include "petscsnes.h" so that we can use SNES solvers. Note that this
5 file automatically includes:
6 petscsys.h - base PETSc routines petscvec.h - vectors
7 petscmat.h - matrices
8 petscis.h - index sets petscksp.h - Krylov subspace methods
9 petscviewer.h - viewers petscpc.h - preconditioners
10 petscksp.h - linear solvers
11 */
12 #include <petscsnes.h>
13
14 extern PetscErrorCode FormJacobian1(SNES, Vec, Mat, Mat, void *);
15 extern PetscErrorCode FormFunction1(SNES, Vec, Vec, void *);
16
main(int argc,char ** argv)17 int main(int argc, char **argv)
18 {
19 SNES snes; /* nonlinear solver context */
20 Vec x; /* solution vector */
21 Mat J; /* Jacobian matrix */
22 PetscInt its;
23 PetscScalar *xx;
24 SNESConvergedReason reason;
25 PetscBool test_ghost = PETSC_FALSE;
26
27 PetscFunctionBeginUser;
28 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
29 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_ghost", &test_ghost, NULL));
30
31 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
32 Create nonlinear solver context
33 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
34 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
35
36 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
37 Create matrix and vector data structures; set corresponding routines
38 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
39 /*
40 Create vectors for solution and nonlinear function
41 */
42 if (test_ghost) {
43 PetscInt gIdx[] = {0, 1};
44 PetscInt nghost = 2;
45
46 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nghost", &nghost, NULL));
47 PetscCall(VecCreateGhost(PETSC_COMM_WORLD, 2, PETSC_DECIDE, PetscMin(nghost, 2), gIdx, &x));
48 } else {
49 PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
50 PetscCall(VecSetSizes(x, 2, PETSC_DECIDE));
51 PetscCall(VecSetFromOptions(x));
52 }
53
54 /*
55 Create Jacobian matrix data structure
56 */
57 PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
58 PetscCall(MatSetSizes(J, 2, 2, PETSC_DECIDE, PETSC_DECIDE));
59 PetscCall(MatSetFromOptions(J));
60 PetscCall(MatSetUp(J));
61
62 /*
63 Set function evaluation routine and vector.
64 */
65 PetscCall(SNESSetFunction(snes, NULL, FormFunction1, &test_ghost));
66
67 /*
68 Set Jacobian matrix data structure and Jacobian evaluation routine
69 */
70 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian1, &test_ghost));
71
72 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73 Customize nonlinear solver; set runtime options
74 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
75 PetscCall(SNESSetFromOptions(snes));
76
77 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78 Evaluate initial guess; then solve nonlinear system
79 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
80 PetscCall(VecGetArray(x, &xx));
81 xx[0] = -1.2;
82 xx[1] = 1.0;
83 PetscCall(VecRestoreArray(x, &xx));
84
85 /*
86 Note: The user should initialize the vector, x, with the initial guess
87 for the nonlinear solver prior to calling SNESSolve(). In particular,
88 to employ an initial guess of zero, the user should explicitly set
89 this vector to zero by calling VecSet().
90 */
91
92 PetscCall(SNESSolve(snes, NULL, x));
93 PetscCall(VecViewFromOptions(x, NULL, "-sol_view"));
94 PetscCall(SNESGetIterationNumber(snes, &its));
95 PetscCall(SNESGetConvergedReason(snes, &reason));
96 /*
97 Some systems computes a residual that is identically zero, thus converging
98 due to FNORM_ABS, others converge due to FNORM_RELATIVE. Here, we only
99 report the reason if the iteration did not converge so that the tests are
100 reproducible.
101 */
102 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s number of SNES iterations = %" PetscInt_FMT "\n", reason > 0 ? "CONVERGED" : SNESConvergedReasons[reason], its));
103
104 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105 Free work space. All PETSc objects should be destroyed when they
106 are no longer needed.
107 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108
109 PetscCall(VecDestroy(&x));
110 PetscCall(MatDestroy(&J));
111 PetscCall(SNESDestroy(&snes));
112 PetscCall(PetscFinalize());
113 return 0;
114 }
115
VecCheckGhosted(Vec X,PetscBool test_rev)116 PetscErrorCode VecCheckGhosted(Vec X, PetscBool test_rev)
117 {
118 PetscFunctionBeginUser;
119 PetscCall(VecGhostUpdateBegin(X, INSERT_VALUES, SCATTER_FORWARD));
120 PetscCall(VecGhostUpdateEnd(X, INSERT_VALUES, SCATTER_FORWARD));
121 if (test_rev) {
122 PetscCall(VecGhostUpdateBegin(X, INSERT_VALUES, SCATTER_REVERSE));
123 PetscCall(VecGhostUpdateEnd(X, INSERT_VALUES, SCATTER_REVERSE));
124 }
125 PetscFunctionReturn(PETSC_SUCCESS);
126 }
127
FormFunction1(SNES snes,Vec x,Vec f,PetscCtx ctx)128 PetscErrorCode FormFunction1(SNES snes, Vec x, Vec f, PetscCtx ctx)
129 {
130 PetscScalar *ff;
131 const PetscScalar *xx;
132
133 PetscFunctionBeginUser;
134 if (*(PetscBool *)ctx) {
135 PetscCall(VecCheckGhosted(x, PETSC_FALSE));
136 PetscCall(VecCheckGhosted(f, PETSC_TRUE));
137 }
138
139 /*
140 Get pointers to vector data.
141 - For default PETSc vectors, VecGetArray() returns a pointer to
142 the data array. Otherwise, the routine is implementation dependent.
143 - You MUST call VecRestoreArray() when you no longer need access to
144 the array.
145 */
146 PetscCall(VecGetArrayRead(x, &xx));
147 PetscCall(VecGetArray(f, &ff));
148
149 /* Compute function */
150 ff[0] = -2.0 + 2.0 * xx[0] + 400.0 * xx[0] * xx[0] * xx[0] - 400.0 * xx[0] * xx[1];
151 ff[1] = -200.0 * xx[0] * xx[0] + 200.0 * xx[1];
152
153 /* Restore vectors */
154 PetscCall(VecRestoreArrayRead(x, &xx));
155 PetscCall(VecRestoreArray(f, &ff));
156 PetscFunctionReturn(PETSC_SUCCESS);
157 }
158
FormJacobian1(SNES snes,Vec x,Mat jac,Mat B,PetscCtx ctx)159 PetscErrorCode FormJacobian1(SNES snes, Vec x, Mat jac, Mat B, PetscCtx ctx)
160 {
161 const PetscScalar *xx;
162 PetscScalar A[4];
163 PetscInt idx[2];
164 PetscMPIInt rank;
165
166 PetscFunctionBeginUser;
167 if (*(PetscBool *)ctx) PetscCall(VecCheckGhosted(x, PETSC_FALSE));
168 /*
169 Get pointer to vector data
170 */
171 PetscCall(VecGetArrayRead(x, &xx));
172
173 /*
174 Compute Jacobian entries and insert into matrix.
175 - Since this is such a small problem, we set all entries for
176 the matrix at once.
177 */
178 A[0] = 2.0 + 1200.0 * xx[0] * xx[0] - 400.0 * xx[1];
179 A[1] = -400.0 * xx[0];
180 A[2] = -400.0 * xx[0];
181 A[3] = 200;
182
183 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)x), &rank));
184 idx[0] = 2 * rank;
185 idx[1] = 2 * rank + 1;
186
187 PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES));
188
189 /*
190 Restore vector
191 */
192 PetscCall(VecRestoreArrayRead(x, &xx));
193
194 /*
195 Assemble matrix
196 */
197 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
198 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
199 if (jac != B) {
200 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
201 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
202 }
203 PetscFunctionReturn(PETSC_SUCCESS);
204 }
205
206 /*TEST
207
208 test:
209 suffix: 1
210 args: -snes_monitor_short -snes_max_it 1000 -sol_view
211 requires: !single
212
213 test:
214 suffix: 2
215 args: -snes_monitor_short -snes_max_it 1000 -snes_type newtontrdc -snes_trdc_use_cauchy false -sol_view
216 requires: !single
217
218 test:
219 suffix: ghosts
220 nsize: {{1 2}}
221 args: -snes_max_it 4 -snes_type {{newtontr newtonls}} -nghost {{0 1 2}} -test_ghost
222 requires: !single
223
224 TEST*/
225