Actual source code: ex1.c

  2: static char help[] = "Newton's method to solve a two-variable system, sequentially.\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:      petsc.h       - base PETSc routines   petscvec.h - vectors
 12:      petscsys.h    - system routines       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: typedef struct {
 20:   Vec         xloc,rloc;    /* local solution, residual vectors */
 21:   VecScatter  scatter;
 22: } AppCtx;

 24: /* 
 25:    User-defined routines
 26: */

 34: int main(int argc,char **argv)
 35: {
 36:   SNES           snes;         /* nonlinear solver context */
 37:   KSP            ksp;          /* linear solver context */
 38:   PC             pc;           /* preconditioner context */
 39:   Vec            x,r;          /* solution, residual vectors */
 40:   Mat            J;            /* Jacobian matrix */
 42:   PetscInt       its;
 43:   PetscMPIInt    size,rank;
 44:   PetscScalar    pfive = .5,*xx;
 45:   PetscTruth     flg;
 46:   AppCtx         user;         /* user-defined work context */
 47:   IS             isglobal,islocal;

 49:   PetscInitialize(&argc,&argv,(char *)0,help);
 50:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 51:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

 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:   if (size > 1){
 70:     VecCreateSeq(PETSC_COMM_SELF,2,&user.xloc);
 71:     VecDuplicate(user.xloc,&user.rloc);

 73:     /* Create the scatter between the global x and local xloc */
 74:     ISCreateStride(MPI_COMM_SELF,2,0,1,&islocal);
 75:     ISCreateStride(MPI_COMM_SELF,2,0,1,&isglobal);
 76:     VecScatterCreate(x,isglobal,user.xloc,islocal,&user.scatter);
 77:     ISDestroy(isglobal);
 78:     ISDestroy(islocal);
 79:   }

 81:   /*
 82:      Create Jacobian matrix data structure
 83:   */
 84:   MatCreate(PETSC_COMM_WORLD,&J);
 85:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
 86:   MatSetFromOptions(J);

 88:   PetscOptionsHasName(PETSC_NULL,"-hard",&flg);
 89:   if (!flg) {
 90:     /* 
 91:      Set function evaluation routine and vector.
 92:     */
 93:     SNESSetFunction(snes,r,FormFunction1,&user);

 95:     /* 
 96:      Set Jacobian matrix data structure and Jacobian evaluation routine
 97:     */
 98:     SNESSetJacobian(snes,J,J,FormJacobian1,PETSC_NULL);
 99:   } else {
100:     if (size != 1) SETERRQ(1,"This case is a uniprocessor example only!");
101:     SNESSetFunction(snes,r,FormFunction2,PETSC_NULL);
102:     SNESSetJacobian(snes,J,J,FormJacobian2,PETSC_NULL);
103:   }

105:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106:      Customize nonlinear solver; set runtime options
107:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108:   /* 
109:      Set linear solver defaults for this problem. By extracting the
110:      KSP, KSP, and PC contexts from the SNES context, we can then
111:      directly call any KSP, KSP, and PC routines to set various options.
112:   */
113:   SNESGetKSP(snes,&ksp);
114:   KSPGetPC(ksp,&pc);
115:   PCSetType(pc,PCNONE);
116:   KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);

118:   /* 
119:      Set SNES/KSP/KSP/PC runtime options, e.g.,
120:          -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
121:      These options will override those specified above as long as
122:      SNESSetFromOptions() is called _after_ any other customization
123:      routines.
124:   */
125:   SNESSetFromOptions(snes);

127:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128:      Evaluate initial guess; then solve nonlinear system
129:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130:   if (!flg) {
131:     VecSet(x,pfive);
132:   } else {
133:     VecGetArray(x,&xx);
134:     xx[0] = 2.0; xx[1] = 3.0;
135:     VecRestoreArray(x,&xx);
136:   }
137:   /*
138:      Note: The user should initialize the vector, x, with the initial guess
139:      for the nonlinear solver prior to calling SNESSolve().  In particular,
140:      to employ an initial guess of zero, the user should explicitly set
141:      this vector to zero by calling VecSet().
142:   */

144:   SNESSolve(snes,PETSC_NULL,x);
145:   SNESGetIterationNumber(snes,&its);
146:   if (flg) {
147:     Vec f;
148:     VecView(x,PETSC_VIEWER_STDOUT_WORLD);
149:     SNESGetFunction(snes,&f,0,0);
150:     VecView(r,PETSC_VIEWER_STDOUT_WORLD);
151:   }

153:   if (!rank){
154:     PetscPrintf(PETSC_COMM_SELF,"number of Newton iterations = %D\n\n",its);
155:   }

157:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158:      Free work space.  All PETSc objects should be destroyed when they
159:      are no longer needed.
160:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

162:   VecDestroy(x); VecDestroy(r);
163:   MatDestroy(J); SNESDestroy(snes);
164:   if (size > 1){
165:     VecDestroy(user.xloc);
166:     VecDestroy(user.rloc);
167:     VecScatterDestroy(user.scatter);
168:   }
169:   PetscFinalize();
170:   return 0;
171: }
172: /* ------------------------------------------------------------------- */
175: /* 
176:    FormFunction1 - Evaluates nonlinear function, F(x).

178:    Input Parameters:
179: .  snes - the SNES context
180: .  x    - input vector
181: .  ctx  - optional user-defined context

183:    Output Parameter:
184: .  f - function vector
185:  */
186: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
187: {
189:   PetscScalar    *xx,*ff;
190:   AppCtx         *user = (AppCtx*)ctx;
191:   Vec            xloc=user->xloc,floc=user->rloc;
192:   VecScatter     scatter=user->scatter;
193:   MPI_Comm       comm;
194:   PetscMPIInt    size,rank;
195:   PetscInt       rstart,rend;

197:   PetscObjectGetComm((PetscObject)snes,&comm);
198:   MPI_Comm_size(comm,&size);
199:   MPI_Comm_rank(comm,&rank);
200:   if (size > 1){
201:     /* 
202:        This is a ridiculous case for testing intermidiate steps from sequential code development to parallel implementation.
203:        (1) scatter x into a sequetial vector;
204:        (2) each process evaluates all values of floc; 
205:        (3) scatter floc back to the parallel f.
206:      */
207:     VecScatterBegin(scatter,x,xloc,INSERT_VALUES,SCATTER_FORWARD);
208:     VecScatterEnd(scatter,x,xloc,INSERT_VALUES,SCATTER_FORWARD);

210:     VecGetOwnershipRange(f,&rstart,&rend);
211:     VecGetArray(xloc,&xx);
212:     VecGetArray(floc,&ff);
213:     ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
214:     ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;
215:     VecRestoreArray(floc,&xx);
216:     VecRestoreArray(xloc,&xx);

218:     VecScatterBegin(scatter,floc,f,INSERT_VALUES,SCATTER_REVERSE);
219:     VecScatterEnd(scatter,floc,f,INSERT_VALUES,SCATTER_REVERSE);
220:   } else {
221:     /*
222:      Get pointers to vector data.
223:        - For default PETSc vectors, VecGetArray() returns a pointer to
224:          the data array.  Otherwise, the routine is implementation dependent.
225:        - You MUST call VecRestoreArray() when you no longer need access to
226:          the array.
227:     */
228:     VecGetArray(x,&xx);
229:     VecGetArray(f,&ff);

231:     /* Compute function */
232:     ff[0] = xx[0]*xx[0] + xx[0]*xx[1] - 3.0;
233:     ff[1] = xx[0]*xx[1] + xx[1]*xx[1] - 6.0;

235:     /* Restore vectors */
236:     VecRestoreArray(x,&xx);
237:     VecRestoreArray(f,&ff);
238:   }
239:   return 0;
240: }
241: /* ------------------------------------------------------------------- */
244: /*
245:    FormJacobian1 - Evaluates Jacobian matrix.

247:    Input Parameters:
248: .  snes - the SNES context
249: .  x - input vector
250: .  dummy - optional user-defined context (not used here)

252:    Output Parameters:
253: .  jac - Jacobian matrix
254: .  B - optionally different preconditioning matrix
255: .  flag - flag indicating matrix structure
256: */
257: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
258: {
259:   PetscScalar    *xx,A[4];
261:   PetscInt       idx[2] = {0,1};

263:   /*
264:      Get pointer to vector data
265:   */
266:   VecGetArray(x,&xx);

268:   /*
269:      Compute Jacobian entries and insert into matrix.
270:       - Since this is such a small problem, we set all entries for
271:         the matrix at once.
272:   */
273:   A[0] = 2.0*xx[0] + xx[1]; A[1] = xx[0];
274:   A[2] = xx[1]; A[3] = xx[0] + 2.0*xx[1];
275:   MatSetValues(*B,2,idx,2,idx,A,INSERT_VALUES);
276:   *flag = SAME_NONZERO_PATTERN;

278:   /*
279:      Restore vector
280:   */
281:   VecRestoreArray(x,&xx);

283:   /* 
284:      Assemble matrix
285:   */
286:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
287:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
288:   if (*jac != *B){
289:     MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
290:     MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
291:   }
292:   return 0;
293: }

295: /* ------------------------------------------------------------------- */
298: PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
299: {
301:   PetscScalar    *xx,*ff;

303:   /*
304:      Get pointers to vector data.
305:        - For default PETSc vectors, VecGetArray() returns a pointer to
306:          the data array.  Otherwise, the routine is implementation dependent.
307:        - You MUST call VecRestoreArray() when you no longer need access to
308:          the array.
309:   */
310:   VecGetArray(x,&xx);
311:   VecGetArray(f,&ff);

313:   /*
314:      Compute function
315:   */
316:   ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
317:   ff[1] = xx[1];

319:   /*
320:      Restore vectors
321:   */
322:   VecRestoreArray(x,&xx);
323:   VecRestoreArray(f,&ff);
324:   return 0;
325: }
326: /* ------------------------------------------------------------------- */
329: PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
330: {
331:   PetscScalar    *xx,A[4];
333:   PetscInt       idx[2] = {0,1};

335:   /*
336:      Get pointer to vector data
337:   */
338:   VecGetArray(x,&xx);

340:   /*
341:      Compute Jacobian entries and insert into matrix.
342:       - Since this is such a small problem, we set all entries for
343:         the matrix at once.
344:   */
345:   A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
346:   A[2] = 0.0;                     A[3] = 1.0;
347:   MatSetValues(*B,2,idx,2,idx,A,INSERT_VALUES);
348:   *flag = SAME_NONZERO_PATTERN;

350:   /*
351:      Restore vector
352:   */
353:   VecRestoreArray(x,&xx);

355:   /*
356:      Assemble matrix
357:   */
358:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
359:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
360:   if (*jac != *B){
361:     MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
362:     MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
363:   }
364:   return 0;
365: }