2: static char help[] = "Newton's method for a two-variable system, sequential.\n\n";

  4: /*
  5:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
  6:    file automatically includes:
  7:      petscsys.h       - base PETSc routines   petscvec.h - vectors
  8:      petscmat.h - matrices
  9:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 10:      petscviewer.h - viewers               petscpc.h  - preconditioners
 11:      petscksp.h   - linear solvers
 12: */

 23: #include <petscsnes.h>

 25: /*
 26:    User-defined routines
 27: */
 28: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat,Mat,void*);
 29: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
 30: extern PetscErrorCode FormJacobian2(SNES,Vec,Mat,Mat,void*);
 31: extern PetscErrorCode FormFunction2(SNES,Vec,Vec,void*);

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

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

 48:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 49:      Create nonlinear solver context
 50:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 51:   SNESCreate(PETSC_COMM_WORLD,&snes);

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

 64:   /*
 65:      Create Jacobian matrix data structure
 66:   */
 67:   MatCreate(PETSC_COMM_WORLD,&J);
 68:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
 69:   MatSetFromOptions(J);
 70:   MatSetUp(J);

 72:   PetscOptionsHasName(NULL,NULL,"-hard",&flg);
 73:   if (!flg) {
 74:     /*
 75:      Set function evaluation routine and vector.
 76:     */
 77:     SNESSetFunction(snes,r,FormFunction1,NULL);

 79:     /*
 80:      Set Jacobian matrix data structure and Jacobian evaluation routine
 81:     */
 82:     SNESSetJacobian(snes,J,J,FormJacobian1,NULL);
 83:   } else {
 84:     SNESSetFunction(snes,r,FormFunction2,NULL);
 85:     SNESSetJacobian(snes,J,J,FormJacobian2,NULL);
 86:   }

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

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

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

127:   SNESSolve(snes,NULL,x);
128:   if (flg) {
129:     Vec f;
130:     VecView(x,PETSC_VIEWER_STDOUT_WORLD);
131:     SNESGetFunction(snes,&f,0,0);
132:     VecView(r,PETSC_VIEWER_STDOUT_WORLD);
133:   }

135:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136:      Free work space.  All PETSc objects should be destroyed when they
137:      are no longer needed.
138:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

140:   VecDestroy(&x)); PetscCall(VecDestroy(&r);
141:   MatDestroy(&J)); PetscCall(SNESDestroy(&snes);
142:   PetscFinalize();
143:   return 0;
144: }
145: /* ------------------------------------------------------------------- */
146: /*
147:    FormFunction1 - Evaluates nonlinear function, F(x).

149:    Input Parameters:
150: .  snes - the SNES context
151: .  x    - input vector
152: .  ctx  - optional user-defined context

154:    Output Parameter:
155: .  f - function vector
156:  */
157: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
158: {
159:   const PetscScalar *xx;
160:   PetscScalar       *ff;

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

172:   /* Compute function */
173:   ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
174:   ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;

176:   /* Restore vectors */
177:   VecRestoreArrayRead(x,&xx);
178:   VecRestoreArray(f,&ff);
179:   return 0;
180: }
181: /* ------------------------------------------------------------------- */
182: /*
183:    FormJacobian1 - Evaluates Jacobian matrix.

185:    Input Parameters:
186: .  snes - the SNES context
187: .  x - input vector
188: .  dummy - optional user-defined context (not used here)

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

201:   /*
202:      Get pointer to vector data
203:   */
204:   VecGetArrayRead(x,&xx);

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

215:   /*
216:      Restore vector
217:   */
218:   VecRestoreArrayRead(x,&xx);

220:   /*
221:      Assemble matrix
222:   */
223:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
224:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
225:   if (jac != B) {
226:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
227:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
228:   }
229:   return 0;
230: }

232: /* ------------------------------------------------------------------- */
233: PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
234: {
235:   const PetscScalar *xx;
236:   PetscScalar       *ff;

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

248:   /*
249:      Compute function
250:   */
251:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
252:   ff[1] = xx[1];

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

268:   /*
269:      Get pointer to vector data
270:   */
271:   VecGetArrayRead(x,&xx);

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

282:   /*
283:      Restore vector
284:   */
285:   VecRestoreArrayRead(x,&xx);

287:   /*
288:      Assemble matrix
289:   */
290:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
291:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
292:   if (jac != B) {
293:     MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
294:     MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
295:   }
296:   return 0;
297: }

299: /*TEST

301:    test:
302:       args: -ksp_gmres_cgs_refinement_type refine_always -snes_monitor_short
303:       requires: !single

305:    test:
306:       suffix: 2
307:       requires: !single
308:       args:  -snes_monitor_short
309:       output_file: output/ex1_1.out

311:    test:
312:       suffix: 3
313:       args: -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab  -snes_monitor_short
314:       requires: !single
315:       output_file: output/ex1_1.out

317:    test:
318:       suffix: 4
319:       args: -ksp_view_solution ascii:ex1_2_sol.tmp::append  -snes_monitor_short
320:       requires: !single
321:       output_file: output/ex1_1.out

323:    test:
324:       suffix: 5
325:       args: -ksp_view_solution ascii:ex1_2_sol.tmp:ascii_matlab:append  -snes_monitor_short
326:       requires: !single
327:       output_file: output/ex1_1.out

329:    test:
330:       suffix: 6
331:       args: -ksp_view_solution ascii:ex1_2_sol.tmp:default:append  -snes_monitor_short
332:       requires: !single
333:       output_file: output/ex1_1.out

335:    test:
336:       suffix: X
337:       args: -ksp_monitor_short -ksp_type gmres -ksp_gmres_krylov_monitor -snes_monitor_short -snes_rtol 1.e-4
338:       requires: !single x

340: TEST*/