Actual source code: ex42.c

petsc-3.11.4 2019-09-28
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*/



 10: /*
 11:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 12:    file automatically includes:
 13:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 14:      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: extern PetscErrorCode FormJacobian1(SNES,Vec,Mat,Mat,void*);
 22: 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);if (ierr) return ierr;

 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 ierr;
110: }
111: /* ------------------------------------------------------------------- */
112: /*
113:    FormFunction1 - Evaluates nonlinear function, F(x).

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

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

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

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

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

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

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

169:   /*
170:      Get pointer to vector data
171:   */
172:   VecGetArrayRead(x,&xx);

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

185:   /*
186:      Restore vector
187:   */
188:   VecRestoreArrayRead(x,&xx);

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



204: /*TEST

206:    test:
207:       args: -snes_monitor_short -snes_max_it 1000
208:       requires: !single

210: TEST*/