Actual source code: ex42.c

petsc-3.3-p7 2013-05-11
  2: static char help[] = "Newton's method to solve a two-variable system that comes from the Rosenbrock function.\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: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
 20: extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);

 24: int main(int argc,char **argv)
 25: {
 26:   SNES           snes;         /* nonlinear solver context */
 27:   Vec            x,r;          /* solution, residual vectors */
 28:   Mat            J;            /* Jacobian matrix */
 30:   PetscInt       its;
 31:   PetscScalar    *xx;
 32:   SNESConvergedReason reason;

 34:   PetscInitialize(&argc,&argv,(char *)0,help);
 35: 
 36:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 37:      Create nonlinear solver context
 38:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 39:   SNESCreate(PETSC_COMM_WORLD,&snes);

 41:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 42:      Create matrix and vector data structures; set corresponding routines
 43:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 44:   /*
 45:      Create vectors for solution and nonlinear function
 46:   */
 47:   VecCreate(PETSC_COMM_WORLD,&x);
 48:   VecSetSizes(x,PETSC_DECIDE,2);
 49:   VecSetFromOptions(x);
 50:   VecDuplicate(x,&r);

 52:   /*
 53:      Create Jacobian matrix data structure
 54:   */
 55:   MatCreate(PETSC_COMM_WORLD,&J);
 56:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2);
 57:   MatSetFromOptions(J);
 58:   MatSetUp(J);

 60:   /* 
 61:      Set function evaluation routine and vector.
 62:   */
 63:   SNESSetFunction(snes,r,FormFunction1,PETSC_NULL);

 65:   /* 
 66:      Set Jacobian matrix data structure and Jacobian evaluation routine
 67:   */
 68:   SNESSetJacobian(snes,J,J,FormJacobian1,PETSC_NULL);

 70:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 71:      Customize nonlinear solver; set runtime options
 72:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 73:   SNESSetFromOptions(snes);

 75:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 76:      Evaluate initial guess; then solve nonlinear system
 77:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 78:   VecGetArray(x,&xx);
 79:   xx[0] = -1.2; xx[1] = 1.0;
 80:   VecRestoreArray(x,&xx);

 82:   /*
 83:      Note: The user should initialize the vector, x, with the initial guess
 84:      for the nonlinear solver prior to calling SNESSolve().  In particular,
 85:      to employ an initial guess of zero, the user should explicitly set
 86:      this vector to zero by calling VecSet().
 87:   */

 89:   SNESSolve(snes,PETSC_NULL,x);
 90:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);
 91:   SNESGetIterationNumber(snes,&its);
 92:   SNESGetConvergedReason(snes,&reason);
 93:   /*
 94:      Some systems computes a residual that is identically zero, thus converging
 95:      due to FNORM_ABS, others converge due to FNORM_RELATIVE.  Here, we only
 96:      report the reason if the iteration did not converge so that the tests are
 97:      reproducible.
 98:   */
 99:   PetscPrintf(PETSC_COMM_WORLD,"%s number of SNES iterations = %D\n\n",reason>0?"CONVERGED":SNESConvergedReasons[reason],its);

101:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102:      Free work space.  All PETSc objects should be destroyed when they
103:      are no longer needed.
104:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

106:   VecDestroy(&x); VecDestroy(&r);
107:   MatDestroy(&J); SNESDestroy(&snes);
108:   PetscFinalize();
109:   return 0;
110: }
111: /* ------------------------------------------------------------------- */
114: /* 
115:    FormFunction1 - Evaluates nonlinear function, F(x).

117:    Input Parameters:
118: .  snes - the SNES context
119: .  x    - input vector
120: .  ctx  - optional user-defined context

122:    Output Parameter:
123: .  f - function vector
124:  */
125: PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
126: {
128:   PetscScalar    *xx,*ff;

130:   /*
131:     Get pointers to vector data.
132:     - For default PETSc vectors, VecGetArray() returns a pointer to
133:     the data array.  Otherwise, the routine is implementation dependent.
134:     - You MUST call VecRestoreArray() when you no longer need access to
135:     the array.
136:   */
137:   VecGetArray(x,&xx);
138:   VecGetArray(f,&ff);

140:   /* Compute function */
141:   ff[0] = -2.0 + 2.0*xx[0] + 400.0*xx[0]*xx[0]*xx[0] - 400.0*xx[0]*xx[1];
142:   ff[1] = -200.0*xx[0]*xx[0] + 200.0*xx[1];

144:   /* Restore vectors */
145:   VecRestoreArray(x,&xx);
146:   VecRestoreArray(f,&ff);
147:   return 0;
148: }
149: /* ------------------------------------------------------------------- */
152: /*
153:    FormJacobian1 - Evaluates Jacobian matrix.

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

160:    Output Parameters:
161: .  jac - Jacobian matrix
162: .  B - optionally different preconditioning matrix
163: .  flag - flag indicating matrix structure
164: */
165: PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *dummy)
166: {
167:   PetscScalar    *xx,A[4];
169:   PetscInt       idx[2] = {0,1};

171:   /*
172:      Get pointer to vector data
173:   */
174:   VecGetArray(x,&xx);

176:   /*
177:      Compute Jacobian entries and insert into matrix.
178:       - Since this is such a small problem, we set all entries for
179:         the matrix at once.
180:   */
181:   A[0] = 2.0 + 1200.0*xx[0]*xx[0] - 400.0*xx[1];
182:   A[1] = -400.0*xx[0];
183:   A[2] = -400.0*xx[0];
184:   A[3] = 200;
185:   MatSetValues(*B,2,idx,2,idx,A,INSERT_VALUES);
186:   *flag = SAME_NONZERO_PATTERN;

188:   /*
189:      Restore vector
190:   */
191:   VecRestoreArray(x,&xx);

193:   /* 
194:      Assemble matrix
195:   */
196:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
197:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
198:   if (*jac != *B){
199:     MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
200:     MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
201:   }
202:   return 0;
203: }