Actual source code: ex42.c
petsc-3.11.4 2019-09-28
2: static char help[] = "Newton's method to solve a two-variable system that comes from the Rosenbrock function.\n\n";
4: /*T
5: Concepts: SNES^basic example
6: T*/
10: /*
11: Include "petscsnes.h" so that we can use SNES solvers. Note that this
12: file automatically includes:
13: petscsys.h - base PETSc routines petscvec.h - vectors
14: petscmat.h - matrices
15: petscis.h - index sets petscksp.h - Krylov subspace methods
16: petscviewer.h - viewers petscpc.h - preconditioners
17: petscksp.h - linear solvers
18: */
19: #include <petscsnes.h>
21: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat,Mat,void*);
22: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
24: int main(int argc,char **argv)
25: {
26: SNES snes; /* nonlinear solver context */
27: Vec x,r; /* solution, residual vectors */
28: Mat J; /* Jacobian matrix */
29: PetscErrorCode ierr;
30: PetscInt its;
31: PetscScalar *xx;
32: SNESConvergedReason reason;
34: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
36: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
37: Create nonlinear solver context
38: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
39: SNESCreate(PETSC_COMM_WORLD,&snes);
41: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
42: Create matrix and vector data structures; set corresponding routines
43: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
44: /*
45: Create vectors for solution and nonlinear function
46: */
47: VecCreate(PETSC_COMM_WORLD,&x);
48: VecSetSizes(x,PETSC_DECIDE,2);
49: VecSetFromOptions(x);
50: VecDuplicate(x,&r);
52: /*
53: Create Jacobian matrix data structure
54: */
55: MatCreate(PETSC_COMM_WORLD,&J);
56: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
57: MatSetFromOptions(J);
58: MatSetUp(J);
60: /*
61: Set function evaluation routine and vector.
62: */
63: SNESSetFunction(snes,r,FormFunction1,NULL);
65: /*
66: Set Jacobian matrix data structure and Jacobian evaluation routine
67: */
68: SNESSetJacobian(snes,J,J,FormJacobian1,NULL);
70: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71: Customize nonlinear solver; set runtime options
72: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
73: SNESSetFromOptions(snes);
75: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76: Evaluate initial guess; then solve nonlinear system
77: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
78: VecGetArray(x,&xx);
79: xx[0] = -1.2; xx[1] = 1.0;
80: VecRestoreArray(x,&xx);
82: /*
83: Note: The user should initialize the vector, x, with the initial guess
84: for the nonlinear solver prior to calling SNESSolve(). In particular,
85: to employ an initial guess of zero, the user should explicitly set
86: this vector to zero by calling VecSet().
87: */
89: SNESSolve(snes,NULL,x);
90: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
91: SNESGetIterationNumber(snes,&its);
92: SNESGetConvergedReason(snes,&reason);
93: /*
94: Some systems computes a residual that is identically zero, thus converging
95: due to FNORM_ABS, others converge due to FNORM_RELATIVE. Here, we only
96: report the reason if the iteration did not converge so that the tests are
97: reproducible.
98: */
99: PetscPrintf(PETSC_COMM_WORLD,"%s number of SNES iterations = %D\n",reason>0 ? "CONVERGED" : (char*) SNESConvergedReasons[reason],its);
101: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: Free work space. All PETSc objects should be destroyed when they
103: are no longer needed.
104: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106: VecDestroy(&x); VecDestroy(&r);
107: MatDestroy(&J); SNESDestroy(&snes);
108: PetscFinalize();
109: return ierr;
110: }
111: /* ------------------------------------------------------------------- */
112: /*
113: FormFunction1 - Evaluates nonlinear function, F(x).
115: Input Parameters:
116: . snes - the SNES context
117: . x - input vector
118: . ctx - optional user-defined context
120: Output Parameter:
121: . f - function vector
122: */
123: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
124: {
125: PetscErrorCode ierr;
126: PetscScalar *ff;
127: const PetscScalar *xx;
129: /*
130: Get pointers to vector data.
131: - For default PETSc vectors, VecGetArray() returns a pointer to
132: the data array. Otherwise, the routine is implementation dependent.
133: - You MUST call VecRestoreArray() when you no longer need access to
134: the array.
135: */
136: VecGetArrayRead(x,&xx);
137: VecGetArray(f,&ff);
139: /* Compute function */
140: ff[0] = -2.0 + 2.0*xx[0] + 400.0*xx[0]*xx[0]*xx[0] - 400.0*xx[0]*xx[1];
141: ff[1] = -200.0*xx[0]*xx[0] + 200.0*xx[1];
143: /* Restore vectors */
144: VecRestoreArrayRead(x,&xx);
145: VecRestoreArray(f,&ff);
146: return 0;
147: }
148: /* ------------------------------------------------------------------- */
149: /*
150: FormJacobian1 - Evaluates Jacobian matrix.
152: Input Parameters:
153: . snes - the SNES context
154: . x - input vector
155: . dummy - optional user-defined context (not used here)
157: Output Parameters:
158: . jac - Jacobian matrix
159: . B - optionally different preconditioning matrix
160: . flag - flag indicating matrix structure
161: */
162: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
163: {
164: const PetscScalar *xx;
165: PetscScalar A[4];
166: PetscErrorCode ierr;
167: PetscInt idx[2] = {0,1};
169: /*
170: Get pointer to vector data
171: */
172: 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;
183: MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES);
185: /*
186: Restore vector
187: */
188: VecRestoreArrayRead(x,&xx);
190: /*
191: Assemble matrix
192: */
193: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
194: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
195: if (jac != B) {
196: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
197: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
198: }
199: return 0;
200: }
204: /*TEST
206: test:
207: args: -snes_monitor_short -snes_max_it 1000
208: requires: !single
210: TEST*/