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