Actual source code: ex42.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  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,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 */
 29:   PetscErrorCode      ierr;
 30:   PetscInt            its;
 31:   PetscScalar         *xx;
 32:   SNESConvergedReason reason;

 34:   PetscInitialize(&argc,&argv,(char*)0,help);

 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,NULL);

 65:   /*
 66:      Set Jacobian matrix data structure and Jacobian evaluation routine
 67:   */
 68:   SNESSetJacobian(snes,J,J,FormJacobian1,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,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",reason>0 ? "CONVERGED" : (char*) 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: {
127:   PetscErrorCode    ierr;
128:   PetscScalar       *ff;
129:   const PetscScalar *xx;

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

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

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

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

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

173:   /*
174:      Get pointer to vector data
175:   */
176:   VecGetArrayRead(x,&xx);

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

189:   /*
190:      Restore vector
191:   */
192:   VecRestoreArrayRead(x,&xx);

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