Actual source code: ex42.c
petsc-3.8.4 2018-03-24
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*/
8: /*
9: Include "petscsnes.h" so that we can use SNES solvers. Note that this
10: file automatically includes:
11: petscsys.h - base PETSc routines petscvec.h - vectors
12: petscmat.h - matrices
13: petscis.h - index sets petscksp.h - Krylov subspace methods
14: petscviewer.h - viewers petscpc.h - preconditioners
15: petscksp.h - linear solvers
16: */
17: #include <petscsnes.h>
19: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat,Mat,void*);
20: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
22: int main(int argc,char **argv)
23: {
24: SNES snes; /* nonlinear solver context */
25: Vec x,r; /* solution, residual vectors */
26: Mat J; /* Jacobian matrix */
27: PetscErrorCode ierr;
28: PetscInt its;
29: PetscScalar *xx;
30: SNESConvergedReason reason;
32: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
34: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
35: Create nonlinear solver context
36: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
37: SNESCreate(PETSC_COMM_WORLD,&snes);
39: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
40: Create matrix and vector data structures; set corresponding routines
41: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
42: /*
43: Create vectors for solution and nonlinear function
44: */
45: VecCreate(PETSC_COMM_WORLD,&x);
46: VecSetSizes(x,PETSC_DECIDE,2);
47: VecSetFromOptions(x);
48: VecDuplicate(x,&r);
50: /*
51: Create Jacobian matrix data structure
52: */
53: MatCreate(PETSC_COMM_WORLD,&J);
54: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
55: MatSetFromOptions(J);
56: MatSetUp(J);
58: /*
59: Set function evaluation routine and vector.
60: */
61: SNESSetFunction(snes,r,FormFunction1,NULL);
63: /*
64: Set Jacobian matrix data structure and Jacobian evaluation routine
65: */
66: SNESSetJacobian(snes,J,J,FormJacobian1,NULL);
68: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: Customize nonlinear solver; set runtime options
70: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
71: SNESSetFromOptions(snes);
73: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74: Evaluate initial guess; then solve nonlinear system
75: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
76: VecGetArray(x,&xx);
77: xx[0] = -1.2; xx[1] = 1.0;
78: VecRestoreArray(x,&xx);
80: /*
81: Note: The user should initialize the vector, x, with the initial guess
82: for the nonlinear solver prior to calling SNESSolve(). In particular,
83: to employ an initial guess of zero, the user should explicitly set
84: this vector to zero by calling VecSet().
85: */
87: SNESSolve(snes,NULL,x);
88: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
89: SNESGetIterationNumber(snes,&its);
90: SNESGetConvergedReason(snes,&reason);
91: /*
92: Some systems computes a residual that is identically zero, thus converging
93: due to FNORM_ABS, others converge due to FNORM_RELATIVE. Here, we only
94: report the reason if the iteration did not converge so that the tests are
95: reproducible.
96: */
97: PetscPrintf(PETSC_COMM_WORLD,"%s number of SNES iterations = %D\n",reason>0 ? "CONVERGED" : (char*) SNESConvergedReasons[reason],its);
99: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: Free work space. All PETSc objects should be destroyed when they
101: are no longer needed.
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104: VecDestroy(&x); VecDestroy(&r);
105: MatDestroy(&J); SNESDestroy(&snes);
106: PetscFinalize();
107: return ierr;
108: }
109: /* ------------------------------------------------------------------- */
110: /*
111: FormFunction1 - Evaluates nonlinear function, F(x).
113: Input Parameters:
114: . snes - the SNES context
115: . x - input vector
116: . ctx - optional user-defined context
118: Output Parameter:
119: . f - function vector
120: */
121: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
122: {
123: PetscErrorCode ierr;
124: PetscScalar *ff;
125: const PetscScalar *xx;
127: /*
128: Get pointers to vector data.
129: - For default PETSc vectors, VecGetArray() returns a pointer to
130: the data array. Otherwise, the routine is implementation dependent.
131: - You MUST call VecRestoreArray() when you no longer need access to
132: the array.
133: */
134: VecGetArrayRead(x,&xx);
135: VecGetArray(f,&ff);
137: /* Compute function */
138: ff[0] = -2.0 + 2.0*xx[0] + 400.0*xx[0]*xx[0]*xx[0] - 400.0*xx[0]*xx[1];
139: ff[1] = -200.0*xx[0]*xx[0] + 200.0*xx[1];
141: /* Restore vectors */
142: VecRestoreArrayRead(x,&xx);
143: VecRestoreArray(f,&ff);
144: return 0;
145: }
146: /* ------------------------------------------------------------------- */
147: /*
148: FormJacobian1 - Evaluates Jacobian matrix.
150: Input Parameters:
151: . snes - the SNES context
152: . x - input vector
153: . dummy - optional user-defined context (not used here)
155: Output Parameters:
156: . jac - Jacobian matrix
157: . B - optionally different preconditioning matrix
158: . flag - flag indicating matrix structure
159: */
160: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
161: {
162: const PetscScalar *xx;
163: PetscScalar A[4];
164: PetscErrorCode ierr;
165: PetscInt idx[2] = {0,1};
167: /*
168: Get pointer to vector data
169: */
170: VecGetArrayRead(x,&xx);
172: /*
173: Compute Jacobian entries and insert into matrix.
174: - Since this is such a small problem, we set all entries for
175: the matrix at once.
176: */
177: A[0] = 2.0 + 1200.0*xx[0]*xx[0] - 400.0*xx[1];
178: A[1] = -400.0*xx[0];
179: A[2] = -400.0*xx[0];
180: A[3] = 200;
181: MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES);
183: /*
184: Restore vector
185: */
186: VecRestoreArrayRead(x,&xx);
188: /*
189: Assemble matrix
190: */
191: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
192: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
193: if (jac != B) {
194: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
195: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
196: }
197: return 0;
198: }