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 uniprocessor example
  6:    Processors: 1
  7: T*/

  9: /* 
 10:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 11:    file automatically includes:
 12:      petsc.h       - base PETSc routines   petscvec.h - vectors
 13:      petscsys.h    - system routines       petscmat.h - matrices
 14:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 15:      petscviewer.h - viewers               petscpc.h  - preconditioners
 16:      petscksp.h   - linear solvers
 17: */
 18:  #include petscsnes.h

 20: /* 
 21:    User-defined routines
 22: */

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

 43:   PetscInitialize(&argc,&argv,(char *)0,help);
 44:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 45:   if (size != 1) SETERRQ(1,"This is a uniprocessor example only!");

 47:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 48:      Create nonlinear solver context
 49:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 51:   SNESCreate(PETSC_COMM_WORLD,&snes);

 53:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 54:      Create matrix and vector data structures; set corresponding routines
 55:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 57:   /*
 58:      Create vectors for solution and nonlinear function
 59:   */
 60:   VecCreateSeq(PETSC_COMM_SELF,2,&x);
 61:   VecDuplicate(x,&r);

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

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

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

 86:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 87:      Customize nonlinear solver; set runtime options
 88:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

126:   SNESSolve(snes,PETSC_NULL,x);
127:   SNESGetIterationNumber(snes,&its);
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:   PetscPrintf(PETSC_COMM_SELF,"number of Newton iterations = %D\n\n",its);

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

142:   VecDestroy(x); VecDestroy(r);
143:   MatDestroy(J); SNESDestroy(snes);

145:   PetscFinalize();
146:   return 0;
147: }
148: /* ------------------------------------------------------------------- */
151: /* 
152:    FormFunction1 - Evaluates nonlinear function, F(x).

154:    Input Parameters:
155: .  snes - the SNES context
156: .  x - input vector
157: .  dummy - optional user-defined context (not used here)

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

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

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

183:   /*
184:      Restore vectors
185:   */
186:   VecRestoreArray(x,&xx);
187:   VecRestoreArray(f,&ff);

189:   return 0;
190: }
191: /* ------------------------------------------------------------------- */
194: /*
195:    FormJacobian1 - Evaluates Jacobian matrix.

197:    Input Parameters:
198: .  snes - the SNES context
199: .  x - input vector
200: .  dummy - optional user-defined context (not used here)

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

213:   /*
214:      Get pointer to vector data
215:   */
216:   VecGetArray(x,&xx);

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

228:   /*
229:      Restore vector
230:   */
231:   VecRestoreArray(x,&xx);

233:   /*
234:      Assemble matrix
235:   */
236:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
237:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);

239:   return 0;
240: }


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

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

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

267:   /*
268:      Restore vectors
269:   */
270:   VecRestoreArray(x,&xx);
271:   VecRestoreArray(f,&ff);

273:   return 0;
274: }
275: /* ------------------------------------------------------------------- */
278: PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
279: {
280:   PetscScalar    *xx,A[4];
282:   PetscInt       idx[2] = {0,1};

284:   /*
285:      Get pointer to vector data
286:   */
287:   VecGetArray(x,&xx);

289:   /*
290:      Compute Jacobian entries and insert into matrix.
291:       - Since this is such a small problem, we set all entries for
292:         the matrix at once.
293:   */
294:   A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
295:   A[2] = 0.0;                     A[3] = 1.0;
296:   MatSetValues(*jac,2,idx,2,idx,A,INSERT_VALUES);
297:   *flag = SAME_NONZERO_PATTERN;

299:   /*
300:      Restore vector
301:   */
302:   VecRestoreArray(x,&xx);

304:   /*
305:      Assemble matrix
306:   */
307:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
308:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);

310:   return 0;
311: }