Actual source code: ex1.c

petsc-3.4.5 2014-06-29
  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: /*
 20:    User-defined routines
 21: */
 22: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
 23: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
 24: extern PetscErrorCode FormJacobian2(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
 25: extern PetscErrorCode FormFunction2(SNES,Vec,Vec,void*);

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

 42:   PetscInitialize(&argc,&argv,(char*)0,help);
 43:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 44:   if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Example is only for sequential runs");

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

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

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

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

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

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

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

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

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

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

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

141:   VecDestroy(&x); VecDestroy(&r);
142:   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: .  ctx  - optional user-defined context

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

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

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

180:   /* Restore vectors */
181:   VecRestoreArrayRead(x,&xx);
182:   VecRestoreArray(f,&ff);
183:   return 0;
184: }
185: /* ------------------------------------------------------------------- */
188: /*
189:    FormJacobian1 - Evaluates Jacobian matrix.

191:    Input Parameters:
192: .  snes - the SNES context
193: .  x - input vector
194: .  dummy - optional user-defined context (not used here)

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

208:   /*
209:      Get pointer to vector data
210:   */
211:   VecGetArrayRead(x,&xx);

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

223:   /*
224:      Restore vector
225:   */
226:   VecRestoreArrayRead(x,&xx);

228:   /*
229:      Assemble matrix
230:   */
231:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
232:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
233:   if (*jac != *B) {
234:     MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
235:     MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
236:   }
237:   return 0;
238: }

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

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

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

265:   /*
266:      Restore vectors
267:   */
268:   VecRestoreArrayRead(x,&xx);
269:   VecRestoreArray(f,&ff);
270:   return 0;
271: }
272: /* ------------------------------------------------------------------- */
275: PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
276: {
277:   const PetscScalar *xx;
278:   PetscScalar       A[4];
279:   PetscErrorCode    ierr;
280:   PetscInt          idx[2] = {0,1};

282:   /*
283:      Get pointer to vector data
284:   */
285:   VecGetArrayRead(x,&xx);

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

297:   /*
298:      Restore vector
299:   */
300:   VecRestoreArrayRead(x,&xx);

302:   /*
303:      Assemble matrix
304:   */
305:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
306:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
307:   if (*jac != *B) {
308:     MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
309:     MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
310:   }
311:   return 0;
312: }