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: */

 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 */
 45:   PetscMPIInt    size;
 46:   PetscScalar    pfive = .5,*xx;
 47:   PetscBool      flg;

 49:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 50:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 51:   if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Example is only for sequential runs");

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

208:   /*
209:      Get pointer to vector data
210:   */
211:   VecGetArrayRead(x,&xx);

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

222:   /*
223:      Restore vector
224:   */
225:   VecRestoreArrayRead(x,&xx);

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

239: /* ------------------------------------------------------------------- */
240: PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
241: {
242:   PetscErrorCode    ierr;
243:   const PetscScalar *xx;
244:   PetscScalar       *ff;

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

256:   /*
257:      Compute function
258:   */
259:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
260:   ff[1] = xx[1];

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

277:   /*
278:      Get pointer to vector data
279:   */
280:   VecGetArrayRead(x,&xx);

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

291:   /*
292:      Restore vector
293:   */
294:   VecRestoreArrayRead(x,&xx);

296:   /*
297:      Assemble matrix
298:   */
299:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
300:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
301:   if (jac != B) {
302:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
303:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
304:   }
305:   return 0;
306: }

308: /*TEST

310:    test:
311:       args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
312:       requires: !single

314:    test:
315:       suffix: 2
316:       requires: !single
317:       args:  -snes_monitor_short
318:       output_file: output/ex1_1.out

320:    test:
321:       suffix: 3
322:       args: -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab  -snes_monitor_short
323:       requires: !single
324:       output_file: output/ex1_1.out

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

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

338:    test:
339:       suffix: 6
340:       args: -ksp_view_solution ascii:ex1_2_sol.tmp:default:append  -snes_monitor_short
341:       requires: !single
342:       output_file: output/ex1_1.out

344:    test:
345:       suffix: X
346:       args: -ksp_monitor_short -ksp_type gmres -ksp_gmres_krylov_monitor -snes_monitor_short -snes_rtol 1.e-4
347:       requires: !single x

349: TEST*/