Actual source code: ex1.c

  1: /*$Id: ex1.c,v 1.26 2001/08/07 03:04:16 balay Exp $*/

  3: static char help[] = "Newton's method to solve a two-variable system, sequentially.\n\n";

  5: /*T
  6:    Concepts: SNES^basic uniprocessor example
  7:    Processors: 1
  8: T*/

 10: /* 
 11:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 12:    file automatically includes:
 13:      petsc.h       - base PETSc routines   petscvec.h - vectors
 14:      petscsys.h    - system routines       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: */
 19:  #include petscsnes.h

 21: /* 
 22:    User-defined routines
 23: */
 24: extern int FormJacobian1(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
 25: extern int FormFunction1(SNES,Vec,Vec,void*);
 26: extern int FormJacobian2(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
 27: extern int FormFunction2(SNES,Vec,Vec,void*);

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

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

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

 50:   SNESCreate(PETSC_COMM_WORLD,&snes);

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

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

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

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

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

 84:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 85:      Customize nonlinear solver; set runtime options
 86:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

124:   SNESSolve(snes,x);
125:   SNESGetIterationNumber(snes,&its);
126:   if (flg) {
127:     Vec f;
128:     VecView(x,PETSC_VIEWER_STDOUT_WORLD);
129:     SNESGetFunction(snes,&f,0,0);
130:     VecView(r,PETSC_VIEWER_STDOUT_WORLD);
131:   }

133:   PetscPrintf(PETSC_COMM_SELF,"number of Newton iterations = %d\n\n",its);

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

140:   VecDestroy(x); VecDestroy(r);
141:   MatDestroy(J); SNESDestroy(snes);

143:   PetscFinalize();
144:   return 0;
145: }
146: /* ------------------------------------------------------------------- */
149: /* 
150:    FormFunction1 - Evaluates nonlinear function, F(x).

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

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

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

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

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

187:   return 0;
188: }
189: /* ------------------------------------------------------------------- */
192: /*
193:    FormJacobian1 - Evaluates Jacobian matrix.

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

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

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

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

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

230:   /*
231:      Assemble matrix
232:   */
233:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
234:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);

236:   return 0;
237: }


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

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

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

264:   /*
265:      Restore vectors
266:   */
267:   VecRestoreArray(x,&xx);
268:   VecRestoreArray(f,&ff);

270:   return 0;
271: }
272: /* ------------------------------------------------------------------- */
275: int FormJacobian2(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
276: {
277:   PetscScalar *xx,A[4];
278:   int    ierr,idx[2] = {0,1};

280:   /*
281:      Get pointer to vector data
282:   */
283:   VecGetArray(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(*jac,2,idx,2,idx,A,INSERT_VALUES);
293:   *flag = SAME_NONZERO_PATTERN;

295:   /*
296:      Restore vector
297:   */
298:   VecRestoreArray(x,&xx);

300:   /*
301:      Assemble matrix
302:   */
303:   MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
304:   MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);

306:   return 0;
307: }