Actual source code: ex1.c


  2: static char help[] = "Newton's method for a two-variable system, sequential.\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: */
This examples solves either
\begin{equation}
F\genfrac{(}{)}{0pt}{}{x_0}{x_1} = \genfrac{(}{)}{0pt}{}{x^2_0 + x_0 x_1 - 3}{x_0 x_1 + x^2_1 - 6}
\end{equation}
or if the {\tt -hard} options is given
\begin{equation}
F\genfrac{(}{)}{0pt}{}{x_0}{x_1} = \genfrac{(}{)}{0pt}{}{\sin(3 x_0) + x_0}{x_1}
\end{equation}
 27: #include <petscsnes.h>

 29: /*
 30:    User-defined routines
 31: */
 32: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat,Mat,void*);
 33: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
 34: extern PetscErrorCode FormJacobian2(SNES,Vec,Mat,Mat,void*);
 35: extern PetscErrorCode FormFunction2(SNES,Vec,Vec,void*);

 37: int main(int argc,char **argv)
 38: {
 39:   SNES           snes;         /* nonlinear solver context */
 40:   KSP            ksp;          /* linear solver context */
 41:   PC             pc;           /* preconditioner context */
 42:   Vec            x,r;          /* solution, residual vectors */
 43:   Mat            J;            /* Jacobian matrix */
 44:   PetscMPIInt    size;
 45:   PetscScalar    pfive = .5,*xx;
 46:   PetscBool      flg;

 48:   PetscInitialize(&argc,&argv,(char*)0,help);
 49:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

 52:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 53:      Create nonlinear solver context
 54:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 55:   SNESCreate(PETSC_COMM_WORLD,&snes);

 57:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 58:      Create matrix and vector data structures; set corresponding routines
 59:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 60:   /*
 61:      Create vectors for solution and nonlinear function
 62:   */
 63:   VecCreate(PETSC_COMM_WORLD,&x);
 64:   VecSetSizes(x,PETSC_DECIDE,2);
 65:   VecSetFromOptions(x);
 66:   VecDuplicate(x,&r);

 68:   /*
 69:      Create Jacobian matrix data structure
 70:   */
 71:   MatCreate(PETSC_COMM_WORLD,&J);
 72:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
 73:   MatSetFromOptions(J);
 74:   MatSetUp(J);

 76:   PetscOptionsHasName(NULL,NULL,"-hard",&flg);
 77:   if (!flg) {
 78:     /*
 79:      Set function evaluation routine and vector.
 80:     */
 81:     SNESSetFunction(snes,r,FormFunction1,NULL);

 83:     /*
 84:      Set Jacobian matrix data structure and Jacobian evaluation routine
 85:     */
 86:     SNESSetJacobian(snes,J,J,FormJacobian1,NULL);
 87:   } else {
 88:     SNESSetFunction(snes,r,FormFunction2,NULL);
 89:     SNESSetJacobian(snes,J,J,FormJacobian2,NULL);
 90:   }

 92:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 93:      Customize nonlinear solver; set runtime options
 94:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 95:   /*
 96:      Set linear solver defaults for this problem. By extracting the
 97:      KSP and PC contexts from the SNES context, we can then
 98:      directly call any KSP and PC routines to set various options.
 99:   */
100:   SNESGetKSP(snes,&ksp);
101:   KSPGetPC(ksp,&pc);
102:   PCSetType(pc,PCNONE);
103:   KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);

105:   /*
106:      Set SNES/KSP/KSP/PC runtime options, e.g.,
107:          -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
108:      These options will override those specified above as long as
109:      SNESSetFromOptions() is called _after_ any other customization
110:      routines.
111:   */
112:   SNESSetFromOptions(snes);

114:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115:      Evaluate initial guess; then solve nonlinear system
116:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117:   if (!flg) {
118:     VecSet(x,pfive);
119:   } else {
120:     VecGetArray(x,&xx);
121:     xx[0] = 2.0; xx[1] = 3.0;
122:     VecRestoreArray(x,&xx);
123:   }
124:   /*
125:      Note: The user should initialize the vector, x, with the initial guess
126:      for the nonlinear solver prior to calling SNESSolve().  In particular,
127:      to employ an initial guess of zero, the user should explicitly set
128:      this vector to zero by calling VecSet().
129:   */

131:   SNESSolve(snes,NULL,x);
132:   if (flg) {
133:     Vec f;
134:     VecView(x,PETSC_VIEWER_STDOUT_WORLD);
135:     SNESGetFunction(snes,&f,0,0);
136:     VecView(r,PETSC_VIEWER_STDOUT_WORLD);
137:   }

139:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140:      Free work space.  All PETSc objects should be destroyed when they
141:      are no longer needed.
142:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

144:   VecDestroy(&x)); PetscCall(VecDestroy(&r);
145:   MatDestroy(&J)); PetscCall(SNESDestroy(&snes);
146:   PetscFinalize();
147:   return 0;
148: }
149: /* ------------------------------------------------------------------- */
150: /*
151:    FormFunction1 - Evaluates nonlinear function, F(x).

153:    Input Parameters:
154: .  snes - the SNES context
155: .  x    - input vector
156: .  ctx  - optional user-defined context

158:    Output Parameter:
159: .  f - function vector
160:  */
161: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
162: {
163:   const PetscScalar *xx;
164:   PetscScalar       *ff;

166:   /*
167:    Get pointers to vector data.
168:       - For default PETSc vectors, VecGetArray() returns a pointer to
169:         the data array.  Otherwise, the routine is implementation dependent.
170:       - You MUST call VecRestoreArray() when you no longer need access to
171:         the array.
172:    */
173:   VecGetArrayRead(x,&xx);
174:   VecGetArray(f,&ff);

176:   /* Compute function */
177:   ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
178:   ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;

180:   /* Restore vectors */
181:   VecRestoreArrayRead(x,&xx);
182:   VecRestoreArray(f,&ff);
183:   return 0;
184: }
185: /* ------------------------------------------------------------------- */
186: /*
187:    FormJacobian1 - Evaluates Jacobian matrix.

189:    Input Parameters:
190: .  snes - the SNES context
191: .  x - input vector
192: .  dummy - optional user-defined context (not used here)

194:    Output Parameters:
195: .  jac - Jacobian matrix
196: .  B - optionally different preconditioning matrix
197: .  flag - flag indicating matrix structure
198: */
199: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
200: {
201:   const PetscScalar *xx;
202:   PetscScalar       A[4];
203:   PetscInt          idx[2] = {0,1};

205:   /*
206:      Get pointer to vector data
207:   */
208:   VecGetArrayRead(x,&xx);

210:   /*
211:      Compute Jacobian entries and insert into matrix.
212:       - Since this is such a small problem, we set all entries for
213:         the matrix at once.
214:   */
215:   A[0]  = 2.0*xx[0] + xx[1]; A[1] = xx[0];
216:   A[2]  = xx[1]; A[3] = xx[0] + 2.0*xx[1];
217:   MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES);

219:   /*
220:      Restore vector
221:   */
222:   VecRestoreArrayRead(x,&xx);

224:   /*
225:      Assemble matrix
226:   */
227:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
228:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
229:   if (jac != B) {
230:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
231:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
232:   }
233:   return 0;
234: }

236: /* ------------------------------------------------------------------- */
237: PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
238: {
239:   const PetscScalar *xx;
240:   PetscScalar       *ff;

242:   /*
243:      Get pointers to vector data.
244:        - For default PETSc vectors, VecGetArray() returns a pointer to
245:          the data array.  Otherwise, the routine is implementation dependent.
246:        - You MUST call VecRestoreArray() when you no longer need access to
247:          the array.
248:   */
249:   VecGetArrayRead(x,&xx);
250:   VecGetArray(f,&ff);

252:   /*
253:      Compute function
254:   */
255:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
256:   ff[1] = xx[1];

258:   /*
259:      Restore vectors
260:   */
261:   VecRestoreArrayRead(x,&xx);
262:   VecRestoreArray(f,&ff);
263:   return 0;
264: }
265: /* ------------------------------------------------------------------- */
266: PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
267: {
268:   const PetscScalar *xx;
269:   PetscScalar       A[4];
270:   PetscInt          idx[2] = {0,1};

272:   /*
273:      Get pointer to vector data
274:   */
275:   VecGetArrayRead(x,&xx);

277:   /*
278:      Compute Jacobian entries and insert into matrix.
279:       - Since this is such a small problem, we set all entries for
280:         the matrix at once.
281:   */
282:   A[0]  = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
283:   A[2]  = 0.0;                     A[3] = 1.0;
284:   MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES);

286:   /*
287:      Restore vector
288:   */
289:   VecRestoreArrayRead(x,&xx);

291:   /*
292:      Assemble matrix
293:   */
294:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
295:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
296:   if (jac != B) {
297:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
298:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
299:   }
300:   return 0;
301: }

303: /*TEST

305:    test:
306:       args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
307:       requires: !single

309:    test:
310:       suffix: 2
311:       requires: !single
312:       args:  -snes_monitor_short
313:       output_file: output/ex1_1.out

315:    test:
316:       suffix: 3
317:       args: -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab  -snes_monitor_short
318:       requires: !single
319:       output_file: output/ex1_1.out

321:    test:
322:       suffix: 4
323:       args: -ksp_view_solution ascii:ex1_2_sol.tmp::append  -snes_monitor_short
324:       requires: !single
325:       output_file: output/ex1_1.out

327:    test:
328:       suffix: 5
329:       args: -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab:append  -snes_monitor_short
330:       requires: !single
331:       output_file: output/ex1_1.out

333:    test:
334:       suffix: 6
335:       args: -ksp_view_solution ascii:ex1_2_sol.tmp:default:append  -snes_monitor_short
336:       requires: !single
337:       output_file: output/ex1_1.out

339:    test:
340:       suffix: X
341:       args: -ksp_monitor_short -ksp_type gmres -ksp_gmres_krylov_monitor -snes_monitor_short -snes_rtol 1.e-4
342:       requires: !single x

344: TEST*/