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 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>
14: extern PetscErrorCode FormJacobian1(SNES, Vec, Mat, Mat, void *);
15: extern PetscErrorCode FormFunction1(SNES, Vec, Vec, void *);
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;
27: PetscFunctionBeginUser;
28: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
29: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_ghost", &test_ghost, NULL));
31: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
32: Create nonlinear solver context
33: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
34: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
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;
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: }
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));
62: /*
63: Set function evaluation routine and vector.
64: */
65: PetscCall(SNESSetFunction(snes, NULL, FormFunction1, &test_ghost));
67: /*
68: Set Jacobian matrix data structure and Jacobian evaluation routine
69: */
70: PetscCall(SNESSetJacobian(snes, J, J, FormJacobian1, &test_ghost));
72: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73: Customize nonlinear solver; set runtime options
74: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
75: PetscCall(SNESSetFromOptions(snes));
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));
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: */
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));
104: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105: Free work space. All PETSc objects should be destroyed when they
106: are no longer needed.
107: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: PetscCall(VecDestroy(&x));
110: PetscCall(MatDestroy(&J));
111: PetscCall(SNESDestroy(&snes));
112: PetscCall(PetscFinalize());
113: return 0;
114: }
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: }
128: PetscErrorCode FormFunction1(SNES snes, Vec x, Vec f, void *ctx)
129: {
130: PetscScalar *ff;
131: const PetscScalar *xx;
133: PetscFunctionBeginUser;
134: if (*(PetscBool *)ctx) {
135: PetscCall(VecCheckGhosted(x, PETSC_FALSE));
136: PetscCall(VecCheckGhosted(f, PETSC_TRUE));
137: }
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));
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];
153: /* Restore vectors */
154: PetscCall(VecRestoreArrayRead(x, &xx));
155: PetscCall(VecRestoreArray(f, &ff));
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
159: PetscErrorCode FormJacobian1(SNES snes, Vec x, Mat jac, Mat B, void *ctx)
160: {
161: const PetscScalar *xx;
162: PetscScalar A[4];
163: PetscInt idx[2];
164: PetscMPIInt rank;
166: PetscFunctionBeginUser;
167: if (*(PetscBool *)ctx) { PetscCall(VecCheckGhosted(x, PETSC_FALSE)); }
168: /*
169: Get pointer to vector data
170: */
171: PetscCall(VecGetArrayRead(x, &xx));
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;
183: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)x), &rank));
184: idx[0] = 2 * rank;
185: idx[1] = 2 * rank + 1;
187: PetscCall(MatSetValues(B, 2, idx, 2, idx, A, INSERT_VALUES));
189: /*
190: Restore vector
191: */
192: PetscCall(VecRestoreArrayRead(x, &xx));
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: }
206: /*TEST
208: test:
209: suffix: 1
210: args: -snes_monitor_short -snes_max_it 1000 -sol_view
211: requires: !single
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
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
224: TEST*/