Actual source code: ex1.c

petsc-3.3-p7 2013-05-11
  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: */
 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: */
 27: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
 28: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
 29: extern PetscErrorCode FormJacobian2(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
 30: extern PetscErrorCode FormFunction2(SNES,Vec,Vec,void*);

 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:   PetscBool      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);
 87:   MatSetUp(J);

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

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

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

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

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

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

154:   PetscPrintf(PETSC_COMM_WORLD,"number of SNES iterations = %D\n",its);

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

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

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

182:    Output Parameter:
183: .  f - function vector
184:  */
185: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
186: {
187:   PetscErrorCode    ierr;
188:   const PetscScalar *xx;
189:   PetscScalar       *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
203:            code development to parallel implementation.
204:        (1) scatter x into a sequetial vector;
205:        (2) each process evaluates all values of floc; 
206:        (3) scatter floc back to the parallel f.
207:      */
208:     VecScatterBegin(scatter,x,xloc,INSERT_VALUES,SCATTER_FORWARD);
209:     VecScatterEnd(scatter,x,xloc,INSERT_VALUES,SCATTER_FORWARD);

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

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

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

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

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

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

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

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

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

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

296: /* ------------------------------------------------------------------- */
299: PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
300: {
302:   PetscScalar    *xx,*ff;

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

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

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

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

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

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

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