2: static char help[] = "Newton's method for a two-variable system, sequential.\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: */

 29: #include <petscsnes.h>

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

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

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

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

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

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

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

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

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

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

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

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


143:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144:      Free work space.  All PETSc objects should be destroyed when they
145:      are no longer needed.
146:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

157:    Input Parameters:
158: .  snes - the SNES context
159: .  x    - input vector
160: .  ctx  - optional user-defined context

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

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

181:   /* Compute function */
182:   ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
183:   ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;

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

194:    Input Parameters:
195: .  snes - the SNES context
196: .  x - input vector
197: .  dummy - optional user-defined context (not used here)

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

211:   /*
212:      Get pointer to vector data
213:   */
214:   VecGetArrayRead(x,&xx);

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

225:   /*
226:      Restore vector
227:   */
228:   VecRestoreArrayRead(x,&xx);

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

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

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

259:   /*
260:      Compute function
261:   */
262:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
263:   ff[1] = xx[1];

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

280:   /*
281:      Get pointer to vector data
282:   */
283:   VecGetArrayRead(x,&xx);

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

294:   /*
295:      Restore vector
296:   */
297:   VecRestoreArrayRead(x,&xx);

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




314: /*TEST

316:    test:
317:       args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
318:       requires: !single

320:    test:
321:       suffix: 2
322:       requires: !single
323:       args:  -snes_monitor_short
324:       output_file: output/ex1_1.out

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

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

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

344:    test:
345:       suffix: 6
346:       args: -ksp_view_solution ascii:ex1_2_sol.tmp:default:append  -snes_monitor_short
347:       requires: !single
348:       output_file: output/ex1_1.out

350:    test:
351:       suffix: X
352:       args: -ksp_monitor_short -ksp_type gmres -ksp_gmres_krylov_monitor -snes_monitor_short -snes_rtol 1.e-4
353:       requires: !single x

355: TEST*/